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))
