Forecasting in Big Data Environments: an Adaptable and Automated Shrinkage Estimation of Neural Networks (AAShNet)

TL;DR

This paper introduces an adaptable, automated shrinkage estimation method for neural networks tailored for high-dimensional, nonlinear forecasting tasks, improving prediction accuracy and interpretability in big data environments.

Contribution

It develops a novel high-dimensional learning approach combining L1 and L2 penalties with gradient-based hyperparameter tuning for neural network forecasting.

Findings

Improved forecast accuracy over existing models.

Effective selection of predictors and network topology.

Pruning of nonlinear components when data is linear.

Abstract

This paper considers improved forecasting in possibly nonlinear dynamic settings, with high-dimension predictors ("big data" environments). To overcome the curse of dimensionality and manage data and model complexity, we examine shrinkage estimation of a back-propagation algorithm of a deep neural net with skip-layer connections. We expressly include both linear and nonlinear components. This is a high-dimensional learning approach including both sparsity L1 and smoothness L2 penalties, allowing high-dimensionality and nonlinearity to be accommodated in one step. This approach selects significant predictors as well as the topology of the neural network. We estimate optimal values of shrinkage hyperparameters by incorporating a gradient-based optimization technique resulting in robust predictions with improved reproducibility. The latter has been an issue in some approaches. This is statistically interpretable and unravels some network structure, commonly left to a black box. An additional advantage is that the nonlinear part tends to get pruned if the underlying process is linear. In an application to forecasting equity returns, the proposed approach captures nonlinear dynamics between equities to enhance forecast performance. It offers an appreciable improvement over current univariate and multivariate models by RMSE and actual portfolio performance.

Figures

Tables

Table 1. Table I:

Portfolio	Return	Sharp Ratio	Ave(RMSE)
AAShNet	10.2%	1.143	0.01558
B&H	9.55%	0.767	-
RW	8.17%	0.770	0.01564
Ridge	5.45%	0.561	0.01613
Lasso	6.87%	0.741	0.01590

Equations

$$y_t={\mathcal{L}}_t+{\mathcal{N}}_t$$

$$y_{t}=\Phi({\bm{x}};{\bm{w}})=\sum_{i\rightarrow k}x_{it}w_{ik}+\sum_{j\rightarrow k}\phi_{j}\bigg{(}\sum_{i\rightarrow j}x_{it}w_{ij}\bigg{)}w_{jk}+\varepsilon_{t},$$

$$w^{new}=w^{old}-\eta \text{leavevmode}\frac{\partial E(\mathbf{w})}{\partial w}$$

$${\mathbf{w}}^{\ast }$$

$$\left\{\begin{array}[]{lr}\Delta=\frac{\partial E({\bm{w}})}{\partial w}+\frac{\partial\Omega({\bm{w},\lambda})}{\partial w}\\ w^{new}=w^{old}\leavevmode\nobreak\ -\leavevmode\nobreak\ \eta\Delta\end{array}\right.$$

$${\mathbf{w}}^{\ast }=\underset{\mathbf{w}}{\mathrm{argmin}}\text{leavevmode}E(\mathbf{w}|\lambda ,X)\text{leavevmode}+\text{leavevmode}\frac{{\lambda }_2}{2}\sum _{i\to k}{w}_{ik}^2\text{leavevmode}+\text{leavevmode}{\lambda }_1(\sum _{i\to j}|w_{ij}|+\sum _{j\to k}|w_{jk}|)$$

$$\left\{\begin{array}[]{lr}w_{ik}^{new}=w_{ik}^{old}-\eta\leavevmode\nobreak\ (\frac{\partial E({w|\lambda,X})}{\partial w_{ik}}+\lambda_{2}w_{ik}^{old})\\ w_{ij}^{new}=w_{ij}^{old}-\eta\leavevmode\nobreak\ (\frac{\partial E({w|\lambda,X})}{\partial w_{ij}}+\lambda_{1}sgn(w_{ij}^{old}))\\ w_{jk}^{new}=w_{jk}^{old}-\eta\leavevmode\nobreak\ (\frac{\partial E({w|\lambda,X})}{\partial w_{jk}}+\lambda_{1}sgn(w_{jk}^{old}))\\ \end{array}\right.$$

$${\lambda }^{\ast }$$

$$\text{s.t.}w(\lambda )$$

$$y_{t}=\Phi({\bm{x}};{\bm{w}})=\alpha\sum_{i\rightarrow k}x_{it}w_{ik}+(1-\alpha)\sum_{j\rightarrow k}\phi_{j}\bigg{(}\sum_{i\rightarrow j}x_{it}w_{ij}\bigg{)}w_{jk}+\varepsilon_{t},$$

$$v_1$$

$$v_2$$

$$v_3$$

$$v_4$$

$$v_5$$

$$y$$

$$\dot{v_1}$$

$$\dot{v_2}$$

$$\dot{v_3}$$

$$\dot{v_4}$$

$$\dot{v_5}$$

$$\dot{y}$$

$$\overline{v_5}$$

$$\overline{v_4}$$

$$\overline{v_3}$$

$$\overline{v_1}$$

$$\overline{v_2}$$

$$\overline{v_1}$$

$$\overline{x_1}$$

$$\overline{x_2}$$

Full text

\usdate

1.5em

Forecasting in Big Data Environments: an Adaptable and Automated Shrinkage Estimation of Neural Networks (AAShNet)

*Ali Habibnia * Department of Economics, Virginia Tech. Email: habibnia@vt.edu   

Esfandiar Maasoumi Department of Economics, Emory University.

ABSTRACT

This paper considers improved forecasting in possibly nonlinear dynamic settings, with high-dimension predictors (“big data” environments). To overcome the curse of dimensionality and manage data and model complexity, we examine shrinkage estimation of a back-propagation algorithm of a neural net with skip-layer connections. We expressly include both linear and nonlinear components. This is a high-dimensional learning approach including both sparsity $L_{1}$ and smoothness $L_{2}$ penalties, allowing high-dimensionality and nonlinearity to be accommodated in one step. This approach selects significant predictors as well as the topology of the neural network. We estimate optimal values of shrinkage hyperparameters by incorporating a gradient-based optimization technique resulting in robust predictions with improved reproducibility. The latter has been an issue in some approaches. This is statistically interpretable and unravels some network structure, commonly left to a black box. An additional advantage is that the nonlinear part tends to get pruned if the underlying process is linear. In an application to forecasting equity returns, the proposed approach captures nonlinear dynamics between equities to enhance forecast performance. It offers an appreciable improvement over current univariate and multivariate models by RMSE and actual portfolio performance.

Key Words: Nonlinear Shrinkage Estimation, Gradient-based Hyperparameter Optimization, High-dimensional Nonlinear Time Series, Neural Networks

JEL classification: C45, C51, C52, C53, C61.

I Introduction

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 (Huber, (2011);Cherkassky et al., (1994); Moody, (1991)). 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 - Pearson, (1901); Eckart and Young, (1936); factor analysis - Spearman, (1904); canonical correlations analysis - Hotelling, (1936), and several others111Cunningham and Ghahramani, (2015) surveyed the literature on linear dimensionality reduction in their work..). Feature selection is accomplished by such methods as Ridge - Hoerl and Kennard, (1970); LASSO - Tibshirani, (1996) and Elastic Net - Zou and Hastie, (2005)).

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 Ng, (2004). 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 Ng, (2004). 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})$.

[TABLE]

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.

II The Model

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 function222Using 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

[TABLE]

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 II.1. 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.(II.1) 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 Bryson et al., (1979) and popularized in the field of artificial neural network research by Werbos, (1988) and Rumelhart et al., (1986). 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 Baydin et al., (2017). 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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

which can be realized by iteratively adjusting the parameters using the updating rules below

[TABLE]

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}$. Lee et al., (2006)

III Gradient-based Hyperparameter Optimization

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 generalization333Generalization 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 (Bergstra et al., 2011b ). Transparent tuning of hyperparameters can be part of an Hyperparameter Optimization (HPO), as an outer loop in training procedures.

The de-facto naïve 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. Bergstra and Bengio, (2012) show empirically and theoretically that random search more efficient than grid search. Statistical techniques such as cross-validation Wahba, (1990), bootstrapping Efron and Tibshirani, (1994), and Bayesian methods MacKay, (1992) 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.444This assumption can be relaxed. see: Chu and Marron, (1991). For a general description of the CV see Geisser, (1975), and Arlot and Celisse, (2010) for a comprehensive review on cross-validation procedures and their applications in different algorithms and frameworks. Several studies such as Rivals and Personnaz, (1999) 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 (Hutter et al., (2011); Bergstra et al., 2011a ; Bergstra et al., (2013); Snoek et al., (2012)), 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 Larsen et al., (1996) and Andersen et al., (1997), emerged in the 1990s. We can distinguish two main approaches of gradient-based optimization: Implicit differentiation and iterative differentiation.

Implicit differentiation, first proposed by Larsen et al., (1996), 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$. Bengio, (2000) 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 Keerthi et al., (2007). Pedregosa, (2016) proposes HOAG which uses inexact gradients, allowing the gradient with respect to hyperparameters to be computed approximately.

In iterative differentiation, first proposed by Domke, (2012), 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ïve 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. Maclaurin et al., (2015) 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:

[TABLE]

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 Goodfellow et al., (2016). We modify our training (Algorithm 1) 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 2 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:

[TABLE]

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.

IV Case Study: Return Prediction

Research into modelling and forecasting financial returns has a long history. Several models are described in Tsay, (2005) and Campbell et al., (1996) 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 Sharpe, (1964) and Lintner, (1965) 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 Fama and French, (2004). This failure is perhaps due to the absence of multiple factors. Arbitrage pricing theory (APT) is a general model proposed by Ross, (1976) 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 Hsieh, (1991); Bollerslev et al., (1994); Brooks, (1996); Cont, (2001).

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 “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. Armstrong and Collopy, (1992), Pesaran and Timmermann, (1995, 2000), Granger and Pesaran, (2000) and Engle and Colacito, (2006) 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 IV.1 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 DeMiguel et al., (2009). 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 I provides evidence for out-of-sample forecasting ability of this model vis-à-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.

V Concluding Remarks

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-à-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.

VI Appendix: Automatic Differentiation

There are three main approaches that computer can work out the derivatives: Numerical, Symbolic and Automatic differentiation. Automatic differentiation refers to a family of procedures to automatically calculate exact derivatives of any function, including program subroutines, with time complexity at most a small constant factor of the time complexity of the original function. It is not inherently ill-conditioned and unstable similar to the numerical method and has much less computational complexity. It also does not suffer from expression swell problem of symbolic differentiation.

AD augments the standard computation with calculation of derivatives whose combination through chain rule gives the derivative for overall composition. AD can be applied on evaluation trace of arbitrary program subroutines which can be more than closed-form functions and are in fact capable of incorporating complex control flows which do not directly alter values. An automatic differentiator takes a code subroutine that computes a function of several independent variables as input and gives as output a code that computes the original function along the gradient of the function with respect to the independent variables. As most of the functions are piece-wise differentiable and control flows not directly interfering with calculations, chain rule can be used repeatedly in such a way that gradients are calculated along intermediate values being computed.

Based on modus operandi of automatic differentiation there can be two implementations of this technique; the forward mode and the reverse mode. We investigate each method, by applying them on the same trivial function $y=f(x_{1},x_{2})=x_{1}x_{2}-cos(x_{1})$ at $(x_{1},x_{2})=(6,3)$.

[TABLE]

In forward mode, we build a Forward Primal Trace of the values propagating through the function and a corresponding Forward Tangent Trace. Eq. VI.1 shows the forward evaluation of primals. Forward primal trace depicts the natural flow of composition. Eq. VI.2 is the corresponding tangent trace for $\dot{y}=\frac{\partial f}{\partial x_{1}}$, that is the rate of change of the function $f$ with respect to the input $x_{1}$. Notice that both traces are evaluated as written, top to bottom. To calculate the derivative with respect to $n$ different parameters, $n$ forward mode differentiations would be needed. This makes the forward-mode very inefficient for deep learning models where the number of parameters may amount to millions.

[TABLE]

The reverse mode works by complementing each intermediate variable $v_{i}$ with an adjoint $\bar{v_{i}}$ representing the sensitivity of output $y$ to changes in $v_{i}$. In reverse mode the code is executed and the trace is stored in memory at first stage. At second stage, the adjoints are calculated in opposite direction of the execution of the original function. The reverse adjoint trace corresponding to Eq. VI.1 is depicted in VI.3.

[TABLE]

For a function $f\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ whose number of operations to be evaluated is denoted by $\textup{ops}(f)$, the complexity of calculating the Jacobian by forward and reverse modes are $n\times c\times\textup{ops}(f)$ and $m\times c\times\textup{ops}(f)$, respectively, where it is guaranteed that $c<6$Griewank and Walther, (2008). That is if $n\gg m$, backward-mode is preferable, although it would have increased memory requirements. And forward mode should be used when the number of dependent variables is greater than the number of independent variables.