[Yuima-commits] r676 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 17 01:53:09 CEST 2018
Author: lorenzo
Date: 2018-09-17 01:53:08 +0200 (Mon, 17 Sep 2018)
New Revision: 676
Modified:
pkg/yuima/R/simulateForPpr.R
Log:
Simulate PPR with feedback in the SDE
Modified: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R 2018-07-18 08:19:47 UTC (rev 675)
+++ pkg/yuima/R/simulateForPpr.R 2018-09-16 23:53:08 UTC (rev 676)
@@ -182,16 +182,441 @@
methodfGn = methodfGn, sampling = sampling,
subsampling = subsampling)
}else{
- object <- aux.simulatPPRROldNew(object, nsim = nsim, seed = seed,
+ if(object at PPR@RegressWithCount){
+ object <-aux.simulatPPRWithCount(object, 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, sampling = sampling,
+ subsampling = subsampling)
+
+ }else{
+
+ object <- aux.simulatPPRROldNew(object, 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, sampling = sampling,
subsampling = subsampling)
+ }
}
return(object)
}
+eulerPPR<-function(xinit,yuima,Initial,Terminal, 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
+
+ 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)
+ }
+
+ dX<-vector(mode="numeric",length(dX_dummy))
+
+ for(i in 1:length(dX_dummy)){
+ dX[i] <- eval(dX_dummy[i], env)
+ }
+ }
+
+ delta <- yuima at sampling@delta
+
+ if(!length(sdeModel at measure.type)){
+
+ 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 = Initial) #LM
+ }else{
+ 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)))
+
+ p.b <- function(t, X=numeric(d.size)){
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], X[i], env)
+ }
+ assign(modeltime, t, env)
+ 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){
+ dW <- rbind( rep(1, n)*delta , dW)
+ }
+
+ JP <- sdeModel at jump.coeff
+ mu.size <- length(JP)
+
+ 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)
+
+ j.size <- length(JP[[1]])
+ tmp <- matrix(0, mu.size, j.size)
+
+ for(i in 1:mu.size){
+ for(j in 1:j.size){
+ tmp[i,j] <- eval(JP[[i]][j],env)
+ }
+ }
+ return(tmp)
+ }
+
+ dZ <- dL
+
+ if(is.null(dim(dZ)))
+ dZ <- matrix(dZ,nrow=1)
+ for(i in 1:n){
+ assign(sdeModel at jump.variable, dZ[,i], env)
+
+ if(sdeModel at J.flag){
+ dZ[,i] <- 1
+ }
+ tmp.j <- p.b.j(t=Initial+(i - 1)*delta, X=dX) # YK
+ if(sum(dim(tmp.j))==2)
+ tmp.j <- as.numeric(tmp.j)
+ dX <- dX + p.b(t=Initial+(i - 1)*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i] # YK
+ X_mat[, i+1] <- dX
+ }
+ #tsX <- ts(data=t(X_mat), deltat=delta, start=yuima at sampling@Initial)
+
+ }
+ return(X_mat)
+}
+
+aux.simulatPPRWithCount<-function(object, 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 = 0.5,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling){
+
+ samp <- sampling
+ Model <- object at model
+ gFun <- object at gFun
+ Kern <- object at Kernel
+ object at sampling <- samp
+ randomGenerator<-object at model@measure$df
+
+ if(missing(increment.W) | is.null(increment.W)){
+
+ if( Model at hurst!=0.5 ){
+
+ grid<-sampling2grid(object at sampling)
+ isregular<-object at sampling@regular
+
+ if((!isregular) || (methodfGn=="Cholesky")){
+ dW<-CholeskyfGn(grid, Model at hurst,Model at noise.number)
+ yuima.warn("Cholesky method for simulating fGn has been used.")
+ } else {
+ dW<-WoodChanfGn(grid, Model at hurst,Model at noise.number)
+ }
+
+ } else {
+
+ # delta<-(Terminal-Initial)/n
+ delta <- samp at delta
+ if(!is.Poisson(Model)){ # if pure CP no need to setup dW
+ dW <- rnorm(samp at n * Model at noise.number, 0, sqrt(delta))
+ dW <- matrix(dW, ncol=samp at n, nrow=Model at noise.number,byrow=TRUE)
+ } else {
+ dW <- matrix(0,ncol=samp at n,nrow=1) # maybe to be fixed
+ }
+ }
+
+ } else {
+ dW <- increment.W
+ }
+ if(missing(xinit)){
+ if(length(object at model@xinit)!=0){
+ xinit<-numeric(length=length(object at model@xinit))
+ for(i in c(1:object at model@equation.number))
+ xinit[i] <- eval(object at model@xinit[i])
+ }else{
+ xinit <- rep(0,object at model@equation.number)
+ object at model@xinit<-xinit
+ }
+ }
+ if(missing(hurst)){
+ hurst<-0.5
+ }
+ if(samp at regular){
+ tForMeas<-samp at delta
+ NumbIncr<-samp at n
+ if(missing(true.parameter)){
+ eval(parse(text= paste0("measureparam$",
+ object at model@time.variable," <- tForMeas",collapse="")))
+ }else{
+ measureparam<-true.parameter[object at model@parameter at measure]
+ eval(parse(text= paste0("measureparam$",
+ object at model@time.variable," <- tForMeas",collapse="")))
+
+ }
+ Noise.L <- t(rand(object = randomGenerator, n=NumbIncr, param=measureparam))
+ rownames(Noise.L)<-Model at solve.variable
+ #dIncr <- apply(cbind(0,Noise.L),1,diff)
+ Noise.count <- Noise.L[object at PPR@counting.var,]
+ Noise.Laux <- Noise.L
+ for(i in c(1:length(object at PPR@counting.var))){
+ Noise.Laux[object at PPR@counting.var[i],]<-0
+ }
+ }
+ myenv<-new.env()
+ par.len <- length(object at PPR@allparam)
+ if(par.len>0){
+ for(i in 1:par.len){
+ pars <- object at PPR@allparam[i]
+
+ for(j in 1:length(true.parameter)){
+ if( is.na(match(pars, names(true.parameter)[j]))!=TRUE){
+ assign(object at PPR@allparam[i], true.parameter[[j]],myenv)
+ }
+ }
+ }
+ }
+ assign("dL",Noise.Laux,myenv)
+
+
+
+ CovariateSim<- eulerPPR(xinit=xinit,yuima=object,dW=dW,
+ Initial=samp at Initial,Terminal=samp at Terminal,
+ env=myenv)
+ rownames(CovariateSim)<- Model at solve.variable
+ assign("info.PPR", object at PPR, myenv)
+
+ dimCov <- length(object at PPR@covariates)
+ if (dimCov>0){
+ for(j in c(1:dimCov)){
+ assign(object at PPR@covariates[j],
+ as.numeric(CovariateSim[object at PPR@covariates[j],1]),
+ envir = myenv)
+ }
+ }
+
+ ExprHaz <- constHazIntPr(g.Fun = object at gFun@formula,
+ Kern.Fun = object at Kernel, covariates = object at PPR@covariates,
+ counting.var = object at PPR@counting.var)$Intens
+
+
+ # Start Simulation PPR
+ compErrHazR4 <- function(samp, Kern,
+ capitalTime, Model,
+ my.env, ExprHaz, Time, dN,
+ Index, pos){
+ assign(Kern at variable.Integral@var.time, Time, envir = my.env)
+
+ assign(Model at time.variable, capitalTime, envir = my.env)
+ l <- 1
+ for(i in c(1:length(Kern at variable.Integral@var.dx)) ){
+ if(any(Kern at variable.Integral@var.dx[i]==my.env$info.PPR at counting.var)){
+ assign(paste0("d",Kern at variable.Integral@var.dx[i]), dN[l,], envir =my.env)
+ l <- l + 1
+ }
+ if(Kern at variable.Integral@var.dx[i]%in%my.env$info.PPR at covariates){
+ assign(paste0("d",Kern at variable.Integral@var.dx[i]),
+ diff(c(0,my.env[[Kern at variable.Integral@var.dx[i]]])) ,
+ envir =my.env)
+ }
+ }
+ condPointIngrid <- samp at grid[[1]]<=my.env$t
+ PointIngridInt <- samp at grid[[1]][condPointIngrid]
+ CondJumpGrid <- PointIngridInt %in% my.env$s
+ assign("CondJumpGrid", CondJumpGrid, envir = my.env)
+ Lambda <- NULL
+ # for(h in c(1:Index)){
+ # Lambdadum <- eval(ExprHaz[[h]], envir = my.env)
+ # Lambda <- rbind(Lambda,Lambdadum)
+ #
+ Lambda <- eval(ExprHaz[[pos]], envir = my.env)
+ # rownames(Lambda) <- my.env$info.PPR at counting.var
+ return(Lambda)
+
+ }
+ dN <- matrix(0,object at gFun@dimension[1],object at gFun@dimension[2])
+ grid <- samp at grid[[1]]
+ const <- -log(runif(gFun at dimension[1]))
+ condMyTR <- const<delta
+ while(any(condMyTR)){
+ if(sum(condMyTR)==0){
+ const <- -log(runif(length(condMyTR)))
+ condMyTR <- const<delta
+ }else{
+ const[condMyTR] <- -log(runif(sum(condMyTR)))
+ condMyTR <- const<delta
+ }
+ }
+ jumpT<-NULL
+ i <- 1
+ dimGrid <-length(grid)
+ cond <- const
+ Index <- gFun at dimension[1]
+ inter_i <- rep(i,Index)
+ noExit<-rep(T,Index)
+ while(any(noExit)){
+
+ for(j in c(1:Index)){
+ HazardRate<-0
+ while(cond[j]>0 && noExit[j]){
+ lambda<-compErrHazR4(samp, Kern, capitalTime=samp at grid[[1]][inter_i[j]],
+ Model, myenv, ExprHaz, Time=jumpT, dN, Index, j)
+ # lambda<-hawkesInt(mu=mu, alpha=alpha, beta=beta,
+ # timet=grid[i], JumpT=jumpT)
+ incrlambda <- lambda*delta
+ HazardRate <- HazardRate+incrlambda
+ cond[j] <- const[j]-HazardRate
+ inter_i[j]<-inter_i[j]+1
+ if(inter_i[j]>=(dimGrid-1)){
+ noExit[j] <- FALSE
+ }
+ if(inter_i[j]<dim(CovariateSim)[2]){
+ 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(CovariateSim[object at PPR@covariates[j],1:inter_i[j]]),
+ envir = myenv)
+ }
+ }
+
+
+ # Line 354 necessary for the development of the code.
+ # cat("\n ", i, grid[i])
+ }
+ }
+ }
+
+ i <- min(inter_i)
+
+ if(any(noExit)){
+ if(i<dim(CovariateSim)[2]){
+ jumpT<-c(jumpT,grid[i])
+
+ if(dim(dN)[2]==1 & all(dN[,1]==0)){
+ dN[i==inter_i,1] <- Noise.count[i-1]
+ Noise.Laux[object at PPR@counting.var,i-1]<-Noise.count[i-1]
+ dumdN <- dN
+ }else{
+ dumdN <- rep(0,Index)
+ dumdN[i==inter_i] <- Noise.count[i-1]
+ Noise.Laux[object at PPR@counting.var,i-1] <- dumdN[i==inter_i]
+ dN <- cbind(dN,dumdN)
+ }
+ cat("\n ", i, grid[i])
+ assign("dL",Noise.Laux,myenv)
+
+ CovariateSim<- eulerPPR(xinit=xinit,yuima=object,dW=dW,
+ Initial=samp at Initial,Terminal=samp at Terminal,
+ env=myenv)
+
+ rownames(CovariateSim)<- Model at solve.variable
+
+
+ const <- -log(runif(object at gFun@dimension[1]))
+ condMyTR <- const<delta
+ while(any(condMyTR)){
+ if(sum(condMyTR)==0){
+ const <- -log(runif(length(condMyTR)))
+ condMyTR <- const<delta
+ }else{
+ const[condMyTR] <- -log(runif(sum(condMyTR)))
+ condMyTR <- const<delta
+ }
+ }
+
+ cond <- const
+ if(all(noExit)){
+ inter_i <- rep(i, Index)
+ }else{
+ if(any(noExit)){
+ inter_i[noExit] <- i
+ inter_i[!noExit] <- samp at n+1
+ }
+ }
+
+ }
+ }
+ }
+ tsX <- ts(data=t(CovariateSim), deltat=delta, start=object at sampling@Initial)
+ object at data <- setData(original.data=tsX)
+ for(i in 1:length(object at data@zoo.data))
+ index(object at data@zoo.data[[i]]) <- object at sampling@grid[[1]] ## to be fixed
+
+ #object at model@hurst <-tmphurst
+
+ if(missing(subsampling))
+ return(object)
+ subsampling(object, subsampling)
+
+}
+
aux.simulatPPRROldNew<-function(object, nsim = nsim, seed = seed,
xinit = xinit, true.parameter = true.parameter,
space.discretized = space.discretized, increment.W = increment.W,
More information about the Yuima-commits
mailing list