The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
95,998 characters
Degrees of Freedom and Information Criteria for the Synthetic Control Method
\lstset{
basicstyle=\ttfamily\small,
numbers=left,
keywordstyle= \color{ blue!70},commentstyle=\color{red!50!green!50!blue!50}
}
\title{Degrees of Freedom and Information Criteria for \\the Synthetic Control
Method\thanks{First draft: December 16, 2019. We thank Alberto Abadie, Dmitry Arkhangelsky, Guanglei Hong, Kaspar Wuthrich, as well as the editor and anonymous referees for insightful comments and questions. The first two authors are in alphabetical order and made equal contributions.}}
\author{
Guillaume A.~Pouliot
\thanks{Department of Economics, Rice University,
6100 Main St, Houston, TX 77005, United States. Email: [email removed].}
\and Zhen Xie
\thanks{Department of Economics, Northwestern University,
2211 Campus Drive, Evanston, IL 60208. Email: [email removed].}
\and Ziyi Liu
\thanks{Haas School of Business, University of California, Berkeley,
2220 Piedmont Ave, Berkeley, CA 94720. Email: [email removed].}
}
\date{\today}
\maketitle
\vspace{-0.2in}
\begin{abstract}
We provide an analytical characterization of the model flexibility of the synthetic control method (SCM) in the familiar form of degrees of freedom. We obtain estimable information criteria, which may be used to circumvent cross-validation when selecting either the tuning parameter in
penalized variants of SCM or the weighting matrix in the SCM with covariates. We assess the impact of car license rationing in Tianjin; while a natural match is available, both it and other donors are noisy, inviting the use of SCM to average over approximately matching donors. The very large number of candidate donors calls for
penalized variants of SCM and we observe that model selection using information criteria outperforms that based on cross-validation.
\end{abstract}
\noindent\textbf{Keywords}: Synthetic Controls, Model Selection, Information Criteria, Degrees
of Freedom, Lagrange Multiplier Theory, Chinese Automotive Industry.\\[4pt]
\newpage
\section{Introduction}
The synthetic control method has become a standard regression tool
in economics, political science, and a handful of other fields. See
\citeA{abadie2021using} for a recent
survey and pedagogical introduction. As such, methodological research
has endeavored to append the synthetic control estimator with standard
regression output. This includes quality of fit criteria, confidence intervals, $p$-values,
etc \cite{abadie2010synthetic,chernozhukov2018t, cattaneo2021prediction}.
This paper endeavors to do just that and delivers the degrees of freedom and information criteria for the synthetic control method. As such, our motivation, albeit somewhat jejune, is commensurately uncontroversial.
Our first contribution is to produce
the degrees of freedom. We find that the degrees
of freedom, or the effective number of estimated parameters, of the synthetic
control method without covariates is \emph{one less than the expected
number of donors having nonzero estimated coefficients}. We
obtain a more general result that covers the case with covariates.
Our motivation for producing this as-of-yet-unavailable statistic
is two-fold. First, the question ``does the synthetic control method
overfit?'' is, as we argue below, non-trivial and important.
The degrees of freedom offer a clear and intuitive answer to that
question.
We find that SCM does \emph{not} overfit in most of the seminal
applications we revisited. It does however overfit in ``high-dimensional" applications such as the one we investigate in this paper.
Second, we produce information criteria for synthetic control methods.
We find that, in applications where the number of donors is large relative to the number of pre-treatment observations,
overfitting
(through model selection flexibility) does arise and regularization
is required.
The short pre-treatment series, relative to the number of donors, makes it such that cross-validation
strategies
can perform poorly. An information
criterion, which assesses the out-of-sample performance of the estimator
evaluated on the entire pre-treatment data, is expected to do better.
Equipped with closed-form expressions for the degrees of freedom
that have sample analogs, we produce just such information criteria
and indeed observe (in simulation and placebo cases) that they outperform
cross-validation in terms of producing accurate counterfactual and treatment effect estimates.
Importantly, the produced information criteria can also be used to select the weighting matrix in the synthetic control problem with
covariates, see problem formulation (\ref{eq:SCwithCOVbeginning})-(\ref{eq:SCwithCOVend}).
\subsubsection*{Degrees of freedom}
Degrees of freedom expressions are inherently interesting for the
synthetic control method. Indeed, consider two of the more striking
--and attractive-- features of the synthetic control method, typically displayed in its standard output. We reproduce in \Cref{figure1}
the output of \citeA{abadie2010synthetic} as an example.
First, the fitted regression coefficient estimates typically exhibit
substantial sparsity; the synthetic control is a linear combination
of a few donors, with many donors getting an estimated weight of zero.
Second, the observed and fitted paths, before treatment, often suggest
a high in-sample fit, as is the case in \Cref{figure1}.
These two features capture our motivation for inquiring into the degrees
of freedom of the synthetic control method. The high degree of sparsity in the estimated regression coefficients suggests the possibility
of extensive implicit model selection. Given that the series being fitted are often short, one may worry that this additional model flexibility brings about overfitting, and that this is in fact what explains the surprising quality of the in-sample fit, hence casting doubt on the
quality of the counterfactual and treatment effect estimates.
\begin{figure}
\centering
\begin{center}
\begin{minipage}[t]{0.30\textwidth}
\centering
\resizebox{1.45\textwidth}{!}{
\begin{tabular}{@{} l r l r @{}}
\toprule
\hline
State & Weight & State & Weight \\
\midrule
Alabama & 0 & Montana & 0 \\
Alaska & -- & Nebraska & 0 \\
Arizona & -- & Nevada & 0.062 \\
Arkansas & 0 & New Hampshire & 0.001 \\
Colorado & 0 & New Jersey & -- \\
Connecticut & 0 & New Mexico & 0 \\
Delaware & 0 & New York & -- \\
District of Columbia & -- & North Carolina & 0 \\
Florida & -- & North Dakota & 0 \\
Georgia & 0 & Ohio & 0 \\
Hawaii & -- & Oklahoma & 0 \\
Idaho & 0 & Oregon & -- \\
Illinois & 0 & Pennsylvania & 0 \\
Indiana & 0.175& Rhode Island & 0 \\
Iowa & 0 & South Carolina & 0.343 \\
Kansas & 0 & South Dakota & 0 \\
Kentucky & 0 & Tennessee & 0 \\
Louisiana& 0 & Texas & 0.182 \\
Maine & 0.236& Utah & 0 \\
Maryland & -- & Vermont & 0 \\
Massachusetts & -- & Virginia & 0 \\
Michigan & -- & Washington & -- \\
Minnesota& 0 & West Virginia & 0 \\
Mississippi & 0 & Wisconsin & 0 \\
Missouri & 0 & Wyoming & 0 \\
\bottomrule
\end{tabular}
}
\end{minipage}
\hfill
\begin{minipage}[t]{0.53\textwidth}
\vspace{-15em}
\centering
\includegraphics[scale=0.64]{graphs/California_tobacco/main.pdf}
\end{minipage}
\end{center}
\caption{\emph{Synthetic control output of California Proposition 99 investigation. The degrees of freedom estimate is 5. The left-hand side and right-hand side panels are, respectively, our replications of Table 2 and Figure 2 of \protect\citeA{abadie2010synthetic}.}}
\label{figure1}
\end{figure}
The question \emph{``does synthetic control overfit?''} thus stands open. On the one hand, the best subset selection method with sample
size and number of active independent variables typical of synthetic
control applications would be expected to overfit (see \Cref{figure2}).
On the other hand, placebo exercises in which one forecasts non-treated
series and compares the forecasted to the realized series often suggest
a reasonable quality of forecast (e.g., \citeNP{abadie2010synthetic}).
This motivates the development of an analytical measure of the model
flexibility of the method, the most compelling of which is, in our
opinion, its degrees of freedom.
The key idea is to produce an estimable expression for the degrees of freedom by applying Stein's Lemma, see Section \ref{section3}. To the best of our knowledge, this was first done by \citeA{meyer2000degrees}. The closest paper to ours is \citeA{zou2007degrees}, which uses Stein's Lemma to compute the degrees of freedom of the lasso (see Figure \ref{figure2}). Their result and the similarities between the lasso and SCM estimation problems inspired this analysis.
\subsubsection*{Information criteria}
The wide applicability of the synthetic control method has brought
applied analysts to consider ``high-dimensional'' applications --i.e.,
with many donors-- even though the training, or pre-treatment, period may be short
relative to the number of donors.
\footnote{For instance, \citeA{cavallo2013catastrophic}
study the causal impact of catastrophic natural disasters on economic
growth, they have 196 donors and fewer than 40 training periods in
each regression. \citeA{bifulco2017using} study the
impact of an education program on district enrollments and graduation
rates, they have 275 donors and 10 training periods. \citeA{bohn2014did}
study the effect of Arizona\textquoteright s 2007 Legal
Arizona Workers Act (LAWA) on the proportion of the state\textquoteright
s population, they have 46 donors and 9 training periods. \citeA{pieters2016effect}
study the effect of democratization on child mortality, they have
24 untreated units and 10 training periods. \citeA{heersink2017disasters}
study the effect of natural disasters on politicians\textquoteright{}
electoral fortunes, they have 100 donors and 29 training periods.
\citeA{peri2019labor} study the effects of the
Mariel Boatlift on wages and employment, they have 43 donors and six
training periods. \citeA{billmeier2013assessing}
study the impact of economic liberalization on the real GDP per capita,
they have 180 donors and 38 training periods.} Such applications
of course raise \emph{a priori} concerns about overfitting,
and these have brought about the development of penalized
methods,
all of which require the specification of a tuning parameter
(\citeNP{abadie2021penalized, athey2021matrix, doudchenko2016balancing}).
Typically, such tuning parameters are estimated by cross-validation,
with different specific algorithms preferred by different authors,
see \Cref{table1}.
\begin{table}[H]\centering
\begin{tabular}{
>{\centering\arraybackslash}p{0.28\linewidth}
>{\centering\arraybackslash}p{0.67\linewidth}
}
\toprule
CV method & Reference\\
\midrule
Rolling window & \citeA{kellogg2020combining}
\\[8pt]
\multirow{3}{*}{\parbox{\linewidth}{\centering Pre-Intervention Holdout Validation on the Treated}} & \citeA{abadie2015comparative}, \citeA{xu2017generalized},\\
& \citeA{abadie2021penalized},\\
& \citeA{ben2021augmented}
\\[8pt]
\multirow{2}{*}{\parbox{\linewidth}{\centering Leave-one-out Validation on the Untreated}} & \citeA{abadie2021penalized}\\
\\
\bottomrule
\end{tabular}
\caption{\emph{Survey of cross-validation procedures suggested in the synthetic
controls literature.}}
\label{table1}
\end{table}
Cross-validation is also called upon in synthetic control applications
to estimate the weighting matrix of the inner problem in the synthetic
control method with covariates. See (\ref{eq:SCwithCOVbeginning})-(\ref{eq:SCwithCOVend})
below for details about the construction of synthetic controls using
covariates.
However, cross-validation is often ill-suited for model selection
with the synthetic control method.
On the one hand, ``pre-intervention holdout validation on the treated''
is a data-splitting procedure; the estimator is trained on the first
half of the pre-treatment data, and the second half is used as a test
set. As such, we may expect it to be severely biased with short pre-treatment series; conceptually, the synthetic control method trained on, say, half of the pre-treatment data may behave very differently from one trained on the whole pre-treatment data. On the other hand, the ``leave-one-out validation on the untreated'' approach treats a donor as a placebo to-be-treated unit, and uses the post-treatment data as a test set.
It ignores the true to-be-treated unit, and averages over the post-treatment
forecast error of a selection of such placebo to-be-treated units.
It relies on the assumption that, in the absence of a treatment effect,
the conditional distribution of any donor --perhaps after some preselection--
given the other donors is the same as that of the to-be-treated unit
given the donors.
This is a strong, additional assumption to make about the underlying data generating process.
As investigated in \citeA{kellogg2020combining}, the pre-intervention holdout approach may produce a misleading estimate
of the prediction mean-squared error. \citeA{kellogg2020combining} instead propose to use rolling-window cross-validation which they find to perform best, at least for the estimator they propose.
As detailed below, these procedures remain data hungry, and we are well motivated to consider out-of-sample error estimation methods that
use information from the entire pre-treatment data without implicitly assuming symmetry between the donor and the to-be-treated unit.
A classical alternative to cross-validation is to compare models using
an information criterion (\citeNP{claeskens2008model}).
Information criteria append to the in-sample loss a penalty for model
flexibility, thus allowing for the comparison of models of different
flexibility while training them on the entire pre-treatment data.
The key point is that while both information criteria (may they be AIC, BIC, or other) and cross-validation techniques estimate out-of-sample fit, information criteria do so using all available data, while cross-validation requires splitting the data. This brings about a tradeoff between the stronger assumptions required for the validity of information criteria and the larger sample required for the accuracy of cross-validation. The two approaches thus being complementary, a complete toolkit for a regression method should include both.
It is worthwhile to note that the aforementioned cross-validation procedures are not only
data-hungry, they also require
tuning from the user. The rolling window and pre-intervention holdout
cross-validation approaches require careful tuning of the size of the
training and test sets. Likewise, for leave-one-out validation on the
untreated, the pool of placebo units (those ``similar enough'' to
the treated unit) must be selected; this can be delicate, as illustrated
in \citeA{abadie2010synthetic}. Meanwhile, the information
criteria approach requires no tuning and our publicly available code
can be used directly.
The challenge in producing information criteria is to give a computable
expression for a penalty term capturing model flexibility. As detailed
in \Cref{section3} below, this exercise is intimately
tied to the estimation of degrees of freedom.
The intellectual history of this approach to producing information criteria can be summarized as follows. \citeA{stein1981estimation} proved the now famous Stein's Lemma, displayed below in (\ref{eq:SteinLemma}), and used it to produce what is now called Stein's unbiased risk estimate (SURE), displayed below in (\ref{eq:homoskedasticSURE}). This implicitly defines a general formulation for degrees of freedom as part of the penalty term in SURE, when the latter is thought of as an information criterion. The literature attributes to \citeA{efron2004estimation} and \citeA{hastie1990generalized} the definition of degrees of freedom in terms of covariance between observed and fitted values, which matches that deduced from the penalty term of SURE, before applying Stein's lemma. \citeA{meyer2000degrees} used Stein's lemma to compute the degrees of freedom for shape restricted regression, setting the stage for other such applications. Amongst those, we find the degrees of freedom of the lasso, a result first given in \citeA{zou2007degrees}, further generalized in \citeA{tibshirani2012degrees}, and which has inspired the analysis herein.
General results producing degrees of freedom estimates for classes of models are given in \citeA{kato2009degrees} and \citeA{chen2020degrees}.
\subsubsection*{Impact of rationing on the demand for cars}
The aforementioned methodological developments are of general interest
but were, for the authors, motivated by the analysis of the market
for new automobiles in China after the introduction of rationing.
We investigate model-specific changes in sales in Tianjin following the introduction of a lottery-auction hybrid for the distribution
of car licenses. In order to build counterfactual --as if there
had been no rationing-- time series for each individual car model,
the synthetic control method turns out to be a natural choice but falls prey to overfitting, and the penalized synthetic control method (\citeNP{abadie2021penalized})
is more reliable. However, cross-validation performs poorly for tuning
parameter selection and we are able to carry out a more robust analysis
using our information criteria instead. We present model-specific
treatment effects of the policy, thus allowing for a detailed study
of the market for new cars in Tianjin.
This is a somewhat novel, or at least uncommon, application of the synthetic
control method. While a natural match is available --the same car
model in an untreated city-- it is very noisy. By instead averaging
over many approximate matches, we can produce a synthetic control
having a much smaller variance than the natural donor does, at a relatively
small cost in bias.
\subsubsection*{Reduced form framework, target of inference, and notation}
The typical notation in the synthetic control and in the degrees of
freedom literature is somewhat at odds. We clarify the connection.
The synthetic control literature uses notation
appealing to the potential outcomes framework with, say, $(Y_{j,t}(0),Y_{j,t}(1))$,
$j=0,1,...,p$, being the potential outcomes under control and treatment
of unit $j$ at time $t$. Some units will remain untreated over time
while some will be treated after some time period $T^{*}$, i.e., $Y_{j,t}(1)$
is observed instead of $Y_{j,t}(0)$ for such $j$'s for $t>T^{*}$.\footnote{This is the canonical framework, scattered treatment times may also
be considered.}
The object of interest is some treatment effect $\tau=Y_{j,t}(1)-Y_{j,t}(0)$ for the --here, unique-- treated unit $j=0$ at some time $t>T^{*}$.
The synthetic control estimator being a constrained least-squares
estimator, its implicit target is the constrained best linear predictor
and, under correct specification of the regression function, the conditional
expectation. Specifically, under correct specification, the population
estimand is $Y_{j,t}(1)-E\left[\left.Y_{j,t}(0)\right|Y_{1,t}(0),\dots,Y_{p,t}(0)\right]$.
Since only a single $Y_{j,t}(1)$ is observed, quality of inference
about $\tau$ is assessed conditionally on $Y_{j,t}(1)$ and is thus
commensurate to quality of inference about $E\left[\left.Y_{j,t}(0)\right|Y_{1,t}(0),\dots,Y_{p,t}(0)\right]$.
Furthermore, since only pre-treatment data is involved in the study
of this object, the notation more typical of the degrees of freedom
literature suffices, and in fact will prove quite handy and natural.
Bridging to the notation we will use, we define the vectors $\mathbf{Y}=\left(Y_{0,1}(0),\dots,Y_{0,T^{*}}(0)\right)^{T}$
, $\mathbf{X}_{j}=\left(Y_{j,1}(0),\dots,Y_{j,T^{*}}(0)\right)^{T}$, $j=1,...,p$, and $\mathbf{X}=\left(\mathbf{X}_{1},\dots,\mathbf{X}_{p}\right)$.
Letting $n=T^{*}$, we recover the conventional $\mathbf{X}\in\mathbb{R}^{n \times p}$.
We use $\mathbf{X}_{J}$, $J\subset\{1,\dots,p\}$, to designate the submatrix
made of the columns designated by $J$. We use $\mathbf{1}_{p}$, $p\in\mathbb{N}$,
to designate the vertical vector of ones of length $p$. For any vector
$\mathbf{w}\in\mathbb{R}^{n}$, we use $\left\Vert \mathbf{w}\right\Vert _{2}$
to indicate its $\ell_{2}$ norm and for any set $\mathcal{B}$, we
use $\left|\mathcal{B}\right|$ to indicate its cardinality.
\bigskip{}
Our methodological contributions can be succinctly summarized as follows.
We produce information criteria for synthetic control methods.
In the case of the penalized synthetic control method, this allows for the selection of the tuning parameter controlling model selection.
In the classical case of the synthetic control method with covariates, this allows for the selection of the weighting matrix.
Importantly, model
selection based on the information criteria circumvents cross-validation,
which we argue and demonstrate can be misleading in practice.
We furthermore produce closed form expressions for the degrees of freedom for the synthetic control method, with or without covariates, as well as for penalized versions.
These have sample analogs and may be reported alongside standard output.
See Figures \ref{figure1} and \ref{fig:one model fits} for a suggestion of modified output.
The remainder of the article is divided as follows. \Cref{section2}
defines the synthetic control method and several penalized
extensions. \Cref{section3} produces the information
criteria and degrees of freedom estimates for the different methods.
\Cref{section4} uses some of the herein developed
methodology to study the impact of license rationing on the sales
of individual car models in Tianjin. \Cref{section5}
discusses and concludes.
\section{The Synthetic Control Method}
\label{section2}
We are interested in the classical synthetic control method, and particularly
interested in its penalized extensions.
\subsection{The Synthetic Control Method Without Covariates}
Given an observed series $\mathbf{Y}\in\mathbb{R}^{n}$ and $p$ observed
series from ``donors'' collected as a matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$,
the vector of optimal donor weights $\hat{\beta} \in \mathbb{R}^{p}$ is the solution
of the optimization problem
\begin{equation}
\min_{\beta \in \mathbb{R}^{p}}\left\Vert \mathbf{Y}-\mathbf{X}\beta\right\Vert^{2} _{2}\label{eq:SCwithoutCOVbeginning}
\end{equation}
subject to
\begin{equation}
\mathbf{1}_{p}^{T}\beta=1,\ \beta\ge0,\label{SCwithoutCOVend}
\end{equation}
where the inequality applies pointwise to vectors and $\mathbf{1}_{p} \in \mathbb{R}^{p}$ is the $p$-tuple whose entries are all 1.
The standard regression output of the synthetic control method contains
a table of regression coefficient estimates, as well as a plot of
the observed series before and after treatment, overlaid with the
treated unit's series fitted by the synthetic control method before treatment,
and forecasted after treatment. Because the forecasted series is a
function of untreated units --called ``donors''-- estimated on
pre-treatment data, it is interpreted as forecasting a non-treated
counterfactual --a ``synthetic control''-- for the treated unit.
The difference between the realized and forecasted series provides
an estimate of the treatment effect.
\subsection{The Synthetic Control Method With Covariates}
In some cases, $n_{\mathrm{cov}}$ ``covariate'' variables are believed to
satisfy --at least approximately-- the same linear relationship
as the aforementioned series, and the coefficient $\beta$ is constrained
to be a synthetic control solution for these covariates. Let the matrix
$\mathbf{D}\in\mathbb{R}^{n_{\mathrm{cov}}\times p}$ and the vector
$\mathbf{Z}\in\mathbb{R}^{n_{\mathrm{cov}}}$ collect the independent
and dependent covariate variables, respectively. The optimal solution
$\hat{\beta}$ is then obtained by solving the two-level
program
\begin{equation}
\min_{\beta}\left\Vert \mathbf{Y}-\mathbf{X}\beta\right\Vert_{2}^2
\label{eq:SCwithCOVbeginning}
\end{equation}
subject to
\begin{equation}
\beta\in \underset{\beta'\in\mathbb{S}}{\arg\min}\ \left\Vert \mathbf{Z}-\mathbf{D}\beta'\right\Vert_{V}^2,
\label{eq:SCwithCOVend}
\end{equation}
where $\mathbb{S}=\left\{ \beta\in\mathbb{R}^p: \mathbf{1}_p^T\beta=1, \beta\geq 0 \right\}$, $\left\Vert \bm{a}\right\Vert _{V}=(\bm{a}^{T}V\bm{a})^{1/2}$ for $V\succ 0$,
and $\left\Vert {\bm a}\right\Vert _{2}=({\bm a}^T{\bm a})^{1/2}$ for conformable vectors ${\bm a}$.
The diagonal matrix $V$ effectively weighs the ``observations''
of the inner regression problem. As such, it may be considered as a hyperparameter in the above formulation and may be selected using an out-of-sample criterion. Indeed, we may pick $V$ by cross-validation, solving at
each iteration the problem (\ref{eq:SCwithCOVbeginning})-(\ref{eq:SCwithCOVend})
with $V$ fixed, and searching for a ``good" $V$ matrix. This is a suggestion of \citeA{abadie2015comparative}. As detailed below, our proposed information criteria will be an attractive alternative to cross-validation.
\subsection{The Penalized Synthetic Control Method}
In settings where there are many donors relative to the number of
pre-treatment time periods, one may be concerned about overfitting.
Indeed, such a relatively large number of donors makes it more likely
to find a linear combination of donors that closely matches the to-be-treated
unit, even though they in fact have little predictive power.
In particular,
a collection of very ``far away'' donors may give a tight in-sample
fit, even though we may have \emph{a priori} knowledge or belief that such
``far away'' donors are more likely to be poor individual matches and are expected
as such to make for a poor synthetic control.
We may therefore want
to limit the choice of donors by relying on the conventional prior that donors more ``similar'' to the to-be-treated unit are more reliable in the specific sense that they offer a
better out-of-sample forecast.
A natural way to implement such a prior belief for frequentist estimation
is to add to the objective function a shrinkage term penalizing coefficients that put more weight on \emph{a priori} bad matches. This was suggested
in \citeA{abadie2010synthetic} and implemented in \citeA{abadie2021penalized}.
The estimator $\hat{\beta}_{\mathrm{pen}}$ is defined as the minimizer
of the penalized synthetic control method (PSCM) problem
\begin{equation}
\min_{\beta}\left\Vert \mathbf{Y}-\mathbf{X}\beta\right\Vert _{2}^{2}+\lambda\sum_{j=1}^{p}\beta_{j}\left\Vert \mathbf{Y}-\mathbf{X}_{j}\right\Vert _{2}^{2}\label{eq:PenalizedCSbeginning}
\end{equation}
subject to
\begin{equation}
\mathbf{1}_{p}^{T}\beta=1,\ \beta\ge0,\label{eq:PenalizedSCend}
\end{equation}
where $\mathbf{X}_{j}$ is the $j^{\mathrm{th}}$ column of $\mathbf{X}$.
Note that $\lambda>0$ must be selected by the user.
\subsection{Constrained Ridge SCM}
A different and more flexible approach to regularizing the synthetic control method is proposed in \citeA{arkhangelsky2021synthetic}. The optimization problem is
\begin{equation}
\min_{\beta_{0},\beta} \ \left\Vert \mathbf{Y}-\beta_{0}\mathbf{1}_{n}-\mathbf{X}\beta\right\Vert _{2}^{2}+\lambda\left\Vert \beta\right\Vert _{2}^{2}
\label{eq:CR objective}
\end{equation}
subject to
\begin{equation}
\mathbf{1}_p^{T}\beta=1,\ \beta\ge0,
\label{eq:CR constraint}
\end{equation}
where $\beta_0\in\mathbb{R}$ and $\lambda > 0$ is a regularization parameter that must be selected by the user.
First, this modified synthetic control method includes an intercept term.
As pointed out in \citeA{doudchenko2016balancing}, it allows for a systematic additive difference between the treated unit and control units, which is an important feature of the standard difference-in-difference strategy. \citeA{ferman2021synthetic} also show that including an intercept term can improve the estimator in terms of bias and variance when the pre-treatment fit is not perfect.
Second, following \citeA{doudchenko2016balancing}, adding an $\ell_2$ regularization penalty can increase the dispersion and ensure the uniqueness of the synthetic control weights.
\subsection{Elastic Net SCM}
Another popular cousin of the SCM is the elastic net variant of \citeA{doudchenko2016balancing}. Their modified problem is
\begin{equation}
\min_{\beta_{0},\beta} \ \left\Vert \mathbf{Y}-\beta_{0}\mathbf{1}_n-\mathbf{X}\beta\right\Vert _{2}^{2}+\lambda_{1}\left\Vert \beta\right\Vert _{1}+\lambda_{2}\left\Vert \beta\right\Vert _{2}^{2},
\label{eq:EN}
\end{equation}
with $\lambda_1, \lambda_2>0$. Remark that the optimization is not subject to the simplex constraint.
\section{Information Criteria and Degrees of Freedom}
\label{section3}
Because the selection of the tuning parameter $\lambda$ in (\ref{eq:PenalizedCSbeginning})-(\ref{eq:PenalizedSCend}), (\ref{eq:CR objective})-(\ref{eq:CR constraint}),
(\ref{eq:EN}), or of the weighting matrix $V$ in (\ref{eq:SCwithCOVbeginning})-(\ref{eq:SCwithCOVend}),
boils down to comparing models of different flexibility, we want to
use a measure of out-of-sample fit as our criterion.
Our preferred such notion of fit is the error between the model's fitted values and
the target regression function, the conditional expectation.
Specifically, we want to select a model, or tuning parameter, which minimizes the risk
\begin{equation}
\mathcal{R}:=E\left\Vert \hat{\mathbf{Y}}-E\left[\mathbf{Y}\left|\mathbf{X}\right.\right]\right\Vert _{2}^{2},
\end{equation}
where $\hat{\mathbf{Y}}\in\mathbb{R}^{n}$ are the fitted values produced
by the candidate model.
We begin by noting that, since the goal is model selection, we will care about the risk only up to a constant. Note that
\begin{equation}
\mathcal{R} =E\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2E\left[\sum_{i=1}^{n}\operatorname*{Cov}\left(\left.\hat{Y}_{i},Y_{i}\right|\mathbf{X}\right)\right]+\mathrm{const} ,
\label{eq:proprisk}
\end{equation}
where ``$\mathrm{const}$" always stands for some constant that does not depend on the fitted values, and is thus irrelevant for assessment of the quality of fit.
We speak of an \emph{information criterion} to refer to estimates of $\mathcal{R}$ based on the decomposition (\ref{eq:proprisk}), i.e., a measure of in-sample fit plus a penalty term for model flexibility.
The information criteria approach requires an estimable penalty term.
To accomplish this, we rely on the modern theory for degrees of freedom (\citeNP{meyer2000degrees}, \citeNP{zou2007degrees}) based on Stein's lemma and the computation of divergences.
Stein's
lemma --and thus the entire degrees of freedom literature-- requires the Gaussian assumption
\begin{equation}
\mathbf{Y}|\mathbf{X}\sim N\left(E\left[\mathbf{Y}|\mathbf{X}\right],\Sigma_{Y|X}\right).\label{GaussianAssumption}
\end{equation}
This assumption is discussed in detail in \Cref{sec:discussion-Gaussian}.
Assuming (\ref{GaussianAssumption}) with diagonal covariance matrix $\Sigma_{Y|X}$ and almost differentiability
(see \Cref{subsectionA4}), Stein's Lemma (\citeNP{stein1981estimation}) states that
\begin{equation}
\operatorname*{Cov}\left(\left.Y_{i},\hat{Y}_{i} \ \right| \ \mathbf{X}\right)=\sigma_{i}^{2}E_{Y|X}\left[\frac{\partial\hat{Y}_{i}}{\partial Y_{i}}\right],\label{eq:SteinLemma}
\end{equation}
where $\sigma_{i}^{2}=V\left(\left.Y_{i}\right|\mathbf{X}\right)$.
This is a simple yet powerful result.
While the left-hand side of (\ref{eq:SteinLemma}) is intractable,\footnote{This is immediate. Consider the favorable case of unbiased forecasts; computation of the conditional covariance requires knowledge of the conditional expectation $E\left[\mathbf{Y}|\mathbf{X}\right]$,
which is the object we are trying to estimate in the first place.}
the right-hand side may be estimated in closed form for the regression models covered in this article.
Under the Gaussian assumption (\ref{GaussianAssumption}), equation
(\ref{eq:proprisk}) is equivalently expressed as
\begin{equation}
\mathcal{R} = E\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2E\left[\sum_{i=1}^{n}\sigma_{i}^{2}E_{Y|X}\left[\frac{\partial\hat{Y}_{i}}{\partial Y_{i}}\right]\right]+\mathrm{const.}\label{genSURE}
\end{equation}
Remarkably, this closed form has a sample analog as long as we can plug in estimates for the divergence and the variance.
To produce a simple plug-in estimate of (\ref{genSURE}), we assume the data is independently and identically distributed and conditionally homoskedastic such that
\begin{equation}
\sigma_{i}^{2}=\sigma^{2},\ \forall \ i, \label{eq:Ahomoskedastic}
\end{equation}
for some $\sigma^{2}$.
Specifically, supposing (\ref{eq:Ahomoskedastic}) holds and omitting the additive constant, expression (\ref{genSURE})
rewrites as
\begin{equation}
\mathrm{IC} := E\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2\sigma^{2} E \left[\mathrm{df}\left(\hat{\mathbf{Y}}\right)\right],\label{eq:homoskedasticSURE}
\end{equation}
where
\begin{equation}
\label{eq:homoskedastic_df}
\mathrm{df}\left(\hat{\mathbf{Y}}\right)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\operatorname*{Cov}\left(\left.Y_{i},\hat{Y}_{i}\right| \mathbf{X} \right)=\mathrm{Tr}\left(E\left[\left.\nabla\hat{\mathbf{Y}}\right|\mathbf{X}\right]\right),
\end{equation}
and the divergence $\nabla\hat{\mathbf{Y}}$ is the Jacobian of $\hat{\mathbf{Y}}$ with respect to $\mathbf{Y}$.
The information criterion $\mathrm{IC}$ is our main estimand of --proportional--
risk. \emph{The key is to develop theory that produces computable
estimates of the divergence $\nabla\hat{\mathbf{Y}}$}. We do so
in \Cref{subsectionA4}
and produce sample analogs of the degrees of freedom \eqref{eq:homoskedastic_df}
in \Cref{sec:info-criteria}.
Recognizing similarities between the lasso and SCM estimation problems, we recuperate the insight of \citeA{zou2007degrees} and use Stein's Lemma to produce an expression for the degrees of freedom and a risk estimate for the SCM that have sample analogs.
We obtained analytical and computational simplifications by assuming that the data was normally distributed, which is not exact in typical applications. However, the resulting expression appears to be a good estimate of model flexibility in non-Gaussian regimes.
Indeed,
both theory and simulations substantiate the claim of robustness against departure from Gaussianity. See \Cref{sec:simulation}.
The conditional homoskedasticity invoked in (\ref{eq:Ahomoskedastic}) can, however, be more challenging.
Under heavy heteroskedasticity, (\ref{eq:homoskedasticSURE}) may not be a good approximation of (\ref{genSURE}). As detailed below, in such regimes, we recommend the less simple but more robust alternative estimator (\ref{eq:robustIC}).
We remark that under homoskedasticity, the risk expression (\ref{eq:homoskedasticSURE})
relates in an immediate fashion to degrees of freedom as well as to
more familiar expressions for information criteria, such as AIC, BIC
or Mallow's $C_{p}$.
The other classical approach for comparing the quality of fit of different models without falling prey to overfitting is to estimate the fit on held-out data; this is referred to as \emph{cross-validation}.
The cross-validation approach is understood to estimate
$\mathcal{CV}:=E\left\Vert \mathbf{Y}_{\mathrm{new}}-\hat{\mathbf{Y}}_{\mathrm{new}}\right\Vert _{2}^{2},$
which is equal, up to a constant, to
\begin{equation}
E\left\Vert \hat{\mathbf{Y}}_{\mathrm{new}}-E\left[\mathbf{Y}\left|\mathbf{X}_{\mathrm{new}}\right.\right]\right\Vert _{2}^{2},
\end{equation}
where $\left(\mathbf{X}_{\mathrm{new}},\mathbf{Y}_{\mathrm{new}}\right)$
is an observation drawn from the same data generating process as,
but not included in, the data set which the model forecasting $\mathbf{Y}_{\mathrm{new}}$
with $\hat{\mathbf{Y}}_{\mathrm{new}}$ was trained on.\footnote{Explicitly, $E\left\Vert \mathbf{Y}_{\mathrm{new}}-\hat{\mathbf{Y}}_{\mathrm{new}}\right\Vert _{2}^{2}=E\left\Vert E\left[\mathbf{Y}\left|\mathbf{X}_{\mathrm{new}}\right.\right]+\varepsilon_{\mathrm{new}}-\hat{\mathbf{Y}}_{\mathrm{new}}\right\Vert _{2}^{2} = E\left\Vert \hat{\mathbf{Y}}_{\mathrm{new}}-E\left[\mathbf{Y}\left|\mathbf{X}_{\mathrm{new}}\right.\right]\right\Vert _{2}^{2}+\mathrm{const},$
because $\varepsilon_{\mathrm{new}}=\mathbf{Y}_{\mathrm{new}}-E\left[\mathbf{Y}\left|\mathbf{X}_{\mathrm{new}}\right.\right]$
is independent of $\hat{\mathbf{Y}}_{\mathrm{new}}$.}
Of course, both information criteria and cross-validation, as they have their own advantages and disadvantages, are standard model selection tools in modern statistics and econometrics.
Information criteria typically require strong assumptions in order to deliver closed-form expressions. Cross-validation typically involves some form of sample-splitting, and as such may be biased or data hungry.
Finally, remark that using Stein's lemma and computable expressions
for the divergence, we obtain an analytical and estimable expression
for $\mathrm{df}(\hat{\mathbf{Y}})$ in order to quantify analytically
the model flexibility of the synthetic control method in typical applications. In \Cref{sec:df-sc}, we carry out such an inquiry and answer our motivating question ``does the synthetic control method overfit?''.
\subsection{The Degrees of Freedom of the Synthetic Control Methods}\label{sec:df-sc}
Once a tractable analytic expression obtains for the divergence, the
risk estimate and degrees of freedom obtain upon verifying regularity
conditions. However, computing the divergence in specific cases can
be non-trivial. Correspondingly, as \citeA{tibshirani2015stein} point out, a small industry of computing divergences has blossomed.\footnote{For instance, \citeA{meyer2000degrees} compute divergences for convex
constrained regression estimators, \citeA{mukherjee2015degrees} compute divergences for reduced rank regressions, \citeA{candes2013unbiased} compute the degrees of freedom of singular value thresholding and
spectral estimators, \citeA{deledalle2012risk} compute divergences for
singular value thresholding, \citeA{mazumder2020computing} compute divergences
for estimators with rank penalty, \citeA{minami2020degrees} gives the degrees
of freedom of estimators with submodular penalties, and \citeA{chen2020degrees} study least-squares estimators with linear penalties and constraints.}
\Cref{subsectionA4} collects the technical derivations of the divergences. With those in hand, the desired degrees of freedom expressions, and thus the information criteria, readily obtain as corollaries.
The degrees of freedom sometimes obtain in rather elegant closed-form expressions.
Recall that our definition of degrees of freedom for a model and its
fitted values $\hat{\mathbf{Y}}$ is
\[
\mathrm{df}\left(\hat{\mathbf{Y}}\right)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\operatorname*{Cov}\left(\left.Y_{i},\hat{Y}_{i}\right|\mathbf{X}\right),
\]
assuming (\ref{eq:Ahomoskedastic}) holds. This is an attractive measurement of model flexibility.
Indeed, define the noise perturbations $\varepsilon_{i}=Y_{i}-E\left[\left.Y_i\right|\mathbf{X}\right]$,
$i=1,...,n$, then the definitional equivalence $\operatorname*{Cov}(Y_{i},\hat{Y}_{i} \ | \ \mathbf{X})=\operatorname*{Cov}(\varepsilon_{i},\hat{Y}_{i} \ | \ \mathbf{X})$
makes explicit that $\mathrm{df}(\hat{\mathbf{Y}})$ is
measuring how much the fitted values are adapting to noise. A familiar
instance of degrees of freedom is the number of regression coefficients
in ordinary least-squares (OLS).
We indeed obtain that statistic back in the Gaussian setting, as a special case of the more general definition. For instance, for OLS with full column rank design matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$,
under assumptions (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}), application of Stein's lemma delivers
the familiar quantity
\[
\mathrm{df}\left(\hat{\mathbf{Y}}_{\mathrm{ols}}\right)=E\left[\mathrm{Tr}\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\right]=p,
\]
where $\hat{\mathbf{Y}}_{\mathrm{ols}}=\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\mathbf{Y}$.
We are able to quantify the model flexibility, or ``effective
number of coefficients'', of the synthetic control method. Consider
the case of the synthetic control method without covariates.
Quite remarkably, we find that in spite of the implicit-- and sometimes extensive-- model selection carried out by the synthetic control method,
the expected model flexibility remains that of a linear regression with sum-to-one constraint using only the selected donors.
In other words, the implicit model selection is carried out at no additional
cost in degrees of freedom.\footnote{There are different ways to develop intuition for that result. For the lasso, \citeA{tibshirani2015stein} informally speak of the penalization on the value of the coefficients perfectly offsetting the model selection flexibility,
thus producing the same degrees of freedom, in expectation, as if
carrying out ordinary least-squares on the selected model. Problem
(\ref{eq:SCwithoutCOVbeginning})-(\ref{SCwithoutCOVend}) has an
equivalent representation as a lasso problem for non-negative coefficients and for a specific tuning parameter, where the coefficients are furthermore constrained to sum to one, thus explaining the one fewer degrees of freedom.}
\begin{proposition} [Degrees of Freedom of SCM Without Covariates]
Let $\tilde{\mathbf{X}}=(\mathbf{X}^\top,\mathbf{1}_p)^\top\in\mathbb{R}^{(n+1)\times p}$, an augmented donor matrix.
Suppose that $\mathbf{Y}|\mathbf{X}$ follows the probability law stipulated in (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}).
Then,
\begin{equation*}
\mathrm{df}\left(\mathbf{X}\hat{\beta}_{\mathrm{sc}}\right)=E_{Y|X}\left[\mathrm{rank}(\tilde{\mathbf{X}}_{\mathcal{A}}) \right]-1,
\end{equation*}
where $\mathcal{A}=\mathcal{A}\left(\mathbf{Y}\right):=\{ j : \hat{\beta}_{\mathrm{sc},j}(\mathbf{Y}) > 0 \}$ is the active support corresponding to a solution $\hat{\beta}_{\mathrm{sc}}(\mathbf{Y})$ of the synthetic control problem (\ref{eq:SCwithoutCOVbeginning})-(\ref{SCwithoutCOVend}).
\label{proposition4}
\end{proposition}
With continuous data, the result is even more directly interpretable
as it can be formulated directly in terms of the number of nonzero
weight coefficients.\footnote{The rank deficient case can still be handled by selecting a ``canonical''
active set $\mathcal{A}^{*}$, for instance that selected by the penalized
estimator of \citeA{abadie2021penalized} for an arbitrarily small but
nonzero tuning parameter for the penalty term.}
Let $\mathcal{H}\left(\mathbf{X}\right)
$ designate the convex hull of the columns of the design matrix.\footnote{Specifically, $\mathcal{H}\left(\mathbf{X}\right)
=\left\{ \sum_{j=1}^p \beta_j\mathbf{X}_j: \beta_j\geq 0 \ \text{for all} \ j, \ \sum_{j=1}^p \beta_j=1 \right\}$.}
\begin{corollary}
Suppose that the conditions of \Cref{proposition4} hold and that the conditioned upon design matrix $\mathbf{X}$ follows a distribution that is absolutely continuous with respect to Lebesgue measure on $\mathbb{R}^{n\times p}$. If
$\mathbf{Y}\notin \mathcal{H}\left(\mathbf{X}\right)$,
then the synthetic control solution is unique. Moreover, with probability one over $\mathbf{X}$, we have
\begin{equation*}
\mathrm{df}\left(\mathbf{X}\hat{\beta}_{\mathrm{sc}}\right)=E_{Y|X}\left|\mathcal{A}\right|-1.
\end{equation*}
\label{corollary:dof}
\end{corollary}
\vspace{-2.5em}
An immediate implication of \Cref{corollary:dof} is that one less than the number of non-zero coefficients is an unbiased estimate of the degrees of freedom of the synthetic control method without covariates.
The finite-sample unbiasedness is an attractive feature of the estimate, as emphasized in \citeA{shen2002adaptive,efron2004estimation,shen2006optimal}.
Furthermore, the estimate is consistent, as is the analogous degrees of freedom estimate for the lasso
\citeA{zou2007degrees}. See \Cref{subsectionA5}.
\Cref{figure2} illustrates the theory on simulated data.\footnote{The simulation was done with $p=40$ donors, $n=200$ observations
for each data set, and each point is produced by averaging over 400
simulations. In order to get more variation in the degrees of freedom
of the synthetic control method, the equality constraint $\mathbf{1}_{p}^{T}\beta=1$
was replaced with $\mathbf{1}_{p}^{T}\beta=a$ for different values
of $a$.} We compare with the famous lasso result (\citeNP{zou2007degrees}), which
states that $\mathrm{df}(\mathbf{X}\hat{\beta}_{\mathrm{Lasso}})=E_{Y|X}\left|\mathcal{A}_{\mathrm{Lasso}}\right|$,
where $\mathcal{A}_{\mathrm{Lasso}}$ is the index set of the active independent variables in the lasso regression.
\begin{figure}
\begin{center}
\includegraphics[scale=0.45]{graphs/theoretical_simulation/plot_df_nocov_big.pdf}
\includegraphics[scale=0.45]{graphs/theoretical_simulation/df_best_subset.pdf}
\end{center}
\caption{\emph{Degrees of freedom of the synthetic control method without covariates
(blue) and of the lasso (green) on the left-hand side, and of the
best subset selection regression on the right-hand side.}}
\label{figure2}
\end{figure}
The right-hand side plot of \Cref{figure2} illustrates the ``real cost'' of unrestricted model selection in terms of degrees of freedom and potential overfitting. The degrees of freedom of the best subset selection method,
with respect to the size of the subset, are presented for comparison
with their synthetic control analog. While a synthetic control fitted
model has degrees of freedom one smaller than the expected number
of active donors, we can see that the explicit model search introduces
much more model flexibility.
Our main interest lies in penalized synthetic control methods.
We are therefore interested in their degrees of
freedom.
\begin{proposition}[Degrees of Freedom of Penalized Synthetic Controls]
Suppose that $\mathbf{Y}|\mathbf{X}$ follows the probability law stipulated in (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}).
Suppose that the conditioned upon design matrix $\mathbf{X}$ follows a distribution that is absolutely continuous with respect to Lebesgue measure on $\mathbb{R}^{n\times p}$.
Then, with probability one over $\mathbf{X}$,
\begin{equation*}
\mathrm{df}\left(\mathbf{X}\hat{\beta}_{\mathrm{pen}}\right)=(1+\lambda)(E_{Y|X}\left[\mathrm{rank}(\mathbf{X}_{\mathcal{A}}) \right]-1),
\end{equation*}
where $\mathrm{rank}(\mathbf{X}_{\mathcal{A}})=\min\{|\mathcal{A}|,n\}$ and $\mathcal{A}=\mathcal{A}\left(\mathbf{Y}\right):=\{ j : \hat{\beta}_{\mathrm{pen},j}(\mathbf{Y}) > 0 \}$ is the active support corresponding to the solution $\hat{\beta}_{\mathrm{pen}}\left(\mathbf{Y}\right)$ of the penalized synthetic control problem (\ref{eq:PenalizedCSbeginning})--(\ref{eq:PenalizedSCend}).
\label{proposition5}
\end{proposition}
When interpreting \Cref{proposition5}, it is important to note that $E|\mathcal{A}|$ is a --typically decreasing-- function of $\lambda$.
It may also be instructive to observe that, if one takes $\lambda$ to be arbitrarily large, one obtains (abstracting away from ties) the one-nearest neighbor matching estimator with $\ell_2$ distance.
\begin{proposition}[Degrees of Freedom of Constrained Ridge SCM]
Suppose $\mathbf{Y}|\mathbf{X}$ follows the probability law stipulated in (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}).
Then, the constrained ridge synthetic control fit $\hat{\mathbf{Y}}_{\mathrm{crsc}}=\mathbf{X}\hat{\beta}_{\mathrm{crsc}}+\hat{\beta}_{\mathrm{crsc},0}\mathbf{1}_n$ has degrees of freedom
\begin{equation}
\mathrm{df}\left(\hat{\mathbf{Y}}_{\mathrm{crsc}}\right)=\operatorname{E}_{Y|X}\left[\sum_{j=1}^{|\mathcal{A}|}
\frac{s_j^2}{s_j^2+\lambda}+\lambda
\cfrac{\mathbf{1}_{|\mathcal{A}|}^T\left( \mathbf{V}\left( \mathbf{S}^2+\lambda\mathbf{I}_{|\mathcal{A}|} \right)^2\mathbf{V}^T \right)^{-1}\mathbf{1}_{|\mathcal{A}|}}{\mathbf{1}_{|\mathcal{A}|}^T\left( \mathbf{V}\left( \mathbf{S}^2+\lambda\mathbf{I}_{|\mathcal{A}|} \right)\mathbf{V}^T \right)^{-1}\mathbf{1}_{|\mathcal{A}|}}\right],
\label{dfcrsc}
\end{equation}
where $\mathcal{A}=\mathcal{A}\left(\mathbf{Y}\right):=\{ j: \hat{\beta}_{\mathrm{crsc},j}(\mathbf{Y}) > 0 \}$ is the active support, the singular value decomposition of $\tilde{\mathbf{X}}_{\mathcal{A}}=(\mathbf{I}_n-\frac{1}{n}\mathbf{1}_n\mathbf{1}_n^T)\mathbf{X}_{\mathcal{A}}\in\mathbb{R}^{n\times|\mathcal{A}|}$ has the form $\tilde{\mathbf{X}}_{\mathcal{A}}=\mathbf{U}\mathbf{S}\mathbf{V}^T$, $\mathbf{U}\in\mathbb{R}^{n\times|\mathcal{A}|}$ and $ \mathbf{V}\in\mathbb{R}^{|\mathcal{A}|\times|\mathcal{A}|}$ are orthogonal matrices, and $\mathbf{S}\in\mathbb{R}^{|\mathcal{A}|\times|\mathcal{A}|}$ is a diagonal matrix with diagonal entries $s_1\geq s_2\geq...\geq s_{|\mathcal{A}|}\geq 0$, the singular values of $\tilde{\mathbf{X}}_{\mathcal{A}}$.
\label{prodfcrsc}
\end{proposition}
\begin{remark}
The change-of-variable from $\mathbf{X}_{\mathcal{A}}$ to $\tilde{\mathbf{X}}_{\mathcal{A}}$ aligns with the observation that demeaning the data before applying the SCM estimator is equivalent to adding an intercept (Footnote 4, \citeNP{ferman2021synthetic}).
\end{remark}
\begin{remark}
We discuss two special cases. As $\lambda\downarrow0$, the constrained ridge SCM approaches the synthetic control estimator with an intercept. Its degrees of freedom is the expected number of nonzero coefficients $E|\mathcal{A}|$, implying that adding the intercept costs one degree of freedom. If we remove the sum-to-one and nonnegative constraints, then the estimator (after demeaning) reduces to the usual ridge regression, whose degrees of freedom is exactly the first term of (\ref{dfcrsc}).
We thus recuperate, as a special case, the result of \citeA[~Equ. 3.50]{hastie2009elements}.
\end{remark}
The elastic net SCM can be expressed as a lasso problem and is treated as such by \citeA{tibshirani2012degrees}.
\begin{proposition}[Degrees of Freedom of Elastic Net SCM, \cite{tibshirani2012degrees}]
Suppose that $\mathbf{Y}|\mathbf{X}$ follows the probability law stipulated in (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}).
Then, the elastic net synthetic control fit $\hat{\mathbf{Y}}_{\mathrm{elast}}=\mathbf{X}\hat{\beta}_{\mathrm{elast}}+\hat{\beta}_{\mathrm{elast},0}\mathbf{1}_n$ has degrees of freedom
\begin{align*}
\mathrm{df}(\hat{\mathbf{Y}}_{\mathrm{elast}})&=E_{Y|X}\left[\mathrm{Tr}\left(\tilde{\mathbf{X}}_{\mathcal{A}}\left(\tilde{\mathbf{X}}_{\mathcal{A}}^{T}\tilde{\mathbf{X}}_{\mathcal{A}}+\lambda_{2}\mathbf{I}_{|\mathcal{A}|}\right)^{-1}\tilde{\mathbf{X}}_{\mathcal{A}}^{T}\right)\right]+1
=E_{Y|X}\left[\sum_{j=1}^{|\mathcal{A}|}
\frac{s_j^2}{s_j^2+\lambda_2}\right]+1,
\end{align*}
where
$\mathcal{A}=\mathcal{A}\left(\mathbf{Y}\right):=\{ j: \hat{\beta}_{\mathrm{elast},j}(\mathbf{Y}) \neq 0 \}$ is the active support and $\tilde{\mathbf{X}}_{\mathcal{A}}$ as well as $\mathbf{S}=\text{diag}(s_1,...,s_{|\mathcal{A}|})$ are defined in \Cref{prodfcrsc}.
\label{propescm}
\end{proposition}
\begin{remark}
The degrees of freedom for elastic net SCM and ridge regression with an intercept share the same expression. However, their computed values generally differ because the active set $\mathcal{A}$ depends on the estimated coefficients, which vary by method.
\end{remark}
In the synthetic control problems with covariates (\ref{eq:SCwithCOVbeginning})-(\ref{eq:SCwithCOVend}), the ``covariates'' play the role of
special observations, or time periods, and the synthetic control coefficient $\beta$ is constrained
to produce the best fit possible on these special observations.
When the inner problem admits multiple solutions, it restrains the model flexibility
without imposing a unique solution, and that is reflected by a reduction in the degrees of
freedom equal to the number of covariates, $n_{\mathrm{cov}}$.
Define the augmented donor matrix $\tilde{\mathbf{X}}=(\mathbf{X}^\top,\mathbf{D}^\top,\mathbf{1}_p)^{\top}\in\mathbb{R}^{(n+n_{\mathrm{cov}}+1)\times p}$
and let $\mathrm{relint}\,\mathcal{H}(\mathbf{D})$ designate the relative (to its affine hull) interior of $\mathcal{H}(\mathbf{D})$.
\begin{proposition} [Degrees of Freedom of SCM With Covariates]
Suppose that, conditionally on $(\mathbf{X},\mathbf{D},\mathbf{Z})$, the outcome $\mathbf{Y}$ follows the probability law stipulated in (\ref{GaussianAssumption}) and (\ref{eq:Ahomoskedastic}).
Suppose that $\mathbf{D}\in\mathbb{R}^{n_{\mathrm{cov}}\times p}$ follows a distribution that is absolutely continuous with respect to Lebesgue measure and that $n_{\mathrm{cov}} \le p$.
If $\mathbf{Z}\in \mathrm{relint}\, \mathcal{H}(\mathbf{D})$, then
\[
\mathrm{df}(\hat{\mathbf{Y}})
=
E_{Y\mid X,D,Z}\!\left[\mathrm{rank}\!\left(\tilde{\mathbf{X}}_{\mathcal{A}}\right)\right]
-
n_{\mathrm{cov}}-1,
\]
with probability one over $\mathbf{D}$, where $\hat{\mathbf{Y}}=\mathbf{X\hat{\beta}(\mathbf{Y})}$, and $\mathcal{A}(\mathbf{Y})=\{j:\hat{\beta}_j(\mathbf{Y})>0\}$ is the active set corresponding to a solution $\hat{\beta}(\mathbf{Y})$ of the synthetic control problem
(\ref{eq:SCwithCOVbeginning})--(\ref{eq:SCwithCOVend}).
\label{proposition7}
\end{proposition}
Analogously to the case without covariates, if the design matrix $\mathbf{X}$ is drawn from a continuous distribution, the degrees of freedom expression can be further simplified.
\begin{corollary}
Assume the conditions of \Cref{proposition7} hold and $(\mathbf{X},\mathbf{D})$ has a distribution that is absolutely continuous with respect to Lebesgue measure. If $\mathbf{Y}\notin \{ \mathbf{X}\beta: \mathbf{D}\beta=\mathbf{Z}, \mathbf{1}^T\beta=1, \beta \geq 0 \}$,
then the synthetic control solution is unique. Moreover, with probability one over $(\mathbf{X},\mathbf{D})$, then
\begin{equation*}
\mathrm{df}\left(\hat{\mathbf{Y}}\right)=E_{Y|X,D,Z}\left|\mathcal{A}\right|-n_{\mathrm{cov}}-1.
\end{equation*}
\label{cor:card_with_cov}
\end{corollary}
\vspace{-2.5em}
In some cases, the covariates uniquely determine the synthetic control coefficient $\beta$.
This happens when the independent variable covariate $\mathbf{Z}$ cannot be reproduced exactly
as a convex combination of the dependent covariate variables $\mathbf{D}$.
In other words, when $\mathbf{Z}$ is not in the convex hull of $\mathbf{D}$, which we denote
$\mathbf{Z} \notin \mathcal{H}(\mathbf{D})$.
When that is the case, the outcome variable $\mathbf{Y}$ does not influence the fit $\hat{\mathbf{Y}}$,
and no overfitting whatsoever arises.
With cross-validation for instance, when that is the case, the in-sample fit in the training
set is of the same magnitude---and equal in expectation---to the out-of-sample fit in the
test set.
This absence of fitting to $\mathbf{Y}$ is likewise captured by the degrees of freedom of the
synthetic control method with covariates, which are zero in that specific case.
\begin{proposition} [Degrees of Freedom of SCM With Covariates]
Suppose all assumptions stated in \Cref{proposition7} hold. If $\mathbf{Z}\notin \mathcal{H}(\mathbf{D})$, then the synthetic
control fitted values $\hat{\mathbf{Y}}$ have degrees of freedom $\mathrm{df}(\hat{\mathbf{Y}})=0$.
\label{prop:df_SCM_cov_notin}
\end{proposition}
\subsection{Information Criteria for the Synthetic Control Methods}\label{sec:info-criteria}
The information criteria estimate, the sample analog to \eqref{eq:homoskedasticSURE},
is
\begin{equation}
\widehat{\mathrm{IC}}:=\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2\hat{\sigma}^{2}\hat{\mathrm{df}}\left(\hat{\mathbf{Y}}\right).
\label{eq:information-criteria-estimate}
\end{equation}
All the degrees of freedom expressions $\mathrm{df}(\hat{\mathbf{Y}})$
derived in \Cref{sec:df-sc} are expectations of observed quantities, the
latter are thus used as natural sample analogs $\hat{\mathrm{df}}(\hat{\mathbf{Y}})$.
The sample variance in (\ref{eq:information-criteria-estimate}) is estimated as
\[
\hat{\sigma}^{2}=\frac{1}{n-\hat{p}}\sum_{i=1}^{n}\hat{\varepsilon}_{i}^{2},
\]
where the fitted residuals are from the unpenalized synthetic control estimate
and $\hat{p} = \hat{\mathrm{df}}(\mathbf{X}\hat{\beta}_{\mathrm{sc}})$.
This is the general approach suggested by \citeA{friedman2001elements}.
If one is concerned that the unpenalized synthetic control method is overfitting, our recommended, conservative approach is to use out-of-sample errors as $\hat{\varepsilon}_i$, for $i$ ranging over a test set in the pre-treatment period.
This will naturally tend to produce an inflated, rather than a deflated, variance estimate.
Because larger values of the tuning parameter in, say, the penalized synthetic control method, effectively shrink the estimate towards the matching estimate and tend to produce an estimate that is more stable, we consider this ``upward bias'' as erring on the side of caution.
See Section \ref{section4} for an example.
The estimated information criterion used to select the tuning parameter
$\lambda$ when implementing the penalized synthetic control method
is
\[
\widehat{\mathrm{IC}}_{\mathrm{pen}}(\lambda):=\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2\hat{\sigma}^{2}\left(1+\lambda\right)\left(\left|\mathcal{A}\right|-1\right),
\]
where, importantly, both $\hat{\mathbf{Y}}$ and $\mathcal{A}$ are functions of $\lambda$.
Likewise, the tuning parameter $\lambda$ in the constrained ridge SCM, the pair $(\lambda_1, \lambda_2)$ in the elastic net SCM, and the weighting matrix $V$ in the SCM with covariates can be selected by minimizing the sample information criteria that append the in-sample loss $\|\mathbf{Y} - \hat{\mathbf{Y}}\|_2^2$ with the corresponding sample analog of the degrees of freedom, as provided in \Cref{prodfcrsc,propescm,proposition7}, respectively.
\subsection{Heteroskedasticity-Robust Information Criteria}\label{sec:discussion-Hete}
The homoskedasticity assumption \eqref{eq:Ahomoskedastic} may appear restrictive to some users and this concern warrants discussion.
We find that the proposed information criteria (\ref{eq:homoskedasticSURE}) performs well in the regime of moderate heteroskedasticity but exhibits substantial bias in simulation settings with high heteroskedasticity, see the left panel in \Cref{fig:robustsim}.
We are thus motivated to produce an alternative information criterion that is more robust to heteroskedasticity.
As detailed in Appendix \ref{subsection:consistency_hc}, one such information criterion for high-heteroskedasticity regimes is
\begin{equation}
\widehat{\mathrm{IC}}_{\mathrm{HR}}=n^{-1}\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+\frac{2}{n-\hat{p}}\sum_{i=1}^{n}\hat{\varepsilon}_{i}^{2}\frac{\partial\hat{Y}_{i}}{\partial Y_{i}}
\label{eq:robustIC}
\end{equation}
where $\hat{p}:=\hat{\mathrm{df}}(\mathbf{X}\hat{\beta}_{\mathrm{sc}})$ is the degrees of freedom estimate of the unpenalized SCM estimator. The second term in (\ref{eq:robustIC}) serves as a consistent estimator for the covariance penalty term in (\ref{genSURE}).
The central and right panels of Figure \ref{fig:robustsim} display the performance of the heteroskedasticity robust information criteria estimate (\ref{eq:robustIC}) under both homoskedasticity and high heteroskedasticity regimes. We see that the heteroskedasticity robust information criteria remains unbiased in the homoskedastic case and approximates well, on average, its population analog even under severe heteroskedasticity.
Our information criteria can also be extended to accommodate a non-diagonal conditional covariance matrix $\Sigma_{Y|X}$, as arises in the presence of serial dependence.
By leveraging tools from heteroskedasticity-and-autocorrelation-robust (HAR) variance estimation, we produce a HAR information criteria that is robust to both conditional heteroskedasticity and serial correlation of unknown form. The formal procedure is detailed in Appendix \ref{section:hac_ic}.
The methodology can also accommodate structured conditional covariance matrices.
Consider \citeA{liu1994siegel}, who gives a multivariate
generalization of Stein's Lemma according to which the information
criteria become
\[
E\left\Vert \mathbf{Y}-\hat{\mathbf{Y}}\right\Vert _{2}^{2}+2E\left[\mathrm{Tr}\left(\Sigma_{Y|X}E\left[\left.\nabla\hat{\mathbf{Y}}\right|\mathbf{X}\right]\right)\right].
\]
The form is more general, but it brings about the challenge of estimating
$\Sigma_{Y|X}$. It may be estimated in different ways.
If one is
willing to make the necessary model specification assumptions,
one may
use a structured conditional covariance model.
\begin{figure}
\begin{center}
\includegraphics[scale=0.25]{graphs/theoretical_simulation/HC_heterDGP_homoIC.pdf} \
\includegraphics[scale=0.25]{graphs/theoretical_simulation/HC_heterDGP_heterIC.pdf} \
\includegraphics[scale=0.25]{graphs/theoretical_simulation/HC_homoDGP_heterIC.pdf}
\end{center}
\vspace{1em}
\footnotesize\textit{} The left-hand side figure presents pairs of population and average estimated IC using (\ref{eq:homoskedasticSURE}) for different samples simulated from a high-heteroskedasticity data generating process. The center and right-hand side figures present pairs of population and average estimated IC using (\ref{eq:robustIC}) for different samples simulated from a high-heteroskedasticity and homoskedastic data generating process, respectively. In the heteroskedastic DGP, the outcome variables are generated using a linear model with covariate $X_{t,i}$ and errors drawn from a normal distribution with mean zero and variance proportional to $\exp{(X_{t,1})}$.
\caption{\emph{Assessment of robustness to heteroskedasticity}}
\label{fig:robustsim}
\end{figure}
\subsection{Discussion of Gaussian Assumption}\label{sec:discussion-Gaussian}
While developing an estimator under Gaussian assumptions is by no
means exceptional, it remains important to ask if such an assumption
can be relaxed, and if we are robust to its misspecification.
First,
recent developments have formalized the intuition that degrees of freedom as herein defined capture the flexibility of a method in a robust manner. Specifically, \citeA{fathi2022relaxing} provide theory formalizing how to assess the bias incurred when using standard forms of SURE even though the observations are not Gaussian.
They argue that the bias can be sufficiently small so as to yield estimates useful for the selection of tuning parameters. Formally, they produce bounds on the bias of the SURE estimate which go to zero as the distribution of the observations approaches a normal distribution.
Second, we can see from the simulations in \Cref{sec:simulation} that while the
semi-synthetic data is decidedly not Gaussian, our estimate of the
degrees of freedom, and thus of out-of-sample performance, is remarkably
robust.
\section{Forecasting Counterfactual Car Sales Under Rationing in Tianjin}
\label{section4}
In Tianjin, on December 16, 2013, the municipal government introduced
a hybrid half-lottery and half-auction system for the procurement
of license plates which heavily rationed the number of licenses
issued.\footnote{The measures to control vehicle purchase and restrict the traffic
in Tianjin were first announced in a press conference by the Tianjin
Municipal People\textquoteright s Government at 7 p.m. of December
15, 2013 (The State Council of the People\textquoteright s Republic
of China, 2013). The controls and restrictions were effective five
hours later on December 16, 2013 at midnight. The Tianjin Municipal
People\textquoteright s Government suspended any new vehicle registration
and vehicle transfer in Tianjin between December 16, 2013 and January
15, 2014 in order to ensure a smooth transition and preparation for
the new rules. For more detailed background information, see \citeA{daljord2021black}.} This was done in order to limit pollution from car emissions, which
had become a public health hazard.
Because the auction allowed wealthier individuals a better chance
of obtaining a license, the rationing induced a change in the population
of car buyers, and therefore a change in the demand for different
models (\citeNP{li2018better}).
This change in demand may be non-trivial. While \citeA{daljord2021black}
document the change in the population of car buyers in Beijing when there is a lottery without auction and the licenses are transacted solely
on a black market, they do not study the impact of rationing
on the demand for specific models. However, manufacturers and policymakers
would naturally be interested in the impact of such a policy on the sales
of car models of different prices and fuel consumption.
Specifically, we are interested in establishing model-specific counterfactual
demand, and thus the impact on sales, of individual car models.
The strategy is to use the city of Shijiazhuang to build counterfactuals.
Shijiazhuang is comparable to Tianjin in that, for instance, it is
geographically close and has a population on the same order of magnitude
(15 vs 11 million). Crucially, Shijiazhuang did not have rationing
and we observe sales data for that city.
We highlight two features of the empirical analysis. The first feature is that we use the synthetic control method not because we are missing a good match, but to combine many good albeit noisy matches so as to attenuate variance while minimizing cost in bias. Bias and variance are traded off according to the tuning parameter of the penalized synthetic control.
The second feature is that we consider the simulation and data analysis in conjunction,
and simulate from a data generating process meant to emulate that
of the observed data. We do this for two reasons. First, the simulation
should be thought of as part of the application; to best carry out model
selection (i.e., picking the tuning parameter), we must first elect a model selection method, and since we do not have enough observations to empirically compare model selection methods on observed data (a very data-hungry procedure), we resort to simulated data. Second, we are interested in the simulation for its own sake; in
order to properly assess how practical the theory is, we want to carry out robustness checks on simulated yet realistic data.
\subsection{Empirical Approach}
An apparently natural approach would be to use matching with a single
match. Indeed, all models investigated are sold in both Shijiazhuang
and Tianjin, and thus the Shijiazhuang sales for, say, a Toyota Highlander,
make a natural counterfactual for the Toyota Highlander sales in Tianjin,
post rationing. However, the relatively small number of sales for any
given model makes the time series relatively noisy and methodology
tackling this issue ought to be considered and compared.
On the one hand, filtering is a natural way to tackle the issue of
noisy donors data. Specifically, we want to average over similar time
series and thus reduce noise at little cost in bias. Intuitively,
the underlying demand patterns for analogous models from different
brands such as, say, a Toyota Highlander and a Honda CR-V, may be very
similar, thus allowing for variance reduction at little cost in bias
when averaging over both series to build a counterfactual series for
the Toyota Highlander. The synthetic control method is specifically aimed at
building a control unit by averaging over multiple possible controls.
We thus consider it as a natural alternative to matching.\footnote{Standard $k$-nearest-neighbor matching ($k>1$) would impose equal weights on potentially distant donors. Kernel-weighted matching could relax this restriction, but in practice requires tuning bandwidth parameters.}
To the best of our knowledge, this is a novel use of the synthetic control method.
The SCM is typically used when an obvious match cannot be found; here,
a natural match can be found, but it is noisy and we prefer to combine
with slightly imperfect matches in order to attenuate the noise.
On the other hand, regularization is a natural way to tackle the issue of a noisy to-be-treated, or outcome variable, particularly when there are many independent variables.
This suggests the use
of a regularized alternative to synthetic controls, such as
penalized synthetic controls (\citeNP{abadie2021penalized}).
Conveniently, penalized methods can themselves provide evidence of overfitting when the unpenalized model is too flexible.
This is explored in Figure \ref{Fig:single_model_control_city}, where we see that regularizing the synthetic control estimator improves its out-of-sample performance (left panel), and we find preliminary evidence that the SURE provides a reasonable estimate of the optimal tuning parameter.
Of course the probability simplex constraint, as well as the covariates if any are used, act as a form of regularization.
We find that the additional regularization induced by the above two methods is beneficial in our application.
Since our methodological premise is that the unpenalized synthetic control method overfits, we use out-of-sample errors to estimate the variance term in the information criteria. Specifically, we fit an unpenalized synthetic control estimator on the first two-thirds of the pre-treatment data, and collect the out-of-sample errors $\hat{\varepsilon}_i$, for $i$ ranging over the last third of the pre-treatment data.
\subsection{Simulation}\label{sec:simulation}
We consider two complementary simulation designs.
First, we consider a Gaussian factor model. In that design, the Gaussian assumption underlying the theory applies, and we can furthermore compute the true risk for comparison.
Second, we keep with the factor model but sample the outcome unit's residuals from their empirical distribution, thus assessing the robustness of the theory and the relative performance of different model selection methods in a more realistic design in which we can still compute the true risk.
\begin{figure}
\begin{center}
\includegraphics[scale=0.30]{graphs/single_model/overfit_LHS.pdf} \ \ \ \ \
\includegraphics[scale=0.30]{graphs/single_model/overfit_RHS.pdf}
\end{center}
\footnotesize\textit{}
Left panel: out-of-sample prediction mean-squared error, computed by first estimating PSCM using pre-treatment Highlander market shares data from Shijiazhuang and then evaluating the mean-squared error in the post-treatment periods. Right panel: estimated SURE computed using the same pre-treatment data. Both are plotted as functions of $\lambda$, with the vertical blue dashed line in each panel indicating the minimand.
\caption{\emph{\label{fig:bootstrap_risk}Prediction MSE and SURE Estimates for Different \(\lambda\) Values in PSCM (Shijiazhuang)}}
\label{Fig:single_model_control_city}
\end{figure}
While our interest is in the quality of fit of the synthetic control
method, we use and estimate a more involved but more flexible model to simulate from. Motivated by the theoretical insight of \citeA{abadie2010synthetic} and by the successful implementation on Current Population
Survey (\citeNP{ferman2021synthetic}), we fit a factor model.
We estimate a factor model on the Shijiazhuang data and simulate from the fitted model. Let $Y_{t}$ denote the sales of the treated unit in month $t$, and let $\mathbf X_{t}=(X_{1t},\dots ,X_{pt})^{T}$ collect sales for the $p$ donor units. We model sales as
\begin{equation}
Y_t=\psi_t^{T}L_0+U_{0t}\quad \mathrm{and}\quad X_{jt}=\psi_t^{T}L_j+U_{jt},
\qquad j=1,\dots,p,\; t=1,\dots,T,
\label{eq:fmy}
\end{equation}
where $\psi_t\in\mathbb{R}^r$ is a vector of common factors, $L_j\in\mathbb{R}^r$ is the factor loading for unit $j$, and $U_{jt}$ is an idiosyncratic error. Let $\mathbf{U}_t=(U_{0t},U_{1t},...,U_{pt})^T$, we normalize that
\begin{equation}
E\left[\psi_{t}\right]=\mathbf{0}_{r},\ E\left[\psi_{t}\psi_{t}^{T}\right]=\mathbf{I}_r,\ E[\mathbf{U}_t]=\mathbf{0}_{p+1}, \ E\left[\mathbf{U}_{t}\mathbf{U}_{t}^{T}\right]=\Sigma,
\label{equ:factor_restriction}
\end{equation}
and assume $E[\psi_t\mathbf{U}_t^T]=0$. We take $\Sigma \in \mathbb{R}^{(p+1)\times (p+1)}$ to be diagonal so the idiosyncratic shocks are uncorrelated across units. To connect the factor structure to synthetic controls, let $L_{-0}=[L_{1},\dots ,L_{p}]\in\mathbb R^{r\times p}$, and impose $L_{0}=L_{-0}\,\beta^{*}$ for some weights $\beta^{*}\geq 0, \mathbf{1}^T\beta^{*}=1$. In implementation, we use $\hat{\beta}$ from the unpenalized synthetic control fit as a plug-in for $\beta^{*}$.
We estimate the factors and loadings from the Shijiazhuang panel using the algorithm of \citeA{xu2017generalized}.\footnote{We model the preprocessed data, which was passed through an MA(3) filter and was demeaned.
We fit separate models for sales and market shares. The extended model with time fixed effects was considered but has higher information criteria \cite{bai2002determining}. }
To generate a full panel, we simulate a total of $T=36$ time periods from the factor model, consisting of $T^*=23$ pre-treatment periods and 13 post-treatment periods, as in the application.
For our first simulation design, we consider the Gaussian factor model. Specifically, we simulate the data independently across $t$ from
\[
\left(\begin{array}{c}
Y_{t}\\
\mathbf{X}_{t}
\end{array}\right)\sim N\left(\mathbf{0}_{p+1},\mathbf{L}\mathbf{L}^{T}+\Sigma\right),
\]
where $\mathbf{L}=\left(L_{0},L_1,...,L_{p}\right)^{T}\in\mathbb{R}^{(p+1)\times r}$. Under this design, $Y_t|\mathbf{X}_t$ is Gaussian, and the conditional expectation $E\left[\left.Y_{t}\right|\mathbf{X}_{t}\right]$
is available in closed form (\ref{equ:expectation_expression}), and the true risk is analytically tractable.
In this simulation design, the assumptions underlying the SURE hold exactly, and the degrees of freedom estimate is unbiased, as illustrated in the top-left panel of Figure \ref{fig:Study-of-Robustness.}.
For the second simulation design, we maintain the Gaussian factor structure for the $p$ donor units, but replace the treated unit's idiosyncratic shocks with draws from the empirical distribution of the fitted residuals. This preserves the common factor structure while inducing non-Gaussianity in $Y_t|\mathbf{X}_t$. Because the donors remain Gaussian and the treated shock remains mean-zero and independent of $\mathbf{X}_t$, the conditional mean $E[Y_t|\mathbf{X}_t]$ remains available in closed-form, and the true risk can still be computed exactly. Appendix \ref{simdetails} provides simulation details.
\begin{figure}
\begin{center}
\includegraphics[scale=0.33]{graphs/empirical_simulation/gaussian_dof.pdf}
\includegraphics[scale=0.33]{graphs/empirical_simulation/empirical_dof.pdf}
\end{center}
\begin{center}
\includegraphics[scale=0.33]{graphs/empirical_simulation/residual_distribution.pdf}
\includegraphics[scale=0.33]{graphs/empirical_simulation/residual_qq.pdf}
\end{center}
\footnotesize\textit{} Estimated
degrees of freedom versus true degrees of freedom for the Gaussian factor
model (top left), Estimated degrees of freedom versus true degrees
of freedom for the Gaussian factor model with empirical residuals (top
right), histogram of empirical residuals (bottom left), quantile-quantile
plot of empirical residuals (bottom right).
\caption{\emph{\label{fig:Study-of-Robustness.}Assessment of robustness to Gaussian assumption}}
\end{figure}
In Figure \ref{fig:Study-of-Robustness.} we study the robustness
of our estimate of degrees of freedom, and hence of the risk, to the
ubiquitous failure of the Gaussian assumption. We see in the
histogram (bottom-left panel) and quantile-quantile plot (bottom-right panel) that the fitted residuals are
decidedly not Gaussian, but the quantile-quantile plot nevertheless
suggests a moderate enough departure from Gaussianity that our estimate
of degrees of freedom may still be accurate enough to be useful. This
is validated by the top right panel of Figure \ref{fig:Study-of-Robustness.},
which shows that the estimated degrees of freedom still line up well with
the truth, in expectation, when we simulate the errors by drawing them from the empirical distribution of the fitted residuals.
While the simulation exercise is intrinsically motivated by the need to assess the robustness of SURE to distributional assumptions, it is likewise motivated by the need to select a model selection method for the application whose data generating process the simulation emulates.
To that end, Figure \ref{fig:Risk vs lambda Gaussian sim} displays typical output
from this simulation. We see in the top-left panel that the simulation
design produces a U-shape for the true risk.
This replicates nicely, albeit somewhat more pronouncedly, the plot of estimated risk presented in Figure \ref{Fig:single_model_control_city}.
\begin{figure}
\begin{center}
\includegraphics[scale=0.30]{graphs/empirical_simulation/onedraw_risk.pdf}
\includegraphics[scale=0.30]{graphs/empirical_simulation/onedraw_sure.pdf}
\end{center}
\begin{center}
\includegraphics[scale=0.30]{graphs/empirical_simulation/onedraw_hcv.pdf}
\includegraphics[scale=0.30]{graphs/empirical_simulation/onedraw_vcv.pdf}
\end{center}
\footnotesize\textbf{} Risk (top-left), estimated SURE (top-right), out-of-sample mean-squared error in horizontal cross-validation (bottom-left) and out-of-sample mean-squared error in vertical cross-validation (bottom-right). All as a function of $\lambda$ for a single, representative realization of the simulation. The dashed lines indicate the $\lambda$ values that minimize each respective criterion.
\caption{\label{fig:Risk vs lambda Gaussian sim}\emph{Risk, SURE, and Cross-Validation MSE Across \(\lambda\) in a Single Simulation Realization}}
\end{figure}
The information criteria based on SURE capture this true pattern, and select approximately the same model.
Horizontal and vertical cross-validations (see Appendix \ref{subsection:tuning_parameter_selection} for formal descriptions) display substantially different
patterns, and select tuning parameters far from the oracle tuning parameter, i.e., the one selected according to the population risk.
\begin{table}[!h]
\begin{center}
\footnotesize
\begin{tabular}{lcccc|cccc}
\toprule
\multicolumn{1}{c}{} & \multicolumn{4}{c|}{Gaussian} & \multicolumn{4}{c}{Empirical} \\
\midrule
& {RMSE} & {RMSE} & {Median} & {Mean} & {RMSE} & {RMSE} & {Median} & {Mean} \\
& {$\hat{\tau}_{1}\times 10^{-2}$} & {$\hat{\tau}_{12}\times 10^{-2}$} & {$\hat{\lambda}$} & {$|\hat{\lambda}-\lambda_{risk}|$} & {$\hat{\tau}_{1}\times 10^{-2}$} & {$\hat{\tau}_{12}\times 10^{-2}$} & {$\hat{\lambda}$} & {$|\hat{\lambda}-\lambda_{risk}|$} \\
\midrule
{Risk} & 10.347 & 36.113 & 0.192 & 0.000 & 10.446 & 36.123 & 0.192 & 0.000 \\
{IC$_{\text{oracle}}$} & 10.789 & 37.498 & 0.133 & 0.909 & 10.764 & 37.512 & 0.150 & 0.882 \\
{$\widehat{\text{IC}}$} & 10.768 & 37.462 & 0.277 & 1.007 & 10.802 & 37.538 & 0.313 & 0.979 \\
{CV-Horizontal} & 10.880 & 37.767 & 0.150 & 1.278 & 10.842 & 37.693 & 0.170 & 1.206 \\
{CV-Vertical} & 10.883 & 38.065 & 0.000 & 0.796 & 10.931 & 38.030 & 0.000 & 0.768 \\
{Rolling Window} & 10.869 & 37.691 & 0.192 & 1.133 & 10.828 & 37.665 & 0.217 & 1.083 \\
\bottomrule
\end{tabular}
\end{center}
\footnotesize\textbf{}
The quantities $\hat{\tau}_{1}$ and $\hat{\tau}_{12}$ are the treatment effect estimates one month and one year after treatment, respectively, using the tuning parameter selected by each procedure. The root mean-squared error (RMSE) is evaluated over 5,000 Monte-Carlo replications.
“Median $\hat\lambda$’’ is the median selected tuning parameter, and \(|\hat\lambda-\lambda_{\text{risk}}|\) is the mean absolute deviation from the oracle penalty that minimizes the population risk. Formal definitions for the selection procedures are provided in Appendix \ref{subsection:tuning_parameter_selection}.
\caption{\emph{ Prediction performance and selected tuning parameters from different selection procedures }}
\label{tab:selection_compare}
\end{table}
In Table \ref{tab:selection_compare}, we consider different methods under the
two aforementioned simulation designs.
The proposed information criteria appear to systematically outperform the vertical and horizontal cross-validation approaches and perform comparably to but marginally better than the rolling-window cross-validation (also detailed in Appendix \ref{subsection:tuning_parameter_selection}), whose better performance amongst cross-validation methods had already been documented
\cite{kellogg2020combining}.
From the simulation exercise, we conclude that the SURE information criterion is more reliable for studying the data at hand and we thus elect it as our model selection method for analyzing the data.
\subsection{Data Analysis}
Having elected the information criteria approach \eqref{eq:information-criteria-estimate} as our model selection method, we may now carry out the regression exercise.
We first investigate in detail the impact of rationing for a single, popular model, the Toyota Highlander.
This procedure can be automated
and allows for the joint analysis of a selection of models, which we carry out subsequently.
In preprocessing, we pass the data through an MA(3) filter. The start date of the time series is January 2012, and the end date is December 2014.
\subsubsection*{Analysis for a single model}
We analyze in isolation the impact of rationing on the demand for
the Toyota Highlander. At a high level, the core task is to build a counterfactual;
the time series of demand for the Highlander had there not been rationing.
As intuited and anticipated above, pure matching approaches deliver
a poor counterfactual for the to-be-treated unit. Whether we match
the Highlander in Tianjin to the Highlander in Shijiazhuang, or to the nearest
model in Shijiazhuang according to the penalty term, so in the $\ell_{2}$ sense,
we get an unconvincing fit even as assessed by the plot of the time
series. See Figure \ref{fig:one model fits}.
\begin{figure}
\begin{center}
\includegraphics[scale=0.30]{graphs/single_model/outcome_sc.pdf} \ \ \ \ \
\includegraphics[scale=0.30]{graphs/single_model/outcome_psc.pdf}
\end{center}
\begin{center}
\includegraphics[scale=0.30]{graphs/single_model/outcome_modelmatch.pdf} \ \ \ \ \
\includegraphics[scale=0.30]{graphs/single_model/outcome_mindis.pdf}
\end{center}
\footnotesize\textbf{} Unpenalized synthetic controls (top-left), penalized synthetic
controls with the tuning parameter selected by SURE
(top-right), matching by car model (bottom-left), and matching by $\ell_2$ distance
(bottom-right). The black line plots the observed outcome series and the pale blue line plots the estimated counterfactual series.
\caption{\label{fig:one model fits}\emph{Observed and Estimated Counterfactual Series for the Highlander in Tianjin.}}
\end{figure}
This motivates the use of a synthetic control, a match that is averaged over multiple donors and is thus expected to have lower variance.\footnote{For a more general discussion of methods blending difference-in-differences and synthetic controls, see \citeA{doudchenko2016balancing}.}
Concerns of model flexibility --there are 76 donors after eliminating
models that were either introduced after the beginning of our time
series or discarded before the end-- and noisy to-be-treated series
motivate the use of penalized synthetic controls (\citeNP{abadie2021penalized}). As is generally the case, two general approaches avail to estimate
the tuning parameter of the penalty term: the \emph{cross-validation}
approach and the\emph{ information criteria} approach.
In the simulation exercise of \Cref{sec:simulation}, we found
that the information criteria approach provided a more accurate estimate
of prediction error and, crucially, produced a model that better predicted
treatment effects, especially at short horizons.
It may be that certain time-invariant variables improve the fit when
used as covariates. Two natural covariates are the stock price of
the model's brand and the average --over time-- of its outcome data.
We consider the penalized synthetic control model with these two covariates,
and select both $\lambda$ and $V$ according to our information criteria.
The selected model put full weight on the covariate constructed as
an average of outcome data, and had a worse information criterion than
the model without covariates. We therefore opt for proceeding
without covariates.
The White test for heteroskedasticity (\citeNP{white1980heteroskedasticity}, \citeNP{breusch1979simple}) produced a
$p$-value of $0.12$.
This suggests a moderate
amount of heteroskedasticity and, in light of our simulation study, leaves us sufficiently confident that the theory will remain by and large reliable.
Remark that the choice of tuning parameter selection method is consequential.
As is well exemplified in the case of the Highlander, two different selection methods can yield substantially different tuning parameters, leading to different
treatment effect estimates. As is immediate from Figure \ref{fig:real data plots 1},
the estimate of prediction mean-squared error is essentially monotone in $\lambda$
according to cross-validation, detecting none of the expected overfitting.
The plot of the estimated risk versus $\lambda$, according
to the information criterion estimate, however presents the U-shape
typical of overfitting scenarios; penalization for smaller values
of $\lambda$ reduces overfitting and decreases the prediction error,
but for too large values of $\lambda$ the regression underfits and
the prediction error increases.
The implied difference in the treatment effect estimate is economically
important. For instance, while the unpenalized synthetic control method (which corresponds to $\lambda=0$) estimates an increase in relative demand of 20\% for the Highlander, the information criteria-based estimate predicts an increase of 36\%.
This difference in output from different model selection methodologies ought to be qualified. The conceptual motivation for penalizing is that we are confident, or have an \emph{a priori} belief, that ``nearby'' donors are ``better'' donors; indeed, that is the motivation for matching in the first place.
A more heavily penalized synthetic control estimate is less prone to use ``far away'' donors whose fluctuations may coincidentally cancel and produce a good in-sample fit without capturing signal.
For instance, the linear combination of a luxury minivan, say the Buick GL8, and an ultra economical hatchback, say the Zotye Z100, may give a better in-sample fit for the Highlander than does the linear combination of the CR-V and Highlander itself but, especially since there are many more donors than training periods, we would suspect that this is an instance of overfitting via model selection.
A larger tuning parameter forces the estimator towards the matching
estimator, and away from perhaps coincidental linear combinations
of donors producing a tighter in-sample fit. In a way that echoes the model selection interpretation with the lasso --which forces the estimated regression
function towards zero or another value considered \emph{a priori}
plausible-- the more penalized estimate is more ``conservative''
in a desirable sense.
When penalizing, to what extent are we relying more on the ``nearby'' donors?
The penalty $\sqrt{\sum_{j=1}^{p}\hat{\beta}_{j}\left\Vert \mathbf{Y}-\mathbf{X}_{j}\right\Vert _{2}^{2}}$,
where $\hat{\beta}_{j}$ is a function of $\lambda$, equals 0.92 when no penalty is applied,
and equals 0.37 for $\lambda_{\mathrm{IC}}^{*}=0.3$. Hence, in the $\ell_{2}$
sense, the relevant donors are ``on average'' two and a half times
farther from the to-be-treated unit in the unpenalized synthetic control estimate.
\begin{figure}
\begin{center}
\includegraphics[scale=0.30]{graphs/single_model/selection_sure.pdf} \ \ \ \ \
\includegraphics[scale=0.30]{graphs/single_model/selection_hcv.pdf}
\end{center}
\begin{center}
\includegraphics[scale=0.30]{graphs/single_model/selection_vcv.pdf} \ \ \ \ \
\includegraphics[scale=0.30]{graphs/single_model/selection_te.pdf}
\end{center}
\footnotesize\textbf{} Estimated information criterion (top-left), horizontal cross-validation loss (top-right), vertical cross-validation loss (bottom-left), estimated one-year treatment effect (bottom-right), each plotted as a function of \(\lambda\).
Vertical lines indicate the tuning parameter selected by different methods: \(\lambda^{\ast}_{\mathrm{VCV}}\) (green), \(\lambda^{\ast}_{\mathrm{IC}}\) (red), and \(\lambda^{\ast}_{\mathrm{HCV}}\) (blue).
\caption{\label{fig:real data plots 1}\emph{Tuning Parameter Selection for PSCM and Estimated Treatment Effects for the Highlander in Tianjin}}
\end{figure}
Our main conclusion from the single model analysis is thus that proportional sales for the Toyota Highlander increased substantially
due to the introduction of rationing.
This is quite interesting. The Highlander is a mid-range car
and it was not \emph{ex ante} obvious that it could be a common choice for auction or lottery winners.
We do not attempt to differentiate
between the two types of winners and leave this for further research.
Note however that such a question may be tackled using the optimal
transport reduced form methodology developed in \citeA{daljord2021black}.
\begin{figure}[!h]
\centering
\includegraphics[scale=0.375]{graphs/multiple_models/histogram_treatment_effects_sales.pdf}
\includegraphics[scale=0.375]{graphs/multiple_models/histogram_treatment_effects_shares.pdf}
\begin{flushleft}
\footnotesize Distribution of one-year treatment effects on the number of units sold (left panel) and
distribution of one-year treatment effects on market shares (right panel).
Kernel densities are superimposed on each histogram.
\end{flushleft}
\vspace{-1em}
\caption{\emph{Distribution of car-model-specific treatment effects one year after rationing}}
\label{fig:distribution_effects}
\end{figure}
\subsubsection*{Joint analysis for multiple models}
We wish to investigate the heterogeneity in treatment effects across
different car models. Did luxury cars indeed fare well under the new quota? What about low-end cars? To accommodate such considerations,
we consider simultaneously the treatment effect estimates of multiple
car models.
\begin{figure}
\begin{center}
\includegraphics[scale=0.37]{graphs/multiple_models/multiple_treatment_effects.pdf}
\end{center}
\caption{\emph{\label{fig:Heterogeneous-treatment-effects}Time Series of Car-Model-Specific Treatment Effect Estimates.}}
\end{figure}
We restrict the analysis to the 78 models that sold at least 2,000 units over the sample window.
We automate the procedure carried out for the Toyota Highlander.
For each model, we estimate the penalized synthetic control model using the tuning parameter selected by the estimated SURE information criterion and recover the full counterfactual path of the model.
We first assess the plausibility of the homoskedasticity assumption underlying our suggested risk approximation (\ref{eq:information-criteria-estimate}).
Applying the White test described above, we fail to reject homoskedasticity at the 5\% level for 47 of the 78 models when the outcome is market share (65 at the 1\% level) and for 41 models when the outcome is the number of units sold (59 at the 1\% level). Although heteroskedasticity cannot be ruled out entirely, it seems too modest to overturn the qualitative conclusions.
With all individual treatment effects in hand, we can produce their
histogram and assess visually the distribution of treatment effects.
Figure~\ref{fig:distribution_effects} plots the empirical distribution of one-year treatment effects across models. The left panel shows effects on sales levels; almost all mass lies below zero.
This is expected since the rationing decreased total sales. The right panel plots effects on market shares.
Here the distribution is skewed; most models experience a slight loss in relative share, but the right tail is longer and thinner, capturing the fact that a group of models gained a larger market share.
Of course this does not indicate which specific cars are being sold
in greater or lesser proportions.
To accomplish that, Figure \ref{fig:Heterogeneous-treatment-effects} displays the treatment effect paths for the 13 models whose cumulative sales exceeded 10,000. All series trend downward, as expected given the sharp contraction in license supply, yet the magnitude of the decline varies markedly: upper-mainstream sedans such as the \emph{Magotan} and \emph{Sagitar} lost far fewer units than budget models such as the \emph{Tengyi C30} or the \emph{Corolla EX}.
Figure \ref{fig:effect-versus-price} relates each model’s one-year market-share treatment effect to its average pre-policy price (MSRP). Higher-priced vehicles experienced larger gains (or smaller losses) in relative share. An OLS fit yields a positive slope of \(1.5\times10^{-6}\),
implying that a ¥100,000 increase in price is associated with a 0.15-percentage-point rise in post-rationing market share.
The empirical analysis thus suggests that mid- to high-priced models fared better under rationing, which is consistent with the fact that rationed plates were allocated —through auction or secondary markets— to higher-income households.
\begin{figure}
\begin{center}
\includegraphics[scale=0.37]{graphs/multiple_models/prices.pdf}
\end{center}
\caption{\emph{\label{fig:effect-versus-price}Treatment effect versus price,
given in proportional sales. The dashed line gives the OLS fit.}}
\end{figure}
This confirms the expected
pattern. More expensive cars tend to see a smaller decrease
in sales. The cars with the most extreme reduction in relative demand
tend to be the cheapest.
\section{Conclusion}
\label{section5}
When used in high-dimensional settings, regression methods are often
appended with a penalty term to avoid overfitting, with the lasso being a popular example.
The synthetic control method is no exception, and penalized extensions have been developed to deal with high-dimensional settings.
Such extensions require model selection, either according to cross-validation or to an information criterion.
While information criteria are commonplace for the lasso, and the default in some packages, such methodology was missing for the synthetic control method.
We have developed novel theory and methodology in order to carry out model selection according to an information criterion.
We have argued that the information criteria approach is more reliable than cross-validation in our application of interest.
The herein developed theory delivers degrees of freedom estimates for the synthetic control method.
These in turn produce reassuring theoretical guarantees that the good in-sample fit of the method in early, seminal applications (e.g., \citeNP{abadie2003economic}, \citeNP{abadie2010synthetic}, \citeNP{abadie2015comparative}) was due to the information content and not the flexibility of the synthetic control model.
In analyzing Chinese car sales data, we find that the synthetic control method can be profitably used for filtering when good but noisy matches are available. In our application, this called for the use of penalized variants of SCM in order to avoid overfitting, and of an information criterion to select the correct tuning parameter. Thus equipped, we were able to produce model specific treatment effect paths for a selection of cars.
\newpage