rjdennis / solvedsge.jl Goto Github PK
View Code? Open in Web Editor NEWA Julia package to solve, simulate, and analyze nonlinear DSGE models.
License: MIT License
A Julia package to solve, simulate, and analyze nonlinear DSGE models.
License: MIT License
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
using SolveDSGE
ERROR: LoadError: LoadError: UndefVarError: FloatingPoint not defined
in include_from_node1(::String) at ./loading.jl:488
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_from_node1(::String) at ./loading.jl:488
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in eval(::Module, ::Any) at ./boot.jl:234
in eval(::Module, ::Any) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in require(::Symbol) at ./loading.jl:415
in require(::Symbol) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
while loading /Users/boris/.julia/v0.5/SolveDSGE/src/types.jl, in expression starting on line 3
while loading /Users/boris/.julia/v0.5/SolveDSGE/src/SolveDSGE.jl, in expression starting on line 3
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
I can't help but ask, can this package solve the occasionally binds model?hh
@JuliaRegistrator register()
How about using CUDA.jl to speed up the calculation in smolyak method? (the speed gain is quite promising according to github.com/ikarib/smolyak
Hi,
I am wondering if you could indicate a reference for the model of Example 3: Three shock New Keynesian Model, please?
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
I am trying to assign parameters to a model, along the lines of model2a.txt in the test directory. The model works fine when parameters are given in the .txt file, but errors when they are assigned. The model and solution file are below, the assignment or not is determined by un-commenting the relevant block in the .txt file, and by commenting/uncommenting the assign_parameters line in the .jl file.
The beginning part of the error message, when assigning. is
The model's variables are now in this order: ["z", "η", "k", "y", "c", "n", "r", "w", "MUC", "MUL"]
The following parameters do not have values assigned: ["α", "β", "δ", "γ", "ρ₁", "σ₁", "ρ₂", "σ₂", "ψ"]
ERROR: LoadError: UndefVarError: pax not defined
Stacktrace:
[1] nlsolve_static_equations(f::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SolveDSGE.va is a missing pax.
I am using SolveDSGE v0.4.6. By the way, the package is very nice and fast, and seems to agree with what I used to get using Dynare and Octave. Any help much appreciated. Thanks, M.
Here are the two files:
.txt file:
`# CK 2 shock model from "Indirect Likelihood Inference"
states:
k, z, η
end
jumps:
y, c, n, r, w, MUC, MUL
end
shocks:
ϵ
u
end
parameters:
α = 0.33
β = 0.99
δ = 0.025
γ = 2.0
ρ₁ = 0.9
σ₁ = 0.02
ρ₂ = 0.7
σ₂ = 0.01
ψ = 3.4174389608451694
end
#parameters:
#α
#β
#δ
#γ
#ρ₁
#σ₁
#ρ₂
#σ₂
#ψ
#end
equations:
MUC = c^(-γ)
MUL = ψexp(η)
r = α * exp(z) * k^(α-1) * n^(1-α)
w = (1-α)exp(z) k^α * n^(-α)
MUC = βMUC(+1) * (1 + r(+1) - δ)
MUL/MUC = w
z(+1) = ρ₁z + σ₁u
η(+1) = ρ₂η + σ₂ϵ
y = exp(z) * (k^α) * (n^(1-α))
k(+1) = y - c + (1-δ)*k
end
`
.jl solution file:
`using SolveDSGE
filename = "CK.txt"
path = joinpath(@DIR,filename)
process_model(path)
processed_filename = "CK_processed.txt"
processed_path = joinpath(@DIR,processed_filename)
dsge = retrieve_processed_model(processed_path)
function ParamsAndSS()
α = 0.33
β = 0.99
δ = 0.025
γ = 2.0
ρ₁ = 0.9
σ₁ = 0.02
ρ₂ = 0.7
σ₂ = 0.01
nss = 1.0/3.0
c1 = ((1/β + δ - 1)/α)^(1/(1-α))
kss = nss/c1
iss = δkss;
yss = kss^α * nss^(1-α)
css = yss - iss;
MUCss = css^(-γ)
rss = α * kss^(α-1) * nss^(1-α)
wss = (1-α) (kss)^α * nss^(-α)
MULss = wss*MUCss
ψ = (css^(-γ)) * (1-α) * (kss^α) * (nss^(-α))
p = [α, β, δ, γ, ρ₁, σ₁, ρ₂, σ₂, ψ]
ss = [0.0, 0.0, kss, yss, css, nss, rss, wss, MUCss, MULss]
return p, ss
end
tol = 1e-8
maxiters = 1000
p, ss = ParamsAndSS()
#dsge = assign_parameters(dsge, p)
ss = compute_steady_state(dsge, ss, tol, maxiters)
scheme = PerturbationScheme(ss, 1.0, "third")
solution = solve_model(dsge, scheme)
`
Hi.
I'm writing a post on solving dynamic economic models in the Julia ecosystem.
I'm interested in 4 broad categories of models:
I have multiple solutions for the deterministic models.
Is it possible to use your package to solve finite horizon models?
@JuliaRegistrator register()
I'm getting an "applicable method may be too new: running in world age X, while current world is Y" error when using the solve_model function
inside another function (which I need in order for instance to solve repeatedly the model for different parameter values).
It seems to be an issue caused by running the retrieve_processed_model
function in global scope or inside a function.
A simple example to replicate the issue is as follows (taken from your model3b
example model):
using SolveDSGE
filename = "model3b.txt"
path3 = joinpath(@__DIR__,filename)
process_model(path3)
processed_filename = "model3b_processed.txt"
processed_path3 = joinpath(@__DIR__,processed_filename)
function testmeplz(processed_path3)
dsge3 = retrieve_processed_model(processed_path3)
x3 = [0.05, 3.05, 0.7, 0.7]
tol = 1e-8
maxiters = 1000
ss3 = compute_steady_state(dsge3,x3,tol,maxiters)
N = PerturbationScheme(ss3,1.0,"first")
soln_fo3 = solve_model(dsge3,N)
return soln_fo3
end
testmeplz(processed_path3)
This simple example generates the same kind of error I get, but for the compute_steady_state
function. When solving my model, I get it with the solve_model
function because I have a customized function that solves the steady state exactly, instead of relying on a solver.
If in the previous example I pasted you move dsge3 = retrieve_processed_model(processed_path3)
outside the testmeplz
function, things work as you would expect. This makes me think the issue has to do with running retrieve_processed_model
in global or local scope.
Do you know how I can solve this issue? Is there something you can do on your side?
Many thanks,
Alberto
@JuliaRegistrator register()
@JuliaRegistrator register()
The consumption saving problem has no steady state (in general).
I'd expect the perturbation methods not to work here, but projection methods should be fine b/c they don't need a ss.
It doesn't seem to be working well:
using SolveDSGE, Plots;
cd("mydir")
# deterministic con-sav
open("Det_CS.txt", "w") do io
println(io,"# cs\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, r= 0.05 \nend\n")
println(io,"equations:")
println(io,"k(+1) = (1+r)*k - c")
println(io,"(c(+1))^σ = β*(1 + r)*(c^σ)")
println(io,"end")
end;
path = joinpath(@__DIR__,"Det_CS.txt")
process_model(path) # The model's variables are now in this order: ["k", "c"]
ngm = retrieve_processed_model(path)
β=0.99; σ=1.0; r=0.05; #r= (1/0.99) - 1 #r=0.05;
g0 = ones(2)
dom = [2.0; 0.01] #domain for state vars. 2*n_s array. #row 1 upper values. row 2 lower values
tf = 1e-8; tv = 1e-6; #inner loop tol #outer loop tol
n_thr = 4; #num threads to use.
###########################################################################
# Projection, ChebyshevSchemeDet
###########################################################################
ng = chebyshev_nodes; #node generator: chebyshev_nodes/chebyshev_extrema
nn1 = [21]; #number of nodes for each state var #vec integers
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
###########################################################################
ng = chebyshev_gauss_lobatto; #node generator: chebyshev_gauss_lobatto/clenshaw_curtis_equidistant
nl1 = 3; #number of layers to be used in the approximation
no1 = 6; #order of Cheby poly
S = SmolyakSchemeDet(g0,ng,nl1,dom[:,:],tf,tv,maxiters)
soln_sa = solve_model(ngm,S,n_thr)
nn1 = [21]; #number of nodes for each state var #vec integers
PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
soln_pl = solve_model(ngm,PL,n_thr)
###########################################################################
# Analysis
###########################################################################
T=20
k0 = [1.5]; # initial state
#c_t = (1-β)*(β^t)*(1+r)^t * a_0 #log-ut
#a_t = (β^t)*(1+r)^t * a_0
sim_d4 = simulate(soln_cha, k0, T)
sim_d5 = simulate(soln_sa, k0, T)
sim_d6 = simulate(soln_pl, k0, T)
# Closed Form Sol
k_cf = zeros(T+1)
k_cf[1] = k0[1]
c_cf = zeros(T)
for t in 1:T
k_cf[t+1] = (β*(1+r))^(t+1) * k0[1]
c_cf[t] = (1-β) * (β*(1+r))^(t+1) * k0[1]
end
cf = [k_cf[1:T] c_cf]
lv = reshape(ngm.variables,1,ngm.number_variables)
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)
#p7 = plot!(1:T,sim_d7', title="Proj. & Pert", lab=lv)
#p8 = plot!(1:T,sim_d8', title="Proj. & Pert", lab=lv)
plot(p4,p5,p6, size=1.25 .* (600, 400))
This returns the wrong solution & the solvers don't give feedback about correctness...
@JuliaRegistrator register()
Consider the deterministic NGM:
k(+1) = (1-δ)*k + k^α - c
1 = β*((c/c(+1))^σ)*(1-δ + α*(k(+1)^(α-1.0)) )
when I use log-transforms (as in Levintal) both Perturbation & Projection methods work:
exp(log(k(+1))) = (1-δ)*exp(log(k)) + exp(log(k))^α - exp(log(c))
1 = β*( ( exp(log(c)) / exp(log(c(+1)))) ^σ)*(1-δ + α*(exp(log(k(+1)))^(α-1.0)) )
when I add r
, Perturbation methods still work but Projection methods give errors
exp(log(k(+1))) = (1-δ)*exp(log(k)) + exp(log(k))^α - exp(log(c))
exp(log(r)) = (1-δ) + α*(exp(log(k))^(α-1.0))
1 = β*( ( exp(log(c)) / exp(log(c(+1)))) ^σ)*exp(log(r(+1)))
Julia code:
using SolveDSGE; cd("mydir");
open("Det_NCGM.txt", "w") do io
println(io,"# ncgm\n")
println(io,"states:\nk \nend\n")
println(io,"jumps:\nc, r \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,"exp(log(k(+1))) = (1-δ)*exp(log(k)) + exp(log(k))^α - exp(log(c))")
println(io,"exp(log(r)) = (1-δ) + α*(exp(log(k))^(α-1.0)) ")
println(io,"1 = β*( ( exp(log(c)) / exp(log(c(+1)))) ^σ)*exp(log(r(+1)))")
println(io,"end")
end;
path = joinpath(@__DIR__,"Det_NCGM.txt")
process_model(path)
ngm = retrieve_processed_model(path)
tol = 1e-8; maxiters = 1000; x = .25*ones(ngm.number_variables);
ss = compute_steady_state(ngm, x, tol, maxiters)
ss = ss.zero
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;
rss= α*(kss)^(α-1);
ss ≈ [kss, css, rss] #true
# 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
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
dom = [sss*2.0; sss*0.01] #domain for state vars. 2*n_s array. #row 1 upper vals. row 2 lower vals
tf = 1e-8; tv = 1e-6; #inner loop tol #outer loop tol
# ChebyshevSchemeDet
ng = chebyshev_nodes; #node generator: chebyshev_nodes/chebyshev_extrema
nn1 = 21 * ones(ngm.number_states) .|> Int64
no1 = 6; #order of Cheby poly
P = ChebyshevSchemeDet(g0,ng,nn1,no1,dom,tf,tv,maxiters)
# Smolyak
ng = chebyshev_gauss_lobatto; #node generator: chebyshev_gauss_lobatto/clenshaw_curtis_equidistant
nl1 = 3; #number of layers to be used in the approximation
no1 = 6; #order of Cheby poly
S = SmolyakSchemeDet(g0,ng,nl1,dom,tf,tv,maxiters)
# PL
nn1 = [21]; #number of nodes for each state var #vec integers
PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
So far everything works.
When I try to solve using Chebyhsev I get:
julia> soln_cha = solve_model(ngm,P)
ERROR: DomainError with -0.040396881204646656:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math .\math.jl:33
[2] log(x::Float64)
@ Base.Math .\special\log.jl:285
[3] (::SolveDSGE.var"#projection_equations#211"{Vector{Float64}, Vector{Vector{Float64}}, Int64, Vector{Float64}, typeof(chebyshev_evaluate)})(f::Vector{Float64}, x::Vector{Float64})
@ SolveDSGE C:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\SolveDSGE_jl_ProjectionPerturbation\AZ\Det_NCGM_processed.txt:87
[4] value!!(obj::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, x::Vector{Float64})
@ NLSolversBase C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:166
[5] value!!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:164 [inlined]
[6] value!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:28 [inlined]
[7] trust_region_(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:174
[8] trust_region (repeats 2 times)
@ C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:235 [inlined]
[9] nlsolve(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#29#31", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:26
[10] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:xtol, :iterations), Tuple{Float64, Int64}}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:52
[11] solve_nonlinear(model::SolveDSGE.REModel{Int64, String}, scheme::ChebyshevSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:649
[12] solve_model(model::SolveDSGE.REModel{Int64, String}, scheme::ChebyshevSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:5124
[13] top-level scope
@ REPL[41754]:1
When I try to solve using Smolyak I get:
julia> soln_sa = solve_model(ngm,S)
ERROR: DomainError with -0.010377766333114088:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math .\math.jl:33
[2] log(x::Float64)
@ Base.Math .\special\log.jl:285
[3] (::SolveDSGE.var"#projection_equations#211"{Vector{Float64}, Vector{Vector{Float64}}, Matrix{Int64}, Vector{Float64}, typeof(smolyak_evaluate)})(f::Vector{Float64}, x::Vector{Float64})
@ SolveDSGE C:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\SolveDSGE_jl_ProjectionPerturbation\AZ\Det_NCGM_processed.txt:87
[4] value!!(obj::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, x::Vector{Float64})
@ NLSolversBase C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:166
[5] value!!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:164 [inlined]
[6] value!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:28 [inlined]
[7] trust_region_(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:174
[8] trust_region (repeats 2 times)
@ C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:235 [inlined]
[9] nlsolve(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#29#31", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:26
[10] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:xtol, :iterations), Tuple{Float64, Int64}}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:52
[11] solve_nonlinear(model::SolveDSGE.REModel{Int64, String}, scheme::SmolyakSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:2392
[12] solve_model(model::SolveDSGE.REModel{Int64, String}, scheme::SmolyakSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:5124
[13] top-level scope
@ REPL[41760]:1
When I try to solve using PL I get:
julia> soln_pl = solve_model(ngm,PL)
ERROR: DomainError with -0.010377766333114088:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math .\math.jl:33
[2] log(x::Float64)
@ Base.Math .\special\log.jl:285
[3] (::SolveDSGE.var"#projection_equations_pl#212"{Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{Float64}, typeof(piecewise_linear_evaluate)})(f::Vector{Float64}, x::Vector{Float64})
@ SolveDSGE C:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\SolveDSGE_jl_ProjectionPerturbation\AZ\Det_NCGM_processed.txt:108
[4] value!!(obj::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, x::Vector{Float64})
@ NLSolversBase C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:166
[5] value!!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:164 [inlined]
[6] value!
@ C:\Users\azevelev\.julia\packages\NLSolversBase\geyh3\src\interface.jl:28 [inlined]
[7] trust_region_(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:174
[8] trust_region (repeats 2 times)
@ C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\solvers\trust_region.jl:235 [inlined]
[9] nlsolve(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#29#31", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:26
[10] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:xtol, :iterations), Tuple{Float64, Int64}}})
@ NLsolve C:\Users\azevelev\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:52
[11] solve_nonlinear(model::SolveDSGE.REModel{Int64, String}, scheme::PiecewiseLinearSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:3737
[12] solve_model(model::SolveDSGE.REModel{Int64, String}, scheme::PiecewiseLinearSchemeDet{Float64, Int64})
@ SolveDSGE C:\Users\azevelev\.julia\packages\SolveDSGE\5G9HJ\src\solution_functions.jl:5124
[13] top-level scope
@ REPL[41764]:1
I've seen these complex argument
errors a few times before using SolveDSGE.jl
...
I would like to ask advice on setting the number of grids for high dimensional problems. I am now dealing with a 6 dimensional Bellman equation and would like to ask how many grids are appropriate?
where i can find "Blanchard_kahn_form" ? when i ran, the code show it is not defined. thank you for your kind help.
@JuliaRegistrator register()
Consider a discrete-time, deterministic firm investment problem w/ adjustment costs:
This problem has a closed form solution (for comparison):
There are 3 regions of legal parameters where capital is growing/constant/shrinking:
(the lower limit for z
is based on the Second Order Condition, the upper limit for z
is based on the transversality condition)
This model has no steady-state except for the knife-edge case z = r + δ
where every initial capital level is the steady state.
The sol is accurate (except Cheby/Smolyak) for params in the decreasing capital region:
The sols are no good for params in the constant capital region (closed form sols in blue/orange):
The sols are no good for params in the growing capital region (closed form sols in blue/orange):
They're not even in the correct direction...
I was hoping one of the benefits of global methods is they can handle models like these.
Can this package (eventually) solve these types of problems? with shocks?
Here is my code:
using SolveDSGE, Plots;
# Param
p=1.00; δ=0.10; θ=2.00; r=0.15;
z= r + δ -.05 # k_t→0 less accurate farther away from r + δ
z= r + δ # k_t→k_0
z= r + δ +.02 # k_t→ ∞
open("Det_INV.txt", "w") do io
println(io,"# inv\n")
println(io,"states:\nk \nend\n")
println(io,"jumps:\ni \nend\n")
println(io,"shocks:\nend\n")
println(io,"parameters:\nδ=0.10, θ=2.00, r=0.15, z=0.20 \nend\n")
println(io,"equations:")
println(io,"k(+1) = (1-δ)*k + i")
println(io,"1 + θ*(i/k - δ) = (1/(1+r))*(z-δ + (θ/2.0)*(i(+1)/k(+1) - δ)^2.0 + 1 + θ*(i(+1)/k(+1) - δ))")
println(io,"end")
end;
path = joinpath(@__DIR__,"Det_INV.txt")
process_model(path)
ngm = retrieve_processed_model(path)
tol = 1e-8; maxiters = 1000;
x = .25*ones(ngm.number_variables)
ss = compute_steady_state(ngm, x, tol, maxiters)
ss = ss.zero
sss = ss[1:ngm.number_states] #ss for states only
# 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)
# ProjectionSchemes: ChebyshevSchemes/SmolyakSchemes/PiecewiseLinearSchemes
g0 = ss; #initial guess
dom = [sss*2.0; sss*0.01] #domain for state vars. 2*n_s array. #row 1 upper vals. row 2 lower vals
tf = 1e-8; #inner loop tol
tv = 1e-6; #outer loop tol
n_thr = 4; #num threads to use.
# Projection, ChebyshevSchemeDet
ng = chebyshev_nodes; #node generator: chebyshev_nodes/chebyshev_extrema
nn1 = 21 * ones(ngm.number_states) .|> Int64
nn2 = 71 * ones(ngm.number_states) .|> Int64
no1 = 6; #order of Cheby poly
P = ChebyshevSchemeDet(g0,ng,nn1,no1,dom,tf,tv,maxiters)
PP = ChebyshevSchemeDet(g0,ng,nn1,4, dom,tf,tv,maxiters)
PPP = ChebyshevSchemeDet(g0,ng,nn2,6, dom,tf,tv,maxiters)
soln_cha = solve_model(ngm,P)
soln_cha = solve_model(ngm,P,n_thr)
soln_chb = solve_model(ngm,PP,n_thr)
soln_chc = solve_model(ngm,PPP,n_thr)
# Smolyak
ng = chebyshev_gauss_lobatto; #node generator: chebyshev_gauss_lobatto/clenshaw_curtis_equidistant
nl1 = 3; #number of layers to be used in the approximation
no1 = 6; #order of Cheby poly
S = SmolyakSchemeDet(g0,ng,nl1,dom,tf,tv,maxiters)
soln_sa = solve_model(ngm,S,n_thr)
# PL
nn1 = [21]; #number of nodes for each state var #vec integers
PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
soln_pl = solve_model(ngm,PL,n_thr)
# Compose
soln_c1 = solve_model(ngm,soln_to,PPP)
soln_c2 = solve_model(ngm,soln_c1,PP)
T = 100 # num Sim pds
ω = 1.0; # distance from steady-state
sss[1]=1.0; # K_0
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)
sim_d7 = simulate(soln_c1, sss*ω, T)
# Closed Form Sol
Q=(1+r)*((1+r*θ)-√θ*√(2*(r+δ-z)+θ*r^2))
gk=(Q/(1+r) -1)*(1/θ)
#
k_cf = zeros(T+1)
k_cf[1] = sss[1]*ω
i_cf = zeros(T)
for t in 1:T
k_cf[t+1] = (k_cf[t]) * (1+gk) #k(+1) =
i_cf[t] = (k_cf[t]) * (δ+gk) #i =
end
cf = [k_cf[1:T] i_cf]
# sim returns k[1:T] NOT k[1:T+1]
lv = reshape(ngm.variables,1,ngm.number_variables)
p = plot(legend=:bottomright);
p0 = plot!(1:T,cf, title="Closed Form", lab="");
p1 = plot(p0, 1:T,sim_d1', title="Pert. 1st Order", lab=lv)
p2 = plot(p0, 1:T,sim_d2', title="Pert. 2nd Order", lab=lv)
p3 = plot(p0, 1:T,sim_d3', title="Pert. 3rd Order", lab=lv)
p4 = plot(p0, 1:T,sim_d4', title="Proj. Cheby 1", lab=lv)
p5 = plot(p0, 1:T,sim_d5', title="Proj. Smol 1", lab=lv)
p6 = plot(p0, 1:T,sim_d6', title="Proj. PL 1", lab=lv)
p7 = plot(p0, 1:T,sim_d7', title="Proj. & Pert", lab=lv)
#p8 = plot!(1:T,sim_d8', title="Proj. & Pert", lab=lv)
plot(p1,p2,p3,p4,p5,p6, size=1.25 .* (600, 400))
@JuliaRegistrator register()
First, thanks for this very nice package! I saw in the commit comment for v0.4 that changes were made to facilitate estimation, and that it is possible to solve the model passing in the parameters. This is something I would like to do, to provide students an example of estimation using simulated moments and related methods. I don't see documentation for how to do this - is any available, or planned? I can probably figure it out looking at the code, but if you're planning to add more details, I'm happy to wait a bit.
If you are interested, I'm planning to move the code for the model at page 451 of the document https://github.com/mcreel/Econometrics/blob/master/econometrics.pdf to Julia, using SolveDSGE.jl, instead of Dynare, This will help a lot to explain methods that need to do simulations at different parameter values.
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
@JuliaRegistrator register()
Hi,
I might have missed it but are there options to inspect the accuracy of a solution?
One common approach is to examine/plot residuals: for example (here)
@JuliaRegistrator register()
Would it be easier to
decision_rule(soln::R1) where {R1 <: SolutionType}
to src/analysis_functions.jl
soln
object in src/solution_functions.jl
soln_pl = solve_model(ngm,PL)
soln_pl.decision_rule(x)
The 2nd approach might involve structures.jl
struct FirstOrderSolutionStoch{T <: Real,S <: Integer} <: PerturbationSolutionStoch
# x(t+1) = hx*x(t) + k*v(t+1)
# y(t) = gx*x(t)
hbar::Union{T,Array{T,1}} # steady state values for predetermined variables
hx::Union{Array{T,2},Array{T,1}} # Transition matrix for predetermined variables
k::Union{Array{T,2},Array{T,1}} # Innovation loading matrix
gbar::Union{T,Array{T,1}} # steady state values for nonpredetermined variables
gx::Union{Array{T,2},Array{T,1}} # Decision rule matrix linking nonpredetermined variables to predetermined variables
sigma::Union{Array{T,2},Array{T,1}} # Innovation variance-covariance matrix
grc::S # Number of eigenvalues greater than cutoff
soln_type::String # "determinate", "indeterminate", or "explosive"
#
decision_rule::Function # policy function
end
julia> filename = "model3b.txt"
"model3b.txt"
julia> path3 = joinpath(@__DIR__,filename)
"/Users/AZevelev/Dropbox/COMPUTATION/Julia/DE_OC/SolveDSGE_ProjectionPerturbation/AZ_1/model3b.txt"
julia> process_model(path3)
The model's variables are now in this order: ["tech", "cap", "conslag", "cons"]
julia> processed_filename = "processed_model5a.txt"
"processed_model5a.txt"
julia> processed_path3 = joinpath(@__DIR__,processed_filename)
"/Users/AZevelev/Dropbox/COMPUTATION/Julia/DE_OC/SolveDSGE_ProjectionPerturbation/AZ_1/processed_model5a.txt"
julia> dsge3 = retrieve_processed_model(processed_path3)
ERROR: could not open file /Users/AZevelev/Dropbox/COMPUTATION/Julia/DE_OC/SolveDSGE_ProjectionPerturbation/AZ_1/processed_model5a_processed.txt
Stacktrace:
[1] include(::Function, ::Module, ::String) at ./Base.jl:380
[2] include at ./Base.jl:368 [inlined]
[3] include at /Users/AZevelev/.julia/packages/SolveDSGE/GMwhJ/src/SolveDSGE.jl:1 [inlined]
[4] retrieve_processed_model(::String) at /Users/AZevelev/.julia/packages/SolveDSGE/GMwhJ/src/parser_functions.jl:692
[5] top-level scope at none:1
@JuliaRegistrator register()
@JuliaRegistrator register()
The tag name "0.2.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @RJDennis
@JuliaRegistrator register()
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:
compute_steady_state()
gives the wrong steady state, so I used the closed form hereHi,
I was trying to execute the examples in my computer, but I got this error:
The model's variables are now in this order: ["z", "k", "c", "muc"]
ERROR: LoadError: StackOverflowError:
in expression starting at /Users/Github/SolveDSGE.jl/examples/example_one/solution_file_1a.jl:45
Something similar happens for the rest of examples, and I am not sure how to solve it.
Any help is appreciated.
I am confused by how this package pre-computes integrals for Chebyshev projection. The code snippet below pre-computes integrals in the package:
function compute_chebyshev_integrals(eps_nodes,eps_weights,order,sigma)
integrals = Array{Float64}(undef,order+1)
for i = 1:(order+1)
integrals[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
end
return integrals
end
My understanding is that this code calculates the Gauss-Hermite quadrature of exp(i * eps)
, where i
is the order of the polynomial approximation. The code then uses these integrals to re-scale the weights on each basis function, as shown below (Lines 733-744 of solution_functions.jl)
for i = 1:length(jumps_approximated)
weights[i] = chebyshev_weights(variables[jumps_approximated[i]],grid,order,domain)
end
scaled_weights = copy(weights)
for i = 1:length(jumps_approximated)
for j = 1:ns
index = [1:ndims(weights[i]);]
index[1],index[j] = index[j],index[1]
scaled_weights[i] = permutedims(integrals[j].*permutedims(weights[i],index),index)
end
end
I double checked the code for chebyshev_weights
, and these weights appear to be the coefficients on the order i
Chebyshev polynomial (please correct me if I'm wrong).
In the Judd, Maliar, and Maliar paper, these are the correct integrals to calculate when the approximating polynomials are ordinary polynomials. To precompute integrals for other polynomial families (and other function classes), the authors provide Proposition 1. By that proposition, I believe the correct integrals should actually be something like
function compute_chebyshev_integrals(eps_nodes,eps_weights,order,sigma,basis_functions)
integrals = Array{Float64}(undef,order+1)
for i = 1:(order+1)
f = basis_functions[i] # basis_functions is a vector of each basis function, which corresponds to the maximal order of the Chebyshev basis functions
integrals[i] = sum(f(sqrt(2)*sigma*eps_nodes).*eps_weights)*pi^(-1/2)
end
return integrals
end
Please let me know if I'm misunderstanding either your code or the paper.
Best,
William
@JuliaRegistrator register()
@JuliaRegistrator register()
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.