reason is that some experimental designs will make the GR estimator invariant
to certain parameter choices. To demonstrate this, consider the setting of a
two-arm completely randomized design with one pretreatment covariate. Let 1
denote the treatment arm and 0 denote the control arm, and let $n_{1}$ and
$n_{0}$ denote the number of treated and control units, respectively. Define
$\bar{x}=\frac{1}{n}\sum_{i}x_{i}=0$ (centered),
$\hat{\bar{x}}_{1}=\frac{1}{n_{1}}\sum_{i}\mathds{R}_{1i}x_{i}$, and
$\hat{\bar{x}}_{0}=\frac{1}{n_{0}}\sum_{i}\mathds{R}_{0i}x_{i}$, and note
$\hat{\bar{x}}_{0}=-\frac{n_{1}}{n_{0}}\hat{\bar{x}}_{1}$. Suppose we use a
separate-slope linear model to adjust
$f^{0}(x_{i},\beta)=\beta_{0}+\beta_{1}x_{i}$ and
$f^{1}(x_{i},\beta)=\beta_{2}+\beta_{3}x_{i}$. The GR estimator for the ATE is

$\displaystyle\begin{split}&(\beta_{2}-\beta_{0})+\frac{1}{n}\sum_{i}x_{i}(\beta_{3}-\beta_{1})-\left(\frac{1}{n_{1}}\sum_{i}\mathds{R}_{1i}(\beta_{2}+x_{i}\beta_{3})-\frac{1}{n_{0}}\sum_{i}\mathds{R}_{0i}(\beta_{0}+x_{i}\beta_{1})\right)\\\
&=0-\hat{\bar{x}}_{1}\beta_{3}+\hat{\bar{x}}_{0}\beta_{1}=-\hat{\bar{x}}_{1}(\beta_{3}+\frac{n_{1}}{n_{0}}\beta_{1}).\end{split}$
(58)

We observe that the intercepts $\beta_{0}$ and $\beta_{2}$ are canceled and
that there are multiple pairs of $(\tilde{\beta}_{1},\tilde{\beta}_{3})$ that
are equivalent to $(\beta_{1},\beta_{3})$, such as $\tilde{\beta}_{1}=0$ and
$\tilde{\beta}_{3}=\beta_{3}+\frac{n_{1}}{n_{0}}\beta_{1}$. For nonlinear
models, the problem is more subtle. For example, for a separate-slope logistic
model, we have the GR estimator

$\displaystyle\begin{split}&\frac{1}{n}\sum_{i}\left(\frac{\exp(\beta_{2}+\beta_{3}x_{i})}{1+\exp(\beta_{2}+\beta_{3}x_{i})}-\frac{\exp(\beta_{0}+\beta_{1}x_{i})}{1+\exp(\beta_{0}+\beta_{1}x_{i})}\right)\\\
&-\frac{1}{n_{1}}\sum_{i}\mathds{R}_{1i}\frac{\exp(\beta_{2}+\beta_{3}x_{i})}{1+\exp(\beta_{2}+\beta_{3}x_{i})}+\frac{1}{n_{0}}\sum_{i}\mathds{R}_{0i}\frac{\exp(\beta_{0}+\beta_{1}x_{i})}{1+\exp(\beta_{0}+\beta_{1}x_{i})}.\end{split}$
(59)

When $\beta_{1}$ and $\beta_{3}$ are distant from 0, the intercepts will have
an impact on the estimator. However, if $\beta_{1}$ and $\beta_{3}$ are close
to 0, the intercepts will cancel, as in the linear case. We note that the
identification problem highlighted here is design specific. For example, this
problem will also arise for pairwise randomized designs but not for Bernoulli
designs.353535Note in a pairwise randomized design, a GR estimator using
linear models is also invariant to variables that have no within pair
variation. In general, there are two sources of the weak/non-identification
problem. The first is weak/non identification of arm-specific intercepts. This
typically happens when a design has the same assignment probabilities for all
units and the numbers of units in all arms are fixed in all possible random
treatment allocations. This is can be checked by inspecting the eigenvalues of
the matrix $\frac{1}{n}\mathbf{1}^{\prime}\mathbf{D}\mathbf{1}$. Small
eigenvalues correspond to possible weak/non identification problems for the
intercepts. Another source of weak/non identification is that the design
induces co-linearity or cancellation of covariates in different treatment
arms. This can be checked by inspecting the eigenvalues of
$\frac{1}{n}{}\mathbf{x}\mathbf{D}\mathbf{x}$. Small eigenvalues correspond to
possible cancellation of covariates among different arms.

To use the Opt-GR estimators, we recommend that researchers inspect the
eigenvalues of $\frac{1}{n}\mathbf{x}\mathbf{D}\mathbf{x}$ before applying the
estimators. If one of the eigenvalues is small, the researcher may wish to
avoid specifying an intercept for the treatment arms or use a model with same-
slope adjustments instead of a model with separate-slope adjustments. As a
default, we recommend that researchers use a same-slope adjustment and choose
intercepts with care by inspecting
$\frac{1}{n}\mathbf{1}^{\prime}\mathbf{D}\mathbf{1}$.

We now compare the three classes of estimators in terms of implementation
difficulty and theoretical performance. The QMLE-GR estimators are the easiest
to implement in practice. The optimization problem is usually convex, and
efficient algorithms are widely available for solving such problems. No-harm
GR estimators are also simple to use. They are based on the QMLE estimates and
estimate the additional multiplicative constant in a straightforward manner.
Opt-GR estimators are the theoretically optimal estimators within the class of
models used for adjustment, but for nonlinear models, their implementation
usually requires solving a nonconvex problem.363636For most typical models
with arm-specific intercepts, such as linear, logit, and Poisson models, the
No-harm GR estimators are dominated by the Opt-GR estimators in terms of
asymptotic variances. This is because the class of models under consideration
does not change after one multiplies the imputation functions by a constant.

Besides the implementation difficulties, our asymptotic analysis does not
inform the finite-sample performance of the estimators. In general, when
$|||\mathbf{D}|||_{2}$ is large relative to the sample size, we expect the No-
harm GR and Opt-GR estimators to exhibit large finite-sample variability. We
examine their performances in the simulation section with various
$\frac{|||\mathbf{D}|||_{2}}{n}$ values.

We end this section by noting that one can hybrid the insights of No-harm GR
and Opt-GR estimators. We introduce here a class of Opt-I GR
estimators.373737The letter I stands for ”imputed”. This class of estimators
is Opt-GR estimators with linear models. Instead of using the full set of
covariates, the linear models have arm-specific intercepts and a single
covariate that is imputed by a QMLE model, for example, an OLS model (Opt-I
OLS GR) or a logit model (Opt-I Logit GR).383838In terms of asymptotic
variance, the Opt-GR estimator with a linear model is more efficient compared
with the Opt-I OLS GR estimator because they both use linear models for
adjustment. A theoretical comparison between the Opt-GR estimator with the
logit model and Opt-I Logit GR estimator is difficult because two estimators
use different classes of models for adjustment. We find this class of
estimators has good finite-sample performance in our simulations. We detail
the constructions below. Recall the definition of the intercept matrix
$\mathbf{1}$ in equation (5) and the stack of imputations $f(\theta)$ defined
before Section 6.

Algorithm 4 Opt-I estimator for
$\gamma_{c,n}=c^{\prime}\frac{1}{n}\mathbf{1}^{\prime}y$

A vector of QMLE coefficient estimates $\hat{\theta}_{n}$ from Algorithm 1, a
column contrast vector $c\in\mathbb{R}^{k}$, imputation functions
$f^{a}(\cdot,\theta)$, $a=1...,k$

  1. 1.

Define
$\hat{X}^{c}=\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})\begin{bmatrix}\mathbf{1}&f(\hat{\theta}_{n})\end{bmatrix}\in\mathbf{R}^{kn\times(p+1)}$.

  2. 2.

Define
$\hat{\beta}=\begin{bmatrix}\hat{\beta}^{1},...,\hat{\beta}^{k+1}\end{bmatrix}=(\hat{X}^{c^{\prime}}\mathbf{D}\hat{X}^{c})^{+}\hat{X}^{c^{\prime}}\mathbf{D}\boldsymbol{\pi}^{-1}\mathds{R}\mathbf{diag}(y)\mathbf{1}c\in\mathbf{R}^{k+1}$

  3. 3.

Compute $\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}$, for
$a=1,...,k$, where

$\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}=\frac{1}{n}\sum_{i=1}^{n}\left(\hat{\beta}^{a}+\hat{\beta}^{k+1}f^{a}(x_{i},\hat{\theta}_{n})\right)+\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\left(y_{ai}-\hat{\beta}^{a}-\hat{\beta}^{k+1}f^{a}(x_{i},\hat{\theta}_{n})\right).$
(60)

  4. 4.

Output
$\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}=\sum_{a=1}^{k}c_{a}\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}$.

The assumptions, asymptotic theory, and variance-bound estimation procedure
for the Opt-I estimators are similar to those of the No-harm GR estimators. We
leave them to Appendix B due to space constraints.

##  7 Examples

In this section, we provide two applications for our results. In Section 7.1,
we check Assumptions 6 and 9 for a two-arm completely randomized design.
Similar strategies can be employed to check Assumptions 6 and 9 for other
designs, such as pairwise randomized designs, cluster randomization designs
with uniform bounded cluster sizes, and stratified randomization designs with
a finite number of strata.393939A similar calculation appeared in Middleton,
2021b  with multiple errors. We reproduce the calculation having corrected
these errors. In Section 7.2, we demonstrate our framework by considering
QMLE-GR, No-harm GR, Opt-GR, and Opt-I GR estimators with linear and logit
models in network experimental designs (Aronow and Samii, (2017)).

###  7.1 Checking Assumption 6, $|||\mathbf{D}|||_{2}=O(1)$, and Assumption 9,
$\sigma_{\max}((\tilde{\mathbf{D}}\otimes\tilde{\mathbf{D}})\circ\mathbf{S})=o(n)$,
for a two-arm completely randomized design

Table 7.1 and Table 7.1 detail information for checking Assumptions 6 and 9 as
well as Assumptions 6′ and 9′ in Appendix J. Recall the setup of the design:
for a sample of size $n$, $n_{t}$ units are randomly selected into the
treatment group, and the rest of the $n_{c}=n-n_{t}$ units are selected into
the control group. We assume $\frac{n_{t}}{n}\to p$ as $n\to\infty$, with
$0<c<p<C<1$. The absolute row sum of the matrix $|||\mathbf{D}|||_{2}$ is
either $\frac{n_{t}}{n_{c}}+\frac{n_{c}}{n_{t}}+1+1=O(1)$ or
$\frac{n_{c}}{n_{t}}+\frac{n_{t}}{n_{c}}+1+1=O(1)$. Thus,
$|||\mathbf{D}|||_{2}=O(1)$. The last column checks Assumption 6′.
To check Assumption 9, we first note that the variance bound being used is the
Neyman bound, for example, see Imbens and Rubin, (2015).404040We check the
conditions for the Neyman bound as it is commonly used in practice. One can
check the conditions for the AS bound in the same way.  We can use the same
argument as above to establish
$|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}=O(1)$.
For the first quantity
$\sigma_{\max}((\tilde{\mathbf{D}}\otimes\tilde{\mathbf{D}})\circ\mathbf{S})$,
we use Lemma 24. Table 7.1 contains information about the slice indexed by
$\mathds{R}_{1i}$, $i=1,...,n$. For information about the slice indexed by
$\mathds{R}_{0i}$, $i=1,...,n$, one can simply exchange the roles of $n_{t}$
and $n_{c}$ in Table 7.1. By symmetry, we see
$\sigma_{\max}((\tilde{\mathbf{D}}\otimes\tilde{\mathbf{D}})\circ\mathbf{S})=O(1)$
using Lemma 24. One can also see by symmetry that
$\frac{1}{n}||\left(\tilde{\mathbf{D}}\otimes\tilde{\mathbf{D}}\right)\circ\mathbf{S}||_{1}=O(1)$,
as is required in A9′.
This calculation can be easily generalized to check Assumptions 6, 6′, 9, and
9′ for cluster randomization designs with uniform bounded cluster sizes and
stratified randomization designs with a finite number of strata.

Table 7.1: Entries for $\mathbf{D}$ for a two-arm completely randomized design

Partition | Entry Type | Count | Entry Value |  $\frac{1}{n}$Count$\scriptstyle\times$(Entry Value)
---|---|---|---|---
$\mathbf{D}_{11}$ | Diagonal | $\scriptstyle n$ | $\scriptstyle\frac{n_{c}}{n_{t}}$ | $\scriptstyle\frac{n_{c}}{n_{t}}=O(1)$
| Off-diagonal | $\scriptstyle n(n-1)$ | $\scriptstyle-\frac{n_{c}}{n_{t}(n-1)}$ | $\scriptstyle-\frac{n_{c}}{n_{t}}=O(1)$
$\mathbf{D}_{12}$ or $\mathbf{D}_{21}$ | Diagonal | $\scriptstyle 2n$ | $\scriptstyle-1$ | $\scriptstyle-2=O(1)$
| Off-diagonal | $\scriptstyle 2n(n-1)$ | $\scriptstyle\frac{1}{(n-1)}$ | $\scriptstyle 2=O(1)$
$\mathbf{D}_{22}$ | Diagonal | $\scriptstyle n$ | $\scriptstyle\frac{n_{t}}{n_{c}}$ | $\scriptstyle\frac{n_{t}}{n_{c}}=\scriptstyle O(1)$
| Off-diagonal | $\scriptstyle n(n-1)$ | $\scriptstyle-\frac{n_{t}}{n_{c}(n-1)}$ | $\scriptstyle-\frac{n_{t}}{n_{c}}=O(1)$

  * •

Note: $n$ is the number of units. $n_{t}$ is the number of units in the
treatment group, and $n_{c}$ is the number of units in the control group. The
matrix $\mathbf{D}$ is divided into four block matrices with
$\mathbf{D}_{ab}$, $a,b=1,2$. The block matrix $\mathbf{D}_{ab}$ consists of
positions between rows $(a-1)n+1$ and $an$ and columns $(b-1)n+1$ and $bn$.
The Entry Type column indicates whether the entries are diagonal or off-
diagonal in the $\mathbf{D}_{ab}$ matrix. The Count column counts the number
of entries of each type. The Entry Value column records the value of the entry
of each type.

Entries for
$(\tilde{\mathbf{D}}^{N}\otimes\tilde{\mathbf{D}}^{N})\circ\mathbf{S}$ for a
two-arm completely randomized design with a fixed $\mathds{R}_{1i}$ jkl
pattern
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{1j},\mathds{R}_{1k}\mathds{R}_{1l})\frac{\tilde{\mathbf{D}}^{N}_{ij}\tilde{\mathbf{D}}^{N}_{kl}}{\pi_{ij}\pi_{lk}}$
$(v_{111})$
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{1j},\mathds{R}_{0k}\mathds{R}_{0l})\frac{\tilde{\mathbf{D}}^{N}_{ij}\tilde{\mathbf{D}}^{N}_{kl}}{\pi_{ij}\pi_{lk}}$
$(v_{100})$ Count O(1)? i=j=k=l $\frac{n^{2}n_{c}}{n_{t}^{3}}$
$-\frac{n^{2}}{n_{t}n_{c}}$ 1 Y j$\not=$i, k=l=i
$-\frac{n^{2}n_{c}}{n_{t}^{3}(n-1)}$ $\frac{n^{2}}{(n-1)n_{c}n_{t}}$ n-1 Y
k$\not=$i, j=l=i $-\frac{n^{2}n_{c}}{n_{t}^{3}(n-1)}$
$\frac{n^{2}}{(n-1)n_{t}n_{c}}$ n-1 Y l$\not=$i, j=k=i
$-\frac{n^{2}n_{c}}{n_{t}^{3}(n-1)}$ $\frac{n^{2}}{(n-1)n_{t}n_{c}}$ n-1 Y
j$=$i, k$=$l $-\frac{n_{c}n^{2}}{(n-1)n_{t}^{3}}$
$\frac{n^{2}}{n_{t}n_{c}(n-1)}$ n-1 Y j$=$k, i$=$l
$\frac{n_{c}(n+n_{t}-1)n^{2}}{n_{t}(n_{t}-1)n_{c}^{2}(n-1)^{2}}$
-$\frac{n^{2}}{n_{t}n_{c}(n-1)^{2}}$ n-1 Y j$=$l, i$=$k
$\frac{n_{c}(n+n_{t}-1)n^{2}}{n_{t}(n_{t}-1)n_{c}^{2}(n-1)^{2}}$
-$\frac{n^{2}}{n_{t}n_{c}(n-1)^{2}}$ n-1 Y j$=$i, k$\not=$l$\not=$i
$\frac{2n_{c}n^{2}}{n_{t}^{3}(n-1)(n-2)}$
$\frac{2n^{2}}{n_{t}n_{c}(n-1)(n-2)}$ (n-1)(n-2) Y k=i, j$\not=$l$\not=$i
$\left[\frac{(n_{t}-2)n(n-1)}{(n-2)n_{t}(n_{t}-1)}-1\right]\frac{n^{2}}{n_{t}^{2}(n-1)^{2}}$
$-\frac{n^{2}}{n_{c}n_{t}(n-1)^{2}}$ (n-1)(n-2) Y l=i, j$\not=$k$\not=$i
$\left[\frac{(n_{t}-2)n(n-1)}{(n-2)n_{t}(n_{t}-1)}-1\right]\frac{n^{2}}{n_{t}^{2}(n-1)^{2}}$
$-\frac{n^{2}}{n_{c}n_{t}(n-1)^{2}}$ (n-1)(n-2) Y i$\not=$j$\not=$k$\not=$l
$\frac{(4n-6)n_{t}(n_{t}-1)-(4n_{t}-6)n(n-1)}{(n-2)(n-3)n_{t}(n_{t}-1)}\frac{n^{2}}{n_{t}^{2}(n-1)^{2}}$
$\frac{(4n-6)n^{2}}{(n-2)(n-3)n_{t}n_{c}(n-1)^{2}}$ (n-1)(n-2)(n-3) Y

  * •

Note: This table computes the entries of the tensor
$(\tilde{\mathbf{D}}^{N}\otimes\tilde{\mathbf{D}}^{N})\circ\mathbf{S}$ for
completely randomized experiments with the Neyman variance bound. WLOG, the
$i$th entry is fixed, and indices $j$, $k$, and $l$ vary. The first column
lists all possible $(i,j,k,l)$ patterns. The second column computes the case
for
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{1j},\mathds{R}_{1k}\mathds{R}_{1l})\frac{\tilde{\mathbf{D}}^{N}_{ij}\tilde{\mathbf{D}}^{N}_{kl}}{\pi_{ij}\pi_{lk}}$.
The third column computes the case for
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{1j},\mathds{R}_{0k}\mathds{R}_{0l})\frac{\tilde{\mathbf{D}}^{N}_{ij}\tilde{\mathbf{D}}^{N}_{kl}}{\pi_{ij}\pi_{lk}}$.
The fourth column counts the number of such patterns with $i$ fixed. The last
column confirms if the sum ($(|v_{111}|+|v_{100}|)$)$*$Count is of the order
$O(1)$. Note that with the Neyman bound, patterns involving
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{0j},\mathds{R}_{1k}\mathds{R}_{0l})$
and
$\mathbf{COV}(\mathds{R}_{1i}\mathds{R}_{1j},\mathds{R}_{1k}\mathds{R}_{0l})$
(the off-diagonal block) are multiplied with zeros, so we do not need to check
them here.

###  7.2 Model Adjustments in Network Experiments

In this section, we examine network experiments proposed in Aronow and Samii,
(2017) using our framework. Components of the experimental design include:

  * •

A finite population $U_{n}$ with units indexed by $i=1,...,n$. Each unit has a
trait vector $\xi_{i}\in\Xi_{n}$ (i.e., network connections) and a
pretreatment covariate vector $x_{i}\in\mathbf{R}^{p}$. Let $\Theta_{n}$
denote the set of traits.

  * •

An experimental design that randomly selects units into $M$ treatment values.
One realization of the assignment vector has the form
$Z=(Z_{1},...,Z_{n})\in\\{0,...,M-1\\}^{n}$. The distribution of the random
assignment vector $Z$, denoted as $P(Z)$, is known. Let
$\Omega_{n}\subset\\{0,...,M-1\\}^{n}$ denote the set of possible random
assignment vectors.

  * •

An exposure mapping that maps the assignment treatment vectors and a unit-
specific trait to an exposure value,
$F_{n}:\Omega_{n}\times\Xi_{n}\to\Delta_{n}$, where $\Delta_{n}$ denotes the
set of possible exposure values. This map is specified by the researcher
depending on the research questions at hand. $\Delta_{n}$ is usually specified
to be a finite set.414141In this setup, it is possible that the exposure
mappings are misspecified. See Aronow and Samii, (2017), Sävje, (2021), Sävje
et al., (2021), and Leung, (2022) for a discussion of estimation and
inferential theories in this context. We will proceed as if the exposure
mapping is correctly specified. Estimation and inference under misspecified
mappings do not immediately follow from our framework, although generalization
is possible under additional assumptions.

For one experiment, researchers randomly draw an assignment vector $Z$ and
observe the scalar outcomes $\\{Y_{i}(Z)\\}_{i=1}^{n}$. Notice that up to this
stage the outcome for unit $i$ has depended on the entire assignment vector.
Consistent estimation under unrestricted interference is deemed virtually
impossible Sävje et al., (2021). One strategy to alleviate the problem is to
restrict the interference patterns using exposure mappings. The potential
outcomes are assumed to be correctly indexed by the exposure value:

######  Assumption 15.

For $i=1,...,n$ and $Z,\tilde{Z}^{\prime}\in\Omega_{n}$,
$Y_{i}(Z)=Y_{i}(Z^{\prime})$ if $F_{n}(Z,\xi_{i})=F_{n}(\tilde{Z},\xi_{i})$.
$\Delta_{n}=\Delta$ for all $n$ and $\Delta$ is a finite set.

Assumption 15 implies that the exposures are ”effective treatments” as defined
in Manski, (2013). The assumption that $\Delta_{n}$ equals a finite set
$\Delta$ is a typical assumption made in the literature Aronow and Samii,
(2017). We note that this experimental setup is very general, and it can be
generalized to other settings in which the exposure mappings are not
necessarily mediated by a network.
Enumerate the element in $\Delta$ as $\\{1,...,k\\}$. With Assumption 15, one
can write the potential outcomes associated with unit $i$ as
$\left(y_{1i},...,y_{ki}\right)\in\mathbf{R}^{k}$. The assignment vector
associated with unit $i$ can be written as
$\left(\mathds{R}_{1i},...,\mathds{R}_{ki}\right)$. This maps the problem back
to our framework. Denote the probability of unit $i$ receiving treatment $k$
as $\pi_{i}(k)$ and the joint probability of units $i$ and $j$ receiving
treatments $k$ and $l$ as $\pi_{ij}(k,l)$. We make two assumptions, as in
Aronow and Samii, (2017), that are sufficient for consistent parameter
estimation, consistent plug-in variance bound estimation, and a local-
dependence central limit theorem. We need to use the concept of a dependency
graph.424242A dependency graph tailored to our setting is can be described as
follows. The dependency graph is a graph indexed by the random variables in
the treatment vector $\mathds{R}1_{\scriptscriptstyle{kn}}$. Consider two
subsets of the random variables $\mathds{R}_{A}$ and $\mathds{R}_{B}$. If
there is no arrow pointing from the variables in set $\mathds{R}_{A}$ to the
variables in set $\mathds{R}_{B}$ and vice versa, the variables in
$\mathds{R}_{A}$ are jointly independent of the variables in $\mathds{R}_{B}$.
For a detailed discussion, see Chen and Shao, (2004) and Chen et al., (2011).

######  Assumption 16.

For a positive $c_{\ref{A:networkprobability},1}\in(0,1)$,
$\pi_{i}(d_{k})>c_{\ref{A:networkprobability},1}>0$ uniformly in $i$ and in
$n$. For a positive $c_{\ref{A:networkprobability},2}\in(0,1)$,
$\pi_{ij}(k,l)>c_{\ref{A:networkprobability},2}>0$ uniformly in $i$ and in $n$
whenever $\pi_{ij}(k,l)\not=0$.

######  Assumption 17.

The dependency graph of the random assignment vector
$\mathds{R}1_{\scriptscriptstyle{kn}}$ has bounded degrees uniformly in $n$.

These two assumptions are commonly made in the literature. With these two
assumptions, one can show that the first-order design matrix $\mathbf{D}$ and
the AS bound matrix have bounded row norms.
$\tilde{\mathbf{D}}^{{\scriptscriptstyle{\textnormal{AS}}}}$ has a bounded row
norm as well. Furthermore, Assumption 9 is also satisfied by Proposition 6.2
in Aronow and Samii, (2017). To make the assumptions concrete, we provide an
example from Section 9 in Aronow and Samii, (2017) below.

######  Example 7.

Consider a situation where we observe $n$ units connected in undirected
networks. Each unit $i$ is associated with the trait $\theta_{i}$, which is
the $i$th row vector of the unnormalized adjacency matrix. The treatment
values are $\\{0,1\\}$, and the treatment assignment vector is denoted as
$Z\in\\{0,1\\}^{n}$ with $Z_{i}$ denoting the treatment assignment of unit
$i$. The exposure mapping is assumed to be:

$f(Z,\theta_{i})=\begin{cases}d_{11}(\text{Direct+Indirect Exposure}):\hskip
11.38109ptZ_{i}I(Z^{\prime}\theta_{i}>0)=1\\\ d_{10}(\text{Isolated Direct
Exposure}):\hskip 11.38109ptZ_{i}I(Z^{\prime}\theta_{i}=0)=1\\\
d_{01}(\text{Indirect Exposure}):\hskip
11.38109pt(1-Z_{i})I(Z^{\prime}\theta_{i}>0)=1\\\ d_{00}(\text{No
Exposure}):\hskip 11.38109pt(1-Z_{i})I(Z^{\prime}\theta_{i}=0)=1\\\
\end{cases}.$ (61)

The design is assumed to be a Bernoulli with probability $\frac{1}{2}$. The
network degree is bounded uniformly in $n$, and there are no isolated nodes.

One can confirm that Assumptions 16 and 17 are satisfied in this setup. The
assumption of bounded degrees is most appropriate, for example, when people
form networks mostly within a community and the number of people in a
community is bounded. Situations like this include schools and villages.434343
In some social network settings, there may exist nodes of very large degrees.
This may cause two problems: 1) the assignment probability becomes too small,
and/or 2) the dependence between exposure maps becomes too strong. In these
cases, one may wish to restrict the parameters of interest to subgroups (e.g.,
people with low degrees) and consider specific exposure mappings for precise
estimates (e.g., the sample average treatment effects condition on the event
that all the high-degree nodes are assigned to treatment).
Recall that $p$ is the dimension of the pretreatment covariates. We consider
the case of a same-slope adjustment and consider two types of imputation
functions: linear models $f^{a}(x_{i},\theta)=\gamma^{a}+x_{i}^{\prime}\beta$,
$a=1,...,k$ and logistic models
$f^{a}(x_{i},\theta)=\frac{\exp(\gamma^{a}+x_{i}^{\prime}\beta)}{1+\exp(\gamma^{a}+x_{i}^{\prime}\beta)}$,
$a=1,...,k$. For the QMLE-GR estimator with the linear model, the finite
population criterion and the sample equivalent are

$\mathcal{L}_{n}^{\textnormal{ln}}(\theta)=\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}(y_{ai}-\gamma^{a}-x_{i}^{\prime}\beta)^{2},\hskip
11.38109pt\hat{\mathcal{L}}_{n}^{\textnormal{ln}}(\theta)=-\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\pi_{ai}}(y_{ai}-\gamma^{a}-x_{i}^{\prime}\beta)^{2},$

where $\theta=(\\{\gamma^{a}\\}_{a=1}^{k},\beta)\in\mathbf{R}^{k+p}$. For the
QMLE-GR estimator with the logistic regression model, they are

$\mathcal{L}_{n}^{\textnormal{lg}}(\theta)=-\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\left[y_{ai}(\gamma^{a}+x_{i}^{\prime}\beta)-\log(1+\exp{(\gamma^{a}+x_{i}^{\prime}\beta)})\right]$

and

$\hat{\mathcal{L}}_{n}^{\textnormal{lg}}(\theta)=-\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\pi_{ai}}\left[y_{ai}(\gamma^{a}+x_{i}^{\prime}\beta)-\log(1+\exp{(\gamma^{a}+x_{i}^{\prime}\beta)})\right],$

where $\theta=(\\{\gamma^{a}\\}_{a=1}^{k},\beta)\in\mathbf{R}^{k+p}$. For the
Opt-GR estimators, we form the criterion using the setup in Section 6.3 with
$\Omega=\mathbf{D}$. We choose the AS bound for a bounding matrix. The plug-in
variance bound estimators are constructed using the formulae in Theorems 4, 5,
and 6. For the Opt-GR estimator with the logistic model, we need more moment
assumptions in order to satisfy Assumption 14-(i)-(b).

######  Assumption 18 (Bounded 8th moments).

For all $n$ and $z_{i}\in\\{x_{1i},...,x_{pi}\\}$,

$\frac{1}{n}\sum_{i=1}^{n}z_{i}^{8}<C_{\ref{A:Bounded8thmoment}}<\infty,$ (62)

where $C_{\ref{A:Bounded8thmoment}}$ is a finite constant.

The following theorems specialize Theorems 4, 5, and 6 in this network
experiment setting.

######  Theorem 7.

For linear models:

  1. 1.

Under Assumptions 1, 2, 16, and 17, Theorem 4 holds for the QMLE-GR estimator.
In addition, under Assumption 12-(ii), Theorem 5 holds for the No-harm-GR
estimator, and, under Assumption 19-(ii), Theorem 8 holds for the Opt-I GR
estimator.

  2. 2.

Under Assumptions 1, 13-(iv), 16, and 17, Theorem 6 holds for the Opt-GR
estimator.

For logistic models:

  1. 1.

Under Assumptions 1, 10-(i), (ii), (vii), 16, 17, and 18, Theorem 4 holds for
the QMLE-GR estimator. In addition, under Assumption 12-(ii), Theorem 5 holds
for the No-harm-GR estimator, and, under Assumption 19-(ii), Theorem 8 holds
for the Opt-I GR estimator.

  2. 2.

Under Assumptions 1, 13-(i), (ii), (iv), 16, 17, and 18, Theorem 6 holds for
the Opt-GR estimator.

######  Remark 11.

With a CLT under local dependence Chen and Shao, (2004), one can construct a
confidence interval using a normal approximation if the asymptotic variance
$\frac{1}{n}c^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(y-f(\theta_{n}))\hskip
1.0pt\mathbf{D}\hskip 1.0pt\mathbf{diag}(y-f(\theta_{n}))\mathbf{1}c$ is
uniformly bounded above and below by positive constants for large $n$. That
is, there exist positive constants $c$ and $C$ such that the asymptotic
variance is bounded between $c$ and $C$ uniformly for large $n$. The upper
bound is satisfied under our assumptions. The lower bound is a typical
assumption. If the assumption on the lower bound is violated, the estimator
will converge to the true parameter with a faster-than-$\sqrt{n}$ rate, and
the confidence interval using a normal approximation may not have the
asymptotically correct coverage.

######  Remark 12.

For the Opt-GR estimators with logit models, the compactness and well-
separation assumptions in A13-(i) is necessary. This is because the gradient
of a logit model contains the factor
$\frac{1}{\exp(\gamma^{a}+x_{i}\beta)}\frac{\exp(\gamma^{a}+x_{i}\beta)}{1+\exp(\gamma^{a}+x_{i}\beta)}$.
This weight will converge to 0 if $\gamma^{a}$ or an entry in $\beta$ escapes
to infinity and, as a result, the first-order condition will approach zero. In
theory, such a problem can be solved, for example, by redefining the criterion
that does not suffer from the weighting problem or explicitly regularizing the
parameters. We leave it as a next step to understand the theoretical and
practical solutions to this problem.

##  8 Examples

In this section, we demonstrate the performance of our method. We study a
network experiment based on the data in Cai et al., (2015). Section 7.1
describes the background and the dataset. Section 7.2 uses the first-order
design matrix to understand the strengths and weaknesses of different designs
in this setting. Section 7.3 describes our imputation and simulation
procedures. Section 7.4 discusses simulation results.

###  8.1 Background and Dataset

Cai et al., (2015) studies the influence of social networks on weather
insurance adoption. The paper’s experimental design features two village-level
randomizations and two household-level randomizations. Each village is first
randomly assigned to receive insurance price variations, and villages that do
not receive price variations are further randomized to receive different
default buying options. Within each village, households are randomly assigned
to four treatment arms: first-round simple, first-round intensive, second-
round simple, and second-round intensive. Households assigned to the second
round are further randomly assigned to receive different information about the
take-up decisions made by households in the first round.

The village-level price variation is used to measure the monetary equivalence
of social network effects. The village-level variation of default buying
options and the household-level variation of take-up information are used to
study the mechanism of social network effects. The household-level
randomization to different treatment arms is used to measure the social
network effects. Our simulation below builds on the household-level
randomization to four treatment arms and abstracts away from the other three
randomizations. The experimental design in the paper is a village-level
stratified randomization. The computer code used for randomization in Cai et
al., (2015) is not immediately available.444444A brief description of the
randomization procedure is given in footnote 8 of Cai et al., (2015). As a
result, we decide to construct the randomization procedures in the simulation
according to our own understanding of the randomization procedures. Any
implications we draw from the simulation results below should not be related
to the original paper.

The experiment is carried out in villages in rural China. The unit of the
experiment is a household, which typically consisted of 4 to 5 people. The
outcome of interest is whether insurance take-up occurs, a binary variable,
which is observed as households’ individual purchase decisions at the end of
each information session. The treatment variables are the treatment statuses
of households and those of the households that are nominated as friends.
Households’ pretreatment covariate data and their social network data are
collected from surveys. The pretreatment covariates include variables such as
demographics, rice production, income, and experiences with natural disasters.
The social network is collected through friend nominations: a social network
survey asked household heads ”to list five close friends, either within or
outside the village, with whom they most frequently discussed rice production
or financial issues.” For more detailed information, please see Section II-B
of Cai et al., (2015).

We briefly describe here the social networks used in the simulation. The
friend nomination graph is a directed graph with edges pointing from the
nominators to the nominees. The network has 4806 households and 41 disjoint
weak components.454545A component in a directed graph is called a weak
component if it is a connected component in the directed graph with directed
edges replaced by non-directed ones. The average outdegree is 3.5, and the
maximum indegree is 18.464646We note that this network is a subset of the
complete network in the dataset. In the complete network, almost all
households nominated 5 friends. For our simulation, we dropped households with
missing household size and rice production area information, which are the
variables used for stratification. The transformed network includes 4806
households and an average degree of 3.5. This number is similar to that of the
network used in the analysis in the paper by Cai et al., (2015), as some
friend households in the dataset are not assigned to any treatment arms. We
plot the second largest component of the network in Figure 1.474747We have
chosen the second largest component because the largest component is too large
to fit on a page to show finer details. The network has clear community
structures, as most friendship ties form within natural villages. We employ
the social network and pretreatment covariate information in the data. There
are 12 exposure mappings considered in the paper, which we include in Table
8.1.484848In Cai et al., (2015), comparing exposures 1 and 2 provides the
effect of the intensive information session. Comparing exposures 1 with 3, 1
with 4, 2 with 8, and 2 with 9 yields information about diffusion not captured
by the friendship network (Column 6 in Table 2 of the paper). Comparing
exposures 1, 3, and 4 with exposures 5, 6, and 7 and exposures 2, 8, and 9
with exposures 10, 11, and 12 gives information about social network effects
(Column 5 in Table 2 of the paper).

A List of Exposure Mappings Considered in Cai et al., (2015) Exposure Index
Description 1 In FRS 2 In FRI 3 In SRS, with no friends in the first round 4
In SRS, with friends in FRS and no friends in FRI 5 In SRS, with friends in
FRI 6 In SRS, with two friends in FRI 7 In SRS, with more than two friends in
FRI 8 In SRI, with no friends in the first round 9 In SRI, with friends in FRS
and no friends in FRI 10 In SRI, with one friend in FRI 11 In SRI, with two
friends in FRI 12 In SRI, with more than two friends in FRI

  * •

The abbreviations stand for the treatment arms. FRS: First Round Simple, FRI:
First Round Intensive, SRS: Second Round Simple, and SRI: Second Round
Intensive.

Figure 1: The Second Largest Component of the Network in the Simulation

Note: Figure 1 plots the second largest component used in the simulation. Each
node represents a household. Two nodes are connected if one of the pair
nominates another as a friend. Natural villages are groups of households and
are the basis of randomization.

###  8.2 Use $|||\mathbf{D}|||_{2}$ to understand different designs

In this section, we compare three designs using our proposed measure
$|||\mathbf{D}|||_{2}$ from Section 4. The purpose of this demonstration is to
show how to use the $|||\mathbf{D}|||_{2}$ measure to understand the relative
strengths and weaknesses of various designs. We note that this is not the only
measure researchers might use to understand different experimental designs.
There are other aspects of experimental designs that the measure
$|||\mathbf{D}|||_{2}$ does not capture, such as covariate balancing, However,
we believe $|||\mathbf{D}|||_{2}$ is a useful value for researchers to examine
to understand their designs.

We consider three experimental designs: 1) a finely stratified natural
village-level randomization, 2) a natural village-level randomization, and 3)
a Bernoulli design. We label the three designs as Design A, Design B and
Design C, respectively.

The finely stratified natural village-level randomization is inspired by the
experimental design in Cai et al., (2015). We partition households in each
natural village into four groups based on their household sizes and rice
production areas. For each natural village, we calculate the medians of these
two variables and classify households into four groups: LL (below median
household sizes and below median rice production areas), LH (below median
household sizes and above median rice production areas), HL (above median
household sizes and below median rice production areas), and HH (above median
household sizes and above median rice production areas). In some villages,
there are strata with fewer than four households. In this case, we merge each
stratum with a stratum of the same type from another village. We choose the
other village such that it belongs to the same connected component of the
friendship network as the natural village of the stratum to be merged.
Households in each stratum are then completely randomized to four treatment
arms. If the number of households is not a multiple of four, there are at most
three remainder households. We assign the first remainder household to Second
Round Intensive, the second remainder household (if it exists) to Second Round
Simple, and the third household (if it exists) to First Round Intensive.

For the natural village-level randomization, we first merge two villages with
fewer than 10 people each with one other village. With each natural village,
we randomly assign households to four arms with proportions of $\frac{1}{10}$
in FRS, $\frac{1}{10}$ in FRI, $\frac{2}{5}$ in SRS, and $\frac{2}{5}$ in SRI.
Practically, we repeat the vector $(4,3,2,1,4,3,4,3,4,3)$ for each village.
When the number of households is not a multiple of 10, we start from the left
of the vector until the remaining households are exhausted.

For Bernoulli designs, we randomly assign households to four treatment arms
with probability (1/4,1/4,1/4,1/4).

Table 8.2 below documents the $|||\mathbf{D}|||_{2}$ measure for comparing
exposure 1 with exposures 2-7 with the three different designs.494949In
general, for each experimental design, there are $12\times 11$
$|||\mathbf{D}|||_{2}$ values to compute when comparing effects between two
arms. For exposition purposes, we only include 6 such measures for simplicity.
One might suspect that the diagonal entries (i.e.,
$\text{{Var}}(\frac{\mathds{R}_{ai}}{{\pi_{ai}}})$) mainly contribute to the
value $|||\mathbf{D}|||_{2}$. To see whether this is true, we report in Table
8.2 the largest eigenvalues but with the diagonal entries set to zero. We
denote these alternative measures as $|||\mathbf{D}^{o}||_{2}$.

Comparing $|||\mathbf{D}|||_{2}$ for different designs Design A Design B
Design C Exposures 1vs2 7.33 13.17 4.03 Exposures 1vs3 $\infty$ 41.73 341.90
Exposures 1vs4 65.52 48.76 37.32 Exposures 1vs5 27.16 41.66 21.50 Exposures
1vs6 $\infty$ $\infty$ 81.07 Exposures 1vs7 $\infty$ $\infty$ 384.91

  * •

This table reports the largest eigenvalues of the first-order design matrices
for comparing effects between exposure 1 and other exposures under various
designs.

Comparing $|||\mathbf{D}^{o}|||_{2}$ for different designs Design A Design B
Design C Exposures 1vs2 6.00 11.01 1.03 Exposures 1vs3 $\infty$ 24.63 191.48
Exposures 1vs4 30.09 28.86 19.90 Exposures 1vs5 8.65 19.01 6.99 Exposures 1vs6
$\infty$ $\infty$ 19.49 Exposures 1vs7 $\infty$ $\infty$ 134.40

  * •

This table reports the largest eigenvalues of the first-order design matrices
(with diagonal entries set to 0) for comparing effects between exposure 1 and
other exposures under various designs.

Each entry in the table documents the $|||\mathbf{D}|||_{2}$ value for
comparing exposure 1 with another exposure under one of the three experimental
designs. We first notice that there are 5 $\infty$ signs in the table. These
$\infty$ signs correspond to cases where some units have 0 assignment
probabilities to some exposures. For example, for exposure 6 in the finely
stratified design, 37 people have 0 assignment probabilities and 80 people
have assignment probabilities smaller than 0.01. For exposure 3 in the finely
stratified design, 37 people have 0 assignment probabilities and 1019 people
have assignment probabilities smaller than 0.01.505050For exposure 7 in the
finely stratified design, 184 people have 0 assignment probabilities and 1847
people have probabilities smaller than 0.01. For exposure 6 in the stratified
design, 32 people have 0 assignment probabilities and 935 people have
probabilities smaller than 0.01. For exposure 7 in the village stratified
design, 724 people have 0 assignment probabilities and 3820 people have
probabilities smaller than 0.01. The 0 and small probabilities are caused by
small strata: if a unit and its friends are all in a small stratum, then the
unit may never experience some exposure states. For example, if a unit and its
5 friends constitute a stratum and at least one unit is assigned to the first
round, then the unit will never experience a state where none of its friends
are in the first round.

Secondly, we observe that the $|||\mathbf{D}^{o}|||_{2}$ values in Table 8.2
are roughly half the size of their counterparts in Table 8.2. This
demonstrates both the diagonal and off-diagonal elements contribute to in the
$|||\mathbf{D}|||_{2}$ values.

Finally, looking at the columns for the Bernoulli design, we know that
estimating the difference of effects between exposure 3 and exposure 1 and the
difference of effects between exposure 7 and exposure 1 could be difficult.
This is intuitive because the event that none of the 3$\sim$5 friends of a
unit are assigned to the first round and the event that more than 2 of a
unit’s friends are assigned to the first round intensive session both have a
low probability. In consequence, we will not have enough samples to accurately
measure these effects.

###  8.3 Simulation Design

We report simulation results for comparisons of exposures 1 and 2, exposures 1
and 3, exposures 1 and 4 with the Bernoulli design, exposures 1 and 4 with the
finely stratified village-level randomization, and exposures 1 and 4 with the
natural village-level randomization. We choose to report results for comparing
exposures 1 and 2, exposures 1 and 3, and exposures 1 and 4 with the Bernoulli
design to examine the performance of our estimators across different values of
$|||\mathbf{D}|||_{2}$. We choose to report results for comparing exposures 1
and 4 with three different designs to examine the performance of our
estimators across different designs with approximately similar
$|||\mathbf{D}|||_{2}$ values.

For each comparison of two exposures, we impute the potential outcomes in two
ways515151We only impute one set of potential comes when comparing exposures 1
and 2 in Design A. Design A with exposures 1 and 2 is a simple Bernoulli
design such that the estimators have the same expected performance with the
two types of imputation strategies discussed below. . In the first simulation
scenario, we impute the potential outcomes using a logistic model
$y_{ai}=\mathbf{I}\\{\beta_{a}+x_{i}^{\prime}\beta>\epsilon_{i}\\},\hskip
2.0pt\epsilon_{i}\sim\text{Logistic}(0,1)$ with exposure-specific intercepts
and the same slope across all arms.525252We fix one realization of the shock
vector $(\epsilon_{i})_{i=1}^{n}$ throughout our simulations. The pretreatment
covariates we use in the simulation are household head gender, household head
age, household size, area of rice production, risk aversion, perceived
disaster probability of next year, and household head literacy.535353The
coefficient values are taken from the estimation result of Table 1 Column 2
Cai et al., (2015). In this scenario and barring finite-sample issues, QMLE-
GR, No-harm-GR, Opt-GR, and Opt-I GR estimators are expected to work similarly
well and show improvement over the baseline HT estimator. In the second
simulation scenario, we impute the potential outcomes such that the Opt-GR
estimators will have a large efficiency gain over the QMLE-GR estimators.
Specifically, for each pair of average potential outcomes, e.g. exposure 1 and
exposure 3, we take the sum of the first 5 eigenvectors corresponding to the
largest eigenvalue of the matrix
$(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})^{\prime}\mathbf{diag}(c^{\prime}\mathbf{1})\mathbf{D}_{13}\mathbf{diag}(c^{\prime}\mathbf{1})(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})-(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{D}\mathbf{x})^{-1}\mathbf{x}^{\prime})^{\prime}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})\mathbf{D}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{D}\mathbf{x})^{-1}\mathbf{x}^{\prime})$,
add 0.5 onto each entry and round them to the nearest integer.545454With
$c=(-1,1)$,
$\frac{1}{n}y^{\prime}(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})^{\prime}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})\mathbf{D}_{13}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})y$
is the asymptotic variance of the WLS estimator for comparing exposures 1 and
3, and
$\frac{1}{n}y^{\prime}(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})^{\prime}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})\mathbf{D}_{13}\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})(\mathbf{i}_{2n}-\mathbf{x}(\mathbf{x}^{\prime}\mathbf{x})^{-1}\mathbf{x}^{\prime})y$
is the asymptotic variance of the Opt-GR estimator with a linear model. Our
choices of the potential-outcome vector make the efficiency gain from using
the Opt-GR estimator large.  This scenario is used to demonstrate the
performance of the No-harm GR, Opt-GR, and Opt-I GR estimators. More details
for data constructions, imputations, and implementations can be found in
Appendix H. We refer to the first simulation scenario as Sim-Impute and the
second simulation scenario Sim-Optimal.555555In conclusion, we compare 9 sets
of simulations (5 exposure comparisons and each with 2 datasets and 1 exposure
comparison with 1 dataset).

We compare the following 11 estimators:

  1. 1.

Horvitz-Thompson estimator (HT),

  2. 2.

Hajek estimator (HA),

  3. 3.

OLS estimator (OLS),

  4. 4.

WLS estimator with inverse probability weighting (WLS),

  5. 5.

No-harm estimator using the imputations of the WLS estimator (NH WLS),

  6. 6.

QMLE-GR estimator with a logit model (Logit),565656We set
$\omega_{ai}=\pi_{ai}$ where $\omega_{ai}$ is defined in equation (44). This
is to mimic the exercise where researchers estimate a logit model without any
weighting.

  7. 7.

No-harm estimator using the imputations of a logit model (NH Logit),

  8. 8.

Opt-GR estimator with a linear model (Opt Linear),

  9. 9.

Opt-GR estimator with a logit model (Opt Logit),

  10. 10.

Opt-I GR estimator with imputations using an OLS model (Opt-I OLS),

  11. 11.

Opt-I GR estimator with imputations using a logit model (Opt-I Logit).

All adjustment models have a separate intercept for each arm and the same
coefficients on the covariates (non-interacted models).

###  8.4 Simulation Results

We report simulation results for comparing 1) exposures 1 and 4 in Design A,
Scenario Sim-Impute, 2) exposures 1 and 4 in Design A, Scenario Sim-Optimal,
3) exposures 1 and 3 in Design C, Scenario Sim-Impute, and 4) exposures 1 and
3 in Design C, Scenario Sim-Optimal. Simulation results for other comparisons
and designs are reported in Appendix I.575757 We say two villages are
connected if one household from one village nominates a household from the
other village as a friend. The use of normal critical value for inference in
Design A and Design B is justified by the local dependence CLT Chen and Shao,
(2004). We consider the setup where all natural villages have bounded sizes
and each natural village is connected to a finite number of other natural
villages.

We note three points from these simulation results. Firstly, when one compares
across all tables, the variances of estimators increase as
$|||\mathbf{D}|||_{2}$ increases. The increase is visible by comparing, for
example, the variances of the WLS estimators in Table 8.4 and Table
8.4.585858The two simulation scenarios use the same set of potential outcomes.
However, this increase is not necessarily linear as a function of
$|||\mathbf{D}|||_{2}$ since the measure $|||\mathbf{D}|||_{2}$ considers only
the worst case scenario, as discussed at the end of Section 4.1.

Secondly, the Opt-GR estimators (Opt Linear and Opt Logit) bring variance
reductions but also face bias-variance trade-offs in the finite sample. For
example, in Table 8.4, the Opt-GR linear estimator has a variance around 26%
percent lower than that of the WLS estimator. However, the Opt-GR estimators
also incur a finite sample bias. The performance of the Opt-GR estimators
deteriorate as $|||\mathbf{D}|||$ gets larger, as can be seen in Table I. This
phenomenon happens because the case considered in Table I is a case where the
imputations are weighted by small assignment probabilities. As a consequence,
the bias and variance in the estimated coefficients are magnified, especially
when the imputation functions are linear and contain multiple covariates.

Thirdly, in our simulations, the Opt-I OLS and Opt-I Logit estimators perform
reasonably well in all cases. The Opt-I OLS and Opt-I Logit estimators are
less efficient compared with the Opt Linear and Opt Logit estimators in terms
of theoretical asymptotic variance, but the loss of efficiency appears to be
small and the Opt-I OLS and Opt-I Logit estimators have better finite sample
performance. Taken together and based on our simulation results, we consider
the Opt-I OLS and Opt-I Logit as viable alternatives for Opt Linear and Opt
Logit in many practical settings.

Simulation Results for Comparing Exposures 1 and 4 in Design A and Scenario
Sim-Impute, $\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}=0.12$ HT HA OLS WLS NH
(WLS) Logit NH (Logit) OC (Linear) OC (Logit) Opt-I (OLS) Opt-I (Logit) Bias2
$\times$ N 0.00 0.00 0.01 0.00 0.01 0.00 0.01 0.42 0.36 0.04 0.04 Variance
$\times$ N 6.49 4.39 4.12 4.41 4.39 4.41 4.39 4.40 4.36 4.37 4.37 Mean Squared
Error $\times$ N 6.49 4.39 4.13 4.41 4.40 4.41 4.40 4.83 4.72 4.41 4.40
Estimated Variance Bound $\times$ N 7.51 5.21 4.89 5.13 5.14 5.13 5.19 5.21
5.26 5.21 5.21 95% Normal CI Coverage 0.96 0.96 0.96 0.96 0.96 0.96 0.96 0.95
0.95 0.96 0.96 Theoretical Asy. Variance $\times$ N 6.53 4.44 4.15 4.44 4.44
4.44 4.44 4.44 4.44 4.44 4.44 Theoretical Asy. Variance Bound $\times$ N 7.50
5.22 4.94 5.22 5.22 5.23 5.23 5.22 5.22 5.22 5.22

  * •

Table 8.4 reports simulation results for comparing exposures 1 and 4 in Design
A in the Sim-Impute scenario. N=4509 is the sample size. The number of
simulations is 3000. HT refers to the Horvitz-Thompson estimator, HA refers to
the Hajek estimator, OLS refers to the ordinary least square estimator, WLS
refers to the inverse-assignment-probability weighted least square estimator,
NH (WLS) refers to the No-harm estimators using the imputations of the WLS
estimator, Logit refers to the QMLE estimator with a logit model, NH (Logit)
refers to the No-harm estimator using the imputations of a logit model, OC
(Linear) refers to the Opt-GR estimator with a linear model and OC (Logit)
refers to the Opt-GR estimator with a logit model, Opt-I OLS refers to the
Opt-I GR estimator with an OLS model, and Opt-I Logit refers to the Opt-I GR
estimator with a logit model.

  * •

The row Bias2 $\times$ N reports the squared biases of the estimators,
multiplied by the sample size. The row Variance $\times$ N reports the
variances of the estimators, multiplied by the sample size. The row Mean
Squared Error $\times$ N reports the MSEs of the estimators, multiplied by the
sample size. The row Estimated Variance Bound $\times$ N reports the averaged
estimates of the variance bound estimators, multiplied by the sample size. The
row 95% Normal CI Coverage reports the empirical coverage rates of nominal 95
percent confidence intervals. The row True Asy. Variance reports the
theoretical asymptotic variances. The row True Asy. Variance Bound reports the
theoretical asymptotic variance bounds.

Simulation Results for Comparing Exposures 1 and 4 in Design A and Scenario
Sim-Optimal, $\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}=0.12$ HT HA OLS WLS NH
(WLS) Logit NH (Logit) OC (Linear) OC (Logit) Opt-I (OLS) Opt-I (Logit) Bias2
$\times$ N 0.00 0.00 0.04 0.00 0.02 0.00 0.01 0.69 0.52 0.05 0.05 Variance
$\times$ N 9.47 5.24 7.43 5.69 5.68 8.27 7.50 4.58 4.24 4.16 4.13 Mean Squared
Error $\times$ N 9.47 5.24 7.46 5.69 5.70 8.27 7.51 5.27 4.77 4.21 4.18
Estimated Variance Bound $\times$ N 10.62 6.20 8.16 6.44 6.45 6.46 8.52 5.09
5.07 5.22 5.15 95% Normal CI Coverage 0.96 0.97 0.96 0.96 0.96 0.92 0.96 0.95
0.96 0.97 0.97 Theoretical Asy. Variance$\times$ N 9.21 5.14 7.26 5.45 5.45
8.36 7.48 3.88 3.68 3.89 3.73 Theoretical Asy. Variance Bound $\times$ N 10.62
6.23 8.19 6.53 6.52 9.38 8.57 5.15 5.00 5.49 5.03

  * •

Table 8.4 reports simulation results for comparing exposures 1 and 4 in Design
A in the Sim-Optimal scenario. N=4509 is the sample size. The number of
simulations is 3000. HT refers to the Horvitz-Thompson estimator, HA refers to
the Hajek estimator, OLS refers to the ordinary least square estimator, WLS
refers to the inverse-assignment-probability weighted least square estimator,
NH (WLS) refers to the No-harm estimators using the imputations of the WLS
estimator, Logit refers to the QMLE estimator with a logit model, NH (Logit)
refers to the No-harm estimator using the imputations of a logit model, OC
(Linear) refers to the Opt-GR estimator with a linear model and OC (Logit)
refers to the Opt-GR estimator with a logit model, Opt-I OLS refers to the
Opt-I GR estimator with an OLS model, and Opt-I Logit refers to the Opt-I GR
estimator with a logit model.

  * •

The row Bias2 $\times$ N reports the squared biases of the estimators,
multiplied by the sample size. The row Variance $\times$ N reports the
variances of the estimators, multiplied by the sample size. The row Mean
Squared Error $\times$ N reports the MSEs of the estimators, multiplied by the
sample size. The row Estimated Variance Bound $\times$ N reports the averaged
estimates of the variance bound estimators, multiplied by the sample size. The
row 95% Normal CI Coverage reports the empirical coverage rates of nominal 95
percent confidence intervals. The row True Asy. Variance reports the
theoretical asymptotic variances. The row True Asy. Variance Bound reports the
theoretical asymptotic variance bounds.

Simulation Results for Comparing Exposures 1 and 3 in Design C and Scenario
Sim-Impute, $\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}=0.27$ HT HA OLS WLS NH
(WLS) Logit NH (Logit) OC (Linear) OC (Logit) Opt-I (OLS) Opt-I (Logit) Bias2
$\times$ N 0.00 0.01 0.33 0.01 0.77 0.01 0.69 19.88 14.51 1.88 1.89 Variance
$\times$ N 32.05 16.50 8.76 16.48 17.16 16.52 16.83 22.96 20.02 17.34 17.38
Mean Squared Error $\times$ N 32.05 16.51 9.09 16.48 17.93 16.53 17.51 42.84
34.53 19.22 19.28 Estimated Variance Bound $\times$ N 35.63 16.95 10.16 15.25
16.25 15.22 17.68 18.80 28.20 18.07 18.08 95% Normal CI Coverage 0.94 0.93
0.95 0.92 0.92 0.93 0.93 0.80 0.90 0.92 0.92 Theoretical Asy. Variance
$\times$ N 31.57 16.12 8.59 16.11 16.11 16.15 16.14 16.10 16.10 16.11 16.11
Theoretical Asy. Variance Bound $\times$ N 35.66 17.64 10.34 17.64 17.64 17.68
17.66 17.63 17.63 17.64 17.64

  * •

Table 8.4 reports simulation results for comparing exposures 1 and 3 in Design
C in the Sim-Impute scenario. N=4509 is the sample size. The number of
simulations is 3000. HT refers to the Horvitz-Thompson estimator, HA refers to
the Hajek estimator, OLS refers to the ordinary least square estimator, WLS
refers to the inverse-assignment-probability weighted least square estimator,
NH (WLS) refers to the No-harm estimators using the imputations of the WLS
estimator, Logit refers to the QMLE estimator with a logit model, NH (Logit)
refers to the No-harm estimator using the imputations of a logit model, OC
(Linear) refers to the Opt-GR estimator with a linear model and OC (Logit)
refers to the Opt-GR estimator with a logit model, Opt-I OLS refers to the
Opt-I GR estimator with an OLS model, and Opt-I Logit refers to the Opt-I GR
estimator with a logit model.

  * •

The row Bias2 $\times$ N reports the squared biases of the estimators,
multiplied by the sample size. The row Variance $\times$ N reports the
variances of the estimators, multiplied by the sample size. The row Mean
Squared Error $\times$ N reports the MSEs of the estimators, multiplied by the
sample size. The row Estimated Variance Bound $\times$ N reports the averaged
estimates of the variance bound estimators, multiplied by the sample size. The
row 95% Normal CI Coverage reports the empirical coverage rates of nominal 95
percent confidence intervals. The row True Asy. Variance reports the
theoretical asymptotic variances. The row True Asy. Variance Bound reports the
theoretical asymptotic variance bounds.

Simulation Results for Comparing Exposures 1 and 3 in Design C and Scenario
Sim-Optimal, $\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}=0.27$ HT HA OLS WLS NH
(WLS) Logit NH (Logit) OC (Linear) OC (Logit) Opt-I (OLS) Opt-I (Logit) Bias2
$\times$ N 0.01 0.00 21.86 0.18 0.90 0.01 0.35 38.87 17.04 5.15 4.77 Variance
$\times$ N 62.49 19.81 13.07 27.49 28.89 49.03 41.95 34.62 27.88 20.44 20.40
Mean Squared Error $\times$ N 62.50 19.81 34.93 27.67 29.79 49.04 42.30 73.50
44.93 25.59 25.17 Estimated Variance Bound $\times$ N 71.40 21.35 14.66 25.55
27.70 49.38 48.15 16.95 21.15 23.57 22.99 95% Normal CI Coverage 0.95 0.94
0.77 0.92 0.92 0.89 0.94 0.68 0.78 0.94 0.94 Theoretical Asy. Variance
$\times$ N 60.71 19.59 13.05 26.46 26.30 50.32 44.56 12.99 11.73 13.09 12.81
Theoretical Asy. Variance Bound $\times$ N 71.10 22.50 15.17 29.89 29.74 55.41
50.20 15.16 13.71 15.47 15.11

  * •

Table 8.4 reports simulation results for comparing exposures 1 and 3 in Design
C in the Sim-Optimal scenario. N=4509 is the sample size. The number of
simulations is 3000. HT refers to the Horvitz-Thompson estimator, HA refers to
the Hajek estimator, OLS refers to the ordinary least square estimator, WLS
refers to the inverse-assignment-probability weighted least square estimator,
NH (WLS) refers to the No-harm estimators using the imputations of the WLS
estimator, Logit refers to the QMLE estimator with a logit model, NH (Logit)
refers to the No-harm estimator using the imputations of a logit model, OC
(Linear) refers to the Opt-GR estimator with a linear model and OC (Logit)
refers to the Opt-GR estimator with a logit model, Opt-I OLS refers to the
Opt-I GR estimator with an OLS model, and Opt-I Logit refers to the Opt-I GR
estimator with a logit model.

  * •

The row Bias2 $\times$ N reports the squared biases of the estimators,
multiplied by the sample size. The row Variance $\times$ N reports the
variances of the estimators, multiplied by the sample size. The row Mean
Squared Error $\times$ N reports the MSEs of the estimators, multiplied by the
sample size. The row Estimated Variance Bound $\times$ N reports the averaged
estimates of the variance bound estimators, multiplied by the sample size. The
row 95% Normal CI Coverage reports the empirical coverage rates of nominal 95
percent confidence intervals. The row True Asy. Variance reports the
theoretical asymptotic variances. The row True Asy. Variance Bound reports the
theoretical asymptotic variance bounds.

Online Appendix

##  Appendix A Additional Results in Section 4

We prove Lemma 4.1 (here as Lemma A.1).

######  Lemma A.1.

Under Assumption 3, the WLS estimator with
$\mathds{W}^{\scriptscriptstyle{\textnormal{WLS}}}=\begin{bmatrix}\mathbf{i}_{k}|\mathbf{0}_{k\times
p}\end{bmatrix}\left(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}\right)^{+}\mathbf{x}^{\prime}\mathbf{m}$
is equivalent to a completely imputed estimator with
$\mathds{Q}=\mathbf{0}_{k\times kn}$ and
$\mathds{P}=\mathbf{x}\left(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}\right)^{+}\mathbf{x}^{\prime}\mathbf{m}$.

###### Proof.

First notice the Completely Imputed estimator with
$\mathds{Q}=\mathbf{0}_{k\times kn}$ and
$\mathds{P}=\mathbf{x}\left(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}\right)^{+}\mathbf{x}^{\prime}\mathbf{m}$
can be written as

$\displaystyle\widehat{\delta}^{\scriptscriptstyle{\textnormal{CI}}}=$
$\displaystyle\mathds{W}^{\scriptscriptstyle{\textnormal{CI}}}\mathds{R}y=\left(\mathbf{1}^{\prime}\mathbf{1}\right)\mathbf{1}^{\prime}\mathbf{x}\left(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}\right)^{-1}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y=\left(\mathbf{1}^{\prime}\mathbf{1}\right)\mathbf{1}^{\prime}\mathbf{x}\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=\left(\mathbf{1}^{\prime}\mathbf{1}\right)\mathbf{1}^{\prime}\widehat{y}$

Note by Assumption 3, we have
$\frac{1}{n}\mathbf{1}^{\prime}\mathbf{x}=\begin{bmatrix}\mathbf{i}_{k}|\mathbf{0}_{k\times
p}\end{bmatrix}\in\mathbf{R}^{k\times(k+p)}.$ Thus
$\widehat{\delta}^{\scriptscriptstyle{\textnormal{CI}}}=\frac{1}{n}\mathbf{1}^{\prime}\mathbf{x}\
\widehat{b}_{n}^{\scriptscriptstyle{\textnormal{WLS}}}=\begin{bmatrix}\mathbf{i}_{k}|\mathbf{0}_{k\times
p}\end{bmatrix}\widehat{b}_{n}^{\scriptscriptstyle{\textnormal{WLS}}}=\hat{\mu}_{n}^{{\scriptscriptstyle{\textnormal{WLS}}}}$,
where $\hat{\mu}_{n}^{{\scriptscriptstyle{\textnormal{WLS}}}}$ is WLS
estimator. This proves Lemma A.1. ∎

##  Appendix B Additional Results in Section 6

This section contains additional results for the Opt-I estimators. Define
$X^{c}=\mathbf{diag}(c^{\prime}\mathbf{1}^{\prime})\begin{bmatrix}\mathbf{1}&f(\theta_{n})\end{bmatrix}\in\mathbf{R}^{kn\times(p+1)}$
and
$X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}=\begin{bmatrix}\mathbf{1}&f(\theta_{n})\end{bmatrix}\in\mathbf{R}^{kn\times(p+1)}$.

######  Assumption 19.

  1. (i)

$|||\mathbf{D}|||_{1}=O(1)$.

  2. (ii)

There exists a positive integer $N$ and a positive constant $c_{\ref{A:OCI}}$
such that

$\frac{1}{n}(X^{c})^{\prime}\mathbf{D}X^{c}\succeq
c_{\ref{A:OCI}}1_{\scriptscriptstyle{k+1}},$ (63)

uniformly for all $n\geq N$.

######  Theorem 8.

Define $\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}$ and
$\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}$ as in
Algorithm 4. Further define $\beta$ as
$\beta=\begin{bmatrix}\beta^{1},...,\beta^{k+1}\end{bmatrix}=\left((X^{c})^{\prime}\mathbf{D}X^{c}\right)^{-1}X^{c}\mathbf{D}y$
and
$\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}=\frac{1}{n}\sum_{i=1}^{n}(\beta^{a}+\beta^{k+1}f^{a}(x_{i},\theta_{n}))+\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}(y_{ai}-\beta^{a}-\beta^{k+1}f^{a}(x_{i},\theta_{n}))$
for $a=1,...,k$, where
$\theta_{n}\equiv\arg\min_{\theta\in\Theta}\mathcal{L}_{n}(\theta)$ and
$\mathcal{L}_{n}(\theta)$ is defined in (44). Under Assumptions 1, 6, 9, 10,
11 and 19, we have i) $\hat{\beta}-\beta=O_{p}(\frac{1}{\sqrt{n}})$, ii)
$\widehat{\mu}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}_{n,a}-\hat{\mu}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}_{n,a}=o_{p}(\frac{1}{\sqrt{n}})$,
and iii)
$\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}}}-\gamma_{n,c}=O_{p}(\frac{1}{\sqrt{n}})$.
Define
$\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}=\sum_{a=1}^{k}c_{a}\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}$.
The variance of
$\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}$
can be expressed as

$\text{{Var}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})=\frac{1}{n^{2}}c^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\beta)\hskip
1.0pt\mathbf{D}\hskip
1.0pt\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\beta)\mathbf{1}c.$

Define the variance bound
$\tilde{\text{{Var}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{NOH}}},{\scriptscriptstyle{\textnormal{L}}}})=\frac{1}{n^{2}}c^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\hat{\beta})\tilde{\mathbf{D}}\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\hat{\beta})\mathbf{1}c\in\mathbf{R}$,
with an identified variance bound matrix $\tilde{\mathbf{D}}$. The plug-in
variance-bound estimator
$\widehat{\tilde{\text{{Var}}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})=\frac{1}{n^{2}}\mathbf{1}^{\prime}\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\beta))\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\mathbf{diag}(y-X^{{\scriptscriptstyle{\textnormal{Opt-I}}}}\beta))\mathbf{1}\in\mathbf{R}$
is consistent:
$n(\hat{\tilde{\text{{Var}}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})-\tilde{\text{{Var}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}))\overset{p}{\to}0$.
If there exists a positive constant $c_{\ref{Thm:OCI}}$ such that
$n\tilde{\text{{Var}}}(\widehat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})\geq
c_{\ref{Thm:OCI}}\mathbf{i}_{k}$ uniformly for large $n$, then
$\frac{\widehat{\tilde{\text{{Var}}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})}{\tilde{\text{{Var}}}(\hat{\gamma}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}})}\overset{p}{\to}1$.

##  Appendix C Auxiliary Lemmas

This section proves several auxiliary lemmas instrumental in the proofs below.
Note Assumptions 1′, 6′ and 9′ are stated in Section J.

######  Lemma C.1 (Well-behaved WLS Design Matrix).

Under Assumption 2 and Assumption 5, there exists a $N$ such that for all
$n\geq N$,
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\pi\mathbf{x}\in\mathbf{R}^{(k+p)\times(k+p)}$
is invertible. Moreover,
$\lambda_{\min}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x})$
is bounded away from 0 uniformly for $n\geq N$.

###### Proof.

Let $\mathbf{M}$ be the limiting matrix in Assumption 2 and set a
$\epsilon<\frac{\lambda_{\min}(\mathbf{M})}{2}$. By Assumption 2, there exists
a $N$ such that for all $n\geq N$,
$||\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x}-\mathbf{M}||_{2}<\epsilon$. We
have
$|\lambda_{\min}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x})-\lambda_{\min}(\mathbf{M})|<\epsilon$
by Theorem 6.3.8 in Horn and Johnson, (2012). This implies
$\lambda_{\min}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x})\geq\frac{\lambda_{\min}(\mathbf{M})}{2}$
for all $n\geq N$.
Since $\mathbf{m}$ and $\pi$ are positive diagonal matrices,
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\pi\mathbf{x}\in\mathbf{R}^{(k+p)\times(k+p)}$
is a positive semidefinte matrix. We only need to check its smallest
eigenvalue. For any $t\in\mathbf{R}^{k+p}$

$\displaystyle
t^{\prime}\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}t$
$\displaystyle\geq
c_{5}\frac{1}{n}||\mathbf{x}t||_{2}^{2}=c_{5}\times(t^{\prime}\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x}t)>\frac{c_{5}\lambda_{\min}(\mathbf{M})}{2}||t||_{2}^{2}>0$

where $c_{5}$ is the constant from Assumption 5. This inequality holds for all
$n\geq N$, proving the statement. ∎

######  Lemma C.2 (Bounded WLS Coefficients).

Let
$b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=(\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{+}(\mathbf{x}\mathbf{m}\boldsymbol{\pi}y)$.
Under Assumptions 1, 2 and 5, or Assumptions 1′, 2 and 5,
$||b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}=O(1)$

###### Proof.

We showed in Lemma C.1 that
$\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}$ is invertible for
large $n$, so we shall assume
$b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=(\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{-1}(\mathbf{x}\mathbf{m}\boldsymbol{\pi}y)$.
For a $\epsilon>0$ and $n$ large,

$\displaystyle\begin{split}||\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y||_{2}^{2}=&\frac{1}{n^{2}}y^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}^{\prime}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y\\\
&\leq\lambda_{\max}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x})\times\frac{1}{n}||\mathbf{m}\pi
y||_{2}^{2}\\\
&\leq\lambda_{\max}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x})\times
C_{5}^{2}\frac{1}{n}||y||_{2}^{2}\\\ &\leq(\lambda_{\max}(M)+\epsilon)\times
C_{5}^{2}\frac{1}{n}||y||_{2}^{2}<\infty,\end{split}$ (64)

where $C_{5}$ is the constant from Assumption 5. For the last inequality we
used Assumption 2 to bound
$\lambda_{\max}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{x})$, and Assumption 1
(or Assumption 1′) to bound $\frac{1}{n}||y||_{2}^{2}$. Then,

$\displaystyle\begin{split}b^{{\scriptscriptstyle{\textnormal{WLS}}}^{\prime}}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}&=(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y)^{\prime}(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{-1}(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{-1}(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y)\\\
&\leq\frac{1}{\lambda_{\min}^{2}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\pi\mathbf{x})}\times||\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y||_{2}^{2}=O(1)\end{split}$

by Lemma C.1 and (64). ∎

Let $z_{n}\in\mathbf{R}^{kn}$. We define the HT estimator:

$\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}z_{n}\in\mathbf{R}^{k}$

######  Lemma C.3 (Rate of the HT Estimators).

If $\frac{1}{n}||z_{n}||_{2}^{2}=O(1)$,
$nV(\widehat{\delta}^{HT})=O(|||\mathbf{D}|||_{2})$.

###### Proof.

Notice we can write:

$\displaystyle\begin{split}\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\mathbf{diag}(z_{n})1_{\scriptscriptstyle{kn}}=\frac{1}{n}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\boldsymbol{\pi}^{-1}\mathds{R}1_{\scriptscriptstyle{kn}}\end{split}$

We have for any $t\in\mathbf{R}^{k}$

$\displaystyle\begin{split}nV(t^{\prime}\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})&=nt^{\prime}\frac{1}{n}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\text{{Var}}(\boldsymbol{\pi}^{-1}\mathds{R}1_{\scriptscriptstyle{kn}})\mathbf{diag}(z_{n})\mathbf{1}t\frac{1}{n}\\\
&=\frac{1}{n}(t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z_{n}))\mathbf{D}(\mathbf{diag}(z_{n})\mathbf{1}t)\leq|||\mathbf{D}|||_{2}\frac{1}{n}||t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})||_{2}^{2}\leq|||\mathbf{D}|||_{2}\frac{1}{n}t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\mathbf{diag}(z_{n})\mathbf{1}t\\\
&\leq|||\mathbf{D}|||_{2}||t||_{2}^{2}\frac{1}{n}\lambda_{\max}(\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\mathbf{diag}(z_{n})\mathbf{1})\leq|||\mathbf{D}|||_{2}||t||_{2}^{2}\frac{1}{n}\mathbf{Tr}(\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\mathbf{diag}(z_{n})\mathbf{1})\\\
&=|||\mathbf{D}|||_{2}||t||_{2}^{2}\frac{1}{n}||z_{n}||_{2}^{2}=||t||_{2}^{2}|||\mathbf{D}|||_{2}\times
O(1).\end{split}$

where the last line holds by our assumption on $\frac{1}{n}||z_{n}||_{2}^{2}$.
We conclude the largest eigenvalue of the positive semi-definite matrix
$V(\hat{\delta}^{HT})$ is of order $O(\frac{|||\mathbf{D}|||_{2}}{n})$.
$||nV(\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})||_{2}\leq\sqrt{k}|||nV(\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})|||_{2}=O(|||\mathbf{D}|||_{2})$,
proving the statement. ∎

######  Lemma C.4.

If $||z_{n}||_{\infty}=O(1)$ for all $n$,
$nV(\delta^{HT}_{c})=O(\frac{1}{n}||\mathbf{D}||_{1})$

###### Proof.

Notice we can write:

$\displaystyle\begin{split}\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\mathbf{diag}(z_{n})1_{\scriptscriptstyle{kn}}=\frac{1}{n}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\boldsymbol{\pi}^{-1}\mathds{R}1_{\scriptscriptstyle{kn}}\end{split}$

We have for any $t\in\mathbf{R}^{k}$

$\displaystyle\begin{split}nV(t^{\prime}\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})&=nt^{\prime}\frac{1}{n}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})\text{{Var}}(\boldsymbol{\pi}^{-1}\mathds{R}1_{\scriptscriptstyle{kn}})\mathbf{diag}(z_{n})\mathbf{1}\frac{1}{n}t\\\
&=\frac{1}{n}(t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z_{n}))\mathbf{D}(\mathbf{diag}(z_{n})\mathbf{1}t)\leq\frac{1}{n}||\mathbf{D}||_{1}\times||t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z_{n})||_{\infty}^{2}\leq\frac{1}{n}||\mathbf{D}||_{1}\times||t||_{\infty}^{2}\times
k||z_{n}||_{\infty}^{2}\\\
&=\frac{1}{n}||\mathbf{D}||_{1}\times||t||_{2}^{2}\times O(1)\end{split}$

where the first and second inequality is algebraic and the last line holds by
our assumption on $||z_{n}||_{\infty}$. We conclude the largest eigenvalue of
the positive semi-definite matrix $V(\hat{\delta}^{HT})$ is of order
$O(\frac{\frac{1}{n}||\mathbf{D}||_{1}}{n})$.
$||nV(\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})||_{2}\leq\sqrt{k}|||nV(\hat{\delta}^{{\scriptscriptstyle{\textnormal{HT}}}})|||_{2}=O(\frac{1}{n}||\mathbf{D}||_{1})$,
proving the statement. ∎

######  Lemma C.5.

Under Assumptions 1, 2 and 5, and if $\frac{|||\mathbf{D}|||_{2}}{n}=o(1)$,

$\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=O_{p}(\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}).$

Under Assumptions 1′, 2, and 5 and if
$\frac{\frac{1}{n}||\mathbf{D}||_{1}}{n}=o(1)$,

$\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=O_{p}(\sqrt{\frac{\frac{1}{n}||\mathbf{D}||_{1}}{n}}).$

###### Proof.

We prove the first claim. The second one is proved analogously. We first show
the ”numerator” vector $\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y$
is $\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}$-consistent for
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}y$.
Notice $\frac{1}{n}\mathbf{x}^{\prime}\mathds{R}y$ is an unbiased estimator
for $\frac{1}{n}\mathbf{x}^{\prime}\boldsymbol{\pi}y$. We only need to show
that the variance is of the order $O(\frac{|||\mathbf{D}|||_{2}}{n})$. Let
$\mathbf{x}_{i}$ be the column vector created from the $i$th column of
$\mathbf{x}$. Then the $i$th element of
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y$ can be written,

$\displaystyle\\{\frac{1}{n}\mathbf{x}\mathbf{m}\mathds{R}y\\}_{i}=$
$\displaystyle\frac{1}{n}\mathbf{x}_{i}^{\prime}\mathbf{m}\mathds{R}y=\frac{1}{n}1_{\scriptscriptstyle{kn}}^{\prime}\mathbf{diag}(\mathbf{x}_{i})\mathbf{m}\mathds{R}y=\frac{1}{n}1_{\scriptscriptstyle{kn}}^{\prime}\mathds{R}\mathbf{m}\mathbf{diag}(\mathbf{x}_{i})y$
$\displaystyle=$ $\displaystyle
1_{\scriptscriptstyle{k}}^{\prime}\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\boldsymbol{\pi}\mathbf{m}\mathbf{diag}(\mathbf{x}_{i})y=1_{\scriptscriptstyle{k}}^{\prime}\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\boldsymbol{\pi}\mathbf{m}(\mathbf{x}_{i}\circ
y)$

Under Assumption 1, we have

$\displaystyle\frac{1}{n}||\boldsymbol{\pi}\mathbf{m}(\mathbf{x}_{i}\circ
y)||_{2}^{2}$ $\displaystyle\leq
C_{5}^{2}\times\frac{1}{n}||(\mathbf{x}_{i}\circ y)||_{2}^{2}\leq
C_{5}^{2}\sqrt{\frac{1}{n}||\mathbf{x}_{i}||_{4}^{4}}\times\sqrt{\frac{1}{n}||y||_{4}^{4}}<\infty,$

where the first inequality is by Assumption 5 and the second line is by the
Cauchy-Schwartz inequality and Assumption 1. Thus by Lemma C.3, we prove the
$\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}$ consistency of the numerator vector.
Thus we have
$\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}||\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y-\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y||_{2}=O_{p}(1)$

Similarly, the $(i,j)$ element of the WLS ”denominator” matrix, can be written
as:

$\displaystyle\\{\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}\\}_{ij}=$
$\displaystyle
1_{\scriptscriptstyle{k}}^{\prime}\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\boldsymbol{\pi}\mathbf{m}(\mathbf{x}_{i}\circ\mathbf{x}_{j}).$

Following the same argument as above, we can show
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}$ is
$\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}$-consistent for
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}$. Note by
using the Frobenius norm and the Weyl’s inequality, the smallest eigenvalue of
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}$ converges to
the smallest eigenvalue of
$\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}$ in
probability. Thus for a small
$0<\epsilon<\frac{1}{2}\lambda_{\min}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x})$
and with Lemma C.1,
$\mathbf{P}(\lambda_{\min}(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x})>\epsilon)\to
1$. Thus we have
$\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}||(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x})^{+}-(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{-1}||_{2}=O_{p}(1)$.

Finally note the algebraic decomposition that for
$\widehat{A},A\in\mathbf{R}^{k_{1}\times k_{2}}$ and
$\widehat{B},B\in\mathbf{R}^{k_{2}\times k_{3}}$

$\widehat{A}\widehat{B}-AB=(\widehat{A}-A)(\widehat{B}-B)+(\widehat{A}-A)B+A(\widehat{B}-B),$

Let
$\widehat{A}=(\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x})^{+}$,
$\widehat{B}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y$,
$A=(\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}\mathbf{x})^{-1}$ and
$B=\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y$, we have:

$\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}}||\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||=O_{p}(1)$

∎

######  Lemma C.6.

Consider a sequence of symmetric matrices
$\widehat{A}_{n}\in\mathbf{R}^{k\times k},n=1,2..,$ and a symmetric matrix
$A_{n}\in\mathbf{R}^{k\times k}$. If for every $t\in\mathbf{R}^{k}$,
$t^{\prime}\widehat{A}_{n}t-t^{\prime}A_{n}t\overset{p}{\to}0$, then
$\widehat{A}_{n}-A_{n}\overset{p}{\to}\mathbf{0}_{k\times k}$.

###### Proof.

Using the standard basis vectors in $\mathbf{R}^{k}$, one can show the
difference in diagonal entries converges in probability to 0 zero. Then
looking at all the two-by-two principal minors, one can show the off diagonal
entries converge in probability to 0 as well. ∎

######  Remark 13.

Lemma C.6 is not true for asymmetric matrices. For example,
$\widehat{A}_{n}=\begin{bmatrix}0&1\\\ -1&0\end{bmatrix}$ and
$\widehat{A}_{n}=\begin{bmatrix}0&\hskip 5.69054pt0\\\ 0&\hskip
5.69054pt0\end{bmatrix}$. $t^{\prime}\widehat{A}_{n}t=t^{\prime}A_{n}t=0$ for
all $t\in\mathbf{R}^{2}$.

We now state a tensor inequality. Consider a fourth-order $n$-dimensional
tensor. We denote it as $\mathbf{A}=[a_{ijkl}]\in\mathbf{R}^{n\times n\times
n\times n}$. We shall understand it as a multi-linear function:

$\mathbf{A}:\mathbf{R}^{n}\times\mathbf{R}^{n}\times\mathbf{R}^{n}\times\mathbf{R}^{n}\to\mathbf{R}$

where for $x,y,z,a\in\mathbf{R}^{n}$,

$\mathbf{A}(w,x,y,z)=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}a_{ijkl}w_{i}x_{j}y_{k}z_{l}$
(65)

Consider the following maximization problem:

$\displaystyle\max_{w,x,y,z}\mathbf{A}(w,x,y,z)$ $\displaystyle\text{ subject
to
}\sum_{i=1}^{n}w_{i}^{4}=\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}^{4}=\sum_{i=1}^{n}z_{i}^{4}=1$

We denote its optimal value as $\sigma_{\max}(\mathbf{A})$. This optimal value
exists because we are optimizing a continuous function on a compact set. Note
further the problem

$\displaystyle\max_{w,x,y,z}|\mathbf{A}(w,x,y,z)|$ $\displaystyle\text{
subject to
}\sum_{i=1}^{n}w_{i}^{4}=\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}^{4}=\sum_{i=1}^{n}z_{i}^{4}=1$

has the same solution as the problem above because we can always take the
negative of one of the vectors.

Note for any vector $w,x,y,z\in\mathbf{R}^{n}$, we have:

$\displaystyle|\mathbf{A}(w,x,y,z)|=||w||_{4}||x||_{4}||y||_{4}||z||_{4}|\mathbf{A}(\frac{w}{||w||_{4}},\frac{x}{||x||_{4}},\frac{y}{||y||_{4}},\frac{z}{||z||_{4}})|\leq\sigma(\mathbf{A})||w||_{4}||x||_{4}||y||_{4}||z||_{4}$

The following lemma bounds $\sigma_{\max}(A)$. We define quantities:

$||\mathbf{A}||_{-i}=\max_{i}\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}|a_{ijkl}|$

and analogously,

$||\mathbf{A}||_{-j}=\max_{j}\sum_{i=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}|a_{ijkl}|$

and similarly for $||A||_{-k}$ and $||A||_{-l}$. We define:

$||\mathbf{A}||_{\infty}=\max\\{||\mathbf{A}||_{-i},||\mathbf{A}||_{-j},||\mathbf{A}||_{-k},||\mathbf{A}||_{-l}\\}$

######  Lemma C.7.

Consider a fourth-order $n-$dimensional tensor
$\mathbf{A}=[a_{ijkl}]\in\mathbf{R}^{n\times n\times n\times n}$.

$\sigma_{\max}(\mathbf{A})\leq||\mathbf{A}||_{\infty}$

###### Proof.

First note if $a_{ijkl}=0$ for all $ijkl$ indices, this inequality is
trivially satisfied. Thus we assume there exists at least one $a_{ijkl}\not=0$
for an $ijkl$ index.
We observe that the objective function (65) is continuous and the feasible set
is compact, so at least one optimal solution exists. Moreover, 0 is not in the
feasible set. As a result, we conclude that a solution for the optimization
problem exists and the local independence constraint qualification is
satisfied at each solution (Nocedal and Wright, (2006), P320). Let
$\\{w^{*},x^{*},y^{*},z^{*}\\}\subset\mathbf{R}^{n}$ denote one of the optimal
solutions. For the vector $w^{*}$, the KKT condition states that there exists
a $\lambda_{w}^{*}$ such that

$\frac{\partial}{\partial
w_{i}}\mathbf{A}(w^{*},x^{*},y^{*},z^{*})-\lambda_{w}^{*}4w_{i}^{*3}=0,\text{
for all }i=1,...,n$

We have

$\frac{\partial}{\partial
w_{i}}\mathbf{A}(w^{*},x^{*},y^{*},z^{*})=\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}a_{ijkl}x^{*}_{j}y^{*}_{k}z^{*}_{l}$

We further have

$\displaystyle
4\sum_{i=1}^{n}w_{i}^{*}(\lambda_{w}^{*}w_{i}^{*3})=\sum_{i=1}^{n}w_{i}^{*}\frac{\partial}{\partial
w_{i}}\mathbf{\mathbf{A}}(w^{*},x^{*},y^{*},z^{*})=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}a_{ijkl}w^{*}_{i}x^{*}_{j}y^{*}_{k}z^{*}_{l}=\mathbf{\mathbf{A}}(w^{*},x^{*},y^{*},z^{*}),$

Given the constraint $\sum_{i=1}^{n}(w_{i}^{*})^{4}=1$, we have the following
equality:

$\lambda_{w}^{*}=\frac{1}{4}\mathbf{A}(w^{*},x^{*},y^{*},z^{*})$

Now notice:

$\displaystyle
4\lambda^{*}_{w}||w^{*}||_{\infty}^{3}=4\lambda^{*}_{w}\max_{i}|w_{i}^{*}|^{3}=\max_{i}|\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}a_{ijkl}x^{*}_{j}y^{*}_{k}z^{*}_{l}|$
$\displaystyle\leq\max_{i}\\{\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}|a_{ijkl}|\\}\times||x^{*}||_{\infty}\times||y^{*}||_{\infty}\times||z^{*}||_{\infty}$
$\displaystyle\leq||\mathbf{A}||_{\infty}\times||x^{*}||_{\infty}\times||y^{*}||_{\infty}\times||z^{*}||_{\infty},$

Notice a similar argument also works for $x^{*}$ and $y^{*}$ and $z^{*}$, and
$\lambda^{*}_{w}=\lambda^{*}_{x}=\lambda^{*}_{y}=\lambda^{*}_{z}=\frac{1}{4}\mathbf{A}(w^{*},x^{*},y^{*},z^{*})$.
Define $\lambda^{*}=\frac{1}{4}\mathbf{A}(w^{*},x^{*},y^{*},z^{*})$. We have
the following four equations:

$\displaystyle 4\lambda^{*}||w^{*}||^{3}_{\infty}$
$\displaystyle\leq||\mathbf{A}||_{\infty}\times||x^{*}||_{\infty}\times||y^{*}||_{\infty}\times||z^{*}||_{\infty}$
$\displaystyle 4\lambda^{*}||x^{*}||^{3}_{\infty}$
$\displaystyle\leq||\mathbf{A}||_{\infty}\times||w^{*}||_{\infty}\times||y^{*}||_{\infty}\times||z^{*}||_{\infty}$
$\displaystyle 4\lambda^{*}||y^{*}||^{3}_{\infty}$
$\displaystyle\leq||\mathbf{A}||_{\infty}\times||w^{*}||_{\infty}\times||x^{*}||_{\infty}\times||z^{*}||_{\infty}$
$\displaystyle 4\lambda^{*}||z^{*}||^{3}_{\infty}$
$\displaystyle\leq||\mathbf{A}||_{\infty}\times||w^{*}||_{\infty}\times||x^{*}||_{\infty}\times||y^{*}||_{\infty}$

Define
$v^{*}=\max\\{||w^{*}||_{\infty},||x^{*}||_{\infty},||y^{*}||_{\infty},||z^{*}||_{\infty}\\}$,
we have the inequality:

$4\lambda^{*}\times(v^{*})^{3}\leq||\mathbf{A}||_{\infty}\times(v^{*})^{3}$

Since $v^{*}\not=0$, We then have the inequality
$\mathbf{A}(w^{*},x^{*},y^{*},z^{*})=4\lambda^{*}\leq||\mathbf{A}||_{\infty}$.
∎

Denote $||z||_{4}^{4}=\sum_{i,j}z^{4}_{ij}$, $||z||_{1}=\sum_{i,j}|z_{ij}|$
and $||z||_{\infty}=\max_{i,j}|z_{ij}|$ for a matrix $z$. For a fourth-order
tensor $\mathbf{A}=\left[a_{ijkl}\right]_{ijkl}$, define
$||\mathbf{A}||_{1}=\sum_{ijkl}|a_{ijkl}|$

######  Lemma C.8 (Infeasible variance estimation).

Let $\tilde{\mathbf{D}}$ be a valid variance bound matrix, and
$\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}$ be the inverse
probability weighted version of the bounding matrix $\tilde{\mathbf{D}}$. Let
$\mathbf{Q}$ denote the fourth-order tensor
$(\tilde{\mathbf{D}}\bigotimes\tilde{\mathbf{D}})\circ\mathbf{S}$, as defined
in Section 5.3. Let $z\in\mathbf{R}^{kn\times k}$.
Consider the estimator:

$\widehat{\tilde{\text{{Var}}}}=\frac{1}{n^{2}}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z\mathbf{1}^{\prime}\in\mathbf{R}^{k\times
k}$

for the quantity

$\tilde{\text{{Var}}}=\frac{1}{n^{2}}\mathbf{1}^{\prime}z\tilde{\mathbf{D}}z\mathbf{1}\in\mathbf{R}^{k\times
k}$

Then

$\widehat{\tilde{\text{{Var}}}}-\tilde{\text{{Var}}}=\sqrt{\frac{1}{n^{3}}\sigma_{\max}(\mathbf{Q})}\times\sqrt{\frac{1}{n}||z||_{4}^{4}}$

or

$\widehat{\tilde{\text{{Var}}}}-\tilde{\text{{Var}}}=\sqrt{\frac{1}{n^{4}}||\mathbf{Q}||_{1}}\times\sqrt{||z||_{\infty}^{4}}$

###### Proof.

For an arbitrary $t\in\mathbf{R}^{k}$, consider the quadratic form
$t^{\prime}(\widehat{\tilde{\text{{Var}}}}-\tilde{\text{{Var}}})t$. We are
interested in upper bounding its convergence rate.

Note $t^{\prime}\widehat{\tilde{\text{{Var}}}}t$ is unbiased for
$t^{\prime}\tilde{\text{{Var}}}t$. To upper bound the convergence rate, we
study its variance:

$\displaystyle\text{{Var}}\bigg{(}t^{\prime}\widehat{\tilde{\text{{Var}}}}t\bigg{)}$
$\displaystyle=\frac{1}{n^{4}}\text{{E}}\left[\left(t^{\prime}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}zt-t^{\prime}z\tilde{\mathbf{D}}zt\right)^{2}\right]$
(66)

Define a vector $a=t^{\prime}z\in\mathbf{R}^{kn}$. Use
$\mathds{R}_{i},i=1,...,kn$ to denote the $ith$ diagonal element of
$\mathds{R}$. Some algebra shows that (66) has the form:

$\frac{1}{n^{4}}\mathbf{Q}(t^{\prime}z,t^{\prime}z,t^{\prime}z,t^{\prime}z)=\frac{1}{n^{4}}\sum_{i,j,k,l=1}^{kn}\mathbf{COV}(\mathds{R}_{i}\mathds{R}_{j},\mathds{R}_{k}\mathds{R}_{l})\frac{\tilde{d}_{ij}\tilde{d}_{kl}}{\pi_{ij}\pi_{kl}}(t^{\prime}z)_{i}(t^{\prime}z)_{j}(t^{\prime}z)_{k}(t^{\prime}z)_{l},$

where $\tilde{d}_{ij}$ denotes the $ij$th entry of the matrix
$\tilde{\mathbf{D}}$, $\pi_{ij}$ denotes the $ij$th entry of the matrix
$\mathbf{p}$ and $(t^{\prime}z)_{i}$ denotes the $i$th entry of the vector
$t^{\prime}z$.

Note $\mathbf{Q}(\cdot,\cdot,\cdot,\cdot)$ is a fourth-order $kn$ dimensional
tensor. Then we have,

$\displaystyle\text{{Var}}\bigg{(}t^{\prime}\widehat{\tilde{\text{{Var}}}}t\bigg{)}=\frac{1}{n^{4}}\text{{E}}\left[\left(t^{\prime}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}zt-t^{\prime}z\tilde{\mathbf{D}}zt\right)^{2}\right]=\frac{1}{n^{4}}\mathbf{Q}(t^{\prime}z,t^{\prime}z,t^{\prime}z,t^{\prime}z)$
$\displaystyle\leq\frac{1}{n^{4}}\sigma_{\max}(\mathbf{Q})||t^{\prime}z||_{4}^{4}\leq\frac{1}{n^{4}}\sigma_{\max}(\mathbf{Q})k^{3}||z||_{4}^{4}\times||t||_{\infty}^{4}\leq\frac{1}{n^{4}}\sigma_{\max}(\mathbf{Q})k^{3}||z||_{4}^{4}\times||t||_{2}^{4}$

Thus we have

$t^{\prime}(\widehat{\tilde{\text{{Var}}}}-\tilde{\text{{Var}}})t=\sqrt{\frac{1}{n^{3}}\sigma_{\max}(\mathbf{Q})\frac{1}{n}||z||_{4}^{4}}\times
k^{3}||t||_{2}^{2}$

Alternatively, we can bound

$\displaystyle\text{{Var}}\bigg{(}t^{\prime}\widehat{\tilde{\text{{Var}}}}t\bigg{)}=\frac{1}{n^{4}}\text{{E}}\left[\left(t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z)\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\mathbf{diag}(z)\mathbf{1}t-t^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(z)\tilde{\mathbf{D}}\mathbf{diag}(z)\mathbf{1}t\right)^{2}\right]$
$\displaystyle\leq$
$\displaystyle||t||_{\infty}^{4}\times||z||_{\infty}^{4}\times\frac{1}{n^{4}}||\mathbf{Q}||_{1}\leq||t||_{2}^{4}\times||z||_{\infty}^{4}\times\frac{1}{n^{4}}||\mathbf{Q}||_{1}$

Thus we have

$t^{\prime}(\widehat{\tilde{\text{{Var}}}}-\tilde{\text{{Var}}})t=||t||_{2}^{2}\sqrt{\frac{1}{n^{4}}||\mathbf{Q}||_{1}\times||z||_{\infty}^{4}}$

∎

######  Lemma C.9.

(Feasible Estimators Converging to Infeasible Estimators) Let
$\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}$ be the inverse
probability weighted version of the bounding matrix $\tilde{\mathbf{D}}$.
Consider plug-in estimator of the form

$\frac{1}{n^{2}}\hat{z}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}\in\mathbf{R}^{k\times
k},$

where $\hat{z}\in\mathbf{R}^{kn}$. We have

$\displaystyle\frac{1}{n^{2}}\hat{z}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}-\frac{1}{n^{2}}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z$
$\displaystyle=$ $\displaystyle
O\left(\frac{1}{n^{2}}||\hat{z}-z||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}+\frac{1}{n^{2}}||\hat{z}-z||_{2}\times||z||_{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}\right)$

or

$\displaystyle\frac{1}{n^{2}}\hat{z}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}-\frac{1}{n^{2}}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z$
$\displaystyle=$ $\displaystyle
O\left(\frac{1}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}^{2}+\frac{1}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}\times||z||_{\infty}\right)$

###### Proof.

We have the following algebraic manipulation:

$\displaystyle\begin{split}&\frac{1}{n^{2}}\hat{z}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}-\frac{1}{n^{2}}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}z=\frac{1}{n^{2}}(\hat{z}-z)^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}+\frac{1}{n^{2}}z^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\hat{z}-\frac{1}{n^{2}}z^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z\\\
=&\frac{1}{n^{2}}(\hat{z}-z)^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)+\frac{1}{n^{2}}(\hat{z}-z)\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z+\frac{1}{n^{2}}z\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)\end{split}$
(67)

Using the $|||\cdot|||_{2}$ norm and the submultiplicativity of a matrix norm,
we have the quantity upper bounded by

$\displaystyle|||\frac{1}{n^{2}}(\hat{z}-z)^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)+\frac{1}{n^{2}}(\hat{z}-z)\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z+\frac{1}{n^{2}}z^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)|||_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n^{2}}|||(\hat{z}-z)\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)|||_{2}+\frac{1}{n^{2}}|||(\hat{z}-z)\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}z|||_{2}+\frac{1}{n^{2}}|||z^{\prime}\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}(\hat{z}-z)|||_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n^{2}}|||\hat{z}-z|||_{2}^{2}\times|||\mathds{R}|||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}+\frac{2}{n^{2}}|||\hat{z}-z|||_{2}\times|||\mathds{R}|||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}\times|||z|||_{2}$
$\displaystyle\leq\frac{1}{n^{2}}||\hat{z}-z||_{2}^{2}\times|||\mathds{R}|||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}+\frac{2}{n^{2}}||\hat{z}-z||_{2}\times||z||_{2}\times|||\mathds{R}|||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}$

Thus we have for (67)

$(\ref{(15)})=O\left(\frac{1}{n^{2}}||\hat{z}-z||_{2}^{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}+\frac{1}{n^{2}}||\hat{z}-z||_{2}\times||z||_{2}\times|||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}|||_{2}\right),$

where we used $|||\mathds{R}|||_{2}=1$ and changed $\frac{2}{n^{2}}$ to
$\frac{1}{n^{2}}$ for simplicity.
Similarly, using the $||\cdot||_{\infty}$ vector norm, we have an alternative
upper bound:

$\frac{1}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}^{2}+\frac{2}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}\times||z||_{\infty}$

Thus we have for (67)

$(\ref{(15)})=O\left(\frac{1}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}^{2}+\frac{1}{n^{2}}||\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}||_{1}\times||\hat{z}-z||_{\infty}\times||z||_{\infty}\right).$

∎

######  Lemma C.10.

Let $\mathbf{\Omega}\in\mathbf{R}^{n\times n}$ be a symmetric matrix,
$y\in\mathbf{R}^{n}$ and $x\in\mathbf{R}^{n}$ be two arbitrary column vectors.
The quantity
$x^{\prime}\mathbf{\Omega}^{\prime}\mathbf{diag}(y)\mathbf{diag}(y)\mathbf{\Omega}x$
can be upper bounded by :

$x^{\prime}\mathbf{\Omega}^{\prime}\mathbf{diag}(y)\mathbf{diag}(y)\mathbf{\Omega}x\leq|||\mathbf{\Omega}|||^{2}_{1}\sqrt{\sum_{i=1}^{n}x_{i}^{4}}\sqrt{\sum_{i=1}^{n}y_{i}^{4}}$

###### Proof.

Notice the stated quantity is the l2 norm of the vector
$x^{\prime}\mathbf{\Omega}^{\prime}\mathbf{diag}(y)\in\mathbf{R}^{n}$. Denote
the $j$th column of $\mathbf{\Omega}$ by $\Omega_{j}$, and the (i,j) th entry
of $\mathbf{\Omega}$ by $\mathbf{\Omega}_{ij}$. The $i$th entry of the vector
$x^{\prime}\mathbf{\Omega}^{\prime}\mathbf{diag}(y)$ is
$x^{\prime}\mathbf{\Omega}_{i}y_{i}$. The l2 norm thus can be written as

$\displaystyle\sum_{i=1}^{n}(x^{\prime}\mathbf{\Omega}_{i}y_{i})^{2}=\sum_{i=1}^{n}(\sum_{j=1}^{n}x_{j}\mathbf{\Omega}_{ji}y_{i})^{2}\leq\sum_{i=1}^{n}(\sum_{j=1}^{n}|x_{j}||\mathbf{\Omega}_{ji}||y_{i}|)^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}(\sum_{j=1}^{n}|x_{j}|\sqrt{|\mathbf{\Omega}_{ji}|}|y_{i}|\sqrt{|\mathbf{\Omega}_{ji}|})^{2}\leq\sum_{i=1}^{n}(\sum_{j=1}^{n}|x_{j}|^{2}|y_{i}|^{2}|\mathbf{\Omega}_{ji}|)\times(\sum_{j=1}^{n}|\mathbf{\Omega}_{ji}|)$
$\displaystyle\leq$
$\displaystyle|||\mathbf{\Omega}|||_{1}\sum_{i=1}^{n}\sum_{j=1}^{n}|x_{j}|^{2}|y_{i}|^{2}|\mathbf{\Omega}_{ji}|$

Now notice the expression
$\sum_{i=1}^{n}\sum_{j=1}^{n}|x_{j}|^{2}|y_{i}|^{2}|\mathbf{\Omega}_{ji}|$ is
the expanded expression for the quadratic form $(x\circ
x)^{\prime}|\mathbf{\Omega}|(y\circ y)$, where $|\mathbf{\Omega}|$ replaces
entries in $\mathbf{\Omega}$ with their absolute values and $\circ$ denotes
the hadamard product. Thus we can continue the inequality

$\displaystyle|||\mathbf{\Omega}|||_{1}\sum_{i=1}^{n}\sum_{j=1}^{n}|x_{j}|^{2}|y_{i}|^{2}|\mathbf{\Omega}_{ij}|\leq|||\mathbf{\Omega}|||_{1}|||(|\mathbf{\Omega}|)|||_{2}\times\sqrt{\sum_{i=1}^{n}x_{i}^{4}}\sqrt{\sum_{i=1}^{n}y_{i}^{4}}$
$\displaystyle\leq$
$\displaystyle|||\mathbf{\Omega}|||_{1}^{2}\sqrt{\sum_{i=1}^{n}x_{i}^{4}}\sqrt{\sum_{i=1}^{n}y_{i}^{4}}.$

For the last line we used the fact for a symmetric matrix,

$|||(|\mathbf{\Omega}|)|||_{2}=\lambda_{\max}(|\mathbf{\Omega}|)\leq|||(\mathbf{\Omega})|||_{1}=|||\mathbf{\Omega}|||_{1},$

where $\lambda_{\max}(|\mathbf{\Omega}|)$ denote the largest absolute
eigenvalue of $|\mathbf{\Omega}|$, ∎

######  Remark 14.

Taking $\mathbf{\Omega}=\begin{bmatrix}1&1&...&1\\\ 1&0&...&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 1&0&...&0\end{bmatrix}$,
$y=e_{1}\in\mathbf{R}^{n}$ and $x=\mathbf{1}\in\mathbf{R}^{n}$, one can show
the upper bound
$|||\mathbf{\Omega}|||_{1}|||(|\mathbf{\Omega}|)|||_{2}\times\sqrt{\sum_{i=1}^{n}x_{i}^{4}}\sqrt{\sum_{i=1}^{n}y_{i}^{4}}$
is tight in terms of the order of $n$.

Let $\Theta$ be a compact set in a finite dimensional Euclidean space, and
$B(\theta,\delta)$ denote a closed ball in $\Theta$ of radius $\delta\geq 0$
(with the $l_{2}$ norm) centered at $\theta$. The following lemma adapts
Theorem 1 Andrews, (1992) to our setting.

######  Lemma C.11.

Let $\Theta$ be a compact set in a finite dimensional Euclidean space and
$\mathbf{P}_{n}$ be the probability measure induced by the random assignments.
Consider a sequence of continuous deterministic functions
$Q_{n}(\cdot):\Theta\to\mathbf{R}$ and continuous stochastic functions
$\hat{Q}_{n}(\cdot):\Theta\to\mathbf{R}$.595959We assume the criterion
functions are continuous to avoid measurability issues. This is condition is
satisfied by all models considered in this paper.  If
$|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|=o_{p}(1)$ pointwise on $\theta\in\Theta$
and the stochastic function $\hat{Q}_{n}(\cdot)-Q_{n}(\theta)$ is uniformly
stochastically equicontinuous: for all $\epsilon>0$, there exists a $\delta>0$
such that

$\displaystyle\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)-\left(\hat{Q}_{n}(\theta^{\prime})-Q_{n}(\theta^{\prime})\right)|>\epsilon)<\epsilon$

uniformly for large $n$. Then
$\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|=o_{p}(1)$.

###### Proof.

Fix a given $\epsilon>0$ and let $\delta$ to be the corresponding radius in
the stochastic equicontinuity condition. Because $\Theta$ is compact, we can
find a finite cover of $\Theta$, $\\{B(\theta_{j},\delta)\\}_{j=1,...,J}$. We
then have:

$\displaystyle\mathbf{P}_{n}(\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|>2\epsilon)\leq$
$\displaystyle\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)-\left(\hat{Q}_{n}(\theta^{\prime})-Q_{n}(\theta^{\prime})\right)|>\epsilon)$
$\displaystyle+\mathbf{P}_{n}(\max_{j\leq
J}|\hat{Q}_{n}(\theta_{j})-Q_{n}(\theta_{j})|>\epsilon)<2\epsilon$

uniformly for large $n$. Then for each $\epsilon>0$, we can find a $n$ and a
$\tilde{\epsilon}<<\epsilon$ such that
$\mathbf{P}_{n}(\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|>\epsilon)\leq\mathbf{P}_{n}(\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|>\tilde{\epsilon})<\tilde{\epsilon}$.
Thus $\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|=o_{p}(1)$ ∎

The following lemma is the standard consistency proof for GMM estimators.

######  Lemma C.12.

Let $\Theta$ be a compact set in a finite dimensional Euclidean space and
$\mathbf{P}_{n}$ be the probability measure induced by the random assignments.
Consider a sequence of continuous deterministic functions
$Q_{n}(\cdot):\Theta\to\mathbf{R}$ and continuous stochastic functions
$\hat{Q}_{n}(\cdot):\Theta\to\mathbf{R}$. Let
$\theta_{n}\equiv\arg\inf_{\theta\in\Theta}Q_{n}(\theta)$ and define
$\hat{\theta}_{n}\equiv\arg\min_{\theta\in\Theta}\hat{Q}_{n}(\theta)$. If
there exists a positive $c$ such that $\inf_{\theta\in\Theta\backslash
B(\theta_{n},\epsilon)}Q_{n}(\theta)-Q_{n}(\theta)>c\epsilon^{2}$ uniformly
for all large $n$ and
$\sup_{\theta\in\Theta}|\hat{Q}_{n}(\theta)-Q_{n}(\theta)|=o_{p}(1)$, then

$\hat{\theta}_{n}-\theta_{n}=o_{p}(1).$

###### Proof.

We have

$\displaystyle\mathbf{P}_{n}(||\hat{\theta}_{n}-\theta_{n}||>2\epsilon)\leq$
$\displaystyle\mathbf{P}_{n}(Q_{n}(\hat{\theta}_{n})-Q_{n}(\theta_{n})>c4\epsilon^{2})$
$\displaystyle\leq$
$\displaystyle\mathbf{P}_{n}(Q_{n}(\hat{\theta}_{n})-Q_{n}(\theta_{n})-\hat{Q}_{n}(\hat{\theta}_{n})+\hat{Q}_{n}(\theta_{n})>c2\epsilon^{2})+\mathbf{P}_{n}(\hat{Q}_{n}(\hat{\theta}_{n})-\hat{Q}_{n}(\theta_{n})>c2\epsilon^{2})$
$\displaystyle\leq$ $\displaystyle
2\mathbf{P}_{n}(\sup_{\theta\in\Theta}|Q_{n}(\theta)-Q_{n}(\theta)|>c\epsilon^{2})\to
0,$

as $n\to\infty$, where
$\mathbf{P}_{n}(\hat{Q}_{n}(\hat{\theta}_{n})-\hat{Q}_{n}(\theta_{n})>c2\epsilon^{2})=0$
because $\hat{\theta}_{n}$ is the minimzer of $\hat{Q}_{n}(\theta)$. ∎

The following Lemma uses notation in Section 6.

######  Lemma C.13.

Define a GR estimator for arm a as
$\hat{\mu}_{n,a}=\frac{1}{n}\sum_{i}f^{a}(x_{i},\hat{\theta}_{n})+\frac{1}{n}\sum_{i}\frac{\mathds{R}_{ai}}{\pi_{ai}}\left(y_{ai}-f^{a}(x_{i},\hat{\theta}_{n})\right)$,
and
$\hat{\mu}_{n,a}^{L}=\frac{1}{n}\sum_{i}f^{a}(x_{i},\theta_{n})+\frac{1}{n}\sum_{i}\frac{\mathds{R}_{ai}}{\pi_{ai}}\left(y_{ai}-f^{a}(x_{i},\theta_{n})\right)$.
If $\hat{\theta}_{n}-\theta_{n}=O_{p}(\frac{1}{\sqrt{n}})$, and there exists a
positive integer $N$ such that the following conditions for $f^{a}$ are
satisfied uniformly for all $n\geq N$:

  1. (i)

$f^{a}(x_{i},\theta)$ is two times differentiable in $\theta$ for all $x_{i}$
values, i=1…,n.

  2. (ii)

There exists a $C_{1}$ and an $\epsilon>0$ such that
$\frac{1}{n}\sum_{a,i}||\nabla_{\theta}f^{a}(x_{i},\theta)||^{2}_{2}<C_{1}$
for all $\theta\in B(\theta_{n},\epsilon)$.

  3. (iii)

There exists a $C_{2}$ and an $\epsilon>0$ such that
$\frac{1}{n}\sum_{a,i}\sup_{\theta\in
B(\theta_{n},\epsilon)}||\nabla_{\theta\theta}f^{a}(x_{i},\theta)||_{1}<C_{2}$,

then
$\frac{1}{n}\sum_{i}(f^{a}(x_{i},\hat{\theta}_{n})-f^{a}(x_{i},\theta_{n}))^{2}=O_{p}(\frac{1}{n})$
and
$\hat{\mu}_{n,a}-\hat{\mu}_{n,a}^{L}=o_{p}(\frac{|||\mathbf{D}|||_{2}}{n})$

###### Proof.

We first show
$\frac{1}{n}\sum_{i}(f^{a}(x_{i},\hat{\theta}_{n})-f^{a}(x_{i},\theta_{n}))^{2}=O_{p}(\frac{1}{n})$.
We first observe two useful facts:

  1. 1.

For any $\theta_{1},\theta_{2},\tilde{\theta}\in B(\theta_{n},\epsilon)$ and
$||\theta_{1}-\theta_{2}||_{2}\leq\delta$, by (ii),

$\frac{1}{n}\sum_{i=1}^{n}(\nabla_{\theta}f^{a}(x_{i},\tilde{\theta})^{\prime}(\theta_{1}-\theta_{2}))^{2}\leq\\{\frac{1}{n}\sum_{i=1}^{n}||\nabla_{\theta}f^{a}(x_{i},\tilde{\theta})||_{2}^{2}\\}\times||\theta_{1}-\theta_{2}||_{2}^{2}\leq
C_{1}\delta^{2}$ (68)

  2. 2.

For any $\theta_{1},\theta_{2}\in B(\theta_{n},\delta)$ with
$\delta<\epsilon$, we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\theta_{2})-f(x_{i},\theta_{1}))^{2}=\frac{2}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\tilde{\theta}_{2})-f^{a}(x_{i},\theta_{1}))\nabla_{\theta}f^{a}(x_{i},\tilde{\theta}_{2})^{\prime}(\theta_{2}-\theta_{1})$
$\displaystyle\leq$ $\displaystyle
2\sqrt{\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\tilde{\theta}_{2})-f^{a}(x_{i},\theta_{1}))^{2}}\times\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\nabla_{\theta}f^{a}(x_{i},\tilde{\theta}_{2})^{\prime}(\theta_{1}-\theta_{2}))^{2}},$

where $\tilde{\theta}_{2}$ is between $\theta_{1}$ and $\theta_{2}$. Take
supremum on both side over $\theta_{2},\theta_{1},\tilde{\theta}_{2}$, we
arrive at:

$\displaystyle\sup_{\theta_{1},\theta_{2}\in
B(\theta_{n},\delta)}\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\theta_{2})-f^{a}(x_{i},\theta_{1}))^{2}$
$\displaystyle\leq$ $\displaystyle 2\sqrt{\sup_{\theta_{1},\theta_{2}\in
B(\theta_{n},\delta)}\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\theta_{2})-f^{a}(x_{i},\theta_{1}))^{2}}\times\sqrt{C_{1}\delta^{2}},$

which yields:

$\sup_{\theta_{1},\theta_{2}\in
B(\theta_{n},\delta)}\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\theta_{2})-f^{a}(x_{i},\theta_{1}))^{2}\leq
4C_{1}\delta^{2}$ (69)

Because $\hat{\theta}_{n}-\theta_{n}=O_{p}(\frac{1}{n})$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\hat{\theta}_{n})-f^{a}(x_{i},\theta_{n}))^{2}\leq$
$\displaystyle\sup_{\theta_{1},\theta_{2}\in
B(\theta_{n},||\hat{\theta}_{n}-\theta_{n}||)}\frac{1}{n}\sum_{i=1}^{n}(f^{a}(x_{i},\theta_{2})-f^{a}(x_{i},\theta_{1}))^{2}\leq
4C_{1}||\hat{\theta}_{n}-\theta_{n}||^{2}=O_{p}(\frac{1}{n}).$

This proves the first claim. Now the second claim follows by:

$\displaystyle\hat{\mu}_{n,a}-\hat{\mu}_{n,a}^{\scriptscriptstyle{\textnormal{L}}}=\frac{1}{n}\sum_{i=1}^{n}\left[f^{a}(x_{i},\hat{\theta}_{n})-f^{a}(x_{i},\theta_{n})\right]+\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}(f^{a}(x_{i},\theta_{n})-f^{a}(x_{i},\hat{\theta}_{n}))$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}f^{a}(x,\theta_{n})^{\prime}(\hat{\theta}_{n}-\theta_{n})+(\hat{\theta}_{n}-\theta_{n})^{\prime}\\{\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta\theta}f^{a}(x,\tilde{\theta}_{n})\\}(\hat{\theta}_{n}-\theta_{n})$
$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta}f^{a}(x,\theta_{n})^{\prime}(\hat{\theta}_{n}-\theta_{n})-(\hat{\theta}_{n}-\theta_{n})^{\prime}\\{\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta\theta}f^{a}(x,\tilde{\theta}_{n})\\}(\hat{\theta}_{n}-\theta_{n})$
$\displaystyle=o_{p}(\frac{1}{\sqrt{n}}),$

where $\tilde{\theta}_{n}$ is between $\hat{\theta}_{n}-\theta_{n}$. The final
line follows by noticing

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}f^{a}(x,\theta_{n})^{\prime}(\hat{\theta}_{n}-\theta_{n})-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta}f^{a}(x_{i},\theta_{n})^{\prime}(\hat{\theta}_{n}-\theta_{n})=-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}-\boldsymbol{\pi}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta}f^{a}(x_{i},\theta_{n})^{\prime}(\hat{\theta}_{n}-\theta_{n})$
$\displaystyle=$ $\displaystyle
O_{p}(\frac{1}{\sqrt{n}})O_{p}(\frac{1}{\sqrt{n}})$

by (ii) and Lemma C.3. Because $\tilde{\theta}_{n}\in B(\theta_{n},\epsilon)$
with probability approaching one, the middle term can be upper bounded by

$\displaystyle\begin{split}||\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}-\pi_{ai}}{\pi_{ai}}\nabla_{\theta\theta}f^{a}(x_{i},\tilde{\theta}_{n})||_{1}\leq\frac{1}{n}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\pi_{ai}}\sup_{\theta\in
B(\theta_{n},\epsilon)}||\nabla_{\theta\theta}f^{a}(x_{i},\theta)||_{1}+\frac{1}{n}\sum_{i=1}^{n}\sup_{\theta\in
B(\theta_{n},\epsilon)}||\nabla_{\theta\theta}f^{a}(x_{i},\theta)||_{1}.\end{split}$
(70)

The upper bound is of order $O_{p}(1)$ by (iii) and the Markov inequality.
Together with $\hat{\theta}_{n}-\theta_{n}=O_{p}(\frac{1}{\sqrt{n}})$, this
implies the $O_{p}(\frac{1}{n})$ rate for the second order remainder term. ∎

##  Appendix D Proofs in Section 3

###  D.1 Proof of Theorem 1

###### Proof.

We have by definition $\text{{E}}[\hat{m}_{n}^{s}]=m_{n}^{s}$ for
$s=1,...,l_{1}$. It can be seen, as in Lemma C.3,
$\text{{Var}}(\hat{m}_{n}^{s})=\frac{1}{n^{2}}\phi^{s}\mathbf{D}\phi^{s}\leq\frac{1}{n}||\phi^{s}||^{2}_{2}\times\frac{1}{n}|||\mathbf{D}|||_{2}$.
Thus we have
$\hat{m}_{n}^{s}-m_{n}^{s}=O_{p}(\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}})=o_{p}(1)$
for $s=1,...,l_{1}$. By Assumption 4-3, $||\hat{\mu}_{n}-\mu_{n}||_{2}\leq
C\sqrt{\sum_{s=1}^{l_{1}}(\hat{m}_{n}^{s}-m_{n}^{s})^{2}}$. Thus
$||\hat{\mu}_{n}-\mu_{n}||_{2}=O_{p}(\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}})$.
∎

###  D.2 Proof of Corollary 1

###### Proof.

We prove only for the GR estimators. Other estimators can be proved
analogously. To verify Assumption 4 for the GR estimators, notice:

$\displaystyle\hat{\mu}^{{\scriptscriptstyle{\textnormal{GR}}}}_{n}$
$\displaystyle=\underbrace{\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}y}_{(1)}-\underbrace{(\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi}^{-1})^{\prime}\mathbf{x})}_{(2)}\times\underbrace{(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x})^{-1}}_{(3)}\times\underbrace{(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y)}_{(4)}$
$\displaystyle=F(k\text{ moments in (1)},k+p\text{ moments in
(2)},\frac{(k+p)(k+p+1)}{2}$ $\displaystyle\text{ moments in (3)},k+p\text{
moments in (4)}\hskip 2.0pt)$

We rewrite the GR estimator using the moment estimator fomulation as
$F(\tilde{y}_{n},\tilde{x}_{n},\tilde{A}_{n},\tilde{b}_{n})$ where
$\tilde{y}_{n}\in\mathbf{R}^{k}$ denotes the k entries corresponding to the k
moments in (1), $\tilde{x}_{n}\in\mathbf{R}^{k}$ denotes the k+p entries
corresponding to the k+p moments in (2),
$\tilde{A}_{n}\in\mathbf{R}^{(k+p)\times(k+p)}$ is a symmetric matrix
corresponding to the $\frac{(k+p)(k+p+1)}{2}\text{ moments in (3)}$ and
$\tilde{b}_{n}\in\mathbf{R}^{(k+p)}$ denotes the k+p entries corresponding to
the k+p moments in (4).606060Some moments in (3) can appear in $\tilde{A}_{n}$
twice. Let $y_{n}=\frac{1}{n}\mathbf{1}^{\prime}y$, $x_{n}=0$,
$A_{n}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}$ and
$b=\frac{1}{n}\mathbf{x}\mathbf{m}\boldsymbol{\pi}y$. By Assumption 2,
Assumption 5 and Lemma C.1, $\lambda_{\min}(A_{n})$ is bounded away from 0
uniformly for $n\geq N$. For a sufficient small $\epsilon$ and by Corollary
6.3.8 in Horn and Johnson, (2012), $\lambda_{\min}(\tilde{A}_{n})$ is
uniformly bounded away from 0 for all
$||\tilde{A}_{n}-A_{n}||_{2}<\epsilon$.616161$||\cdot||_{2}$ is the Frobenious
norm of a matrix. For $\tilde{y}_{n},\tilde{x}_{n},\tilde{A}_{n}$ and
$\tilde{b}_{n}$ that are sufficiently close to ${y}_{n},0,A_{n}$ and $b_{n}$,
we have:

$\displaystyle||F(\tilde{y}_{n},\tilde{x}_{n},\tilde{A}_{n},\tilde{b}_{n})-F({y}_{n},0,A_{n},b_{n})||_{2}=||y_{n}-\tilde{y}_{n}-\tilde{x}_{n}\tilde{A}_{n}^{-1}\tilde{b}_{n}||_{2}$
$\displaystyle\leq||y_{n}-\tilde{y}_{n}||_{2}+||\tilde{x}_{n}\tilde{A}_{n}^{-1}\tilde{b}_{n}||\leq||y_{n}-\tilde{y}_{n}||_{2}+||\tilde{x}_{n}||_{2}\times||\tilde{A}_{n}^{-1}\tilde{b}_{n}||_{2}\leq||y_{n}-\tilde{y}_{n}||_{2}+C||\tilde{x}_{n}||_{2},$

where $C$ is a constant independent of $n$ and we use Assumption 2, Assumption
5, Lemma C.1 and Lemma C.2 to bound $||A_{n}^{-1}\tilde{b}_{n}||_{2}$ by a
constant C independent of n. Assumption 4-(iii) is satisfied by Assumption 1.
That $\frac{|||\mathbf{D}|||_{2}}{n}$ is of order o(1) is satisfied by
Assumption 6. ∎

###  D.3 Proof of Lemma 4.2

###### Proof.

The CI estimators converge to
$\frac{1}{n}\mathbf{1}\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}$.
The difference from average potential outcomes is:

$\frac{1}{n}\mathbf{1}^{\prime}y-\frac{1}{n}\mathbf{1}^{\prime}\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}$
(71)

Note

$\mathbf{1}^{\prime}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathbf{m}^{-1}\mathbf{m}\boldsymbol{\pi}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})$

By the definition of WLS,

$\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=0$

Thus equation (71) is 0 if each column of
$\mathbf{m}^{-1}\boldsymbol{\pi}^{-1}\mathbf{1}$ is in the column space of
$\mathbf{x}$, which is the condition stated in Lemma 4.2. ∎

###  D.4 Proof of Lemma 4.3

###### Proof.

We observe that
$\widehat{\mu}^{{\scriptscriptstyle{\textnormal{GR}}}}_{n}=\widehat{\mu}^{\scriptscriptstyle{\textnormal{CI}}}_{n}$
if

$\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}(y-\mathbf{x}\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=0$
(72)

By the definition of sample WLS,

$\mathbf{x}^{\prime}\mathbf{m}\mathds{R}(y-\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}})=0$

Thus equation (72) is 0 if each column
$\mathds{R}\boldsymbol{\pi}^{-1}\mathbf{1}$ is in the column space of
$\mathds{R}\mathbf{m}\mathbf{x}$, which is the condition stated in Lemma 4.3.
∎

###  D.5 Proof of Lemma 4.4

###### Proof.

The MI estimators converge to $\frac{1}{n}\mathbf{1}^{\prime}\pi
y+\frac{1}{n}\mathbf{1}^{\prime}(\mathbf{i}-\boldsymbol{\pi})\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}$.
The difference from average potential outcomes is:

$\displaystyle\begin{split}\frac{1}{n}\mathbf{1}^{\prime}y-\frac{1}{n}\mathbf{1}^{\prime}\pi
y-\frac{1}{n}\mathbf{1}^{\prime}(\mathbf{i}-\boldsymbol{\pi})\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=\frac{1}{n}\mathbf{1}^{\prime}(\mathbf{i}-\boldsymbol{\pi})(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\end{split}$
(73)

Note

$\displaystyle\mathbf{1}^{\prime}(\mathbf{i}-\boldsymbol{\pi})(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=\mathbf{1}^{\prime}(\mathbf{i}-\boldsymbol{\pi})\boldsymbol{\pi}^{-1}\mathbf{m}^{-1}\mathbf{m}\boldsymbol{\pi}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=\mathbf{1}^{\prime}(\boldsymbol{\pi}^{-1}-\mathbf{i})\mathbf{m}^{-1}\mathbf{m}\boldsymbol{\pi}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})$

By the definition of WLS,

$\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=0$

Thus equation (22) is 0 if each column of
$\mathbf{m}^{-1}(\boldsymbol{\pi}^{-1}-\mathbf{i})\mathbf{1}$ is in the column
space of $\mathbf{x}$, which is the condition stated in Lemma 4.4. ∎

###  D.6 Proof of Lemma 4.5

###### Proof.

We observe that
$\widehat{\mu}^{{\scriptscriptstyle{\textnormal{GR}}}}_{n}=\widehat{\mu}^{\scriptscriptstyle{\textnormal{MI}}}_{n}$
if

$\mathbf{1}^{\prime}(\boldsymbol{\pi}^{-1}-\mathbf{i})\mathds{R}(y-\mathbf{x}\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})=0$
(74)

By the definition of sample WLS,

$\mathbf{x}^{\prime}\mathbf{m}\mathds{R}(y-\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}})=0$

Thus equation (74) is 0 if each column
$\mathds{R}(\boldsymbol{\pi}^{-1}-\mathbf{i})\mathbf{1}$ is in the column
space of $\mathds{R}\mathbf{m}\mathbf{x}$, which is the condition stated in
Lemma 4.5. ∎

##  Appendix E Proofs in Section 4

###  E.1 Proof of Theorem 2

###### Proof.

By an argument similar to Lemma C.3, we have
$\hat{m}_{n}^{s}-m_{n}^{s}=O_{p}(\sqrt{\frac{|||\mathbf{D}|||_{2}}{n}})=o_{p}(1)$
for $s=1,...,l_{1}$. By Assumption 8, we have, with probability approaching
one,
$||\hat{\mu}_{n}-\mu_{n}^{\scriptscriptstyle{\textnormal{L}}}||_{2}=O_{p}(\frac{|||\mathbf{D}|||_{2}}{n})$.

For any linearized estimator it has the form:

$\mu_{n}+dF_{m_{n}}(\hat{m}_{n,r}-m_{n,r})\in\mathbf{R}^{k}$

Note the $s$th entry of the vector $\hat{m}_{n,r}-m_{n,r}$ has the form
$\frac{1}{n}1_{\scriptscriptstyle
kn}^{\prime}(\mathds{R}-\boldsymbol{\pi})\phi^{s}$. Thus the random term can
be written as:

$\displaystyle
dF_{m_{n}}(\hat{m}_{n,r}-m_{n,r})=\sum_{s=1}^{l_{1}}dF_{m_{n}}^{s}(\hat{m}^{s}_{n}-m^{s}_{n})=\sum_{s=1}^{l_{1}}dF_{m_{n}}^{s}\frac{1}{n}1_{\scriptscriptstyle
kn}^{\prime}\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi})\phi^{s}$
$\displaystyle=\sum_{s=1}^{l_{1}}dF_{m_{n}}^{s}\frac{1}{n}1_{\scriptscriptstyle
kn}^{\prime}\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi})\mathbf{diag}(\phi^{s})1_{\scriptscriptstyle
kn}=\frac{1}{n}\sum_{s=1}^{l_{1}}dF_{m_{n}}^{s}1_{\scriptscriptstyle
kn}^{\prime}\mathbf{diag}(\phi^{s})\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi})1_{\scriptscriptstyle
kn}$
$\displaystyle=\underbrace{\left(\frac{1}{n}\sum_{s=1}^{l_{1}}dF_{m_{n}}^{s}1_{\scriptscriptstyle
kn}^{\prime}\mathbf{diag}(\phi^{s})\right)}_{z^{\prime}}\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi})1_{\scriptscriptstyle
kn}$

Notice that $z$ involves only fixed quantities, thus we have:

$\text{{Var}}(dF_{m_{n}}(\hat{m}_{n,r}-m_{n,r}))=z^{\prime}\text{{Var}}(\boldsymbol{\pi}^{-1}(\mathds{R}-\boldsymbol{\pi})1_{\scriptscriptstyle
kn})z=z^{\prime}\mathbf{D}z$

∎

###  E.2 Proof of Corollary 2

###### Proof.

We prove only for the GR estimators. The Theorem can be proved in an analogous
way for other estimators. Recall the form of a GR estimator:

$\hat{\mu}_{n}^{{\scriptscriptstyle{\textnormal{GR}}}}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}y-\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}+\frac{1}{n}\mathbf{1}^{\prime}\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n},$

where
$\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}=(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x})^{+}(\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y)$
with a diagonal weighting matrix $\mathbf{m}$ with strictly positive entries.
The linear expansion of the GR estimator can be shown to have the form
$\frac{1}{n}\mathbf{1}^{\prime}y+\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}(\mathds{R}-\pi)(y-\mathbf{x}^{\prime}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})$.

Note the estimators depends on five sets of moments,
$\hat{m}_{y,n}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}y$,
$\hat{m}_{x,n}=\frac{1}{n}\mathbf{1}^{\prime}\boldsymbol{\pi}^{-1}\mathds{R}\mathbf{x}$,
$m_{x,n}=\frac{1}{n}\mathbf{1}^{\prime}\mathbf{x}$,
$\hat{m}_{xy,n}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}y$ and
$\hat{m}_{xx,n}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\mathds{R}\mathbf{x}$.
We also define $m_{y,n}=\frac{1}{n}\mathbf{1}^{\prime}y$,
$m_{xx,n}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}\mathbf{x}$
and $m_{xy,n}=\frac{1}{n}\mathbf{x}^{\prime}\mathbf{m}\boldsymbol{\pi}y$.
Consider moments $\tilde{m}_{xx,n}$ and $\tilde{m}_{xy,n}$ that are in a
sufficiently small local neighborhood of $m_{xx,n}$ and $m_{xy,n}$.
$\tilde{m}_{xx,n}$ is invertible because $m_{xx,n}$ is invertible by
Assumption 2, Assumption 5 and Lemma C.1, and the eigenvalues are continuous
functions of entries. Note we have the bound:

$\displaystyle\tilde{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}$
$\displaystyle=\tilde{m}_{xx,n}^{-1}\tilde{m}_{xy,n}-m_{xx,n}^{-1}m_{xy,n}$
$\displaystyle=(\tilde{m}_{xx,n}^{-1}-m_{xx,n}^{-1})\tilde{m}_{xy,n}+(m_{xx,n})^{-1}(\tilde{m}_{xy,n}-m_{xy,n})$
$\displaystyle=\tilde{m}_{xx,n}^{-1}\left(m_{xx,n}-\tilde{m}_{xx,n}\right)m_{xx,n}^{-1}\tilde{m}_{xy,n}+(m_{xx,n})^{-1}(\tilde{m}_{xy,n}-m_{xy,n})$

Then

$||\tilde{b}_{n}^{{\scriptscriptstyle{\textnormal{WLS}}}}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}\leq
C(||m_{xx,n}-\tilde{m}_{xx,n}||_{2}\times+||\tilde{m}_{xy,n}-m_{xy,n}||_{2}),$

where $C$ is a constant independent of $n$ for large $n$ by by Assumption 2,
Assumption 5 and Lemma C.1. The first $||\cdot||_{2}$ norm on the right-hand
side denotes the matrix Frobenius norm and the second $||\cdot||_{2}$ norm
denotes the standard vector norm. Thus for a local linear approximation we
have, for moments $\tilde{m}_{y}$, $\tilde{m}_{x}$, $\tilde{m}_{xx,n}$ and
$\tilde{m}_{xy,n}$ in a sufficiently small local neighborhood of $m_{y}$,
$m_{x}$, $m_{xx,n}$ and $m_{xy,n}$,

$\displaystyle F(\tilde{m}_{n})=$ $\displaystyle
m_{x,n}\tilde{m}_{xx,n}^{-1}\tilde{m}_{xy,n}+\tilde{m}_{y,n}-\tilde{m}_{x,n}\tilde{m}_{xx,n}^{-1}\tilde{m}_{xy,n}$
$\displaystyle F(m_{n})=$ $\displaystyle m_{y,n}$ $\displaystyle
dF_{m_{n}}(\tilde{m}_{n}^{r}-m_{n}^{r})=$
$\displaystyle(\tilde{m}_{y,n}-m_{y,n})-(\tilde{m}_{x,n}-m_{x,n})m_{xx,n}^{-1}m_{xy,n}$
$\displaystyle||F(\tilde{m}_{n})-F(m_{n})-dF_{m_{n}}(\tilde{m}_{n}^{r}-m_{n}^{r})||_{2}$
$\displaystyle\leq$
$\displaystyle||-(m_{x,n}-\tilde{m}_{x,n})\tilde{m}_{xx,n}^{-1}\tilde{m}_{xy,n}-(m_{x,n}-\tilde{m}_{x,n})m_{xx,n}^{-1}m_{xy,n}||_{2}\leq||(m_{x,n}-\tilde{m}_{x,n})(\tilde{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})||_{2}$
$\displaystyle\lesssim$
$\displaystyle||m_{x,n}-\tilde{m}_{x,n}||_{2}^{2}+||\tilde{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}^{2}\lesssim(||m_{x,n}-\tilde{m}_{x,n}||_{2}^{2}+||m_{xx,n}-\tilde{m}_{xx,n}||^{2}_{2}+||\tilde{m}_{xy,n}-m_{xy,n}||^{2}_{2}),$

where the constant is independent of $n$ by our assumption. This verifies
Assumption 8 for the GR estimators.
∎

###  E.3 Proof of Theorem 3 and Theorem 3’

###### Proof.

They are proved by using Lemma C.8 and Lemma C.9. ∎

###  E.4 Proof of Corollary 3 and Corollary 3’

We prove only for the GR estimators. The Corollary can be proved in an
analogous way for other estimators. The plug-in variance bound estimator for a
GR estimator is:

$\hat{\tilde{V}}(\hat{\mu}^{{\scriptscriptstyle{\textnormal{HT}}}})=\frac{1}{n}\mathbf{1}^{\prime}\mathbf{diag}(y-\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\mathds{R}\tilde{\mathbf{D}}_{{}{/}}{}_{\scriptscriptstyle\mathbf{p}}\mathds{R}\mathbf{diag}(y-\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\mathbf{1}\frac{1}{n}$

Using the notation in Theorem 3, we identify
$z=\mathbf{1}^{\prime}\mathbf{diag}(y-\mathbf{x}b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\in\mathbf{R}^{kn\times
k}$ and
$\hat{z}=\mathbf{1}^{\prime}\mathbf{diag}(y-\mathbf{x}\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\in\mathbf{R}^{kn\times
k}$. Note
$\mathbf{x}(\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})\in\mathbf{R}^{kn}$.
$||\cdot||_{2}$ denotes the Frobenius norm if applied to a matrix and the
$l_{2}$ vector norm if applied to a vector. We have

$\displaystyle\frac{1}{n}||\hat{z}-z||_{2}^{2}=\frac{1}{n}||\mathbf{1}^{\prime}\mathbf{diag}(\mathbf{x}(\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}))||_{2}^{2}=\frac{1}{n}||\mathbf{x}(\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n})||_{2}^{2}\leq\frac{1}{n}\lambda_{\max}(\mathbf{x}^{\prime}\mathbf{x})||\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}||\mathbf{x}||_{2}^{2}||\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}^{2}=O_{p}(\frac{|||\mathbf{D}|||_{2}}{n}),$

where $\frac{1}{n}||\mathbf{x}||_{2}^{2}=O(1)$ by Assumption 1 and
$||\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}^{2}=O_{p}(\frac{|||\mathbf{D}|||_{2}}{n})$
by Lemma C.5. $\frac{1}{n}||z||_{2}^{2}=O(1)$ by Assumption 1 and Lemma C.2.
Similarly under Assumption 1′ and by Lemma C.5,

$||\hat{z}-z||_{\infty}=||\mathbf{x}(\hat{b}-b)||_{\infty}\leq
C\times||\hat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_{n}||_{2}=O_{p}(\sqrt{\frac{\frac{1}{n}||\mathbf{D}||_{1}}{n}}),$

and $||z||_{\infty}=O(1)$ by Assumption 1′ and Lemma C.2.

##  Appendix F Proof in Section 6

###  F.1 Proof of Theorem 4

We first check the stochastic equicontinuity for large $n$.Notice by A10-(iv),
A6 and Lemma C.3 we have pointwise convergence for the criterion function: for
each $\theta\in\Theta$

$\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\pi_{ai}}\omega_{ai}g^{a}(y_{ai},x_{i},\theta)-\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\omega_{ai}g^{a}(y_{ai},x_{i},\theta)=o_{p}(1).$

By A10-(v), $Q_{n}(\theta)$ is continuous uniformly over $\theta\in\Theta$ and
for large $n$. Then,

$\displaystyle\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}\left|\hat{\mathcal{L}}_{n}(\theta)-\mathcal{L}_{n}(\theta)-(\hat{\mathcal{L}}_{n}(\theta^{\prime})-\mathcal{L}_{n}(\theta^{\prime})|\right]>3\epsilon)$
$\displaystyle\leq\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}\left|\hat{\mathcal{L}}_{n}(\theta)-\hat{\mathcal{L}}_{n}(\theta^{\prime})|\right]>2\epsilon)+\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}\left|\mathcal{L}_{n}(\theta)-\mathcal{L}_{n}(\theta^{\prime})|\right]>\epsilon),$

First note the second term is a degenerate probability event: it happens with
probability 1 or 0. For the first term, we take a $2\delta$-covering
$\mathcal{N}_{2\delta}$ of $\Theta$. Note $\mathcal{N}_{2\delta}$ is a finite
set by A10-(i). By the triangle inequality, for each $\theta$ there exists a
$\tilde{\theta}\in\mathcal{N}_{2\delta}$ such that
$d(\theta^{\prime},\tilde{\theta})<2\delta$ for all $\theta^{\prime}\in
B(\theta,\delta)$. Thus we can bound the first term:

$\begin{split}&\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}\left|\hat{\mathcal{L}}_{n}(\theta)-\hat{\mathcal{L}}_{n}(\theta^{\prime})|\right]>2\epsilon)\\\
=&\mathbf{P}_{n}(\sup_{\theta\in\Theta}\sup_{\theta^{\prime}\in
B(\theta,\delta)}\left|\hat{\mathcal{L}}_{n}(\theta)-\hat{\mathcal{L}}_{n}(\tilde{\theta})+\hat{\mathcal{L}}_{n}(\tilde{\theta})-\hat{\mathcal{L}}_{n}(\theta^{\prime})|\right]>2\epsilon)\\\
\leq&2\mathbf{P}_{n}(\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\pi_{ai}}D(y_{ai},x_{i})h(2\delta)>\epsilon)\leq
2\frac{\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}D(y_{ai},x_{i})h(2\delta)}{\epsilon}\end{split}$
(75)

where the last line is by the Markov inequality. By setting $\delta$ small
enough, we prove the desired inequality with A10-(v). With Lemma C.11 and
Lemma C.12, we conclude $\hat{\theta}_{n}-\theta_{n}=o_{p}(1)$.

By A10-(ii), $\hat{\theta}_{n}$ is in the interior of $\Theta$ with
probability approaching one. Hence by the first order condition,
$\hat{\theta}_{n}$ satisfies:

$\nabla_{\theta}\hat{\mathcal{L}}_{n}(\hat{\theta}_{n})=\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta}g^{a}(y_{ai},x_{i},\hat{\theta}_{n})=0$

A Taylor expansion around $\theta=\theta_{n}$ yields:

$\displaystyle\begin{split}0=&\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta}g^{a}(y_{ai},x_{i},\theta_{n})\\\
+&\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta\theta}g^{a}(y_{ai},x_{i},\theta_{n})(\hat{\theta}_{n}-\theta_{n})+o(||\hat{\theta}_{n}-\theta_{n}||)\end{split}$
(76)

The last term is justified by bounding the higher order reminder terms as
follows. Take the $t$th entry of the gradient
$\nabla_{\theta}\hat{\mathcal{L}}_{n}$. This row corresponds to the partial
derivative of $\hat{\mathcal{L}}_{n}$ with respect to the $t$th parameter
$\theta_{t}$. Its Taylor expansion has the form:

$\displaystyle 0=$
$\displaystyle\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\frac{\partial}{\partial\theta_{t}}g^{a}(y_{ai},x_{i},\theta_{n})+\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta}^{\prime}\frac{\partial}{\partial\theta_{t}}g^{a}(y_{ai},x_{i},\theta_{n})(\hat{\theta}_{n}-\theta_{n})$
$\displaystyle+$
$\displaystyle\frac{1}{n}(\hat{\theta}_{n}-\theta_{n})^{\prime}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta\theta}\frac{\partial}{\partial\theta_{t}}g^{a}(y_{ai},x_{i},\tilde{\theta}_{n})(\hat{\theta}_{n}-\theta_{n})$

where $\tilde{\theta}_{n}$ is between $\theta_{n}$ and $\hat{\theta}_{n}$. As
$\hat{\theta}_{n}$ enters $B(\theta_{n},\epsilon)$ with probability one, we
have for the remainder term,

$\displaystyle\begin{split}&|||\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta\theta}\frac{\partial}{\partial\theta_{t}}g^{a}(y_{ai},x_{i},\tilde{\theta}_{n})|||_{2}\leq\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}|||\nabla_{\theta\theta}\frac{\partial}{\partial\theta_{t}}g^{a}(y_{ai},x_{i},\tilde{\theta}_{n})|||_{2}\\\
\leq&\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\sup_{\theta\in
B(\theta_{n},\epsilon)}||\nabla_{\theta\theta\theta}g^{a}(y_{ai},x_{i},\tilde{\theta}_{n})||_{1}=O_{p}(1)\end{split}$
(77)

where for the $O_{p}(1)$ argument we used the Markov inequality, A10-(ix),
A10-(viii) and the equivalence of the $l_{1}$ and $l_{2}$ norm. Also by
A10-(ix), A10-(vi) and Lemma C.3,
$\frac{1}{n}\sum_{a,i}\omega_{ai}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta\theta}g^{a}(y_{ai},x_{i},\theta_{n})$
converges in probability to
$\frac{1}{n}\sum_{a,i}\omega_{ai}\nabla_{\theta\theta}g^{a}(y_{ai},x_{i},\theta_{n})$.
Hence it is invertible with probability approaching one. Thus using the above
argument and the fact $\hat{\theta}_{n}-\theta_{n}=o_{p}(1)$, we rearrange
(76) to get to

$\displaystyle\begin{split}&(\mathbf{i}_{s}+o_{p}(1))(\hat{\theta}_{n}-\theta_{n})=-(\frac{1}{n}\sum_{a,i}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta\theta}g^{a}(y_{ai},x_{i},\theta_{n}))^{-1}\frac{1}{n}\sum_{a=1}^{k}\sum_{i=1}^{n}\frac{\mathds{R}_{ai}}{\boldsymbol{\pi}_{ai}}\omega_{ai}\nabla_{\theta}g^{a}(y_{ai},x_{i},\theta_{n})\\\
=&O_{p}(\frac{1}{\sqrt{n}}),\end{split}$ (78)

where for the last line we used the FOC
$\frac{1}{n}\sum_{i=1}^{n}\omega_{ai}\nabla_{\theta}g^{a}(y_{ai},x_{i},\theta_{n})=0$,
A10-(vi), A10-(vii) and Lemma C.3.

The remaining conclusion follows from Assumption 6, Assumption 9, Assumption
11, Lemma C.13, Lemma C.8 and Lemma C.9.

###  F.2 Proof of Theorem 5

We first note
$\alpha_{n}^{c}=\frac{\frac{1}{n}c^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(y)\mathbf{D}\mathbf{diag}(f({\theta}_{n})\mathbf{1}c}{\frac{1}{n}c^{\prime}\mathbf{1}^{\prime}\mathbf{diag}(f(\theta_{n}))\mathbf{D}\mathbf{diag}(f(\theta_{n}))\mathbf{1}c}$
is uniformly bounded above for large $n$. The nominator,