Useful example for the readme: comparison w/ closed form for deterministic NCGM

[azev77]

using SolveDSGE, Plots;
cd("my_directory")

# write .txt model in Julia
# deterministic Neoclassical Growth Model 
open("Det_NCGM.txt", "w") do io
    println(io,"# ncgm\n")
    println(io,"states:\nk \nend\n")
    println(io,"jumps:\nc \nend\n")
    println(io,"shocks:\nend\n")
    println(io,"parameters:\nβ=0.99, σ=1.0, δ=1.0, α=.3 \nend\n")
    println(io,"equations:")
    println(io,"k(+1) = (1-δ)*k + k^α - c")
    println(io,"(c(+1))^σ = β*(c^σ)*(1-δ + α*(k(+1)^(α-1.0)) )")
    println(io,"end")
end;

###########################################################################
path = joinpath(@__DIR__,"Det_NCGM.txt")
process_model(path)
dsge = retrieve_processed_model(path)
x = ones(dsge.number_variables)  #Guess ["k", "c"]
tol = 1e-8; maxiters = 1_000; #Int(1e+3);
ss = compute_steady_state(dsge,x,tol,maxiters) # wrong ss

kss=(((1/.99) + 1 -1)/.3)^(1/(.3-1.0))
css = kss^.3 - kss
ss = [kss, css]
ss = [kss, css]

###########################################################################
# Perturbation, Order = 1,2,3
###########################################################################
cutoff = 1.0; # Perturbation HP 
N   = PerturbationScheme(ss,cutoff,"first")
NN  = PerturbationScheme(ss,cutoff,"second")
NNN = PerturbationScheme(ss,cutoff,"third")

soln_fo = solve_model(dsge,N)
soln_so = solve_model(dsge,NN)
soln_to = solve_model(dsge,NNN)

###########################################################################
# Projection, ChebyshevSchemeDet # not working...
###########################################################################

###########################################################################
# Analysis 
###########################################################################
T = 100; # num Sim pds 
sss = ss[1:dsge.number_states]; #ss for states only 

sim_d1 = simulate(soln_fo, sss*.8, T)
sim_d2 = simulate(soln_so, sss*.8, T)
sim_d3 = simulate(soln_to, sss*.8, T)
# Closed Form Sol 
k_cf    = zeros(T+1) 
k_cf[1] = sss[1]*.8
c_cf    = zeros(T)
for t in 1:T
    k_cf[t+1] = (.3*.99) * (k_cf[t]^.3)  #k(+1) = (.3*.99) * (k^.3)
    c_cf[t]   = (1-.3*.99) * (k_cf[t]^.3) #c     = (1-.3*.99) * (k^.3)
end

p1 = plot()
plot!(k_cf[2:T+1] , lab="k_t closed form")
plot!(c_cf        , lab="c_t closed form")
plot!(1:T, sim_d1')
plot!(1:T, sim_d2')
plot!(1:T, sim_d3')

sim_d1 = simulate(soln_fo, sss*1.2, T)
sim_d2 = simulate(soln_so, sss*1.2, T)
sim_d3 = simulate(soln_to, sss*1.2, T)
k_cf    = zeros(T+1)
k_cf[1] = sss[1]*1.2
c_cf    = zeros(T)
for t in 1:T
    k_cf[t+1] = (.3*.99) * (k_cf[t]^.3)
    c_cf[t]   = (1-.3*.99) * (k_cf[t]^.3)
end
p1 = plot()
plot!(k_cf[2:T+1] , lab="k_t closed form")
plot!(c_cf        , lab="c_t closed form")
plot!(1:T, sim_d1')
plot!(1:T, sim_d2')
plot!(1:T, sim_d3')

Start 20% below ss:

Start 20% above ss:

A few issues:

1. compute_steady_state() gives the wrong steady state, so I used the closed form here
2. the projection methods don't seem to work for this problem, not sure why,,,
3. At some point it would be great if we could automatically plot the policy fcns (then I can compare w/ sol from vfi etc)

[RJDennis]

I didn't have a problem with solving the model using Chebyshev polynomials or with the labeling issue in the plots (your second and third issues). See the code below. Regarding the steady state, your right that the starting values you used didn't produce the correct steady state. It may be the the max-iteration limit was reached in which case I can have the function issue a warning. This is something I always planned to do, but never got around to.

ngm = retrieve_processed_model(processed_path)

tol = 1e-8
maxiters = 1000
x = [0.2,0.4]
ss = compute_steady_state(ngm, x, tol, maxiters)

cutoff = 1.0
N   = PerturbationScheme(ss,cutoff,"first")
NN  = PerturbationScheme(ss,cutoff,"second")
NNN = PerturbationScheme(ss,cutoff,"third")

soln_fo = solve_model(ngm, N)
soln_so = solve_model(ngm, NN)
soln_to = solve_model(ngm, NNN)

P = ChebyshevSchemeDet(ss,chebyshev_nodes,[21],6,[0.22; 0.13],tol,1e-6,maxiters)

soln_nla = solve_model(ngm,P)
soln_nlb = solve_model(ngm,soln_to,P)

T = 100 # num Sim pds 

sss = ss[1:ngm.number_states] #ss for states only 

sim_d1   = simulate(soln_fo, sss*.8, T)
sim_d2   = simulate(soln_so, sss*.8, T)
sim_d3   = simulate(soln_to, sss*.8, T)
sim_proj = simulate(soln_nla, sss*.8, T)

using Plots

p1 = plot(1:T,sim_d1',title = "First Order Pert.",labels   = reshape(ngm.variables,1,ngm.number_variables))
p2 = plot(1:T,sim_d2',title = "Second Order Pert.",labels   = reshape(ngm.variables,1,ngm.number_variables))
p3 = plot(1:T,sim_d3',title = "Third Order Pert.",labels   = reshape(ngm.variables,1,ngm.number_variables))
p4 = plot(1:T,sim_proj',title = "Proj. (Chebyshev poly)",labels = reshape(ngm.variables,1,ngm.number_variables))

plot(p1,p2,p3,p4)

[azev77]

Thanks. My third issue was about being able to automatically plot the policy function given by each solution.

Example for the first order perturbation:

#closed form for capital
p_cf(x) = (α*β) * (x^α)

# x(t+1) = hx*x(t)          # First Order 
p_fo(x; soln=soln_fo) = soln.hx[1] * (x-soln.hbar[1]) + soln.hbar[1]

plot(legend=:topleft);
plot!(0.01:.01:2*kss, p_cf, lab="Closed form")
plot!(0.01:.01:2*kss, p_fo, lab="Pert 1st")

[RJDennis]

I misunderstood. From the figures you showed I thought the problem was with the labelling. I suppose what you have in mind is plotting specific slices through the policy function, allowing one state variale to vary at a time. I can produce a function that will return the values for the decision variables for any given state which may help with what you want, but I don't want to make any plotting package a dependency.

[azev77]

Sorry I wasn't clear, I just plotted what I meant in my comment above.

My question wasn't about plotting per-se. (I agree it's a bad idea to add dependencies.)
It was about automatically returning the policy function to the user.
The way I manually created: p_fo(x; soln=soln_fo) = soln.hx[1] * (x-soln.hbar[1]) + soln.hbar[1]

You prob can re-use the policy function later in your package when simulating etc

[azev77]

This can show off the potential benefits of global vs local methods.
A simple example can easily illustrate what kinds of method are less/more accurate far away from the steady state etc..
(might also be used as a good starting point for vfi/pfi methods etc)

# closed form for capital: k_{t+1} = f(k_t)
p_cf(x) = (α*β) * (x^α)

# x(t+1) = hx*x(t)          #+ k*v(t+1)
#   y(t) = gx*x(t)
# x(t+1) - xss = hx*(x(t) - xss)          
p_fo(x; s=soln_fo) = s.hx[1] * (x-s.hbar[1]) + s.hbar[1]

#using LinearAlgebra
# x(t+1) = hx*x(t) + (1/2)*[kron(I,x(t))]'hxx*[kron(I,x(t))]
function p_so(x; s=soln_so)
    nx = length(s.hbar)
    hxx = Matrix(reshape(s.hxx',nx*nx,nx)')
    #
    a1 = s.hx*(x-s.hbar[1])
    a2 = s.hx*(0.) + (1/2)*hxx*kron(x-s.hbar[1],x-s.hbar[1])
    a3 = a1 + a2 + s.hbar  
    return a3[1]  
end

# ChebyshevSolutionDet
using ChebyshevApprox
function p_ch(x; soln=soln_cha)
    w = chebyshev_weights_threaded(soln.variables[1],soln.nodes,soln.order,soln.domain)
    return chebyshev_evaluate(w,[x],soln.order,soln.domain)  
end

function p_chc(x; soln=soln_chc)
    w = chebyshev_weights_threaded(soln.variables[1],soln.nodes,soln.order,soln.domain)
    return chebyshev_evaluate(w,[x],soln.order,soln.domain)  
end


plot(legend=:topleft);
plot!(0.01:.01:2*kss, p_cf, lab="Closed form")
plot!(0.01:.01:2*kss, p_fo, lab="Pert 1st")
plot!(0.01:.01:2*kss, p_so, lab="Pert 2nd")
plot!(0.01:.01:2*kss, p_ch, lab="Proj Cheby")
plot!(0.01:.01:2*kss, p_chc,lab="Proj Cheby")

[RJDennis]

I assume you have set the domain you used for the Chebyshev polynomials appropriately.

[azev77]

Here's what I did:

using SolveDSGE, Plots;
cd("mydir")

# deterministic Neoclassical Growth Model 
open("Det_NCGM.txt", "w") do io
    println(io,"# ncgm\n")
    println(io,"states:\nk \nend\n")
    println(io,"jumps:\nc \nend\n")
    println(io,"shocks:\nend\n")
    println(io,"parameters:\nβ=0.99, σ=1.0, δ=1.0, α=.3 \nend\n")
    println(io,"equations:")
    println(io,"k(+1) = (1-δ)*k + k^α - c")
    println(io,"(c(+1))^σ = β*(c^σ)*(1-δ + α*(k(+1)^(α-1.0)) )")
    println(io,"end")
end;

###########################################################################
path = joinpath(@__DIR__,"Det_NCGM.txt")
process_model(path)
ngm = retrieve_processed_model(path)

###########################################################################
# Steady State 
###########################################################################
tol = 1e-8
maxiters = 1000
x = [0.25,0.25] 
#x = [0.,0.] Error. 
#x = [0.5,0.5] #Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
ss = compute_steady_state(ngm, x, tol, maxiters)
sss = ss[1:ngm.number_states] #ss for states only 

#ss closed form compare. 
β=0.99; σ=1.0; δ=1.0; α=.3;
kss=(((1/β) + 1 -δ)/α)^(1/(α-1.0)); 
css = kss^α - kss;
ss ≈ [kss, css]                 #true

# 
n_thr = 4; #num threads to use. 
dom = [sss*2.0; sss*0.01]
tf = 1e-8; tv = 1e-6; #inner loop tol #outer loop tol

###########################################################################
# Perturbation, Order = 1,2,3
###########################################################################
g0 = ss;          # point about which to perturb the model (generally the steady state) 
cutoff = 1.0      # param that separates unstable from stable eigenvalues
# order = "first" # order of the perturbation
N   = PerturbationScheme(g0,cutoff,"first")   # Klein (2000)
NN  = PerturbationScheme(g0,cutoff,"second")  # Gomme and Klein (2011)
NNN = PerturbationScheme(g0,cutoff,"third")   # Binning (2013) w refinement from Levintal (2017)

soln_fo = solve_model(ngm, N)
soln_so = solve_model(ngm, NN)
soln_to = solve_model(ngm, NNN)

###########################################################################
# Projection, ChebyshevSchemeDet
###########################################################################
g0 = ss; #initial guess 
ng = chebyshev_nodes; #node generator: chebyshev_nodes/chebyshev_extrema
nn1 = [21]; #number of nodes for each state var #vec integers 
#nq = 9; #number quadrature points. Int. #Don't need if deterministic!!!!
no1 = 6; #order of Cheby poly 

P    = ChebyshevSchemeDet(g0,ng,nn1,no1,dom,tf,tv,maxiters)
soln_cha = solve_model(ngm,P,n_thr)

###########################################################################
# Projection, SmolyakSchemeDet
# Projection, PiecewiseLinearSchemeStoch
# Compose 
###########################################################################
g0 = ss; #initial guess 
ng = chebyshev_gauss_lobatto; #node generator: chebyshev_gauss_lobatto/clenshaw_curtis_equidistant
nl1 = 3; #number of layers to be used in the approximation
#nq = 9; #number quadrature points. Int. #Don't need if deterministic!!!!
no1 = 6; #order of Cheby poly 

S = SmolyakSchemeDet(g0,ng,nl1,dom[:,:],tf,tv,maxiters)
soln_sa = solve_model(ngm,S,n_thr)

# PL 
g0 = ss; #initial guess 
nn1 = [21]; #number of nodes for each state var #vec integers 
#nq = 9; #number quadrature points. Int. #Don't need if deterministic!!!!

PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
soln_pl = solve_model(ngm,PL,n_thr)

###########################################################################
# Analysis
###########################################################################
T=10 # num Sim pds 
ω = 1.0 - 0.90; # start 90% below steady-state

sim_d1   = simulate(soln_fo, sss*ω,  T)
sim_d2   = simulate(soln_so, sss*ω,  T)
sim_d3   = simulate(soln_to, sss*ω,  T)
sim_d4   = simulate(soln_cha, sss*ω, T)
sim_d5   = simulate(soln_sa, sss*ω,  T)
sim_d6   = simulate(soln_pl, sss*ω,  T)
# Closed Form Sol 
k_cf    = zeros(T+1) 
k_cf[1] = sss[1]*ω
c_cf    = zeros(T)
for t in 1:T
    k_cf[t+1] = (α*β) * (k_cf[t]^α)  #k(+1) = (.3*.99) * (k^.3)
    c_cf[t]   = (1-α*β) * (k_cf[t]^α) #c     = (1-.3*.99) * (k^.3)
end
cf = [k_cf[1:T] c_cf]
# sim returns k[1:T]    NOT     k[1:T+1]

lv = reshape(ngm.variables,1,ngm.number_variables)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p1 = plot!(1:T,sim_d1', title="Pert. 1st Order", lab=lv)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p2 = plot!(1:T,sim_d2', title="Pert. 2nd Order", lab=lv)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p3 = plot!(1:T,sim_d3', title="Pert. 3rd Order", lab=lv)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p4 = plot!(1:T,sim_d4', title="Proj. Cheby 1",   lab=lv)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p5 = plot!(1:T,sim_d5', title="Proj. Smol  1",   lab=lv)

p  = plot();
p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
p6 = plot!(1:T,sim_d6', title="Proj. PL    1",   lab=lv)

plot(p1,p2,p3,p4,p5,p6, size=1.25 .* (600, 400))

Observe: 3rd order is closer to closed form than 2nd order, & 2nd closed than 1st.
I'd expect Levintal's 5th order to be even closer...
Projection methods here are very impressive!

Now plot the policy functions over the same domain as used to solve model:

# closed form for capital: k_{t+1} = f(k_t)
p_cf(x) = (α*β) * (x^α)

p_fo(x; s=soln_fo) = s.hx[1] * (x-s.hbar[1]) + s.hbar[1]

function p_so(x; s=soln_so)
    nx = length(s.hbar)
    hxx = Matrix(reshape(s.hxx',nx*nx,nx)')
    #
    a1 = s.hx*(x-s.hbar[1])
    a2 = s.hx*(0.) + (1/2)*hxx*kron(x-s.hbar[1],x-s.hbar[1])
    a3 = a1 + a2 + s.hbar  
    return a3[1]  
end

using ChebyshevApprox
function p_ch(x; soln=soln_cha)
    w = chebyshev_weights_threaded(soln.variables[1],soln.nodes,soln.order,soln.domain)
    return chebyshev_evaluate(w,[x],soln.order,soln.domain)  
end

s_gr = sss[1]*0.01:.01:sss[1]*2.0
plot(legend=:topleft, size=1.25 .* (600, 400) );
plot!(s_gr, p_cf, lab="Closed form")
plot!(s_gr, p_fo, lab="Pert 1st")
plot!(s_gr, p_so, lab="Pert 2nd")
plot!(s_gr, p_ch, lab="Proj Cheby")

Now plot policy functions w/ upper bound higher than domain used to solve model:

s_gr = sss[1]*0.01:.01:1.0
plot(legend=:topleft, size=1.25 .* (600, 400) );
plot!(s_gr, p_cf, lab="Closed form")
plot!(s_gr, p_fo, lab="Pert 1st")
plot!(s_gr, p_so, lab="Pert 2nd")
plot!(s_gr, x->max(p_ch(x),0.05), lab="Proj Cheby")

[azev77]

I'll record some unsolicited, non-urgent comments here (obviously ignore the ones you want...)

• Rendahl 2017 applies Linear Time Iteration to two very simple difference equations
  I feel like these extremely simple remind users that the pkg just solves systems of (non-linear, stochastic) difference equations & might have an appeal to a wider range of users
• Having the solvers return the entire policy function in the solution object might simplify your code.
  For example simulating is "just" iterating on the policy function, instead of coding a separate simulation method for each solver.
  It could also make it easier to compute Euler equation errors etc..
• In principal you can use almost any function approx method for projection. (That's what I learned from JFV.)
  Have you considered using neural nets from Flux.jl?
  Flux can potentially allow users to solve very high dimensional problems, very quickly.
  E.g. Flux is used to solve DEs.
  @mcreel have you tried anything like this?
  For example, DsgeSolves.jl uses neural networks for projection...
• There is a small inconsistency seen in the above example:
  P = ChebyshevSchemeDet(g0,ng,nn1,no1,dom,tf,tv,maxiters)
S = SmolyakSchemeDet(g0,ng,nl1,dom[:,:],tf,tv,maxiters)
PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
  For the Smolyak, the domain needs to be converted from a vector to matrix dom[:,:]

[RJDennis]

A few quick reactions, first, as I mentioned above, I can see the usefulness of generating function for the decision rules and law-of-motions for the states for users, and that I can use such functions in my code to reduce some repetition, so I will do that. Second, I have imposed some restrictions on what can be used for function approximation through the way that I have done the quadrature. You seem to be focused on deterministic models so you don't see the benefits to how the quadrature in SolveDSGE.jl is done. Something like Flux.jl could be used for deterministic models, but I think that I would need to change the way the quadrature is done to use it for stochastic models. Third, I'll sort out the vector/matrix thing for the Smolyak case, just for consistency.

[RJDennis]

Oh, and computing Euler equation errors is on the todo list.

[azev77]

My guess is the standardized decision rule might make it easier to compute generic EE errors.
For the deterministic NCGM above:

function ee(x, pol)
    xp  = pol(x)
    xpp = pol(xp)
    #
    c  = (1-δ)*x  + x^α  - xp
    cp = (1-δ)*xp + xp^α - xpp 
    #
    err = β*(c^σ)*(1-δ + α*(xp^(α-1.0)) )  - (cp)^σ
    return err 
end

s_gr = sss[1]*0.01:.01:sss[1]*2.25
plot(legend=:bottomright, size=1.25 .* (600, 400) );
plot!(s_gr, ee.(s_gr, p_cf), lab="Closed form")
plot!(s_gr, ee.(s_gr, p_fo), lab="Pert 1st")
plot!(s_gr, ee.(s_gr, p_so), lab="Pert 2nd")
plot!(s_gr, ee.(s_gr, p_ch), lab="Proj Cheby")