[Yuima-commits] r33 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 25 15:13:29 CET 2009
Author: abrouste
Date: 2009-11-25 15:13:28 +0100 (Wed, 25 Nov 2009)
New Revision: 33
Modified:
pkg/yuima/R/AllClasses.R
pkg/yuima/R/CholeskyfGn.R
pkg/yuima/R/sampling2grid.R
pkg/yuima/R/simulate.R
pkg/yuima/R/yuima.model.R
pkg/yuima/R/yuima.sampling.R
Log:
yuima at sampling class up + simulate function with fgn methods
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/AllClasses.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -51,11 +51,13 @@
# sampling is now empty, but should give informations on the sampling
# type, rate, deltas, etc.
-setClass("yuima.sampling", representation(Terminal = "numeric",
+setClass("yuima.sampling", representation(Initial = "numeric",
+ Terminal = "numeric",
division = "numeric",
- Initial = "numeric",
+ delta = "numeric",
grid = "numeric",
- random = "logical"
+ random = "ANY",
+ regular = "logical"
)
)
Modified: pkg/yuima/R/CholeskyfGn.R
===================================================================
--- pkg/yuima/R/CholeskyfGn.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/CholeskyfGn.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -1,5 +1,5 @@
CholeskyfGn <-
-function(mesh, H)
+function(mesh, H, dim)
{
##--------------------------------------------------------------------------
@@ -9,7 +9,8 @@
##--------------------------------------------------------------------------
## Input : mesh : mesh grid where the fBm is evaluated
-## H : self-similarity parameter
+## H : self-similarity parameter
+## dim : valued in R^dim
##
##
## Output : simulation of a standard fractional Brownian noise
@@ -22,10 +23,12 @@
## Taille memoire N^2
## -------------------------------------------------------------------------
- N<-length(mesh)-1 #j'ai retirer 1 # N+1 is the size of the fGn sample to be simulated
+N<-length(mesh)-2 # N+1 is the size of the fGn sample to be simulated
+fGn<-matrix(0,dim,N+1)
+ for (k in 1:dim){
matcov <- matrix(0,N+1,N+1) # Covariance matrix of the fGn
H2 <- 2*H
-
+
for (i in (1:(N+1))) {
j <- i:(N+1)
matcov[i, j]<- 0.5 * (abs(mesh[i+1]-mesh[j])^H2 + abs(mesh[i]-mesh[j+1])^H2 - abs(mesh[i] - mesh[j])^H2-abs(mesh[i+1] - mesh[j+1])^H2)
@@ -33,6 +36,8 @@
}
L <- chol(matcov)
Z <- rnorm(N+1)
- fGn <- t(L) %*% Z
+ fGn[k,] <- t(L) %*% Z
+ }
+
return(fGn)
}
\ No newline at end of file
Modified: pkg/yuima/R/sampling2grid.R
===================================================================
--- pkg/yuima/R/sampling2grid.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/sampling2grid.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -10,8 +10,8 @@
#Check if this grid is random
if (samp at random==FALSE){
-#Regular deterministic grid
- return(seq(samp at Initial,samp at Terminal,length=samp at division))
+#Regular deterministic grid
+ return(seq(samp at Initial[1],samp at Terminal[1],length=samp at division[1]+1)) #Attention! do not treat multifrequency grid
}
return(1) #To do: Random grid
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/simulate.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -10,12 +10,12 @@
##:: function simulate
##:: solves SDE and returns result
setGeneric("simulate",
- function(object, nsim, seed, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL)
+ function(object, nsim, seed, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL, methodfGn="Cholesky")
standardGeneric("simulate")
)
setMethod("simulate", "yuima",
- function(object, nsim=1, seed=NULL, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL){
+ function(object, nsim=1, seed=NULL, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL,methodfGn="Cholesky"){
##:: errors checks
@@ -130,12 +130,18 @@
if(missing(increment.W)){
if( sdeModel at hurst!=0.5 ){
grid<-sampling2grid(yuima at sampling)
- dW<-CholeskyfGn(grid, sdeModel at hurst)
- dW <- t(matrix(dW, nrow=division, ncol=r.size))
+
+ if(methodfGn=="Cholesky"){
+ dW<-CholeskyfGn(grid, sdeModel at hurst,r.size)
+ }else{
+ cat("\nNot done presently\n")
+ return(NULL)
+ }
+
} else {
delta<-Terminal/division
dW <- rnorm(division * r.size, 0, sqrt(delta))
- dW <- t(matrix(dW, nrow=division, ncol=r.size))
+ dW <- matrix(dW, ncol=division, nrow=r.size,byrow=TRUE)
}
} else {
dW <- increment.W
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/yuima.model.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -55,7 +55,7 @@
## set yuima model from SDE
setModel <- function(drift=NULL,
diffusion=NULL,
-hurst=0.5, # as.numeric(NULL),
+ hurst=0.5,
jump.coeff=character(),
measure=list(),
measure.type=character(),
Modified: pkg/yuima/R/yuima.sampling.R
===================================================================
--- pkg/yuima/R/yuima.sampling.R 2009-11-24 22:23:48 UTC (rev 32)
+++ pkg/yuima/R/yuima.sampling.R 2009-11-25 14:13:28 UTC (rev 33)
@@ -5,35 +5,77 @@
setMethod("initialize", "yuima.sampling",
- function(.Object, Terminal, division, Initial, grid, random){
- eqn <- length(Terminal)
+ function(.Object, Initial, Terminal, division, delta, grid, random, regular){
+
+ if(length(grid)>0){
+
+ testInitial<-(min(grid)==Initial)
+ testTerminal<-(max(grid)==Terminal)
+ testdivision<-(abs(division-diff(range(grid))/mean(diff(grid))+1)<10^(-10))
+ testdelta<-(abs(delta-mean(diff(grid)))<10^(-10))
+ testregular<-all(abs(diff(diff(grid)))<10^(-10))
+
+ if(!testInitial){
+ cat("\n Start time has been set with the grid \n")
+ }
+ if(!testTerminal){
+ cat("\n Terminal time has been set with the grid \n")
+ }
+ if(!testdivision){
+ cat("\n Division has been set with the grid \n")
+ }
+ if(!testdelta){
+ cat("\n delta has been set with the grid \n")
+ }
+ if(testregular){
+ .Object at division <- diff(range(grid))/mean(diff(grid))+1
+ .Object at delta <- mean(diff(grid))
+ .Object at regular <- TRUE
+ }else{
+ .Object at division <- length(grid)-1
+ .Object at delta <- as.numeric(NULL)
+ .Object at regular <- FALSE
+ }
+ .Object at grid <- grid
+ .Object at Initial <- min(grid)
+ .Object at Terminal <- max(grid)
+ .Object at random <- random
+ }else{
+ # There is no grid
+ eqn <- length(Terminal)
if(length(Terminal)==length(division)){
- .Object at Terminal <- Terminal
+ .Object at Initial <- Initial
+ .Object at Terminal <- Terminal
.Object at division <- division
- .Object at Initial <- Initial
+ .Object at delta <- delta
.Object at grid <- grid
- .Object at random <- random
- }else if(length(Terminal)==1){
- .Object at division <- division
- .Object at Terminal <- rep(Terminal, length(division))
- .Object at Initial <- Initial
+ .Object at random <- random
+ .Object at regular <- regular
+ }else if(length(Terminal)==1){
+ .Object at Initial <- Initial
+ .Object at Terminal <- rep(Terminal, length(division))
+ .Object at division <- division
+ .Object at delta <- delta
.Object at grid <- grid
- .Object at random <- random
+ .Object at random <- random
+ .Object at regular <- regular
}else if(length(division)==1){
- .Object at Terminal <- Terminal
+ .Object at Initial <- Initial
+ .Object at Terminal <- Terminal
.Object at division <- rep(division, length(Terminal))
- .Object at Initial <- Initial
+ .Object at delta <- delta
.Object at grid <- grid
- .Object at random <- random
+ .Object at random <- random
+ .Object at regular <- regular
}else{
cat("\nDimension missmatch.\n")
return(NULL)
- }
+ }}
return(.Object)
})
-setSampling <- function(Terminal, division, Initial=0, grid=as.numeric(NULL), random=FALSE){
- return(new("yuima.sampling", Terminal=Terminal, division=division, Initial=Initial, grid=grid, random=random))
+setSampling <- function(Initial=0, Terminal=1, division=100, delta=0.1, grid=as.numeric(NULL), random=FALSE){
+ return(new("yuima.sampling", Initial=Initial, Terminal=Terminal, division=division, delta=delta, grid=grid, random=random, regular=TRUE))
}
# accessors
More information about the Yuima-commits
mailing list