Deep Neural Network Estimation in Panel Data Models
In this paper we study neural networks and their approximating power in panel
data models. We provide asymptotic guarantees on deep feed-forward neural
network estimation of the conditional mean, building on the work of Farrell et
al. (2021), and explore latent patterns in the cross-section. We use the
proposed estimators to forecast the progression of new COVID-19 cases across
the G7 countries during the pandemic. We find significant forecasting gains
over both linear panel and nonlinear time series models. Containment or
lockdown policies, as instigated at the national-level by governments, are
found to have out-of-sample predictive power for new COVID-19 cases. We
illustrate how the use of partial derivatives can help open the "black-box" of
neural networks and facilitate semi-structural analysis: school and workplace
closures are found to have been effective policies at restricting the
progression of the pandemic across the G7 countries. But our methods illustrate
significant heterogeneity and time-variation in the effectiveness of specific
containment policies.
Ilias Chronopoulos, Katerina Chrysikou, George Kapetanios, James Mitchell, Aristeidis Raftapostolos
# Deep Neural Network Estimation in Panel Data Models+

###### Abstract

In this paper we study neural networks and their approximating power in panel data models. We provide asymptotic guarantees on deep feed-forward neural network estimation of the conditional mean, building on the work of Farrell et al. (2021), and explore latent patterns in the cross-section. We use the proposed estimators to forecast the progression of new COVID-19 cases across the G7 countries during the pandemic. We find significant forecasting gains over both linear panel and nonlinear time series models. Containment or lockdown policies, as instigated at the national-level by governments, are found to have out-of-sample predictive power for new COVID-19 cases. We illustrate how the use of partial derivatives can help open the "black-box" of neural networks and facilitate semi-structural analysis: school and workplace closures are found to have been effective policies at restricting the progression of the pandemic across the G7 countries. But our methods illustrate significant heterogeneity and time-variation in the effectiveness of specific containment policies.

JEL codes: C33, C45.

Keywords: Machine learning, neural networks, panel data, nonlinearity, forecasting, COVID-19, policy interventions.
Introduction

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 Pesaran (2015)) 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, Balestra and Nerlove (1966), the fixed effects (within) estimator, see, for example, Mundlak (1961, 1978), and the Swamy (1970) estimator. Alternative ways to model individual heterogeneity in linear models are found in Hsiao (1974, 1975), with a thorough discussion in Hsiao and Pesaran (2004, 2008) and Part VI of Pesaran (2015).

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 Fernandez-Val and Weidner (2016), who adapt the analytical and jackknife bias correction methods introduced in Hahn and Newey (2004) to nonlinear models with additive or interactive individual and time effects, and Chen et al. (2021) who address estimation and inference in general nonlinear models using iterative estimation. Hacooglu Hoke and Kapetanios (2021) provide an approach for estimation and inference in nonlinear conditional mean panel data models in the presence of cross-sectional dependence. Jochmans (2017) develops the asymptotic properties of GMM estimators for models with two-way multiplicative fixed effects, while Charbonneau (2013) 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, Hornik (1991), Hornik et al. (1989), Gallant and White (1992), and Park and Sandberg (1991). 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, Kapetanios and Blake (2010)).

More recent work by Liang and Srikant (2016) and Yarotsky (2017, 2018) considers feedforward neural networks as approximations for complex functions that accommodate multiple layers, provided sufficiently many hidden neurons and layers are available. Other examples, like Bartlett et al. (2019), provide the theoretical framework for neural network estimation, while Schmidt-Hieber (2020) 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 Farrell et al. (2021), building on the work of Yarotsky (2017) and Bartlett et al. (2019), studies deep neural networks and considers their use for semi-parametric inference.

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_{ii}|\mathbf{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 Farrell et al. (2021), 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 Hastie et al. (2022) and Kelly et al. (2022), as well as Remark 5 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, Athey and Imbens (2017), Wager and Athey (2018), Belloni et al. (2014), Joseph (2019), Chronopoulos et al. (2023), and Kapetanios et al. (2023).

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 2 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 3, 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 4. Section 5 concludes. We relegate to the online appendix additional forecasting results, data summaries, and further discussion of the prediction evaluation tests used.

## 2 Theoretical considerations: the deep neural panel data model

Let \(y_{it}\) be the observation for the \(i^{th}\) cross-sectional unit at time \(t\) generated by the following panel data model:

\[E\left(y_{it}|\mathbf{x}_{it}\right)=\widetilde{h}_{i}\left(\mathbf{x}_{it}\right), \quad i=1,\ldots,N,\;t=1,\ldots,T, \tag{1}\]

where \(\left\{\mathbf{x}_{it}\right\}=\left\{\left(x_{t,1},\ldots,x_{t,p}\right)^{ \prime}\right\}\) 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:

\[y_{it}=\widetilde{h}_{i}\left(\mathbf{x}_{it}\right)+\varepsilon_{it}, \tag{2}\]

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(\mathbf{x}_{it}\right)\) can be decomposed as follows:

\[\widetilde{h}_{i}\left(\mathbf{x}_{it}\right)=h\left(\mathbf{x}_{it}\right)+h_{i} \left(\mathbf{x}_{it}\right), \tag{3}\]

where the function \(h(\cdot)\) is the common component of the model, and is our main focus of interest, and \(h_{i}\left(\mathbf{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:

\[y_{it}=\mathbf{x}_{it}^{\prime}\mathbf{\beta}_{i}+\varepsilon_{it}=\mathbf{x}_{it}^{\prime} \mathbf{\beta}+\mathbf{x}_{it}^{\prime}\mathbf{\eta}_{i}+\varepsilon_{it}, \tag{4}\]

where \(E(\mathbf{\eta}_{i}|\mathbf{x}_{it},\varepsilon_{it})=0\). Equation (4) 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(\mathbf{\cdot};\mathbf{\theta}\right)\). Here, the functional form is known up to the parameter vector \(\mathbf{\theta}\), which is a vector of ancillary parameters, such as network weights and biases. More details on the choice of \(g\left(\mathbf{\cdot};\mathbf{\theta}\right)\) and the role of various neural network parameters will be provided in Section 3 below. Therefore, we parameterize (3) by proposing the following panel model:

\[y_{it}=g\left(\mathbf{x}_{it};\mathbf{\theta}^{0}\right)+g\left(\mathbf{x}_{it};\mathbf{ \theta}_{i}^{0}\right)+\varepsilon_{it}, \tag{5}\]

where \(\mathbf{\theta}^{0}\)and \(\mathbf{\theta}_{i}^{0}\) denote the values of the parameters that best approximate \(h\) and \(h_{i}\), respectively (see (3)), in a sense to be defined below.

It is useful to draw some parallels between (4) and (5). We note that most multi-layer neural network architectures have a final linear layer given by:

\[g\left(\mathbf{x}_{it};\mathbf{\theta}^{0}\right)=\mathbf{\theta}_{L}^{0\prime}\mathbf{f} \left(\mathbf{x}_{it}\right),\]

where \(\mathbf{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 \(\mathbf{f}\) and \(\mathbf{f}_{i}\), of the form:

\[y_{it}=\mathbf{\theta}_{L}^{0\prime}\mathbf{f}\left(\mathbf{x}_{it}\right)+\mathbf{\theta}_{i,L}^{0\prime}\mathbf{f}_{i}\left(\mathbf{x}_{it}\right)+\varepsilon_{it}, \tag{6}\]

which is reminiscent of (4) and thus provides a clear rationale for our nonlinear extension of it.

Furthermore, it provides a rationale for thinking that \(\mathbf{\theta}_{i}^{0}\) plays a similar role to the idiosyncratic coefficients, \(\mathbf{\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 \(\mathbf{x}_{it}\) contains regressors from other cross-sectional units, allowing for complex nonlinear interactions across units. In the limit, each unit could have \(\left(\mathbf{x}_{1t},\ldots,\mathbf{x}_{Nt}\right)\) as the regressor vector.

Next, we consider the conditions needed to identify \(h(\cdot)\). We require certain definitions. First, we define \(\mathbf{\varepsilon}_{it}\equiv y_{it}-h\left(\mathbf{x}_{it}\right)-h_{i}\left(\mathbf{x }_{it}\right)\) and \(u_{it}=h_{i}\left(\mathbf{x}_{it}\right)+\mathbf{\varepsilon}_{it}\), where the latter is in analogy to the usual composite error term for the linear heterogeneous panel model, given as \(\mathbf{x}_{it}^{\prime}\mathbf{\eta}_{i}+\varepsilon_{it}\). The assumption below generalizes the usual identification assumption on \(\mathbf{\eta}_{i}\) made in linear heterogeneous panel models.
**Assumption 1**: _For all \(i=1,\ldots,N,\ t=1,\ldots,T\)_

1. _For some positive constant_ \(C\)_, we assume that_ \(h(\cdot)\) _in (_7_) below is bounded, such that_ \(\|h\|_{\infty}\leq C.\)__\(h_{i}\) _is bounded similarly to_ \(h\)_._
2. \(\{u_{it}\}\) _is independent and bounded across_ \(i\)_._
3. \(E[u_{it}]h\left(\mathbf{x}_{it}\right)]=0.\)__

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 (6) as it is similar, in functionality, to assuming that \(E(\mathbf{\eta}_{i}|\mathbf{x}_{it})=0\).

Next we align our discussion with Farrell et al. (2021). The overall goal of neural network estimation in Farrell et al. (2021) is to estimate an unknown smooth function \(h(\cdot)\) that maps covariates, \(\mathbf{X}\), to an outcome \((T\times N)\) matrix \(\mathbf{Y}\), by minimizing a loss function \(g_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{\theta}\right)\) with respect to the parameterization \(\mathbf{\theta}\) of a neural network function \(g\left(\mathbf{\cdot};\mathbf{\theta}^{0}\right)\). Formally,

\[h=\operatorname*{arg\,min}_{\mathbf{\theta}}E\left[g_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{ \theta}\right)\right]. \tag{7}\]

This is a minimization of a population quantity and assumes that the true function \(h(\cdot)\) is the unique solution of (7). Note that while Farrell et al. (2021) do not specify a true function \(h(\cdot)\), we take a further step and assume that the true functions in (3) coincide with the unique solutions of (7). We do not specify \(\mathbf{Y}\) and \(\mathbf{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_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{\theta}\right)\) in order for our results to hold. We require the following assumptions:

**Assumption 2**: _For some constant \(C_{g_{*}}>0\), we assume that \(g_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{\theta}\right)\) satisfies:_

\[\sup_{\mathbf{Y},\mathbf{X}}|g_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{\theta}_{1}\right)-g_{*} \left(\mathbf{Y},\mathbf{X};\mathbf{\theta}_{2}\right)|\leq C_{g_{*}}\left\|\mathbf{\theta}_{ 1}-\mathbf{\theta}_{2}\right\|,\ \text{for Frobenius norm}\ \left\|\!\left\cdot\right\|\!\right|.\]

**Assumption 3**: _Consider a Holder space \(\mathcal{W}^{b,\infty}([-1,1]^{d})\), with \(b=1,2,\ldots\), where \(\mathcal{W}^{b,\infty}([-1,1]^{d})\) is the space of functions on \([-1,1]^{d}\) in \(L^{\infty}\), along with their weak derivatives. Recall \(h\) in (7), we assume that \(h\) lies within \(\mathcal{W}^{b,\infty}([-1,1]^{d})\), with a norm in \(\mathcal{W}^{b,\infty}([-1,1]^{d})\):_

\[\left\|h\right\|_{\mathcal{W}^{b,\infty}([-1,1]^{d})}=\max_{\mathbf{a}:|\mathbf{a}| \leq b}\operatorname*{ess\,sup}_{x\in[-1,1]^{d}}\left|D^{a}h(x)\right|, \tag{8}\]

_where \(\mathbf{a}=(a_{1},a_{2},\ldots,a_{d})\in[-1,1]^{d}\), \(|\mathbf{a}|=a_{1}+a_{2}+\cdots+a_{d}\) and \(D^{a}h\) is the corresponding weak derivative._
**Remark 1**: _In Assumption 3 we state a smoothness assumption following existing theoretical results; see, for example, Farrell et al. (2021) and Yarotsky (2017, 2018). A more detailed discussion of Holder-Sobolev and Besov spaces is available in Gine and Nickl (2016). Assumption 3 holds for both \(h(\cdot)\) and \(h_{i}(\cdot)\), \(i=1,\ldots,N.\)_

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_{*}\left(\mathbf{Y},\mathbf{X};\mathbf{\theta}\right)=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t =1}^{T}\left(y_{it}-g(\mathbf{x}_{it};\mathbf{\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 Farrell et al. (2021) and the discussion below. Such networks approximate smooth functions well, as shown in Yarotsky (2017, 2018).

A further assumption is required on the processes \(\{\mathbf{x}_{t,k}\}\) and \(\{\varepsilon_{it}\}\), for some \(k=1,\ldots,p\), where \(p\) is the number of covariates.

**Assumption 4**: _We assume the following:_

1. _The rows of_ \(\mathbf{X}_{t}\) _are_ i.i.d. _realizations from a Gaussian distribution whose_ \(p\)_-dimensional inner product matrix_ \(\mathbf{\Sigma}\) _has a strictly positive minimum eigenvalue, such that_ \(\Lambda_{\min}^{2}>0\) _and_ \(\Lambda_{\min}^{-2}=O(1).\)__
2. _The rows of the error term_ \(\varepsilon_{it}\) _are are_ i.i.d. _realizations from a Gaussian distribution, such that_ \(\varepsilon_{it}\sim N(0,\sigma_{\varepsilon}I_{N})\)_._
3. \(\varepsilon_{it}\) _and_ \(\mathbf{x}_{it}\) _are mutually independent._

**Remark 2**: _Assumption 4 states that the covariate process and the errors are continuous and have all their existing moments, while being mutually independent. These strict assumptions are standard in the neural network literature and useful in order to continue with the analysis on a more simplified basis. These assumptions are strict. But it is reasonable to conjecture that similar results to those given below would hold under weaker conditions._

Having rewritten the loss function, as squared error loss, with respect to the re-parameterized panel model in (5), we construct a pooled-type nonlinear estimator, \(\widehat{\mathbf{\theta}}\), such that:

\[\widehat{\mathbf{\theta}}=\operatorname*{arg\,min}_{\mathbf{\theta}}g_{*}\left(\mathbf{Y },\mathbf{X};\mathbf{\theta}\right)=\operatorname*{arg\,min}_{\mathbf{\theta}\in\mathbb{R }^{d}}\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left[y_{it}-g(\mathbf{x}_{it};\mathbf{ \theta})\right]^{2}, \tag{9}\]
which obeys Assumption 2. Therefore, our estimator of \(h\left(\mathbf{x}_{it}\right)\) is given by \(g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})\). Then, we proceed to estimate \(h_{i}\left(\mathbf{x}_{it}\right)\) by \(g(\mathbf{x}_{it};\widehat{\mathbf{\theta}}_{i})\), where \(\widehat{\mathbf{\theta}}_{i}\) is given by:

\[\widehat{\mathbf{\theta}}_{i}=\operatorname*{arg\,min}_{\mathbf{\theta}_{i}\in\mathbb{ R}^{d}}\frac{1}{T}\sum_{t=1}^{T}\left[y_{it}-g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})-g (\mathbf{x}_{it};\mathbf{\theta}_{i})\right]^{2}, \tag{10}\]

for each \(i\), given \(\widehat{\mathbf{\theta}}\) from (9).

Next, we argue that the estimation in (9) can effectively separate \(h\left(\mathbf{x}_{it}\right)\) from \(h_{i}\left(\mathbf{x}_{it}\right)\) and that the unit-wise second step estimation in (10) can retrieve \(h_{i}\left(\mathbf{x}_{it}\right)\). We do this by noting the following. Consider the loss function for \(\widehat{\mathbf{\theta}}\) in (9). We have:

\[\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left[y_{it}-g(\mathbf{x}_{it} ;\mathbf{\theta})\right]^{2} =\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left[\left(h\left(\mathbf{ x}_{it}\right)-g(\mathbf{x}_{it};\mathbf{\theta})\right)+h_{i}\left(\mathbf{x}_{it} \right)+\varepsilon_{it}\right]^{2}\] \[=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left[\left(h\left(\mathbf{ x}_{it}\right)-g(\mathbf{x}_{it};\mathbf{\theta})\right)\right]^{2}+\frac{1}{NT}\sum_{i=1}^{N} \sum_{t=1}^{T}\varepsilon_{it}^{2}\] \[+\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}h_{i}\left(\mathbf{x}_{it} \right)^{2}+\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(h\left(\mathbf{x}_{it} \right)-g(\mathbf{x}_{it};\mathbf{\theta})\right)h_{i}\left(\mathbf{x}_{it}\right) \tag{11}\] \[+\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(h\left(\mathbf{x}_{it} \right)-g(\mathbf{x}_{it};\mathbf{\theta})\right)\varepsilon_{it}+\frac{1}{NT}\sum_{i= 1}^{N}\sum_{t=1}^{T}h_{i}\left(\mathbf{x}_{it}\right)\varepsilon_{it}=\sum_{j=1}^{ 6}A_{j}.\]

Under Assumptions 1-4, 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 (3), \(A_{4}\) is \(O_{p}(N^{-1/2})\). Then it immediately follows that the loss function is minimized when \(\mathbf{\theta}=\mathbf{\theta}^{0}\), in view of our identification assumption in (7). It therefore follows that \(\widehat{\mathbf{\theta}}\rightarrow^{p}\mathbf{\theta}^{0}\) and \(g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})\rightarrow^{p}g(\mathbf{x}_{it};\mathbf{\theta}^ {0})=h\left(\mathbf{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:

\[g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})=\operatorname*{arg\,min}_{\mathbf{\theta}\in \mathbb{R}^{d}}\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(\mathbf{x}_{it};\mathbf{\theta})\right]^{2}=\operatorname* {arg\,min}_{\mathbf{\theta}\in\mathbb{R}^{d}}\frac{1}{N}\sum_{i=1}^{N}\left[\bar{y }_{i}-\bar{g}_{i}(\mathbf{\theta})\right]^{2} \tag{12}\]

and the associated model with a composite error is:

\[\bar{y}_{i}=\bar{g}_{i}(\mathbf{\theta})+u_{i},\quad i=1,\ldots,N, \tag{13}\]

where:

\[u_{i}=\frac{1}{T}\sum_{t=1}^{T}g(\mathbf{x}_{it};\mathbf{\theta}_{i})+ \frac{1}{T}\sum_{t=1}^{T}\varepsilon_{it}=\bar{g}_{i}(\mathbf{\theta}_{i})+\bar{ \varepsilon}_{i}. \tag{14}\]

Note that \(u_{i}\) obeys Assumption 1.2. Moreover, this setting corresponds to that of Theorem 1
in Farrell et al. (2021) enabling the use of the rates derived in this theorem. This analysis is summarized and extended in the following proposition:

**Proposition 1**: _Suppose Assumptions 1-4 hold. Let \(g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})\) be the deep network estimator defined in (12). Then, for some \(\psi<1/2\), the following holds:_

\[\sup_{i,t}\left\|g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})-h\left(\mathbf{x}_{it} \right)\right\|_{2}^{2}=O_{P}(N^{-\psi}). \tag{15}\]

The proof of Proposition 1 follows from the proof of Theorem 1 in Farrell et al. (2021), using the arguments made above the Proposition to recast our panel framework into the one of Farrell et al. (2021), by separately identifying \(h\) and \(h_{i}\). In Proposition 1, we use the results from Theorem 1 of Farrell et al. (2021) to obtain an asymptotic rate of convergence for the error in (15). 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 Farrell et al. (2021). 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 Bartlett et al. (2019) and in Lemma 6 of Farrell et al. (2021). Fast convergence of (15) 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 \(\widehat{\mathbf{\theta}}_{i}\) by minimizing the loss over \(\mathbf{\theta}_{i}\): \(L_{i}=\frac{1}{T}\sum_{t=1}^{T}[y_{it}-g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})-g (\mathbf{x}_{it};\mathbf{\theta}_{i})]^{2}.\) Given the rate in Proposition 1, it immediately follows that \(g\left(\mathbf{x}_{it};\mathbf{\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}\).

**Remark 3**: _Before concluding, it is of interest to consider whether an idiosyncratic component, \(g(\mathbf{x}_{it};\mathbf{\theta}_{i}^{0})\), is needed, in addition to the common component, \(g(\mathbf{x}_{it};\mathbf{\theta}^{0})\). This could be tested by a nonlinear version of a poolability test. One way to proceed is by fitting only the common component and then determining whether the residuals, \(\widehat{u}_{it}=y_{it}-g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})\), can be further explained by unit-wise neural network regressions. Again, one way to do this is by constructing unit-wise \(R^{2}\) statistics. This is intuitive, if we recall the quasi linear representation given by (6). One can regress \(\widehat{u}_{it}\) on \(\mathbf{f}_{i}(\mathbf{x}_{it})\) to obtain such \(R^{2}\) statistics. Then the null hypothesis that \(\widetilde{h}_{i}\left(\mathbf{x}_{it}\right)=h(\mathbf{x}_{it})\), can be tested using the test statistic:_

\[P=\frac{1}{\widehat{\sigma}\sqrt{N}}\sum_{i=1}^{N}\left(TR_{i}^{2}-m\right),\]

_where \(\widehat{\sigma}^{2}=\frac{1}{N}\sum_{i=1}^{N}\left(TR_{i}^{2}-m\right)^{2}\) and an appropriate centering factor, \(m\), needs to be chosen. This could be the dimension of \(\mathbf{f}_{i}\left(\mathbf{x}_{it}\right)\), although care needs to be taken, given that \(\mathbf{f}_{i}\left(\mathbf{x}_{it}\right)\) will contain estimated parameters. One way to resolve this issue may be to estimate the neural networks over a different time period to that used to run the unit-wise regressions of \(\widehat{u}_{it}\) on \(\mathbf{f}_{i}(\mathbf{x}_{it})\). Then, under our assumptions, including that of cross sectional independence, and for an appropriate choice of \(m\), \(P\) is asymptotically standard normal under the null hypothesis._
_Further exploration of this test is of interest. However, a full and rigorous analysis is beyond the scope of the current paper._

## 3 Implementation considerations

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.

### Neural network construction

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 2. 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 (9). 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)^{ \prime}}+\boldsymbol{b}^{(l)^{\prime}})\), while the final output of the network is \(\boldsymbol{g}^{(L)}=\boldsymbol{g}^{(L-1)}\boldsymbol{W}^{(L)^{\prime}}+ \boldsymbol{b}^{(L)^{\prime}}\) 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:

(16)
where \(\mathbf{W}^{(l)}\) is a \(M^{(l)}\times M^{(l-1)}\) matrix of weights, \(\mathbf{b}^{(l)}\) is a \(M^{(l)}\times N\) matrix of biases at layer \(l\), with \(\mathbf{b}^{(1)}=\mathbf{0}\). Notice that at \(l=1\), the dimensions of \(\mathbf{W}^{(1)}\) are \(M^{(1)}\times p\) and of \(\mathbf{b}^{(1)}\) are \(M^{(1)}\times N\). At the final layer, that is, at \(l=L\), the dimensions of \(\mathbf{W}^{(L)}\) are \(1\times M^{(L-1)}\) and of \(\mathbf{b}^{(L)}\) are \(1\times N\).

Note that throughout the paper we use \(\mathbf{\theta}\) to denote a stacked vector containing all ancillary trainable parameters affiliated with the network estimation, as defined below:

\[\mathbf{\theta}=\left(\operatorname{vec}\left(\mathbf{W}^{(1)^{\prime}}\right),\ldots,\operatorname{vec}\left(\mathbf{W}^{(L)^{\prime}}\right),\operatorname{vec} \left(\mathbf{b}^{(1)^{\prime}}\right),\operatorname{vec}\left(\mathbf{b}^{(2)^{ \prime}}\right),\ldots,\mathbf{b}^{(L)^{\prime}}\right)^{\prime}. \tag{17}\]

We define the overall number of parameters as \(d=|\mathbf{\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 \((\mathbf{W}^{(l)},\mathbf{b}^{(l)})\) are calculated through back-propagation (using the chain-rule) to train the network.

**Remark 4**: _The exact (composition) structure described in (16) holds for a subclass of feed-forward neural networks, specifically that one that refers to fully connected layers (the one being consecutive to the other) but has no other connections. Each layer has a number of hidden units that are of the same order of magnitude. This architecture is the most commonly used in empirical research, and is often referred to as a Multi-layer Perceptron (MLP). Furthermore, the exact structure in (16) does not hold generally for any feed-forward neural network._

The specific choice of the network architecture is crucial and affects the complexity and the approximating power of \(g\left(\mathbf{\cdot};\mathbf{\theta}\right)\) in (16). 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 (16), in Figure 1.

The neural network in Figure 1 consists of two inputs \(\mathbf{X}_{t}\in\mathbb{R}^{N\times p}\), \(\mathbf{X}_{t}=(\mathbf{x}_{t}^{(1)},\mathbf{x}_{t}^{(2)})\), in particular, \(p=2\), where \(\mathbf{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 \(\widehat{\mathbf{y}}_{t}\). Between the inputs and output \((\mathbf{X}_{t},\widehat{\mathbf{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.1

Figure 1: Illustration of a _feed-forward neural network_ with two input matrices \((\mathbf{x}_{t}^{(1)},\mathbf{x}_{t}^{(2)})^{\prime}\), two layers, \(L=2\), \(5\) nodes, \(M=5\), eighteen connections, \(W=14\), and one (fitted) output \(\widehat{\mathbf{y}}_{t}.\) The inputs are illustrated with a white circle, the neurons with grey circles, the output with a black circle.
Notice that the illustration in Figure 1 can correspond to a nonlinear pooled-type estimation of \(g(\mathbf{X}_{t};\mathbf{\theta}^{0})\) in (5), where we use (9) to obtain \(\widehat{\mathbf{\theta}}\), defined in (17), with input \(\mathbf{X}_{t}=(\mathbf{x}_{t}^{(1)},\ldots,\mathbf{x}_{t}^{(p)})\), and output \(\widehat{\mathbf{y}}_{t}\in\mathbb{R}^{N\times 1}.\) The remainder of (5), \(g(\mathbf{X}_{t};\mathbf{\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(\mathbf{x}_{it};\widehat{\mathbf{\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.

### Implementation and regularization

In this section we discuss some operational implementation aspects required for the estimation of the panel neural network estimators proposed in Section 2. We focus discussion on the following panel estimator, \(g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})\), obtained from the optimization of (12):

\[g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})=\operatorname*{arg\,min}_{\mathbf{\theta}\in \mathbb{R}^{d}}\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(\mathbf{x}_{it};\mathbf{\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 Goodfellow et al. (2016) and Gu et al. (2020, 2021). In order to be able to minimize (12) and obtain a feasible solution for the panel estimator \(g(\mathbf{x}_{it};\widehat{\mathbf{\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_{*}(\mathbf{y},\mathbf{X};\mathbf{\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, \(\mathbf{\theta}\). The penalized estimator based on the LASSO is obtained as the solution to the following minimization problem:

\[g(\mathbf{x}_{it};\widehat{\mathbf{\theta}})^{\text{LASSO}}=\operatorname*{arg\,min} _{\mathbf{\theta}\in\mathbb{R}^{d}}\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(\mathbf{x}_{it};\mathbf{\theta})\right]^{2} +\lambda\left\|\mathbf{\theta}\right\|_{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, Zhang et al. (2021). 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 Gu et al. (2020) for a detailed discussion of early stopping and dropout.

Batch normalization, proposed by Ioffe and Szegedy (2015), 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.

**Remark 5**: _In this paper, we consider both penalized and non-penalized estimation. While using the latter might seem problematic due to the large number of parameters that needs to be estimated, we find, in our empirical work, that this is not necessarily the case. This is not surprising. Recent work in the statistical and machine learning literature highlights what is known as the double descent effect. For linear regressions, this relates to the use of generalized inverses to construct least squares estimators, when the number of variables, \(p\), exceeds the number of observations, \(T.\) Such estimators work better either when \(p\) is small (and standard matrix inversion can be used) or when \(p\) is much larger than \(T.\) Then the quality of the performance of an estimator is implicitly measured in terms of the "bias-variance trade-off," where an optimal performance resides at the lowest reported bias and variance of the corresponding model (either linear or nonlinear). While it is widely accepted that the "bias-variance trade-off" function resembles a U-shaped curve, it has been observed, for example see Belkin et al. (2019) and Hastie et al. (2022), that beyond the interpolation limit the test loss descends again, hence the name "double-descent." To understand why, note that such estimators implicitly impose penalization by using generalized inverses and so choose the parameter vector with the smallest norm, among all admissible vectors; see Hastie et al. (2022). So once \(p\) is much larger than \(T\) such a selection becomes more consequential, as many more candidate vectors are admissible. This linear effect is also present for neural network estimation given the connection to linear models highlighted in Section 2, as discussed in detail in Hastie et al. (2022) and Kelly et al. (2022)._

### Cross-validation

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 \([0.01,0.001]\). 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, Gu et al. (2020).

To select the trainable parameters, \(\mathbf{\theta}\), and the hyperparameters discussed above, we follow Gu et al. (2020, 2021) 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, \(\mathbf{\theta}\), given a specific set of hyperparameters; the _validation_ sub-sample, which is used to tune the different hyperparameters given \(\widehat{\mathbf{\theta}}\) from the _training_ sub-sample2; 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^{*}\) denote the total sample size for the specific window, then the _training_ sub-sample consists of \(\lfloor 0.8T^{*}\rfloor\), the validation sub-sample consists of \(\lfloor 0.2T^{*}\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.

Footnote 2: Note that while \(\widehat{\mathbf{\theta}}\) is used in the tuning of the hyperparameters, it is only estimated at the _training_ sub-sample.

### Optimization

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 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 Kingma and Ba (2014)3, which is a more efficient version of SGD. Finally, we set the number of epochs to \(5,000\) and use early stopping following Gu et al. (2020) to mitigate potential overfitting.

Footnote 3: ADAM is using estimates for the first and second moments of the gradient to calculate the learning rate.

## 4 Empirical analysis: forecasting new COVID-19 cases

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.

### The COVID-19 data and the Oxford stringency index

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 \(\mathbf{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. Knutson et al. (2023), Mathieu et al. (2021), and Caporale et al. (2022) 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 \(\mathbf{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, Gu et al. (2021)) and rank-normalize all of our variables into the \([0,1]\) interval as follows:

\[\widetilde{\mathbf{x}}_{i}=\frac{\mathbf{x}_{i}-\min(\mathbf{x}_{i})}{\max(\mathbf{x}_{i})-\min( \mathbf{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 2 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.

Figure 2: The Oxford stringency index and new COVID-19 cases per 100K of the population
#### 4.1.1 Out-of-sample forecasting design

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 (5), the \(h-\)day-ahead forecast of new COVID-19 cases per 100K is:

\[\widehat{y}_{i,t+h}|\mathcal{F}_{t}=\widehat{g}\left(\mathbf{x}_{it};\mathbf{\theta}^{ \star}\right)+\widehat{g}\left(\mathbf{x}_{it};\mathbf{\theta}_{i}^{\star}\right), \tag{18}\]

where \(\widehat{g}\left(\mathbf{x}_{it+h};\mathbf{\theta}^{\star}\right)\) denotes the corresponding fit of the pooled network, \(\widehat{g}\left(\mathbf{x}_{it+h};\mathbf{\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\), \(\mathbf{\theta}^{\star}\) denotes the optimal weights obtained from the CV for the deep pooled model, and \(\mathbf{\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(\mathbf{x}_{it};\mathbf{\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 \(\mathbf{x}_{it}\), as per (18)). To ease the computational burden, given that we re-estimate the model and use CV (as discussed in Section 3) 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 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 Canova and Ciccarelli (2009) and specify the \(i^{th}\) equation of the PVAR with \(q\) lags as:

\[y_{it}=A_{1i}\boldsymbol{Y}_{t-1}+\cdots+A_{qi}\boldsymbol{Y}_{t-q}+\epsilon_{ it},\epsilon_{it}\sim\text{i.i.d.}N\left(0,\sigma_{i}\right), \tag{19}\]

where \(A_{ji}\) for \(j=1,\ldots,q\) are coefficient matrices, we have dropped the intercept for notational simplicity, \(\boldsymbol{Y}_{t}=(z_{1t}^{\prime},\ldots,z_{Nt}^{\prime})^{\prime}\), and \(z_{it}=(y_{it},\boldsymbol{x}_{it})^{\prime}\). We set \(q=28\). We estimate the PVAR by OLS and compute \(h\)-day-ahead forecasts of \(y_{it+h}\) from (19) via iteration.

#### 4.1.2 Forecast evaluation

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):

\[\text{RMSE}_{i}=\sqrt{\frac{1}{T}\underset{t=1}{\overset{T}{\sum}}\left(y_{ i,t+h}-\widehat{y}_{i,t+h}|\mathcal{F}_{t}\right)^{2}},\quad i=1,\ldots,N.\]

We use the Diebold and Mariano (1995) (DM) test to test whether differences in forecast accuracy across models are statistically significant. We follow Harvey et al. (1997) and use their small-sample adjustment.

Table 1 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 B.1-B.3 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 2 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 1, 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 2 shows that after an additional two weeks policy was effective in all G7 countries, except Italy and the US.

Table 3 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 3 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 Giacomini and Rossi (2010) 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.

\begin{table}
\begin{tabular}{l l l l l l l l} \hline \hline  & Canada & France & Germany & Italy & Japan & UK & US \\ \hline  & & & \(h=7\) & & & & \\ \hline Deep pooled & 1.084 & 0.928 & 1.033 & 1.229 & 0.878 & 1.043 & 1.321 \\ Deep idiosyncratic & 1.015 & 0.897 & 1.038 & 1.152 & 0.836 & 1.309 & 1.329 \\ Deep time-series & 0.847 & 0.946 & 1.076 & 0.881 & 1.044 & 0.949 & 1.040 \\ PVAR(28) & 0.875\({}^{***}\) & 0.868\({}^{***}\) & 0.933 & 0.850\({}^{*}\) & 0.792\({}^{**}\) & 0.880\({}^{*}\) & 0.596\({}^{**}\) \\ \hline  & & & \(h=14\) & & & & \\ \hline Deep pooled & 1.012 & 0.907\({}^{*}\) & 0.896 & 1.154 & 0.861 & 0.952 & 1.225 \\ Deep idiosyncratic & 0.967 & 0.868\({}^{*}\) & 0.913 & 1.087 & 0.841 & 1.210 & 1.191 \\ Deep time-series & 0.910 & 0.915 & 1.085 & 0.920 & 1.045 & 0.920 & 0.960 \\ PVAR(28) & 0.880\({}^{***}\) & 0.856\({}^{***}\) & 0.930 & 0.855\({}^{*}\) & 0.773\({}^{**}\) & 0.891\({}^{*}\) & 0.603\({}^{**}\) \\ \hline  & & & \(h=21\) & & & & \\ \hline Deep pooled & 0.953 & 0.879\({}^{*}\) & 0.864\({}^{*}\) & 1.117 & 0.855\({}^{*}\) & 0.957 & 1.078 \\ Deep idiosyncratic & 0.903 & 0.838\({}^{**}\) & 0.874 & 0.989\({}^{*}\) & 0.850 & 1.124 & 1.059 \\ Deep time-series & 0.951 & 0.878 & 1.094 & 0.923 & 1.033 & 0.936 & 0.947 \\ PVAR(28) & 0.887\({}^{***}\) & 0.839\({}^{***}\) & 0.928 & 0.854\({}^{*}\) & 0.794\({}^{**}\) & 0.889\({}^{*}\) & 0.607\({}^{***}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: RMSE ratios, comparing the forecast accuracy of each respective model with and without the aggregate Oxford stringency index at 7, 14, and 21 days-ahead. Ratios \(<1\) indicate superior predictive ability for the model with the stringency index. For a description of the 4 forecasting models, see the notes to Table 1. \(*\), \(**\), and \(***\) denote rejection of the null hypothesis of equality of forecast mean squared errors with and without the aggregate Oxford stringency index at the 10%, 5%, and 1% levels of significance, respectively, using the modified Diebold and Mariano (1995) test with the Harvey et al. (1997) adjustment.
### Policy effectiveness

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

Footnote 4: We prefer the use of partial derivatives over Shapley additive explanation values, as proposed by Lundberg and Lee (2017), since derivatives tend to be less noisy (see, for example, Chronopoulos et al. (2023)) 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 Aas et al. (2021), who also discuss solutions. However, these are computationally intensive, potentially still quite poor approximations, and not appropriate for large sets of variables.

\begin{table}
\begin{tabular}{l l l l l l l l} \hline \hline  & Canada & France & Germany & Italy & Japan & UK & US \\ \hline \multicolumn{8}{c}{\(h=7\)} \\ \hline Deep pooled & 0.941 & 0.834\({}^{**}\) & 0.852\({}^{**}\) & 0.848 & 0.852\({}^{*}\) & 0.842\({}^{**}\) & 0.845\({}^{**}\) \\ Deep idiosyncratic & 0.866 & 0.861 & 0.913 & 0.861 & 0.875 & 1.059 & 0.947 \\ Deep time-series & 0.906 & 1.063 & 1.071 & 0.909 & 1.067 & 0.901 & 1.075 \\ PVAR(28) & 0.840 & 0.614\({}^{***}\) & 0.992 & 0.791 & 0.807 & 0.850 & 0.739 \\ \hline \multicolumn{8}{c}{\(h=14\)} \\ \hline Deep pooled & 0.921 & 0.893 & 0.868\({}^{**}\) & 0.843\({}^{**}\) & 0.860\({}^{*}\) & 0.887\({}^{*}\) & 0.908\({}^{*}\) \\ Deep idiosyncratic & 0.796\({}^{*}\) & 0.900 & 0.955 & 0.822\({}^{*}\) & 0.877 & 1.021 & 0.925 \\ Deep time-series & 0.934 & 1.079 & 1.066 & 0.917 & 1.096 & 0.890 & 1.045 \\ PVAR(28) & 0.826 & 0.579\({}^{***}\) & 1.043 & 0.783 & 0.802 & 0.840 & 0.732\({}^{*}\) \\ \hline \multicolumn{8}{c}{\(h=21\)} \\ \hline Deep pooled & 0.900 & 0.884\({}^{**}\) & 0.845\({}^{**}\) & 0.825\({}^{***}\) & 0.843\({}^{**}\) & 0.931 & 0.887\({}^{***}\) \\ Deep idiosyncratic & 0.772\({}^{**}\) & 0.869 & 0.991 & 0.787\({}^{**}\) & 0.891 & 0.953 & 0.864 \\ Deep time-series & 0.994 & 1.090 & 1.118 & 0.919 & 1.084 & 0.906 & 1.009 \\ PVAR(28) & 0.830 & 0.557\({}^{***}\) & 1.078 & 0.785 & 0.846 & 0.834 & 0.730\({}^{**}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: RMSE ratios, comparing the forecast accuracy of each respective model with and without the disaggregate Oxford stringency index at 7, 14, and 21 days-ahead. Ratios \(<\!1\) indicate superior predictive ability for the model with the stringency index. For a description of the 4 forecasting models, see the notes to Table 1. \(*\), \(**\), and \(***\) denote rejection of the null hypothesis of equality of forecast mean squared errors with and without the disaggregate Oxford stringency index at the 10%, 5%, and 1% levels of significance, respectively, using the modified Diebold and Mariano (1995) test with the Harvey et al. (1997) adjustment.
Footnote 4: We use the notation of [21], which is the same as the one used in [21], but it is not clear that the \(\mathbf{X}_{t}\) is not a \(\mathbf{X}_{t}\)-invariant.

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, Goodfellow et al. (2016) 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 \(\mathbf{X}_{t}\in\mathbb{R}^{N\times p}\), where \(\mathbf{X}_{t}=(\mathbf{x}_{t}^{(1)},\ldots,\mathbf{x}_{t}^{(p)})\). Then for some \(i=1,\ldots,N\), \(j=1,\ldots,p\) and \(t=1,\ldots,T\), the partial derivatives of \(g(\mathbf{X}_{t};\widehat{\theta})\) with respect to the \(j^{th}\) characteristic in \(\mathbf{X}_{t}\) are:

\[d_{i\,j,\,t}=\frac{\partial g\left(\mathbf{X}_{t};\widehat{\mathbf{\theta}}\right)}{ \partial x_{i,\,j,\,t-h}}, \tag{20}\]

where \(g(\mathbf{X}_{t};\widehat{\mathbf{\theta}})\) is the function (see Section 2) that approximates the number of new cases per-100K across the \(i\) different countries, in our case the G7 countries. We assess the partial derivatives across time since, following Kapetanios (2007), we expect them to vary due to the inherent nonlinearity of the neural network.

In our work we present the partial derivatives, defined in (20), 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 Kapetanios et al. (2023) 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 3, 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.5 There are three features that we draw out from Figure 3. 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 Giacomini and Rossi (2010) fluctuation tests reported in the online appendix, policy in Japan is again seen in Figure 3 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 2). Containment policies in Italy tended, relative to the other countries, to have a more muted effect.

Given the evidence from Table 3 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 A.2-A.3), we should in any case not over interpret these partial derivatives.

Figure 4 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 3, we again see evidence across countries that the effects of school and university closures were far

Figure 3: Partial derivatives: The effects of policy (as measured by the Oxford stringency index) on new COVID-19 cases 7, 14, 21, and 28 days after the policy change.
more effective prior to January 2022. Thereafter, the effects are much more modest.

Turning to Figure 5, 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.

Figure 4: Partial derivatives: The effects of school and university closures on new COVID-19 cases 7, 14, 21, and 28 days after the policy change

Figure 5: Partial derivatives: The effects of workplace closures on new COVID-19 cases 7, 14, 21, and 28 days after the policy change
Conclusion

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 Farrell et al. (2021).

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.
