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.
83,308 characters · 13 sections · 86 citation commands
Deep Neural Network Estimation in Panel Data Models
Panel data models are widely used in economics and finance. They combine both cross-sectional and time series data. One important advantage of panel data over time-series methods (see, for example, Chapters 26 and 28 of pesaran2015time) is their ability to control for unobserved heterogeneity both in the temporal and longitudinal dimensions. One can then approximate this latent individual heterogeneity through identifiable effects that are otherwise non-detectable in traditional time-series data sets. There are several ways to model and control for individual heterogeneity in linear panel data models: the random effects estimator, see, for example, balestra1966pooling, the fixed effects (within) estimator, see, for example, mundlak1961empirical, mundlak1978pooling, and the swamy1970efficient estimator. Alternative ways to model individual heterogeneity in linear models are found in hsiao1974statistical, hsiao1975some, with a thorough discussion in hsiao2004random, hsiao2008random and Part VI of pesaran2015time.
The work summarized above focuses on linear heterogeneous panel data models. However, the importance of nonlinearity has attracted increased interest in the literature. Notable contributions are fernandez2016individual, who adapt the analytical and jackknife bias correction methods introduced in hahn2004jackknife to nonlinear models with additive or interactive individual and time effects, and CHEN2021296 who address estimation and inference in general nonlinear models using iterative estimation. haciouglu2021common provide an approach for estimation and inference in nonlinear conditional mean panel data models in the presence of cross-sectional dependence. jochmans2017two develops the asymptotic properties of GMM estimators for models with two-way multiplicative fixed effects, while charbonneau2013multiple considers a logit conditional maximum likelihood approach to investigate whether existing panel methods for eliminating a single fixed effect can be modified to eliminate multiple fixed effects.
In this paper we also focus on the estimation of nonlinear panels. We propose the use of a novel machine learning (ML) panel data estimator based on neural networks. To help delineate the contributions of this paper and the empirical application that we consider, we first provide a high-level summary of the current literature on ML.
Statistical ML is a major interdisciplinary research area. In the last decade, ML methods have been incorporated, in various forms, across the natural, social, medical, and economic sciences, leading to significant research outputs. There are two main reasons for such widespread adoption. Firstly, ML methods and specifically neural networks, the focus of this paper, have been found to exhibit outstanding empirical performance when forecasting, specifically with high-dimensional data sets. Secondly, they have great capacity to uncover potentially unknown and both highly complicated and nonlinear relationships in the data. In conjunction with increased availability of high-dimensional data sets, and policymakers' understandable desire for accurate forecasts, considerable attention has been paid to ML.
Studies have shown that feed-forward neural networks can approximate any continuous function of several real variables arbitrarily well; see, for example, hornik1991approximation, hornik1989multilayer, galant1992learning, and par1991. Other nonparametric approaches, for example, splines, wavelets, the Fourier basis, as well as simple polynomial approximations, have the universal approximation property, based on the Stone--Weierstrass theorem. However, it has been convincingly argued that neural networks outperform them in prediction (see, for example, bla2010).
More recent work by liang2016deep and YAROTSKY2017103, yarotsky2018optimal considers feed-forward neural networks as approximations for complex functions that accommodate multiple layers, provided sufficiently many hidden neurons and layers are available. Other examples, like bartlett2019nearly, provide the theoretical framework for neural network estimation, while Hieber2019 focuses on the adaptation property of neural networks, showing that they can strictly improve on classical methods. If the unknown target function is a composition of simpler functions, then the composition-based deep net estimator is superior to estimators that do not use compositions. Lastly, recent work of farrell2021deep, building on the work of YAROTSKY2017103 and bartlett2019nearly, studies deep neural networks and considers their use for semi-parametric inference. \newline In this paper we focus on nonlinear panel data models, where the source of nonlinearity lies in the conditional mean. Our contribution to the literature is as follows. We propose a ML estimator of the conditional mean, $E(y_{it}|\boldsymbol{x}_{it})$, based on neural networks and explore the idea of heterogeneity in a nonlinear panel model by allowing the conditional mean to have a panel -- common nonlinear component -- as well as a nonlinear idiosyncratic component. We base our theoretical results mainly on farrell2021deep, expanding their contribution to a panel data framework. We also find evidence of the double descent effect, whereby complex models can perform well without the need for explicit regularization (see hastie2022surprises and kelly2022virtue, as well as Remark (ref) below).
We use the new deep panel data models to forecast the transmission of new COVID-19 cases during the pandemic across a number of countries. We consider the G7 countries. In contrast to theoretical epidemiological models, that may be specified incorrectly, our proposed neural network models are flexible reduced-form models. They let the data determine the path of new infections over time, by modeling this path as dependent on the lagged levels of the number of infections. By comparing the models against a deep (nonlinear) time-series model, that does not aim to exploit cross-country dependencies, we test whether there are benefits when forecasting new COVID-19 cases to pooling data across countries. We find that there clearly are. Importantly, our model also captures the nonlinear features of a pandemic, particularly in its early waves.
Neural networks have great capacity to approximate complicated nonlinear functions and have been found to forecast well. But they are frequently criticized as non-interpretable (of being a “black box”), since they do not offer simple summaries of relationships in the data. Recently, there have been a number of papers that try to make ML output interpretable; see, for example, athey2017state, wager2018estimation, belloni2014inference, joseph2019shapley, chronopoulos2023forecasting, and kap2023.
In this paper, given the many but contrasting (across time and countries) containment or social-distancing policies instigated to moderate the path of the COVID-19 pandemic, we use our model to shed light on the relative effectiveness -- across time and across the G7 countries -- of these policies at lowering the number of new COVID-19 cases. We do so by exploring how the use of partial derivatives, calculated from the output of our proposed neural network, can help examine the effectiveness of policy. We examine the derivatives over time and find that some, but not all, containment policies were effective at lowering new COVID-19 cases. These policies tended to be more effective two to three weeks after the policy change. There is also considerable heterogeneity across countries in the effectiveness of these policies. Policy, as a whole, was somewhat less effective in Italy, and was more effective in Japan in late summer 2022, later than in the other G7 countries.
The remainder of the paper proceeds as follows. In Section (ref) we introduce our main theoretical results: we discuss non-asymptotic bounds for a (potentially heterogeneous) neural network panel estimator based on a quadratic loss function. In Section (ref), we discuss both methodological and implementation aspects of the proposed methodology. We undertake the modeling and forecasting of new COVID-19 cases, and the assessment of the effectiveness of containment policies, in Section (ref). Section (ref) concludes. We relegate to the online appendix additional forecasting results, data summaries, and further discussion of the prediction evaluation tests used.
Let $y_{it}$ be the observation for the $i^{th}$ cross-sectional unit at time $t$ generated by the following panel data model:
where $\{\boldsymbol{x}_{it}\}=\{(x_{t,1},\ldots,x_{t,p})^{\prime}\}$ is a $p$-dimensional vector of regressors, belonging to unit $i$, and $\widetilde {h}_{i}(\cdot)$ are unknown functions that will be approximated with neural networks. Throughout, we abstract from unconditional mean considerations, for simplicity, by assuming $E(y_{it})=0$. This can be achieved by simple unit-by-unit demeaning of the dependent variable. Therefore, the model we entertain is given by:
where $\varepsilon_{it}$ is an error term. Next, we provide a crucial decomposition to justify the use of a panel structure. We assume that $\widetilde{h}_{i}\left( \boldsymbol{x}_{it}\right)$ can be decomposed as follows:
where the function $h(\cdot)$ is the common component of the model, and is our main focus of interest, and $h_{i}\left( \boldsymbol{x}_{it}\right) $ are idiosyncratic components that will also be approximated with neural networks. Assumptions needed for the identification of $h(\cdot)$ will be given below. The main motivation for this decomposition is the familiar linear heterogeneous panel data model, that takes the form:
where $E(\boldsymbol{\eta}_{i}|\boldsymbol{x}_{it},\varepsilon_{it})=0$. Equation (ref) allows coefficients to vary across individual units. We wish to consider and analyze a nonlinear extension of this heterogeneous panel data model.
The next step of our proposal involves approximating $h(\cdot)$ and $h_{i}(\cdot)$ with neural network functional parameterizations, given by $g\left( \boldsymbol{\cdot};\boldsymbol{\theta}\right) $. Here, the functional form is known up to the parameter vector $\boldsymbol{\theta}$, which is a vector of ancillary parameters, such as network weights and biases. More details on the choice of $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$ and the role of various neural network parameters will be provided in Section (ref) below. Therefore, we parameterize (ref) by proposing the following panel model:
where $\boldsymbol{\theta}^{0}$and $\boldsymbol{\theta}_{i}^{0}$ denote the values of the parameters that best approximate $h$ and $h_{i}$, respectively (see (ref)), in a sense to be defined below.
It is useful to draw some parallels between (ref) and (ref). We note that most multi-layer neural network architectures have a final linear layer given by:
where $\boldsymbol{f}$ is a vector of known functions that form part of the neural network architecture and $L$ denotes the number of network layers. Then, it follows that we have a linear representation, in $\boldsymbol{f}$ and $\boldsymbol{f}_{i}$, of the form:
which is reminiscent of (ref) and thus provides a clear rationale for our nonlinear extension of it.
Furthermore, it provides a rationale for thinking that $\boldsymbol{\theta}_{i}^{0}$ plays a similar role to the idiosyncratic coefficients, $\boldsymbol{\eta}_{i}$, of the linear model. Of course, one can use a different network architecture for the panel and idiosyncratic components, but for simplicity we keep the same structure. The model above encompasses a variety of nonlinear specifications. It is also worth emphasizing that the dimension of the regressor vector could be very large. So it is conceivable that each $\boldsymbol{x}_{it}$ contains regressors from other cross-sectional units, allowing for complex nonlinear interactions across units. In the limit, each unit could have $(\boldsymbol{\boldsymbol{x}}_{1t} ,\ldots,\boldsymbol{\boldsymbol{x}}_{Nt})$ as the regressor vector.
Next, we consider the conditions needed to identify $h(\cdot)$. We require certain definitions. First, we define $\boldsymbol{\varepsilon}_{it}\equiv y_{it}-h\left( \boldsymbol{x}_{it}\right) -h_{i}\left( \boldsymbol{x} _{it}\right) $ and $u_{it}=h_{i}\left( \boldsymbol{x}_{it}\right) +\boldsymbol{\varepsilon}_{it}$, where the latter is in analogy to the usual composite error term for the linear heterogeneous panel model, given as $\boldsymbol{x}_{it}^{\prime}\boldsymbol{\eta}_{i}+\varepsilon_{it}$. The assumption below generalizes the usual identification assumption on $\boldsymbol{\eta}_{i}$ made in linear heterogeneous panel models.
This assumption enables separation of $h(\cdot)$ and $h_{i}(\cdot)$ when panel pooled estimation is carried out. For neural network estimation, much stricter assumptions will be needed. In particular, the second part of the assumption is justified in view of our later assumption that $h(\cdot)$ and $h_{i}(\cdot)$ can be well approximated by neural network architectures and the linear aspect of neural networks discussed in (ref) as it is similar, in functionality, to assuming that $E(\boldsymbol{\eta}_{i}|\boldsymbol{x}_{it})=0$.
Next we align our discussion with farrell2021deep. The overall goal of neural network estimation in farrell2021deep is to estimate an unknown smooth function $h(\cdot)$ that maps covariates, $\boldsymbol{X}$, to an outcome ($T\times N$) matrix $\boldsymbol{Y}$, by minimizing a loss function $g_{\ast}\left( \boldsymbol{Y},\boldsymbol{X};\boldsymbol{\theta }\right) $ with respect to the parameterization $\boldsymbol{\theta}$ of a neural network function $g\left( \boldsymbol{\cdot}; \boldsymbol{\theta}^{0}\right) $. Formally,
This is a minimization of a population quantity and assumes that the true function $h(\cdot)$ is the unique solution of (ref). Note that while farrell2021deep do not specify a true function $h(\cdot)$, we take a further step and assume that the true functions in (ref) coincide with the unique solutions of (ref). We do not specify \ $\boldsymbol{Y}$ and $\boldsymbol{X}$ further, since we will apply this general estimation strategy both to get a panel-based estimate of $h(\cdot)$ and estimates of $h_{i}(\cdot)$ via unit-specific estimation.
For now, we present further sufficient general conditions on $h(\cdot)$ and $g_{\ast}\left( \boldsymbol{Y},\boldsymbol{X};\boldsymbol{\theta}\right)$ in order for our results to hold. We require the following assumptions:
In our case, and in what follows, we specialize the general framework above by using a squared error loss function that, for the panel setting, becomes: \[ g_{\ast}\left(\boldsymbol{Y},\boldsymbol{X};\boldsymbol{\theta}\right)=\frac{1}{NT}\sum _{i=1}^{N}\sum_{t=1}^{T}\left( {y}_{it}-g(\boldsymbol{x}_{it} ;\boldsymbol{\theta})\right) ^{2}. \] In our analysis we use feed-forward neural networks architectures with rectified linear unit (ReLU) activation functions and weights that are unbounded following farrell2021deep and the discussion below. Such networks approximate smooth functions well, as shown in YAROTSKY2017103, yarotsky2018optimal.
A further assumption is required on the processes $\{ \boldsymbol{x}_{t, k} \}$ and $\{ {\varepsilon}_{it} \}$, for some $k=1,\ldots, p$, where $p$ is the number of covariates.
Having rewritten the loss function, as squared error loss, with respect to the re-parameterized panel model in (ref), we construct a pooled-type nonlinear estimator, ${\boldsymbol{\widehat{\theta}}}$, such that:
which obeys Assumption (ref). Therefore, our estimator of $h\left( \boldsymbol{x}_{it}\right) $ is given by $g(\boldsymbol{x}_{it} ;\boldsymbol{\widehat{\theta}})$. Then, we proceed to estimate $h_{i}\left( \boldsymbol{x}_{it}\right) $ by $g(\boldsymbol{x}_{it} ;\boldsymbol{\widehat{\theta}_{i}})$, where $\boldsymbol{\widehat{\theta}_{i}}$ is given by:
for each $i$, given $\boldsymbol{\widehat{\theta}}$ from (ref). \
Next, we argue that the estimation in (ref) can effectively separate $h\left( \boldsymbol{x}_{it}\right) $ from $h_{i}\left( \boldsymbol{x} _{it}\right) $ and that the unit-wise second step estimation in (ref) can retrieve $h_{i}\left( \boldsymbol{x}_{it}\right) $. We do this by noting the following. Consider the loss function for $\boldsymbol{\widehat{\theta}}$ in (ref). We have:
Under Assumptions (ref)--(ref), terms $A_{2}$ and $A_{3}$ converge in probability to positive limits, while $A_{5}$ and $A_{6}$ converge, in probability to zero, and in fact, furthermore, both are $O_{p}((NT)^{-1/2})$. Additionally, under (ref), $A_{4}$ is $O_{p}(N^{-1/2} )$. Then it immediately follows that \ the loss function is minimized when $\boldsymbol{\theta=\theta}^{0}$, in view of our identification assumption in (ref). It therefore follows that $\boldsymbol{\widehat{\theta} }\rightarrow^{p}\boldsymbol{\theta}^{0}$ and $g(\boldsymbol{x}_{it} ;\boldsymbol{\widehat{\theta}})\rightarrow^{p}g(\boldsymbol{x}_{it} ;\boldsymbol{\theta}^{0})=h\left( \boldsymbol{x}_{it}\right) $. This proves that the best pooled panel neural network approximation coincides with the true panel function.
Next, we can consider a closely related and, in fact, asymptotically equivalent minimization problem given by:
and the associated model with a composite error is:
where:
Note that $u_{i}$ obeys Assumption (ref).(ref). Moreover, this setting corresponds to that of Theorem 1 in farrell2021deep enabling the use of the rates derived in this theorem. This analysis is summarized and extended in the following proposition:
The proof of Proposition (ref) follows from the proof of Theorem 1 in farrell2021deep, using the arguments made above the Proposition to recast our panel framework into the one of farrell2021deep, by separately identifying $h$ and $h_i$. In Proposition (ref), we use the results from Theorem 1 of farrell2021deep to obtain an asymptotic rate of convergence for the error in (ref). It is clear that this rate of convergence is not optimal, since $\psi<1/2$. We have provided a simplified result compared to Theorem 1 of farrell2021deep. Refinements related to factors, such as the depth and width of the neural network used, can be obtained. These are also discussed in Theorem 6 of bartlett2019nearly and in Lemma 6 of farrell2021deep. Fast convergence of (ref) depends on the trade-off between the number of neurons and layers, and more specifically on the parameterization of their relationship, that controls the approximating power of the network.
We note that, in addition, one can obtain consistency for $\boldsymbol{\widehat{\theta}} _{i}$ by minimizing the loss over $\boldsymbol{\theta}_{i}$: $L_{i}=\frac{1} {T}\sum_{t=1}^{T}[y_{it}-g(\boldsymbol{x}_{it};\boldsymbol{\widehat{\theta} })-g(\boldsymbol{x}_{it};\boldsymbol{\theta}_{i})]^{2}.$ Given the rate in Proposition (ref), it immediately follows that $g\left( \boldsymbol{x}_{it};\boldsymbol{\theta}_{i}^{0}\right) $ can be consistently estimated at rate $T^{-\psi}$, as long as $T=o(N^{\xi})$ for some $\xi<1$, given that then the uniform rate in Proposition 1 is faster than $T^{-\psi}$.
In this section we provide details on implementation of the proposed nonlinear estimators. First, we summarize the overall neural network construction which relates to the choice of the neural network's architecture and can be summarized by the functional parameterization $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$, used in the approximation of $h(\cdot)$. We limit our attention to the construction of $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$, since it directly applies to $g\left(\boldsymbol{\cdot};\boldsymbol{\theta_i}\right)$. Then we illustrate how regularization can be applied in the context of the proposed estimators. Finally, we discuss both the cross-validation exercise used to select the different parameters and hyperparameters of the corresponding network and optimization algorithm.
We focus on the construction of the feed-forward neural network functional parameterization, $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$, used to approximate $h(\cdot)$ in Section (ref). The feed-forward architecture consists of: an input layer, where the covariates are introduced given an initial set of weights to the inner (hidden) part of the network; the hidden layers, where a number of computational nodes are collected in each hidden layer and nonlinear transformations on the (weighted) covariates occur; and the output layer that gives the final predictions and a choice for the activation function $\sigma(x):\mathbb{R}\rightarrow\mathbb{R}$ that is applied element-wise. The architecture is feed-forward, since in each of the hidden layers there exist several interconnected neurons that allow information to flow from one layer to the other, but only in one direction. The connections between layers correspond to weights.
We use $L$ to define the total number of hidden layers and $M^{(l)}$, $l=1, \ldots, L$ to define the total number of neurons at the $l^{th}$ layer. $L$ and $M^{(l)}$ are measures for the depth and width of the neural network, respectively. We use the ReLU activation function, $\sigma_l(\boldsymbol{X}_{t}):=\max(\boldsymbol{X}_{t}, 0)$, where $\boldsymbol{X}_{t}$ is a $N\times p $ matrix of characteristics for $ t=1,\ldots, T$; $l=1, \ldots,L-1$ and a linear activation function for $l=L$. The activation functions are applied elementwise. To explain the exact computation of the outcome of the feed-forward neural network, we focus on the pooled-type estimator in (ref). We assume that the widths (the number of neurons), $M^{(l)}$, and depth (the number of hidden layers), $L$, of the network are constant positive numbers.
Each of the neurons undergoes a computation similar to the linear combination received in each hidden layer $l$: $\boldsymbol{g}^{(l)} = \sigma_{l} (\boldsymbol{g}^{(l-1)}\boldsymbol{W}^{(l)'} + \boldsymbol{b}^{(l)'}) $, while the final output of the network is $\boldsymbol{g}^{(L)} =\boldsymbol{g}^{(L-1)}\boldsymbol{W}^{(L)'} + \boldsymbol{b}^{(L)'}$ and $\boldsymbol{g}^{(0)} = \boldsymbol{X}_{t}$. We can then define for some $t=1,\ldots, T,$ $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$ as:
where $\boldsymbol{W}^{(l)} $ is a $M^{(l)} \times M^{(l-1)}$ matrix of weights, $\boldsymbol{b}^{(l)}$ is a $M^{(l)} \times N$ matrix of biases at layer $l$, with $\boldsymbol{b}^{(1)} = \boldsymbol{0}$. Notice that at $l=1$, the dimensions of $\boldsymbol{W}^{(1)} $ are $M^{(1)}\times p$ and of $\boldsymbol{b}^{(1)}$ are $M^{(1)} \times N$. At the final layer, that is, at $l=L$, the dimensions of $\boldsymbol{W}^{(L)} $ are $1 \times M^{(L-1)}$ and of $\boldsymbol{b}^{(L)}$ are $1 \times N$.
Note that throughout the paper we use $\boldsymbol{\theta}$ to denote a stacked vector containing all ancillary trainable parameters affiliated with the network estimation, as defined below:
We define the overall number of parameters as $d=|\boldsymbol{\theta}|$. The optimization of the neural network proceeds in a forward fashion (from the input layer, that is, $l=1$, to the output $l=L$) and layer-by-layer through an optimizer, for example, a version of stochastic gradient descent (SGD), where the gradients of the parameters $(\boldsymbol{W}^{(l)}, \boldsymbol{b}^{(l)})$ are calculated through back-propagation (using the chain-rule) to train the network.
The specific choice of the network architecture is crucial and affects the complexity and the approximating power of $g\left(\boldsymbol{\cdot};\boldsymbol{\theta}\right)$ in (ref). Our analysis involves primarily theoretical arguments that are widely applicable in feed-forward neural networks when we deal with panel data. We present an example of a feed-forward neural network, based on (ref), in Figure (ref).
The neural network in Figure (ref) consists of two inputs $\boldsymbol{X}_{t}\in\mathbb{R}^{N\times p}$, $\boldsymbol{X}_{t} = ( \boldsymbol{x}_{t}^{(1)} , \boldsymbol{x}_{t}^{(2)} )$, in particular, $p=2$, where $\boldsymbol{x}_{t}^{(j)}$ is a $N\times 1$ vector of one characteristic at $t=1,\ldots, T$ for some $j=1,2$, and one fitted output $\boldsymbol{\widehat{y}}_{t}$. Between the inputs and output $(\boldsymbol{X}_{t}, \boldsymbol{\widehat{y}}_{t})^{\prime}$, are $M$ hidden computational nodes/neurons, in particular $M=5$. The neurons are connected directly forming an acyclic graph which specifies a fixed architecture.\footnote{ The network in Figure (ref) can be used to optimize (ref) and also (ref), but now used for each $i= 1,\ldots,N$. }
Notice that the illustration in Figure (ref) can correspond to a nonlinear pooled-type estimation of $g( \boldsymbol{X} _{t};\boldsymbol{\theta}^{0}) $ in (ref), where we use (ref) to obtain $\boldsymbol{\widehat{\theta}}$, defined in (ref), with input $\boldsymbol{X}_{t} = (\boldsymbol{x}^{(1)} _{t}, \ldots, \boldsymbol{x}^{(p)}_{t} )$, and output $\boldsymbol{\widehat{y}}_{t}\in\mathbb{R}^{N\times 1}. $ The remainder of (ref), $g( \boldsymbol{X}_{t};\boldsymbol{\theta}_{i}^{0}), \; i=1,\ldots,N, \; t=1,\ldots, T $ can be described conceptually as the heterogeneous component, that differs cross-sectionally. One can obtain the estimate of this heterogeneous component following the same steps as those used to obtain $g( \boldsymbol{x}_{it};\boldsymbol{\widehat{\theta}}) $, with the only difference that now the feed-forward neural network is estimated unit-wise, similar to the logic of a fixed effects estimator for linear panel data models.
In this section we discuss some operational implementation aspects required for the estimation of the panel neural network estimators proposed in Section (ref). We focus discussion on the following panel estimator, ${g}(\boldsymbol{x}_{it}; \boldsymbol{\widehat{ \theta}})$, obtained from the optimization of (ref): \[ g(\boldsymbol{x}_{it};\boldsymbol{\widehat{\theta}})= \underset{_{\boldsymbol{\theta}\in\mathbb{R}^{d}}}{\arg\min} \, \frac{1}{N}\sum_{i=1}^{N}\left[ \frac{1}{T}\sum_{t=1}^{T}y_{it}-\frac{1 }{T}\sum_{t=1}^{T}g(\boldsymbol{x}_{it};\boldsymbol{\theta })\right] ^{2}. \]
This nonlinear panel estimator, and generally neural network estimators, have many significant advantages over traditional panel models, mainly summarized in their great capacity at approximating highly nonlinear and complicated associations between variables and outstanding forecasting performance; see, for example, the discussion in goodfellow2016deep and gu2020empirical, gu2020autoencoder. In order to be able to minimize (ref) and obtain a feasible solution for the panel estimator ${g}(\boldsymbol{x}_{it}; \boldsymbol{\widehat{ \theta}})$, we need to choose the overall architecture of the neural network. Following the discussion above, this reduces to choices for the total number of layers $L$, total number of neurons $M^{(l)}$, at each $l=1,\ldots, L$ layers, a loss function $g_{\ast}(\boldsymbol{y},\boldsymbol{X;\theta})$, which in this paper is taken to be the MSE loss, an updating rule for the weights (learning rate, $\gamma$) during optimization, and the optimization algorithm itself, typically taken to be some variant of SGD.
However, neural networks tend to overfit, which can lead to a severe deterioration in their (forecasting) performance. A common empirical solution to this is to impose a penalty on the trainable parameters of the neural network, $\boldsymbol{\theta}$. The penalized estimator based on the LASSO is obtained as the solution to the following minimization problem: \[ g(\boldsymbol{x}_{it};\boldsymbol{\widehat{\theta}})^{\textsc{LASSO}}= \underset{_{\boldsymbol{\theta}\in\mathbb{R}^{d}}}{\arg\min} \, \frac{1}{N}\sum_{i=1}^{N}\left[ \frac{1}{T}\sum_{t=1}^{T}y_{it}-\frac{1 }{T}\sum_{t=1}^{T}g(\boldsymbol{x}_{it};\boldsymbol{\theta })\right] ^{2} + \lambda \left \Vert \boldsymbol{\theta} \right \Vert_{1}, \] where $\lambda$ is the regularization parameter. Note that while explicit regularization improves empirical solutions of neural networks estimators under low signal-to-noise ratios, its role is not clear theoretically, since there are cases where simpler SGD solutions present similar solutions; see, for example, zhang2021understanding. Other commonly used regularization techniques frequently employed empirically to assist in the estimation of neural networks, relate to batch normalization, early stopping, and dropout. We succinctly discuss batch normalization below, given its importance because of the cross-sectional aspect of our estimator. We refer the reader to gu2020empirical for a detailed discussion of early stopping and dropout.
Batch normalization, proposed by ioffe2015batch, is a technique used to control the variability of the covariates across different regions of the network and datasets. It is used to address the issue of internal covariate shift, where inputs of hidden layers may follow different distributions than their counterparts in the validation sample. This is a prevalent issue when fitting, in particular, deep neural networks. Effectively, batch normalization cross-sectionally demeans and standardizes the variance of the batch inputs.
The cross-validation (CV) scheme consists of choices on the overall architecture of the neural network: the total number of layers ($L$), neurons ($M$), the learning rate ($\gamma$) of SGD, the batch size, dropout rate, level of regularization $(\lambda)$, and a choice on the activation functions.
Regarding the choice on the activation functions, we use ReLU for the hidden layers and a linear function for the output layer. We tune the learning rate of the optimizer, $\gamma$, from five discrete values in the interval $\left[0.01, 0.001\right]$. We tune the depth and width of the neural networks using the following grids, $[1,3,5,10,15]$ and [5, 10, 15, 20, 30], respectively. Hence the choice between deep or shallow learning is completely data-driven, as it is selected from the CV scheme. We set the batch size to 14. For the tuning of the regularization parameter, $\lambda$, used for LASSO penalisation, we use the following grid $c\sqrt{\log{p}/NT}$, where $c=[0.001, 0.01, 0.1, 0.5, 1, 5, 10]$. We also use dropout regularization, where the dropout probability is up to $10$ percent; see, for example, gu2020empirical.
To select the trainable parameters, $\boldsymbol{\theta}$, and the hyperparameters discussed above, we follow gu2020empirical, gu2020autoencoder and divide our data into three disjoint time periods that maintain the temporal ordering of the data: the training sub-sample, which is used to estimate the parameters of the model, $\boldsymbol{\theta}$, given a specific set of hyperparameters; the validation sub-sample, which is used to tune the different hyperparameters given $\boldsymbol{\widehat{\theta}}$ from the training sub-sample\footnote{Note that while $\boldsymbol{\widehat{\theta}}$ is used in the tuning of the hyperparameters, it is only estimated at the training sub-sample.}; and, finally, the testing sub-sample which is truly out-of-sample and is used to evaluate our nonlinear models' forecasting performance. As discussed in detail below, our forecasting exercise is recursive, based on an expanding window size. Hence, at each expanding window, we need to use the train-validation-split of the sample and estimate the relevant parameters and tune the hyperparameters. At each expanding window, let $T^\ast$ denote the total sample size for the specific window, then the training sub-sample consists of $\lfloor 0.8T^\ast \rfloor$, the validation sub-sample consists of $\lfloor 0.2 T^\ast \rfloor - c$, and finally, the testing sub-sample consists of 7, 14, or 21 observations depending on the forecast horizon, $h$, respectively. $c$ is chosen so that the testing sub-sample always has $h$ observations and $\lfloor \cdot \rfloor$ stands for the floor function.
The estimation of neural networks is generally a computational cumbersome optimization problem due to nonlinearities and non-convexities. The most commonly used solution utilizes SGD to train a neural network. SGD uses a batch of a specific size, that is, a small subset of the data at each \textcolor{black}{epoch (iteration)} of the optimization to evaluate the gradient, to alleviate the computation hurdle. The step of the derivative at each epoch is controlled by the learning rate, $\gamma$. We use the adaptive moment estimation algorithm (ADAM) proposed by kingma2014adam\footnote{ADAM is using estimates for the first and second moments of the gradient to calculate the learning rate.}, which is a more efficient version of SGD. Finally, we set the number of epochs to $5,000$ and use early stopping following gu2020empirical to mitigate potential overfitting.
In this section, after introducing the data, we examine the predictive ability of the proposed model(s) for forecasting the daily path of new COVID-19 cases across the G7 countries. We compare the forecasting results from our new models against two restricted alternatives: a neural network without a cross-sectional dimension and a linear panel data VAR (PVAR). Comparison against these alternatives lets us examine the importance of firstly modeling the panel dimension and secondly of allowing for nonlinearities. To assess the out-of-sample Granger causality of pandemic-induced lockdown policies on the spread of COVID-19, we compare the forecasting performance of our models with and without measures of the stringency of government-imposed containment and lockdown policies. Such (non-pharmaceutical) policies were differentially adopted by many countries from March 2020, including the G7, to reduce the spread of COVID-19. Then, we discuss how partial derivatives can be used to help interpret the output of the deep panel models. They can be used to help assess the efficacy of the different containment policy measures taken by individual countries to contain the spread of COVID-19.
Our interest is modeling and forecasting, at a daily frequency, reports of new COVID-19 cases per 100K of the population over the sample April 2020 through December 2022 for the G7 countries. We source these data from the World Health Organization coronavirus dashboard.
As $\boldsymbol{x}_{it}$ variables, for each country, $i$, at day, $t$, we consider a set of 7 lagged COVID-19 related indicators, as well as lags of new cases per-100K (our $y_{it}$ variable). For parsimony, we confine attention to lags at 7, 14, 21, and 28 days. These 7 variables, plus lags of the dependent variable, may all have explanatory power for $y_{it}$. The 7 variables (all reported per 100K of the population) comprise: new deaths, the reproduction rate, new tests, the share of COVID-19 tests that are positive measured as a rolling 7-day average (this is the inverse of tests per case), the number of people vaccinated, the number of people fully vaccinated, and the number of total boosters. knutson2022estimating, mathieu2021global, and CAPORALE202277 also consider such COVID-19 related variables, given that they all likely relate (contemporaneously or at a lag) to the number of new COVID-19 cases.
To assess the role of containment policies in explaining and forecasting the spread of new COVID-19 cases, we then consider specifications that augment the aforementioned set of $\boldsymbol{x}_{it}$ variables by adding in a measure or measures of the stringency of the government response to COVID-19. Specifically, we use the government response stringency index, as compiled by the Oxford Coronavirus Government Response Tracker (OxCGRT). This index is a composite measure based on 9 response indicators, namely: school closures, workplace closures, the cancellation of public events, restriction on gatherings, public transport closing, requirements to stay at home, movement restriction, restrictions on international travel, and public information campaigns. Throughout the pandemic the Oxford stringency index was a widely consulted measure of policy. Since the Oxford index is an aggregation of 9 indicators, with the weights subjectively chosen by Oxford, we also experiment with forecasting when the underlying 9 disaggregates enter individually into our models, so that, in effect, we objectively use the data to weight the disaggregates. Note that we always consider the lagged effects of policy changes on new COVID-19 cases, mitigating endogeneity concerns that, for example, stricter lockdown policies follow increases in new COVID-19 cases.
Throughout, $t$ corresponds to a day and we use a trailing seven-day rolling average to smooth the data. The cross-sectional dimension of our panel is $p=36$ when we consider the aggregate stringency index (as published by Oxford) and $p=68$ when we consider the disaggregated stringency index. We further follow the literature (see, for example, gu2020autoencoder) and rank-normalize all of our variables into the $[0, 1]$ interval as follows: \[ \widetilde{\boldsymbol{x}}_i = \frac{\boldsymbol{x}_i - \min(\boldsymbol{x}_i)}{\max(\boldsymbol{x}_i)-\min(\boldsymbol{x}_i)}, \quad i=1,\ldots, N. \] This normalization minimizes the influence of severely outlying observations stemming from covariate distributions that may have significant departures from normality, a common feature of COVID-19 data, especially at the beginning of the pandemic.
The online Data Appendix provides additional data details. Figure (ref) presents the aggregate stringency index and plots new COVID-19 cases per 100K of the population through our sample period. This figure shows that there are apparent commonalities across countries, both in the stringency of policy and the evolution of new COVID-19 cases. But there are differences too, with Japan standing out as having looser containment policies than the other countries during mid-2020 and then experiencing a later spike in new COVID-19 cases in summer 2022. Thus, it remains an empirical question whether forecasting new COVID-19 cases is improved by pooling information across countries.
We recursively produce forecasts of $y_{it}$ -- new COVID-19 cases -- by estimating our set of models using expanding estimation windows and evaluate these forecasts over the out-of-sample period February 6, 2021 through December 24, 2022. Given ((ref)), the $h-$day-ahead forecast of new COVID-19 cases per 100K is:
where $ \widehat{g}\left( \boldsymbol{x}_{it+h};\boldsymbol{\theta }^{\star}\right)$ denotes the corresponding fit of the pooled network, $\widehat{g}\left( \boldsymbol{x}_{it+h};\boldsymbol{\theta }_{i}^{\star}\right) $ denotes the unit-by-unit fit of the network, and $\widehat{y}_{i,t+h} |\mathcal{F}_t$ denotes the deep idiosyncratic forecast. $ \mathcal{F}_t$ denotes the information set up to time $t$, for some $t=1,\ldots, T$, $ \boldsymbol{\theta }^{\star}$ denotes the optimal weights obtained from the CV for the deep pooled model, and $\boldsymbol{\theta }_{i}^{\star}$ denotes the optimal weights obtained from the CV for the deep idiosyncratic model. We compare this forecast against, what we call, the “deep pooled” forecast that sets $\widehat{g}\left( \boldsymbol{x}_{it};\boldsymbol{\theta }_{i}^{\star}\right) =0$.
We recursively compute $h=7$, $h=14$, and $h=21$ day-ahead forecasts using an expanding estimation window (relating $y_{i,t+h}$ to $\boldsymbol{x_{it}}$, as per ((ref))). To ease the computational burden, given that we re-estimate the model and use CV (as discussed in Section (ref)) at each window, we increase the size of the estimation windows in increments of 7 days. We now summarize how estimation and forecasting works for $h=7$ (forecasting at the longer horizons proceeds analogously): We first estimate our models using daily data from April 1, 2020 through January 30, 2021 ($T^{0}=305$) and produce forecasts 7 days-ahead. Then we estimate from April 1, 2020 through February 6, 2021 ($T^{1}=312$) and again produce forecasts 7 days-ahead. We carry on this process until we finally estimate our models over the sample April 1, 2020 through December 17, 2021 ($T^{700}=991$) producing forecasts 7 days-ahead. This results in an out-of-sample sample size of $700$ days. We do not consider forecasting earlier than 7 days-ahead, given that the incubation period of COVID-19 is typically around one week, so that we should not expect policy changes to have effects within one week. During the first wave of the pandemic, many governments revised their policy measures to restrict the virus once a week, which also helps rationalize our choice of forecast horizons. Forecasts for longer horizons, $h$, are obtained similarly.
To test if and how our proposed deep neural network panel data models confer forecasting gains, we compare them against two benchmarks that switch off firstly panel (cross-country) interactions and secondly nonlinear effects. We do so by estimating: (i) a “deep time-series” model that is identical to our deep neural network panel data model but is estimated separately for each country; and (ii) a panel VAR (PVAR) model that does allow for cross-country interactions, but assumes linearity in terms of how $x_{it}$ affects $y_{it+h}$. Testing our model against these two special cases isolates whether it is allowing for cross-country interaction and/or for nonlinearity that is advantageous.
We follow canova2009panel and specify the $i^{th}$ equation of the PVAR with $q$ lags as:
where $A_{ji}$ for $j=1,\ldots, q$ are coefficient matrices, we have dropped the intercept for notational simplicity, $\boldsymbol{Y}_t= ({z}_{1t}', \ldots, {z}_{Nt}')'$, and $z_{it}= ({y}_{it}, \boldsymbol{x}_{it})'$. We set $q= 28$. We estimate the PVAR by OLS and compute $h$-day-ahead forecasts of $y_{it+h}$ from ((ref)) via iteration.
In this section we evaluate the forecasting performance of the proposed nonlinear panel estimator(s) relative to the two benchmark models, namely, the linear PVAR(28), and the deep time-series neural network. We then examine whether the inclusion of policy related variables affects forecast accuracy. Specifically, to test for out-of-sample Granger causality of the policy measures adopted by governments to contain the spread of COVID-19, we compare the forecast accuracy of all of our models with and without the aggregate and disaggregate Oxford stringency indexes.
We evaluate the accuracy of the forecasts of new COVID-19 cases using the root mean squared forecast error (RMSE):
We use the Dieb95 {(DM)} test to test whether differences in forecast accuracy across models are statistically significant. We follow harvey1997testing and use their small-sample adjustment.
Table (ref) compares the accuracy of our two deep nonlinear models against the two benchmarks when we do not include the stringency-based measures of policy and instead focus on predicting new COVID cases using lags of new COVID cases and the other 7 COVID-related measures. The results are striking. Both deep nonlinear panel models provide significant forecasting gains over both the linear PVAR(28) model and the deep time-series neural network at all three forecast horizons. This shows the importance of both the panel dimension and nonlinearites in forecasting the daily path of new COVID-19 cases across the G7 countries. Of the two deep models, the deep pooled estimator delivers, for all 7 countries, more accurate forecasts than the deep idiosyncratic model. Simpler models often work better when forecasting and this appears to be the case here too: allowing for additional country-specific effects in our deep pooled model hinders out-of-sample forecasting performance. Tables (ref)-(ref) in the online appendix show that the forecasting gains, of the deep models against the time series model, are statistically significant. This evidences that the gains from modeling and forecasting new COVID-19 cases come from pooling data (in a nonlinear manner) across the G7 countries.
We next test whether the containment or lockdown policies, imposed at the national-level, help forecast new COVID-19 cases. If the policies were effective, conditioning on them should deliver more accurate forecasts. Table (ref) presents the relative RMSE ratios for each of the four forecasting models when estimating including and excluding the aggregate stringency index. We see that focusing on the deep models, given their higher accuracy as seen in Table (ref), policy as measured by the aggregate stringency index was only effective in France and Japan at 7 days. In the other 5 countries, the RMSE ratios are greater than unity, indicating that better forecasts of new COVID-19 cases are made without the stringency index. Interestingly, for the less accurate deep time-series and PVAR models, policy appears to have been more effective. But consistent with it taking time for policy changes to affect the path of the pandemic, Table (ref) shows that after an additional two weeks policy was effective in all G7 countries, except Italy and the US.
Table (ref) then tests whether the Oxford stringency data have more value-added when forecasting if we let the models decide how much weight to attach to each of the 9 components (policy levers) in the aggregate stringency index. The fact that the RMSE ratios, for the preferred deep models, are now less than unity across all 7 countries indicates that policy was effective after all: but it is important to let the data determine what policies matter in which country. Table (ref) indicates that at $h=7$ days policy was least effective in Canada and Italy, as while policy interventions still affect new COVID-19 cases, unlike in the other G7 countries, these effects are not statistically significant. However, again demonstrating that policy changes take time to have impact, policy has a larger effect after another week (at $h=14$ days), as in both Canada and Italy the relative RMSE ratios are lower at 14 days than at 7 days.
In the online appendix, we provide additional checks on the forecasting performance of our models. We show that the forecasting gains from the models conditioning on the disaggregate stringency index are often stronger in the first half of our out-of-sample window, when in absolute terms the forecasting errors were higher as COVID-19 infection rates were higher and more volatile. Analysis also indicates that the gains of our deep pooled models, over the linear PVAR model, were higher during these earlier waves of COVID-19. This is consistent with the pandemic exhibiting highly nonlinear features in its earlier waves, before vaccinations and other immunities helped restrain the spread of COVID-19. The fluctuation test of giacomini2010forecast is used to show that policy in Italy and Japan proved to be effective later than in the other G7 countries: it is only by the fall of 2022 that we see policy having a marked effect on forecast accuracy. We also present results with the lasso penalization and discuss the observed double descent pattern that our deep models forecast better without any penalty.
A common critique of ML algorithms is their putative trade-off between accuracy and interpretability. The output of a highly complicated ML model, such as a deep neural network of the sort we consider, may fit the data well in-sample and even, as we find, out-of-sample. But the model itself is often hard to interpret. In this section, we illustrate how the use of partial derivatives provides one way to assess the impact of covariates. We focus on examination of the effects of changes in policy, as measured by the aggregate and the disaggregate stringency indexes, on the transmission of new COVID-19 cases.
The use of partial derivatives to interpret model output is, of course, common practice in econometrics, ranging from the simple linear regression model to impulse response analysis. In this section, we show how partial derivatives can be used in deep neural networks to interpret highly nonlinear relationships between covariates and the dependent variable.\footnote{We prefer the use of partial derivatives over Shapley additive explanation values, as proposed by lundberg2017unified, since derivatives tend to be less noisy (see, for example, chronopoulos2023forecasting) and computationally less expensive to compute. Perhaps, though, the biggest disadvantage is the set of implicit assumptions, used in the operational construction of Shapley values. A major one is the assumption that inputs are statistically independent. This is discussed in AIonshap, who also discuss solutions. However, these are computationally intensive, potentially still quite poor approximations, and not appropriate for large sets of inputs. While partial derivatives (as well as coefficients in linear models) have similar issues, as discussed in pes2014a, these issues are both more transparent in nature, and, as discussed in pes2014a, far easier to address.}
While our deep neural networks are highly nonlinear, their solution/output via SGD optimization methods, can be treated as differentiable functions, as the majority of activation functions are differentiable. In this paper, we consider the case of ReLU, which is not differentiable at zero, whereas it is at every other point of $\mathbb{R}$. From a computational standpoint, the gradient descent, heuristically, works well enough to treat it as a differentiable function. Furthermore, goodfellow2016deep argue that this issue is negligible and machine learning softwares are prone to rounding errors, making them very unlikely to compute the gradient at a singularity point. Note that even in this extreme case, both SGD and ADAM, will use the right sub-gradient at zero.
Let the matrix of characteristics be denoted $\boldsymbol{X}_t \in \mathbb{R}^{N\times p}$, where $\boldsymbol{X}_t = (\boldsymbol{x}^{(1)}_{t}, \ldots, \boldsymbol{x}^{(p)}_{t})$. Then for some $i=1, \ldots, N$, $j=1, \ldots, p$ and $t=1, \ldots, T$, the partial derivatives of $g(\boldsymbol{X}_t; \widehat{\theta})$ with respect to the $j^{th}$ characteristic in $\boldsymbol{X}_{t}$ are:
where ${g}(\boldsymbol{X}_t; {\boldsymbol{\widehat{\theta}}})$ is the function (see Section (ref)) that approximates the number of new cases per-$100$K across the $i$ different countries, in our case the G7 countries. We assess the partial derivatives across time since, following kapetanios2007measuring, we expect them to vary due to the inherent nonlinearity of the neural network.
In our work we present the partial derivatives, defined in (ref), without adding confidence bands around them to assess statistical significance. The reason for this is that there is currently no rigorous technology in the literature to produce these, especially in the case of penalized estimation. However, recent work by kap2023 uses a bootstrap approach to construct confidence bands around partial derivatives. A full modification of this work for use in panel models is an interesting and promising avenue to proceed, but is left for future research.
In Figure (ref), we present the partial derivatives with respect to the aggregate stringency index at horizons $h\in \{7, 14, 21, 28 \}.$ Thereby we evaluate the dynamic effectiveness of the stringency policies adopted across the G7 countries.\footnote{ We present the partial derivatives as sixty-day moving averages to smooth out noise.} There are three features that we draw out from Figure (ref). First, policy is more effective at containing the spread of COVID-19 after 7 days. Stronger and more negative effects of increases in stringency are seen after 7 days. Secondly, with the exception of in Japan, policy was most effective in the late fall of 2021 and in early 2022, at the time of the highly contagious Omicron variant. The dynamic effects of policy are, on average, much weaker in the second half of our sample. This is consistent with higher vaccination rates meaning that from mid-2021 (non-immunization) policies became less effective at restraining the spread of new COVID-19 cases. Thirdly, there is considerable cross-country variation in the effectiveness of policy. As referenced above when summarizing the giacomini2010forecast fluctuation tests reported in the online appendix, policy in Japan is again seen in Figure (ref) to have been most effective in late-summer 2022, consistent with COVID-19 cases peaking later in Japan than in the other countries (see Figure (ref)). Containment policies in Italy tended, relative to the other countries, to have a more muted effect.
Given the evidence from Table (ref) that the disaggregated stringency index confers additional forecasting gains relative to the aggregate index, we next look at the partial derivatives with respect to the 9 components of the Oxford index. This way we aim to shed light on the effectiveness of specific policy measures. We focus on the effects of school and university closings and of workplace closings, since of the 9 components of the Oxford stringency index these tend to be the specific policies associated with the largest marginal effects. Results for the other policy measures are provided in the online appendix. Given the high degree of correlation between the different policy measures (see online Tables (ref)-(ref)), we should in any case not over interpret these partial derivatives.
Figure (ref) shows that over time (as $h$ increases from 7 to 28 days) the effects of school and university closings had an increasingly strong effect. For most countries, as expected, these effects are negative: the closures lead to a fall in new COVID-19 cases. These negative effects are especially strong in Italy. But in the UK, the effects are not so clean cut, with the closings appearing to have a positive effect during the early stages of COVID-19. As in Figure (ref), we again see evidence across countries that the effects of school and university closures were far more effective prior to January 2022. Thereafter, the effects are much more modest.
Turning to Figure (ref), we see that while workplace closures tended to have a negative effect on COVID-19 soon after the policy change, in particular in Germany and the UK, thereafter the effects are more uncertain and variable across countries. This can be attributed not just to difficulties in isolating the direct effects of one policy change versus another (related) one, but because in the intervening period there were likely additional and perhaps offsetting changes.
This paper proposes a nonlinear panel data estimator of the conditional mean based on neural networks. We explore heterogeneity and latent patterns in the cross-section, and derive an estimator to account for these patterns. Furthermore, we provide asymptotic arguments for the proposed methodology building on the work of farrell2021deep.
We use the proposed estimators to forecast, in a simulated out-of-sample experiment, the progression of the COVID-19 pandemic across the G7 countries. We find significant forecasting gains over both linear panel data models and time series neural networks. Containment or lockdown policies, as instigated at the national-level, are found to have out-of-sample predictive power for the spread of new COVID-19 cases. Using partial derivatives to help interpret the panel neural networks, we find considerable heterogeneity and time-variation in the effectiveness of specific containment policies.