[Yuima-commits] r581 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 25 11:55:48 CET 2017
Author: kyuta
Date: 2017-01-25 11:55:48 +0100 (Wed, 25 Jan 2017)
New Revision: 581
Added:
pkg/yuima/src/euler.c
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/sim.euler.R
Log:
modified sim.euler.R and added euler.c to implement the Euler-Maruyama scheme by the C code (only the diffusion case)
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-01-24 18:04:37 UTC (rev 580)
+++ pkg/yuima/DESCRIPTION 2017-01-25 10:55:48 UTC (rev 581)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.5.5
+Version: 1.5.6
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2017-01-24 18:04:37 UTC (rev 580)
+++ pkg/yuima/NEWS 2017-01-25 10:55:48 UTC (rev 581)
@@ -41,5 +41,6 @@
2016/01/14: modified rng.R, rng.Rd, adaBayes.Rd, qmle.Rd, setModel.Rd, spectralcov.Rd
2016/05/26: added rpts and rnts in rng.R and the corresponding c language file
2016/07/08: fixed some bugs in llag.R and cce_functions.c
-2016/10/04: modified setMultiModel.R, sim.euler.R and yuima.model.R to generate nts and pts process
-2016/12/16: added rGIG, rGH, dGIG and dGH in rng.R and the corresponding c language file YU
\ No newline at end of file
+2016/10/04: modified setMultiModel.R, sim.euler.R and yuima.model.R to generate nts and pts process
+2016/12/16: added rGIG, rGH, dGIG and dGH in rng.R and the corresponding c language file YU
+2017/01/25: modified sim.euler.R and added euler.c to implement the Euler-Maruyama scheme by the C code (only the diffusion case)
\ No newline at end of file
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2017-01-24 18:04:37 UTC (rev 580)
+++ pkg/yuima/R/sim.euler.R 2017-01-25 10:55:48 UTC (rev 581)
@@ -1,374 +1,385 @@
-euler<-function(xinit,yuima,dW,env){
-
- sdeModel<-yuima at model
-
- modelstate <- sdeModel at solve.variable
- modeltime <- sdeModel at time.variable
- V0 <- sdeModel at drift
- V <- sdeModel at diffusion
- r.size <- sdeModel at noise.number
- d.size <- sdeModel at equation.number
- Terminal <- yuima at sampling@Terminal[1]
- Initial <- yuima at sampling@Initial[1]
-
- n <- yuima at sampling@n[1]
- dL <- env$dL
-
-# dX <- xinit
-
- # 06/11 xinit is an expression: the structure is equal to that of V0
- if(length(unique(as.character(xinit)))==1 &&
- is.numeric(tryCatch(eval(xinit[1],env),error=function(...) FALSE))){
- dX_dummy<-xinit[1]
- dummy.val<-eval(dX_dummy, env)
- if(length(dummy.val)==1){dummy.val<-rep(dummy.val,length(xinit))}
- for(i in 1:length(modelstate)){
- assign(modelstate[i],dummy.val[i] ,env)
- }
- dX<-vector(mode="numeric",length(dX_dummy))
-
- for(i in 1:length(xinit)){
- dX[i] <- dummy.val[i]
- }
- }else{
- dX_dummy <- xinit
- if(length(modelstate)==length(dX_dummy)){
- for(i in 1:length(modelstate)) {
- if(is.numeric(tryCatch(eval(dX_dummy[i],env),error=function(...) FALSE))){
- assign(modelstate[i], eval(dX_dummy[i], env),env)
- }else{
- assign(modelstate[i], 0, env)
- }
- }
- }else{
- yuima.warn("the number of model states do not match the number of initial conditions")
- return(NULL)
- }
-
- # 06/11 we need a initial variable for X_0
- dX<-vector(mode="numeric",length(dX_dummy))
-
- for(i in 1:length(dX_dummy)){
- dX[i] <- eval(dX_dummy[i], env)
- }
- }
-##:: set time step
- # delta <- Terminal/n
- delta <- yuima at sampling@delta
-
-##:: check if DRIFT and/or DIFFUSION has values
- has.drift <- sum(as.character(sdeModel at drift) != "(0)")
- var.in.diff <- is.logical(any(match(unlist(lapply(sdeModel at diffusion, all.vars)), sdeModel at state.variable)))
- #print(is.Poisson(sdeModel))
-
- ##:: function to calculate coefficients of dW(including drift term)
- ##:: common used in Wiener and CP
- p.b <- function(t, X=numeric(d.size)){
- ##:: assign names of variables
- for(i in 1:length(modelstate)){
- assign(modelstate[i], X[i], env)
- }
- assign(modeltime, t, env)
- ##:: solve diffusion term
- if(has.drift){
- tmp <- matrix(0, d.size, r.size+1)
- for(i in 1:d.size){
- tmp[i,1] <- eval(V0[i], env)
- for(j in 1:r.size){
- tmp[i,j+1] <- eval(V[[i]][j],env)
- }
- }
- } else { ##:: no drift term (faster)
- tmp <- matrix(0, d.size, r.size)
- if(!is.Poisson(sdeModel)){ # we do not need to evaluate diffusion
- for(i in 1:d.size){
- for(j in 1:r.size){
- tmp[i,j] <- eval(V[[i]][j],env)
- } # for j
- } # foh i
- } # !is.Poisson
- } # else
- return(tmp)
- }
-
- X_mat <- matrix(0, d.size, n+1)
- X_mat[,1] <- dX
-
- if(has.drift){ ##:: consider drift term to be one of the diffusion term(dW=1)
- dW <- rbind( rep(1, n)*delta , dW)
- }
-
-
- if(!length(sdeModel at measure.type)){ ##:: Wiener Proc
- ##:: using Euler-Maruyama method
-
- if(var.in.diff & (!is.Poisson(sdeModel))){ ##:: diffusions have state variables and it is not Poisson
- ##:: calcurate difference eq.
- for( i in 1:n){
- # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i]
- dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] # LM
- X_mat[,i+1] <- dX
- }
- }else{ ##:: diffusions have no state variables (not use p.b(). faster)
- sde.tics <- seq(0, Terminal, length=(n+1))
- sde.tics <- sde.tics[2:length(sde.tics)]
-
- X_mat[, 1] <- dX
-
- ##:: assign names of variables
- for(i in 1:length(modelstate)){
- assign(modelstate[i], dX[i])
- }
- assign(modeltime, sde.tics)
- t.size <- length(sde.tics)
-
- ##:: solve diffusion term
- if(has.drift){
- pbdata <- matrix(0, d.size*(r.size+1), t.size)
- for(i in 1:d.size){
- pbdata[(i-1)*(r.size+1)+1, ] <- eval(V0[i], env)
- for(j in 1:r.size){
- pbdata[(i-1)*(r.size+1)+j+1, ] <- eval(V[[i]][j], env)
- }
- }
- dim(pbdata)<-(c(r.size+1, d.size*t.size))
- }else{
- pbdata <- matrix(0, d.size*r.size, t.size)
- if(!is.Poisson(sdeModel)){
- for(i in 1:d.size){
- for(j in 1:r.size){
- pbdata[(i-1)*r.size+j, ] <- eval(V[[i]][j], env)
- } # for j
- } # for i
- } # !is.Poisson
- dim(pbdata)<-(c(r.size, d.size*t.size))
- } # else
-
- pbdata <- t(pbdata)
-
- ##:: calcurate difference eq.
- for( i in 1:n){
- if(!is.Poisson(sdeModel))
- dX <- dX + pbdata[((i-1)*d.size+1):(i*d.size), ] %*% dW[, i]
- X_mat[, i+1] <- dX
- }
- }
- #tsX <- ts(data=t(X_mat), deltat=delta , start=0)
- tsX <- ts(data=t(X_mat), deltat=delta , start = yuima at sampling@Initial) #LM
- }else{ ##:: Levy
- JP <- sdeModel at jump.coeff
- mu.size <- length(JP)
- # cat("\n Levy\n")
- ##:: function to solve c(x,z)
- p.b.j <- function(t, X=numeric(d.size)){
- for(i in 1:length(modelstate)){
- assign(modelstate[i], X[i], env)
- }
- assign(modeltime, t, env)
- # tmp <- numeric(d.size)
- j.size <- length(JP[[1]])
- tmp <- matrix(0, mu.size, j.size)
- # cat("\n inside\n")
- #print(dim(tmp))
- for(i in 1:mu.size){
- for(j in 1:j.size){
- tmp[i,j] <- eval(JP[[i]][j],env)
- }
- # tmp[i] <- eval(JP[i], env)
- }
- return(tmp)
- }
- # print(ls(env))
-
- ### WHY I AM DOING THIS?
- # tmp <- matrix(0, d.size, r.size)
- #
- #for(i in 1:d.size){
- # for(j in 1:r.size){
- # cat("\n here\n")
- # tmp[i,j] <- eval(V[[i]][j],env)
- # } # for j
- # }
- ###
-
- if(sdeModel at measure.type == "CP" ){ ##:: Compound-Poisson type
-
- ##:: delete 2010/09/13 for simulate func bug fix by s.h
- ## eta0 <- eval(sdeModel at measure$intensity)
-
- ##:: add 2010/09/13 for simulate func bug fix by s.h
- eta0 <- eval(sdeModel at measure$intensity, env) ## intensity param
-
- ##:: get lambda from nu()
- tmp.expr <- function(my.x){
- assign(sdeModel at jump.variable,my.x)
- return(eval(sdeModel at measure$df$expr))
- }
- #lambda <- integrate(sdeModel at measure$df$func, 0, Inf)$value * eta0
- #lambda <- integrate(tmp.expr, 0, Inf)$value * eta0 ##bug:2013/10/28
-
- dummyList<-as.list(env)
- # print(str(dummyList))
- #print(str(idx.dummy))
- lgth.meas<-length(yuima at model@parameter at measure)
- if(lgth.meas>1){
- for(i in c(2:lgth.meas)){
- idx.dummy<-yuima at model@parameter at measure[i]
- #print(i)
- #print(yuima at model@parameter at measure[i])
- assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
- }
- }
-
-
- lambda <- integrate(tmp.expr, -Inf, Inf)$value * eta0
-
- ##:: lambda = nu() (p6)
- # N_sharp <- rpois(1,Terminal*eta0) ##:: Po(Ne)
- N_sharp <- rpois(1,(Terminal-Initial)*eta0) ##:: Po(Ne)
- if(N_sharp == 0){
- JAMP <- FALSE
- }else{
- JAMP <- TRUE
- Uj <- sort( runif(N_sharp, Initial, Terminal) )
- # Uj <- sort( runif(N_sharp, 0, Terminal) )
- ij <- NULL
- for(i in 1:length(Uj)){
- Min <- min(which(Initial+ c(1:n)*delta > Uj[i]))
- # Min <- min(which(c(1:n)*delta > Uj[i]))
- ij <- c(ij, Min)
- }
- }
-
- ##:: make expression to create iid rand J
- if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel at measure$df$expr)){
- ##:: e.g. dnorm(z,1,1) -> rnorm(mu.size*N_sharp,1,1)
- F <- suppressWarnings(parse(text=gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel at measure$df$expr, perl=TRUE)))
- }else{
- stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
- }
-
- ##:: delete 2010/09/13 for simulate func bug fix by s.h
- ## randJ <- eval(F) ## this expression is evaluated locally not in the yuimaEnv
-
- ##:: add 2010/09/13 for simulate func bug fix by s.h
- F.env <- new.env(parent=env)
- assign("mu.size", mu.size, envir=F.env)
- assign("N_sharp", N_sharp, envir=F.env)
-
- randJ <- eval(F, F.env) ## this expression is evaluated in the F.env
-
- j <- 1
- for(i in 1:n){
- if(JAMP==FALSE || sum(i==ij)==0){
- Pi <- 0
- }else{
- if(is.null(dL)){
- J <- eta0*randJ[j]/lambda
- j <- j+1
- ##cat(paste(J,"\n"))
- ##Pi <- zeta(dX, J)
- assign(sdeModel at jump.variable, J, env)
-
- if(sdeModel at J.flag){
- J <- 1
- }
-
- # Pi <- p.b.j(t=i*delta,X=dX) * J #LM
- Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) * J
- }else{# we add this part since we allow the user to specify the increment of CP LM 05/02/2015
- # Pi <- p.b.j(t=i*delta,X=dX) %*% dL[, i] #LM
- Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) %*% dL[, i]
- }
- ##Pi <- p.b.j(t=i*delta, X=dX)
- }
- # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi # LM
- dX <- dX + p.b(t=yuima at sampling@Initial + i*delta, X=dX) %*% dW[, i] + Pi
- X_mat[, i+1] <- dX
- }
- # tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
- tsX <- ts(data=t(X_mat), deltat=delta, start=yuima at sampling@Initial)
- ##::end CP
- }else if(sdeModel at measure.type=="code"){ ##:: code type
- ##:: Jump terms
- code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
- args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
- #print(args)
- dZ <- switch(code,
- rNIG=paste("rNIG(n, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta, ", args[6],")"),
- rIG=paste("rIG(n,", args[2], "*delta, ", args[3], ")"),
- rgamma=paste("rgamma(n, ", args[2], "*delta, ", args[3], ")"),
- rbgamma=paste("rbgamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], "*delta, ", args[5], ")"),
-## rngamma=paste("rngamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
- rvgamma=paste("rvgamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta,", args[6],")"),
-## rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
- rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], "*delta^(1/",args[2],"), ", args[5], "*delta)"),
- rpts=paste("rpts(n, ", args[2], ", ", args[3], "*delta,", args[4],")"),
- rnts=paste("rnts(n, ", args[2], ", ", args[3], "*delta,", args[4], ", ", args[5], ", ", args[6],"*delta,", args[7], ")")
- )
- ## added "rpts" and "rnts" by YU (2016/10/4)
- dummyList<-as.list(env)
- #print(str(dummyList))
- lgth.meas<-length(yuima at model@parameter at measure)
- #print(lgth.meas)
- if(lgth.meas!=0){
- for(i in c(1:lgth.meas)){
- #print(i)
- #print(yuima at model@parameter at measure[i])
- idx.dummy<-yuima at model@parameter at measure[i]
- #print(str(idx.dummy))
- assign(idx.dummy,dummyList[[idx.dummy]])
- #print(str(idx.dummy))
- #print(str(dummyList[[idx.dummy]]))
- #print(get(idx.dummy))
- }
- }
-
- if(is.null(dZ)){ ##:: "otherwise"
- cat(paste("Code \"", code, "\" not supported yet.\n", sep=""))
- return(NULL)
- }
- if(!is.null(dL))
- dZ <- dL
- else
- dZ <- eval(parse(text=dZ))
- ##:: calcurate difference eq.
- #print(str(dZ))
- if(is.null(dim(dZ)))
- dZ <- matrix(dZ,nrow=1)
- # print(dim(dZ))
- # print(str(sdeModel at jump.variable))
- for(i in 1:n){
- assign(sdeModel at jump.variable, dZ[,i], env)
-
- if(sdeModel at J.flag){
- dZ[,i] <- 1
- }
- # cat("\np.b.j call\n")
- # tmp.j <- p.b.j(t=i*delta, X=dX) #LM
- tmp.j <- p.b.j(t=yuima at sampling@Initial+i*delta, X=dX)
- #print(str(tmp.j))
- #cat("\np.b.j cback and dZ\n")
- # print(str(dZ[,i]))
- # print(sum(dim(tmp.j)))
- if(sum(dim(tmp.j))==2)
- tmp.j <- as.numeric(tmp.j)
- #print(str(tmp.j))
- #print(str(p.b(t = i * delta, X = dX) %*% dW[, i]))
- # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i] #LM
- dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i]
- X_mat[, i+1] <- dX
- }
- # tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
- tsX <- ts(data=t(X_mat), deltat=delta, start=yuima at sampling@Initial)
- ##::end code
- }else{
- cat(paste("Type \"", sdeModel at measure.type, "\" not supported yet.\n", sep=""))
- return(NULL)
- }
- }##::end levy
- yuimaData <- setData(original.data=tsX)
- return(yuimaData)
-
-
-}
+euler<-function(xinit,yuima,dW,env){
+
+ sdeModel<-yuima at model
+
+ modelstate <- sdeModel at solve.variable
+ modeltime <- sdeModel at time.variable
+ V0 <- sdeModel at drift
+ V <- sdeModel at diffusion
+ r.size <- sdeModel at noise.number
+ d.size <- sdeModel at equation.number
+ Terminal <- yuima at sampling@Terminal[1]
+ Initial <- yuima at sampling@Initial[1]
+
+ n <- yuima at sampling@n[1]
+ dL <- env$dL
+
+# dX <- xinit
+
+ # 06/11 xinit is an expression: the structure is equal to that of V0
+ if(length(unique(as.character(xinit)))==1 &&
+ is.numeric(tryCatch(eval(xinit[1],env),error=function(...) FALSE))){
+ dX_dummy<-xinit[1]
+ dummy.val<-eval(dX_dummy, env)
+ if(length(dummy.val)==1){dummy.val<-rep(dummy.val,length(xinit))}
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i],dummy.val[i] ,env)
+ }
+ dX<-vector(mode="numeric",length(dX_dummy))
+
+ for(i in 1:length(xinit)){
+ dX[i] <- dummy.val[i]
+ }
+ }else{
+ dX_dummy <- xinit
+ if(length(modelstate)==length(dX_dummy)){
+ for(i in 1:length(modelstate)) {
+ if(is.numeric(tryCatch(eval(dX_dummy[i],env),error=function(...) FALSE))){
+ assign(modelstate[i], eval(dX_dummy[i], env),env)
+ }else{
+ assign(modelstate[i], 0, env)
+ }
+ }
+ }else{
+ yuima.warn("the number of model states do not match the number of initial conditions")
+ return(NULL)
+ }
+
+ # 06/11 we need a initial variable for X_0
+ dX<-vector(mode="numeric",length(dX_dummy))
+
+ for(i in 1:length(dX_dummy)){
+ dX[i] <- eval(dX_dummy[i], env)
+ }
+ }
+##:: set time step
+ # delta <- Terminal/n
+ delta <- yuima at sampling@delta
+
+##:: check if DRIFT and/or DIFFUSION has values
+ has.drift <- sum(as.character(sdeModel at drift) != "(0)")
+ var.in.diff <- is.logical(any(match(unlist(lapply(sdeModel at diffusion, all.vars)), sdeModel at state.variable)))
+ #print(is.Poisson(sdeModel))
+
+ ##:: function to calculate coefficients of dW(including drift term)
+ ##:: common used in Wiener and CP
+ p.b <- function(t, X=numeric(d.size)){
+ ##:: assign names of variables
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], X[i], env)
+ }
+ assign(modeltime, t, env)
+ ##:: solve diffusion term
+ if(has.drift){
+ tmp <- matrix(0, d.size, r.size+1)
+ for(i in 1:d.size){
+ tmp[i,1] <- eval(V0[i], env)
+ for(j in 1:r.size){
+ tmp[i,j+1] <- eval(V[[i]][j],env)
+ }
+ }
+ } else { ##:: no drift term (faster)
+ tmp <- matrix(0, d.size, r.size)
+ if(!is.Poisson(sdeModel)){ # we do not need to evaluate diffusion
+ for(i in 1:d.size){
+ for(j in 1:r.size){
+ tmp[i,j] <- eval(V[[i]][j],env)
+ } # for j
+ } # foh i
+ } # !is.Poisson
+ } # else
+ return(tmp)
+ }
+
+ X_mat <- matrix(0, d.size, n+1)
+ X_mat[,1] <- dX
+
+ if(has.drift){ ##:: consider drift term to be one of the diffusion term(dW=1)
+ dW <- rbind( rep(1, n)*delta , dW)
+ }
+
+
+ if(!length(sdeModel at measure.type)){ ##:: Wiener Proc
+ ##:: using Euler-Maruyama method
+
+ if(0){ # old version (Jan 25, 2017)
+ if(var.in.diff & (!is.Poisson(sdeModel))){ ##:: diffusions have state variables and it is not Poisson
+ ##:: calcurate difference eq.
+ for( i in 1:n){
+ # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i]
+ dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] # LM
+ X_mat[,i+1] <- dX
+ }
+ }else{ ##:: diffusions have no state variables (not use p.b(). faster)
+ sde.tics <- seq(0, Terminal, length=(n+1))
+ sde.tics <- sde.tics[2:length(sde.tics)]
+
+ X_mat[, 1] <- dX
+
+ ##:: assign names of variables
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], dX[i])
+ }
+ assign(modeltime, sde.tics)
+ t.size <- length(sde.tics)
+
+ ##:: solve diffusion term
+ if(has.drift){
+ pbdata <- matrix(0, d.size*(r.size+1), t.size)
+ for(i in 1:d.size){
+ pbdata[(i-1)*(r.size+1)+1, ] <- eval(V0[i], env)
+ for(j in 1:r.size){
+ pbdata[(i-1)*(r.size+1)+j+1, ] <- eval(V[[i]][j], env)
+ }
+ }
+ dim(pbdata)<-(c(r.size+1, d.size*t.size))
+ }else{
+ pbdata <- matrix(0, d.size*r.size, t.size)
+ if(!is.Poisson(sdeModel)){
+ for(i in 1:d.size){
+ for(j in 1:r.size){
+ pbdata[(i-1)*r.size+j, ] <- eval(V[[i]][j], env)
+ } # for j
+ } # for i
+ } # !is.Poisson
+ dim(pbdata)<-(c(r.size, d.size*t.size))
+ } # else
+
+ pbdata <- t(pbdata)
+
+ ##:: calcurate difference eq.
+ for( i in 1:n){
+ if(!is.Poisson(sdeModel))
+ dX <- dX + pbdata[((i-1)*d.size+1):(i*d.size), ] %*% dW[, i]
+ X_mat[, i+1] <- dX
+ }
+ }
+ }
+
+ # new version (Jan 25, 2017)
+ b <- parse(text=paste("c(",paste(as.character(V0),collapse=","),")"))
+ vecV <- parse(text=paste("c(",paste(as.character(unlist(V)),collapse=","),")"))
+
+ X_mat <- .Call("euler", dX, Initial, as.integer(r.size),
+ rep(1, n) * delta, dW, modeltime, modelstate, quote(eval(b, env)),
+ quote(eval(vecV, env)), env, new.env())
+
+ #tsX <- ts(data=t(X_mat), deltat=delta , start=0)
+ tsX <- ts(data=t(X_mat), deltat=delta , start = yuima at sampling@Initial) #LM
+ }else{ ##:: Levy
+ JP <- sdeModel at jump.coeff
+ mu.size <- length(JP)
+ # cat("\n Levy\n")
+ ##:: function to solve c(x,z)
+ p.b.j <- function(t, X=numeric(d.size)){
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], X[i], env)
+ }
+ assign(modeltime, t, env)
+ # tmp <- numeric(d.size)
+ j.size <- length(JP[[1]])
+ tmp <- matrix(0, mu.size, j.size)
+ # cat("\n inside\n")
+ #print(dim(tmp))
+ for(i in 1:mu.size){
+ for(j in 1:j.size){
+ tmp[i,j] <- eval(JP[[i]][j],env)
+ }
+ # tmp[i] <- eval(JP[i], env)
+ }
+ return(tmp)
+ }
+ # print(ls(env))
+
+ ### WHY I AM DOING THIS?
+ # tmp <- matrix(0, d.size, r.size)
+ #
+ #for(i in 1:d.size){
+ # for(j in 1:r.size){
+ # cat("\n here\n")
+ # tmp[i,j] <- eval(V[[i]][j],env)
+ # } # for j
+ # }
+ ###
+
+ if(sdeModel at measure.type == "CP" ){ ##:: Compound-Poisson type
+
+ ##:: delete 2010/09/13 for simulate func bug fix by s.h
+ ## eta0 <- eval(sdeModel at measure$intensity)
+
+ ##:: add 2010/09/13 for simulate func bug fix by s.h
+ eta0 <- eval(sdeModel at measure$intensity, env) ## intensity param
+
+ ##:: get lambda from nu()
+ tmp.expr <- function(my.x){
+ assign(sdeModel at jump.variable,my.x)
+ return(eval(sdeModel at measure$df$expr))
+ }
+ #lambda <- integrate(sdeModel at measure$df$func, 0, Inf)$value * eta0
+ #lambda <- integrate(tmp.expr, 0, Inf)$value * eta0 ##bug:2013/10/28
+
+ dummyList<-as.list(env)
+ # print(str(dummyList))
+ #print(str(idx.dummy))
+ lgth.meas<-length(yuima at model@parameter at measure)
+ if(lgth.meas>1){
+ for(i in c(2:lgth.meas)){
+ idx.dummy<-yuima at model@parameter at measure[i]
+ #print(i)
+ #print(yuima at model@parameter at measure[i])
+ assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
+ }
+ }
+
+
+ lambda <- integrate(tmp.expr, -Inf, Inf)$value * eta0
+
+ ##:: lambda = nu() (p6)
+ # N_sharp <- rpois(1,Terminal*eta0) ##:: Po(Ne)
+ N_sharp <- rpois(1,(Terminal-Initial)*eta0) ##:: Po(Ne)
+ if(N_sharp == 0){
+ JAMP <- FALSE
+ }else{
+ JAMP <- TRUE
+ Uj <- sort( runif(N_sharp, Initial, Terminal) )
+ # Uj <- sort( runif(N_sharp, 0, Terminal) )
+ ij <- NULL
+ for(i in 1:length(Uj)){
+ Min <- min(which(Initial+ c(1:n)*delta > Uj[i]))
+ # Min <- min(which(c(1:n)*delta > Uj[i]))
+ ij <- c(ij, Min)
+ }
+ }
+
+ ##:: make expression to create iid rand J
+ if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel at measure$df$expr)){
+ ##:: e.g. dnorm(z,1,1) -> rnorm(mu.size*N_sharp,1,1)
+ F <- suppressWarnings(parse(text=gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel at measure$df$expr, perl=TRUE)))
+ }else{
+ stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
+ }
+
+ ##:: delete 2010/09/13 for simulate func bug fix by s.h
+ ## randJ <- eval(F) ## this expression is evaluated locally not in the yuimaEnv
+
+ ##:: add 2010/09/13 for simulate func bug fix by s.h
+ F.env <- new.env(parent=env)
+ assign("mu.size", mu.size, envir=F.env)
+ assign("N_sharp", N_sharp, envir=F.env)
+
+ randJ <- eval(F, F.env) ## this expression is evaluated in the F.env
+
+ j <- 1
+ for(i in 1:n){
+ if(JAMP==FALSE || sum(i==ij)==0){
+ Pi <- 0
+ }else{
+ if(is.null(dL)){
+ J <- eta0*randJ[j]/lambda
+ j <- j+1
+ ##cat(paste(J,"\n"))
+ ##Pi <- zeta(dX, J)
+ assign(sdeModel at jump.variable, J, env)
+
+ if(sdeModel at J.flag){
+ J <- 1
+ }
+
+ # Pi <- p.b.j(t=i*delta,X=dX) * J #LM
+ Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) * J
+ }else{# we add this part since we allow the user to specify the increment of CP LM 05/02/2015
+ # Pi <- p.b.j(t=i*delta,X=dX) %*% dL[, i] #LM
+ Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) %*% dL[, i]
+ }
+ ##Pi <- p.b.j(t=i*delta, X=dX)
+ }
+ # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi # LM
+ dX <- dX + p.b(t=yuima at sampling@Initial + i*delta, X=dX) %*% dW[, i] + Pi
+ X_mat[, i+1] <- dX
+ }
+ # tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
+ tsX <- ts(data=t(X_mat), deltat=delta, start=yuima at sampling@Initial)
+ ##::end CP
+ }else if(sdeModel at measure.type=="code"){ ##:: code type
+ ##:: Jump terms
+ code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
+ args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
+ #print(args)
+ dZ <- switch(code,
+ rNIG=paste("rNIG(n, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta, ", args[6],")"),
+ rIG=paste("rIG(n,", args[2], "*delta, ", args[3], ")"),
+ rgamma=paste("rgamma(n, ", args[2], "*delta, ", args[3], ")"),
+ rbgamma=paste("rbgamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], "*delta, ", args[5], ")"),
+## rngamma=paste("rngamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
+ rvgamma=paste("rvgamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta,", args[6],")"),
+## rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
+ rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], "*delta^(1/",args[2],"), ", args[5], "*delta)"),
+ rpts=paste("rpts(n, ", args[2], ", ", args[3], "*delta,", args[4],")"),
+ rnts=paste("rnts(n, ", args[2], ", ", args[3], "*delta,", args[4], ", ", args[5], ", ", args[6],"*delta,", args[7], ")")
+ )
+ ## added "rpts" and "rnts" by YU (2016/10/4)
+ dummyList<-as.list(env)
+ #print(str(dummyList))
+ lgth.meas<-length(yuima at model@parameter at measure)
+ #print(lgth.meas)
+ if(lgth.meas!=0){
+ for(i in c(1:lgth.meas)){
+ #print(i)
+ #print(yuima at model@parameter at measure[i])
+ idx.dummy<-yuima at model@parameter at measure[i]
+ #print(str(idx.dummy))
+ assign(idx.dummy,dummyList[[idx.dummy]])
+ #print(str(idx.dummy))
+ #print(str(dummyList[[idx.dummy]]))
+ #print(get(idx.dummy))
+ }
+ }
+
+ if(is.null(dZ)){ ##:: "otherwise"
+ cat(paste("Code \"", code, "\" not supported yet.\n", sep=""))
+ return(NULL)
+ }
+ if(!is.null(dL))
+ dZ <- dL
+ else
+ dZ <- eval(parse(text=dZ))
+ ##:: calcurate difference eq.
+ #print(str(dZ))
+ if(is.null(dim(dZ)))
+ dZ <- matrix(dZ,nrow=1)
+ # print(dim(dZ))
+ # print(str(sdeModel at jump.variable))
+ for(i in 1:n){
+ assign(sdeModel at jump.variable, dZ[,i], env)
+
+ if(sdeModel at J.flag){
+ dZ[,i] <- 1
+ }
+ # cat("\np.b.j call\n")
+ # tmp.j <- p.b.j(t=i*delta, X=dX) #LM
+ tmp.j <- p.b.j(t=yuima at sampling@Initial+i*delta, X=dX)
+ #print(str(tmp.j))
+ #cat("\np.b.j cback and dZ\n")
+ # print(str(dZ[,i]))
+ # print(sum(dim(tmp.j)))
+ if(sum(dim(tmp.j))==2)
+ tmp.j <- as.numeric(tmp.j)
+ #print(str(tmp.j))
+ #print(str(p.b(t = i * delta, X = dX) %*% dW[, i]))
+ # dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i] #LM
+ dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i]
+ X_mat[, i+1] <- dX
+ }
+ # tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
+ tsX <- ts(data=t(X_mat), deltat=delta, start=yuima at sampling@Initial)
+ ##::end code
+ }else{
+ cat(paste("Type \"", sdeModel at measure.type, "\" not supported yet.\n", sep=""))
+ return(NULL)
+ }
+ }##::end levy
+ yuimaData <- setData(original.data=tsX)
+ return(yuimaData)
+
+
+}
Added: pkg/yuima/src/euler.c
===================================================================
--- pkg/yuima/src/euler.c (rev 0)
+++ pkg/yuima/src/euler.c 2017-01-25 10:55:48 UTC (rev 581)
@@ -0,0 +1,83 @@
+//
+// euler.c
+//
+//
+// Created by Yuta Koike on 2016/01/23.
+//
+//
+
+#include <R.h>
+#include <Rmath.h>
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+SEXP euler(SEXP x0, SEXP t0, SEXP R, SEXP dt, SEXP dW, SEXP modeltime, SEXP modelstate,
+ SEXP drift, SEXP diffusion, SEXP env, SEXP rho){
+
+ int i, j, k, n, d, r;
+ double *rdt, *rdW, *rX, *rx0, *b, *sigma;
+ SEXP X, xpar, tpar;
+
+ PROTECT(x0 = AS_NUMERIC(x0));
+ rx0 = REAL(x0);
+ PROTECT(dt = AS_NUMERIC(dt));
+ rdt = REAL(dt);
+ PROTECT(dW = AS_NUMERIC(dW));
+ rdW = REAL(dW);
+
+ d = length(x0);
+ n = length(dt);
+
+ r = *INTEGER(R);
+
+ /* PROTECT(X = allocVector(REALSXP, n + 1));
+ rX = REAL(X);
+ rX[0] = *REAL(x0);*/
+ PROTECT(X = allocMatrix(REALSXP, d, n + 1));
+ rX = REAL(X);
+
+ for (j = 0; j < d; j++) {
+ rX[j] = rx0[j];
+ }
+
+ PROTECT(tpar = allocVector(REALSXP, 1));
+ PROTECT(xpar = allocVector(REALSXP, 1));
+
+ PROTECT(t0 = AS_NUMERIC(t0));
+ REAL(tpar)[0] = REAL(t0)[0]; /* initial time */
+
+ for (i = 0; i < n; i++) {
+
+ /* assign the current variables */
+ defineVar(installChar(STRING_ELT(modeltime, 0)), tpar, env);
+
+ for (j = 0; j < d; j++) {
+ /*xpar = allocVector(REALSXP, 1);*/
+ REAL(xpar)[0] = rX[j + i * d];
+ defineVar(installChar(STRING_ELT(modelstate, j)), duplicate(xpar), env);
+ /*defineVar(installChar(STRING_ELT(modelstate, j)), xpar, env);*/
+ }
+
+ /*defineVar(install("env"), env, rho);*/
+
+ /* evaluate coefficients */
+ b = REAL(eval(drift, rho));
+ sigma = REAL(eval(diffusion, rho));
+
+ for (j = 0; j < d; j++) {
+ rX[j + (i + 1) * d] = rX[j + i * d] + b[j] * rdt[i];
+ for (k = 0; k < r; k++) {
+ rX[j + (i + 1) * d] += sigma[k + j * r] * rdW[k + i * r];
+ }
+ }
+
+ /*defineVar(install("t"), tpar, rho);*/
+ /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * rdt[i] + *REAL(eval(diffusion, rho)) * rdW[i];*/
+ /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * REAL(dt)[i] + *REAL(eval(diffusion, rho)) * REAL(dW)[i];*/
+
+ REAL(tpar)[0] += rdt[i];
+ }
+
+ UNPROTECT(7);
+ return(X);
+}
More information about the Yuima-commits
mailing list