Demo

Sequential Monte Carlo

Sequential Monte Carlo methods tackle filtering problems without beyond the linear Gaussian cases approachable by a Kalman Filter. At the heart of this module is a simple framework for defining state space models beyond the limitations of a few cases. Moreover, I introduce a class of SMC algorithms that not only solve the filtering problem of latent state inference, but also parameter learning in a fully Bayesian framework.

Why Julia

Julia is a just in time (JIT) compiled language which generates machine code on the fly, all under a legible syntax which is deeply rooted in computational maths. The resulting script is not only fast, but reads similarly to pseudocode from an applied math journal.

Getting Started

Like other probablistic programming languages, SequentialMonteCarlo.jl begins with a model. There are various models already included, but this package also provides the flexibility to create your own. For starters, let's define a univariate log stochastic volatility model and simulate T=100 periods.

sv_mod = StateSpaceModel(
    StochasticVolatility(μ=-1.0,ρ=0.9,σ=0.2),
    (1,1)
)

T = 100
x,y = simulate(sv_mod,T)

This construction weighs heavily on a concept called multiple dispatch; where the same function can be defined differently given various types of inputs. Interestingly, this concept allows us to define a general function like transition(model::SSM,x) that has multiple meanings given a model. Thus allowing us to perform estimation over the entire parent class of models.

Particle Filtering

Suppose we want to infer the latent process in our log volatility model, specifically in a way that is tailored for online filtering. In SequentialMonteCarlo.jl, we can do this by employing a simple bootstrap filter; since, at each step, we generate a cloud of xs, we can collect summary statistics.

# define summary statistics
probs = [0.25,0.50,0.75]
xq    = fill(zeros(Float64,3),T)
logZ  = 0.0

# initialize filter
x,w,logμ = bootstrap_filter(sv_mod,y[1],N=1024)
xq[1]    = quantile(x,probs)
logZ    += logμ

# subsequent iterations
for t in 1:T
    logμ,w,ess = bootstrap_filter!(x,w,y[t],sv_mod)
    @printf("t = %4d\tess = %4.3f",t,ess)

    # update summary statistics
    xq[t] = quantile(x,probs)
    logZ += logμ
end

Since these naive methods are building blocks to more involved algorithms, function calls are made at each time step. Furthermore, running bootstrap_filter!(...) actually modifies the original memory of the input values, notably x and w, which performs updates without allocating a new return value. Instead the function just returns a collection of summary statistics over the updated cloud x.

If the user instead wants the log marginal likelihood p(y[1:T]) then the call is simplified by the following line of code:

logZ = log_likelihood(1024,y[1:T],model)

The construction of these particle filters is not perfect on its own, since it is meant specifically to perform joint estimation. However, they are fully functional and meant to be flexible enough when called in high volumes. Something that LowLevelParticleFilters.jl actually fails to do.

Joint Inference

Suppose a given state space model has an uncertain construction such that the model parameters are unobserved. Now the problem becomes twofold: what can we infer about the latent states as well as the parameters?

To solve this problem, I introduce two main algorithms: SMC² and density tempered SMC.

Running these algorithms requires a prior for the parameter space as well as a model constructor as we defined above with lg_mod(). Below we define lg_prior such that length(lg_prior) == length(θ) where θ is the input to the function lg_mod. Subsequently, this model specifically defines inference such that θ = [A,Q,R] taking B and ``

# model constructor
lg_mod(θ) = StateSpaceModel(
    LinearGaussian(A=θ[1],B=1.0,Q=θ[2],R=θ[3]),
    (1,1)
)

# prior definition
lg_prior = product_distribution([
    TruncatedNormal(0,1,-1,1),
    LogNormal(),
    LogNormal()
])

Upon declaration of the prior and model constructor, we define a generic SMC sampler with 513 parameter particles, 1024 state particles, and ESS threshold of 0.5, and 3 MCMC steps. To demonstrate, we can run density tempered SMC like so...

lg_dt_smc = SMC(512,1024,lg_mod,lg_prior,3,0.5)
density_tempered(lg_dt_smc,y)

For an online algorithm like SMC², we treat it similarly to the particle filter by using mutating functions to change properties of lg_smc² as time progresses.

lg_smc² = SMC(512,1024,lg_mod,lg_prior,3,0.5)
smc²(lg_smc²,y)

for t in 2:T
    smc²!(lg_smc²,y,t)
end

Data Generating Process

Now assume we have a simple stochastic volatility model in which we want to perform joint inference over latent states x[1:t] and model parameters θ = [μ,ρ,σ] given observations y[1:t]. For review, the data generating process can be summarized by the following state transition and observation processes.

x[t] ~ Normal(μ+ρ*(μ-x[t-1]),σ)
y[t] ~ Normal(0,exp(0.5*x[t]))

For demonstration, we simulated T = 1000 periods given μ = -1.0, ρ = 0.9, and σ = 0.2. The set up for SMC² is exactly like it was before using 1024 state particles, 512 parameter particles, and 3 MCMC steps.

sv_prior = product_distribution([
    Normal(0.0,2.0),
    Uniform(-1.0,1.0),
    Gamma(1.0,1.0)
])

sv_mod(θ) = StateSpaceModel(
    StochasticVolatility(θ...),
    (1,1)
)

sv_smc² = SMC(512,512,sv_mod,sv_prior,3,0.5)

For the addition of a visual, I implemented a function plot_histograms which returns the weighted histograms given the weighted sample produced by the algorithm. To specifically define the iteration, we run the following.

smc²(sv_smc²,y)
plt = @animate for t in 2:T
    smc²!(sv_smc²,y,t)
    plot_histograms(sv_smc²,["μ","ρ","σ"])
    plot!(plot_title=string("\n",t))
end every 10

This estimation also works for density tempered SMC, although I omit the plots for the sake of brevity. The following code summarizes this procedure:

sv_dt_smc = SMC(512,1024,sv_mod,sv_prior,3,0.5)
density_tempered(sv_dt_smc,y)

Application

Using consumer expenditures data from FRED, we can replicate the model used in Stock & Watson for inflation de-trending. Specifically, they use the unobserved component stochastic volatility model, which I have outlined below using notation consistent with this report.

x[t] ~ Normal(x[t-1],exp(0.5*σx[t-1]))
y[t] ~ Normal(x[t],exp(0.5*σy[t]))

σx[t] ~ Normal(σx[t-1],γ)
σy[t] ~ Normal(σy[t-1],γ)

For this model, we define the generator ucsv_mod(θ) where θ = [γ,x[0],(σx[0],σy[0])] using the following specification.

ucsv_mod(θ) = StateSpaceModel(
    UCSV(θ[1],θ[2],(θ[3],θ[4])),
    (3,1)
)

ucsv_prior = product_distribution([
    Uniform(0.0,1.0),   # γ
    Normal(3.0,2.0),    # x[0]
    Uniform(0.0,2.0),   # σx[0]
    Uniform(0.0,2.0)    # σy[0]
])

To deploy SMC² for online estimation, we use ucsv_smc² = SMC(8192,512,ucsv_mod,ucsv_prior,5,0.5), and iterate as we demonstrated before. After iteration concludes, we produce the following graphs for trend and cycle respectively: