DiffusionRjgqd (version 0.1.1)

BiJGQD.mcmc: MCMC Inference on Bivariate Jump Generalized Quadratic Diffusions (2D JGQDs).

Description

BiJGQD.mcmc() uses parametrised coefficients (provided by the user as R-functions) to construct a C++ program in real time that allows the user to perform Bayesian inference on the resulting diffusion model. Given a set of starting parameters and other input parameters, a MCMC chain is returned for further analysis. The user may specify any model within the JGQD framework by defining parametrised functions giving the form of the coefficients of the model. That is, BiJGQD.density generates approximate transitional densities for a class of bivariate jump diffusion processes with SDE:

BivJumpDiff1.png

where

BivJumpDiff4.png

BivJumpDiff3.png

and

BivJumpDiff2.png

describes a bivariate Poisson process with jump matrix:

BivJumpDiff5.png

and intensity vector

BivJumpDiff6.png

Usage

BiJGQD.mcmc(X, time, mesh = 10, theta, sds, updates = 10, burns = min(round(updates/2),25000), exclude = NULL, plot.chain = TRUE, RK.order = 4, wrt = FALSE, Tag = NA, Dtype = "Saddlepoint", Jdist = "MVNormal", Jtype ='Add', adapt = 0, print.output = TRUE, decode = TRUE, palette = 'mono')

Arguments

X

A matrix containing rows of data points to be modelled. Although observations are allowed to be non-equidistant, observations in both dimensions are assumed to occur at the same time epochs (i.e. time gives the time signature for both dimensions).

time

A vector containing the time epochs at which observations were made.

mesh

The number of mesh points in the time discretization.

theta

The parameter vector of the process. theta are taken as the starting values of the MCMC chain and gives the dimension of the parameter vector used to calculate the DIC. Care should be taken to ensure that each element in theta is in fact used within the coefficient-functions, otherwise redundant parameters will be counted in the calculation of the DIC.

sds

Proposal distribution standard deviations. That is, for the i-th parameter the proposal distribution is ~ Normal(...,sds[i]^2).

updates

The number of MCMC updates/iterations to perform (including burn-in).

burns

The number of MCMC updates/iterations to burn.

exclude

Vector indicating which transitions to exclude from the analysis. Default = NULL.

plot.chain

If = TRUE (default), a trace plot of the MCMC chain will be made along with a trace of the acceptance rate.

RK.order

The order of the Runge-Kutta solver used to approximate the trajectories of cumulants. Must be 4 (default) or 10.

Tag

Tag can be used to name (tag) an MCMC run e.g. Tag='Run_1'

Dtype

The density approximant to use. Valid types are "Saddlepoint" (default), "Edgeworth" or "Normal".

Jdist

Valid entries are 'MVNormal' (currently).

Jtype

Valid types are 1 or 2.

adapt

For development purposes.

wrt

If TRUE a .cpp file will be written to the current directory. For bug report diagnostics.

print.output

If TRUE information about the model and algorithm is printed to the console.

decode

Should the algorithm estimate jump detection probabilities? Default value is TRUE.

palette

Colour palette for drawing trace plots. Default palette = 'mono', otherwise a qualitative palette will be used.

Value

Syntactical jargon

Synt. [1]: The coefficients of the 2D JGQD may be parameterized using the reserved variable theta. For example: a00 <- function(t){theta[1]*(theta[2]+sin(2*pi*(t-theta[3])))}. Synt. [2]: Due to syntactical differences between R and C++ special functions have to be used when terms that depend on t. When the function cannot be separated in to terms that contain a single t, the prod(a,b) function must be used. For example: a00 <- function(t){0.1*(10+0.2*sin(2*pi*t)+0.3*prod(sqrt(t),1+cos(3*pi*t)))}. Here sqrt(t)*cos(3*pi*t) constitutes the product of two terms that cannot be written i.t.o. a single t. To circumvent this isue, one may use the prod(a,b) function. Synt. [3]: Similarly, the ^ - operator is not overloaded in C++. Instead the pow(x,p) function may be used to calculate x^p. For example sin(2*pi*t)^3 in: a00 <- function(t){0.1*(10+0.2*pow(sin(2*pi*t),3))}.

Warning

Warning [1]: The parameter mesh is used to discretize the transition horizons between successive observations. It is thus important to ensure that mesh is not too small when large time differences are present in time. Check output for max(dt) and divide by mesh.

References

Updates available on GitHub at https://github.com/eta21.

Daniels, H.E. 1954 Saddlepoint approximations in statistics. Ann. Math. Stat., 25:631--650.

Eddelbuettel, D. and Romain, F. 2011 Rcpp: Seamless R and C++ integration. Journal of Statistical Software, 40(8):1--18,. URL http://www.jstatsoft.org/v40/i08/.

Eddelbuettel, D. 2013 Seamless R and C++ Integration with Rcpp. New York: Springer. ISBN 978-1-4614-6867-7.

Eddelbuettel, D. and Sanderson, C. 2014 Rcpparmadillo: Accelerating R with high-performance C++ linear algebra. Computational Statistics and Data Analysis, 71:1054--1063. URL http://dx.doi.org/10.1016/j.csda.2013.02.005.

Feagin, T. 2007 A tenth-order Runge-Kutta method with error estimate. In Proceedings of the IAENG Conf. on Scientifc Computing.

Varughese, M.M. 2013 Parameter estimation for multivariate diffusion systems. Comput. Stat. Data An., 57:417--428.

See Also

JGQD.remove and JGQD.mcmc.

Examples

#===============================================================================

 # For detailed notes and examples on how to use the BiJGQD.mcmc() function, see
 # the following vignette:

 RShowDoc('Part_4_Likelihood_Inference',type='html','DiffusionRjgqd')

#===============================================================================