# Nonlinear impact estimation in spatial autoregressive models

Jean-Sauveur Ay, Julie Le Gallo, Kassoum Ayouba

► To cite this version:

Jean-Sauveur Ay, Julie Le Gallo, Kassoum Ayouba. Nonlinear impact estimation in spatial autoregressive models. Economics Letters, 2017, 163, pp.59-64. 10.1016/j.econlet.2017.11.031. hal-01662244

HAL Id: hal-01662244

https://hal.science/hal-01662244v1

Submitted on 16 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Nonlinear impact estimation in spatial autoregressive
models

Jean-Sauveur Ayª, Kassoum Ayoubaª,*, Julie Le Galloª

ªCESAER, AgroSup Dijon, INRA, Université de Bourgogne Franche-Comté, 21000 Dijon,
France

Abstract

This paper extends the literature on the calculation and interpretation of im-
pacts for spatial autoregressive models. Using a Bayesian framework, we show
how the individual direct and indirect impacts associated with an exogenous
variable introduced in a nonlinear way in such models can be computed, the-
oretically and empirically. Rather than averaging the individual impacts, we
suggest to graphically analyze them along with their confidence intervals cal-
culated from Markov chain Monte Carlo (MCMC). We also explicitly derive
the form of the gap between individual impacts in the spatial autoregressive
model and the corresponding model without a spatial lag and show, in our ap-
plication on the Boston dataset, that it is higher for spatially highly connected
observations.

Keywords: Spatial econometrics, marginal impacts, spline, Markov chain
Monte Carlo

JEL: C11, C21

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: jean-sauveur.ay@inra.fr (Jean-Sauveur Ay),
ayouba.kassoum@inra.fr (Kassoum Ayouba), julie.le-gallo@agrosupdijon.fr (Julie Le
Gallo)</footnote>
Nonlinear impact estimation in spatial autoregressive
models

Abstract

This paper extends the literature on the calculation and interpretation of im-
pacts for spatial autoregressive models. Using a Bayesian framework, we show
how the individual direct and indirect impacts associated with an exogenous
variable introduced in a nonlinear way in such models can be computed, the-
oretically and empirically. Rather than averaging the individual impacts, we
suggest to graphically analyze them along with their confidence intervals cal-
culated from Markov chain Monte Carlo (MCMC). We also explicitly derive
the form of the gap between individual impacts in the spatial autoregressive
model and the corresponding model without a spatial lag and show, in our ap-
plication on the Boston dataset, that it is higher for spatially highly connected
observations.

Keywords: Spatial econometrics, marginal impacts, spline, Markov chain
Monte Carlo

JEL: C11, C21
# 1. Introduction

Spatial autoregressive (SAR) models are now widely used to analyze spa-
tial economic interactions in various applications. For example, such a model
is appropriate for the housing market since housing prices depend on prices of
recently sold neighboring homes (Anselin and Lozano-Gracia, 2008). This de-
pendence structure comes from the fact that sellers presumably use information
on neighboring homes to determine the asking price. The use of a SAR model,
besides providing a richer characterization of the market, has important impli-
cations for impacts calculations (for early statements of this issue, see Anselin
and Le Gallo (2006), Kim et al. (2003), Kelejian et al. (2006)). Indeed, their
computation and interpretation is less straightforward than in standard multi-
ple linear a-spatial regression models: any change in an explanatory variable for
a given observation not only affects the observation itself (direct impact) but
also affects all other observations indirectly (indirect impact).

LeSage and Pace (2009) show how to theoretically derive these marginal im-
pacts in SAR models. For *n* spatial observations, they obtain a *n* × *n* matrix
of impacts for one exogenous variable. In order to have a compact represen-
tation of these impacts, they propose to report one direct impact equal to the
average of the diagonal elements of the matrix of marginal impacts and one
indirect impact equals to the average row sums of the non-diagonal elements of
that matrix. However, when the exogenous variable of interest is introduced in
a nonlinear way in the SAR model (e.g. in the form of polynomial or splines
function), averaging the impacts in such a way is irrelevant.

In this paper, we extend the work of LeSage and Pace (2009) on impacts
computation and estimation when the exogenous variable of interest appears in
a nonlinear way in a SAR model. We also derive the form of the gap between
impacts in the spatial autoregressive model and the corresponding model with-
out a spatial lag and show in our application that it is higher for spatially highly
connected observations.

The remainder of the paper is organized as follows. The second section
presents the theoretical derivation of direct and indirect impacts associated with
an exogenous variable introduced in a nonlinear way in a SAR model. The
third section introduces the estimation strategy which is applied on the well-
known Boston housing dataset described in the fourth Section. The obtained
empirical results are presented in the fifth Section. The sixth section concludes
and comments on possible extensions.
## 2. Impacts in theory

Consider the following spatial autoregressive model<fnref n="1" />

$$y = \rho Wy + X\beta + f(z) + \epsilon \quad (1)$$

where $y$ is the $(n \times 1)$ dependent variable that exhibits variation across spatial observations, $X$ is the $(n \times k)$ matrix of linear explanatory variables including an intercept term, with the associated parameters $\beta$ contained in a $(k \times 1)$ vector, $z$ is the $(n \times 1)$ variable the impact of which on $y$ is nonlinear, and $W$ is a specified constant $(n \times n)$ spatial weight matrix with the usual assumptions. $z$, our variable of interest is additively separable from the other $X$ to simplify notations. Our approach can easily be extended to many variables introduced non-linearly with interaction effects between them. We assume that each term of the disturbance vector $\epsilon$ of dimension $(n \times 1)$ is normally and identically distributed with zero mean and variance $\sigma^2$. The scalar $\rho$ measures the strength of the spatial dependence. $f(\cdot)$ is a linear-in-parameters function, for instance a polynomial function of degree $p$: $f(z) = \sum_{j=1}^{p} \gamma_j z^j$, a spline<fnref n="2" /> function of order $p$ and $q$ knots $t_l, t_{l+1}, \dots, t_q$: $f(z) = \sum_{j=1}^{p} \gamma_j z^j + \sum_{l=1}^{q} \delta_l (z - t_l)_+^p$, a B-spline function of order $p$ and $q$ knots: $f(z) = \sum_{j=1}^{p+q} \gamma_j B_j(z)^3$

Equation (1) can be rewritten in the following reduced form

$$y = V(W)X\beta + V(W)f(z) + V(W)\epsilon \quad (2)$$

with

$$V(W) = (I_n - \rho W)^{-1} = \begin{bmatrix} V_{11} & V_{12} & V_{13} & \cdots & V_{1n} \\ V_{21} & V_{22} & V_{23} & \cdots & V_{2n} \\ V_{31} & V_{32} & V_{33} & \cdots & V_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ V_{n1} & V_{n2} & V_{n3} & \cdots & V_{nn} \end{bmatrix}$$

and $I_n$ the identity matrix of order $n$. Note that we assumed that the matrix $I_n - \rho W$ is not singular to reach the reduced form in Equation (2).

In this paper, we are interested in the estimation of the partial derivative of $y$ with respect to changes in the nonlinear variable $z$ in our SAR model. In models containing a spatial lag, the measure of the partial derivative of the dependent variable with respect to an explanatory variable is less straightforward than in standard linear models. Indeed the standard linear regression interpretation of coefficient estimates ($\hat{\beta}_q = \frac{\partial \hat{y}}{\partial x_q}$) as partial derivatives no longer holds in SAR

<footnote n="1">¹Our approach can also be applied on spatial Durbin model.</footnote>

<footnote n="2">²See Hastie and Tibshirani (1990) for details.</footnote>

<footnote n="3">³Note that the B-spline function (B stands for “basis”) is an extension of the spline function. The $B_j, j = 1, \dots, p+q$ are defined recursively to avoid numerical instability.</footnote>
model since the matrix of explanatory variable is transformed by the $n \times n$ in-
verse matrix $V(W)$. In such a model, any change to an explanatory variable
for a given observation affects the dependent variable of the observation itself
(direct impact) and potentially the dependent variable of all other observations
(indirect impact) through $V(W)$. We elaborate on this observation to derive
the impacts for the variable $z$ introduced in a nonlinear way.

Starting from the reduced form in Equation (2), the matrix of responses to
a change of the nonlinear variable $z$ on $y$ is given by

$$ \begin{bmatrix} \frac{\partial y_1}{\partial z_1} & \frac{\partial y_1}{\partial z_2} & \cdots & \frac{\partial y_1}{\partial z_n} \\ \frac{\partial y_2}{\partial z_1} & \frac{\partial y_2}{\partial z_2} & \cdots & \frac{\partial y_2}{\partial z_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial y_n}{\partial z_1} & \frac{\partial y_n}{\partial z_2} & \cdots & \frac{\partial y_n}{\partial z_n} \end{bmatrix} = \begin{bmatrix} V_{11}f_z(z_1) & V_{12}f_z(z_2) & \cdots & V_{1n}f_z(z_n) \\ V_{21}f_z(z_1) & V_{22}f_z(z_2) & \cdots & V_{2n}f_z(z_n) \\ \vdots & \vdots & \ddots & \vdots \\ V_{n1}f_z(z_1) & V_{n2}f_z(z_2) & \cdots & V_{nn}f_z(z_n) \end{bmatrix} $$

with $f_z(z)$ the derivative of $f(z)$.

Since we have an $n \times n$ matrix of impacts, the challenge is to find a way
to compactly report them. In the case of an exogenous variable that enters
linearly, LeSage and Pace (2009) suggest to compute the average of the main
diagonal elements and the average of the off-diagonal elements of the impacts
matrix to obtain a summary measure of the direct impact and the indirect
impact respectively. This method of summarizing impacts is irrelevant when
the effect of a variable as $z$ can be positive and negative at different parts of it
support. Therefore, we propose to plot the individual impacts (i.e. impacts of
each observation) along with their confidence intervals. The individual direct
(IDI) and total (ITI) impacts of the nonlinear variable $z$ on $y_i$ are given by<fnref n="4" />

$$ IDI_i = \frac{\partial y_i}{\partial z_i} = V_{ii}f_z(z_i) \qquad (3) $$

$$ ITI_i = \frac{\partial y_i}{\partial z} = \sum_{j=1}^{n} V_{ij}f_z(z_j) \qquad (4) $$

From these expressions, we can reach the individual indirect impact (III) as

<footnote n="4">⁴Note that our approach is also appropriate to the case of an interacted explanatory variable, e.g., $y = \rho Wy + z\beta + xz\gamma + \epsilon$, where the individual direct and total impact would be respectively $\frac{\partial y_i}{\partial z_i} = V_{ii}(\beta + x_i\gamma)$ and $\frac{\partial y_i}{\partial z} = \sum_{j=1}^{n} V_{ij}(\beta + x_j\gamma)$. This logically means that the values of the variable $x$ at each observation produce observation-level impact estimates of interest. To assess the role played by $x$ in determining observation-level effects, we could plot these impacts estimates for sorted values of $x$. We thank the reviewer for pointing out this.</footnote>
follows

$$III_i = ITI_i - IDI_i \quad (5)$$

### 3. Estimation strategy

To get the individual impacts for z, we proceed in two steps. We first estimate the parameters of the SAR model, and then in a second step use them to compute $V(W)$ and the derivative of $f(z)$, i.e. $f_z(z)$ for each observation i. Note that, to estimate $V(W) = (I_n - \rho W)^{-1}$, we have followed an approach suggested by LeSage and Pace (2009, page 111), based on the Leontief expansion of $V(W)$, which is computationally more efficient than the brute force approach (i.e. calculating $V(W)$ per se) for large n: $V(W) = (I_n - \rho W)^{-1} \approx I_n + \rho W + \rho^2 W^2 + ... + \rho^q W^q$. We have chosen q = 250 to approximate $V(W)$ in our application.<fnref n="5" />

Let us begin with the estimation of the SAR model. We follow LeSage and Pace (2009) who have shown how to estimate SAR models using a Bayesian methodology. We briefly recall the approach in the next paragraph.<fnref n="6" />

The Bayesian model or alternatively the posterior distribution of the model parameters up to a constant term is given by the product of the likelihood of the available data (y, X and z) and the prior distributions of the parameters of the model. To get the mean and variance for the model parameters, an option is to numerically integrate the posterior distribution. Unfortunately, the form of the posterior distribution is often complex, leading to intractable numerical integration. To circumvent this problem, we rely on the methodology known as Markov Chain Monte Carlo (MCMC) to estimate the SAR model parameters. The general idea behind MCMC is that rather than work with the posterior density of our parameters, the same goal could be achieved by analyzing a large random sample from the posterior distribution. Specifically, we use the *Metropolis within Gibbs sampling* since our posterior distribution can be split in three posterior distributions: one of unknown form for $\rho$ and two for $\sigma^2$ and the parameters associated to the linear and nonlinear variables following a Gamma inverse and a Normal distributions respectively.

Once we have estimated the parameter of our model, to get the individual impacts, we need the derivative of $f(z)$, i.e. $f_z(z)$. A straightforward way of getting it is to analytically compute the derivative of $f(z)$. For instance, when $f(z)$ is a polynomial function of degree p, its derivative is given by: $f_{z_i}(z) =$

<footnote n="5">⁵The trace-based approach proposed by LeSage and Pace (2009, pp. 114-115) to approximate the direct and indirect impacts (defined respectively as the average of the diagonal elements of the matrix of marginal impacts and the average row sums of the non-diagonal elements of that matrix) does not apply to our case since we are not interested in aggregated measure of impacts.</footnote>

<footnote n="6">⁶A full presentation of the procedure can be found in LeSage and Pace (2009, pp. 124-153).</footnote>
$\sum_{j=1}^{p} j\gamma_j z^{j-1}$. However, in some cases (e.g. B-splines functions), using the analytical derivative of $f(z)$ may not be that simple and it is better to proceed differently. In this paper, we make use of simulations to get the derivatives. Starting from the SAR model, the derivative of $f(z)$ for observation $i$ is given by

$$f_z(z_i) = \lim_{\xi \to 0} \frac{\hat{y}_{i2} - \hat{y}_{i1}}{\xi} \qquad (6)$$

where $\hat{y}_{i1}$ and $\hat{y}_{i2}$ are predictions of $y_i$ from the model in Equation (1) using $z$ and $z + \xi$ respectively. All the estimated parameters except those related to $f(z)$ are set their mean when calculating $\hat{y}_{i1}$ and $\hat{y}_{i2}$.

## 4. Data

We use the Boston dataset for the application. It is available to the R community as part of the library spdep (Bivand et al. (2013), Bivand and Piras (2015)). The dataset reports for each of the units of observation (506 in total), the (corrected) median value of owner-occupied homes in the Boston area (CMEDV), together with several variables which might help to explain its variation such as: per-capita crime rate (CRIM) by town, nitric oxides concentration (NOX), proportions of non-retail business acres per town (INDUS), proportions of units built prior to 1940 (AGE), weighted distances to five Boston employment centers (DIS), index of accessibility to highways (RAD), property-tax rate (TAX) per $10,000, pupil-teacher ratios by town (PTRATIO); proportion of blacks (B) and percentage of the lower status of the population (LSTAT).

## 5. Results and discussion

We apply the methodology described above to Boston house prices. Specifically, we estimate the SAR specification in Equation (1), where $y$ is the median value of owner-occupied homes (CMEDV), $X$ the vector of the following linear variables: TAX, INDUS, CRIM and AGE. The nonlinear variable $z$ is represented by DIS. $f$ is a cubic B-spline function ($f(DIST) = \sum_{j=1}^{3} \gamma_j B_j(DIST))^7$ used to capture the potential nonlinear effect of DIS on the dependent variable. The spatial weights matrix $W$ is constructed using the sphere of influence approach (see Bivand (2010, page 66) for details).

### 5.1. Results of the Bayesian estimation

We use the following starting values for the parameters of our Bayesian model: 0.5 for $\sigma^2$, 0.1 for $\rho$ and 0.5 for each of the remaining parameters (MCMC

<footnote n="7">⁷Since we do not have an a priori on the number and the position of knots on the support of $z$, we set $q = 0$.</footnote>
starting values). We generate a total of 10,000 simulations, and then discard the first 1,000 as a “burn-in” – to mitigate startup effects. The Metropolis-Hastings algorithm used to draw values for $\rho$ perform very well: the acceptance probability is 48%.<fnref n="8" /> The chains of retained simulations show no sign of nonstationarity (see supplemental Figure A online). Note that we provide some more formal diagnostics to the visual (graphical) convergence check. Specifically, we monitor the convergence of our chains of parameters using diagnostics suggested by Geweke (1992) and Gelman and Rubin (1992). They provide no evidence against convergence for our chains.

Table 1 presents point estimates (i.e. posterior means), associated measures of uncertainty (i.e. posterior standard deviations) and the upper and lower bound of a 95% highest posterior density interval<fnref n="9" /> (columns labeled Q2.5% and Q97.5%) for the parameters of our model.

<table><caption>Table 1: Parameter estimates</caption><thead><tr><th></th><th>MEAN</th><th>STD</th><th>Q2.5</th><th>Q97.5</th></tr></thead><tbody><tr><td>INTERCEPT</td><td>22.8258</td><td>2.3160</td><td>18.2907</td><td>27.4274</td></tr><tr><td>B<sub>1</sub>(DIST)</td><td>-13.4145</td><td>3.4085</td><td>-20.1183</td><td>-6.7841</td></tr><tr><td>B<sub>2</sub>(DIST)</td><td>-9.3020</td><td>3.2576</td><td>-15.7103</td><td>-2.8356</td></tr><tr><td>B<sub>3</sub>(DIST)</td><td>-11.8411</td><td>4.1098</td><td>-19.7626</td><td>-3.9209</td></tr><tr><td>CRIM</td><td>-0.1101</td><td>0.0342</td><td>-0.1781</td><td>-0.0432</td></tr><tr><td>TAX</td><td>-0.0051</td><td>0.0022</td><td>-0.0094</td><td>-0.0007</td></tr><tr><td>AGE</td><td>-0.0595</td><td>0.0131</td><td>-0.0847</td><td>-0.0336</td></tr><tr><td>INDUS</td><td>-0.3057</td><td>0.0603</td><td>-0.4241</td><td>-0.1863</td></tr><tr><td>σ<sup>2</sup></td><td>25.9011</td><td>1.6983</td><td>22.7676</td><td>29.4221</td></tr><tr><td>ρ</td><td>0.6890</td><td>0.0290</td><td>0.6285</td><td>0.7448</td></tr></tbody></table>

We wish to compare the Bayesian MCMC sampling procedure and the maximum likelihood method results.<fnref n="10" /> To this end, the maximum likelihood estimates for the SAR model and sample data are presented in Table 2 alongside Bayesian estimates based on MCMC sampling scheme. We also report the t-statistics for the two estimations methods as in LeSage and Pace (2009, pp. 143-145). We observe that all estimates and t-statistics are nearly identical, suggesting they would produce similar inferences.

<footnote n="8">⁸LeSage and Pace (2009) state that the optimal acceptance rate (i.e., the one which minimizes the autocorrelations across the sample values) is between 40% and 60%.</footnote>

<footnote n="9">⁹This coverage region is analogous to the 95% confidence intervals.</footnote>

<footnote n="10">¹⁰The maximum likelihood estimation results are available in supplemental Table A online.</footnote>
<table><caption>Table 2: Comparison of SAR model estimates</caption><thead><tr><th rowspan="2"></th><th colspan="2">Max. Likelihood</th><th colspan="2">Bayesian M-H</th></tr><tr><th>Estimates</th><th>t-stat</th><th>Estimates</th><th>t-stat</th></tr></thead><tbody><tr><td>INTERCEPT</td><td>22.5830</td><td>9.1977</td><td>22.8258</td><td>9.8556</td></tr><tr><td>B<sub>1</sub>(DIST)</td><td>-13.3724</td><td>-3.9837</td><td>-13.4145</td><td>-3.9355</td></tr><tr><td>B<sub>2</sub>(DIST)</td><td>-9.2372</td><td>-2.8515</td><td>-9.3020</td><td>-2.8554</td></tr><tr><td>B<sub>3</sub>(DIST)</td><td>-11.7986</td><td>-2.9195</td><td>-11.8411</td><td>-2.8811</td></tr><tr><td>CRIM</td><td>-0.1086</td><td>-3.2452</td><td>-0.1101</td><td>-3.2144</td></tr><tr><td>TAX</td><td>-0.0050</td><td>-2.2547</td><td>-0.0051</td><td>-2.2828</td></tr><tr><td>AGE</td><td>-0.0590</td><td>-4.4745</td><td>-0.0595</td><td>-4.5202</td></tr><tr><td>INDUS</td><td>-0.3029</td><td>-5.0132</td><td>-0.3057</td><td>-5.0650</td></tr><tr><td>σ<sup>2</sup></td><td>25.2490</td><td></td><td>25.9011</td><td></td></tr><tr><td>ρ</td><td>0.6949</td><td>22.805</td><td>0.6890</td><td>23.6795</td></tr></tbody></table>

## 5.2. Impacts estimation

For each of the 10,000 iterations, we compute the individual direct, total and indirect impacts of the variable *DIS* introduced in a nonlinear way in the specification. Let us focus on the individual direct impacts to clearly get how they are calculated. For observation *i*, we compute for each iteration the individual direct impact for the corresponding value of *DIS* using the parameter estimates and the derivative of the nonlinear function *f*, with equation (6). We end up with a chain of individual direct impacts corresponding to the number of iterations. As above, we discard the first 1,000 draws. On the remaining chain, we calculate the mean and the 2.5% and 97.5% percentiles which correspond respectively to the estimate of the individual direct impacts and its 95% confidence interval. The procedure is repeated for each of the 506 observations of *DIS*. The estimation of the individual total and indirect impacts follows the same procedure.

Figure 1 plots for each value of *DIS*, the direct and indirect impacts along with their confidence intervals. In the top panels of Figure 1, we plot the direct and indirect impacts of *DIS* for each observation *i*. We have also applied a LOESS regression to these two sets of data point to get some more readable curves of impacts (bottom panels)
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of direct and indirect individual impacts</figcaption></figure>

Two important observations follow from Figure 1. First, we observe that the impacts (direct and indirect) are negative, suggesting that Boston employment centers are attractive: the closer a house is to city center, the higher its price and the price of its neighbors. Second, the marginal impacts are nonlinear: the more *DIS* increases, the less are the magnitude of the impacts. This means that
an increase of *DIS* of one mile has a higher impact on values of homes close to the center compared to the others. The non-linearity of the impacts implies that taking their average is indeed irrelevant in our case, and may lead to erroneous conclusions.

To better see the nonlinear effect of *DIS*, we directly plot it against the dependent variable *CMED*.<fnref n="11" /> For that purpose, for the value of *DIS* for observation *i*, we calculate for each iteration its prediction using the parameter estimates. Note that we use the unbiased predictor proposed by Kelejian and Prucha (2007, page 367) which is compatible with our SAR model in Equation (2): the prediction is based on the information set $\Lambda = \{Z, W\}$ where Z is the vector of linear and nonlinear exogenous variables. To get rid of the effect of the other exogenous variable and focus only on *DIS*, we set as in Fox (1987) their values to their mean. From there, we follow the same procedure as in the calculation of marginal impacts to end with an estimate of the effect of distance and the 95% confidence interval. The procedure is repeated for each of the 506 values of *DIS*. The obtained plot (see Figure 2) clearly shows the nonlinear effect of *DIS* on the endogenous variable *CMED*.

<footnote n="11">¹¹The coefficients associated to splines variables in Table 1 are not directly interpretable.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Impact of DIS</figcaption></figure>

### 5.3. Structure of W and impacts in nonlinear models

In this section, we analyze the link between the spatial structure of the observations and the form of the impacts in nonlinear models in the context of the application. To that purpose, let us begin with the individual direct impacts. From Equation (3), we know that when there is no spatial autoregressive term (i.e. when the model is a-spatial), the individual direct impact is $f_z(z_i)$. Therefore, the difference between spatial and a-spatial models individual direct impacts is the following

$$ GAP_{i}^{Direct} = f_{z}(z_{i}) [V_{ii} - 1] \qquad (7) $$

Two observations emerge from the expression in Equation (7). First, the gap between the two individual direct impacts depends on the structure of spatial connections between observations ($w_{ij}$, $i, j = 1, ..., n$). Second, from the Leontief expansion ($V(W) = [I_n - \rho W]^{-1} = I_n + \rho W + \rho^2 W^2 + ...$), the elements of $V(W)$ are always positive, and even greater than one for the diagonal terms ($V_{ii}$), since the estimate of the spatial autoregressive coefficient ($\rho$) is positive. Therefore the sign of the gap between the individual impacts in the spatial and the a-spatial models is of the sign of $f_z(z_i)$. As a consequence, $f_z(z_i)$ is the upper
bound (when $f_z(z_i)$ is negative) or the lower bound (when $f_z(z_i)$ is positive) estimation of the individual impacts in our spatial autoregressive model.

Let us move now on individual total impacts. From Equation (4), we can deduce that for an a-spatial model the individual total impact is $f_z(z_i)$, which is the same as the individual direct impact.<fnref n="12" /> Therefore, the difference between spatial and a-spatial individual total impacts is the following

$$GAP_{i}^{Total} = f_{z}(z_{i}) [V_{ii} - 1] + \sum_{j=1, j \neq i}^{n} V_{ij}f_{z}(z_{j}) \quad (8)$$

We note that the gap between the two individual total impacts, as well, depends on the structure of spatial connections between observations. Unlike for individual direct impacts, we can not say anything about the sign of $GAP_{i}^{Total}$.

To confirm these observations, we plot for each observation i the direct and total impacts against $f_z(z_i)$, which is the total impact for the a-spatial model. The procedure is repeated for each observation in the Boston dataset and the obtained results are plotted in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Impact in the a-spatial model and in the SAR model</figcaption></figure>

The analysis of Figure 3 confirms our previous observations. First, the

<footnote n="12">¹²This makes sense as there is no individual indirect impact in a-spatial models.</footnote>
derivative $f_z(z_i)$, indeed appears as the upper bound estimation of the direct impact. It is even the upper bound estimation of the total impact, since all the $f_z(z_i)$ are negative in our application.<fnref n="13" /> This means in our case that the estimated total impacts in the a-spatial model are overestimated if the data generating process is a SAR model. Second, the gap between the two curves is increasing as we move towards the center (low distances). Since the observations close to the center are often the most connected (dense zones), that suggests that the more connected is the observation, the more is the gap (i.e. the gap depends on spatial connections among observations).

To confirm this intuition, we focus on the structure of the spatial weights matrix. Specifically, we calculate for each observation, its number of neighbors defined as the number of houses in a 1,000 meters buffer. Then we plot for each observation the gap between $f_z(z_i)$ and the total impact against its number of neighbors. The result is presented in Figure 5 and clearly shows that central observations are associated with higher gap. This result is interesting as it states that the more an observation is linked to the others, the more important is the bias if spatial autocorrelation is wrongly omitted. In the Boston housing market, our result states that median housing prices decrease more rapidly for houses located in dense areas than for the houses located in sparse areas.

This result has another important implication: it confirms that averaging the impacts may overlook important heterogeneity even if the exogenous variable of interest is linear. Indeed, for a linear exogenous variable $x_k$ with the associated parameter $\beta_k$, the individual total impact is $\beta_k \sum_{j=1}^n V_{ij}$. If $\beta_k$ is for example positive, considering only the average leads to the underestimation of the impact of $x_k$ for observations in central areas and the overestimation in sparse areas. Therefore, in specific applications where considering the heterogeneity of impacts is important, for instance to design place-based policies, it is better to take advantage of the “noise” introduced by spatial autocorrelation between observations rather than smoothing the impacts by averaging them, even for linear exogenous variables.

<footnote n="13">¹³When $f_z(z_i) \le 0$ for $i = 1, ..., n$, $GAP_i^{Total}$ is $\le 0$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Relation between the number of neighbours and the GAP</figcaption></figure>

6. Conclusion

We provide in this paper a framework for the estimation of impacts associ-
ated with an exogenous variable introduced in a non linear way in spatial autore-
gressive model. We also show that instead of averaging the individual impacts,
one should consider to graphically analyze them along with their confidence
interval. Indeed, averaging the individual impacts, besides being inappropriate
for nonlinear exogenous variable, smooths spatial interaction effects which may
be of interest in spatial autoregressive models.
References

Anselin, L., Le Gallo, J., 2006. Interpolation of air quality measures in hedonic house price models: spatial aspects. Spatial Economic Analysis 1, 31–52.
doi:http://dx.doi.org/10.1080/17421770600661337.

Anselin, L., Lozano-Gracia, N., 2008. Errors in variables and spatial effects in hedonic house price models of ambient air quality. Empirical economics 34, 5–34.
doi:https://doi.org/10.1007/s00181-007-0152-3.

Bivand, R., Hauke, J., Kossowski, T., 2013. Computing the jacobian in gaussian spatial autoregressive models: An illustrated comparison of available methods. Geographical Analysis 45, 150–179.
URL: http://www.jstatsoft.org/v63/i18/, doi:http://dx.doi.org/10.1111/gean.12008.

Bivand, R., Piras, G., 2015. Comparing implementations of estimation methods for spatial econometrics. Journal of Statistical Software 63, 1–36.
URL: http://www.jstatsoft.org/v63/i18/.

Bivand, R.S., 2010. Spatial econometric functions in r, in: Handbook of Applied Spatial Analysis. Fischer M., Getis A. (eds). Springer, Berlin, Heidelberg.
doi:https://doi.org/10.1007/978-3-642-03647-7_4.

Fox, J., 1987. Effect displays for generalized linear models. Sociological methodology 17, 347–361.
doi:https://doi.org/10.2307/271037.

Gelman, A., Rubin, D., 1992. Inference From iterative simulation using multiple sequences, in: Bayesian Statistics. J. Bernardo, J. Berger, A. Dawid and A. Smith, eds, Clarendon Press, Oxford, pp. 457-511.

Geweke, J., 1992. Evaluating the accuracy of sampling-based approaches to calculating posterior moments, in: Bayesian Statistics. J. Bernardo, J. Berger, A. Dawid and A. Smith, eds, Clarendon Press, Oxford, pp. 169-194.

Hastie, T., Tibshirani, R., 1990. Generalized additive models. Wiley Online Library.
doi:http://dx.doi.org/10.1002/0471667196.ess0297.pub2.

Kelejian, H.H., Prucha, I.R., 2007. The relative efficiencies of various predictors in spatial econometric models containing spatial lags. Regional Science and Urban Economics 37, 363–374.
doi:https://doi.org/10.1016/j.regsciurbeco.2006.11.005.

Kelejian, H.H., Tavlas, G.S., Hondroyiannis, G., 2006. A spatial modelling approach to contagion among emerging economies. Open economies review 17, 423–441.
doi:https://doi.org/10.1007/s11079-006-0357-7.

Kim, C.W., Phipps, T.T., Anselin, L., 2003. Measuring the benefits of air quality improvement: a spatial hedonic approach. Journal of environmental economics and management 45, 24–39.
doi:https://doi.org/10.1016/S0095-0696(02)00013-X.

LeSage, J., Pace, R.K., 2009. Introduction to Spatial Econometrics. CRC Press.
Online Supplementary Materials

Graphical convergence check

<figure><img src="image_6.png" /><figcaption>Figure A: MCMC chains</figcaption></figure>
Maximum likelihood estimation results

<table><caption>Table A: SAR estimation</caption><thead><tr><th></th><th><em>Dependent variable:</em></th></tr></thead><tbody><tr><td><em>B</em><sub>1</sub>(DIST)</td><td>-13.372<sup>***</sup><br>(3.357)</td></tr><tr><td><em>B</em><sub>2</sub>(DIST)</td><td>-9.237<sup>***</sup><br>(3.239)</td></tr><tr><td><em>B</em><sub>3</sub>(DIST)</td><td>-11.799<sup>***</sup><br>(4.041)</td></tr><tr><td>CRIM</td><td>-0.109<sup>***</sup><br>(0.033)</td></tr><tr><td>TAX</td><td>-0.005<sup>**</sup><br>(0.002)</td></tr><tr><td>AGE</td><td>-0.059<sup>***</sup><br>(0.013)</td></tr><tr><td>INDUS</td><td>-0.303<sup>***</sup><br>(0.060)</td></tr><tr><td>INTERCEPT</td><td>22.583<sup>***</sup><br>(2.455)</td></tr><tr><td>Observations</td><td>506</td></tr><tr><td>Log Likelihood</td><td>-1,574.987</td></tr><tr><td><em>ρ</em></td><td>0.69493</td></tr><tr><td><em>σ</em><sup>2</sup></td><td>25.249</td></tr><tr><td>Akaike Inf. Crit.</td><td>3,169.974</td></tr><tr><td>Wald Test</td><td>520.068<sup>***</sup> (df = 1)</td></tr><tr><td>LR Test</td><td>303.378<sup>***</sup> (df = 1)</td></tr></tbody></table>

Note: *p<0.1; **p<0.05; ***p<0.01