R: MLEs for Generalized Quadratic Diffusion Models (GQDs).
GQD.mle
R Documentation
MLEs for Generalized Quadratic Diffusion Models (GQDs).
Description
GQD.mle() constructs a C++ program in real time that allows the user to perform maximum likelihood inference on scalar GQDs. Given a set of starting parameters, the maximum likelihood estimates are calculated via minimization of minus the likelihood function via the built-in R-function optim.
GQD.mcmc() performs inference using the Metropolis-Hastings algorithm for jump diffusions of the form:
Time series (vector) of discretely observed points of the process of interest. These may be non-equidistant observations (see time).
time
A vector of time-stamps associated with each observation in X.
mesh
The number mesh points between any two given data points.
theta
The set of starting parameters for the optimization routine.
control
List of control variables to be used by optim.
method
Method to be used by optim.
exclude
Vector indicating which transitions to exclude from the analysis. Default = NULL.
RK.order
The order of the Runge-Kutta scheme used. Value must be 4 or (default) 10.
Dtype
Character string indicating the type of density approximation (see details) to use. Types: 'Saddlepoint', 'Normal', 'Gamma', 'InvGamma' and 'Beta' are supported (default = 'Saddlepoint').
Trunc
Vector of length 2 containing the cumulant truncation order and the density truncation order respectively. May take on values 4, 6 and 8 with the constraint that Trunc[1] >= Trunc[2]. Default is c(4,4).
P
Normalization parameter indicating the number of points to use when normalizing members of the Pearson system (see details)
alpha
Normalization parameter controlig the mesh concentration when normalizing members of the Pearson system (see details). Increasing alpha decreases concentration around the mean and vice versa (default alpha = 0).
lower,upper
Lower and upper bounds for the normalization range.
Tag
Tag can be used to name (tag) an MCMC run e.g. Tag='Run_1'
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.
Value
opt
The output from optim.
model.info
A list of variables pertaining to inference calculations.
model.info$elapsed.time
The run-time, in h/m/s format,of the MCMC procedure (excluding compile time).
model.info$time.homogeneous
‘No’ if the model has time-homogeneous coefficients and ‘Yes’ otherwise.
model.info$p
The dimension of theta.
Syntactical jargon
Synt. [1]: The coefficients of the GQD may be parameterized using the reserved variable theta. For example:
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:
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:
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.
Warning [2]: Note that minus the likelihood is minimized, as such the optim output (hessian) needs to be adjusted accordingly if used for calculating confidence intervals. Furthermore, GQD.mle may be temperamental under certain conditions
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
GQD.remove, GQD.mcmc, BiGQD.mcmc, BiGQD.mle, GQD.passage and GQD.TIpassage.
Examples
#===============================================================================
# Simulate a time inhomogeneous diffusion.
#-------------------------------------------------------------------------------
data(SDEsim1)
attach(SDEsim1)
par(mfrow=c(1,1))
expr1=expression(dX[t]==2*(5+3*sin(0.5*pi*t)-X[t])*dt+0.5*sqrt(X[t])*dW[t])
plot(Xt~time,type='l',col='blue',xlab='Time (t)',ylab=expression(X[t]),main=expr1)
#------------------------------------------------------------------------------
# Define coefficients of the process.
#------------------------------------------------------------------------------
GQD.remove()
G0 <- function(t){theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))}
G1 <- function(t){-theta[1]}
Q0 <- function(t){theta[4]*theta[4]}
theta.start <- c(1,1,1,1) # Starting values for the chain
mesh.points <- 10 # Number of mesh points
m1 <- GQD.mle(Xt,time,mesh=mesh.points,theta=theta.start)
GQD.remove()
G0 <- function(t){theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))}
G1 <- function(t){-theta[1]}
Q1 <- function(t){theta[4]*theta[4]}
theta.start <- c(1,1,1,1) # Starting values for the chain
mesh.points <- 10 # Number of mesh points
m2 <- GQD.mle(Xt,time,mesh=mesh.points,theta=theta.start)
# Check estimates:
GQD.estimates(m1)
GQD.estimates(m2)
# Compare models:
GQD.aic(list(m1,m2))
#===============================================================================
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(DiffusionRgqd)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/DiffusionRgqd/GQD.mle.Rd_%03d_medium.png", width=480, height=480)
> ### Name: GQD.mle
> ### Title: MLEs for Generalized Quadratic Diffusion Models (GQDs).
> ### Aliases: GQD.mle
> ### Keywords: syntax C++ maximum likelihood
>
> ### ** Examples
>
> ## No test:
> #===============================================================================
> # Simulate a time inhomogeneous diffusion.
> #-------------------------------------------------------------------------------
>
> data(SDEsim1)
> attach(SDEsim1)
> par(mfrow=c(1,1))
> expr1=expression(dX[t]==2*(5+3*sin(0.5*pi*t)-X[t])*dt+0.5*sqrt(X[t])*dW[t])
> plot(Xt~time,type='l',col='blue',xlab='Time (t)',ylab=expression(X[t]),main=expr1)
>
> #------------------------------------------------------------------------------
> # Define coefficients of the process.
> #------------------------------------------------------------------------------
>
> GQD.remove()
[1] "Removed : NA "
> G0 <- function(t){theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))}
> G1 <- function(t){-theta[1]}
> Q0 <- function(t){theta[4]*theta[4]}
>
> theta.start <- c(1,1,1,1) # Starting values for the chain
> mesh.points <- 10 # Number of mesh points
>
> m1 <- GQD.mle(Xt,time,mesh=mesh.points,theta=theta.start)
Compiling C++ code. Please wait.
================================================================
Generalized Ornstein-Uhlenbeck
================================================================
_____________________ Drift Coefficients _______________________
G0 : theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))
G1 : -theta[1]
G2
___________________ Diffusion Coefficients _____________________
Q0 : theta[4]*theta[4]
Q1
Q2
_______________________ Model/Chain Info _______________________
Time Homogeneous : No
Data Resolution : Homogeneous: dt=0.25
# Removed Transits. : None
Density approx. : 2nd Ord. Truncation + Std Normal Dist.
Elapsed time : 00:00:00
... ... ... ... ... ... ... ... ... ... ...
dim(theta) : 4
----------------------------------------------------------------
>
> GQD.remove()
[1] "Removed : G0 G1 Q0"
>
> G0 <- function(t){theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))}
> G1 <- function(t){-theta[1]}
> Q1 <- function(t){theta[4]*theta[4]}
>
> theta.start <- c(1,1,1,1) # Starting values for the chain
> mesh.points <- 10 # Number of mesh points
>
> m2 <- GQD.mle(Xt,time,mesh=mesh.points,theta=theta.start)
Compiling C++ code. Please wait.
================================================================
Generalized Quadratic Diffusion (GQD)
================================================================
_____________________ Drift Coefficients _______________________
G0 : theta[1]*(theta[2]+theta[3]*sin(0.25*pi*t))
G1 : -theta[1]
G2
___________________ Diffusion Coefficients _____________________
Q0
Q1 : theta[4]*theta[4]
Q2
_______________________ Model/Chain Info _______________________
Time Homogeneous : No
Data Resolution : Homogeneous: dt=0.25
# Removed Transits. : None
Density approx. : 4 Ord. Truncation +4th Ord. Saddlepoint Appr.
Elapsed time : 00:00:00
... ... ... ... ... ... ... ... ... ... ...
dim(theta) : 4
----------------------------------------------------------------
>
> # Check estimates:
> GQD.estimates(m1)
Estimate Lower_95 Upper_95
theta[1] 2.087 1.805 2.369
theta[2] 5.016 4.907 5.124
theta[3] 2.898 2.737 3.059
theta[4] 1.138 1.052 1.223
> GQD.estimates(m2)
Estimate Lower_95 Upper_95
theta[1] 2.036 1.773 2.299
theta[2] 5.014 4.906 5.121
theta[3] 2.911 2.756 3.066
theta[4] 0.498 0.460 0.535
>
> # Compare models:
> GQD.aic(list(m1,m2))
Convergence p min.likelihood AIC BIC N
Model 1 0 4 246.5406 501.0812 517.0570 401
Model 2 0 4 225.0633 [=] 458.1267 [=] 474.1025 401
>
> #===============================================================================
> ## End(No test)
>
>
>
>
>
> dev.off()
null device
1
>