Appendix D — Structure of pre-built functions

Code

using Distributions
include("../src/loadData.jl")

D.1 Primary function

The provided code is built around a user-written function, denoted here bootstrapFilter(θ, Y, opts) . This function should encode the user’s entire state-space model by implementing a bootstrap filter, including initialisation and the projection-weighting-resampling steps.

This function should accept:

A parameter vector θ
Sorted dataframe of data Y
Options dictionary opts

and should return a tuple of the form (X, W), where:

X contains the particle values (typically in the form of an \(N \times T \times S\) vector where \(S\) is the number of hidden-states)
W is a \(N \times T\) matrix of observation weights

D.1.1 Example

Only minor modifications need to be made to the example in Section 5.1.1 in order to satisfy these requirements:

Code

σ = 0.1 # Model parameters
nzdata = loadData("NZCOVID") # Dataframe containing model data
opts = Dict() # A dictionary of parameter values
opts["T"] = length(nzdata.Ct)
opts["N"] = 10000
opts["pR0"] = Uniform(0, 10)

function runSimpleModel(σ, nzdata::DataFrame, opts::Dict)

    # Extract frequently used options
    T = opts["T"]
    N = opts["N"]

    # Initialise output matrices
    X = zeros(N, T)
    W = zeros(N, T)

    # Sample from initial distribution
    X[:,1] = rand.(opts["pR0"], N)

    # Run the filter
    for tt = 2:T

        # Project according to the state-space model
        X[:,tt] = exp.(rand.(Normal.(log.(X[:,tt-1]), σ)))

        # Weight according to the observation model
        Λ = sum(nzdata.Ct[tt-1:-1:1] .* ω[1:tt-1])
        W[:,tt] = pdf.(Poisson.(X[:,tt] .* Λ), nzdata.Ct[tt])

        # Resample
        inds = wsample(1:N, W[:,tt], N; replace=true)
        X[:, max(tt - L, 1):tt] = X[inds, max(tt - L, 1):tt]

    end

    return(X, W)

end

D.2 Pre-built functions

We provide a collection of functions that accept bootstrapFilter (and also θ, Y, and opts, depending on the function) as an argument:

Function	Description
`estimateLoglik()`	Calculates and returns the log-likelihood estimate (Equation 5.5). Includes optional argument `ignoreerror` (default false), that returns \(-\infty\) if the model returns an error.
`simplePMMH()`	A simple implementation of the PMMH algorithm. Returns a matrix of accepted values `C` and a matrix of sampled values with the corresponding likelihood estimate and accept/reject decision `OUT`.
`simplePMMHMulti()`	A multithreaded wrapper of `simplePMMH()`.
`PMMH()`	A (slightly) more comprehensive and adaptive implementation of the PMMH algorithm. Returns the same arguments, as well as a diagnostics dataframe.
`PMMHMulti()`	A multithreaded wrapper of `PMMH()`.