[Yuima-commits] r759 - pkg/yuima/R

[noreply at r-forge.r-project.org]

Author: lorenzo
Date: 2021-10-09 16:49:51 +0200 (Sat, 09 Oct 2021)
New Revision: 759

Modified:
   pkg/yuima/R/qmleLevy.R
Log:
added variance-covariance matrix to qmleLevy output

Modified: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R	2021-09-28 19:51:15 UTC (rev 758)
+++ pkg/yuima/R/qmleLevy.R	2021-10-09 14:49:51 UTC (rev 759)
@@ -10,6 +10,7 @@
                    Est.Incr = "NoIncr",
                    aggregation = TRUE)
 {
+  oldyuima<-yuima #line1
   myjumpname <- yuima at model@jump.variable
   mymeasureparam <- yuima at model@parameter at measure
   if(!(Est.Incr %in% c("NoIncr","Incr","IncrPar")))
@@ -180,6 +181,7 @@
     diffupper <- upper[1:length(yuima at model@parameter at diffusion)]
     difflower <- lower[1:length(yuima at model@parameter at diffusion)]
     fres <- qmle(yuima=yuima,start=diffstart,lower=difflower,upper=diffupper,rcpp = TRUE,joint = FALSE,method = "L-BFGS-B")
+    DiffHessian<- fres at details$hessian #182 
     
     DIPAR <- yuima at model@parameter at diffusion
     DIFFUSION <- yuima at model@diffusion
@@ -236,6 +238,7 @@
     yuima at model@parameter at diffusion <- character(0)
     
     sres<-qmle(yuima=yuima,start=dristart,lower=drilower,upper=driupper,rcpp = TRUE,method = "L-BFGS-B")
+    DriftHessian <- sres at details$hessian #239
     
     if((length(ovp) == 0) && (third)){
       yuima at model@diffusion <- DIFFUSION
@@ -291,6 +294,12 @@
                      vcov = vcov0, min = min0, details = details0, minuslogl = minusquasilogl,
                      method = sres at method, nobs=nobs0, model=sdeModel)
       # res <- list(first = fres at coef, second = sres at coef)}
+      if(length(oldyuima at data@zoo.data)==1){
+        myGamhat <- matrix(0,length(coef),length(coef))
+        myGamhat[1:dim(DiffHessian)[1],1:dim(DiffHessian)[2]]<-DiffHessian
+        myGamhat[dim(DiffHessian)[1]+1:dim(DriftHessian)[1],dim(DiffHessian)[1]+1:dim(DriftHessian)[2]]<-DriftHessian#293
+        myGamhat<-myGamhat/oldyuima at sampling@Terminal
+      }
     }else{
       yuima.stop("third estimation may be theoretical invalid under the presence of an overlapping parameter.")
     }
@@ -386,6 +395,11 @@
     # }
     nova<-sqrt((jump.term)^2)
     resi<-(1/(nova[1:(s.size-1)]))*(inc[1:(s.size-1)]-h*drif.term[1:(s.size-1)])
+    if(length(oldyuima at data@zoo.data)==1){
+      coefSigdiff<- 1/h*sum(resi^4)  #389 resi
+      coefDriftSig <- 1/h*sum(resi^3)
+    }
+    
     if(aggregation){
       Ter <- yuima at sampling@Terminal
       ures <- numeric(floor(Ter))
@@ -502,6 +516,64 @@
     Incr.Lev <- setData(original.data=Incr.Lev)
   }
   
+  if(length(oldyuima at data@zoo.data)==1){
+    mydiff<-oldyuima at model@jump.coeff[[1]]
+    mydiffDer <-deriv(mydiff,oldyuima at model@parameter at jump)
+    myenvdiff<- new.env()
+    if(length(oldyuima at model@parameter at jump)>=1){
+      for(i in c(1:length(oldyuima at model@parameter at jump))){
+        assign(value=coef[oldyuima at model@parameter at jump[i]],x=oldyuima at model@parameter at jump[i],envir=myenvdiff)
+      }
+    }
+    EvalPartDiff <- Vectorize(FUN= function(myenvdiff,mydiffDer, data){
+      assign(x=oldyuima at model@solve.variable, value=data,envir = myenvdiff)
+      return(attr(eval(mydiffDer, envir=myenvdiff),"gradient"))},vectorize.args = "data") 
+    
+    DiffJumpCoeff<-EvalPartDiff(myenvdiff,mydiffDer, data=pX)
+    if(!is.matrix(DiffJumpCoeff)){
+      sigmadiffhat<- as.matrix(sum(DiffJumpCoeff^2/jump.term[1:(oldyuima at sampling@n-1)]^2)/(oldyuima at sampling@n))*coefSigdiff
+      DiffJumpCoeff<- t(DiffJumpCoeff)
+    }else{
+      sigmadiffhat<- matrix(0,dim(DiffJumpCoeff)[1],dim(DiffJumpCoeff)[1])
+      for(t in c(1:dim(DiffJumpCoeff)[2])){
+        sigmadiffhat <-sigmadiffhat+as.matrix(DiffJumpCoeff[,t])%*%DiffJumpCoeff[,t]/jump.term[t]^2
+      }
+      sigmadiffhat<-  sigmadiffhat/(oldyuima at sampling@n)*coefSigdiff
+    }
+    
+    mydrift<-oldyuima at model@drift[[1]]
+    mydriftDer <-deriv(mydrift,oldyuima at model@parameter at drift)
+    
+    myenvdrift<- new.env()
+    if(length(oldyuima at model@parameter at drift)>=1){
+      for(i in c(1:length(oldyuima at model@parameter at drift))){
+        assign(value=coef[oldyuima at model@parameter at drift[i]],x=oldyuima at model@parameter at drift[i],envir=myenvdrift)
+      }
+    }
+    
+    DriftDerCoeff<-EvalPartDiff(myenvdrift,mydriftDer, data=pX)
+    if(!is.matrix(DriftDerCoeff)){
+      sigmadrifthat<- as.matrix(sum(DriftDerCoeff^2/jump.term[1:(oldyuima at sampling@n-1)]^2)/(oldyuima at sampling@n))*coefDriftSig
+      DriftDerCoeff <- t(DriftDerCoeff)
+    }else{
+      sigmadrifthat<- matrix(0,dim(DriftDerCoeff)[1],dim(DriftDerCoeff)[1])
+      for(t in c(1:dim(DriftDerCoeff)[2])){
+        sigmadrifthat <-sigmadrifthat+as.matrix(DriftDerCoeff[,t])%*%DriftDerCoeff[,t]/jump.term[t]^2
+      }
+      sigmadrifthat<- sigmadrifthat/(oldyuima at sampling@n)
+    }
+    sigmadriftdiff <- matrix(0, dim(sigmadrifthat)[2],dim(sigmadiffhat)[1])
+    for(t in c(1:dim(DriftDerCoeff)[2]))
+      sigmadriftdiff<-sigmadriftdiff+DriftDerCoeff[,t]%*%t(DiffJumpCoeff[,t])
+    
+    sigmadriftdiff<-sigmadriftdiff/oldyuima at sampling@n
+    
+    MatSigmaHat <- rbind(cbind(sigmadiffhat,t(sigmadriftdiff)),cbind(sigmadriftdiff,sigmadrifthat))
+    
+    res at coef<-res at coef[c(oldyuima at model@parameter at jump,oldyuima at model@parameter at drift)]
+    res at vcov<-solve(t(myGamhat)%*%solve(MatSigmaHat)%*%myGamhat*oldyuima at sampling@Terminal)
+  }
+  
   if(Est.Incr == "Incr"){