Inference for Multivariate Stochastic Differential Equations with msde

Martin Lysy, JunYong Tong, Nigel Delaney

2019-01-11

Introduction

A \(d\)-dimensional stochastic differential equation (SDE) \(\boldsymbol{Y}_t = (Y_{1t},\ldots,Y_{dt})\) is written as \[ \mathrm{d}\boldsymbol{Y}_t = \boldsymbol{\Lambda}_{\boldsymbol{\theta}}(\boldsymbol{Y}_t)\,\mathrm{d}t + \boldsymbol{\Sigma}_{\boldsymbol{\theta}}(\boldsymbol{Y}_t)^{1/2}\,\mathrm{d}\boldsymbol{B}_t, \] where \(\boldsymbol{\Lambda}_{\boldsymbol{\theta}}(\boldsymbol{y})\) and \(\boldsymbol{\Sigma}_{\boldsymbol{\theta}}(\boldsymbol{y})\) are the drift and diffusion functions, and \(\boldsymbol{B}_t = (B_{1t},\ldots,B_{dt})\) is \(d\)-dimensional Brownian motion. The msde package implements a Markov Chain Monte Carlo (MCMC) algorithm to sample from the posterior distribution \(p(\boldsymbol{\theta}\mid \boldsymbol{Y})\) of the parameters given discrete observations \(\boldsymbol{Y}= (\boldsymbol{Y}_{0},\ldots,\boldsymbol{Y}_{N})\) recorded at times \(t_{0},\ldots,t_{N}\), with some of the \(d\) components of \(\boldsymbol{Y}_t\) possibly latent. To do this efficiently, msde requires on-the-fly C++ compiling of user-specified models. Instructions for setting up R to compile C++ code are provided in the Installation section.

Creating an sde.model object

The SDE model used throughout this vignette is the so-called Lotka-Volterra predator-prey model. Let \(H_t\) and \(L_t\) denote the number of Hare and Lynx at time \(t\) coexisting in a given habitat. The Lotka-Volterra SDE describing the interactions between these two animal populations is given by (Golightly and Wilkinson 2010): \[ \begin{bmatrix} \mathrm{d} H_t \\ \mathrm{d} L_t \end{bmatrix} = \begin{bmatrix} \alpha H_t - \beta H_tL_t \\ \beta H_tL_t - \gamma L_t \end{bmatrix}\, \mathrm{d} t + \begin{bmatrix} \alpha H_t + \beta H_tL_t & -\beta H_tL_t \\ -\beta H_tL_t & \beta H_tL_t + \gamma L_t\end{bmatrix}^{1/2} \begin{bmatrix} \mathrm{d} B_{1t} \\ \mathrm{d} B_{2t} \end{bmatrix}. \] Thus we have \(d = 2\), \(\boldsymbol{Y}_t = (H_t, L_t)\), and \(\boldsymbol{\theta}= (\alpha, \beta, \gamma)\).

The sdeModel class definition

In order to build this model in C++, we create a header file lotvolModel.h containing the class definition for an sdeModel object. The basic structure of this class is given below:

The meaning of each class member is as follows:

Thus the whole file lotvolModel.h is given below:

The additions to the previous code sections are:

  1. The header include guards (#ifndef/#define/#endif).
  2. The sdeModel:: identifier is prepended to the class member definitions when these are written outside of the class declaration itself.
  3. The inline keyword before the class member definitions, both of which ensure that only one instance of these functions is passed to the C++ compiler.

Compiling and checking the sde.model object

One the sdeModel class is created as the C++ level, it is compiled in R using the following commands:

## Loading required package: msde

Before using the model for inference, it is useful to make sure that the C++ entrypoints are error-free. To facilitate this, msde provides R wrappers to the internal C++ drift, diffusion, and validator functions, which can then be checked against R versions as follows:

# helper functions

# random matrix of size nreps x length(x) from vector x
jit.vec <- function(x, nreps) {
  apply(t(replicate(n = nreps, expr = x, simplify = "matrix")), 2, jitter)
}
# maximum absolute and relative error between two arrays
max.diff <- function(x1, x2) {
  c(abs = max(abs(x1-x2)), rel = max(abs(x1-x2)/max(abs(x1), 1e-8)))
}

# R sde functions

# drift and diffusion
lv.drift <- function(x, theta) {
  dr <- c(theta[1]*x[1] - theta[2]*x[1]*x[2], # alpha * H - beta * H*L
          theta[2]*x[1]*x[2] - theta[3]*x[2]) # beta * H*L - gamma * L
  dr
}
lv.diff <- function(x, theta) {
  df <- matrix(NA, 2, 2)
  df[1,1] <- theta[1]*x[1] + theta[2]*x[1]*x[2] # alpha * H + beta * H*L
  df[1,2] <- -theta[2]*x[1]*x[2] # -beta * H*L
  df[2,1] <- df[1,2] # -beta * H*L
  df[2,2] <- theta[2]*x[1]*x[2] + theta[3]*x[2] # beta * H*L + gamma * L
  chol(df) # always use sd scale in R
}

# validators
lv.valid.data <- function(x, theta) all(x > 0)
lv.valid.params <- function(theta) all(theta > 0)

# generate some test values
nreps <- 12
x0 <- c(H = 71, L = 79)
theta0 <- c(alpha = .5, beta = .0025, gamma = .3)
X <- jit.vec(x0, nreps)
Theta <- jit.vec(theta0, nreps)

# drift and diffusion check

# R versions
dr.R <- matrix(NA, nreps, lvmod$ndims) # drift
df.R <- matrix(NA, nreps, lvmod$ndims^2) # diffusion
for(ii in 1:nreps) {
  dr.R[ii,] <- lv.drift(x = X[ii,], theta = Theta[ii,])
  # flattens diffusion matrix into a row
  df.R[ii,] <- c(lv.diff(x = X[ii,], theta = Theta[ii,]))
}

# C++ versions
dr.cpp <- sde.drift(model = lvmod, x = X, theta = Theta)
df.cpp <- sde.diff(model = lvmod, x = X, theta = Theta)

# compare
max.diff(dr.R, dr.cpp)
## abs rel 
##   0   0
## abs rel 
##   0   0
##     x theta 
##  TRUE  TRUE

Simulating trajectories from the Lotka-Volterra model

The basis for both simulation and inference with SDEs is the Euler-Maruyama approximation (Maruyama 1955), which states that over a small time interval \(\Delta t\), the (intractable) transition density of the SDE can be approximated by \[ \boldsymbol{Y}_{t+\Delta t} \mid \boldsymbol{Y}_t \approx \mathcal{N}\Big(\boldsymbol{Y}_t + \boldsymbol{\Lambda}_\boldsymbol{\theta}(\boldsymbol{Y}_t)\Delta t, \boldsymbol{\Sigma}_\boldsymbol{\theta}(\boldsymbol{Y}_t)\Delta t\Big), \] with convergence to the true SDE dynamics as \(\Delta t\to 0\).

In order to simulate data from the Lotka-Volterra SDE model, we use the function sde.sim. Here we’ll generate \(N = 50\) observations of the process with initial values \(\boldsymbol{Y}_0 = (71, 79)\), and parameter values \(\boldsymbol{\theta}= (.5, .0025,.3)\), with time between observations of \(\Delta t= 1\) year. The dt.sim argument to sde.sim specifies the internal observation time used by the Euler-Maruyama approximation.

# simulation parameters
theta0 <- c(alpha = .5, beta = .0025, gamma = .3) # true parameter values
x0 <- c(H = 71, L = 79) # initial SDE values
N <- 50 # number of observations
dT <- 1 # time between observations (years)

# simulate data
lvsim <- sde.sim(model = lvmod, x0 = x0, theta = theta0,
                 nobs = N-1, # N-1 steps forward
                 dt = dT,
                 dt.sim = dT/100) # internal observation time
## Number of normal draws required: 4900
## Running simulation...
## Execution time: 0 seconds.
## Bad Draws = 0
# plot data
Xobs <- rbind(c(x0), lvsim$data) # include first observation
tseq <- (1:N-1)*dT # observation times
clrs <- c("black", "red")
par(mar = c(4, 4, 1, 0)+.1)
plot(x = 0, type = "n", xlim = range(tseq), ylim = range(Xobs),
     xlab = "Time (years)", ylab = "Population")
lines(tseq, Xobs[,"H"], type = "o", pch = 16, col = clrs[1])
lines(tseq, Xobs[,"L"], type = "o", pch = 16, col = clrs[2])
legend("topleft", legend = c("Hare", "Lynx"), fill = clrs)

Inference for multivariate SDE models

Parameter inference is conducted used a well-known data augmentation scheme due to Pedersen (1995). Assume that, as above, the SDE observations are evenly spaced with interobservation time \(\Delta t\). For any integer \(m > 0\), let \(\boldsymbol{Y}_{(m)} = (\boldsymbol{Y}_{m,0},\ldots,\boldsymbol{Y}_{m,Nm})\) denote the value of the SDE at equally spaced intervals of \(\Delta t_m = \Delta t/m\). Thus, \(\boldsymbol{Y}_{(1)} = \boldsymbol{Y}\) corresponds to the observed data, and for \(m > 1\), we have \(\boldsymbol{Y}_{m,nm} = \boldsymbol{Y}_n\). Thus we shall refer to the “missing data” as \(\boldsymbol{Y}_{\mathrm{miss}}= \boldsymbol{Y}_{(m)} \setminus \boldsymbol{Y}\). The Euler-Maruyama approximation to the complete likelihood is \[ \mathcal L(\boldsymbol{\theta}\mid \boldsymbol{Y}_{(m)}) = \prod_{n=0}^{Nm-1} \varphi\Big(\boldsymbol{Y}_{m,n+1} \mid \boldsymbol{Y}_{m,n} + \boldsymbol{\Lambda}_\boldsymbol{\theta}(\boldsymbol{Y}_{m,n})\Delta t_m, \boldsymbol{\Sigma}_\boldsymbol{\theta}(\boldsymbol{Y}_{m,n}\Delta t)\Big), \] where \(\varphi(\boldsymbol{y}\mid \mathbf \mu, \mathbf \Sigma)\) is the PDF of \(\boldsymbol{y}\sim \mathcal{N}(\mathbf \mu, \mathbf \Sigma)\). The Bayesian data augmentation scheme then consists of chosing a prior \(\pi(\boldsymbol{\theta})\) and sampling from the posterior distribution \[ p_m(\boldsymbol{\theta}\mid \boldsymbol{Y}) = \int \mathcal L(\boldsymbol{\theta}\mid \boldsymbol{Y}_{(m)}) \times \pi(\boldsymbol{\theta}) \, \mathrm{d}\boldsymbol{Y}_{\mathrm{miss}}. \] As \(m \to \infty\) this approximate posterior converges to the true SDE posterior \(p(\boldsymbol{\theta}\mid \boldsymbol{Y})\).

Posterior sampling from the Euler-Maruyama posterior is accomplished with the function sde.post. In the following example we use \(m = 1\), i.e., there is no missing data. We’ll use a Lebesgue prior \(\pi(\boldsymbol{\theta}) \propto 1\); more information on the default and custom prior specifications can be found in the following sections.

# initialize the posterior sampler
init <- sde.init(model = lvmod, x = Xobs, dt = dT,
                 m = 1, theta = c(.1, .1, .1))

nsamples <- 2e4
burn <- 2e3
lvpost <- sde.post(model = lvmod, init = init,
                   hyper = NULL, #prior specification
                   nsamples = nsamples, burn = burn)
## Output size:
## params = 60000
## Running posterior sampler...
## Execution time: 0.16 seconds.
## alpha accept: 44.3%
## beta accept: 44.7%
## gamma accept: 44%
# posterior histograms
tnames <- expression(alpha, beta, gamma)
par(mfrow = c(1,3))
for(ii in 1:lvmod$nparams) {
  hist(lvpost$params[,ii], breaks = 25, freq = FALSE,
       xlab = tnames[ii],
       main = parse(text = paste0("p[1](", tnames[ii], "*\" | \"*bold(Y))")))
  # superimpose true parameter value
  abline(v = theta0[ii], lwd = 4, lty = 2)
}

Missing data specification with sde.init

In the example above there was no missing data, i.e., \(m = 1\) and both components of the SDE are observed at each time \(t_{0},\ldots,t_{n}\). In order to refine the Euler-Maruyama approximation, we simply pass a larger value of \(m\) to sde.init:

We can also assume that only the first \(q < d\) components of the SDE are observed, with the last \(d-q\) being latent. In this case with \(q = 1\), the lynx population would be unobserved, and this is specified with the nvar.obs argument:

Note that the initial data x must still be supplied as an \((N+1) \times d\) matrix, with the missing values corresponding to initial values for the MCMC sampler.

Default prior specification

Since msde allows for some of the \(q\) components of \(\boldsymbol{Y}_t\) to be latent, a prior must be specified not only for \(\boldsymbol{\theta}\) but also for the latent variables in the initial observation \(\boldsymbol{Y}_0\).

In the example above, we assumed a Lebesgue prior \(\pi(\boldsymbol{\theta}, \boldsymbol{Y}_0) \propto 1\), with the restriction that \(\boldsymbol{\theta}, \boldsymbol{Y}_0 > 0\) (as specified in the sdeModel class definition via the isValidData and isValidParams validators).

Default prior

The default prior in msde is a multivariate normal, for which the (fixed) hyper-parameters are supplied via the hyper argument to sde.post. The hyper argument can either be NULL, or a list with elements mu and Sigma. These consist of a named vector named matrix specifying the mean and variance of the named elements. Unnamed elements are given a Lebesgue prior. So for example, posterior inference for the dataset above with \(m = 1\), latent variable \(L\), and prior distribution \[ L_0, \alpha, \gamma \stackrel{\mathrm{iid}}{\sim}\mathcal{N}(1, 1), \qquad \pi(H_0, \beta \mid H_0, \alpha, \gamma) \propto 1 \] is obtained as follows:

## Output size:
## params = 60000
## data = 2e+05
## Running posterior sampler...
## Execution time: 0.45 seconds.
## Bridge accept: min = 0.00909%, avg = 44.5%
## alpha accept: 44.1%
## beta accept: 43.5%
## gamma accept: 44.1%
## L0 accept: 43.9%

Custom prior specification

msde provides a two-stage mechanism for specifying user-defined priors. This is illustrated below with a simple log-normal prior on \(\boldsymbol{\eta}= (\alpha, \gamma, \beta, L_0) = (\eta _{1},\ldots,\eta _{4})\), namely \[ \pi(\boldsymbol{\eta}) \iff \log(\eta_i) \stackrel{\mathrm{iid}}{\sim}\mathcal{N}(\mu_i, \sigma_i^2), \] where the hyperparameters are \(\boldsymbol{\mu }= (\mu _{1},\ldots,\mu _{4})\) and \(\boldsymbol{\sigma }= (\sigma _{1},\ldots,\sigma _{4})\). Ultimately, the hyperparameters will be passed to sde.post as a two-element list, e.g.,

hyper <- list(mu = c(0,0,0,0), sigma = c(1,1,1,1))

The sdePrior class definition

The first step of the prior specification is to define it at the C++ level, through the sdPrior class. The class corresponding to the log-normal prior above is defined in the file lotvolPrior.h pasted below.

The meaning of each class member is as follows:

Formatting the R input to the C++ code

C++ is much less forgiving that R when it comes to accepting incorrect inputs. Thus, if a user accidently passed the hyperparameters to sde.post as e.g.,

at best this would cause garbage MCMC output and at worst, the R session to terminate abruptly. For this reason msde provides an input-checking mechanism, which can also used to format the hyperparameters into the list-of-numeric-vectors input expected by the C++ code. This is done by passing an appropriate input checking function to sde.make.model through the argument hyper.check. This argument accepts a function with the exact signature of the example below. In this example, hyper must be a list with elements mu and sigma, which are either:

  1. scalars, in which case each is replicated four times.
  2. vectors of length four, in which case the order is determined by \(\boldsymbol{\eta}= (\alpha, \beta, \gamma, L_0)\).

Compiling and checking the prior

Now we are ready to create the sde.model object:

We can also test the C++ implementation of the prior against one written in R, using the msde function sde.prior.

##          abs          rel 
## 2.842171e-14 8.830490e-17

Installation

The msde package requires on-the-fly C++ compiling of user-specified models which is handled through the R package Rcpp (Eddelbuettel and François 2011).

Enable C++ compiling for R

To make sure the C++ compiler is set up correctly, install the Rcpp package and from within R run the following:

If the code compiles and outputs 8.6 then the C++ compiler is interfaced with R correctly.

Optimize settings for the C++ compiler (optional)

It’s possible to speed up msde by a reasonable amount by passing a few flags to the C++ compiler. This can be done by creating a Makevars file (or Makevars.win on Windows). To do this, find your home folder by running the R command Sys.getenv("HOME"), and in that folder create a subfolder called .R containing the Makevars file (if it doesn’t exist already). Now add the following lines to th .R/Makevars file:

The first two options make the C++ code faster and the third uses the clang++ compiler instead of g++, which has better error messages (and is the default compiler on OS X).

Enable OpenMP support (optional)

On a multicore machine, msde can parallelize some of its computations with OpenMP directives. This can’t be done through R on Windows and is supported by default on recent versions of g++/clang++ on Linux. For OS X, the default version of clang++ does not support OpenMP but it is supported by that of the LLVM Project. This can be installed through Homebrew. After installing Homebrew using the instructions from the previous link, in Terminal run brew install llvm. Then have R use LLVM’s clang++ compiler by setting the following in .R/Makevars:

Note that these compiler directives alone will not enable OpenMP support. This happens at compile time by linking against -fopenmp, which is done internally by msde.

References

Eddelbuettel, D. and François, R., 2011. Rcpp: Seamless R and C++ integration. Journal of Statistical Software, 40 (8), 1–18.

Golightly, A. and Wilkinson, D.J., 2010. Discussion of ‘Particle Markov chain Monte Carlo methods’ by Christophe Andrieu, Arnaud Doucet, Roman Holenstein. JRSS B, 59 (2), 341–357.

Maruyama, G., 1955. Continuous Markov processes and stochastic equations. Rendiconti del Circolo Matematico di Palermo, 4 (1), 48–90.

Pedersen, A.R., 1995. A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations. Scandinavian Journal of Statistics, 22 (1), 55–71.