Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
158,146 characters · 16 sections · 80 citation commands
\noindentAbstract
One simple, and often very effective, way to attenuate the impact of nuisance parameters on maximum likelihood estimation of a parameter of interest is to recenter the profile score for that parameter. We apply this general principle to the quasi-maximum likelihood estimator (QMLE) of the autoregressive parameter $\lambda$ in a spatial autoregression. The resulting estimator for $\lambda$ has better finite sample properties compared to the QMLE for $\lambda$, especially in the presence of a large number of covariates. It can also solve the incidental parameter problem that arises, for example, in social interaction models with network fixed effects, or in spatial panel models with individual or time fixed effects. However, spatial autoregressions present specific challenges for this type of adjustment, because recentering the profile score may cause the adjusted estimate to be outside the usual parameter space for $\lambda$. Conditions for this to happen are given, and implications are discussed. For inference, we propose confidence intervals based on a Lugannani--Rice approximation to the distribution of the adjusted QMLE of $\lambda$. Based on our simulations, the coverage properties of these intervals are excellent even in models with a large number of covariates.
Among the several difficulties posed by nuisance parameters, a fundamental problem in a frequentist framework is that the profile likelihood function for a parameter of interest is typically not a genuine likelihood, in the sense that it does not correspond to the density of any observable random variable. To tackle this issue, a number of modified profile likelihoods have been proposed LaskarKing2001, Pace2006. Such modified profile likelihoods are genuine likelihoods only in special cases, but they can often be interpreted as approximations to a genuine likelihood. In practice, they tend to perform better than the profile likelihood, particularly when there is little information in the data about the nuisance parameters, which is likely to happen, for example, if the number of nuisance parameters is large relative to the sample size. The better performance of modified profile likelihoods is not necessarily captured by first-order asymptotic theory. Indeed, if the number of nuisance parameters does not depend on the sample size, modified profile likelihoods usually produce estimators that are first-order asymptotically equivalent to the maximum likelihood estimator (MLE). On the other hand, if the number of nuisance parameters increases with the sample size, a modified profile likelihood may be preferable even in terms of first-order asymptotic properties Jiang96,Sartori2003, Arellano2007.
One consequence of a profile likelihood not being a genuine likelihood is that the expectation of the profile score is generally nonzero, which means that setting the profile score to zero does not provide an unbiased estimating equation. A simple modified profile likelihood is therefore obtained by centering the profile score. We refer to the modified profile likelihood obtained in this way as the adjusted profile likelihood, and to the associated MLE as the adjusted MLE. The principle underlying this type of adjustment has been motivated from various perspectives Neyman1948, Conniffe1987, McCullagh1990, and has been applied to several statistical models Macaskill1993, DhaeneJoch2015.
The present paper is concerned with the adjusted MLE of the autoregressive parameter $\lambda$ in a spatial autoregression. Reliable estimation of $\lambda$ is important in many applications in economics, as well as in other fields. For example, in social interaction models $\lambda$ captures the endogenous effect, the assessment of which may be crucial for policy purposes Mansky93,Moffitt01,LeeLiuLin2010. More precisely, we will consider the quasi-MLE (QMLE), that is, the MLE obtained by maximizing the Gaussian likelihood, but without assuming that the error distribution is truly Gaussian. The literature on estimating (cross-sectional or panel) spatial autoregressions is now very large, the QMLE being perhaps the most popular estimator. An early reference for the QMLE in spatial autoregressions is Ord1975, whereas a rigorous first-order asymptotic analysis of the QMLE was given only much later, in an influential paper by Lee2004, under conditions that have become standard in the literature. Recently, higher-order approximations to the distribution of the QMLE of $\lambda$ have become available Robinson2015, and, motivated by the fact that the QMLE of $\lambda$ can suffer from substantial bias, a number of studies have suggested bias reduction procedures BaoUllah2007, Bao2013, Yang2015. In fact, part of the bias in the QMLE of a parameter of interest can be attributed to the bias in the profile likelihood estimating equation, so centering the profile score can itself be interpreted as a bias reduction technique. LeeYu2010, while discussing bias reduction in spatial panel data models, state explicitly that correction methods based on modifying the score \textquotedblleft might be possible and would be of interest in future research\textquotedblright\ (see their footnote 24). Yu2015 show that, in the absence of incidental parameters, the adjusted QMLE of $\lambda$ is first-order asymptotically equivalent to the QMLE of $\lambda$. Yu2015 also derive the second-order bias of the adjusted QMLE, and compare by simulation the adjusted QMLE with other bias reduction techniques. Prior to that article, Durban2000 had provided a preliminary investigation of the adjusted QMLE in a class of models that include spatial autoregressions. Notwithstanding the last two papers, the general distributional properties of the adjusted QMLE of $\lambda$ remain unclear. Indeed, Yu2015 conclude their analysis by suggesting that the adjusted QMLE of $\lambda$ \textquotedblleft deserves a deep study in future\textquotedblright.
Our main contributions are as follows. First, on studying the properties of the adjusted profile likelihood, we find that the distributions of the QMLE of $\lambda$ and its adjusted version may be supported on different intervals. This is due to the fact that, in spatial autoregressions, the parameter $\lambda$ is usually restricted to a certain interval containing the origin. Such a restriction is incorporated in the QMLE (which, strictly speaking, should therefore be referred to as a restricted QMLE), but may be violated once the profile score is recentered. We discuss the implications of the supports being different, and give conditions for this to happen. Second, for inference about $\lambda$, we propose confidence intervals based on the Lugannani80 saddlepoint approximation to the cdf of the adjusted QMLE. Contrary to the commonly employed Wald confidence intervals, these confidence intervals have excellent coverage properties even when the dimension of the nuisance parameter is large. Third, we consider the social interaction model of LeeLiuLin2010, and show that the adjusted QMLE solves the incidental parameter problem due to the network fixed effects without requiring the condition on $W$ (row normalization) that is needed for the estimator in LeeLiuLin2010. Fourth, we compare the QMLE and the adjusted QMLE by Monte Carlo simulation, and find evidence that the latter is preferable in a variety of circumstances.
The rest of the paper is structured as follows. Section (ref) introduces the spatial autoregressive (SAR) model and the QMLE of $\lambda$. Section (ref) discusses the properties of the adjusted profile likelihood for $\lambda$, and introduces the confidence intervals based on the adjusted QMLE. Section (ref) briefly considers the spatial error model, which is less popular than the SAR model in economic applications, but provides important motivation for the adjusted ML procedure. Section (ref) contains simulation evidence on the performance of the adjusted QMLE and associated confidence intervals. Section (ref) discusses extensions. Appendix (ref) contains some auxiliary results needed for the proofs, which can be found in Appendix (ref). Various technical materials related to the paper can be found in the online Supplement.
Throughout the paper, all vectors and matrices are real valued unless otherwise indicated. The null space of a matrix $A$ is denoted by $\operatorname{null}(A)$, the column space by $\operatorname{col}(A),$ and the orthogonal complement of $\operatorname{col}(A)$ by $\operatorname{col} ^{\perp}(A)$. Also, $M_{A}$ denotes the orthogonal projector onto $\operatorname{col}^{\perp}(A)$ ($M_{A}\coloneqq I_{n}-A(A^{\prime}A)^{-1}A^{\prime}$ if $A$ has full column rank). Finally, $\mu_{\mathbb{R}^{n}}$ denotes the Lebesgue measure on $\mathbb{R}^{n}$, and \textquotedblleft a.s.\textquotedblright\ stands for almost surely, with respect to $\mu_{\mathbb{R}^{n}}$.
We consider the spatial autoregressive (SAR) model
where $y$ is the $n\times1$ vector of observed random variables, $\lambda$ is a scalar parameter, $W$ is a spatial weights matrix, $X$ is an $n\times k$ matrix of regressors with full column rank and with $k\leq n-2$, $\beta \in\mathbb{R}^{k}$, $\sigma$ is a positive scale parameter, and $\varepsilon$ is a zero mean $n\times1$ random vector. For simplicity, we take $X$ and $W$ to be non-stochastic, and we assume that $W$ is completely known. Alternatively, we could allow $X$ and $W$ to be random, and interpret the analysis as conditional on them, provided that they are independent of $\varepsilon$. Some of the columns of $X$ may be spatial lags of some other columns, to allow for the estimation of, in the terminology of social network analysis, contextual effects. When there is no $X$, the model is called a pure SAR model. We allow for equation ((ref)) to also represent a spatial panel data model, or a model in which individuals are divided in several networks. In both cases, $W$ is a block diagonal matrix, with the number of blocks being given by, respectively, the number of time points and the number of networks. Additive fixed effects along one or both of the panel dimensions, or network fixed effects, can be added. For the purpose of estimation, fixed effects are treated as parameters, and hence can be included in $\beta$. Throughout the paper we assume that $W$ has at least one (real) negative eigenvalue and at least one (real) positive eigenvalue. This assumption is virtually always satisfied in applications, especially because the diagonal entries of $W$ are usually set to zero. The smallest real eigenvalue of $W$ is denoted by $\omega_{\min}$, and the largest real eigenvalue is normalized, without loss of generality, to 1.
On rewriting equation ((ref)) as $S(\lambda)y=X\beta+\sigma\varepsilon$, where $S(\lambda)\coloneqq I_{n}-\lambda W$, it is clear that in order for $y$ to be uniquely determined, given $X$\ and $\varepsilon$, it is necessary that $S(\lambda)$ is nonsingular. This requires $\lambda\neq\omega^{-1}$, for any nonzero real eigenvalue $\omega$ of $W$ (nonreal complex eigenvalues of $W$ do not need to be considered here, because $\lambda$ is assumed to be real, and $\omega^{-1}$ is real if and only $\omega$ is). In both applications and theoretical studies, the parameter space for $\lambda$ is usually restricted much further, namely to the largest interval containing the origin in which $S(\lambda)$ is nonsingular, that is, \[ \Lambda\coloneqq (\omega_{\min}^{-1},1), \] or a subset thereof (possibly independent of $n$) such as $(-1,1)$. Without such restrictions the models are believed to be too erratic to be useful in practice, and $\lambda$ is difficult to interpret.\footnote{From a large sample perspective, the restriction to $(-1,1)$ is often imposed, along with a uniform boundeness condition on $W$, to guarantee that the variances of the $y_{i}$'s do not explode as $n$ grows. Also, note that $\Lambda$ and the entries of $W$ are allowed to depend on $n$, although this is not emphasized in our notation.}
The following assumption is required to rule out some pathological combinations of $W$ and $X$.
To provide some intuition for Assumption (ref), we note that the condition $M_{X}(\omega I_{n}-W)=0$ is equivalent to $\mathrm{\operatorname{col}}(\omega I_{n}-W)\subseteq\operatorname{col}(X)$, and we distinguish two cases. First, if Assumption (ref) is violated for the eigenvalue $\omega=0$ (i.e., $\mathrm{\operatorname{col}} (W)\subseteq\operatorname{col}(X)$), it is evident from equation ((ref)) that estimating $\beta$ and $\lambda$ separately must be problematic. Second, if Assumption (ref) is violated for some eigenvalue $\omega\neq0$, then for any $y\in\mathbb{R}^{n}$ it is possible to find a $\beta\in \mathbb{R}^{k}$ such that $S(\omega^{-1})y=X\beta$, meaning that a SAR model with $\lambda=\omega^{-1}$ can provide perfect fit for any $y$. It is clear that in this case any sensible inferential procedure should suggest that $\lambda=\omega^{-1}$ and $\sigma=0$, for any $y$, and whatever the true values of $\lambda$ and $\sigma$ are.\footnote{For the specific case of the QMLE, see part (i) of Lemma (ref) in Section (ref) of the Supplement. Section (ref) of the Supplement also contains further technical remarks about Assumption (ref).}
The following example provides a simple illustration of Assumption (ref).
\begin{example2} [Group interaction]There are $R\geq1$ groups of individuals. Individuals interact uniformly within their group, and do not interact across different groups. If each group has the same size, say $m>1$, this type of interaction can be represented by the block-diagonal weights matrix
where $B_{m}\coloneqq \frac{1}{m-1}\mathopen{}\mathclose\bgroup\originalleft( \iota_{m}\iota_{m}^{\prime}-I_{m}\aftergroup\egroup\originalright) $, with $\iota_{m}$ an $m\times1$ vector of all ones. See Lee2007b and HillierMartellosio2016 for theoretical studies of this model, and Carrell2013 and boucher2014 for recent applications. For matrix ((ref)), one can easily verify that $\omega_{\min}=-\frac{1}{m-1}$ and $\operatorname{col}(\omega_{\min}I_{n}-W)=\operatorname{col}(I_{R}\otimes \iota_{m}).$ Noting that $I_{R}\otimes\iota_{m}$ is the design matrix of the group fixed effects, it follows that, in a SAR model with weights matrix ((ref)), Assumption (ref) is violated (for $\omega=\omega _{\min}$) whenever group fixed effects are included in the model (recall that the fixed effects are treated as parameters to be estimated and hence are included in $\beta$). It should be noted that the presence of group intercepts does not cause a violation of Assumption (ref) when the model is unbalanced, that is, not all groups have the same size. \end{example2}
Indeed, it is well known that the model of Example (ref) suffers from an identifiability problem. Specifically, Lee2007b and Bramoulle2009 show that the parameters of a SAR model with weights matrix ((ref)), group fixed effects, and contextual effects are not identifiable after removal of the group fixed effects by a within transformation. It is easily verified that, in this model, the identifiability problem occurs even without contextual effects.
We now define the QMLE of the parameters in model ((ref)). By quasi-likelihood we mean the likelihood that would prevail under the condition $\varepsilon\sim\mathrm{N}(0,I_{n})$. Omitting additive constants, the quasi-log-likelihood is
for any $\lambda$ such that $S(\lambda)$ is nonsingular. To avoid tedious repetitions, in the remainder of the paper we will often omit the \textquotedblleft quasi-\textquotedblright\ in front of \textquotedblleft log-likelihood\textquotedblright.
The QMLE in most common use is the maximizer of $l(\beta,\sigma^{2},\lambda)$ under the condition that $\lambda$ is in $\Lambda$ (or in a subset thereof). That is, the QMLE is \[ (\hat{\beta}_{\mathrm{ML}},\hat{\sigma}_{\mathrm{ML}}^{2},\hat{\lambda }_{\mathrm{ML}})=\underset{\beta\in\mathbb{R}^{k},\hspace{0.1667em}\sigma ^{2}>0,\hspace{0.1667em}\lambda\in\Lambda}{\operatorname*{argmax}}l(\beta,\sigma^{2} ,\lambda). \] Maximization with respect to $\beta$ and $\sigma^{2}$ gives $\hat{\beta }_{\mathrm{ML}}(\lambda)\coloneqq (X^{\prime}X)^{-1}X^{\prime}S(\lambda)y$ and $\hat{\sigma}_{\mathrm{ML}}^{2}(\lambda)\coloneqq \frac{1}{n}y^{\prime}S^{\prime }(\lambda)M_{X}S(\lambda)y$. The corresponding profile, or concentrated, log-likelihood for $\lambda$ is, again omitting additive constants,
The QMLE of $\lambda$ can be equivalently defined as
The function $l(\lambda)$ is a.s.\ well defined (see Section (ref) in the Supplement), and, clearly, is continuously differentiable whenever it is well defined. Also, it is easy to see that, under Assumption (ref) , $l(\lambda)$ a.s.\ goes to $-\infty$ at each real zero of $\det(S(\lambda))$ Hillier2017. Thus, $l(\lambda)$ has a.s.\ at least one critical point corresponding to a maximum in any interval between two consecutive real zeros of $\det\mathopen{}\mathclose\bgroup\originalleft( S(\lambda)\aftergroup\egroup\originalright) $. We also define the unrestricted QMLE $\hat{\lambda}_{\mathrm{uML}}\coloneqq \operatorname*{argmax}_{\lambda\in\Lambda_{u}}l(\lambda )$, where $\Lambda_{u}\coloneqq \{\lambda\in\mathbb{R}:\det(S(\lambda))\neq0\}$. The estimators $\hat{\lambda}_{\mathrm{ML}}$ and $\hat{\lambda}_{\mathrm{uML}}$ are different, because the global maximum of $l(\lambda)$ over $\Lambda_{u}$ is not necessarily in $\Lambda$, even when the true value of $\lambda$ is in $\Lambda$.
A profile likelihood for a parameter of interest does not take into account the sampling variability associated to the estimation of the nuisance parameters, and hence, as mentioned in the introduction, is generally not a genuine likelihood. As a consequence, a profile score estimating equation is generally biased, which is likely to induce bias in the QMLE of the parameter of interest. The adjusted QMLE solves the estimating equation obtained by recentering the profile score. We now apply this general adjustment to the estimation of $(\sigma^{2},\lambda)$ in the SAR model. The score for $(\sigma^{2},\lambda)$ is centered assuming only that $\mathrm{E} (\varepsilon)=0$ and $\mathrm{var}(\varepsilon)=I_{n}$. Note that we treat only $\beta$ as the nuisance parameter, not $(\beta,\sigma^{2})$, because an adjusted estimator of $\sigma^{2}$ is required for estimation of $\lambda $.\footnote{Treating $\sigma^{2}$ as a nuisance parameter too, and consequently recentering the score for $\lambda$ only (i.e., the score associated to the log-likelihood ((ref))) would not produce an adjusted estimator for $\sigma^{2}$. Also, the (exact) recentering of $s(\lambda)$ would require stronger assumptions than what required for the recentering of $s(\sigma^{2},\lambda)$; see Section (ref) of the Supplement.}
On concentrating just $\beta$ out of the Gaussian log-likelihood ((ref)), the profile log-likelihood for $(\sigma^{2},\lambda)$ is
with profile score
where $G(\lambda)\coloneqq WS^{-1}(\lambda)$. We now compute the expectation of $s(\sigma^{2},\lambda)$ under the SAR model $y=\lambda Wy+X\beta +\sigma\varepsilon$. Assuming that $\mathrm{E}(\varepsilon)=0$ and $\mathrm{var}(\varepsilon)=I_{n}$, we have\footnote{If the matrices $X$ or $W$ were stochastic, but independent of $\varepsilon$, we would be conditioning on them here.}
For a pure model, $\mathrm{E}(s(\sigma^{2},\lambda))=0$ for any $\lambda$ such that $S(\lambda)$ is nonsingular, meaning that the estimating equation $s(\sigma^{2},\lambda)=0$ is unbiased.\footnote{Depending on the sample size and on the structure of $W$, $\hat{\lambda}_{\mathrm{ML}}$ can be considerably biased even in the pure case, despite the fact that the estimating equation $s(\sigma^{2},\lambda)=0$ is unbiased in that case. The adjustment studied in this paper is not designed to reduce this type of bias in $\hat{\lambda }_{\mathrm{ML}}$. Instead, it aims at reducing the bias due to the nuisance parameter $\beta$.} When regressors are present, however, the unaccounted variability in the estimation of $\beta$ causes the estimating equation $s(\sigma^{2},\lambda)=0$ to be biased. Note that the expectation ((ref)) does not depend on $\beta$, so recentering the score is straightforward. The adjusted profile score for $(\sigma^{2},\lambda)$, defined as $s_{\mathrm{a} }(\sigma^{2},\lambda)\coloneqq s(\sigma^{2},\lambda)-\mathrm{E}(s(\sigma^{2} ,\lambda))$, is
Setting $s_{\mathrm{a}1}(\sigma^{2},\lambda)=0$ gives $\hat{\sigma }_{\mathrm{aML}}^{2}(\lambda)\coloneqq \frac{n}{n-k}\hat{\sigma}_{\mathrm{ML}} ^{2}(\lambda)$. That is, recentering the score automatically delivers the usual degrees of freedom correction for $\hat{\sigma}_{\mathrm{ML}} ^{2}(\lambda)$. The adjusted QMLE for $\lambda$ must then be a zero of
or, which is a.s.\ the same, must solve the estimating equation
where \[ R(\lambda)\coloneqq M_{X}\mathopen{}\mathclose\bgroup\originalleft( G(\lambda)-\frac{\mathrm{tr}(M_{X}G(\lambda))} {n-k}I_{n}\aftergroup\egroup\originalright) . \]
We emphasize that, by construction, the estimating equation ((ref)) is exactly unbiased provided that $\mathrm{E}(\varepsilon)=0$ and $\mathrm{var}(\varepsilon)=I_{n}$ (no further distributional assumptions are required).
Given the adjusted profile score $s_{\mathrm{a}}(\sigma^{2},\lambda)$, one can define the function with gradient equal to $s_{\mathrm{a}}(\sigma^{2} ,\lambda)$, which we refer to as the adjusted likelihood for $(\sigma ^{2},\lambda)$, denoted by $l_{\mathrm{a}}(\sigma^{2},\lambda)$. Letting $\operatorname{Re}\mathopen{}\mathclose\bgroup\originalleft[ \cdot\aftergroup\egroup\originalright] $ denote the real part of a complex number, the following result gives a closed form expression for $l_{\mathrm{a} }(\sigma^{2},\lambda)$.
From expression ((ref)), we immediately obtain the adjusted likelihood for $\lambda$ only,
(note that this is the likelihood with score $s_{\mathrm{a}2}(\lambda)$). Equations ((ref)) and ((ref)) may be useful for graphical or optimization purposes, but are not used for the analytical results that follow in the paper. Indeed, even the results that (for simplicity and to aid intuition) are stated in terms of the adjusted likelihood, are proved using only the expression for the adjusted score, and could equally be formulated in terms of the adjusted score.
So far, the adjusted QMLE of $\lambda$ has been introduced as a zero of $s_{\mathrm{a}2}(\lambda)$, but of course $s_{\mathrm{a}2}(\lambda)$ may have many (real) zeros. The question therefore arises as to how exactly the adjusted QMLE should be defined. Recall from Section (ref) that $\hat{\lambda}_{\mathrm{ML}}$ is defined as the maximizer of $l(\lambda)$ over $\Lambda$. One may therefore be tempted to define the adjusted QMLE of $\lambda$ as the maximizer of $l_{\mathrm{a}}(\lambda)$ over $\Lambda$. However, we shall see that $s_{\mathrm{a}2}(\lambda)$ may have no zeros in $\Lambda$ (or, equivalently, the adjusted profile likelihood $l_{\mathrm{a} }(\lambda)$ may have no maximum on $\Lambda$), so it makes sense to define the adjusted QMLE on an interval larger than $\Lambda$, $\Lambda_{\mathrm{a}}$ say. By analogy with the unadjusted case, we shall define $\Lambda _{\mathrm{a}}$ to be the shortest open interval containing the origin with the property that $l_{\mathrm{a}}(\lambda)\rightarrow-\infty$ a.s.\ at both extremes of $\Lambda_{\mathrm{a}}.$ But first, in order to fully understand the problem, we need to study the behavior of the profile score $s_{\mathrm{a} 2}(\lambda)$, or, equivalently, of $l_{\mathrm{a}}(\lambda)$, near the zeros of $\det(S(\lambda))$.
Perhaps unexpectedly, the profile score $s(\sigma^{2},\lambda)$ and its adjusted version $s_{\mathrm{a}}(\sigma^{2},\lambda)$ may have very different behavior as $\lambda$ approaches a real zero of $\det(S(\lambda))$. Obviously, different behavior of $s(\sigma^{2},\lambda)$ and $s_{\mathrm{a}}(\sigma ^{2},\lambda)$ implies different behavior of the functions $l(\sigma ^{2},\lambda)$ and $l_{\mathrm{a}}(\sigma^{2},\lambda)$, and hence of their profile versions $l(\lambda)$ and $l_{\mathrm{a}}(\lambda)$. For simplicity, we state the results in this section in terms of the univariate functions $l(\lambda)$ and $l_{\mathrm{a}}(\lambda)$. We shall see that the different behavior of $l(\lambda)$ and $l_{\mathrm{a}}(\lambda)$ accounts for important differences in the properties of the associated estimators for $\lambda$.
To start with, note that the analysis is straightforward for the pure model. In that case, we trivially have $l_{\mathrm{a}}(\lambda)=l(\lambda)$, and plugging $\hat{\sigma}_{\mathrm{ML}}^{2}(\lambda)=\frac{1}{n}y^{\prime }S^{\prime}(\lambda)S(\lambda)y$ in equation ((ref)) reveals that $l(\lambda)$ a.s.\ approaches $-\infty$ near any real zero of $\det (S(\lambda))$. The presence of regressors complicates the analysis. We shall confine attention to semisimple eigenvalues of $W$. An eigenvalue is said to be semisimple if its algebraic and geometric multiplicities are equal Meyer2000. While it simplifies the analysis considerably, the restriction to semisimple eigenvalues does not imply a great loss of generality. For example, all eigenvalues of a diagonalizable matrix are semisimple, and any simple eigenvalue is semisimple (an eigenvalue being simple if it has algebraic multiplicity equal to one). The behavior of $l_{\mathrm{a}} (\lambda)$ close to the eigenvalue $\omega=1$ is often particularly important, and that eigenvalue is in most cases semisimple in applications. For example, it is semisimple when $W$ is row stochastic (even if $W$ is not diagonalizable, and the algebraic multiplicity of $\omega=1$ is larger than one), or when $W$ is irreducible .\footnote{A row stochastic matrix is a square nonnegative matrix whose row sums are all equal to 1. Irreducibility can be defined in terms of the graph of a matrix as follows. Let the graph of an $n\times n$ matrix $A$ be the directed graph on $n$ vertices in which there is an an edge from vertex $i$ to vertex $j$ if and only $A(i,j)\neq0$. Also, call a graph strongly connected if there is a sequence of directed edges from any vertex $i$ to any vertex $j$. Then, $A$ is irreducible if and only if the graph of $A$ is strongly connected Meyer2000.}
The simplification afforded by the restriction to semisimple eigenvalues is that we can express the conditions for $l_{\mathrm{a}}(\lambda)$ to diverge or be bounded near a real zero of $\det(S(\lambda))$ in terms of a projector onto an eigenspace of $W$. Without the restriction to semisimple eigenvalues, the conditions would need to be stated in terms of projections onto generalized eigenspaces. If the eigenvalue $\omega$ of $W$ is semisimple, then $\mathrm{null}\mathopen{}\mathclose\bgroup\originalleft( W-\omega I_{n}\aftergroup\egroup\originalright) $ (the eigenspace of $W$ associated to the eigenvalue $\omega$) and $\operatorname{col}\mathopen{}\mathclose\bgroup\originalleft( W-\omega I_{n}\aftergroup\egroup\originalright) $ are complementary subspaces of $\mathbb{C}^{n}$, and therefore we can define a (unique) projector, denoted by $Q_{\omega}$, onto $\mathrm{null}\mathopen{}\mathclose\bgroup\originalleft( W-\omega I_{n}\aftergroup\egroup\originalright) $ along $\operatorname{col} \mathopen{}\mathclose\bgroup\originalleft( W-\omega I_{n}\aftergroup\egroup\originalright) $.\footnote{The matrix $Q_{\omega}$ is a \textquotedblleft spectral projector\textquotedblright\ of $W$; see, for instance, Meyer2000, where explicit general representations for such projectors can also be found. Particular cases are given in the proof of Lemma (ref) and in the proof of Lemma (ref) in the Supplement.} Recalling that the zeros of $\det(S(\lambda))$ are the reciprocals of the nonzero eigenvalues of $W$, we can prove the following result.
We can now compare the behavior of the adjusted profile log-likelihood $l_{\mathrm{a}}(\lambda)$ near the real zeros of of $\det(S(\lambda))$, as established by Theorem (ref), with the behavior of the unadjusted profile log-likelihood $l(\lambda)$ near those points. Recall from Section (ref) that $\lim_{\lambda\rightarrow\omega^{-1}} l(\lambda)=-\infty$ a.s., for any nonzero real eigenvalue $\omega$ of $W$ (under Assumption (ref)). Thus, $l(\lambda)$ and its adjusted version $l_{\mathrm{a}}(\lambda)$ have the same behavior near a point $\omega^{-1}$, for a semisimple nonzero real eigenvalue $\omega$, only in case (i) of Theorem (ref). In case (ii), $l_{\mathrm{a}}(\lambda)$ can be extended to a function that is a.s.\ continuous at $\lambda=\omega^{-1}$. This can be achieved by extending the domain of $l_{\mathrm{a}}(\lambda)$ to include $\omega^{-1}$, and setting $l_{\mathrm{a}}(\omega^{-1})\coloneqq \lim _{\lambda\rightarrow\omega^{-1}}l_{\mathrm{a}}(\lambda)$. From now on, when we say that $l_{\mathrm{a}}(\lambda)$ is continuous at $\lambda=\omega^{-1}$ we implicitly assume that this extension has been performed. In case (iii) of Theorem (ref), $l_{\mathrm{a}}(\lambda)$ is unbounded from above near $\omega^{-1}$. From unreported numerical experiments, it appears that $\mathrm{tr}(M_{X}Q_{\omega})<0$ is an extremely rare occurrence for pairs $(W,X)$ which are likely to be encountered in applications.\footnote{For example, for all simulation designs in Section (ref), we find that $l_{\mathrm{a}}(\lambda)$ is always bounded from above on $\Lambda_{\mathrm{a}}$. A complete understanding of when $l_{\mathrm{a}}(\lambda)$ may be unbounded from above would be of interest, but is left for future research. Note that the fact that $l_{\mathrm{a}}(\lambda)$ can be unbounded from above is not surprising, since $l(\lambda)$ itself can be unbounded from above (see Lemma (ref) in the Supplement).} We shall also see shortly that $\mathrm{tr}(M_{X}Q_{\omega})<0$ cannot occur if $W$ is symmetric.
Ruling out the pathological case (iii), it is useful to try and understand which of cases (i) and (ii) in Theorem (ref) is likely to occur in applications. At first sight, the condition $\mathrm{tr} (M_{X}Q_{\omega})=0$ in case (ii) may look very restrictive. Indeed, Lemma (ref) in Appendix (ref) establishes that $\mathrm{tr}(M_{X}Q_{\omega})\neq0$ for generic, in the measure theoretic sense, full column rank $X$ (and for any fixed $W$). However, $X$ typically contains an intercept (or group intercepts), and this implies that $\mathrm{tr}(M_{X}Q_{\omega})=0$ occurs, at $\omega=1$, for a very large class of matrices $W$ used in practice. To see why this is the case, the next lemma provides a condition for $\mathrm{tr}(M_{X}Q_{\omega})=0$ in terms of the eigenspace $\mathrm{null}(W-\omega I_{n})$ (the eigenspace of $W$ associated to the eigenvalue $\omega$).
Part (i) of Lemma (ref) establishes that $\mathrm{null}(W-\omega I_{n})\subseteq\operatorname{col}(X)$ is sufficient for case (ii) of Theorem (ref) to apply, and hence for $l(\lambda)$ and $l_{\mathrm{a} }(\lambda)$ to have different behaviour near $\lambda=\omega^{-1}$ (for any semisimple nonzero real eigenvalue $\omega$ of $W$, and provided that Assumption (ref) holds). It turns out that this sufficient condition is very often satisfied in applications. Two examples are given next, the second one being a generalization of the first.
\begin{example2} [Row stochastic and irreducible weights matrix]In applications of spatial autoregressions, $W$ is often row stochastic and irreducible (cf. footnote (ref)), and an intercept is included in the regression. By the Perron--Frobenius Theorem Horn1985, $\omega=1$ is a simple (and hence semisimple) eigenvalue of $W$, and the associated eigenspace $\mathrm{null} (W-I_{n})$ is spanned by a vector of identical entries, and therefore is in $\operatorname{col}(X)$. It follows, under Assumption (ref), that $l(\lambda)$ a.s.\ approaches $-\infty$ as $\lambda\rightarrow1$, while $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous at $\lambda=1$. \end{example2}
\begin{example2} [Block diagonal weights matrix]Example (ref) generalizes immediately to the case when $W$ is a block diagonal matrix whose blocks are row stochastic and irreducible matrices, of, say, size $m_{r}\times m_{r}$. That is, using direct sum notation, $W=\bigoplus_{r=1}^{R}W_{r}$. This situation arises, for instance, in a social interaction model on $R$ networks LeeLiuLin2010, or in a spatial panel model where individuals are followed over time LeeYu2010. When $W$ has this structure, the eigenspace $\mathrm{null}(W-I_{n})$ is spanned by the columns of the (network or time) fixed effects matrix $\bigoplus_{r=1} ^{R}\iota_{m_{r}}$, and therefore is in $\operatorname{col}(X)$ as long as the regressions contains those fixed effects. In that case, and provided that Assumption (ref) holds, $l(\lambda)$ a.s.\ approaches $-\infty$ as $\lambda\rightarrow1$, while $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous at $\lambda=1$. \end{example2}
According to part (ii) of Lemma (ref), when $W$ is symmetric, the condition $\mathrm{null}(W-\omega I_{n})\subseteq\operatorname{col}(X)$ is also necessary for $l_{\mathrm{a}}(\lambda)$ to be a.s.\ continuous at $\lambda=\omega^{-1}$. Thus, for symmetric $W$, Theorem (ref) reduces to the simple statement that $\lim_{\lambda\rightarrow\omega^{-1} }l_{\mathrm{a}}(\lambda)$ is a.s.\ bounded if $\mathrm{null}(W-\omega I_{n})\subseteq\operatorname{col}(X)$, $-\infty$ otherwise, for any semisimple nonzero real eigenvalue $\omega$.
We end this section by providing a graphical comparison of $l(\lambda)$ and $l_{\mathrm{a}}(\lambda)$. Consider a SAR model with weights matrix $W$ equal to the row normalized adjacency matrix of an Erd{\H{o}}s-R{\'{e}}nyi $G(n,p)$ graph Erdos1959. The $G(n,p)$ graph is a random graph on $n$ vertices, with an edge between any two vertices being present with probability $p$, independently of every other edge. Suppose that the regression contains an intercept, and that, for simplicity, the graph is connected. Then, since $W$ is row-stochastic and irreducible, $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous at $\lambda=1$ and a.s.\ approaches $-\infty$ as $\lambda$ approaches any other singularity of $S(\lambda)$. Figure (ref) displays $l(\lambda)$ and $l_{\mathrm{a}}(\lambda)$, for one random draw of $G(n,p)$ and for one random draw from the intercept-only model $y=.5Wy+\iota _{n}+\varepsilon$, with $\varepsilon\sim\mathrm{N}(0,I_{n})$. The log-likelihood functions are plotted for $\lambda\in(0,\omega_{3}^{-1})$, where $\omega_{3}$ is the third largest eigenvalue of $W$, and $l_{\mathrm{a} }(\lambda)$ has been lowered so that the two likelihoods have the same maximum value.\footnote{For the particular random draw of $\varepsilon$ underlying Figure (ref), $\hat{\lambda}_{\mathrm{ML}}$ and its adjusted version $\hat{\lambda}_{\mathrm{a}\mathrm{ML}}$ (defined as in Section (ref) below) are about 0.478 and 0.506, respectively. Over $10^{6}$ draws from $\varepsilon\sim\mathrm{N}(0,I_{n})$ (and for the same draw of $G(n,p)$ used for Figure (ref)), empirical bias and RMSE are $-0.039$ and $0.128$ for $\hat{\lambda}_{\mathrm{ML}}$ and $-0.007$ and $0.126$ for $\hat{\lambda}_{\mathrm{a}\mathrm{ML}}$.} Note that, contrary to $l(\lambda)$, $l_{\mathrm{a}}(\lambda)$ does not go to $-\infty$ as $\lambda\rightarrow1$.
Theorem (ref) establishes that the adjusted profile log-likelihood $l_{\mathrm{a}}(\lambda)$ may, in contrast to $l(\lambda)$, be a.s.\ continuous at the extremes of the parameter space $\Lambda$. As a consequence, there is no guarantee that $l_{\mathrm{a}}(\lambda)$ has a maximum over $\Lambda$, which suggests that $l_{\mathrm{a}}(\lambda)$ should be maximized over a larger set, $\Lambda_{\mathrm{a}}$ say. As anticipated in Section (ref), it is natural to define $\Lambda _{\mathrm{a}}$ as the shortest open interval containing the origin with the property that $l_{\mathrm{a}}(\lambda)\rightarrow-\infty$ a.s.\ at both extremes of $\Lambda_{\mathrm{a}}$. For example, in the case of Figure (ref), $\Lambda=(-1.195,1)$ and $\Lambda_{\mathrm{a}} =(-1.195,1.178).$ Note that the extremes of $\Lambda_{\mathrm{a}}$ must always be zeros of $\det(S(\lambda))$, because $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous between consecutive real zeros of $\det(S(\lambda ))$.\footnote{The fact that $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous between consecutive real zeros of $\det(S(\lambda))$ also means that, ignoring the pathological cases in which $l_{\mathrm{a}}(\lambda)$ is a.s.\ unbounded from above, $\Lambda_{\mathrm{a}}$ is the smallest open set containing the origin on which $l_{\mathrm{a}}(\lambda)$ is guaranteed to have a maximum a.s.} Hence, our definition of $\Lambda_{\mathrm{a}}$ requires the following assumption.
It is clear from Section (ref) that Assumption (ref) can be violated only in very special cases. In those cases, one could take the left (resp., right) endpoint of $\Lambda _{\mathrm{a}}$ to be $-\infty$ (resp., $+\infty$), but we refrain from doing this, for simplicity.
It is natural to define the adjusted QMLE of $\lambda$ as the maximizer of $l_{\mathrm{a}}(\lambda)$ over $\Lambda_{\mathrm{a}}$, that is, \[ \hat{\lambda}_{\mathrm{aML}}\coloneqq \operatorname*{argmax}_{\lambda\in\Lambda_{\mathrm{a}} }l_{\mathrm{a}}(\lambda). \]
Maximization of $l_{\mathrm{a}}(\lambda)$ over a subset of $\Lambda _{\mathrm{a}}$ would yield a censored version of $\hat{\lambda}_{\mathrm{aML} }$. In particular, let $\bar{\Lambda}$ be $\Lambda$ augmented with one of its endpoints if $l_{\mathrm{a}}(\lambda)$ is bounded near that endpoint (augmented with both endpoints if $l_{\mathrm{a}}(\lambda)$ is bounded near both endpoints), and let $\bar{\lambda}_{\mathrm{aML}}\coloneqq \operatorname*{argmax}_{\lambda \in\bar{\Lambda}}l_{\mathrm{a}}(\lambda)$ be the estimator that is obtained by maximizing $l_{\mathrm{a}}(\lambda)$ over $\bar{\Lambda}$. Since $\Lambda\subseteq\Lambda_{\mathrm{a}}$, $\bar{\lambda}_{\mathrm{aML}}$ is a censored version of $\hat{\lambda}_{\mathrm{aML}}$.
Note that the set $\Lambda_{\mathrm{a}}$, contrary to $\Lambda$, may depend on $X$ (because, by Theorem (ref), whether or not $l_{\mathrm{a} }(\lambda)\rightarrow-\infty$ a.s.\ at some zero of $\det(S(\lambda))$ depends on $X$). For a fixed $X$, $\Lambda_{\mathrm{a}}=\Lambda$ if $l_{\mathrm{a}}(\lambda)\rightarrow-\infty$ a.s.\ as $\lambda$ approaches the extremes of $\Lambda$, and $\Lambda_{\mathrm{a}}\supset\Lambda$ otherwise. But it is also possible to compare $\Lambda$ and $\Lambda_{\mathrm{a}}$ for generic $X$ (in the measure theoretic sense), or better, since the model typically contains an intercept, for a generic matrix $X$ containing an intercept. The following example is important.
\begin{example2} Suppose $W$ is row-stochastic and irreducible, and that an intercept is included in the model. Let $X=(\iota_{n},\widetilde{X})$ and $\widetilde{\mathcal{X}}\coloneqq \{\widetilde{X}\in\mathbb{R}^{n\times(k-1)} :\mathrm{rank}(X)=k\}$. We know from Example (ref) that, in this case, $l_{\mathrm{a}}(\lambda)$ is a.s.\ continuous at $\lambda=1$, for all $\widetilde{X}\in\widetilde{\mathcal{X}}$. Consider now an arbitrary semisimple nonzero real eigenvalue $\omega\neq1$ of $W$. By Lemma (ref) in Appendix (ref) and Theorem (ref), $l_{\mathrm{a}}(\lambda)$ is a.s.\ unbounded from below or from above near $\omega^{-1}$ for $\mu_{\mathbb{R}^{n\times(k-1)}}$-almost every $\widetilde{X}\in\widetilde{\mathcal{X}}$. Assuming that $\omega_{\min}$ and the second largest (positive) eigenvalue of $W$, denoted by $\omega_{2}$, are semisimple, it follows that $\Lambda_{\mathrm{a}}=(\omega_{\min} ^{-1},\omega_{2}^{-1})$, for $\mu_{\mathbb{R}^{n\times(k-1)}}$-almost every $\widetilde{X}\in\widetilde{\mathcal{X}}\setminus\mathopen{}\mathclose\bgroup\originalleft\{ \mathcal{P} _{\omega_{\min}} \cup\mathcal{P}_{\omega_{2}}\aftergroup\egroup\originalright\} ,$ where $\mathcal{P} _{\omega}$ is the set of pathological $\widetilde{X}$ such that $l_{\mathrm{a} }(\lambda)$ is a.s.\ unbounded from above near $\omega^{-1}$ (i.e., $\mathcal{P}_{\omega}\coloneqq \{\widetilde{X}\in\mathbb{R}^{n\times(k-1)} :\mathrm{rank}(X)=k$ and $\mathrm{tr}(M_{X}Q_{\omega})<0\}$). As discussed earlier, $\mathcal{P} _{\omega}$ is expected to be very small or even $\mu_{\mathbb{R} ^{n\times(k-1)}}$-null in cases of interest in applications. By Lemma (ref), $\mathcal{P}_{\omega}$ is empty if $W$ is symmetric, so in that case $\Lambda_{\mathrm{a}}=(\omega_{\min}^{-1} ,\omega_{2}^{-1})$ for $\mu_{\mathbb{R}^{n\times(k-1)}}$-almost every $\widetilde{X}\in\widetilde{\mathcal{X}}$. \end{example2}
Three remarks about the result in Example (ref) that $\Lambda_{\mathrm{a}}=(\omega_{\min}^{-1},\omega_{2}^{-1})$, for generic non-pathological $\widetilde{X}$ are in order. First, note that $\omega _{2}^{-1}>1$ and recall that $\Lambda=(\omega_{\min}^{-1},1)$, independently of $X$. Thus, how much larger $\Lambda_{\mathrm{a}}$ is compared to $\Lambda$ depends only on the eigenvalue gap $1-\omega_{2}$. There is considerable evidence in the graph theory literature that the eigenvalue gap $1-\omega_{2}$ tends to be large when the graph underlying $W$ has good connectivity and randomness properties Brouwer11, and this is indeed something we will come upon in our Monte Carlo experiments later. Second, on replacing the intercept $\iota_{n}$ with group intercepts $\bigoplus_{r=1}^{R}\iota_{m_{r}}$, as in Example (ref), the result in Example (ref) generalizes immediately to the case when $W$ is block diagonal with row stochastic and irreducible blocks. Third, the argument used in Example (ref) can also be applied to compare $\Lambda_{\mathrm{a}}$ and $\Lambda$ for weights matrices that are not row-stochastic and irreducible (or are not block diagonal with row-stochastic and irreducible blocks). To do this, note that the result $\Lambda _{\mathrm{a}}=(\omega_{\min}^{-1},\omega_{2}^{-1})$ for generic non-pathological $\widetilde{X}$ in Example (ref) arises because, in that case, $\iota_{n}\in\operatorname{col}(X)$ and $\iota_{n}$ spans the eigenspace $\mathrm{null}(W-I_{n})$, so that the condition in Part (i) of Lemma (ref) is satisfied. When $W$ is not row-stochastic and irreducible, such a special interaction between $\operatorname{col}(X)$ and $W$ will typically not occur, and as a result $\Lambda_{\mathrm{a}} =\Lambda=(\omega_{\min}^{-1},1)$ for generic non-pathological $\widetilde{X}$. This is most easily seen when $W$ is symmetric. In that case, by part (ii) of Lemma (ref) and Theorem (ref), $\Lambda _{\mathrm{a}}$ is different from $\Lambda$ only if $\mathrm{null}(W-I_{n})$ or $\mathrm{null}(W-\omega_{\min}I_{n})$ are in $\operatorname{col}(X)$. In general there is no reason why $\operatorname{col}(X)$ should contain those eigenspaces.
The original motivation for seeking an unbiased estimating equation is, of course, bias reduction, or more generally, improved inference on $\lambda$ (and possibly $\sigma^{2}).$ The simulation evidence reported in Section (ref) below is unequivocal that the hoped-for bias reduction is certainly achieved, and without detriment to the mean squared error. However, there is one caveat that must be mentioned here, which we discuss next.
We have seen that the intervals $\Lambda\ $and $\Lambda_{\mathrm{a}}$ on which the distributions of $\hat{\lambda}_{\mathrm{ML}}$ and $\hat{\lambda }_{\mathrm{aML}}$ are supported are different in many cases of interest. That is, there are circumstances in which $\hat{\lambda}_{\mathrm{aML}}$ lies outside $\Lambda.$ This raises the question of whether, from the point of view of interpreting the parameter $\lambda$, that outcome is acceptable.\footnote{The situation is somewhat similar to, but more subtle than, that in which the MLE for a variance (or covariance matrix) satisfies the expected non-negativity (or positive definiteness) requirement, but a bias-corrected version of it does not. For example, subtraction of the estimated first (bias) term in an asymptotic expansion for the expectation of the MLE can have this undesirable outcome. We also note that several other estimators of $\lambda$, for example IV, GMM, or indirect inference estimators, may have support larger than $\Lambda$ Kelejian98,Lee2007a,Kyriacou2017.} Ultimately, this will depend on the context, but censoring the estimator to ensure that it lies in $\Lambda$ will\ clearly entail sacrifice in terms of bias. The extent of the censoring would usually be greater the larger is $\lambda$ in absolute value, and would also depend on characteristics of $W$ such as its sparseness (see the Monte Carlo simulations in Section (ref)). But, presumably, the greater the censoring, the more sacrifice there will be in terms of bias-reduction. In fact, simulations reported in Section (ref) of the Supplement suggest that, when $\hat{\lambda }_{\mathrm{aML}}$ is outside $\Lambda$, often the unrestricted maximizer of $l(\lambda)$ is also outside $\Lambda$ (i.e., $\hat{\lambda}_{\mathrm{ML}} \neq\hat{\lambda}_{\mathrm{uML}}$). This suggests that, whenever $\Lambda_{\mathrm{a}}\neq\Lambda$, $\hat{\lambda}_{\mathrm{aML}}$ should be regarded as a modification to the QMLE that maximizes $l(\lambda)$ over $\Lambda_{\mathrm{a}}$, not over $\Lambda$.
Confidence intervals for a parameter of interest based on the QMLE may perform poorly if the data contains little information about the nuisance parameters. This may be the case, for instance, when the number of nuisance parameters is large relative to the sample size. The theory presented in the paper so far allows the construction of confidence intervals for $\lambda$ with accurate coverage even in the case of little information about $\beta$.
We say that a differentiable function is single-peaked on an open interval if, on that interval, it has a maximum and no other stationary points corresponding to minima or maxima (so the function is non-decreasing to the left of the peak, and non-increasing to the right of the peak). We shall see later in this subsection that $l_{\mathrm{a}}(\lambda)$ is single-peaked on $\Lambda_{\mathrm{a}}$ quite generally (Proposition (ref)), but before doing that we show how single-peakedness can be used to construct confidence intervals for $\lambda.$ If $l_{\mathrm{a}}(\lambda)$ is single-peaked on $\Lambda_{\mathrm{a}}$, then the cdf of $\hat{\lambda }_{\mathrm{aML}}$ admits the representation
for any $z\in\Lambda_{\mathrm{a}}$.\footnote{The notation $\Pr(\hat{\lambda }_{\mathrm{aML}}\leq z;\beta,\sigma^{2},\lambda)$ emphasizes that, of course, we are assuming that $y$ is generated by the SAR model ((ref)). But it is worth remarking that, for a general random vector $y$, we would still have $\Pr(\hat{\lambda}_{\mathrm{aML}}\leq z)=\Pr(y^{\prime}S^{\prime }(z)R(z)S(z)y\leq0)$; all that is required is that the likelihood ((ref)) is used for estimation.} That is, the cdf of $\hat{\lambda }_{\mathrm{aML}}$ at $z$ equals the cdf of the quadratic form $y^{\prime }S^{\prime}(z)R(z)S(z)y$ at zero. Thus, one can use some approximation to the cdf (at zero) of $y^{\prime}S^{\prime}(z)R(z)S(z)y$ to obtain an approximation to the cdf (at $z$) of $\hat{\lambda}_{\mathrm{aML}}$. Since the cumulant generating function of a quadratic form is available in closed form, at least under normality, one natural candidate is the Lugannani--Rice saddlepoint approximation Lugannani80. In fact, using a saddlepoint approximation to the cdf of a score to recover the cdf of the corresponding estimator had already been investigated in Daniels1983 Butler2007. The approximate cdf of $\hat {\lambda}_{\mathrm{aML}}$ can then be inverted to obtain confidence intervals for $\lambda$. This approach to constructing confidence intervals has been applied by Hillier2017 to the (unadjusted) QMLE of $\lambda$. Whilst revising the present paper we discovered that essentially the same approach had previously been suggested by Paige2009 for general quadratic estimating equations, and then by Jeganathan2015 specifically for some spatial models.
Let $\widetilde{\Pr}(\hat{\lambda}_{\mathrm{aML}}\leq z;\beta,\sigma ^{2},\lambda)$ denote the approximation to the cdf of $\hat{\lambda }_{\mathrm{aML}}$, obtained by the Lugannani--Rice formula (see Section (ref) of the Supplement for details). The parameters $\beta$ and $\sigma^{2}$ in the approximation can be replaced with their QMLEs given $\lambda$, $\hat{\beta}_{\mathrm{ML}}(\lambda)$ and $\hat{\sigma }_{\mathrm{aML}}^{2}(\lambda)$, to give the approximation $\widetilde{\Pr }(\hat{\lambda}_{\mathrm{aML}}\leq z;\lambda)\coloneqq \widetilde{\Pr}(\hat{\lambda }_{\mathrm{aML}}\leq z;\hat{\beta}_{\mathrm{ML}}(\lambda),\hat{\sigma }_{\mathrm{aML}}^{2}(\lambda),\lambda)$. Confidence intervals for $\lambda$ based on $\hat{\lambda}_{\mathrm{aML}}$ can then be constructed by inverting $\widetilde{\Pr}(\hat{\lambda}_{\mathrm{aML}}\leq z;\lambda)$. More specifically, replace $z$ with the observed value of $\hat{\lambda }_{\mathrm{aML}}$, say $\hat{\lambda}_{\mathrm{aML}}^{\mathrm{obs}}$, and let $\lambda_{1}\coloneqq \inf\{\lambda:\widetilde{\Pr}(\hat{\lambda}_{\mathrm{aML}} \leq\hat{\lambda}_{\mathrm{aML}}^{\mathrm{obs}};\lambda)=1-\alpha_{1}\}$ and $\lambda_{2}=\sup\{\lambda:$ $\widetilde{\Pr}(\hat{\lambda}_{\mathrm{aML}} \leq\hat{\lambda}_{\mathrm{aML}}^{\mathrm{obs}};\lambda)=\alpha_{2}\}$. Then $(\lambda_{1},\lambda_{2})$ is an approximate $\mathopen{}\mathclose\bgroup\originalleft( 1-\alpha_{1}-\alpha _{2}\aftergroup\egroup\originalright) \%$ two-sided confidence interval for $\lambda$.\footnote{Note that the set $\{\lambda:\alpha_{1}\leq\widetilde{\Pr}(\hat{\lambda }_{\mathrm{aML}}\leq z;\lambda)\leq\alpha_{2}\}$ is in general a union of intervals. It is an interval if $\widetilde{\Pr}(\hat{\lambda}_{\mathrm{aML} }\leq z;\lambda)$ is monotonic in $\lambda$. It seems reasonable to expect this monotonicity to hold in many cases, at least approximately.}
The Lugannani--Rice approximation we employ is, as is frequently the case, a normal based one. That is, it is constructed on the basis of the cumulant generating function of $y^{\prime}S^{\prime}(z)R(z)S(z)y$ that obtains if $\varepsilon\sim\mathrm{N}(0,I_{n})$. It is important to emphasize, however, that we intend the approximation to be used generally. Indeed, there is considerable evidence in the literature that a normal-based Lugannani--Rice approximation is typically very accurate for the distribution of a statistic that has a limiting normal distribution. There is also evidence that a normal-based Lugannani--Rice approximation works generally very well for the distribution of a quadratic form in nonnormal variables as long as we are far from the lower tail of the distribution Wood1993. This seems to be relevant in our case as we are only interested in the cdf of $y^{\prime }S^{\prime}(z)R(z)S(z)y$ at $0$, and $R(z)$ is indefinite.
The key condition for representation ((ref)) to hold is that $l_{\mathrm{a}}(\lambda)$ is single-peaked on $\Lambda_{\mathrm{a}}$. We now discuss this condition. The following very mild assumption allows us to make use of the results derived earlier for semisimple eigenvalues.\footnote{When $\Lambda_{\mathrm{a}}=\Lambda$, there is nothing to assume, because in that case $\Lambda$ does not contain any zero of $\det(S(\lambda))$. When $\Lambda_{\mathrm{a}}\supset\Lambda$, in most cases of practical interest the only zero of $\det(S(\lambda))$ in $\Lambda_{\mathrm{a}}$ is the one corresponding to the eigenvalue $1$, which, as pointed out earlier, is virtually always semisimple in applications.}
It turns out that if $W$ is symmetric, condition ((ref)) is satisfied for any $X$ (see Lemma (ref) in Appendix (ref)), and hence $\Lambda_{\mathrm{a}}$ is a.s. single-peaked in that case. If $W$ is not symmetric, the condition depends on both $W$ and $X$.\footnote{It is interesting to note that, for the unadjusted profile likelihood $l(\lambda)$, single-peakedness over $\Lambda$ holds if $n\mathrm{tr}(G^{2}(\lambda ))>\mathopen{}\mathclose\bgroup\originalleft[ \mathrm{tr}(G(\lambda))\aftergroup\egroup\originalright] ^{2}$ for all $\lambda\in\Lambda$, a condition that does not depend on $X$ Hillier2017.} Extensive numerical experimentation with nonsymmetric $W$ suggests that for the vast majority of pairs $(W,X)$ likely to be met in applications, $l_{\mathrm{a}}(\lambda)$ is a.s.\ single-peaked on $\Lambda_{\mathrm{a}}$. And, for pairs $(W,X)$ such that condition ((ref)) is not satisfied, the probability (for some distribution of $y$ that is absolutely continuous w.r.t. $\mu_{\mathbb{R}^{n}}$) that $l_{\mathrm{a}}(\lambda)$ is multi-peaked is generally very small. Note that the right hand side of representation ((ref)) is likely to provide a good approximation to the cdf of $\hat{\lambda}_{\mathrm{aML}}$ whenever the probability that $l_{\mathrm{a} }(\lambda)$ is multi-peaked is nonzero but small. The performance of the saddlepoint confidence intervals is assessed by numerical simulation in Section (ref). Naturally, the confidence intervals can be inverted to obtain hypothesis tests on $\lambda$, but we do not investigate the power properties of such tests here.
In this section we apply the score adjustment to a social interaction model with network fixed effects and contextual effects Lee2007b, Bramoulle2009, LeeLiuLin2010, Lin2015, Fortin2015. There are $R$ networks, with network $r$ having $m_{r}$ individuals. The model is
where $W_{r}$ is the weights matrix of network $r$, $\alpha_{r}$ is an unobserved network fixed effect, $\widetilde{X}_{r}$ is an $m_{r}\times \tilde{k}$ matrix of regressors, $\gamma$ and $\delta$ are $\tilde{k}\times1$ parameters. In terms of equation ((ref)), $y=(y_{1}^{\prime} ,\ldots,y_{R}^{\prime})^{\prime}$, $W=\bigoplus_{r=1}^{R}W_{r}$, $\beta =(\gamma^{\prime},\delta^{\prime},\alpha_{1},\ldots,\alpha_{R})^{\prime}$, $\varepsilon=(\varepsilon_{1}^{\prime},\ldots,\varepsilon_{R}^{\prime})^{\prime} $, and $X=(\widetilde{X},W\widetilde{X},\bigoplus_{r=1}^{R}\iota_{m_{r}})$, with $\widetilde{X}\coloneqq (\widetilde{X}_{1}^{\prime},\ldots,\widetilde{X}_{R} ^{\prime})^{\prime}$. We assume that $X$ has full column rank.
To avoid the incidental parameter problem that arises in the case of many networks of fixed dimensions, inference in this model usually proceeds by first eliminating the network fixed effects. Define the orthogonal projector $M_{r}\coloneqq I_{m_{r}}-\frac{1}{m_{r}}\iota_{m_{r}}\iota_{m_{r}}^{\prime}$, and let $F_{r}$ be the $m_{r}\times\mathopen{}\mathclose\bgroup\originalleft( m_{r}-1\aftergroup\egroup\originalright) $ matrix of orthonormal eigenvectors of $M_{r}$ corresponding to the eigenvalue $1$, so that $F_{r}^{\prime}F_{r}=I_{m_{r}-1}$, $F_{r}F_{r}^{\prime}=M_{r}$, and $F_{r}^{\prime}\iota_{m_{r}}=0$. In a likelihood framework, the standard procedure is that proposed by LeeLiuLin2010, which eliminates the fixed effects by premultiplying each equation in ((ref)) by $F_{r}^{\prime}$. Under the condition that each weights matrix $W_{r}$ has all row sums equal to one (i.e., $W_{r}\iota_{m_{r}}=\iota_{m_{r}}$ for all $r=1,\ldots,R$), we have $F_{r}^{\prime}W_{r}=F_{r}^{\prime}W_{r}F_{r} F_{r}^{\prime}$, and therefore the transformed model is
where $y_{r}^{\ast}\coloneqq F_{r}^{\prime}y_{r}$, $W_{r}^{\ast}\coloneqq F_{r}^{\prime} W_{r}F_{r},$ $\widetilde{X}_{r}^{\ast}\coloneqq F_{r}^{\prime}\widetilde{X}_{r}$, and $\varepsilon_{r}^{\ast}\coloneqq F_{r}^{\prime}\varepsilon_{r}$. Note that $\mathrm{var}(\varepsilon_{r}^{\ast})=I_{m-1}$ if $\mathrm{var}(\varepsilon _{r})=I_{m}$. We denote the profile (quasi) log-likelihood for $(\sigma ^{2},\lambda)$ based on the transformed model ((ref)) by $l_{\mathrm{LLL}}(\sigma^{2},\lambda)$. If the condition $W_{r}\iota_{m_{r} }=\iota_{m_{r}}$, for each $r=1,\ldots,R$, is not satisfied, the transformation by $F_{r}^{\prime}$ does not yield a reduced form and hence a likelihood. In that case, only non-likelihood procedures, such as the GMM of LiuLee2010, are available.
The score adjustment discussed in the present paper provides an alternative likelihood solution to the large-$R$ incidental parameter problem. It provides a large-$R$ consistent estimator of all model parameters (see Remark (ref)), like the LeeLiuLin2010 procedure, but, contrary to LeeLiuLin2010, it does not require the constant row sums condition. The score adjustment is performed exactly as for the general SAR model, treating the $\alpha_{r}$'s as parameters, profiling out the parameter $\beta =(\gamma^{\prime},\delta^{\prime},\alpha_{1},\ldots,\alpha_{R})^{\prime}$ as in equation ((ref)), and recalling that the expectation ((ref)) does not depend on $\beta$.
The next proposition shows that the score adjustment method is equivalent to the LeeLiuLin2010 method, if the latter applies (i.e., $W_{r} \iota_{m_{r}}=\iota_{m_{r}}$ for all $r$) and there are no covariates (i.e., the model is $y_{r}=\lambda W_{r}y_{r}+\alpha_{r}\iota_{m_{r}}+\sigma \varepsilon_{r}$, $r=1,\ldots,R$).
When covariates are present, the estimators of $\sigma^{2}$ and $\lambda$ (and hence of $\beta$) obtained by the score adjustment method are different from the ones obtained by the LeeLiuLin2010 method. Both methods solve the large-$R$ incidental parameter problem, but, in addition to not requiring the constant row sum condition, the score adjustment approach also implements a correction that deals with the nuisance parameters $\gamma$ and $\delta$. The two methods are compared by simulation in Section (ref) .\footnote{Our focus is on finite sample results. From the point of view of first-order asymptotics, we expect the two estimators to be equivalent under the conditions in LeeLiuLin2010. Such conditions include, in particular, the fact that $\tilde{k}$ does not vary with $n$; if $\tilde{k}$ increased with $n$, the adjusted QMLE might have an advantage even from the point of view of first-order asymptotics.}
The spatial error model
is considerably less popular than the SAR model ((ref)) in economic applications, but, as we shall see, provides important motivation for the score adjustment considered in this paper. We now briefly consider this model, leaving all details to Section (ref) of the Supplement. It is convenient to generalize the error structure in model ((ref)) to $A(\theta)u=\sigma\varepsilon$, where $A(\theta)$ is a square matrix that is invertible for any value of the parameter $\theta$ in a subset $\Theta$ of $\mathbb{R}^{p}$. In the case of the spatial error model, we may take $A(\theta)=S(\lambda).$
The profile (quasi) log-likelihood for $(\sigma^{2},\theta)$ based on the assumption $\varepsilon\sim\mathrm{N}(0,I_{n})$ is (up to an additive constant)
where $U(\theta)\coloneqq A^{\prime}(\theta)M_{A(\theta)X}A(\theta)$. As for the case of the SAR model, the score associated to ((ref)) can be (exactly) recentered assuming only that $\mathrm{E}(\varepsilon)=0$ and $\mathrm{var}(\varepsilon)=I_{n}$. The corresponding adjusted likelihood is
Remarkably, and contrary to the case of a SAR model, $l_{\mathrm{a}} (\sigma^{2},\theta)$ is a genuine likelihood. Hence, the corresponding score provides an unbiased and information unbiased estimating equation (under the assumptions that $\mathrm{E}(\varepsilon)=0$ and $\mathrm{var}(\varepsilon )=I_{n}$). In fact, $l_{\mathrm{a}}(\sigma^{2},\theta)$ is equivalent to the likelihood used for restricted ML (REML) estimation, which has been shown to be quite generally preferable to ML estimation in a regression model with correlated errors Thompson1962, Patterson1971, RahmanKing1997.
Maximization of ((ref)) for fixed $\theta$ gives $\hat{\sigma}_{\mathrm{a}\mathrm{ML}}^{2}(\theta)=\frac{1}{n-k}y^{\prime }U(\theta)y$, and thus the profile adjusted likelihood for $\theta$ only is
To fully understand the effect of the score adjustment, it is useful to relate the likelihood $l_{\mathrm{a}}(\theta)$ to invariance properties of the model. Assuming that the distribution of $\varepsilon$ does not depend on the parameters $\beta$ and $\sigma$ (and that the parameters are identifiable), it is straightforward to check that the model is invariant under the group $\mathcal{G}_{X}$ of transformations $y\rightarrow\kappa y+X\delta$ in the sample space, for any $\kappa>0$, any $\delta\in\mathbb{R}^{k}$, and for a fixed $X$ Lehmann2005. The likelihood $l_{\mathrm{a}}(\theta)$ corresponds to the density of a maximal invariant under the group $\mathcal{G}_{X}$ (and $l_{\mathrm{a} }(\sigma^{2},\theta)$ corresponds to the density of a maximal invariant under the group of transformations $y\rightarrow y+X\delta$, $\delta\in \mathbb{R}^{k}$). That is, using the adjusted likelihood $l_{\mathrm{a} }(\theta)$ corresponds to imposing that inference should be invariant with respect to the group under which the model itself is invariant, as advocated by the \textquotedblleft principle of invariance\textquotedblright.
We do not give details here, but it can be shown that, as $\lambda$ approaches any real zero of $\det\mathopen{}\mathclose\bgroup\originalleft( S(\lambda)\aftergroup\egroup\originalright) $, $l(\lambda)\rightarrow -\infty$ a.s., whereas $l_{\mathrm{a}}(\lambda)$ may either a.s. approach $-\infty$ or be a.s.\ bounded in a spatial error model.\footnote{Given the correspondence between $l_{\mathrm{a}}(\lambda)$ and the density of a maximal invariant under the group $\mathcal{G}_{X}$, the fact that $l_{\mathrm{a} }(\lambda)$ may be bounded as $\lambda\rightarrow\omega^{-1}$, for some nonzero real eigenvalue $\omega$ of $W$, provides an explanation of why the limiting power of an invariant test for $\lambda=0$ may be strictly in $(0,1)$ as $\lambda\rightarrow\omega^{-1}$ Martellosio2010. On the other hand, if $l_{\mathrm{a}}(\lambda)\rightarrow-\infty$, then the power of an invariant test for $\lambda=0$ may approach 0 or 1 as $\lambda\rightarrow\omega^{-1}$ depending on the location of the critical region.} Thus, as might have been expected, in the spatial error model the profile log-likelihood and its adjusted version behave very much as in the SAR model.\footnote{There is one important difference though. Contrary to the case of the SAR model, in the spatial error model $l_{\mathrm{a}}(\lambda)$ can never approach $+\infty$ a.s.\ as $\lambda$ approaches a real zero of $\det\mathopen{}\mathclose\bgroup\originalleft( S(\lambda)\aftergroup\egroup\originalright) $. This is because $l_{\mathrm{a}}(\lambda)$ is a genuine likelihood, so a.s.\ unboundedness from above of $l_{\mathrm{a}}(\lambda)$ as $\lambda$ approaches a zero of $\det\mathopen{}\mathclose\bgroup\originalleft( S(\lambda)\aftergroup\egroup\originalright) $ would imply the existence of a density that diverges to +$\infty$ almost everywhere on $\mathbb{R}^{n}$ as $\lambda$ approaches that zero.} So, in particular, it makes sense to define the adjusted QMLE of $\lambda$ on a support different from $\Lambda$, as in the SAR model. More generally, this implies that when the covariance parameter $\theta$ is restricted to a certain set, the REML estimator of $\theta$ does not necessarily respect that restriction, something that, to the best of our knowledge, has not been noted before in the literature.
We conduct Monte Carlo experiments to investigate the performance of the adjusted MLE of $\lambda$ in the SAR model. First, we consider the case of a single network, and then the case of the multiple networks model of Section (ref). The number of replications is $10^{6}$ in all experiments.
A Watts-Strogatz random graph is formed by rewiring with probability $p$ each link in a $h$-ahead $h$-behind circular matrix Watts1998.\footnote{A $h$-ahead $h$-behind circular matrix $A$ is the adjacency matrix of a graph made of $n$ vertices on a circle such that each vertex is linked to its $h$ nearest neighbors on each side. That is, $A(i,j)=1$ if $0<\mathopen{}\mathclose\bgroup\originalleft\vert i-j\aftergroup\egroup\originalright\vert \operatorname{mod}(n-1-h)\leq h$), and $A(i,j)=0$ otherwise.} The extreme cases $p=0$ and $p=1$ correspond, respectively, to a $h$-ahead $h$-behind circular matrix and to a Erd{\H{o}}s-R{\'{e}}nyi random graph. Watts-Strogatz graphs are popular in social network analysis because, for relatively small values of $p$, they can reproduce important characteristics of many real world networks, including high clustering and small distances between most nodes. In our first numerical experiment, we take $W$ in model ((ref)) to be the normalized adjacency matrix of (one realization of) a Watts-Strogatz graph. Normalization of $W$ is either a row normalization or a spectral normalization (by spectral normalization we mean that the matrix is rescaled by its spectral radius). Note that the two normalizations are equivalent when $p=0$, due to the fact that a $h$-ahead $h$-behind circular matrix has constant row sums. The matrix $X$ contains an intercept, $\tilde {k}$ regressors, and, to allow for contextual effects, the spatially lagged version of those regressors. That is, $X=(\iota_{n},\widetilde{X} ,W\widetilde{X})$, where $\widetilde{X}$ is $n\times\tilde{k}$. The matrix $\widetilde{X}$ is drawn in each repetition; half of its columns are drawn from independent $\mathrm{N}(0,I_{n})$ distributions, half from independent uniform distributions on $[0,1]$.\footnote{Given our theoretical framework, it could be argued that a design with $\widetilde{X}$ fixed across repetitions would be more appropriate. However, we prefer to vary $\widetilde{X}$ across repetitions, in an attempt to investigate the performance of the adjusted QMLE in an environment where $\widetilde{X}$ is random.} For $p>0$, $W$ is drawn once and then kept constant across repetitions. We set $\beta=\iota _{2\tilde{k}+1}$, $\sigma=1$, and $n=200$. The errors $\varepsilon_{i}$ are generated from independent standard normal distributions.
In the present context, $\hat{\lambda}_{\mathrm{ML}}$ is, subject to regularity conditions, consistent and asymptotically normal as $n\rightarrow \infty$ Lee2004, and $\hat{\lambda}_{\mathrm{ML}}$ and $\hat{\lambda }_{\mathrm{aML}}$ are first-order asymptotically equivalent Yu2015. Table (ref) compares the finite sample performance of $\hat{\lambda}_{\mathrm{ML}}$ and $\hat{\lambda}_{\mathrm{aML}}$ for a range of values of $p$, $h$, and $\lambda$, when $\tilde{k}=2$. The columns headed by $\hat{\lambda}_{\mathrm{ML}}$ and $\hat{\lambda}_{\mathrm{aML}}$ give the bias(s.d.) of the estimators. $\Delta\%$ stands for percentage change (of the absolute bias or RMSE) from $\hat{\lambda}_{\mathrm{ML}}$ to $\hat{\lambda }_{\mathrm{aML}}$. For the case of row normalization, we also report $\omega_{2}^{-1}$, the percentage of times when $\hat{\lambda}_{\mathrm{aML} }>1$, denoted by $\%(\hat{\lambda}_{\mathrm{aML}}>1)$, and bias and s.d. of $\bar{\lambda}_{\mathrm{aML}}$ (which, recall, is the estimator that is set to 1 when $\hat{\lambda}_{\mathrm{aML}}>1$). According to the results in Section (ref), in the present setting $\Lambda_{\mathrm{a}}$ must be the same in each repetition, with its right endpoint being equal to $\omega_{2}^{-1}$ when $W$ is row normalized, equal to 1 when $W$ is spectrally normalized (and $p>0$). Thus, $\hat{\lambda}_{\mathrm{aML}}$ may be greater than $1$ when $W$ is row normalized, while $\hat{\lambda }_{\mathrm{aML}}<1$ in all repetitions when $W$ is spectrally normalized. The results in Table (ref) show that, when $p=0$, $\hat{\lambda}_{\mathrm{aML}}$ provides a significant improvement compared to $\hat{\lambda}_{\mathrm{ML}}$ both in terms of bias and RMSE. As $p$ increases, the reduction in bias afforded by $\hat{\lambda}_{\mathrm{aML}}$ increases but at the cost of greater variability, for both normalizations. The bias of $\hat{\lambda}_{\mathrm{aML}}$ is small, unless $h$ is large and $p$ is small. Note that $\omega_{2}^{-1}$ increases with $p$ and $h$ (cf. the paragraph after Example (ref)), and $\%(\hat{\lambda }_{\mathrm{aML}}>1)$ is nondecreasing in $p$, $h,$ and $\lambda$, and can be very large when $p$, $h,$ and $\lambda$ are large. As $\%(\hat{\lambda }_{\mathrm{aML}}>1)$ increases, $\bar{\lambda}_{\mathrm{aML}}$ becomes less biased, and more variable, compared to $\hat{\lambda}_{\mathrm{aML}}$. \enlargethispage{\baselineskip} \enlargethispage{\baselineskip}
Table (ref) analyzes the impact of the number of regressors. We take $\lambda=.5$ and $h=5$. For these values of $\lambda$ and $h$, $\mathrm{Pr}(\hat{\lambda}_{\mathrm{aML}}>1)$ is either zero or negligible in the row normalized case, so, contrary to the previous table, we do not report $\omega_{2}^{-1}$, $\%(\hat{\lambda}_{\mathrm{aML}}>1)$, and $\bar{\lambda }_{\mathrm{aML}}$. As expected, the bias reduction afforded by $\hat{\lambda }_{\mathrm{aML}}$ increases as $\tilde{k}$ increases. As $\tilde{k}$ increases, the relative performance of $\hat{\lambda}_{\mathrm{aML}}$, compared to $\hat{\lambda}_{\mathrm{ML}}$, improves also in terms of RMSE if $p$ is not large.
\thispagestyle{empty}\newgeometry{left=3cm, right=1.5cm, top=2cm, bottom=1cm}
\restoregeometry
We now generate data according to model ((ref)). For simplicity, we consider a balanced case: each of the $R$ networks has the same number, $m$, of individuals. The matrices $W_{r}$ and $\widetilde{X}$ are generated as in Section (ref). That is, each $W_{r}$ is the normalized adjacency matrix of (a realization of) a Watts-Strogatz graph, and the regressors are drawn in each repetition, $\tilde{k}/2$ of them from a standard normal distribution, and the other $\tilde{k}/2$ from a $\mathrm{Unif}(0,1)$ distribution. Errors and the fixed effects $\alpha_{r}$ are drawn (in each repetition) from $\mathrm{N}(0,1)$ distributions, and we set $\gamma =\delta=\iota_{\tilde{k}}$ and $\sigma=1.$\footnote{The simulated model is one in which fixed effects and regressors are random and independent of each other, and $W$ is non-stochastic. Note that we could allow for some correlation between fixed effects and regressors, but this is not needed for our purposes.}
Table (ref) compares $\hat{\lambda}_{\mathrm{aML}}$ to the estimator $\hat{\lambda}_{\mathrm{LLL}}$ obtained by maximizing the LeeLiuLin2010 likelihood $l_{\mathrm{LLL}}(\sigma^{2},\lambda)$, for a range of values of $R$, $m$, and $\tilde{k}$.\footnote{Some results for larger $R$ are given in Table (ref) in the Supplement.} We choose $\lambda=0.5$, and set the parameters of the Watts-Strogatz random graph to $h=5$ and $p=0.2$. When $W$ is row normalized, $\hat{\lambda}_{\mathrm{aML}}$ performs similarly to $\hat{\lambda}_{\mathrm{LLL}}$.\footnote{Note that, as for the design in Section (ref), in the present setting the set $\Lambda_{\mathrm{a}}$ is the same in all repetitions.} In fact, it has lower bias (but note that the bias of $\hat{\lambda}_{\mathrm{LLL}}$ is already very small in most of the cases considered in the table) and slightly slower RMSE. As expected, the relative performance of $\hat{\lambda}_{\mathrm{LLL}}$ improves as $\tilde{k}$ increases.\footnote{Similarly to the cross-sectional case, the reduction in RMSE is higher for lower values of $p$ (for example, in the extreme case $p=0$, $\Delta\%\text{RMSE}$ is -4.417, -14.707, -24.278 when $\tilde{k}$ is $2,6,10$, respectively, and $R=10$, $m=20$).} When $W$ is spectrally normalized, $\hat{\lambda}_{\mathrm{LLL}}$ cannot be obtained, so only results for $\hat{\lambda}_{\mathrm{aML}}$ are reported. The simulation results show that $\hat{\lambda}_{\mathrm{aML}}$ performs very satisfactorily even in that case.
Table (ref) reports empirical coverages of Wald confidence intervals based on first-order asymptotic normality of $\hat{\lambda }_{\mathrm{LLL}}$ (in the columns headed by $\mathrm{W}_{\mathrm{LLL}}$), empirical coverages of Wald confidence intervals based on first-order asymptotic normality of $\hat{\lambda}_{\mathrm{aML}}$ (in the columns headed by $\mathrm{W}_{\mathrm{aML}}$), and empirical coverages of saddlepoint confidence intervals based on $\hat{\lambda}_{\mathrm{aML}}$ (in the columns headed by $\mathrm{s}_{\mathrm{aML}}$). Since $\hat{\lambda}_{\mathrm{LLL}}$ is not available when $W$ is spectrally normalized, we only report the case of row normalized $W$.\footnote{When $W$ is spectrally normalized, coverages of the confidence intervals based on $\hat{\lambda}_{\mathrm{aML}}$ are similar to the case of row normalization.} The nominal size is 95%. We consider equi-tailed two-sided confidence intervals, and right-sided confidence intervals of the form $(-\infty,\lambda_{U})$, where $\lambda_{U}$ is a suitably selected upper end-point.\footnote{The 95% Wald confidence intervals are $\hat{\lambda}\pm1.96\sqrt{\hat{v}}$ (two-sided) and $(-\infty ,\hat{\lambda}+1.645\sqrt{\hat{v}})$ (right-sided), where $\hat{\lambda}$ is either $\hat{\lambda}_{\mathrm{LLL}}$ or $\hat{\lambda}_{\mathrm{aMLE}}$, and $\hat{v}$ denotes the asymptotic variance given in Proposition 6.1 of LeeLiuLin2010 and evaluated at the LLL or aMLE estimates of $\lambda,\beta,\sigma^{2}$.} The errors $\varepsilon_{ri}$ are generated independently from either (a) a standard normal distribution, (b) a gamma distribution with shape parameter 1 and scale parameter 1, demeaned by the population mean. Mean, variance, skewness, and kurtosis are $0,1,0,3$ in case (a) and $0,1,2,9$ in case (b). The Wald confidence intervals based on $\hat{\lambda}_{\mathrm{aML}}$ offer an improvement over Wald confidence intervals based on $\hat{\lambda}_{\mathrm{LLL}}$ when $\tilde{k}$ is not too small, and particularly in the case of right sided confidence intervals. The Lugannani--Rice approximation delivers a further improvement, and indeed the coverages of the saddlepoint confidence intervals based on $\hat{\lambda }_{\mathrm{aML}}$ are excellent in all cases considered in the table, even under the gamma distribution. Conversely, the empirical coverage of the Wald confidence intervals based on $\hat{\lambda}_{\mathrm{LLL}}$ is acceptable in the two-sided case when $R=m=30$, but quickly deteriorates as $\tilde{k}$ increases or $n$ decreases, and is considerably worse in the right-sided case. Of course, one minus a coverage in the table gives the size of the test for $\lambda=0$ obtained by inverting the confidence interval. Table (ref) in the Supplement reports coverages under other severely non-normal distributions; again, the saddlepoint confidence intervals based on the adjusted QMLE are very accurate in all cases considered.
Recentering the profile score for a parameter of interest is one possible way to deal with nuisance parameters. In this paper, we have applied this general principle to the estimation of the autoregressive parameter $\lambda$ in a spatial autoregression. The resulting adjusted QMLE for $\lambda$ successfully reduces the bias in the QMLE, provides confidence intervals with excellent coverage properties (and hence tests for $\lambda$ with excellent size properties) even when the dimension of the nuisance parameter is large, and is as straightforward to compute as the original QMLE. The adjusted QMLE can also solve the incidental parameter problem that occur, for example, in social network models with network fixed effects. However, due to the fact that the parameter space for $\lambda$ is usually restricted to a certain interval, the spatial autoregressive setting presents challenges for the score adjustment procedure that do not arise in other models. Namely, the distributions of the QMLE and of its adjusted version can be supported on different intervals, which means that a comparison between the two estimators is not straightforward. Our simulations suggest that the adjusted QMLE generally performs better than the QMLE in terms of RMSE, particularly in models with a large number of covariates.
This paper has focused on a simple version of a spatial autoregressive model. In empirical applications, it is typically desirable to extend the model in various directions. For example, one may want to allow for endogeneity coming from $X$ or $W$, or for model errors that are subject to a spatial autoregressive structure themselves. Such extensions would not preclude the use of the score adjustment procedure, but would typically imply that the expectation of the profile score for $\lambda$ depends on nuisance parameters. In that case, as discussed in McCullagh1990, the expectation would need to be obtained numerically, rather than analytically, and the resulting estimating equation would be unbiased up to some order, rather than being exactly unbiased.
It is also worth mentioning that it should be possible to use the score adjustment procedure in conjunction with modifications to the QMLE that allow for unknown heteroskedasticity LiuYang2015. Finally, the adjusted QMLE should be effective also in models where the number of regressors $k$ is increasing with the sample size GuptaRobinson2016. If $k$ increases sufficiently quickly with $n$, then the adjusted QMLE should have advantages, with respect to the QMLE, even in terms of first-order asymptotics.