[Yuima-commits] r214 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 13 10:59:27 CET 2012
Author: kyuta
Date: 2012-12-13 10:59:27 +0100 (Thu, 13 Dec 2012)
New Revision: 214
Added:
pkg/yuima/R/noisy.sampling.R
pkg/yuima/man/noisy.sampling.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/R/bns.test.R
pkg/yuima/R/cce.R
pkg/yuima/R/llag.R
pkg/yuima/R/mpv.R
pkg/yuima/man/bns.test.Rd
pkg/yuima/man/cce.Rd
pkg/yuima/man/llag.Rd
pkg/yuima/man/mpv.Rd
Log:
add noisy.sampling.R, noisy.sampling.Rd
modify bns.test.R, cce.R, llag.R, mpv.R, bns.test.Rd, cce.Rd, llag.Rd, mpv.Rd
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2012-10-14 09:37:37 UTC (rev 213)
+++ pkg/yuima/DESCRIPTION 2012-12-13 09:59:27 UTC (rev 214)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.197
-Date: 2012-10-05
+Version: 0.1.198
+Date: 2012-12-13
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2012-10-14 09:37:37 UTC (rev 213)
+++ pkg/yuima/NAMESPACE 2012-12-13 09:59:27 UTC (rev 214)
@@ -71,6 +71,7 @@
export(cce)
export(llag)
export(poisson.random.sampling)
+export(noisy.sampling)
export(mpv)
export(bns.test)
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2012-10-14 09:37:37 UTC (rev 213)
+++ pkg/yuima/NEWS 2012-12-13 09:59:27 UTC (rev 214)
@@ -1 +1,3 @@
-2012/10/05: add mpv.R, bns.test.R, mpv.Rd, bns.test.Rd
\ No newline at end of file
+2012/10/05: add mpv.R, bns.test.R, mpv.Rd, bns.test.Rd
+2012/12/13: add noisy.sampling.R, noisy.sampling.Rd
+ modify bns.test.R, cce.R, llag.R, mpv.R, bns.test.Rd, cce.Rd, llag.Rd, mpv.Rd
\ No newline at end of file
Modified: pkg/yuima/R/bns.test.R
===================================================================
--- pkg/yuima/R/bns.test.R 2012-10-14 09:37:37 UTC (rev 213)
+++ pkg/yuima/R/bns.test.R 2012-12-13 09:59:27 UTC (rev 214)
@@ -18,7 +18,7 @@
bns.stat_standard <- function(yuima,r=rep(1,4)){
- JV <- mpv(yuima,2)-mpv(yuima)
+ JV <- mpv(yuima,r=2)-mpv(yuima,r=c(1,1))
IQ <- mpv(yuima,r)
return(JV/sqrt((pi^2/4+pi-5)*IQ))
@@ -26,8 +26,8 @@
bns.stat_ratio <- function(yuima,r=rep(1,4),adj=TRUE){
- bpv <- mpv(yuima)
- RJ <- 1-bpv/mpv(yuima,2)
+ bpv <- mpv(yuima,r=c(1,1))
+ RJ <- 1-bpv/mpv(yuima,r=2)
avar <- mpv(yuima,r)/bpv^2
if(adj){
@@ -39,8 +39,8 @@
bns.stat_log <- function(yuima,r=rep(1,4),adj=TRUE){
- bpv <- mpv(yuima)
- RJ <- log(mpv(yuima,2)/bpv)
+ bpv <- mpv(yuima,r=c(1,1))
+ RJ <- log(mpv(yuima,r=2)/bpv)
avar <- mpv(yuima,r)/bpv^2
if(adj){
Modified: pkg/yuima/R/cce.R
===================================================================
--- pkg/yuima/R/cce.R 2012-10-14 09:37:37 UTC (rev 213)
+++ pkg/yuima/R/cce.R 2012-12-13 09:59:27 UTC (rev 214)
@@ -1,89 +1,2254 @@
-# cce() Cumulative Covariance Estimator
-
-#
-# CumulativeCovarianceEstimator
-#
-
-# returns a matrix of var-cov estimates
-
-setGeneric( "cce", function(x) standardGeneric("cce") )
-
-setMethod("cce", "yuima", function(x) cce(x at data) )
-
-setMethod("cce", "yuima.data", function(x) {
-
- data <- get.zoo.data(x)
- n.series <- length(data)
-#if(n.series <2)
-# stop("Please provide at least 2-dimensional time series")
-
-# 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]]))
-
-# NA data must be skipped
- idt <- which(is.na(ser.X[[i]]))
- if(length(idt>0)){
- ser.X[[i]] <- (ser.X[[i]])[-idt]
- ser.times[[i]] <- (ser.times[[i]])[-idt]
- }
- if(length(ser.X[[i]])<2) {
- stop("length of data (w/o NA) must be more than 1")
- }
-# 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{
- 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
- }
- }
- }else{
- cmat[i,j] <- sum(ser.diffX[[i]]^2)
- }
- cmat[i,j] <- cmat[j,i]
- }
-
- }
- if(n.series >1){
- sdmat <- diag(sqrt(diag(cmat)))
- cormat <- solve(sdmat) %*% cmat %*% solve(sdmat)
- }else{
- cormat <- as.matrix(1)
- }
- return( list(covmat=cmat,cormat=cormat) )
-})
+
+# 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()
+
+ 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]
+ }
+
+ }
+
+ result <- vector(d.size, mode="list")
+
+ for(d in 1:d.size){
+ result[[d]] <- data[[d]][ser.samplings[[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")
+
+ 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))
+ }
+
+ }
+
+ result <- vector(d.size, mode="list")
+
+ for(d in 1:d.size){
+ result[[d]] <- data[[d]][ser.samplings[[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)
+
+ mu <- c()
+ w <- c()
+ q <- c()
+ r <- c()
+
+ 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
+ }else{
+ q[1] <- mu0
+ }
+ r[1] <- 1
+ }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] <- 1
+ if(xtime[1]<ytime[w0]){
+ r[1] <- w0-1
+ }else{
+ r[1] <- w0
+ }
+ }else{
+ q[1] <- 1
+ r[1] <- 1
+ }
+
+ 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
+ }else{
+ q[i+1] <- mu[i]
+ }
+ 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
+ }else{
+ r[i+1] <- w[i]
+ }
+ }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
+ #}
+
+ }
+
+ return(list(xg=as.numeric(x)[mu[1:(i-1)]],
+ xl=as.numeric(x)[q[1:(i-1)]],
+ ygamma=as.numeric(y)[w[1:(i-1)]],
+ ylambda=as.numeric(y)[r[1:(i-1)]],
+ num.data=i-1))
+
+}
+
+
+##############################################################
+# 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
+ }
+ }
+
+ return(mean(result))
+ #return(result/frequency)
+}
+
+Omega_BNHLS <- function(zdata,sec=120,utime){
+
+ #q <- ceiling(sec/mean(diff(as.numeric(time(zdata))*utime)))
+ q <- ceiling(sec*(length(zdata)-1)/utime)
+ 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+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
+ obj <- abs(diff(as.numeric(x)))
+
+ result <- (pi/2)*obj[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=3){
+
+ 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)
+
+ #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(sqrt(n)),3)
+
+ coef2 <- coef^2
+ rho <- double(n)
+
+ 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)*
+ rollapply(x[-n],width=K,FUN="BPV",align="left")/(K-2)
+ }else{
+ #rho <- coef2*(mad(diff(as.numeric(x)))/0.6745)^2
+ rho <- coef2*(mad(diff(x))/0.6745)^2
+ }
+
+ result[[d]] <- rho[-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{
+ 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
+ }
+ }
+ }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,utime){
+
+ 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)
+
+ 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
+
+ Num <- ncol(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)
+
+ 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)
+ 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,utime){
+
+ n.series <- length(data)
+
+ if(missing(theta)&&missing(kn))
+ theta <- selectParam.pavg(data,utime=utime,a.theta=7585/1161216,
+ b.theta=151/20160,c.theta=1/24)
+
+ 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.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] <- 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)
+ tmp <- barX[-length(barX)]
+ cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)/(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.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] <- 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)
+ 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
+
+ }
+ }
+ }
+ }
+ }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.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[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)
+ }
+ }
+ }
+
+ 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.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] <- cmat[i,j]/(psi.kn^2)
+ cmat[j,i] <- cmat[i,j]
+
+ }else{
+
+ diffX <- diff(as.numeric(data[[i]]))
+
+ # pre-averaging
+ kni <- min(kn[i,i],length(data[[i]])) # kn must be less than the number of the observations
+ weight <- sapply((1:(kni-1))/kni,g)
+ psi.kn <- sum(weight)
+
+ barX <- rollapplyr(diffX,width=kni-1,FUN="%*%",weight)
+ tmp <- barX[-length(barX)]
+ cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kni-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
+
+ }
+ }
+ }
+ }else{# non-refreshing, non-cwise
+
+ # 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.numX <- integer(n.series)
+ ser.barX <- vector(n.series, mode="list")
+
+ 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]] )
+ ser.numX[i] <- length(ser.diffX[[i]])
+ }
+
+ # pre-averaging
+ if(missing(kn)){
+ kn <- min(max(ceiling(mean(theta)*sqrt(sum(ser.numX))),2),
+ ser.numX+1)
+ }
+ kn <- kn[1]
+
+ weight <- sapply((1:(kn-1))/kn,g)
+ psi.kn <- sum(weight)
+
+ for(i in 1:n.series){
+ ser.barX[[i]] <- rollapplyr(ser.diffX[[i]],width=kn-1,FUN="%*%",weight)
+ }
+
+ ser.num.barX <- sapply(ser.barX,"length")-1
+
+ # core part of cce
+
+ cmat <- matrix(0, n.series, n.series) # cov
+ 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[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 <- cmat/(psi.kn^2)
+
+ }
+ }
+
+ return(cmat)
+}
+
+
+################################################################
+
+# previous tick Two Scales realized CoVariance
+## data: a list of zoos c.two: a postive number
+
+TSCV <- function(data,K,c.two,J=1,adj=TRUE,utime){
+
+ X <- do.call("rbind",lapply(refresh_sampling(data),"as.numeric"))
+ N <- ncol(X)
+
+ if(missing(K)){
+ if(missing(c.two)) c.two <- selectParam.TSCV(data,utime=utime)
+ K <- ceiling(mean(c.two)*N^(2/3))
+ }
+
+ scale <- (N-K+1)*J/(K*(N-J+1))
+
+ diffX1 <- apply(X,1,"diff",lag=K)
+ diffX2 <- apply(X,1,"diff",lag=J)
+
+ cmat <- t(diffX1)%*%diffX1/K-scale*t(diffX2)%*%diffX2/J
+ if(adj) cmat <- cmat/(1-scale)
+
+ return(cmat)
+}
+
+################################################################
+
+# Generalized multiscale estimator
+## data: a list of zoos c.multi: a postive number
+
+GME <- function(data,c.multi,utime){
+
+ d.size <- length(data)
+
+ if(missing(c.multi)) c.multi <- selectParam.GME(data,utime=utime)
+
+ c.multi <- matrix(c.multi,d.size,d.size)
+ c.multi <- (c.multi+t(c.multi))/2
+
+ cmat <- matrix(0,d.size,d.size)
+ for(i in 1:d.size){
+ for(j in i:d.size){
+
+ if(i!=j){
+
+ sdata <- Bibsynchro(data[[i]],data[[j]])
+
+ N <- sdata$num.data
+
+ M <- ceiling(c.multi[i,j]*sqrt(N))
+ M <- min(c(M,N)) # M must be smaller than N
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 214
More information about the Yuima-commits
mailing list