# Topological superfluids in two-dimensional Fermi gas with Rashba spin-orbit
coupling

Xiaosen Yang  Department of Physics, Jiangsu University, Zhenjiang, 212013, P.
R. China  Ho-Kin Tang  Shenzhen JL Computational Science And Applied Research
Institute, Shenzhen, 518131, P. R. China  Noah Fan Qi Yuan  Shenzhen JL
Computational Science And Applied Research Institute, Shenzhen, 518131, P. R.
China  Beibing Huang  Department of Physics, Yancheng Institute of Technology,
Yancheng, 224051, P. R. China  Guangcan Guo  CAS Key Lab of Quantum
Information, University of Science and Technology of China, Hefei, 230026,
P.R. China and Synergetic Innovation Center of Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei, 230026, P.R.
China  Ming Gong  CAS Key Lab of Quantum Information, University of Science
and Technology of China, Hefei, 230026, P.R. China and Synergetic Innovation
Center of Quantum Information and Quantum Physics, University of Science and
Technology of China, Hefei, 230026, P.R. China<EMAIL_ADDRESS>Hai-Qing
Lin  Shenzhen JL Computational Science And Applied Research Institute,
Shenzhen, 518131, P. R. China  Beijing Computational Science Research Center,
Beijing, 100084, P. R. China

###### Abstract

The realization of spin-orbit coupling (SOC) in ultracold atoms has triggered
an intensive exploring of topological superfluids in the degenerate Fermi
gases based on mean-field theory, which has not yet been reported in
experiments. Here, we demonstrate the topological phase transitions in the
system via the numerically exact quantum Monte Carlo method. Without prior
assumptions, our unbiased real-space calculation shows that spin-orbit
coupling can stabilize an unconventional pairing in the weak SOC regime, in
which the Fulde-Ferrell-Larkin-Ovchinnikov pairing coexists with the Bardeen-
Cooper-Schrieffer pairing. Furthermore, we use the jumps in the spin
polarization at the time-reversal invariant momenta to qualify the topological
phase transition, where we find the critical exponent deviated from the mean-
field theory. Our results pave the way for the searching of unconventional
pairing and topological superfluids with degenerate Fermi gases.

## Introduction

Figure 1: Competition between the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO)
pairing and the Bardeen-Cooper-Schrieffer (BCS) pairing: The figure shows pair
momentum distribution $n_{pair}(\textbf{k})$ for the spin-orbit coupling
$\lambda=0,0.18,0.4$ (from left to right) with attractive interaction $U=3.5$,
Zeeman field $h=1.7$, inverse temperature $\beta=10$ and filling factor
$n=0.375$. The FFLO-type pairing with nonzero momentum peaks at pair momentum
distribution while the BCS-type pairing with the peak at the zero momentum.
The FFLO+BCS pairing has both the zero momentum and nonzero momentum peaks at
pair momentum distribution. With increasing spin-orbit coupling, the type of
pairing goes from FFLO, FFLO+BCS, to BCS. In the inset, we illustrate the
Fermi surface for the pairing in different phases with the arrows pointing to
the direction of $\langle S_{k}\rangle$. We simulate the system in a $8\times
16$ periodic supercell.

Spin-orbit coupling (SOC) has triggered extensive concerns in condensed
matters and ultracold atoms for its essential role in many exotic phenomena,
especially for the realization of topological phases such as topological
insulator [1, 2], topological superconductors [2, 3], spin Hall effects [4,
5], Weyl semimetal [6, 7, 8], Floquet topological phase [9, 10, 11] and so on.
In ultracold atoms, SOC can be realized both in fermionic [12, 13] and bosonic
[14, 15] systems by Raman couplings. These milestone breakthroughs have shed
the light on exploring various novel phases with different topology and
pairings in experiments [5, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

The ultracold atoms also provide a clean and ideal platform to study the
unconventional pairings and related phenomena in the Fermi gases [26]. In
solid materials, the Fermi energy is much larger than the Zeeman splitting and
SOC, thus only the Bardeen-Cooper-Schrieffer (BCS) pairing is favored in most
of the materials. As a result, while the Fulde-Ferrell-Larkin-Ovchinnikov
(FFLO) [27, 28] superconductors with finite momentum pairings have been
predicted through the competition between magnetism and pairing energies in
the 1960s, it is hard to be realized in solid materials. The degenerate Fermi
gases can be used to explore this kind of pairing for the comparable energies
between the Fermi energy, Zeeman field and SOC.

By the tuning of the Fermi surfaces from the SOC and in-plane Zeeman field,
which breaks the inversion symmetry [29], the FFLO state can be realized,
yielding the so-called topological FFLO superfluids [30, 31]. The Zeeman field
is found to drive the topological phase transitions within the superfluid
phase [32, 33, 34, 35, 36]. The competition between FFLO and BCS pairing near
the pairing-normal phase transition has been studied within mean field theory
[37, 38], while results deep in the pairing phase at low temperatures are
still needed. In the presence of SOC, the transition from the FFLO to BCS
phase becomes continuous in one dimensional systems [39, 40]. However, the
finite temperature phase transition only occurs in high dimensions, where no
concrete evidence of the continuous behavior has been found.

In this manuscript, we extend the mean-field study of the superfluid phase to
the numerically exact quantum Monte Carlo (QMC) method [41], in which the
pairing correlation is calculated and the topological invariant in terms of
quantized $\pi$ Berry phase on the spin-non-degenerate Fermi surface is
applied. The spin-orbit coupled lattice model with attractive interactions
have been investigated using the auxiliary field QMC method [42, 43, 44, 45],
in which no infamous sign problem happens in the spin-balanced case [46].
These previous studies show that SOC can optimize the superconducting
temperature, leads to rich pairing structure like the mixture of spin-singlet
and spin-triplet [42, 43], suppresses charge density wave order near the half-
filling, and even contributes to the topological signature of edge current
[44]. However, when we include the Zeeman field, the sign problem restricts
the range of parameter that we can explore, but we can still study FFLO by
selecting parameters with acceptable sign [47].

Our results show that the ratio of Zeeman field $h$ and SOC strength $\lambda$
play an important role in determining the pairing properties, from the FFLO
pairing in the small $\lambda/h$ limit to the BCS pairing in the large
$\lambda/h$ limit. We present a strong evidence that the FFLO+BCS phase is
stabilized between two limits. In large SOC regime, we characterize the
topological properties of the superfluids by spin polarization at time-
reversal-invariant momenta. We investigate the phase diagram and find the
critical exponent associated with the scaling is different from the result
given by the mean-field study. The topological superfluids have not yet been
realized in experiments. Our results provide strong evidences for the
realization of topological superfluids with the current available interactions
in experiments, thus pave the way for the realization of these exotic
superfluids with degenerate Fermi gases using the state-of-art setups.

## Model

We start with the Rashba spin-orbit coupled Fermi gas loaded on a 2D square
optical lattice which can be described by the following fermionic Hubbard
Hamiltonian:

$\displaystyle H$ $\displaystyle=$
$\displaystyle-t\sum_{\langle{\bm{i},\bm{j}}\rangle}c^{\dagger}_{\bm{i}s}c_{\bm{j}s}+i\lambda\sum_{\langle{\bm{i},\bm{j}}\rangle}c_{\bm{i}s}^{\dagger}(\bm{\hat{e}}_{\bm{i},\bm{j}}\times\bm{\sigma})^{s,s^{\prime}}_{z}c_{\bm{j}s^{\prime}}-U\sum_{\bm{i}}n_{\bm{i}\uparrow}n_{\bm{i}\downarrow}-\mu\sum_{\bm{i}}n_{\bm{i}}+h\sum_{\bm{i}}c^{\dagger}_{\bm{i}s}\sigma^{s,s^{\prime}}_{z}c_{\bm{i}s^{\prime}},$

where $c^{\dagger}_{\bm{i}s}$ and $c_{\bm{i}s}$ are the creation and
annihilation operators for fermions at site $\bm{r_{i}}$ with spin
$s=(\uparrow,\downarrow)$. $n_{\bm{i}}=\sum_{s}n_{\bm{i}s}$ is the total
density operator in real space on site $\bm{r_{i}}$ with
$n_{\bm{i}s}=c^{\dagger}_{\bm{i}s}c_{\bm{i}s}$. The hopping amplitude between
the nearest-neighbor sites $\langle{\bm{i},\bm{j}}\rangle$ is $t$, which we
take equal to unity to set the energy scale. $\bm{\hat{e}}_{\bm{i},\bm{j}}$ is
the unit vector connecting sites $\bm{r_{i}}$ and $\bm{r_{j}}$.
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are the Pauli matrices.
$\lambda$ and $U(U>0)$ stand for Rashba SOC strength and the strength of on-
site attractive interaction. $\mu$ and $h$ are the chemical potential and
Zeeman field, respectively. The desired filling factor
$n=\sum_{\bm{i}s}\langle n_{\bm{i}s}\rangle/N$ on an $N$-site periodic
supercell lattice is obtained by tuning the chemical potential $\mu$, where
$\langle\dots\rangle$ is the thermodynamic average as defined in Eq. (5) in
Method section. Without interaction, the dispersion relation of square lattice
model is $\epsilon_{k}=-2t(\cos{k_{x}}+\cos{k_{y}})-\mu$. At low filling, the
model shows similar dispersion relation to the Fermi gas with
$\epsilon_{k}=k^{2}/2m-\mu^{\prime}$ with $\mu^{\prime}=\mu+2t$ and $m=1/2t$.
We study the regime away from half filling to detect the unconventional
pairing. In the current experiments of Fermi gas in the 2D optical lattices
[48, 49], $t\sim 0.2-1$ kHz using $U\sim(7-8)t$; and the coldest temperature
achieved is about $k_{B}T\sim 150-250$ Hz ($k_{B}$ is the Boltzman constant),
which corresponds to $\beta t\sim 4$. Thus in our manuscript, we focus on the
regime with $\beta t\sim 4-12$.

## Results

### Competition between FFLO pairing and BCS pairing

Figure 2: The FFLO+BCS pairing response to the change of the Zeeman field $h$
and the filling $n$. (keeping $\lambda$ fixed at 0.18 as in Fig. 1(b)) We plot
the pair momentum distribution $n_{pair}$ along the $k_{y}$ direction with
$k_{x}=0$. (a) With decreasing $h$ from $1.7$ to $1.5$, the ratio of
$h/\lambda$ drops, the BCS pairing peaked at $k=(0,0)$ appears (dotted blue
line). (b) Keeping the ratio of $h/\lambda$ the same, FFLO+BCS pairing
persists for increasing $n$ from 0.375 to 0.41 (dashed red line). However, if
the filling $n$ is small, only BCS pairing happens because of the band
structure (dotted blue line).

To investigate the effect of SOC on the pairing properties, we measure the
pair momentum distribution as the Fourier transform of the pairing correlation
[39, 50, 51, 52]:

$\displaystyle
n_{pair}(\bm{k})=\sum_{i,j}e^{i\bm{k}\cdot(\bm{r_{i}}-\bm{r_{j}})}\langle\Delta^{\dagger}_{\bm{i}}\Delta_{\bm{j}}\rangle,$
(1)

where
$\Delta_{\bm{i}}^{\dagger}=c_{\bm{i}\uparrow}^{\dagger}c_{\bm{i}\downarrow}^{\dagger}$
creates a singlet pair at the site $\bm{r_{i}}$ and $\bm{k}$ is the momentum.
This distribution has a simple interpretation. For BCS pairing, it exhibits
peaks at zero momentum; while for FFLO pairing between
$|\bm{k},\uparrow\rangle$ and $|-\bm{k}+\bm{Q},\downarrow\rangle$, it exhibits
peaks at $\pm\bm{Q}$.

The pair momentum distribution is shown in Fig. 1 for various SOC with
$U=3.5$, $h=1.7$, $\beta=10$ and filling factor $n=0.375$. In the absence of
SOC, the Zeeman field lifts the degeneracy of spin. The pairing occur between
different Fermi circles with different spin and hence with finite momentum.
The pair momentum distribution is different from the balanced case and has a
circular ridge shape. In Fig. 1(a), the peak of the distribution is at the rim
of radius $|\bm{Q}|=|\bm{k}_{F\downarrow}|-|\bm{k}_{F\uparrow}|$, where
$\bm{k}_{F\uparrow}$ and $\bm{k}_{F\downarrow}$ are the Fermi momenta of
spin-$\uparrow$ and -$\downarrow$ Fermi circles, respectively. This nonzero
momentum peak indicates the stability of the FFLO pairing in the system [53].

As shown in Fig. 1(b), in the presence of SOC, the FFLO peak drops slightly.
Interestingly, a zero momentum peak emerges in the distribution when the
strength of SOC increases. Thus, the pair momentum distribution has a novel
form: the circular ridge and the zero momentum peak coexist in the momentum
space. This distribution is different from the FFLO and BCS, but contains the
signatures of both two type pairing. This indicates the pairing of BCS and
FFLO is compatible in a small range of SOC strength, corresponding to the
FFLO-BCS pairing. In Fig. 1(c), as the SOC increases further, the zero
momentum peak will ascend further until the finite circular ridge disappears
in the distribution, resembling BCS pairing. We find that BCS pairing
dominates for the moderate and strong strength of SOC. With increasing SOC,
the pairing undergoes a transition from FFLO to BCS through FFLO-BCS.

The exotic pairing arises from the dramatic change of spin texture induced by
competition between SOC and Zeeman field [37, 38]. The insets of Fig. 1
illustrate how the spin texture in different Fermi circles changes with SOC as
Zeeman field is fixed nonzero. Without SOC, the pairing occurs between the
spin-$\uparrow$ and spin-$\downarrow$ from different bands, resulting in the
finite momentum pairing of FFLO. With increasing SOC, two helical bands
consist of the mixture of the spin-$\uparrow$ and spin-$\downarrow$ fermions.
The FFLO pairing still occurs between two helical bands as the interband
pairing, as $z$-component of the spin persists to be polarized by Zeeman
field. However, BCS pairing arises from the $xy$-component of the spin within
the same band as the intraband pairing. Both types of pairing coexist in the
FFLO+BCS phase. Further increasing SOC, the $xy$-component of the spin becomes
dominant on the Fermi surface of the helical bands, so the intraband BCS
pairing dominates.

To investigate the stability of the exotic FFLO+BCS pairing, we show how the
Zeeman field and filling factor change the $n_{pair}$. In Fig. 2(a), with
decreasing the Zeeman field, the BCS pairing is enhanced and the zero momentum
peak becomes larger. With increasing the Zeeman field, the total pairing is
suppressed, and the FFLO peak momentum becomes larger. In Fig. 2(b), when the
filling factor decreases, the FFLO pairing is suppressed. For the small
filling, only the lower chirality band is occupied. Only the intraband BCS
pairing is allowed and the interband FFLO pairing is prohibited.

### Topologically nontrivial superfluids and the topological phase transitions

Figure 3: Phase diagram for fermions on square lattice with Rashba spin-orbit
coupling and Zeeman field. There is three topological superfluid phases with
nonzero winding number $W=\pm 1,-2$. Increasing the Zeeman field, there is
topological phase transition from topologically trivial superfluid ($W=0$)
into topologically nontrivial superfluid phases ($W\neq 0$). The on-site
attraction strength is $U=4$ and the spin-orbit coupling strength is
$\lambda=1$. We classify the topological phase according to the spin
polarization $\delta n(\textbf{k})=\langle
n_{\downarrow}(\textbf{k})\rangle-\langle n_{\uparrow}(\textbf{k})\rangle$,
which has jumps at time-reversal invariant momenta through the topological
phase transition. The jumps are notated as ($\star$) in (b)-(e). The inverse
temperature is set as $\beta=12$.

The system can realize the topological nontrivial superfluid states which
break the time-reversal symmetry. The topological nontrivial phase is chiral
superfluid phase, which can be characterized by winding number defined from
the Green’s function at all frequencies[54, 4]:

$\displaystyle W=\frac{1}{24\pi^{2}}\int
d^{2}kd\omega\mathrm{Tr}[\epsilon^{\mu\nu\rho}(\mathcal{G}\partial_{\mu}\mathcal{G}^{-1})(\mathcal{G}\partial_{\nu}\mathcal{G}^{-1})(\mathcal{G}\partial_{\rho}\mathcal{G}^{-1})],$
(2)

where the winding number $W$ is zero for the topologically trivial phase and
an nonzero integer for the topologically nontrivial phase. The winding number
changes when the system undergoes a topological phase transition. The winding
number can be simplified by the ‘topological Hamiltonian’, which depends only
on the zero-frequency Green’s function [55, 56]. But it is difficult to
extract the anomalous Green’s function from the QMC simulation. Thus, we do
not measure the winding number directly, instead we make use of the spin
polarization at the time reversal invariant (TRI) momenta [57, 58], since they
are easy to measure in both QMC simulations and cold atom experiments. For a
square lattice system, the winding number can be determined by the spin
polarization at four TRI momenta $\textbf{K}=(0,0),(0,\pi),(\pi,0),(\pi,\pi)$
as follows

$\displaystyle W$ $\displaystyle=$
$\displaystyle\sum_{\textbf{K}}\cos(K_{x})\cos(K_{y})\delta n(\textbf{K}),$
(3)

where $\delta
n(\textbf{K})=n_{\downarrow}(\textbf{K})-n_{\uparrow}(\textbf{K})$ is the spin
polarization at the TRI momentum K. Prior to these results, we have verified
that without SOC (setting $\lambda=0$), the jump of spin difference at the TRI
momenta is impossible, as expected for the standard BCS superfluids with
$W=0$. Meanwhile, without Zeeman field (see Fig. 3 (a) at $h=0$), the jump of
spin difference is also impossible, yielding a trivial superfluid with $W=0$.
There are also several special conditions: (i) The system is always in a
normal phase without cold atoms ($n=0$); (ii) This model possesses particle-
hole symmetry, thus the physics of $n$ and $(2-n)$ yields the same phase
diagram with $W\rightarrow W$; (iii) This model is also symmetric about
$h\rightarrow-h$ (with $W\rightarrow-W$) and $\lambda\rightarrow-\lambda$
(with $W\rightarrow W$).

Fig. 3 shows the phase diagram in the $n-h$ plane with $U=4$ and $\lambda=1$.
The topological invariants of the superfluid phases are determined by the spin
polarization at the four TRI momenta $\delta n(\textbf{K})$. Without the
Zeeman field, the system preserves the time-reversal symmetry (TRS). Then, the
superfuild phase is topologically trivial. When we introduce the Zeeman field,
the TRS is broken and the superfluid phase can be topologically nontrivial.

For small Zeeman field, the spin polarization $\delta n(\textbf{K})$ at four
TRI momenta are also zero and the superfluid phase is topologically trivial
with zero winding number ($W=0$). With increasing Zeeman field, the spin
polarization at the four TRI momenta jumps from $0$ into $1$ and the
superfluid becomes topologically nontrivial with nonzero topological
invariant. We mark the unit spin polarization at the TRI momenta as stars in
Fig. 3 (b)-(e). For large filling, the spin polarization remains unchanged at
$[(0,0),(\pi,\pi)]$, while has a jump from $0$ into $1$ at $[(0,\pi),(\pi,0)]$
by increasing the Zeeman field as shown in Fig. 3 (b). Therefore, the phase is
topologically nontrivial with nonzero winding number ($W=-2$). By increasing
the Zeeman field further, the spin polarization jumps at momentum $(0,0)$ as
shown in Fig. 3 (c). Then, the topological invariant changes from $-2$ to $-1$
at the topological phase transition boundary. For intermediate filling, the
spin polarization jumps from $0$ into $1$ at $(0,0)$ and remains unchange at
other TRI momenta as increasing of the Zeeman field. Thus, the winding number
of the topological superfluid is $W=1$ as shown in Fig. 3 (e). For the
topologically nontrivial phases, the robust topological edge states also can
stabilizes at the boundary of the system. The three topological superfluid
phases have distinct different behavior in momentum distribution and can be
easily determined in cold atom experiments. The pairing is easily to be broken
by the Zeeman field and the phase is normal state with negligible pair
correlation with $\mathrm{max}[n_{pair}]<0.1$ for the statistical error of the
QMC method.

We also compare the QMC phase diagram with the phase diagram obtained from the
mean-field approach (See Supplementary Material [59]). There are some
significant differences, outlined as following. (A) In the mean-field theory,
the order parameter is destroyed when $h>3.8$; however, in the QMC simulation,
the superfluids can be survived even when $h>5$. At half filling, a much
larger Zeeman field is required to destroy the pairing. (B) At low filling
($n\sim 0$), a finite Zeeman field ($h\sim 2$) is required to destroy the
pairing, while from QMC, the critical Zeeman field is proportional to the
filling factor $n$. These quantitative difference may imply the limitation of
the mean-field theory in the understanding of the topological superfluids in
degenerate Fermi gases; and the qualitatively similarities between these two
approaches shed light on the realization of the topological superfluids based
on the current experimental setups.

Figure 4: Topological phase transition and its associated scaling We find the
magnetization, equivalent to $\delta n(0,0)$, as a functions of Zeeman field
for inverse temperature $\beta=4,8,12$ with $U=4$, $\lambda=1$ and $n=0.4$.
The polarization increases from 0 into 1 as the increasing of the Zeeman
field. The vertical lines with $\delta n(0,0)=0.5$ defines the phase boundary
$h_{c}=1.64$, $1.72$ and $1.73$ for respective temperatures. The inset shows
the universal scaling law of spin polarization using Eq. 4, by which we obtain
the critical exponent $\nu=4/3\pm 0.1$. The Green curve is the scaling
function.

To probe scaling near the phase transition, we measure the spin polarization
$m\equiv\delta n(0,0)$ at zero momentum as a function of Zeeman field for
different temperature $\beta=1/k_{B}T=4,8,12$ with filling factor $n=0.4$ as
shown in Fig. 4. The polarization $m$ approaches $1$ across the phase
transition to the topologically nontrivial phase with nonzero winding number
($W=1$). In the mean field framework, the polarization is
$m=\frac{1}{2}\left[1+\tanh(\frac{\beta(h-h_{c})}{2})\right]$. Using the QMC
result, the spin polarization as shown in insert of Fig. 4, can be fitted by

$\displaystyle m=\delta
n(0,0)=f(\beta^{1/\nu}(h-h_{c})/2)=\left[1+\tanh\left(\beta^{1/\nu}(h-h_{c})/2\right)\right]/2,$
(4)

where we find critical exponent $\nu=4/3\pm 0.1$. In the limit of zero
temperature($\beta\rightarrow\infty$), the sharp change of polarization
becomes the discontinuous jump and the critical Zeeman field is $h_{c}=1.74\pm
0.08$.

## Conclusion

We have investigated the unconventional pairing and the topological phase
transition in the spin-imbalanced spin-orbit coupled Fermi gases by unbiased
QMC simulation. We study the system at an experimentally accessible
temperature. We find the SOC can enrich the properties of the pairing and
stabilize the exotic FFLO+BCS pairing. The interplay between FFLO pairing and
BCS pairing results from the change of spin texture on the Fermi surface due
to the competition between the Zeeman field and the spin-orbit coupling. We
develop the scheme to probe the topological phase transition in the QMC method
by characterizing the jump of the spin polarization at time-reversal invariant
momenta, in which the spin polarization can be measured in the ultracold-atom
experiment directly. We find that with increasing the Zeeman field,
topological phase transitions between topologically trivial superfluid phase
and nontrivial phases occur, and we give the phase diagram.

The 2D spin-imbalanced spin-orbit Fermi gas could be further explored with the
numerically exact QMC method. One interesting direction would be to develop
the method to probe the topological signature such as edge current or other
unconventional orders and to further explore the topological aspect of
FFLO+BCS phase. Moreover, it would also be beneficial to explore
unconventional pairing in the presence of in-plane magnetic fields, and phase
diagrams of other related systems. Our results lay the theoretical foundation
for the understanding of these topological superfluids, and also for the
detection of these interesting physics.

## Method

Quantum Monte Carlo method is a numerically exact method for investigating
strongly correlated systems. It is especially useful in 2D, as the system size
and dimensionality cause problems in other techniques like exact
diagonalization and density matrix renormalization group method. In this
manuscript, we use the finite temperature determinant quantum Monte-Carlo
(DQMC), whose correlation functions are given by

$\displaystyle\langle\hat{O}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{Z}Tr\big{[}\hat{O}e^{\beta
H}\big{]}=\frac{1}{Z}\int\mathcal{D}[\phi_{i,\tau}]e^{-S[\phi_{i,\tau}]}\hat{O}[\phi_{i,\tau}]\,,$
(5)

, where $\hat{O}$ is the measurable operator, $\beta$ is the inverse
temperature, and $S$ is the action with $\phi_{i,\tau}$ the set of auxiliary
fields in both spatial dimension and time dimension (please refer to
literature for details [60]).

The main step in DQMC is the matrix propagation using the propagators from
interacting part $e^{-\Delta\tau H_{int}}$ and non-interacting part
$e^{-\Delta\tau H_{non}}$ of the Hamiltonian, where $-\Delta\tau$ is the time
interval for Trotter’s decomposition. We use a Hubbard-Stratonovich
transformation to convert the interaction term into a non-interacting term
coupled to an discrete auxiliary field. This transformation enables us to
treat Hubbard models with attractive electron interactions. If the interaction
is local, $e^{-\Delta\tau H_{int}}$ is usually a diagonal matrix after the
Hubbard-Stratanovich transformation. We use the discrete Hubbard-Stratanovich
transformation, where the auxiliary field is of spin +1 or -1.

To study SOC systems using DQMC [42, 43, 44], we have made an advancement in
the QMC methods by considering a matrix representing both spin-$\uparrow$ and
spin-$\downarrow$. Without SOC, the propagating matrix of non-interacting part
is a $N\times N$ matrix, where $N$ is the number of sites. The matrix is used
to do the propagation for spin-$\uparrow$ part and spin-$\downarrow$ part
separately. With SOC, the size of the propagator in the simulation becomes
$2N\times 2N$, as the spin-flip terms destroy the block-diagonal property in
the matrix. Therefore, the computation time requires in the SOC algorithm is
around 4 times more than the case without SOC, as the majority of time is
spent in the matrix multiplication with time complexity of $(2N)^{3}$ and
$N^{3}\times 2$ in the case of with and without SOC respectively.

Negative sign problem is an unsolved problem in the QMC method [61, 62]. The
negative weight appears in some Monte Carlo samples, leading to the
exponential growth of computational time required in the case of serious sign
problem. It limits the usage of QMC method to specific types of systems. In
the case of spin-balanced SOC system with attractive interaction, it has been
shown that no sign problem exists in the system [46]. However, in the case of
spin-imbalanced system, the sign problem arises. To handle the problem, we use
the large-scale finite temperature DQMC to get the results in some regime of
parameter space. In this manuscript, we carry out the simulation in $8\times
16$ lattice with inverse temperature $\beta=4-12$, which corresponds to the
experimental setting and ensures the convergence of our result to the ground
state. Note that the temperature $\beta=4$ has already achieved in experiments
in the lattice model several years ago [48].

## Acknowledgements

H.-Q.L. acknowledges the financial support from NSAF (No. U1930402) and NSFC
(No. 11734002). X.Y. and B.H. acknowledge the support from NSFC (No.11504143,
No.11547047). M.G. acknowledges the financial support from the NYTTP (No.
KJ2030000001) and the NNSFC (No.GG2470000101).

## Author contributions statement

X.Y. and H.K.T. contributed equally to this work. X.Y., H.K.T., M.G. and
H.Q.L. conceived the project. X.Y. and H.K.T. conducted the calculations. All
authors took part in analyzing the results and wrote the manuscript.

## Data Availability

The data and code supporting the findings of this study are available from the
corresponding authors upon request.

## References

  * [1] Hasan, M. Z. & Kane, C. L.  Colloquium : Topological insulators.  _Rev. Mod. Phys._ 82, 3045–3067, DOI: 10.1103/RevModPhys.82.3045 (2010).
  * [2] Qi, X.-L. & Zhang, S.-C.  Topological insulators and superconductors.  _Rev. Mod. Phys._ 83, 1057–1110, DOI: 10.1103/RevModPhys.83.1057 (2011).
  * [3] Sato, M. & Ando, Y.  Topological superconductors: a review.  _Reports on Progress in Physics_ 80, 076501, DOI: 10.1088/1361-6633/aa6ac7 (2017).
  * [4] Qi, X.-L., Wu, Y.-S. & Zhang, S.-C.  Topological quantization of the spin hall effect in two-dimensional paramagnetic semiconductors.  _Phys. Rev. B_ 74, 085308, DOI: 10.1103/PhysRevB.74.085308 (2006).
  * [5] Zhu, S.-L., Fu, H., Wu, C.-J., Zhang, S.-C. & Duan, L.-M.  Spin hall effects for cold atoms in a light-induced gauge potential.  _Phys. Rev. Lett._ 97, 240401, DOI: 10.1103/PhysRevLett.97.240401 (2006).
  * [6] Burkov, A. A. & Balents, L.  Weyl semimetal in a topological insulator multilayer.  _Phys. Rev. Lett._ 107, 127205, DOI: 10.1103/PhysRevLett.107.127205 (2011).
  * [7] Weng, H., Fang, C., Fang, Z., Bernevig, B. A. & Dai, X.  Weyl semimetal phase in noncentrosymmetric transition-metal monophosphides.  _Phys. Rev. X_ 5, 011029, DOI: 10.1103/PhysRevX.5.011029 (2015).
  * [8] Armitage, N. P., Mele, E. J. & Vishwanath, A.  Weyl and dirac semimetals in three-dimensional solids.  _Rev. Mod. Phys._ 90, 015001, DOI: 10.1103/RevModPhys.90.015001 (2018).
  * [9] Lindner, N. H., Refael, G. & Galitski, V.  Floquet topological insulator in semiconductor quantum wells.  _Nature Physics_ 7, 490, DOI: 10.1038/nphys1926 (2011).
  * [10] Rudner, M. S., Lindner, N. H., Berg, E. & Levin, M.  Anomalous edge states and the bulk-edge correspondence for periodically driven two-dimensional systems.  _Phys. Rev. X_ 3, 031005, DOI: 10.1103/PhysRevX.3.031005 (2013).
  * [11] Xu, S.-Y. _et al._ Discovery of a weyl fermion semimetal and topological fermi arcs.  _Science_ 349, 613–617, DOI: 10.1126/science.aaa9297 (2015).
  * [12] Wang, P. _et al._ Spin-orbit coupled degenerate fermi gases.  _Phys. Rev. Lett._ 109, 095301, DOI: 10.1103/PhysRevLett.109.095301 (2012).
  * [13] Cheuk, L. W. _et al._ Spin-injection spectroscopy of a spin-orbit coupled fermi gas.  _Phys. Rev. Lett._ 109, 095302, DOI: 10.1103/PhysRevLett.109.095302 (2012).
  * [14] Lin, Y.-J. _et al._ Bose-einstein condensate in a uniform light-induced vector potential.  _Phys. Rev. Lett._ 102, 130401, DOI: 10.1103/PhysRevLett.102.130401 (2009).
  * [15] Lin, Y.-J., Compton, R. L., Jimenez-Garcia, K., Porto, J. V. & Spielman, I. B.  Synthetic magnetic fields for ultracold neutral atoms.  _Nature_ 462, 628–632, DOI: 10.1038/nature08609 (2009).
  * [16] Goldman, N. _et al._ Non-abelian optical lattices: Anomalous quantum hall effect and dirac fermions.  _Phys. Rev. Lett._ 103, 035301, DOI: 10.1103/PhysRevLett.103.035301 (2009).
  * [17] Yi, W. & Guo, G.-C.  Phase separation in a polarized fermi gas with spin-orbit coupling.  _Phys. Rev. A_ 84, 031608, DOI: 10.1103/PhysRevA.84.031608 (2011).
  * [18] Ozawa, T., Pitaevskii, L. P. & Stringari, S.  Supercurrent and dynamical instability of spin-orbit-coupled ultracold bose gases.  _Phys. Rev. A_ 87, 063610, DOI: 10.1103/PhysRevA.87.063610 (2013).
  * [19] Banerjee, S., Erten, O. & Randeria, M.  Ferromagnetic exchange, spin-orbit coupling and spiral magnetism at the laalo_3/srtio_3 interface.  _Nat. Phys._ 9, 626–630, DOI: 10.1038/nphys2702 (2013).
  * [20] Bert, J. A. _et al._ Direct imaging of the coexistence of ferromagnetism and superconductivity at the laalo3/srtio3 interface.  _Nat. Phys._ 7, 767–771, DOI: 10.1038/nphys2079 (2011).
  * [21] Cole, W. S., Zhang, S., Paramekanti, A. & Trivedi, N.  Bose-hubbard models with synthetic spin-orbit coupling: Mott insulators, spin textures, and superfluidity.  _Phys. Rev. Lett._ 109, 085302, DOI: 10.1103/PhysRevLett.109.085302 (2012).
  * [22] Li, Y., Pitaevskii, L. P. & Stringari, S.  Quantum tricriticality and phase transitions in spin-orbit coupled bose-einstein condensates.  _Phys. Rev. Lett._ 108, 225301, DOI: 10.1103/PhysRevLett.108.225301 (2012).
  * [23] Wang, C., Gao, C., Jian, C.-M. & Zhai, H.  Spin-orbit coupled spinor bose-einstein condensates.  _Phys. Rev. Lett._ 105, 160403, DOI: 10.1103/PhysRevLett.105.160403 (2010).
  * [24] Yu, Z.-Q. & Zhai, H.  Spin-orbit coupled fermi gases across a feshbach resonance.  _Phys. Rev. Lett._ 107, 195305, DOI: 10.1103/PhysRevLett.107.195305 (2011).
  * [25] Zhai, H.  Spin-orbit coupled quantum gases.  _Int. J. Mod. Phys. B_ 26, 1230001, DOI: 10.1142/S0217979212300010 (2012).
  * [26] Mitra, D. _et al._ Quantum gas microscopy of an attractive Fermi–Hubbard system.  _Nat. Phys._ 14, 173–177, DOI: 10.1038/nature08482 (2018).
  * [27] Fulde, P. & Ferrell, R. A.  Superconductivity in a strong spin-exchange field.  _Phys. Rev._ 135, A550–A563, DOI: 10.1103/PhysRev.135.A550 (1964).
  * [28] Larkin, A. I. & Ovchinnikov, Y. N.  Inhomogeneous state of superconductors.  _Sov. Phys. JETP_ 20, 762 (1965).
  * [29] Zheng, Z. _et al._ FFLO Superfluids in 2D Spin-Orbit Coupled Fermi Gases.  _Sci. Rep._ 4, 6535, DOI: 10.1038/srep06535 (2014).
  * [30] Huang, B., Yang, X., Xu, N., Zhou, J. & Gong, M.  Dynamical instability with respect to finite-momentum pairing in quenched bcs superconducting phases.  _Phys. Rev. B_ 99, 014517, DOI: 10.1103/PhysRevB.99.014517 (2019).
  * [31] Qu, C. _et al._ Topological superfluids with finite-momentum pairing and majorana fermions.  _Nat. Commun._ 4, 2710, DOI: 10.1038/ncomms3710 (2013).
  * [32] Sato, M., Takahashi, Y. & Fujimoto, S.  Non-abelian topological order in $s$-wave superfluids of ultracold fermionic atoms.  _Phys. Rev. Lett._ 103, 020401, DOI: 10.1103/PhysRevLett.103.020401 (2009).
  * [33] Gong, M., Tewari, S. & Zhang, C.  Bcs-bec crossover and topological phase transition in 3d spin-orbit coupled degenerate fermi gases.  _Phys. Rev. Lett._ 107, 195303, DOI: 10.1103/PhysRevLett.107.195303 (2011).
  * [34] Gong, M., Chen, G., Jia, S. & Zhang, C.  Searching for majorana fermions in 2d spin-orbit coupled fermi superfluids at finite temperature.  _Phys. Rev. Lett._ 109, 105302, DOI: 10.1103/PhysRevLett.109.105302 (2012).
  * [35] Dong, Y., Dong, L., Gong, M. & Pu, H.  Dynamical phases in quenched spin-orbit-coupled degenerate Fermi gas.  _Nature Communications_ 6, 6103, DOI: 10.1038/ncomms7103 (2015).
  * [36] Huang, B., Chan, C. F. & Gong, M.  Large chern-number topological superfluids in a coupled-layer system.  _Phys. Rev. B_ 91, 134512, DOI: 10.1103/PhysRevB.91.134512 (2015).
  * [37] Zwicknagl, G., Jahns, S. & Fulde, P.  Critical magnetic field of ultra-thin superconducting films and interfaces.  _Journal of the Physical Society of Japan_ 86, 083701, DOI: 10.7566/JPSJ.86.083701 (2017).
  * [38] Yuan, N. F. Q. & Fu, L.  Topological metals and finite-momentum superconductors.  _Proc. Natl. Acad. Sci. U.S.A._ 118, DOI: 10.1073/pnas.2019063118 (2021).
  * [39] Liang, J. _et al._ Unconventional pairings of spin-orbit coupled attractive degenerate fermi gas in a one-dimensional optical lattice.  _Sci. Rep._ 5, 14863, DOI: 10.1038/srep14863 (2015).
  * [40] Roy, M. S. & Kumar, M.  Fulde-ferrel-larkin-ovchinnikov phase in one dimensional fermi gas with attractive interactions and transverse spin-orbit coupling (2021).  ArXiv:2108.03314.
  * [41] Becca, F. & Sorella, S.  _Quantum Monte Carlo approaches for correlated systems_ (Cambridge University Press, 2017).
  * [42] Tang, H. K., Yang, X., Sun, J. & Lin, H. Q.  Berezinskii-kosterlitz-thoules phase transition of spin-orbit coupled fermi gas in optical lattice.  _EPL_ 107, 40003, DOI: 10.1209/0295-5075/107/40003 (2014).
  * [43] Shi, H., Rosenberg, P., Chiesa, S. & Zhang, S.  Rashba spin-orbit coupling, strong interactions, and the bcs-bec crossover in the ground state of the two-dimensional fermi gas.  _Phys. Rev. Lett._ 117, 040401, DOI: 10.1103/PhysRevLett.117.040401 (2016).
  * [44] Rosenberg, P., Shi, H. & Zhang, S.  Ultracold atoms in a square lattice with spin-orbit coupling: Charge order, superfluidity, and topological signatures.  _Phys. Rev. Lett._ 119, 265301, DOI: 10.1103/PhysRevLett.119.265301 (2017).
  * [45] Rosenberg, P., Shi, H. & Zhang, S.  Accurate computations of rashba spin-orbit coupling in interacting systems: From the fermi gas to real materials.  _Journal of Physics and Chemistry of Solids_ 128, 161 – 168, DOI: 10.1016/j.jpcs.2017.12.026 (2019).
  * [46] Wu, C. & Zhang, S.-C.  Sufficient condition for absence of the sign problem in the fermionic quantum monte carlo algorithm.  _Phys. Rev. B_ 71, 155115, DOI: 10.1103/PhysRevB.71.155115 (2005).
  * [47] Wolak, M. J., Grémaud, B., Scalettar, R. T. & Batrouni, G. G.  Pairing in a two-dimensional fermi gas with population imbalance.  _Phys. Rev. A_ 86, 023630, DOI: 10.1103/PhysRevA.86.023630 (2012).
  * [48] mazurenko, A. _et al._ A cold-atom fermi–hubbard antiferromagnet.  _Nature_ 545, 462, DOI: 10.1038/nature22362 (2017).
  * [49] Gell, M., Wurz1, N., Samland, J., Chan, C. F. & Höhl, M.  Competing magnetic orders in a bilayer hubbard model with ultracold atoms.  _Nature_ 589, 40, DOI: 10.1038/s41586-020-03058-x (2020).
  * [50] Feiguin, A. E. & Heidrich-Meisner, F.  Pairing states of a polarized fermi gas trapped in a one-dimensional optical lattice.  _Phys. Rev. B_ 76, 220508, DOI: 10.1103/PhysRevB.76.220508 (2007).
  * [51] Rizzi, M. _et al._ Fulde-Ferrell-Larkin-Ovchinnikov pairing in one-dimensional optical lattices.  _Phys. Rev. B_ 77, 245105, DOI: 10.1103/PhysRevB.77.245105 (2008).
  * [52] Tezuka, M. & Ueda, M.  Density-matrix renormalization group study of trapped imbalanced fermi condensates.  _Phys. Rev. Lett._ 100, 110403, DOI: 10.1103/PhysRevLett.100.110403 (2008).
  * [53] Wolak, M. J., Grémaud, B., Scalettar, R. T. & Batrouni, G. G.  Pairing in a two-dimensional fermi gas with population imbalance.  _Phys. Rev. A_ 86, 023630, DOI: 10.1103/PhysRevA.86.023630 (2012).
  * [54] Volovik, G. E.  _The Universe in a Helium Droplet_ (Oxford University Press,USA, 2003).
  * [55] Wang, Z. & Zhang, S.-C.  Simplified topological invariants for interacting insulators.  _Phys. Rev. X_ 2, 031008, DOI: 10.1103/PhysRevX.2.031008 (2012).
  * [56] Wang, Z. & Zhang, S.-C.  Strongly correlated topological superconductors and topological phase transitions via green’s function.  _Phys. Rev. B_ 86, 165116, DOI: 10.1103/PhysRevB.86.165116 (2012).
  * [57] Luo, X., Zhou, K., Liu, W., Liang, Z. & Zhang, Z.  Fidelity susceptibility and topological phase transition of a two-dimensional spin-orbit-coupled fermi superfluid.  _Phys. Rev. A_ 89, 043612, DOI: 10.1103/PhysRevA.89.043612 (2014).
  * [58] Zhang, W. & Yi, W.  Topological fulde–ferrell–larkin–ovchinnikov states in spin–orbit-coupled fermi gases.  _Nat. Commun._ 4, 2711, DOI: 10.1038/ncomms3711 (2013).
  * [59] See Supplemental Material for details of the mean field studies.
  * [60] Santos, R. R. d.  Introduction to quantum monte carlo simulations for fermionic systems.  _Brazilian Journal of Physics_ 33, 36–54, DOI: 10.1590/S0103-97332003000100003 (2003).
  * [61] Troyer, M. & Wiese, U.-J.  Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations.  _Phys. Rev. Lett._ 94, 170201, DOI: 10.1103/PhysRevLett.94.170201 (2005).
  * [62] Li, Z.-X. & Yao, H.  Sign-problem-free fermionic quantum monte carlo: Developments and applications.  _Annual Review of Condensed Matter Physics_ 10, 337–356, DOI: 10.1146/annurev-conmatphys-033117-054307 (2019).


# Optimal selection of the number of control units in kNN algorithm to
estimate average treatment effects

Andrés Ramírez-Hassan Department of Economics, Universidad EAFIT, Medellín,
Colombia; email<EMAIL_ADDRESS>Raquel Vargas-Correa Department of
Economics, Universidad EAFIT, Medellín, Colombia; email<EMAIL_ADDRESS>Gustavo Garcia Department of Economics, Universidad EAFIT, Medellín, Colombia;
email<EMAIL_ADDRESS>Daniel Londoño Department of Economics,
Universidad EAFIT, Medellín, Colombia; email<EMAIL_ADDRESS>
###### Abstract

We propose a simple approach to optimally select the number of control units
in $k$ nearest neighbors (kNN) algorithm focusing in minimizing the mean
squared error for the average treatment effects. Our approach is non-
parametric where confidence intervals for the treatment effects were
calculated using asymptotic results with bias correction. Simulation exercises
show that our approach gets relative small mean squared errors, and a balance
between confidence intervals length and type I error. We analyzed the average
treatment effects on treated (ATET) of participation in 401(k) plans on
accumulated net financial assets confirming significant effects on amount and
positive probability of net asset. Our optimal $k$ selection produces
significant narrower ATET confidence intervals compared with common practice
of using $k=1$.

JEL Classification: C14, C18, O18, R42.
Keywords: Average treatment effects, k nearest neighbors, optimal number of
control units, 401(k).

##  1 Introduction

Matching estimators are commonly used to estimate average treatment effects
when treatment assignment is independent of outcome conditional on confounding
variables (Cameron and Trivedi,, 2005, p. 871). However, it seems that optimal
selection of the number of matching units ($k$) is an open question as most of
the applied literature fixes $k$ small without any optimal criterion, often
$k=1$ (Abadie and Imbens,, 2016). Choice of $k$ is a trade-off between bias
and variance. Small $k$ implies small bias, but high variance, and vice versa.
So, the aim of this paper is to propose a simple approach to optimally select
$k$. In particular, we focus on the $k$ nearest-neighbors (kNN) algorithm when
used to estimate treatment effects: average treatment effect (ATE), and
average treatment effect on treated (ATET).

Optimal selection of $k$ in kNN algorithm is done by cross validation when
prediction or classification are the objectives (Hastie et al.,, 2009, p.
241). This is easily performed as the statistical object of interest is
observable in-sample. On the other hand, average treatment effects are not
observable anywhere. So, we follow similar arguments to Athey et al., (2015)
proposing an unbiased function of observable variables to optimally choose $k$
through cross validation when the main inferential concern is treatment
effects.

We apply our proposal to analyze ATET due to enrollment in 401(k) on amount
and probability of positive accumulated net financial assets. We found
significant ATET whose point estimates are approximately $15 K and 19% given
an optimal $k^{*}=19$ in both cases.

After this brief introduction, Section 2 shows our proposal to optimally
select $k$ in the kNN algorithm. Section 3 displays results from some
simulations exercises, and Section 4 shows the results of our application
regarding 401(k) enrollment. Section 5 ends with some concluding remarks.

##  2 Optimal selection of $k$ in kNN

We will focus on the average treatment effect on treated as we will use this
in our application. However, results associated with the average treatment
effects are in the Appendix (see subsection 6.1). We build the matching units
using $k$ nearest neighbors (kNN) with replacement, where $k$ is selected
minimizing the mean squared error between a conditional unbiased estimator of
the average treatment effect on the treated (ATET) and a matching estimator of
the ATET.

We consider a framework where there is a random binary treatment
($D_{i}=\left\\{0,1\right\\}$) that is independent of the outcome variable
($Y_{i}$) conditional on independent observable variables
($\bm{X}_{i}\subset\mathbb{R}^{P}$), $i=1,2,\dots,N$. The average treatment
effect on treated,

$\Delta^{T}=\mathbb{E}\left[(Y_{i}(1)-Y_{i}(0))|D_{i}=1,\bm{X}_{i}=\bm{x}_{i}\right],$
(1)

where $\bm{X}_{i}$ are control variables, $\left\\{Y_{i}(0),Y_{i}(1)\right\\}$
are the potential outcomes under different states, untreated and treated,

$\displaystyle Y_{i}=Y_{i}(D_{i})=\begin{Bmatrix}Y_{i}(0),&D_{i}=0\\\
Y_{i}(1),&D_{i}=1\end{Bmatrix}.$

However, we do not simultaneously observe the same individual under both
states. So, we propose to use the matching approach to build the synthetic
controls, the general formula is (Cameron and Trivedi,, 2005, p. 875),

$\hat{\Delta}^{TM}=\frac{1}{N_{1}}\sum_{i\in\left\\{D_{i}=1\right\\}}\left[Y_{i}(1)-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x})}Y_{j}(0)\right],$ (2)

where
$A_{k}^{i}(\bm{x})=\left\\{j:||\bm{X}_{i}-\bm{X}_{j}||<||\bm{X}_{i}-\bm{X}_{l}||,j\in\left\\{D_{j}=0\right\\},l\in\left\\{D_{l}=0\right\\}\right\\}$
is the set of the $k$ closest untreated units to treated unit $i$ in terms of
the Euclidean norm in covariates, $j,l=1,2,\dots,N_{0}$, $||\cdot||$ is the
Euclidean norm, $N_{0}$ is the sample size of untreated individuals
($D_{i}=0$), $N_{1}$ is the sample size of treated individuals ($D_{i}=1$),
and $k$ is the cardinality of $A_{k}^{i}(\bm{x})$, that is, the number of
controls for treated unit,

$|A_{k}^{i}(\bm{x})|=\sum_{j\in
A_{k}^{i}(\bm{x})}\mathbbm{1}\left\\{||\bm{X}_{i}-\bm{X}_{j}||<||\bm{X}_{i}-\bm{X}_{l}||\right\\}=k.$

Therefore, the synthetic control for treated individual $i$ is built using an
average of untreated nearest individuals (neighbors).

The choice of $k$ is a trade-off between bias and variance. Selecting just one
neighbor minimizes bias, this is equal to zero if there is exact match
($k=1,\bm{x_{i}}=\bm{x_{j}}$), but implies high variability. On the other
hand, a large amount of neighbors decreases variance but increases bias. In
general, $\bm{x_{i}}=\bm{x_{j}}$ is an event of probability 0 for continuous
covariates, then Abadie and Imbens, (2006) show that matching estimators have
an asymptotic bias. This asymptotic bias can be ignored if
$N_{1}^{P_{c}/2}=O(N_{0})$ where $P_{c}$ is the number of continuous
covariates ($P_{c}/2>1$). The bias of the ATET is

$B^{TM}=\frac{1}{N_{1}}\sum_{i\in\left\\{D_{i}=1\right\\}}\left[\mu_{0}(\bm{x}_{i})-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x})}\mu_{0}(\bm{x}_{j})\right],$ (3)

where $\mu_{0}(\bm{x}_{i})=\mathbb{E}\left[Y|\bm{X}=\bm{x},D=0\right]$. This
can be estimated non-parametrically using a series expansion estimator to
obtain $\hat{B}^{TM}$ (Abadie and Imbens,, 2011).

Therefore by equations 2 and 3, the bias-corrected matching estimator is

$\hat{\Delta}^{BC}=\hat{\Delta}^{TM}-\hat{B}^{TM}.$ (4)

Abadie and Imbens, (2006) show under suitable assumptions (see Assumptions 1,
2’, 3’ and 4 in their paper) that

$(V^{E,T}+V^{\Delta(\bm{X}),T})^{-1/2}\sqrt{N_{1}}(\hat{\Delta}^{TM}-B^{TM}-\Delta^{T})\xrightarrow[]{d}N(0,1),$

where
$V^{E,T}=\frac{1}{N_{1}}\sum_{i=1}^{N}\left(D_{i}+(1-D_{i})\frac{J_{k}^{i}}{k}\right)^{2}\sigma^{2}(\bm{X}_{i},D_{i})$,
$J_{k}^{i}=\sum_{j=1}^{N}\mathbbm{1}\left\\{i\in A_{k}^{j}\right\\}$ is the
number of times unit $i$ is used as a match given $k$ matches per unit,
$\sigma^{2}(\bm{X}_{i},D_{i})=\mathbb{V}(Y|\bm{X}_{i}=\bm{x}_{i},D_{i}=d_{i})$,
which can be consistently estimated by
$\hat{\sigma}^{2}(\bm{X}_{i},D_{i})=\frac{k}{k+1}\left(Y_{i}-\frac{1}{k}\sum_{l\in
L_{k}^{i}}\right)$, $L_{k}^{i}(\bm{x})$ is the set of the $k$ closest units to
unit $i$ in terms of the Euclidean norm in covariates such that have the same
value for the treatment variable, and
$V^{\Delta(\bm{X}),T}=\mathbb{E}\left[(\left[\mu_{1}(\bm{X}_{i})-\mu_{0}(\bm{X}_{i})\right]-{\Delta}^{T})^{2}|\bm{D}=1\right]$
is the variance of the conditional average treatment effect on treated.

In an ideal situation, we select $k$ in the regression nearest neighbor
framework (kNN) minimizing the mean squared error of the average treatment
effect on treated (equation 1), that is,

$\operatorname*{arg\,min}_{k}\mathbb{E}\left({\hat{\Delta}^{BC}}-{\Delta}^{T}\right)^{2}.$
(5)

However, the average treatment effect on treated ($\Delta^{T}$) is not
observed. Therefore, we follow similar arguments to Athey et al., (2015)
proposing an unbiased function of potentially observable variables to define
the ATET in program 5.

We have the following assumptions:

######  Assumption 1.

Stable unit treatment value assumption (SUTVA)

######  Assumption 2.

Ignorability (unconfoundedness)

$D_{i}\perp Y_{i}(0)|\bm{X}_{i}=\bm{x}_{i}.$

######  Assumption 3.

Overlap (matching)

$0\leq P(D_{i}=1|\bm{X}_{i}=\bm{x}_{i})<1.$

######  Remark 1.

Assumption 1 implies not spillover, interaction or general equilibrium effects
(Rubin,, 1978). Assumption 2 establishes that treatment assignment ignores
untreated outcome given control variables (Rubin,, 1978; Rosenbaum and Rubin,,
1983). Assumption 3 says that it is necessary to have overlap in subsamples of
$\bm{X}_{i}$, that is, for each treated individual it is necessary to have an
analogous untreated (synthetic control) to identify the treatment effect on
treated (Rosenbaum and Rubin,, 1983).

We assume conditional independent and identically distributed random sample
drawn from a population. This implicitly implies Assumption 1 (SUTVA). Then,
we define the individual treatment effect on treated

${\Delta_{i}^{T*}}=\frac{Y_{i}(D_{i}-P(\bm{X}_{i}))}{P\left[D=1\right](1-P(\bm{X}_{i}))},$
(6)

where $P\left[D=1\right]$ is the marginal probability of treatment, which is
the same for every individual given random sampling, and
$P(\bm{X}_{i})=P\left[D_{i}=1|\bm{X}_{i}=\bm{x}\right]$ is the conditional
probability of treatment given regressors $\bm{X}_{i}$, that is, the
propensity score (Rosenbaum and Rubin,, 1983).

######  Proposition 1.

Assuming a conditional independent and identically distributed random sample
such that Assumptions 2 and 3 are satisfied, then:

$\mathbb{E}\left[{\Delta_{i}^{T*}}\big{|}P(\bm{X}_{i}),\bm{X}_{i}=\bm{x}\right]=\Delta^{T}.$

###### Proof.

We omit conditioning on $P\left[D=1\right]$ and $\bm{X}_{i}=\bm{x}$ to
simplify notation.

$\displaystyle\mathbb{E}\left[{\Delta_{i}^{T*}}\big{|}P\left[D=1\right],\bm{X}_{i}=\bm{x}\right]=$
$\displaystyle\mathbb{E}\left[\frac{Y_{i}(D_{i}-P(\bm{X}_{i}))}{P\left[D=1\right](1-P(\bm{X}_{i}))}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{P\left[D=1\right](1-P(\bm{X}_{i}))}\left\\{(1-P(\bm{X}_{i}))\mathbb{E}\left[D_{i}(Y_{i}(1)-Y_{i}(0))\right]\right.$
$\displaystyle\left.+\mathbb{E}\left[Y_{i}(0)(D_{i}-P(\bm{X}_{i}))\right]\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{P\left[D=1\right](1-P(\bm{X}_{i}))}\left\\{(1-P(\bm{X}_{i}))\mathbb{E}\left[\mathbb{E}\left[D_{i}(Y_{i}(1)-Y_{i}(0))\big{|}D_{i}\right]\right]\right.$
$\displaystyle\left.+\mathbb{E}\left[D_{i}\right]\mathbb{E}\left[Y_{i}(0)\right]-P(\bm{X}_{i})\mathbb{E}\left[Y_{i}(0)\right]\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{P\left[D=1\right](1-P(\bm{X}_{i}))}\left\\{(1-P(\bm{X}_{i}))P\left[D=1\right]\mathbb{E}\left[(Y_{i}(1)-Y_{i}(0))\big{|}D_{i}=1\right]\right.$
$\displaystyle\left.+\mathbb{E}\left[Y_{i}(0)\right](\mathbb{E}\left[D_{i}\right]-P(\bm{X}_{i}))\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[(Y_{i}(1)-Y_{i}(0))\big{|}D_{i}=1\right]$

The first equality is by definition, the second is after some simply algebra
and taking into account that $Y_{i}=D_{i}Y_{i}(1)+(1-D_{i})Y_{i}(0)$ and
$D_{i}Y_{i}=D_{i}Y_{i}(1)+D_{i}Y_{i}(0)-D_{i}Y_{i}(0)$ given
$D_{i}=\left\\{0,1\right\\}$. The third is an application of the law of
iterative expectations for the first term, and the Assumption 2 (ignorability)
for the second term. The fourth equality takes again into account that
$D_{i}=\left\\{0,1\right\\}$, and $P\left[D_{i}=1\right]=P\left[D=1\right]$
(random sampling). Finally, we take into account that
$\mathbb{E}\left[D_{i}|\bm{X}_{i}=\bm{x}_{i}\right]=P\left[D_{i}=1|\bm{X}_{i}=\bm{x}\right]=P(\bm{X}_{i})$.
Notice that Assumption 3 is required to have
$P\left[D=1\right](1-P(\bm{X}_{i}))\neq 0$.

######  Remark 2.

Observe that
$\mathbb{E}\left[Y_{i}(0)D_{i}\right]=P\left[D=1\right]\mathbb{E}\left[Y_{i}(0)|D_{i}=1\right]$.
Therefore, we also obtain the required statement, if

$\mathbb{E}\left[Y_{i}(0)|D_{i}=1\right]=\mathbb{E}\left[Y_{i}(0)|D_{i}=0\right]=\mathbb{E}\left[Y_{i}(0)\right]$

conditional on $\bm{X}_{i}=\bm{x}_{i}$ is satisfied, that is, the conditional
mean assumption on control group is satisfied (Angrist and Krueger,, 1999, p.
1316). Outcome of untreated individuals does not determine participation.

∎

We require estimators for $P\left[D=1\right]$ and $P(\bm{X}_{i})$ to estimate
$\Delta_{i}^{T*}$. We use $\frac{1}{N}\sum_{i=1}^{N}D_{i}$, $N=N_{0}+N_{1}$,
as estimator for the marginal probability of treatment, and a logit model to
estimate the conditional probability of treatment. Obviously, these estimators
can be changed. Notice that this introduces an extra source of variability as
we do not observe $P\left[D=1\right]$ and $P(\bm{X}_{i})$ in observational
data sets. This kind of issue is discussed by Athey et al., (2015) who point
out that transformations like equation 7 does not optimally use all available
information, this may mean extra variability. Additionally, Abadie and Imbens,
(2016) show that ATE using matching based on the estimated propensity score is
more efficient than matching based on the true propensity score. On the other
hand, ATET using the estimated propensity score can be more o less efficient
compared to the true propensity score. We performed some simulation exercises
that show better inferential performance based on the estimated propensity
score (available upon authors request). In addition, as our ATET (ATE) is
asymptotically biased corrected, the mean squared error decomposition shows
that variance explains the most part (75% on average), as expected. However,
it remains a proportion to finite sample bias effect (25% on average).

Give that ${\Delta_{i}^{T*}}$ in equation 8 is a conditional unbiased
estimator for $\Delta^{T}$ (Proposition 1), we select $k$ in the kNN
regression framework such that

$\operatorname*{arg\,min}_{{k}}\frac{1}{G}\sum_{g=1}^{G}\left\\{\left(\hat{\Delta}_{g}^{BC}-\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}\widehat{\Delta_{i,g}^{T*}}\right)^{2}\right\\},$
(7)

where

$\widehat{\Delta_{i,g}^{T*}}=\left(\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}D_{i,g}^{Test}\right)^{-1}\left(\frac{Y_{i,g}^{Test}(D_{i,g}^{Test}-\widehat{P(\bm{X}_{i,g}^{Test})})}{1-\widehat{P(\bm{X}_{i,g}^{Test}})}\right),$

$\hat{\Delta}_{g}^{TM}=\frac{1}{N_{1,g}}\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\left[Y_{i,g}^{Test}(1)-\frac{1}{K}\sum_{j\in
A_{K}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)\right],$

$\hat{B}_{g}^{TM}=\frac{1}{N_{1,g}}\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\left[\hat{\mu}_{0}(\bm{x}_{i})_{g}^{Test}-\frac{1}{K}\sum_{j\in
A_{K}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{0}(\bm{x}_{j})_{g}^{Train}\right],$

and

$\hat{\Delta}_{g}^{BC}=\hat{\Delta}_{g}^{TM}-\hat{B}_{g}^{TM}$

where $G$ is the number of groups, $N_{g}$ is the sample size of each group,
and $N_{1,g}$ is the sample size of the treated group in group $g$,
$g=1,2,\dots,G$. Superscripts $Train$ and $Test$ refer to train and test data
sets (see below).

Notice that $\widehat{\Delta_{i,g}^{T*}}$ is the sample version of the
conditional unbiased estimator of the ATET $\Delta_{i}^{T*}$, so its average
is a consistent estimator for ATET. On the other hand, $\hat{\Delta}_{g}^{BC}$
is bias-corrected version of the matching estimator (Abadie and Imbens,,
2006).

kNN regression is a standard approach in the machine and statistical learning
communities (Hastie et al.,, 2009, p. 14). So, we follow $k$-fold cross-
validation to select the optimal $k$ (Hastie et al.,, 2009, p. 241). In
particular, we randomly split the data set in $G$ roughly equal-sized groups
such that we keep the proportion between treated and untreated individuals in
each group, then we calculate
$\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}\widehat{\Delta_{i,g}^{T*}}$,
$g=1,2,\dots,G$. We fit the model using $G-1$ groups, that is, we fix $k$, and
identify the $k$ nearest neighbors in the control group at the training data
set in terms of the Euclidean norm to individual $i$ in the treated group at
the $G$-th left group (test data set), that is,

$|A_{k}^{i}(\bm{x}_{g}^{Test})|=\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\mathbbm{1}\left\\{||\bm{X}_{i,g}^{Test}-\bm{X}_{j}^{Train}||<||\bm{X}_{i,g}^{Test}-\bm{X}_{l}^{Train}||\right\\}=k,$

then we estimate $\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)$, and obtain
$\hat{\Delta}_{i,g}^{TM}=Y_{i,g}^{Test}(1)-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)$ and
$\hat{\Delta}_{i,g}^{BC}=\hat{\mu}_{0}(\bm{x}_{i})_{g}^{Test}-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{0}(\bm{x}_{j})_{g}^{Test}$ for
$i\in\left\\{D_{i,g}^{Test}=1\right\\}$. So, we can calculate the mean squared
error for the $G$-th left group. This is done for $g=1,2,\dots,G$ and
different values of $k$ such that we average over the $G$ prediction errors
given $k$, and select the $k$ which has the minimum average error (see program
7).

Given an optimal $k=k^{*}$, we use equation 4 to estimate the average
treatment effect on treated using the whole sample size. As we said, this
estimator converges in distribution to a normal distribution, so we use the
asymptotic result to build confidence intervals.

Algorithm A1 summarizes our methodological proposal:

Algorithm A1 Optimal $k$ in kNN to estimate ATET

1:Calculate $\frac{1}{N}\sum_{i=1}^{N}D_{i}$.

2:for $i=1,2,\dots,N$ do

  * •

Fit a logit model where the dependent variable is $D_{i}$ to get
$\widehat{P(\bm{X}_{i})}$.

  * •

Calculate the $\widehat{\Delta_{i}^{T*}}$, the sample version of equation 8.

3:end for

4:Split randomly the data set in $G$ roughly equal-sized groups such that the
proportion between treated and untreated in each group is keep the same as in
the original data set.

5:Let us set the training data set as $G-1$ groups, and the test data set the
$g$-th left group.

6:for $k=1,2,\dots,K$ do

7: for $g=1,2,\dots,G$ do

  * •

Calculate $\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}\widehat{\Delta_{i,g}^{T*}}$,
where $N_{g}$ is the sample size of $g$-th group (test data set).

  * •

Find the $k$ nearest neighbors of $\bm{X}_{i,g}$ for treated individuals of
the $g$th group ($D_{i,g}=1$ in the test data set) among the untreated
individuals in the training data set ($D_{j,h\notin g}=0$ in the remaining
$G-1$ groups).

  * •

Calculate $\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)$ and $\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{0}(\bm{x}_{j})_{g}^{Train}$.

  * •

Obtain $\hat{\Delta}_{i,g}^{TM}=Y_{i,g}^{Test}(1)-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)$ and
$\hat{B}_{i,g}^{TM}=\hat{\mu}_{0}(\bm{x}_{i})_{g}^{Test}-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{0}(\bm{x}_{j})_{g}^{Train}$ for
$i\in\left\\{D_{i,g}^{Test}=1\right\\}$.

  * •

Estimate
$\hat{\Delta}_{g}^{TM}=\frac{1}{N_{1,g}}\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\hat{\Delta}_{i,g}^{TM}$,
$\hat{B}_{g}^{TM}=\frac{1}{N_{1,g}}\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\hat{B}_{i,g}^{TM}$
and $\hat{\Delta}_{g}^{BC}=\hat{\Delta}_{g}^{TM}-\hat{B}_{g}^{TM}$

8: end for  • Calculate the mean squared error:
$\frac{1}{G}\sum_{g=1}^{G}\left\\{\left(\hat{\Delta}_{g}^{BC}-\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}\widehat{\Delta_{i,g}^{T*}}\right)^{2}\right\\}$

9:end for

10:Select $k$ that minimizes the mean squared error ($k^{*}$).

11:for $i=1,2,\dots,N_{1}$ do

  * •

Given $k^{*}$ and the whole data set, find the $k^{*}$ nearest neighbors in
the untreated group to each individual in the treated group.

  * •

Calculate $\hat{\Delta}_{i}^{BC}=Y_{i}(1)-\frac{1}{k^{*}}\sum_{j\in
A_{k^{*}}^{i}(\bm{x})}Y_{j}(0)$

12:end for

13:Use $\hat{\Delta}^{BC}$ as a measure of central tendency of the ATET at
neighborhood level.

14:Use theorems 4 and 7 in Abadie and Imbens, (2006) to build confidence
intervals of the ATET.

##  3 A Monte Carlos study

We set exactly the same simulation setting as Otsu and Rai, (2017) who
proposed a weighted bootstrap to perform inference of matching estimators for
ATE and ATET as it is well known that naive bootstrap produces invalid
inference due to failing to reproduce the distribution of the number of times
each unit is used as a match (Abadie and Imbens,, 2008).

We show the simulation setting for exposition purposes:

$\displaystyle Y_{i}(1)=$
$\displaystyle\tau+m_{j}(||\bm{X}_{i}||)+\epsilon_{i},$ $\displaystyle
Y_{i}(0)=$ $\displaystyle m_{j}(||\bm{X}_{i}||)+\epsilon_{i},$ $\displaystyle
D_{i}=$ $\displaystyle\mathbbm{1}\left\\{P(\bm{X}_{i})\geq v_{i}\right\\},\
v_{i}\sim U\left[0,1\right],$ $\displaystyle P(\bm{X}_{i})=$
$\displaystyle\gamma_{1}+\gamma_{2}||\bm{X}_{i}||,\
X_{i}=\left[X_{i1},\dots,X_{iP}\right]^{\top},$ $\displaystyle X_{ij}=$
$\displaystyle\psi_{i}|\zeta_{ij}|/||\zeta_{i}||,\ j=1,2,\dots,P,$
$\displaystyle\psi_{i}\sim U\left[0,1\right],\ \zeta_{i}\sim
N(\bm{0},\bm{I}_{P}),\ \epsilon_{i}\sim N(0,0.2^{2}),$

where $\epsilon_{i}$, $v_{i}$, $\psi_{i}$ and $\zeta_{i}$ are mutually
independent. The ATET as well as ATE is $\tau=0.5$, $\gamma_{1}=0.15$,
$\gamma=0.7$, $P=5$ and there are six different curves for
$m_{j}(||\bm{X}_{i}||)$ (see Otsu and Rai, (2017) for details):

$\displaystyle m_{1}(||\bm{X}_{i}||)=$ $\displaystyle 0.15+0.7||\bm{X}_{i}||,$
$\displaystyle m_{2}(||\bm{X}_{i}||)=$ $\displaystyle
0.1+0.5||\bm{X}_{i}||+0.5\exp\left(-200(||\bm{X}_{i}||-0.7)^{2}\right),$
$\displaystyle m_{3}(||\bm{X}_{i}||)=$ $\displaystyle
0.8-2(||\bm{X}_{i}||-0.9)^{2}-5(||\bm{X}_{i}||-0.7)^{3}-10(||\bm{X}_{i}||-0.6)^{10},$
$\displaystyle m_{4}(||\bm{X}_{i}||)=$ $\displaystyle
0.2+(1-||\bm{X}_{i}||)^{0.5}-0.6(0.9-||\bm{X}_{i}||)^{2},$ $\displaystyle
m_{5}(||\bm{X}_{i}||)=$ $\displaystyle
0.2+(1-||\bm{X}_{i}||)^{0.5}-0.6(0.9-||\bm{X}_{i}||)^{2}-0.1||\bm{X}_{i}||\cos(30||\bm{X}_{i}||),$
$\displaystyle m_{6}(||\bm{X}_{i}||)=$ $\displaystyle
0.4+0.25\sin(8||\bm{X}_{i}||-5)+0.4\exp\left(-16(4||\bm{X}_{i}||-2.5)^{2}\right).$

In addition, we perform other simulation exercises where the outcome variable
is binary. In particular,

$\displaystyle Y_{i}(1)$ $\displaystyle=\begin{Bmatrix}1,&\tau+\beta
m_{j}(||\bm{X}_{i}||)+\epsilon_{i}>0\\\ 0,&\tau+\beta
m_{j}(||\bm{X}_{i}||)+\epsilon_{i}\leq 0\end{Bmatrix},$ $\displaystyle
Y_{i}(0)$ $\displaystyle=\begin{Bmatrix}1,&\beta
m_{j}(||\bm{X}_{i}||)+\epsilon_{i}>0\\\ 0,&\beta
m_{j}(||\bm{X}_{i}||)+\epsilon_{i}\leq 0\end{Bmatrix},$

where $\epsilon_{i}\sim LG(0,1)$, $LG$ denotes a logistic distribution,
$\beta=0.5$, $\tau=0.5$, and
$P(D_{i}=1)=P(\gamma_{1}+\gamma_{2}||\bm{X}_{i}||>v_{i})=F_{LG}(\gamma_{1}+\gamma_{2}||\bm{X}_{i}||)$,
$\gamma_{1}=0.15$, $\gamma_{2}=0.4$, $v_{i}\sim LG(0,1)$ and $F_{LG}$ is the
distribution function. Other components of the specification are as in the
previous setting.

Observe that in this setting
$ATE_{i}=\mathbbm{E}(Y_{i}(1)-Y_{i}(0)|\bm{X}_{i}=\bm{x}_{i},Z_{i}=z_{i})=P(Y_{i}(1)=1|\bm{X}_{i}=\bm{x}_{i},Z_{i}=z_{i})-P(Y_{i}(0)=0|\bm{X}_{i}=\bm{x}_{i},Z_{i}=z_{i})=F_{LG}(\tau+\beta
m_{j}(||\bm{X}_{i}||))-F_{LG}(\beta m_{j}(||\bm{X}_{i}||))$, which has the
same analytical expression for the ATET, but conditional on the relevant
sample ($D_{i}=1$).

Our results are based on a sample size equal to 100 (Otsu and Rai,,
2017).111We perform other simulation exercises using larger sample sizes.
Outcomes show similar results. Available upon authors request. We have 1,000
replications for each simulation exercise.

First, we present the results for the average treatment effects on treated for
the continuous outcome variable.222Results for the average treatment effects
are in the Appendix, subsection 6.2. Outcomes are similar to ATET. Table 1
displays the mean relative squared errors (MRSE),

$MRSE=\frac{1}{S_{k}}\sum_{s=1}^{S_{k}}\left(\frac{\hat{\Delta}^{BC}_{s}(k)-ATET}{\hat{\Delta}^{BC}_{s}(k^{*})-ATET}\right)^{2},\
k^{*}\neq k,k=1,2,\dots 20,$

where $\hat{\Delta}^{BC}_{s}(k)$ and $\hat{\Delta}^{BC}_{s}(k^{*})$ are the
ATET estimates using any $k$ and the optimal $k=k^{*}$, and $S_{k}$ is the
number of times that $k$ was not optimal.

We see in Table 1 that there are relative large advantages using our optimal
$k^{*}$ as all figures are larger than 1.

Table 1: Mean relative squared errors continuous outcome: Average treatment effects on treated $k$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$
---|---|---|---|---|---|---
1 | 479,945.8 | 361.6 | 71.9 | 74.9 | 64,319.7 | 5,308.1
2 | 14,005.9 | 194.9 | 90.9 | 54.2 | 462.4 | 24,887.9
3 | 658.6 | 126.7 | 129.6 | 10.4 | 18,965.7 | 7,391.0
4 | 25,471.9 | 85.2 | 118.5 | 14.5 | 11.4 | 92.9
5 | 805.4 | 43.1 | 121.7 | 11.0 | 2,366.4 | 101.1
6 | 9,758.6 | 19.9 | 138.3 | 19.5 | 7,236.0 | 290.2
7 | 15,678.7 | 20.7 | 142.0 | 21.0 | 29,098.9 | 171.5
8 | 1,519.8 | 16.6 | 177.4 | 12.9 | 27,464.4 | 183.4
9 | 29,276.6 | 14.4 | 193.9 | 10.2 | 23,308.7 | 336.1
10 | 10,104.8 | 19.2 | 220.3 | 14.0 | 38,362.6 | 574.2
11 | 1,616.9 | 19.8 | 249.1 | 15.1 | 41,302.2 | 709.1
12 | 6,571.2 | 22.3 | 251.4 | 28.6 | 38,928.3 | 1,318.5
13 | 23,269.0 | 25.5 | 336.5 | 47.6 | 40,001.6 | 1,474.0
14 | 14,743.3 | 27.5 | 400.7 | 38.0 | 46,703.4 | 2,351.0
15 | 10,691.7 | 29.2 | 395.8 | 55.7 | 42,344.9 | 3,035.4
16 | 23,648.5 | 33.4 | 424.8 | 92.7 | 44,138.8 | 3,372.0
17 | 19,488.5 | 41.6 | 459.2 | 126.6 | 45,173.9 | 4,027.6
18 | 20,453.9 | 44.2 | 502.1 | 172.5 | 41,846.3 | 3,577.2
19 | 26,292.6 | 44.3 | 486.7 | 214.3 | 50,218.6 | 3,886.7
20 | 29,382.8 | 45.0 | 547.4 | 263.2 | 51,969.0 | 2,973.6
$MRSE=\frac{1}{S_{k}}\sum_{s=1}^{S_{k}}\left(\frac{\hat{\Delta}^{BC}_{s}(k)-ATET}{\hat{\Delta}^{BC}_{s}(k^{*})-ATET}\right)^{2},\
k^{*}\neq k,k=1,2,\dots 20$ where
$\hat{\Delta}^{BC}_{s}(k)$ and $\hat{\Delta}^{BC}_{s}(k^{*})$ are the ATET
estimates using any $k$ and
the optimal $k=k^{*}$, and $S_{k}$ is the number of times that $k$ was not
optimal. It seems that there are relative large advantages using
our optimal $k^{*}$.

Figure 1 shows the average over 1,000 replications of the ratios between the
length of the 95% confidence intervals associated with any number of neighbors
($k$, $x$ axis) and the optimal selection using our proposal. This is done for
the six data generating settings ($m_{j}(||\bm{X}_{i}||),j=1,2,\dots 6$). We
see in this figure that a small $k$ implies relative wide intervals, whereas
large $k$ implies narrow intervals. This is expected as small $k$ implies high
variability compared with large $k$. Observe that intervals lengths decreases
very fast, and stabilize around $k=10$. On the other hand, we see in Figure 2
that the type I error due to rejecting the null hypothesis $ATET=\tau$, given
a significance level equal to 5%, is the lowest using $k=1$. Using $k$
minimizing the mean squared error gives a relatively small type I error
compared with other $k$ selection.

Figure 3 displays the histograms of the optimal $k$ in our six simulation
settings. We see that optimal selection is dominated by small $k$. This could
be evidence to support $k=1$, which is widely used in many applications.
Although, it seems that $k=\left\\{2,3\right\\}$ are the most relevant, and
large values of $k$ cannot be discarded, especially when using large sample
sizes (we performed other simulations exercises that show this).

Figure 1: 95% confidence interval length ratio continuous outcome: Average
treatment effects on treated. Small $k$ implies relative wide intervals,
whereas large $k$ implies narrow intervals. The dotted line is optimal
selection of $k$ minimizing mean squared error. By construction has a ratio
equal to 1. Figure 2: Type I errors (5% nominal size) continuous outcome:
Average treatment effects on treated. Optimal $k$ gives a relatively small
type I error compared with other $k$ selection. $k=1$ always gets the best
nominal size. Figure 3:  Optimal $k$ frequencies continuous outcome: Average
treatment effects on treated. It seems that small $k$ dominates optimal $k$
choice. Although, large $k$ cannot be discarded.

Results of the ATET for our simulation setting using a binary outcome are
qualitatively similar to the continuous outcome. Table 2 show the mean
relative squared errors. Again, all figures are larger than 1, which means
that using a non-optimal $k$ implies larger estimation errors.

Table 2: Mean relative squared errors binary outcome: Average treatment effects on treated $k$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$
---|---|---|---|---|---|---
1 | 168.8 | 143.8 | 185.1 | 222.2 | 8,684.6 | 2,581.9
2 | 103.1 | 95.7 | 72.8 | 52.1 | 51.5 | 1,342.0
3 | 54.6 | 54.3 | 69.7 | 56.1 | 56.1 | 737.1
4 | 56.4 | 21.3 | 23.9 | 17.9 | 114.1 | 441.2
5 | 50.0 | 23.1 | 30.6 | 27.0 | 481.3 | 478.2
6 | 37.2 | 13.7 | 15.6 | 20.9 | 1,460.7 | 466.2
7 | 51.7 | 13.1 | 18.3 | 24.3 | 1,228.6 | 256.8
8 | 18.3 | 18.5 | 26.5 | 22.3 | 1,917.7 | 230.9
9 | 34.8 | 17.0 | 26.1 | 27.0 | 2,316.8 | 317.2
10 | 86.7 | 26.4 | 34.6 | 36.3 | 2,720.1 | 463.3
11 | 178.4 | 34.6 | 44.4 | 62.4 | 2,198.2 | 352.9
12 | 188.5 | 36.4 | 39.2 | 54.6 | 2,082.1 | 314.2
13 | 246.8 | 36.0 | 38.8 | 60.4 | 3,024.4 | 363.5
14 | 327.1 | 41.5 | 42.7 | 68.9 | 3,977.1 | 517.1
15 | 383.3 | 47.5 | 47.5 | 88.9 | 5,120.8 | 709.2
16 | 316.1 | 50.2 | 51.8 | 92.2 | 7,136.7 | 888.7
17 | 292.2 | 51.1 | 46.4 | 98.8 | 8,163.2 | 1,222.5
18 | 274.7 | 54.6 | 54.4 | 110.8 | 8,501.3 | 1,175.7
19 | 254.6 | 55.6 | 60.3 | 110.4 | 8,730.9 | 1,105.3
20 | 248.2 | 61.1 | 64.3 | 124.8 | 9,135.4 | 1,050.8
$MRSE=\frac{1}{S_{k}}\sum_{s=1}^{S_{k}}\left(\frac{\hat{\Delta}^{BC}_{s}(k)-ATET}{\hat{\Delta}^{BC}_{s}(k^{*})-ATET}\right)^{2},\
k^{*}\neq k,k=1,2,\dots 20$
where $\hat{\Delta}^{BC}_{s}(k)$ and $\hat{\Delta}^{BC}_{s}(k^{*})$ are the
ATET estimates using
any $k$ and the optimal $k=k^{*}$, and $S_{k}$ is the number of times
that $k$ was not optimal. It seems that there are relative
large advantages using an optimal $k$.

Figure 4 shows the relative 95% confidence intervals length ratios associated
with $k$ and $k^{*}$. All six simulation settings look very similar with wider
confidence intervals when $k$ is small, and narrower intervals when $k$ is
large. There is stabilization around $k=10$. Figure 5 shows type I errors
using 5% as significance level. Larger $k$ implies larger errors, and using
$k^{*}$ implies small type I errors in general, except when $k$ is very small.
Small type I error can be explained by means of Figure 6 where we see that
optimal selection of $k$ is also based on small $k$. Here
$k^{*}=\left\\{2,3,4,5\right\\}$ are particularly relevant.

In general, we can observe that our proposal implies relatively smaller
estimation errors achieving a compromise between confidence interval length
and type I error having as a result of small optimal $k$ most of the times.

Figure 4: 95% confidence interval length ratio binary outcome: Average
treatment effects on treated. Small $k$ implies relative wide intervals,
whereas large $k$ implies narrow intervals. The dotted line is optimal
selection of $k$ minimizing mean squared error. By construction has a ratio
equal to 1. Figure 5: Type I errors (5% nominal size) binary outcome: Average
treatment effects on treated. Optimal $k$ gives a relatively small type I
error compared with other $k$ selection. $k=1$ always gets the best nominal
size. Figure 6:  Optimal $k$ frequencies binary outcome: Average treatment
effects on treated. It seems that small $k$ dominates optimal $k$ choice.
Although, large $k$ cannot be discarded.

##  4 Empirical example

To illustrate our proposal we estimate the average treatment effect on treated
of 401(k) participation on accumulated net financial assets. 401(k) is a
retirement plan where savings contributions are provided by employers deduced
from employees’ payment before taxes such that taxes on capital gains are
avoided. Employees should achieve some eligibility criteria such that the plan
does not have universal applicability in principle. The impact evaluation of
401(k) on net financial assets has been done previously by Benjamin, (2003);
Chernozhukov and Hansen, (2004); Conley et al., (2012); Chernozhukov et al.,
(2017) among others. We overcome the lack of random assignment following
Chernozhukov et al., (2017) arguments. Enrollment in a 401(k) plan is
exogenous when controlling for observable variables which drive employment
decisions, such as income, when the plan initially became available.
Therefore, our dataset is from the 1991 Survey of Income and Program
Participation composed by 9,915 observations.

We control for the same set of variables as Chernozhukov and Hansen, (2004).
In particular, age (linear and squared), education (years), household size,
income (seven levels), and binary variables indicating: defined benefit
pension, home ownership, marital status (married), participation in IRA
(individual retirement account) and two-earners status (both household heads
contribute to household income). Details can be found in Benjamin, (2003).

Table 3 shows unconditional mean difference tests of our outcome variables:
accumulated net financial assets, and a binary variable indicating positive
accumulated net financial assets. Employees enrolled in 401(k) plans have
$27.4 more accumulated net financial assets on average than employees who are
not enrolled. The unconditional probability of positive net financial assets
is 31% on average higher for the former. We can see that there are
unconditional statistically significant differences.

Table 3: Descriptive statistics: Outcome variables Response | Mean treated | Mean untreated | Statistic
---|---|---|---
Continuous | $38.26K (79.09) | $10.89K (55.26) | 16.28∗∗
Binary | 0.86 (0.35) | 0.55 (0.5) | 35.27 8∗∗
Standard error in parenthesis. The test of the difference of two proportions
uses a normally
distributed test statistic calculated as
$z=\frac{\hat{p}_{1}-\hat{p}_{2}}{(\hat{p}_{p}(1-\hat{p}_{p})(1/n_{1}+1/n2))^{0.5}}$,
$\hat{p}_{p}=\frac{x_{1}+x_{2}}{n_{1}+n_{2}}$, $x_{l}$ and $n_{l}$ are
the number of successes and sample sizes in each group, $l=1,2$. The test of
the difference
of two means uses a Student’s t calculated as
$t=\frac{\bar{x}_{1}-\bar{x}_{2}}{(s_{x_{1}}^{2}/n_{1}+s_{x_{2}}^{2}/n_{2})^{0.5}}$.
1 Null hypothesis: means are equal. Critical value at 5% level of significance
is 1.96
∗∗ Rejection of null hypothesis at 5% significance level.

Figure 7 shows the ATET in the accumulated net financial assets using our
optimal selection, $k^{*}=19$, and $k=1$, as a reference. We see similar ATET
point estimates, which are approximately $ 15K (similar outcomes are found by
Conley et al., (2012) under the exogeneity assumption using a Bayesian
approach for a local average treatment effects, LATE). However, the 95%
confidence interval using $k^{*}$ is ($10.9K,$18.9K), whereas $k=1$ interval
is ($9.1K,$19.9K), that is 35% larger. Results regarding the probability of
positive accumulated net financial assets can be seen in Figure 8. Both,
$k^{*}=19$ and $k=1$ indicate statistically significant ATET whose point
estimate is approximately 19.9% and 18.7% respectively. The 95% confidence
intervals are (16.3%, 21.3%) and (17.7%, 21.3%) for $k=1$ and $k^{*}=19$. This
means that using $k=1$ implies a 38.9% wider interval than using $k^{*}=19$.
Observe that both exercises indicate the same optimal $k$. This means that
diagnostic balance tests are based on the same control groups. Figure 9 shows
the love plot indicating that control groups based on $k^{*}=19$ satisfy mean
balancing conditions.

Figure 7: Average treatment effects on treated due to 401(k) participation on
accumulated net financial assets: Point estimates with $k=1$ (gray dot) and
optimal $k^{*}$ (black dot), and 95% confidence intervals. Both, $k=1$ and
$k^{*}$ indicate 5% significant positive ATET of approximately $15K. 95%
confidence intervals using $k=1$ are 28% wider compared to $k^{*}$. Figure 8:
Average treatment effects on treated due to 401(k) participation on the
probability of positive accumulated net financial assets: Point estimates with
$k=1$ (gray dot) and optimal $k^{*}$ (black dot), and 95% confidence
intervals. Both, $k=1$ and $k^{*}$ indicate 5% significant positive ATET of
approximately 19%. 95% confidence intervals using $k=1$ are 28% wider compared
to $k^{*}$. Figure 9: Mean balancing condition for control variables used
calculating ATET: Love plot using $k^{*}=19$. It seems that mean balancing
conditions are satisfied as statistical tests (black dots) are inside critical
levels (dashed red lines).

##  5 Conclusions and final remarks

We present a simple approach to obtain an optimal selection of the number of
control units ($k$) to build the “synthetic” or “counterfactual” unit for each
individual in the k nearest neighbors algorithm. Our approach is based on a
simple unbiased reference estimator (individual treatment effect on treated
for ATET and individual treatment effect for ATE) such that we select $k$
minimizing the mean squared error with respect to the reference estimator.
Therefore, the optimal $k=k^{*}$ achieves a balance between confidence
interval length and type I error.

Our application suggests that $k^{*}$ can be relatively large compared to the
common practice of using $k=1$. This implies significant reductions regarding
confidence interval lengths. Although, both choices show 5% statistically
significant effects of 401(k) enrollment in the probability of positive
accumulated net financial assets, and its amount.

A final remark is that using as reference estimator for treatments effects the
one obtained using $k=1$ is an interesting alternative. We obtained same
qualitative results, that is, a balance between interval length and type I
error, and a considerable reduction in MRSE. Simulations exercises confirm
this.

## References

  * Abadie and Imbens, (2006) Abadie, A. and Imbens, G. (2006).  Large sample properties of matching estimators for average treatment effects.  Econometrica, 74(1):235–267.
  * Abadie and Imbens, (2008) Abadie, A. and Imbens, G. (2008).  On the failure of the bootstrap for matching estimators.  Econometrica, 76(6):1537––1557.
  * Abadie and Imbens, (2011) Abadie, A. and Imbens, G. (2011).  Bias-corrected matching estimators for average treatment effects.  Journal of Business & Economic Statistics, 29(1):1–11.
  * Abadie and Imbens, (2016) Abadie, A. and Imbens, G. (2016).  Matching on the estimated propensity score.  Econometrica, 84(2):781––807.
  * Angrist and Krueger, (1999) Angrist, J. and Krueger, A. (1999).  Empirical strategies in labor economics.  In Ashenfelter, O. and Card, D., editors, Handbook of Labor Economics, volume 3. ELSEVIER.
  * Athey et al., (2015) Athey, S., Imbens, G., and Ramachandra, V. (2015).  Machine learning methods for estimating heterogeneous causal effects.  Technical report, Stanford University.
  * Benjamin, (2003) Benjamin, D. (2003).  Does 401(k) eligibility increase savings? Evidence from propensity score subclassification.  Journal of Public Economics, 87:1259–1290.
  * Cameron and Trivedi, (2005) Cameron, C. and Trivedi, K. (2005).  Microeconometrics: Methods and Applications.  Cambridge University Press.
  * Chernozhukov et al., (2017) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2017).  Double/debiased machine learning for treatment and causal parameters.  The Econometrics Journal, 21:1–68.
  * Chernozhukov and Hansen, (2004) Chernozhukov, V. and Hansen, C. (2004).  The effects of 401 (k) participation on the wealth distribution: an instrumental quantile regression analysis.  Review of Economics and statistics, 86(3):735–751.
  * Conley et al., (2012) Conley, T., Hansen, C., and Rossi, P. (2012).  Plausibly exogenous.  The Review of Economics and Statistics, 94(1):260–272.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning Data Mining, Inference, and Prediction.  Springer, second edition edition.
  * Otsu and Rai, (2017) Otsu, T. and Rai, Y. (2017).  Bootstrap inference of matching estimators for average treatment effects.  Journal of the American Statistical Association, 112(520):1720–1732.
  * Rosenbaum and Rubin, (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55.
  * Rubin, (1978) Rubin, D. (1978).  Bayesian inference for causal effects: The role of randomization.  Annals of Statistics, 6:34–58.

##  6 Appendix

###  6.1 Optimal selection of $k$ in kNN: Average treatment effects

The average treatment effect is

$\Delta^{ATE}=\mathbb{E}\left[(Y_{i}(1)-Y_{i}(0))|\bm{X}_{i}=\bm{x}_{i}\right].$

A matching estimator for the ATE is given by

$\displaystyle\hat{\Delta}^{ATE}=\frac{1}{N}\sum_{i=1}^{N}\left[\hat{Y}_{i}(1)-\hat{Y}_{i}(0)\right],$

where

$\displaystyle\hat{Y}_{i}(1)$
$\displaystyle=\begin{Bmatrix}\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x})}Y_{j}(0),&D_{i}=0\\\ Y_{i},&D_{i}=1\end{Bmatrix},$

$\displaystyle\hat{Y}_{i}(0)$ $\displaystyle=\begin{Bmatrix}Y_{i},&D_{i}=0\\\
\frac{1}{k}\sum_{j\in A_{k}^{i}(\bm{x})}Y_{j}(1),&D_{i}=1\end{Bmatrix},$

$A_{k}^{i}(\bm{x})$ defines the relevant control group conditional on being
treated or untreated.

The bias of the ATE is

$B^{ATE}=\frac{1}{N}\left\\{\sum_{i\in\left\\{D_{i}=1\right\\}}\left[\mu_{0}(\bm{x}_{i})-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x})}\mu_{0}(\bm{x}_{j})\right]-\sum_{i\in\left\\{D_{i}=0\right\\}}\left[\mu_{1}(\bm{x}_{i})-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x})}\mu_{1}(\bm{x}_{j})\right]\right\\},$

where $\mu_{0}(\bm{x}_{i})=\mathbb{E}\left[Y|\bm{X}=\bm{x},D=0\right]$ and
$\mu_{1}(\bm{x}_{i})=\mathbb{E}\left[Y|\bm{X}=\bm{x},D=1\right]$. These can be
estimated non-parametrically using a series expansion estimator to obtain
$\hat{B}^{TM}$ (Abadie and Imbens,, 2011).

Therefore, the bias-corrected matching estimator is

$\hat{\Delta}^{BCATE}=\hat{\Delta}^{ATE}-\hat{B}^{ATE}.$

Theorem 4(i) in Abadie and Imbens, (2006) shows asymptotic distribution
convergence results for the bias corrected average treatment effect, and
section 3.2 in their paper has required variance expressions (marginal and
conditional).

The individual treatment effect is

${\Delta_{i}^{ATE*}}=\frac{Y_{i}(D_{i}-P(\bm{X}_{i}))}{P(\bm{X}_{i})(1-P(\bm{X}_{i}))},$
(8)

such that
$\mathbb{E}\left[{\Delta_{i}^{ATE*}}\big{|}\bm{X}_{i}=\bm{x}\right]=\Delta^{ATE}$
(Athey et al.,, 2015).

Therefore, the optimal selection of $k$ is given by solving the program

$\operatorname*{arg\,min}_{{k}}\frac{1}{G}\sum_{g=1}^{G}\left\\{\left(\hat{\Delta}_{g}^{BCATE}-\frac{1}{N_{g}}\sum_{i=1}^{N_{g}}\widehat{\Delta_{i,g}^{ATE*}}\right)^{2}\right\\},$

where

$\widehat{\Delta_{i,g}^{ATE*}}=\left(\frac{Y_{i,g}^{Test}(D_{i,g}^{Test}-\widehat{P(\bm{X}_{i,g}^{Test})})}{\widehat{P(\bm{X}_{i,g}^{Test})}(1-\widehat{P(\bm{X}_{i,g}^{Test})})})\right),$

$\displaystyle\hat{\Delta}_{g}^{ATE}=$
$\displaystyle\frac{1}{N_{g}}\left\\{\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\left[Y_{i,g}^{Test}(1)-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(0)\right]\right.$
$\displaystyle\left.-\sum_{i\in\left\\{D_{i,g}^{Test}=0\right\\}}\left[Y_{i,g}^{Test}(0)-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}Y_{j,g}^{Train}(1)\right]\right\\},$

$\displaystyle\hat{B}_{g}^{ATE}=$
$\displaystyle\frac{1}{N_{g}}\left\\{\sum_{i\in\left\\{D_{i,g}^{Test}=1\right\\}}\left[\hat{\mu}_{0}(\bm{x}_{i})_{g}^{Test}-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{0}(\bm{x}_{j})_{g}^{Train}\right]\right.$
$\displaystyle-\left.\sum_{i\in\left\\{D_{i,g}^{Test}=0\right\\}}\left[\hat{\mu}_{1}(\bm{x}_{i})_{g}^{Test}-\frac{1}{k}\sum_{j\in
A_{k}^{i}(\bm{x}_{g}^{Test})}\hat{\mu}_{1}(\bm{x}_{j})_{g}^{Train}\right]\right\\},$

and

$\hat{\Delta}_{g}^{BCATE}=\hat{\Delta}_{g}^{ATE}-\hat{B}_{g}^{ATE}.$

An algorithm similar to Algorithm A1 can be implemented to obtain an optimal
$k$ for the ATE given these definitions.

###  6.2 Simulation outcomes: Average treatment effects

Table 4: Mean relative squared errors continuous outcome: Average treatment effects $k$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$
---|---|---|---|---|---|---
1 | 709.3 | 4,511.1 | 42.0 | 790.9 | 586.1 | 249.2
2 | 12.8 | 31.1 | 97.7 | 403.6 | 200.8 | 63.2
3 | 9.3 | 162.6 | 104.7 | 329.2 | 272.2 | 85.2
4 | 20.4 | 805.9 | 26.0 | 17.5 | 54.9 | 49.1
5 | 7.5 | 569.2 | 19.5 | 34.7 | 87.3 | 41.1
6 | 11.4 | 382.9 | 25.5 | 39.5 | 539.8 | 27.3
7 | 14.2 | 522.7 | 29.1 | 100.4 | 949.7 | 41.5
8 | 11.8 | 855.1 | 31.7 | 171.5 | 1,055.5 | 46.0
9 | 7.6 | 495.1 | 34.4 | 563.4 | 1,104.3 | 37.2
10 | 9.8 | 697.9 | 40.9 | 856.8 | 946.1 | 55.8
11 | 8.4 | 483.9 | 47.7 | 888.2 | 949.8 | 82.0
12 | 12.8 | 372.7 | 52.9 | 1,050.7 | 1,073.7 | 71.3
13 | 12.4 | 298.0 | 49.5 | 1,557.7 | 1,240.0 | 65.3
14 | 10.7 | 360.9 | 56.7 | 1,759.5 | 1,445.1 | 49.8
15 | 9.1 | 371.0 | 67.7 | 1,818.8 | 1,742.4 | 37.8
16 | 14.7 | 429.1 | 68.6 | 1,756.0 | 1,915.3 | 38.3
17 | 17.8 | 445.8 | 90.8 | 1,746.5 | 2,096.7 | 57.8
18 | 10.9 | 382.9 | 104.9 | 1,782.8 | 2,252.2 | 84.7
19 | 9.2 | 386.4 | 146.7 | 1,623.2 | 2,586.7 | 155.0
20 | 9.4 | 233.2 | 184.2 | 1,739.2 | 2,722.4 | 212.9
$MRSE=\frac{1}{S_{k}}\sum_{s=1}^{S_{k}}\left(\frac{\hat{\Delta}^{BC}_{s}(k)-ATE}{\hat{\Delta}^{BC}_{s}(k^{*})-ATE}\right)^{2},\
k^{*}\neq k,k=1,2,\dots 20$
where $\hat{\Delta}^{BC}_{s}(k)$ and $\hat{\Delta}^{BC}_{s}(k^{*})$ are the
ATE estimates using any $k$
and the optimal $k=k^{*}$, and $S_{k}$ is the number of times that $k$ was
not optimal. It seems that there are relative large advantages
using our optimal $k^{*}$.
Figure 10: 95% confidence interval length ratio continuous outcome: Average treatment effects. Small $k$ implies relative wide intervals, whereas large $k$ implies narrow intervals. The dotted line is the optimal selection of $k$ minimizing mean squared error. By construction has a ratio equal to 1. Figure 11: Type I errors continuous outcome: Average treatment effects. The optimal $k$ gives an average type I error compared with other $k$ selection. $k=1$ always gets the best nominal size Table 5: Mean relative squared errors binary outcome: Average treatment effects $k$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$
---|---|---|---|---|---|---
1 | 73.5 | 123.3 | 191.2 | 275.5 | 194.6 | 1,444.2
2 | 14.5 | 43.5 | 838.5 | 209.1 | 25.3 | 187,820.3
3 | 12.7 | 41.6 | 906.2 | 266.3 | 42.3 | 378,910.0
4 | 35.3 | 37.7 | 27.5 | 433.5 | 106.5 | 508,679.5
5 | 15.7 | 28.7 | 26.1 | 308.9 | 149.3 | 519,137.8
6 | 12.5 | 30.2 | 22.5 | 511.0 | 150.7 | 626,520.8
7 | 12.8 | 49.5 | 38.6 | 754.6 | 148.2 | 739,373.9
8 | 9.4 | 81.3 | 40.0 | 1,114.0 | 182.0 | 826,213.6
9 | 19.0 | 59.8 | 31.3 | 1,048.2 | 133.6 | 755,273.2
10 | 12.7 | 79.5 | 23.2 | 1,199.5 | 134.1 | 764,009.2
11 | 13.9 | 44.8 | 26.8 | 988.9 | 140.7 | 786,963.7
12 | 16.5 | 44.5 | 30.9 | 1,140.5 | 131.7 | 715,959.8
13 | 24.7 | 33.5 | 18.0 | 1,375.7 | 139.4 | 727,003.1
14 | 18.5 | 25.1 | 25.9 | 1,362.9 | 157.4 | 697,692.1
15 | 23.3 | 20.4 | 48.2 | 1,213.5 | 194.6 | 711,163.3
16 | 26.3 | 15.1 | 13.0 | 1,198.5 | 205.3 | 738,156.2
17 | 23.8 | 15.0 | 14.2 | 1,269.1 | 196.7 | 695,079.4
18 | 17.9 | 14.6 | 11.9 | 1,325.0 | 220.8 | 661,129.6
19 | 18.2 | 14.7 | 12.0 | 1,267.1 | 235.7 | 615,817.0
20 | 19.4 | 16.2 | 24.4 | 1,381.8 | 254.8 | 647,296.8
$MRSE=\frac{1}{S_{k}}\sum_{s=1}^{S_{k}}\left(\frac{\hat{\Delta}^{BC}_{s}(k)-ATE}{\hat{\Delta}^{BC}_{s}(k^{*})-ATE}\right)^{2},\
k^{*}\neq k,k=1,2,\dots 20$
where $\hat{\Delta}^{BC}_{s}(k)$ and $\hat{\Delta}^{BC}_{s}(k^{*})$ are the
ATE estimates using
any $k$ and the optimal $k=k^{*}$, and $S_{k}$ is the number of times
that $k$ was not optimal. It seems that there are relative
large advantages using an optimal $k$.
Figure 12: 95% confidence interval length ratio binary outcome: Average
treatment effects. Small $k$ implies relative wide intervals, whereas large
$k$ implies narrow intervals. The dotted line is optimal selection of $k$
minimizing mean squared error. By construction has a ratio equal to 1. Figure
13: Type I errors (5% nominal size) binary outcome: Average treatment effects.
The optimal $k$ gives a relatively small type I error compared with other $k$
selection. $k=1$ always gets the best nominal size.