[Yuima-commits] r146 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 8 02:28:53 CEST 2011
Author: hinohide
Date: 2011-04-08 02:28:53 +0200 (Fri, 08 Apr 2011)
New Revision: 146
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/asymptotic_term.R
pkg/yuima/R/simFunctional.R
pkg/yuima/man/asymptotic_term.Rd
Log:
asymptotic_term is extended to multivariate cases
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2011-02-17 16:16:43 UTC (rev 145)
+++ pkg/yuima/DESCRIPTION 2011-04-08 00:28:53 UTC (rev 146)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.184
-Date: 2011-02-17
+Version: 0.1.185
+Date: 2011-04-08
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2011-02-17 16:16:43 UTC (rev 145)
+++ pkg/yuima/NEWS 2011-04-08 00:28:53 UTC (rev 146)
@@ -1,3 +1,6 @@
+Change in Version 0.1.185
+ o asymptotic_term is now possible to handle multivariate cases.
+
Change in Version 0.1.180
o a minor bug on cce is fixed.
Change in Version 0.1.179
Modified: pkg/yuima/R/asymptotic_term.R
===================================================================
--- pkg/yuima/R/asymptotic_term.R 2011-02-17 16:16:43 UTC (rev 145)
+++ pkg/yuima/R/asymptotic_term.R 2011-04-08 00:28:53 UTC (rev 146)
@@ -19,8 +19,10 @@
#e <- gete(yuima at functional)
assign(expand.var, gete(yuima at functional))
- Terminal <- yuima at sampling@Terminal
- division <- yuima at sampling@n
+## Terminal <- yuima at sampling@Terminal
+ Terminal <- yuima at sampling@Terminal[1]
+## division <- yuima at sampling@n
+ division <- yuima at sampling@n[1]
xinit <- getxinit(yuima at functional)
state <- yuima at model@state.variable
V0 <- yuima at model@drift
@@ -29,6 +31,7 @@
d.size <- yuima at model@equation.number
k.size <- length(F)
+ print("compute X.t0")
X.t0 <- Get.X.t0(yuima, expand.var=expand.var)
delta <- deltat(X.t0)
@@ -55,7 +58,7 @@
return(tmp)
}
- # function to solve Y_{t} (between (13.9) and (13.10)) using runge kutta method
+ # function to solve Y_{t} (between (13.9) and (13.10)) using runge kutta method. Y_{t} is GL(d) valued (matrices)
Get.Y <- function(){
## init
dt <- Terminal/division
@@ -82,12 +85,13 @@
for(t in 1:division){
## Xt
for(i in 1:d.size){
- assign(state[i],X.t0[t,i])
+ assign(state[i],X.t0[t,i]) ## state[i] is x_i, for example.
}
## k1
for(i in 1:(d.size*d.size)){
assign(Ystate[i],Yt[i])
}
+
for(i in 1:(d.size*d.size)){
k1[i] <- dt*eval(F[i])
}
@@ -291,7 +295,7 @@
return(mu)
}
- # function to calculate a_{s}^{alpha} in thesis p5
+ # function to calculate a_{s}^{alpha} in bookchapter p5
funca <- function(e=0){
#init
division <- nrow(X.t0)
@@ -307,8 +311,8 @@
# prepare expression of derivatives
for(k in 1:k.size){
- dxf[k] <- deriv(f[[1]][k],state) #expression of derived f0 by x
- dxF[k] <- deriv(F[k],state) #expression of derived F by x
+ dxf[k] <- deriv(f[[1]][k],state) #expression of d f0/dx
+ dxF[k] <- deriv(F[k],state) #expression of d F/dx
def[[k]] <- list()
for(r in 2:(r.size+1)){
def[[k]][[r-1]] <- deriv(f[[r]][k],"e") #expression of derived fa by e
@@ -597,7 +601,7 @@
# if(require(adapt)){
if(require(cubature)){
#tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi0)$value
- tmp <- adaptIntegrate(gz_pi0, lower=rep(min,k.size),upper=rep(max,k.size))$integral
+ tmp <- adaptIntegrate(gz_pi0, lower=rep(min,k.size),upper=rep(max,k.size))$integral
} else {
tmp <- NA
}
@@ -1262,6 +1266,7 @@
get.d1.term<- function(){
lambda.max <- max(eigen(Sigma)$values)
## get g(z)*pi1(z)
+
gz_pi1 <- function(z){
tmp <- g(z) * get.pi1(z)
return( tmp )
@@ -1269,18 +1274,47 @@
## integrate
if( k.size ==1){ # use 'integrate()'
tmp <- integrate(gz_pi1,-Inf,Inf)$value
- }else if( 2 <= k.size || k.size <= 20 ){ # use 'adapt()' to solve multi-dim integration8
+ }else if(2 <= k.size || k.size <= 20 ){ # use sampling approximation for solving multi-dim integration. 2010/11/13
+ set.seed(123)
max <- 10*lambda.max
min <- -10*lambda.max
+ tmp.x <- seq(min,max,length=100)
+ my.x <- NULL
+ for(k in 1:k.size){
+ my.x <- rbind(my.x,tmp.x)
+ }
+ est.points <- my.x[sample(seq(1,nrow(my.x),by=1),nrow(my.x),rep=FALSE),sample(seq(1,ncol(my.x),by=1),ncol(my.x),rep=FALSE)]
+ tmp <- 0
+ for(i in 1:ncol(est.points)){
+ tmp <- tmp + gz_pi1(est.points[,i])
+ }
+ tmp <- tmp/ncol(est.points)
+## }else if( 2 <= k.size || k.size <= 20 ){ # use 'cubatur()' to solve multi-dim integration.
+ }else if( (2 <= k.size || k.size <= 20) && FALSE ){ # use 'cubatur()' to solve multi-dim integration.
+ max <- 10*lambda.max
+## max <- 10*lambda.max/k.size
+ min <- -10*lambda.max
+## min <- -10*lambda.max/k.size
# if(require(adapt)){
if(require(cubature)){
# tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi1)$value
- tmp <- adaptIntegrate(gz_pi1,lower=rep(min,k.size),upper=rep(max,k.size))$integral
+ print("Messages bellow are for debug...")
+ print(date())
+## tmp <- adaptIntegrate(gz_pi1,lower=rep(min,k.size),upper=rep(max,k.size))$integral
+ tmp <- adaptIntegrate(gz_pi1,lower=rep(min,k.size),upper=rep(max,k.size),tol=1e-5,maxEval=500/k.size)
+ print(date())
+ print("tolerance")
+ print(tmp$error)
+ print("function evaluated times")
+ print(tmp$functionEvaluations)
+ print("return code: if it is not 0, error occured in integration (may be does not converge)")
+ print(tmp$returnCode)
+ tmp <- tmp$integral
} else {
tmp <- NA
}
}else{
- stop("length k is too big.")
+ stop("length k is too long.")
}
return(tmp)
}
@@ -1361,10 +1395,9 @@
yuima.warn("Get variables ...")
Y <- Get.Y()
mu <- funcmu()
- aMat <- funca()
+ aMat <- funca() ## ¤³¤³¤Ç¥¨¥é¡¼ 2010/11/24, TBC
Sigma <- funcsigma()
-
# calculate each variables shown in p.9 and
# prepare for trapezium integration
Modified: pkg/yuima/R/simFunctional.R
===================================================================
--- pkg/yuima/R/simFunctional.R 2011-02-17 16:16:43 UTC (rev 145)
+++ pkg/yuima/R/simFunctional.R 2011-04-08 00:28:53 UTC (rev 146)
@@ -4,7 +4,7 @@
##:: fix bug 07/23
assign(expand.var, e)
- division <- yuima at sampling@n #number of observed time
+ division <- yuima at sampling@n[1] #number of observed time
F <- getF(yuima at functional)
d.size <- yuima at model@equation.number
k.size <- length(F)
@@ -26,8 +26,9 @@
funcf <- function(yuima,X,e, expand.var){
##:: fix bug 07/23
assign(expand.var, e)
-
- division <- yuima at sampling@n #number of observed time
+
+ ## division <- yuima at sampling@n #number of observed time
+ division <- yuima at sampling@n[1]
F <- getF(yuima at functional)
f <- getf(yuima at functional)
r.size <- yuima at model@noise.number
@@ -35,6 +36,7 @@
k.size <- length(F)
modelstate <- yuima at model@state.variable
resf <- array(0,dim=c(k.size,division,(r.size+1))) #value of f0,f1,~,fr. size:array[k.size,division,r.size+1]
+
for(r in 1:(r.size+1)){
for(t in 1:division){
@@ -54,42 +56,41 @@
# function to calculate Fe in (13.2).
# core function of 'simFunctional'
funcFe. <- function(yuima,X,e,expand.var="e"){
- ##:: fix bug 07/23
- assign(expand.var, e) ## expand.var
+
+ ##:: fix bug 07/23
+ assign(expand.var, e) ## expand.var
+ F <- getF(yuima at functional)
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ k.size <- length(F)
+ modelstate <- yuima at model@state.variable
+ ## division <- yuima at sampling@n
+ division <- yuima at sampling@n[1] ## 2010/11/13 modified. Currently, division must be the same for each path
+ Terminal <- yuima at sampling@Terminal
+ delta <- Terminal/division #length between observed times
+ dw <- matrix(0,division-1,r.size+1) #Wr(t) size:matrix[division,r.size+1]
+
+ dw[,1] <- rep(Terminal/(division-1),length=division-1) #W0(t)=t
- F <- getF(yuima at functional)
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- k.size <- length(F)
- modelstate <- yuima at model@state.variable
- division <- yuima at sampling@n
- Terminal <- yuima at sampling@Terminal
- delta <- Terminal/division #length between observed times
- dw <- matrix(0,division-1,r.size+1) #Wr(t) size:matrix[division,r.size+1]
-
- dw[,1] <- rep(Terminal/(division-1),length=division-1) #W0(t)=t
-
- for(r in 2:(r.size+1)){
- tmp <- rnorm(division,0,sqrt(delta)) #calculate Wr(t)
- dw[,r] <- tmp[2:division]-tmp[1:(division-1)] #calculate dWr(t)
- }
-
- resF <- funcF(yuima,X,e,expand.var=expand.var) #calculate F with X,e. size:vector[k.size]
- resf <- funcf(yuima,X,e,expand.var=expand.var) #calculate f with X,e. size:array[k.size,division,r.size+1]
-
- Fe <- numeric(k.size)
- for(k in 1:k.size){
- Fe[k] <- sum(resf[k,1:division-1,]*dw)+resF[k] #calculate Fe using resF and resf as (13.2).
- }
- return(Fe)
- }
+ for(r in 2:(r.size+1)){
+ tmp <- rnorm(division,0,sqrt(delta)) #calculate Wr(t)
+ dw[,r] <- tmp[2:division]-tmp[1:(division-1)] #calculate dWr(t)
+ }
+
+ resF <- funcF(yuima,X,e,expand.var=expand.var) #calculate F with X,e. size:vector[k.size]
+ resf <- funcf(yuima,X,e,expand.var=expand.var) #calculate f with X,e. size:array[k.size,division,r.size+1] ## ¤¤¤Þ, ¤³¤³¤Ç¥¨¥é¡¼ 2010/11/13
+ Fe <- numeric(k.size)
+ for(k in 1:k.size){
+ Fe[k] <- sum(resf[k,1:division-1,]*dw)+resF[k] #calculate Fe using resF and resf as (13.2).
+ }
+ return(Fe)
+}
# Get.X.t0
# function to calculate X(t=t0) using runge kutta method
Get.X.t0 <- function(yuima, expand.var="e"){
##:: fix bug 07/23
assign(expand.var,0) ## epsilon=0
-
r.size <- yuima at model@noise.number
d.size <- yuima at model@equation.number
k.size <- length(getF(yuima at functional))
@@ -97,27 +98,71 @@
V0 <- yuima at model@drift
V <- yuima at model@diffusion
- Terminal <- yuima at sampling@Terminal
- division <- yuima at sampling@n
+ Terminal <- yuima at sampling@Terminal[1]
+## division <- yuima at sampling@n
+ division <- yuima at sampling@n[1] ## 2010/11/13
##:: fix bug 07/23
#pars <- #yuima at model@parameter at all[1] #epsilon
-
+
xinit <- getxinit(yuima at functional)
+
## init
dt <- Terminal/division
X <- xinit
Xt <- xinit
+ ## k <- matrix(numeric(d.size^2), nrow=d.size,ncol=d.size)
+## k1 <- matrix(numeric(d.size^2), nrow=d.size,ncol=d.size)
+## k2 <- matrix(numeric(d.size^2), nrow=d.size,ncol=d.size)
+## k3 <- matrix(numeric(d.size^2), nrow=d.size,ncol=d.size)
+## k4 <- matrix(numeric(d.size^2), nrow=d.size,ncol=d.size)
k <- numeric(d.size)
k1 <- numeric(d.size)
k2 <- numeric(d.size)
k3 <- numeric(d.size)
k4 <- numeric(d.size)
-
+
##:: fix bug 07/23
#assign(pars,0) ## epsilon=0
+ ## ## runge kutta
+## for(t in 1:division){
+## ## k1
+## for(i in 1:d.size){
+## assign(modelstate[i],Xt[i])
+## }
+## for(i in 1:d.size){
+## k1[,i] <- dt*eval(V0[i])
+## }
+## ## k2
+## for(i in 1:d.size){
+## assign(modelstate[i],Xt[i]+k1[i]/2)
+## }
+## for(i in 1:d.size){
+## k2[,i] <- dt*eval(V0[i])
+## }
+## ## k3
+## for(i in 1:d.size){
+## assign(modelstate[i],Xt[i]+k2[i]/2)
+## }
+## for(i in 1:d.size){
+## k3[,i] <- dt*eval(V0[i])
+## }
+## ## k4
+## for(i in 1:d.size){
+## assign(modelstate[i],Xt[i]+k3[i])
+## }
+
+## for(i in 1:d.size){
+## k4[,i] <- dt*eval(V0[i])
+## }
+## ## V0(X(t+dt))
+## k <- (k1+k2+k2+k3+k3+k4)/6
+## Xt <- Xt+k
+## X <- rbind(X,Xt)
+## }
+
## runge kutta
for(t in 1:division){
## k1
@@ -153,10 +198,11 @@
Xt <- Xt+k
X <- rbind(X,Xt)
}
+
## return matrix : (division+1)*d
rownames(X) <- NULL
colnames(X) <- modelstate
- return(ts(X,deltat=dt,start=0))
+ return(ts(X,deltat=dt[1],start=0))
}
# simFunctional
@@ -196,7 +242,6 @@
##:: fix bug 07/23
X.t0 <- Get.X.t0(yuima, expand.var=expand.var)
F0 <- funcFe.(yuima, X.t0, 0, expand.var=expand.var)
-
return(F0)
})
Modified: pkg/yuima/man/asymptotic_term.Rd
===================================================================
--- pkg/yuima/man/asymptotic_term.Rd 2011-02-17 16:16:43 UTC (rev 145)
+++ pkg/yuima/man/asymptotic_term.Rd 2011-04-08 00:28:53 UTC (rev 146)
@@ -53,6 +53,38 @@
asymp <- asymptotic_term(yuima, block=10, rho,g)
asymp
sum(asymp$d0 + e * asymp$d1)
+
+
+### An example of multivariate case: Heston model
+## a <- 1;C <- 1;d <- 10;R<-.1
+## diff.matrix <- matrix( c("x1*sqrt(x2)*e", "e*R*sqrt(x2)",0,"sqrt(x2*(1-R^2))*e"), 2,2)
+## model <- setModel(drift = c("a*x1","C*(10-x2)"), diffusion = diff.matrix,solve.variable=c("x1","x2"),state.variable=c("x1","x2"))
+## call option is evaluated by averating
+## max{ (1/T)*int_0^T Xt^e dt, 0}, the first argument is the functional of interest:
+##
+## Terminal <- 1
+## xinit <- c(1,1)
+##
+## f <- list( c(expression(0), expression(0)), c(expression(0), expression(0)) , c(expression(0), expression(0)) )
+## F <- expression(x1,x2)
+##
+## division <- 1000
+## e <- .3
+##
+## yuima <- setYuima(model = model, sampling = setSampling(Terminal=Terminal, n=division))
+## yuima <- setFunctional( yuima, f=f,F=F, xinit=xinit,e=e)
+##
+## rho <- expression(x1)
+## F0 <- F0(yuima)
+## get_ge <- function(x){
+## return( max(x[1],0))
+## }
+## g <- function(x) get_ge(x)
+## set.seed(123)
+## asymp <- asymptotic_term(yuima, block=10, rho,g)
+## sum(asymp$d0 + e * asymp$d1)
+
+
}
More information about the Yuima-commits
mailing list