[ANN] Update to SolveDSGE

RJDennis · March 7, 2020, 2:21pm

Hi All,

I’ve recently released an update to SolveDSGE. The update allows the model that is to be solved to be stored separately in a file that gets read and processed prior to solution. In addition, models can be solved using first-, second-, or third-order perturbation, or by projection methods. The solution obtained from one method can easily be used as the initialization for another, allowing accurate solutions to be obtained more quickly. There is a user guide to help you get started and examples of how to use SolveDSGE are provided in the examples folder. Existing code using the package should be unaffected by the update.

Albert_Zevelev · March 7, 2020, 7:21pm

Thanks this is really useful!
I’m having trouble plotting the transition path from an initial state using perturbation, but it is beautiful w/ projection.

Consider a deterministic growth model stored in “model.txt”:

states:
cap
end

jumps:
cons
end

shocks:
end

parameters:
betta = 0.99
sigma = 1.1
delta = 0.025
alpha = 0.30
end

equations:
cap(+1) = (1.0 - delta)*cap + cap^alpha - cons
cons^(-sigma) = betta*cons(+1)^(-sigma)*(1.0 - delta + alpha*cap(+1)^(alpha - 1.0))
end

Now Julia:

using SolveDSGE
filename = "model.txt";
path = joinpath(@__DIR__,filename);
dsge = get_model(path)

x = [3.05, 0.7] #cap/cons
tol = 1e-8; maxiters = 1000;
ss = compute_steady_state(dsge,x,tol,maxiters)

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

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

#Projection Cheby
P = ChebyshevSchemeDet(ss,chebyshev_nodes,[11],4,[23.0; 19.5],tol,1e-6,maxiters)
PP = ChebyshevSchemeDet(ss,chebyshev_nodes,[21],4,[23.0; 19.5],tol,1e-6,maxiters)
PPP = ChebyshevSchemeDet(ss,chebyshev_nodes,[31],5,[23.0; 19.5],tol,1e-6,maxiters)

soln_nla = solve_model(dsge,P)
soln_nlb = solve_model(dsge,soln_so,PP)
soln_nlc = solve_model(dsge,soln_nlb,PPP)

#Projection Smolyak #can't solve right now.
#Proj PiecewiseLinear  #can't solve right now

sim_data = simulate(soln_nla, [ ss[1]/2 ],100) #capital below ss
sim_data = simulate(soln_nla, [ ss[1]*2 ],100) #capital above ss

sim_data = simulate(soln_to, [ ss[1]/2 ],100)
sim_data = simulate(soln_to, [ ss[1]*2 ],100)

plot(1:T, sim_data', label = ["cap" "cons"])
plot!(ss, seriestype = :hline, label=false)

Plotting one at a time.
Projection methods look good, but Perturbation methods start from random initial points.

Also, is there any way to enter model info/equations using Julia ?
It would be really convenient & Julia makes it easy to modify parameters etc…
Thanks again!!!

Albert_Zevelev · March 7, 2020, 7:36pm

Also, is it possible to use your package to solve finite horizon models?

RJDennis · March 7, 2020, 9:44pm

I’ve just pushed an update that I believe fixes this problem. You can access it from master or wait until the registration request is merged — probably an hour or so.

RJDennis · March 7, 2020, 9:45pm

Sorry, I haven’t added any functionality for solving finite horizon models with this update.

Albert_Zevelev · March 8, 2020, 7:01pm

@RJDennis the deterministic simulations look great now!
I’m having some issues w/ simulations for stochastic models.

#Stochastic NCGM
states:
k, z
end

jumps:
c
end

shocks:
ε
end

parameters:
β = 0.99
σ = 1.1
δ = 0.025
α = 0.30
ρ = 0.8
σ_ε = 0.01 
end

equations:
k(+1)  = (1.0 - δ)*k + exp(z)*k^α - c
c^(-σ) = β*c(+1)^(-σ)*(1.0 - δ + α*exp(z(+1))*k(+1)^(α - 1.0))
z(+1)  = ρ*z + σ_ε*ε
end

Julia program:

filename = "StochasticNCGM.txt";
path = joinpath(@__DIR__,filename); m = get_model(path);
xg = [0.05, 3.05, 0.7]; tol = 1e-8; maxiters = 1000;
ss = compute_steady_state(m,xg,tol,maxiters); #Steady state #zss, kss, css

Third order perturbation gives an error message when simulating

P3 = PerturbationScheme(ss,1.0,"third");
s_P3 = solve_model(m,P3)
sim_data = simulate(s_P3, x0b, T)

UndefVarError: hx not defined

P2 = PerturbationScheme(ss,1.0,"second");
s_P2 = solve_model(m,P2);
C1 = ChebyshevSchemeStoch(ss,chebyshev_nodes,[21,21],9,3,[0.07 23.0; -0.07 19.5],tol,1e-6,maxiters);
s_C1 = solve_model(m,C1);
L1 = SmolyakSchemeStoch(ss,chebyshev_gauss_lobatto,9,4,[0.07 23.0; -0.07 19.5],tol,1e-6,maxiters);
s_L1 = solve_model(m,L1);

T=100; x0a= ss[1:2]*2;
#
sim_data = simulate(s_P2, x0a, T)
plot(1:T, sim_data', label = hcat(m.variables...))
plot!(ss, seriestype = :hline, label=false)
#
sim_data = simulate(s_C1, x0a, T)
plot(1:T, sim_data', label = hcat(m.variables...))
plot!(ss, seriestype = :hline, label=false)
#
sim_data = simulate(s_L1, x0a, T)
plot(1:T, sim_data', label = hcat(m.variables...))
plot!(ss, seriestype = :hline, label=false)

#
T_irf=50; ix_shock=1; n_MC= 10000;
#
pos_imps, neg_imps = impulses(s_P2,T_irf,ix_shock,n_MC)
plot(1:50, pos_imps', title = "Positive IRF", label = hcat(m.variables...))
plot(1:50, neg_imps', title = "Negative IRF", label = hcat(m.variables...))
#
pos_imps, neg_imps = impulses(s_C1,T_irf,ix_shock,n_MC)
plot(1:50, pos_imps', title = "Positive IRF", label = hcat(m.variables...))
plot(1:50, neg_imps', title = "Negative IRF", label = hcat(m.variables...))

The IRFs look great but the simulations don’t seem to start above steady state.

Albert_Zevelev · March 8, 2020, 7:29pm

I think some of this is due to the default: burn_in = 1000 currently hard-coded in all stochastic simulations.

1 It could be nice to keep burn_in = 1000 as the default but also to allow the user to customize.
For example:

g(x; burn_in=1000) = x + burn_in
g(5) #1005
g(5, burn_in=3) #8

Here:

function simulate(soln::R,initial_state::Array{T,1},n::S; burn_in=1000)

@sglyon does something similar: (CLMMJuliaPythonMatlab/sim_resids.jl at master · sglyon/CLMMJuliaPythonMatlab · GitHub) setting default nburn=200

2 Similarly it could be nice to allow user to customize the seed or possibly generate multiple simulations

function simulate(soln::R,initial_state::Array{T,1},n::S; burn_in=1000, seed=123456)

3 currently the only option for generating shocks is randn this can also be customized or maybe even allow a different distribution for each shock?

I submitted a PR for burn_in. This won’t cause any breaks bc the default burn_in remains the same. Likewise the default arguments for simulate remain: simulate(soln,initial_state,n)

RJDennis · March 8, 2020, 10:27pm

I’ve pushed an update. I removed burn_in from the simulate function. That way the user can generate their own burn_in period through their choice of how many observations to simulate and then just remove however many initial periods they want. And I made the initial seed an optional argument. I knew that me hard-coding the length of the burn_in period and setting the seed would raise questions and need to be changed even before I made the update, but I ran out of energy.

At this point I have kept randn for the random number generator, primarily because normality is used in the quadrature when solving the models. I’ll give some thought to this and maybe change it and allow some customization at a later stage. I wasn’t able to replicate the error with the third-order perturbation. Let me know if you encounter it again.

Two little things: You are probably just doing it for illustration, but the initial states you are using are generally outside the domain used to solve the model, which could easily lead to inaccuracies or errors. Finally, your suggestions are most welcome, but some of them might be better raised as issues on the github page rather than in general discourse.

jmatovu · April 1, 2020, 6:47am

I am having challenges with get_model function. Cant locate it in the directory and hence cant get files. Where could this function be?