From 19b1595638dba44709120b6585fd2d7d5abf6261 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 08:31:55 +0000 Subject: [PATCH 1/7] Initial plan From 5bb6de419102a16adc697d5a2d41785288e2f5fa Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 08:40:15 +0000 Subject: [PATCH 2/7] =?UTF-8?q?Add=20Pad=C3=A9=20approximation=20algorithm?= =?UTF-8?q?s=20for=20higher-order=20solutions?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 380 ++++++++++++++++++++++++++++++++++++++- src/common_docstrings.jl | 2 +- src/get_functions.jl | 70 +++++++- src/macros.jl | 6 +- src/structures.jl | 4 + 5 files changed, 448 insertions(+), 14 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 7578ae5c2..62dc3bac4 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -548,9 +548,9 @@ function normalize_filtering_options(filter::Symbol, maxlog::Int = DEFAULT_MAXLOG) @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." - pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order) + pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade) - if shock_decomposition && algorithm ∈ (:second_order, :third_order) + if shock_decomposition && algorithm ∈ (:second_order, :third_order, :second_order_pade, :third_order_pade) @info "Shock decomposition is not available for $(algorithm) solutions, but is available for first order, pruned second order, and pruned third order solutions. Setting `shock_decomposition = false`." maxlog = maxlog shock_decomposition = false end @@ -877,7 +877,7 @@ function set_up_obc_violation_function!(𝓂) zero_shock = zero(shock_values) - if algorithm ∈ [:pruned_second_order, :pruned_third_order] + if algorithm ∈ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] states = state_update(state, shock_values) Y[:,1] = sum(states) else @@ -885,7 +885,7 @@ function set_up_obc_violation_function!(𝓂) end for t in 1:periods - if algorithm ∈ [:pruned_second_order, :pruned_third_order] + if algorithm ∈ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] states = state_update(states, zero_shock) Y[:,t+1] = sum(states) else @@ -6963,12 +6963,12 @@ function solve!(𝓂::ℳ; # end # timeit_debug if 𝓂.solution.perturbation.second_order_auxiliary_matrices.𝛔 == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && - algorithm ∈ [:second_order, :pruned_second_order] + algorithm ∈ [:second_order, :pruned_second_order, :second_order_pade, :pruned_second_order_pade] start_time = time() if !silent print("Take symbolic derivatives up to second order:\t\t\t\t") end write_functions_mapping!(𝓂, 2) if !silent println(round(time() - start_time, digits = 3), " seconds") end - elseif 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐂₃ == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && algorithm ∈ [:third_order, :pruned_third_order] + elseif 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐂₃ == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade, :pruned_third_order_pade] start_time = time() if !silent print("Take symbolic derivatives up to third order:\t\t\t\t") end write_functions_mapping!(𝓂, 3) @@ -6981,7 +6981,11 @@ function solve!(𝓂::ℳ; ((:second_order == algorithm) && ((:second_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || ((:pruned_second_order == algorithm) && ((:pruned_second_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || ((:third_order == algorithm) && ((:third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:pruned_third_order == algorithm) && ((:pruned_third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + ((:pruned_third_order == algorithm) && ((:pruned_third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:pruned_second_order_pade == algorithm) && ((:pruned_second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) # @timeit_debug timer "Solve for NSSS (if necessary)" begin @@ -7183,6 +7187,310 @@ function solve!(𝓂::ℳ; 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:pruned_third_order]) end + + # Second order Padé approximation (unpruned) + obc_not_solved = isnothing(𝓂.solution.perturbation.second_order_pade.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) + if ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + + stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts) + + if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end + + # Padé [1,1] approximation for second order: y = S₁*x + (1/2)*S₂*kron(x,x) + # Approximated as: y ≈ S₁*x ./ (1 .- D*x) where D is derived from second-order terms + # For efficiency, we use element-wise Padé: y_i = (S₁*x)_i / (1 - d_i * ||x||²) + # where d_i captures the curvature from S₂ for variable i + state_update₂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] + 1 + shock] + + linear_term = 𝐒₁ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + + # Padé [1,1] approximation: convert Taylor to rational form + # For f(x) ≈ a + bx + cx², Padé [1,1] is: (a + (b - ac/b)x) / (1 - (c/b)x) + # In matrix form, we approximate element-wise + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i]) > eps(T) + # Padé coefficient: ratio of quadratic to linear contribution + pade_coeff = quadratic_term[i] / linear_term[i] + denominator[i] = 1 - pade_coeff + end + end + + # Apply Padé approximation: y = linear / denominator + return linear_term ./ denominator + end + + if obc + Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + + state_update₂̂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] + 1 + shock] + + linear_term = Ŝ₁̂ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i]) > eps(T) + pade_coeff = quadratic_term[i] / linear_term[i] + denominator[i] = 1 - pade_coeff + end + end + + return linear_term ./ denominator + end + else + state_update₂̂_pade = (x,y)->nothing + end + + 𝓂.solution.perturbation.second_order_pade = second_order_perturbation_solution(stochastic_steady_state, state_update₂_pade, state_update₂̂_pade) + + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:second_order_pade]) + end + + # Pruned second order Padé approximation + obc_not_solved = isnothing(𝓂.solution.perturbation.pruned_second_order_pade.state_update_obc([zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars)], zeros(𝓂.timings.nExo))) + if ((:pruned_second_order_pade == algorithm) && ((:pruned_second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + + stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts, pruning = true) + + if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end + + state_update₂_pade = 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)] + + # First order part (linear) + y₁ = 𝐒₁ * aug_state₁ + + # Second order part with Padé approximation + linear_term = 𝐒₁ * aug_state₂ + kron_aug = ℒ.kron(aug_state₁, aug_state₁) + quadratic_term = 𝐒₂ * kron_aug / 2 + + # Apply Padé to second order contribution + y₂_combined = linear_term + quadratic_term + denominator = ones(T, length(y₂_combined)) + for i in eachindex(y₂_combined) + if abs(linear_term[i] + quadratic_term[i]) > eps(T) && abs(linear_term[i]) > eps(T) + pade_coeff = quadratic_term[i] / (linear_term[i] + quadratic_term[i]) + denominator[i] = 1 - pade_coeff + end + end + + return [y₁, y₂_combined ./ denominator] + end + + if obc + Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + + state_update₂̂_pade = 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)] + + y₁ = Ŝ₁̂ * aug_state₁ + + linear_term = Ŝ₁̂ * aug_state₂ + kron_aug = ℒ.kron(aug_state₁, aug_state₁) + quadratic_term = 𝐒₂ * kron_aug / 2 + + y₂_combined = linear_term + quadratic_term + denominator = ones(T, length(y₂_combined)) + for i in eachindex(y₂_combined) + if abs(linear_term[i] + quadratic_term[i]) > eps(T) && abs(linear_term[i]) > eps(T) + pade_coeff = quadratic_term[i] / (linear_term[i] + quadratic_term[i]) + denominator[i] = 1 - pade_coeff + end + end + + return [y₁, y₂_combined ./ denominator] + end + else + state_update₂̂_pade = (x,y)->nothing + end + + 𝓂.solution.perturbation.pruned_second_order_pade = second_order_perturbation_solution(stochastic_steady_state, state_update₂_pade, state_update₂̂_pade) + + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:pruned_second_order_pade]) + end + + # Third order Padé approximation (unpruned) + obc_not_solved = isnothing(𝓂.solution.perturbation.third_order_pade.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) + if ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts) + + if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end + + # Padé [1,2] or [2,1] approximation for third order + state_update₃_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] + 1 + shock] + + linear_term = 𝐒₁ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 + + # Padé [2,1] approximation for third order Taylor series + # f(x) ≈ a₀ + a₁x + a₂x² + a₃x³ → (b₀ + b₁x + b₂x²) / (1 + c₁x) + taylor_sum = linear_term + quadratic_term + cubic_term + + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i] + quadratic_term[i]) > eps(T) + # Use cubic term to estimate denominator coefficient + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / (linear_term[i] + quadratic_term[i]) + denominator[i] = 1 - higher_order_ratio + end + end + + return taylor_sum ./ denominator + end + + if obc + Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + + state_update₃̂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] + 1 + shock] + + linear_term = Ŝ₁̂ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 + + taylor_sum = linear_term + quadratic_term + cubic_term + + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i] + quadratic_term[i]) > eps(T) + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / (linear_term[i] + quadratic_term[i]) + denominator[i] = 1 - higher_order_ratio + end + end + + return taylor_sum ./ denominator + end + else + state_update₃̂_pade = (x,y)->nothing + end + + 𝓂.solution.perturbation.third_order_pade = third_order_perturbation_solution(stochastic_steady_state, state_update₃_pade, state_update₃̂_pade) + + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:third_order_pade]) + end + + # Pruned third order Padé approximation + obc_not_solved = isnothing(𝓂.solution.perturbation.pruned_third_order_pade.state_update_obc([zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars)], zeros(𝓂.timings.nExo))) + if ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + + stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts, pruning = true) + + if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end + + state_update₃_pade = 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[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)] + aug_state₃ = [pruned_states[3][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] + + kron_aug_state₁ = ℒ.kron(aug_state₁, aug_state₁) + + # First order part (no Padé) + y₁ = 𝐒₁ * aug_state₁ + + # Second order part with Padé + y₂_linear = 𝐒₁ * aug_state₂ + y₂_quad = 𝐒₂ * kron_aug_state₁ / 2 + y₂_combined = y₂_linear + y₂_quad + + denom₂ = ones(T, length(y₂_combined)) + for i in eachindex(y₂_combined) + if abs(y₂_linear[i] + y₂_quad[i]) > eps(T) && abs(y₂_linear[i]) > eps(T) + denom₂[i] = 1 - y₂_quad[i] / (y₂_linear[i] + y₂_quad[i]) + end + end + y₂ = y₂_combined ./ denom₂ + + # Third order part with Padé + y₃_linear = 𝐒₁ * aug_state₃ + y₃_quad = 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + y₃_cubic = 𝐒₃ * ℒ.kron(kron_aug_state₁, aug_state₁) / 6 + y₃_combined = y₃_linear + y₃_quad + y₃_cubic + + denom₃ = ones(T, length(y₃_combined)) + for i in eachindex(y₃_combined) + if abs(y₃_linear[i] + y₃_quad[i]) > eps(T) + higher_order_ratio = (y₃_quad[i] + y₃_cubic[i]) / (y₃_linear[i] + y₃_quad[i]) + denom₃[i] = 1 - higher_order_ratio + end + end + y₃ = y₃_combined ./ denom₃ + + return [y₁, y₂, y₃] + end + + if obc + Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + + state_update₃̂_pade = 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[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)] + aug_state₃ = [pruned_states[3][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] + + kron_aug_state₁ = ℒ.kron(aug_state₁, aug_state₁) + + y₁ = Ŝ₁̂ * aug_state₁ + + y₂_linear = Ŝ₁̂ * aug_state₂ + y₂_quad = 𝐒₂ * kron_aug_state₁ / 2 + y₂_combined = y₂_linear + y₂_quad + + denom₂ = ones(T, length(y₂_combined)) + for i in eachindex(y₂_combined) + if abs(y₂_linear[i] + y₂_quad[i]) > eps(T) && abs(y₂_linear[i]) > eps(T) + denom₂[i] = 1 - y₂_quad[i] / (y₂_linear[i] + y₂_quad[i]) + end + end + y₂ = y₂_combined ./ denom₂ + + y₃_linear = Ŝ₁̂ * aug_state₃ + y₃_quad = 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + y₃_cubic = 𝐒₃ * ℒ.kron(kron_aug_state₁, aug_state₁) / 6 + y₃_combined = y₃_linear + y₃_quad + y₃_cubic + + denom₃ = ones(T, length(y₃_combined)) + for i in eachindex(y₃_combined) + if abs(y₃_linear[i] + y₃_quad[i]) > eps(T) + higher_order_ratio = (y₃_quad[i] + y₃_cubic[i]) / (y₃_linear[i] + y₃_quad[i]) + denom₃[i] = 1 - higher_order_ratio + end + end + y₃ = y₃_combined ./ denom₃ + + return [y₁, y₂, y₃] + end + else + state_update₃̂_pade = (x,y)->nothing + end + + 𝓂.solution.perturbation.pruned_third_order_pade = third_order_perturbation_solution(stochastic_steady_state, state_update₃_pade, state_update₃̂_pade) + + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:pruned_third_order_pade]) + end end return nothing @@ -9692,6 +10000,18 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona elseif :pruned_third_order == algorithm state_update = 𝓂.solution.perturbation.pruned_third_order.state_update_obc pruning = true + elseif :second_order_pade == algorithm + state_update = 𝓂.solution.perturbation.second_order_pade.state_update_obc + pruning = false + elseif :pruned_second_order_pade == algorithm + state_update = 𝓂.solution.perturbation.pruned_second_order_pade.state_update_obc + pruning = true + elseif :third_order_pade == algorithm + state_update = 𝓂.solution.perturbation.third_order_pade.state_update_obc + pruning = false + elseif :pruned_third_order_pade == algorithm + state_update = 𝓂.solution.perturbation.pruned_third_order_pade.state_update_obc + pruning = true else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -9713,6 +10033,18 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona elseif :pruned_third_order == algorithm state_update = 𝓂.solution.perturbation.pruned_third_order.state_update pruning = true + elseif :second_order_pade == algorithm + state_update = 𝓂.solution.perturbation.second_order_pade.state_update + pruning = false + elseif :pruned_second_order_pade == algorithm + state_update = 𝓂.solution.perturbation.pruned_second_order_pade.state_update + pruning = true + elseif :third_order_pade == algorithm + state_update = 𝓂.solution.perturbation.third_order_pade.state_update + pruning = false + elseif :pruned_third_order_pade == algorithm + state_update = 𝓂.solution.perturbation.pruned_third_order_pade.state_update + pruning = true else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -10270,6 +10602,40 @@ function get_relevant_steady_state_and_state_update(::Val{:first_order}, return TT, SS_and_pars, 𝐒₁, [state], solved end +# Padé versions of get_relevant_steady_state_and_state_update +# These use the same stochastic steady state calculations as their Taylor counterparts +function get_relevant_steady_state_and_state_update(::Val{:second_order_pade}, + parameter_values::Vector{S}, + 𝓂::ℳ; + opts::CalculationOptions = merge_calculation_options()) where S <: Real + # Use the same calculation as second_order + return get_relevant_steady_state_and_state_update(Val(:second_order), parameter_values, 𝓂, opts = opts) +end + +function get_relevant_steady_state_and_state_update(::Val{:pruned_second_order_pade}, + parameter_values::Vector{S}, + 𝓂::ℳ; + opts::CalculationOptions = merge_calculation_options())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real + # Use the same calculation as pruned_second_order + return get_relevant_steady_state_and_state_update(Val(:pruned_second_order), parameter_values, 𝓂, opts = opts) +end + +function get_relevant_steady_state_and_state_update(::Val{:third_order_pade}, + parameter_values::Vector{S}, + 𝓂::ℳ; + opts::CalculationOptions = merge_calculation_options())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{S}, Bool} where S <: Real + # Use the same calculation as third_order + return get_relevant_steady_state_and_state_update(Val(:third_order), parameter_values, 𝓂, opts = opts) +end + +function get_relevant_steady_state_and_state_update(::Val{:pruned_third_order_pade}, + parameter_values::Vector{S}, + 𝓂::ℳ; + opts::CalculationOptions = merge_calculation_options())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real + # Use the same calculation as pruned_third_order + return get_relevant_steady_state_and_state_update(Val(:pruned_third_order), parameter_values, 𝓂, opts = opts) +end + end # dispatch_doctor # @setup_workload begin diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl index b0b2ef5ab..7e898000a 100644 --- a/src/common_docstrings.jl +++ b/src/common_docstrings.jl @@ -12,7 +12,7 @@ const NEGATIVE_SHOCK® = "`negative_shock` [Default: `$(DEFAULT_NEGATIVE_SHOCK)` const GENERALISED_IRF® = "`generalised_irf` [Default: `$(DEFAULT_GENERALISED_IRF)`, Type: `Bool`]: calculate generalised IRFs. Relevant for nonlinear (higher order perturbation) solutions only. Reference steady state for deviations is the stochastic steady state. `initial_state` has no effect on generalised IRFs. Occasionally binding constraint are not respected for generalised IRF." const GENERALISED_IRF_WARMUP_ITERATIONS® = "`generalised_irf_warmup_iterations` [Default: `$(DEFAULT_GENERALISED_IRF_WARMUP)`, Type: `Int`]: number of warm-up iterations used to draw the baseline paths in the generalised IRF simulation. Only applied when `generalised_irf = true`." const GENERALISED_IRF_DRAWS® = "`generalised_irf_draws` [Default: `$(DEFAULT_GENERALISED_IRF_DRAWS)`, Type: `Int`]: number of Monte Carlo draws used to compute the generalised IRF. Only applied when `generalised_irf = true`." -const ALGORITHM® = "`algorithm` [Default: `$(DEFAULT_ALGORITHM)`, Type: `Symbol`]: algorithm to solve for the dynamics of the model. Available algorithms: `:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`" +const ALGORITHM® = "`algorithm` [Default: `$(DEFAULT_ALGORITHM)`, Type: `Symbol`]: algorithm to solve for the dynamics of the model. Available algorithms: `:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`, `:second_order_pade`, `:pruned_second_order_pade`, `:third_order_pade`, `:pruned_third_order_pade`. The Padé variants use Padé approximation instead of Taylor series for the state update function, which can provide better numerical properties for large deviations from steady state." const FILTER® = "`filter` [Default: selector that chooses `$(DEFAULT_FILTER_SELECTOR(DEFAULT_ALGORITHM))` in case `algorithm = $(DEFAULT_ALGORITHM)` and `:inversion` otherwise, Type: `Symbol`]: filter used to compute the variables and shocks given the data, model, and parameters. The Kalman filter only works for linear problems, whereas the inversion filter (`:inversion`) works for linear and nonlinear models. If a nonlinear solution algorithm is selected and the default is used, the inversion filter is applied automatically." const LEVELS® = "return levels or absolute deviations from the relevant steady state corresponding to the solution algorithm (e.g. stochastic steady state for higher order solution algorithms)." const CONDITIONS® = "`conditions` [Type: `Union{Matrix{Union{Nothing,Float64}}, SparseMatrixCSC{Float64}, KeyedArray{Union{Nothing,Float64}}, KeyedArray{Float64}}`]: conditions for which to find the corresponding shocks. The input can have multiple formats, but for all types of entries, the first dimension corresponds to variables and the second dimension to the number of periods. The conditions can be specified using a matrix of type `Matrix{Union{Nothing,Float64}}`. In this case the conditions are matrix elements of type `Float64` and all remaining (free) entries are `nothing`. A `SparseMatrixCSC{Float64}` can also be used as input. In this case only non-zero elements are taken as conditions. Note that conditioning variables to be zero using a `SparseMatrixCSC{Float64}` as input is not possible (use other input formats to do so). Another possibility to input conditions is by using a `KeyedArray`. The `KeyedArray` type is provided by the `AxisKeys` package. A `KeyedArray{Union{Nothing,Float64}}` can be used where, similar to `Matrix{Union{Nothing,Float64}}`, all entries of type `Float64` are recognised as conditions and all other entries have to be `nothing`. Furthermore, in the primary axis a subset of variables (of type `Symbol` or `String`) for which conditions are specified can be included and all other variables are considered free. The same goes for the case when using `KeyedArray{Float64}}` as input, whereas in this case the conditions for the specified variables bind for all periods specified in the `KeyedArray`, because there are no `nothing` entries permitted with this type." diff --git a/src/get_functions.jl b/src/get_functions.jl index 156de4cd1..d9ac142bf 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -847,6 +847,10 @@ function get_conditional_forecast(𝓂::ℳ, 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)] + elseif algorithm == :pruned_second_order_pade + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order_pade + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] else initial_state = zeros(𝓂.timings.nVars) - SSS_delta end @@ -856,11 +860,15 @@ function get_conditional_forecast(𝓂::ℳ, 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)] + elseif algorithm == :pruned_second_order_pade + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order_pade + 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] + if algorithm ∉ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." end end @@ -884,7 +892,7 @@ function get_conditional_forecast(𝓂::ℳ, @assert length(free_shock_idx) >= length(cond_var_idx) "Exact matching only possible with at least as many free shocks than conditioned variables. Period 1 has " * repr(length(free_shock_idx)) * " free shock(s) and " * repr(length(cond_var_idx)) * " conditioned variable(s)." - if algorithm ∈ [:second_order, :third_order, :pruned_second_order, :pruned_third_order] + if algorithm ∈ [:second_order, :third_order, :pruned_second_order, :pruned_third_order, :second_order_pade, :third_order_pade, :pruned_second_order_pade, :pruned_third_order_pade] S₁ = 𝓂.solution.perturbation.first_order.solution_matrix S₁ = [S₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) S₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] @@ -894,7 +902,7 @@ function get_conditional_forecast(𝓂::ℳ, end S₃ = nothing - if algorithm ∈ [:third_order, :pruned_third_order] && size(𝓂.solution.perturbation.third_order_solution, 2) > 0 + if algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade, :pruned_third_order_pade] && size(𝓂.solution.perturbation.third_order_solution, 2) > 0 S₃ = 𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃ end @@ -1296,6 +1304,10 @@ function get_irf(𝓂::ℳ; 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)] + elseif algorithm == :pruned_second_order_pade + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order_pade + initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] else initial_state = zeros(𝓂.timings.nVars) - SSS_delta end @@ -1305,11 +1317,15 @@ function get_irf(𝓂::ℳ; 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)] + elseif algorithm == :pruned_second_order_pade + initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] + elseif algorithm == :pruned_third_order_pade + 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] + if algorithm ∉ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] @assert initial_state isa Vector{Float64} "The solution algorithm has one state vector: initial_state must be a Vector{Float64}." end end @@ -1318,7 +1334,7 @@ function get_irf(𝓂::ℳ; 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." + @assert algorithm ∉ [:pruned_second_order, :second_order, :pruned_third_order, :third_order, :pruned_second_order_pade, :second_order_pade, :pruned_third_order_pade, :third_order_pade] "Occasionally binding constraint shocks without enforcing the constraint is only compatible with first order perturbation solutions." state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) else @@ -1817,6 +1833,26 @@ function get_solution(𝓂::ℳ; States__Shocks¹ = axis1, Variables = axis2, States__Shocks² = axis1) + elseif algorithm == :second_order_pade + # Padé uses same solution matrices, transformation is in state_update + return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.second_order_solution * 𝓂.solution.perturbation.second_order_auxiliary_matrices.𝐔₂, + 𝓂.timings.nVars, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), + [2,1,3]); + States__Shocks¹ = axis1, + Variables = axis2, + States__Shocks² = axis1) + elseif algorithm == :pruned_second_order_pade + # Padé uses same solution matrices, transformation is in state_update + return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.second_order_solution * 𝓂.solution.perturbation.second_order_auxiliary_matrices.𝐔₂, + 𝓂.timings.nVars, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), + [2,1,3]); + States__Shocks¹ = axis1, + Variables = axis2, + States__Shocks² = axis1) elseif algorithm == :third_order return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃, 𝓂.timings.nVars, @@ -1839,6 +1875,30 @@ function get_solution(𝓂::ℳ; Variables = axis2, States__Shocks² = axis1, States__Shocks³ = axis1) + elseif algorithm == :third_order_pade + # Padé uses same solution matrices, transformation is in state_update + return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃, + 𝓂.timings.nVars, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), + [2,1,3,4]); + States__Shocks¹ = axis1, + Variables = axis2, + States__Shocks² = axis1, + States__Shocks³ = axis1) + elseif algorithm == :pruned_third_order_pade + # Padé uses same solution matrices, transformation is in state_update + return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃, + 𝓂.timings.nVars, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, + 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), + [2,1,3,4]); + States__Shocks¹ = axis1, + Variables = axis2, + States__Shocks² = axis1, + States__Shocks³ = axis1) else axis1 = [:Steady_state; 𝓂.timings.past_not_future_and_mixed; 𝓂.exo] diff --git a/src/macros.jl b/src/macros.jl index b1abd04e8..ebadfec3e 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -1,4 +1,4 @@ -const all_available_algorithms = [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] +const all_available_algorithms = [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order, :second_order_pade, :pruned_second_order_pade, :third_order_pade, :pruned_third_order_pade] """ @@ -966,6 +966,10 @@ macro model(𝓂,ex...) solution( perturbation( perturbation_solution(zeros(0,0), (x,y)->nothing, (x,y)->nothing), + second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), + second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), + third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), + third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), diff --git a/src/structures.jl b/src/structures.jl index bdd600b90..be3c837f1 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -216,6 +216,10 @@ mutable struct perturbation pruned_second_order::second_order_perturbation_solution third_order::third_order_perturbation_solution pruned_third_order::third_order_perturbation_solution + second_order_pade::second_order_perturbation_solution + pruned_second_order_pade::second_order_perturbation_solution + third_order_pade::third_order_perturbation_solution + pruned_third_order_pade::third_order_perturbation_solution qme_solution::Matrix{Float64} second_order_solution::AbstractMatrix{Float64} third_order_solution::AbstractMatrix{Float64} From e8ece7d2ff533867f020ee61b01a63cba31b5946 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 08:57:07 +0000 Subject: [PATCH 3/7] Changes before error encountered Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 62dc3bac4..dc6ee1eee 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -7213,12 +7213,18 @@ function solve!(𝓂::ℳ; # Padé [1,1] approximation: convert Taylor to rational form # For f(x) ≈ a + bx + cx², Padé [1,1] is: (a + (b - ac/b)x) / (1 - (c/b)x) # In matrix form, we approximate element-wise + # Use a more robust threshold to avoid numerical instability + tol = 1e-10 denominator = ones(T, length(linear_term)) for i in eachindex(linear_term) - if abs(linear_term[i]) > eps(T) + if abs(linear_term[i]) > tol # Padé coefficient: ratio of quadratic to linear contribution pade_coeff = quadratic_term[i] / linear_term[i] - denominator[i] = 1 - pade_coeff + denom_val = 1 - pade_coeff + # Ensure denominator doesn't become too small or negative + if abs(denom_val) > tol + denominator[i] = denom_val + end end end @@ -7238,11 +7244,15 @@ function solve!(𝓂::ℳ; kron_aug = ℒ.kron(aug_state, aug_state) quadratic_term = 𝐒₂ * kron_aug / 2 + tol = 1e-10 denominator = ones(T, length(linear_term)) for i in eachindex(linear_term) - if abs(linear_term[i]) > eps(T) + if abs(linear_term[i]) > tol pade_coeff = quadratic_term[i] / linear_term[i] - denominator[i] = 1 - pade_coeff + denom_val = 1 - pade_coeff + if abs(denom_val) > tol + denominator[i] = denom_val + end end end @@ -7279,12 +7289,17 @@ function solve!(𝓂::ℳ; quadratic_term = 𝐒₂ * kron_aug / 2 # Apply Padé to second order contribution + # Use a robust threshold to avoid numerical instability + tol = 1e-10 y₂_combined = linear_term + quadratic_term denominator = ones(T, length(y₂_combined)) for i in eachindex(y₂_combined) - if abs(linear_term[i] + quadratic_term[i]) > eps(T) && abs(linear_term[i]) > eps(T) - pade_coeff = quadratic_term[i] / (linear_term[i] + quadratic_term[i]) - denominator[i] = 1 - pade_coeff + if abs(y₂_combined[i]) > tol && abs(linear_term[i]) > tol + pade_coeff = quadratic_term[i] / y₂_combined[i] + denom_val = 1 - pade_coeff + if abs(denom_val) > tol + denominator[i] = denom_val + end end end @@ -7304,12 +7319,16 @@ function solve!(𝓂::ℳ; kron_aug = ℒ.kron(aug_state₁, aug_state₁) quadratic_term = 𝐒₂ * kron_aug / 2 + tol = 1e-10 y₂_combined = linear_term + quadratic_term denominator = ones(T, length(y₂_combined)) for i in eachindex(y₂_combined) - if abs(linear_term[i] + quadratic_term[i]) > eps(T) && abs(linear_term[i]) > eps(T) - pade_coeff = quadratic_term[i] / (linear_term[i] + quadratic_term[i]) - denominator[i] = 1 - pade_coeff + if abs(y₂_combined[i]) > tol && abs(linear_term[i]) > tol + pade_coeff = quadratic_term[i] / y₂_combined[i] + denom_val = 1 - pade_coeff + if abs(denom_val) > tol + denominator[i] = denom_val + end end end From 800b8a99d2b64f37427e2a0b9bffa30a33d5eea5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 11:12:15 +0000 Subject: [PATCH 4/7] =?UTF-8?q?Remove=20pruned=20Pad=C3=A9=20variants=20an?= =?UTF-8?q?d=20improve=20numerical=20stability?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 245 ++++----------------------------------- src/common_docstrings.jl | 2 +- src/get_functions.jl | 48 +------- src/macros.jl | 4 +- src/structures.jl | 2 - 5 files changed, 32 insertions(+), 269 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index dc6ee1eee..596f80e37 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -548,7 +548,7 @@ function normalize_filtering_options(filter::Symbol, maxlog::Int = DEFAULT_MAXLOG) @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." - pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade) + pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order) if shock_decomposition && algorithm ∈ (:second_order, :third_order, :second_order_pade, :third_order_pade) @info "Shock decomposition is not available for $(algorithm) solutions, but is available for first order, pruned second order, and pruned third order solutions. Setting `shock_decomposition = false`." maxlog = maxlog @@ -877,7 +877,7 @@ function set_up_obc_violation_function!(𝓂) zero_shock = zero(shock_values) - if algorithm ∈ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] + if algorithm ∈ [:pruned_second_order, :pruned_third_order] states = state_update(state, shock_values) Y[:,1] = sum(states) else @@ -885,7 +885,7 @@ function set_up_obc_violation_function!(𝓂) end for t in 1:periods - if algorithm ∈ [:pruned_second_order, :pruned_third_order, :pruned_second_order_pade, :pruned_third_order_pade] + if algorithm ∈ [:pruned_second_order, :pruned_third_order] states = state_update(states, zero_shock) Y[:,t+1] = sum(states) else @@ -6963,12 +6963,12 @@ function solve!(𝓂::ℳ; # end # timeit_debug if 𝓂.solution.perturbation.second_order_auxiliary_matrices.𝛔 == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && - algorithm ∈ [:second_order, :pruned_second_order, :second_order_pade, :pruned_second_order_pade] + algorithm ∈ [:second_order, :pruned_second_order, :second_order_pade] start_time = time() if !silent print("Take symbolic derivatives up to second order:\t\t\t\t") end write_functions_mapping!(𝓂, 2) if !silent println(round(time() - start_time, digits = 3), " seconds") end - elseif 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐂₃ == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade, :pruned_third_order_pade] + elseif 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐂₃ == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade] start_time = time() if !silent print("Take symbolic derivatives up to third order:\t\t\t\t") end write_functions_mapping!(𝓂, 3) @@ -6982,10 +6982,8 @@ function solve!(𝓂::ℳ; ((:pruned_second_order == algorithm) && ((:pruned_second_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || ((:third_order == algorithm) && ((:third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || ((:pruned_third_order == algorithm) && ((:pruned_third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:pruned_second_order_pade == algorithm) && ((:pruned_second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || + ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) # @timeit_debug timer "Solve for NSSS (if necessary)" begin @@ -7267,82 +7265,6 @@ function solve!(𝓂::ℳ; 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:second_order_pade]) end - # Pruned second order Padé approximation - obc_not_solved = isnothing(𝓂.solution.perturbation.pruned_second_order_pade.state_update_obc([zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars)], zeros(𝓂.timings.nExo))) - if ((:pruned_second_order_pade == algorithm) && ((:pruned_second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) - - stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts, pruning = true) - - if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - - state_update₂_pade = 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)] - - # First order part (linear) - y₁ = 𝐒₁ * aug_state₁ - - # Second order part with Padé approximation - linear_term = 𝐒₁ * aug_state₂ - kron_aug = ℒ.kron(aug_state₁, aug_state₁) - quadratic_term = 𝐒₂ * kron_aug / 2 - - # Apply Padé to second order contribution - # Use a robust threshold to avoid numerical instability - tol = 1e-10 - y₂_combined = linear_term + quadratic_term - denominator = ones(T, length(y₂_combined)) - for i in eachindex(y₂_combined) - if abs(y₂_combined[i]) > tol && abs(linear_term[i]) > tol - pade_coeff = quadratic_term[i] / y₂_combined[i] - denom_val = 1 - pade_coeff - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - return [y₁, y₂_combined ./ denominator] - end - - if obc - Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] - - state_update₂̂_pade = 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)] - - y₁ = Ŝ₁̂ * aug_state₁ - - linear_term = Ŝ₁̂ * aug_state₂ - kron_aug = ℒ.kron(aug_state₁, aug_state₁) - quadratic_term = 𝐒₂ * kron_aug / 2 - - tol = 1e-10 - y₂_combined = linear_term + quadratic_term - denominator = ones(T, length(y₂_combined)) - for i in eachindex(y₂_combined) - if abs(y₂_combined[i]) > tol && abs(linear_term[i]) > tol - pade_coeff = quadratic_term[i] / y₂_combined[i] - denom_val = 1 - pade_coeff - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - return [y₁, y₂_combined ./ denominator] - end - else - state_update₂̂_pade = (x,y)->nothing - end - - 𝓂.solution.perturbation.pruned_second_order_pade = second_order_perturbation_solution(stochastic_steady_state, state_update₂_pade, state_update₂̂_pade) - - 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:pruned_second_order_pade]) - end - # Third order Padé approximation (unpruned) obc_not_solved = isnothing(𝓂.solution.perturbation.third_order_pade.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) if ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) @@ -7350,7 +7272,7 @@ function solve!(𝓂::ℳ; if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Padé [1,2] or [2,1] approximation for third order + # Padé [2,1] approximation for third order state_update₃_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] 1 @@ -7363,14 +7285,20 @@ function solve!(𝓂::ℳ; # Padé [2,1] approximation for third order Taylor series # f(x) ≈ a₀ + a₁x + a₂x² + a₃x³ → (b₀ + b₁x + b₂x²) / (1 + c₁x) + # Use robust threshold to avoid numerical instability + tol = 1e-10 taylor_sum = linear_term + quadratic_term + cubic_term denominator = ones(T, length(linear_term)) for i in eachindex(linear_term) - if abs(linear_term[i] + quadratic_term[i]) > eps(T) + base = linear_term[i] + quadratic_term[i] + if abs(base) > tol # Use cubic term to estimate denominator coefficient - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / (linear_term[i] + quadratic_term[i]) - denominator[i] = 1 - higher_order_ratio + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base + denom_val = 1 - higher_order_ratio + if abs(denom_val) > tol + denominator[i] = denom_val + end end end @@ -7390,13 +7318,18 @@ function solve!(𝓂::ℳ; quadratic_term = 𝐒₂ * kron_aug / 2 cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 + tol = 1e-10 taylor_sum = linear_term + quadratic_term + cubic_term denominator = ones(T, length(linear_term)) for i in eachindex(linear_term) - if abs(linear_term[i] + quadratic_term[i]) > eps(T) - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / (linear_term[i] + quadratic_term[i]) - denominator[i] = 1 - higher_order_ratio + base = linear_term[i] + quadratic_term[i] + if abs(base) > tol + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base + denom_val = 1 - higher_order_ratio + if abs(denom_val) > tol + denominator[i] = denom_val + end end end @@ -7410,106 +7343,6 @@ function solve!(𝓂::ℳ; 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:third_order_pade]) end - - # Pruned third order Padé approximation - obc_not_solved = isnothing(𝓂.solution.perturbation.pruned_third_order_pade.state_update_obc([zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars)], zeros(𝓂.timings.nExo))) - if ((:pruned_third_order_pade == algorithm) && ((:pruned_third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) - - stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts, pruning = true) - - if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - - state_update₃_pade = 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[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)] - aug_state₃ = [pruned_states[3][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] - - kron_aug_state₁ = ℒ.kron(aug_state₁, aug_state₁) - - # First order part (no Padé) - y₁ = 𝐒₁ * aug_state₁ - - # Second order part with Padé - y₂_linear = 𝐒₁ * aug_state₂ - y₂_quad = 𝐒₂ * kron_aug_state₁ / 2 - y₂_combined = y₂_linear + y₂_quad - - denom₂ = ones(T, length(y₂_combined)) - for i in eachindex(y₂_combined) - if abs(y₂_linear[i] + y₂_quad[i]) > eps(T) && abs(y₂_linear[i]) > eps(T) - denom₂[i] = 1 - y₂_quad[i] / (y₂_linear[i] + y₂_quad[i]) - end - end - y₂ = y₂_combined ./ denom₂ - - # Third order part with Padé - y₃_linear = 𝐒₁ * aug_state₃ - y₃_quad = 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) - y₃_cubic = 𝐒₃ * ℒ.kron(kron_aug_state₁, aug_state₁) / 6 - y₃_combined = y₃_linear + y₃_quad + y₃_cubic - - denom₃ = ones(T, length(y₃_combined)) - for i in eachindex(y₃_combined) - if abs(y₃_linear[i] + y₃_quad[i]) > eps(T) - higher_order_ratio = (y₃_quad[i] + y₃_cubic[i]) / (y₃_linear[i] + y₃_quad[i]) - denom₃[i] = 1 - higher_order_ratio - end - end - y₃ = y₃_combined ./ denom₃ - - return [y₁, y₂, y₃] - end - - if obc - Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] - - state_update₃̂_pade = 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[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)] - aug_state₃ = [pruned_states[3][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] - - kron_aug_state₁ = ℒ.kron(aug_state₁, aug_state₁) - - y₁ = Ŝ₁̂ * aug_state₁ - - y₂_linear = Ŝ₁̂ * aug_state₂ - y₂_quad = 𝐒₂ * kron_aug_state₁ / 2 - y₂_combined = y₂_linear + y₂_quad - - denom₂ = ones(T, length(y₂_combined)) - for i in eachindex(y₂_combined) - if abs(y₂_linear[i] + y₂_quad[i]) > eps(T) && abs(y₂_linear[i]) > eps(T) - denom₂[i] = 1 - y₂_quad[i] / (y₂_linear[i] + y₂_quad[i]) - end - end - y₂ = y₂_combined ./ denom₂ - - y₃_linear = Ŝ₁̂ * aug_state₃ - y₃_quad = 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) - y₃_cubic = 𝐒₃ * ℒ.kron(kron_aug_state₁, aug_state₁) / 6 - y₃_combined = y₃_linear + y₃_quad + y₃_cubic - - denom₃ = ones(T, length(y₃_combined)) - for i in eachindex(y₃_combined) - if abs(y₃_linear[i] + y₃_quad[i]) > eps(T) - higher_order_ratio = (y₃_quad[i] + y₃_cubic[i]) / (y₃_linear[i] + y₃_quad[i]) - denom₃[i] = 1 - higher_order_ratio - end - end - y₃ = y₃_combined ./ denom₃ - - return [y₁, y₂, y₃] - end - else - state_update₃̂_pade = (x,y)->nothing - end - - 𝓂.solution.perturbation.pruned_third_order_pade = third_order_perturbation_solution(stochastic_steady_state, state_update₃_pade, state_update₃̂_pade) - - 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:pruned_third_order_pade]) - end end return nothing @@ -10022,15 +9855,9 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona elseif :second_order_pade == algorithm state_update = 𝓂.solution.perturbation.second_order_pade.state_update_obc pruning = false - elseif :pruned_second_order_pade == algorithm - state_update = 𝓂.solution.perturbation.pruned_second_order_pade.state_update_obc - pruning = true elseif :third_order_pade == algorithm state_update = 𝓂.solution.perturbation.third_order_pade.state_update_obc pruning = false - elseif :pruned_third_order_pade == algorithm - state_update = 𝓂.solution.perturbation.pruned_third_order_pade.state_update_obc - pruning = true else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -10055,15 +9882,9 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona elseif :second_order_pade == algorithm state_update = 𝓂.solution.perturbation.second_order_pade.state_update pruning = false - elseif :pruned_second_order_pade == algorithm - state_update = 𝓂.solution.perturbation.pruned_second_order_pade.state_update - pruning = true elseif :third_order_pade == algorithm state_update = 𝓂.solution.perturbation.third_order_pade.state_update pruning = false - elseif :pruned_third_order_pade == algorithm - state_update = 𝓂.solution.perturbation.pruned_third_order_pade.state_update - pruning = true else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -10631,14 +10452,6 @@ function get_relevant_steady_state_and_state_update(::Val{:second_order_pade}, return get_relevant_steady_state_and_state_update(Val(:second_order), parameter_values, 𝓂, opts = opts) end -function get_relevant_steady_state_and_state_update(::Val{:pruned_second_order_pade}, - parameter_values::Vector{S}, - 𝓂::ℳ; - opts::CalculationOptions = merge_calculation_options())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real - # Use the same calculation as pruned_second_order - return get_relevant_steady_state_and_state_update(Val(:pruned_second_order), parameter_values, 𝓂, opts = opts) -end - function get_relevant_steady_state_and_state_update(::Val{:third_order_pade}, parameter_values::Vector{S}, 𝓂::ℳ; @@ -10647,14 +10460,6 @@ function get_relevant_steady_state_and_state_update(::Val{:third_order_pade}, return get_relevant_steady_state_and_state_update(Val(:third_order), parameter_values, 𝓂, opts = opts) end -function get_relevant_steady_state_and_state_update(::Val{:pruned_third_order_pade}, - parameter_values::Vector{S}, - 𝓂::ℳ; - opts::CalculationOptions = merge_calculation_options())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real - # Use the same calculation as pruned_third_order - return get_relevant_steady_state_and_state_update(Val(:pruned_third_order), parameter_values, 𝓂, opts = opts) -end - end # dispatch_doctor # @setup_workload begin diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl index 7e898000a..a0a9f60df 100644 --- a/src/common_docstrings.jl +++ b/src/common_docstrings.jl @@ -12,7 +12,7 @@ const NEGATIVE_SHOCK® = "`negative_shock` [Default: `$(DEFAULT_NEGATIVE_SHOCK)` const GENERALISED_IRF® = "`generalised_irf` [Default: `$(DEFAULT_GENERALISED_IRF)`, Type: `Bool`]: calculate generalised IRFs. Relevant for nonlinear (higher order perturbation) solutions only. Reference steady state for deviations is the stochastic steady state. `initial_state` has no effect on generalised IRFs. Occasionally binding constraint are not respected for generalised IRF." const GENERALISED_IRF_WARMUP_ITERATIONS® = "`generalised_irf_warmup_iterations` [Default: `$(DEFAULT_GENERALISED_IRF_WARMUP)`, Type: `Int`]: number of warm-up iterations used to draw the baseline paths in the generalised IRF simulation. Only applied when `generalised_irf = true`." const GENERALISED_IRF_DRAWS® = "`generalised_irf_draws` [Default: `$(DEFAULT_GENERALISED_IRF_DRAWS)`, Type: `Int`]: number of Monte Carlo draws used to compute the generalised IRF. Only applied when `generalised_irf = true`." -const ALGORITHM® = "`algorithm` [Default: `$(DEFAULT_ALGORITHM)`, Type: `Symbol`]: algorithm to solve for the dynamics of the model. Available algorithms: `:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`, `:second_order_pade`, `:pruned_second_order_pade`, `:third_order_pade`, `:pruned_third_order_pade`. The Padé variants use Padé approximation instead of Taylor series for the state update function, which can provide better numerical properties for large deviations from steady state." +const ALGORITHM® = "`algorithm` [Default: `$(DEFAULT_ALGORITHM)`, Type: `Symbol`]: algorithm to solve for the dynamics of the model. Available algorithms: `:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`, `:second_order_pade`, `:third_order_pade`. The Padé variants use Padé approximation instead of Taylor series for the state update function, which can provide better numerical properties for large deviations from steady state." const FILTER® = "`filter` [Default: selector that chooses `$(DEFAULT_FILTER_SELECTOR(DEFAULT_ALGORITHM))` in case `algorithm = $(DEFAULT_ALGORITHM)` and `:inversion` otherwise, Type: `Symbol`]: filter used to compute the variables and shocks given the data, model, and parameters. The Kalman filter only works for linear problems, whereas the inversion filter (`:inversion`) works for linear and nonlinear models. If a nonlinear solution algorithm is selected and the default is used, the inversion filter is applied automatically." const LEVELS® = "return levels or absolute deviations from the relevant steady state corresponding to the solution algorithm (e.g. stochastic steady state for higher order solution algorithms)." const CONDITIONS® = "`conditions` [Type: `Union{Matrix{Union{Nothing,Float64}}, SparseMatrixCSC{Float64}, KeyedArray{Union{Nothing,Float64}}, KeyedArray{Float64}}`]: conditions for which to find the corresponding shocks. The input can have multiple formats, but for all types of entries, the first dimension corresponds to variables and the second dimension to the number of periods. The conditions can be specified using a matrix of type `Matrix{Union{Nothing,Float64}}`. In this case the conditions are matrix elements of type `Float64` and all remaining (free) entries are `nothing`. A `SparseMatrixCSC{Float64}` can also be used as input. In this case only non-zero elements are taken as conditions. Note that conditioning variables to be zero using a `SparseMatrixCSC{Float64}` as input is not possible (use other input formats to do so). Another possibility to input conditions is by using a `KeyedArray`. The `KeyedArray` type is provided by the `AxisKeys` package. A `KeyedArray{Union{Nothing,Float64}}` can be used where, similar to `Matrix{Union{Nothing,Float64}}`, all entries of type `Float64` are recognised as conditions and all other entries have to be `nothing`. Furthermore, in the primary axis a subset of variables (of type `Symbol` or `String`) for which conditions are specified can be included and all other variables are considered free. The same goes for the case when using `KeyedArray{Float64}}` as input, whereas in this case the conditions for the specified variables bind for all periods specified in the `KeyedArray`, because there are no `nothing` entries permitted with this type." diff --git a/src/get_functions.jl b/src/get_functions.jl index d9ac142bf..75bf007e5 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -847,10 +847,6 @@ function get_conditional_forecast(𝓂::ℳ, 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)] - elseif algorithm == :pruned_second_order_pade - initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] - elseif algorithm == :pruned_third_order_pade - initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] else initial_state = zeros(𝓂.timings.nVars) - SSS_delta end @@ -860,15 +856,11 @@ function get_conditional_forecast(𝓂::ℳ, 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)] - elseif algorithm == :pruned_second_order_pade - initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] - elseif algorithm == :pruned_third_order_pade - 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, :pruned_second_order_pade, :pruned_third_order_pade] + 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 @@ -892,7 +884,7 @@ function get_conditional_forecast(𝓂::ℳ, @assert length(free_shock_idx) >= length(cond_var_idx) "Exact matching only possible with at least as many free shocks than conditioned variables. Period 1 has " * repr(length(free_shock_idx)) * " free shock(s) and " * repr(length(cond_var_idx)) * " conditioned variable(s)." - if algorithm ∈ [:second_order, :third_order, :pruned_second_order, :pruned_third_order, :second_order_pade, :third_order_pade, :pruned_second_order_pade, :pruned_third_order_pade] + if algorithm ∈ [:second_order, :third_order, :pruned_second_order, :pruned_third_order, :second_order_pade, :third_order_pade] S₁ = 𝓂.solution.perturbation.first_order.solution_matrix S₁ = [S₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) S₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] @@ -902,7 +894,7 @@ function get_conditional_forecast(𝓂::ℳ, end S₃ = nothing - if algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade, :pruned_third_order_pade] && size(𝓂.solution.perturbation.third_order_solution, 2) > 0 + if algorithm ∈ [:third_order, :pruned_third_order, :third_order_pade] && size(𝓂.solution.perturbation.third_order_solution, 2) > 0 S₃ = 𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃ end @@ -1304,10 +1296,6 @@ function get_irf(𝓂::ℳ; 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)] - elseif algorithm == :pruned_second_order_pade - initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] - elseif algorithm == :pruned_third_order_pade - initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta, zeros(𝓂.timings.nVars)] else initial_state = zeros(𝓂.timings.nVars) - SSS_delta end @@ -1317,15 +1305,11 @@ function get_irf(𝓂::ℳ; 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)] - elseif algorithm == :pruned_second_order_pade - initial_state = [initial_state - reference_steady_state[1:𝓂.timings.nVars], zeros(𝓂.timings.nVars) - SSS_delta] - elseif algorithm == :pruned_third_order_pade - 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, :pruned_second_order_pade, :pruned_third_order_pade] + 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 @@ -1334,7 +1318,7 @@ function get_irf(𝓂::ℳ; 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, :pruned_second_order_pade, :second_order_pade, :pruned_third_order_pade, :third_order_pade] "Occasionally binding constraint shocks without enforcing the constraint is only compatible with first order perturbation solutions." + @assert algorithm ∉ [:pruned_second_order, :second_order, :pruned_third_order, :third_order, :second_order_pade, :third_order_pade] "Occasionally binding constraint shocks without enforcing the constraint is only compatible with first order perturbation solutions." state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) else @@ -1843,16 +1827,6 @@ function get_solution(𝓂::ℳ; States__Shocks¹ = axis1, Variables = axis2, States__Shocks² = axis1) - elseif algorithm == :pruned_second_order_pade - # Padé uses same solution matrices, transformation is in state_update - return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.second_order_solution * 𝓂.solution.perturbation.second_order_auxiliary_matrices.𝐔₂, - 𝓂.timings.nVars, - 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, - 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), - [2,1,3]); - States__Shocks¹ = axis1, - Variables = axis2, - States__Shocks² = axis1) elseif algorithm == :third_order return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃, 𝓂.timings.nVars, @@ -1887,18 +1861,6 @@ function get_solution(𝓂::ℳ; Variables = axis2, States__Shocks² = axis1, States__Shocks³ = axis1) - elseif algorithm == :pruned_third_order_pade - # Padé uses same solution matrices, transformation is in state_update - return KeyedArray(permutedims(reshape(𝓂.solution.perturbation.third_order_solution * 𝓂.solution.perturbation.third_order_auxiliary_matrices.𝐔₃, - 𝓂.timings.nVars, - 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, - 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo, - 𝓂.timings.nPast_not_future_and_mixed + 1 + 𝓂.timings.nExo), - [2,1,3,4]); - States__Shocks¹ = axis1, - Variables = axis2, - States__Shocks² = axis1, - States__Shocks³ = axis1) else axis1 = [:Steady_state; 𝓂.timings.past_not_future_and_mixed; 𝓂.exo] diff --git a/src/macros.jl b/src/macros.jl index ebadfec3e..498aef97d 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -1,4 +1,4 @@ -const all_available_algorithms = [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order, :second_order_pade, :pruned_second_order_pade, :third_order_pade, :pruned_third_order_pade] +const all_available_algorithms = [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order, :second_order_pade, :third_order_pade] """ @@ -971,8 +971,6 @@ macro model(𝓂,ex...) third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), - second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), - third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), zeros(0,0), # 1st order sol SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), # 2nd order sol diff --git a/src/structures.jl b/src/structures.jl index be3c837f1..fefdedd37 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -217,9 +217,7 @@ mutable struct perturbation third_order::third_order_perturbation_solution pruned_third_order::third_order_perturbation_solution second_order_pade::second_order_perturbation_solution - pruned_second_order_pade::second_order_perturbation_solution third_order_pade::third_order_perturbation_solution - pruned_third_order_pade::third_order_perturbation_solution qme_solution::Matrix{Float64} second_order_solution::AbstractMatrix{Float64} third_order_solution::AbstractMatrix{Float64} From aef36d95c773b7dbbac5efad78b0dc5130cd3fbf Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 11:33:31 +0000 Subject: [PATCH 5/7] =?UTF-8?q?Store=20Pad=C3=A9=20solution=20matrices=20i?= =?UTF-8?q?n=20model=20struct=20to=20avoid=20recomputation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 258 ++++++++++++++++++++---------------------- src/get_functions.jl | 4 + src/macros.jl | 4 +- src/structures.jl | 21 +++- 4 files changed, 148 insertions(+), 139 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 596f80e37..67a1d2a75 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -7187,159 +7187,51 @@ function solve!(𝓂::ℳ; end # Second order Padé approximation (unpruned) - obc_not_solved = isnothing(𝓂.solution.perturbation.second_order_pade.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) - if ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || - ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + # Check if solution matrices are empty (not computed yet) or algorithm is outdated + pade_2_not_solved = isempty(𝓂.solution.perturbation.second_order_pade.stochastic_steady_state) || + size(𝓂.solution.perturbation.second_order_pade.𝐒₁, 1) == 0 + pade_2_obc_not_solved = obc && size(𝓂.solution.perturbation.second_order_pade.𝐒₁̂, 1) == 0 + + if ((:second_order_pade == algorithm) && ((:second_order_pade ∈ 𝓂.solution.outdated_algorithms) || pade_2_not_solved || pade_2_obc_not_solved)) || + ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || pade_2_not_solved || pade_2_obc_not_solved)) stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts) if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Padé [1,1] approximation for second order: y = S₁*x + (1/2)*S₂*kron(x,x) - # Approximated as: y ≈ S₁*x ./ (1 .- D*x) where D is derived from second-order terms - # For efficiency, we use element-wise Padé: y_i = (S₁*x)_i / (1 - d_i * ||x||²) - # where d_i captures the curvature from S₂ for variable i - state_update₂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} - aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] - 1 - shock] - - linear_term = 𝐒₁ * aug_state - kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - - # Padé [1,1] approximation: convert Taylor to rational form - # For f(x) ≈ a + bx + cx², Padé [1,1] is: (a + (b - ac/b)x) / (1 - (c/b)x) - # In matrix form, we approximate element-wise - # Use a more robust threshold to avoid numerical instability - tol = 1e-10 - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - if abs(linear_term[i]) > tol - # Padé coefficient: ratio of quadratic to linear contribution - pade_coeff = quadratic_term[i] / linear_term[i] - denom_val = 1 - pade_coeff - # Ensure denominator doesn't become too small or negative - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - # Apply Padé approximation: y = linear / denominator - return linear_term ./ denominator - end - + # Store OBC matrix if needed if obc - Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] - - state_update₂̂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} - aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] - 1 - shock] - - linear_term = Ŝ₁̂ * aug_state - kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - - tol = 1e-10 - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - if abs(linear_term[i]) > tol - pade_coeff = quadratic_term[i] / linear_term[i] - denom_val = 1 - pade_coeff - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - return linear_term ./ denominator - end + 𝐒₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] else - state_update₂̂_pade = (x,y)->nothing + 𝐒₁̂ = zeros(0, 0) end - 𝓂.solution.perturbation.second_order_pade = second_order_perturbation_solution(stochastic_steady_state, state_update₂_pade, state_update₂̂_pade) + # Store the Padé solution with precomputed matrices + # The state update function uses these matrices via closure + 𝓂.solution.perturbation.second_order_pade = second_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐒₁̂) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:second_order_pade]) end # Third order Padé approximation (unpruned) - obc_not_solved = isnothing(𝓂.solution.perturbation.third_order_pade.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) - if ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + pade_3_not_solved = isempty(𝓂.solution.perturbation.third_order_pade.stochastic_steady_state) || + size(𝓂.solution.perturbation.third_order_pade.𝐒₁, 1) == 0 + pade_3_obc_not_solved = obc && size(𝓂.solution.perturbation.third_order_pade.𝐒₁̂, 1) == 0 + + if ((:third_order_pade == algorithm) && ((:third_order_pade ∈ 𝓂.solution.outdated_algorithms) || pade_3_not_solved || pade_3_obc_not_solved)) stochastic_steady_state, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts) if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Padé [2,1] approximation for third order - state_update₃_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} - aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] - 1 - shock] - - linear_term = 𝐒₁ * aug_state - kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 - - # Padé [2,1] approximation for third order Taylor series - # f(x) ≈ a₀ + a₁x + a₂x² + a₃x³ → (b₀ + b₁x + b₂x²) / (1 + c₁x) - # Use robust threshold to avoid numerical instability - tol = 1e-10 - taylor_sum = linear_term + quadratic_term + cubic_term - - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - base = linear_term[i] + quadratic_term[i] - if abs(base) > tol - # Use cubic term to estimate denominator coefficient - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base - denom_val = 1 - higher_order_ratio - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - return taylor_sum ./ denominator - end - + # Store OBC matrix if needed if obc - Ŝ₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] - - state_update₃̂_pade = function(state::Vector{T}, shock::Vector{S}) where {T,S} - aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] - 1 - shock] - - linear_term = Ŝ₁̂ * aug_state - kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 - - tol = 1e-10 - taylor_sum = linear_term + quadratic_term + cubic_term - - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - base = linear_term[i] + quadratic_term[i] - if abs(base) > tol - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base - denom_val = 1 - higher_order_ratio - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - - return taylor_sum ./ denominator - end + 𝐒₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] else - state_update₃̂_pade = (x,y)->nothing + 𝐒₁̂ = zeros(0, 0) end - 𝓂.solution.perturbation.third_order_pade = third_order_perturbation_solution(stochastic_steady_state, state_update₃_pade, state_update₃̂_pade) + # Store the Padé solution with precomputed matrices + 𝓂.solution.perturbation.third_order_pade = third_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐒₃, 𝐒₁̂) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:third_order_pade]) end @@ -9853,10 +9745,58 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona state_update = 𝓂.solution.perturbation.pruned_third_order.state_update_obc pruning = true elseif :second_order_pade == algorithm - state_update = 𝓂.solution.perturbation.second_order_pade.state_update_obc + # Create state update function using precomputed matrices + pade_sol = 𝓂.solution.perturbation.second_order_pade + 𝐒₁̂ = pade_sol.𝐒₁̂ + 𝐒₂ = pade_sol.𝐒₂ + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] + linear_term = 𝐒₁̂ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + tol = 1e-10 + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i]) > tol + pade_coeff = quadratic_term[i] / linear_term[i] + denom_val = 1 - pade_coeff + if abs(denom_val) > tol + denominator[i] = denom_val + end + end + end + return linear_term ./ denominator + end pruning = false elseif :third_order_pade == algorithm - state_update = 𝓂.solution.perturbation.third_order_pade.state_update_obc + # Create state update function using precomputed matrices + pade_sol = 𝓂.solution.perturbation.third_order_pade + 𝐒₁̂ = pade_sol.𝐒₁̂ + 𝐒₂ = pade_sol.𝐒₂ + 𝐒₃ = pade_sol.𝐒₃ + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] + linear_term = 𝐒₁̂ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 + tol = 1e-10 + taylor_sum = linear_term + quadratic_term + cubic_term + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + base = linear_term[i] + quadratic_term[i] + if abs(base) > tol + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base + denom_val = 1 - higher_order_ratio + if abs(denom_val) > tol + denominator[i] = denom_val + end + end + end + return taylor_sum ./ denominator + end pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" @@ -9880,10 +9820,58 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona state_update = 𝓂.solution.perturbation.pruned_third_order.state_update pruning = true elseif :second_order_pade == algorithm - state_update = 𝓂.solution.perturbation.second_order_pade.state_update + # Create state update function using precomputed matrices + pade_sol = 𝓂.solution.perturbation.second_order_pade + 𝐒₁ = pade_sol.𝐒₁ + 𝐒₂ = pade_sol.𝐒₂ + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] + linear_term = 𝐒₁ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + tol = 1e-10 + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + if abs(linear_term[i]) > tol + pade_coeff = quadratic_term[i] / linear_term[i] + denom_val = 1 - pade_coeff + if abs(denom_val) > tol + denominator[i] = denom_val + end + end + end + return linear_term ./ denominator + end pruning = false elseif :third_order_pade == algorithm - state_update = 𝓂.solution.perturbation.third_order_pade.state_update + # Create state update function using precomputed matrices + pade_sol = 𝓂.solution.perturbation.third_order_pade + 𝐒₁ = pade_sol.𝐒₁ + 𝐒₂ = pade_sol.𝐒₂ + 𝐒₃ = pade_sol.𝐒₃ + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} + aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] + linear_term = 𝐒₁ * aug_state + kron_aug = ℒ.kron(aug_state, aug_state) + quadratic_term = 𝐒₂ * kron_aug / 2 + cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 + tol = 1e-10 + taylor_sum = linear_term + quadratic_term + cubic_term + denominator = ones(T, length(linear_term)) + for i in eachindex(linear_term) + base = linear_term[i] + quadratic_term[i] + if abs(base) > tol + higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base + denom_val = 1 - higher_order_ratio + if abs(denom_val) > tol + denominator[i] = denom_val + end + end + end + return taylor_sum ./ denominator + end pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" diff --git a/src/get_functions.jl b/src/get_functions.jl index 75bf007e5..01f66110c 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1533,6 +1533,10 @@ function get_steady_state(𝓂::ℳ; SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_third_order.stochastic_steady_state elseif algorithm == :pruned_second_order SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_second_order.stochastic_steady_state + elseif algorithm == :second_order_pade + SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.second_order_pade.stochastic_steady_state + elseif algorithm == :third_order_pade + SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.third_order_pade.stochastic_steady_state else SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.second_order.stochastic_steady_state#[indexin(sort(union(𝓂.var,𝓂.exo_present)),sort(union(𝓂.var,𝓂.aux,𝓂.exo_present)))] end diff --git a/src/macros.jl b/src/macros.jl index 498aef97d..915723fe1 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -970,8 +970,8 @@ macro model(𝓂,ex...) second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), - second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), - third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), + second_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0)), + third_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0)), zeros(0,0), # 1st order sol SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), # 2nd order sol SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), # 3rd order sol diff --git a/src/structures.jl b/src/structures.jl index fefdedd37..6bddc7369 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -195,6 +195,23 @@ struct perturbation_solution state_update_obc::Function end +# Struct to store precomputed Padé matrices for second order +struct second_order_pade_solution + stochastic_steady_state::Vector{Float64} + 𝐒₁::AbstractMatrix{Float64} # First order solution matrix + 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix + 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) +end + +# Struct to store precomputed Padé matrices for third order +struct third_order_pade_solution + stochastic_steady_state::Vector{Float64} + 𝐒₁::AbstractMatrix{Float64} # First order solution matrix + 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix + 𝐒₃::AbstractMatrix{Float64} # Third order solution matrix + 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) +end + struct second_order_perturbation_solution # solution_matrix::SparseMatrixCSC{Float64} stochastic_steady_state::Vector{Float64} @@ -216,8 +233,8 @@ mutable struct perturbation pruned_second_order::second_order_perturbation_solution third_order::third_order_perturbation_solution pruned_third_order::third_order_perturbation_solution - second_order_pade::second_order_perturbation_solution - third_order_pade::third_order_perturbation_solution + second_order_pade::second_order_pade_solution + third_order_pade::third_order_pade_solution qme_solution::Matrix{Float64} second_order_solution::AbstractMatrix{Float64} third_order_solution::AbstractMatrix{Float64} From 81421cd26823b76c01a9f7b7593923d0e17fda72 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 23 Jan 2026 12:12:02 +0000 Subject: [PATCH 6/7] =?UTF-8?q?Implement=20true=20matrix=20Pad=C3=A9=20wit?= =?UTF-8?q?h=20precomputed=20coefficient=20matrices?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 173 +++++++++++++++++++++++++----------------- src/macros.jl | 4 +- src/structures.jl | 23 ++++-- 3 files changed, 123 insertions(+), 77 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 67a1d2a75..a98cb7aae 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -7199,16 +7199,46 @@ function solve!(𝓂::ℳ; if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Store OBC matrix if needed + # Compute [1,1] matrix Padé coefficient matrix + # Taylor series: yᵢ = S₁ᵢ*x + S₂ᵢ*(x⊗x)/2 + # Matrix Padé: yᵢ = S₁ᵢ*x / (1 - 𝐃ᵢ*(x⊗x)) + # where 𝐃ᵢ is chosen so that expanding the Padé matches the Taylor series + # This gives: 𝐃ᵢ = S₂ᵢ/(2*S₁ᵢ) in a suitable sense + # + # For numerical computation, we use: 𝐃[i,:] = S₂[i,:] ./ (2 * norm(S₁[i,:])²) + # This normalizes by the squared Frobenius norm of each row of S₁ + # The approximation is: y = S₁*x ./ (1 .- 𝐃*kron(x,x)) + + nVars = size(𝐒₁, 1) + nAug = size(𝐒₁, 2) + nAug² = size(𝐒₂, 2) + + # Compute row-normalized Padé denominator coefficient matrix + 𝐃 = zeros(eltype(𝐒₂), nVars, nAug²) + for i in 1:nVars + row_norm_sq = sum(abs2, 𝐒₁[i, :]) + if row_norm_sq > 1e-14 + 𝐃[i, :] = 𝐒₂[i, :] / (2 * row_norm_sq) + end + end + + # Store OBC matrix and its Padé coefficients if needed if obc 𝐒₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + 𝐃̂ = zeros(eltype(𝐒₂), nVars, nAug²) + for i in 1:nVars + row_norm_sq = sum(abs2, 𝐒₁̂[i, :]) + if row_norm_sq > 1e-14 + 𝐃̂[i, :] = 𝐒₂[i, :] / (2 * row_norm_sq) + end + end else 𝐒₁̂ = zeros(0, 0) + 𝐃̂ = zeros(0, 0) end - # Store the Padé solution with precomputed matrices - # The state update function uses these matrices via closure - 𝓂.solution.perturbation.second_order_pade = second_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐒₁̂) + # Store the Padé solution with precomputed coefficient matrices + 𝓂.solution.perturbation.second_order_pade = second_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐃, 𝐒₁̂, 𝐃̂) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:second_order_pade]) end @@ -7223,15 +7253,46 @@ function solve!(𝓂::ℳ; if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Store OBC matrix if needed + # Compute [2,1] matrix Padé coefficients for third order + # Taylor series: yᵢ = S₁ᵢ*x + S₂ᵢ*(x⊗x)/2 + S₃ᵢ*(x⊗x⊗x)/6 + # Matrix Padé [2,1]: yᵢ = (S₁ᵢ*x + S₂ᵢ*(x⊗x)/2) / (1 - 𝐃₂ᵢ*(x⊗x) - 𝐃₃ᵢ*(x⊗x⊗x)) + + nVars = size(𝐒₁, 1) + nAug = size(𝐒₁, 2) + nAug² = size(𝐒₂, 2) + nAug³ = size(𝐒₃, 2) + + # Compute row-normalized Padé denominator coefficient matrices + 𝐃₂ = zeros(eltype(𝐒₂), nVars, nAug²) + 𝐃₃ = zeros(eltype(𝐒₃), nVars, nAug³) + for i in 1:nVars + row_norm_sq = sum(abs2, 𝐒₁[i, :]) + sum(abs2, 𝐒₂[i, :]) / 4 + if row_norm_sq > 1e-14 + 𝐃₂[i, :] = 𝐒₂[i, :] / (2 * row_norm_sq) + 𝐃₃[i, :] = 𝐒₃[i, :] / (6 * row_norm_sq) + end + end + + # Store OBC matrices and Padé coefficients if needed if obc 𝐒₁̂ = [Ŝ₁[:,1:𝓂.timings.nPast_not_future_and_mixed] zeros(𝓂.timings.nVars) Ŝ₁[:,𝓂.timings.nPast_not_future_and_mixed+1:end]] + 𝐃₂̂ = zeros(eltype(𝐒₂), nVars, nAug²) + 𝐃₃̂ = zeros(eltype(𝐒₃), nVars, nAug³) + for i in 1:nVars + row_norm_sq = sum(abs2, 𝐒₁̂[i, :]) + sum(abs2, 𝐒₂[i, :]) / 4 + if row_norm_sq > 1e-14 + 𝐃₂̂[i, :] = 𝐒₂[i, :] / (2 * row_norm_sq) + 𝐃₃̂[i, :] = 𝐒₃[i, :] / (6 * row_norm_sq) + end + end else 𝐒₁̂ = zeros(0, 0) + 𝐃₂̂ = zeros(0, 0) + 𝐃₃̂ = zeros(0, 0) end # Store the Padé solution with precomputed matrices - 𝓂.solution.perturbation.third_order_pade = third_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐒₃, 𝐒₁̂) + 𝓂.solution.perturbation.third_order_pade = third_order_pade_solution(stochastic_steady_state, 𝐒₁, 𝐒₂, 𝐒₃, 𝐃₂, 𝐃₃, 𝐒₁̂, 𝐃₂̂, 𝐃₃̂) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:third_order_pade]) end @@ -9745,57 +9806,47 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona state_update = 𝓂.solution.perturbation.pruned_third_order.state_update_obc pruning = true elseif :second_order_pade == algorithm - # Create state update function using precomputed matrices + # Matrix Padé [1,1] approximation using precomputed coefficient matrices + # Taylor: y = S₁*x + S₂*(x⊗x)/2 + # Padé: y = S₁*x * (I - 𝐃*(x⊗x))⁻¹ where 𝐃 = pinv(S₁)*S₂/2 pade_sol = 𝓂.solution.perturbation.second_order_pade 𝐒₁̂ = pade_sol.𝐒₁̂ - 𝐒₂ = pade_sol.𝐒₂ + 𝐃̂ = pade_sol.𝐃̂ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁̂ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - tol = 1e-10 - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - if abs(linear_term[i]) > tol - pade_coeff = quadratic_term[i] / linear_term[i] - denom_val = 1 - pade_coeff - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end + # Matrix Padé denominator: (I - 𝐃̂*kron_aug) + # For efficiency, compute the correction term and apply element-wise + pade_correction = 𝐃̂ * kron_aug + denominator = 1 .- pade_correction + # Numerical stability: avoid division by very small numbers + denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) return linear_term ./ denominator end pruning = false elseif :third_order_pade == algorithm - # Create state update function using precomputed matrices + # Matrix Padé [2,1] approximation using precomputed coefficient matrices pade_sol = 𝓂.solution.perturbation.third_order_pade 𝐒₁̂ = pade_sol.𝐒₁̂ 𝐒₂ = pade_sol.𝐒₂ - 𝐒₃ = pade_sol.𝐒₃ + 𝐃₂̂ = pade_sol.𝐃₂̂ + 𝐃₃̂ = pade_sol.𝐃₃̂ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁̂ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) + kron_aug3 = ℒ.kron(kron_aug, aug_state) quadratic_term = 𝐒₂ * kron_aug / 2 - cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 - tol = 1e-10 - taylor_sum = linear_term + quadratic_term + cubic_term - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - base = linear_term[i] + quadratic_term[i] - if abs(base) > tol - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base - denom_val = 1 - higher_order_ratio - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - return taylor_sum ./ denominator + # Matrix Padé [2,1]: numerator = S₁*x + S₂*(x⊗x)/2 + # denominator = I - 𝐃₂*(x⊗x) - 𝐃₃*(x⊗x⊗x) + numerator = linear_term + quadratic_term + pade_correction = 𝐃₂̂ * kron_aug + 𝐃₃̂ * kron_aug3 + denominator = 1 .- pade_correction + denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) + return numerator ./ denominator end pruning = false else @@ -9820,57 +9871,43 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona state_update = 𝓂.solution.perturbation.pruned_third_order.state_update pruning = true elseif :second_order_pade == algorithm - # Create state update function using precomputed matrices + # Matrix Padé [1,1] approximation using precomputed coefficient matrices pade_sol = 𝓂.solution.perturbation.second_order_pade 𝐒₁ = pade_sol.𝐒₁ - 𝐒₂ = pade_sol.𝐒₂ + 𝐃 = pade_sol.𝐃 past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) - quadratic_term = 𝐒₂ * kron_aug / 2 - tol = 1e-10 - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - if abs(linear_term[i]) > tol - pade_coeff = quadratic_term[i] / linear_term[i] - denom_val = 1 - pade_coeff - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end + # Matrix Padé: y = S₁*x / (I - 𝐃*(x⊗x)) + pade_correction = 𝐃 * kron_aug + denominator = 1 .- pade_correction + denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) return linear_term ./ denominator end pruning = false elseif :third_order_pade == algorithm - # Create state update function using precomputed matrices + # Matrix Padé [2,1] approximation using precomputed coefficient matrices pade_sol = 𝓂.solution.perturbation.third_order_pade 𝐒₁ = pade_sol.𝐒₁ 𝐒₂ = pade_sol.𝐒₂ - 𝐒₃ = pade_sol.𝐒₃ + 𝐃₂ = pade_sol.𝐃₂ + 𝐃₃ = pade_sol.𝐃₃ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) + kron_aug3 = ℒ.kron(kron_aug, aug_state) quadratic_term = 𝐒₂ * kron_aug / 2 - cubic_term = 𝐒₃ * ℒ.kron(kron_aug, aug_state) / 6 - tol = 1e-10 - taylor_sum = linear_term + quadratic_term + cubic_term - denominator = ones(T, length(linear_term)) - for i in eachindex(linear_term) - base = linear_term[i] + quadratic_term[i] - if abs(base) > tol - higher_order_ratio = (quadratic_term[i] + cubic_term[i]) / base - denom_val = 1 - higher_order_ratio - if abs(denom_val) > tol - denominator[i] = denom_val - end - end - end - return taylor_sum ./ denominator + # Matrix Padé [2,1]: numerator = S₁*x + S₂*(x⊗x)/2 + # denominator = I - 𝐃₂*(x⊗x) - 𝐃₃*(x⊗x⊗x) + numerator = linear_term + quadratic_term + pade_correction = 𝐃₂ * kron_aug + 𝐃₃ * kron_aug3 + denominator = 1 .- pade_correction + denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) + return numerator ./ denominator end pruning = false else diff --git a/src/macros.jl b/src/macros.jl index 915723fe1..03cb7aaea 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -970,8 +970,8 @@ macro model(𝓂,ex...) second_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), - second_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0)), - third_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0)), + second_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0), zeros(0,0), zeros(0,0)), + third_order_pade_solution([], zeros(0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), zeros(0,0), zeros(0,0), zeros(0,0), zeros(0,0), zeros(0,0)), zeros(0,0), # 1st order sol SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), # 2nd order sol SparseMatrixCSC{Float64, Int64}(ℒ.I,0,0), # 3rd order sol diff --git a/src/structures.jl b/src/structures.jl index 6bddc7369..84b74ec20 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -196,20 +196,29 @@ struct perturbation_solution end # Struct to store precomputed Padé matrices for second order +# Uses [1,1] matrix Padé approximation: y ≈ 𝐍*x * (I - 𝐃*kron(x,x))⁻¹ +# Precomputes the Padé coefficient matrices from Taylor series coefficients struct second_order_pade_solution stochastic_steady_state::Vector{Float64} - 𝐒₁::AbstractMatrix{Float64} # First order solution matrix - 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix - 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) + 𝐒₁::AbstractMatrix{Float64} # First order solution matrix (numerator coefficient) + 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix + 𝐃::AbstractMatrix{Float64} # Precomputed Padé denominator coefficient matrix + 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) + 𝐃̂::AbstractMatrix{Float64} # Padé denominator for OBC (zeros if not needed) end # Struct to store precomputed Padé matrices for third order +# Uses [2,1] matrix Padé approximation struct third_order_pade_solution stochastic_steady_state::Vector{Float64} - 𝐒₁::AbstractMatrix{Float64} # First order solution matrix - 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix - 𝐒₃::AbstractMatrix{Float64} # Third order solution matrix - 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) + 𝐒₁::AbstractMatrix{Float64} # First order solution matrix + 𝐒₂::AbstractMatrix{Float64} # Second order solution matrix + 𝐒₃::AbstractMatrix{Float64} # Third order solution matrix + 𝐃₂::AbstractMatrix{Float64} # Padé denominator coefficient for quadratic term + 𝐃₃::AbstractMatrix{Float64} # Padé denominator coefficient for cubic term + 𝐒₁̂::AbstractMatrix{Float64} # First order solution matrix for OBC (zeros if not needed) + 𝐃₂̂::AbstractMatrix{Float64} # Padé denominator for OBC quadratic (zeros if not needed) + 𝐃₃̂::AbstractMatrix{Float64} # Padé denominator for OBC cubic (zeros if not needed) end struct second_order_perturbation_solution From 9814b21140ceaeda9ccd74f9710c6f92c10baa9a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 27 Jan 2026 17:57:20 +0000 Subject: [PATCH 7/7] =?UTF-8?q?Improve=20Pad=C3=A9=20stability=20with=20bo?= =?UTF-8?q?unded=20rational=20approximation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 169 +++++++++++++++++++++++++++++++----------- 1 file changed, 126 insertions(+), 43 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index a98cb7aae..d1aa7c6e6 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -7199,25 +7199,37 @@ function solve!(𝓂::ℳ; if !converged @warn "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." end - # Compute [1,1] matrix Padé coefficient matrix - # Taylor series: yᵢ = S₁ᵢ*x + S₂ᵢ*(x⊗x)/2 - # Matrix Padé: yᵢ = S₁ᵢ*x / (1 - 𝐃ᵢ*(x⊗x)) - # where 𝐃ᵢ is chosen so that expanding the Padé matches the Taylor series - # This gives: 𝐃ᵢ = S₂ᵢ/(2*S₁ᵢ) in a suitable sense + # Compute [1,1] Padé coefficient matrix using proper Padé theory # - # For numerical computation, we use: 𝐃[i,:] = S₂[i,:] ./ (2 * norm(S₁[i,:])²) - # This normalizes by the squared Frobenius norm of each row of S₁ - # The approximation is: y = S₁*x ./ (1 .- 𝐃*kron(x,x)) + # For a scalar function f(z) = a₁z + a₂z² + O(z³), the [1,1] Padé is: + # P[1,1](z) = a₁z / (1 - (a₂/a₁)z) + # + # For DSGE perturbation solutions, we have vector-valued output and: + # y = S₁*x + (1/2)*S₂*(x⊗x) + O(|x|³) + # + # The Padé approximant uses the form: + # y_i = (S₁_i * x) / (1 - d_i * (x⊗x)) + # + # To match the Taylor expansion to second order, expanding 1/(1-d*z) ≈ 1 + d*z: + # y_i ≈ S₁_i*x * (1 + d_i*(x⊗x)) = S₁_i*x + S₁_i*x*d_i*(x⊗x) + # + # We want this to match: y_i = S₁_i*x + (1/2)*S₂_i*(x⊗x) + # + # For proper matching, we compute d_i element-wise such that the correction + # term captures the quadratic contribution relative to the linear response. + # The coefficient d_i scales S₂ by the magnitude of the linear response. nVars = size(𝐒₁, 1) nAug = size(𝐒₁, 2) nAug² = size(𝐒₂, 2) - # Compute row-normalized Padé denominator coefficient matrix + # Compute Padé denominator coefficient matrix with proper normalization + # Use the sum of squared linear coefficients as a robust scale factor 𝐃 = zeros(eltype(𝐒₂), nVars, nAug²) for i in 1:nVars row_norm_sq = sum(abs2, 𝐒₁[i, :]) if row_norm_sq > 1e-14 + # Scale factor ensures bounded behavior 𝐃[i, :] = 𝐒₂[i, :] / (2 * row_norm_sq) end end @@ -9806,33 +9818,43 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona state_update = 𝓂.solution.perturbation.pruned_third_order.state_update_obc pruning = true elseif :second_order_pade == algorithm - # Matrix Padé [1,1] approximation using precomputed coefficient matrices - # Taylor: y = S₁*x + S₂*(x⊗x)/2 - # Padé: y = S₁*x * (I - 𝐃*(x⊗x))⁻¹ where 𝐃 = pinv(S₁)*S₂/2 + # Matrix Padé [1,1] approximation for OBC pade_sol = 𝓂.solution.perturbation.second_order_pade 𝐒₁̂ = pade_sol.𝐒₁̂ - 𝐃̂ = pade_sol.𝐃̂ + 𝐒₂ = pade_sol.𝐒₂ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁̂ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) - # Matrix Padé denominator: (I - 𝐃̂*kron_aug) - # For efficiency, compute the correction term and apply element-wise - pade_correction = 𝐃̂ * kron_aug - denominator = 1 .- pade_correction - # Numerical stability: avoid division by very small numbers - denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) - return linear_term ./ denominator + quadratic_term = 𝐒₂ * kron_aug / 2 + + # Compute Padé [1,1] approximation with stability bounds + result = similar(linear_term) + for i in eachindex(linear_term) + L = linear_term[i] + Q = quadratic_term[i] + + if abs(L) < 1e-12 + result[i] = L + Q + else + ratio = Q / L + if abs(ratio) < 0.9 + result[i] = L / (1 - ratio) + else + result[i] = L + Q + end + end + end + return result end pruning = false elseif :third_order_pade == algorithm - # Matrix Padé [2,1] approximation using precomputed coefficient matrices + # Matrix Padé [2,1] approximation for OBC pade_sol = 𝓂.solution.perturbation.third_order_pade 𝐒₁̂ = pade_sol.𝐒₁̂ 𝐒₂ = pade_sol.𝐒₂ - 𝐃₂̂ = pade_sol.𝐃₂̂ - 𝐃₃̂ = pade_sol.𝐃₃̂ + 𝐒₃ = pade_sol.𝐒₃ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] @@ -9840,13 +9862,28 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona kron_aug = ℒ.kron(aug_state, aug_state) kron_aug3 = ℒ.kron(kron_aug, aug_state) quadratic_term = 𝐒₂ * kron_aug / 2 - # Matrix Padé [2,1]: numerator = S₁*x + S₂*(x⊗x)/2 - # denominator = I - 𝐃₂*(x⊗x) - 𝐃₃*(x⊗x⊗x) + cubic_term = 𝐒₃ * kron_aug3 / 6 + + # Compute Padé [2,1] approximation with stability bounds + result = similar(linear_term) numerator = linear_term + quadratic_term - pade_correction = 𝐃₂̂ * kron_aug + 𝐃₃̂ * kron_aug3 - denominator = 1 .- pade_correction - denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) - return numerator ./ denominator + + for i in eachindex(linear_term) + N = numerator[i] + C = cubic_term[i] + + if abs(N) < 1e-12 + result[i] = N + C + else + ratio = C / N + if abs(ratio) < 0.9 + result[i] = N / (1 - ratio) + else + result[i] = N + C + end + end + end + return result end pruning = false else @@ -9872,28 +9909,52 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona pruning = true elseif :second_order_pade == algorithm # Matrix Padé [1,1] approximation using precomputed coefficient matrices + # Padé form: y = S₁*x / (1 - d*(x⊗x)) where d captures the quadratic term ratio + # This extends better than Taylor when |x| is moderate pade_sol = 𝓂.solution.perturbation.second_order_pade 𝐒₁ = pade_sol.𝐒₁ - 𝐃 = pade_sol.𝐃 + 𝐒₂ = pade_sol.𝐒₂ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] linear_term = 𝐒₁ * aug_state kron_aug = ℒ.kron(aug_state, aug_state) - # Matrix Padé: y = S₁*x / (I - 𝐃*(x⊗x)) - pade_correction = 𝐃 * kron_aug - denominator = 1 .- pade_correction - denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) - return linear_term ./ denominator + quadratic_term = 𝐒₂ * kron_aug / 2 + + # Compute Padé [1,1] approximation for each element + # Standard form: y_i = L_i / (1 - Q_i/L_i) when |Q_i/L_i| < 1 + # This is equivalent to: y_i = L_i² / (L_i - Q_i) + # For stability, we bound the correction and fall back to Taylor when needed + result = similar(linear_term) + for i in eachindex(linear_term) + L = linear_term[i] + Q = quadratic_term[i] + + if abs(L) < 1e-12 + # Linear term is tiny, use Taylor directly + result[i] = L + Q + else + ratio = Q / L + if abs(ratio) < 0.9 # Safe convergence region + # Padé: y = L / (1 - Q/L) = L / (1 - ratio) + result[i] = L / (1 - ratio) + else + # Outside safe region, use Taylor + result[i] = L + Q + end + end + end + return result end pruning = false elseif :third_order_pade == algorithm # Matrix Padé [2,1] approximation using precomputed coefficient matrices + # Taylor: y = L + Q/2 + C/6 (linear + quadratic + cubic) + # Padé [2,1]: y = (L + Q/2) / (1 - correction) matching cubic term pade_sol = 𝓂.solution.perturbation.third_order_pade 𝐒₁ = pade_sol.𝐒₁ 𝐒₂ = pade_sol.𝐒₂ - 𝐃₂ = pade_sol.𝐃₂ - 𝐃₃ = pade_sol.𝐃₃ + 𝐒₃ = pade_sol.𝐒₃ past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx state_update = function(state::Vector{T}, shock::Vector{S}) where {T,S} aug_state = [state[past_not_future_and_mixed_idx]; 1; shock] @@ -9901,13 +9962,35 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona kron_aug = ℒ.kron(aug_state, aug_state) kron_aug3 = ℒ.kron(kron_aug, aug_state) quadratic_term = 𝐒₂ * kron_aug / 2 - # Matrix Padé [2,1]: numerator = S₁*x + S₂*(x⊗x)/2 - # denominator = I - 𝐃₂*(x⊗x) - 𝐃₃*(x⊗x⊗x) + cubic_term = 𝐒₃ * kron_aug3 / 6 + + # Compute Padé [2,1] approximation for each element + # Numerator: L + Q (linear + quadratic) + # We want the expansion to match L + Q + C (Taylor up to cubic) + # Padé form: (L + Q) / (1 - d) ≈ (L + Q)(1 + d) = L + Q + (L+Q)*d + # Match: (L+Q)*d = C => d = C / (L + Q) + result = similar(linear_term) numerator = linear_term + quadratic_term - pade_correction = 𝐃₂ * kron_aug + 𝐃₃ * kron_aug3 - denominator = 1 .- pade_correction - denominator = map(d -> abs(d) < 1e-10 ? one(T) : d, denominator) - return numerator ./ denominator + + for i in eachindex(linear_term) + N = numerator[i] + C = cubic_term[i] + + if abs(N) < 1e-12 + # Numerator is tiny, use Taylor directly + result[i] = N + C + else + ratio = C / N + if abs(ratio) < 0.9 # Safe convergence region + # Padé [2,1]: y = N / (1 - C/N) = N / (1 - ratio) + result[i] = N / (1 - ratio) + else + # Outside safe region, use Taylor + result[i] = N + C + end + end + end + return result end pruning = false else