R/asymptotic_term_second.R

In yuima: The YUIMA Project Package for SDEs

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@model)) stop("model object is missing!")
	if(is.null(yuima@sampling)) stop("sampling object is missing!")
	if(is.null(yuima@functional)) stop("functional object is missing!")

	f <- getf(yuima@functional)
	F <- getF(yuima@functional)

	##:: fix bug 07/23
	#e <- gete(yuima@functional)
	assign(expand.var, gete(yuima@functional))

	Terminal <- yuima@sampling@Terminal[1]
	division <- yuima@sampling@n[1]
	xinit <- getxinit(yuima@functional)
	state <- yuima@model@state.variable
	V0 <- yuima@model@drift
	V <- yuima@model@diffusion
	r.size <- yuima@model@noise.number
	d.size <- yuima@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@model@parameter@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 <- double(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 <- double(k.size)
		}else{
			n7 <- length(get_x_F)
			n8 <- sum(get_x_F == 0)

			n <- n7 - n8

			if(n == 0){
				mu2_1 <- double(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 <- double(k.size)
			}else{
				mu2_2 <- double(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))
		third.tmp <- matrix(0,k.size,d.size)

		n1 <- length(get_x_f0)
		n2 <- sum(get_x_f0 == 0)

		n <- n1 - n2

		third[,,division] <- 0

## modified 12/06

		if(n != 0){
		  for(s in (division-1):1){
		    third.tmp <- third.tmp + matrix(get_x_f0[,,s],k.size,d.size) %*%
				     tmpY[,,s] * delta

		    third[,,s] <- third.tmp %*% matrix(get_Y_e_V[,,s],d.size,r.size)
		  }
		}

		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
	## 

	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  'cubature' 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, lowerLimit=lower,upperLimit=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 <- double(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 <- double(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 <- double(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
	## 

	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 <- Terminal/(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, Sigma=Sigma)
	return(terms)
})