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.
diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl
index 702cc0567..22f941f44 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,24 @@ 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 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 +652,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 +720,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 +768,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 +783,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 +859,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 +877,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 +903,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 +996,9 @@ 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_shocks => filter_free_shocks,
+                           :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,
@@ -1020,10 +1095,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 +1152,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 +1230,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 +1313,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 +1396,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 +1481,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 +1526,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 +1542,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 +1623,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 +1649,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 +1732,9 @@ 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_shocks => filter_free_shocks,
+                           :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,
@@ -1605,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
@@ -1670,7 +1884,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,
+         :filter_free_shocks, :measurement_error_std])
 
     annotate_tol_diff!(annotate_diff_input, model_estimates_active_plot_container)
     
@@ -1734,6 +1949,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)
 
@@ -1768,7 +1995,7 @@ function plot_model_estimates!(𝓂::ℳ,
             end
         end
         
-        if not_zero_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.
@@ -1920,6 +2147,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 +2156,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 +2192,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/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/src/get_functions.jl b/src/get_functions.jl
index ed46623b7..4c82559ce 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{S},
+                                state::AbstractVector{<:AbstractVector{StateT}},
+                                shocks::AbstractMatrix{T},
+                                nT::Int,
+                                past_idx::Vector{Int},
+                                n_warm::Int,
+                                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)
+    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{S1},<:AbstractMatrix{S2}},
+                                state::AbstractVector{StateT},
+                                shocks::AbstractMatrix{T},
+                                nT::Int,
+                                past_idx::Vector{Int},
+                                n_warm::Int,
+                                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)
+    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{S1},<:AbstractMatrix{S2},<:AbstractMatrix{S3}},
+                                state::AbstractVector{StateT},
+                                shocks::AbstractMatrix{T},
+                                nT::Int,
+                                past_idx::Vector{Int},
+                                n_warm::Int,
+                                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)
+    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{S1},<:AbstractMatrix{S2}},
+                                state::AbstractVector{<:AbstractVector{StateT}},
+                                shocks::AbstractMatrix{T},
+                                nT::Int,
+                                past_idx::Vector{Int},
+                                n_warm::Int,
+                                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)
+    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{S1},<:AbstractMatrix{S2},<:AbstractMatrix{S3}},
+                                state::AbstractVector{<:AbstractVector{StateT}},
+                                shocks::AbstractMatrix{T},
+                                nT::Int,
+                                past_idx::Vector{Int},
+                                n_warm::Int,
+                                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)
+    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)::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)
+    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/functionality_tests.jl b/test/functionality_tests.jl
index 666f99aae..62ebd90e1 100644
--- a/test/functionality_tests.jl
+++ b/test/functionality_tests.jl
@@ -516,6 +516,251 @@ function functionality_test(m, m2; algorithm = :first_order, plots = true)
                                             forecast_periods = forecast_periods)
                 end
             end
+
+        # 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)
+
+        clear_solution_caches!(m, algorithm)
+        plot_model_estimates(m, data, zero_shocks;
+                             algorithm = algorithm,
+                             data_in_levels = false,
+                             presample_periods = 0)
+
+        clear_solution_caches!(m, algorithm)
+        plot_model_estimates(m, data, zero_shocks;
+                             algorithm = algorithm,
+                             data_in_levels = false,
+                             presample_periods = 3)
+
+        # 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)
+
+        # 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)
+
+        # 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)
+
+        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)
+        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
+
+        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 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)
+        plot_model_estimates!(m, data, zero_shocks, vector_std;
+                              algorithm = algorithm,
+                              data_in_levels = false)
+
+        # 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)
+
+        clear_solution_caches!(m, algorithm)
+        plot_model_estimates!(m, data, zero_shocks, matrix_std;
+                              algorithm = algorithm,
+                              data_in_levels = false)
+
+        # 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)
+        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 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)
+
+        clear_solution_caches!(m, algorithm)
+        plot_model_estimates!(m, data, zero_shocks, matrix_std;
+                              algorithm = algorithm,
+                              data_in_levels = false,
+                              forecast_periods = 6)
+
+        # 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))
+
+        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))
+
+        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
+        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 — no std
+        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 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 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 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 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 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 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 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)
+
+        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
 
         @testset "plot_solution" begin