# Co-Optimization of Design and Fabrication Plans for Carpentry: Supplemental
Material

Haisen Zhao<EMAIL_ADDRESS>University of Washington and Shandong
University ,  Max Willsey<EMAIL_ADDRESS>University of Washington
,  Amy Zhu<EMAIL_ADDRESS>University of Washington ,  Chandrakana
Nandi<EMAIL_ADDRESS>University of Washington ,  Zachary Tatlock
<EMAIL_ADDRESS>University of Washington ,  Justin Solomon
<EMAIL_ADDRESS>Massachusetts Institute of Technology and  Adriana Schulz
<EMAIL_ADDRESS>University of Washington

††journal: TOG

##  1\. Technique details

This section provides additional details of our ICEE algorithm for extracting
a Pareto front where each solution $s$ represents a (design, fabrication)
pair.

###  1.1. Generating Design Variants

Design variants can be generated manually by a user or automatically. In our
paper, we implement an automatic method, which includes three steps: (1)
detecting all the assembly connectors between two neighboring parts; (2)
enumerating all candidate connector variants for each connect; (3) generating
design variants by selecting different connector variation of each connector,
as shown in Figure 4 of the main paper.

###  1.2. Fabrication arrangement generation

Given a specific design variation, we use a fabrication arrangement generation
algorithm to update the BOP E-graph such that the E-graph encodes more
fabrication arrangements. This algorithm uses two heuristic described as
Section 4.3.3.

For an input design variation $d_{i}$ which consists of parts $p_{j},0...n$,
this algorithm first groups a library of stock lumbers by their dimensions,
e.g., lumber $2^{\prime\prime}\times 4^{\prime\prime}\times
24^{\prime\prime}$, lumber $2^{\prime\prime}\times 4^{\prime\prime}\times
48^{\prime\prime}$, and lumber $2^{\prime\prime}\times 4^{\prime\prime}\times
96^{\prime\prime}$ are grouped together. For the parts assigned to this group,
this algorithm starts by packing all the parts on the largest stock lumber of
the current group. The packing process will be terminated when (1) the current
stock lumber is maximally packed or (2) all parts are packed. In the first
case, the packing process can continue by switching to another stock lumber
(also the largest one) until the second termination condition is reached. Such
a layout process is called a full Traversal which traverses all of the parts
in a specific order. Prior work [Wu et al., 2019] has used a similar approach
but used the number of Traversals as the termination criterion which prevents
them to control the amount of fabrication arrangements to generate with this
heuristic-driven method.

###  1.3. Cost Metrics

This section describes how we compute the three costs: material usage
($f_{c}$), cutting precision ($f_{p}$), and fabrication time ($f_{t}$). Our
formulae are updated versions of the ones used by Wu et al. [2019]. Our key
improvements are 1) we include a quantitative evaluation method. Each cost
metric is associated with a meaningful unit: material usage ($f_{c}$) is in
dollars, cutting precision ($f_{p}$) is in inches, fabrication time ($f_{t}$)
is in minutes. 2) we evaluate stock load and unload time to make the result
fabrication time much more reasonable.

#### Material Cost

We compute material cost as:

$\displaystyle f_{c}=\sum_{i=1}^{n}p_{i}.$

where $p_{i}$, the price of the $i^{th}$ piece of stock is estimated based on
costs from standard US vendors [McMASTER-CARR, 2021] as shown in Table 1, and
$n$ is the total number of pieces of lumber used.

Stock | Dimension | Material Cost ($)
---|---|---
$2^{\prime\prime}\times 2^{\prime\prime}$ | $24^{\prime\prime}$ | 3.0
$2^{\prime\prime}\times 2^{\prime\prime}$ | $48^{\prime\prime}$ | 5.5
$2^{\prime\prime}\times 2^{\prime\prime}$ | $96^{\prime\prime}$ | 10.0
$2^{\prime\prime}\times 4^{\prime\prime}$ | $24^{\prime\prime}$ | 3.0
$2^{\prime\prime}\times 4^{\prime\prime}$ | $48^{\prime\prime}$ | 5.5
$2^{\prime\prime}\times 4^{\prime\prime}$ | $96^{\prime\prime}$ | 10.0
$4^{\prime\prime}\times 4^{\prime\prime}$ | $24^{\prime\prime}$ | 7.5
$4^{\prime\prime}\times 4^{\prime\prime}$ | $48^{\prime\prime}$ | 13.75
$4^{\prime\prime}\times 4^{\prime\prime}$ | $96^{\prime\prime}$ | 25.0
$2^{\prime\prime}\times 8^{\prime\prime}$ | $24^{\prime\prime}$ | 7.5
$2^{\prime\prime}\times 8^{\prime\prime}$ | $48^{\prime\prime}$ | 13.75
$2^{\prime\prime}\times 8^{\prime\prime}$ | $96^{\prime\prime}$ | 25.0
$1/2^{\prime\prime}$ | $12^{\prime\prime}\times 20^{\prime\prime}$ | 5.5
$1/2^{\prime\prime}$ | $24^{\prime\prime}\times 20^{\prime\prime}$ | 10.0
$1/2^{\prime\prime}$ | $48^{\prime\prime}\times 36^{\prime\prime}$ | 30.0
3/4 ′′ | $12^{\prime\prime}\times 20^{\prime\prime}$ | 7.0
3/4 ′′ | $24^{\prime\prime}\times 20^{\prime\prime}$ | 12.0
3/4 ′′ | $48^{\prime\prime}\times 36^{\prime\prime}$ | 32.0
Table 1. Prices of stocks

Tool | Full setup (s) | Partial setup (s) | Op Time
---|---|---|---
Lumber | Plywood
Chopsaw | 60 | 60 | 15 | 1 s
Bandsaw | 20 | 90 | N/A | 1 inch/s
Jigsaw | 30 | 60 | N/A | 1 inch/s
Tracksaw | 180 | 180 | 75 | 4.5 inch/s
Drill | 20 | 20 | N/A | 0.1 inch (depth)/s
Table 2. Fabrication times for different tools. All times are converted to
minutes in our results.

Tool | Operation Error
---|---
Chopsaw |  $\nicefrac{{1}}{{64}}^{th}$ of an inch
Bandsaw |  $\nicefrac{{1}}{{16}}^{th}$ of an inch
Jigsaw |  $\nicefrac{{3}}{{16}}^{th}$ of an inch
Tracksaw |  $\nicefrac{{1}}{{32}}^{nd}$ of an inch
Drill |  $\nicefrac{{1}}{{32}}^{nd}$ of an inch
Table 3. Error per cut for each tool.

Stock | F-Load (s) | P-Load (s) | F-Unload (s) | P-Unload (s)
---|---|---|---|---
$2^{\prime\prime}\times 2^{\prime\prime}\times 4^{\prime\prime}$ | 10 | 1 | 5 | 1
$2^{\prime\prime}\times 2^{\prime\prime}\times 48^{\prime\prime}$ | 20 | 2 | 8 | 2
$2^{\prime\prime}\times 2^{\prime\prime}\times 96^{\prime\prime}$ | 40 | 3 | 15 | 2
$2^{\prime\prime}\times 4^{\prime\prime}\times 24^{\prime\prime}$ | 10 | 1 | 5 | 1
$2^{\prime\prime}\times 4^{\prime\prime}\times 48^{\prime\prime}$ | 20 | 2 | 8 | 2
$2^{\prime\prime}\times 4^{\prime\prime}\times 96^{\prime\prime}$ | 40 | 4 | 15 | 2
$4^{\prime\prime}\times 4^{\prime\prime}\times 24^{\prime\prime}$ | 15 | 2 | 5 | 1
$4^{\prime\prime}\times 4^{\prime\prime}\times 48^{\prime\prime}$ | 30 | 4 | 10 | 2
$4^{\prime\prime}\times 4^{\prime\prime}\times 96^{\prime\prime}$ | 60 | 6 | 20 | 3
$2^{\prime\prime}\times 8^{\prime\prime}\times 24^{\prime\prime}$ | 15 | 2 | 5 | 1
$2^{\prime\prime}\times 8^{\prime\prime}\times 48^{\prime\prime}$ | 30 | 4 | 10 | 2
$2^{\prime\prime}\times 8^{\prime\prime}\times 96^{\prime\prime}$ | 60 | 6 | 20 | 3
$1/2^{\prime\prime}\times 12^{\prime\prime}\times 20^{\prime\prime}$ | 30 | 3 | 10 | 2
$1/2^{\prime\prime}\times 24^{\prime\prime}\times 20^{\prime\prime}$ | 50 | 5 | 15 | 2
$1/2^{\prime\prime}\times 48^{\prime\prime}\times 36^{\prime\prime}$ | 100 | 10 | 20 | 2
3/4 ${}^{\prime\prime}\times 12^{\prime\prime}\times 20^{\prime\prime}$ | 30 | 3 | 10 | 2
3/4 ${}^{\prime\prime}\times 24^{\prime\prime}\times 20^{\prime\prime}$ | 50 | 5 | 15 | 2
3/4 ${}^{\prime\prime}\times 48^{\prime\prime}\times 36^{\prime\prime}$ | 100 | 10 | 20 | 2
Table 4. Loading and unloading time of stocks. All times are converted to
minutes in our results.

#### Time

We asked an expert carpenter to assign a fabrication time to each tool (Table
3) and estimate a stock loading and unloading time for each piece of wood
stock (Table 4). They reported the time taken for (1) full setup of the tool,
(2) partial setup when applicable (a partial setup in one where only some
parameters are changed), (3) full stock load (unload) time (setting up a piece
of wood stock to the tool workspace) (4) partial stock load (unload) time
(setting up a piece of wood stock to the tool workspace where stocks are
stacked) and (5) performing a single operation (e.g., cut). Out of all the
used tools, tracksaw has the most elaborate setup process which is indicated
by its long setup times (both full and partial). Bandsaw and jigsaw are set to
different full setup time while cutting lumber and plywood sheet. Only chopsaw
and tracksaw allow partial setups, the remaining tools do not. The operation
time for a chopsaw is constant — all cuts take a second. For all other tools,
the operation time is based on the length cut per second. For example, the
operation time for a tracksaw cut of length $l^{\prime\prime}$ is $l/4.5$
seconds.

We refer to a cut as ”partial” if it requires only a partial setup. For
example, if the $i^{th}$ cut is a partial cut on a chopsaw, the time taken for
this cut would be $15s$ for partial setup and $1s$ for the cutting operation
leading to a total of $16s$ for the cut to be completed. The fabrication time,
$f_{t}$ is therefore computed as:

$\displaystyle f_{t}=\sum_{i=1}^{k}(s_{i}+w_{i}+o_{i})$

where $s_{i}$, $w_{i}$ and $o_{i}$ are the setup time, stock load and unload
time, and operation time for the $i^{th}$ cut respectively, and $k$ is the
total number of cuts. $s_{i}$ and $o_{i}$ are computed based on Table 3.
$w_{i}$ is computed based on Table 4. For a cut with stacking, $w_{i}$ is
measured as the sum of a full load and unload time for the first piece of
stock, and multiple partial load and unload times for the following stocks.

#### Precision

To measure the precision of a fabrication plan, we compute (1) measurement
error, $\epsilon$, and (2) operation error, $p$. We use
$m=1/16^{\prime\prime}$ as the minimum measurement that can be made in any of
the currently supported tools. The measurement error for length, $m^{\prime}$
is computed as:

$\displaystyle\epsilon=min(m^{\prime}\%m,m-(m^{\prime}\%m))$

which, intuitively, is the residual length that cannot be measured by our
tools. For operation error, we again asked an expert to estimate the error per
operation for each tool (Table 3). The error for the $i^{th}$ cut is therefore
$\epsilon_{i}+p_{i}$. $f_{p}$ is then computed as:

$\displaystyle f_{p}=\sum_{i=1}^{k}(\epsilon_{i}+p_{i})$

where $k$ is the total number of cuts. Lower values of $f_{p}$ indicate higher
precision.

####  1.3.1. Mixed-material implementation

By introducing a new material, we must accommodate the cost of using this new
material in our metrics. Based on an expert’s input, we set material price of
metal to be 20 times that of wood. All of the same tools can be used to cut
metal sheets, using a specific blade. Setup time is the same, but the
execution times is set to 10 times slower, and stock loading and loading time
is 5 times slower. Finally, the jigsaw’s precision error is set to 2 times
that of wood. Modulo the cost difference, we apply the same process to
evaluate the fabrication cost.

###  1.4. Pareto front extraction

In Section 4.3.4 of the main paper, we propose two methods to speed up the
Pareto front extraction. This section provides more details of the two
methods. For an atomic e-node which has not been previously optimized, we
simply try a maximum, $P$, different orders of cuts (each cutting order is
evaluated using the cost metrics in Section 1.3), then select the cutting
order with minimum precision error ($f_{p}$) and the one with minimum
fabrication time ($f_{t}$).

To apply the branch and bound technique, we need to define the lower and upper
bound for the precision error ($f_{p}$) and the time cost ($f_{t}$). A term
$\mathcal{T}$’s cutting order can be initialized by taking the cutting order
used in its contained atomic e-node. We first select the cutting order with
minimized precision of each atomic e-node, then evaluate its precision error
defined as the precision upper bound of the term. The time upper bound is set
by taking the cutting order with minimized time of each atomic e-node.  The
precision (time) lower bound is defined as the evaluated precision (time)
cost, ignoring the dependency relationship between cuts belongs to the same
stock.  The dependency relationship indicates the order of a cut with respect
to other cuts. In other words, the time lower bound is computed by the
assumption that each cut is independent of the other.

If the precision (or) lower bound is not dominated by the Pareto front of all
computed solutions $\mathcal{S}$, we run an optimization that uses the upper
bound as a starting point. For the optimization, we randomly flip the order of
some cuts until $t$ iterations. We set $t$ as 20 in our experiments.

##  2\. Results and Discussion

In this section, we provide additional results of our ICEE pipeline for
exploring the space of design variations and fabrication plans of Figure 7.
Table 5 shows detailed comparison between our pipeline and the no design
exploration pipeline. At the end of this supplemental material, for each
model, we demonstrate some representative design variations and fabrication
plans extracted with our ICEE pipeline and the baseline method.

Model | $HV_{base}$ | $HV_{our}$ | Ref. Point | $M_{base}(\$)$ | $P_{base}(^{\prime\prime})$ | $T_{base}(m)$ | $M_{our}(\$)$ | $P_{our}(^{\prime\prime})$ | $T_{our}(m)$
---|---|---|---|---|---|---|---|---|---
Frame | 881570 | 902017 | (100.0,100.0) | 10.00 | N/A | 1.98 | 8.50 | N/A | 1.37
L-Frame | 796637 | 797655 | (100.0,100.0) | 18.50 | N/A | 2.15 | 18.50 | N/A | 2.05
A-bookcase | 698414 | 702256 | (100.0,100.0,100.0) | 23.00 | 0.16 | 8.82 | 23.00 | 0.13 | 8.53
S-Chair | 683887 | 698016 | (100.0,100.0) | 25.50 | N/A | 7.92 | 25.50 | N/A | 5.92
Table | 833919 | 849705 | (100.0,100.0) | 20.00 | N/A | 7.18 | 20.00 | N/A | 5.78
F-Cube | 817648 | 825935 | (100.0,100.0) | 15.50 | N/A | 3.08 | 15.50 | N/A | 2.12
Window | 754863 | 777711 | (100.0,100.0) | 20.00 | N/A | 4.93 | 20.00 | N/A | 2.38
Bench | 355824 | 369260 | (100.0,100.0) | 55.50 | N/A | 17.03 | 53.00 | N/A | 20.38
A-Chair | 596346 | 578014 | (100.0,100.0) | 35.50 | N/A | 6.25 | 35.50 | N/A | 9.52
F-Pot | 24155184 | 24264553 | (300.0,300.0,300.0) | 13.00 | 0.29 | 19.13 | 13.00 | 0.26 | 17.91
Z-Table | 19716012 | 20386788 | (300.0,300.0,300.0) | 55.50 | 0.35 | 30.69 | 53.00 | 0.28 | 24.28
Loom | 20469248 | 21040283 | (300.0,300.0,300.0) | 27.50 | 1.47 | 48.26 | 25.50 | 0.58 | 43.86
J-Gym | 13925991 | 15657249 | (300.0,300.0,300.0) | 96.00 | 1.82 | 72.02 | 89.00 | 0.62 | 51.38
D-Chair | 543771 | 539241 | (100.0,100.0) | 35.50 | N/A | 15.00 | 35.50 | N/A | 16.20
Bookcase | 20266354 | 21957242 | (300.0,300.0,300.0) | 40.00 | 0.86 | 37.53 | 30.00 | 0.27 | 24.71
Dresser | 21237296 | 22866170 | (300.0,300.0,300.0) | 30.00 | 0.61 | 36.48 | 30.00 | 0.14 | 15.42
Table 5.  Performance comparison between our method and the baseline method
(no design exploration pipeline). In this table, we first report the
hypervolume value of our method ($HV_{our}$) and the baseline method
($HV_{base}$). The reference points are also listed used for the hyper-volume
computation. We also report the minimal material usage, the minimal cutting
precision and the minimal fabrication time of the Pareto fronts of baseline
method ($M_{base},P_{base},T_{base}$) and the Pareto fronts of our pipeline
($M_{our},P_{our},T_{our}$). The cutting precision of some models is not
reported for that precision cost has not been taken into the objective
metrics. We use bold font to indicate the hypervolume and these fabrication
costs where we are better than the baseline method.

###  2.1. Carpentry Compiler Parameters

During the comparison between our pipeline and the Carpentry Compiler pipeline
[Wu et al., 2019], we use the default parameter setting used in their
experiments, the number of Traversal $T=50$, the number of top e-nodes $n=10$,
the maximal different orders of cuts $P=25$, the population size during
E-graph extraction as $120$, the probability of crossover and mutation
$p_{c}=0.95,p_{m}=0.1$.

###  2.2. Running Time Analysis

On average, our pipeline spends 4.2% of its runtime in design extraction, 5.8%
of its time generating packings, 64.5% of its time assigning cutting orders,
and 21.6% of its time in applying genetic algorithms to extract Pareto front
of terms. This translates to 10.0% of its runtime in the expansion and
contraction phase, 86.1% in the extraction phase, and the other 3.9%
everywhere else.

###  2.3. Comparison with baseline

Figure 1. Pareto fronts generated from our pipeline and baseline method. Model | $HV_{base}$ | $HV_{our}$ | Ref. Point
---|---|---|---
Frame | 901997 | 902017 | (100.0,100.0)
J-Gym | 15036794 | 15657249 | (300.0,300.0,300.0)
L-Frame | 797573 | 797655 | (100.0,100.0)
Table | 751206 | 849705 | (100.0,100.0)
Window | 778158 | 777711 | (100.0,100.0)
Table 6. Hypervolume result of the performance validation experiment. In this
table, we report the hypervolume of the Pareto fronts computed from our method
($HV_{our}$) and the ”baseline” method ($HV_{base}$) . ”Baseline” indicates
extracting the Pareto front fabrication plans fro each design variation
explored by our method independently with the Carpentry Compiler pipeline [Wu
et al., 2019]. The reference points used for hypervolume computation are also
reported.

As shown in Figure 1, we cannot find the exact same front due to convergence
and stochasticity. In order to make this comparison, we have tuned the
parameters to be as close as possible. This means we do not find the exact
same hypervolume indicated in Table 6. However, hypervolume is unintuitive to
compare, so we believe presenting plots for comparison is far more compelling.

Figure 2.  Pareto fronts of the A-Chair model with different alpha, (a) shows
the result of a small alpha ($0.5$), (b) shows the result of the default alpha
($0.75$), (c) shows the result of a bigger alpha ($0.95$). The Pareto fronts
generated from no design exploration pipeline stay the same in (a,b,c). Pareto
fronts of the J-Gym model with tuning parameters of the no design exploration
pipeline, (d) show results of the default parameters, (e) show the Pareto
fronts with tuning parameters of the no design exploration pipeline. The
Pareto fronts generated from our pipeline stay the same in (d,e).

###  2.4. Convergence

Due to the combinatorial nature of the problem, we cannot guarantee that we
discover the true Pareto front. In what follows we discuss how this limitation
can affect the results of our evaluation.

####  2.4.1. Increasing the Design Space

An implication of the intractable search is that we have a finite amount of
resources to spend searching the space. Realistically, this means we can only
explore a subset of the design space. Consider the design of the “Window” in
Figure 7. If we permit this model to have angled joints instead of just 90
degree joints, the space of design variants becomes much larger. Ideally, the
result for this new window variant should be at least as optimal as the result
for the original window. However, our approach was able to explore only a
portion of the vast design space ($|\mathcal{D}|=1520052$, compared to
$|\mathcal{D}|=10463$ for the simpler window model) and found a worse
landscape of solutions (Figure 3). For models with a large design space, if
the initial design variants picked by our algorithm are too far from the
optimal ones, then our algorithm may never reach those optimal points.

Figure 3. Example showing a limitation of our approach: the star dots indicate
the Pareto fronts with allowing the Window model to have arbitrarily angled
connectors; the colored dots indicate ones without arbitrarily angled
connectors. Both results are generated with the default parameter of our ICEE
pipeline.

To avoid similar potential pitfalls, we might be able to allocate our
resources to effectively direct the search. In particular, one thing we can do
to navigate this space more precisely is to expose a parameter, alpha.

$\alpha$ is a parameter that roughly defines the tradeoff between breadth
(exploring many designs) and depth (exploring variations within a small number
of designs). For our experiments, we set $\alpha$ to a default value of
$0.75$. However, in cases where the initial design is believed to be optimal
or close to optimal, the user may increase $\alpha$ to encourage deeper
exploration of the initial and similar designs. In cases where our baseline of
optimizing the fabrication for a single design outperforms exploring many
designs in the default configuration, we have found that increasing $\alpha$
to $0.95$ makes the performance at least as good again.

This note about parameter tuning reveals interesting implications. In
particular, our new algorithm does not always fully dominate fabrication-plan-
only exploration with default parameters. There are models where the input
design is already fairly optimized and doesn’t offer much opportunity for
improvement. Take a model that is simple or is composed of many instances of
the same shape. The Adirondack chair is such an example. It would be clear
that an optimal packing arrangement would minimize material, cuts, and error
by stacking cuts; design variations would only deviate from the simplicity of
the fabrication plan and there are no improvements possible. The baseline
method has the advantage with these models, because it spends its search time
deeply exploring fabrication plans, while ICEE must also spend time searching
across design variations. Tuning the parameter alpha can help us perform as
well as the baseline method, as seen in Figure 2 (a) (b) (c).

From these examples, we see that comparing approaches while using parameter
tuning is trickier and subtler.

####  2.4.2. Comparisons with [Wu et al., 2019]

Notably, comparisons with the approaches not employing design exploration are
also susceptible to the imprecision in convergence. In previous sections, we
ran each tool with its default parameters. Tuning the parameters for each
approach would change the results and make the difference in Pareto fronts
difficult to assess. Neither approach discovers the true Pareto front, so an
absolute comparison depends on the input parameters. Fortunately, we avoid the
brunt of this pitfall: the selling point of our design space exploration
approach is not about eking out marginal improvements over the baseline, but
about exploring a completely different space.

Though we may not search any one design as thoroughly when using default
parameters, our approach is able to do something the baseline cannot: discover
optimizations to the input design when they exist. Consider elaborate models
where there are many ways to pack and fabricate any design, and due to the
complexity and number of parts, the designer has diminishing intuition about
how to achieve some desired point in the space of tradeoffs.

An example that illustrates this complexity is the Jungle Gym model. The model
uses a mix of 1D and 2D wood and has a number of components. From Figure 2 (d)
(e), we see that the initial design was suboptimal, and just by exploring
several different design variations, we find a number of plans that save
greatly on our three objectives. Even the expert plan is completely dominated
by the pareto front, indicating that good designs are not always readily
apparent.

Our tool takes much longer to finish running as it explores the design space
jointly with the fabrication space. Because this is a synthesis problem, one
might ask if allowing the baseline to run for longer would result in better
fabrication plans for the initial design. We chose not to attempt a comparison
of the two algorithms for two reasons. Although it is possible that allowing
the fabrication-space-only algorithm to run longer will help it find lower-
cost plans, there is no parameter that directly controls how long the
algorithm runs. The second reason is that we are supplying a new approach that
explores a different solution space, and targets a different problem.
Controlling for the running time of the algorithm would still not create an
apples-to-apples comparison.

Percentage improvement (%)
---
Model | Carpenter price ($\$/hour$)
0 | 10 | 20 | 40 | 80 | 160 | 240 | 400
Frame | 15 | 19 | 21 | 20 | 15 | 10 | 12 | 16
L-Frame | 0 | 1 | 2 | 2 | 1 | 3 | 1 | 0
A-bookcase | 0 | 1 | 2 | 6 | 7 | 6 | 6 | 5
S-Chair | 0 | 1 | 2 | 2 | 6 | 9 | 11 | 14
Table | 0 | 1 | 2 | 3 | 4 | 7 | 10 | 12
F-Cube | 0 | 4 | 5 | 3 | 3 | 5 | 6 | 8
Window | 0 | 2 | 4 | 12 | 20 | 26 | 28 | 31
Bench | 5 | 4 | 4 | 4 | 4 | 5 | 1 | -4
A-Chair | 0 | -2 | -4 | -4 | -3 | -4 | -9 | -16
F-Pot | 0 | 3 | 4 | 5 | 7 | 9 | 9 | 8
Z-Table | 5 | 5 | 6 | 8 | 11 | 11 | 11 | 12
Loom | 7 | 3 | 1 | 0 | 3 | 3 | 4 | 5
J-Gym | 7 | 7 | 7 | 8 | 11 | 17 | 20 | 23
D-Chair | 0 | 1 | 1 | 2 | 3 | 2 | 0 | -2
Bookcase | 25 | 16 | 10 | 7 | 8 | 12 | 14 | 18
Dresser | 0 | 2 | 4 | 6 | 8 | 25 | 33 | 42
Table 7. The percent improvement of minimum-cost (after scalarization) plans
of the baseline compared to the minimum-cost plans of design space
exploration, when scalarized at different prices.

###  2.5. Scalarization Results

Table 7 shows the scalarization of some of the tradeoffs found on the Pareto
front with our method. Note that these results are all for the wood models
with our default cost metric. When time isn’t worth a lot, plans with low
material cost dominate. When time is worth more much than materials, low time
become cheapest. We observe a variety of fabrication plans being effective at
different points.

Figure 4. Design variants and fabrication plans extracted with our ICEE
algorithm (Frame) Figure 5. Design variants and fabrication plans extracted
with our ICEE algorithm (A-Chair) Figure 6. Design variants and fabrication
plans extracted with our ICEE algorithm (Table) Figure 7. Design variants and
fabrication plans extracted with our ICEE algorithm (Window) Figure 8. Design
variants and fabrication plans extracted with our ICEE algorithm (L-Frame)
Figure 9. Design variants and fabrication plans extracted with our ICEE
algorithm (D-Chair) Figure 10. Design variants and fabrication plans extracted
with our ICEE algorithm (Bench) Figure 11. Design variants and fabrication
plans extracted with our ICEE algorithm (F-Cube) Figure 12. Design variants
and fabrication plans extracted with our ICEE algorithm (J-Gym) Figure 13.
Design variants and fabrication plans extracted with our ICEE algorithm
(Bookcase) Figure 14. Design variants and fabrication plans extracted with our
ICEE algorithm (A-Bookcase) Figure 15. Design variants and fabrication plans
extracted with our ICEE algorithm (Dresser) Figure 16. Design variants and
fabrication plans extracted with our ICEE algorithm (F-Pot) Figure 17. Design
variants and fabrication plans extracted with our ICEE algorithm (Z-Table)
Figure 18. Design variants and fabrication plans extracted with our ICEE
algorithm (Loom)

## References

  * [1]
  * McMASTER-CARR [2021] McMASTER-CARR. 2021\.  McMASTER-CARR.  https://www.mcmaster.com/.
  * Wu et al. [2019] Chenming Wu, Haisen Zhao, Chandrakana Nandi, Jeffrey I Lipton, Zachary Tatlock, and Adriana Schulz. 2019\.  Carpentry compiler.  _ACM Transactions on Graphics (TOG)_ 38, 6 (2019), 1–14.


Synthetic Matching Control Method∗ ††∗We are very grateful to thank for
invaluable comments from Kaspar Wüthrich.  †††Institute of Science and
Technology for Brain-inspired Intelligence, Fudan University, Shanghai, China;
E-mail<EMAIL_ADDRESS>

Rong J.B. Zhu†

Abstract. Estimating weights in the synthetic control method involves an
optimization procedure that simultaneously selects and aligns control units in
order to closely match the treated unit. However, this simultaneous selection
and alignment of control units may lead to a loss of efficiency in the
synthetic control method. Another concern arising from the aforementioned
procedure is its susceptibility to under-fitting due to imperfect pre-
treatment fit. It is not uncommon for the linear combination, using
nonnegative weights, of pre-treatment period outcomes for the control units to
inadequately approximate the pre-treatment outcomes for the treated unit. To
address both of these issues, this paper proposes a simple and effective
method called _Synthetic Matching Control_ (SMC). The SMC method begins by
performing the univariate linear regression to establish a proper match
between the pre-treatment periods of the control units and the treated unit.
Subsequently, a SMC estimator is obtained by synthesizing (taking a weighted
average) the matched controls. To determine the weights in the synthesis
procedure, we propose an approach that utilizes a criterion of unbiased risk
estimator. Theoretically, we show that the synthesis way is asymptotically
optimal in the sense of achieving the lowest possible squared error. Extensive
numerical experiments highlight the advantages of the SMC method.

Keywords: Synthetic Control, Treatment Effects, Panel Data, Ensemble

##  1 Introduction

The synthetic control (SC) method is a popular approach of evaluating the
effects of policy changes. It allows estimation of the impact of a treatment
on a single unit in panel data settings with a modest number of control units
and with many pre-treatment periods (Abadie and Gardeazabal, 2003; Abadie et
al., 2010). The key idea under the SC method is to construct a weighted
average of control units, known as a synthetic control, that matches the
treated unit’s pre-treatment outcomes. The estimated impact is then calculated
as the difference in post-treatment outcomes between the treated unit and the
synthetic control. See Abadie (2021) for recent reviews.

The SC method utilizes a constrained optimization to solve for weights,
typically resulting in sparse weights where only a few control units have non-
zero weights (Abadie and L’Hour, 2021). This estimation process can be seen as
an automatic procedure of simultaneously selecting and aligning control units
in order to closely match the treated unit. However, the simultaneous
selection and alignment of control units can lead to a loss of efficiency of
the SC method. This procedure also arises another concern is that it is
susceptible to under-fitting due to imperfect pre-treatment fit. The method
requires that the synthetic control’s pre-treatment outcomes closely match the
pre-treatment outcomes for the treated unit (Abadie et al., 2015). This
requirement is often too stringent for using synthetic control alone due to
interpolation bias, as discussed in Kellogg et al. (2021). It is not uncommon
to encounter situations where a linear combination, using nonnegative weights,
of the pre-treatment period outcomes for the control units fails to accurately
approximate to the pre-treatment outcomes of the treated unit.

In this article, we present a straightforward yet effective method called
_Synthetic Matching Control_ (SMC) to address these issues. SMC begins with
the pre-treatment fit for each control units that match the treated unit. It
employs univariate linear regression to better fit on pre-treatment periods,
thereby obtaining a matched characteristic for the control unit. The pre-
treatment fit provides SMC with the flexibility, allowing for the inclusion of
negative coefficients on certain units, thereby reducing interpolation bias.
Subsequently, SMC synthesizes the pre-treatment fits of all control units to
generate an SMC estimator. Our proposed synthesis method builds on the
literature on model averaging (Hansen, 2007), which we study here for
synthetic controls. Like the SC method, SMC can control extrapolation bias by
applying the synthesis method. When the pre-treatment fit on certain control
units is good for the SMC estimator, their corresponding weights are far away
from zero, indicating a heavy reliance on these units. Conversely, if the
weights are zero, it signifies that SMC does not depend on those units.

We conduct both detailed simulation studies and an empirical study of the
economic costs of conflicts in Basque, Spain, to shed light on when the SMC
method performs well. We find evidence that SMC has lower mean-squared
prediction error than alternatives in these studies. Choosing weights in SMC
is more flexible than alternatives so that it can reduce the extrapolation
bias of synthetic estimators with restrictions.

The article is organized as follows. Section 1.1 briefly reviews related work.
Section 2 introduces the set-up and the SC method. Section 3 presents the SMC
method, which includes Section 3.1 on the pre-treatment fit for each control
unit, Section 3.2 on the synthesis method, and Section 3.3 on an
asymptotically optimal weighting method for the SMC method. Two extensions are
considered: the case when units are more than time periods in Section 4.1, and
the incorporation of auxiliary covariates in Section 4.2. Section 5 reports on
extensive simulation studies as well as an application to the Basque dataset.
Finally, Section 6 discusses the limitations of the method and some possible
directions for further research.

###  1.1 Related Work

This article is closely related to the studies that investigate the SC
estimator when the pre-treatment fit is imperfect. One approach to addressing
the issue is to relax the restriction that the weights are nonnegative.
Doudchenko and Imbens (2016) argues that negative weights would be beneficial
in many settings and proposes adding an intercept into the SC problem.
Similarly, Amjad et al. (2018) proposes a denoising algorithm that combines
negative weights with a preprocessing step, and Ferman and Pinto (2021) also
argues that a demeaned version of the SC method is already efficient. Another
approach is to use an outcome model for reducing the imperfect fit. Powell
(2018) allows for extrapolation by constructing the SC unit based on the
fitted values on unit-specific time periods. Ben-Michael et al. (2021)
proposes the augmented synthetic control method, which uses an outcome model
to estimate bias resulting from imperfect pre-treatment fit and de-biases the
original SC estimator. The SMC method relates these estimators in the sense of
addressing the issue of imperfect pre-treatment fit. However, it differs from
them in intent. The concept of SMC aims to mitigate interpolation bias through
unit matching while concurrently reducing extrapolation bias by ensembling all
matched units. We delve these aspects further in Section 3.2.

Several related articles have addressed the challenge of dealing with datasets
that include too many control units, leading to that the solution of the SC
estimator is not unique (Abadie et al., 2015). Robbins et al. (2017) and
Abadie and L’Hour (2021) adapt the original SC proposal to incorporate a
penalty on the weights into the SC optimization problem. Gobillon and Magnac
(2016) makes use of dimension reduction strategies to improve the estimator’s
performance. Doudchenko and Imbens (2016) suggests selecting the set of best
controls by restricting the number of controls allowed to be different from
zero using an $l_{0}$-penalty on the weights. While the SMC method does not
employ a penalty to tackle the problem of an excessive number of unit, the
criterion used in a penalty-style fashion stems from constructing an unbiased
estimator of the risk associated with the synthetic estimator. In situations
where the number of units is large with respect to the number of time periods,
a preprocessing step of screening units can extend the SMC method.

Our article is related to Kellogg et al. (2021), which proposes the matching
and synthetic control (MASC) estimator by using a weighted average of the SC
and matching estimators to balance interpolation and extrapolation bias. Our
SMC method differs from MASC in several ways. First, SMC does not use the
matching estimator and instead considers pre-treatment fit of each control
unit that matches the treated unit. Second, while the MASC estimator is a
weighted average between the SC estimator and the matching estimator, the SMC
estimator is a synthetic estimator that incorporates all matched controls.
Finally, the methods also differ in terms of how the weights are chosen. For
the MASC estimator, the only weight is chosen by the cross-validation method.
In contrast, the SMC method involves solving multiple weights, and we employ
the unbiased risk estimator criterion, as developed in model averaging by
Hansen (2007), to determine these weights.

Our article is also related to Athey et al. (2019) and Viviano and Bradic
(2022) that have considered the benefits of model averaging in the context of
synthetic control. Athey et al. (2019) combines several regularized SC and
matrix completion estimators developed in Doudchenko and Imbens (2016) and
Athey et al. (2021), while Viviano and Bradic (2022) combines a large number
of estimators from the machine learning literature. In contrast, the purpose
of the weighting scheme in our SMC method is to ensemble all matched units and
to mitigate the risk associated with the resulting synthetic estimator. We
utilize the pre-treatment fitted estimators from each control unit as the
estimators to be averaged. This averaging process helps alleviate the risk
inherent in the synthetic estimator. This approach distinguishes our method
from previous methods, which involves averaging various types of estimators or
combining multiple estimators from the machine learning literature.

Finally, in addition to SC-style weighting strategies, there have been
articles that directly use outcome modeling approaches. These include the
panel data approach in Hsiao et al. (2012), the generalized synthetic control
method in Xu (2017), the matrix completion method in Athey et al. (2021), and
the synthetic difference-in-differences method in Arkhangelsky et al. (2021).
In this article we focus on studying the synthetic control framework.

##  2 Synthetic Control Method

We consider the canonical SC panel data setting with $j=1,\cdots,J+1$ units
observed for $t=1,\cdots,T$ time periods. We restrict attention to the case
where a single unit receives treatment, and follow the convention that the
first one $j=1$ is treated and that the remaining $J$ ones are control units.
Let $T_{0}$ be the number of pre-intervention periods, with $1\leq T_{0}<T$.
We adopt the potential outcomes framework (Neyman 1923); the potential
outcomes for unit $j$ in period $t$ under control and treatment are
$Y_{jt}(0)$ and $Y_{jt}(1)$, respectively. Thus, the observed outcomes

$\displaystyle Y_{1t}$ $\displaystyle=\begin{cases}Y_{1t}(0)\text{ if }t\leq
T_{0},\\\ Y_{1t}(1)\text{ if }t>T_{0};\end{cases}$ $\displaystyle Y_{jt}$
$\displaystyle=Y_{jt}(0)\text{ for }j=2,\cdots,J+1.$

We define the treated potential outcome as

$Y_{jt}(1)=Y_{jt}(0)+\tau_{jt},$

where $\tau_{jt}$ is the effect of the intervention for unit $j$ at time $t$.
Since the first unit is treated, the key estimand of interest is
$\tau_{1T}=Y_{1t}(1)-Y_{1t}(0)$. We separate the counterfactual outcome
$Y_{1t}(0)$ into a model component $\mu_{1t}=\mathbb{E}\left[Y_{1t}(0)\right]$
plus an error term $\epsilon_{1t}$:

$\displaystyle Y_{1t}(0)$ $\displaystyle=\mu_{1t}+\epsilon_{1t},\ \
t=1,\cdots,T,$ (1)

where $\epsilon_{1t}$ is a zero mean error term with variance
$\sigma_{t}^{2}=\mathbb{E}\left[\epsilon_{1t}^{2}\right]$. Here we allow for
heteroskedasticity.

Let $\mathbf{y}_{1}$ be a $(T_{0}\times 1)$ vector of pre-intervention
characteristics of the treated unit that we aim to match as closely as
possible, and $\mathbf{Y}_{0}$ be $(T_{0}\times J)$ matrix that contains the
same variables for the control units. A synthetic control is defined as a
weighted average of the control units. Let
${\boldsymbol{w}}=(w_{2},\cdots,w_{J+1})^{\top}$ be the weight vector in the
unit simplex in $\mathbb{R}^{J}$:

${\cal H}_{sc}=\left\\{w_{j}\in[0,1]:\sum_{j=2}^{J+1}w_{j}=1\right\\}.$

In the SC method, the weight vector ${\boldsymbol{w}}$ is chosen to solve the
following optimization problem:

$\displaystyle\tilde{{\boldsymbol{w}}}^{\text{sc}}$
$\displaystyle=\arg\min_{{\boldsymbol{w}}\in{\cal
H}_{sc}}\|\mathbf{y}_{1}-\mathbf{Y}_{0}{\boldsymbol{w}}\|_{\mathbf{V}},$ (2)

where
$\|\mathbf{y}_{1}-\mathbf{Y}_{0}{\boldsymbol{w}}\|_{\mathbf{V}}=\sqrt{(\mathbf{y}_{1}-\mathbf{Y}_{0}{\boldsymbol{w}})^{\top}\mathbf{V}(\mathbf{y}_{1}-\mathbf{Y}_{0}{\boldsymbol{w}})}$
with some symmetric and positive semidefinite matrix $\mathbf{V}$.

The introduction of $\mathbf{V}$ is a crucial step in the SC estimator as it
helps to reduce the mean square error of the estimator. The $\mathbf{V}$
matrix is used to apply a linear transformation on the variables in
$\mathbf{Y}_{0}$ and $\mathbf{y}_{1}$ based on their predictive power on the
outcome. The optimal choice of $\mathbf{V}$ assigns weights that minimize the
mean square error of the synthetic control estimator. A common way is to
choose positive definite and diagonal matrices for $\mathbf{V}$, which results
in the minimum mean squared prediction error of the outcome variable for the
pre-intervention periods (Abadie and Gardeazabal, 2003; Abadie et al., 2010).

Then a synthetic control estimator is constructed by

$\displaystyle\hat{Y}_{1t}(0)$
$\displaystyle=\sum_{j=2}^{J+1}\tilde{w}_{j}^{\text{sc}}Y_{jt},\ \
t\in\\{1,\cdots,T\\}.$

The treatment effect $\tau_{1t}$ is given by the comparison between the
outcome for the treated unit and the outcome for the synthetic control
estimator at time $t\in\\{1,\cdots,T\\}$:

$\hat{\tau}_{1t}=Y_{1t}-\hat{Y}_{1t}(0)=Y_{1t}-\sum_{j=2}^{J+1}\tilde{w}_{j}^{\text{sc}}Y_{jt}.$

The weights $\tilde{{\boldsymbol{w}}}^{\text{sc}}$ in the SC estimator are
typically sparse, meaning that they are only non-zero for a few control units
(Abadie and L’Hour, 2021). This feature is considered as an attractive
property since it provides a way for experts to use their knowledge to
evaluate the plausibility of the resulting estimates (Abadie, 2021).

In the SC method, the optimization problem (2) involves simultaneously
pursuing two objectives: matching and synthesis. The unit simplex serves as
the constraint for the parameter ${\boldsymbol{w}}$ in (2), ensuring that the
estimate $\tilde{{\boldsymbol{w}}}^{\text{sc}}$ often exhibits sparsity. This
sparsity indicates a procedure of selecting control units that are close to
the treated unit. However, the SC method does not account for the potential
fit relation between the treated unit and the control units, which can result
in a loss of accuracy in matching the control units. Furthermore, the unit
simplex ensures that the weights in $\tilde{{\boldsymbol{w}}}^{\text{sc}}$ sum
up to 1, representing a weighted average of the selected control units. As
discussed in Kellogg et al. (2021), the synthesis strategy in the SC method
aims to minimize extrapolation bias but may be susceptible to interpolation
bias.

##  3 Synthetic Matching Control Method

###  3.1 Unit Matching

The goal of unit matching is to establish a correspondence between each
control unit and the treated unit, aiming to minimize the distance between
them. This process entails conducting a univariate regression analysis where
the treated unit is regressed on the control unit. By doing so, we can
estimate the counterfactual outcome for each control unit based on the pre-
intervention regression fit.

We approximate $\mu_{1t}$ in (1) by control unit $j$ in the following working
model

$\displaystyle\mu_{1t}$ $\displaystyle=\theta_{j}Y_{jt}+b_{jt}.$ (3)

where $b_{jt}$ represents the approximation error of $\mu_{1t}$ by
$\theta_{j}Y_{jt}$. We estimate $\theta_{j}\in\mathbb{R}$ using the
unconstrained univariate ordinary least squares method. Without loss of
generality (see Appendix A), we assume that $\mathbf{y}_{1}$ and
$\mathbf{Y}_{0}$ are centered, i.e., $\mathbf{1}^{\top}\mathbf{y}_{1}=0$ and
$\mathbf{1}^{\top}\mathbf{Y}_{0}=\boldsymbol{0}$. We construct the matched
control $\theta_{j}\mathbf{y}_{j}$, where $\theta_{j}$ is obtained by
minimizing the simple Euclidean distance:
$\min_{\theta_{j}\in\mathbb{R}}\|\mathbf{y}_{1}-\theta_{j}\mathbf{y}_{j}\|_{\mathbf{V}}^{2}$.
In practice, we use the procedure employed in the SC method to choose
$\mathbf{V}$ (Abadie and Gardeazabal, 2003; Abadie et al., 2010, 2015).
Because $\mathbf{V}$ can be absorbed into $\mathbf{y}_{1}$ and
$\mathbf{y}_{j}$, without loss of generality we simply rewrite the
minimization as

$\min_{\theta_{j}\in\mathbb{R}}\|\mathbf{y}_{1}-\theta_{j}\mathbf{y}_{j}\|^{2}.$

Using the regression method is to mimic the behavior of the treated unit
before the intervention as closely as possible. It follows the least squares
estimator

$\hat{\theta}_{j}=\left(\mathbf{y}_{j}^{\top}\mathbf{y}_{j}\right)^{-1}\mathbf{y}_{j}^{\top}\mathbf{y}_{1}.$

Consequently, the value of the control unit $j$ based on the matched
characteristic is given by

$\tilde{\mathbf{y}}_{j}=\hat{\theta}_{j}\mathbf{y}_{j}.$

###  3.2 Synthetic Matching Control

We utilize the synthetic method by assigning a weight for the matching
estimator for each control unit. Let
$\boldsymbol{w}=(w_{2},\cdots,w_{J+1})^{\top}$ be the weight vector in the
simplex in $\mathbb{R}^{J}$:

$\mathcal{H}_{J}=\left\\{w_{i}\in[0,1]:j=2,\cdots,J+1\right\\}.$

Note we do not require that $\sum_{j=2}^{J+1}w_{j}=1$. Then the synthetic
matched control estimator is

$\hat{\mathbf{y}}_{1}(\boldsymbol{w})=\sum_{j=2}^{J+1}w_{j}\hat{\theta}_{j}\mathbf{y}_{j}.$
(4)

Let us call it _“Synthetic Matching Control”_ (SMC). In contrast to SC, the
weight in SMC is separate from the regression coefficient. The weight
represents the degree to which the matched characteristic is considered in the
synthesis process, while the regression coefficient reflects the relationship
between the treated and control units. Similar to SC, SMC leverages the
synthetic method to control extrapolation error by assigning more weight to
the matched characteristics that demonstrate higher similarity.

Denote $\boldsymbol{\mu}_{1}=(\mu_{11},\cdots,\mu_{1T_{0}})^{\top}$ and
$\boldsymbol{\epsilon}_{1}=(\epsilon_{11},\cdots,\epsilon_{1T_{0}})^{\top}$.
We decompose the error as follows:

$\displaystyle\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\boldsymbol{\mu}_{1}$
$\displaystyle=\sum_{j=2}^{J+1}w_{j}\hat{\theta}_{j}\mathbf{y}_{j}-\sum_{j=2}^{J+1}w_{j}\mathbf{y}_{1}+\sum_{j=2}^{J+1}w_{j}\mathbf{y}_{1}-\boldsymbol{\mu}_{1}$
$\displaystyle=\sum_{j=2}^{J+1}w_{j}(\hat{\theta}_{j}\mathbf{y}_{j}-\mathbf{y}_{1})+(\sum_{j=2}^{J+1}w_{j}-1)\boldsymbol{\mu}_{1}+\sum_{j=2}^{J+1}w_{j}\boldsymbol{\epsilon}_{1}.$
(5)

Eqn. (3.2) illustrates two components within the error: the first term
represents the interpolation error, resulting from the pre-treatment fitting
errors $\hat{\theta}_{j}\mathbf{y}_{j}-\mathbf{y}_{1}$ with weights $w_{j}$;
the second term is extrapolation error
$(\sum_{j=2}^{J+1}w_{j}-1)\boldsymbol{\mu}_{1}+\sum_{j=2}^{J+1}w_{j}\boldsymbol{\epsilon}_{1}$.
For the SC estimator
$\tilde{\mathbf{y}}_{1}^{\text{sc}}(\boldsymbol{w})=\sum_{j=2}^{J+1}w_{j}\mathbf{y}_{j}$,
the corresponding error is decomposed as

$\displaystyle\tilde{\mathbf{y}}_{1}^{\text{sc}}(\boldsymbol{w})-\boldsymbol{\mu}_{1}$
$\displaystyle=\sum_{j=2}^{J+1}w_{j}(\mathbf{y}_{j}-\mathbf{y}_{1})+(\sum_{j=2}^{J+1}w_{j}-1)\boldsymbol{\mu}_{1}+\sum_{j=2}^{J+1}w_{j}\boldsymbol{\epsilon}_{1}.$
(6)

Comparing (3.2) with (6), we observe two advantages of SMC over SC. First, the
process of unit matching reduces the interpolation error since
$\hat{\theta}_{j}\mathbf{y}_{j}$ is the best linear prediction based on
control unit $j$. This helps to improve the accuracy of the estimated
counterfactual outcome. Second, the constraint $\sum_{j=2}^{J+1}w_{j}=1$ in
the SC method is not necessarily aimed at minimizing the error. Although the
term $(\sum_{j=2}^{J+1}w_{j}-1)\boldsymbol{\mu}_{1}$ disappears when
$\sum_{j=2}^{J+1}w_{j}=1$, it does not guarantee that other terms become
small. In the next section, we adopt the mean-squared error as a measure to
minimize in order to determine the weights ${\boldsymbol{w}}$. This objective
provides a quantitative criterion to optimize the SMC estimator and improve
its overall performance.

###  3.3 Calculation of $\boldsymbol{w}$

From Section 3.1, the prediction of $\boldsymbol{\mu}_{1}$ on the control unit
$j$ is $\tilde{\mathbf{y}}_{j}=\hat{\theta}_{j}\mathbf{y}_{j}$. We rewrite it
as $\tilde{\mathbf{y}}_{j}=\mathbf{H}_{j}\mathbf{y}_{1}$, where
$\mathbf{H}_{j}=\mathbf{y}_{j}\left(\mathbf{y}_{j}^{\top}\mathbf{y}_{j}\right)^{-1}\mathbf{y}_{j}^{\top}$,
implying that (4) is rewritten as

$\hat{\mathbf{y}}_{1}(\boldsymbol{w})=\sum\nolimits_{j=2}^{J+1}w_{j}\mathbf{H}_{j}\mathbf{y}_{1}.$

Define
$\mathbf{H}(\boldsymbol{w})=\sum\nolimits_{j=2}^{J+1}w_{j}\mathbf{H}_{j}$ with
$\ell_{jt}$ as the $t$-th diagonal element of $\mathbf{H}_{j}$, and let
$R(\boldsymbol{w})=\mathbb{E}\left[\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\boldsymbol{\mu}_{1}\|^{2}\right]$
denote the risk. We have that

$\displaystyle\mathbb{E}\left[\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\mathbf{y}_{1}\|^{2}\right]-R(\boldsymbol{w})$
$\displaystyle=\mathbb{E}\left[\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\mathbf{y}_{1}\|^{2}-\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\boldsymbol{\mu}_{1}\|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\mathbf{\epsilon}^{\top}\mathbf{\epsilon}-2\mathbf{\epsilon}^{\top}(\mathbf{H}(\boldsymbol{w})\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{1}+\mathbf{H}(\boldsymbol{w})\mathbf{\epsilon})\right]$
$\displaystyle=\sum\nolimits_{t=1}^{T_{0}}\sigma_{t}^{2}-2\sum\nolimits_{j=2}^{J+1}w_{j}\sum\nolimits_{t=1}^{T_{0}}\sigma_{t}^{2}\ell_{jt}.$
(7)

Eqn. (3.3) demonstrates that the expression

$\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\mathbf{y}_{1}\|^{2}+2\sum\nolimits_{j=2}^{J+1}w_{j}\sum\nolimits_{t=1}^{T_{0}}\sigma_{t}^{2}\ell_{jt}-\sum\nolimits_{t=1}^{T_{0}}\sigma_{t}^{2}$

serves as an unbiased estimator of $R(\boldsymbol{w})$. This motivates the
utilization of the following criterion to obtain ${\boldsymbol{w}}$:

${\cal
C}_{0}(\boldsymbol{w})=\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\mathbf{y}_{1}\|^{2}+2\sum\nolimits_{j=2}^{J+1}\sigma_{j}^{2}w_{j}.$

Here $\sigma_{j}^{2}=\sum\nolimits_{t=1}^{T_{0}}\sigma_{t}^{2}\ell_{jt}$. It
is worth noting that $\sum\nolimits_{t=1}^{T_{0}}\ell_{jt}=1$, in this
context, $\sigma_{j}^{2}$ is defined as a weighted average of $\sigma_{t}^{2}$
over $t=1,\cdots,T_{0}$, where the weights are represented by $\ell_{jt}$.

Replacing $\sigma_{j}^{2}$ by
$\hat{\sigma}^{2}=\|\mathbf{y}_{1}-\mathbf{Y}_{0}[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})]^{-1}\mathbf{Y}_{0}^{\top}\mathbf{y}_{1}\|^{2}$,
where $\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})$ denotes the diagonal
matrix formed by the diagonal elements of
$\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0}$, we propose an approximate Mallows’
$C_{p}$ criterion

${\cal
C}(\boldsymbol{w})=\|\hat{\mathbf{y}}_{1}(\boldsymbol{w})-\mathbf{y}_{1}\|^{2}+2\hat{\sigma}^{2}\sum\nolimits_{j=2}^{J+1}w_{j}.$
(8)

From (8), the weight vector is obtained as

$\hat{\boldsymbol{w}}=\mathop{\arg\min}\limits_{\boldsymbol{w}\in\mathcal{H}_{J}}{{\cal
C}}(\boldsymbol{w}).$

With $\hat{{\boldsymbol{w}}}$ into (4), we obtain the SMC estimator

$\hat{Y}_{1t}(0)=\sum_{j=2}^{J+1}\hat{w}_{j}\hat{\theta}_{j}Y_{jt}\ \
t\in\\{1,\cdots,T\\}.$ (9)

We summarize the procedure of obtaining the SMC estimator as Algorithm 1.

In (9), the SMC estimator is represented as a linear weighting estimator of
the outcomes of control units $Y_{jt}$, similar to the SC estimator. The
weights $\hat{w}_{j}\hat{\theta}_{j}$ can be understood as the solution to a
penalized synthetic control problem. It is important to note that these
weights consist of two components: $\hat{w}_{j}$ and $\hat{\theta}_{j}$. This
formulation allows for negative weights and enables extrapolation beyond the
convex hull of the control units. In unit matching alone, the estimator
$\hat{\theta}_{j}$ allows for arbitrarily weights even when there is no
correlation between the treated unit and the control unit $j$. On the
contrary, by imposing the constraint of the convex hull of the matched control
units, the sum of weights is penalized directly. This constraint effectively
controls the amount of extrapolation error.

Obtain
$\hat{\theta}_{j}=\left(\mathbf{y}_{j}^{\top}\mathbf{y}_{j}\right)^{-1}\mathbf{y}_{j}^{\top}\tilde{{\bf
y}}_{1}$ for $j\in\\{2,\cdots,J+1\\}$.

Solve $\hat{{\boldsymbol{w}}}$ by
$\hat{{\boldsymbol{w}}}=\mathop{\arg\min}\nolimits_{{\boldsymbol{w}}\in\mathcal{H}_{J}}\left\\{\|{\bf
y}_{1}-w_{j}\hat{\theta}_{j}{\bf
y}_{j}\|^{2}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}\right\\}.$

Obtain $\hat{Y}_{1t}(0)=\sum_{j=2}^{J+1}\hat{w}_{j}\hat{\theta}_{j}Y_{jt}$,
$t\in\\{1,\cdots,T\\}$.

Algorithm 1 The SMC estimator

We shall provide a justification: the synthetic estimator with the weights
$\hat{{\boldsymbol{w}}}$ asymptotically achieve the minimum loss of the
infeasible best possible synthetic estimator.

######  Theorem 1.

Denote
$L({\boldsymbol{w}})=\|\hat{\mathbf{y}}_{1}({\boldsymbol{w}})-\boldsymbol{\mu}_{1}\|^{2}$.
Assume that (1) $\max_{t}\mathbb{E}\left[\epsilon_{t}^{4}\right]\leq
c_{1}<\infty$ for some constant $c_{1}$, (2)
$\|\boldsymbol{\mu}_{1}\|^{2}/T_{0}\leq c_{2}<\infty$ for some constant
$c_{2}$, and (3)
$J^{-1}\|\boldsymbol{\mu}_{1}-\mathbf{Y}_{0}[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})]^{-1}\mathbf{Y}_{0}^{\top}\boldsymbol{\mu}_{1}\|^{2}\rightarrow\infty$
as $T_{0}\rightarrow\infty$, then

$\displaystyle\frac{L(\hat{{\boldsymbol{w}}})}{\inf_{{\boldsymbol{w}}\in\mathcal{H}_{J}}L({\boldsymbol{w}})}\rightarrow
1,$

in probability, as $T_{0}\to\infty$.

###### Proof.

The technical proofs are given in Appendix B. ∎

Theorem 1 demonstrates the asymptotic optimality of the proposed method,
showing that that the squared error obtained by $\hat{{\boldsymbol{w}}}$ is
asymptotically equivalent to the infeasible optimal weight vector. The
asymptotic optimality is commonly observed in statistical problems, such as
model selection (Li, 1987) and model averaging (Hansen, 2007). This result
justifies that the SMC estimator is asymptotically optimal in the class of
synthetic estimators where the weight vector ${\boldsymbol{w}}$ is restricted
to the set $\mathcal{H}_{J}$.

The conditions of $\max_{t}\mathbb{E}\left[\epsilon_{t}^{4}\right]\leq c_{1}$
and $\|\boldsymbol{\mu}_{1}\|^{2}/T_{0}\leq c_{2}$ are quite mild since they
only require bounded fourth moments of errors and that
$\|\boldsymbol{\mu}_{1}\|^{2}=O(T_{0})$, respectively. The key condition,
$J^{-1}\|\boldsymbol{\mu}_{1}-\mathbf{Y}_{0}[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})]^{-1}\mathbf{Y}_{0}^{\top}\boldsymbol{\mu}_{1}\|^{2}\rightarrow\infty$
as $T_{0}\rightarrow\infty$, means that the squared model approximation error
is large relative to the number of control units. This condition described is
typically considered to be mild in the context of the synthetic control
problem. This is because achieving a perfect approximation through univariate
regression on a simple control unit is rarely attainable.

##  4 Extensions

In this section, we consider two elaborations to the basic setup. First, we
extend it to cases where units are more than time periods. Second, we extend
it by incorporating auxiliary covariates.

###  4.1 Screening Units When They are Too Many

We extend the application of the SMC method to cases where $J\geq T_{0}$ or
$J\approx T_{0}$. To accomplish this, we propose a practical procedure that
involves screening the units using the sure independent ranking and screening
(SIRS) method (Zhu et al., 2011) to reduce the number of units. In high-
dimensional statistics, Theorems 2 and 3 in Zhu et al. (2011) indicate that
SIRS can reduce the dimensionality without losing any active variables with a
probability approaching one. We prefer SIRS over the original sure
independence screening proposed by Fan and Lv (2008) because it allows us to
assume that no linear candidate model is correct. In model (3), any working
model is an approximation of the expected counterfactual value $\mu_{1t}$.

For applying SIRS into the control units, we assume that $\mu_{1t}$ depends
only on some of the control units, called as active units, in this study. SIRS
screens the units based on the magnitude of the following statistics instead
of the marginal correlation,

$\tilde{\eta}_{j}=\frac{1}{T_{0}}\sum_{t=1}^{T_{0}}\left\\{\frac{1}{T_{0}}\sum_{l=1}^{T_{0}}Y_{jt}I_{(-\infty,\
Y_{1t})}(Y_{1l})\right\\}^{2},\quad\text{for}\quad j=2,\cdots,J+1.$

Derivation and interpretation of this statistics can be found in Zhu et al.
(2011). We use the statistics $\tilde{\eta}_{j}$ for screening units, then
obtain a set that involves any activate units. Once the screened units are
reduced, we perform the SMC method on these units to obtain the estimator. We
summarize it as Algorithm 2 below.

Step 1: Screen units by the SIRS method to get the subset $\mathbf{Y}_{s}$
from $\mathbf{Y}_{0}$;

Step 2: Perform Algorithm 1 on $\mathbf{y}_{1}$ and $\mathbf{Y}_{s}$.

Algorithm 2 The SMC estimator when control units are too many

###  4.2 Incorporating auxiliary covariates

We have focused on matching pre-treatment values of the outcome variable. In
practice, we typically observe a set of auxiliary covariates as well. For
example, in the study of Proposition 99, Abadie et al. (2010) considers the
following covariates: average retail price of cigarettes, per capita state
personal income, per capita beer consumption, and the percentage of the
population age 15–24.

It is natural to incorporate auxiliary covariates in applying the SMC method.
For unit $j$, denote $\mathbf{x}_{j}$ as a $(p\times 1)$ vector of observed
covariates that are not affected by the intervention. Let
$\mathbf{X}=(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{J+1})$.
Analogous to the SC method (Abadie et al., 2010), We define the augmented
$(T_{d}\times 1)$, where $T_{d}=T_{0}+p$, vector of pre-intervention
characteristics for the treated unit
$\mathbf{z}_{1}=(\mathbf{y}_{1}^{\top},\mathbf{x}_{1})^{\top}\in\mathbb{R}^{T_{d}}$.
Similarly, $\mathbf{Z}_{0}$ is a $(T_{d}\times J_{0})$ matrix that contains
the same variables for the control units.

We apply Algorithm 1 on $\mathbf{z}_{1}$ and $\mathbf{Z}_{0}$ to obtain
$\hat{w}_{j}^{(\mathbf{z})}$ and $\hat{\theta}_{j}^{(\mathbf{z})}$, and then
obtain the SMC estimator

$\hat{Y}_{1t}(0)=\sum_{j=2}^{J+1}\hat{w}_{j}^{(\mathbf{z})}\hat{\theta}_{j}^{(\mathbf{z})}Y_{jt}\
\ t\in\\{1,\cdots,T\\}.$ (10)

We summarize it as Algorithm 3 below.

Step 1: Combine $\mathbf{y}_{1}$ with $\mathbf{x}_{1}$ to obtain
$\mathbf{z}_{1}$; Similarly combine $\mathbf{Y}_{0}$ with $\mathbf{X}_{0}$ to
obtain $\mathbf{Z}_{0}$;

Step 2: Perform Algorithm 1 on $\mathbf{z}_{1}$ and $\mathbf{Z}_{0}$, and then
obtain $\hat{w}_{j}^{(\mathbf{z})}$ and $\hat{\theta}_{j}^{(\mathbf{z})}$;

Step 3: Obtain the SMC estimator according to (10).

Algorithm 3 The SMC estimator when auxiliary covariates are incorporated

##  5 Empirical Studies

In this section, we conduct extensive Monte Carlo simulation studies to assess
the performance of various methods, finding where and how the SMC estimator
performs compared to existing estimators, and subsequently we perform an
empirical analyse on a real dataset to examine the behavior of the SMC method.

###  5.1 Simulation Studies

Now we investigate the finite sample performance of alternative estimators in
two simulation experiments, one using a generative model and the other using a
working model. We compare several representative synthetic estimators,
including: (a) the original SC method (SC) in Abadie et al. (2010), (b) the
de-meaned SC (dSC) in Ferman and Pinto (2021), (c) the augmented SC Method
(ASC) in Ben-Michael et al. (2021), (d) the matching and SC method (MASC) in
Kellogg et al. (2021), (e) OLS in Hsiao et al. (2012), and (f) the constrained
lasso (lasso) in Chernozhukov et al. (2021).

Generative Model. In this experiment, all units are generated according to the
factor model as follows

$Y_{jt}(0)=\alpha_{t}+\lambda_{j}F_{t}+\epsilon_{jt},$

where time specific terms $\alpha_{t}\sim\mathcal{N}(0,1)$, and unobserved
factors $F_{t}\sim\mathcal{N}(0,1)$. Similarly as Li (2020), we consider three
sets of factor loadings $\boldsymbol{\lambda}$. (1)
$\boldsymbol{\lambda}_{1}$: $\lambda_{j}=1$ for $j=1,\cdots,7$, and
$\lambda_{j}=0$ for $j=8,\cdots,J+1$; (2) $\boldsymbol{\lambda}_{2}$:
$\lambda_{1}=3$, $\lambda_{j}=1$ for $j=2,\cdots,7$, and $\lambda_{j}=0$ for
$j=8,\cdots,J+1$; (3) $\boldsymbol{\lambda}_{3}$: $\lambda_{1}=3$,
$\lambda_{j}=1$ for $j=2,\cdots,J+1$; For (1), all nonzero factor loadings are
set to be ones so that both treated and control units with nonzero loadings
are drawn from a common distribution. In contrast, for (2), treated and
control units are drawn from two heterogeneous distributions, since loadings
for the treated unit all equal to 3, and the control units’ nonzero loadings
all equal to 1. The setting of $\boldsymbol{\lambda}_{3}$ is similar to
$\boldsymbol{\lambda}_{2}$, except that control units’ loadings all equal to
1. We set $T=50$ with $T_{0}=40$ and $J=20$. The errors
$\epsilon_{jt}\sim\mathcal{N}(0,\sigma^{2})$, where we use three values of
$\sigma$ values, 1, 0.5, and 0.1, to investigate the impact of $\sigma$.

To evaluate each estimator, we compute the mean squared prediction error
(MSPE), which is defined as
$\text{MSPE}=(T-T_{0})^{-1}\sum\limits_{t=T_{0}+1}^{T}\|\hat{Y}_{1,t}(0)-Y_{1,t}(0)\|^{2}$,
by calculating the average loss across 500 simulations. The results are
reported in Table 1. For the homoscedastic setting of
$\boldsymbol{\lambda}_{1}$, the SC, dSC, ASC, and MASC methods perform
similarly as lasso, and are better than SMC. For the heteroscedastic settings
of $\boldsymbol{\lambda}_{2}$ and $\boldsymbol{\lambda}_{3}$, SMC is better
than other methods except for the case of $\boldsymbol{\lambda}_{2}$ and
$\sigma=0.1$, where it is slightly worse than ASC. Especially for the dense
setting $\boldsymbol{\lambda}_{3}$, the SC, dSC, ASC, and MASC methods work as
badly as lasso, while SMC works pretty well. Meanwhile, OLS works well in this
setting, but is still worse than SMC.

Comparing the results across various $\sigma$ values, we find that the above
observations hold true. Notably, we observe that the MSPE values of both the
SMC and OLS estimators approach zero as $\sigma$ decreases from 1 to 0.1
across all three scenarios of $\boldsymbol{\lambda}$. However, we do not
observe this convergence in the SC, dSC, MASC, and lasso methods for the
heteroscedastic settings, nor in ASC for the setting of
$\boldsymbol{\lambda}_{3}$.

In order to assess the impact of $(T,J)$, we consider three scenarios: $(T,J)$
takes on values of $(100,20)$, $(50,10)$, and $(50,50)$. We set $T_{0}=T-10$.
The results are reported in Table 2. When comparing the cases of
$(T,J)=(100,20)$, $(50,10)$ or $(50,50)$ with the case of $(T,J)=(50,20)$ in
Table 1, we find that these methods consistently demonstrate robust
performance across different values of $(T,J)$. In the
$\boldsymbol{\lambda}_{1}$ setting, SMC performs similarly to the SC, dSC,
ASC, MASC, and lasso methods. However, it outperforms the other methods in the
$\boldsymbol{\lambda}_{2}$ and $\boldsymbol{\lambda}_{3}$ settings. It is
noteworthy that in the case of $(50,50)$, we have $T_{0}<J$. For this case, we
employ the SMC estimator by applying Algorithm 2 with SIRS preprocessing on
screening units as described in Section 4.1. These results obtained in the
$(50,50)$ case highlight the effectiveness of the SIRS preprocessing step on
screening units within the extended SMC method.

Table 1: The post-period MSPE of alternative estimators. $(\boldsymbol{\lambda},\sigma)$ | SC | dSC | ASC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---|---
$(\boldsymbol{\lambda}_{1},1)$ | 1.521 | 1.493 | 1.474 | 1.562 | 2.566 | 1.564 | 1.578
$(\boldsymbol{\lambda}_{1},0.5)$ | 0.285 | 0.290 | 0.279 | 0.412 | 0.477 | 0.295 | 0.336
$(\boldsymbol{\lambda}_{1},0.1)$ | 0.014 | 0.015 | 0.014 | 0.021 | 0.031 | 0.015 | 0.017
$(\boldsymbol{\lambda}_{2},1)$ | 6.572 | 6.805 | 3.561 | 6.656 | 3.827 | 6.839 | 2.528
$(\boldsymbol{\lambda}_{2},0.5)$ | 5.758 | 5.871 | 0.900 | 5.828 | 1.061 | 5.809 | 0.742
$(\boldsymbol{\lambda}_{2},0.1)$ | 5.246 | 5.364 | 0.024 | 5.253 | 0.031 | 5.347 | 0.021
$(\boldsymbol{\lambda}_{3},1)$ | 19.80 | 19.97 | 21.87 | 19.94 | 2.752 | 19.48 | 2.010
$(\boldsymbol{\lambda}_{3},0.5)$ | 18.10 | 18.29 | 18.81 | 18.24 | 0.721 | 18.22 | 0.489
$(\boldsymbol{\lambda}_{3},0.1)$ | 17.89 | 17.98 | 19.12 | 17.96 | 0.034 | 18.01 | 0.021
Table 2: The post-period MSPE of alternative estimators under various $(T,J)$. $(T,J)$ | $(\boldsymbol{\lambda},\sigma)$ | SC | dSC | ASC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---|---|---
| $(\boldsymbol{\lambda}_{1},1)$ | 1.221 | 1.238 | 1.214 | 1.291 | 1.444 | 1.260 | 1.291
| $(\boldsymbol{\lambda}_{1},0.5)$ | 0.314 | 0.321 | 0.309 | 0.338 | 0.393 | 0.321 | 0.331
| $(\boldsymbol{\lambda}_{1},0.1)$ | 0.013 | 0.013 | 0.013 | 0.014 | 0.015 | 0.013 | 0.014
| $(\boldsymbol{\lambda}_{2},1)$ | 20.316 | 20.252 | 3.223 | 20.133 | 3.528 | 20.259 | 2.952
(100,20) | $(\boldsymbol{\lambda}_{2},0.5)$ | 19.524 | 19.349 | 0.823 | 19.447 | 0.896 | 19.298 | 0.716
| $(\boldsymbol{\lambda}_{2},0.1)$ | 18.947 | 18.879 | 0.037 | 18.932 | 0.029 | 18.850 | 0.027
| $(\boldsymbol{\lambda}_{3},1)$ | 14.557 | 21.023 | 15.304 | 14.098 | 1.813 | 20.462 | 1.730
| $(\boldsymbol{\lambda}_{3},0.5)$ | 12.967 | 19.572 | 14.023 | 12.926 | 0.455 | 19.478 | 0.423
| $(\boldsymbol{\lambda}_{3},0.1)$ | 12.734 | 19.179 | 13.593 | 12.732 | 0.018 | 19.165 | 0.017
| $(\boldsymbol{\lambda}_{1},1)$ | 1.166 | 1.202 | 1.152 | 1.251 | 1.527 | 1.213 | 1.204
| $(\boldsymbol{\lambda}_{1},0.5)$ | 0.360 | 0.378 | 0.351 | 0.403 | 0.416 | 0.371 | 0.375
| $(\boldsymbol{\lambda}_{1},0.1)$ | 0.013 | 0.013 | 0.012 | 0.014 | 0.015 | 0.012 | 0.013
| $(\boldsymbol{\lambda}_{2},1)$ | 26.252 | 26.847 | 4.592 | 26.931 | 3.716 | 27.192 | 3.246
(50,10) | $(\boldsymbol{\lambda}_{2},0.5)$ | 25.219 | 26.184 | 1.035 | 25.263 | 1.033 | 26.292 | 0.952
| $(\boldsymbol{\lambda}_{2},0.1)$ | 25.755 | 26.596 | 0.038 | 25.674 | 0.032 | 26.588 | 0.028
| $(\boldsymbol{\lambda}_{3},1)$ | 13.669 | 15.136 | 14.827 | 13.458 | 2.683 | 15.000 | 2.531
| $(\boldsymbol{\lambda}_{3},0.5)$ | 11.900 | 13.765 | 13.155 | 11.820 | 0.718 | 13.839 | 0.689
| $(\boldsymbol{\lambda}_{3},0.1)$ | 11.543 | 13.792 | 12.797 | 11.557 | 0.022 | 13.760 | 0.021
| $(\boldsymbol{\lambda}_{1},1)$ | 1.375 | 1.428 | 1.368 | 1.477 | 5.867 | 1.383 | 1.482
| $(\boldsymbol{\lambda}_{1},0.5)$ | 0.334 | 0.335 | 0.335 | 0.358 | 1.541 | 0.337 | 0.349
| $(\boldsymbol{\lambda}_{1},0.1)$ | 0.012 | 0.013 | 0.012 | 0.012 | 0.055 | 0.012 | 0.013
| $(\boldsymbol{\lambda}_{2},1)$ | 11.753 | 10.580 | 2.366 | 11.825 | 12.150 | 11.114 | 2.906
(50,50) | $(\boldsymbol{\lambda}_{2},0.5)$ | 11.591 | 10.436 | 0.917 | 11.728 | 2.562 | 10.153 | 0.851
| $(\boldsymbol{\lambda}_{2},0.1)$ | 11.472 | 10.417 | 0.040 | 11.428 | 0.121 | 10.286 | 0.030
| $(\boldsymbol{\lambda}_{3},1)$ | 42.628 | 50.944 | 44.088 | 42.713 | 6.288 | 51.190 | 2.645
| $(\boldsymbol{\lambda}_{3},0.5)$ | 41.936 | 50.807 | 43.524 | 41.814 | 1.644 | 50.473 | 0.610
| $(\boldsymbol{\lambda}_{3},0.1)$ | 41.520 | 50.186 | 42.391 | 41.572 | 0.059 | 50.225 | 0.024

Working Model. Our second simulation is based on the following working model:

$\mathbf{y}_{1}=\mathbf{Y}_{0}\boldsymbol{\theta}+{\boldsymbol{\epsilon}}.$

Here rows of the $T$ by $J$ predictor matrix $\mathbf{Y}_{0}$ is generated
from a multivariate normal distribution with mean $\mathbf{0}$ and covariance
matrix $\boldsymbol{\Sigma}$ with the $(i,j)$th element
$\sigma_{ij}=\rho^{|i-j|}$, where the value of $\rho$ is set to 0.8. As the
first simulation does, $(T,J)$ is set to $(50,20)$, where $T_{0}=40$. The true
coefficients $\boldsymbol{\theta}=(\theta_{1},\cdots,\theta_{J})^{\top}$ are
set to $(c/7,\cdots,c/7,0,\cdots,0)^{\top}$, where the first 7 coefficients
are non-zero with $c>0$. We control the summation of these coefficients by
setting various $c$: 0.5, 1, and 2. When $c=1$, the SC method is expected to
be perfect. For $c=0.5$ or 2, there is gap between the restriction that the
weights sum up to one and the real approximating model. The errors
${\boldsymbol{\epsilon}}$ is considered to be the heteroscedastic setting of
$\epsilon_{t}\sim N(0,\sigma_{t}^{2})$ with
$\sigma_{t}^{2}=\sigma_{0}^{2}\frac{\|\mathbf{y}^{t}\|^{2}}{J+\|\mathbf{y}^{t}\|^{2}}$,
where $\mathbf{y}^{t}$ means the $t$-th row of $\mathbf{Y}_{0}$. Here the
value of $\sigma_{0}^{2}$ is selected to control the coefficient of
determination $R^{2}$ to vary among 0.4, 0.6, and 0.8.

The results in the working model setting are reported in Table 3, which shows
the post-period MSPE values of alternative estimators. For the case where
$c=1$, the SC, dSC, ASC, MASC and lasso methods perform similarly, but achieve
a lower MSPE than SMC, which aligns with our expectations. However, in cases
where $c=0.5$, SMC achieves a lower risk than the SC, dSC, ASC, and MASC
methods, and performs similarly as lasso. Additionally, for cases where $c=2$,
SMC achieves the best among all methods. These observations remain consistent
across various $R^{2}$ values.

Table 3: The post-period MSPE of alternative estimators in the working model setting. $c$ | $R^{2}$ | SC | dSC | ASC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---|---|---
| 0.4 | 1.1059 | 1.1279 | 1.1227 | 1.1280 | 1.8637 | 1.0634 | 1.1615
1 | 0.6 | 0.4941 | 0.5091 | 0.4971 | 0.5744 | 1.0746 | 0.5089 | 0.5349
| 0.8 | 0.1947 | 0.2039 | 0.1918 | 0.2115 | 0.3898 | 0.2045 | 0.2270
| 0.4 | 0.2237 | 0.2327 | 0.2268 | 0.2349 | 0.2922 | 0.1664 | 0.1566
0.5 | 0.6 | 0.1440 | 0.1489 | 0.1451 | 0.1510 | 0.1026 | 0.0768 | 0.0702
| 0.8 | 0.1110 | 0.1141 | 0.1118 | 0.1265 | 0.0492 | 0.0431 | 0.0272
| 0.4 | 4.9207 | 5.1163 | 4.5520 | 5.0719 | 7.5369 | 5.1716 | 4.3806
2 | 0.6 | 2.7837 | 2.8202 | 2.7010 | 2.7860 | 3.6843 | 2.8342 | 2.4110
| 0.8 | 1.3589 | 1.3843 | 1.1473 | 1.3118 | 1.3252 | 1.3425 | 0.8332

###  5.2 The Basque dataset

We study the effect of terrorism on per capita GDP in Basque, Spain. The
Basque dataset is from Abadie and Gardeazabal (2003). It consists of per
capita GDP of 17 regions in Spain from 1955 to 1997, and 12 other covariates
of each region over the same time interval, representing education,
investment, sectional shares, and population density in each region. We
incorporate auxiliary covariates which include averages for the 13
characteristics from 1960 to 1969, and scale each covariate so that it has
equal variance of outcomes. In this study, the treated unit is the Basque
Country, and the treatment is the onset of separatist terrorism, which begins
in 1970.

Placebo Analyses. Similar to Abadie and Gardeazabal (2003), we conduct a
placebo study to compare alternative estimators in the real data. We perform
placebo analyses on each region, excluding Basque, as the placebo region. We
calculate the mean squared prediction error (MSPE) for each region by taking
the differences between its actual and fitted outcome paths in each of the
post-period years (1970-1997), squaring these differences and then averaging
them among these years. The results of our analysis are presented in Table 4,
which shows that, on average, SMC tends to have the lowest MSPE. In addition,
we include the pre-period fit of these estimator in Table 6 of Appendix to
further demonstrate their performance. Interestingly, we observe that SMC does
not exhibit the best pre-period fit on average (it is the second best), while
ASC demonstrates the best pre-period fit on average. This observation suggests
that SMC is less prone to over-fitting compared to ASC.

Table 4: Performance (MSPE) of alternative estimators in the placebo study. region | SC | dSC | ASC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---|---
Andalucia | 0.41 | 0.15 | 0.17 | 0.32 | 5.18 | 0.13 | 0.22
Aragon | 0.03 | 0.12 | 0.06 | 0.04 | 0.03 | 0.02 | 0.01
Asturias | 0.71 | 0.56 | 3.40 | 0.73 | 0.30 | 0.44 | 0.67
Baleares | 2.12 | 3.68 | 1.24 | 2.12 | 2.51 | 4.73 | 0.38
Canarias | 0.07 | 0.10 | 0.35 | 0.07 | 0.96 | 0.02 | 0.25
Cantabria | 0.37 | 0.65 | 0.90 | 0.34 | 1.87 | 0.56 | 0.12
Leon | 0.01 | 0.12 | 0.08 | 0.01 | 0.15 | 0.13 | 0.13
Mancha | 0.07 | 0.02 | 0.04 | 0.04 | 0.49 | 0.39 | 0.04
Cataluna | 0.44 | 0.03 | 0.14 | 0.44 | 0.73 | 1.33 | 0.04
Valenciana | 0.15 | 0.14 | 0.09 | 0.04 | 0.08 | 0.03 | 0.31
Extremadura | 0.74 | 0.06 | 0.17 | 0.74 | 0.99 | 0.63 | 0.32
Galicia | 0.01 | 0.02 | 0.04 | 0.01 | 0.07 | 0.01 | 0.11
Madrid | 0.11 | 0.38 | 3.75 | 0.11 | 0.48 | 0.16 | 0.48
Murcia | 0.21 | 0.27 | 0.07 | 0.20 | 1.18 | 0.03 | 0.51
Navarra | 0.04 | 0.04 | 0.03 | 0.05 | 0.04 | 0.12 | 0.03
Rioja | 0.04 | 0.03 | 0.08 | 0.19 | 0.09 | 0.30 | 0.04
average | 0.35 | 0.40 | 0.66 | 0.34 | 0.95 | 0.56 | 0.23

Synthetic Basque. We estimate the effect of exposure to terrorism on GDP per
capita in Basque, Spain. We present the GDP per capita for both Basque and its
synthetic control, generated using the SC and SMC methods, in Figure 1. We
also show the difference between the actual values and the synthetic values.
Due to saving space and for ease of comparison, we just report the comparison
between SMC and the original SC method, ignoring other estimators. We observe
that SMC performs better in tracking the post-intervention trend. To gain a
better understanding of the SMC method, we inspect the comprehensive weight
$\theta_{j}w_{j}$ assigned to each control unit $j$ and compare it with the
weights obtained from other methods. The results are reported in Table 5. The
SC, dSC, and MASC methods assign positive weights to just two control units.
The lasso method assigns positive weights to four control units. Comparing to
them, SMC assigns non-zero weight to four units. It is important to note that
OLS has no zero weight as it is not constrained. This suggests that, compared
with other methods with constraints, SMC’s greater flexibility may bring
greater predictive power.

Figure 1: Study of the Basque Country. Left: Actual and counterfactual per capita GDP of the Basque Country. Right: The difference of actual and counterfactual values (ATT). Table 5: Estimates of weights for alternative estimators in the Basque study. The comprehensive weights $\theta_{j}w_{j}$ and the coefficients are reported for SMC and OLS, respectively. region | SC | dSC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---
Andalucia | 0 | 0 | 0 | 0.217 | -0.112 | 0
Aragon | 0 | 0 | 0 | -3.059 | 0 | 0
Asturias | 0 | 0 | 0 | 1.460 | 0 | 0
Baleares | 0 | 0.582 | 0 | -0.365 | 0.365 | 0
Canarias | 0 | 0 | 0 | -0.245 | 0 | 0
Cantabria | 0 | 0 | 0 | 0.048 | 0.048 | 0.270
Leon | 0 | 0 | 0 | 0.080 | 0 | 0
Mancha | 0 | 0 | 0 | 0.861 | 0 | 0
Cataluna | 0.851 | 0 | 0.851 | -0.174 | 0.174 | 0
Valenciana | 0 | 0 | 0 | 1.766 | 0.065 | 0.058
Extremadura | 0 | 0 | 0 | -0.401 | 0 | 0
Galicia | 0 | 0 | 0 | -0.396 | 0 | 0
Madrid | 0.149 | 0.418 | 0.149 | 0.533 | 0 | 0.342
Murcia | 0 | 0 | 0 | -0.547 | 0 | 0
Navarra | 0 | 0 | 0 | 1.302 | 0.234 | 0
Rioja | 0 | 0 | 0 | 0.673 | 0 | 0.567
Intercept | - | -0.335 | - | -2.580 | 1.097 | -0.192

##  6 Discussion

This paper makes three contributions: (1) We propose a simple and effective
method, synthetic matching control, by synthesizing the matched controls. (2)
We determine the weights by minimizing the unbiased risk estimate criterion.
We demonstrate that this method is asymptotically optimal, achieving the
lowest possible squared error. (3) We expand the method in several domains,
including conducting inference, extending it to cases where control units is
more than time periods, and incorporating auxiliary covariates.

There are several potential directions for future work. First, we focus on the
simple linear regression to reduce the issue of imperfect fit. Consequently,
the synthesis method relies on the simple linear regression. However, the SMC
method may be applicable for more complex data structures, such as discrete,
count, or hierarchical outcomes. Therefore, extending the method to broader
regression models is an interesting direction. Second, for settings with
multiple treated units, we can fit SMC separately for each treated unit, as in
Abadie (2021). However, this approach brings a loss of efficiency due to the
correlation of treated units. Therefore, efficiently extending the method to
multiple treated units is a worthy direction. Third, a set of auxiliary
covariates is available, we pool the auxiliary covariates and the outcomes
togother to conduct the SMC method. However, this approach may bring extra
risk when the linear approximation relation in the covariates is different
from that in the outcomes. Therefore, exploring ways to incorporate auxiliary
covariates into the SMC method while minimizing such risks is another worthy
problem for future research. Finally, extending it to more complicated
situations, such as staggered adoption where units take up the treatment at
different times (Ben-Michael et al., 2022), is another challenging direction.

## References

  * Abadie (2021) Abadie, A. (2021), “Using Synthetic Controls: Feasibility, Data Requirements, and Methodological Aspects,” Journal of Economic Literature, 391–425.
  * Abadie et al. (2010) Abadie, A., Diamond, A., and Hainmueller, J. (2010), “Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of Cali- fornia’s Tobacco Control Program,” Journal of the American Statistical Association, 493–505.
  * Abadie et al. (2015) — (2015), “Comparative Politics and the Synthetic Control Method,” American Journal of Political Science, 495–510.
  * Abadie and Gardeazabal (2003) Abadie, A. and Gardeazabal, J. (2003), “The Economic Costs of Conflict: A Case Study of the Basque Country,” The American Economic Review, 113–132.
  * Abadie and L’Hour (2021) Abadie, A. and L’Hour, J. (2021), “A Penalized Synthetic Control Estimator for Disaggregated Data,” Journal of the American Statistical Association, 536, 1817–1834.
  * Amjad et al. (2018) Amjad, M., Shah, D., and Shen, D. (2018), “Robust Synthetic Control,” The Journal of Machine Learning Research, 19, 802–852.
  * Arkhangelsky et al. (2021) Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., and Wager, S. (2021), “Synthetic Difference in Differences,” American Economic Review, 111, 4088–4118.
  * Athey et al. (2021) Athey, S., Bayati, M., Doudchenko, N., Imbens, G., and Khosravi, K. (2021), “Matrix Completion Methods for Causal Panel Data Models,” Journal of the American Statistical Association, 116, 1716–1730.
  * Athey et al. (2019) Athey, S., Bayati, M., Imbens, G., and Qu, Z. (2019), “Ensemble Methods for Causal Effects in Panel Data Settings,” AEA Papers and Proceedings, 109, 65–70.
  * Ben-Michael et al. (2021) Ben-Michael, E., Feller, A., and othstein, J. (2021), “The Augmented Synthetic Control Method,” Journal of the American Statistical Association, 536, 1789–1803.
  * Ben-Michael et al. (2022) Ben-Michael, E., Feller, A., and Rothstein, J. (2022), “Synthetic Controls With Staggered Adoption,” Journal of the Royal Statistical Society, Series B, 84, 351–381.
  * Chernozhukov et al. (2021) Chernozhukov, V., Wüthrich, K., and Zhu, Y. (2021), “An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls,” Journal of the American Statistical Association, 536, 1849–1864.
  * Doudchenko and Imbens (2016) Doudchenko, N. and Imbens, G. W. (2016), “Balancing, Regression, Difference-In-Differences and Synthetic Control Methods: A Synthesis,” Working Paper 22791, National Bureau of Economic Research.
  * Fan and Lv (2008) Fan, J. and Lv, J. (2008), “Sure Independence Screening for Ultrahigh Dimensional Feature Space (with discussion),” Journal of the Royal Statistical Society, Series B, 70, 849–911.
  * Ferman and Pinto (2021) Ferman, B. and Pinto, C. (2021), “Synthetic Controls with Imperfect Pretreatment Fit,” Quantitative Economics, 12, 1197–1221.
  * Gobillon and Magnac (2016) Gobillon, L. and Magnac, T. (2016), “Regional Policy Evaluation: Inter- active Fixed Effects and Synthetic Controls,” Review of Economics and Statistics, 98, 535–551.
  * Hansen (2007) Hansen, B. (2007), “Least Squares Model Averaging,” Econometrica, 75, 1175–1189.
  * Hsiao et al. (2012) Hsiao, C., Ching, S., and Wan, K. (2012), “A Panel Data Approach for Program Evaluation: Measuring the Benefits of Political and Economic Integration of Hong Kong With Mainland China,” Journal of Applied Econometrics, 27, 705–740.
  * Kellogg et al. (2021) Kellogg, M., Mogstad, M., Guillaume A. Pouliot, G., and Torgovitsky, A. (2021), “Combining Matching and Synthetic Control to Tradeoff Biases From Extrapolation and Interpolation,” Journal of the American Statistical Association, 536, 1804–1814.
  * Li (2020) Li, K. (2020), “Statistical Inference for Average Treatment Effects Estimated by Synthetic Control Methods,” Journal of the American Statistical Association, 115, 2068–2083.
  * Li (1987) Li, K.-C. (1987), “Asymptotic optimality for $C_{p}$, $C_{l}$, Cross-validation and Generalized Cross-validation: Discrete index set,” The Annals of Statistics, 15, 958–975.
  * Powell (2018) Powell, D. (2018), “Imperfect Synthetic Controls: Did the Massachusetts Health Care Reform Save Lives?” Tech. rep., Santa Monica, CA: RAND Corporation.
  * Robbins et al. (2017) Robbins, M., Saunders, J., and Kilmer, B. (2017), “A Framework for Synthetic Control Methods With High-Dimensional, Micro-Level Data: Evaluating a Neighborhood-Specific Crime Intervention,” Journal of the American Statistical Association, 112, 109–126.
  * Viviano and Bradic (2022) Viviano, D. and Bradic, J. (2022), “Synthetic Learner: Model-free Inference on Treatments over Time,” Journal of Econometrics.
  * Whittle (1960) Whittle, P. (1960), “Bounds for the Moments of Linear and Quadratic Forms in Independent Variables,” Theory of Probability and its Applications, 5, 302–305.
  * Xu (2017) Xu, Y. (2017), “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models,” Political Analysis, 57–76.
  * Zhu et al. (2011) Zhu, L., Li, L., Li, R., and Zhu, L. (2011), “Model-free Feature Screening for Ultrahigh-dimensional Data,” Journal of the American Statistical Association, 106, 1464–1475.

##  Appendix A Centralization Without Loss of Generality

No intercept is used in the matching estimator since the assumption that
$\mathbf{y}_{1}$ and $\mathbf{y}_{j}$ are centered. If $\mathbf{y}_{1}$ and
$\mathbf{y}_{j}$ are not centered, we should add an intercept term
$b\in\mathbb{R}$ in the matching estimator, which is constructed by

$\displaystyle\mathbf{y}_{1}$ $\displaystyle=b+\theta_{j}\mathbf{y}_{j}.$ (11)

The parameters $b$ and $\theta_{j}$ are obtained by minimizing the simple
Euclidean distance:

$\min_{b\in\mathbb{R},\theta\in\mathbb{R}}\|\mathbf{y}_{1}-b-\theta_{j}\mathbf{y}_{j}\|.$

Denote $\tilde{\mathbf{y}}_{j}=\mathbf{y}_{j}-\bar{y}_{j}\mathbf{1}$ and
$\tilde{\mathbf{y}}_{1}=\mathbf{y}_{1}-\bar{y}_{1}\mathbf{1}$. The
minimization above implies that

$\hat{\theta}_{j}=(\tilde{\mathbf{y}}_{j}^{\top}\tilde{\mathbf{y}}_{j})^{-1}\tilde{\mathbf{y}}_{j}^{\top}\tilde{\mathbf{y}}_{1},\quad\tilde{b}_{j}=\bar{y}_{1}-\bar{\mathbf{y}}_{j}^{\top}\hat{\theta}_{j}.$

The Synthetic Matched Control estimator is

$\hat{\mathbf{y}}_{1}(\boldsymbol{w})=\sum_{j=2}^{J+1}w_{j}(\bar{y}_{1}+\tilde{\mathbf{y}}_{j}\hat{\theta}_{j})=\bar{y}_{1}+\sum_{j=2}^{J+1}w_{j}\tilde{\mathbf{y}}_{j}\hat{\theta}_{j}.$
(12)

Comparing (12) with (4), it is without loss of generality that we assume that
$\mathbf{y}_{1}$ and $\mathbf{y}_{j}$ are centered.

##  Appendix B Proof of Theorem 1

The solution
$\hat{\boldsymbol{\theta}}=(\hat{\theta}_{2},\cdots,\hat{\theta}_{J+1})^{\top}$
satisfies:

$\displaystyle\hat{\boldsymbol{\theta}}=$
$\displaystyle\left[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})\right]^{-1}\mathbf{Y}_{0}^{\top}\mathbf{y}_{1}.$
(13)

Inserting $\hat{\boldsymbol{\theta}}$ into the criterion
$\mathcal{C}_{c}({\boldsymbol{w}})$, we have

$\displaystyle\mathcal{C}_{c}({\boldsymbol{w}})$
$\displaystyle=\left\|\mathbf{y}_{1}-\mathbf{Y}_{0}\mathbf{W}\left[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})\right]^{-1}\mathbf{Y}_{0}^{\top}\mathbf{y}_{1}\right\|^{2}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}.$
(14)

By the above analysis on Algorithm 1, we just need to prove the asymptotic
optimality by analyzing (14). In the theorem, we have the following
conditions:

$\displaystyle(C1)\quad\max_{t}\mathbb{E}\left[\epsilon_{t}^{4}\right]\leq
c_{1}<\infty\text{ for some constant }c_{1}.$ (15)
$\displaystyle(C2)\quad\|\boldsymbol{\mu}_{1}\|^{2}/T_{0}\leq
c_{2}<\infty\text{ for some constant }c_{2}.$ (16) $\displaystyle(C3)\quad
J^{-1}\|\boldsymbol{\mu}_{1}-\mathbf{Y}_{0}[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})]^{-1}\mathbf{Y}_{0}^{\top}\boldsymbol{\mu}_{1}\|^{2}\rightarrow\infty\text{
as }T_{0}\rightarrow\infty.$ (17)

Let
$\mathbf{M}({\boldsymbol{w}})=\mathbf{Y}_{0}\mathbf{W}\left[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})\right]^{-1}\mathbf{Y}_{0}^{\top}$
and ${\bf A}({\boldsymbol{w}})={\bf I}-\mathbf{M}({\boldsymbol{w}})$. We have

$\displaystyle\hat{\bf
y}_{1}({\boldsymbol{w}})=\mathbf{M}({\boldsymbol{w}}){\bf y}_{1},$ (18)
$\displaystyle L({\boldsymbol{w}})=\|\hat{\bf
y}_{1}({\boldsymbol{w}})-\boldsymbol{\mu}_{1}\|^{2}=\|\mathbf{M}({\boldsymbol{w}}){\bf
y}_{1}-\boldsymbol{\mu}_{1}\|^{2}.$ (19)

From (19), a simple calculation follows

$\displaystyle R({\boldsymbol{w}})=E\\{L({\boldsymbol{w}})\\}=\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|^{2}+\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}\geq\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|^{2}.$ (20)

Thus, combing (20) and Condition (17), we have

$\displaystyle\inf_{w\in\mathcal{H}_{J}}J^{-1}R({\boldsymbol{w}})\rightarrow\infty.$
(21)

Denote
${\bf\Sigma}=\text{diag}\\{\sigma_{1}^{2},\cdots,\sigma_{T_{0}}^{2}\\}$. From
(14) and (19), we have

$\displaystyle\mathcal{C}_{c}({\boldsymbol{w}})-L({\boldsymbol{w}})$
$\displaystyle=\|\hat{\bf
y}_{1}({\boldsymbol{w}})-\boldsymbol{\mu}_{1}-{\boldsymbol{\epsilon}}\|^{2}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}-\left\|\mathbf{M}({\boldsymbol{w}}){\bf
y}_{1}-\boldsymbol{\mu}_{1}\right\|^{2}$ $\displaystyle=-2\left\\{\hat{\bf
y}_{1}({\boldsymbol{w}})-\boldsymbol{\mu}_{1}\right\\}^{\top}{\boldsymbol{\epsilon}}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}+\left\|{\boldsymbol{\epsilon}}\right\|^{2}$
$\displaystyle=-2{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\bf
y}_{1}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}+2\boldsymbol{\mu}_{1}^{\top}{\boldsymbol{\epsilon}}+\left\|{\boldsymbol{\epsilon}}\right\|^{2}$
$\displaystyle=-2{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}})\boldsymbol{\mu}_{1}-2{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}+2\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}+2\boldsymbol{\mu}_{1}^{\top}{\boldsymbol{\epsilon}}+\left\|{\boldsymbol{\epsilon}}\right\|^{2},$
(22)

where the last three terms in the last line do not involve ${\boldsymbol{w}}$.
From (19) and (20), we have

$\displaystyle R({\boldsymbol{w}})-L({\boldsymbol{w}})=$
$\displaystyle\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|^{2}+\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}-\left\|\mathbf{M}({\boldsymbol{w}}){\bf
y}_{1}-\boldsymbol{\mu}_{1}\right\|^{2}$ $\displaystyle=$
$\displaystyle\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|^{2}+\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}-\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}-\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}\right\|^{2}$
$\displaystyle=$
$\displaystyle-\left\|\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}\right\|^{2}+\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}+2{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}.$ (23)

Thus, to prove the asymptotic optimality, We need only to verify that

$\begin{array}[]{lll}&(\text{i})\ \sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}R^{-1}({\boldsymbol{w}})\left|{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right|=o_{P}(1),\\\
&(\text{ii})\ \sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}R^{-1}({\boldsymbol{w}})[\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}-{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}]=o_{P}(1),\\\
&(\text{iii})\ \sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}R^{-1}({\boldsymbol{w}})\left\|\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}\right\|^{2}=o_{P}(1),\\\
&(\text{iv})\ \sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}R^{-1}({\boldsymbol{w}})\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}=o_{P}(1),\\\
&(\text{v})\ \sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}R^{-1}({\boldsymbol{w}})\left|{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right|=o_{P}(1).\end{array}$

Let
$\mathbf{H}=\mathbf{Y}_{0}\left[\text{diag}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})\right]^{-1}\mathbf{Y}_{0}^{\top}$.
A direct simplification yields

$\displaystyle\mathbf{M}({\boldsymbol{w}})$ $\displaystyle\leq\mathbf{H},$
(24)

where the inequality is in Loewner ordering. Given Condition (15), applying
Theorem 2 of Whittle (1960) leads to that, for some constant $0<\eta<\infty$,

$\displaystyle\mathbb{E}\left[\left({\boldsymbol{\epsilon}}^{\top}\mathbf{H}{\boldsymbol{\epsilon}}-\mathbb{E}\left[{\boldsymbol{\epsilon}}^{\top}\mathbf{H}{\boldsymbol{\epsilon}}\right]\right)^{2}\right]\leq\eta(\hbox{trace}\\{\mathbf{H}{\bf\Sigma}\mathbf{H}({\boldsymbol{w}})\\}).$
(25)

Thus, we have

$\displaystyle\mathbf{\epsilon}^{\top}\mathbf{H}\mathbf{\epsilon}=$
$\displaystyle
E(\mathbf{\epsilon}^{\top}\mathbf{H}\mathbf{\epsilon})+O_{P}\bigg{\\{}\sqrt{Var(\mathbf{\epsilon}^{\top}\mathbf{H}\mathbf{\epsilon})}\bigg{\\}}$
$\displaystyle=$
$\displaystyle\hbox{trace}\\{{\bf\Sigma}\mathbf{H}\\}+O_{P}\bigg{\\{}\sqrt{\hbox{trace}\\{{\bf\Sigma}\mathbf{H}\\}}\bigg{\\}}$
$\displaystyle=$ $\displaystyle O_{p}(J).$ (26)

From (21) and (B), we have

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{{\boldsymbol{\epsilon}}^{\top}\mathbf{H}{\boldsymbol{\epsilon}}}{R({\boldsymbol{w}})}\leq\frac{O_{P}(J)}{R({\boldsymbol{w}})}=o_{P}(1).$
(27)

From (21) and (27), we have

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left|{\boldsymbol{\epsilon}}^{\top}\mathbf{H}\mathbf{A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right|}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left\|\mathbf{H}{\boldsymbol{\epsilon}}\right\|\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left\|\mathbf{H}{\boldsymbol{\epsilon}}\right\|}{R^{1/2}({\boldsymbol{w}})}=o_{P}(1).$

Thus, (i) is proved.

From Condition (16), we have

$\displaystyle\hat{\sigma}^{2}=\frac{\|\mathbf{y}_{1}-\mathbf{Y}_{0}(\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0})^{-1}\mathbf{Y}_{0}^{\top}\mathbf{y}_{1}\|^{2}}{T_{0}-J}=O_{P}(1).$
(28)

It follows that

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{|\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}-{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}|}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\hat{\sigma}^{2}{\boldsymbol{w}}^{\top}\mathbf{1}}{R({\boldsymbol{w}})}+\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{{\boldsymbol{\epsilon}}^{\top}\mathbf{H}{\boldsymbol{\epsilon}}}{R({\boldsymbol{w}})}=o_{P}(1),$

where the last step is from (21), (27), and (28). Thus, (ii) is proved.

From (24), we have

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\|\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}\|^{2}}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{{\boldsymbol{\epsilon}}^{\top}\mathbf{H}{\boldsymbol{\epsilon}}}{R({\boldsymbol{w}})}=o_{P}(1),$

where the last step is from (21). Thus, (iii) is proved.

From (24), we have

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\hbox{trace}\left\\{\mathbf{M}({\boldsymbol{w}}){\bf\Sigma}\mathbf{M}^{\top}({\boldsymbol{w}})\right\\}}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\hbox{trace}\\{\mathbf{H}{\bf\Sigma}\\}}{R({\boldsymbol{w}})}=o_{P}(1),$

where the last step is from (21) and Condition (15). Thus, (iv) is proved.

From (20) and (24), we have

$\displaystyle\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left|{\boldsymbol{\epsilon}}^{\top}\mathbf{M}({\boldsymbol{w}}){\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right|}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left\|\mathbf{M}({\boldsymbol{w}}){\boldsymbol{\epsilon}}\right\|\left\|{\bf
A}({\boldsymbol{w}})\boldsymbol{\mu}_{1}\right\|}{R({\boldsymbol{w}})}\leq\sup\limits_{{\boldsymbol{w}}\in{\cal
H}_{J}}\frac{\left\|\mathbf{H}{\boldsymbol{\epsilon}}\right\|}{R^{1/2}({\boldsymbol{w}})}=o_{P}(1),$

where the last step is from (27). It means that (vi) is proved. Therefore, the
theorem is proved.

Table 6: The pre-preiod fit of alternative estimators in the placebo study. The fit by the residual sum of squares (RSS) is measured.  region | SC | dSC | ASC | MASC | OLS | lasso | SMC
---|---|---|---|---|---|---|---
Andalucia | 1.26e-3 | 5.01e-2 | 3.96e-3 | 2.44e-3 | 1.76e-2 | 8.17e-2 | 5.87e-5
Aragon | 4.12e-4 | 1.95e-4 | 1.35e-4 | 7.52e-4 | 6.46e-5 | 1.15e-3 | 6.39e-4
Asturias | 9.15e-5 | 4.96e-5 | 9.76e-3 | 1.09e-4 | 1.34e-3 | 2.42e-3 | 4.81e-5
Baleares | 9.52e-2 | 3.04e-1 | 4.73e-3 | 9.52e-2 | 5.35e-3 | 6.2e-1 | 1.64e-2
Canarias | 1.32e-3 | 2.78e-3 | 1.16e-3 | 1.32e-3 | 2.17e-3 | 3.8e-2 | 1.35e-3
Cantabria | 1.09e-3 | 4.36e-4 | 4.57e-3 | 5.67e-4 | 4.52e-3 | 9.89e-4 | 1.39e-4
Leon | 4.60e-4 | 8.42e-4 | 8.76e-5 | 5.79e-4 | 3.18e-4 | 6.26e-3 | 8.72e-3
Mancha | 4.14e-3 | 3.72e-4 | 5.96e-4 | 8.4e-3 | 1.05e-3 | 9.41e-3 | 2.33e-3
Cataluna | 1.36e-2 | 1.83e-4 | 8.95e-4 | 1.36e-2 | 1.39e-3 | 1.04e-1 | 9.36e-5
Valenciana | 9.81e-4 | 6.96e-4 | 7.42e-4 | 4.51e-3 | 1.45e-4 | 2.48e-3 | 3.13e-3
Extremadura | 1.15e-1 | 1.08e-2 | 1.40e-3 | 1.15e-1 | 2.04e-3 | 8.53e-2 | 4.65e-4
Galicia | 3.01e-4 | 3.45e-4 | 3.07e-5 | 3.01e-4 | 8.51e-5 | 1.22e-2 | 3.03e-4
Madrid | 7.21e-1 | 5.61e-2 | 2.25e-2 | 7.21e-1 | 2.14e-2 | 1.19e-1 | 2.34e-1
Murcia | 1.37e-3 | 5.33e-4 | 1.44e-3 | 2.63e-3 | 3.11e-3 | 9.07e-3 | 1.24e-2
Navarra | 2.82e-4 | 3.01e-4 | 6.16e-5 | 5.09e-4 | 1.63e-4 | 1.13e-2 | 2.36e-4
Rioja | 7.17e-4 | 8.40e-4 | 1.30e-4 | 3.11e-3 | 4.95e-4 | 1.38e-2 | 3.75e-4
average | 5.98e-2 | 2.68e-2 | 3.27e-3 | 6.06e-2 | 3.83e-2 | 6.98e-2 | 1.75e-2