Closed azev77 closed 3 years ago
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)
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")
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.
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
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")
I assume you have set the domain you used for the Chebyshev polynomials appropriately.
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")
I'll record some unsolicited, non-urgent comments here (obviously ignore the ones you want...)
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[:,:]
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.
Oh, and computing Euler equation errors is on the todo list.
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")
Start 20% below ss: Start 20% above ss:
A few issues:
compute_steady_state()
gives the wrong steady state, so I used the closed form here