[Yuima-commits] r457 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jul 9 22:57:55 CEST 2016
Author: lorenzo
Date: 2016-07-09 22:57:55 +0200 (Sat, 09 Jul 2016)
New Revision: 457
Added:
pkg/yuima/R/PointProcessClasses.R
pkg/yuima/R/setPpr.R
pkg/yuima/R/simulateForPpr.R
pkg/yuima/man/info.Ppr.Rd
pkg/yuima/man/setHawkes.Rd
pkg/yuima/man/setPpr.Rd
pkg/yuima/man/yuima.Hawkes.Rd
pkg/yuima/man/yuima.Ppr.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/setHawkes.R
pkg/yuima/R/simulateForMapsIntegralAndOperator.R
pkg/yuima/R/simulateMultiProcess.R
Log:
New Methods and Classes for Point Process are now available
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-07-08 04:44:22 UTC (rev 456)
+++ pkg/yuima/DESCRIPTION 2016-07-09 20:57:55 UTC (rev 457)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.1.2
+Version: 1.1.3
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2016-07-08 04:44:22 UTC (rev 456)
+++ pkg/yuima/NAMESPACE 2016-07-09 20:57:55 UTC (rev 457)
@@ -83,7 +83,10 @@
"variable.Integral",
"Integrand",
"Integral.sde",
-"yuima.Integral"
+"yuima.Integral",
+"info.Ppr",
+"yuima.Ppr",
+"yuima.Hawkes"
)
exportMethods(
@@ -126,6 +129,8 @@
export(setCharacteristic)
export(setCarma)
export(setPoisson)
+export(setPpr)
+export(setHawkes)
export(dconst)
export(rconst)
Added: pkg/yuima/R/PointProcessClasses.R
===================================================================
--- pkg/yuima/R/PointProcessClasses.R (rev 0)
+++ pkg/yuima/R/PointProcessClasses.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -0,0 +1,73 @@
+setClass("info.Ppr",
+ representation(allparam = "character",
+ allparamPpr = "character",
+ common ="character",
+ counting.var = "character",
+ var.dx = "character",
+ upper.var = "character",
+ lower.var = "character",
+ covariates = "character",
+ var.dt = "character",
+ additional.info = "character",
+ Info.measure = "list")
+ )
+
+setClass("yuima.Ppr",
+ representation(Ppr = "info.Ppr",
+ gFun = "info.Output",
+ Kernel = "Integral.sde"),
+ contains="yuima"
+)
+
+setMethod("initialize",
+ "info.Ppr",
+ function(.Object,
+ allparam = character(),
+ allparamPpr = character(),
+ common = character(),
+ counting.var = character(),
+ var.dx = character(),
+ upper.var = character(),
+ lower.var = character(),
+ covariates = character(),
+ var.dt = character(),
+ additional.info = character(),
+ Info.measure = list()){
+ .Object at allparam <- allparam
+ .Object at allparamPpr <- allparamPpr
+ .Object at common <- common
+ .Object at counting.var <- counting.var
+ .Object at var.dx <- var.dx
+ .Object at upper.var <- upper.var
+ .Object at lower.var <- lower.var
+ .Object at covariates <- covariates
+ .Object at var.dt <- var.dt
+ .Object at additional.info <- additional.info
+ .Object at Info.measure <- Info.measure
+ return(.Object)
+ }
+)
+
+setMethod("initialize",
+ "yuima.Ppr",
+ function(.Object,
+ Ppr = new("info.Ppr"),
+ gFun = new("info.Output"),
+ Kernel = new("Integral.sde"),
+ yuima = new("yuima")){
+ #.Object at param <- param
+ .Object at Ppr <- Ppr
+ .Object at gFun <- gFun
+ .Object at Kernel <- Kernel
+ .Object at data <- yuima at data
+ .Object at model <- yuima at model
+ .Object at sampling <- yuima at sampling
+ .Object at characteristic <- yuima at characteristic
+ .Object at functional <- yuima at functional
+ return(.Object)
+ }
+)
+
+setClass("yuima.Hawkes",
+ contains="yuima.Ppr"
+)
Modified: pkg/yuima/R/setHawkes.R
===================================================================
--- pkg/yuima/R/setHawkes.R 2016-07-08 04:44:22 UTC (rev 456)
+++ pkg/yuima/R/setHawkes.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -1,49 +1,42 @@
-# setHawkes <- function(yuima, counting.var="N", gFun, Kernel,
-# var.dx, var.dt = "s", lambda.var = "lambda",
-# lower.var="0", upper.var = "t",
-# nrow =1 ,ncol=1){
-#
-# g <- yuima:::setMaps(func = gFun, yuima = yuima,
-# nrow = nrow, ncol = ncol)
-#
-# yuimadum <- yuima
-# yuimadum at time.variable <- var.dt
-#
-# HawkesType <- FALSE
-# if(counting.var %in% var.dx){
-# HawkesType <- TRUE
-# }
-# if(!HawkesType){
-# Integral <- yuima:::setIntegral(yuima=yuimadum,
-# integrand = Kernel, var.dx = var.dx,
-# lower.var = lower.var, upper.var = upper.var,
-# out.var = "", nrow = nrow, ncol = ncol)
-# }else{
-# Integral <- yuima:::setIntegral(yuima=yuimadum,
-# integrand = Kernel, var.dx = var.dx,
-# lower.var = lower.var, upper.var = upper.var,
-# out.var = "", nrow = nrow, ncol = ncol, type ="")
-# }
-# if(g at Output@dimension[1]!=Integral$dimIntegrand[1]){
-# yuima.stop("dimension gFun and kernel mismatch")
-# }
-#
-#
-# allparam <-unique(c(yuima at parameter@all, g at param@allparamMap,
-# Integral$param$IntegrandParam))
-# common <- unique(c(g at param@common, Integral$param$common))
-# paramHawkes <- list(allparam = allparam, common = common,
-# gFun = g at param@allparamMap,
-# Kern = Integral$param$IntegrandParam)
-#
-# # IntPpr<- yuima:::setIntegral(yuima=yuimadum,
-# # integrand = Kernel, var.dx = "N",
-# # lower.var = lower.var, upper.var = upper.var,
-# # out.var = "", nrow = nrow, ncol = ncol)
-#
-# return(list(Count.Proc = counting.var,
-# gFun = list(param=g at param, output=g at Output),
-# Kernel = Integral, paramHawkes = paramHawkes,
-# model = yuima, SelfEx = HawkesType))
-#
-# }
+setHawkes <- function(lower.var="0", upper.var = "t", var.dt = "s",
+ process = "N", dimension = 1, intensity = "lambda",
+ ExpKernParm1="c", ExpKernParm2 ="a",
+ const = "nu", measure = NULL, measure.type = NULL){
+
+ PROCESS <- paste0(process,c(1:dimension))
+ leng <- length(PROCESS)
+
+ mod1 <- setModel(drift = rep("0",leng),
+ diffusion = matrix("0",leng,leng),
+ jump.coeff = diag("1",leng,leng),
+ measure = measure, measure.type = measure.type,
+ solve.variable = PROCESS)
+
+ INTENSITY <- as.list(paste0(intensity,c(1:dimension)))
+
+ gFun <- paste0(const,c(1:dimension))
+
+ Ccoeff<-as.character(MatrCoeff(ExpKernParm1, dimension))
+ Acoeff<-as.character(MatrCoeff(ExpKernParm2, dimension))
+ Kernelpar<-c(Acoeff,Ccoeff)
+
+ Kernel<- matrix(paste0(Ccoeff,"*exp(-",Acoeff,"*(","t-",var.dt,"))"),dimension, dimension)
+
+ res <- aux.setPpr(yuima = mod1, counting.var=PROCESS,
+ gFun, Kernel,
+ var.dx = PROCESS, var.dt = var.dt, lambda.var = INTENSITY,
+ lower.var=lower.var, upper.var = upper.var,
+ nrow =dimension ,ncol=dimension, general = FALSE)
+
+ return(res)
+}
+
+MatrCoeff<-function(lett, dimension){
+ c1<-paste0(lett,c(1:dimension))
+
+ cMatrix<-matrix(NA,dimension,dimension)
+ for(i in c(1:dimension)){
+ cMatrix[i,]<-paste0(c1[i],c(1:dimension))
+ }
+ return(cMatrix)
+}
Added: pkg/yuima/R/setPpr.R
===================================================================
--- pkg/yuima/R/setPpr.R (rev 0)
+++ pkg/yuima/R/setPpr.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -0,0 +1,115 @@
+setPpr <- function(yuima, counting.var="N", gFun, Kernel,
+ var.dx, var.dt = "s", lambda.var = "lambda",
+ lower.var="0", upper.var = "t",
+ nrow =1 ,ncol=1){
+
+ general <- TRUE
+ ret <- aux.setPpr(yuima, counting.var, gFun,
+ Kernel, var.dx, var.dt, lambda.var,
+ lower.var="0", upper.var = "t",
+ nrow =1 ,ncol=1,general =general)
+ return(ret)
+}
+
+
+aux.setPpr <-function(yuima, counting.var="N", gFun, Kernel,
+ var.dx, var.dt = "s", lambda.var = "lambda",
+ lower.var="0", upper.var = "t",
+ nrow =1 ,ncol=1, general){
+ g <- setMap(func = gFun, yuima = yuima,
+ nrow = nrow, ncol = ncol)
+
+ yuimadum <- yuima
+ yuimadum at time.variable <- var.dt
+
+ HawkesType <- FALSE
+ if(counting.var %in% var.dx){
+ HawkesType <- TRUE
+ }
+ if(!HawkesType){
+ Integral <- setIntegral(yuima=yuimadum,
+ integrand = Kernel, var.dx = var.dx,
+ lower.var = lower.var, upper.var = upper.var,
+ out.var = "", nrow = nrow, ncol = ncol)
+ }else{
+ Integral <- setIntegral(yuima=yuimadum,
+ integrand = Kernel, var.dx = var.dx,
+ lower.var = lower.var, upper.var = upper.var,
+ out.var = "", nrow = nrow, ncol = ncol)
+ }
+ if(g at Output@dimension[1]!=Integral at Integral@Integrand at dimIntegrand[1]){
+ yuima.stop("dimension gFun and kernel mismatch")
+ }
+
+
+ allparam <- unique(c(yuima at parameter@all, g at Output@param at allparamMap,
+ Integral at Integral@param.Integral at Integrandparam))
+ common <- unique(c(g at Output@param at common,
+ Integral at Integral@param.Integral at common))
+ paramHawkes <- list(allparam = allparam, common = common,
+ gFun = g at Output@param at allparamMap,
+ Kern = Integral at Integral@param.Integral at Integrandparam)
+
+ # IntPpr<- yuima:::setIntegral(yuima=yuimadum,
+ # integrand y= Kernel, var.dx = "N",
+ # lower.var = lower.var, upper.var = upper.var,
+ # out.var = "", nrow = nrow, ncol = ncol)
+
+ # return(list(Count.Proc = counting.var,
+ # gFun = list(param=g at Output@param, output=g at Output),
+ # Kernel = Integral, paramHawkes = paramHawkes,
+ # model = yuima, SelfEx = HawkesType))
+ yuima1 <- setYuima(model =yuima)
+ type <- yuima at measure.type
+ if(type == "code"){
+ measure <- list(df = as.character(yuima at measure$df$expr))
+ }
+ if(type == "CP"){
+ measure <- list(intensity = as.character(yuima at measure$intensity),
+ df= as.character(yuima at measure$df$expr))
+ }
+ if(general){
+ covariates<-character()
+ if(sum(!(counting.var==yuima at solve.variable))!=0){
+ covariates<-yuima at solve.variable[!(counting.var==yuima at solve.variable)]
+ }
+ Ppr <- new("info.Ppr",
+ allparam = paramHawkes$allparam,
+ allparamPpr = unique(c(paramHawkes$gFun,paramHawkes$Kern)),
+ common = paramHawkes$common,
+ counting.var = counting.var,
+ var.dx = var.dx,
+ upper.var = upper.var,
+ lower.var = lower.var,
+ covariates = covariates,
+ var.dt = var.dt,
+ additional.info = lambda.var,
+ Info.measure = list(type=type,measure=measure))
+
+
+ ret <- new("yuima.Ppr", Ppr = Ppr,
+ gFun = g at Output,
+ Kernel = Integral at Integral,
+ yuima = yuima1)
+ }else{
+
+ Ppr <- new("info.Ppr",
+ allparam = paramHawkes$allparam,
+ allparamPpr = unique(c(paramHawkes$gFun,paramHawkes$Kern)),
+ common = paramHawkes$common,
+ counting.var = counting.var,
+ var.dx = var.dx,
+ upper.var = upper.var,
+ lower.var = lower.var,
+ covariates = character(),
+ var.dt = var.dt,
+ additional.info = "Exponential Hawkes",
+ Info.measure = list(type=type,measure=measure))
+ ret <- new("yuima.Hawkes", Ppr = Ppr,
+ gFun = g at Output,
+ Kernel = Integral at Integral,
+ yuima = yuima1)
+
+ }
+ return(ret)
+}
Modified: pkg/yuima/R/simulateForMapsIntegralAndOperator.R
===================================================================
--- pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-07-08 04:44:22 UTC (rev 456)
+++ pkg/yuima/R/simulateForMapsIntegralAndOperator.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -46,10 +46,12 @@
assign(nam[i],par[i])
}
my.data<-eval(object at Output@formula[[1]])
- for(i in c(2:prod(object at Output@dimension))){
- my.data<-cbind(my.data,
- eval(object at Output@formula[[i]]))
- }
+ if(prod(object at Output@dimension)>1){
+ for(i in c(2:prod(object at Output@dimension))){
+ my.data<-cbind(my.data,
+ eval(object at Output@formula[[i]]))
+ }
+ }
names(my.data)<-object at Output@param at out.var
data1 <- setData(my.data)
Added: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R (rev 0)
+++ pkg/yuima/R/simulateForPpr.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -0,0 +1,487 @@
+setMethod("simulate", "yuima.Hawkes",
+ 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"
+ ...){
+ res <- aux.simulatHawkes(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(res)
+ }
+)
+
+aux.simulatHawkes<- function(object, nsim, seed,
+ xinit, true.parameter, space.discretized, increment.W,
+ increment.L, method, hurst, methodfGn, sampling, subsampling){
+
+# object at Kernel@param.Integral at allparam
+# simOzaki.aux(gFun=object at gFun@formula,a,cCoeff, Time, numJump)
+ ret <-NULL
+}
+
+setMethod("simulate", "yuima.Ppr",
+ 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"
+ ...){
+ res <- aux.simulatPpr(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(res)
+ }
+)
+
+aux.simulatPpr<- 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 = hurst,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling){
+ Time <- sampling at Terminal
+ numbVardx <- length(object at Ppr@var.dx)
+ numbCountVar <- length(object at Ppr@counting.var)
+ U <- runif(numbCountVar)
+
+ my.env<- new.env()
+
+ true.parameter <- unlist(true.parameter)
+
+ if(!all(names(true.parameter)==object at Ppr@allparam)){
+ yuima.stop("true.parameters mismatch the model parameters")
+ }
+ for(i in c(1:length(object at Ppr@allparam))){
+ assign(object at Ppr@allparam[i],true.parameter[object at Ppr@allparam[i]], envir = my.env)
+ }
+
+ assign("t",object at gFun@param at time.var, envir = my.env)
+
+
+ nameu <- object at gFun@param at time.var
+ assign("dt",sampling at delta, envir = my.env)
+
+ if(is.null(increment.W)){
+ dimW <- length(object at model@diffusion[[1]])
+ W <- matrix(rnorm(dimW*sampling at n,mean=0,sd= sqrt(sampling at delta)),nrow=dimW,ncol=sampling at n)
+ }
+ Condcovariate <- TRUE
+ if(is.null(increment.L)){
+ dimL <- length(object at model@jump.coeff[[1]])
+ L <- matrix(0,nrow=dimL,ncol=sampling at n)
+ Condcovariate <- FALSE
+ # if(length(object at Ppr@covariates)!=0)
+ # Condcovariate <- TRUE
+ cond <- !(object at model@solve.variable %in% object at Ppr@counting.var)
+ if(any(cond)){
+ Condcovariate <- TRUE
+ }
+ dimMd <- length(object at model@solve.variable)
+ dumMod <- setModel(drift = rep("0",dimMd),
+ diffusion = matrix("0",dimMd,1),
+ jump.coeff = diag("1",dimMd,dimMd),
+ measure = object at Ppr@Info.measure$measure,
+ measure.type = object at Ppr@Info.measure$type,
+ solve.variable = object at model@solve.variable)
+ if(length(object at model@parameter at measure)!=0){
+ simMod <- simulate(object = dumMod,
+ true.parameter = true.parameter[object at model@parameter at measure],
+ sampling = sampling)
+ }else{
+ simMod <- simulate(object = dumMod,
+ sampling = sampling)
+ }
+
+ L <- t(diff(simMod at data@original.data))
+ }
+
+ assign("Condcovariate",Condcovariate, envir = my.env)
+ assign("W", W, envir = my.env)
+
+
+ rownames(L)<- object at model@solve.variable
+
+ assign("L", L, envir = my.env)
+
+ assign("All.labKern",object at Kernel@variable.Integral,envir = my.env)
+ assign("All.labgFun",object at gFun@param,envir = my.env)
+
+
+
+ Fun1 <- function(u,env){
+ part <- seq(0,u,by=env$dt)
+ env$t<-part[-length(part)]
+ if(Condcovariate){
+ yuima<- object at model
+ for(i in c(1:length(object at Ppr@covariates))){
+ assign(object at Ppr@covariates[i],
+ eval(yuima at xinit[yuima at solve.variable==object at Ppr@covariates[i]],
+ envir = env), envir = env)
+ }
+ if(u!=0){
+ # Mat<-matrix(0,length(yuima at solve.variable),length(env$t)+1)
+ # for(i in c(1:length(yuima at solve.variable))){
+ # Mat[i,1] = eval(yuima at xinit[i],envir = env)
+ # }
+ Linc <- env$L[,c(1:(length(part)-1))]
+ # Linc[yuima at solve.variable!=object at Ppr@covariates,]<-matrix(0,
+ # sum(yuima at solve.variable!=object at Ppr@covariates), dim(Linc)[2])
+ Linc[yuima at solve.variable!=object at Ppr@covariates,] <- 0
+ DumUnderlMod <- simulate(yuima, true.parameter = true.parameter,
+ increment.L = env$L[,c(1:(length(part)-1))],
+ sampling = setSampling(Terminal = u, n= (length(part)-1)))
+
+
+ for(i in c(1:length(object at Ppr@covariates))){
+ VariableDum <- DumUnderlMod at data@original.data[,yuima at solve.variable==object at Ppr@covariates[i]]
+ assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
+ }
+ }
+ }
+ (log(env$U)+sum(eval(env$gFun,envir = env)*env$dt))^2
+ }
+
+ Fun2 <- function(u,env){
+ u <- max(env$old_u,u)
+ dumpart <- seq(0,env$old_u, by=env$dt)
+ part <- seq(env$old_u,u,by=env$dt)
+ t_k <- env$t
+ env$t<-part[-length(part)]
+ if(u>=sampling at Terminal){
+ # Think a better solution
+ my.env$utrue<-u
+ return(0)
+ }
+ if(Condcovariate){
+ LevIncr <- env$L[, length(dumpart)+c(1:(length(env$t)))]
+ LevIncr[object at Ppr@counting.var,]<-0
+ yuima<- object at model
+ xinit<- numeric(length(object at Ppr@covariates))
+ names(xinit)<- object at Ppr@covariates
+ for(i in c(1:length(object at Ppr@covariates))){
+ xinit[i] <- env[[object at Ppr@covariates[i]]]
+ }
+
+ xinitCount <- numeric(length(object at Ppr@counting.var))
+ names(xinitCount) <- object at Ppr@counting.var
+ for(i in c(1:length(xinitCount))){
+ xinitCount[i] <- tail(env[[object at Ppr@counting.var[i]]],n = 1)
+ }
+ xinit <- c(xinit,xinitCount)
+ if(part[length(part)]-part[1]!=0){
+ DumVarCov <- simulate(yuima,
+ true.parameter = true.parameter,
+ increment.L = LevIncr,
+ sampling = setSampling(Terminal = (part[length(part)]-part[1]),
+ n = dim(LevIncr)[2]),
+ xinit=xinit[yuima at solve.variable])
+ for(i in c(1:length(object at Ppr@covariates))){
+ VariableDum <- DumVarCov at data@original.data[,yuima at solve.variable==object at Ppr@covariates[i]]
+ assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
+ }
+ }else{
+ for(i in c(1:length(object at Ppr@covariates))){
+ VariableDum <- xinit[yuima at solve.variable==object at Ppr@covariates[i]]
+ assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
+ }
+ }
+ #Insert Here simulation Covariate
+ }
+ integG <-sum(eval(env$gFun,envir = env)*env$dt)
+ env$s <- unique(c(env$s,t_k))[-length(env$s)]
+ dumt <- env$t
+ num <- length(env$Kern)
+ integKer <- 0
+ for(j in c(1:length(dumt))){
+ env$t <- dumt[j]
+ dumKernInt <- 0
+ for(i in c(1:num)){
+ lab.dx <- env$All.labKern at var.dx[i]
+ dumKernInt <- dumKernInt+sum(eval(env$Kern,envir=env)*diff(eval(env[[lab.dx]])))
+ }
+ integKer <- integKer + dumKernInt
+ }
+ NewTerm <- 0
+ if(env$Condcovariate){
+ ## Insert Her
+ }
+ my.env$utrue<-u
+ (log(env$U)+ integG + integKer+NewTerm)^2
+ }
+
+
+ u <- numeric(length = numbCountVar)
+ names(u) <- object at Ppr@counting.var
+ for(i in c(1:numbCountVar)){
+ assign("gFun", object at gFun@formula[[i]], envir=my.env)
+ assign("U",runif(1),envir = my.env)
+ u[i]<- as.numeric(optim(0,Fun1,env=my.env)$par)
+ }
+
+ t_1 <- min(u)
+
+
+
+ if(t_1>Time){
+ yuima.stop("No jump occurs in the considered time interval.
+ Increasing Terminal in setSampling is suggested")
+ }
+
+ condt1<- u%in%t_1
+
+ namesContVarJump <- names(u[condt1])
+
+ JUMP <- matrix(0,nrow=numbCountVar,ncol=sampling at n)
+
+ rownames(JUMP)<- object at Ppr@counting.var
+ pos<-sum(sampling at grid[[1]][-1]<=t_1)
+ t_1 <- sampling at grid[[1]][-1][pos]
+ recordTime<-c(0,t_1)
+ pos0<-0
+
+ JUMP[namesContVarJump, pos] <- L[namesContVarJump, pos]
+ ntot <- sampling at n
+ dL <- L
+ dL[object at Ppr@counting.var,c((pos0+1):pos)]<-JUMP[object at Ppr@counting.var,c((pos0+1):pos)]
+
+ X_mat <- matrix(0, length(object at model@solve.variable),
+ ntot)
+ rownames(X_mat) <- object at model@solve.variable
+
+ dummyX <- simulate(object at model, true.parameter = true.parameter,
+ increment.W = if(is.matrix(W[,1:pos])){W[,1:pos]}else{t(as.matrix(W[,1:pos]))},
+ increment.L = if(is.matrix(dL[,1:pos])){dL[,1:pos]}else{t(as.matrix(dL[,1:pos]))},
+ sampling = setSampling(Terminal = t_1,
+ n = t_1/sampling at delta))
+ X_mat[,1:pos] <- t(dummyX at data@original.data)[,-1]
+
+ t_jump <- t_1
+ if(length(object at Kernel@variable.Integral at var.dx)==1){
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 1, diff)
+ }else{
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 2, diff)
+ }
+
+
+
+ Index <- matrix(c(1:prod(object at Kernel@Integrand at dimIntegrand)),
+ nrow = object at Kernel@Integrand at dimIntegrand[1],
+ ncol = object at Kernel@Integrand at dimIntegrand[2])
+
+ assign(object at Kernel@variable.Integral at var.time,
+ sampling at grid[[1]][c((pos0+1):(pos))],
+ envir = my.env)
+
+ assign(object at gFun@param at time.var, t_1, envir = my.env)
+ for(i in c(1:object at Kernel@Integrand at dimIntegrand[2])){
+ assign(object at Kernel@variable.Integral at var.dx[i],
+ as.numeric(Comulat.dx[,i]),
+ envir = my.env)
+ }
+ KernDum <- list()
+ for(i in c(1:object at Kernel@Integrand at dimIntegrand[1])){
+ dumKern <- expression()
+ for(j in c(1:object at Kernel@Integrand at dimIntegrand[2])){
+ id <- as.numeric(Index[i,j])
+ dumKern <- c(dumKern,object at Kernel@Integrand at IntegrandList[[id]])
+
+ }
+ KernDum[[i]] <- dumKern
+ }
+
+
+ udumm <- numeric(length = numbCountVar)
+ names(udumm) <- object at Kernel@variable.Integral at var.dx
+
+ assign("L",dL,envir = my.env)
+ pos0 <- pos
+ assign("pos0", pos, envir = my.env)
+ assign("old_u",t_1, envir = my.env)
+
+ while(t_jump<Time){
+
+
+ oldt_1<-t_1
+ for(i in c(1:numbCountVar)){
+ assign("gFun", object at gFun@formula[[i]], envir=my.env)
+ assign("Kern", KernDum[[i]], envir=my.env)
+ my.env$utrue<-0
+ while(my.env$utrue<oldt_1){
+ assign("U",runif(1),envir = my.env)
+ optim((t_1+2*my.env$dt),Fun2,method = "Nelder-Mead",
+ env=my.env)$par
+ u[i] <- as.numeric(my.env$utrue)
+ }
+ }
+
+ t_1 <- min(u)
+
+ condt1<- u%in%t_1
+ namesContVarJump <- names(u[condt1])
+
+ mypos<-sum(sampling at grid[[1]][-1]<=t_1)
+ if((pos0+1)<mypos){
+ pos<-sum(sampling at grid[[1]][-1]<=t_1)
+ t_jump<- t_1
+ t_1 <- sampling at grid[[1]][-1][pos]
+ recordTime<-c(recordTime,t_1)
+
+
+ #if(t_1!=sampling at Terminal){
+
+ pos <- min(pos,dim(L)[2])
+ JUMP[namesContVarJump, pos] <- L[namesContVarJump, pos]
+ dL[object at Ppr@counting.var,c((pos0+1):pos)]<-JUMP[object at Ppr@counting.var,c((pos0+1):pos)]
+ aa<-setSampling(Terminal = (t_1-my.env$old_u),
+ n = length((pos0+1):pos))
+ dummyX <- simulate(object at model, true.parameter = true.parameter,
+ increment.W = if(is.matrix(W[,(pos0+1):pos])){W[,(pos0+1):pos]}else{t(as.matrix(W[,(pos0+1):pos]))},
+ increment.L = if(is.matrix(dL[,(pos0+1):pos])){dL[,(pos0+1):pos]}else{t(as.matrix(dL[,(pos0+1):pos]))},
+ sampling = aa,
+ xinit=X_mat[,(pos0)])
+ X_mat[,(pos0+1):pos] <- t(dummyX at data@original.data)[,-1]
+ if(length(object at Kernel@variable.Integral at var.dx)==1){
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 1, diff)
+ }else{
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 2, diff)
+ }
+ if(!is.matrix(Comulat.dx)){
+ Comulat.dx <-t(as.matrix(Comulat.dx))
+ }
+
+ Index <- matrix(c(1:prod(object at Kernel@Integrand at dimIntegrand)),
+ nrow = object at Kernel@Integrand at dimIntegrand[1],
+ ncol = object at Kernel@Integrand at dimIntegrand[2])
+ assign(object at Kernel@variable.Integral at var.time,
+ sampling at grid[[1]][c((pos0+1):(pos))],
+ envir = my.env)
+ assign(object at gFun@param at time.var, t_1, envir = my.env)
+ for(i in c(1:object at Kernel@Integrand at dimIntegrand[2])){
+
+ assign(object at Kernel@variable.Integral at var.dx[i],
+ as.numeric(Comulat.dx[,i]),
+ envir = my.env)
+ }
+ pos0<-pos
+ assign("pos0", pos, envir = my.env)
+ assign("old_u",t_1, envir = my.env)
+
+ #}
+ }
+ assign("L",dL,envir = my.env)
+ }
+ X_mat[namesContVarJump,pos]<-X_mat[namesContVarJump,pos]
+ res.dum <- list(X_mat=X_mat,timeJump = recordTime, grid=sampling)
+
+ solve.variable <-unique(c(object at model@solve.variable))
+ N.VarPPr<-length(solve.variable)
+
+ dummy.mod <- setModel(drift=rep("0",N.VarPPr),
+ diffusion = NULL, jump.coeff = diag(rep("1",N.VarPPr)),
+ measure = object at Ppr@Info.measure$measure,
+ measure.type = object at Ppr@Info.measure$type,
+ solve.variable = solve.variable, xinit=c(object at model@xinit))
+
+ mynewincr <- if(is.matrix(res.dum$X_mat)){t(as.matrix(apply(cbind(0,res.dum$X_mat),1,diff)))}else{apply(cbind(0,res.dum$X_mat),1,diff)}
+
+ interResMod <- simulate(object = dummy.mod,
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
+
+ resGfun<-new("yuima.Output",
+ Output = object at gFun,
+ yuima=setYuima(model=dummy.mod,sampling = sampling))
+
+ interResGfun <- simulate(object = resGfun,
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
+ dummyObject <- object at Kernel
+ dummyObject at variable.Integral@out.var <-object at Ppr@additional.info
+ resInt <- new("yuima.Integral",
+ Integral = dummyObject,
+ yuima = setYuima(model=dummy.mod,sampling = sampling))
+
+ interResInt <- simulate(object = resInt,
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
+ DataIntensity <- interResGfun at data@original.data + interResInt at data@original.data
+ InterMDia<-zoo(interResMod at data@original.data, order.by = index(DataIntensity))
+ Alldata <-merge(InterMDia,DataIntensity)
+ colnames(Alldata)<-c(solve.variable,object at Ppr@additional.info)
+ # for(i in c(1:N.VarPPr)){
+ # assign(solve.variable[i],interRes at data@original.data[,i],envir=my.env)
+ # }
+ # dummy<-NULL
+ # for(t in c(1:length(object at Ppr@additional.info))){
+ # dummy <-eval(object at gFun)
+ # assign(object at Ppr@additional.info[[]])
+ # }
+ object at data<-setData(Alldata)
+ return(object)
+ }
+
+
+
+
+# simOzaki.aux<-function(gFun,a,cCoeff, Time, numJump){
+# t_k<-0
+# N<-0
+# S<-1
+#
+# T_k<-c(t_k)
+#
+# N_k<-c(N)
+# U<-runif(1)
+# t_k <- -log(U)/gFun
+# if(t_k<Time){
+# T_k<-c(T_k,t_k)
+# N<-N+numJump
+# N_k<-c(N_k, N)
+# }
+# while(t_k<=Time){
+# U<-runif(1)
+# optim.env<-new.env()
+# assign("U",U,envir=optim.env)
+# assign("t_k",t_k,envir=optim.env)
+# assign("c",cCoeff,envir=optim.env)
+# assign("a",a,envir=optim.env)
+# assign("S",S,envir=optim.env)
+# assign("gFun",gFun,envir=optim.env)
+#
+# min<-function(u,env){
+# U<-env$U
+# t_k<-env$t_k
+# c<-env$c
+# a<-env$a
+# S<-env$S
+# gFun<-env$gFun
+# y<-(log(U)+gFun*(u-t_k)+c/a*S*(1-exp(-a*(u-t_k))))^2
+# }
+# y<-optim(par=t_k,min, env=optim.env )$par
+# S<- exp(-a*(y-t_k))*S+1
+# t_k<-y
+# T_k<-c(T_k,t_k)
+# N<-N+numJump
+# N_k<-c(N_k, N)
+# }
+# return(list(T_k=T_k,N_k=N_k))
+# }
Modified: pkg/yuima/R/simulateMultiProcess.R
===================================================================
--- pkg/yuima/R/simulateMultiProcess.R 2016-07-08 04:44:22 UTC (rev 456)
+++ pkg/yuima/R/simulateMultiProcess.R 2016-07-09 20:57:55 UTC (rev 457)
@@ -218,6 +218,7 @@
if(sdeModel at measure.type=="CP"){
intens <- as.character(sdeModel at measure$intensity)
dens <- as.character(sdeModel at measure$df$expr)
+
dumCP <- setPoisson(intensity = intens, df = dens,
dimension = length(sdeModel at jump.coeff[[1]]))
dummSamp <- yuima at sampling
@@ -329,7 +330,10 @@
# 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)
Added: pkg/yuima/man/info.Ppr.Rd
===================================================================
--- pkg/yuima/man/info.Ppr.Rd (rev 0)
+++ pkg/yuima/man/info.Ppr.Rd 2016-07-09 20:57:55 UTC (rev 457)
@@ -0,0 +1,11 @@
+\name{info.Ppr}
+\docType{class}
+\alias{info.Ppr}
+\alias{info.Ppr-class}
+\alias{initialize,info.Ppr-method}
+
+
+\title{Class for information about Point Process}
+\description{
+Auxiliar class for definition of an object of class \code{\link{yuima.Ppr}} and \code{\link{yuima.Hawkes}}. see the documentation for more details.
+}
Added: pkg/yuima/man/setHawkes.Rd
===================================================================
--- pkg/yuima/man/setHawkes.Rd (rev 0)
+++ pkg/yuima/man/setHawkes.Rd 2016-07-09 20:57:55 UTC (rev 457)
@@ -0,0 +1,43 @@
+\name{setHawkes}
+\alias{setHawkes}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Constructor of Hawkes model}
+\description{
+...
+}
+\usage{
+setHawkes(lower.var = "0", upper.var = "t", var.dt = "s",
+ process = "N", dimension = 1, intensity = "lambda",
+ ExpKernParm1 = "c", ExpKernParm2 = "a", const = "nu",
+ measure = NULL, measure.type = NULL)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{lower.var}{...}
+ \item{upper.var}{...}
+ \item{var.dt}{...}
+ \item{process}{...}
+ \item{dimension}{...}
+ \item{intensity}{...}
+ \item{ExpKernParm1}{...}
+ \item{ExpKernParm2}{...}
+ \item{const}{...}
+ \item{measure}{...}
+ \item{measure.type}{...}
+}
+\details{...}
+\value{...
+}
+\references{
+...}
+\author{
+...}
+\note{
+...}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{...}
+%\examples{}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 457
More information about the Yuima-commits
mailing list