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
solvedsge.jl's People
Contributors
![femtocleaner[bot] avatar](https://avatars.githubusercontent.com/in/4123?v=4 "femtocleaner[bot]")
Stargazers
Watchers
Forkers
dingliumath mmkuang msaidf thomasbourany mostlyeconmix tkelman tranquilhero pgnepal binhpham79 jngod2011 wenddymacro amgalanb dawievlill pvfonseca shizelong1985 azev77 xinzhicaishu pbrazval allisterh marsteslasolvedsge.jl's Issues
Register v0.4.4
@JuliaRegistrator register()
Register version 0.4.14
@JuliaRegistrator register()
Register version 0.3.10
@JuliaRegistrator register()
Register version 0.4.0
@JuliaRegistrator register()
Register version 0.3.8
@JuliaRegistrator register()
Register version 0.4.14
@JuliaRegistrator register()
Error in julia 0.5
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
Register version 0.3.3
@JuliaRegistrator register()
Register version 0.3.7
@JuliaRegistrator register()
Register v0.4.7
@JuliaRegistrator register()
Register version 0.3.1
@JuliaRegistrator register()
occasionally binds
I can't help but ask, can this package solve the occasionally binds model?hh
Register v0.4.10
@JuliaRegistrator register()
using CUDA to alleviate high-dimensional problem in projection method
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
Example 3
Hi,
I am wondering if you could indicate a reference for the model of Example 3: Three shock New Keynesian Model, please?
Register version 0.3.2
@JuliaRegistrator register()
Register version 0.4.11
@JuliaRegistrator register()
Register v0.4.9
@JuliaRegistrator register()
Error when assigning parameters
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
Use this to verify steady stats
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)
`
Example with Finite Horizon
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?
Register v0.4.5
@JuliaRegistrator register()
Running solve_model inside a function
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
Register version 0.3.9
@JuliaRegistrator register()
Register v0.4.8
@JuliaRegistrator register()
Models w/o steady states: Consumption Saving Problems
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...
Register version 0.4.1
@JuliaRegistrator register()
nested task error: DomainError with -0.01... Try log(Complex(x)).
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]:1When 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]:1When 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]:1I've seen these complex argument errors a few times before using SolveDSGE.jl...
Grids setting
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?
Blanchard_kahn_form not defined
where i can find "Blanchard_kahn_form" ? when i ran, the code show it is not defined. thank you for your kind help.
Register version 0.3.12
@JuliaRegistrator register()
Firm investment: Another model w/o a steady state
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))Register version 0.3.11
@JuliaRegistrator register()
Question about estimation
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.
Register v0.4.3
@JuliaRegistrator register()
Register version 0.3.6
@JuliaRegistrator register()
Register version 0.4.13
@JuliaRegistrator register()
Register 0.3.9
@JuliaRegistrator register()
Feature Request: Options to Inspect Solution Accuracy
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)
Register version 0.3
@JuliaRegistrator register()
Policy function design
Would it be easier to
- add a function
decision_rule(soln::R1) where {R1 <: SolutionType}tosrc/analysis_functions.jl
I think this is the approach you're considering - include it as part of the
solnobject insrc/solution_functions.jl
In this way, for example:
soln_pl = solve_model(ngm,PL)
The policy function could besoln_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
endExamples 1-3 typo
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:1Register v0.4.6
@JuliaRegistrator register()
Register v0.4.2
@JuliaRegistrator register()
Error tagging new release
The tag name "0.2.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @RJDennis
Register version 0.4.12
@JuliaRegistrator register()
Useful example for the readme: comparison w/ closed form for deterministic NCGM
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 here- the projection methods don't seem to work for this problem, not sure why,,,
- At some point it would be great if we could automatically plot the policy fcns (then I can compare w/ sol from vfi etc)
ERROR: LoadError: StackOverflowError:
Hi,
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.
Pre-computation of integrals for Chebyshev projection
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
Register version 0.3.4
@JuliaRegistrator register()
Register version v0.3.5
@JuliaRegistrator register()
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.