[Yuima-commits] r197 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 10 16:49:29 CET 2012
Author: kamatani
Date: 2012-01-10 16:49:28 +0100 (Tue, 10 Jan 2012)
New Revision: 197
Added:
pkg/yuima/R/asymptotic_term_third.R
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/asymptotic_term_second.R
Log:
Fixed bugs of asymptotc term
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2011-12-14 14:57:29 UTC (rev 196)
+++ pkg/yuima/DESCRIPTION 2012-01-10 15:49:28 UTC (rev 197)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.1936
+Version: 0.1.194
Date: 2011-09-21
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Modified: pkg/yuima/R/asymptotic_term_second.R
===================================================================
--- pkg/yuima/R/asymptotic_term_second.R 2011-12-14 14:57:29 UTC (rev 196)
+++ pkg/yuima/R/asymptotic_term_second.R 2012-01-10 15:49:28 UTC (rev 197)
@@ -1,1067 +1,1072 @@
-
-yuima.warn <- function(x){
- cat(sprintf("\nYUIMA: %s\n", x))
-}
-
-# in this source we note formulae like latex
-
-
-setGeneric("asymptotic_term",
- function(yuima,block=100, rho, g, expand.var="e")
- standardGeneric("asymptotic_term"))
-
-setMethod("asymptotic_term",signature(yuima="yuima"),
- function(yuima,block=100, rho, g, expand.var="e"){
-
- if(is.null(yuima at model)) stop("model object is missing!")
- if(is.null(yuima at sampling)) stop("sampling object is missing!")
- if(is.null(yuima at functional)) stop("functional object is missing!")
-
- f <- getf(yuima at functional)
- F <- getF(yuima at functional)
-
- ##:: fix bug 07/23
- #e <- gete(yuima at functional)
- assign(expand.var, gete(yuima at functional))
-
- Terminal <- yuima at sampling@Terminal[1]
- division <- yuima at sampling@n[1]
- xinit <- getxinit(yuima at functional)
- state <- yuima at model@state.variable
- V0 <- yuima at model@drift
- V <- yuima at model@diffusion
- r.size <- yuima at model@noise.number
- 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)
-
- ##:: fix bug 07/23
- pars <- expand.var #yuima at model@parameter at all[1] #epsilon
-
- # fix bug 20110628
- if(k.size==1){
- G <- function(x){
- n <- length(x)
- res <- numeric(0)
- for(i in 1:n){
- res <- append(res,g(x[i]))
- }
- return(res)
- }
- }else{
- G <- function(x) return(g(x))
- }
-
-
- # function to return symbolic derivatives of myfunc by mystate(multi-state)
-
- Derivation.vector <- function(myfunc,mystate,dim1,dim2){
-
- tmp <- vector(dim1*dim2,mode="expression")
-
- for(i in 1:dim1){
- for(j in 1:dim2){
- tmp[(i-1)*dim2+j] <- parse(text=deparse(D(myfunc[i],mystate[j])))
- }
- }
-
- return(tmp)
- }
-
- # function to return symbolic derivatives of myfunc by mystate(single state)
-
- Derivation.scalar <- function(myfunc,mystate,dim){
-
- tmp <- vector(dim,mode="expression")
-
- for(i in 1:dim){
- tmp[i] <- parse(text=deparse(D(myfunc[i],mystate)))
- }
- return(tmp)
- }
-
-
- # 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(env){
-
- ## init
- dt <- env$delta
- assign(pars,0) ## epsilon=0
- Yinit <- as.vector(diag(d.size))
- Yt <- Yinit
- Y <- Yinit
- k <- numeric(d.size*d.size)
- k1 <- numeric(d.size*d.size)
- k2 <- numeric(d.size*d.size)
- k3 <- numeric(d.size*d.size)
- k4 <- numeric(d.size*d.size)
- Ystate <- paste("y",1:(d.size*d.size),sep="")
-
- F <- NULL
- F.n <- vector(d.size,mode="expression")
- for(n in 1:d.size){
- for(i in 1:d.size){
- F.tmp <- dx.drift[((i-1)*d.size+1):(i*d.size)]
-
- F.n[i] <- parse(text=paste(paste("(",F.tmp,")","*",
- Ystate[((n-1)*d.size+1):(n*d.size)],sep=""),
- collapse="+"))
- }
- F <- append(F,F.n)
- }
-
- ## runge kutta
- for(t in 1:(division-1)){
-
- ## Xt
- for(i in 1:d.size){
- 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])
- }
-
- ## k2
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k1[i]/2)
- }
-
- for(i in 1:(d.size*d.size)){
- k2[i] <- dt*eval(F[i])
- }
-
- ## k3
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k2[i]/2)
- }
-
- for(i in 1:(d.size*d.size)){
- k3[i] <- dt*eval(F[i])
- }
-
- ## k4
-##modified
- ## Xt+1
- for(i in 1:d.size){
- assign(state[i],X.t0[t+1,i])
- }
-##
-
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k3[i])
- }
-
- for(i in 1:(d.size*d.size)){
- k4[i] <- dt*eval(F[i])
- }
-
- ## F(Y(t+dt))
- k <- (k1+k2+k2+k3+k3+k4)/6
- Yt <- Yt+k
- Y <- rbind(Y,Yt)
- }
-
- ## return matrix : (division+1)*(d.size*d.size)
- rownames(Y) <- NULL
- colnames(Y) <- Ystate
- return(ts(Y,deltat=dt,start=0))
- }
-
-
- #l*t
-
- e_F <- function(X.t0,tmp.F,env){
-
- d.size <- env$d.size
- k.size <- env$k.size
- division <- env$division
-
- result <- matrix(0,k.size,division)
-
- assign(pars[1],0)
-
- de.F <- list()
-
- for(k in 1:k.size){
- de.F[[k]] <- parse(text=deparse(D(tmp.F[k],pars[1])))
- }
-
- for(t in 1:division){
- for(d in 1:d.size){
- assign(state[d],X.t0[t,d])
- }
-
- for(k in 1:k.size){
- result[k,t] <- eval(de.F[[k]])
- }
- }
-
- return(result)
- }
-
-
- # function to calculate mu in thesis p5
- funcmu <- function(e=0,env){
-
- d.size <- env$d.size
- k.size <- env$k.size
- division <- env$division
- delta <- env$delta
-
- n1 <- length(get_e_f0)
- n2 <- sum(get_e_f0 == 0)
-
- n3 <- length(get_e_F)
- n4 <- sum(get_e_F == 0)
-
- n <- (n1 - n2 != 0) * (n3 - n4 != 0)
-
- if(n == 0){
- mu1 <- real(k.size)
- }else{
- mu1 <- apply(get_e_f0,1,sum) * delta + get_e_F[,division]
- }
-
- n5 <- length(get_Y_e_V0)
- n6 <- sum(get_Y_e_V0 == 0)
-
- n <- n5 - n6
-
- if(n == 0){
- mu2 <- real(k.size)
- }else{
- n7 <- length(get_x_F)
- n8 <- sum(get_x_F == 0)
-
- n <- n7 - n8
-
- if(n == 0){
- mu2_1 <- real(k.size)
- }else{
- mu2_1 <- matrix(get_x_F[,,division],k.size,d.size) %*%
- tmpY[,,d.size] %*% matrix(apply(get_Y_e_V0,1,sum),d.size,1) *
- delta
- }
-
- n9 <- length(get_x_f0)
- n10 <- sum(get_x_f0 == 0)
-
- n <- n7 - n8
-
- if(n == 0){
- mu2_2 <- real(k.size)
- }else{
- mu2_2 <- real(k.size)
-
- for(l in 1:k.size){
- for(i in 1:d.size){
- for(j in 1:d.size){
-
- tmp1 <- I0(get_Y_e_V0[j,],env)
- tmp2 <- get_x_f0[l,i,] * tmpY[i,j,] * tmp1
-
- mu2_2[l] <- mu2_2[l] + I0(tmp2,env)
- }
- }
- }
- }
-
- mu2 <- mu2_1 + mu2_2
- }
-
- mu <- mu1 + mu2
-
- return(mu)
- }
-
-
- # function to calculate a_{s}^{alpha} in bookchapter p5
-
- #k*r*t
-
- funca <- function(e=0,env=env){
-
- d.size <- env$d.size
- r.size <- env$r.size
- k.size <- env$k.size
- division <- env$division
- delta <- env$delta
-
- first <- get_e_f
-
- second <- array(0,dim=c(k.size,r.size,division))
-
- second.tmp <- matrix(get_x_F[,,division],k.size,d.size) %*%
- tmpY[,,division]
-
- for(t in 1:division){
- second[,,t] <- second.tmp %*% matrix(get_Y_e_V[,,t],d.size,r.size)
- }
-
- third <- array(0,dim=c(k.size,r.size,division))
-
- n1 <- length(get_x_f0)
- n2 <- sum(get_x_f0 == 0)
-
- n <- n1 - n2
-
- third[,,division] <- 0
-
- if(n != 0){
- for(s in (division-1):1){
- third[,,s] <- third[,,s+1] + matrix(get_x_f0[,,s],k.size,d.size) %*%
- tmpY[,,s] %*% matrix(get_Y_e_V[,,s],d.size,r.size) * delta
- }
- }
-
- result <- first + second + third
-
- return(result) #size:array[k.size,r.size,division]
- }
-
-
- # function to calculate sigma in thesis p5
- # require: aMat
- funcsigma <- function(e=0,env=env){
-
- k.size <- env$k.size
- division <- env$division
- delta <- env$delta
-
- sigma <- matrix(0,k.size,k.size) #size:matrix[k.size,k.size]
-
- for(t in 1:(division-1)){
- a.t <- matrix(aMat[,,t],k.size,r.size)
- sigma <- sigma+(a.t%*%t(a.t)) * delta
- }
-
- # Singularity check
- if(any(eigen(sigma)$value<=0.0001)){
- yuima.warn("Eigen value of covariance matrix in very small.")
- }
-
- return(sigma)
- }
-
-
- ## integrate start:1 end:t number to integrate:block
- # because integration at all 0:T takes too much time,
- # we sample only 'block' points and use trapezium rule to integrate
-
- make.range.for.trapezium.fomula <- function(t,block){
-
- if(t/block <= 1){ # block >= t : just use all points
- range <- c(1:t)
- }else{ # make array which includes points to use
- range <- as.integer( (c(0:block) * (t/block))+1)
-
- if( range[block+1] < t){
- range[block+2] <- t
- }else if( range[block+1] > t){
- range[block+1] <- t
- }
- }
- return(range)
- }
-
-
- ## multi dimension gausian distribusion
-
- normal <- function(x,mu,Sigma){
-
- if(length(x)!=length(mu)){
- print("Error:length of x != one of mu")
- }
-
- dimension <- length(x)
- tmp <- 1/((sqrt(2*pi))^dimension * sqrt(det(Sigma))) * exp(-1/2 * t((x-mu)) %*% solve(Sigma) %*% (x-mu) )
- return(tmp)
- }
-
-
- ## get d0
- ## required library(adapt)
-
- get.d0.term <- function(){
-
- ## get g(z)*pi0(z)
-
-##modified
-# gz_pi0 <- function(z){
-# return( G(z) * H0 *normal(z,mu=mu,Sigma=Sigma))
-# }
-
- gz_pi0 <- function(z){
- tmpz <- matA %*% z
- return(G(tmpz) * H0 * normal(tmpz,mu=mu,Sigma=Sigma)) #det(matA) = 1
- }
-##
-
- gz_pi02 <- function(z){
- return(G(z) * H0 * dnorm(z,mean=mu,sd=sqrt(Sigma)))
- }
-
- ## integrate
- if( k.size ==1){ # use 'integrate' if k=1
-
- tmp <- integrate(gz_pi02,mu-7*sqrt(Sigma),mu+7*sqrt(Sigma))$value
-
- }else if( 2 <= k.size || k.size <= 20 ){ # use library 'adapt' to solve multi-dimentional integration
-
- lambda <- eigen(Sigma)$values
- matA <- eigen(Sigma)$vector
-
- lower <- c()
- upper <- c()
-
- for(k in 1:k.size){
- max <- 7 * sqrt(lambda[k])
- lower[k] <- - max
- upper[k] <- max
- }
-
- if(require(cubature)){
- tmp <- adaptIntegrate(gz_pi0, lower=lower,upper=upper)$integral
- } else {
- tmp <- NA
- }
-
- }else{
- stop("length k is too big.")
- }
- return(tmp)
- }
-
-
- ###############################################################################
- # following funcs are part of d1 term
- # because they are finally integrated at 'get.d1.term()',
- # these funcs are called over and over again.
- # so, we use trapezium rule for integration to save time and memory.
-
- # these funcs almost calculate each formulas including trapezium integration.
- # see each formulas in thesis to know what these funcs do.
-
- # some funcs do alternative calculation at k=1.
- # it depends on 'integrate()' function
- ###############################################################################
-
-
- #a:l, b:l*j*r*t/j*r*t, c:l*r*t
-
- #l
-
- ft1_1 <- function(a1,env){
-
- k.size <- env$k.size
-
- result <- a1
-
- return(result)
- }
-
-
- #first:l*k*k, second:l
-
- ft1_2 <- function(b1,env){
-
- d.size <- env$d.size
- k.size <- env$k.size
- block <- env$block
-
- first <- array(0,dim=c(k.size,k.size,k.size))
- second <- real(k.size)
-
- for(l in 1:k.size){
- for(j in 1:d.size){
- tmp <- I_12(b1[[1]][l,j,,],b1[[2]][j,,],env)
-
- first[l,,] <- first[l,,] + tmp$first[,,block]
- second <- second + tmp$second[block]
- }
- }
-
- return(list(first=first,second=second))
- }
-
-
- #l*k
-
- ft1_3 <- function(c1,env){
-
- k.size <- env$k.size
- block <- env$block
-
- result <- matrix(0,k.size,k.size)
-
- for(l in 1:k.size){
- tmp <- I_1(c1[l,,],env)
-
- result[l,] <- tmp[,block]
- }
-
- return(result)
- }
-
-
- F_tilde1__1 <- function(get_F_tilde1,env){
-
- k.size <- env$k.size
-
- temp <- list()
- temp$first <- array(0,dim=c(k.size,k.size,k.size))
- temp$second <- matrix(0,k.size,k.size)
- temp$third <- real(k.size)
-
-
- calc.range <- c(1:3)
-
- tmp1 <- get_F_tilde1$result1
-
- n1 <- length(tmp1)
- n2 <- sum(tmp1 == 0)
-
- if(n1 == n2){
- calc.range <- calc.range[calc.range != 1]
- }
-
- tmp2 <- get_F_tilde1$result2[[1]]
-
- n3 <- length(tmp2)
- n4 <- sum(tmp2 == 0)
-
- tmp3 <- get_F_tilde1$result2[[2]]
-
- n5 <- length(tmp3)
- n6 <- sum(tmp3 == 0)
-
- n <- (n3 - n4 != 0) * (n5 - n6 != 0)
-
- if(n == 0){
- calc.range <- calc.range[calc.range != 2]
- }
-
- tmp3 <- get_F_tilde1$result3
-
- n7 <- length(tmp3)
- n8 <- sum(tmp3 == 0)
-
- if(n7 == n8){
- calc.range <- calc.range[calc.range != 3]
- }
-
-
- for(i in calc.range){
-
- tmp <- switch(i,"a","b","c")
-
- result <- switch(tmp,"a"=ft1_1(get_F_tilde1[[i]],env),
- "b"=ft1_2(get_F_tilde1[[i]],env),
- "c"=ft1_3(get_F_tilde1[[i]],env))
-
- nlist <- length(result)
-
- if(nlist != 2){
-
- tmp1 <- length(dim(result)) #2 or NULL
-
- if(tmp1 == 0){
- tmp1 <- 3
- }
-
- temp[[tmp1]] <- temp[[tmp1]] + result
-
- }else{
-
- temp[[1]] <- temp[[1]] + result[[1]]
-
- temp[[3]] <- temp[[3]] + result[[2]]
-
- }
- }
-
- first <- temp$first
- second <- temp$second
- third <- temp$third
-
- return(list(first=first,second=second,third=third))
- }
-
-
- F_tilde1_x <- function(l,x){
-
- first <- matrix(get_F_tilde1__1$first[l,,],
- k.size,k.size)
-
- result1 <- 0
-
- for(k1 in 1:k.size){
- for(k2 in 1:k.size){
- result1 <- result1 + first[k1,k2] * x[k1] * x[k2]
- }
- }
-
- second <- get_F_tilde1__1$second[l,]
-
- result2 <- 0
-
- for(k1 in 1:k.size){
- result2 <- result2 + second[k1] * x[k1]
- }
-
- result3 <- get_F_tilde1__1$third[l]
-
- result <- result1 + result2 + result3
-
- return(result)
- }
-
-
- #h:d*block
-
- #first:numeric(1), second:k.size*block
-
- Di_bar <- function(h,env){
-
- block <- env$block
-
- first <- real(1)
-
- tmp4 <- matrix(0,r.size,block)
- tmp6 <- matrix(0,r.size,block)
-
- tmp5 <- matrix(0,r.size,block)
-
- for(i in 1:d.size){
- tmp1 <- h[i,] * get_Y_D[i,]
- first <- first + I0(tmp1,env)[block]
-
- for(j in 1:d.size){
- tmp2 <- h[i,] * tmpY[i,j,]
- tmp3 <- I0(tmp2,env)
- tmp4 <- tmp4 + tmp3[block] * get_Y_e_V[j,,]
-
- for(r in 1:r.size){
- tmp5[r,] <- tmp3 * get_Y_e_V[j,r,]
- }
- tmp6 <- tmp6 + tmp5
- }
- }
-
- tmp7 <- tmp4 - tmp6
- second <- I_1(tmp7,env)
-
- return(list(first=first,second=second))
- }
-
-
- Di_bar_x <- function(h,x){
-
- tmp1 <- Di_bar(h,env)
-
- tmp2 <- tmp1$second
-
- result <- tmp1$first + I_1_x(x,tmp2,env)
-
- return(result)
- }
-
-
- get.P2 <- function(z){
-
- block <- env$block
-
- if(k.size==1){
- First <- Di_bar_x(di.rho,z)
- Second <- rep(de.rho %*% Diff[block,] * delta ,length(z))
- }else{
- First <- Di_bar_x(di.rho,z)
- Second <- de.rho %*% Diff[block,] * delta
- }
-
- tmp <- First + Second
-
- return(tmp)
- }
-
-
- # now calculate pi1 using funcs above
-
- get.pi1 <- function(z){
-
- First <- 0
- z.tilda <- z - mu
-
- if(k.size ==1){
- zlen <- length(z)
- result <- c()
-
- for(i in 1:zlen){
- First <- (F_tilde1_x(1,z.tilda[i]+delta) *
- dnorm(z[i]+delta,mean=mu,sd=sqrt(Sigma)) -
- F_tilde1_x(1,z.tilda[i]) *
- dnorm(z[i],mean=mu,sd=sqrt(Sigma)))/delta
-
- Second <- get.P2(z.tilda[i]) * dnorm(z[i],mean=mu,sd=sqrt(Sigma))
-
- result[i] <- First + Second
- }
- }else{
- tmp <- numeric(k.size)
-
- for(k in 1:k.size){
- dif <- numeric(k.size)
- dif[k] <- dif[k] + delta
- z.delta <- z.tilda + dif
- tmp[k] <- (F_tilde1_x(k,z.delta) * normal(z.delta,numeric(k.size),Sigma) -
- F_tilde1_x(k,z.tilda) * normal(z.tilda,numeric(k.size),Sigma))/delta
- }
-
- First <- sum(tmp)
-
- Second <- get.P2(z.tilda) * normal(z.tilda,numeric(k.size),Sigma)
-
- result <- First + Second
- }
-
- pi1 <- - H0 * result
- return(pi1)
- }
-
- # calculate d1
- ## required library(adapt)
-
- get.d1.term<- function(){
-
- ## get g(z)*pi1(z)
-
- gz_pi1 <- function(z){
- tmp <- G(z) * get.pi1(z)
- return( tmp )
- }
-
- ## integrate
- if(k.size == 1){ # use 'integrate()'
- tmp <- integrate(gz_pi1,mu-7*sqrt(Sigma),mu+7*sqrt(Sigma))$value
-
- }else if(2 <= k.size || k.size <= 20){
-
- lambda <- eigen(Sigma)$values
-
- matA <- eigen(Sigma)$vector
-
- gz_pi1 <- function(z){
- tmpz <- matA %*% z
- tmp <- G(tmpz) * get.pi1(tmpz) #det(matA) = 1
- return( tmp )
- }
-
- my.x <- matrix(0,k.size,20^k.size)
- dt <- 1
-
- for(k in 1:k.size){
- max <- 7 * sqrt(lambda[k])
- min <- -7 * sqrt(lambda[k])
- tmp.x <- seq(min,max,length=20)
- dt <- dt * (tmp.x[2] - tmp.x[1])
- my.x[k,] <- rep(tmp.x,each=20^(k.size-k),times=20^(k-1))
- }
-
- tmp <- 0
-
- for(i in 1:20^k.size){
- tmp <- tmp + gz_pi1(my.x[,i])
- }
-
- tmp <- tmp * dt
-
- }else{
- stop("length k is too long.")
- }
- return(tmp)
- }
-
-
- # 'rho' is a given function of (X,e) (p.7 formula(13.19))
-
- # d(rho)/di
- get.di.rho <- function(){
-
- assign(pars[1],0)
- di.rho <- numeric(d.size)
- tmp <- matrix(0,d.size,block+1)
-
- for(i in 1:d.size){
- di.rho[i] <- deriv(rho,state[i])
- }
-
- for(t in 1:(block+1)){
- for(i in 1:d.size){
- assign(state[i],X.t0[t,i])
- }
- for(j in 1:d.size){
- tmp[j,t]<- attr(eval(di.rho[j]),"grad")
- }
- }
- return(tmp)
- }
-
-
- # d(rho)/de
- get.de.rho <- function(){
-
- assign(pars[1],0)
- tmp <- matrix(0,1,block+1)
- de.rho <- deriv(rho,pars[1])
-
- for(t in 1:(block+1)){
- for(i in 1:d.size){
- assign(state[i],X.t0[t,i])
- }
- tmp[1,t]<- attr(eval(de.rho),"grad")
- }
- return(tmp)
- }
-
-
- # H_{t}^{e} at t=0 (see formula (13.19))
- # use trapezium integration
- get.H0 <- function(){
-
- assign(pars[1],0)
- tmp <- matrix(0,1,division-1)
-
- for(t in 1:(division-1)){
- for(i in 1:d.size){
- assign(state[i],X.t0[t,i])
- }
- tmp[1,t] <- eval(rho)
- }
- H0 <- exp(-(sum(tmp) * delta))
- return(as.double(H0) )
- }
-
-
- #################################################
- # Here is an entry point of 'asymptotic_term()' #
- #################################################
-
- ## initialization part
-
- division <- nrow(X.t0)
-# delta <- T/(division - 1)
-
- # make expressions of derivation of V0
- dx.drift <- Derivation.vector(V0,state,d.size,d.size)
- de.drift <- Derivation.scalar(V0,pars,d.size)
- dede.drift <- Derivation.scalar(de.drift,pars,d.size)
-
- # make expressions of derivation of V
- dx.diffusion <- as.list(NULL)
- for(i in 1:d.size){
- dx.diffusion[i] <- list(Derivation.vector(V[[i]],state,d.size,r.size))
- }
-
- de.diffusion <- as.list(NULL)
- for(i in 1:d.size){
- de.diffusion[i] <- list(Derivation.scalar(V[[i]],pars,r.size))
- }
-
- dede.diffusion <- as.list(NULL)
- for(i in 1:d.size){
- dede.diffusion[i] <- list(Derivation.scalar(de.diffusion[[i]],pars,r.size))
- }
- dxdx.drift <- list()
- dxdx.diffusion <- list()
-
- for(i1 in 1:d.size){
- dxdx.drift[[i1]] <- list()
- dxdx.diffusion[[i1]] <- list()
-
- for(i2 in 1:d.size){
- dxdx.drift[[i1]][[i2]] <- list()
- dxdx.diffusion[[i1]][[i2]] <- list()
-
- for(i in 1:d.size){
- tmp1 <- parse(text=deparse(D(V0[i],state[i2])))
- dxdx.drift[[i1]][[i2]][i] <- parse(text=deparse(D(tmp1,state[i1])))
- dxdx.diffusion[[i1]][[i2]][[i]] <- list()
-
- for(r in 1:r.size){
- tmp2 <- parse(text=deparse(D(V[[i]][r],state[i2])))
- dxdx.diffusion[[i1]][[i2]][[i]][r] <- parse(text=deparse(D(tmp2,state[i1])))
- }
- }
- }
- }
-
- dxde.drift <- list()
- dxde.diffusion <- list()
-
- for(i1 in 1:d.size){
- dxde.drift[[i1]] <- list()
- dxde.diffusion[[i1]] <- list()
-
- for(i in 1:d.size){
- tmp1 <- parse(text=deparse(D(V0[i],pars[1])))
- dxde.drift[[i1]][i] <- parse(text=deparse(D(tmp1,state[i1])))
- dxde.diffusion[[i1]][[i]] <- list()
-
- for(r in 1:r.size){
- tmp2 <- parse(text=deparse(D(V[[i]][r],pars[1])))
- dxde.diffusion[[i1]][[i]][r] <- parse(text=deparse(D(tmp2,state[i1])))
- }
- }
- }
-
-
- env <- new.env()
- env$d.size <- d.size
- env$r.size <- r.size
- env$k.size <- k.size
- env$division <- division
- env$delta <- delta
- env$state <- state
- env$pars <- pars
- env$block <- division
- env$my.range <- c(1:division)
-
- yuima.warn("Get variables ...")
- Y <- Get.Y(env=env)
-
- tmpY <- array(0,dim=c(d.size,d.size,division))
- invY <- array(0,dim=c(d.size,d.size,division))
-
- for(t in 1:division){
- tmpY[,,t] <- matrix(Y[t,],d.size,d.size)
- invY[,,t] <- solve(tmpY[,,t])
- }
-
- env$invY <- invY
-
- get_e_f0 <- e_f0(X.t0,f,env)
- get_e_f <- e_f(X.t0,f,env)
- get_x_f0 <- x_f0(X.t0,f,env)
- get_e_F <- e_F(X.t0,F,env)
- get_x_F <- x_F(X.t0,F,env)
-
- get_Y_e_V0 <- Y_e_V0(X.t0,de.drift,env)
- get_Y_e_V <- Y_e_V(X.t0,de.diffusion,env)
-
- mu <- funcmu(env=env)
- aMat <- funca(env=env) ## ¤³¤³¤Ç¥¨¥é¡¼ 2010/11/24, TBC
- Sigma <- funcsigma(env=env)
-
- invSigma <- solve(Sigma)
-
- di.rho <- get.di.rho()
- de.rho <- get.de.rho()
- H0 <- get.H0()
-
- yuima.warn("Done.")
-
- ## calculation
- yuima.warn("Calculating d0 ...")
- d0 <- get.d0.term()
- yuima.warn("Done\n")
-
-
- yuima.warn("Calculating d1 term ...")
-
- # prepare for trapezium integration
-
- my.range <- make.range.for.trapezium.fomula(division,block)
- my.diff <- my.range[2:(block+1)] - my.range[1:block]
-
- Diff <- matrix(0,block+1,block+1)
-
- for(t in 2:length(my.range)){
- for(s in 1:t){
- if( s==1 || t==s ){
- if(s==1){
- Diff[t,s] <- my.diff[s]
- }else if(t==s){
- Diff[t,s] <- my.diff[s-1]
- }
- }else{
- Diff[t,s] <- my.diff[s-1] + my.diff[s]
- }
- }
- }
-
- Diff <- Diff / 2 # trapezium integrations need "/2"
-
-
- env$aMat <- aMat
- env$mu <- mu
- env$Sigma <- Sigma
- env$invSigma <- invSigma
- env$invY <- array(invY[,,my.range],dim=c(d.size,d.size,block+1))
- env$block <- block + 1
- env$my.range <- my.range
- env$Diff <- Diff
-
- tmpY <- array(tmpY[,,my.range],dim=c(d.size,d.size,block+1))
-
- get_Y_e_V0 <- Y_e_V0(X.t0,de.drift,env)
- get_Y_e_V <- Y_e_V(X.t0,de.diffusion,env)
- get_Y_D <- Y_D(X.t0,tmpY,get_Y_e_V0,env)
- get_Y_x1_x2_V0 <- Y_x1_x2_V0(X.t0,dxdx.drift,env)
- get_Y_x1_x2_V <- Y_x1_x2_V(X.t0,dxdx.diffusion,env)
- get_Y_x_e_V0 <- Y_x_e_V0(X.t0,dxde.drift,env)
- get_Y_x_e_V <- Y_x_e_V (X.t0,dxde.diffusion,env)
- get_Y_e_e_V0 <- Y_e_e_V0(X.t0,dede.drift,env)
- get_Y_e_e_V <- Y_e_e_V(X.t0,dede.diffusion,env)
-
- get_D0_t <- D0_t(tmpY, get_Y_D, get_Y_e_V)
- get_e_t <- e_t(tmpY, get_Y_e_V, get_Y_D, get_Y_x1_x2_V0, get_Y_x_e_V0, get_Y_e_e_V0, env)
- get_U_t <- U_t(tmpY, get_Y_D, get_Y_x1_x2_V0, get_Y_x_e_V0, env)
- get_U_hat_t <- U_hat_t(tmpY, get_Y_e_V, get_Y_D, get_Y_x1_x2_V0, get_Y_x_e_V, env)
- get_E0_t <- E0_t(tmpY, get_Y_e_V, get_Y_x1_x2_V0, get_Y_e_e_V, get_e_t, get_U_t, get_U_hat_t, env)
-
- get_e_f0 <- e_f0(X.t0,f,env)
- get_e_f <- e_f(X.t0,f,env)
- get_x_f0 <- x_f0(X.t0,f,env)
- get_x1_x2_f0 <- x1_x2_f0(X.t0,f,env)
- get_x1_x2_f <- x1_x2_f(X.t0,f,env)
- get_x_e_f0 <- x_e_f0(X.t0,f,env)
- get_x_e_f <- x_e_f(X.t0,f,env)
- get_e_e_f0 <- e_e_f0(X.t0,f,env)
- get_e_e_f <- e_e_f(X.t0,f,env)
-
- get_F_t <- F_t (tmpY, get_Y_e_V, get_Y_D, get_x1_x2_f0, get_x_e_f0, get_e_e_f0, env)
- get_W_t <- W_t (tmpY, get_Y_D, get_x1_x2_f0, get_x_e_f0, env)
- get_W_hat_t <- W_hat_t (tmpY, get_Y_e_V, get_x1_x2_f0, get_x_e_f, env)
-
- get_F_tilde1_1 <- F_tilde1_1(tmpY, get_Y_e_V, get_x1_x2_f0, get_e_e_f, get_F_t, get_W_t, get_W_hat_t, env)
- get_F_tilde1_2 <- F_tilde1_2(get_E0_t,get_x_f0,env)
-
- get_x_F <- x_F(X.t0,F,env)
- get_x1_x2_F <- x1_x2_F(X.t0,F,env)
- get_x_e_F <- x_e_F(X.t0,F,env)
- get_e_e_F <- e_e_F(X.t0,F,env)
-
- get_F_tilde1_3 <- F_tilde1_3(get_E0_t,get_x_F,env)
- get_F_tilde1_4 <- F_tilde1_4(tmpY, get_Y_e_V, get_Y_D, get_x_F, get_x1_x2_F, get_x_e_F, get_e_e_F, env)
- get_F_tilde1 <- F_tilde1(get_F_tilde1_1, get_F_tilde1_2, get_F_tilde1_3, get_F_tilde1_4)
-
- get_F_tilde1__1 <- F_tilde1__1(get_F_tilde1,env)
-
- d1 <- get.d1.term()
- yuima.warn("Done\n")
-
- terms <- list(d0=d0, d1=d1)
- return(terms)
-})
-
-
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 197
More information about the Yuima-commits
mailing list