pdSpecEst2D {pdSpecEst}

Intrinsic wavelet HPD time-varying spectral estimation

Description

pdSpecEst2D calculates a (d,d)-dimensional HPD wavelet-denoised time-varying spectral matrix estimator by applying the following steps to an initial noisy HPD time-varying spectral estimate (obtained with e.g., pdPgram2D):

1. a forward intrinsic AI wavelet transform, with WavTransf2D,

2. (tree-structured) thresholding of the wavelet coefficients, with pdCART,

3. an inverse intrinsic AI wavelet transform, with InvWavTransf2D.

The complete estimation procedure is described in more detail in Chapter 5 of (Chau 2018).

Usage

pdSpecEst2D(P, order = c(3, 3), metric = "Riemannian", alpha = 1,
  return_val = "f", ...)

Arguments

P

a (d,d,n1,n2)-dimensional array of HPD matrices corresponding to a rectangular surface of (d,d)-dimensional HPD matrices of size n_1 \times n_2, with n_1 = 2^{J_1} and n_2 = 2^{J_2} for some J_1, J_2 > 0.

order

a 2-dimensional numeric vector (1,1) \le order \le (9,9) corresponding to the marginal orders of the intrinsic 2D AI refinement scheme, defaults to order = c(3, 3).

metric

the metric that the space of HPD matrices is equipped with. The default choice is "Riemannian", but this can also be one of: "logEuclidean", "Cholesky", "rootEuclidean" or "Euclidean". See also the Details section below.

alpha

an optional tuning parameter in the wavelet thresholding procedure. The penalty (or sparsity) parameter in the tree-structured wavelet thresholding procedure in pdCART is set to alpha times the estimated universal threshold, defaults to alpha = 1.

return_val

an optional argument that specifies whether the denoised spectral estimator is returned or not. See the Details section below.

...

additional arguments for internal use.

Details

The input array P corresponds to an initial noisy HPD time-varying spectral estimate of the (d, d)-dimensional spectral matrix at a time-frequency grid of size m_1 \times m_2, with m_1, m_2 dyadic numbers. This can be e.g., a multitaper HPD time-varying periodogram given as output by the function pdPgram2D.
P is transformed to the wavelet domain by the function WavTransf2D, which applies an intrinsic 2D AI wavelet transform based on a metric specified by the user. The noise is removed by tree-structured thresholding of the wavelet coefficients based on the trace of the whitened coefficients with pdCART by minimization of a complexity penalized residual sum of squares (CPRESS) criterion via the fast tree-pruning algorithm in (Donoho 1997). The penalty (i.e., sparsity) parameter in the optimization procedure is set equal to alpha times the universal threshold, where the noise variance of the traces of the whitened wavelet coefficients are determined from the finest wavelet scale. See Chapter 5 of (Chau 2018) for further details.
The function computes the forward and inverse intrinsic 2D AI wavelet transform in the space of HPD matrices equipped with one of the following metrics: (i) the affine-invariant Riemannian metric (default) as detailed in e.g., (Bhatia 2009)[Chapter 6] or (Pennec et al. 2006); (ii) the log-Euclidean metric, the Euclidean inner product between matrix logarithms; (iii) the Cholesky metric, the Euclidean inner product between Cholesky decompositions; (iv) the Euclidean metric; or (v) the root-Euclidean metric. The default choice of metric (affine-invariant Riemannian) satisfies several useful properties not shared by the other metrics, see (Chau and von Sachs 2019) or (Chau 2018) for more details. Note that this comes at the cost of increased computation time in comparison to one of the other metrics.
If return_val = 'f' the thresholded wavelet coefficients are transformed back to the time-frequency domain by the inverse intrinsic 2D AI wavelet transform via InvWavTransf2D, returning the wavelet-denoised HPD time-varying spectral estimate.

Value

The function returns a list with the following five components:

f

a (d,d,m1,m2)-dimensional array of HPD matrices, corresponding to the HPD wavelet-denoised estimate on the same resolution grid of size m_1 \times m_2 as specified by the input array P. If return_val != 'f', the inverse wavelet transform of the thresholded wavelet coefficients is not computed and f is set equal to NULL.

D

the 2D pyramid of threshold wavelet coefficients. This is a list of arrays, where each array contains the rectangular grid (d,d)-dimensional thresholded wavelet coefficients from the coarsest wavelet scale j = 0 up to the finest wavelet scale j = jmax.

M0

a numeric array containing the midpoint(s) at the coarsest scale j = 0 in the 2D midpoint pyramid.

tree.weights

a list of logical values specifying which coefficients to keep, with each list component corresponding to an individual wavelet scale starting from the coarsest wavelet scale j = 0.

D.raw

the 2D pyramid of non-thresholded wavelet coefficients in the same format as the component $D.

References

Bhatia R (2009). Positive Definite Matrices. Princeton University Press, New Jersey.

Chau J (2018). Advances in Spectral Analysis for Multivariate, Nonstationary and Replicated Time Series. phdthesis, Universite catholique de Louvain.

Chau J, von Sachs R (2019). “Intrinsic wavelet regression for curves of Hermitian positive definite matrices.” Journal of the American Statistical Association. doi: 10.1080/01621459.2019.1700129.

Donoho D (1997). “CART and best-ortho-basis: a connection.” The Annals of Statistics, 25(5), 1870–1911.

Pennec X, Fillard P, Ayache N (2006). “A Riemannian framework for tensor computing.” International Journal of Computer Vision, 66(1), 41–66.

See Also

pdPgram2D, WavTransf2D, InvWavTransf2D, pdCART

Examples

## Not run: 
P <- rExamples2D(c(2^6, 2^6), 2, example = "tvar")$P
f <- pdSpecEst2D(P)

## End(Not run)