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

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

Author: eguchi
Date: 2018-10-30 03:17:37 +0100 (Tue, 30 Oct 2018)
New Revision: 682

Modified:
   pkg/yuima/R/IC.R
Log:
modified

Modified: pkg/yuima/R/IC.R
===================================================================
--- pkg/yuima/R/IC.R	2018-10-30 02:16:59 UTC (rev 681)
+++ pkg/yuima/R/IC.R	2018-10-30 02:17:37 UTC (rev 682)
@@ -1,62 +1,67 @@
 ## information criteria
 
-IC <- function(yuima, data, start, lower, upper, joint = FALSE, rcpp = FALSE,...){
-  if(is(yuima,"yuima") == TRUE){
-    
+IC <- function(yuima, data = NULL, start, lower, upper, joint = FALSE, rcpp = FALSE,...){
+  if(missing(yuima)) stop("yuima object is missing.")
+  
+  if(!is(yuima,"yuima")) stop("This function is for yuima-class.")
+  
+  state.num <- length(yuima at model@state.variable)
+  if(is.null(yuima at data@zoo.data) == TRUE){
     if(is.matrix(data) == FALSE){
+      n <- length(data)
       sub.zoo.data <- list(zoo(x = data, order.by = yuima at sampling@grid[[1]]))
       names(sub.zoo.data)[1] <- "Series 1"
     }else{
       sub.zoo.data <- list()
+      if(ncol(data)-state.num != 0){
+        data <- t(data)
+      }
+      n <- nrow(data)
       for(i in 1:ncol(data)){
         sub.zoo.data <- c(sub.zoo.data, list(zoo(x = data[,i], order.by = yuima at sampling@grid[[1]])))
         names(sub.zoo.data)[i] <- paste("Series", i)
       }
     }
     yuima at data@zoo.data <- sub.zoo.data
-    alpha <- yuima at model@parameter at drift
-    beta <- yuima at model@parameter at diffusion
-    Terminal <- yuima at sampling@Terminal
-    n <- yuima at sampling@n
-    
-    para.num.init  <- match(yuima at model@parameter at all, names(start))
-    para.num.low  <- match(yuima at model@parameter at all, names(lower))
-    para.num.upp  <- match(yuima at model@parameter at all, names(upper))
-    para.start <- NULL
-    para.lower <- NULL
-    para.upper <- NULL
-    for(j in 1:length(yuima at model@parameter at all)){
-      para.start <- c(para.start, list(start[[para.num.init[j]]]))
-      para.lower <- c(para.lower, list(lower[[para.num.low[j]]]))
-      para.upper <- c(para.upper, list(upper[[para.num.upp[j]]]))
-    }
-    names(para.start) <- yuima at model@parameter at all
-    names(para.lower) <- yuima at model@parameter at all
-    names(para.upper) <- yuima at model@parameter at all
-
-    mle <- qmle(yuima, start = para.start, lower = para.lower, upper = para.upper, method = "L-BFGS-B", joint = joint, rcpp = rcpp)
-    hess <- list(mle at details$hessian)
-    hess.diff <- subset(hess[[1]], rownames(hess[[1]])%in%beta, select=beta)
-    hess.drif <- subset(hess[[1]], rownames(hess[[1]])%in%alpha, select=alpha)
-    
-    esti <- coef(mle)
-    names(esti) <- c(beta, alpha)
-    cic <- summary(mle)@m2logL+2*(length(alpha)+length(beta))
-    bic <- summary(mle)@m2logL+length(alpha)*log(Terminal[1])+length(beta)*log(n[1])
-    if(det(hess.diff) > 0 && det(hess.drif) > 0){
-      qbic <- summary(mle)@m2logL+log(det(hess.diff))+log(det(hess.drif))
-    }else{
-      warning(cat("det(hessian of diffusion coefficient)<=0 or det(hessian of drift coefficient)<=0 \n"))
-      qbic <- summary(mle)@m2logL+length(alpha)*log(Terminal)+length(beta)*log(n)
-    }
-    
-    final.res <- list(par = esti, BIC = bic, QBIC = qbic, CIC = cic)
-    return(final.res)
   }else{
-    yuima.stop("This function is for yuima-class.")
+    n <- yuima at sampling@n[1]
   }
+  n <- n-1
+  #alpha <- yuima at model@parameter at drift
+  #beta <- yuima at model@parameter at diffusion
+  Terminal <- yuima at sampling@Terminal[1]
+  
+  para.num.init  <- match(yuima at model@parameter at all, names(start))
+  para.num.low  <- match(yuima at model@parameter at all, names(lower))
+  para.num.upp  <- match(yuima at model@parameter at all, names(upper))
+  para.start <- NULL
+  para.lower <- NULL
+  para.upper <- NULL
+  for(j in 1:length(yuima at model@parameter at all)){
+    para.start <- c(para.start, list(start[[para.num.init[j]]]))
+    para.lower <- c(para.lower, list(lower[[para.num.low[j]]]))
+    para.upper <- c(para.upper, list(upper[[para.num.upp[j]]]))
+  }
+  names(para.start) <- yuima at model@parameter at all
+  names(para.lower) <- yuima at model@parameter at all
+  names(para.upper) <- yuima at model@parameter at all
+  
+  mle <- qmle(yuima, start = para.start, lower = para.lower, upper = para.upper, method = "L-BFGS-B", joint = joint, rcpp = rcpp)
+  hess <- list(mle at details$hessian)
+  hess.diff <- subset(hess[[1]], rownames(hess[[1]])%in%yuima at model@parameter at diffusion, select=yuima at model@parameter at diffusion)
+  hess.drif <- subset(hess[[1]], rownames(hess[[1]])%in%yuima at model@parameter at drift, select=yuima at model@parameter at drift)
+  
+  esti <- coef(mle)
+  names(esti) <- c(yuima at model@parameter at diffusion, yuima at model@parameter at drift)
+  cic <- summary(mle)@m2logL+2*(length(yuima at model@parameter at drift)+length(yuima at model@parameter at diffusion))
+  bic <- summary(mle)@m2logL+length(yuima at model@parameter at drift)*log(Terminal)+length(yuima at model@parameter at diffusion)*log(n)
+  if(det(hess.diff) > 0 && det(hess.drif) > 0){
+    qbic <- summary(mle)@m2logL+log(det(hess.diff))+log(det(hess.drif))
+  }else{
+    qbic <- summary(mle)@m2logL+length(yuima at model@parameter at drift)*log(Terminal)+length(yuima at model@parameter at diffusion)*log(n)
+  }
+  
+  final.res <- list(par = esti, BIC = bic, QBIC = qbic, CIC = cic)
+  return(final.res)
 }
 
-
-
-