[Yuima-commits] r244 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 10 16:41:32 CEST 2013
Author: abrouste
Date: 2013-10-10 16:41:32 +0200 (Thu, 10 Oct 2013)
New Revision: 244
Added:
pkg/yuima/R/mmfrac.R
Modified:
pkg/yuima/R/AllClasses.R
pkg/yuima/R/qgv.R
pkg/yuima/R/simulate.R
pkg/yuima/R/zzz.R
Log:
Fractional estimation
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2013-07-18 10:42:25 UTC (rev 243)
+++ pkg/yuima/R/AllClasses.R 2013-10-10 14:41:32 UTC (rev 244)
@@ -16,7 +16,7 @@
# Class 'yuima.model'
setClass("yuima.model",representation(drift="expression",
diffusion="list",
- hurst="numeric",
+ hurst="ANY",
jump.coeff="expression",
measure="list",
measure.type="character",
Added: pkg/yuima/R/mmfrac.R
===================================================================
--- pkg/yuima/R/mmfrac.R (rev 0)
+++ pkg/yuima/R/mmfrac.R 2013-10-10 14:41:32 UTC (rev 244)
@@ -0,0 +1,84 @@
+##estimate theta2 by LSE
+
+mmfrac <- function(yuima,...){
+
+ call <- match.call()
+
+ if(missing(yuima)){
+ yuima.stop("yuima object is missing.")
+ }
+
+ if(class(yuima)!="yuima"){
+ yuima.stop("an object of class yuima is needed.")
+ }
+
+ Ddiffx0 <- D(yuima at model@diffusion[[1]],yuima at model@state.variable) == 0
+ Ddifft0 <- D(yuima at model@diffusion[[1]],yuima at model@time.variable) == 0
+
+ bconst<-FALSE
+ Hconst<-FALSE
+
+ diffnoparam<-length(yuima at model@parameter at diffusion)==0
+
+ if (Ddiffx0 & Ddifft0 & diffnoparam){
+ if(eval(yuima at model@diffusion[[1]])==0){
+ yuima.stop("the model has no fractional part.")
+ }
+ }
+
+ if (length(yuima at model@parameter at drift)!=1){
+ yuima.stop("the drift is malformed.")
+ }
+
+
+ dx <- D(yuima at model@drift,yuima at model@state.variable)
+ dbx <- D(yuima at model@diffusion[[1]],yuima at model@state.variable)==0
+ dbt <- D(yuima at model@diffusion[[1]],yuima at model@time.variable)==0
+ bxx <- D(dx,yuima at model@state.variable)==0
+ bmix <- D(dx,yuima at model@time.variable)==0
+
+ isfOU <- (bxx || bmix) && (dbx && dbt)
+
+ if (!isfOU){yuima.stop("estimation not available for this model")}
+
+ process<-yuima at data@zoo.data[[1]]
+
+
+ T<-yuima at sampling@Terminal
+ est<-qgv(yuima,...)
+
+ H<-est$coeff[1]
+ sigma<-est$coeff[2]
+
+ if ((!is.na(H))&(!is.na(sigma))){
+
+ theta<-(2*mean(process^2)/(sigma^2*gamma(2*H+1)))^(-1/2/H)
+ sH<-sqrt((4*H-1)*(1+(gamma(1-4*H)*gamma(4*H-1))/(gamma(2-2*H)*gamma(2*H))))
+ sdtheta<- sqrt(theta)*(sH/(2*H))/sqrt(T)
+
+ }else{
+ yuima.warn("Diffusion estimation not available, can not estimate the drift parameter")
+ }
+
+ x<-c(est$coeff,theta)
+ names(x)<-c(names(est$coeff),yuima at model@parameter at drift)
+ sdx<-matrix(,3,3)
+ diag(sdx)<-c(diag(est$vcov),sdtheta)
+ colnames(sdx)<-names(x)
+ rownames(sdx)<-names(x)
+
+ obj <- list(coefficients=x,vcov=sdx,call=call)
+ class(obj) <- "mmfrac"
+ return(obj)
+
+}
+
+print.mmfrac<-function(x,...){
+print(x)
+}
+
+
+}
+
+
+
Modified: pkg/yuima/R/qgv.R
===================================================================
--- pkg/yuima/R/qgv.R 2013-07-18 10:42:25 UTC (rev 243)
+++ pkg/yuima/R/qgv.R 2013-10-10 14:41:32 UTC (rev 244)
@@ -1,18 +1,85 @@
-#Estimating local Hölder exponent with the constant.
+##:: function qgv
+##:: Estimating the local Hölder exponent of the path and the constant
+qgv<- function(yuima, filter.type="Daubechies", order=2, a=NULL, ...){
-qgv<- function(yuima, filter.type="Daubechies", order=2, with.constant.est=TRUE, a=NULL, ...){
-
call <- match.call()
if(missing(yuima)){
yuima.stop("yuima object is missing.")
}
+
+ if(class(yuima)!="yuima"){
+ yuima.stop("an object of class yuima is needed.")
+ }
+
-#Test for order to be integer
-
- if (missing(a)){
+ Ddiffx0 <- D(yuima at model@diffusion[[1]],yuima at model@state.variable) == 0
+ Ddifft0 <- D(yuima at model@diffusion[[1]],yuima at model@time.variable) == 0
+
+ bconst<-FALSE
+ Hconst<-FALSE
+
+ diffnoparam<-length(yuima at model@parameter at diffusion)==0
+
+ if (Ddiffx0 & Ddifft0 & diffnoparam){
+
+ if(eval(yuima at model@diffusion[[1]])==0){
+ yuima.stop("the model has no fractional part.")
+ }
+
+ bconst<-TRUE
+ }
+
+ y<-yuima at model@drift
+ dx <- D(y,yuima at model@state.variable)
+ bt <- D(y,yuima at model@time.variable)==0
+ bxx <- D(dx,yuima at model@state.variable)==0
+
+
+ dbx <- D(yuima at model@diffusion[[1]],yuima at model@state.variable)==0
+ dbt <- D(yuima at model@diffusion[[1]],yuima at model@time.variable)==0
+
+
+ isfOU<-(bxx && bt) && (dbx && dbt)
+ if(isfOU){yuima.warn("It is fOU")} #To be deleted
+
+
+ if (!is.na(yuima at model@hurst)){
+ yuima.warn("No Hurst exponent will be estimated")
+ Hconst<-TRUE
+ }
+
+
+ H <- yuima at model@hurst
+ sdH<-NA
+
+ if(Hconst & bconst){
+ x<-c(H,eval(yuima at model@diffusion[[1]]))
+ names(x)<-c("hurst","const")
+ sdx<-diag(c(NA,NA))
+ colnames(sdx)<-names(x)
+ rownames(sdx)<-names(x)
+ sdx[2,1] <- sdx[1,2] <- NA
+ obj <- list(coefficients=x,vcov=sdx,call=call)
+ class(obj) <- "qgv"
+ return(obj)
+ }
+
+ isregular<-yuima at sampling@regular
+
+ if (!isregular){
+ yuima.stop("qgv method is only working for regular grid.")
+ }
+
+ if (!(order %in% 1:10)){
+ yuima.warn("Classical filter implement are of order ranged in [1,10], order have been fixed to 2.")
+ order=2
+ }
+
+ if (missing(a)){
+
if (filter.type=="Daubechies"){
if (order==2){
@@ -27,16 +94,25 @@
-.8365163037378077,
.2241438680420134,
.1294095225512603)
+ order=2
}
}else if (filter.type=="Classical"){
mesh<-0:order
a=(-1)^(1-mesh)/2^order*choose(order,mesh)
+ }else{
+ yuima.warn("No such type of filter. Filter have been fixed to Daubechies' filter of order 2.")
+ a<-1/sqrt(2)*c(.4829629131445341,
+ -.8365163037378077,
+ .2241438680420134,
+ .1294095225512603)
+ order=2
+ }
+
+
}
- }
-
L<-length(a)
a2<-rep(0,2*L)
a2[seq(1,2*L,by=2)]<-a
@@ -44,26 +120,108 @@
process<-yuima at data@zoo.data[[1]]
N<-length(process)
-
#Computation of the generalized quadratic variations
V1<-sum(filter(process,a)^2,na.rm=TRUE)
V2<-sum(filter(process,a2)^2,na.rm=TRUE)
-
- H<-1/2*log2(V2/V1)
-
- if (H>1){H<-0.999} #Over-estimation of H
- if (H<0){H<-0.001}
-
-
- #Compute the estimation of the constant C.
-
- delta<-yuima at sampling@delta
-
- constfilt<-sum(a%*%t(a)*abs(matrix(0:(L-1),L,L)-matrix(0:(L-1),L,L,byrow=TRUE))^(2*H))
-
- C<- sqrt(-2*V1/(N-L)/(constfilt*delta^(2*H)))
+
+ if(!Hconst){
+ H<-1/2*log2(V2/V1)
+ }
+
+ nconst <- "const"
+
+ sdC<-NA
+ C <- NA
+ if(isfOU){
+ if( bconst||(H>=1)||(H<=0)){
+ if(diffnoparam){
+ C<-eval(yuima at model@diffusion[[1]])
+ }
+ }else{
+ #Compute the estimation of the constant C.
+ delta<-yuima at sampling@delta
+ constfilt<-sum(a%*%t(a)*abs(matrix(0:(L-1),L,L)-matrix(0:(L-1),L,L,byrow=TRUE))^(2*H))
+ C<- sqrt(-2*V1/(N-L)/(constfilt*delta^(2*H)))
+ nconst<-as.character(yuima at model@diffusion[[1]])
+ }
+ } else {
+ nconst<-as.character(yuima at model@diffusion[[1]])
+ }
+
+
+ if(isfOU){
+
+ #Compute the standard error
+ infty<-100
+
+ C11<-rep(0,2*infty+1)
+ C11bis<-rep(0,2*infty+1)
+ C12<-rep(0,2*infty+1)
+ C22<-rep(0,2*infty+1)
+ l<-order+1
+
+ for (i in (-infty:infty)){
+
+ for (q in 0:l){
+ for (r in 0:l){
+ C11[i+infty+1]<-C11[i+infty+1]+a[q+1]*a[r+1]*abs(q-r+i)^(2*H)
+ }
+ }
+
+ for (q in 0:l){
+ for (r in 0:l){
+ C12[i+infty+1]<-C12[i+infty+1]+a[q+1]*a[r+1]*abs(2*q-r+i)^(2*H)
+ }
+ }
+
+ for (q in 0:l){
+ for (r in 0:l){
+ C22[i+infty+1]<-C22[i+infty+1]+a[q+1]*a[r+1]*abs(2*q-2*r+i)^(2*H)
+ }
+ }
+
+ for (q in 0:(2*l)){
+ for (r in 0:(2*l)){
+ C11bis[i+infty+1]<-C11bis[i+infty+1]+a2[q+1]*a2[r+1]*abs(q-r+i)^(2*H)
+ }
+ }
+
+ }
+
+ rho11<-1/2*sum((C11/C11[infty+1])^2)
+ rho11dil<-1/2*sum((C11bis/C11bis[infty+1])^2)
+
+ rho12<-1/2*sum((C12/2^H/C11[infty+1])^2)
+ rho22<-1/2*sum((C22/2^H/2^H/C11[infty+1])^2)
+
+ sigma1<-1/(log(2))^2*(rho11+rho11dil-2*rho12)
+
+ if((!Hconst)&(H>0)){
+ sdH<-sqrt(sigma1)/sqrt(N)
+ yuima.warn("sdH computed") #to be deleted
+
+ }
+
+ if ((H>0)&(H<1)&(!bconst)){
+ sigma2<-C^2/2*sigma1
+ sdC<-sqrt(sigma2)/sqrt(N)*log(N)
+ }
- return(list(C=C,H=H))
+ }
+
+x<-c(H,C)
+names(x)<-c("hurst",nconst)
+sdx<-diag(c(sdH,sdC))
+colnames(sdx)<-names(x)
+rownames(sdx)<-names(x)
+sdx[2,1] <- sdx[1,2] <- NA
+obj <- list(coefficients=x,vcov=sdx,call=call)
+class(obj) <- "qgv"
+return(obj)
+
}
+print.qgv<-function(x,...){
+print(x)
+}
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2013-07-18 10:42:25 UTC (rev 243)
+++ pkg/yuima/R/simulate.R 2013-10-10 14:41:32 UTC (rev 244)
@@ -13,7 +13,7 @@
setGeneric("simulate",
function(object, nsim=1, seed=NULL, xinit, true.parameter, space.discretized=FALSE,
- increment.W=NULL, increment.L=NULL, methodfGn="WoodChan",
+ 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),
@@ -26,7 +26,7 @@
setMethod("simulate", "yuima.model",
function(object, nsim=1, seed=NULL, xinit, true.parameter,
space.discretized=FALSE, increment.W=NULL, increment.L=NULL,
- methodfGn="WoodChan",
+ hurst, methodfGn="WoodChan",
sampling, subsampling,
#Initial = 0, Terminal = 1, n = 100, delta,
# grid, random = FALSE, sdelta=as.numeric(NULL),
@@ -42,19 +42,20 @@
} else {
tmpsamp <- sampling
}
- tmpyuima <- setYuima(model=object, sampling=tmpsamp)
+ tmpyuima <- setYuima(model=object, sampling=tmpsamp)
+
out <- simulate(tmpyuima, nsim=nsim, seed=seed, xinit=xinit,
true.parameter=true.parameter,
space.discretized=space.discretized,
increment.W=increment.W, increment.L=increment.L,
- methodfGn=methodfGn, subsampling=subsampling)
+ hurst=hurst,methodfGn=methodfGn, subsampling=subsampling)
return(out)
})
setMethod("simulate", "yuima",
function(object, nsim=1, seed=NULL, xinit, true.parameter,
space.discretized=FALSE, increment.W=NULL, increment.L=NULL,
- methodfGn="WoodChan",
+ hurst,methodfGn="WoodChan",
sampling, subsampling,
Initial = 0, Terminal = 1, n = 100, delta,
grid = as.numeric(NULL), random = FALSE, sdelta=as.numeric(NULL),
@@ -66,9 +67,20 @@
yuima <- object
if(missing(yuima)){
- yuima.warn("yuima object is missing.")
- return(NULL)
+ yuima.warn("yuima object is missing.")
+ return(NULL)
}
+
+ tmphurst<-yuima at model@hurst
+
+ if(!missing(hurst)){
+ yuima at model@hurst=hurst
+ }
+
+ if (is.na(yuima at model@hurst)){
+ yuima.warn("Specify the hurst parameter.")
+ return(NULL)
+ }
tmpsamp <- NULL
if(is.null(yuima at sampling)){
@@ -178,6 +190,8 @@
if(space.discretized){
##:: using Space-discretized Euler-Maruyama method
yuima at data <- space.discretized(xinit, yuima, yuimaEnv)
+
+ yuima at model@hurst<-tmphurst
return(yuima)
}
@@ -215,6 +229,8 @@
for(i in 1:length(yuima at data@zoo.data))
index(yuima at data@zoo.data[[i]]) <- yuima at sampling@grid[[1]] ## to be fixed
yuima at model@xinit <- xinit
+ yuima at model@hurst <-tmphurst
+
if(missing(subsampling))
return(yuima)
subsampling(yuima, subsampling)
Modified: pkg/yuima/R/zzz.R
===================================================================
--- pkg/yuima/R/zzz.R 2013-07-18 10:42:25 UTC (rev 243)
+++ pkg/yuima/R/zzz.R 2013-10-10 14:41:32 UTC (rev 244)
@@ -2,12 +2,17 @@
#.noGenerics <- TRUE
-.onLoad <- function(libname, pkgname)
+.onAttach <- function(libname, pkgname)
{
# require(methods)
# require(zoo)
- packageStartupMessage(sprintf("\n\n%s\nThis is YUIMA Project package.\nNon stable release. Use at your own risk!!!\nFunctions and documentation may change in the final release\n%s\n\n",paste(rep("#",60),collapse=""),paste(rep("#",60),collapse="")))
+ packageStartupMessage(rep("#",60))
+ packageStartupMessage("This is YUIMA Project package.")
+ packageStartupMessage("Non stable release. Use at your own risk!!!")
+ packageStartupMessage("Functions and documentation may change in the final release")
+ packageStartupMessage(rep("#",60))
+
# require(KernSmooth, quietly=TRUE)
# library.dynam("yuima", pkgname, libname)
}
More information about the Yuima-commits
mailing list