Wavelet variances

[stephaneguerrier]

Hi guys,

Once Issue #4 has been addressed the next part is to compute various statistics from the wavelet decomposition, here are the main ones:

• Classical wavelet variance: this is the classical estimator proposed by Percival. This is already implemented in the gmwm package. This should be straightforward.
• Robust wavelet variance: this robust estimator @robertomolinari and I proposed. This is also implemented on the gmwm package. However, it might be good to double check that and in particular the covariance matrix. @robertomolinari: could you add some details here?
• Spatial (robust) wavelet variance: This should identical to the first two points. The only difference is that wavelet variance should be computed on a vectorized version of the spatial wavelet decomposition. If the first two points are addressed this should be simple. @robertomolinari : could you add some details here (if needed!) Thanks!
• Cross-covariance: This is the covariance between two wavelet transforms. @HaotianXu: could you please add some details here? Thanks!

Note that for all the quantities described above we will need to compute the point estimate together with its (estimated) variance.

[HaotianXu]

Regarding Cross-covariance, I attached a .r file containing a function I wrote.
The Wavelet Cross-covariance (only consider the cross covariance of wc with lag = 0) can be calculated with following steps:

• apply modwt (or dwt) to each univariate TS
• for each pair of 2 univariate TSs, calculate their Cross-covariance (Cov(the i-th wc of the first TS at level j, the i-th wc of the second TS at level j)) at every level j (j = 1, ... J)

Apart from the wccv, the function also returns the variance of each wccv and its 95% CI. The variance is calculated based on the spectrum density (of the cross product of wc from different TS) evaluate at 0 (spectrum0 function from library(coda)). spectrum0() may give some warning messages when level j is large (less number of wc).

###################################################
library(coda)

Compute WCCV (based on the modwt function in GMWM package)

--------------------------------------------

INPUT:

Xt = time series matrix with num.ts number of columns and N number of rows

OUTPUT:

A list contains:

wccv: empirical wavelet cross-covariance basd on Haar MODWT

wccv.cov: variance of empirical wavelet cross-covariance

ci_low: low bound of the 95% CI of empirical wavelet cross-covariance

ci_high: high bound of the 95% CI of empirical wavelet cross-covariance

wccv = function(Xt){
num.ts = ncol(Xt)
N = nrow(Xt)
J = floor(log2(N)) - 1
wccv.mat = matrix(NA, num.ts*(num.ts-1)/2, J)
cov.mat = matrix(NA, num.ts*(num.ts-1)/2, J)
ci.low.mat = matrix(NA, num.ts*(num.ts-1)/2, J)
ci.high.mat = matrix(NA, num.ts*(num.ts-1)/2, J)
c = 0
for(i in 1:(num.ts-1)){
coe1 = modwt(Xt[,i])
for(j in (i+1):num.ts){
c = c+1
coe2 = modwt(Xt[,j])
for(k in 1:J){
wccv.mat[c,k] = mean((unlist(coe1[k])) * (unlist(coe2[k])))

    cov.mat[c,k] = 1/(N - 2^k + 1) * spectrum0((unlist(coe1[[k]])*unlist(coe2[[k]])), max.freq = 0.5, order = 2, max.length = 130)$spec
    ci.low.mat[c,k] =  wccv.mat[c,k] + qnorm(0.025) * sqrt(cov.mat[c,k])
    ci.high.mat[c,k] = wccv.mat[c,k] + qnorm(1-0.025) * sqrt(cov.mat[c,k])
  }
}

}
list(wccv = t(wccv.mat), wccv.cov = t(cov.mat), ci_low = t(ci.low.mat), ci_high = t(ci.high.mat))
}