The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
68,577 characters
Saddlepoint approximations for spatial panel data models
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if11
{
\title{\bf Saddlepoint approximations for \\ spatial panel data models}
\author{Chaonan Jiang\thanks{
Research Center for Statistics and Geneva School of Economics and Management, University of Geneva, Blv. Pont
d'Arve 40, 1211 Geneva, Switzerland, \url{[email removed]}}, \hspace{0.2cm}
Davide La Vecchia \thanks{Research Center for Statistics and Geneva School of Economics and Management, University of Geneva, Blv. Pont
d'Arve 40, 1211 Geneva, Switzerland, \url{[email removed]}}, \hspace{0.2cm} Elvezio Ronchetti \thanks{Research Center for Statistics
and Geneva School of Economics and Management, University of Geneva, Blv. Pont d'Arve 40,
1211 Geneva, Switzerland, \url{[email removed]}}, \hspace{0.2cm} Olivier Scaillet
\thanks{Geneva Finance Research Institute, Geneva School of Economics and Management, University of Geneva and Swiss Finance Institute, Blv. Pont d'Arve 40,
1211 Geneva, Switzerland, \url{[email removed]}}}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Saddlepoint approximations for \\ spatial panel data models}
\end{center}
\medskip
} \fi
\bigskip
\begin{abstract}
We develop new higher-order asymptotic techniques for the Gaussian maximum likelihood
estimator in a spatial panel data model, with fixed effects, time-varying covariates, and spatially correlated errors. Our saddlepoint density and tail area approximation
feature \textit{relative error} of order $O(1/(n(T-1)))$ with $n$ being the cross-sectional dimension and $T$ the time-series dimension. The main theoretical tool is the tilted-Edgeworth technique in a non-identically distributed setting. The density
approximation is always non-negative, does not need resampling, and is accurate in the tails.
Monte Carlo experiments on density approximation and testing in the presence of nuisance parameters illustrate the good performance of our approximation over first-order asymptotics and Edgeworth expansions. An empirical application to the investment-saving relationship in OECD (Organisation for Economic Co-operation and Development) countries shows disagreement between testing results based on first-order asymptotics and saddlepoint techniques.
\end{abstract}
\noindent
{\it Keywords:} Higher-order asymptotics, investment-saving, random field, tail area.
\vfill
\newpage
\spacingset{1.5}
\addtolength{\textheight}{.5in}
\newpage
\section{Introduction} \label{Sec: Intro}
Accounting for spatial dependence is of interest both from an applied and a theoretical point of view. Indeed,
panel data with spatial cross-sectional interaction enable empirical researchers to take into account the temporal dimension and, at the same time, control for the
spatial dependence. From a theoretical point of view, the special features of panel data with spatial effects
present the challenge to develop new methodological tools.
Much of the machinery for conducting statistical inference on panel data models has been established under
the simplifying assumption of cross-sectional independence. This assumption may be
inadequate in many cases. For instance,
correlation across spatial data comes typically from competition, spillovers, or aggregation. The presence of such a
correlation might be anticipated in observable variables and/or in the unobserved disturbances in a statistical model and
ignoring it can have adverse effects on routinely-applied inferential procedures.
See, e.g., \citet{GG10}, \citet{R12}, \citet{C15}, \citet{CW15}, and recently \citet{CMW19}
for book-length discussions in the statistical literature. In the econometric literature, see, e.g.,
\citet{KKP07}, \citet{LY10}, \citet{RR14_ET}, \citet{RR15}, and, for book-length presentations, \citet[Ch. 13]{BaltagiBook}, \citet{An13}, and \citet{KP17}.
Different nonparametric, semiparametric, and parametric approaches have been proposed to incorporate cross-sectional
dependence in panel data models. The choice on the modeling approach depends on the time series
($T$) and cross-sectional ($n$) dimensions.
A nonparametric approach
is only feasible when $T$ is large relative to $n$. In other situations, typically when $T$ is very small (e.g., $T=2$) and $n$ is large, semiparametric models have been employed, including time varying regressors (namely factor models)
and spatial autoregressive component, when information on spatial distances is available.
Least squares and quasi maximum-likelihood estimator represent the main popular tools for estimation
within this setting.
When both
$T$ and $n$ are small, the parametric approach is the sensible choice and
(Gaussian) likelihood-based procedures are applied to define the maximum likelihood estimator (MLE).
{There is a vast literature on the MLE for spatial autoregressive models, an early reference being \citet{O75}. The derivation of the first-order asymptotics is available in the econometric literature; we refer to the seminal paper by \citet{L04}. For the class of spatial autoregressive processes, with fixed effects, time-varying covariates, and spatially correlated errors that we consider in this paper, the first-order asymptotic results for the Gaussian MLE are available in \citet{LY10}, where the authors derive asymptotic approximations (the exact finite-sample distribution
being intractable), when the cross-sectional dimension $n$ is large and $T$ is finite or large.}
{The main issue related to first-order asymptotic approximations is that, when $n$ is not very large, such approximations may be unreliable: alternatives are highly recommended. \citet{BU07} provide analytic formulae for the second-order bias and mean squared error of the MLE for the spatial parameter $\lambda$, in a Gaussian model. \citet{B13} and \citet{Ya05} extend these approximations to include also exogenous explanatory variables, which remain valid also when the process is not Gaussian. \citet{RR14_EJ,RR14_ET} develop Edgeworth-improved tests for no
spatial correlation in spatial autoregressive models for pure cross-sectional data based on least squares estimation and Lagrange
multiplier tests. Moreover, \citet{RR15} work on the concentrated likelihood and derive an Edgeworth expansion for MLE of $\lambda$
in the setting of a first-order spatial autoregressive panel data model, with fixed effects and without covariates. \citet{HM18} (see their \S 6) and \citet{HM20} (see their \S 3.5) propose saddlepoint approximations for the profile likelihood estimator of $\lambda$.}
Resampling methods are also available alternatives to improve on the first-order asymptotics,
achieving higher-order asymptotic refinements in terms of {\it absolute} error.
However, it requires
either a bias correction or an asymptotically pivotal statistics;
see \citet{Hall92} and \citet{Horowitz01}
in the i.i.d.\ setting. To the best of our knowledge, for spatial panel models considered in this paper, such results are not
available.
The aim of this paper is to introduce \textit{saddlepoint approximations}
for parametric spatial autoregressive panel
data models with fixed effects and time-varying covariates. They overcome the problems mentioned above by means
of the tilted-Edgeworth technique. For general references on saddlepoint approximations in the i.i.d.\ setting,
see the seminal paper of \cite{D54} and the book-length presentations of
\cite{FR90}, \cite{Jensen95}, \cite{K06},
and \cite{brazzaleetal2007}. For a result about testing on spatial dependence, see \citet{T02}, and for developments
in time series models, see \citet{LR19}.
We remark that we could cast the methodology of this paper into the framework of statistical analysis of random fields on a network graph, where
the underlying, known, network graph describes the spatial structure of the stochastic process; see e.g. \citet{K09} Ch. 8 for a book-length introduction. In \S 2, we briefly comment on this approach. For the ease-of-reference to the extant econometric literature,
in the rest of the paper, we prefer to stick to the econometric notation and terminology of spatial panel data models.
The paper is organized as follows. In \S\ref{Sec: Motivation}, we provide a
motivating example.
\S\ref{Sec: Setting} defines the general model setting and the estimation method, whereas the detailed methodology is presented in
\S\ref{Sec: Methodology}.
In particular in \S\ref{Sec: Links_econometrics} we provide a detailed discussion about the connections with the econometric literature.
Algorithms and computational aspects are
discussed in \S\ref{Sec: algm}.
\S\ref{Sec: MC} provides numerical comparison with other methods and, in \S\ref{Sec: test_nuisance_parameters}, we tackle the problem of testing in the presence of nuisance parameters.
In \S\ref{Sec: Real}, we present an empirical application. The online
Supplementary Material contains detailed derivations, technical appendices and additional numerical experiments.
\vspace{-0.5cm}
\section{Motivating example} \label{Sec: Motivation}
We motivate our research by a Monte Carlo (MC) exercise illustrating the low accuracy of the routinely applied first-order asymptotics in the setting of spatial panel data model. We consider the model:
\begin{equation}
\begin{aligned}
Y_{nt} &= \lambda_0 W_n Y_{nt} + X_{nt}\beta_0 + c_{n0} + E_{nt},& \\
E_{nt} &=\rho_0 M_n E_{nt} +V_{nt}, & \quad t=1,....,T, &
\end{aligned}
\label{Eq: GeneralY}
\end{equation}
where $Y_{nt} = (y_{1t},y_{2t},...,y_{nt})$, $X_{nt}$ is an $n \times k$ matrix of non stochastic time-varying regressors, $c_{n0}$ is an $n \times 1$ vector of fixed effects, and $V_{nt} =(v_{1t},v_{2t},..,v_{nt})'$ are $n \times 1$ vectors with $v_{it} \sim \mathcal{N}(0,\sigma_0^2)$, i.i.d. across $i$ and $t$. The matrices $W_n$ and $M_n$ are weighting matrices
describing the spatial dynamics. Following the literature, we label this model SARAR(1,1) to emphasize the spatial dependence in both the response variable $Y_{nt}$ and the error $E_{nt}$.
As in the MC example in \citet{LY10} p.\ 172, we generate samples from (\ref{Eq: GeneralY}) using $\theta_0 = (\beta_0,\lambda_0,\rho_0,\sigma_0^2)' = (1.0,0.2,0.5,1)'$, $T=5$, and $k=4$ covariates. The quantities $X_{nt}$, $c_{n0}$ and $V_{nt}$ are generated from independent standard normal distributions and, as it is customary in the econometric literature, we set $W_n = M_n$, where the off-diagonal elements are different from zero, while the diagonal elements are all zero. We consider two sample sizes: $n=24$ (small sample) and $n=100$ (moderate/large sample). The choice of $n=24$ is related to the empirical data analysis that we consider in \S \ref{Sec: Real}, where we apply the model in (\ref{Eq: GeneralY}) to conduct inference on the investment-saving relation for the 24 OECD (Organisation for Economic Co-operation and Development) countries. Similar sample sizes arise in many
real data analyses, where panel datasets contain a limited number of cross-sectional units, e.g.,
sampling can be expensive and/or time consuming, as it is typically the case in field studies.
As it is customary in the statistical/econometric software, we illustrate the inference issues related to the use of the first-order asymptotics by means of three different spatial weight matrices: Rook matrix, Queen matrix, and Queen matrix with torus. In Figure \ref{Fig: Wn1}, we display the geometry of $Y_{nt}$ as implied by each considered spatial matrix: the plots highlight that different matrices imply different spatial relations. For instance, we see that the Rook matrix implies fewer links than the Queen matrix. Indeed, the Rook criterion defines neighbours by the existence of a common edge between two spatial units, whilst the Queen criterion is less rigid and defines neighbours as spatial units sharing an edge or a vertex. Besides, we may interpret $\{Y_{nt}\}$ as a $n$-dimensional random field on the network graph which describes the known underlying spatial structure. Then, $W_n$ represents the weighted adjacency matrix (in the spatial econometrics literature, $W_n$ is called contiguity matrix). In Figure \ref{Fig: Wn1}, we display the geometry of a random field on a regular lattice (undirected graph).
In the real data example of \S \ref{Sec: Real}, we consider a random field over a manifold (a sphere), providing two additional examples for $W_n$.
\begin{figure}[htbp]
\begin{center}
\begin{tabular}{l}
\hspace{2.cm} Rook \hspace{3.75cm} Queen \hspace{3.8cm} Queen torus \\
\includegraphics[width=0.25\textwidth, height=0.25\textheight]{Wn24_Rook.eps}\hspace{1.4cm}
\includegraphics[width=0.25\textwidth, height=0.25\textheight]{Wn24_Queen.eps}\hspace{1.5cm}
\includegraphics[width=0.25\textwidth, height=0.25\textheight]{Wn24_QueenTorus.eps}\\
\end{tabular}
\caption{Different types of neighboring structure for $Y_{nt}$, as implied by different types of $W_n$ matrix, for $n=24$.}
\label{Fig: Wn1}
\end{center}
\end{figure}
To illustrate the inferential issues entailed by the use of first-order asymptotics, for each type of $W_n$, we generate a sample of $n$ observations. Since $c_{n0}$ creates an incidental parameter issue, we
eliminate it by the standard differentiation procedure, and for each MC run we estimate the model parameter $\theta$ using the transformation approach of \citet{LY10}, with maximum likelihood estimation method; we refer to the R package \texttt{spml} for implementation details. We set the MC size to 5000.
We illustrate graphically the behavior of the first-order asymptotic theory in finite sample by comparing the distribution of $\hat\lambda$ to the Gaussian asymptotic distribution (see \S \ref{Sec: M_fun}
for details). Via QQ-plot analysis, Figure \ref{fig1} shows that the Gaussian approximation can be either too thin or too thick in the tails with respect to the ``exact'' distribution. For instance, when $n=24$ and $W_n$ is rook, the Gaussian quantiles are larger than the ``exact'' ones in the left tail, while we observe the opposite phenomenon in the right tail. Similar considerations hold for the other types of $W_n$. The more complex is
the geometry of $W_n$ (e.g., $W_n$ has Queen structure) the more pronounced are the departures from the Gaussian. For $n=100$, and $W_n$ Rook, the MLE displays a distribution which is in line with the Gaussian one (see bottom left panel). However, when $W_n$ becomes more complex (e.g., Queen with torus),
larger departures in the tails are still evident. In Appendix D.1, we illustrate that similar conclusions are available also for the simpler SAR(1) model:
\begin{equation}
\begin{aligned}
Y_{nt} &=& \lambda_0 W_n Y_{nt} + c_{n0} + V_{nt}, \quad {\text{for}} \quad t=1,2,
\end{aligned}
\label{Eq: RR}
\end{equation}
where $\theta_0=(\lambda_0,\sigma_0^2)'$. More generally, unreported results suggest that, in the considered SARAR setting, the ``exact'' and the asymptotic distribution, as well as the saddlepoint approximation, agree for the considered types of $W_n$, when $n\geq250$. \\
\begin{figure}[htbp]
\begin{center}
\begin{tabular}{l} \hspace{3cm} Rook \hspace{3.1cm} Queen \hspace{2cm} Queen torus \\
\begin{turn}
{90} \hspace{2.1cm} $n = 24$
\end{turn}
\includegraphics[width=0.9\textwidth, height=0.175\textheight]{spml_n24T5_M.eps} \\
\begin{turn}
{90} \hspace{2.1cm} $n=100$
\end{turn}
\includegraphics[width=0.9\textwidth, height=0.175\textheight]{spml_n100T5_M.eps}
\end{tabular}
\caption{SARAR(1,1) model: QQ-plot vs normal of the MLE $\hat\lambda$, for different sample sizes ($n=24$ and $n=100$), $\lambda_0=0.2$, and different types of $W_n$ matrix.}
\label{fig1}
\end{center}
\end{figure}
\section{Model setting and estimation method} \label{Sec: Setting}
Let us consider a random field described by the {SARAR(1,1) model in (\ref{Eq: GeneralY}).
{We label by $P_{\theta_0}\in \mathcal{P}$, with $\theta_0 \in \Theta \subset \mathbb{R}^d$, the actual underlying distribution, which is characterized by $\theta_0 = (\beta_0,\lambda_0,\rho_0,\sigma_0^2)'$, the true parameter value.} The matrix $W_n$ is an $n \times n$ nonstochastic spatial weight matrix that generates the spatial dependence
on $y_{it}$ among cross sectional units. The matrix $X_{nt}$ is an $n \times k$ matrix of non stochastic time varying regressors, and $c_{n0}$ is an $n \times
1$ vector of fixed effects. Similarly, $M_n$ is an $n \times n$ spatial weight matrix for the disturbances --- quite often $W_n = M_n$. Moreover, we define $S_n(\lambda)= I_n - \lambda W_n$, and analogously
$R_n(\rho) = I_n - \rho M_n$.
The vector $c_{n0}$ introduces an incidental parameter problem; see \citet{LY10} and \citet{RR15}. To cope with this issue, we follow the standard approach, and we transform the model in order to derive consistent estimator for the model parameter $\theta=(\beta',\lambda,\rho,\sigma^2)'$ and $\theta\in \Theta \subset \mathbb{R}^{d}$. To achieve the goal, we first eliminate the individual effects by the deviation from the time-mean operator $J_T = (I_T - \frac{1}{T} l_T l'_T)$, where
$I_T$ is the $T \times T$ identity matrix, and $l_T = (1,...,1)'$, namely the $T\times 1$ vector of ones. Without creating linear dependence in the resulting
disturbances, we adopt the transformation introduced by \citet{LY10}.
First, let the orthonormal eigenvector matrix of $J_T$ be
$[F_{T,T-1}, \frac{1}{\sqrt{T}} l_T]$, {where $[\cdot]$ represents a matrix horizontal concatenation and} $F_{T,T-1}$ is the $T \times (T-1)$ submatrix
corresponding to the unit eigenvalues.
Then, for any $n \times T$ matrix $[Z_{n1}, ... , Z_{nT}]$, we
define the transformed $n\times (T - 1)$ matrix
$
[Z^*_{n1}, ... , Z^*_{nT}]=[Z_{n1}, ... , Z_{nT}]F_{T,T-1}.
$
Similarly, we define $X^*_{nt} = [X^*_{nt,1},X^*_{nt,2},...,X^*_{nt,k}]$. Thus, we transform the model in (\ref{Eq: GeneralY}) and we obtain:
\begin{equation}
\begin{aligned}
Y^*_{nt} &= \lambda_0 W_n Y^*_{nt} + X^*_{nt}\beta_0 + E^*_{nt}, &\\
E^*_{nt} &= \rho_0 M_n E^*_{nt} +V^*_{nt}, &\quad t=1,2,...,T.&
\end{aligned}
\label{Eq: GeneralY_Transf}
\end{equation}
Since $
\left(V^{*'}_{n1},...,V^{*'}_{n(T-1)}\right)' = \left(F'_{T,T-1} \otimes I_n\right)\left(V'_{n1},...,V'_{n(T-1)}\right)',
$
and the $v_{it}$ are i.i.d., we have
$$
\mathbb{E}\left[\left(V^{*'}_{n1},...,V^{*'}_{n(T-1)}\right)'\left(V^{*'}_{n1},...,V^{*'}_{n(T-1)}\right)\right] = \sigma_0^2 I_{n(T-1)},
$$
where $\mathbb{E}[\cdot]$ represents the expectation taken w.r.t. $P_{\theta_0}$.
{Now, the Gaussian assumption on the innovation terms implies that $v^{*}_{it}$ are independent for all $i$ and $t$; without this
assumption, they would be simply uncorrelated. See \citet{LY10} p. 167.}
Thus, defining $\zeta =(\beta',\lambda,\rho)'$, the log-likelihood is:
\begin{eqnarray*}
\ln L_{n,T}(\theta) = \ell_{n,T}(\theta) &=& - \frac{n(T-1)}{2} \ln (2\pi\sigma^2) + (T-1) [\ln \vert S_n(\lambda) \vert + \ln \vert R_n(\rho) \vert] \nonumber \\
& & \qquad - \frac{1}{2\sigma^2} \sum_{t=1}^{T-1} V^{*'}_{nt}(\zeta)V^*_{nt}(\zeta),
\label{Eq: LogL1}
\end{eqnarray*}
where $V^*_{nt}(\zeta) = R_n(\rho) [S_n(\lambda) Y^*_{nt} - X^*_{nt} \beta].$ {As remarked in \citet{LY10}, the function $L_{n,T}$ has a conditional likelihood interpretation: it is the likelihood conditional on the time average $\sum_{t=1}^{T} Y_{nt}/T$, which is a sufficient statistic for $c_{n0}$, under normality.}
We rewrite $\ell_{n,T}(\theta)$ in terms of a quadratic form in $\tilde{V}_{nt}(\zeta)$ as:
\begin{eqnarray}
\ell_{n,T}(\theta) &=& - \frac{n(T-1)}{2} \ln (2\pi\sigma^2) + (T-1) [\ln \vert S_n(\lambda) \vert + \ln \vert R_n(\rho) \vert] \nonumber \\
& & - \qquad \frac{1}{2\sigma^2} \sum_{t=1}^{T} \tilde{V}'_{nt}(\zeta)\tilde{V}_{nt}(\zeta),
\label{Eq: LogL2}
\end{eqnarray}
where
$
\tilde{V}_{nt}(\zeta) = R_n(\rho) [S_n(\lambda) \tilde{Y}_{nt} - \tilde{X}_{nt} \beta],
$
with
\begin{equation}
\tilde{Y}_{nt} = {Y}_{nt} - \sum_{t=1}^{T} Y_{nt}/T, \quad \tilde{X}_{nt} = X_{nt} - \sum_{t=1}^{T} X_{nt}/T. \label{Eq. transf}
\end{equation}
The MLE $\hat\theta_{n,T}$ for $\theta$ is
an $M$-estimator obtained by solving
$
\hat\theta_{n,T}=\text{arg}\max_{\theta \in \Theta} \ell_{n,T}(\theta).
$
It implies the system of estimating equations:
\begin{equation}
0=\frac{\partial\ell_{n,T}(\hat\theta_{n,T})}{\partial {\theta} } = \sum_{t=1}^{T} (T-1)^{-1}\psi_{nt}(\hat\theta_{n,T}),
\label{Eq: M_est}
\end{equation}
{ where $\psi_{nt}$
is the likelihood score function}
\begin{equation} \label{Eq: score_text}
\psi_{nt}(\theta) = \left( \begin{array}{cc}
\frac{T-1}{\sigma^2} (R_n(\rho) \tilde{X}_{nt})' \tilde{V}_{nt}(\zeta) \\
\frac{T-1}{\sigma^2}\left( \left(\ddot{G}_n\ddot{X}_{nt}\beta\right)' \tilde{V}_{nt}(\zeta) + \tilde{V}_{nt}^{'}\ddot{G}_{n}^{'}\tilde{V}_{nt}\right)-\frac{(T-1)^{2}}{T}\text{tr}({G}_n(\lambda))
\\
\frac{T-1}{\sigma^2} (H_n(\rho) \tilde{V}_{nt}(\zeta))' \tilde{V}_{nt}(\zeta) - \frac{(T-1)^{2}}{T}\text{tr}(H_n(\rho)) \\
\frac{T-1}{2\sigma^4} \left(\tilde{V}'_{nt}(\zeta)\tilde{V}_{nt}(\zeta) - \frac{n(T-1)}{T} \sigma^2 \right)
\end{array} \right),
\end{equation}
where $ G_n(\lambda) = W_n S_n^{-1}, \quad H_n(\rho) = M_n R_n^{-1}, \quad \ddot{G}_n(\lambda) = R_nG_nR_n^{-1}$, and $\ddot{X}_{nt} =R_n \tilde{X}_{nt}.$
\section{Methodology}
\label{Sec: Methodology}
We assume $ n \gg T$, so we deal with so-called \textit{micro panels}: in the econometric literature, this type of data typically involve annual records covering a short time span for each individual. Within this setting for $T$ being fixed,
the standard asymptotic arguments rely crucially on the number $n$
of individuals tending to infinity; see \citet{LY10}. In contrast, in our development,
we consider small sample cross-sectional asymptotics (\citet{FR90}), and we still leave $T$ fixed (possibly small). However,
we will keep $T$ in the notation of normalizing factors to demonstrate the
improved rate of convergence that would result if $T \to \infty$ or it is large. The derivation of our higher-order techniques relies on three steps: (i) defining a second-order asymptotic (von Mises) expansion for the MLE, see \S\ref{Sec: M_fun}; (ii) identifying the corresponding $U$-statistic, see \S\ref{Sec: IIordervM};
(iii) deriving the Edgeworth expansion for the $U$-statistic as in \citet{BGVZ86} and deriving the saddlepoint
density by means of the tilted-Edgeworth technique, see \S \ref{Sec: Ustat} and \S\ref{Sec: sadd}.
Similar approaches are available in the standard setting of i.i.d.\ random variables in \citet{ER86}, \citet{BNSC89}, and \citet{GR96}.
\subsection{The $M$-functional related to the MLE and its first-order asymptotics} \label{Sec: M_fun}
Let us first define the $M$-functional related to the MLE. To this end, we remark that the likelihood score function in (\ref{Eq: score_text}) is a vector in $\mathbb{R}^{d}$, and each $l$-th element of
this vector, for $l=1,...,d$, is a sum of $n$ terms. In what follows, for $i=1,...,n$, we denote by $\psi_{i,t,l}(\theta)$ the $i$-th term, at time $t$, of this sum for the $l$-th component of the score.
To specify $\psi_{i,t,l}(\theta)$, we set
$R_n(\rho) =\left (r_1^{'}(\rho), r_2^{'}(\rho), \cdots, r_n^{'}(\rho)\right)',$
$$\tilde{X}_{nt}=\left[\tilde{X}_{nt,1},\tilde{X}_{nt,2}, \cdots, \tilde{X}_{nt,k}\right],$$
$\tilde{V}_{nt}(\zeta) =\left (\tilde{v}_{1t}(\zeta), \tilde{v}_{2t}(\zeta), \cdots, \tilde{v}_{nt}(\zeta)\right)'$ and $H_n(\rho) =\left (h_1^{'}(\rho), h_2^{'}(\rho), \cdots, h_n^{'}(\rho)\right)'$, where $r_i(\rho)$ and $h_i(\rho)$ are the $i_{th}$ row of $R_n(\rho)$ and $H_n(\rho)$, $g_{ii}$ and $h_{ii}$ are $i_{th}$ element of the diagonal of $G_n(\lambda)$ and $H_n(\rho)$, respectively. Then, from (\ref{Eq: score_text}), it follows
\begin{equation}
\psi_{i,t}(\theta)=\left( \begin{array}{c} \psi_{i,t,1}(\theta),\\ \psi_{i,t,2}(\theta) \\ \vdots \\ \psi_{i,t,d}(\theta)
\end{array} \right)_{d\times 1}
= \left( \begin{array}{cc}
\frac{T-1}{\sigma^2} r_i(\rho)\tilde{X}_{nt,1}\tilde{v}_{it}(\zeta) \\
\frac{T-1}{\sigma^2} r_i(\rho)\tilde{X}_{nt,2}\tilde{v}_{it}(\zeta)\\
\vdots \\
\frac{T-1}{\sigma^2} r_i(\rho)\tilde{X}_{nt,k}\tilde{v}_{it}(\zeta) \\
\frac{T-1}{\sigma^2}r_i(\rho) \left(G_n\tilde{X}_{nt} \beta +G_nR_n^{-1}(\rho)\tilde{V}_{nt}(\zeta)\right) \tilde{v}_{it}(\zeta) - \frac{(T-1)^{2}}{T}g_{ii} \\
\frac{T-1}{\sigma^2} h_i(\rho)\tilde{V}_{nt}(\zeta) \tilde{v}_{it}(\zeta) - \frac{ (T-1)^{2}}{T}h_{ii}\\
\frac{T-1}{2\sigma^4} \left(\tilde{v}_{it}(\zeta)^2 - \frac{T-1}{T} \sigma^2 \right)
\end{array} \right)_{d\times 1} . \label{Eq: s_it}
\end{equation}
Thus, for every $t=1,2,...,T$, we have
$\displaystyle
\psi_{nt}(\theta)=\left(\sum_{i=1}^{n} \psi_{i,t,1}(\theta),..., \sum_{i=1}^{n} \psi_{i,t,d}(\theta)\right)',
$
and, from (\ref{Eq: M_est}), it follows that the MLE is the solution to
\begin{equation}
\frac{1}{n} \sum_{t=1}^{T} \left(\sum_{i=1}^{n} (T-1)^{-1}\psi_{i,t,1}(\hat\theta_{n,T}),..., \sum_{i=1}^{n} (T-1)^{-1}\psi_{i,t,d}(\hat\theta_{n,T})\right)' = 0. \label{Eq. MestPn}
\end{equation}
The $M$-functional $\vartheta$ related to the MLE
is implicitly defined as
the unique functional root of:
\begin{equation}
\mathbb{E}\left\{\sum_{t=1}^{T} \left({T-1}\right)^{-1} \psi_{nt}\left[\vartheta(P_{\theta_0})\right] \right\}=0, \label{Eq: Mfunctional}
\end{equation}
{or equivalently via the asymptotic maximization $\theta_0 = \text{arg}\max_{\theta \in \Theta} \mathbb{E}[\ell_{n,T}(\theta_0)]$; see e.g., \citet{L04}. In what follows, we write $\theta_0 = \vartheta(P_{\theta_0})$ to emphasize the dependence of the functional on the measure $P_{\theta_0}$.}
{The finite sample version of the $M$-functional in (\ref{Eq: Mfunctional}) is the $M$-estimator defined in (\ref{Eq. MestPn}), or equivalently via the finite sample maximization $\hat\theta_{n,T}=\text{arg}\max_{\theta \in \Theta} \ell_{n,T}(\theta)$.
In what follows, we write $\hat\theta_{n,T}= \vartheta(P_{n,T})$, where $P_{n,T}$ is the measure associated to the $n$-dimensional sample}. We can check the uniqueness of the M-estimator on a case-by-case basis, using Assumption A (see below) and working on the Gaussian log-likelihood. For instance, in the case of the SAR model, we can compute the second derivative of $\ell_{n,T}$ w.r.t. $\lambda$ and check that $\ell_{n,T}$ is a concave function, admitting a unique maximizer. Alternatively, we can solve the estimating equations implied by first-order conditions related to $\ell_{n,T}$ resorting on a one-step procedure and using for instance the GMM estimator (see \citet{LY10} and reference therein) as a preliminary estimator; for a book-length description of one-step procedure; see, e.g., \citet{vdW98} Ch. 5.
In what follows, we set
$m:=n(T-1)$, with $m\rightarrow\infty$, as $n\rightarrow\infty$. Then, we introduce
\textbf{Assumption A.}{\it
\begin{enumerate}[(i)]
\item The elements $\omega_{n,ij}$ of $W_n$ and the elements $m_{n,ij}$ of $M_n$
in (\ref{Eq: GeneralY}) are at most of order $\tilde{h}_n^{-1}$, denoted by $O(1/\tilde{h}_n)$, uniformly in all i,j, where the rate sequence $\{\tilde{h}_n\}$ is bounded, and $\tilde{h}_n$ is bounded away from zero for all $n$.
As a normalization, we have $\omega_{n,ii} = m_{n,ii}=0$, for all i.
\item $n$ diverges, while $T\geq 2$ and it is finite.
\item Assumptions 2-5 and Assumption 7 in \citet{LY10} are satisfied.
\item Denote $C_n = \ddot{G}_n-{n}^{-1}{tr(\ddot{G}_n})I_n$ and $D_n = H_n- {n}^{-1} {tr(H_n)}I_n$ where $\ddot{G}_n = R_nG_nR_n^{-1}$ and $H_n=M_nR_n^{-1}$. Then $C_n^s = C_n+C_n^{'}$ and $D_n^s= D_n+D_n^{'}$. The limit of ${n^{-2}}\left[ tr(C_n^sC_n^s)tr(D_n^sD_n^s)-tr^2(C_n^sD_n^s)\right]$ is strictly positive as
$n \rightarrow \infty$.
\end{enumerate}
}
Assumptions A$(i)$
characterizes the behavior of $W_n$ and $M_n$ in terms of $n$, and $W_n$ and $M_n$ are row-normalized. It means $\omega_{n,ij}= d_{ij}/\sum_{j=1}^{n}d_{ij}$, where $d_{ij}$ is the spatial distance of the $i-${th} and the $j-${th} units in some (characteristic) space. For each $i$, the weight $\omega_{n,ij}$ defines an average of neighboring values. In what follows, we consider spatial weight matrices (like, e.g., Rook and Queen) such that $\sum_{j=1}^n d_{ij}=O(\tilde{h}_n)$ uniformly in $i$ and the row-normalized weight matrix satisfies Assumption A$(i)$; see e.g. \citet{L04}. For instance, $W_n$ as Rook creates a square tessellation with $\tilde{h}_n=4$ for the inner fields on the chessboard, and $\tilde{h}_n=2$ and $\tilde{h}_n=3$ for the corner and border fields, respectively.
Assumption A$(ii)$ defines the asymptotic scheme of our theoretical development, in which
we consider $n$ cross-sectional units and we leave $T$ fixed. Assumption A$(iii)$ refers to \citet{LY10}, who develop the first-order asymptotic theory.
All $W_n$, $M_n$, $S_n^{-1}(\lambda)$, $R_n^{-1}(\rho)$ are uniformly bounded by Assumption A$(iv)$ which guarantees the convergence of the asymptotic variance, see below.
Assumption A$(iv)$ states the identification conditions of the model and the conditions for the nonsingularity of the limit of the information matrix. In particular, it implies that the $(d\times d)$-matrix
\begin{equation}
M_{i,T}(\psi,P_{\theta _{0}}) = \mathbb{E}\left[-(T-1)^{-1}\sum_{t=1}^T {\partial \psi_{i,t}(\theta)}/{\partial\theta}\Big\vert_{\theta = \theta_0}\right]
\label{Eq: Mit}
\end{equation}
is non-singular. {Under Assumption A$(i)$-A$(iv)$, Theorem 1 part(ii) in \citet{LY10} shows
that
$\displaystyle
\lim_{n \rightarrow \infty} \hat\theta_{n,T}
= \theta_0 $.} Furthermore, Theorem 2 point (ii) in \citet{LY10} implies, as $n \to \infty$, that
$\hat\theta_{n,T}$ satisfies
$
\sqrt{m} (\hat\theta_{n,T}- \theta_0) \overset{\mathcal{D}}{\rightarrow} \mathcal{N} \left( 0, \Sigma^{-1}_{0,T} \right),$ and $\Sigma_{0,T}
= \text{plim}_{n \rightarrow \infty} \Sigma_{0,n,T}
$. The operator $\text{plim}$ stands for the limit in probability and the expression of $ \Sigma_{0,n,T}$ is available in the online Supplementary Material (see Appendix B). The first-order asymptotics is obtained letting $n\rightarrow \infty$;
there is no need for $T\rightarrow \infty$ to obtain a consistent and asymptotically normal $M$-estimator.
\vspace{-0.4cm}
\subsection{Second-order von Mises expansion} \label{Sec: IIordervM}
To define a higher-order density approximation to the finite-sample density of the MLE, we need to derive its higher-order asymptotic expansion, making use of \\
\textbf{Assumption B.} {\it
\begin{enumerate}[(i)]
\item
${\partial^2 \psi_{i,t,l}(\theta)}/{\partial\theta \partial\theta'}$ exists at $\theta=\theta_0$, for every $i=1,..,n$, $t=1,...,T$ and $l=1,..,d$. \\
\item
The $d\times d$-matrix
$\mathbb{E} \left[(T-1)^{-1} \sum_{t=1}^{T} {\partial^2 \psi_{i,t,l}(\theta)}/{\partial\theta \partial\theta'}\Big\vert_
{\theta=\theta_0} \right ]$
is positive semi-definite, for every $l=1,..,d$.
\end{enumerate}
}
Then, we state the following
\begin{lemma} \label{Lemma_IIvM}
Let the MLE be defined as in (\ref{Eq: M_est}). Under Assumptions A-B, the following expansion holds:
\begin{equation} \label{IIorder_txt}
\vartheta(P_{n,T})-\vartheta(P_{\theta _{0}})=\frac{1}{n} \sum_{i=1}^n IF_{i,T} (\psi,P_{\theta _{0}}) +\frac{1}{2n^2} \sum_{i=1}^n \sum_{j=1}^n \varphi_{i,j,T}(\psi, P_{\theta _{0}}) + O_P(m^{-3/2}),
\end{equation}
where
\begin{equation} \label{IF_it}
IF_{i,T}(\psi,P_{\theta _{0}})=M_{i,T}^{-1}(\psi,P_{\theta _{0}}) (T-1)^{-1}\sum_{t=1}^T \psi_{i,t}(\theta_{0}),
\end{equation}
and
\begin{eqnarray} \label{kerneldue_txt}
\varphi_{i,j, T}(\psi,P_{\theta_0})
&=& IF_{i,T}(\psi, P_{\theta_0})+IF_{j,T}(\psi, P_{\theta_0})
+ M_{i,T}^{-1}(\psi, P_{\theta_0})\Gamma_{i,j,T}(\psi,P_{\theta_0}) \notag \\
&&+M_{i,T}^{-1}(\psi,P_{\theta _{0}})\left\{(T-1)^{-1}\sum_{t=1}^{T}\frac{\partial \psi_{j,t}(\theta)}{\partial\theta}\Big\vert_{\theta_0}IF_{i,T}(\psi,P_{\theta_0}) \notag \right. \\
&& \left.+ (T-1)^{-1}\sum_{t=1}^{T} \frac{\partial \psi_{i,t}(\theta)}{\partial\theta}\Big\vert_{\theta=\theta_0} IF_{j,T}(\psi, P_{\theta_0})\right\},
\end{eqnarray}
where
\begin{equation}
\Gamma_{i,j,T}(\psi,P_{\theta_0})'=\left(
\begin{array}{ccc}
IF'_{j,T}(\psi,P_{\theta_0})& \mathbb E\left[\sum_{t=1}^{T} \frac{\partial^2 \psi_{i,t,1}(\theta_0)}{\partial\theta \partial\theta'}\Big\vert_
{\theta=\theta_0} \right
]&
IF_{i,T}(\psi,P_{\theta_0}) \\
\vdots \\
IF'_{j,T}(\psi ,P_{\theta_0}) & \mathbb{E}\left[\sum_{t=1}^{T} \frac{\partial^2 \psi_{i,t,d}(\theta_0)}{\partial\theta \partial\theta'}
\Big\vert_{\theta=\theta_0} \right
] &
IF_{i,T}(\psi,P_{\theta_0})
\end{array}
\right)\ ,
\label{gammader}
\end{equation}
and $M_{i,T}(\psi, P_{\theta _{0}})$ is defined by (\ref{Eq: Mit}).
\end{lemma}
In (\ref{IIorder_txt}), we interpret the quantities $ IF_{i,T}(\psi, P_{\theta _{0}}) $, the first-order von Mises kernel, and $ \varphi_{i,j,T}(\psi, P_{\theta _{0}})$, the second-order von Mises kernel, as functional derivatives of the $M$-functional related to the MLE;
see \citet{Fernholz2001}. Specifically, the first term, of order $m^{-1} \propto n^{-1}$, is the Influence Function (IF) and represents the standard tool applied to derive the first-order (Gaussian) asymptotic theory of the MLE; see, e.g., \citet{vdW98} and \citet{BaltagiBook} for a book-length introduction. The second term in (\ref{IIorder_txt}), of order $m^{-2} \propto n^{-2}$, plays a pivotal role in our derivation of higher-order approximation.
\vspace{-0.5cm}
\subsection{Approximation via $U$-statistic} \label{Sec: Ustat}
The result of Lemma \ref{Lemma_IIvM} together with the chain rule define a second-order asymptotic expansion for a real-valued function of the MLE,
such as a component of $\vartheta(P_{n,T})$ or a linear contrast. In Lemma \ref{Lemma_IIvM_q}, we show that
we can write the asymptotic expansion in terms of a $U$-statistic of order two. To this end, we introduce the following assumption.
\textbf{Assumption C.} \\
{\it Let
$q$ be a function from $\mathbb{R}^{d}$ to $\mathbb{R}$, which has continuous and nonzero
gradient at $\theta = \theta_0$ and continuous second derivative at $\theta = \theta_0$. }
Then, we have
\begin{lemma} \label{Lemma_IIvM_q}
Under Assumptions A-C, the following expansion holds:
$$
q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})] = \frac{2}{n(n-1)}
\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} h_{i,j,T}\left(\psi,P_{\theta_0} \right) + O_P(m^{-3/2}),
$$
where
\begin{eqnarray} \label{Eq: h_ijT}
h_{i,j,T}\left(\psi,P_{\theta_0} \right) &=& g_{i,T}\left(\psi,P_{\theta_0} \right)+ g_{j,T}\left(\psi,P_{\theta_0} \right)+\gamma_{i,j,T}(\psi,P_{\theta_0}) \nonumber \\
&=& \frac{1}{2} \left\{ IF'_{i,T}(\psi,P_{\theta_0}) +IF'_{j,T}(\psi,P_{\theta_0}) + \varphi'_{i,j, T}(\psi,P_{\theta_0})\right\}\frac{\partial q(\vartheta)}{\partial\vartheta}\Big\vert
_{\theta=\theta_0} \nonumber \\
&+& \frac{1}{2} IF'_{i,T}(\psi,P_{\theta_0}) \frac{\partial^2 q(\vartheta)}{\partial\vartheta\partial\vartheta'}\Big\vert
_{\theta=\theta_0} IF_{j,T}(\psi,P_{\theta_0}) ,
\end{eqnarray}
with
\begin{equation} \label{Eq: g_iT}
g_{i,T}(\psi,P_{\theta_0}) = \frac{1}{2} \left( IF'_{i,T}(\psi,P_{\theta_0}) \frac{\partial q(\vartheta)}{\partial\vartheta}\Big\vert_{\theta=\theta_0} \right),
\end{equation}
\begin{equation}
\gamma_{i,j,T}(\psi,P_{\theta_0})=\frac{1}{2}
\left( \varphi'_{i,j, T}(\psi,P_{\theta_0})\frac{\partial q(\vartheta)}{\partial\vartheta}\Big\vert
_{\theta=\theta_0} + IF'_{i,T}(\psi,P_{\theta_0}) \frac{\partial^2 q(\vartheta)}{\partial\vartheta\partial\vartheta'}\Big\vert
_{\theta=\theta_0} IF_{j,T}(\psi,P_{\theta_0}) \right). \label{Eq: gamma_ijT}
\end{equation}
\end{lemma}
The function $q$ may select, e.g., a single component of the vector $\theta_0$. In many empirical applications, the most interesting parameter is the spatial correlation coefficient $\lambda_0$, and the null hypothesis is zero correlation versus the alternative hypothesis of positive spatial correlation---the aim being to check whether there is a contagion effect.
\vspace{-0.3cm}
\subsection{Higher-order asymptotics} \label{Sec: sadd}
Making use of Lemma \ref{Lemma_IIvM} and Lemma \ref{Lemma_IIvM_q}, we derive
the Edgeworth and the saddlepoint approximation
to the distribution of a real-valued function $q$ of the MLE.
Let $f_{n,T}(z)$ be the true density of $q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})]$ at the point
$z \in \mathcal{A}$, where $\mathcal{A}$ is a compact subset of $\mathbb{R}^d$.
Our derivation of the saddlepoint density approximation to $f_{n,T}(z)$ is based on the
tilted-Edgeworth expansion for $U$-statistics of order two. With this regard, a remark is in order.
From (\ref{Eq: s_it}), we see that
the terms in the random vector $\psi_{nt}(\theta_0)$ depend on the rows of the weight matrix $W_n(\rho)$ and $M_n(\lambda)$. As
a consequence, these terms are independent but not identically distributed random variables, and we need to derive
the Edgeworth expansion for our $U$-statistic taking into account this aspect.
To this end, we approximate
the cumulant generating function (c.g.f.) of our $U$-statistic by
\textit{summing (in $i$ and $j$) the (approximate) c.g.f.} of each $h_{i,j,T}$ kernel.
This is an extension of the derivation by
\citet{BGVZ86} for i.i.d. random variables.
To elaborate further, we introduce
\textbf{Assumption D.} \\ {\it
Suppose that there exist positive numbers $\delta$, $\delta_1$, $C$ and positive and continuous functions $\chi_{j}$: $(0,\infty)\to (0,\infty)$, $j=1,2,$ satisfying $\lim_{z\to\infty} \chi_1(z)=0$, $\lim_{z\to\infty} \chi_2(z)\geq \delta_1 >0$, and a real number $\alpha$ such that $\alpha \geq 2+\delta>2$,
\begin{enumerate}[(i)]
\item $ \mathbb{E}\left[\vert \gamma_{i,j,T}(\psi,P_{\theta_0})\vert^{\alpha}\right]<C$ for any $i$ and $j$, $1\leq i< j\leq n$,
\item $\mathbb{E}\left[ g_{i,T}(\psi,P_{\theta_0})^{4}1_{[z,\infty)}(\vert g_{i,T}(\psi,P_{\theta_0})\vert)\right]<\chi_1(z)$ for all $z>0$ and any $i$, $1\leq i\leq n$,
\item $\Big\vert\mathbb{E}\left[e^{\iota \nu g_{i,T}(\psi,P_{\theta_0})}\right]\Big\vert\leq 1-\chi_2(z) <1$ for all $z>0$ and any $i$, $1\leq i\leq n$ and $\iota^2=-1$,
\item $\vert \vert M_{i,T}(\psi,P_{\theta _{0}})- M_{j,T}(\psi,P_{\theta _{0}}) \vert \vert=O(n^{-1})$ uniformly in $\lambda$ and $\rho$.
\end{enumerate} }
A few comments are in order. Assumptions D$(i)$-$(iii)$ are similar to the technical assumptions in \citet{BGVZ86} p.
1465 and 1477. However, there are some differences between our assumptions and theirs. Indeed, to
take into account the non identical distribution of $\psi_{i,t}$ and $\psi_{j,t}$, for $i \neq j$,
we consider the first- and second-order von Mises kernels for each $i$ (as in D$(i)$-$(iii)$).
It is different from \citet{BGVZ86}: compare, e.g., our D$(ii)$ to their Eq. (1.17).
D$(iv)$ is not considered in \citet{BGVZ86}: it is a peculiar assumption needed for our higher-order asymptotics (the technical aspects are available in
Lemma A.1 and its proof in Appendix A). {In Appendix D.2,
we illustrate that, in the case of the SAR(1) model, the validity of D($iv$) is related to more primitive expressions involving the entries of (some powers of) $W_n$.} For other models, one should derive such primitive expressions on a case-by-case basis. For the sake of generality, here we provide an intuition on D($iv$). Let us consider two different locations $i$ and $j$.
From (\ref{Eq: Mit}), we see that D$(iv)$
imposes a structure on the information available at different locations. Indeed, $ M_{i,T}(\psi,P_{\theta _{0}})$ and $ M_{j,T}(\psi,P_{\theta _{0}})$ contribute to the asymptotic variance of the MLE. Since $ M_{i,T}(\psi,P_{\theta _{0}}) $ is related to the information available at the $i$-th location, D$(iv)$ essentially assumes that
there exists an informative content which is common to location $i$ and $j$, whilst the ({Frobenious} norm of the) information content specific to each location is of order $O(n^{-1})$.
\color{black}
\begin{proposition} \label{Prop_EDG}
Under Assumptions A-D, the Edgeworth expansion $\Lambda_m(z)$ for the c.d.f. $F_m$ of $\sigma_{n,T}^{-1}\{q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})]\}$ is
\begin{eqnarray}
\Lambda_m(z) &=& \Phi(z) - \phi(z)\left\{n^{-1/2} \frac{\kappa_{n,T}^{(3)}}{3!} (z^2-1) + n^{-1} \frac{\kappa_{n,T}^{(4)}}{4!} (z^3-3z) + n^{-1} \frac{\kappa_{n,T}^{(3)}}{72} (z^5-10z^2+15z) \right\} \nonumber \\% \notag \\
\label{Eq: EDG_CDF}
\end{eqnarray}
where $z \in \mathcal{A}$, $\sigma_{n,T}$ is the standard deviation of $q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})]$, $\Phi(z)$ and $\phi(z)$ are the c.d.f. and p.d.f of a standard normal r.v. respectively, $\kappa_{n,T}^{(3)}n^{-1/2}$ and $\kappa_{n,T}^{(4)}n^{-1}$ are the approximate third and fourth cumulants of $\sigma_{n,T}^{-1}\{q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})]\}$, as defined in (A.15) and (A.18), respectively.
Then
$\sup_z|F_m(z)-\Lambda_m(z)| =o(m^{-1}).$
\end{proposition}
In addition, we can get the saddlepoint density approximation by exponentially tilting the Edgeworth expansion.
\begin{proposition} \label{Prop_SAD}
Under Assumption A-D, the saddlepoint density approximation to the density of $q[\vartheta(P_{n,T})]-q[\vartheta(P_{\theta_0})]$ at the point $z \in \mathcal{A}$ is
\begin{equation}
p_{n,T}(z)=\left[ \frac{n}{2\pi \tilde{\mathcal{K}}_{n,T}^{\prime \prime }(\nu)}
\right] ^{1/2}\exp \left\{ n \left[\tilde{\mathcal{K}}_{n,T}(\nu)-\nu z\right]\right\},
\label{Eq: SAD_panel}
\end{equation}
with relative error of order $O(m^{-1})$, $\nu:=\nu(z)$ is the saddlepoint defined by
\begin{equation}
\tilde{\mathcal{K}}_{n,T}^{\prime}(\nu)=z,
\label{Eq. sadd}
\end{equation}
the function $\tilde{\mathcal{K}}_{n,T}$ is the approximate c.g.f. of
$\sqrt{n}(q[\vartheta(P_{n,T})]- q[\vartheta(P_{\theta_0})]) $, as defined in (A.42),
while $\tilde{\mathcal{K}}_{n,T}^{\prime}$ and $\tilde{\mathcal{K}}_{n,T}^{\prime \prime}$ represent the first and second derivative of $\tilde{\mathcal{K}}_{n,T}$, respectively. Moreover,
\begin{eqnarray}
P\left\{ q[\vartheta(P_{n,T})] - q[\vartheta(P_{\theta_0})] > z \right\} & = & \left[ 1- \Phi(r) + \phi(r) \left( \frac{1}{c}- \frac{1}{r} \right) \right] \left[1+O(m^{-1})\right],
\label{LR_pvalue}
\end{eqnarray}
$$c=\nu \left[ \tilde{\mathcal{K}}_{n,T}^{\prime \prime}(\nu )\right] ^{1/2} \quad \text{and} \quad
r=\text{sgn}(\nu )\left\{ 2n\left[ \nu z-\tilde{\mathcal{K}}_{n,T}(\nu )\right] \right\}^{1/2}.$$
\end{proposition}
The proofs of these Propositions are available in Appendix A.
They rely on a argument similar to the one applied in the proof of \citet{F82} for the derivation of a saddlepoint
density approximation of multivariate $M$-estimators,
and in \citet{GR96} for general statistics. Following \citet{D80}, we can further normalize
$p_{n,T}$ to obtain a proper density by dividing the right hand side of (\ref{Eq: SAD_panel}) by its integral with respect to $z$. This normalization typically improves
even further the accuracy of the approximation.
\subsection{Links with the econometric literature}
\label{Sec: Links_econometrics}
The expansions in Proposition \ref{Prop_EDG} and Proposition \ref{Prop_SAD} are connected with the results on higher-order expansions available in the spatial econometric literature, as cited in \S \ref{Sec: Intro}. However, some key differences are worth a mention.
(i) The Edgeworth expansion in \citet{RR15} is for the concentrated MLE of $\lambda$ and it is based on a higher-order Taylor expansion of the concentrated likelihood score; see also \citet{RR14_EJ,RR14_ET, RR15}, \citet{HM18}, and \citet{HM20}. In contrast, our method is based on a von Mises expansion of the MLE functional of the whole model parameter and we resort on a marginalization procedure to obtain the saddlepoint density approximation of the parameter(s) of interest.
Therefore, Lemma \ref{Lemma_IIvM} and Lemma \ref{Lemma_IIvM_q} give generality and flexibility to our approach: not only we may focus on $\lambda$, but also, e.g., on $\rho$ (which contains information on the spatial dependence of the innovation terms) and/or on $\beta$ (which convey information on the significance of the time-varying covariates).
(ii) Our saddlepoint approximation is more general than the Edgeworth-based approximations available in the econometric literature, since we work with a larger class of models, which includes the model in \citet{RR15} as a special case.
(iii) Although the inference (e.g., testing)
derived using the Edgeworth expansion improves on the standard first-order asymptotics,
it is well-known (see, e.g., \citet{FR90}) that, in general, this technique provides a
good approximation in the center of the distribution, but can be inaccurate in the tails, where the Edgeworth expansion
can even become negative. It can lead to inaccurate approximations. Our saddlepoint approximation
is a density-like object and is always nonnegative.
(iv) Our saddlepoint approximation yields a tail-area approximation via a Lugannani-Rice type formula. A similar result is not available for the Edgeworth expansion of the concentrated MLE derived in \citet{RR15}. Recently, \citet{HM20} studied the adjusted profile likelihood estimation method and obtained a result similar to our tail-area approximation. Their formula is derived for the spatial autoregressive model with covariates. However, they do not prove the higher-order properties of their approximation.
In Proposition \ref{Prop_EDG}, we prove that our saddlepoint density approximation
features \textit{relative error} of order $O(1/(n(T-1)))$. This has to be contrasted with the extant Edgeworth expansion, which entails an absolute error of lower order---more precisely, the error order is $o((nT)^{-1/2})$, when the entries of the spatial matrix are $O(1)$; see Eq. (2.15) in \citet{RR15}. Achieving a small relative error is appealing in
tail areas where the probabilities are small.
(v) In the comparison with the bootstrap, our methodology
does not need resampling. Moreover, it does neither require
bias correction, nor any studentization.
\section{Computational aspects} \label{Sec: algm}
Most of the quantities related to the saddlepoint density approximation $p_{n,T}$ and the tail area in (\ref{LR_pvalue}) are available in closed-form.
In Appendix C.1, we provide an algorithm (see Algorithm 1) in which we itemize the main computational steps needed to implement the saddlepoint tail area approximation,
for a given transformation $q$
and for a given reference parameter $\theta_0$---it is, e.g., the parameter characterizing the null hypothesis in a simple hypothesis testing, where the tail area
probability is an approximate $p$-value.
\vspace{-0.5cm}
\section{Comparisons with other approximations and testing in the presence of nuisance parameters} \label{Sec: MC}
\vspace{-0.25cm}
We compare the performance of our saddlepoint approximations to other routinely-applied asymptotic techniques. To start with, we consider the SAR(1) model, where $\lambda$ is the only unknown parameter. Then, we move to the SARAR(1,1) model, where we illustrate how to take care of nuisance parameters. We use the same setting as in \S \ref{Sec: Motivation}; we refer to the online Supplementary Material (Appendix D) for more details and for additional results.
\subsection{Comparisons with other asymptotic techniques} \label{Sec: AsyvsSad}
\textit{Saddlepoint vs first-order asymptotics.} For the SAR(1) model, we analyse the behaviour of the MLE of $\lambda_0$,
whose PP-plots are available in Figure \ref{Fig: PP}. For each type of
$W_n$, for $n=24$ and $n=100$, the plots show that the saddlepoint approximation is closer to the ``exact'' probability than the first-order asymptotics approximation.
For $W_n$ Rook, the saddlepoint approximation improves on the routinely-applied first-order asymptotics. In Figure
\ref{Fig: PP}, the accuracy gains are evident also for $W_n$ Queen and Queen with torus, where the first-order asymptotic theory displays large errors essentially over the whole support (specially in the tails). On the contrary,
the saddlepoint approximation is close to the 45-degree line.
\begin{figure}[hbt!]
\begin{center}
\begin{tabular}{ccc}
Rook & Queen & Queen torus \\
\begin{turn}
{90} \hspace{2.5cm} $n = 24$
\end{turn}
\includegraphics[width=0.3\textwidth, height=0.265\textheight]{PPplot_n24WnRook.eps} &
\includegraphics[width=0.3\textwidth, height=0.265\textheight]{PPplot_n24WnQueen} &
\includegraphics[width=0.3\textwidth, height=0.265\textheight]{PPplot_n24WnQueenTorus.eps}
\\
\end{tabular}
\caption{SAR(1) model: PP-plots for saddlepoint (continuous line) vs asymptotic normal (dotted line) probability approximation, for the MLE $\hat\lambda$, for $n=24$ and $n=100$, $\lambda_0=0.2$, and different $W_n$.}
\label{Fig: PP}
\end{center}
\end{figure}
\textit{Saddlepoint vs Edgeworth expansion (testing simple hypotheses).} The Edgeworth expansion derived in Proposition \ref{Prop_EDG} represents the natural
alternative to the saddlepoint approximation since it is fully analytic.
To gain insights into the different behavior of the saddlepoint and Edgeworth approximations,
we investigate the size of a hypothesis test based on the approximations.
We set $n=24$ and we assume that $\sigma^2$ is known and equal to one.
We consider the simple null hypothesis $H_0$: $\lambda_0=0$ for a one-sided test of zero against positive values of
spatial correlation. We use 25,000 replications of $\hat\lambda_{n,T}$ to get the empirical estimate $\hat F_0$ of the c.d.f. $F_0$ of the estimator under the null hypothesis. We use the generic notation $G$ for the c.d.f. of one of the Edgeworth, or saddlepoint approximations, under the null hypothesis. For the sake of completeness, we also display the results for the Gaussian (first-order) approximation. The empirical rejection probabilities $\hat \alpha = 1-\hat F_0(G^{-1}(1-\alpha))$ are shown in Figure \ref{Fig: Size} for nominal size $\alpha$ ranging from 1\% to 10\%, and correspond to an estimated size. We have overrejection when we are above the 45-degree line. We observe strong size distortions for the asymptotic and Edgeworth approximations as expected from the previous results.
The saddlepoint approximation exhibits only mild size distortions. For example, we get an estimated size $\hat \alpha$ of 11.72\%, 7.36\%, 5.70\%, for the Normal, Edgeworth, and saddlepoint approximations, for a nominal size of 5\%.
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.5\textwidth, height=0.35\textheight]{Size_n24WnRook.eps}
\caption{SAR(1) model: Estimated $\hat \alpha$ versus nominal size $\alpha$ between 1\% and 10\% under saddlepoint (continuous line), Edgeworth (dotted line with diamonds) and first-order asymptotic approximation (dotted line). $W_n$ is Rook, $n=24$ and $\lambda_0=0.0$.}
\label{Fig: Size}
\end{center}
\end{figure}
\textit{Saddlepoint vs parametric bootstrap.}
The parametric bootstrap represents a (computer-based) competitor, commonly applied in statistics and econometrics.
To compare our saddlepoint approximation to the one obtained by bootstrap, we consider different numbers of bootstrap repetitions, labeled as $B$: we use $B=499$ and $B=999$. For space constraints, in Figure \ref{Fig: FB}, we display the results for $B=499$ (similar plots are available for $B=999$) showing the functional boxplots (as obtained iterating the procedure 100 times) of the bootstrap approximated density, for sample
size $n =24$ and for $W_n$ is Queen.
\begin{figure}[hbt!]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.485\textwidth, height=0.35\textheight]{FB1_Densities_n24WnQueen_50000.eps} &
\includegraphics[width=0.485\textwidth, height=0.35\textheight]{FB1_Densities_n24WnQueen_50000_RightTail.eps}
\end{tabular}
\caption{SAR(1) model. Left panel: Density plots for saddlepoint (continuous line) vs the functional boxplot of the parametric bootstrap probability approximation to the exact density (as expressed by the histogram and obtained using
MC with size 25000), for the MLE $\hat\lambda$ and $W_n$ is Queen. Sample size is $n=24$, while $\lambda_0=0.2$. Right panel: zoom on the right tail. In each plot, we display the functional central curve (dotted line with crosses), the $1_{st}$ and $3_{rd}$ functional quartile (two-dash lines).
}
\label{Fig: FB}
\end{center}
\end{figure}
To visualize the variability entailed by the bootstrap, we display the first and third quartile curves (two-dash lines) and the median functional curve (dotted line with crosses); for details about functional boxplots, we refer to \citet{Sun11} and to R routine \texttt{fbplot}. We notice that, while the bootstrap median functional curve (representing a typical bootstrap density approximation) is close to the actual density (as represented by the histogram), the range between the quartile curves illustrates that the bootstrap approximation has a variability.
Clearly, the variability depends on $B$: the larger is $B$, the smaller is the variability. However, larger values of $B$ entail bigger computational costs: when $B=499$, the bootstrap is almost as fast as the saddlepoint density approximation ({computation time
about 7 minutes, on a 2.3 GHz Intel Core i5 processor}),
but for $B=999$, it is three times slower. We refer to Appendix D.5 for additional numerical results.
\subsection{Testing in the presence of nuisance parameters}
\label{Sec: test_nuisance_parameters}
\subsubsection{Saddlepoint test for composite hypotheses} \label{Sec: testcomp}
Our saddlepoint density and/or tail approximations are helpful for testing simple hypotheses about $\theta_0$; see \S \ref{Sec: AsyvsSad}.
Another interesting case suggested by the Associate Editor and an anonymous referee that has a strong practical relevance is related to testing a
composite null hypothesis. It is a problem which is different from the one considered so far in the paper, because it raises the issue of dealing with nuisance parameters.
To tackle this problem, several possibilities are available.
For instance, we may
fix the nuisance parameters at the MLE estimates. Alternatively, we may consider to use
the (re-centered) profile estimators, as suggested, e.g., in \citet{HM18} and \citet{HM20}.
Combined with the saddlepoint density in (\ref{Eq: SAD_panel}), these techniques yield a ready solution to the nuisance parameter problem. In our numerical experience (see Appendix D.6 for an experiment about the SAR(1)), these solutions may preserve
reasonable accuracy in some cases.
Nevertheless, the main theoretical drawback related
to the use of MLE values for the nuisance parameter(s) is that it would not guarantee that
the second-order properties derived in the previous sections still hold.
To cope with this issue, we propose to build on
\citet{RRY03}, who derive a saddlepoint test statistic which
takes into account explicitly the nuisance parameters, while preserving relative error in normal region. We feel this test statistic represents the natural candidate within our setting: it shares the same spirit as our saddlepoint density approximation and it is derived going through steps which are similar to ours.
The paper by \citet{RRY03} defines the test statistic in the i.i.d.\ setting, while \citet{Lo09} and \citet{CR10} extend it to the non-i.i.d.\ data setting.
Let us consider a SARAR model whose parameter is
$\theta = (\theta_{10}, \theta_2)'$,
where $\theta_{10}$ is specified by the null composite hypothesis: typically, the null concerns $\lambda$ only, while $\theta_2$ contains all the nuisance parameters. More specifically, the parameter is
$\theta = (\lambda, \beta, \rho, \sigma^2)'$ and the general function
$q(\theta)$ used in the previous sections is simply
$q(\theta) = \lambda$. Thus, we have the composite hypothesis:
\begin{equation}
\mathcal{H}_0: \lambda=\lambda_0 = 0 \quad \text{vs} \quad \mathcal{H}_1: \lambda > 0, \label{TestComp}
\end{equation}
where $\theta= (\lambda,\theta_2)'$, with $\theta_{10}=\lambda_0$ and $\theta_2=(\beta, \rho, \sigma^2)'$. Then, we define the test statistic
\begin{equation}
{SAD}_n(\hat{\lambda}) = 2n \
\underset{\theta_2}{\inf} \ \underset{\nu}{\sup} \ -\mathcal{K}_\psi(\nu,\hat\lambda,\theta_2). \label{Eq. h}
\end{equation}
The function $\mathcal{K}_\psi(\nu,(\lambda,\theta_2))$ is the c.g.f. of the estimating function:
\begin{equation}
\mathcal{K}_\psi(\nu,\hat\lambda,\theta_2) = n^{-1}\sum_{i = 1}^{n}\ln E_{P_{(\lambda_{0},\theta_2)}} \exp(\nu^{T}\psi_{i}^{(T)}(\hat{\lambda},\theta_2)), \label{Eq. cgfpsi}
\end{equation}
where $\psi_{i}^{(T)}(\lambda,\theta_2):=\sum_{t=1}^{T}(T-1)^{-1}\psi_{i,t}(\lambda,\theta_2)$ and $\psi_{i,t}$ is as in (4.1). The c.g.f. $\mathcal{K}_\psi$ has a role analogous to the one of the c.g.f. of the $U$-statistic, that we derived in \S \ref{Sec: Methodology}. We highlight that the expected value in (\ref{Eq. cgfpsi}) is taken w.r.t. the probability $P_{(\lambda_{0},\theta_2)}$, where $\lambda_0$ is specified by the null, while the nuisance parameters are not fixed: the infimum over $\theta_2$ takes care of the nuisance parameters.
In our inference procedure, we have that
$\hat{\theta}_{n,T}=(\hat\lambda,\hat\theta_2)'$ is the solution to
$
\sum_{i=1}^{n} \psi_{i}^{(T)}(\lambda,\theta_2) = 0.
$
Under the null hypothesis, the test statistic ${SAD}_n(\hat{\lambda}) $ is asymptotically $\chi_1^2$ distributed with a relative error of order $O(m^{-1})$ in the normal region.
\subsubsection{Implementation aspects}
To implement the test (\ref{Eq. h}) for the problem (\ref{TestComp}), in Appendix C.2, we propose an algorithm (see Algorithm 2) and itemize the main steps needed to compute the test statistic.
\subsubsection{Numerical results}
Let us work with a SARAR(1,1) model, having no covariates and known variance $\sigma^2=1$ and $n=24$. It implies that
$\theta=(\lambda,\rho)'$ and we consider the problem in (\ref{TestComp}), with $\rho$ being the nuisance
parameter. We set three different values $\rho=0.25, 0.5, 0.75$ to analyze numerically the impact that
the spatial dependence in the innovation term has on $SAD_n$. We study the behaviour of the Wald
test, as obtained using the first-order asymptotic theory and making use of the expression
of the asymptotic variance as available in Appendix B. We compare the Wald test to $SAD_n$--to implement (\ref{Eq. h}) we make use of the R routine \texttt{nlm}. We consider
two types of spatial matrix $W_n$, the Rook and the Queen, and we set $W_n\equiv M_n$. Both test statistics are asymptotically $\chi_1^2$ distributed
under the null hypothesis. To compare them in small samples, we first obtain the 95th and 97.5th quantile of each test statistic; then we compute the corresponding probability as obtained using the $\chi_1^2$. We display the results in Table \ref{TabNuisance}. We see that the Wald test has severe size distortion. For instance, for $\rho=0.25$, we observe a relative error of about $30\%$, for the quantile of $95\%$, when $W_n$ is Rook, while the saddlepoint test entails a relative error of about $1.8\%$. Looking at the performance of ${SAD}_n$, we see that it is uniformly more accurate than the Wald test: considering all cases, we observe a maximal relative error of about $2\%$, for the quantile of $95\%$, when $\rho=0.75$ and $W_n$ is Queen; in the same setting, the Wald test entails a relative error of about $24\%$. Moreover, the size is fairly constant for the different values of $\rho$: it illustrates that the test statistic takes care correctly of the nuisance parameter.
\begin{table}[h!]
\begin{center}
\begin{tabular}{ccccccccccccccccccccccc}
\toprule
& & \multicolumn{3}{l}{$\rho=0.25$} & & \multicolumn{3}{l}{$\rho=0.5$} & & \multicolumn{3}{l}{$\rho=0.75$} \\
& & $95.00\%$ & $97.50\%$ & & & $95.00\%$ & $97.50\%$ & & & $95.00\%$ & $97.50\%$\\
\cmidrule{3-5} \cmidrule{7-9} \cmidrule{11-13}\\
Rook & \\
& \text{Wald} & $66.08\%$ & $89.89\%$ && & $98.33\%$ & $99.41\%$ && & $99.99\%$ & $99.99\%$ && \\
& ${SAD}_n$ & $96.71\%$ & $97.18\%$ && & $96.66\%$ & $97.18\%$ && & $95.55\%$ & $96.04\%$ &&\\
Queen & \\
& \text{Wald} & $72.71\%$ & $80.20\%$ && & $90.48\%$ & $98.00\%$ && & $98.36\%$ & $99.22\%$ && \\
& ${SAD}_n$ & $94.79\%$ & $96.94\%$ && & $95.90\%$ & $98.20\%$ && & $96.94\%$ & $97.49\%$ && \\
\bottomrule
\end{tabular}
\end{center}
\caption{Wald and $SAD_n$ test for the problem (\ref{TestComp})--spatial dependence in the presence of a nuisance parameter, in a SARAR(1,1) model--with no covariates and known variance. The quantiles are obtained using 100 repetitions of each test statistic.}
\label{TabNuisance}
\end{table}
\section{Empirical application} \label{Sec: Real}
\citet{FH80} document empirically that domestic saving rate in a country has a positive correlation with the domestic investment rate. It contrasts with the understanding that, if capital is perfectly mobile between countries, most of any incremental saving is invested to get the highest return regardless of any locations, and that such correlation should actually vanish. \citet{DE10} suggest to use spatial modeling since several papers challenge these findings but under the strong assumption that investment rates are independent across countries.
Such an assumption might influence the conclusions of applied statial economics.
In this empirical exercise, we investigate the presence of spatial autocorrelation in the investment-saving relationship. We consider investment and saving rates for 24 OECD countries between 1960 and 2000 (41 years). Because of macroeconomic reasons (deregulating financial markets), we divide the whole period into shorter sub-periods: 1960-1970, 1971-1985 and 1986-2000, as advocated by \citet{DE10}.
Since the cross-sectional size is only $n=24$, the asymptotics may suffer from size distortion as documented in \S\ref{Sec: MC}.
Therefore, we resort on a saddlepoint test to investigate whether or not there are inferential issues (coming from finite sample distortions and nuisance parameters) in the use of the first-order asymptotic theory. In line with the econometric literature, we specify the following SARAR(1,1) model for the three sub-periods:
\begin{equation}
\begin{aligned}
\text{ Inv}_{nt} &= \lambda_0 W_n \text{ Inv}_{nt} + \beta_0 \text{Sav}_{nt} + c_{n0} + E_{nt},& \\
E_{nt} &=\rho_0 M_n E_{nt} +V_{nt}, & \quad t=1,2, \cdots, T &
\end{aligned}
\label{Eq: Inv}
\end{equation}
where $\text{Inv}_{nt}$ is the $n \times 1$ vector of investment rates for all countries and $\text{Sav}_{nt}$ is the $n \times 1$ vector of saving rates. Each element $v_{it}$ in $V_{nt}$ is i.i.d across $i$ and $t$, having Gaussian distribution with zero mean and variance $\sigma_0^2$. $c_{n0}$ is the vector of fixed effects.
We assume $W_n=M_n$ and adopt two different weight matrices as in \citet{DE10}. The first one is based on the inverse distance. Each element $\omega_{ij}$ in $W_n $ is $d_{ij}^{-1}$, where $d_{ij}$ is the arc distance between capitals of countries $i$ and $j$. The second is the binary seven nearest neighbors (7NN) weight matrix. More precisely, $\omega_{ij}$=1, if $d_{ij} \leq d_{i}$ and $i \neq j$. Otherwise, $\omega_{ij}=0$, where $d_i$ is the $7_{th}$ order smallest arc-distance between countries $i$ and $j$ such that each country $i$ has exactly 7 neighbors. Both weight matrices are row-normalized.
We estimate the parameters using the MLE described in \S\ref{Sec: Setting}. Table \ref{Table: Estimates} gathers the point estimates (and their standard errors) that agree with the magnitudes found by \citet{DE10}. To investigate the validity of the model (\ref{Eq: Inv}), we test for spatial dependence, working on $\lambda=0$ and/or $\rho=0$. Specifically, our aim is to detect if and in which period(s) the inference yielded by the first-order asymptotic theory differs from the inference obtained using our saddlepoint test. With this goal, in
Table \ref{Table: Pvalues} we provide the $p$-values for testing (at the $5\%$ level) three different composite hypotheses: in the first row, we consider the problem of testing for $\lambda=0$; in the second row, we test for $\rho=0$; in the third row, we test for $\lambda=\rho=0$. To perform the tests, we consider the routinely-applied Wald test (as obtained using the first-order asymptotic approximation, ASY) and the saddlepoint test ($SAD_n$}). In each testing procedure, we treat the parameters not specified by the null hypothesis as nuisance parameters. In the $SAD_n$ test, we take care of the nuisance as indicated in (\ref{Eq. h}), while in the ASY test we simply plug-in the MLE estimates for the nuisance parameters---as it is customary in the econometric software based on the first-order asymptotic theory.
In the period 60-70, both ASY and $SAD_n$ yield the same inference, for both the considered types of weight matrix, with conventional significance levels.
The other sub-periods display some discrepancies between the inference obtained via ASY and via $SAD_n$. We do not want to discuss all discrepancies but only briefly comment on some key differences---we highlights the corresponding values in Table \ref{Table: Pvalues}. In the sub-period 71-85 under 7NN $W_n$, the saddlepoint test finds no evidence against no spatial dependence in the investing rates across countries, and vice-versa for the asymptotic approximation. Moreover, the ASY test does not find evidence against $\rho=0$, while the $SAD_n$ test rejects this composite hypothesis. Thus, the $SAD_n$ test indicates a spillover through the contemporary shocks between countries. This spillover goes through the innovations, i.e., through the unexpected part of the model dynamics, a finding not documentable when one relies on the first-order asymptotic theory. This results suggests that a test statistic designed to perform well in small samples and in the presence of nuisance parameters is able to document spatial dependence in the disturbances $E_{nt}$. Some differences are detectable also in the sub-period 86-00, under the inverse distance matrix.
\begin{table}[htb!]
\begin{center}
\begin{tabular}{cccccccc}
\toprule
& \multicolumn{3}{l}{Weight matrix: inverse distance} & & \multicolumn{3}{l}{Weight matrix: 7 nearest neighbours} \\
\cmidrule{2-4} \cmidrule{6-8}
&1960-1970& 1971-1985& 1986-2000& & 1960-1970& 1971-1985& 1986-2000\\
\midrule
$\beta$ & 0.935(0.05)
&0.638(0.04)
& 0.356(0.07)
& & 0.932(0.05) &0.633(0.04) & 0.368(0.07) \\
$\lambda$ & 0.004(0.10)&0.381(0.11)& 0.430(0.30)
& & -0.016(0.09)& 0.340(0.10) &0.437(0.18) \\
$\rho$ &-0.305(0.22)& 0.334(0.16)& 0.222(0.40)
& & -0.219(0.19)&0.258(0.15) &0.025(0.28)
\\
\bottomrule
\end{tabular}
\caption{SARAR(1,1) model: Maximum likelihood estimates of Parameters $\beta$, $\lambda$, $\rho$. Standard errors are between brackets.}
\label{Table: Estimates}
\end{center}
\end{table}
\begin{table}[hbt!]
\begin{center}
\begin{tabular}{ccccccccc}
\toprule
& &\multicolumn{3}{l}{Weight matrix: inverse distance} & & \multicolumn{3}{l}{Weight matrix: 7 nearest neighbours} \\
\cmidrule{3-5} \cmidrule{7-9}
& &1960-1970& 1971-1985& 1986-2000& & 1960-1970& 1971-1985& 1986-2000\\
\midrule
\multirow {2}{*}{$\lambda=0$}&$SAD_n$ & 1.0000 &0.0096 &\textbf{0.0000} & & 0.9998& \textbf{0.2248}& 0.0000 \\
&ASY &1.0000 &0.0116 &\underline{0.5679}& &0.9987 & \underline{0.0130} &0.1123\\
\multirow {2}{*}{$\rho = 0$}&$SAD_n$ & 0.1134 &0.0024 &0.1217
& & 0.3232& \textbf{0.0403}& 0.9993 \\
&ASY&0.5890 &0.2261& 0.9578 & & 0.7101& \underline{0.3898} &0.9998\\
\multirow {2}{*}{$\lambda =\rho=0$}&$SAD_n$ &0.1414& 0.0000& 0.0000
& & 0.2603& 0.0000 &0.0000
\\
&ASY &0.4615& 0.0000 &0.0000& &0.5042& 0.0000 &0.0000\\
\bottomrule
\end{tabular}
\caption{SARAR(1,1) model: $p$-values of Saddlepoint ($SAD_n$) and Wald (ASY) tests for several composite hypotheses.}
\label{Table: Pvalues}
\end{center}
\end{table}
\begin{center}
{\large\bf SUPPLEMENTARY MATERIAL}
\end{center}
The online supplementary material includes proofs,
lengthy analytical derivations and additional numerical results for the SAR(1) model. All the codes and data are available in our Github repository.
{\small
\bibliographystyle{myapalike1}
\bibliography{biblio_Panel}}