JGQD.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 jump diffusion model. Given a set of starting parameters, a MCMC chain is returned for further analysis.
The structure of the model is predefined and coefficients may be provided for models nested within the generalized quadratic diffusion framework.JGQD.mcmc() performs inference using the Metropolis-Hastings algorithm for jump diffusions of the form:
ScalarEqn1.png
where
ScalarEqn2.png
ScalarEqn3.png
and
ScalarEqn4.png
describes a Poisson process with jumps of the form:
ScalarEqn6.png
arriving with intensity
ScalarEqn5.png
subject to a jump distribition of the form:
ScalarEqn7.png
JGQD.mcmc(X, time, mesh = 10, theta, sds, updates = 1000, burns = min(round(updates/2), 25000), Jtype = "Add", Jdist = "Normal", Dtype = "Saddlepoint", RK.order = 4, exclude = NULL, plot.chain = TRUE, wrt = FALSE, Tag = NA, factorize = TRUE, print.output = TRUE, decode = TRUE, palette = 'mono')
time
).
X
.
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
[i]^2)
burns
values are omitted from the inference, although the entire chain is returned.
NULL
.TRUE
(default), a trace plot is made of the resulting MCMC chain (see details).
'Saddlepoint'
is supported in the current version of the software.Tag
can be used to name (tag) an MCMC run e.g. Tag='Run_1'
TRUE
a .cpp file will be written to the current directory. For bug report diagnostics.TRUE
, model information is printed to the console.TRUE
. palette = 'mono'
, otherwise a qualitative palette will be used.theta
. For example: G0 <- 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: G0 <- 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: G0 <- function(t){0.1*(10+0.2*pow(sin(2*pi*t),3))}
. JGQD.mcmc()
operates by searching the workspace for functions with names that match the coefficients of the predefined stochastic differential equation. Only the required coefficients need to be specified e.g. G0(t)
,G1(t)
and Q0(t)
for an Ornstein-Uhlenbeck model. Unspecified coefficients are ignored. When a new model is to be defined, the current model may be removed from the workspace by using the JGQD.remove
function, after which the new coefficients may be supplied.
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.
JGQD.remove
, BiJGQD.mcmc
.
#===============================================================================
# For detailed notes and examples on how to use the JGQD.mcmc() function, see
# the following vignette:
RShowDoc('Part_4_Likelihood_Inference',type='html','DiffusionRjgqd')
#===============================================================================
Run the code above in your browser using DataLab