Skip to content

Commit

Permalink
got 1 period optim going
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorek1
thorek1 committed Oct 17, 2023
1 parent b7fc564 commit d1f1488
Show file tree
Hide file tree
Showing 5 changed files with 500 additions and 28 deletions.
2 changes: 2 additions & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ DynarePreprocessor_jll = "23afba7c-24e5-5ee2-bc2c-b42e07f0492a"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -30,6 +31,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
StatusSwitchingQP = "73a7d86c-03e3-42fb-a0f5-70ef6cda08f4"
Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
Expand Down
37 changes: 28 additions & 9 deletions src/MacroModelling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ import SpecialFunctions: erfcinv, erfc
import SymPyPythonCall as SPyPyC
import Symbolics
import ForwardDiff as ℱ
import JuMP: AffExpr
# import Zygote
import SparseArrays: SparseMatrixCSC, SparseVector, AbstractSparseArray#, sparse, spzeros, droptol!, sparsevec, spdiagm, findnz#, sparse!
import LinearAlgebra as ℒ
Expand Down Expand Up @@ -124,6 +125,24 @@ Base.show(io::IO, 𝓂::ℳ) = println(io,



function parse_obc_shock_bounds(expr::Expr)
# Determine the order of the shock and bound in the expression
shock_first = isa(expr.args[2], Symbol)

# Extract the shock and bound from the expression
shock = shock_first ? expr.args[2] : expr.args[3]
bound_expr = shock_first ? expr.args[3] : expr.args[2]

# Evaluate the bound expression to get a numerical value
bound = eval(bound_expr) |> Float64

# Determine whether the bound is a lower or upper bound based on the comparison operator and order
is_upper_bound = (expr.args[1] in (:<, :≤) && shock_first) || (expr.args[1] in (:>, :≥) && !shock_first)

return shock, is_upper_bound, bound
end



function mat_mult_kron(A::AbstractArray{T},B::AbstractArray{T},C::AbstractArray{T}; tol::AbstractFloat = eps()) where T <: Real
n_rowB = size(B,1)
Expand Down Expand Up @@ -601,7 +620,7 @@ function match_pattern(strings::Union{Set,Vector}, pattern::Regex)
return filter(r -> match(pattern, string(r)) !== nothing, strings)
end

function parse_occasionally_binding_constraints(equations_block; max_obc_shift::Int = 20)
function parse_occasionally_binding_constraints(equations_block; max_obc_shift::Int = 10)
eqs = []
condition_list = []

Expand Down Expand Up @@ -646,7 +665,7 @@ function parse_occasionally_binding_constraints(equations_block; max_obc_shift::
x.args[1] ∈ obc_shocks ?
begin
obc_shock = intersect([x.args[1]], obc_shocks)[1]
obc_shifts = [Expr(:ref,Meta.parse(string(obc_shock) * "ᵒᵇᶜ⁽⁻"*super(string(i))*"⁾"),:(x-$i)) for i in 0:max_obc_shift]
obc_shifts = [Expr(:ref,Meta.parse(string(obc_shock) * "ᵒᵇᶜ⁽⁻"*super(string(i))*"⁾"),i > 0 ? :(x - $i) : :x) for i in 0:max_obc_shift]
Expr(:call,:+, x, obc_shifts...)
end :
x :
Expand Down Expand Up @@ -2675,7 +2694,7 @@ function solve!(𝓂::ℳ;

@assert solved "Could not find stable first order solution."

state_update₁ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end
state_update₁ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end

𝓂.solution.perturbation.first_order = perturbation_solution(sol_mat, state_update₁)
𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:riccati, :first_order])
Expand All @@ -2694,7 +2713,7 @@ function solve!(𝓂::ℳ;

@assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1."

state_update₂ = function(state::Vector{Float64}, shock::Vector{Float64})
state_update₂ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}})
aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx]
1
shock]
Expand All @@ -2715,7 +2734,7 @@ function solve!(𝓂::ℳ;

@assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1."

state_update₂ = function(pruned_states::Vector{Vector{Float64}}, shock::Vector{Float64})
state_update₂ = function(pruned_states::Vector{Vector{<: Union{Real,AffExpr}}}, shock::Vector{<: Union{Real,AffExpr}})
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)]

Expand All @@ -2736,7 +2755,7 @@ function solve!(𝓂::ℳ;

@assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1."

state_update₃ = function(state::Vector{Float64}, shock::Vector{Float64})
state_update₃ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}})
aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx]
1
shock]
Expand All @@ -2754,7 +2773,7 @@ function solve!(𝓂::ℳ;

@assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1."

state_update₃ = function(pruned_states::Vector{Vector{Float64}}, shock::Vector{Float64})
state_update₃ = function(pruned_states::Vector{Vector{<: Union{Real,AffExpr}}}, shock::Vector{<: Union{Real,AffExpr}})
aug_state₁ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 1; shock]
aug_state₁̂ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 0; shock]
aug_state₂ = [pruned_states[2][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)]
Expand Down Expand Up @@ -2786,7 +2805,7 @@ function solve!(𝓂::ℳ;

sol_mat, converged = calculate_quadratic_iteration_solution(∇₁; T = 𝓂.timings)

state_update₁ₜ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end
state_update₁ₜ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end

𝓂.solution.perturbation.quadratic_iteration = perturbation_solution(sol_mat, state_update₁ₜ)
𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:quadratic_iteration, :binder_pesaran])
Expand All @@ -2807,7 +2826,7 @@ function solve!(𝓂::ℳ;

sol_mat = calculate_linear_time_iteration_solution(∇₁; T = 𝓂.timings)

state_update₁ₜ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end
state_update₁ₜ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end

𝓂.solution.perturbation.linear_time_iteration = perturbation_solution(sol_mat, state_update₁ₜ)
𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:linear_time_iteration])
Expand Down
10 changes: 10 additions & 0 deletions src/macros.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,14 @@ macro model(𝓂,ex...)

model_ex, condition_list = parse_occasionally_binding_constraints(model_ex)

obc_shock_bounds = Tuple{Symbol, Bool, Float64}[]

for c in condition_list
if c isa Expr
push!(obc_shock_bounds, parse_obc_shock_bounds(c))
end
end

# write down dynamic equations and add auxilliary variables for leads and lags > 1
for (i,arg) in enumerate(model_ex.args)
if isa(arg,Expr)
Expand Down Expand Up @@ -849,6 +857,8 @@ macro model(𝓂,ex...)

$T,

$obc_shock_bounds,

solution(
perturbation( perturbation_solution(SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), x->x),
perturbation_solution(SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), x->x),
Expand Down
3 changes: 3 additions & 0 deletions src/structures.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -349,6 +349,9 @@ mutable struct ℳ
model_third_order_derivatives::Vector{Function}

timings::timings

obc_shock_bounds::Vector{Tuple{Symbol, Bool, Float64}}

solution::solution
# symbolics::symbolics

Expand Down
Loading

0 comments on commit d1f1488

Please sign in to comment.