Skip to content

Commit

Permalink
obc shock param
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorek1
thorek1 committed Nov 12, 2023
1 parent 92e5bc0 commit 8f90a6c
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 23 deletions.
27 changes: 19 additions & 8 deletions src/MacroModelling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -329,10 +329,21 @@ function set_up_obc_violation_function!(𝓂)
state_update,
reference_steady_state,
𝓂,
algorithm,
periods,
shock_values)
T = 𝓂.timings

# if algorithm ∈ [:first_order, :riccati, :linear_time_iteration, :quadratic_iteration]
# Ytype = JuMP.AffExpr
# elseif algorithm ∈ [:pruned_second_order, :second_order]
# Ytype = JuMP.QuadExpr
# elseif algorithm ∈ [:pruned_third_order, :third_order]
# Ytype = JuMP.NonlinearExpr
# end

# Y = zeros(Ytype, T.nVars, periods+2)

Y = zeros(JuMP.AffExpr, T.nVars, periods+2)

shock_values[contains.(string.(T.exo),"α΅’α΅‡αΆœ")] .= x
Expand Down Expand Up @@ -1058,10 +1069,10 @@ function parse_occasionally_binding_constraints(equations_block; max_obc_shift::
# push!(eqs, :($(obc) = $(Expr(:ref, obc.args[1], -1)) * 0.3 + $(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(max_obc_shift)) * "⁾"), 0))))
push!(eqs, :($(obc) = $(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(max_obc_shift)) * "⁾"), 0))))

push!(eqs, :($(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻⁰⁾"), 0)) = $(Expr(:ref, Meta.parse(string(obc.args[1]) * "⁽" * super(string(max_obc_shift)) * "⁾"), :x))))
push!(eqs, :($(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻⁰⁾"), 0)) = activeα΅’α΅‡αΆœshocks * $(Expr(:ref, Meta.parse(string(obc.args[1]) * "⁽" * super(string(max_obc_shift)) * "⁾"), :x))))

for i in 1:max_obc_shift
push!(eqs, :($(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(i)) * "⁾"), 0)) = $(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(i-1)) * "⁾"), -1)) + $(Expr(:ref, Meta.parse(string(obc.args[1]) * "⁽" * super(string(max_obc_shift-i)) * "⁾"), :x))))
push!(eqs, :($(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(i)) * "⁾"), 0)) = $(Expr(:ref, Meta.parse(string(obc.args[1]) * "ᴸ⁽⁻" * super(string(i-1)) * "⁾"), -1)) + activeα΅’α΅‡αΆœshocks * $(Expr(:ref, Meta.parse(string(obc.args[1]) * "⁽" * super(string(max_obc_shift-i)) * "⁾"), :x))))
end
end

Expand Down Expand Up @@ -4602,14 +4613,14 @@ function irf(state_update::Function,
end
end
else
past_states, past_shocks, solved = obc_state_update(initial_state, zero(shock_history[:,1]), shock_history[:,1])
past_states, past_shocks, solved = obc_state_update(initial_state, zero(shock_history[:,1]), shock_history[:,1], state_update, algorithm)

always_solved = solved
if !solved
@warn "No solution at iteration 1"
else
for t in 2:periods
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shock_history[:,t])
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shock_history[:,t], state_update, algorithm)

if !solved @warn "No solution at iteration $t" end

Expand Down Expand Up @@ -4657,14 +4668,14 @@ function irf(state_update::Function,
end
end
else
past_states, past_shocks, solved = obc_state_update(initial_state, shck, shck)
past_states, past_shocks, solved = obc_state_update(initial_state, shck, shck, state_update, algorithm)

always_solved = solved
if !solved
@warn "No solution at iteration 1"
else
for t in 2:periods
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shck)
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shck, state_update, algorithm)

if !solved @warn "No solution at iteration $t" end
always_solved = always_solved && solved
Expand Down Expand Up @@ -4712,15 +4723,15 @@ function irf(state_update::Function,
end
end
else
past_states, past_shocks, solved = obc_state_update(initial_state, zero(shock_history[:,1]), shock_history[:,1])
past_states, past_shocks, solved = obc_state_update(initial_state, zero(shock_history[:,1]), shock_history[:,1], state_update, algorithm)

always_solved = solved

if !solved
@warn "No solution at iteration 1"
else
for t in 2:periods
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shock_history[:,t])
past_states, past_shocks, solved = obc_state_update(past_states, past_shocks, shock_history[:,t], state_update, algorithm)

if !solved @warn "No solution at iteration $t" end

Expand Down
22 changes: 15 additions & 7 deletions src/get_functions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -841,8 +841,6 @@ function get_irf(𝓂::β„³;

@assert !(shocks == :none && generalised_irf) "Cannot compute generalised IRFs for model without shocks."

state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂)

reference_steady_state, solution_error = 𝓂.solution.outdated_NSSS ? 𝓂.SS_solve_func(𝓂.parameter_values, 𝓂, verbose) : (copy(𝓂.solution.non_stochastic_steady_state), eps())

if algorithm == :second_order
Expand Down Expand Up @@ -879,6 +877,14 @@ function get_irf(𝓂::β„³;
occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0
end

if occasionally_binding_constraints #&&
@assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraints only compatible with first order perturbation solutions."

solve!(𝓂, parameters = :activeα΅’α΅‡αΆœshocks => 1, verbose = false, dynamics = true, algorithm = algorithm)
end

state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂)

if generalised_irf
girfs = girf(state_update,
SSS_delta,
Expand All @@ -894,11 +900,9 @@ function get_irf(𝓂::β„³;
return girfs
else
if occasionally_binding_constraints
function obc_state_update(past_states::Vector{R}, past_shocks::Vector{R}, present_shocks::Vector{R}) where R <: Float64
function obc_state_update(past_states::Vector{R}, past_shocks::Vector{R}, present_shocks::Vector{R}, state_update::Function, algorithm::Symbol) where R <: Float64
unconditional_forecast_horizon = 40

state_update = 𝓂.solution.perturbation.first_order.state_update

reference_steady_state = 𝓂.solution.non_stochastic_steady_state

obc_shock_idx = contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")
Expand All @@ -907,7 +911,7 @@ function get_irf(𝓂::β„³;

num_shocks = sum(obc_shock_idx)Γ·periods_per_shock

constraints_violated = any(JuMP.value.(𝓂.obc_violation_function(zeros(num_shocks*periods_per_shock), past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1]) .> 1e-12)
constraints_violated = any(JuMP.value.(𝓂.obc_violation_function(zeros(num_shocks*periods_per_shock), past_states, past_shocks, state_update, reference_steady_state, 𝓂, algorithm, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1]) .> 1e-12)

if constraints_violated
# Find shocks fulfilling constraint
Expand Down Expand Up @@ -948,7 +952,7 @@ function get_irf(𝓂::β„³;

JuMP.@objective(model, Min, x' * β„’.I * x)

JuMP.@constraint(model, 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1] .<= 0)
JuMP.@constraint(model, 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, algorithm, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1] .<= 0)

JuMP.optimize!(model)

Expand Down Expand Up @@ -994,6 +998,10 @@ function get_irf(𝓂::β„³;
negative_shock = negative_shock)
end

if occasionally_binding_constraints #&& algorithm ∈ [:pruned_second_order, :second_order, :pruned_third_order, :third_order]
solve!(𝓂, parameters = :activeα΅’α΅‡αΆœshocks => 0, verbose = false, dynamics = true, algorithm = algorithm)
end

return irfs
end
end
Expand Down
5 changes: 5 additions & 0 deletions src/macros.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1351,6 +1351,11 @@ macro parameters(𝓂,ex...)
return quote
mod = @__MODULE__

if any(contains.(string.(mod.$𝓂.var), "α΅’α΅‡αΆœ"))
push!($calib_parameters, :activeα΅’α΅‡αΆœshocks)
push!($calib_values, 0)
end

calib_parameters, calib_values = expand_indices($calib_parameters, $calib_values, [mod.$𝓂.parameters_in_equations; mod.$𝓂.var])
calib_eq_parameters, calib_equations_list, ss_calib_list, par_calib_list = expand_calibration_equations($calib_eq_parameters, $calib_equations_list, $ss_calib_list, $par_calib_list, [mod.$𝓂.parameters_in_equations; mod.$𝓂.var])
calib_parameters_no_var, calib_equations_no_var_list = expand_indices($calib_parameters_no_var, $calib_equations_no_var_list, [mod.$𝓂.parameters_in_equations; mod.$𝓂.var])
Expand Down
24 changes: 16 additions & 8 deletions src/plotting.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,8 +378,6 @@ function plot_irf(𝓂::β„³;

solve!(𝓂, parameters = parameters, verbose = verbose, dynamics = true, algorithm = algorithm)

state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂)

NSSS, solution_error = 𝓂.solution.outdated_NSSS ? 𝓂.SS_solve_func(𝓂.parameter_values, 𝓂, verbose) : (𝓂.solution.non_stochastic_steady_state, eps())

full_SS = sort(union(𝓂.var,𝓂.aux,𝓂.exo_present))
Expand Down Expand Up @@ -434,13 +432,21 @@ function plot_irf(𝓂::β„³;
variables = variables isa String_input ? variables .|> Meta.parse .|> replace_indices : variables

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
end

if occasionally_binding_constraints #&&
@assert algorithm βˆ‰ [:pruned_second_order, :second_order, :pruned_third_order, :third_order] "Occasionally binding constraints only compatible with first order perturbation solutions."

solve!(𝓂, parameters = :activeα΅’α΅‡αΆœshocks => 1, verbose = false, dynamics = true, algorithm = algorithm)
end

state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂)

if generalised_irf
Y = girf(state_update,
SSS_delta,
Expand All @@ -455,11 +461,9 @@ function plot_irf(𝓂::β„³;
negative_shock = negative_shock)#, warmup_periods::Int = 100, draws::Int = 50, iterations_to_steady_state::Int = 500)
else
if occasionally_binding_constraints
function obc_state_update(past_states::Vector{R}, past_shocks::Vector{R}, present_shocks::Vector{R}) where R <: Float64
function obc_state_update(past_states::Vector{R}, past_shocks::Vector{R}, present_shocks::Vector{R}, state_update::Function, algorithm::Symbol) where R <: Float64
unconditional_forecast_horizon = 40

state_update = 𝓂.solution.perturbation.first_order.state_update

reference_steady_state = 𝓂.solution.non_stochastic_steady_state

obc_shock_idx = contains.(string.(𝓂.timings.exo),"α΅’α΅‡αΆœ")
Expand All @@ -468,7 +472,7 @@ function plot_irf(𝓂::β„³;

num_shocks = sum(obc_shock_idx)Γ·periods_per_shock

constraints_violated = any(JuMP.value.(𝓂.obc_violation_function(zeros(num_shocks*periods_per_shock), past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1]) .> 1e-12)
constraints_violated = any(JuMP.value.(𝓂.obc_violation_function(zeros(num_shocks*periods_per_shock), past_states, past_shocks, state_update, reference_steady_state, 𝓂, algorithm, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1]) .> 1e-12)

if constraints_violated
# Find shocks fulfilling constraint
Expand Down Expand Up @@ -509,7 +513,7 @@ function plot_irf(𝓂::β„³;

JuMP.@objective(model, Min, x' * β„’.I * x)

JuMP.@constraint(model, 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1] .<= 0)
JuMP.@constraint(model, 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, algorithm, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[1] .<= 0)

JuMP.optimize!(model)

Expand Down Expand Up @@ -556,6 +560,10 @@ function plot_irf(𝓂::β„³;
end
end

if occasionally_binding_constraints #&& algorithm ∈ [:pruned_second_order, :second_order, :pruned_third_order, :third_order]
solve!(𝓂, parameters = :activeα΅’α΅‡αΆœshocks => 0, verbose = false, dynamics = true, algorithm = algorithm)
end

if shocks isa KeyedArray{Float64} || shocks isa Matrix{Float64}
periods += size(shocks)[2]
end
Expand Down

0 comments on commit 8f90a6c

Please sign in to comment.