adding macro @Dynare and new functions

Manifest.toml

+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.0"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[Parsers]]
+deps = ["Dates", "Test"]
+git-tree-sha1 = "75d07cb840c300084634b4991761886d0d762724"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "1.0.1"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

Project.toml

@@ -2,3 +2,7 @@ name = "Dynare"
 uuid = "5203de40-99df-439e-afbc-014de65cb9ef"
 authors = ["michel "]
 version = "0.1.0"
+
+[deps]
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/Dynare.jl

 module Dynare
 
 include("model.jl")
+export get_abc, get_de
 include("parser/DynareParser.jl")
+export parser, get_jacobian_at_steadystate!
 include("parser/DynarePreprocessor.jl")
+export dynare_preprocess
+include("steady_state/SteadyState.jl")
+export steady_state!
+include("Utils.jl")
+export get_power_deriv
+
+export @dynare
+
+macro dynare(modfilename, args...)
+    if  occursin(r"\.mod$", modfilename)
+        modname = modfilename[1:length(modfilename)-4]
+    else
+        modname = modfilename
+    end
+    dynare_preprocess(modname*".mod", args)
+    context = Dynare.Context()
+    parser(modname, context)
+    return context
+end
 
-export parser, dynare_preprocess
 end # module

src/model.jl

@@ -2,9 +2,11 @@ using LinearAlgebra
 export Model, get_de, get_abc, inverse_order_of_dynare_decision_rule
 
 struct Model
-    endo_nbr
-    current_exogenous_nbr
+    endogenous_nbr
+    exogenous_nbr
     lagged_exogenous_nbr
+    exogenous_deterministic_nbr
+    parameter_nbr
     lead_lag_incidence
     n_static
     n_fwrd
@@ -61,7 +63,7 @@ struct Model
     steady_state!
 end
 
-function Model(modfilename, endo_nbr, lead_lag_incidence, current_exogenous_nbr, lagged_exogenous_nbr)
+function Model(modfilename, endo_nbr, lead_lag_incidence, exogenous_nbr, lagged_exogenous_nbr, exogenous_deterministic_nbr, parameter_nbr)
     i_static = findall((lead_lag_incidence[1,:] .== 0) .& (lead_lag_incidence[3,:] .== 0))
     p_static = lead_lag_incidence[2,i_static]
     i_dyn = findall((lead_lag_incidence[1,:] .> 0) .| (lead_lag_incidence[3,:] .> 0))
@@ -98,7 +100,7 @@ function Model(modfilename, endo_nbr, lead_lag_incidence, current_exogenous_nbr,
     i_cur_both = findall(lead_lag_incidence[2,i_both] .> 0)
     n_cur_both = length(i_cur_both)
     p_cur_both = lead_lag_incidence[2,i_both[i_cur_both]]
-    icolsD = [1:n_cur_bkwrd; n_bkwrd+n_both .+ (1:(n_fwrd+n_both))]
+    icolsD = [1:n_cur_bkwrd; n_bkwrd + n_both .+ (1:(n_fwrd+n_both))]
     jcolsD = [p_cur_bkwrd; p_fwrd; p_both_f]
     # derivatives of current values of variables that are both
     # forward and backward are included in the E matrix
@@ -111,14 +113,18 @@ function Model(modfilename, endo_nbr, lead_lag_incidence, current_exogenous_nbr,
     DErows2 = (n_dyn-n_both) .+ (1:n_both)
     gx_rows = n_bkwrd .+ (1:(n_fwrd+n_both))
     hx_rows = 1:(n_bkwrd + n_both)
-    i_current_exogenous = maximum(lead_lag_incidence) .+ (1:current_exogenous_nbr)
+    i_current_exogenous = maximum(lead_lag_incidence) .+ (1:exogenous_nbr)
     i_lagged_exogenous = 0:-1
     serially_correlated_exogenous = false
     dynamic! = load_dynare_function(modfilename*"Dynamic.jl")
     static! = load_dynare_function(modfilename*"Static.jl")
-    steady_state! = load_dynare_function(modfilename*"SteadyState2.jl")
-    Model(endo_nbr, current_exogenous_nbr, lagged_exogenous_nbr,
-          lead_lag_incidence, n_static, n_fwrd, n_bkwrd, n_both,
+    if isfile(modfilename*"SteadyState2.jl")
+        steady_state! = load_dynare_function(modfilename*"SteadyState2.jl")
+    else
+        steady_state! = Nothing
+    end
+    Model(endo_nbr, exogenous_nbr, lagged_exogenous_nbr, exogenous_deterministic_nbr,
+          parameter_nbr, lead_lag_incidence, n_static, n_fwrd, n_bkwrd, n_both,
           n_states, DErows1, DErows2, n_dyn, i_static, i_dyn, i_bkwrd,
           i_bkwrd_b, i_bkwrd_ns, i_fwrd, i_fwrd_b, i_fwrd_ns, i_both,
           p_static, p_bkwrd, p_bkwrd_b, p_fwrd, p_fwrd_b, p_both_b,
@@ -132,34 +138,6 @@ function Model(modfilename, endo_nbr, lead_lag_incidence, current_exogenous_nbr,
           steady_state!)   
 end
 
-Model(modfilename, endo_nbr, lli, current_exogenous_nbr) =
-    Model(modfilename, endo_nbr, lli, current_exogenous_nbr, 0)
-    
-function get_de(jacobian,model)
-    n1 = size(model.DErows1,1)
-    n2 = model.n_dyn - n1;
-    d = zeros(model.n_dyn,model.n_dyn)
-    e = zeros(model.n_dyn,model.n_dyn)
-    d[1:n1,model.icolsD] = jacobian[:,model.jcolsD]
-    e[1:n1,model.icolsE] = -jacobian[:,model.jcolsE]
-    u = Matrix{Float64}(I, n2, n2)                                    
-    d[model.DErows2,model.colsUD] = u
-    e[model.DErows2,model.colsUE] = u
-    return d, e
-end
-
-function get_abc(model::Model,jacobian::Array{Float64})
-    i_rows = model.n_static+1:model.endo_nbr
-    n = length(i_rows)
-    a = zeros(Float64,n,n)
-    b = zeros(Float64,n,n)
-    c = zeros(Float64,n,n)
-    a[:,model.i_bkwrd_ns] = view(jacobian,i_rows,model.p_bkwrd_b)
-    b[:,model.i_current_ns]  = view(jacobian,i_rows,model.p_current_ns)
-    c[:,model.i_fwrd_ns]  = view(jacobian,i_rows,model.p_fwrd_b)
-    return a, b, c
-end
-
 function inverse_order_of_dynare_decision_rule(m::Model)
     inverse_order_var = Vector{Int64}(undef, m.endo_nbr)
     for i = 1:m.n_static
@@ -189,8 +167,7 @@ end
 function load_dynare_function(filename)
     file = readlines(filename)
     # drop using Utils
-    deleteat!(file, 6)
-    str = join(file, "\n")
-    return eval(Meta.parse(str))
+    file[6] = "using Dynare: get_power_deriv"
+    return eval(Meta.parse(join(file, "\n")))
 end
 

src/parser/DynareParser.jl

+using FastLapackInterface
+using FastLapackInterface.LinSolveAlgo
 using JSON
+using LinearRationalExpectations
 
 @enum SymbolType Endogenous Exogenous ExogenousDeterministic Parameter DynareFunction
 
@@ -12,20 +15,61 @@ end
 struct Options
 end
 
+struct ModelResults
+    endogenous_steady_state::Vector{Float64}
+    exogenous_steady_state::Vector{Float64}
+    exogenous_deterministic_steady_state::Vector{Float64}
+end
+
 struct Results
+    model_results::Vector{ModelResults}
 end
 
-struct Context
+mutable struct Work
+    params::Vector{Float64}
+    residuals::Vector{Float64}
+    temporary_values::Vector{Float64}
+    dynamic_variables::Vector{Float64}
+    jacobian::Matrix{Float64}
+    qr_jacobian::Matrix{Float64}
+end
+
+function work_allocate(w::Work, m::Model)
+    if length(w.params) == 0
+        resize!(w.params, m.parameter_nbr)
+    end
+    if length(w.residuals) == 0
+        resize!(w.residuals, m.endogenous_nbr)
+    end
+    if length(w.temporary_values) == 0
+        resize!(w.temporary_values, sum(m.dynamic!.tmp_nbr))
+    end
+    ncol = m.n_bkwrd + m.n_current + m.n_fwrd + 2*m.n_both
+    if length(w.dynamic_variables) == 0
+        resize!(w.dynamic_variables, ncol)
+    end
+    ncol1 = ncol + m.exogenous_nbr
+    if length(w.jacobian) == 0
+        w.jacobian = Matrix{Float64}(undef, m.endogenous_nbr, ncol1)
+    end
+    if length(w.qr_jacobian) == 0
+        w.qr_jacobian = Matrix{Float64}(undef, m.endogenous_nbr, ncol1)
+    end
+end          
+
+mutable struct Context
     symboltable::Dict{String, Symbol}
     models::Vector{Model}
     options::Dict{String, Any}
     results::Results
+    work::Work
     function Context()
         symboltable = Dict()
         models = Vector{Model}(undef,1)
         options = Dict()
-        results = Results()
-        new(symboltable, models, options, results)
+        results = Results([ModelResults([],[],[])])
+        work = Work([], [], [], [], Matrix{Float64}(undef, 0, 0), Matrix{Float64}(undef, 0, 0))
+        new(symboltable, models, options, results, work)
     end
 end
 
@@ -48,16 +92,22 @@ function parser(modfilename, context::Context)
     exo_nbr = set_symbol_table!(context.symboltable, modeljson["exogenous"], Exogenous)
     exo_det_nbr = set_symbol_table!(context.symboltable, modeljson["exogenous_deterministic"], ExogenousDeterministic)
     param_nbr = set_symbol_table!(context.symboltable, modeljson["parameters"], Parameter)
-    params = Vector{Float64}(undef, param_nbr)
     Sigma_e = zeros(exo_nbr, exo_nbr)
     model_info = get_model_info(modeljson["model_info"])
     context.models[1] = Model(modfilename,
                               endo_nbr,
                               model_info.lead_lag_incidence,
-                              exo_nbr)
+                              exo_nbr,
+                              0,
+                              exo_det_nbr,
+                              param_nbr)
+    context.results = Results([ModelResults(Vector{Float64}(undef, endo_nbr),
+                                    Vector{Float64}(undef, exo_nbr),
+                                            Vector{Float64}(undef, exo_det_nbr))])
+    work_allocate(context.work, context.models[1])
     for field in modeljson["statements"]
         if field["statementName"] == "param_init"
-            initialize_parameter!(params, field, context.symboltable)
+            initialize_parameter!(context.work.params, field, context.symboltable)
         elseif field["statementName"] == "native"
             native_statement(field)
         elseif field["statementName"] == "initval"
@@ -184,9 +234,57 @@ end
 function check(field)
 end
 
-function compute_stoch_simul(context); end;
-function compute_prefect_foresight_setup(context); end;
-function compute_perfect_foresight_solver(context); end;
+function compute_stoch_simul(context)
+    m = context.models[1]
+    results = context.results
+    if context.options["stoch_simul"]["dr_cycle_reduction"]
+        algo = "CR"
+    else
+        algo = "GS"
+    end
+    ws = LinearRationalExpectationsWs(algo,
+                                      m.endogenous_nbr,
+                                      m.exogenous_nbr,
+                                      m.exogenous_deterministic_nbr,
+                                      m.i_fwrd_b,
+                                      m.i_current,
+                                      m.i_bkwrd_b,
+                                      m.i_both,
+                                      m.i_static)
+    LinearRationalExpectations.remove_static!(jacobian, ws)
+    if algo == "GS"
+        LinearRationalExpectations.get_de!(ws, jacobian)
+    else
+        LinearRationalExpectations.get_abc!(ws, jacobian)
+    end
+    LinearRationalExpectations.first_order_solver!(results, algo, jacobian, options, ws)
+    println(results)
+end
 
+function compute_prefect_foresight_setup(context); end
+function compute_perfect_foresight_solver(context); end
 
+function get_dynamic_variables!(y::Vector{Float64}, steadystate::Vector{Float64}, lli::Matrix{Int64})
+    for i = 1:size(lli,2)
+        value = steadystate[i]
+        for j = 1:size(lli,1)
+            k = lli[j, i]
+            if k > 0
+                y[k] = value
+            end
+        end
+    end
+end
 
+function get_jacobian_at_steadystate!(work::Work, steadystate, exogenous, m::Model, period::Int64)
+    lli = m.lead_lag_incidence
+    get_dynamic_variables!(work.dynamic_variables, steadystate, lli)
+    m.dynamic!.dynamic!(work.temporary_values,
+                        work.residuals,
+                        work.jacobian,
+                        work.dynamic_variables,
+                        exogenous,
+                        work.params,
+                        steadystate,
+                        period)  
+end

src/parser/DynarePreprocessor.jl

-DYNARE_ROOT = "/data/projects/dynare/git/preprocessor/src/dynare_m"
-
-dynare_preprocess(modfilename) = run(`$DYNARE_ROOT $modfilename language=julia output=third json=compute`) 
+DYNARE_BINARY = "/data/projects/dynare/git/preprocessor/src/dynare_m"
 
+function dynare_preprocess(modfilename, args)
+    dynare_args = [basename(modfilename), "language=julia", "output=third", "json=compute"]
+    offset = 0
+    for a in args
+        if occursin(r"^output=", a)
+            deleteat!(dynare_args, 3)
+            offset = 1
+        elseif occursin(r"^json=", a)
+            deleteat!(dynare_args, 4 - offset)
+        end
+    end
+    append!(dynare_args, args)
+    println(dynare_args)
+    current_directory = pwd()
+    directory = dirname(modfilename)
+    cd(directory)
+    run(`$DYNARE_BINARY $dynare_args`)
+    cd(current_directory)
+end

test/test_DynareParser.jl

@@ -5,3 +5,5 @@ context = Dynare.Context()
 
 parser("models/example1/example1", context)
 println(context)
+steady_state!(context)
+println(context.results)