NLsolve.jl: It seems like the solver does nothing


#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:

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.