diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md index b8b0d1d63..55b83062d 100644 --- a/docs/src/unfinished_docs/todo.md +++ b/docs/src/unfinished_docs/todo.md @@ -64,6 +64,9 @@ - [ ] FastDifferentiation is faster in taking derivatives and more efficient in writing functions but does not support custom functions (e.g. normlogpdf) - [ ] fix this inference errors for large fuctions. they are slow. fix derivatives in general. - [ ] 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 - [ ] put write_derivatives_function and lock structure inside function - [ ] take apart solve_matrix_equation for various cases - [ ] try static arrays in KF diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 3643e4c82..fcdc0d27f 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -503,6 +503,175 @@ function minimize_distance_to_conditions(X::Vector{S}, p)::S where S end +# function match_conditions(res::Vector{S}, X::Vector{S}, jac::Matrix{S}, p) where S +# conditions, state_update, shocks, cond_var_idx, free_shock_idx, state, 𝒷 = p + +# if length(jac) > 0 +# jac .= π’œ.jacobian(𝒷(), xx -> begin +# shocks[free_shock_idx] .= xx +# return abs2.(conditions[cond_var_idx] - state_update(state, convert(typeof(xx), shocks))[cond_var_idx]) +# end, X)[1]' +# end + +# shocks[free_shock_idx] .= X + +# 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} + if state isa Vector{T} + pruning = false + else + pruning = true + end + + if length(grad) > 0 + grad .= π’œ.gradient(𝒷(), xx -> begin + state_copy = deepcopy(state) + + XX = reshape(xx, length(X) Γ· warmup_iters, warmup_iters) + + for i in 1:warmup_iters + state_copy = state_update(state_copy, XX[:,i]) + end + + return precision_factor .* sum(abs2, data - (pruning ? sum(state_copy) : state_copy)[cond_var_idx]) + end, X)[1] + end + + state_copy = deepcopy(state) + + XX = reshape(X, length(X) Γ· warmup_iters, warmup_iters) + + for i in 1:warmup_iters + state_copy = state_update(state_copy, XX[:,i]) + end + + return precision_factor .* sum(abs2, data - (pruning ? sum(state_copy) : state_copy)[cond_var_idx]) +end + + + + +function match_initial_data!(res::Vector{S}, X::Vector{S}, jac::Matrix{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} + if state isa Vector{T} + pruning = false + else + pruning = true + end + + if length(jac) > 0 + jac .= π’œ.jacobian(𝒷(), xx -> begin + state_copy = deepcopy(state) + + XX = reshape(xx, length(X) Γ· warmup_iters, warmup_iters) + + for i in 1:warmup_iters + state_copy = state_update(state_copy, XX[:,i]) + end + + return precision_factor .* abs.(data - (pruning ? sum(state_copy) : state_copy)[cond_var_idx]) + end, X)[1]' + end + + if length(res) > 0 + state_copy = deepcopy(state) + + XX = reshape(X, length(X) Γ· warmup_iters, warmup_iters) + + for i in 1:warmup_iters + state_copy = state_update(state_copy, XX[:,i]) + end + + res .= abs.(data - (pruning ? sum(state_copy) : state_copy)[cond_var_idx]) + end +end + + + +function minimize_distance_to_initial_data(X::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, pruning::Bool)::S where {S, T} + state_copy = deepcopy(state) + + XX = reshape(X, length(X) Γ· warmup_iters, warmup_iters) + + for i in 1:warmup_iters + state_copy = state_update(state_copy, XX[:,i]) + end + + return precision_factor .* sum(abs2, data - (pruning ? sum(state_copy) : state_copy)[cond_var_idx]) +end + + + + + +function minimize_distance_to_data(X::Vector{S}, Data::Vector{T}, State::Union{Vector{T},Vector{Vector{T}}}, state_update::Function, cond_var_idx::Vector{Union{Nothing, Int64}}, precision_factor::Float64, pruning::Bool)::S where {S, T} + return precision_factor .* sum(abs2, Data - (pruning ? sum(state_update(State, X)) : state_update(State, X))[cond_var_idx]) +end + + +function minimize_distance_to_data!(X::Vector{S}, grad::Vector{S}, Data::Vector{T}, State::Union{Vector{T},Vector{Vector{T}}}, state_update::Function, cond_var_idx::Vector{Union{Nothing, Int64}}, precision_factor::Float64) where {S, T} + if State isa Vector{T} + pruning = false + else + pruning = true + end + + if length(grad) > 0 + grad .= π’œ.gradient(𝒷(), xx -> precision_factor .* sum(abs2, Data - (pruning ? sum(state_update(State, xx)) : state_update(State, xx))[cond_var_idx]), X)[1] + end + + return precision_factor .* sum(abs2, Data - (pruning ? sum(state_update(State, X)) : state_update(State, X))[cond_var_idx]) +end + + + +function match_data_sequence!(res::Vector{S}, X::Vector{S}, jac::Matrix{S}, Data::Vector{T}, State::Union{Vector{T},Vector{Vector{T}}}, state_update::Function, cond_var_idx::Vector{Union{Nothing, Int64}}, precision_factor::Float64) where {S, T} + if State isa Vector{T} + pruning = false + else + pruning = true + end + + if length(jac) > 0 + jac .= π’œ.jacobian(𝒷(), xx -> precision_factor .* abs.(Data - (pruning ? sum(state_update(State, xx)) : state_update(State, xx))[cond_var_idx]), X)[1]' + end + + if length(res) > 0 + res .= precision_factor .* abs.(Data - (pruning ? sum(state_update(State, X)) : state_update(State, X))[cond_var_idx]) + end +end + + function set_up_obc_violation_function!(𝓂) present_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚€β‚Ž$"))) @@ -7554,8 +7723,7 @@ function irf(state_update::Function, periods::Int = 40, shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, variables::Union{Symbol_input,String_input} = :all, - shock_size::Real = 1, - negative_shock::Bool = false)::Union{KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{String}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{String}}}} + negative_shock::Bool = false) pruning = initial_state isa Vector{Vector{Float64}} @@ -7589,7 +7757,7 @@ function irf(state_update::Function, shock_idx = parse_shocks_input_to_index(shocks,T) end - var_idx = parse_variables_input_to_index(variables, T) |> sort + var_idx = parse_variables_input_to_index(variables, T) axis1 = T.var[var_idx] @@ -7601,7 +7769,7 @@ function irf(state_update::Function, always_solved = true if shocks == :simulate - shock_history = randn(T.nExo,periods) * shock_size + shock_history = randn(T.nExo,periods) shock_history[contains.(string.(T.exo),"α΅’α΅‡αΆœ"),:] .= 0 @@ -7648,9 +7816,9 @@ function irf(state_update::Function, Y = zeros(T.nVars,periods,length(shock_idx)) for (i,ii) in enumerate(shock_idx) - if shocks βˆ‰ [:simulate, :none] && shocks isa Union{Symbol_input,String_input} + if shocks != :simulate && shocks isa Union{Symbol_input,String_input} shock_history = zeros(T.nExo,periods) - shock_history[ii,1] = negative_shock ? -shock_size : shock_size + shock_history[ii,1] = negative_shock ? -1 : 1 end past_states = initial_state @@ -7688,6 +7856,185 @@ 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) +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) + + periods = size(shock_history, 2) + + Y = zeros(T.nVars, periods, length(shock_idx)) + + pruning = initial_state isa Vector{Vector{Float64}} + + periods_of_change = steady_states_and_state_update |> keys |> collect |> sort + + initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] + + 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) + + reference_steady_state, __, SSS_delta, state_update = steady_states_and_state_update[p] + + if k > 1 + Ξ”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 + initial_state_copy[i][1] += Ξ”NSSS + # initial_state_copy[i][2] += Ξ”NSSS#Ξ”SSS_delta + if length(initial_state_copy[i]) > 3 + initial_state_copy[i][3] += Ξ”NSSS + end + else + 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) + end + end + end + + return Y +end + + + +function girf(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, + warmup_periods::Int = 100, + draws::Int = 50) + + 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 + 1) + + 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) + + Y = zeros(T.nVars, periods + 1, length(shock_idx)) + + for (i,ii) in enumerate(shock_idx) + initial_state_copy = deepcopy(initial_state) + + for draw in 1:draws + initial_state_copyΒ² = deepcopy(initial_state_copy) + + for i in 1:warmup_periods + initial_state_copyΒ² = state_update(initial_state_copyΒ², randn(T.nExo)) + end + + Y₁ = zeros(T.nVars, periods + 1) + Yβ‚‚ = zeros(T.nVars, periods + 1) + + baseline_noise = randn(T.nExo) + + 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 + + if pruning + initial_state_copyΒ² = state_update(initial_state_copyΒ², baseline_noise) + + initial_state₁ = deepcopy(initial_state_copyΒ²) + initial_stateβ‚‚ = deepcopy(initial_state_copyΒ²) + + Y₁[:,1] = initial_state_copyΒ² |> sum + Yβ‚‚[:,1] = initial_state_copyΒ² |> sum + else + Y₁[:,1] = state_update(initial_state_copyΒ², baseline_noise) + Yβ‚‚[:,1] = state_update(initial_state_copyΒ², baseline_noise) + end + + for t in 1:periods + baseline_noise = randn(T.nExo) + + if pruning + initial_state₁ = state_update(initial_state₁, baseline_noise) + initial_stateβ‚‚ = state_update(initial_stateβ‚‚, baseline_noise + shock_history[:,t]) + + Y₁[:,t+1] = initial_state₁ |> sum + Yβ‚‚[:,t+1] = initial_stateβ‚‚ |> sum + else + Y₁[:,t+1] = state_update(Y₁[:,t],baseline_noise) + Yβ‚‚[:,t+1] = state_update(Yβ‚‚[:,t],baseline_noise + shock_history[:,t]) + end + end + + Y[:,:,i] += Yβ‚‚ - Y₁ + end + Y[:,:,i] /= draws + end + + 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 + + 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,2:end,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) +end + + + + function irf(state_update::Function, initial_state::Union{Vector{Vector{Float64}},Vector{Float64}}, level::Vector{Float64}, diff --git a/src/get_functions.jl b/src/get_functions.jl index baf25f97f..b8aa708c8 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1140,7 +1140,7 @@ 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, negative_shock::Bool = false, @@ -1177,42 +1177,98 @@ 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 # end # timeit_debug @@ -1261,16 +1317,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 # @timeit_debug timer "Calculate IRFs" begin girfs = girf(state_update, @@ -1366,17 +1436,14 @@ function get_irf(𝓂::β„³; negative_shock = negative_shock) # end # timeit_debug else - # @timeit_debug timer "Calculate IRFs" begin - 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, - shock_size = shock_size, - negative_shock = negative_shock) - # end # timeit_debug + 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/src/plotting.jl b/src/plotting.jl index 953795ede..db0e78dce 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -142,7 +142,6 @@ function plot_model_estimates(𝓂::β„³, delete!(attributes_redux, :framestyle) - # write_parameters_input!(𝓂, parameters, verbose = verbose) @assert filter ∈ [:kalman, :inversion] "Currently only the kalman filter (:kalman) for linear models and the inversion filter (:inversion) for linear and nonlinear models are supported." @@ -485,7 +484,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, @@ -528,45 +527,108 @@ function plot_irf(𝓂::β„³; 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) |> sort - 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 + 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 + solve!(𝓂, parameters = parameters, opts = opts, dynamics = true, algorithm = algorithm, obc = occasionally_binding_constraints || obc_shocks_included) reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm, opts = opts) @@ -588,7 +650,7 @@ 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 if algorithm βˆ‰ [:pruned_second_order, :pruned_third_order] @@ -596,18 +658,39 @@ function plot_irf(𝓂::β„³; 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, @@ -739,15 +822,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, - shock_size = shock_size, - 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 @@ -777,27 +861,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) diff --git a/test/test_break_points.jl b/test/test_break_points.jl new file mode 100644 index 000000000..7bb721bb3 --- /dev/null +++ b/test/test_break_points.jl @@ -0,0 +1,675 @@ + +using MacroModelling, StatsPlots + +include("../models/Smets_Wouters_2003.jl") + +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)] + +# 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])) + + +include("../models/Smets_Wouters_2007.jl") + +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 = [1,5,30]) + +irfs = get_irf(Smets_Wouters_2007, parameters = break_points, shocks = :none, levels = true, algorithm = :pruned_second_order) + +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)) + + + + + +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 +T = 𝓂.timings +periods = 100 +algorithm = :first_order + +# parameters = nothing + +variables = :all_excluding_obc +shocks = :none # :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 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_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_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),"α΅’α΅‡αΆœ")]) != []) +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) + +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 + 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) + +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 + + +steady_states_and_state_update = Dict{Int, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Function}}() + +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 + + +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) + + +Y = zeros(T.nVars, periods, length(shock_idx)) + +pruning = initial_state isa Vector{Vector{Float64}} + +periods_of_change = steady_states_and_state_update |> keys |> collect |> sort + +initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] + +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) + + println(concerned_periods) + reference_steady_state, NSSS, 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] + + if pruning + for j in initial_state_copy[i] + j += Ξ”SS + end + else + initial_state_copy[i] += Ξ”SS + 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] + end + end +end + + +Y[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] + +(Y .+ reference_steady_state)[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] + +reference_steady_states[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:] + + + +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 + + + + + + + + +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) + +# reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm) + + +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 + 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 1 ∈ concerned_periods + initial_state_copy = [deepcopy(initial_state) for _ in shock_idx] + end + + 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 + + +(Y .+ reference_steady_states)[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:,:] + +reference_steady_states[indexin([:dinve,:pinfobs,:robs,:dwobs,:labobs],T.var),:] + + + + + +# 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) + + 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,:] + + + + + + + + + + + +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 +