[Yuima-commits] r259 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 30 04:55:05 CET 2013
Author: kyuta
Date: 2013-10-30 04:55:05 +0100 (Wed, 30 Oct 2013)
New Revision: 259
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/cce.R
pkg/yuima/src/cce_functions.c
Log:
modify cce.R, cce_functions.c
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2013-10-28 01:47:48 UTC (rev 258)
+++ pkg/yuima/DESCRIPTION 2013-10-30 03:55:05 UTC (rev 259)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.213
-Date: 2013-10-28
+Version: 0.1.214
+Date: 2013-10-30
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2013-10-28 01:47:48 UTC (rev 258)
+++ pkg/yuima/NEWS 2013-10-30 03:55:05 UTC (rev 259)
@@ -13,4 +13,5 @@
2013/04/14: modify bns.test.Rd, mpv.Rd
2013/10/28: add cce_functions.c
modify cce.R, llag.R, sim.euler.R, bns.test.Rd, cce.Rd, llag.Rd, mpv.Rd, noisy.sampling.Rd
-2013/10/28: modify llag.R
\ No newline at end of file
+2013/10/28: modify llag.R
+2013/10/30: modify cce.R, cce_functions.c
\ No newline at end of file
Modified: pkg/yuima/R/cce.R
===================================================================
--- pkg/yuima/R/cce.R 2013-10-28 01:47:48 UTC (rev 258)
+++ pkg/yuima/R/cce.R 2013-10-30 03:55:05 UTC (rev 259)
@@ -1,2601 +1,2630 @@
-
-# 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)
-
- 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))
-
- 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 <- (26/35)*2
- b.theta <- (12/5)*(NS^2+3*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 <- max(floor(mean(theta)*Num^(1/2+delta)),2)
-
- 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
- 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(kn[i,j],sapply(ser.X,"length")) # kn must be less than the numbers of the observations
- weight <- sapply((1:(knij-1))/knij,g)
- #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
- psi.kn <- sum(weight)
-
- # 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",
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 259
More information about the Yuima-commits
mailing list