Getting Started with pfR

The Big Picture

In pf, the c++ library, each model/algorithm pair gets its own class. Writing this class’s methods defines a time series model, and the base class you choose to inherit from selects the particle filter that will be used. This typically requires writing two files: a c++ header file that declares the class, and a c++ source file that implements all the model specifics.

In pfR, the R package, there is one extra step. In addition to the header and source file, there is some additional c++ code that helps you to use everything from within an R session.

The beauty of this is that pf files are usable in pfR, and the majority of pfR code is usable in a c++ program. The downside to this is that R programmers will need a passing familiarity with c++–fortunately most of the c++ code is boilerplate, and pfR has a function that will write most of it for you.

Examples

Before we do any specific example, we’ll have to define a state space model, and choose some notation.

For \(t =1, \ldots, T\), let \(\{y_t\}\) be your observed time series, \(\{x_t\}\) be some set of hidden states, and (optional) \(\{z_t\}\) be a sequence of inputs you have at your disposal to help better predict each \(y_t\).

Defining a state-space model requires defining three kinds of distributions:

In the diagram below, circle nodes correspond with observed variables, square nodes represent hidden/latent variables, and arrows demonstrate how conditional dependence relationships (the dotted lines convey that they are optional).

Example 1: Univariate Stochastic Volatility with Leverage

Here’s a state space model from Harvey and Shephard (1996) Yu (2005) that can help us understand a segment of univariate financial returns data \(y_1, \ldots y_T\).

\[ \begin{gather*} y_t = \exp(x_t/2)\epsilon_t, \hspace{10mm} t = 1, \ldots, T \\ x_{t+1} = \mu + \phi(x_t - \mu) + \eta_t, \hspace{10mm} t = 1, \ldots, T-1 \\ x_1 \sim \mathcal{N}(\mu, \sigma^2 / (1-\phi^2)) \\ \begin{bmatrix} \epsilon_t \\ \eta_t \end{bmatrix} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}) \end{gather*} \]

\(\mathcal{N}(\cdot, \cdot)\) denotes a (possibly multivariate) Gaussian distribution, \(-1 < \rho < 1\) is the leverage parameter (typically negative for equities), and

\[ \boldsymbol{\Sigma} = \begin{bmatrix}1 & \rho\sigma \\ \rho\sigma & \sigma^2 \end{bmatrix}. \]

The interpretation of \(x_t\) is the (logarithm of) the volatility of the returns and our covariate will be \(z_t = y_{t-1}\). Writing this model in our notation yields

\(f(x_t \mid x_{t-1}, z_t, \theta)\) is univariate normal with mean \(\mu + \phi(x_t - \mu) + \rho\sigma y_{t-1}\exp[-x_t/2]\) and variance \(\sigma^2(1 - \rho^2)\).

Coding It Up

Suppose that the algorithm we want is a bootstrap filter with covariates. Our next step is writing all the c++ code for this model and this algorithm. Here is the bit of code that will do most of the work for you. The first argument is the prefix for all of your filenames and will be used throughout the rest of your work.

library(pfr)
createPFCPPTemplates("svol_leverage", "BSWC", fileDir = "./")

This will save three files to your desktop–svol_leverage_bswc.h, svol_leverage_bswc.cpp, and svol_leverage_bswc_export.cpp–and then open them in an RStudio session. You will notice some TODOs in the files. Consider these as requests for you to fill in the details of your model.

After you are finished editing them, they should like something like this, this, and this.

Compiling The Code

To compile the code you just wrote, simply run

svol_lev <- buildModelFuncs(".", "svol_leverage")

svol_lev is an Rcpp Module Eddelbuettel (2013) object that holds all the functions that allow you to use your model. You can call these functions like this:

svol_lev$svol_leverage_bswc_approx_LL(rnorm(100), c(.9, 0.0, 1.0, -.2))
svol_lev$svol_leverage_bswc_approx_filt(rnorm(100), c(.9, 0.0, 1.0, -.2))

Example 2: A Factor Multivariate Stochastic Volatility Model

Here’s another state space model from Chib, Omori, and Asai (2009) that can handle multivariate returns data \(\mathbf{y}_1, \mathbf{y}_2, \ldots\).

\[ \begin{gather*} \mathbf{y}_t = \mathbf{B} \mathbf{f}_t + \mathbf{V}_t^{1/2} \boldsymbol{\epsilon}_t, \hspace{10mm} \boldsymbol{\epsilon}_t \sim \mathcal{N}_p(\boldsymbol{0}, \mathbf{I}) \\ \mathbf{f}_t = \mathbf{D}_t^{1/2} \boldsymbol{\gamma}_t, \hspace{10mm} \boldsymbol{\gamma}_t \sim \mathcal{N}_q(\boldsymbol{0}, \mathbf{I}) \\ \mathbf{x}_{t+1} = \boldsymbol{\mu} + \boldsymbol{\Phi} (\mathbf{x}_{t} - \boldsymbol{\mu}) + \boldsymbol{\eta}_t, \hspace{10mm} \boldsymbol{\eta}_t \sim \mathcal{N}_{p+q}(\boldsymbol{0}, \boldsymbol{\Sigma} ) \\\\ \mathbf{x}_1 \sim \mathcal{N}(\boldsymbol{\mu}, \text{diag}(\sigma^2_1/(1-\phi_1^2), \ldots, \sigma^2_{p+q}/(1-\phi_{p+q}^2) ) ) \\ \end{gather*} \]

where \[ \begin{gather*} \mathbf{V}_t = \text{diag}[\exp(x_{1,t}), \ldots, \exp(x_{p,t}) ], \\ \mathbf{D}_t = \text{diag}[\exp(x_{p+1,t}), \ldots, \exp(x_{p+q,t}) ], \\ \boldsymbol{\Phi} = \text{diag}(\phi_1, \ldots, \phi_{p+q}), \\ \boldsymbol{\Sigma} = \text{diag}(\sigma^2_1, \ldots, \sigma^2_{p+q}) , \end{gather*} \]

\(p\) is the number of observations, and \(q\) is the number of factors.

Suppose you want to use an auxiliary particle filter Pitt and Shephard (1999) for this model. The c++ file templates are generated with

createPFCPPTemplates("mean_factor_msvol", "APF", fileDir = "./")

The finished files look like the ones location here, here and here. That code is compiled with

fsvol <- buildModelFuncs(".", "mean_factor_msvol")

and the (approximate) log-likelihood and filtering functions are run with

paramEst <- c(1.22, #B
              -1.17, -.01, .08, #mu
              .85, .84, .81, #Phi
              3.7, 3.4, 3.49) #Sigma
fakeData <- matrix(rnorm(50), ncol=2)
fsvol$mean_factor_msvol_apf_approx_LL(fakeData, paramEst)
fsvol$mean_factor_msvol_apf_approx_filt(fakeData, paramEst)

References

Chib, Siddhartha, Yasuhiro Omori, and Manabu Asai. 2009. “Multivariate Stochastic Volatility.” CIRJE F-Series, 365–400.
Eddelbuettel, Dirk. 2013. Seamless R and C++ Integration with Rcpp. New York: Springer. https://doi.org/10.1007/978-1-4614-6868-4.
Harvey, Andrew C., and Neil Shephard. 1996. “Estimation of an Asymmetric Stochastic Volatility Model for Asset Returns.” Journal of Business & Economic Statistics 14 (4): 429–34. https://doi.org/10.1080/07350015.1996.10524672.
Pitt, Michael K., and Neil Shephard. 1999. “Filtering via Simulation: Auxiliary Particle Filters.” Journal of the American Statistical Association 94 (446): 590–99. http://www.jstor.org/stable/2670179.
Yu, Jun. 2005. “On Leverage in a Stochastic Volatility Model.” Journal of Econometrics 127 (2): 165–78. https://doi.org/https://doi.org/10.1016/j.jeconom.2004.08.002.