[Yuima-commits] r30 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 24 19:32:20 CET 2009
Author: abrouste
Date: 2009-11-24 19:32:20 +0100 (Tue, 24 Nov 2009)
New Revision: 30
Added:
pkg/yuima/R/CholeskyfGn.R
pkg/yuima/R/sampling2grid.R
Modified:
pkg/yuima/R/AllClasses.R
pkg/yuima/R/simulate.R
pkg/yuima/R/yuima.model.R
pkg/yuima/R/yuima.sampling.R
Log:
add fractional hurst index in class + simulate with fgN
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2009-11-24 18:17:47 UTC (rev 29)
+++ pkg/yuima/R/AllClasses.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -52,7 +52,10 @@
# type, rate, deltas, etc.
setClass("yuima.sampling", representation(Terminal = "numeric",
- division = "numeric"
+ division = "numeric",
+ Initial = "numeric",
+ grid = "numeric",
+ random = "logical"
)
)
Added: pkg/yuima/R/CholeskyfGn.R
===================================================================
--- pkg/yuima/R/CholeskyfGn.R (rev 0)
+++ pkg/yuima/R/CholeskyfGn.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -0,0 +1,38 @@
+CholeskyfGn <-
+function(mesh, H)
+{
+
+##--------------------------------------------------------------------------
+## Author : Alexandre Brouste
+## Project: Yuima and Fieldsim
+##--------------------------------------------------------------------------
+
+##--------------------------------------------------------------------------
+## Input : mesh : mesh grid where the fBm is evaluated
+## H : self-similarity parameter
+##
+##
+## Output : simulation of a standard fractional Brownian noise
+## for the mesh grid evaluation by Choleki s
+## decomposition of the covariance matrix of the fGn.
+##--------------------------------------------------------------------------
+
+##---------------------------------------------------------
+## Complexity O(N^3) via method chol
+## Taille memoire N^2
+## -------------------------------------------------------------------------
+
+ N<-length(mesh)-1 #j'ai retirer 1 # N+1 is the size of the fGn sample to be simulated
+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)
+ matcov[j, i]<- matcov[i, j]
+ }
+ L <- chol(matcov)
+ Z <- rnorm(N+1)
+ fGn <- t(L) %*% Z
+ return(fGn)
+}
\ No newline at end of file
Property changes on: pkg/yuima/R/CholeskyfGn.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/yuima/R/sampling2grid.R
===================================================================
--- pkg/yuima/R/sampling2grid.R (rev 0)
+++ pkg/yuima/R/sampling2grid.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -0,0 +1,19 @@
+sampling2grid <- function(samp){
+
+#Check if there is already a grid
+ if (length(samp at grid)>0){
+
+#It is the general case where there is a grid
+#including irregular deterministic case
+ return(samp at grid)
+ }
+
+#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))
+ }
+
+ return(1) #To do: Random grid
+}
+
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2009-11-24 18:17:47 UTC (rev 29)
+++ pkg/yuima/R/simulate.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -126,15 +126,21 @@
##:: using Euler-Maruyama method
delta <- Terminal/division
- ##:: Diffusion terms
- if( missing(increment.W)){
-
- dW <- rnorm(division * r.size, 0, sqrt(delta))
- dW <- t(matrix(dW, nrow=division, ncol=r.size))
- }else{
- dW <- increment.W
- }
-
+ if (missing(increment.W)){
+ if (sdeModel at hurst!=1/2){
+ grid<-sampling2grid(yuima at sampling)
+ dW<-CholeskyfGn(grid,sdeModel at hurst)
+ dW <- t(matrix(dW, nrow=division, ncol=r.size))
+ }else{
+ delta<-Terminal/division
+ dW <- rnorm(division * r.size, 0, sqrt(delta))
+ dW <- t(matrix(dW, nrow=division, ncol=r.size))
+ }
+
+ }else{
+ dW <- increment.W
+ }
+
yuima at data <- euler(xinit, yuima, dW, yuimaEnv)
return(yuima)
})
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2009-11-24 18:17:47 UTC (rev 29)
+++ pkg/yuima/R/yuima.model.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -19,6 +19,7 @@
function(.Object,
drift,
diffusion,
+ hurst,
jump.coeff,
measure,
measure.type,
@@ -33,6 +34,7 @@
xinit){
.Object at drift <- drift
.Object at diffusion <- diffusion
+ .Object at hurst <- hurst
.Object at jump.coeff <- jump.coeff
.Object at measure <- measure
.Object at measure.type <- measure.type
@@ -53,6 +55,7 @@
## set yuima model from SDE
setModel <- function(drift=NULL,
diffusion=NULL,
+ hurst=as.numeric(NULL),
jump.coeff=character(),
measure=list(),
measure.type=character(),
@@ -379,6 +382,7 @@
tmp <- new("yuima.model",
drift= DRIFT,
diffusion= DIFFUSION,
+ hurst=as.numeric(hurst),
jump.coeff=JUMP,
measure= MEASURE,
measure.type= measure.type,
Modified: pkg/yuima/R/yuima.sampling.R
===================================================================
--- pkg/yuima/R/yuima.sampling.R 2009-11-24 18:17:47 UTC (rev 29)
+++ pkg/yuima/R/yuima.sampling.R 2009-11-24 18:32:20 UTC (rev 30)
@@ -5,17 +5,26 @@
setMethod("initialize", "yuima.sampling",
- function(.Object, Terminal, division){
+ function(.Object, Terminal, division, Initial, grid, random){
eqn <- length(Terminal)
if(length(Terminal)==length(division)){
.Object at Terminal <- Terminal
.Object at division <- division
- }else if(length(Termianl)==1){
+ .Object at Initial <- Initial
+ .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 grid <- grid
+ .Object at random <- random
}else if(length(division)==1){
.Object at Terminal <- Terminal
.Object at division <- rep(division, length(Terminal))
+ .Object at Initial <- Initial
+ .Object at grid <- grid
+ .Object at random <- random
}else{
cat("\nDimension missmatch.\n")
return(NULL)
@@ -23,8 +32,8 @@
return(.Object)
})
-setSampling <- function(Terminal, division){
- return(new("yuima.sampling", Terminal=Terminal, division=division))
+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))
}
# accessors
More information about the Yuima-commits
mailing list