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.
118,255 characters · 16 sections · 89 citation commands
Causal inference and policy evaluation without a control group$^*$
Keywords: potential outcomes framework, counterfactual forecasting, machine learning, short panels, panel cross-validation, educational inequality.
JEL-Codes: C18, C21, C53, I24.
---------------------
\newgeometry{left=2.25cm,right=2.25cm,top=2.5cm,bottom=2.5cm}
\frenchspacing Today, the econometric toolbox for causal panel analysis features many alternative approaches for estimating causal effects, such as difference-in-differences Card:Krueger:1994, the synthetic control method Abadie:Diamond:Hainmueller:2010, two-way Angrist:Pischke:2009 and interactive Bai:2009 fixed effects models, and matrix completion methods Athey:Bayati:Doudchenko:Imbens:Khosravi:2021. However, all of these popular methodologies depend on a critical requirement: the availability of untreated units forming a credible control group, without which they cannot be applied. This poses a significant empirical challenge in observational studies, as there are at least two relevant cases in which a suitable control group does not exist: (i) when the treatment simultaneously affects all units, such as in the context of a large-scale shock or a nationwide program with universal participation Duflo:2017; (ii) when only a subgroup of units receives the treatment, but the set of untreated units cannot form a valid control group due to violations of the no-interference assumption Cox:1958.\footnote{Spillovers are common in many empirical settings where the observations are correlated in time, space, or both, and ignoring them can lead to very misleading inferences Sobel:2006. In principle, spillovers can be modeled and accounted for at the cost of additional and often restrictive assumptions. However, in practice, most methods just postulate the absence of interference.} Under these challenging but not uncommon circumstances, the use of standard causal panel data methods is precluded, leaving a gap in the econometric toolkit of empirical researchers.
\frenchspacing We fill this gap by introducing the Machine Learning Control Method (MLCM), a new estimator based on flexible counterfactual forecasting via machine learning (ML). The MLCM is a versatile technique that can leverage any available supervised ML algorithm and enables the identification and estimation of several policy-relevant causal parameters---including individual, average, and conditional average treatment effects (CATEs)---in practically relevant environments without a control group. It is flexible in terms of data structure, as it can be applied across a wide variety of panel settings, including very short panels and contexts with staggered adoption. The intuition behind the MLCM is straightforward: when you cannot rely on untreated units to build the counterfactual, you can forecast it. Specifically, use supervised ML techniques trained and tuned on the pre-treatment data to forecast the evolution of post-treatment outcomes in the absence of the treatment. Then, estimate treatment effects as the difference between observed and forecasted post-treatment outcomes.
Our approach builds upon and bridges three different methodological currents: causal panel data methods, time series forecasting, and causal ML. From the literature on causal panel analysis, we adopt the conceptual and identification framework based on the potential outcomes model Rubin:1974. Since earlier contributions Ashenfelter:Card:1985, Card:1990 to its more recent developments Arkhangelsky:Athey:Hirshberg:Imbens:Wager:2021, Roth:Santanna:Bilinski:Poe:2023, this entire literature revolves around addressing the “fundamental problem of causal inference” Holland:1986 in an observational panel data setting where, for each unit, it is impossible to observe both potential outcomes at once. Therefore, causal inference can be viewed as a missing data problem, with the key goal being the imputation of missing potential outcomes for the treated units Imbens:Rubin:2015. Most widely used causal panel data methods achieve this imputation by leveraging untreated units to construct a control group. Our key contribution to this literature is to establish identification conditions in the absence of a control group and to introduce a novel estimator of causal effects that imputes missing potential outcomes without relying on untreated units.
The time series literature features some forecasting approaches that researchers have employed for counterfactual building in scenarios without a control group.\footnote{A rich tradition also exists for forecasting with panel data (e.g., arellano2011nonlinear, baltagi2008econometric, liu2020forecasting). However, our interest lies in counterfactual forecasting. In this regard, the Rubin causal framework is absent from this tradition, which focuses on different estimands and non-causal settings compared to the recent causal panel data literature arkhangelsky2024causal.} Some of them are straightforward: the mean method estimates the counterfactual “no-treatment scenario” of the dependent variable $\operatorname{Y}_{i}$ as the average of its pre-treatment values, while the naive method considers the last pre-treatment value of $\operatorname{Y}_{i}$ as the counterfactual Hyndman:Athanasopoulos:2021. A more sophisticated forecasting approach is the interrupted time series (ITS) analysis, first introduced by Box:Tiao:1975. ITS analysis generally fits regression models or autoregressive integrated moving average (ARIMA) models to the entire time series of observed data; this is done after postulating a structure on the intervention effect (e.g., constant level shift). By construction, ITS delivers a single average effect for the whole post-intervention period. More recent studies Brodersen:Gallusser:Koehler:Remy:Scott:2015, Chernozhukov:Wuthrich:Zhu:2021,Menchetti:Cipollini:Mealli:2023, rambachan2021common consider estimation and inference with other time series techniques (such as Bayesian Structural Time Series, ARIMA models, or impulse response functions) in time series settings where no direct control group may be available. From a panel data perspective, these methods have significant limitations. First, they are designed for cases with a single treated series. Second, they often impose restrictive assumptions and functional forms. Third, they can be implemented only when many pre-treatment periods are available. However, in the traditional panel data literature, the number of units is much larger than the number of time periods arkhangelsky2024causal, and the time dimension is typically substantially smaller than in time series settings, rendering these methods mostly inapplicable. We draw from this literature the idea of forecasting counterfactuals by exploiting the temporal information. We go beyond it by proposing a more flexible, ML-powered method explicitly designed for counterfactual forecasting with panel data.
In recent years, a new methodological literature at the intersection between causal inference with panel data and ML has evolved. The matrix completion methods developed by Athey:Bayati:Doudchenko:Imbens:Khosravi:2021 use nuclear norm regularization and observed control outcomes to impute missing elements in the counterfactual matrix of untreated unit/period combinations. semenova2017estimation propose an approach for estimating CATEs characterized by high-dimensional parameters in both homogeneous cross-sectional and unit-heterogeneous dynamic panel data settings, leveraging ML techniques for improved inference. Artificial control methods Carvalho:Masini:Medeiros:2018,Masini:Medeiros:2021 and synthetic learners Viviano:Bradic:2023 harness supervised ML algorithms to predict counterfactuals and assess treatment effects in a panel setting with a single treated unit, many untreated units, and a long pre-intervention window. All these methods need a set of units assumed to be completely unaffected by the treatment to build a counterfactual scenario. Therefore, none of these causal ML techniques can be applied in the challenging econometric setting we focus on. We share with this literature the idea of harnessing the power of ML in the service of causal inference rather than for pure prediction, exploiting the fact that counterfactual building is ultimately a predictive task Varian:2016. We extend it by developing a new causal ML estimator that estimates individual treatment effects without depending on the no-interference assumption or leveraging unaffected units.
We begin by thoroughly discussing and formalizing identification without a control group in the potential outcomes framework. We then propose an adaptive estimator based on a horse-race competition between several different ML algorithms, a novel cross-validation (CV) procedure for model assessment and selection tailored for panel data, and block-bootstrapping for inference. The MLCM comes with a set of diagnostic and placebo tests and is characterized by a high level of generality: it delivers individual treatment effects that can either be the object of interest or aggregated into several policy-relevant causal estimands. Since our approach allows for unrestricted treatment effect heterogeneity, we also propose an automated search for heterogeneity based on an easy-to-interpret regression tree.
To showcase the potential of the MLCM, we present an extensive simulation study, a replication of the minimum wage application in callaway2021difference without using control units, and an empirical application in which we investigate the effects of the COVID-19 pandemic on educational inequality in Italian Local Labor Markets (LLMs). We find that the pandemic led to a generalized drop in students’ performance, which is particularly pronounced in LLMs that before COVID-19 were characterized by higher levels of unemployment and inequality and lower educational attainment.
At the time of writing, we are aware of only one other method for assessing causality in panel settings without using control units.\footnote{There are also some empirical studies in energy economics that, without proposing a new formal method, have applied counterfactual prediction using ML algorithms on very high-frequency (hourly) electricity data abrell2022effective, jarvis2022private, prest2023rcts.} In independent work subsequent to ours, Botosaru:Giacomini:Weidner:2023 propose estimating average treatment effects in the absence of a control group by regressing pre-treatment outcomes on known basis functions of time. Compared to the MLCM, this method relies on different assumptions, such as absence of interference and the validity of a central limit theorem, and imposes functional form restrictions. Overall, we view the two methods as complementary, depending on the assumptions one is willing to make and on the characteristics of the available data.\footnote{For instance, the estimator by Botosaru:Giacomini:Weidner:2023 can be employed in data-scarce environments where only outcome data are available, or when the researcher is comfortable relying on parametric assumptions. In contrast, the MLCM is more flexible and better suited for data-rich or high-dimensional settings where information on many covariates is available. This context is particularly relevant when the primary goal is to uncover heterogeneity in treatment effects associated with differences in observed characteristics.}
The possibility of moving beyond the reliance on untreated units paves the way to evaluating many potential treatments---such as universal policies, large-scale shocks, and programs engendering interactions between units---that, due to the lack of a valid comparison group, have so far been under-explored in empirical studies. To answer these causal questions, interested researchers can implement the MLCM on real-world datasets using the companion R package MachineControl.\footnote{The latest version of the package is available on GitHub at \href{https://github.com/FMenchetti/MachineControl}{\textcolor{blue}{this link}}.}
The rest of this paper is organized as follows. Section (ref) outlines the causal framework by defining the identification assumptions, the causal estimands, and the estimators. Section (ref) describes the implementation process of the MLCM and presents the simulation study. Section (ref) illustrates the empirical application, while Section (ref) concludes.
In this section, we present the causal framework for an observational short-panel setup where the intervention is a simultaneous policy change or shock that affects directly or indirectly the entire statistical population under scrutiny.\footnote{To use language from arkhangelsky2024causal, in this benchmark case, we consider a thin panel matrix, with a large number of cross-sectional units and a small number of time periods. However, the framework easily accommodates the use of long panels and fat and square matrices. In addition, the MLCM can be applied to both balanced and unbalanced panels, as well as in block-assignment and staggered adoption settings where one cannot credibly rely on the plausibility of the no-interference assumption (i.e., outcomes of never-treated or not-yet-treated units are contaminated by spillovers).} We start by defining the assumptions and the corresponding causal estimands within the potential outcomes framework. We then provide the proof that such estimands can be identified under those assumptions. Finally, we introduce novel estimators for scenarios without a suitable control group. To make the notation (and the method) as general as possible, the definitions given in this section assume the presence of some contemporaneous covariates unaffected by the intervention. Including these covariates in the MLCM can potentially mitigate certain assumptions. However, as we will discuss later, caution must be exercised in their use and selection.
Denote with $\operatorname{Y}_{i,t}$ the outcome of unit $i$ at time $t$, and let $\operatorname{W}_{i,t} \in \{0,1\}$ be a random variable describing the treatment assignment of unit $i = 1, \dots, N$ at time $t = 1, \dots,t_0,\dots, T$, where $1$ indicates the treatment, $0$ indicates control, and $t_0$ denotes the intervention date.\footnote{The word “treatment” is commonly used in the context of randomized controlled trials. As we are dealing with an observational study, we use the terms “treatment”, “intervention”, “policy”, and “shock” interchangeably.} As we are focusing on a single and simultaneous intervention, we can then write $\operatorname{W}_{i,t} = 0$ for all $i$ and $t \leq t_0$, $\operatorname{W}_{i,t} = 1$ for all $i$ and $t > t_0$, therefore, $t_0+1$ is the first post-intervention period. In other words, all units are treated simultaneously and remain treated thereafter. Under the potential outcomes framework, for each unit, we can define a potential outcome under $\operatorname{W}_{i,t} = 0$ and a different potential outcome under $\operatorname{W}_{i,t} = 1$.
We now introduce the assumptions that are needed to define, identify, and estimate causal effects in an observational panel setting without controls, offering novel theoretical insights such as identification conditions for non-linear, multi-step-ahead causal effects. We emphasize that such assumptions are not particularly restrictive, as the MLCM is designed to be highly flexible. Our method, whose implementation is detailed in Section (ref), starts by running several ML algorithms on pre-treatment data and then selects the one producing the most accurate forecasts. Given that it is not possible to know in advance which ML method will be chosen, the notation is kept as general as possible.\footnote{We believe this is an advantage over most existing approaches in the causal ML literature: while Carvalho:Masini:Medeiros:2018 and Masini:Medeiros:2021 focus on LASSO and Wager:Athey:2018 adopt tree-based methods, the MLCM is flexible and can easily adapt to different data-generating processes. A notable exception in this context is the synthetic learner proposed by Viviano:Bradic:2023. In contrast to their technique, our method is based on a horse-race competition between several alternative ML algorithms, whereas the synthetic learner relies on an ensemble procedure that combines various estimators. }
The first assumption contributes to define potential outcomes and establishes the crucial link between them and the treatment assignment. We rely on a less stringent version of the Stable Unit Treatment Value Assumption (SUTVA) Rubin:1974, Imbens:Rubin:2015, retaining only its second part: the treatment is the same for all the units. This is an a priori assumption that the potential outcomes do not depend on different treatment intensities or the mechanism used to assign the treatment.
While we maintain this no-multiple-versions-of-treatment assumption, we drop the first part of SUTVA, namely, the no-interference assumption Cox:1958.\footnote{Following Imbens:Rubin:2015, we consider the case with general equilibrium effects as a scenario in which there are widespread violations of the no-interference assumption.} We consider this as one of the main advantages of our approach because, while it has long been known that violations and failures of the no-interference assumption can lead to misleading inferences in many social science settings Sobel:2006, this strong assumption is rarely questioned or tested in practice Chiu:Lan:Liu:Ziyi:Xu:2023. This is a key departure from most evaluation methods and it is important to delve into its implications. SUTVA-related interference among units can be of two types: (1) from treated to control units, and (2) among treated units themselves.\footnote{Even if the treatment is the same for all units, residual spillover effects may arise due to interference among the treated units. Consider, for example, a study investigating excess mortality from COVID-19. While the pandemic affects all areas, those with fewer intensive care beds may need to send patients to neighboring hospitals, prompting an additional (indirect) rise in mortality rates there also. Spillovers due to individuals’ characteristics within the same treatment group are discussed in Ogburn:VanderWeele:2014. } First, we completely circumvent pitfalls regarding Type-1 interference, since our counterfactual scenario is generated without relying on untreated units. Second, we do not postulate the lack of interference among treated units and allow for any possible Type-2 interference. We avoid relying on the no-interference-among-treated assumption because of the likely presence of interference across both the temporal and cross-sectional dimensions in social science applications Xu:2023. Consequently, our focus is on estimating the total effect of the treatment, which encompasses both the direct effect on a unit from the treatment it receives and the indirect effects arising from the spillover and general equilibrium effects. We claim that in real-world scenarios with likely interaction between units, focusing on the total effect of the treatment is a sensible choice as it captures the actual impact on each unit under the realized treatment assignment mechanism. Under Assumption (ref), we can index the potential outcomes to the common treatment path: $\operatorname{Y}_{i,t}(1)$ indicates the potential outcome under the treatment, and $\operatorname{Y}_{i,t}(0)$ is the counterfactual outcome.
The next assumption contributes to the definition of causal effects.
This is not a restrictive assumption, since it holds after a suitable transformation of the data. For example, we can start from a multiplicative effect and then apply a logarithmic transformation to find an additive structure. Since inference is conducted using block-bootstrap, if it is necessary to express the effect in its original form, the inverse transformation can be applied to the bootstrap distribution, so as to recover the original multiplicative effect and its confidence interval. In conjunction with Assumption (ref) (see below), this assumption implies that we effectively adopt a Partially Linear Model, akin to seminal causal machine learning methods chernozhukov2018double, Wager:Athey:2018. The Partially Linear Model strikes a balance between structure and flexibility: the causal-effect component of the model remains simple and interpretable, while the untreated potential outcome can exhibit almost arbitrary complexity chernozhukov2024applied.
Assumption (ref) rules out the possibility of anticipatory effects of the intervention on both the outcome and the covariates Abadie:Diamond:Hainmueller:2010, Borusyak:Jaravel:Spiess:2021.
This assumption implies that in the pre-intervention period, the observed outcome corresponds to the potential outcome absent the policy, i.e., $\operatorname{Y}_{i,t} \equiv \operatorname{Y}_{i,t}(0)$ and that $\mathbf{X}_{i,t} \equiv \mathbf{X}_{i,t}(0)$ throughout all the analysis period. Additionally, Assumption (ref) precludes any anticipatory actions by agents, implying that they cannot alter or manipulate their outcomes prior to receiving the treatment. This assumption can be empirically tested by verifying whether pre-treatment effects are, on average, zero. Finally, post-treatment exogeneity of the covariates is not essential for identification, as in many applications, it is more realistic to use only lagged covariates. Thus, part (iii) of Assumption (ref) is relevant only in cases where post-treatment values of some covariates are needed to control for post-treatment confounders and these covariates can be considered as exogenous Brodersen:Gallusser:Koehler:Remy:Scott:2015.\footnote{Although motivating covariates’ choice is often sufficient, it can be also verified by testing for the presence of treatment effects on each covariate: those that are significantly impacted by the intervention must be removed from the model.}
The next identification assumption pertains to the potential outcomes model, similar to Carvalho:Masini:Medeiros:2018 and Masini:Medeiros:2021, but instead of using control units, we rely on a dynamic specification based on lagged outcomes.
Notice that Equation ((ref)) assumes poolability, i.e., homogeneous intercept and slope for all units, meaning that once we account for a highly predictive set of covariates and include temporal dynamics through lagged outcomes, the residual component is essentially random noise.\footnote{For panel data studies with small T and large N, it is usual to pool the observations baltagi2008econometric.} We posit this assumption can be met in practice by using a subset of very predictive covariates out of a larger initial set selected ex ante on the basis of domain knowledge. This approach ensures that major sources of heterogeneity among units are effectively addressed. In addition, we can provide empirical support for this assumption by examining the fit of the pooled model during the pre-intervention period. A good fit, indicated, for instance, by low Mean Squared Error (MSE), would suggest that the pooled model adequately captures the variations in the data. In cases where the researcher also wants to account for time-invariant, group-level unobserved heterogeneity (e.g., units’ fixed effects), the covariate set could be augmented by adopting recent approaches based on sufficient representations for categorical variables Johannemann:Hadad:Athey:Wager:2019.
We remark that $h(\cdot)$ can be either a linear or a non-linear function. In the linear case, no additional assumptions are needed for the identification of causal effects. In the non-linear case, no further assumptions are required if the effect is defined only a single step ahead, as in our empirical application. Conversely, identification of non-linear multi-step-ahead causal effects is notoriously difficult and has not yet been fully explored in the literature (for time series in non-causal settings, see Chen:Yang:Hafner:2004). Therefore, our paper also contributes to the formal identification of causal effects in this challenging context. Specifically, when $h(\cdot)$ is non-linear and the focus is on multi-step-ahead causal effects, we propose modeling the potential outcomes in the post-intervention period as follows.
Notice that Equation ((ref)) does not impose any particularly restrictive assumption on the potential outcomes. Indeed, it merely formalizes the dependence structure between potential outcomes and their conditional expectations already implied by Equation ((ref)), e.g., for $k = 2$ and $p = q = 0$, we have $\operatorname{Y}_{i,t_0+2}(0) = h(\operatorname{Y}_{i,t_0+1|t_0}(0) + \epsilon_{i,t_0+1}, \mathbf{X}_{i,t_0+2}) + \epsilon_{i,t_0+2} = g_2(\operatorname{Y}_{t_0+1|t_0}(0), \operatorname{Y}_{i,t_0}, \mathbf{X}_{i,t_0+2}, \mathbf{X}_{i,t_0+1}) + \epsilon_{i,t_0+2}$.
The panel CV procedure---thoroughly described in Section (ref)---which we propose to estimate the pre-intervention potential outcome model in Equation ((ref)), will reveal whether the winner of the horse-race is a linear or non-linear model. If a linear model is selected, Equation ((ref)) remains valid for all $t = 1, \dots, T$. Otherwise, we can use the non-linear specification in Equation ((ref)) to identify and estimate post-intervention counterfactual outcomes after $t_0+1$.
Finally, the potential outcome model defined by Equation ((ref)) does not assume stationarity of the outcome variable. In similar contexts, studies often adopt a weaker assumption of conditional stationarity DHaultfoeuille:Hoderlein:Sasaki:2023, Hoderlein:White:2012. Formally, let $f_{\operatorname{Y}_{i,t}(0)|\operatorname{Y}_{i,t-1}^{(p)}(0),\mathbf{X}_{i,t}^{(q)}}$ be the conditional density of the potential outcomes under control. Conditional stationarity assumes that this density remains unchanged one step ahead, i.e., $f_{\operatorname{Y}_{i,t}(0)|\operatorname{Y}_{i,t-1}^{(p)}(0),\mathbf{X}_{i,t}^{(q)}} = f_{\operatorname{Y}_{i,t+1}(0)|\operatorname{Y}_{i,t}^{(p)}(0),\mathbf{X}_{i,t+1}^{(q)}}$. In other words, its functional form, which matches the form of the error term, is time-invariant, although the conditional mean can still be time-varying.\footnote{Consider this simple (non-linear) example: $Y_t = \sin{Y_{t-1}^2} + \varepsilon_t$ where $\varepsilon_t \sim f(0,\sigma^2_{\epsilon})$. Under conditional stationarity, we have that $Y_t|Y_{t-1} \sim f(\sin{Y_{t-1}^2}, \sigma^2_{\epsilon})$, $Y_{t+1}|Y_{t} \sim f(\sin{Y_{t}^2}, \sigma^2_{\epsilon})$ and so on. Notice that while the conditional mean is time-varying, the form of the distribution $f(\cdot)$ remains the same. More generally, in our setting of interest, a short panel with small T and large N, non-stationarity is not an issue of particular concern baltagi2008econometric. Nevertheless, as we will show below, the simulation results reported in Supplemental Appendix (ref) demonstrate that even in the case of an explosive autoregressive coefficient (i.e., when the coefficient of the past lag of the outcome $Y_{t-1}$ is set at $\phi = 1.2$), both the bias and the coverage rates of the MLCM are in line with the results obtained for stationary processes (i.e., when $\phi = 0.8$).} This means that the error term would be distributed in the same way during the post-intervention period. In our case, as we explain in Section (ref), we target the causal effect at the individual level and specifically its expected value conditional on past information. Therefore, for identification, estimation, and inference, we do not require the error distribution to remain unchanged after the intervention. Instead, it is sufficient that the error term remains uncorrelated and has a mean of zero, a condition already embedded in Assumptions (ref) and (ref). In other words, we only require that, in the absence of the policy, the potential outcomes model would continue to be correctly specified. This implies the absence of other unforecastable shocks or policies affecting the outcome of interest in the post-intervention window. While it is not possible to explicitly test for this, we can assess how well the model fits the pre-intervention data via the rigorous panel CV procedure (see Section (ref)).
We introduce novel estimands that could be of general interest for observational studies based on panel data in the absence of control units. The causal estimands discussed here pertain to Case (i) defined above, namely, when the treatment affects all, or most, of the available units, simultaneously. In Supplemental Appendix (ref), we also define causal estimands for the Case (ii) scenario with violations of the no-interference assumption.
We begin with the unconditional individual causal effect estimand.
Note that this estimand arises naturally as a consequence of Assumption ((ref)). However, identification, estimation, and inference of this causal effect would require stringent assumptions (e.g., strong stationarity of the potential outcome distribution). Therefore, our focus is on its conditional mean.
This estimand emphasizes that, because the potential outcomes are temporally related, the causal effect of the policy should be defined conditionally on past information. Furthermore, in many situations, the full distribution of $\Delta_{i,t_0+k}$ is not of interest unless we specifically target causal effects at the median or other quantiles. Typically, the mean effect is the primary concern, making $\tau_{i,t_0+k}$ a more suitable summary measure of policy effects. The following theorem establishes the identification conditions for the individual effect $\tau_{i,t_0+k}$.
The following estimands are unit-averages of the individual effects defined above. For the reasons previously stated, we focus on averaging the conditional individual effects, which can be considered the panel analog of the ATE. We also frame them from a finite-sample perspective, as the MLCM is designed for statistical and econometric settings where the treatment affects, either directly or indirectly, the entire population under study.\footnote{For further details on the difference between the finite-sample and superpopulation perspectives in causal inference, see Imbens:Rubin:2015. In observational panel settings, only Rambachan:Roth:2023 defined causal estimands when the entire population is observed, but their approach relies on control units. }
Note that in this benchmark scenario (Case (i) where only treated units are available), the term $\tau_{t_0+k}$ corresponds to the Average Treatment effect on the Treated (ATT).
The subsequent estimand measures the average individual effect within a subset of units with the same values of selected covariates. This is commonly referred to as the CATE, and it is of particular interest when there is reason to believe that the intervention has produced heterogeneous effects on different subpopulations of units as defined by their distinct characteristics. Following existing literature on CATE in high-dimensional settings with a mix of discrete and continuous covariates Fan:Hsu:Lieli:Zhang:2022, Knaus:Lechner:Strittmatter:2021, Chernozhukov:Demirer:Duflo:Fernandez:2018, we focus on its low-dimensional summary called “group-average treatment effect”.
Note that the set of candidate conditioning variables $G_{i,t}$ used to determine CATEs does not necessarily have to match those utilized for counterfactual forecasting. This is because the variables that predict outcomes can, and often do, differ at least partially from those that predict treatment effect heterogeneity. We also remark that in a general panel setting with more than one post-treatment periods, the above definitions imply the existence of vectors representing estimated average effects, i.e., $ (\tau_{t_0+1}, \dots, \tau_T )$ and $( \tau_{t_0+1}(g), \dots, \tau_T (g))$. Furthermore, researchers are sometimes also interested in temporal aggregations of such effects.
We now introduce the following Theorem.
We now introduce estimators of the causal quantities defined in Section (ref). Note that imputing future counterfactual outcomes based on past information and covariates is essentially a forecasting problem. Therefore, in the following definitions, the estimator $\widehat{\operatorname{Y}}_{t_0+k|t_0}$ represents the $k$-step-ahead forecast based on the model trained in the pre-intervention period. We first define an estimator for the individual effect at time $t_0+1$.
A detailed description of the estimation algorithm for $h(\cdot)$ will be given in Section (ref). Although our empirical application only evaluates causal effects at time $t_0+1$, one might also be interested in estimating effects over multiple post-intervention periods, as in our simulations and replication study. There are two main strategies for generating multi-step-ahead forecasts Chevillon:2007, Forecasting:2022: the recursive approach, where each forecast is defined and estimated using previous forecasts, and the direct approach, where forecasts are produced by estimating separate models for each forecast horizon. In the case of linear models, the recursive approach yields more efficient parameter estimates, especially when the model is correctly specified and for a long forecasting horizon Pesaran:Pick:Timmermann:2011. However, with non-linear models, the recursive approach is known to be asymptotically biased, making the direct approach often preferable Bontempi:Taieb:2011. One criticism of the direct forecasting approach is that it assumes independence between the forecasts, as each one is based on a separate model. Hybrid forecasts, instead, represent a combination of the recursive and direct approaches, effectively integrating the strengths of both strategies Forecasting:2022. Like the direct approach, hybrid strategies require estimating separate models for each forecasting horizon. However, similar to the recursive approach, each model also includes the direct forecast from the previous step. This method offers full flexibility by allowing different models for each time horizon, while still accounting for the dependency between direct forecasts. This is exactly what Equation ((ref)) captures: each $\operatorname{Y}_{i,t_0+k-1|t_0}(0), \dots, \operatorname{Y}_{i,t_0+1|t_0}(0)$ is a direct forecast (i.e., the expected potential outcome up to $k-1$ steps ahead from $t_0$), and $\operatorname{Y}_{i,t_0+k}(0)$ is then modeled as a flexible function of these direct forecasts, along with observed past lags and covariates.
Therefore, depending on whether the potential outcome model is linear or non-linear, the next definition outlines two different estimators for multi-step-ahead evaluations.
A detailed description of the estimation algorithm for $g_k(\cdot)$ will be given in Section (ref). To clarify how multi-step-ahead forecasts are computed and to highlight the differences between the linear and non-linear case, we provide the following example.
Building on Equations ((ref)-(ref)), the next definition summarizes the ATE and CATE estimators under the MLCM.
Inference on the estimated causal effects defined above is conducted using block-bootstrap. Refer to Supplemental Appendix (ref) for a detailed description of the block-bootstrap algorithms used to derive confidence intervals for the ATE and CATEs. Finally, we stress that deriving the theoretical properties of the above-defined estimators would deviate from the intended nature of the methodology, since the MLCM can be implemented with any supervised ML algorithm, and it is unknown in advance which one will outperform the others.\footnote{To use Breiman's words: “Nowhere is it written on a stone tablet what kind of model should be used to solve problems involving data." breiman2001statistical.} Instead, in Section (ref), we perform a simulation study showing that the MLCM can achieve forecast unbiasedness and is able to detect causal effects even in challenging econometric environments with non-linearities and irrelevant covariates included in the model specification.
Supervised ML techniques primarily aim to minimize the out-of-sample prediction error, generalizing well on unseen data. The degree of flexibility is the result of a trade-off: increased flexibility can enhance in-sample fit but may diminish out-of-sample fit due to overfitting. ML algorithms tackle this trade-off by relying on empirical tuning to choose the optimal level of complexity.
The standard ML approach is to randomly split the sample into two sets, containing, for instance, 2/3 and 1/3 of observations. One then uses the first set to train ML algorithms (training set) and the second to test them (testing set). This introduces a “firewall” principle: none of the data involved in generating the prediction function is used to evaluate it Mullainathan:Spiess:2017. The out-of-sample performance of the model on the unseen (held-out) data of the testing set can be considered a reliable measure of the “true” performance on future data. In order to solve the bias-variance trade-off and prevent overfitting, one can rely on automatic tuning using tools such as random k-fold CV on the training sample to select the best-performing values of the tuning parameters in terms of an a priori defined metric, such as the MSE.
We depart from this standard ML routine and reorient it towards the counterfactual forecasting goal. First, we do not randomly split the data, but we train, tune, and evaluate the models only on the pre-treatment data (Design Stage); then, we use the final selected model to forecast counterfactual post-treatment outcomes. In this forecasting perspective, the unseen data on which the ML models must generalize well are not the outcomes of different units (as in typical out-of-sample prediction tasks), but future observations of the outcome for the same set of units employed to train the models. Stated differently, the aim is to make ML models learn as best as possible the pre-intervention outcome trend for each treated unit, so as to predict the best possible counterfactual outcome under the no-treatment scenario. The key implication of this unconventional ML setup is the shift in focus: the primary concern becomes ensuring unbiasedness in forecasts, rather than focusing only on forecast accuracy.
Second, and related, we do not carry out hyperparameter tuning and model selection with random k-fold CV. The panel dimension creates an additional challenge regarding how to implement CV, because standard CV does not account for the temporal structure of the data arkhangelsky2024causal. To address this challenge, we introduce a resampling technique suited for forecasting tasks on panel data—panel CV—which is described below.
Finally, a key concern in ML regards the trade-off between accuracy and interpretability. Such a trade-off is relevant when ML is used for tasks that take into consideration transparency aspects. In principle, our method can be used with any supervised ML routine, including black-box techniques like deep neural networks. However, maintaining some degree of transparency can be important, because higher interpretability of the estimated counterfactuals bolsters the credibility of the proposed approach Abadie:2021. In practice, users should decide on the basis of a comparative assessment across a mix of models characterized by different layers of complexity by balancing any improvements in performance from complex models against the loss of interpretability associated with their use.
The MLCM is rooted in the idea of developing a ML forecasting model that can closely reproduce the outcome trajectories of treated units in the pre-intervention window, so that any post-intervention divergence between the observed and forecasted outcomes can be attributed to the treatment under the identification assumptions. To achieve this, the implementation of the MLCM requires ten empirical steps, which are divided into the Design Stage and the Analysis Stage. The full process is summarized in Box (ref) and described in detail below.
In the Design Stage, the first step involves the selection of supervised ML algorithms that will play the horse-race of performance testing on the pre-treatment data. This competition among multiple methods is conceptually analogous to the way the best ML learner is selected in the double/debiased ML method chernozhukov2018double.\footnote{The MLCM can leverage any supervised ML algorithm, including deep learning techniques, some of which, such as RNNs and transformers, are explicitly designed to learn temporal dynamics; however, they cannot be employed in the panel settings we focus on (small T, large N), as they require a large number of time periods to be applicable. It is also possible to stack several different ML algorithms and form complex ensemble learners, but this would come at the cost of a substantial loss in transparency.}
In the second step, we recommend deploying some tweaks involving feature pre-selection and engineering drawing from consolidated practices in applied predictive modeling Kuhn:Johnson:2013. A crucial aspect of this pre-processing is the strategy for selecting predictors. There are two main schools of thought: those who advocate for purely data-driven selection argue that one should build a dataset as large as possible, and then let the algorithm autonomously decide which variables matter for the prediction task. Others stress the importance of subject matter knowledge: the researcher should select ex ante the relevant predictors, and then feed only those to the algorithm. The rationale is that subject matter knowledge can separate meaningful from irrelevant information, eliminating detrimental noise and enhancing the underlying signal Kuhn:Johnson:2013. We propose a hybrid approach: build a large initial dataset on the basis of domain knowledge, then adopt preliminary and data-driven variable selection criteria to drop non-informative predictors.
As outlined in the causal framework, counterfactual forecasting is carried out by using an information set mainly comprising lagged values of outcomes and covariates. However, determining the optimal number of lags to include is an empirical question. We recommend including at least two lagged values of both outcomes and covariates, and then using a data-driven approach to select a subsample of the most relevant features. The latter step allows for a reduction of the risk of overfitting and degradation of performances as well as facilitating interpretability. Finally, feature engineering and data pre-processing matter too because how the predictors enter into the model is also important Kuhn:Johnson:2013.
To carry out model selection and validation on the pre-treatment sample, we propose a panel CV approach. Using an alternative CV procedure is necessary since ML methods do not natively handle longitudinal data and are designed for predicting rather than forecasting.
Our panel CV approach adapts time series CV based on expanding training windows to a panel setting.\footnote{Here we use an expanding window (as we are in a short-panel setting), but the procedure can also be implemented using a rolling window approach, which might be more appropriate in specific cases (e.g., with long panels or non-stationarity of the outcome).} The intuition is provided in Figure (ref) and constitutes an adaptation from Hyndman:Athanasopoulos:2021.
We start from a set of candidate algorithms $\mathcal{H} = (h_1, \dots, h_J)$, each having its own set of parameters $\theta^{(j)}$ and hyperparameters $\gamma^{(j)}$. For example, in the empirical application, we use four learners: LASSO, PLS, random forest and stochastic gradient boosting. For each learner, we carry out panel CV as detailed in Algorithm (ref). In short, the ML algorithms are trained repeatedly using an expanding window approach, with hyperparameters tuned at each step to minimize forecast errors. This sequential procedure ensures that, for each unit, there are no “future” observations in the training set and no “past” observations in the validation set. Note that the same procedure can be applied to forecast post-treatment values of important covariates affected by the intervention. We can assume that such covariates follow Equation ((ref)) and use the ML algorithms and the panel CV routine to forecast counterfactual values of $\widehat{\mathbf{X}}_{i, t_0+k|t_0}(0)$, which we can then incorporate into the prediction of $\widehat{\operatorname{Y}}_{i,t_0+k|t_0}$ as in liu2020forecasting.
To provide evidence about internal validity, we suggest running diagnostic checks (e.g., showing that the distribution of the pre-treatment forecasting errors is approximately Gaussian and centered around zero) as well as placebo tests, which have become a key device for assessing the credibility of research designs in observational settings eggers2024placebo. Following Liu:Wang:Xu:2022, panel placebo tests can be implemented by hiding one or more periods of observations right before the onset of the treatment and using a model trained on the rest of the pre-treatment periods to predict the untreated outcomes of the held-out period(s). If the identifying assumptions are valid, the differences between the observed and forecast outcomes in those periods should be close to zero.\footnote{If large discrepancies, i.e., forecasting errors, are found at this stage, this would indicate pre-intervention shocks or policy changes occurring between pre-treatment periods. These factors should be accounted for in the model to bridge the gap between observed and forecasted pre-treatment outcomes.} Given that we estimate unit-level treatment effects, we are also able to test whether most unit-level placebo differences are close to zero. The Design Stage thus ends with a battery of performance, diagnostic, and placebo tests.
The Analysis Stage starts with final model selection and training: on the basis of the comparative performance assessment in the Design Stage, pick the best-performing model, re-train it on the full pre-treatment sample (using the hyperparameter values obtained in the Design Stage), then use it to forecast counterfactual outcomes in the post-intervention period. Once that is done, treatment effects for each unit are given by the difference between the post-treatment observed data and the corresponding ML-generated counterfactual forecasts. The ATE is the average of the individual effects.
When dealing with a single post-intervention period, step 7 of the Analysis Stage is straightforward. However, when performing a multi-step-ahead forecast using a data-driven ML routine, one must take special care to avoid including post-treatment outcomes in the prediction. Multi-step-ahead forecast becomes more challenging when the selected ML algorithm is non-linear, as we cannot simply plug-in the forecasted value $\widehat{\operatorname{Y}}_{t_0+1|t_0}(0)$ into the already estimated model. As explained in Section (ref), one solution is to re-estimate the potential outcome model for values of $k \geq 2$. Importantly, at the end of the panel CV procedure, we know the selected ML algorithm and can adjust our forecasting strategy accordingly. Specifically, we propose using Equation ((ref)) to model post-intervention potential outcomes, as this approach has proven effective in identifying non-linear, multi-step-ahead causal effects, and estimate them with Equation ((ref)). From a practical standpoint, this requires repeating the panel CV for each step of the forecasting horizon. The detailed procedure is described in Algorithm (ref).
During the counterfactual forecasting step, especially when opting for more complex and less transparent techniques, users should consider leveraging innovations in Interpretable Machine Learning and Explainable Artificial Intelligence molnar2020interpretable, such as model-agnostic Shapley Additive Explanations (SHAP) lundberg2017unified, to enhance explainability and transparency of the counterfactual building process Abadie:2021. \frenchspacing
After estimating individual effects, data-driven CATEs can be computed via a regression tree analysis. Specifically, this approach uses the estimated treatment effects as the outcome variable, regresses them on many potentially associated variables, and lets the algorithm pick the main predictors and their critical thresholds. The resulting tree reports the group-average treatment effects for all units in each terminal node. As for causal trees and forests Athey:Imbens:2016, Wager:Athey:2018, this data-driven search for heterogeneity of causal effects removes a major degree of discretion because the researcher can only select the set of covariates that can be used by the tree to build the subgroups. However, the two approaches differ regarding both purpose and implementation: our data-driven technique is a post-estimation approach aimed at automatically recovering and visualizing the relevant heterogeneity dimensions, while causal trees and forests are counterfactual methods for the direct estimation of heterogeneous treatment effects, which they achieve by leveraging control units. Moreover, we are not interested in the out-of-sample performance of the regression tree, but in retrieving CATEs for the entire population of interest. To this end, there is no need to split the sample into training and testing sets or to prune the tree by adjusting the complexity parameter. However, it may be necessary to set a pre-specified minimum node size to preserve the interpretability of the resulting tree.
Finally, standard errors and confidence intervals for the ATE and CATEs are estimated through the block-bootstrap approach described in Supplemental Appendix (ref).\footnote{In staggered adoption settings, uncertainty quantification with block-bootstrapping can be extended to other estimands, such as group-time average treatment effects. See the replication exercise in Supplemental Appendix (ref).}
To investigate the performance of the MLCM in detecting average treatment effects in panel datasets, we performed an extensive simulation study using different data-generating processes (linear and non-linear) and various combinations of pre- and post-intervention periods. We generated $1,000$ panel datasets, each consisting of $400$ units and $T = 7$ or $12$ time periods, according to the following two models,
where: the error term $\epsilon_{i,t}$ is generated from a Normal distribution with standard deviation $\sigma_{\epsilon} = 2$, i.e., $\epsilon_{i,t} \sim N(0,2)$, $\phi = 0.8$ is the autoregressive coefficient, $\boldsymbol{\beta} = (\beta^{(1)}, \dots, \beta^{(11)})'$ is a $m \times 1$ vector of coefficients and $\mathbf{X}_{i,t-1} = \left(\operatorname{X}_{i,t-1}^{(1)}, \dots, \operatorname{X}_{i,t-1}^{(11)} \right)$ is a $1 \times m$ vector of $11$ predictors, both continuous and categorical, also containing interaction terms and correlated regressors. We allowed the covariates to vary in time and across units in the dataset by adding, respectively, a random term $\nu_t$ (which, as will become clear later, varies for different covariates) and a random term $u_i \sim N(1,1)$. The latter term adds variability (and thus, heterogeneity) between the units. In particular, the covariates are generated as follows: $\operatorname{X}_{i,t}^{(1)} = 0.1 t+u_i+\nu_t^{(1)} , \nu_t^{(1)} \sim N(0,1)$, $\operatorname{X}_{i,t}^{(2)} = 0.1 t+u_i+\nu_t^{(2)} , \nu_t^{(2)} \sim N(0,0.2)$, $\operatorname{X}_{i,t}^{(3,4,5)} = u_i+\nu_t^{(3,4,5)} , \nu_t^{(3,4,5)} \sim MVN(0, \Sigma)$, $\operatorname{X}_{i,t}^{(6)} = u_i-\nu_t^{(6)}, \nu_t^{(6)} \sim N(0,1)$, $\operatorname{X}_{i,t}^{(7)} = \left(0.1 t + \nu_t^{(1)} \right)^2 + u_i + \nu_t^{(7)}, \nu_t^{(7)} \sim N(0,0.2)$, $\operatorname{X}_{i,t}^{(8)} \in \{0,1\}$, $\operatorname{X}_{i,t}^{(9)} \in \{1,2,3\}$, $\operatorname{X}_{i,t}^{(10)} = \operatorname{X}_{i,t}^{(3)} \cdot \operatorname{X}_{i,t}^{(9)}$, $\operatorname{X}_{i,t}^{(11)} = \operatorname{X}_{i,t}^{(2)} \cdot \operatorname{X}_{i,t}^{(8)}$.\footnote{The variance-covariance matrix of the multivariate normal distribution used to generate covariates 3--5 is set to $\Sigma =
$. Notice that, to put the ML algorithms under further stress, the covariance between $\operatorname{X}_{i,t}^{(3)}$ and $\operatorname{X}_{i,t}^{(5)}$ is $0.7$ so the two variables are highly correlated but $\operatorname{X}_{i,t}^{(3)}$ is ten times more important (in terms of coefficient) than $\operatorname{X}_{i,t}^{(5)}$.} These choices regarding the data-generating process for the covariates included in our simulation study are driven by the goal of mirroring real-world empirical settings where the methodology might be applied. In addition, to better reflect a typical real-world scenario where only a subset of covariates is relevant, we set certain coefficients to zero: specifically, $\beta^{(1)}$, $\beta^{(8)}$, and $\beta^{(9)}$. The remaining coefficients are generated as follows: $\beta^{(2)} = \beta^{(6)} = \beta^{(10)} = 2$, $\beta^{(3)} = \beta^{(7)} = 1$, $\beta^{(4)} = 2.5$, $\beta^{(5)} = 0.1$ and $\beta^{(11)} = 1.5$. We also remark that, for computational reasons, in this simulation study we focus on a low-dimensional set of covariates. However, the MLCM is particularly well-suited for data-rich environments and sparse settings with many more covariates. Therefore, the simulation results provided here should be interpreted as conservative evidence regarding the performance of the MLCM. Finally, we assume exogeneity of the predictors in both the pre- and post-intervention periods (third part of Assumption (ref)). Table (ref) below provides an overview of the generated datasets.
At $t_0+1$ (corresponding to T = $5$ in Table (ref)), we included a fictional intervention that increases the outcome for each unit by $2$ standard deviations at time $t_0+1$, $1.5$ standard deviations at time $t_0+2$, and $1$ standard deviation at time $t_0+3$. In other words, the intervention has a decreasing effect over time. We made this choice because adding a unit-specific component in the covariates generates heterogeneity; therefore, the scale of each $\operatorname{Y}_i$ varies across the $i$’s. We measured the performance of the MLCM in terms of both the bias of the estimated effect from the true impact and the interval coverage. Also note that, under this setup, the number of pre-intervention periods is $t_0 = 4$ when $T = 7$ and $t_0 = 9$ when $T = 12$. As estimators, we employ four different ML algorithms - LASSO, Partial Least Squares, random forest, and stochastic gradient boosting.
The results are summarized in Table (ref) below and show that, across all post-intervention horizons, the MLCM achieves a very low bias. This holds true for both linear and non-linear model specifications. The bias tends to decrease when the number of pre-intervention time periods increases, as more information is present in the data. Nevertheless, the bias at the shortest time period is still very low, which reinforces our belief that the MLCM can be effectively used for short panels. For instance, under the linear specification and $t_0 = 4$, the relative bias increases from $0.2\%$ measured at $t_0+1$ to $0.9\%$ measured at $t_0+3$. However, by increasing the number of pre-intervention time points, the bias drops at $0.4\%$ even at $t_0+3$. We also observe that the relative bias under a linear model specification is much lower than that under the non-linear specification, which was expected since non-linearities are typically more challenging to detect. In any case, the bias under a non-linear model specification is still low (the maximum bias is $4.9\%$ and is measured in the most challenging scenario, i.e., at the third time horizon when the pre-intervention series is very short, $t_0 = 4$).
Similar observations can be made for the interval coverage, which was estimated based on $1,000$ block-bootstrap iterations. At the first time horizon after the intervention, the interval coverage is close or equal to the nominal $95\%$ level for both linear and non-linear model specifications. Then, it tends to slightly decrease as we move further away from the intervention, but remains close to the nominal level for both linear and non-linear processes.
Overall, these numerical studies demonstrate that the MLCM can achieve forecast unbiasedness and that the coverage of the estimator is high even when considering short pre-intervention windows. The simulation results also suggest that when we are interested in the estimation of causal effects $k$-step ahead from the treatment, having a greater number of pre-intervention periods further mitigates the bias, both under linear and non-linear data-generating processes.
In Supplemental Appendix (ref), we provide many additional results that show that these key findings remain largely unchanged if we consider alternative data-generating processes for both the outcome variable and the covariates, and if we reduce the number of available units. Notably, when $\phi$ is set to $1.2$ (cf. Table (ref)), rather than to $0.8$ as in Table (ref), i.e., when moving from a stationary to a non-stationary data-generating process with an explosive autoregressive coefficient for the lag of the outcome variable, the MLCM exhibits very similar performance.
The main empirical illustration provided below is an original study on the educational effects of the COVID-19 crisis, thus focusing on a simultaneous treatment affecting all available units. In Supplemental Appendix (ref), to demonstrate how to leverage the MLCM in settings where there are untreated units, but there may be violations of the no-interference assumption, we revisit the application on minimum wage with staggered treatment adoption in callaway2021difference without using their control group. Both empirical illustrations rely on short panel datasets.
In the aftermath of the COVID-19 pandemic, the inequality legacy of this unprecedented crisis has emerged as an issue of great policy relevance. The available evidence on income inequality documents a decrease driven by short-run and temporary government compensation policies, which were mostly targeted at the poorest segments of the population Stantcheva:2022. Concerning education, instead, recent micro-level evidence Agostinelli:Doepke:Sorrenti:Zilibotti:2022, Carlana:Ferrara:Lopez:2023 shows that school closures had a large, persistent, and unequal effect on learning. The educational gaps caused by the school closures may, in turn, permanently affect the lifetime income possibilities of the current generation of students, with vast repercussions on future inequalities Werner:Woessmann:2023. Therefore, it is through the human capital channel that the inequality effects of the pandemic may eventually appear in the medium and long run. It follows that to anticipate longer-term consequences on income inequality and territorial disparities, it is necessary to gauge the magnitude of educational losses and their distribution across areas.
To our knowledge, there is no granular evidence on the geography of the education effects of the COVID-19 crisis. This is mainly due to identification challenges caused by the sudden spread of the pandemic across the world, which resulted in the absence of a suitable untreated group. Italy is an important case study, as it ranks among the hardest-hit countries and was the first Western country to impose a strict nationwide lockdown. During the lockdown that began on March 9, 2020, schools across the entire national territory were closed and remained so until the end of the school year. This resulted in Italy ranking among the OECD countries with the highest number of weeks of school closures and distance learning Battisti:Maggio:2023. In this setting, we leverage the MLCM with non-staggered adoption, given that the treatment---the COVID-19 shock---simultaneously affected all units.\footnote{We refer to the treatment variable as COVID-19 “shock”, which encompasses both the pandemic and the containment policies, as we aim to capture the total effect of the pandemic, i.e., its direct and indirect effects on educational outcomes. In line with previous literature, we consider the restrictive measures, including the lockdown and the school closures, as a manifestation of the pandemic, not as a separate treatment. COVID-related learning losses can originate from many pandemic-induced channels, such as school closures and distance learning, prolonged absence from school due to one or more infections, absence due to parents’ concerns about their children’s health, and many other transmission mechanisms. At the aggregate level, only the overall impact of the pandemic remains discernible.}
We employ LLM yearly data covering the period from 2013 to 2020. We cover all Italian LLMs except some of the smallest ones for which the education data are unavailable due to privacy protection, for a total of $579$ LLMs ($95\%$ of Italian LLMs and over $99\%$ of the Italian population). The dependent variable is the standardized math test score of fifth-grade students referred to the school year started in 2020. Following the approach by Carlana:Ferrara:Lopez:2023, the math score is standardized with respect to its pre-pandemic mean. This way, the detected treatment effects can be interpreted as deviations from the pre-COVID baseline. We focus on younger students for two reasons: i) learning disruptions earlier in life typically have longer-lasting and more severe effects, so younger children may have been more heavily impacted Stantcheva:2022; ii) primary schools were excluded from the heterogeneous school closures involved in the tier system of regional restrictions implemented by the Italian government since the onset of the second wave (November 2020), so the treatment is the same for all units (Assumption (ref)). The policy relevance of the math outcome is illustrated by the fact that the most recent report by the OECD Program for International Student Assessment (PISA), published in December 2023, reported that math scores dropped globally in the last few years, making headlines around the world.\footnote{See \href{https://www.nytimes.com/2023/12/05/us/math-scores-pandemic-pisa.html}{\textcolor{blue}{here}} for coverage of the news from the New York Times.}
The initial pre-treatment information set includes over $150$ variables (see Table (ref) in Supplemental Appendix (ref) for a detailed description of the variables). In this set of covariates, we included the first three lags of all the predictors as covariates and three lags of the outcome variable. This implies that we collapse the original 2013–2020 dataset into a dataset covering the period 2017–2020.
The implementation process for this application is reported in Box (ref), while the outcomes of the empirical analysis are reported in Section (ref).
We use a data-driven approach to restrict the information set. More specifically, we follow Athey:Wager:2019 and Basu:Kumbier:Brown:Yu:2018 and apply a pilot random forest on the pre-treatment data to pick up a subset of the most important predictors according to the importance ranking produced by the forest.\footnote{For this preliminary operation, we use default hyperparameter settings that typically perform well with random forests Athey:Wager:2019. See Figure (ref) for the variable importance ranking of the preliminary forest.} In order to select the precise number of relevant predictors in a data-driven manner, we include this number as an additional parameter in the subsequent panel CV routine which we employ to tune the hyperparameters of all the selected ML algorithms (LASSO, Partial Least Squares, stochastic gradient boosting, and random forest ).\footnote{For all algorithms, we tune the subset of most predictive variables (10, 20 or 30) to be included in the analysis according to the preliminary forest. For LASSO, we tune the penalty parameter (the values from 0.1 to 0.9 with an increase of 0.1 in each run); for PLS, we tune the number of components (all possible values from 1 to the number of variables). For boosting, we tune the number of trees (1,000 or 2,000), the maximum depth of each tree (1 or 2), the minimum number of observations in terminal nodes (5 or 10) and the learning rate (0.001, 0.002 or 0.005); for random forest, we tune the number of variables randomly sampled as candidates at each split (1/2, 1/3 or 1/4 of the number of variables), and use a fixed number of 1,000 trees. All candidate hyperparameter values are tested on each testing set.}
Note that in this application, due to the ubiquitous nature of this shock, we do not invoke the third part of Assumption (ref) (exogeneity of post-treatment covariates) and only include pre-treatment values of time-varying variables and time-invariant predictors in the information set. Assumption (ref) is satisfied because the pandemic was completely unanticipated. Regarding Assumption (ref), given the unprecedented disruptions brought about by the pandemic and based on institutional and subject matter knowledge, we deem it plausible that any divergence from the forecasted counterfactual in educational outcomes can be attributed exclusively to the COVID-19 shock, and not to other concomitant shocks or policies.\footnote{There were no educational reforms at the primary school level that could have influenced the standardized math test scores of fifth-grade students in the period under scrutiny.} Finally, Assumption (ref) is not needed for identification, as we only estimate one-step-ahead causal effects here.
Table (ref) below reports the average performance—in terms of panel CV MSE—of the four selected ML algorithms in forecasting the standardized math test score of fifth-grade students across the pre-treatment test sets (2018 and 2019).\footnote{A replication notebook of the main results reported in this section is publicly available \href{https://marclet.github.io/MLCM-Replication-Notebook/}{\textcolor{blue}{here}}.}
The best-performing MLCM version is the one using the random forest with twenty variables and 1/3 of them, i.e., six, as candidates at each split. Overall, the fully non-linear random forest and boosting fare better than the linear ML methods, suggesting the presence of non-linearities in the data-generating process. Moreover, the non-linear methods outperform LASSO and Partial Least Squares in all cases, and their performance is stronger in the presence of very few variables and less deteriorated by the inclusion of variables with poor predictive power (cf. Table (ref) in Supplemental Appendix (ref)). See Table (ref) in Supplemental Appendix (ref) for a detailed description of the 20 variables with the highest importance score attributed by the preliminary forest).
Panel CV on pre-treatment periods implicitly enables in-time placebo analysis along the lines of Bertrand:Duflo:Mullainathan:2004 and Liu:Wang:Xu:2022. In-time placebos are performed on the same pool of treated units where we “fake” that the treatment occurred at time $t_0-1$ and use only information up to $t_0-1$ to forecast the counterfactuals at time $t_0$. Since we know what the real values at time $t_0$ are, we can shift the treatment date artificially backward in time to assess the difference with the main estimates and evaluate the forecasting accuracy and unbiasedness of our estimator. As reported in the table, the placebo ATEs are indistinguishable from zero and statistically insignificant both in 2018 (ATE: 0.0554, with 95% confidence intervals (–0.1182; 0.2765) and 2019 (ATE: 0.0254, with 95% confidence intervals (–0.0273; 0.0728).\footnote{The counterfactual estimates exhibit substantial similarity regardless of whether a preliminary random forest or an alternative preliminary procedure is employed, such as boosting or selecting covariates based solely on the highest absolute correlation coefficients (e.g., the correlation between the counterfactual estimates and those obtained using default hyperparameter values is 0.9985). Detailed results are available upon request.}
In addition, the unit-level placebo map reported in Figure (ref) depicts the temporal average (for the years 2018 and 2019) individual `treatment' effects estimated with the best-performing ML algorithm, random forest, and shows that almost all LLMs exhibit no trace of significant differences between the forecasted and observed math scores in the years before the COVID-19 pandemic. In particular, only 2% of the LLMs report a difference lower than –1 standard deviation.\footnote{Figure (ref) in Supplemental Appendix (ref) presents the placebo maps for the single pre-treatment years (2018 and 2019). In these instances as well, few LLMs exhibit a difference lower than –1 standard deviation between the forecasted and observed math scores prior to the pandemic (4.3% of LLMs in 2018 and 5.7% of LLMs in 2019). Furthermore, the LLMs with extreme values tend to be the smallest ones, which are also the most volatile. }
Finally, Figure (ref) reports a diagnostic test for the performance of the MLCM with random forest. The figure illustrates that the distribution of the placebo forecasting errors is approximately normal and centered around zero, not only for the temporal average but also for the single pre-treatment years. This demonstrates that the ATE estimator is nearly unbiased even when applied on very short panels and without relying on the use of post-treatment information.
Figure (ref) shows the local effects of the pandemic on math scores of fifth-grade students across Italian LLMs and reports the ATE. These estimates come from the best-performing technique of the Design Stage, the MLCM with random forest.\footnote{Figure (ref) in Supplemental Appendix (ref) provides SHAP values for this model for the counterfactual 2020 forecasts.}
The first insight is that the COVID-19 shock engendered major learning losses in Italy. The ATE points to a decrease in the math score of –0.7608 standard deviation (SD) with respect to the pre-pandemic average, and this effect is strongly statistically significant (95% confidence intervals: –0.7976; –0.6884). This result is in line with micro-level studies documenting large negative impacts on educational outcomes in Italy Battisti:Maggio:2023, Carlana:Ferrara:Lopez:2023.
As shown above, some of the smallest LLMs (the most volatile) were imprecisely estimated in the placebo exercise. Since equal weight is attributed to all treated units, the presence of these LLMs might affect the accuracy of the ATE estimate. To address this potential concern, we recommend a sensitivity analysis that re-estimates the ATE following the exclusion of units with the poorest fit in the pre-treatment period. Here, we removed the top 1%, 2%, and 5% of the most imprecisely estimated LLMs in the placebo tests. The resultant ATE estimates were -0.7500, -0.7568, and -0.7341, respectively—each closely aligning with the ATE estimate of -0.7608.
The second key finding is that this average impact masks considerable heterogeneity across the Italian territory. Some LLMs experienced much larger drops in students’ performances. In particular, the MLCM detects learning losses of more than one standard deviation in clusters of local economies mostly located in the South of Italy (Campania, Apulia, Molise, Calabria, Sicily, and Sardinia), and drops in math scores greater than two SD in scattered areas predominantly concentrated in Central and Southern regions.
Figure (ref) presents data-driven CATEs estimated with the regression tree analysis.\footnote{For this analysis, we selected a set of LLM-level variables including, among others, pre-pandemic variables such as the per capita income, the unemployment rate, the Gini index, the share of university graduates, two proxies of access to distance learning that capture the digital divide across areas, and excess deaths registered during the first wave of the pandemic Cerqua:DiStefano:Letta:Miccoli:2021. The full list can be found in Table (ref) in Supplemental Appendix (ref).} The algorithm selects only three predictors to construct the tree. The analysis reveals that the most substantial learning losses (1.22 SD below the pre-COVID baseline) occurred in LLMs characterized by an unemployment rate equal to or above 10.19%, a percentage of university graduates below 26%, and a Gini index equal to or above 0.42. Conversely, the territories with an above-average treatment effect (–0.56 SD) exhibit an unemployment rate below 10.19% and a percentage of university graduates equal to or above 23%. Finally, all the CATEs reported in the tree are statistically significant.
Since areas with higher unemployment, greater inequality, and lower education rates at baseline have experienced disproportionate impacts, the pre-existing gaps across the country have been amplified. If such gaps persist, it is likely that they will act as a catalyst for future economic inequality, leading to widened territorial disparities in the medium and long term. In summary, we can anticipate that, in the absence of counterbalancing policy efforts, the pandemic will exacerbate long-standing inequalities that predate COVID-19.
Identifying causal effects is challenging without a credible control group. We overcome this challenge by proposing a new method based on counterfactual forecasting with machine learning. The MLCM can employ any off-the-shelf ML algorithm to estimate policy-relevant causal parameters in a wide variety of panel settings, including very short panels, without relying on untreated units. To illustrate its potential, we presented numerical studies, a replication of the minimum wage application in callaway2021difference without using untreated units, and an empirical analysis on the effects of the COVID-19 crisis on education in Italy. The companion R package \href{https://github.com/FMenchetti/MachineControl}{\textcolor{blue}{MachineControl}} provides an easy-to-use implementation of the proposed approach. Large shocks, international economic policies, nationwide policy changes, and regional programs engendering substantial interaction between units, are all real-world scenarios where a control group may not exist. In such cases, researchers can harness the MLCM, which complements the existing econometric toolbox for causal inference and policy evaluation.