NLsolve.jl: It seems like the solver does nothing

New to Julia
1

After I had some problems in the beginning, I sort of gave up. Now I want to give it another try and it would be really nice to get to the solution of the problem this time.

I have a basic root-finding problem. It is an economics problem and broadly speaking I want to know what “optimal” consumption is depending on two state variables, capital and productivity.

Here is the code:

# This example: Neoclassical Growth Model
using QuantEcon

# Parameters
global γ = 2;    # Risk aversion
global α = 0.3;  # Capital share
global β = 0.98; # Discount factor

# First, create a grid using the QuantEcon toolbox
K_grid = collect(linspace(0.1, 10, 5))
z_grid = collect(linspace(0.9, 1.1, 3))
p_grid = [1/3 1/3 1/3] # iid shocks and they are equally probable

grid=gridmake(K_grid,z_grid)

# First, Rootfinding and linear policy function for c.


using NLsolve

f! = function(α_N,real_res,grid) # we want to find roots for f
  global γ,α,β
  K_next    = grid[:,2].*grid[:,1].^α - (grid*α_N)
  next_grid = gridmake(K_next,z_grid)
  c_next    = (next_grid*α_N)
  RHS       = β*reshape(c_next.^(-γ).*α.*(next_grid[:,1]).^(α-1).*next_grid[:,2],15,3)*p_grid'
  res       = (grid*α_N).^(-γ) - RHS
  real_res  = grid'*res
end

g!(α_N,real_res) = f!(α_N,real_res,grid)
guess = [0.01, 0.1]
res = nlsolve(g!,guess)
g!(guess,[0,0])

The output from the second to last line is:

grafik

nlsolve claims that the residual f(x) is already very small and convergence has been achieved. However, if I execute g!(guess,[0,0]), I get the following:

grafik

So the residual is by no means close to zero.

Does anyone know how I can go about it?

2

It doesn’t look like you modify α_N?

3

I don’t quite understand. I am trying to find the α_N that solves the equation and thus minimizes real_res.

4

Oh, it’s flipped from what I remember. But re-read the docs:

https://github.com/JuliaNLSolvers/NLsolve.jl#simple-example

You don’t output the residual. Rather, you should update the residual inplace. The easiest way to make the change here is to do real_res .= grid'*res. Of course, you can make that a lot better by making it not allocate, like At_mul_B!(real_res,grid,res) I think is the call for that here.

5

Thanks, it worked. The working code is now:

# Simple implementation of several colocation algorithms to solve
# models globally that are otherwise a bit clunky to solve.

# This example: Neoclassical Growth Model
using QuantEcon

# Parameters
global γ = 2;    # Risk aversion
global α = 0.3;  # Capital share
global β = 0.98; # Discount factor

# First, create a grid using the QuantEcon toolbox
K_grid = collect(linspace(0.1, 3, 5))
z_grid = collect(linspace(0.9, 1.1, 3))
p_grid = [1/3 1/3 1/3] # iid shocks and they are equally probable

grid=gridmake(K_grid,z_grid)

# First, Rootfinding and linear policy function for c.


using NLsolve
α_N = [α*β/(1+α*β) 0]'
real_res = [10.0 10.0]'
f! = function(α_N,real_res,grid) # we want to find roots for f
  global γ,α,β
  @show α_N # To see progress
  K_next    = grid[:,2].*grid[:,1].^α - (grid*α_N)
  next_grid = gridmake(K_next,z_grid)
  c_next    = next_grid*α_N
  RHS       = β*reshape(c_next.^(-γ).*α.*(next_grid[:,1]).^(α-1).*next_grid[:,2],15,3)*p_grid'
  res       = (grid*α_N).^(-γ) - RHS
  real_res  .= grid'*res
end

g!(α_N,real_res) = f!(α_N,real_res,grid)

guess = [α*β/(1+α*β) 0.3]'
res = nlsolve(g!,guess)

# Solution
α_N=res.zero
6

Hey all,

quick question regarding your final code: When I re-run it I get the error message:
DomainError:
Exponentiation yielding a complex result requires a complex argument

do you have any idea why that is?

Thanks

7

I will look at it in the following week and let you know then. I also realized that the code is not working anymore, but didn’t look into it too much so far.

8

This error message arises because c_next probably takes negative values. Taking roots of negative values yields complex numbers, which cause problems in the other mathematical operations. One way to go about it is to calculate RHS only for positive values of c_next and set res to high values where c_next is negative. The same might be true with respect to K_next of course.