diff --git a/Project.toml b/Project.toml index 6b1ad0cf..14a2b697 100644 --- a/Project.toml +++ b/Project.toml @@ -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" @@ -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" diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 46e75248..375f92b9 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -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 β„’ @@ -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...) @@ -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 diff --git a/src/get_functions.jl b/src/get_functions.jl index 82ff05c4..cd14f58d 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -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 diff --git a/src/plotting.jl b/src/plotting.jl index 7f7c87f4..8477d961 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -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