[Yuima-commits] r636 - in pkg/yuima: . R inst/ybook
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 2 16:12:53 CET 2018
Author: iacus
Date: 2018-02-02 16:12:53 +0100 (Fri, 02 Feb 2018)
New Revision: 636
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qgv.R
pkg/yuima/inst/ybook/chapter2.R
pkg/yuima/inst/ybook/chapter4.R
Log:
updated for book
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2018-01-16 13:38:47 UTC (rev 635)
+++ pkg/yuima/DESCRIPTION 2018-02-02 15:12:53 UTC (rev 636)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.7.6
+Version: 1.7.7
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2)
Author: YUIMA Project Team
Modified: pkg/yuima/R/qgv.R
===================================================================
--- pkg/yuima/R/qgv.R 2018-01-16 13:38:47 UTC (rev 635)
+++ pkg/yuima/R/qgv.R 2018-02-02 15:12:53 UTC (rev 636)
@@ -143,7 +143,7 @@
nconst<-as.character(yuima at model@diffusion[[1]])
}
} else {
- nconst<-as.character(yuima at model@diffusion[[1]])
+ nconst<- as.character(yuima at model@diffusion[[1]])
}
@@ -205,6 +205,7 @@
}
+nconst <- gsub("[()]", "", nconst)
x<-c(H,C)
names(x)<-c("hurst",nconst)
sdx<-diag(c(sdH,sdC))
Modified: pkg/yuima/inst/ybook/chapter2.R
===================================================================
--- pkg/yuima/inst/ybook/chapter2.R 2018-01-16 13:38:47 UTC (rev 635)
+++ pkg/yuima/inst/ybook/chapter2.R 2018-02-02 15:12:53 UTC (rev 636)
@@ -391,7 +391,8 @@
coef(fit)
## ------------------------------------------------------------------------
-getSymbols("DEXUSEU", src="FRED", to="2016-12-31")
+getSymbols("DEXUSEU", src="FRED")
+DEXUSEU <- DEXUSEU["/2016"]
head(DEXUSEU)
meanCIR <- mean(DEXUSEU, na.rm=TRUE)
meanCIR
@@ -478,21 +479,19 @@
## ------------------------------------------------------------------------
library(quantmod)
Delta <- 1/252
-getSymbols("DEXUSEU", src="FRED", to = "2016-12-31")
-meanCIR <- mean(DEXUSEU, na.rm=TRUE)
+getSymbols("DEXUSEU", src="FRED")
+DEXUSEU <- DEXUSEU["/2016"]
+USEU <- setData(na.omit(DEXUSEU), delta=Delta)
+meanCIR <- mean(get.zoo.data(USEU)[[1]])
gBm <- setModel(drift="mu*x", diffusion="sigma*x")
-mod <- setYuima(model=gBm, data=setData(na.omit(DEXUSEU),
- delta=Delta))
+mod <- setYuima(model=gBm, data=USEU)
cir1 <- setModel(drift="theta1-theta2*x", diffusion="sigma*sqrt(x)")
cir2 <- setModel(drift="kappa*(mu-x)", diffusion="sigma*sqrt(x)")
ckls <- setModel(drift="theta1-theta2*x", diffusion="sigma*x^gamma")
-mod1 <- setYuima(model=cir1, data=setData(na.omit(DEXUSEU),
- delta=Delta))
-mod2 <- setYuima(model=cir2, data=setData(na.omit(DEXUSEU),
- delta=Delta))
-mod3 <- setYuima(model=ckls, data=setData(na.omit(DEXUSEU),
- delta=Delta))
- gBm.fit <- qmle(mod, start=list(mu=1, sigma=1),
+mod1 <- setYuima(model=cir1, data=USEU)
+mod2 <- setYuima(model=cir2, data=USEU)
+mod3 <- setYuima(model=ckls, data=USEU)
+gBm.fit <- qmle(mod, start=list(mu=1, sigma=1),
lower=list(mu=0.1, sigma=0.1),
upper=list(mu=100, sigma=10))
cir1.fit <- qmle(mod1, start=list(theta1=1, theta2=1, sigma=0.5),
@@ -546,10 +545,10 @@
tmp <- AIC(gBm.fit,cir1.fit,cir2.fit,ckls.fit)
## ------------------------------------------------------------------------
-alpha <- c("1-mu11*S1+mu12*S2","2+mu21*S1-mu22*S2")
-beta <- matrix(c("s1*S1","s2*S1", "-s3*S2","s4*S2"),2,2)
-mod.est <- setModel(drift=alpha, diffusion=beta,
- solve.var=c("S1","S2"),state.variable=c("S1","S2"))
+a <- c("1-mu11*X1+mu12*X2","2+mu21*X1-mu22*X2")
+b <- matrix(c("s1*X1","s2*X1", "-s3*X2","s4*X2"),2,2)
+mod.est <- setModel(drift=a, diffusion=b,
+ solve.var=c("X1","X2"),state.variable=c("X1","X2"))
truep <- list(mu11=.9, mu12=0, mu21=0, mu22=0.7,
s1=.3, s2=0,s3=0,s4=.2)
low <- list(mu11=1e-5, mu12=1e-5, mu21=1e-5, mu22=1e-5,
@@ -573,10 +572,10 @@
fit1 <- qmle(X, start=truep, lower=low, upper=upp, method="L-BFGS-B")
## ------------------------------------------------------------------------
-alpha <- c("1-mu11*S1","2-mu22*S2")
-beta <- matrix(c("s1*S1","0", "0","s4*S2"),2,2)
-mod.est2 <- setModel(drift=alpha, diffusion=beta,
- solve.var=c("S1","S2"),state.variable=c("S1","S2"))
+a <- c("1-mu11*X1","2-mu22*X2")
+b <- matrix(c("s1*X1","0", "0","s4*X2"),2,2)
+mod.est2 <- setModel(drift=a, diffusion=b,
+ solve.var=c("X1","X2"),state.variable=c("X1","X2"))
truep <- list(mu11=.9, mu22=0.7, s1=.3,s4=.2)
low <- list(mu11=1e-5, mu22=1e-5, s1=1e-5, s4=1e-5)
upp <- list(mu11=2, mu22=2, s1=1, s4=1)
@@ -766,7 +765,7 @@
pdf("figures/plot-returns.pdf",width=9,height=4)
par(mar=c(4,4,2,1))
plot(X,main="log returns of AAPL")
-abline(v=cp$tau, lty=3)
+abline(v=cp$tau, lty=3,lwd=2,col="red")
dev.off()
## ----plot-cpoint-aapl,echo=FALSE,results='hide'--------------------------
@@ -827,7 +826,7 @@
p2 <- 0.3
newsamp <- setSampling(random=list(rdist=c(
function(x) rexp(x, rate=p1*n/Terminal),
- function(x) rexp(x, rate=p1*n/Terminal))) )
+ function(x) rexp(x, rate=p2*n/Terminal))) )
## ------------------------------------------------------------------------
Y <- subsampling(X, sampling=newsamp)
@@ -869,7 +868,7 @@
p2 <- 0.3
newsamp <- setSampling(random=list(rdist=c(
function(x) rexp(x, rate=p1*n/Terminal),
- function(x) rexp(x, rate=p1*n/Terminal))) )
+ function(x) rexp(x, rate=p2*n/Terminal))) )
Y <- subsampling(yuima, sampling = newsamp)
## ----cceplot3,fig.keep='none'--------------------------------------------
Modified: pkg/yuima/inst/ybook/chapter4.R
===================================================================
--- pkg/yuima/inst/ybook/chapter4.R 2018-01-16 13:38:47 UTC (rev 635)
+++ pkg/yuima/inst/ybook/chapter4.R 2018-02-02 15:12:53 UTC (rev 636)
@@ -150,7 +150,8 @@
set.seed(127)
x <- rIG(100000,delta,gamma)
hist(x,xlim=c(0,2),ylim=c(0,2),breaks=100,freq=FALSE)
-curve(dIG(x,delta,gamma),add=TRUE,col="red")
+curve(dIG(x,delta,gamma),add=TRUE,col="red",
+ from=min(x), to=max(x), n=500)
mean(x)
var(x)
@@ -158,7 +159,8 @@
pdf("figures/plot-dIG.pdf",width=9,height=4)
par(mar=c(4,4,2,1))
hist(x,xlim=c(0,2),ylim=c(0,2),breaks=100,freq=FALSE)
-curve(dIG(x,delta,gamma),add=TRUE,col="red")
+curve(dIG(x,delta,gamma),add=TRUE,col="red",
+ from=min(x), to=max(x), n=500)
dev.off()
## ----fig.keep='none'-----------------------------------------------------
@@ -198,7 +200,8 @@
## hist(IGsimdata,xlim=c(0,2), ylim=c(0,2), breaks=100, freq=FALSE,
## main=expression(paste("Distribution of ",X[1],
## " and Density of IG(1,2)")))
-## curve(dIG(x,delta,gamma),add=TRUE,col="red")
+## curve(dIG(x,delta,gamma),add=TRUE,col="red",
+## from = 0.001, to = 5, n=500)
## ----plot-IGprocd,echo=FALSE,results='hide'------------------------------
pdf("figures/plot-IGprocd.pdf",width=9,height=4)
@@ -214,7 +217,8 @@
}
hist(IGsimdata,xlim=c(0,2), ylim=c(0,2), breaks=100, freq=FALSE,
main=expression(paste("Distribution of ",X[1]," and Density of IG(1,2)")))
-curve(dIG(x,delta,gamma),add=TRUE,col="red")
+curve(dIG(x,delta,gamma),add=TRUE,col="red",
+ from = 0.001, to = 5, n=500)
dev.off()
## ----fig.keep='none'-----------------------------------------------------
@@ -463,6 +467,151 @@
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