## 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)
## 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)
## 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)
## 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)
## 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)
## 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)
## 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)
\[ \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)
## $joint
## theta1 theta0
## 1.9739888 0.9926288
\[ \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 \(NIG\left(\delta,0,\delta t,0,1\right)\) with \(\delta=10\). True value \((\theta_{0,0},\theta_{1,0})=(1,2)\).
# Model construction
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="rNIG(z,10,0,10,0)"))
n<-100000; tp<-0.1; N<-n*tp; T<-100; hn<-T/N
sam<-setSampling(Terminal = T, n=n) ## set sampling scheme
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 = 1, true.parameter = true,sampling = sam,
subsampling = subsam) ## generate a path
plot(yuima, col="blue")
# dev.copy2pdf(file="qmle-nig.pdf", width=8, height=5) ## figure output (scale: inch)
upper1<-list(theta0 = 4, theta1 = 4)
lower1<-list(theta0 = 0.5, theta1 = 1)
start1<-list(theta0 = runif(1,0.5,4), theta1 = runif(1,1,4))
qmleLevy(yuima, start=start1, lower=lower1, upper=upper1, joint = FALSE, third = FALSE)
## $first
## theta1
## 2.000248
##
## $second
## theta0
## 0.9820445
\[ \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
\[ \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