diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 6e474788..4a199c10 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -1365,9 +1365,14 @@ function calculate_jacobian(parameters::Vector{<: Number}, 𝓂::ℳ) par = ComponentVector( vcat(parameters,calibrated_parameters),Axis(vcat(𝓂.parameters,𝓂.calibration_equations_parameters))) SS = ComponentVector(non_stochastic_steady_state, Axis(sort(union(𝓂.exo_present,𝓂.var)))) - SS_past = SS[[indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_past))...] - SS_present = SS[[indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_present))...] - SS_future = SS[[indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_future))...] + past_idx = [indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_past = length(past_idx) > 0 ? SS[past_idx] : zeros(0) #; zeros(length(𝓂.exo_past))...] + + present_idx = [indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_present = length(present_idx) > 0 ? SS[present_idx] : zeros(0)#; zeros(length(𝓂.exo_present))...] + + future_idx = [indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_future = length(future_idx) > 0 ? SS[future_idx] : zeros(0)#; zeros(length(𝓂.exo_future))...] shocks_ss = zeros(length(𝓂.exo)) @@ -1383,10 +1388,15 @@ function calculate_hessian(parameters::Vector{<: Number}, 𝓂::ℳ) par = ComponentVector( vcat(parameters,calibrated_parameters),Axis(vcat(𝓂.parameters,𝓂.calibration_equations_parameters))) SS = ComponentVector(non_stochastic_steady_state, Axis(𝓂.var)) - - SS_past = SS[[indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_past))...] - SS_present = SS[[indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_present))...] - SS_future = SS[[indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_future))...] + + past_idx = [indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_past = length(past_idx) > 0 ? SS[past_idx] : zeros(0) #; zeros(length(𝓂.exo_past))...] + + present_idx = [indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_present = length(present_idx) > 0 ? SS[present_idx] : zeros(0)#; zeros(length(𝓂.exo_present))...] + + future_idx = [indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_future = length(future_idx) > 0 ? SS[future_idx] : zeros(0)#; zeros(length(𝓂.exo_future))...] shocks_ss = zeros(length(𝓂.exo)) @@ -1404,10 +1414,15 @@ function calculate_third_order_derivatives(parameters::Vector{<: Number}, 𝓂:: par = ComponentVector( vcat(parameters,calibrated_parameters),Axis(vcat(𝓂.parameters,𝓂.calibration_equations_parameters))) SS = ComponentVector(non_stochastic_steady_state, Axis(𝓂.var)) - - SS_past = SS[[indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_past))...] - SS_present = SS[[indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_present))...] - SS_future = SS[[indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...]]#; zeros(length(𝓂.exo_future))...] + + past_idx = [indexin(sort([𝓂.var_past; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_past,𝓂.exo_past))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_past = length(past_idx) > 0 ? SS[past_idx] : zeros(0) #; zeros(length(𝓂.exo_past))...] + + present_idx = [indexin(sort([𝓂.var_present; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_present,𝓂.exo_present))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_present = length(present_idx) > 0 ? SS[present_idx] : zeros(0)#; zeros(length(𝓂.exo_present))...] + + future_idx = [indexin(sort([𝓂.var_future; map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾|ᴸ⁽[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), union(𝓂.aux_future,𝓂.exo_future))]), sort(union(𝓂.var,𝓂.exo_present)))...] + SS_future = length(future_idx) > 0 ? SS[future_idx] : zeros(0)#; zeros(length(𝓂.exo_future))...] shocks_ss = zeros(length(𝓂.exo)) @@ -1868,10 +1883,11 @@ function irf(state_update::Function, initial_state::Vector{Float64}, T::timings; return KeyedArray(Y[var_idx,:,:]; Variables = T.var[var_idx], Period = 1:periods, Shock = [:simulate]) elseif shocks == :none Y = zeros(T.nVars,periods,1) - Y[:,1,1] = state_update(initial_state,[0.0]) + shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : [0.0] + Y[:,1,1] = state_update(initial_state,shck) for t in 1:periods-1 - Y[:,t+1,1] = state_update(Y[:,t,1],[0.0]) + Y[:,t+1,1] = state_update(Y[:,t,1],shck) end return KeyedArray(Y[var_idx,:,:]; Variables = T.var[var_idx], Period = 1:periods, Shock = [:none]) @@ -1993,7 +2009,7 @@ function parse_shocks_input_to_index(shocks::Symbol_input, T::timings) if shocks == :all shock_idx = 1:T.nExo elseif shocks == :none - shock_idx = 1:T.nExo + shock_idx = 1 elseif shocks == :simulate shock_idx = 1 elseif shocks isa Matrix{Symbol} diff --git a/src/get_functions.jl b/src/get_functions.jl index 277c2a4b..017d9743 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -62,6 +62,8 @@ function get_irf(𝓂::ℳ, sol_mat = calculate_first_order_solution(jacc; T = 𝓂.timings) state_update = function(state::Vector, shock::Vector) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end + + shocks = 𝓂.timings.nExo == 0 ? :none : shocks shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) @@ -181,6 +183,12 @@ function get_irf(𝓂::ℳ; var_idx = parse_variables_input_to_index(variables, 𝓂.timings) end + shocks = 𝓂.timings.nExo == 0 ? :none : shocks + + if shocks == :none && generalised_irf + @error "Cannot compute generalised IRFs for model without shocks." + end + if generalised_irf girfs = girf(state_update, 𝓂.timings; diff --git a/src/plotting.jl b/src/plotting.jl index a565f0a6..e29eb154 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -76,6 +76,8 @@ function plot(𝓂::ℳ; end init_state = initial_state == [0.0] ? zeros(𝓂.timings.nVars) : initial_state - collect(get_non_stochastic_steady_state_internal(𝓂)) + + shocks = 𝓂.timings.nExo == 0 ? :none : shocks shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) @@ -206,7 +208,6 @@ function plot(𝓂::ℳ; # elseif length(pp) > 1 if length(pp) > 0 - shock_string = ": " * string(𝓂.timings.exo[shock_idx[shock]]) if shocks == :simulate shock_string = ": simulate all" @@ -215,6 +216,7 @@ function plot(𝓂::ℳ; shock_string = "" shock_name = "no_shock" else + shock_string = ": " * string(𝓂.timings.exo[shock_idx[shock]]) shock_name = string(𝓂.timings.exo[shock_idx[shock]]) end diff --git a/test/runtests.jl b/test/runtests.jl index fa47aaef..c66f91d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,34 @@ using Random include("test_standalone_function.jl") end +@testset "Model without shocks" begin + @model m begin + K[0] = (1 - δ) * K[-1] + I[0] + Z[0] = (1 - ρ) * μ + ρ * Z[-1] + I[1] = ((ρ + δ - Z[0])/(1 - δ)) + ((1 + ρ)/(1 - δ)) * I[0] + end + + @parameters m begin + ρ = 0.05 + δ = 0.10 + μ = .17 + σ = .2 + end + + m_ss = get_steady_state(m) + @test isapprox(m_ss(:,:Steady_state),[1/7.5,1/.75,.17],rtol = eps(Float32)) + + m_sol = get_solution(m) + @test isapprox(m_sol(:,:K),[1/.75,.9,.04975124378109454],rtol = eps(Float32)) + + init = m_ss(:,:Steady_state) |> collect + init[2] *= 1.5 + get_irf(m, initial_state = init, shocks = :none) + + plot_irf(m, initial_state = init, shocks = :none) + @test true +end + @testset "Distribution functions, general and SS" begin @model RBC_CME begin