Skip to content

Commit

Permalink
implement complementarity
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorek1
thorek1 committed Nov 11, 2023
1 parent c7ceda1 commit 2771550
Show file tree
Hide file tree
Showing 4 changed files with 142 additions and 110 deletions.
3 changes: 2 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
Expand Down Expand Up @@ -58,8 +59,8 @@ Krylov = "^0.9"
LaTeXStrings = "^1"
LinearOperators = "^2"
MacroTools = "^0.5"
MadNLP = "^0.7"
MatrixEquations = "^2"
NLopt = "^1"
PrecompileTools = "^1"
RecursiveFactorization = "^0.2"
Reexport = "^1"
Expand Down
51 changes: 32 additions & 19 deletions src/MacroModelling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ import SymPyPythonCall as SPyPyC
import Symbolics
import ForwardDiff as ℱ
import JuMP
import NLopt # MadNLP doesnt support noninear constraints # StatusSwitchingQP not reliable
import MadNLP # NLopt # MadNLP doesnt support noninear constraints # StatusSwitchingQP not reliable
# import Zygote
import SparseArrays: SparseMatrixCSC, SparseVector, AbstractSparseArray#, sparse, spzeros, droptol!, sparsevec, spdiagm, findnz#, sparse!
import LinearAlgebra as ℒ
Expand Down Expand Up @@ -355,7 +355,7 @@ function set_up_obc_violation_function!(𝓂)
$(steady_state...)
$(steady_state_obc...)

constraint_values = Vector{JuMP.AffExpr}[]
constraint_values = Vector[]
shock_sign_indicators = Bool[]

$(𝓂.obc_violation_equations...)
Expand Down Expand Up @@ -415,24 +415,37 @@ function write_obc_violation_equations(𝓂)

maximisation = contains(string(plchldr), "⁺")

if dyn_1
if maximisation
push!(cond1, :(push!(constraint_values, $(x.args[3].args[2]))))
# push!(cond2, :(push!(constraint_values, $(x.args[3].args[2]))))
else
push!(cond1, :(push!(constraint_values, -$(x.args[3].args[2]))))
# push!(cond2, :(push!(constraint_values, -$(x.args[3].args[2])))) # RBC
end
end
# if dyn_1
# if maximisation
# push!(cond1, :(push!(constraint_values, $(x.args[3].args[2]))))
# # push!(cond2, :(push!(constraint_values, $(x.args[3].args[2]))))
# else
# push!(cond1, :(push!(constraint_values, -$(x.args[3].args[2]))))
# # push!(cond2, :(push!(constraint_values, -$(x.args[3].args[2])))) # RBC
# end
# end

# if dyn_2
# if maximisation
# push!(cond1, :(push!(constraint_values, $(x.args[3].args[3]))))
# # push!(cond2, :(push!(constraint_values, $(x.args[3].args[3])))) # testmax
# else
# push!(cond1, :(push!(constraint_values, -$(x.args[3].args[3]))))
# # push!(cond2, :(push!(constraint_values, -$(x.args[3].args[3])))) # RBC
# end
# end


if dyn_2
if maximisation
push!(cond1, :(push!(constraint_values, $(x.args[3].args[3]))))
# push!(cond2, :(push!(constraint_values, $(x.args[3].args[3])))) # testmax
else
push!(cond1, :(push!(constraint_values, -$(x.args[3].args[3]))))
# push!(cond2, :(push!(constraint_values, -$(x.args[3].args[3])))) # RBC
end
if maximisation
push!(cond1, :(push!(constraint_values, [sum($(x.args[3].args[2]) .* $(x.args[3].args[3]))])))
push!(cond1, :(push!(constraint_values, $(x.args[3].args[2]))))
push!(cond1, :(push!(constraint_values, $(x.args[3].args[3]))))
# push!(cond1, :(push!(constraint_values, max.($(x.args[3].args[2]), $(x.args[3].args[3])))))
else
push!(cond1, :(push!(constraint_values, [sum($(x.args[3].args[2]) .* $(x.args[3].args[3]))])))
push!(cond1, :(push!(constraint_values, -$(x.args[3].args[2]))))
push!(cond1, :(push!(constraint_values, -$(x.args[3].args[3]))))
# push!(cond1, :(push!(constraint_values, min.($(x.args[3].args[2]), $(x.args[3].args[3])))))
end

if maximisation
Expand Down
99 changes: 54 additions & 45 deletions src/get_functions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -907,54 +907,63 @@ function get_irf(𝓂::ℳ;

num_shocks = sum(obc_shock_idx)÷periods_per_shock

# Find shocks fulfilling constraint
# model = JuMP.Model(MadNLP.Optimizer)
model = JuMP.Model(NLopt.Optimizer)
JuMP.set_attribute(model, "algorithm", :LD_SLSQP)
# JuMP.set_attribute(model, "algorithm", :LD_MMA)

JuMP.set_silent(model)

# JuMP.set_attribute(model, "tol", 1e-12)

# Create the variables over the full set of indices first.
JuMP.@variable(model, x[1:num_shocks*periods_per_shock])

# Now loop through obc_shock_bounds to set the bounds on these variables.
# maxmin_indicators = 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2]
# for (idx, v) in enumerate(maxmin_indicators)
# idxs = (idx - 1) * periods_per_shock + 1:idx * periods_per_shock
# if v
# # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# JuMP.set_upper_bound.(x[idxs], 0)
# # JuMP.set_lower_bound.(x[idxs], 0)
# else
# # JuMP.set_upper_bound.(x[idxs], 0)
# JuMP.set_lower_bound.(x[idxs], 0)
# end
# # # else
# # # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# # # JuMP.set_lower_bound.(x[idxs], 0)
# # # else
# # # JuMP.set_upper_bound.(x[idxs], 0)
# # # end
# # # end
# end

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.optimize!(model)
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)

solved = JuMP.termination_status(model) ∈ [JuMP.OPTIMAL,JuMP.LOCALLY_SOLVED]

# precision = JuMP.objective_value(model)

# if precision > eps(Float32) @warn "Bounds enforced up to reduced precision: $precision" end # I need the dual value (constraints). this relates to the shock size
if constraints_violated
# Find shocks fulfilling constraint
model = JuMP.Model(MadNLP.Optimizer)
# model = JuMP.Model(NLopt.Optimizer)
# JuMP.set_attribute(model, "algorithm", :LD_SLSQP)
# JuMP.set_attribute(model, "algorithm", :AUGLAG)
# JuMP.set_attribute(model, "local_optimizer", :LD_LBFGS)
# JuMP.set_attribute(model, "algorithm", :LD_MMA)

JuMP.set_silent(model)

# JuMP.set_attribute(model, "tol", 1e-12)

# Create the variables over the full set of indices first.
JuMP.@variable(model, x[1:num_shocks*periods_per_shock])

# Now loop through obc_shock_bounds to set the bounds on these variables.
# maxmin_indicators = 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2]
# for (idx, v) in enumerate(maxmin_indicators)
# idxs = (idx - 1) * periods_per_shock + 1:idx * periods_per_shock
# if v
# # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# JuMP.set_upper_bound.(x[idxs], 0)
# # JuMP.set_lower_bound.(x[idxs], 0)
# else
# # JuMP.set_upper_bound.(x[idxs], 0)
# JuMP.set_lower_bound.(x[idxs], 0)
# end
# # # else
# # # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# # # JuMP.set_lower_bound.(x[idxs], 0)
# # # else
# # # JuMP.set_upper_bound.(x[idxs], 0)
# # # end
# # # end
# end

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.optimize!(model)

solved = JuMP.termination_status(model) ∈ [JuMP.OPTIMAL,JuMP.LOCALLY_SOLVED]

# precision = JuMP.objective_value(model)

# if precision > eps(Float32) @warn "Bounds enforced up to reduced precision: $precision" end # I need the dual value (constraints). this relates to the shock size

present_shocks[contains.(string.(𝓂.timings.exo),"ᵒᵇᶜ")] .= JuMP.value.(x)
else
solved = true
end

present_states = state_update(past_states,JuMP.value.(past_shocks))
present_shocks[contains.(string.(𝓂.timings.exo),"ᵒᵇᶜ")] .= JuMP.value.(x)

return present_states, present_shocks, solved
end
Expand Down
99 changes: 54 additions & 45 deletions src/plotting.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -468,54 +468,63 @@ function plot_irf(𝓂::ℳ;

num_shocks = sum(obc_shock_idx)÷periods_per_shock

# Find shocks fulfilling constraint
# model = JuMP.Model(MadNLP.Optimizer)
model = JuMP.Model(NLopt.Optimizer)
JuMP.set_attribute(model, "algorithm", :LD_SLSQP)
# JuMP.set_attribute(model, "algorithm", :LD_MMA)

JuMP.set_silent(model)

# JuMP.set_attribute(model, "tol", 1e-12)

# Create the variables over the full set of indices first.
JuMP.@variable(model, x[1:num_shocks*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)

# Now loop through obc_shock_bounds to set the bounds on these variables.
# maxmin_indicators = 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2]
# for (idx, v) in enumerate(maxmin_indicators)
# idxs = (idx - 1) * periods_per_shock + 1:idx * periods_per_shock
# if v
# # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# JuMP.set_upper_bound.(x[idxs], 0)
# # JuMP.set_lower_bound.(x[idxs], 0)
# else
# # JuMP.set_upper_bound.(x[idxs], 0)
# JuMP.set_lower_bound.(x[idxs], 0)
# end
# # # else
# # # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# # # JuMP.set_lower_bound.(x[idxs], 0)
# # # else
# # # JuMP.set_upper_bound.(x[idxs], 0)
# # # end
# # # end
# end

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.optimize!(model)

solved = JuMP.termination_status(model) ∈ [JuMP.OPTIMAL,JuMP.LOCALLY_SOLVED]

# precision = JuMP.objective_value(model)

# if precision > eps(Float32) @warn "Bounds enforced up to reduced precision: $precision" end # I need the dual value (constraints). this relates to the shock size
if constraints_violated
# Find shocks fulfilling constraint
model = JuMP.Model(MadNLP.Optimizer)
# model = JuMP.Model(NLopt.Optimizer)
# JuMP.set_attribute(model, "algorithm", :LD_SLSQP)
# JuMP.set_attribute(model, "algorithm", :AUGLAG)
# JuMP.set_attribute(model, "local_optimizer", :LD_LBFGS)
# JuMP.set_attribute(model, "algorithm", :LD_MMA)

JuMP.set_silent(model)

# JuMP.set_attribute(model, "tol", 1e-12)

# Create the variables over the full set of indices first.
JuMP.@variable(model, x[1:num_shocks*periods_per_shock])

# Now loop through obc_shock_bounds to set the bounds on these variables.
# maxmin_indicators = 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2]
# for (idx, v) in enumerate(maxmin_indicators)
# idxs = (idx - 1) * periods_per_shock + 1:idx * periods_per_shock
# if v
# # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# JuMP.set_upper_bound.(x[idxs], 0)
# # JuMP.set_lower_bound.(x[idxs], 0)
# else
# # JuMP.set_upper_bound.(x[idxs], 0)
# JuMP.set_lower_bound.(x[idxs], 0)
# end
# # # else
# # # if 𝓂.obc_violation_function(x, past_states, past_shocks, state_update, reference_steady_state, 𝓂, unconditional_forecast_horizon, JuMP.AffExpr.(present_shocks))[2][idx]
# # # JuMP.set_lower_bound.(x[idxs], 0)
# # # else
# # # JuMP.set_upper_bound.(x[idxs], 0)
# # # end
# # # end
# end

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.optimize!(model)

solved = JuMP.termination_status(model) ∈ [JuMP.OPTIMAL,JuMP.LOCALLY_SOLVED]

# precision = JuMP.objective_value(model)

# if precision > eps(Float32) @warn "Bounds enforced up to reduced precision: $precision" end # I need the dual value (constraints). this relates to the shock size

present_shocks[contains.(string.(𝓂.timings.exo),"ᵒᵇᶜ")] .= JuMP.value.(x)
else
solved = true
end

present_states = state_update(past_states,JuMP.value.(past_shocks))
present_shocks[contains.(string.(𝓂.timings.exo),"ᵒᵇᶜ")] .= JuMP.value.(x)

return present_states, present_shocks, solved
end
Expand Down

0 comments on commit 2771550

Please sign in to comment.