[Yuima-commits] r43 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Dec 12 03:01:09 CET 2009
Author: hinohide
Date: 2009-12-12 03:01:08 +0100 (Sat, 12 Dec 2009)
New Revision: 43
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/AllClasses.R
pkg/yuima/R/sim.euler.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/simulate.Rd
pkg/yuima/man/yuima.model-class.Rd
Log:
Modified interface for Levy process modeling
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/DESCRIPTION 2009-12-12 02:01:08 UTC (rev 43)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.74
-Date: 2009-12-09
+Version: 0.0.75
+Date: 2009-12-12
Depends: methods, zoo
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/R/AllClasses.R 2009-12-12 02:01:08 UTC (rev 43)
@@ -28,7 +28,8 @@
equation.number="numeric",
dimension="numeric",
solve.variable="character",
- xinit="numeric"
+ xinit="numeric",
+ J.flag="logical"
)
)
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/R/sim.euler.R 2009-12-12 02:01:08 UTC (rev 43)
@@ -125,7 +125,7 @@
}
return(tmp)
}
-
+
if(sdeModel at measure.type == "CP"){ ##:: Compound-Poisson type
eta0 <- eval(sdeModel at measure$intensity)
##:: get lambda from nu()
@@ -146,13 +146,13 @@
}
##:: make expression to create iid rand J
- if(grep("^[dexp|dnorm|dgamma]", sdeModel at measure$df$expr)) {
+ if(grep("^[dexp|dnorm|dgamma]", 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.")
}
- randJ <- eval(F) ## this expression is evaluated locally not in the yuimaEnv
+ randJ <- eval(F) ## this expression is evaluated locally not in the yuimaEnv
j <- 1
for(i in 1:division){
if(JAMP==FALSE || sum(i==ij)==0){
@@ -161,16 +161,21 @@
J <- eta0*randJ[j]/lambda
j <- j+1
##cat(paste(J,"\n"))
- ##Pi <- zeta(dX,J)
+ ##Pi <- zeta(dX, J)
assign(sdeModel at jump.variable, J, env)
- ##Pi <- p.b.j(t=i*delta,X=dX) %*% J
- Pi <- p.b.j(t=i*delta, X=dX)
+
+ if(sdeModel at J.flag){
+ J <- 1
+ }
+
+ Pi <- p.b.j(t=i*delta,X=dX) * J
+ ##Pi <- p.b.j(t=i*delta, X=dX)
}
dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi
X_mat[, i+1] <- dX
}
tsX <- ts(data=t(X_mat), deltat=delta, start=0)
-
+ ##::end CP
}else if(sdeModel at measure.type=="code"){ ##:: code type
##:: Jump terms
code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
@@ -194,18 +199,22 @@
##:: calcurate difference eq.
for(i in 1:division){
- assign(sdeModel at jump.variable, dZ[i])
+ assign(sdeModel at jump.variable, dZ[i], env)
+
+ if(sdeModel at J.flag){
+ dZ[i] <- 1
+ }
+
dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +p.b.j(t=i*delta, X=dX) * dZ[i]
X_mat[, i+1] <- dX
}
tsX <- ts(data=t(X_mat), deltat=delta, start=0)
-
+ ##::end code
}else{
cat(paste("Type \"", sdeModel at measure.type, "\" not supported yet.\n", sep=""))
return(NULL)
}
- }
-
+ }##::end levy
yuimaData <- setData(original.data=tsX)
return(yuimaData)
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/R/yuima.model.R 2009-12-12 02:01:08 UTC (rev 43)
@@ -31,7 +31,8 @@
equation.number,
dimension,
solve.variable,
- xinit){
+ xinit,
+ J.flag){
.Object at drift <- drift
.Object at diffusion <- diffusion
.Object at hurst <- hurst
@@ -47,6 +48,7 @@
.Object at dimension <- dimension
.Object at solve.variable <- solve.variable
.Object at xinit <- xinit
+ .Object at J.flag <- J.flag
return(.Object)
})
@@ -149,7 +151,7 @@
measure.par <- unique( c( all.vars(MEASURE$intensity), all.vars(MEASURE$df$expr) ) )
##measure.par$intensity <- unique(all.vars(MEASURE$intensity))
##::end check df name ####################
-
+ ##::end CP
}else if(measure.type=="code"){ ##::code
if(length(measure)!=1){
cat(paste("\nlength of measure must be one on type", measure.type, ".\n"))
@@ -187,7 +189,7 @@
measure.par <- unique(all.vars(MEASURE$df$expr))
##::end check df name ####################
-
+ ##::end code
}else if(measure.type=="density"){ ##::density
if(length(measure)!=1){
cat(paste("\nlength of measure must be one on type", measure.type, ".\n"))
@@ -227,7 +229,7 @@
measure.par <- unique(all.vars(MEASURE[[names(measure)]]$expr))
##::end check df name ####################
-
+ ##::end density
}else{ ##::else
cat(paste("\nmeasure type", measure.type, "isn't supported.\n"))
return(NULL)
@@ -345,7 +347,11 @@
}
##:: get parameters in jump expression
+ J.flag <- FALSE
jump.par <- unique(all.vars(JUMP))
+ if(length(na.omit(match(jump.par, jump.variable)))){
+ J.flag <- TRUE
+ }
jump.idx <- as.numeric(na.omit(match(c(state.variable, time.variable, jump.variable, solve.variable), jump.par)))
if(length(jump.idx)>0){
jump.par <- jump.par[-jump.idx]
@@ -356,7 +362,7 @@
if(length(measure.idx)>0){
measure.par <- measure.par[-measure.idx]
}
-
+
##:: order parameters for 'yuima.pars'
##id1 <- which(diff.par %in% drift.par)
##id2 <- which(drift.par %in% diff.par)
@@ -401,7 +407,7 @@
length(tmppar at measure)
),
solve.variable= solve.variable,
- xinit= xinit)
-
+ xinit= xinit,
+ J.flag <- J.flag)
return(tmp)
}
Modified: pkg/yuima/man/simulate.Rd
===================================================================
--- pkg/yuima/man/simulate.Rd 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/man/simulate.Rd 2009-12-12 02:01:08 UTC (rev 43)
@@ -28,10 +28,6 @@
It simulates solutions of stochastic differential equations with Gaussian noise,
fractional Gaussian noise awith/without jumps.
-
-By default, \code{methodfGn = "Cholesky"} uses the exact method based on the
-Cholesky decomposition. This method may not work, i.e. the computation is
-stopped, for too long grids due to memory constraints.
}
\value{
\item{yuima}{A multi-dimensional \code{yuima} object(time series object).}
@@ -67,22 +63,11 @@
ou1 <- setYuima(model=mod1, sampling=samp)
# Solve SDEs using Euler-Maruyama method.
-ou1 <- simulate(ou1, xinit=1, true.p = list(theta=-1))
+ou1 <- simulate(ou, xinit=1, true.p = list(theta=-1))
plot(ou1)
-# Path-simulation for 1-dim diffusion process with fractional noise.
-# dXt = theta*Xt*dt + dWt^h
-mod2 <- setModel(drift="theta*y", diffusion=1, solve.variable=c("y"), hurst=0.3)
-str(mod2)
-ou2 <- setYuima(model=mod2, sampling=samp)
-# Solve SDEs using Euler-Maruyama method.
-ou2 <- simulate(ou2, xinit=1, true.p = list(theta=-1))
-plot(ou2)
-
-
-
# A multi-dimensional (correlated) diffusion process.
# To describe the following model:
# X=(X1,X2,X3); dXt = U(t,Xt)dt + V(t)dWt
@@ -210,8 +195,9 @@
##:: sample for Levy process ("CP" type)
+## specify the jump term as c(x,t)dz
obj.model <- setModel(drift=c("-theta*x"), diffusion="sigma",
-jump.coeff="z", measure=list(intensity="1", df=list("dnorm(z, 0, 1)")),
+jump.coeff="1", measure=list(intensity="1", df=list("dnorm(z, 0, 1)")),
measure.type="CP", solve.variable="x")
##:: Parameters
@@ -228,15 +214,25 @@
obj.sampling <- setSampling(Terminal=T, division=division)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
X <- simulate(obj.yuima, xinit=xinit, true.parameter=list(theta=theta, sigma=sigma))
+X11()
plot(X)
+##:: sample for Levy process ("CP" type)
+## specify the jump term as c(x,t,z)
+obj.model <- setModel(drift=c("-theta*x"), diffusion="sigma",
+jump.coeff="z",J.flag=TRUE, measure=list(intensity="1", df=list("dnorm(z, 0, 1)")),
+measure.type="CP", solve.variable="x")
+
+
+
##:: sample for Levy process ("code" type)
## dX_{t} = -x dt + dZ_t
obj.model <- setModel(drift="-x", xinit=1, jump.coeff="1", measure.type="code", measure=list(df="rIG(z, 1, 0.1)"))
obj.sampling <- setSampling(Terminal=10, division=10000)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
result <- simulate(obj.yuima)
+X11()
plot(result)
##:: sample for multidimensional Levy process ("code" type)
Modified: pkg/yuima/man/yuima.model-class.Rd
===================================================================
--- pkg/yuima/man/yuima.model-class.Rd 2009-12-09 05:08:24 UTC (rev 42)
+++ pkg/yuima/man/yuima.model-class.Rd 2009-12-12 02:01:08 UTC (rev 43)
@@ -43,6 +43,7 @@
the stochastic differential equation has to be solved.}
\item{\code{xinit}:}{contains the initial value of the stochastic
differential equation.}
+ \item{\code{J.flag}:}{wheather jump.coeff include jump.variable.}
}
}
\author{The YUIMA Project Team}
More information about the Yuima-commits
mailing list