# Generation of Levy processes

## Compound Poisson process

## Example: Compound Poisson jumps
set.seed(234)
# Compound Poisson process is characterized by intensiry and jump law.
mod.cpp <- setPoisson(intensity = "lambda", df = list("dnorm(z,mu,sigma)"))

T <- 10 # terminal sampling time
samp <- setSampling(Initial = 0, Terminal = T, n = 100*T)
cpp <- setYuima(model = mod.cpp, sampling = samp)

param <- list(lambda = 10, mu = 0, sigma = 2)

# str(cpp) # to see the structure of cpp.

X.cpp <- simulate(cpp, xinit = 0, true.parameter = param)
plot(X.cpp, col = "blue", main="Compound Poisson process")
grid()

# dev.copy2pdf(file="cpp.pdf", width=8, height=5) # figure output (scale: inch)

## Inverse Gaussian subordinator

## Example: IG-Levy process (subordinator)
set.seed(234)
T <- 1
n <- 1000 # "n" may be set here.
jump <- "1" # the jump coefficient, here unity.

mod.ig <- setModel(drift = "0", diffusion="0", jump.coeff = jump, solve.variable = "x",
state.variable = "x", measure.type = "code", measure = list(df="rIG(z,3,3)"))

samp <- setSampling(Terminal=T, n=n)
ig <- setYuima(model=mod.ig, sampling=samp)

# str(ig)

X.ig <- simulate(ig, xinit=0)
plot(X.ig, col="blue", main="IG-Levy process (Subordinator)")
grid()

# dev.copy2pdf(file="ig.pdf", width=8, height=5) ## figure output (scale: inch)

## Normal inverse Gaussian (NIG) Levy process

## Example: NIG-Levy process
set.seed(234)
T <- 1; n <- 1000; jump <- "1"

mod.nig <- setModel(drift = "0", diffusion="0", jump.coeff = jump,
solve.variable = "x", state.variable = "x",
measure.type = "code", measure = list(df="rNIG(z,3,0,3,0)"))

samp <- setSampling(Terminal=T, n=n)
nig <- setYuima(model=mod.nig, sampling=samp)

# str(nig)

X.nig <- simulate(nig, xinit=0)
plot(X.nig, col="blue", main="NIG-Levy process")
grid()

# dev.copy2pdf(file="nig.pdf", width=8, height=5) # scale: inch

## Checking by histogram:
n <- 100000
X.nig <- rNIG(n,5,-0.5,5,0)           # NIG(5,-0.5,5,0)-rng
pdf <- function(x) dNIG(x,5,-0.5,5,0) # NIG(5,-0.5,5,0)-pdf
hist(X.nig, prob=T, breaks=50, col="lightblue",
xlim=c(-5,4), ylim=c(0,0.45), xlab="", ylab="Empirical density",
main="NIG-rng histgram checking")
curve(pdf, add = TRUE, col="red", lwd = 2)  # target density plot
grid()

# dev.copy2pdf(file="nig-hist.pdf", width=8, height=5) ## figure output (scale: inch)

# Generation of Levy driven SDE

## Diffusion with compound Poisson jumps

## Example: Diffusion with compound Poisson jumps
set.seed(234)
T <- 1; n <- 1000; jump <- "-1"

mod.dcppj <- setModel(drift=c("sin(x)-theta*x"), diffusion="sigma", jump.coeff=jump,
measure=list(intensity="lambda",df=list("dgamma(z,3,3)")),
measure.type="CP",solve.variable="x")
samp <- setSampling(Terminal=T, n=n)
dcppj <- setYuima(model=mod.dcppj, sampling=samp)

param <- list(lambda = 10, theta = 1, sigma = 2)

# str(dcppj)

X.dcppj <- simulate(dcppj, xinit=0, true.parameter=param)
plot(X.dcppj, col="blue", main="Diffusion with compound Poisson jumps")
grid()

# dev.copy2pdf(file="dcppj.pdf", width=8, height=5) ## figure output (scale: inch)

## Geometric NIG Levy process

## Example: Geometric NIG-Levy process
set.seed(123)
T <- 1; n <- 1000; drift <- "0"; diff <- "0"; jump <- "x"

mod.gnig <- setModel(drift = drift, diffusion=diff, jump.coeff = jump,
solve.variable = "x", state.variable = "x",
measure.type = "code", measure = list(df="rNIG(z,5,0,5,0)"))

samp <- setSampling(Terminal=T, n=n)
gnig <- setYuima(model=mod.gnig, sampling=samp)

# str(gnig)

X.gnig <- simulate(gnig, xinit=1)
plot(X.gnig, col="blue", main="Geometric NIG-Levy process")
grid()

# dev.copy2pdf(file="gnig.pdf", width=8, height=5) ## figure output (scale: inch)

## Stable (Infinite variance) SDE

## Example: Non-linear stable-Levy SDE
set.seed(234)
T <- 1; n <- 2000; drift <- "-x/sqrt(1+x^2)"; jump <- "1"

mod.st <- setModel(drift = drift, diffusion="0", jump.coeff = jump,
solve.variable = "x", state.variable = "x",
measure.type = "code", measure = list(df="rstable(z,1.7,0,1,0)"))

samp <- setSampling(Terminal=T, n=n)
st <- setYuima(model=mod.st, sampling=samp)

# str(st)

X.st <- simulate(st, xinit=0)
plot(X.st, col="blue", main="Stable-Levy SDE")
grid()

# dev.copy2pdf(file="st.pdf", width=8, height=5) ## figure output (scale: inch)

## 2-dimensional SDE

## Example: 2-dim. NIG SDE
set.seed(345)
T <- 1; n <- 1000; x0 <- c(0,0)

# parameter values
beta <- c(0.5,-0.3); mu <- c(0,0); delta0 <- 1; alpha <- 5
Lambda <- matrix(c(1,0,0,1),2,2)

drif <- c("-2*x1","0.3*x1-1/sqrt(1+x2^2)")
jump <- matrix(c("1/sqrt(x1^2+1)","0","-0.5","1"),2,2)

mod.nigsde2 <- setModel(drift=drif, xinit=x0, jump.coeff=jump,
solve.variable=c("x1","x2"), measure.type="code",
measure=list(df="rNIG(z,alpha,beta,delta0,mu,Lambda)"))

samp <- setSampling(Terminal=T, n=n)
nigsde2 <- setYuima(model=mod.nigsde2, sampling=samp)

# str(nigsde2)

X.nigsde2 <- simulate(nigsde2, true.par=list(alpha=alpha, beta=beta,
delta0=delta0, mu=mu, Lambda=Lambda))
plot(X.nigsde2, col="blue", main="2-dimensional NIG-SDE")

# dev.copy2pdf(file="nigsde2.pdf", width=8, height=5) ## figure output (scale: inch)

# qmleLevy

## Bilateral-gamma driven SDE

\begin{aligned} dX_{t}&=-\theta_{0}X_{t}dt+\frac{\theta_{1}}{\sqrt{1+X_{t}^{2}}}dZ_{t}, \end{aligned} where the driving noise $$Z_{t}$$ obayes $$bGamma(t,\sqrt{2},t,\sqrt{2})$$. True value $$(\theta_{0,0},\theta_{1,0})=(1,2)$$.

# Model construction: drift and jump coefficients
dri<-"-theta0*x"
jum<-"theta1/(1+x^2)^(1/2)"

yuima<-setModel(drift = dri, jump.coeff = jum, solve.variable = "x",
state.variable = "x", measure.type = "code",
measure = list(df="rbgamma(z,1,sqrt(2),1,sqrt(2))"))

## set subsampling rate and terminal time, and stepsize h.
n<-100000; tp<-0.1; N<-n*tp; T<-100; hn<-T/N

## set sampling scheme
sam<-setSampling(Terminal = T, n=n)
subsam<-setSampling(Terminal = T, n=N)
yuima<-setYuima(model = yuima, sampling = sam)

set.seed(123)
true<-list(theta0 = 1, theta1 = 2) ## true values
yuima<-simulate(yuima, xinit = 0, true.parameter = true,sampling = sam,
subsampling = subsam) ## generate a path
plot(yuima, col="blue")

# dev.copy2pdf(file="qmle-bg.pdf", width=8, height=5) ## figure output (scale: inch)

## set upper and lower bounds
upper<-list(theta0 = 40, theta1 = 40)
lower<-list(theta0 = 0.5, theta1 = 1)
start<-list(theta0 = runif(1,0.5,40), theta1 = runif(1,1,40)) ## set initial values

qmleLevy(yuima, start=start, lower=lower, upper=upper, joint = TRUE)
##   theta1
## 2.000248
##
## second ## theta0 ## 0.9820445 ## 2-dimensional variance-gamma driven SDE \begin{aligned} d\begin{pmatrix}X_{1,t}\\X_{2,t}\end{pmatrix}=\begin{pmatrix}1-\theta_0X_{1,t}-X_{2,t}\\-\theta_1X_{2,t}\end{pmatrix}dt+\begin{pmatrix} \frac{\theta_2}{1+X_{1,t}^2}+1& 0\\1 & 1\end{pmatrix} dZ_t, \end{aligned} where the driving noise $$Z_t$$ obeys $$rvgamma\left(\frac{1}{2}t,1,\begin{pmatrix}0\\0\\\end{pmatrix},\begin{pmatrix}0\\0\\\end{pmatrix},\begin{pmatrix}1&0\\0&1\\\end{pmatrix}\right)$$. The true value is $$(\theta_0,\theta_1,\theta_2)=(1,2,3)$$. ## set parameter values in noise lambda<-1/2 alpha<-1 beta<-c(0,0) mu<-c(0,0) Lambda<-matrix(c(1,0,0,1),2,2) ## set coefficients dri<-c("1-theta0*x1-x2","-theta1*x2") jum<-matrix(c("theta2/(1+x1^2)+1","0","1","1"),2,2) yuima <- setModel(drift=dri, jump.coeff=jum, solve.variable=c("x1","x2"), state.variable = c("x1","x2"), measure.type="code", measure=list(df="rvgamma(z, lambda, alpha, beta, mu, Lambda)")) n<-100000; tp<-0.1; N<-floor(n*tp); T<-100; hn<-T/N sam<-setSampling(Terminal = T, n=n) subsam<-setSampling(Terminal = T, n=N) yuima<-setYuima(model = yuima, sampling = sam) true<-list(theta0 = 1,theta1 = 2,theta2 = 3,lambda=lambda, alpha=alpha, beta=beta, mu=mu, Lambda=Lambda) upper<-list(theta0 = 4, theta1 = 4, theta2 = 5) lower<-list(theta0 = 0.5, theta1 = 1, theta2 = 1) set.seed(123) yuima<-simulate(yuima, xinit = c(0,0), true.parameter = true, sampling = sam, subsampling = subsam) ## generate a path plot(yuima, col="blue") # dev.copy2pdf(file="qmle-2vg.pdf", width=8, height=5) ## figure output (scale: inch) start<-list(theta0 = runif(1,0.5,4), theta1 = runif(1,1,4), theta2 = runif(1,1,5)) qmleLevy(yuima,start=start,lower=lower,upper=upper,joint = FALSE,third = TRUE)  ##first
##   theta2
## 2.608618
##
## $second ## theta0 theta1 ## 0.9527039 1.9516423 ## ##$third
##   theta2
## 2.598297

## 2-dimensional NIG driven SDE

\begin{aligned} d\begin{pmatrix}X_{1,t}\\X_{2,t}\end{pmatrix}=\begin{pmatrix}1-\theta_0X_{1,t}\\-\theta_1X_{2,t}\end{pmatrix}dt+\begin{pmatrix} \exp\left(-\frac{\theta_2}{1+X_{1,t}^2}\right)& 0\\1 & \exp\left(-\frac{\theta_3}{\sqrt{1+X_{2,t}^2}}\right) \end{pmatrix} dZ_t, \end{aligned} where the driving noise $$Z_t$$ obeys $$NIG_2\left(\frac{1}{\sqrt{\pi}}t,1, \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&0\\0&1\end{pmatrix}\right)$$. The true value is $$(\theta_0,\theta_1,\theta_2,\theta_3)=(1,2,3,4)$$.

## set parameters in noise
alpha<-0.5
a<-1/sqrt(pi)
b<-1
beta<-c(0,0)
mu<-c(0,0)
Lambda<-matrix(c(1,0,0,1),2,2)

dri<-c("1-theta0*x1","-theta1*x2")
jum<-matrix(c("exp(-theta2/(1+x1^2))","0","1",
"exp(-theta3/sqrt(1+x2^2))"),2,2) ## set coefficients

yuima <- setModel(drift=dri, jump.coeff=jum,
solve.variable=c("x1","x2"), state.variable = c("x1","x2"),
measure.type="code",
measure=list(df="rnts(z,alpha,a,b,beta,mu,Lambda)"))

n<-100000; tp<-0.1; N<-floor(n*tp); T<-100; hn<-T/N

sam<-setSampling(Terminal = T, n=n)
subsam<-setSampling(Terminal = T, n=N)
yuima<-setYuima(model = yuima, sampling = sam)

true<-list(theta0 = 1,theta1 = 2,theta2 = 3,theta3 = 4,alpha=alpha, a=a, b=b,
beta=beta, mu=mu, Lambda=Lambda)
upper<-list(theta0 = 4, theta1 = 4, theta2 = 5, theta3 = 5)
lower<-list(theta0 = 0.5, theta1 = 1, theta2 = 1, theta3 = 1)

set.seed(567)
Y<-simulate(yuima, xinit = c(0,0), true.parameter = true,
sampling = sam, subsampling = subsam)
plot(Y, col="blue")

# dev.copy2pdf(file="qmle-2vg.pdf", width=8, height=5) ## figure output (scale: inch)

start<-list(theta0 = runif(1,0.5,4), theta1 = runif(1,1,4),
theta2 = runif(1,1,5), theta3 = runif(1,1,5))
qmleLevy(Y,start=start,lower=lower,upper=upper,joint = FALSE, third = TRUE) 
## $first ## theta2 theta3 ## 2.773936 4.020123 ## ##$second
##    theta0    theta1
## 0.9628266 1.9959077
##
## \$third
##   theta2   theta3
## 2.894978 4.024005