From 9ca893ca1da526f0623f06df068b82e17060273a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Hermes=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Mon, 13 Jun 2016 11:58:43 +0200
Subject: [PATCH] Added routines for steady state computation (julia).
---
julia/DynareModel.jl | 3 +-
julia/DynareOptions.jl | 2 +-
julia/DynareOutput.jl | 5 +--
julia/SteadyState.jl | 96 ++++++++++++++++++++++++++++++++++++++++
preprocessor/ModFile.cc | 7 +--
tests/julia/rbc/test1.jl | 81 ++++++++++++++++++++++++++++++++-
tests/julia/rbc/test2.jl | 4 +-
7 files changed, 187 insertions(+), 11 deletions(-)
create mode 100644 julia/SteadyState.jl
diff --git a/julia/DynareModel.jl b/julia/DynareModel.jl
index a85442c9d2..f89cb9effd 100644
--- a/julia/DynareModel.jl
+++ b/julia/DynareModel.jl
@@ -18,8 +18,7 @@ module DynareModel
# along with Dynare. If not, see <http://www.gnu.org/licenses/>.
##
-
-export Endo, Exo, ExoDet, Param, dynare_model
+export Model, Endo, Exo, ExoDet, Param, dynare_model
abstract Atom
diff --git a/julia/DynareOptions.jl b/julia/DynareOptions.jl
index bdcca5c162..c5cc5e8881 100644
--- a/julia/DynareOptions.jl
+++ b/julia/DynareOptions.jl
@@ -19,7 +19,7 @@ module DynareOptions
##
-export dynare_options
+export Options, dynare_options
type Options
dynare_version::ASCIIString
diff --git a/julia/DynareOutput.jl b/julia/DynareOutput.jl
index 81a4132f75..b3ca48a23a 100644
--- a/julia/DynareOutput.jl
+++ b/julia/DynareOutput.jl
@@ -18,8 +18,7 @@ module DynareOutput
# along with Dynare. If not, see <http://www.gnu.org/licenses/>.
##
-
-export dynare_output
+export Ouput, dynare_output
type Output
dynare_version::ASCIIString
@@ -28,7 +27,7 @@ type Output
end
function dynare_output()
- return Output("", # dynare_version
+ return Output("", # dynare_version
Array(Float64, 0), # steady_state
Array(Float64, 0) # exo_steady_state
)
diff --git a/julia/SteadyState.jl b/julia/SteadyState.jl
new file mode 100644
index 0000000000..8d628563b6
--- /dev/null
+++ b/julia/SteadyState.jl
@@ -0,0 +1,96 @@
+module SteadyState
+
+##
+ # Copyright (C) 2016 Dynare Team
+ #
+ # This file is part of Dynare.
+ #
+ # Dynare is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU General Public License as published by
+ # the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # Dynare is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU General Public License
+ # along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+##
+
+using NLsolve
+
+import DynareModel.Model
+import DynareOutput.Output
+
+export steady, steady!
+export steady_state, steady_state!
+
+function steady(model::Model, oo::Output)
+ if model.analytical_steady_state || model.user_written_analytical_steady_state
+ steadystate = zeros(length(model.endo))
+ model.steady_state(steadystate, oo.exo_steady_state, model.params)
+ return steadystate
+ else
+ error("You have to provide a closed form solution for the steady state, or declare a guess\nfor the steady state as a third input argument.")
+ end
+end
+
+function steady!(model::Model, oo::Output)
+ if model.analytical_steady_state || model.user_written_analytical_steady_state
+ model.steady_state(oo.steady_state, oo.exo_steady_state, model.params)
+ return
+ else
+ error("You have to provide a closed form solution for the steady state, or declare a guess\nfor the steady state as a third input argument.")
+ end
+end
+
+function steady(model::Model, oo::Output, yinit::Vector{Float64})
+ ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
+ r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
+ if converged(r)
+ return r.zero
+ else
+ return fill(nan(Float64), length(yinit))
+ end
+end
+
+function steady!(model::Model, oo::Output, yinit::Vector{Float64})
+ ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
+ r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
+ if converged(r)
+ oo.steady_state = r.zero
+ else
+ oo.steady_state = fill(nan(Float64), length(yinit))
+ end
+end
+
+function steady_state(model::Model, oo::Output)
+ ys = steady(model, oo)
+ display_steady_state(model, oo, ys)
+end
+
+function steady_state!(model::Model, oo::Output)
+ steady!(model, oo)
+ display_steady_state(model, oo, oo.steady_state)
+end
+
+function display_steady_state(model::Model, oo::Output, ys::Vector{Float64})
+ println("\n\nSTEADY STATE:\n")
+ for i = 1:length(model.endo)
+ println(string(model.endo[i].name, " = ", ys[i]))
+ end
+end
+
+function issteadystate(model::Model, oo::Output, ys::Vector{Float64})
+ residuals = zeros(Float64, length(ys))
+ compute_static_model_residuals!(model, oo, ys, residuals)
+ return maximum(abs(residuals))<1e-6
+end
+
+function compute_static_model_residuals!(model::Model, oo::Output, ys::Vector{Float64}, residuals::Vector{Float64})
+ model.static(ys, oo.exo_steady_state, model.params, residuals)
+end
+
+end
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 2b2c684bdb..05c7bce25f 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -1104,11 +1104,12 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
<< "#" << endl
<< "# NB: this file was automatically generated by Dynare" << endl
<< "# from " << basename << ".mod" << endl
- << "#" << endl
+ << "#" << endl << endl
<< "using DynareModel" << endl
<< "using DynareOptions" << endl
- << "using DynareOutput" << endl
+ << "using DynareOutput" << endl << endl
<< "using Utils" << endl
+ << "using SteadyState" << endl << endl
<< "using " << basename << "Static" << endl
<< "using " << basename << "Dynamic" << endl
<< "if isfile(\"" << basename << "SteadyState.jl" "\")" << endl
@@ -1197,7 +1198,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
<< " using " << basename << "DynamicParamsDerivs" << endl
<< " model_.dynamic_params_derivs = " << basename << "DynamicParamsDerivs.params_derivs" << endl
<< "end" << endl
- << "end" << endl;
+ << "end" << endl;
jlOutputFile.close();
cout << "done" << endl;
}
diff --git a/tests/julia/rbc/test1.jl b/tests/julia/rbc/test1.jl
index 8f7ba07ca5..c4579fa377 100644
--- a/tests/julia/rbc/test1.jl
+++ b/tests/julia/rbc/test1.jl
@@ -10,4 +10,83 @@ importall Dynare
@dynare "rbc1.mod"
-print(model_.fname)
+importall SteadyState
+
+# First call to the steady state routine (analytical)
+@time SteadyState.steady!(model_, oo_)
+
+# First call to the steady state routine (analytical)
+@time SteadyState.steady!(model_, oo_)
+
+paramsinit = copy(model_.params);
+
+yinit = deepcopy(oo_.steady_state)
+
+y_init = copy(yinit)
+y_init[1] = 1.1*yinit[1]
+y_init[2] = 0.9*yinit[2]
+
+# First call to the steady state routine (numerical)
+println("First call to the numerical steady state routine")
+@time SteadyState.steady!(model_, oo_, yinit)
+
+# Check results
+@assert maximum(abs(oo_.steady_state-yinit))<1e-9
+
+y_init = copy(yinit)
+yinit[1] = 1.1*yinit[1]
+yinit[2] = 0.9*yinit[2]
+
+# Second call to the steady state routine (numerical)
+println("Second call to the numerical steady state routine")
+@time SteadyState.steady!(model_, oo_, yinit)
+
+# change alpha
+println("Change α")
+model_.params[4] = max(min(1.0, model_.params[4]*1.2), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change delta
+println("Change δ")
+model_.params[6] = max(min(1.0, model_.params[6]*1.2), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change beta
+println("Change β")
+model_.params[1] = max(min(1-1e-6, model_.params[1]*0.99), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change tau
+println("Change τ")
+model_.params[3] /= 1.5
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change Epsilon
+println("Change ϵ")
+model_.params[5] *= 1.5
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
diff --git a/tests/julia/rbc/test2.jl b/tests/julia/rbc/test2.jl
index aa7a9566bd..f9eae9ebfa 100644
--- a/tests/julia/rbc/test2.jl
+++ b/tests/julia/rbc/test2.jl
@@ -11,4 +11,6 @@ importall Dynare
@dynare "rbc2.mod"
-print(model_.fname)
+importall SteadyState
+
+steady!(model_, oo_)
--
GitLab