EconBase
← Back to paper

The Harmonic Synthetic Control Method

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.

132,410 characters · 31 sections · 58 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

The Harmonic Synthetic Control Method

abstractSynthetic control methods can produce misleading counterfactual predictions when outcome series contain unit-specific stochastic trends, a common feature of nonstationary macroeconomic data. Existing remedies, such as pre-filtering or differencing, reduce spurious matching but may discard shared nonstationary variation that helps estimate donor weights. We propose Harmonic Synthetic Control (HSC), which replaces this binary choice with a soft allocation mechanism. HSC jointly estimates donor weights and a treated-unit-specific smooth residual component, then extrapolates this component into post-treatment periods using a time-series forecaster. A tuning parameter, selected by rolling-origin cross-validation, governs the division between donor matching and forecasting. As it varies, HSC continuously interpolates between synthetic control applied to differenced outcomes and synthetic control applied to raw outcomes with an intercept or trend. We provide a spectral interpretation showing how HSC downweights low-frequency residual components in donor matching and assigns them to the forecasting branch. A prediction-error decomposition separates weight-estimation distortion from residual-forecasting error. Monte Carlo exercises show that HSC adapts across regimes, performing well when stochastic trends are predominantly common or idiosyncratic, while estimators fixed to one regime can fail in the other. \noindentKeywords: synthetic control, nonstationarity, spurious regression, causal inference, frequency domain

\thispagestyle{empty} \doublespacing

\setcounter{page}{1} \abovedisplayskip=5pt \belowdisplayskip=5pt

Introduction

Synthetic control methods construct counterfactuals for treated units by finding weighted combinations of untreated donors that match the treated unit's pre-treatment outcome path abadie2003ecc,abadie2010synthetic, abadie2015comparative. The logic is that if a weighted combination of donors can reproduce the treated unit's outcomes before treatment, the same combination should approximate what the treated unit's outcomes would have been in the absence of treatment. A difficulty arises when outcome series are nonstationary: the pre-treatment fit that synthetic control exploits may be spurious. Specifically, when units contain unit-specific stochastic trends, a convex combination of donors can closely track the treated unit's pre-treatment path through coincidental co-movement rather than shared structure, producing in-sample fit that breaks down out of sample and leading to biased estimation and distorted inference phillips1986understanding,masini2021counterfactual, masini2022counterfactual,shi2025synthetic.

Existing responses to this problem face a tradeoff over whether to preprocess the data, such as detrending or differencing. These transformations turn nonstationary time series into stationary ones, thereby alleviating the spurious matching risk. However, these transformations also discard nonstationary variation that the treated unit potentially shares with donor units, which is the main source of identifying variation for donor-weight estimation in synthetic control ferman2021synthetic,abadie2021using. Other variants of synthetic control, such as augmented or bias-corrected extensions, can reduce the residual imbalance left by imperfect pre-treatment fit ben2021augmented,arkhangelsky2021synthetic, but because they still begin from weights chosen to match raw pre-treatment outcomes, they can inherit the same underlying weight distortion caused by spurious matching.

We formalize this tradeoff through a conceptual distinction. In macroeconomic panel data, nonstationarity typically takes the form of persistent stochastic components, such as random walks or other integrated processes, whose variance grows without bound over time. We categorize this stochastic trend variation into two sources. By a shared stochastic trend, we mean a stochastic trend component whose innovations are shared across units; such a trend moves the treated unit and donors together, possibly with unit-specific responses. Synthetic control is designed to exploit this shared structure, using a weighted combination of donors to approximate the treated unit's trajectory. By an idiosyncratic stochastic trend, we mean a stochastic trend component whose realizations are unit-specific and do not generate stable comovement across units, even though they may appear correlated by chance in any finite sample. This is the source of the spurious matching problem. We use stochastic trend as an umbrella term for both.

figure[figure omitted — 1,001 chars of source]

Based on these concepts, we decompose untreated potential outcomes of the treated unit and all control donors into three components (Figure (ref)). The first is a component $L$, possibly low-rank and driven by shared latent factors whose loadings vary across units. These factors may include both stochastic trend and short-run variations.\footnote{We avoid the labels stationary and nonstationary because what matters for our analysis is bounded long-run variance, not strict stationarity. By short-run we mean variation with bounded long-run variance, which includes strictly stationary processes as well as bounded-variance nonstationary processes.} The second is idiosyncratic short-run noise $\varepsilon$. It does not fundamentally distort donor-weight estimation: it averages out over long pre-treatment windows under standard moment conditions. The third is the idiosyncratic stochastic trend $R$, which is not governed by the shared factor structure. Unlike $\varepsilon$, $R$ can severely distort donor-weight estimation: independent stochastic trends produce realized correlations whose magnitude does not vanish as $T_0$ grows, so longer pre-treatment windows do not help phillips1986understanding. Whether the stochastic trend variation in the treated unit's series is mostly shared(captured by $L$) or idiosyncratic (captured by $R$) is generally unknown to the researcher ex ante. In the first case, preprocessing removes the variation that synthetic control would otherwise use for donor-weight estimation. In the second case, allowing stochastic trends to enter weight estimation creates the risk of spurious matching. This tradeoff motivates us to propose a solution that adapts to both regimes.

This paper proposes harmonic synthetic control (HSC), which replaces this binary choice of preprocessing with a soft, data-driven allocation. HSC jointly estimates donor weights and a treated-unit-specific smooth component $E$ that absorbs idiosyncratic trend-like variation not reproducible by a convex combination of donors. The roughness of $E$ is controlled by a penalty $\|D_q E\|_2^2$, where $D_q$ is the $q$th-order difference operator, with $q \in \{1, 2\}$ specifying the order of smoothness. The division of labor between donor matching and the smooth component is governed by a single tuning parameter $\rho \in [0,1]$, selected by rolling-origin cross-validation. When $\rho$ is close to $1$, HSC imposes a high penalty on the roughness of $E$ and recovers synthetic control with an intercept or with an intercept plus a linear trend, depending on whether $q=1$ or $q=2$. When $\rho$ is close to $0$, $E$ absorbs the entire discrepancy between the treated unit and the convex combination of donors, and the weight estimation in HSC approaches synthetic control applied to first or second differences of the raw outcomes (for $q=1$ or $q=2$, respectively). At intermediate $\rho$, HSC continuously interpolates between these two endpoints. In post-treatment periods, the smooth component is forecast by a time series forecaster, and the counterfactual is constructed by adding the forecast of $E$ to the donor matching component. HSC does not attempt to disentangle which portion of the stochastic trend variation is shared and which is idiosyncratic; instead, cross-validation selects the allocation between donor matching and the smooth component that yields the best out-of-sample predictive performance.

This soft allocation has a clean spectral interpretation. Any pre-treatment series can be decomposed into components at different frequencies: low-frequency content carries slowly varying (trend-like) variation, while high-frequency content carries short-run variation. The HSC weight estimation problem down-weights low-frequency components and amplifies high-frequency components, thus alleviating the spurious matching risk. The treated-unit-specific component $E$ absorbs the low-frequency residual variation that remains after donor matching, and the post-treatment forecast of $E$ extrapolates this low-frequency, trend-like component forward in time. Equivalently, HSC can be understood as applying synthetic control to a soft spectral transformation of the raw data that interpolates between two extremes: at $\rho = 0$, the transformation reduces to $q$-th order differencing, which strongly suppresses low-frequency content; at $\rho = 1$, the transformation removes only the null-space component (constants for $q=1$, constants plus linear trends for $q=2$), leaving other frequencies unchanged. The method's name reflects the form of this interpolation: at each frequency, the spectral gain of the HSC transformation is a weighted harmonic mean of the gains at the two endpoints, with weights $1-\rho$ and $\rho$.

We develop an envelope bound on the HSC counterfactual's prediction error. The error decomposes into a weight-estimation term and a forecasting term, each depending on $\rho$. The weight-estimation term reflects that the researcher observes $Y = L + R + \varepsilon$ rather than the shared component $L$ alone; the HSC-estimated donor weights therefore deviate from the oracle weights constructed from $L$ alone. This term captures the tradeoff between the risk of spurious matching at large $\rho$ and the downweighting of useful low-frequency variation in $L$ at small $\rho$. The forecasting term captures the error that would remain even with oracle weights; depending on the prediction quality of the time series forecaster, this term can be monotonically increasing, monotonically decreasing, or non-monotonic in $\rho$. Together these two terms determine the tradeoff that cross-validation aims to balance in the choice of $\rho$.

The synthetic control literature has produced many methods that match donors to the treated unit's pre-treatment outcomes, differing in weight constraints, bias-correction strategies, and temporal aggregation doudchenko2016balancing,ben2021augmented,arkhangelsky2021synthetic,sun2024temporal. A growing subset of this literature addresses synthetic control specifically under nonstationarity. masini2021counterfactual,masini2022counterfactual show that counterfactual estimation requires a cointegrating relationship between treated and donor units; without it, estimated effects diverge and inference suffers severe size distortion. harvey2021cointegration model the common stochastic trend explicitly and propose stationarity tests on the pre-treatment difference as a diagnostic for donor selection. shi2025synthetic decompose outcomes into trend and cycle via the Hamilton filter and restrict donor matching to the cyclical residual. These contributions diagnose the nonstationarity problem or propose hard filters that remove persistent variation before matching. HSC introduces the smooth component $E$ as a new degree of freedom; rather than diagnosing nonstationarity or applying a fixed filter, HSC offers a soft, data-driven alternative that retains shared stochastic trend while mitigating idiosyncratic spurious matching.

The estimator's mechanism draws on a different technical toolkit than is standard in the synthetic control literature. The roughness penalty $\|D_q E\|_2^2$ connects HSC to the Whittaker--Henderson smoothing framework whittaker1922new,henderson1924new, the Hodrick--Prescott filter hodrick1997postwar, and penalized spline formulations eilers1996flexible. The spectral decomposition of the HSC metric, in which the tuning parameter $\rho$ acts as a frequency-dependent gain function on the pre-treatment residual, imports ideas from spectral analysis in time series into the synthetic control setting. This connection provides both interpretive clarity ($\rho$ governs a soft spectral partition between what is matched cross-sectionally and what is forecasted univariately) and computational tractability, since the profiled HSC objective reduces to a standard constrained quadratic program.

We illustrate the method on the canonical Hong Kong example, the path of per-capita GDP after the 1997 return of Hong Kong to Chinese sovereignty, studied by hsiao2012panel and shi2025synthetic. Cross-validation selects an interior allocation rather than either preprocessing extreme, confirming that the data prefer a soft partition between donor matching and the smooth component. HSC distributes donor weight broadly across the control pool, whereas level-matching and filter-based competitors either concentrate weight on a few donors or extrapolate the treated unit's own trend and overshoot the observed series. On a rolling-origin out-of-sample criterion HSC is the most accurate estimator among those we compare, and this ranking is stable across the cross-validation horizon and across two very different donor-selection philosophies. The application thus reproduces, in real data, the soft-allocation behavior that the theory and the Monte Carlo evidence predict.

The remainder of the paper is organized as follows. Section (ref) formalizes the outcome decomposition, establishes notation, and characterizes the two failure modes (spurious donor matching and over-filtering) that motivate the need for a soft allocation mechanism. Section (ref) introduces the HSC estimator, derives its profiled representation, and constructs the forecast operator that extrapolates the smooth component into post-treatment periods. Section (ref) develops the spectral interpretation, showing that the tuning parameter $\rho$ acts as a frequency-dependent gain function, and describes the cross-validation procedure for selecting $\rho$. Section (ref) presents the prediction-error decomposition into weight-estimation and forecasting terms, and discusses the trade-off under the selection of $\rho$. Section (ref) reports Monte Carlo evidence for the HSC estimator. Section (ref) applies HSC to the 1997 Hong Kong handover and compares it with established alternatives. Section (ref) concludes.

A Tradeoff between Spurious Donor Matching and Over-Filtering

This section formalizes the allocation problem described in the Introduction. We first establish notation for the synthetic control setting, then provide the formal $L + R + \varepsilon$ decomposition and characterize two failure modes: spurious donor matching when the idiosyncratic stochastic trend $R$ dominates, and over-filtering when shared stochastic trend in $L$ is discarded. A simulated illustration shows that existing methods commit to one regime or the other.

Setup and notation

We observe a panel of $N_0+1$ units over $T=T_0+T_{\mathrm{post}}$ periods. Unit $i=1$ is the treated unit; units $i=2,\ldots,N_0+1$ form the donor pool. Treatment is imposed at the end of period $T_0$, so the pre-treatment window is $t=1,\ldots,T_0$ and the post-treatment window is $t=T_0+1,\ldots,T$. Let $Y_{it}(0)$ denote the untreated potential outcome for unit $i$ at time $t$. Under the standard no-anticipation assumption, we observe $Y_{1t}=Y_{1t}(0)$ for $t\le T_0$. Let $X_t=(Y_{2t},\ldots,Y_{N_0+1,t})'$ denote the $N_0\times 1$ vector of donor outcomes at time $t$. We write $Y_{\mathrm{pre}}=(Y_{1,1},\ldots,Y_{1,T_0})'$ and $X_{\mathrm{pre}}=(X_1,\ldots,X_{T_0})'$ for the $T_0\times 1$ and $T_0\times N_0$ pre-treatment arrays. $Y_{\mathrm{post}}$ and $X_{\mathrm{post}}$ are defined similarly.

The synthetic control estimator constructs a counterfactual for the treated unit as a weighted combination of donors. The weights $\hat\omega$ are chosen from the simplex $\Delta=\{\omega\in\mathbb{R}^{N_0}:\omega\ge 0,\;\mathbf{1}'\omega=1\}$ by minimizing the pre-treatment sum of squared residuals:

equation*[equation* omitted — 120 chars of source]

and the counterfactual at horizon $h\ge 1$ is $\hat{Y}_{1,T_0+h}(0) = X_{T_0+h}'\hat\omega$. A common simplex-constrained extension adds an intercept $\hat\alpha$ to absorb a constant level shift between the treated unit and the weighted donors. The weights and intercept are estimated jointly:

equation*[equation* omitted — 180 chars of source]

which is equivalent to matching on demeaned pre-treatment outcomes doudchenko2016balancing,ferman2021synthetic. The counterfactual becomes $\hat{Y}_{1,T_0+h}(0) = X_{T_0+h}'\hat\omega + \hat\alpha$.

Both formulations share the same basic structure: donor weights are chosen to minimize a pre-treatment loss computed on the observed outcome series, and the counterfactual extrapolates those weights into the post-treatment window.

Decomposing untreated outcomes and the allocation problem

We now formalize the decomposition introduced informally in Section (ref) (Figure (ref)). We decompose untreated potential outcomes into three components:

equation[equation omitted — 75 chars of source]

Here $L_{it} = \Lambda_i'F_t$ is a shared low-rank component, where $F_t$ collects latent factors and $\Lambda_i$ is the corresponding vector of unit-specific loadings. The factors $F_t$ may contain both short-run and stochastic trend movements. The defining property of $L_{it}$ is that it is shared: when a convex combination of donors matches the treated unit's loadings $\Lambda_1$, it reproduces $L_{1t}$ in every period. By contrast, $R_{it}$ and $\varepsilon_{it}$ are idiosyncratic. The term $R_{it}$ denotes an idiosyncratic stochastic trend, whose long-run variance grows without bound, whereas $\varepsilon_{it}$ denotes idiosyncratic short-run noise with bounded long-run variance. Unlike the shared component, neither $R_{it}$ nor $\varepsilon_{it}$ is governed by a shared factor structure. As a result, observed outcomes may mask the low-rank structure in $L_{it}$, with the main difficulty arising from the idiosyncratic stochastic trend component $R_{it}$.

This decomposition refines the outcome framework of arkhangelsky2021synthetic by separating idiosyncratic stochastic trends from short-run idiosyncratic noise. In arkhangelsky2021synthetic's notation, untreated outcomes are represented by a systematic component $\mathcal{L}_{it}$ and an idiosyncratic error $\epsilon_{it}$. Their asymptotic analysis allows $\mathcal{L}$ to be an approximately low-rank systematic matrix, without imposing a fixed known rank, while their Assumption 1 restricts the rows $\epsilon_{i\cdot}$ of the noise matrix to be i.i.d.\ Gaussian vectors with covariance matrix $\Sigma\in\mathbb R^{T\times T}$ whose eigenvalues are bounded and bounded away from zero. This condition permits temporal dependence and some forms of nonstationary covariance heterogeneity. In the growing-$T$ asymptotics under which Assumption 1 is imposed, however, it excludes integrated unit-specific stochastic trends: if $\epsilon_{it}$ were a random walk with innovation variance $\sigma^2$, then $\operatorname{Cov}(\epsilon_{is},\epsilon_{it})=\sigma^2\min\{s,t\}$, so the largest eigenvalue of $\Sigma$ grows on the order of $T^2$ and is not bounded uniformly in $T$, violating the bounded-eigenvalue requirement of Assumption 1. Our decomposition therefore writes the idiosyncratic component as $R_{it}+\varepsilon_{it}$, where $R_{it}$ captures the unit-specific stochastic trend excluded by the bounded-eigenvalue condition and $\varepsilon_{it}$ denotes the remaining short-run noise. This split makes explicit the type of persistent idiosyncratic variation that we argue is central to the spurious matching problem.

A parallel restriction appears in ferman2021synthetic, who study synthetic control in the fixed-$N$, large-$T_0$ regime under a linear factor model with common factors, unit-specific loadings, and idiosyncratic shocks $\epsilon_{it}$. Their Assumption 4 requires the pre-treatment moments to converge in probability to non-stochastic constants: $T_0^{-1}\sum_t \epsilon_{it}^2 \xrightarrow{p} \sigma_\epsilon^2$, with analogous conditions on the common factors and on cross-products of factors and noise. The substantive content matches arkhangelsky2021synthetic's Assumption 1: the idiosyncratic noise must have bounded long-run variance with well-behaved pre-treatment sample moments. An idiosyncratic random walk in $\epsilon_{it}$ violates Assumption 4 in the same way it violates the bounded-eigenvalue condition.

As discussed in Section (ref), the fitting criterion in Section (ref) operates on observed outcomes and therefore does not distinguish among $L_{it}$, $R_{it}$, and $\varepsilon_{it}$. The resulting allocation problem, deciding how much stochastic trend variation to attribute to the shared component $L_{it}$ versus the idiosyncratic stochastic trend $R_{it}$, leads to two failure modes formalized in the following subsections.

The risk of spurious donor matching

When $R_{it}$ contributes substantially to $Y_{it}(0)$, the pre-treatment optimization can achieve a close fit by assigning weight to donors whose independent persistent movements happen to co-move with the treated unit over the observed window. This is a manifestation of spurious regression in the synthetic control setting: the in-sample fit may appear excellent, yet it is driven by coincidental trending behavior rather than a shared factor structure, and therefore does not extend beyond the pre-treatment period. shi2025synthetic formalize this as the spurious synthetic control problem, noting that “even if a country's GDP can be closely approximated by a weighted average of others over a given period, such a fit may arise purely from coincidental trending behavior.” The simplex constraints $\omega\ge 0$, $\mathbf{1}'\omega=1$ narrow the feasible set but do not prevent spurious fit: independent random walks can still produce close pre-treatment matches within the simplex.

It is worth emphasizing that the central failure is weight distortion, not merely imperfect fit. When $R_{it}$ is quantitatively important, the optimizer is drawn toward donors whose idiosyncratic trends happen to track $R_{1t}$ in the pre-treatment period, pulling the weights away from the combination that would best reproduce the shared component $L_{it}$. This distinction between weight distortion and lack of fit is important because several recent proposals, including augmented synthetic control ben2021augmented and synthetic difference-in-differences arkhangelsky2021synthetic, improve on the standard synthetic control by augmenting it with an outcome model that corrects for the residual imbalance left by imperfect pre-treatment fit. Because these methods still begin from the donor weights chosen to match pre-treatment outcome levels, they can inherit the same underlying weight-distortion channel: the bias-correction step operates conditional on already-distorted weights, so it can at best mitigate but does not fully undo this distortion.

More fundamentally, the existing literature on synthetic control with nonstationary data identifies cointegration between the treated unit and the synthetic control as the key condition separating valid donor matching from spurious fit. In the regression-based framework studied by masini2021counterfactual, counterfactual estimation requires the data-generating process to admit a cointegrating relationship involving the treated unit; without such a relationship, the regression is spurious. masini2022counterfactual show that in this case the estimated treatment effect diverges and that ignoring the nonstationary nature of the data leads to severe over-rejection of the null hypothesis of no effect. harvey2021cointegration reach a similar conclusion from a structural time series perspective: they model the shared stochastic trend explicitly and argue that a synthetic control is valid when the target and control series share a common stochastic trend, proposing stationarity tests on the pre-treatment difference as a diagnostic for donor selection.

In the language of Section (ref), cointegration between the treated unit and the synthetic control therefore requires that the stochastic trend content of their difference be eliminated not only in the pre-treatment period but also out of sample. When the stochastic trend variation is entirely driven by shared factors in $L_{it}$, a cointegrating relationship among the units is guaranteed by the factor structure, and the donor weights that minimize the pre-treatment criterion can recover the common component both in pre- and post-treatment periods. When the treated unit also contains a quantitatively important idiosyncratic stochastic trend $R_{1t}$, the cointegrating relationship involving the treated unit may not exist, and the critiques from the cited literature then apply directly.

The cost of hard filtering

The previous subsection showed that level-based matching is vulnerable to spurious donor matching when $R_{it}$ is large. Two natural responses have been proposed, both of which remove stochastic trend variation before constructing the synthetic control.

The first is explicit pre-filtering. shi2025synthetic propose decomposing the outcome into a trend and a cyclical component using the Hamilton filter hamilton2018you, forecasting the treated unit's trend from its own lagged values, and restricting donor matching to the stationary cyclical residual. Their strategy is designed to eliminate spurious matching from idiosyncratic stochastic trend by construction under the maintained trend-cycle decomposition: the Hamilton filter removes the trend component, and donor matching is then applied to the resulting cyclical residual rather than to the raw persistent series.

The second is differencing. If $R_{it}$ is a random walk, first-differencing yields $\Delta Y_{it}=\Delta L_{it}+\Delta R_{it}+\Delta\varepsilon_{it}$, where $\Delta R_{it}$ is now stationary. This logic underlies the concern raised by masini2022counterfactual that nonstationary data can produce spurious counterfactual estimates, motivating transformations such as differencing before applying counterfactual methods, and is related to abadie2021using's (abadie2021using) observation that matching in differences can help when levels are not credibly matchable.

Both strategies address the spurious matching problem, but they share a common cost: neither distinguishes between $R_{it}$ and the stochastic trend part of $L_{it}$. If the common factor $F_t$ is a random walk with heterogeneous loadings $\Lambda_i$, then first-differencing removes the resulting stochastic trend in $\Lambda_i'F_t$ entirely. ferman2021synthetic show that when common factors include diverging nonstationary components, these components dominate the pre-treatment fitting criterion as $T_0$ grows, providing the primary source of identifying variation for donor weight estimation. Mechanically removing them discards an important signal for donor matching. What remains after filtering or differencing is a stationary factor structure with potentially much lower signal-to-noise ratio. The resulting donor weights can become substantially more imprecise and perform poorly in recovering the shared component of $L_{it}$. abadie2021using indeed notes that matching in first differences can inflate the variance attributable to noise, potentially inducing an increase in bias.

Illustration of the tradeoff

The preceding two subsections identified a tension between two failure modes: spurious donor matching when $R_{it}$ dominates, and over-filtering when shared stochastic trend variation in $L_{it}$ is discarded. The relative severity of these two risks depends on the magnitude of $R_{it}$ relative to $L_{it}$.

Figure (ref) illustrates the tradeoff with a simulated example. The data-generating process is $Y_{it}(0)=\Lambda_i' F_t+\kappa\cdot R_{it}+\varepsilon_{it}$, where $F_t$ is a shared random-walk factor with heterogeneous loadings ($\Lambda_i \sim N(1,1)$), $R_{it}$ are idiosyncratic random walks, and $\varepsilon_{it}$ is iid noise with $\sigma_\varepsilon=2$. The common factor structure is the same across the two rows; the only difference is the importance of the idiosyncratic stochastic trend, governed by $\kappa$. The top row sets $\kappa=0$ (only shared stochastic trend); the bottom row sets $\kappa=2$ (an idiosyncratic stochastic trend added on top of the same common structure). Each column applies a different estimator. The first is synthetic control in levels with an intercept, which constructs the counterfactual as $\hat{Y}_{1,T_0+h}(0) = X_{T_0+h}'\hat\omega^{\mathrm{lev}} + \hat\alpha$, where $\hat\alpha$ absorbs a constant shift between the treated unit and the weighted donors. The second is synthetic control in first differences, which estimates donor weights on the differenced data $\Delta Y_{it}$ and anchors the counterfactual at the last pre-treatment observation: $\hat{Y}_{1,T_0+h}(0) = X_{T_0+h}'\hat\omega^{\mathrm{dif}} + (Y_{1,T_0} - X_{T_0}'\hat\omega^{\mathrm{dif}})$.

figure[figure omitted — 970 chars of source]

In the top-left panel, synthetic control with intercept captures the common factor structure and tracks the treated unit closely in both the pre- and post-treatment periods (RMSE $=2.2$). In the top-right panel, first-differencing weakens the dominant common signal; matching on differenced data in the presence of large stationary noise ($\sigma_\varepsilon=2$) produces weights that recover the common structure less precisely, and the reconstructed level path has a larger discrepancy (RMSE $=6.0$). In the bottom-left panel, synthetic control with intercept shows a large performance discrepancy between pre- and post- treatment periods, suggesting an instance of spurious donor matching (RMSE $=15.3$). In the bottom-right panel, first-differencing synthetic control avoids the spurious match and produces a more accurate counterfactual (RMSE $=4.4$).

In practice, the researcher does not know how important the idiosyncratic stochastic trend component is relative to the shared component. A method that commits fully to either level matching or hard filtering will fail in one regime or the other. This motivates the need for a soft allocation mechanism rather than a binary choice. The goal of this paper is to provide a synthetic control estimator that lets the data determine how much of the stochastic trend variation in the treated unit should be allocated between donor matching and a smooth treated-specific component.

Harmonic Synthetic Control

Section (ref) showed that synthetic control's pre-treatment fitting criterion does not distinguish between the shared low-rank component $L_{it}$ and the idiosyncratic stochastic trend $R_{it}$; level matching and hard filtering each fail in one regime or the other. To bridge these two regimes, we propose the harmonic synthetic control (HSC) estimator. The construction proceeds in three steps: weight estimation for donor matching, the treated-unit-specific time series forecaster, and the construction of the counterfactual. HSC does not attempt to disentangle $L_{it}$ and $R_{it}$ in the data-generating process. Instead, it seeks the optimal allocation between modeling stochastic trend variation as shared low-rank structure and modeling it as idiosyncratic stochastic trends. A single tuning parameter $\rho \in [0,1]$ governs such an allocation.

Donor weights of HSC

We use $q\in\{1,2\}$ to denote the smoothness order. The difference operator $D_q$ is the $q$th-order difference operator: $D_1$ is the $(T_0-1)\times T_0$ first-difference operator with rows $(D_1 x)_t = x_{t+1}-x_t$, and $D_2$ is the $(T_0-2)\times T_0$ second-difference operator with rows $(D_2 x)_t =x_{t+2}-2x_{t+1}+x_t$.\footnote{Explicitly, $D_1$ has $-1$ on the main diagonal and $+1$ on the first superdiagonal; $D_2$ has entries $+1, -2, +1$ on three consecutive diagonals.} Let $K_q := D_q'D_q$, which is symmetric positive semidefinite with null space $\mathrm{Null}(K_q)$ equal to $\mathrm{span}\{\mathbf{1}_{T_0}\}$ for $q=1$ and $\mathrm{span}\{\mathbf{1}_{T_0}, t_\mathrm{pre}\}$ for $q=2$, where $t_\mathrm{pre}:=(1,\dots,T_0)'$. Recall that $X_\mathrm{pre}\in\ensuremath{\mathbb{R}}^{T_0\times N_0}$ denotes the pre-treatment donor outcome matrix, $Y_{\mathrm{pre}}\in\ensuremath{\mathbb{R}}^{T_0}$ is the treated unit's pre-treatment outcome vector, and $\Delta_{N_0}:=\{\omega\in\ensuremath{\mathbb{R}}^{N_0}:\omega\ge 0,\;\mathbf{1}'\omega=1\}$ is the unit simplex.

definition[Harmonic Synthetic Control, $\rho\in(0,1)$] For $\rho\in(0,1)$, the HSC estimator jointly solves \begin{equation} \bigl(\hat\omega(\rho,q),\,\hat E_\mathrm{pre}(\rho,q)\bigr) \;\in\; \mathop{\rm argmin}_{\omega\in\Delta_{N_0},\; E\in\ensuremath{\mathbb{R}}^{T_0}} \left\{ \frac{1}{\rho}\bigl\|Y_{\mathrm{pre}}-X_\mathrm{pre}\omega - E\bigr\|_2^2 \;+\;\frac{1}{1-\rho}\|D_q E\|_2^2 \;+\;\zeta^2 T_0\|\omega\|_2^2 \right\}, \end{equation} where $\zeta>0$ is a ridge regularization parameter.\footnote{Following arkhangelsky2021synthetic, we set the default value of $\zeta$ as $T_\mathrm{post}^{1/4}\hat\sigma$, where $\hat\sigma$ is the standard deviation of all elements of the first-differenced donor matrix $D_1 X_\mathrm{pre}$. This is the single-treated-unit specialization of the formula $\zeta = (N_{\mathrm{tr}}T_\mathrm{post})^{1/4}\hat\sigma$ used for multiple treated units.} The first term penalizes the discrepancy between the treated unit and the donor-weighted combination after removing the latent smooth component $E$; the second penalizes the roughness of $E$ through the $q$th-difference penalty $\|D_q E\|_2^2 = E'K_q E$; the third is a ridge penalty that helps stabilize the donor weights.

$D_q$ determines what counts as a smooth component. When $q=1$, the penalty $\|D_1E\|_2^2=\sum_{t=1}^{T_0-1}(E_{t+1}-E_t)^2$ penalizes local changes in level, so smoother paths are those that vary less from one period to the next. Constant vectors, which are in the null space of $K_1$, receive no penalty. When $q=2$, the penalty $\|D_2E\|_2^2=\sum_{t=1}^{T_0-2}(E_{t+2}-2E_{t+1}+E_t)^2$ penalizes local changes in slope, so smoother paths are those with less curvature. Intercept-plus-linear-trend components, which are in the null space of $K_2$, receive no penalty. Thus, for $q=1$ the smooth component is encouraged to be locally flat in levels, whereas for $q=2$ it is encouraged to be locally linear in time.

The tuning parameter $\rho$ governs the relative incentives in the joint optimization over $(\omega,E)$. Conditional on donor weights $\omega$, the optimizer chooses $E$ to balance fidelity to the residual $Y_{\mathrm{pre}}-X_\mathrm{pre}\omega$ against the smoothness restriction imposed by $\|D_qE\|_2^2$. When $\rho$ is small, the fit term $\frac{1}{\rho}\|Y_{\mathrm{pre}}-X_\mathrm{pre}\omega-E\|_2^2$ receives relatively large weight compared with the roughness penalty. Therefore, conditional on a given $\omega$, the optimizer is more willing to let $E$ track a larger and potentially rougher portion of the residual, where roughness is measured by $D_q$. When $\rho$ is large, the roughness penalty $\frac{1}{1-\rho}\|D_qE\|_2^2$ receives relatively large weight, so conditional on a given $\omega$, the optimizer is forced to choose a smoother $E$, where smoothness is in the sense measured by $D_q$.

Because $\omega$ and $E$ are chosen jointly, the tuning parameter $\rho$ does not act on a fixed pre-treatment discrepancy. Instead, it determines how the joint optimizer splits the treated series $Y_\mathrm{pre}$ between the donor-matched component $X_\mathrm{pre}\hat\omega$ and the treated-unit-specific smooth component $\hat{E}$, with larger $\rho$ forcing $\hat{E}$ to be smoother.

Profiling and the HSC metric

The joint formulation in Definition (ref) is useful for intuition, but the estimator becomes more transparent and computationally tractable after profiling out the treated-unit-specific smooth component \(E\). This yields an equivalent weight-estimation problem in which the pre-treatment residual is measured under a \(\rho\)- and \(q\)-dependent metric.\footnote{We use the term “metric” informally throughout: \(W_{\rho,q}\) is symmetric positive semidefinite (not positive definite), and it annihilates \(\mathrm{Null}(K_q)\), so the quadratic form \(r'W_{\rho,q}r\) is a seminorm rather than a norm. Residual components in \(\mathrm{Null}(K_q)\) contribute zero to the HSC criterion and are handled separately through the smooth component \(\hat E_\mathrm{pre}\); see the remark after Proposition (ref).}

For any candidate donor weight vector \(\omega\in\Delta_{N_0}\), define the pre-treatment discrepancy/residual as

equation[equation omitted — 95 chars of source]

Fix \(\rho\in(0,1)\). For each \(\omega\), the inner minimization over \(E\) in (ref) is a strictly convex quadratic program: \[ \min_{E\in\ensuremath{\mathbb{R}}^{T_0}} \left\{ \frac{1}{\rho}\|r_\mathrm{pre}(\omega)-E\|_2^2 +\frac{1}{1-\rho}\|D_qE\|_2^2 \right\}. \] To express its solution, let

equation[equation omitted — 123 chars of source]
proposition[Profiled representation] Fix \(q\in\{1,2\}\) and \(\rho\in(0,1)\). For every \(\omega\in\Delta_{N_0}\), the inner problem in \(E\) has the unique minimizer \begin{equation} \hat E_\mathrm{pre}(\omega;\rho,q)=S_{\rho,q}\,r_\mathrm{pre}(\omega). \end{equation} Substituting this optimizer back into the criterion yields the profiled objective \begin{equation} \hat\omega(\rho,q) \;\in\; \mathop{\rm argmin}_{\omega\in\Delta_{N_0}} \left\{ r_\mathrm{pre}(\omega)'W_{\rho,q}r_\mathrm{pre}(\omega) +\zeta^2T_0\|\omega\|_2^2 \right\}, \end{equation} where \begin{equation} W_{\rho,q}:=\frac{1}{\rho}(I_{T_0}-S_{\rho,q}), \qquad \rho\in(0,1). \end{equation} Equivalently, the fitted treated-unit-specific smooth component is \begin{equation} \hat E_\mathrm{pre}(\rho,q) = S_{\rho,q}\bigl(Y_{\mathrm{pre}}-X_\mathrm{pre}\hat\omega(\rho,q)\bigr). \end{equation}

The proof is deferred to Appendix (ref). Proposition (ref) shows that donor weight estimation in HSC can be understood through the metric $W_{\rho,q}$, which re-weights the pre-treatment discrepancy $r_\mathrm{pre}$ in a standard ridge-regularized quadratic program on the simplex. The operator $S_{\rho,q}$ acts as a smoother that extracts the smooth part of $\hat r_\mathrm{pre}$ as the smooth component $\hat E_\mathrm{pre}$.

We now show that this family of metrics extends continuously to the boundary cases $\rho=0$ and $\rho=1$. Define

align[align omitted — 154 chars of source]

where $P_{0,q}$ denotes the orthogonal projector onto $\mathrm{Null}(K_q)$.

proposition[Continuous extension of the HSC metric] For each fixed \(q\in\{1,2\}\): \begin{enumerate}[label=(\roman*)] • The interior families satisfy \[ \lim_{\rho\downarrow 0}S_{\rho,q}=I_{T_0}, \qquad \lim_{\rho\uparrow 1}S_{\rho,q}=P_{0,q}, \] and \[ \lim_{\rho\downarrow 0}W_{\rho,q}=K_q, \qquad \lim_{\rho\uparrow 1}W_{\rho,q}=I_{T_0}-P_{0,q}. \] Hence the boundary definitions (ref) and (ref) are the unique continuous extensions of \(\{S_{\rho,q}\}\) and \(\{W_{\rho,q}\}\) from \((0,1)\) to \([0,1]\). • \(W_{\rho,q}\) is symmetric positive semidefinite for every \(\rho\in[0,1]\). \end{enumerate}

The proof is deferred to Appendix (ref). The key implication of Proposition (ref) is that the HSC weight problem extends continuously from the interior \(\rho\in(0,1)\) to the endpoint cases \(\rho=0\) and \(\rho=1\), where the two boundary values correspond to familiar special cases.

At \(\rho=0\), \(S_{0,q}=I_{T_0}\) means the smoother assigns the entire residual to the treated-unit-specific “smooth” component; correspondingly, \(W_{0,q}=K_q=D_q'D_q\), so the metric measures only the \(q\)th-order roughness of the residual: \[ r_\mathrm{pre}(\omega)'W_{0,q}r_\mathrm{pre}(\omega) = \|D_qr_\mathrm{pre}(\omega)\|_2^2. \] Hence the profiled objective becomes

equation[equation omitted — 206 chars of source]

that is, synthetic control applied to the \(q\)th-differenced outcomes with a ridge penalty.

At \(\rho=1\), because \(W_{1,q}=I_{T_0}-P_{0,q}\), the profiled criterion penalizes only the component of the residual orthogonal to \(\mathrm{Null}(K_q)\). Using the projection identity \(r'(I_{T_0}-P_{0,q})r=\min_{\gamma}\|r-Z_{0,q}\gamma\|_2^2\), where \(Z_{0,q}\) is a basis of \(\mathrm{Null}(K_q)\), the profiled objective at \(\rho=1\) is equivalent to the simultaneous least-squares fit of donor weights and unregularized null-space coefficients \(\gamma\). Thus, for \(q=1\), since \(\mathrm{Null}(K_1)=\mathrm{span}\{\mathbf 1_{T_0}\}\), the endpoint estimator is synthetic control with an intercept; for \(q=2\), since \(\mathrm{Null}(K_2)=\mathrm{span}\{\mathbf 1_{T_0},t_\mathrm{pre}\}\), the endpoint estimator is synthetic control with an intercept and a linear trend in time.

By Proposition (ref), null-space components of the residual are not tuned by $\rho$. For every $v \in \mathrm{Null}(K_q)$ and every $\rho \in (0,1)$, $(I_{T_0} + \lambda_\rho K_q)v = v$, so $S_{\rho,q}v = v$ and $W_{\rho,q}v = 0$; the same identities hold at $\rho \in \{0,1\}$ by definitions. Thus, a constant component of the residual (for $q=1$) or any intercept-plus-linear-trend component (for $q=2$) is always absorbed entirely into $\hat E_\mathrm{pre}(\rho,q)$ and contributes nothing to the donor-matching criterion. The tuning parameter $\rho$ reallocates only the residual components outside $\mathrm{Null}(K_q)$ between donor matching and the smooth component; null-space components are always assigned to the smooth branch.

With these properties at $\rho=0$ and $\rho=1$, we can extend the definition of the HSC weight estimation to $\rho \in [0,1]$.

definition[Harmonic Synthetic Control, \texorpdfstring{$\rho\in[0,1]$}{rho in [0,1]}] For \(\rho\in[0,1]\) and \(q\in\{1,2\}\), define the HSC weight estimator by \begin{equation} \hat\omega(\rho,q) \;\in\; \mathop{\rm argmin}_{\omega\in\Delta_{N_0}} \left\{ r_\mathrm{pre}(\omega)'W_{\rho,q}r_\mathrm{pre}(\omega) +\zeta^2T_0\|\omega\|_2^2 \right\}, \end{equation} where \(W_{\rho,q}\) is the continuously extended metric family from Proposition (ref). The associated treated-unit-specific smooth component is \begin{equation} \hat E_\mathrm{pre}(\rho,q):=S_{\rho,q}\,\hat r_\mathrm{pre}(\rho,q), \qquad \hat r_\mathrm{pre}(\rho,q):=Y_{\mathrm{pre}}-X_\mathrm{pre}\hat\omega(\rho,q), \end{equation} with \(S_{\rho,q}\) continuously extended to \([0,1]\) as in Proposition (ref).

This profiled Definition (ref) is the HSC weight estimator used throughout the remainder of the paper. Section (ref) discusses the interpretation of this metric from the spectral perspective and shows how $W_{\rho,q}$ amplifies or downweights components at different frequencies.

Forecast operator and the HSC counterfactual

The previous subsection develops the HSC weight estimator $\hat\omega(\rho,q)$ and shows that it arises from minimizing the pre-treatment residual norm under the metric $W_{\rho,q}$. However, the donor weights alone do not define the full counterfactual. The profiled representation produces a smooth component $\hat E_\mathrm{pre}(\rho,q)$, which captures the smooth part of the donor matching residual. To construct a post-treatment counterfactual, this smooth component is extrapolated to post-treatment periods. This term serves a similar function to the bias-correction term in augmented synthetic control ben2021augmented.

To see why a forecasting step is needed, consider the decomposition implied by HSC. For each $(\rho,q)$, the pre-treatment residual $\hat r_\mathrm{pre}(\rho,q):=Y_{\mathrm{pre}}-X_\mathrm{pre}\hat\omega(\rho,q)$ is split by the smoother $S_{\rho,q}$ into two components: \[ \hat r_\mathrm{pre}(\rho,q) = \underbrace{S_{\rho,q}\hat r_\mathrm{pre}(\rho,q)}_{\displaystyle\hat E_\mathrm{pre}(\rho,q)} \;+\; \underbrace{(I_{T_0}-S_{\rho,q})\hat r_\mathrm{pre}(\rho,q)}_{\displaystyle\hat u_\mathrm{pre}(\rho,q)}. \] The first term, $\hat E_\mathrm{pre}(\rho,q)$, is the treated-unit-specific smooth component: the smooth part of the residual that cannot be explained by the donors. The second term, $\hat u_\mathrm{pre}(\rho,q)$, is the rough remainder, containing the higher-roughness components of the residual. Note that any component of the residual in $\mathrm{Null}(K_q)$ (a constant for $q=1$, or an intercept-plus-linear-trend component for $q=2$) is contained in $\hat E_\mathrm{pre}(\rho,q)$.

The treated unit's pre-treatment outcome therefore admits the three-part decomposition.

equation[equation omitted — 264 chars of source]

In the post-treatment period, the donor-matched component extends to $X_\mathrm{post}\hat\omega(\rho,q)$ using observed donor outcomes. The rough remainder is the part of the pre-treatment residual that is suppressed by the smoother $S_{\rho,q}$; it is treated as noise and is not carried forward. The smooth component, however, represents a treated-unit-specific component that, if present in the pre-treatment period, is likely to persist. Discarding it would leave a predictable bias in the counterfactual. This motivates the introduction of a forecast operator to extrapolate $\hat E_\mathrm{pre}(\rho,q)$ into post-treatment periods.

Let $G_q\in\ensuremath{\mathbb{R}}^{T_\mathrm{post}\times T_0}$ be a deterministic linear forecast operator. The HSC counterfactual is defined as

equation[equation omitted — 145 chars of source]

The counterfactual thus combines two branches: a donor-matching component, $X_\mathrm{post}\hat\omega(\rho,q)$, that uses cross-sectional information to predict the treated unit's counterfactual, and a time series forecasting component, $G_q\hat E_\mathrm{pre}(\rho,q)$, that uses temporal information in the pre-treatment residual to extrapolate the smooth component. Writing $\Pi_{\rho,q}:=G_qS_{\rho,q}\in\ensuremath{\mathbb{R}}^{T_\mathrm{post}\times T_0}$ for the composed forecast-smoother, the counterfactual can also be expressed as

equation[equation omitted — 157 chars of source]

The forecast operator $G_q$ is treated as a fixed linear map throughout the analysis of donor weights and counterfactual construction. In practice, $G_q$ may be estimated from the pre-treatment data. We impose one structural requirement on $G_q$.

definition[Admissible forecast operator] A linear operator $G_q:\ensuremath{\mathbb{R}}^{T_0}{\rightarrow}\ensuremath{\mathbb{R}}^{T_\mathrm{post}}$ is an admissible forecast operator of order $q$ if, for every polynomial $p$ of degree less than $q$, \begin{equation} G_q\,\bigl(p(t)\bigr)_{t=1}^{T_0} \;=\; \bigl(p(t)\bigr)_{t=T_0+1}^{T_0+T_\mathrm{post}}. \end{equation}

For $q=1$, the requirement reduces to $G_1\,\mathbf{1}_{T_0} = \mathbf{1}_{T_\mathrm{post}}$: constants are continued as constants. For $q=2$, it adds $G_2\,(1,2,\ldots,T_0)' = (T_0+1,\ldots,T_0+T_\mathrm{post})'$: linear trends are continued as linear trends. These are the only cases we use in what follows.

Definition (ref) ensures that the forecast operator extrapolates the components which the roughness penalty leaves unpenalized, namely those in $\mathrm{Null}(K_q)$, to the post-treatment period without distortion. The components that pass mechanically into $\hat E$ (constants for $q=1$; constants and linear trends for $q=2$) are thereby continued unchanged. Since $S_{\rho,q}$ also leaves $\mathrm{Null}(K_q)$ untouched for every $\rho\in[0,1]$, the composed operator $\Pi_{\rho,q} = G_q S_{\rho,q}$ inherits the same continuation property whenever $G_q$ is admissible.

Definition (ref) is a design requirement on the forecast operator, not a property that generic forecasters automatically satisfy. The simplest admissible operators are constant extrapolation (for $q=1$) and linear extrapolation (for $q=2$). Constant extrapolation defines $(G_1^{\mathrm{const}}x)_h := x_{T_0}$ for $h=1,\dots,T_\mathrm{post}$: each post-treatment period receives the last pre-treatment value. Then $G_1^{\mathrm{const}}\mathbf{1}_{T_0}=\mathbf{1}_{T_\mathrm{post}}$. Linear extrapolation defines $(G_2^{\mathrm{lin}}x)_h := x_{T_0}+h\cdot(x_{T_0}-x_{T_0-1})$ for $h=1,\dots,T_\mathrm{post}$. Then $G_2^{\mathrm{lin}}\mathbf{1}_{T_0}=\mathbf{1}_{T_\mathrm{post}}$ and $G_2^{\mathrm{lin}}t_\mathrm{pre}=t_\mathrm{post}$.

In practice, one may wish to use a richer forecasting model, for example an autoregressive or ARIMA-type specification, to extrapolate $\hat E_\mathrm{pre}(\rho,q)$. Because such procedures involve parameter estimation, the resulting forecast operator $\hat G_q$ is data-dependent and need not be admissible on its own. A convenient way to enforce the requirement is to separate the null-space and non-null-space components of the input.

Recall that $P_{0,q}$ denotes the orthogonal projector onto $\mathrm{Null}(K_q)$, and $P_{\perp,q} := I_{T_0} - P_{0,q}$. Definition (ref) constrains a forecast operator only through its action on $\mathrm{Null}(K_q)$. Let $G_q^{\mathrm{null}}$ denote the unique admissible action on the null space: it has the property that, for $v\in\mathrm{Null}(K_q)$ equal to $\bigl(p(t)\bigr)_{t=1}^{T_0}$ for the unique polynomial $p$ of degree less than $q$, $G_q^{\mathrm{null}} v = \bigl(p(t)\bigr)_{t=T_0+1}^{T_0+T_\mathrm{post}}$. Construct the forecast operator \[ \widetilde G_q := G_q^{\mathrm{null}} P_{0,q} + \hat G_q P_{\perp,q}. \] For every $z \in \mathrm{Null}(K_q)$, $\widetilde G_q z = G_q^{\mathrm{null}} z$, so Definition (ref) holds by construction. For every $z \in \mathrm{Null}(K_q)^\perp$, $\widetilde G_q z = \hat G_q z$, so the data-driven forecaster retains full flexibility on the $\mathrm{Null}(K_q)^{\perp}$ component. For components in $\mathrm{Null}(K_q)$, $G_q^{\mathrm{null}}$ simply lets the constant ($q=1$) or the linear trend ($q=2$) persist into the post-treatment periods. These are the carry-forward $G_1^{\mathrm{const}}$ and the two-point linear extension $G_2^{\mathrm{lin}}$ defined above, respectively.

The tuning parameter $\rho$ determines the content of the smooth component $\hat E_\mathrm{pre}(\rho,q)$ and therefore what the forecast operator $G_q$ is asked to extrapolate. At $\rho=1$, the smoother $S_{1,q} = P_{0,q}$ extracts only the null-space component of the residual, and $G_q$ continues these unpenalized components into the post-treatment period. For interior values $\rho \in (0,1)$, the smooth component absorbs additional variation beyond the null space, and $G_q$ extrapolates these additional components. At $\rho=0$, the smoother $S_{0,q} = I_{T_0}$ assigns the entire residual to $\hat E_\mathrm{pre}$, so $G_q$ carries the full pre-treatment discrepancy forward. The choice of $\rho$ thus calibrates the burden on $G_q$: minimal at $\rho=1$, where extrapolation amounts to leaving null-space components unchanged, and maximal at $\rho=0$, where the entire residual must be forecast.

Spectral Interpretation and Tuning

The previous section defined the HSC estimator and derived its profiled representation. We now develop the complementary perspective on this allocation. The metric $W_{\rho,q}$ acts as a frequency-dependent gain function, so $\rho$ determines which frequency components of $Y_\mathrm{pre}$ enter the donor-matching criterion and which are diverted to the time series forecaster. We use cross-validation to select $\rho$ by minimizing out-of-sample prediction error. We then illustrate the resulting adaptation using the simulated data introduced in Section (ref).

Spectral decomposition of the HSC metric

The previous section showed that $\rho$ controls a smooth allocation of $Y_\mathrm{pre}$ between donor matching and the time series forecaster. We now examine the complementary question: how does $\rho$ shape the donor-matching criterion? The spectral decomposition of $W_{\rho,q}$ provides a precise answer and reveals that the key difference between HSC and existing synthetic control methods lies in how each method weights different frequency components in the optimization criterion.

Since the penalty matrix $K_q = D_q^{\top} D_q$ is symmetric positive semidefinite, $\mathbb{R}^{T_0}$ admits an orthonormal basis of eigenvectors of $K_q$, denoted $v_{1,q}, \ldots, v_{T_0,q}$, with corresponding eigenvalues $0 \leq \mu_{1,q} \leq \cdots \leq \mu_{T_0,q}$. Any vector of length $T_0$ can be written as a linear combination of these eigenvectors. Each eigenvalue $\mu_{j,q}$ measures the roughness of the corresponding basis function as assessed by the $q$th-difference penalty: $\mu_{j,q} = v_{j,q}' K_q v_{j,q} = \|D_q v_{j,q}\|_2^2$. The first $q$ eigenvectors span the null space $\mathrm{Null}(K_q)$, corresponding to constants when $q=1$ and constants and linear trends when $q=2$, with $\mu_{j,q}=0$. As the index $j$ increases beyond $q$, the eigenvectors oscillate progressively more rapidly and $\mu_{j,q}$ grows.\footnote{For $q=1$, the eigenvalues admit the closed form $\mu_{j,1} = 4\sin^2\!\bigl((j-1)\pi/(2T_0)\bigr)$ for $j=1,\dots,T_0$. For $q=2$, the eigenvalue spectrum grows more steeply, as illustrated in Figure (ref)(a).} Small $\mu_{j,q}$ corresponds to low-frequency variation, such as slow-moving trends and long cycles, while large $\mu_{j,q}$ corresponds to high-frequency variation, including rapid, short-run oscillations.

For a residual vector $r_{\mathrm{pre}}(\omega)$, define its spectral coordinates by $\tilde r_j(\omega) := v_{j,q}' r_{\mathrm{pre}}(\omega)$, the coefficient on the $j$-th eigenvector. One can thus write $r_{\mathrm{pre}}(\omega) = \sum_{j=1}^{T_0} \tilde r_j(\omega) v_{j,q}$. In these coordinates, the smoother $S_{\rho,q}$ and the profiled metric $W_{\rho,q}$ each act componentwise through scalar functions of the eigenvalue $\mu$:

equation[equation omitted — 203 chars of source]
equation[equation omitted — 225 chars of source]

The shrinkage function $s_q(\mu;\rho)$ is decreasing in $\mu$: smoother components (small $\mu$) in $r_{\mathrm{pre}}(\omega)$ are less shrunk toward zero than rougher components (larger $\mu$). Components in $\mathrm{Null}(K_q)$ ($\mu=0$) survive intact. The weight function $w_q(\mu;\rho)$ is increasing in $\mu$: rougher components in $r_{\mathrm{pre}}(\omega)$ receive more emphasis in the donor-matching problem. Components in $\mathrm{Null}(K_q)$ ($\mu=0$) are completely excluded from the donor-matching problem. From this spectral perspective, HSC therefore routes low-frequency discrepancy primarily to the time series branch and high-frequency discrepancy primarily to the donor-matching branch.

Figure (ref) displays these two functions and their consequences for the simulated pre-treatment series introduced in Section (ref), with $q=1$ on the left and $q=2$ on the right. We highlight three features of the spectral decomposition.

First, Panel (a) plots the shrinkage function $s_q(\mu;\rho)$ for five values of $\rho$. This function determines how much of each spectral coordinate of the pre-treatment residual is retained in the smooth component $\hat E_\mathrm{pre}(\rho,q)$. For $\rho \in (0,1)$, each curve equals one at $\mu=0$ (null-space components are fully retained) and decays smoothly toward zero as $\mu$ increases, with smaller $\rho$ producing slower decay. At $\rho=0$, the shrinkage function equals one everywhere: the entire residual is assigned to $\hat E_\mathrm{pre}$ and extrapolated by $G_q$. At $\rho=1$, the function collapses to the null-space indicator, $s_q(\mu;1) = \mathbf{1}\{\mu=0\}$. Intermediate $\rho$ values interpolate smoothly between these extremes.

Second, Panel (b) plots the complementary weight function $w_q(\mu;\rho)$, which determines how much of each spectral coordinate enters the donor-matching criterion (ref). Every curve passes through $(\mu,w) = (1,1)$: components with $\mu > 1$ are amplified relative to uniform weighting, components with $\mu < 1$ are shrunk. At $\rho=0$, $w_q(\mu;0) = \mu$ (the 45-degree line), so the profiled metric reduces to $K_q$ and donors are matched entirely in $q$th differences. At $\rho=1$, $w_q(\mu;1) = \mathbf{1}\{\mu>0\}$ (the horizontal line), so all components outside $\mathrm{Null}(K_q)$ are retained intact in the donor-matching criterion, while the null-space directions are absorbed by $\hat E_\mathrm{pre}$. The $\rho=1$ endpoint thus coincides with SC with intercept at $q=1$ and with SC with intercept-plus-linear-trend at $q=2$.

Third, the eigenvalue ranges of $K_1$ and $K_2$ differ substantially. For $q=1$, the eigenvalues of $K_1$ lie in $[0,4]$ when $T_0=80$, whereas for $q=2$ the eigenvalues of $K_2$ span a much wider range (up to roughly $16$). Under the same $\rho$, the wider spectrum therefore yields more extreme amplification of rough components in the sense of $K_2$ than in the sense of $K_1$ in the donor-matching criterion.

figure[figure omitted — 2,085 chars of source]

Because $W_{\rho,q}$ is symmetric positive semidefinite, it admits a symmetric square root $\widetilde C_{\rho,q} := W_{\rho,q}^{1/2} = V_q \mathop{\rm diag}\!\bigl(\sqrt{w_q(\mu_{j,q};\rho)}\bigr) V_q'$, where $V_q := (v_{1,q}, \ldots, v_{T_0,q})$ is the orthogonal matrix of eigenvectors introduced above. The donor-matching problem can equivalently be written as \[ \hat\omega(\rho,q) \in \arg\min_{\omega \in \Delta_{N_0}} \left\{ \|\widetilde C_{\rho,q}\,r_{\mathrm{pre}}(\omega)\|_2^2 + \zeta^2 T_0 \|\omega\|_2^2 \right\}. \] The transformed series $\widetilde C_{\rho,q}\,Y_\mathrm{pre}$ makes the effect of $\rho$ on the data visible. Panel (c) of Figure (ref) plots this transformation for the treated unit from the shared + idiosyncratic stochastic trend regime ($\kappa=2$) of Figure (ref), at five values of $\rho$.

In the $q=1$ panel, at $\rho=1$, the transformation projects out the constant direction and retains all other components with equal weight; the result resembles a demeaned version of $Y_\mathrm{pre}$. At $\rho=0$, the transformation reduces to $K_1^{1/2}$, which has the same quadratic form as the first-difference operator ($\|K_1^{1/2} r\|_2^2 = \|D_1 r\|_2^2$): slow trends are strongly attenuated and period-to-period oscillations are amplified. At $\rho \in \{0.5, 0.75, 0.95\}$, low-frequency components are progressively attenuated while high-frequency components are increasingly amplified. The $q=2$ panel exhibits the same pattern with sharper contrast. At $\rho=1$, the transformation projects out both the constant and linear-trend directions. At $\rho=0$, it applies $K_2^{1/2}$, which has the same quadratic form as the second-difference operator: the transformed series is visibly rougher than the raw data because rapid oscillations that contribute little to the original amplitude are greatly amplified.

The name “harmonic” reflects the spectral structure of HSC. The weight $w_q(\mu;\rho) = \mu/\bigl((1-\rho) + \rho\mu\bigr)$ admits the decomposition \[ \frac{1}{w_q(\mu;\rho)} \;=\; (1-\rho)\cdot\frac{1}{\mu} \;+\; \rho\cdot 1, \] which identifies $w_q(\mu;\rho)$ as the weighted harmonic mean of $\mu$ and $1$ with weights $(1-\rho)$ and $\rho$. The two endpoints recover the arguments themselves: $w_q(\mu;0) = \mu$ and $w_q(\mu;1) = \mathbf{1}\{\mu>0\}$. This structure mirrors the coefficients $1/\rho$ and $1/(1-\rho)$ on the smoothness and residual terms in the primal objective (ref). Because the harmonic mean is dominated by the smaller of its arguments, $w_q(\mu;\rho)$ is small whenever either $\mu$ is small (a low-frequency, near-null-space component) or $\rho$ is small (the user has shifted emphasis toward $\hat E_\mathrm{pre}$), giving the metric its characteristic soft cutoff.

Selection of $\rho$

The central tuning decision in HSC is how to allocate low-frequency components of the pre-treatment discrepancy between donor matching and the smooth component $\hat E_\mathrm{pre}$. We select $\rho$ by rolling-origin cross-validation using only pre-treatment data of the treated unit:\footnote{Each unit may carry its own idiosyncratic stochastic trend, so the $\rho$ that minimizes prediction error on a control unit need not be appropriate for the treated unit. We therefore restrict cross-validation to the treated unit's own pre-treatment history.} at each fold, HSC is fitted on a shortened training window and used to predict held-out pre-treatment outcomes, and the $\rho$ minimizing average prediction error is selected. This mimics the forecasting task the estimator faces at the treatment date.

The cross-validation procedure requires the researcher to specify four inputs: the smoothness order $q \in \{1,2\}$, the forecast operator $G_q$, a forecast horizon $h \ge 1$, and a number of folds $L \ge 1$. The rolling forecast origins are then $k_\ell = T_0 - h - L + \ell$ for $\ell = 1, \ldots, L$; by construction, $k_L + h = T_0$, so all validation windows lie within the pre-treatment period and no post-treatment data are used. At each origin $k_\ell$, the training window is $\{1,\dots,k_\ell\}$ and the validation window is $\{k_\ell+1,\dots,k_\ell+h\}$. For each candidate $\rho$ and each origin $k_\ell$, the procedure (i) re-estimates the HSC weights $\hat\omega(\rho,q)$ and the smooth component $\hat E_\mathrm{pre}(\rho,q)$ on the training window, (ii) refits any data-dependent parameters of $G_q$ on the training window (for instance, the coefficients of an ARIMA specification), (iii) constructs the counterfactual forecast $\hat Y_{k_\ell+s}(0;\rho,q)$ for the validation window by combining the donor-matched component with the $G_q$-extrapolated smooth component, and (iv) computes the squared prediction errors against the actual outcomes $Y_{k_\ell+1},\dots,Y_{k_\ell+h}$. Because all quantities are re-estimated within each fold, the validation error is an out-of-sample measure of the joint performance of donor matching and time series extrapolation. The cross-validation criterion is

equation[equation omitted — 155 chars of source]

where $\hat Y_{k_\ell+s}(0;\rho,q)$ denotes the HSC counterfactual prediction for period $k_\ell+s$ constructed from the training window $\{1,\dots,k_\ell\}$.

One practical consideration in constructing the candidate grid is worth noting. The equivalent smoothing parameter $\lambda_\rho = \rho/(1-\rho)$ is a convex function of $\rho$ that increases slowly near $\rho=0$ but diverges as $\rho \uparrow 1$. Consequently, a uniformly spaced grid in $\rho$ induces a nonuniform grid in $\lambda_\rho$, with increasingly coarse coverage on the $\lambda_\rho$ scale near $\rho=1$. When the cross-validation curve exhibits dramatic change near $\rho=1$, a convenient alternative is to construct the grid on a logarithmic scale in $\lambda_\rho$ and map back via $\rho = \lambda_\rho/(1+\lambda_\rho)$, while retaining the boundary $\rho=0$ and $\rho=1$ explicitly.

Beyond the choice of $\rho$, the same procedure can also be used to select $(\rho, q)$ jointly by minimizing $\mathrm{CV}(\rho, q)$ over the pair. If multiple forecast operators are under consideration, the grid extends further to $(\rho, q, G_q)$ triples, comparing, for example, constant extrapolation against an ARIMA specification for each $(\rho, q)$ combination. Joint selection remains computationally inexpensive because the quadratic program for each $(\rho, q, k_\ell)$ combination is fast to solve.\footnote{Our implementation uses Gurobi.}

An illustrative example

We illustrate HSC on the two simulated regimes from Section (ref). Two configurations vary both the smoothness order $q$ and the forecast operator $G_q$ to show the roles of all three tuning choices ($\rho$, $q$, and $G_q$). In both configurations, $\rho$ is selected by the rolling-origin cross-validation procedure of Section (ref) with one-step-ahead horizon ($h=1$), 10 folds ($L=10$), and a grid of 21 equally spaced values in $[0,1]$.

\paragraph*{Configuration 1: $q=1$, constant extrapolation.} The first configuration sets $q=1$ and uses the constant-extrapolation operator $G_1^{\mathrm{const}}$, the most conservative admissible forecaster: each post-treatment period receives the last pre-treatment value of $\hat E_\mathrm{pre}$.

Figure (ref) displays the results. Each row corresponds to one regime: shared stochastic trend ($\kappa=0$, top) and shared + idiosyncratic stochastic trend ($\kappa=2$, bottom). The three columns show the cross-validation curve, the counterfactual fit, and the decomposition of the fitted counterfactual into its donor-matched and time series components.

In the regime with only shared stochastic trend, cross-validation selects $\hat\rho=1$, the endpoint corresponding to SC with intercept. This is consistent with the spectral interpretation of Section (ref): when the stochastic trend factor is shared across all units, the donor pool can reproduce the treated unit's low-frequency dynamics, so HSC lets the donor-matching branch carry most of the load in modeling the treated series. The decomposition panel confirms this allocation: the donor component $X_\mathrm{post}\hat\omega$ (green) almost coincides with the counterfactual $\hat Y_\mathrm{post}(0)$ (blue), while the extrapolated smooth component $G_1^{\mathrm{const}}\hat E_\mathrm{pre}$ (red) is a flat horizontal line.

In the regime with both shared and idiosyncratic stochastic trends, cross-validation selects $\hat\rho=0$. The donor pool can no longer reproduce the treated unit's idiosyncratic stochastic trend, so cross-validation finds that matching entirely in first differences, with the time series branch carrying more of the prediction, yields the best out-of-sample predictions. The decomposition panel reflects this: the extrapolated smooth component $G_1^{\mathrm{const}}\hat E_\mathrm{pre}$ now carries a large level correction that accounts for the drift accumulated up to $T_0$. In this regime, HSC reduces to ridge-regularized synthetic control in first differences, anchored at $T_0$, as shown in Section (ref).

figure[figure omitted — 1,638 chars of source]

The contrast between $\hat\rho=1$ and $\hat\rho=0$ provides a sharp demonstration of the adaptive mechanism of HSC: the same estimator with the same tuning procedure selects the two extreme endpoints of the $\rho$ continuum, recovering pure level matching when the stochastic trend is shared and pure differencing when the drift is partly idiosyncratic, precisely the regime-dependent behavior that Section (ref) argues is needed.

\paragraph*{Configuration 2: $q=2$, ARIMA forecast.}

The second configuration sets $q=2$ and uses an ARIMA(1,1,0) model as the data-driven forecaster $\hat G_q$ within the null-space separation construction of Section (ref). Definition (ref) is satisfied by design: the null-space component of $\hat E_\mathrm{pre}$ (which now includes both an intercept and a linear trend) is extrapolated by the canonical linear forecaster $G_2^{\mathrm{lin}}$, while the non-null-space component is extrapolated by an ARIMA model.

Figure (ref) displays the results in the same format as Figure (ref). In the shared stochastic trend regime, cross-validation again selects $\hat\rho=1$. As in Configuration 1, the donor pool suffices to match the treated unit's dynamics, and the time series forecaster plays a minimal role. The post-treatment RMSE (2.7) is comparable to Configuration 1 (2.6).

In the regime with both shared and idiosyncratic stochastic trends, cross-validation selects $\hat\rho=0.33$. The decomposition panel reveals that the smooth component $\hat E_\mathrm{pre}$ exhibits a clear downward trend in the pre-treatment period, and the ARIMA forecaster extrapolates this trend into the post-treatment window. This trending extrapolation is visible in the red dotted line, which continues to decline after $T_0$ rather than remaining flat as in Figure (ref). As a result, less of the idiosyncratic drift needs to be absorbed by the level shift alone, and the post-treatment RMSE improves from 4.3 to 3.4.\footnote{The small RMSE difference (4.3 vs.\ 4.4) between Configuration 1 at $\hat\rho=0$ and the unregularized first-differenced SC of Figure (ref) reflects the ridge term $\zeta^2 T_0 \|\omega\|_2^2$.}

figure[figure omitted — 1,220 chars of source]

This comparison highlights the complementary roles of the three tuning choices. The tuning parameter $\rho$ controls how much of the pre-treatment discrepancy is allocated to each branch. The forecast operator $G_q$ determines how the smooth component is extrapolated. The smoothness order $q$ determines what counts as “smooth”. Across both configurations, the cross-validation procedure adapts $\rho$ to the data-generating regime.

Prediction-Error Decomposition

Having established the HSC estimator and its spectral interpretation, we now develop a formal decomposition of the counterfactual prediction error. This decomposition separates the prediction error into a weight-estimation component, which captures the discrepancy introduced by estimating donor weights from observed outcomes rather than from the underlying shared structure, and a forecasting component, which captures the prediction error that would remain even with oracle weights. Each component depends on $\rho$, and their interplay provides a formal decomposition that helps interpret the adaptive $\rho$-selection illustrated in Section (ref).

Throughout this section we fix $q\in\{1,2\}$ and suppress the $q$-subscript when no ambiguity arises, writing $K=K_q$, $S_\rho=S_{\rho,q}$, $W_\rho=W_{\rho,q}$, $P_0=P_{0,q}$, and $P_\perp=I_{T_0}-P_{0,q}$. We use the null-space-separated forecast operator $\widetilde G_q=G_q^{\mathrm{null}}P_0+\hat G_q P_\perp$ constructed in Section (ref) and write $\Pi_\rho:=\widetilde G_q S_\rho$ for the composed operator. Since $S_\rho$ leaves $\mathrm{Null}(K)$ components intact and $\widetilde G_q$ continues $\mathrm{Null}(K)$ components to post-treatment periods (Definition (ref)), the composed operator inherits the null-space continuation property:

equation[equation omitted — 114 chars of source]

Outcome model and oracle benchmark

We connect the prediction error to the data-generating structure discussed in Section (ref). Recall from Section (ref) the decomposition $Y_{it}(0)=L_{it}+R_{it}+\varepsilon_{it}$, where $L_{it}$ denotes the shared component that provides the signal for donor matching, $R_{it}$ is an idiosyncratic stochastic trend, and $\varepsilon_{it}$ is idiosyncratic short-run noise. For the prediction-error analysis, we work with the two-component grouping

equation[equation omitted — 83 chars of source]

where $\mathcal{R}_{it}:=R_{it}+\varepsilon_{it}$ collects all components not explained by the shared component.\footnote{The grouping $\mathcal{R}_{it}=R_{it}+\varepsilon_{it}$ is adopted for analytical convenience: the prediction-error decomposition distinguishes between the shared component $L$ and everything else, regardless of whether that remainder is short-run noise or a stochastic trend. The persistence structure of $\mathcal{R}_{it}$ matters for the spurious matching channel of Term A (Section (ref)).}

Define donor and treated stacks for $L$ and $\mathcal{R}$ analogously to $X_\mathrm{pre}$ and $Y_\mathrm{pre}$: \[ X_\mathrm{pre} = L_{0,\mathrm{pre}}+\mathcal{R}_{0,\mathrm{pre}}, \quad Y_\mathrm{pre} = L_{1,\mathrm{pre}}+\mathcal{R}_{1,\mathrm{pre}}, \qquad X_\mathrm{post} = L_{0,\mathrm{post}}+\mathcal{R}_{0,\mathrm{post}}, \quad Y_\mathrm{post}(0) = L_{1,\mathrm{post}}+\mathcal{R}_{1,\mathrm{post}}. \]

We compare the HSC estimator against an oracle that observes the signal $L$ directly and solves the $\rho=1$ HSC program on $L$ (with the null space projected out):\footnote{The oracle uses the same ridge parameter $\zeta$ as the HSC estimator. This ensures that the oracle benchmark is not trivially superior due to different regularization, so the decomposition isolates the effects of observing $Y$ rather than $L$ and of using $W_\rho$ rather than $P_\perp$.}

equation[equation omitted — 227 chars of source]

where $P_\perp$ projects onto $\mathrm{Null}(K)^\perp$, removing an intercept ($q=1$) or an intercept and linear trend ($q=2$). This is the same quadratic program as the $\rho=1$ endpoint of HSC (Definition (ref) with $W_{1,q}=P_\perp$), but on the shared component $L$ rather than the observed outcome $Y=L+\mathcal{R}$. The oracle is treated as a fixed benchmark throughout the analysis.

The prediction-error decomposition

Recall from (ref) that the HSC counterfactual at any $\rho$ takes the form $\hat Y_\mathrm{post}(0;\rho,q) = X_\mathrm{post}\hat\omega(\rho,q) + \Pi_\rho\,r_\mathrm{pre}(\hat\omega(\rho,q))$, where $r_\mathrm{pre}(\omega):=Y_\mathrm{pre}-X_\mathrm{pre}\omega$ is the pre-treatment residual. We decompose the prediction error by introducing an oracle predictor that replaces the estimated weights with $\omega^\mathrm{oracle}$ but retains the same $\rho$-dependent smoother and forecaster:

equation[equation omitted — 182 chars of source]

where $r_\mathrm{pre}^\mathrm{oracle}:=Y_\mathrm{pre}-X_\mathrm{pre}\omega^\mathrm{oracle}$ is the oracle pre-treatment residual. Adding and subtracting $\widetilde Y_\mathrm{post}^\mathrm{oracle}(0;\rho)$ in the prediction error yields:

proposition[Prediction-error decomposition] For every $\rho\in[0,1]$, \begin{equation} \underbrace{Y_\mathrm{post}(0)-\hat Y_\mathrm{post}(0;\rho,q)}_{prediction error} \;=\; \underbrace{\bigl(\widetilde Y_\mathrm{post}^\mathrm{oracle}(0;\rho) -\hat Y_\mathrm{post}(0;\rho,q)\bigr)}_{\mathrm{Term A}(\rho) \ (weight estimation)} \;+\; \underbrace{\bigl(Y_\mathrm{post}(0) -\widetilde Y_\mathrm{post}^\mathrm{oracle}(0;\rho)\bigr)}_{\mathrm{Term B}(\rho) \ (forecasting)}. \end{equation} The two terms admit the closed forms \begin{align} \mathrm{Term A}(\rho) &= \bigl(X_\mathrm{post}-\Pi_\rho X_\mathrm{pre}\bigr) \bigl(\omega^\mathrm{oracle}-\hat\omega(\rho,q)\bigr), \\ \mathrm{Term B}(\rho) &= \bigl(Y_\mathrm{post}(0)-X_\mathrm{post}\omega^\mathrm{oracle}\bigr) \;-\; \Pi_\rho\bigl(Y_\mathrm{pre}-X_\mathrm{pre}\omega^\mathrm{oracle}\bigr). \end{align}

The derivation given in Appendix (ref) uses only the linearity of $r_\mathrm{pre}(\omega)$ in $\omega$ and the linearity of $\Pi_\rho$.

Term A isolates the cost of using the estimated weights $\hat\omega(\rho,q)$ rather than the oracle weights $\omega^\mathrm{oracle}$ in post-treatment periods. The weight discrepancy $\omega^\mathrm{oracle}-\hat\omega(\rho,q)$ is transferred into prediction error through the donor-forecast residual matrix $C_\rho:=X_\mathrm{post}-\Pi_\rho X_\mathrm{pre}$, whose column $j$ is the forecast residual when the composed operator $\Pi_\rho$ extrapolates donor $j$ from the pre-treatment to the post-treatment period. The subtracted term $\Pi_\rho X_\mathrm{pre}$ appears because the donor weights enter the counterfactual twice, directly through the post-treatment donor block $X_\mathrm{post}\omega$ and indirectly through the pre-treatment residual $r_\mathrm{pre}(\omega)=Y_\mathrm{pre}-X_\mathrm{pre}\omega$ that $\Pi_\rho$ extrapolates forward. Term A therefore depends on $\rho$ through both the estimated weights and $C_\rho$.

Term B is the prediction error that would remain even if the oracle weights were available. It measures how accurately the composed operator $\Pi_\rho$ extrapolates, from the pre-treatment to the post-treatment period, the discrepancy between the treated unit and its oracle synthetic counterfactual, where the oracle weights are those defined by the shared component $L$. Term B depends on $\rho$ through the smoother $S_\rho$ embedded in $\Pi_\rho$, which controls how much of the pre-treatment oracle residual is passed to the forecaster and how much is discarded.

Weight-estimation error: three channels of distortion

We now characterize the Term A cost by identifying three distinct channels through which the weight discrepancy $\omega^\mathrm{oracle}-\hat\omega(\rho,q)$ arises.

The estimated weights $\hat\omega(\rho,q)$ and the oracle weights $\omega^\mathrm{oracle}$ minimize different objectives over the same constraint set $\Delta_{N_0}$. The profiled HSC objective (Definition (ref)) is

equation[equation omitted — 140 chars of source]

while the oracle objective (ref) is

equation[equation omitted — 168 chars of source]

These two objectives share the same ridge penalty but differ in two ways. First, the metric: HSC evaluates the pre-treatment residual under $W_\rho$, whereas the oracle uses $P_\perp$. The two coincide only at $\rho=1$ (where $W_\rho=P_\perp$). Second, the data: HSC observes $X_\mathrm{pre}=L_{0,\mathrm{pre}}+\mathcal{R}_{0,\mathrm{pre}}$ and $Y_\mathrm{pre}=L_{1,\mathrm{pre}}+\mathcal{R}_{1,\mathrm{pre}}$, whereas the oracle operates on the signal $L$ alone. The weight discrepancy $\omega^\mathrm{oracle}-\hat\omega(\rho,q)$ reflects the combined effect of these two differences.

To separate the contributions of $L$ and $\mathcal{R}$ to the weight discrepancy, we decompose the oracle residual by component. For each $Z\in\{L,\mathcal{R}\}$, define $e^Z_\mathrm{pre}:=Z_{1,\mathrm{pre}}-Z_{0,\mathrm{pre}}\omega^\mathrm{oracle}\in\ensuremath{\mathbb{R}}^{T_0}$, so that $r_\mathrm{pre}^\mathrm{oracle}=e^L_\mathrm{pre}+e^\mathcal{R}_\mathrm{pre}$ by additivity (ref). The signal residual $e^L_\mathrm{pre}$ contains a null-space component $P_0\,e^L_\mathrm{pre}\in\mathrm{Null}(K)$ that the oracle objective does not use, and a complement $P_\perp\,e^L_\mathrm{pre}\in\mathrm{Null}(K)^\perp$ that the oracle directly minimizes. Since both $W_\rho$ and $P_\perp$ annihilate $\mathrm{Null}(K)$ components, we have $(W_\rho-P_\perp)\,P_\perp\,e^L_\mathrm{pre}=(W_\rho-P_\perp)\,e^L_\mathrm{pre}$ and $W_\rho\,P_\perp\,e^L_\mathrm{pre}=W_\rho\,e^L_\mathrm{pre}$.

The following proposition bounds $\|\mathrm{Term~A}(\rho)\|_2$ by three interpretable components that correspond to the metric difference, the data difference, and their interaction. Define the bound components

align[align omitted — 389 chars of source]

where $\|v\|_{Q_\rho^{-1}}:=\sqrt{v'Q_\rho^{-1}v}$,\; $\|v\|_{W_\rho}:=\sqrt{v'W_\rho v}$, and

equation[equation omitted — 116 chars of source]

is the Hessian of $F_\rho$ (up to the constant $2T_0$).

proposition[Weight-estimation error envelope] For every $\rho\in[0,1]$, \begin{equation} \bigl\|\mathrm{Term A}(\rho)\bigr\|_2 \;\le\; \mathcal{P}_\rho \bigl[ A_1(\rho)+A_2(\rho)+A_3(\rho) \bigr], \end{equation} where $\mathcal{P}_\rho :=\|(X_\mathrm{post}-\Pi_\rho X_\mathrm{pre})Q_\rho^{-1/2}\|_{\mathrm{op}}$ is the transfer multiplier.

Proposition (ref) bounds the weight-estimation error with four ingredients, each with a distinct role. The terms $A_1(\rho)$, $A_2(\rho)$, and $A_3(\rho)$ are three sources of the weight gap $\omega^\mathrm{oracle}-\hat\omega(\rho,q)$, each tied to a specific distortion mechanism. The matrix $Q_\rho$, the ridge-regularized Gram matrix of the filtered donor series $W_\rho^{1/2}X_\mathrm{pre}$, summarizes how sharply the filtered donor pool distinguishes different weight-reallocation patterns; the dual norm $\|\cdot\|_{Q_\rho^{-1}}$ through which $A_1$ and $A_2$ are measured reflects this geometry. The multiplier $\mathcal{P}_\rho$ converts the weight gap into a prediction error and is the only place in the Term A envelope where the forecast operator $\widetilde G_q$ enters. We unpack each ingredient in turn, then synthesize them into two opposing forces that shape Term A in $\rho$.

The metric distortion $A_1(\rho)$ measures the contribution to the weight gap $\omega^\mathrm{oracle}-\hat\omega$ that arises because HSC and the oracle use different metrics to evaluate the pre-treatment residual: HSC uses $W_\rho$, whereas the oracle uses $P_\perp$. Thus $A_1$ is the structural price HSC pays for evaluating the pre-treatment residual under $W_\rho$ rather than the oracle's $P_\perp$, which leaves the weight criterion no longer aligned with the oracle's identification target. Two limiting cases eliminate this price entirely. First, at $\rho=1$, $W_\rho=P_\perp$ and $A_1$ vanishes because the two metrics coincide. Second, $A_1$ vanishes for all $\rho$ whenever $P_\perp\,e^L_\mathrm{pre}=0$, that is, whenever the treated unit's signal $L_{1,\mathrm{pre}}$ can be written as a convex combination of donor signals $L_{0,\mathrm{pre}}\omega^\mathrm{oracle}$ plus a $\mathrm{Null}(K)$ component.\footnote{Concretely, $P_\perp\,e^L_\mathrm{pre}=0$ means that the oracle achieves zero residual-fit term in (ref): the mismatch in $L$ lies entirely in $\mathrm{Null}(K)$ (a level shift for $q=1$, an affine trend for $q=2$). Both $W_\rho$ and $P_\perp$ assign zero weight to $\mathrm{Null}(K)$, so the choice of metric is inconsequential.} Beyond these two cases, the metric discrepancy $W_\rho-P_\perp$ widens as $\rho$ decreases from $1$, which pushes $A_1$ upward. However, because $A_1$ is measured in the $Q_\rho^{-1}$ norm, which also depends on $\rho$, the net behavior of $A_1(\rho)$ need not be monotonic. Notably, $A_1$ can be nonzero even when $\mathcal{R}=0$, which is the sense in which it is a purely structural channel.

The interaction $A_2(\rho)$ captures how idiosyncratic components in the donor units $\mathcal{R}_{0,\mathrm{pre}}$ can perturb weight estimation when the signal $L$ is not perfectly matched. The oracle benchmark is defined using the shared component $L$ only, and the HSC objective is evaluated on the observed donor outcomes $X_{\mathrm{pre}}=L_{0,\mathrm{pre}}+\mathcal{R}_{0,\mathrm{pre}}$. As a result, $\mathcal{R}_{0,\mathrm{pre}}$ can accidentally align with $W_\rho\,e^L_{\mathrm{pre}}$ and appear to help reduce the remaining signal mismatch in sample, thereby shifting $\hat\omega(\rho)$ away from $\omega^\mathrm{oracle}$. This term is large when $\mathcal{R}_{0,\mathrm{pre}}$ has a substantial projection onto $W_\rho\,e^L_{\mathrm{pre}}$, and it vanishes whenever $P_\perp\,e^L_\mathrm{pre}=0$.

The spurious matching $A_3(\rho)$ measures how tempted the HSC weight criterion is to chase idiosyncratic components in the donors. Unlike $A_1$ and $A_2$, this channel does not require any mismatch in $L$. The oracle weights $\omega^\mathrm{oracle}$ are chosen to fit $L$ alone. The $\mathcal{R}$-mismatch $e^\mathcal{R}_\mathrm{pre}=\mathcal{R}_{1,\mathrm{pre}}-\mathcal{R}_{0,\mathrm{pre}}\omega^\mathrm{oracle}$ is untouched by the oracle and can be large. The HSC weight criterion thus has an incentive to deviate from $\omega^\mathrm{oracle}$ to absorb it, pulling weights toward donors whose idiosyncratic components happen to co-move with $\mathcal{R}_{1,\mathrm{pre}}$.

The sensitivity of $A_3$ to $\rho$ depends critically on whether $\mathcal{R}$ is short-run noise or a stochastic trend. At $\rho=1$, $W_\rho=P_\perp$, and the weight criterion applies no spectral down-weighting beyond removing the null-space component. If $\mathcal{R}$ contains a random-walk component, $\|e^\mathcal{R}_\mathrm{pre}\|_{P_\perp}$ grows at rate $O_p(T_0)$, so $A_3(1)=O_p\!\bigl(T_0^{1/2}\bigr){\rightarrow}\infty$, which reflects the spurious regression phenomenon discussed in Section (ref). At $\rho=0$, the metric reduces to $W_0=K=D_q'D_q$, so the weight criterion operates on the $q$th differences of $e^\mathcal{R}_\mathrm{pre}$; differencing renders the random-walk component stationary and yields $A_3(0)=O_p(1)$, thereby controlling the spurious channel. For interior values $\rho\in(0,1)$, the spectral weight $w(\mu;\rho)=\mu/\bigl((1-\rho)+\rho\mu\bigr)$ interpolates smoothly between the two extremes: larger $\rho$ retains more low-frequency energy in the weight criterion and is therefore more vulnerable to spurious matching when $\mathcal{R}$ contains a stochastic trend, whereas smaller $\rho$ down-weights low-frequency components more aggressively, at the cost of the metric distortion already discussed for $A_1(\rho)$. When $e^\mathcal{R}_\mathrm{pre}$ contains only short-run noise, either because $\mathcal{R}$ itself is stationary or because the treated unit and the control donors form a cointegrated system under the oracle weights so that the stochastic-trend components cancel, $A_3$ is $O_p(1)$ at both endpoints.

The operative quantity is $\lambda_{\max}(Q_\rho^{-1})$, the largest eigenvalue of $Q_\rho^{-1}$. It measures how weakly the HSC criterion identifies the donor weights: a large value means the filtered donors are nearly collinear in some direction, so small score discrepancies are amplified through the dual norm $\|\cdot\|_{Q_\rho^{-1}}$, inflating $A_1$ and $A_2$. How $\lambda_{\max}(Q_\rho^{-1})$ changes with $\rho$ can be non-monotonic. At $\rho=0$, the metric is $K_q$, so low-frequency directions are already most strongly down-weighted while high-frequency directions are amplified. Moving $\rho$ away from zero gradually restores weight on low-frequency components and reduces the amplification of high-frequency components. Depending on the frequency composition of the control donors, these two effects can make $\lambda_{\max}(Q_\rho^{-1})$ peak at an interior $\rho\in(0,1)$. The ridge floor $\zeta^2$ in $Q_\rho$ guarantees that $\lambda_{\max}(Q_\rho^{-1})\le 1/\zeta^2$ even when the filtered donors are nearly collinear; the no-ridge comparison in Appendix (ref) shows that without this floor $\lambda_{\max}(Q_\rho^{-1})$ can grow dramatically.

The transfer multiplier $\mathcal{P}_\rho=\|(X_\mathrm{post}-\Pi_\rho X_\mathrm{pre})Q_\rho^{-1/2}\|_\mathrm{op}$ depends on $\rho$ through two distinct mechanisms. The first is the inverse curvature $Q_\rho^{-1/2}$, which reflects the same identification geometry that shapes $A_1$ and $A_2$. The second is $C_\rho$; its dependence on $\rho$ runs through the smoother $S_\rho$, which determines how much of each donor's pre-treatment path is passed to the forecaster.

Taken together, the envelope $\mathcal{P}_\rho[A_1+A_2+A_3]$ is governed by two opposing forces in $\rho$. At large $\rho$ the dominant cost is spurious matching: the $A_3$ channel grows when $\mathcal{R}$ carries stochastic trends, formalizing the spurious donor matching risk of Section (ref). At small $\rho$ the dominant cost is identification loss: the metric gap $W_\rho-P_\perp$ widens, possibly inflating $A_1$, and the filtered donor design $W_\rho^{1/2}X_\mathrm{pre}$ sheds low-frequency variation, inflating $\lambda_{\max}(Q_\rho^{-1})$. This is the over-filtering cost of Section (ref), and it is most severe when the donor series are dominated by low-frequency variation, as is typical for macroeconomic data. The net shape of the envelope in $\rho$ is therefore non-monotonic in general, and depends on the size and structure of the signal mismatch $e^L_\mathrm{pre}$, the persistence of $\mathcal{R}$, the frequency composition of the donor design $X_\mathrm{pre}$, and the forecast operator $\widetilde G_q$.

Forecasting error

Term B is the prediction error that would remain even if the HSC weights coincided with the oracle weights. Its size depends on how accurately the composed operator $\Pi_\rho=\widetilde G_q S_\rho$ extrapolates the oracle pre-treatment residual into the post-treatment window. Because $S_\rho$ varies with $\rho$ while $\widetilde G_q$ does not, the $\rho$-dependence of Term B is governed entirely by the smoother and by what it forwards to the forecaster.

Recall the oracle pre-treatment residual $r_\mathrm{pre}^\mathrm{oracle}$ from Section (ref), and define its post-treatment counterpart:

equation[equation omitted — 320 chars of source]

The null-space component is the exception under $\Pi_\rho$. By the null-space continuation property (ref), $\Pi_\rho v = G_q^{\mathrm{null}}v$ for every $v\in\mathrm{Null}(K)$ regardless of $\rho$: the smoother leaves such a component intact and the forecaster then continues it by $G_q^{\mathrm{null}}$, in both steps independently of $\rho$. The null-space content of the oracle residual therefore contributes a fixed offset to the post-treatment prediction at every $\rho$ and plays no role in the $\rho$-dependent tradeoff. Accordingly, define

equation[equation omitted — 336 chars of source]

Both vectors are the oracle residual with the same null-space content $P_0\,r_\mathrm{pre}^\mathrm{oracle}$ removed: directly in the pre-period, and through its canonical continuation $G_q^{\mathrm{null}}\,P_0\, r_\mathrm{pre}^\mathrm{oracle}$ in the post-period. With this common adjustment, $\mathrm{Term~B}(\rho)$ takes the form

equation[equation omitted — 137 chars of source]

with the derivation given in Appendix (ref). Both $\eta_\mathrm{post}^\mathrm{oracle}$ and $\eta_\mathrm{pre}^\mathrm{oracle}$ are $\rho$-independent; all $\rho$-dependence in $\mathrm{Term~B}(\rho)$ enters through the smoother $S_\rho$.

At the endpoint $\rho=0$, $S_0=I$ and $\mathrm{Term~B}(0)=\eta_\mathrm{post}^\mathrm{oracle}-\widetilde G_q\, \eta_\mathrm{pre}^\mathrm{oracle}$. At the other endpoint $\rho=1$, $S_1=P_0$, and since $\eta_\mathrm{pre}^\mathrm{oracle}\in\mathrm{Null}(K)^\perp$ by construction, $S_1\,\eta_\mathrm{pre}^\mathrm{oracle}=0$ and therefore $\mathrm{Term~B}(1)=\eta_\mathrm{post}^\mathrm{oracle}$. At $\rho=1$ the forecaster makes no contribution to Term B beyond the canonical null-space continuation already absorbed into $\eta_\mathrm{post}^\mathrm{oracle}$, so the choice of $\widetilde G_q$ is irrelevant at this endpoint. Between the endpoints, as $\rho$ increases from $0$ to $1$, the smoothed input $S_\rho\,\eta_\mathrm{pre}^\mathrm{oracle}$ decreases monotonically toward zero in the spectral sense of Section (ref), so $\rho$ is the dial that controls how much of $\eta_\mathrm{pre}^\mathrm{oracle}$ reaches the forecaster.

How $\rho$ affects $\|\mathrm{Term~B}\|$ then depends on how well the raw forecast $\widetilde G_q\,\eta_\mathrm{pre}^\mathrm{oracle}$ tracks $\eta_\mathrm{post}^\mathrm{oracle}$. When $\widetilde G_q\,\eta_\mathrm{pre}^\mathrm{oracle}$ already predicts $\eta_\mathrm{post}^\mathrm{oracle}$ well, the raw forecast is useful and $\rho=0$ is preferred. When $\widetilde G_q\,\eta_\mathrm{pre}^\mathrm{oracle}$ over-extrapolates the noisier part of the residual, smoothing its input first improves the forecast and an interior $\rho>0$ is preferred. When $\widetilde G_q\,\eta_\mathrm{pre}^\mathrm{oracle}$ is far from $\eta_\mathrm{post}^\mathrm{oracle}$, the forecaster is harmful, and $\rho=1$, which discards its input entirely and leaves $\mathrm{Term~B}$ equal to $\eta_\mathrm{post}^\mathrm{oracle}$, is preferred. This logic suggests that a longer post-treatment window favors a larger $\rho$: the time series forecaster becomes less reliable at distant horizons, which pushes the preferred regime toward stronger regularization or full suppression of the time series forecaster.

Implications for tuning

Sections (ref) and (ref) characterized the two errors that HSC trades off as $\rho$ varies. We now collect their implications for the choices a practitioner makes. HSC exposes four such choices, and the decomposition shows they play distinct roles. The tuning parameter $\rho$ allocates the pre-treatment residual between donor matching and the time series forecaster. The forecast operator $\widetilde G_q$ and the smoothness order $q$ are structural: $\widetilde G_q$ fixes how the non-null residual is extrapolated, and $q$ fixes what counts as smooth and how the null space is continued. The cross-validation horizon $h$ does not change the estimator; it determines which of the two errors the cross-validation criterion weighs most heavily. We take $\rho$ first, then $\widetilde G_q$, $q$, and $h$.

$\rho$ is the allocation lever, and the tradeoff it controls is two-sided. At large $\rho$ the weight criterion retains the low-frequency content of the pre-treatment residual, so when $\mathcal{R}$ carries a stochastic trend the spurious matching channel $A_3$ grows and Term A rises. At small $\rho$ the criterion is restricted toward high-frequency content: the metric distortion $A_1$ widens and the filtered donor design sheds the low-frequency variation that identifies the weights, so Term A rises again through identification loss. The $\rho$-shape of Term B is governed instead by how accurately the composed operator extrapolates the oracle residual: when the raw forecast tracks $\eta_\mathrm{post}^\mathrm{oracle}$ well a small $\rho$ is preferred, and when it does not a large $\rho$, which suppresses the forecaster's input, is preferred. Neither error is monotone in $\rho$, and their sum has no general optimum; the best $\rho$ depends on the data-generating regime.

The forecast operator $\widetilde G_q$ is the lever common to both errors. It enters Term A only through $C_\rho$ inside the transfer multiplier $\mathcal{P}_\rho$, and it drives Term B directly through $\widetilde G_q\,S_\rho\,\eta_\mathrm{pre}^\mathrm{oracle}$. A forecaster that extrapolates the oracle residual well lowers the Term B floor. Its effect on the multiplier $\mathcal{P}_\rho$ is separate: $\mathcal{P}_\rho$ depends on the forecaster only through the donor forecast residuals $C_\rho=X_\mathrm{post}-\Pi_\rho X_\mathrm{pre}$, and the same operator can behave differently on $\eta_\mathrm{pre}^\mathrm{oracle}$ than on the donor paths, so the two effects need not move together. In the Monte Carlo study, the constant carry-forward and the ARIMA$(1,1,0)$ forecasters both perform well.\footnote{In both cases the rule is the data-driven component $\hat G_q$ of the construction in Definition (ref): it is applied to the non-null part $P_{\perp,q}r$ of the pre-treatment residual $r$, while the null-space part $P_{0,q}r$ is continued by the canonical $G_q^{\mathrm{null}}$. When $\hat G_q$ is the constant carry-forward, the two parts recombine in closed form. For $q=1$, the continued mean plus the carried-forward demeaned residual equals the last entry $r_{T_0}$ held constant, so $\widetilde G_1$ coincides with the constant carry-forward applied directly to the raw pre-treatment residual. For $q=2$, the composed forecast at horizon $h$ is $r_{T_0}+\hat\beta\,h$, where $\hat\beta$ is the slope of the line fitted to the pre-treatment residual over $\mathrm{Null}(K_2)$; equivalently, the fitted linear trend is extrapolated with its level re-anchored to the last entry $r_{T_0}$. Under the ARIMA$(1,1,0)$ rule, the null-space part is continued in the same way, while the non-null part is forecast by the ARIMA model.}

The smoothness order $q$ is a structural choice: it fixes the penalty $K_q$, the null space $\mathrm{Null}(K_q)$, and the canonical continuation $G_q^{\mathrm{null}}$, and the decomposition makes its role precise. At $\rho=0$ the metric is $K_q=D_q'D_q$, so $A_3(0)\propto\|D_q\,e^\mathcal{R}_\mathrm{pre}\|_2$ is controlled only if $q$ is large enough that $D_q$ stationarizes the idiosyncratic component of $\mathcal{R}$. An $I(1)$ idiosyncratic stochastic trend is stationarized by $q=1$, whereas an $I(2)$ idiosyncratic stochastic trend is not and requires $q=2$. Raising $q$ to $2$ also enlarges the null space and allows a more flexible specification: since $\mathrm{Null}(K_1)\subset\mathrm{Null}(K_2)$, any approximately affine gap between the treated unit's signal and its donor combination is then absorbed at no metric-distortion cost. These benefits come with two costs. First, the same enlargement forces $G_q^{\mathrm{null}}$ to continue that affine direction: for $q=2$ it extends the line fitted to $P_{0,2}\,r_\mathrm{pre}^\mathrm{oracle}$, and by the null-space continuation property (ref) this continuation sits inside $\mathrm{Term~B}$ at every $\rho$ and can carry an extrapolation bias that grows with the post-treatment window $T_\mathrm{post}$, whereas the $q=1$ continuation carries a level forward and leaves a floor that is flat in the horizon. Second, the spectrum of $K_2$ is much wider than that of $K_1$, which makes the amplification to the high-frequency components much more significant. The metric gap $W_{\rho,q}-P_{\perp,q}$, and hence $A_1$, rises more steeply as $\rho$ falls from $1$.

The cross-validation horizon $h$ does not alter the estimator but selects which error the criterion minimizes. At a short horizon the forecaster's extrapolation bias is typically small, so the criterion is dominated by Term A. At a long horizon the extrapolation bias accumulates and Term B can dominate; the criterion then rewards a large $\hat\rho$ mechanically, which suppresses the time series forecaster. It is worth noting that a larger $h$ leaves less pre-treatment data for cross-validation. In practice, when the pre-treatment window is short, the researcher must balance the post-treatment horizon of interest against the amount of pre-treatment data retained for cross-validation.

These choices are not independent. The cross-validation criterion of Section (ref) selects $\rho$, and optionally $q$ and $\widetilde G_q$ jointly. In practice the structural choices can be guided by what is known about the application. Set $q$ to the smallest order that plausibly stationarizes the raw data. Choose $\widetilde G_q$ conservatively unless the pre-treatment data give clear evidence that a richer forecaster predicts better. Then let cross-validation at the policy-relevant horizon $h$ select $\rho$.

Monte Carlo Evidence

Sections (ref)--(ref) motivate HSC as a soft allocation mechanism, develop its spectral interpretation, and derive a prediction-error decomposition that separates donor matching from residual extrapolation. This section reports a Monte Carlo study that evaluates HSC's finite-sample performance against standard synthetic control estimators and documents how the cross-validated tuning parameter $\hat\rho$ adapts to the underlying data-generating regime. Full details of the data-generating process appear in Appendix (ref).

Design

The data-generating process retains the additive structure used throughout the paper. Untreated potential outcomes follow

equation[equation omitted — 114 chars of source]

where $L_{j,t}=\sum_{k=1}^{3}\Lambda_{j,k}F_{k,t}$ is a low-rank component built from three factors $F_{k,t}$ (one random walk, one ARIMA$(1,1,0)$, one stationary AR(1)), $\kappa \mathcal{E}_{j,t}$ is a unit-specific ARIMA$(1,1,0)$ component whose innovations interpolate between a common shock and an idiosyncratic shock by $\sqrt{\rho_u}u_t^{\mathrm{c}}+\sqrt{1-\rho_u}u_{j,t}^{\mathrm{i}}$, $\varepsilon_{j,t}\sim\mathcal{N}(0,1)$ is stationary noise, $\alpha_j$ is a unit fixed effect, and $\delta_t$ is a time fixed effect. The factor paths $F_{k,t}$, the idiosyncratic component $\mathcal{E}_{j,t}$, the noise $\varepsilon_{j,t}$, and the time fixed effects $\delta_t$ are redrawn in every replication; the factor loadings $\Lambda_{j,k}$ and the unit fixed effects $\alpha_j$ are drawn once and held fixed across replications. The treated unit's loadings are constructed as a sparse Dirichlet convex combination of donor loadings, placing the treated unit inside the donor convex hull. Two scalar parameters govern the persistence structure: $\kappa\in\{0,0.5,1,2\}$ controls the amplitude of the unit-specific stochastic trend, and $\rho_u\in\{0,0.5,1\}$ controls how much of that persistence is shared across units. We replicate every $(\kappa,\rho_u)$ cell $R=500$ times with $T_0=200$ pre-treatment periods, $T_\mathrm{post}=20$ post-treatment periods, and $N_0=50$ donors. The treated unit receives no treatment effect, so the post-period RMSE between the estimated counterfactual and the untreated potential outcome measures predictive accuracy.

We compare five baseline synthetic control estimators against HSC. The baselines are plain SC abadie2010synthetic, SC with an intercept (SC-INT, doudchenko2016balancing), synthetic difference-in-differences arkhangelsky2021synthetic, and two variants of the synthetic business-cycle estimator of shi2025synthetic that differ in the pre-treatment filter used to extract the cyclical component (SBCA-ARIMA and SBCA-Hamilton). Plain SC matches in levels and is therefore vulnerable to spurious matching. SC-INT removes a unit-specific level shift with an intercept, and SDID constructs its unit weights from a level-matching problem with an intercept and a ridge penalty; both still match the residual variation in levels. The SBCA family applies a pre-treatment filter to separate trend from cycle, matches donors on the cycle, and extrapolates the treated trend independently.

For HSC we evaluate four time series forecasters that all satisfy Definition (ref): in each case the forecaster is applied only to the non-null-space component, while the null-space component is continued by the canonical $G_q^{\mathrm{null}}$. The last_constant forecaster carries the last fitted value of the non-null-space component forward as a constant; under $q=1$ this, combined with the canonical constant continuation of the null space, recovers carrying the last fitted value of the residual forward, and under $q=2$ it adds a constant offset to the null-space linear extension. The arima110 forecaster fits an ARIMA$(1,1,0)$ to the non-null-space component, the correctly specified model for the DGP's idiosyncratic ARIMA$(1,1,0)$ stochastic trend. The ar forecaster fits a stationary AR(4) to the non-null-space component, whose forecasts mean-revert. The hamilton forecaster forecasts the non-null-space component with the $h$-step-ahead regression of hamilton2018you. Each forecaster is evaluated at both smoothness orders $q\in\{1,2\}$. The cross-validation horizon is fixed at $h=1$ in this section; Appendix (ref) examines the effect of choosing $h=20$ instead.

HSC ties or improves on baselines across regimes

Figure (ref) reports the post-period RMSE pooled across the 20 post-treatment periods for each of the 12 $(\kappa,\rho_u)$ cells. Within each panel we show the five baseline estimators as single bars, and we show each HSC forecaster as two bars side by side: the lighter bar reports the $q=1$ result and the darker bar reports $q=2$.

figure[figure omitted — 828 chars of source]

Three patterns are visible. First, when the component $\mathcal{E}_{j,t}$ is absent or shared across units, all four HSC forecasters tie SC-INT and SDID and substantially improve on plain SC and the SBCA family. The top row of the figure ($\kappa=0$) and the right column ($\rho_u=1$) display this behavior: HSC's pooled RMSE is within a few percent of SC-INT and SDID, plain SC sits well above the others because it cannot remove the heterogeneous unit intercepts, and the SBCA filters strip away common variation that the other methods can match. Second, when the $\mathcal{E}_{j,t}$ component is present and partially or fully idiosyncratic ($\kappa\ge 0.5$, $\rho_u\le 0.5$), HSC with $q=1$ delivers the lowest pooled RMSE in most cells. The margin between HSC and SC-INT or SDID grows with $\kappa$ and is largest in the corner with the most idiosyncratic drift. The SBCA family performs well when $\mathcal{E}_{j,t}$ is completely idiosyncratic but becomes worse when $\mathcal{E}_{j,t}$ is partially shared. Third, the difference in performance between time series forecasters is visible here; the last_constant and arima110 forecasters perform better than the ar and hamilton forecasters. The $q=1$ HSC also performs better than the $q=2$ configuration. These differences in configurations reflect the design of the DGP, as the idiosyncratic stochastic trend is an ARIMA(1,1,0) model and $q=1$ suffices to control the spurious matching. Appendix (ref) shows that this pooled ranking holds horizon by horizon for the strongest configurations, the constant carry-forward and the $q=1$ ARIMA$(1,1,0)$ forecaster, and Appendix (ref) attributes the pooled advantage to a variance reduction that more than offsets a small bias penalty.

Cross-validation adapts to the regime

The argument in Section (ref) predicts that the optimal $\rho$ depends on the data: when the stochastic trend is shared across units, level matching identifies donor weights well and $\rho$ near one is optimal; when the stochastic trend is idiosyncratic, the donor pool cannot reproduce it and filtering it out by pushing $\rho$ toward zero is preferred. Figure (ref) reports the distribution of the cross-validated $\hat\rho$ across the same $(\kappa,\rho_u)$ grid for the four HSC forecasters at both smoothness orders.

figure[figure omitted — 612 chars of source]

The cross-validated selection behaves as the theory predicts. When the idiosyncratic stochastic trend is absent ($\kappa=0$, the top row of the figure or $\rho_u=1$, the right column of the figure) the distribution of $\hat\rho$ almost piles up at one for every forecaster and both smoothness orders: with no idiosyncratic stochastic trend component to filter out, the CV objective rewards matching in levels. When the stochastic trend is present and at least partially unit-specific ($\kappa\ge 1$, $\rho_u\le 0.5$, bottom-left region of the figure), $\hat\rho$ shifts downward: medians fall to around $0.5$ at $\kappa=1$ and near zero at $\kappa=2$ with $\rho_u=0$, with substantial dispersion across replications. The shift happens for all four forecasters and at both smoothness orders, confirming that the rolling-origin CV identifies the correct allocation between donor matching and time series forecaster directly from the data. Two further properties of the cross-validated allocation are deferred to the appendix: Appendix (ref) reports how $\hat\rho$ and post-period RMSE shift when the CV horizon $h$ is extended from one to twenty.

Empirical Application: The 1997 Handover of Hong Kong

We illustrate HSC on the study of the per-capita GDP after the 1997 return of Hong Kong to Chinese sovereignty. The example was introduced by hsiao2012panel and revisited by shi2025synthetic. hsiao2012panel difference the data to a stationary growth-rate outcome and select a small set of geographically and economically proximate donors, including mainland China and Hong Kong's Asian trading partners, by an information criterion, fitting an unrestricted regression of the treated series on the selected donors. shi2025synthetic instead work with the nonstationary level of annual per-capita GDP, restrict the donor pool to developed economies with comparable long-run growth, and impose a hard separation between a treated-unit trend forecast from Hong Kong's own history and a donor-matched cyclical component. Throughout we use the sign convention $\hat\tau_t = Y_{1t}-\hat Y_{1t}(0)$, so a negative value means observed Hong Kong GDP lies below the estimated no-handover counterfactual. The main text uses the annual data of shi2025synthetic so that the comparison with the most closely related estimator is exact; Appendix (ref) reports robustness to the cross-validation horizon and to the hsiao2012panel geographic-neighbour donor pool.

Data and cross-validated configuration

The panel is the one assembled by shi2025synthetic: annual real per-capita GDP for Hong Kong and eleven developed donor economies, including Australia, Austria, Canada, Denmark, France, Germany, Italy, Korea, the Netherlands, New Zealand, and the United States over $1961$--$2003$. The treatment year is $1997$, giving $T_0=36$ pre-treatment years ($1961$--$1996$) and $T_\mathrm{post}=7$ post-treatment years ($1997$--$2003$); the United Kingdom and mainland China are excluded as parties directly involved in the handover, and economies exposed to the $1997$--$98$ Asian financial crisis or with heterogeneous welfare-state structures are excluded, following shi2025synthetic. The panel-data implementation of hsiao2012panel, which instead selects geographically proximate Asian donors by an information criterion on differenced data, is examined as a robustness check in Appendix (ref).

We evaluate HSC under four configurations: the last_constant and ARIMA$(1,1,0)$ forecasters, each at roughness orders $q\in\{1,2\}$, which are the configurations that performed well in the Monte Carlo study of Section (ref). The tuning parameter $\hat\rho$ is selected by rolling-origin cross-validation with one-step-ahead horizon ($h=1$), a $21$-point $\rho$-grid, and the SDID-style ridge $\zeta=T_\mathrm{post}^{1/4}\hat\sigma_{\Delta X}$. Figure (ref) reports the cross-validated mean squared prediction error along the $\rho$-grid for the four configurations. The cross-validation selects an interior optimum: the best configuration is ARIMA$(1,1,0)$ at $q=1$ with $\hat\rho=0.11$.

figure[figure omitted — 873 chars of source]

Counterfactual comparison

Figure (ref) overlays the cross-validation-selected HSC counterfactual with those of plain synthetic control abadie2010synthetic, synthetic control with an intercept doudchenko2016balancing, synthetic difference-in-differences arkhangelsky2021synthetic, and the synthetic business-cycle estimator of shi2025synthetic (SBCA-Hamilton), together with observed Hong Kong GDP. The estimators fall into three groups. The HSC counterfactual tracks observed Hong Kong closely throughout the post-treatment window, reaching about \$30{,}000 by $2003$ against an observed \$28{,}100, an implied $\hat\tau_{2003}\approx-\$1{,}900$. SC, SC-INT, and SDID drift moderately above the observed series. SBCA-Hamilton diverges sharply: because it forecasts Hong Kong's post-$1997$ trend by a recursive linear projection of its own pre-$1997$ history, and Hong Kong's pre-handover growth was unusually steep, that projection rises to roughly \$36{,}100 by $2003$, an implausibly large effect.

figure[figure omitted — 935 chars of source]

Donor-weight diversification

Figure (ref) compares the donor weights that HSC, SDID, SC-INT, and SBCA-Hamilton assign across the eleven donors. The contrast is stark. HSC distributes weight broadly across all eleven economies, with no single weight exceeding $0.19$ (the largest are Korea $0.18$, Germany $0.14$, the United States $0.13$, and Italy $0.11$). SC-INT collapses onto a corner solution, placing $0.91$ on the United States and $0.09$ on Korea. SBCA-Hamilton concentrates on four donors (Italy $0.43$, Germany $0.25$, Korea $0.18$, the United States $0.09$). SDID is intermediate: its ridge penalty de-concentrates the weights relative to SC-INT. The United States weight falls from $0.91$ to $0.56$ and mass spreads to Denmark ($0.23$), Korea, and Germany.

figure[figure omitted — 780 chars of source]

Out-of-sample accuracy

Pre-treatment fit cannot discriminate among these estimators, because each minimizes a different in-sample criterion. We therefore evaluate every method by the same rolling-origin, one-step-ahead cross-validated MSPE used to select $\hat\rho$. HSC is the most accurate method by a wide margin: at $h=1$ the selected configuration attains a CV-MSPE of $4.8\times10^{5}$, and all four HSC configurations ($4.8$--$5.4\times10^{5}$) fall below every competing estimator---SBCA-Hamilton ($1.2\times10^{6}$), SDID ($1.5\times10^{6}$), SC-INT ($3.6\times10^{6}$), and plain SC ($9.0\times10^{6}$). HSC thus improves on the synthetic business-cycle estimator by roughly a factor of two and a half, and on the level-matching estimators by one to two orders of magnitude, on a criterion that uses only pre-treatment data. Appendix (ref) shows that this ranking is preserved when the cross-validation horizon is lengthened to $h=4$ and when the donor pool is replaced by the geographic-neighbour pool of hsiao2012panel.

Conclusion

Harmonic synthetic control (HSC) addresses counterfactual estimation when untreated outcomes may contain both shared and idiosyncratic stochastic trends, a regime that the researcher cannot reliably distinguish ex ante. Instead of committing in advance to matching in the raw level or to differencing before matching, HSC introduces a treated-unit-specific smooth component and a single tuning parameter that rolling-origin cross-validation uses to allocate predictive responsibility between donor matching and time series forecaster. The spectral interpretation, the prediction-error decomposition, the Monte Carlo evidence, and the Hong Kong application all point to the same conclusion: a soft, data-driven allocation is more robust and can adapt to different regimes.

The present paper develops and evaluates the HSC point estimator; formal uncertainty quantification is the natural next step. A prediction interval for the HSC counterfactual should combine donor-weight estimation uncertainty with out-of-sample forecast-error calibration for the smooth component, extending the synthetic-control prediction-interval framework of cattaneo2021prediction to the soft-allocation setting. Because that construction rests on additional assumptions beyond those required for the point estimator, we leave it to future works.

\onehalfspacing