From 6ebfe20b3aae0670d5fabb668c5d9b9400f81791 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Tue, 23 Apr 2024 19:22:28 +0200 Subject: [PATCH 1/7] first partial attempt (incomplete) --- docs/src/unfinished_docs/todo.md | 2 + test/test_break_points.jl | 330 +++++++++++++++++++++++++++++++ 2 files changed, 332 insertions(+) create mode 100644 test/test_break_points.jl diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md index 771821728..7c9150fe2 100644 --- a/docs/src/unfinished_docs/todo.md +++ b/docs/src/unfinished_docs/todo.md @@ -4,6 +4,8 @@ - [ ] ss transition by entering new parameters at given periods - [ ] check downgrade tests +- [ ] break point estimation, obc, shock decomp +- [ ] higher ordshock decomps - [ ] figure out why PG and IS return basically the prior - [ ] allow external functions to calculate the steady state (and hand it over via SS or get_loglikelihood function) - need to use the check function for implicit derivatives and cannot use it to get him a guess from which he can use internal solver going forward - [ ] go through custom SS solver once more and try to find parameters and logic that achieves best results diff --git a/test/test_break_points.jl b/test/test_break_points.jl new file mode 100644 index 000000000..1eb9b7f61 --- /dev/null +++ b/test/test_break_points.jl @@ -0,0 +1,330 @@ +using MacroModelling + +include("../models/Smets_Wouters_2007.jl") + +irfs = get_irf(Smets_Wouters_2007) + +break_points = KeyedArray([.4,.3,.25]', Variable =[:ctrend], Time = [1,15,30]) + + +# translate this to the combined parameter vector +SS() + + +import MacroModelling: get_relevant_steady_states, parse_shocks_input_to_index, parse_algorithm_to_state_update, obc_objective_optim_fun, obc_constraint_optim_fun, String_input, Symbol_input, timings + +𝓂 = Smets_Wouters_2007 +periods = 40 +algorithm = :first_order + +parameters = nothing + +variables = :all_excluding_obc +shocks = :all_excluding_obc +negative_shock = false +generalised_irf = false +initial_state = [0.0] +levels = false +ignore_obc = false +verbose = false + +shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks + +shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks + +shocks = 𝓂.timings.nExo == 0 ? :none : shocks + +@assert !(shocks == :none && generalised_irf) "Cannot compute generalised IRFs for model without shocks." + +stochastic_model = length(𝓂.timings.exo) > 0 + +obc_model = length(𝓂.obc_violation_equations) > 0 + +if shocks isa Matrix{Float64} + @assert size(shocks)[1] == 𝓂.timings.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." + + periods += size(shocks)[2] + + shock_history = zeros(𝓂.timings.nExo, periods) + + shock_history[:,1:size(shocks)[2]] = shocks + + shock_idx = 1 + + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:]) > 1e-10 +elseif shocks isa KeyedArray{Float64} + shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) + + periods += size(shocks)[2] + + @assert length(setdiff(shock_input, 𝓂.timings.exo)) == 0 "Provided shocks which are not part of the model." + + shock_history = zeros(𝓂.timings.nExo, periods + 1) + + shock_history[indexin(shock_input,𝓂.timings.exo),1:size(shocks)[2]] = shocks + + shock_idx = 1 + + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks(intersect(𝓂.timings.exo,axiskeys(shocks,1)),:)) > 1e-10 +else + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) +end + +if ignore_obc + occasionally_binding_constraints = false +else + occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 +end + + +irfs = irf(state_update, + initial_state, + levels ? reference_steady_state + SSS_delta : SSS_delta, + 𝓂.timings; + periods = periods, + shocks = shocks, + variables = variables, + negative_shock = negative_shock) + + + + +var_idx = parse_variables_input_to_index(variables, T) + +axis1 = T.var[var_idx] + +if any(x -> contains(string(x), "β—–"), axis1) + axis1_decomposed = decompose_name.(axis1) + axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] +end + + + + +if shocks == :simulate + axis2 = [:simulate] +elseif shocks == :none + axis2 = [:none] +else + axis2 = shocks isa Union{Symbol_input,String_input} ? + shock_idx isa Int ? + [T.exo[shock_idx]] : + T.exo[shock_idx] : + [:Shock_matrix] + + if any(x -> contains(string(x), "β—–"), axis2) + axis2_decomposed = decompose_name.(axis2) + axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] + end +end + + +function irff(state_update::Function, + initial_state::Union{Vector{Vector{Float64}},Vector{Float64},Matrix{Float64}}, + level::Vector{Float64}, + T::timings; + periods::Int = 40, + shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, + variables::Union{Symbol_input,String_input} = :all, + negative_shock::Bool = false) + + pruning = initial_state isa Vector{Vector{Float64}} + + shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks + + shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks + + if shocks isa Matrix{Float64} + @assert size(shocks)[1] == T.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." + + # periods += size(shocks)[2] + + shock_history = zeros(T.nExo, periods) + + shock_history[:,1:size(shocks)[2]] = shocks + + shock_idx = 1 + elseif shocks isa KeyedArray{Float64} + shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) + + # periods += size(shocks)[2] + + @assert length(setdiff(shock_input, T.exo)) == 0 "Provided shocks which are not part of the model." + + shock_history = zeros(T.nExo, periods) + + shock_history[indexin(shock_input,T.exo),1:size(shocks)[2]] = shocks + + shock_idx = 1 + else + shock_idx = parse_shocks_input_to_index(shocks,T) + end + + + if shocks == :simulate + shock_history = randn(T.nExo,periods) + + shock_history[contains.(string.(T.exo),"α΅’α΅‡αΆœ"),:] .= 0 + + Y = zeros(T.nVars,periods,1) + + initial_state = state_update(initial_state,shock_history[:,1]) + + Y[:,1,1] = pruning ? sum(initial_state) : initial_state + + for t in 1:periods-1 + initial_state = state_update(initial_state,shock_history[:,t+1]) + + Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state + end + elseif shocks == :none + Y = zeros(T.nVars,periods,1) + + shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : zeros(T.nExo) + + initial_state = state_update(initial_state, shck) + + Y[:,1,1] = pruning ? sum(initial_state) : initial_state + + for t in 1:periods-1 + initial_state = state_update(initial_state, shck) + + Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state + end + else + Y = zeros(T.nVars,periods,length(shock_idx)) + + for (i,ii) in enumerate(shock_idx) + initial_state_copy = deepcopy(initial_state) + + if shocks != :simulate && shocks isa Union{Symbol_input,String_input} + shock_history = zeros(T.nExo,periods) + shock_history[ii,1] = negative_shock ? -1 : 1 + end + + initial_state_copy = state_update(initial_state_copy, shock_history[:,1]) + + Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + + for t in 1:periods-1 + initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1]) + + Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + end + end + end + + return Y +end + + + +break_points +break_points = KeyedArray([nothing .3 .25 +.77 NaN .9], Variable =[:ctrend, :constepinf], Time = [1,15,30]) + + +function transform_break_points(break_points) + dict = Dict{Int, Dict{Symbol, Float64}}() + + for i in 1:size(break_points,2) + time = axiskeys(break_points,2)[i] + + if !haskey(dict, time) + dict[time] = Dict{Symbol, Float64}() + end + + for k in 1:size(break_points,1) + parameter = axiskeys(break_points,1)[k] + value = break_points[k,i] + if !(isnothing(value) || !isfinite(value)) + dict[time][parameter] = value + end + end + end + + # if !haskey(dict, 0) + # dict[0] = Dict{Symbol, Float64}() + # end + + return dict +end + +break_points_dict = transform_break_points(break_points) +break_points_dict + +previous_parameter_values = deepcopy(𝓂.parameter_values) + +periods_of_change = break_points_dict |> keys |> collect |> sort + + +solve!(𝓂, parameters = break_points_dict[0], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + + +for (i,p) in enumerate(periods_of_change) + concerned_periods = (i == 1 ? 1 : periods_of_change[i-1]):periods_of_change[i] + periods_with_these_parameters = periods_of_change[i] - (i == 1 ? 1 : periods_of_change[i-1]) + + solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + + reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + + unspecified_initial_state = initial_state == [0.0] + + if unspecified_initial_state + if algorithm == :pruned_second_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + else + initial_state = zeros(𝓂.timings.nVars) - SSS_delta + end + else + if initial_state isa Vector{Float64} + if algorithm == :pruned_second_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + else + initial_state = initial_state - NSSS + end + else + if algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] + @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." + end + end + end + + if occasionally_binding_constraints + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + elseif obc_shocks_included + @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." + + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + else + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + end + + irfs1 = irff(state_update, + initial_state, + levels ? reference_steady_state + SSS_delta : SSS_delta, + 𝓂.timings; + periods = p, + shocks = shocks, + variables = variables, + negative_shock = negative_shock) + +end + +irfs1 = irff(state_update, + initial_state, + levels ? reference_steady_state + SSS_delta : SSS_delta, + 𝓂.timings; + periods = periods, + shocks = shocks, + variables = variables, + negative_shock = negative_shock) + +irfs1[:,1,:] \ No newline at end of file From f70e5ebaee3fa6b3920071663c5f215a7654285a Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 24 Apr 2024 14:11:54 +0200 Subject: [PATCH 2/7] Update Turing version in Project.toml to 0.30, 0.31 --- Project.toml | 2 +- docs/src/unfinished_docs/todo.md | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 59f1b2c70..6fbb17e62 100644 --- a/Project.toml +++ b/Project.toml @@ -93,7 +93,7 @@ SymPyPythonCall = "^0.2" Symbolics = "^5" Test = "^1" ThreadedSparseArrays = "^0.2.3" -Turing = "0.29, 0.30" +Turing = "0.30, 0.31" Unicode = "^1" Zygote = "^0.6" julia = "1.8" diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md index 7c9150fe2..3d340115f 100644 --- a/docs/src/unfinished_docs/todo.md +++ b/docs/src/unfinished_docs/todo.md @@ -4,6 +4,7 @@ - [ ] ss transition by entering new parameters at given periods - [ ] check downgrade tests +- [ ] obc and conditional forecasting should be the same problem. see if you can use the simple matrix linear algebra from cond fcst in obc - [ ] break point estimation, obc, shock decomp - [ ] higher ordshock decomps - [ ] figure out why PG and IS return basically the prior From 9c27e4ecb8c1dcf471a3ea9a7898d6319c3fda12 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 24 Apr 2024 17:07:01 +0200 Subject: [PATCH 3/7] first working prototype --- src/MacroModelling.jl | 310 +++++++++------------ src/get_functions.jl | 146 +++++++--- test/test_break_points.jl | 560 +++++++++++++++++++++++++++++--------- 3 files changed, 670 insertions(+), 346 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 5d9982e0b..24ccf1ca7 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -393,8 +393,33 @@ end # res .= abs2.(conditions[cond_var_idx] - state_update(state, convert(typeof(X), shocks))[cond_var_idx]) # end +transform_break_points(break_points::Dict{Int, Dict{Symbol, Float64}}) = break_points +function transform_break_points(break_points::KeyedArray{Float64})::Dict{Int, Dict{Symbol, Float64}} + dict = Dict{Int, Dict{Symbol, Float64}}() + + for i in 1:size(break_points,2) + time = axiskeys(break_points,2)[i] + + if !haskey(dict, time) + dict[time] = Dict{Symbol, Float64}() + end + + for k in 1:size(break_points,1) + parameter = axiskeys(break_points,1)[k] + value = break_points[k,i] + if !(isnothing(value) || !isfinite(value)) + dict[time][parameter] = value + end + end + end + + if !haskey(dict, 0) + dict[0] = Dict{Symbol, Float64}() + end + return dict +end function minimize_distance_to_initial_data!(X::Vector{S}, grad::Vector{S}, data::Vector{T}, state::Union{Vector{T},Vector{Vector{T}}}, state_update::Function, warmup_iters::Int, cond_var_idx::Vector{Union{Nothing, Int64}}, precision_factor::Float64) where {S, T} @@ -6299,125 +6324,54 @@ end -function irf(state_update::Function, - initial_state::Union{Vector{Vector{Float64}},Vector{Float64}}, - level::Vector{Float64}, - T::timings; - periods::Int = 40, - shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, - variables::Union{Symbol_input,String_input} = :all, - negative_shock::Bool = false) - - pruning = initial_state isa Vector{Vector{Float64}} - - shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks - - shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks - - if shocks isa Matrix{Float64} - @assert size(shocks)[1] == T.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." - - # periods += size(shocks)[2] - - shock_history = zeros(T.nExo, periods) - - shock_history[:,1:size(shocks)[2]] = shocks - - shock_idx = 1 - elseif shocks isa KeyedArray{Float64} - shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) - - # periods += size(shocks)[2] - - @assert length(setdiff(shock_input, T.exo)) == 0 "Provided shocks which are not part of the model." - - shock_history = zeros(T.nExo, periods) - - shock_history[indexin(shock_input,T.exo),1:size(shocks)[2]] = shocks - - shock_idx = 1 - else - shock_idx = parse_shocks_input_to_index(shocks,T) - end - - var_idx = parse_variables_input_to_index(variables, T) - - axis1 = T.var[var_idx] - - if any(x -> contains(string(x), "β—–"), axis1) - axis1_decomposed = decompose_name.(axis1) - axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] - end - - if shocks == :simulate - shock_history = randn(T.nExo,periods) - - shock_history[contains.(string.(T.exo),"α΅’α΅‡αΆœ"),:] .= 0 - - Y = zeros(T.nVars,periods,1) - - initial_state = state_update(initial_state,shock_history[:,1]) - - Y[:,1,1] = pruning ? sum(initial_state) : initial_state +# function irf(steady_states_and_state_update::Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}, +# initial_state::Union{Vector{Vector{Float64}},Vector{Float64}}, +# shock_history::Matrix{Float64}, +# T::timings) +function irf(steady_states_and_state_update::Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}, + initial_state::Union{Vector{Vector{Float64}},Vector{Float64}}, + shock_history::Array{Float64,3}, + shock_idx::Vector{Int}, + levels::Bool, + T::timings) - for t in 1:periods-1 - initial_state = state_update(initial_state,shock_history[:,t+1]) + periods = size(shock_history, 2) - Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state - end + Y = zeros(T.nVars, periods, length(shock_idx)) - return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:simulate]) - elseif shocks == :none - Y = zeros(T.nVars,periods,1) + pruning = initial_state isa Vector{Vector{Float64}} - shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : zeros(T.nExo) + periods_of_change = steady_states_and_state_update |> keys |> collect |> sort - initial_state = state_update(initial_state, shck) + initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] - Y[:,1,1] = pruning ? sum(initial_state) : initial_state + for i in shock_idx + for (k,p) in enumerate(periods_of_change) + concerned_periods = (k == 1 ? 1 : periods_of_change[k]):(k == length(periods_of_change) ? periods : periods_of_change[k+1] - 1) - for t in 1:periods-1 - initial_state = state_update(initial_state, shck) + reference_steady_state, __, SSS_delta, state_update = steady_states_and_state_update[p] - Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state - end + if k > 1 + Ξ”SS = steady_states_and_state_update[periods_of_change[k]][2] - steady_states_and_state_update[periods_of_change[k-1]][2] - return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) - else - Y = zeros(T.nVars,periods,length(shock_idx)) - - for (i,ii) in enumerate(shock_idx) - initial_state_copy = deepcopy(initial_state) - - if shocks != :simulate && shocks isa Union{Symbol_input,String_input} - shock_history = zeros(T.nExo,periods) - shock_history[ii,1] = negative_shock ? -1 : 1 + if pruning + for j in initial_state_copy[i] + j += Ξ”SS + end + else + initial_state_copy[i] += Ξ”SS + end end - initial_state_copy = state_update(initial_state_copy, shock_history[:,1]) - - Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + for t in concerned_periods + initial_state_copy[i] = state_update(initial_state_copy[i], shock_history[:,t,i]) - for t in 1:periods-1 - initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1]) - - Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + Y[:,t,i] = pruning ? sum(initial_state_copy[i]) : initial_state_copy[i] .+ (levels ? reference_steady_state + SSS_delta : SSS_delta) end end - - axis2 = shocks isa Union{Symbol_input,String_input} ? - shock_idx isa Int ? - [T.exo[shock_idx]] : - T.exo[shock_idx] : - [:Shock_matrix] - - if any(x -> contains(string(x), "β—–"), axis2) - axis2_decomposed = decompose_name.(axis2) - axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] - end - - return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) end + + return Y end @@ -6589,9 +6543,9 @@ function parse_shocks_input_to_index(shocks::Union{Symbol_input,String_input}, T elseif shocks == :all_excluding_obc shock_idx = findall(.!contains.(string.(T.exo),"α΅’α΅‡αΆœ")) elseif shocks == :none - shock_idx = 1 + shock_idx = [1] elseif shocks == :simulate - shock_idx = 1 + shock_idx = [1] elseif shocks isa Matrix{Symbol} if length(setdiff(shocks,T.exo)) > 0 @warn "Following shocks are not part of the model: " * join(string.(setdiff(shocks,T.exo)),", ") @@ -7964,80 +7918,80 @@ function filter_and_smooth(𝓂::β„³, end -if VERSION >= v"1.9" - @setup_workload begin - # Putting some things in `setup` can reduce the size of the - # precompile file and potentially make loading faster. - @model FS2000 precompile = true begin - dA[0] = exp(gam + z_e_a * e_a[x]) - log(m[0]) = (1 - rho) * log(mst) + rho * log(m[-1]) + z_e_m * e_m[x] - - P[0] / (c[1] * P[1] * m[0]) + bet * P[1] * (alp * exp( - alp * (gam + log(e[1]))) * k[0] ^ (alp - 1) * n[1] ^ (1 - alp) + (1 - del) * exp( - (gam + log(e[1])))) / (c[2] * P[2] * m[1])=0 - W[0] = l[0] / n[0] - - (psi / (1 - psi)) * (c[0] * P[0] / (1 - n[0])) + l[0] / n[0] = 0 - R[0] = P[0] * (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ ( - alp) / W[0] - 1 / (c[0] * P[0]) - bet * P[0] * (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ (1 - alp) / (m[0] * l[0] * c[1] * P[1]) = 0 - c[0] + k[0] = exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ (1 - alp) + (1 - del) * exp( - (gam + z_e_a * e_a[x])) * k[-1] - P[0] * c[0] = m[0] - m[0] - 1 + d[0] = l[0] - e[0] = exp(z_e_a * e_a[x]) - y[0] = k[-1] ^ alp * n[0] ^ (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) - gy_obs[0] = dA[0] * y[0] / y[-1] - gp_obs[0] = (P[0] / P[-1]) * m[-1] / dA[0] - log_gy_obs[0] = log(gy_obs[0]) - log_gp_obs[0] = log(gp_obs[0]) - end - - @parameters FS2000 silent = true precompile = true begin - alp = 0.356 - bet = 0.993 - gam = 0.0085 - mst = 1.0002 - rho = 0.129 - psi = 0.65 - del = 0.01 - z_e_a = 0.035449 - z_e_m = 0.008862 - end - - ENV["GKSwstype"] = "nul" - - @compile_workload begin - # all calls in this block will be precompiled, regardless of whether - # they belong to your package or not (on Julia 1.8 and higher) - @model RBC precompile = true begin - 1 / c[0] = (0.95 / c[1]) * (Ξ± * exp(z[1]) * k[0]^(Ξ± - 1) + (1 - Ξ΄)) - c[0] + k[0] = (1 - Ξ΄) * k[-1] + exp(z[0]) * k[-1]^Ξ± - z[0] = 0.2 * z[-1] + 0.01 * eps_z[x] - end - - @parameters RBC silent = true precompile = true begin - Ξ΄ = 0.02 - Ξ± = 0.5 - end +# if VERSION >= v"1.9" +# @setup_workload begin +# # Putting some things in `setup` can reduce the size of the +# # precompile file and potentially make loading faster. +# @model FS2000 precompile = true begin +# dA[0] = exp(gam + z_e_a * e_a[x]) +# log(m[0]) = (1 - rho) * log(mst) + rho * log(m[-1]) + z_e_m * e_m[x] +# - P[0] / (c[1] * P[1] * m[0]) + bet * P[1] * (alp * exp( - alp * (gam + log(e[1]))) * k[0] ^ (alp - 1) * n[1] ^ (1 - alp) + (1 - del) * exp( - (gam + log(e[1])))) / (c[2] * P[2] * m[1])=0 +# W[0] = l[0] / n[0] +# - (psi / (1 - psi)) * (c[0] * P[0] / (1 - n[0])) + l[0] / n[0] = 0 +# R[0] = P[0] * (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ ( - alp) / W[0] +# 1 / (c[0] * P[0]) - bet * P[0] * (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ (1 - alp) / (m[0] * l[0] * c[1] * P[1]) = 0 +# c[0] + k[0] = exp( - alp * (gam + z_e_a * e_a[x])) * k[-1] ^ alp * n[0] ^ (1 - alp) + (1 - del) * exp( - (gam + z_e_a * e_a[x])) * k[-1] +# P[0] * c[0] = m[0] +# m[0] - 1 + d[0] = l[0] +# e[0] = exp(z_e_a * e_a[x]) +# y[0] = k[-1] ^ alp * n[0] ^ (1 - alp) * exp( - alp * (gam + z_e_a * e_a[x])) +# gy_obs[0] = dA[0] * y[0] / y[-1] +# gp_obs[0] = (P[0] / P[-1]) * m[-1] / dA[0] +# log_gy_obs[0] = log(gy_obs[0]) +# log_gp_obs[0] = log(gp_obs[0]) +# end - get_SS(FS2000) - get_SS(FS2000, parameters = :alp => 0.36) - get_solution(FS2000) - get_solution(FS2000, parameters = :alp => 0.35) - get_standard_deviation(FS2000) - get_correlation(FS2000) - get_autocorrelation(FS2000) - get_variance_decomposition(FS2000) - get_conditional_variance_decomposition(FS2000) - get_irf(FS2000) +# @parameters FS2000 silent = true precompile = true begin +# alp = 0.356 +# bet = 0.993 +# gam = 0.0085 +# mst = 1.0002 +# rho = 0.129 +# psi = 0.65 +# del = 0.01 +# z_e_a = 0.035449 +# z_e_m = 0.008862 +# end + +# ENV["GKSwstype"] = "nul" + +# @compile_workload begin +# # all calls in this block will be precompiled, regardless of whether +# # they belong to your package or not (on Julia 1.8 and higher) +# @model RBC precompile = true begin +# 1 / c[0] = (0.95 / c[1]) * (Ξ± * exp(z[1]) * k[0]^(Ξ± - 1) + (1 - Ξ΄)) +# c[0] + k[0] = (1 - Ξ΄) * k[-1] + exp(z[0]) * k[-1]^Ξ± +# z[0] = 0.2 * z[-1] + 0.01 * eps_z[x] +# end - data = simulate(FS2000)([:c,:k],:,:simulate) - get_loglikelihood(FS2000, data, FS2000.parameter_values) - get_mean(FS2000, silent = true) - # get_SSS(FS2000, silent = true) - # get_SSS(FS2000, algorithm = :third_order, silent = true) +# @parameters RBC silent = true precompile = true begin +# Ξ΄ = 0.02 +# Ξ± = 0.5 +# end - # import StatsPlots - # plot_irf(FS2000) - # plot_solution(FS2000,:k) # fix warning when there is no sensitivity and all values are the same. triggers: no strict ticks found... - # plot_conditional_variance_decomposition(FS2000) - end - end -end +# get_SS(FS2000) +# get_SS(FS2000, parameters = :alp => 0.36) +# get_solution(FS2000) +# get_solution(FS2000, parameters = :alp => 0.35) +# get_standard_deviation(FS2000) +# get_correlation(FS2000) +# get_autocorrelation(FS2000) +# get_variance_decomposition(FS2000) +# get_conditional_variance_decomposition(FS2000) +# get_irf(FS2000) + +# data = simulate(FS2000)([:c,:k],:,:simulate) +# get_loglikelihood(FS2000, data, FS2000.parameter_values) +# get_mean(FS2000, silent = true) +# # get_SSS(FS2000, silent = true) +# # get_SSS(FS2000, algorithm = :third_order, silent = true) + +# # import StatsPlots +# # plot_irf(FS2000) +# # plot_solution(FS2000,:k) # fix warning when there is no sensitivity and all values are the same. triggers: no strict ticks found... +# # plot_conditional_variance_decomposition(FS2000) +# end +# end +# end end diff --git a/src/get_functions.jl b/src/get_functions.jl index 7ca0af3b8..1aec1d473 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -962,9 +962,9 @@ And data, 4Γ—40Γ—1 Array{Float64, 3}: function get_irf(𝓂::β„³; periods::Int = 40, algorithm::Symbol = :first_order, - parameters::ParameterType = nothing, + parameters::Union{ParameterType, KeyedArray{Float64}, Dict{Int, Dict{Symbol, Real}}} = nothing, variables::Union{Symbol_input,String_input} = :all_excluding_obc, - shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all_excluding_obc, + shocks::Union{Symbol_input, String_input, Matrix{Float64}, KeyedArray{Float64}} = :all_excluding_obc, negative_shock::Bool = false, generalised_irf::Bool = false, initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = [0.0], @@ -984,46 +984,100 @@ function get_irf(𝓂::β„³; obc_model = length(𝓂.obc_violation_equations) > 0 + if ignore_obc + occasionally_binding_constraints = false + else + occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 + end + + if parameters isa ParameterType + break_points_dict = Dict{Int, Dict{Symbol, Float64}}(0 => Dict{Symbol, Float64}()) + periods_of_change = [0] + else + break_points_dict = transform_break_points(parameters) + periods_of_change = break_points_dict |> keys |> collect |> sort + end + if shocks isa Matrix{Float64} @assert size(shocks)[1] == 𝓂.timings.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." - + periods += size(shocks)[2] - - shock_history = zeros(𝓂.timings.nExo, periods) - - shock_history[:,1:size(shocks)[2]] = shocks - - shock_idx = 1 - + + shock_idx = [1] + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[:, 1:size(shocks)[2], 1] = shocks + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:]) > 1e-10 elseif shocks isa KeyedArray{Float64} shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) - + periods += size(shocks)[2] - + @assert length(setdiff(shock_input, 𝓂.timings.exo)) == 0 "Provided shocks which are not part of the model." - - shock_history = zeros(𝓂.timings.nExo, periods + 1) - - shock_history[indexin(shock_input,𝓂.timings.exo),1:size(shocks)[2]] = shocks - - shock_idx = 1 - + + shock_idx = [1] + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[indexin(shock_input,𝓂.timings.exo), 1:size(shocks)[2], 1] = shocks + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks(intersect(𝓂.timings.exo,axiskeys(shocks,1)),:)) > 1e-10 - else + elseif shocks == :simulate + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:] .= 0 + + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) + elseif shocks == :none + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) + elseif shocks isa Union{Symbol_input,String_input} shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + for (i,ii) in enumerate(shock_idx) + shock_history[ii,1,i] = negative_shock ? -1 : 1 + end obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) end + + solve!(𝓂, parameters = parameters isa ParameterType ? parameters : nothing, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) - if ignore_obc - occasionally_binding_constraints = false + var_idx = parse_variables_input_to_index(variables, 𝓂.timings) + + axis1 = 𝓂.timings.var[var_idx] + + if any(x -> contains(string(x), "β—–"), axis1) + axis1_decomposed = decompose_name.(axis1) + axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] + end + + if shocks == :simulate + axis2 = [:simulate] + elseif shocks == :none + axis2 = [:none] else - occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 + axis2 = shocks isa Union{Symbol_input,String_input} ? + shock_idx isa Int ? + [𝓂.timings.exo[shock_idx]] : + 𝓂.timings.exo[shock_idx] : + [:Shock_matrix] + + if any(x -> contains(string(x), "β—–"), axis2) + axis2_decomposed = decompose_name.(axis2) + axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] + end end - solve!(𝓂, parameters = parameters, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) - reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) unspecified_initial_state = initial_state == [0.0] @@ -1052,16 +1106,30 @@ function get_irf(𝓂::β„³; end end - if occasionally_binding_constraints - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) - elseif obc_shocks_included - @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." + steady_states_and_state_update = Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}() - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) - else - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + pre_break_point_parameters = deepcopy(𝓂.parameter_values) + + for p in periods_of_change + solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + + if occasionally_binding_constraints + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + elseif obc_shocks_included + @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." + + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + else + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + end + + reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + + steady_states_and_state_update[p] = (reference_steady_state, NSSS, SSS_delta, state_update) end + solve!(𝓂, parameters = pre_break_point_parameters, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + if generalised_irf girfs = girf(state_update, initial_state, @@ -1150,14 +1218,14 @@ function get_irf(𝓂::β„³; variables = variables, negative_shock = negative_shock) else - irfs = irf(state_update, + raw_irfs = irf(steady_states_and_state_update, initial_state, - levels ? reference_steady_state + SSS_delta : SSS_delta, - 𝓂.timings; - periods = periods, - shocks = shocks, - variables = variables, - negative_shock = negative_shock) + shock_history, + shock_idx, + levels, + 𝓂.timings) + + irfs = KeyedArray(raw_irfs[var_idx,:,:]; Variables = axis1, Periods = 1:(periods + periods_of_change[end]), Shocks = axis2) end return irfs diff --git a/test/test_break_points.jl b/test/test_break_points.jl index 1eb9b7f61..7115cc217 100644 --- a/test/test_break_points.jl +++ b/test/test_break_points.jl @@ -1,26 +1,70 @@ + using MacroModelling +function transform_break_points(break_points) + dict = Dict{Int, Dict{Symbol, Float64}}() + + for i in 1:size(break_points,2) + time = axiskeys(break_points,2)[i] + + if !haskey(dict, time) + dict[time] = Dict{Symbol, Float64}() + end + + for k in 1:size(break_points,1) + parameter = axiskeys(break_points,1)[k] + value = break_points[k,i] + if !(isnothing(value) || !isfinite(value)) + dict[time][parameter] = value + end + end + end + + if !haskey(dict, 0) + dict[0] = Dict{Symbol, Float64}() + end + + return dict +end + include("../models/Smets_Wouters_2007.jl") irfs = get_irf(Smets_Wouters_2007) -break_points = KeyedArray([.4,.3,.25]', Variable =[:ctrend], Time = [1,15,30]) +irfs = get_irf(Smets_Wouters_2007, shocks = :none, levels = true) +get_parameters(Smets_Wouters_2007, values = true) + +break_points = KeyedArray([.4,.3,.25]', Variable =[:ctrend], Time = [3,15,30]) + +irfs = get_irf(Smets_Wouters_2007, parameters = break_points, shocks = :none, levels = true) + +starting_vals = get_irf(Smets_Wouters_2007, shocks = :none, levels = true, periods = 1) +irfalt = get_irf(Smets_Wouters_2007, parameters = :ctrend => .4, shocks = :none, levels = true, initial_state = vec(starting_vals)) -# translate this to the combined parameter vector -SS() -import MacroModelling: get_relevant_steady_states, parse_shocks_input_to_index, parse_algorithm_to_state_update, obc_objective_optim_fun, obc_constraint_optim_fun, String_input, Symbol_input, timings + +parameters = KeyedArray( + [nothing .3 .25 + .77 NaN .9], + Variable = [:ctrend, :constepinf], + Time = [2, 15, 30]) + + + +# translate this to the combined parameter vector +import MacroModelling: get_relevant_steady_states, parse_shocks_input_to_index, parse_algorithm_to_state_update, obc_objective_optim_fun, obc_constraint_optim_fun, String_input, Symbol_input, timings, parse_variables_input_to_index, ParameterType 𝓂 = Smets_Wouters_2007 -periods = 40 +T = 𝓂.timings +periods = 100 algorithm = :first_order -parameters = nothing +# parameters = nothing variables = :all_excluding_obc -shocks = :all_excluding_obc +shocks = :none # :all_excluding_obc negative_shock = false generalised_irf = false initial_state = [0.0] @@ -28,6 +72,10 @@ levels = false ignore_obc = false verbose = false + + + + shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks @@ -40,16 +88,30 @@ stochastic_model = length(𝓂.timings.exo) > 0 obc_model = length(𝓂.obc_violation_equations) > 0 +if ignore_obc + occasionally_binding_constraints = false +else + occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 +end + +if parameters isa ParameterType + break_points_dict = Dict{Int, Dict{Symbol, Float64}}(0 => Dict{Symbol, Float64}()) + periods_of_change = [0] +else + break_points_dict = transform_break_points(parameters) + periods_of_change = break_points_dict |> keys |> collect |> sort +end + if shocks isa Matrix{Float64} @assert size(shocks)[1] == 𝓂.timings.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." periods += size(shocks)[2] - shock_history = zeros(𝓂.timings.nExo, periods) - - shock_history[:,1:size(shocks)[2]] = shocks + shock_idx = [1] + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) - shock_idx = 1 + shock_history[:, 1:size(shocks)[2], 1] = shocks obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:]) > 1e-10 elseif shocks isa KeyedArray{Float64} @@ -59,41 +121,44 @@ elseif shocks isa KeyedArray{Float64} @assert length(setdiff(shock_input, 𝓂.timings.exo)) == 0 "Provided shocks which are not part of the model." - shock_history = zeros(𝓂.timings.nExo, periods + 1) + shock_idx = [1] - shock_history[indexin(shock_input,𝓂.timings.exo),1:size(shocks)[2]] = shocks + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) - shock_idx = 1 + shock_history[indexin(shock_input,𝓂.timings.exo), 1:size(shocks)[2], 1] = shocks obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks(intersect(𝓂.timings.exo,axiskeys(shocks,1)),:)) > 1e-10 -else +elseif shocks == :simulate shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:] .= 0 + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) -end +elseif shocks == :none + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) -if ignore_obc - occasionally_binding_constraints = false -else - occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 -end + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) +elseif shocks isa Union{Symbol_input,String_input} + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) -irfs = irf(state_update, - initial_state, - levels ? reference_steady_state + SSS_delta : SSS_delta, - 𝓂.timings; - periods = periods, - shocks = shocks, - variables = variables, - negative_shock = negative_shock) + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + for (i,ii) in enumerate(shock_idx) + shock_history[ii,1,i] = negative_shock ? -1 : 1 + end + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) +end +solve!(𝓂, parameters = parameters isa ParameterType ? parameters : nothing, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) -var_idx = parse_variables_input_to_index(variables, T) +var_idx = parse_variables_input_to_index(variables, 𝓂.timings) -axis1 = T.var[var_idx] +axis1 = 𝓂.timings.var[var_idx] if any(x -> contains(string(x), "β—–"), axis1) axis1_decomposed = decompose_name.(axis1) @@ -101,8 +166,6 @@ if any(x -> contains(string(x), "β—–"), axis1) end - - if shocks == :simulate axis2 = [:simulate] elseif shocks == :none @@ -110,8 +173,8 @@ elseif shocks == :none else axis2 = shocks isa Union{Symbol_input,String_input} ? shock_idx isa Int ? - [T.exo[shock_idx]] : - T.exo[shock_idx] : + [𝓂.timings.exo[shock_idx]] : + 𝓂.timings.exo[shock_idx] : [:Shock_matrix] if any(x -> contains(string(x), "β—–"), axis2) @@ -120,154 +183,284 @@ else end end +reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) -function irff(state_update::Function, - initial_state::Union{Vector{Vector{Float64}},Vector{Float64},Matrix{Float64}}, - level::Vector{Float64}, - T::timings; - periods::Int = 40, - shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, - variables::Union{Symbol_input,String_input} = :all, - negative_shock::Bool = false) - - pruning = initial_state isa Vector{Vector{Float64}} +unspecified_initial_state = initial_state == [0.0] - shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks +if unspecified_initial_state + if algorithm == :pruned_second_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + else + initial_state = zeros(𝓂.timings.nVars) - SSS_delta + end +else + if initial_state isa Vector{Float64} + if algorithm == :pruned_second_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + else + initial_state = initial_state - NSSS + end + else + if algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] + @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." + end + end +end - shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks - if shocks isa Matrix{Float64} - @assert size(shocks)[1] == T.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." +steady_states_and_state_update = Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}() - # periods += size(shocks)[2] +for p in periods_of_change + solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) - shock_history = zeros(T.nExo, periods) + if occasionally_binding_constraints + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + elseif obc_shocks_included + @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." - shock_history[:,1:size(shocks)[2]] = shocks + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + else + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + end - shock_idx = 1 - elseif shocks isa KeyedArray{Float64} - shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) + reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) - # periods += size(shocks)[2] + steady_states_and_state_update[p] = (reference_steady_state, NSSS, SSS_delta, state_update) +end - @assert length(setdiff(shock_input, T.exo)) == 0 "Provided shocks which are not part of the model." - - shock_history = zeros(T.nExo, periods) - shock_history[indexin(shock_input,T.exo),1:size(shocks)[2]] = shocks +periods += periods_of_change[end] +# function irf(steady_states_and_state_update::Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}, +# initial_state::Union{Vector{Vector{Float64}},Vector{Float64}}, +# shock_history::Matrix{Float64}, +# T::timings; +# periods::Int = 40, +# negative_shock::Bool = false) - shock_idx = 1 - else - shock_idx = parse_shocks_input_to_index(shocks,T) - end +Y = zeros(T.nVars, periods, length(shock_idx)) - if shocks == :simulate - shock_history = randn(T.nExo,periods) +pruning = initial_state isa Vector{Vector{Float64}} - shock_history[contains.(string.(T.exo),"α΅’α΅‡αΆœ"),:] .= 0 +periods_of_change = steady_states_and_state_update |> keys |> collect |> sort - Y = zeros(T.nVars,periods,1) +initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] - initial_state = state_update(initial_state,shock_history[:,1]) +for i in shock_idx + for (k,p) in enumerate(periods_of_change) + concerned_periods = (k == 1 ? 1 : periods_of_change[k]):(k == length(periods_of_change) ? periods : periods_of_change[k+1] - 1) - Y[:,1,1] = pruning ? sum(initial_state) : initial_state + println(concerned_periods) + reference_steady_state, NSSS, SSS_delta, state_update = steady_states_and_state_update[p] - for t in 1:periods-1 - initial_state = state_update(initial_state,shock_history[:,t+1]) + if k > 1 + Ξ”SS = steady_states_and_state_update[periods_of_change[k]][2] - steady_states_and_state_update[periods_of_change[k-1]][2] - Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state + if pruning + for j in initial_state_copy[i] + j += Ξ”SS + end + else + initial_state_copy[i] += Ξ”SS + end end - elseif shocks == :none - Y = zeros(T.nVars,periods,1) - shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : zeros(T.nExo) + for t in concerned_periods + initial_state_copy[i] = state_update(initial_state_copy[i], shock_history[:,t,i]) - initial_state = state_update(initial_state, shck) + Y[:,t,i] = pruning ? sum(initial_state_copy[i]) : initial_state_copy[i] + end + end +end - Y[:,1,1] = pruning ? sum(initial_state) : initial_state - for t in 1:periods-1 - initial_state = state_update(initial_state, shck) +Y[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] - Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state - end - else - Y = zeros(T.nVars,periods,length(shock_idx)) +(Y .+ reference_steady_state)[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] - for (i,ii) in enumerate(shock_idx) - initial_state_copy = deepcopy(initial_state) - - if shocks != :simulate && shocks isa Union{Symbol_input,String_input} - shock_history = zeros(T.nExo,periods) - shock_history[ii,1] = negative_shock ? -1 : 1 - end +reference_steady_states[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:] - initial_state_copy = state_update(initial_state_copy, shock_history[:,1]) - Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy - for t in 1:periods-1 - initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1]) +Y = zeros(T.nVars,periods,length(shock_idx)) - Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy - end - end +for (i,ii) in enumerate(shock_idx) + initial_state_copy = deepcopy(initial_state) + + if shocks != :simulate && shocks isa Union{Symbol_input,String_input} + shock_history = zeros(T.nExo,periods) + shock_history[ii,1] = negative_shock ? -1 : 1 end - return Y + initial_state_copy = state_update(initial_state_copy, shock_history[:,1]) + + Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + + for t in 1:periods-1 + initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1]) + + Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + end end -break_points -break_points = KeyedArray([nothing .3 .25 -.77 NaN .9], Variable =[:ctrend, :constepinf], Time = [1,15,30]) -function transform_break_points(break_points) - dict = Dict{Int, Dict{Symbol, Float64}}() + + + +get_parameters(Smets_Wouters_2007, values = true) + +# solve!(𝓂, parameters = break_points_dict[30], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) - for i in 1:size(break_points,2) - time = axiskeys(break_points,2)[i] +# reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) - if !haskey(dict, time) - dict[time] = Dict{Symbol, Float64}() + +initial_state = [0.0] + +pruning = initial_state isa Vector{Vector{Float64}} + +Y = zeros(T.nVars,periods + periods_of_change[end],length(shock_idx)); + +reference_steady_states = zeros(T.nVars,periods + periods_of_change[end]); + +for (i,p) in enumerate(periods_of_change) +# i = 1 +# p = periods_of_change[i] + + concerned_periods = (i == 1 ? 1 : periods_of_change[i]):(i == length(periods_of_change) ? p + periods : periods_of_change[i+1] - 1) + # periods_with_these_parameters = (i == length(periods_of_change) ? p + periods : periods_of_change[i+1]) - (i == 1 ? 1 : periods_of_change[i]) + + # if periods_with_these_parameters == 0 continue end + + solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + + reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + + reference_steady_states[:,concerned_periods] .= reference_steady_state + + Ξ”reference_steady_state = reference_steady_state - reference_steady_stateꜜ + Ξ”NSSS = NSSS - NSSSꜜ + Ξ”SSS_delta = SSS_delta - SSS_deltaꜜ + + # println(maximum(abs, Ξ”reference_steady_state), maximum(abs, Ξ”NSSS), maximum(abs, Ξ”SSS_delta)) + + unspecified_initial_state = initial_state == [0.0] + + if unspecified_initial_state + if algorithm == :pruned_second_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] + Ξ”initial_state = [Ξ”reference_steady_state[1:𝓂.timings.nVars], Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + Ξ”initial_state = [Ξ”reference_steady_state[1:𝓂.timings.nVars], Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta, Ξ”reference_steady_state[1:𝓂.timings.nVars]] + else + initial_state = zeros(𝓂.timings.nVars) - SSS_delta + Ξ”initial_state = Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta end - - for k in 1:size(break_points,1) - parameter = axiskeys(break_points,1)[k] - value = break_points[k,i] - if !(isnothing(value) || !isfinite(value)) - dict[time][parameter] = value + else + if initial_state isa Vector{Float64} + if algorithm == :pruned_second_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] + Ξ”initial_state = [Ξ”reference_steady_state[1:𝓂.timings.nVars], Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta] + elseif algorithm == :pruned_third_order + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] + Ξ”initial_state = [Ξ”reference_steady_state[1:𝓂.timings.nVars], Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta, Ξ”reference_steady_state[1:𝓂.timings.nVars]] + else + initial_state = initial_state - NSSS + Ξ”initial_state = Ξ”reference_steady_state[1:𝓂.timings.nVars] - Ξ”SSS_delta + end + else + if algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] + @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." end end end + + if occasionally_binding_constraints + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + elseif obc_shocks_included + @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." + + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + else + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + end - # if !haskey(dict, 0) - # dict[0] = Dict{Symbol, Float64}() - # end + if 1 ∈ concerned_periods + initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] + end - return dict + initial_state_copy = [s - Ξ”initial_state for s in initial_state_copy] + + for (i,ii) in enumerate(shock_idx) + for t in concerned_periods .- 1 + initial_state_copy[i] = state_update(initial_state_copy[i], shock_history[:,t+1,i]) + + Y[:,t+1,i] = pruning ? sum(initial_state_copy[i]) : initial_state_copy[i] + end + end + reference_steady_stateꜜ, NSSSꜜ, SSS_deltaꜜ = reference_steady_state, NSSS, SSS_delta + # println("concerned_periods: ", concerned_periods) + # println("periods_with_these_parameters: ", periods_with_these_parameters) + # println("parameters: ", break_points_dict[p]) end -break_points_dict = transform_break_points(break_points) -break_points_dict -previous_parameter_values = deepcopy(𝓂.parameter_values) +(Y .+ reference_steady_states)[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] -periods_of_change = break_points_dict |> keys |> collect |> sort +reference_steady_states[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:] -solve!(𝓂, parameters = break_points_dict[0], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) -for (i,p) in enumerate(periods_of_change) - concerned_periods = (i == 1 ? 1 : periods_of_change[i-1]):periods_of_change[i] - periods_with_these_parameters = periods_of_change[i] - (i == 1 ? 1 : periods_of_change[i-1]) - solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) +# check against alternative way + +solve!(𝓂, parameters = previous_parameter_values, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + +reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + +iirrff = get_irf(𝓂, initial_state = reference_steady_state[1:𝓂.timings.nVars], parameters = break_points_dict[2], shocks = :none) + +iirrff([:dinve,:pinfobs,:robs,:dwobs,:labobs],:,:) + +Y[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] + + + +solve!(𝓂, parameters = previous_parameter_values, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + + + + +solve!(𝓂, parameters = previous_parameter_values, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + +solve!(𝓂, parameters = break_points_dict[1], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + +i = 2 + +concerned_periods = (i == 1 ? 1 : periods_of_change[i-1]):periods_of_change[i] +periods_with_these_parameters = periods_of_change[i] - (i == 1 ? 1 : periods_of_change[i-1]) + +# for (i,p) in enumerate(periods_of_change) +# concerned_periods = (i == 1 ? 1 : periods_of_change[i-1]):periods_of_change[i] +# periods_with_these_parameters = periods_of_change[i+1] - (i == 1 ? 1 : periods_of_change[i]) +# if p == 1 +# pars = break_points_dict[p] +# else +# pars = nothing +# end + +# write a function that goes from breakpoint to breakpoint, solves the model with the new parameters, and then computes the IRFs for as many periods until the next breakpoint + + + solve!(𝓂, parameters = pars, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) @@ -316,7 +509,7 @@ for (i,p) in enumerate(periods_of_change) variables = variables, negative_shock = negative_shock) -end +# end irfs1 = irff(state_update, initial_state, @@ -327,4 +520,113 @@ irfs1 = irff(state_update, variables = variables, negative_shock = negative_shock) -irfs1[:,1,:] \ No newline at end of file +irfs1[:,1,:] + + + + + + + + + + + +function irff(state_update::Function, + initial_state::Union{Vector{Vector{Float64}},Vector{Float64},Matrix{Float64}}, + level::Vector{Float64}, + T::timings; + periods::Int = 40, + shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, + variables::Union{Symbol_input,String_input} = :all, + negative_shock::Bool = false) + + pruning = initial_state isa Vector{Vector{Float64}} + + shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks + + shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks + + if shocks isa Matrix{Float64} + @assert size(shocks)[1] == T.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." + + # periods += size(shocks)[2] + + shock_history = zeros(T.nExo, periods) + + shock_history[:,1:size(shocks)[2]] = shocks + + shock_idx = [1] + elseif shocks isa KeyedArray{Float64} + shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) + + # periods += size(shocks)[2] + + @assert length(setdiff(shock_input, T.exo)) == 0 "Provided shocks which are not part of the model." + + shock_history = zeros(T.nExo, periods) + + shock_history[indexin(shock_input,T.exo),1:size(shocks)[2]] = shocks + + shock_idx = [1] + else + shock_idx = parse_shocks_input_to_index(shocks,T) + end + + + if shocks == :simulate + shock_history = randn(T.nExo,periods) + + shock_history[contains.(string.(T.exo),"α΅’α΅‡αΆœ"),:] .= 0 + + Y = zeros(T.nVars,periods,1) + + initial_state = state_update(initial_state,shock_history[:,1]) + + Y[:,1,1] = pruning ? sum(initial_state) : initial_state + + for t in 1:periods-1 + initial_state = state_update(initial_state,shock_history[:,t+1]) + + Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state + end + elseif shocks == :none + Y = zeros(T.nVars,periods,1) + + shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : zeros(T.nExo) + + initial_state = state_update(initial_state, shck) + + Y[:,1,1] = pruning ? sum(initial_state) : initial_state + + for t in 1:periods-1 + initial_state = state_update(initial_state, shck) + + Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state + end + else + Y = zeros(T.nVars,periods,length(shock_idx)) + + for (i,ii) in enumerate(shock_idx) + initial_state_copy = deepcopy(initial_state) + + if shocks != :simulate && shocks isa Union{Symbol_input,String_input} + shock_history = zeros(T.nExo,periods) + shock_history[ii,1] = negative_shock ? -1 : 1 + end + + initial_state_copy = state_update(initial_state_copy, shock_history[:,1]) + + Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + + for t in 1:periods-1 + initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1]) + + Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy + end + end + end + + return Y +end + From 584d856ab8ddb10b0bb88279f5060deb7a77ce4c Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 24 Apr 2024 22:43:33 +0200 Subject: [PATCH 4/7] fixed visualisation issues with plots --- src/MacroModelling.jl | 15 ++-- src/plotting.jl | 200 ++++++++++++++++++++++++++++++------------ 2 files changed, 154 insertions(+), 61 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 24ccf1ca7..e55c34750 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -6352,21 +6352,24 @@ function irf(steady_states_and_state_update::Dict{Int, Tuple{Vector{Float64}, Ve reference_steady_state, __, SSS_delta, state_update = steady_states_and_state_update[p] if k > 1 - Ξ”SS = steady_states_and_state_update[periods_of_change[k]][2] - steady_states_and_state_update[periods_of_change[k-1]][2] - + Ξ”NSSS = steady_states_and_state_update[periods_of_change[k-1]][2] - steady_states_and_state_update[periods_of_change[k]][2] + Ξ”SSS_delta = steady_states_and_state_update[periods_of_change[k-1]][3] - steady_states_and_state_update[periods_of_change[k]][3] + if pruning - for j in initial_state_copy[i] - j += Ξ”SS + initial_state_copy[i][1] += Ξ”NSSS + initial_state_copy[i][2] -= Ξ”SSS_delta + if length(initial_state_copy[i]) > 3 + initial_state_copy[i][3] += Ξ”NSSS end else - initial_state_copy[i] += Ξ”SS + initial_state_copy[i] += Ξ”NSSS end end for t in concerned_periods initial_state_copy[i] = state_update(initial_state_copy[i], shock_history[:,t,i]) - Y[:,t,i] = pruning ? sum(initial_state_copy[i]) : initial_state_copy[i] .+ (levels ? reference_steady_state + SSS_delta : SSS_delta) + Y[:,t,i] = (pruning ? sum(initial_state_copy[i]) : initial_state_copy[i]) .+ (levels ? reference_steady_state + SSS_delta : SSS_delta) end end end diff --git a/src/plotting.jl b/src/plotting.jl index 63175cd05..3a8c9095d 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -104,7 +104,7 @@ function plot_model_estimates(𝓂::β„³, guidefont = 8, legendfontsize = 8, tickfontsize = 8, - framestyle = :box) + framestyle = :semi) # write_parameters_input!(𝓂, parameters, verbose = verbose) @@ -378,7 +378,7 @@ function plot_irf(𝓂::β„³; periods::Int = 40, shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all_excluding_obc, variables::Union{Symbol_input,String_input} = :all_excluding_auxilliary_and_obc, - parameters::ParameterType = nothing, + parameters::Union{ParameterType, KeyedArray{Float64}, Dict{Int, Dict{Symbol, Real}}} = nothing, show_plots::Bool = true, save_plots::Bool = false, save_plots_format::Symbol = :pdf, @@ -399,51 +399,113 @@ function plot_irf(𝓂::β„³; guidefont = 8, legendfontsize = 8, tickfontsize = 8, - framestyle = :box) + framestyle = :semi) shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks shocks = shocks isa String_input ? shocks .|> Meta.parse .|> replace_indices : shocks - + shocks = 𝓂.timings.nExo == 0 ? :none : shocks + @assert !(shocks == :none && generalised_irf) "Cannot compute generalised IRFs for model without shocks." + stochastic_model = length(𝓂.timings.exo) > 0 obc_model = length(𝓂.obc_violation_equations) > 0 - if shocks isa Matrix{Float64} - @assert size(shocks)[1] == 𝓂.timings.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." + if ignore_obc + occasionally_binding_constraints = false + else + occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 + end - shock_idx = 1 + if parameters isa ParameterType + break_points_dict = Dict{Int, Dict{Symbol, Float64}}(0 => Dict{Symbol, Float64}()) + periods_of_change = [0] + else + break_points_dict = transform_break_points(parameters) + periods_of_change = break_points_dict |> keys |> collect |> sort + end + if shocks isa Matrix{Float64} + @assert size(shocks)[1] == 𝓂.timings.nExo "Number of rows of provided shock matrix does not correspond to number of shocks. Please provide matrix with as many rows as there are shocks in the model." + + periods += size(shocks)[2] + + shock_idx = [1] + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[:, 1:size(shocks)[2], 1] = shocks + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:]) > 1e-10 elseif shocks isa KeyedArray{Float64} - shock_idx = 1 - - obc_shocks = 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")] - - obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks(intersect(obc_shocks, axiskeys(shocks,1)),:)) > 1e-10 - else + shock_input = map(x->Symbol(replace(string(x),"β‚β‚“β‚Ž" => "")),axiskeys(shocks)[1]) + + periods += size(shocks)[2] + + @assert length(setdiff(shock_input, 𝓂.timings.exo)) == 0 "Provided shocks which are not part of the model." + + shock_idx = [1] + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[indexin(shock_input,𝓂.timings.exo), 1:size(shocks)[2], 1] = shocks + + obc_shocks_included = stochastic_model && obc_model && sum(abs2,shocks(intersect(𝓂.timings.exo,axiskeys(shocks,1)),:)) > 1e-10 + elseif shocks == :simulate shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) - + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + shock_history[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ"),:] .= 0 + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) - end + elseif shocks == :none + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) + elseif shocks isa Union{Symbol_input,String_input} + shock_idx = parse_shocks_input_to_index(shocks,𝓂.timings) + + shock_history = zeros(𝓂.timings.nExo, periods + periods_of_change[end], length(shock_idx)) + + for (i,ii) in enumerate(shock_idx) + shock_history[ii,1,i] = negative_shock ? -1 : 1 + end - if shocks isa KeyedArray{Float64} || shocks isa Matrix{Float64} - periods = max(periods, size(shocks)[2]) + obc_shocks_included = stochastic_model && obc_model && (intersect((((shock_idx isa Vector) || (shock_idx isa UnitRange)) && (length(shock_idx) > 0)) ? 𝓂.timings.exo[shock_idx] : [𝓂.timings.exo[shock_idx]], 𝓂.timings.exo[contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")]) != []) end - - variables = variables isa String_input ? variables .|> Meta.parse .|> replace_indices : variables + + solve!(𝓂, parameters = parameters isa ParameterType ? parameters : nothing, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) var_idx = parse_variables_input_to_index(variables, 𝓂.timings) - if ignore_obc - occasionally_binding_constraints = false - else - occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 + axis1 = 𝓂.timings.var[var_idx] + + if any(x -> contains(string(x), "β—–"), axis1) + axis1_decomposed = decompose_name.(axis1) + axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] end - solve!(𝓂, parameters = parameters, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + if shocks == :simulate + axis2 = [:simulate] + elseif shocks == :none + axis2 = [:none] + else + axis2 = shocks isa Union{Symbol_input,String_input} ? + shock_idx isa Int ? + [𝓂.timings.exo[shock_idx]] : + 𝓂.timings.exo[shock_idx] : + [:Shock_matrix] + + if any(x -> contains(string(x), "β—–"), axis2) + axis2_decomposed = decompose_name.(axis2) + axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] + end + end reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) @@ -464,24 +526,47 @@ function plot_irf(𝓂::β„³; elseif algorithm == :pruned_third_order initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] else - initial_state = initial_state - reference_steady_state[1:𝓂.timings.nVars] + initial_state = initial_state - NSSS end else - @assert algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] && initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." + if algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] + @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." + end end end + + steady_states_and_state_update = Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}() + + pre_break_point_parameters = deepcopy(𝓂.parameter_values) + reference_steady_states = zeros(𝓂.timings.nVars ,periods + periods_of_change[end]) - if occasionally_binding_constraints - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) - elseif obc_shocks_included - @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks without enforcing the constraint is only compatible with first order perturbation solutions." + for (k,p) in enumerate(periods_of_change) + solve!(𝓂, parameters = break_points_dict[p], verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) - else - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) - end + if occasionally_binding_constraints + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + elseif obc_shocks_included + @assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraint shocks witout enforcing the constraint is only compatible with first order perturbation solutions." + + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) + else + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + end + + reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + + steady_states_and_state_update[p] = (reference_steady_state, NSSS, SSS_delta, state_update) + concerned_periods = (k == 1 ? 1 : periods_of_change[k]):(k == length(periods_of_change) ? periods + periods_of_change[end] : periods_of_change[k+1] - 1) + + for t in concerned_periods + reference_steady_states[:,t] .= reference_steady_state + end + end + + solve!(𝓂, parameters = pre_break_point_parameters, verbose = verbose, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) + if generalised_irf Y = girf(state_update, initial_state, @@ -611,14 +696,16 @@ function plot_irf(𝓂::β„³; variables = variables, negative_shock = negative_shock) .+ SSS_delta[var_idx] else - Y = irf(state_update, - initial_state, - zeros(𝓂.timings.nVars), - 𝓂.timings; - periods = periods, - shocks = shocks, - variables = variables, - negative_shock = negative_shock) .+ SSS_delta[var_idx] + levels = false + + raw_Y = irf(steady_states_and_state_update, + initial_state, + shock_history, + shock_idx, + levels, + 𝓂.timings) + + Y = KeyedArray(raw_Y[var_idx,:,:]; Variables = axis1, Periods = 1:(periods + periods_of_change[end]), Shocks = axis2) end end @@ -648,27 +735,30 @@ function plot_irf(𝓂::β„³; end for i in 1:length(var_idx) - SS = reference_steady_state[var_idx[i]] + SS = reference_steady_states[var_idx[i],:] - can_dual_axis = gr_back && all((Y[i,:,shock] .+ SS) .> eps(Float32)) && (SS > eps(Float32)) + can_dual_axis = gr_back && all((Y[i,:,shock] .+ SS) .> eps(Float32)) && SS[1] .> eps(Float32) # && periods_of_change == 1 if !(all(isapprox.(Y[i,:,shock],0,atol = eps(Float32)))) push!(pp,begin - StatsPlots.plot(Y[i,:,shock] .+ SS, - title = replace_indices_in_symbol(𝓂.timings.var[var_idx[i]]), + StatsPlots.plot(SS, ylabel = "Level", - label = "") + label = "", + title = replace_indices_in_symbol(𝓂.timings.var[var_idx[i]]), + color = :black) if can_dual_axis StatsPlots.plot!(StatsPlots.twinx(), - 100*((Y[i,:,shock] .+ SS) ./ SS .- 1), + 100*(SS .- SS[1]) ./ SS[1], ylabel = LaTeXStrings.L"\% \Delta", - label = "") + label = "", + color = :black) end - StatsPlots.hline!(can_dual_axis ? [SS 0] : [SS], - color = :black, - label = "") + StatsPlots.plot!(can_dual_axis ? [collect(Y[i,:,shock] .+ SS) collect(100*((Y[i,:,shock] .+ SS) ./ SS[1] .- 1))] : [Y[i,:,shock] .+ SS], + label = "", + grid = :all, + color = 1) end) if !(plot_count % plots_per_page == 0) @@ -850,7 +940,7 @@ function plot_conditional_variance_decomposition(𝓂::β„³; guidefont = 8, legendfontsize = 8, tickfontsize = 8, - framestyle = :box) + framestyle = :semi) fevds = get_conditional_variance_decomposition(𝓂, periods = 1:periods, @@ -1028,7 +1118,7 @@ function plot_solution(𝓂::β„³, guidefont = 8, legendfontsize = 8, tickfontsize = 8, - framestyle = :box) + framestyle = :semi) state = state isa Symbol ? state : state |> Meta.parse |> replace_indices @@ -1378,7 +1468,7 @@ function plot_conditional_forecast(𝓂::β„³, guidefont = 8, legendfontsize = 8, tickfontsize = 8, - framestyle = :box) + framestyle = :semi) conditions = conditions isa KeyedArray ? axiskeys(conditions,1) isa Vector{String} ? rekey(conditions, 1 => axiskeys(conditions,1) .|> Meta.parse .|> replace_indices) : conditions : conditions From a276e4d1651ef0cf27cf361ac1956595fcc4e0da Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 24 Apr 2024 23:31:50 +0200 Subject: [PATCH 5/7] removed pruned higher order effect on stochastic vol states --- src/MacroModelling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index e55c34750..9435753c2 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -6357,7 +6357,7 @@ function irf(steady_states_and_state_update::Dict{Int, Tuple{Vector{Float64}, Ve if pruning initial_state_copy[i][1] += Ξ”NSSS - initial_state_copy[i][2] -= Ξ”SSS_delta + # initial_state_copy[i][2] += Ξ”NSSS#Ξ”SSS_delta if length(initial_state_copy[i]) > 3 initial_state_copy[i][3] += Ξ”NSSS end From 0c7ae32c787d27f53375989d5a9c3aa6b6584461 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Thu, 25 Apr 2024 00:25:53 +0200 Subject: [PATCH 6/7] Fix typo in conditions_in_levels variable name --- src/plotting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plotting.jl b/src/plotting.jl index 3a8c9095d..eca275e64 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -1378,7 +1378,7 @@ The left axis shows the level, and the right axis the deviation from the (non) s - `periods` [Default: `40`, Type: `Int`]: the total number of periods is the sum of the argument provided here and the maximum of periods of the shocks or conditions argument. - $PARAMETERS - $VARIABLES -`conditions_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the conditions are provided in levels. If `true` the input to the conditions argument will have the non stochastic steady state substracted. +- `conditions_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the conditions are provided in levels. If `true` the input to the conditions argument will have the non stochastic steady state substracted. - $LEVELS - $ALGORITHM - `show_plots` [Default: `true`, Type: `Bool`]: show plots. Separate plots per shocks and varibles depending on number of variables and `plots_per_page`. From 153fd747394bd5de284419ed08e97e2dedfc39e5 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Thu, 25 Apr 2024 00:40:53 +0200 Subject: [PATCH 7/7] Add StatsPlots package and update break_points in test_break_points.jl --- test/test_break_points.jl | 93 ++++++++++++++++++++++++++++----------- 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/test/test_break_points.jl b/test/test_break_points.jl index 7115cc217..7bb721bb3 100644 --- a/test/test_break_points.jl +++ b/test/test_break_points.jl @@ -1,31 +1,72 @@ -using MacroModelling +using MacroModelling, StatsPlots -function transform_break_points(break_points) - dict = Dict{Int, Dict{Symbol, Float64}}() - - for i in 1:size(break_points,2) - time = axiskeys(break_points,2)[i] +include("../models/Smets_Wouters_2003.jl") - if !haskey(dict, time) - dict[time] = Dict{Symbol, Float64}() - end - - for k in 1:size(break_points,1) - parameter = axiskeys(break_points,1)[k] - value = break_points[k,i] - if !(isnothing(value) || !isfinite(value)) - dict[time][parameter] = value - end - end - end +SSS(Smets_Wouters_2003, parameters = [:std_scaling_factor => 50, :alpha => .3], algorithm = :pruned_second_order) + +# break_points = KeyedArray([.35,.33,.34]', Variable = [:alpha], Time = [3,15,30]) + +break_points = KeyedArray([.37]', Variable = [:alpha], Time = [50]) + +# SSS(Smets_Wouters_2003, algorithm = :pruned_second_order) + +# SS(Smets_Wouters_2003, parameters = [:std_scaling_factor => 10, :alpha => .3]) + + +irfs = get_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, levels = true, algorithm = :pruned_second_order, periods = 40) + + +plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :pruned_third_order, periods = 80) + +plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :third_order, periods = 80) + +plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :pruned_second_order, periods = 80, variables = :q_f) + +plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :pruned_second_order, periods = 80) + + +plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :second_order, periods = 80) + + + +# plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, algorithm = :second_order, periods = 400) + +# plot_irf(Smets_Wouters_2003, parameters = break_points, shocks = :none, periods = 400) + + +# SSS_delta_old = NSSS_old - reference_steady_state_old +# SSS_delta_new = NSSS_new - reference_steady_state_new + +# Ξ”SSS_delta = SSS_delta_old - SSS_delta_new + + +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta_old] +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta_old - Ξ”SSS_delta] +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta_old - (SSS_delta_old - SSS_delta_new)] +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta_old - SSS_delta_old + SSS_delta_new] +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - (SSS_delta + SSS_delta2)] + +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - NSSS + reference_steady_state] +# initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - ((NSSS - reference_steady_state) + (NSSS2 - reference_steady_state2))] + +# state_update = function(pruned_states::Vector{Vector{T}}, shock::Vector{S}) where {T,S} +# aug_state₁ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 1; shock] +# aug_stateβ‚‚ = [pruned_states[2][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] - if !haskey(dict, 0) - dict[0] = Dict{Symbol, Float64}() - end +# return [𝐒₁ * aug_state₁, 𝐒₁ * aug_stateβ‚‚ + 𝐒₂ * β„’.kron(aug_state₁, aug_state₁) / 2] # strictly following Andreasen et al. (2018) +# end + +SSS(Smets_Wouters_2003, parameters = [:std_scaling_factor => 30, :alpha => .35], algorithm = :second_order) +SSS(Smets_Wouters_2003, parameters = [:std_scaling_factor => 30, :alpha => .35], algorithm = :pruned_second_order) + + +# get_parameters(Smets_Wouters_2003, values = true) + +irfsalt = get_irf(Smets_Wouters_2003, parameters = nothing, shocks = :none, levels = true, algorithm = :pruned_second_order, initial_state = collect(irfs[:,1,1])) + +plot_irf(Smets_Wouters_2003, shocks = :none, algorithm = :pruned_second_order, initial_state = collect(irfs[:,1,1])) - return dict -end include("../models/Smets_Wouters_2007.jl") @@ -34,9 +75,11 @@ irfs = get_irf(Smets_Wouters_2007) irfs = get_irf(Smets_Wouters_2007, shocks = :none, levels = true) get_parameters(Smets_Wouters_2007, values = true) -break_points = KeyedArray([.4,.3,.25]', Variable =[:ctrend], Time = [3,15,30]) +break_points = KeyedArray([.4,.3,.25]', Variable =[:ctrend], Time = [1,5,30]) + +irfs = get_irf(Smets_Wouters_2007, parameters = break_points, shocks = :none, levels = true, algorithm = :pruned_second_order) -irfs = get_irf(Smets_Wouters_2007, parameters = break_points, shocks = :none, levels = true) +SSS(Smets_Wouters_2007, algorithm = :pruned_second_order) starting_vals = get_irf(Smets_Wouters_2007, shocks = :none, levels = true, periods = 1) irfalt = get_irf(Smets_Wouters_2007, parameters = :ctrend => .4, shocks = :none, levels = true, initial_state = vec(starting_vals))