Nothing
R/simulateMultiProcess.R
In yuima: The YUIMA Project Package for SDEs
Defines functions
Multi.Euler
simCode
aux.simulate.multimodel
#
setMethod("simulate", "yuima.multimodel",
function(object, nsim=1, seed=NULL, xinit, true.parameter,
space.discretized=FALSE, increment.W=NULL, increment.L=NULL, method="euler",
hurst, methodfGn="WoodChan",
sampling, subsampling,
#Initial = 0, Terminal = 1, n = 100, delta,
# grid, random = FALSE, sdelta=as.numeric(NULL),
# sgrid=as.numeric(NULL), interpolation="none"
...){
tmpsamp <- NULL
if(missing(sampling)){
tmpsamp <- setSampling(...)
# tmpsamp <- setSampling(Initial = Initial, Terminal = Terminal, n = n,
# delta = delta, grid = grid, random = random, sdelta=sdelta,
# sgrid=sgrid, interpolation=interpolation)
} else {
tmpsamp <- sampling
}
tmpyuima <- setYuima(model=object, sampling=tmpsamp)
out <- simulate(tmpyuima, nsim=nsim, seed=seed, xinit=xinit,
true.parameter=true.parameter,
space.discretized=space.discretized,
increment.W=increment.W, increment.L=increment.L,
method=method,
hurst=hurst,methodfGn=methodfGn, subsampling=subsampling)
return(out)
})
aux.simulate.multimodel<-function(object, nsim, seed, xinit, true.parameter,
space.discretized, increment.W, increment.L,method,
hurst=0.5,methodfGn,
sampling, subsampling,
Initial, Terminal, n, delta,
grid, random, sdelta,
sgrid, interpolation){
##:: errors checks
if(is(object@model@measure$df,"yuima.law")&& is.null(increment.L)){
samp<- object@sampling
randomGenerator<-object@model@measure$df
if(samp@regular){
tForMeas<-samp@delta
NumbIncr<-samp@n[1]
if(missing(true.parameter)){
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}else{
measureparam<-true.parameter[object@model@parameter@measure]
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}
Noise<- rand(object = randomGenerator, n=NumbIncr, param=measureparam)
}else{
# Just For Irregular Grid
tForMeas<-diff(samp@grid[[1]]-samp@Initial)
my.InternalFunforLevy<-function(tForMeas,
randomGenerator,
true.parameter,object){
if(missing(true.parameter)){
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}else{
measureparam<-true.parameter[object@model@parameter@measure]
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}
Noise<- rand(object = randomGenerator, n=samp@n[1], param=measureparam)
}
Noise<-sapply(X=tForMeas, FUN=my.InternalFunforLevy,
randomGenerator=randomGenerator,
true.parameter=true.parameter,
object=object)
}
increment.L=Noise
if(is(object@model@measure$df,"yuima.law")&&!is.null(increment.L)){
dummy<-object
dummy@model@measure$df <- expression()
if(missing(xinit)){
res<- aux.simulate.multimodel(object=dummy, nsim, xinit=object@model@xinit ,seed, true.parameter,
space.discretized, increment.W,
increment.L = t(increment.L), method, hurst, methodfGn,
sampling=object@sampling)
}else{
res<- aux.simulate.multimodel(object=dummy, nsim, xinit, seed, true.parameter,
space.discretized, increment.W,
increment.L = increment.L, method, hurst, methodfGn,
sampling=object@sampling)
}
res@model <- object@model
if(missing(subsampling))
return(res)
return(subsampling(res, subsampling))
}
}
##:1: error on yuima model
yuima <- object
samp <- object@sampling
if(missing(yuima)){
yuima.warn("yuima object is missing.")
return(NULL)
}
tmphurst<-yuima@model@hurst
if(!missing(hurst)){
yuima@model@hurst=hurst
}
if (is.na(yuima@model@hurst)){
yuima.warn("Specify the hurst parameter.")
return(NULL)
}
tmpsamp <- NULL
if(is.null(yuima@sampling)){
if(missing(sampling)){
tmpsamp <- setSampling(Initial = Initial,
Terminal = Terminal, n = n, delta = delta,
grid = grid, random = random, sdelta=sdelta,
sgrid=sgrid, interpolation=interpolation)
} else {
tmpsamp <- sampling
}
} else {
tmpsamp <- yuima@sampling
}
yuima@sampling <- tmpsamp
sdeModel <- yuima@model
Terminal <- yuima@sampling@Terminal[1]
n <- yuima@sampling@n[1]
r.size <- sdeModel@noise.number
d.size <- sdeModel@equation.number
##:2: error on xinit
if(missing(xinit)){
xinit <- sdeModel@xinit
} else {
if(length(xinit) != d.size){
if(length(xinit)==1){
xinit <- rep(xinit, d.size)
} else {
yuima.warn("Dimension of xinit variables missmatch.")
return(NULL)
}
}
}
xinit <- as.expression(xinit) # force xinit to be an expression
par.len <- length(sdeModel@parameter@all)
if(missing(true.parameter) & par.len>0){
true.parameter <- vector(par.len, mode="list")
for(i in 1:par.len)
true.parameter[[i]] <- 0
names(true.parameter) <- sdeModel@parameter@all
}
yuimaEnv <- new.env()
if(par.len>0){
for(i in 1:par.len){
pars <- sdeModel@parameter@all[i]
for(j in 1:length(true.parameter)){
if( is.na(match(pars, names(true.parameter)[j]))!=TRUE){
assign(sdeModel@parameter@all[i], true.parameter[[j]],yuimaEnv)
}
}
#assign(sdeModel@parameter@all[i], true.parameter[[i]], yuimaEnv)
}
}
if(space.discretized){
if(r.size>1){
warning("Space-discretized EM cannot be used for multi-dimentional models. Use standard method.")
space.discretized <- FALSE
}
if(length(sdeModel@jump.coeff)){
warning("Space-discretized EM is for only Wiener Proc. Use standard method.")
space.discretized <- FALSE
}
}
##:: Error check for increment specified version.
if(!missing(increment.W) & !is.null(increment.W)){
if(space.discretized == TRUE){
yuima.warn("Parameter increment must be invalid if space.discretized=TRUE.")
return(NULL)
}else if(dim(increment.W)[1] != r.size){
yuima.warn("Length of increment's row must be same as yuima@model@noise.number.")
return(NULL)
}else if(dim(increment.W)[2] != n){
yuima.warn("Length of increment's column must be same as sampling@n[1].")
return(NULL)
}
}
##:: Error check for increment specified version.
if(!missing(increment.L) & !is.null(increment.L)){
if(space.discretized == TRUE){
yuima.warn("Parameter increment must be invalid if space.discretized=TRUE.")
return(NULL)
}else if(dim(increment.L)[1] != length(yuima@model@jump.coeff[[1]]) ){ #r.size){
yuima.warn("Length of increment's row must be same as yuima@model@noise.number.")
return(NULL)
}else if(dim(increment.L)[2] != n){
yuima.warn("Length of increment's column must be same as sampling@n[1].")
return(NULL)
}
}
yuimaEnv$dL <- increment.L
if(space.discretized){
##:: using Space-discretized Euler-Maruyama method
yuima@data <- space.discretized(xinit, yuima, yuimaEnv)
yuima@model@hurst<-tmphurst
return(yuima)
}
##:: using Euler-Maruyama method
delta <- samp@delta
if(missing(increment.W) | is.null(increment.W)){
if( sdeModel@hurst!=0.5 ){
grid<-sampling2grid(yuima@sampling)
isregular<-yuima@sampling@regular
if((!isregular) || (methodfGn=="Cholesky")){
dW<-CholeskyfGn(grid, sdeModel@hurst,r.size)
yuima.warn("Cholesky method for simulating fGn has been used.")
} else {
dW<-WoodChanfGn(grid, sdeModel@hurst,r.size)
}
} else {
delta<-samp@delta
if(!is.Poisson(sdeModel)){ # if pure CP no need to setup dW
dW <- rnorm(n * r.size, 0, sqrt(delta))
dW <- matrix(dW, ncol=n, nrow=r.size,byrow=TRUE)
} else {
dW <- matrix(0,ncol=n,nrow=1) # maybe to be fixed
}
}
} else {
dW <- increment.W
}
# if(is.Poisson(sdeModel)){
# yuima@data <- simCP(xinit, yuima, yuimaEnv)
# } else {
# yuima@data <- euler(xinit, yuima, dW, yuimaEnv)
# }
if(is(sdeModel,"yuima.multimodel")&&!is(sdeModel@measure$df,"yuima.law")){
if(length(sdeModel@measure.type)==1){
if(sdeModel@measure.type=="CP"){
intens <- as.character(sdeModel@measure$intensity)
dens <- as.character(sdeModel@measure$df$expr)
dumCP <- setPoisson(intensity = intens, df = dens,
dimension = length(sdeModel@jump.coeff[[1]]))
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = dummSamp@Terminal,
n = dummSamp@n)
traj <- simulate(object = dumCP,
sampling = samp,
true.parameter = true.parameter)
Incr.levy <- diff(traj@data@zoo.data[[1]])
if(length(traj@data@zoo.data)>1){
for(i in c(2:length(traj@data@zoo.data))){
Incr.levy<-cbind(Incr.levy,diff(traj@data@zoo.data[[i]]))
}
}
}else{
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = dummSamp@Terminal,
n = dummSamp@n)
xinitCode <- yuima@model@xinit
dimJumpCoeff <- length(yuima@model@jump.coeff[[1]])
dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
Dumsolve.variable<-paste0("MyLevyDum",c(1:dimJumpCoeff))
if(!is(sdeModel@measure$df,"yuima.law")){
LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
df = as.character(sdeModel@measure$df$expr),
measure.type = sdeModel@measure.type,
solve.variable = Dumsolve.variable)
}else{
LevyMod <- setModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
measure = sdeModel@measure,
measure.type = sdeModel@measure.type,
solve.variable = Dumsolve.variable)
}
yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
yuimaLevy@model@dimension <- dimJumpCoeff
traj<- simCode(xinit=xinitCode,yuima = yuimaLevy, env=yuimaEnv)
Incr.levy <- diff(traj@data@zoo.data[[1]])
if(length(traj@data@zoo.data)>1){
for(i in c(2:length(traj@data@zoo.data))){
Incr.levy<-cbind(Incr.levy,diff(traj@data@zoo.data[[i]]))
}
}
#yuima.stop("code multivariate Levy will be implemented as soon as possible")
}
}else{
if(any(sdeModel@measure.type=="CP")){
intens <- as.character(sdeModel@measure$intensity)
dens <- as.character(sdeModel@measure$df$expr)
# If we consider independence between CP and the Other Levy
# we have:
numbLev <- length(sdeModel@measure.type)
posCPindex <- c(1:numbLev)[sdeModel@measure.type%in%"CP"]
CPmeasureComp <- paste0(dens,"[,c(",toString(posCPindex),")]")
intens <- as.character(sdeModel@measure$intensity)
dumCP <- setPoisson(intensity = intens, df = CPmeasureComp,
dimension = length(posCPindex))
# Simulation CP part
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = unique(dummSamp@Terminal),
n = unique(dummSamp@n))
trajCP <- simulate(object = dumCP, sampling = samp,
true.parameter = true.parameter)
dimJumpCoeff <- length(yuima@model@measure.type)
dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
Dumsolve.variable <- paste0("MyLevyDum",c(1:dimJumpCoeff))
dummy.measure.code <- as.character(sdeModel@measure$df$expr)
LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
df = dummy.measure.code,
measure.type = "code",
solve.variable = Dumsolve.variable)
yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
yuimaLevy@model@dimension <- dimJumpCoeff
trajcode<- simCode(xinit=rep("0",length=dimJumpCoeff),
yuima = yuimaLevy, env=yuimaEnv)
countCP <- 0
countcode <- 0
if(yuima@model@measure.type[1]=="CP"){
Incr.levy <- as.matrix(as.numeric(diff(trajCP@data@zoo.data[[1]])))
countcode <- countcode+1
}else{
if(yuima@model@measure.type[1]=="code"){
Incr.levy <- as.matrix(as.numeric(diff(trajcode@data@zoo.data[[1]])))
countCP <- countCP+1
}
}
if(length(yuima@model@measure.type)>1){
for(i in c(2:length(yuima@model@measure.type))){
if(yuima@model@measure.type[i]=="CP"){
Incr.levy<-cbind(Incr.levy,as.numeric(diff(trajCP@data@zoo.data[[(i-countCP)]])))
countcode <- countcode+1
}else{
if(yuima@model@measure.type[i]=="code"){
Incr.levy <- cbind(Incr.levy,as.numeric(diff(trajcode@data@zoo.data[[i]])))
countCP <- countCP+1
}
}
}
}
# yuima.stop("Levy with CP and/or code")
}
}
}
if(!is.null(increment.L))
Incr.levy<-t(increment.L)
assign("dL",t(Incr.levy),envir=yuimaEnv)
sim <- Multi.Euler(xinit,yuima,dW,env=yuimaEnv)
yuima@data@zoo.data<-as.list(numeric(length=length(sim@zoo.data))) #LM nov2016
# yuima@data@zoo.data<-sim@zoo.data
for(i in 1:length(yuima@data@zoo.data)){
yuima@data@zoo.data[[i]]<-sim@zoo.data[[i]]
index(yuima@data@zoo.data[[i]]) <- yuima@sampling@grid[[1]]
}## to be fixed
yuima@data@original.data <- sim@original.data
yuima@model@xinit <- xinit
yuima@model@hurst <-tmphurst
if(missing(subsampling))
return(yuima)
subsampling(yuima, subsampling)
}
simCode <- function(xinit,yuima,env){
sdeModel<-yuima@model
modelstate <- sdeModel@solve.variable
modeltime <- sdeModel@time.variable
Terminal <- yuima@sampling@Terminal[1]
Initial <- yuima@sampling@Initial[1]
dimension <- yuima@model@dimension
dummy.val <- numeric(dimension)
if(length(xinit) != dimension)
xinit <- rep(xinit, dimension)[1:dimension]
if(length(unique(as.character(xinit)))==1 &&
is.numeric(tryCatch(eval(xinit[1],envir=env),error=function(...) FALSE))){
dX_dummy<-xinit[1]
dummy.val<-eval(dX_dummy, envir=env)
if(length(dummy.val)==1){
dummy.val<-rep(dummy.val,dimension)
}
for(i in 1:length(modelstate)){
assign(modelstate[i],dummy.val[i] ,envir=env)
}
} else {
for(i in 1:dimension){
dummy.val[i] <- eval(xinit[i], envir=env)
}
}
### Simulation of CP using Lewis' method
##:: Levy
JP <- eval(sdeModel@jump.coeff[[1]], envir=env)
mu.size <- length(JP)
# print(str(JP))
#assign(sdeModel@measure$intensity, env) ## intensity param
# .CPintensity <- function(.t) {
# assign(modeltime, .t, envir=env)
# eval(sdeModel@measure$intensity, envir=env)
# }
dummyList<-as.list(env)
lgth.meas<-length(yuima@model@parameter@measure)
if(lgth.meas>1){
for(i in c(2:lgth.meas)){
idx.dummy<-yuima@model@parameter@measure[i]
assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
}
}
# we use Lewis' acceptance/rejection method
#if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel@measure$df$expr)){
##:: e.g. dnorm(z,1,1) -> rnorm(mu.size*N_sharp,1,1)
#gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel@measure$df$expr, perl=TRUE)
#} else{
#stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
#}
if(!is(sdeModel@measure$df,"yuima.law")){
dumStringMeas <- toString(sdeModel@measure$df$expr)
dumStringMeas1 <- substr(x=dumStringMeas, start=2,stop=nchar(x = dumStringMeas))
dumStringMeas2 <- paste0("r",dumStringMeas1)
tmpMeas2 <- strsplit(x=dumStringMeas2,split="")
posMeas2 <- match("(" , tmpMeas2[[1]])[1]
dumStringMeas3 <- substr(x=dumStringMeas2, start=1,stop=(posMeas2-1))
a<-deparse(args(eval(parse(text=dumStringMeas3))))[1]
b<-gsub("^function (.+?)","(",a)
b1 <- substr(x=b,start =1, stop=(nchar(b)-1))
FinalMeasRandn<-paste0(dumStringMeas3,b1)
dummyvarMaes <- all.vars(parse(text=FinalMeasRandn))
posDum<- match(c(sdeModel@jump.variable,sdeModel@parameter@measure),dummyvarMaes)
if(length(posDum)+1!=length(dummyvarMaes)){
yuima.stop("too input variables in the random number function")
}
deltaVar <- dummyvarMaes[-posDum]
# ell <- optimize(f=.CPintensity, interval=c(Initial, Terminal), maximum = TRUE)$objective
# ellMax <- ell * 1.01
F <- suppressWarnings(parse(text=gsub("^r(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", parse(text=FinalMeasRandn), perl=TRUE)))
F.env <- new.env(parent=env)
N_sharp <- unique(yuima@sampling@n)
TrueDelta <- unique(yuima@sampling@delta)
assign(deltaVar, TrueDelta, envir=F.env)
assign("mu.size", mu.size, envir=F.env)
assign("N_sharp", N_sharp, envir=F.env)
randJ <- eval(F, envir=F.env) ## this expression is evaluated in the F.env
randJ <- rbind(dummy.val, randJ)
}else{
TrueDelta <- unique(yuima@sampling@delta)
randJ<- env$dL
}
RANDLevy <- apply(randJ,2,cumsum)
tsX <- zoo(x=RANDLevy)
yuimaData <- setYuima(data=setData(tsX, delta=TrueDelta))
#yuimaData <- subsampling(yuimaData, sampling=yuima@sampling)
return(yuimaData)
}
Multi.Euler<-function(xinit,yuima,dW,env){
sdeModel<-yuima@model
modelstate <- sdeModel@solve.variable
modeltime <- sdeModel@time.variable
V0 <- sdeModel@drift
V <- sdeModel@diffusion
r.size <- sdeModel@noise.number
d.size <- sdeModel@equation.number
Terminal <- yuima@sampling@Terminal[1]
n <- yuima@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
##:: check if DRIFT and/or DIFFUSION has values
has.drift <- sum(as.character(sdeModel@drift) != "(0)")
var.in.diff <- is.logical(any(match(unlist(lapply(sdeModel@diffusion, all.vars)), sdeModel@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@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]
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)
}else{ ##:: Levy
JP <- sdeModel@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@measure.type == "CP" ){ ##:: Compound-Poisson type
#
# ##:: delete 2010/09/13 for simulate func bug fix by s.h
# ## eta0 <- eval(sdeModel@measure$intensity)
#
# ##:: add 2010/09/13 for simulate func bug fix by s.h
# eta0 <- eval(sdeModel@measure$intensity, env) ## intensity param
#
# ##:: get lambda from nu()
# tmp.expr <- function(my.x){
# assign(sdeModel@jump.variable,my.x)
# return(eval(sdeModel@measure$df$expr))
# }
# #lambda <- integrate(sdeModel@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@model@parameter@measure)
# if(lgth.meas>1){
# for(i in c(2:lgth.meas)){
# idx.dummy<-yuima@model@parameter@measure[i]
# #print(i)
# #print(yuima@model@parameter@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)
# if(N_sharp == 0){
# JAMP <- FALSE
# }else{
# JAMP <- TRUE
# Uj <- sort( runif(N_sharp, 0, Terminal) )
# ij <- NULL
# for(i in 1:length(Uj)){
# 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@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@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@jump.variable, J, env)
#
# if(sdeModel@J.flag){
# J <- 1
# }
#
# Pi <- p.b.j(t=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]
# }
# ##Pi <- p.b.j(t=i*delta, X=dX)
# }
# dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi
# X_mat[, i+1] <- dX
# }
# tsX <- ts(data=t(X_mat), deltat=delta, start=0)
# ##::end CP
# }else if(sdeModel@measure.type=="code"){ ##:: code type
# ##:: Jump terms
# code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel@measure$df$expr, perl=TRUE))
# args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel@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], ")"),
# rngamma=paste("rngamma(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)")
# )
# dummyList<-as.list(env)
# #print(str(dummyList))
# lgth.meas<-length(yuima@model@parameter@measure)
# #print(lgth.meas)
# if(lgth.meas!=0){
# for(i in c(1:lgth.meas)){
# #print(i)
# #print(yuima@model@parameter@measure[i])
# idx.dummy<-yuima@model@parameter@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@jump.variable))
for(i in 1:n){
assign(sdeModel@jump.variable, dZ[,i], env)
if(sdeModel@J.flag){
dZ[,i] <- 1
}
# cat("\np.b.j call\n")
tmp.j <- p.b.j(t=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)))
#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]
X_mat[, i+1] <- dX
}
tsX <- ts(data=t(X_mat), deltat=delta, start=0)
##::end code
# }else{
# cat(paste("Type \"", sdeModel@measure.type, "\" not supported yet.\n", sep=""))
# return(NULL)
# }
}##::end levy
yuimaData <- setData(original.data=tsX)
return(yuimaData)
}
Try the yuima package in your browser
Any scripts or data that you put into this service are public.
yuima documentation built on Nov. 14, 2022, 3:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Multi.Euler simCode aux.simulate.multimodel
#
setMethod("simulate", "yuima.multimodel",
function(object, nsim=1, seed=NULL, xinit, true.parameter,
space.discretized=FALSE, increment.W=NULL, increment.L=NULL, method="euler",
hurst, methodfGn="WoodChan",
sampling, subsampling,
#Initial = 0, Terminal = 1, n = 100, delta,
# grid, random = FALSE, sdelta=as.numeric(NULL),
# sgrid=as.numeric(NULL), interpolation="none"
...){
tmpsamp <- NULL
if(missing(sampling)){
tmpsamp <- setSampling(...)
# tmpsamp <- setSampling(Initial = Initial, Terminal = Terminal, n = n,
# delta = delta, grid = grid, random = random, sdelta=sdelta,
# sgrid=sgrid, interpolation=interpolation)
} else {
tmpsamp <- sampling
}
tmpyuima <- setYuima(model=object, sampling=tmpsamp)
out <- simulate(tmpyuima, nsim=nsim, seed=seed, xinit=xinit,
true.parameter=true.parameter,
space.discretized=space.discretized,
increment.W=increment.W, increment.L=increment.L,
method=method,
hurst=hurst,methodfGn=methodfGn, subsampling=subsampling)
return(out)
})
aux.simulate.multimodel<-function(object, nsim, seed, xinit, true.parameter,
space.discretized, increment.W, increment.L,method,
hurst=0.5,methodfGn,
sampling, subsampling,
Initial, Terminal, n, delta,
grid, random, sdelta,
sgrid, interpolation){
##:: errors checks
if(is(object@model@measure$df,"yuima.law")&& is.null(increment.L)){
samp<- object@sampling
randomGenerator<-object@model@measure$df
if(samp@regular){
tForMeas<-samp@delta
NumbIncr<-samp@n[1]
if(missing(true.parameter)){
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}else{
measureparam<-true.parameter[object@model@parameter@measure]
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}
Noise<- rand(object = randomGenerator, n=NumbIncr, param=measureparam)
}else{
# Just For Irregular Grid
tForMeas<-diff(samp@grid[[1]]-samp@Initial)
my.InternalFunforLevy<-function(tForMeas,
randomGenerator,
true.parameter,object){
if(missing(true.parameter)){
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}else{
measureparam<-true.parameter[object@model@parameter@measure]
eval(parse(text= paste0("measureparam$",
object@model@time.variable," <- tForMeas",collapse="")))
}
Noise<- rand(object = randomGenerator, n=samp@n[1], param=measureparam)
}
Noise<-sapply(X=tForMeas, FUN=my.InternalFunforLevy,
randomGenerator=randomGenerator,
true.parameter=true.parameter,
object=object)
}
increment.L=Noise
if(is(object@model@measure$df,"yuima.law")&&!is.null(increment.L)){
dummy<-object
dummy@model@measure$df <- expression()
if(missing(xinit)){
res<- aux.simulate.multimodel(object=dummy, nsim, xinit=object@model@xinit ,seed, true.parameter,
space.discretized, increment.W,
increment.L = t(increment.L), method, hurst, methodfGn,
sampling=object@sampling)
}else{
res<- aux.simulate.multimodel(object=dummy, nsim, xinit, seed, true.parameter,
space.discretized, increment.W,
increment.L = increment.L, method, hurst, methodfGn,
sampling=object@sampling)
}
res@model <- object@model
if(missing(subsampling))
return(res)
return(subsampling(res, subsampling))
}
}
##:1: error on yuima model
yuima <- object
samp <- object@sampling
if(missing(yuima)){
yuima.warn("yuima object is missing.")
return(NULL)
}
tmphurst<-yuima@model@hurst
if(!missing(hurst)){
yuima@model@hurst=hurst
}
if (is.na(yuima@model@hurst)){
yuima.warn("Specify the hurst parameter.")
return(NULL)
}
tmpsamp <- NULL
if(is.null(yuima@sampling)){
if(missing(sampling)){
tmpsamp <- setSampling(Initial = Initial,
Terminal = Terminal, n = n, delta = delta,
grid = grid, random = random, sdelta=sdelta,
sgrid=sgrid, interpolation=interpolation)
} else {
tmpsamp <- sampling
}
} else {
tmpsamp <- yuima@sampling
}
yuima@sampling <- tmpsamp
sdeModel <- yuima@model
Terminal <- yuima@sampling@Terminal[1]
n <- yuima@sampling@n[1]
r.size <- sdeModel@noise.number
d.size <- sdeModel@equation.number
##:2: error on xinit
if(missing(xinit)){
xinit <- sdeModel@xinit
} else {
if(length(xinit) != d.size){
if(length(xinit)==1){
xinit <- rep(xinit, d.size)
} else {
yuima.warn("Dimension of xinit variables missmatch.")
return(NULL)
}
}
}
xinit <- as.expression(xinit) # force xinit to be an expression
par.len <- length(sdeModel@parameter@all)
if(missing(true.parameter) & par.len>0){
true.parameter <- vector(par.len, mode="list")
for(i in 1:par.len)
true.parameter[[i]] <- 0
names(true.parameter) <- sdeModel@parameter@all
}
yuimaEnv <- new.env()
if(par.len>0){
for(i in 1:par.len){
pars <- sdeModel@parameter@all[i]
for(j in 1:length(true.parameter)){
if( is.na(match(pars, names(true.parameter)[j]))!=TRUE){
assign(sdeModel@parameter@all[i], true.parameter[[j]],yuimaEnv)
}
}
#assign(sdeModel@parameter@all[i], true.parameter[[i]], yuimaEnv)
}
}
if(space.discretized){
if(r.size>1){
warning("Space-discretized EM cannot be used for multi-dimentional models. Use standard method.")
space.discretized <- FALSE
}
if(length(sdeModel@jump.coeff)){
warning("Space-discretized EM is for only Wiener Proc. Use standard method.")
space.discretized <- FALSE
}
}
##:: Error check for increment specified version.
if(!missing(increment.W) & !is.null(increment.W)){
if(space.discretized == TRUE){
yuima.warn("Parameter increment must be invalid if space.discretized=TRUE.")
return(NULL)
}else if(dim(increment.W)[1] != r.size){
yuima.warn("Length of increment's row must be same as yuima@model@noise.number.")
return(NULL)
}else if(dim(increment.W)[2] != n){
yuima.warn("Length of increment's column must be same as sampling@n[1].")
return(NULL)
}
}
##:: Error check for increment specified version.
if(!missing(increment.L) & !is.null(increment.L)){
if(space.discretized == TRUE){
yuima.warn("Parameter increment must be invalid if space.discretized=TRUE.")
return(NULL)
}else if(dim(increment.L)[1] != length(yuima@model@jump.coeff[[1]]) ){ #r.size){
yuima.warn("Length of increment's row must be same as yuima@model@noise.number.")
return(NULL)
}else if(dim(increment.L)[2] != n){
yuima.warn("Length of increment's column must be same as sampling@n[1].")
return(NULL)
}
}
yuimaEnv$dL <- increment.L
if(space.discretized){
##:: using Space-discretized Euler-Maruyama method
yuima@data <- space.discretized(xinit, yuima, yuimaEnv)
yuima@model@hurst<-tmphurst
return(yuima)
}
##:: using Euler-Maruyama method
delta <- samp@delta
if(missing(increment.W) | is.null(increment.W)){
if( sdeModel@hurst!=0.5 ){
grid<-sampling2grid(yuima@sampling)
isregular<-yuima@sampling@regular
if((!isregular) || (methodfGn=="Cholesky")){
dW<-CholeskyfGn(grid, sdeModel@hurst,r.size)
yuima.warn("Cholesky method for simulating fGn has been used.")
} else {
dW<-WoodChanfGn(grid, sdeModel@hurst,r.size)
}
} else {
delta<-samp@delta
if(!is.Poisson(sdeModel)){ # if pure CP no need to setup dW
dW <- rnorm(n * r.size, 0, sqrt(delta))
dW <- matrix(dW, ncol=n, nrow=r.size,byrow=TRUE)
} else {
dW <- matrix(0,ncol=n,nrow=1) # maybe to be fixed
}
}
} else {
dW <- increment.W
}
# if(is.Poisson(sdeModel)){
# yuima@data <- simCP(xinit, yuima, yuimaEnv)
# } else {
# yuima@data <- euler(xinit, yuima, dW, yuimaEnv)
# }
if(is(sdeModel,"yuima.multimodel")&&!is(sdeModel@measure$df,"yuima.law")){
if(length(sdeModel@measure.type)==1){
if(sdeModel@measure.type=="CP"){
intens <- as.character(sdeModel@measure$intensity)
dens <- as.character(sdeModel@measure$df$expr)
dumCP <- setPoisson(intensity = intens, df = dens,
dimension = length(sdeModel@jump.coeff[[1]]))
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = dummSamp@Terminal,
n = dummSamp@n)
traj <- simulate(object = dumCP,
sampling = samp,
true.parameter = true.parameter)
Incr.levy <- diff(traj@data@zoo.data[[1]])
if(length(traj@data@zoo.data)>1){
for(i in c(2:length(traj@data@zoo.data))){
Incr.levy<-cbind(Incr.levy,diff(traj@data@zoo.data[[i]]))
}
}
}else{
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = dummSamp@Terminal,
n = dummSamp@n)
xinitCode <- yuima@model@xinit
dimJumpCoeff <- length(yuima@model@jump.coeff[[1]])
dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
Dumsolve.variable<-paste0("MyLevyDum",c(1:dimJumpCoeff))
if(!is(sdeModel@measure$df,"yuima.law")){
LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
df = as.character(sdeModel@measure$df$expr),
measure.type = sdeModel@measure.type,
solve.variable = Dumsolve.variable)
}else{
LevyMod <- setModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
measure = sdeModel@measure,
measure.type = sdeModel@measure.type,
solve.variable = Dumsolve.variable)
}
yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
yuimaLevy@model@dimension <- dimJumpCoeff
traj<- simCode(xinit=xinitCode,yuima = yuimaLevy, env=yuimaEnv)
Incr.levy <- diff(traj@data@zoo.data[[1]])
if(length(traj@data@zoo.data)>1){
for(i in c(2:length(traj@data@zoo.data))){
Incr.levy<-cbind(Incr.levy,diff(traj@data@zoo.data[[i]]))
}
}
#yuima.stop("code multivariate Levy will be implemented as soon as possible")
}
}else{
if(any(sdeModel@measure.type=="CP")){
intens <- as.character(sdeModel@measure$intensity)
dens <- as.character(sdeModel@measure$df$expr)
# If we consider independence between CP and the Other Levy
# we have:
numbLev <- length(sdeModel@measure.type)
posCPindex <- c(1:numbLev)[sdeModel@measure.type%in%"CP"]
CPmeasureComp <- paste0(dens,"[,c(",toString(posCPindex),")]")
intens <- as.character(sdeModel@measure$intensity)
dumCP <- setPoisson(intensity = intens, df = CPmeasureComp,
dimension = length(posCPindex))
# Simulation CP part
dummSamp <- yuima@sampling
samp <- setSampling(Initial = dummSamp@Initial,
Terminal = unique(dummSamp@Terminal),
n = unique(dummSamp@n))
trajCP <- simulate(object = dumCP, sampling = samp,
true.parameter = true.parameter)
dimJumpCoeff <- length(yuima@model@measure.type)
dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
Dumsolve.variable <- paste0("MyLevyDum",c(1:dimJumpCoeff))
dummy.measure.code <- as.character(sdeModel@measure$df$expr)
LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
df = dummy.measure.code,
measure.type = "code",
solve.variable = Dumsolve.variable)
yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
yuimaLevy@model@dimension <- dimJumpCoeff
trajcode<- simCode(xinit=rep("0",length=dimJumpCoeff),
yuima = yuimaLevy, env=yuimaEnv)
countCP <- 0
countcode <- 0
if(yuima@model@measure.type[1]=="CP"){
Incr.levy <- as.matrix(as.numeric(diff(trajCP@data@zoo.data[[1]])))
countcode <- countcode+1
}else{
if(yuima@model@measure.type[1]=="code"){
Incr.levy <- as.matrix(as.numeric(diff(trajcode@data@zoo.data[[1]])))
countCP <- countCP+1
}
}
if(length(yuima@model@measure.type)>1){
for(i in c(2:length(yuima@model@measure.type))){
if(yuima@model@measure.type[i]=="CP"){
Incr.levy<-cbind(Incr.levy,as.numeric(diff(trajCP@data@zoo.data[[(i-countCP)]])))
countcode <- countcode+1
}else{
if(yuima@model@measure.type[i]=="code"){
Incr.levy <- cbind(Incr.levy,as.numeric(diff(trajcode@data@zoo.data[[i]])))
countCP <- countCP+1
}
}
}
}
# yuima.stop("Levy with CP and/or code")
}
}
}
if(!is.null(increment.L))
Incr.levy<-t(increment.L)
assign("dL",t(Incr.levy),envir=yuimaEnv)
sim <- Multi.Euler(xinit,yuima,dW,env=yuimaEnv)
yuima@data@zoo.data<-as.list(numeric(length=length(sim@zoo.data))) #LM nov2016
# yuima@data@zoo.data<-sim@zoo.data
for(i in 1:length(yuima@data@zoo.data)){
yuima@data@zoo.data[[i]]<-sim@zoo.data[[i]]
index(yuima@data@zoo.data[[i]]) <- yuima@sampling@grid[[1]]
}## to be fixed
yuima@data@original.data <- sim@original.data
yuima@model@xinit <- xinit
yuima@model@hurst <-tmphurst
if(missing(subsampling))
return(yuima)
subsampling(yuima, subsampling)
}
simCode <- function(xinit,yuima,env){
sdeModel<-yuima@model
modelstate <- sdeModel@solve.variable
modeltime <- sdeModel@time.variable
Terminal <- yuima@sampling@Terminal[1]
Initial <- yuima@sampling@Initial[1]
dimension <- yuima@model@dimension
dummy.val <- numeric(dimension)
if(length(xinit) != dimension)
xinit <- rep(xinit, dimension)[1:dimension]
if(length(unique(as.character(xinit)))==1 &&
is.numeric(tryCatch(eval(xinit[1],envir=env),error=function(...) FALSE))){
dX_dummy<-xinit[1]
dummy.val<-eval(dX_dummy, envir=env)
if(length(dummy.val)==1){
dummy.val<-rep(dummy.val,dimension)
}
for(i in 1:length(modelstate)){
assign(modelstate[i],dummy.val[i] ,envir=env)
}
} else {
for(i in 1:dimension){
dummy.val[i] <- eval(xinit[i], envir=env)
}
}
### Simulation of CP using Lewis' method
##:: Levy
JP <- eval(sdeModel@jump.coeff[[1]], envir=env)
mu.size <- length(JP)
# print(str(JP))
#assign(sdeModel@measure$intensity, env) ## intensity param
# .CPintensity <- function(.t) {
# assign(modeltime, .t, envir=env)
# eval(sdeModel@measure$intensity, envir=env)
# }
dummyList<-as.list(env)
lgth.meas<-length(yuima@model@parameter@measure)
if(lgth.meas>1){
for(i in c(2:lgth.meas)){
idx.dummy<-yuima@model@parameter@measure[i]
assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
}
}
# we use Lewis' acceptance/rejection method
#if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel@measure$df$expr)){
##:: e.g. dnorm(z,1,1) -> rnorm(mu.size*N_sharp,1,1)
#gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel@measure$df$expr, perl=TRUE)
#} else{
#stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
#}
if(!is(sdeModel@measure$df,"yuima.law")){
dumStringMeas <- toString(sdeModel@measure$df$expr)
dumStringMeas1 <- substr(x=dumStringMeas, start=2,stop=nchar(x = dumStringMeas))
dumStringMeas2 <- paste0("r",dumStringMeas1)
tmpMeas2 <- strsplit(x=dumStringMeas2,split="")
posMeas2 <- match("(" , tmpMeas2[[1]])[1]
dumStringMeas3 <- substr(x=dumStringMeas2, start=1,stop=(posMeas2-1))
a<-deparse(args(eval(parse(text=dumStringMeas3))))[1]
b<-gsub("^function (.+?)","(",a)
b1 <- substr(x=b,start =1, stop=(nchar(b)-1))
FinalMeasRandn<-paste0(dumStringMeas3,b1)
dummyvarMaes <- all.vars(parse(text=FinalMeasRandn))
posDum<- match(c(sdeModel@jump.variable,sdeModel@parameter@measure),dummyvarMaes)
if(length(posDum)+1!=length(dummyvarMaes)){
yuima.stop("too input variables in the random number function")
}
deltaVar <- dummyvarMaes[-posDum]
# ell <- optimize(f=.CPintensity, interval=c(Initial, Terminal), maximum = TRUE)$objective
# ellMax <- ell * 1.01
F <- suppressWarnings(parse(text=gsub("^r(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", parse(text=FinalMeasRandn), perl=TRUE)))
F.env <- new.env(parent=env)
N_sharp <- unique(yuima@sampling@n)
TrueDelta <- unique(yuima@sampling@delta)
assign(deltaVar, TrueDelta, envir=F.env)
assign("mu.size", mu.size, envir=F.env)
assign("N_sharp", N_sharp, envir=F.env)
randJ <- eval(F, envir=F.env) ## this expression is evaluated in the F.env
randJ <- rbind(dummy.val, randJ)
}else{
TrueDelta <- unique(yuima@sampling@delta)
randJ<- env$dL
}
RANDLevy <- apply(randJ,2,cumsum)
tsX <- zoo(x=RANDLevy)
yuimaData <- setYuima(data=setData(tsX, delta=TrueDelta))
#yuimaData <- subsampling(yuimaData, sampling=yuima@sampling)
return(yuimaData)
}
Multi.Euler<-function(xinit,yuima,dW,env){
sdeModel<-yuima@model
modelstate <- sdeModel@solve.variable
modeltime <- sdeModel@time.variable
V0 <- sdeModel@drift
V <- sdeModel@diffusion
r.size <- sdeModel@noise.number
d.size <- sdeModel@equation.number
Terminal <- yuima@sampling@Terminal[1]
n <- yuima@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
##:: check if DRIFT and/or DIFFUSION has values
has.drift <- sum(as.character(sdeModel@drift) != "(0)")
var.in.diff <- is.logical(any(match(unlist(lapply(sdeModel@diffusion, all.vars)), sdeModel@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@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]
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)
}else{ ##:: Levy
JP <- sdeModel@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@measure.type == "CP" ){ ##:: Compound-Poisson type
#
# ##:: delete 2010/09/13 for simulate func bug fix by s.h
# ## eta0 <- eval(sdeModel@measure$intensity)
#
# ##:: add 2010/09/13 for simulate func bug fix by s.h
# eta0 <- eval(sdeModel@measure$intensity, env) ## intensity param
#
# ##:: get lambda from nu()
# tmp.expr <- function(my.x){
# assign(sdeModel@jump.variable,my.x)
# return(eval(sdeModel@measure$df$expr))
# }
# #lambda <- integrate(sdeModel@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@model@parameter@measure)
# if(lgth.meas>1){
# for(i in c(2:lgth.meas)){
# idx.dummy<-yuima@model@parameter@measure[i]
# #print(i)
# #print(yuima@model@parameter@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)
# if(N_sharp == 0){
# JAMP <- FALSE
# }else{
# JAMP <- TRUE
# Uj <- sort( runif(N_sharp, 0, Terminal) )
# ij <- NULL
# for(i in 1:length(Uj)){
# 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@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@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@jump.variable, J, env)
#
# if(sdeModel@J.flag){
# J <- 1
# }
#
# Pi <- p.b.j(t=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]
# }
# ##Pi <- p.b.j(t=i*delta, X=dX)
# }
# dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi
# X_mat[, i+1] <- dX
# }
# tsX <- ts(data=t(X_mat), deltat=delta, start=0)
# ##::end CP
# }else if(sdeModel@measure.type=="code"){ ##:: code type
# ##:: Jump terms
# code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel@measure$df$expr, perl=TRUE))
# args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel@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], ")"),
# rngamma=paste("rngamma(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)")
# )
# dummyList<-as.list(env)
# #print(str(dummyList))
# lgth.meas<-length(yuima@model@parameter@measure)
# #print(lgth.meas)
# if(lgth.meas!=0){
# for(i in c(1:lgth.meas)){
# #print(i)
# #print(yuima@model@parameter@measure[i])
# idx.dummy<-yuima@model@parameter@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@jump.variable))
for(i in 1:n){
assign(sdeModel@jump.variable, dZ[,i], env)
if(sdeModel@J.flag){
dZ[,i] <- 1
}
# cat("\np.b.j call\n")
tmp.j <- p.b.j(t=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)))
#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]
X_mat[, i+1] <- dX
}
tsX <- ts(data=t(X_mat), deltat=delta, start=0)
##::end code
# }else{
# cat(paste("Type \"", sdeModel@measure.type, "\" not supported yet.\n", sep=""))
# return(NULL)
# }
}##::end levy
yuimaData <- setData(original.data=tsX)
return(yuimaData)
}
Try the yuima package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.