[Yuima-commits] r654 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 30 16:03:41 CEST 2018
Author: lorenzo
Date: 2018-05-30 16:03:41 +0200 (Wed, 30 May 2018)
New Revision: 654
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/simulateForPpr.R
Log:
Simulation PPR covariates kernel
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2018-05-23 04:55:27 UTC (rev 653)
+++ pkg/yuima/DESCRIPTION 2018-05-30 14:03:41 UTC (rev 654)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.8.1
+Version: 1.8.2
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2)
Author: YUIMA Project Team
Modified: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R 2018-05-23 04:55:27 UTC (rev 653)
+++ pkg/yuima/R/simulateForPpr.R 2018-05-30 14:03:41 UTC (rev 654)
@@ -52,13 +52,79 @@
}
)
-constHazIntPr <- function(g.Fun , Kern.Fun){
+constHazIntPr <- function(g.Fun , Kern.Fun, covariates, counting.var){
numb.Int <- length(g.Fun)
Int.Intens <- list()
for(i in c(1:numb.Int)){
dum.g <- as.character(g.Fun[i])
+ dum.g <- gsub("(", "( ", fixed=TRUE,x = dum.g)
+ dum.g <- gsub(")", " )", fixed=TRUE,x = dum.g)
+
+ if(length(counting.var)>0){
+ for(i in c(1:length(counting.var))){
+ my.countOld <- paste0(counting.var[i] ," ")
+ #my.countNew <- paste0("as.numeric(", counting.var[i] ,")")
+ my.countNew <- paste0("(", counting.var[i] ,")")
+ dum.g <- gsub(my.countOld, my.countNew, x = dum.g, fixed=TRUE)
+ my.countOld <- paste0(counting.var[i] ,"[",Kern.Fun at variable.Integral@upper.var,"]")
+
+ my.countNew <- paste0("(", counting.var[i] ,")")
+ dum.g <- gsub(my.countOld, my.countNew, x = dum.g, fixed=TRUE)
+ }
+ }
+ if(length(covariates)>0){
+ for(i in c(1:length(covariates))){
+ my.covarOld <- paste0(covariates[i] ," ")
+ my.covarNew <- covariates[i]
+ dum.g <- gsub(my.covarOld, my.covarNew, x = dum.g, fixed=TRUE)
+ my.covarOld <- paste0(covariates[i] ,"[",Kern.Fun at variable.Integral@upper.var,"]")
+ my.covarNew <- covariates[i]
+ dum.g <- gsub(my.covarOld, my.covarNew, x = dum.g, fixed=TRUE)
+ }
+ }
+
dum.g <- paste("tail(",dum.g,", n=1L)")
dum.Ker <- as.character(Kern.Fun at Integrand@IntegrandList[[i]])
+ dum.Ker <- gsub("(", "( ", fixed=TRUE,x = dum.Ker)
+ dum.Ker <- gsub(")", " )", fixed=TRUE,x = dum.Ker)
+
+ if(length(counting.var)>0){
+ for(i in c(1:length(counting.var))){
+ my.countOld <- paste0(counting.var[i] ," ")
+ my.countNew <- paste0( counting.var[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@var.time ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+ my.countOld <- paste0(counting.var[i] ,"[",Kern.Fun at variable.Integral@upper.var,"]")
+ my.countNew <- paste0( counting.var[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@upper.var ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ my.countOld <- paste0(counting.var[i] ,"[",Kern.Fun at variable.Integral@var.time,"]")
+ my.countNew <- paste0(counting.var[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@var.time ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ }
+ }
+ if(length(covariates)>0){
+ for(i in c(1:length(covariates))){
+
+ my.countOld <- paste0(covariates[i] ," ")
+ my.countNew <- paste0( covariates[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@var.time ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ my.countOld <- paste0(covariates[i] ,"[",Kern.Fun at variable.Integral@upper.var,"]")
+ my.countNew <- paste0( covariates[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@upper.var ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ my.countOld <- paste0(covariates[i] ,"[",Kern.Fun at variable.Integral@var.time,"]")
+ my.countNew <- paste0( covariates[i] ,
+ "[ as.character( ",Kern.Fun at variable.Integral@var.time ," ) ]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+ }
+ }
dif.dx <- paste("d",Kern.Fun at variable.Integral@var.dx, sep="")
dum.Ker <- paste(dum.Ker,dif.dx, sep = "*")
if(length(dum.Ker)>1){
@@ -73,6 +139,9 @@
res <- list(Intens = Int.Intens)
}
+
+
+
aux.simulatPPR<- function(object, nsim = nsim, seed = seed,
xinit = xinit, true.parameter = true.parameter,
space.discretized = space.discretized, increment.W = increment.W,
@@ -104,6 +173,110 @@
increment.L = increment.L, method = method, hurst = 0.5,
methodfGn = methodfGn, sampling = sampling,
subsampling = subsampling){
+ myhawkesP <- function(simMod, Kern,
+ samp, Model, my.env, ExprHaz,
+ Time, dN){
+ noExit<-TRUE
+ const <- -log(runif(1))
+ delta <- samp at delta
+ grid <- samp at grid[[1]]
+ while(const<delta){
+ const <- -log(runif(1))
+ }
+ jumpT<-NULL
+ i <- 1
+ dimGrid <-length(grid)
+ cond <- const
+ allconst <- NULL
+ allcond <- NULL
+ allhaz <- NULL
+ while(noExit){
+ HazardRate<-0
+ while(cond>0 && noExit){
+ #lastJump <- tail(jumpT,n=1L)
+ lambda<-compErrHazR2(simMod, Kern, capitalTime=samp at grid[[1]][i], Model, my.env, ExprHaz,
+ Time=jumpT, dN)
+ # lambda<-hawkesInt(mu=mu, alpha=alpha, beta=beta,
+ # timet=grid[i], JumpT=jumpT)
+ incrlambda <- lambda*delta
+ HazardRate <- HazardRate+incrlambda
+ cond <- const-HazardRate
+ i<-i+1
+ if(i>=(dimGrid-1)){
+ noExit <- FALSE
+ }
+ if(i<dim(simMod at data@original.data)[1]){
+ dimCov <- length(object at PPR@covariates)
+
+ if (dimCov>0){
+ for(j in c(1:dimCov)){
+ my.covdata <- simMod at data@original.data[1:i,object at PPR@covariates[j]]
+ names(my.covdata) <-simMod at sampling@grid[[1]][1:i]
+
+ assign(object at PPR@covariates[j],
+ my.covdata,
+ envir = my.env)
+
+ # assign(object at PPR@covariates[j],
+ # as.numeric(simMod at data@original.data[1:i,object at PPR@covariates[j]]),
+ # envir = my.env)
+ }
+ }
+
+
+ # Line 354 necessary for the development of the code.
+ # cat("\n ", i, grid[i])
+ }
+ }
+ if(i<dim(simMod at data@original.data)[1]){
+ jumpT<-c(jumpT,grid[i])
+ # if(i==7001){
+ # cat("\n",noExit)
+ # }
+ if(dN[1]==0){
+ dN <- 1
+ }else{
+ dN <- c(dN,1)
+ }
+ names(dN)<-jumpT
+ allhaz <- c(allhaz,HazardRate)
+ allcond <- c(allcond,cond)
+ cond <- const
+ allconst <- c(allconst, const)
+ const <- -log(runif(1))
+ while(const<delta){
+ const <- -log(runif(1))
+ }
+ }
+ }
+ return(list(jumpT=jumpT,allcond=allcond,allconst=allconst, allhaz=allhaz))
+ }
+
+
+ compErrHazR2 <- function(simMod, Kern,
+ capitalTime, Model, my.env, ExprHaz,
+ Time, dN){
+ # dummyLambda <- numeric(length=(samp at n+1))
+ if(length(Kern at variable.Integral@var.dx)==1){
+ # MyPos <- sum(samp at grid[[1]]<=tail(Time,n=1L))
+ assign(Kern at variable.Integral@var.time, Time, envir = my.env)
+ # cond <- -log(cost)-sum(dummyLambda)*samp at delta
+
+ assign(Model at time.variable, capitalTime, envir = my.env)
+ assign(paste0("d",Kern at variable.Integral@var.dx), dN, envir =my.env)
+
+ # condPointIngrid <- simMod at sampling@grid[[1]]<=my.env$t
+ # PointIngridInt <- simMod at sampling@grid[[1]][condPointIngrid]
+ # CondJumpGrid <- PointIngridInt %in% my.env$s
+ # assign("CondJumpGrid", CondJumpGrid, envir = my.env)
+
+ Lambda <- eval(ExprHaz[[1]], envir=my.env)
+ return(Lambda)
+ }else{
+ return(NULL)
+ }
+ }
+
# We need an expression for the evaluation of the hazard
if(missing(hurst)){
@@ -262,7 +435,8 @@
Data.tot <- as.matrix(simMod at data@original.data)
ExprHaz <- constHazIntPr(g.Fun = object at gFun@formula,
- Kern.Fun = object at Kernel)$Intens
+ Kern.Fun = object at Kernel, covariates = object at PPR@covariates,
+ counting.var = object at PPR@counting.var)$Intens
if(length(ExprHaz)==1){
@@ -278,104 +452,30 @@
dimCov <- length(object at PPR@covariates)
if (dimCov>0){
for(j in c(1:dimCov)){
- assign(object at PPR@covariates[j],
- as.numeric(simMod at data@original.data[1,object at PPR@covariates[j]]),
+ # assign(object at PPR@covariates[j],
+ # as.numeric(simMod at data@original.data[1,object at PPR@covariates[j]]),
+ # envir = my.env)
+ my.covdata <- simMod at data@original.data[1,object at PPR@covariates[j]]
+ names(my.covdata) <-simMod at sampling@grid[[1]][1]
+
+ assign(object at PPR@covariates[j],
+ my.covdata,
envir = my.env)
+
}
}
- compErrHazR2 <- function(simMod, Kern,
- capitalTime, Model, my.env, ExprHaz,
- Time, dN){
- # dummyLambda <- numeric(length=(samp at n+1))
- if(length(Kern at variable.Integral@var.dx)==1){
- # MyPos <- sum(samp at grid[[1]]<=tail(Time,n=1L))
- assign(Kern at variable.Integral@var.time, Time, envir = my.env)
- # cond <- -log(cost)-sum(dummyLambda)*samp at delta
+
- assign(Model at time.variable, capitalTime, envir = my.env)
- assign(paste0("d",Kern at variable.Integral@var.dx), dN, envir =my.env)
- Lambda <- eval(ExprHaz[[1]], envir=my.env)
- return(Lambda)
- }else{
- return(NULL)
- }
- }
-
IntensityProc <- 0
#set.seed(1)
dN<-0
- myhawkesP <- function(simMod, Kern,
- samp, Model, my.env, ExprHaz,
- Time){
- noExit<-TRUE
- const <- -log(runif(1))
- delta <- samp at delta
- grid <- samp at grid[[1]]
- while(const<delta){
- const <- -log(runif(1))
- }
- jumpT<-0
- i <- 1
- dimGrid <-length(grid)
- cond <- const
- allconst <- NULL
- allcond <- NULL
- allhaz <- NULL
- while(noExit){
- HazardRate<-0
- while(cond>0 && noExit){
- #lastJump <- tail(jumpT,n=1L)
- lambda<-compErrHazR2(simMod, Kern, capitalTime=samp at grid[[1]][i], Model, my.env, ExprHaz,
- Time=jumpT, dN)
- # lambda<-hawkesInt(mu=mu, alpha=alpha, beta=beta,
- # timet=grid[i], JumpT=jumpT)
- incrlambda <- lambda*delta
- HazardRate <- HazardRate+incrlambda
- cond <- const-HazardRate
- i<-i+1
- if(i>=(dimGrid-1)){
- noExit <- FALSE
- }
- if(i<dim(simMod at data@original.data)[1]){
- dimCov <- length(object at PPR@covariates)
-
- if (dimCov>0){
- for(j in c(1:dimCov)){
- assign(object at PPR@covariates[j],
- as.numeric(simMod at data@original.data[1:i,object at PPR@covariates[j]]),
- envir = my.env)
- }
- }
-
+
- # Line 354 necessary for the development of the code.
- # cat("\n ", i, grid[i])
- }
- }
- if(i<dim(simMod at data@original.data)[1]){
- jumpT<-c(jumpT,grid[i])
- # if(i==7001){
- # cat("\n",noExit)
- # }
- dN<-c(dN,1)
- allhaz <- c(allhaz,HazardRate)
- allcond <- c(allcond,cond)
- cond <- const
- allconst <- c(allconst, const)
- const <- -log(runif(1))
- while(const<delta){
- const <- -log(runif(1))
- }
- }
- }
- return(list(jumpT=jumpT,allcond=allcond,allconst=allconst, allhaz=allhaz))
- }
-
prova1<-myhawkesP(simMod, Kern,
samp, Model, my.env, ExprHaz,
- Time)
+ Time, dN)
Time<-prova1$jumpT
# return(Time)
@@ -507,7 +607,8 @@
# # cat(sprintf("\n%.5f ", tail(Time,n=1L)))
# # }
# }
- cond <- samp at grid[[1]][-1] %in% Time[-1]
+ #cond <- samp at grid[[1]][-1] %in% Time[-1]
+ cond <- samp at grid[[1]][-1] %in% Time
countVar <- Model at solve.variable %in% object at PPR@counting.var
increment.L[!cond, countVar]<-0
if(missing(xinit)){
More information about the Yuima-commits
mailing list