[Yuima-commits] r260 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 20 12:50:04 CET 2013
Author: iacus
Date: 2013-11-20 12:50:03 +0100 (Wed, 20 Nov 2013)
New Revision: 260
Added:
pkg/yuima/R/setCarma.R
pkg/yuima/man/setCarma.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/adaBayes.R
pkg/yuima/R/asymptotic_term_second.R
pkg/yuima/R/asymptotic_term_third.R
pkg/yuima/R/poisson.random.sampling.R
pkg/yuima/R/sim.euler.R
pkg/yuima/R/sim.euler.space.discretized.R
pkg/yuima/R/simulate.R
pkg/yuima/R/toLatex.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/model.parameter-class.Rd
Log:
added setCarma and modifications to yuima class
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/DESCRIPTION 2013-11-20 11:50:03 UTC (rev 260)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.214
-Date: 2013-10-30
+Version: 0.1.215
+Date: 2013-11-20
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/NAMESPACE 2013-11-20 11:50:03 UTC (rev 260)
@@ -58,6 +58,7 @@
export(setData)
export(setSampling)
export(setCharacteristic)
+export(setCarma)
export(dim)
export(length)
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/AllClasses.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -1,98 +1,101 @@
-# Class Definitions
-# This source MUST be loaded first
-
-# Class 'yuima.pars'
-
-# parameter object included in 'yuima.model'
-setClass("model.parameter",representation(all="character",
- common="character",
- diffusion="character",
- drift="character",
- jump="character",
- measure="character"
- )
- )
-
-# Class 'yuima.model'
-setClass("yuima.model",representation(drift="expression",
- diffusion="list",
- hurst="ANY",
- jump.coeff="expression",
- measure="list",
- measure.type="character",
- parameter="model.parameter",
- state.variable="character",
- jump.variable="character",
- time.variable="character",
- noise.number="numeric",
- equation.number="numeric",
- dimension="numeric",
- solve.variable="character",
- xinit="numeric",
- J.flag="logical"
- )
- )
-
-# Class 'yuima.data'
-
-# we want yuimaS4 to use any class of data as input
-# the original data will be stored in OrigData
-# we convert these objects internally to "zoo" object
-# in the future, we may want to use more flexible
-# classes
-
-setClass("yuima.data", representation(original.data = "ANY",
- zoo.data = "ANY"
- )
- )
-
-
-# Class 'yuima.sampling'
-
-# sampling is now empty, but should give informations on the sampling
-# type, rate, deltas, etc.
-
-setClass("yuima.sampling", representation(Initial = "numeric",
- Terminal = "numeric",
- n = "numeric",
- delta = "numeric",
- grid = "ANY",
- random = "ANY",
- regular = "logical",
- sdelta = "numeric",
- sgrid = "ANY",
- oindex = "ANY",
- interpolation = "character"
- )
- )
-
-# Class 'yuima.functional'
-
-# functional model used in 'asymptotic term' procedure
-
-setClass("yuima.functional", representation(F = "ANY",
- f = "list",
- xinit = "numeric",
- e = "numeric"
- )
- )
-
-
-# Class 'yuima'
-
-# this is the principal class of yuima project. It may contain up to
-# three slots for now: the data, the model and the sampling
-
-setClass("yuima.characteristic", representation(equation.number = "numeric",
- time.scale = "numeric"
- )
- )
-
-
-setClass("yuima", representation(data = "yuima.data",
- model = "yuima.model",
- sampling = "yuima.sampling",
- characteristic = "yuima.characteristic",
- functional = "yuima.functional"
- )
- )
+# Class Definitions
+# This source MUST be loaded first
+
+# Class 'yuima.pars'
+
+# parameter object included in 'yuima.model'
+setClass("model.parameter",representation(all="character",
+ common="character",
+ diffusion="character",
+ drift="character",
+ jump="character",
+ measure="character",
+# Insert parameters for starting conditions
+ xinit="character"
+ )
+ )
+
+# Class 'yuima.model'
+setClass("yuima.model",representation(drift="expression",
+ diffusion="list",
+ hurst="ANY",
+ jump.coeff="expression",
+ measure="list",
+ measure.type="character",
+ parameter="model.parameter",
+ state.variable="character",
+ jump.variable="character",
+ time.variable="character",
+ noise.number="numeric",
+ equation.number="numeric",
+ dimension="numeric",
+ solve.variable="character",
+# xinit="numeric",
+ xinit="expression",
+ J.flag="logical"
+ )
+ )
+
+# Class 'yuima.data'
+
+# we want yuimaS4 to use any class of data as input
+# the original data will be stored in OrigData
+# we convert these objects internally to "zoo" object
+# in the future, we may want to use more flexible
+# classes
+
+setClass("yuima.data", representation(original.data = "ANY",
+ zoo.data = "ANY"
+ )
+ )
+
+
+# Class 'yuima.sampling'
+
+# sampling is now empty, but should give informations on the sampling
+# type, rate, deltas, etc.
+
+setClass("yuima.sampling", representation(Initial = "numeric",
+ Terminal = "numeric",
+ n = "numeric",
+ delta = "numeric",
+ grid = "ANY",
+ random = "ANY",
+ regular = "logical",
+ sdelta = "numeric",
+ sgrid = "ANY",
+ oindex = "ANY",
+ interpolation = "character"
+ )
+ )
+
+# Class 'yuima.functional'
+
+# functional model used in 'asymptotic term' procedure
+
+setClass("yuima.functional", representation(F = "ANY",
+ f = "list",
+ xinit = "numeric",
+ e = "numeric"
+ )
+ )
+
+
+# Class 'yuima'
+
+# this is the principal class of yuima project. It may contain up to
+# three slots for now: the data, the model and the sampling
+
+setClass("yuima.characteristic", representation(equation.number = "numeric",
+ time.scale = "numeric"
+ )
+ )
+
+
+setClass("yuima", representation(data = "yuima.data",
+ model = "yuima.model",
+ sampling = "yuima.sampling",
+ characteristic = "yuima.characteristic",
+ functional = "yuima.functional"
+ )
+ )
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/adaBayes.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -144,13 +144,13 @@
n <- length(yuima)[1]
env <- new.env()
- assign("X", as.matrix(onezoo(yuima)), env=env)
- assign("deltaX", matrix(0, n-1, d.size), env=env)
- assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), env=env)
+ assign("X", as.matrix(onezoo(yuima)), envir=env)
+ assign("deltaX", matrix(0, n-1, d.size), envir=env)
+ assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
for(t in 1:(n-1))
env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
- assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
+ assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
mle <- qmle(yuima, "start"=start, "lower"=lower,"upper"=upper, "method"="L-BFGS-B")
integ <- function(idx.fixed=NULL,f=f,start=start,par=NULL,hessian=FALSE,upper,lower){
Modified: pkg/yuima/R/asymptotic_term_second.R
===================================================================
--- pkg/yuima/R/asymptotic_term_second.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/asymptotic_term_second.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -436,7 +436,7 @@
}
if(require(cubature)){
- tmp <- adaptIntegrate(gz_pi0, lower=lower,upper=upper)$integral
+ tmp <- adaptIntegrate(gz_pi0, lowerLimit=lower,upperLimit=upper)$integral
} else {
tmp <- NA
}
Modified: pkg/yuima/R/asymptotic_term_third.R
===================================================================
--- pkg/yuima/R/asymptotic_term_third.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/asymptotic_term_third.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -434,7 +434,7 @@
}
if(require(cubature)){
- tmp <- adaptIntegrate(gz_pi0, lower=lower,upper=upper)$integral
+ tmp <- adaptIntegrate(gz_pi0, lowerLimit=lower,upperLimit=upper)$integral
} else {
tmp <- NA
}
Modified: pkg/yuima/R/poisson.random.sampling.R
===================================================================
--- pkg/yuima/R/poisson.random.sampling.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/poisson.random.sampling.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -22,11 +22,11 @@
Elength <- T*n*rate[i]
## make random numbers following exponential distribution
- Time <- diffinv(rexp(Elength,r=deltat*n*rate[i]))
+ Time <- diffinv(rexp(Elength,rate=deltat*n*rate[i]))
# make up a deficit
while(Time[length(Time)]< Num){
# adding random numbers (by 10% of Elength)
- Time2 <- diffinv(rexp(trunc(Elength/10+1),r=deltat*n*rate[i]))+Time[length(Time)]
+ Time2 <- diffinv(rexp(trunc(Elength/10+1),rate=deltat*n*rate[i]))+Time[length(Time)]
# restrain duplication
Time <- append(Time,Time2[-1])
}
Added: pkg/yuima/R/setCarma.R
===================================================================
--- pkg/yuima/R/setCarma.R (rev 0)
+++ pkg/yuima/R/setCarma.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -0,0 +1,291 @@
+# Carma_Model<-setClass("Carma_Model",
+# slots = c(Cogarch_Model_Log="logical",
+# Under_Lev="yuima.model"),
+# prototype=list(Cogarch_Model_Log = FALSE,
+# Under_Lev=NULL),
+# contains= "yuima")
+#ma.par, ar.par, lin.par=NULL is.null
+#state Variable consistente con Yuima.????
+# Inserire Solve Variable???
+# ... che mi serve per passare gli stessi parametri di setModel
+# per i ... guardare qmle
+# call<-matchcall()
+# mydots guardare
+setCarma<-function(p,q,loc.par=NULL,ar.par="beta",ma.par="alpha",lin.par=NULL,Carma.var="v",Latent.var="x", ...){
+# We use the same parametrization as in Brockwell[2000]
+
+# mydots$Carma.var= V
+# mydots$Latent.var= X ?????
+ # The Carma model is defined by two equations:
+ # \begin{eqnarray}
+ # V_{t}&=&\alpha_{0}+a'Y_{t-}\\
+ # dY_{t}&=& BY_{t-}dt+ e dL\\
+ # \end{eqnarray}
+ # p is the number of the moving average parameters \alpha
+ # q is the number of the autoregressive coefficient \beta
+
+#Default parameters
+
+ call <- match.call()
+ quadratic_variation<-FALSE
+
+ mydots <- as.list(call)[-(1:2)]
+
+
+# hurst<-0.5
+# jump.coeff <- character()
+# measure <- list()
+# measure.type <- character()
+# state.variable <- "Y"
+# jump.variable <- "z"
+# time.variable <- "t"
+# mydots$xinit<- NULL
+ if (is.null(mydots$hurst)){
+ mydots$hurst<-0.5
+ }
+
+ if(is.null(mydots$time.variable)){
+ mydots$time.variable<-"t"
+ }
+
+ if(is.null(mydots$jump.variable)){
+ mydots$jump.variable<-"z"
+ }
+
+# if(is.null(mydots$Carma.var)){
+# Carma.var<-"V"
+# } else{
+# Carma.var<-mydots$Carma.var
+# }
+
+# if(is.null(mydots$Latent.var)){
+# Latent.var<-"X"
+# } else{
+# Latent.var<-mydots$Latent.var
+# }
+
+ if(is.null(mydots$xinit)){
+ if(is.null(mydots$XinExpr)){
+ mydots$xinit<-as.character(0*c(1:p))
+ }else{
+ if(mydots$XinExpr==TRUE){
+ Int.Var<-paste(Latent.var,"0",sep="")
+ mydots$xinit<-paste(Int.Var,c(0:(p-1)),sep="")
+ }
+ }
+ } else{
+ dummy<-as.character(mydots$xinit)
+ mydots$xinit<-dummy[-1]
+ }
+
+ if(p<q){
+ yuima.stop("order of AR must be larger than MA order")
+ }
+
+ beta_coeff0<-paste("-",ar.par,sep="")
+ beta_coeff<-paste(beta_coeff0,p:1,sep="")
+ coeff_alpha<-c(paste(ma.par,0:q,sep=""),as.character(matrix(0,1,p-q-1)))
+ fin_alp<-length(coeff_alpha)
+ # We built the drift condition
+
+ Y_coeff<-paste(Latent.var,0:(p-1),sep="")
+ fin_Y<-length(Y_coeff)
+ V1<-paste(coeff_alpha,Y_coeff,sep="*")
+ V2<-paste(V1,collapse="+")
+# alpha0<-paste(ma.par,0,sep="")
+ if(is.null(loc.par)){
+ V<-paste("(",V2,")",collapse="")
+ } else {
+ Vt<-paste(loc.par,V2,sep="+")
+ V<-paste("(",Vt,")",collapse="")
+ }
+ drift_last_cond<-paste(paste(beta_coeff,Y_coeff,sep="*"),collapse="")
+ # Drift condition for the dV_{t}
+
+ drift_first_cond_1<-c(paste(coeff_alpha[-fin_alp],Y_coeff[-1],sep="*"))
+ drift_first_cond_2<-paste(drift_first_cond_1,collapse="+")
+ drift_first_cond_a<-paste("(",drift_last_cond,")",sep="")
+ drift_first_cond_b<-paste(coeff_alpha[fin_alp],drift_first_cond_a,sep="*")
+ drift_first_cond<-paste(drift_first_cond_2,drift_first_cond_b,sep="+")
+
+ if(length(Y_coeff)>1)
+ {drift_Carma<-c(drift_first_cond,Y_coeff[2:length(Y_coeff)],drift_last_cond)}else
+ {drift_Carma<-c(drift_first_cond,drift_last_cond)}
+ # We need to consider three different situations
+
+ if(is.null(mydots$jump.coeff) & is.null(mydots$measure) &
+ is.null(mydots$measure.type) & quadratic_variation==FALSE){
+ # The Carma model is driven by a Brwonian motion
+ if (is.null(lin.par)){
+ diffusion_Carma<-matrix(c(coeff_alpha[fin_alp],as.character(matrix(0,(p-1),1)),"1"),(p+1),1)
+ # Latent.var<-Y_coeff
+ Model_Carma1<-setModel(drift=drift_Carma,
+ diffusion=diffusion_Carma,
+ hurst=mydots$hurst,
+ state.variable=c(Carma.var,Y_coeff),
+ solve.variable=c(Carma.var,Y_coeff),
+ xinit=c(V,mydots$xinit))
+ if(length(Model_Carma1)==0){
+ stop("Yuima model was not built")
+ } else {
+ return(Model_Carma1)
+ }
+ } else{
+ if(ma.par==lin.par){
+ first_term<-paste(coeff_alpha[fin_alp],V,sep="*")
+ diffusion_Carma<-matrix(c(first_term,as.character(matrix(0,(p-1),1)),V),(p+1),1)
+ Model_Carma1<-setModel(drift=drift_Carma,
+ diffusion=diffusion_Carma,
+ hurst=mydots$hurst,
+ state.variable=c(Carma.var,Y_coeff),
+ solve.variable=c(Carma.var,Y_coeff),
+ xinit=c(V,mydots$xinit))
+ return(Model_Carma1)
+ }else{
+# coeff_gamma<-c(paste(lin.par,1:p,sep=""),as.character(matrix(0,1,p-q)))
+ coeff_gamma<-c(paste(lin.par,1:p,sep=""))
+ Gamma1<-paste(coeff_gamma,Y_coeff,sep="*")
+ Gamma2<-paste(Gamma1,collapse="+")
+ gamma0<-paste(lin.par,0,sep="")
+ Gammat<-paste(gamma0,Gamma2,sep="+")
+ Gamma<-paste("(",Gammat,")",collapse="")
+ first_term<-paste(coeff_alpha[fin_alp],Gamma,sep="*")
+
+ diffusion_Carma<-matrix(c(first_term,as.character(matrix(0,(p-1),1)),Gamma),(p+1),1)
+ Model_Carma1<-setModel(drift=drift_Carma,
+ diffusion=diffusion_Carma,
+ hurst=mydots$hurst,
+ state.variable=c(Carma.var,Y_coeff),
+ solve.variable=c(Carma.var,Y_coeff),
+ xinit=c(V,mydots$xinit))
+
+ return(Model_Carma1)
+ }
+
+ }
+ } else {
+ if( is.null(mydots$jump.coeff) & is.null(mydots$measure) &
+ is.null(mydots$measure.type) & is.null(lin.par) &
+ quadratic_variation==TRUE){
+
+ stop("The CoGarch model requires a Carma process driven by the discrete part of the quadratic covariation:
+ You Must specify the Levy Measure")
+
+ } else {
+
+ if(quadratic_variation==FALSE & is.null(lin.par)){
+ # warning("At the moment, we need specify the underlying L?vy directly")
+ # diffusion_Carma<-matrix(c(coeff_alpha[fin_alp],as.character(matrix(0,(q-1),1)),"1"),(q+1),1)
+ # Model_Carma1<-setModel(drift=drift_Carma, diffusion=diffusion_Carma,
+ # hurst=hurst, state.variable=c(V,Y_coeff), solve.variable=c(V,Y_coeff))
+ # under_Lev1<-setModel(drift="0",diffusion="0",jump.coeff="1" ,
+ # measure=measure ,measure.type=measure.type ,
+ # jump.variable=jump.variable , time.variable=time.variable)
+ # if(length(Model_Carma1)==0){
+ # stop("Yuima model was not built")
+ # } else {
+ # Model_Carma<-Carma_Model()
+ # Model_Carma at model <- Model_Carma1
+ # Model_Carma at Cogarch_Model_Log <- Cogarch_Model
+ # Model_Carma at Under_Lev <-under_Lev1
+ # return(Model_Carma)
+ # }
+
+ # LM 27/09 We use a modified
+ # setModel that allows us to build a sde where the slot model at jump.coeff is an vector
+
+ # jump_Carma<-matrix(c(coeff_alpha[fin_alp],as.character(matrix(0,(q-1),1)),"1"),(q+1),1)
+ jump_Carma<-c(coeff_alpha[fin_alp],as.character(matrix(0,(p-1),1)),"1")
+ Model_Carma<-setModel(drift=drift_Carma,
+ diffusion = NULL,
+ hurst=mydots$hurst,
+ jump.coeff=jump_Carma,
+ measure=eval(mydots$measure),
+ measure.type=mydots$measure.type,
+ jump.variable=mydots$jump.variable,
+ time.variable=mydots$time.variable,
+ state.variable=c(Carma.var,Y_coeff),
+ solve.variable=c(Carma.var,Y_coeff),
+ xinit=c(V,mydots$xinit))
+ return(Model_Carma)
+ } else {
+ if (quadratic_variation==FALSE ){
+ # Selecting Quadratic_Variation==FALSE and specifying the Heteroskedatic.param in the model,
+ # The coefficient of the error term is a vector in which the last element is an affine linear function
+ # of the vector space Y
+
+ # We have to consider two different cases:
+ # a) The last component of the error term is $V_{t-}=\alpha_{0}+a'Y_{t-}$. Usually
+ # this condition is used for the definition of the COGARCH(p,q) introduced in Brockwell and Davis and
+ # b) The last component of the error term is a linear function of the state variable $Y_{t}$
+ # different of the observable variable V.
+ if(ma.par==lin.par){
+ jump_first_comp<-paste(coeff_alpha[fin_alp],V,sep="*")
+ jump_Carma<-c(jump_first_comp,as.character(matrix(0,(p-1),1)),V)
+ }else{
+# coeff_gamma<-c(paste(lin.par,1:p,sep=""),as.character(matrix(0,1,q-p)))
+ coeff_gamma<-c(paste(lin.par,1:p,sep=""))
+ Gamma1<-paste(coeff_gamma,Y_coeff,sep="*")
+ Gamma2<-paste(Gamma1,collapse="+")
+ gamma0<-paste(lin.par,0,sep="")
+ Gammat<-paste(gamma0,Gamma2,sep="+")
+ Gamma<-paste("(",Gammat,")",collapse="")
+ jump_first_comp<-paste(coeff_alpha[fin_alp],Gamma,sep="*")
+ jump_Carma<-c(jump_first_comp,as.character(matrix(0,(p-1),1)),Gamma)
+ }
+
+# Model_Carma<-setModel(drift=drift_Carma,
+# diffusion =NULL,
+# hurst=0.5,
+# jump.coeff=jump_Carma,
+# measure=eval(mydots$measure),
+# measure.type=mydots$measure.type,
+# jump.variable=mydots$jump.variable,
+# time.variable=mydots$time.variable,
+# state.variable=c(Carma.var,Y_coeff),
+# solve.variable=c(Carma.var,Y_coeff),
+# c(V,mydots$xinit))
+# return(Model_Carma)
+ }
+
+ Model_Carma<-setModel(drift=drift_Carma,
+ diffusion =NULL,
+ hurst=mydots$hurst,
+ jump.coeff=jump_Carma,
+ measure=eval(mydots$measure),
+ measure.type=mydots$measure.type,
+ jump.variable=mydots$jump.variable,
+ time.variable=mydots$time.variable,
+ state.variable=c(Carma.var,Y_coeff),
+ solve.variable=c(Carma.var,Y_coeff),
+ c(V,mydots$xinit))
+ return(Model_Carma)
+ if(quadratic_variation==TRUE){
+#
+ stop("Work in Progress: Implementation of CARMA model for CoGarch.
+ We need the increments of Quadratic Variation")
+#
+# diffusion_first_cond<-paste(coeff_alpha[fin_alp],V,sep="*")
+# diffusion_Carma<-matrix(c(diffusion_first_cond,as.character(matrix(0,(q-1),1)),Vt),(q+1),1)
+# # At the present time, Yuima does not support Multi - Jumps
+# Model_Carma1<-setModel(drift=drift_Carma, diffusion=diffusion_Carma,
+# hurst=hurst, state.variable=c(V,Y_coeff), solve.variable=c(V,Y_coeff))
+# under_Lev1<-setModel(drift="0",diffusion="0",jump.coeff="1" ,
+# measure=measure ,measure.type=measure.type ,
+# jump.variable=jump.variable , time.variable=time.variable)
+# if(length(Model_Carma1)==0){
+# stop("Yuima model was not built")
+# } else {
+# Model_Carma<-Carma_Model()
+# Model_Carma at model <- Model_Carma1
+# Model_Carma at Cogarch_Model_Log <- Cogarch_Model
+# Model_Carma at Under_Lev <-under_Lev1
+# return(Model_Carma)
+# }
+#
+ }
+ }
+ }
+
+ }
+}
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/sim.euler.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -12,8 +12,30 @@
Terminal <- yuima at sampling@Terminal[1]
n <- yuima at sampling@n[1]
- dX <- xinit
+# dX <- xinit
+
+ # 06/11 xinit is an expression: the structure is equal to that of V0
+ 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)
+ }
+ # 06/11 we need a initial variable for X_0
+ dX<-vector(mode="numeric",length(dX_dummy))
+ for(i in 1:length(dX_dummy)){
+ dX[i] <- eval(dX_dummy[i], env)
+ }
+
+
+
##:: set time step
delta <- Terminal/n
Modified: pkg/yuima/R/sim.euler.space.discretized.R
===================================================================
--- pkg/yuima/R/sim.euler.space.discretized.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/sim.euler.space.discretized.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -1,106 +1,126 @@
-space.discretized<-function(xinit,yuima, env){
-
-
-##:: initialize state variable
- 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]
- n <- yuima at sampling@n[1]
- dX <- xinit
-
-##:: set time step
- delta <- Terminal/n
-
-
-
-
-##:: using Space-discretized Euler-Maruyama method
-
-
-
-##:: function for approximation of function G
-gfunc <- function(x){
- c0 <- 10
- c1 <- 10
- ret <- rep(0, length(x))
- idx <- which(x < 1/c0)
- ret[idx] <- 1
-
- idx <- which(1/c0 <= x)
- ret[idx] <- 1-pnorm(x[idx])
- for(i in 1:length(idx)){
- n <- 1:floor(c1/x[idx[i]])
- ret[idx[i]] <- 4 * (ret[idx[i]] - sum( pnorm((4*n+1)*x[idx[i]]) - pnorm((4*n-1)*x[idx[i]]) ))
- }
-
- idx <- which(1 < ret)
- ret[idx] <- 1
- return(ret)
-}
-
-
-
- dxx <- 0.0001
- xx <- seq(0, 1.7, dxx)
-
- ##:: approximate function G(gg)
- gg <- gfunc(xx)
- appfunc <- suppressWarnings(approxfun(gg, xx))
-
- ##:: calculate inverse of G
- unif.a <- runif(n*2)
- inv.a <- pmin(qnorm(1 - unif.a/4), appfunc(unif.a), na.rm=TRUE)
-
- ##:: make random time steps
- ep <- sqrt(delta)
- dTW <- (ep/inv.a)^2
- time_idx <- cumsum(dTW) ##:: time index should be attached
- div_sd <- min(which(time_idx > Terminal)) ##:: cut by time=1
- time_idx <- time_idx[1:div_sd]
-
- ##:: add diffusion term
- dTW <- rbind(dTW[1:div_sd],
- t(matrix( (rbinom(div_sd*r.size, 1, 0.5)*2-1) * ep,
- nrow=div_sd,
- ncol=r.size)
- )
- )
-
- X_mat <- matrix(0, d.size, div_sd+1)
- X_mat[,1] <- dX
-
- ##:: function to calculate coefficients of dTW
- p.b <- function(t, X=numeric(d.size)){
- ##:: assign names of variables
- for(i in 1:length(modelstate)){
- assign(modelstate[i], X[i], env)
- }
- assign(modeltime, t, env)
- 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)
- }
- }
- return(tmp)
- }
- ##:: calcurate difference equation
- for(i in 1:div_sd){
- dX <- dX + p.b(t=time_idx[i], X=dX) %*% dTW[,i]
- X_mat[,i+1] <- dX
- }
- ##tsX <- ts(data=t(X_mat), deltat=delta , start=0)
- ##:: output zoo data
- zooX <- zoo(x=t(X_mat), order.by=c(0, time_idx))
- yuimaData <- setData(original.data=zooX)
- return(yuimaData)
-
-}
+space.discretized<-function(xinit,yuima, env){
+
+
+##:: initialize state variable
+ 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]
+ n <- yuima at sampling@n[1]
+# dX <- xinit
+
+ 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)
+ }
+ # 20/11 we need a initial variable for X_0
+ dX<-vector(mode="numeric",length(dX_dummy))
+
+ for(i in 1:length(dX_dummy)){
+ dX[i] <- eval(dX_dummy[i], env)
+ }
+
+
+##:: set time step
+ delta <- Terminal/n
+
+
+
+
+##:: using Space-discretized Euler-Maruyama method
+
+
+
+##:: function for approximation of function G
+gfunc <- function(x){
+ c0 <- 10
+ c1 <- 10
+ ret <- rep(0, length(x))
+ idx <- which(x < 1/c0)
+ ret[idx] <- 1
+
+ idx <- which(1/c0 <= x)
+ ret[idx] <- 1-pnorm(x[idx])
+ for(i in 1:length(idx)){
+ n <- 1:floor(c1/x[idx[i]])
+ ret[idx[i]] <- 4 * (ret[idx[i]] - sum( pnorm((4*n+1)*x[idx[i]]) - pnorm((4*n-1)*x[idx[i]]) ))
+ }
+
+ idx <- which(1 < ret)
+ ret[idx] <- 1
+ return(ret)
+}
+
+
+
+ dxx <- 0.0001
+ xx <- seq(0, 1.7, dxx)
+
+ ##:: approximate function G(gg)
+ gg <- gfunc(xx)
+ appfunc <- suppressWarnings(approxfun(gg, xx))
+
+ ##:: calculate inverse of G
+ unif.a <- runif(n*2)
+ inv.a <- pmin(qnorm(1 - unif.a/4), appfunc(unif.a), na.rm=TRUE)
+
+ ##:: make random time steps
+ ep <- sqrt(delta)
+ dTW <- (ep/inv.a)^2
+ time_idx <- cumsum(dTW) ##:: time index should be attached
+ div_sd <- min(which(time_idx > Terminal)) ##:: cut by time=1
+ time_idx <- time_idx[1:div_sd]
+
+ ##:: add diffusion term
+ dTW <- rbind(dTW[1:div_sd],
+ t(matrix( (rbinom(div_sd*r.size, 1, 0.5)*2-1) * ep,
+ nrow=div_sd,
+ ncol=r.size)
+ )
+ )
+
+ X_mat <- matrix(0, d.size, div_sd+1)
+ X_mat[,1] <- dX
+
+ ##:: function to calculate coefficients of dTW
+ p.b <- function(t, X=numeric(d.size)){
+ ##:: assign names of variables
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], X[i], env)
+ }
+ assign(modeltime, t, env)
+ 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)
+ }
+ }
+ return(tmp)
+ }
+ ##:: calcurate difference equation
+ for(i in 1:div_sd){
+ dX <- dX + p.b(t=time_idx[i], X=dX) %*% dTW[,i]
+ X_mat[,i+1] <- dX
+ }
+ ##tsX <- ts(data=t(X_mat), deltat=delta , start=0)
+ ##:: output zoo data
+ zooX <- zoo(x=t(X_mat), order.by=c(0, time_idx))
+ yuimaData <- setData(original.data=zooX)
+ return(yuimaData)
+
+}
\ No newline at end of file
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2013-10-30 03:55:05 UTC (rev 259)
+++ pkg/yuima/R/simulate.R 2013-11-20 11:50:03 UTC (rev 260)
@@ -1,238 +1,240 @@
-## We have splitted the simulate function into blocks to allow for future
-## methods to be added. S.M.I. & A.B.
-## Interface to simulate() changed to match the S3 generic function in the
-## package stats
-## added an environment to let R find the proper values defined in the main
-## body of the function, which in turn calls different simulation methods
-## All new simulation methods should look into the yuimaEnv for local variables
-## when they need to "eval" R expressions
-
-##:: function simulate
-##:: solves SDE and returns result
-subsampling <- function(x,y) return(x)
-
-setGeneric("simulate",
- function(object, nsim=1, seed=NULL, xinit, true.parameter, space.discretized=FALSE,
- increment.W=NULL, increment.L=NULL, hurst, methodfGn="WoodChan",
- sampling=sampling, subsampling=subsampling, ...
-# Initial = 0, Terminal = 1, n = 100, delta,
-# grid = as.numeric(NULL), random = FALSE, sdelta=as.numeric(NULL),
-# sgrid=as.numeric(NULL), interpolation="none"
- )
- standardGeneric("simulate")
- )
-
-
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 260
More information about the Yuima-commits
mailing list