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.
136,364 characters · 51 sections · 69 citation commands
Causal Forecasting in Panel Data: A Two-Way Synthetic Forecasting Approach
Estimating causal effects in panel data is a central problem in policy evaluation. Existing methods largely address retrospective questions of the form: what would have happened to a target unit under a different intervention during the observed panel? In many applications, however, decision-makers face prospective questions: what will happen to a target unit under an intervention it has not yet experienced, beyond the observed panel? This article develops a framework for answering such causal forecasting questions by integrating the retrospective counterfactual logic of synthetic-controls-based approaches with the extrapolative structure of multivariate time-series forecasting. Building on the latent factor models that justify unit-side regressions in synthetic controls, we impose low-rank temporal structure on the latent time factors to identify prospective causal forecast estimands. We operationalize this strategy through the Two-Way Synthetic Forecasting estimator, or TWSF, which learns cross-unit relationships from pre-treatment outcomes and combines them with a time-series model learned from treated donor trajectories under the intervention of interest. Under suitable conditions, we establish finite-sample forecasting error bounds that imply pointwise consistency and introduce an orthogonalized correction that yields asymptotic normality and thus enables pointwise inference. We extend the framework to fixed multi-step forecasting horizons through both direct and recursive procedures, each of which inherits analogous pointwise guarantees. We corroborate the theory with simulation studies and illustrate the practical utility of TWSF by studying the public-health impact of opening NFL stadiums during the 2020 season.
Keywords: synthetic controls; synthetic interventions; multivariate singular spectrum analysis; debiased machine learning; pointwise inference
In September 2020, at the height of the COVID-19 pandemic, the National Football League (NFL) season was set to begin. Each team faced a difficult and highly localized policy decision: should they admit fans into their stadiums, or should they keep their stadiums closed? These decisions were made in coordination with local communities and public-health authorities, and they reflected local guidelines, disease prevalence, risk tolerance, and operational constraints. Teams that admitted fans adopted a range of mitigation measures, including staggered entry, masking requirements, health questionnaires, staff temperature checks, compliance officers, modified concessions, social distancing in seating and lines, mobile ticketing, enhanced cleaning protocols, public-health messaging, and capacity limits mac21a. Ultimately, some teams, such as the Atlanta Falcons, allowed fans to attend games, while others, such as the Seattle Seahawks, remained closed.
From a research perspective, this setting provided a collection of natural experiments and raised several causal questions about the public-health consequences of opening stadiums. The first question is perhaps the most familiar: {\em what would have happened to a city that opened its stadium had it instead kept the stadium closed?} For example, what would COVID-19 case rates in Atlanta have looked like had the Falcons not admitted fans?
To formalize this question, let $Y_{it}(d) \in \Rb$ denote the potential case rate for city $i$ at time $t$ under policy $d$, where $d=0$ denotes keeping the stadium closed and $d=1$ denotes opening the stadium to fans. Let $Y_{it}$ denote the observed case rate. Prior to the start of the NFL season on September 10, 2020, all cities were operated under $d=0$. We refer to this interval as the {\em pre-treatment} period and denote its length by $T_0$. After the season began, some cities opened their stadiums while others kept them closed. We refer to this horizon as the {\em post-treatment} period and denote its length as $T_1$. Let $\Ic_0$ represent the set of cities whose teams remained closed and $\Ic_1$ the set of cities whose teams opened their stadiums. Thus, for $t \le T_0$, we observe $Y_{it} = Y_{it}(0)$ for every city $i$, whereas for $t > T_0$, we observe $Y_{it} = Y_{it}(1)$ for $i \in \Ic_1$ and $Y_{it} = Y_{it}(0)$ for $i \in \Ic_0$. Notably, our formulation rules out interference and dynamic effects. With this notation, the first question asks for the missing counterfactual trajectory $Y_{it}(0)$ for $i \in \Ic_1$ and $T_0 < t \le T \coloneqq T_0 + T_1$. That is, for a treated city that opened its stadium, what would its case rates have been during the observed post-treatment period had it remained closed?
The synthetic controls (SC) methodology provides an elegant answer to this question abadie1, abadie2. Taking $i$ to be Atlanta, the synthetic control method constructs a “synthetic Atlanta” as a weighted composition of control cities that kept their stadiums closed. The weights are chosen so that the synthetic city closely tracks Atlanta during the pre-treatment period:
where $\mathcal{W}$ is a constraint set. The counterfactual path for Atlanta under closure is then estimated by
Comparing the observed Atlanta trajectory, $Y_{it} = Y_{it}(1)$, to its synthetic closed-stadium trajectory, $\widehat{Y}_{it}(0)$, provides an estimate of the effect of opening the stadium. Applying this type of analysis, nfl_pnas found little evidence that opening NFL stadiums substantially increased local COVID-19 case counts, suggesting that, for the 2020 NFL season, the benefits of providing a tightly controlled outdoor spectating environment, counterbalanced the risks associated with opening.
While this first analysis is valuable, it leaves a complementary question unanswered: {\em what would have happened to a city that kept its stadium closed had it instead opened?} For example, what would COVID-19 case rates in Seattle have looked like had the Seahawks admitted fans? This second question reverses the direction of the missing counterfactual. Instead of estimating $Y_{it}(0)$ for treated cities $i \in \Ic_1$, we now seek $Y_{it}(1)$ for $i \in \Ic_0$ during $T_0 < t \le T$. The synthetic interventions (SI) methodology addresses this question by extending the SC-logic to multiple intervention states synth_iv. Taking $i$ as Seattle, the SI estimator constructs a “synthetic Seattle” from cities that opened their stadiums, whereby the weights are again chosen to match Seattle during the pre-treatment period:
The counterfactual open-stadium trajectory is then estimated as
The gap between $\widehat{Y}_{it}(1)$ and the observed closed-stadium trajectory, $Y_{it} = Y_{it}(0)$, estimates what would have happened had a control city opened its stadium during the observed post-treatment period.
Together, SC and SI answer two natural questions. SC asks what would have happened to treated cities had they remained untreated. SI asks what would have happened to control cities had they instead received treatment. Both methods, however, look {\em backward}: they reconstruct missing potential outcomes during a period that has already been observed for at least one donor group. Now consider a decision-maker in a control city midway through the 2020 NFL season. Up to that point, the city has kept its stadium closed, while other cities have already opened. The decision-maker may certainly ask the retrospective question: what would our case rates have been had we opened at the start of the season? SI is designed for precisely this kind of analysis. But the decision-maker’s more immediate policy question is prospective: {\em what will happen if we open for the next game?}
This question differs in kind from the previous two. It does not ask for a missing potential outcome during a period that has already been observed for other cities. Instead, it asks for a counterfactual forecast beyond the observed panel: $Y_{it}(1)$ for $i \in \Ic_0$ and $t > T$. For a control city that has remained closed through time $T$, we want to forecast what would happen if it switches to treatment after $T$. The challenge is that no outcomes are observed beyond $T$, for either treated or control cities.
As currently formulated, neither SC nor SI can directly answer this question. Both methods rely on donor-unit outcomes during the target period of interest. SC estimates a treated unit’s missing control trajectory by using observed control-unit outcomes over the same post-treatment period. SI estimates a control unit’s missing treated trajectory by using observed treated-unit outcomes over the same post-treatment period. Once the target period lies beyond the observed data, the donor outcomes needed to construct these counterfactuals are themselves unobserved. A natural alternative is to turn to time-series forecasting. However, traditional time-series methods answer yet another question. They forecast what will happen if the current regime continues: what will happen to Atlanta if it continues operating under its observed open-stadium policy, or what will happen to Seattle if it continues operating under its observed closed-stadium policy? They are not, by themselves, designed to forecast what will happen if a unit switches to a treatment condition it has not yet experienced.
This distinction motivates the central question of this article: {\em Can we forecast beyond the observed panel to predict what will happen to a control unit under a treatment it has not yet experienced?} Equivalently, can we estimate $Y_{it}(1)$ for a control unit $i \in \Ic_0$ and a future time $t > T$, using only observations up to time $T$?
This article proposes a framework for causal forecasting in panel data. The framework combines the retrospective counterfactual logic of SC and SI with the prospective structure of time-series forecasting. Accordingly, our target estimand is not a retrospective treatment effect during an already observed post-treatment period. Instead, for a control unit that has not yet received treatment, we aim to forecast its future potential outcome under treatment beyond the observed panel.
Our approach builds on the latent-factor foundations of SC and SI. Classical SC methods are commonly motivated by matrix factor models in which potential outcomes depend on latent unit and time factors. SI extends this idea to settings with multiple treatment states by modeling potential outcomes through latent unit, time, and intervention factors, with the control condition serving as the baseline intervention. We build upon this structure by imposing time-series dynamics on the latent time factors. This added structure allows us to use the same observation patterns that arise in SC and SI designs, while expanding the class of questions from retrospective counterfactual reconstruction to prospective causal forecasting.
Under this model, we establish an identification result for causal forecasts of the kind described above. To operationalize this strategy, we introduce the Two-Way Synthetic Forecasting (TWSF) estimator, which combines two sources of structure. First, it learns cross-unit relationships from the pre-treatment period, as in SI. Second, it learns temporal dynamics from the post-treatment trajectories of treated units, as in multivariate time-series forecasting. These two components are then combined to forecast the treated potential outcome of a control unit beyond the observed horizon.
We prove that TWSF is pointwise consistent under suitable signal-to-noise, rank, and spectral conditions. We then introduce an orthogonalized correction that removes first-order sensitivity to the estimated unit and time weights, yielding an asymptotically normal estimator and valid pointwise inference for one-step-ahead causal forecasts. We further extend the method to fixed multi-step forecasting horizons through both direct and recursive procedures, showing that analogous pointwise inferential guarantees continue to hold.
Finally, we support our statistical claims through simulation studies and revisit the NFL stadium-opening study of nfl_pnas. In the empirical application, we ask not only what would have happened during the observed season had cities made different opening decisions, but also what would have happened had closed-stadium cities switched to opening at a future decision point using only information available up to that time. The empirical findings are broadly consistent with the original NFL analysis: opening stadiums does not appear to generate a substantial increase in local COVID-19 case rates, although the timing of the opening decision can matter.
This article connects two literatures that have largely developed separately: causal inference with panel data and time-series forecasting. The former is designed to estimate missing potential outcomes, but typically within an observed post-treatment period. The latter is designed to extrapolate beyond the observed horizon, but usually under the continuation of an already observed regime. We combine these perspectives in order to forecast a future counterfactual outcome for a unit that has not yet experienced the treatment.
Causal inference with panel data has been studied extensively. SC, introduced in the seminal works of abadie1, abadie2, is a prominent tool with many methodological variants hcw, sc_enet, Xu2016GeneralizedSC, li_bell, RSC, arco, li2020, abadie_hour, victor21. Related extensions adapt SC to staggered adoption eli_stagger and leverage multiple outcomes mRSC, sun25. Particularly relevant is SI, which generalizes SC to accommodate multiple interventions synth_iv. Broadly speaking, these methods are connected through low-rank factor models and matrix-completion perspectives bai2019matrix, fernandez2020low, CAHAN2023113. Additionally, most theoretical guarantees in this literature are established with global error metrics, such as the Frobenius norm, though recent work has developed entrywise guarantees and pointwise inference Choi20092024, yan2024entrywiseinferencemissingpanel. Our work shares this pointwise inferential perspective, but targets a prospective causal forecast beyond the observed panel.
A prominent work on the relationship between unit-side regressions, such as SC, and time-side regressions from the unconfoundedness literature rubin_rosenbaum, imbens_wooldridge is athey2021matrix. Along this thread, sameroot clarifies how these approaches coincide in estimation while differing in inference. Operationally, augmented SC asc and synthetic difference-in-differences sdid, combine unit- and time-side adjustments into doubly robust estimators. Recent approaches have also incorporated temporal structure into panel models. These include neural-network methods combined with SC synbeats, Bayesian SC with autoregressive factors Pang_Liu_Xu_2022, multitask Gaussian-process models for dynamic treatment effects eli23, chen23, and recurrent or continuous-time neural methods for counterfactual response modeling lim18, seedat22, vand23. These methods offer flexible ways to model time-varying potential outcomes, but do not conduct forecasting and often lack theoretical guarantees.
A close methodological precursor is the line of work connecting matrix estimation with time-series forecasting ssa, mssa. In particular, multivariate singular spectrum analysis (mSSA) forecasts the continuation of a multivariate time series under an observed regime by exploiting low-dimensional structure in trajectory matrices mssa; notably, it does not address counterfactual regime switches. Nevertheless, this perspective informs our treatment of time dynamics. The closest work to ours is focus, which also studies causal forecasting in panel data using a low-rank factor model and provides pointwise estimation and inference guarantees. Their setting is complementary. Observationally, they require target units to have previously experienced treatment, whereas we study target units that remain untreated throughout the observed panel. From a modeling standpoint, they impose a stationary autoregressive structure on the temporal factors, while we impose a low-rank temporal structure.
Section (ref) introduces the causal forecasting panel framework and defines the one-step-ahead causal estimand. Section (ref) presents the TWSF estimator. Section (ref) establishes pointwise consistency and finite-sample forecasting error bounds. Section (ref) introduces the orthogonalized TWSF estimator and proves pointwise asymptotic normality, enabling valid inference. Section (ref) extends the framework to fixed multi-step forecasting horizons using direct and recursive procedures. Section (ref) presents simulation studies. Section (ref) revisits the NFL stadium-opening application. Section (ref) concludes this article. The proofs are collected in the Appendix.
For any positive integer $a$, let $[a] = \{1, \dots, a\}$. For a vector $\boldsymbol{v} \in \Rb^a$, let $\| \boldsymbol{v} \|_p$ denote its $\ell_p$-norm. We define the inner product between vectors $\boldsymbol{u}, \boldsymbol{v} \in \Rb^a$ as $\langle \boldsymbol{u}, \boldsymbol{v} \rangle = \boldsymbol{u}^\top \boldsymbol{v} = \sum_{\ell=1}^a u_a v_a$. For a matrix $\bX \in \Rb^{a \times b}$, we denote its operator and Frobenius norms as $\| \bX \|_\text{op}$ and $\| \bX \|_F$. Let $\| \boldsymbol{x} \|_{\psi_2}$ denote the sub-Gaussian norm of a random vector $\boldsymbol{x}$. Denote by $\boldsymbol{I}$ the identity matrix and by $\bX^\dagger$ the Moore-Penrose pseudoinverse of $\bX$. Convergence in probability and distribution are denoted by $\xrightarrow{p}$ and $\rightsquigarrow$, respectively.
We abstract away the motivating example in Section (ref), while continuing to use the NFL study as a concrete guide. Consider a panel data of observations on $N \ge 1$ units across $T \ge 1$ time periods. To articulate our causal objectives, we adopt the potential outcomes framework attributed to neyman and rubin1976. For each unit $i \in [N]$ and time period $t \in [T]$, let $Y_{it}(0)$ and $Y_{it}(1)$ denote the potential outcomes under control and treatment, respectively. Let $D_{it} \in \{0,1\}$ denote the treatment indicator. The observed outcome is
The observation law of (ref) encodes the {\em stable unit treatment value assumption}, which implicitly rules out spillover (network) effects between units imbens_rubin_2015. Moreover, it highlights the central challenge of causal inference: for any unit-time pair, only one potential outcome can be observed.
We focus on an observation pattern that captures the NFL setting in Section (ref). As in the canonical SC setup, all $N$ units are observed under control during the first $T_0$ time periods, which we call the pre-treatment period. During the following $T_1$ time periods, which we call the post-treatment period, the $N$th unit remains under control, while all other units receive treatment. For a visualization, see Figure (ref). Formally, for all $i \in [N]$ and $t \le T_0$, $D_{it} = 0$. For $T_0 < t \le T$, where $T \coloneqq T_0 + T_1$,
Accordingly, the control set is $\Ic_0 = \{N\}$ and the treated donor pool is $\Ic_1 = [N-1]$. For compactness, let $\bD \in \{0,1\}^{N \times T}$ collect the treatment assignments and $\boldsymbol{Y} \in \Rb^{N \times T}$ collect the observed outcomes. We reemphasize that the observed panel ends at time $T$. Our objective, however, is prospective: we aim to forecast the treated potential outcome of the control unit beyond the observed panel.
To make this forecasting problem identifiable, we impose structure on the potential outcomes. Specifically, we adopt the SI factor model introduced in synth_iv.
The key feature of Assumption (ref) is that the unit factors $\boldsymbol{u}_i$ are invariant across time and intervention states. This invariance is what enables relationships learned among units under one intervention to be transported to outcomes under another intervention. In particular, weights learned from pre-treatment outcomes under control can be used to reconstruct treated potential outcomes, provided the latent unit structure is stable. More flexible models are possible. For example, one could allow the unit factor $\boldsymbol{u}_i$ to depend on the state $d$. However, doing so would substantially increase model complexity and would weaken the cross-intervention link that allows us to extend SC reasoning to the prospective forecasting problem considered in this article. Model (ref) can also be viewed as a multi-intervention extension of the multiplicative tensor factor models that underpin the tensor completion literature tensor_missing, gandy, anandkumar2014tensor, barak2015noisy. From this perspective, Assumption (ref) imposes a natural low-dimensional structure that allows time-intervention factors to vary flexibly while keeping unit factors stable.
Under the factor model in (ref), we define our one-step-ahead causal forecasting estimand as
In words, (ref) is the expected potential outcome for the control unit $N$ under treatment one period beyond the observed panel. Returning to the NFL example, suppose a policy-maker in Tennessee is deciding whether to admit fans at the next home game on October 13, 2020. By then, the policymaker has roughly one month of case-count data from cities whose NFL teams had already opened their stadiums. The relevant policy question is whether Tennessee should continue its closed-stadium policy or switch to an open-stadium policy. In this context, $\theta$ represents Tennessee's expected COVID-19 case count on October 14 if Tennessee admits fans on October 13.
For now, we focus on this one-step-ahead estimand. Section (ref) extends the object to a fixed multi-step forecasting horizon.
A policy-maker is unlikely to make such a decision arbitrarily. Instead, the decision will likely depend on factors related to the county, trends, and anticipated case trajectories under each policy. We therefore allow treatment assignment to be confounded, provided the confounding is mediated by the latent factors. Let $\Ec \coloneqq \{\boldsymbol{u}_i, \boldsymbol{v}_t(d), D_{it}: i \in [N], t \in [T+1], d \in \{0,1\} \}$ denote the information set comprising the latent factors and intervention assignments.
Assumption (ref) says that, conditional on the latent structure, the idiosyncratic noise is mean zero. Together with Assumption (ref), it implies that potential outcomes are mean independent of treatment assignment once the latent factors are held fixed. Thus, the latent factors play the role of unobserved confounders, analogous to observed covariates under a classical selection-on-observables assumption.
This type of latent ignorability is aligned with assumptions used in the causal panel data literature athey2021matrix, asc, synth_iv and in related contexts kallus2018causal. As with any assumption about unobserved confounding, Assumption (ref) is not directly testable and should be assessed using domain knowledge about the data-generating process.
Assumption (ref) requires the latent factor of the control unit to lie within the linear span of the latent factors of the treated donor units. This condition is what allows unit $N$ to be reconstructed as a weighted combination of units in $\Ic_1$. From a linear-algebraic perspective, Assumption (ref) is closely related to the factor structure in Assumption (ref). When the factor model is low rank, with $r \ll \min\{N, 2T\}$, the latent unit factors are necessarily linearly dependent, rendering Assumption (ref) a natural consequence. Even so, low rank alone does not guarantee that the donor units span $\boldsymbol{u}_N$ as pathological configurations can still arise in which the donor pool does not cover the target unit. We therefore state Assumption (ref) separately, because it directly governs the identifiability of unit-level causal forecasts.
In practical terms, Assumption (ref) requires the donor pool $\Ic_1$ to be sufficiently large and rich. The condition is thus analogous to common support in standard causal inference, where treated and control units must overlap in their covariate distributions. Although the latent factors are unobserved, the plausibility of Assumption (ref) can be assessed indirectly through pre-treatment fit diagnostics, as is standard in the SC literature abadie_survey.
The following proposition, akin to synth_iv, formalizes the spatial reconstruction implied by our assumptions thus far.
Proposition (ref) shows that the target unit $N$ can be reconstructed from the treated donor units $\Ic_1$. Namely, Proposition (ref)(b) provides a way to learn the relevant unit weights during the pre-treatment period, when all units are observed under control. However, Proposition (ref)(a) still involves donor outcomes at time $T+1$, which lies beyond the observed panel. Since no outcomes are observed after time $T$, Proposition (ref)(a) alone is not enough to recover $\theta$. We therefore impose additional structure on the treated latent time factors.
To extrapolate beyond the observed panel, we endow the latent time factors under treatment with additional structure. Assumption (ref) implies that there are $r$ canonical latent temporal components under treatment. For each component $a \in [r]$, define the time series $g_{a}(t) \coloneqq V_{ta} (1)$ for $t \in \mathbb{Z}$. For any positive integers $m, n$ and shift $s \in \Zb$, we define the Hankel matrix generated by $g_a$ as $\bH(g_a; m, n, s) \in \Rb^{m \times n}$ with entries $H_{ij}(g_a; m, n, s) \coloneqq g_a(s + i + j - 2)$ for $i \in [m]$ and $j \in [n]$. Equivalently,
Motivated by mssa, our temporal assumption is that each latent treated time factor admits a low-rank Hankel representation.
Assumption (ref) is flexible enough to capture a rich family of time series dynamics, including any finite sum of harmonic-polynomial components, differentiable periodic functions, and H\"older continuous functions mssa. Hence, the low-rank Hankel structure provides natural finite-dimensional representations or approximations for many time series with trend and periodicity. Spectral representations of stationary processes further motivate this structure: sample paths of many stationary processes can be expressed through weighted combinations of harmonics, with richer representations obtained as the number of harmonics increases shumway2011time. The class of time series with low-rank Hankel structure is also closed under component-wise addition and multiplication mssa.
We next translate this temporal structure into a form that can be used for estimation. Fix a Page length $L \in \mathbb{Z}_+$ satisfying $L \le T_1$, and suppose $B \coloneqq T_1/L$ is an integer. For each donor $j \in \Ic_1$, define the treated signal $f_j(t) = \langle \boldsymbol{u}_j, \boldsymbol{v}_t(1) \rangle$ for $t \in \Zb$. The corresponding Page matrix over the post-treatment period is
Stack the donor Page matrices as $\bar{\bP}_L = [\bP_L(f_1) \cdots \bP_L(f_{N-1})] \in \Rb^{L \times B N_1}$. Notably, each Page matrix can be obtained by selecting the appropriate columns from the corresponding Hankel matrix. Therefore, the rank of the Hankel matrix controls the rank of the Page matrix. The next result, attributed to mssa, makes this precise.
Proposition (ref) asserts that the stacked Page matrix for the treated donor signals is low rank, with rank bounded by the product of the factor-model rank $r$ and the Hankel rank $G$. This low-rank temporal structure yields a forecasting relation analogous to the spatial reconstruction in Proposition (ref).
Proposition (ref) is similar in spirit to mssa. Proposition (ref)(a) expresses each donor unit's treated outcome at time $T+1$ as a common linear function of its recent treated outcomes. Proposition (ref)(b) demonstrates how this same temporal relationship can be learned from the post-treatment period, when all donor units in $\Ic_1$ are observed under treatment. Operationally, it suggests forming the stacked Page matrix over the donor units and regressing the final row of each Page block on the preceding $L-1$ rows.
Combining the spatial reconstruction in Proposition (ref) with the temporal reconstruction in Proposition (ref) yields the following identification result.
Theorem (ref) shows that the causal forecasting estimand $\theta$ can be represented using only donor-pool outcomes observed up to time $T$. The unknown quantities enabling this representation are the spatial weights $\boldsymbol{\beta}$ and temporal weights $\boldsymbol{\alpha}$. Proposition (ref)(b) indicates that $\boldsymbol{\beta}$ can be learned from pre-treatment outcomes, while Proposition (ref)(b) indicates that $\boldsymbol{\alpha}$ can be learned from post-treatment donor trajectories.
Importantly, Theorem (ref) allows treatment assignment to depend on latent factors. Since the latent factors determine the conditional mean of the potential outcomes, treatment assignment may depend on the signal itself. In the NFL setting, this means that cities may choose whether to open or close stadiums based on how they expect case counts to evolve under each policy. The framework therefore accommodates a broad class of unobserved confounding mechanisms while still identifying a prospective causal forecast.
Theorem (ref) suggests a constructive estimation strategy. To estimate the causal forecast, we need to learn two sets of weights. The first set captures cross-unit structure: it reconstructs the target control unit from the treated donor units using pre-treatment outcomes. The second set captures temporal structure: it forecasts future treated outcomes from recent treated histories using the post-treatment trajectories of the donor units. We call the resulting procedure the Two-Way Synthetic Forecasting estimator, or TWSF.
Let $L \in \mathbb{Z}_+$ be an algorithmic hyperparameter. Define $K \coloneqq L - 1$, $B \coloneqq T_1 / L \in \mathbb{N}$, and $M \coloneqq (B-1)N_1$. To learn the spatial structure, collect the pre-treatment outcomes of the target unit $N$ as
and collect the pre-treatment outcomes of the treated donor units $j \in \Ic_1$ as
To learn the temporal structure, encode the post-treatment observations for each $j \in \Ic_1$ into a training Page matrix. Specifically, define
The final row of this matrix contains the one-step-ahead responses within each Page block, while the first $K$ rows contains the corresponding lagged predictors. Stacking these training Page matrices across $\Ic_1$ gives
Here, $\bZ_\texttt{lag} \in \Rb^{K \times M}$ contains the first $K$ rows of the stacked Page matrix, while $\boldsymbol{z}_{\texttt{next}} \in \Rb^M$ contains the final row. The training Page matrix deliberately excludes the final $L$ post-treatment observations as they are reserved for forecasting the target outcome at time $T+1$. Accordingly, collect the final $K$ observed lags for the treated donor units in the matrix
A key ingredient of TWSF is a spectral de-noising step based on hard singular value thresholding (HSVT). For any matrix $\bA \in \Rb^{m \times n}$, write its singular value decomposition (SVD) as
For any $k \le \min\{m, n\}$, define the rank-$k$ approximation of $\bA$ as
The TWSF algorithm proceeds in four steps.
To interpret TWSF, return to the example in which the target unit $N$ is Tennessee. The vector $\widehat{\boldsymbol{\beta}}$ describes the weighted combination of treated donor cities whose pre-treatment case-count trajectories best reconstruct Tennessee’s pre-treatment trajectory. In contrast, the vector $\widehat{\boldsymbol{\alpha}}$ describes the weighted combination of recent treated lags that best predicts the next treated outcome among cities whose NFL teams had already opened their stadiums. The procedure to acquire $\widehat{\boldsymbol{\beta}}$ is motivated by SI, while the the strategy to learn $\widehat{\boldsymbol{\alpha}}$ follows mSSA.
TWSF admits two interpretations. In the first interpretation, TWSF begins with an SI-step via $\widehat{\boldsymbol{\beta}}$ to impute Tennessee's counterfactual trajectory had it opened its stadium during the final $K$ observed post-treatment periods, characterized by $\bW^\top \widehat{\boldsymbol{\beta}} \in \Rb^K$. The estimator then treats this imputed trajectory as the recent treated history of Tennessee and applies the forecasting weights $\widehat{\boldsymbol{\alpha}}$ to project one period ahead, yielding $\widehat{\theta}$. In the second interpretation, TWSF reverses the order. It starts with an mSSA-step via $\widehat{\boldsymbol{\alpha}}$ to forecast the time $T+1$ treated outcome for each treated donor unit, collected in $\bW \widehat{\boldsymbol{\alpha}} \in \Rb^{N_1}$. It then applies the SI-weights $\widehat{\boldsymbol{\beta}}$ to these donor forecasts, producing the same estimate $\widehat{\theta}$. Thus, \texttt{TWSF} can be read either as “impute, then forecast” or as ”forecast, then impute.” In both cases, the estimate arrives at the same bilinear form.
This bilinear structure resembles the form of doubly robust estimators such as synthetic difference-in-differences sdid and augumented SC asc, which also combine unit- and time-side adjustments. The interpretation, however, is different. In those methods, time weights are used to balance pre-treatment and post-treatment periods in order to sharpen a treatment-effect contrast. In TWSF, time weights are forecasting coefficients learned from the low-rank Hankel structure of the treated latent time factors.
The de-noising steps in (ref) and (ref), followed by the regression steps in (ref) and (ref), are instances of principal component regression (PCR). PCR first projects the observed design matrix onto its leading singular components and then performs least squares using the resulting low-rank approximation. This formulation is natural under Assumption (ref). The observed matrices $\boldsymbol{Y}_{\Ic_1, \texttt{pre}}$ and $\bZ_{\texttt{lag}}$ are contaminated versions of their conditional means, $\mathbb{E}[\boldsymbol{Y}_{\Ic_1, \texttt{pre}} \mid \Ec]$ and $\mathbb{E}[\bZ_\texttt{lag} \mid \Ec]$, which have low-rank structure. The leading singular values of the observed matrices therefore capture the underlying signal, while smaller singular directions primarily reflect noise. PCR uses this spectral separation to discard noisy directions before estimating the spatial and temporal weights.
In practice, the true ranks of the signal matrices are unknown. Several principled strategies can be used to choose the spectral thresholds $k_y$ and $k_z$. Standard cross-validation procedures provide one option. Other data-driven alternatives include selecting the number of principal components needed to explain a target fraction of spectral energy, or using universal singular-value thresholding rules that retain only singular values above a data-dependent cutoff, such as the procedures proposed in donoho14 and Chatterjee15.
We next establish statistical guarantees for the TWSF estimator. Throughout this section, recall the definitions $B \coloneqq T_1 / L \in \mathbb{N}$, $K \coloneqq L - 1$, and $M \coloneqq (B-1)N_1$.
Before turning to estimation error, we address an important subtlety. The TWSF estimator uses PCR to estimate the coefficient vectors $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$. However, if the population design matrices are rank deficient, these coefficient vectors need not be uniquely identifiable. For example, if $\mathbb{E}[\boldsymbol{Y}_{\Ic_1, \texttt{pre}} \mid \Ec]$ is rank deficient, then the pre-treatment observations cannot identify $\boldsymbol{\beta}$ itself. They identify only the component of $\boldsymbol{\beta}$ lying in the row space of the population design matrix. An analogous issue arises for $\boldsymbol{\alpha}$. The following assumption ensures that these recoverable components are sufficient for estimating the causal forecast.
Assumption (ref) requires the out-of-sample forecasting block $\mathbb{E}[\bW \mid \Ec]$ to lie within the linear span of the in-sample covariates used to learn the spatial and temporal weights. The first inclusion says that the unit-side variation in $\mathbb{E}[\bW \mid \Ec]$ is represented in the pre-treatment donor design. The second says that the time-side variation in $\mathbb{E}[\bW \mid \Ec]$ is represented in the lagged treated-donor design. This condition allows the estimable projections of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ to identify the same causal forecast as the original weights, which we formalize next.
Corollary (ref) reformulates Theorem (ref) in terms of the recoverable components $\boldsymbol{\alpha}^*$ and $\boldsymbol{\beta}^*$. This is essential for estimation because these projected quantities are uniquely defined and identifiable from the population regression problems. In this sense, Assumptions (ref) and (ref) are complementary. Assumption (ref) guarantees that a valid spatial representation exists, while Assumption (ref) ensures that the component of this representation recoverable from the pre-treatment design is sufficient for forecasting. Similarly, Assumptions (ref) and (ref) together ensure that the relevant temporal forecasting relationship is both available and attainable.
More broadly, Assumption (ref) is a form of generalization condition. It requires the out-of-sample covariates in $\mathbb{E}[\bW \mid \Ec]$ to remain within the linear span of the in-sample covariates used for learning. This plays a role analogous to covariate-distribution similarity in standard statistical learning, but it is expressed in a linear-algebraic form suited to panel settings with treatment-induced distribution shifts.
We now impose additional conditions that allow us to bound the forecasting error of TWSF.
The precise bound $[-1,1]$ is without loss of generality and can be extended to $[a,b]$ for any $a, b \in \Rb$.
Assumption (ref) allows the latent temporal factors to be correlated and to follow rich time-series dynamics, but it requires the idiosyncratic shocks to be independent. This can be restrictive, but it is useful for a first finite-sample analysis of an estimator of the TWSF form. Similar independence structures are used in early analyses of SC, SI, and mSSA analysis abadie2, synth_iv, mssa, and in related work on causal forecasting such as focus, albeit they consider a different temporal factor structure. Extending the theory to more general dependent noise structures is an important direction for future work.
Assumption (ref) requires the nonzero singular values of the population design matrices to be sufficiently strong and well balanced. This condition ensures that the low-rank signal can be reliably separated from noise. Although potentially restrictive, it formalizes the idea that the signal should not be overwhelmed by idiosyncratic variation. It is therefore analogous to spectral gap, pervasiveness, and beta-min conditions used in factor models, matrix completion, and high-dimensional regression chamberlainfactor, bai2019matrix, fan2018eigenvector, beta_min. However, unlike traditional factor-based methods, our aim is not to recover the latent factors themselves but to exploit their presence and separation from noise to guide estimation. Empirically, the plausibility of this condition can be assessed through spectral diagnostics such as scree plots and estimated condition numbers.
We now state finite-sample guarantees for the TWSF estimator. Throughout, the notation $\lesssim$ suppresses constants that do not depend on the model dimensions.
The first result controls the estimation errors of the PCR subroutines used to recover $\boldsymbol{\beta}^*$ and $\boldsymbol{\alpha}^*$.
Proposition (ref) shows that the spatial and temporal PCR subroutines consistently recover the identifiable components of the population weights under suitable conditions, discussed below. These rates improve upon the best known PCR bounds, as established in pcr_jmlr. To see this, take $N_1 \asymp T_0 \asymp d$ and $\sigma, r_y = \Oc(1)$. Then, ignoring log factors, pcr_jmlr yields an error rate of $\Oc(d^{-3/4})$, whereas Proposition (ref) yields an error rate of $\Oc(d^{-1})$. This improvement is potentially useful beyond the present setting, since the related SI and mSSA estimators studied in synth_iv and mssa also rely on PCR subroutines and establish rates that are aligned with pcr_jmlr.
It is important to highlight that Proposition (ref) hinges on several rank conditions. To begin, note that (ref) can be interpreted as a signal-to-noise requirement. Standard concentration results imply that the operator-norm noise in $\boldsymbol{Y}_{\Ic_1, \texttt{pre}}$ scales as $\sqrt{N_1} + \sqrt{T_0}$, up to logarithmic factors. Assumption (ref) ensures that the smallest nonzero singular value of the population signal $\mathbb{E}[\boldsymbol{Y}_{\Ic_1, \texttt{pre}} \mid \Ec]$ is sufficiently large relative to this noise level. The same reasoning applies to $\bZ_\texttt{lag}$. Moreover, the proposition assumes oracle knowledge of the true ranks $r_y$ and $r_z$. A formal treatment of rank misspecification is left for future work, but in practice data-driven methods such as the thresholding scheme of donoho14 performs remarkably well, and overestimating the rank has been shown to be less harmful than underestimating it pcr_jmlr.
We now translate the parameter estimation bounds into a finite-sample error bound for the causal forecast.
Theorem (ref) provides a pointwise causal forecasting guarantee under unobserved confounding. To interpret the rate, consider a balanced regime in which $\sigma = \Oc(1)$ and $N_1 \asymp T_0 \asymp K \asymp M \asymp d$. If $r_y \asymp r_z \asymp r_d$, then
Hence, if $r_y, r_z = \Oc(1)$, then the bound simplifies as $R_{\text{HP}} = \Oc(d^{-1/2} \log(d)) = o(1)$, so the TWSF estimator is pointwise consistent for $\theta$.
\setcounter{assumption}{0}
Theorem (ref) shows that the TWSF estimator can achieve vanishing causal forecast error under suitable spectral and rank conditions. However, this high-probability error bound does not by itself provide a tractable limiting distribution to enable inference. The challenge is that the estimator depends on two estimated nuisance objects, $\boldsymbol{\beta}^*$ and $\boldsymbol{\alpha}^*$. Namely, a first-order expansion of the estimator contains terms such as
These terms are linear in the PCR estimation errors, and are not automatically negligible relative to the stochastic fluctuation of the estimator, and their distribution is arduous to characterize directly. Therefore, although the vanilla TWSF estimator is natural for prediction, we introduce a bias-correction based on orthogonalization to enable inference.
Motivated by the semiparametric and debiased machine-learning literature on locally robust scores Newey1994, ChernozhukovEtAl2018DML, ChernozhukovEtAl2022LocallyRobust, ChernozhukovNeweySingh2022Automatic, ChernozhukovNeweySingh2022Riesz, HirshbergWager2021, we construct an orthogonalized TWSF estimator using Riesz representers associated with the two residual moment equations used to estimate the nuisance parameters. For compactness, write $\bar{\bX} = \mathbb{E}[\bX \mid \Ec]$ for any random object $\bX$. Consider an augmented population score
The two residual moments vanish under Corollary (ref) with Propositions (ref) and (ref), so (ref) holds for any choice of $(\boldsymbol{q}_\alpha, \boldsymbol{q}_\beta)$. We choose these representers so that the score is first-order insensitive to perturbations in $(\boldsymbol{\alpha}^*, \boldsymbol{\beta}^*)$. Differentiating with respect to the nuisance parameters yield the orthogonality requirements
Assumption (ref) ensures that these equations are feasible. Indeed, $\bar{\bW} \boldsymbol{\alpha}^* \in \text{col}(\bar{\bW}) \subseteq \text{col}(\bar{\boldsymbol{Y}}_{\Ic_1, \texttt{pre}})$, and thus,
Analogously, $\bar{\bW}^\top \boldsymbol{\beta}^* \in \text{row}(\bar{\bW}) \subseteq \text{row}(\bar{\bZ}_\texttt{lag}^\top)$, and hence,
We therefore define the population Riesz representers as the minimum $\ell_2$-norm solutions
These representers do not change the population target since the residual moments vanish. Their role is instead to cancel the first-order sensitivity of the score to nuisance-estimation error induced by $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\beta}}$.
The orthogonalized TWSF estimator replaces the population matrices in the Riesz equations by rank-truncated estimates and evaluates the corresponding empirical orthogonal score.
The first term in (ref) is the vanilla TWSF estimator. The second and third terms are residual corrections from the unit-side and time-side regressions. These corrections remove the leading effect of estimating $\boldsymbol{\beta}^*$ and $\boldsymbol{\alpha}^*$, enabling asymptotic normality.
For inference, we strengthen the previous noise and spectral conditions of Assumptions (ref) and (ref).
Assumption (ref) specializes Assumption (ref) to homoskedastic Gaussian noise. It is not essential in principle, but it avoids additional Lyapunov-type conditions and allows us to state a cleaner asymptotic normality result.
Assumption (ref) extends the well-balanced spectrum condition to the forecasting block $\bW$. This ensures that the Riesz representers can be estimated accurately enough for valid inference.
For the remainder of this article, define
Theorem (ref) shows that the orthogonalized TWSF estimator is asymptotically normal around the causal forecast $\theta$, provided the remainder term $\Psi$ is negligible relative to the leading standard deviation $\upsilon$. To interpret the condition, consider the following balanced regime.
Under this regime, Theorem (ref) states
If, in addition, we impose the lower bound $\| \boldsymbol{q}^*_\beta \|_2 + \| \boldsymbol{q}^*_\alpha \|_2 \gtrsim d^{-1/2}$, then $\upsilon \gtrsim d^{-1/2}$ and thus,
Hence, when $r_d \gg \log(d)$, the condition $\Psi/ \upsilon = o(1)$ is satisfied if $r_d = o(d^{1/3})$. When $r \lesssim \log(d)$, it is sufficient that $r_d \log^2(d) = o(d)$. In the particularly simple case where the ranks remain bounded, the condition follows immediately.
To implement Theorem (ref) for conducting inference, we need to estimate the leading variance $\upsilon^2$.
Consider the plug-in estimator
where $\widehat{\sigma}^2$ is any estimate of $\sigma^2$.
Proposition (ref) shows that the plug-in variance estimator is consistent and the studentized error is asymptotically normal under the same general scaling regimes that make the orthogonalized remainder negligible. Namely, in the balanced regime of Definition (ref), together with $\| \boldsymbol{q}^*_\beta \|_2 + \| \boldsymbol{q}^*_\alpha \|_2 \gtrsim d^{-1/2}$, we have $\upsilon^2/\sigma^2 \gtrsim d^{-1}$. The bound in Proposition (ref) then implies
Thus, when $r_d \gg \log(d)$, the plug-in variance estimator is consistent under the same rank scaling $r_d = o(d^{1/3})$ that appears in Theorem (ref).
Collectively, Theorem (ref) and Propositions (ref) justify pointwise confidence intervals for $\theta$. For any significance level $a \in (0,1)$, an asymptotic $(1-a) \times 100\%$ confidence interval is
where $z_{1-a/2}$ is the upper $a/2$ quantile of the standard normal distribution. This result contributes to the growing literature on inference in panel data and SC-settings by providing pointwise inferential guarantees for prospective causal forecasts victor18, victor21, Choi20092024, yan2024entrywiseinferencemissingpanel.
Proposition (ref) requires a consistent estimator of $\sigma^2$. We consider two variance estimators based on the unit- and time-sided PCR residuals:
In turn, these can also be combined into a pooled estimator:
The following proposition controls the estimation errors of these variance estimators.
Proposition (ref) shows that the proposed residual-based variance estimators are consistent under mild scaling conditions. In the balanced regime of Definition (ref),
Thus, consistency of the variance estimators follows whenever $r_d = o(d)$. Notably, the estimators in (ref) and (ref) are similar in spirit to the residual variance estimators used in synth_iv. However, Proposition (ref) provides explicit high-probability guarantees for their accuracy, which are absent in synth_iv.
So far, our focus has been one-step-ahead causal forecasting: for a control unit that has remained untreated through the observed panel $T$, we looked to forecast its treated potential outcome at time $T+1$. We now extend the framework to a fixed multi-step horizon. For a fixed integer $h \ge 1$, define the horizon-specific causal estimands
For the analysis, let $\Ec_h \coloneqq \Ec \cup \left\{ \boldsymbol{v}_{T+\ell}(1): \ell \in [h] \right\}$ denote the enlarged information set containing the latent variables through the forecast horizon $T+h$. Throughout this section, we treat $h$ as fixed while the panel dimensions grow.
There are two natural ways to extend one-step TWSF to a fixed multi-step horizon. The first is a {\em direct} approach: for each horizon $\ell$, we learn a separate $\ell$-step-ahead forecasting rule from the treated donor trajectories. The second is a {\em recursive} approach: we learn the one-step-ahead forecasting rule once, then iterate it forward to obtain forecasts at horizons $\ell \in [h]$. These two strategies mirror a classical distinction in time-series forecasting, but here they are embedded inside the causal forecasting structure of TWSF.
\setcounter{assumption}{0} \setcounter{cor}{0} \setcounter{theorem}{0}
The direct approach treats each forecast horizon as its own supervised learning problem. Instead of estimating a one-step temporal rule and applying it repeatedly, we construct Page blocks whose first $K$ rows serve as predictors and whose next $h$ rows serve as horizon-specific responses. This allows the estimator to learn a separate temporal coefficient vector for each lead time $\ell \in [h]$. The appeal of this strategy is that it avoids compounding one-step forecast errors. Its cost is that each training block must be longer, so fewer non-overlapping Page blocks are available for estimation.
Fix a lag length $K$, and define $L_h \coloneqq K + h$, $B_h \coloneqq T_1 / L_h \in \Nb$, and $M_h \coloneqq (B_h-1)N_1$. For each treated donor unit $j \in \Ic_1$, define the direct training Page matrix by
Stacking these matrices across treated donor units gives
Here, $\bZ_{\texttt{dir}\text{-}\texttt{lag}} \in \Rb^{K \times M_h}$ contains the first $K$ rows of the stacked Page matrix, while $\boldsymbol{z}_{\texttt{dir}\text{-}\texttt{next}, \ell} \in \Rb^{M_h}$ for $\ell \in [h]$ contains the row for the $\ell$-step-ahead response. The final observed donor histories used for forecasting are collected in $\bW$, as defined in (ref). The direct orthogonalized TWSF estimator proceeds as follows.
The first term in (ref) is the direct multi-step analogue of the vanilla TWSF forecast. The second and third terms are the unit- and time-side orthogonalization corrections. The interpretation is straightforward. The vector $\bW^\top \widehat{\boldsymbol{\beta}} \in \Rb^K$ remains the imputed treated history of unit $N$ over the final $K$ observed periods. The coefficient vector $\widehat{\boldsymbol{\alpha}}^\texttt{dir}_\ell$ maps this treated history directly to the $\ell$-step-ahead outcome. Thus, the direct method estimates $h$ separate forecasting rules, all of which share the same synthetic unit weights $\widehat{\boldsymbol{\beta}}$.
The key novelty relative to the one-step estimator is the construction of a longer Page block of length $K+h$. In the one-step case, the final row of each Page block is the next-period response. In the direct multi-step case, the final $h$ rows are treated as separate future responses. When $h=1$, the direct estimator reduces to the original one-step orthogonalized TWSF estimator of (ref).
The following theorem shows that the direct estimator inherits the pointwise inferential guarantees of the one-step estimator.
Theorem (ref) is essentially a horizon-by-horizon extension of Theorem (ref) and Proposition (ref). Once the direct Page matrix has been constructed, each horizon $\ell$ has the same structure as the original one-step problem. The only substantive change is that the effective number of temporal training samples becomes $M_h$, reflecting the longer Page blocks required by the direct method.
Under the balanced regime of Definition (ref), the same scaling conditions that justify Theorem (ref) and Proposition (ref) also justify Theorem (ref), up to constants involving $h$. Thus, for each fixed horizon $\ell$, one obtains a pointwise confidence interval analogous to (ref):
These intervals are pointwise in $\ell$. A simultaneous confidence band for the full $h$-period forecast path would require a joint Gaussian approximation or a multiple-testing correction, which we leave for future work.
\setcounter{assumption}{0} \setcounter{cor}{0} \setcounter{theorem}{0}
The recursive approach takes the opposite view. Rather than learning a separate forecasting rule for each horizon, it learns the one-step-ahead temporal rule once and iterates that rule forward. This approach is more parsimonious and uses the same Page construction as the original one-step estimator. It is therefore attractive when the post-treatment period is not long enough to support horizon-specific direct regressions. Its cost is that errors in the one-step rule may propagate through the recursion, so the analysis must account for how uncertainty in the one-step coefficient affects all future horizons.
For any vector $\boldsymbol{x} \in \Rb^K$, define the companion map
This matrix updates a $K$-lag state by shifting the lags forward and appending the one-step linear forecast. The recursive orthogonalized TWSF estimator proceeds as follows.
The recursive estimator can be interpreted as follows. The vector $\bW^\top \widehat{\boldsymbol{\beta}}$ is the imputed treated $K$-lag state of the target unit at time $T$. The companion matrix $\boldsymbol{\Pi}(\widehat{\boldsymbol{\alpha}})$ advances this state by one period. Applying the companion map $\ell$ times yields the coefficient vector $\widehat{\boldsymbol{\alpha}}^\texttt{rec}_\ell$, which maps the current imputed treated history directly to the $\ell$-step-ahead forecast.
The main novelty relative to the one-step estimator is that the temporal coefficient is now a nonlinear function of the one-step coefficient $\widehat{\boldsymbol{\alpha}}$. Consequently, the time-side orthogonalization must account for this nonlinearity. This is why the recursive Riesz representer uses the Jacobian $\boldsymbol{J}_\ell(\widehat{\boldsymbol{\alpha}})$. The Jacobian captures how a first-order perturbation in the one-step temporal coefficient propagates to the $\ell$-step recursive coefficient.
The recursive method requires a strengthened recoverability condition.
When $h=1$, Assumption (ref) reduces to the recoverability condition of Assumption (ref) used in the one-step analysis. For larger fixed $h$, it requires not only the final observed lag block $\mathbb{E}[\bW \mid \Ec]$, but also its recursively shifted versions, to remain in the row span of the one-step temporal training design. This is the linear-algebraic price of using one estimated one-step rule to generate multiple future forecasts. Accordingly, it is sensible to expect that Assumption (ref) is more plausible for reasonably lengthed horizons $h$.
Theorem (ref) shows that the recursive estimator also admits valid pointwise inference for each fixed horizon. Relative to Theorem (ref) and Proposition (ref), the new feature is the multiplicative factor $\mathfrak{C}_{h,\alpha}$, which captures the propagation of one-step estimation error through the recursive companion map. For fixed $h$, this factor remains controlled under the same balanced scaling regimes used earlier. Thus, in the balanced regime of Definition (ref), the same rank and signal-to-noise conditions that justify the one-step orthogonalized estimator continue to justify recursive multi-step inference. Theorem (ref) also clarifies the limitation of the recursive approach. The analysis is for fixed $h$. If the forecast horizon were allowed to grow with the panel dimensions, then the stability of the companion recursion would become central, and $\mathfrak{C}_{h,\alpha}$ could no longer be treated as a harmless constant.
For each $\ell \in [h]$, Theorem (ref) supports the pointwise confidence interval
As with the direct approach, these intervals are pointwise in the forecast horizon. A simultaneous confidence band for the entire forecast path would require controlling the joint distribution of the leading Gaussian terms across $\ell \in [h]$, together with the corresponding remainders.
We conduct a Monte Carlo study to evaluate the finite-sample behavior of the direct and recursive orthogonalized TWSF estimators developed in Section (ref). The simulations focus on two aspects of performance: estimation accuracy and pointwise confidence interval coverage. The design is chosen to mirror the balanced fixed-rank regime used to interpret the theoretical results, while allowing the dimension of the panel to grow.
We first describe the simulation setup common to the accuracy and coverage evaluations.
We consider forecast horizons $h \in \{1, 5, 10\}$. In the case $h=1$, the direct and recursive estimators coincide and reduce to the original one-step orthogonalized TWSF estimator introduced in Section (ref). For $h > 1$, the direct estimator learns horizon-specific forecasting rules, while the recursive estimator learns a one-step forecasting rule and iterates it forward.
The simulations are indexed by $d \in \{25, 50, 75, 100, 150, 250\}$. For each $d$, we set $N_1 = T_0 = K =d$, and choose $T_1 \coloneqq 4 (K + 10)$. For the direct estimator at horizon $h$, we use Page blocks of length $K+h$. For the recursive estimator, we use the one-step Page blocks of length $K+1$. In both cases, incomplete trailing Page blocks are discarded.
For each value of $d$, we generate ten independent latent-factor realizations. Conditional on each latent realization, we generate one hundred independent noise realizations, giving $R = 10 \times 100 = 1000$ replications per design point $d$. This nested design reduces sensitivity to any particular latent-factor draw while still averaging over a large number of idiosyncratic noise realizations.
Potential outcomes are generated according to the latent factor model of Assumption (ref) and the observation pattern reflects the setup in Section (ref): prior to $T_0$, all units are observed under control; over the next $T_1$ periods, the donor units are exposed to treatment while the target unit $N$ remains under control. Our target estimand is then defined in accordance with (ref).
For each donor $j \in \Ic_1$, we generate $\boldsymbol{u}_j \coloneqq [1, \xi_{j2}, \xi_{j3}, \xi_{j4}]^\top$ for $j \in \Ic_1$, where $\xi_{j2}, \xi_{j3}, \xi_{j4}$ are standard normal variables. To enforce Assumption (ref), we draw a random subset $\Sc \subset \Ic_1$ with $| \Sc | = 8$, draw convex weights $(\lambda_j)_{j \in \Sc}$ from a symmetric Dirichlet distribution, and set $\boldsymbol{u}_N = \sum_{j \in \Sc} \lambda_j \boldsymbol{u}_j$.
The latent time factors are finite mixtures of harmonics. Let $\boldsymbol{b}(t)$ collect one very low-frequency sine/cosine pair and several medium- and long-period sine/cosine pairs. We then set $\boldsymbol{v}_t(0) = \bA_0 \boldsymbol{b}(t)$ and $\boldsymbol{v}_t(1) = \bA_1 \boldsymbol{b}(t)$, where $\bA_0, \bA_1 \in \Rb^{4 \times 8}$ are fixed loading matrices. These matrices encode how the primitive harmonic components enter the four latent time factors under control and treatment. Each row corresponds to a latent time factor, while each column corresponds to a harmonic basis component; the entries determine the amplitude, sign, and cross-factor mixing of the corresponding harmonic. The first two columns correspond to a very low-frequency harmonic. We set these entries to zero in $\bA_0$, so the control dynamics contain no systematic low-frequency drift, and allow them to be nonzero in $\bA_1$, so the treated dynamics exhibit a smooth upward drift over the observed window. The remaining columns generate medium- and long-period oscillations under both intervention states, with different amplitudes and mixtures. This creates a nontrivial treatment-control contrast while remaining compatible with the low-rank Hankel structure in Assumption (ref).
After forming the raw signals, we scale them so that $\max_{i,t,d} | \langle \boldsymbol{u}_i, \boldsymbol{v}_t(d) \rangle | \le 0.8$, upholding Assumption (ref). We then perturb the signal with independent mean-zero Gaussian noise with $\sigma = 0.15$, following Assumptions (ref) and (ref). The noise level is non-negligible relative to the entrywise signal, but does not overwhelm the leading signal spectra. The loading matrices and scaling are calibrated so that the population design and forecasting blocks have the intended ranks and well-separated nonzero spectra, consistent with Assumption (ref), throughout the simulation grid.
We provide the direct and recursive orthogonalized TWSF estimators the oracle ranks and lag length. This allows the simulation to focus on the finite-sample behavior predicted by Theorems (ref) and (ref) without introducing rank- and lag-selection errors. We evaluate the performance of these estimators via several metrics. For method $m \in $\{dir, rec\}, replication $s \in [R]$, and lead $\ell \in [h]$, denote the forecast error as $e_\ell^{(s,m)} \coloneqq \tilde{\theta}_\ell^{(s,m)} - \theta^{(s)}_\ell$. We report three error metrics:
which translate as the average terminal bias, terminal root-mean-squared-error (RMSE), and path RMSE. Terminal bias and terminal RMSE evaluate centering and accuracy at the same forecast point $T+h$ used for pointwise coverage. Path RMSE summarizes accuracy over the full $h$-period forecast trajectory. We do not use average path bias as a primary metric because signed errors can cancel across horizons, whereas path RMSE provides a more stable measure of full-trajectory performance. Finally, we evaluate terminal pointwise 90% coverage:
where $\widehat{\upsilon}_h^{(m)}$ is the plug-in standard error from the corresponding direct or recursive approaches defined as in (ref) and (ref). We also report the average length of these pointwise intervals,
We reemphasize that these confidence intervals are pointwise at $T+h$ and should not be interpreted as simultaneous confidence bands for the entire forecast path.
Figure (ref) reports the estimation error metrics as the dimension $d$ increases. The direct estimator is shown in blue and the recursive estimator in red, with shaded bands indicating one standard error computed from the nested design. For $h=1$, the two estimators perfectly match, as expected from the one-step equivalence. For $h = 5$ and $h=10$, the estimation errors generally decline as $d$ increases, consistent with the balanced fixed-rank scaling of the theory. The terminal bias is small and tends toward zero, while both terminal RMSE and path RMSE decrease with the panel dimension. For $h=1$, the two RMSE metrics coincide and hence is not reported separately.
Figure (ref) reports empirical pointwise coverage at the 90% nominal level. Coverage generally improves as $d$ increases and approaches the nominal target for the larger designs. The smallest designs exhibit noticeable undercoverage, especially for the direct estimator. This is consistent with the fact that the direct method uses longer Page blocks, leaving fewer effective temporal training samples in small designs. The recursive estimator is overall more stable in these simulations, particularly at smaller dimensions, reflecting its more parsimonious one-step training design. At larger dimensions, the two methods become more similar. Average interval lengths decrease with $d$, indicating that improved coverage is not driven by excessively wide intervals, but rather by improved centering and variance calibration.
Overall, the simulation results support the main theoretical predictions. In the balanced fixed-rank regime, estimation error decreases as the dimensions grow, and the plug-in pointwise confidence intervals become increasingly well calibrated. For longer horizons, the results illustrate the expected tradeoff between the direct and recursive strategies: the direct estimator learns horizon-specific forecasting rules but has fewer temporal training samples, whereas the recursive estimator is more parsimonious but propagates uncertainty through the companion recursion.
We revisit the NFL stadium-opening study of nfl_pnas, introduced in Section (ref). Whereas nfl_pnas studied the retrospective question of what would have happened to cities that opened their stadiums had they instead remained closed, we study a prospective version of the problem: what would have happened to a city that had kept its stadium closed if it had opened to fans at a future decision date?
During the 2020 NFL season, each team worked with local officials and public-health authorities to decide whether to admit fans or keep its stadium closed. These decisions were highly localized, reflecting local COVID-19 prevalence, public-health guidance, stadium logistics, and community risk tolerance. Following nfl_pnas, we use daily county-level COVID-19 case-count data from the New York Times NYT:21. The data are available beginning in April 2020, giving us a long pre-season period over which to learn cross-city relationships. Our observational units are the counties containing NFL stadiums; for readability, we refer to these units by their corresponding NFL city names. As in nfl_pnas, we omit Arizona because the relevant fan-origin county data were not available. After consolidating multi-team markets such as New York and Los Angeles, this leaves $N=29$ NFL cities. Of these, 18 admitted fans at some point during the 2020 regular season, while 11 kept their stadiums closed throughout the period we study. The first opening dates are given in Table (ref).
This adoption pattern is staggered, whereas the theoretical framework in Section (ref) assumes a simplified observation pattern with a common treatment date. We therefore adapt the empirical design to approximate the setting required by TWSF. This is similar in spirit to the empirical strategy used in synth_iv, where a staggered empirical setting is mapped into the framework of the estimator. A formal extension of TWSF to staggered treatment timing is an important direction for future work.
Throughout the empirical analysis, we use the six earliest open-stadium cities that first admitted fans during September as the treated donor pool: $\Ic_1 = $ \{Kansas City,Jacksonville,Cleveland,Dallas,Indianapolis,Miami\}, giving a modest $N_1 = 6$. As shown in Table (ref), these cities had all opened by September 20. We exclude Denver, even though it opened on September 27, because including it would shorten the window available for learning treated time dynamics. With this donor pool, we begin learning temporal relationships on September 21, the first day after all donor cities had admitted fans. For each target city and opening date $\tau$, we use observations through $\tau$ and forecast the following $h=14$ days. We include the opening date itself in the observed history because reported COVID-19 case counts on the day of a game cannot plausibly reflect infections generated by that same game. Thus, for a target city $i$, the horizon-specific estimand is $\theta_{i, \ell}(\tau) = \mathbb{E}[Y_{i, \tau+\ell}(1) \mid \Ec_{\tau+14}]$, $\ell \in [14]$. The 14-day horizon is chosen because any effect of a stadium opening on reported COVID-19 cases would be expected to appear only after a biological and reporting delay. PCR tests are more informative several days after exposure, symptoms often appear around day five, and later cases may arise up to roughly two weeks after exposure covid1, covid2. In all analyses, the unit-side weights $\widehat{\boldsymbol{\beta}}$ are learned from April 1 through September 10, the date of the first stadium opening. The time-side weights $\widehat{\boldsymbol{\alpha}}$ are learned from September 21 through the relevant opening date $\tau$. No case-count data during the forecast horizon $\tau+1, \dots, \tau+14$ are used in estimation.
Because the forecast horizon is $h=14$, the direct the direct multi-step estimator of Section (ref) is not feasible for the early target dates: valid Page blocks of length $K+h$ cannot be constructed with the available treated-donor time series. We therefore use the recursive orthogonalized TWSF estimator of Section (ref) and report pointwise 90% confidence intervals based on (ref).
The main counterfactual quantities of interest are not directly observable. We therefore begin with a validation exercise in which the target outcomes are eventually observed. The goal is to assess whether TWSF can forecast treated trajectories for cities that opened later in the season, using only information that would have been available before those cities opened. The validation targets are the treated cities that opened after the early donor pool $\Ic_1$. For each target city, we pretend that the observed panel ends on the day of its first open home game. For example, Carolina opened on October 4, so for Carolina we use data only through October 4 and forecast case counts from October 5 through October 18. The realized case counts over this two-week horizon are withheld during estimation and used only for validation. Although these observations are subject to idiosyncratic shocks and therefore imperfect, they remain a practical target for validation.
As aforementioned, we use the same donor pool $\Ic_1$ for all validation targets. Thus, the validation exercise evaluates whether the six earliest open-stadium cities contain enough cross-sectional and temporal information to forecast the later open-stadium cities. The length of the temporal training window varies across targets. For Carolina, Cincinnati, Houston, and Tampa Bay, the time-side model is learned from September 21 through October 4, giving the shortest treated-donor training window. For later targets such as Baltimore and Washington, the time-side model is learned over a substantially longer period. This variation provides an informal check on whether longer treated-donor histories improve the quality of the temporal forecast.
We select the algorithmic hyperparameters by cross-validation using root-mean-squared error. Table (ref) reports the selected lag length and ranks used in the recursive orthogonalized TWSF estimator. The validation results are visualized in Figure (ref). The observed case counts are shown as solid black curves, while the recursive orthogonalized TWSF forecasts are shown as dash-dotted red curves with pointwise 90% confidence intervals as the shaded bands. Overall, the forecasts track the realized trajectories closely. Out of the 11 validation targets, three cities---Carolina, Cincinnati, and Pittsburgh---have trajectories that are not fully contained within the pointwise confidence intervals. The remaining cities are well aligned with the observed trajectories. This provides encouraging evidence that the estimator can combine pre-season cross-city relationships with post-opening treated-donor dynamics to forecast future treated outcomes. The validation results also suggest that forecast quality improves as the treated-donor training window grows. Carolina and Cincinnati opened on October 4, giving the shortest window for learning $\widehat{\boldsymbol{\alpha}}$, and both are among the harder cases. Later-opening cities generally exhibit closer agreement between forecasts and observations. This pattern is consistent with the theoretical role of the temporal sample size in Proposition (ref) and Theorem (ref).
Crucially, this validation exercise should not be interpreted as a test of all identifying assumptions. It does not verify the latent factor model, rule out interference, or account for dynamic treatment histories, and thus all results should be interpreted with these limitations in mind. With that said, our validation study still offers a useful empirical diagnostic: when applied prospectively to cities that later opened their stadiums, the method produces forecasts that are generally closer to the realized treated trajectories.
We now turn to the main counterfactual exercise. For each of the 11 cities that kept their stadiums closed, we ask what would have happened had the city opened at a plausible decision date during the 2020 season. Specifically, for each control city, we set the hypothetical opening date as the first home game after October 1. This yields: October 4: \{Chicago, Detroit, Las Vegas, Los Angeles, San Francisco\}, October 11: \{New York, Seattle\}, October 18: \{Minnesota, New England\}, November 1:\{Buffalo, Green Bay\}. We then forecast the city’s open-stadium potential outcome path over the following 14 days.
For each city $i$ and hypothetical opening date $\tau$, we estimate $\widehat{Y}_{i, \tau+\ell}(1)$, $\ell \in [14]$, using only observations available through $\tau$. We compare this estimated open-stadium trajectory to the observed closed-stadium trajectory, $Y_{i, \tau+\ell} \coloneqq Y_{i, \tau+\ell}(0)$. Examining the forecast contrast, $\widehat{Y}_{i, \tau+\ell}(1) - Y_{i, \tau+\ell}$, sheds insight into the public-health impact of opening the stadium, whereby positive values suggest that opening is forecasted to produce more cases than remaining closed. We use the same donor pool $\Ic_1$ and recursive orthogonalized TWSF estimator as in the validation exercise of Section (ref) to produce our counterfactual forecasts. For each hypothetical opening date, we use the hyperparameters selected by cross-validation for the corresponding validation date in Table (ref).
The counterfactual results are shown in Figure (ref). For most control cities, the estimated open-stadium trajectory remains close to the observed closed-stadium trajectory, and the observed trajectory is generally contained within the pointwise 90% confidence intervals. Three cities---Buffalo, Green Bay, and Seattle---exhibit more noticeable deviations. Seattle is largely covered for most of the two-week horizon but deviates near the end. Green Bay is a known problematic case in nfl_pnas, which notes that a reporting error at the Wisconsin Department of Health Services generated a large idiosyncratic shock in the relevant pre-intervention period green_bay. Buffalo is also difficult to interpret because Erie County experienced a pronounced fall 2020 COVID-19 surge and became subject to New York’s micro-cluster restrictions at the start of November 2020 buffalo1, buffalo2.
Aside from these exceptions, the counterfactual estimates are broadly aligned with the observed closed-stadium trajectories. In fact, except for New England, the estimated open-stadium trajectories are generally below the observed closed-stadium trajectories. We do not interpret this as evidence that opening stadiums reduced cases. Rather, the more conservative interpretation is that the forecasts provide little evidence of a substantial increase in local case counts from opening stadiums under the mitigation protocols used during the 2020 NFL season. This conclusion is consistent with nfl_pnas, which found little evidence that opening stadiums to fans increased local case counts. One possible explanation is that strict mitigation and crowd-management protocols offset the risks associated with opening. Fans attending a controlled outdoor event with masking, distancing, capacity restrictions, and staff oversight may not have generated substantially more transmission than the alternative activities they would otherwise have pursued. At the same time, the counterfactual estimates vary across hypothetical opening dates. This is substantively plausible. The epidemic environment changed rapidly during Fall 2020, and the temporal model is trained on different treated-donor histories depending on the decision date. Thus, while the overall conclusion is that opening does not appear to induce a large increase in case counts, the timing of the opening decision can matter.
Several qualifications are important. First, the confidence intervals reported above are pointwise in the forecast horizon and should not be interpreted as simultaneous confidence bands for the full 14-day path. Second, the analysis treats stadium opening as binary and does not explicitly model attendance levels, differences in compliance or enforcement, or repeated home games. Finally, the findings are specific to the 2020 NFL season, outdoor stadium environments, and the mitigation protocols in place before the emergence of later, more transmissible variants.
With these caveats in mind, the case study illustrates a practical value of TWSF as a prospective causal forecasting tool. The method does not eliminate the usual concerns that accompany observational causal analysis, but it provides a disciplined way to combine cross-unit information, treated-donor time dynamics, and uncertainty quantification when decision-makers must evaluate a policy change before observing its consequences.
This article develops a framework for causal forecasting in panel data. The central problem is prospective rather than retrospective: for a control unit that has not yet received treatment, we seek to forecast its future potential outcome under treatment using only the panel observed up to the decision time. This setting falls between two classical paradigms: the SC and SI logic provides a foundation to estimate counterfactual outcomes within an observed post-treatment period, while multivariate time-series methods such as mSSA enable one to forecast future outcomes under the continuation of an observed regime. The framework developed here combines these ideas to address counterfactual regime switches beyond the observed panel. Our approach builds upon the latent factor structure underlying SI, and adds temporal dynamics through a low-rank Hankel representation of the treated latent time factors. This structure yields an identification result that expresses the causal forecast through two sets of weights: unit weights learned from the pre-treatment period \'a la SI and temporal forecasting weights learned from treated donor trajectories \'a la mSSA, culminating in the \texttt{TWSF} estimator.
We establish finite-sample pointwise error bounds and consistency for TWSF under suitable conditions, and then introduce an orthogonalized version of the estimator that enables asymptotic normality and valid pointwise inference. We further extend the framework to fixed multi-step forecasting horizons using both direct and recursive procedures, and show that analogous pointwise inferential guarantees continue to hold. The simulation results support the theoretical findings: estimation error generally decreases with the panel dimensions, and empirical coverage approaches the nominal level as the sample size grows. The NFL stadium-opening application illustrates the practical value of the method.
Several directions remain open. First, the theory currently relies on independent idiosyncratic noise, with Gaussian noise imposed for the inference results. Extending the analysis to allow serial dependence, cross-sectional dependence, heteroskedasticity, or heavier-tailed shocks would make the framework more broadly applicable. Second, the theoretical observation pattern assumes a common treatment date, whereas many empirical settings, including the NFL application, involve staggered adoption. Developing a formal staggered-adoption version of TWSF is therefore an important next step. Third, different temporal structures for the latent factors, such as stationary autoregressive dynamics of the kind studied in focus, may be preferable in some applications and could lead to alternative causal forecasting estimators. Fourth, the unit-side and time-side regressions could be adapted using alternative regression formulations. Finally, the current multi-step theory treats the forecast horizon as fixed and provides pointwise inference. Extending the results to growing horizons or to simultaneous confidence bands for the entire forecast path could further strengthen the framework for sequential policy decision-making.