3.2.1 Tests for Univariate and Multivariate Normality

It is known [²¹] that PCA (implicitly) and FA (explicitly) require a normally distributed multivariate sample in order to produce completely reliable results, i.e., they will only produce uncorrelated and independent components if the sample data have no higher order statistics beyond the variance.

Thus, if the samples do not fulfill these conditions, we will be prompted to use a more suitable technique such as ICA to uncover the underlying sources in a non-Gaussian sample. Therefore, we first tested the univariate normality (UVN) of each individual series, since ICA requires that not more than one of the observed signals (the returns on equities) be non-Gaussian.

Tables 1 to 4 present the descriptive statistics up to the fourth moment of the four databases used in this study. We can observe that the skewness and the kurtosis of practically all the stocks differs from those of the Gaussian distribution.

We also carried out the Jarque-Bera test for UVN on the four databases, rejecting the null hypothesis of normality at 5% of probability for all the stocks in the daily databases, but not rejecting it for only one stock in the weekly databases that was normally distributed. The last two columns of the Tables 1 to 4 present the results of the Jarque-Bera test.

We used two classical alternatives for assessing the multivariate normality (MVN) tests: the Mardia [³³] and the Henze-Zirkler [¹⁸] MVN tests. Mardia’s test is based on the multivariate skewness and kurtosis of the sample. Henze-Zirkler’s (H-Z) test considers a measure of the distance between the characteristic function of the MVN and the empirical one, where the computed statistic will be lognormally distributed, if the data is multivariate normal. Both techniques have shown very good performance in measuring the MVN against other classic and newer alternatives, as [³⁴] remark in their study.

We performed two tests following the accepted criterion of applying more than one MVN test when assessing this property of a sample.⁵ Our results with both tests reject the null hypothesis of MVN at 5% of probability for all the databases. Tables 5 and 6 present the results of Mardia’s and H-Z’s tests, respectively.

We extended this analysis by making an experiment concerning the horizon of Mardia’s test, i.e., we ran the test using different numbers of observations so as to check the multivariate normality in different scenarios. The results showed that from 101 observations on, inclusive, the sample is non-Gaussian according to the three statistics.

On the basis of the foregoing results⁶, we cannot accept as completely reliable the outcomes of techniques assuming the multivariate normality of data such as PCA and FA; thus, we are led to the application of more suitable techniques like ICA. In fact, this part of our investigation represents an important, but in most cases ignored, aspect in empiric studies that uses classic multivariate techniques to extract the pervasive factors; since in many cases the MVN is assumed but not tested, the results and conclusions may be flawed.

In addition, the assumption done in the ICA models, is that the third and fourth moments differ significantly from the values of a Gaussian distribution.

In addition, the tests of normality are based on checking this assumption. In particular the non-linearities used for the implementation of the experiments in this paper, guaranteed the presence of high order interactions from the Taylor expansion, and therefore the presence of moments of all orders.

3.2.2 Estimation of the ICA Model

In order to estimate the ICA model in expression (2), we used the ICASSO methodology [²⁰], which is based on the FastICA algorithm [²²]⁷. According to the foregoing authors, the FastICA algorithm is based on a fixed-point iteration scheme for finding the local extrema of the objective functions. The basic iteration for the vector w for each IC obtained by this method is:

where the nonlinearity g can be almost any smooth function such as:

and g’ is the derivative of g(.).⁸

The final vector gives one of the ICs as a linear combination in y = w^T z. The specific resulting algorithm depends both on the estimation principle used and the approach selected to estimate several numbers of ICs, i.e., the nonlinearity and the decorrelation method chosen. In [²¹], the authors state that by setting the options, nonlinearity tanh (hyperbolic tangent) and symmetric approach, one can obtain a good estimation of the ICA model; this would be equivalent to performing the three estimation approaches at the same time.

In addition, the positive kurtosis obtained in the multivariate normality tests leads us to use the hyperbolic tangent function.

Furthermore, as reported in [¹⁴], the best trade-off for estimating the ICA model, from statistical performance and computational load perspectives, is represented by the FastICA algorithm with symmetric orthogonalization and tanh nonlinearity estimation. In our study we followed these specifications.

The election of the ideal number of ICs to estimate still represents an unsolved problem.

Although in ICA literature we can find diverse criteria to determine this number, in most cases it is actually chosen by trial and error without any theoretical basis. One alternative is to reduce the number of dimensions in the whitening pre-processing stage, considering some criteria from among those used in PCA or FA, and to estimate the same number of ICs. For the sake of comparison with our previous study, we use the same test window, which ranges from two to nine components.⁹

As stated by [²⁰], one problem that the ICA estimation presents is that the reliability of the estimated ICs is not known since the results are stochastic, i.e., they might be dissimilar in different runs of the algorithm.

Thus, the results of a single run of the FastICA algorithm could not be completely trusted and an additional analysis of the reliability of the estimation should be performed. In this context, reliability has two aspects the algorithmic and the statistical. According to the former authors, ICASSO methodology represents an alternative for dealing with this problem, since it ensures the algorithmic and statistical stability and reliability of the estimated components by running the FastICA algorithm many times, using different initial conditions and/or a differently bootstrapped data set.

Following [²⁰], ICASSO first runs the FastICA algorithm M times on data set X=[x1,x2,…xN], composed of N samples of k vectors; then, ICASSO forms clusters with the ICs produced in each run according to their similarity. Mutual similarities between estimates are computed, using the absolute value of their linear correlation coefficient as the measure of similarity:

These elements form the similarity matrix, which can be obtained by:

where, Σ is the covariance matrix of dataset x, and W^ is the estimates of demixing matrices W^i from each run i=1,2,…,M gathered in a single matrix:

According to [¹⁹], reliable estimates of ICs correspond to tight clusters, since they agglomerate estimates generated by many runs of the algorithm which are similar, even when the initial values and datasets for the estimation have been changed. Conversely, estimates which do not belong to any cluster are considered unreliable estimates. The centrotype of each cluster is considered a more reliable estimate than that generated by any single run.

Besides the previously declared parameters for FastICA, there are some additional parameters to set when using ICASSO, such as the resampling mode, number of resampling cycles (M) and number of clusters (L). In order to ensure both statistical and algorithmic reliability, in our study we used both resampling modes, i.e., each time the dataset was bootstrapped and the initial conditions of the algorithm were randomized. We used the default number of resampling cycles fixed by the software, i.e., 30, and we set the number of clusters according to the number of ICs (m) estimated in each experiment in order to obtain squared mixing (A) and demixing (W) matrices.

The demixing matrix (W) computed by ICASSO corresponds to the centrotypes of each cluster as well, representing a more reliable estimate than that produced by a single run of FastICA; however, they are not strictly orthogonalized. In the context of our research where we need to obtain orthogonalized ICs, we will have to make an orthogonalization procedure in a later step.

Consequently, we first took the demixing matrix (W) produced by ICASSO, then we computed the mixing matrix:

and the matrix of independent components or sources:

3.2.3 Ranking and Orthogonalization of the Independent Components

The ICA model does a decomposition by means of a criterion related to statistical independence, which does not allow to order in a natural way the components and thus the residual. The criterion presented in this section is one criterion that has sense in the application at hand. In contrast with the case of linear regression or PCA, where the driving noise is easy to identify, because it is a residual obtained after the components of maximum variance are determined, in the case of ICA such an interpretation will not be natural. Because of this, in the literature about ICA it is not clearly specified the difference between the components and the residual, and therefore the results are usually presented as a complete projection in the space statistically independent components.

Then, next we ordered the independent components in terms of their explained variability by means of the criterion proposed by [¹²]. This criterion ranks the ICs according to the amount of variance of the stocks that explains each one of them, thus we obtain a ranked matrix of independent components (S^r), as well as sorted mixing (A^r) and demixing matrices (W^r).

Finally, we orthogonalized the matrix of ICs by means of the following process of transformation:

where V is a transformation matrix to decorrelate the matrix of sorted independent components, and S^o represents the matrix of orthogonalized ICs.

3.2.4 Extraction of Underlying Systematic Risk Factors Via ICA

In each one of the four databases, we computed eight multifactor models in order to extract a window from two to nine independent components. Then, we proceeded to reconstruct the original variables according to the generation process of expression (1), including the inverse of the transformation matrix V in order to orthogonalize the mixing matrix A as well:

The reproduced values were very similar to the observed series for greater part of the equities in all the datasets, which indicates that the generative multifactor model performed by ICA was effective. However, stocks such as GMODELO, CEMEX, SORIANA and GCARSO were not very well reconstructed, especially in the cases of daily returns and excesses, due to the high volatility they presented during the studied period. To save space, we only present the line plots for the first five stocks appearing in the returns and excesses observed and reproduced from each database.

Figures 1 to 4 present the results of the case when we extracted nine underlying factors; the reconstruction performance is evident.¹⁰ An interesting fact of the ICA algorithm is that it captures the global interaction between stocks, independently of the non-stationarity of the joint behavior. That is, the required assumption in the model is that there are independent sources that are mixed by a matrix W.

Note: Logarithmic returns of the first five stocks observed in each database and their respective reconstructions using the estimated ICA model. Stock symbols of the stocks presented appear above each line plots.

Fig. 1. Line plots of the observed and reproduced stocks 

If the matrix does not change, the ICA algorithm will give an estimation, and therefore, given that the matrix does not change, it will impute the components of volatility to some of the non-observable factors.

3.2.5 Independence Test

In order to test the independence of the computed ICs, we ran the Hilbert-Schmidt Independence Criterion (HSIC) test [¹⁵]¹¹, which tests whether random variables X and Y are independent based on a sample of observed pairs (x_i, y_i). The results of our independence tests confirmed the statistical independence, between each pair of components estimated from the weekly and daily databases.

3.2.6 Econometric Contrast

We carried out an econometric contrast under a statistical approach to the Arbitrage Pricing Theory (APT) using the underlying systematic risk factors extracted via ICA. The APT’s pricing equation is expressed as follows:

In the same outline that in [²⁸], λ₀ represents the riskless interest rate, λ_k the risk premium for each kind of systematic risk factor, and β_k the exposures to each type of systematic risk. We tested the former expression by way of an average cross-section methodology estimating the coefficients by ordinary least squares (OLS) in the following regression model:

We used again the two-stage methodology for the econometric contrast of the APT used in our aforementioned study [²⁸], which is explained as follows: In the first stage, we estimated the betas to be used in expression 18 from the scores of the extracted factor. In the second stage, we estimated the lambdas. In the first stage we estimated the betas by regressing the factor scores obtained by ICA as a cross-section on the returns and excesses. In order to improve the efficiency of the parameter estimates and to eliminate autocorrelation in the error terms of the regressions, we used weighted least squares (WLS) to estimate the entire system of equations at the same time.

The results of the regressions in the four databases were very good, producing, in almost all cases, statistically significant parameters, high values of the R² coefficients and results in the Durbin-Watson test of autocorrelation, which lead us to the non-rejection of the null hypothesis of no-autocorrelation. In the second stage we estimated the lambdas or risk premia in expression 17 by regressing the betas obtained in the first stage as a cross-section on the average returns and excesses, using ordinary least squares (OLS).

In order to avoid the econometric problems of heteroskedasticity and autocorrelation in the residuals of the model estimated through OLS, we corrected it by means of the Newey-West heteroskedasticity and autocorrelation consistent covariance estimates (HEC). Additionally, we verified the normality in the residuals by carrying out the Jarque-Bera test of normality.

In order to accept the APT pricing model, we require the statistical significance of at least one parameter lambda different from λ₀, and the equality of the independent term to its theoretic value, i.e., the average returns, in the models expressed in returns:

and zero, in the models expressed in excesses of the riskless interest rate:

We used Wald’s test to confirm these equalities.

In Table 7, we present a summary of the results of the econometric contrast for the four databases. In general, the results of the explanation power, the adjusted R-squared (R²*), the statistical significance of the multivariate test (F), and the Jarque-Bera normality test of the residuals are very good in almost all the contrasted models. The univariate tests for the individual statistical significance of the parameters (statistic t) priced from one to five factors exclusive of λ₀ in the weekly and daily databases, thus giving evidence in favor of the APT in 27 models.

Nevertheless, only four models fulfilled both the statistical significance of the parameters and the equality of the independent term to its theoretic value, in addition to the fulfilment of normality in the residuals.

The referred models appear marked in Table 7, where we used the same methodology of presentation and analysis of the results as in our preceding paper [²⁸].