# Estimates spectrum of multivariate time series 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`). 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 [REMOVE_ME]$$ f_{y}(\lambda) = \boldsymbol \beta' f_{\mathbf{X}}(\lambda) \boldsymbol \beta \geq 0, [REMOVE_ME_2]$$ 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)`. ```r mvspectrum( series, method = c("mvspec", "pgram", "pspectrum", "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: use `"pspectrum"` for `pspectrum`; 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 `pspectrum` or `mvspec` (e.g., `taper`) - `mvspectrum.output`: an object of class `"mvspectrum"` representing the multivariate spectrum of $\mathbf{X}_t$ (not necessarily `normalize`d). - `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. ## Returns `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`). It also has an `"normalized"` attribute, 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`. ## 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`). 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)`. ## 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 `method`s in `mvspectrum` (this is highly recommended to obtain more reasonable/stable estimates). ## Examples ```r 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 ## Not run: ss <- mvspectrum(XX[, 1], method="pspectrum", smoothing = TRUE) mvspectrum(XX, normalize = TRUE) ## End(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 = "pgram")) * 2 sum(mvspectrum(xx, normalize = FALSE, method = "mvspec")) * 2 ## Not run: sum(mvspectrum(xx, normalize = FALSE, method = "pspectrum")) * 2 ## End(Not run) ## Not run: 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 = "mvspec") mvspectrum2wcov(normalize_mvspectrum(S.U)) ## End(Not run) XX <- matrix(rnorm(1000), ncol = 2) SS <- mvspectrum(XX, "mvspec") ss1 <- mvspectrum(XX[, 1], "mvspec") 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, "mvspec") ss.yy.comb <- spectrum_of_linear_combination(SS, beta.tmp) ss.yy <- mvspectrum(yy, "mvspec") plot(ss.yy, log = TRUE) # using plot.mvspectrum() lines(ss.yy.comb, col = "red", lty = 1, lwd = 2) ``` ## References See References in `spectrum`, `pspectrum`, `mvspec`.

mvspectrum function