Movatterモバイル変換

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 andvignette("msde-exmodels") for several example models.

Author(s)

Maintainer: Martin Lysymlysy@uwaterloo.ca

Authors:

• Feiyu Zhu

• JunYong Tong

Other contributors:

• Trevor Kitt [contributor]

• Nigel Delaney [contributor]

Examples

# Posterior inference for Heston's model# compile modelhfile <- 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 versionhmod <- sde.examples("hest")# Simulate dataX0 <- 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/252nobs <- 1000hest.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 observationsinit <- sde.init(model = hmod, x = hest.sim$data,                 dt = hest.sim$dt, theta = theta)# Initialize posterior sampling argumentnsamples <- 1e4burn <- 1e3hyper <- NULL # flat priorhest.post <- sde.post(model = hmod, init = init, hyper = hyper,                      nsamples = nsamples, burn = burn)# plot the histogram for the sampled parameterspar(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

Thep-dimensional multivariate Ornstein-Uhlenbeck (mOU) processY_t = (Y_{1t}, \ldots, Y_{dt}) satisfies the SDE

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

whereB_t = (B_{1t}, \ldots, B_{pt}) isp-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,

whereY_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 ofsde.post() as in the example.

Value

Scalar value of the loglikelihood. See Details.

Examples

# bivariate OU modelbmod <- sde.examples("biou")# simulate some data# true parameter valuesGamma0 <- .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 scaletheta0 <- c(Gamma0, Lambda0, Psi0[c(1,3,4)])names(theta0) <- bmod$param.names# initial valueY0 <- rnorm(2)names(Y0) <- bmod$data.names# simulationdT <- runif(1, max = .1) # time stepnObs <- 10bsim <- sde.sim(bmod, x0 = Y0, theta = theta0,                dt = dT, dt.sim = dT, nobs = nObs)YObs <- bsim$data# inference via MCMCbinit <- sde.init(bmod, x = YObs, dt = dT, theta = theta0,                  nvar.obs = 1) # second component is unobserved# only Lambda1 is unknownfixed.params <- rep(TRUE, bmod$nparams)names(fixed.params) <- bmod$param.namesfixed.params["Lambda1"] <- FALSE# prior on (Lambda1, Y_0)hyper <- list(mu = c(0,0), Sigma = diag(2))names(hyper$mu) <- bmod$data.namesdimnames(hyper$Sigma) <- rep(list(bmod$data.names), 2)# posterior samplingnsamples <- 1e5burn <- 1e3bpost <- sde.post(bmod, binit, hyper = hyper,                  fixed.params = fixed.params,                  nsamples = nsamples, burn = burn)L1.mcmc <- bpost$params[,"Lambda1"]# analytic posteriorL1.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 densityL1.Kalman <- exp(L1.loglik - max(L1.loglik))L1.Kalman <- L1.Kalman/sum(L1.Kalman)/(L1.seq[2]-L1.seq[1])# comparehist(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 insde.prior() andsde.post(). This default prior is multivariate normal on the elements of⁠(theta, x0)⁠ specified by each ofnames(mu),rownames(Sigma), andcolnames(Sigma). The remaining components are given Lebesgue priors, or a full Lebesgue prior ifhyper == NULL. If the names ofmu andSigma 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 intheta.

xId

The index of the corresponding variables inx0.

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 withndims^2 columns containing the diffusion function evaluated atx andtheta. 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 modelhmod <- sde.examples("hest")#'# single inputtheta <- 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 inputsnreps <- 10Theta <- apply(t(replicate(nreps, theta)), 2, jitter)X0 <- apply(t(replicate(nreps, x0)), 2, jitter)sde.diff(model = hmod, x = X0, theta = Theta)#'# mixed inputssde.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 withndims columns containing the drift function evaluated atx andtheta. If either input contains invalid SDE data or parameters an error is thrown.

Examples

# load Heston's modelhmod <- sde.examples("hest")# single inputx0 <- 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 inputsnreps <- 10Theta <- 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 withOpenMP disabled. A full description of the example models can be found in the package vignette; to view it runvignette("msde-exmodels").

Value

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

See Also

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

Examples

# Heston's modelhmod <- sde.examples("hest") # load pre-compiled model# inspect model's C++ codehfile <- sde.examples("hest", file.only = TRUE)cat(readLines(hfile), sep = "\n")## Not run: # compile it from scratchparam.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 tosde.post().

Usage

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

Arguments

Value

Ansde.init object, corresponding to a list with elements:

data

An⁠ncomp x ndims⁠ matrix of complete data, wherencomp = 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 ofdata. Note thatnvar.obs.m[(i-1)*m+1] == nvar.obs[ii], and thatnvar.obs.m[i-1] == 0 ifi is not a multiple ofm.

params

Parameter initial values.

Examples

# load Heston's modelhmod <- sde.examples("hest")# generate some observed datanObs <- 5x0 <- c(X = log(1000), Z = 0.1)X0 <- apply(t(replicate(nObs, x0)), 2, jitter)dT <- .6theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)# no missing datasde.init(model = hmod, x = X0, dt = dT, theta = theta)# all but endpoint volatilities are missingsde.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 observationsm <- 3 # divide each observation interval into m equally spaced timepointssde.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 ofx and/or first dimension oftheta. If input contains invalid data or parameters an error is thrown.

Examples

# load Heston's modelhmod <- sde.examples("hest")# Simulate datanreps <- 10nobs <- 100theta <- 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/252hsim <- sde.sim(model = hmod, x0 = X0, theta = Theta,                dt = dT, dt.sim = dT/10, nobs = nobs, nreps = nreps)# single parameter, single datasde.loglik(model = hmod, x = hsim$data[,,1], dt = dT, theta = theta)# multiple parameters, single datasde.loglik(model = hmod, x = hsim$data[,,1], dt = dT, theta = Theta)# multiple parameters, multiple datasde.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

Ansde.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 withOpenMP.

See Also

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

Examples

# header (C++) file for Heston's modelhfile <- sde.examples("hest", file.only = TRUE)cat(readLines(hfile), sep = "\n")# compile the modelparam.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 lengthnparams + ndims, or a named vector containing the elements defined bysde.init$nvar.obs.m[1] (the missing variables in the first SDE observation) andfixed.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 stepn 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 stepn, 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}).

Whenadapt is not logical, it is a list with elementsmax andrate, such thatdelta(n) = min(max, 1/n^rate). These elements can be scalars or vectors in the same manner asmwg.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 elementsisamples,icomp, andidims 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 bydata.out.

init

Thesde.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 thedata output.

mwg.sd

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

hyper

The hyperparameter specification.

loglik

Ifloglik.out == TRUE, the vector ofnsamples complete data loglikelihoods calculated at each posterior sample.

last.iter

A list with elementsdata andparams giving the last MCMC sample. Useful for resuming the MCMC from that point.

last.miss

Iflast.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 modelhmod <- sde.examples("hest") # load pre-compiled model# Simulate dataX0 <- 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/252nobs <- 1000hest.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 observationsinit <- sde.init(model = hmod, x = hest.sim$data,                 dt = hest.sim$dt, theta = theta)# Initialize posterior sampling argumentnsamples <- 1e4burn <- 1e3hyper <- NULL # flat priorhest.post <- sde.post(model = hmod, init = init, hyper = hyper,                      nsamples = nsamples, burn = burn)# plot the histogram for the sampled parameterspar(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 thesdePrior class. At theR level, thehyper.check argument ofsde.make.model() is a function with argumentshyper,param.names,data.names used to converthyper into a list ofNULL or double-vectors which get passed on to the⁠C++⁠ code. This function can also be used to throwR-level errors to protect the⁠C++⁠ code from invalid inputs, as is done for the default prior inmvn.hyper.check(). For a full example see the "Custom Prior" section invignette("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 parametersrv.names <- c("alpha","gamma","rho")mu <- rnorm(3)Sigma <- crossprod(matrix(rnorm(9),3,3))names(mu) <- rv.namescolnames(Sigma) <- rv.namesrownames(Sigma) <- rv.nameshyper <- list(mu = mu, Sigma = Sigma)# Simulate datanreps <- 10theta <- 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 withY_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 thevariance scale. At each step, a while-loop is used until a valid SDE draw is produced. The simulation algorithm terminates afternreps trajectories are drawn or once a total ofmax.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 modelhmod <- sde.examples("hest")# initial valuesx0 <- c(X = log(1000), Z = 0.1)theta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)# simulate datadT <- 1/252nobs <- 2000burn <- 500hsim <- 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^2hmod <- sde.examples("hest") # load modeltheta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)# valid datax0 <- c(X = log(1000), Z = 0.1)sde.valid.data(model = hmod, x = x0, theta = theta)# invalid datax0 <- c(X = log(1000), Z = -0.1)sde.valid.data(model = hmod, x = x0, theta = theta)# valid parameterstheta <- c(alpha = 0.1, gamma = 1, beta = 0.8, sigma = 0.6, rho = -0.8)sde.valid.params(model = hmod, theta = theta)# invalid parameterstheta <- c(alpha = 0.1, gamma = -4, beta = 0.8, sigma = 0.6, rho = -0.8)sde.valid.params(model = hmod, theta = theta)