Particle filter - #310
Open
thorek1 wants to merge 25 commits into
Open
Conversation
… error Add `filter = :particle` to `get_loglikelihood` for all perturbation orders (:first_order through :pruned_third_order), integrating out the structural shocks by Monte Carlo. Three variants selectable via `particle_filter_algorithm`: :bootstrap (Dynare-style sequential-importance-resampling), :auxiliary (Pitt-Shephard), and :tempered (Herbst-Schorfheide). Includes a resampling suite (systematic/stratified/multinomial/residual), adaptive ESS resampling, Lyapunov-based initial cloud, and a seeded `rng` for reproducibility. Generalize measurement error (`measurement_error_std`) to the filter-based likelihood path, including the Kalman filter (forward pass in src and the ForwardDiff extension); it reduces exactly to the previous behaviour when zero. The particle filter is a stochastic, non-differentiable estimator and errors clearly under reverse-mode AD (Zygote/Mooncake) and ForwardDiff, steering users to gradient-free samplers (Pigeons, nested sampling). Validated that all three variants reproduce the exact Kalman log-likelihood on a small RBC model and on the Smets-Wouters (2007) linear model / US data, up to the expected Var/2 finite-particle bias. New test set `particle_filter` (test_particle_filter.jl + test_particle_filter_sw07.jl), wired into CI. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HhhiqdZNkzbGZKJ8oSeoTp
Rewrite the first-order (linear) path of all three particle-filter variants (bootstrap, auxiliary, tempered) for speed. Particles are stored as the columns of an nVars×N matrix and the whole swarm is propagated with two BLAS gemm calls (Xₜ = A·Xₜ₋₁ + B·Eₜ) instead of N type-unstable, per-call-allocating closure invocations. Particle pools are double-buffered so resampling and the tempered Metropolis mutation run fully in place; the mutation caches Base = A·Anc once per stage and only recomputes the batched shock term B·Eprop per proposal. Resampling uses preallocated index/cumulative-weight buffers (in-place `*_resample_indices!`), following LowLevelParticleFilters.jl; inverse measurement-error variances and the log-normaliser are cached. On the Smets-Wouters (2007) linear model (nVars=40, N=10k-20k) this cuts a tempered-filter likelihood evaluation from ~192.7M allocations / 11.75 GiB to ~630 allocations / ~20 MiB (≈300,000× fewer allocations), and speeds up the filters by roughly 11× (bootstrap/auxiliary) and 5-6× (tempered). All variants still reproduce the exact Kalman likelihood. Higher orders continue to use the generic methods. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HhhiqdZNkzbGZKJ8oSeoTp
Extend the performance work to the nonlinear orders (:second_order through
:pruned_third_order) for all three variants. Replace the type-unstable,
per-call-allocating `state_update` closure with typed in-place transitions
(`nonpruned_state_update_{2,3}!`, `pf_pruned_{2nd,3rd}!`) that gather the
augmented state with explicit loops and apply the solution matrices via
`kron!`/`mul!` with preallocated scratch. Particle pools are double-buffered and
built with a concretely-typed, `Val`-dispatched initialiser; the per-period loop
runs behind a function barrier so the large kwarg method body does not lose the
pool element type. A `Core.Box` (from a captured-and-reassigned pool variable in
the tempered filter) is avoided by keeping the captured cloud read-only and
swapping separate locals.
On Smets-Wouters (2007), pruned second order, N=4000 this cuts allocations from
~6.25M (bootstrap) / ~55M (tempered) to ~66k / ~153k, with GC time near zero.
All higher-order likelihood values are bit-identical to before. Values, missing
data, resampling schemes and error guards all verified; test_particle_filter
(47) and the SW07 test (7) pass.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HhhiqdZNkzbGZKJ8oSeoTp
Standalone benchmark (benchmark/particle_filter_llpf_comparison.jl) comparing MacroModelling's particle filters against a bootstrap ParticleFilter from LowLevelParticleFilters.jl on the same first-order DSGE state space, cross-checked against the exact Kalman likelihood. LLPF/Distributions are not package deps; the script header explains running it in a throwaway environment. Findings (RBC 2-obs / SW07 7-obs, N=20000): on the well-conditioned RBC all three (MacroModelling bootstrap, MacroModelling tempered, LLPF bootstrap) agree with the Kalman value; MacroModelling's specialised filter is ~2x faster than LLPF on RBC and ~15x faster on SW07 (0.25s vs 3.7s per evaluation), and its bootstrap tracks the Kalman likelihood closely on SW07 where the generic LLPF setup (rank-deficient DSGE process noise mapped to a jittered MvNormal) drifts. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HhhiqdZNkzbGZKJ8oSeoTp
Benchmark Results
Benchmark PlotsA plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR. |
Collapse `particle_filter_algorithm` into the `filter` argument: the filter is now fully identified by one symbol — `:bootstrap_particle`, `:auxiliary_particle` or `:tempered_particle` (with `:particle` kept as an alias for the bootstrap filter). A filter registry in `default_options.jl` drives validation and the internal variant dispatch. Prefix the particle-specific options so they read unambiguously alongside the other filter options: `particle_resampling`, `particle_resampling_threshold`, `particle_initial_state_scaling` and `particle_rng`. The filter bodies bind the historical short names once so the optimised hot loops are untouched. `measurement_error_std` now defaults to `:auto`, which resolves per filter: no measurement error for the Kalman and inversion filters (their previous behaviour), and 10% of each observable's sample standard deviation for the particle filters, which are degenerate without it. Raise the default particle count to 10_000, which keeps an SW07-sized problem accurate to a couple of log-likelihood points in well under a second per evaluation. Annotate the filter internals: why resampling is needed at all and how the four schemes trade variance against cost, what each stage and buffer of the bootstrap recursion does, and an intuitive account of what the auxiliary filter's look-ahead buys and why dividing the preview back out keeps it unbiased. The LLPF benchmark now takes its initial-state covariance from the package's own Lyapunov solver (via `particle_initial_state_covariance`) instead of a hand-rolled doubling loop, so both filters start from an identical prior. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Add a `measurement_error_covariance` keyword to `get_loglikelihood` for correlated measurement error, superseding `measurement_error_std` when supplied and validated as square, symmetric and positive definite. Nothing in a filter requires H to be diagonal — only H⁻¹ and log det H are ever needed. The Kalman filter forms its innovation covariance as a matrix anyway, so it now accepts an arbitrary H, in both the dense and the missing-data recursion (where the observed sub-block H[idx, idx] is used) and in the ForwardDiff extension. A diagonal covariance reproduces the per-observable standard deviations exactly. The particle filters keep the diagonal fast path: their inner loop is an elementwise quadratic form, and correlated measurement error in a DSGE is more naturally written into the model itself — measurement-error processes in the observation equations move the correlation into the state transition and leave H diagonal. An off-diagonal covariance is therefore rejected with a message pointing at both alternatives. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
New `docs/src/filters.md`, registered in the page tree, covering all filter options in one place: the state-space setup they share, a comparison table and a short decision rule, then the maths and assumptions of each filter — the Kalman recursion and why it is exact and differentiable, what the inversion filter inverts and why that rules out measurement error and requires as many shocks as observables, and the predict/weight/resample skeleton of the particle filters. Explains the properties that matter in practice: why the particle likelihood is unbiased but its logarithm is biased downward by about Var/2 (so it sits below the Kalman value and converges from below), why measurement error is required at all, what each particle variant buys, and how the resampling schemes trade variance against cost. Adds a section on how the filters map into one another — the particle filters converge to the Kalman likelihood on a linear model, the inversion and particle filters make opposite measurement-error assumptions, and correlated measurement error can be moved into the model's observation equations so that every filter can handle it. Replace the Dynare reference in the estimation tutorial with the primary sources (Gordon, Salmond & Smith 1993; Fernández-Villaverde & Rubio-Ramírez 2007; Pitt & Shephard 1999; Herbst & Schorfheide 2019) and update it to the new filter names. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…ints Add a `filter_data_with_model` method for the particle filters, so `get_model_estimates`, `get_estimated_variables`, `get_estimated_shocks` and the estimate plots work with `filter = :bootstrap_particle` and friends. It returns the filtered moments of the particle cloud: the weighted mean of the states, the weighted spread as standard deviations, and the weighted mean of the drawn shocks. All three particle variants target the same filtering distribution — they differ only in how efficiently they estimate the likelihood — so the moments come from the standard predict/weight/resample recursion whichever variant is selected. Two honest limitations, both signalled: these are filtered rather than smoothed estimates (a particle smoother is a different algorithm), and a linear shock decomposition does not exist for a nonlinear filter, so `decomposition` is returned as zeros with an informational message. Thread the particle options (`measurement_error_std`, `n_particles`, `particle_resampling`, `particle_resampling_threshold`, `particle_initial_state_scaling`, `particle_rng`) through those entry points and the two plotting functions that take a `filter`. Tested against the Kalman estimates on a linear model, where the filtered particle paths track the smoothed Kalman paths closely, plus nonlinear orders and the combined estimates entry point. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #310 +/- ##
==========================================
- Coverage 88.05% 88.03% -0.02%
==========================================
Files 35 37 +2
Lines 30243 32306 +2063
==========================================
+ Hits 26629 28440 +1811
- Misses 3614 3866 +252 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Implement smoothing for the particle filters. `smooth = true` now returns E[xₜ | y₁..T] instead of E[xₜ | y₁..ₜ] for `get_model_estimates`, `get_estimated_variables`, `get_estimated_shocks` and the estimate plots. The method is fixed-interval smoothing along the filter's genealogy: every particle surviving at T carries the ancestral line that produced it, and those lines are draws from p(x₁..T | y₁..T), so averaging them with the terminal weights gives the smoothed moments. The textbook backward-kernel smoother is not usable here: it reweights by p(xₜ₊₁ | xₜ), which for a DSGE with fewer shocks than states is a Dirac on a lower-dimensional manifold and hence undefined. The known limitation (path degeneracy of the ancestral lines) and the memory cost are documented. `normalize_filtering_options` now permits smoothing for the Kalman and particle filters and only disables it for the inversion filter. Verified on a linear model: smoothing more than halves the distance of the particle estimates to the Durbin-Koopman smoother (0.020 vs 0.047 filtered). Finish the argument cleanup: the last `particle_filter_algorithm` mention (a stale header comment in particle.jl) is gone, so the filter is identified by a single symbol everywhere. Replace the remaining Dynare references in the filter sources with the primary works — Gordon, Salmond & Smith (1993) and Fernández-Villaverde & Rubio-Ramírez (2007) for the bootstrap filter. All six filter-taking entry points now expose the same particle options, and a `measurement_error_std` passed to an entry point that cannot use it (the Kalman smoother path) reports that it is ignored instead of dropping it silently. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Shock decomposition works for the particle filters after all — it needs a shock
path, and the smoother now supplies one, so the same attribution the inversion
filter uses applies. At first order the contributions are additive and the split
is exact. At pruned second and third order they are not additive, which is what
the Aumann-Shapley (marginal contribution) attribution is for, so those orders
reuse `aumann_shapley_shock_decomposition_pruned_{2nd,3rd}_order!` from
inversion.jl under `marginal_contribution = true`. Non-pruned second/third order
have no decomposition at any filter, matching the existing gating. The previous
blanket "not defined for particle filters" claim was wrong and is gone.
One subtlety this surfaced: the Aumann-Shapley routine checks that the
contributions reproduce the supplied `variables`, but a smoothed *mean* is not a
model trajectory (averaging does not commute with a nonlinear transition), so
passing it in left a closure error the routine tried to remove by refining its
quadrature until it ran past its hand-coded node limit. The pruned decomposition
therefore attributes the trajectory implied by the smoothed shocks — the same
object the inversion filter decomposes — which closes exactly.
Also fixes a duplicated keyword block in `plot_model_estimates` introduced when
the particle options were threaded into the plotting extension, which made the
function uncallable with any arguments. Adds plotting coverage (Kalman, filtered
and smoothed particle filters, and a pruned second-order run with the
Aumann-Shapley decomposition) so the extension is exercised by the particle
filter test set rather than only by the plots jobs.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Implement the sequential pruned attribution so `marginal_contribution = false` works too, rather than reporting that Aumann-Shapley is required. It runs one trajectory per shock with only that shock switched on plus one with all of them: each single-shock path is that shock's contribution, the all-shock path minus their sum is the interaction the nonlinearity creates, and the remainder goes to the residual. This reproduces the inversion filter's column layout, so both attributions now match it exactly — nExo+3 columns for the sequential split and nExo+2 for the Aumann-Shapley split, which distributes the interaction across the shocks instead of isolating it. Decomposition works off whichever shock path the run produced, so it is available for the filtered estimates (`smooth = false`) as well as the smoothed ones (`smooth = true`); the two give different decompositions because the underlying shock estimates differ. Tests cover both attributions at both pruned orders under both settings, and check that at first order — where the split is additive — the shocks explain the bulk of the movement and the residual is left carrying the initial-state contribution. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…tics Reviewer feedback: the argument is a covariance, not a standard deviation, and it had three different names across the code (`measurement_error_std`, `measurement_error_variances`, `measurement_error_covariance`). - one user-facing kwarg `measurement_error`, subsuming the separate covariance argument: scalar = common variance, vector = per-observable variances, matrix = full covariance. The pre-existing positional `measurement_error_std` on the filter-free `get_loglikelihood` keeps its name, since it genuinely is a standard deviation. - kernels take `measurement_error` too, and `:auto` is resolved once at the user-facing layer (`resolve_measurement_error`) so no kernel sees a Symbol. - the particle filters now accept a correlated covariance: `DenseMeasurementError` caches a Cholesky factor of H restricted to each missing-data pattern. A diagonal matrix is reduced to the variance vector so the elementwise fast path is unchanged. - `on_failure_loglikelihood` defaults to -1e6 for the particle filters instead of -Inf; a stochastic failure should reject a proposal, not kill a chain. - `get_estimated_variable_standard_deviations` gains `algorithm`/`filter` and the particle kwargs, reporting the cloud spread at any perturbation order. - particle-filter keyword docs move to `common_docstrings.jl` as shared constants; `filters.md` and the estimation tutorial follow the new semantics. - answer the review questions inline: why the Kalman filter needs H at all (particle.jl header), why the initial covariance is first order at every perturbation order, and why the inversion filter's smoother is a no-op rather than unsupported (its filtered estimate is already the smoothed one). Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Reviewer feedback: the buffers were allocated per call. Inside a sampler the
likelihood is evaluated thousands of times at identical dimensions, so they
now live in `𝓂.workspaces.particle`, lazily sized by
`ensure_particle_workspace!(nVars, nExo, n_particles)` in the same style as the
Kalman and Lyapunov workspaces.
Six nVars×N and three nExo×N matrices plus the per-particle vectors cover the
simultaneous needs of all three first-order kernels (the tempered one needs the
most: ancestors, states, Metropolis proposals, and a swap partner for each).
Repeat evaluations at N = 20,000 now allocate 18 KiB instead of ~5 MB.
The higher-order kernels keep per-call pools: their particles are
`Vector{Vector}` whose element layout depends on the pruning order, so they do
not share a fixed buffer shape.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Reviewer feedback: the tempering arguments should be reachable from every
function that takes `filter`. They now are — but adding them as inert kwargs
would have been worse than not having them, so `filter_data_with_model` learns
the tempered recursion rather than only the bootstrap one.
`:bootstrap_particle` and `:auxiliary_particle` genuinely need no distinction
here: the auxiliary filter's look-ahead proposal changes the variance of the
*likelihood* estimate, not the cloud it leaves behind. `:tempered_particle`
does change the cloud — within each period it bridges from the prior to the
full measurement density, resampling and rejuvenating the shocks by
random-walk Metropolis at each stage — so the estimates and the smoother both
benefit from the extra distinct support points.
The smoother composes the within-period resampling maps with the
end-of-period ones (`within[t]` then `parent[t-1]`) so the genealogy still
walks correctly. Verified: on the linear model the tempered smoothed path is
closer to the Kalman smoother than the filtered path, and the tempering
controls demonstrably change the result.
Also documents why forward-filtering backward-smoothing is not available:
re-pairing a stored x_{t+1} with a different ancestor needs an ε solving
g(x_t^i, ε) = x_{t+1}, which is overdetermined when there are fewer shocks
than states, so every backward weight is zero. FFBS would need a
kernel-regularised transition; the genealogy smoother is exact as written, and
tempering is the lever against its path degeneracy.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…urement error Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
The first cut looked up the Cholesky factor through a `Dict` keyed on the observed-row pattern, once per particle per period, and solved against a view of the scratch buffer. Both allocate: 87 MB per evaluation at N = 20,000 over 40 periods. The pattern is constant inside a period, so a one-entry memo guarded by an allocation-free comparison skips the dictionary entirely, and the triangular solve is a hand-written forward substitution over a plain matrix — for the handful of observables a DSGE has, faster than dispatching to BLAS and with no view to heap-allocate. `vᵀH⁻¹v = ‖L⁻¹v‖²`, so the substitution produces the solve and the quadratic form in one pass. 87 MB → 20 KiB, same likelihood. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
The change of variables y -> eps contributes -log|det Z| (Z = CB) to the loglikelihood, but `logabsdets` was placed inside the -1/2 alongside the quadratic form, so it contributed -1/2 log|det Z|. This is not a harmless constant: log|det Z| carries the shock standard deviations, so halving it halves the likelihood's penalty against large shocks. On an AR(1) the profile likelihood peaks at sig^2 = 2Q/T instead of Q/T — estimated shock standard deviations inflate by exactly sqrt(2), confirmed numerically at 1.4151 vs sqrt(2) = 1.4142. The inversion filter is the default for every nonlinear algorithm, so this affected higher-order estimation throughout. Verified on an AR(1) observed directly, where the conditional likelihood is closed form: Kalman and inversion now both return 959.0753, matching exactly; before, inversion returned 270.6024, short by exactly half the Jacobian. Fixed in all ten primal branches (first order, pruned 2nd/3rd, missing-data variants) and the ten copies in the rrules, together with every hand-written Jacobian cotangent, which simply doubles when the primal term goes from -1/2 L to -L (sign conventions are preserved by scaling). The gradient cross-checks (ForwardDiff vs Zygote vs FiniteDifferences, all five orders including the under-identified case) stay at 264/264, so the primal and the pullbacks moved together. Adds test/test_inversion_filter_likelihood.jl, which pins the *level* against closed-form answers rather than only checking finiteness or cross-path agreement — the bug was invisible to both, since every path shared it. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
The sharp statement of how the filters relate is that the inversion filter *is* the Kalman filter started at P₁ = BB'. The inversion filter assumes x₀ is known exactly, so the only uncertainty about x₁ is that period's shocks; given at least as many observables as shocks the update then drives the posterior covariance to exactly zero and it stays there, so the two agree period by period rather than merely asymptotically. Verified to machine precision: identical on the small RBC, and 1.6e-9 apart on Smets-Wouters (2007) with 7 shocks, 7 observables and 184 periods. The ergodic prior is a genuinely different starting point — 489 log points away on SW07 — so the test also asserts that gap, guarding against passing for the trivial reason that the initial covariance does not matter. The particle-filter side of the same statement was already covered: the initial cloud is drawn with the ergodic covariance, Var(x₁) = AΣA' + BB' = Σ, which is `initial_covariance = :theoretical`. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
The particle filters' `initial_covariance` is Var(x₀) — the cloud is drawn around the initial state and *then* propagated — while the Kalman filter's argument of the same name is P₁ = Var(x₁), the first predicted state. They correspond as P₁ = A·Var(x₀)·A' + BB'. The distinction is invisible at the `:theoretical` default, because the ergodic covariance is the fixed point of exactly that map and so is carried to itself. That is why passing `:theoretical` to both filters lines them up, and why the difference went unnoticed. It bites as soon as an explicit matrix is supplied: reproducing a Kalman run with P₁ = BB' (i.e. the inversion filter) needs a *zero* matrix here, not BB'. Verified on SW07 at both ends of the correspondence, with measurement error and 20,000 particles: Var(x₀)=0 vs Kalman P₁=BB' agree to 1.55 log points, and Var(x₀)=BB' vs Kalman P₁=A BB' A'+BB' to 3.34 — both within Monte-Carlo error and on the expected (downward) side of it. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…ee path filters.md explained each filter but not the two inputs that cut across all of them, nor the conditions under which they are the same filter. Adds: - "Measurement error and the initial covariance": what each is and, more to the point, how they differ — H is noise on the observation and acts every period forever, P₁ is a prior on the state and decays. The three jobs H does (genuine measurement error, stochastic singularity, misspecification) are separated, along with the fourth, purely computational one for the particle filters. Includes the timing-convention warning (particle filters take Var(x₀), the Kalman filter Var(x₁)) and the interaction: H > 0 forces a strictly positive steady-state P, which is why the inversion filter's P ≡ 0 and measurement error are contradictory rather than merely unimplemented. - "When are they the same filter?": the equivalences as a table with their conditions. The sharp one is that the inversion filter *is* the Kalman filter started at P₁ = BB' — exact period by period given n_y ≥ n_ε and H = 0. - "More shocks than observables": spells out the implicit assumption the inversion filter makes there. Minimum norm is the conditional mean, so the score stays a proper density; what fails is the clamp P ≡ 0, which needs rank(CB) = n_ε. The filter then understates the innovation covariance, by an amount governed by how much of the unidentified subspace propagates into the next period's observables — measurable, and cheap, from a first-order Kalman recursion. - "The filter-free likelihood": the fourth option, which does not filter at all but treats the shocks as parameters. Notes why it has no initial covariance, why measurement error is mandatory there, and why its argument keeps the `_std` suffix (a matrix means per-period standard deviations, not a covariance). Verified by building the page: 4 tables, the admonition and both cross- references render, with no pipe collision in the table maths. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…alence `get_loglikelihood` accepted `initial_covariance` but the estimate entry points did not — they always used the ergodic covariance, with no way to override it. Adds the keyword to `get_shock_decomposition`, `get_estimated_shocks`, `get_estimated_variables`, `get_estimated_variable_standard_deviations` and the two estimate plots, threading it through `filter_data_with_model` into `filter_and_smooth`. It is forwarded only to the Kalman and particle filters; the inversion filter has no state covariance, and says so rather than silently ignoring a supplied matrix. The same pass also fixes two kwarg blocks that had been missed when the tempering controls were threaded through. With that in place the equivalence can be checked where it actually bites — on the paths, not on one scalar. With P₁ = BB' the Kalman gain is BB'C'(CBB'C')⁻¹ = B Z⁺, which *is* the inversion filter's state recursion, so the two must track the same states and shocks. On SW07 they agree to ~1e-10 across all variables and all 184 periods, against 0.99 relative deviation under the ergodic prior. One subtlety, now documented and tested: the states match the Kalman *smoothed* estimates, not the filtered ones. Seven observations do not pin all forty model variables contemporaneously, though the full sample does — checkable directly, since under P₁ = BB' the smoothed dispersion collapses to ~3.6e-5 while the filtered dispersion does not. Exact identification of the state is the inversion filter's assumption, so the smoothed estimates are what it reproduces. The shocks match under both. filters.md gains a "What you get by default" table and the state-equivalence section, plus a note that `initial_covariance` is expressed in different bases on the likelihood and estimate paths. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
Above first order there is no Kalman filter to compare against, so the
higher-order particle machinery — the pruned and non-pruned second- and
third-order transitions, their Kronecker scratch, the pruned Vector{Vector}
particle layout — had nothing exact to be checked against. Two references
close that gap.
1. A *linear* model has no higher-order solution terms, so every perturbation
order describes the same system and the exact answer is the Kalman
likelihood at every order. Smets-Wouters (2007) log-linearised is such a
model, and the tests assert that premise rather than assume it: the
inversion filter is deterministic, so it must return an identical value at
all five orders (it does, to 1e-10). Deviations from the Kalman value are
then -2.7 / -0.3 / -0.3 at first, pruned second and second order — all
within Monte-Carlo error and all on the expected downward side. Pruned and
non-pruned agree exactly, which is itself a check on the pruning code.
2. On a genuinely nonlinear model the reference is the inversion filter. As
H -> 0 the measurement density collapses onto the change of variables
y -> eps, so p(y|x) -> N(eps_hat;0,I)/|det Z|, which is the inversion
filter's per-period term; a degenerate initial cloud matches its other
assumption, that x0 is known. The test asserts the *direction* rather than a
tolerance, because the observed gap is non-monotonic (-35.9, -1.6, -3.2 as
the measurement-error variance falls through 1e-4, 1e-5, 1e-6). That is not
a defect: shrinking H is exactly what degenerates the importance weights, so
the approach stalls at a floor set by particle noise.
Third order on forty variables costs minutes per evaluation — the first version
of this took 19 minutes — so the third-order sweep runs on a small linear model
instead, where all five orders finish in seconds. SW07 keeps first and second
order, where it is cheap and the dimensionality is the point.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019YdYHujHfAMxP3jt9kwzzT
…elling.jl into particle-filter
* Keep higher-order state updates compressed * Cache higher-order power contractions * Fix compressed higher-order CI regressions * Cache invariant compressed selector maps * Fix compressed higher-order review paths * Correct review verification count * Fix ForwardDiff higher-order gradient paths * Rewrite the particle filters and add a guided (conditionally optimal) variant The estimates path was ~500x more allocating than the likelihood path and no particle filter gave shock estimates that survived a change of RNG seed on a Smets-Wouters-sized problem. Both are addressed here. Performance. The filters are rebuilt around a batched particle cloud: nVars x N matrices, one per pruned state component, so a period costs a handful of BLAS gemm calls instead of N gemv calls, and the estimates path shares that machinery instead of duplicating it in a closure that dispatched dynamically per particle. The block loop is split across Julia's threads. On SW07 at pruned second order get_estimated_shocks went from 12.06 s / 4900 MiB / 319M allocations to 0.35 s / 19 MiB / 36k for the bootstrap filter, and 92.8 s to 8.4 s for the tempered one. src/filter/particle.jl is shorter than before. Correctness. Batched compressed-Kronecker kernels match the reference vector kernels; the batched transition matches the model's own state_update at all five perturbation orders; on a linear model every order and every variant reproduces the Kalman likelihood, ordered by efficiency. Missing data, smoothing and both shock decompositions verified. Tempered filter. The Metropolis mutation is now preconditioned by the stage's own Gaussian covariance and its scale adapts towards 25% acceptance; a fixed isotropic step cannot suit a target that contracts as phi rises and is anisotropic across shocks. Its defaults move to tempering_target_ratio 1.5 and 4 mutation steps, which measurably improve both the estimates and the likelihood per unit of compute (see the comments in default_options.jl for the numbers). New filter :guided_particle. With about as many shocks as observables and a small measurement error, the observation nearly determines the shock given the ancestor, and that conditional is available in closed form. Drawing from it makes the importance weight identically constant when the transition is linear in the shock. The proposal's centre is refined by Gauss-Newton steps against the true residual, and the filter then anneals from the proposal to the exact conditional along q^(1-beta) pi^beta, which costs 1.2 stages per period on average and only spends more where the proposal is poor. On the euro-area data it needs eight to twelve times less compute than the tempered filter for equal accuracy or better, on both the estimates and the likelihood, at first and pruned second order. :particle now aliases :guided_particle rather than :bootstrap_particle, and tempering_mh_steps defaults to 2 for the guided filter and 4 for the others via a filter selector. This changes results for existing callers passing :particle. Also: the filters now warn instead of silently returning noise when the cloud has degenerated, and both AGENT_PROGRESS.md and tasks/todo.md record the measurements behind each choice, including the approaches that were tried and rejected (full adaptation, proposal over-dispersion, deferred within-period resampling). Verified with test/test_particle_filter.jl (142/142; its plotting testset needs a libgobject this container lacks) and a focused check covering all four variants at every perturbation order. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fa11GzkuyTSc6xjoEP17br * delete Agent progress * Fix higher-order inversion-filter pullbacks against finite differences The compressed-Kronecker migration left several reverse-mode paths on the uncompressed convention, or reusing buffers that were still live. Weights. d/dx compressed_kron²_power(aug(x)) is 2·compressed_kron²(aug, ∂aug/∂x) and the cubic analogue carries a 3, so the 1/2 and 1/6 in front of the 𝐒₂ and 𝐒₃ terms cancel down to 1 and 1/2 in the pullback. The non-pruned stochastic steady-state pullbacks kept the forward weights, which is the main regression here: the SSS correction is itself second-order small, so this read as 0.15 % on ∂state but cancelled up to ~390 % in the likelihood gradient, and it also reached get_steady_state's parameter Jacobian. The same off-by-a-factor appears in the cubic third-order warmup term, in accumulate_cubic_kron_jacobian_pullback! (whose helper is already the exact VJP), and in the second-order warmup, where a pair VJP was double-counting the 1/2. Aliasing. compressed_kron{²,³}_power_vjp! overwrite their output while the identity variants accumulate; three call sites passed a shared accumulator and so dropped the terms added before them. Separately, the dense pruned third-order path used the ∂kronstate¹⁻_vol cotangent as forward scratch, and the third-order with-missing backward loop read a workspace buffer still holding the last forward step's value. KKT blocks. The dense pruned second- and third-order shock-solver cotangents were assembled by uncompressed ℒ.kron!, which crashed on any model with nExo ≠ 2. Rewritten in the compressed pair/triple bases. Forward side. third_order_warmup_observation_and_jacobian still built its shock-state mixed terms uncompressed, so it disagreed with its own pullback; compressed_pair_hessian! sets both triangles, hence the factor 2 the shock solver's Newton step was missing; the ForwardDiff Dual overloads needed the same compressed treatment. FS2000 parameter gradients, Mooncake against central_fdm(5,1), worst relative error over T = 1/5/40 and all five algorithms: 8.5e-08, against 3.9 before. test_inversion_filter_gradients 264/264, test_higher_order_1 7869/7869, test_missing_data 128/128, test_rrule_robustness 57/57. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Address particle-filter review comments Naming. The bridging controls were named after tempering but are shared with the guided filter, so tempering_* becomes particle_* and the two per-variant values keep their variant's name (DEFAULT_GUIDED_MH_STEPS, DEFAULT_TEMPERED_MH_STEPS). The convention is now stated at the top of the defaults block. Magic numbers that were inline in particle.jl — scratch budget, block sizes, Metropolis adaptation targets and bounds, the low-ESS warning threshold — move to named defaults. propagate_block! loses its algorithm if/elseif chain in favour of dispatch on the perturbation order. build_guided_proposal reuses its factorisation across the Gauss-Newton refinement instead of re-solving. Comments that leaned on jargon are rewritten to say what the code does, and the defaults' commentary is cut back to the conclusions with the supporting tables dropped. Statistics was declared as a test dependency without a compat entry, which Aqua flags. JET could not correlate get_loglikelihood's measurement-error guard with the separate `if` that uses it, so the narrowing moves to the use site. Deletes the tasks/ scratch files, which were working notes rather than sources. test_particle_filter, Aqua and the JET hot-path suite all pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Correct the particle-filter guidance in the filters guide The page told readers to fall back to :tempered_particle "when the model is strongly nonlinear", which its own pruned-second-order benchmark contradicts, and to prefer it for estimates. Both were right when written, before the guided filter gained its own Metropolis sweeps; neither is now. Replaces the vague trigger with the assumption that actually binds — the guided proposal is built from the first-order shock impact on observables, so it wants the observation roughly linear in the shock — with the diagnostic that fires when it fails (the post-bridge ESS warning) and the price of switching (~10x). The estimates recommendation flips to :guided_particle, which the page's own table has at 0.078 seed dispersion against 0.093 at an eighth of the cost. The "estimates versus likelihoods" section now attributes the fix to mutation rather than to tempering, since both filters mutate. Also repairs a broken @ref. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Address review comments on the compressed higher-order work Correctness questions raised in review, both answered with finite differences and now pinned by tests. The compressed-power derivative identities are d/dx compressed_kron²_power(x) = 2·compressed_kron²(x, dx) and d/dx compressed_kron³_power(x) = 3·compressed_kron³(x, x, dx), so the forward Taylor weights 1/2 and 1/6 become 1 and 1/2 in every Jacobian built on them. That is why the third-order warmup recursion divides the cubic term by 2 where the uncompressed form divided by 6, and where the factor-of-3 change in occasionally_binding_constraints comes from. Checked against central differences at 1e-10; the alternatives are off by exactly 2 and 3. The OBC Jacobian was also checked end to end against differences of its own output (7.9e-10 at second order, 3.9e-10 at third). Structure. The Aumann-Shapley historical shock decomposition belongs to neither filter that uses it, so it moves to src/filter/decomposition.jl. The particle filter's three column-wise Kronecker kernels were a second implementation of what perturbation/solution.jl already does per vector; they become compressed_kron{²,³}_power_columns!/compressed_kron²_columns! next to the kernels they delegate to. The compressed_kron test set was orphaned — no CI job ever ran it — and now runs as part of basic. Performance, each measured rather than assumed. - The cubic shock/state index sets and their row maps were rebuilt per call (a sort plus a binary search per row) though they depend only on the model's dimensions. Memoised on those dimensions. - The KKT pair-block cotangents read the target column back on the right of a `.+=`, which materialises it: 240 B per call at nExo = 2 up to 328 kB at nExo = 40, and 1.6-5.8x slower than two `ger!`s. Rewritten as preallocated rank-1 updates. End to end it is a small share (0.5 MB of 150 MB on Smets-Wouters at nExo = 7) but it grows with nExo. - `A_mat .-= 2 .* Matrix(I(n_x))` in the ForwardDiff KKT blocks materialised a dense identity for a diagonal update. Now a loop. - `@simd` on the compressed kernels was tried and is 1.6-2.6x *slower* — the body's branch is perfectly predicted and `@simd` trades it for a masked store. Recorded so nobody retries it. - Densifying 𝐒₂/𝐒₃ for the particle filters was questioned. In the compressed basis they are 28-87% dense on every model measured, and dense `gemm` beats sparse over all of that range (3.2-7.9x with one BLAS thread, up to 26x with four); the crossover is near 1% density. Measurements are in the comment. The `filtered` keyword on solve_stochastic_steady_state_newton had no caller passing `false` — both live call sites pre-slice and pass `true` — so it and its dead branch are gone from all four methods. Documentation and naming. The filters guide now maps every knob to the algorithm step it controls, for the tempered and the guided filter alike. The bridging options are shared by both filters, not tempered-only as their docstrings implied, and now say so; `filter`'s docstring notes that the inversion filter is the fastest nonlinear option, that it can fail outright rather than degrade, and that its smoothed and filtered estimates coincide. The guided filter's internal `particle_mh_steps` fallback said 4 where the public selector gives it 2. Semicolon-joined statements are split one per line, cryptic locals in the particle filters are named for what they hold, and the de-indented lines the compressed migration left in inversion.jl are re-indented. Comments answering the rest: how the Laplace approximation feeds the guided proposal and why the reported likelihood does not depend on it, how the guided bridge differs from the tempered one, why the mutating filters keep last period's cloud when the others do not, and why the Gauss-Newton step is fine with fewer observables than shocks. test_compressed_kron 239/239 (including the new argument-order, derivative-weight and column-wise sets), test_inversion_filter_gradients 264/264, test_jet_hot_paths 302/302, test_filter_free_gradients 247/247, test_initial_state 187/187, test_particle_filter 146/146, test_missing_data 128/128, test_rrule_robustness 57/57, test_inversion_filter_likelihood 7/7. test_higher_order_1 was not run to completion locally; CI covers it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Take the cubic index sets from the model's constants The model already cached these. `third_order_indices` carries shock_state_state_idxs/_rows and shock_shock_state_idxs/_rows, filled once by ensure_conditional_forecast_constants!, and the inversion filter loops already read them; the memoised compressed_cubic_shock_maps added in the previous commit was a second, global, dict-and-lock cache of the same four vectors built a different way. It is gone. What is left is one construction — compressed_shock_state_state_index_map and its shock-shock sibling — used both to fill the model's constants and by the few pullback entry points that can be called without the model in scope, so the two cannot drift. The joint-warmup solver and its observation/Jacobian routine now take the sets as keyword arguments, mirroring what the rrules already did, and all ten call sites in the filters and their pullbacks pass the model's cached vectors down. The remaining fallback is lazy, where the old code hit the dict even when the caller had supplied everything. Dead code the sweep turned up, all of it left behind by the compressed migration: - second_order_warmup_observation_and_jacobian built the *cubic* index sets and never read them. That predates the memoisation, which only made the dead work cheap. - third_order_indices.I_exo2, a sparse I(nExo^2) built for every third-order model and read by nothing. The compressed shock pair is nExo(nExo+1)/2. - `II = I(n_exo_pair)` in all six inversion filter functions, assigned and never used, two of them allocating a sparse identity. The compressed path replaced kron(II, state_vol) with compressed_triple_state_to_pair!. `II` is still live in rrules.jl, where the uses are real. - compressed_shock_shock_state_indices, which had no callers. Also answered in a comment above DEFAULT_GUIDED_NEWTON_STEPS: how the guided proposal's Gauss-Newton step relates to the inversion filter. It is the same problem with the shock prior left in — the inversion filter solves r(eps) = 0, this maximises -0.5||eps||^2 - 0.5 r(eps)' H^-1 r(eps) — and the two coincide as the measurement error vanishes, since M = I + Bo'H^-1 Bo -> Bo'H^-1 Bo and K -> pinv(Bo). Checked: ||K r - pinv(Bo) r|| falls as O(sigma^2), 2e-1 / 2e-3 / 2e-7 at sigma^2 = 1e-2 / 1e-4 / 1e-8, for square Bo and both rectangular shapes. The comment also records why it does not inherit the inversion filter's hard failure (M is positive definite by construction, so the solve succeeds at any shape; the filter's own bail-out is cloud degeneracy, a different condition) and why running the inversion filter per particle instead would give up the shared factorisation that makes the proposal cheap. New testset pins the two builders against the sorted set plus the binary-search row map at four shapes. test_inversion_filter_gradients 264/264 (this is the suite that drives the third-order joint-warmup pullbacks), test_inversion_filter_likelihood 7/7, test_compressed_kron clean including the new set. Function-name diff on inversion.jl against the previous commit: nothing missing, nothing added. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Thore Kockerols <Thore.Kockerols@ecb.europa.eu> Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.