[Yuima-commits] r333 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 28 21:43:13 CEST 2014
Author: iacus
Date: 2014-09-28 21:43:13 +0200 (Sun, 28 Sep 2014)
New Revision: 333
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/AllClasses.R
pkg/yuima/R/qmle.R
pkg/yuima/R/rng.R
pkg/yuima/R/setPoisson.R
pkg/yuima/R/sim.euler.R
pkg/yuima/R/simCP.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/setPoisson.Rd
pkg/yuima/man/simulate.Rd
Log:
added multidim levy
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/DESCRIPTION 2014-09-28 19:43:13 UTC (rev 333)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.33
-Date: 2014-09-23
+Version: 1.0.34
+Date: 2014-09-29
Depends: methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/AllClasses.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -19,7 +19,8 @@
setClass("yuima.model",representation(drift="expression",
diffusion="list",
hurst="ANY",
- jump.coeff="expression",
+ jump.coeff="list",
+#jump.coeff="expression",
measure="list",
measure.type="character",
parameter="model.parameter",
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/qmle.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -87,7 +87,7 @@
for(r in 1:r.size){
#for(d.tmp in 1:d){
if(d.size==1){
- measure[1,r,] <- eval(JUMP[r],envir=tmp.env)
+ measure[1,r,] <- eval(JUMP[[r]],envir=tmp.env)
}else{
for(d in 1:d.size){
measure[d,r,] <- eval(JUMP[[d]][r],envir=tmp.env)
Modified: pkg/yuima/R/rng.R
===================================================================
--- pkg/yuima/R/rng.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/rng.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -28,7 +28,7 @@
if(length(mu)!=length(beta)){
stop("Error: wrong input dimension.")
}
- if(missing(Lambda)){
+ if(missing(Lambda) | is.na(Lambda)){
## univariate case
if(length(mu)!=1 || length(beta)!=1){
stop("Error: wrong input dimension.")
@@ -175,11 +175,12 @@
rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0,Lambda){
## Error check
+ #print(match.call())
if(length(mu)!=length(beta)){
stop("Error: wrong input dimension.")
}
- if(missing(Lambda)){
+ if(missing(Lambda) | is.na(Lambda)){
## univariate case
gamma <- sqrt(alpha^2 - beta^2)
if(gamma <0){
@@ -215,6 +216,7 @@
z <- mu + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% beta,x),length(beta),x)+t(matrix(rep(sqrt(tau),length(beta)),x,length(beta))) * (sqrt.L %*% t(matrix(eta,x,length(beta))))
X <- z
+ # print(str(X))
return(X)
}
Modified: pkg/yuima/R/setPoisson.R
===================================================================
--- pkg/yuima/R/setPoisson.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/setPoisson.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -37,7 +37,7 @@
return(.Object)
})
-setPoisson <- function(intensity=1, df=NULL, scale=1, ...){
+setPoisson <- function(intensity=1, df=NULL, scale=1, dimension=1, ...){
poisson.model <- setModel(drift="0", diffusion="0",
jump.coeff=as.expression(scale),
@@ -57,8 +57,8 @@
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,
+ equation.number = dimension,
+ dimension = dimension,
solve.variable=poisson.model at solve.variable,
xinit=poisson.model at xinit,
J.flag = poisson.model at J.flag
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/sim.euler.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -1,5 +1,4 @@
euler<-function(xinit,yuima,dW,env){
-
sdeModel<-yuima at model
@@ -96,7 +95,7 @@
}
- if(!length(sdeModel at jump.coeff)){ ##:: Wiener Proc
+ if(!length(sdeModel at measure.type)){ ##:: Wiener Proc
##:: using Euler-Maruyama method
if(var.in.diff & (!is.Poisson(sdeModel))){ ##:: diffusions have state variables and it is not Poisson
@@ -154,20 +153,39 @@
}else{ ##:: Levy
JP <- sdeModel at jump.coeff
mu.size <- length(JP)
-
+ # cat("\n Levy\n")
##:: function to solve c(x,z)
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)
- tmp <- numeric(d.size)
- for(i in 1:d.size){
- tmp[i] <- eval(JP[i], env)
+ # tmp <- numeric(d.size)
+ j.size <- length(JP[[1]])
+ tmp <- matrix(0, mu.size, j.size)
+ # cat("\n inside\n")
+ #print(dim(tmp))
+ for(i in 1:mu.size){
+ for(j in 1:j.size){
+ tmp[i,j] <- eval(JP[[i]][j],env)
+ }
+ # tmp[i] <- eval(JP[i], env)
}
return(tmp)
}
+ # print(ls(env))
+ ### WHY I AM DOING THIS?
+ # tmp <- matrix(0, d.size, r.size)
+ #
+ #for(i in 1:d.size){
+ # for(j in 1:r.size){
+ # cat("\n here\n")
+ # tmp[i,j] <- eval(V[[i]][j],env)
+ # } # for j
+ # }
+ ###
+
if(sdeModel at measure.type == "CP"){ ##:: Compound-Poisson type
##:: delete 2010/09/13 for simulate func bug fix by s.h
@@ -185,10 +203,14 @@
#lambda <- integrate(tmp.expr, 0, Inf)$value * eta0 ##bug:2013/10/28
dummyList<-as.list(env)
+ # print(str(dummyList))
+ #print(str(idx.dummy))
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]
+ #print(i)
+ #print(yuima at model@parameter at measure[i])
assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
}
}
@@ -254,22 +276,31 @@
##:: Jump terms
code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
+ #print(args)
dZ <- switch(code,
- rNIG=paste("rNIG(n, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta)"),
+ rNIG=paste("rNIG(n, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta, ", args[6],")"),
rIG=paste("rIG(n,", args[2], "*delta, ", args[3], ")"),
rgamma=paste("rgamma(n, ", args[2], "*delta, ", args[3], ")"),
rbgamma=paste("rbgamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], "*delta, ", args[5], ")"),
## rngamma=paste("rngamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
- rngamma=paste("rngamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta)"),
+ rngamma=paste("rngamma(n, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta,", args[6],")"),
## rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], "*delta^(1/",args[2],"), ", args[5], "*delta)")
)
dummyList<-as.list(env)
+ #print(str(dummyList))
lgth.meas<-length(yuima at model@parameter at measure)
+ #print(lgth.meas)
if(lgth.meas!=0){
for(i in c(1:lgth.meas)){
+ #print(i)
+ #print(yuima at model@parameter at measure[i])
idx.dummy<-yuima at model@parameter at measure[i]
- assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
+ #print(str(idx.dummy))
+ assign(idx.dummy,dummyList[[idx.dummy]])
+ #print(str(idx.dummy))
+ #print(str(dummyList[[idx.dummy]]))
+ #print(get(idx.dummy))
}
}
@@ -279,15 +310,28 @@
}
dZ <- eval(parse(text=dZ))
##:: calcurate difference eq.
-
+ #print(str(dZ))
+ if(is.null(dim(dZ)))
+ dZ <- matrix(dZ,nrow=1)
+ # print(dim(dZ))
+ # print(str(sdeModel at jump.variable))
for(i in 1:n){
- assign(sdeModel at jump.variable, dZ[i], env)
+ assign(sdeModel at jump.variable, dZ[,i], env)
if(sdeModel at J.flag){
- dZ[i] <- 1
+ dZ[,i] <- 1
}
-
- dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +p.b.j(t=i*delta, X=dX) * dZ[i]
+ # cat("\np.b.j call\n")
+ tmp.j <- p.b.j(t=i*delta, X=dX)
+ #print(str(tmp.j))
+ #cat("\np.b.j cback and dZ\n")
+ # print(str(dZ[,i]))
+ # print(sum(dim(tmp.j)))
+ if(sum(dim(tmp.j))==2)
+ tmp.j <- as.numeric(tmp.j)
+ #print(str(tmp.j))
+ #print(str(p.b(t = i * delta, X = dX) %*% dW[, i]))
+ dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i]
X_mat[, i+1] <- dX
}
tsX <- ts(data=t(X_mat), deltat=delta, start=0)
Modified: pkg/yuima/R/simCP.R
===================================================================
--- pkg/yuima/R/simCP.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/simCP.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -7,27 +7,34 @@
modeltime <- sdeModel at time.variable
Terminal <- yuima at sampling@Terminal[1]
Initial <- yuima at sampling@Initial[1]
-
+ dimension <- yuima at model@dimension
+ dummy.val <- numeric(dimension)
+ if(length(xinit) != dimension)
+ xinit <- rep(xinit, dimension)[1:dimension]
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))
+ dummy.val<-rep(dummy.val,dimension)
}
for(i in 1:length(modelstate)){
assign(modelstate[i],dummy.val[i] ,envir=env)
}
- }
+ } else {
+ for(i in 1:dimension){
+ dummy.val[i] <- eval(xinit[i], envir=env)
+ }
+
+ }
-
### Simulation of CP using Lewis' method
##:: Levy
- JP <- eval(sdeModel at jump.coeff, envir=env)
+ JP <- eval(sdeModel at jump.coeff[[1]], envir=env)
mu.size <- length(JP)
-
+ # print(str(JP))
#assign(sdeModel at measure$intensity, env) ## intensity param
.CPintensity <- function(.t) {
@@ -48,12 +55,12 @@
# we use Lewis' acceptance/rejection method
- if(grep("^[dexp|dnorm|dgamma|dconst]", 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 dconst, dexp, dnorm and dgamma yet.")
- }
+ #} 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
@@ -75,10 +82,10 @@
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)
+ randJ <- as.matrix(randJ, ncol=yuima at dimension)
+
+ randJ <- rbind(dummy.val, randJ)
+ CP <- apply(randJ,2,cumsum)
tsX <- zoo(x=CP, order.by=E)
yuimaData <- setYuima(data=setData(tsX))
yuimaData <- subsampling(yuimaData, sampling=yuima at sampling)
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/R/yuima.model.R 2014-09-28 19:43:13 UTC (rev 333)
@@ -81,7 +81,8 @@
drift = expression() ,
diffusion = list() ,
hurst = 0.5,
- jump.coeff = expression(),
+ jump.coeff = list(),
+#jump.coeff = expression(),
measure=list(),
measure.type=character(),
parameter = new("model.parameter"),
@@ -161,23 +162,22 @@
# return(NULL)
# }
- }else if(measure.type=="CP"){ ##::CP
- if(length(measure)!=2){
- yuima.warn(paste("length of measure must be two on type", measure.type, "."))
- return(NULL)
- }
-
- if(!is.list(measure)){
- measure <- list(intensity=measure[1], df=measure[2])
- }else{
- if(length(measure[[1]])!=1 || length(measure[[2]])!=1){
- yuima.warn("multi dimentional jump term is not supported yet.")
- return(NULL)
- }
+ } else if(measure.type=="CP"){ ##::CP
+ # if(length(measure)!=2){
+ # yuima.warn(paste("length of measure must be two on type", measure.type, "."))
+ # return(NULL)
+ #}
+ if(!is.list(measure)){
+ measure <- list(intensity=measure[1], df=measure[2],dimension=measure[3])
+ } else {
+ #if(length(measure[[1]])!=1 || length(measure[[2]])!=1){
+ # yuima.warn("multi dimentional jump term is not supported yet.")
+ # return(NULL)
+ # }
##::naming measure list ########
tmpc <- names(measure)
if(is.null(tmpc)){
- names(measure) <- c("intensity", "df")
+ names(measure) <- c("intensity", "df","dimension")
}else{
whichint <- match("intensity", tmpc)
whichdf <- match("df", tmpc)
@@ -220,7 +220,8 @@
##measure.par$intensity <- unique(all.vars(MEASURE$intensity))
##::end check df name ####################
##::end CP
- }else if(measure.type=="code"){ ##::code
+
+ } else if(measure.type=="code"){ ##::code
if(length(measure)!=1){
yuima.warn(paste("length of measure must be one on type", measure.type, "."))
return(NULL)
@@ -308,7 +309,6 @@
##::end measure and jump term #####################################
-
##:: check for errors and reform values
if(any(time.variable %in% state.variable)){
yuima.warn("time and state(s) variable must be different.")
@@ -399,13 +399,8 @@
return(res)
}
- if (length(jump.coeff)==0){
- JUMP <- parse(text=jump.coeff)
- }else{
- JUMP <- vector(n.eqn, mode="expression")
- }
-
- ##:: make expressions of drifts and diffusions
+
+ ##:: make expressions of drifts and diffusions and jump
for(i in 1:n.eqn){
DRIFT[i] <- pre.proc(loc.drift[i,])
# 22/11 Simplify expressions
@@ -423,14 +418,57 @@
DIFFUSION[[i]][j] <- yuima.Simplifyobj(expr)
}
#22/11
- if (length(JUMP)>0){
- JUMP[i] <- parse(text=jump.coeff[i])
- JUMP[i] <- yuima.Simplifyobj(JUMP[i])
- }
+
+#if (length(JUMP)>0){
+# JUMP[i] <- parse(text=jump.coeff[i])
+# JUMP[i] <- yuima.Simplifyobj(JUMP[i])
+#}
+
}
+
+
+#print(length(jump.coeff))
+#if (length(jump.coeff)==0){
+# JUMP <- list(parse(text=jump.coeff))
+#}else{
+# # JUMP <- vector(n.eqn, mode="expression")
+# JUMP <- vector(n.eqn, mode="list")
+#}
+
+if(length(jump.coeff)==0){
+ JUMP <- list()
+} else {
+ if(length(jump.coeff)==1 & !is.matrix(jump.coeff)){ # is a scalar
+ expr <- parse(text=jump.coeff)
+ if(length(expr)==0){
+ expr <- expression(0) # expr must have something
+ }
+ JUMP <- list(yuima.Simplifyobj(expr))
+ } else { # must be matrix, n.col = dimension of Lévy noise
+ jump.coeff <- as.matrix(jump.coeff)
+ c.j <- ncol(jump.coeff)
+ r.j <- nrow(jump.coeff)
+ #print(c.j)
+ #print(r.j)
+ #print(jump.coeff)
+ JUMP <- vector(r.j, mode="list")
+ for(i in 1:r.j){
+ for(j in 1:c.j){
+ #cat(sprintf("\ni=%d,j=%d\n",i,j))
+ expr <- parse(text=jump.coeff[i,j])
+ if(length(expr)==0){
+ expr <- expression(0) # expr must have something
+ }
+ JUMP[[i]][j] <- yuima.Simplifyobj(expr)
+ }
+ }
+ }
+}
+#print(str(JUMP))
+
#
-
+
##:: get parameters in drift expression
drift.par <- unique(all.vars(DRIFT))
# Modification starting point 6/11
Modified: pkg/yuima/man/setPoisson.Rd
===================================================================
--- pkg/yuima/man/setPoisson.Rd 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/man/setPoisson.Rd 2014-09-28 19:43:13 UTC (rev 333)
@@ -15,13 +15,14 @@
}
\usage{
-setPoisson(intensity = 1, df = NULL, scale = 1, ...)
+setPoisson(intensity = 1, df = NULL, scale = 1, dimension=1, ...)
}
\arguments{
\item{intensity}{either and expression or a numerical value representing the intensity function of the Poisson process Nt.}
\item{df}{is the density of jump random variables Y.}
\item{scale}{this is the scaling factor \code{c}.}
+ \item{dimension}{this is the dimension of the jump component.}
\item{...}{passed to \code{\link{setModel}}}
}
Modified: pkg/yuima/man/simulate.Rd
===================================================================
--- pkg/yuima/man/simulate.Rd 2014-09-26 03:07:51 UTC (rev 332)
+++ pkg/yuima/man/simulate.Rd 2014-09-28 19:43:13 UTC (rev 333)
@@ -137,7 +137,8 @@
# A non-negative process (CIR process)
# dXt= a*(c-y)*dt + b*sqrt(Xt)*dWt
sq <- function(x){y = 0;if(x>0){y = sqrt(x);};return(y);}
- model<- setModel(drift="0.8*(0.2-x)",diffusion="0.5*sq(x)",solve.variable=c("x"))
+ model<- setModel(drift="0.8*(0.2-x)",
+ diffusion="0.5*sq(x)",solve.variable=c("x"))
T<-10
n<-1000
sampling <- setSampling(Terminal=T,n=n)
@@ -226,7 +227,8 @@
my.dW <- t(matrix(my.dW, nrow=n, ncol=yuima.obj at model@noise.number))
## Solve SDEs using Euler-Maruyama method.
-yuima.obj.path <- simulate(yuima.obj, space.discretized=FALSE, increment.W=my.dW)
+yuima.obj.path <- simulate(yuima.obj, space.discretized=FALSE,
+ increment.W=my.dW)
if( !is.null(yuima.obj.path) ){
dev.new()
# x11()
@@ -290,18 +292,24 @@
## dX = (theta - A X)dt + dZ,
## theta=(theta_1, theta_2) = c(1,.5)
## A=[a_ij], a_11 = 2, a_12 = 1, a_21 = 1, a_22=2
-## xinit <- c(1,1)
-## beta <- c(.1,.1)
-## mu <- c(0,0)
-## Lambda <- matrix(c(1,0,0,1),2,2)
-## obj.model <- setModel(drift=c("1 - 2*x1-x2",".5-x1-2*x2"), xinit=c(1,1),
-## solve.variable=c("x1","x2"), jump.coeff="1", measure.type="code",
-## measure=list(df="rIG(z, alpha=1, beta=beta, mu=mu, Lambda=Lambda)"))
-## obj.sampling <- setSampling(Terminal=10, n=10000)
-## obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
-## result <- simulate(obj.yuima)
-## plot(result)
+require(yuima)
+x0 <- c(1,1)
+beta <- c(.1,.1)
+mu <- c(0,0)
+delta0 <- 1
+alpha <- 1
+Lambda <- matrix(c(1,0,0,1),2,2)
+cc <- matrix(c(1,0,0,1),2,2)
+obj.model <- setModel(drift=c("1 - 2*x1-x2",".5-x1-2*x2"), xinit=x0,
+solve.variable=c("x1","x2"), jump.coeff=cc, measure.type="code",
+ measure=list(df="rNIG(z, alpha, beta, delta0, mu, Lambda)"))
+obj.sampling <- setSampling(Terminal=10, n=10000)
+obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
+result <- simulate(obj.yuima,true.par=list( alpha=alpha,
+ beta=beta, delta0=delta0, mu=mu, Lambda=Lambda))
+plot(result)
+
# Path-simulation for a Carma(p=2,q=1) model driven by a Brownian motion:
carma1<-setCarma(p=2,q=1)
str(carma1)
More information about the Yuima-commits
mailing list