# Doubly robust nearest neighbors in factor models

###### Abstract

In this technical note, we introduce an improved variant of nearest neighbors
for counterfactual inference in panel data settings where multiple units are
assigned multiple treatments over multiple time points, each sampled with
constant probabilities. We call this estimator a doubly robust nearest
neighbor estimator and provide a high probability non-asymptotic error bound
for the mean parameter corresponding to each unit at each time. Our guarantee
shows that the doubly robust estimator provides a (near-)quadratic improvement
in the error compared to nearest neighbor estimators analyzed in prior work
for these settings.

By Raaz Dwivedi1,2, Katherine Tian1, Sabina Tomkins3
Predrag Klasnja3, Susan Murphy1 and Devavrat Shah2

_Harvard University 1, Massachusetts Institute of Technology2 and University
of Michigan3_

##  1 Introduction

We consider a panel data setting with $N+1$ users that are assigned multiple
binary treatments taking value in $\\{0,1\\}$ over $T+1$ time points.111We use
$N+1$ and $T+1$ instead of $N$ and $T$ for theoretical convenience later on.
We use $A_{i,t}$ to denote the treatment assigned to unit $i\in[N]$ at time
$t\in[T]$ and $Y_{i,t}$ to denote the outcome observed for unit $i$ after the
treatment was assigned at time $t$. Overall, we assume the following
observation model:

$\textstyle
Y_{i,t}^{(a)}=\theta_{i,t}^{(a)}+\varepsilon_{i,t}\quad\text{for}\quad
a\in\mathcal{A}\quad\quad\text{and}\quad Y_{i,t}=Y_{i,t}^{(A_{i,t})},$ (1)

where $A_{i,t}$ are assigned by drawing Bernoulli random variables with
probability $p$ in an i.i.d manner. Here for each treatment $a\in\\{0,1\\}$,
the tuple $(Y_{i,t}^{(a)},\theta_{i,t}^{(a)})$ denotes the counterfactual
outcome (also referred to as the potential outcome) and its mean under
treatment $A_{i,t}=a$; $\varepsilon_{i,t}$ denotes exogenous and i.i.d zero
mean noise with variance $\sigma^{2}$ and $|\varepsilon_{i,t}|\leq
c_{\varepsilon}$ almost surely. Our model 1 follows the classical Neyman-Rubin
causal modeling set-up [7, 8] and to write it we make two assumptions: (i) The
_consistency assumption_ , also known as the stable unit treatment value
assumption (SUTVA), i.e., the observed outcome is equal to the potential
outcome and (ii) _No delayed or spill-over effect_ , i.e., the counterfactuals
at time $t$ do not depend on any past treatment. Our goal is to estimate
$\theta_{i,t}^{(a)}$ and provide an error guarantee for each estimate.
Throughout this work, we focus on the parameters for $a=1$ but the same
machinery can be used to derive results for the parameters for $a=0$. We drop
superscript $(a)$ (used to denote treatments) going forward and use the
simpler notation $\theta_{i,t}=\theta_{i,t}^{(1)}$ and
$\Theta=\\{\theta_{i,t}\\}_{i,t\in[N+1]\times[T+1]}$ where
$[K]\triangleq\\{1,\ldots,K\\}$.

It is well known that without any assumptions, estimating the mean parameters
is an impossible task due to more unknowns than the (noisy) knowns. To make
progress, we assume a latent factor model for the mean parameters that has
also been used in several prior works [3, 5, 1, 2]. In particular, we assume
that the mean counterfactual parameters $\Theta$ are a function of two sets of
$d$-dimensional latent factors—$\mathcal{U}\triangleq\\{u_{i},i\in[N+1]\\}$
and $\mathcal{V}\triangleq\\{v_{t},t\in[T+1]\\}$ and given these factors, they
are given by

$\displaystyle\theta_{i,t}=\\!\langle u_{i},v_{t}\rangle\quad\text{for}\quad
i\in[N\\!+\\!1],t\in[T\\!+\\!1],$ (2)

and assume that $\mathopen{}\mathclose{{}\left\|{u_{i}}}\right\|_{2}\leq
c_{u}$ and $\mathopen{}\mathclose{{}\left\|{v_{t}}}\right\|_{2}\leq c_{v}$
almost surely. Here, we assume a linear factorization for simplicity and later
on discuss the conditions under which our results can be extended to a non-
linear factor model of the form $\theta_{i,t}=f(u_{i},v_{t})$ for general non-
linear function $f$ satisfying some regularity conditions.

##  2 Algorithm

We start with an intuitive construction of the doubly robust nearest neighbors
(DR-NN) estimate for a given tuple $(i,t)$ in Sec. 2.1 followed by a formal
description in Sec. 2.2.

###  2.1 An intuitive construction of doubly robust nearest neighbors estimate

The notation, and definitions in this section should be treated as informal
versions of the estimates defined formally in Sec. 2.2. We abuse notation to
provide a parallel description behind a constructive intuition in building the
improved nearest neighbor estimate.

Here we consider an idealized setting with no noise so that
$\theta_{i,t}=Y_{i,t}$ and we assume that we observe all but one outcome
$Y_{i,t}$, and our goal is to estimate the unobserved $\theta_{i,t}$. There
are two natural variants of nearest neighbors to estimate $\theta_{i,t}$.

#### Unit-based estimate

In the first variant, that we call the unit-based estimate, one can compare
the counterfactuals of unit $i$ with all other units across time and identify
the set $\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ such that $u_{j}\approx
u_{i}$ for each $j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$. We compute
the unit-based estimate for $\theta_{i,t}$ using a simple average over the
counterfactuals of these neighboring units at time $t$:

$\displaystyle\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}=\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}Y_{j,t}}{S^{\textrm{unit}}_{i,t,\eta_{1}}}=\mathopen{}\mathclose{{}\left\langle\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}u_{j}}{S^{\textrm{unit}}_{i,t,\eta_{1}}},v_{t}}\right\rangle\triangleq\mathopen{}\mathclose{{}\left\langle\widehat{u}_{i},v_{t}}\right\rangle.$
(3)

The error of this estimate can be computed as

$\displaystyle\theta_{i,t}-\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}=\mathopen{}\mathclose{{}\left\langle
u_{i}-\widehat{u}_{i},v_{t}}\right\rangle\quad\Longrightarrow\quad\|{\theta_{i,t}-\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}}\|_{2}\leq\|{u_{i}-\widehat{u}_{i}}\|_{2}\|{v_{t}}\|_{2}.$
(4)

#### Time-based estimate

The second variant is the time-based estimate, where the roles of units and
time are swapped. In particular, we compare the counterfactuals at time $t$
with all other time points across units and identify the set
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ such that $v_{t^{\prime}}\approx
v_{t}$ for each $t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$. We
compute the time-based estimate for $\theta_{i,t}$ using a simple average over
the counterfactuals of unit $i$ at these neighboring time points:

$\displaystyle\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}=\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}Y_{i,t^{\prime}}}{S^{\textrm{time}}_{i,t,\eta_{2}}}=\mathopen{}\mathclose{{}\left\langle
u_{i},\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}v_{t^{\prime}}}{S^{\textrm{time}}_{i,t,\eta_{2}}}}\right\rangle\triangleq\mathopen{}\mathclose{{}\left\langle
u_{i},\widehat{v}_{t}}\right\rangle,$ (5)

for which the error can be computed as

$\displaystyle\theta_{i,t}-\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}=\mathopen{}\mathclose{{}\left\langle
u_{i},v_{t}-\widehat{v}_{t}}\right\rangle\quad\Longrightarrow\quad\|{\theta_{i,t}-\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}}\|_{2}\leq\|{u_{i}}\|_{2}\|{v_{t}-\widehat{v}_{t}}\|_{2}.$
(6)

#### A doubly-robust estimate

It is natural to ask if we can have best of both worlds, i.e., can we combine
the two estimates in a meaningful way to obtain an estimate
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ such that

$\displaystyle\theta_{i,t}-\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}=\mathopen{}\mathclose{{}\left\langle
u_{i}-\widehat{u}_{i},v_{t}-\widehat{v}_{t}}\right\rangle\quad\text{and}\quad\|{\theta_{i,t}-\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}}\|_{2}\leq\|{u_{i}-\widehat{u}_{i}}\|_{2}\|{v_{t}-\widehat{v}_{t}}\|_{2}?$
(7)

If we could, then the estimate
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ is doubly-robust in
the sense that as long as either of $\widehat{u}_{i}$ or $\widehat{v}_{t}$
provide a good estimate for $u_{i}$ or $v_{t}$ respectively, we would obtain a
good estimate for $\theta_{i,t}$. That is,
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ would be doubly
robust to estimating the two latent factors $(u_{i},v_{t})$. Moreover, the
error of the improved estimate is a product of the error in estimating
$u_{i}$, and $v_{t}$ from the displays 4 and 6, so that
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ would exhibit an
error better than both $\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}$ and
$\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}$. To find an expression for
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$, we expand the RHS of
7, and rearrange terms to find that

$\displaystyle\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$
$\displaystyle=\mathopen{}\mathclose{{}\left\langle
u_{i},\widehat{v}_{t}}\right\rangle+\mathopen{}\mathclose{{}\left\langle\widehat{u}_{i},v_{t}}\right\rangle-\mathopen{}\mathclose{{}\left\langle\widehat{u}_{i},\widehat{v}_{t}}\right\rangle$
(8)
$\displaystyle=\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}Y_{i,t^{\prime}}}{S^{\textrm{time}}_{i,t,\eta_{2}}}+\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}Y_{i,t^{\prime}}}{S^{\textrm{time}}_{i,t,\eta_{2}}}-\mathopen{}\mathclose{{}\left\langle\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}u_{j}}{S^{\textrm{unit}}_{i,t,\eta_{1}}},\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}v_{t^{\prime}}}{S^{\textrm{time}}_{i,t,\eta_{2}}}}\right\rangle$
(9)
$\displaystyle=\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}(Y_{i,t^{\prime}}+Y_{j,t}-Y_{j,t^{\prime}})}{S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}},$
(10)

which can be easily computed from the observed data. Of course, the
construction here was informal, and several key steps remain to be formalized,
which we now turn our attention to.

###  2.2 A formal description of the doubly robust nearest neighbors estimate

We begin with the two vanilla variants of nearest neighbors followed by the
doubly robust variant.

#### Unit-based estimate

Given a threshold $\eta_{1}$, define the unit neighbors
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ as

$\displaystyle\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}\triangleq\\!\mathopen{}\mathclose{{}\left\\{j\neq
i\\!\ \big{|}\
\\!\widehat{\rho}^{\textrm{unit}}_{i,t}(j)\\!\leq\\!\eta_{1}}\right\\}\,\
\text{with}\
\,\widehat{\rho}^{\textrm{unit}}_{i,t}(j)\\!\triangleq\\!\displaystyle\frac{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}(Y_{i,t^{\prime}}\\!-\\!Y_{j,t^{\prime}})^{2}}{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}}\ \text{for}\ j\neq i.$ (11)

In words, $\widehat{\rho}^{\textrm{unit}}_{i,t}(j)$ denotes the mean squared
distance of observed outcomes for unit $i$ and $j$ at time points other than
time $t$ when both units were in the control treatment; and
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ denotes the units that are within
$\eta_{1}$ distance from unit $i$. Intuitively,
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ denotes the units that are similar
to the unit of interest $i$. The _unit-based_ nearest neighbor estimate is
defined as a simple average of the outcomes at time $t$ for the units in
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ and in control at time $t$:

$\displaystyle\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}\triangleq\displaystyle\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}Y_{j,t}}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}}.$
(12)

if the denominator is not $0$; else the estimate is given by

$\displaystyle\quad\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}=Y_{i,t}\quad\text{if}\quad
A_{i,t}=1,\quad\text{else}\quad\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}=\frac{\sum_{j=1}^{N}A_{j,t}Y_{j,t}}{\sum_{j=1}^{N}A_{j,t}}.$
(13)

#### Time-based estimate

Given the threshold $\eta$, define the time neighbors
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ as

$\displaystyle\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}\triangleq\\!\mathopen{}\mathclose{{}\left\\{t^{\prime}\neq
t\\!\ \big{|}\
\\!\widehat{\rho}^{\textrm{time}}_{i,t}(t^{\prime})\\!\leq\\!\eta_{2}}\right\\}\,\
\text{with}\
\,\widehat{\rho}^{\textrm{time}}_{i,t}(t^{\prime})\\!\triangleq\\!\displaystyle\frac{\sum_{j\neq
i}A_{j,t}A_{j,t}(Y_{j,t}\\!-\\!Y_{j,t^{\prime}})^{2}}{\sum_{j\neq
i}A_{j,t}A_{j,t^{\prime}}}\ \text{for}\ t^{\prime}\neq t.$ (14)

In words, $\widehat{\rho}^{\textrm{time}}_{i,t}(t^{\prime})$ denotes the mean
squared distance of observed outcomes for time $t$ and $t^{\prime}$ across
units other than $i$ which were under the control at both the time points; and
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ denotes the time points that are
within $\eta_{2}$ distance from time $t$. Intuitively,
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ denotes the time points that are
similar to the time point of interest $t$. The _time-based_ nearest neighbor
estimate is defined as a simple average of the outcomes for unit $i$ at those
time points in $\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ at which it was in
control:

$\displaystyle\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}\triangleq\displaystyle\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}Y_{i,t^{\prime}}}{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}}.$
(15)

if the denominator is not $0$; else the estimate is given by the time analog
of 13:

$\displaystyle\quad\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}=Y_{i,t}\quad\text{if}\quad
A_{i,t}=1,\quad\text{else}\quad\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}=\frac{\sum_{t^{\prime}=1}^{T}A_{i,t^{\prime}}Y_{i,t^{\prime}}}{\sum_{t^{\prime}=1}^{T}A_{i,t^{\prime}}}.$
(16)

#### A doubly robust estimate

We now define the _doubly robust estimate_ studied in this work. With the
definitions 11 and 14 in place, and using the notation
$\boldsymbol{\eta}\triangleq(\eta_{1},\eta_{2})$, the doubly robust estimate
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ can be defined as

$\displaystyle\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}\triangleq\displaystyle\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}^{c}A_{j,t}^{c}A_{j,t^{\prime}}^{c}(Y_{i,t^{\prime}}\\!+\\!Y_{j,t}\\!-\\!Y_{j,t^{\prime}})}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}^{c}A_{j,t}^{c}A_{j,t^{\prime}}^{c}},$
(17)

when the denominator is not zero. Else, we can choose the estimates

$\displaystyle\quad\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}=Y_{i,t}\
\text{if}\ A_{i,t}=1,\ \ \ \text{else}\ \
\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}=\frac{\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}A_{i,t^{\prime}}^{c}A_{j,t}^{c}A_{j,t^{\prime}}^{c}(Y_{i,t^{\prime}}\\!+\\!Y_{j,t}\\!-\\!Y_{j,t^{\prime}})}{\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}A_{i,t^{\prime}}^{c}A_{j,t}^{c}A_{j,t^{\prime}}^{c}}$
(18)

if the denominator in the latter expression is non-zero, else either of the
estimates 13 and 16. In words, the doubly robust estimate makes use of data
from both similar units and similar time points. Each term in the numerator of
definition 17 is made of the sum of (i) the outcome of unit $i$ at a similar
time point $t^{\prime}$, and (ii) the outcome of a similar unit $j$ at time
$t$, (iii) with an offset term of the outcome of unit $j$ and time
$t^{\prime}$. Thus to be able to compute such a term, we need both the units
$i$ and $j$ to be in control at time $t^{\prime}$, and unit $j$ to be in
control at time $t$.

######  Remark 1 (Sample split for theoretical analysis).

To prove our results, we assume that the doubly robust estimate was
constructed using a data-split. In particular, we use do a $2\times 2$ data
split chosen uniformly at random: (i) Two equal splits of time points
$\mathopen{}\mathclose{{}\left\\{t^{\prime}\in[T]:t^{\prime}\neq t}\right\\}$
into $\mathbf{T}_{1}$ and $\mathbf{T}_{2}$, and (ii) two equal splits of units
$\mathopen{}\mathclose{{}\left\\{j\in[N]:j\neq i}\right\\}$ into
$\mathbf{N}_{1}$ and $\mathbf{N}_{2}$. We use time points $\mathbf{T}_{1}$ to
compute the unit distances $\widehat{\rho}^{\textrm{unit}}_{i,t}(j)$ for
$j\in\mathbf{N}_{1}$, and identify the unit-neighbors only in the set
$\mathbf{N}_{1}$, i.e.,
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}\subset\mathbf{N}_{1}$. Next, we use
units $\mathbf{N}_{2}$ to compute the time distances
$\widehat{\rho}^{\textrm{time}}_{i,t}(t^{\prime})$ for
$t^{\prime}\in\mathbf{T}_{2}$, and identify the time-neighbors only in the set
$\mathcal{T}_{2}$, i.e.,
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}\subset\mathbf{T}_{2}$. On the one
hand, this data split inflates the variance of our estimator. On the other
hand, it removes the (hard to analyze) bias that might arise due to data-reuse
in estimating the two neighbors sets
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}$ and
$\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$ without any data split.

###  2.3 Examples

Our results (established in the sequel) for the nearest neighbor algorithm
rely on the number of neighbors apriori available for a given unit and time.
Here to simplify notation, we drop the superscript unit and time, and we
understand the indices $(i,j)$ are reserved for units, and $(t,t^{\prime})$
are reserved for time points:

$\displaystyle\mathbf{S}_{i,r}$
$\displaystyle\\!\triangleq\\!\mathopen{}\mathclose{{}\left\\{j\neq i\
\big{|}\
(u_{j}\\!-\\!u_{i})^{\top}\\!\Sigma_{\textrm{time}}(u_{j}\\!-\\!u_{i})\\!\leq\\!r}\right\\},\quad\quad\qquad
S_{i,r}\triangleq|\mathbf{S}_{i,r}|,$ (19) $\displaystyle\mathbf{S}_{t,r}$
$\displaystyle\\!\triangleq\\!\mathopen{}\mathclose{{}\left\\{t^{\prime}\neq
t\ \big{|}\
(v_{t^{\prime}}\\!-\\!v_{t})^{\top}\\!\Sigma_{\textrm{unit}}(v_{t^{\prime}}\\!-\\!v_{t})\leq
r}\right\\},\quad\text{and}\quad S_{t,r}\triangleq|\mathbf{S}_{t,r}|$ (20)

where

$\displaystyle\Sigma_{\textrm{unit}}\triangleq\frac{1}{N}\sum_{i=1}^{N}u_{i}u_{i}^{\top}\quad\text{and}\quad\Sigma_{\textrm{time}}\triangleq\frac{1}{T}\sum_{t=1}^{T}v_{t}v_{t}^{\top}.$
(21)

We also define

$\displaystyle\lambda_{\textrm{unit}}=\lambda_{\min}(\Sigma_{\textrm{unit}})\quad\text{and}\quad\lambda_{\textrm{time}}=\lambda_{\min}(\Sigma_{\textrm{time}})$
(22)

We discuss two distinct examples for generating mechanisms on the latent
factors such that one can characterize the size of the sets defined in 19 and
20. In particular, we consider two settings for distributions on
$\mathcal{P}_{\textrm{unit}}$ and $\mathcal{P}_{\textrm{time}}$. We assume
that the distributions are such that the population level covariances satisfy
$\mathbb{E}[u_{i}u_{i}^{\top}]=\mathbf{I}_{d}$ and
$\mathbb{E}[v_{t}v_{t}^{\top}]=\mathbf{I}_{d}$:

######  Example (Discrete unit factors).

In this setting, we consider $\mathcal{P}_{\textrm{unit}}$ and
$\mathcal{P}_{\textrm{time}}$ as uniform distributions over
$M_{\textrm{unit}}$ and $M_{\textrm{time}}$ number of distinct $d$-dimensional
vectors. When
$u_{i}\stackrel{{\scriptstyle\textrm{i.i.d}}}{{\sim}}\mathcal{P}_{\textrm{unit}}$
and
$v_{t}\stackrel{{\scriptstyle\textrm{i.i.d}}}{{\sim}}\mathcal{P}_{\textrm{time}}$,
then we can show that $S_{i,r}\geq\frac{N}{2M_{\textrm{unit}}}$ for
$r=\Omega(T^{-\frac{1}{2}})$, and $S_{t,r}\geq\frac{T}{2M_{\textrm{time}}}$
for $r=\Omega(N^{-\frac{1}{2}})$. Moreover, standard matrix concentration
inequalities yield that $\lambda_{\textrm{unit}},\lambda_{\textrm{time}}\asymp
1$ for $N,T$ large enough.

######  Example (Continuous unit factors).

In this setting, we consider $\mathcal{P}_{\textrm{unit}}$ and
$\mathcal{P}_{\textrm{time}}$ as uniform distributions over
$M_{\textrm{unit}}$ and $M_{\textrm{time}}$ a bounded set in $\mathbb{R}^{d}$
(say $[-2,2]^{d}$). In this case, with
$u_{i}\stackrel{{\scriptstyle\textrm{i.i.d}}}{{\sim}}\mathcal{P}_{\textrm{unit}}$
and
$v_{t}\stackrel{{\scriptstyle\textrm{i.i.d}}}{{\sim}}\mathcal{P}_{\textrm{time}}$,
we can show that $S_{i,r}=\Theta(Nr^{\frac{d}{2}})$ for
$r=\Omega(T^{-\frac{1}{2}})$, and $S_{t,r}=\Theta(Tr^{\frac{d}{2}})$ for
$r=\Omega(N^{-\frac{1}{2}})$. Moreover, standard matrix concentration
inequalities yield that $\lambda_{\textrm{unit}},\lambda_{\textrm{time}}\asymp
1$ for $N,T$ large enough.

The two examples represent two different regimes for the latent factors:
Example Example can be interpreted as a parametric model with degree of
freedom $M$, and Example Example denotes a non-parametric model in $d$
dimensions.

##  3 Main results

In this section, we present our main result, a unit$\times$time wise non-
asymptotic guarantee for the doubly-robust nearest neighbor estimate and then
discuss its consequences.

Given $\eta$ and $\delta$, we define

$\displaystyle
e_{\delta}\\!\triangleq\\!94(c_{v}c_{u}\\!+\\!c_{\varepsilon})^{2}\sqrt{\log(36NTd/\delta)}$
(23)

Recall the definitions 19 and 20 of the sets $\mathbf{S}_{i,r}$ and
$\mathbf{S}_{t,r}$. To state the next result, we assume that the
hyperparameters $\eta_{1},\eta_{2}$ are such that

$\displaystyle
ce_{\delta}\mathopen{}\mathclose{{}\left((S_{i,\eta_{1}^{\prime}})^{-\frac{1}{2}}+(S_{t,\eta_{2}^{\prime}})^{-\frac{1}{2}}}\right)$
$\displaystyle\leq p\quad\text{where}\quad$ (24a)
$\displaystyle\eta_{1}^{\prime}\triangleq\eta_{1}\\!-\\!2\sigma^{2}\\!-\\!\frac{e_{\delta}}{p\sqrt{T}},$
$\displaystyle\quad\text{and}\quad\eta_{2}^{\prime}\triangleq\eta_{2}\\!-\\!2\sigma^{2}\\!-\\!\frac{e_{\delta}}{p\sqrt{N}}.$
(24b)

We are now ready to state our main result.

######  Theorem 3.1 (Error bound for
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$).

Given a fixed $\delta\in(0,1)$, suppose that the hyperparameter
$\boldsymbol{\eta}=(\eta_{1},\eta_{2})$ satisfies the regularity condition
24a. Then there exists universal constants $c,c^{\prime}$ such that
conditional on the latent factors $\mathcal{U},\mathcal{V}$, the doubly robust
estimate $\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ with sample
split (Rem. 1) satisfies

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}\\!-\\!\theta_{i,t})^{2}$
$\displaystyle\\!\leq\\!\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}+\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}},\quad\text{where}\quad$
(25a) $\displaystyle\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}$
$\displaystyle\\!\triangleq\\!\frac{2}{\lambda_{\textrm{unit}}\lambda_{\textrm{time}}}\\!\bigg{(}\\!\eta_{1}^{\prime}\\!+\\!\\!\frac{2e_{\delta}}{p\sqrt{T}}\\!\bigg{)}\\!\\!\bigg{(}\\!\eta_{2}^{\prime}\\!+\\!\\!\frac{2e_{\delta}}{p\sqrt{N}}\\!\bigg{)},\quad\text{and}\quad$
(25b) $\displaystyle\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}}$
$\displaystyle\\!\triangleq\\!\frac{c\log(NT/\delta)\sigma^{2}}{p}\\!\bigg{(}\\!\frac{1}{S_{i,\eta_{1}^{\prime}}}\\!\\!+\\!\frac{1}{S_{t,\eta_{2}^{\prime}}}\\!\\!+\\!\frac{1}{p^{2}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}\\!}\bigg{)},$
(25c) with probability at least $\displaystyle
1\\!-\\!2\delta\\!-\\!2S_{i,\eta_{1}^{\prime}}e^{-c^{\prime}p^{2}S_{i,\eta_{1}^{\prime}}}\\!-\\!2S_{t,\eta_{2}^{\prime}}e^{-c^{\prime}p^{2}S_{t,\eta_{2}^{\prime}}}\\!-\\!e^{-c^{\prime}p^{3}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}}.$
(25d)

See Sec. 4 for its proof where we also provide a tighter bound (the bound
above states worse logarithmic factors for brevity). Thm. 3.1 provides a
generic error bound on our doubly robust estimate for $\theta_{i,t}$
conditional on $\mathcal{U},\mathcal{V}$, with high probability over the
sampling mechanism and the noise variables, and decomposes the error into the
bias term $\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}$ 25b, and the variance
term $\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}}$ 25c.

###  3.1 Improvement for examples

Before turning to a detailed discussion of how this leads to an improved
guarantee over the vanilla NN variants, we first unpack the error via the next
corollary for the settings in Examples Example and Example. For clarity, we
assume $T=N$, assume $p=\Omega(N^{-\beta})$ for some $\beta\in[0,1]$—a larger
$\beta$ implies more missingness, treat $\delta$ as a constant, and track
dependency of the error only on $(N,T)$ in terms of $\widetilde{O}$ notation
that hides logarithmic dependencies, and constants that depend on
$(d,\sigma^{2},c_{u},c_{v},c_{\varepsilon})$.

######  Corollary 1 (Error improvement for examples).

Consider a setting satisfying the assumptions of Thm. 3.1 with $T=N$ and
$p=\Theta(N^{-\beta})$ and let $e_{p,\delta,N}=\frac{e_{\delta}}{p\sqrt{N}}$.
Then the following statements hold.

  1. (a)

For the setting in Example Example with
$M_{\textrm{unit}}=M_{\textrm{time}}=M$ and $N$ large enough such that
$e_{p,\delta,N}=o(\min_{u\neq
u^{\prime}}\mathopen{}\mathclose{{}\left\|{u\\!-\\!u^{\prime}}}\right\|^{2}_{2})$,
then for the choice $\eta_{1}=\eta_{2}=2\sigma^{2}+2e_{p,\delta,N}$, with
probability at least
$1-2\delta-2Ne^{-N^{1-2\beta}/M}-e^{-N^{2-3\beta}/M^{2}}$, we have

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}-\theta_{i,t})^{2}=\widetilde{\mathcal{O}}\mathopen{}\mathclose{{}\left(\frac{M}{N^{1-\beta}}\mathopen{}\mathclose{{}\left(1+\frac{M}{N^{1-2\beta}}}\right)}\right).$
(26)

  2. (b)

For the setting in Example Example when $\beta=0$ with the choice
$\eta_{1}=\eta_{2}=2\sigma^{2}+e_{p,\delta,N}+cN^{-\frac{2}{d+4}}$, with
probability at least $1-2\delta-cNe^{-cN^{4/(d+4)}}$, we have

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}-\theta_{i,t})^{2}=\widetilde{\mathcal{O}}(N^{-\frac{4}{d+4}}).$
(27)

See Sec. 5 for the proof, where we also provide the error expressions for part
(b) with a general $\beta\leq\frac{1}{2}$. For a constant $p$ ($\beta=0$)
Corollary 1 provides the following error for the DR-NN estimate:

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}-\theta_{i,t})^{2}=\begin{cases}\widetilde{\mathcal{O}}(M/N)\
\ \ \text{for \lx@cref{creftype~refnum}{example:finite}, and}\ \\\\[5.69054pt]
\widetilde{\mathcal{O}}(N^{-\frac{4}{d+4}})\ \text{for
\lx@cref{creftype~refnum}{example:continuous}}\ .\end{cases}$ (28)

Notably $M/N$ resembles a parametric rate, while $N^{-\frac{4}{d+4}}$ denotes
a non-parameteric rate that suffers from the curse of dimensionality. In both
settings, 28 establishes the dominance of our DR-NN estimate compared to the
unit or time-based NN estimate as the error scaling for
$\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}$ (which is same as that for
$\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}$ due to symmetry in
assumptions) as proven in prior work [4, Thm. 1], [3, Cor. 1]:222Both works do
not provide the stated bound explicitly: [4, Thm. 1] derives it implicitly in
their proofs, and [3, Thm. 3, Cor. 1] establishes the bound for an adaptively
sampled setting that provides the stated bound in our setting as a special
case (see [3, Rem. 5]).,

$\displaystyle(\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}-\theta_{i,t})^{2}=\begin{cases}\widetilde{\mathcal{O}}(N^{-\frac{1}{2}}+M/N)\
\ \text{for \lx@cref{creftype~refnum}{example:finite} and}\ \\\\[5.69054pt]
\widetilde{\mathcal{O}}(N^{-\frac{2}{d+2}})\quad\qquad\ \text{for
\lx@cref{creftype~refnum}{example:continuous}}\ .\end{cases}$ (29)

Putting 28 and 29 together, we conclude that compared to the unit- or time-
based NN estimate, DR-NN estimate provides a quadratic and near-quadratic in
$N$ improvement in the error for discrete factors and continuous factors
respectively.

We note that the assumption of uniformity in Examples Example and Example is
for simplifying the presentation of Corollary 1, and our proofs can be easily
adapted to establish error bounds that scale with the minimum value of the
probability mass function instead of $M$ with discrete factors, and with the
lower bound on density for continuous factors. With non-identity covariance
matrices $\Sigma_{\textrm{unit}},\Sigma_{\textrm{time}}$, the results involve
constants that scale with dimension $d$, the condition number, and the minimum
eigenvalue of the matrices—so that the prior conclusions continue to hold when
all these quantities are treated as constants as $N$ and $T$ grow. Finally,
when the factor probability mass functions or density functions are not lower
bounded, using arguments from from the [5, proof of Cor. 3] and in a manner
analogous to [3, Rem. 3], we can conclude that it remains lower bounded with
high probability over the draw of $u_{i}$ and $v_{t}$ so that suitably
modified Thm. 3.1 and 1 can be interpreted as a high probability error bound
on a uniformly drawn estimate from all the $N\times T$ estimates.

###  3.2 Discussion about error improvement in general settings

We now provide a detailed discussion on how Thm. 3.1 establishes two favorable
properties of the DR-NN compared to the unit- or time-based NN.

We omit $(i,t)$ in our notation and use the shorthand:

$\displaystyle\mathrm{Err}_{\eta_{\textrm{DR}}}^{\textrm{DR}}\triangleq\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}+\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}},$
(30)

where $\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}$ and
$\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}}$ respectively denote the bias
and variance of the DR-NN estimate from Thm. 3.1. Prior work [4, Thm. 1], [3,
Thm. 3] (also see Footnote 2) establishes that

$\displaystyle(\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}\\!-\\!\theta_{i,t})^{2}$
$\displaystyle\leq\underbrace{\frac{2}{\lambda_{\textrm{time}}}\\!\bigg{(}\\!\eta_{1}^{\prime}\\!+\\!\\!\frac{ce_{\delta}}{p\sqrt{T}}\\!\bigg{)}}_{\triangleq\mathbb{B}^{\textrm{unit}}_{\eta_{1}}}+\underbrace{\frac{c^{\prime}e_{\delta}\sigma^{2}}{pS_{i,\eta_{1}^{\prime}}}}_{\triangleq\mathbb{V}^{\textrm{unit}}_{\eta_{1}}}\triangleq\mathrm{Err}_{\eta_{1}}^{\textrm{unit}},\quad\text{and}\quad$
(31a)
$\displaystyle(\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}\\!-\\!\theta_{i,t})^{2}$
$\displaystyle\leq\underbrace{\frac{2}{\lambda_{\textrm{unit}}}\\!\bigg{(}\\!\eta_{2}^{\prime}\\!+\\!\\!\frac{ce_{\delta}}{p\sqrt{N}}\\!\bigg{)}}_{\triangleq\mathbb{B}^{\textrm{time}}_{\eta_{2}}}+\underbrace{\frac{c^{\prime}e_{\delta}\sigma^{2}}{pS_{t,\eta_{2}^{\prime}}}}_{\triangleq\mathbb{V}^{\textrm{time}}_{\eta_{2}}}\triangleq\mathrm{Err}_{\eta_{2}}^{\textrm{time}},$
(31b)

with high probability. (For simplicity we treat $p$ as a constant going
forward.) Comparing 31a and 31b with 25, we find that up to logarithmic
factors, the bias and variance of DR-NN estimator satisfies

$\displaystyle\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}\approx\mathbb{B}^{\textrm{unit}}_{\eta_{1}}\mathbb{B}^{\textrm{time}}_{\eta_{2}}\quad\text{and}\quad\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}}\approx\mathbb{V}^{\textrm{unit}}_{\eta_{1}}+\mathbb{V}^{\textrm{time}}_{\eta_{2}}+\textrm{higher-
order-term},$ (32)

i.e., roughly speaking the bias and variance of
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ are respectively
given by the product of biases, and sum of variances of the two vanilla NN
estimates ($\widehat{\theta}_{i,t,\eta_{1}}^{\textrm{unit}}$ and
$\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}$). Note that
$\mathbb{B}^{\textrm{unit}}_{\eta_{1}}$ and
$\mathbb{V}^{\textrm{unit}}_{\eta_{1}}$ exhibit a trade-off with $\eta_{1}$:
As $\eta_{1}$ decreases to $2\sigma^{2}$ (and consequently $\eta_{1}^{\prime}$
to $0$), the $\mathbb{B}^{\textrm{unit}}_{\eta_{1}}$ 31a decreases to
$\mathcal{O}(T^{-\frac{1}{2}})$, but $\mathbb{V}^{\textrm{unit}}_{\eta_{1}}$
increases as $S_{i,\eta_{1}^{\prime}}$ decreases to $0$ 19. The threshold
$\eta_{2}$ characterizes a similar trade-off between
$\mathbb{B}^{\textrm{time}}_{\eta_{2}}$ and
$\mathbb{V}^{\textrm{time}}_{\eta_{2}}$. Putting these observations together,
yield the following two main conclusions.

  1. 1.

Doubly robust to heterogeneity in user and time factors. First, we consider
the setting where one of the time and unit neigbors estimate provides trivial
error. When one of the factors, say time factors, are too heterogeneous so
that we are forced to choose $\eta_{2}$ large so that
$\eta_{2}-2\sigma^{2}=\Omega(1)$ and $S_{t,\eta_{2}^{\prime}}=\Omega(T)$ (note
that $S_{t,\eta_{2}^{\prime}}=T$ as $\eta_{2}\to\infty$). In such a case, we
have $\mathbb{B}^{\textrm{time}}_{\eta_{2}}=\mathcal{O}(1)$, and
$\mathbb{V}^{\textrm{time}}_{\eta_{2}}=\widetilde{\mathcal{O}}(T^{-1})$ and as
a result $\mathrm{Err}_{\eta_{2}}^{\textrm{time}}=\Omega(1)$, i.e., the time
neighbors estimate provides a trivial error and is not useful. However, for
the DR-NN estimate with $\boldsymbol{\eta}=(\eta_{1},\eta_{2})$, we have

$\displaystyle\mathbb{B}^{\textrm{DR}}_{\boldsymbol{\eta}}\precsim\mathbb{B}^{\textrm{unit}}_{\eta_{1}},\quad\text{and}\quad\mathbb{V}^{\textrm{DR}}_{\boldsymbol{\eta}}\precsim\mathbb{V}^{\textrm{unit}}_{\eta_{1}}\quad\Longrightarrow\quad\mathrm{Err}_{\boldsymbol{\eta}}^{\textrm{DR}}\precsim\mathrm{Err}_{\eta_{1}}^{\textrm{unit}}.$
(33)

Similarly, when $\eta_{1}$ has to be chosen too large, we can find that
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ enjoys roughly the
same error guarantee as $\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}$.
Putting the two pieces together, we can write

$\displaystyle\mathrm{Err}_{\boldsymbol{\eta}}^{\textrm{DR}}\precsim\min\mathopen{}\mathclose{{}\left\\{\mathrm{Err}_{\eta_{1}}^{\textrm{unit}},\,\mathrm{Err}_{\eta_{2}}^{\textrm{time}}}\right\\}.$
(34)

In other words, $\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$ is
doubly robust to the heterogeneity in the user and time factors—as long as one
of them provides enough good neighbors (either for unit $i$ or time $t$), our
DR estimator provides a non-trivial error guarantee.

  2. 2.

Improved error rates when unit and time neighbors provide vanishing errors. We
now consider the setting where both time and unit neigbors estimate provides
vanishing error. The bias-variance trade-off 31a implies that the optimal
choice of $\eta_{1}$ satisfies
$\mathbb{B}^{\textrm{unit}}_{\eta_{1}}\asymp\mathbb{V}^{\textrm{unit}}_{\eta_{1}}$;
denote this choice by $\eta_{\textrm{unit}}$. Similarly, denote the optimal
choice of $\eta_{2}$ for $\widehat{\theta}_{i,t,\eta_{2}}^{\textrm{time}}$ as
$\eta_{\textrm{time}}$. In the regime of optimal vanishing error, we have

$\displaystyle\mathrm{Err}_{\eta_{\textrm{unit}}}^{\textrm{unit}}\asymp\mathbb{B}^{\textrm{unit}}_{\eta_{\textrm{unit}}}\asymp\mathbb{V}^{\textrm{unit}}_{\eta_{\textrm{unit}}}\ll
1,\quad\mathrm{Err}_{\eta_{\textrm{time}}}^{\textrm{time}}\asymp\mathbb{B}^{\textrm{time}}_{\eta_{\textrm{time}}}\asymp\mathbb{V}^{\textrm{time}}_{\eta_{\textrm{time}}}\ll
1.$ (35)

These rough calculations immediately imply that the choice
$\boldsymbol{\eta}=(\eta_{\textrm{unit}},\eta_{\textrm{time}})$ is sub-optimal
for $\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$, since the error
for this estimate at such a choice is dominated by the variance term:

$\displaystyle\mathbb{B}^{\textrm{unit}}_{\eta_{\textrm{unit}}}\mathbb{B}^{\textrm{time}}_{\eta_{\textrm{time}}}\asymp\mathbb{V}^{\textrm{unit}}_{\eta_{\textrm{unit}}}\mathbb{V}^{\textrm{time}}_{\eta_{\textrm{time}}}\ll\mathbb{V}^{\textrm{unit}}_{\eta_{\textrm{unit}}}+\mathbb{V}^{\textrm{time}}_{\eta_{\textrm{time}}}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:dr_decomp}}}{{\Longrightarrow}}\mathrm{Err}_{(\eta_{\textrm{unit}},\eta_{\textrm{time}})}^{\textrm{DR}}=\mathcal{O}\mathopen{}\mathclose{{}\left(\mathbb{V}^{\textrm{unit}}_{\eta_{\textrm{unit}}}+\mathbb{V}^{\textrm{time}}_{\eta_{\textrm{time}}}}\right).$
(36)

In contrast, the optimal choices would satisfy
$\mathbb{B}^{\textrm{unit}}_{\eta_{1}}\mathbb{B}^{\textrm{time}}_{\eta_{2}}\asymp\mathbb{V}^{\textrm{unit}}_{\eta_{1}}+\mathbb{V}^{\textrm{time}}_{\eta_{2}}$.
Under weak regularity conditions, as we increase the value of
$(\eta_{1},\eta_{2})$ (in each coordinate) from
$\eta_{\textrm{unit}},\eta_{\textrm{time}}$ the bias terms in 36 increase, and
the variance terms decrease yielding that, at the optimum choice for DR-NN
estimate—denoted by
$\boldsymbol{\eta}_{\textrm{DR}}\triangleq(\eta_{1,\textrm{DR}},\eta_{2,\textrm{DR}})$,
we have
$\mathbb{V}^{\textrm{unit}}_{\eta_{1,\textrm{DR}}}\ll\mathbb{V}^{\textrm{unit}}_{\eta_{\textrm{unit}}}$,
and
$\mathbb{V}^{\textrm{time}}_{\eta_{2,\textrm{DR}}}\ll\mathbb{V}^{\textrm{time}}_{\eta_{\textrm{time}}}$,
and consequently

$\displaystyle\mathrm{Err}_{\boldsymbol{\eta}_{\textrm{DR}}}^{\textrm{DR}}\asymp\mathbb{B}^{\textrm{unit}}_{\eta_{1,\textrm{DR}}}\mathbb{B}^{\textrm{time}}_{\eta_{2,\textrm{DR}}}\asymp\mathbb{V}^{\textrm{unit}}_{\eta_{1,\textrm{DR}}}+\mathbb{V}^{\textrm{time}}_{\eta_{2,\textrm{DR}}}\ll\min\mathopen{}\mathclose{{}\left\\{\mathrm{Err}_{\eta_{\textrm{unit}}}^{\textrm{unit}},\mathrm{Err}_{\eta_{\textrm{time}}}^{\textrm{time}}}\right\\}.$
(37)

That is, when both hyperparameters are tuned, the error of DR-NN estimate
would provide a significant improvement over the error obtained by the either
of the unit or time neighbors estimate.

######  Remark 2 (No neighbors yields vacuous guarantee).

Thm. 3.1 applies only to the estimate 17 and is meaningful when there are
enough number of unit and time neighbors. It does not provide any guarantee
for the estimate 18. If any of the problem parameters $(N,T,p)$ or the
algorithm parameter $\eta$ is too small, the terms $S_{i,\eta_{1}^{\prime}}$
or $S_{t,\eta_{2}^{\prime}}$ might become $0$ so that the guarantee from Thm.
3.1 would be vacuous.

######  Remark 3 (Non-linear factor model).

Thm. 3.1 can be adapted for a Lipschitz non-linear factor model. In
particular, if instead of 2, the mean parameters $\Theta$ satisfy

$\displaystyle\theta_{i,t}\\!=f(u_{i},v_{t})$ (38)

for a twice continuously-differentiable function $f$ that is also
$L$-Lipschitz with $L$-Lipschitz gradients, then our doubly robust NN would
once again provide an improvement in the error over vanilla NN estimates. The
key intuition behind this improvement can be described via the following
taylor series decomposition. Let $\nabla_{u}f$ and $\nabla_{v}f$ denote the
partial derivatives of $f$ with respect to $u$ and $v$ respectively. Similarly
let $\nabla^{2}_{u}f,\nabla^{2}_{v}f,\nabla^{2}_{u,v}f$ denote the suitable
second-order partial derivatives. Then for an arbitrary set of vectors
$\\{u_{i},v_{t},u_{j},v_{t}\\}$ in the domain of $f$, standard algebra with
Taylor’s theorem yields that

$\displaystyle\mathopen{}\mathclose{{}\left(f(u_{j},v_{t})+f(u_{i},v_{t^{\prime}})-f(u_{j},v_{t^{\prime}})}\right)-f(u_{i},v_{t})$
(39)
$\displaystyle=(u_{j}-u_{i})^{\top}\mathopen{}\mathclose{{}\left(\nabla_{u}^{2}f(u^{\prime},v_{t})-\nabla_{u}^{2}f(u^{\prime\prime},v^{\prime\prime})}\right)(u_{j}-u_{i})$
(40)
$\displaystyle\qquad+(v_{t}^{\prime}-v_{t})^{\top}\mathopen{}\mathclose{{}\left(\nabla_{v}^{2}f(u_{i},v^{\prime})-\nabla_{v}^{2}f(u^{\prime\prime},v^{\prime\prime})}\right)(v_{t}^{\prime}-v_{t})$
(41)
$\displaystyle\qquad-2(u_{j}-u_{i})^{\top}\nabla_{u,v}^{2}f(u^{\prime\prime},v^{\prime\prime})(v_{t}^{\prime}-v_{t}),$
(42)

for some $u^{\prime},v^{\prime},u^{\prime\prime},v^{\prime\prime}$. Noting
that
$\sup_{u,v}\mathopen{}\mathclose{{}\left\|{\nabla_{u,v}^{2}f(u,v)}}\right\|_{\mathrm{op}}\leq
L$ due to the $L$-Lipschitzness of the gradients, we find that

$\displaystyle\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left(f(u_{j},v_{t})+f(u_{i},v_{t^{\prime}})-f(u_{j},v_{t^{\prime}})}\right)-f(u_{i},v_{t})}\right|\leq
2L(\|{u_{i}-u_{j}}\|_{2}+\|{v_{t}-v_{t^{\prime}}}\|_{2})^{2}$ (43)

This intuition serves as the non-linear analog of 7 for using the DR-NN
estimator as defined in 17. Rem. 4 highlights how the proof can be modified
for this setting to derive a suitable variant of Thm. 3.1.

##  4 Proof of Thm. 3.1: Error bound for
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}$

We make use of the following high probability bound for binomial random
variables throughout our proofs:

######  Lemma 4.1.

Given $X_{i}\sim\mathrm{Bernoulli}(p_{i})$ for $i\in[N]$, we have

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{N}X_{i}\leq\frac{\sum_{i=1}^{N}p_{i}}{2}}\right)\leq\exp\mathopen{}\mathclose{{}\left(-\frac{\sum_{i=1}^{N}p_{i}}{8}}\right).$
(44)

###### Proof.

The standard Binomial-Chernoff bound [6, Thm. 4.2] for independent bernoulli
variables $\mathopen{}\mathclose{{}\left\\{X_{i}}\right\\}$ with
$\mu\triangleq\sum_{i=1}^{N}\mathbb{E}[X_{i}]=\sum_{i=1}^{N}p_{i}$ is given by

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{N}X_{i}\leq(1-\delta)\mu}\right)\leq\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\cdot\delta^{2}\mu}\right)\quad\text{for}\quad\delta\in(0,1).$
(45)

Substituting $\delta=\frac{1}{2}$ yields the claim. ∎

###  4.1 Basic decomposition

Define the shorthands $M_{i,j,t,t^{\prime}}\triangleq
A_{i,t^{\prime}}A_{j,t}A_{j,t^{\prime}}$,

$\displaystyle\widehat{\theta}$
$\displaystyle\triangleq\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}(\theta_{i,t^{\prime}}\\!+\\!\theta_{j,t}\\!-\\!\theta_{j,t^{\prime}})}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}},\quad\text{and}\quad$
(46) $\displaystyle\widehat{\varepsilon}$
$\displaystyle\triangleq\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}(\varepsilon_{i,t^{\prime}}\\!+\\!\varepsilon_{j,t}\\!-\\!\varepsilon_{j,t^{\prime}})}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}}.$
(47)

In the proof, we would include the probability of an event that implies

$\displaystyle\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}>0,$
(48)

so we assume this condition going forward. Then by definition 17 and 1, we
have
$\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}=\widehat{\theta}+\widehat{\varepsilon}$,
so that

$\displaystyle(\theta_{i,t}-\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}})^{2}$
$\displaystyle=(\theta_{i,t}-\widehat{\theta})^{2}+\widehat{\varepsilon}^{2}+2(\theta_{i,t}-\widehat{\theta})\widehat{\varepsilon}\leq
2(\theta_{i,t}-\widehat{\theta})^{2}+2\widehat{\varepsilon}^{2}.$ (49)

We claim that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[(\theta_{i,t}\\!-\\!\widehat{\theta})^{2}\leq\frac{1}{\lambda_{\textrm{time}}\lambda_{\textrm{unit}}}\mathopen{}\mathclose{{}\left(\eta_{1}^{\prime}+\frac{2e_{\delta}}{p\sqrt{T}}}\right)\mathopen{}\mathclose{{}\left(\eta_{2}^{\prime}+\frac{2e_{\delta}}{p\sqrt{N}}}\right)}\right]\geq
1-\delta$ (50)

and

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\widehat{\varepsilon}^{2}\leq\frac{c\log(NT/\delta)\sigma^{2}}{p}\\!\bigg{(}\frac{1}{S_{i,\eta_{1}^{\prime}}}\\!+\\!\frac{1}{S_{t,\eta_{2}^{\prime}}}\\!+\\!\frac{1}{p^{2}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}}\bigg{)}}\right]$
(51) $\displaystyle\geq
1\\!-\\!\delta\\!-\\!2S_{i,\eta_{1}^{\prime}}e^{-c^{\prime}p^{2}S_{i,\eta_{1}^{\prime}}}\\!-\\!2S_{t,\eta_{2}^{\prime}}e^{-c^{\prime}p^{2}S_{t,\eta_{2}^{\prime}}}\\!-\\!e^{-c^{\prime}p^{3}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}},$
(52)

which in turn implies the claimed result in the theorem.

We now prove each of the bounds separately.

###  4.2 Proof of the bound 51 on variance

Define

$\displaystyle\widehat{\varepsilon}_{1}$
$\displaystyle\triangleq\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\varepsilon_{i,t^{\prime}}A_{i,t^{\prime}}\mathopen{}\mathclose{{}\left(\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}}\right)}{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\mathopen{}\mathclose{{}\left(\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}}\right)},$
(53) $\displaystyle\widehat{\varepsilon}_{2}$
$\displaystyle\triangleq\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\varepsilon_{j,t}A_{j,t}\mathopen{}\mathclose{{}\left(\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}A_{j,t^{\prime}}}\right)}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}\mathopen{}\mathclose{{}\left(\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}A_{j,t^{\prime}}}\right)},$
(54) $\displaystyle\quad\text{and}\quad\widehat{\varepsilon}_{3}$
$\displaystyle\triangleq\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{j,t}A_{i,t^{\prime}}A_{j,t^{\prime}}\varepsilon_{j,t^{\prime}}}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{j,t}A_{i,t^{\prime}}A_{j,t^{\prime}}}.$
(55)

Then by definition 47, we find that

$\displaystyle\widehat{\varepsilon}=\widehat{\varepsilon}_{1}+\widehat{\varepsilon}_{2}-\widehat{\varepsilon}_{3}\quad\Longrightarrow\quad\widehat{\varepsilon}^{2}\leq
3(\widehat{\varepsilon}_{1}^{2}+\widehat{\varepsilon}_{2}^{2}+\widehat{\varepsilon}_{3}^{2}).$
(56)

Let
$\mathfrak{A}\triangleq\\{A_{j,t^{\prime}},A_{j,t},A_{i,t^{\prime}},j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}\\}$.
Then note that $\mathfrak{A}$ is independent of the noise
$\mathfrak{Z}\triangleq\\{\varepsilon_{j,t^{\prime}},\varepsilon_{j,t},\varepsilon_{i,t^{\prime}},j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}\\}$
due to the data-split and the fact that the noise is exogenously generated.
Consequently, conditioned on $\mathfrak{A}$, using Berstein’s inequality (and
ignoring the second higher order term),333If $|X_{i}-\mu_{i}|<b$ a.s.,
$\mathrm{Var}(X_{i})\leq\sigma^{2}$, then
$\mathbb{P}(|\frac{1}{n}|\sum_{i=1}^{n}(X_{i}-\mu_{i})|\geq t)|\leq
2\exp(-nt^{2}/(\sigma^{2}+bt/3))$. we find that

$\displaystyle\widehat{\varepsilon}_{1}^{2}$ $\displaystyle\leq
c\sigma^{2}\log(6S^{\textrm{time}}_{i,t,\eta_{2}}/\delta)\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}\big{)}^{2}}{\big{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}\big{)}\big{)}^{2}},$
(57) $\displaystyle\widehat{\varepsilon}_{2}^{2}$ $\displaystyle\leq
c\sigma^{2}\log(6S^{\textrm{unit}}_{i,t,\eta_{1}}/\delta)\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}\big{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}A_{j,t^{\prime}}\big{)}^{2}}{\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}\big{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{m,t^{\prime}}A_{j,t^{\prime}}\big{)}\big{)}^{2}},$
(58) $\displaystyle\widehat{\varepsilon}_{3}^{2}$
$\displaystyle\leq\frac{c\sigma^{2}\log(6S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}/\delta)}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{j,t}A_{i,t^{\prime}}A_{j,t^{\prime}}},$
(59)

with probability $1-\delta$. Moreover note that the variables
$A_{j,t^{\prime}},A_{j,t},A_{i,t^{\prime}}$ in $\mathfrak{A}$ are independent
of the sets
$\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}$
due to the data-split, and hence Lem. 4.1 implies that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}}{p^{2}S^{\textrm{unit}}_{i,t,\eta_{1}}}\in\mathopen{}\mathclose{{}\left[\frac{1}{2},\frac{3}{2}}\right]|\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\right]$
$\displaystyle\geq
1-2\exp\mathopen{}\mathclose{{}\left(-\frac{p^{2}S^{\textrm{unit}}_{i,t,\eta_{1}}}{8}}\right)\
\text{for any $t^{\prime}$,}\ $ (60)
$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}}{pS^{\textrm{time}}_{i,t,\eta_{2}}}\geq\frac{1}{2}|\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\right]$
$\displaystyle\geq
1-\exp\mathopen{}\mathclose{{}\left(-\frac{pS^{\textrm{time}}_{i,t,\eta_{2}}}{8}}\right),$
(61)
$\displaystyle\quad\text{and}\quad\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{j,t}A_{i,t^{\prime}}A_{j,t^{\prime}}}{p^{3}S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}}\geq\frac{1}{2}|\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\right]$
$\displaystyle\geq
1-\exp\mathopen{}\mathclose{{}\left(-\frac{p^{3}S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}}{8}}\right).$
(62)

Using the probability bound 60, and a standard union bound argument, we find
that

$\displaystyle\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}\big{)}^{2}}{\bigg{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}A_{j,t^{\prime}}\big{)}\bigg{)}^{2}}$
$\displaystyle\leq\frac{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\big{(}\frac{3}{2}p^{2}S^{\textrm{unit}}_{i,t,\eta_{1}}\big{)}^{2}}{\bigg{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}\frac{1}{2}p^{2}S^{\textrm{unit}}_{i,t,\eta_{1}}\bigg{)}^{2}}$
(63)
$\displaystyle=\frac{9}{\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}}\leq\frac{18}{pS^{\textrm{time}}_{i,t,\eta_{2}}},$
(64)

with probability at least

$\displaystyle
1-\exp\mathopen{}\mathclose{{}\left(-\frac{pS^{\textrm{time}}_{i,t,\eta_{2}}}{8}}\right)-2S^{\textrm{unit}}_{i,t,\eta_{1}}\exp\mathopen{}\mathclose{{}\left(-\frac{p^{2}S^{\textrm{unit}}_{i,t,\eta_{1}}}{8}}\right).$
(65)

A similar computation yields that

$\displaystyle\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}\big{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{i,t^{\prime}}A_{j,t^{\prime}}\big{)}^{2}}{\big{(}\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}A_{j,t}\big{(}\sum_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}A_{m,t^{\prime}}A_{j,t^{\prime}}\big{)}\big{)}^{2}}\leq\frac{18}{pS^{\textrm{unit}}_{i,t,\eta_{1}}}$
(66)

with probability at least

$\displaystyle
1-\exp\mathopen{}\mathclose{{}\left(-\frac{pS^{\textrm{unit}}_{i,t,\eta_{1}}}{8}}\right)-2S^{\textrm{time}}_{i,t,\eta_{2}}\exp\mathopen{}\mathclose{{}\left(-\frac{p^{2}S^{\textrm{time}}_{i,t,\eta_{2}}}{8}}\right).$
(67)

Putting together the pieces, we find that conditional on
$\mathcal{U},\mathcal{V}$, we have

$\displaystyle\widehat{\varepsilon}_{1}^{2}\leq\frac{18c\sigma^{2}\log(\frac{6S^{\textrm{time}}_{i,t,\eta_{2}}}{\delta})}{pS^{\textrm{time}}_{i,t,\eta_{2}}},\quad\widehat{\varepsilon}_{2}^{2}\leq\frac{18c\sigma^{2}\log(\frac{6S^{\textrm{unit}}_{i,t,\eta_{1}}}{\delta})}{pS^{\textrm{unit}}_{i,t,\eta_{1}}},\quad\text{and}\quad\widehat{\varepsilon}_{3}^{2}\leq\frac{2c\sigma^{2}\log(6S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}/\delta)}{p^{3}S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}}$
(68)

with probability at least

$\displaystyle
1-\delta-e^{-{pS^{\textrm{unit}}_{i,t,\eta_{1}}}/{8}}-2S^{\textrm{time}}_{i,t,\eta_{2}}e^{-{p^{2}S^{\textrm{time}}_{i,t,\eta_{2}}}/{8}}-e^{-\frac{pS^{\textrm{unit}}_{i,t,\eta_{1}}}{8}}-2S^{\textrm{time}}_{i,t,\eta_{2}}e^{-{p^{2}S^{\textrm{time}}_{i,t,\eta_{2}}}/{8}}-e^{-p^{3}S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}/8},$
(69)

which can be further simplified (by lower bounding) as

$\displaystyle
1\\!-\\!\delta\\!-\\!(2S^{\textrm{time}}_{i,t,\eta_{2}}\\!+\\!1)e^{-{p^{2}S^{\textrm{time}}_{i,t,\eta_{2}}}/{8}}\\!-\\!(2S^{\textrm{time}}_{i,t,\eta_{2}}\\!+\\!1)e^{-{p^{2}S^{\textrm{time}}_{i,t,\eta_{2}}}/{8}}\\!-\\!e^{-p^{3}S^{\textrm{unit}}_{i,t,\eta_{1}}S^{\textrm{time}}_{i,t,\eta_{2}}/8},$
(70)

Next, we note that under event $\mathcal{E}_{\textrm{dist}}$, we have

$\displaystyle\mathbf{S}_{i,\eta_{1}^{\prime}}\subseteq\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}\quad\text{and}\quad\mathbf{S}_{t,\eta_{2}^{\prime}}\subseteq\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}},$
(71)

and since the data split was done using Bernoulli$(\frac{1}{2})$ flips, we
obtain that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[S^{\textrm{unit}}_{i,t,\eta_{1}}\geq\frac{S_{i,\eta_{1}^{\prime}}}{4}}\right]\geq
1-\exp\mathopen{}\mathclose{{}\left(-\frac{S_{t,\eta_{2}^{\prime}}}{16}}\right)\quad\text{and}\quad\mathbb{P}\mathopen{}\mathclose{{}\left[S^{\textrm{time}}_{i,t,\eta_{2}}\geq\frac{S_{t,\eta_{2}^{\prime}}}{4}}\right]\geq
1-\exp\mathopen{}\mathclose{{}\left(-\frac{S_{t,\eta_{2}^{\prime}}}{16}}\right).$
(72)

Putting the pieces together, we conclude that

$\displaystyle\widehat{\varepsilon}^{2}\leq\frac{c\log(NT/\delta)\sigma^{2}}{p}\\!\bigg{(}\frac{1}{S_{i,\eta_{1}^{\prime}}}\\!+\\!\frac{1}{S_{t,\eta_{2}^{\prime}}}\\!+\\!\frac{1}{p^{2}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}}\bigg{)}$
(73)

with probability at least

$\displaystyle
1\\!-\\!\delta\\!-\\!2S_{i,\eta_{1}^{\prime}}e^{-c^{\prime}p^{2}S_{i,\eta_{1}^{\prime}}}\\!-\\!2S_{t,\eta_{2}^{\prime}}e^{-c^{\prime}p^{2}S_{t,\eta_{2}^{\prime}}}\\!-\\!e^{-c^{\prime}p^{3}S_{i,\eta_{1}^{\prime}}S_{t,\eta_{2}^{\prime}}},$
(74)

for some universal $c,c^{\prime}>0$.

###  4.3 Proof of the bound 50 on bias

We have

$\displaystyle\theta_{i,t}-\widehat{\theta}$
$\displaystyle=\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}(\theta_{i,t}-(\theta_{i,t^{\prime}}\\!+\\!\theta_{j,t}\\!-\\!\theta_{j,t^{\prime}}))}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}}$
(75)
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:bilinear}}}{{=}}\frac{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}\mathopen{}\mathclose{{}\left\langle
u_{i}-u_{j},v_{t}-v_{t^{\prime}}}\right\rangle}{\sum_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}},t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}M_{i,j,t,t^{\prime}}}$
(76)

so that

$\displaystyle|\theta_{i,t}\\!-\\!\widehat{\theta}|\\!\leq\\!\max_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\max_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}|\mathopen{}\mathclose{{}\left\langle
u_{i}\\!-\\!u_{j},v_{t}\\!-\\!v_{t^{\prime}}}\right\rangle|\leq\max_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\mathopen{}\mathclose{{}\left\|{u_{i}\\!-\\!u_{j}}}\right\|_{2}\max_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\mathopen{}\mathclose{{}\left\|{v_{t}\\!-\\!v_{t^{\prime}}}}\right\|_{2}.$
(77)

Recall the notations 21 and 23. Then the next lemma shows that the distances
between units, and time concentrates around suitably defined “population-
level” quantities:

######  Lemma 4.2.

Under the setting of Thm. 3.1, for any fixed $\delta\in(0,1]$, the event
$\mathcal{E}_{\textrm{dist}}$

$\displaystyle\sup_{j\neq
i}|\widehat{\rho}^{\textrm{unit}}_{i,t}(j)\\!-\\!\rho^{\textrm{unit}}_{i,t}(j)|\stackrel{{\scriptstyle(a)}}{{\leq}}\frac{e_{\delta}}{p\sqrt{T}}\quad\text{and}\quad\sup_{t^{\prime}\neq
t}|\widehat{\rho}^{\textrm{time}}_{i,t}(t^{\prime})-\rho^{\textrm{time}}_{i,t}(t^{\prime})|\stackrel{{\scriptstyle(b)}}{{\leq}}\frac{e_{\delta}}{p\sqrt{N}},$
(78) $\displaystyle\ \text{where}\
\rho^{\textrm{unit}}_{i,t}(j)\triangleq\mathopen{}\mathclose{{}\left\|{u_{i}\\!-\\!u_{j}}}\right\|^{2}_{\Sigma_{\textrm{time}}}\\!+\\!2\sigma^{2}\quad\text{and}\quad\rho^{\textrm{time}}_{i,t}(t^{\prime})\triangleq\mathopen{}\mathclose{{}\left\|{v_{t}\\!-\\!v_{t^{\prime}}}}\right\|^{2}_{\Sigma_{\textrm{unit}}}\\!+\\!2\sigma^{2},$
(79)

satisfies $\mathbb{P}(\mathcal{E}_{\textrm{dist}}|\mathcal{U},\mathcal{V})\geq
1-\delta$, where $e_{\delta}$ is defined in 23.

We now apply Lem. 4.2 to finish our proof. Note that the definitions 21 and 22
imply that

$\displaystyle\mathopen{}\mathclose{{}\left\|{u_{i}-u_{j}}}\right\|_{2}^{2}\leq\frac{\rho^{\textrm{unit}}_{i,t}(j)-2\sigma^{2}}{\lambda_{\min}(\Sigma_{\textrm{time}})},$
(80)

and hence on the event $\mathcal{E}_{\textrm{dist}}$, we obtain that

$\displaystyle\max_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\\!\\!\mathopen{}\mathclose{{}\left\|{u_{i}-u_{j}}}\right\|_{2}^{2}\leq\max_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\frac{\rho^{\textrm{unit}}_{i,t}(j)-2\sigma^{2}}{\lambda_{\min}(\Sigma_{\textrm{time}})}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:event_dist}}}{{\leq}}\frac{\widehat{\rho}^{\textrm{unit}}_{i,t}(j)-2\sigma^{2}+\frac{e_{\delta}}{p\sqrt{T}}}{\lambda_{\textrm{time}}}$
(81)
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:reliable_nbr}}}{{\leq}}\frac{\eta_{1}-2\sigma^{2}+\frac{e_{\delta}}{p\sqrt{T}}}{\lambda_{\textrm{time}}}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:thresh_prime}}}{{=}}\frac{1}{\lambda_{\textrm{time}}}\mathopen{}\mathclose{{}\left(\eta_{1}^{\prime}+\frac{2e_{\delta}}{p\sqrt{T}}}\right).$
(82)

Similarly, we also find that

$\displaystyle\max_{j\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\\!\\!\mathopen{}\mathclose{{}\left\|{v_{t}-v_{t^{\prime}}}}\right\|_{2}^{2}\leq\frac{1}{\lambda_{\textrm{unit}}}\mathopen{}\mathclose{{}\left(\eta_{2}^{\prime}+\frac{2e_{\delta}}{p\sqrt{N}}}\right),$
(83)

on the event $\mathcal{E}_{\textrm{dist}}$. Putting these bounds together with
77, we conclude that

$\displaystyle(\theta_{i,t}\\!-\\!\widehat{\theta})^{2}\leq\frac{1}{\lambda_{\textrm{time}}\lambda_{\textrm{unit}}}\mathopen{}\mathclose{{}\left(\eta_{1}^{\prime}+\frac{2e_{\delta}}{p\sqrt{T}}}\right)\mathopen{}\mathclose{{}\left(\eta_{2}^{\prime}+\frac{2e_{\delta}}{p\sqrt{N}}}\right)$
(84)

with probability $1-\delta$ and we are done.

######  Remark 4.

When $\theta_{i,t}$ is given by a non-linear function as in Rem. 3, then the
calculations in 43 can be used to show that 77 would in fact be replaced by

$\displaystyle|\theta_{i,t}\\!-\\!\widehat{\theta}|\\!\leq\\!2L(\max_{j\in\mathbf{S}^{\textrm{unit}}_{i,t,\eta_{1}}}\mathopen{}\mathclose{{}\left\|{u_{i}\\!-\\!u_{j}}}\right\|_{2}+\max_{t^{\prime}\in\mathbf{S}^{\textrm{time}}_{i,t,\eta_{2}}}\mathopen{}\mathclose{{}\left\|{v_{t}\\!-\\!v_{t^{\prime}}}}\right\|_{2})^{2}.$
(85)

Then a suitable adaptation of the discussion in [3, Sec. 4.1] can be used to
translate the tail bounds from Lem. 4.2 (which in turn can be established
using a variant of [3, Lem. 1] under the additional assumption that the latent
factors are drawn i.i.d) to a bound on the distance between the latent
factors.

###  4.4 Proof of Lem. 4.2

We establish that the bound stated in 78(a) holds with probability at least
$1-\delta/2$. Proof of 78(b) can be derived in an analogous manner so that the
claim 78 holds with probability at least $1-\delta$ as desired.

Now to establish the bound 78(a), first we discuss our arguments without any
data-split. Define the shorthands

$\displaystyle V_{i,j}$
$\displaystyle\triangleq\mathopen{}\mathclose{{}\left[A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}}\right]_{t^{\prime}\neq
t}\in\mathbb{R}^{d\times
T},\quad\widetilde{\varepsilon}_{i,j}=\mathopen{}\mathclose{{}\left[\varepsilon_{i,t^{\prime}}-\varepsilon_{j,t^{\prime}}}\right]_{t^{\prime}\neq
t}\in\mathbb{R}^{T},\quad\widetilde{\varepsilon}^{(2)}_{i,j}\triangleq\mathopen{}\mathclose{{}\left[(\varepsilon_{i,t^{\prime}}-\varepsilon_{j,t^{\prime}})^{2}}\right]_{t^{\prime}\neq
t}\in\mathbb{R}^{T}$ (86) $\displaystyle\widehat{\Sigma}_{i,j}$
$\displaystyle\triangleq\frac{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}}{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}}\in\mathbb{R}^{d\times
d},\quad\widetilde{A}_{i,j}=\mathopen{}\mathclose{{}\left[A_{i,t^{\prime}}A_{j,t^{\prime}}}\right]_{t^{\prime}\neq
t}\in\\{0,1\\}^{T},\quad\text{and}\quad T_{i,j}\triangleq\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}.$ (87)

Then we also have that

$\displaystyle\widehat{\Sigma}_{i,j}=\frac{1}{T_{i,j}}V_{i,j}V_{i,j}^{\top}.$
(88)

We have

$\displaystyle\widehat{\rho}^{\textrm{unit}}_{i,t}(j)=\frac{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}(Y_{i,t^{\prime}}\\!-\\!Y_{j,t^{\prime}})^{2}}{\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}}$
$\displaystyle=\frac{1}{T_{i,j}}\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}(Y_{i,t^{\prime}}\\!-\\!Y_{j,t^{\prime}})^{2}$
(89)
$\displaystyle=\frac{1}{T_{i,j}}\bigg{[}(u_{i}-u_{j})^{\top}\mathopen{}\mathclose{{}\left(\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}}\right)(u_{i}-u_{j})$
(90) $\displaystyle\quad\qquad+\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}(\varepsilon_{i,t^{\prime}}-\varepsilon_{j,t^{\prime}})^{2}$
(91) $\displaystyle\quad\qquad+2(u_{i}-u_{j})^{\top}\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}(\varepsilon_{i,t^{\prime}}-\varepsilon_{j,t^{\prime}})\bigg{]}$
(92)
$\displaystyle=(u_{i}\\!-\\!u_{j})^{\top}\widehat{\Sigma}_{i,j}(u_{i}\\!-\\!u_{j})\\!+\\!\frac{1}{T_{i,j}}\mathopen{}\mathclose{{}\left\langle\widetilde{A}_{i,j},\widetilde{\varepsilon}^{(2)}_{i,j}}\right\rangle$
(93)
$\displaystyle\qquad\quad\\!+\\!\mathopen{}\mathclose{{}\left\langle\widetilde{\varepsilon}_{i,j},\frac{1}{T_{i,j}}V_{i,j}(u_{i}\\!-\\!u_{j})}\right\rangle$
(94)

We claim that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{\Sigma_{\textrm{time}}-\widehat{\Sigma}_{i,j}}}\right\|_{\mathrm{op}}\leq\frac{7c_{v}^{2}\sqrt{\log(9d/\delta)}}{p\sqrt{T}}\
\text{and}\ T_{i,j}\geq\frac{p^{2}T}{2}|\mathcal{V}}\right]\geq
1-\frac{\delta}{3}.$ (95)

if $p^{2}T\geq 12\log(9/\delta)$, which in turn is implied by the condition
24a. We defer the proof of this claim to the end of this section.

Note that $\widetilde{\varepsilon}^{(2)}_{i,j}-2\sigma^{2}\mathbf{1}$ is a
mean zero sub-Gaussian random vector with parameter $4c_{\varepsilon}^{2}$ and
$\widetilde{\varepsilon}_{i,j}$ is a mean zero sub-Gaussian random vector with
parameter $2c_{\varepsilon}$. Consequently, we have

$\textstyle\mathbb{P}\mathopen{}\mathclose{{}\left[|\mathopen{}\mathclose{{}\left\langle\widetilde{A}_{i,j},\widetilde{\varepsilon}^{(2)}_{i,j}}\right\rangle-2\sigma^{2}T_{i,j}|\leq
8c_{\varepsilon}^{2}\|{\widetilde{A}_{i,j}}\|_{2}\sqrt{\log(\frac{3}{\delta})}|\widetilde{A}_{i,j}}\right]$
$\textstyle\geq 1-\frac{\delta}{3},\ \text{and}\ $ (96)
$\textstyle\mathbb{P}\mathopen{}\mathclose{{}\left[|\mathopen{}\mathclose{{}\left\langle\widetilde{\varepsilon}_{i,j},\frac{1}{T_{i,j}}V_{i,j}(u_{i}\\!-\\!u_{j})}\right\rangle|\leq
4c_{\varepsilon}\frac{1}{T_{i,j}}\|{V_{i,j}(u_{i}-u_{j})}\|_{2}\sqrt{\log(\frac{3}{\delta})}|V_{i,j},\mathcal{U}}\right]$
$\textstyle\geq 1-\frac{\delta}{3}.$ (97)

Noting that $\|{\widetilde{A}_{i,j}}\|_{2}=\sqrt{T_{i,j}}$ and

$\displaystyle\frac{1}{T_{i,j}}\|{V_{i,j}(u_{i}-u_{j})}\|_{2}\leq\frac{1}{T_{i,j}}\mathopen{}\mathclose{{}\left\|{V_{i,j}}}\right\|_{\mathrm{op}}\mathopen{}\mathclose{{}\left\|{u_{i}-u_{j}}}\right\|_{2}\leq
2c_{u}\sqrt{\frac{\lambda_{\max}(\widehat{\Sigma}_{i,j})}{T_{i,j}}},$ (98)

respectively, we find that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\frac{1}{T_{i,j}}\mathopen{}\mathclose{{}\left\langle\widetilde{A}_{i,j},\widetilde{\varepsilon}^{(2)}_{i,j}}\right\rangle-2\sigma^{2}}\right|\leq
8c_{\varepsilon}^{2}\sqrt{\frac{\log(\frac{3}{\delta})}{T_{i,j}}}|\widetilde{A}_{i,j}}\right]$
$\displaystyle\geq 1-\frac{\delta}{3},\ \text{and}\ $ (99)
$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left\langle\widetilde{\varepsilon}_{i,j},\frac{1}{T_{i,j}}V_{i,j}(u_{i}\\!-\\!u_{j})}\right\rangle}\right|\leq
8c_{\varepsilon}c_{u}\sqrt{\frac{\lambda_{\max}(\widehat{\Sigma}_{i,j})\cdot\log(\frac{3}{\delta})}{T_{i,j}}}|V_{i,j},\mathcal{U}}\right]$
$\displaystyle\geq 1-\frac{\delta}{3}.$ (100)

Moreover, under 95, we also have
$\lambda_{\max}(\widehat{\Sigma}_{i,j})\leq\lambda_{\max}(\Sigma_{\textrm{time}})+c_{v}^{2}\leq
2c_{v}^{2}$ when $p^{2}T\geq 49\log(9d/\delta)$, which is implied by the
condition 24a. Putting the pieces together, we have

$\displaystyle|\widehat{\rho}^{\textrm{unit}}_{i,t}(j)-\rho^{\textrm{unit}}_{i,t}(j)|$
$\displaystyle\leq\bigg{|}(u_{i}\\!-\\!u_{j})^{\top}\widehat{\Sigma}_{i,j}(u_{i}\\!-\\!u_{j})\\!+\\!\frac{1}{T_{i,j}}\mathopen{}\mathclose{{}\left\langle\widetilde{A}_{i,j},\widetilde{\varepsilon}^{(2)}_{i,j}}\right\rangle\\!+\\!\mathopen{}\mathclose{{}\left\langle\widetilde{\varepsilon}_{i,j},\frac{1}{T_{i,j}}V_{i,j}(u_{i}\\!-\\!u_{j})}\right\rangle$
(101)
$\displaystyle\qquad\qquad-\mathopen{}\mathclose{{}\left\|{u_{i}\\!-\\!u_{j}}}\right\|^{2}_{\Sigma_{\textrm{time}}}\\!+\\!2\sigma^{2}\bigg{|}$
(102)
$\displaystyle\leq\mathopen{}\mathclose{{}\left|(u_{i}\\!-\\!u_{j})^{\top}(\widehat{\Sigma}_{i,j}-\Sigma_{\textrm{time}})(u_{i}\\!-\\!u_{j})}\right|\\!+\\!\mathopen{}\mathclose{{}\left|\frac{1}{T_{i,j}}\mathopen{}\mathclose{{}\left\langle\widetilde{A}_{i,j},\widetilde{\varepsilon}^{(2)}_{i,j}}\right\rangle-2\sigma^{2}}\right|$
(103)
$\displaystyle\qquad\qquad\\!+\\!\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left\langle\widetilde{\varepsilon}_{i,j},\frac{1}{T_{i,j}}V_{i,j}(u_{i}\\!-\\!u_{j})}\right\rangle}\right|$
(104)
$\displaystyle\leq\frac{28c_{u}^{2}c_{v}^{2}\sqrt{\log(9d/\delta)}+(12c_{\varepsilon}^{2}+16c_{u}c_{v}c_{\varepsilon})\sqrt{\log(\frac{3}{\delta})}}{p\sqrt{T}}$
(105)
$\displaystyle\leq\frac{56(c_{u}c_{v}+c_{\varepsilon})^{2}\sqrt{\log(9d/\delta)}}{p\sqrt{T}}$
(106)

conditional on $\mathcal{U},\mathcal{V}$, with probability at least
$1-\delta$.

Now we take into account the data-split. Since the data split is done using
independent Bernoulli coin flips, essentially same arguments work except that
we have to replace $T$ by $T/2$ and also have to account for the error between
$\Sigma_{\textrm{time}}$ and the matrix
$\frac{\sum_{t^{\prime}=1}^{T}B_{t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}}{\sum_{t^{\prime}=1}^{T}B_{t^{\prime}}}$
where $B_{t^{\prime}}$ are Bernoulli$(\frac{1}{2})$ flips. To that end, we can
simply apply 95 with $p^{2}$ replaced by $\frac{1}{2}$ so that we need to add
another error term
$\frac{14c_{u}^{2}c_{v}^{2}\sqrt{\log(18d/\delta)}}{\sqrt{T}}$ (so that the
appropriate tail bound holds with probability $1-\delta/2$). Adding this term,
replacing $\delta$ by $\delta/2$ and $T$ by $T/2$ in the error bounds above,
we conclude that

$\displaystyle|\widehat{\rho}^{\textrm{unit}}_{i,t}(j)-\rho^{\textrm{unit}}_{i,t}(j)|\leq\frac{94(c_{u}c_{v}+c_{\varepsilon})^{2}\sqrt{\log(18d/\delta)}}{p\sqrt{T}}$
(107)

conditional on $\mathcal{U},\mathcal{V}$, with probability at least
$1-\delta$.

Taking a union bound over all $j\neq i$, replacing $\delta$ by $\delta/(4N)$
and recalling the definition 23 of $e_{\delta}$ yields that the bound 78(a)
holds with probability at least $1-\delta/2$ as claimed in the beginning of
this proof.

#### Proof of 95

Our proof makes use of the standard Matrix Chernoff tail bounds. In
particular, we apply [9, Thm. 5.1.1] for the sequence
$\\{A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}\\}_{t^{\prime}\neq
t}$. Note that

$\displaystyle
A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}\succeq
0,\quad\|{A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}}\|_{\mathrm{op}}\leq
c_{v}^{2},\quad\text{and}\quad\mathbb{E}[\sum_{t^{\prime}\neq
t}A_{i,t^{\prime}}A_{j,t^{\prime}}v_{t^{\prime}}v_{t^{\prime}}^{\top}|\mathcal{V}]=p^{2}T\Sigma_{\textrm{time}}.$
(108)

Let $\lambda_{1}\triangleq\lambda_{\max}(\Sigma_{\textrm{time}})$. Then
$\mu_{\max}\triangleq\lambda_{\max}(p^{2}T\Sigma_{\textrm{time}})=p^{2}T\lambda_{1}$.
Moreover,
$\lambda_{\min}(p^{2}T\Sigma_{\textrm{time}})=p^{2}T\lambda_{\textrm{time}}$.
Substituting $\mu_{\max}=p^{2}T\lambda_{1}$,
$\mu_{\min}=p^{2}T\lambda_{\textrm{time}}$ and $R=c_{v}^{2}$ in [9, Thm.
5.1.1], we find that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\lambda_{\max}(V_{i,j}V_{i,j}^{\top})\geq(1+\epsilon)p^{2}T\lambda_{1}|\mathcal{V}}\right]$
$\displaystyle\leq
d\mathopen{}\mathclose{{}\left(\frac{e^{\epsilon}}{(1+\epsilon)^{1+\epsilon}}}\right)^{p^{2}T\lambda_{1}/c_{v}^{2}}\stackrel{{\scriptstyle(i)}}{{\leq}}de^{-\epsilon^{2}p^{2}T\lambda_{1}/(3c_{v}^{2})}\
\text{and}\ $ (109)
$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\lambda_{\min}(V_{i,j}V_{i,j}^{\top})\leq(1-\epsilon)p^{2}T\lambda_{\textrm{time}}|\mathcal{V}}\right]$
$\displaystyle\leq
d\mathopen{}\mathclose{{}\left(\frac{e^{-\epsilon}}{(1-\epsilon)^{1-\epsilon}}}\right)^{p^{2}T\lambda_{\textrm{time}}/c_{v}^{2}}\stackrel{{\scriptstyle(i)}}{{\leq}}de^{-\epsilon^{2}p^{2}T\lambda_{\textrm{time}}/(2c_{v}^{2})}$
(110)

for $\epsilon\in[0,1]$, where steps (i) and (ii) follow from the following
facts:

$\displaystyle\frac{e^{\epsilon}}{(1+\epsilon)^{1+\epsilon}}\leq
e^{-{\epsilon^{2}}/{3}}\quad\text{and}\quad\frac{e^{-\epsilon}}{(1-\epsilon)^{1-\epsilon}}\leq
e^{-{\epsilon^{2}}/{2}}\quad\text{for}\quad\epsilon\in[0,1].$ (111)

Moreover, [6, Thm. 4.1] implies that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[T_{i,j}\geq(1+\epsilon)p^{2}T}\right]\leq\mathopen{}\mathclose{{}\left(\frac{e^{\epsilon}}{(1+\epsilon)^{1+\epsilon}}}\right)^{p^{2}T}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:exp_fun_bounds}}}{{\leq}}e^{-\epsilon^{2}p^{2}T/3}.$
(112)

Setting all the tail bounds equal to $\delta/9$, we conclude that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\|{V_{i,j}V_{i,j}^{\top}-p^{2}T\Sigma_{\textrm{time}}}\|_{\mathrm{op}}\leq
c_{v}^{2}p\sqrt{3T\log(9d/\delta)}}\right]$ $\displaystyle\geq 1-2\delta/9\
\text{and}\ $ (113)
$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[|T_{i,j}-p^{2}T|\leq
p\sqrt{3T\log(9/\delta)}}\right]$ $\displaystyle\geq 1-\delta/9,$ (114)

where we have also used the fact that
$\lambda_{\textrm{time}}\leq\lambda_{1}\leq c_{v}^{2}$ to simplify the
constants in 113. We can directly verify that $T_{i,j}\geq p^{2}T/2$ under the
event in the display 114 if $p^{2}T\geq 12\log(9/\delta)$. Next, we note that

$\displaystyle\widehat{\Sigma}_{i,j}-\Sigma_{\textrm{time}}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:sig_V_reln}}}{{=}}\frac{1}{T_{i,j}}V_{i,j}V_{i,j}^{\top}-\Sigma_{\textrm{time}}=\frac{V_{i,j}V_{i,j}^{\top}-p^{2}T\Sigma_{\textrm{time}}}{T_{i,j}}+\frac{T_{i,j}-p^{2}T}{T_{i,j}}\Sigma_{\textrm{time}}.$
(115)

Noting that $\|{\Sigma_{\textrm{time}}}\|_{\mathrm{op}}=\lambda_{1}\leq
c_{v}^{2}$, under the events in the display 113 and 114, we find that

$\displaystyle\|{\widehat{\Sigma}_{i,j}-\Sigma_{\textrm{time}}}\|_{\mathrm{op}}$
$\displaystyle\leq\frac{\|{V_{i,j}V_{i,j}^{\top}-p^{2}T\Sigma_{\textrm{time}}}\|_{\mathrm{op}}}{T_{i,j}}+\frac{|T_{i,j}-p^{2}T|}{T_{i,j}}\|{\Sigma_{\textrm{time}}}\|_{\mathrm{op}}$
(116)
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\frac{c_{v}^{2}p\sqrt{3T\log(9d/\delta)}}{p^{2}T/2}+\frac{p\sqrt{3T\log(9/\delta)}}{p^{2}T/2}c_{v}^{2}$
(117) $\displaystyle\leq\frac{7c_{v}^{2}\sqrt{\log(9d/\delta)}}{p\sqrt{T}},$
(118)

where the inequality (i) holds with probability at least $1-\delta/3$ given
the arguments above. The proof is now complete.

##  5 Proof of Corollary 1: Error improvement for examples

We prove each part separately.

#### Proof of part (a)

We have
$\eta_{2}^{\prime}+2\frac{e_{\delta}}{p\sqrt{T}}=\eta_{1}^{\prime}+2\frac{e_{\delta}}{p\sqrt{N}}=4e_{p,\delta,N}$
and hence as noted in Example Example, we have
$S_{i,\eta_{1}^{\prime}}\geq\frac{N}{2M_{\textrm{unit}}}=\frac{N}{2M}$ and
$S_{t,\eta_{2}^{\prime}}\geq\frac{T}{2M_{\textrm{time}}}=\frac{N}{2M}$ since
$\eta_{1}^{\prime},\eta_{2}^{\prime}=\Omega(N^{-\frac{1}{2}})$. Substituting
these quantities in the bound from Thm. 3.1, we find that

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}-\theta_{i,t})^{2}\leq
c\frac{e_{\delta}^{2}}{p^{2}N}+\frac{c\log(N^{2}/\delta)\sigma^{2}}{p}\\!\bigg{(}\\!\frac{2M}{N}+\frac{M^{2}}{p^{2}N^{2}}\bigg{)}=\widetilde{\mathcal{O}}\mathopen{}\mathclose{{}\left(\frac{M}{N^{1-\beta}}+\frac{M^{2}}{N^{2-3\beta}}}\right),$
(119)

with the claimed probability. Our result follows.

#### Proof of part (b)

We first derive a general bound for $p=\Theta(N^{-\beta})$ the choice
$\eta_{1}=\eta_{2}=2\sigma^{2}+e_{p,\delta,N}+r$ for
$r=\Omega(e_{p,\delta,N})$ so that
$\eta_{1}^{\prime}=\eta_{2}^{\prime}=\Theta(r)$. As noted in Example Example,
we then have

$\displaystyle
S_{i,\eta_{1}^{\prime}}=\Theta(Nr^{\frac{d}{2}})\quad\text{and}\quad
S_{t,\eta_{2}^{\prime}}=\Theta(Nr^{\frac{d}{2}}).$ (120)

Substituting these quantities in the bound from Thm. 3.1, we find that

$\displaystyle(\widehat{\theta}_{i,t,\boldsymbol{\eta}}^{\textrm{DR}}-\theta_{i,t})^{2}$
$\displaystyle\leq
cr^{2}+c\frac{re_{\delta}}{p\sqrt{N}}+c\frac{e_{\delta}^{2}}{p^{2}N}+\frac{c}{pNr^{\frac{d}{2}}}+\frac{c}{p^{3}N^{2}r^{d}}$
(121)
$\displaystyle=\widetilde{\mathcal{O}}\mathopen{}\mathclose{{}\left(r^{2}+\frac{r}{N^{\frac{1-2\beta}{2}}}+\frac{1}{N^{1-\beta}r^{\frac{d}{2}}}+\frac{1}{N^{2-3\beta}r^{d}}}\right).$
(122)

with probability at least
$1-2\delta-2Ne^{-cN^{1-2\beta}r^{d/2}}-e^{-cN^{2-3\beta}r^{d}}$. Let $f(r)$
denote the function inside the parentheses in the last display. We have

$\displaystyle r=N^{-\frac{1-2\beta}{2}}\quad\Longrightarrow\quad g(r)$
$\displaystyle=\mathcal{O}(\max\\{N^{-(1-2\beta)},\ N^{(1-2\beta)d-4},\
N^{(\frac{1}{2}-\beta)d-(1-3\beta)}\\})$ (123) $\displaystyle
r=N^{-\frac{2(1-\beta)}{d+4}}\quad\Longrightarrow\quad g(r)$
$\displaystyle=\mathcal{O}(\max\\{N^{-\frac{4(1-\beta)}{d+4}},\
N^{-\frac{4+\frac{d}{2}-(6+d)\beta}{d+4}},\
N^{-\frac{8-\beta(d+12)}{d+4}}\\})$ (124) $\displaystyle
r=N^{-\frac{2-3\beta}{d+2}}\quad\Longrightarrow\quad g(r)$
$\displaystyle=\mathcal{O}(\max\\{N^{-\frac{2(2-3\beta)}{d+2}},\
N^{-\frac{\frac{d+5}{2}(1-2\beta)+\frac{1}{2}}{d+2}},\
N^{-\frac{2(1-\beta)+\frac{d\beta}{2}}{d+2}}\\}).$ (125)

For any given choice of $\beta\in[0,\frac{1}{2}]$, we can choose an $r$ that
obtains the minimum value for $g$ to obtain a suitable bound. For $\beta=0$,
this choice is $r=N^{-\frac{2(1-\beta)}{d+4}}$ with $\beta=0$, which in fact
yields the claimed result in the corollary.

[Acknowledgments] The authors thank Avi Feller and James M. Robins for their
helpful comments.

## References

  * [1] [author] Agarwal, AnishA., Shah, DevavratD. and Shen, DennisD. (2021). Synthetic interventions. arXiv preprint arXiv:2006.07691.
  * [2] [author] Arkhangelsky, DmitryD., Athey, SusanS., Hirshberg, David AD. A., Imbens, Guido WG. W. and Wager, StefanS. (2019). Synthetic difference in differences Technical Report, National Bureau of Economic Research.
  * [3] [author] Dwivedi, RaazR., Murphy, SusanS. and Shah, DevavratD. (2022). Counterfactual inference for sequential experiments. arXiv preprint arXiv:2202.06891.
  * [4] [author] Lee, Christina E.C. E., Li, YihuaY., Shah, DevavratD. and Song, DogyoonD. (2016). Blind Regression: Nonparametric Regression for Latent Variable Models via Collaborative Filtering. In Advances in Neural Information Processing Systems 29 2155–2163.
  * [5] [author] Li, YihuaY., Shah, DevavratD., Song, DogyoonD. and Yu, Christina LeeC. L. (2019). Nearest neighbors for matrix estimation interpreted as blind regression for latent variable model. IEEE Transactions on Information Theory 66 1760–1784.
  * [6] [author] Motwani, RajeevR. and Raghavan, PrabhakarP. (1995). Randomized algorithms. Cambridge university press.
  * [7] [author] Neyman, JerzyJ. (1923). Sur les applications de la theorie des probabilites aux experiences agricoles: Essai des principes. Master’s Thesis.
  * [8] [author] Rubin, Donald B.D. B. (1976). Inference and Missing Data. Biometrika 63 581–592.
  * [9] [author] Tropp, Joel AJ. A. (2012). User-friendly tools for random matrices: An introduction. Technical Report.