From 7f813c0427109a81edb6f1e69c098f28347b33d5 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Wed, 2 Aug 2023 10:42:42 -0400 Subject: [PATCH 01/10] fix wasted memory creating dense zeros before conversion to sparse --- src/model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/model.jl b/src/model.jl index 61721b9f..bbb7d26c 100644 --- a/src/model.jl +++ b/src/model.jl @@ -29,7 +29,7 @@ struct DynamicWs ] for i = 2:order push!(derivatives, - zeros(endogenous_nbr, + spzeros(endogenous_nbr, (3*endogenous_nbr + exogenous_nbr)^i) ) end From b5f1278708e6947932709be4d1ab6f256fc1f4a9 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Wed, 2 Aug 2023 13:46:20 -0400 Subject: [PATCH 02/10] keep F derivatives as sparse in Dynare --- src/perturbations.jl | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/src/perturbations.jl b/src/perturbations.jl index da5f2043..c4b716be 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -518,26 +518,11 @@ function compute_stoch_simul!( model, ) m = model - k = vcat( - m.i_bkwrd_b, - m.endogenous_nbr .+ model.i_current, - 2*m.endogenous_nbr .+ model.i_fwrd_b, - 3*m.endogenous_nbr .+ collect(1:m.exogenous_nbr) - ) - n = 3*m.endogenous_nbr + m.exogenous_nbr - nc = m.n_bkwrd + 2*m.n_both + m.n_current + m.n_fwrd + m.exogenous_nbr - F1 = zeros(m.endogenous_nbr, nc) - @views F1 .= ws.derivatives[1][:, k] - kk = vec(reshape(1:n*n, n, n)[k, k]) - - sp = ws.derivatives[2] - F2 = Matrix(ws.derivatives[2])[:, kk] - F = [F1, F2] G = Vector{Matrix{Float64}}(undef, 0) push!(G, context.results.model_results[1].linearrationalexpectations.g1) push!(G, zeros(model.endogenous_nbr, (model.n_states + model.exogenous_nbr + 1)^2)) - ws = KOrderPerturbations.KOrderWs(model.endogenous_nbr, + solverWs = KOrderPerturbations.KOrderWs(model.endogenous_nbr, model.n_fwrd + model.n_both, model.n_states, model.n_current, @@ -549,10 +534,10 @@ function compute_stoch_simul!( order) moments = [0, vec(model.Sigma_e)] KOrderPerturbations.k_order_solution!(G, - F, + ws.derivatives, moments, order, - ws) + solverWs) copy!(results.solution_derivatives, G) end From 25d3357e6a63a52e0184a7c604c3ac870077049a Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Thu, 3 Aug 2023 13:15:58 -0400 Subject: [PATCH 03/10] fix filename typo in test_scenario.jl --- test/test_scenario.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_scenario.jl b/test/test_scenario.jl index 733cf7d7..c3151edd 100644 --- a/test/test_scenario.jl +++ b/test/test_scenario.jl @@ -5,7 +5,7 @@ using JLD2 using LinearAlgebra using Test -@dynare "models/example1ss/example1ss" +@dynare "models/example1pf/example1pf_conditional" context = load("models/example1pf/example1pf_conditional/output/example1pf_conditional.jld2", "context") function make_f_J(context, scenario, periods, infoperiod) @@ -305,7 +305,7 @@ end @testset "example1pf_conditional" begin context = @dynare "models/example1pf/example1pf_conditional" - e = AxisArrayTables.data(simulation(:e)) + e = AxisArrayTables.data(simulation(:e)) u = AxisArrayTables.data(simulation(:u)) y = AxisArrayTables.data(simulation(:y)) c1 = AxisArrayTables.data(simulation(:c, simnbr=1)) From 10d2b7cb474270ff2fcf614eeb58ac692fa5d697 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Fri, 11 Aug 2023 18:37:08 -0400 Subject: [PATCH 04/10] WIP: IRFS for 2nd order simulation --- src/perturbations.jl | 74 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/src/perturbations.jl b/src/perturbations.jl index c4b716be..8b67d167 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -442,6 +442,19 @@ function stoch_simul_core!(context::Context, ws::DynamicWs, options::StochSimulO filename, ) end + # TODO: we should probably have a block like this + # else if options.order == 2 && options.irf > 0 + # irfs2(context, options.irf) + # path = "$(context.modfileinfo.modfilepath)/graphs/" + # mkpath(path) + # filename = "$(path)/irfs2" + # plot_irfs( + # context.results.model_results[1].irfs, + # model, + # context.symboltable, + # filename, + # ) + # end if (periods = options.periods) > 0 steadystate = results.trends.endogenous_steady_state linear_trend = results.trends.endogenous_linear_trend @@ -538,11 +551,72 @@ function compute_stoch_simul!( moments, order, solverWs) + # FOR DEBUGGING: single simulation + KOrderPerturbations.simulate_run(G, 100, solverWs) + + irfs2(G, 100, context, solverWs) + path = "$(context.modfileinfo.modfilepath)/graphs/" + filename = "$(path)/irfs2" + + # for development: delete all files that starts with "irfs2" + filter(x -> startswith(x, "irfs2"), readdir("$(path)", join=true)) .|> x -> rm(x) + + display(context.results.model_results[1].irfs) + plot_irfs(context.results.model_results[1].irfs, model, context.symboltable, filename) + copy!(results.solution_derivatives, G) end end +function irfs2(GD, periods, context, solverWs::KOrderWs) + model = context.models[1] + results = context.results.model_results[1] + endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)] + exogenous_names = [Symbol(n) for n in get_exogenous_longname(context.symboltable)] + + nshock = solverWs.nshock + u_shock = [zeros(nshock) for _ in 1:periods] + u_shock_det = [zeros(nshock) for _ in 1:periods] + + gy1 = GD[1][:, 1] + y0 = zeros(size(gy1)[1]) + + n = length(y0) + simWs = SimulateWs(GD, n, solverWs) + + for exo_var in 1:model.exogenous_nbr + sigma = diag(model.Sigma_e)[exo_var] + det_shock = sqrt(sigma) + + delta_sims::Vector{Array{Array{Float64}}} = [] + for n_sim in 1:100 + for i in 2:periods + random = randn() * det_shock^2 + u_shock[i][exo_var] = random + end + + u_shock_det = deepcopy(u_shock) + u_shock_det[2][exo_var] += det_shock + simWs = SimulateWs(GD, n, solverWs) + rand_sim = simulate(GD, y0, u_shock, periods, simWs) + rand_det_sim = simulate(GD, y0, u_shock_det, periods, simWs) + + delta_sim = rand_det_sim - rand_sim + + push!(delta_sims, delta_sim) + end + tdf_vec = mean(delta_sims) + tdf_mat = reduce(vcat,transpose.(tdf_vec)) + + tdf = AxisArrayTable(transpose(tdf_mat'), Undated(1):Undated(periods), endogenous_names) + + results.irfs[exogenous_names[exo_var]] = tdf + end + +end + + function compute_first_order_solution!( context::Context, endogenous::AbstractVector{Float64}, From f33f10965ca2a2810a5620e3e2ec43e342becb24 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Fri, 11 Aug 2023 18:47:01 -0400 Subject: [PATCH 05/10] WIP modify tests to execute orders 1 and 2 --- test/models/example1/example1.mod | 3 ++- test/models/example2/example2.mod | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/test/models/example1/example1.mod b/test/models/example1/example1.mod index 14f02928..0af7773b 100644 --- a/test/models/example1/example1.mod +++ b/test/models/example1/example1.mod @@ -52,6 +52,7 @@ end; check; -stoch_simul(dr=cycle_reduction, order=1); +stoch_simul(dr=cycle_reduction, order=1, irf=100); +stoch_simul(dr=cycle_reduction, order=2, irf=100); diff --git a/test/models/example2/example2.mod b/test/models/example2/example2.mod index 3f80fac7..9f489821 100644 --- a/test/models/example2/example2.mod +++ b/test/models/example2/example2.mod @@ -54,5 +54,5 @@ end; check; -stoch_simul(order=1); - +stoch_simul(order=1, irf=100); +stoch_simul(order=2, irf=100); From e1e9b8f2bb1ba1389997b58cef4358847402349b Mon Sep 17 00:00:00 2001 From: MichelJuillard Date: Mon, 14 Aug 2023 10:20:58 +0200 Subject: [PATCH 06/10] comment out 2nd order IRFs set deterministic ut0 shock --- src/perturbations.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/perturbations.jl b/src/perturbations.jl index 8b67d167..af729e2c 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -408,6 +408,7 @@ function localapproximation!(; context::Context=context, ncol = m.n_bkwrd + m.n_current + m.n_fwrd + 2 * m.n_both tmp_nbr = m.dynamic_tmp_nbr ws = DynamicWs(context, order=options.order) + @show options.order stoch_simul_core!(context, ws, options) end @@ -552,18 +553,19 @@ function compute_stoch_simul!( order, solverWs) # FOR DEBUGGING: single simulation - KOrderPerturbations.simulate_run(G, 100, solverWs) - + ut0 = sqrt(model.Sigma_e[1,1]) + KOrderPerturbations.simulate_run(G, ut0, 100, solverWs) + #= irfs2(G, 100, context, solverWs) - path = "$(context.modfileinfo.modfilepath)/graphs/" - filename = "$(path)/irfs2" + #path = "$(context.modfileinfo.modfilepath)/graphs/" + #filename = "$(path)/irfs2" # for development: delete all files that starts with "irfs2" - filter(x -> startswith(x, "irfs2"), readdir("$(path)", join=true)) .|> x -> rm(x) + #filter(x -> startswith(x, "irfs2"), readdir("$(path)", join=true)) .|> x -> rm(x) display(context.results.model_results[1].irfs) plot_irfs(context.results.model_results[1].irfs, model, context.symboltable, filename) - + =# copy!(results.solution_derivatives, G) end From 90fe2719ee0bf09cc4fbaa95f71e57f0916ca068 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Wed, 16 Aug 2023 11:36:39 -0400 Subject: [PATCH 07/10] WIP IRFS --- src/perturbations.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/perturbations.jl b/src/perturbations.jl index af729e2c..953075a2 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -446,6 +446,7 @@ function stoch_simul_core!(context::Context, ws::DynamicWs, options::StochSimulO # TODO: we should probably have a block like this # else if options.order == 2 && options.irf > 0 # irfs2(context, options.irf) + irfs2(context, options.irf, solverWs) # path = "$(context.modfileinfo.modfilepath)/graphs/" # mkpath(path) # filename = "$(path)/irfs2" @@ -555,25 +556,25 @@ function compute_stoch_simul!( # FOR DEBUGGING: single simulation ut0 = sqrt(model.Sigma_e[1,1]) KOrderPerturbations.simulate_run(G, ut0, 100, solverWs) - #= - irfs2(G, 100, context, solverWs) - #path = "$(context.modfileinfo.modfilepath)/graphs/" - #filename = "$(path)/irfs2" + copy!(results.solution_derivatives, G) + irfs2(context, 100, solverWs) + path = "$(context.modfileinfo.modfilepath)/graphs/" + filename = "$(path)/irfs2" # for development: delete all files that starts with "irfs2" #filter(x -> startswith(x, "irfs2"), readdir("$(path)", join=true)) .|> x -> rm(x) display(context.results.model_results[1].irfs) - plot_irfs(context.results.model_results[1].irfs, model, context.symboltable, filename) - =# - copy!(results.solution_derivatives, G) + # plot_irfs(context.results.model_results[1].irfs, model, context.symboltable, filename) + end - end -function irfs2(GD, periods, context, solverWs::KOrderWs) +function irfs2(context, periods, solverWs::KOrderWs) model = context.models[1] results = context.results.model_results[1] + GD = results.solution_derivatives + endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)] exogenous_names = [Symbol(n) for n in get_exogenous_longname(context.symboltable)] From effbd74e677573a7477ed85f154ff89092b0ab0b Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Thu, 17 Aug 2023 14:52:20 -0400 Subject: [PATCH 08/10] Polish up IRFS 2nd order feature and integration in Dynare - remove any debugging code - unify the function signature between `irfs` and `irfs2` - call `irfs2` properly from `stoch_simul_core!` --- src/perturbations.jl | 83 +++++++++++++++++--------------------------- 1 file changed, 31 insertions(+), 52 deletions(-) diff --git a/src/perturbations.jl b/src/perturbations.jl index 953075a2..30156c7e 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -427,36 +427,28 @@ function stoch_simul_core!(context::Context, ws::DynamicWs, options::StochSimulO display_stoch_simul2(context, options) end end - if options.order == 1 && options.irf > 0 - exogenous_names = get_exogenous_longname(context.symboltable) - n = model.endogenous_nbr - m = model.exogenous_nbr - y = Dict{Symbol,AxisArrayTable} - irfs!(context, options.irf) + + if options.irf > 0 path = "$(context.modfileinfo.modfilepath)/graphs/" mkpath(path) filename = "$(path)/irfs" + + if options.order == 1 + irfs!(context, options.irf) + elseif options.order == 2 + irfs2!(context, options.irf) + end + plot_irfs( context.results.model_results[1].irfs, model, context.symboltable, filename, ) + display(context.results.model_results[1].irfs[:e]) end - # TODO: we should probably have a block like this - # else if options.order == 2 && options.irf > 0 - # irfs2(context, options.irf) - irfs2(context, options.irf, solverWs) - # path = "$(context.modfileinfo.modfilepath)/graphs/" - # mkpath(path) - # filename = "$(path)/irfs2" - # plot_irfs( - # context.results.model_results[1].irfs, - # model, - # context.symboltable, - # filename, - # ) - # end + + if (periods = options.periods) > 0 steadystate = results.trends.endogenous_steady_state linear_trend = results.trends.endogenous_linear_trend @@ -553,67 +545,54 @@ function compute_stoch_simul!( moments, order, solverWs) - # FOR DEBUGGING: single simulation - ut0 = sqrt(model.Sigma_e[1,1]) - KOrderPerturbations.simulate_run(G, ut0, 100, solverWs) copy!(results.solution_derivatives, G) - irfs2(context, 100, solverWs) - path = "$(context.modfileinfo.modfilepath)/graphs/" - filename = "$(path)/irfs2" - - # for development: delete all files that starts with "irfs2" - #filter(x -> startswith(x, "irfs2"), readdir("$(path)", join=true)) .|> x -> rm(x) - - display(context.results.model_results[1].irfs) - # plot_irfs(context.results.model_results[1].irfs, model, context.symboltable, filename) - end end -function irfs2(context, periods, solverWs::KOrderWs) +function irfs2!(context, periods) model = context.models[1] results = context.results.model_results[1] GD = results.solution_derivatives + exo_nbr = model.exogenous_nbr + state_index = model.i_bkwrd_b + endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)] exogenous_names = [Symbol(n) for n in get_exogenous_longname(context.symboltable)] - - nshock = solverWs.nshock - u_shock = [zeros(nshock) for _ in 1:periods] - u_shock_det = [zeros(nshock) for _ in 1:periods] gy1 = GD[1][:, 1] - y0 = zeros(size(gy1)[1]) - - n = length(y0) - simWs = SimulateWs(GD, n, solverWs) + n = size(gy1)[1] + y0 = zeros(n) - for exo_var in 1:model.exogenous_nbr + simWs = SimulateWs(GD, n, state_index, model.exogenous_nbr) + + for exo_var in 1:exo_nbr + u_shock = [zeros(exo_nbr) for _ in 1:periods] + sigma = diag(model.Sigma_e)[exo_var] det_shock = sqrt(sigma) delta_sims::Vector{Array{Array{Float64}}} = [] for n_sim in 1:100 - for i in 2:periods + for i in 1:periods random = randn() * det_shock^2 u_shock[i][exo_var] = random end - - u_shock_det = deepcopy(u_shock) - u_shock_det[2][exo_var] += det_shock - simWs = SimulateWs(GD, n, solverWs) rand_sim = simulate(GD, y0, u_shock, periods, simWs) + + u_shock_det = deepcopy(u_shock) + u_shock_det[1][exo_var] += det_shock rand_det_sim = simulate(GD, y0, u_shock_det, periods, simWs) delta_sim = rand_det_sim - rand_sim push!(delta_sims, delta_sim) end - tdf_vec = mean(delta_sims) - tdf_mat = reduce(vcat,transpose.(tdf_vec)) - tdf = AxisArrayTable(transpose(tdf_mat'), Undated(1):Undated(periods), endogenous_names) - + mean_irf = mean(delta_sims) + mean_irf_matrix = reduce(vcat, transpose.(mean_irf)) + + tdf = AxisArrayTable(mean_irf_matrix, Undated(1):Undated(periods), endogenous_names) results.irfs[exogenous_names[exo_var]] = tdf end From c4f49cdc7308edb1eed29c5df87ccccf3394d1c8 Mon Sep 17 00:00:00 2001 From: MichelJuillard Date: Fri, 18 Aug 2023 21:23:09 +0200 Subject: [PATCH 09/10] new version for irfs2 --- src/perturbations.jl | 60 +++++++++++++++---------------- test/models/example1/example1.mod | 3 +- 2 files changed, 32 insertions(+), 31 deletions(-) diff --git a/src/perturbations.jl b/src/perturbations.jl index 30156c7e..9d367a4b 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -408,7 +408,6 @@ function localapproximation!(; context::Context=context, ncol = m.n_bkwrd + m.n_current + m.n_fwrd + 2 * m.n_both tmp_nbr = m.dynamic_tmp_nbr ws = DynamicWs(context, order=options.order) - @show options.order stoch_simul_core!(context, ws, options) end @@ -427,7 +426,6 @@ function stoch_simul_core!(context::Context, ws::DynamicWs, options::StochSimulO display_stoch_simul2(context, options) end end - if options.irf > 0 path = "$(context.modfileinfo.modfilepath)/graphs/" mkpath(path) @@ -549,7 +547,7 @@ function compute_stoch_simul!( end end -function irfs2!(context, periods) +function irfs2!(context, periods; burning = 100) model = context.models[1] results = context.results.model_results[1] GD = results.solution_derivatives @@ -557,40 +555,38 @@ function irfs2!(context, periods) exo_nbr = model.exogenous_nbr state_index = model.i_bkwrd_b + model.Sigma_e + active_exogenous = findall(diag(model.Sigma_e) .> 0) + nx = length(active_exogenous) endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)] exogenous_names = [Symbol(n) for n in get_exogenous_longname(context.symboltable)] - + exogenous_names = view(exogenous_names, active_exogenous) + Sigma_e = view(model.Sigma_e, active_exogenous, active_exogenous) + C = transpose(cholesky(Sigma_e).U) gy1 = GD[1][:, 1] n = size(gy1)[1] y0 = zeros(n) simWs = SimulateWs(GD, n, state_index, model.exogenous_nbr) - for exo_var in 1:exo_nbr - u_shock = [zeros(exo_nbr) for _ in 1:periods] - - sigma = diag(model.Sigma_e)[exo_var] - det_shock = sqrt(sigma) + for (i, exo_var) in enumerate(active_exogenous) + u_shock = [zeros(exo_nbr) for _ in 1:periods + burning] + det_shock = C[:, i] - delta_sims::Vector{Array{Array{Float64}}} = [] - for n_sim in 1:100 - for i in 1:periods - random = randn() * det_shock^2 - u_shock[i][exo_var] = random + mean_sims = [zeros(model.endogenous_nbr) for _ in 1:periods] + replic = 1000 + for n_sim in 1:replic + for i in 1:periods + burning + u_shock[i][active_exogenous] = C*randn(nx) end - rand_sim = simulate(GD, y0, u_shock, periods, simWs) - + rand_sim = simulate(GD, y0, u_shock, periods + burning, simWs) u_shock_det = deepcopy(u_shock) - u_shock_det[1][exo_var] += det_shock - rand_det_sim = simulate(GD, y0, u_shock_det, periods, simWs) - - delta_sim = rand_det_sim - rand_sim - - push!(delta_sims, delta_sim) + @views u_shock_det[burning + 1] .+= det_shock + rand_det_sim = simulate(GD, y0, u_shock_det, periods + burning, simWs) + @views mean_sims .+= (rand_det_sim[burning + 1:end] .- rand_sim[burning + 1:end])/replic end - - mean_irf = mean(delta_sims) - mean_irf_matrix = reduce(vcat, transpose.(mean_irf)) + + mean_irf_matrix = reduce(vcat, transpose.(mean_sims)) tdf = AxisArrayTable(mean_irf_matrix, Undated(1):Undated(periods), endogenous_names) results.irfs[exogenous_names[exo_var]] = tdf @@ -637,19 +633,23 @@ end function irfs!(context, periods) model = context.models[1] + model.Sigma_e results = context.results.model_results[1] + active_exogenous = findall(diag(model.Sigma_e) .> 0) endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)] exogenous_names = [Symbol(n) for n in get_exogenous_longname(context.symboltable)] - C = cholesky(model.Sigma_e + 1e-14 * I) - x = zeros(model.exogenous_nbr) + exogenous_names = view(exogenous_names, active_exogenous) + Sigma_e = view(model.Sigma_e, active_exogenous, active_exogenous) + C = transpose(cholesky(model.Sigma_e).U) + x = zeros(length(active_exogenous)) A = zeros(model.endogenous_nbr, model.endogenous_nbr) B = zeros(model.endogenous_nbr, model.exogenous_nbr) make_A_B!(A, B, model, results) - for i = 1:model.exogenous_nbr + for (i, j) in enumerate(active_exogenous) fill!(x, 0) - x[i] = robustsqrt(model.Sigma_e[i, i]) + @views x .= C[:, j] yy = Matrix{Float64}(undef, size(A, 1), periods) - mul!(view(yy, :, 1), B, x) + @views mul!(view(yy, :, 1), B[:, active_exogenous], x) for j = 2:periods mul!(view(yy, :, j), A, view(yy, :, j - 1)) end diff --git a/test/models/example1/example1.mod b/test/models/example1/example1.mod index 0af7773b..def8fdf4 100644 --- a/test/models/example1/example1.mod +++ b/test/models/example1/example1.mod @@ -47,12 +47,13 @@ shocks; var e; stderr 0.009; var u; stderr 0.009; //var e, u = phi*0.009*0.009; -var e, u = 0.009*0.009; +var e, u = 0.1*0.009*0.009; end; check; stoch_simul(dr=cycle_reduction, order=1, irf=100); +irfs1 = deepcopy(context.results.model_results[1].irfs) stoch_simul(dr=cycle_reduction, order=2, irf=100); From 6f7c63e5b128c0a46184136e99fe7182ee9332fa Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Wed, 23 Aug 2023 15:46:57 -0400 Subject: [PATCH 10/10] add display of moments and correlation statistics at second order except of function names, most of this code is probably generalization to k-order solutions --- src/perturbations.jl | 76 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/src/perturbations.jl b/src/perturbations.jl index 9d367a4b..ae361edc 100644 --- a/src/perturbations.jl +++ b/src/perturbations.jl @@ -77,6 +77,41 @@ function display_stoch_simul2(context::Context, options::StochSimulOptions) LRE_results = results.linearrationalexpectations stationary_variables = LRE_results.stationary_variables display_solution_function2(results.solution_derivatives, endogenous_names, exogenous_names, m) + simulation_results = long_second_order_simulation(context) + display_mean_sd_variance2(simulation_results, endogenous_names, m) + display_correlation2(simulation_results, endogenous_names, m) + display_autocorrelation2(simulation_results, endogenous_names, m, options) +end + +function long_second_order_simulation(context; periods = 100_000, burning = 100) + model = context.models[1] + results = context.results.model_results[1] + GD = results.solution_derivatives + + exo_nbr = model.exogenous_nbr + original_endo_nbr = model.original_endogenous_nbr + state_index = model.i_bkwrd_b + + gy1 = GD[1][:, 1] + y0 = zeros(size(gy1, 1)) + + active_exogenous = findall(diag(model.Sigma_e) .> 0) + Sigma_e = view(model.Sigma_e, active_exogenous, active_exogenous) + C = transpose(cholesky(Sigma_e).U) + n_active = length(active_exogenous) + + u_shock = [zeros(exo_nbr) for _ in 1:periods] + for i in 1:periods + u_shock[i][active_exogenous] = C*randn(n_active) + end + + simWs = SimulateWs(GD, length(y0), state_index, model.exogenous_nbr) + simulation_result_vec = simulate(GD, y0, u_shock, periods, simWs) + simulation_result_mat = stack(simulation_result_vec)' + + # ignore burning periods and useless endogenous variables + simulation_result_clean = simulation_result_mat[burning:end, 1:original_endo_nbr] + return simulation_result_clean end function display_solution_function( @@ -223,6 +258,20 @@ function display_mean_sd_variance( dynare_table(data, title) end +function display_mean_sd_variance2(sim_results, endogenous_names, model::Model) + original_endo_nbr = model.original_endogenous_nbr + + title = "SIMULATED MOMENTS" + data = Matrix{Any}(undef, original_endo_nbr + 1, 4) + data[1, :] = ["VARIABLE", "MEAN", "STD. DEV.", "VARIANCE"] + data[2:end, 1] = endogenous_names[1:original_endo_nbr] + data[2:end, 2] .= map(mean, eachcol(sim_results)) + data[2:end, 3] .= map(std, eachcol(sim_results)) + data[2:end, 4] .= map(var, eachcol(sim_results)) + println("\n") + dynare_table(data, title) +end + function display_variance_decomposition( LREresults::LinearRationalExpectationsResults, endogenous_names::AbstractVector{String}, @@ -318,6 +367,19 @@ function display_correlation( dynare_table(data, title) end +function display_correlation2(sim_results, endogenous_names, model::Model) + original_endo_nbr = model.original_endogenous_nbr + + title = "SIMULATED CORRELATION MATRIX" + data = Matrix{Any}(undef, original_endo_nbr + 1, original_endo_nbr + 1) + data[1, 1] = "" + data[1, 2:end] = endogenous_names[1:original_endo_nbr] |> permutedims + data[2:end, 1] = endogenous_names[1:original_endo_nbr] + data[2:end, 2:end] = cor(sim_results) + println("\n") + dynare_table(data, title) +end + function display_autocorrelation( LREresults::LinearRationalExpectationsResults, endogenous_names::AbstractVector{String}, @@ -370,6 +432,20 @@ function display_autocorrelation( dynare_table(data, title) end +function display_autocorrelation2(sim_results, endogenous_names, model, options::StochSimulOptions) + original_endo_nbr = model.original_endogenous_nbr + lags = [i for i in 1:options.nar] + + title = "SIMULATED AUTOCORRELATION COEFFICIENTS" + data = Matrix{Any}(undef, original_endo_nbr + 1, options.nar + 1) + data[1, 1] = "" + data[2:end, 1] = endogenous_names[1:original_endo_nbr] + data[1, 2:end] = lags + data[2:end, 2:end] = autocor(sim_results, lags)' + println("\n") + dynare_table(data, title) +end + function make_A_B!( A::Matrix{Float64}, B::Matrix{Float64},