msde: Bayesian Inference for Multivariate Stochastic Differential Equations

Description

Implements an MCMC sampler for the posterior distribution of arbitrary time-homogeneous multivariate stochastic differential equation (SDE) models with possibly latent components. The package provides a simple entry point to integrate user-defined models directly with the sampler's C++ code, and parallelizes large portions of the calculations when compiled with 'OpenMP'.

Details

See package vignettes; vignette("msde-quicktut") for a tutorial and vignette("msde-exmodels") for several example models.

Author(s)

Maintainer: Martin Lysy mlysy@uwaterloo.ca

Authors:

• Feiyu Zhu

• JunYong Tong

Other contributors:

• Trevor Kitt [contributor]

• Nigel Delaney [contributor]

Examples

# Posterior inference for Heston's model

# compile model
hfile <- sde.examples("hest", file.only = TRUE)
param.names <- c("alpha", "gamma", "beta", "sigma", "rho")
data.names <- c("X", "Z")
hmod <- sde.make.model(ModelFile = hfile,
                       param.names = param.names,
                       data.names = data.names)
# or simply load pre-compiled version
hmod <- sde.examples("hest")

# Simulate data
X0 <- c(X = log(1000), Z = 0.1)
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
dT <- 1/252
nobs <- 1000
hest.sim <- sde.sim(model = hmod, x0 = X0, theta = theta,
                    dt = dT, dt.sim = dT/10, nobs = nobs)

# initialize MCMC sampler
# both components observed, no missing data between observations
init <- sde.init(model = hmod, x = hest.sim$data,
                 dt = hest.sim$dt, theta = theta)

# Initialize posterior sampling argument
nsamples <- 1e4
burn <- 1e3
hyper <- NULL # flat prior
hest.post <- sde.post(model = hmod, init = init, hyper = hyper,
                      nsamples = nsamples, burn = burn)

# plot the histogram for the sampled parameters
par(mfrow = c(2,3))
for(ii in 1:length(hmod$param.names)) {
  hist(hest.post$params[,ii],breaks=100, freq = FALSE,
       main = parse(text = hmod$param.names[ii]), xlab = "")
}

Loglikelihood for multivariate Ornstein-Uhlenbeck process.

Description

Computes the exact Euler loglikelihood for any amount of missing data using a Kalman filter.

Usage

mou.loglik(X, dt, nvar.obs, Gamma, Lambda, Phi, mu0, Sigma0)

Arguments

Details

The p-dimensional multivariate Ornstein-Uhlenbeck (mOU) process Y_t = (Y_{1t}, \ldots, Y_{dt}) satisfies the SDE

dY_t = (\Gamma Y_t + \Lambda)dt + \Phi^{1/2} dB_t,

where B_t = (B_{1t}, \ldots, B_{pt}) is p-dimensional Brownian motion. Its Euler discretization is of the form

Y_{n+1} = Y_n + (\Gamma Y_n + \Lambda) \Delta_n + \Phi^{1/2} \Delta B_n,

where Y_n = Y(t_n), \Delta_n = t_{n+1} - t_n and

\Delta B_n = B(t_{n+1}) - B(t_n) \stackrel{\textnormal{ind}}{\sim} \mathcal N(0, \Delta_n).

Thus, Y_0, \ldots, Y_N is multivariate normal Markov chain for which the marginal distribution of any subset of timepoints and/or components can be efficiently calculated using the Kalman filter. This can be used to check the MCMC output of sde.post() as in the example.

Value

Scalar value of the loglikelihood. See Details.

Examples

# bivariate OU model
bmod <- sde.examples("biou")

# simulate some data

# true parameter values
Gamma0 <- .1 * crossprod(matrix(rnorm(4),2,2))
Lambda0 <- rnorm(2)
Phi0 <- crossprod(matrix(rnorm(4),2,2))
Psi0 <- chol(Phi0) # precompiled model uses the Cholesky scale
theta0 <- c(Gamma0, Lambda0, Psi0[c(1,3,4)])
names(theta0) <- bmod$param.names
# initial value
Y0 <- rnorm(2)
names(Y0) <- bmod$data.names

# simulation
dT <- runif(1, max = .1) # time step
nObs <- 10
bsim <- sde.sim(bmod, x0 = Y0, theta = theta0,
                dt = dT, dt.sim = dT, nobs = nObs)
YObs <- bsim$data

# inference via MCMC
binit <- sde.init(bmod, x = YObs, dt = dT, theta = theta0,
                  nvar.obs = 1) # second component is unobserved
# only Lambda1 is unknown
fixed.params <- rep(TRUE, bmod$nparams)
names(fixed.params) <- bmod$param.names
fixed.params["Lambda1"] <- FALSE
# prior on (Lambda1, Y_0)
hyper <- list(mu = c(0,0), Sigma = diag(2))
names(hyper$mu) <- bmod$data.names
dimnames(hyper$Sigma) <- rep(list(bmod$data.names), 2)

# posterior sampling
nsamples <- 1e5
burn <- 1e3
bpost <- sde.post(bmod, binit, hyper = hyper,
                  fixed.params = fixed.params,
                  nsamples = nsamples, burn = burn)
L1.mcmc <- bpost$params[,"Lambda1"]

# analytic posterior
L1.seq <- seq(min(L1.mcmc), max(L1.mcmc), len = 500)
L1.loglik <- sapply(L1.seq, function(l1) {
  lambda <- Lambda0
  lambda[1] <- l1
  mou.loglik(X = YObs, dt = dT, nvar.obs = 1,
             Gamma = Gamma0, Lambda = lambda, Phi = Phi0,
             mu0 = hyper$mu, Sigma0 = hyper$Sigma)
})
# normalize density
L1.Kalman <- exp(L1.loglik - max(L1.loglik))
L1.Kalman <- L1.Kalman/sum(L1.Kalman)/(L1.seq[2]-L1.seq[1])

# compare
hist(L1.mcmc, breaks = 100, freq = FALSE,
     main = expression(p(Lambda[1]*" | "*bold(Y)[1])),
     xlab = expression(Lambda[1]))
lines(L1.seq, L1.Kalman, col = "red")
legend("topright", legend = c("Analytic", "MCMC"),
       pch = c(NA, 22), lty = c(1, NA), col = c("red", "black"))

Argument checking for the default multivariate normal prior.

Description

Argument checking for the default multivariate normal prior.

Usage

mvn.hyper.check(hyper, param.names, data.names)

Arguments

Details

This function is not meant to be called directly by the user, but rather to parse the hyper-parameters of a default multivariate normal prior distribution to be passed to the C++ code in sde.prior() and sde.post(). This default prior is multivariate normal on the elements of ⁠(theta, x0)⁠ specified by each of names(mu), rownames(Sigma), and colnames(Sigma). The remaining components are given Lebesgue priors, or a full Lebesgue prior if hyper == NULL. If the names of mu and Sigma are inconsistent an error is thrown.

Value

A list with the following elements:

mean

The mean vector.

cholSd

The upper upper Cholesky factor of the variance matrix.

thetaId

The index of the corresponding variables in theta.

xId

The index of the corresponding variables in x0.

SDE diffusion function.

Description

Computes the SDE model's diffusion function given data and parameter values.

Usage

sde.diff(model, x, theta)

Arguments

Value

A matrix with ndims^2 columns containing the diffusion function evaluated at x and theta. Each row corresponds to the upper triangular Cholesky factor of the diffusion matrix. If either input contains invalid SDE data or parameters an error is thrown.

Examples

# load Heston's model
hmod <- sde.examples("hest")
#'
# single input
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
x0 <- c(X = log(1000), Z = 0.1)
sde.diff(model = hmod, x = x0, theta = theta)
#'
# multiple inputs
nreps <- 10
Theta <- apply(t(replicate(nreps, theta)), 2, jitter)
X0 <- apply(t(replicate(nreps, x0)), 2, jitter)
sde.diff(model = hmod, x = X0, theta = Theta)
#'
# mixed inputs
sde.diff(model = hmod, x = x0, theta = Theta)

SDE drift function.

Description

Computes the SDE model's drift function given data and parameter values.

Usage

sde.drift(model, x, theta)

Arguments

Value

A matrix with ndims columns containing the drift function evaluated at x and theta. If either input contains invalid SDE data or parameters an error is thrown.

Examples

# load Heston's model
hmod <- sde.examples("hest")

# single input
x0 <- c(X = log(1000), Z = 0.1)
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
sde.drift(model = hmod, x = x0, theta = theta)

# multiple inputs
nreps <- 10
Theta <- apply(t(replicate(nreps,theta)),2,jitter)
X0 <- apply(t(replicate(nreps,x0)),2,jitter)
sde.drift(model = hmod, x = X0, theta = Theta)

Example SDE models.

Description

Provides sample ⁠C++⁠ code for several SDE models.

Usage

sde.examples(
  model = c("hest", "pgnet", "lotvol", "biou", "eou"),
  file.only = FALSE
)

Arguments

Details

All pre-compiled models are with the default prior and with OpenMP disabled. A full description of the example models can be found in the package vignette; to view it run vignette("msde-exmodels").

Value

An sde.model object, or the path to the C++ model header file.

See Also

sde.make.model() for sde.model objects, mvn.hyper.check() for specification of the default prior.

Examples

# Heston's model
hmod <- sde.examples("hest") # load pre-compiled model

# inspect model's C++ code
hfile <- sde.examples("hest", file.only = TRUE)
cat(readLines(hfile), sep = "\n")

## Not run: 
# compile it from scratch
param.names <- c("alpha", "gamma", "beta", "sigma", "rho")
data.names <- c("X", "Z")
hmod <- sde.make.model(ModelFile = hfile,
                       param.names = param.names,
                       data.names = data.names)

## End(Not run)

MCMC initialization.

Description

Specifies the observed SDE data, interobservation times, initial parameter and missing data values to be supplied to sde.post().

Usage

sde.init(model, x, dt, m = 1, nvar.obs, theta)

Arguments

Value

An sde.init object, corresponding to a list with elements:

data

An ⁠ncomp x ndims⁠ matrix of complete data, where ncomp = N_m = m * (nobs-1)+1.

dt.m

The complete data interobservation time, dt_m = dt/m.

nvar.obs.m

The number of variables observed per row of data. Note that nvar.obs.m[(i-1)*m+1] == nvar.obs[ii], and that nvar.obs.m[i-1] == 0 if i is not a multiple of m.

params

Parameter initial values.

Examples

# load Heston's model
hmod <- sde.examples("hest")

# generate some observed data
nObs <- 5
x0 <- c(X = log(1000), Z = 0.1)
X0 <- apply(t(replicate(nObs, x0)), 2, jitter)
dT <- .6
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)

# no missing data
sde.init(model = hmod, x = X0, dt = dT, theta = theta)

# all but endpoint volatilities are missing
sde.init(model = hmod, x = X0, dt = dT, m = 1,
         nvar.obs = c(2, rep(1, nObs-2), 2), theta = theta)

# all volatilities missing,
# two completely missing SDE timepoints between observations
m <- 3 # divide each observation interval into m equally spaced timepoints
sde.init(model = hmod, x = X0, dt = dT,
         m = m, nvar.obs = 1, theta = theta)

SDE loglikelihood function.

Description

Evaluates the loglikelihood function given SDE data and parameter values.

Usage

sde.loglik(model, x, dt, theta, ncores = 1)

Arguments

Value

A vector of loglikelihood evaluations, of the same length as the third dimension of x and/or first dimension of theta. If input contains invalid data or parameters an error is thrown.

Examples

# load Heston's model
hmod <- sde.examples("hest")

# Simulate data
nreps <- 10
nobs <- 100
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
Theta <- apply(t(replicate(nreps, theta)), 2, jitter)
x0 <- c(X = log(1000), Z = 0.1)
X0 <- apply(t(replicate(nreps,x0)), 2, jitter)
dT <- 1/252
hsim <- sde.sim(model = hmod, x0 = X0, theta = Theta,
                dt = dT, dt.sim = dT/10, nobs = nobs, nreps = nreps)

# single parameter, single data
sde.loglik(model = hmod, x = hsim$data[,,1], dt = dT, theta = theta)
# multiple parameters, single data
sde.loglik(model = hmod, x = hsim$data[,,1], dt = dT, theta = Theta)
# multiple parameters, multiple data
sde.loglik(model = hmod, x = hsim$data, dt = dT, theta = Theta)

Create an SDE model object.

Description

Compiles the C++ code for various SDE-related algorithms and makes the routines available within R.

Usage

sde.make.model(
  ModelFile,
  PriorFile = "default",
  data.names,
  param.names,
  hyper.check,
  OpenMP = FALSE,
  ...
)

Arguments

Value

An sde.model object, consisting of a list with the following elements:

ptr

Pointer to C++ sde object (sdeRobj) implementing the member functions: drift/diffusion, data/parameter validators, loglikelihood, prior distribution, forward simulation, MCMC algorithm for Bayesian inference.

ndims, nparams

The number of SDE components and parameters.

data.names, param.names

The names of the SDE components and parameters.

omp

A logical flag for whether or not the model was compiled for multicore functionality with OpenMP.

See Also

sde.drift(), sde.diff(), sde.valid(), sde.loglik(), sde.prior(), sde.sim(), sde.post().

Examples

# header (C++) file for Heston's model
hfile <- sde.examples("hest", file.only = TRUE)
cat(readLines(hfile), sep = "\n")


# compile the model
param.names <- c("alpha", "gamma", "beta", "sigma", "rho")
data.names <- c("X", "Z")
hmod <- sde.make.model(ModelFile = hfile,
                       param.names = param.names,
                       data.names = data.names)

hmod

MCMC sampler for the SDE posterior.

Description

A Metropolis-within-Gibbs sampler for the Euler-Maruyama approximation to the true posterior density.

Usage

sde.post(
  model,
  init,
  hyper,
  nsamples,
  burn,
  mwg.sd = NULL,
  adapt = TRUE,
  loglik.out = FALSE,
  last.miss.out = FALSE,
  update.data = TRUE,
  data.out,
  update.params = TRUE,
  fixed.params,
  ncores = 1,
  verbose = TRUE
)

Arguments

Details

The Metropolis-within-Gibbs (MWG) jump sizes can be specified as a scalar, a vector or length nparams + ndims, or a named vector containing the elements defined by sde.init$nvar.obs.m[1] (the missing variables in the first SDE observation) and fixed.params (the SDE parameters which are not held fixed). The default jump sizes for each MWG random variable are ⁠.25 * |initial_value|⁠ when ⁠|initial_value| > 0⁠, and 1 otherwise.

adapt == TRUE implements an adaptive MCMC proposal by Roberts and Rosenthal (2009). At step n of the MCMC, the jump size of each MWG random variable is increased or decreased by \delta(n), depending on whether the cumulative acceptance rate is above or below the optimal value of 0.44. If \sigma_n is the size of the jump at step n, then the next jump size is determined by

\log(\sigma_{n+1}) = \log(\sigma_n) \pm \delta(n), \qquad \delta(n) = \min(.01, 1/n^{1/2}).

When adapt is not logical, it is a list with elements max and rate, such that delta(n) = min(max, 1/n^rate). These elements can be scalars or vectors in the same manner as mwg.sd.

For SDE models with thousands of latent variables, data.out can be used to thin the MCMC missing data output. An integer vector or scalar returns specific or evenly-spaced posterior samples from the ⁠ncomp x ndims⁠ complete data matrix. A list with elements isamples, icomp, and idims determines which samples, time points, and SDE variables to return. The first of these can be a scalar or vector with the same meaning as before.

Value

A list with elements:

params

An ⁠nsamples x nparams⁠ matrix of posterior parameter draws.

data

A 3-d array of posterior missing data draws, for which the output dimensions are specified by data.out.

init

The sde.init object which initialized the sampler.

data.out

A list of three integer vectors specifying which timepoints, variables, and MCMC iterations correspond to the values in the data output.

mwg.sd

A named vector of Metropolis-within-Gibbs standard devations used at the last posterior iteration.

hyper

The hyperparameter specification.

loglik

If loglik.out == TRUE, the vector of nsamples complete data loglikelihoods calculated at each posterior sample.

last.iter

A list with elements data and params giving the last MCMC sample. Useful for resuming the MCMC from that point.

last.miss

If last.miss.out == TRUE, an ⁠nsamples x nmissN⁠ matrix of all posterior draws for the missing data in the final observation. Useful for SDE forecasting at future timepoints.

accept

A named list of acceptance rates for the various components of the MCMC sampler.

References

Roberts, G.O. and Rosenthal, J.S. "Examples of adaptive MCMC." Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. http://www.probability.ca/jeff/ftpdir/adaptex.pdf.

Examples

# Posterior inference for Heston's model
hmod <- sde.examples("hest") # load pre-compiled model

# Simulate data
X0 <- c(X = log(1000), Z = 0.1)
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
dT <- 1/252
nobs <- 1000
hest.sim <- sde.sim(model = hmod, x0 = X0, theta = theta,
                    dt = dT, dt.sim = dT/10, nobs = nobs)

# initialize MCMC sampler
# both components observed, no missing data between observations
init <- sde.init(model = hmod, x = hest.sim$data,
                 dt = hest.sim$dt, theta = theta)

# Initialize posterior sampling argument
nsamples <- 1e4
burn <- 1e3
hyper <- NULL # flat prior
hest.post <- sde.post(model = hmod, init = init, hyper = hyper,
                      nsamples = nsamples, burn = burn)

# plot the histogram for the sampled parameters
par(mfrow = c(2,3))
for(ii in 1:length(hmod$param.names)) {
  hist(hest.post$params[,ii],breaks=100, freq = FALSE,
       main = parse(text = hmod$param.names[ii]), xlab = "")
}

SDE prior function.

Description

Evaluates the SDE prior given data, parameter, and hyperparameter values.

Usage

sde.prior(model, theta, x, hyper)

Arguments

Details

The prior is constructed at the ⁠C++⁠ level by defining a function (i.e., public member)

double logPrior(double *theta, double *x)

within the sdePrior class. At the R level, the hyper.check argument of sde.make.model() is a function with arguments hyper, param.names, data.names used to convert hyper into a list of NULL or double-vectors which get passed on to the ⁠C++⁠ code. This function can also be used to throw R-level errors to protect the ⁠C++⁠ code from invalid inputs, as is done for the default prior in mvn.hyper.check(). For a full example see the "Custom Prior" section in vignette("msde-quicktut").

Value

A vector of log-prior densities evaluated at the inputs.

Examples

hmod <- sde.examples("hest") # load Heston's model

# setting prior for 3 parameters
rv.names <- c("alpha","gamma","rho")
mu <- rnorm(3)
Sigma <- crossprod(matrix(rnorm(9),3,3))
names(mu) <- rv.names
colnames(Sigma) <- rv.names
rownames(Sigma) <- rv.names
hyper <- list(mu = mu, Sigma = Sigma)

# Simulate data
nreps <- 10
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
x0 <- c(X = log(1000), Z = 0.1)
Theta <- apply(t(replicate(nreps,theta)),2,jitter)
X0 <- apply(t(replicate(nreps,x0)),2,jitter)

sde.prior(model = hmod, x = X0, theta = Theta, hyper = hyper)

Simulation of multivariate SDE trajectories.

Description

Simulates a discretized Euler-Maruyama approximation to the true SDE trajectory.

Usage

sde.sim(
  model,
  x0,
  theta,
  dt,
  dt.sim,
  nobs,
  burn = 0,
  nreps = 1,
  max.bad.draws = 5000,
  verbose = TRUE
)

Arguments

Details

The simulation algorithm is a Markov process with Y_0 = x_0 and

Y_{t+1} \sim \mathcal{N}(Y_t + \mathrm{dr}(Y_t, \theta) dt_{\mathrm{sim}}, \mathrm{df}(Y_t, \theta) dt_{\mathrm{sim}}),

where \mathrm{dr}(y, \theta) is the SDE drift function and \mathrm{df}(y, \theta) is the diffusion function on the variance scale. At each step, a while-loop is used until a valid SDE draw is produced. The simulation algorithm terminates after nreps trajectories are drawn or once a total of max.bad.draws are reached.

Value

A list with elements:

data

An array of size ⁠nobs x ndims x nreps⁠ containing the simulated SDE trajectories.

params

The vector or matrix of parameter values used to generate the data.

⁠dt, dt.sim⁠

The actual and internal interobservation times.

nbad

The total number of bad draws.

Examples

# load pre-compiled model
hmod <- sde.examples("hest")

# initial values
x0 <- c(X = log(1000), Z = 0.1)
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)

# simulate data
dT <- 1/252
nobs <- 2000
burn <- 500
hsim <- sde.sim(model = hmod, x0 = x0, theta = theta,
                dt = dT, dt.sim = dT/10,
                nobs = nobs, burn = burn)

par(mfrow = c(1,2))
plot(hsim$data[,"X"], type = "l", xlab = "Time", ylab = "",
     main = expression(X[t]))
plot(hsim$data[,"Z"], type = "l", xlab = "Time", ylab = "",
     main = expression(Z[t]))

SDE data and parameter validators.

Description

Checks whether input SDE data and parameters are valid.

Usage

sde.valid.data(model, x, theta)

sde.valid.params(model, theta)

Arguments

Value

A logical scalar or vector indicating whether the given data/parameter pair is valid.

Examples

# Heston's model
# valid data is: Z > 0
# valid parameters are: gamma, sigma > 0, |rho| < 1, beta > .5 * sigma^2
hmod <- sde.examples("hest") # load model

theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)

# valid data
x0 <- c(X = log(1000), Z = 0.1)
sde.valid.data(model = hmod, x = x0, theta = theta)

# invalid data
x0 <- c(X = log(1000), Z = -0.1)
sde.valid.data(model = hmod, x = x0, theta = theta)

# valid parameters
theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)
sde.valid.params(model = hmod, theta = theta)

# invalid parameters
theta <- c(alpha = 0.1, gamma = -4, beta = 0.8, sigma = 0.6, rho = -0.8)
sde.valid.params(model = hmod, theta = theta)