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:

[x] Classical wavelet variance: this is the classical estimator proposed by Percival. This is already implemented in the gmwm package. This should be straightforward.
[x] 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!
[x] 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.

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)) }

SMAC-Group / wv

Wavelet variances #5