The Ordinary Least Eigenvalues Estimator

The Ordinary Least Eigenvalues

Introduction

Linear regression models with individual-specific effects are widely used to fit data with network structures. Such linear models were used to explain trade flows between countries (AndersonWincoop2003, Fally2015), to fit matched employer-employee data (AKM, Bonhommeetal2019), or to study teacher effect on student performance (Jacksonetal2014), to mention a few examples. In applications, these linear regression models are most often used with a particular specification of the individual-specific effect. A popular specification consists of including the individual effects as additive error terms. This is the approach taken for instance in AKM and Jacksonetal2014. Broadly, three types of estimators are used under this specification. The two way fixed effects estimator (AKM) exploits the additivity of the model to eliminate the individual effects. After double differencing, the initial model is turned into a regular linear regression model (free of the individual effects), and estimators are obtained by least squares on the transformed model. When the data is non bipartite with a number $N$ of agents (respectively, when it is bipartite, with $N$ and $M$ agents on each side), the two way fixed estimator of the slope parameters converges at the optimal rate of $N$ (resp. $\sqrt{NM}$). The two way fixed estimator comes with a significant caveat: the slope parameters on any agent-specific observable covariates disappear in the double differencing process, in the same way as the individual effects. Those can be recovered in a second stage by ordinary least squares if we further assume the individual effects to be exogenous with respect to the additive observable. The second stage OLS estimators for the slope parameters on the additive covariates converges at a $\sqrt{N}$ rate.

A second approach appeals to the standard OLS estimator (e.g. Rose2004, FafchampsGubet2007). In the dydadic linear regression setting, the OLS estimator is in general $\sqrt{N}$ consistent for all the parameters. Given that for some covariates the two way fixed effects estimator can provide $N$-consistent estimators, the OLS estimator is severely inefficient. Other approaches consist of estimating a fixed effects model, by regressing the output variable on the covariates, individual indicators and interaction of individual interactions. As such, this procedure suffers from the incidental parameter problem, since the number of parameters to be estimated grows faster than the sample size. The incidental parameter problem is usually resolved through a dimension reduction, for instance by “grouping" individuals into a small number of categories that share the same individual effect, or simimilarly by assuming the individual effects are drawn from some discrete distribution (BonhommeManresa2015, Bonhommeetal2019).

This paper proposes an $N$-consistent estimator for the slope parameters under an exogeneity assumption on the individual (unobserved) effects. We exploit the matrix structure of network data and identify the individual effects' footprint on the spectrum of the output matrix. We then correct for the unobservables'effect on the spectrum. The correction proposed in this paper consists of tweaking the standard minimization problem that yields the OLS estimator in a way that mechanically attenuates the effect of the unobservables on the explained variable's spectrum.

This paper contributes to the literature on the dyadic linear regression model in a few ways. First, we provide a rate optimal estimator for all slope parameters. In the process, we provide valid estimators for the individual effects and overcome the incidental parameter issue without any appeal to discretization. This is an improvement over existing estimation procedures that are either rate inferior (e.g. OLS), or that only provide rate optimal estimators for a subset of slope parameters (e.g. the two-way fixed effects). This comes at the cost of the exogeneity assumption the paper imposes on the individual effects. Second, to the best of my knowledge, this is the first time tools from the random matrix literature are used for inference. These tools allow for the analysis of a new set of estimators that are based on a matrix representation of the regression model. These new estimators are severely non linear and are defined as extrema of globally ill-behaved objective functions. We show that in a simple sub-specification, our estimator is semi-parametrically efficient, unlike OLS or any other linear estimator. This suggests that achieving semi parametric efficiency in models with dyadic dependence might require the recourse to non standard estimators.

The next section introduces the setup and lays out the main intuitions leading up to the definition of the new estimator. Section (ref) discusses the estimator's theoretical properties and numerical implementation. Section (ref) proposes estimators for the asymptotic bias and covariance matrix. Finally, section (ref) shows the results from Monte Carlo simulations. All proof are differed to the end of the paper (section (ref)).

The eigenvalue corrected OLS

Consider the model:

for all $i\neq j$, where $A_i$'s are i.i.d centered random variables with finite fourth moments. The $V_{ij}$'s are i.i.d centered square integrable random variables with $V_{ij}=V_{ji}$, $\beta_0:=(\beta_{0,1},.., \beta_{0,L})$ is the parameter of interest and $\gamma \geq 0$ and $\delta \in \{ -1, +1\}$ are unknown nuisance parameters.\footnote{$\gamma$ is set to be positive because the model (ref) can also be re-expressed as$Y_{ij}= X_{ij}\beta +(-\gamma)((-A_i)+(-A_j))+ \delta\times (-A_i) \times (-A_j)+V_{ij}$. The sign of $\gamma$ is not identified.} The covariates $X$ are such that for all $i,j,l$ such that $i\neq j$: $X_{ij,l}= X_{ji, l}=\phi(X_i, X_j, W_{ij})$, for some (unknown) function $\phi$, i.i.d random variables $X_i$ and $i.i.d.$ variables $W_{ij}$. By convention, $X_{ii,l}=0$ for all $i$ and $l$ by convention.

The model in equation ((ref)) can also be re-expressed:

which reduces the study of the model (ref) to that of the model:

where the errors $U_i:=\gamma+A_i$ are no longer assumed to be centered. All the slope parameters remain unchanged as you move from model ((ref)) to ((ref)) (or ((ref))), only the intercept is altered by the correction term“$-\delta \gamma^2$" in equation ((ref)). Therefore, any “good" estimators for the parameters of the model ((ref)) also provide good estimators for the parameters in models ((ref)) and ((ref)), except perhaps for their intercepts. For reasons that will soon become clear, the rest of the paper focuses on the model ((ref)). This paper focuses on the estimation of the slope parameter and does not propose an improved estimator for the intercept.

Let $N$ be the sample size (number of nodes or agents $i$). Denote $Y$ and $V$ the $N\times N$ matrices with entries $Y_{ij}$, $V_{ij}$ and $X_l$ the matrix with entries $X_{ij,l}$ for every $l=1..L$. $Y$ and $X_l$'s diagonal entries are equal to zero. $V$'s $i$th diagonal term is equal to $\delta(E(U_1^2)-U_i^2)$. Finally, stack the individual random effects into a vector denoted $U$. This allows for the formulation of model (ref) in a compact matrix form :

Looking at model ((ref)) through the lens of the matrix formulation ((ref)) allows for novel interpretations of classical estimators. It also gives access to potentially interesting estimators that are based on standard matrix functions. To illustrate, consider the ordinary least squares estimator on model (ref), defined by $$ \hat{\mu}_{OLS}:= \arg \min_{\mu \in \mathbb{R}^L} \sum_{i,j}\left(Y_{ij}-\sum_{l=1}^L \mu_{l} X_{ij,l}\right)^2$$ Under the formulation ((ref)), $ \hat{\mu}_{OLS}$ can also be expressed:

where, for any $N\times N$ matrix $M$, $Trace(M)$ denotes $M$'s trace, $\lambda_1(M)\geq \lambda_2(M)\geq \dots \lambda_N(M)$ are $M$'s eigenvalues ranked from largest to smallest and $M( . )$ is the matrix valued function that takes an element $\mu \in \mathbb{R}^L$ and returns the matrix

Equation ((ref)) indicates that the OLS estimator can also be defined as a minimizer of the average squared eigenvalues of the matrix $M(\mu)$. Let's examine the distribution of $M(\mu)$'s eigenvalues for values of $\mu$ that are “close" to the true value $\mu_0$, assuming $\delta=1$ (the treatment for $\delta=-1$ is similar). Begin with the value $\mu=\mu_0$, that is, let's look at the distribution of the eigenvalues of the matrix $UU'+V$. \footnote{ We ignore the effect of the matrix $E(U_1^2)I_N$ in the discussion that follows. $E(U_1^2)I_N$ simply shifts all eigenvalues by the same quantity $E(U_1^2)$. The shift size will turn out to be of a low order of magnitude compared to the bulk of $UU'+V$'s eigenvalues and its effect will be negligible anyways.} Figure (ref) shows the histogram of the eigenvalues of the simulated matrix $\frac{1}{\sqrt{N}}\left( UU'+V\right)$, where the $U$'s and $V$'s are i.i.d standard normal and the sample size is set to $N=2000$.

The histogram in figure (ref) shows two distinct parts: to the left, a block of eigenvalues concentrated between values $\sim -2$ and $\sim +2$, and a single eigenvalue, further to the right, at around value $\sim 46$. After proper rescaling (and ignoring the single eigenvalue to the left for the rescaled histogram to fit on a page) the block of eigenvalues to the left has the shape of a semi-circle as shown in figure (ref) .

To rationalize the shape of the histogram (ref), let's examine the eigenvalues of each of the terms composing $M(\mu_0)$. The matrix $UU'$ is of rank 1, its unique non null eigenvalue is equal to $U'U=\sum_i U_i^2$ which is of the same order as $N E(U_1^2)$ when $N$ is large enough.

Figure (ref) shows the histogram of $V$'s eigenvalues. The two histograms in (ref) and (ref) are seemingly identical. Only $M(\mu_0)$'s outlier eigenvalue (the one approximately equal to 46) is absent from $V$'s histogram. This should come as no surprise: the matrix $M(\mu_0)$ is a rank 1 deformation of $V$. The impact of rank 1 deformations on the eigenvalues of the original matrix ($V$ here) is well studied (e.g. Bunch1978). Because $UU'$s unique eigenvalue is positive, modifying $V$ through $UU'$ shifts all of $V$'s eigenvalues upwards such that $V$'s eigenvalues are interlaced with $V+UU'$'s, that is, for $i=2, ..., N$: \[\lambda_{i}(V) \leq \lambda_i\left(V+UU'\right)\leq \lambda_{i-1} (V) \] and \[ \lambda_{1} (V) \leq \lambda_1\left(V+UU'\right)\] Provided that $V$'s eigenvalues (rescaled by $\frac{1}{\sqrt{N}}$) are concentrated roughly between -2 and 2, then the inequalities above predict that $V+UU'$'s $N-1$ smallest eigenvalues will be only shifted by a small amount, which explains why the figures (ref) and (ref) are not visually distinguishible.

The semi-circle in figure (ref) is reminiscent of Weigner's semi-circle law in the random matrix literature (see for instance FlorentAntti2018). Weigner's law states that the empirical distribution of the eigenvalues of a random symmetric matrix with centered square integrable entries “converges" (in a sense that is made precise below) to a distribution with a semi-circular probability density function . Formally:\\

Theorem [FlorentAntti2018]. Let $W_N$ be a sequence of $N\times N$ symmetric matrices with i.i.d centered square integrable entries below the diagonal with $E(W_{N, 12}^2)=\sigma^2$, and i.i.d centered square integrable entries on the diagonal. Let $Z_N$ be a rescaling of $W_N$: $Z_N:=\frac{1}{\sqrt{N}}W_N$. Let $\mathcal{S}$ be the probability distribution characterized by the the probability density function $f$: $$ f(x) := \left\{

\right. $$ Then for any measurable set $I$: $$ \frac{\sum_{i=1}^N \mathbbm{1}\{ \lambda_i(Z_N)\in I\}}{N}\rightarrow_p \mathcal{S}(I) \\$$ }

In particular, Furedi1981 show that $V$'s largest eigenvalue is of order $\sqrt{N}$ with probability approaching 1 as $N$ grows.

These observations combined suggest the following rough interpretation of the histogram (ref): $M(\mu_0)$'s $N-1$ smallest eigenvalues are of order $\sqrt{N}$ and are "very close" to $V$'s eigenvalues, whereas the largest eigenvalue is due to the $UU'$ deformation and is of order $N$.

Let's extend these intuitions to values of $\mu$ that are different from the true parameter $\mu_0$. If $\mu$ is too far from $\mu_0$, then the term $ \sum_{l=1}^L (\mu_{0,l} - \mu_l) X_{l}$ in equation ((ref)) can become dominant and dwarf the contributions of $V$ and $UU'$ in $M(\mu)$'s eigenvalue distribution. In the other extreme, when the candidate $\mu$ is“very close" to $\mu_0$, then the contribution of the covariates' term becomes negligible and we obtain a histogram that is similar to the one in figure (ref).

The values of $\mu$ that are abberantly far from $\mu_0$ lead to the eigenvalues of $ \sum_{l=1}^L (\mu_{0,l} - \mu_l) X_{l}$ being of a higher than order $\sqrt{N}$. Subsequently, they are easy to eliminate as they produce a histogram that is grossly different from the one in figure (ref). However, this rough discrimination strategy will be ineffective for values of $\mu$ that return a term $\sum_{l=1}^L (\mu_{0,l} - \mu_l) X_{l}$ of order $\sqrt{N}$ or lower. In any case, for model ((ref)), the OLS estimator is known to be $\sqrt{N}$ consistent in general (See for instance Menzel2021 or section 4 in Graham2020). The following lemma shows that any $\sqrt{N}$ estimator is in fact “close enough" for our purposes.\footnote{Because the intercept is shifted by $-\delta \gamma^2$ when we move from the original model ((ref)) to model ((ref)), the OLS estimator of the intercept would need to be corrected to account for the shift. That is done in proposition (ref) and relegated to the appendix.}

The lemma implies that any initial estimator $\tilde{\mu}$ that is $\sqrt{N}$ consistent - like the OLS estimator - would yield a covariates' term such that $\sum_{l=1}^L (\mu_{0,l} - \tilde{\mu}_l) X_{l}=O_p(\sqrt{N})$. It would produce an eigenvalue histogram for $M(\tilde{\mu})$ that is similar to figure (ref) with one outlier eigenvalue of order $N$, due to the rank 1 modification $UU'$, and a cloud of eigenvalues that are of a smaller order $\sqrt{N}$ but that need not form a semi-circle this time.

Provided that the candidate $\mu$ is close enough to $\mu_0$, the largest eigenvalue of $M(\mu)$ is, at least up to a first order approximation, closely tied to the error term $UU'$. Notice that, absent the $UU'$ from the model (ref) (or the random effects $A_i$ and $A_j$ from the model ((ref))), we would be back to the standard linear regression model with i.i.d. and exogenous noise $V_{ij}$. In that case, we know that OLS is efficient, and since the sample size is $\frac{N(N-1)}{2}$, the rate of convergence of the OLS estimator would be $N$, rather than $\sqrt{N}$ under models ((ref)) or ((ref)).

An appealing idea is then to modify the objective function in the matrix form definition of the OLS in ((ref)) to remove the contribution of the random effects. Following the intuition laid down so far, this can for instance be done by removing $M(\mu)$'s largest eigenvalue from the sum of squared errors before minimizing. The new estimator would be a solution to the minimization problem

First, let's show that the problem ((ref)) admits at least a solution. Let $f_N$ the function $$f_N: \mu \rightarrow \left( \sum_{i\neq j} X_{ij}'X_{ij}-\sum_{i\neq j, k\neq i,j} \nu_i({\mu}) \nu_j({\mu}) X_{jk}'X_{ik} \right)^{-1} \left( \sum_{i\neq j}X_{ij}'Y_{ij} -\sum_{i\neq j, k\neq i,j} \nu_i({\mu}) \nu_j({\mu}) X_{jk}'Y_{ik} \right)$$

The condition on $E(X_{12}'X_{12})$ is standard in the classical least squares theory and insures that in the population, none of the regressors is a linear combination of the others (see for instance Wooldridge2010, chapter 4).

The optimization problem ((ref)) involves functions that are in general not smooth. It is not solvable in closed form. In addition to guaranteeing the existence of a solution, lemma (ref) provides a practical tool to study the behavior of estimators obtained through the optimization problem ((ref)). Intuitively, equation ((ref)) is a first order condition of a minimization problem that is equivalent to ((ref)). Let $\mu^* \in \arg \min_{\mu \in \mathbb{R}^L} \sum_{i=2}^N \lambda_i \left(M(\mu)\right)^2 $ and note {

}

The last equality allows for the expression of ${\mu}^*$ as the minimizer of a smooth and convex function over $\mathbb{R}^L$ (in fact, strictly convex with probability 1, when $N$ is large enough): $$\mu \rightarrow \sum_{i\neq j} \left(Y_{ij}-\sum_{l=1}^L \mu_{l} X_{ij,l}\right)^2- \sum_{i\neq j,k\neq i,j}\nu_i({\mu}^*) \nu_j ({\mu}^*) \left(Y_{ik}-\sum_{l=1}^L \mu_{l} X_{ik,l}\right)\left(Y_{kj}-\sum_{l=1}^L \mu_{l} X_{kj,l}\right)$$ the first order condition results in the fixed point problem ((ref)). The proof in section (ref) closely follows this sketch.

Lemma (ref) does not guarantee the uniqueness of the solution to the minimization problem ((ref)). The iteration process just described, when it converges, could converge to one of many potential fixed point of ((ref)) (solutions to ((ref))). Additionally, the iteration process could be explosive, leading the iterations to diverge rather than approach one of the fixed points. The function $f_N$ is generally ill-behaved. In general, it is neither convex, nor quasi-convex, nor differentiable. Figure (ref) illustrates $f_N$'s behavior for the simplest model nested in model ((ref)): $Y_{ij}=\mu_0 +U_iU_j+V_{ij}$ for $\sigma_U=\sigma_V=1$, $\mu_0=1$ and for $N=100$. In this example, $f_N$ is convex between $\approx -0.5$ and $\approx 2$, it has a point of inflexion, smoothly switching convexity at $\approx -0.5$. $f_N$ is not derivable at $\approx 2$. However, $f_N$ has a unique minimum (on the interval displayed in figure (ref)), that is close to the true parameter $\mu=\mu_0=1$. The figure also points to the direction that the results in the sequel will follow: I show that with high probability, $f_N$ is well behaved in a shrinking neighborhood of $\mu_0$ . $\mu_0$ being unknown, knowledge of a good enough first stage estimator will be essential throughout the paper. In particular, we study the estimator defined defined in ((ref)) by studying single successive iterations on the fixed point problem ((ref)). It turns out that when the iteration process is initiated with a good first stage estimator, e.g. OLS, it converges to a fixed point or a minimizer (formal statements are presented in corollary (ref) in the following section). In fact, the main reason why the errors $U$ are required to be exogenous throughout the paper is that without exogeneity, no valid first stage estimators are available to the best of my knowledge.

The estimator(s) studied in this paper are obtained by iterating equation ((ref)), that is, by plugging some “reasonable" initial candidate estimator in the right hand side of ((ref)) to obtain what we show is a more precise estimator on the left hand side, then iterating this process as needed until the true fixed point distribution is achieved. We show that when $\gamma$ in model ((ref)) is equal to zero, a single iteration starting with the ordinary least squares (or any $\sqrt{N}$-consistent first stage) is enough to achieve full efficiency. This is reminiscent of “one-step theorems" in standard cross-section models where, knowing that efficiency is achieved by an estimator minimizing some objective function, full efficiency can also be achieved through a single iteration of a gradient descent algorithm applied to the objective function provided that the algorithm is initiated with an appropriate first stage estimator (see NeweyMcFadden1994 section 3.4 for a detailed discussion). As in the standard cross-section setting, iterating beyond the first iteration has no first order effect.

For non-null $\gamma$ however, any $\sqrt{N}$-consistent first stage can only return a $\sqrt{N}-$ consistent estimator following any finite number of iterations, even though the asymptotic variance of the generated estimators decays exponentially in the number of iterations. For this case, the next section lays out a general method to obtain a rate $N$ estimator.

The estimator's asymptotic behavior

The first result examines a single iteration of the fixed point problem ((ref)).

A detailed proof is presented in section (ref). Proposition (ref) shows that $\left( E(X_{12}X_{12}')-\frac{E(U_1)^2}{E(U_1^2)}E(X_{12}X_{23}')\right)$ is invertible, insuring that $K$ is well defined.

Equation ((ref)) describes how the distribution of the single iteration estimator relates to the first stage estimator's. An immediate corollary of proposition ((ref)) is that the single iteration estimator $\hat{\mu}$ is consistent and converges to $\mu_0$ at least as at a rate of ${\sqrt{N}}$. Also, up to a first order approximation, the single stage estimator depends linearly on the initial $\tilde{\mu}$.

Whether the iteration process improves the quality of estimation depends on the matrix $K$. When $E(U_1)=0$ in proposition (ref) (i.e. when $\gamma=0$ in model ((ref))), the matrix $K$ is null and equation ((ref)) becomes $$\hat{\mu}-\mu_0=O_p\left( \frac{1}{N}\right).$$ After a single iteration, we are able to achieve the optimal rate of convergence $N$. Unfortunately, proposition (ref) does not provide the asymptotic distribution of $\hat{\mu}$ or the effect iterations have beyond the first iteration. To answer both these questions, we need to zoom into the $O_p\left( \frac{1}{\sqrt{N}}\right)$ term in equation ((ref)) and determine how it depends on the first stage estimator $\tilde{\mu}$ and/or how it behaves asymptotically. The next proposition and its proof in section (ref) address this case.

A “one step theorem" applies, one iteration is enough to achieve full efficiency. The argument proving the efficiency of $\hat{\mu}$ in proposition (ref) is simple: consider the alternative model $ Y_{ij}=\sum_{l=1}^L \mu_{0,l} X_{ij,l}+ V_{ij}=X_{ij}\mu_0+V_{ij}$ is $i.i.d.$ with $i.i.d.$ errors $V_{ij}$. In this model, the ordinary least squares estimator is known to be efficient and asymptotically normal, with asymptotic covariance matrix $2\sigma_v^2 E(X_{12}X_{12}')^{-1}$ - the same asymptotic distribution as in ((ref)) (see for instance Chamberlain1987 or Newey1990). Given that our model of interest ((ref)) is noisier than the alternative model, the following corollary holds.

When $K\neq 0$, the one step theorem no longer applies. After any finite number of iterations, the new estimator is still $\sqrt{N}$-consistent. To understand the role of $K$ when $K\neq 0$, consider the simple case where we have a single regressor ($L=1$). $K$ becomes a scalar and when $|K|<1$, $\hat{\mu}$ is to a first order closer to $\mu_0$ than $\tilde{\mu}$. If the first stage estimator is asymptotically normal (the standard ordinary least squares estimator for example) with an asymptotic variance of $\tilde{\sigma}^2$, then $\hat{\mu}$ is normally distributed with variance $K^2 \tilde{\sigma}^2 <\tilde{\sigma}^2$. Moreover, as we iterate, the variance decays exponentially in the number of iterations. Conversely, if $|K|>1$, iterations produce noisier estimators, and the variance explodes exponentially with the number of iterations. Finally, if $|K|=1$, then the new estimator is asymptotically equivalent to the first stage estimator, iteration is neither useful nor harmful.

Simplify further, and assume that the single regressor is in fact just a constant $X_{ij}=1$, that is, we are interested in estimating the mean of $Y_{ij}$. The constant $K$ becomes $K=\frac{E(U_1)^2}{E(U_1^2)}$ which is positive and strictly smaller than 1 (since by assumption $\sigma_U^2>0$), the iterations improve estimation quality.

When $L>1$, $K$ is a matrix. Rather than comparing $K$ to 1, the relevant comparison is now between $K$ and $I_L$ - the identity matrix of dimension $L$ - in the partial order on symmetric matrices. When $K^2> I_L$, that is, when $K$'s eigenvalues are all larger than 1 in absolute value, the successive iterations follow an explosive path of covariance matrices. The conclusions are similar to the univariate setting in the two cases: $K^2< I_L$ or $K^2=I_L$. In the multivariate case however, these three cases are not exhaustive, since $>$ here is only a partial order. Fortunately, the next proposition shows that the only possibile case, given our assumptions, is in fact $0<K<I_L$.

Together, the propositions (ref) and (ref) imply that given a $\sqrt{N}$-consistent initial estimator and a fixed $\epsilon>0$, we iterate the process described in the equation ((ref)) to obtain a new $\sqrt{N}$-consistent estimator with a variance that is smaller than $\epsilon$ (or $\epsilon I_L$ in the multivariate case). This strongly suggests that an estimator with a faster than $\sqrt{N}$ rate of convergence exists. In fact, using a simple trick, the propositions (ref) and (ref) provide a rate $N$ (a rate optimal) estimator.\footnote{A similar idea is used for instance to construct a “Generalized Jackknife estimator" (e.g. PowellStockStoker1989,CattaneoCrumpJansson2013). In the context of proposition (ref), the term $\sqrt{N}(\tilde{\mu}-\mu_0)$ is eliminated by taking the convex combination, in the same fashion that the bias is removed in the generalized jacknife by taking a convex combination of estimator with the same bias. }

Another corollary of proposition (ref) is that, if $f_N$ has a fixed point $\hat{\mu}^*$ that is $\sqrt{N}$-consistent, then equation ((ref)) yields: $$(I-K)\sqrt{N}(\hat{\mu}^*-\mu_0)=O_p\left(\frac{1}{\sqrt{N}}\right) $$ so $\hat{\mu}^*$ is in fact $N$-consistent. Proposition (ref) is silent about the exact asymptotic distribution of $\hat{\mu}^*$ and about its existence.

To establish the existence of a fixed point, notice that equation ((ref)) has the flavor of Taylor expansion, where the matrix $K$ would represent a gradient. Because the matrix $K$ has a spectral radius that is smaller than 1 (proposition (ref)), then $f_N$ must be contracting in a local sense. Then (a variation on) the Banach fixed point theorem should prove existence. This intuition is the main idea for the proof for the next corollary.

So $ \hat{\mu}^*$ exists with probability approaching 1 and is rate optimal. It is left to determine its asymptotic distribution. We need to compute a higher order term in the expansion ((ref)) of proposition (ref). That is the purpose of proposition (ref).

Because of the presence of the residual $R_N$ in equation ((ref)), the new expansion is fundamentally different from the previous one (((ref)) in proposition (ref)). The effect of an iteration on the estimation quality is now ambiguous and depends on how the first stage estimator $\tilde{\mu}$ relates to the residual $R_N$. Even if they were independent, it is not clear whether iteration improves estimation. Unfortunately, even though proposition (ref) provides the asymptotic distribution of $R_N$, that is not enough to fully characterize the distribution of the single iteration estimator. For that, we would need the joint distribution of the first stage $\Tilde{\mu}$ and $R_N$, which is challenging even for a single iteration. However, we can see that because $K$ has a smaller than one spectrum (by proposition (ref)), as we iterate, the contribution of the initial (first stage or input) estimator fades away. Intuitively, starting with some $N-$ consistent first stage $\tilde{\mu}(=: \hat{\mu}_0)$, from ((ref)):

So the limit distribution (when $m$ approaches infinity) should not depend on the initial estimator $\tilde{\mu}$. Corollary (ref) formalizes these thoughts.

Notice that $\hat{\mu}^*$ is asymptotically biased. In the next section, we offer a correction to this bias by proposing a consistent estimator for the bias term.

The corollary show that if we initiate a sequence $\hat{\mu}_0:=\tilde{\mu}$, for some initial $\sqrt{N}$-consistent estimator $\tilde{\mu}$, and then we iterate “infinitely many" $\hat{\mu}_{m+1}:=f_N(\hat{\mu}_m)$ as in corollary (ref), then with high probability $\hat{\mu}_m$ approaches a fixed point $\hat{\mu}^*$. As is standard in numerical optimization methods, “infinitely many" repetitions can in practice be read as “sufficiently many repetitions". None of the results so far in this paper provides any guidance regarding how many repetitions are enough. In fact, one of proposition (ref)'s corollaries can be concerning: equation ((ref)) establishes that if we initiate with a $\sqrt{N}$ consistent estimator, then we can only hope the iteration process to return $\sqrt{N}$-consistent estimators if we stop after a finite number of iterations. Therefore, from equation ((ref)), we get a sense of what a lower bound on the number of iterations should be, and it is rather massive. The number of iterations should be a diverging function of the sample size $N$ for us to have any hope to escape the $\sqrt{N}$ rate of convergence. How fast the number of iterations grows with $N$ will have an effect on the rate of convergence of the final estimator, but it is hard to tell what the proper order of magnitude is. It is even less clear what the rate of convergence would be if the number of iterations is indexed on some stoppage criterion on the value of the objective function, as is usually the case in standard numerical optimization algorithms. In simulations, the question of the number of iterations does not seem to be problematic. The standard optimization methods deliver distributions that are in line with the predictions of the asymptotic results presented so far, in particular the asymptotic distribution of $\hat{\mu}^*$ in corollary (ref).

Fortunately, the propositions (ref) and (ref) can be put to use differently to extract an estimator that is asymptotically equivalent to the minimizer $\hat{\mu}^*$. The alternative estimator requires exactly 2 iterations over the function $f_N$ and is therefore numerically more efficient. Using the alternative estimator, we can circumvent the concerns we highlighted around the number of iterations that are sufficient to achieve the desired asymptotic distribution.

An equivalent estimator

First, assume that the matrix $K$ is observed. Let $\Tilde{\mu}$ be an initial $\sqrt{N}$ consistent estimator and let $\hat{\mu}_1$ be the estimator returned in the equation ((ref)) after a single iteration. Write $\check{\mu}_1^*:=G\hat{\mu}_1+(I_L-G)\tilde{\mu}$, for some fixed $L\times L$ matrix $G$ and where $I_L$ is the identity matrix of dimension $L$. We will choose the matrix $G$ so that $\check{\mu}_1^*$ converges at rate $N$. Write

choosing $G$ such that $I_L - G(I_L-K)=0$ - i.e. $G=(I_L-K)^{-1}$ - yields a rate $N$ estimator. Note that by the proposition (ref), $I_L-K$ is invertible and $G$ is well defined.

In practice, the matrix $K$ is not observed. Instead, it needs to be estimated and plugged in to generate an estimator for $G$. Assume we have a consistent estimator $\hat{K}$ for $K$. Define $\hat{G}:=(I_L-\hat{K})^{-1}$ and $\check{\mu}_1:= \hat{G}\hat{\mu}+(I_L-\hat{G})\tilde{\mu}$. As for $\check{\mu}^*$

If $\hat{K}$ is a $\sqrt{N}$ - consistent estimator for $K$, that is, if $\hat{K}-K=O_p\left(\frac{1}{\sqrt{N}}\right)$, then the new estimator $\check{\mu}$ is rate optimal. The following proposition offers an example of a $\sqrt{N}$ consistent estimator for $K$.

Proposition (ref) allows for the construction of an estimator that is rate optimal. However, studying the asymptotic distribution of $\check{\mu}_1$ defined in ((ref)) is challenging. It requires that we determine the joint asymptotic distribution of $\hat{K}$, $\sqrt{N}(\tilde{\mu}-\mu_0)$ and the residual of order $O_p\left(\frac{1}{\sqrt{N}}\right)$ in equation ((ref)). However, as for the study of the fixed point $\hat{\mu}^*$, as we iterate, the effect of first stage estimator fades away. Rather than iterating here again, we use the same linear combination trick that allows us again to achieve the “infinite iterations" distribution using one iteration only.

Let $\hat{\mu}_2:=f_N(\check{\mu}_1)$ and define $\check{\mu}_2:= \hat{G}\hat{\mu}_2+(I_L-\hat{G})\check{\mu}_1$. Following the steps in equation ((ref)),

the last equality is a consequence of proposition (ref). This proves that $\check{\mu}_2$ is asymptotically equivalent to $\hat{\mu}^*$, the fixed point studied through corollary (ref).

To summerize, the alternative estimation procedure follows these steps:

Inference and bias correction

To be able to do inference on the (asymptotically equivalent) estimators presented in the previous section. We need to

{ $$\sigma_V^2 (I-K)^{-1 }\left(E(X_{12}X_{12}')- \frac{E(U_1)^2}{E(U_1^2)}E(X_{12}X_{32}')\right)^{-1} \Sigma \left(E(X_{12}X_{12}')- \frac{E(U_1)^2}{E(U_1^2)}E(X_{12}X_{32}')\right)^{-1}(I-K)^{-1 } $$} (c.f. proposition (ref))

To provide a consistent estimator for the bias term, we can use $\hat{K}$ defined in proposition (ref) as a consistent estimator for $K$. The matrices $ E(X_{12}X_{12}')$, $E(X_{12}X_{32}')$ and the vector $E(X_{12})$ can be estimated through their sample analogues. $\delta$ and the moments of $U$ are what is left to be estimated. Assume that $\delta=1$, section 2 explained how the eigenvector corresponding to the largest eigenvalue of $M(\tilde{\mu})$ is a good approximation to the normalized vector $\frac{U}{|| U||_2}$. Moreover, the largest eigenvalue informs about $U'U$, the norm of the vector $U$. Combining both, we can recover an estimator for $U$. When $\delta$ is -1, then we reason in terms of the largest eigenvalue in absolute value, and its corresponding eigenvalue. The difference when $\delta=-1$ is that the corresponding eigenvalue in fact estimates $-\frac{U}{||U||_2}$ rather than $\frac{U}{||U||_2}$ and a sign correction is necessary. The sign of $\delta$ is, with probability approaching 1, the sign of the largest eigenvalue in absolute value. These ideas are formalized through lemma (ref).

Lemma (ref) leads to the following consistent estimator for the bias term.

Finally, to obtain an estimator for the covariance matrix, it is left to provide a consistent estimator for the variance $\sigma_V^2$. Let's go back to the model ((ref)) $$ Y_{ij}=\sum_{l=1}^L \beta_{0,l} X_{ij,l}+ \delta U_i U_j+V_{ij} $$ Denote $\epsilon_{ij}:= \delta U_i U_j- \delta E(U_1)^2+V_{ij} $. First, observe that $$ \sigma_\epsilon^2:=Var(\epsilon_{ij})=E(U_1^2)^2-E(U_1)^4+\sigma_V^2 $$ Therefore, given the estimators for $E(U_1)$ and $E(U_1^2)$ provided in lemma (ref), any estimator $\hat{\sigma}_\epsilon^2$ that is consistent for $\sigma_\epsilon^2$, we obtain a consistent estimator $\hat{\sigma}_V^2:= \hat{\sigma}_\epsilon^2- \left( \frac{\hat{U}_i^2}{N}\right)^2+ \left( \frac{\hat{U}_i}{N}\right)^4$, where $\hat{U}_i$ are defined in lemma (ref). Notice that the sub-sample of the data $\{ (Y_{2k-1, 2k}, X_{2k-1, 2k} ), k=1...[N/2]\}$ is i.i.d. with i.i.d residuals $\epsilon_{2k-1, 2k}$.

Using the i.i.d. sub-sample, standard estimators for $\sigma_\epsilon^2$ are available from the cross-sectional OLS literature. Following for instance Wooldridge2010 (Section 4.2.2.), let $\hat{\epsilon}_{2k-1,2k}$ for $k=1...[N/2]$ be the residuals OLS fit of $Y_{2k-1, 2k}$ on $X_{2k-1, 2k}$ (and an intercept, when $X_{2k-1, 2k}$ does not include an intercept). A consistent estimator for $\epsilon$ is

To summarize:

Simulation study

I run $S=10000$ simulations on each of the 4 following designs, with a network of $N=100$ nodes in each simulation.

for each of the two designs $X\sim Unif(0,1)$, $\beta_{0,1}=\beta_{0,2}=E(A_1^2)=E(V_{12}^2)=1$.\footnote{I also generated simulations with $X\sim \mathcal{N}(0,1)$ or $X\sim 1+ \mathcal{N}(0,1)$ and the outcomes are similar.} The histograms for the estimated slope parameters $\beta_{0,2}$ are in the figures (ref) to (ref). In each graph, we show the histogram for the OLS estimator (in bleu) on the original model ((ref)) as a benchmark, the estimator $\hat{\mu}_{EIG}$ defined in this paper in green. The OLS estimator is semi-parametrically efficient in the model without individual effects, $Y_{ij}:=\beta_{0,1}+\beta_{0,2}X_{ij}+V_{ij}$ as a “gold standard" in orange, it is also estimated for each of the simulations and displayed in orange in the figures (ref) to (ref) as an oracle estimator. The estimators for the intercepts are not shown since the slope parameter are our concern in this paper. As discussed in the introduction, our estimator is $N-$ consistent for $\beta_{0,1}-\delta \gamma=1-1=0$ rather than for $\beta_{0,1}=1$. The term $\delta \gamma $ can't be estimated at a higher rate than $\sqrt{N}$. Any estimator for $\beta_{0,1}$ based on our estimator and an estimated correction for $\delta \gamma $ would only yield a $\sqrt{N}$-consistent estimator, even though $\beta_{0,1}-\delta \gamma$ is estimated at rate $N$.

The first two histograms (figures (ref) and (ref)) confirm the result in proposition (ref). The histogram for the eigenvalue-corrected estimator (in green) is close to the oracle (orange). On both histogram, the OLS estimator (blue) seems to have a larger variance. In fact, the OLS estimator has a non standard asymptotic distribution (cf. Menzel2021) and its distribution is slightly skewed to the left. The skew is not visible in figures (ref) and (ref), because the variance of $A$ is not large enough (see figure (ref) for a version of figure (ref) with a $Var(A)=100$ and where the skew is now obvious on the OLS estimator, whereas the eigenvalue corrected estimator is unaffected).

Figures (ref) and (ref) show that the histogram of the OLS estimator (blue) is much less concentrated than the eigenvalue-corrected estimator (green). This reflects the prediction of corollary (ref). The eigenvector corrected estimator is itself less efficient than the oracle (orange), but is rate optimal.

Extensions

The idea that this paper explores finds its origins in a simple observation: the matrix $ UU' + V$ has a single order $N$ eigenvalue that is due to the interaction between (re-centered) individual effects $U_i \times U_j$. A natural extension, when the model has higher order interactions, as in the model $$Y_{ij}=X_{ij}\mu_0+U_iU_j+U_i^2U_j^2+ V_{ij}$$ for instance. Then, similarly, the matrix of residuals would have two largest eigenvalues of order $N$ and $N-2$ eigenvalues of order at most $\sqrt{N}$. An estimator of the same flavor as the eigenvector corrected estimator discussed in this paper would consist of removing the two largest eigenvalues from the objective function in equation ((ref)), rather than just the largest eigenvalue.

More generally, if the model of interest is instead $$ Y_{ij}=X_{ij}\mu_0+g(U_i,U_j)+ V_{ij}$$ for some unknown but well-behaved and symmetric function $g$, under some smoothness assumption, $g$ has a spectral representation: $$g(x,y)=\sum_{k=1}^{+\infty}c_k \phi_k(x) \phi_k(y)$$ for some $\mathcal{L}^2$-orthonormal basis $(\phi_k)$ and a decaying sequence $c_k$ of real scalars. Because $c_k$ vanishes to 0, correcting for a certain number of first terms in the expansion of $g$ by removing the corresponding number of largest eigenvalues seems like an interesting extension of the estimator the current paper concerned itself with.

Proofs and intermediary results

This section details the proofs of all the results in the paper. It begins by showing how the OLS estimator of the intercept (in mode ((ref))) can be adjusted to obtain a $\sqrt{N}$- consistent estimator of the modified estimator (in model (ref)). Then we provide the technical ingredients (propositions (ref) and (ref)) that our main results heavily rely on.

Adjustment to the intercept

On the distribution of the largest eigenvalue

Proof of Lemma (ref)

Assume $X$ satisfies the lemma's assumptions. Let $\lambda$ be one of $X$'s eigenvalues and let $x$ be a corresponding eigenvector. Then

where $||.||_2$ designates the Euclidean norm for vectors and the spectral norm for matrices. Hence $$ |\lambda|\leq ||X||_2$$ but we know that the spectral normal is smaller than the Forbenius norm for any matrix. Therefore: $$ |\lambda| \leq \sqrt{\sum_{i,j} X_{ij}^2}$$ It is left to show that $\sum_{i,j} X_{ij}^2=O_p(N^2)$. Decompose $$\sum_{ij}X_{ij}^2=\sum_{ij} E(X_{ij}^2|X_i, X_j)+\sum_{ij}X_{ij}^2-E(X_{ij}^2|X_i, X_j)$$ First, by a U-statistic law of large numbers (e.g. theorem 3.1.3. in Ustattheory), $\sum_{ij} E(X_{ij}^2|X_i, X_j)= O_p(N^2)$. For the second term in the decomposition, it is enough to note: $$ Var\left(\frac{1}{N^2}\sum_{ij}X_{ij}^2 \bigg| (X_i)_{i=1}^\infty \right)\rightarrow 0, \; \mbox{ almost surely.}$$

Proof of lemma (ref)

Proof of propositions (ref) and (ref)

Proof of proposition (ref)

Proof of proposition (ref)

Proof of corollary (ref)

Proof of proposition (ref)

Proof of lemma (ref)