Particle filter

.github/workflows/ci.yml

@@ -173,6 +173,10 @@ jobs:
             os: ubuntu-latest
             arch: x64
             test_set: "system_prior_estimation"
+          - version: '1'
+            os: ubuntu-latest
+            arch: x64
+            test_set: "particle_filter"
     steps:
       - uses: actions/checkout@v7
       - uses: julia-actions/setup-julia@v2

Project.toml

@@ -103,6 +103,7 @@ RuntimeGeneratedFunctions = "0.5"
 Showoff = "1"
 SparseArrays = "1"
 SpecialFunctions = "2"
+Statistics = "1"
 StatsPlots = "0.15"
 Subscripts = "0.1.3"
 Suppressor = "0.2"
@@ -135,10 +136,11 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Pigeons = "0eb8d820-af6a-4919-95ae-11206f830c31"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["ADTypes", "Aqua", "BenchmarkTools", "CondaPkg", "PythonCall", "JET", "Dates", "DelimitedFiles", "DifferentiationInterface", "DynamicPPL", "ForwardDiff", "Mooncake", "FlexiChains", "MCMCChains", "LineSearches", "Optim", "Test", "Turing", "Pigeons", "FiniteDifferences", "StatsPlots", "Preferences", "Zygote"]
+test = ["ADTypes", "Aqua", "BenchmarkTools", "CondaPkg", "PythonCall", "JET", "Dates", "DelimitedFiles", "Statistics", "DifferentiationInterface", "DynamicPPL", "ForwardDiff", "Mooncake", "FlexiChains", "MCMCChains", "LineSearches", "Optim", "Test", "Turing", "Pigeons", "FiniteDifferences", "StatsPlots", "Preferences", "Zygote"]

benchmark/particle_filter_llpf_comparison.jl

@@ -0,0 +1,104 @@
+# Benchmark: MacroModelling's particle filter vs LowLevelParticleFilters.jl
+#
+# Compares the log-likelihood and per-evaluation wall time of MacroModelling's
+# built-in particle filters (`filter = :particle`) against a bootstrap
+# `ParticleFilter` from LowLevelParticleFilters.jl, on the same first-order
+# (linear) DSGE state space, cross-checked against the exact Kalman likelihood.
+#
+# LowLevelParticleFilters and Distributions are NOT dependencies of the package;
+# run this in a throwaway environment, e.g.
+#
+#   julia --project=/tmp/pfbench -e '
+#       using Pkg; Pkg.develop(path="."); Pkg.add(["LowLevelParticleFilters","Distributions","AxisKeys","DelimitedFiles"])'
+#   julia --project=/tmp/pfbench benchmark/particle_filter_llpf_comparison.jl
+#
+# The DSGE structure is mapped onto LLPF's generic interface as
+#   xₜ = A·xₜ₋₁ + wₜ,   wₜ ~ N(0, B Bᵀ)      (rank-deficient process noise)
+#   yₜ = xₜ[observables] + vₜ,  vₜ ~ N(0, H)
+# with the initial cloud drawn from the ergodic covariance Σ (discrete Lyapunov).
+
+using MacroModelling
+using LowLevelParticleFilters
+using Distributions
+using Random, DelimitedFiles, AxisKeys
+import Statistics
+import LinearAlgebra as ℒ
+import MacroModelling: get_relevant_steady_state_and_state_update,
+                       particle_initial_state_covariance, merge_calculation_options
+
+# Pull the exact first-order state space (deviation form) the package filter uses.
+# The initial-state covariance comes from the package's own Lyapunov solver via
+# `particle_initial_state_covariance`, i.e. the very routine the particle filter
+# uses to spread its initial cloud, so both filters start from the same prior.
+function extract_linear(m, data_levels, observables, params)
+    constants, SS_and_pars, 𝐒, _, _ = get_relevant_steady_state_and_state_update(Val(:first_order), params, m)
+    T = constants.post_model_macro
+    nVars = T.nVars; nPast = T.nPast_not_future_and_mixed
+    ssnames = constants.post_complete_parameters.SS_and_pars_names
+    obs_idx = convert(Vector{Int}, indexin(observables, ssnames))
+    A = zeros(nVars, nVars); A[:, T.past_not_future_and_mixed_idx] .= 𝐒[:, 1:nPast]
+    B = Matrix(𝐒[:, nPast+1:end])
+    dev = collect(data_levels) .- SS_and_pars[obs_idx]
+    Σ, _ = particle_initial_state_covariance(m, T, merge_calculation_options(), :theoretical)
+    return A, B, obs_idx, Σ, dev, nVars
+end
+
+function bench_model(name, m, data, observables, me; N = 20000, nseed = 8)
+    params = m.parameter_values
+    kal = get_loglikelihood(m, data(observables), params; filter = :kalman,
+                            presample_periods = 0, initial_covariance = :theoretical,
+                            measurement_error = me .^ 2)
+    println("\n==== $name (N=$N) ====")
+    println("Kalman+ME = ", round(kal, digits = 3))
+
+    for pf_filter in (:bootstrap_particle, :tempered_particle)
+        Nn = pf_filter == :tempered_particle ? N ÷ 3 : N
+        t0 = time()
+        lls = [get_loglikelihood(m, data(observables), params; filter = pf_filter,
+                    algorithm = :first_order, presample_periods = 0, initial_covariance = :theoretical,
+                    measurement_error = me .^ 2,
+                    n_particles = Nn, particle_rng = Random.Xoshiro(s)) for s in 1:nseed]
+        println("MacroModelling ", rpad(String(pf_filter), 19), " N=$Nn  mean=", round(Statistics.mean(lls), digits = 2),
+                "  std=", round(Statistics.std(lls), digits = 2), "  time/run=", round((time() - t0) / nseed, digits = 3), "s")
+    end
+
+    A, B, obs_idx, Σ, dev, nVars = extract_linear(m, data(observables), observables, params)
+    me_var = (me isa AbstractVector ? collect(me) : fill(me, length(observables))) .^ 2
+    nObs = length(observables); nT = size(dev, 2)
+    df = MvNormal(zeros(nVars), ℒ.Symmetric(B * B') + 1e-10ℒ.I)
+    dg = MvNormal(zeros(nObs), ℒ.Diagonal(me_var))
+    d0 = MvNormal(zeros(nVars), ℒ.Symmetric(Σ) + 1e-10ℒ.I)
+    u = [Float64[] for _ in 1:nT]
+    y = [collect(dev[:, t]) for t in 1:nT]
+    t0 = time()
+    lls = map(1:nseed) do s
+        Random.seed!(s)
+        loglik(ParticleFilter(N, (x, u, p, t) -> A * x, (x, u, p, t) -> x[obs_idx], df, dg, d0), u, y)
+    end
+    println("LLPF           ", rpad("bootstrap", 10), " N=$N  mean=", round(Statistics.mean(lls), digits = 2),
+            "  std=", round(Statistics.std(lls), digits = 2), "  time/run=", round((time() - t0) / nseed, digits = 3), "s")
+end
+
+@model RBC2 begin
+    1 / c[0] = (β / c[1]) * (α * exp(z[1]) * k[0]^(α - 1) + (1 - δ))
+    c[0] + k[0] = (1 - δ) * k[-1] + q[0]
+    q[0] = exp(z[0]) * k[-1]^α * exp(g[0])
+    z[0] = ρz * z[-1] + std_z * eps_z[x]
+    g[0] = ρg * g[-1] + std_g * eps_g[x]
+end
+@parameters RBC2 begin
+    std_z = 0.01; std_g = 0.01; ρz = 0.4; ρg = 0.6; δ = 0.02; α = 0.5; β = 0.95
+end
+Random.seed!(12345)
+data_rbc = MacroModelling.simulate(RBC2, periods = 40)([:c, :q], :, :simulate)
+bench_model("RBC (2 obs)", RBC2, data_rbc, [:c, :q], 0.002; N = 20000, nseed = 8)
+
+dat, header = readdlm(joinpath(@__DIR__, "..", "test", "data", "usmodel.csv"), ',', header = true)
+dat = Float64.(dat); csv = vec(Symbol.(strip.(header)))
+dsw = KeyedArray(dat', Variable = csv, Time = axes(dat, 1))([:dy, :dc, :dinve, :labobs, :pinfobs, :dw, :robs], 47:230)
+obs_sw = [:dy, :dc, :dinve, :labobs, :pinfobs, :dwobs, :robs]
+dsw = rekey(dsw, :Variable => obs_sw)
+include(joinpath(@__DIR__, "..", "models", "Smets_Wouters_2007_linear.jl"))
+SS(Smets_Wouters_2007_linear, parameters = [:crhoms => 0.01, :crhopinf => 0.01, :crhow => 0.01, :cmap => 0.01, :cmaw => 0.01])
+bench_model("SW07 (7 obs)", Smets_Wouters_2007_linear, dsw, obs_sw,
+            2.0 .* [Statistics.std(collect(dsw(o))) for o in obs_sw]; N = 20000, nseed = 6)

docs/make.jl

@@ -48,6 +48,7 @@ makedocs(
             "Variance Decomposition" => "plot_conditional_variance_decomposition.md",
             "Model Estimates" => "plot_model_estimates.md",
         ],
+        "Filters" => "filters.md",
         "Steady State" => "steady_state.md",
         "Shapley decompositions (higher order)" => "shapley_decompositions.md",
         "Speed Benchmarks" => "speed.md",