[Yuima-commits] r623 - in pkg/yuima: . inst/ybook
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 3 18:02:47 CEST 2017
Author: iacus
Date: 2017-09-03 18:02:47 +0200 (Sun, 03 Sep 2017)
New Revision: 623
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/inst/ybook/chapter1.R
pkg/yuima/inst/ybook/chapter2.R
pkg/yuima/inst/ybook/chapter3.R
pkg/yuima/inst/ybook/chapter4.R
Log:
changed ybook scripts
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-08-16 08:28:45 UTC (rev 622)
+++ pkg/yuima/DESCRIPTION 2017-09-03 16:02:47 UTC (rev 623)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.6.9
+Version: 1.7.0
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/inst/ybook/chapter1.R
===================================================================
--- pkg/yuima/inst/ybook/chapter1.R 2017-08-16 08:28:45 UTC (rev 622)
+++ pkg/yuima/inst/ybook/chapter1.R 2017-09-03 16:02:47 UTC (rev 623)
@@ -11,6 +11,9 @@
options(width=55)
options(continue=" ")
require(yuima)
+Rver <- paste(version$major,version$minor, collapse="",sep=".")
+YUIMAver <- as.character(read.dcf(file=system.file("DESCRIPTION", package="yuima"),
+ fields="Version"))
## ----eval=FALSE----------------------------------------------------------
## ybook(1)
@@ -186,7 +189,7 @@
## ------------------------------------------------------------------------
sol <- c("x1","x2") # variable for numerical solution
b <- c("-theta*x1","-x1-gamma*x2") # drift vector
-s <- matrix(c("1","x1","0","delta","x2","0"),2,3) # diff. mat.
+s <- matrix(c("1","x1","0","beta","x2","0"),2,3) # diff. mat.
mymod <- setModel(drift = b, diffusion = s, solve.variable = sol)
## ------------------------------------------------------------------------
@@ -201,7 +204,7 @@
## ------------------------------------------------------------------------
set.seed(123)
X2 <- simulate(mymod, sampling=samp,
- true.param=list(theta=1,gamma=1,delta=1))
+ true.param=list(theta=1,gamma=1,beta=1))
X2
## ----results='hide'------------------------------------------------------
@@ -260,10 +263,10 @@
set.seed(123)
Y.sub <- simulate(mymod,sampling=setSampling(delta=0.001,n=1000),
subsampling=setSampling(delta=0.01,n=100),
- true.par=list(theta=1,delta=1,gamma=1))
+ true.par=list(theta=1,beta=1,gamma=1))
set.seed(123)
Y <- simulate(mymod, sampling=setSampling(delta=0.001,n=1000),
- true.par=list(theta=1,delta=1,gamma=1))
+ true.par=list(theta=1,beta=1,gamma=1))
plot(Y, plot.type="single")
points(get.zoo.data(Y.sub)[[1]],col="red")
points(get.zoo.data(Y.sub)[[2]],col="green",pch=18)
@@ -293,26 +296,21 @@
## my.yuima <- setYuima(data=setData(X), model=mod)
## ----echo=TRUE,results='hide',tidy=TRUE----------------------------------
-data <- read.csv("http://chart.yahoo.com/table.csv?s=IBM&g=d&x=.csv",
- stringsAsFactor=FALSE)
+require(quantmod)
+getSymbols("IBM", to = "2017-07-31")
## ----echo=TRUE-----------------------------------------------------------
-head(data)
+str(IBM)
-## ----results='hide',tidy=FALSE-------------------------------------------
-tmp <- zoo(data$Close, order.by=as.Date(data$Date,
- format="%Y-%m-%d"))
-str(tmp)
+## ----echo=TRUE-----------------------------------------------------------
+head(IBM)
-## ----echo=FALSE----------------------------------------------------------
-writeLines(strwrap(capture.output(str(tmp)),width=60))
-
## ------------------------------------------------------------------------
-x <- setYuima(data=setData(tmp))
+x <- setYuima(data=setData(IBM$IBM.Close))
str(x at data)
## ------------------------------------------------------------------------
-y <- setYuima(data=setData(tmp, delta=1/252))
+y <- setYuima(data=setData(IBM$IBM.Close, delta=1/252))
str(y at data)
## ----setData,fig.keep='none'---------------------------------------------
@@ -332,11 +330,11 @@
## ----quantmod,message=FALSE----------------------------------------------
library(quantmod)
-getSymbols("IBM")
+getSymbols("IBM", to = "2017-07-31")
attr(IBM, "src")
## ------------------------------------------------------------------------
-getSymbols("IBM",src="google")
+getSymbols("IBM", to = "2017-07-31", src="google")
attr(IBM, "src")
## ------------------------------------------------------------------------
@@ -346,22 +344,7 @@
attr(EURUSD, "src")
str(EURUSD)
-## ----results="hide",message=FALSE----------------------------------------
-require(fImport)
-X <- yahooSeries("IBM")
-str(X)
-
-## ----echo=FALSE----------------------------------------------------------
-writeLines(strwrap(capture.output(str(X)),width=60))
-
## ----results='hide'------------------------------------------------------
-X <- yahooImport("IBM")
-str(X)
-
-## ----echo=FALSE----------------------------------------------------------
-writeLines(strwrap(capture.output(str(X)),width=60))
-
-## ----results='hide'------------------------------------------------------
library(tseries)
x <- get.hist.quote("IBM")
str(x)
@@ -369,14 +352,6 @@
## ----echo=FALSE----------------------------------------------------------
writeLines(strwrap(capture.output(str(x)),width=60))
-## ----results='hide'------------------------------------------------------
-x <- get.hist.quote(instrument = "EUR/USD", provider = "oanda",
- start = Sys.Date() - 300)
-str(x)
-
-## ----echo=FALSE----------------------------------------------------------
-writeLines(strwrap(capture.output(str(x)),width=60))
-
## ------------------------------------------------------------------------
mydat <- get.zoo.data(y)[[1]]
str(mydat)
@@ -434,7 +409,7 @@
## ------------------------------------------------------------------------
Xreg <- zooreg(some.data, start = c(1964, 2), frequency = 12)
-time(Xreg)[1:12]
+time(Xreg)
## ------------------------------------------------------------------------
Y <- as.ts(X)
@@ -659,6 +634,19 @@
quotes[ (time(quotes) >= initial) & (time(quotes)<= terminal), 4:9]
## ------------------------------------------------------------------------
+getSymbols("IBM", from="2015-01-01", to = "2016-12-31")
+str(IBM)
+
+## ------------------------------------------------------------------------
+IBM["2015-01","IBM.Close"]
+
+## ------------------------------------------------------------------------
+IBM["2016-02-11/2016-03-05","IBM.Close"]
+
+## ------------------------------------------------------------------------
+IBM["/2015-02-11","IBM.Close"]
+
+## ------------------------------------------------------------------------
mod2 <- setModel(drift = "-mu*x", diffusion = "1/(1+x^gamma)")
mod2
Modified: pkg/yuima/inst/ybook/chapter2.R
===================================================================
--- pkg/yuima/inst/ybook/chapter2.R 2017-08-16 08:28:45 UTC (rev 622)
+++ pkg/yuima/inst/ybook/chapter2.R 2017-09-03 16:02:47 UTC (rev 623)
@@ -61,6 +61,7 @@
plot(x1, main="mod1, xinit=1")
plot(x2, main="mod1, xinit=expression(rnorm(1))")
plot(x3, main="mod2, xinit=3")
+par(mfrow=c(1,1))
## ----echo=FALSE,results='hide'-------------------------------------------
pdf("figures/plot-mod1diff2.pdf",width=9,height=4)
@@ -96,7 +97,7 @@
## ------------------------------------------------------------------------
hyper1 <- setModel( diff="sigma",
- drift="(sigma/2)^2*(beta-alpha*((x-mu)/(sqrt(delta^2+(x-mu)^2))))")
+ drift="(sigma^2/2)*(beta-alpha*((x-mu)/(sqrt(delta^2+(x-mu)^2))))")
## ------------------------------------------------------------------------
hyper1
@@ -104,7 +105,7 @@
## ------------------------------------------------------------------------
hyper2 <- setModel(drift="0",
- diffusion = "sigma*alpha*exp(0.5*sqrt(delta^2+(x-mu)^2)-
+ diffusion = "sigma*exp(0.5*alpha*sqrt(delta^2+(x-mu)^2)-
beta*(x-mu))")
## ------------------------------------------------------------------------
@@ -271,7 +272,7 @@
lower <- list(theta1=0,theta2=0)
upper <- list(theta1=1,theta2=1)
bayes1 <- adaBayes(yuima, start=param.init, prior=prior,
- lower=lower,upper=upper)
+ lower=lower,upper=upper, method="nomcmc")
## ----echo=TRUE-----------------------------------------------------------
coef(summary(bayes1))
@@ -390,7 +391,7 @@
coef(fit)
## ------------------------------------------------------------------------
-getSymbols("DEXUSEU", src="FRED")
+getSymbols("DEXUSEU", src="FRED", to="2016-12-31")
head(DEXUSEU)
meanCIR <- mean(DEXUSEU, na.rm=TRUE)
meanCIR
@@ -447,43 +448,37 @@
model<- setModel(drift="t1*(t2-x)",diffusion="t3")
## ------------------------------------------------------------------------
-T<-10
-n<-1000
-sampling <- setSampling(Terminal=T,n=n)
+T<-300
+n<-3000
+sampling <- setSampling(Terminal=T, n=n)
yuima<-setYuima(model=model, sampling=sampling)
-h0 <- list(t1=0.3, t2=1, t3=0.25)
-h1 <- list(t1=0.3, t2=0.2, t3=0.1)
+h00 <- list(t1=0.3, t2=1, t3=0.25)
+h01 <- list(t1=0.3, t2=0.2, t3=0.1)
set.seed(123)
-X <- simulate(yuima, xinit=1, true=h0)
+X <- simulate(yuima, xinit=1, true.par=h00)
## ------------------------------------------------------------------------
phi1 <- function(x) 1-x+x*log(x)
## ------------------------------------------------------------------------
-phi.test(X, H0=h0, H1=h1,phi=phi1)
-
-## ------------------------------------------------------------------------
-phi.test(X, H0=h0, phi=phi1, start=h0,
+phi.test(X, H0=h00, phi=phi1, start=h00,
lower=list(t1=0.1, t2=0.1, t3=0.1),
upper=list(t1=2,t2=2,t3=2),method="L-BFGS-B")
-## ----echo=FALSE----------------------------------------------------------
-pval <- phi.test(X, H0=h0, phi=phi1, start=h0,
+## ----echo=FALSE,results='hide'-------------------------------------------
+pval <- phi.test(X, H0=h00, phi=phi1, start=h00,
lower=list(t1=0.1, t2=0.1, t3=0.1),
upper=list(t1=2,t2=2,t3=2),method="L-BFGS-B")$pvalue
## ------------------------------------------------------------------------
-phi.test(X, H0=h1, phi=phi1, start=h0,
+phi.test(X, H0=h01, phi=phi1, start=h00,
lower=list(t1=0.1, t2=0.1, t3=0.1),
upper=list(t1=2,t2=2,t3=2),method="L-BFGS-B")
## ------------------------------------------------------------------------
-phi.test(X, H0=h0, H1=h1)
-
-## ------------------------------------------------------------------------
library(quantmod)
Delta <- 1/252
-getSymbols("DEXUSEU", src="FRED")
+getSymbols("DEXUSEU", src="FRED", to = "2016-12-31")
meanCIR <- mean(DEXUSEU, na.rm=TRUE)
gBm <- setModel(drift="mu*x", diffusion="sigma*x")
mod <- setYuima(model=gBm, data=setData(na.omit(DEXUSEU),
@@ -765,7 +760,7 @@
## ----fig.keep='none'-----------------------------------------------------
X <- diff(log(get.zoo.data(mod)[[1]]))
plot(X)
-abline(v=cp$tau, lty=3)
+abline(v=cp$tau, lty=3,lwd=2,col="red")
## ----plot-returns,echo=FALSE, fig.keep='none',results='hide'-------------
pdf("figures/plot-returns.pdf",width=9,height=4)
@@ -803,7 +798,6 @@
## ------------------------------------------------------------------------
set.seed(123)
-
Terminal <- 1
n <- 1000
# Cumulative Covariance
@@ -908,13 +902,11 @@
# intentionally displace the second time series
data1 <- get.zoo.data(yuima)[[1]]
-
data2 <- get.zoo.data(yuima)[[2]]
time2 <- time( data2 )
theta2 <- 0.05 # the lag of x2 behind x1
stime2 <- time2 + theta2
time(data2) <- stime2
-
data3 <- get.zoo.data(yuima)[[3]]
time3 <- time( data3 )
theta3 <- 0.12 # the lag of x3 behind x1
@@ -937,6 +929,13 @@
llag(yuima)
llag(syuima)
+## ----plot-shifted-ci,echo=FALSE, fig.keep='none',results='hide'----------
+pdf("figures/plot-shifted-ci.pdf",width=9,height=5)
+par(mar=c(4,5,2,1))
+par(mfrow=c(1,3))
+llag(syuima,plot=TRUE,ci=TRUE)
+dev.off()
+
## ------------------------------------------------------------------------
data2 <- get.zoo.data(yuima)[[2]]
time2 <- time( data2 )
Modified: pkg/yuima/inst/ybook/chapter3.R
===================================================================
--- pkg/yuima/inst/ybook/chapter3.R 2017-08-16 08:28:45 UTC (rev 622)
+++ pkg/yuima/inst/ybook/chapter3.R 2017-09-03 16:02:47 UTC (rev 623)
@@ -64,20 +64,20 @@
poisson4
## ----mod5,fig.keep='none'------------------------------------------------
-mod5 <- setPoisson(intensity="alpha+lambda*t",
+mod5 <- setPoisson(intensity="alpha+beta*t",
df=list("dnorm(z,mu,sigma)"))
set.seed(123)
poisson5 <- simulate(mod5, sampling=samp,
- true.par=list(alpha=1,lambda=.5,mu=0, sigma=2))
+ true.par=list(alpha=2,beta=.5,mu=0, sigma=2))
plot(poisson5)
f <- function(t,alpha,beta) alpha + beta*t
-curve(f(x,alpha=2,beta=0.3)-20,0,30,add=TRUE,col="red",lty=3)
+curve(f(x,alpha=2,beta=0.5)-20,0,30,add=TRUE,col="red",lty=3,lwd=2)
## ----plot-poi5,echo=FALSE,results='hide'---------------------------------
pdf("figures/plot-poi5.pdf",width=9,height=4)
par(mar=c(4,4,1,1))
plot(poisson5,type="S")
-curve(f(x,alpha=2,beta=0.3)-20,0,30,add=TRUE,col="red",lty=3)
+curve(f(x,alpha=2,beta=0.5)-20,0,30,add=TRUE,col="red",lty=3,lwd=2)
dev.off()
## ----mod6,fig.keep='none'------------------------------------------------
@@ -88,13 +88,13 @@
true.par=list(theta=1.5,mu=0, sigma=2))
plot(poisson6)
f <- function(t,theta) theta*t^(theta-1)
-curve(f(x,theta=1.5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,theta=1.5),0,30,add=TRUE,col="red",lty=3,lwd=2)
## ----plot-poi6,echo=FALSE,results='hide'---------------------------------
pdf("figures/plot-poi6.pdf",width=9,height=4)
par(mar=c(4,4,1,1))
plot(poisson6,type="S")
-curve(f(x,theta=1.5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,theta=1.5),0,30,add=TRUE,col="red",lty=3,lwd=2)
dev.off()
## ----mod7,fig.keep='none'------------------------------------------------
@@ -105,13 +105,13 @@
true.par=list(lambda=.2,beta=10,gamma=1))
plot(poisson7)
f <- function(t,beta,lambda) beta*exp(-lambda*t)
-curve(f(x,beta=10,lambda=0.5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,beta=10,lambda=0.2),0,30,add=TRUE,col="red",lty=3,lwd=2)
## ----plot-poi7,echo=FALSE,results='hide'---------------------------------
pdf("figures/plot-poi7.pdf",width=9,height=4)
par(mar=c(4,4,1,1))
plot(poisson7,type="S")
-curve(f(x,beta=10,lambda=0.5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,beta=10,lambda=0.2),0,30,add=TRUE,col="red",lty=3,lwd=2)
dev.off()
## ----mod8,fig.keep='none'------------------------------------------------
@@ -123,13 +123,13 @@
plot(poisson8)
f <- function(t,a,omega,phi,lambda) 0.5*a*(1+cos(omega*t+phi))+lambda
curve(f(x,a=2,omega=0.5,phi=3.14,lambda=5),0,30,add=TRUE,
- col="red",lty=3)
+ col="red",lty=3,lwd=2)
## ----plot-poi8,echo=FALSE,results='hide'---------------------------------
pdf("figures/plot-poi8.pdf",width=9,height=4)
par(mar=c(4,4,1,1))
plot(poisson8,type="S")
-curve(f(x,a=2,omega=0.5,phi=3.14,lambda=5),0,30,add=TRUE, col="red",lty=3)
+curve(f(x,a=2,omega=0.5,phi=3.14,lambda=5),0,30,add=TRUE, col="red",lty=3,lwd=2)
dev.off()
## ----fig.keep='none'-----------------------------------------------------
@@ -140,13 +140,13 @@
true.par=list(a=1,theta=0.5,lambda=5,mu=0,sigma=1))
plot(poisson9)
f <- function(t,a,theta,lambda) a*cos(theta*t)+lambda
-curve(f(x,a=1,theta=0.5,lambda=5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,a=1,theta=0.5,lambda=5),0,30,add=TRUE,col="red",lty=3,lwd=2)
## ----plot-poi9,echo=FALSE,results='hide'---------------------------------
pdf("figures/plot-poi9.pdf",width=9,height=4)
par(mar=c(4,4,1,1))
plot(poisson9,type="S")
-curve(f(x,a=1,theta=0.5,lambda=5),0,30,add=TRUE,col="red",lty=3)
+curve(f(x,a=1,theta=0.5,lambda=5),0,30,add=TRUE,col="red",lty=3,lwd=2)
dev.off()
## ----mod10,fig.keep='none'-----------------------------------------------
@@ -189,7 +189,9 @@
Lambda <- matrix(c(1,0,0,1),2,2)
t(rNIG(n,alpha=alpha,beta=beta,delta=delta0,mu=mu,Lambda=Lambda))
}
-d2DNIG <- function(n,alpha){
+# the next fake density plays no role in simulation
+# but it is needed for model specification
+d2DNIG <- function(n,alpha){
rep(0,2)
}
Modified: pkg/yuima/inst/ybook/chapter4.R
===================================================================
--- pkg/yuima/inst/ybook/chapter4.R 2017-08-16 08:28:45 UTC (rev 622)
+++ pkg/yuima/inst/ybook/chapter4.R 2017-09-03 16:02:47 UTC (rev 623)
@@ -18,15 +18,16 @@
sigma <- 1
lambda <- 10
samp <- setSampling(Terminal=10, n=1000)
-mod10 <- setPoisson(intensity="lambda", df=list("dnorm(z,mu,sigma)"))
-y10 <- simulate(mod10,sampling=samp,
+mod10b <- setPoisson(intensity="lambda", df=list("dnorm(z,mu,sigma)"))
+y10b <- simulate(mod10b,sampling=samp,
true.par=list(lambda=lambda,mu=0.1, sigma=2))
-y10
+y10b
## ----fig.keep='none'-----------------------------------------------------
BGmodel <- setModel(drift="0", xinit="0", jump.coeff="1",
- measure.type="code", measure=list(df="rbgamma(z, delta.minus=2,
- gamma.minus=0.6, delta.plus=1.4, gamma.plus=0.3)"))
+ measure.type="code", measure=list(df="rbgamma(z, delta.plus=1.4,
+ gamma.plus=0.3, delta.minus=2,
+ gamma.minus=0.6)"))
n <- 1000
samp <- setSampling(Terminal=1, n=n)
BGyuima <- setYuima(model=BGmodel, sampling=samp)
@@ -285,7 +286,7 @@
## }
## hist(VGPsimdata,xlim=c(-7,10),ylim=c(0,0.22),breaks=100,freq=FALSE,
## main=expression(paste("Distribution of ",X[1.8],
-## " and Density of NG")))
+## " and Density of VG")))
## curve(dvgamma(x,lambda*T,alpha,beta,mu*T),add=TRUE,col="red")
## ----plot-VGPproc2,echo=FALSE,results='hide'-----------------------------
@@ -302,7 +303,7 @@
VGPsimdata <- c(VGPsimdata,as.numeric(x1[1]))
}
hist(VGPsimdata,xlim=c(-7,10),ylim=c(0,0.22),breaks=100,freq=FALSE,
- main=expression(paste("Distribution of ",X[1.8], " and Density of NG")))
+ main=expression(paste("Distribution of ",X[1.8], " and Density of VG")))
curve(dvgamma(x,lambda*T,alpha,beta,mu*T),add=TRUE,col="red")
dev.off()
@@ -318,7 +319,7 @@
normal.rn <- rnorm(n,0,1)
iv.rn <- rIG(n,delta*T,gamma)
z <- mu*T+beta*iv.rn+sqrt(iv.rn)*normal.rn
-title <- expression(paste(NIG[1.8],
+title <- expression(paste(NIGP[1.8],
" built by subordination (green) and rNIG (white)"))
nig.rn <- rNIG(n,alpha,beta,delta*T,mu*T)
hist(z,xlim=c(-1,10),ylim=c(0,0.61),breaks=100, freq=FALSE,
@@ -432,16 +433,16 @@
main=expression(Y[t]),probability=TRUE,col="red")
## experiment by convolution
nrep <- 3000000
-X1 <- rnts(nrep,alpha,delta*t,gamma,beta,mu*t,Lambda)
+Xt <- rnts(nrep,alpha,delta*t,gamma,beta,mu*t,Lambda)
X05 <- rnts(nrep,alpha,delta*t/2,gamma,beta,mu*t/2,Lambda)
X05.prime <- rnts(nrep,alpha,delta*t/2,gamma,beta,mu*t/2,Lambda)
Xsum <- X05+X05.prime
-hist(X1,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
- main=expression(X[1]),probability=TRUE)
+hist(Xt,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
+ main=expression(X[t]),probability=TRUE)
hist(Xsum,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
main=expression(paste(X[t/2]+X[t/2],"'")),
probability=TRUE,col="red")
-ks.test(X1,Xsum)
+ks.test(Xt,Xsum)
## ----plot-NTSPproc,echo=FALSE,results='hide'-----------------------------
pdf("figures/plot-NTSPproc.pdf",width=9,height=4)
@@ -451,162 +452,17 @@
main=expression(X[t]),probability=TRUE)
hist(y,xlim=c(-3,3),ylim=c(0,1.2),breaks=200,
main=expression(Y[t]),probability=TRUE,col="red")
-hist(X1,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
-main=expression(X[1]),probability=TRUE)
+hist(Xt,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
+main=expression(X[t]),probability=TRUE)
hist(Xsum,xlim=c(-3,3),ylim=c(0,1.2),breaks=300,
main=expression(paste(X[t/2]+X[t/2],"'")),
probability=TRUE,col="red")
dev.off()
-rm(X1)
+rm(Xt)
rm(X05)
rm(X05.prime)
rm(Xsum)
-## ----eval=FALSE----------------------------------------------------------
-## alpha <- 0.5
-## beta <- -0.4
-## sigma <- 0.7
-## gamma <- 0.5
-## n <- 1000
-## T <- 1.8
-## ASmodel <- setModel(drift=0, jump.coeff=1, measure.type="code",
-## measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-## samp <- setSampling(Terminal=T, n=n)
-## ASyuima <- setYuima(model=ASmodel, sampling=samp)
-## set.seed(129)
-## for (i in 1:10) {
-## result <- simulate(ASyuima, true.par=list(alpha=alpha,
-## beta=beta,sigma=sigma,gamma=gamma))
-## plot(result,xlim=c(0,T),ylim=c(-40,10),col=i,
-## main=expression(paste("Paths of stable process (",
-## alpha==0.5,",",beta==-0.4,")")),par(new=T))
-## }
-##
-## #param2
-## alpha <- 1
-## beta <- -0.4
-## sigma <- 0.7
-## gamma <- 0.5
-## AS2model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-## measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-## AS2yuima <- setYuima(model=AS2model, sampling=samp)
-## for (i in 1:10) {
-## result <- simulate(AS2yuima, true.par=list(alpha=alpha,
-## beta=beta,sigma=sigma,gamma=gamma))
-## plot(result,xlim=c(0,T),ylim=c(-5,5),col=i,
-## main=expression(paste("Paths of stable process (",
-## alpha==1,",",beta==-0.4,")")),par(new=T))
-## }
-##
-## #param3
-## alpha <- 1
-## beta <- 0.4
-## sigma <- 0.7
-## gamma <- 0.5
-## AS3model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-## measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-## AS3yuima <- setYuima(model=AS3model, sampling=samp)
-## for (i in 1:10) {
-## result <- simulate(AS3yuima, true.par=list(alpha=alpha,
-## beta=beta,sigma=sigma,gamma=gamma))
-## plot(result,xlim=c(0,T),ylim=c(-5,5),col=i,
-## main=expression(paste("Paths of stable process (",
-## alpha==1,",",beta==0.4,")")),par(new=T))
-## }
-##
-## #param4
-## alpha <- 1.5
-## beta <- 0.4
-## sigma <- 0.7
-## gamma <- 0.5
-## AS4model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-## measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-## AS4yuima <- setYuima(model=AS4model, sampling=samp)
-## for (i in 1:10) {
-## result <- simulate(AS4yuima, true.par=list(alpha=alpha,
-## beta=beta, sigma=sigma,gamma=gamma))
-## plot(result,xlim=c(0,T),ylim=c(-3,5),col=i,
-## main=expression(paste("Paths of stable process (",
-## alpha==1.5,",",beta==0.4,")")),par(new=T))
-## }
-
-## ----plot-ASproc,echo=FALSE,results='hide'-------------------------------
-pdf("figures/plot-ASproc1.pdf",width=9,height=4)
-par(mar=c(4,4,2,1))
-alpha <- 0.5
-beta <- -0.4
-sigma <- 0.7
-gamma <- 0.5
-n <- 1000
-T <- 1.8
-ASmodel <- setModel(drift=0, jump.coeff=1, measure.type="code",
-measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-samp <- setSampling(Terminal=T, n=n)
-ASyuima <- setYuima(model=ASmodel, sampling=samp)
-set.seed(129)
-for (i in 1:10) {
-result <- simulate(ASyuima, true.par=list(alpha=alpha,
-beta=beta,sigma=sigma,gamma=gamma))
-plot(result,xlim=c(0,T),ylim=c(-40,10),col=i,
-main=expression(paste("Paths of stable process (",
-alpha==0.5,",",beta==-0.4,")")),par(new=T))
-}
-dev.off()
-pdf("figures/plot-ASproc2.pdf",width=9,height=4)
-par(mar=c(4,4,2,1))
-#param2
-alpha <- 1
-beta <- -0.4
-sigma <- 0.7
-gamma <- 0.5
-AS2model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-AS2yuima <- setYuima(model=AS2model, sampling=samp)
-for (i in 1:10) {
-result <- simulate(AS2yuima, true.par=list(alpha=alpha,
-beta=beta,sigma=sigma,gamma=gamma))
-plot(result,xlim=c(0,T),ylim=c(-5,5),col=i,
-main=expression(paste("Paths of stable process (",
-alpha==1,",",beta==-0.4,")")),par(new=T))
-}
-dev.off()
-pdf("figures/plot-ASproc3.pdf",width=9,height=4)
-par(mar=c(4,4,2,1))
-#param3
-alpha <- 1
-beta <- 0.4
-sigma <- 0.7
-gamma <- 0.5
-AS3model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-AS3yuima <- setYuima(model=AS3model, sampling=samp)
-for (i in 1:10) {
-result <- simulate(AS3yuima, true.par=list(alpha=alpha,
-beta=beta,sigma=sigma,gamma=gamma))
-plot(result,xlim=c(0,T),ylim=c(-5,5),col=i,
-main=expression(paste("Paths of stable process (",
-alpha==1,",",beta==0.4,")")),par(new=T))
-}
-dev.off()
-pdf("figures/plot-ASproc4.pdf",width=9,height=4)
-par(mar=c(4,4,2,1))
-#param4
-alpha <- 1.5
-beta <- 0.4
-sigma <- 0.7
-gamma <- 0.5
-AS4model <- setModel(drift=0, jump.coeff=1, measure.type="code",
-measure=list(df="rstable(z,alpha,beta,sigma,gamma)"))
-AS4yuima <- setYuima(model=AS4model, sampling=samp)
-for (i in 1:10) {
-result <- simulate(AS4yuima, true.par=list(alpha=alpha,beta=beta,
-sigma=sigma,gamma=gamma))
-plot(result,xlim=c(0,T),ylim=c(-3,5),col=i,
-main=expression(paste("Paths of stable process (",
-alpha==1.5,",",beta==0.4,")")),par(new=T))
-}
-dev.off()
-
## ----fig.keep='none'-----------------------------------------------------
modJump <- setModel(drift = c("-theta*x"), diffusion = "sigma",
jump.coeff=c("gamma+x/sqrt(1+x^2)"),
More information about the Yuima-commits
mailing list