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.
68,615 characters · 17 sections · 84 citation commands
Saddlepoint approximations for spatial panel data models
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Higher-order asymptotics, investment-saving, random field, tail area.
\spacingset{1.5} {.5in}
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., GG10, R12, C15, CW15, and recently CMW19 for book-length discussions in the statistical literature. In the econometric literature, see, e.g., KKP07, LY10, RR14_ET, RR15, and, for book-length presentations, BaltagiBook, An13, and 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 O75. The derivation of the first-order asymptotics is available in the econometric literature; we refer to the seminal paper by 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 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. 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. B13 and Ya05 extend these approximations to include also exogenous explanatory variables, which remain valid also when the process is not Gaussian. 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, 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. HM18 (see their \S 6) and 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 Hall92 and 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 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 D54 and the book-length presentations of FR90, Jensen95, K06, and brazzaleetal2007. For a result about testing on spatial dependence, see T02, and for developments in time series models, see 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. 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), we provide a motivating example. \S(ref) defines the general model setting and the estimation method, whereas the detailed methodology is presented in \S(ref). In particular in \S(ref) we provide a detailed discussion about the connections with the econometric literature. Algorithms and computational aspects are discussed in \S(ref). \S(ref) provides numerical comparison with other methods and, in \S(ref), we tackle the problem of testing in the presence of nuisance parameters. In \S(ref), we present an empirical application. The online Supplementary Material contains detailed derivations, technical appendices and additional numerical experiments.
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:
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 LY10 p.\ 172, we generate samples from ((ref)) 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), where we apply the model in ((ref)) 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), 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), we display the geometry of a random field on a regular lattice (undirected graph). In the real data example of \S (ref), we consider a random field over a manifold (a sphere), providing two additional examples for $W_n$.
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 LY10, with maximum likelihood estimation method; we refer to the R package 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) for details). Via QQ-plot analysis, Figure (ref) 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:
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$. \\
Let us consider a random field described by the {SARAR(1,1) model in ((ref)). {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 LY10 and 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 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)) and we obtain:
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 LY10 p. 167.} Thus, defining $\zeta =(\beta',\lambda,\rho)'$, the log-likelihood is:
where $V^*_{nt}(\zeta) = R_n(\rho) [S_n(\lambda) Y^*_{nt} - X^*_{nt} \beta].$ {As remarked in 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:
where $ \tilde{V}_{nt}(\zeta) = R_n(\rho) [S_n(\lambda) \tilde{Y}_{nt} - \tilde{X}_{nt} \beta], $ with
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:
{ where $\psi_{nt}$ is the likelihood score function}
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}.$
We assume $ n \gg T$, so we deal with so-called 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 LY10. In contrast, in our development, we consider small sample cross-sectional asymptotics (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); (ii) identifying the corresponding $U$-statistic, see \S(ref); (iii) deriving the Edgeworth expansion for the $U$-statistic as in BGVZ86 and deriving the saddlepoint density by means of the tilted-Edgeworth technique, see \S (ref) and \S(ref). Similar approaches are available in the standard setting of i.i.d.\ random variables in ER86, BNSC89, and GR96.
Let us first define the $M$-functional related to the MLE. To this end, we remark that the likelihood score function in ((ref)) 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)), it follows
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)), it follows that the MLE is the solution to
The $M$-functional $\vartheta$ related to the MLE is implicitly defined as the unique functional root of:
{or equivalently via the asymptotic maximization $\theta_0 = \text{arg}\max_{\theta \in \Theta} \mathbb{E}[\ell_{n,T}(\theta_0)]$; see e.g., 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)) is the $M$-estimator defined in ((ref)), 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 LY10 and reference therein) as a preliminary estimator; for a book-length description of one-step procedure; see, e.g., vdW98 Ch. 5.
In what follows, we set $m:=n(T-1)$, with $m\rightarrow\infty$, as $n\rightarrow\infty$. Then, we introduce
Assumption A.{\it
}
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. 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 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
is non-singular. {Under Assumption A$(i)$-A$(iv)$, Theorem 1 part(ii) in LY10 shows that $\displaystyle \lim_{n \rightarrow \infty} \hat\theta_{n,T} = \theta_0 $.} Furthermore, Theorem 2 point (ii) in 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.
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 \\ Assumption B. {\it
} Then, we state the following
In ((ref)), 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 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., vdW98 and BaltagiBook for a book-length introduction. The second term in ((ref)), of order $m^{-2} \propto n^{-2}$, plays a pivotal role in our derivation of higher-order approximation.
The result of Lemma (ref) 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), 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.
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
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.
Making use of Lemma (ref) and Lemma (ref), 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)), 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 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 BGVZ86 for i.i.d. random variables. To elaborate further, we introduce
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$,
}
A few comments are in order. Assumptions D$(i)$-$(iii)$ are similar to the technical assumptions in 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 BGVZ86: compare, e.g., our D$(ii)$ to their Eq. (1.17). D$(iv)$ is not considered in 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)), 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}
In addition, we can get the saddlepoint density approximation by exponentially tilting the Edgeworth expansion.
The proofs of these Propositions are available in Appendix A. They rely on a argument similar to the one applied in the proof of F82 for the derivation of a saddlepoint density approximation of multivariate $M$-estimators, and in GR96 for general statistics. Following D80, we can further normalize $p_{n,T}$ to obtain a proper density by dividing the right hand side of ((ref)) by its integral with respect to $z$. This normalization typically improves even further the accuracy of the approximation.
The expansions in Proposition (ref) and Proposition (ref) are connected with the results on higher-order expansions available in the spatial econometric literature, as cited in \S (ref). However, some key differences are worth a mention.
(i) The Edgeworth expansion in 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 RR14_EJ,RR14_ET, RR15, HM18, and 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) and Lemma (ref) 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 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., 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 RR15. Recently, 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), we prove that our saddlepoint density approximation features 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 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.
Most of the quantities related to the saddlepoint density approximation $p_{n,T}$ and the tail area in ((ref)) 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.
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); we refer to the online Supplementary Material (Appendix D) for more details and for additional results.
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). 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), 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.
Saddlepoint vs Edgeworth expansion (testing simple hypotheses). The Edgeworth expansion derived in Proposition (ref) 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) 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%.
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), 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.
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 Sun11 and to R routine 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.
Our saddlepoint density and/or tail approximations are helpful for testing simple hypotheses about $\theta_0$; see \S (ref). 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 HM18 and HM20. Combined with the saddlepoint density in ((ref)), 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 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 RRY03 defines the test statistic in the i.i.d.\ setting, while Lo09 and 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:
where $\theta= (\lambda,\theta_2)'$, with $\theta_{10}=\lambda_0$ and $\theta_2=(\beta, \rho, \sigma^2)'$. Then, we define the test statistic
The function $\mathcal{K}_\psi(\nu,(\lambda,\theta_2))$ is the c.g.f. of the estimating function:
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). We highlight that the expected value in ((ref)) 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.
To implement the test ((ref)) for the problem ((ref)), in Appendix C.2, we propose an algorithm (see Algorithm 2) and itemize the main steps needed to compute the test statistic.
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)), 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)) we make use of the R routine 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). 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.
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. 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 DE10. Since the cross-sectional size is only $n=24$, the asymptotics may suffer from size distortion as documented in \S(ref). 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:
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 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). Table (ref) gathers the point estimates (and their standard errors) that agree with the magnitudes found by DE10. To investigate the validity of the model ((ref)), 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) 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)), 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). 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.
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.