Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
82,128 characters · 15 sections · 0 citation commands
Understanding fluctuations through Multivariate Circulant Singular Spectrum Analysis
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Block circulant matrices, double diagonalization, co-movement, eigenstructure, time series.
\spacingset{1.45}
Signal extraction is an old problem across various disciplines, but typical decomposition procedures differ depending on the field. For instance, engineers work with amplitude and frequency modulated (AM-FM) signals. They usually deal with high frequencies and use different methods to extract possibly highly non-linear fluctuations, such as those based on the Hilbert Transform (Huang et al., 1998; Gianfelici et al., 2007; and Biagetti et al., 2015). However, these methods were not designed to work with low frequencies. On the other hand, economists need to extract long and medium-term information, like the trend and cycle, free of seasonality, for policy analysis. In general, they use parametric models either based on ARIMA formulations (see, e.g. Hillmer and Tiao, 1982) or in the state space framework (see e.g. Harvey, 1989) that may not work well with short or highly nonlinear time series.
Singular Spectrum Analysis (SSA) is another alternative for signal extraction based on subspace algorithms (see, for instance, the surveys by Ghil et al., 2002, or Golyandina and Zhigljavsky, 2013). After choosing a window length $L$, SSA builds a related trajectory matrix by putting together lagged pieces of the original time series and performs its Singular Value Decomposition. Different alternatives of SSA, such as Basic and Toeplitz SSA, need to identify the frequencies associated with the estimated components after they have been extracted. Within this framework Circulant Singular Spectrum Analysis (CiSSA), is a nonparametric procedure developed in Bógalo et al. (2021) for univariate time series, that allows for an automated matching between the extracted signals and their frequency of oscillation.
The multivariate extension of Singular Spectrum Analysis, M-SSA, (Broomhead and King, 1986b) appears simultaneously with univariate SSA (Broomhead and King, 1986a; Fraedrich, 1986). M-SSA, as is the case with its univariate counterpart, also needs to identify the frequencies of the reconstructed components and, as Plaut and Vautard (1994) observe, is able to extract patterns both along time and across time series. While originally only applied to climatology, since the work of Ghil et al. (2002) who compared different versions, the application of this technique has been extended to a wide range of disciplines: biometrics (see, e.g., Lee et al., 2014), seismic activity (see, e.g., Cheng et al., 2019), geolocalization (see, e.g., Gruszczynska et al., 2017), business cycles and economics (see, e.g., Carvalho (de) and Rua, 2017; Hassani et al., 2013; Silva et al., 2018) and medical diagnostic (see, e.g., Jain et al., 2020).
Multivariate Circulant Singular Spectrum Analysis (M-CiSSA) is a novel unified framework for multivariate signal extraction. It builds a new trajectory matrix, different to previous M-SSA, and a related block circulant matrix of second order moments that allows us to compute the cross-spectral density matrices at different frequencies in a straightforward way. It is based on the properties of the eigen-structure of a block circulant matrix related to the second order moments of the vector of time series and their lags. Also, further diagonalization of these matrices enables us to obtain the eigenvectors and eigenvalues that provide the principal components and their contribution to the total variability within a specific frequency. Both aspects, block frequency identification and decomposition within a block, are new contributions in this setup. As a consequence, we can prove a uniqueness theorem that states the reconstruction of the univariate components by the sum of the multivariate subcomponents per frequency.
M-CiSSA offers several advantages. Firstly, eigenvalues and eigenvectors of each of these blocks resulting from the block diagonalization contain all the variability of the corresponding frequencies. This enables us to understand co-movements at different frequencies and even the cyclical position among variables. Secondly, the uniqueness theorem is new to SSA and is very useful to further understand the formation of the individual cycles in terms of the multivariate common drivers. Thirdly, the solution proposed can be applied to any range of frequencies and to very different types of time series.
All in all, our approach provides a comprehensive framework for jointly analyzing fluctuations, solving, in a unified way, the problem of extracting the underlying components of a set of time series, disentangling their sources of variation and assessing their relative cyclical position at each frequency. In the particular case of common information and without the need for fitting a factor model, M-CiSSA also identifies common sources of variation and assigns them to a particular frequency. It can also serve as a tool for denoising the latent signals.
We illustrate the performance of this technique through two examples. First, we use a synthetic example of two series, each of them generated as the sum of various latent signals including amplitude and modulated ones. We show that M-CiSSA is able to clearly separate the underlying signals. We complete our illustrations with the analysis of a real data set of Primary Commodity Energy Prices in order to characterize the main latent signals (trends, cycles) driving their fluctuations. We found co-movements in the long and medium term for all series except US natural gas, and identified the additional decoupling of coals and Japanese natural gas only at the medium-term cyclical frequency. Finally, we have discovered which prices can be held or sustained by forces outside the common evolution of the markets.
The remainder of the paper is structured as follows. Section 2 presents the proposed methodology. Section 3 shows the two data applications. Finally, we draw our conclusions in section 4.
In this section we will present our new proposal that generalizes the univariate CiSSA (Bógalo et al. 2021). The theory behind CiSSA starts, as any other SSA procedure, with the transformation of the original series into a related trajectory matrix, its decomposition into elementary units, and further reconstruction to regain the original time series dimensions. The novelty is that the decomposition is based on a circulant second order matrix that guarantees an explicit expression for the eigenvalues and eigenvectors that relates them to frequencies. This perfect match between frequencies and the eigen-structure of the matrix allows us for an automated identification of the signals in the time series. The multivariate M-CiSSA follows a similar approach that substitutes the circulant matrix by a block circulant matrix that allows us to block diagonalize and match each block with a frequency. While each block contains all the information related to the multivariate frequency, further diagonalization within the block will identify the main drivers that explain the fluctuations.
For a better understanding of the new proposal, we will first review the univariate algorithm. Afterwards, we will introduce the steps of the multivariate version highlighting the novelties and generalizations that were needed compared to the univariate case. Finally, we will present the auxiliary results that are the pillars of our new proposal.
In what follows, let $\left\{ x_{t}\right\}_{t=1} ^{T}$ be a realization of a stochastic process $\left\{ x_{t}\right\}$ $t \in \cal T$ and $\mathbf{x}=(x_{1},...,x_{T})^{\prime}$, where the prime denotes transpose and $L$ a positive integer, called the window length, such that $1<L<T/2$.\footnote{See that we use the same notation for the stochastic process and for the observed time series. When needed, it will be explicitly clarified in the main text which one we are referring to.} The algorithm works in four steps\footnote{For additional details, see Bógalo et al. (2021).}:
1st step: Embedding. Select the window length $L$ and convert the univariate time series into a matrix by defining a trajectory matrix $\mathbf{X}$, $N=T-L+1$, as follows
where $\mathbf{x}_{j}=(x_{j},...,x_{j+L-1})^{\prime }$ indicates the $L \times 1$ vector with origin at time $j$.
2nd step: Decomposition Compute the second moments
and define
that are the elements of the circulant matrix
and are also needed to compute its eigenvalues
The associated eigenvectors are given by $k=1,...,L$
where $u_{k,j}=\exp \left( -i2\pi (j-1)\frac{k-1}{L}\right)$.
Form the elementary matrices of rank 1, $\mathbf{X}_{k}=\mathbf{u}_{k}\mathbf{u}_{k}^*\mathbf{X}$, where the superindex $*$ denotes the conjugate transpose.
3rd step: Grouping. Based on the following relationship between the eigenvalues and the spectral density function of the data $f$ (see Lancaster, 1969)
associate the $k$-th eigenvalue and corresponding eigenvector to the frequency $w_k=\frac{k-1}{L},k=1,...,L.$ Group the elementary matrices $\mathbf{X}_{k}$ into $G$ disjoint groups. For each group $j=1,...,G$, we can define the matrix $\mathbf{X}_{I_{j}}=\sum_{k\in I_j} \mathbf{X} _k$ summing all the elementary matrices within this group.
4th step: Reconstruction. Convert each matrix from step 3 into a time series of the original dimension, denoted as $\widetilde{\mathbf{x}}^{(j)}=( \widetilde{x}_{1}^{(j)},...,\widetilde{x}_{T}^{(j)})^{\prime }$ by diagonal averaging. Denoting by $\widetilde{x}_{r,s}$ the elements of the matrix $\mathbf{X}_{I_{j}}$, the reconstruction is done by averaging the elements of this matrix over its antidiagonals,
and constitutes the extracted signals associated to a particular frequency or range of frequencies.
Note that with CiSSA the analyst just needs to choose the desired frequencies beforehand to obtain the related signals. This is different from previous attempts to automate SSA where the components are first estimated and in a second step associated to a frequency, see Ghil and Mo (1991), Vautard et al. (1992), Alexandrov and Golyandina (2005), Alexandrov (2009), Alonso and Salgado (2008), Bilancia and Campobasso (2010), Arteche and García-Enríquez (2017), Carvalho (de) and Rúa (2017) and recently, Golyandina and Zhornikova (2023), among others.
The goal we pursue with M-CiSSA is more ambitious as it aims to understand the formation of fluctuations by frequency, uncovering their commonalities and specificities, and also the relative cyclical position of the signals. The M-CiSSA solution is based on the block diagonalizaton of a block circulant matrix related to the second moments of the data. Because of the properties of block circulant matrices, there is an exact match that associates a block with a specific known frequency. Further diagonalization within each block explains the cross-section variability by frequency and, therefore, allows us to understand how series are related at each frequency. A new insight of the multivariate problem is the so-called multichannel analysis or the knowledge of how each series contributes to the formation of the signal. We will extract the unobserved signals from several channels or time series in contrast to the univariate analysis where the trend and different cycles were only fed by one series or channel. When the different time series are measures of a particular variable in different points of a grid, this gives a space-time decomposition. The eigenvalue problem to be solved is similar to that of the univariate case but replacing each data point by a vector of $M$ time series. This rises another issue and is how each channel or univariate time series contributes to the formation of the different unobserved signals. In this sense, each element of the eigenvectors in univariate analysis is now replaced by a segment or subvector of size $M$. Putting together the elements coming from each time series, we can find the contribution of each one of them to the formation of the unobserved signals. We will make explicit in our proposal how to identify the "pieces" of eigenvectors corresponding to each channel or time series.
So, in what follows, we will first introduce the new trajectory matrix; then, we will present the four steps of the M-CiSSA procedure and, finally, we will justify the new algorithm.
We start by defining the new trajectory matrix. Let $ \mathbf{x}_t=\left(x_t^{\left(1\right)},\cdots,x_t^{\left(M\right)}\right)'$, be an $M$-dimensional stationary and, for the sake of simplicity, zero mean stochastic process $t \in \cal T$ with autocovariance matrices $ \mathbf{\Gamma}_k= E[\mathbf{x}_{t+k}\mathbf{x}_t^{\prime}], k=0,\pm1,\cdots,\pm\left(L-1\right)$. Consider now its realization of length $T$. Given a window length $L<T/2$ we construct a trajectory matrix, of dimensions $LM\times N$ with $ N=T-L+1$ as,
Notice that this construction of the trajectory matrix differs from that of other SSA multivariate algorithms\footnote{The trajectory matrix in other SSA multivariate algorithms stacks the univariate trajectory matrices instead of considering each value of the matrix as a vector of time series.}. The consideration of this alternative trajectory matrix is a relevant issue, since its matrix of second order moments is block Toeplitz and will lay the foundation for the construction of the related block circulant matrix.
The four steps of the new M-CiSSA algorithm present several novelties. First, on the embedding step we will use the new trajectory matrix that we have previously defined. This will allow us to define the blocks of the elementary matrices associated with each time series in the decomposition step. The diagonalization will be done in two steps: block diagonalization, first, followed by a further diagonalization within each block. Putting together the pieces of eigenvectors corresponding to the same time series is easy in our case as the element of the subvector corresponding to series $i$ occupies always the same position $i$ within the subvector and allows us to easily identify the contribution of each channel or time series to the formation of the unobserved components. As with univariate CiSSA, the grouping will be automated and we extend the results of the univariate case, where we match frequencies with eigenstructure, to the multivariate setup. As an auxiliary result, we need to show that the new proposed circulant matrices are asymptotically equivalent to the traditional variance-covariance matrices. The reconstruction step is done as usual in SSA algorithms. In what follows we are going to, first, present the 4 steps of the algorithm and, afterwards, to show the results that allow us to introduce the aforementioned novelties
The four steps of the M-CiSSA algorithm are:
1st step: Embedding.
Form the big trajectory matrix $\mathbf{X}$ as in ((ref)).
2nd step: Double decomposition.
Compute the sample second moment matrices $\hat{\mathbf{\Gamma}}_k=\frac{1}{T-k}\sum_{t=1}^{T-k}\mathbf{x}_{t+k}\mathbf{x}_t^\prime$, $k=1,\cdots,L-1$ and from them, the linear combinations
Build the block circulant matrix $\mathbf{S}_\mathbf{C}$ given by
Diagonalize the $LM \times LM$ matrix $\mathbf{S}_\mathbf{C}$ in two steps.
Notice that the Fourier unitary matrix defined in ((ref)) has complex values. Therefore, we would like to make a change of basis so that the new set of eigenvectors that diagonalize $\mathbf{S}_C$, are all real. Proposition 2 (see Section 2.4.2), allows to substitute $\mathbf{V}$ by a real orthonormal basis $\mathbf{\tilde{V}}$ that diagonalizes $ \mathbf{S}_C$ with elements $\widetilde{\mathbf{V}}=\left[{\widetilde{\mathbf{v}}}_1|\cdots|{\widetilde{\mathbf{v}}}_{LM}\right]$, being ${\widetilde{\mathbf{v}}}_j={\widetilde{\mathbf{v}}}_{k,m}$ with $j=\left(k-1\right)M+m\quad, \forall \: k=1,\cdots,L$ and $m=1,\cdots,M$ defined by
where $\mathcal{R}_{\mathbf{v}}$ denotes the real part of the vector $\mathbf{v}$.
Therefore, we can compute the elementary matrix for the {\it m-th} subcomponent at frequency $\omega_k$ as
Consequently, the trajectory matrix can be decomposed as
and the contribution to the total variability of the elementary matrix $\mathbf{X}_{k,m}$ is the ratio
The piece of ${{\mathbf{\tilde{v}}}_{k,m}}$ corresponding to the series $(i)$ is denoted by $\mathbf{\tilde{v}}_{k,m}^{\left( i \right)}$ and is given by
where ${{\mathbf{1}}_{M,i}}$ is a vertical vector of length $M$ with a $1$ in the {\it i-th} position and $0$ elsewhere. Therefore, the elementary matrix of the series $(i)$ for the {\it m-th} subcomponent of the frequency $\omega_k$ is obtained as
In addition, the participation index of the series $(i)$ within the {\it m-th} subcomponent of the frequency $\omega_k$ is calculated, from ((ref)), by the following expression
3rd step: Grouping.
The spectral density function of the vector process $\mathbf{x}_t$ is symmetric and, therefore, $\mathbf{\hat{F}}_{k}^{{}}=\mathbf{\hat{F}}_{L+2-k}^{\prime}$ for $k=2,\cdots ,\left\lfloor \tfrac{L+1}{2} \right\rfloor$. Then $\hat{\mathbf{D}}_k=\hat{\mathbf{D}}_{L+2-k}$ and, consequently, the subcomponents $\mathbf{w}_{k,m}=\mathbf{X}'\widetilde{\mathbf{v}}_{k,m}$ and $\mathbf{w}_{L+2-k,m}=\mathbf{X}'\widetilde{\mathbf{v}}_{L+2-k,m}$ are harmonics of the same frequency. This motivates the creation of elementary pairs per subcomponent and frequency $B_{k,m}=\left\{\left(k,m\right),\left(L+2-k,m\right)\right\}$ for $k=2,\cdots,\left\lfloor \tfrac{L+1}{2} \right\rfloor$ except $B_{1,m}=\left\{\left(1,m\right)\right\}$ and occasionally $B_{\frac{L}{2}+1,m}=\left\{\left(\frac{L}{2}+1,m\right)\right\}$ if $L$ is even. The matrices corresponding to the pairs $B_{k,m}$ are given by the sum of two elementary matrices per subcomponent and frequency
Therefore, both the matrices associated with the elementary pairs $B_{k,m}$ by subcomponent and frequency, and the oscillatory components obtained from them are previously identified with a determined frequency as in the univariate case.
Under the assumption of separability, we define $G$ disjoint groups of the elementary pairs per subcomponent and frequency. The resulting matrix for each of the disjoint groups is defined as the sum of the matrices associated with the pairs $B_{k,m}$ included. If $I_j=\left\{ B_{k_{j_1},m_{j_1}}, ..., B_{k_{j_q},m_{j_q}}\right\}$, $j=1,...,G$ is each disjoint group of $j_q$ pairs $B_{k,m}$ with $ 1\le j_q\le \left\lfloor \tfrac{L+1}{2} \right\rfloor M$, then the matrix $\mathbf{X}_{I_j}$ from group $I_j$ is calculated as the sum of the corresponding matrices defined by ((ref)), $ \mathbf{X}_{I_j}=\sum_{B_{k,m}\in I_j}\mathbf{X}_{B_{k,m}}$ and the trajectory matrix $\mathbf{X}$ can be decomposed as the sum of these group matrices
If we wish to extract the oscillation at frequency $\frac{k-1}{L}$, the appropriate group is $B_k=\left\{B_{k,m}\; \forall\: m, 1\le m\le M\right\} $ being the resulting matrix $ \mathbf{X}_{B_k}=\sum_{m=1}^{M}\mathbf{X}_{B_{k,m}}$
4th step: Reconstruction.
Finally, each $L\times\ N$ matrix ${{\mathbf{X}}_{{{I}_{j}}}}$ of vectors $M\times 1$ from the previous step is transformed into a new time series of length $T$ by diagonal averaging, producing the reconstructed time series $\widetilde{\mathbf{x}}_{I_j,t}=\left({\widetilde{x}}_{I_j,t}^{\left(1\right)},\cdots,\widetilde{x}_{I_j,t}^{\left(M\right)}\right)'$. If $\mathbf{x}_{r,s}^{I_j}$ are the vector elements of the matrix $\mathbf{X}_{I_j}$, then the values of the reconstructed vector time series $ \widetilde{\mathbf{x}}_{I_j,t}$, with $L<N$ are calculated as in Vautard et al. (1992) but the formula is adapted to vectors to “hankelize” the matrix $\mathbf{X}_{I_j}$ with the operator $\text{H}(\centerdot)$ as follows:
The reconstructed multivariate times series resulting from $B_{k,m}$ are called elementary reconstructed vector time series by subcomponent $m$ and frequency $k$.
The pseudo-code with the logical sequence of the details of the four steps of M-CiSSA is provided in Algorithm 1.
\spacingset{1}
\spacingset{1.45}
The M-CiSSA algorithm described so far requires stationary time series. However, it is straightforward to show that it can also be applied to non-stationary time series. Other versions of multivariate SSA, Basic SSA (Broomhead and King, 1986b) and Toeplitz SSA (Plaut and Vautard, 1994), are also implemented on non-stationary time series. In the case of Circulant SSA, Bógalo et al. (2021) prove its validity for non-stationary univariate time series approximating the discontinuities of the spectrum by means of a pseudo-spectrum. The same approximation can be applied in the case of $M$ series instead of just one. More recently, there are examples of applying alternative versions of SSA to non-stationary multivariate time series (see, e.g., Groth et al., 2011; Carvalho (de) and Rua, 2017). In a related but different context, Peña and Yohai (2016) also use principal components and lagged principal components with non-stationary time series.
Consider the sequence of lagged cross variance-covariance matrices of the population as a function of the window length $L$ and denote it by $\mathbf{T}_L$. Each matrix in that sequence is an $L\times L$ block Toeplitz matrix with $M\times M$ blocks resulting in a $LM\times LM$ Hermitian matrix. Let,
It is well known that the sequence $\left\{ \mathbf{\Gamma}_k \right\} _{k\in\mathbb{Z}}$ can be generated as
where $\omega\in\left[0,1\right]$ is the frequency in cycles per unit of time, $\text{i}=\sqrt{-1}$ is the imaginary unit and $\mathbf{F}\left(\omega\right)$ is the spectral density matrix of the stochastic process of the vector time series $\mathbf{x}_t$, that is, the ”matrix” Fourier series given by
The sequence $\left\{\mathbf{\Gamma}_k\right\}_{k\in\mathbb{Z}}$ are the Fourier coefficients of the matrix-valued function $\mathbf{F}\left(\omega\right)$. As a consequence, the continuous and $2\pi$-periodic matrix-valued function $\mathbf{F}\left(\omega\right)$ of a real variable is the generating function or symbol of the matrix $\mathbf{T}_L\left(\mathbf{F}\right)=\mathbf{T}_L$ that originates the sequence of block Toeplitz matrices that we denote $\left\{\mathbf{T}_L\left(\mathbf{F}\right)\right\}$.
With this approach, we do not have a closed formula for the eigenvalues and eigenvectors of the Toeplitz matrix of second moments given in ((ref)). This problem can be solved, if instead of block Toeplitz matrices, we use block circulant matrices.
Let $\mathbf{C}_L$ be an $L\times L$ block circulant matrix with blocks $M\times M$, that is, $\mathbf{C}_L$ is an $LM\times LM$ matrix of the form
where $\mathbf{\Omega}_k\in\mathbb{C}^{M\times M},\; k=0,1,\cdots,L-1$. We say that $\mathbf{C}_L$ is block circulant since each block row is built from a right shift of the blocks of the previous block row. Each block of $\mathbf{C}_L$ can be generated according to Gutiérrez-Gutiérrez and Crespo (2008) as
being $\mathbf{F}\left(\omega\right)$ the generating function of the sequence $\left\{\mathbf{T}_L\left(\mathbf{F}\right)\right\}$ that also generates the sequence of block circulant matrices $\left\{\mathbf{C}_L\left(\mathbf{F}\right)\right\}$. Moreover, as Gutiérrez-Gutiérrez and Crespo (2008), Lemma 6.1, show the two block matrices sequences $\left\{\mathbf{T}_L\left(\mathbf{F}\right)\right\}$ and $\left\{\mathbf{C}_L\left(\mathbf{F}\right)\right\}$ are asymptotically equivalent as $L\rightarrow \infty $, $\mathbf{T}_{L}(\mathbf{F})$ $\sim $ $\mathbf{C}_{L}(\mathbf{F})$, in the sense that both matrices have bounded eigenvalues and $\underset{L\rightarrow \infty }{\lim }\frac{ \left\Vert \mathbf{T}_{L}(\mathbf{F})-\mathbf{C}_{L}(\mathbf{F})\right\Vert_{F}}{\sqrt{L }}=0$, where $\left\Vert \text{\textperiodcentered }\right\Vert_{F}$ is the Frobenius norm.
The advantage of using the block circulant matrix $\mathbf{C}_L$ instead of the block Toeplitz matrix $\mathbf{T}_{L}$ is that the former can be block diagonalized, while the later cannot. However, in order to build the block matrix ((ref)), we should either know the matrix function $\mathbf{F}$ or the infinite sequence $\left\{\mathbf{\Gamma}_k\right\}_{k\in\mathbb{Z}}$. We realize that, in practice, we will have a finite number of second order matrices so, in order to make this approach operational, we generalize the results given by Pearl (1973) for the continuous and $2\pi$-periodic scalar function $f$ to the continuous and $2\pi$-periodic matrix function $\mathbf{F}$. In particular, similarly to the proposal put forward by Pearl (1973) for $1 \times 1 $ blocks, we suggest using as $M \times M$ blocks of the first row in $\mathbf{C}_L$ given by ((ref)):
This way of defining the block circulant matrix $\mathbf{C}_L$ is associated with the continuous and $2\pi$-periodic matrix function ${\widetilde{\mathbf{F}}}$,
that generates the sequence of block circulant matrices $\left\{\mathbf{C}_L\left(\widetilde{\mathbf{F}}\right)\right\}$. Theorem (ref) shows the asymptotic equivalence between the sequences $\left\{\mathbf{T}_L\left(\mathbf{F}\right)\right\}$ and $\left\{\mathbf{C}_L\left(\widetilde{\mathbf{F}}\right)\right\}$, denoted as $\mathbf{T}_{L}(\mathbf{F})$ $\sim $ $\mathbf{C}_{L}(\widetilde{\mathbf{F}})$. This theorem is the basis of our proposed algorithm since we will use the later sequences in our Multivariate Circulant SSA proposal, which we have called M-CiSSA.
Once the block circulant matrix is defined, we check the properties of its diagonalization. In particular, we see that the diagonal blocks are associated with known frequencies and that further diagonalizing within the blocks is possible. As a consequence, these results allow us to understand fluctuations by frequency and within frequency.
To do so, we start by showing how the block circulant matrix sequence $\left\{\mathbf{C}_L\left(\mathbf{F}\right)\right\}$ is diagonalized. For any matrix valued function $\mathbf{F}:\left[0,1\right]\rightarrow\mathbb{C}^{M\times M}$ which is continuous and $2\pi$-periodic, the block circulant matrix $\mathbf{C}_L\left(\mathbf{F}\right)$ is characterized according to Gutiérrez-Gutiérrez and Crespo (2008) by a block diagonalization given by
where $\mathbf{U}_L$ is given by ((ref)). Each block $\mathbf{F}_k=\mathbf{F}\left(\frac{k-1}{L}\right)$, $k=1,\cdots,L$ represents the cross spectral density matrix of the multivariate stochastic process $\mathbf{x}_t$ for the frequency $\omega_k=\frac{k-1}{L}, k=1,\cdots,L$ and can be unitarily diagonalized. In this way, we obtain $\mathbf{F}_k=\mathbf{E}_k\mathbf{D}_k\mathbf{E}_k^*$ with $\mathbf{E}_k\in\mathbb{C}^{M\times M}$ and $\mathbf{E}_k\mathbf{E}_k^*=\mathbf{E}_k^*\mathbf{E}_k=\mathbf{I}_M$, where $\mathbf{E}_k=\left[\mathbf{e}_{k,1}|\cdots|\mathbf{e}_{k,M}\right]$ contains the eigenvectors, and the diagonal matrix $\mathbf{D}_k=\operatorname{diag}{\left(\lambda_{k,1},\cdots,\lambda_{k,M}\right)}$ contains the ordered eigenvalues $\lambda_{k,1}\geq\cdots\geq\lambda_{k,M}\geq0$ of $\mathbf{F}_k$. Therefore, the unitary diagonalization of the Hermitian matrix $\mathbf{C}_L\left(\mathbf{F}\right)$ is given by $\mathbf{C}_L\left(\mathbf{F}\right)=\mathbf{VD}\mathbf{V}^*$ with
where $\mathbf{E}=\operatorname{diag}{\left(\mathbf{E}_1,\cdots,\mathbf{E}_L\right)}$ and $\mathbf{D}=\operatorname{diag}{\left(\mathbf{D}_1,\cdots,\mathbf{D}_L\right)}$. As a consequence, there are $M$ eigenvectors associated with each frequency $\omega_k=\frac{k-1}{L}$. The {\it j-th} eigenvector $\mathbf{v}_j,\; j=1,...,LM$ of the matrix $\mathbf{C}_L\left(\mathbf{F}\right)$ is given by
for $k=1,\cdots,L$ and $m=1,\cdots, M$ where $\mathbf{u}_k$ is the {\it k-th} column of the Fourier unitary matrix $\mathbf{U}_L$ and $\mathbf{e}_{k,m}$ is the {\it m-th} eigenvector of the cross spectral density matrix $\mathbf{F}_k$. Notice that $ \mathbf{F}$ is symmetric with respect to the frequency $ \frac{1}{2}$ as deduced from ((ref)). This means $\mathbf{F}_k=\mathbf{F}_{L+2-k}^{\prime};\; k=2,\cdots,\left\lfloor\frac{L+1}{2}\right\rfloor$. So the corresponding eigenvectors are conjugated $ \mathbf{E}_k=\overline{\mathbf{E}}_{L+2-k}$ and the associated eigenvalues are equal $ \mathbf{D}_k=\mathbf{D}_{L+2-k}$. Proposition (ref) states how to orthogonally diagonalize the circulant matrix $\mathbf{C}_L\left(\mathbf{F}\right)$.
Given the {\it i-th} variable of the set of time series, its corresponding signal for a particular frequency can be estimated in two different ways: in a univariate framework (or M-CiSSA with $M=1$) or within the multivariate setup. In this section, we prove that the univariate signal is the result of adding up the $M$ subcomponents within a given frequency for that series. This is relevant, for instance, to understand the formation of the cycles associated to each time series. While we could extract the cycle within the univariate framework, in doing so within the multivariate setup allows us to decompose the univariate cycle as sum of subcomponents that reflect the relation among the different variables. In particular, we are able to disentangle which part of the univariate cycle is common to the other variables and which part is specific or idyosincratic and not shared with the rest of the variables.
The univariate estimation in CiSSA of the oscillatory component of the {\it i-th} series for each frequency ${{\omega }_{k}}=\tfrac{k-1}{L}$, $k=1,\cdots ,L$, originates a single time series or component associated with the elementary matrix by frequency $\mathbf{X}_{k}^{\left( i \right)}$. However, the estimation of that same oscillatory component by M-CiSSA produces $M$ series or subcomponents respectively associated with the $ M$ elementary matrices by subcomponent and frequency $\mathbf{X}_{k,m}^{\left( i \right)}$. The sum of these matrices, $\sum\limits_{m=1}^{M}{\mathbf{X}_{k,m}^{\left( i \right)}}$, originates the estimation in M-CiSSA of the oscillatory component of the {\it i-th} series at frequency ${{\omega }_{k}}$. Theorem (ref) proves that the two signals are the same.
This result allows us to use M-CiSSA as a way to de-noise the extracted signals and also to obtain the common spectral signals and estimate co-movements. Regarding de-noising, M-CiSSA allows us to separate the signal from an over imposed colored noise, as defined by Allen and Robertson (1996), by selecting a reduced number of components associated with the non-null eigenvalues estimated for the frequency $\omega_k$. In general, to estimate the signal of a harmonic, it will be enough to select a small number of subcomponents associated with their higher eigenvalues that will characterize both the amplitude and the dating of the oscillatory components. In this sense, the subcomponents help to extract the common spectral signals as well as the co-movements in order to analyze their characterization (procyclical or anticyclical) and cyclical position (leading, coincident, or lagging) of the oscillatory components as Groth et al (2011) observed. Therefore, these subcomponents describe the formation of the oscillatory components in a multivariate setup.
Theorem (ref) also proves the empirical result by Plaut and Vautard (1994) that, in M-SSA, an oscillatory pair does not explain all the variability due to a harmonic. Classical versions of M-SSA only consider the highest eigenvalues and omit information related to that harmonic. This information may hold great interest for economic analysis because it allows us to observe, for each of the series, the different gaps for the same oscillatory component.
The following synthetic example shows the performance of M-CiSSA with various types of signals. We consider two ($M=2$) signals $x_{1}(t)$ and $x_{2}(t)$ each one of them generated as the sum of a linear trend $T_{i}(t)$ plus a signal modulated in amplitude $S_{i}(t)$ plus another oscillatory component modulated both in amplitude and frequency $Y_{i}(t)$, such that $x_{i}(t)=T_{i}(t)+S_{i}(t)+Y_{i}(t)$, $i=1,2$. The trend has positive slope for the first series $T_{1}(t)=0.5t$ and negative for the second one $T_{2}(t)=-0.25t$. The AM components are generated as $S_{1}(t)=A_{S}(t)sin(\omega_{S}t)$ and $S_{2}(t)=A_{S}(t)sin(\omega_{S}t-\pi/2)$ being out of phase $\pi/2$. Finally, the AM-FM component are generated as $Y_{i}(t)=A_{Y}(t)sin(\omega_{Y_{i,a}} t+ \omega_{Y_{i,b}}\frac{t^2}{2T})$, $i=1,2$. The modulated amplitudes are generated as $A_{S}(t)=2+0.3cos(\omega_{A,S} t)$ and $A_{Y}(t)=1+0.1cos(\omega_{A,Y} t)$. Notice that the frequency modulated signals show linearly increasing frequencies in time, $\omega_{Y_{i}}=\omega_{Y_{i,a}} + \omega_{Y_{i,b}}\frac{t}{2T}$, $i=1,2$. The chosen values for the frequencies are $f_{S}=125Hz$, $f_{A,S}=1Hz$, $f_{A,Y}=5Hz$, $f_{Y_{1,a}}=50Hz$, $f_{Y_{2,a}}=180Hz$ and $f_{Y_{1,b}}=f_{Y_{2,b}}=40Hz$, being $\omega=2\pi f$. Similar synthetic signals have been used, for example, in Biagietti et al. (2015) and Bógalo et al. (2021). The sampling frequency is 1000$Hz$ and the signals are observed for 10 seconds, therefore the number of observations is 10000. The left panels in Figure (ref), show the simulated time series (first row) and the AM and AM-FM signals (second and third rows respectively) for a span of 2 seconds.
We choose $L=200$ and perform the block eigendecomposition by frequency that contains the spectral information associated to the frequency $w_k=\frac{k-1}{L}, k=1, ..., 200$. Each block is characterized by the 2 eigenvalues $\hat{\lambda}_{k,m}, m=1,2$ of the diagonal matrix $\hat{\textbf{D}}_k, k=1, ..., 200$ that also define the contribution to the total variability of the bivariate system as in ((ref)). The top right panel of Figure 1 shows the estimation of the spectral density of the bivariate system by the first and second (dynamic) eigenvalues, measured in dB, for each frequency, where in the x-axis we have represented both, the values of $k$ and the equivalent normalized frequencies $w_k$ to highlight the automated identification that M-CiSSA provides.
First, it can be seen that in this particular example the first eigenvalue dominates the second one for every frequency. Therefore, it would be enough to just analyze the first eigenvalue, instead of the trace for each block, to identify the most relevant frequencies of fluctuation. Notice that all the values of the second eigenvalue are negative (the scale is logarithmic) and, therefore, its information content is negligible. Second, it clearly shows a peak at the zero frequency $(k=1)$ that captures the linear trend and another one at the normalized frequency of 0.125 ($k=26$) that corresponds to $S_{i}(t), i=1,2$ and the 2 "plateau" corresponding to the frequency modulated signals $Y_1(t)$ and $Y_2(t)$. The first "plateau" goes from $k=11$ to $k=19$ and the second one from $k=33$ to $k=41$, that translated into frequencies correspond to the normalized frequencies between 0.05 to 0.07 for the first one and frequencies between 0.18 and 0.20 for the second one. Therefore, M-CiSSA will capture the modulated frequency by adding components of adjacent frequencies. These components are the more relevant in the bivariate system and they account for 66.25% of the total variability.
Once the analysis of the block diagonalization is complete, we study the within blocks diagonalization. Our goal in this second step is to isolate and reconstruct the latent signals as well as to know how much each of the observed series or channels contributes to that reconstruction. To illustrate that M-CiSSA is accurate to reproduce AM signals, the middle row of Figure (ref) shows the synthetic signals $S_1(t)$ and $S_2(t)$ on the left column and their reconstruction on the right one, for a period of 2 seconds. The two signals are reconstructed using, first, ((ref)) for $k=26$ and its symmetric value in $L+2-k=176$ which correspond to the frequency of 125Hz; and then, the antidiagonal averaging given in ((ref)). The plot in the middle column, right panel shows the two columns of the reconstructed series ${{\mathbf{\tilde{x}}}_{{{I}_{j}},t}}$ for the same time span of 2 seconds. For the AM-FM signals the same procedure is applied for $k=11$ to $k=19$ and their symmetric counterparts (left, bottom panel) where we show the high fluctuations of the first plateau for the reconstruction of $Y_1(t)$, being negligible those linked to the reconstruction of $Y_2(t)$; and $k=33$ to $k=41$ (right, bottom panel) and, again, their symmetric counterparts, where the variation in $Y_2(t)$ is shown\footnote{The reconstruction of the trend follows the same methodology with $k=0$ and exactly reproduces the generated ones. Results are available from the authors upon request.}.
All in all, we have seen that M-CiSSA is able to capture latent components of very different nature (non-stationary, common, idiosyncratic, modulated in amplitude, with time varying frequency), that can be out of phase, reconstructing the sources of variation in a multivariate context. The introduced circulant matrices of second moments allows the match between frequencies and eigenvalues, providing an automated identification of the extracted signals.
We now apply M-CiSSA to the Primary Commodity Energy Prices published by the International Monetary Fund (IMF). Energy commodities accounted for 40.9% of the world trade between 2014 and 2016 and are central to competitiveness in industry, notably influencing consumers by affecting their energy consumption patterns and total expenditure. The evolution of energy prices impacts on other non-energy primary commodities (Kratschell and Schmidt, 2017), exchange rates (Xu et al., 2019), and inflation (Garratt and Petrella, 2019), among others. Energy prices not only affect economic activity but they are also influenced by it (as Alquist et al., 2019, and Kilian and Zhou, 2018, show for commodity prices in general). Also real price shifts may affect the commodity demand and result in one of the causes of fluctuations of the business cycle. This relationship suggests that besides analyzing the long-run behaviour for policy issues, studying the cyclical frequency of energy prices is also of paramount importance. A second issue is related to decoupling among different energy prices, as on a theoretical and empirical basis, oil and natural gas are close substitutes in the long run. In this regard, some authors have found that oil prices drove US gas prices, and also that oil prices led the co-movement between European and North American natural gas prices (see, e.g., Brown and Yucel, 2008, 2009). Despite this evidence, US oil and gas prices seem to have decoupled since 2009 triggering a new discussion on whether or not this is permanent (see, e.g., Erdos, 2012; Nick and Thoenes, 2014; Zhang and Ji, 2018). Since target policies should differ depending on whether decoupling amongst markets occurs in the long run or at medium frequency, we have tackled this problem by addressing the particular frequencies at which markets might be decoupled.
We have applied M-CiSSA to the multivariate analysis of the monthly Primary Commodity Energy Prices included in the category ENERGY by the IMF. We have analyzed the sample that goes from January 1992 until December 2022 ($T=372$ observations). ENERGY comprises a set of nine commodity prices: Australian and South African Coal (COALAU, COALSA); Brent, Dubai and West Texas Intermediate Crude Oil (OILBRE, OILDUB, OILWTI); European, Indonesian and US Natural Gas (NGASEU, NGASJP, NGASUS) and Propane (PROPANE). The original prices in US\$ have been deflated by the US Consumer Price Index to turn them into real terms\footnote{Consideration of the data directly in nominal terms yields the same conclusions. Results are available upon request.} and they have been transformed into index numbers (2016=100) for homogeneity. Figure (ref) shows their graphs.
Before applying CiSSA and M-CiSSA, we only need to choose one parameter, the window length $L$. Due to the monthly periodicity, the first consideration is that $L$ should be multiple of 12. On the other hand, $L$ must also be a multiple of the cycle periods to be analyzed as recommended by Golyandina and Zhigljavsky (2013). In this regard, as economic business cycles expand approximately every 8 years, we have chosen $L=12 \times 8 =96$ and built the trajectory matrix and the associated block circulant matrix comprising 96 blocks of size $9 \times 9$. After diagonalizing this matrix, each resulting block is associated with a frequency $w_k=\frac{k-1}{L}$ and within each block we can further diagonalize it and understand the fluctuations within each frequency. In particular, $k=1$ will correspond to the trend, $k=2$ to $6$ and their symmetric counterparts $96$ to $92$, to the business cycle.
Table (ref) shows the accumulated contribution to the total variability of each frequency, computed as the ratio in ((ref)) but summing all the eigenvalues corresponding to the same block over the total sum of eigenvalues. The accumulated contribution shows that the trend is the most informative signal capturing 36.7% of the total variability. Furthermore, 96-month (8 year) cycles explain 22.3% of the total variability. The 48-month (4 year) cycles account for 16.0% of the variability. Therefore, by analyzing trends and 8 and 4-year cycles, we can explain 75.3% of the variability of energy primary commodity prices.
Information about the trend is represented by the first $9 \times 9$ block of the approximation to the cross spectral density. Remember that, according to the uniqueness theorem (Theorem (ref)), the univariate estimation in CiSSA of an oscillatory component at each frequency is the sum of the 9 subcomponents associated with this frequency in the multivariate M-CiSSA setup. If the number of subcomponents required to pick up practically the total variability of the long-run behaviour of the individual series, is less than the total number of series, then there are common long-run trends. Disentangling what is common and what is idiosyncratic to each series at each frequency is one of the contributions of M-CiSSA. Table (ref) shows that the variability of the trend is mainly explained by the first three eigenvalues of this first block, as they account for 99.7%.
Figure (ref) shows the estimated trends for every energy commodity price and the reconstruction from the three first subcomponents. There is great similarity between the univariate trend and the multivariate one reconstructed by the 3 subcomponents, therefore understanding the flutctuations implied by the eigenvectors will characterize the long-run drivers of the 9 series.
Looking at the individual series, we can see in the top panel of Table (ref) that the 3 subcomponents explain almost all the variability of the trend of each variable, being enough to reconstruct the long run behaviour of each series (the minimum variability explained is 98.7% for NGASJP).
Additionally, the analysis of the eigenvectors also helps to understand the construction of the common forces of the trend. From equation ((ref)), the eigenvectors corresponding to the first block $k=1$ are just vectors of ones multiplied by their own constant, and they capture the changing level of the series. Each eigenvector can be subdivided in the subcomponents associated to each time series. Table (ref) shows the relative weight of each variable\footnote{The relative weights are computed as $100\times$ the sum of the squares of the components of the eigenvector associated to each variable. Recall that eigenvectors have modulus 1.} in the first, second and third eigenvectors. The loadings\footnote{For brevity, the loadings are not shown but are available from the authors.} for the first eigenvector are all positive with a small contribution of NGASUS that becomes the main driver in the second eigenvector. Therefore, our first conclusion is that the long-run behaviour of NGASUS is decoupled from the rest of the energy prices. The third eigenvector gives positive weight to all oils and negative to coals and natural gases in Europe and the USA, therefore, it separates the oil market from the coal and natural gas markets.
The second and third components in terms of relevance are the 96-month (8-year) and 48-month (4-year) cycles that account for the 22.6% and 16.0% of the total variability (Table (ref)).
Regarding the 96-month cycle, it is represented by its corresponding $9 \times 9$ block for $k=2$ and $k=96$ of the cross power spectral density. Table (ref) shows that the variability within this frequency can be approximated (98.1%) by the sum of the first 3 eigenvectors. Table (ref) shows the high approximation, over 95.4% in all the series, of the sum of these 3 subcomponents. Therefore, understanding the formation of the corresponding subcomponents allows to understand the common drivers of the 96-month cycle.
Table (ref) shows the relative weights of each of the three eigenvectors. It can be seen that almost half of the total weight corresponds to NGASEU in the first eigenvector, while in the second one, it is NGASUS the main driver. In the third eigenvector oil (specially OILBRE and OILDUB) and coal (COALAU and COALSA) gain more relevance compared to the previous subcomponents.
The left panel of Figure (ref) shows the segments of the first, second and third eigenvectors corresponding to each commodity. The conclusions of Table (ref) still hold at the look of the amplitude of the waves of the different eigenvector segments, however the analysis can be extended to understand the nature of the different cycles. Regarding the first subcomponent (top-left graph) that explains (Table (ref)) 66.8% of the 96-month variability, we can see that all the segments share the same minima and maxima. However, differences appear in the second and third subcomponents that jointly explain more than 30% of the variability. The most striking effect is the decoupling in both graphs (left panel, middle and bottom graphs) of NGASUS (light blue line) as it has its own cycle in the second and third graphs, corresponding to additional subcomponents, different from the bulk of cycles that we see in the first graph (top panel, left column).
A similar analysis can be performed to understand the formation of the 4-year cycle. Again, 3 subcomponents are needed in order to explain almost all the variability, $98.3\%$ (See Table (ref)). And for the 4 year cycle associated to each variable, the sum of these 3 subcomponents explain more than $95.9\%$ (see Table (ref)). As it can be seen in Table (ref), the first subcomponent in the 4-year cycle is dominated by NGASEU, while the second one is mixture of NGASUS and PROPANE. Finally, the last subcomponent is mainly dominated only by NGASUS. (Figure (ref), right panel), we see that all prices are aligned for the first set of subvectors (top graph), showing NGASEU more amplitude in the cyclical 4-year wave for the first subcomponent, while in the second subcomponent (right panel, middle graph) PROPANE and NGASUS are separated from NGASEU. Finally, the third subcomponent (right panel, bottom graph) NGASUS presents a specific movement opposed to the remaining variables.
All in all, we have seen how 8 and 4-year cycles can be explained in a multivariate setup with M-CiSSA.
We have introduced a new nonparametric methodology, M-CiSSA, that enables us to identify common fluctuations within a group of time series by frequency. It is also useful to deseasonalize and denoise the extracted signals. M-CiSSA is based on the properties of block circulant matrices and its diagonalization that links the eigen-structure of this matrix and the multivariate spectral density at different frequencies. Further diagonalization within each block of the spectral density matrix enables us to understand the multivariate information contained for each frequency.
We have proved that the uniqueness of each extracted component per frequency, that is, the extracted components in a univariate fashion coincide with the sum of all the subcomponents for each frequency extracted in a multivariate way. The added value of the multivariate approach is the ability to disentangle what is common from and what is idiosyncratic. This cannot be done using the univariate approach because there is no information regarding what is common. We depart from the abundant literature of dynamic factor models in several ways: firstly, we are able to separate common and idiosyncratic fluctuations by frequency; secondly, we can also discover phase shifts among the extracted signals for each series. Additionally, the multivariate approach is useful for further denoise the univariate extracted signals. We have illustrated the ability of M-CiSSA to separate different underlying signals through a synthetic example. We have also applied M-CiSSA to Primary Energy Commodity Prices. On the one hand, understanding the formation of Primary Energy Commodity Prices is key for medium and long-term energy policy issues. In particular, we have discovered the nature of decoupling at different frequencies and that markets can have a very different behaviour in the medium and long run. This is precisely one novelty of our approach: we can ascertain at what frequencies markets are decoupled. The comparison between the results in the univariate and multivariate cases also enables us to discover which prices are held or sustained by forces outside the global evolution of the markets. In particular, in the long run, we found that although coals are not decoupled, their prices are higher than suggested by the global forces of the markets. We also found that the long-run prices for natural gas in Europe exhibit a lower level than suggested, probably due to the European policy of installing liquid natural gas reserves.
Finally, our approach is quite general and can be applied to many disciplines involving time series analysis. Among other things, it can be used to separate long, medium, and short-run analysis, cyclical analysis, forecasting scenarios, deseasonalizing and denoising time series, and extracting the common from the idiosyncratic signals by frequency due to the ability of the procedure to extract signals at the desired frequencies specified by the user.