Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
37,121 characters · 5 sections · 55 citation commands
Forecasting in Big Data Environments: an Adaptable and Automated Shrinkage Estimation of Neural Networks (AAShNet)
\setlist{noitemsep} \onehalfspacing
\thispagestyle{empty}
\centerline{ABSTRACT}
{\bf Key Words:} Nonlinear Shrinkage Estimation, Gradient-based Hyperparameter Optimization, High-dimensional Nonlinear Time Series, Neural Networks \\ {\bf JEL classification:} C45, C51, C52, C53, C61.
An important step in designing modern predictive models is to cope with high-dimensional data, presenting large numbers of (cor)related variables and complex properties. “Big data" is both an increase in the number of samples collected over time, and an increase in the number of potential explanatory variables and predictors. When dimension grows, the specificities of high-dimensional spaces and data must then be taken into account in the design of predictive models. While this is valid in general, its importance is heightened when using nonlinear tools such as artificial neural networks. Most nonlinear models involve more parameters than the dimension of the data space which may result in a lack of identifiability, lead to instability, and overfitting (huber2011;fromstatisticstoneuralnetworkstheoryandpatternrecognitionapplications1994; moody1991). Selection of significant predictors, and model complexity are the key tasks of designing accurate predictive models in data-rich environments.
Feature extraction and feature selection are broadly the two main approaches to dimensionality reduction. Extraction transforms the original features into a lower dimensional space preserving all its fundamentals. Feature selection methods select a small subset of the original features without a transformation. Extraction methods include principal component analysis - pearson1901; eckart1936; factor analysis - spearman1904; canonical correlations analysis - hotelling1936, and several others\footnotemark\footnotetext{cunningham2015 surveyed the literature on linear dimensionality reduction in their work.}.). Feature selection is accomplished by such methods as Ridge - hoerl1970; LASSO - tibshirani1996 and Elastic Net - zou2005).
In this work, our main focus is on feature selection techniques. We apply shrinkage approaches (usually referred to as regularization in machine learning literature). We embed feature selection in the backpropagation algorithm as part of its overall operation. Accordingly, we extend our loss function to include $L_1$ norm for the weights of the dense network, and $L_2$ norm for the weights in the skip-layer. The dense network corresponds to a multilayer neural network, whereas the skip-layer denotes the direct connection from each of the input variables to each of the output variables, which is similar to a linear regression model.
Shrinkage is an implicitly embedded feature selection. It is an example of model selection since only a subset of variables contributes to the final predictor. It has frequently been observed that $L_1$ shrinkage produces many zero parameters, leading to some features being dropped and a sparse model. Only those parameters whose impact on the empirical risk is considerable appear in the fitted model ng2004. Shrinkage is a proper means of controlling complexity in the nonlinear component. From an optimization point of view we have a neural network learned/estimated by LASSO. This prevents hidden units from getting stuck near zero and/or exploding weights.
Simultaneously, we employ the $L_2$ shrinkage on the skip-layer connections (linear part of the model), in order to penalize groups of parameters, and encourage the sum of the squares of the parameters to be small. Therefore we will not drop specific features from linear component, making it possible to interpret the marginal impact of predictors on the target variable. It is worth mentioning that the linear part of the model can be interpreted as a Ridge regression.
There are other benefits to shrinkage/regularization. Empirically, penalizing the magnitude of network parameters is also a way to reduce overfitting and to increase prediction accuracy ng2004. This is especially true in the state-of-art models, such as deep learning models with large number of parameters. Our proposed algorithm combines the neural network’s advantage of describing the nonlinear process with the superior accuracy of feature selection that is provided by a penalized loss function that combines $L_1$ and $L_2$ norms.
Many studies have suggested neural networks as a promising alternative to linear regression models. Empirical evidence on out-of-sample forecasting performance is, however, mixed. It is challenging to determine linear or nonlinear components. Linearity tests do often suggest that real world series are rarely purely linear or nonlinear.
We consider the possibility that the series $(y_t)$ contain both a linear component, $(\mathcal{L}_t)$, and a nonlinear component $(\mathcal{N}_t)$.
Neural network alone is not best suited to handle both linear and nonlinear components, especially when the linear component is superior to the nonlinear component.
Two different approaches to model and forecast series with both linear and nonlinear patterns are available. The first approach is a two step methodology to combine linear time series models and neural network models. In this approach, the first step residuals are obtained from the fitted linear model $\hat{e_t}= y_t - \hat{\mathcal{L}_t}$. In the second step a nonlinear model (e.g., GARCH, neural nets) is trained on the residuals of the first step. In principle, this “hybrid” two step approach can provide superior predictions when both the linear and neural network model are well specified. In practice, however, two types of model specification errors are introduced without an ability to assess their mutual impact.
The alternative approach that we are proposing in this paper models both linear and nonlinear components adoptively. It is based on a neural network with skip-layer connections including both linear and nonlinear structures.
The rest of the paper is organized as follows. Section II provides the basic framework of the proposed model. In Section III we investigate proper estimation of shrinkage hyperparameters and introduce gradient-based techniques based on reverse-mode automatic differentiation (RMAD) to accomplish this. Section IV presents an application to US financial returns. Section V contains some concluding remarks.
In this study, we examine a feedforward neural network with one hidden layer, known as a dense network. Neural network models can be seen as generalizations of linear models, when one allows direct connections from the input variables to the output layer with a linear transfer function\footnotemark \footnotetext{Using linear function for the output unit activation function (in conjunction with nonlinear activations amongst the hidden units) allows the network to perform a powerful form of nonlinear regression. So, the network can predict continuous target values using a linear combination of signals that arise from one layer of nonlinear transformations of the input.}, that we refer to as the skip-layer. The model is expressed as
where $\Phi$ describes the network by a vector function. We associate subscript $i$ with the input layer, subscript $j$ with the hidden layer, and subscript $k$ with the output layer. $x_{it} = (x_{1t},x_{2t},...,x_{mn})$ is the value of the $i$th input node, which can be a constant input representing biases, a matrix of lagged values of $y_t$ and some exogenous variables. $\phi_j(.)$ and $J$ are activation functions and number of neurons used at the hidden layer. A single-hidden-layer neural network with skip-layer connections is shown in Figure (ref). A network with only one hidden layer and skip-layer connections has three sets of weights: those for direct connections between the inputs and the output ($w_{ik}$), those connecting the inputs to the hidden layer ($w_{ij}$), and those connecting the output of the hidden layer to the final output layer($w_{jk}$).
First term in Eq.(ref) represents a linear regression term. The second term, denoting the dense network of the two layers, hidden and output, is usually referred to as a multi-layer perceptron in the literature. It has been shown to be able to perform well with nonlinear complex data. A greater capacity of the dense network, compared to the skip-layer, is realized by stacking two layers, enabling it to model more complex data. A differentiable nonlinear activation function $\phi$ is used in the hidden units. $\varepsilon_t$ is a random disturbance term which captures all other factors influencing $y$ than the $x$. A linear component term moves the model in the linear direction. This aids statistical interpretation and unravels the structure behind the network, otherwise left to a black box. This simultaneous approach has the advantage, when we apply shrinkage techniques to estimate network parameters for an essentially linear process, of pruning the hidden neurons.
Estimation of network elementary parameters based on prediction error minimisation is known as training/learning. The most common cost/risk function is the mean squared prediction error (MSE), $E = \frac{1}{n}\sum_{t=1}^n (y_t - \hat{y}_t)^2$. Given target values $y_t$ and network estimated outputs $\hat{y}_t$ error functions are obtained for each parameter set, followed by tuning of the parameters.
The error surface becomes increasingly complicated with the number of input variables and network parameters. It is common to employ the conventional feed-forward neural network, trained with the popular and revolutionary gradient-descent-type algorithm known as backpropagation. The backpropagation algorithm was first introduced by bryson1979 and popularized in the field of artificial neural network research by werbos1988 and rumelhart1986. Error function's sensitivity to network parameters is assessed via Gradient Descent optimization. Gradient is normally defined as the first order derivative of the error function with respect to each of the model parameters. Working out the gradients can be performed in a completely mechanical way known as Automatic Differentiation baydin2017automatic. AD employs the Jacobian matrix of gradients for each parameter $w_i$ to identify directions that decrease the height of the error surface (see Appendix). In fact backpropagation is only a specific case of reverse-mode AD that is applied to an objective function errors as functions of model parameters. The weight adjustment is given by
Where the constant $\eta$ is the learning rate (step size) for updating elementary parameters, its value falls between zero and one. By iteratively repeating this mechanism, the network can be trained in a way that converges to the optima. The set of new elementary parameters are repeatedly presented to the network until the error value is minimized. Around the optimum point, all the elements of the gradient would be very small, leading to tiny changes in new parameters.
We add the $L_1$ and $L_2$ penalties in training our modelto the loss function $\widetilde{E}(.)$, the original MSE. The following optimization problem is used for training:
where the regularization term $\Omega(\bm {w, \lambda})$ is a combination of the L1 norm and the L2 norm of the parameter vector. $\lambda$ sets the impact of shrinkage on the loss, with larger values resulting in more penalization. Using the regularized objective causes the training procedure to be inclined to smaller parameter values; unless larger parameters considerably improve the original error value (MSE). Assuming a fixed $\lambda$, to learn $w^*$, we only need to include the derivative of $\Omega(\bm {w, \lambda})$ in our derivatives:
Where $\Delta$ is the gradient of the regularized loss function. $\lambda > 0$ is proportional to complexity of the model but is not a parameter that appears in the model. It is a hyperparameter. In the next section, we explain the impact of hyperparameters and elaborate on our procedure for tuning them.
We employ $L_1$ and $L_2$ shrinkage on the parameters of the dense network and skip-layer, respectively; as is depicted by following optimization problem:
which can be realized by iteratively adjusting the parameters using the updating rules below
Where $\lambda_1$ and $\lambda_2$ are non-negative values known as shrinkage hyperparameters. $L_1$ sparsity norm and $L_2$ smoothing norm are two closely related regularizers that can be used to impose a penalty on the complexity of the model that is to be learned. Shrinkage estimation of the model can be seen as an implementation of Occam's razor, introducing a controllable trade-off between fitting data and model complexity, enabling us to have models of less complexity with adequate generalization capability. Regularization in neural networks limits the magnitude of network parameters by adding a penalty for weights to the model error function. In this study, $L_2$ shrinkage penalizes parameters in skip-layer connections by adding sum of their squared values to the error term. $L_1$ shrinkage penalizes parameters in the dense network to encourage the topology of the learned network to be sparse. The relative importance of the compromise between finding small weights and minimizing the original risk function depends on the size of $\lambda$.
To use $L_2$ shrinkage, we add a $\lambda_2 w$ term to the gradient as the derivative of $w^2$ is $2w$. $L_2$ shrinkage works with all forms of learning algorithms, but does not provide implicit feature selection. The derivative of the absolute value of $w$ is $w/|w|$, however $L_1$ norm is not differentiable at zero and hence poses a problem for gradient-based methods.
The problem can be solved using the exact gradient, which is discontinuous at zero. We can also solve the problem by the smooth approximation approach which will allow us to use gradient descent. To smooth out the $L_1$ norm using an approximation, we use $\sqrt{w^2 +\epsilon}$ in place of $|w|$ , where $\epsilon$ is a smoothing parameter which can also be interpreted as a sort of sparsity parameter. When $\epsilon$ is large compared to $w$, the expression $w + \epsilon$ is dominated by $\epsilon$ and taking the squared root yields approximately $\sqrt{\epsilon}$. lee2006
The major drawback of shrinkage is that it introduces additional hyperparameters. In practice we have two set of parameters: model elementary parameters (network weights and biases), and learning algorithm hyperparameters (magnitude of $L_1$ and $L_2$ penalties, and learning rate). We would ideally like to determine these hyperparameters to get optimal generalization\footnotemark \footnotetext{Generalization means building a model on one set of training data and hope that it makes effective predictions on a different set of test data.}. As opposed to elementary parameters, these hyperparamters cannot be directly trained by the data. Whereas the elementary parameters specify how to transform the input data into the desired output, the hyperparameters define how our model and algorithm are actually structured.
The performance and robustness of neural networks relies to a large extent on hyperparameters. Tuning these hyperparameters not only makes the investigation of methods difficult, but also hinders reproducibility (bergstra2011algorithms). Transparent tuning of hyperparameters can be part of an Hyperparameter Optimization (HPO), as an outer loop in training procedures.
The de-facto na\"ive approach of searching through combinations of potential values of hypergradients and choosing the one that performed the best (a.k.a. grid search) is very time-consuming and becomes quickly infeasible as the dimension of hyperparameter space grows. In many practical applications manually searching the space of hyperparameter settings is tedious and tends to lead to unsatisfactory outcomes. bergstra2012 show empirically and theoretically that random search more efficient than grid search. Statistical techniques such as cross-validation wahba1990, bootstrapping efron1994, and Bayesian methods mackay1992 can also assist in determining hyperparameters.
HPO must be guided by some performance metric, typically measured by cross-validation (CV) on the training set, or evaluation on a held-out validation set. The rationale behind CV is to split the data into the training samples used for learning the algorithm, and the validation samples (one or several folds) for estimating the risk of each algorithm and for evaluation of its performance. CV consists of averaging several hold-out estimators (folds) of the risk corresponding to different splits of the data, and selecting the algorithm with the smallest estimated risk. Within each fold, hyperparameters are fixed and we only estimate model elementary parameters. The validation samples play the role of new unseen data as long as the data are i.i.d.\footnotemark \footnotetext{This assumption can be relaxed. see: chu1991.} For a general description of the CV see geisser1975, and arlot2010 for a comprehensive review on cross-validation procedures and their applications in different algorithms and frameworks. Several studies such as rivals1999 show cases in which CV performance is less than satisfactory.
Recently, automated approaches for estimation of hyperparameters have been proposed which can provide substantial improvements and transparency. Although one may also “hyperparameterize” certain discrete choices in design of the model (e.g. number of hidden units), we focus only on the continuous hyperparameters in this work. There are a number of gradient-free automated optimization methods (hutter2011; bergstra2011; bergstra2013; snoek2012), all of which rely on multiple complete training runs with varied fixed hyperparameters. Hyperparameters are chosen to optimize the validation loss after complete training of the model parameters.
Gradient-based HPO approaches, proposed by larsen1996 and andersen1997, emerged in the 1990s. We can distinguish two main approaches of gradient-based optimization: Implicit differentiation and iterative differentiation.
Implicit differentiation, first proposed by larsen1996, computes the derivative of the cost $L_{valid}$ with respect to $\lambda$ based on the observation that, under some regularity conditions, the implicit function theorem can be applied in order to calculate the gradients of the loss function. In particular, the cost function is assumed to smooth and converge to local minima. The inner optimization $w(\lambda) \in \textup{argmin}_{w}L_{train}$ can be characterized by the implicit equation $\nabla_w L_{train} = 0$. bengio2000gradient derived the gradients for unconstrained cost function and applied the algorithm to $L2$ shrinkage for linear regression. The method has also been used to find kernel parameters of Support Vector Machines keerthi2007efficient. pedregosa2016hyperparameter proposes HOAG which uses inexact gradients, allowing the gradient with respect to hyperparameters to be computed approximately.
In iterative differentiation, first proposed by domke2012generic, the gradient for hyperparameters are calculated by differentiating each iteration of the inner optimization loop and using the chain rule to aggregate the results. However, the problem with this reverse-mode approach is that one must retain the entire history of elementary parameter updates, making a na\"ive implementation impractical due to memory constraints. Reverse-mode differentiation requires intermediate variables to be maintained in the memory for the reverse pass and evaluation of validation loss needs hundreds or thousands of inner optimization iterations. maclaurin2015 later extended this for setting of stochastic gradient descent via reverse mode automatic differentiation of validation loss.The burden of storing the entire training trajectory $w_1,\cdots,w_T$ is avoided by an algorithm that exactly reverses SGD with momentum to compute gradients with respect to all training parameters, only using a relatively small memory footprint, making a solution feasible for large-scale big data machine learning problems.
We defined the updating rule for elementary parameters as $w_{t+1} = w_t - \eta \nabla L_{train}$ where $L_{train} = \widetilde{E}(w_t|\lambda, X_{train})$ is the regularized loss value on train data. To calculate hypergradients we rely on the unregularized loss function, that is $L_{valid} = E(w_t|\lambda, X_{valid})$, as the actual generalization performance of the model, on unseen data points, does not directly depend on regularizers; otherwise the model with no regularization would be always selected:
There are cases where SGD can become very slow. The method of momentum is designed to accelerate learning, especially in the face of high curvature, small but consistent gradients, or noisy gradients Goodfellow2016. We modify our training (Algorithm (ref)) to include a velocity variable $v$ storing the momentum by calculating exponentially decaying moving average of past gradients.
where $\gamma_t$ is the momentum decay rate. The training procedure starts with elementary parameters velocity $v_1=0$ and $w_1$ and ends with $v_T$ and $w_T = w_{T-1} + \eta_{T-1}v_{T-1}$. Algorithm (ref) is then used to calculate the gradients of validation loss with regard to the hyperparameters.
The velocity $v_t$ is needed to reverse the path, otherwise without momentum, $g_t$ and $\eta_t$ alone would not be able to recover $w_{t-1}$. Notice that the loss of information caused by finite precision arithmetic in computers leads to failure of this algorithm. For this reason, we need to store the bits lost in $v_t$ when multiplied by $\gamma_{t}$.
Given this powerful gradient-based mechanism for finding hyperparameters, a natural extension to our model is to introduce a hyperparameter $\alpha$ denoting the contribution of skip-layer and dense-network in producing predictions with higher generalization. That is to say, our model can be reformulated as:
where $\alpha$ assumes a value between zero and one. Appreciating that the skip-layer and the dense network have unbalanced effect on the outcome, one can see how this may result in faster convergence of training procedure. More importantly, $\alpha$ can be interpreted as the activation of skip-layer and dense network and can point to linearity or nonlinearity components.
Research into modelling and forecasting financial returns has a long history. Several models are described in tsay2005 and campbell1996 that attempt to explain return time series using linear combinations of one or more financial market factors. The most widely studied single factor model is the capital asset pricing model (CAPM) of sharpe1964 and lintner1965 that relates the expected return of equities to the expected rate of return on a market index (such as the Standard and Poor’s 500 Index). The empirical performance of CAPM is poor as it cannot explain the behaviour of asset returns, see fama2004. This failure is perhaps due to the absence of multiple factors. Arbitrage pricing theory (APT) is a general model proposed by ross1976 to account for these deficiencies. APT presents a linear approximate model of expected asset returns based on an unknown number of macroeconomic “factors” or market indices. The relationship between the factors and historical returns is routinely determined linearly.
Return time series present characteristics such as comovement, nonlinearity, non-Gausianity (skewness and heavy tails), volatility clustering and leverage effect. This makes the modelling task very challenging, see hsieh1991; bollerslev1994; brooks1996; cont2001.
The data are daily returns of $ m =$ 418 equities on the S&P 500 index from 03.01.2006 through 28.09.2018, for a total of 3208 observations. The initial sample 03.01.2006 - 28.09.2017 is used for estimation (training), with $T =$ 2957 in-sample size. The holdout sample period 01.10.2017 - 28.09.2018 (251 observations) is employed to examine the models' out-of-sample forecasting performance. 1-step (here one day) ahead forecasts of targets ($ \hat y_{it+1|t} $) are based on a rolling estimation window. Parameter estimates are updated every five steps.
We believe accounting for comovements between financial returns is important in forecasting returns. Consequently, the lags of other equities are included as predictors for any return series. We examine the nonlinear high-dimensional forecasting model described in the prior sections (AAShNet model) as well as several competing models and benchmarks.
We compare our proposed model with a benchmark, the sample mean of ${\bm y}_t$ over the in-sample window, as the 1-step ahead forecast. This corresponds to assuming the log daily price of follows a random walk (RW) with drift. It is almost equivalent to the \enquote{zero forecast} when the in-sample window is large enough. Furthermore, a buy-and-hold (B&H) strategy in the market portfolio (S&P 500 Index) has been considered as another benchmark. To understand whether allowing nonlinearity improves portfolio performance we examine the AAShNet algorithm (with Ridge and Lasso) optimized by cross-validation.
Since predictability of financial returns has major consequences for financial decision making, the model with minimal forecast error is deemed optimal. However, the model with minimum forecast error does not necessarily guarantee profit maximization, the primary objective of financial decision makers. armstrong1992, pesaran1995,pesaran2000, granger2000 and engle2006 argue that a forecast evaluation criterion should be related to decision making and judge predictability of financial returns in terms of portfolio simulation. More specifically, a trading (portfolio) simulation approach assumes that all competing models are applied with stock market virtual investment decisions, and out-of-sample portfolio performances are used to evaluate the predictability of alternative models.
Consequently, this paper examines both statistical and portfolio performance measures (the out-of-sample RMSE and the portfolio performance during the out-of-sample period). Figure (ref) illustrates portfolio excess returns for the out-of-sample period for the proposed model (ASShNet) against competing approaches. We randomly selected 50 stocks out of 418 stocks to construct the portfolio. However, the forecast of each selected stock is based on the lags of all 418 equities.
Consider a passive, equally weighted $(1/M)$ portfolios with short selling. This portfolio is known to be a very stringent benchmark that many optimization models fail to outperform (see demiguel2009. We compute the portfolio’s out-of-sample excess returns and volatility as well as the Sharpe ratio. Sharpe ratio measures risk-adjusted returns, a portfolio with a greater Sharpe ratio offers greater returns for the same risk. If a portfolio with lower Sharpe ratio has returned better over a time period than another portfolio with a higher ratio, the risk of losing by investing in the former fund will be higher.
The proposed penalized neural net behaves noticeably better in this empirical analysis. Table (ref) provides evidence for out-of-sample forecasting ability of this model vis-\`a-vis competing approaches in terms of the Sharp ratio. AAShNet also offers an appreciable improvement over linear shrinkage models and benchmarks based on RMSE and actual portfolio performance. In the Ridge and Lasso regressions, the best model is selected by cross-validation. We perform generalized cross-validation, which is an efficient leave-one-out cross-validation.
AAShNet produces higher returns (10.24%) at the end of the out-of-sample period, with a Sharpe ratio of (1.143) that is superior to alternative models. This indicates that significantly improved forecast is obtained by modelling nonlinear dynamics among variables. One should note that Random Walk with drift and AR(1) are special cases of shrinkage models and AAShNet when there is no dependence on other equities.
Forecasting with many predictors has received a good deal of attention in recent years. Shrinkage methods are one of the most common approaches for forecasting with many predictors. Such methods have generally ignored nonlinear dynamic relations among predictors and the target variable.
In this study, we suggested an Adaptable and Automated Shrinkage Estimation of Neural Networks (AAShNet). We explained how skip-layer connections move the model in the right direction when the data contains both linear and nonlinear components. To overcome the curse of dimensionality and to manage model complexity, we penalized the model loss function with $L_1$ and $L_2$ norms. Setting the size of shrinkage is still an open question. Recent studies have proposed automated approaches for estimation of algorithm hyperparameters. We employed the gradient-based automated approaches which treat shrinkage hyperparameters in the same manner as the network weights during training, and simultaneously optimize both sets of parameters.
The empirical application to forecasting daily returns of equities in the S&P 500 index from 2006 to 2018 provides support for the out-of-sample forecasting ability of AAShNet algorithm vis-\`a-vis some competing approaches, both in terms of statistical criteria and trading simulation performance. Our empirical results encourage further research toward other possible applications of the proposed model.