Functional Singular Spectrum Analysis

TL;DR

This paper introduces functional SSA, a novel method for analyzing functional time series by combining functional data analysis with univariate SSA, demonstrated through simulations and real-world applications.

Contribution

The paper presents a new functional SSA method that extends traditional SSA to functional data, offering an alternative to MSSA and dFPCA with practical implementation tools.

Findings

Functional SSA outperforms MSSA in certain scenarios.

The method provides a competitive alternative to dFPCA.

An R package and web app facilitate practical use.

Abstract

In this paper, we introduce a new extension of the Singular Spectrum Analysis (SSA) called functional SSA to analyze functional time series. The new methodology is developed by integrating ideas from functional data analysis and univariate SSA. We explore the advantages of the functional SSA in terms of simulation results and two real data applications. We compare the proposed approach with Multivariate SSA (MSSA) and dynamic Functional Principal Component Analysis (dFPCA). The results suggest that further improvement to MSSA is possible, and the new method provides an attractive alternative to the dFPCA approach that is used for analyzing correlated functions. We implement the proposed technique to an application of remote sensing data and a call center dataset. We have also developed an efficient and user-friendly R package and a shiny web application to allow interactive exploration of the results.

Tables

Table 1. Table 1 : The mean of RMSE for 1000 generation of the simulated model by dFPCA, MSSA and FSSA approaches.

Model	$\omega$	N	dFPCA		MSSA				FSSA
				$L=20$	$L=30$	$L=40$		$L=20$	$L=30$	$L=40$
Setting 1 GWN	0.00	50	0.092	0.026	0.022	0.019		0.010	0.011	0.014
		100	0.088	0.024	0.020	0.018		0.008	0.007	0.007
		150	0.074	0.024	0.020	0.017		0.007	0.006	0.006
		200	0.071	0.024	0.019	0.017		0.007	0.006	0.006
	0.10	50	0.062	0.028	0.024	0.022		0.009	0.011	0.014
		100	0.045	0.027	0.023	0.020		0.006	0.006	0.007
		150	0.037	0.027	0.022	0.019		0.005	0.005	0.005
		200	0.033	0.026	0.022	0.019		0.005	0.005	0.005
	0.25	50	0.046	0.028	0.024	0.022		0.009	0.011	0.014
		100	0.038	0.027	0.022	0.020		0.006	0.006	0.007
		150	0.032	0.027	0.022	0.019		0.005	0.005	0.005
		200	0.030	0.026	0.022	0.019		0.005	0.005	0.005
Setting 2 ${\Vert \Psi \Vert }_{𝒮}^2=0$	0.00	50	1.008	0.251	0.246	0.270		0.244	0.248	0.284
		100	0.993	0.205	0.185	0.179		0.190	0.176	0.174
		150	0.974	0.192	0.166	0.153		0.175	0.155	0.145
		200	0.960	0.186	0.158	0.143		0.168	0.145	0.134
	0.10	50	0.636	0.220	0.211	0.235		0.204	0.216	0.264
		100	0.633	0.193	0.165	0.154		0.158	0.144	0.143
		150	0.633	0.188	0.158	0.141		0.149	0.130	0.120
		200	0.632	0.186	0.154	0.136		0.144	0.123	0.112
	0.25	50	0.636	0.223	0.214	0.235		0.206	0.218	0.263
		100	0.633	0.194	0.167	0.156		0.160	0.146	0.144
		150	0.633	0.188	0.157	0.140		0.148	0.129	0.120
		200	0.632	0.185	0.153	0.136		0.143	0.122	0.111
Setting 3 ${\Vert \Psi \Vert }_{𝒮}^2=0.5$	0.00	50	0.841	0.316	0.310	0.350		0.305	0.311	0.364
		100	0.818	0.271	0.240	0.231		0.250	0.228	0.224
		150	0.806	0.262	0.223	0.203		0.238	0.207	0.192
		200	0.801	0.258	0.216	0.194		0.233	0.198	0.180
	0.10	50	0.696	0.277	0.267	0.300		0.255	0.270	0.330
		100	0.709	0.241	0.207	0.193		0.198	0.180	0.178
		150	0.710	0.235	0.197	0.176		0.187	0.162	0.151
		200	0.710	0.232	0.192	0.170		0.180	0.153	0.139
	0.25	50	0.707	0.214	0.205	0.225		0.198	0.209	0.254
		100	0.707	0.187	0.160	0.149		0.153	0.139	0.138
		150	0.705	0.181	0.151	0.135		0.141	0.123	0.114
		200	0.705	0.178	0.148	0.130		0.136	0.116	0.106
Setting 4 ${\Vert \Psi \Vert }_{𝒮}^2=0.9$	0.00	50	0.926	0.486	0.477	0.535		0.464	0.476	0.553
		100	0.915	0.441	0.391	0.373		0.404	0.368	0.358
		150	0.901	0.432	0.371	0.339		0.391	0.341	0.318
		200	0.823	0.432	0.363	0.326		0.388	0.331	0.301
	0.10	50	0.830	0.341	0.330	0.384		0.313	0.331	0.413
		100	0.831	0.286	0.245	0.229		0.235	0.212	0.210
		150	0.832	0.279	0.233	0.208		0.220	0.190	0.176
		200	0.829	0.275	0.227	0.200		0.213	0.179	0.163
	0.25	50	0.829	0.210	0.200	0.226		0.195	0.205	0.255
		100	0.829	0.175	0.149	0.138		0.140	0.129	0.127
		150	0.829	0.169	0.141	0.125		0.129	0.112	0.105
		200	0.833	0.167	0.138	0.121		0.124	0.106	0.097

Table 2. Table 2 : P-values of the multivariate trace periodicity test of Hörmann et al., ( 2018 ) and the bootstrap trend test (Grosjean and Ibanez,, 2018 ) on four FTS ( y t ​ ( s ) subscript 𝑦 𝑡 𝑠 y_{t}(s) , R 1 subscript 𝑅 1 R_{1} , R 2 subscript 𝑅 2 R_{2} , R 3 subscript 𝑅 3 R_{3} ).

FTS	periodicity test	trend test
$y_t\left(s\right)$	0.00	0.02
$R_1$	0.00	0.03
$R_2$	0.45	0.03
$R_3$	0.30	0.07

Table 3. Table 3 : P-values of the multivariate trace periodicity test of Hörmann et al., ( 2018 ) on six sets of FTS ( y t ​ ( s ) subscript 𝑦 𝑡 𝑠 y_{t}(s) , R 1 subscript 𝑅 1 R_{1} , R 2 subscript 𝑅 2 R_{2} , R 3 subscript 𝑅 3 R_{3} , R 4 subscript 𝑅 4 R_{4} , R 5 subscript 𝑅 5 R_{5} ) for the periods of length 7 7 7 and 30 30 30 days.

FTS	d=7	d=30
$y_t\left(s\right)$	0	0.86
$R_1$	0	0.96
$R_2$	0	0.97
$R_3$	0	0.30
$R_4$	1	0.00
$R_5$	1	0.17

Equations

$$\mathbf{X}=\mathcal{T}\left({\mathbf{y}}_N\right)=\left[{\mathbf{x}}_1,\dots ,{\mathbf{x}}_K\right],$$

$$\mathbf{X}=\sum _{i=1}^L\sqrt{{\lambda }_i}{\mathbf{u}}_i{\mathbf{v}}_i^{\top }=\sum _{i=1}^L{\mathbf{X}}_i.$$

$${\mathbf{X}}_i={\mathbf{u}}_i{\mathbf{u}}_i^{\top }\mathbf{X}=\left[({\mathbf{u}}_i\otimes {\mathbf{u}}_i){\mathbf{x}}_1,\dots ,({\mathbf{u}}_i\otimes {\mathbf{u}}_i){\mathbf{x}}_K\right],$$

$$\mathbf{X}={\mathbf{X}}_{I_1}+{\mathbf{X}}_{I_2}+\cdots +{\mathbf{X}}_{I_m}.$$

$${\mathbf{y}}_N={\tilde{\mathbf{y}}}_1+{\tilde{\mathbf{y}}}_2+\dots +{\tilde{\mathbf{y}}}_m,$$

$${\rho }^{(w)}({\text{y}}_N^{(1)},{\text{y}}_N^{(2)})=\frac{{⟨{\text{y}}_N^{(1)},{\text{y}}_N^{(2)}⟩}_w}{\sqrt{{⟨{\text{y}}_N^{(1)},{\text{y}}_N^{(1)}⟩}_w}\sqrt{{⟨{\text{y}}_N^{(2)},{\text{y}}_N^{(2)}⟩}_w}},$$

$$\text{mathbfcal}Z\mathbf{a}=\left(\begin{array}{c}{\sum }_{j=1}^K{a}_j{z}_{1,j}\\ ⋮\\ {\sum }_{j=1}^K{a}_j{z}_{L,j}\end{array}\right),z_{i,j}\in \mathbb{H},\text{and}\mathbf{a}=(a_1,\dots ,a_K)\in {\mathbb{R}}^K.$$

$$⟨\text{mathbfcal}{Z}_1,\text{mathbfcal}{Z}_2{⟩}_{\mathcal{F}}:=\sum _{i=1}^L\sum _{j=1}^K⟨z_{i,j}^{(1)},z_{i,j}^{(2)}⟩.$$

$$\sum _{i=1}^N\Vert x_i-\bar{x}{\Vert }^2\le \sum _{i=1}^N\Vert x_i-y{\Vert }^2,$$

$${\mathbf{x}}_j(s):={\left(\begin{array}{c}{y}_j(s),y_{j+1}(s),\dots ,y_{j+L-1}(s)\end{array}\right)}^{\top },j=1,\dots ,K,$$

$$\text{mathbfcal}X\mathbf{a}:=\sum _{j=1}^K{a}_j{\mathbf{x}}_j=\left(\begin{array}{c}{\sum }_{j=1}^K{a}_j{y}_j\\ {\sum }_{j=1}^K{a}_j{y}_{j+1}\\ ⋮\\ {\sum }_{j=1}^K{a}_j{y}_{j+L-1}\end{array}\right),\mathbf{a}={\left(a_1,\dots ,a_K\right)}^{\top }\in {\mathbb{R}}^K.$$

$$\mathbf{X}(s)=\left[\begin{array}{c}{\mathbf{x}}_1(s),\dots ,{\mathbf{x}}_K(s)\end{array}\right].$$

$$\text{mathbfcal}{X}^{\ast }\mathbf{z}=\left(\begin{array}{c}{\sum }_{i=1}^L⟨y_i,z_i⟩\\ {\sum }_{i=1}^L⟨y_{i+1},z_i⟩\\ ⋮\\ {\sum }_{i=1}^L⟨y_{i+K-1},z_i⟩\end{array}\right),\mathbf{z}={\left(z_1,\dots ,z_L\right)}^{\top }\in {\mathbb{H}}^L,$$

$$\text{mathbfcal}S\mathbf{z}$$

$$c_{i,j}(s,u):=\sum _{k=1}^K{y}_{i+k-1}(s)y_{j+k-1}(u),\text{for}s,u\in [0,1].$$

$$\text{mathbfcal}S\mathbf{z}(u)$$

$$\text{mathbfcal}S{\mathbf{\psi }}_i={\lambda }_i{\mathbf{\psi }}_i,\text{and}{\lambda }_i⟶0\text{as}i⟶\infty .$$

$$\text{mathbfcal}S=\sum _{i=1}^{\infty }{\lambda }_i{\mathbf{\psi }}_i\otimes {\mathbf{\psi }}_i.$$

$$\text{C}(s,u)=\sum _{i=1}^{\infty }{\lambda }_i{\mathbf{\psi }}_i(s){\mathbf{\psi }}_i^{\top }(u).$$

$$\text{mathbfcal}{X}_i\mathbf{a}:=\sum _{j=1}^K{a}_j({\mathbf{\psi }}_i\otimes {\mathbf{\psi }}_i){\mathbf{x}}_j=({\mathbf{\psi }}_i\otimes {\mathbf{\psi }}_i)\sum _{j=1}^K{a}_j{\mathbf{x}}_j.$$

$${\mathbf{X}}_i(s)$$

$$=\left[\begin{array}{c}({\mathbf{\psi }}_i\otimes {\mathbf{\psi }}_i){\mathbf{x}}_1(s),\dots ,({\mathbf{\psi }}_i\otimes {\mathbf{\psi }}_i){\mathbf{x}}_K(s)\end{array}\right].$$

$$\text{mathbfcal}X=\sum _{i=1}^{\infty }\text{mathbfcal}{X}_i.$$

$${\mathcal{S}}^{†}:=\text{mathbfcal}{X}^{\ast }\text{mathbfcal}X=\left[\begin{array}{ccc}{\sum }_{i=1}^L⟨y_i,y_i⟩& \cdots & {\sum }_{i=1}^L⟨y_i,y_{i+K-1}⟩\\ ⋮& \ddots & ⋮\\ {\sum }_{i=1}^L⟨y_{i+K-1},y_i⟩& \cdots & {\sum }_{i=1}^L⟨y_{i+K-1},y_{i+K-1}⟩\end{array}\right].$$

$$\text{mathbfcal}X=\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\text{v}}_i\otimes {\mathbf{\psi }}_i,$$

$$\text{mathbfcal}X\mathbf{a}=\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}⟨{\text{v}}_i,\mathbf{a}{⟩}_{{\mathbb{R}}^K}{\mathbf{\psi }}_i,$$

$$\text{mathbfcal}X=\text{mathbfcal}{X}_{I_1}+\text{mathbfcal}{X}_{I_2}+\cdots +\text{mathbfcal}{X}_{I_m}.$$

$$\Vert \text{mathbfcal}{X}_{I_q}-\text{mathbfcal}{\tilde{X}}_{I_q}{\Vert }_{\mathcal{F}}^2\le \Vert \text{mathbfcal}{X}_{I_q}-\widetilde{\text{mathbfcal}Z}{\Vert }_{\mathcal{F}}^2,\text{for any}\widetilde{\text{mathbfcal}Z}\in {\mathbb{H}}_H^{L\times K}.$$

$${\tilde{x}}_{i,j}^q=\frac{1}{n_s}\sum _{(l,k):l+k=s}{x}_{l,k}^q,$$

$${⟨{\text{y}}_N^{(1)},{\text{y}}_N^{(2)}⟩}_w=\sum _{i=1}^N{w}_i⟨y_i^{(1)},y_i^{(2)}⟩,$$

$${⟨{\text{y}}_N^{(1)},{\text{y}}_N^{(2)}⟩}_w=0.$$

Full text

Keywords: Functional Time Series, Hilbert Space, Singular Spectrum Analysis, SVD

1 Introduction

One of the popular approaches in the decomposition of time series is accomplished using the rates of change. In this approach, the observed time series is partitioned (decomposed) into informative trends plus potential seasonal (cyclical) and noise (irregular) components. Aligned with this principle, Singular Spectrum Analysis (SSA) is a model-free procedure that is commonly used as a nonparametric technique in analyzing the time series. SSA is intrinsically motivated as an exploratory and model building tool rather than a confirmatory procedure (Golyandina et al.,, 2001). SSA does not require restrictive assumptions such as stationarity, linearity, and normality. It can be used for a wide range of purposes such as trend and periodic component detection and extraction, smoothing, forecasting, change-point detection, gap filling, causality and so on; (see, e.g. Golyandina et al.,, 2001; Moskvina and Zhigljavsky,, 2003; Kondrashov et al.,, 2010; Golyandina and Osipov,, 2007; Mohammad and Nishida,, 2011; Mahmoudvand and Rodrigues,, 2016; Rodrigues and Mahmoudvand,, 2016).

The implementation of SSA over time series is similar to that of Principal Components Analysis (PCA) of multivariate data. In contrast to PCA, which is applied to a data matrix with independent rows, SSA is applied to a time series. It provides a representation of the given time series in terms of eigentriples of a so-called trajectory matrix (Alexandrov,, 2009).

Up to this day, many studies have been published with extensions and applications of SSA. Extensions to a multivariate model as well as to a two-dimensional setting can be found, e.g., in Golyandina and Zhigljavsky, (2013); Golyandina et al., (2018); Hassani and Mahmoudvand, (2018) and references therein. In the regular SSA, we assume that the observation at each time point is scalar, vector or array. As a matter of interest, one may consider a series of curves observed over time, and use the basics of Hilbert space in the functional data analysis (FDA) framework to introduce the concept of functional SSA (FSSA).

While the research in FDA has grown extensively in recent years, there have been relatively few contributions dealing with functional time series (FTS); see, e.g., Hörmann and Kokoszka, (2012) and Bosq, (2000). Although most of the current FTS approaches focus on a parametric fit for inferences and forecasting, there exist some other approaches in the literature that extend FPCA to incorporate the temporal correlation of FTS. For instance, Hörmann et al., (2015) introduced dynamic FPCA (dFPCA) to analyze FTS. This approach assumes the strong assumption of stationarity which is not generally held in practice. It would be of interest to non-parametrically decompose a nonstationary FTS to reveal the respective trends plus seasonal and irregular components in an appropriate manner. Consistent with this approach, and as a first step, Fraiman et al., (2014) introduced a new concept of trends for the FTS. Furthermore, Hörmann et al., (2018) considered the periodic components for the FTS and derived several procedures to test the periodicity using frequency domain analysis. To the best of our knowledge, existing studies mainly focus on detecting rather than extracting interpretable components.

Since one of the primary missions of SSA is to extract trends and periodic components of a regular (non-functional) time series, it would be rational to establish a similar elegant nonparametric procedure to extract such components in FTS. In this paper we use the basics of SSA and multivariate functional PCA (MFPCA), introduced in Happ and Greven, (2018); Chiou et al., (2014), to develop FSSA. In a nutshell, the core of SSA is to use PCA on the variables being lagged versions of a single time series. Since a lagged vectors of FTS forms a multivariate functional variable, we use the theory of MFPCA to develop the FSSA procedure. The new methodology, FSSA, not only can serve as a nonparametric dimension reduction tool to decompose the functional time series; it can also be used as a visualization tool to illustrate the concept of seasonality and periodicity in the functional space over time.

In order to depict the idea of our approach and to show its utility, consider the following motivating example involving a real dataset which is described in detail in the supplementary materials. This data provides the intraday number of calls to a call center, during different times of the days for one year. The associated $365$ curves is represented in an overlapping pattern in Figure 1 (left). In Figure 1 (right) we investigate the pattern among weekdays and weekend days. As we can see, the intraday patterns of weekends (Friday and Saturday) are significantly different from workdays while workdays seem to have similar patterns. Investigators used variants of FPCA to analyze the call center data in literature (Shen and Huang,, 2005; Huang et al.,, 2008; Maadooliat et al.,, 2015). For illustration purpose, we compare the results of the proposed method (FSSA) and dFPCA on this dataset. Figure 2 (top) presents the projection of the data into the first four FPCs of the dFPCA obtained from the freqdom.fda package in R (Hörmann et al.,, 2015). We used seven different colors to differentiate between different days of a week. As one may observe, visually, there is no clear separation in either one of the FPC graphs in the top row. Although this may not be surprising, as the purpose of PCA of any type is to reduce the dimensionality, and not necessary decompose the data into regular trends, periodic and irregular components. In contrast, the grouping results that we obtained using the FSSA on the call center data are given in Figure 2 (bottom). It can be seen that the functional behavior of seven days of a week can be well-distinguished, visually, using either one of the last two groups (groups 3 and 4).

The rest of the paper is organized as follows. Section 2 reviews the core of SSA for completeness. Section  3 presents the theoretical foundations and some properties of the proposed method (FSSA), and Section 4 provides implementation details. Section 5.1 reports simulation results to illustrate the use of the proposed approach in analyzing FTS, and to compare it with MSSA and dFPCA. Application to a real data example on remote sensing is given in Section 5.2. Section 6 provides some discussions and concluding remarks.

2 General scheme of SSA

As we mentioned in Section 1, SSA can be used for many purposes. However, as we intend to introduce the functional version of SSA for decomposing FTS, we review a general scheme of SSA to perform time series decomposition in this section.

2.1 Univariate SSA

Throughout this section, we consider $y_{i}$’s are elements of Euclidean space $\mathbb{R}$. Suppose that ${\bf y}_{N}=(y_{1},y_{2},\ldots,y_{N})^{\top}$ is a realization of size $N$ from a time series. The basic SSA algorithm consists of four steps: Embedding, Decomposition, Grouping, and Reconstruction.

Step 1. Embedding

This step generates a multivariate object by tracking a moving window of size $L$ over the original time series, where $L$ is called window length parameter and $1<L<N$. Embedding can be regarded as a mapping operator $\mathcal{T}$ that transfers the series ${\bf y}_{N}$ into a so-called trajectory matrix ${\bf X}$ of dimension $L\times K$, defines by

[TABLE]

where $K=N-L+1$ and ${\boldsymbol{x}}_{j}=(y_{j},y_{j+1},\ldots,y_{j+L-1})^{\top}$, for $j=1,\ldots,K$, are called lagged vectors. Note that the trajectory matrix ${\bf X}$, is a Hankel matrix, which means that all the elements along the anti-diagonals are equal. Indeed, the embedding operator $\mathcal{T}$ is a one-to-one mapping from $\mathbb{R}^{N}$ into $\mathbb{R}^{L\times K}_{H}\subseteq\mathbb{R}^{L\times K}$, where $\mathbb{R}^{L\times K}_{H}$ is the set of all $L\times K$ Hankel matrices.

Step 2. Decomposition

In this step, the singular value decomposition (SVD) for the trajectory matrix is computed as:

[TABLE]

where $\sqrt{\lambda_{i}}$ is the $i^{th}$ singular value of $\bf X,\ {\bf u}_{i}$ and {\bf v}_{i}\ are the associated (orthonormal) left and right singular vectors, and ${\bf X}_{i}=\sqrt{\lambda_{i}}{\bf u}_{i}{\bf v}_{i}^{\top}$ is called the respective elementary matrix. Note that ${\bf u}_{i}$ is an eigenvector of ${\bf XX}^{\top}$ corresponding to the eigenvalue $\lambda_{i}$. Moreover, it yields that

[TABLE]

where, in this section, $\otimes$ denotes the outer (tensor) product of two vectors.

Step 3. Grouping

Consider a partition of the set of indices $\{1,2,\ldots,r\}$, where $r$ is the rank of the matrix $\bf{X}$, into $m$ disjoint subsets $\{{I_{1}},{I_{2}},\ldots,{I_{m}}\}$. For any positive integer $q$, i.e. $1\leq q\leq m$, the matrix ${\bf X}_{I_{q}}$ is defined as ${\bf X}_{I_{q}}=\sum_{i\in{I_{q}}}{\bf X}_{i}$. Thus, by the expansion (2) we have the *grouped matrix decomposition *

[TABLE]

Each group in (4) should correspond to a component in time series decomposition. These components can be considered as trend, cycle, seasonal, noise, etc.

Step 4 Reconstruction

Finally, the resulting matrices ${\bf X}_{I_{q}}$ in (4), are transformed back into the form of the original series ${\bf y}_{N}$ by an inverse operator $\mathcal{T}^{-1}$. In order to do this, first, it is necessary that each matrix ${\bf X}_{I_{q}}$ to be approximated by a matrix in $\mathbb{R}^{L\times K}_{H}$. This approximation is performed optimally in the sense of orthogonal projection of ${\bf X}_{I_{q}}$ on $\mathbb{R}^{L\times K}_{H}$ with respect to the Frobenius norm. Denote this projection by $\Pi:\mathbb{R}^{L\times K}\rightarrow\mathbb{R}^{L\times K}_{H}$. It is shown that the projection $\Pi$ is the averaging of the matrix elements over the antidiagonals $i+j=const$. By combining the results of this step and (4), we obtain the final decomposition of the series in the form of

[TABLE]

where $\tilde{\bf y}_{q}=\mathcal{T}^{-1}\Pi({\bf X}_{I_{q}})$, for $q=1,\ldots,m$.

The above algorithm can be extended to perform Multivariate SSA (MSSA) for analyzing multivariate time series. The only difference is in defining the trajectory matrix which can be defined by stacking univariate trajectory matrices horizontally or vertically (Hassani and Mahmoudvand,, 2018).

It is well known that SSA does not require restrictive assumptions; however, it is ideal to have a time series with separable components. Therefore, we present tools to measure the separability of components in the next subsection.

2.2 Separability

Let $\textbf{y}_{N}^{(i)}=\left\{y_{1}^{(i)},\ldots,y_{N}^{(i)}\right\}$, for $i=1,2$ be two time series and consider an additive model as $\textbf{y}_{N}=\textbf{y}_{N}^{(1)}+\textbf{y}_{N}^{(2)}$. The series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ are called separable when each lagged vector of $\textbf{y}_{N}^{(1)}$ is orthogonal to the lagged vectors of $\textbf{y}_{N}^{(2)}$. To measure the degree of separability between two time series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$, Golyandina et al., (2001) introduced the so-called w-correlation

[TABLE]

where, $\left\langle\textbf{y}_{N}^{(t)},\textbf{y}_{N}^{(s)}\right\rangle_{w}=\sum_{i=1}^{N}w_{i}y_{i}^{(t)}y_{i}^{(s)}$ for $t,s=1,2$ and $w_{i}=\min\{i,L,N-i+1\}$. We call two series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ w-orthogonal if $\rho^{(w)}(\textbf{y}_{N}^{(1)},\textbf{y}_{N}^{(2)})=0$ for approriate values of $L$ (see the next subsection for more details). Note that $\tilde{\bf y}_{q}$, $q=1,\ldots,m$, is the reconstructed component produced by the group $I_{q}$, and the matrix of {\bm{\rho}}^{(w)}=\bigl{\{}\rho^{(w)}(\tilde{\bf y}_{i},\tilde{\bf y}_{j})\bigr{\}}_{i,j=1}^{m} is called w-correlation matrix.

2.3 Parameter Selection

There are two basic parameters in SSA procedure; window length ($L$) and grouping parameters. Choosing improper values for these parameters yields an incomplete reconstruction and misleading results in subsequent analysis. In spite of the importance of choosing $L$ and grouping parameters for SSA, no ideal solution has been yet proposed. A thorough review of the problem shows that the optimal choice of the parameters depends on the intrinsic structure of the data and the purposes of the study (Golyandina et al.,, 2001; Golyandina and Zhigljavsky,, 2013). However, there are several recommendations and rules that work well for a wide range of scenarios. It is recommended to select the window length parameter, $L$, to be a large integer that is multiple of the periodicities of the time series, but not larger than $\frac{N}{2}$.

In addition, there are several utilities for effective grouping. These tools include analyzing the periodogram, paired plot of the singular vectors, scree plot of the singular values, and w-correlation plot; see Golyandina et al., (2001) for more details.

3 Theoretical Foundations of FSSA

We start this section with the mathematical foundations that are used to develop the functional SSA procedure. From hereafter, we consider $\textbf{y}_{N}=\left(y_{1},\ldots,y_{N}\right)^{\top}$ is a FTS of length $N$. This means that each elements $y_{i}:[0,1]\rightarrow\mathbb{R}$ belongs to $\mathbb{H}:=\mathcal{L}^{2}([0,1])$, the space of square integrable real functions defined on the interval $[0,1]$. Here, the space $\mathbb{H}$ is a Hilbert space, equipped with inner product $\left\langle x,y\right\rangle=\int_{0}^{1}x(s)y(s)ds$. For a given positive integer $k$, the space $\mathbb{H}^{k}$ denotes the Cartesian product of $k$ copies of $\mathbb{H}$; i.e. for an element ${\boldsymbol{x}}\in\mathbb{H}^{k}$, it has the form ${\boldsymbol{x}}(s)=\begin{pmatrix}x_{1}(s),x_{2}(s),\ldots,x_{k}(s)\end{pmatrix}^{\top}$, where $x_{i}\in\mathbb{H}$, and $s\in[0,1]$. Then $\mathbb{H}^{k}$ is a Hilbert space equipped with the inner product $\langle\boldsymbol{x},\boldsymbol{y}\rangle_{\mathbb{H}^{k}}=\sum_{i=1}^{k}\langle x_{i},y_{i}\rangle$. The norms will be denoted by $\|\cdot\|$ and $\|\cdot\|_{\mathbb{H}^{k}}$ in the spaces $\mathbb{H}$ and $\mathbb{H}^{k}$, respectively. For $x\in\mathbb{H}_{1}$, and $y\in\mathbb{H}_{2}$, where $\mathbb{H}_{1}$ and $\mathbb{H}_{2}$ are two Hilbert spaces, we define the tensor(outer) product corresponding to the operator $x\otimes y:\mathbb{H}_{1}\rightarrow\mathbb{H}_{2}$, as ($x\otimes y)h:=\langle x,h\rangle y$, where $h\in\mathbb{H}_{1}$.

For positive integers $L\ \text{and}\ K$, we denote $\mathbb{H}^{L\times K}$ as the linear space spanned by operators $\mathbfcal{Z}:\mathbb{R}^{K}\rightarrow\mathbb{H}^{L}$, specified by $[{z}_{i,j}]_{i=1,\ldots,L}^{j=1,\ldots,K}$ where

[TABLE]

We call an operator $\tilde{\mathbfcal{Z}}=[\tilde{{z}}_{i,j}]\in\mathbb{H}^{L\times K}$ Hankel if $\|\tilde{{z}}_{i,j}-g_{s}\|=0$, for some $g_{s}\in\mathbb{H}$, where $s=i+j$. The space of such Hankel operators will be denoted $\mathbb{H}_{H}^{L\times K}$. For two given operators $\mathbfcal{Z}_{1}=[\textbf{z}^{(1)}_{i,j}]$ and $\mathbfcal{Z}_{2}=[\textbf{z}^{(2)}_{i,j}]$ in $\mathbb{H}^{L\times K}$, define

[TABLE]

It follows immediately that $\langle\cdot,\cdot\rangle_{\mathcal{F}}$, defines an inner product on $\mathbb{H}^{L\times K}$. We will call it Frobenius inner product of two operators in $\mathbb{H}^{L\times K}$. The associated Frobenius norm is $\|\mathbfcal{Z}\|_{\mathcal{F}}=\sqrt{\langle\mathbfcal{Z},\mathbfcal{Z}\rangle_{\mathcal{F}}}$. Before discussing the FSSA algorithm, here we present a lemma that will be used in the last step of the proposed algorithm. Note that the Proofs for all lemmas, theorems and propositions are given in the supplementary materials.

Lemma 3.1**.**

Let $x_{i},\ i=1,\ldots,N$ be elements of the Hilbert space $\mathbb{H}$. If $\bar{x}=\frac{1}{N}\sum_{i=1}^{N}x_{i}$, then

[TABLE]

for all $y\in\mathbb{H}$.

3.1 FSSA algorithm

For an integer $1\leq L\leq{N}/{2}$, let $K=N-L+1$ and define a set of multivariate functional vectors in $\mathbb{H}^{L}$ by

[TABLE]

where ${\boldsymbol{x}}_{j}$’s denote the functional $L-$lagged vectors. The following algorithm provides the FSSA results in four steps.

Step 1. Embedding

Define the operator $\mathbfcal{X}:\mathbb{R}^{K}\rightarrow\mathbb{H}^{L}$ with

[TABLE]

We call $\mathbfcal{X}$ the trajectory operator. It is easy to see that $\mathbfcal{X}=\mathcal{T}\textbf{y}_{N}$, where $\mathcal{T}$ is an operator from $\mathbb{H}^{N}$ to $\mathbb{H}_{H}^{L\times K}$. Evaluating $\mathbfcal{X}\boldsymbol{a}$ at a given point $s\in[0,1]$ is same as the matrix product ${\bf X}(s)\boldsymbol{a}$, where ${\bf X}(s)$ is an $L\times K$ Hankel matrix given by

[TABLE]

Proposition 3.1**.**

The operator $\mathbfcal{X}$ is a bounded linear operator. If we define $\mathbfcal{X}^{*}:\mathbb{H}^{L}\rightarrow\mathbb{R}^{K}$, given by

[TABLE]

then $\mathbfcal{X}^{*}$ is an adjoint operator for $\mathbfcal{X}$.

Step 2. Decomposition

Define the operator $\bm{\mathcal{S}}:\mathbb{H}^{L}\rightarrow\mathbb{H}^{L}$ by $\bm{\mathcal{S}}:=\mathbfcal{X}\mathbfcal{X}^{*}$. Therefore, for given ${\boldsymbol{z}}\in\mathbb{H}^{L}$ it implies that

[TABLE]

Here, the operator $\bm{\mathcal{S}}$ can be also considered as an $L\times L$ matrix with the operator entries $\mathcal{S}_{i,j}:\mathbb{H}\rightarrow\mathbb{H}$, given by $\mathcal{S}_{i,j}=\sum_{l=1}^{K}y_{i+l-1}\otimes y_{j+l-1}$, where $i,j=1,\ldots,L$. Note that, the operator $\mathcal{S}_{i,j}$ defines an integral operator on $\mathbb{H}$, associated to the kernel

[TABLE]

Let us define ${\bf C}:[0,1]\times[0,1]\rightarrow\mathbb{R}^{L\times L}$ to be a kernel matrix with the elements $\{c_{i,j}\}$. Note that ${\bf C}(s,u)={\bf X}(s){\bf X}(u)^{\top}$. It is easy to show that the associated integral operator of ${\bf C}$ is $\mathbfcal{S}$, i.e,

[TABLE]

Proposition 3.2**.**

The operator $\bm{\mathcal{S}}$ defined in (12) is a linear, self-adjoint, positive definite, bounded, continuous and compact operator.

By the results of the Proposition 3.2 and the Hilbert-Schmidt Theorem (e.g. Simon, (1980), Thm. VI.16), it follows that there exists an orthonormal basis system $\left\{\boldsymbol{\psi}_{i},\ i\in\mathbb{N}\right\}$ of $\mathbb{H}^{L}$ such that

[TABLE]

Furthermore, using the Spectral Theorem (e.g. Werner, (2006), Thm. VI.3.2.) implies

[TABLE]

Since the kernel $\textbf{C}(s,u)$ is continuous, it admits the expansion

[TABLE]

This result is known as multivariate Mercer’s Theorem, (see e.g. Happ and Greven, (2018) Prop. 3). For any positive $i$, define an operator $\mathbfcal{X}_{i}:\mathbb{R}^{K}\rightarrow\mathbb{H}^{L}$, given by

[TABLE]

We call $\mathbfcal{X}_{i}$ an elementary operator. Note that $\mathbfcal{X}_{i}\in\mathbb{H}^{L\times K}$. Evaluating $\mathbfcal{X}_{i}\boldsymbol{a}$ at a given point $s\in[0,1]$ is equivalent to the matrix product ${\bf X}_{i}(s)\boldsymbol{a}$, where ${\bf X}_{i}(s)$ is an $L\times K$ matrix given by

[TABLE]

Note that, ${\bf X}_{i}(s)$’s can be considered as functional extension of the elementary matrices defined in (3), where ${\bf X}_{i}(s)$ is projecting columns of ${\bf X}(s)$ into a spaced spanned by $\boldsymbol{\psi_{i}}(s)$.

Proposition 3.3**.**

The elementary operators $\mathbfcal{X}_{i}$’s are bounded operators of rank one. Furthermore $\mathbfcal{X}_{i}$’s decompose the trajectory operator $\mathbfcal{X}$ as

[TABLE]

The next theorem provides the SVD of the trajectory operator $\mathbfcal{X}$, to obtain the associated eigentriples ($\sqrt{\lambda_{i}},\boldsymbol{v}_{i},\boldsymbol{\psi}_{i}$) in the decomposition step.

Theorem 3.1**.**

Let ${\{\boldsymbol{\psi}_{i}\}_{i=1}^{\infty}}$ and ${\{\lambda_{i}\}_{i=1}^{\infty}}$ be the eigenelements of $\bm{\mathcal{S}}$ given in (15), and

[TABLE]

Then, the SVD of the trajectory operator $\mathbfcal{X}$ can be written as

[TABLE]

where $\textbf{v}_{i}=\begin{pmatrix}\frac{\langle\boldsymbol{\psi}_{i},{\boldsymbol{x}}_{1}\rangle_{\mathbb{H}^{L}}}{\sqrt{\lambda_{i}}},\ldots,\frac{\langle\boldsymbol{\psi}_{i},{\boldsymbol{x}}_{K}\rangle_{\mathbb{H}^{L}}}{\sqrt{\lambda_{i}}}\end{pmatrix}^{\top}$. Furthermore for any $\boldsymbol{a}\in\mathbb{R}^{K}$, using (22) we have

[TABLE]

where

• i)

$\{\lambda_{i}\}_{i=1}^{\infty}$* is the set of non-ascending eigenvalues of $\bm{\mathcal{S}}^{\dagger},$ and*

• ii)

$\boldsymbol{v}_{i}$’s are the associated orthonormal eigenvectors of $\bm{\mathcal{S}}^{\dagger},$ satisfying $\mathbfcal{X}\boldsymbol{v}_{i}=\sqrt{\lambda_{i}}\boldsymbol{\psi}_{i}$.

We refer to $\sqrt{\lambda_{i}}$ as singular value, $\boldsymbol{\psi}_{i}$ as left singular function, and $\boldsymbol{v}_{i}$ as right singular vector associated to the $i^{th}$ component of the trajectory operator $\mathbfcal{X}$. The singular vectors, $\boldsymbol{v}_{i}$’s, can be used to produce paired plots similar to the ones in the SSA literature (Golyandina et al.,, 2001).

Step 3. Grouping

The grouping step is the procedure of rearranging and partitioning the elementary operators ${\mathbfcal X}_{i}$’s in (20). To do this, we mimic the approaches used in step 3 of Section 2 for the univariate SSA and implement the equivalent functional version of those in Haghbin et al., (2019). Note that, in practice, we use a finite set of elementary operators, and one can consider a partition $\{{I_{1}},{I_{2}},\ldots,{I_{m}}\}$ of the set of indices such that we have the expansion

[TABLE]

Step 4. Reconstruction

At this step, for any given $q$ ($1\leq q\leq m$), we would like to use $\mathcal{T}^{-1}:\mathbb{H}_{H}^{L\times K}\rightarrow\mathbb{H}^{N}$ to transform back each operator $\mathbfcal{X}_{I_{q}}$ in (24) to $\tilde{\bf y}_{N}^{q}$, and hence construct a functional version of the decomposition given in (5). But since $\mathbfcal{X}_{I_{q}}\in\mathbb{H}^{L\times K}$, first it is necessary to project $\mathbfcal{X}_{I_{q}}$ to $\mathbb{H}_{H}^{L\times K}$. Note that $\mathbb{H}_{H}^{L\times K}$ is a closed subspace of $\mathbb{H}^{L\times K}$, therefore by Projection Theorem, there exist a unique $\mathbfcal{\tilde{X}}_{I_{q}}\in\mathbb{H}_{H}^{L\times K}$ such that

[TABLE]

To specify $\mathbfcal{\tilde{X}}_{I_{q}}$, we denote the elements of $\mathbfcal{X}_{I_{q}}$ and $\mathbfcal{\tilde{X}}_{I_{q}}$ by $[x_{i,j}^{q}]$ and $[\tilde{x}_{i,j}^{q}]$, respectively. Using Lemma 3.1, it is easy to extend the diagonal averaging approach given by Golyandina et al., (2001) to $\mathbb{H}^{L\times K}$ and obtain $\tilde{x}_{i,j}^{q}$’s as following:

[TABLE]

where $s=i+j$ and $n_{s}$ stands for the number of $(l,k)$ pairs such that $l+k=s$. Denote this projection by $\Pi_{\mathbb{H}}:\mathbb{H}^{L\times K}\rightarrow\mathbb{H}^{L\times K}_{H}$, and set $\mathbfcal{\tilde{X}}_{I_{q}}=\Pi_{\mathbb{H}}\mathbfcal{X}_{I_{q}}$. Now we can define $\tilde{\bf y}_{N}^{q}=\mathcal{T}^{-1}\tilde{\mathbfcal{X}}_{I_{q}}$, and reconstruct the functional time series.

3.2 Separability

Let $\textbf{y}_{N}=\textbf{y}_{N}^{(1)}+\textbf{y}_{N}^{(2)}$, where $\textbf{y}_{N}^{(i)}=\{y_{1}^{(i)},\ldots,y_{N}^{(i)}\}$, $i=1,2$, are FTS. Using a fixed window length $L$, for each series $\textbf{y}_{N}^{(i)}$, denote $\{{\boldsymbol{x}}_{k}^{(i)}\}_{k=1}^{K}$ as a sequence of functional lagged vectors, and $\mathcal{L}^{(i)}$ as the linear space spanned by $\{{\boldsymbol{x}}_{k}^{(i)}\}_{k=1}^{K}$. Analogous to univariate SSA, separability of the series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ is equivalent to $\mathcal{L}^{(1)}\bot\ \mathcal{L}^{(2)}$, which is same as $\langle{\boldsymbol{x}}_{k}^{(1)},{\boldsymbol{x}}_{k^{\prime}}^{(2)}\rangle_{\mathbb{H}^{L}}=0$, for all $k,k^{\prime}=1,\ldots,K$. Furthermore, a necessary condition for separability can be defined based on w-correlation measure. To do this, consider the weighted inner product of two series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ as

[TABLE]

where $w_{i}=\min\{i,L,N-i+1\}$. We call the series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ w-orthogonal if

[TABLE]

Theorem 3.2**.**

If the series $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ are separable, then they are w-orthogonal.

Also, to quantify the degrees of separability of two FTS, the functional version of the w-correlation measure can be obtained by replacing the new definition of the weighted inner product (26) into (6).

4 Implementation Strategy

In practice, functional data are being recorded discretely and then converted to functional objects using proper smoothing techniques. We refer to Ramsay and Silverman, (2007) for more details on preprocessing the raw data. Let $\left\{\nu_{i}\right\}_{i\in\mathbb{N}}$ be a known basis system (not necessarily orthogonal) of the space $\mathbb{H}$. Each functional observation in $\mathbb{H}$ can be projected into subspace ${\mathbb{H}}_{d}:=sp\left\{\nu_{i}\right\}$, where $d$ can be determined by variety of techniques (e.g. cross-validation). Therefore, each $y_{j}\in\mathbb{H}_{d}$ is uniquely represented by

[TABLE]

Let us define quotient sequence, $q_{k}$, and reminder sequence, $r_{k}$, by

[TABLE]

Note that for any given $k$ ($1\leq k\leq Ld$), one may use (29) to determine $q_{k}$ and $r_{k}$ uniquely, so these sequences are well defined. Now, consider the objects ${\boldsymbol{\phi}}_{k}\in\mathbb{H}_{d}^{L}$, as a vector of length $L$ with all coordinates are zero except $r_{k}$-th, which is $\nu_{q_{k}}$.

Lemma 4.1**.**

The sequence $\{{\boldsymbol{\phi}}_{k}\}_{k=1}^{Ld}$ is a basis system for $\mathbb{H}_{d}^{L}$, where $\mathbb{H}_{d}^{L}$ is the Cartesian product of $L$ copies of $\mathbb{H}_{d}$.

Using Lemma 4.1, each element ${\boldsymbol{z}}\in\mathbb{H}_{d}^{L}$ admits a unique representation

[TABLE]

where ${\mathbf{c}_{\boldsymbol{z}}}=\left(\langle{\boldsymbol{z}},\widetilde{{\boldsymbol{\phi}}}_{1}\rangle_{\mathbb{H}^{L}},\ldots,\langle{\boldsymbol{z}},\widetilde{{\boldsymbol{\phi}}}_{Ld}\rangle_{\mathbb{H}^{L}}\right)^{\top}$ corresponds to ${\boldsymbol{z}}$, ${\mathbfcal{F}}:\mathbb{R}^{Ld}\rightarrow\mathbb{H}_{d}^{L}$ belongs to $\mathbb{H}^{Ld\times L}_{d}$, and $\widetilde{{\boldsymbol{\phi}}}_{i}$ is the dual basis of ${\boldsymbol{\phi}}_{i}$. Note that, in the special case, when $\nu_{i}$’s are orthonormal (so ${\boldsymbol{\phi}}_{i}$’s are), $\widetilde{{\boldsymbol{\phi}}}_{i}={\boldsymbol{\phi}}_{i}$ (see Christensen,, 1995, for more details). Applying the linear operator $\mathbfcal{S}$, defined in (12), on (30) implies

[TABLE]

where ${\mathbf{S}}^{\top}=\left[\langle{\mathbfcal{S}}{\boldsymbol{\phi}}_{i},\widetilde{{\boldsymbol{\phi}}}_{j}\rangle_{\mathbb{H}^{L}}\right]_{i,j=1}^{Ld}$. We call ${\mathbf{S}}$ the corresponding matrix of ${\mathbfcal{S}}$.

Lemma 4.2**.**

Let ${\boldsymbol{z}}$ be a functional object in $\mathbb{H}_{d}^{L}$, then ${\mathbfcal{S}}\boldsymbol{z}=\lambda{\boldsymbol{z}}$ if and only if ${\mathbf{S}}{\mathbf{c}_{\boldsymbol{z}}}=\lambda{\mathbf{c}_{\boldsymbol{z}}}$.

Theorem 4.1**.**

Suppose the Gram matrix $\mathbf{G}:=\left[\langle{\boldsymbol{\phi}}_{i},{\boldsymbol{\phi}}_{j}\rangle_{\mathbb{H}^{L}}\right]_{i,j=1}^{Ld}$. Then the following holds:

• (i)

$\mathbf{G}=\left[\delta_{r_{i},r_{j}}\langle\nu_{q_{i}},\nu_{q_{j}}\rangle\right]_{i,j=1}^{Ld}$.

• (ii)

$\left\langle\mathbfcal{S}{\boldsymbol{\phi}}_{i},{\boldsymbol{\phi}}_{j}\right\rangle_{\mathbb{H}^{L}}=\sum_{m=1}^{K}\langle y_{r_{i}+m-1},\nu_{q_{i}}\rangle\langle y_{r_{j}+m-1},\nu_{q_{j}}\rangle$.

• (iii)

${\mathbf{S}}_{0}:=\left[\langle{\mathbfcal{S}}{\boldsymbol{\phi}}_{i},{\boldsymbol{\phi}}_{j}\rangle_{\mathbb{H}^{L}}\right]_{i,j=1}^{Ld}$* is a symmetric matrix.*

• (iv)

${\mathbf{S}}=\mathbf{G}^{-1}\mathbf{S}_{0}$.

Now we have the recipes to proceed with the following algorithm and obtain the eigenfunctions of $\mathbfcal{S}$, $\boldsymbol{\psi}_{i}$’s, used in the decomposition step. For a given set of basis $\{\nu_{i}\}_{i=1}^{d}$, and a FTS, $\textbf{y}_{N}$:

• •

Use Theorem 4.1 to compute the matrices $\mathbf{G},\ \mathbf{S}_{0}$, and $\mathbf{S}$.

• •

Use the eigendecomposition of $\mathbf{S}$ to obtain eigenpairs $(\lambda_{i},{\mathbf{c}_{{\boldsymbol{\psi}}_{i}}})$ for $i=1,\ldots,Ld$.

• •

Use (30) and Lemma 4.2 to obtain $\boldsymbol{\psi_{i}}$’s, eigenfunctions of $\mathbfcal{S}$.

Now, one can use ${{\boldsymbol{\psi}}_{i}}$’s to decompose the FTS to elementary operators $\mathbfcal{X}_{i}$’s. Note that, in practice, we represent the elementary operators in matrix form. Therefore, one may observe that the equivalent functional elementary matrix ${\bf X}_{i}(s)$, given in (3.1), is just a projection of the functional lagged vectors $\{{\boldsymbol{x}}_{j}(s)\}_{j=1}^{K}$, given in (8), onto ${{\boldsymbol{\psi}}_{i}}$. Furthermore, in the diagonal averaging step, we can incorporate the averaging over the associate basis coefficients of $x_{l,k}^{q}$ in (25) to obtain the respective basis coefficients for $\tilde{x}_{i,j}^{q}$. For more details see Haghbin et al., (2019).

5 Numerical study

In this section, first, we present a simulation study to elaborate the use of the FSSA compared with dFPCA and MSSA under different scenarios. To do so, we utilize the implementation of the proposed model that is available as an R package named Rfssa in the CRAN repository (Haghbin et al.,, 2019). We also use the freqdom.fda (Hörmann et al.,, 2015) and Rssa (Golyandina et al.,, 2015) packages to obtain the dFPCA and MSSA results. In the second subsection, we analyze a remote sensing data using Rfssa and provide some visualization tools that come handy in the grouping step.

We developed a shiny app, included in the Rfssa package, also available at https://fssa.shinyapps.io/fssa/, to demonstrate and reproduce different aspects of the simulation setup. Furthermore, it can be used to compare the results of dFPCA, MSSA, and FSSA on the remote sensing data, call center dataset or any other FTS, provided by the end-user.

5.1 Simulation study

For the simulation setup, consider the functional time series of lengths $N=50,100,150$ and $200$ which are observed in $n=100$ fixed equidistant discrete points on $[0,1]$ from the following model:

[TABLE]

A cubic B-spline basis functions with 15 degrees of freedom is used to convert $\{Y_{t}(s_{i})\}$’s into smooth (continuous) functional curves. In this model, $m_{t}(s)$ is considered to be a periodic component defined as

[TABLE]

where $\omega$ is the model frequency with three different values ($\omega=0,0.1$ and $0.25$).

The $X_{t}(s)$ in (32) is a stochastic term that is generated under four different settings with an increasing trend in complexity. In the first setting, we consider $\{X_{t}(s_{i}),t=1,\ldots,N\ \text{and}\ i=1,\ldots,n\}$ are drawn from an independent Gaussian White Noise (GWN) process with zero mean and standard deviation equal to $0.1$. It is expected to obtain an acceptable performance from FPCA for reconstructing the FTS in the first setting as intuitively FPCA outperforms under this ideal framework (see Maadooliat et al.,, 2015, for more details).

In the remaining three settings, the $\{X_{t}(s)\}$ processes are simulated from a functional autoregressive model of order 1, $FAR(1)$, defined by

[TABLE]

where $\Psi$ is an integral operator with a parabolic kernel as follow

[TABLE]

The constant $\gamma_{0}$ is chosen such that the Hilbert-Schmidt norm defined by

[TABLE]

acquires the values $\|\Psi\|_{\mathcal{S}}^{2}=0,0.5$, and $0.9$, for the remaining three settings, respectively. In these settings, the white noise terms $\varepsilon_{t}(s)$ are considered as independent trajectories of the standard Brownian motion over the interval $[0,1]$. It is worth to note that as we increase the Hilbert-Schmidt norm, $\|\Psi\|_{\mathcal{S}}^{2}$, in the FAR(1) models, the dependency structure of consecutive FTS gets more twisted, and we expect it would be more challenging to reconstruct the true structures, $\{X_{t}(s)\}$.

To compare the performance of FSSA and MSSA, we further consider three window length parameters ($L=20,30$ and $40$) in our simulation setup. For the sake of consistency in all of the reconstruction procedures (dFPCA, MSSA and FSSA), we use the first two leading eigen-components. As a measure of goodness of fit, we use the Root Mean Square Error (RMSE) defined as:

[TABLE]

where $\hat{Y}_{t}(s_{i})$ is the FTS reconstructed by each method. We repeat each setting $1000$ times and report the mean of the RMSE’s in Table 1.

By comparing the results in Table 1, it can be seen that FSSA outperform dFPCA in different scenarios. This may not be surprising, as the main task of dFPCA is dimension reduction. Except for the first setting, MSSA also outperforms dFPCA significantly. Furthermore, FSSA performs better than MSSA in most of the cases except the case where the length of the FTS is small ($N=50$) and the window size, $L$, is getting closer to $N$. However, it is clear that FSSA is the optimal method for reconstructing the longer FTS ($N\geq 100$).

For all methods, RMSE decreases as the length of the series increases. For two smaller frequencies ($\omega=0\ \textrm{and}\ 0.1$), the average of RMSE increases as the noise structure becomes more complex in settings 1 through 4 while it decreases for $\omega=0.25$. This might be happening due to the unpredicted cross-correlation of the functional noise structures and the periodic form of FTS.

The efficiency of MSSA and FSSA for different window lengths ($L$), time series lengths ($N$) and frequencies ($\omega$), the ratio of RMSE of MSSA to FSSA is examined in Figure 3. The overall pattern confirms the improvement in RMSE for FSSA as the time series get longer (larger $N$). Overall, as $L$ is increasing, the pattern of ratio of RMSE’s remains unchanged. Although as the window length becomes larger, either the improvement diminishes for longer FTS, or disappears (or reverses) for smaller $N$. It is also worth to note that in setting 1 (GWN), based on the right panel of Figure 3 and Table 1, the FSSA dominates the other two methods in all combinations of parameters with a better efficiency scale.

5.2 Application to Remote Sensing Data

Tracking changes in vegetation over time has become of interest to researchers that want to preserve wildlife. One technique used to track how much vegetation is present in a region is through field surveys. The problem with this technique is that it is difficult to implement especially in low population areas (Panuju and Trisasongko,, 2012). The use of remote sensing data in this context is preferred in order to detect man-made or natural changes in vegetation. One source of remote sensing data is from NASA’s MODerate-resolution Imaging Spectroradiometer (MODIS) satellite which provides images, twice daily, of regions around the globe at varying spatial resolutions (Tuck et al.,, 2014). Normalized Difference Vegetation Index (NDVI) is a commonly used pixel-wise index in MODIS satellite images. The NDVI values are bounded between zero and one, where index values that are closer to one indicate that more vegetation is present and smaller values indicate the absence of vegetation (Tuck et al.,, 2014).

Many studies have used the spectral NDVI measure and its variants in order to remotely track changes in vegetation over time (Lambin,, 1999). The temporal average, and temporal variability of NDVI images have been used as explanatory variables for the number of different types of vegetation present in many regions (Tuck et al.,, 2014). Also, Panuju and Trisasongko, (2012) used the maximum value of NDVI images, taken of the Jambi province in Indonesia, to form a time series. The resulted time series was then analyzed using an X12-ARIMA model in order to identify trend and seasonal changes in woody vegetation (Panuju and Trisasongko,, 2012). While these statistics (e.g. maximum and average) have seen some success in tracking the changes of the NDVI images, it may fail capturing the distribution of the vegetation. Therefore, one may seek for a more comprehensive measure (e.g. a FTS) to describe the distribution of the NDVI images.

Here we use $448$ NDVI images taken in $16$ days increments between February $18$, $2002$ to July $28$, $2019$. The images have a spatial resolution of $250$ meter and are of a square region just outside of the city of Jambi, Indonesia between $103.61^{\circ}E-103.68^{\circ}E$ and $1.60^{\circ}S-1.67^{\circ}S$. Figure 4 shows the respective NDVI images taken on June 10, 2002 and June 10, 2019. It is interesting to note that, although the respective NDVI images are not similar, the sample means of the NDVI values are very close (differ by only about 0.0032 unit). As it is shown in Figure 4, the kernel density estimates (KDEs) are much more informative in representing the distribution of the NDVI images.

We follow the Silverman’s rule of thumb (Silverman,, 1986) for bandwidth selection and obtain the KDEs. Then we project the results onto a functional space spanned by a cubic B-spline basis, selected via the GCV criterion. Figure 5(A) shows the projected KDEs onto the B-spline space. We pass the results ($448$ FTS) as input to the FSSA algorithm with the lag, $L=45$, and study the behavior of the NDVI images over almost two decade, using the Rfssa package.

According to the subplots in Figures 5(B-C: the singular values and w-correlation plots), a suitable partition would be grouping the first and fourth components separately, plus the second and third components jointly ($G=\{1,2$–$3,4\}$). The remaining subplots in Figures 5(D-F: the right singular vectors and left singular functions) indicate the first component captures mean behavior, the second and third capture annual behavior, and the fourth captures trend.

An important goal of analyzing the NDVI data is to investigate the existence of a temporal trend in the NDVI images over time. It is interesting to see that FSSA can distinguish and separate the overall mean structure (top-left subfigures in Figures 5(D-F)) and the trend pattern (bottom-right subfigures in Figures 5(D-F)) in two different components (the first and fourth components). The trend component is causing an interesting change-point behavior after almost a decade (bottom-right subfigures in Figures 5(D-F) and Figure 6(C)). We would like to report that SSA algorithm with the lag, $L=45$, cannot separate this trend structure, and the overall mean pattern. Therefore, SSA combines these two (the overall mean and the trend components) into one component (the SSA results are omitted due to space constraints).

Next, we build the reconstruction of the NDVI images using the grouping suggested by the FSSA exploratory plots, $G$, as shown in Figure 6. It is clear that Figure 6(A) shows the reconstructed overall mean, where it does not change over time. Figure 6(B) provides the annual harmonic behavior, and one may observe the change-point behavior after a decade in Figure 6(C). The last subfigure, Figure 6(D), presents the sum of the trend and overall mean components.

It would be of interest to confirm the properties of the reconstructed groups using some rigorous statistical procedures. In order to do this, we provide the multivariate trace periodicity test of Hörmann et al., (2018) to test for the annual seasonal effect of FTS, and a bootstrap procedure to test existence of a trend (Grosjean and Ibanez,, 2018) over the time series of the mean of the coefficients associated to the B-spline basis. We obtain the results of these two tests (periodicity and trend) on four sets of FTS (original signal $y_{t}(s)$, $R_{1}$, $R_{2}$, $R_{3}$), where $R_{i}$ represents residual curves obtained via removing the reconstructed FTS by the first $i$ groups, from the original signal $y_{t}(s)$. Table 2 provides the p-values of the tests for the periodicity and the trend. It is clear that periodicity test captures the annual pattern for $y_{t}(s)$ and $R_{1}$ that contain the seasonal components (p-value=0). Also the trend test is significant for all of the FTS, except $R_{3}$ that does not contain the third group (fourth component).

While the FSSA procedure was able to separate out the year long seasonal pattern, it was also able to detect the less obvious trend component present in the NDVI data; which indicates a loss of vegetation over the previous 20 years. It was determined that between 2001 and 2015, grassland and shrubs in the Jambi region accounted for a large amount of the vegetation lost due to controlled burning for human use (Prasetyo et al.,, 2016). The region specified in images that we study here is primarily upland agriculture and grass and it appears to have been a hot spot for controlled fires (Prasetyo et al.,, 2016).

6 Discussion

In this paper, we constructed the FSSA procedure by incorporating the FDA techniques in basic SSA via MFPCA. The contribution of the proposed model is to provide practitioners with some tools to utilize the advantages of SSA in FTS. Accordingly, the researchers can analyze functional sequences (e.g., time series, longitudinal, or spatial data) via FSSA. Alternatively one may approach the problem using MSSA, given that the data points are measured in fixed, regular grid points over time.

As for the ease of use, an efficient and user-friendly R implementation of FSSA is developed in the Rfssa package. Furthermore, a shiny web application is also included in the package, and it is available at https://fssa.shinyapps.io/fssa/ for reproducing the results of this paper or analyzing any other FTS.

Supplementary Materials

Application to call center dataset

To illustrate the advantages of FSSA, especially its main capability in extracting different functional components (i.e. trend, harmonic and noise), we explore the call center dataset analyzed in Maadooliat et al., (2015). This dataset provides the number of calls to a call center per 6 minutes intervals, between January 1 through December 31, 1999. Suppose $y_{t}(s_{i}),\ t=1,\ldots 365,\ i=1,\ldots,240$, is the square root of number of calls during the time interval $s_{i}$ on day $t$. Figure 1 (left) shows the projection of the $\{y_{t}(s_{i})\}$’s (vectors of length $240$) into a functional space spanned by a cubic B-spline using GCV criterion.

An important goal of analyzing the call center data is to investigate the existence of periodic behaviors (e.g., weekly or monthly). Figure 1 (right) visually confirms the existence of a strong weekly pattern in the dataset. Since one cannot visually confirm the presence of a monthly behavior using similar graphs, it would be interesting to show that FSSA can provide tools and machinery to extract such weaker signals.

In order to capture both monthly and weekly pattern by FSSA, first, we choose window length as multiple of $7$ and close to $30$, i.e., $L=28$. Then, we provide several plots using Rfssa package for grouping the components (Figure 7). These plots are the functional form (analogy) of the ones commonly used in the SSA literature (Golyandina et al.,, 2001). As it can be seen in the plot of leading singular values (scree plot), the first singular value is relatively large, and there exists three evident pairs with almost equal leading singular values correspond to the three components. The w-correlation plot suggests partitioning the eigentriples into five groups: $1$, $2-3$, $4-5$, $6-7$, $8-9$ and the remainder that does not seem to contain any strong signal. Considering the remaining plots (right singular vectors and pairs of singular vectors), one can see the eigentriple pairs $2-3$, $4-5$ and $6-7$ are related to a one-week periodicity with frequencies $1/7$, $2/7$ and $3/7$, while the last group, eigentriple pair $8-9$, describes a weak monthly cycle. These groups can reproduce the reconstructed FTS. The functional components associated with the first four groups are presented in Figure 2 (bottom) from left to right. Furthermore, some creative visualization tools that are implemented in the Rfssa package can be used to extract within/between days patterns for the call center data by employing the estimated left singular functions (Figure 8). It is worth to mention that in Figure 8 (right), there are 28 curves associated with all nine singular functions (graphs). One may note that for the first singular function, all curves resemble a similar pattern respective to the main trend. Furthermore, singular functions $2-7$ contain seven distinguish patterns that each consists of four curves whereas $28$ distinct curves construct singular functions $8-9$.

For further clarification, we provide the multivariate trace periodicity test of Hörmann et al., (2018) on six sets of FTS (original signal $y_{t}(s)$, $R_{1}$, $R_{2}$, $R_{3}$, $R_{4}$, $R_{5}$), where $R_{i}$ represents residual curves obtained via removing the reconstructed FTS by the first $i$ groups, from the original signal $y_{t}(s)$. Table 3 provides the p-values of the test for the periods of length $7$ and $30$ days (p-values for testing the weekly and monthly patterns). It is clear that periodicity test captures the weekly pattern for $y_{t}(s)$, $R_{1}$, $R_{2}$ and $R_{3}$ that contain either all or part of the weekly components (p-value=0). After subtracting the functional mean of all curves (first eigentriple) and the weekly components (eigentriple 2-7), the monthly pattern in $R_{4}$ is not weak anymore, and the periodicity test can capture the monthly cycle in $R_{4}$ (p-value=0). Finally, $R_{5}$ is the remainder of the signal after removing all of the weekly and monthly components, and that’s why the associated p-values in the last row are not significant anymore.

Proof of theorems and propositions

Proof of Lemma.

3.1

[TABLE]

∎

Proof of Prop.

3.1 For given ${\boldsymbol{a}}\in\mathbb{R}^{K}$ and ${\boldsymbol{z}}\in\mathbb{H}^{L}$ we have

[TABLE]

∎

Proof of Prop.

3.2 The linearity follows by bilinear form of the inner product. To prove self-adjoint property, for given ${\boldsymbol{z},w}\in\mathbb{H}^{L}$, (12) gives

[TABLE]

which implies $\bm{\mathcal{S}}$ is self-adjoint. Moreover, we have

[TABLE]

which implies that $\bm{\mathcal{S}}$ is positive definite. To prove boundedness, let $C=\sum_{j=1}^{K}\|{\boldsymbol{x}}_{j}\|_{\mathbb{H}^{L}}^{2}\geq 0$. Then

[TABLE]

As an immediate result of boundedness one can show the continuity. If ${\boldsymbol{z}}\in\mathbb{H}^{L}$ and ${\boldsymbol{z}}_{n}$ is a sequence in $\mathbb{H}^{L}$ such that $\|{\boldsymbol{z}}_{n}-{\boldsymbol{z}}\|_{\mathbb{H}^{L}}\rightarrow 0$ as $n\rightarrow\infty$, then we have

[TABLE]

In other words, $\bm{\mathcal{S}}{\boldsymbol{z}}_{n}\rightarrow\bm{\mathcal{S}}{\boldsymbol{z}}$, which proves the continuity of $\bm{\mathcal{S}}$. To prove compactness, define $B:=\max_{i=1,\ldots,N}\|{y}_{i}\|^{2}$, and let $\{{\boldsymbol{f}}_{n}\}$ be a weak-null sequence in $\mathbb{H}^{L}$, i.e.

[TABLE]

for all $\boldsymbol{a}\in{\mathbb{H}^{L}}$. Then using (13) and (14), with some efforts, we have

[TABLE]

Therefore, $\bm{\mathcal{S}}{\boldsymbol{f}}_{n}\rightarrow 0$ as $n\rightarrow 0$, which implies the compactness (e.g. Weidmann, (1980) Thm 6.3). ∎

Proof of Prop.

3.3 For any $i\in\mathbb{N}$, define \textbf{s}_{i}:=\bigl{(}\langle\boldsymbol{\psi}_{i},\boldsymbol{x}_{1}\rangle_{\mathbb{H}^{L}},\langle\boldsymbol{\psi}_{i},\boldsymbol{x}_{2}\rangle_{\mathbb{H}^{L}},\allowbreak\ldots,\langle\boldsymbol{\psi}_{i},\boldsymbol{x}_{K}\rangle_{\mathbb{H}^{L}}\bigr{)}^{\top}. Then using the definition of the operator $\mathbfcal{X}_{i}$, for a given ${\boldsymbol{a}}\in\mathbb{R}^{K}$, we have

[TABLE]

which yields that $\mathbfcal{X}_{i}$ is a bounded operator of rank one (e.g. Weidmann, (1980) Thm. 6.1). Moreover, since $\left\{\boldsymbol{\psi}_{i},\ i\in\mathbb{N}\right\}$ is an orthonormal basis system for $\mathbb{H}^{L}$, for all $\boldsymbol{a}\in\mathbb{R}^{K}$ we have

[TABLE]

∎

Proof of Thm.

3.1 Note that $\bm{\mathcal{S}}^{\dagger}:\mathbb{R}^{K}\rightarrow\mathbb{R}^{K}$ while $\bm{\mathcal{S}}:\mathbb{H}^{L}\rightarrow\mathbb{H}^{L}.$ Furthermore,

[TABLE]

To prove orthonormality, we have

[TABLE]

Moreover, using (9) we have

[TABLE]

∎

Proof of Thm.

3.2 The Equation (26) gives

[TABLE]

Hence, separability of $\textbf{y}_{N}^{(1)}$ and $\textbf{y}_{N}^{(2)}$ implies $\langle{\boldsymbol{x}}_{k}^{(1)},{\boldsymbol{x}}_{k}^{(2)}\rangle_{\mathbb{H}^{L}}=0$ for all $k=1,\ldots,K$, which complete the proof. ∎

Proof of Lemma.

4.1 The proof will be divided into two steps. In the first step it will be shown that $\mathbb{H}_{d}^{L}=sp\{{\boldsymbol{\phi}}_{1},\ldots,{\boldsymbol{\phi}}_{Ld}\}$. In the second step it will be proved that ${\boldsymbol{\phi}}_{1},\ldots,{\boldsymbol{\phi}}_{Ld}$ are linearly independent. Suppose that ${\boldsymbol{z}}=\left(z_{1},\ldots,z_{L}\right)^{\top}$, where $z_{i}\in\mathbb{H}_{d}$. By definition, each elements of ${\boldsymbol{z}}$ admits the basis expansions $z_{j}=\sum_{i=1}^{d}b_{i,j}\nu_{i},\ j=1,\ldots,L$. Therefore

[TABLE]

which implies the first step. To prove linear independency, if $\sum_{k=1}^{Ld}a_{k}{\boldsymbol{\phi}}_{k}={\boldsymbol{0}}$ then,

[TABLE]

[TABLE]

This means $\sum_{j=1}^{d}a_{(j-1)L+i}\nu_{j}=0$ for all $i=1,\ldots,L$ and consequently $a_{(j-1)L+i}=0$ for $j=1,\ldots,d$, since $\{\nu_{i}\}_{i=1}^{d}$ are linear independent. ∎

Proof of Thm.

4.1 The proof (i) is clear, and (iii) is straightforward from (ii). To prove the part (ii), using (13) gives

[TABLE]

Finally, to prove (iv), the Equation (30) gives,

[TABLE]

The Gram matrix $\mathbf{G}$ is Hermitian and nonsingular, since the basis ${\boldsymbol{\phi}}_{i}$ are linearly independent. Let $\mathbf{G}^{-1}:=\left[h_{i,j}\right]_{i,j=1}^{Ld}$ be the inverse of $\mathbf{G}$. The Equation (36) gives

[TABLE]

which implies

[TABLE]

Applying this yields

[TABLE]

which implies $\mathbf{S}=\mathbf{G}^{-1}\mathbf{S}_{0}$.

∎

R-package for FSSA routine

Rfssa package includes the code and a shiny app, https://fssa.shinyapps.io/fssa/, to perform the FSSA algorithm, and present the illustrative figures described in this work. The package also contains the NDVI remote sensing and the call center datasets used as examples in the article. (available in CRAN)