Compiles the C++ code for various SDE-related algorithms and makes the routines available within R.
sde.make.model( ModelFile, PriorFile = "default", data.names, param.names, hyper.check, OpenMP = FALSE, ... )
| ModelFile | Path to the header file where the SDE model is defined. |
|---|---|
| PriorFile | Path to the header file where the SDE prior is defined. See |
| data.names | Vector of names for the SDE components. Defaults to |
| param.names | Vector of names for the SDE parameters. Defaults to |
| hyper.check | A function with arguments |
| OpenMP | Logical; whether the model is compiled with |
| ... | additional arguments to |
An sde.model object, consisting of a list with the following elements:
ptrPointer 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, nparamsThe number of SDE components and parameters.
data.names, param.namesThe names of the SDE components and parameters.
ompA logical flag for whether or not the model was compiled for multicore functionality with OpenMP.
# header (C++) file for Heston's model hfile <- sde.examples("hest", file.only = TRUE) cat(readLines(hfile), sep = "\n")#> /// @file hestModel.h #> #> #ifndef hestModel_h #> #define hestModel_h 1 #> #> /// SDE model class for Heston's stochastic volatility model. #> /// #> /// The model is given by #> /// ``` #> /// dXt = (alpha - .125 * Zt^2)dt + .5 * Zt dB_Xt #> /// dZt = (beta/Zt - .5*gamma * Zt)dt + sigma * dB_Zt #> /// cor(B_Xt, B_Zt) = rho #> /// ``` #> /// #> /// The data vector is `x = (X, Z)` and the parameter vector is `theta = (alpha, gamma, beta, sigma, rho)`. #> class sdeModel { #> public: #> static const int nParams = 5; ///< Number of model parameters. #> static const int nDims = 2; ///< Number of SDE dimensions. #> static const bool sdDiff = true; ///< Diffusion is on the standard deviation scale. #> static const bool diagDiff = false; ///< Diffusion is not diagonal. #> /// SDE drift function. #> void sdeDr(double *dr, double *x, double *theta); #> /// SDE diffusion function. #> void sdeDf(double *df, double *x, double *theta); #> /// SDE data validator. #> bool isValidData(double *x, double *theta); #> /// SDE parameter validator. #> bool isValidParams(double *theta); #> }; #> #> /// @param[out] dr Array into which to store the calculated drift. #> /// @param[in] x Array of SDE components at a given time point. #> /// @param[in] theta Array of SDE parameters. #> inline void sdeModel::sdeDr(double *dr, double *x, double *theta) { #> dr[0] = (theta[0] - .125 * x[1]*x[1]); // x #> dr[1] = (theta[2]/x[1] - .5 * theta[1]*x[1]); // z #> return; #> } #> #> /// @param[out] df Array into which to store the calculated diffusion matrix. #> /// @param[in] x Array of SDE components at a given time point. #> /// @param[in] theta Array of SDE parameters. #> inline void sdeModel::sdeDf(double *df, double *x, double *theta) { #> df[0] = .5 * x[1]; #> df[2] = theta[3]; #> df[3] = sqrt(1.0-theta[4]*theta[4]) * df[2]; #> df[2] *= theta[4]; #> return; #> } #> #> /// @param[in] x Array of SDE components at a given time point. #> /// @param[in] theta Array of SDE parameters. #> /// #> /// @return Whether or not the SDE data `x` is valid. In this case we must have `Zt > 0`. #> inline bool sdeModel::isValidData(double *x, double *theta) { #> return(x[1] > 0.0); #> } #> #> /// @param[in] theta Array of SDE parameters. #> /// #> /// @return Whether or not the SDE parameters `theta` are valid. In this case we must have `gamma, sigma > 0`, `beta > sigma^2/2`, and `|rho| < 1`. #> inline bool sdeModel::isValidParams(double *theta) { #> bool isValid; #> isValid = (theta[1] > 0) && (theta[3] > 0); #> isValid = isValid && (-1.0 < theta[4]) && (1.0 > theta[4]); #> isValid = isValid && (theta[2] > 0.5 * theta[3] * theta[3]); #> return(isValid); #> } #> #> #endif# \donttest{ # 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#> $cobj #> C++ object <0x7ff96cd8d390> of class 'msde_10cd64c5b27d' <0x7ff972cb7300> #> #> $ndims #> [1] 2 #> #> $nparams #> [1] 5 #> #> $data.names #> [1] "X" "Z" #> #> $param.names #> [1] "alpha" "gamma" "beta" "sigma" "rho" #> #> $hyper.check #> function (hyper, param.names, data.names) #> { #> nparams <- length(param.names) #> ndims <- length(data.names) #> var.names <- c(param.names, data.names) #> prior.names <- c("mu", "Sigma") #> if (!is.null(hyper)) { #> if (is.null(names(hyper)) || !identical(sort(names(hyper)), #> sort(prior.names))) { #> stop("hyper must be NULL or a list with elements mu and Sigma.") #> } #> mu.names <- names(hyper$mu) #> Sigma.rnames <- rownames(hyper$Sigma) #> Sigma.cnames <- colnames(hyper$Sigma) #> if (!identical(mu.names, Sigma.rnames) || !identical(mu.names, #> Sigma.cnames)) { #> stop("names(mu), rownames(Sigma), and colnames(Sigma) are not consistent.") #> } #> if (!is.valid.vars(mu.names, c(param.names, data.names))) { #> stop("names(mu) must be a unique subset of param.names and data.names.") #> } #> var.id <- sapply(mu.names, function(x) which(x == var.names)) #> var.ord <- order(var.id) #> mu <- hyper$mu[var.ord] #> Sigma <- hyper$Sigma[var.ord, var.ord] #> var.id <- sort(var.id) #> theta.id <- var.id[var.id <= nparams] #> x.id <- var.id[var.id > nparams] - nparams #> prior.args <- list(mean = mu, cholSd = chol(Sigma), thetaId = theta.id, #> xId = x.id) #> prior.args <- lapply(prior.args, function(x) { #> if (length(x) == 0) #> NULL #> else as.double(x) #> }) #> } #> else { #> prior.args <- list(NULL) #> } #> prior.args #> } #> <bytecode: 0x7ff973976dd0> #> <environment: namespace:msde> #> #> $omp #> [1] FALSE #> #> attr(,"class") #> [1] "sde.model"# }