[Yuima-commits] r326 - in pkg/yuima: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sat Sep 6 12:05:57 CEST 2014


Author: lorenzo
Date: 2014-09-06 12:05:56 +0200 (Sat, 06 Sep 2014)
New Revision: 326

Modified:
   pkg/yuima/DESCRIPTION
   pkg/yuima/R/AllClasses.R
   pkg/yuima/R/qmle.R
   pkg/yuima/R/yuima.R
   pkg/yuima/man/yuima.carma.qmle-class.Rd
Log:
Summary fo qmle.carma

Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION	2014-09-03 17:30:59 UTC (rev 325)
+++ pkg/yuima/DESCRIPTION	2014-09-06 10:05:56 UTC (rev 326)
@@ -1,8 +1,8 @@
 Package: yuima
 Type: Package
 Title: The YUIMA Project package for SDEs
-Version: 1.0.28
-Date: 2014-09-04
+Version: 1.0.29
+Date: 2014-09-06
 Depends: methods, zoo, stats4, utils, expm, cubature, mvtnorm
 Author: YUIMA Project Team
 Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>

Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R	2014-09-03 17:30:59 UTC (rev 325)
+++ pkg/yuima/R/AllClasses.R	2014-09-06 10:05:56 UTC (rev 326)
@@ -123,7 +123,8 @@
 
 # Class yuima.carma.qmle
 setClass("yuima.carma.qmle",representation(Incr.Lev = "ANY",
-                                           model = "yuima.carma"
+                                           model = "yuima.carma",
+                                           logL.Incr = "ANY"
                                            ),
                             contains="mle"
          )
@@ -141,10 +142,13 @@
 contains="mle"
 )
 
-setClass("yuima.carma.qmle",representation(Incr.Lev = "ANY",
-model = "yuima.carma"
-),
-contains="mle"
+setClass("summary.yuima.carma.qmle",representation(MeanI = "ANY",
+                                                   SdI = "ANY",
+                                                   logLI = "ANY",
+                                                   TypeI = "ANY",
+                                                   NumbI = "ANY",
+                                                   StatI ="ANY"),
+contains="summary.mle"
 )
 
 setClass("summary.yuima.CP.qmle",

Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R	2014-09-03 17:30:59 UTC (rev 325)
+++ pkg/yuima/R/qmle.R	2014-09-06 10:05:56 UTC (rev 326)
@@ -141,7 +141,7 @@
   # In this case we use a two step procedure:
   # First) The Coefficient are obtained by QMLE computed using the Kalman Filter.
   # Second) Using the result in Brockwell, Davis and Yang (2007) we retrieve 
-  # the underlying Levy. The estimated increments are used to find the Lévy parameters.
+  # the underlying Levy. The estimated increments are used to find the L?vy parameters.
 
 #   if(is(yuima at model, "yuima.carma")){
 #     yuima.warm("two step procedure for carma(p,q)")
@@ -521,7 +521,7 @@
               new.start <- start[-c(idx.fixed,idx.measure)] # considering only initial guess for
             }
 
-            if(length(new.start)>1){ #Â?multidimensional optim # Adjust lower for no negative Noise
+            if(length(new.start)>1){ #??multidimensional optim # Adjust lower for no negative Noise
                 if(is.CARMA(yuima) && (NoNeg.Noise==TRUE))
                     if(mean.noise %in% names(lower)){lower[mean.noise]<-10^-7}
 				oout <- optim(new.start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
@@ -919,10 +919,14 @@
     min.jump <- 0
     
     
-    if(length(c(diff.par,drift.par))>0)
-	  min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env)
+    if(length(c(diff.par,drift.par))>0 & !is.CARMA(yuima)){ # LM 04/09/14
+	    min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env)
+    }else{
+      if(length(c(diff.par,drift.par))>0 & is.CARMA(yuima)){
+        min.diff <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+      }
+    }
     
-    
     if(length(c(measure.par))>0 & !is.CARMA(yuima))
         min.jump <-   minusquasipsi(yuima=yuima, param=mycoef[measure.par], print=print, env=env)
     
@@ -960,7 +964,7 @@
     if( Est.Incr=="Carma.IncPar" || Est.Incr=="Carma.Inc" ){
     final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef), 
                    vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl, 
-                    method = method, nobs=yuima.nobs)
+                    method = method, nobs=yuima.nobs, logL.Incr = NULL)
     }else{
       if(Est.Incr=="Carma.Par"){
       final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef), 
@@ -1032,7 +1036,7 @@
      carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef), 
                           vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl, 
                           method = method, Incr.Lev = inc.levy.fin,
-                          model = yuima at model, nobs=yuima.nobs)
+                          model = yuima at model, nobs=yuima.nobs, logL.Incr = NULL)
      return(carma_final_res)
    }
    
@@ -1087,7 +1091,7 @@
           carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef), 
                                vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl, 
                                method = method, Incr.Lev = inc.levy,
-                               model = yuima at model)
+                               model = yuima at model, logL.Incr = NULL)
           return(carma_final_res)
         }
         
@@ -1189,7 +1193,7 @@
           carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef), 
                                vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl, 
                                method = method, Incr.Lev = inc.levy,
-                               model = yuima at model)
+                               model = yuima at model, logL.Incr = NULL)
           return(carma_final_res)
         }
         if(aggregation==TRUE){
@@ -1313,7 +1317,8 @@
       carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(coef), 
                      vcov = cov, min = min, details = oout, minuslogl = minusquasilogl, 
                      method = method, Incr.Lev = inc.levy.fin,
-                           model = yuima at model, nobs=yuima.nobs)
+                           model = yuima at model, nobs=yuima.nobs,
+                     logL.Incr = tryCatch(-result.Lev$value,error=function(theta){NULL}))
     }else{
       if(Est.Incr=="Carma.Par"){
         carma_final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(coef), 
@@ -2009,20 +2014,6 @@
 
 
 
-
-# Plot Method for yuima.carma.qmle
-setMethod("plot",signature(x="yuima.carma.qmle"),
-          function(x, ...){
-            Time<-index(x at Incr.Lev)
-            Incr.L<-coredata(x at Incr.Lev)
-            if(is.complex(Incr.L)){
-              yuima.warn("Complex increments. We plot only the real part")
-              Incr.L<-Re(Incr.L)
-            }
-            plot(x=Time,y=Incr.L, ...)
-          }
-)
-
 setMethod("summary", "yuima.qmle",
 function (object, ...)
 {
@@ -2109,6 +2100,68 @@
 )
 # Utilities for estimation of levy in continuous arma model
 
+setMethod("summary", "yuima.carma.qmle",
+          function (object, ...)
+          {
+            cmat <- cbind(Estimate = object at coef, `Std. Error` = sqrt(diag(object at vcov)))
+            m2logL <- 2 * object at min
+            data<-Re(coredata(object at Incr.Lev))
+            data<- data[!is.na(data)]
+            
+            tmp <- new("summary.yuima.carma.qmle", call = object at call, coef = cmat,
+                       m2logL = m2logL, 
+                       MeanI = mean(data),
+                       SdI = sd(data),
+                       logLI = object at logL.Incr,
+                       TypeI = object at model@measure.type,
+                       NumbI = length(data),
+                       StatI =summary(data)
+            )
+            tmp
+          }
+)
+
+setMethod("show", "summary.yuima.carma.qmle",
+          function (object)
+          {
+            
+            cat("Two Stage Quasi-Maximum likelihood estimation\n\nCall:\n")
+            print(object at call)
+            cat("\nCoefficients:\n")
+            print(coef(object))
+            cat("\n-2 log L:", object at m2logL, "\n")
+            
+            cat(sprintf("\n\nNumber of increments: %d\n",object at NumbI))
+            cat(sprintf("\nAverage of increments: %f\n",object at MeanI))
+            cat(sprintf("\nStandard Dev. of increments: %f\n",object at SdI))
+            if(!is.null(object at logLI)){
+              cat(sprintf("\n\n-2 log L of increments: %f\n",-2*object at logLI))
+            }
+            cat("\nSummary statistics for increments:\n")
+            print(object at StatI)
+            cat("\n")
+          }
+)
+
+
+
+  # Plot Method for yuima.carma.qmle
+setMethod("plot",signature(x="yuima.carma.qmle"),
+          function(x, ...){
+            Time<-index(x at Incr.Lev)
+            Incr.L<-coredata(x at Incr.Lev)
+            if(is.complex(Incr.L)){
+              yuima.warn("Complex increments. We plot only the real part")
+              Incr.L<-Re(Incr.L)
+            }
+            plot(x=Time,y=Incr.L, ...)
+          }
+)
+
+
+
+
+
 #Density code for compound poisson
 
 #CPN
@@ -2571,7 +2624,7 @@
         paramLev[1]<-paramLev[1]/dt
       }
   }
-  results<-list(estLevpar=paramLev,covLev=covLev)
+  results<-list(estLevpar=paramLev,covLev=covLev, value=firs.prob$value)
   return(results)
 }
 

Modified: pkg/yuima/R/yuima.R
===================================================================
--- pkg/yuima/R/yuima.R	2014-09-03 17:30:59 UTC (rev 325)
+++ pkg/yuima/R/yuima.R	2014-09-06 10:05:56 UTC (rev 326)
@@ -243,6 +243,13 @@
 # setter
 setYuima <-
   function(data=NULL, model=NULL, sampling=NULL, characteristic=NULL, functional=NULL){
+    if(is.CARMA(model)&& !is.null(data)){
+      if(is.null(dim(data at original.data))){
+        data<-setData(zoo(x=matrix(as.numeric(data at original.data),length(data at original.data),
+                                   (model at info@p+1)), order.by=time(data at zoo.data[[1]])))  
+      
+      } 
+    }# LM 05/09/14 
     return(new("yuima", data=data, model=model, sampling=sampling, characteristic=characteristic,functional=functional))
   }
 

Modified: pkg/yuima/man/yuima.carma.qmle-class.Rd
===================================================================
--- pkg/yuima/man/yuima.carma.qmle-class.Rd	2014-09-03 17:30:59 UTC (rev 325)
+++ pkg/yuima/man/yuima.carma.qmle-class.Rd	2014-09-06 10:05:56 UTC (rev 326)
@@ -14,6 +14,7 @@
   \describe{
     \item{\code{Incr.Lev}:}{is an object of class \code{\link{zoo}} that contains the estimated increments of the noise obtained using \code{\link{CarmaNoise}}.}
     \item{\code{model}:}{is an object of of \code{\link{yuima.carma-class}}.}
+    \item{\code{logL.Incr}:}{is an object of class \code{numeric} that contains the value of the log-likelihood for estimated Levy increments.}
     \item{\code{call}:}{is an object of class \code{language}. }
     \item{\code{coef}:}{is an object of class \code{numeric} that contains estimated parameters.}
     \item{\code{fullcoef}:}{is an object of class \code{numeric} that contains estimated and fixed parameters.}



More information about the Yuima-commits mailing list