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

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

Author: lorenzo
Date: 2021-10-11 21:39:50 +0200 (Mon, 11 Oct 2021)
New Revision: 760

Modified:
   pkg/yuima/R/qmleLevy.R
Log:
Standard errors qmleLevy.R

Modified: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R	2021-10-09 14:49:51 UTC (rev 759)
+++ pkg/yuima/R/qmleLevy.R	2021-10-11 19:39:50 UTC (rev 760)
@@ -604,12 +604,92 @@
     upperjump <- upper0[lev.names]
     esti <- optim(fn = minusloglik, lower = lowerjump, upper = upperjump, 
                   par = para, method = "L-BFGS-B")
-    #optimHess(par=esti$par, fn=minusloglik)
+
+    HessianEta <- optimHess(par=esti$par, fn=minusloglik)
     res at coef<-c(res at coef,esti$par)
     res at fullcoef[names(para)]<-esti$par
-    res at vcov<-cbind(res at vcov,matrix(NA,ncol=length(esti$par),nrow=dim(res at vcov)[1]))
+    if(length(oldyuima at data@zoo.data)==1 & is(oldyuima at model@measure$df, "yuima.law")){
+      if(!aggregation){
+        Ter <- yuima at sampling@Terminal
+        ures <- numeric(floor(Ter))
+        for(l in 1:floor(Ter)){
+          ures[l] <- sum(resi[(floor((l-1)/h)):(floor(l/h)-1)]) 
+        } 
+      }
+      mypar<-res at coef[oldyuima at model@parameter at measure]
+      mypar[oldyuima at model@measure$df at time.var]<-1
+      fdataeta<-dens(oldyuima at model@measure$df,x=ures,mypar)# f(eps, eta)
+      del <- 10^-3
+      fdatadeltaeta<- dens(oldyuima at model@measure$df,x=ures+del,mypar) # f(eps + delta, eta)
+      dummy <- t(rep(1,length(mypar[oldyuima at model@measure$df at param.measure])))
+      myeta<- as.matrix(mypar[oldyuima at model@measure$df at param.measure])%*%dummy
+      myetapert <- myeta+del*diag(rep(1,dim(myeta)[1])) 
+      fdataetadelta<-sapply(X=1:dim(myeta)[1],FUN = function(X){
+        par<-myetapert[,X]
+        par[oldyuima at model@measure$df at time.var]<-1
+        dens(oldyuima at model@measure$df,x=ures,par)
+      }
+      )# f(eps, eta+delta)
+      fdatadeltaetadelta<-sapply(X=1:dim(myeta)[1],FUN = function(X){
+        par<-myetapert[,X]
+        par[oldyuima at model@measure$df at time.var]<-1
+        dens(oldyuima at model@measure$df,x=ures+del,par)
+      }
+      ) # f(eps +deta, eta+delta) 
+      term1<-1/(fdataeta)
+      term2 <- fdatadeltaeta*term1%*%dummy
+      
+      term2 <- fdatadeltaetadelta-term2*fdataetadelta
+      mixpartial<-t((as.matrix(term1)%*%dummy)/del^2*term2/oldyuima at sampling@Terminal)
+      
+      
+      # construction of b_i 
+      # DiffJumpCoeff, DriftDerCoeff, jump.term length(resi)
+      # step1 <- t(DiffJumpCoeff)%*%DiffJumpCoeff
+      DerMeta <- 1/del*(fdataetadelta - as.matrix(fdataeta)%*%rep(1,dim(fdataetadelta)[2]))*(as.matrix(term1)%*%rep(1,dim(fdataetadelta)[2]))
+      SigmaEta <- t(DerMeta)%*%DerMeta/oldyuima at sampling@Terminal
+      
+      SigmaEtaAlpha<- 1/oldyuima at sampling@n*DriftDerCoeff%*%(t(DiffJumpCoeff)/(as.matrix(jump.term[-length(jump.term)]^2)%*%rep(1,dim(DiffJumpCoeff)[1])) )
+      SigmaEtaAlpha <- SigmaEtaAlpha*sum(resi^3)/oldyuima at sampling@delta
+      
+      b_i <- matrix(0,floor(Ter),length(c(oldyuima at model@parameter at drift, oldyuima at model@parameter at jump)))
+      Coef1 <- matrix(0,floor(Ter),length( oldyuima at model@parameter at jump))
+      Coef2 <- matrix(0,floor(Ter),length( oldyuima at model@parameter at drift))
+      
+      for(l in 1:floor(Ter)){
+        pos <- (floor((l-1)/h)):(floor(l/h)-1)
+        if(length(oldyuima at model@parameter at jump)==1){
+          b_i[l,1:length(oldyuima at model@parameter at jump)] <- sum(-DiffJumpCoeff[,pos]*resi[pos]/jump.term[pos])
+          Coef1[l,]<-sum(DiffJumpCoeff[,pos]/jump.term[pos]*(resi[pos]^2-h))
+        }else{
+          interm <- as.matrix(resi[pos]/jump.term[pos])
+          b_i[l,1:length(oldyuima at model@parameter at jump)] <-  -t(DiffJumpCoeff[,pos]%*%interm)
+          interm2 <- as.matrix((resi[pos]^2-h)/jump.term[pos])
+          Coef1[l,] <-DiffJumpCoeff[,pos]%*%interm2
+        }
+        if(length(oldyuima at model@parameter at drift)==1){
+          b_i[l,1:length(oldyuima at model@parameter at drift)+length(oldyuima at model@parameter at jump)] <- -sum(h* DriftDerCoeff[,pos]%*%jump.term[pos])
+          Coef2[l,]<-sum(DriftDerCoeff[,pos]/jump.term[pos]*(resi[pos]))
+        }else{
+          b_i[l,1:length(oldyuima at model@parameter at drift)+length(oldyuima at model@parameter at jump)] <--h* DriftDerCoeff[,pos]%*%jump.term[pos]
+          interm3 <- as.matrix((resi[pos])/jump.term[pos])
+          Coef2[l,] <-DriftDerCoeff[,pos]%*%interm3
+        }
+      }
+      MatrUnder <- mixpartial%*%b_i
+      I_n <- cbind(rbind(myGamhat,MatrUnder),rbind(matrix(0, dim(myGamhat)[1],dim(HessianEta)[2]),HessianEta/oldyuima at sampling@Terminal))
+      SigmaGammaEta <- t(DerMeta)%*%Coef1/Ter
+      SigmaAlphaEta <- t(DerMeta)%*%Coef2/Ter
+      # dim(fdataetadelta), length(fdataeta),  length(term1)
+      dum <- cbind(SigmaGammaEta , SigmaAlphaEta)
+      MatSigmaHat1<- rbind(cbind(MatSigmaHat,t(dum)),cbind(dum,SigmaEta)) 
+      InvIn<- solve(I_n)
+      res at vcov <-InvIn%*%MatSigmaHat1%*%t(InvIn)/Ter
+    }else{
+      res at vcov<-cbind(res at vcov,matrix(NA,ncol=length(esti$par),nrow=dim(res at vcov)[1]))
+    }
     colnames(res at vcov)<-names(res at fullcoef)
-    res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
+    #res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
     rownames(res at vcov)<-names(res at fullcoef)
     res at min<-c(res at min,esti$value)
     res at nobs<-c(res at nobs,length(Incr.Lev at zoo.data[[1]]))