diff --git a/julia/v0.5/Dynare.jl b/julia/v0.5/Dynare.jl new file mode 100644 index 0000000000000000000000000000000000000000..82efd64f2d672b17931f8b69fa72cdfeefa732bd --- /dev/null +++ b/julia/v0.5/Dynare.jl @@ -0,0 +1,57 @@ +module Dynare + +## + # Copyright (C) 2015-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/>. +## + +export @compile, @dynare + +function compile(modfile) + # Add cd to path if not already there + if isempty(findin([pwd()], LOAD_PATH)) + unshift!(LOAD_PATH, pwd()) + end + # Process modfile + println(string("Using ", WORD_SIZE, "-bit preprocessor")) + preprocessor = string(dirname(@__FILE__()), "/preprocessor", WORD_SIZE, "/dynare_m") + run(`$preprocessor $modfile language=julia output=dynamic`) +end + +macro dynare(modfiles...) + ex = Expr(:toplevel) + if length(modfiles)>1 + for modfile in modfiles + eval(:(compile($modfile))) + basename = split(modfile, ".mod"; keep=false) + push!(ex.args, Expr(:import, symbol(basename[1]))) + end + else + eval(:(compile($modfiles))) + basename = split(modfiles[1], ".mod"; keep=false) + push!(ex.args, Expr(:importall, symbol(basename[1]))) + end + return ex +end + +macro compile(modfiles...) + for modfile in modfiles + eval(:(compile($modfile))) + end +end + +end diff --git a/julia/v0.5/DynareModel.jl b/julia/v0.5/DynareModel.jl new file mode 100644 index 0000000000000000000000000000000000000000..6c4b6110e938a654fcd3b2740dc712db8b81775e --- /dev/null +++ b/julia/v0.5/DynareModel.jl @@ -0,0 +1,181 @@ +module DynareModel +## + # Copyright (C) 2015 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/>. +## + +export Model, Endo, Exo, ExoDet, Param, dynare_model + +abstract Atom + +immutable Endo <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable Exo <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable ExoDet <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable Param <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable AuxVars + endo_index::Int + var_type::Int + orig_index::Int + orig_lead_lag::Int + eq_nbr::Int + orig_expr::UTF8String +end + +immutable PredVars + index::Int +end + +immutable ObsVars + index::Int +end + +immutable DetShocks + exo_det::Int + exo_id::Int + multiplicative::Bool + periods::Vector{Int} + value::Float64 +end + +immutable EquationTag + eq_nbr::Int + name::UTF8String + value::UTF8String +end + +type Model + fname::ASCIIString + dname::ASCIIString + dynare_version::ASCIIString + endo::Vector{Endo} + exo::Vector{Exo} + exo_det::Vector{ExoDet} + param::Vector{Param} + aux_vars::Vector{AuxVars} + pred_vars::Vector{Int} + obs_vars::Vector{Int} + orig_endo_nbr::Int + orig_eq_nbr::Int + eq_nbr::Int + ramsey_eq_nbr::Int + det_shocks::Vector{DetShocks} + nstatic::Int + nfwrd::Int + npred::Int + nboth::Int + nsfwrd::Int + nspred::Int + ndynamic::Int + maximum_lag::Int + maximum_lead::Int + maximum_endo_lag::Int + maximum_endo_lead::Int + maximum_exo_lag::Int + maximum_exo_lead::Int + lead_lag_incidence::Matrix{Int} + nnzderivatives::Vector{Int} + analytical_steady_state::Bool + user_written_analytical_steady_state::Bool + static_and_dynamic_models_differ::Bool + equation_tags::Vector{UTF8String} + exo_names_orig_ord::Vector{Int} + sigma_e::Matrix{Float64} + correlation_matrix::Matrix{Float64} + h::Matrix{Float64} + correlation_matrix_me::Matrix{Float64} + sigma_e_is_diagonal::Bool + params::Vector{Float64} + static::Function + static_params_derivs::Function + dynamic::Function + dynamic_params_derivs::Function + first_derivatives::Function + steady_state::Function +end + +function dynare_model() + return Model("", # fname + "", # dname + "", # dynare_version + Array(Endo,0), # endo + Array(Exo,0), # exo + Array(ExoDet,0), # exo_det + Array(Param,0), # param + Array(AuxVars,0), # aux_vars + Array(Int,0), # pred_vars + Array(Int,0), # obs_vars + 0, # orig_endo_nbr + 0, # orig_eq_nbr + 0, # eq_nbr + 0, # ramsey_eq_nbr + Array(DetShocks,0), # det_shocks + 0, # nstatic + 0, # nfwrd + 0, # npred + 0, # nboth + 0, # nsfwrd + 0, # nspred + 0, # ndynamic + 0, # maximum_lag + 0, # maximum_lead + 0, # maximum_endo_lag + 0, # maximum_endo_lead + 0, # maximum_exo_lag + 0, # maximum_exo_lead + Array(Int, 3, 0), # lead_lag_incidence + zeros(Int, 3), # nnzderivatives + false, # analytical_steady_state + false, # user_written_analytical_steady_state + false, # static_and_dynamic_models_differ + Array(ASCIIString,0), # equation_tags + Array(Int64,1), # exo_names_orig_ord + Array(Float64, 0, 0), # sigma_e (Cov matrix of the structural innovations) + Array(Float64, 0, 0), # correlation_matrix (Corr matrix of the structural innovations) + Array(Float64, 0, 0), # h (Cov matrix of the measurement errors) + Array(Float64, 0, 0), # correlation_matrix_me (Cov matrix of the measurement errors) + true, # sigma_e_is_diagonal + Array(Float64, 0), # params + function()end, # static + function()end, # static_params_derivs + function()end, # dynamic + function()end, # dynamic_params_derivs + function()end, # first_derivatives + function()end # steady_state + ) +end + +end diff --git a/julia/v0.5/DynareOptions.jl b/julia/v0.5/DynareOptions.jl new file mode 100644 index 0000000000000000000000000000000000000000..b25bdc27cd0e09ce9a6da651d3aaac6ddadb5871 --- /dev/null +++ b/julia/v0.5/DynareOptions.jl @@ -0,0 +1,51 @@ +module DynareOptions +## + # Copyright (C) 2015 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/>. +## + +export Options, dynare_options + +type PFMSolver + maxit::Int + periods::Int + tolx::Float64 + tolf::Float64 +end + +function pfmsolver_set_defaults() + return PFMSolver(500, # maxit (Maximum number of iterations in Newton algorithm) + 400, # periods (Number of periods to return to the steady state) + 1e-6, # tolx (Tolerance criterion on the paths for the endogenous variables) + 1e-6 # tolf (Tolerance criterion on the stacked non linear equations) + ) +end + +type Options + dynare_version::ASCIIString + linear::Bool + pfmsolver::PFMSolver +end + +function dynare_options() + return Options("", # dynare_version + false, # linear + pfmsolver_set_defaults() # pfmsolver + ) +end + +end diff --git a/julia/v0.5/DynareOutput.jl b/julia/v0.5/DynareOutput.jl new file mode 100644 index 0000000000000000000000000000000000000000..b3ca48a23a3ea399b792d37b881774376375f8e9 --- /dev/null +++ b/julia/v0.5/DynareOutput.jl @@ -0,0 +1,36 @@ +module DynareOutput +## + # Copyright (C) 2015 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/>. +## + +export Ouput, dynare_output + +type Output + dynare_version::ASCIIString + steady_state::Vector{Float64} + exo_steady_state::Vector{Float64} +end + +function dynare_output() + return Output("", # dynare_version + Array(Float64, 0), # steady_state + Array(Float64, 0) # exo_steady_state + ) +end + +end diff --git a/julia/v0.5/PerfectForesightModelSolver.jl b/julia/v0.5/PerfectForesightModelSolver.jl new file mode 100644 index 0000000000000000000000000000000000000000..6a4cf042f2556545e02954910fdb8fca8c56fa01 --- /dev/null +++ b/julia/v0.5/PerfectForesightModelSolver.jl @@ -0,0 +1,134 @@ +module PerfectForesightModelSolver + +## + # 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/>. +## + +import DynareModel.Model +import DynareOutput.Output +import DynareOptions.Options + +export simulate_perfect_foresight_model! + +function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64}, exogenousvariables::Matrix{Float64}, steadystate::Vector{Float64}, model::Model, options::Options) + + lead_lag_incidence = model.lead_lag_incidence + + nyp = countnz(lead_lag_incidence[1,:]) + ny0 = countnz(lead_lag_incidence[2,:]) + nyf = countnz(lead_lag_incidence[3,:]) + + ny = length(model.endo) + nd = nyp+ny0+nyf + + periods = options.pfmsolver.periods + params = model.params + + tmp = lead_lag_incidence[2:3,:]' + i_cols_A1 = find(tmp) + i_cols_1 = tmp[i_cols_A1] + + tmp = lead_lag_incidence[1:2,:]' + i_cols_AT = find(tmp) + i_cols_T = tmp[i_cols_AT] + + tmp = lead_lag_incidence[2,:]' + i_cols_A0 = find(tmp) + i_cols_0 = tmp[i_cols_A0] + + i_cols_j = collect(1:nd) + i_upd = ny+collect(1:periods*ny) + + Y = vec(endogenousvariables) + z = Y[find(lead_lag_incidence')] + + jacobian = zeros(Float64, ny, length(z)+length(model.exo)) + residuals = zeros(Float64, ny) + + println("\nMODEL SIMULATION:\n") + + rd = zeros(Float64, periods*ny) + + iA = zeros(Int64, periods*model.nnzderivatives[1]) + jA = zeros(Int64, periods*model.nnzderivatives[1]) + vA = zeros(Float64, periods*model.nnzderivatives[1]) + + convergence = false + iteration = 0 + + while !convergence + iteration += 1 + i_rows = collect(1:ny) + i_cols_A = find(lead_lag_incidence') + i_cols = i_cols_A + m = 0 + for it = 2:(periods+1) + model.dynamic(Y[i_cols], exogenousvariables, params, steadystate, it, residuals, jacobian) + if it==(periods+1) & it==2 + (r, c, v) = findnz(jacobian[:,i_cols_0]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A0[c] + vA[k] = v + elseif it==(periods+1) + (r, c, v) = findnz(jacobian[:,i_cols_T]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A[i_cols_T[c]] + vA[k] = v + elseif it==2 + (r, c, v) = findnz(jacobian[:,i_cols_1]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A1[c] + vA[k] = v + else + (r, c, v) = findnz(jacobian[:,i_cols_j]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A[c] + vA[k] = v + end + m += length(v) + rd[i_rows] = residuals + i_rows += ny + i_cols += ny + if it>2 + i_cols_A += ny + end + end + err = maximum(abs(rd)) + println("Iter. ", iteration, "\t err. ", round(err, 12)) + if err<options.pfmsolver.tolf + iteration -= 1 + convergence = true + end + A = sparse(iA[1:m], jA[1:m], vA[1:m]) + @time dy = A\rd + Y[i_upd] -= dy + if maximum(abs(dy))<options.pfmsolver.tolx + convergence = true + end + end + if convergence + println("\nPFM solver converged in ", iteration, " iterations!\n") + endogenousvariables = reshape(Y, ny, periods+2) + end +end + +end diff --git a/julia/v0.5/PerfectForesightModelSolver2.jl b/julia/v0.5/PerfectForesightModelSolver2.jl new file mode 100644 index 0000000000000000000000000000000000000000..46929c5fc3d2fc0cdeab2d67a069f97341425b0b --- /dev/null +++ b/julia/v0.5/PerfectForesightModelSolver2.jl @@ -0,0 +1,185 @@ +module PerfectForesightModelSolver2 + +## + # 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/>. +## + +import DynareModel.Model +import DynareOutput.Output +import DynareOptions.Options + +export simulate_perfect_foresight_model! + +function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64}, exogenousvariables::Matrix{Float64}, steadystate::Vector{Float64}, model::Model, options::Options) + + lead_lag_incidence = model.lead_lag_incidence + + nyp = countnz(lead_lag_incidence[1,:]) + ny0 = countnz(lead_lag_incidence[2,:]) + nyf = countnz(lead_lag_incidence[3,:]) + + ny = length(model.endo) + nd = nyp+ny0+nyf + + periods = options.pfmsolver.periods + params = model.params + + tmp = lead_lag_incidence[2:3,:]' + i_cols_A1 = find(tmp) + i_cols_1 = tmp[i_cols_A1] + + tmp = lead_lag_incidence[1:2,:]' + i_cols_AT = find(tmp) + i_cols_T = tmp[i_cols_AT] + + tmp = lead_lag_incidence[2,:]' + i_cols_A0 = find(tmp) + i_cols_0 = tmp[i_cols_A0] + + i_cols_j = collect(1:nd) + i_upd = ny+collect(1:periods*ny) + + Y = vec(endogenousvariables) + z = Y[find(lead_lag_incidence')] + + jacobian = zeros(Float64, ny, length(z)+length(model.exo)) + residuals = zeros(Float64, ny) + + println("\nMODEL SIMULATION:\n") + + rd = zeros(Float64, periods*ny) + A = spzeros(periods*ny,periods*ny) + convergence = false + iteration = 0 + + ws = PerfectForesightModelWS(periods,model) + + while !convergence + iteration += 1 + perfect_foresight_model!(Y,exogenousvariables,model,steadystate,ws,rd) + err = maximum(abs(rd)) + println("Iter. ", iteration, "\t err. ", round(err, 12)) + if err<options.pfmsolver.tolf + iteration -= 1 + convergence = true + break + end + A = perfect_foresight_model!(Y,exogenousvariables,model,steadystate,ws,A) + @time dy = A\rd + Y[i_upd] -= dy + if maximum(abs(dy))<options.pfmsolver.tolx + convergence = true + end + end + if convergence + println("\nPFM solver converged in ", iteration, " iterations!\n") + endogenousvariables = reshape(Y, ny, periods+2) + end +end + +type PerfectForesightModelWS + periods::Int64 + iA::Array{Int64,1} + jA::Array{Int64,1} + vA::Array{Float64,1} + iP::Array{Int64,1} + jP::Array{Int64,1} + vP::Array{Float64,1} + function PerfectForesightModelWS(n,model) + periods = n + nzd = model.nnzderivatives[1] + iA = zeros(Int64,n*nzd) + jA = zeros(Int64,n*nzd) + vA = zeros(Float64,n*nzd) + iP = zeros(Int64,nzd) + jP = zeros(Int64,nzd) + vP = zeros(Float64,nzd) + new(periods,iA,jA,vA,iP,jP,vP) + end +end + +function perfect_foresight_model!(Y::Array{Float64,1},exogenousvariables::Array{Float64,2},model::Model,steadystate::Array{Float64,1},ws::PerfectForesightModelWS,residuals::Array{Float64,1}) + periods = ws.periods + ny = length(model.endo) + i_rows = collect(1:ny) + i_cols = find(model.lead_lag_incidence') + m = 0 + for it = 2:(periods+1) + res = sub(residuals,m+1:m+model.eq_nbr) + Yview = sub(Y,i_cols) + model.dynamic(Yview, exogenousvariables, model.params, steadystate, it, res) + m += model.eq_nbr + i_cols += ny + end +end + +function perfect_foresight_model!(Y::Array{Float64,1},exogenousvariables::Array{Float64,2},model::Model,steadystate::Array{Float64,1},ws::PerfectForesightModelWS,A::SparseMatrixCSC{Float64,Int64}) + periods = ws.periods + ny = length(model.endo) + i_rows = collect(1:ny) + i_cols = find(model.lead_lag_incidence') + m = 0 + offset_r = 0 + offset_c = -ny + nzd = 0 + for it = 2:(periods+1) + Yview = sub(Y,i_cols) + if it == 2 + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, ws.iP, ws.jP, ws.vP) + nzd = count(i->(0 .< i .<= 3*ny),ws.jP) + k1 = 1 + for k = 1:nzd + if ws.jP[k] > ny + ws.iA[k1] = ws.iP[k] + ws.jA[k1] = ws.jP[k] + offset_c + ws.vA[k1] = ws.vP[k] + k1 += 1 + end + end + m = k1 - 1 + elseif it==(periods+1) + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, ws.iP, ws.jP, ws.vP) + for k=1:nzd + if ws.jP[k] <= 2*ny + m += 1 + ws.iA[m] = ws.iP[k] + offset_r + ws.jA[m] = ws.jP[k] + offset_c + ws.vA[m] = ws.vP[k] + end + end + else + I = sub(ws.iA,m+1:m+nzd) + J = sub(ws.jA,m+1:m+nzd) + V = sub(ws.vA,m+1:m+nzd) + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, I,J,V) + for k=m+1:m+nzd + ws.iA[k] += offset_r + ws.jA[k] += offset_c + end + m += nzd + end + offset_r += ny + offset_c += ny + i_cols += ny + end + A = sparse(ws.iA[1:m], ws.jA[1:m], ws.vA[1:m]) + assert(size(A)==(ny*periods,ny*periods)) + return A +end + +end diff --git a/julia/v0.5/SteadyState.jl b/julia/v0.5/SteadyState.jl new file mode 100644 index 0000000000000000000000000000000000000000..8d628563b6074ebbe5079189176e2bdd69ee567c --- /dev/null +++ b/julia/v0.5/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/julia/v0.5/Utils.jl b/julia/v0.5/Utils.jl new file mode 100644 index 0000000000000000000000000000000000000000..6fd7eff5316a3f97410331793983eed1a5404f6b --- /dev/null +++ b/julia/v0.5/Utils.jl @@ -0,0 +1,36 @@ +module Utils +## + # Copyright (C) 2015 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/>. +## + +export get_power_deriv + +function get_power_deriv(x::Float64, p::Real, k::Int) + if abs(x)<1e-12 && p>0 && k>p && typeof(p)==Int + dxp = .0 + else + dxp = x^(p-k) + for i = 0:k-1 + dxp *= p + p -= 1 + end + end + return dxp +end + +end diff --git a/julia/v0.5/preprocessor64/dynare_m b/julia/v0.5/preprocessor64/dynare_m new file mode 120000 index 0000000000000000000000000000000000000000..9d36d77010cefc22462bc097d09c06cfeebf3fdc --- /dev/null +++ b/julia/v0.5/preprocessor64/dynare_m @@ -0,0 +1 @@ +/home/michel/dynare/git/master/preprocessor/dynare_m \ No newline at end of file diff --git a/julia/v0.5/v0.5/Dynare.jl b/julia/v0.5/v0.5/Dynare.jl new file mode 100644 index 0000000000000000000000000000000000000000..82efd64f2d672b17931f8b69fa72cdfeefa732bd --- /dev/null +++ b/julia/v0.5/v0.5/Dynare.jl @@ -0,0 +1,57 @@ +module Dynare + +## + # Copyright (C) 2015-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/>. +## + +export @compile, @dynare + +function compile(modfile) + # Add cd to path if not already there + if isempty(findin([pwd()], LOAD_PATH)) + unshift!(LOAD_PATH, pwd()) + end + # Process modfile + println(string("Using ", WORD_SIZE, "-bit preprocessor")) + preprocessor = string(dirname(@__FILE__()), "/preprocessor", WORD_SIZE, "/dynare_m") + run(`$preprocessor $modfile language=julia output=dynamic`) +end + +macro dynare(modfiles...) + ex = Expr(:toplevel) + if length(modfiles)>1 + for modfile in modfiles + eval(:(compile($modfile))) + basename = split(modfile, ".mod"; keep=false) + push!(ex.args, Expr(:import, symbol(basename[1]))) + end + else + eval(:(compile($modfiles))) + basename = split(modfiles[1], ".mod"; keep=false) + push!(ex.args, Expr(:importall, symbol(basename[1]))) + end + return ex +end + +macro compile(modfiles...) + for modfile in modfiles + eval(:(compile($modfile))) + end +end + +end diff --git a/julia/v0.5/v0.5/DynareModel.jl b/julia/v0.5/v0.5/DynareModel.jl new file mode 100644 index 0000000000000000000000000000000000000000..6c4b6110e938a654fcd3b2740dc712db8b81775e --- /dev/null +++ b/julia/v0.5/v0.5/DynareModel.jl @@ -0,0 +1,181 @@ +module DynareModel +## + # Copyright (C) 2015 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/>. +## + +export Model, Endo, Exo, ExoDet, Param, dynare_model + +abstract Atom + +immutable Endo <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable Exo <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable ExoDet <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable Param <: Atom + name::UTF8String + tex_name::UTF8String + long_name::UTF8String +end + +immutable AuxVars + endo_index::Int + var_type::Int + orig_index::Int + orig_lead_lag::Int + eq_nbr::Int + orig_expr::UTF8String +end + +immutable PredVars + index::Int +end + +immutable ObsVars + index::Int +end + +immutable DetShocks + exo_det::Int + exo_id::Int + multiplicative::Bool + periods::Vector{Int} + value::Float64 +end + +immutable EquationTag + eq_nbr::Int + name::UTF8String + value::UTF8String +end + +type Model + fname::ASCIIString + dname::ASCIIString + dynare_version::ASCIIString + endo::Vector{Endo} + exo::Vector{Exo} + exo_det::Vector{ExoDet} + param::Vector{Param} + aux_vars::Vector{AuxVars} + pred_vars::Vector{Int} + obs_vars::Vector{Int} + orig_endo_nbr::Int + orig_eq_nbr::Int + eq_nbr::Int + ramsey_eq_nbr::Int + det_shocks::Vector{DetShocks} + nstatic::Int + nfwrd::Int + npred::Int + nboth::Int + nsfwrd::Int + nspred::Int + ndynamic::Int + maximum_lag::Int + maximum_lead::Int + maximum_endo_lag::Int + maximum_endo_lead::Int + maximum_exo_lag::Int + maximum_exo_lead::Int + lead_lag_incidence::Matrix{Int} + nnzderivatives::Vector{Int} + analytical_steady_state::Bool + user_written_analytical_steady_state::Bool + static_and_dynamic_models_differ::Bool + equation_tags::Vector{UTF8String} + exo_names_orig_ord::Vector{Int} + sigma_e::Matrix{Float64} + correlation_matrix::Matrix{Float64} + h::Matrix{Float64} + correlation_matrix_me::Matrix{Float64} + sigma_e_is_diagonal::Bool + params::Vector{Float64} + static::Function + static_params_derivs::Function + dynamic::Function + dynamic_params_derivs::Function + first_derivatives::Function + steady_state::Function +end + +function dynare_model() + return Model("", # fname + "", # dname + "", # dynare_version + Array(Endo,0), # endo + Array(Exo,0), # exo + Array(ExoDet,0), # exo_det + Array(Param,0), # param + Array(AuxVars,0), # aux_vars + Array(Int,0), # pred_vars + Array(Int,0), # obs_vars + 0, # orig_endo_nbr + 0, # orig_eq_nbr + 0, # eq_nbr + 0, # ramsey_eq_nbr + Array(DetShocks,0), # det_shocks + 0, # nstatic + 0, # nfwrd + 0, # npred + 0, # nboth + 0, # nsfwrd + 0, # nspred + 0, # ndynamic + 0, # maximum_lag + 0, # maximum_lead + 0, # maximum_endo_lag + 0, # maximum_endo_lead + 0, # maximum_exo_lag + 0, # maximum_exo_lead + Array(Int, 3, 0), # lead_lag_incidence + zeros(Int, 3), # nnzderivatives + false, # analytical_steady_state + false, # user_written_analytical_steady_state + false, # static_and_dynamic_models_differ + Array(ASCIIString,0), # equation_tags + Array(Int64,1), # exo_names_orig_ord + Array(Float64, 0, 0), # sigma_e (Cov matrix of the structural innovations) + Array(Float64, 0, 0), # correlation_matrix (Corr matrix of the structural innovations) + Array(Float64, 0, 0), # h (Cov matrix of the measurement errors) + Array(Float64, 0, 0), # correlation_matrix_me (Cov matrix of the measurement errors) + true, # sigma_e_is_diagonal + Array(Float64, 0), # params + function()end, # static + function()end, # static_params_derivs + function()end, # dynamic + function()end, # dynamic_params_derivs + function()end, # first_derivatives + function()end # steady_state + ) +end + +end diff --git a/julia/v0.5/v0.5/DynareOptions.jl b/julia/v0.5/v0.5/DynareOptions.jl new file mode 100644 index 0000000000000000000000000000000000000000..b25bdc27cd0e09ce9a6da651d3aaac6ddadb5871 --- /dev/null +++ b/julia/v0.5/v0.5/DynareOptions.jl @@ -0,0 +1,51 @@ +module DynareOptions +## + # Copyright (C) 2015 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/>. +## + +export Options, dynare_options + +type PFMSolver + maxit::Int + periods::Int + tolx::Float64 + tolf::Float64 +end + +function pfmsolver_set_defaults() + return PFMSolver(500, # maxit (Maximum number of iterations in Newton algorithm) + 400, # periods (Number of periods to return to the steady state) + 1e-6, # tolx (Tolerance criterion on the paths for the endogenous variables) + 1e-6 # tolf (Tolerance criterion on the stacked non linear equations) + ) +end + +type Options + dynare_version::ASCIIString + linear::Bool + pfmsolver::PFMSolver +end + +function dynare_options() + return Options("", # dynare_version + false, # linear + pfmsolver_set_defaults() # pfmsolver + ) +end + +end diff --git a/julia/v0.5/v0.5/DynareOutput.jl b/julia/v0.5/v0.5/DynareOutput.jl new file mode 100644 index 0000000000000000000000000000000000000000..b3ca48a23a3ea399b792d37b881774376375f8e9 --- /dev/null +++ b/julia/v0.5/v0.5/DynareOutput.jl @@ -0,0 +1,36 @@ +module DynareOutput +## + # Copyright (C) 2015 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/>. +## + +export Ouput, dynare_output + +type Output + dynare_version::ASCIIString + steady_state::Vector{Float64} + exo_steady_state::Vector{Float64} +end + +function dynare_output() + return Output("", # dynare_version + Array(Float64, 0), # steady_state + Array(Float64, 0) # exo_steady_state + ) +end + +end diff --git a/julia/v0.5/v0.5/PerfectForesightModelSolver.jl b/julia/v0.5/v0.5/PerfectForesightModelSolver.jl new file mode 100644 index 0000000000000000000000000000000000000000..6a4cf042f2556545e02954910fdb8fca8c56fa01 --- /dev/null +++ b/julia/v0.5/v0.5/PerfectForesightModelSolver.jl @@ -0,0 +1,134 @@ +module PerfectForesightModelSolver + +## + # 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/>. +## + +import DynareModel.Model +import DynareOutput.Output +import DynareOptions.Options + +export simulate_perfect_foresight_model! + +function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64}, exogenousvariables::Matrix{Float64}, steadystate::Vector{Float64}, model::Model, options::Options) + + lead_lag_incidence = model.lead_lag_incidence + + nyp = countnz(lead_lag_incidence[1,:]) + ny0 = countnz(lead_lag_incidence[2,:]) + nyf = countnz(lead_lag_incidence[3,:]) + + ny = length(model.endo) + nd = nyp+ny0+nyf + + periods = options.pfmsolver.periods + params = model.params + + tmp = lead_lag_incidence[2:3,:]' + i_cols_A1 = find(tmp) + i_cols_1 = tmp[i_cols_A1] + + tmp = lead_lag_incidence[1:2,:]' + i_cols_AT = find(tmp) + i_cols_T = tmp[i_cols_AT] + + tmp = lead_lag_incidence[2,:]' + i_cols_A0 = find(tmp) + i_cols_0 = tmp[i_cols_A0] + + i_cols_j = collect(1:nd) + i_upd = ny+collect(1:periods*ny) + + Y = vec(endogenousvariables) + z = Y[find(lead_lag_incidence')] + + jacobian = zeros(Float64, ny, length(z)+length(model.exo)) + residuals = zeros(Float64, ny) + + println("\nMODEL SIMULATION:\n") + + rd = zeros(Float64, periods*ny) + + iA = zeros(Int64, periods*model.nnzderivatives[1]) + jA = zeros(Int64, periods*model.nnzderivatives[1]) + vA = zeros(Float64, periods*model.nnzderivatives[1]) + + convergence = false + iteration = 0 + + while !convergence + iteration += 1 + i_rows = collect(1:ny) + i_cols_A = find(lead_lag_incidence') + i_cols = i_cols_A + m = 0 + for it = 2:(periods+1) + model.dynamic(Y[i_cols], exogenousvariables, params, steadystate, it, residuals, jacobian) + if it==(periods+1) & it==2 + (r, c, v) = findnz(jacobian[:,i_cols_0]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A0[c] + vA[k] = v + elseif it==(periods+1) + (r, c, v) = findnz(jacobian[:,i_cols_T]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A[i_cols_T[c]] + vA[k] = v + elseif it==2 + (r, c, v) = findnz(jacobian[:,i_cols_1]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A1[c] + vA[k] = v + else + (r, c, v) = findnz(jacobian[:,i_cols_j]) + k = collect(1:length(v))+m + iA[k] = i_rows[r] + jA[k] = i_cols_A[c] + vA[k] = v + end + m += length(v) + rd[i_rows] = residuals + i_rows += ny + i_cols += ny + if it>2 + i_cols_A += ny + end + end + err = maximum(abs(rd)) + println("Iter. ", iteration, "\t err. ", round(err, 12)) + if err<options.pfmsolver.tolf + iteration -= 1 + convergence = true + end + A = sparse(iA[1:m], jA[1:m], vA[1:m]) + @time dy = A\rd + Y[i_upd] -= dy + if maximum(abs(dy))<options.pfmsolver.tolx + convergence = true + end + end + if convergence + println("\nPFM solver converged in ", iteration, " iterations!\n") + endogenousvariables = reshape(Y, ny, periods+2) + end +end + +end diff --git a/julia/v0.5/v0.5/PerfectForesightModelSolver2.jl b/julia/v0.5/v0.5/PerfectForesightModelSolver2.jl new file mode 100644 index 0000000000000000000000000000000000000000..46929c5fc3d2fc0cdeab2d67a069f97341425b0b --- /dev/null +++ b/julia/v0.5/v0.5/PerfectForesightModelSolver2.jl @@ -0,0 +1,185 @@ +module PerfectForesightModelSolver2 + +## + # 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/>. +## + +import DynareModel.Model +import DynareOutput.Output +import DynareOptions.Options + +export simulate_perfect_foresight_model! + +function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64}, exogenousvariables::Matrix{Float64}, steadystate::Vector{Float64}, model::Model, options::Options) + + lead_lag_incidence = model.lead_lag_incidence + + nyp = countnz(lead_lag_incidence[1,:]) + ny0 = countnz(lead_lag_incidence[2,:]) + nyf = countnz(lead_lag_incidence[3,:]) + + ny = length(model.endo) + nd = nyp+ny0+nyf + + periods = options.pfmsolver.periods + params = model.params + + tmp = lead_lag_incidence[2:3,:]' + i_cols_A1 = find(tmp) + i_cols_1 = tmp[i_cols_A1] + + tmp = lead_lag_incidence[1:2,:]' + i_cols_AT = find(tmp) + i_cols_T = tmp[i_cols_AT] + + tmp = lead_lag_incidence[2,:]' + i_cols_A0 = find(tmp) + i_cols_0 = tmp[i_cols_A0] + + i_cols_j = collect(1:nd) + i_upd = ny+collect(1:periods*ny) + + Y = vec(endogenousvariables) + z = Y[find(lead_lag_incidence')] + + jacobian = zeros(Float64, ny, length(z)+length(model.exo)) + residuals = zeros(Float64, ny) + + println("\nMODEL SIMULATION:\n") + + rd = zeros(Float64, periods*ny) + A = spzeros(periods*ny,periods*ny) + convergence = false + iteration = 0 + + ws = PerfectForesightModelWS(periods,model) + + while !convergence + iteration += 1 + perfect_foresight_model!(Y,exogenousvariables,model,steadystate,ws,rd) + err = maximum(abs(rd)) + println("Iter. ", iteration, "\t err. ", round(err, 12)) + if err<options.pfmsolver.tolf + iteration -= 1 + convergence = true + break + end + A = perfect_foresight_model!(Y,exogenousvariables,model,steadystate,ws,A) + @time dy = A\rd + Y[i_upd] -= dy + if maximum(abs(dy))<options.pfmsolver.tolx + convergence = true + end + end + if convergence + println("\nPFM solver converged in ", iteration, " iterations!\n") + endogenousvariables = reshape(Y, ny, periods+2) + end +end + +type PerfectForesightModelWS + periods::Int64 + iA::Array{Int64,1} + jA::Array{Int64,1} + vA::Array{Float64,1} + iP::Array{Int64,1} + jP::Array{Int64,1} + vP::Array{Float64,1} + function PerfectForesightModelWS(n,model) + periods = n + nzd = model.nnzderivatives[1] + iA = zeros(Int64,n*nzd) + jA = zeros(Int64,n*nzd) + vA = zeros(Float64,n*nzd) + iP = zeros(Int64,nzd) + jP = zeros(Int64,nzd) + vP = zeros(Float64,nzd) + new(periods,iA,jA,vA,iP,jP,vP) + end +end + +function perfect_foresight_model!(Y::Array{Float64,1},exogenousvariables::Array{Float64,2},model::Model,steadystate::Array{Float64,1},ws::PerfectForesightModelWS,residuals::Array{Float64,1}) + periods = ws.periods + ny = length(model.endo) + i_rows = collect(1:ny) + i_cols = find(model.lead_lag_incidence') + m = 0 + for it = 2:(periods+1) + res = sub(residuals,m+1:m+model.eq_nbr) + Yview = sub(Y,i_cols) + model.dynamic(Yview, exogenousvariables, model.params, steadystate, it, res) + m += model.eq_nbr + i_cols += ny + end +end + +function perfect_foresight_model!(Y::Array{Float64,1},exogenousvariables::Array{Float64,2},model::Model,steadystate::Array{Float64,1},ws::PerfectForesightModelWS,A::SparseMatrixCSC{Float64,Int64}) + periods = ws.periods + ny = length(model.endo) + i_rows = collect(1:ny) + i_cols = find(model.lead_lag_incidence') + m = 0 + offset_r = 0 + offset_c = -ny + nzd = 0 + for it = 2:(periods+1) + Yview = sub(Y,i_cols) + if it == 2 + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, ws.iP, ws.jP, ws.vP) + nzd = count(i->(0 .< i .<= 3*ny),ws.jP) + k1 = 1 + for k = 1:nzd + if ws.jP[k] > ny + ws.iA[k1] = ws.iP[k] + ws.jA[k1] = ws.jP[k] + offset_c + ws.vA[k1] = ws.vP[k] + k1 += 1 + end + end + m = k1 - 1 + elseif it==(periods+1) + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, ws.iP, ws.jP, ws.vP) + for k=1:nzd + if ws.jP[k] <= 2*ny + m += 1 + ws.iA[m] = ws.iP[k] + offset_r + ws.jA[m] = ws.jP[k] + offset_c + ws.vA[m] = ws.vP[k] + end + end + else + I = sub(ws.iA,m+1:m+nzd) + J = sub(ws.jA,m+1:m+nzd) + V = sub(ws.vA,m+1:m+nzd) + model.first_derivatives(Y[i_cols], exogenousvariables, model.params, steadystate, it, I,J,V) + for k=m+1:m+nzd + ws.iA[k] += offset_r + ws.jA[k] += offset_c + end + m += nzd + end + offset_r += ny + offset_c += ny + i_cols += ny + end + A = sparse(ws.iA[1:m], ws.jA[1:m], ws.vA[1:m]) + assert(size(A)==(ny*periods,ny*periods)) + return A +end + +end diff --git a/julia/v0.5/v0.5/SteadyState.jl b/julia/v0.5/v0.5/SteadyState.jl new file mode 100644 index 0000000000000000000000000000000000000000..8d628563b6074ebbe5079189176e2bdd69ee567c --- /dev/null +++ b/julia/v0.5/v0.5/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/julia/v0.5/v0.5/Utils.jl b/julia/v0.5/v0.5/Utils.jl new file mode 100644 index 0000000000000000000000000000000000000000..6fd7eff5316a3f97410331793983eed1a5404f6b --- /dev/null +++ b/julia/v0.5/v0.5/Utils.jl @@ -0,0 +1,36 @@ +module Utils +## + # Copyright (C) 2015 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/>. +## + +export get_power_deriv + +function get_power_deriv(x::Float64, p::Real, k::Int) + if abs(x)<1e-12 && p>0 && k>p && typeof(p)==Int + dxp = .0 + else + dxp = x^(p-k) + for i = 0:k-1 + dxp *= p + p -= 1 + end + end + return dxp +end + +end