From c7fbfbf3220a913e11ad768f8681548d2abfb38a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 31 Dec 2025 16:31:03 +0000 Subject: [PATCH 01/40] Initial plan From bb3f131de4ae0f918db11d5eb807d529a6fd3984 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 31 Dec 2025 17:16:25 +0000 Subject: [PATCH 02/40] Fix convert_to_ss_equation return type and add simplify methods for non-Expr types Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index aa32cae2e..200ca0722 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -2936,6 +2936,11 @@ end Max = max Min = min +# Simplify methods for non-Expr types (pass through) +simplify(ex::Symbol)::Symbol = ex +simplify(ex::Int)::Int = ex +simplify(ex::Float64)::Float64 = ex + function simplify(ex::Expr)::Union{Expr,Symbol,Int} ex_ss = convert_to_ss_equation(ex) @@ -2953,7 +2958,7 @@ function simplify(ex::Expr)::Union{Expr,Symbol,Int} x, parsed) end -function convert_to_ss_equation(eq::Expr)::Expr +function convert_to_ss_equation(eq::Expr)::Union{Expr, Symbol, Int} postwalk(x -> x isa Expr ? x.head == :(=) ? From 5e793c26bd784f020d6f5d599ed130d28b82e6d8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 31 Dec 2025 17:24:20 +0000 Subject: [PATCH 03/40] Add type checks in macros.jl for backward looking models with exp/log/etc functions Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/macros.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/macros.jl b/src/macros.jl index c7c0f8ed1..37a775690 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -333,14 +333,14 @@ macro model(𝓂,ex...) bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 1e12)) : (eps(), 1e12) x end : - x.args[2].head == :ref ? + x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 1e12)) : (eps(), 1e12) x end : x : - x.args[2].head == :call ? # nonnegative expressions + x.args[2] isa Expr && x.args[2].head == :call ? # nonnegative expressions begin if precompile replacement = x.args[2] @@ -381,14 +381,14 @@ macro model(𝓂,ex...) bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 1e12)) : (eps(), 1e12) x end : - x.args[2].head == :ref ? + x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 1e12)) : (eps(), 1e12) x end : x : - x.args[2].head == :call ? # nonnegative expressions + x.args[2] isa Expr && x.args[2].head == :call ? # nonnegative expressions begin if precompile replacement = x.args[2] @@ -425,14 +425,14 @@ macro model(𝓂,ex...) bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 1-eps())) : (eps(), 1-eps()) x end : - x.args[2].head == :ref ? + x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 1-eps())) : (eps(), 1-eps()) x end : x : - x.args[2].head == :call ? # nonnegative expressions + x.args[2] isa Expr && x.args[2].head == :call ? # nonnegative expressions begin if precompile replacement = x.args[2] @@ -469,14 +469,14 @@ macro model(𝓂,ex...) bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], -1e12), min(bounds[x.args[2]][2], 600)) : (-1e12, 600) x end : - x.args[2].head == :ref ? + x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # have exp terms bound so they dont go to Inf begin bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], -1e12), min(bounds[x.args[2].args[1]][2], 600)) : (-1e12, 600) x end : x : - x.args[2].head == :call ? # nonnegative expressions + x.args[2] isa Expr && x.args[2].head == :call ? # nonnegative expressions begin if precompile replacement = x.args[2] @@ -513,14 +513,14 @@ macro model(𝓂,ex...) bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 2-eps())) : (eps(), 2-eps()) x end : - x.args[2].head == :ref ? + x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin - bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 2-eps())) : (eps(), 2-eps()) + bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]] ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 2-eps())) : (eps(), 2-eps()) x end : x : - x.args[2].head == :call ? # nonnegative expressions + x.args[2] isa Expr && x.args[2].head == :call ? # nonnegative expressions begin if precompile replacement = x.args[2] From 1efac1359c403ff2f89135abaf7d6c5cca5da7c7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 31 Dec 2025 17:42:42 +0000 Subject: [PATCH 04/40] =?UTF-8?q?Fix=20syntax=20error=20in=20macros.jl:=20?= =?UTF-8?q?])=20=E2=86=92=20)=20in=20haskey=20call?= 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/macros.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/macros.jl b/src/macros.jl index 37a775690..c5b1f5f8e 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -516,7 +516,7 @@ macro model(𝓂,ex...) x.args[2] isa Expr && x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin - bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]] ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 2-eps())) : (eps(), 2-eps()) + bounds[x.args[2].args[1]] = haskey(bounds, x.args[2].args[1]) ? (max(bounds[x.args[2].args[1]][1], eps()), min(bounds[x.args[2].args[1]][2], 2-eps())) : (eps(), 2-eps()) x end : x : From db479ab771f9fe45fb379552b270df0bc2d922bc Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 31 Dec 2025 18:17:39 +0000 Subject: [PATCH 05/40] Add test file for backward looking models (VAR and Solow growth) Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- test/runtests.jl | 2 + test/test_backward_looking_models.jl | 93 ++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) create mode 100644 test/test_backward_looking_models.jl diff --git a/test/runtests.jl b/test/runtests.jl index 23667f45e..9e742a6f7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -793,6 +793,8 @@ if test_set == "basic" plots = false # test_higher_order = false + include("test_backward_looking_models.jl") + @testset verbose = true "Provide parameters later" begin include("models/Backus_Kehoe_Kydland_1992.jl") diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl new file mode 100644 index 000000000..46195dc3b --- /dev/null +++ b/test/test_backward_looking_models.jl @@ -0,0 +1,93 @@ +# Tests for purely backward looking models (nFuture_not_past_and_mixed = 0) +# These models have no forward-looking variables and include VARs and nonlinear growth models + +@testset verbose = true "Backward looking models" begin + + @testset "Linear backward looking model (VAR)" begin + # Two-variable VAR model + @model VAR2 begin + y[0] = a11 * y[-1] + a12 * x[-1] + sigma_y * eps_y[x] + x[0] = a21 * y[-1] + a22 * x[-1] + sigma_x * eps_x[x] + end + + @parameters VAR2 begin + a11 = 0.5 + a12 = 0.3 + a21 = 0.2 + a22 = 0.4 + sigma_y = 0.1 + sigma_x = 0.1 + end + + # Verify it's a backward looking model + @test VAR2.timings.nFuture_not_past_and_mixed == 0 + + # Test steady state (should be zero for this linear model) + ss = get_steady_state(VAR2) + @test isapprox(ss(:y, :Steady_state), 0.0, atol=1e-10) + @test isapprox(ss(:x, :Steady_state), 0.0, atol=1e-10) + + # Test solution + sol = get_solution(VAR2) + @test !isnothing(sol) + + # Test IRFs + irf = get_irf(VAR2) + @test size(irf, 1) == 2 # 2 variables + @test size(irf, 3) == 2 # 2 shocks + + # Test simulation + sim = simulate(VAR2) + @test size(sim, 1) == 2 # 2 variables + + VAR2 = nothing + end + + @testset "Nonlinear backward looking model (Solow growth)" begin + # Solow growth model - purely backward looking nonlinear model + @model SolowGrowth begin + k[0] = (1 - delta) * k[-1] + s * k[-1]^alpha * exp(z[-1]) + z[0] = rho * z[-1] + sigma * eps_z[x] + end + + @parameters SolowGrowth begin + alpha = 0.3 + delta = 0.1 + s = 0.2 + rho = 0.9 + sigma = 0.01 + end + + # Verify it's a backward looking model + @test SolowGrowth.timings.nFuture_not_past_and_mixed == 0 + + # Test steady state exists and is positive for capital + ss = get_steady_state(SolowGrowth) + k_ss = ss(:k, :Steady_state) + z_ss = ss(:z, :Steady_state) + + @test k_ss > 0 # Capital should be positive at steady state + @test isapprox(z_ss, 0.0, atol=1e-10) # z steady state should be 0 + + # Verify steady state: k_ss = (s/delta)^(1/(1-alpha)) + k_ss_analytical = (0.2/0.1)^(1/(1-0.3)) + @test isapprox(k_ss, k_ss_analytical, rtol=1e-6) + + # Test solution + sol = get_solution(SolowGrowth) + @test !isnothing(sol) + + # Test IRFs + irf = get_irf(SolowGrowth) + @test size(irf, 1) == 2 # 2 variables (k, z) + @test size(irf, 3) == 1 # 1 shock + + # Test simulation + sim = simulate(SolowGrowth) + @test size(sim, 1) == 2 # 2 variables + + SolowGrowth = nothing + end + + GC.gc() +end From 884b1fa7e10e415337fe764693c127b7410c01b8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 20:49:50 +0000 Subject: [PATCH 06/40] Add newton algorithm for backward looking models with residual and jacobian functions Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 279 ++++++++++++++++++++++++++- src/common_docstrings.jl | 2 +- src/get_functions.jl | 6 +- src/macros.jl | 3 +- src/structures.jl | 3 + test/test_backward_looking_models.jl | 27 ++- 6 files changed, 308 insertions(+), 12 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 200ca0722..2e197bdd4 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -2906,7 +2906,7 @@ function get_relevant_steady_states(𝓂::ℳ, end relevant_SS = get_steady_state(𝓂, algorithm = algorithm, - stochastic = algorithm != :first_order, + stochastic = algorithm ∉ [:first_order, :newton], return_variables_only = true, derivatives = false, verbose = opts.verbose, @@ -6795,7 +6795,8 @@ 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))) || + ((:newton == algorithm) && ((:newton ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) # @timeit_debug timer "Solve for NSSS (if necessary)" begin @@ -6858,6 +6859,12 @@ function solve!(𝓂::ℳ; 𝓂.solution.non_stochastic_steady_state = SS_and_pars 𝓂.solution.outdated_NSSS = solution_error > opts.tol.NSSS_acceptance_tol + + # Generate newton simulation functions for backward looking models + if 𝓂.timings.nFuture_not_past_and_mixed == 0 + write_newton_simulation_functions!(𝓂) + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:newton]) + end end obc_not_solved = isnothing(𝓂.solution.perturbation.second_order.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) @@ -7997,6 +8004,131 @@ function write_auxiliary_indices!(𝓂::ℳ) return nothing end + +""" + write_newton_simulation_functions!(𝓂::ℳ) + +Generate residual and jacobian functions for newton-based simulation of backward looking models. +The residual function evaluates F(x_t, x_{t-1}, e_t) = 0 +The jacobian function computes ∂F/∂x_t (derivatives w.r.t. present values) + +These functions are only generated for backward looking models (nFuture_not_past_and_mixed = 0). +""" +function write_newton_simulation_functions!(𝓂::ℳ; + cse = true, + skipzeros = true) + # Only generate for backward looking models + if 𝓂.timings.nFuture_not_past_and_mixed > 0 + return nothing + end + + # Get all variable lists from dynamic equations + dyn_var_present_list = collect(reduce(union, 𝓂.dyn_present_list)) + dyn_var_past_list = collect(reduce(union, 𝓂.dyn_past_list)) + dyn_exo_list = collect(reduce(union, 𝓂.dyn_exo_list)) + dyn_ss_list = Symbol.(string.(collect(reduce(union, 𝓂.dyn_ss_list))) .* "₍ₛₛ₎") + + present = map(x -> Symbol(replace(string(x), r"₍₀₎" => "")), string.(dyn_var_present_list)) + past = map(x -> Symbol(replace(string(x), r"₍₋₁₎" => "")), string.(dyn_var_past_list)) + exo = map(x -> Symbol(replace(string(x), r"₍ₓ₎" => "")), string.(dyn_exo_list)) + stst = map(x -> Symbol(replace(string(x), r"₍ₛₛ₎" => "")), string.(dyn_ss_list)) + + # Sort variables + vars_raw = vcat(dyn_var_present_list[indexin(sort(present), present)], + dyn_var_past_list[indexin(sort(past), past)], + dyn_exo_list[indexin(sort(exo), exo)]) + + pars_ext = vcat(𝓂.parameters, 𝓂.calibration_equations_parameters) + parameters_and_SS = vcat(pars_ext, dyn_ss_list[indexin(sort(stst), stst)]) + + np = length(parameters_and_SS) + nv = length(vars_raw) + nps = length(𝓂.parameters) + + # Number of present, past, and shock variables + n_present = length(dyn_var_present_list) + n_past = length(dyn_var_past_list) + n_exo = length(dyn_exo_list) + nVars = 𝓂.timings.nVars + + # Create symbolic variables: 𝔙 = [present; past; shocks] + Symbolics.@variables 𝔓[1:np] 𝔙[1:nv] + + parameter_dict = Dict{Symbol, Symbol}() + back_to_array_dict = Dict{Symbolics.Num, Symbolics.Num}() + + for (i, v) in enumerate(parameters_and_SS) + push!(parameter_dict, v => :($(Symbol("𝔓_$i")))) + push!(back_to_array_dict, Symbolics.parse_expr_to_symbolic(:($(Symbol("𝔓_$i"))), @__MODULE__) => 𝔓[i]) + end + + for (i, v) in enumerate(vars_raw) + push!(parameter_dict, v => :($(Symbol("𝔙_$i")))) + push!(back_to_array_dict, Symbolics.parse_expr_to_symbolic(:($(Symbol("𝔙_$i"))), @__MODULE__) => 𝔙[i]) + end + + # Handle calibration equations + calib_vars = Symbol[] + calib_expr = [] + for v in 𝓂.calibration_equations_no_var + push!(calib_vars, v.args[1]) + push!(calib_expr, v.args[2]) + end + + calib_replacements = Dict{Symbol, Any}() + for (i, x) in enumerate(calib_vars) + replacement = Dict(x => calib_expr[i]) + for ii in i+1:length(calib_vars) + calib_expr[ii] = replace_symbols(calib_expr[ii], replacement) + end + push!(calib_replacements, x => calib_expr[i]) + end + + # Parse dynamic equations + dyn_equations = 𝓂.dyn_equations |> + x -> replace_symbols.(x, Ref(calib_replacements)) |> + x -> replace_symbols.(x, Ref(parameter_dict)) |> + x -> Symbolics.parse_expr_to_symbolic.(x, Ref(@__MODULE__)) |> + x -> Symbolics.substitute.(x, Ref(back_to_array_dict)) + + # Convert to vector for build_function + dyn_equations_vec = Symbolics.Num.(dyn_equations) + + # Build residual function: F(present, past, shocks; parameters) = 0 + _, residual_func = Symbolics.build_function(dyn_equations_vec, 𝔓, 𝔙, + cse = cse, + expression_module = @__MODULE__, + expression = Val(false))::Tuple{<:Function, <:Function} + + # Compute jacobian w.r.t. present variables only (first n_present elements of 𝔙) + present_vars = 𝔙[1:n_present] + ∇_present = Symbolics.sparsejacobian(dyn_equations_vec, collect(present_vars)) + + # Convert to matrix form + ∇_present_mat = convert(Matrix, ∇_present) + + # Build jacobian function + _, jacobian_func = Symbolics.build_function(∇_present_mat, 𝔓, 𝔙, + cse = cse, + skipzeros = skipzeros, + expression_module = @__MODULE__, + expression = Val(false))::Tuple{<:Function, <:Function} + + # Create buffers + residual_buffer = zeros(Float64, length(dyn_equations)) + jacobian_buffer = zeros(Float64, size(∇_present_mat)) + + # Store in model + 𝓂.newton_simulation_functions = ( + residual_func = residual_func, + jacobian_func = jacobian_func, + residual_buffer = residual_buffer, + jacobian_buffer = jacobian_buffer + ) + + return nothing +end + write_parameters_input!(𝓂::ℳ, parameters::Nothing; verbose::Bool = true) = return parameters write_parameters_input!(𝓂::ℳ, parameters::Pair{Symbol,Float64}; verbose::Bool = true) = write_parameters_input!(𝓂::ℳ, OrderedDict(parameters), verbose = verbose) write_parameters_input!(𝓂::ℳ, parameters::Pair{S,Float64}; verbose::Bool = true) where S <: AbstractString = write_parameters_input!(𝓂::ℳ, OrderedDict{Symbol,Float64}(parameters[1] |> Meta.parse |> replace_indices => parameters[2]), verbose = verbose) @@ -9506,6 +9638,11 @@ 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 :newton == algorithm + @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "Newton algorithm is only available for backward looking models (nFuture_not_past_and_mixed = 0)." + # For OBC, fall back to first order for now + state_update = 𝓂.solution.perturbation.first_order.state_update_obc + pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -9527,6 +9664,11 @@ 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 :newton == algorithm + @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "Newton algorithm is only available for backward looking models (nFuture_not_past_and_mixed = 0)." + # Create newton-based state update function + state_update = create_newton_state_update(𝓂) + pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" state_update = (x,y)->nothing @@ -9537,6 +9679,139 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona return state_update, pruning end + +""" + create_newton_state_update(𝓂::ℳ) + +Create a state update function that uses Newton's method to solve for present values +given past values and shocks. Only for backward looking models. +""" +function create_newton_state_update(𝓂::ℳ) + # Get the newton simulation functions + residual_func = 𝓂.newton_simulation_functions.residual_func + jacobian_func = 𝓂.newton_simulation_functions.jacobian_func + + # Get steady state and parameters + SS_and_pars = 𝓂.solution.non_stochastic_steady_state + parameters = 𝓂.parameter_values + + # Get variable indices + dyn_var_present_list = collect(reduce(union, 𝓂.dyn_present_list)) + dyn_var_past_list = collect(reduce(union, 𝓂.dyn_past_list)) + dyn_exo_list = collect(reduce(union, 𝓂.dyn_exo_list)) + + present = map(x -> Symbol(replace(string(x), r"₍₀₎" => "")), string.(dyn_var_present_list)) + past = map(x -> Symbol(replace(string(x), r"₍₋₁₎" => "")), string.(dyn_var_past_list)) + exo = map(x -> Symbol(replace(string(x), r"₍ₓ₎" => "")), string.(dyn_exo_list)) + + n_present = length(present) + n_past = length(past) + n_exo = length(exo) + nVars = 𝓂.timings.nVars + + # Get SS and parameter info + dyn_ss_list = Symbol.(string.(collect(reduce(union, 𝓂.dyn_ss_list))) .* "₍ₛₛ₎") + stst = map(x -> Symbol(replace(string(x), r"₍ₛₛ₎" => "")), string.(dyn_ss_list)) + + pars_ext = vcat(𝓂.parameters, 𝓂.calibration_equations_parameters) + parameters_and_SS_syms = vcat(pars_ext, dyn_ss_list[indexin(sort(stst), stst)]) + + SS_and_pars_names = vcat(Symbol.(string.(sort(union(𝓂.var, 𝓂.exo_past, 𝓂.exo_future)))), 𝓂.calibration_equations_parameters) + + # Get indices for mapping + present_sorted = sort(present) + past_sorted = sort(past) + exo_sorted = sort(exo) + + present_idx_in_SS = indexin(present_sorted, SS_and_pars_names) + past_idx_in_SS = indexin(past_sorted, SS_and_pars_names) + + # Get indices for past_not_future_and_mixed (used in state vector) + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + + # Create the state update function + function state_update_newton(state::Vector{T}, shock::Vector{S}) where {T, S} + # Get steady state values for parameters and SS variables + # Parameters are in the first part, SS values come after + pars_values = 𝓂.parameter_values + ss_values = SS_and_pars[indexin(sort(stst), SS_and_pars_names)] + params_and_ss = vcat(pars_values, ss_values) + + # Initial guess: steady state values for present variables + present_guess = SS_and_pars[present_idx_in_SS] + + # Get past values from state + past_values = state[past_not_future_and_mixed_idx][indexin(past_sorted, 𝓂.timings.past_not_future_and_mixed[indexin(past_sorted, 𝓂.timings.past_not_future_and_mixed)])] + + # Reorder past values to match the sorted order in the dynamic equations + # past_values should be ordered same as past_sorted + past_from_state = zeros(n_past) + for (i, p) in enumerate(past_sorted) + idx_in_state = findfirst(x -> x == p, 𝓂.timings.past_not_future_and_mixed) + if idx_in_state !== nothing + past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + end + end + + # Get shock values - reorder to match sorted order + shock_values = zeros(n_exo) + for (i, e) in enumerate(exo_sorted) + idx_in_shock = findfirst(x -> x == e, 𝓂.timings.exo) + if idx_in_shock !== nothing + shock_values[i] = shock[idx_in_shock] + end + end + + # Combine: 𝔙 = [present; past; shocks] + vars = vcat(present_guess, past_from_state, shock_values) + + # Newton iterations + max_iter = 50 + tol = 1e-10 + + residual = zeros(length(present_guess)) + jacobian = zeros(length(present_guess), n_present) + + for iter in 1:max_iter + # Update vars with current present guess + vars[1:n_present] = present_guess + + # Evaluate residual and jacobian + residual_func(residual, params_and_ss, vars) + jacobian_func(jacobian, params_and_ss, vars) + + # Check convergence + if ℒ.norm(residual) < tol + break + end + + # Newton step: present_new = present_old - J^{-1} * F + try + Δ = jacobian \ residual + present_guess = present_guess - Δ + catch + # If jacobian is singular, break + break + end + end + + # Construct full state vector in the correct order + result = copy(SS_and_pars[1:nVars]) # Start with SS values + + # Fill in the present values we just computed + for (i, p) in enumerate(present_sorted) + idx = findfirst(x -> x == p, 𝓂.timings.var) + if idx !== nothing + result[idx] = present_guess[i] + end + end + + return result + end + + return state_update_newton +end + @stable default_mode = "disable" begin function find_variables_to_exclude(𝓂::ℳ, observables::Vector{Symbol}) diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl index 9238fc4e7..af6be3e33 100644 --- a/src/common_docstrings.jl +++ b/src/common_docstrings.jl @@ -11,7 +11,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`, `:newton` (for backward looking models only - uses Newton's method to solve the nonlinear equations directly)" 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 cc8eea449..9e44e01bf 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1447,12 +1447,12 @@ function get_steady_state(𝓂::ℳ; sylvester_algorithm³ = (isa(sylvester_algorithm, Symbol) || length(sylvester_algorithm) < 2) ? :bicgstab : sylvester_algorithm[2]) if stochastic - if algorithm == :first_order + if algorithm ∈ [:first_order, :newton] @info "Stochastic steady state requested but algorithm is $algorithm. Setting `algorithm = :second_order`." maxlog = DEFAULT_MAXLOG algorithm = :second_order end else - if algorithm != :first_order + if algorithm ∉ [:first_order, :newton] @info "Non-stochastic steady state requested but algorithm is $algorithm. Setting `stochastic = true`." maxlog = DEFAULT_MAXLOG stochastic = true end @@ -1741,7 +1741,7 @@ function get_solution(𝓂::ℳ; silent = silent, algorithm = algorithm) - if algorithm == :first_order + if algorithm ∈ [:first_order, :newton] solution_matrix = 𝓂.solution.perturbation.first_order.solution_matrix end diff --git a/src/macros.jl b/src/macros.jl index c5b1f5f8e..ce08c0132 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, :newton] """ @@ -913,6 +913,7 @@ macro model(𝓂,ex...) (zeros(0,0), x->x), # third_order_derivatives (zeros(0,0), x->x), # third_order_derivatives_parameters (zeros(0,0), x->x), # third_order_derivatives_SS_and_pars + (residual_func = (x,y,z)->nothing, jacobian_func = (x,y,z)->nothing, residual_buffer = Float64[], jacobian_buffer = zeros(0,0)), # newton_simulation_functions # (x->x, SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0), 𝒟.prepare_jacobian(x->x, 𝒟.AutoForwardDiff(), [0]), SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0)), # third_order_derivatives # ([], SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0)), # model_jacobian # ([], Int[], zeros(1,1)), # model_jacobian diff --git a/src/structures.jl b/src/structures.jl index f918b6a37..b887ccfaf 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -426,6 +426,9 @@ mutable struct ℳ third_order_derivatives_parameters::Tuple{AbstractMatrix{<: Real},Function} third_order_derivatives_SS_and_pars::Tuple{AbstractMatrix{<: Real},Function} + # Newton simulation functions for backward looking models + newton_simulation_functions::NamedTuple{(:residual_func, :jacobian_func, :residual_buffer, :jacobian_buffer), Tuple{Function, Function, Vector{Float64}, Matrix{Float64}}} + # model_jacobian::Tuple{Vector{Function}, SparseMatrixCSC{Float64}} # model_jacobian::Tuple{Vector{Function}, Vector{Int}, Matrix{<: Real}} # # model_jacobian_parameters::Function diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 46195dc3b..e30a7e1c1 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -31,15 +31,23 @@ sol = get_solution(VAR2) @test !isnothing(sol) - # Test IRFs - irf = get_irf(VAR2) - @test size(irf, 1) == 2 # 2 variables - @test size(irf, 3) == 2 # 2 shocks + # Test IRFs with first_order algorithm + irf_first = get_irf(VAR2, algorithm = :first_order) + @test size(irf_first, 1) == 2 # 2 variables + @test size(irf_first, 3) == 2 # 2 shocks # Test simulation sim = simulate(VAR2) @test size(sim, 1) == 2 # 2 variables + # Test newton algorithm - should produce same results as first_order for linear model + irf_newton = get_irf(VAR2, algorithm = :newton) + @test size(irf_newton, 1) == 2 # 2 variables + @test size(irf_newton, 3) == 2 # 2 shocks + + # For linear models, first_order and newton should give approximately the same results + @test isapprox(collect(irf_first), collect(irf_newton), rtol=1e-6) + VAR2 = nothing end @@ -77,7 +85,7 @@ sol = get_solution(SolowGrowth) @test !isnothing(sol) - # Test IRFs + # Test IRFs with first_order algorithm irf = get_irf(SolowGrowth) @test size(irf, 1) == 2 # 2 variables (k, z) @test size(irf, 3) == 1 # 1 shock @@ -86,6 +94,15 @@ sim = simulate(SolowGrowth) @test size(sim, 1) == 2 # 2 variables + # Test newton algorithm for nonlinear backward looking model + irf_newton = get_irf(SolowGrowth, algorithm = :newton) + @test size(irf_newton, 1) == 2 # 2 variables (k, z) + @test size(irf_newton, 3) == 1 # 1 shock + + # Test simulation with newton + sim_newton = simulate(SolowGrowth, algorithm = :newton) + @test size(sim_newton, 1) == 2 # 2 variables + SolowGrowth = nothing end From 1359ddbc51294f1df070a488a584439ba25f74ad Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 21:48:31 +0000 Subject: [PATCH 07/40] Fix newton state_update to add SS before evaluation and subtract after Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 2e197bdd4..967c8dd4b 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -9729,6 +9729,11 @@ function create_newton_state_update(𝓂::ℳ) # Get indices for past_not_future_and_mixed (used in state vector) past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + # Get steady state values for present and past variables + SS_present = SS_and_pars[present_idx_in_SS] + SS_past = SS_and_pars[past_idx_in_SS] + SS_vars = SS_and_pars[1:nVars] + # Create the state update function function state_update_newton(state::Vector{T}, shock::Vector{S}) where {T, S} # Get steady state values for parameters and SS variables @@ -9737,19 +9742,17 @@ function create_newton_state_update(𝓂::ℳ) ss_values = SS_and_pars[indexin(sort(stst), SS_and_pars_names)] params_and_ss = vcat(pars_values, ss_values) - # Initial guess: steady state values for present variables - present_guess = SS_and_pars[present_idx_in_SS] - - # Get past values from state - past_values = state[past_not_future_and_mixed_idx][indexin(past_sorted, 𝓂.timings.past_not_future_and_mixed[indexin(past_sorted, 𝓂.timings.past_not_future_and_mixed)])] + # Initial guess: steady state values for present variables (in levels) + present_guess = copy(SS_present) - # Reorder past values to match the sorted order in the dynamic equations - # past_values should be ordered same as past_sorted + # Get past values from state (which are in deviations from SS) + # Add steady state to convert to levels past_from_state = zeros(n_past) for (i, p) in enumerate(past_sorted) idx_in_state = findfirst(x -> x == p, 𝓂.timings.past_not_future_and_mixed) if idx_in_state !== nothing - past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + # state is in deviations, add SS to get levels + past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] end end @@ -9762,7 +9765,7 @@ function create_newton_state_update(𝓂::ℳ) end end - # Combine: 𝔙 = [present; past; shocks] + # Combine: 𝔙 = [present; past; shocks] (all in levels for present/past) vars = vcat(present_guess, past_from_state, shock_values) # Newton iterations @@ -9773,10 +9776,10 @@ function create_newton_state_update(𝓂::ℳ) jacobian = zeros(length(present_guess), n_present) for iter in 1:max_iter - # Update vars with current present guess + # Update vars with current present guess (in levels) vars[1:n_present] = present_guess - # Evaluate residual and jacobian + # Evaluate residual and jacobian (with values in levels) residual_func(residual, params_and_ss, vars) jacobian_func(jacobian, params_and_ss, vars) @@ -9795,18 +9798,19 @@ function create_newton_state_update(𝓂::ℳ) end end - # Construct full state vector in the correct order - result = copy(SS_and_pars[1:nVars]) # Start with SS values + # Construct full state vector in the correct order (in levels first) + result_levels = copy(SS_vars) - # Fill in the present values we just computed + # Fill in the present values we just computed (in levels) for (i, p) in enumerate(present_sorted) idx = findfirst(x -> x == p, 𝓂.timings.var) if idx !== nothing - result[idx] = present_guess[i] + result_levels[idx] = present_guess[i] end end - return result + # Subtract steady state to return deviations from SS + return result_levels - SS_vars end return state_update_newton From 209427f6db898c3d77da0982c5650a4f92559b6e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 22:12:37 +0000 Subject: [PATCH 08/40] Add simulate_newton function for backward looking models with custom initial states Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 65 +++++++++++---- src/get_functions.jl | 114 +++++++++++++++++++++++++++ test/test_backward_looking_models.jl | 40 ++++++++++ 3 files changed, 202 insertions(+), 17 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 967c8dd4b..67d8a2631 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -189,7 +189,7 @@ export plot_irfs!, plot_irf!, plot_IRF!, plot_girf!, plot_simulations!, plot_sim export Normal, Beta, Cauchy, Gamma, InverseGamma -export get_irfs, get_irf, get_IRF, simulate, get_simulation, get_simulations, get_girf +export get_irfs, get_irf, get_IRF, simulate, get_simulation, get_simulations, get_girf, simulate_newton export get_conditional_forecast export get_solution, get_first_order_solution, get_perturbation_solution, get_second_order_solution, get_third_order_solution export get_steady_state, get_SS, get_ss, get_non_stochastic_steady_state, get_stochastic_steady_state, get_SSS, steady_state, SS, SSS, ss, sss @@ -9621,7 +9621,7 @@ end # dispatch_doctor # return [𝐒₁ * aug_state₁̃, 𝐒₁ * aug_state₂̃ + 𝐒₂ * kron_aug_state₁ / 2, 𝐒₁ * aug_state₃̃ + 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + 𝐒₃ * ℒ.kron(kron_aug_state₁,aug_state₁) / 6] # end -function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasionally_binding_constraints::Bool)::Tuple{Function, Bool} +function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasionally_binding_constraints::Bool; levels::Bool = false)::Tuple{Function, Bool} if occasionally_binding_constraints if algorithm == :first_order state_update = 𝓂.solution.perturbation.first_order.state_update_obc @@ -9666,8 +9666,8 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona pruning = true elseif :newton == algorithm @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "Newton algorithm is only available for backward looking models (nFuture_not_past_and_mixed = 0)." - # Create newton-based state update function - state_update = create_newton_state_update(𝓂) + # Create newton-based state update function with levels option + state_update = create_newton_state_update(𝓂, levels = levels) pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" @@ -9681,12 +9681,15 @@ end """ - create_newton_state_update(𝓂::ℳ) + create_newton_state_update(𝓂::ℳ; levels::Bool = false) Create a state update function that uses Newton's method to solve for present values given past values and shocks. Only for backward looking models. + +When `levels = false` (default): Input state is in deviations from SS, output is in deviations. +When `levels = true`: Input state is in levels, output is in levels. Use this for models without a steady state. """ -function create_newton_state_update(𝓂::ℳ) +function create_newton_state_update(𝓂::ℳ; levels::Bool = false) # Get the newton simulation functions residual_func = 𝓂.newton_simulation_functions.residual_func jacobian_func = 𝓂.newton_simulation_functions.jacobian_func @@ -9742,18 +9745,37 @@ function create_newton_state_update(𝓂::ℳ) ss_values = SS_and_pars[indexin(sort(stst), SS_and_pars_names)] params_and_ss = vcat(pars_values, ss_values) - # Initial guess: steady state values for present variables (in levels) - present_guess = copy(SS_present) - - # Get past values from state (which are in deviations from SS) - # Add steady state to convert to levels + # Get past values from state past_from_state = zeros(n_past) for (i, p) in enumerate(past_sorted) idx_in_state = findfirst(x -> x == p, 𝓂.timings.past_not_future_and_mixed) if idx_in_state !== nothing - # state is in deviations, add SS to get levels - past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] + if levels + # state is already in levels, use directly + past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + else + # state is in deviations, add SS to get levels + past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] + end + end + end + + # Initial guess: use past values as starting point for present (in levels) + # For models without SS, this is more robust than using SS values + if levels + present_guess = zeros(n_present) + for (i, p) in enumerate(present_sorted) + idx_in_past = findfirst(x -> x == p, past_sorted) + if idx_in_past !== nothing + present_guess[i] = past_from_state[idx_in_past] + else + # For variables not in past, use 1.0 as default (common for capital, etc.) + present_guess[i] = 1.0 + end end + else + # Use steady state values for present variables (in levels) + present_guess = copy(SS_present) end # Get shock values - reorder to match sorted order @@ -9798,8 +9820,12 @@ function create_newton_state_update(𝓂::ℳ) end end - # Construct full state vector in the correct order (in levels first) - result_levels = copy(SS_vars) + # Construct full state vector in the correct order (in levels) + if levels + result_levels = zeros(nVars) + else + result_levels = copy(SS_vars) + end # Fill in the present values we just computed (in levels) for (i, p) in enumerate(present_sorted) @@ -9809,8 +9835,13 @@ function create_newton_state_update(𝓂::ℳ) end end - # Subtract steady state to return deviations from SS - return result_levels - SS_vars + if levels + # Return in levels (no SS subtraction) + return result_levels + else + # Subtract steady state to return deviations from SS + return result_levels - SS_vars + end end return state_update_newton diff --git a/src/get_functions.jl b/src/get_functions.jl index 9e44e01bf..117f81cab 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1367,6 +1367,120 @@ Wrapper for [`get_irf`](@ref) with `generalised_irf = true`. get_girf(𝓂::ℳ; kwargs...) = get_irf(𝓂; kwargs..., generalised_irf = true) +""" +$(SIGNATURES) +Simulate a backward looking model forward from an initial state using Newton's method. + +This function is designed for backward looking models (models with no forward-looking variables) +that may not have a steady state or stable growth path. It iterates the system forward using +Newton's method to solve the nonlinear equations at each time step. + +# Arguments +- `𝓂`: The model object +- `initial_state`: A `Vector{Float64}` or `KeyedArray` with initial values for state variables in levels. + For `KeyedArray`, keys should be variable names (Symbols). + +# Keyword Arguments +- `periods` [Default: `40`, Type: `Int`]: Number of periods to simulate +- `shocks` [Default: `:simulate`, Type: `Union{Symbol, Matrix{Float64}, KeyedArray{Float64}}`]: + Shock values. Use `:simulate` for random shocks, `:none` for no shocks, or provide a matrix/KeyedArray. +- `parameters` [Default: `nothing`]: Parameter values to use +- `variables` [Default: `:all_excluding_obc`]: Variables to return +- `verbose` [Default: `false`]: Print progress information + +# Returns +- `KeyedArray` with variables in rows and periods in columns, values in levels. + +# Example +```julia +# Model without steady state (random walk) +@model RandomWalk begin + x[0] = x[-1] + sigma * eps_x[x] +end + +@parameters RandomWalk begin + sigma = 0.1 +end + +# Simulate from initial state x = 5.0 +result = simulate_newton(RandomWalk, [5.0], periods = 20) +``` +""" +function simulate_newton(𝓂::ℳ, + initial_state::Union{Vector{Float64}, KeyedArray{Float64}}; + periods::Int = DEFAULT_PERIODS, + shocks::Union{Symbol_input, String_input, Matrix{Float64}, KeyedArray{Float64}} = :simulate, + parameters::ParameterType = nothing, + variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, + verbose::Bool = DEFAULT_VERBOSE, + tol::Tolerances = Tolerances()) + + @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "simulate_newton is only available for backward looking models (nFuture_not_past_and_mixed = 0)." + + opts = merge_calculation_options(tol = tol, verbose = verbose) + + # Solve the model to get newton functions + solve!(𝓂, + parameters = parameters, + opts = opts, + dynamics = true, + algorithm = :newton) + + # Process shocks + shocks_processed, negative_shock, shock_size, periods, _, _ = process_shocks_input(shocks, false, 1.0, periods, 𝓂) + + # Generate shock history + if shocks_processed == :simulate + shock_history = randn(𝓂.timings.nExo, periods) + elseif shocks_processed == :none + shock_history = zeros(𝓂.timings.nExo, periods) + elseif shocks_processed isa Matrix + shock_history = shocks_processed + periods = size(shock_history, 2) + else + shock_history = zeros(𝓂.timings.nExo, periods) + end + + # Convert KeyedArray initial_state to Vector if needed + if initial_state isa KeyedArray + state_vec = zeros(𝓂.timings.nVars) + for (var, val) in zip(axiskeys(initial_state, 1), initial_state) + idx = findfirst(x -> x == var, 𝓂.timings.var) + if idx !== nothing + state_vec[idx] = val + end + end + initial_state = state_vec + end + + @assert length(initial_state) == 𝓂.timings.nVars "Initial state must have $(𝓂.timings.nVars) elements (one for each variable: $(𝓂.timings.var))" + + # Create newton state update in levels mode + state_update = create_newton_state_update(𝓂, levels = true) + + # Simulate forward + Y = zeros(𝓂.timings.nVars, periods) + current_state = copy(initial_state) + + for t in 1:periods + current_state = state_update(current_state, shock_history[:, t]) + Y[:, t] = current_state + end + + # Process variables + var_idx = parse_variables_input_to_index(variables, 𝓂.timings) |> sort + + axis1 = 𝓂.timings.var[var_idx] + + if any(x -> contains(string(x), "◖"), axis1) + axis1_decomposed = decompose_name.(axis1) + axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] + end + + return KeyedArray(Y[var_idx, :]; Variables = axis1, Periods = 1:periods) +end + + diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index e30a7e1c1..bad79474f 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -106,5 +106,45 @@ SolowGrowth = nothing end + @testset "Backward looking model with custom initial state (simulate_newton)" begin + # Test with Solow growth model from custom initial capital + @model SolowGrowth2 begin + k[0] = (1 - delta) * k[-1] + s * k[-1]^alpha * exp(z[-1]) + z[0] = rho * z[-1] + sigma * eps_z[x] + end + + @parameters SolowGrowth2 begin + alpha = 0.3 + delta = 0.1 + s = 0.2 + rho = 0.9 + sigma = 0.01 + end + + # Verify it's a backward looking model + @test SolowGrowth2.timings.nFuture_not_past_and_mixed == 0 + + # Start from initial capital below steady state + k_ss = (0.2/0.1)^(1/(1-0.3)) + initial_k = k_ss * 0.5 # Start at 50% of steady state capital + initial_state_solow = [initial_k, 0.0] # [k, z] + + result_solow = simulate_newton(SolowGrowth2, initial_state_solow, periods = 10, shocks = :none) + + @test size(result_solow, 1) == 2 # 2 variables + @test size(result_solow, 2) == 10 # 10 periods + + # Capital should increase over time (converging toward steady state) + @test result_solow(:k, 1) > initial_k # First period k should increase + @test result_solow(:k, 10) > result_solow(:k, 1) # Later periods should be higher + + # Test simulation with shocks + result_with_shocks = simulate_newton(SolowGrowth2, initial_state_solow, periods = 10, shocks = :simulate) + @test size(result_with_shocks, 1) == 2 + @test size(result_with_shocks, 2) == 10 + + SolowGrowth2 = nothing + end + GC.gc() end From bd7128cbe0bd9dcce1027a5ec58d45fad4837977 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 22:49:18 +0000 Subject: [PATCH 09/40] Integrate initial_state in levels into existing API for backward looking models with newton Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 2 +- src/get_functions.jl | 194 ++++++++++----------------- test/test_backward_looking_models.jl | 78 +++++++++-- 3 files changed, 142 insertions(+), 132 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 67d8a2631..451c60563 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -189,7 +189,7 @@ export plot_irfs!, plot_irf!, plot_IRF!, plot_girf!, plot_simulations!, plot_sim export Normal, Beta, Cauchy, Gamma, InverseGamma -export get_irfs, get_irf, get_IRF, simulate, get_simulation, get_simulations, get_girf, simulate_newton +export get_irfs, get_irf, get_IRF, simulate, get_simulation, get_simulations, get_girf export get_conditional_forecast export get_solution, get_first_order_solution, get_perturbation_solution, get_second_order_solution, get_third_order_solution export get_steady_state, get_SS, get_ss, get_non_stochastic_steady_state, get_stochastic_steady_state, get_SSS, steady_state, SS, SSS, ss, sss diff --git a/src/get_functions.jl b/src/get_functions.jl index 117f81cab..3bd308049 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1244,6 +1244,66 @@ function get_irf(𝓂::ℳ; # end # timeit_debug + # Check if model is backward looking and has no valid steady state or is explosive + is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 + unspecified_initial_state = initial_state == [0.0] + + # For backward looking models, check steady state validity and if model is explosive + no_valid_steady_state = false + is_explosive = false + + if is_backward_looking + # Try to get steady state and check if it's valid + param_vals = isnothing(parameters) ? 𝓂.parameter_values : parameters isa AbstractDict ? [parameters[k] for k in 𝓂.parameters] : parameters + SS_and_pars, (solution_error, _) = get_NSSS_and_parameters(𝓂, Float64.(param_vals), opts = opts) + + # Check if steady state solution failed or contains invalid values + ss_invalid = solution_error > tol.NSSS_acceptance_tol || + isnan(solution_error) || + any(isnan, SS_and_pars) || + any(isinf, SS_and_pars) + + if ss_invalid + no_valid_steady_state = true + else + # Steady state found - check if model is explosive by doing one newton iteration + # First, ensure newton functions are generated + solve!(𝓂, parameters = parameters, opts = opts, dynamics = true, algorithm = :newton) + + # Create newton state update function + state_update_check = create_newton_state_update(𝓂, levels = true) + + # Get steady state values for variables only + SS_vars = SS_and_pars[1:𝓂.timings.nVars] + + # Do one iteration with no shocks starting from steady state + zero_shocks = zeros(𝓂.timings.nExo) + next_state = state_update_check(SS_vars, zero_shocks) + + # If next state differs from SS, model is explosive + is_explosive = !isapprox(next_state, SS_vars, rtol = 1e-8) + end + + no_valid_steady_state = no_valid_steady_state || is_explosive + end + + # For backward looking models without valid steady state or explosive: + # - algorithm must be :newton + # - initial_state must be provided (in levels) + if is_backward_looking && no_valid_steady_state + if algorithm != :newton + @assert false "Model is backward looking with no valid steady state (or is explosive). Use algorithm = :newton and provide initial_state in levels." + end + if unspecified_initial_state + @assert false "Model is backward looking with no valid steady state (or is explosive). Provide initial_state in levels." + end + end + + # Determine if we should use levels mode: + # - backward looking model with newton algorithm and provided initial_state + # - OR backward looking model with no valid steady state (already checked algorithm = :newton and initial_state provided above) + use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state + # @timeit_debug timer "Solve model" begin solve!(𝓂, @@ -1261,10 +1321,12 @@ function get_irf(𝓂::ℳ; reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm, opts = opts) # end # timeit_debug - - unspecified_initial_state = initial_state == [0.0] - - if unspecified_initial_state + + if use_levels_mode + @assert initial_state isa Vector{Float64} "initial_state must be a Vector{Float64} for backward looking models with newton algorithm." + @assert length(initial_state) == 𝓂.timings.nVars "initial_state must have $(𝓂.timings.nVars) elements (one for each variable: $(𝓂.timings.var))" + # initial_state is already in levels, don't modify it + elseif unspecified_initial_state if algorithm == :pruned_second_order initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] elseif algorithm == :pruned_third_order @@ -1295,10 +1357,16 @@ function get_irf(𝓂::ℳ; state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) else - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false, levels = use_levels_mode) end - level = levels ? reference_steady_state + SSS_delta : SSS_delta + # When using levels mode, output is always in levels + level = (levels || use_levels_mode) ? reference_steady_state + SSS_delta : SSS_delta + + # For use_levels_mode, the level adjustment is different since state_update returns levels directly + if use_levels_mode + level = zeros(𝓂.timings.nVars) # No adjustment needed, state_update returns levels + end responses = compute_irf_responses(𝓂, state_update, @@ -1367,120 +1435,6 @@ Wrapper for [`get_irf`](@ref) with `generalised_irf = true`. get_girf(𝓂::ℳ; kwargs...) = get_irf(𝓂; kwargs..., generalised_irf = true) -""" -$(SIGNATURES) -Simulate a backward looking model forward from an initial state using Newton's method. - -This function is designed for backward looking models (models with no forward-looking variables) -that may not have a steady state or stable growth path. It iterates the system forward using -Newton's method to solve the nonlinear equations at each time step. - -# Arguments -- `𝓂`: The model object -- `initial_state`: A `Vector{Float64}` or `KeyedArray` with initial values for state variables in levels. - For `KeyedArray`, keys should be variable names (Symbols). - -# Keyword Arguments -- `periods` [Default: `40`, Type: `Int`]: Number of periods to simulate -- `shocks` [Default: `:simulate`, Type: `Union{Symbol, Matrix{Float64}, KeyedArray{Float64}}`]: - Shock values. Use `:simulate` for random shocks, `:none` for no shocks, or provide a matrix/KeyedArray. -- `parameters` [Default: `nothing`]: Parameter values to use -- `variables` [Default: `:all_excluding_obc`]: Variables to return -- `verbose` [Default: `false`]: Print progress information - -# Returns -- `KeyedArray` with variables in rows and periods in columns, values in levels. - -# Example -```julia -# Model without steady state (random walk) -@model RandomWalk begin - x[0] = x[-1] + sigma * eps_x[x] -end - -@parameters RandomWalk begin - sigma = 0.1 -end - -# Simulate from initial state x = 5.0 -result = simulate_newton(RandomWalk, [5.0], periods = 20) -``` -""" -function simulate_newton(𝓂::ℳ, - initial_state::Union{Vector{Float64}, KeyedArray{Float64}}; - periods::Int = DEFAULT_PERIODS, - shocks::Union{Symbol_input, String_input, Matrix{Float64}, KeyedArray{Float64}} = :simulate, - parameters::ParameterType = nothing, - variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, - verbose::Bool = DEFAULT_VERBOSE, - tol::Tolerances = Tolerances()) - - @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "simulate_newton is only available for backward looking models (nFuture_not_past_and_mixed = 0)." - - opts = merge_calculation_options(tol = tol, verbose = verbose) - - # Solve the model to get newton functions - solve!(𝓂, - parameters = parameters, - opts = opts, - dynamics = true, - algorithm = :newton) - - # Process shocks - shocks_processed, negative_shock, shock_size, periods, _, _ = process_shocks_input(shocks, false, 1.0, periods, 𝓂) - - # Generate shock history - if shocks_processed == :simulate - shock_history = randn(𝓂.timings.nExo, periods) - elseif shocks_processed == :none - shock_history = zeros(𝓂.timings.nExo, periods) - elseif shocks_processed isa Matrix - shock_history = shocks_processed - periods = size(shock_history, 2) - else - shock_history = zeros(𝓂.timings.nExo, periods) - end - - # Convert KeyedArray initial_state to Vector if needed - if initial_state isa KeyedArray - state_vec = zeros(𝓂.timings.nVars) - for (var, val) in zip(axiskeys(initial_state, 1), initial_state) - idx = findfirst(x -> x == var, 𝓂.timings.var) - if idx !== nothing - state_vec[idx] = val - end - end - initial_state = state_vec - end - - @assert length(initial_state) == 𝓂.timings.nVars "Initial state must have $(𝓂.timings.nVars) elements (one for each variable: $(𝓂.timings.var))" - - # Create newton state update in levels mode - state_update = create_newton_state_update(𝓂, levels = true) - - # Simulate forward - Y = zeros(𝓂.timings.nVars, periods) - current_state = copy(initial_state) - - for t in 1:periods - current_state = state_update(current_state, shock_history[:, t]) - Y[:, t] = current_state - end - - # Process variables - var_idx = parse_variables_input_to_index(variables, 𝓂.timings) |> sort - - axis1 = 𝓂.timings.var[var_idx] - - if any(x -> contains(string(x), "◖"), axis1) - axis1_decomposed = decompose_name.(axis1) - axis1 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis1_decomposed] - end - - return KeyedArray(Y[var_idx, :]; Variables = axis1, Periods = 1:periods) -end - - diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index bad79474f..978f7a1cd 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -106,7 +106,7 @@ SolowGrowth = nothing end - @testset "Backward looking model with custom initial state (simulate_newton)" begin + @testset "Backward looking model with initial_state in levels" begin # Test with Solow growth model from custom initial capital @model SolowGrowth2 begin k[0] = (1 - delta) * k[-1] + s * k[-1]^alpha * exp(z[-1]) @@ -129,22 +129,78 @@ initial_k = k_ss * 0.5 # Start at 50% of steady state capital initial_state_solow = [initial_k, 0.0] # [k, z] - result_solow = simulate_newton(SolowGrowth2, initial_state_solow, periods = 10, shocks = :none) + # For backward looking models with newton algorithm, initial_state is in levels + result = get_irf(SolowGrowth2, + algorithm = :newton, + initial_state = initial_state_solow, + shocks = :none, + periods = 10) - @test size(result_solow, 1) == 2 # 2 variables - @test size(result_solow, 2) == 10 # 10 periods + @test size(result, 1) == 2 # 2 variables + @test size(result, 2) == 10 # 10 periods # Capital should increase over time (converging toward steady state) - @test result_solow(:k, 1) > initial_k # First period k should increase - @test result_solow(:k, 10) > result_solow(:k, 1) # Later periods should be higher - - # Test simulation with shocks - result_with_shocks = simulate_newton(SolowGrowth2, initial_state_solow, periods = 10, shocks = :simulate) - @test size(result_with_shocks, 1) == 2 - @test size(result_with_shocks, 2) == 10 + @test result(:k, 1) > initial_k # First period k should increase + @test result(:k, 10) > result(:k, 1) # Later periods should be higher SolowGrowth2 = nothing end + @testset "Model without analytical steady state (explosive growth)" begin + # Explosive growth model - no stable steady state + # y[0] = (1 + g) * y[-1] + a * z[-1] where g > 0 leads to explosive growth + # This model doesn't have a finite steady state - it grows forever + @model ExplosiveGrowth begin + y[0] = (1 + g) * y[-1] + a * z[-1] + z[0] = rho * z[-1] + sigma * eps_z[x] + end + + @parameters ExplosiveGrowth begin + g = 0.02 # 2% growth rate - no convergence to steady state + a = 0.5 + rho = 0.9 + sigma = 0.01 + end + + # Verify it's a backward looking model + @test ExplosiveGrowth.timings.nFuture_not_past_and_mixed == 0 + + # For explosive backward looking models: + # - must use algorithm = :newton + # - must provide initial_state (in levels) + initial_y = 100.0 # Start at y = 100 + initial_state = [initial_y, 0.0] # [y, z] + + result = get_irf(ExplosiveGrowth, + algorithm = :newton, + initial_state = initial_state, + shocks = :none, + periods = 10) + + @test size(result, 1) == 2 # 2 variables + @test size(result, 2) == 10 # 10 periods + + # Access values properly (result is 3D with shocks dimension) + y_values = result(:y, :, :none) + z_values = result(:z, :, :none) + + # y should grow each period at rate g (with z = 0) + @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 + @test y_values[2] > y_values[1] # Continues growing + @test y_values[10] > y_values[1] # Still growing at period 10 + + # z should stay at 0 with no shocks + @test isapprox(z_values[1], 0.0, atol=1e-10) + @test isapprox(z_values[10], 0.0, atol=1e-10) + + # Test that calling without initial_state on explosive model throws error + @test_throws AssertionError get_irf(ExplosiveGrowth, algorithm = :newton) + + # Test that calling with wrong algorithm on explosive model throws error + @test_throws AssertionError get_irf(ExplosiveGrowth, algorithm = :first_order, initial_state = initial_state) + + ExplosiveGrowth = nothing + end + GC.gc() end From 8bd65fe1d8532bb8ff7a9ea6f9b99f139c5282a4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 23:50:38 +0000 Subject: [PATCH 10/40] Add unit root model support and fix SS equation parsing for backward looking models Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 17 ++++++++-- src/macros.jl | 11 +++++-- test/test_backward_looking_models.jl | 48 ++++++++++++++++++++++++++++ 3 files changed, 71 insertions(+), 5 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index 3bd308049..5c3b48b04 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1251,6 +1251,7 @@ function get_irf(𝓂::ℳ; # For backward looking models, check steady state validity and if model is explosive no_valid_steady_state = false is_explosive = false + has_unit_root = false if is_backward_looking # Try to get steady state and check if it's valid @@ -1282,6 +1283,16 @@ function get_irf(𝓂::ℳ; # If next state differs from SS, model is explosive is_explosive = !isapprox(next_state, SS_vars, rtol = 1e-8) + + # For unit root models with valid SS: the model stays at SS without shocks + # but deviates permanently with shocks (like random walk) + # These are NOT explosive - they have valid SS that the newton solver found + # Unit root detection: check if the first-order solution matrix has eigenvalue = 1 + if !is_explosive && algorithm == :newton + # For unit root models, we can use SS as reference + # Deviations will be computed against no-shock reference path + has_unit_root = true # Mark as unit root model for later processing + end end no_valid_steady_state = no_valid_steady_state || is_explosive @@ -1300,9 +1311,9 @@ function get_irf(𝓂::ℳ; end # Determine if we should use levels mode: - # - backward looking model with newton algorithm and provided initial_state - # - OR backward looking model with no valid steady state (already checked algorithm = :newton and initial_state provided above) - use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state + # - backward looking model with newton algorithm and provided initial_state (for explosive models) + # - For unit root models with valid SS and newton algorithm, we can work in deviations from no-shock path + use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state && !has_unit_root # @timeit_debug timer "Solve model" begin diff --git a/src/macros.jl b/src/macros.jl index ce08c0132..3555e6258 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -647,13 +647,20 @@ macro model(𝓂,ex...) if precompile ss_aux_equation = unblock(prs_ex) else - ss_aux_equation = simplify(unblock(prs_ex)) + simplified_eq = simplify(unblock(prs_ex)) + # If equation simplifies to a constant (Int/Float), keep the original unsimplified form + # This preserves variable dependencies for the incidence matrix + if simplified_eq isa Int || simplified_eq isa Float64 + ss_aux_equation = unblock(prs_ex) + else + ss_aux_equation = simplified_eq + end end end if ss_aux_equation isa Symbol push!(ss_aux_equations, Expr(:call,:-,ss_aux_equation,0)) - else#if !(ss_aux_equation isa Int) + else # println(eq) push!(ss_aux_equations, ss_aux_equation) end diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 978f7a1cd..6838cc076 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -202,5 +202,53 @@ ExplosiveGrowth = nothing end + @testset "Unit root model (AR(1) with rho=1, has valid SS at zero)" begin + # Unit root model using AR(1) formulation with rho=1 + # y[0] = rho_y * y[-1] + sigma * eps_y[x] with rho_y = 1.0 + # This has a valid steady state at y = 0 but has a unit root + # Written this way (instead of y[0] = y[-1] + e[x]) so SymPy can solve SS + @model UnitRootAR1 begin + y[0] = rho_y * y[-1] + sigma * eps_y[x] + z[0] = rho_z * z[-1] + sigma_z * eps_z[x] + end + + @parameters UnitRootAR1 begin + rho_y = 1.0 # Unit root - permanent shocks + rho_z = 0.5 # Stationary + sigma = 0.1 + sigma_z = 0.1 + end + + # Verify it's a backward looking model + @test UnitRootAR1.timings.nFuture_not_past_and_mixed == 0 + + # Should find steady state at zero + ss = get_steady_state(UnitRootAR1) + @test isapprox(ss(:y, :Steady_state), 0.0, atol=1e-10) + @test isapprox(ss(:z, :Steady_state), 0.0, atol=1e-10) + + # For unit root models with valid SS, newton algorithm should work WITHOUT requiring initial_state + # Deviations are computed against the no-shock reference path + irf_newton = get_irf(UnitRootAR1, algorithm = :newton) + @test size(irf_newton, 1) == 2 # 2 variables + @test size(irf_newton, 3) == 2 # 2 shocks + + # The IRF for y to eps_y shock should show permanent effect (unit root) + # First period deviation should be equal to sigma (shock size = 1) + @test isapprox(irf_newton(:y, 1, :eps_y), 0.1, rtol=1e-6) # sigma * 1 + # Effect should persist (unit root with rho=1) + @test isapprox(irf_newton(:y, 10, :eps_y), 0.1, rtol=1e-6) # Still 0.1 + + # For z with rho = 0.5, effect should decay + @test isapprox(irf_newton(:z, 1, :eps_z), 0.1, rtol=1e-6) # sigma_z * 1 + @test irf_newton(:z, 10, :eps_z) < irf_newton(:z, 1, :eps_z) # Decays + + # First order and newton should give same results for linear model + irf_first = get_irf(UnitRootAR1, algorithm = :first_order) + @test isapprox(collect(irf_first), collect(irf_newton), rtol=1e-6) + + UnitRootAR1 = nothing + end + GC.gc() end From 61d8d604ca1da1236c5ecfc5330cb928c4954cdb Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 1 Jan 2026 23:58:55 +0000 Subject: [PATCH 11/40] Refactor backward looking model detection: rename has_unit_root to has_valid_ss_for_reference and fix tests Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 21 ++++------ test/test_backward_looking_models.jl | 61 ++++++++++++++-------------- 2 files changed, 38 insertions(+), 44 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index 5c3b48b04..44bd5a69b 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1251,7 +1251,6 @@ function get_irf(𝓂::ℳ; # For backward looking models, check steady state validity and if model is explosive no_valid_steady_state = false is_explosive = false - has_unit_root = false if is_backward_looking # Try to get steady state and check if it's valid @@ -1283,21 +1282,15 @@ function get_irf(𝓂::ℳ; # If next state differs from SS, model is explosive is_explosive = !isapprox(next_state, SS_vars, rtol = 1e-8) - - # For unit root models with valid SS: the model stays at SS without shocks - # but deviates permanently with shocks (like random walk) - # These are NOT explosive - they have valid SS that the newton solver found - # Unit root detection: check if the first-order solution matrix has eigenvalue = 1 - if !is_explosive && algorithm == :newton - # For unit root models, we can use SS as reference - # Deviations will be computed against no-shock reference path - has_unit_root = true # Mark as unit root model for later processing - end end no_valid_steady_state = no_valid_steady_state || is_explosive end + # For backward looking models: determine if we have a valid SS to use as reference + # If SS is valid (not explosive), we can work in deviations from SS + has_valid_ss_for_reference = is_backward_looking && !no_valid_steady_state + # For backward looking models without valid steady state or explosive: # - algorithm must be :newton # - initial_state must be provided (in levels) @@ -1311,9 +1304,9 @@ function get_irf(𝓂::ℳ; end # Determine if we should use levels mode: - # - backward looking model with newton algorithm and provided initial_state (for explosive models) - # - For unit root models with valid SS and newton algorithm, we can work in deviations from no-shock path - use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state && !has_unit_root + # - backward looking model with newton algorithm and provided initial_state (for explosive models without valid SS) + # - For backward looking models with valid SS and newton algorithm, we work in deviations from SS (not levels) + use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state && !has_valid_ss_for_reference # @timeit_debug timer "Solve model" begin diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 6838cc076..2e4fc2989 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -146,60 +146,61 @@ SolowGrowth2 = nothing end - @testset "Model without analytical steady state (explosive growth)" begin - # Explosive growth model - no stable steady state - # y[0] = (1 + g) * y[-1] + a * z[-1] where g > 0 leads to explosive growth - # This model doesn't have a finite steady state - it grows forever - @model ExplosiveGrowth begin + @testset "Model with unstable eigenvalue but valid SS at zero" begin + # Model y[0] = (1 + g) * y[-1] with g > 0 + # Has eigenvalue > 1 (unstable) but valid SS at y = 0 + # From SS, iterating stays at SS. From any other point, it explodes. + @model UnstableButValidSS begin y[0] = (1 + g) * y[-1] + a * z[-1] z[0] = rho * z[-1] + sigma * eps_z[x] end - @parameters ExplosiveGrowth begin - g = 0.02 # 2% growth rate - no convergence to steady state + @parameters UnstableButValidSS begin + g = 0.02 # 2% growth rate - eigenvalue > 1 a = 0.5 rho = 0.9 sigma = 0.01 end # Verify it's a backward looking model - @test ExplosiveGrowth.timings.nFuture_not_past_and_mixed == 0 + @test UnstableButValidSS.timings.nFuture_not_past_and_mixed == 0 - # For explosive backward looking models: - # - must use algorithm = :newton - # - must provide initial_state (in levels) + # This model HAS a valid SS at y=0, z=0 + ss = get_steady_state(UnstableButValidSS) + @test isapprox(ss(:y, :Steady_state), 0.0, atol=1e-10) + @test isapprox(ss(:z, :Steady_state), 0.0, atol=1e-10) + + # Since SS is valid, newton works WITHOUT requiring initial_state + # (starts from SS=0 and computes IRFs in deviations) + irf_newton = get_irf(UnstableButValidSS, algorithm = :newton) + @test size(irf_newton, 1) == 2 # 2 variables + @test size(irf_newton, 3) == 1 # 1 shock + + # IRF from SS should show the effect of a shock + # y responds to z shock through the 'a' parameter + @test irf_newton(:y, 1, :eps_z) != 0.0 # y responds to z shock + + # Can also provide initial_state in levels for simulation from non-SS point initial_y = 100.0 # Start at y = 100 - initial_state = [initial_y, 0.0] # [y, z] + initial_state_levels = [initial_y, 0.0] # [y, z] - result = get_irf(ExplosiveGrowth, + result = get_irf(UnstableButValidSS, algorithm = :newton, - initial_state = initial_state, + initial_state = initial_state_levels, shocks = :none, periods = 10) @test size(result, 1) == 2 # 2 variables @test size(result, 2) == 10 # 10 periods - # Access values properly (result is 3D with shocks dimension) - y_values = result(:y, :, :none) - z_values = result(:z, :, :none) - # y should grow each period at rate g (with z = 0) + # Since initial_state is in levels and model has valid SS, + # result is in levels (initial_state_levels provided) + y_values = result(:y, :, :none) @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 @test y_values[2] > y_values[1] # Continues growing - @test y_values[10] > y_values[1] # Still growing at period 10 - - # z should stay at 0 with no shocks - @test isapprox(z_values[1], 0.0, atol=1e-10) - @test isapprox(z_values[10], 0.0, atol=1e-10) - - # Test that calling without initial_state on explosive model throws error - @test_throws AssertionError get_irf(ExplosiveGrowth, algorithm = :newton) - - # Test that calling with wrong algorithm on explosive model throws error - @test_throws AssertionError get_irf(ExplosiveGrowth, algorithm = :first_order, initial_state = initial_state) - ExplosiveGrowth = nothing + UnstableButValidSS = nothing end @testset "Unit root model (AR(1) with rho=1, has valid SS at zero)" begin From 9135814e1ff9820194c01a9d91f214be2ebde000 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 07:11:12 +0000 Subject: [PATCH 12/40] Add conditional forecasting support for newton algorithm (partial implementation) Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 22 +++++--- src/get_functions.jl | 76 ++++++++++++++++++++++++++++ test/test_backward_looking_models.jl | 44 ++++++++++++++++ 3 files changed, 134 insertions(+), 8 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 451c60563..f7685b769 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -845,9 +845,11 @@ end function minimize_distance_to_conditions(X::Vector{S}, p)::S where S Conditions, State_update, Shocks, Cond_var_idx, Free_shock_idx, State, Pruning, precision_factor = p - Shocks[Free_shock_idx] .= X + # Create a copy of Shocks with the correct element type to support ForwardDiff + Shocks_typed = convert(Vector{S}, copy(Shocks)) + Shocks_typed[Free_shock_idx] .= X - new_State = State_update(State, convert(typeof(X), Shocks)) + new_State = State_update(State, Shocks_typed) cond_vars = Pruning ? sum(new_State) : new_State @@ -9779,7 +9781,9 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) end # Get shock values - reorder to match sorted order - shock_values = zeros(n_exo) + # Use similar type as shock to support ForwardDiff dual numbers + shock_values = similar(shock, n_exo) + fill!(shock_values, zero(eltype(shock))) for (i, e) in enumerate(exo_sorted) idx_in_shock = findfirst(x -> x == e, 𝓂.timings.exo) if idx_in_shock !== nothing @@ -9788,14 +9792,16 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) end # Combine: 𝔙 = [present; past; shocks] (all in levels for present/past) - vars = vcat(present_guess, past_from_state, shock_values) + # Convert present_guess and past_from_state to match shock type for ForwardDiff + ShockType = eltype(shock) + vars = vcat(convert(Vector{ShockType}, present_guess), convert(Vector{ShockType}, past_from_state), shock_values) # Newton iterations max_iter = 50 tol = 1e-10 - residual = zeros(length(present_guess)) - jacobian = zeros(length(present_guess), n_present) + residual = zeros(ShockType, length(present_guess)) + jacobian = zeros(ShockType, length(present_guess), n_present) for iter in 1:max_iter # Update vars with current present guess (in levels) @@ -9822,9 +9828,9 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) # Construct full state vector in the correct order (in levels) if levels - result_levels = zeros(nVars) + result_levels = zeros(ShockType, nVars) else - result_levels = copy(SS_vars) + result_levels = convert(Vector{ShockType}, copy(SS_vars)) end # Fill in the present values we just computed (in levels) diff --git a/src/get_functions.jl b/src/get_functions.jl index 44bd5a69b..43ece25e5 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -986,6 +986,82 @@ function get_conditional_forecast(𝓂::ℳ, shocks[free_shock_idx,i] = CC \ (conditions[cond_var_idx,i] - state_update(Y[:,i-1], Float64[shocks[:,i]...])[cond_var_idx]) + Y[:,i] = state_update(Y[:,i-1], Float64[shocks[:,i]...]) + end + elseif algorithm == :newton + # For backward looking models with newton algorithm + # Use optimization to find minimum norm shocks that satisfy conditions + precision_factor = 1.0 + + p = (conditions[:,1], state_update, shocks[:,1], cond_var_idx, free_shock_idx, initial_state, false, precision_factor) + + res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), + zeros(length(free_shock_idx)), + Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3)), + Optim.Options(f_abstol = eps(), g_tol= 1e-30); + autodiff = :forward) end + + matched = Optim.minimum(res) < 1e-12 + + if !matched + res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), + zeros(length(free_shock_idx)), + Optim.LBFGS(), + Optim.Options(f_abstol = eps(), g_tol= 1e-30); + autodiff = :forward) end + + matched = Optim.minimum(res) < 1e-12 + end + + @assert matched "Numerical stabiltiy issues for restrictions in period 1." + + x = Optim.minimizer(res) + + shocks[free_shock_idx,1] .= x + + Y[:,1] = state_update(initial_state, Float64[shocks[:,1]...]) + + for i in 2:size(conditions,2) + cond_var_idx = findall(conditions[:,i] .!= nothing) + + if conditions_in_levels + conditions[cond_var_idx,i] .-= reference_steady_state[cond_var_idx] + SSS_delta[cond_var_idx] + else + conditions[cond_var_idx,i] .-= SSS_delta[cond_var_idx] + end + + free_shock_idx = findall(shocks[:,i] .== nothing) + + shocks[free_shock_idx,i] .= 0 + + @assert length(free_shock_idx) >= length(cond_var_idx) "Exact matching only possible with at least as many free shocks than conditioned variables. Period " * repr(i) * " has " * repr(length(free_shock_idx)) * " free shock(s) and " * repr(length(cond_var_idx)) * " conditioned variable(s)." + + p = (conditions[:,i], state_update, shocks[:,i], cond_var_idx, free_shock_idx, Y[:,i-1], false, precision_factor) + + res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), + zeros(length(free_shock_idx)), + Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3)), + Optim.Options(f_abstol = eps(), g_tol= 1e-30); + autodiff = :forward) end + + matched = Optim.minimum(res) < 1e-12 + + if !matched + res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), + zeros(length(free_shock_idx)), + Optim.LBFGS(), + Optim.Options(f_abstol = eps(), g_tol= 1e-30); + autodiff = :forward) end + + matched = Optim.minimum(res) < 1e-12 + end + + @assert matched "Numerical stabiltiy issues for restrictions in period $i." + + x = Optim.minimizer(res) + + shocks[free_shock_idx,i] .= x + Y[:,i] = state_update(Y[:,i-1], Float64[shocks[:,i]...]) end end diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 2e4fc2989..90ed46b1f 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -251,5 +251,49 @@ UnitRootAR1 = nothing end + @testset "Conditional forecasting with newton algorithm" begin + # Test conditional forecasting for backward looking models with newton + @model VAR_cond begin + y[0] = a11 * y[-1] + a12 * x[-1] + sigma_y * eps_y[x] + x[0] = a21 * y[-1] + a22 * x[-1] + sigma_x * eps_x[x] + end + + @parameters VAR_cond begin + a11 = 0.5 + a12 = 0.3 + a21 = 0.2 + a22 = 0.4 + sigma_y = 0.1 + sigma_x = 0.1 + end + + # Verify it's a backward looking model + @test VAR_cond.timings.nFuture_not_past_and_mixed == 0 + + # Test conditional forecast with newton - condition y to be 0.05 in period 1 + conditions = KeyedArray(Matrix{Union{Nothing,Float64}}(nothing, 2, 2), + Variables = VAR_cond.var, + Periods = 1:2) + conditions[:y, 1] = 0.05 # Condition y in period 1 + conditions[:x, 2] = 0.02 # Condition x in period 2 + + # Use newton algorithm for conditional forecast + cf_newton = get_conditional_forecast(VAR_cond, conditions, algorithm = :newton) + + # Check that conditions are met + @test isapprox(cf_newton(:y, 1), 0.05, rtol=1e-6) + @test isapprox(cf_newton(:x, 2), 0.02, rtol=1e-6) + + # Compare with first_order conditional forecast - should give similar results for linear model + cf_first = get_conditional_forecast(VAR_cond, conditions, algorithm = :first_order) + @test isapprox(cf_first(:y, 1), 0.05, rtol=1e-6) + @test isapprox(cf_first(:x, 2), 0.02, rtol=1e-6) + + # Shocks should be similar for linear models + @test isapprox(cf_newton(:eps_y₍ₓ₎, 1), cf_first(:eps_y₍ₓ₎, 1), rtol=1e-4) + + VAR_cond = nothing + end + GC.gc() end From 44dcdbf77018c4cd305442d9f1f463ca8ee5544e Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 11:20:54 +0100 Subject: [PATCH 13/40] Refactor backward looking model handling: streamline state updates and assertions for initial states --- src/MacroModelling.jl | 6 ++---- src/get_functions.jl | 38 ++++++++------------------------------ src/macros.jl | 11 ++--------- 3 files changed, 12 insertions(+), 43 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index f7685b769..44779ac3d 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -845,11 +845,9 @@ end function minimize_distance_to_conditions(X::Vector{S}, p)::S where S Conditions, State_update, Shocks, Cond_var_idx, Free_shock_idx, State, Pruning, precision_factor = p - # Create a copy of Shocks with the correct element type to support ForwardDiff - Shocks_typed = convert(Vector{S}, copy(Shocks)) - Shocks_typed[Free_shock_idx] .= X + Shocks[Free_shock_idx] .= X - new_State = State_update(State, Shocks_typed) + new_State = State_update(State, convert(typeof(X), Shocks)) cond_vars = Pruning ? sum(new_State) : new_State diff --git a/src/get_functions.jl b/src/get_functions.jl index 43ece25e5..b0ed07281 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1359,31 +1359,17 @@ function get_irf(𝓂::ℳ; # If next state differs from SS, model is explosive is_explosive = !isapprox(next_state, SS_vars, rtol = 1e-8) end - - no_valid_steady_state = no_valid_steady_state || is_explosive end - # For backward looking models: determine if we have a valid SS to use as reference - # If SS is valid (not explosive), we can work in deviations from SS - has_valid_ss_for_reference = is_backward_looking && !no_valid_steady_state - # For backward looking models without valid steady state or explosive: # - algorithm must be :newton # - initial_state must be provided (in levels) if is_backward_looking && no_valid_steady_state - if algorithm != :newton - @assert false "Model is backward looking with no valid steady state (or is explosive). Use algorithm = :newton and provide initial_state in levels." - end - if unspecified_initial_state - @assert false "Model is backward looking with no valid steady state (or is explosive). Provide initial_state in levels." - end + @assert algorithm == :newton "Model is backward looking with no valid steady state (or is explosive). Use algorithm = :newton and provide initial_state in levels." + + @assert !unspecified_initial_state "Model is backward looking with no valid steady state (or is explosive). Provide initial_state in levels." end - # Determine if we should use levels mode: - # - backward looking model with newton algorithm and provided initial_state (for explosive models without valid SS) - # - For backward looking models with valid SS and newton algorithm, we work in deviations from SS (not levels) - use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state && !has_valid_ss_for_reference - # @timeit_debug timer "Solve model" begin solve!(𝓂, @@ -1401,12 +1387,10 @@ function get_irf(𝓂::ℳ; reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm, opts = opts) # end # timeit_debug - - if use_levels_mode - @assert initial_state isa Vector{Float64} "initial_state must be a Vector{Float64} for backward looking models with newton algorithm." - @assert length(initial_state) == 𝓂.timings.nVars "initial_state must have $(𝓂.timings.nVars) elements (one for each variable: $(𝓂.timings.var))" - # initial_state is already in levels, don't modify it - elseif unspecified_initial_state + + unspecified_initial_state = initial_state == [0.0] + + if unspecified_initial_state if algorithm == :pruned_second_order initial_state = [zeros(𝓂.timings.nVars), zeros(𝓂.timings.nVars) - SSS_delta] elseif algorithm == :pruned_third_order @@ -1440,13 +1424,7 @@ function get_irf(𝓂::ℳ; state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false, levels = use_levels_mode) end - # When using levels mode, output is always in levels - level = (levels || use_levels_mode) ? reference_steady_state + SSS_delta : SSS_delta - - # For use_levels_mode, the level adjustment is different since state_update returns levels directly - if use_levels_mode - level = zeros(𝓂.timings.nVars) # No adjustment needed, state_update returns levels - end + level = levels ? reference_steady_state + SSS_delta : SSS_delta responses = compute_irf_responses(𝓂, state_update, diff --git a/src/macros.jl b/src/macros.jl index 3555e6258..ce08c0132 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -647,20 +647,13 @@ macro model(𝓂,ex...) if precompile ss_aux_equation = unblock(prs_ex) else - simplified_eq = simplify(unblock(prs_ex)) - # If equation simplifies to a constant (Int/Float), keep the original unsimplified form - # This preserves variable dependencies for the incidence matrix - if simplified_eq isa Int || simplified_eq isa Float64 - ss_aux_equation = unblock(prs_ex) - else - ss_aux_equation = simplified_eq - end + ss_aux_equation = simplify(unblock(prs_ex)) end end if ss_aux_equation isa Symbol push!(ss_aux_equations, Expr(:call,:-,ss_aux_equation,0)) - else + else#if !(ss_aux_equation isa Int) # println(eq) push!(ss_aux_equations, ss_aux_equation) end From 537a2427f5176cc99a64a57e58ffdc615f52db62 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 10:27:43 +0000 Subject: [PATCH 14/40] Refactor: Add backward_looking struct to solution, remove newton_simulation_functions from model Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 91 ++++++++++++++++++------------------------- src/macros.jl | 2 +- src/structures.jl | 14 +++++-- 3 files changed, 50 insertions(+), 57 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 44779ac3d..ea5869c95 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -8118,12 +8118,31 @@ function write_newton_simulation_functions!(𝓂::ℳ; residual_buffer = zeros(Float64, length(dyn_equations)) jacobian_buffer = zeros(Float64, size(∇_present_mat)) - # Store in model - 𝓂.newton_simulation_functions = ( - residual_func = residual_func, - jacobian_func = jacobian_func, - residual_buffer = residual_buffer, - jacobian_buffer = jacobian_buffer + # Compute jacobian w.r.t. shocks for conditional forecasting + shock_vars = 𝔙[n_present + n_past + 1 : n_present + n_past + n_exo] + ∇_shocks = Symbolics.sparsejacobian(dyn_equations_vec, collect(shock_vars)) + ∇_shocks_mat = convert(Matrix, ∇_shocks) + + _, jacobian_shock_func = Symbolics.build_function(∇_shocks_mat, 𝔓, 𝔙, + cse = cse, + skipzeros = skipzeros, + expression_module = @__MODULE__, + expression = Val(false))::Tuple{<:Function, <:Function} + + jacobian_shock_buffer = zeros(Float64, size(∇_shocks_mat)) + + # Create newton state update function + state_update = create_newton_state_update(𝓂, residual_func, jacobian_func, residual_buffer, jacobian_buffer) + + # Store in model's solution.backward_looking struct + 𝓂.solution.backward_looking = backward_looking_solution( + state_update, + residual_func, + jacobian_func, + jacobian_shock_func, + residual_buffer, + jacobian_buffer, + jacobian_shock_buffer ) return nothing @@ -9666,8 +9685,8 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona pruning = true elseif :newton == algorithm @assert 𝓂.timings.nFuture_not_past_and_mixed == 0 "Newton algorithm is only available for backward looking models (nFuture_not_past_and_mixed = 0)." - # Create newton-based state update function with levels option - state_update = create_newton_state_update(𝓂, levels = levels) + # Use the stored backward_looking state update function + state_update = 𝓂.solution.backward_looking.state_update pruning = false else # @assert false "Provided algorithm not valid. Valid algorithm: $all_available_algorithms" @@ -9681,19 +9700,15 @@ end """ - create_newton_state_update(𝓂::ℳ; levels::Bool = false) + create_newton_state_update(𝓂::ℳ, residual_func, jacobian_func, residual_buffer, jacobian_buffer) Create a state update function that uses Newton's method to solve for present values given past values and shocks. Only for backward looking models. -When `levels = false` (default): Input state is in deviations from SS, output is in deviations. -When `levels = true`: Input state is in levels, output is in levels. Use this for models without a steady state. +Works in deviations from SS: Input state is in deviations from SS, output is in deviations. """ -function create_newton_state_update(𝓂::ℳ; levels::Bool = false) - # Get the newton simulation functions - residual_func = 𝓂.newton_simulation_functions.residual_func - jacobian_func = 𝓂.newton_simulation_functions.jacobian_func - +function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian_func::Function, + residual_buffer::Vector{Float64}, jacobian_buffer::Matrix{Float64}) # Get steady state and parameters SS_and_pars = 𝓂.solution.non_stochastic_steady_state parameters = 𝓂.parameter_values @@ -9745,38 +9760,17 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) ss_values = SS_and_pars[indexin(sort(stst), SS_and_pars_names)] params_and_ss = vcat(pars_values, ss_values) - # Get past values from state + # Get past values from state (state is in deviations, add SS to get levels) past_from_state = zeros(n_past) for (i, p) in enumerate(past_sorted) idx_in_state = findfirst(x -> x == p, 𝓂.timings.past_not_future_and_mixed) if idx_in_state !== nothing - if levels - # state is already in levels, use directly - past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] - else - # state is in deviations, add SS to get levels - past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] - end + past_from_state[i] = state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] end end - # Initial guess: use past values as starting point for present (in levels) - # For models without SS, this is more robust than using SS values - if levels - present_guess = zeros(n_present) - for (i, p) in enumerate(present_sorted) - idx_in_past = findfirst(x -> x == p, past_sorted) - if idx_in_past !== nothing - present_guess[i] = past_from_state[idx_in_past] - else - # For variables not in past, use 1.0 as default (common for capital, etc.) - present_guess[i] = 1.0 - end - end - else - # Use steady state values for present variables (in levels) - present_guess = copy(SS_present) - end + # Initial guess: use steady state values for present variables (in levels) + present_guess = copy(SS_present) # Get shock values - reorder to match sorted order # Use similar type as shock to support ForwardDiff dual numbers @@ -9825,11 +9819,7 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) end # Construct full state vector in the correct order (in levels) - if levels - result_levels = zeros(ShockType, nVars) - else - result_levels = convert(Vector{ShockType}, copy(SS_vars)) - end + result_levels = convert(Vector{ShockType}, copy(SS_vars)) # Fill in the present values we just computed (in levels) for (i, p) in enumerate(present_sorted) @@ -9839,13 +9829,8 @@ function create_newton_state_update(𝓂::ℳ; levels::Bool = false) end end - if levels - # Return in levels (no SS subtraction) - return result_levels - else - # Subtract steady state to return deviations from SS - return result_levels - SS_vars - end + # Subtract steady state to return deviations from SS + return result_levels - SS_vars end return state_update_newton diff --git a/src/macros.jl b/src/macros.jl index ce08c0132..14bf4786d 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -913,7 +913,6 @@ macro model(𝓂,ex...) (zeros(0,0), x->x), # third_order_derivatives (zeros(0,0), x->x), # third_order_derivatives_parameters (zeros(0,0), x->x), # third_order_derivatives_SS_and_pars - (residual_func = (x,y,z)->nothing, jacobian_func = (x,y,z)->nothing, residual_buffer = Float64[], jacobian_buffer = zeros(0,0)), # newton_simulation_functions # (x->x, SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0), 𝒟.prepare_jacobian(x->x, 𝒟.AutoForwardDiff(), [0]), SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0)), # third_order_derivatives # ([], SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0)), # model_jacobian # ([], Int[], zeros(1,1)), # model_jacobian @@ -975,6 +974,7 @@ macro model(𝓂,ex...) second_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)), third_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),Dict{Vector{Int}, Int}(),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)) ), + backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0)), Float64[], # Set([:first_order]), Set(all_available_algorithms), diff --git a/src/structures.jl b/src/structures.jl index b887ccfaf..f96871f1c 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -240,8 +240,19 @@ struct ss_solve_block extended_ss_problem::function_and_jacobian end +struct backward_looking_solution + state_update::Function + residual_func::Function + jacobian_func::Function + jacobian_shock_func::Function # Jacobian w.r.t. shocks for conditional forecasting + residual_buffer::Vector{Float64} + jacobian_buffer::Matrix{Float64} + jacobian_shock_buffer::Matrix{Float64} +end + mutable struct solution perturbation::perturbation + backward_looking::backward_looking_solution non_stochastic_steady_state::Vector{Float64} # algorithms::Set{Symbol} outdated_algorithms::Set{Symbol} @@ -426,9 +437,6 @@ mutable struct ℳ third_order_derivatives_parameters::Tuple{AbstractMatrix{<: Real},Function} third_order_derivatives_SS_and_pars::Tuple{AbstractMatrix{<: Real},Function} - # Newton simulation functions for backward looking models - newton_simulation_functions::NamedTuple{(:residual_func, :jacobian_func, :residual_buffer, :jacobian_buffer), Tuple{Function, Function, Vector{Float64}, Matrix{Float64}}} - # model_jacobian::Tuple{Vector{Function}, SparseMatrixCSC{Float64}} # model_jacobian::Tuple{Vector{Function}, Vector{Int}, Matrix{<: Real}} # # model_jacobian_parameters::Function From 2b5f30313158d5e30f7a0f46164fd6878b4044c4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 10:37:57 +0000 Subject: [PATCH 15/40] Fix get_irf to use parse_algorithm_to_state_update for newton and define use_levels_mode Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index b0ed07281..ee3082a75 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1324,6 +1324,10 @@ function get_irf(𝓂::ℳ; is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 unspecified_initial_state = initial_state == [0.0] + # Determine if we should use levels mode for newton + # - Use levels mode when algorithm is newton AND initial_state is provided AND model is backward looking + use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state + # For backward looking models, check steady state validity and if model is explosive no_valid_steady_state = false is_explosive = false @@ -1346,18 +1350,19 @@ function get_irf(𝓂::ℳ; # First, ensure newton functions are generated solve!(𝓂, parameters = parameters, opts = opts, dynamics = true, algorithm = :newton) - # Create newton state update function - state_update_check = create_newton_state_update(𝓂, levels = true) + # Use parse_algorithm_to_state_update to get the newton state update function + state_update_check, _ = parse_algorithm_to_state_update(:newton, 𝓂, false, levels = false) # Get steady state values for variables only SS_vars = SS_and_pars[1:𝓂.timings.nVars] - # Do one iteration with no shocks starting from steady state + # Do one iteration with no shocks starting from steady state (in deviations = zeros) zero_shocks = zeros(𝓂.timings.nExo) - next_state = state_update_check(SS_vars, zero_shocks) + zero_state = zeros(𝓂.timings.nVars) + next_state = state_update_check(zero_state, zero_shocks) - # If next state differs from SS, model is explosive - is_explosive = !isapprox(next_state, SS_vars, rtol = 1e-8) + # If next state differs from zero (SS in deviations), model is explosive + is_explosive = !isapprox(next_state, zero_state, rtol = 1e-8) end end From 89d57f4c999ae1a493832b54cfc010c5c73cb989 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 10:45:14 +0000 Subject: [PATCH 16/40] Remove unused use_levels_mode and levels parameter from parse_algorithm_to_state_update Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 2 +- src/get_functions.jl | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index ea5869c95..f6b4338da 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -9640,7 +9640,7 @@ end # dispatch_doctor # return [𝐒₁ * aug_state₁̃, 𝐒₁ * aug_state₂̃ + 𝐒₂ * kron_aug_state₁ / 2, 𝐒₁ * aug_state₃̃ + 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + 𝐒₃ * ℒ.kron(kron_aug_state₁,aug_state₁) / 6] # end -function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasionally_binding_constraints::Bool; levels::Bool = false)::Tuple{Function, Bool} +function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasionally_binding_constraints::Bool)::Tuple{Function, Bool} if occasionally_binding_constraints if algorithm == :first_order state_update = 𝓂.solution.perturbation.first_order.state_update_obc diff --git a/src/get_functions.jl b/src/get_functions.jl index ee3082a75..4025df90f 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1324,10 +1324,6 @@ function get_irf(𝓂::ℳ; is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 unspecified_initial_state = initial_state == [0.0] - # Determine if we should use levels mode for newton - # - Use levels mode when algorithm is newton AND initial_state is provided AND model is backward looking - use_levels_mode = is_backward_looking && algorithm == :newton && !unspecified_initial_state - # For backward looking models, check steady state validity and if model is explosive no_valid_steady_state = false is_explosive = false @@ -1351,7 +1347,7 @@ function get_irf(𝓂::ℳ; solve!(𝓂, parameters = parameters, opts = opts, dynamics = true, algorithm = :newton) # Use parse_algorithm_to_state_update to get the newton state update function - state_update_check, _ = parse_algorithm_to_state_update(:newton, 𝓂, false, levels = false) + state_update_check, _ = parse_algorithm_to_state_update(:newton, 𝓂, false) # Get steady state values for variables only SS_vars = SS_and_pars[1:𝓂.timings.nVars] @@ -1426,7 +1422,7 @@ function get_irf(𝓂::ℳ; state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, true) else - state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false, levels = use_levels_mode) + state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) end level = levels ? reference_steady_state + SSS_delta : SSS_delta From d20b6f29c554d1285183bcb8f2f3b5bc4d0a6f69 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 10:50:55 +0000 Subject: [PATCH 17/40] Implement LagrangeNewton for conditional forecasting with newton algorithm Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 139 +++++++++++++++++++++++-------------------- 1 file changed, 74 insertions(+), 65 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index 4025df90f..616e0933f 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -990,79 +990,88 @@ function get_conditional_forecast(𝓂::ℳ, end elseif algorithm == :newton # For backward looking models with newton algorithm - # Use optimization to find minimum norm shocks that satisfy conditions - precision_factor = 1.0 + # Use LagrangeNewton formulation to find minimum norm shocks that satisfy conditions + # Solve: min ||e||² s.t. y[cond_var_idx] = conditions where y = state_update(past, e) - p = (conditions[:,1], state_update, shocks[:,1], cond_var_idx, free_shock_idx, initial_state, false, precision_factor) - - res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), - zeros(length(free_shock_idx)), - Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3)), - Optim.Options(f_abstol = eps(), g_tol= 1e-30); - autodiff = :forward) end - - matched = Optim.minimum(res) < 1e-12 - - if !matched - res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), - zeros(length(free_shock_idx)), - Optim.LBFGS(), - Optim.Options(f_abstol = eps(), g_tol= 1e-30); - autodiff = :forward) end - - matched = Optim.minimum(res) < 1e-12 + # For linear backward looking models, the shock jacobian is constant + # Get the shock impact matrix from the solution + shock_impact = 𝓂.solution.perturbation.first_order.solution_matrix[:, 𝓂.timings.nPast_not_future_and_mixed+1:end] + + # LagrangeNewton: solve for minimum norm shocks + # The system for period 1: + # target = y[cond_var_idx] where y = state_update(initial_state, shocks) + # For linear model: target = (A * initial_state + B * shocks)[cond_var_idx] + # Minimum norm solution: shocks = B'[cond_var_idx,:] * (B[cond_var_idx,:] * B'[cond_var_idx,:])^{-1} * (target - A * initial_state)[cond_var_idx] + + # Get the state-independent part (from initial_state) + y_no_shock = state_update(initial_state, zeros(length(𝓂.exo))) + + # Extract shock jacobian for conditioned variables + B_cond = shock_impact[cond_var_idx, free_shock_idx] + + # Target residual (what we need from shocks) + target_residual = Float64[conditions[idx, 1] for idx in cond_var_idx] - y_no_shock[cond_var_idx] + + # Minimum norm solution using pseudoinverse + # x = B' * (B * B')^{-1} * target_residual + BBt = B_cond * B_cond' + + if length(cond_var_idx) == length(free_shock_idx) + # Square system - solve directly + x = B_cond \ target_residual + else + # Underdetermined - minimum norm solution + x = B_cond' * (BBt \ target_residual) end - - @assert matched "Numerical stabiltiy issues for restrictions in period 1." - - x = Optim.minimizer(res) - - shocks[free_shock_idx,1] .= x - - Y[:,1] = state_update(initial_state, Float64[shocks[:,1]...]) - - for i in 2:size(conditions,2) - cond_var_idx = findall(conditions[:,i] .!= nothing) + + shocks[free_shock_idx, 1] .= x + + Y[:, 1] = state_update(initial_state, Float64[shocks[:, 1]...]) + + # Verify solution + matched = maximum(abs.(Y[cond_var_idx, 1] - Float64[conditions[idx, 1] for idx in cond_var_idx])) < 1e-10 + @assert matched "Numerical stability issues for restrictions in period 1." + + for i in 2:size(conditions, 2) + cond_var_idx = findall(conditions[:, i] .!= nothing) if conditions_in_levels - conditions[cond_var_idx,i] .-= reference_steady_state[cond_var_idx] + SSS_delta[cond_var_idx] + conditions[cond_var_idx, i] .-= reference_steady_state[cond_var_idx] + SSS_delta[cond_var_idx] else - conditions[cond_var_idx,i] .-= SSS_delta[cond_var_idx] + conditions[cond_var_idx, i] .-= SSS_delta[cond_var_idx] end - - free_shock_idx = findall(shocks[:,i] .== nothing) - - shocks[free_shock_idx,i] .= 0 - + + free_shock_idx = findall(shocks[:, i] .== nothing) + + shocks[free_shock_idx, i] .= 0 + @assert length(free_shock_idx) >= length(cond_var_idx) "Exact matching only possible with at least as many free shocks than conditioned variables. Period " * repr(i) * " has " * repr(length(free_shock_idx)) * " free shock(s) and " * repr(length(cond_var_idx)) * " conditioned variable(s)." - - p = (conditions[:,i], state_update, shocks[:,i], cond_var_idx, free_shock_idx, Y[:,i-1], false, precision_factor) - - res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), - zeros(length(free_shock_idx)), - Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3)), - Optim.Options(f_abstol = eps(), g_tol= 1e-30); - autodiff = :forward) end - - matched = Optim.minimum(res) < 1e-12 - - if !matched - res = @suppress begin Optim.optimize(x -> minimize_distance_to_conditions(x, p), - zeros(length(free_shock_idx)), - Optim.LBFGS(), - Optim.Options(f_abstol = eps(), g_tol= 1e-30); - autodiff = :forward) end - - matched = Optim.minimum(res) < 1e-12 + + # Get the state-independent part for this period + y_no_shock = state_update(Y[:, i-1], zeros(length(𝓂.exo))) + + # Extract shock jacobian for conditioned variables + B_cond = shock_impact[cond_var_idx, free_shock_idx] + + # Target residual + target_residual = Float64[conditions[idx, i] for idx in cond_var_idx] - y_no_shock[cond_var_idx] + + # Minimum norm solution + BBt = B_cond * B_cond' + + if length(cond_var_idx) == length(free_shock_idx) + x = B_cond \ target_residual + else + x = B_cond' * (BBt \ target_residual) end - - @assert matched "Numerical stabiltiy issues for restrictions in period $i." - - x = Optim.minimizer(res) - - shocks[free_shock_idx,i] .= x - - Y[:,i] = state_update(Y[:,i-1], Float64[shocks[:,i]...]) + + shocks[free_shock_idx, i] .= x + + Y[:, i] = state_update(Y[:, i-1], Float64[shocks[:, i]...]) + + # Verify solution + matched = maximum(abs.(Y[cond_var_idx, i] - Float64[conditions[idx, i] for idx in cond_var_idx])) < 1e-10 + @assert matched "Numerical stability issues for restrictions in period $i." end end From dd21265fa988cc48f031f147ac0cf71379779782 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 11:03:29 +0000 Subject: [PATCH 18/40] Use Newton iterations with finite differences for nonlinear conditional forecasting Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 132 +++++++++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 50 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index 616e0933f..f5af7cc7a 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -993,44 +993,85 @@ function get_conditional_forecast(𝓂::ℳ, # Use LagrangeNewton formulation to find minimum norm shocks that satisfy conditions # Solve: min ||e||² s.t. y[cond_var_idx] = conditions where y = state_update(past, e) - # For linear backward looking models, the shock jacobian is constant - # Get the shock impact matrix from the solution - shock_impact = 𝓂.solution.perturbation.first_order.solution_matrix[:, 𝓂.timings.nPast_not_future_and_mixed+1:end] - - # LagrangeNewton: solve for minimum norm shocks - # The system for period 1: - # target = y[cond_var_idx] where y = state_update(initial_state, shocks) - # For linear model: target = (A * initial_state + B * shocks)[cond_var_idx] - # Minimum norm solution: shocks = B'[cond_var_idx,:] * (B[cond_var_idx,:] * B'[cond_var_idx,:])^{-1} * (target - A * initial_state)[cond_var_idx] - - # Get the state-independent part (from initial_state) - y_no_shock = state_update(initial_state, zeros(length(𝓂.exo))) - - # Extract shock jacobian for conditioned variables - B_cond = shock_impact[cond_var_idx, free_shock_idx] - - # Target residual (what we need from shocks) - target_residual = Float64[conditions[idx, 1] for idx in cond_var_idx] - y_no_shock[cond_var_idx] - - # Minimum norm solution using pseudoinverse - # x = B' * (B * B')^{-1} * target_residual - BBt = B_cond * B_cond' + # Helper function to find shocks using Newton iterations with nonlinear jacobian + function find_newton_shocks(past_state::Vector{Float64}, + target::Vector{Float64}, + cond_idx::Vector{Int}, + free_idx::Vector{Int}, + fixed_shocks::Vector{Float64}) + n_free = length(free_idx) + n_cond = length(cond_idx) + + # Initial guess for free shocks (zeros) + e_free = zeros(n_free) + + max_iter = 100 + tol = 1e-10 + + for iter in 1:max_iter + # Build full shock vector + e_full = copy(fixed_shocks) + e_full[free_idx] .= e_free + + # Evaluate model at current shock guess + y = state_update(past_state, e_full) + + # Compute residual: y[cond_idx] - target + r = y[cond_idx] - target + + # Check convergence + if ℒ.norm(r) < tol + break + end + + # Get shock jacobian ∂y/∂e using finite differences + # (since jacobian_shock_func gives ∂F/∂e not ∂y/∂e directly) + ε = 1e-7 + J = zeros(n_cond, n_free) + for j in 1:n_free + e_plus = copy(e_full) + e_plus[free_idx[j]] += ε + y_plus = state_update(past_state, e_plus) + J[:, j] = (y_plus[cond_idx] - y[cond_idx]) / ε + end + + # Minimum norm update: Δe = J' * (J * J')^{-1} * r + if n_cond == n_free + # Square system + Δe = J \ r + else + # Underdetermined - minimum norm solution + JJt = J * J' + Δe = J' * (JJt \ r) + end + + e_free = e_free - Δe + + if ℒ.norm(Δe) < tol + break + end + end + + return e_free + end - if length(cond_var_idx) == length(free_shock_idx) - # Square system - solve directly - x = B_cond \ target_residual - else - # Underdetermined - minimum norm solution - x = B_cond' * (BBt \ target_residual) + # Period 1 + target_1 = Float64[conditions[idx, 1] for idx in cond_var_idx] + fixed_shocks_1 = zeros(length(𝓂.exo)) + for (idx, val) in enumerate(shocks[:, 1]) + if val !== nothing + fixed_shocks_1[idx] = val + end end - shocks[free_shock_idx, 1] .= x + e_free = find_newton_shocks(initial_state, target_1, cond_var_idx, free_shock_idx, fixed_shocks_1) + shocks[free_shock_idx, 1] .= e_free Y[:, 1] = state_update(initial_state, Float64[shocks[:, 1]...]) # Verify solution - matched = maximum(abs.(Y[cond_var_idx, 1] - Float64[conditions[idx, 1] for idx in cond_var_idx])) < 1e-10 - @assert matched "Numerical stability issues for restrictions in period 1." + matched = maximum(abs.(Y[cond_var_idx, 1] - target_1)) < 1e-8 + @assert matched "Numerical stability issues for restrictions in period 1. Max error: $(maximum(abs.(Y[cond_var_idx, 1] - target_1)))" for i in 2:size(conditions, 2) cond_var_idx = findall(conditions[:, i] .!= nothing) @@ -1047,31 +1088,22 @@ 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 " * repr(i) * " has " * repr(length(free_shock_idx)) * " free shock(s) and " * repr(length(cond_var_idx)) * " conditioned variable(s)." - # Get the state-independent part for this period - y_no_shock = state_update(Y[:, i-1], zeros(length(𝓂.exo))) - - # Extract shock jacobian for conditioned variables - B_cond = shock_impact[cond_var_idx, free_shock_idx] - - # Target residual - target_residual = Float64[conditions[idx, i] for idx in cond_var_idx] - y_no_shock[cond_var_idx] - - # Minimum norm solution - BBt = B_cond * B_cond' - - if length(cond_var_idx) == length(free_shock_idx) - x = B_cond \ target_residual - else - x = B_cond' * (BBt \ target_residual) + target_i = Float64[conditions[idx, i] for idx in cond_var_idx] + fixed_shocks_i = zeros(length(𝓂.exo)) + for (idx, val) in enumerate(shocks[:, i]) + if val !== nothing + fixed_shocks_i[idx] = val + end end - shocks[free_shock_idx, i] .= x + e_free = find_newton_shocks(Y[:, i-1], target_i, cond_var_idx, free_shock_idx, fixed_shocks_i) + shocks[free_shock_idx, i] .= e_free Y[:, i] = state_update(Y[:, i-1], Float64[shocks[:, i]...]) # Verify solution - matched = maximum(abs.(Y[cond_var_idx, i] - Float64[conditions[idx, i] for idx in cond_var_idx])) < 1e-10 - @assert matched "Numerical stability issues for restrictions in period $i." + matched = maximum(abs.(Y[cond_var_idx, i] - target_i)) < 1e-8 + @assert matched "Numerical stability issues for restrictions in period $i. Max error: $(maximum(abs.(Y[cond_var_idx, i] - target_i)))" end end From cbff4e6f118a075d82e618cea869d6e06622076e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 11:14:27 +0000 Subject: [PATCH 19/40] Use analytical jacobian from jacobian_shock_func for conditional forecasting instead of finite differences Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/get_functions.jl | 120 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 109 insertions(+), 11 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index f5af7cc7a..aff311559 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -993,7 +993,110 @@ function get_conditional_forecast(𝓂::ℳ, # Use LagrangeNewton formulation to find minimum norm shocks that satisfy conditions # Solve: min ||e||² s.t. y[cond_var_idx] = conditions where y = state_update(past, e) - # Helper function to find shocks using Newton iterations with nonlinear jacobian + # Get analytical jacobian functions from backward_looking solution + residual_func = 𝓂.solution.backward_looking.residual_func + jacobian_func = 𝓂.solution.backward_looking.jacobian_func + jacobian_shock_func = 𝓂.solution.backward_looking.jacobian_shock_func + + # Get SS and parameter info for evaluating jacobians + SS_and_pars = 𝓂.solution.non_stochastic_steady_state + + dyn_var_present_list = collect(reduce(union, 𝓂.dyn_present_list)) + dyn_var_past_list = collect(reduce(union, 𝓂.dyn_past_list)) + dyn_exo_list = collect(reduce(union, 𝓂.dyn_exo_list)) + dyn_ss_list = Symbol.(string.(collect(reduce(union, 𝓂.dyn_ss_list))) .* "₍ₛₛ₎") + + present = map(x -> Symbol(replace(string(x), r"₍₀₎" => "")), string.(dyn_var_present_list)) + past = map(x -> Symbol(replace(string(x), r"₍₋₁₎" => "")), string.(dyn_var_past_list)) + exo = map(x -> Symbol(replace(string(x), r"₍ₓ₎" => "")), string.(dyn_exo_list)) + stst = map(x -> Symbol(replace(string(x), r"₍ₛₛ₎" => "")), string.(dyn_ss_list)) + + n_present = length(present) + n_past = length(past) + n_exo = length(exo) + nVars = 𝓂.timings.nVars + + pars_ext = vcat(𝓂.parameters, 𝓂.calibration_equations_parameters) + SS_and_pars_names = vcat(Symbol.(string.(sort(union(𝓂.var, 𝓂.exo_past, 𝓂.exo_future)))), 𝓂.calibration_equations_parameters) + + present_sorted = sort(present) + past_sorted = sort(past) + exo_sorted = sort(exo) + + present_idx_in_SS = indexin(present_sorted, SS_and_pars_names) + past_idx_in_SS = indexin(past_sorted, SS_and_pars_names) + + SS_present = SS_and_pars[present_idx_in_SS] + SS_past = SS_and_pars[past_idx_in_SS] + SS_vars = SS_and_pars[1:nVars] + + pars_values = 𝓂.parameter_values + ss_values = SS_and_pars[indexin(sort(stst), SS_and_pars_names)] + params_and_ss = vcat(pars_values, ss_values) + + past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx + + # Helper function to compute solution jacobian ∂y/∂e using implicit function theorem + # F(y, y_past, e) = 0 => ∂y/∂e = -(∂F/∂y)^{-1} * (∂F/∂e) + function compute_solution_shock_jacobian(past_state::Vector{Float64}, y_current::Vector{Float64}, e_current::Vector{Float64}) + # Build vars vector in levels: [present; past; shocks] + past_from_state = zeros(n_past) + for (i, p) in enumerate(past_sorted) + idx_in_state = findfirst(x -> x == p, 𝓂.timings.past_not_future_and_mixed) + if idx_in_state !== nothing + past_from_state[i] = past_state[past_not_future_and_mixed_idx][idx_in_state] + SS_past[i] + end + end + + # Get present values in levels from y_current (deviations) + SS + present_levels = zeros(n_present) + for (i, p) in enumerate(present_sorted) + idx = findfirst(x -> x == p, 𝓂.timings.var) + if idx !== nothing + present_levels[i] = y_current[idx] + SS_present[i] + end + end + + # Reorder shocks to match sorted order + shock_values = zeros(n_exo) + for (i, e) in enumerate(exo_sorted) + idx_in_shock = findfirst(x -> x == e, 𝓂.timings.exo) + if idx_in_shock !== nothing + shock_values[i] = e_current[idx_in_shock] + end + end + + vars = vcat(present_levels, past_from_state, shock_values) + + # Evaluate analytical jacobians + J_y = zeros(n_present, n_present) + J_e = zeros(n_present, n_exo) + jacobian_func(J_y, params_and_ss, vars) + jacobian_shock_func(J_e, params_and_ss, vars) + + # ∂y/∂e = -(∂F/∂y)^{-1} * (∂F/∂e) + # This gives derivatives in the sorted order; we need to map back to model order + dydt_sorted = -(J_y \ J_e) + + # Map from sorted present vars to model's var order + # and from sorted exo to model's exo order + dydt = zeros(nVars, length(𝓂.timings.exo)) + for (i, p) in enumerate(present_sorted) + var_idx = findfirst(x -> x == p, 𝓂.timings.var) + if var_idx !== nothing + for (j, e) in enumerate(exo_sorted) + exo_idx = findfirst(x -> x == e, 𝓂.timings.exo) + if exo_idx !== nothing + dydt[var_idx, exo_idx] = dydt_sorted[i, j] + end + end + end + end + + return dydt + end + + # Helper function to find shocks using Newton iterations with analytical jacobian function find_newton_shocks(past_state::Vector{Float64}, target::Vector{Float64}, cond_idx::Vector{Int}, @@ -1024,16 +1127,11 @@ function get_conditional_forecast(𝓂::ℳ, break end - # Get shock jacobian ∂y/∂e using finite differences - # (since jacobian_shock_func gives ∂F/∂e not ∂y/∂e directly) - ε = 1e-7 - J = zeros(n_cond, n_free) - for j in 1:n_free - e_plus = copy(e_full) - e_plus[free_idx[j]] += ε - y_plus = state_update(past_state, e_plus) - J[:, j] = (y_plus[cond_idx] - y[cond_idx]) / ε - end + # Get shock jacobian ∂y/∂e using analytical formula via implicit function theorem + dydt_full = compute_solution_shock_jacobian(past_state, y, e_full) + + # Extract submatrix for conditioned vars and free shocks + J = dydt_full[cond_idx, free_idx] # Minimum norm update: Δe = J' * (J * J')^{-1} * r if n_cond == n_free From 78f5dcda02ce141fee4e5ee7b58073c654fc611f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 13:29:25 +0000 Subject: [PATCH 20/40] Address review: rename jacobian_func to jacobian_state_func, move newton code to separate condition Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 24 +++++------ src/get_functions.jl | 97 ++++++++++++++++++++++++++----------------- src/structures.jl | 2 +- 3 files changed, 71 insertions(+), 52 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index f6b4338da..5c6766cb2 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -6859,12 +6859,12 @@ function solve!(𝓂::ℳ; 𝓂.solution.non_stochastic_steady_state = SS_and_pars 𝓂.solution.outdated_NSSS = solution_error > opts.tol.NSSS_acceptance_tol - - # Generate newton simulation functions for backward looking models - if 𝓂.timings.nFuture_not_past_and_mixed == 0 - write_newton_simulation_functions!(𝓂) - 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:newton]) - end + end + + # Generate newton simulation functions for backward looking models + if (algorithm == :newton) && (:newton ∈ 𝓂.solution.outdated_algorithms) && (𝓂.timings.nFuture_not_past_and_mixed == 0) + write_newton_simulation_functions!(𝓂) + 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:newton]) end obc_not_solved = isnothing(𝓂.solution.perturbation.second_order.state_update_obc(zeros(𝓂.timings.nVars), zeros(𝓂.timings.nExo))) @@ -8108,7 +8108,7 @@ function write_newton_simulation_functions!(𝓂::ℳ; ∇_present_mat = convert(Matrix, ∇_present) # Build jacobian function - _, jacobian_func = Symbolics.build_function(∇_present_mat, 𝔓, 𝔙, + _, jacobian_state_func = Symbolics.build_function(∇_present_mat, 𝔓, 𝔙, cse = cse, skipzeros = skipzeros, expression_module = @__MODULE__, @@ -8132,13 +8132,13 @@ function write_newton_simulation_functions!(𝓂::ℳ; jacobian_shock_buffer = zeros(Float64, size(∇_shocks_mat)) # Create newton state update function - state_update = create_newton_state_update(𝓂, residual_func, jacobian_func, residual_buffer, jacobian_buffer) + state_update = create_newton_state_update(𝓂, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer) # Store in model's solution.backward_looking struct 𝓂.solution.backward_looking = backward_looking_solution( state_update, residual_func, - jacobian_func, + jacobian_state_func, jacobian_shock_func, residual_buffer, jacobian_buffer, @@ -9700,14 +9700,14 @@ end """ - create_newton_state_update(𝓂::ℳ, residual_func, jacobian_func, residual_buffer, jacobian_buffer) + create_newton_state_update(𝓂::ℳ, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer) Create a state update function that uses Newton's method to solve for present values given past values and shocks. Only for backward looking models. Works in deviations from SS: Input state is in deviations from SS, output is in deviations. """ -function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian_func::Function, +function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian_state_func::Function, residual_buffer::Vector{Float64}, jacobian_buffer::Matrix{Float64}) # Get steady state and parameters SS_and_pars = 𝓂.solution.non_stochastic_steady_state @@ -9801,7 +9801,7 @@ function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian # Evaluate residual and jacobian (with values in levels) residual_func(residual, params_and_ss, vars) - jacobian_func(jacobian, params_and_ss, vars) + jacobian_state_func(jacobian, params_and_ss, vars) # Check convergence if ℒ.norm(residual) < tol diff --git a/src/get_functions.jl b/src/get_functions.jl index aff311559..9e70659b8 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -995,7 +995,7 @@ function get_conditional_forecast(𝓂::ℳ, # Get analytical jacobian functions from backward_looking solution residual_func = 𝓂.solution.backward_looking.residual_func - jacobian_func = 𝓂.solution.backward_looking.jacobian_func + jacobian_state_func = 𝓂.solution.backward_looking.jacobian_state_func jacobian_shock_func = 𝓂.solution.backward_looking.jacobian_shock_func # Get SS and parameter info for evaluating jacobians @@ -1036,9 +1036,10 @@ function get_conditional_forecast(𝓂::ℳ, past_not_future_and_mixed_idx = 𝓂.timings.past_not_future_and_mixed_idx - # Helper function to compute solution jacobian ∂y/∂e using implicit function theorem + # Helper function to compute Jacobian J = ∂y[cond]/∂e[free] directly using implicit function theorem # F(y, y_past, e) = 0 => ∂y/∂e = -(∂F/∂y)^{-1} * (∂F/∂e) - function compute_solution_shock_jacobian(past_state::Vector{Float64}, y_current::Vector{Float64}, e_current::Vector{Float64}) + # Returns J and also J_y (for algebraic simplification in underdetermined case) + function compute_jacobian_and_update(past_state::Vector{Float64}, y_current::Vector{Float64}, e_current::Vector{Float64}, r::Vector{Float64}, cond_idx::Vector{Int}, free_idx::Vector{Int}) # Build vars vector in levels: [present; past; shocks] past_from_state = zeros(n_past) for (i, p) in enumerate(past_sorted) @@ -1069,31 +1070,63 @@ function get_conditional_forecast(𝓂::ℳ, vars = vcat(present_levels, past_from_state, shock_values) # Evaluate analytical jacobians - J_y = zeros(n_present, n_present) - J_e = zeros(n_present, n_exo) - jacobian_func(J_y, params_and_ss, vars) - jacobian_shock_func(J_e, params_and_ss, vars) + J_y_full = zeros(n_present, n_present) + J_e_full = zeros(n_present, n_exo) + jacobian_state_func(J_y_full, params_and_ss, vars) + jacobian_shock_func(J_e_full, params_and_ss, vars) - # ∂y/∂e = -(∂F/∂y)^{-1} * (∂F/∂e) - # This gives derivatives in the sorted order; we need to map back to model order - dydt_sorted = -(J_y \ J_e) + # Build index maps: cond_idx (model order) -> sorted order for present vars + # free_idx (model order) -> sorted order for exo vars + cond_idx_sorted = Int[] + for c in cond_idx + var_name = 𝓂.timings.var[c] + sorted_idx = findfirst(x -> x == var_name, present_sorted) + if sorted_idx !== nothing + push!(cond_idx_sorted, sorted_idx) + end + end - # Map from sorted present vars to model's var order - # and from sorted exo to model's exo order - dydt = zeros(nVars, length(𝓂.timings.exo)) - for (i, p) in enumerate(present_sorted) - var_idx = findfirst(x -> x == p, 𝓂.timings.var) - if var_idx !== nothing - for (j, e) in enumerate(exo_sorted) - exo_idx = findfirst(x -> x == e, 𝓂.timings.exo) - if exo_idx !== nothing - dydt[var_idx, exo_idx] = dydt_sorted[i, j] - end - end + free_idx_sorted = Int[] + for f in free_idx + exo_name = 𝓂.timings.exo[f] + sorted_idx = findfirst(x -> x == exo_name, exo_sorted) + if sorted_idx !== nothing + push!(free_idx_sorted, sorted_idx) end end - return dydt + n_cond = length(cond_idx) + n_free = length(free_idx) + + # Extract submatrices + J_y = J_y_full[cond_idx_sorted, :] # n_cond x n_present + J_e = J_e_full[:, free_idx_sorted] # n_present x n_free + + # Algebraic simplification: instead of J = -(J_y_full \ J_e_full)[cond_sorted, free_sorted] + # and then Δe = J \ r or J' * (J*J')^{-1} * r + # + # For square case (n_cond == n_free == n_present): + # J \ r where J = -(J_y_full \ J_e_full) + # => solve J_y_full * J * J_e_full^{-1} = -I for J, then J \ r + # => Δe = -J_e_full \ (J_y_full \ r) when system is square + # + # For general case: + # We need the full J matrix to compute minimum norm solution + + # Compute J = -(J_y_full \ J_e_full)[cond_idx_sorted, free_idx_sorted] + dydt_sorted = -(J_y_full \ J_e_full) + J = dydt_sorted[cond_idx_sorted, free_idx_sorted] + + if n_cond == n_free + # Square system + Δe = J \ r + else + # Underdetermined - minimum norm solution: Δe = J' * (J * J')^{-1} * r + JJt = J * J' + Δe = J' * (JJt \ r) + end + + return Δe end # Helper function to find shocks using Newton iterations with analytical jacobian @@ -1103,7 +1136,6 @@ function get_conditional_forecast(𝓂::ℳ, free_idx::Vector{Int}, fixed_shocks::Vector{Float64}) n_free = length(free_idx) - n_cond = length(cond_idx) # Initial guess for free shocks (zeros) e_free = zeros(n_free) @@ -1127,21 +1159,8 @@ function get_conditional_forecast(𝓂::ℳ, break end - # Get shock jacobian ∂y/∂e using analytical formula via implicit function theorem - dydt_full = compute_solution_shock_jacobian(past_state, y, e_full) - - # Extract submatrix for conditioned vars and free shocks - J = dydt_full[cond_idx, free_idx] - - # Minimum norm update: Δe = J' * (J * J')^{-1} * r - if n_cond == n_free - # Square system - Δe = J \ r - else - # Underdetermined - minimum norm solution - JJt = J * J' - Δe = J' * (JJt \ r) - end + # Compute Newton update using analytical jacobian and implicit function theorem + Δe = compute_jacobian_and_update(past_state, y, e_full, r, cond_idx, free_idx) e_free = e_free - Δe diff --git a/src/structures.jl b/src/structures.jl index f96871f1c..ca42446f5 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -243,7 +243,7 @@ end struct backward_looking_solution state_update::Function residual_func::Function - jacobian_func::Function + jacobian_state_func::Function jacobian_shock_func::Function # Jacobian w.r.t. shocks for conditional forecasting residual_buffer::Vector{Float64} jacobian_buffer::Matrix{Float64} From 97548483e929485ebe3898f3e4df1c5ce1d6cdbd Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 15:30:22 +0100 Subject: [PATCH 21/40] Remove redundant condition for newton algorithm in solve! function --- src/MacroModelling.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 5c6766cb2..2b392a24c 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -6795,8 +6795,7 @@ 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))) || - ((:newton == algorithm) && ((:newton ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) + ((:pruned_third_order == algorithm) && ((:pruned_third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) # @timeit_debug timer "Solve for NSSS (if necessary)" begin From b2cbed8cd551fe9335594d33da85de0ab8b3ffc4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 15:08:31 +0000 Subject: [PATCH 22/40] Address review: Use LinearSolve caching for Newton, skip find_newton_shocks when no conditions Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 29 +++++++++++++++++++---------- src/get_functions.jl | 8 ++++++++ src/macros.jl | 2 +- src/structures.jl | 1 + 4 files changed, 29 insertions(+), 11 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 2b392a24c..c172c83c3 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -8117,6 +8117,10 @@ function write_newton_simulation_functions!(𝓂::ℳ; residual_buffer = zeros(Float64, length(dyn_equations)) jacobian_buffer = zeros(Float64, size(∇_present_mat)) + # Create LinearSolve cache for Newton iterations + prob = 𝒮.LinearProblem(jacobian_buffer, residual_buffer) + lu_buffer = 𝒮.init(prob, 𝒮.LUFactorization(), verbose = isdefined(𝒮, :LinearVerbosity) ? 𝒮.LinearVerbosity(𝒮.SciMLLogging.Minimal()) : false) + # Compute jacobian w.r.t. shocks for conditional forecasting shock_vars = 𝔙[n_present + n_past + 1 : n_present + n_past + n_exo] ∇_shocks = Symbolics.sparsejacobian(dyn_equations_vec, collect(shock_vars)) @@ -8131,7 +8135,7 @@ function write_newton_simulation_functions!(𝓂::ℳ; jacobian_shock_buffer = zeros(Float64, size(∇_shocks_mat)) # Create newton state update function - state_update = create_newton_state_update(𝓂, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer) + state_update = create_newton_state_update(𝓂, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer, lu_buffer) # Store in model's solution.backward_looking struct 𝓂.solution.backward_looking = backward_looking_solution( @@ -8141,7 +8145,8 @@ function write_newton_simulation_functions!(𝓂::ℳ; jacobian_shock_func, residual_buffer, jacobian_buffer, - jacobian_shock_buffer + jacobian_shock_buffer, + lu_buffer ) return nothing @@ -9699,7 +9704,7 @@ end """ - create_newton_state_update(𝓂::ℳ, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer) + create_newton_state_update(𝓂::ℳ, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer, lu_buffer) Create a state update function that uses Newton's method to solve for present values given past values and shocks. Only for backward looking models. @@ -9707,7 +9712,8 @@ given past values and shocks. Only for backward looking models. Works in deviations from SS: Input state is in deviations from SS, output is in deviations. """ function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian_state_func::Function, - residual_buffer::Vector{Float64}, jacobian_buffer::Matrix{Float64}) + residual_buffer::Vector{Float64}, jacobian_buffer::Matrix{Float64}, + lu_buffer::𝒮.LinearCache) # Get steady state and parameters SS_and_pars = 𝓂.solution.non_stochastic_steady_state parameters = 𝓂.parameter_values @@ -9807,14 +9813,17 @@ function create_newton_state_update(𝓂::ℳ, residual_func::Function, jacobian break end - # Newton step: present_new = present_old - J^{-1} * F - try - Δ = jacobian \ residual - present_guess = present_guess - Δ - catch - # If jacobian is singular, break + # Newton step using LinearSolve: present_new = present_old - J^{-1} * F + lu_buffer.A = jacobian + lu_buffer.b = residual + 𝒮.solve!(lu_buffer) + + if !isfinite(sum(lu_buffer.u)) + # If solution is not finite (singular jacobian), break break end + + present_guess = present_guess - lu_buffer.u end # Construct full state vector in the correct order (in levels) diff --git a/src/get_functions.jl b/src/get_functions.jl index 9e70659b8..f78841264 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1193,6 +1193,14 @@ function get_conditional_forecast(𝓂::ℳ, for i in 2:size(conditions, 2) cond_var_idx = findall(conditions[:, i] .!= nothing) + # Skip find_newton_shocks if no conditions in this period - propagate state with zero shocks + if isempty(cond_var_idx) + free_shock_idx = findall(shocks[:, i] .== nothing) + shocks[free_shock_idx, i] .= 0 + Y[:, i] = state_update(Y[:, i-1], Float64[shocks[:, i]...]) + continue + end + if conditions_in_levels conditions[cond_var_idx, i] .-= reference_steady_state[cond_var_idx] + SSS_delta[cond_var_idx] else diff --git a/src/macros.jl b/src/macros.jl index 14bf4786d..768923ec5 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -974,7 +974,7 @@ macro model(𝓂,ex...) second_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)), third_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),Dict{Vector{Int}, Int}(),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)) ), - backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0)), + backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization())), Float64[], # Set([:first_order]), Set(all_available_algorithms), diff --git a/src/structures.jl b/src/structures.jl index ca42446f5..0cd07f3ea 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -248,6 +248,7 @@ struct backward_looking_solution residual_buffer::Vector{Float64} jacobian_buffer::Matrix{Float64} jacobian_shock_buffer::Matrix{Float64} + lu_buffer::𝒮.LinearCache end mutable struct solution From 16aedf48c3c7fad0f52e73de45b35fc0187683a7 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 17:41:16 +0100 Subject: [PATCH 23/40] Refactor get_conditional_forecast to use nVars for Y initialization instead of solution_matrix size --- src/get_functions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index f78841264..baa917140 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -848,7 +848,7 @@ function get_conditional_forecast(𝓂::ℳ, var_idx = parse_variables_input_to_index(variables, 𝓂.timings) |> sort - Y = zeros(size(𝓂.solution.perturbation.first_order.solution_matrix,1),periods) + Y = zeros(𝓂.timings.nVars, periods) cond_var_idx = findall(conditions[:,1] .!= nothing) From 4f66cdad2887d99131fa5f355324e0a588e7559e Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 17:46:19 +0100 Subject: [PATCH 24/40] Fix backward model validation: update steady state checks and assertions --- src/get_functions.jl | 28 +++------------------------- 1 file changed, 3 insertions(+), 25 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index baa917140..8f84a08f8 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1504,37 +1504,15 @@ function get_irf(𝓂::ℳ; isnan(solution_error) || any(isnan, SS_and_pars) || any(isinf, SS_and_pars) - - if ss_invalid - no_valid_steady_state = true - else - # Steady state found - check if model is explosive by doing one newton iteration - # First, ensure newton functions are generated - solve!(𝓂, parameters = parameters, opts = opts, dynamics = true, algorithm = :newton) - - # Use parse_algorithm_to_state_update to get the newton state update function - state_update_check, _ = parse_algorithm_to_state_update(:newton, 𝓂, false) - - # Get steady state values for variables only - SS_vars = SS_and_pars[1:𝓂.timings.nVars] - - # Do one iteration with no shocks starting from steady state (in deviations = zeros) - zero_shocks = zeros(𝓂.timings.nExo) - zero_state = zeros(𝓂.timings.nVars) - next_state = state_update_check(zero_state, zero_shocks) - - # If next state differs from zero (SS in deviations), model is explosive - is_explosive = !isapprox(next_state, zero_state, rtol = 1e-8) - end end # For backward looking models without valid steady state or explosive: # - algorithm must be :newton # - initial_state must be provided (in levels) - if is_backward_looking && no_valid_steady_state - @assert algorithm == :newton "Model is backward looking with no valid steady state (or is explosive). Use algorithm = :newton and provide initial_state in levels." + if is_backward_looking && ss_invalid + @assert algorithm == :newton "Model is backward looking with no valid steady state. Use algorithm = :newton and provide initial_state in levels." - @assert !unspecified_initial_state "Model is backward looking with no valid steady state (or is explosive). Provide initial_state in levels." + @assert !unspecified_initial_state "Model is backward looking with no valid steady state. Provide initial_state in levels." end # @timeit_debug timer "Solve model" begin From f8e042a19671b0452859df6f250499685959ce37 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 16:52:41 +0000 Subject: [PATCH 25/40] Address review: Use LinearSolve caching for conditional forecasting, let solve! handle parameter parsing Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 17 ++++++++- src/get_functions.jl | 85 +++++++++++++++++++++---------------------- src/macros.jl | 2 +- src/structures.jl | 5 +++ 4 files changed, 63 insertions(+), 46 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index c172c83c3..4092d9902 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -8134,6 +8134,17 @@ function write_newton_simulation_functions!(𝓂::ℳ; jacobian_shock_buffer = zeros(Float64, size(∇_shocks_mat)) + # Create buffers for conditional forecasting (fixed size based on model dimensions) + n_present = length(present_vars) + n_exo = length(shock_vars) + jac_state_buffer = zeros(Float64, n_present, n_present) + jac_shock_buffer = zeros(Float64, n_present, n_exo) + dydt_buffer = zeros(Float64, n_present, n_exo) + + # Create LinearSolve cache for J_y_full \ J_e_full operations in conditional forecasting + jac_state_prob = 𝒮.LinearProblem(jac_state_buffer, jac_shock_buffer) + lu_jac_state_cache = 𝒮.init(jac_state_prob, 𝒮.LUFactorization(), verbose = isdefined(𝒮, :LinearVerbosity) ? 𝒮.LinearVerbosity(𝒮.SciMLLogging.Minimal()) : false) + # Create newton state update function state_update = create_newton_state_update(𝓂, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer, lu_buffer) @@ -8146,7 +8157,11 @@ function write_newton_simulation_functions!(𝓂::ℳ; residual_buffer, jacobian_buffer, jacobian_shock_buffer, - lu_buffer + lu_buffer, + jac_state_buffer, + jac_shock_buffer, + dydt_buffer, + lu_jac_state_cache ) return nothing diff --git a/src/get_functions.jl b/src/get_functions.jl index 8f84a08f8..fb12efa35 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -998,6 +998,12 @@ function get_conditional_forecast(𝓂::ℳ, jacobian_state_func = 𝓂.solution.backward_looking.jacobian_state_func jacobian_shock_func = 𝓂.solution.backward_looking.jacobian_shock_func + # Get buffers and caches from backward_looking solution + jac_state_buffer = 𝓂.solution.backward_looking.jac_state_buffer + jac_shock_buffer = 𝓂.solution.backward_looking.jac_shock_buffer + dydt_buffer = 𝓂.solution.backward_looking.dydt_buffer + lu_jac_state_cache = 𝓂.solution.backward_looking.lu_jac_state_cache + # Get SS and parameter info for evaluating jacobians SS_and_pars = 𝓂.solution.non_stochastic_steady_state @@ -1069,11 +1075,11 @@ function get_conditional_forecast(𝓂::ℳ, vars = vcat(present_levels, past_from_state, shock_values) - # Evaluate analytical jacobians - J_y_full = zeros(n_present, n_present) - J_e_full = zeros(n_present, n_exo) - jacobian_state_func(J_y_full, params_and_ss, vars) - jacobian_shock_func(J_e_full, params_and_ss, vars) + # Evaluate analytical jacobians using buffers + fill!(jac_state_buffer, 0.0) + fill!(jac_shock_buffer, 0.0) + jacobian_state_func(jac_state_buffer, params_and_ss, vars) + jacobian_shock_func(jac_shock_buffer, params_and_ss, vars) # Build index maps: cond_idx (model order) -> sorted order for present vars # free_idx (model order) -> sorted order for exo vars @@ -1098,32 +1104,26 @@ function get_conditional_forecast(𝓂::ℳ, n_cond = length(cond_idx) n_free = length(free_idx) - # Extract submatrices - J_y = J_y_full[cond_idx_sorted, :] # n_cond x n_present - J_e = J_e_full[:, free_idx_sorted] # n_present x n_free - - # Algebraic simplification: instead of J = -(J_y_full \ J_e_full)[cond_sorted, free_sorted] - # and then Δe = J \ r or J' * (J*J')^{-1} * r - # - # For square case (n_cond == n_free == n_present): - # J \ r where J = -(J_y_full \ J_e_full) - # => solve J_y_full * J * J_e_full^{-1} = -I for J, then J \ r - # => Δe = -J_e_full \ (J_y_full \ r) when system is square - # - # For general case: - # We need the full J matrix to compute minimum norm solution + # Compute J = -(J_y_full \ J_e_full)[cond_idx_sorted, free_idx_sorted] using LinearSolve cache + # Solve J_y_full * X = J_e_full for X (store in dydt_buffer) + lu_jac_state_cache.A = jac_state_buffer + lu_jac_state_cache.b = jac_shock_buffer + 𝒮.solve!(lu_jac_state_cache) + dydt_buffer .= -lu_jac_state_cache.u - # Compute J = -(J_y_full \ J_e_full)[cond_idx_sorted, free_idx_sorted] - dydt_sorted = -(J_y_full \ J_e_full) - J = dydt_sorted[cond_idx_sorted, free_idx_sorted] + J = dydt_buffer[cond_idx_sorted, free_idx_sorted] if n_cond == n_free - # Square system - Δe = J \ r + # Square system - use direct linear solve + J_prob = 𝒮.LinearProblem(J, r) + J_sol = 𝒮.solve(J_prob, 𝒮.LUFactorization()) + Δe = J_sol.u else # Underdetermined - minimum norm solution: Δe = J' * (J * J')^{-1} * r JJt = J * J' - Δe = J' * (JJt \ r) + JJt_prob = 𝒮.LinearProblem(JJt, r) + JJt_sol = 𝒮.solve(JJt_prob, 𝒮.LUFactorization()) + Δe = J' * JJt_sol.u end return Δe @@ -1486,18 +1486,27 @@ function get_irf(𝓂::ℳ; # end # timeit_debug - # Check if model is backward looking and has no valid steady state or is explosive + # Check if model is backward looking is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 unspecified_initial_state = initial_state == [0.0] - # For backward looking models, check steady state validity and if model is explosive - no_valid_steady_state = false - is_explosive = false + # @timeit_debug timer "Solve model" begin + + solve!(𝓂, + parameters = parameters, + opts = opts, + dynamics = true, + algorithm = algorithm, + # timer = timer, + obc = occasionally_binding_constraints || obc_shocks_included) + + # end # timeit_debug + # For backward looking models, check steady state validity after solve! has processed parameters + ss_invalid = false if is_backward_looking - # Try to get steady state and check if it's valid - param_vals = isnothing(parameters) ? 𝓂.parameter_values : parameters isa AbstractDict ? [parameters[k] for k in 𝓂.parameters] : parameters - SS_and_pars, (solution_error, _) = get_NSSS_and_parameters(𝓂, Float64.(param_vals), opts = opts) + # Get steady state using already-updated parameters + SS_and_pars, (solution_error, _) = get_NSSS_and_parameters(𝓂, 𝓂.parameter_values, opts = opts) # Check if steady state solution failed or contains invalid values ss_invalid = solution_error > tol.NSSS_acceptance_tol || @@ -1514,18 +1523,6 @@ function get_irf(𝓂::ℳ; @assert !unspecified_initial_state "Model is backward looking with no valid steady state. Provide initial_state in levels." end - - # @timeit_debug timer "Solve model" begin - - solve!(𝓂, - parameters = parameters, - opts = opts, - dynamics = true, - algorithm = algorithm, - # timer = timer, - obc = occasionally_binding_constraints || obc_shocks_included) - - # end # timeit_debug # @timeit_debug timer "Get relevant steady state" begin diff --git a/src/macros.jl b/src/macros.jl index 768923ec5..0a8ec1036 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -974,7 +974,7 @@ macro model(𝓂,ex...) second_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)), third_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),Dict{Vector{Int}, Int}(),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)) ), - backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization())), + backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization()), zeros(0,0), zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization())), Float64[], # Set([:first_order]), Set(all_available_algorithms), diff --git a/src/structures.jl b/src/structures.jl index 0cd07f3ea..cc60fca91 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -249,6 +249,11 @@ struct backward_looking_solution jacobian_buffer::Matrix{Float64} jacobian_shock_buffer::Matrix{Float64} lu_buffer::𝒮.LinearCache + # Buffers for conditional forecasting + jac_state_buffer::Matrix{Float64} # For J_y_full in conditional forecasting + jac_shock_buffer::Matrix{Float64} # For J_e_full in conditional forecasting + dydt_buffer::Matrix{Float64} # For J_y_full \ J_e_full result + lu_jac_state_cache::𝒮.LinearCache # Cache for J_y_full \ b operations end mutable struct solution From 68bb2c7645139c71dcd12449c460336989c09861 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 17:15:05 +0000 Subject: [PATCH 26/40] Remove LinearSolve from matrix operations in conditional forecasting, keep buffers only Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 11 ++--------- src/get_functions.jl | 23 +++++++---------------- src/macros.jl | 2 +- src/structures.jl | 2 -- 4 files changed, 10 insertions(+), 28 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 4092d9902..0ea014cb9 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -8138,12 +8138,7 @@ function write_newton_simulation_functions!(𝓂::ℳ; n_present = length(present_vars) n_exo = length(shock_vars) jac_state_buffer = zeros(Float64, n_present, n_present) - jac_shock_buffer = zeros(Float64, n_present, n_exo) - dydt_buffer = zeros(Float64, n_present, n_exo) - - # Create LinearSolve cache for J_y_full \ J_e_full operations in conditional forecasting - jac_state_prob = 𝒮.LinearProblem(jac_state_buffer, jac_shock_buffer) - lu_jac_state_cache = 𝒮.init(jac_state_prob, 𝒮.LUFactorization(), verbose = isdefined(𝒮, :LinearVerbosity) ? 𝒮.LinearVerbosity(𝒮.SciMLLogging.Minimal()) : false) + jac_shock_buffer_cond = zeros(Float64, n_present, n_exo) # Create newton state update function state_update = create_newton_state_update(𝓂, residual_func, jacobian_state_func, residual_buffer, jacobian_buffer, lu_buffer) @@ -8159,9 +8154,7 @@ function write_newton_simulation_functions!(𝓂::ℳ; jacobian_shock_buffer, lu_buffer, jac_state_buffer, - jac_shock_buffer, - dydt_buffer, - lu_jac_state_cache + jac_shock_buffer_cond ) return nothing diff --git a/src/get_functions.jl b/src/get_functions.jl index fb12efa35..1a9c52aa4 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -998,11 +998,9 @@ function get_conditional_forecast(𝓂::ℳ, jacobian_state_func = 𝓂.solution.backward_looking.jacobian_state_func jacobian_shock_func = 𝓂.solution.backward_looking.jacobian_shock_func - # Get buffers and caches from backward_looking solution + # Get buffers from backward_looking solution jac_state_buffer = 𝓂.solution.backward_looking.jac_state_buffer jac_shock_buffer = 𝓂.solution.backward_looking.jac_shock_buffer - dydt_buffer = 𝓂.solution.backward_looking.dydt_buffer - lu_jac_state_cache = 𝓂.solution.backward_looking.lu_jac_state_cache # Get SS and parameter info for evaluating jacobians SS_and_pars = 𝓂.solution.non_stochastic_steady_state @@ -1104,26 +1102,19 @@ function get_conditional_forecast(𝓂::ℳ, n_cond = length(cond_idx) n_free = length(free_idx) - # Compute J = -(J_y_full \ J_e_full)[cond_idx_sorted, free_idx_sorted] using LinearSolve cache - # Solve J_y_full * X = J_e_full for X (store in dydt_buffer) - lu_jac_state_cache.A = jac_state_buffer - lu_jac_state_cache.b = jac_shock_buffer - 𝒮.solve!(lu_jac_state_cache) - dydt_buffer .= -lu_jac_state_cache.u + # Compute J = -(J_y_full \ J_e_full)[cond_idx_sorted, free_idx_sorted] using standard \ operator + # Solve J_y_full * X = J_e_full for X + dydt = -(jac_state_buffer \ jac_shock_buffer) - J = dydt_buffer[cond_idx_sorted, free_idx_sorted] + J = dydt[cond_idx_sorted, free_idx_sorted] if n_cond == n_free # Square system - use direct linear solve - J_prob = 𝒮.LinearProblem(J, r) - J_sol = 𝒮.solve(J_prob, 𝒮.LUFactorization()) - Δe = J_sol.u + Δe = J \ r else # Underdetermined - minimum norm solution: Δe = J' * (J * J')^{-1} * r JJt = J * J' - JJt_prob = 𝒮.LinearProblem(JJt, r) - JJt_sol = 𝒮.solve(JJt_prob, 𝒮.LUFactorization()) - Δe = J' * JJt_sol.u + Δe = J' * (JJt \ r) end return Δe diff --git a/src/macros.jl b/src/macros.jl index 0a8ec1036..35822b9b6 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -974,7 +974,7 @@ macro model(𝓂,ex...) second_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)), third_order_auxiliary_matrices(SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),Dict{Vector{Int}, Int}(),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0),SparseMatrixCSC{Int, Int64}(ℒ.I,0,0)) ), - backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization()), zeros(0,0), zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization())), + backward_looking_solution((x,y)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, (x,y,z)->nothing, Float64[], zeros(0,0), zeros(0,0), 𝒮.init(𝒮.LinearProblem(zeros(1,1), zeros(1)), 𝒮.LUFactorization()), zeros(0,0), zeros(0,0)), Float64[], # Set([:first_order]), Set(all_available_algorithms), diff --git a/src/structures.jl b/src/structures.jl index cc60fca91..c4f4de8d8 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -252,8 +252,6 @@ struct backward_looking_solution # Buffers for conditional forecasting jac_state_buffer::Matrix{Float64} # For J_y_full in conditional forecasting jac_shock_buffer::Matrix{Float64} # For J_e_full in conditional forecasting - dydt_buffer::Matrix{Float64} # For J_y_full \ J_e_full result - lu_jac_state_cache::𝒮.LinearCache # Cache for J_y_full \ b operations end mutable struct solution From 516c8d931115651448db313438cdf4aa3c6319a7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 17:19:06 +0000 Subject: [PATCH 27/40] Add baseline_path support for backward looking models with initial_state Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 36 ++++++++++++++++++++++++++++-------- src/get_functions.jl | 24 +++++++++++++++++++++++- 2 files changed, 51 insertions(+), 9 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 0ea014cb9..5d9d59681 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -8699,7 +8699,8 @@ function compute_irf_responses(𝓂::ℳ, generalised_irf_warmup_iterations::Int, generalised_irf_draws::Int, enforce_obc::Bool, - algorithm::Symbol) + algorithm::Symbol, + baseline_path::Union{Nothing, Matrix{Float64}} = nothing) if enforce_obc function obc_state_update(present_states, present_shocks::Vector{R}, state_update::Function) where R <: Float64 @@ -8783,7 +8784,8 @@ function compute_irf_responses(𝓂::ℳ, variables = variables, negative_shock = negative_shock, warmup_periods = generalised_irf_warmup_iterations, - draws = generalised_irf_draws) + draws = generalised_irf_draws, + baseline_path = baseline_path) else return irf(state_update, initial_state, @@ -8793,7 +8795,8 @@ function compute_irf_responses(𝓂::ℳ, shocks = shocks, shock_size = shock_size, variables = variables, - negative_shock = negative_shock) + negative_shock = negative_shock, + baseline_path = baseline_path) end end end @@ -8949,7 +8952,8 @@ function irf(state_update::Function, shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = :all, variables::Union{Symbol_input,String_input} = :all, shock_size::Real = 1, - negative_shock::Bool = false)::Union{KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{String}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{String}}}} + negative_shock::Bool = false, + baseline_path::Union{Nothing, Matrix{Float64}} = nothing)::Union{KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{String}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{String}}}} pruning = initial_state isa Vector{Vector{Float64}} @@ -9025,7 +9029,12 @@ function irf(state_update::Function, Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state end - return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) + # Use baseline_path if provided, otherwise use constant level + if baseline_path !== nothing + return KeyedArray(Y[var_idx,:,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) + else + return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) + end else Y = zeros(T.nVars,periods,length(shock_idx)) @@ -9059,7 +9068,12 @@ function irf(state_update::Function, axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] end - return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + # Use baseline_path if provided, otherwise use constant level + if baseline_path !== nothing + return KeyedArray(Y[var_idx,:,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + else + return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + end end end @@ -9075,7 +9089,8 @@ function girf(state_update::Function, shock_size::Real = 1, negative_shock::Bool = false, warmup_periods::Int = 100, - draws::Int = 50)::Union{KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{String}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{String}}}} + draws::Int = 50, + baseline_path::Union{Nothing, Matrix{Float64}} = nothing)::Union{KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{String}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{String},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{Symbol}}}, KeyedArray{Float64, 3, NamedDimsArray{(:Variables, :Periods, :Shocks), Float64, 3, Array{Float64, 3}}, Tuple{Vector{Symbol},UnitRange{Int},Vector{String}}}} pruning = initial_state isa Vector{Vector{Float64}} @@ -9241,7 +9256,12 @@ function girf(state_update::Function, axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] end - return KeyedArray(Y[var_idx,2:end,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + # Use baseline_path if provided, otherwise use constant level + if baseline_path !== nothing + return KeyedArray(Y[var_idx,2:end,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + else + return KeyedArray(Y[var_idx,2:end,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + end end diff --git a/src/get_functions.jl b/src/get_functions.jl index 1a9c52aa4..a131f3f1d 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1558,6 +1558,27 @@ function get_irf(𝓂::ℳ; end level = levels ? reference_steady_state + SSS_delta : SSS_delta + + # For backward looking models with initial_state provided, compute baseline path + # The baseline path is the no-shock forward iteration from initial_state + # This will be used for computing deviations instead of a constant SS level + baseline_path = nothing + if is_backward_looking && !unspecified_initial_state && algorithm == :newton + # Compute the no-shock baseline path in deviations + nVars = 𝓂.timings.nVars + baseline_path = zeros(nVars, periods) + zero_shocks = zeros(length(𝓂.timings.exo)) + + # Start from the initial_state (already converted to deviations above) + baseline_state = copy(initial_state) + for t in 1:periods + baseline_state = state_update(baseline_state, zero_shocks) + baseline_path[:, t] = baseline_state + end + + # Add the level offset (SSS_delta for deviations mode, full level for levels mode) + baseline_path .+= level + end responses = compute_irf_responses(𝓂, state_update, @@ -1572,7 +1593,8 @@ function get_irf(𝓂::ℳ; generalised_irf_warmup_iterations = generalised_irf_warmup_iterations, generalised_irf_draws = generalised_irf_draws, enforce_obc = occasionally_binding_constraints, - algorithm = algorithm) + algorithm = algorithm, + baseline_path = baseline_path) return responses From 129368119b157f7504ad62a3efcf7cd443bafab9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 17:28:19 +0000 Subject: [PATCH 28/40] Fix baseline path deviations for explosive models and add comprehensive tests Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- src/MacroModelling.jl | 15 +++++---- src/get_functions.jl | 13 ++++---- test/test_backward_looking_models.jl | 48 ++++++++++++++++++++++++++++ 3 files changed, 64 insertions(+), 12 deletions(-) diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 5d9d59681..348707b31 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -9029,9 +9029,10 @@ function irf(state_update::Function, Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state end - # Use baseline_path if provided, otherwise use constant level + # Use baseline_path if provided: return (Y - baseline_path + level) = shock effect + level offset + # baseline_path is in deviations, so Y - baseline_path gives pure shock effect if baseline_path !== nothing - return KeyedArray(Y[var_idx,:,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) + return KeyedArray(Y[var_idx,:,:] .- baseline_path[var_idx,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) else return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none]) end @@ -9068,9 +9069,10 @@ function irf(state_update::Function, axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] end - # Use baseline_path if provided, otherwise use constant level + # Use baseline_path if provided: return (Y - baseline_path + level) = shock effect + level offset + # baseline_path is in deviations, so Y - baseline_path gives pure shock effect if baseline_path !== nothing - return KeyedArray(Y[var_idx,:,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + return KeyedArray(Y[var_idx,:,:] .- baseline_path[var_idx,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) else return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) end @@ -9256,9 +9258,10 @@ function girf(state_update::Function, axis2 = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in axis2_decomposed] end - # Use baseline_path if provided, otherwise use constant level + # Use baseline_path if provided: return (Y - baseline_path + level) = shock effect + level offset + # baseline_path is in deviations, so Y - baseline_path gives pure shock effect if baseline_path !== nothing - return KeyedArray(Y[var_idx,2:end,:] .+ baseline_path[var_idx,:]; Variables = axis1, Periods = 1:periods, Shocks = axis2) + return KeyedArray(Y[var_idx,2:end,:] .- baseline_path[var_idx,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) else return KeyedArray(Y[var_idx,2:end,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = axis2) end diff --git a/src/get_functions.jl b/src/get_functions.jl index a131f3f1d..e468580b1 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1561,10 +1561,11 @@ function get_irf(𝓂::ℳ; # For backward looking models with initial_state provided, compute baseline path # The baseline path is the no-shock forward iteration from initial_state - # This will be used for computing deviations instead of a constant SS level + # When levels=false: IRFs will show deviations from this baseline (shock effect) + # When levels=true: baseline_path is not used (return simulation in levels) baseline_path = nothing - if is_backward_looking && !unspecified_initial_state && algorithm == :newton - # Compute the no-shock baseline path in deviations + if is_backward_looking && !unspecified_initial_state && algorithm == :newton && !levels + # Compute the no-shock baseline path in deviations from NSSS nVars = 𝓂.timings.nVars baseline_path = zeros(nVars, periods) zero_shocks = zeros(length(𝓂.timings.exo)) @@ -1575,9 +1576,9 @@ function get_irf(𝓂::ℳ; baseline_state = state_update(baseline_state, zero_shocks) baseline_path[:, t] = baseline_state end - - # Add the level offset (SSS_delta for deviations mode, full level for levels mode) - baseline_path .+= level + # baseline_path is now in deviations from NSSS + # In irf(), we compute: Y - baseline_path + level + # which gives: simulation_deviations - baseline_deviations + SSS_delta = shock_effect end responses = compute_irf_responses(𝓂, diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 90ed46b1f..355366244 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -200,6 +200,54 @@ @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 @test y_values[2] > y_values[1] # Continues growing + # Test the key identity: levels=false deviation + baseline_levels = levels=true + # This verifies that baseline path deviations work correctly for explosive models + + # Get IRF with shocks in deviations mode (levels=false is default) + irf_dev = get_irf(UnstableButValidSS, + algorithm = :newton, + initial_state = initial_state_levels, + periods = 10) + + # Get IRF in levels mode + irf_lev = get_irf(UnstableButValidSS, + algorithm = :newton, + initial_state = initial_state_levels, + levels = true, + periods = 10) + + # Get baseline path (no shocks) in levels + baseline_lev = get_irf(UnstableButValidSS, + algorithm = :newton, + initial_state = initial_state_levels, + shocks = :none, + levels = true, + periods = 10) + + # Identity: irf_dev + baseline_lev ≈ irf_lev + # (deviation from baseline + baseline in levels = simulation in levels) + @test isapprox(collect(irf_dev(:y, :, :eps_z)) .+ collect(baseline_lev(:y, :, :none)), + collect(irf_lev(:y, :, :eps_z)), rtol=1e-10) + @test isapprox(collect(irf_dev(:z, :, :eps_z)) .+ collect(baseline_lev(:z, :, :none)), + collect(irf_lev(:z, :, :eps_z)), rtol=1e-10) + + # Test conditional forecasting with non-trivial starting point + conditions = KeyedArray(Matrix{Union{Nothing,Float64}}(nothing, 2, 2), + Variables = UnstableButValidSS.var, + Periods = 1:2) + conditions[:y, 1] = initial_y * 1.1 # Condition y to be 10% above initial in period 1 + conditions[:z, 2] = 0.05 # Condition z in period 2 + + cf_result = get_conditional_forecast(UnstableButValidSS, + conditions, + algorithm = :newton, + initial_state = initial_state_levels, + periods = 5) + + # Check that conditions are satisfied + @test isapprox(cf_result(:y, 1), initial_y * 1.1, rtol=1e-4) + @test isapprox(cf_result(:z, 2), 0.05, rtol=1e-4) + UnstableButValidSS = nothing end From 165b068b3f3da24c6edabd1a102e6cd3d97e506d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 18:09:28 +0000 Subject: [PATCH 29/40] Fix backward looking model tests: proper 3D indexing, correct expectations for y response timing Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- test/test_backward_looking_models.jl | 49 ++++++++++++++++------------ 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 355366244..9f63e6163 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -130,26 +130,30 @@ initial_state_solow = [initial_k, 0.0] # [k, z] # For backward looking models with newton algorithm, initial_state is in levels + # shocks = :none returns result with 3rd dimension as [:none] result = get_irf(SolowGrowth2, algorithm = :newton, initial_state = initial_state_solow, shocks = :none, + levels = true, # Request levels output periods = 10) @test size(result, 1) == 2 # 2 variables @test size(result, 2) == 10 # 10 periods # Capital should increase over time (converging toward steady state) - @test result(:k, 1) > initial_k # First period k should increase - @test result(:k, 10) > result(:k, 1) # Later periods should be higher + # Access with 3 indices since result is 3D (Vars, Periods, Shocks) + @test result(:k, 1, :none) > initial_k # First period k should increase + @test result(:k, 10, :none) > result(:k, 1, :none) # Later periods should be higher SolowGrowth2 = nothing end @testset "Model with unstable eigenvalue but valid SS at zero" begin - # Model y[0] = (1 + g) * y[-1] with g > 0 + # Model y[0] = (1 + g) * y[-1] + a * z[-1] with g > 0 # Has eigenvalue > 1 (unstable) but valid SS at y = 0 - # From SS, iterating stays at SS. From any other point, it explodes. + # Note: y depends on z[-1], not z[0], so y response to eps_z shock + # appears only in period 2+ (delayed by one period) @model UnstableButValidSS begin y[0] = (1 + g) * y[-1] + a * z[-1] z[0] = rho * z[-1] + sigma * eps_z[x] @@ -177,8 +181,9 @@ @test size(irf_newton, 3) == 1 # 1 shock # IRF from SS should show the effect of a shock - # y responds to z shock through the 'a' parameter - @test irf_newton(:y, 1, :eps_z) != 0.0 # y responds to z shock + # y depends on z[-1], so y responds to z shock in period 2 (delayed by one period) + @test irf_newton(:y, 2, :eps_z) != 0.0 # y responds to z shock in period 2 + @test irf_newton(:y, 1, :eps_z) == 0.0 # y has no immediate response (depends on z[-1]) # Can also provide initial_state in levels for simulation from non-SS point initial_y = 100.0 # Start at y = 100 @@ -188,14 +193,14 @@ algorithm = :newton, initial_state = initial_state_levels, shocks = :none, + levels = true, # Request levels output periods = 10) @test size(result, 1) == 2 # 2 variables @test size(result, 2) == 10 # 10 periods # y should grow each period at rate g (with z = 0) - # Since initial_state is in levels and model has valid SS, - # result is in levels (initial_state_levels provided) + # With levels = true, result is in levels y_values = result(:y, :, :none) @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 @test y_values[2] > y_values[1] # Continues growing @@ -232,11 +237,14 @@ collect(irf_lev(:z, :, :eps_z)), rtol=1e-10) # Test conditional forecasting with non-trivial starting point - conditions = KeyedArray(Matrix{Union{Nothing,Float64}}(nothing, 2, 2), - Variables = UnstableButValidSS.var, - Periods = 1:2) - conditions[:y, 1] = initial_y * 1.1 # Condition y to be 10% above initial in period 1 - conditions[:z, 2] = 0.05 # Condition z in period 2 + # Note: y[1] = (1+g)*y[-1] + a*z[-1] depends only on predetermined values, + # so we can only condition z (which has shock eps_z) in period 1. + # For period 2, y depends on z[1] which is affected by eps_z[1], so both can be conditioned. + conditions = Matrix{Union{Nothing,Float64}}(nothing, 2, 2) + z_idx = findfirst(x -> x == :z, UnstableButValidSS.var) + y_idx = findfirst(x -> x == :y, UnstableButValidSS.var) + conditions[z_idx, 1] = 0.05 # Condition z in period 1 (affected by eps_z shock) + conditions[z_idx, 2] = 0.04 # Condition z in period 2 cf_result = get_conditional_forecast(UnstableButValidSS, conditions, @@ -245,8 +253,8 @@ periods = 5) # Check that conditions are satisfied - @test isapprox(cf_result(:y, 1), initial_y * 1.1, rtol=1e-4) - @test isapprox(cf_result(:z, 2), 0.05, rtol=1e-4) + @test isapprox(cf_result(:z, 1), 0.05, rtol=1e-4) + @test isapprox(cf_result(:z, 2), 0.04, rtol=1e-4) UnstableButValidSS = nothing end @@ -319,11 +327,12 @@ @test VAR_cond.timings.nFuture_not_past_and_mixed == 0 # Test conditional forecast with newton - condition y to be 0.05 in period 1 - conditions = KeyedArray(Matrix{Union{Nothing,Float64}}(nothing, 2, 2), - Variables = VAR_cond.var, - Periods = 1:2) - conditions[:y, 1] = 0.05 # Condition y in period 1 - conditions[:x, 2] = 0.02 # Condition x in period 2 + # Use Matrix directly instead of KeyedArray to avoid indexing issues + conditions = Matrix{Union{Nothing,Float64}}(nothing, 2, 2) + y_idx = findfirst(x -> x == :y, VAR_cond.var) + x_idx = findfirst(x -> x == :x, VAR_cond.var) + conditions[y_idx, 1] = 0.05 # Condition y in period 1 + conditions[x_idx, 2] = 0.02 # Condition x in period 2 # Use newton algorithm for conditional forecast cf_newton = get_conditional_forecast(VAR_cond, conditions, algorithm = :newton) From 769bf9a7bca66e9606a3103fb217e0dd0c09e721 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 20:32:41 +0000 Subject: [PATCH 30/40] Add deviations_from parameter, default algorithm=newton for backward looking models, add tests Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 6 ++- src/common_docstrings.jl | 3 +- src/default_options.jl | 5 +++ src/get_functions.jl | 19 +++++--- test/test_backward_looking_models.jl | 66 ++++++++++++++++++---------- 5 files changed, 67 insertions(+), 32 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 9c60d34c1..c5de88e76 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -3,7 +3,7 @@ module StatsPlotsExt using MacroModelling import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, merge_calculation_options, MODEL®, DATA®, PARAMETERS®, ALGORITHM®, FILTER®, VARIABLES®, SMOOTH®, SHOW_PLOTS®, SAVE_PLOTS®, SAVE_PLOTS_NAME®, SAVE_PLOTS_FORMAT®, SAVE_PLOTS_PATH®, PLOTS_PER_PAGE®, MAX_ELEMENTS_PER_LEGENDS_ROW®, EXTRA_LEGEND_SPACE®, PLOT_ATTRIBUTES®, QME®, SYLVESTER®, LYAPUNOV®, TOLERANCES®, VERBOSE®, DATA_IN_LEVELS®, PERIODS®, SHOCKS®, SHOCK_SIZE®, NEGATIVE_SHOCK®, GENERALISED_IRF®, GENERALISED_IRF_WARMUP_ITERATIONS®, CONDITIONS_IN_LEVELS®, GENERALISED_IRF_DRAWS®, INITIAL_STATE®, IGNORE_OBC®, CONDITIONS®, SHOCK_CONDITIONS®, LEVELS®, LABEL®, RENAME_DICTIONARY®, parse_shocks_input_to_index, parse_variables_input_to_index, replace_indices, replace_indices_special, filter_data_with_model, get_relevant_steady_states, replace_indices_in_symbol, parse_algorithm_to_state_update, girf, decompose_name, obc_objective_optim_fun, obc_constraint_optim_fun, compute_irf_responses, process_ignore_obc_flag, adjust_generalised_irf_flag, process_shocks_input, normalize_filtering_options, infer_step -import MacroModelling: DEFAULT_ALGORITHM, DEFAULT_FILTER_SELECTOR, DEFAULT_WARMUP_ITERATIONS, DEFAULT_VARIABLES_EXCLUDING_OBC, DEFAULT_SHOCK_SELECTION, DEFAULT_PRESAMPLE_PERIODS, DEFAULT_DATA_IN_LEVELS, DEFAULT_SHOCK_DECOMPOSITION_SELECTOR, DEFAULT_SMOOTH_SELECTOR, DEFAULT_LABEL, DEFAULT_SHOW_PLOTS, DEFAULT_SAVE_PLOTS, DEFAULT_SAVE_PLOTS_FORMAT, DEFAULT_SAVE_PLOTS_PATH, DEFAULT_PLOTS_PER_PAGE_SMALL, DEFAULT_TRANSPARENCY, DEFAULT_MAX_ELEMENTS_PER_LEGEND_ROW, DEFAULT_EXTRA_LEGEND_SPACE, DEFAULT_VERBOSE, DEFAULT_QME_ALGORITHM, DEFAULT_SYLVESTER_SELECTOR, DEFAULT_SYLVESTER_THRESHOLD, DEFAULT_LARGE_SYLVESTER_ALGORITHM, DEFAULT_SYLVESTER_ALGORITHM, DEFAULT_LYAPUNOV_ALGORITHM, DEFAULT_PLOT_ATTRIBUTES, DEFAULT_ARGS_AND_KWARGS_NAMES, DEFAULT_PLOTS_PER_PAGE_LARGE, DEFAULT_SHOCKS_EXCLUDING_OBC, DEFAULT_VARIABLES_EXCLUDING_AUX_AND_OBC, DEFAULT_PERIODS, DEFAULT_SHOCK_SIZE, DEFAULT_NEGATIVE_SHOCK, DEFAULT_GENERALISED_IRF, DEFAULT_GENERALISED_IRF_WARMUP, DEFAULT_GENERALISED_IRF_DRAWS, DEFAULT_INITIAL_STATE, DEFAULT_IGNORE_OBC, DEFAULT_PLOT_TYPE, DEFAULT_CONDITIONS_IN_LEVELS, DEFAULT_SIGMA_RANGE, DEFAULT_FONT_SIZE, DEFAULT_VARIABLE_SELECTION, DEFAULT_FORECAST_PERIODS +import MacroModelling: DEFAULT_ALGORITHM, DEFAULT_FILTER_SELECTOR, DEFAULT_WARMUP_ITERATIONS, DEFAULT_VARIABLES_EXCLUDING_OBC, DEFAULT_SHOCK_SELECTION, DEFAULT_PRESAMPLE_PERIODS, DEFAULT_DATA_IN_LEVELS, DEFAULT_SHOCK_DECOMPOSITION_SELECTOR, DEFAULT_SMOOTH_SELECTOR, DEFAULT_LABEL, DEFAULT_SHOW_PLOTS, DEFAULT_SAVE_PLOTS, DEFAULT_SAVE_PLOTS_FORMAT, DEFAULT_SAVE_PLOTS_PATH, DEFAULT_PLOTS_PER_PAGE_SMALL, DEFAULT_TRANSPARENCY, DEFAULT_MAX_ELEMENTS_PER_LEGEND_ROW, DEFAULT_EXTRA_LEGEND_SPACE, DEFAULT_VERBOSE, DEFAULT_QME_ALGORITHM, DEFAULT_SYLVESTER_SELECTOR, DEFAULT_SYLVESTER_THRESHOLD, DEFAULT_LARGE_SYLVESTER_ALGORITHM, DEFAULT_SYLVESTER_ALGORITHM, DEFAULT_LYAPUNOV_ALGORITHM, DEFAULT_PLOT_ATTRIBUTES, DEFAULT_ARGS_AND_KWARGS_NAMES, DEFAULT_PLOTS_PER_PAGE_LARGE, DEFAULT_SHOCKS_EXCLUDING_OBC, DEFAULT_VARIABLES_EXCLUDING_AUX_AND_OBC, DEFAULT_PERIODS, DEFAULT_SHOCK_SIZE, DEFAULT_NEGATIVE_SHOCK, DEFAULT_GENERALISED_IRF, DEFAULT_GENERALISED_IRF_WARMUP, DEFAULT_GENERALISED_IRF_DRAWS, DEFAULT_INITIAL_STATE, DEFAULT_IGNORE_OBC, DEFAULT_PLOT_TYPE, DEFAULT_CONDITIONS_IN_LEVELS, DEFAULT_SIGMA_RANGE, DEFAULT_FONT_SIZE, DEFAULT_VARIABLE_SELECTION, DEFAULT_FORECAST_PERIODS, DEFAULT_ALGORITHM_BACKWARD_LOOKING, DEFAULT_DEVIATIONS_FROM import DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF, TYPEDSIGNATURES, TYPEDFIELDS import LaTeXStrings @@ -1802,7 +1802,7 @@ function plot_irf(𝓂::ℳ; save_plots_name::Union{String, Symbol} = "irf", save_plots_path::String = DEFAULT_SAVE_PLOTS_PATH, plots_per_page::Int = DEFAULT_PLOTS_PER_PAGE_LARGE, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, shock_size::Real = DEFAULT_SHOCK_SIZE, negative_shock::Bool = DEFAULT_NEGATIVE_SHOCK, generalised_irf::Bool = DEFAULT_GENERALISED_IRF, @@ -1810,6 +1810,8 @@ function plot_irf(𝓂::ℳ; generalised_irf_draws::Int = DEFAULT_GENERALISED_IRF_DRAWS, initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = DEFAULT_INITIAL_STATE, ignore_obc::Bool = DEFAULT_IGNORE_OBC, + plot_baseline::Bool = 𝓂.timings.nFuture_not_past_and_mixed == 0, + deviations_from::Symbol = DEFAULT_DEVIATIONS_FROM, rename_dictionary::AbstractDict{<:Union{Symbol, String}, <:Union{Symbol, String}} = Dict{Symbol, String}(), plot_attributes::Dict = Dict(), verbose::Bool = DEFAULT_VERBOSE, diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl index af6be3e33..0d464b89a 100644 --- a/src/common_docstrings.jl +++ b/src/common_docstrings.jl @@ -42,4 +42,5 @@ const INITIAL_STATE® = "`initial_state` [Default: `$(DEFAULT_INITIAL_STATE)`, T const INITIAL_STATE®1 = "`initial_state` [Default: `$(DEFAULT_INITIAL_STATE)`, Type: `Vector{Float64}`]: The initial state defines the starting point for the model (in levels, not deviations). The state includes all variables as well as exogenous variables in leads or lags if present. `get_irf(𝓂, shocks = :none, variables = :all, periods = 1)` returns a `KeyedArray` with all variables. The `KeyedArray` type is provided by the `AxisKeys` package." const LABEL® = "`label` [Type: `Union{Real, String, Symbol}`]: label to attribute to this function call in the plots. The default is the number of previous function calls since the last call to the function version with ! + 1." const RENAME_DICTIONARY® = "`rename_dictionary` [Default: `Dict()`, Type: `Dict{Symbol, String}`]: dictionary mapping variable or shock symbols to custom display names in plots. For example: `Dict(:dinve => \"Investment growth\", :c => \"Consumption\")`. Variables/shocks not in the dictionary will use their default names." -const CONDITIONS_IN_LEVELS® = "`conditions_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the conditions are provided in levels. If `true` the input to the conditions argument will have the relevant steady state subtracted (non-stochastic or stochastic steady state depending on the solution algorithm)." \ No newline at end of file +const CONDITIONS_IN_LEVELS® = "`conditions_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the conditions are provided in levels. If `true` the input to the conditions argument will have the relevant steady state subtracted (non-stochastic or stochastic steady state depending on the solution algorithm)." +const DEVIATIONS_FROM® = "`deviations_from` [Default: `$(DEFAULT_DEVIATIONS_FROM)`, Type: `Symbol`]: reference point for computing deviations when `levels = false`. Options are `:steady_state` (deviations from the relevant steady state, the default behaviour) or `:baseline` (deviations from the no-shock path starting from `initial_state` - useful for backward looking models with explosive dynamics). Ignored when `levels = true`." \ No newline at end of file diff --git a/src/default_options.jl b/src/default_options.jl index 16473aa82..acce1c809 100644 --- a/src/default_options.jl +++ b/src/default_options.jl @@ -3,6 +3,11 @@ # General algorithm and filtering defaults const DEFAULT_ALGORITHM = :first_order const DEFAULT_ALGORITHM_SELECTOR = stochastic -> stochastic ? :second_order : :first_order +# For backward looking models, default to newton algorithm +const DEFAULT_ALGORITHM_BACKWARD_LOOKING = :newton + +# Deviations mode for IRF output +const DEFAULT_DEVIATIONS_FROM = :steady_state # :steady_state or :baseline const DEFAULT_FILTER_SELECTOR = algorithm -> algorithm == :first_order ? :kalman : :inversion const DEFAULT_SHOCK_DECOMPOSITION_SELECTOR = algorithm -> algorithm ∉ (:second_order, :third_order) const DEFAULT_SMOOTH_SELECTOR = filter -> filter == :kalman diff --git a/src/get_functions.jl b/src/get_functions.jl index e468580b1..37c5f2a52 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1439,10 +1439,11 @@ And data, 4×40×1 Array{Float64, 3}: """ function get_irf(𝓂::ℳ; periods::Int = DEFAULT_PERIODS, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, parameters::ParameterType = nothing, variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = DEFAULT_SHOCKS_EXCLUDING_OBC, + deviations_from::Symbol = DEFAULT_DEVIATIONS_FROM, negative_shock::Bool = DEFAULT_NEGATIVE_SHOCK, generalised_irf::Bool = DEFAULT_GENERALISED_IRF, generalised_irf_warmup_iterations::Int = DEFAULT_GENERALISED_IRF_WARMUP, @@ -1514,6 +1515,14 @@ function get_irf(𝓂::ℳ; @assert !unspecified_initial_state "Model is backward looking with no valid steady state. Provide initial_state in levels." end + + # Validate deviations_from parameter + @assert deviations_from ∈ [:steady_state, :baseline] "deviations_from must be either :steady_state or :baseline" + + # If levels = true, ignore deviations_from (levels are levels) + if levels && deviations_from == :baseline + @warn "deviations_from = :baseline is ignored when levels = true (returning levels)" + end # @timeit_debug timer "Get relevant steady state" begin @@ -1559,12 +1568,12 @@ function get_irf(𝓂::ℳ; level = levels ? reference_steady_state + SSS_delta : SSS_delta - # For backward looking models with initial_state provided, compute baseline path + # For backward looking models with deviations_from = :baseline, compute baseline path # The baseline path is the no-shock forward iteration from initial_state - # When levels=false: IRFs will show deviations from this baseline (shock effect) - # When levels=true: baseline_path is not used (return simulation in levels) + # When deviations_from = :baseline and levels = false: IRFs will show deviations from this baseline (shock effect) + # When levels = true: baseline_path is not used (return simulation in levels) baseline_path = nothing - if is_backward_looking && !unspecified_initial_state && algorithm == :newton && !levels + if !levels && deviations_from == :baseline && is_backward_looking && algorithm == :newton # Compute the no-shock baseline path in deviations from NSSS nVars = 𝓂.timings.nVars baseline_path = zeros(nVars, periods) diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 9f63e6163..b0c1781e4 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -31,7 +31,12 @@ sol = get_solution(VAR2) @test !isnothing(sol) - # Test IRFs with first_order algorithm + # Test IRFs - default algorithm should be :newton for backward looking models + irf_default = get_irf(VAR2) + @test size(irf_default, 1) == 2 # 2 variables + @test size(irf_default, 3) == 2 # 2 shocks + + # Test IRFs with first_order algorithm (explicitly specified) irf_first = get_irf(VAR2, algorithm = :first_order) @test size(irf_first, 1) == 2 # 2 variables @test size(irf_first, 3) == 2 # 2 shocks @@ -40,13 +45,8 @@ sim = simulate(VAR2) @test size(sim, 1) == 2 # 2 variables - # Test newton algorithm - should produce same results as first_order for linear model - irf_newton = get_irf(VAR2, algorithm = :newton) - @test size(irf_newton, 1) == 2 # 2 variables - @test size(irf_newton, 3) == 2 # 2 shocks - # For linear models, first_order and newton should give approximately the same results - @test isapprox(collect(irf_first), collect(irf_newton), rtol=1e-6) + @test isapprox(collect(irf_first), collect(irf_default), rtol=1e-6) VAR2 = nothing end @@ -85,7 +85,7 @@ sol = get_solution(SolowGrowth) @test !isnothing(sol) - # Test IRFs with first_order algorithm + # Test IRFs - default algorithm for backward looking is :newton irf = get_irf(SolowGrowth) @test size(irf, 1) == 2 # 2 variables (k, z) @test size(irf, 3) == 1 # 1 shock @@ -94,15 +94,29 @@ sim = simulate(SolowGrowth) @test size(sim, 1) == 2 # 2 variables - # Test newton algorithm for nonlinear backward looking model - irf_newton = get_irf(SolowGrowth, algorithm = :newton) - @test size(irf_newton, 1) == 2 # 2 variables (k, z) - @test size(irf_newton, 3) == 1 # 1 shock + # Test first_order algorithm for comparison + irf_first = get_irf(SolowGrowth, algorithm = :first_order) + @test size(irf_first, 1) == 2 # 2 variables (k, z) + @test size(irf_first, 3) == 1 # 1 shock # Test simulation with newton sim_newton = simulate(SolowGrowth, algorithm = :newton) @test size(sim_newton, 1) == 2 # 2 variables + # Test conditional forecasting with newton + # Condition z to be 0.02 in period 1 + conditions = Matrix{Union{Nothing,Float64}}(nothing, 2, 2) + k_idx = findfirst(x -> x == :k, SolowGrowth.var) + z_idx = findfirst(x -> x == :z, SolowGrowth.var) + conditions[z_idx, 1] = 0.02 # Condition z in period 1 + conditions[z_idx, 2] = 0.01 # Condition z in period 2 + + cf_newton = get_conditional_forecast(SolowGrowth, conditions, algorithm = :newton) + + # Check that conditions are met + @test isapprox(cf_newton(:z, 1), 0.02, rtol=1e-4) + @test isapprox(cf_newton(:z, 2), 0.01, rtol=1e-4) + SolowGrowth = nothing end @@ -146,6 +160,9 @@ @test result(:k, 1, :none) > initial_k # First period k should increase @test result(:k, 10, :none) > result(:k, 1, :none) # Later periods should be higher + # Test plot_irf works (just check it doesn't error) + # plot_irf(SolowGrowth2, algorithm = :newton, show_plots = false) + SolowGrowth2 = nothing end @@ -176,7 +193,8 @@ # Since SS is valid, newton works WITHOUT requiring initial_state # (starts from SS=0 and computes IRFs in deviations) - irf_newton = get_irf(UnstableButValidSS, algorithm = :newton) + # Default algorithm for backward looking models is now :newton + irf_newton = get_irf(UnstableButValidSS) @test size(irf_newton, 1) == 2 # 2 variables @test size(irf_newton, 3) == 1 # 1 shock @@ -205,13 +223,12 @@ @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 @test y_values[2] > y_values[1] # Continues growing - # Test the key identity: levels=false deviation + baseline_levels = levels=true - # This verifies that baseline path deviations work correctly for explosive models - - # Get IRF with shocks in deviations mode (levels=false is default) - irf_dev = get_irf(UnstableButValidSS, + # Test deviations_from = :baseline + # Get IRF with shocks in deviations from baseline mode + irf_baseline = get_irf(UnstableButValidSS, algorithm = :newton, initial_state = initial_state_levels, + deviations_from = :baseline, periods = 10) # Get IRF in levels mode @@ -229,11 +246,11 @@ levels = true, periods = 10) - # Identity: irf_dev + baseline_lev ≈ irf_lev + # Identity: irf_baseline + baseline_lev ≈ irf_lev # (deviation from baseline + baseline in levels = simulation in levels) - @test isapprox(collect(irf_dev(:y, :, :eps_z)) .+ collect(baseline_lev(:y, :, :none)), + @test isapprox(collect(irf_baseline(:y, :, :eps_z)) .+ collect(baseline_lev(:y, :, :none)), collect(irf_lev(:y, :, :eps_z)), rtol=1e-10) - @test isapprox(collect(irf_dev(:z, :, :eps_z)) .+ collect(baseline_lev(:z, :, :none)), + @test isapprox(collect(irf_baseline(:z, :, :eps_z)) .+ collect(baseline_lev(:z, :, :none)), collect(irf_lev(:z, :, :eps_z)), rtol=1e-10) # Test conditional forecasting with non-trivial starting point @@ -286,7 +303,8 @@ # For unit root models with valid SS, newton algorithm should work WITHOUT requiring initial_state # Deviations are computed against the no-shock reference path - irf_newton = get_irf(UnitRootAR1, algorithm = :newton) + # Default algorithm for backward looking is :newton + irf_newton = get_irf(UnitRootAR1) @test size(irf_newton, 1) == 2 # 2 variables @test size(irf_newton, 3) == 2 # 2 shocks @@ -334,8 +352,8 @@ conditions[y_idx, 1] = 0.05 # Condition y in period 1 conditions[x_idx, 2] = 0.02 # Condition x in period 2 - # Use newton algorithm for conditional forecast - cf_newton = get_conditional_forecast(VAR_cond, conditions, algorithm = :newton) + # Use newton algorithm (default for backward looking) for conditional forecast + cf_newton = get_conditional_forecast(VAR_cond, conditions) # Check that conditions are met @test isapprox(cf_newton(:y, 1), 0.05, rtol=1e-6) From bdd90fdc759dc4b3e91d8d6e3fd4dd8585204e9b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 20:38:44 +0000 Subject: [PATCH 31/40] Add docstrings for deviations_from and plot_baseline parameters Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 2 ++ src/get_functions.jl | 1 + 2 files changed, 3 insertions(+) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index c5de88e76..ee695c7a7 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -1750,6 +1750,8 @@ If the model contains occasionally binding constraints and `ignore_obc = false` - $GENERALISED_IRF_DRAWS® - $INITIAL_STATE® - $IGNORE_OBC® +- `plot_baseline` [Default: `true` for backward looking models, `false` otherwise, Type: `Bool`]: plot the baseline path (no-shock trajectory from `initial_state`) as a solid black line instead of the reference steady state. Useful for backward looking models with explosive dynamics. +- `deviations_from` [Default: `:steady_state`, Type: `Symbol`]: reference point for deviations when not in levels mode. Options: `:steady_state` (deviations from relevant steady state) or `:baseline` (deviations from no-shock path starting from `initial_state`). - `label` [Default: `1`, Type: `Union{Real, String, Symbol}`]: label to attribute to this function call in the plots. - $SHOW_PLOTS® - $SAVE_PLOTS® diff --git a/src/get_functions.jl b/src/get_functions.jl index 37c5f2a52..f8fad3c7a 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1386,6 +1386,7 @@ If the model contains occasionally binding constraints and `ignore_obc = false` - $PARAMETERS® - $(VARIABLES®(DEFAULT_VARIABLES_EXCLUDING_OBC)) - $SHOCKS® +- $DEVIATIONS_FROM® - $NEGATIVE_SHOCK® - $GENERALISED_IRF® - $GENERALISED_IRF_WARMUP_ITERATIONS® From 90ab1587eb629252226994166e38fe4cf0e7bc5e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 20:40:12 +0000 Subject: [PATCH 32/40] Remove commented-out plot_irf test Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- test/test_backward_looking_models.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index b0c1781e4..50b5af13d 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -160,9 +160,6 @@ @test result(:k, 1, :none) > initial_k # First period k should increase @test result(:k, 10, :none) > result(:k, 1, :none) # Later periods should be higher - # Test plot_irf works (just check it doesn't error) - # plot_irf(SolowGrowth2, algorithm = :newton, show_plots = false) - SolowGrowth2 = nothing end From 394b6f00f5b3b19351bec9e2ef59dd2296adc8d0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 20:48:26 +0000 Subject: [PATCH 33/40] Rename plot_baseline to reference parameter in plot_irf Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index ee695c7a7..baf5b7fd5 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -1750,7 +1750,7 @@ If the model contains occasionally binding constraints and `ignore_obc = false` - $GENERALISED_IRF_DRAWS® - $INITIAL_STATE® - $IGNORE_OBC® -- `plot_baseline` [Default: `true` for backward looking models, `false` otherwise, Type: `Bool`]: plot the baseline path (no-shock trajectory from `initial_state`) as a solid black line instead of the reference steady state. Useful for backward looking models with explosive dynamics. +- `reference` [Default: `:baseline` for backward looking models, `:steady_state` otherwise, Type: `Symbol`]: reference path for plotting. Options: `:steady_state` (horizontal line at steady state) or `:baseline` (no-shock trajectory from `initial_state`). Useful for backward looking models with explosive dynamics. - `deviations_from` [Default: `:steady_state`, Type: `Symbol`]: reference point for deviations when not in levels mode. Options: `:steady_state` (deviations from relevant steady state) or `:baseline` (deviations from no-shock path starting from `initial_state`). - `label` [Default: `1`, Type: `Union{Real, String, Symbol}`]: label to attribute to this function call in the plots. - $SHOW_PLOTS® @@ -1812,7 +1812,7 @@ function plot_irf(𝓂::ℳ; generalised_irf_draws::Int = DEFAULT_GENERALISED_IRF_DRAWS, initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = DEFAULT_INITIAL_STATE, ignore_obc::Bool = DEFAULT_IGNORE_OBC, - plot_baseline::Bool = 𝓂.timings.nFuture_not_past_and_mixed == 0, + reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, deviations_from::Symbol = DEFAULT_DEVIATIONS_FROM, rename_dictionary::AbstractDict{<:Union{Symbol, String}, <:Union{Symbol, String}} = Dict{Symbol, String}(), plot_attributes::Dict = Dict(), From d4a0adec3843514118eec8cffc52b3845be558c2 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 22:03:35 +0100 Subject: [PATCH 34/40] rename to reference --- ext/StatsPlotsExt.jl | 6 +++--- src/common_docstrings.jl | 2 +- src/default_options.jl | 2 +- src/get_functions.jl | 20 ++++++++++---------- test/test_backward_looking_models.jl | 20 ++++++++++---------- 5 files changed, 25 insertions(+), 25 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index baf5b7fd5..609e8cae8 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -3,7 +3,7 @@ module StatsPlotsExt using MacroModelling import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, merge_calculation_options, MODEL®, DATA®, PARAMETERS®, ALGORITHM®, FILTER®, VARIABLES®, SMOOTH®, SHOW_PLOTS®, SAVE_PLOTS®, SAVE_PLOTS_NAME®, SAVE_PLOTS_FORMAT®, SAVE_PLOTS_PATH®, PLOTS_PER_PAGE®, MAX_ELEMENTS_PER_LEGENDS_ROW®, EXTRA_LEGEND_SPACE®, PLOT_ATTRIBUTES®, QME®, SYLVESTER®, LYAPUNOV®, TOLERANCES®, VERBOSE®, DATA_IN_LEVELS®, PERIODS®, SHOCKS®, SHOCK_SIZE®, NEGATIVE_SHOCK®, GENERALISED_IRF®, GENERALISED_IRF_WARMUP_ITERATIONS®, CONDITIONS_IN_LEVELS®, GENERALISED_IRF_DRAWS®, INITIAL_STATE®, IGNORE_OBC®, CONDITIONS®, SHOCK_CONDITIONS®, LEVELS®, LABEL®, RENAME_DICTIONARY®, parse_shocks_input_to_index, parse_variables_input_to_index, replace_indices, replace_indices_special, filter_data_with_model, get_relevant_steady_states, replace_indices_in_symbol, parse_algorithm_to_state_update, girf, decompose_name, obc_objective_optim_fun, obc_constraint_optim_fun, compute_irf_responses, process_ignore_obc_flag, adjust_generalised_irf_flag, process_shocks_input, normalize_filtering_options, infer_step -import MacroModelling: DEFAULT_ALGORITHM, DEFAULT_FILTER_SELECTOR, DEFAULT_WARMUP_ITERATIONS, DEFAULT_VARIABLES_EXCLUDING_OBC, DEFAULT_SHOCK_SELECTION, DEFAULT_PRESAMPLE_PERIODS, DEFAULT_DATA_IN_LEVELS, DEFAULT_SHOCK_DECOMPOSITION_SELECTOR, DEFAULT_SMOOTH_SELECTOR, DEFAULT_LABEL, DEFAULT_SHOW_PLOTS, DEFAULT_SAVE_PLOTS, DEFAULT_SAVE_PLOTS_FORMAT, DEFAULT_SAVE_PLOTS_PATH, DEFAULT_PLOTS_PER_PAGE_SMALL, DEFAULT_TRANSPARENCY, DEFAULT_MAX_ELEMENTS_PER_LEGEND_ROW, DEFAULT_EXTRA_LEGEND_SPACE, DEFAULT_VERBOSE, DEFAULT_QME_ALGORITHM, DEFAULT_SYLVESTER_SELECTOR, DEFAULT_SYLVESTER_THRESHOLD, DEFAULT_LARGE_SYLVESTER_ALGORITHM, DEFAULT_SYLVESTER_ALGORITHM, DEFAULT_LYAPUNOV_ALGORITHM, DEFAULT_PLOT_ATTRIBUTES, DEFAULT_ARGS_AND_KWARGS_NAMES, DEFAULT_PLOTS_PER_PAGE_LARGE, DEFAULT_SHOCKS_EXCLUDING_OBC, DEFAULT_VARIABLES_EXCLUDING_AUX_AND_OBC, DEFAULT_PERIODS, DEFAULT_SHOCK_SIZE, DEFAULT_NEGATIVE_SHOCK, DEFAULT_GENERALISED_IRF, DEFAULT_GENERALISED_IRF_WARMUP, DEFAULT_GENERALISED_IRF_DRAWS, DEFAULT_INITIAL_STATE, DEFAULT_IGNORE_OBC, DEFAULT_PLOT_TYPE, DEFAULT_CONDITIONS_IN_LEVELS, DEFAULT_SIGMA_RANGE, DEFAULT_FONT_SIZE, DEFAULT_VARIABLE_SELECTION, DEFAULT_FORECAST_PERIODS, DEFAULT_ALGORITHM_BACKWARD_LOOKING, DEFAULT_DEVIATIONS_FROM +import MacroModelling: DEFAULT_ALGORITHM, DEFAULT_FILTER_SELECTOR, DEFAULT_WARMUP_ITERATIONS, DEFAULT_VARIABLES_EXCLUDING_OBC, DEFAULT_SHOCK_SELECTION, DEFAULT_PRESAMPLE_PERIODS, DEFAULT_DATA_IN_LEVELS, DEFAULT_SHOCK_DECOMPOSITION_SELECTOR, DEFAULT_SMOOTH_SELECTOR, DEFAULT_LABEL, DEFAULT_SHOW_PLOTS, DEFAULT_SAVE_PLOTS, DEFAULT_SAVE_PLOTS_FORMAT, DEFAULT_SAVE_PLOTS_PATH, DEFAULT_PLOTS_PER_PAGE_SMALL, DEFAULT_TRANSPARENCY, DEFAULT_MAX_ELEMENTS_PER_LEGEND_ROW, DEFAULT_EXTRA_LEGEND_SPACE, DEFAULT_VERBOSE, DEFAULT_QME_ALGORITHM, DEFAULT_SYLVESTER_SELECTOR, DEFAULT_SYLVESTER_THRESHOLD, DEFAULT_LARGE_SYLVESTER_ALGORITHM, DEFAULT_SYLVESTER_ALGORITHM, DEFAULT_LYAPUNOV_ALGORITHM, DEFAULT_PLOT_ATTRIBUTES, DEFAULT_ARGS_AND_KWARGS_NAMES, DEFAULT_PLOTS_PER_PAGE_LARGE, DEFAULT_SHOCKS_EXCLUDING_OBC, DEFAULT_VARIABLES_EXCLUDING_AUX_AND_OBC, DEFAULT_PERIODS, DEFAULT_SHOCK_SIZE, DEFAULT_NEGATIVE_SHOCK, DEFAULT_GENERALISED_IRF, DEFAULT_GENERALISED_IRF_WARMUP, DEFAULT_GENERALISED_IRF_DRAWS, DEFAULT_INITIAL_STATE, DEFAULT_IGNORE_OBC, DEFAULT_PLOT_TYPE, DEFAULT_CONDITIONS_IN_LEVELS, DEFAULT_SIGMA_RANGE, DEFAULT_FONT_SIZE, DEFAULT_VARIABLE_SELECTION, DEFAULT_FORECAST_PERIODS, DEFAULT_ALGORITHM_BACKWARD_LOOKING, DEFAULT_REFERENCE import DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF, TYPEDSIGNATURES, TYPEDFIELDS import LaTeXStrings @@ -1751,7 +1751,7 @@ If the model contains occasionally binding constraints and `ignore_obc = false` - $INITIAL_STATE® - $IGNORE_OBC® - `reference` [Default: `:baseline` for backward looking models, `:steady_state` otherwise, Type: `Symbol`]: reference path for plotting. Options: `:steady_state` (horizontal line at steady state) or `:baseline` (no-shock trajectory from `initial_state`). Useful for backward looking models with explosive dynamics. -- `deviations_from` [Default: `:steady_state`, Type: `Symbol`]: reference point for deviations when not in levels mode. Options: `:steady_state` (deviations from relevant steady state) or `:baseline` (deviations from no-shock path starting from `initial_state`). +- `reference` [Default: `:steady_state`, Type: `Symbol`]: reference point for deviations when not in levels mode. Options: `:steady_state` (deviations from relevant steady state) or `:baseline` (deviations from no-shock path starting from `initial_state`). - `label` [Default: `1`, Type: `Union{Real, String, Symbol}`]: label to attribute to this function call in the plots. - $SHOW_PLOTS® - $SAVE_PLOTS® @@ -1813,7 +1813,7 @@ function plot_irf(𝓂::ℳ; initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = DEFAULT_INITIAL_STATE, ignore_obc::Bool = DEFAULT_IGNORE_OBC, reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, - deviations_from::Symbol = DEFAULT_DEVIATIONS_FROM, + reference::Symbol = DEFAULT_REFERENCE, rename_dictionary::AbstractDict{<:Union{Symbol, String}, <:Union{Symbol, String}} = Dict{Symbol, String}(), plot_attributes::Dict = Dict(), verbose::Bool = DEFAULT_VERBOSE, diff --git a/src/common_docstrings.jl b/src/common_docstrings.jl index 0d464b89a..8d3a712af 100644 --- a/src/common_docstrings.jl +++ b/src/common_docstrings.jl @@ -43,4 +43,4 @@ const INITIAL_STATE®1 = "`initial_state` [Default: `$(DEFAULT_INITIAL_STATE)`, const LABEL® = "`label` [Type: `Union{Real, String, Symbol}`]: label to attribute to this function call in the plots. The default is the number of previous function calls since the last call to the function version with ! + 1." const RENAME_DICTIONARY® = "`rename_dictionary` [Default: `Dict()`, Type: `Dict{Symbol, String}`]: dictionary mapping variable or shock symbols to custom display names in plots. For example: `Dict(:dinve => \"Investment growth\", :c => \"Consumption\")`. Variables/shocks not in the dictionary will use their default names." const CONDITIONS_IN_LEVELS® = "`conditions_in_levels` [Default: `true`, Type: `Bool`]: indicator whether the conditions are provided in levels. If `true` the input to the conditions argument will have the relevant steady state subtracted (non-stochastic or stochastic steady state depending on the solution algorithm)." -const DEVIATIONS_FROM® = "`deviations_from` [Default: `$(DEFAULT_DEVIATIONS_FROM)`, Type: `Symbol`]: reference point for computing deviations when `levels = false`. Options are `:steady_state` (deviations from the relevant steady state, the default behaviour) or `:baseline` (deviations from the no-shock path starting from `initial_state` - useful for backward looking models with explosive dynamics). Ignored when `levels = true`." \ No newline at end of file +const REFERENCE® = "`reference` [Default: `$(DEFAULT_REFERENCE)`, Type: `Symbol`]: reference point for computing deviations when `levels = false`. Options are `:steady_state` (deviations from the relevant steady state, the default behaviour) or `:baseline` (deviations from the no-shock path starting from `initial_state` - useful for backward looking models with explosive dynamics). Ignored when `levels = true`." \ No newline at end of file diff --git a/src/default_options.jl b/src/default_options.jl index acce1c809..3822fcd0f 100644 --- a/src/default_options.jl +++ b/src/default_options.jl @@ -7,7 +7,7 @@ const DEFAULT_ALGORITHM_SELECTOR = stochastic -> stochastic ? :second_order : :f const DEFAULT_ALGORITHM_BACKWARD_LOOKING = :newton # Deviations mode for IRF output -const DEFAULT_DEVIATIONS_FROM = :steady_state # :steady_state or :baseline +const DEFAULT_REFERENCE = :steady_state # :steady_state or :baseline const DEFAULT_FILTER_SELECTOR = algorithm -> algorithm == :first_order ? :kalman : :inversion const DEFAULT_SHOCK_DECOMPOSITION_SELECTOR = algorithm -> algorithm ∉ (:second_order, :third_order) const DEFAULT_SMOOTH_SELECTOR = filter -> filter == :kalman diff --git a/src/get_functions.jl b/src/get_functions.jl index f8fad3c7a..bfad6a6c6 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -1386,7 +1386,7 @@ If the model contains occasionally binding constraints and `ignore_obc = false` - $PARAMETERS® - $(VARIABLES®(DEFAULT_VARIABLES_EXCLUDING_OBC)) - $SHOCKS® -- $DEVIATIONS_FROM® +- $REFERENCE® - $NEGATIVE_SHOCK® - $GENERALISED_IRF® - $GENERALISED_IRF_WARMUP_ITERATIONS® @@ -1444,7 +1444,7 @@ function get_irf(𝓂::ℳ; parameters::ParameterType = nothing, variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, shocks::Union{Symbol_input,String_input,Matrix{Float64},KeyedArray{Float64}} = DEFAULT_SHOCKS_EXCLUDING_OBC, - deviations_from::Symbol = DEFAULT_DEVIATIONS_FROM, + reference::Symbol = DEFAULT_REFERENCE, negative_shock::Bool = DEFAULT_NEGATIVE_SHOCK, generalised_irf::Bool = DEFAULT_GENERALISED_IRF, generalised_irf_warmup_iterations::Int = DEFAULT_GENERALISED_IRF_WARMUP, @@ -1517,12 +1517,12 @@ function get_irf(𝓂::ℳ; @assert !unspecified_initial_state "Model is backward looking with no valid steady state. Provide initial_state in levels." end - # Validate deviations_from parameter - @assert deviations_from ∈ [:steady_state, :baseline] "deviations_from must be either :steady_state or :baseline" + # Validate reference parameter + @assert reference ∈ [:steady_state, :baseline] "reference must be either :steady_state or :baseline" - # If levels = true, ignore deviations_from (levels are levels) - if levels && deviations_from == :baseline - @warn "deviations_from = :baseline is ignored when levels = true (returning levels)" + # If levels = true, ignore reference (levels are levels) + if levels && reference == :baseline + @warn "reference = :baseline is ignored when levels = true (returning levels)" end # @timeit_debug timer "Get relevant steady state" begin @@ -1569,12 +1569,12 @@ function get_irf(𝓂::ℳ; level = levels ? reference_steady_state + SSS_delta : SSS_delta - # For backward looking models with deviations_from = :baseline, compute baseline path + # For backward looking models with reference = :baseline, compute baseline path # The baseline path is the no-shock forward iteration from initial_state - # When deviations_from = :baseline and levels = false: IRFs will show deviations from this baseline (shock effect) + # When reference = :baseline and levels = false: IRFs will show deviations from this baseline (shock effect) # When levels = true: baseline_path is not used (return simulation in levels) baseline_path = nothing - if !levels && deviations_from == :baseline && is_backward_looking && algorithm == :newton + if !levels && reference == :baseline && is_backward_looking && algorithm == :newton # Compute the no-shock baseline path in deviations from NSSS nVars = 𝓂.timings.nVars baseline_path = zeros(nVars, periods) diff --git a/test/test_backward_looking_models.jl b/test/test_backward_looking_models.jl index 50b5af13d..f438b9f2e 100644 --- a/test/test_backward_looking_models.jl +++ b/test/test_backward_looking_models.jl @@ -100,7 +100,7 @@ @test size(irf_first, 3) == 1 # 1 shock # Test simulation with newton - sim_newton = simulate(SolowGrowth, algorithm = :newton) + sim_newton = simulate(SolowGrowth) @test size(sim_newton, 1) == 2 # 2 variables # Test conditional forecasting with newton @@ -111,7 +111,7 @@ conditions[z_idx, 1] = 0.02 # Condition z in period 1 conditions[z_idx, 2] = 0.01 # Condition z in period 2 - cf_newton = get_conditional_forecast(SolowGrowth, conditions, algorithm = :newton) + cf_newton = get_conditional_forecast(SolowGrowth, conditions) # Check that conditions are met @test isapprox(cf_newton(:z, 1), 0.02, rtol=1e-4) @@ -146,7 +146,7 @@ # For backward looking models with newton algorithm, initial_state is in levels # shocks = :none returns result with 3rd dimension as [:none] result = get_irf(SolowGrowth2, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_solow, shocks = :none, levels = true, # Request levels output @@ -205,7 +205,7 @@ initial_state_levels = [initial_y, 0.0] # [y, z] result = get_irf(UnstableButValidSS, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_levels, shocks = :none, levels = true, # Request levels output @@ -220,24 +220,24 @@ @test y_values[1] ≈ initial_y * (1 + 0.02) rtol=1e-6 # y_1 = (1+g) * y_0 @test y_values[2] > y_values[1] # Continues growing - # Test deviations_from = :baseline + # Test reference = :baseline # Get IRF with shocks in deviations from baseline mode irf_baseline = get_irf(UnstableButValidSS, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_levels, - deviations_from = :baseline, + reference = :baseline, periods = 10) # Get IRF in levels mode irf_lev = get_irf(UnstableButValidSS, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_levels, levels = true, periods = 10) # Get baseline path (no shocks) in levels baseline_lev = get_irf(UnstableButValidSS, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_levels, shocks = :none, levels = true, @@ -262,7 +262,7 @@ cf_result = get_conditional_forecast(UnstableButValidSS, conditions, - algorithm = :newton, + # algorithm = :newton, initial_state = initial_state_levels, periods = 5) From 594803863a08d5dcfdc7702567c961a06de50d19 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 21:24:58 +0000 Subject: [PATCH 35/40] Default algorithm for backward looking models is :newton in all functions Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 7 +++---- src/get_functions.jl | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 609e8cae8..51b82121e 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -1813,7 +1813,6 @@ function plot_irf(𝓂::ℳ; initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = DEFAULT_INITIAL_STATE, ignore_obc::Bool = DEFAULT_IGNORE_OBC, reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, - reference::Symbol = DEFAULT_REFERENCE, rename_dictionary::AbstractDict{<:Union{Symbol, String}, <:Union{Symbol, String}} = Dict{Symbol, String}(), plot_attributes::Dict = Dict(), verbose::Bool = DEFAULT_VERBOSE, @@ -2488,7 +2487,7 @@ function plot_irf!(𝓂::ℳ; save_plots_name::Union{String, Symbol} = "irf", save_plots_path::String = DEFAULT_SAVE_PLOTS_PATH, plots_per_page::Int = DEFAULT_PLOTS_PER_PAGE_SMALL, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, shock_size::Real = DEFAULT_SHOCK_SIZE, negative_shock::Bool = DEFAULT_NEGATIVE_SHOCK, generalised_irf::Bool = DEFAULT_GENERALISED_IRF, @@ -4800,7 +4799,7 @@ function plot_conditional_forecast(𝓂::ℳ, parameters::ParameterType = nothing, variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, conditions_in_levels::Bool = DEFAULT_CONDITIONS_IN_LEVELS, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, label::Union{Real, String, Symbol} = DEFAULT_LABEL, show_plots::Bool = DEFAULT_SHOW_PLOTS, save_plots::Bool = DEFAULT_SAVE_PLOTS, @@ -5256,7 +5255,7 @@ function plot_conditional_forecast!(𝓂::ℳ, parameters::ParameterType = nothing, variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, conditions_in_levels::Bool = DEFAULT_CONDITIONS_IN_LEVELS, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, label::Union{Real, String, Symbol} = length(conditional_forecast_active_plot_container) + 1, show_plots::Bool = DEFAULT_SHOW_PLOTS, save_plots::Bool = DEFAULT_SAVE_PLOTS, diff --git a/src/get_functions.jl b/src/get_functions.jl index bfad6a6c6..160ddb9f1 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -737,7 +737,7 @@ function get_conditional_forecast(𝓂::ℳ, parameters::ParameterType = nothing, variables::Union{Symbol_input,String_input} = DEFAULT_VARIABLES_EXCLUDING_OBC, conditions_in_levels::Bool = DEFAULT_CONDITIONS_IN_LEVELS, - algorithm::Symbol = DEFAULT_ALGORITHM, + algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, levels::Bool = false, verbose::Bool = DEFAULT_VERBOSE, tol::Tolerances = Tolerances(), From f23d759a9b81a8cf57764feceeff5c9cf226dd45 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 21:54:14 +0000 Subject: [PATCH 36/40] Add reference argument to plot_irf, plot_irf!, plot_conditional_forecast, plot_conditional_forecast! Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 78 +++++++++++++++++++++++++++++++------------- 1 file changed, 56 insertions(+), 22 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 51b82121e..d8235bc7e 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -17,6 +17,7 @@ import Showoff import DataStructures: OrderedSet import SparseArrays: SparseMatrixCSC import NLopt +import Statistics: mean using DispatchDoctor import MacroModelling: plot_irfs, plot_irf, plot_IRF, plot_simulations, plot_simulation, plot_solution, plot_girf, plot_conditional_forecast, plot_conditional_variance_decomposition, plot_forecast_error_variance_decomposition, plot_fevd, plot_model_estimates, plot_shock_decomposition, plotlyjs_backend, gr_backend, compare_args_and_kwargs, get_irf @@ -1893,6 +1894,25 @@ function plot_irf(𝓂::ℳ; state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂, false) end + # Compute baseline path when reference = :baseline for backward looking models + is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 + baseline_path_for_plot = nothing + + if reference == :baseline && is_backward_looking && algorithm == :newton + nVars = 𝓂.timings.nVars + baseline_path_for_plot = zeros(nVars, periods_extended) + zero_shocks = zeros(length(𝓂.timings.exo)) + + # Compute baseline in deviations from NSSS + baseline_state = copy(initial_state) + for t in 1:periods_extended + baseline_state = state_update(baseline_state, zero_shocks) + baseline_path_for_plot[:, t] = baseline_state + end + # Convert to levels for plotting: baseline_path_for_plot + NSSS + baseline_path_for_plot = baseline_path_for_plot .+ NSSS[1:nVars] + end + level = zeros(𝓂.timings.nVars) Y = compute_irf_responses(𝓂, @@ -2027,11 +2047,14 @@ function plot_irf(𝓂::ℳ; for (i,v) in enumerate(var_idx) SS = reference_steady_state[v] + + # Get baseline path for this variable if using baseline reference + var_baseline_path = isnothing(baseline_path_for_plot) ? nothing : baseline_path_for_plot[v, :] if !(all(isapprox.(Y[i,:,shock],0,atol = eps(Float32)))) variable_name = variable_names_display[i] - push!(pp, standard_subplot(Y[i,:,shock], SS, variable_name, gr_back, pal = pal)) + push!(pp, standard_subplot(Y[i,:,shock], SS, variable_name, gr_back, pal = pal, baseline_path = var_baseline_path)) if !(plot_count % plots_per_page == 0) plot_count += 1 @@ -2113,43 +2136,51 @@ function standard_subplot(irf_data::AbstractVector{S}, variable_name::R, gr_back::Bool; pal::StatsPlots.ColorPalette = StatsPlots.palette(:auto), - xvals = 1:length(irf_data)) where {S <: AbstractFloat, R <: Union{String, Symbol}} - finite_vals = filter(isfinite, irf_data) - can_dual_axis = gr_back && !isempty(finite_vals) && all((finite_vals .+ steady_state) .> eps(Float32)) && (steady_state > eps(Float32)) + xvals = 1:length(irf_data), + baseline_path::Union{Nothing, AbstractVector{S}} = nothing) where {S <: AbstractFloat, R <: Union{String, Symbol}} + # If baseline_path is provided, use it; otherwise use steady_state for reference line + use_baseline = !isnothing(baseline_path) + + if use_baseline + # baseline_path is in levels - use it for reference + reference_line = baseline_path + # irf_data is in deviations from baseline, so we add baseline_path to get levels + plot_data = irf_data .+ baseline_path + else + reference_line = fill(steady_state, length(irf_data)) + plot_data = irf_data .+ steady_state + end + + finite_vals = filter(isfinite, plot_data) + ref_for_dual_axis = use_baseline ? mean(filter(isfinite, baseline_path)) : steady_state + can_dual_axis = gr_back && !isempty(finite_vals) && all(finite_vals .> eps(Float32)) && (ref_for_dual_axis > eps(Float32)) xrotation = length(string(xvals[1])) > 5 ? 30 : 0 p = StatsPlots.plot(xvals, - irf_data .+ steady_state, + plot_data, title = variable_name, ylabel = "Level", xrotation = xrotation, color = pal[1], label = "") - - StatsPlots.hline!([steady_state], + + # Plot reference line (either baseline path or horizontal steady state) + if use_baseline + StatsPlots.plot!(p, xvals, reference_line, color = :black, label = "") + else + StatsPlots.hline!([steady_state], + color = :black, + label = "") + end lo, hi = StatsPlots.ylims(p) - # if !(xvals isa UnitRange) - # low = 1 - # high = length(irf_data) - - # # Compute nice ticks on the shifted range - # ticks_shifted, _ = StatsPlots.optimize_ticks(low, high, k_min = 4, k_max = 6) - - # ticks_shifted = Int.(ceil.(ticks_shifted)) - - # labels = xvals[ticks_shifted] - - # StatsPlots.plot!(xticks = (ticks_shifted, labels)) - # end - if can_dual_axis StatsPlots.plot!(StatsPlots.twinx(), - ylims = (100 * (lo / steady_state - 1), 100 * (hi / steady_state - 1)), + ylims = (100 * (lo / ref_for_dual_axis - 1), 100 * (hi / ref_for_dual_axis - 1)), xrotation = xrotation, ylabel = LaTeXStrings.L"\% \Delta") end @@ -2496,6 +2527,7 @@ function plot_irf!(𝓂::ℳ; initial_state::Union{Vector{Vector{Float64}},Vector{Float64}} = DEFAULT_INITIAL_STATE, ignore_obc::Bool = DEFAULT_IGNORE_OBC, plot_type::Symbol = DEFAULT_PLOT_TYPE, + reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, rename_dictionary::AbstractDict{<:Union{Symbol, String}, <:Union{Symbol, String}} = Dict{Symbol, String}(), plot_attributes::Dict = Dict(), transparency::Float64 = DEFAULT_TRANSPARENCY, @@ -4801,6 +4833,7 @@ function plot_conditional_forecast(𝓂::ℳ, conditions_in_levels::Bool = DEFAULT_CONDITIONS_IN_LEVELS, algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, label::Union{Real, String, Symbol} = DEFAULT_LABEL, + reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, show_plots::Bool = DEFAULT_SHOW_PLOTS, save_plots::Bool = DEFAULT_SAVE_PLOTS, save_plots_format::Symbol = DEFAULT_SAVE_PLOTS_FORMAT, @@ -5257,6 +5290,7 @@ function plot_conditional_forecast!(𝓂::ℳ, conditions_in_levels::Bool = DEFAULT_CONDITIONS_IN_LEVELS, algorithm::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? DEFAULT_ALGORITHM_BACKWARD_LOOKING : DEFAULT_ALGORITHM, label::Union{Real, String, Symbol} = length(conditional_forecast_active_plot_container) + 1, + reference::Symbol = 𝓂.timings.nFuture_not_past_and_mixed == 0 ? :baseline : :steady_state, show_plots::Bool = DEFAULT_SHOW_PLOTS, save_plots::Bool = DEFAULT_SAVE_PLOTS, save_plots_format::Symbol = DEFAULT_SAVE_PLOTS_FORMAT, From ef331872c092e745839b85f8066583322915bf48 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 22:31:58 +0000 Subject: [PATCH 37/40] Remove Statistics dependency by replacing mean with sum/length Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index d8235bc7e..cb88db525 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -17,7 +17,7 @@ import Showoff import DataStructures: OrderedSet import SparseArrays: SparseMatrixCSC import NLopt -import Statistics: mean + using DispatchDoctor import MacroModelling: plot_irfs, plot_irf, plot_IRF, plot_simulations, plot_simulation, plot_solution, plot_girf, plot_conditional_forecast, plot_conditional_variance_decomposition, plot_forecast_error_variance_decomposition, plot_fevd, plot_model_estimates, plot_shock_decomposition, plotlyjs_backend, gr_backend, compare_args_and_kwargs, get_irf @@ -2152,7 +2152,8 @@ function standard_subplot(irf_data::AbstractVector{S}, end finite_vals = filter(isfinite, plot_data) - ref_for_dual_axis = use_baseline ? mean(filter(isfinite, baseline_path)) : steady_state + finite_baseline_vals = filter(isfinite, baseline_path) + ref_for_dual_axis = use_baseline ? (sum(finite_baseline_vals) / length(finite_baseline_vals)) : steady_state can_dual_axis = gr_back && !isempty(finite_vals) && all(finite_vals .> eps(Float32)) && (ref_for_dual_axis > eps(Float32)) xrotation = length(string(xvals[1])) > 5 ? 30 : 0 From 4ccfa83045cd70060a049701d46682ebffda6543 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 23:41:53 +0100 Subject: [PATCH 38/40] Refactor reference line calculation in standard_subplot for improved clarity --- ext/StatsPlotsExt.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index cb88db525..e5a0acaea 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -2146,14 +2146,15 @@ function standard_subplot(irf_data::AbstractVector{S}, reference_line = baseline_path # irf_data is in deviations from baseline, so we add baseline_path to get levels plot_data = irf_data .+ baseline_path + finite_baseline_vals = filter(isfinite, baseline_path) + ref_for_dual_axis = (sum(finite_baseline_vals) / length(finite_baseline_vals)) else reference_line = fill(steady_state, length(irf_data)) plot_data = irf_data .+ steady_state + ref_for_dual_axis = steady_state end finite_vals = filter(isfinite, plot_data) - finite_baseline_vals = filter(isfinite, baseline_path) - ref_for_dual_axis = use_baseline ? (sum(finite_baseline_vals) / length(finite_baseline_vals)) : steady_state can_dual_axis = gr_back && !isempty(finite_vals) && all(finite_vals .> eps(Float32)) && (ref_for_dual_axis > eps(Float32)) xrotation = length(string(xvals[1])) > 5 ? 30 : 0 From b27951703a1caa8863ac489a1c4aa8915bd13ba8 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Fri, 2 Jan 2026 23:46:16 +0100 Subject: [PATCH 39/40] Fix standard_subplot to use finite baseline values for reference and adjust plot_data calculation --- ext/StatsPlotsExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index e5a0acaea..2787f801a 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -2145,9 +2145,9 @@ function standard_subplot(irf_data::AbstractVector{S}, # baseline_path is in levels - use it for reference reference_line = baseline_path # irf_data is in deviations from baseline, so we add baseline_path to get levels - plot_data = irf_data .+ baseline_path + plot_data = irf_data # .+ baseline_path finite_baseline_vals = filter(isfinite, baseline_path) - ref_for_dual_axis = (sum(finite_baseline_vals) / length(finite_baseline_vals)) + ref_for_dual_axis = first(finite_baseline_vals) else reference_line = fill(steady_state, length(irf_data)) plot_data = irf_data .+ steady_state From 7b76497dfc94c6ee45f21e3e7c1769550d49713c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 2 Jan 2026 22:52:48 +0000 Subject: [PATCH 40/40] Use initial state as reference for dual axis instead of first period value Co-authored-by: thorek1 <13523097+thorek1@users.noreply.github.com> --- ext/StatsPlotsExt.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 2787f801a..4d1f7c1b9 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -1897,12 +1897,16 @@ function plot_irf(𝓂::ℳ; # Compute baseline path when reference = :baseline for backward looking models is_backward_looking = 𝓂.timings.nFuture_not_past_and_mixed == 0 baseline_path_for_plot = nothing + initial_state_levels = nothing # Store initial state in levels for dual axis reference if reference == :baseline && is_backward_looking && algorithm == :newton nVars = 𝓂.timings.nVars baseline_path_for_plot = zeros(nVars, periods_extended) zero_shocks = zeros(length(𝓂.timings.exo)) + # Store initial state in levels for dual axis reference + initial_state_levels = initial_state .+ NSSS[1:nVars] + # Compute baseline in deviations from NSSS baseline_state = copy(initial_state) for t in 1:periods_extended @@ -2050,11 +2054,13 @@ function plot_irf(𝓂::ℳ; # Get baseline path for this variable if using baseline reference var_baseline_path = isnothing(baseline_path_for_plot) ? nothing : baseline_path_for_plot[v, :] + # Get initial state value for this variable (in levels) for dual axis reference + var_initial_value = isnothing(initial_state_levels) ? nothing : initial_state_levels[v] if !(all(isapprox.(Y[i,:,shock],0,atol = eps(Float32)))) variable_name = variable_names_display[i] - push!(pp, standard_subplot(Y[i,:,shock], SS, variable_name, gr_back, pal = pal, baseline_path = var_baseline_path)) + push!(pp, standard_subplot(Y[i,:,shock], SS, variable_name, gr_back, pal = pal, baseline_path = var_baseline_path, initial_value = var_initial_value)) if !(plot_count % plots_per_page == 0) plot_count += 1 @@ -2137,7 +2143,8 @@ function standard_subplot(irf_data::AbstractVector{S}, gr_back::Bool; pal::StatsPlots.ColorPalette = StatsPlots.palette(:auto), xvals = 1:length(irf_data), - baseline_path::Union{Nothing, AbstractVector{S}} = nothing) where {S <: AbstractFloat, R <: Union{String, Symbol}} + baseline_path::Union{Nothing, AbstractVector{S}} = nothing, + initial_value::Union{Nothing, S} = nothing) where {S <: AbstractFloat, R <: Union{String, Symbol}} # If baseline_path is provided, use it; otherwise use steady_state for reference line use_baseline = !isnothing(baseline_path) @@ -2146,12 +2153,12 @@ function standard_subplot(irf_data::AbstractVector{S}, reference_line = baseline_path # irf_data is in deviations from baseline, so we add baseline_path to get levels plot_data = irf_data # .+ baseline_path - finite_baseline_vals = filter(isfinite, baseline_path) - ref_for_dual_axis = first(finite_baseline_vals) + # Use initial_value as reference for dual axis if provided, otherwise use steady_state + ref_for_dual_axis = isnothing(initial_value) ? steady_state : initial_value else reference_line = fill(steady_state, length(irf_data)) plot_data = irf_data .+ steady_state - ref_for_dual_axis = steady_state + ref_for_dual_axis = steady_state end finite_vals = filter(isfinite, plot_data)