[Yuima-commits] r328 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 23 10:20:27 CEST 2014
Author: iacus
Date: 2014-09-23 10:20:27 +0200 (Tue, 23 Sep 2014)
New Revision: 328
Added:
pkg/yuima/R/setPoisson.R
pkg/yuima/R/simCP.R
Modified:
pkg/yuima/R/AllClasses.R
pkg/yuima/R/qmle.R
pkg/yuima/R/sim.euler.R
pkg/yuima/R/simulate.R
pkg/yuima/R/yuima.R
pkg/yuima/R/yuima.model.R
Log:
poisson simulator
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/AllClasses.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -56,7 +56,10 @@
representation(info="carma.info"),
contains="yuima.model")
+# Class Compound Poisson
+setClass("yuima.poisson", contains="yuima.model")
+
# Class 'yuima.data'
# we want yuimaS4 to use any class of data as input
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/qmle.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -105,6 +105,14 @@
### also, I am using the same interface of optim to specify upper and lower bounds
### S.M.I. 22/06/2010
+is.Poisson <- function(obj){
+ if(is(obj,"yuima"))
+ return(is(obj at model, "yuima.poisson"))
+ if(is(obj,"yuima.model"))
+ return(is(obj, "yuima.poisson"))
+ return(FALSE)
+}
+
is.CARMA <- function(obj){
if(is(obj,"yuima"))
return(is(obj at model, "yuima.carma"))
Added: pkg/yuima/R/setPoisson.R
===================================================================
--- pkg/yuima/R/setPoisson.R (rev 0)
+++ pkg/yuima/R/setPoisson.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -0,0 +1,72 @@
+
+
+setMethod("initialize", "yuima.poisson",
+ function(.Object,
+ drift = expression() ,
+ diffusion = list() ,
+ hurst = 0.5,
+ jump.coeff = expression(),
+ measure=list(),
+ measure.type=character(),
+ parameter = new("model.parameter"),
+ state.variable = "x",
+ jump.variable = "z",
+ time.variable = "t",
+ noise.number = numeric(),
+ equation.number = numeric(),
+ dimension = numeric(),
+ solve.variable = character(),
+ xinit = expression(),
+ J.flag = logical()){
+ .Object at drift <- drift
+ .Object at diffusion <- diffusion
+ .Object at hurst <- 0.5
+ .Object at jump.coeff <- jump.coeff
+ .Object at measure <- measure
+ .Object at measure.type <- measure.type
+ .Object at parameter <- parameter
+ .Object at state.variable <- state.variable
+ .Object at jump.variable <- jump.variable
+ .Object at time.variable <- time.variable
+ .Object at noise.number <- noise.number
+ .Object at equation.number <- equation.number
+ .Object at dimension <- dimension
+ .Object at solve.variable <- solve.variable
+ .Object at xinit <- xinit
+ .Object at J.flag <- J.flag
+ return(.Object)
+ })
+
+setPoisson <- function(intensity=1, df=NULL, scale=1, ...){
+
+ poisson.model <- setModel(drift="0", diffusion="0",
+ jump.coeff=as.expression(scale),
+ measure=list(intensity=intensity, df=df),
+ measure.type="CP",...)
+
+
+ CPoisson <- new("yuima.poisson",
+ drift=poisson.model at drift,
+ diffusion = poisson.model at diffusion,
+ hurst=poisson.model at hurst,
+ jump.coeff=poisson.model at jump.coeff,
+ measure=poisson.model at measure,
+ measure.type=poisson.model at measure.type,
+ parameter=poisson.model at parameter,
+ state.variable=poisson.model at state.variable,
+ jump.variable=poisson.model at jump.variable,
+ time.variable=poisson.model at time.variable,
+ noise.number = poisson.model at noise.number,
+ equation.number = poisson.model at equation.number,
+ dimension = poisson.model at dimension,
+ solve.variable=poisson.model at solve.variable,
+ xinit=poisson.model at xinit,
+ J.flag = poisson.model at J.flag
+ )
+
+ return(CPoisson)
+
+}
+
+dconst <- function(x,k=1) k*(x==k)
+rconst <- function(n,k=1) rep(k,n)
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/sim.euler.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -56,7 +56,8 @@
##:: check if DRIFT and/or DIFFUSION has values
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)))
-
+ #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)){
@@ -74,14 +75,16 @@
tmp[i,j+1] <- eval(V[[i]][j],env)
}
}
- }else{ ##:: no drift term (faster)
- tmp <- matrix(0, d.size, r.size)
- for(i in 1:d.size){
- for(j in 1:r.size){
- tmp[i,j] <- 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)
}
@@ -96,7 +99,7 @@
if(!length(sdeModel at jump.coeff)){ ##:: Wiener Proc
##:: using Euler-Maruyama method
- if(var.in.diff){ ##:: diffusions have state variables
+ 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]
@@ -127,19 +130,22 @@
dim(pbdata)<-(c(r.size+1, d.size*t.size))
}else{
pbdata <- matrix(0, d.size*r.size, t.size)
- for(i in 1:d.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)
- }
- }
+ 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){
- dX <- dX + pbdata[((i-1)*d.size+1):(i*d.size), ] %*% dW[, i]
+ if(!is.Poisson(sdeModel))
+ dX <- dX + pbdata[((i-1)*d.size+1):(i*d.size), ] %*% dW[, i]
X_mat[, i+1] <- dX
}
}
@@ -205,11 +211,11 @@
}
##:: make expression to create iid rand J
- if(grep("^[dexp|dnorm|dgamma]", sdeModel at measure$df$expr)){
+ if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel at 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 at measure$df$expr, perl=TRUE)))
}else{
- stop("Sorry. CP only supports dexp, dnorm and dgamma yet.")
+ stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
}
##:: delete 2010/09/13 for simulate func bug fix by s.h
Added: pkg/yuima/R/simCP.R
===================================================================
--- pkg/yuima/R/simCP.R (rev 0)
+++ pkg/yuima/R/simCP.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -0,0 +1,86 @@
+simCP<-function(xinit,yuima,env){
+
+
+ sdeModel<-yuima at model
+
+ modelstate <- sdeModel at solve.variable
+ modeltime <- sdeModel at time.variable
+ Terminal <- yuima at sampling@Terminal[1]
+ Initial <- yuima at sampling@Initial[1]
+
+ 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,length(xinit))
+ }
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i],dummy.val[i] ,envir=env)
+ }
+
+ }
+
+
+### Simulation of CP using Lewis' method
+
+##:: Levy
+ JP <- eval(sdeModel at jump.coeff, envir=env)
+ mu.size <- length(JP)
+
+
+ #assign(sdeModel at measure$intensity, env) ## intensity param
+ .CPintensity <- function(.t) {
+ assign(modeltime, .t, envir=env)
+ eval(sdeModel at measure$intensity, envir=env)
+ }
+
+
+ dummyList<-as.list(env)
+
+ lgth.meas<-length(yuima at model@parameter at measure)
+ if(lgth.meas>1){
+ for(i in c(2:lgth.meas)){
+ idx.dummy<-yuima at model@parameter at measure[i]
+ assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
+ }
+ }
+
+ # we use Lewis' acceptance/rejection method
+
+ if(grep("^[dexp|dnorm|dgamma|dconst]", sdeModel at 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 at measure$df$expr, perl=TRUE)))
+ } else{
+ stop("Sorry. CP only supports dconst, dexp, dnorm and dgamma yet.")
+ }
+
+ ell <- optimize(f=.CPintensity, interval=c(Initial, Terminal), maximum = TRUE)$objective
+ ellMax <- ell * 1.01
+
+
+ time <- Initial
+ E <- Initial
+
+ while(time < Terminal) {
+ time <- time - 1/ellMax * log(runif(1))
+ if(runif(1) < .CPintensity(time)/ellMax)
+ E <- c(E, time)
+ }
+ N_sharp <- length(E)-1
+
+ 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, envir=F.env) ## this expression is evaluated in the F.env
+
+ randJ <- JP[1]*randJ
+
+
+ CP <- cumsum(c(eval(xinit[1],envir=env), randJ))
+ # print(CP)
+ tsX <- zoo(x=CP, order.by=E)
+ yuimaData <- setYuima(data=setData(tsX))
+ yuimaData <- subsampling(yuimaData, sampling=yuima at sampling)
+ return(yuimaData at data)
+}
Property changes on: pkg/yuima/R/simCP.R
___________________________________________________________________
Added: svn:executable
+ *
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/simulate.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -218,15 +218,23 @@
} else {
delta<-Terminal/n
- dW <- rnorm(n * r.size, 0, sqrt(delta))
- dW <- matrix(dW, ncol=n, nrow=r.size,byrow=TRUE)
+ 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
}
- yuima at data <- euler(xinit, yuima, dW, yuimaEnv)
+ if(is.Poisson(sdeModel)){
+ yuima at data <- simCP(xinit, yuima, yuimaEnv)
+ } else {
+ yuima at data <- euler(xinit, yuima, dW, yuimaEnv)
+ }
for(i in 1:length(yuima at data@zoo.data))
index(yuima at data@zoo.data[[i]]) <- yuima at sampling@grid[[1]] ## to be fixed
Modified: pkg/yuima/R/yuima.R
===================================================================
--- pkg/yuima/R/yuima.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/yuima.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -293,6 +293,7 @@
is.fracdiff <- FALSE
is.jumpdiff <- FALSE
is.carma <- FALSE
+ is.poisson <- is.Poisson(mod)
if(length(mod at drift)>0) has.drift <- TRUE
if(length(mod at diffusion)>0) has.diff <- TRUE
@@ -316,7 +317,7 @@
is.carma <- TRUE
if( is.wienerdiff | is.fracdiff | is.jumpdiff ){
- if( is.wienerdiff & ! is.carma){
+ if( is.wienerdiff & ! is.carma & !is.poisson){
cat("\nDiffusion process")
if( is.fracdiff ){
if(!is.na(mod at hurst)){
@@ -330,17 +331,20 @@
}
if(is.carma)
cat(sprintf("\nCarma process p=%d, q=%d", mod at info@p, mod at info@q))
-
+ if(is.poisson)
+ cat("\nCompound Poisson process")
+
if( is.jumpdiff ){
- if( is.wienerdiff | is.carma ){
+ if( (is.wienerdiff | is.carma) & !is.poisson ){
cat(" with Levy jumps")
} else {
+ if(!is.poisson)
cat("Levy jump process")
}
}
cat(sprintf("\nNumber of equations: %d", mod at equation.number))
- if(is.wienerdiff | is.fracdiff)
+ if((is.wienerdiff | is.fracdiff) & !is.poisson)
cat(sprintf("\nNumber of Wiener noises: %d", mod at noise.number))
if(is.jumpdiff)
cat(sprintf("\nNumber of Levy noises: %d", 1))
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2014-09-08 10:18:54 UTC (rev 327)
+++ pkg/yuima/R/yuima.model.R 2014-09-23 08:20:27 UTC (rev 328)
@@ -139,7 +139,7 @@
##::end initialize objects ########
##::make type function list ####
- CPlist <- c("dnorm", "dgamma", "dexp")
+ CPlist <- c("dnorm", "dgamma", "dexp", "dconst")
codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rngamma", "rstable")
##::end make type function list ####
More information about the Yuima-commits
mailing list