Likelihood Inference for Non-Linear Multivariate Jump Diffusions with State Dependent Intensity
Example replication materials
Etienne A.D. Pienaar
February, 2016
Readme
This file provides replication materials for examples and analysis conducted in the paper “Likelihood Inference for Non-Linear Multivariate Jump Diffusions with State Dependent Intensity”, which develops the methodology on which the DiffusionRjgqd package is built. As such, this file may not be useful as a stand-alone document, and should be read in conjunction with the paper. Note further:
- For users that find the interface difficult to understand it may be useful to read material on the ralated package: DiffusionRgqd and the introductory vignettes therein.
In order to save time and keep this document readable, we reffrain from producing output for some of the examples.
Some of the examples provided here are given in greater detail within the package vignettes (start here, or skip to the further reading section).
Approximating the transitional density of a jump diffusion
A scalar example (Section 3.3)
library(DiffusionRjgqd)
# Define the jump diffusion in DiffusionRjgqd syntax:
JGQD.remove()## [1] "Removed : NA "G0=function(t){2*5}
G1=function(t){-2}
Q1=function(t){1*(1+0.4*sin(pi*t))^2}
# Set the transition rates and state values for CTMC:
lambda_1 = 1
lambda_2 = 3
beta_1 = 0.25
beta_2 = 1
Jmu = function(t){1}
Jsig = function(t){0.25}
Lam0 = function(t){lambda_1*(beta_2+beta_1*exp(-(beta_1+beta_2)*t))/(beta_1+beta_2)+lambda_2*beta_1/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*t))}
# Define the jump diffusion coefficients for simulation:
mu = function(x,t){G0(t)+G1(t)*x}
sigma = function(x,t){sqrt(Q1(t)*x)}
j = function(x,z){z}
simulate=function(x0=4,N=10000,pts =c(1,2,3,4,5),init.state=1)
{
d = 0 # Time index
delta = 1/1000 # Step size
tt = seq(0,5,delta) # Equispaced time points
X = rep(x0,N) # Initialize state vector
states= rep(c(lambda_1,lambda_2)[init.state],N) # Initialize intensity process
# Storage for record of the evolution of the moments:
moments = matrix(0,4,length(tt))
cumulants = moments
moments[,1] = x0^{1:4}
cumulants[1,1] = x0
# Storage of snapshots of the simulated trajectories:
L = list()
count =1
for(i in 2:length(tt))
{
X=X+mu(X,d)*delta+sigma(X,d)*rnorm(N,sd=sqrt(delta))
# Simulate the occurance of a jump event:
events = (1-exp(-Lam0(d)*delta)>runif(N))
wh=which(events)
# For those trajectories that events do occur, simulate a jumps
# half a step forward:
if(any(events))
{
X[wh]=X[wh]+j(X[wh],rnorm(length(wh),Jmu(d),Jsig(d)))
}
d=d+delta
# Calculate the prob of a change for a given state:
probs1=beta_1/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*delta))
probs2=beta_2/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*delta))
# For a given state, simulate transitions (or remain in a given state):
which1=which(states==lambda_1)
which2=which(states==lambda_2)
if(length(which1)!=0)
{
whh1.2=which(runif(length(which1))<probs1)
if(length(whh1.2)!=0)
{
states[which1][whh1.2]=lambda_2
}
}
if(length(which2)!=0)
{
whh2.2=which(runif(length(which2))<probs2)
if(length(whh2.2)!=0)
{
states[which2][whh2.2]=lambda_1
}
}
# Record the moments of the process:
S1=sum(X);S2=sum(X^2);S3=sum(X^3);S4=sum(X^4)
moments[,i]=c(S1,S2,S3,S4)/N
cumulants[1,i] = S1/N
cumulants[2,i] = 1/(N*(N-1))*(N*S2-S1^2)
cumulants[3,i] = 1/(N*(N-1)*(N-2))*(2*S1^3-3*N*S1*S2+N^2*S3)
cumulants[4,i] = (1/(N*(N-1)*(N-2)*(N-3))*(-6*S1^4+12*N*S1^2*S2-3*N*(N-1)*S2^2
-4*N*(N+1)*S1*S3+N^2*(N+1)*S4))
# Take snapshots at pts[count]:
if(sum(pts==round(d,3))!=0)
{
L[[count]] = hist(X,plot=F,breaks=25)
count=count+1
}
}
return(list(moments=moments,cumulants=cumulants,time=tt,hists=L,pts=pts))
}
# Simulate the process:
res.sim=simulate()
# Calculate the approximate transition density:
res = JGQD.density(4,seq(2,14,1/10),0,5,1/100)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 1*(1+0.4*sin(pi*t))^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : lambda_1*(beta_2+beta_1*exp(-(beta_1+beta_2)*t))/(beta_1+beta_2)+lambda_2*beta_1/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*t))
## Lam1
## Lam2
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.25
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= # Compare moments:
exprs =c(expression(log(m[1](t))),expression(log(m[2](t))),expression(log(m[3](t)))
,expression(log(m[4](t))))
par(mfrow=c(1,1))
for(i in 1:4)
{
plot(log(res.sim$moments[i,])~res.sim$time,type='l',ylab = exprs[i],xlab='Time (t)'
,lwd=2,col="#BBCCEE")
lines(log(res$moments[i,])~res$time,lty='dashed',col='#222299',lwd=2)
}
legend('bottomright',lty=c('solid','dashed'),col=c('#BBCCEE','#222299'),
legend=c('Simulated','Moment eqns'),lwd=2,bty='n')
library(rgl)
um = rbind(
c(-0.7762160, 0.6299503, 0.02552053, 0),
c(-0.2308982, -0.3217098, 0.91825324, 0),
c(0.5866640, 0.7068703, 0.39517075, 0),
c(0.0000000, 0.0000000, 0.00000000, 1))
r3dDefaults$userMatrix =um
open3d(windowRect=c(0,0,360*1.5,360*1.5)+30,zoom=6/8)## wgl
## 1 persp3d(res$Xt,res$time,pmax(pmin(res$density,1),0),col='white',alpha=0.5,box=F,
xlab='Xt',ylab='Time',zlab='')
surface3d(res$Xt,res$time[-c(1,201:501)],pmax(pmin(res$density[,2:200],1),0),col='white')
cols=colorRampPalette(c("red", "yellow"))
library(colorspace)
colpal=function(n){rev(sequential_hcl(n,power=0.8,l=c(40,100)))}
for(i in 2:5)
{
h1 =res.sim$hists[[i]]
y=rep(h1$density,each=2)
x=c(rbind(h1$breaks[-length(h1$breaks)],h1$breaks[-1]))
hd=cbind(0,y,0)
tt=res.sim$pts[i]
surface3d(x,c(tt-0.0001,tt,tt+0.0001),hd,col=colpal(5)[i],alpha=1)
lines3d(res$Xt,tt,res$density[,i*100],col='black',lwd=2)
}
rgl.snapshot('temp.png')
library(png)
imag = readPNG(paste0(getwd(),'/temp.png'))
par(mfrow=c(1,1))
plot(1:2, type='n', main="", xlab="", ylab="",axes=FALSE)
lim <- par()
rasterImage(imag, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
# Generate a comparison for starting in the high intensity regime
Lam0 = function(t){lambda_1*beta_2/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*t))+lambda_2*(beta_2+beta_1*exp(-(beta_1+beta_2)*t))/(beta_1+beta_2)}
# Simulate the process:
res2.sim=simulate(init.state=2)
# Calculate the approximate transition density:
res2 = JGQD.density(4,seq(2,14,1/10),0,5,1/100)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 1*(1+0.4*sin(pi*t))^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : lambda_1*beta_2/(beta_1+beta_2)*(1-exp(-(beta_1+beta_2)*t))+lambda_2*(beta_2+beta_1*exp(-(beta_1+beta_2)*t))/(beta_1+beta_2)
## Lam1
## Lam2
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.25
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= persp3d(res$Xt,res$time,pmax(pmin(res$density,1),0),col='white',alpha=0.9,box=F,
xlab='Xt',ylab='Time',zlab='')
surface3d(res$Xt,res$time[-c(1,201:501)],pmax(pmin(res$density[,2:200],1),0),col='white')
surface3d(res2$Xt,res2$time,pmax(pmin(res2$density,1),0),col='lightblue',alpha=0.9)
surface3d(res2$Xt,res2$time[-c(1,201:501)],pmax(pmin(res2$density[,2:200],1),0),col='lightblue')
rgl.snapshot('temp.png')
library(png)
imag = readPNG(paste0(getwd(),'/temp.png'))
plot(1:2, type='n', main="", xlab="", ylab="",axes=FALSE)
lim <- par()
rasterImage(imag, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
Short time-scale dynamics and the excess factorization (Section 3.4)
library(DiffusionRjgqd)
# Define the process in DiffusionRjgqd syntax:
JGQD.remove()## [1] "Removed : G0 G1 Q1 Lam0 Jmu Jsig"G0=function(t){2*5}
G1=function(t){-2}
Q1=function(t){0.25}
Jmu = function(t){1.0}
Jsig = function(t){0.5}
Lam0 = function(t){0}
Lam1 = function(t){0.5*(1+sin(3*pi*t))}
Lam2 = function(t){0.1*(1+cos(3*pi*t))}
# Define the jump diffusion coefficients for simulation:
mu = function(x,t){G0(t)+G1(t)*x}
sigma = function(x,t){sqrt(Q1(t)*x)}
j = function(x,z){z}
lambd = function(x,t){(Lam0(t)+Lam1(t)*x+Lam2(t)*x^2)}
simulate=function(x0=4,TT=1.15,N=50000,pts =c(1,2,3,4,5),brks=25)
{
d=0 # Time index
delta=1/2000 # Step size
tt=seq(0,TT,delta) # Equispaced points on [0,TT]
X=rep(x0,N) # Initialize the state vector
isjump = rep(0,N) # Used for counting the number of jumps
probs=matrix(1,3,length(tt)) # Used to store probabilities
# Storage of snapshots of the simulated trajectories:
L = list()
count =1
for(i in 2:length(tt))
{
# Simulate the diffuse part:
X=X+mu(X,d)*delta+sigma(X,d)*rnorm(N,sd=sqrt(delta))
# Simulate the occurance of a jump event:
events = (1-exp(-lambd(X,d)*delta)>runif(N))
wh=which(events)
whn = which(!events)
# For those trajectories that events do occur, simulate a jumps:
if(any(events))
{
X[wh]=X[wh]+j(X[wh],rnorm(length(wh),Jmu(d),Jsig(d)))
}
d=d+delta
isjump = isjump +events
probs[1,i] = sum(isjump==0)/N
probs[2,i] = sum(isjump==1)/N
probs[3,i] = sum(isjump==2)/N
# Take snapshots at pts[count]:
if(sum(pts==round(d,4))!=0)
{
L[[count]] = hist(X,plot=F,breaks=brks)
count=count+1
}
}
return(list(probs=probs,time =tt,X=X,hists=L,pts=pts))
}
res2a=simulate()
resa=JGQD.density(Xs=4,Xt=seq(0,20,1/10),s=0,t=1.15,delt=1/100)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 0.25
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : 0
## Lam1 : 0.5*(1+sin(3*pi*t))
## Lam2 : 0.1*(1+cos(3*pi*t))
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.5
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= Lam1=function(t){0.2}
Lam2=function(t){0}
res2b=simulate()
resb=JGQD.density(Xs=4,Xt=seq(0,20,1/10),s=0,t=1.15,delt=1/100)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 0.25
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : 0
## Lam1 : 0.2
## Lam2 : 0
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.5
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= pp1 = resa$zero_jump_prob
pp2 = resb$zero_jump_prob
epr1 = expression(X[t])
epr2 = expression(d/dt*log(P(N[t]-N[s]==0))==E[X](lambda(X[t],t)))
# Plot the evolution of the 0 jump probabilities
plot(pp1~resa$time,type='l',ylim=c(0,1),lwd=1,col='black',
main='Probability of observing 0 jumps',ylab=expression(P(N[t]-N[0] == 0)),
xlab='t')#,axes=F)
lines(res2a$probs[1,]~res2a$time,col='blue',lty='dashed',lwd=2)
lines(pp2~resa$time,col='black',lty='solid',lwd=1)
lines(res2b$probs[1,]~res2b$time,col='blue',lty='dashed',lwd=2)
abline(h=c(0.0,0.5,0.5),lty='dotted',col='grey')
legend('topright',lty=c('solid','dashed'),lwd=c(1,2),col=c('black','blue'),
legend=c('ODE soln.','Simulated'),bty='n')
text(0.5,0.7,label=expression(lambda(X[t],t) == 0.2*X[t]),pos=4,cex=0.85)
text(0.25,0.25,label=expression(lambda(X[t],t) ==
0.5*(1+sin(3*pi*t))*X[t]+0.1*(1+cos(3*pi*t))*X[t]^2),pos=4,cex=0.85)
# Evaluate the transition densities over short transition horizons:
TT=0.02
Lam1 = function(t){0.5*(1+sin(3*pi*t))}
Lam2 = function(t){0.1*(1+cos(3*pi*t))}
res2a=simulate(TT=TT,pts=seq(0.5*TT,TT,length=4))
resa=JGQD.density(Xs=4,Xt=seq(0,10,1/100),s=0,t=TT,delt=1/200,factorize=TRUE)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 0.25
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : 0
## Lam1 : 0.5*(1+sin(3*pi*t))
## Lam2 : 0.1*(1+cos(3*pi*t))
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.5
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= hist(res2a$X,freq=F,col='gray74',breaks=50,main=paste0('Transition density at t = ',
TT),xlab=epr1,xlim=c(3.5,8),border='white')
lines(resa$density[,TT*200]~resa$Xt,col='darkblue',lwd=2)
legend('topright',lty=c('solid','solid'),lwd=c(2,5),col=c('darkblue','gray74'),
legend=c('Approximation','Simulated'),bty='n')
box()
TT=0.1
Lam1=function(t){0.2}
Lam2=function(t){0}
res2b=simulate(TT=TT)
resb=JGQD.density(Xs=4,Xt=seq(0,10,1/50),s=0,t=TT,delt=1/100,factorize=TRUE)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 0.25
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : 0
## Lam1 : 0.2
## Lam2 : 0
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.5
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= hist(res2b$X,freq=F,col='gray74',breaks=50,main=paste0('Transition density at t = ',
TT),xlab=epr1,xlim=c(3,8),border='white')
lines(resb$density[,TT*100]~resb$Xt,col='darkblue',lwd=2)
legend('topright',lty=c('solid','solid'),lwd=c(2,5),col=c('darkblue','gray74'),
legend=c('Approximation','Simulated'),bty='n')
box()
TT=0.1
Lam1=function(t){0.2}
Lam2=function(t){0}
Q1 =function(t){0.25}
Jsig = function(t){0.05}
res.sim=simulate(TT=TT,pts=seq(0.04,0.1,length=3),brks=55)
res =JGQD.density(Xs=4,Xt=seq(3,7,1/50),s=0,t=TT,delt=1/200,factorize=TRUE)##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 2*5
## G1 : -2
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : 0.25
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : 0
## Lam1 : 0.2
## Lam2 : 0
## ........................... Jumps ..............................
## Normal
## Jmu : 1
## Jsig : 0.05
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## ================================================================= library(rgl)
um = rbind(
c(-0.7762160, 0.6299503, 0.02552053, 0),
c(-0.2308982, -0.3217098, 0.91825324, 0),
c(0.5866640, 0.7068703, 0.39517075, 0),
c(0.0000000, 0.0000000, 0.00000000, 1))
r3dDefaults$userMatrix =um
open3d(windowRect=c(0,0,360*1.5,360*1.5)+30,zoom=6/8)## wgl
## 2 persp3d(res$Xt,res$time,pmax(pmin(res$density,3.5),0),col='white',alpha=0.4,box=F,
xlab='Xt',ylab='Time',zlab='')
surface3d(res$Xt,res$time[1:8],pmax(pmin(res$density[,1:8],3.5),0),col='white')
cols=colorRampPalette(c("red", "yellow"))
library(colorspace)
colpal=function(n){rev(sequential_hcl(n,power=0.8,l=c(40,100)))}
for(i in 1:length(res.sim$pts))
{
h1 =res.sim$hists[[i]]
y=rep(h1$density,each=2)
x=c(rbind(h1$breaks[-length(h1$breaks)],h1$breaks[-1]))
hd=cbind(0,y,0)
tt=res.sim$pts[i]
surface3d(x,c(tt-0.0001,tt,tt+0.0001),hd,col=colpal(5)[i+1],alpha=1)
lines3d(res$Xt,tt,res$density[,tt*200],col='black',lwd=2)
}
rgl.snapshot('temp.png')
library(png)
imag = readPNG(paste0(getwd(),'/temp.png'))
par(mfrow=c(1,1))
plot(1:2, type='n', main="", xlab="", ylab="",axes=FALSE)
lim <- par()
rasterImage(imag, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
A bivariate CIR process with a time inhomogeneous jump mechanism (Section 3.4)
See also:
JGQD.remove()## [1] "Removed : G0 G1 Q1 Lam0 Lam1 Lam2 Jmu Jsig" theta=c(0.5,5,-0.1,0.4,0.4,6,-0.1,0.3,1,0.75,0.75,0.75,0.75)
a00 = function(t){theta[1]*theta[2]}
a10 = function(t){-theta[1]}
a01 = function(t){theta[3]}
c10 = function(t){theta[4]^2}
b00 = function(t){theta[5]*theta[6]}
b01 = function(t){-theta[5]}
b10 = function(t){theta[7]}
f01 = function(t){theta[8]^2}
Lam00= function(t){theta[9]}
Jmu1=function(t){theta[10]*(1+sin(2*pi*t))}
Jmu2=function(t){theta[11]*(1+sin(2*pi*t))}
Jsig11=function(t){theta[12]^2*(1+0.8*sin(2*pi*t))^2}
Jsig22=function(t){theta[13]^2*(1+0.8*sin(2*pi*t))^2}
xx=seq(3,11,1/10)
yy=seq(3,11,1/10)
res= BiJGQD.density(7,7,xx,yy,0,1,1/100,Dtype='Saddlepoint')##
## ================================================================
## GENERALIZED QUADRATIC DIFFUSON
## ================================================================
## _____________________ Drift Coefficients _______________________
## a00 : theta[1]*theta[2]
## a10 : -theta[1]
## a01 : theta[3]
## ... ... ... ... ... ... ... ... ... ... ...
## b00 : theta[5]*theta[6]
## b10 : theta[7]
## b01 : -theta[5]
## ___________________ Diffusion Coefficients _____________________
## c10 : theta[4]^2
## ... ... ... ... ... ... ... ... ... ... ...
## ... ... ... ... ... ... ... ... ... ... ...
## ... ... ... ... ... ... ... ... ... ... ...
## f01 : theta[8]^2
## _______________________ Jump Components ________________________
## ......................... Intensity ............................
## Lam00 : theta[9]
## ........................... Jumps ..............................
## Jmu1 : theta[10]*(1+sin(2*pi*t))
## Jmu2 : theta[11]*(1+sin(2*pi*t))
## Jsig11 : theta[12]^2*(1+0.8*sin(2*pi*t))^2
## Jsig22 : theta[13]^2*(1+0.8*sin(2*pi*t))^2
## ================================================================= # Simulate the process
mux = function(x,y,t){a00(t)+a10(t)*x+a01(t)*y}
sigmax = function(x,y,t){sqrt(c10(t)*x)}
muy = function(x,y,t){b00(t)+b10(t)*x+b01(t)*y}
sigmay = function(x,y,t){sqrt(f01(t)*y)}
lambda1 = function(x,y,t){Lam00(t)}
lambda2 = function(x,y,t){rep(0,length(x))}
j11 = function(x,y,z){z}
j12 = function(x,y,z){z}
j21 = function(x,y,z){z}
j22 = function(x,y,z){z}
simulate=function(x0=7,y0=7,N=10000,TT=5,delta=1/1000,pts,brks=30,plt=F)
{
library(colorspace)
colpal=function(n){rev(sequential_hcl(n,power=0.8,l=c(40,100)))}
d=0 # Time index
tt=seq(0,TT,delta) # Time sequance
X=rep(x0,N) # Initialize state vectors
Y=rep(y0,N)
x.traj = rep(x0,length(tt))
y.traj = rep(y0,length(tt))
x.jump = rep(0,length(tt))
y.jump = rep(0,length(tt))
# Storage for histogram snapshots:
count = 1
L1=list()
L2=list()
evts = rep(0,N)
for(i in 2:length(tt))
{
X=X+mux(X,Y,d)*delta+sigmax(X,Y,d)*rnorm(N,sd=sqrt(delta))
Y=Y+muy(X,Y,d)*delta+sigmay(X,Y,d)*rnorm(N,sd=sqrt(delta))
events1 = (lambda1(X,Y,d)*delta>runif(N))
if(any(events1))
{
wh=which(events1)
evts[wh]=evts[wh]+1
X[wh]=X[wh]+j11(X[wh],Y[wh],rnorm(length(wh),Jmu1(d),sqrt(Jsig11(d))))
Y[wh]=Y[wh]+j21(X[wh],Y[wh],rnorm(length(wh),Jmu2(d),sqrt(Jsig22(d))))
}
events2 = (lambda2(X,Y,d)*delta>runif(N))
d=d+delta
if(sum(round(pts,3)==round(d,3))!=0)
{
if(plt)
{
expr1 = expression(X_t)
expr2 = expression(Y_t)
color.palette=colorRampPalette(c('green','blue','red'))
filled.contour(res$Xt,res$Yt,res$density[,,i],
main=paste0('Transition Density \n (t = ',round(d,2),')'),
color.palette=colpal,
nlevels=41,xlab=expression(X[t]),ylab=expression(Y[t]),plot.axes=
{
# Add simulated trajectories
points(Y~X,pch=c(20,3)[(evts>0)+1],col=c('black','red')[(evts>0)+1],cex=c(0.9,0.6)[(evts>0)+1])
if(any(events2))
{
wh=which(events2)
segments(xpreee[wh],ypreee[wh],X[wh],Y[wh],col='gray')
}
axis(1);axis(2);
# Add a legend
legend('topright',col=c('black','red'),pch=c(20,3),
legend=c('Simulated Trajectories','Jumped'))
yy=contourLines(res$Xt,res$Yt,res$density[,,i],levels=seq(0.01,0.1,length=10))
if(length(yy)>0)
{
for(j in 1:length(yy))
{
lines(yy[[j]])
}
}
})
}
L1[[count]] = hist(X,plot=F,breaks=brks)
L2[[count]] = hist(Y,plot=F,breaks=brks)
count=count+1
#savePlot(paste0('BiExampleTD',count,'.pdf'),type='pdf')
}
}
return(list(time=tt,histsx=L1,histsy=L2,pts=pts))
}
sim=simulate(7,7,N=200,TT=0.75,delta=1/100,plt=T,pts=c(0.13,0.28,0.38,0.51,0.63,0.75))
Comparison to existing methods
See also:
Brownian motion with drift: Short transition horizon approximations (Section 4.1)
rm(list=ls(all=TRUE))
#---------------------------------------------------------------------------
# Part 1: True density, series approximation and saddlepoint comparison
#---------------------------------------------------------------------------
library(DiffusionRjgqd)
JGQD.remove()## [1] "Removed : NA "theta1=c(0.1,0.5,0.5,0.2,0.25)
theta2=c(0.0,0.25,0.5,0.0,0.25)
theta3=c(-0.1,0.75,0.5,0.5,0.1)
pars =rbind(theta1,theta2,theta3)
# Define the jump diffusion using the DiffusionRjgqd syntax:
G0=function(t){theta[1]}
Q0=function(t){theta[2]^2}
# State dependent intensity:
Lam0 = function(t){theta[3]}
Jmu = function(t){theta[4]}
Jsig = function(t){theta[5]}
true.density=function(x,y,t,theta,order =1)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
mu2 = theta[4]
sig2 = theta[5]
dens=exp(-lam*t)*dnorm(y,x+mu*t,sqrt(sig^2*t))
for(i in 1:order)
{
dens= dens +exp(-lam*t)*(lam*t)^i/factorial(i)*dnorm(y,x+mu*t+i*mu2,sqrt(sig^2*t+i*sig2^2))
}
return(list(density=dens,Xt=y))
}
herm.density=function(x,y,t,theta,type =1)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
mu2 = theta[4]
sig2 = theta[5]
if(type == 1)
{
dens=1/sqrt(2*pi*sig^2*t)*exp(-(y-x)^2/(2*sig^2*t)+mu/sig^2*(y-x))*(1-(mu^2/(2*sig^2)+lam)*t)+lam*t/sqrt(2*pi*sig2^2)*exp(-(y-x-mu2)^2/(2*sig2^2))
}
if(type ==2)
{
dens=dnorm(y,x+mu*t,sqrt(sig^2*t))*(1-lam*t)+lam*t*dnorm(y,x+mu2,sqrt(sig2^2))
}
return(list(density=dens,Xt=y))
}
#================================ Short time =======================================
ind=10
for(i in 1:3)
{
theta = pars[i,]
res_1 = JGQD.density(0,seq(-1,1,1/100),0,1,1/100,factorize=T,Jdist='Normal',
Dtype='Normal.A')
res_2 = true.density(0,res_1$Xt,ind*1/100,theta,order=10)
res_3 = herm.density(0,res_1$Xt,ind*1/100,theta,type=1)
expr1 = expression(f[SPT]^{(4)})
expr2 = expression(f[true]^{(10)})
expr3 = expression(f[ser])
eprs =c(expr1,expr2,expr3)
ltys = c('solid','dotted','solid')
cols = c('blue' ,'black','#5AAE61')
lwds = c(1,2,1)
pchs= c(1,NA,3)
thin=seq(1,length(res_1$Xt),10)
plot((res_2$density)~res_2$Xt,type='l',col=cols[2],lty=ltys[2],lwd=lwds[2],
main=paste0('Brownian motion with drift (t = ',round(ind*1/100,2),')'),xlab=expression(X[t]),
ylab='Density',las=1)
lines((res_1$density[,ind])~res_1$Xt,col=cols[1],lty=ltys[1],lwd=lwds[1])
points((res_1$density[thin,ind])~res_1$Xt[thin],col=cols[1],pch=pchs[1],cex=1)
lines((res_3$density)~res_3$Xt,col=cols[3],lty=ltys[3],lwd=lwds[3])
points((res_3$density[thin])~res_3$Xt[thin],col=cols[3],pch=pchs[3],cex=0.9)
legend('topright',legend=eprs,lty=ltys,lwd=lwds,col=cols,pch=pchs,bty='n',cex=1.25,pt.cex=1)
}##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Normal
## Jmu : theta[4]
## Jsig : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Normal
## Jmu : theta[4]
## Jsig : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Normal
## Jmu : theta[4]
## Jsig : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
#================================ Error at short ===========================
par(mfrow=c(1,1))
ind=10
lams= seq(0,2,1/10)
mxs1 = rep(0,length(lams))
mxs2 = rep(0,length(lams))
expr1 = expression(f[SPT]^{(4)})
expr2 = expression(f[ser])
eprs =c(expr1,expr2)
ltys = c('solid','solid')
cols = c('blue' ,'#5AAE61')
lwds = c(1,1)
pchs= c(1,3)
plot(mxs1~lams,type='n',ylim=c(0,0.05),xlab=expression(lambda),ylab=expression(max(abs(hat(f)-f[true]^{(10)}))),main='Max abs. error for varying intensity - Par. set 2.')
for(j in 2){
for(i in 1:length(lams))
{
theta = pars[j,]
theta[3] = lams[i]
res_1 = JGQD.density(0,seq(-1.5,1.5,1/100),0,0.25,1/100,factorize=T,Jdist='Normal',
Dtype='Normal.A',print.output=FALSE)
res_5 = true.density(0,res_1$Xt,ind*1/100,theta,order=10)
res_2 = herm.density(0,res_1$Xt,ind*1/100,theta,type=1)
mxs1[i] = max(abs(res_5$density-res_1$density[,ind]))
mxs2[i] = max(abs(res_5$density-res_2$density))
}
sq=seq(1,length(lams),5)
lines(mxs1~lams,col=cols[1],lwd=lwds[1],lty=ltys[1],pch=pchs[1])
points(mxs1[sq]~lams[sq],col=cols[1],lwd=lwds[1],lty=ltys[1],pch=pchs[1])
lines(mxs2~lams,col=cols[2],lwd=lwds[2],lty=ltys[2])
points(mxs2[sq]~lams[sq],col=cols[2],lwd=lwds[2],lty=ltys[2],pch=pchs[2])
}## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================
## =================================================================legend('topleft',legend=eprs,lty=c(ltys),lwd=c(lwds,1,1,1),
col=c(cols),bty='n',cex=1.25,pt.cex=1,pch=pchs)
#---------------------------------------------------------------------------
# Part 2: Simulated, convolution, series approximation and saddlepoint
# comparison for Exponential and Laplace jumps
#---------------------------------------------------------------------------
rm(list=ls(all=TRUE))
library(DiffusionRjgqd)
JGQD.remove()## [1] "Removed : NA "Xs = 1
tau = 0.1
theta1=c(-0.5,0.75,0.5,1);
theta2=c(0,0.25 ,0.5,2);
theta3=c(0.5,0.5,2,5);
theta4=c(0,0.75,1,3);
pars =rbind(theta1,theta2,theta3,theta4)
# Define the jump diffusion using the DiffusionRjgqd syntax:
G0=function(t){theta[1]}
Q0=function(t){theta[2]^2}
# State dependent intensity:
Lam0 = function(t){theta[3]}
Jlam = function(t){1/theta[4]}
sim=function(tau,theta)
{
N= 100000
X=rep(Xs,N)
dt=1/2000
tt = seq(0,tau,dt)
for(i in 2:length(tt))
{
X = X +theta[1]*dt+theta[2]*rnorm(N,0,sqrt(dt))
events = (theta[3]*dt>runif(N))
if(any(events))
{
wh=which(events)
X[wh]=X[wh]+rexp(length(wh),theta[4])
}
}
return(list(X=X))
}
herm.density=function(x,y,t,theta,type =1)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
lam.z = theta[4]
if(type == 1)
{
dens=1/sqrt(2*pi*sig^2*t)*exp(-(y-x)^2/(2*sig^2*t)+mu/sig^2*(y-x))*(1-(mu^2/(2*sig^2)+lam)*t)+lam*t*dexp(y-x,lam.z)
}
if(type ==2)
{
dens=dnorm(y,x+mu*t,sqrt(sig^2*t))*(1-lam*t)+lam*t*dexp(y-x,lam.z)
}
return(list(density=dens,Xt=y))
}
conv.density=function(x,y,t,theta)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
lam.z = theta[4]
dens = dnorm(y,x+mu*t,sqrt(sig^2*t))*(1-lam*t)+lam*t*lam.z*exp((y-x-mu*t-lam.z*sig^2*t)^2/(2*sig^2*t)-(y-x-mu*t)^2/(2*sig^2*t))*(1-pnorm(0,y-x-mu*t-lam.z*sig^2*t,sqrt(sig^2*t)))
return(list(density=dens,Xt=y))
}
for(i in 1:4)
{
theta=pars[i,]
res_1 = JGQD.density(Xs,seq(0,2.5,1/100),0,1,1/100,factorize=T,Jdist='Exponential',
Dtype='Normal.A')
res_2 = herm.density(Xs,res_1$Xt,tau,theta,type=1)
res_3 = conv.density(Xs,res_1$Xt,tau,theta)
res.sim = sim(tau,theta)
expr1 = expression(f[SPT]^{(4)})
expr2 = expression(f[ser])
expr3 = expression(tilde(f))
eprs =c(expr1,expr2,expr3)
ltys = c('solid','solid','dashed')
cols = c('blue','#5AAE61','black')
lwds = c(1,1,1)
pchs= c(1,3,NA)
thin=seq(1,length(res_1$Xt),20)
hist(res.sim$X,freq=FALSE,border='white',main=paste0('Brownian motion / Exponential jumps (t = ',tau,')'),
breaks=seq(-25,25,1/20),xlab=expression(X[t]),ylab='Density',las=1,xlim=range(res_1$Xt),ylim=c(0,max(res_1$density[,tau*100])*1.2),,col='gray74')
lines((res_1$density[,tau*100])~res_1$Xt,type='l',col=cols[1],lty=ltys[1],lwd=lwds[1])
points((res_1$density[thin,tau*100])~res_1$Xt[thin],col=cols[1],pch=pchs[1],cex=0.9)
lines((res_2$density)~res_2$Xt,col=cols[2],lty=ltys[2],lwd=lwds[2])
points((res_2$density[thin])~res_2$Xt[thin],col=cols[2],pch=pchs[2],cex=0.9)
lines((res_3$density)~res_3$Xt,col=cols[3],lty=ltys[3],lwd=lwds[3])
legend('topright',legend=eprs,lty=ltys,lwd=lwds,col=cols,pch=pchs,bty='n',cex=1.25,pt.cex=1)
box()
}##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/theta[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/theta[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/theta[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/theta[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
rm(list=ls(all=TRUE))
library("smoothmest")## Warning: package 'smoothmest' was built under R version 3.2.4## Loading required package: MASSJGQD.remove()## [1] "Removed : NA "Xs = 1
tau = 0.1
theta1=c(-0.5,0.5 ,1,0.75,0.75);
theta2=c(0,0.25,1,0,0.2);
theta3=c(1.0 ,0.5 ,1,0.5,0.2);
theta4=c(0 ,0.75,1,1,1);
pars =rbind(theta1,theta2,theta3,theta4)
# Define the jump diffusion using the DiffusionRjgqd syntax:
G0=function(t){theta[1]}
Q0=function(t){theta[2]^2}
# State dependent intensity:
Lam0 = function(t){theta[3]}
Ja = function(t){theta[4]}
Jb = function(t){theta[5]}
sim=function(tau,theta)
{
N= 100000
X=rep(Xs,N)
dt=1/2000
tt = seq(0,tau,dt)
for(i in 2:length(tt))
{
X = X +theta[1]*dt+theta[2]*rnorm(N,0,sqrt(dt))
events = (1-exp(-theta[3]*dt)>runif(N))
if(any(events))
{
wh=which(events)
X[wh]=X[wh]+rdoublex(length(wh),theta[4],theta[5])
}
}
return(list(X=X))
}
herm.density=function(x,y,t,theta,type =1)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
mu.z = theta[4]
b.z = theta[5]
if(type == 1)
{
dens=1/sqrt(2*pi*sig^2*t)*exp(-(y-x)^2/(2*sig^2*t)+mu/sig^2*(y-x))*(1-(mu^2/(2*sig^2)+lam)*t)+lam*t*ddoublex(y-x,mu.z,b.z)
}
if(type ==2)
{
dens=dnorm(y,x+mu*t,sqrt(sig^2*t))*(1-lam*t)+lam*t*ddoublex(y-x,mu.z,b.z)
}
return(list(density=dens,Xt=y))
}
conv.density=function(x,y,t,theta)
{
mu = theta[1]
sig = theta[2]
lam = theta[3]
mu.z = theta[4]
b.z = theta[5]
term1 = dnorm(y,x+mu*t,sqrt(sig^2*t))*(1-lam*t)
term2 = lam*t/(2*b.z)*exp((y-x-mu*t-b.z^-1*sig^2*t)^2/(2*sig^2*t)-(y-x-mu*t)^2/(2*sig^2*t)+mu.z/b.z)*(1-pnorm(mu.z,(y-x-mu*t-b.z^2*sig^2*t),sqrt(sig^2*t),lower.tail=T))
term3 = lam*t/(2*b.z)*exp((y-x-mu*t+b.z^-1*sig^2*t)^2/(2*sig^2*t)-(y-x-mu*t)^2/(2*sig^2*t)-mu.z/b.z)*pnorm(mu.z,(y-x-mu*t+b.z^2*sig^2*t),sqrt(sig^2*t),lower.tail=T)
return(list(density=term1+term2+term3,Xt=y))
}
for(i in 1:4)
{
theta=pars[i,]
res_1 = JGQD.density(Xs,seq(0,3,1/100),0,1,1/100,factorize=T,Jdist='Laplace',
Dtype='Normal.A')
res_2 = herm.density(Xs,res_1$Xt,tau,theta,type=1)
res_3 = conv.density(Xs,res_1$Xt,tau,theta)
res.sim = sim(tau,theta)
expr1 = expression(f[SPT]^{(4)})
expr2 = expression(f[ser])
expr3 = expression(tilde(f))
eprs =c(expr1,expr2,expr3)
ltys = c('solid','solid','dashed')
cols = c('blue','#5AAE61','black')
lwds = c(1,1,1)
pchs= c(1,3,NA)
thin=seq(1,length(res_1$Xt),20)
hist(res.sim$X,freq=FALSE,border='white',main=paste0('Brownian motion / Laplace jumps (t = ',tau,')'),
breaks=seq(-15,15,1/20),xlab=expression(X[t]),ylab='Density',las=1,xlim=range(res_1$Xt),ylim=c(0,max(res_1$density[,tau*100])*1.1),col='gray74')
lines((res_1$density[,tau*100])~res_1$Xt,type='l',col=cols[1],lty=ltys[1],lwd=lwds[1])
points((res_1$density[thin,tau*100])~res_1$Xt[thin],col=cols[1],pch=pchs[1],cex=1,lwd=lwds[1])
lines((res_2$density)~res_2$Xt,col=cols[2],lty=ltys[2],lwd=lwds[2])
points((res_2$density[thin])~res_2$Xt[thin],col=cols[2],pch=pchs[2],cex=1,lwd=lwds[2])
lines((res_3$density)~res_3$Xt,col=cols[3],lty=ltys[3],lwd=lwds[3])
legend('topright',legend=eprs,lty=ltys,lwd=lwds,col=cols,pch=pchs,bty='n',cex=1.25,pt.cex=1)
box()
}##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : theta[4]
## Jb : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : theta[4]
## Jb : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : theta[4]
## Jb : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : theta[1]
## G1
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : theta[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : theta[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : theta[4]
## Jb : theta[5]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
State dependent models: Characteristic equations (Section 4.2)
Note: the characteristic_function_to_density () was adapted from an excellent post by Vincent Zoonekynd.
rm(list=ls(all=TRUE))
#=========================================================
# Part 1: BAJD
# Evaluate the transition density at time tau =0.25 using
# the characteristic function
#=========================================================
library(DiffusionRjgqd)
# Some peripheral parameters:
theta1 = c(1 ,5 ,0.5 ,3 ,0.5 ); xlim1 =c(3,6)
theta2 = c(0.5 ,5 ,0.25 ,5 ,1 ) ; xlim2 =c(3,6)
theta3 = c(1 ,5 ,0.5 ,4 ,1.5) ; xlim3 =c(3,6)
theta4 = c(0.75,5 ,0.25 ,6 ,3 ) ; xlim4 =c(2,9)
pars = rbind(theta1,theta2,theta3,theta4)
xs = rbind(xlim1,xlim2,xlim3,xlim4)
taus = c(0.1,0.2,0.25,5)
x0 = 4
sim=function(tau,theta)
{
N= 25000
X=rep(x0,N)
dt=1/2000
tt = seq(0,tau,dt)
for(i in 2:length(tt))
{
X = X +theta[1]*(theta[2]-X)*dt+theta[3]*sqrt(X)*rnorm(N,0,sqrt(dt))
events = (1-exp(-theta[5]*dt)>runif(N))
if(any(events))
{
wh=which(events)
X[wh]=X[wh]+rexp(length(wh),theta[4])
}
}
return(list(X=X))
}
# From CF to the appr. density:
characteristic_function_to_density = function(phi,n,a,b)
{
i = 0:(n-1)
dx = (b-a)/n
x = seq(a,b-dx,length =n)
dt = 2*pi / ( n * dx )
c = -n/2 * dt
d = n/2 * dt
t = c + i * dt
phi_t = phi(t)
X = exp(-(0+1i)*i*dt*a)*phi_t
Y = fft(X)
density = dt /(2*pi)*exp(-(0+1i)*c*x)*Y
list(x = x,density = Re(density))
}
# The characteristic function:
f =function(t,theta =th)
{
a = theta[1]
b = theta[2]
sigma = theta[3]
nu = theta[4]
cc = theta[5]
u=(0+1i)*t
term1 = (1-(sigma^2*u)/(2*a)*(1-exp(-a*tau)))^(-2*a*b/sigma^2)
term2 = ((nu*(1-sigma^2*u/(2*a))-(1-sigma^2*nu/(2*a))*u*exp(-a*tau))/(nu-u))^(cc/(a-sigma^2*nu/2))
term3 = exp(x0*u*exp(-a*tau)/(1-sigma^2/(2*a)*u*(1-exp(-a*tau))))
rs=term1*term2*term3
return(rs)
}
JGQD.remove()## [1] "Removed : NA "G0=function(t){th[1]*th[2]}
G1=function(t){-th[1]}
Q1=function(t){th[3]^2}
Lam0 =function(t){th[5]}
Jlam =function(t){1/th[4]}
for(i in 1:4)
{
th = pars[i,]
tau =taus[i]
#res.sim = sim(tau,th)
d <-characteristic_function_to_density(f,2^10,1,10)
res_1=JGQD.density(4,d$x,0,tau,delt=1/100,Jdist='Exponential',factorize=TRUE)#,Dtype='Gamma.Saddle')
res_2=JGQD.density(4,d$x,0,tau,delt=1/100,Jdist='Exponential',factorize=FALSE)#,Dtype='Gamma.Saddle')
xlims= xs[i,]
ylims= c(0,2.8)
cols= c('blue','#BBCCEE','black')
ltys= c('solid','dashed','solid')
lwds = c(1,2,1)
pchs= c(1,NA,NA)
thin=seq(1,length(res_1$Xt),25)
eprs =c('Saddlepoint (Factorized)',
'Saddlepoint (Unfactorized)',
'Fourier')
#hist(res.sim$X,breaks=seq(0,25,1/10),freq=F,border='white',col='white',
plot(d$x,(d$density),type='n',axes=T ,col=cols[1],lty=ltys[1],lwd=lwds[1],
xlim=xlims, ylab='Density',xlab=expression(X[t]),
main=paste0('BAJD transition density (t = ',tau,')'))
lines(d$density~d$x ,col=cols[3],lty=ltys[3],lwd=lwds[3])
lines(res_1$density[,100*tau]~res_1$Xt ,col=cols[1],lty=ltys[1],lwd=lwds[1])
points(res_1$density[thin,100*tau]~res_1$Xt[thin],col=cols[1],pch=pchs[1])
lines (res_2$density[,100*tau]~res_2$Xt,col=cols[2],lty=ltys[2],lwd=lwds[2])
points(res_2$density[thin,100*tau]~res_2$Xt[thin],col=cols[2],pch=pchs[2])
legend('topright',legend=eprs,lty=ltys,lwd=lwds,col=cols,pch=pchs,bty='n')
}##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : th[1]*th[2]
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0
## Q1 : th[3]^2
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[5]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Exponential
## Jlam : 1/th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
rm(list=ls(all=TRUE))
par(mfrow=c(1,1))
#=========================================================
# Part 2: OU - type jump diffusion with Laplace jumps
# Evaluate the transition density at time tau =0.15 using
# the characteristic function
#=========================================================
library(DiffusionRjgqd)
# Some peripheral parameters:
theta1 = c(1 ,1.2 ,0.5,0.5) ;xlim1 =c(-3,2)
theta2 = c(0.5 ,0.5 ,0.2,0.1 ) ;xlim2 =c(-2,0.5)
theta3 = c(1 ,0.25 ,1 ,0.1 ) ;xlim3 =c(-1.5,0)
theta4 = c(0.75,0.5 ,1 ,0.2) ;xlim4 =c(-2,2)
pars = rbind(theta1,theta2,theta3,theta4)
xs = rbind(xlim1,xlim2,xlim3,xlim4)
taus = c(0.1,0.2,0.25,5)
x0 = -1
# From CF to the appr. density:
characteristic_function_to_density = function(phi,n,a,b)
{
i = 0:(n-1)
dx = (b-a)/n
x = seq(a,b-dx,length =n)
dt = 2*pi / ( n * dx )
c = -n/2 * dt
d = n/2 * dt
t = c + i * dt
phi_t = phi(t)
X = exp(-(0+1i)*i*dt*a)*phi_t
Y = fft(X)
density = dt /(2*pi)*exp(-(0+1i)*c*x)*Y
list(x = x,density = Re(density))
}
# The characteristic function:
f =function(u,theta =th)
{
a = theta[1]
sigma = theta[2]
lam = theta[3]
sigb = theta[4]
u=(0+1i)*u
term1 = exp(u*x0*exp(-a*tau)-u^2*sigma^2/(4*a)*(exp(-2*a*tau)-1))
term2 = ((1-exp(-2*a*tau)*u^2*sigb^2)/(1-u^2*sigb^2))^{lam/(2*a)}
rs=term1*term2
return(rs)
}
JGQD.remove()## [1] "Removed : NA "G0=function(t){0}
G1=function(t){-th[1]}
Q0=function(t){th[2]^2}
Lam0=function(t){th[3]}
Ja =function(t){0}
Jb =function(t){th[4]}
for(i in 1:4)
{
tau=taus[i]
th = pars[i,]
xlims = xs[i,]
d <-characteristic_function_to_density(f,2^10,-5,5)
res_1=JGQD.density(x0,d$x,0,tau,delt=1/100,Jdist='Laplace',factorize=TRUE)
res_2=JGQD.density(x0,d$x,0,tau,delt=1/100,Jdist='Laplace',factorize=FALSE)
ylims= c(0,3.1)
cols= c('blue','#BBCCEE','black')
ltys= c('solid','dashed','solid')
lwds = c(1,2,2)
pchs= c(1,NA,NA)
thin=seq(1,length(res_1$Xt),25)
eprs =c(
'Saddlepoint (Factorized)',
'Saddlepoint (Unfactorized)',
'Fourier')
plot(d$x,(d$density),type='l',axes=T ,col=cols[3],lty=ltys[3],lwd=lwds[3],
xlim=xlims, ylab='Density',xlab=expression(X[t]),
main=paste0('OU transition density (t = ',tau,')'))
lines(res_1$density[,100*tau]~res_1$Xt ,col=cols[1],lty=ltys[1],lwd=lwds[1])
points(res_1$density[thin,100*tau]~res_1$Xt[thin],col=cols[1],pch=pchs[1])
lines (res_2$density[,100*tau]~res_2$Xt,col=cols[2],lty=ltys[2],lwd=lwds[2])
points(res_2$density[thin,100*tau]~res_2$Xt[thin],col=cols[2],pch=pchs[2])
legend('topright',legend=eprs,lty=ltys,lwd=lwds,col=cols,pch=pchs,bty='n')
}##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Normal.A
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
##
## ================================================================
## Jump Generalized Quadratic Diffusion (JGQD)
## ================================================================
## _____________________ Drift Coefficients _______________________
## G0 : 0
## G1 : -th[1]
## G2
## ___________________ Diffusion Coefficients _____________________
## Q0 : th[2]^2
## Q1
## Q2
## _______________________ Jump Mechanism _________________________
## ......................... Intensity ............................
## Lam0 : th[3]
## Lam1
## Lam2
## ........................... Jumps ..............................
## Laplace
## Ja : 0
## Jb : th[4]
## __________________ Distribution Approximant ____________________
## Density approx. : Saddlepoint
## Trunc. Order : 8
## Dens. Order : 4
## =================================================================
Likelihood inference
Exact vs. approximate likelihoods (Section 5.1)
rm(list=ls(all=TRUE))
#=========================================================
# Inference: BM with Normal jumps
#=========================================================
set.seed(1)
library(DiffusionRjgqd)
thetas = c(1 ,1 ,1 ,1 ,1)*1
sim=function(x0,theta)
{
delt = 1/1000
tt= seq(0,25,delt)
N = length(tt)
M = 300
X = matrix(x0,M,N)
for(i in 2:N)
{
X[,i] = X[,i-1] +theta[1]*delt+theta[2]*rnorm(M,0,sqrt(delt))
events = (theta[3]*delt>runif(M))
if(any(events))
{
wh =which(events)
zz=rnorm(length(wh),theta[4],theta[5])
X[wh,i]=X[wh,i]+zz
}
}
ttt= seq(1,N,by=100)
return(list(Xt=X[,ttt],time=tt[ttt]))
}
x=sim(0,thetas)
JGQD.remove()## [1] "Removed : NA "G0=function(t){theta[1]}
Q0=function(t){theta[2]*theta[2]}
Lam0 =function(t){theta[3]}
Jmu =function(t){theta[4]}
Jsig =function(t){theta[5]}
updates =1000
burns =500
theta1 = c(1 ,5 ,0.25 ,3 ,1)+2
sds = c(0.1,0.05,0.02,0.2,0.1)/10
res2 = JGQD.mcmc(x$Xt[1,],x$time,20,theta1,sds,updates,burns,Jdist='Normal',Jtype='Add',print.output = FALSE, plot.chain = FALSE)## =================================================================par(mfrow=c(1,1))
X=x$Xt
N=length(X[1,])
Xt = X[1,-1]
Xs = X[1,-N]
tau = diff(x$time)[1]
mesh2=10
delt=tau/mesh2
like=function(theta)
{
x = Xs; y = Xt; t = tau; order = 10
mu = theta[1]
sig = theta[2]
lam = theta[3]
mu2 = theta[4]
sig2 = theta[5]
dens=exp(-lam*t)*dnorm(y,x+mu*t,sqrt(sig^2*t))
for(ii in 1:order)
{
dens= dens +exp(-lam*t)*(lam*t)^ii/factorial(ii)*dnorm(y,x+mu*t+ii*mu2,sqrt(sig^2*t+ii*sig2^2))
}
return(-sum(log(dens)))
}
like2 = function(theta)
{
rs=solver(Xs,Xt,c(0,theta),mesh2,delt,N-1,x$time[-N],100,0.1,1,10,4)
return(-sum(rs$like))
}
checkf=function(x){any(round(x,2)==0)||any(round(x,2)>=10)}
Nsims = 250
plot(x$Xt[1,]~x$time,type='n',ylim=range(x$Xt),main='Simulated trajectories',xlab='Time (t)',ylab=expression(X[t]))
for(i in 1:200)
{
points(x$Xt[i,]~x$time,pch='.',col='gray54',cex=1)
}
legend('bottomright',legend='Observations',pch='.',col='gray54',bty='n',cex=1,pt.cex=5)
pars0 = matrix(0,5,Nsims)
likes0 = rep(0, Nsims)
pars1 = matrix(0,5,Nsims)
likes1 = rep(0, Nsims)
k=1
j=1
while(k <=Nsims)
{
N=length(X[1,])
Xt = X[j,-1]
Xs = X[j,-N]
tau = diff(x$time)[1]
res0= nlminb(rep(1,5),like , lower = rep(0,5), upper = rep(10,5))
res1= nlminb(rep(1,5),like2, lower = rep(0,5), upper = rep(10,5))
check =!(checkf(res0$par)||checkf(res1$par))
if(check)
{
pars0[,k] = res0$par
likes0[k] = res0$objective
pars1[,k] = res1$par
likes1[k] = res1$objective
k=k+1
}
j=j+1
}
leave.out = c(1)
LIMS=c(-0.7,0.7)*2
expr=expression(hat(mu)[x] - mu[x]^{true},
hat(sigma)[x]- sigma[x]^{true},
hat(lambda) - lambda^{true},
hat(mu)[z] - mu[z]^{true},
hat(sigma)[z]- sigma[z]^{true})
expr2=expression(L[true](theta,D[S] ))
expr3=expression(L[GQD](theta,D[S]))
dt=c(1/10,1/25,1/10,1/10,1/10)*2
for(i in 1:5)
{
h0 =hist(pars1[i,-leave.out]-thetas[i],plot=F,breaks=seq(-3,3,dt[i]))
h1 =hist(pars1[i,-leave.out]-thetas[i],plot=F,breaks=seq(-3,3,dt[i]))
hist(pars0[i,-leave.out]-thetas[i],freq=F,border='white',col='gray84',lwd=1,xlab=expr[i],main=paste0('Deviation of MLE from true parameter'),ylim=c(0,max(c(h0$density,h1$density))),xlim=LIMS,breaks=seq(-3,3,dt[i]))#breaks=seq(-2,2,dt[i])#xlim=range(pars0[i,-leave.out]-thetas[i]))
hist(pars1[i,-leave.out]-thetas[i],add=T,breaks=seq(-3,3,dt[i]),col=rgb(70/255,130/255,180/255,0.25),freq=F,border='white')
legend('topright',bty='n',col=c('gray74',rgb(70/255,130/255,180/255,0.25)),legend=c('True likelihood','Approx. likelihood'),lty=c('solid','solid'),lwd=c(5,5))
abline(v=0,lty='dashed')
box()
}
h0= hist(likes0[-leave.out],freq=F,border='white',col='gray84',main='Likelihood at MLE',xlab='log-likelihood at max.')
h1 =hist(likes1[-leave.out],add=T,col=rgb(70/255,130/255,180/255,0.25),freq=F,border='white')
legend('topright',bty='n',col=c('gray74',rgb(70/255,130/255,180/255,0.25)),legend=c('True likelihood','Approx. likelihood'),lty=c('solid','solid'),lwd=c(5,5))
box()
Inference for non-linear processes (Section 5.2)
rm(list=ls(all=TRUE))
#=========================================================
# Inference: CIR Process With Stochastic Intensity
#=========================================================
set.seed(200707881)
library(DiffusionRjgqd)
library("smoothmest")
thetas = c(1 ,5 ,0.1 ,0.2, 0.5 ,0.5)
sim=function(x0,theta)
{
delt = 1/1000
tt= seq(0,25,delt)
N = length(tt)
X = rep(x0,N)
J = c()
Times = c()
StatesA = c()
StatesB = c()
k=0
dd = 0
for(i in 2:N)
{
X[i] = X[i-1] +theta[1]*(theta[2]-X[i-1])*delt+theta[3]*sqrt(X[i-1])*rnorm(1,0,sqrt(delt))
events = (theta[4]*X[i-1]*delt>runif(1))
if(events)
{
k=k+1
zz=rnorm(1,theta[5],theta[6])
J[k] = zz
Times[k] = tt[i]
StatesA[k] = X[i]
X[i]=X[i]+zz
StatesB[k] = X[i]
}
}
ttt= seq(1,N,by=100)
return(list(Xt=X[ttt],time=tt[ttt],Times=Times,Jumps = J,StatesA = StatesA, StatesB = StatesB,n = length(J)))
}
x=sim(5,thetas)
extract.diff = function(x)
{
n1 = length(x$Xt)
n2 = x$n
X= rep(0,n1+2*n2)
compx = c(x$Xt,x$StatesA)
compt = c(x$time,x$Times)
od=order(compt)
time=compt[od]
X = compx[od]
return(list(X=X,time=time))
}
res = extract.diff(x)
N. =length(res$X)
Xs. = res$X[-N.]
Xt. = res$X[-1]
d = diff(res$time)
skips = which(round(d,3)!=0.100)
d[skips] = 0.1
CIR.likelihood=function(theta)
{
kappa= theta[1]
beta = theta[2]
sig = theta[3]
cc = 2*kappa/((sig^2)*(1-exp(-kappa*d)))
u = cc*Xs.*exp(-kappa*d)
v = cc*Xt.
qqq = 2*kappa*beta/(sig^2)-1
like=(log(cc)+(-u-v)+(qqq/2)*log(v/u) +(log(besselI(2*sqrt(u*v),qqq,expon.scaled=T))+2*sqrt(u*v)))
return(-sum(like[-skips]))
}
Jump.likelihood = function(th){-sum(dnorm(x$Jumps,th[1],th[2],log=TRUE))}
mm =function(X,dtau,a = 1,b = 5){(X-b)*(1-exp(-a*dtau))/a+b*dtau}
Pois.likelihood = function(th)
{
dd = -th[1]*mm(x$StatesB[-x$n],diff(x$Times))+log(th[1])
-sum((dd))
}
JGQD.remove()
G0=function(t){theta[1]*theta[2]}
G1=function(t){-theta[1]}
Q1=function(t){theta[3]*theta[3]}
Lam1 =function(t){theta[4]}
Jmu =function(t){theta[5]}
Jsig =function(t){theta[6]}
updates =50
burns =5
theta1 = c(1 ,10 ,1 ,1 ,1,1)*0.5
sds = c(0.1,0.05,0.02,0.2,0.1,0.1)/20
res2 = JGQD.mcmc(x$Xt,x$time,20,theta1,sds,updates,burns,Jdist='Normal',Jtype='Add',plot.chain = FALSE, print.output = FALSE)
X=x$Xt
N=length(x$X)
Xt = X[-1]
Xs = X[-N]
tau = diff(x$time)[1]
mesh2=10
delt=tau/mesh2
JD.likelihood = function(theta)
{
rs=solver(Xs,Xt,c(0,theta),mesh2,delt,N-1,x$time[-N],100,0.1,1,10,4)
return(-sum((rs$like)))
}
checkf=function(x){any(round(x,2)==0.05)||any(round(x,2)>=10)}
Nsims = 5
pars0 = matrix(0,6,Nsims)
pars1 = matrix(0,6,Nsims)
likes1 = rep(0, Nsims)
k=1
j=1
while(k <=Nsims)
{
x=sim(5,thetas)
res = extract.diff(x)
N. = length(res$X)
Xs. = res$X[-N.]
Xt. = res$X[-1]
d = diff(res$time)
skips = which(round(d,3)!=0.100)
d[skips] = 0.5
Diff.res = nlminb(runif(3),CIR.likelihood ,lower=rep(0.05,3),upper=rep(10,3))
X=x$Xt
Xt = X[-1]
Xs = X[-N]
tau = diff(x$time)[1]
mesh2=10
delt=tau/mesh2
JD.res = nlminb(rep(1,6),JD.likelihood ,lower=rep(0.05,6),upper=rep(10,6))
Jump.res = c(mean(x$Jumps),sd(x$Jumps))
Pois.res = nlminb(runif(1),Pois.likelihood,lower=rep(0.05,1),upper=rep(10,1))
check =!(checkf(JD.res$par))
print(k)
if(check)
{
pars1[,k] = JD.res$par
likes1[k] = JD.res$objective
k=k+1
pars0[,k] = c(Diff.res$par, Pois.res$par,Jump.res)
}
j=j+1
}
leave.out = c(1)
LIMS=c(-1,1)*0.5
expr=expression(hat(alpha)[x] - alpha[x]^{true},
hat(beta)[x] - beta[x]^{true},
hat(sigma)[x] - sigma[x]^{true},
hat(lambda) - lambda^{true},
hat(mu)[z]- mu[z]^{true},
hat(sigma)[z]- sigma[z]^{true})
expr2=expression(L[true](theta,D[S] ))
expr3=expression(L[GQD](theta,D[S]))
par(mfrow=c(1,1))
for(i in 1:6)
{
ddd= density((pars0[i,-leave.out]-thetas[i])/1,bw ='SJ')
val = (pars1[i,-leave.out]-thetas[i])/1
h1 =hist(val,plot=F)
hist(val,freq=F,border='white',col='gray84',lwd=1,xlab=expr[i],main=paste0('Deviation of MLE from true parameter'),ylim=c(0,max(c(ddd$y,h1$density)))
)
h2 =hist((pars0[i,-leave.out]-thetas[i]),add=T,col=rgb(70/255,130/255,180/255,0.25),freq=F,border='white',breaks=seq(-10,10,diff(h1$breaks)[1]))
legend('topright',bty='n',col=c('gray74',rgb(70/255,130/255,180/255,0.25)),legend=c('Jump diff. est.','Direct est.'),lty=c('solid','solid'),lwd=c(5,5))
abline(v=0,lty='dashed')
box()
}Estimate the parameters of a simulated bivariate process (Section 5.3)
data(JSDEsim2)
attach(JSDEsim2)
data(JSDEsim3)
attach(JSDEsim3)
#------------------------------------------------------------------------------
# Define parameterized coefficients of the process, and set up starting
# parameters.
# True model: dX_t = 0.5(2+Y_t-X_t)dt+0.1sqrt{X_tY_t}dB_t +dP_t^1
# dX_t = 1(5-Y_t)dt+0.1sqrt{X_t}dW_t +dP_t^2
# where dP_t^1 = z_tdN_t, dP_t^1 = z_tdN_t describes a Poisson
# process with intensity:
# lambda(X_t,Y_t) = 1
# and
# {z_1,z_2}' ~ Bivariate Normal({0.5,0.5}',diag({0.5,0.5}'))
#------------------------------------------------------------------------------
par(mfrow=c(1,1))
plot(Xt~time,type='l',col='#BBCCEE',ylim=c(-3,13),xlim=c(0,60),axes=F,main='Simulated Trajectory',xlab = 'Time',ylab ='X_t')
lines(Yt~time,type='l',col='#222299')
axis(1,at=seq(0,50,5))
axis(1,at=seq(0,50,5/5),tcl=-0.2,labels=NA)
axis(2,at=seq(-3,13,1))
axis(2,at=seq(-3,13,1/5),tcl=-0.2,labels=NA)
lines(Xjumps[Xjumps!=0]~Jtime[Xjumps!=0],type='h',col='#BBCCEE')
lines(Yjumps[Yjumps!=0]~Jtime[Yjumps!=0],type='h',col='#222299')
mx=mean(Xjumps[Xjumps!=0])
sx=sd(Xjumps[Xjumps!=0])
my=mean(Yjumps[Yjumps!=0])
sy=sd(Yjumps[Yjumps!=0])
legend('topright',lty=c(1,2),col=c('#BBCCEE','#222299'),legend=c(expr1,expr2),cex=0.9)
segments(50,-5,50,9,lty='dotted')
xx=seq(-3,3,1/10)
yy=dnorm(xx,mx,sx)
yy = (yy-min(yy))/(max(yy)-min(yy))*9+51*1
lines(xx~yy,col='#BBCCEE')
yy=dnorm(xx,my,sy)
yy = (yy-min(yy))/(max(yy)-min(yy))*9+51*1
lines(xx~yy,col='#222299')
text(55,10.0,substitute(hat(theta)[8]==a,list(a=round(mx,2))),cex=0.8)
text(55,9.0,substitute(hat(theta)[9]==a,list(a=round(sx,2))),cex=0.8)
text(55,8.0,substitute(hat(theta)[10]==a,list(a=round(my,2))),cex=0.8)
text(55,7.0,substitute(hat(theta)[11]==a,list(a=round(sy,2))),cex=0.8)
abline(h=0,lty='dotted')
X=cbind(Xt,Yt)
# Define the model:
JGQD.remove()
a00 <-function(t){theta[1]*theta[2]}
a10 <-function(t){-theta[1]}
a01 <-function(t){theta[1]}
c11 <-function(t){theta[3]*theta[3]}
b00 <-function(t){theta[4]*theta[5]}
b01 <-function(t){-theta[4]}
f01 <-function(t){theta[6]*theta[6]}
# Constant intensity
Lam00= function(t){theta[7]}
# Normal jumps:
Jmu1 <-function(t){theta[8]}
Jmu2 <-function(t){theta[9]}
Jsig11 <-function(t){theta[10]*theta[10]}
Jsig22 <-function(t){theta[11]*theta[11]}
# Some starting parameters:
theta <-c(rep(1,11))
sds <-c(0.08,0.22,0.01,0.04,0.16,0.01,0.10,0.07,0.09,0.05,0.09)/2
burns <-10000
updates<-50000
res <-BiJGQD.mcmc(X,time,mesh=10,theta,sds,updates,burns=burns)
JGQD.estimates(res,thin=200,burns)[,1:3]Jump diffusion models for Google equity volatility (VXGOG) (Section 6)
See also:
rm(list=ls(all=TRUE))
library(DiffusionRgqd)
library(DiffusionRjgqd)
library(Quandl)
library(MCMCpack)
library(coda)
# Source data for the Google equity volatility:
quandldata1 <- Quandl("CBOE/VXGOG", collapse="daily",
start_date="2010-03-11",end_date="2016-01-01", type="raw")
Vt <- rev(quandldata1[,names(quandldata1)=='Close'])
time1 <-rev(quandldata1[,names(quandldata1)=='Date'])
# Make a plot of the data
plot(Vt~time1,type='l',col='#222299',ylim=c(10,55),main='Google Equity Volatility (VXGOG)',xlab = 'Time',ylab ='Volatility %',lwd=1)
X <- Vt
time <- cumsum(c(0,diff(as.Date(time1))*(1/365)))
updates = 110000
burns = 10000
#===========================================================================
# Jump free models first
#===========================================================================
pfnc = function(theta)
{
dinvgamma(theta[2], 0.001, 0.001)*
dnorm(theta[2], 25, 5)*
dinvgamma(theta[3]^2, 0.001, 0.001)
}
GQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q0 = function(t){theta[3]*theta[3]}
priors=function(theta){pfnc(theta)}
theta = c(100, 50, sqrt(30))
sds = c(2.29, 0.88, 0.72)
model_1 = GQD.mcmc(X,time,10,theta,sds,updates,burns)
GQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q1 = function(t){theta[3]*theta[3]}
priors=function(theta){pfnc(theta)}
theta = c(100, 50, sqrt(30))
sds = c(2.63, 0.96, 0.14)
model_2 = GQD.mcmc(X,time,10,theta,sds,updates,burns)
GQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
priors=function(theta){pfnc(theta)}
theta = c(100, 50, sqrt(30))
sds = c(2.32, 0.86, 0.02)
model_3 = GQD.mcmc(X,time,10,theta,sds,updates,burns)
#===========================================================================
# Jump models
#===========================================================================
pfnc = function(theta)
{
dinvgamma(theta[2] , 0.001, 0.001)*
dnorm(theta[2], 25, 5)*
dinvgamma(theta[3]^2, 0.001, 0.001)*
dinvgamma(theta[4] , 0.001, 0.001)
}
JGQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q1 = function(t){theta[3]*theta[3]}
Lam0 = function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig = function(t){theta[6]}
priors =function(theta){pfnc(theta)}
theta = c(50,25,sqrt(50),100,5,20)
sds = c(1.60,0.90,0.13,3.40,0.46,0.46)/1.5
model_4 = JGQD.mcmc(X,time,10,theta,sds,updates,burns)
JGQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam0 = function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig = function(t){theta[6]}
priors = function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),100,5,20)
sds = c(1.53,0.92,0.02,3.04,0.59,0.66)/1.5
model_5 = JGQD.mcmc(X,time,10,theta,sds,updates,burns)
JGQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam1 = function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig = function(t){theta[6]}
priors = function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),4,5,20)
sds = c(1.53,0.92,0.02,0.34,0.59,0.66)/1.5
model_6 = JGQD.mcmc(X,time,10,theta,sds,updates,burns)
JGQD.remove()
G0 = function(t){theta[1]*theta[2]}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam0= function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig= function(t){theta[6]}
priors=function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),100,0.5,0.5)
sds = c(1.51,1.69,0.02,3.04,0.01,0.01)/1.5
model_7 = JGQD.mcmc(X,time,10,theta,sds,updates,burns, Jtype='Mult')
pfnc = function(theta)
{
rs =dgamma(theta[2] , 0.001, 0.001)*
dnorm(theta[2], 25, 5)*
dgamma(1/theta[3]^2, 0.001, 0.001)*
dgamma(theta[4], 0.001, 0.001)*
dgamma(theta[7], 0.001, 0.001)*
dbeta(theta[8],0.5,0.5)
if(is.na(rs)){return(0.000000000001)}
return(rs)
}
JGQD.remove()
G0 = function(t){theta[1]*theta[2]+theta[1]*theta[7]*sin(8*pi*t+(theta[8]-0.5)*2*pi)}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam0= function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig= function(t){theta[6]}
priors=function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),50,5,10,20,0.5)
sds = c(1.45,0.84,0.02,3.72,0.57,0.53,0.87,0.03)/1.5
model_8 = JGQD.mcmc(X,time,10,theta,sds,updates,burns,decode=F)
JGQD.remove()
G0 = function(t){theta[1]*theta[2]+theta[1]*theta[7]*sin(8*pi*t+(theta[8]-0.5)*2*pi)}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam1= function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig= function(t){theta[6]}
priors=function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),2,20,20,20,0.5)
sds = c(1.87,0.96,0.02,0.17,0.39,0.40,0.71,0.03)/1.5
model_9 = JGQD.mcmc(X,time,5,theta,sds,updates,burns)
JGQD.remove()
G0 = function(t){theta[1]*theta[2]+theta[1]*theta[7]*sin(8*pi*t+(theta[8]-0.5)*2*pi)}
G1 = function(t){-theta[1]}
Q2 = function(t){theta[3]*theta[3]}
Lam0= function(t){theta[4]}
Jmu = function(t){theta[5]}
Jsig = function(t){theta[6]}
priors=function(theta){pfnc(theta)}
theta = c(50,25,sqrt(9),100,0.1,0.1,50,0.5)
sds = c(1.65,0.79,0.02,4.29,0.01,0.01,0.94,0.03)/1.5
model_10 = JGQD.mcmc(X,time,10,theta,sds,updates,burns, Jtype='Mult')
S = list(model_1, model_2, model_3, model_4, model_5, model_6, model_7, model_8, model_9,model_10)
JGQD.dic(S)Further reading
- Methodology paper and replication materials.
- Related packages:
- More Vignettes:
browseVignettes('DiffusionRjgqd')