From cc624d91f39c54763c26771c43bfd9d712e7e337 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Mon, 8 Jun 2026 21:11:44 +0000 Subject: [PATCH 1/7] filter free plotting attempt --- ext/StatsPlotsExt.jl | 294 ++++++++++++++++++++++++++++++++++++++----- src/get_functions.jl | 241 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 + 3 files changed, 508 insertions(+), 29 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 702cc0567..add77bcc5 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -2,7 +2,7 @@ module StatsPlotsExt using MacroModelling -import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, NsssTolerances, SolverTolerances, 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®, STEADY_STATE_FUNCTION®, 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, normalize_presample_periods, trim_informative_sample, adjust_initial_state, SteadyStateFunctionType +import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, NsssTolerances, SolverTolerances, 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®, STEADY_STATE_FUNCTION®, parse_shocks_input_to_index, parse_variables_input_to_index, replace_indices, replace_indices_special, filter_data_with_model, filter_free_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, normalize_presample_periods, trim_informative_sample, adjust_initial_state, SteadyStateFunctionType import MacroModelling: DEFAULT_CACHING, DEFAULT_USE_WORKSPACES, 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 DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF, TYPEDSIGNATURES, TYPEDFIELDS import LaTeXStrings @@ -364,6 +364,30 @@ function build_extended_palette(attributes_redux::Dict; total_pal_len::Int = 100 mapreduce(x -> StatsPlots.coloralpha.(orig_pal, alpha_reduction_factor ^ x), vcat, 0:(total_pal_len ÷ length(orig_pal)) - 1) |> StatsPlots.palette end +function plot_observation_std_band!(p, + xvals, + observations::AbstractVector, + std_values::AbstractVector, + color; + alpha::Float64 = 0.18) + StatsPlots.plot!(p, + xvals, + observations, + ribbon = std_values, + fillalpha = alpha, + linealpha = 0, + linewidth = 0, + label = "", + color = color) + return p +end + +function has_observed_series(variable_idx::Int, obs_idx::AbstractVector{<:Integer}, data_in_deviations, periods) + obs_pos = findfirst(==(variable_idx), obs_idx) + isnothing(obs_pos) && return false + return any(isfinite, vec(data_in_deviations[obs_pos, periods])) +end + function process_rename_dictionary(rename_dictionary::AbstractDict, 𝓂::ℳ) relevant_keys = [k for k in keys(rename_dictionary) if (k isa String ? replace_indices(k) : k) in vcat(𝓂.constants.post_model_macro.var, 𝓂.constants.post_model_macro.exo)] |> sort processed = Any[] @@ -634,6 +658,8 @@ If occasionally binding constraints are present in the model, they are not taken # Arguments - $MODEL® - $DATA® +- `filter_free_shocks` [optional positional, Type: `AbstractMatrix{<:Real}`]: latent structural shock path from filter-free estimation. The matrix must have one row per exogenous shock and one column per data period, plus `max(warmup_iterations - 1, 0)` leading warmup columns when warmup is used. When supplied, the model-implied endogenous path is forward-simulated from these shocks. +- `measurement_error_std` [optional positional, Type: `Real`, `AbstractVector`, or `AbstractMatrix`]: standard deviation(s) for a transparent [+σ,-σ] band around observations. Vectors must have one entry per observable; matrices must have dimensions `(n_observables, n_periods)`. # Keyword Arguments - $PARAMETERS® - $STEADY_STATE_FUNCTION® @@ -700,6 +726,21 @@ simulation = simulate(RBC_CME) plot_model_estimates(RBC_CME, simulation([:k],:,:simulate)) ``` """ +function plot_model_estimates(𝓂::ℳ, + data::KeyedArray, + filter_free_shocks::AbstractMatrix{<:Real}; + kwargs...) + return plot_model_estimates(𝓂, data; filter_free_shocks = filter_free_shocks, kwargs...) +end + +function plot_model_estimates(𝓂::ℳ, + data::KeyedArray, + filter_free_shocks::AbstractMatrix{<:Real}, + measurement_error_std::Union{Real,AbstractVector{<:Real},AbstractMatrix{<:Real}}; + kwargs...) + return plot_model_estimates(𝓂, data; filter_free_shocks = filter_free_shocks, measurement_error_std = measurement_error_std, kwargs...) +end + function plot_model_estimates(𝓂::ℳ, data::KeyedArray; parameters::ParameterType = nothing, @@ -733,7 +774,10 @@ function plot_model_estimates(𝓂::ℳ, sylvester_algorithm::Union{Symbol,Vector{Symbol},Tuple{Symbol,Vararg{Symbol}}} = DEFAULT_SYLVESTER_SELECTOR(𝓂), lyapunov_algorithm::Symbol = DEFAULT_LYAPUNOV_ALGORITHM, caching::Bool = DEFAULT_CACHING, - use_workspaces::Bool = DEFAULT_USE_WORKSPACES) + use_workspaces::Bool = DEFAULT_USE_WORKSPACES, + filter_free_shocks::Union{Nothing,AbstractMatrix{<:Real}} = nothing, + measurement_error_std::Union{Nothing,Real,AbstractVector{<:Real},AbstractMatrix{<:Real}} = nothing, + initial_state = DEFAULT_INITIAL_STATE) # @nospecialize # reduce compile time if !caching invalidate_cache_validity!(𝓂) end @@ -745,14 +789,26 @@ function plot_model_estimates(𝓂::ℳ, sylvester_algorithm² = isa(sylvester_algorithm, Symbol) ? sylvester_algorithm : sylvester_algorithm[1], sylvester_algorithm³ = (isa(sylvester_algorithm, Symbol) || length(sylvester_algorithm) < 2) ? sum(k * (k + 1) ÷ 2 for k in 1:𝓂.constants.post_model_macro.nPast_not_future_and_mixed + 1 + 𝓂.constants.post_model_macro.nExo) > DEFAULT_SYLVESTER_THRESHOLD ? DEFAULT_LARGE_SYLVESTER_ALGORITHM : DEFAULT_SYLVESTER_ALGORITHM : sylvester_algorithm[2], lyapunov_algorithm = lyapunov_algorithm) - warn_irrelevant_tol(tol, algorithm; needs_covariance = filter == :kalman) + warn_irrelevant_tol(tol, algorithm; needs_covariance = isnothing(filter_free_shocks) && filter == :kalman) gr_back, attributes, attributes_redux = setup_plot_attributes(plot_attributes) # write_parameters_input!(𝓂, parameters, verbose = verbose) - filter, smooth, algorithm, shock_decomposition, pruning, warmup_iterations = normalize_filtering_options(filter, smooth, algorithm, shock_decomposition, warmup_iterations) + filter_free_plot = !isnothing(filter_free_shocks) + if filter_free_plot + @assert algorithm ∈ [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] "`plot_model_estimates` with positional shocks only supports perturbation algorithms (`:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`)." + if shock_decomposition + @info "Shock decomposition is not available when plotting a positional filter-free shock path. Setting `shock_decomposition = false`." maxlog = 3 + shock_decomposition = false + end + smooth = false + pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order) + filter = :filter_free + else + filter, smooth, algorithm, shock_decomposition, pruning, warmup_iterations = normalize_filtering_options(filter, smooth, algorithm, shock_decomposition, warmup_iterations) + end if marginal_contribution && shock_decomposition && !pruning @info "`marginal_contribution = true` is only meaningful for pruned higher-order solutions (`:pruned_second_order`, `:pruned_third_order`). Setting `marginal_contribution = false` for `algorithm = $(algorithm)`." maxlog = 3 @@ -809,7 +865,7 @@ function plot_model_estimates(𝓂::ℳ, legend_columns = 1 - legend_items = length(shock_idx) + 3 + pruning + (forecast_periods > 0 ? 1 : 0) + legend_items = length(shock_idx) + 3 + pruning + (forecast_periods > 0 ? 1 : 0) + (filter_free_plot && !isnothing(measurement_error_std) ? 1 : 0) max_columns = min(legend_items, max_elements_per_legend_row) @@ -827,8 +883,22 @@ function plot_model_estimates(𝓂::ℳ, data_in_deviations = MacroModelling.missing_data_to_nan(data) end - data_in_deviations, _, _, informative_periods = trim_informative_sample(data_in_deviations; - require_informative_periods = true) + measurement_error_std_to_plot = nothing + if filter_free_plot + filter_free_output = filter_free_data_with_model(𝓂, + data_in_deviations, + filter_free_shocks, + measurement_error_std, + initial_state; + algorithm = algorithm, + warmup_iterations = warmup_iterations, + opts = opts) + data_in_deviations = filter_free_output.data_in_deviations + informative_periods = filter_free_output.period_range + else + data_in_deviations, _, _, informative_periods = trim_informative_sample(data_in_deviations; + require_informative_periods = true) + end presample_periods = normalize_presample_periods(presample_periods, size(data_in_deviations, 2)) x_axis = axiskeys(data,2)[informative_periods] @@ -839,10 +909,18 @@ function plot_model_estimates(𝓂::ℳ, x_axis = x_axis[periods] - extra_kw = mc ? (; marginal_contribution = true) : NamedTuple() - variables_to_plot, shocks_to_plot, standard_deviations, decomposition = filter_data_with_model(𝓂, data_in_deviations, Val(algorithm), Val(filter), warmup_iterations = warmup_iterations, smooth = smooth, opts = opts; extra_kw...) + if filter_free_plot + variables_to_plot = filter_free_output.variables + shocks_to_plot = filter_free_output.shocks + standard_deviations = zeros(0, 0) + decomposition = zeros(size(variables_to_plot, 1), size(shocks_to_plot, 1) + 2, size(variables_to_plot, 2)) + measurement_error_std_to_plot = filter_free_output.measurement_error_std + else + extra_kw = mc ? (; marginal_contribution = true) : NamedTuple() + variables_to_plot, shocks_to_plot, standard_deviations, decomposition = filter_data_with_model(𝓂, data_in_deviations, Val(algorithm), Val(filter), warmup_iterations = warmup_iterations, smooth = smooth, opts = opts; extra_kw...) + end - if is_pruned + if is_pruned && !filter_free_plot if mc decomposition[:, end - 1, :] .+= SSS_delta else @@ -924,6 +1002,8 @@ function plot_model_estimates(𝓂::ℳ, :decomposition => decomposition, :variables_to_plot => variables_to_plot[var_idx, :], :data_in_deviations => data_in_deviations, + :measurement_error_std => measurement_error_std_to_plot, + :filter_free_plot => filter_free_plot, :shocks_to_plot => shocks_to_plot, :reference_steady_state => reference_steady_state[var_idx], :variable_names => variable_names_display, @@ -945,7 +1025,8 @@ function plot_model_estimates(𝓂::ℳ, plot_count = 1 for v in var_idx - if all(isapprox.(variables_to_plot[v, periods], 0, atol = eps(Float32))) + if all(isapprox.(variables_to_plot[v, periods], 0, atol = eps(Float32))) && + !has_observed_series(v, obs_idx, data_in_deviations, periods) n_subplots -= 1 end end @@ -977,7 +1058,8 @@ function plot_model_estimates(𝓂::ℳ, continue end else - if !(all(isapprox.(variables_to_plot[var_idx[i],periods], 0, atol = eps(Float32)))) + if !(all(isapprox.(variables_to_plot[var_idx[i],periods], 0, atol = eps(Float32))) && + !has_observed_series(var_idx[i], obs_idx, data_in_deviations, periods)) SS = reference_steady_state[var_idx[i]] if shock_decomposition @@ -1020,10 +1102,23 @@ function plot_model_estimates(𝓂::ℳ, color_total = estimate_color) if var_idx[i] ∈ obs_idx # Pad data with NaN for forecast period + obs_row = indexin([var_idx[i]],obs_idx) data_padded = if forecast_periods > 0 - vcat(vec(data_in_deviations[indexin([var_idx[i]],obs_idx),periods]), fill(NaN, forecast_periods)) + vcat(vec(data_in_deviations[obs_row,periods]), fill(NaN, forecast_periods)) else - vec(data_in_deviations[indexin([var_idx[i]],obs_idx),periods]) + vec(data_in_deviations[obs_row,periods]) + end + if !isnothing(measurement_error_std_to_plot) + std_padded = if forecast_periods > 0 + vcat(vec(measurement_error_std_to_plot[obs_row, periods]), fill(NaN, forecast_periods)) + else + vec(measurement_error_std_to_plot[obs_row, periods]) + end + plot_observation_std_band!(p, + extended_x_axis, + shock_decomposition ? data_padded : data_padded .+ SS, + std_padded, + shock_decomposition ? data_color : pal[2]) end StatsPlots.plot!(p, # extended_x_axis, @@ -1064,10 +1159,23 @@ function plot_model_estimates(𝓂::ℳ, if var_idx[i] ∈ obs_idx # Pad data with NaN for forecast period + obs_row = indexin([var_idx[i]],obs_idx) data_padded = if forecast_periods > 0 - vcat(vec(data_in_deviations[indexin([var_idx[i]],obs_idx),periods]), fill(NaN, forecast_periods)) + vcat(vec(data_in_deviations[obs_row,periods]), fill(NaN, forecast_periods)) else - vec(data_in_deviations[indexin([var_idx[i]],obs_idx),periods]) + vec(data_in_deviations[obs_row,periods]) + end + if !isnothing(measurement_error_std_to_plot) + std_padded = if forecast_periods > 0 + vcat(vec(measurement_error_std_to_plot[obs_row, periods]), fill(NaN, forecast_periods)) + else + vec(measurement_error_std_to_plot[obs_row, periods]) + end + plot_observation_std_band!(p, + extended_x_axis, + shock_decomposition ? data_padded : data_padded .+ SS, + std_padded, + shock_decomposition ? data_color : pal[2]) end StatsPlots.plot!(p, extended_x_axis, @@ -1129,6 +1237,18 @@ function plot_model_estimates(𝓂::ℳ, label = "Data", color = shock_decomposition ? data_color : pal[2]) + if !isnothing(measurement_error_std_to_plot) + StatsPlots.plot!(pl, + [NaN], [NaN], + seriestype = :scatter, + label = "[+σ,-σ]", + markershape = :rect, + markersize = 8, + markeralpha = 0.18, + markercolor = shock_decomposition ? data_color : pal[2], + markerstrokewidth = 0) + end + if shock_decomposition additional_labels_prefix = ["Initial value"] additional_labels_suffix = pruning && !mc ? ["Nonlinearities"] : String[] @@ -1200,6 +1320,18 @@ function plot_model_estimates(𝓂::ℳ, label = "Data", color = shock_decomposition ? data_color : pal[2]) + if !isnothing(measurement_error_std_to_plot) + StatsPlots.plot!(pl, + [NaN], [NaN], + seriestype = :scatter, + label = "[+σ,-σ]", + markershape = :rect, + markersize = 8, + markeralpha = 0.18, + markercolor = shock_decomposition ? data_color : pal[2], + markerstrokewidth = 0) + end + if shock_decomposition additional_labels_prefix = ["Initial value"] @@ -1271,6 +1403,8 @@ This function shares most of the signature and functionality of [`plot_model_est # Arguments - $MODEL® - $DATA® +- `filter_free_shocks` [optional positional, Type: `AbstractMatrix{<:Real}`]: latent structural shock path from filter-free estimation. The matrix must have one row per exogenous shock and one column per data period, plus `max(warmup_iterations - 1, 0)` leading warmup columns when warmup is used. When supplied, the model-implied endogenous path is forward-simulated from these shocks. +- `measurement_error_std` [optional positional, Type: `Real`, `AbstractVector`, or `AbstractMatrix`]: standard deviation(s) for a transparent [+σ,-σ] band around observations. Vectors must have one entry per observable; matrices must have dimensions `(n_observables, n_periods)`. # Keyword Arguments - $PARAMETERS® - $STEADY_STATE_FUNCTION® @@ -1354,6 +1488,21 @@ plot_model_estimates(RBC_CME, simulation([:k],:,:simulate)) plot_model_estimates!(RBC_CME, simulation([:k],:,:simulate), parameters = :beta => .99) ``` """ +function plot_model_estimates!(𝓂::ℳ, + data::KeyedArray, + filter_free_shocks::AbstractMatrix{<:Real}; + kwargs...) + return plot_model_estimates!(𝓂, data; filter_free_shocks = filter_free_shocks, kwargs...) +end + +function plot_model_estimates!(𝓂::ℳ, + data::KeyedArray, + filter_free_shocks::AbstractMatrix{<:Real}, + measurement_error_std::Union{Real,AbstractVector{<:Real},AbstractMatrix{<:Real}}; + kwargs...) + return plot_model_estimates!(𝓂, data; filter_free_shocks = filter_free_shocks, measurement_error_std = measurement_error_std, kwargs...) +end + function plot_model_estimates!(𝓂::ℳ, data::KeyedArray; parameters::ParameterType = nothing, @@ -1384,7 +1533,10 @@ function plot_model_estimates!(𝓂::ℳ, sylvester_algorithm::Union{Symbol,Vector{Symbol},Tuple{Symbol,Vararg{Symbol}}} = DEFAULT_SYLVESTER_SELECTOR(𝓂), lyapunov_algorithm::Symbol = DEFAULT_LYAPUNOV_ALGORITHM, caching::Bool = DEFAULT_CACHING, - use_workspaces::Bool = DEFAULT_USE_WORKSPACES) + use_workspaces::Bool = DEFAULT_USE_WORKSPACES, + filter_free_shocks::Union{Nothing,AbstractMatrix{<:Real}} = nothing, + measurement_error_std::Union{Nothing,Real,AbstractVector{<:Real},AbstractMatrix{<:Real}} = nothing, + initial_state = DEFAULT_INITIAL_STATE) # @nospecialize # reduce compile time if !caching invalidate_cache_validity!(𝓂) end @@ -1397,13 +1549,21 @@ function plot_model_estimates!(𝓂::ℳ, sylvester_algorithm³ = (isa(sylvester_algorithm, Symbol) || length(sylvester_algorithm) < 2) ? sum(k * (k + 1) ÷ 2 for k in 1:𝓂.constants.post_model_macro.nPast_not_future_and_mixed + 1 + 𝓂.constants.post_model_macro.nExo) > DEFAULT_SYLVESTER_THRESHOLD ? DEFAULT_LARGE_SYLVESTER_ALGORITHM : DEFAULT_SYLVESTER_ALGORITHM : sylvester_algorithm[2], lyapunov_algorithm = lyapunov_algorithm) - warn_irrelevant_tol(tol, algorithm; needs_covariance = filter == :kalman) + warn_irrelevant_tol(tol, algorithm; needs_covariance = isnothing(filter_free_shocks) && filter == :kalman) gr_back, attributes, attributes_redux = setup_plot_attributes(plot_attributes) # write_parameters_input!(𝓂, parameters, verbose = verbose) - filter, smooth, algorithm, _, pruning, warmup_iterations = normalize_filtering_options(filter, smooth, algorithm, false, warmup_iterations) + filter_free_plot = !isnothing(filter_free_shocks) + if filter_free_plot + @assert algorithm ∈ [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] "`plot_model_estimates!` with positional shocks only supports perturbation algorithms (`:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`)." + smooth = false + pruning = algorithm ∈ (:pruned_second_order, :pruned_third_order) + filter = :filter_free + else + filter, smooth, algorithm, _, pruning, warmup_iterations = normalize_filtering_options(filter, smooth, algorithm, false, warmup_iterations) + end solve!(𝓂, parameters = parameters, @@ -1470,8 +1630,22 @@ function plot_model_estimates!(𝓂::ℳ, data_in_deviations = MacroModelling.missing_data_to_nan(data) end - data_in_deviations, _, _, informative_periods = trim_informative_sample(data_in_deviations; - require_informative_periods = true) + measurement_error_std_to_plot = nothing + if filter_free_plot + filter_free_output = filter_free_data_with_model(𝓂, + data_in_deviations, + filter_free_shocks, + measurement_error_std, + initial_state; + algorithm = algorithm, + warmup_iterations = warmup_iterations, + opts = opts) + data_in_deviations = filter_free_output.data_in_deviations + informative_periods = filter_free_output.period_range + else + data_in_deviations, _, _, informative_periods = trim_informative_sample(data_in_deviations; + require_informative_periods = true) + end presample_periods = normalize_presample_periods(presample_periods, size(data_in_deviations, 2)) x_axis = axiskeys(data,2)[informative_periods] @@ -1482,9 +1656,17 @@ function plot_model_estimates!(𝓂::ℳ, x_axis = x_axis[periods] - variables_to_plot, shocks_to_plot, standard_deviations, decomposition = filter_data_with_model(𝓂, data_in_deviations, Val(algorithm), Val(filter), warmup_iterations = warmup_iterations, smooth = smooth, opts = opts) + if filter_free_plot + variables_to_plot = filter_free_output.variables + shocks_to_plot = filter_free_output.shocks + standard_deviations = zeros(0, 0) + decomposition = zeros(size(variables_to_plot, 1), size(shocks_to_plot, 1) + 2, size(variables_to_plot, 2)) + measurement_error_std_to_plot = filter_free_output.measurement_error_std + else + variables_to_plot, shocks_to_plot, standard_deviations, decomposition = filter_data_with_model(𝓂, data_in_deviations, Val(algorithm), Val(filter), warmup_iterations = warmup_iterations, smooth = smooth, opts = opts) + end - if pruning + if pruning && !filter_free_plot decomposition[:,1:(end - 2 - pruning),:] .+= SSS_delta decomposition[:,end - 2,:] .-= SSS_delta * (size(decomposition,2) - 4) end @@ -1557,6 +1739,8 @@ function plot_model_estimates!(𝓂::ℳ, :decomposition => decomposition, :variables_to_plot => variables_to_plot[var_idx, :], :data_in_deviations => data_in_deviations, + :measurement_error_std => measurement_error_std_to_plot, + :filter_free_plot => filter_free_plot, :shocks_to_plot => shocks_to_plot, :reference_steady_state => reference_steady_state[var_idx], :variable_names => variable_names_display, @@ -1670,7 +1854,8 @@ function plot_model_estimates!(𝓂::ℳ, :tol, :label, :shocks, :shock_names, :variables, :variable_names, - :rename_dictionary, :forecast_periods, :forecast_data, :extended_x_axis]) + :rename_dictionary, :forecast_periods, :forecast_data, :extended_x_axis, + :measurement_error_std, :filter_free_plot]) annotate_tol_diff!(annotate_diff_input, model_estimates_active_plot_container) @@ -1734,6 +1919,18 @@ function plot_model_estimates!(𝓂::ℳ, color = data_color) end + if any(k -> !isnothing(k[:measurement_error_std]), model_estimates_active_plot_container) + StatsPlots.plot!(legend_plot, + [NaN], [NaN], + seriestype = :scatter, + label = "[+σ,-σ]", + markershape = :rect, + markersize = 8, + markeralpha = 0.18, + markercolor = data_color, + markerstrokewidth = 0) + end + sort!(joint_shocks, by = normalize_superscript) sort!(joint_variables, by = normalize_superscript) @@ -1751,6 +1948,7 @@ function plot_model_estimates!(𝓂::ℳ, for var in joint_variables not_zero_anywhere = false + has_data_anywhere = false for k in model_estimates_active_plot_container var_idx = findfirst(==(var), apply_custom_name.(k[:variable_names], Ref(Dict(k[:rename_dictionary])))) @@ -1759,16 +1957,26 @@ function plot_model_estimates!(𝓂::ℳ, if isnothing(var_idx) || not_zero_anywhere # If the variable or shock is not present in the current plot_container, # we skip this iteration. + if !isnothing(var_idx) + obs_axis = collect(axiskeys(k[:data],1)) + obs_symbols = obs_axis isa String_input ? obs_axis .|> Meta.parse .|> replace_indices : obs_axis + obs_symbols_display = [replace_indices_in_symbol.(apply_custom_name(v, Dict(k[:rename_dictionary]))) for v in obs_symbols] + has_data_anywhere = has_data_anywhere || (var ∈ string.(obs_symbols_display)) + end continue else if any(.!isapprox.(k[:variables_to_plot][var_idx, periods], 0, atol = eps(Float32))) not_zero_anywhere = not_zero_anywhere || true # break # If any irf data is not approximately zero, we set the flag to true. end + obs_axis = collect(axiskeys(k[:data],1)) + obs_symbols = obs_axis isa String_input ? obs_axis .|> Meta.parse .|> replace_indices : obs_axis + obs_symbols_display = [replace_indices_in_symbol.(apply_custom_name(v, Dict(k[:rename_dictionary]))) for v in obs_symbols] + has_data_anywhere = has_data_anywhere || (var ∈ string.(obs_symbols_display)) end end - if not_zero_anywhere + if not_zero_anywhere || has_data_anywhere push!(joint_non_zero_variables, var) else # If all irf data for this variable and shock is approximately zero, we skip this subplot. @@ -1920,6 +2128,7 @@ function plot_model_estimates!(𝓂::ℳ, var_indx = findfirst(==(var), apply_custom_name.(k[:variable_names], Ref(Dict(k[:rename_dictionary])))) if var ∈ string.(obs_symbols_display) && !isnothing(var_indx) + obs_row = indexin([var], string.(obs_symbols_display)) if common_axis == [] idx = 1:length(k[:x_axis]) else @@ -1928,12 +2137,24 @@ function plot_model_estimates!(𝓂::ℳ, # Use extended_combined_x_axis length for padding data_in_deviations = fill(NaN, length(extended_combined_x_axis)) - data_in_deviations[idx] = k[:data_in_deviations][indexin([var], string.(obs_symbols_display)), periods] + data_in_deviations[idx] = k[:data_in_deviations][obs_row, periods] # data_in_deviations[idx][1:k[:presample_periods]] .= NaN + plotted_observations = data_in_deviations .+ k[:reference_steady_state][var_indx] + + if !isnothing(k[:measurement_error_std]) + std_values = fill(NaN, length(extended_combined_x_axis)) + std_values[idx] = k[:measurement_error_std][obs_row, periods] + plot_observation_std_band!(p, + extended_combined_x_axis, + plotted_observations, + std_values, + pal[length(model_estimates_active_plot_container) + i]) + end + StatsPlots.plot!(p, extended_combined_x_axis, - data_in_deviations .+ k[:reference_steady_state][var_indx], + plotted_observations, label = "", color = pal[length(model_estimates_active_plot_container) + i] ) @@ -1952,15 +2173,30 @@ function plot_model_estimates!(𝓂::ℳ, var_indx = findfirst(==(var), apply_custom_name.(k[:variable_names], Ref(Dict(k[:rename_dictionary])))) if var ∈ string.(obs_symbols_display) && !isnothing(var_indx) + obs_row = indexin([var], string.(obs_symbols_display)) # Use extended_combined_x_axis length for padding data_in_deviations_padded = fill(NaN, length(extended_combined_x_axis)) - data_vals = k[:data_in_deviations][indexin([var], string.(obs_symbols_display)),:] + data_vals = k[:data_in_deviations][obs_row,:] data_vals[1:k[:presample_periods]] .= NaN data_in_deviations_padded[1:length(combined_x_axis)] = data_vals[periods] + + plotted_observations = data_in_deviations_padded .+ k[:reference_steady_state][var_indx] + + if !isnothing(k[:measurement_error_std]) + std_values_padded = fill(NaN, length(extended_combined_x_axis)) + std_vals = k[:measurement_error_std][obs_row, :] + std_vals[1:k[:presample_periods]] .= NaN + std_values_padded[1:length(combined_x_axis)] = std_vals[periods] + plot_observation_std_band!(p, + extended_combined_x_axis, + plotted_observations, + std_values_padded, + data_color) + end StatsPlots.plot!(p, extended_combined_x_axis, - data_in_deviations_padded .+ k[:reference_steady_state][var_indx], + plotted_observations, label = "", color = data_color ) diff --git a/src/get_functions.jl b/src/get_functions.jl index ed46623b7..c975efed8 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -4832,6 +4832,247 @@ period_me_std(me_std::AbstractVector, idx::AbstractVector{Int}, ::Int) = view(me period_me_std(me_std::AbstractMatrix, idx::AbstractVector{Int}, t::Int) = view(me_std, idx, t) +function materialize_filter_free_measurement_error_std(me_std::Nothing, + n_obs::Int, + nT_input::Int, + period_range::UnitRange{Int}, + ::Type{R}) where R <: Real + return nothing +end + +function materialize_filter_free_measurement_error_std(me_std::Real, + n_obs::Int, + nT_input::Int, + period_range::UnitRange{Int}, + ::Type{R}) where R <: Real + @assert isfinite(me_std) && me_std > zero(me_std) "`measurement_error_std` must be finite and positive." + return fill(convert(R, me_std), n_obs, length(period_range)) +end + +function materialize_filter_free_measurement_error_std(me_std::AbstractVector{<:Real}, + n_obs::Int, + nT_input::Int, + period_range::UnitRange{Int}, + ::Type{R}) where R <: Real + @assert length(me_std) == n_obs "`measurement_error_std` vector must have one entry per observable (got $(length(me_std)), expected $n_obs)." + @assert all(x -> isfinite(x) && x > zero(x), me_std) "`measurement_error_std` entries must be finite and positive." + return repeat(reshape(convert(Vector{R}, me_std), :, 1), 1, length(period_range)) +end + +function materialize_filter_free_measurement_error_std(me_std::AbstractMatrix{<:Real}, + n_obs::Int, + nT_input::Int, + period_range::UnitRange{Int}, + ::Type{R}) where R <: Real + @assert size(me_std) == (n_obs, nT_input) "`measurement_error_std` matrix must have dimensions (n_observables, n_periods) = ($n_obs, $nT_input); got $(size(me_std))." + @assert all(x -> isfinite(x) && x > zero(x), me_std) "`measurement_error_std` entries must be finite and positive." + return convert(Matrix{R}, me_std[:, period_range]) +end + + +function filter_free_state_path(::Val{:first_order}, + 𝐒::AbstractMatrix, + state::AbstractVector{<:AbstractVector{<:Real}}, + shocks::AbstractMatrix{T}, + nT::Int, + past_idx::Vector{Int}, + n_warm::Int, + nVars::Int) where T <: Real + 𝐒₁ = 𝐒 + R = promote_type(eltype(𝐒₁), eltype(shocks), eltype(state[1])) + variables = zeros(R, nVars, nT) + cur_state = convert(Vector{R}, state[1]) + + for t in 1:n_warm + ϵ = view(shocks, :, t) + cur_state = 𝐒₁ * vcat(cur_state[past_idx], ϵ) + end + + for t in 1:nT + ϵ = view(shocks, :, n_warm + t) + cur_state = 𝐒₁ * vcat(cur_state[past_idx], ϵ) + variables[:, t] = cur_state + end + + return variables +end + +function filter_free_state_path(::Val{:second_order}, + 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix}, + state::AbstractVector{<:Real}, + shocks::AbstractMatrix{T}, + nT::Int, + past_idx::Vector{Int}, + n_warm::Int, + nVars::Int) where T <: Real + 𝐒₁, 𝐒₂ = 𝐒 + R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(shocks), eltype(state)) + variables = zeros(R, nVars, nT) + cur_state = convert(Vector{R}, state) + + for t in 1:n_warm + ϵ = view(shocks, :, t) + aug = vcat(cur_state[past_idx], one(R), ϵ) + cur_state = 𝐒₁ * aug + 𝐒₂ * ℒ.kron(aug, aug) / R(2) + end + + for t in 1:nT + ϵ = view(shocks, :, n_warm + t) + aug = vcat(cur_state[past_idx], one(R), ϵ) + cur_state = 𝐒₁ * aug + 𝐒₂ * ℒ.kron(aug, aug) / R(2) + variables[:, t] = cur_state + end + + return variables +end + +function filter_free_state_path(::Val{:third_order}, + 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}, + state::AbstractVector{<:Real}, + shocks::AbstractMatrix{T}, + nT::Int, + past_idx::Vector{Int}, + n_warm::Int, + nVars::Int) where T <: Real + 𝐒₁, 𝐒₂, 𝐒₃ = 𝐒 + R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(𝐒₃), eltype(shocks), eltype(state)) + variables = zeros(R, nVars, nT) + cur_state = convert(Vector{R}, state) + + for t in 1:n_warm + ϵ = view(shocks, :, t) + aug = vcat(cur_state[past_idx], one(R), ϵ) + kaug = ℒ.kron(aug, aug) + cur_state = 𝐒₁ * aug + 𝐒₂ * kaug / R(2) + 𝐒₃ * ℒ.kron(kaug, aug) / R(6) + end + + for t in 1:nT + ϵ = view(shocks, :, n_warm + t) + aug = vcat(cur_state[past_idx], one(R), ϵ) + kaug = ℒ.kron(aug, aug) + cur_state = 𝐒₁ * aug + 𝐒₂ * kaug / R(2) + 𝐒₃ * ℒ.kron(kaug, aug) / R(6) + variables[:, t] = cur_state + end + + return variables +end + +function filter_free_state_path(::Val{:pruned_second_order}, + 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix}, + state::AbstractVector{<:AbstractVector{<:Real}}, + shocks::AbstractMatrix{T}, + nT::Int, + past_idx::Vector{Int}, + n_warm::Int, + nVars::Int) where T <: Real + 𝐒₁, 𝐒₂ = 𝐒 + R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(shocks), eltype(state[1]), eltype(state[2])) + variables = zeros(R, nVars, nT) + cur_state = [convert(Vector{R}, state[1]), convert(Vector{R}, state[2])] + + for t in 1:n_warm + cur_state = pruned_second_order_state_update(cur_state, collect(view(shocks, :, t)), past_idx, nVars, 𝐒₁, 𝐒₂) + end + + for t in 1:nT + cur_state = pruned_second_order_state_update(cur_state, collect(view(shocks, :, n_warm + t)), past_idx, nVars, 𝐒₁, 𝐒₂) + variables[:, t] = cur_state[1] + cur_state[2] + end + + return variables +end + +function filter_free_state_path(::Val{:pruned_third_order}, + 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}, + state::AbstractVector{<:AbstractVector{<:Real}}, + shocks::AbstractMatrix{T}, + nT::Int, + past_idx::Vector{Int}, + n_warm::Int, + nVars::Int) where T <: Real + 𝐒₁, 𝐒₂, 𝐒₃ = 𝐒 + R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(𝐒₃), eltype(shocks), eltype(state[1]), eltype(state[2]), eltype(state[3])) + variables = zeros(R, nVars, nT) + cur_state = [convert(Vector{R}, state[1]), convert(Vector{R}, state[2]), convert(Vector{R}, state[3])] + + for t in 1:n_warm + cur_state = pruned_third_order_state_update(cur_state, collect(view(shocks, :, t)), past_idx, nVars, 𝐒₁, 𝐒₂, 𝐒₃) + end + + for t in 1:nT + cur_state = pruned_third_order_state_update(cur_state, collect(view(shocks, :, n_warm + t)), past_idx, nVars, 𝐒₁, 𝐒₂, 𝐒₃) + variables[:, t] = cur_state[1] + cur_state[2] + cur_state[3] + end + + return variables +end + + +function filter_free_data_with_model(𝓂::ℳ, + data_in_deviations::KeyedArray, + shocks::AbstractMatrix{T}, + measurement_error_std::Union{Nothing,Real,AbstractVector{<:Real},AbstractMatrix{<:Real}} = nothing, + initial_state = DEFAULT_INITIAL_STATE; + algorithm::Symbol = :second_order, + warmup_iterations::Int = DEFAULT_WARMUP_ITERATIONS, + opts::CalculationOptions = merge_calculation_options()) where T <: Real + @assert algorithm ∈ [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] "`plot_model_estimates` with positional shocks only supports perturbation algorithms (`:first_order`, `:second_order`, `:pruned_second_order`, `:third_order`, `:pruned_third_order`)." + @assert warmup_iterations >= 0 "`warmup_iterations` must be non-negative." + + observables = get_and_check_observables(𝓂.constants.post_model_macro, data_in_deviations) + raw_data = missing_data_to_nan(collect(data_in_deviations)) + trimmed_data, _, _, period_range = trim_informative_sample(raw_data; require_informative_periods = true) + + n_obs = length(observables) + nT_input = size(data_in_deviations, 2) + n_warm = max(warmup_iterations - 1, 0) + nT_total = nT_input + n_warm + nT = size(trimmed_data, 2) + nExo = 𝓂.constants.post_model_macro.nExo + + @assert size(shocks, 1) == nExo "`shocks` must have one row per exogenous shock (got $(size(shocks, 1)), expected $nExo)." + @assert size(shocks, 2) == nT_total "`shocks` must have $(nT_total) columns: $nT_input data periods plus $n_warm filter-free warmup shock columns (got $(size(shocks, 2)))." + + constants_obj, SS_and_pars, 𝐒, state, solved = get_relevant_steady_state_and_state_update(Val(algorithm), 𝓂.parameter_values, 𝓂, opts = opts, estimation = true) + @assert solved "No solution for these parameters." + + nVars = 𝓂.constants.post_model_macro.nVars + if initial_state isa AbstractVector{<:Real} + if length(initial_state) == nVars + state_shift = state isa AbstractVector{<:AbstractVector{<:Real}} ? (length(state) == 1 ? zero(state[1]) : -state[2]) : -state + state = adjust_initial_state(initial_state, algorithm, nVars, state_shift, SS_and_pars[1:nVars]) + if algorithm == :first_order + state = [state] + end + end + elseif !isempty(initial_state) + if state isa AbstractVector{<:AbstractVector{<:Real}} + R_state = promote_type(eltype(eltype(state)), eltype(initial_state[1])) + state = [convert(Vector{R_state}, i <= length(initial_state) ? initial_state[i] : state[i]) for i in eachindex(state)] + else + R_state = promote_type(eltype(state), eltype(initial_state[1])) + state = convert(Vector{R_state}, initial_state[1]) + end + end + + visible_cols = isempty(period_range) ? Int[] : n_warm .+ collect(period_range) + aligned_shocks = shocks[:, vcat(1:n_warm, visible_cols)] + filter_free_surface = algorithm == :first_order ? 𝐒 : Tuple(𝐒) + variables = filter_free_state_path(Val(algorithm), filter_free_surface, state, aligned_shocks, nT, 𝓂.constants.post_model_macro.past_not_future_and_mixed_idx, n_warm, nVars) + visible_shocks = aligned_shocks[:, n_warm + 1:end] + R = promote_type(eltype(variables), eltype(shocks), Float64) + aligned_measurement_error_std = materialize_filter_free_measurement_error_std(measurement_error_std, n_obs, nT_input, period_range, R) + names = AxisKeys.NamedDims.dimnames(data_in_deviations) + trimmed_keyed_data = KeyedArray(trimmed_data; NamedTuple{names}((axiskeys(data_in_deviations, 1), axiskeys(data_in_deviations, 2)[period_range]))...) + + return (variables = variables, + shocks = visible_shocks, + measurement_error_std = aligned_measurement_error_std, + data_in_deviations = trimmed_keyed_data, + period_range = period_range) +end + + function filter_free_reduction_indices(past_idx::Vector{Int}, obs_indices::Vector{Int}) kept_rows = sort(unique(vcat(past_idx, obs_indices))) past_in_kept_raw = indexin(past_idx, kept_rows) diff --git a/test/runtests.jl b/test/runtests.jl index 9c6e1447d..80ceb978f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,6 +49,8 @@ elseif test_set == "plots_4" include("test_plots_4.jl") elseif test_set == "plots_5" include("test_plots_5.jl") +elseif test_set == "filter_free_plotting" + include("test_filter_free_plotting.jl") elseif test_set == "higher_order_1" include("test_higher_order_1.jl") elseif test_set == "higher_order_2" From bdd275a2750a6ac42675499e8512f20cfc41e906 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Mon, 8 Jun 2026 21:16:05 +0000 Subject: [PATCH 2/7] Add filter-free shock paths to plot_model_estimates documentation --- docs/src/plot_model_estimates.md | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/docs/src/plot_model_estimates.md b/docs/src/plot_model_estimates.md index 41c581035..a35e63e46 100644 --- a/docs/src/plot_model_estimates.md +++ b/docs/src/plot_model_estimates.md @@ -73,6 +73,28 @@ plot_shock_decomposition(FS2000, data) This produces the same output as `plot_model_estimates` with `shock_decomposition = true`, which is the default setting for first order, pruned second order, and pruned third order solution algorithms. +## Filter-Free Shock Paths + +Filter-free estimates sample the latent structural shocks directly. These shock paths can be passed as a positional argument after the data: + +```julia +filter_free_shocks = zeros(length(get_shocks(FS2000)), size(data, 2)) + +plot_model_estimates(FS2000, data, filter_free_shocks; + algorithm = :second_order) +``` + +The shock path must have one row per exogenous shock and one column per data period. If `warmup_iterations > 1`, prepend `warmup_iterations - 1` additional shock columns, matching the filter-free `get_loglikelihood` convention. The model-implied endogenous path is forward-simulated from the supplied shocks and shown next to the observations. The existing keyword `shocks = :all/:none/...` still selects which shock subplots are displayed. + +An optional fourth positional argument adds a transparent +/- 1 standard-deviation band around the observations: + +```julia +plot_model_estimates(FS2000, data, filter_free_shocks, 0.01; + algorithm = :second_order) +``` + +The standard-deviation input can be a scalar, a vector with one entry per observable, or a matrix with dimensions `(n_observables, n_periods)`. Matrix inputs allow the band width to vary by observable and period. The same positional forms are supported by `plot_model_estimates!` for overlays. + ## Compare Model Estimates with `plot_model_estimates!` The `plot_model_estimates!` function (note the exclamation mark `!`) adds additional model estimates to an existing plot created with `plot_model_estimates`, enabling direct comparison between different scenarios. Any input argument that affects the model's output (such as datasets, solution algorithm, parameter values, filtering methods, or smoothing options) can be varied to compare how these changes influence the estimates. See the respective subsections below (e.g., [Data](#data-required), [Filter](#filter), [Solution Algorithm](#solution-algorithm), [Parameter Values](#parameter-values)) for details on specific arguments. From 54c41896ded0630a0015576d776e6fdcf526aa74 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Mon, 8 Jun 2026 21:31:25 +0000 Subject: [PATCH 3/7] Add filter-free shock tests to functionality_tests and remove filter_free_plotting from runtests --- test/functionality_tests.jl | 285 ++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 - 2 files changed, 285 insertions(+), 2 deletions(-) diff --git a/test/functionality_tests.jl b/test/functionality_tests.jl index 666f99aae..91663782a 100644 --- a/test/functionality_tests.jl +++ b/test/functionality_tests.jl @@ -516,6 +516,291 @@ function functionality_test(m, m2; algorithm = :first_order, plots = true) forecast_periods = forecast_periods) end end + + # Test filter-free positional shocks + if algorithm in [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] + nExo = length(m.constants.post_model_macro.exo) + nT = size(data, 2) + zero_shocks = zeros(nExo, nT) + + obs_axis = collect(axiskeys(data_in_levels, 1)) + n_obs = length(obs_axis) + scalar_std = 0.01 + vector_std = fill(scalar_std, n_obs) + matrix_std = fill(scalar_std, n_obs, nT) + matrix_std[:, 2:2:end] .= 0.02 + + for presample_periods in [0, 3] + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, + algorithm = algorithm, + data_in_levels = false, + presample_periods = presample_periods) + end + + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, zero_shocks, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) + + for presample_periods in [0, 3] + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data, zero_shocks, + algorithm = algorithm, + data_in_levels = false, + presample_periods = presample_periods) + end + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, zero_shocks, + algorithm = algorithm, + data_in_levels = true) + + # filter-free + scalar measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = false) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data, zero_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = false) + + # filter-free + vector measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, vector_std, + algorithm = algorithm, + data_in_levels = false) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data, zero_shocks, vector_std, + algorithm = algorithm, + data_in_levels = false) + + # filter-free + matrix measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = false) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data, zero_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = false) + + # filter-free + measurement_error_std + data_in_levels + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, zero_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, zero_shocks, vector_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, zero_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true) + + # filter-free + measurement_error_std + forecast_periods + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data, zero_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data, zero_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) + end + + # Test Kalman-filtered shocks as positional shocks (first_order only) + if algorithm == :first_order + obs_axis = collect(axiskeys(data_in_levels, 1)) + n_obs = length(obs_axis) + nT = size(data_in_levels, 2) + scalar_std = 0.01 + vector_std = fill(scalar_std, n_obs) + matrix_std = fill(scalar_std, n_obs, nT) + matrix_std[:, 2:2:end] .= 0.02 + + clear_solution_caches!(m, algorithm) + + smoothed_shocks = Matrix(get_estimated_shocks(m, data_in_levels, + algorithm = algorithm, + filter = :kalman, + smooth = true, + data_in_levels = true)) + + clear_solution_caches!(m, algorithm) + + filtered_shocks = Matrix(get_estimated_shocks(m, data_in_levels, + algorithm = algorithm, + filter = :kalman, + smooth = false, + data_in_levels = true)) + + # kalman smoothed + data_in_levels + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, + algorithm = algorithm, + data_in_levels = true) + + # kalman smoothed + forecast_periods + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) + + # kalman filtered + data_in_levels + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, filtered_shocks, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, filtered_shocks, + algorithm = algorithm, + data_in_levels = true) + + # kalman smoothed + scalar measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman smoothed + vector measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, vector_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, vector_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman smoothed + matrix measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman filtered + scalar measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, filtered_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, filtered_shocks, scalar_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman filtered + vector measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, filtered_shocks, vector_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, filtered_shocks, vector_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman filtered + matrix measurement_error_std + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, filtered_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, filtered_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true) + + # kalman smoothed + matrix measurement_error_std + forecast_periods + clear_solution_caches!(m, algorithm) + + plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) + + clear_solution_caches!(m, algorithm) + + plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std, + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) + end end @testset "plot_solution" begin diff --git a/test/runtests.jl b/test/runtests.jl index 80ceb978f..9c6e1447d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,8 +49,6 @@ elseif test_set == "plots_4" include("test_plots_4.jl") elseif test_set == "plots_5" include("test_plots_5.jl") -elseif test_set == "filter_free_plotting" - include("test_filter_free_plotting.jl") elseif test_set == "higher_order_1" include("test_higher_order_1.jl") elseif test_set == "higher_order_2" From 824babf1ce216da0cbcbbf290e454ed899ca0418 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Mon, 8 Jun 2026 21:42:57 +0000 Subject: [PATCH 4/7] Remove unused has_observed_series function and simplify plot_model_estimates logic --- ext/StatsPlotsExt.jl | 25 +++---------------------- 1 file changed, 3 insertions(+), 22 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index add77bcc5..4b67eef02 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -382,12 +382,6 @@ function plot_observation_std_band!(p, return p end -function has_observed_series(variable_idx::Int, obs_idx::AbstractVector{<:Integer}, data_in_deviations, periods) - obs_pos = findfirst(==(variable_idx), obs_idx) - isnothing(obs_pos) && return false - return any(isfinite, vec(data_in_deviations[obs_pos, periods])) -end - function process_rename_dictionary(rename_dictionary::AbstractDict, 𝓂::ℳ) relevant_keys = [k for k in keys(rename_dictionary) if (k isa String ? replace_indices(k) : k) in vcat(𝓂.constants.post_model_macro.var, 𝓂.constants.post_model_macro.exo)] |> sort processed = Any[] @@ -1025,8 +1019,7 @@ function plot_model_estimates(𝓂::ℳ, plot_count = 1 for v in var_idx - if all(isapprox.(variables_to_plot[v, periods], 0, atol = eps(Float32))) && - !has_observed_series(v, obs_idx, data_in_deviations, periods) + if all(isapprox.(variables_to_plot[v, periods], 0, atol = eps(Float32))) n_subplots -= 1 end end @@ -1058,8 +1051,7 @@ function plot_model_estimates(𝓂::ℳ, continue end else - if !(all(isapprox.(variables_to_plot[var_idx[i],periods], 0, atol = eps(Float32))) && - !has_observed_series(var_idx[i], obs_idx, data_in_deviations, periods)) + if !(all(isapprox.(variables_to_plot[var_idx[i],periods], 0, atol = eps(Float32)))) SS = reference_steady_state[var_idx[i]] if shock_decomposition @@ -1948,7 +1940,6 @@ function plot_model_estimates!(𝓂::ℳ, for var in joint_variables not_zero_anywhere = false - has_data_anywhere = false for k in model_estimates_active_plot_container var_idx = findfirst(==(var), apply_custom_name.(k[:variable_names], Ref(Dict(k[:rename_dictionary])))) @@ -1957,26 +1948,16 @@ function plot_model_estimates!(𝓂::ℳ, if isnothing(var_idx) || not_zero_anywhere # If the variable or shock is not present in the current plot_container, # we skip this iteration. - if !isnothing(var_idx) - obs_axis = collect(axiskeys(k[:data],1)) - obs_symbols = obs_axis isa String_input ? obs_axis .|> Meta.parse .|> replace_indices : obs_axis - obs_symbols_display = [replace_indices_in_symbol.(apply_custom_name(v, Dict(k[:rename_dictionary]))) for v in obs_symbols] - has_data_anywhere = has_data_anywhere || (var ∈ string.(obs_symbols_display)) - end continue else if any(.!isapprox.(k[:variables_to_plot][var_idx, periods], 0, atol = eps(Float32))) not_zero_anywhere = not_zero_anywhere || true # break # If any irf data is not approximately zero, we set the flag to true. end - obs_axis = collect(axiskeys(k[:data],1)) - obs_symbols = obs_axis isa String_input ? obs_axis .|> Meta.parse .|> replace_indices : obs_axis - obs_symbols_display = [replace_indices_in_symbol.(apply_custom_name(v, Dict(k[:rename_dictionary]))) for v in obs_symbols] - has_data_anywhere = has_data_anywhere || (var ∈ string.(obs_symbols_display)) end end - if not_zero_anywhere || has_data_anywhere + if not_zero_anywhere push!(joint_non_zero_variables, var) else # If all irf data for this variable and shock is approximately zero, we skip this subplot. From a3cf81d68bc2f6ae4edf5a7d1862b274dfd2c648 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Tue, 9 Jun 2026 15:43:12 +0000 Subject: [PATCH 5/7] fix ! calls and rework tests for plot model estimates --- ext/StatsPlotsExt.jl | 40 +++- src/default_options.jl | 2 + test/functionality_tests.jl | 450 ++++++++++++++++-------------------- 3 files changed, 246 insertions(+), 246 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 4b67eef02..22f941f44 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -997,6 +997,7 @@ function plot_model_estimates(𝓂::ℳ, :variables_to_plot => variables_to_plot[var_idx, :], :data_in_deviations => data_in_deviations, :measurement_error_std => measurement_error_std_to_plot, + :filter_free_shocks => filter_free_shocks, :filter_free_plot => filter_free_plot, :shocks_to_plot => shocks_to_plot, :reference_steady_state => reference_steady_state[var_idx], @@ -1732,6 +1733,7 @@ function plot_model_estimates!(𝓂::ℳ, :variables_to_plot => variables_to_plot[var_idx, :], :data_in_deviations => data_in_deviations, :measurement_error_std => measurement_error_std_to_plot, + :filter_free_shocks => filter_free_shocks, :filter_free_plot => filter_free_plot, :shocks_to_plot => shocks_to_plot, :reference_steady_state => reference_steady_state[var_idx], @@ -1781,6 +1783,42 @@ function plot_model_estimates!(𝓂::ℳ, push!(annotate_diff_input, "Data" => ["#$i" for i in data_idx]) end + filter_free_shocks_idx = Int[] + + if haskey(diffdict, :filter_free_shocks) + unique_ffs = unique(map(x -> isnothing(x) ? nothing : collect(x), diffdict[:filter_free_shocks])) + + for init in diffdict[:filter_free_shocks] + normalized = isnothing(init) ? nothing : collect(init) + for (i,u) in enumerate(unique_ffs) + if u == normalized + push!(filter_free_shocks_idx,i) + continue + end + end + end + + push!(annotate_diff_input, "Filter-free shocks" => ["#$i" for i in filter_free_shocks_idx]) + end + + measurement_error_std_idx = Int[] + + if haskey(diffdict, :measurement_error_std) + unique_mes = unique(map(x -> isnothing(x) ? nothing : collect(x), diffdict[:measurement_error_std])) + + for init in diffdict[:measurement_error_std] + normalized = isnothing(init) ? nothing : collect(init) + for (i,u) in enumerate(unique_mes) + if u == normalized + push!(measurement_error_std_idx,i) + continue + end + end + end + + push!(annotate_diff_input, "Obs. error std" => ["#$i" for i in measurement_error_std_idx]) + end + annotate_rename_dict_diff!(annotate_diff_input, diffdict) # Determine common and combined x axis @@ -1847,7 +1885,7 @@ function plot_model_estimates!(𝓂::ℳ, :shocks, :shock_names, :variables, :variable_names, :rename_dictionary, :forecast_periods, :forecast_data, :extended_x_axis, - :measurement_error_std, :filter_free_plot]) + :filter_free_shocks, :measurement_error_std]) annotate_tol_diff!(annotate_diff_input, model_estimates_active_plot_container) diff --git a/src/default_options.jl b/src/default_options.jl index 8ccc061f4..81416838d 100644 --- a/src/default_options.jl +++ b/src/default_options.jl @@ -109,6 +109,8 @@ const DEFAULT_ARGS_AND_KWARGS_NAMES = Dict( :quadratic_matrix_equation_algorithm => "Quadratic Matrix Equation Algorithm", :sylvester_algorithm => "Sylvester Algorithm", :lyapunov_algorithm => "Lyapunov Algorithm", + :filter_free_shocks => "Filter-free shocks", + :measurement_error_std => "Obs. error std", ) # Turing distribution wrapper defaults diff --git a/test/functionality_tests.jl b/test/functionality_tests.jl index 91663782a..1c7000970 100644 --- a/test/functionality_tests.jl +++ b/test/functionality_tests.jl @@ -516,292 +516,252 @@ function functionality_test(m, m2; algorithm = :first_order, plots = true) forecast_periods = forecast_periods) end end + end - # Test filter-free positional shocks - if algorithm in [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order] - nExo = length(m.constants.post_model_macro.exo) - nT = size(data, 2) - zero_shocks = zeros(nExo, nT) - - obs_axis = collect(axiskeys(data_in_levels, 1)) - n_obs = length(obs_axis) - scalar_std = 0.01 - vector_std = fill(scalar_std, n_obs) - matrix_std = fill(scalar_std, n_obs, nT) - matrix_std[:, 2:2:end] .= 0.02 - - for presample_periods in [0, 3] - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, - algorithm = algorithm, - data_in_levels = false, - presample_periods = presample_periods) - end - - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data_in_levels, zero_shocks, - algorithm = algorithm, - data_in_levels = true) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, - algorithm = algorithm, - data_in_levels = false, - forecast_periods = 6) - - for presample_periods in [0, 3] - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data, zero_shocks, - algorithm = algorithm, - data_in_levels = false, - presample_periods = presample_periods) - end - - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data_in_levels, zero_shocks, - algorithm = algorithm, - data_in_levels = true) - - # filter-free + scalar measurement_error_std - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = false) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data, zero_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = false) - - # filter-free + vector measurement_error_std - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, vector_std, - algorithm = algorithm, - data_in_levels = false) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data, zero_shocks, vector_std, - algorithm = algorithm, - data_in_levels = false) - - # filter-free + matrix measurement_error_std - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = false) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data, zero_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = false) - - # filter-free + measurement_error_std + data_in_levels - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data_in_levels, zero_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = true) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data_in_levels, zero_shocks, vector_std, - algorithm = algorithm, - data_in_levels = true) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data_in_levels, zero_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true) - - # filter-free + measurement_error_std + forecast_periods - clear_solution_caches!(m, algorithm) - - plot_model_estimates(m, data, zero_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = false, - forecast_periods = 6) - - clear_solution_caches!(m, algorithm) - - plot_model_estimates!(m, data, zero_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = false, - forecast_periods = 6) - end - - # Test Kalman-filtered shocks as positional shocks (first_order only) - if algorithm == :first_order - obs_axis = collect(axiskeys(data_in_levels, 1)) - n_obs = length(obs_axis) - nT = size(data_in_levels, 2) - scalar_std = 0.01 - vector_std = fill(scalar_std, n_obs) - matrix_std = fill(scalar_std, n_obs, nT) - matrix_std[:, 2:2:end] .= 0.02 - - clear_solution_caches!(m, algorithm) - - smoothed_shocks = Matrix(get_estimated_shocks(m, data_in_levels, - algorithm = algorithm, - filter = :kalman, - smooth = true, - data_in_levels = true)) - - clear_solution_caches!(m, algorithm) - - filtered_shocks = Matrix(get_estimated_shocks(m, data_in_levels, - algorithm = algorithm, - filter = :kalman, - smooth = false, - data_in_levels = true)) + # filter-free positional shocks — zero shocks, data (deviations), no std + nExo = length(m.constants.post_model_macro.exo) + nT = size(data, 2) + zero_shocks = zeros(nExo, nT) - # kalman smoothed + data_in_levels - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks; + algorithm = algorithm, + data_in_levels = false, + presample_periods = 0) - plot_model_estimates(m, data_in_levels, smoothed_shocks, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks; + algorithm = algorithm, + data_in_levels = false, + presample_periods = 3) - clear_solution_caches!(m, algorithm) + # filter-free positional shocks — zero shocks, data_in_levels, no std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, zero_shocks; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates!(m, data_in_levels, smoothed_shocks, - algorithm = algorithm, - data_in_levels = true) + # filter-free positional shocks — zero shocks, forecast + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks; + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) - # kalman smoothed + forecast_periods - clear_solution_caches!(m, algorithm) + # filter-free positional shocks — zero shocks, plot_model_estimates! + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks; + algorithm = algorithm, + data_in_levels = false, + presample_periods = 0) - plot_model_estimates(m, data_in_levels, smoothed_shocks, - algorithm = algorithm, - data_in_levels = true, - forecast_periods = 6) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks; + algorithm = algorithm, + data_in_levels = false, + presample_periods = 3) - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, zero_shocks; + algorithm = algorithm, + data_in_levels = true) + + # filter-free positional shocks — zero shocks, scalar std + obs_axis = collect(axiskeys(data_in_levels, 1)) + n_obs = length(obs_axis) + scalar_std = 0.01 + vector_std = fill(scalar_std, n_obs) + matrix_std = fill(scalar_std, n_obs, nT) + matrix_std[:, 2:2:end] .= 0.02 - plot_model_estimates!(m, data_in_levels, smoothed_shocks, - algorithm = algorithm, - data_in_levels = true, - forecast_periods = 6) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = false) - # kalman filtered + data_in_levels - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = false) - plot_model_estimates(m, data_in_levels, filtered_shocks, - algorithm = algorithm, - data_in_levels = true) + # filter-free positional shocks — zero shocks, vector std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks, vector_std; + algorithm = algorithm, + data_in_levels = false) - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks, vector_std; + algorithm = algorithm, + data_in_levels = false) - plot_model_estimates!(m, data_in_levels, filtered_shocks, - algorithm = algorithm, - data_in_levels = true) + # filter-free positional shocks — zero shocks, matrix std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = false) - # kalman smoothed + scalar measurement_error_std - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = false) - plot_model_estimates(m, data_in_levels, smoothed_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = true) + # filter-free positional shocks — zero shocks, std, data_in_levels + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, zero_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, zero_shocks, vector_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates!(m, data_in_levels, smoothed_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, zero_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true) - # kalman smoothed + vector measurement_error_std - clear_solution_caches!(m, algorithm) + # filter-free positional shocks — zero shocks, matrix std, forecast + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data, zero_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) - plot_model_estimates(m, data_in_levels, smoothed_shocks, vector_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data, zero_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = false, + forecast_periods = 6) - clear_solution_caches!(m, algorithm) + # Kalman-filtered shocks — smoothed, no std + clear_solution_caches!(m, algorithm) + smoothed_shocks = Matrix(get_estimated_shocks(m, data_in_levels, + algorithm = algorithm, + filter = :kalman, + smooth = true, + data_in_levels = true)) - plot_model_estimates!(m, data_in_levels, smoothed_shocks, vector_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + filtered_shocks = Matrix(get_estimated_shocks(m, data_in_levels, + algorithm = algorithm, + filter = :kalman, + smooth = false, + data_in_levels = true)) - # kalman smoothed + matrix measurement_error_std - clear_solution_caches!(m, algorithm) + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + # Kalman smoothed — forecast + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks; + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) - plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks; + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) - # kalman filtered + scalar measurement_error_std - clear_solution_caches!(m, algorithm) + # Kalman filtered — no std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, filtered_shocks; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates(m, data_in_levels, filtered_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, filtered_shocks; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + # Kalman smoothed — scalar std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates!(m, data_in_levels, filtered_shocks, scalar_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = true) - # kalman filtered + vector measurement_error_std - clear_solution_caches!(m, algorithm) + # Kalman smoothed — vector std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks, vector_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates(m, data_in_levels, filtered_shocks, vector_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks, vector_std; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + # Kalman smoothed — matrix std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates!(m, data_in_levels, filtered_shocks, vector_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true) - # kalman filtered + matrix measurement_error_std - clear_solution_caches!(m, algorithm) + # Kalman filtered — scalar std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, filtered_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates(m, data_in_levels, filtered_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, filtered_shocks, scalar_std; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + # Kalman filtered — vector std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, filtered_shocks, vector_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates!(m, data_in_levels, filtered_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, filtered_shocks, vector_std; + algorithm = algorithm, + data_in_levels = true) - # kalman smoothed + matrix measurement_error_std + forecast_periods - clear_solution_caches!(m, algorithm) + # Kalman filtered — matrix std + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, filtered_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true) - plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true, - forecast_periods = 6) + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, filtered_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true) - clear_solution_caches!(m, algorithm) + # Kalman smoothed — matrix std, forecast + clear_solution_caches!(m, algorithm) + plot_model_estimates(m, data_in_levels, smoothed_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) - plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std, - algorithm = algorithm, - data_in_levels = true, - forecast_periods = 6) - end - end + clear_solution_caches!(m, algorithm) + plot_model_estimates!(m, data_in_levels, smoothed_shocks, matrix_std; + algorithm = algorithm, + data_in_levels = true, + forecast_periods = 6) @testset "plot_solution" begin states = vcat(get_state_variables(m), m.constants.post_model_macro.past_not_future_and_mixed) From 09441fe16a75c0a814adce0a2b371c0db03e2ede Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Tue, 9 Jun 2026 21:36:45 +0000 Subject: [PATCH 6/7] mak data in scope for plot tets --- test/functionality_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/functionality_tests.jl b/test/functionality_tests.jl index 1c7000970..62ebd90e1 100644 --- a/test/functionality_tests.jl +++ b/test/functionality_tests.jl @@ -516,7 +516,6 @@ function functionality_test(m, m2; algorithm = :first_order, plots = true) forecast_periods = forecast_periods) end end - end # filter-free positional shocks — zero shocks, data (deviations), no std nExo = length(m.constants.post_model_macro.exo) @@ -762,6 +761,7 @@ function functionality_test(m, m2; algorithm = :first_order, plots = true) algorithm = algorithm, data_in_levels = true, forecast_periods = 6) + end @testset "plot_solution" begin states = vcat(get_state_variables(m), m.constants.post_model_macro.past_not_future_and_mixed) From bb8ca7da7839c2138be7ede088c9bf0927fc982e Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Wed, 10 Jun 2026 10:13:50 +0000 Subject: [PATCH 7/7] Refactor filter_free_state_path function signatures for improved type handling --- src/get_functions.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/get_functions.jl b/src/get_functions.jl index c975efed8..4c82559ce 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -4871,13 +4871,13 @@ end function filter_free_state_path(::Val{:first_order}, - 𝐒::AbstractMatrix, - state::AbstractVector{<:AbstractVector{<:Real}}, + 𝐒::AbstractMatrix{S}, + state::AbstractVector{<:AbstractVector{StateT}}, shocks::AbstractMatrix{T}, nT::Int, past_idx::Vector{Int}, n_warm::Int, - nVars::Int) where T <: Real + nVars::Int)::Matrix{promote_type(S, T, StateT)} where {S <: Real, T <: Real, StateT <: Real} 𝐒₁ = 𝐒 R = promote_type(eltype(𝐒₁), eltype(shocks), eltype(state[1])) variables = zeros(R, nVars, nT) @@ -4898,13 +4898,13 @@ function filter_free_state_path(::Val{:first_order}, end function filter_free_state_path(::Val{:second_order}, - 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix}, - state::AbstractVector{<:Real}, + 𝐒::Tuple{<:AbstractMatrix{S1},<:AbstractMatrix{S2}}, + state::AbstractVector{StateT}, shocks::AbstractMatrix{T}, nT::Int, past_idx::Vector{Int}, n_warm::Int, - nVars::Int) where T <: Real + nVars::Int)::Matrix{promote_type(S1, S2, T, StateT)} where {S1 <: Real, S2 <: Real, T <: Real, StateT <: Real} 𝐒₁, 𝐒₂ = 𝐒 R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(shocks), eltype(state)) variables = zeros(R, nVars, nT) @@ -4927,13 +4927,13 @@ function filter_free_state_path(::Val{:second_order}, end function filter_free_state_path(::Val{:third_order}, - 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}, - state::AbstractVector{<:Real}, + 𝐒::Tuple{<:AbstractMatrix{S1},<:AbstractMatrix{S2},<:AbstractMatrix{S3}}, + state::AbstractVector{StateT}, shocks::AbstractMatrix{T}, nT::Int, past_idx::Vector{Int}, n_warm::Int, - nVars::Int) where T <: Real + nVars::Int)::Matrix{promote_type(S1, S2, S3, T, StateT)} where {S1 <: Real, S2 <: Real, S3 <: Real, T <: Real, StateT <: Real} 𝐒₁, 𝐒₂, 𝐒₃ = 𝐒 R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(𝐒₃), eltype(shocks), eltype(state)) variables = zeros(R, nVars, nT) @@ -4958,13 +4958,13 @@ function filter_free_state_path(::Val{:third_order}, end function filter_free_state_path(::Val{:pruned_second_order}, - 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix}, - state::AbstractVector{<:AbstractVector{<:Real}}, + 𝐒::Tuple{<:AbstractMatrix{S1},<:AbstractMatrix{S2}}, + state::AbstractVector{<:AbstractVector{StateT}}, shocks::AbstractMatrix{T}, nT::Int, past_idx::Vector{Int}, n_warm::Int, - nVars::Int) where T <: Real + nVars::Int)::Matrix{promote_type(S1, S2, T, StateT)} where {S1 <: Real, S2 <: Real, T <: Real, StateT <: Real} 𝐒₁, 𝐒₂ = 𝐒 R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(shocks), eltype(state[1]), eltype(state[2])) variables = zeros(R, nVars, nT) @@ -4983,13 +4983,13 @@ function filter_free_state_path(::Val{:pruned_second_order}, end function filter_free_state_path(::Val{:pruned_third_order}, - 𝐒::Tuple{<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}, - state::AbstractVector{<:AbstractVector{<:Real}}, + 𝐒::Tuple{<:AbstractMatrix{S1},<:AbstractMatrix{S2},<:AbstractMatrix{S3}}, + state::AbstractVector{<:AbstractVector{StateT}}, shocks::AbstractMatrix{T}, nT::Int, past_idx::Vector{Int}, n_warm::Int, - nVars::Int) where T <: Real + nVars::Int)::Matrix{promote_type(S1, S2, S3, T, StateT)} where {S1 <: Real, S2 <: Real, S3 <: Real, T <: Real, StateT <: Real} 𝐒₁, 𝐒₂, 𝐒₃ = 𝐒 R = promote_type(eltype(𝐒₁), eltype(𝐒₂), eltype(𝐒₃), eltype(shocks), eltype(state[1]), eltype(state[2]), eltype(state[3])) variables = zeros(R, nVars, nT) @@ -5058,7 +5058,7 @@ function filter_free_data_with_model(𝓂::ℳ, visible_cols = isempty(period_range) ? Int[] : n_warm .+ collect(period_range) aligned_shocks = shocks[:, vcat(1:n_warm, visible_cols)] filter_free_surface = algorithm == :first_order ? 𝐒 : Tuple(𝐒) - variables = filter_free_state_path(Val(algorithm), filter_free_surface, state, aligned_shocks, nT, 𝓂.constants.post_model_macro.past_not_future_and_mixed_idx, n_warm, nVars) + variables = filter_free_state_path(Val(algorithm), filter_free_surface, state, aligned_shocks, nT, 𝓂.constants.post_model_macro.past_not_future_and_mixed_idx, n_warm, nVars)::Matrix{promote_type(T, Float64)} visible_shocks = aligned_shocks[:, n_warm + 1:end] R = promote_type(eltype(variables), eltype(shocks), Float64) aligned_measurement_error_std = materialize_filter_free_measurement_error_std(measurement_error_std, n_obs, nT_input, period_range, R)