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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Apr 12 09:28:21 CEST 2018


Author: kyuta
Date: 2018-04-12 09:28:21 +0200 (Thu, 12 Apr 2018)
New Revision: 639

Modified:
   pkg/yuima/DESCRIPTION
   pkg/yuima/NAMESPACE
   pkg/yuima/NEWS
   pkg/yuima/R/cce.R
   pkg/yuima/man/cce.Rd
Log:
a bug in cce is fixed

Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION	2018-04-03 12:20:50 UTC (rev 638)
+++ pkg/yuima/DESCRIPTION	2018-04-12 07:28:21 UTC (rev 639)
@@ -1,7 +1,7 @@
 Package: yuima
 Type: Package
 Title: The YUIMA Project Package for SDEs
-Version: 1.7.7
+Version: 1.7.8
 Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
 Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2)
 Author: YUIMA Project Team

Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE	2018-04-03 12:20:50 UTC (rev 638)
+++ pkg/yuima/NAMESPACE	2018-04-12 07:28:21 UTC (rev 639)
@@ -61,6 +61,7 @@
 importFrom(stats, start)
 importFrom(utils, str)
 importFrom(stats, sd)
+importFrom("stats", cov2cor) # added by YK on Apr. 12, 2018
 
 
 exportClasses("yuima",

Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS	2018-04-03 12:20:50 UTC (rev 638)
+++ pkg/yuima/NEWS	2018-04-12 07:28:21 UTC (rev 639)
@@ -55,4 +55,5 @@
             modified llag.R, bns.test.Rd, llag.Rd, mllag.Rd, hyavar.Rd, cce.Rd, cce_functions.c
 2018/01/16: a bug in computation of asymptotic variances in llag is fixed
             the default value of tol of llag is changed to 1e-7
-            modified llag.R, llag.Rd, cce_functions.c 
\ No newline at end of file
+            modified llag.R, llag.Rd, cce_functions.c 
+2018/04/12: a bug in cce is fixed
\ No newline at end of file

Modified: pkg/yuima/R/cce.R
===================================================================
--- pkg/yuima/R/cce.R	2018-04-03 12:20:50 UTC (rev 638)
+++ pkg/yuima/R/cce.R	2018-04-12 07:28:21 UTC (rev 639)
@@ -1,2718 +1,2731 @@
-
-# cce() Cumulative Covariance Estimator
-
-#
-# CumulativeCovarianceEstimator
-#
-
-# returns a matrix of var-cov estimates
-
-########################################################
-#         functions for the synchronization
-########################################################
-
-# refresh time
-## data: a list of zoos
-
-refresh_sampling <- function(data){
-  
-  d.size <- length(data)
-  
-  if(d.size>1){
-    
-    #ser.times <- vector(d.size, mode="list")
-    #ser.times <- lapply(data,time)
-    ser.times <- lapply(lapply(data,"time"),"as.numeric")
-    ser.lengths <- sapply(data,"length")
-    ser.samplings <- vector(d.size, mode="list")
-    #refresh.times <- c()
-    
-    if(0){
-      tmp <- sapply(ser.times,"[",1)    
-      refresh.times[1] <- max(tmp)
-      
-      I <- rep(1,d.size)
-      
-      for(d in 1:d.size){
-        #ser.samplings[[d]][1] <- max(which(ser.times[[d]]<=refresh.times[1]))
-        while(!(ser.times[[d]][I[d]+1]>refresh.times[1])){
-          I[d] <- I[d]+1
-          if((I[d]+1)>ser.lengths[d]){
-            break
-          }
-        }
-        ser.samplings[[d]][1] <- I[d]
-      }
-      
-      i <- 1
-      
-      #while(all(sapply(ser.samplings,"[",i)<ser.lengths)){
-      while(all(I<ser.lengths)){
-        
-        J <- I
-        tmp <- double(d.size)
-        
-        for(d in 1:d.size){
-          #tmp[d] <- ser.times[[d]][min(which(ser.times[[d]]>refresh.times[i]))
-          repeat{
-            J[d] <- J[d] + 1
-            tmp[d] <- ser.times[[d]][J[d]]
-            if((!(J[d]<ser.lengths))||(tmp[d]>refresh.times[i])){
-              break
-            }
-          }
-        }
-        
-        i <- i+1
-        
-        refresh.times[i] <- max(tmp)
-        
-        for(d in 1:d.size){
-          #ser.samplings[[d]][i] <- max(which(ser.times[[d]]<=refresh.times[i]))
-          while(!(ser.times[[d]][I[d]+1]>refresh.times[i])){
-            I[d] <- I[d]+1
-            if((I[d]+1)>ser.lengths[d]){
-              break
-            }
-          }
-          ser.samplings[[d]][i] <- I[d]
-        }
-        
-      }
-    }
-    
-    MinL <- min(ser.lengths)
-    
-    refresh.times <- double(MinL)
-    refresh.times[1] <- max(sapply(ser.times,"[",1))
-    
-    #idx <- matrix(.C("refreshsampling",
-    #                 as.integer(d.size),
-    #                 integer(d.size),
-    #                 as.double(unlist(ser.times)),
-    #                 as.double(refresh.times),
-    #                 as.integer(append(ser.lengths,0,after=0)),
-    #                 min(sapply(ser.times,FUN="tail",n=1)),
-    #                 as.integer(MinL),
-    #                 result=integer(d.size*MinL))$result,ncol=d.size)
-    idx <- matrix(.C("refreshsampling",
-                     as.integer(d.size),
-                     integer(d.size),
-                     as.double(unlist(ser.times)),
-                     as.double(refresh.times),
-                     as.integer(ser.lengths),
-                     as.integer(diffinv(ser.lengths[-d.size],xi=0)),
-                     min(sapply(ser.times,FUN="tail",n=1)),
-                     as.integer(MinL),
-                     result=integer(d.size*MinL))$result,ncol=d.size)
-    
-    result <- vector(d.size, mode="list")
-    
-    for(d in 1:d.size){
-      #result[[d]] <- data[[d]][ser.samplings[[d]]]
-      result[[d]] <- data[[d]][idx[,d]]
-    }
-    
-    return(result)
-    
-  }else{
-    return(data)
-  }
-  
-}
-
-# refresh sampling for PHY
-## data: a list of zoos
-
-refresh_sampling.PHY <- function(data){
-  
-  d.size <- length(data)
-  
-  if(d.size>1){
-    
-    ser.times <- lapply(lapply(data,"time"),"as.numeric")
-    ser.lengths <- sapply(data,"length")    
-    #refresh.times <- max(sapply(ser.times,"[",1))
-    ser.samplings <- vector(d.size,mode="list")
-    
-    if(0){
-      for(d in 1:d.size){
-        ser.samplings[[d]][1] <- 1
-      }
-      
-      I <- rep(1,d.size)
-      i <- 1
-      
-      while(all(I<ser.lengths)){
-        
-        flag <- integer(d.size)
-        
-        for(d in 1:d.size){
-          while(I[d]<ser.lengths[d]){
-            I[d] <- I[d]+1
-            if(ser.times[[d]][I[d]]>refresh.times[i]){
-              flag[d] <- 1
-              ser.samplings[[d]][i+1] <- I[d]
-              break
-            }
-          }
-        }
-        
-        if(any(flag<rep(1,d.size))){
-          break
-        }else{
-          i <- i+1
-          tmp <- double(d.size)
-          for(d in 1:d.size){
-            tmp[d] <- ser.times[[d]][ser.samplings[[d]][i]]
-          }
-          refresh.times <- append(refresh.times,max(tmp))
-        }
-        
-      }
-    }
-    
-    MinL <- min(ser.lengths)
-    
-    refresh.times <- double(MinL)
-    refresh.times[1] <- max(sapply(ser.times,"[",1))
-    
-    #obj <- .C("refreshsamplingphy",
-    #          as.integer(d.size),
-    #          integer(d.size),
-    #          as.double(unlist(ser.times)),
-    #          rtimes=as.double(refresh.times),
-    #          as.integer(append(ser.lengths,0,after=0)),
-    #          min(sapply(ser.times,FUN="tail",n=1)),
-    #          as.integer(MinL),
-    #          Samplings=integer(d.size*(MinL+1)),
-    #          rNum=integer(1))
-    obj <- .C("refreshsamplingphy",
-              as.integer(d.size),
-              integer(d.size),
-              as.double(unlist(ser.times)),
-              rtimes=as.double(refresh.times),
-              as.integer(ser.lengths),
-              as.integer(diffinv(ser.lengths[-d.size],xi=0)),
-              min(sapply(ser.times,FUN="tail",n=1)),
-              as.integer(MinL),
-              Samplings=integer(d.size*(MinL+1)),
-              rNum=integer(1))
-    
-    refresh.times <- obj$rtimes[1:obj$rNum]
-    idx <- matrix(obj$Samplings,ncol=d.size)
-    
-    result <- vector(d.size, mode="list")
-    
-    for(d in 1:d.size){
-      #result[[d]] <- data[[d]][ser.samplings[[d]]]
-      result[[d]] <- data[[d]][idx[,d]]
-    }
-    
-    return(list(data=result,refresh.times=refresh.times))
-    
-  }else{
-    return(list(data=data,refresh.times=as.numeric(time(data[[1]]))))
-  }
-  
-}
-
-# Bibinger's synchronization algorithm
-## x: a zoo  y: a zoo
-
-Bibsynchro <- function(x,y){
-  
-  xtime <- as.numeric(time(x))
-  ytime <- as.numeric(time(y))
-  
-  xlength <- length(xtime)
-  ylength <- length(ytime)
-  N.max <- max(xlength,ylength)
-  
-  if(0){
-    mu <- integer(N.max)
-    w <- integer(N.max)
-    q <- integer(N.max)
-    r <- integer(N.max)
-    
-    if(xtime[1]<ytime[1]){
-      I <- 1
-      while(xtime[I]<ytime[1]){
-        I <- I+1
-        if(!(I<xlength)){
-          break
-        }
-      }
-      #mu0 <- min(which(ytime[1]<=xtime))
-      mu0 <- I
-      if(ytime[1]<xtime[mu0]){
-        #q[1] <- mu0-1
-        q[1] <- mu0
-      }else{
-        #q[1] <- mu0
-        q[1] <- mu0+1
-      }
-      r[1] <- 2
-    }else if(xtime[1]>ytime[1]){
-      I <- 1
-      while(xtime[I]<ytime[1]){
-        I <- I+1
-        if(!(I<xlength)){
-          break
-        }
-      }
-      #w0 <- min(which(xtime[1]<=ytime))
-      w0 <- I
-      q[1] <- 2
-      if(xtime[1]<ytime[w0]){
-        #r[1] <- w0-1
-        r[1] <- w0
-      }else{
-        #r[1] <- w0
-        r[1] <- w0+1
-      }
-    }else{
-      q[1] <- 2
-      r[1] <- 2
-    }
-    
-    i <- 1
-    
-    repeat{
-      #while((q[i]<xlength)&&(r[i]<ylength)){
-      if(xtime[q[i]]<ytime[r[i]]){
-        #tmp <- which(ytime[r[i]]<=xtime)
-        #if(identical(all.equal(tmp,integer(0)),TRUE)){
-        #  break
-        #}
-        #mu[i] <- min(tmp)
-        if(xtime[xlength]<ytime[r[i]]){
-          break
-        }
-        I <- q[i]
-        while(xtime[I]<ytime[r[i]]){
-          I <- I+1
-        }
-        mu[i] <- I
-        w[i] <- r[i]
-        if(ytime[r[i]]<xtime[mu[i]]){
-          #q[i+1] <- mu[i]-1
-          q[i+1] <- mu[i]
-        }else{
-          #q[i+1] <- mu[i]
-          q[i+1] <- mu[i]+1
-        }
-        r[i+1] <- r[i]+1
-      }else if(xtime[q[i]]>ytime[r[i]]){
-        #tmp <- which(xtime[q[i]]<=ytime)
-        #if(identical(all.equal(tmp,integer(0)),TRUE)){
-        #  break
-        #}
-        if(xtime[q[i]]>ytime[ylength]){
-          break
-        }
-        mu[i] <- q[i]
-        #w[i] <- min(tmp)
-        I <- r[i]
-        while(xtime[q[i]]>ytime[I]){
-          I <- I+1
-        }
-        w[i] <- I
-        q[i+1] <- q[i]+1
-        if(xtime[q[i]]<ytime[w[i]]){
-          #r[i+1] <- w[i]-1
-          r[i+1] <- w[i]
-        }else{
-          #r[i+1] <- w[i]
-          r[i+1] <- w[i]+1
-        }
-      }else{
-        mu[i] <- q[i]
-        w[i] <- r[i]
-        q[i+1] <- q[i]+1
-        r[i+1] <- r[i]+1
-      }
-      
-      i <- i+1
-      
-      if((q[i]>=xlength)||(r[i]>=ylength)){
-        break
-      }
-      
-    }
-    
-    mu[i] <- q[i]
-    w[i] <- r[i]
-  }
-  
-  sdata <- .C("bibsynchro",
-              as.double(xtime),
-              as.double(ytime),
-              as.integer(xlength),
-              as.integer(ylength),
-              mu=integer(N.max),
-              w=integer(N.max),
-              q=integer(N.max),
-              r=integer(N.max),
-              Num=integer(1))
-  
-  Num <- sdata$Num
-  
-  #return(list(xg=as.numeric(x)[mu[1:i]],
-  #            xl=as.numeric(x)[q[1:i]-1],
-  #            ygamma=as.numeric(y)[w[1:i]],
-  #            ylambda=as.numeric(y)[r[1:i]-1],
-  #            num.data=i))
-  return(list(xg=as.numeric(x)[sdata$mu[1:Num]+1],
-              xl=as.numeric(x)[sdata$q[1:Num]],
-              ygamma=as.numeric(y)[sdata$w[1:Num]+1],
-              ylambda=as.numeric(y)[sdata$r[1:Num]],
-              num.data=Num))
-}
-
-
-##############################################################
-#         functions for tuning parameter
-##############################################################
-
-# Barndorff-Nielsen et al. (2009)
-
-RV.sparse <- function(zdata,frequency=1200,utime){
-  
-  znum <- as.numeric(zdata)
-  ztime <- as.numeric(time(zdata))*utime
-  n.size <- length(zdata)
-  end <- ztime[n.size]
-  
-  grid <- seq(ztime[1],end,by=frequency)
-  n.sparse <- length(grid)
-  
-  #result <- double(frequency)
-  #result <- 0
-  #I <- rep(1,n.sparse)
-  
-  #for(t in 1:frequency){
-  #  for(i in 1:n.sparse){
-  #    while((ztime[I[i]+1]<=grid[i])&&(I[i]<n.size)){
-  #      I[i] <- I[i]+1
-  #    }
-  #  }
-  #  result[t] <- sum(diff(znum[I])^2)
-    #result <- result+sum(diff(znum[I])^2)
-  #  grid <- grid+rep(1,n.sparse)
-  #  if(grid[n.sparse]>end){
-  #    grid <- grid[-n.sparse]
-  #    I <- I[-n.sparse]
-  #    n.sparse <- n.sparse-1
-  #  }
-  #}
-  
-  K <- floor(end-grid[n.sparse]) + 1
-  
-  zmat <- matrix(.C("ctsubsampling",
-                    as.double(znum),
-                    as.double(ztime),
-                    as.integer(frequency),
-                    as.integer(n.sparse),
-                    as.integer(n.size),
-                    as.double(grid),
-                    result=double(frequency*n.sparse))$result,
-                 n.sparse,frequency)
-  
-  result <- double(frequency)
-  result[1:K] <- colSums(diff(zmat[,1:K])^2)
-  result[-(1:K)] <- colSums(diff(zmat[-n.sparse,-(1:K)])^2)
-  
-  return(mean(result))
-  #return(result/frequency)
-  #return(znum[I])
-}
-
-Omega_BNHLS <- function(zdata,sec=120,utime){
-  
-  #q <- ceiling(sec/mean(diff(as.numeric(time(zdata))*utime)))
-  #q <- ceiling(sec*(length(zdata)-1)/utime)
-  ztime <- as.numeric(time(zdata))
-  q <- ceiling(sec*(length(zdata)-1)/(utime*(tail(ztime,n=1)-head(ztime,n=1))))
-  obj <- diff(as.numeric(zdata),lag=q)
-  n <- length(obj)
-  
-  #result <- 0
-  result <- double(q)
-  
-  for(i in 1:q){
-    tmp <- obj[seq(i,n,by=q)]
-    n.tmp <- sum(tmp!=0)
-    if(n.tmp>0){
-      result[i] <- sum(tmp^2)/(2*n.tmp)
-    }
-  }
-  
-  #return(result/q)
-  return(mean(result))
-}
-
-NSratio_BNHLS <- function(zdata,frequency=1200,sec=120,utime){
-  
-  IV <- RV.sparse(zdata,frequency,utime)
-  Omega <- Omega_BNHLS(zdata,sec,utime)
-  
-  return(Omega/IV)
-}
-
-# Pre-averaging estimator
-
-selectParam.pavg <- function(data,utime,frequency=1200,sec=120,
-                             a.theta,b.theta,c.theta){
-  
-  if(missing(utime)) utime <- ifelse(is.numeric(time(data[[1]])),23400,1)
-  
-  NS <- sapply(data,FUN=NSratio_BNHLS,frequency=frequency,sec=sec,utime=utime)
-  coef <- (b.theta+sqrt(b.theta^2+3*a.theta*c.theta))/a.theta
-  
-  return(sqrt(coef*NS))
-}
-
-selectParam.TSCV <- function(data,utime,frequency=1200,sec=120){
-  
-  if(missing(utime)) utime <- ifelse(is.numeric(time(data[[1]])),23400,1)
-  
-  NS <- sapply(data,FUN=NSratio_BNHLS,frequency=frequency,sec=sec,
-               utime=utime)
-  
-  return((12*NS)^(2/3))
-}
-
-selectParam.GME <- function(data,utime,frequency=1200,sec=120){
-  
-  if(missing(utime)) utime <- ifelse(is.numeric(time(data[[1]])),23400,1)
-  
-  NS <- sapply(data,FUN=NSratio_BNHLS,frequency=frequency,sec=sec,
-               utime=utime)
-  
-  a.theta <- 52/35
-  #b.theta <- (12/5)*(NS^2+3*NS)
-  b.theta <- (24/5)*(NS^2+2*NS)
-  c.theta <- 48*NS^2
-  
-  return(sqrt((b.theta+sqrt(b.theta^2+12*a.theta*c.theta))/(2*a.theta)))
-}
-
-selectParam.RK <- function(data,utime,frequency=1200,sec=120,cstar=3.5134){
-  
-  if(missing(utime)) utime <- ifelse(is.numeric(time(data[[1]])),23400,1)
-  
-  NS <- sapply(data,FUN=NSratio_BNHLS,frequency=frequency,sec=sec,
-               utime=utime)
-  
-  return(cstar*(NS^(2/5)))
-}
-
-if(0){
-selectParam_BNHLS <- function(yuima,method,utime=1,frequency=1200,sec=120,kappa=7585/1161216,
-                              kappabar=151/20160,kappatilde=1/24,cstar=3.51){
-  data <- get.zoo.data(yuima)
-  
-  switch(method,
-         "PHY"=selectParam.PHY(data,utime,frequency,sec,kappa,kappabar,kappatilde),
-         "GME"=selectParam.GME(data,utime,frequency,sec),
-         "RK"=selectParam.RK(data,utime,frequency,sec,cstar),
-         "PTHY"=selectParam.PHY(data,utime,frequency,sec,kappa,kappabar,kappatilde))
-  
-}
-}
-
-###############################################################
-#      functions for selecting the threshold functions
-###############################################################
-
-# local universal thresholding method
-
-## Bipower variation
-### x: a zoo lag: a positive integer
-
-BPV <- function(x,lag=1){
-  
-  n <- length(x)-1
-  dt <- diff(as.numeric(time(x)))
-  obj <- abs(diff(as.numeric(x)))/sqrt(dt)
-  
-  result <- (pi/2)*(obj[1:(n-lag)]*dt[1:(n-lag)])%*%obj[(1+lag):n]
-  
-  return(result)
-  
-}
-
-## local univaersal thresholding method
-### data: a list of zoos  coef: a positive number
-
-local_univ_threshold <- function(data,coef=5,eps=0.1){
-  
-  d.size <- length(data)
-  
-  result <- vector(d.size,mode="list") 
-  
-  for(d in 1:d.size){
-    
-    x <- data[[d]]
-    #x <- as.numeric(data[[d]])
-    n <- length(x)-1
-    dt <- diff(as.numeric(time(x)))
-    obj <- abs(diff(as.numeric(x)))/sqrt(dt)
-    #dx <- diff(as.numeric(x))
-    
-    #xtime <- time(x)
-    #xtime <- time(data[[d]])
-    #if(is.numeric(xtime)){
-    #  r <- max(diff(xtime))
-    #}else{
-    #  xtime <- as.numeric(xtime)
-      #r <- max(diff(xtime))/(length(x)-1)
-    #  r <- max(diff(xtime))/(tail(xtime,n=1)-head(xtime,n=1))
-    #}
-    #K <- ceiling(sqrt(1/r))
-    #K <- max(ceiling(sqrt(1/r)),3)
-    K <- max(ceiling(n^0.5),3)
-    
-    coef2 <- coef^2
-    rho <- double(n+1)
-    
-    #tmp <- coef2*sum(dx^2)*n^(eps-1)
-    #tmp <- double(n+1)
-    #tmp[-(1:(K-1))] <- coef2*n^eps*rollmean(dx^2,k=K-1,align="left")
-    #tmp[1:(K-1)] <- tmp[K]
-    #dx[dx^2>tmp[-(n+1)]] <- 0
-    
-    if(K<n){
-      #rho[1:K] <- coef2*(mad(diff(as.numeric(x)[1:(K+1)]))/0.6745)^2
-      #rho[1:K] <- coef2*2*log(K)*(mad(diff(x[1:(K+1)]))/0.6745)^2
-      #for(i in (K+1):n){
-      #  rho[i] <- coef2*2*log(length(x))*BPV(x[(i-K):(i-1)])/(K-2)
-      #}
-      #rho[-(1:K)] <- coef2*2*log(n)^(1+eps)*
-      # #rollapply(x[-n],width=K,FUN=BPV,align="left")/(K-2)
-      rho[-(1:(K-1))] <- coef2*n^eps*(pi/2)*
-        rollmean((obj*dt)[-n]*obj[-1],k=K-2,align="left")
-      #rho[-(1:(K-1))] <- coef2*n^eps*rollmean(dx^2,k=K-1,align="left")
-      rho[1:(K-1)] <- rho[K]
-    }else{
-      #rho <- coef2*(mad(diff(as.numeric(x)))/0.6745)^2
-      #rho <- coef2*(mad(diff(x))/0.6745)^2
-      #rho <- coef2*2*log(n)^(1+eps)*BPV(x)
-      rho <- coef2*n^(eps-1)*BPV(x)
-      #rho <- coef2*sum(dx^2)*n^(eps-1)
-    }
-    
-    result[[d]] <- rho[-(n+1)]
-    
-  }
-  
-  return(result)
-  
-}
-
-#################################################################
-#       functions for calculating the estimators
-#################################################################
-
-# Hayashi-Yoshida estimator
-## data: a list of zoos
-
-HY <- function(data) {  
-  
-  n.series <- length(data)
-  
-  # allocate memory
-  ser.X <- vector(n.series, mode="list")     # data in 'x'
-  ser.times <- vector(n.series, mode="list") # time index in 'x'
-  ser.diffX <- vector(n.series, mode="list") # difference of data
-  
-  for(i in 1:n.series){
-    # set data and time index
-    ser.X[[i]] <- as.numeric(data[[i]]) # we need to transform data into numeric to avoid problem with duplicated indexes below
-    ser.times[[i]] <- as.numeric(time(data[[i]]))
-    # set difference of the data 
-    ser.diffX[[i]] <- diff( ser.X[[i]] )
-  }
-  
-  # core part of cce
-  
-  cmat <- matrix(0, n.series, n.series)  # cov
-  for(i in 1:n.series){
-    for(j in i:n.series){ 
-      #I <- rep(1,n.series)
-      #Checking Starting Point
-      #repeat{
-      #while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))){
-      #  if(ser.times[[i]][I[i]] >= ser.times[[j]][I[j]+1]){
-      #    I[j] <- I[j]+1   
-      #  }else if(ser.times[[i]][I[i]+1] <= ser.times[[j]][I[j]]){
-      #    I[i] <- I[i]+1   
-      #  }else{
-      #    break
-      #  }
-      #}
-      
-      
-      #Main Component
-      if(i!=j){
-        #while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))) {
-        #  cmat[j,i] <- cmat[j,i] + (ser.diffX[[i]])[I[i]]*(ser.diffX[[j]])[I[j]]
-        #  if(ser.times[[i]][I[i]+1]>ser.times[[j]][I[j]+1]){
-        #    I[j] <- I[j] + 1
-        #  }else if(ser.times[[i]][I[i]+1]<ser.times[[j]][I[j]+1]){
-        #    I[i] <- I[i] +1
-        #  }else{
-        #    I[i] <- I[i]+1
-        #    I[j] <- I[j]+1
-        #  }
-        #}
-        cmat[j,i] <- .C("HayashiYoshida",as.integer(length(ser.times[[i]])),
-                        as.integer(length(ser.times[[j]])),as.double(ser.times[[i]]),
-                        as.double(ser.times[[j]]),as.double(ser.diffX[[i]]),
-                        as.double(ser.diffX[[j]]),value=double(1))$value
-      }else{
-        cmat[i,j] <- sum(ser.diffX[[i]]^2)
-      }
-      cmat[i,j] <- cmat[j,i]  
-    }
-    
-  }
-  
-  return(cmat)
-}
-
-
-#########################################################
-
-# Modulated realized covariance (Cristensen et al.(2010))
-## data: a list of zoos  theta: a positive number
-## g: a real-valued function on [0,1] (a weight function)
-## delta: a positive number
-
-MRC <- function(data,theta,kn,g,delta=0,adj=TRUE){
-  
-  n.series <- length(data)
-  
-  #if(missing(theta)&&missing(kn))
-  #  theta <- selectParam.pavg(data,utime=utime,a.theta=151/80640,
-  #                            b.theta=1/96,c.theta=1/6)
-  if(missing(theta)) theta <- 1
-  
-  cmat <- matrix(0, n.series, n.series)  # cov
-  
-  # synchronize the data and take the differences
-  #diffX <- lapply(lapply(refresh_sampling(data),"as.numeric"),"diff")
-  #diffX <- do.call("rbind",diffX) # transform to matrix
-  diffX <- diff(do.call("cbind",lapply(refresh_sampling(data),"as.numeric")))
-  
-  #Num <- ncol(diffX)
-  Num <- nrow(diffX)
-  
-  if(missing(kn)) kn <- floor(mean(theta)*Num^(1/2+delta))
-  
-  kn <- min(max(kn,2),Num+1)
-  
-  weight <- sapply((1:(kn-1))/kn,g)
-  psi2.kn <- sum(weight^2)
-  
-  # pre-averaging
-  #myfunc <- function(dx)rollapplyr(dx,width=kn-1,FUN="%*%",weight)
-  #barX <- apply(diffX,1,FUN=myfunc)
-  barX <- filter(diffX,filter=rev(weight),method="c",sides=1)[(kn-1):Num,]
-  
-  cmat <- (Num/(Num-kn+2))*t(barX)%*%barX/psi2.kn
-  
-  if(delta==0){
-    psi1.kn <- kn^2*sum(diff(c(0,weight,0))^2)
-    scale <- psi1.kn/(theta^2*psi2.kn*2*Num)
-    #cmat <- cmat-scale*diffX%*%t(diffX)
-    cmat <- cmat-scale*t(diffX)%*%diffX
-    if(adj) cmat <- cmat/(1-scale)
-  }
-  
-  return(cmat)
-}
-
-#########################################################
-
-# Pre-averaged Hayashi-Yoshida estimator
-## data: a list of zoos  theta: a positive number
-## g: a real-valued function on [0,1] (a weight function)
-## refreshing: a logical value (if TRUE we use the refreshed data)
-
-PHY <- function(data,theta,kn,g,refreshing=TRUE,cwise=TRUE){
-  
-  n.series <- length(data)
-  
-  #if(missing(theta)&&missing(kn))
-  #  theta <- selectParam.pavg(data,a.theta=7585/1161216,
-  #                            b.theta=151/20160,c.theta=1/24)
-  if(missing(theta)) theta <- 0.15
-    
-  cmat <- matrix(0, n.series, n.series)  # cov
-  
-  if(refreshing){
-    
-    if(cwise){
-      
-      if(missing(kn)){# refreshing,cwise,missing kn
-        
-        theta <- matrix(theta,n.series,n.series) 
-        theta <- (theta+t(theta))/2
-        
-        for(i in 1:n.series){
-          for(j in i:n.series){ 
-            if(i!=j){
-              
-              resample <- refresh_sampling.PHY(list(data[[i]],data[[j]]))
-              dat <- resample$data
-              ser.numX <- length(resample$refresh.times)-1
-              
-              # set data and time index
-              ser.X <- lapply(dat,"as.numeric")
-              ser.times <- lapply(lapply(dat,"time"),"as.numeric")
-              
-              # set difference of the data
-              ser.diffX <- lapply(ser.X,"diff")
-              
-              # pre-averaging
-              kn <- min(max(ceiling(theta[i,j]*sqrt(ser.numX)),2),ser.numX)
-              weight <- sapply((1:(kn-1))/kn,g)
-              psi.kn <- sum(weight)
-              
-              #ser.barX <- list(rollapplyr(ser.diffX[[1]],width=kn-1,FUN="%*%",weight),
-              #                 rollapplyr(ser.diffX[[2]],width=kn-1,FUN="%*%",weight))
-              ser.barX <- list(filter(ser.diffX[[1]],rev(weight),method="c",
-                                      sides=1)[(kn-1):length(ser.diffX[[1]])],
-                               filter(ser.diffX[[2]],rev(weight),method="c",
-                                      sides=1)[(kn-1):length(ser.diffX[[2]])])
-              
-              ser.num.barX <- sapply(ser.barX,"length")-1
-              
-              # core part of cce
-              #start <- kn+1
-              #end <- 1
-              #for(k in 1:ser.num.barX[1]){
-              #  while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-kn)<ser.num.barX[2])){
-              #    start <- start + 1
-              #  }
-              #  while((ser.times[[1]][k+kn]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
-              #    end <- end + 1
-              #  }
-              #  cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-kn):end])
-              #}
-              cmat[i,j] <- .C("pHayashiYoshida",
-                              as.integer(kn),
-                              as.integer(ser.num.barX[1]),
-                              as.integer(ser.num.barX[2]),
-                              as.double(ser.times[[1]]),
-                              as.double(ser.times[[2]]),
-                              as.double(ser.barX[[1]]),
-                              as.double(ser.barX[[2]]),
-                              value=double(1))$value
-              
-              cmat[i,j] <- cmat[i,j]/(psi.kn^2)
-              cmat[j,i] <- cmat[i,j]
-              
-            }else{
-              
-              diffX <- diff(as.numeric(data[[i]]))
-              
-              kn <- min(max(ceiling(theta[i,i]*sqrt(length(diffX))),2),
-                        length(data[[i]]))
-              weight <- sapply((1:(kn-1))/kn,g)
-              psi.kn <- sum(weight)
-              
-              #barX <- rollapplyr(diffX,width=kn-1,FUN="%*%",weight)
-              barX <- filter(diffX,rev(weight),method="c",
-                             sides=1)[(kn-1):length(diffX)]
-              tmp <- barX[-length(barX)]
-              #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
-              cmat[i,j] <- tmp%*%rollsum(c(double(kn-1),tmp,double(kn-1)),
-                                         k=2*(kn-1)+1)/(psi.kn)^2
-              
-            }
-          }
-        }
-      }else{# refreshing,cwise,not missing kn
-        
-        kn <- matrix(kn,n.series,n.series)
-        
-        for(i in 1:n.series){
-          for(j in i:n.series){
-            
-            if(i!=j){
-              
-              resample <- refresh_sampling.PHY(list(data[[i]],data[[j]]))
-              dat <- resample$data
-              ser.numX <- length(resample$refresh.times)-1
-              
-              # set data and time index
-              ser.X <- lapply(dat,"as.numeric")
-              ser.times <- lapply(lapply(dat,"time"),"as.numeric")
-              
-              # set difference of the data
-              ser.diffX <- lapply(ser.X,"diff")
-              
-              # pre-averaging
-              knij <- min(max(kn[i,j],2),ser.numX+1)
-              weight <- sapply((1:(knij-1))/knij,g)
-              psi.kn <- sum(weight)
-              
-              #ser.barX <- list(rollapplyr(ser.diffX[[1]],width=knij-1,FUN="%*%",weight),
-              #                 rollapplyr(ser.diffX[[2]],width=knij-1,FUN="%*%",weight))
-              ser.barX <- list(filter(ser.diffX[[1]],rev(weight),method="c",
-                                      sides=1)[(knij-1):length(ser.diffX[[1]])],
-                               filter(ser.diffX[[2]],rev(weight),method="c",
-                                      sides=1)[(knij-1):length(ser.diffX[[2]])])
-              
-              ser.num.barX <- sapply(ser.barX,"length")-1
-              
-              # core part of cce
-              #start <- knij+1
-              #end <- 1
-              #for(k in 1:ser.num.barX[1]){
-              #  while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-knij)<ser.num.barX[2])){
-              #    start <- start + 1
-              #  }
-              #  while((ser.times[[1]][k+knij]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
-              #    end <- end + 1
-              #  }
-              #  cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-knij):end])
-              #}
-              cmat[i,j] <- .C("pHayashiYoshida",
-                              as.integer(knij),
-                              as.integer(ser.num.barX[1]),
-                              as.integer(ser.num.barX[2]),
-                              as.double(ser.times[[1]]),
-                              as.double(ser.times[[2]]),
-                              as.double(ser.barX[[1]]),
-                              as.double(ser.barX[[2]]),
-                              value=double(1))$value
-              
-              cmat[i,j] <- cmat[i,j]/(psi.kn^2)
-              cmat[j,i] <- cmat[i,j]
-              
-            }else{
-              
-              diffX <- diff(as.numeric(data[[i]]))
-              # pre-averaging
-              kni <- min(max(kn[i,i],2),length(data[[i]]))
-              weight <- sapply((1:(kni-1))/kni,g)
-              psi.kn <- sum(weight)
-              
-              #barX <- rollapplyr(diffX,width=kni-1,FUN="%*%",weight)
-              barX <- filter(diffX,rev(weight),method="c",
-                             sides=1)[(kni-1):length(diffX)]
-              tmp <- barX[-length(barX)]
-              
-              # core part of cce
-              #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kni-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
-              cmat[i,j] <- tmp%*%rollsum(c(double(kni-1),tmp,double(kni-1)),
-                                         k=2*(kni-1)+1)/(psi.kn)^2
-            }
-          }
-        }
-      }
-    }else{# refreshing,non-cwise
-      
-      # synchronize the data
-      resample <- refresh_sampling.PHY(data)
-      data <- resample$data
-      ser.numX <- length(resample$refresh.times)-1
-      
-      # if missing kn, we choose it following Barndorff-Nielsen et al. (2011)
-      if(missing(kn)){
-        kn <- min(max(ceiling(mean(theta)*sqrt(ser.numX)),2),ser.numX+1)
-      }
-      kn <- kn[1]
-      
-      weight <- sapply((1:(kn-1))/kn,g)
-      
-      # allocate memory
-      ser.X <- vector(n.series, mode="list")     # data in 'x'
-      ser.times <- vector(n.series, mode="list") # time index in 'x'
-      ser.diffX <- vector(n.series, mode="list") # difference of data
-      ser.barX <- vector(n.series, mode="list")  # pre-averaged data
-      ser.num.barX <- integer(n.series)         
-      
-      for(i in 1:n.series){
-        # set data and time index
-        ser.X[[i]] <- as.numeric(data[[i]]) # we need to transform data into numeric to avoid problem with duplicated indexes below
-        ser.times[[i]] <- as.numeric(time(data[[i]]))
-        
-        # set difference of the data 
-        ser.diffX[[i]] <- diff( ser.X[[i]] )
-        
-        # pre-averaging
-        #ser.barX[[i]] <- rollapplyr(ser.diffX[[i]],width=kn-1,FUN="%*%",weight)
-        ser.barX[[i]] <- filter(ser.diffX[[i]],rev(weight),method="c",
-                                sides=1)[(kn-1):length(ser.diffX[[i]])]
-        ser.num.barX[i] <- length(ser.barX[[i]])-1
-      }
-      
-      # core part of cce
-      
-      for(i in 1:n.series){
-        for(j in i:n.series){ 
-          if(i!=j){
-            #start <- kn+1
-            #end <- 1
-            #for(k in 1:ser.num.barX[i]){
-            #  while(!(ser.times[[i]][k]<ser.times[[j]][start])&&((start-kn)<ser.num.barX[j])){
-            #    start <- start + 1
-            #  }
-            #  while((ser.times[[i]][k+kn]>ser.times[[j]][end+1])&&(end<ser.num.barX[j])){
-            #    end <- end + 1
-            #  }
-            #  cmat[i,j] <- cmat[i,j] + ser.barX[[i]][k]*sum(ser.barX[[j]][(start-kn):end])
-            #}
-            cmat[i,j] <- .C("pHayashiYoshida",
-                            as.integer(kn),
-                            as.integer(ser.num.barX[i]),
-                            as.integer(ser.num.barX[j]),
-                            as.double(ser.times[[i]]),
-                            as.double(ser.times[[j]]),
-                            as.double(ser.barX[[i]]),
-                            as.double(ser.barX[[j]]),
-                            value=double(1))$value
-            cmat[j,i] <- cmat[i,j]
-          }else{
-            tmp <- ser.barX[[i]][1:ser.num.barX[i]]
-            #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)
-            cmat[i,j] <- tmp%*%rollsum(c(double(kn-1),tmp,double(kn-1)),
-                                       k=2*(kn-1)+1) 
-          }
-        }
-      }
-      
-      psi.kn <- sum(weight)
-      cmat <- cmat/(psi.kn^2)
-      
-    }
-  }else{# non-refreshing
-    
-    if(cwise){# non-refreshing,cwise
-      
-      theta <- matrix(theta,n.series,n.series) 
-      theta <- (theta+t(theta))/2
-      
-      if(missing(kn)){
-        ntmp <- matrix(sapply(data,"length"),n.series,n.series)
-        ntmp <- ntmp+t(ntmp)
-        diag(ntmp) <- diag(ntmp)/2
-        kn <- ceiling(theta*sqrt(ntmp))
-        kn[kn<2] <- 2 # kn must be larger than 1
-      }
-      kn <- matrix(kn,n.series,n.series)
-      
-      for(i in 1:n.series){
-        for(j in i:n.series){ 
-          if(i!=j){
-            
-            dat <- list(data[[i]],data[[j]])
-            
-            # set data and time index
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/yuima -r 639


More information about the Yuima-commits mailing list