ForeCA (version 0.2.6)

mvspectrum: Estimates spectrum of multivariate time series

Description

The spectrum of a multivariate time series is a matrix-valued function of the frequency \(\lambda \in [-\pi, \pi]\), which is symmetric/Hermitian around \(\lambda = 0\).

mvspectrum estimates it and returns a 3D array of dimension \(num.freqs \times K \times K\). Since the spectrum is symmetric/Hermitian around \(\lambda = 0\) it is sufficient to store only positive frequencies. In the implementation in this package we thus usually consider only positive frequencies (omitting \(0\)); num.freqs refers to the number of positive frequencies only.

normalize_mvspectrum normalizes the spectrum such that it adds up to \(0.5\) over all positive frequencies (by symmetry it will add up to 1 over the whole range -- thus the name normalize).

For a \(K\)-dimensional time series it adds up to a Hermitian \(K \times K\) matrix with 0.5 in the diagonal and imaginary elements (real parts equal to \(0\)) in the off-diagonal. Since it is Hermitian the mvspectrum will add up to the identity matrix over the whole range of frequencies, since the off-diagonal elements are purely imaginary (real part equals 0) and thus add up to 0.

check_mvspectrum_normalized checks if the spectrum is normalized (see normalize_mvspectrum for the requirements).

mvpgram computes the multivariate periodogram estimate using bare-bone multivariate fft (mvfft). Please use mvspectrum(..., method = 'pgram') instead of mvpgram directly.

This function is merely included to have one method that does not require the astsa nor the sapa R packages. However, it is strongly encouraged to install either one of them to get (much) better estimates. See Details.

get_spectrum_from_mvspectrum extracts the spectrum of one time series from an "mvspectrum" object by taking the i-th diagonal entry for each frequency.

spectrum_of_linear_combination computes the spectrum of the linear combination \(\mathbf{y}_t = \mathbf{X}_t \boldsymbol \beta\) of \(K\) time series \(\mathbf{X}_t\). This can be efficiently computed by the quadratic form $$ f_{y}(\lambda) = \boldsymbol \beta' f_{\mathbf{X}}(\lambda) \boldsymbol \beta \geq 0, $$ for each \(\lambda\). This holds for any \(\boldsymbol \beta\) (even \(\boldsymbol \beta = \boldsymbol 0\) -- not only for \(||\boldsymbol \beta ||_2 = 1\). For \(\boldsymbol \beta = \boldsymbol e_i\) (the i-th basis vector) this is equivalent to get_spectrum_from_mvspectrum(..., which = i).

Usage

mvspectrum(
  series,
  method = c("pgram", "multitaper", "direct", "wosa", "mvspec", "ar"),
  normalize = FALSE,
  smoothing = FALSE,
  ...
)
normalize_mvspectrum(mvspectrum.output)
check_mvspectrum_normalized(f.U, check.attribute.only = TRUE)
mvpgram(series)
get_spectrum_from_mvspectrum(
  mvspectrum.output,
  which = seq_len(dim(mvspectrum.output)[2])
)
spectrum_of_linear_combination(mvspectrum.output, beta)

Arguments

series

a \(T \times K\) array with T observations from the \(K\)-dimensional time series \(\mathbf{X}_t\). Can be a matrix, data.frame, or a multivariate ts object.

method

string; method for spectrum estimation; see method argument in SDF (in the sapa package); use "mvspec" to use mvspec (astsa package); or use "pgram" to use spec.pgram.

normalize

logical; if TRUE the spectrum will be normalized (see Value below for details).

smoothing

logical; if TRUE the spectrum will be smoothed with a nonparametric estimate using gam and an exponential family (with link = log). Only works for univariate spectrum. The smoothing parameter is chosen automatically using generalized cross-validation (see gam for details). Default: FALSE.

…

additional arguments passed to SDF or mvspec (e.g., taper)

mvspectrum.output

an object of class "mvspectrum" representing the multivariate spectrum of \(\mathbf{X}_t\) (not necessarily normalized).

f.U

multivariate spectrum of class 'mvspectrum' with normalize = TRUE.

check.attribute.only

logical; if TRUE it checks the attribute only. This is much faster (it just needs to look up one attribute value), but it might not surface silent bugs. For sake of performance the package uses the attribute version by default. However, for testing/debugging the full computational version can be used.

which

integer(s); the spectrum of which series whould be extracted. By default, it returns all univariate spectra as a matrix (frequencies in rows).

beta

numeric; vector \(\boldsymbol \beta\) that defines the linear combination.

Value

mvspectrum returns a 3D array of dimension \(num.freqs \times K \times K\), where

• num.freqs is the number of frequencies

• K is the number of series (columns in series).

Note that it also has an attribute "normalized" which is FALSE if normalize = FALSE; otherwise TRUE. See normalize_mvspectrum for details.

normalize_mvspectrum returns a normalized spectrum over positive frequencies, which:

univariate:

adds up to \(0.5\),

multivariate:

adds up to Hermitian \(K \times K\) matrix with 0.5 in the diagonal and purely imaginary elements in the off-diagonal.

check_mvspectrum_normalized throws an error if spectrum is not normalized correctly.

get_spectrum_from_mvspectrum returns either a matrix of all univariate spectra, or one single column (if which is specified.)

spectrum_of_linear_combination returns a vector with length equal to the number of rows of mvspectrum.output.

Details

For an orthonormal time series \(\mathbf{U}_t\) the raw periodogram adds up to \(I_K\) over all (negative and positive) frequencies. Since we only consider positive frequencies, the normalized multivariate spectrum should add up to \(0.5 \cdot I_K\) plus a Hermitian imaginary matrix (which will add up to zero when combined with its symmetric counterpart.) As we often use non-parametric smoothing for less variance, the spectrum estimates do not satisfy this identity exactly. normalize_mvspectrum thus adjust the estimates so they satisfy it again exactly.

mvpgram has no options for improving spectrum estimation whatsoever. It thus yields very noisy (in fact, inconsistent) estimates of the multivariate spectrum \(f_{\mathbf{X}}(\lambda)\). If you want to obtain better estimates then please use other methods in mvspectrum (this is highly recommended to obtain more reasonable/stable estimates).

References

See References in spectrum, SDF, mvspec.

Examples

# NOT RUN {
set.seed(1)
XX <- cbind(rnorm(100), arima.sim(n = 100, list(ar = 0.9)))
ss3d <- mvspectrum(XX)
dim(ss3d)

ss3d[2,,] # at omega_1; in general complex-valued, but Hermitian
identical(ss3d[2,,], Conj(t(ss3d[2,,]))) # is Hermitian

ss <- mvspectrum(XX[, 1], smoothing = TRUE)

# }
# NOT RUN {
  mvspectrum(XX, normalize = TRUE)
# }
# NOT RUN {
ss <- mvspectrum(whiten(XX)$U, normalize = TRUE)

xx <- scale(rnorm(100), center = TRUE, scale = FALSE)
var(xx)
sum(mvspectrum(xx, normalize = FALSE, method = "direct")) * 2
sum(mvspectrum(xx, normalize = FALSE, method = "wosa")) * 2


xx <- scale(rnorm(100), center = TRUE, scale = FALSE)
ss <- mvspectrum(xx)
ss.n <- normalize_mvspectrum(ss)
sum(ss.n)
# multivariate
UU <- whiten(matrix(rnorm(40), ncol = 2))$U
S.U <- mvspectrum(UU, method = "wosa")
mvspectrum2wcov(normalize_mvspectrum(S.U))


XX <- matrix(rnorm(1000), ncol = 2)
SS <- mvspectrum(XX, "direct")
ss1 <- mvspectrum(XX[, 1], "direct")

SS.1 <- get_spectrum_from_mvspectrum(SS, 1)
plot.default(ss1, SS.1)
abline(0, 1, lty = 2, col = "blue")


XX <- matrix(arima.sim(n = 1000, list(ar = 0.9)), ncol = 4)
beta.tmp <- rbind(1, -1, 2, 0)
yy <- XX %*% beta.tmp

SS <- mvspectrum(XX, "wosa")
ss.yy.comb <- spectrum_of_linear_combination(SS, beta.tmp)
ss.yy <- mvspectrum(yy, "wosa")

plot(ss.yy, log = TRUE) # using plot.mvspectrum()
lines(ss.yy.comb, col = "red", lty = 1, lwd = 2) 

# }