[Yuima-commits] r511 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 8 20:25:18 CET 2016
Author: lorenzo
Date: 2016-11-08 20:25:17 +0100 (Tue, 08 Nov 2016)
New Revision: 511
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/AuxMethodforPPR.R
pkg/yuima/R/setMultiModel.R
pkg/yuima/R/simulateForMapsIntegralAndOperator.R
pkg/yuima/R/simulateForPpr.R
pkg/yuima/R/simulateMultiProcess.R
Log:
Ppr Model
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/DESCRIPTION 2016-11-08 19:25:17 UTC (rev 511)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.3.6
+Version: 1.3.7
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/R/AuxMethodforPPR.R 2016-11-08 19:25:17 UTC (rev 511)
@@ -21,14 +21,16 @@
}
}
for(j in c(1:length(yuimaPpr at Ppr@counting.var))){
- cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@counting.var[[j]]
+ #cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@counting.var[[j]]
+ cond<-yuimaPpr at model@solve.variable%in%yuimaPpr at Ppr@counting.var[[j]]
assign(yuimaPpr at Ppr@counting.var[[j]],
as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
envir=envPpr[[i]])
}
for(j in c(1:length(yuimaPpr at Ppr@var.dx))){
- cond<-c(colnames(yuimaPpr at data@original.data),yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
+ #cond<-c(colnames(yuimaPpr at data@original.data),yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
+ cond<-c(yuimaPpr at model@solve.variable,yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
if(any(cond[-length(cond)])){
assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
diff(as.numeric(yuimaPpr at data@original.data[1:(i+1),cond[-length(cond)]])),
@@ -36,7 +38,7 @@
}
if(tail(cond,n=1L)){
assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
- as.numeric(diff(Time[1:(i+1),cond[-length(cond)]])),
+ as.numeric(diff(Time[1:(i+1)])),
envir=envPpr[[i]])
}
}
@@ -58,6 +60,11 @@
}
assign("Integrator",Integrator,envir=envPpr[[length(envPpr)]])
assign("Nlamb",length(yuimaPpr at Ppr@counting.var),envir=envPpr[[length(envPpr)]])
+ cond <- yuimaPpr at model@solve.variable %in% yuimaPpr at Ppr@counting.var
+ #assign("CountVar",yuimaPpr at data@original.data[,cond],envir=envPpr[[length(envPpr)]])
+ DummyCountVar <- yuimaPpr at data@original.data[,cond]
+ condDummyCountVar <- t(apply(as.matrix(DummyCountVar),FUN = "diff",MARGIN = 2))!=0
+ assign("CountVar", condDummyCountVar, envir=envPpr[[length(envPpr)]])
out<-NULL
param1<-unlist(parLambda)
@@ -233,7 +240,8 @@
# nrow=1,byrow=TRUE)
IntKer<- matrix(unlist(lapply(envPpr,myfun3,Kern=Kern)),
nrow=lastEnv$Nlamb)
- lambda <- gFunVect+cbind(0,IntKer)
+ # lambda <- gFunVect+cbind(0,IntKer)
+ lambda <- gFunVect+IntKer
time <- (c(lastEnv$s,lastEnv$t[1]))
if(!logLikelihood){
Intensity <- zoo(t(lambda), order.by = time)
@@ -241,14 +249,29 @@
}
dn <- dim(lambda)
if(dn[1]==1){
- logLiklihood2 <- -sum(lambda[,-1]*diff(time)[1])
- logLiklihood1 <- sum(log(lambda[,-1])*lastEnv$Integrator)
+ # logLiklihood2 <- -sum(lambda[,-1]*diff(time)[1])
+ # logLiklihood1 <- sum(log(lambda[,-1])*lastEnv$Integrator)
+ logLiklihood2 <- -sum(lambda*diff(time)[1])
+ # CountProc <- diff(as.numeric(lastEnv$CountVar))!=0
+ logLiklihood1 <- sum(log(lambda)*lastEnv$CountVar)
+
+
}else{
- logLiklihood2 <- -rowSums(lambda[,-1]*diff(time)[1])
- logLiklihood1 <- rowSums(log(lambda[,-1])*lastEnv$Integrator)
+ # logLiklihood2 <- -rowSums(lambda[,-1]*diff(time)[1])
+ # logLiklihood1 <- rowSums(log(lambda[,-1])*lastEnv$Integrator)
+ logLiklihood2 <- -rowSums(lambda*diff(time)[1])
+ #cond <- t(apply(as.matrix(lastEnv$CountVar),FUN = "diff",MARGIN = 2))!=0
+ logLiklihood1 <- rowSums(log(lambda)*lastEnv$CountVar)
}
+ if(is.nan(logLiklihood1)){
+ logLiklihood1 <- -10^10
+ }
+ if(is.nan(logLiklihood2)){
+ logLiklihood2 <- -10^10
+ }
minusLoglik <- -sum(logLiklihood2+logLiklihood1)
- #cat(sprintf("\n%.5f",minusLoglik))
+ # cat(sprintf("\n%.5f",minusLoglik))
+ # cat(sprintf("\n%.5f",param))
return(minusLoglik)
}
Modified: pkg/yuima/R/setMultiModel.R
===================================================================
--- pkg/yuima/R/setMultiModel.R 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/R/setMultiModel.R 2016-11-08 19:25:17 UTC (rev 511)
@@ -136,7 +136,8 @@
if(existsFunction(measurefunc)){
tmp.measure$df$func[[1]] <- eval(parse(text=measurefunc))
}else{
- yuima.stop("density function for jump must be specified")
+ if(measurefunc!="yuima.law")
+ yuima.stop("density function for jump must be specified")
}
@@ -327,10 +328,10 @@
# n.eqn3 <- n.eqn1
# }
- if(n.eqn1 != n.eqn2 || n.eqn1 != n.eqn3){
- yuima.warn("Malformed model, number of equations in the drift and diffusion do not match.")
- return(NULL)
- }
+ if(n.eqn1 != n.eqn2 || n.eqn1 != n.eqn3){
+ yuima.warn("Malformed model, number of equations in the drift and diffusion do not match.")
+ return(NULL)
+ }
n.eqn <- n.eqn1
if(is.null(xinit)){
Modified: pkg/yuima/R/simulateForMapsIntegralAndOperator.R
===================================================================
--- pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-11-08 19:25:17 UTC (rev 511)
@@ -143,48 +143,50 @@
nrow = info.int at dimIntegrand[1], ncol = info.int at dimIntegrand[2])
my.fun <- function(my.list, my.env, i){
dum <- eval(my.list,envir = my.env)
- if(length(dum)==1){
- return(rep(dum,i))
- }else{
+ #if(length(dum)==1){
+ # return(rep(dum,i))
+ #}else{
return(dum[1:i])
- }
+ #}
}
- res<-matrix(0,info.int at dimIntegrand[1],(length(time)-1))
+ # res<-matrix(0,info.int at dimIntegrand[1],(length(time)-1))
+ # element <- matrix(0, nrow =info.int at dimIntegrand[1], ncol = 1)
+ #
+ # for(i in c(1:(length(time)-1))){
+ # assign(info.var at upper.var,time[i+1],envir=my.env)
+ # Inter2 <- lapply(info.int at IntegrandList,
+ # FUN = my.fun, my.env = my.env, i = i)
+ # for(j in c(1:info.int at dimIntegrand[1])){
+ # element[j,] <- M.dX[,c(1:i)]%*%matrix(unlist(Inter2[PosInMatr[j,]]),
+ # ncol = info.int at dimIntegrand[2])
+ # }
+ # res[,i] <- element
+ # }
+ LengTime<-(length(time)-1)
+ my.listenv<-as.list(c(1:LengTime))
+ names(my.listenv)<- rep("i",LengTime)
+ globalMyEnv <-new.env()
+ globalMyEnv$time <- time
+ globalMyEnv$my.env <- my.env
element <- matrix(0, nrow =info.int at dimIntegrand[1], ncol = 1)
-
- for(i in c(1:(length(time)-1))){
- assign(info.var at upper.var,time[i+1],envir=my.env)
- Inter2 <- lapply(info.int at IntegrandList,
- FUN = my.fun, my.env = my.env, i = i)
- for(j in c(1:info.int at dimIntegrand[1])){
- element[j,] <- sum(diag(matrix(unlist(Inter2[PosInMatr[j,]]),
- ncol = info.int at dimIntegrand[2])%*%M.dX[,c(1:i)]))
- }
- res[,i] <- element
- }
- # LengTime<-(length(time)-1)
- # my.listenv<-as.list(c(1:LengTime))
- # names(my.listenv)<- rep("i",LengTime)
- # globalMyEnv <-new.env()
- # globalMyEnv$time <- time
- # globalMyEnv$my.env <- my.env
- # my.listenv2<-lapply(X=my.listenv,
- # globalMyEnv=globalMyEnv,
- # FUN = function(X,globalMyEnv){
- # assign(info.var at upper.var,globalMyEnv$time[X+1],
- # envir= globalMyEnv$my.env)
- # Inter2 <- lapply(info.int at IntegrandList,
- # FUN = my.fun, my.env = globalMyEnv$my.env,
- # i = X)
- # for(j in c(1:info.int at dimIntegrand[1])){
- # element[j,] <- sum(diag(matrix(unlist(Inter2[PosInMatr[j,]]),
- # ncol = info.int at dimIntegrand[2])%*%M.dX[,c(1:X)]))
- # }
- #
- # list(element)
- # })
- # res<-as.numeric(unlist(my.listenv2))
- # res<-matrix(res,info.int at dimIntegrand[1],(length(time)-1))
+ my.listenv2<-lapply(X=my.listenv,
+ globalMyEnv=globalMyEnv,
+ FUN = function(X,globalMyEnv){
+ assign(info.var at upper.var,globalMyEnv$time[X+1],
+ envir= globalMyEnv$my.env)
+ assign(object at model@time.variable,globalMyEnv$time[c(1:X)],
+ envir= globalMyEnv$my.env)
+ Inter2 <- lapply(info.int at IntegrandList,
+ FUN = my.fun, my.env = globalMyEnv$my.env,
+ i = X)
+ for(j in c(1:info.int at dimIntegrand[1])){
+ element[j,] <- M.dX[,c(1:X)]%*%matrix(unlist(Inter2[PosInMatr[j,]]),
+ ncol = info.int at dimIntegrand[2])
+ }
+ list(element)
+ })
+ res<-as.numeric(unlist(my.listenv2))
+ res<-matrix(res,info.int at dimIntegrand[1],(length(time)-1))
res <- cbind(0,res)
rownames(res) <- object at Integral@variable.Integral at out.var
my.data <- zoo(x = t(res), order.by = time)
Modified: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/R/simulateForPpr.R 2016-11-08 19:25:17 UTC (rev 511)
@@ -125,12 +125,19 @@
auxIntMy <- matrix(auxIntMy, Kern at Integrand@dimIntegrand[1],
Kern at Integrand@dimIntegrand[2], byrow=T)
-
- auxInt <- setIntegral(yuima = Model,
+ if(object at Kernel@variable.Integral at var.dx==object at Kernel@variable.Integral at var.time){
+ auxInt <- setIntegral(yuima = Model,
integrand = auxIntMy,
var.dx = Model at time.variable,
upper.var = dummyUpperTime,
lower.var = Kern at variable.Integral@lower.var)
+ }else{
+ auxInt <- setIntegral(yuima = Model,
+ integrand = auxIntMy,
+ var.dx =object at Kernel@variable.Integral at var.dx ,
+ upper.var = dummyUpperTime,
+ lower.var = Kern at variable.Integral@lower.var)
+ }
randomGenerator<-object at model@measure$df
if(samp at regular){
tForMeas<-samp at delta
@@ -146,6 +153,11 @@
}
Noise.L <- t(rand(object = randomGenerator, n=NumbIncr, param=measureparam))
Noise.W <- t(rnorm(NumbIncr, 0,tForMeas))
+ if(length(object at model@diffusion[[1]])>1){
+ for(i in c(2:length(object at model@diffusion[[1]]))){
+ Noise.W <- rbind(Noise.W, rnorm(NumbIncr, 0,tForMeas))
+ }
+ }
if(missing(xinit)){
simg <- simulate(object = auxg, true.parameter = true.parameter[auxg at Output@param at allparam],
sampling = samp, hurst = hurst,
@@ -169,8 +181,10 @@
globPos <- c(globPos,Pos)
}
}
+ globPos <- unique(globPos)
+ globPos <- globPos[(globPos<=samp at n)]
NewNoise.L <- Noise.L
- cod <- object at Ppr@counting.var%in%Model at solve.variable
+ cod <-Model at solve.variable%in%object at Ppr@counting.var
NeWNoise.W<-Noise.W
NeWNoise.W[cod,] <- 0
NewNoise.L[cod,] <- 0
Modified: pkg/yuima/R/simulateMultiProcess.R
===================================================================
--- pkg/yuima/R/simulateMultiProcess.R 2016-11-06 14:24:33 UTC (rev 510)
+++ pkg/yuima/R/simulateMultiProcess.R 2016-11-08 19:25:17 UTC (rev 511)
@@ -273,7 +273,7 @@
# yuima at data <- euler(xinit, yuima, dW, yuimaEnv)
# }
- if(is(sdeModel,"yuima.multimodel")){
+ if(is(sdeModel,"yuima.multimodel")&&!is(sdeModel at measure$df,"yuima.law")){
if(length(sdeModel at measure.type)==1){
if(sdeModel at measure.type=="CP"){
intens <- as.character(sdeModel at measure$intensity)
@@ -305,12 +305,21 @@
dimJumpCoeff <- length(yuima at model@jump.coeff[[1]])
dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
Dumsolve.variable<-paste0("MyLevyDum",c(1:dimJumpCoeff))
- LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
+ if(!is(sdeModel at measure$df,"yuima.law")){
+ LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
diffusion = NULL,
jump.coeff = dumjumpCoeff,
df = as.character(sdeModel at measure$df$expr),
measure.type = sdeModel at measure.type,
solve.variable = Dumsolve.variable)
+ }else{
+ LevyMod <- setModel(drift=rep("0",dimJumpCoeff),
+ diffusion = NULL,
+ jump.coeff = dumjumpCoeff,
+ measure = sdeModel at measure,
+ measure.type = sdeModel at measure.type,
+ solve.variable = Dumsolve.variable)
+ }
yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
yuimaLevy at model@dimension <- dimJumpCoeff
@@ -481,38 +490,42 @@
#} else{
#stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
#}
- dumStringMeas <- toString(sdeModel at 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 at jump.variable,sdeModel at 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)))
+ if(!is(sdeModel at measure$df,"yuima.law")){
+ dumStringMeas <- toString(sdeModel at 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 at jump.variable,sdeModel at 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 at sampling@n)
+ F.env <- new.env(parent=env)
+ N_sharp <- unique(yuima at sampling@n)
+ TrueDelta <- unique(yuima at 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 at 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)
+ randJ<- env$dL
+ }
RANDLevy <- apply(randJ,2,cumsum)
tsX <- zoo(x=RANDLevy)
+
yuimaData <- setYuima(data=setData(tsX, delta=TrueDelta))
#yuimaData <- subsampling(yuimaData, sampling=yuima at sampling)
return(yuimaData)
More information about the Yuima-commits
mailing list