[Yuima-commits] r159 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 22 09:13:48 CEST 2011
Author: kamatani
Date: 2011-07-22 09:13:48 +0200 (Fri, 22 Jul 2011)
New Revision: 159
Added:
pkg/yuima/R/asymptotic_term_third_function.R
Log:
Add third order expansion
Added: pkg/yuima/R/asymptotic_term_third_function.R
===================================================================
--- pkg/yuima/R/asymptotic_term_third_function.R (rev 0)
+++ pkg/yuima/R/asymptotic_term_third_function.R 2011-07-22 07:13:48 UTC (rev 159)
@@ -0,0 +1,9671 @@
+
+ ###########################
+ # third expansion #
+ ###########################
+
+ #bi:r.size*division #aMat:k.size*r.size*division
+ #Sigma:k.size*k.size
+
+ #env:d.size,r.size,k.size,division,delta,state,pars,aMat,mu,Sigma,
+ # invSigma,invY,block,my.range,Diff
+
+ #first:k*k*k*t, second:k*t
+
+ I_123 <- function(b1,b2,b3,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ delta <- env$delta
+ aMat <- env$aMat
+ invSigma <- env$invSigma
+ block <- env$block
+ my.range <- env$my.range
+ Diff <- env$Diff
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) * (n5 - n6 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,block))
+ second <- matrix(0,k.size,block)
+
+ return(list(first=first,second=second))
+ }
+
+ b1a <- matrix(0,k.size,block)
+ b2a <- matrix(0,k.size,block)
+ b3a <- matrix(0,k.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ b2 <- matrix(b2,1,block)
+ b3 <- matrix(b3,1,block)
+ }
+
+ for(t in 1:block){
+ b1a[,t] <- matrix(b1[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b2a[,t] <- matrix(b2[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b3a[,t] <- matrix(b3[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ }
+
+ first <- array(0,dim=c(k.size,k.size,k.size,block))
+
+ tmp1 <- matrix(0,k.size,block)
+ tmp2 <- array(0,dim=c(k.size,k.size,block))
+ tmp3 <- array(0,dim=c(k.size,k.size,k.size,block))
+
+ for(t in 2:block){
+ tmp1[,t] <- b3a %*% Diff[t,] * delta
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(t in 2:block){
+ tmp2[k1,k2,t] <- (b2a[k1,] * tmp1[k2,]) %*%
+ Diff[t,] * delta
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(t in 2:block){
+ tmp3[k1,k2,k3,t] <- (b1a[k1,] * tmp2[k2,k3,]) %*%
+ Diff[t,] * delta
+ }
+ }
+ }
+ }
+
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(t in 1:block){
+ first[k1,k2,,t] <- matrix(tmp3[k1,k2,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ first[k1,,k3,t] <- matrix(first[k1,,k3,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ first[,k2,k3,t] <- matrix(first[,k2,k3,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+
+
+ second1 <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ second1[,t] <- second1[,t] + tmp3[k1,k2,,t] * invSigma[k1,k2]
+ }
+ }
+ second1[,t] <- matrix(second1[,t],1,k.size) %*% invSigma
+ }
+
+ second2 <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k3 in 1:k.size){
+ second2[,t] <- second2[,t] + tmp3[k1,,k3,t] * invSigma[k1,k3]
+ }
+ }
+ second2[,t] <- matrix(second2[,t],1,k.size) %*% invSigma
+ }
+
+ second3 <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ second3[,t] <- second3[,t] + tmp3[,k2,k3,t] * invSigma[k2,k3]
+ }
+ }
+ second3[,t] <- matrix(second3[,t],1,k.size) %*% invSigma
+ }
+
+ second <- - second1 - second2 - second3
+ return(list(first=first,second=second))
+ }
+
+
+ I_123_x <- function(x,get_I_123,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ first <- array(get_I_123$first[,,,block],dim=c(k.size,k.size,k.size))
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ result1 <- result1 + first[k1,k2,k3] *
+ x[k1] * x[k2] * x[k3]
+ }
+ }
+ }
+
+ second <- get_I_123$second[,block]
+
+ result2 <- matrix(x,1,k.size) %*%
+ matrix(second,k.size,1)
+
+ result <- result1 + result2
+
+ return(result)
+ }
+
+
+ #first:k*k*k*k*t, second:k*k*t, third:t
+
+ I_1234 <- function(b1,b2,b3,b4,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ delta <- env$delta
+ aMat <- env$aMat
+ invSigma <- env$invSigma
+ block <- env$block
+ my.range <- env$my.range
+ Diff <- env$Diff
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n7 <- length(b4)
+ n8 <- sum(b4 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) *
+ (n5 - n6 != 0) * (n7 - n8 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,k.size,block))
+ second <- array(0,dim=c(k.size,k.size,block))
+ third <- real(block)
+
+ return(list(first=first,second=second,third=third))
+ }
+
+ b1a <- matrix(0,k.size,block)
+ b2a <- matrix(0,k.size,block)
+ b3a <- matrix(0,k.size,block)
+ b4a <- matrix(0,k.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ b2 <- matrix(b2,1,block)
+ b3 <- matrix(b3,1,block)
+ b4 <- matrix(b4,1,block)
+ }
+
+ for(t in 1:block){
+ b1a[,t] <- matrix(b1[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b2a[,t] <- matrix(b2[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b3a[,t] <- matrix(b3[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b4a[,t] <- matrix(b4[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ }
+
+ first <- array(0,dim=c(k.size,k.size,k.size,k.size,block))
+
+ tmp1 <- matrix(0,k.size,block)
+ tmp2 <- array(0,dim=c(k.size,k.size,block))
+ tmp3 <- array(0,dim=c(k.size,k.size,k.size,block))
+ tmp4 <- array(0,dim=c(k.size,k.size,k.size,k.size,block))
+
+ for(t in 2:block){
+ tmp1[,t] <- b4a %*% Diff[t,] * delta
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(t in 2:block){
+ tmp2[k1,k2,t] <- (b3a[k1,] * tmp1[k2,]) %*%
+ Diff[t,] * delta
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(t in 2:block){
+ tmp3[k1,k2,k3,t] <- (b2a[k1,] * tmp2[k2,k3,]) %*%
+ Diff[t,] * delta
+ }
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ for(t in 2:block){
+ tmp4[k1,k2,k3,k4,t] <- (b1a[k1,] * tmp3[k2,k3,k4,]) %*%
+ Diff[t,]* delta
+ }
+ }
+ }
+ }
+ }
+
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ first[k1,k2,k3,,t] <- matrix(tmp4[k1,k2,k3,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ first[k1,k2,,k4,t] <- matrix(first[k1,k2,,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ first[k1,,k3,k4,t] <- matrix(first[k1,,k3,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ first[,k2,k3,k4,t] <- matrix(first[,k2,k3,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+ }
+ }
+
+
+ second1 <- array(0,dim=c(k.size,k.size,block))
+
+ tmp5 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ tmp5[,,t] <- tmp5[,,t] + tmp4[k1,k2,,,t] *
+ invSigma[k1,k2]
+ }
+ }
+ }
+
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ second1[k3,,t] <- matrix(tmp5[k3,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ second1[,k4,t] <- matrix(second1[,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second2 <- array(0,dim=c(k.size,k.size,block))
+
+ tmp6 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k3 in 1:k.size){
+ tmp6[,,t] <- tmp6[,,t] + tmp4[k1,,k3,,t] *
+ invSigma[k1,k3]
+ }
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(t in 1:block){
+ second2[k2,,t] <- matrix(tmp6[k2,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ second2[,k4,t] <- matrix(second2[,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second3 <- array(0,dim=c(k.size,k.size,block))
+
+ tmp7 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k4 in 1:k.size){
+ tmp7[,,t] <- tmp7[,,t] + tmp4[k1,,,k4,t] *
+ invSigma[k1,k4]
+ }
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(t in 1:block){
+ second3[k2,,t] <- matrix(tmp7[k2,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ second3[,k3,t] <- matrix(second3[,k3,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second4 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ second4[,,t] <- second4[,,t] + tmp4[,k2,k3,,t] *
+ invSigma[k2,k3]
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(t in 1:block){
+ second4[k1,,t] <- matrix(second4[k1,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k4 in 1:k.size){
+ for(t in 1:block){
+ second4[,k4,t] <- matrix(second4[,k4,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second5 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k2 in 1:k.size){
+ for(k4 in 1:k.size){
+ second5[,,t] <- second5[,,t] + tmp4[,k2,,k4,t] *
+ invSigma[k2,k4]
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(t in 1:block){
+ second5[k1,,t] <- matrix(second5[k1,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k3 in 1:k.size){
+ for(t in 1:block){
+ second5[,k3,t] <- matrix(second5[,k3,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second6 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 1:block){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ second6[,,t] <- second6[,,t] + tmp4[,,k3,k4,t] *
+ invSigma[k3,k4]
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(t in 1:block){
+ second6[k1,,t] <- matrix(second6[k1,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(t in 1:block){
+ second6[,k2,t] <- matrix(second6[,k2,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ second <- - second1 - second2 - second3 - second4 - second5 - second6
+
+
+ third1 <- real(block)
+
+ for(t in 1:block){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ third1[t] <- third1[t] + tmp5[k3,k4,t] *
+ invSigma[k3,k4]
+ }
+ }
+ }
+
+
+ third2 <- real(block)
+
+ for(t in 1:block){
+ for(k2 in 1:k.size){
+ for(k4 in 1:k.size){
+ third2[t] <- third2[t] + tmp6[k2,k4,t] *
+ invSigma[k2,k4]
+ }
+ }
+ }
+
+ third3 <- real(block)
+
+
+ d.size <- env$d.size
+ r.size <- env$r.size
+ k.size <- env$k.size
+ state <- env$state
+ pars <- env$pars
+ block <- env$block
+ my.range <- env$my.range
+
+ for(t in 1:block){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ third3[t] <- third3[t] + tmp6[k2,k3,t] *
+ invSigma[k2,k3]
+ }
+ }
+ }
+
+ third <- third1 + third2 + third3
+
+ return(list(first=first,second=second,third=third))
+ }
+
+
+ I_1234_x <- function(x,get_I_1234,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ first <- array(get_I_1234$first[,,,,block],
+ dim=c(k.size,k.size,k.size,k.size))
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(k3 in 1:k.size){
+ for(k4 in 1:k.size){
+ result1 <- result1 + first[k1,k2,k3,k4] *
+ x[k1] * x[k2] * x[k3] * x[k4]
+ }
+ }
+ }
+ }
+
+ second <- matrix(get_I_1234$second[,,block],k.size,k.size)
+
+ result2 <- t(matrix(x,k.size,1)) %*% second %*%
+ matrix(x,k.size,1)
+
+ third <- get_I_1234$third[block]
+
+ result <- result1 + result2 + third
+
+ return(result)
+ }
+
+
+ b1b2 <- function(b1,b2,env){
+
+ r.size <- env$r.size
+ block <- env$block
+
+ result <- real(block)
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0)
+
+ if(n == 0){
+ return(result)
+ }
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ b2 <- matrix(b2,1,block)
+ }
+
+ for(t in 1:block){
+ result[t] <- matrix(b1[,t],1,r.size) %*%
+ matrix(b2[,t],r.size,1)
+ }
+
+ return(result)
+ }
+
+
+ b1_b2 <- function(b1,b2,env){
+
+ delta <- env$delta
+ block <- env$block
+ Diff <- env$Diff
+
+ tmp <- b1b2(b1,b2,env)
+
+ result <- real(block)
+
+ for(t in 1:block){
+ result[t] <- tmp %*% Diff[t,] * delta
+ }
+
+ return(result)
+ }
+
+
+ I0 <- function(f,env){
+
+ delta <- env$delta
+ block <- env$block
+ Diff <- env$Diff
+
+ result <- real(block)
+
+ n1 <- length(f)
+ n2 <- sum(f == 0)
+
+ n <- n1 - n2
+
+ if(n == 0){
+ return(result)
+ }
+
+ for(t in 1:block){
+ result[t] <- f %*% Diff[t,] * delta
+ }
+
+ return(result)
+ }
+
+
+ #k*t
+
+ I_1 <- function(b1,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ delta <- env$delta
+ aMat <- env$aMat
+ invSigma <- env$invSigma
+ block <- env$block
+ my.range <- env$my.range
+ Diff <- env$Diff
+
+ result <- matrix(0,k.size,block)
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ if(n1 == n2){
+ return(result)
+ }
+
+ b1a <- matrix(0,k.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ }
+
+ for(t in 1:block){
+ b1a[,t] <- matrix(b1[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ }
+
+ for(t in 2:block){
+ result[,t] <- t(invSigma) %*% b1a %*% Diff[t,] * delta
+ }
+
+ return(result)
+ }
+
+
+ I_1_x <- function(x,get_I_1,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ tmp <- get_I_1[,block]
+
+ result <- matrix(tmp,1,k.size) %*% matrix(x,k.size,1)
+
+ return(result)
+ }
+
+
+ #first:k*k*t, second:t
+
+ I_12 <- function(b1,b2,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ delta <- env$delta
+ aMat <- env$aMat
+ Sigma <- env$Sigma
+ invSigma <- env$invSigma
+ block <- env$block
+ my.range <- env$my.range
+ Diff <- env$Diff
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,block))
+ second <- real(block)
+
+ return(list(first=first,second=second))
+ }
+
+ b1a <- matrix(0,k.size,block)
+ b2a <- matrix(0,k.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ b2 <- matrix(b2,1,block)
+ }
+
+ for(t in 1:block){
+ b1a[,t] <- matrix(b1[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ b2a[,t] <- matrix(b2[,t],1,r.size) %*%
+ t(matrix(aMat[,,my.range[t]],k.size,r.size))
+ }
+
+ first <- array(0,dim=c(k.size,k.size,block))
+
+ tmp1 <- matrix(0,k.size,block)
+ tmp2 <- array(0,dim=c(k.size,k.size,block))
+
+ for(t in 2:block){
+ tmp1[,t] <- b2a %*% Diff[t,] * delta
+ }
+
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ for(t in 2:block){
+ tmp2[k1,k2,t] <- (b1a[k1,] * tmp1[k2,]) %*%
+ Diff[t,] * delta
+ }
+ }
+ }
+
+ for(k1 in 1:k.size){
+ for(t in 1:block){
+ first[k1,,t] <- matrix(tmp2[k1,,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+ for(k2 in 1:k.size){
+ for(t in 1:block){
+ first[,k2,t] <- matrix(first[,k2,t],1,k.size) %*%
+ invSigma
+ }
+ }
+
+
+ second <- real(block)
+
+ for(t in 1:block){
+ for(k1 in 1:k.size){
+ for(k2 in 1:k.size){
+ second[t] <- second[t] + tmp2[k1,k2,t] * invSigma[k1,k2]
+ }
+ }
+ }
+
+ second <- - second
+
+ return(list(first=first,second=second))
+ }
+
+
+ I_12_x <- function(x,get_I_12,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ first <- matrix(get_I_12$first[,,block],k.size,k.size)
+ second <- get_I_12$second[block]
+
+ result <- matrix(x,1,k.size) %*% first %*%
+ matrix(x,k.size,1) + second
+
+ return(result)
+ }
+
+
+ #first:k*k*t, second:t
+
+ I_1_2 <- function(b1,b2,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,block))
+ second <- real(block)
+
+ return(list(first=first,second=second))
+ }
+
+ first.tmp <- I_12(b1,b2,env)
+ second.tmp <- I_12(b2,b1,env)
+ third.tmp <- b1_b2(b1,b2,env)
+
+ first <- first.tmp$first + second.tmp$first
+ second <- first.tmp$second + second.tmp$second + third.tmp
+
+ return(list(first=first,second=second))
+ }
+
+
+ I_1_2_x <- function(x,get_I_1_2,env){
+
+ result <- I_12_x(x,get_I_1_2,env)
+
+ return(result)
+ }
+
+
+ #first:k*k*t, second:t
+
+ I_1_23 <- function(b1,b2,b3,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ block <- env$block
+ my.range <- env$my.range
+ Diff <- env$Diff
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) * (n5 - n6 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,block))
+ second <- matrix(0,k.size,block)
+
+ return(list(first=first,second=second))
+ }
+
+ first.tmp <- I_123(b1,b2,b3,env)
+ second.tmp <- I_123(b2,b1,b3,env)
+ third.tmp <- I_123(b2,b3,b1,env)
+
+ tmp1 <- b1_b2(b1,b3,env)
+ tmp2 <- matrix(0,r.size,block)
+
+ if(r.size == 1){
+ b2 <- matrix(b2,1,block)
+ }
+
+ for(t in 1:block){
+ tmp2[,t] <- tmp1[t] * b2[,t]
+ }
+
+ fourth.tmp <- I_1(tmp2,env)
+
+ tmp3 <- b1b2(b1,b2,env)
+ tmp4 <- I_1(b3,env)
+ tmp5 <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ tmp5[,t] <- tmp3[t] * tmp4[,t]
+ }
+
+ fifth.tmp <- matrix(0,k.size,block)
+
+ for(k in 1:k.size){
+ fifth.tmp[k,] <- I0(tmp5[k,],env)
+ }
+
+ first <- first.tmp$first + second.tmp$first + third.tmp$first
+ second <- first.tmp$second + second.tmp$second + third.tmp$second +
+ fourth.tmp + fifth.tmp
+
+ return(list(first=first,second=second))
+ }
+
+
+ I_1_23_x <- function(x,get_I_1_23,env){
+
+ result <- I_123_x(x,get_I_1_23,env)
+
+ return(result)
+ }
+
+
+ #first:k*k*t, second:t
+
+ I_12_3 <- function(b1,b2,b3,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ block <- env$block
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) * (n5 - n6 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,block))
+ second <- matrix(0,k.size,block)
+
+ return(list(first=first,second=second))
+ }
+
+ first.tmp <- I_123(b1,b2,b3,env)
+ second.tmp <- I_123(b1,b3,b2,env)
+
+ tmp1 <- b1_b2(b2,b3,env)
+ tmp2 <- matrix(0,r.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ }
+
+ for(t in 1:block){
+ tmp2[,t] <- tmp1[t] * b1[,t]
+ }
+
+ third.tmp <- I_1(tmp2,env)
+
+ first <- first.tmp$first + second.tmp$first
+ second <- first.tmp$second + second.tmp$second + third.tmp
+
+ return(list(first=first,second=second))
+ }
+
+
+ I_12_3_x <- function(x,get_I_12_3,env){
+
+ result <- I_123_x(x,get_I_12_3,env)
+
+ return(result)
+ }
+
+
+ #first:k*k*t, second:t
+
+ I_1_2_3 <- function(b1,b2,b3,env){
+
+ k.size <- env$k.size
+ block <- env$block
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) * (n5 - n6 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,block))
+ second <- matrix(0,k.size,block)
+
+ return(list(first=first,second=second))
+ }
+
+ first.tmp <- I_123(b1,b2,b3,env)
+ second.tmp <- I_123(b2,b1,b3,env)
+ third.tmp <- I_123(b2,b3,b1,env)
+ fourth.tmp <- I_123(b1,b3,b2,env)
+ fifth.tmp <- I_123(b3,b1,b2,env)
+ sixth.tmp <- I_123(b3,b2,b1,env)
+
+ tmp1 <- b1_b2(b1,b3,env)
+ tmp2 <- I_1(b2,env)
+ seventh.tmp <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ seventh.tmp[,t] <- tmp1[t] * tmp2[,t]
+ }
+
+ tmp3 <- b1_b2(b1,b2,env)
+ tmp4 <- I_1(b3,env)
+ eighth.tmp <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ eighth.tmp[,t] <- tmp3[t] * tmp4[,t]
+ }
+
+ tmp5 <- b1_b2(b2,b3,env)
+ tmp6 <- I_1(b1,env)
+ ninth.tmp <- matrix(0,k.size,block)
+
+ for(t in 1:block){
+ ninth.tmp[,t] <- tmp5[t] * tmp6[,t]
+ }
+
+ first <- first.tmp$first + second.tmp$first + third.tmp$first +
+ fourth.tmp$first + fifth.tmp$first + sixth.tmp$first
+
+ second <- first.tmp$second + second.tmp$second + third.tmp$second +
+ fourth.tmp$second + fifth.tmp$second + sixth.tmp$second +
+ seventh.tmp + eighth.tmp + ninth.tmp
+
+ return(list(first=first,second=second))
+ }
+
+
+ I_1_2_3_x <- function(x,get_I_1_2_3,env){
+
+ result <- I_123_x(x,get_I_1_2_3,env)
+
+ return(result)
+ }
+
+
+ #first:k*k*k*k*t, second:k*k*t, third:t
+
+ I_12_34 <- function(b1,b2,b3,b4,env){
+
+ r.size <- env$r.size
+ k.size <- env$k.size
+ block <- env$block
+
+ n1 <- length(b1)
+ n2 <- sum(b1 == 0)
+
+ n3 <- length(b2)
+ n4 <- sum(b2 == 0)
+
+ n5 <- length(b3)
+ n6 <- sum(b3 == 0)
+
+ n7 <- length(b4)
+ n8 <- sum(b4 == 0)
+
+ n <- (n1 - n2 != 0) * (n3 - n4 != 0) *
+ (n5 - n6 != 0) * (n7 - n8 != 0)
+
+ if(n == 0){
+ first <- array(0,dim=c(k.size,k.size,k.size,k.size,block))
+ second <- array(0,dim=c(k.size,k.size,block))
+ third <- real(block)
+
+ return(list(first=first,second=second,third=third))
+ }
+
+ first.tmp <- I_1234(b1,b2,b3,b4,env)
+ second.tmp <- I_1234(b1,b3,b4,b2,env)
+ third.tmp <- I_1234(b1,b3,b2,b4,env)
+
+ tmp1 <- b1_b2(b2,b4,env)
+ tmp2 <- matrix(0,r.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ b2 <- matrix(b2,1,block)
+ b3 <- matrix(b3,1,block)
+ b4 <- matrix(b4,1,block)
+ }
+
+ for(t in 1:block){
+ tmp2[,t] <- tmp1[t] * b3[,t]
+ }
+
+ fourth.tmp <- I_12(b1,tmp2,env)
+
+ tmp3 <- b1_b2(b2,b3,env)
+ tmp4 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp4[,t] <- tmp3[t] * b1[,t]
+ }
+
+ fifth.tmp <- I_12(tmp4,b4,env)
+
+ tmp5 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp5[,t] <- tmp3[t] * b4[,t]
+ }
+
+ sixth.tmp <- I_12(b1,tmp5,env)
+
+ seventh.tmp <- I_1234(b3,b4,b1,b2,env)
+ eighth.tmp <- I_1234(b3,b1,b2,b4,env)
+ ninth.tmp <- I_1234(b3,b1,b4,b2,env)
+
+ tmp6 <- b1_b2(b4,b2,env)
+ tmp7 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp7[,t] <- tmp6[t] * b1[,t]
+ }
+
+ tenth.tmp <- I_12(b3,tmp7,env)
+
+ tmp8 <- b1_b2(b4,b1,env)
+ tmp9 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp9[,t] <- tmp8[t] * b3[,t]
+ }
+
+ eleventh.tmp <- I_12(tmp9,b2,env)
+
+ tmp10 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp10[,t] <- tmp8[t] * b2[,t]
+ }
+
+ twelfth.tmp <- I_12(b3,tmp10,env)
+
+ tmp11 <- b1_b2(b1,b3,env)
+ tmp12 <- I_12(b2,b4,env)
+ thirteenth.tmp <- tmp12
+
+ for(t in 1:block){
+ thirteenth.tmp$first[,,t] <- tmp11[t] * tmp12$first[,,t]
+ thirteenth.tmp$second[t] <- tmp11[t] * tmp12$second[t]
+ }
+
+ tmp13 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp13[,t] <- tmp11[t] * b2[,t]
+ }
+
+ fourteenth.tmp <- I_12(tmp13,b4,env)
+
+ tmp14 <- I_12(b4,b2,env)
+ fifteenth.tmp <- tmp14
+
+ for(t in 1:block){
+ fifteenth.tmp$first[,,t] <- tmp11[t] * tmp14$first[,,t]
+ fifteenth.tmp$second[t] <- tmp11[t] * tmp14$second[t]
+ }
+
+ tmp15 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp15[,t] <- tmp11[t] * b4[,t]
+ }
+
+ sixteenth.tmp <- I_12(tmp15,b2,env)
+
+ tmp16 <- b1_b2(b4,b2,env)
+ tmp17 <- matrix(0,r.size,block)
+
+ for(t in 1:block){
+ tmp17[,t] <- tmp16[t] * b3[,t]
+ }
+
+ tmp18 <- b1b2(b1,tmp17,env)
+
+ seventeenth.tmp <- I0(tmp18,env)
+
+ first <- first.tmp$first + second.tmp$first + third.tmp$first +
+ seventh.tmp$first + eighth.tmp$first + ninth.tmp$first
+
+ second <- first.tmp$second + second.tmp$second + third.tmp$second +
+ fourth.tmp$first + fifth.tmp$first - sixth.tmp$first +
+ seventh.tmp$second + eighth.tmp$second + ninth.tmp$second +
+ tenth.tmp$first + eleventh.tmp$first - twelfth.tmp$first +
+ thirteenth.tmp$first - fourteenth.tmp$first +
+ fifteenth.tmp$first - sixteenth.tmp$first
+
+ third <- first.tmp$third + second.tmp$third + third.tmp$third +
+ fourth.tmp$second + fifth.tmp$second - sixth.tmp$second +
+ seventh.tmp$third + eighth.tmp$third + ninth.tmp$third +
+ tenth.tmp$second + eleventh.tmp$second - twelfth.tmp$second +
+ thirteenth.tmp$second - fourteenth.tmp$second +
+ fifteenth.tmp$second - sixteenth.tmp$second +
+ seventeenth.tmp
+
+ return(list(first=first,second=second,third=third))
+ }
+
+
+ I_12_34_x <- function(x,get_I_12_34,env){
+
+ result <- I_1234_x(x,get_I_12_34,env)
+
+ return(result)
+ }
+
+
+ #first:k*k*k*k*t, second:k*k*t, third:t
+
+ I_1_2_34 <- function(b1,b2,b3,b4,env){
+
+ division <- env$division
+
+ first.tmp <- I_12_34(b1,b2,b3,b4,env)
+ second.tmp <- I_12_34(b2,b1,b3,b4,env)
+
+ tmp1 <- b1_b2(b1,b2,env)
+ tmp2 <- I_12(b3,b4,env)
+ third.tmp <- tmp2
+
+ for(t in 1:division){
+ third.tmp$first[,,t] <- tmp1[t] * tmp2$first[,,t]
+ third.tmp$second[t] <- tmp1[t] * tmp2$second[t]
+ }
+
+ first <- first.tmp$first + second.tmp$first
+ second <- first.tmp$second + second.tmp$second + third.tmp$first
+ third <- first.tmp$third + second.tmp$third + third.tmp$second
+
+ return(list(first=first,second=second,third=third))
+ }
+
+
+ I_1_2_34_x <- function(x,get_I_1_2_34,env){
+
+ result <- I_1234_x(x,get_I_1_2_34,env)
+
+ return(result)
+ }
+
+
+ I_1_2_3_4 <- function(b1,b2,b3,b4){
+
+
+ }
+
+
+##p16 Lemma3
+
+ I0_a <- function(b1,c1,env){
+
+ r.size <- env$r.size
+ block <- env$block
+
+ first.coef <- I0(c1,env)
+ first <- b1
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ }
+
+ second <- matrix(0,r.size,block)
+
+ for(r in 1:r.size){
+ second[r,] <- first.coef * b1[r,]
+ }
+
+ return(list(first.coef=first.coef,first=first,second=second))
+ }
+
+
+ I0_b <- function(b1,b2,c1,env){
+
+ r.size <- env$r.size
+ block <- env$block
+
+ first.coef <- I0(c1,env)
+ first <- list()
+
+ first[[1]] <- b1
+ first[[2]] <- b2
+
+ second <- list()
+
+ second[[1]] <- matrix(0,r.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ }
+
+ for(r in 1:r.size){
+ second[[1]][r,] <- first.coef * b1[r,]
+ }
+
+ second[[2]] <- b2
+
+ return(list(first.coef=first.coef,first=first,second=second))
+ }
+
+
+ I0_c <- function(b1,b2,c1,c2,env){
+
+ r.size <- env$r.size
+ block <- env$block
+
+ tmp1 <- I0(c2,env)
+
+ tmp2 <- c1 * tmp1
+
+ first.coef <- I0(tmp2,env)
+ first <- list()
+
+ first[[1]] <- b1
+ first[[2]] <- b2
+
+ second <- list()
+
+ second[[1]] <- matrix(0,r.size,block)
+
+ if(r.size == 1){
+ b1 <- matrix(b1,1,block)
+ }
+
+ for(r in 1:r.size){
+ second[[1]][r,] <- first.coef * b1[r,]
+ }
+
+ second[[2]] <- b2
+
+ third.coef <- I0(c1,env)
+ third <- list()
+
+ third[[1]] <- matrix(0,r.size,block)
+
+ for(r in 1:r.size){
+ third[[1]][r,] <- tmp1 * b1[r,]
+ }
+
+ third[[2]] <- b2
+
+ fourth <- list()
+
+ fourth[[1]] <- matrix(0,r.size,block)
+
+ for(r in 1:r.size){
+ fourth[[1]][r,] <- third.coef * tmp1 * b1[r,]
+ }
+
+ fourth[[2]] <- b2
+
+ return(list(first.coef=first.coef,first=first,second=second,
+ third.coef=third.coef,third=third,fourth=fourth))
+ }
+
+
+ #d*t
+
+ Y_e_V0 <- function(X.t0,de.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- matrix(0,d.size,block)
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+ for(d in 1:d.size){
+ tmp[d] <- eval(de.drift[[d]])
+ }
+ result[,t] <- invY[,,t] %*% tmp
+ }
+
+ return(result)
+ }
+
+
+ #d*r*t
+
+ Y_e_V <- function(X.t0,de.diffusion,env){
+
+ d.size <- env$d.size
+ r.size <- env$r.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,r.size,block))
+
+ tmp <- matrix(0,d.size,r.size)
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i in 1:d.size){
+ for(r in 1:r.size){
+ tmp[i,r] <- eval(de.diffusion[[i]][[r]])
+ }
+ }
+ result[,,t] <- invY[,,t] %*% tmp
+ }
+
+ return(result)
+ }
+
+
+ #d*t
+
+ Y_D <- function(X.t0,tmpY,get_Y_e_V0,env){
+
+ d.size <- env$d.size
+ delta <- env$delta
+ block <- env$block
+ Diff <- env$Diff
+
+ result <- matrix(0,d.size,block)
+
+ for(i in 1:d.size){
+ for(t in 2:block){
+ result[i,] <- get_Y_e_V0[i,] %*% Diff[t,] * delta
+ }
+ }
+
+ for(t in 2:block){
+ result[,t] <- tmpY[,,t] %*% result[,t]
+ }
+
+ return(result)
+ }
+
+
+ #d*d*d*t
+
+ Y_x1_x2_V0 <- function(X.t0,dxdx.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,d.size,block))
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ for(i2 in 1:d.size){
+ tmp.drift <- dxdx.drift[[i1]][[i2]]
+
+ for(d in 1:d.size){
+ tmp[d] <- eval(tmp.drift[[d]])
+ }
+
+ result[i1,i2,,t] <- invY[,,t] %*% tmp
+ }
+ }
+ }
+
+ return(result)
+ }
+
+
+ #d*d*d*r*t
+
+ Y_x1_x2_V <- function(X.t0,dxdx.diffusion,env){
+
+ d.size <- env$d.size
+ r.size <- env$r.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,d.size,r.size,block))
+
+ tmp <- matrix(0,d.size,r.size)
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ for(i2 in 1:d.size){
+ tmp.diffusion <- dxdx.diffusion[[i1]][[i2]]
+
+ for(i in 1:d.size){
+ for(r in 1:r.size){
+ tmp[i,r] <- eval(tmp.diffusion[[i]][[r]])
+ }
+ }
+ result[i1,i2,,,t] <- invY[,,t] %*% tmp
+ }
+ }
+ }
+
+ return(result)
+ }
+
+
+ #d*d*t
+
+ Y_x_e_V0 <- function(X.t0,dxde.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,block))
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ tmp.drift <- dxde.drift[[i1]]
+
+ for(d in 1:d.size){
+ tmp[d] <- eval(tmp.drift[[d]])
+ }
+
+ result[i1,,t] <- invY[,,t] %*% tmp
+ }
+ }
+
+ return(result)
+ }
+
+
+ #d*d*r*t
+
+ Y_x_e_V <- function(X.t0,dxde.diffusion,env){
+
+ d.size <- env$d.size
+ r.size <- env$r.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,r.size,block))
+
+ tmp <- matrix(0,d.size,r.size)
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ tmp.diffusion <- dxde.diffusion[[i1]]
+
+ for(i in 1:d.size){
+ for(r in 1:r.size){
+ tmp[i,r] <- eval(tmp.diffusion[[i]][[r]])
+ }
+ }
+ result[i1,,,t] <- invY[,,t] %*% tmp
+ }
+ }
+
+ return(result)
+ }
+
+
+ #d*t
+
+ Y_e_e_V0 <- function(X.t0,dede.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- matrix(0,d.size,block)
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(d in 1:d.size){
+ tmp[d] <- eval(dede.drift[[d]])
+ }
+
+ result[,t] <- invY[,,t] %*% tmp
+ }
+
+ return(result)
+ }
+
+
+ #d*r*t
+
+ Y_e_e_V <- function(X.t0,dede.diffusion,env){
+
+ d.size <- env$d.size
+ r.size <- env$r.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,r.size,block))
+
+ tmp <- matrix(0,d.size,r.size)
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i in 1:d.size){
+ for(r in 1:r.size){
+ tmp[i,r] <- eval(dede.diffusion[[i]][[r]])
+ }
+ }
+ result[,,t] <- invY[,,t] %*% tmp
+ }
+
+ return(result)
+ }
+
+
+ #d*d*d*d*t
+
+ Y_x1_x2_x3_V0 <- function(X.t0,dxdxdx.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,d.size,d.size,block))
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ for(i2 in 1:d.size){
+ for(i3 in 1:d.size){
+ tmp.drift <- dxdxdx.drift[[i1]][[i2]][[i3]]
+
+ for(d in 1:d.size){
+ tmp[d] <- eval(tmp.drift[[d]])
+ }
+
+ result[i1,i2,i3,,t] <- invY[,,t] %*% tmp
+ }
+ }
+ }
+ }
+
+ return(result)
+ }
+
+
+ #d*d*d*t
+
+ Y_x1_x2_e_V0 <- function(X.t0,dxdxde.drift,env){
+
+ d.size <- env$d.size
+ state <- env$state
+ pars <- env$pars
+ invY <- env$invY
+ block <- env$block
+ my.range <- env$my.range
+
+ result <- array(0,dim=c(d.size,d.size,d.size,block))
+
+ tmp <- c()
+
+ assign(pars[1],0)
+
+ for(t in 1:block){
+ for(d in 1:d.size){
+ assign(state[d],X.t0[my.range[t],d])
+ }
+
+ for(i1 in 1:d.size){
+ for(i2 in 1:d.size){
+ tmp.drift <- dxdxde.drift[[i1]][[i2]]
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 159
More information about the Yuima-commits
mailing list