The decomposition and forecasting of mutual investment funds using singular spectrum analysis

entropy Article The Decomposition and Forecasting of Mutual Investment Funds Using Singular Spectrum Analysis Paulo Canas Rodrigues 1,2,* , Jonatha Pimentel 1and Patrick Messala 1 and Mohammad Kazemi 3 1Department of Statistics, Federal University of Bahia, 40170-110 Salvador, Brazil; [email protected] (J.P.); [email protected] (P.M.) 2CAST, Faculty of Information Technology and Communication Sciences, Tampere University, FI-33014 Tampere, Finland 3Department of Statistics, Faculty of Mathematical Sciences, Shahrood University of Technology, P.O. Box 3619995161 Shahroud, Iran; [email protected] *Correspondence: [email protected] or [email protected] Received: 7 November 2019; Accepted: 7 January 2020; Published: 9 January 2020   Abstract: Singular spectrum analysis (SSA) is a non-parametric method that breaks down a time series into a set of components that can be interpreted and grouped as trend, periodicity, and noise, emphasizing the separability of the underlying components and separate periodicities that occur at different time scales. The original time series can be recovered by summing all components. However, only the components associated to the signal should be considered for the reconstruction of the noise-free time series and to conduct forecasts. When the time series data has the presence of outliers, SSA and other classic parametric and non-parametric methods might result in misleading conclusions and robust methodologies should be used. In this paper we consider the use of two robust SSA algorithms for model fit and one for model forecasting. The classic SSA model, the robust SSA alternatives, and the autoregressive integrated moving average (ARIMA) model are compared in terms of computational time and accuracy for model fit and model forecast, using a simulation example and time series data from the quotas and returns of six mutual investment funds. When outliers are present in the data, the simulation study shows that the robust SSA algorithms outperform the classical ARIMA and SSA models. Keywords: singular spectrum analysis; robust singular spectrum analysis; time series forecasting; mutual investment funds 1. Introduction Mutual investment funds provide management services to institutional and individual investors, besides great liquidity for financial investments made in them and low transactional costs [ 1 , 2 ]. These funds can be of fixed or variable income and allow to diversify the assets while reducing unsystematic risk. Fixed income mutual investment funds are of low risk, whereas variable-income mutual investment funds vary in terms of risk but also in terms of returns. In this study, we were interested in analyzing the quotas and returns of six of the largest Brazilian based mutual investment funds—three purely based on stocks: (i) Alaska Black, (ii) APEX Long Biased, and (iii) Brasil Capital; and three balanced funds (usually combining a stock component, a bond component, and sometimes a money market component in a single portfolio): (iv) ADAM Strategy, (v) Gavea Macro, and (vi) SPX Nimitz. A natural framework for analyzing mutual investment funds, due to its underlying structure, is a time series method. Entropy 2020,22, 83; doi:10.3390/e22010083 www.mdpi.com/journal/entropy Entropy 2020,22, 83 2 of 22 Singular spectrum analysis (SSA) is a powerful non-parametric technique for time series analysis and forecasting, which incorporates elements of classical time series analysis, multivariate statistics, and matrix algebra. Its main aim is to decompose the original time series into a set of components that can be interpreted as trend components, seasonal components, and noise components [ 3 – 6 ]. SSA has proven both wide usefulness and applicability across many applications [ 7 – 17 ], being that its scope of application ranges from parameter estimation to time series filtering, synchronization analysis, and forecasting [18]. The SSA methodology for model fit can be summarized in four steps: (i) embedding, which maps the original univariate time series into a trajectory matrix; (ii) singular value decomposition (SVD), which helps decomposing the trajectory matrix into the sum of rank-one matrices; (iii) eigentriple grouping, which helps deciding which of the components are associated to the signal and which are associated to the noise; and (iv) diagonal averaging, which maps the rank-one matrices, associated to the signal, back to time series that can be interpreted as trend, seasonal, or other meaningful components. SSA results and interpretation, similarly to many other classical time series methods, can be sensitive to data contamination with outliers [ 19 , 20 ]. In those cases, even a small percentage of outliers can make a big difference on the results for model fit and model forecast. Very few attempts have been made in order to access the effect of the presence of outliers in the data while conducting a SSA. One study [ 21 , 22 ] presented some preliminary results on the effect of outliers in singular spectrum analysis, and [23] made a first attempt to robustify the SSA by considering an SVD based on a robust L1norm [24] instead of the L2norm used in the classical algorithm, which they used for model fit. In this paper we go one step further than [ 23 ] and propose a new robust algorithm for SSA that considers the SVD based on the Huber function [ 25 ]. Moreover, we propose two robust SSA forecasting algorithms, one based on the the L1 norm and another based on the Huber function. Comparisons are made between the classical SSA algorithm, the robust SSA algorithm based on the L1 norm (RLSSA), the robust SSA algorithm based on the Huber function (RHSSA), and the classical autoregressive integrated moving average (ARIMA) model, in terms of computational time and accuracy for model fit and model forecast. These comparisons for decomposing and forecasting time series were done by considering a simulation example and the six mutual investment funds mentioned above. The rest of this paper is organized as follows. Section 2provides the materials and methods containing the data description, a brief introduction to the ARIMA and SSA methodologies, and the details of the proposed robust SSA algorithm that uses the SVD based on the Huber function. Section 3 presents the results and discussion, wherein the ARIMA, SSA, and robust SSA algorithms are compared in terms of model fit and model forecast, using the six mutual investment funds and the simulation example. The paper closes in Section 4, wherein some conclusions are drawn. 2. Materials and Methods 2.1. Data In this paper we consider a dataset that includes daily observations of six mutual investment funds, three based purely on stocks and three balanced funds: Stock funds •Alaska Black: 3 January 2017–30 August 2019 (N= 666 observations). •APEX Long Biased: 15 April 2013–30 August 2019 (N= 1604 observations). •Brasil Capital: 27 August 2012–30 August 2019 (N= 1760 observations). Balanced funds •ADAM Strategy: 29 April 2016–30 August 2019 (N= 838 observations). •Gavea Macro: 30 June 2008–30 August 2019 (N= 2809 observations). •SPX Nimitz: 01 December 2010–30 August 2019 (N= 2199 observations). The datasets were collected from https://infofundos.com.br/carteira. Entropy 2020,22, 83 3 of 22 2.2. ARIMA Model The autoregressive integrated moving average (ARIMA) models are among the most widely used techniques for time series analysis and forecasting. Such a model depends on three parameters: p is the number of lagged observations in the model, i.e., the autoregressive (AR) order; d is the number of times that the original observations are differenced, i.e., the integrated (I) degree; and q is the size of the moving average window, i.e., the order of the moving average (MA) [ 26 ]. This parametric model can then be written as ARIMA(p , d , q) , with p , d , and q non-negative integers. Given a time series YN=y1, . . . , yN, the ARIMA(p,d,q)model can be written as: (1−φ1B1−···−φpBp)(1−B)dyt=c+ (1+θ1B1+···+θqBq)εt, (1) where φ1 , . . . , φp are the parameters or coefficients of the p autoregressive terms; B is the time lag operator, or backward shift, which is a linear operator denoted by Bk such that Lkyt=yt−k , t∈Z ; yt is the observation at the time point t ; c=µ( 1 −φ1−···−φp) ; µ is the mean of ( 1 −B)dyt ; β1 , . . . , βq are the parameters or coefficients of the q moving average terms; and εt is an error term, usually white noise with variance σ2. Alternatively, the model can be written as: (1−φ1B−···−φpBp)(1−B)d(yt−µtd/d!) = (1+θ1B+···+θqBq)εt, (2) which is the parametization used in the “arima” function of the software R [27]. 2.3. Singular Spectrum Analysis Singular spectrum analysis is a non-parametric technique for model fit and model forecasting that decomposes a time series into a number of components that are summed and interpreted as trend, periodicity, and noise. Similarly to many other time series techniques, SSA can be used for solving a wide range of problems, some of the most relevant being its ability to smooth the original time series, and to separate the signal (i.e., trend and oscillatory components with different amplitudes) from the noise components. Therefore, SSA can be used to analyze and reconstruct smoother noise-free time series that can then be used for model forecasting. SSA is divided into two interconnected stages: decomposition and reconstruction of the time series. These stages are divided into two sets each, forming a total of four steps: embedding, singular value decomposition (SVD), grouping, and diagonal averaging. The complete algorithm for model fit is described in the following sub-section. Further details can be found in, e.g., [5,6,28]. 2.3.1. Decomposition In the first stage, the (univariate) time series is converted into a high-dimensional matrix called a trajectory matrix, which is then decomposed into the sum of rank-one matrices based on the SVD. (1) Embedding: Consider a non-zero time series YN={y1 , . . . , yn} with size N> 2. Let L( 1 <L<N) be an integer value called window length and K an integer such that the trajectory matrix includes all values; i.e., K=N−L+ 1. The embedding step is achieved by mapping the original time series into a sequence of Kvectors with length L: Yi= (yi,··· ,yi+L−1)T, 1 ≤i≤K. (3) Entropy 2020,22, 83 4 of 22 Then, the trajectory matrix X , that includes the vectors Yi , i= 1, . . . , K , in its columns can be written as: X= [Y1,··· ,YK]=(yij)LK i,j=1=      y1y2··· yK y2y3··· yK+1 . . .. . ..... . . yLyL+1··· yN      . (4) (2) Singular value decomposition: Let S=XXT , U1 , . . . , UL be the eigenvectors of S , and λ1≥ ··· ≥ λL , its corresponding eigenvalues. If d is the number of non-null eigenvalues of S , and considering Vi=XTUi√λi , we can decompose the trajectory matrix Xas: X= d ∑ i=1 Xi= d ∑ i=1pλiUiVT i. (5) The decomposition stage can be accomplished either by the eigendecomposition of XTX or by the SVD of X ( X=UDVT , D=diag(√λ1 , . . . , √λd) ). A comparison between both decompositions can be found in [29]. 2.3.2. Reconstruction In the second stage, after a separating signal from noise components, a diagonal averaging procedure is conducted in the matrices associated to the signal resulting into the sum of time series components that can then be interpreted as trend or oscillatory components: (1) Eigentriple grouping: This step consists of identifying the first r eigentriples associated with the signal and discarding the d−r eigentriples associated with the noise. Formally, let I=1, . . . , r and Ic=r+1, . . . , d . The goal of this step is to choose Isuch that the trajectory matrix can be written as: XI=∑ i∈IpλiUiVT i+e, (6) where eis the noise term. The number of eigentriples to conduct the reconstruction is often decided based on w-correlations. We shall say that two series Y(1) and Y(2) are approximately separable if all correlations between the rows and the columns of the corresponding trajectory matrices obtained from series Y(1) and Y(2) are close to zero. In [ 5 ] they considered other characteristics of the quality of separability; namely, the weighted correlation or w -correlation, which is a natural measure of deviation of two series Y(1) T and Y(2) Tfrom w-orthogonality: ρ(w) 12 =Y(1) T,Y(2) Tw kY(1) TkwkY(1) Tkw , (7) where kY(i) Tkw=rY(i) T,Y(i) Tw , i= 1, 2, and Y(1) T,Y(2) Tw=∑T t=1wty(1) ty(2) t with wt= min {t,L,T−t+1} . If the absolute value of the w -correlation is small, the two series are almost w -orthogonal. If the absolute value of the w -correlation is large, the series are far from being w -orthogonal and are, therefore, badly separable. Further explanation and intuition about this measure can be found in [5,28]. Other proposals for this choice were proposed by, e.g., [30,31]. (2) Diagonal averaging: Entropy 2020,22, 83 5 of 22 In this step, using anti-diagonal averaging on the matrices included in XI , the noise-free time series is reconstructed. First, the approximate trajectory matrix XI is transformed into a Hankel matrix. Let As={(l , k):l+k=s , 1 ≤l≤L , 1 ≤k≤K} and # (As) be the number of elements in As . The element e xij of the new Hankel matrix e Xis given by: e xij =∑ (l,k)∈As xlk #(As). (8) Next, the Hankel matrix e XI is transformed into a new series of dimension N , and the original time series YNcan be approximated by: e yi=(e xi1for i=1, . . . , L, e xLj for i=L+1, . . . , N,(9) where j=i−L+1. The reconstructed noise-fee time series can then be used for out-of-sample forecasting. 2.4. Robust SSA Despite knowing that SSA has shown to be superior to traditional model-based methods in many applications, the singular value decomposition (second step of the SSA algorithm) is highly sensitive to data contamination with outliers. Very few studies were made in order to access effects of outliers in SSA and to generalize this methodology [ 21 , 22 ]. A first attempt to robustify the SSA by considering an SVD based on a robust L1 norm [ 24 ] instead of the L2 norm used in the classical algorithm, was proposed by [ 23 ]. That robust generalization was compared with the classical SSA algorithm for model fit by these authors. In this subsection we review that robust SSA algorithm proposed by [ 23 ] and propose a new robust algorithm for SSA that considers the SVD based on the Huber function [ 25 ] and also propose an algorithm for robust SSA model forecasting. While the robust algorithms based on the L1 norm are very popular, they have difficulties in handling heavy tail outliers. The robust algorithms based on the Huber function combine the sum of squares loss and the least absolute deviation loss, that is, a quadratic on small errors, but grows linearly for large errors. As a result, the Huber loss function is not only more robust against outliers but also more adaptive for different types of data [ 32 ]. Further details and comparisons between the L1 and Huber loss functions, among others, can be found in [ 33 ]. The R source code is available upon request from the first author of this paper. 2.4.1. Robust SSA Based on the L1Norm The robust SSA algorithm proposed by [ 23 ] replaces the classical SVD based on the least squares L2 norm, by the robust SVD algorithm based on the L1 norm [ 24 ]. This robust SVD is performed iteratively, starting with an initial estimate of the first left singular vector U1 and leading to an outlier-resistant approach that also allows for missing data. The robust SVD based on the L1 norm is implemented under the function “robustSVD()” from the R package “pcaMethods”. 2.4.2. Robust SSA based on the Huber Function Here we propose a new alternative to robustify the SSA algorithm, where the least squares SVD in the step two is replaced by the robust SVD based on the Huber function [ 25 ]. The Huber loss function [34] can be defined as: Lδ(a) = (1 2a2if |a| ≤ δ δ|a|− 1 2δif |a|>δ, (10) Entropy 2020,22, 83 6 of 22 where δ is a parameter that controls the robustness level, and a smaller value of δ usually leads to more robust estimation. The robust SVD based on the Huber function is a special case of robust regularized SVD and can be obtained with the function “RobRSVD” of the “RobRSVD” R package, in the following way: RobRSVD (data, rough = TRUE, uspar = 0, vspar = 0). In this R implementation, the authors consider δ= 1.345, the value commonly used in robust regression that produces 95% efficiency for normal errors [ 35 ]. However, numerical studies suggested that the RobRSVD function is not very sensitive to the choice of δ[25]. More details about this robust SVD can be found in [25]. 2.5. Robust SSA Forecasting Algorithm The standard recurrent SSA forecasting algorithm assumes that a given observation can be written as a linear combination of the L− 1 previous observations [ 5 , 6 , 30 ]. The coefficients of those linear combinations in the classical SSA forecasting algorithm are obtained based on the left singular vectors, U , of the trajectory matrix X . This is valid for SSA because of the orthogonality of the vectors in U and of the full rank decomposition of X , which is not the case for the robust SVD algorithms because of their construction and specific properties. To overcome this limitation for the robust SSA algorithms and to be able to obtain out-of-sample forecasts using a robust SSA algorithm, a three stages approach can be conducted: (i) Use the robust SSA algorithm to obtain a robust approximation for the signal in the trajectory matrix; i.e., conduct the two stages of the robust SSA algorithms, decomposition (using the robust SVD algorithm) and reconstruction, to obtain the noise free (i.e., the signal) trajectory matrix e X; (ii) Apply the standard SVD to the matrix e X obtained in (i) and obtain U∇ j , the vector of the first L− 1 components of Uj and πj , the last component of the vector Uj , j= 1, ··· , r . Then, we can write the coefficient vector baas ba= (baL−1,··· ,ba1)0=1 1−γ2 r ∑ j=1 πjU∇ j, (11) where γ2=∑r j=1π2 j. (iii) The h-steps-ahead out-of-sample recurrent robust SSA forecasts byN+1 , . . . , byN+h , can be obtained as b yt=   ˜ yt, for t=1, ··· ,N ∑L−1 j=1bajb yt−j, for t=N+1, ··· ,N+h(12) where ˜ y1 , . . . , ˜ yN , are the fitted values for the reconstructed time series, as obtained from the robust SSA algorithm in (i). 2.6. Accuracy Measures There are several methods and measures for assessing model accuracy based on the behavior of model errors. Here, there are two types of errors: •Sample errors, called tuning errors; •Out-of-sample errors, called forecast errors. Typically, the root mean squared error (RMSE) is used as a criterion for accessing the precision of a model. The RMSE to investigate the quality of the model fit can be written as: RMSE =v u u t1 N N ∑ t=1 (yt−e yt)2, (13) Entropy 2020,22, 83 7 of 22 where yt are the observed values and e yt the fitted values by the considered model/algorithm (i.e., ARIMA, SSA, robust SSA). To investigate the forecasting accuracy, let us assume that the last g observations are used as a reference (i.e., as test set). Let N0=N−h−g . The RMSE to investigate the quality of the forecasting model can be written as: RMSE =v u u t1 g N ∑ t=N0+h+1 (yt−e yt)2, (14) where ytare the last gobserved values and e ytthe respective h-steps-ahead forecast values. 3. Results and Discussion In this section, comparisons are made between the classical ARIMA model, the classical SSA algorithm, and the robust SSA algorithms, in terms of computational time and accuracy for model fit and model forecast. These comparisons for decomposing and forecasting time series are done by considering a simulation example and the time series of six mutual investment funds. Table 1shows the descriptive statistics for the six mutual investment funds, including the minimum, maximum, and mean returns, being clear that Alaska Black is the fund that shows the largest variation and with the highest mean daily return. On the other end there are Gavea Macro and SPX Nimitz, which show the smallest variations among the considered funds, and low mean returns. In addition to the descriptive measures, Figure 1shows the behavior of the six investment funds over time. From these plots, it is possible to observe that all funds have an overall growing tendency, with similar patterns for Gavea Macro and SPX Nimitz. Brasil Capital Gávea Macro SPX Nimitz ADAM Estrategy Alaska Black APEX Long Biased 2014 2016 2018 2020 2010 2015 2020 2012 2014 2016 2018 2020 2017 2018 2019 2017 2018 2019 2014 2016 2018 1.0 1.5 2.0 2.5 1.0 1.5 2.0 2.5 3.0 1.0 1.5 2.0 2.5 100 150 200 250 300 1.0 1.1 1.2 1.3 1.4 1.5 1.0 1.5 2.0 2.5 3.0 Figure 1. Time series for the returns of the six mutual investment funds, ADAM Strategy, Alaska Black, APEX Long Biased, Brasil Capital, Gávea Macro and SPX Nimitz, from left to right and from top to bottom. The vertical axes show the quota values; i.e., the total net assets of a fund divided by the total number of quotas existing. Entropy 2020,22, 83 8 of 22 Table 1. Descriptive measures for returns of the six mutual investment funds. Investment Fund Minimum Mean Maximum Standard deviation ADAM Strategy −6.26% 0.05% 1.63% 0.0045% Alaska Black −29.62% 0.16% 9.80% 0.0240% APEX Long Biased −8.60% 0.07% 3.72% 0.0085% Brasil Capital −7.55% 0.07% 3.42% 0.0094% Gavea Macro −2.22% 0.04% 2.36% 0.0033% SPX Nimitz −1.92% 0.05% 1.42% 0.0030% 3.1. Model Fit The models/algorithms under comparison for model fit are: (i) ARIMA, (ii) SSA, (iii) robust SSA based on the L1norm (RLSSA), and (iv) robust SSA based on the Huber function (RHSSA). The parameters of the ARIMA model for each of the six mutual investment funds were estimated with the function “auto.arima” from the R package “forecast” [36]. For the SSA and robust SSA algorithms, there are two choices to be made by the researcher: (i) the window length L ; and (ii) the number of eigentriples used for reconstruction r . Three values of L were chosen for each time series, as defined in Table 2— L1=N/ 20, L2=N/ 2, and Lp —being the Lp obtained from the periodogram, based on the largest cycle for each time series [ 37 ] (i.e., about one trimester for ADAM Strategy, one semester for Alaska Black, one year for APEX Long Biased, one quadrimeter for Brasil Capital, one quadrimeter for Gavea Macro, and one quadrimester for SPX Nimitz), and N being the time series length. The choice of the number of eigentriples used for reconstruction r , for each of the considered window lengths and each of the time series, was done by taking into consideration the the w-correlations among components [ 5 ]. Figure 2shows the w-correlation matrices for each of the six mutual investment funds, considering an window length L=N/ 20, and Figure A1 of the appendix shows the w-correlation matrices for each of the six mutual investment funds, considering an window length L=N/ 2. The w-correlation matrices can be obtained with the function “wcor” of the R package “Rssa” [ 38 ] and the number of eigentriples r should be chosen in order to maximize the separability between signal and noise components; i.e., maximize the w-correlation among signal components, maximize the w-correlation among noise components, and minimize the w-correlation between signal and noise components. A summary of the number of eigentriples used for the reconstruction of each time series for each of the window length considered can be seen in Table 2. Since one of the objectives in SSA is to decompose the original time series into interpretable components such as trend and seasonality, plus the noise component that is then discarded, Figure 3 shows the original time series for the Alaska Black mutual investment fund, its trend component (sum of individual trend components), its seasonal component (sum of individual seasonal components), and its residuals (sum of the remaining components associated to noise), considering an window length L=N/ 20 = 33 and r= 12 eigentriples for reconstruction. Similar SSA decompositions for ADAM Strategy, APEX Long Biased, Brasil Capital, ADAM Strategy, Gavea Macro, and SPX Nimitz—considering the values of window length L1 and r1 eigentriples used for reconstruction, as defined in Table 2—can be found in Figures A2–A6 of the appendix, respectively. Entropy 2020,22, 83 9 of 22 Figure 2. W-correlation matrices for each of the six mutual investment funds, ADAM Strategy, Alaska Black, APEX Long Biased, Brasil Capital, Gávea Macro and SPX Nimitz, from left to right and from top to bottom, considering an window length L=N/20. Figure 3. Decomposition of the original time series for the Alaska Black mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 33 and r=12 eigentriples for reconstruction. Entropy 2020,22, 83 16 of 22 Table 12. Mean of the root mean square errors for model forecasting ( M= 1, 5, and 10 steps-ahead), computed for each of the four models, ARIMA, SSA, robust SSA based on the L1 norm, and robust SSA based on the Huber function, for the simulated data, based on 100 runs, using L=24 and r=5. M% of Method Cont. Shift ARIMA SSA RLSSA 0% - 1.685 0.125 0.245 5% yi+2 0.843 0.475 0.330 M = 1 10% yi+2 0.793 0.596 0.426 5% yi×5 3.960 8.461 0.358 10% yi×5 4.359 9.692 0.652 0% - 1.631 0.122 0.222 5% yi+2 0.984 0.475 0.307 M = 5 10% yi+2 0.768 0.586 0.413 5% yi×5 3.789 538.447 0.323 10% yi×5 3.853 17.670 0.720 0% - 1.381 0.127 0.244 5% yi+2 1.320 0.601 0.358 M = 10 10% yi+2 1.148 0.698 0.474 5% yi×5 3.486 22.695 * 4.015 10% yi×5 3.694 622.783 2.320 * 10% trimed mean. The mean value is 1.566 ×106. 4. Conclusions In this paper we considered the problem of model fit and model forecasting in time series. In particular, we analyzed six mutual investment funds. Following up on [ 23 ], who proposed a robust SSA algorithm by replacing the standard least squares SVD by a robust SVD algorithm based on the L1 norm [ 24 ] for model fit, we proposed another robust SSA algorithm where the robust SVD based on the Huber function is considered [ 25 ]. Moreover, we propose a forecasting strategy for the robust SSA algorithms, based on the linear recurrent SSA forecasting algorithm. Comparisons were made between the classical SSA algorithm, the robust SSA algorithms, and the classical ARIMA model, both in terms of computational time and accuracy for model fit and model forecast. Those comparisons were made by using daily observations of six mutual investment funds, and a synthetic data set where the time series were contaminated with outlying observations. For model fit of the six mutual investment funds, the best results were obtained for the SSA model when the window length L was set to be equal to the length of the time series divided by 20, or when the window length is defined as the length of the largest cycle in the time series. The ARIMA model and the robust SSA algorithms alternated for the second best performance. For model forecasting of the six mutual investment funds, the best overall performance was obtained for the classic SSA model considering a lower value for the window length, L=N/ 20 or Lp , followed closely by the ARIMA model and the robust SSA algorithm based on the L1norm. Based on the similarity between the results from the classic SSA model and the robust SSA algorithms, both for model fit and model forecasting, one may assume that the time series data from the six mutual investment funds had no or little data contamination. To access that hypothesis and to better illustrate the usefulness of the robust SSA algorithms, using a scenario with known and controlled outliers, a simulation study and its results were presented in this article. For both mild and and more extreme contamination scenarios, the robust SSA algorithms clearly outperformed the classical AMMI and SSA models, both for model fit and for model forecasting. Another important advantage of the robust SSA algorithms, because of their use of the robust SVD, is that they allow for missing values. Entropy 2020,22, 83 17 of 22 In terms of computational time, the SSA model gives the best performance, the robust algorithms being the most time consuming. A possible future development to reduce the computational time in the robust SSA algorithms is to consider a similar strategy as in [ 39 ], where a randomized SVD algorithm was used to speed up the SSA algorithm. The usefulness of the proposed approach, regarding the forecasting case, can be assessed based on forecasting competitions (e.g., [40]) or large scale forecasting studies (see, e.g., [41]). The methodology and results presented in this paper are of great generality and can be applied to other time series applications. Author Contributions: Conceptualization, P.C.R.; Formal analysis, P.C.R., J.P. and P.M.; Methodology, P.C.R. and M.K.; Software, P.C.R., J.P., P.M. and M.K.; Supervision, P.C.R.; Visualization, J.P. and P.M.; Writing—original draft, P.C.R., J.P., P.M. and M.K.; Writing—review and editing, P.C.R., J.P. and M.K. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Acknowledgments: The authors thank the associate editor and three anonymous reviewers for providing helpful suggestions which contributed to the improvement of the paper. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: ARIMA autoregressive integrated moving average SSA singular spectrum analysis SVD singular value decomposition RHSSA robust SSA algorithm based on the Huber function RLSSA robust SSA algorithm based on the L1norm RMSE root mean squared error Appendix A Figure A1. W-correlation matrices for each of the six mutual investment funds, ADAM Strategy, Alaska Black, APEX Long Biased, Brasil Capital, Gávea Macro and SPX Nimitz, from left to right and from top to bottom, considering an window length L=N/2. Entropy 2020,22, 83 18 of 22 Figure A2. Decomposition of the original time series for the ADAM Strategy mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 41 and r=17 eigentriples used for reconstruction. Figure A3. Decomposition of the original time series for the APEX Long Biased mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 80 and r=14 eigentriples used for reconstruction. Entropy 2020,22, 83 19 of 22 Figure A4. Decomposition of the original time series for the Brasil Capital mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 88 and r=12 eigentriples used for reconstruction. Figure A5. Decomposition of the original time series for the Gavea Macro mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 140 and r=12 eigentriples used for reconstruction. Entropy 2020,22, 83 20 of 22 Figure A6. Decomposition of the original time series for the SPX Nimitz mutual investment fund (top panel), with a trend component (sum of individual trend components, second panel), a seasonal component (sum of individual seasonal components, third panel), and a residual (sum of the remaining components associated to noise, bottom panel), considering an window length L=N/ 20 = 109 and r=8 eigentriples used for reconstruction. Appendix B A second synthetic dataset was obtained by generating random values from the following function and then transforming them into a time series: f(t) = cos (2πwt +φ)+e,t=1, ..., 100, with w= 3 / 8, φ=π/ 8 and e the noise generated from the N( 0, 0.1 ) (right-hand side of Figure A7). A total of 100 simulated time series were considered. 5% Additive Outliers 5% Multiplicative Outliers Synthetic Data 0 25 50 75 100 0 25 50 75 100 0 25 50 75 100 −1.0 −0.5 0.0 0.5 1.0 −1 0 1 2 3 −1 0 1 2 3 Figure A7. Synthetic data without contamination ( right ), data with 5% additive outliers ( left ), and data with 5% multiplicative outliers ( center ). The vertical axes show the simulated value and the horizontal axes show the index of the simulated observation. The data contamination was done in the same manner as described before. An example of 5% additive outliers scenario can be found on the left-hand plot of Figure A7, and an example of 5% multiplicative outliers scenario can be found on the central plot of Figure A7. The results for the root mean square errors for model fit, computed for each of the four models, ARIMA, SSA, robust SSA based on the L1norm, and robust SSA based on the Huber function, can be found in Table A1. Entropy 2020,22, 83 21 of 22 Table A1. Mean of the root mean square errors for model fit, computed for each of the four models, ARIMA, SSA, robust SSA based on the L1 norm, and robust SSA based on the Huber function, for the simulated data, based on 100 runs, using L=24 and r=2. % of Data Contamination Shift ARIMA SSA RLSSA RHSSA 0% - 0.1045 0.0097 0.0099 0.0104 2% yi+2 0.277 0.071 0.058 0.019 5% yi+2 0.351 0.113 0.096 0.032 10% yi+2 0.465 0.161 0.197 0.055 2% yi×5 0.279 0.108 0.026 0.018 5% yi×5 0.386 0.193 0.052 0.040 10% yi×5 0.484 0.338 0.075 0.098 References 1. Varga, G.; Wengert, M. A industria de fundos de investimentos no Brasil. Rev. Econ. Adm. 2011 ,10, 66–109. [CrossRef] 2. Maestri, C.O.N.M.; Malaquias, R.F. Exposition to factors of the investment funds market in Brazil. Rev. Contab. Financ. 2017,28, 61–76. [CrossRef] 3. Broomhead, D.S.; King, G.P. Extracting qualitative dynamics from experimental data. Phys. D Nonlinear Phenom. 1986,20, 217–236. [CrossRef] 4. Fraedrich, K. Estimating the Dimensions of Weather and Climate Attractors. J. Atmos. Sci. 1986 ,43, 419–432. [CrossRef] 5. Golyandina, N.; Nekrutkin, V.; Zhigljavsky, A. Analysis of Time Series Structure: SSA and Related Techniques; Chapman & Hall/CRC: New York, NY, USA, 2001. 6. Golyandina, N.; Zhigljavsky, A. Singular Spectrum Analysis for Time Series; Springer Science and Business Media: Berlin/Heidelberger, Germany, 2013. 7. Hassani, H. Singular spectrum analysis: Methodology and comparison. J. Data Sci. 2007,5, 239–257. 8. Hassani, H.; Zhigljavsky, A. Singular spectrum analysis: methodology and application to economics data. J. Syst. Sci. Complex. 2009,22, 372–394. [CrossRef] 9. Mahmoudvand, R.; Alehosseini, F.; Rodrigues, P.C. Forecasting mortality rate by singular spectrum analysis. RevStat-Stat. J. 2015,13, 193–206. 10. Mahmoudvand, R.; Rodrigues, P.C. Missing value imputation in time series using singular spectrum analysis. Int. J. Energy Stat. 2016,4, 1650005. [CrossRef] 11. Groth, A.; Ghil, M. Synchronization of world economic activity. Chaos: An Interdisciplinary. J. Nonlinear Sci. 2017,27, 127002. 12. Mahmoudvand, R.; Konstantinides, D.; Rodrigues, P.C. Forecasting mortality rate by multivariate singular spectrum analysis. Appl. Stoch. Models Bus. Ind. 2017,33, 717–732. [CrossRef] 13. Zabalza, J.; Qing, C.; Yuen, P.; Sun, G.; Zhao, H.; Ren, J. Fast implementation of two-dimensional singular spectrum analysis for effective data classification in hyperspectral imaging. J. Frankl. Inst. 2018 ,355, 1733–1751. [CrossRef] 14. Mahmoudvand, R.; Rodrigues, P.C.; Yarmohammadi, M. Forecasting daily exchange rates: A comparison between SSA and MSSA. RevStat-Stat. J. 2019,17, 599–616. 15. Mahmoudvand, R.; Rodrigues, P.C. Predicting the Brexit outcome using singular spectrum analysis. J. Comput. Stat. Model. 2019,1, 9–15. 16. Ge, M.; Lv, Y.; Zhang, Y.; Yi, C.; Ma, Y. An effective bearing fault diagnosis technique via local robust principal component analysis and multi-scale permutation entropy. Entropy 2019,21, 959. [CrossRef] 17. Sulandari, W.; Subanar; Lee, M.H.; Rodrigues, P.C. Indonesian electricity load forecasting using singular spectrum analysis. Energy 2020,190, 116408. [CrossRef] 18. Mahmoudvand, R.; Rodrigues, P.C. Prediction intervals for the vector SSA forecasting algorithm in a median based singular spectrum analysis. Comput. Math. Methods 2020. [CrossRef] 19. Reisen, V.A.; Molinares, F.F. Robust estimation in time series with long and short memory properties. Ann. Math. Inform. 2012,39, 207–224. Entropy 2020,22, 83 22 of 22 20. Rodrigues, P.C.; Monteiro, A.; Lourenço, V.M. A Robust additive main effects and multiplicative interaction model for the analysis of genotype-by-environment data. Bioinformatics 2016,32, 58–66. 21. Hassani, H.; Mahmoudvand, R.; Omer, H.N.; Silva, E.S. A preliminary investigation into the effect of outlier(s) on singular spectrum analysis. Fluct. Noise Lett. 2014,13, 1450029. [CrossRef] 22. Rodrigues, P.C.; Mahmoudvand, R. Correlation analysis in contaminated data by singular spectrum analysis. Qual. Reliab. Eng. Int. 2016,32, 2127–2137. [CrossRef] 23. Rodrigues, P.C.; Lourenço, V.M.; Mahmoudvand, R. A robust approach to singular spectrum analysis. Qual. Reliab. Eng. Int. 2018,34, 1437–1447. [CrossRef] 24. Hawkins, D.M.; Liu, L.; Young, S. Robust singular value decomposition. Natl. Inst. Stat. Sci. 2001 ,122, 1–12. 25. Zhang, L.; Shen, H.; Huang, J.Z. Robust regularized singular value decomposition with application to mortality data. Ann. Appl. Stat. 2013,7, 1540–1561. [CrossRef] 26. Brockwell, P.J.; Davis, R.A. Introduction to Time Series and Forecasting; Springer: New York, NY, USA, 1996. 27. Ripley, B.D. Time Series in R 1.5.0. R News, 2/2, 2–7. Available online: https://www.r-project.org/doc/ Rnews/Rnews_2002-2.pdf (accessed on 6 January 2020). 28. Rodrigues, P.C.; Mahmoudvand, R. The benefits of multivariate singular spectrum analysis over the univariate version. J. Frankl. Inst. 2018,355, 544–564. [CrossRef] 29. Ghil, M.; Allen, M.R.; Dettinger, M.D.; Ide, K.; Kondrashov, D.; Mann, M.E.; Robertson, A.W.; Saunders, A.; Tian, Y.; Varadi, F.; et al. Advanced spectral methods for climate time series. Rev. Geophys. 2002 ,40, 3.1–3.41. [CrossRef] 30. Mahmoudvand, R.; Rodrigues, P.C. A new parsimonious recurrent forecasting model in singular spectrum analysis. J. Forecast. 2018,37, 191–200. [CrossRef] 31. Rodrigues, P.C.; Mahmoudvand, R. A new approach for the vector forecast algorithm in singular spectrum analysis. Commun. Stat. Simul. Comput. 2020. [CrossRef] 32. Wen, Q.; Gao, J.; Song, X.; Sun, L.; Tan, J. RobustTrend: A Huber loss with a combined first and second order difference regularization for time series trend filtering. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, Macao, China, 10–16 August 2019; pp. 3856–3862. 33. Bouwmans, T.; Aybat, N.S.; Zahzah, E. Handbook of Robust Low-Rank and Sparse Matrix Decomposition: Applications in Image and Video Processing; CRC Press: New York, NY, USA, 2016. 34. Huber, P.J. Robust estimation of a location parameter. Ann. Math. Stat. 1964,35, 73–101. [CrossRef] 35. Huber, P.J.; Ronchetti, E.M. Robust Statistics; Wiley: Hoboken, NJ, USA, 2009. 36. Hyndman, R.J.; Khandakar, Y. Automatic time series forecasting: The forecast package for R. J. Stat. Softw. 2008,26, 1–22. 37. de Carvalho, M.; Rua, A. Real-Time Nowcasting the US Output Gap: Singular Spectrum Analysis at Work. Int. J. Forecast. 2017,33, 185–198. [CrossRef] 38. Golyandina, N.; Korobeynikov, A.; Shlemov, A.; Usevich, K. Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package. J. Stat. Softw. 2015 ,67. Available online: https://www.jstatsoft. org/article/view/v067i02 (accessed on 6 January 2020). [CrossRef] 39. Rodrigues, P.C.; Tuy, P.G.S.E.; Mahmoudvand, R. Randomized singular spectrum analysis for long time series. J. Stat. Comput. Simul. 2018,88, 1921–1935. [CrossRef] 40. Hyndman, R.J. A brief history of forecasting competitions. Int. J. Forecast. 2020,36, 7–14. [CrossRef] 41. Papacharalampous, G.; Tyralis, H.; Koutsoyiannis, D. Comparison of stochastic and machine learning methods for multi-step ahead forecasting of hydrological processes. Stoch. Environ. Res. Risk Assess. 2019 , 33, 481–514. [CrossRef] c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).