additional files

src/data.jl

+function identify_period_frequency(period)::Periods.Frequency
+    period = uppercase(period)
+    if 'Y' in period
+        frequency = Periods.Year
+    elseif 'A' in period
+        frequency = Periods.Year
+    elseif 'S' in period
+        frequency = Period.Semester
+    elseif 'H' in period
+        frequency = Periods.Semester
+    elseif 'Q' in period
+        frequency = Periods.Quarter
+    elseif 'M' in period
+        frequency = Periods.Month
+    elseif 'W' in period
+        frequency = Periods.Week
+    elseif 'D' in period
+        frequency = Periods.Day
+    elseif (cdash = count("-", period)) == 2
+        frequency = Periods.Day
+    elseif cdash == 1
+        frequency = Periods.Month
+    else
+        throw(ErrorException)
+    end
+end
+
+function get_data(filename, variables, options)
+    df = DataFrame(CSV.File(filename))
+    if uppercase(names(df)[1]) in ["COLUMN1", "DATE", "DATES",
+                                   "PERIOD", "PERIODS"]
+        frequency = identify_period_frequency(uppercase(df[1, 1]))
+    else
+        frequency = Periods.Undated
+    end
+    
+    ny = length(variables)
+
+    start = get(options, "first_obs", 1)
+    last = get(options, "last_obs", size(df, 1) - start + 1) 
+    nobs = last - start + 1
+    Y = Matrix{Union{Missing,Float64}}(undef, ny, nobs)
+    for (i, v) in enumerate(variables)
+        Y[i, :] .= df[start:last, v]
+    end
+    return Y
+end
+

src/graphics.jl

+function plot(variables::Vector{Vector{Float64}};
+              bar_variable = [],
+              first_period = 1,
+              plot_title = "Smoothed value",
+              plot_legend = (),
+              plot_filename = "",
+              plot_legend_position = :topright)
+    local myplot
+    # deal first with bar variable
+    if length(bar_variable) > 0
+        variables = pushfirst!(variables, bar_variable)
+    end
+    for (i, v) in enumerate(variables)
+        if length(plot_legend) > 0
+            thislabel = plot_legend[i]
+        else
+            thislabel = ""
+            plot_legend_position = false
+        end
+        if i == 1
+            if length(bar_variable) > 0
+                len = length(bar_variable)
+                xb = collect(range(first_period, length=len)) 
+                myplot = Plots.bar(xb, bar_variable, label = thislabel,
+                                   legend=plot_legend_position,
+                                   title = plot_title)
+                twinx()
+            else
+                x = collect(range(first_period, length=length(v)))
+                myplot = Plots.plot(x, v, label = thislabel,
+                                    legend=plot_legend_position,
+                                    title = plot_title,
+                                    linewidth = 3)
+            end
+        else
+            x = collect(range(first_period, length=length(v)))
+            myplot = Plots.plot!(x, v,
+                                 label = thislabel,
+                                 linewidth = 3)
+        end
+    end
+    Plots.display(myplot)
+    if length(plot_filename) > 0
+        Plots.savefig(plot_filename)
+    end
+end
+
+function plot(variable::Vector{Float64};
+              bar_variable = [],
+              first_period = 1,
+              plot_title = "Smoothed value",
+              plot_legend = (),
+              plot_filename = "",
+              plot_legend_position = :topright)
+    local myplot
+    # deal first with bar variable
+    if length(plot_legend) > 0
+        thislabel = plot_legend
+    else
+        thislabel = ""
+        plot_legend_position = false
+    end
+    if length(bar_variable) > 0
+        len = length(bar_variable)
+        xb = collect(range(first_period, length=len))
+        myplot = Plots.bar(xb, bar_variable, label = thislabel[1],
+                           legend=plot_legend_position,
+                           title = plot_title)
+        x = collect(range(first_period, length=length(variable)))
+        lims = Plots.ignorenan_extrema(myplot[1].attr[:yaxis])
+        m, M = extrema(variable)
+        tb = (lims[2] - lims[1])/(M-m)
+        variable = tb*(variable .- m) .+ lims[1]
+        Plots.plot!(x, variable, label = thislabel[2],
+                             title = plot_title,
+                    linewidth = 3)
+        myplot = twinx()
+        @show (m, M)
+        myplot = Plots.plot!(myplot, ylims = (m, M))
+    else
+        x = collect(range(first_period, length=length(variable)))
+        Plots.plot(x, variable, label = thislabel[1],
+                   legend=plot_legend_position,
+                   title = plot_title,
+                   linewidth = 3)
+    end
+    Plots.display(myplot)
+    if length(plot_filename) > 0
+        Plots.savefig(plot_filename)
+    end
+end

src/initialization.jl

+function histval!(context, field)
+    symboltable = context.symboltable
+    m = context.models[1]
+    histval = zeros(m.orig_maximum_lag, m.endogenous_nbr)
+    for v in field["vals"]
+        k = symboltable[v["name"]].orderintype
+        l = m.orig_maximum_lag - v["lag"]
+        histval[l, k] = Meta.parse(v["value"])
+    end
+    context.work.histval = histval
+end
+
+function param_init!(context, field)
+    params = context.work.params;
+    symboltable = context.symboltable
+    s = symboltable[field["name"]]
+    k = s.orderintype
+    params[k] = eval(Meta.parse(field["value"]))
+end
+
+function initval!(context, field)
+end
+
+function shocks!(context, field)
+    Sigma = context.models[1].Sigma_e
+    symboltable = context.symboltable
+    set_variance!(Sigma, field["variance"], symboltable)
+    set_stderr!(Sigma, field["stderr"], symboltable)
+    set_covariance!(Sigma, field["covariance"], symboltable)
+    set_correlation!(Sigma, field["correlation"], symboltable)
+end
+
+function set_variance!(Sigma, variance, symboltable)
+    for v in variance
+        k =  symboltable[v["name"]].orderintype
+        Sigma[k, k] = eval(Meta.parse(v["variance"]))
+    end
+end
+
+function set_stderr!(Sigma, stderr, symboltable)
+    for s in stderr
+        k =  symboltable[s["name"]].orderintype
+        x = eval(Meta.parse(s["stderr"]))
+        Sigma[k, k] = x*x
+    end
+end
+
+function set_covariance!(Sigma, covariance, symboltable)
+    for c in covariance
+        k1 =  symboltable[c["name"]].orderintype
+        k2 =  symboltable[c["name2"]].orderintype
+        Sigma[k1, k2] = eval(Meta.parse(c["covariance"]))
+        Sigma[k2, k1] = Sigma[k1, k2]
+    end
+end
+
+function set_correlation!(Sigma, correlation, symboltable)
+    for c in correlation
+        k1 =  symboltable[c["name"]].orderintype
+        k2 =  symboltable[c["name2"]].orderintype
+        corr = eval(Meta.parse(c["correlation"]))
+        Sigma[k2, k1] = sqrt(Sigma[k1, k1]*Sigma[k2, k2])*corr
+    end
+end

src/optimal_policy.jl

+function planner_objective!(context, field)
+end

src/perturbations.jl

+function display_stoch_simul(x, title, context)
+    endogenous_names = get_endogenous_longname(context.symboltable)
+    emptyrow = ["" for _= 1:size(x,1)]
+    column_header = []
+    #    map(x -> push!(column_header, string(x, "\U0209C")), endogenous_names)
+    map(x -> push!(column_header, "$(x)_t"), endogenous_names)
+    row_header = [""]
+    map(x -> push!(row_header, "ϕ($x)"), endogenous_names[context.models[1].i_bkwrd_b])
+    map(x -> push!(row_header, "$(x)_t"), get_exogenous_longname(context.symboltable))
+    data = hcat(row_header,
+                vcat(reshape(column_header, 1, length(column_header)),
+                     x))
+    # Note: ϕ(x) = x_{t-1} - \bar x
+    #    note = string("Note: ϕ(x) = x\U0209C\U0208B\U02081 - ", "\U00305", "x")
+    note = string("Note: ϕ(x) = x_{t-1} - steady_state(x)")
+    println("\n")
+    dynare_table(data, title, column_header, row_header, note)
+end
+
+function make_A_B!(A, B, model, results)
+    vA = view(A, :, model.i_bkwrd_b)
+    vA .= results.linearrationalexpectations.g1_1
+    B .= results.linearrationalexpectations.g1_2
+end
+
+function stoch_simul!(context, field)
+    model = context.models[1]
+    options = context.options
+    results = context.results.model_results[1]
+    work = context.work
+    options["stoch_simul"] = Dict()
+    copy!(options["stoch_simul"], field["options"])
+    compute_stoch_simul!(context)
+    compute_variance!(context)
+    x = results.linearrationalexpectations.g1
+    vx = view(x, :, 1:size(x, 2) - 1)
+    display_stoch_simul(vx', "Coefficients of approximate solution function", context)
+    if (periods = get(options["stoch_simul"], "periods", 0)) > 0
+        simulresults = Matrix{Float64}(undef, periods + 1, model.endogenous_nbr)
+        histval = work.histval
+        if size(histval, 1) == 0
+            y0 = results.endogenous_steady_state
+        else
+            y0 = view(histval, 1, :)
+        end
+        C = cholesky(model.Sigma_e)
+        x = vcat(zeros(1, model.exogenous_nbr),
+                 randn(periods, model.exogenous_nbr)*C.U)
+        c = results.endogenous_steady_state
+        A = zeros(model.endogenous_nbr, model.endogenous_nbr)
+        B = zeros(model.endogenous_nbr, model.exogenous_nbr)
+        make_A_B!(A, B, model, results)
+        simul_first_order!(simulresults, y0, x, c, A, B, periods)
+        if work.model_has_trend
+            simulresults .+= transpose(results.endogenous_trend[:,1]) .+ collect(0:size(simulresults, 1)-1) * transpose(results.endogenous_trend[:, 2]) 
+        end
+        first_period = get(options["stoch_simul"], "first_periods", 1)
+        endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)]
+        data = TimeDataFrame(simulresults, Periods.Undated, first_period, endogenous_names)
+        push!(results.simulations, Simulation("", "stoch_simul", options["stoch_simul"], data))
+    end
+end
+
+function check(context, field)
+end
+
+function compute_stoch_simul!(context)
+    m = context.models[1]
+    results = context.results.model_results[1]
+    options = context.options["stoch_simul"]
+    options["cyclic_reduction"] = Dict()
+    options["generalized_schur"] = Dict()
+    work = context.work
+    Base.invokelatest(steady_state!, context)
+    fill!(results.trends.exogenous_steady_state, 0.0)
+    get_jacobian!(work,
+                  results.trends.endogenous_steady_state,
+                  results.trends.exogenous_steady_state, 
+                  results.trends.endogenous_steady_state,
+                  m,
+                  2)
+    if isnothing(get(options,"dr_cycle_reduction", nothing))
+        algo = "GS"
+    else
+        algo = "CR"
+    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!(work.jacobian, ws)
+    if algo == "GS"
+        LinearRationalExpectations.get_de!(ws, work.jacobian)
+        options["generalized_schur"]["criterium"] = 1 + 1e-6
+    else
+        LinearRationalExpectations.get_abc!(ws, work.jacobian)
+        options["cyclic_reduction"]["tol"] = 1e-8
+    end
+    LinearRationalExpectations.first_order_solver!(results.linearrationalexpectations,
+                                                   algo,
+                                                   work.jacobian,
+                                                   options,
+                                                   ws)
+end
+
+function compute_variance!(context)
+    m = context.models[1]
+    results = context.results.model_results[1]
+    Σy = results.endogenous_variance
+    A = results.linearrationalexpectations.gs1
+    B1 = zeros(m.n_states, m.exogenous_nbr)
+    g1_1 = results.linearrationalexpectations.g1_1
+    g1_2 = results.linearrationalexpectations.g1_2
+    vr1 = view(g1_2, m.i_bkwrd_b, :)
+    B1 .= vr1
+    ws = LyapdWs(m.n_states)
+    Σ = zeros(m.n_states, m.n_states)
+    tmp = zeros(m.n_states, m.exogenous_nbr)
+    mul!(tmp, B1, m.Sigma_e)
+    B = zeros(m.n_states, m.n_states)
+    mul!(B, tmp, B1')
+    extended_lyapd!(Σ, A, B, ws)
+    stationary_variables = results.stationary_variables
+    fill!(stationary_variables, true)
+    state_stationary_variables = view(stationary_variables, m.i_bkwrd_b)
+    nonstate_stationary_variables = view(stationary_variables, m.i_non_states)
+    if any(ws.nonstationary_variables)
+        fill!(Σy, NaN)
+        state_stationary_variables .= .!ws.nonstationary_variables
+        state_stationary_nbr = count(state_stationary_variables)
+        vr3 = view(results.linearrationalexpectations.g1_1, m.i_non_states, :)
+        for i = 1:(m.endogenous_nbr - m.n_states)
+            for j = 1:m.n_states
+                if ws.nonstationary_variables[j] && abs(vr3[i, j]) > 1e-10
+                    nonstate_stationary_variables[j] = false
+                    break
+                end
+            end
+        end
+        nonstate_stationary_nbr = count(nonstate_stationary_variables)
+    else
+        state_stationary_nbr = m.n_states 
+        nonstate_stationary_nbr = m.endogenous_nbr - m.n_states
+    end
+    # state / state
+    stationary_nbr = state_stationary_nbr + nonstate_stationary_nbr
+    A2 = zeros(nonstate_stationary_nbr, state_stationary_nbr)
+    vtmp = view(results.linearrationalexpectations.g1_1, m.i_non_states, :)
+    A2 .= view(vtmp, nonstate_stationary_variables, state_stationary_variables)
+    vtmp = view(Σy, m.i_bkwrd_b, m.i_bkwrd_b)
+    vΣy = view(vtmp, state_stationary_variables, state_stationary_variables)
+    vΣ = view(Σ, state_stationary_variables, state_stationary_variables) 
+    vΣy .= vΣ
+    # endogenous / nonstate
+    n_non_states = m.endogenous_nbr - m.n_states
+    B2 = zeros(nonstate_stationary_nbr, m.exogenous_nbr)
+    vtmp = view(results.linearrationalexpectations.g1_2, m.i_non_states, :)
+    vr2 = view(vtmp, nonstate_stationary_variables, :)
+    B2 .= vr2
+    Σ = zeros(stationary_nbr, nonstate_stationary_nbr)
+    vg1 = view(results.linearrationalexpectations.g1_1, stationary_variables, state_stationary_variables)
+    tmp1 = zeros(stationary_nbr, state_stationary_nbr)
+    mul!(tmp1, vg1, vΣy)
+    vg1 = view(results.linearrationalexpectations.g1_1, m.i_non_states, :)
+    vg11 = view(vg1, nonstate_stationary_variables, state_stationary_variables)
+    mul!(Σ, tmp1, transpose(vg11))
+    tmp2 = zeros(stationary_nbr, m.exogenous_nbr)
+    vg2 = view(results.linearrationalexpectations.g1_2, stationary_variables, :)
+    mul!(tmp2, vg2, m.Sigma_e)
+    mul!(Σ, tmp2, transpose(B2), 1.0, 1.0)
+    vtmp = view(Σy, :, m.i_non_states) 
+    vΣy = view(vtmp, stationary_variables, nonstate_stationary_variables)
+    vΣy .= Σ
+    # nonstate / state
+    vtmp1 = view(Σy, m.i_non_states, m.i_bkwrd_b)
+    vΣy1 = view(vtmp1, nonstate_stationary_variables, state_stationary_variables)
+    vtmp2 = view(Σy, m.i_bkwrd_b, m.i_non_states)
+    vΣy2 = view(vtmp2, state_stationary_variables, nonstate_stationary_variables)
+    vΣy1 .= transpose(vΣy2)
+end
+    
+function simul_first_order!(results, initial_values, x, c, A, B, periods)
+    r_1 = view(results, 1, :)
+    r_1 .= initial_values .- c
+    for t = 2:periods + 1
+        r = view(results, t, :)
+        e = view(x, t, :)
+        mul!(r, B, e)
+        mul!(r, A, r_1, 1.0, 1.0)
+        r_1 .+=  c
+        r_1 = r
+    end
+    r_1 .+= c
+end

test/models/cgg/cgg_ramsey.mod

+var y pie r;
+varexo e_y e_pie;
+
+parameters delta sigma alpha kappa;
+
+delta =  0.44;
+kappa =  0.18;
+alpha =  0.48;
+sigma = -0.06;
+
+model(linear);
+y  = delta * y(-1)  + (1-delta)*y(+1)+sigma *(r - pie(+1)) + e_y; 
+pie  =   alpha * pie(-1) + (1-alpha) * pie(+1) + kappa*y + e_pie;
+end;
+
+shocks;
+var e_y;
+stderr 0.63;
+var e_pie;
+stderr 0.4;
+end;
+
+planner_objective pie^2 + y^2;
+
+ramsey_model(planner_discount=1.0);
+
+stoch_simul(irf=0);
+