R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Generate the Transition Density of a Scalar Jump Generalized...

JGQD.density

R Documentation

Generate the Transition Density of a Scalar Jump Generalized Quadratic Diffusion (GQD).

Description

JGQD.density() approximates the transition density of a scalar generalized quadratic diffusion model (GQD). Given an initial value for the diffusion, Xs, the approximation is evaluated for all Xt at equispaced time-nodes given by splitting [s, t] into subintervals of length delt. JGQD.density() approximates transitional densities of jump diffusions of the form:

where

and

describes a Poisson process with jumps of the form:

arriving with intensity

subject to a jump distribition of the form:

Usage

JGQD.density(Xs = 4, Xt = seq(5, 8, 1/10), s = 0, t = 5, delt =1/100,
             Jdist = "Normal", Jtype = "Add", Dtype = "Saddlepoint",
             Trunc = c(8, 4), factorize = FALSE, factor.type = "Diffusion",
              beta, print.output = TRUE, eval.density = TRUE)

Arguments

`Xs`	Initial value of the process at time s.
`Xt`	Vector of values at which the transition density is to be evaluated over the trajectory of the transition density from time s to t.
`s`	The starting time of the process.
`t`	The time horizon up to and including which the transitional density is evaluated.
`delt`	Size of the time increments at which successive evaluations are made.
`Dtype`	Character string indicating the type of density approximation (see details) to use. Types: `'Saddlepoint'` and `'Edgeworth'` 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 and 8 with the constraint that `Trunc[1] >= Trunc[2]`. Default is `c(4,4)`.
`Jdist`	Valid entries are 'Normal', 'Exponential', 'Gamma' or 'Laplace'.
`Jtype`	Valid types are 'Add' or 'Mult'.
`factorize`	Should factorization be used (default = TRUE).
`factor.type`	Can be used to envoke a special factorization in order to evaluate Hawkes processes nested within the JGQD framework.
`beta`	Variable used for a special case jump structure (for research purposes).
`print.output`	If `TRUE`, model information is printed to the console.
`eval.density`	If `TRUE`, the density is evaluated in addition to calculating the moment eqns.

Details

vdfs

JGQD.density constructs an approximate transition density for a class of quadratic diffusion models. This is done by first evaluating the trajectory of the cumulants/moments of the diffusion numerically as the solution of a system of ordinary differential equations over a time horizon [s,t] split into equi-distant points delt units apart. Subsequently, the resulting cumulants/moments are carried into a density approximant (by default, a saddlepoint approximation) in order to evaluate the transtion surface.

Value

`density`	A matrix giving the density over the spatio-temporal mesh whose vertices are defined by paired permutations of the elements of `X_t` and `time`
`Xt`	A vector of points defining the state space at which the density was evaluated(recycled from input).
`time`	A vector of time points at which the density was evaluated.
`cumulants`	A matrix giving the cumulants of the diffusion. Row i gives the i-th cumulant.
`moments`	A matrix giving the moments of the diffusion. Row i gives the i-th cumulant.
`mesh`	A matrix giving the mesh used for normalization of the density.

Interface

DiffusionRjgqd uses a functional interface whereby th coefficients of a jump diffusion is defined by functions in the current workspace. By defining time-dependent functions with names that match the coefficients of the desired diffusion, DiffusionRjgqd reads the workspace and prepares the appropriate algorithm.

diffusiontable,width=100

In the case of jump diffusions, additional coefficients are required for the jump mechanism as well. Intensity coefficients and jump distributions, along with their corresponding R-names, are given in the tables below.

Intensity:

intensitytable

Jump distributions:

jumptable

Warning

Warning [1]: The system of ODEs that dictate the evolution of the cumulants do so approximately. Thus, although it is unlikely such cases will be encountered in inferential contexts, it is worth checking (by simulation) whether cumulants accurately replicate those of the target jump GQD. Furthermore, it may in some cases occur that the cumulants are indeed accurate whilst the density approximation fails. This can again be verified by simulation after which alternate density approximants may be specified through the variable Dtype.

Warning [2]: The parameter delt is also used as the stepsize for solving a system of ordinary differential equations (ODEs) that govern the evolution of the cumulants of the diffusion. As such delt is required to be small for highly non-linear models in order to ensure sufficient accuracy.

Author(s)

Etienne A.D. Pienaar: etiannead@gmail.com

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.

Examples



#===============================================================================
# Compare transition densities and mean trajectories of a non-linear jump
# diffusion for various jump distributions.
#-------------------------------------------------------------------------------
rm(list=ls(all=TRUE))
library(DiffusionRjgqd)

JGQD.remove()
# Define the jump diffusion using the DiffusionRjgqd syntax:
G1=function(t){0.2*5+0.1*sin(2*pi*t)}
G2=function(t){-0.2}
Q1=function(t){0.2}

# State dependent intensity:
Lam0 = function(t){1}
Lam1    = function(t){0.1}

# Normally distributed jumps: N(1,0.2)
Jmu    = function(t){1.0}
Jsig   = function(t){0.2}
# Normal distribution is the default:
res_1  = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,factorize=FALSE)


# Gamma distributed jumps: Gamma(0.5,0.5)
Jalpha = function(t){0.5}
Jbeta  = function(t){0.5}
# Jdist sets the jump distribution type:
res_2  = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Gamma',factorize=FALSE)

# Laplace jump parameters: Laplace(0.5*(sin(pi*t)),0.2)
Ja    = function(t){0.5*sin(pi*t)}
Jb    = function(t){0.2}
res_3 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)

par(mfrow=c(2,2))
persp(x=res_1$Xt,y=res_1$time,z=pmin(res_1$density,0.8),col=3,xlab='State (X_t)',ylab='Time (t)',
          zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
          main ='Normally distributed jumps',zlim=c(0,0.8))

persp(x=res_2$Xt,y=res_2$time,z=pmin(res_2$density,0.8),col=4,xlab='State (X_t)',ylab='Time (t)',
          zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
          main ='Gamma distributed jumps',zlim=c(0,0.8))

persp(x=res_3$Xt,y=res_3$time,z=pmin(res_3$density,0.8),col=5,xlab='State (X_t)',ylab='Time (t)',
          zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
          main ='Laplace distributed jumps',zlim=c(0,0.8))
plot(res_1$moments[1,]~res_1$time,type='n',main='Mean trajectories',ylab='E[X_t]',xlab='Time (t)')
lines(res_1$moments[1,]~res_1$time,col=1+2)
lines(res_2$moments[1,]~res_2$time,col=1+3)
lines(res_3$moments[1,]~res_3$time,col=1+4)

#===============================================================================
# Compare mean trajectories and zero-jump probabilities of a non-linear jump
# diffusion for various jump intensities.
#-------------------------------------------------------------------------------
JGQD.remove()
# Define the jump diffusion using the DiffusionRjgqd syntax:
G1=function(t){0.2*6}
G2=function(t){-0.2}
Q0=function(t){1.2}

# Laplace jump parameters: Laplace(05,0.2)
Ja    = function(t){0.5}
Jb    = function(t){0.2}


# Constant intensity:
Lam0 = function(t){1}
res_1 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)

# State dependent intensity:
Lam0 = function(t){0}
Lam1    = function(t){0.1}
res_2 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)

# State dependent, inhomogeneous intensity:
Lam0 = function(t){0}
Lam1    = function(t){0.1*(1+sin(5*pi*t))}
res_3 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)


par(mfrow=c(1,2))

plot(res_1$moments[1,]~res_1$time,type='n',main='Mean trajectories',ylab='E[X_t]',xlab='Time (t)')
lines(res_1$moments[1,]~res_1$time,col=2)
lines(res_2$moments[1,]~res_2$time,col=3)
lines(res_3$moments[1,]~res_3$time,col=4)

plot(res_1$zero_jump_prob~res_1$time,type='n',main=expression(P(N_t ==0)),ylab='Probability',
     xlab='Time (t)',ylim=c(0,1))
lines(res_1$zero_jump_prob~res_1$time,col=2)
lines(res_2$zero_jump_prob~res_2$time,col=3)
lines(res_3$zero_jump_prob~res_3$time,col=4)


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

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(DiffusionRjgqd)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/DiffusionRjgqd/JGQD.density.Rd_%03d_medium.png", width=480, height=480)
> ### Name: JGQD.density
> ### Title: Generate the Transition Density of a Scalar Jump Generalized
> ###   Quadratic Diffusion (GQD).
> ### Aliases: JGQD.density
> ### Keywords: transition density moments
> 
> ### ** Examples
> 
> ## No test: 
> 
> #===============================================================================
> # Compare transition densities and mean trajectories of a non-linear jump
> # diffusion for various jump distributions.
> #-------------------------------------------------------------------------------
> rm(list=ls(all=TRUE))
> library(DiffusionRjgqd)
> 
> JGQD.remove()
[1] "Removed :  NA "
> # Define the jump diffusion using the DiffusionRjgqd syntax:
> G1=function(t){0.2*5+0.1*sin(2*pi*t)}
> G2=function(t){-0.2}
> Q1=function(t){0.2}
> 
> # State dependent intensity:
> Lam0 = function(t){1}
> Lam1    = function(t){0.1}
> 
> # Normally distributed jumps: N(1,0.2)
> Jmu    = function(t){1.0}
> Jsig   = function(t){0.2}
> # Normal distribution is the default:
> res_1  = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*5+0.1*sin(2*pi*t)                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0                                                              
 Q1 : 0.2                                                        
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 1                                                        
 Lam1 : 0.1                                                      
 Lam2                                                            
 ........................... Jumps ..............................
 Normal                                                          
 Jmu : 1                                                         
 Jsig : 0.2                                                      
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
=================================================================
> 
> 
> # Gamma distributed jumps: Gamma(0.5,0.5)
> Jalpha = function(t){0.5}
> Jbeta  = function(t){0.5}
> # Jdist sets the jump distribution type:
> res_2  = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Gamma',factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*5+0.1*sin(2*pi*t)                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0                                                              
 Q1 : 0.2                                                        
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 1                                                        
 Lam1 : 0.1                                                      
 Lam2                                                            
 ........................... Jumps ..............................
 Gamma                                                           
 Jalpha : 0.5                                                    
 Jbeta : 0.5                                                     
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
=================================================================
> 
> # Laplace jump parameters: Laplace(0.5*(sin(pi*t)),0.2)
> Ja    = function(t){0.5*sin(pi*t)}
> Jb    = function(t){0.2}
> res_3 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*5+0.1*sin(2*pi*t)                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0                                                              
 Q1 : 0.2                                                        
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 1                                                        
 Lam1 : 0.1                                                      
 Lam2                                                            
 ........................... Jumps ..............................
 Laplace                                                         
 Ja : 0.5*sin(pi*t)                                              
 Jb : 0.2                                                        
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
[1] "yeah"
=================================================================
> 
> par(mfrow=c(2,2))
> persp(x=res_1$Xt,y=res_1$time,z=pmin(res_1$density,0.8),col=3,xlab='State (X_t)',ylab='Time (t)',
+           zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
+           main ='Normally distributed jumps',zlim=c(0,0.8))
> 
> persp(x=res_2$Xt,y=res_2$time,z=pmin(res_2$density,0.8),col=4,xlab='State (X_t)',ylab='Time (t)',
+           zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
+           main ='Gamma distributed jumps',zlim=c(0,0.8))
> 
> persp(x=res_3$Xt,y=res_3$time,z=pmin(res_3$density,0.8),col=5,xlab='State (X_t)',ylab='Time (t)',
+           zlab='Density f(X_t|X_s)',border=NA,shade=0.5,theta=145,
+           main ='Laplace distributed jumps',zlim=c(0,0.8))
> plot(res_1$moments[1,]~res_1$time,type='n',main='Mean trajectories',ylab='E[X_t]',xlab='Time (t)')
> lines(res_1$moments[1,]~res_1$time,col=1+2)
> lines(res_2$moments[1,]~res_2$time,col=1+3)
> lines(res_3$moments[1,]~res_3$time,col=1+4)
> 
> #===============================================================================
> # Compare mean trajectories and zero-jump probabilities of a non-linear jump
> # diffusion for various jump intensities.
> #-------------------------------------------------------------------------------
> JGQD.remove()
[1] "Removed :  G1 G2 Q1 Lam0 Lam1 Jmu Jsig Jalpha Jbeta Ja Jb"
> # Define the jump diffusion using the DiffusionRjgqd syntax:
> G1=function(t){0.2*6}
> G2=function(t){-0.2}
> Q0=function(t){1.2}
> 
> # Laplace jump parameters: Laplace(05,0.2)
> Ja    = function(t){0.5}
> Jb    = function(t){0.2}
> 
> 
> # Constant intensity:
> Lam0 = function(t){1}
> res_1 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*6                                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0 : 1.2                                                        
 Q1                                                              
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 1                                                        
 Lam1                                                            
 Lam2                                                            
 ........................... Jumps ..............................
 Laplace                                                         
 Ja : 0.5                                                        
 Jb : 0.2                                                        
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
[1] "yeah"
=================================================================
> 
> # State dependent intensity:
> Lam0 = function(t){0}
> Lam1    = function(t){0.1}
> res_2 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*6                                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0 : 1.2                                                        
 Q1                                                              
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 0                                                        
 Lam1 : 0.1                                                      
 Lam2                                                            
 ........................... Jumps ..............................
 Laplace                                                         
 Ja : 0.5                                                        
 Jb : 0.2                                                        
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
[1] "yeah"
=================================================================
> 
> # State dependent, inhomogeneous intensity:
> Lam0 = function(t){0}
> Lam1    = function(t){0.1*(1+sin(5*pi*t))}
> res_3 = JGQD.density(4,seq(2,10,1/10),0,2.5,1/100,Jdist = 'Laplace',factorize=FALSE)
                                                                 
 ================================================================
            Jump Generalized Quadratic Diffusion (JGQD)          
 ================================================================
 _____________________ Drift Coefficients _______________________
 G0                                                              
 G1 : 0.2*6                                                      
 G2 : -0.2                                                       
 ___________________ Diffusion Coefficients _____________________
 Q0 : 1.2                                                        
 Q1                                                              
 Q2                                                              
 _______________________ Jump Mechanism _________________________
 ......................... Intensity ............................
 Lam0 : 0                                                        
 Lam1 : 0.1*(1+sin(5*pi*t))                                      
 Lam2                                                            
 ........................... Jumps ..............................
 Laplace                                                         
 Ja : 0.5                                                        
 Jb : 0.2                                                        
 __________________ Distribution Approximant ____________________
 Density approx. : Saddlepoint                                   
 Trunc. Order    : 8                                             
 Dens.  Order    : 4                                             
[1] "yeah"
=================================================================
> 
> 
> par(mfrow=c(1,2))
> 
> plot(res_1$moments[1,]~res_1$time,type='n',main='Mean trajectories',ylab='E[X_t]',xlab='Time (t)')
> lines(res_1$moments[1,]~res_1$time,col=2)
> lines(res_2$moments[1,]~res_2$time,col=3)
> lines(res_3$moments[1,]~res_3$time,col=4)
> 
> plot(res_1$zero_jump_prob~res_1$time,type='n',main=expression(P(N_t ==0)),ylab='Probability',
+      xlab='Time (t)',ylim=c(0,1))
> lines(res_1$zero_jump_prob~res_1$time,col=2)
> lines(res_2$zero_jump_prob~res_2$time,col=3)
> lines(res_3$zero_jump_prob~res_3$time,col=4)
> 
> 
> #===============================================================================
> ## End(No test)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>