Revisiting the Spatial Autoregressive Exponential Model for Counts and Other Nonnegative Variables, with Application to the Knowledge Production Function

Abstract

This paper proposes a two-step pseudo-maximum likelihood estimator of a spatial autoregressive exponential model for counts and other nonnegative variables; it is particularly useful for dealing with zeros. It considers a model specification allowing us to easily determine the direct and indirect partial effects of explanatory variables (spatial spillovers and externalities). A simulation study shows that this method generally behaves better in terms of bias and root mean square error than existing procedures. An empirical example estimating a knowledge production function for the NUTS II European regions is analyzed. Results show that there is spatial dependence between regions on the creation of innovation, where regions less able to transform R&D expenses into innovation benefit from knowledge spatial spillovers through indirect effects. It is also concluded that the socioeconomic environment is important and that, unlike public R&D institutions, private companies are efficient at knowledge production.

1. Introduction

Many empirical applications with spatial data concern the modeling of counts and other nonnegative response variables. Examples are the modeling of trade flows, migration flows, patenting citation and patent creation, number of crashes, firm location and firm birth, number of new patients contracting a given disease, etc. Conventional practices opt to logarithmically transform the dependent variable in order to apply the well-known spatial linear models. This is the approach followed in [1] to model the interregional trade of goods at the NUTS3 level in Spain, in [2] to explain labor migration flows in China, and in [3] to investigate the effect of intraregional labor mobility in the production of knowledge in Europe, to give just a few examples. However, Silva and Tenreyro [4], in the context of cross-sectional data, note that modeling logarithmically transformed variables with a linear model may lead to bias in estimation when heteroscedasticity is present, or to distortions in parameter estimates caused by the need to add a constant to zero observations. The authors propose using the Poisson pseudo-maximum likelihood (PPML) estimator of the model for the untransformed variables as an alternative to ordinary least squares (OLS) of the loglinear model.

Spatial autoregressive models are popular to address spatial dependence. Elhorst [5] discusses the relevance of such models in recent applied spatial econometrics. One reason is that they quantify indirect spatial spillovers, as is pointed in [6,7]. While linear spatial autoregressive models are widely used in the literature, nonlinear spatial autoregressive models, namely models for counts or other nonnegative variables, are not so popular because of their complexity in estimation and in the derivation of marginal effects and spatial externalities (see, for example, Ref. [8], containing a review of the state of the art in spatial econometrics).

A spatial autoregressive Poisson model for counts was first introduced by Besag [9]. His specification includes the spatial lagged count variable in the exponential function that gives the Poisson conditional mean. However, this approach has a severe limitation in that it accommodates only negative spatial dependence to prevent the count process from being explosive. Attempts to overcome this limitation were made in [10], introducing a Winsorized Poisson, and [11], using a spatial filter that uses judiciously selected eigenvectors as regressors. However, these procedures may be computationally demanding and/or lack interpretability, namely because marginal effects and spatial externalities are hard or impossible to calculate given that the model is not invertible and, consequently, there is no reduced form for the dependent variable as a function of the so-called Leontief inverse or spatial multiplier ${(I-\rho W)}^{-1}$ (with $\rho$ the spatial autocorrelation coefficient and $W$ a spatial weighting matrix).

A promising approach is introduced by [12] to model counts that are Poisson distributed, considering that their conditional mean is a function of the logarithm of the spatially weighted conditional means of the neighbors (instead of the spatial lagged count variable as in [9]), in order to have invertibility, and call this model the SAR Poisson. This specification is easy to estimate, allowing for a simple calculation of marginal effects and spatial spillovers. However, because it involves the calculation of the logarithm of the count outcomes in estimation, it has to deal with observations that are zero.

Other approaches are based on Bayesian hierarchical modeling and Markov random fields, considering conditional autoregressive schemes for spatial errors like in [13,14]. These models are typically estimated by Bayesian Markov Chain Monte Carlo (MCMC), which is computationally very burdensome and does not allow for the calculation of spatial externalities. See [15] for an application of a Poisson Bayesian hierarchical model to patent citation flows, estimated with MCMC methods. A comprehensive review of spatial econometric models for count data is given in [16].

This paper proposes a PPML estimation procedure for spatial autoregressive exponential regression (SAR-E) that circumvents the problem of dealing with observations that are zero. The model specification is based on the idea of [12]. In a simulation study, the estimation procedure introduced in this work shows better properties in finite samples than the existing procedures, especially when the spatial autocorrelation coefficient is not close to one. The determination of partial effects is emphasized and the indirect effects are deduced, which enhances the understanding of regional linkages.

This work suggests applying this approach to model counts (not necessarily Poisson distributed) or, in general, other nonnegative variables, extending the applicability of [4] to spatial data.

To illustrate the usefulness of the approach introduced, an empirical study that estimates a knowledge production function (KPF) and knowledge spillovers across European regions is carried out, modeling the number of new patents per million inhabitants. Despite the existence of many applications that estimate KPFs, to our knowledge, only [17] uses a nonlinear SAR approach, namely a Bayesian spatial Tobit regression estimated by MCMC. The one presented here is the first that estimates a spatial autoregressive exponential model for counts, which circumvents the problem of regions with zero patents and simultaneously determines knowledge externalities across regions.

2. The SAR-E Regression

2.1. Model Specification and Partial Effects

This work proposes to model an outcome that is a count or other nonnegative variable, showing spatial dependence, by the following spatial autoregressive exponential specification of the conditional mean (SAR-E regression), which is based in the spatial lag model of counts of [12]:

$$E\left(y|X\right)=\mathsf{\mu }=\exp \left(\rho W\log (\mathsf{\mu })+X\mathsf{\beta }\right),$$

(1)

where $y$ is the vector with observations of the dependent variable for n spatial locations, $\mathsf{\mu }$ is the vector with a conditional mean of $y$, $X$ is a matrix with observations of k explanatory variables for n spatial regions, W is a spatial weighting matrix, and $\mathsf{\beta }$ and $\rho$ are unknown coefficients to be estimated. Observe that, according to Equation (1), the conditional mean of one location $i$, ${\mathsf{\mu }}_i$, is determined as a function of the characteristics of location $i$ through the observed values for the explanatory variables and of a weighted average of the conditional mean of neighboring locations.

Equation (1) serves three purposes. Firstly, it expresses the conditional expectation of a nonnegative and, in particular, a count variable. Observe that count variables are often assumed to have a Poisson distribution, whose expected conditional mean is an exponential function of a set of explanatory variables. Secondly, it incorporates the spatial dependence of the data by means of an autoregressive term, extending the well-known SAR or Spatial Lag linear model to the nonlinear context. Finally, it is invertible, which allows us to easily calculate the partial effects of variables and, in particular, analyze global spatial interactions between regions with the identification of spatial spillovers and externalities.

The reduced form of Equation (1) is:

$$E\left(y|X\right)=\mathsf{\mu }=\exp \left[{\left(I-\rho W\right)}^{-1}X\mathsf{\beta }\right].$$

(2)

The partial effects of the explanatory variables are deduced from Equation (2) leading to,

$$\frac{\partial \mathsf{\mu }}{\partial {x^{\prime }}_k}=\left[\begin{array}{ccc}\frac{\partial {\mu }_1}{\partial x_{1k}}& \cdots & \frac{\partial {\mu }_1}{\partial x_{nk}}\\ ⋮& \ddots & ⋮\\ \frac{\partial {\mu }_n}{\partial x_{1k}}& \cdots & \frac{\partial {\mu }_n}{\partial x_{nk}}\end{array}\right]={\beta }_k{\mathsf{\mu }}^{DIAG}{\left(I-\rho W\right)}^{-1},$$

(3)

where ${\mathsf{\mu }}^{DIAG}$ is a diagonal matrix of order n with elements ${\mu }_i$. Observe that Equation (3) is a $n\times n$ matrix of partial effects, where the elements in the main diagonal are direct effects of the kth explanatory while the off-diagonal elements are indirect effects.

Considering that $A={\left(I-\rho W\right)}^{-1}$, then the direct partial effects in region i are equal to

$$\frac{\partial {\mu }_i}{\partial x_{ik}}={\beta }_k{a}_{ii}{\mu }_i\hspace{1em}i=1,\dots ,n,$$

(4)

and give the expected outcome of a given location due to a one-unit variation of the kth explanatory in the same location.

When the spatial weighting matrix is row-normalized, indirect effects can be divided into spill-in and spill-out spatial spillovers. The spill-in spillover measures the cumulative sum of spatial spillovers that location i receives from all neighboring locations—that is, the sum of expected impacts in the outcome of location i due to a one-unit variation of the kth explanatory in each neighboring location j—and can be calculated as follows:

$$spill-i{n}_i={\beta }_k{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{a}_{ij}{\mu }_i}\hspace{1em}i=1,\dots ,n,$$

(5)

which is the cumulative sum of all off-diagonal elements in row i of Equation (3).

The spill-out spillover effect is the sum of all spatial spillovers that location i transfers to neighboring locations—that is, the sum of the expected impact in the outcome of each location j neighbor from i when the kth explanatory in location i varies by one unit and is equal to

$$spill-ou{t}_i={\beta }_k{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne 1\end{array}}^n{a}_{ji}{\mu }_j}\hspace{1em}i=1,\dots ,n,$$

(6)

or, equivalently, the sum of all off-diagonal elements in column i of Equation (3).

Each region has a direct, a spill-in and a spill-out partial effects. The values analyzed in empirical applications are usually the average of each of these effects over all regions, constituting, respectively, the average direct partial effect, the average spill-in spatial spillover (Aspill-in), and the average spill-out spatial spillover (Aspill-out), that is,

$$Aspill-in=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}spill-i{n}_i}\hspace{1em}Aspill-out=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}spill-ou{t}_i}$$

2.2. Estimation

When the dependent variable is a count with a Poisson distribution, the full information maximum likelihood (FIML) estimator of the reduced form in Equation (2) is derived in [12]. For variables that are not Poisson distributed, which includes some types of counts and other general nonnegative variables, those results can be used in a Poisson pseudo-maximum likelihood (PPML) context, assuming that the conditional mean is correctly specified according to Equation (1). Since the seminal work of [18], PPML has become popular for model estimation because it extends the technique of maximum likelihood to situations where the conditional distribution of the outcome does not need to be specified, but its conditional expectation has to be an exponential function of a linear index. The idea is to use the Poisson probability function to build the likelihood function, even if the outcome is not Poisson distributed, requiring that its expectation coincides with the expectation of a Poisson-distributed variable. As a misspecified distribution was used to define the likelihood function, the covariance matrix of the estimator and, in particular, standard errors, need to be estimated with a robust estimator. Silva and Tenreyro [4] disseminated PPML to estimate the gravity model, which is a particular case of exponential regression.

The PPML approach proposed here results in estimating the unknown coefficients by FIML and the respective standard errors with a robust estimator to safeguard them from variance misspecification, like in situations where there is overdispersion. However, the authors of [12] report severe difficulties in obtaining numerical solutions for FIML estimates. Therefore, they recommend instead a limited information maximum likelihood (LIML) two-step estimator. The first step delivers an estimate for the unknown variable $W\log (\mathsf{\mu })$ obtained with an OLS regression of $W\log (y)$ in the set of regressors $X,WX$ and $W^{2}X$. In the second step, a Poisson regression is performed with regressors $\widehat{W\log (y)}$ and $X$. An expression for the second stage adjusted covariance matrix is given in [12].

This paper proposes a two-step procedure that extends and refines the method described above in two ways. First, it suggests extending the estimation for a pseudo-maximum likelihood framework in order to encompass the modeling of a vast set of outcomes. This approach requires additional care in the estimation of the covariance matrix in the second step. Second, it proposes a different estimation procedure for the first step that circumvents the problem of observations that are zeros. Therefore, the following two-step PPML procedure to estimate the SAR-E model in Equation (1) is recommended:

• Run a PPML regression of $y$ on $X,WX$ and $W^{2}X$ and calculate the predicted values $\widehat{y}$.
• Run a PPML regression of $y$ on $W\log (\widehat{y})$ and $X$.

Observe that the second step of [12] uses $\widehat{W\log (y)}$, the fitted values for the variable $W\log (y)$, while the second step of the proposed method uses $W\log (\widehat{y})$, which is a transformation of the fitted values of variables $y$, $\widehat{y}$.

Standard errors in the second step should take into consideration the pseudo-maximum likelihood framework where the Poisson variance may be misspecified, and should account for the sampling variation in the regressor $W\log (\widehat{y}).$ To overcome these issues, the use of bootstrap standard errors is recommended in the second step. This procedure is easy to implement because it requires only software with a command for Poisson regression and bootstrap standard errors, like STATA [19]. (For non-negative outcomes other than counts, we advise using the command “glm” in STATA with the option “family(Poisson)” instead of “Poisson”.)

3. Simulation Study

In the simulation study, the two-step estimator introduced in Section 2 with a first step based on maximum likelihood (SAR-PPML 1stStep-ML) is compared with the two-step estimator of [12], with the first step being an OLS regression (SAR-PPML 1stStep-OLS) and aspatial PPML. For SAR-PPML 1stStep-OLS, when calculations need the logarithm of the outcome, an ad hoc constant equal to 1 was added to observations that were zeros. Simulations were performed with R [20].

3.1. Simulation Design

The simulation design closely follows that in [12], which is most closely related to other spatial econometric experimental designs, such as those in [21,22,23].

The random dependent variable is generated as $y_i\sim \mathrm{Poisson}({\mu }_i)$ with ${\mu }_i=\exp (A_{i}X\mathsf{\beta })$ where $A_i$ is the ith row of ${(I-\rho W)}^{-1}$. The design matrix includes two covariates, $X_1$ and $X_2$, where the first was randomly generated from a normal distribution, with mean and variance equal to 1 and 2, respectively. Since econometric studies usually incorporate a mix of continuous and dummy variables, following [4], $X_2$ is a dummy variable randomly generated from the Bernoulli distribution with a mean equal to 0.5.

The study considers three alternative spatial weight matrices. They were calculated using the same two-step procedure found in other spatial econometrics simulation studies (see, e.g., [24]). First, n space units are randomly drawn within in the unit square. Secondly, a matrix W0 is constructed given a criterion, and normalized by rows, so that the sum of all elements in each row is 1. In the present study, two different criteria were used, resulting in three alternative spatial weighting matrices. W1 and W3 intend to replicate matrices generated with a contiguity criterion, with neighbors chosen based on the nearest neighbor distance, fixing for W1 that each unit has seven neighbors (the seven units that are closest), while for W3 the number of neighbors is four, which is close to the average number of neighbors observed in the empirical study included in the next section. On the other hand, W2 is created based on an inverse distance criterion, using the Euclidean distance between each unit. The matrix W2 is said to be denser than the matrix W1, since W2 contains more nonzero entries, and matrix W1 is denser than matrix W3.

Monte Carlo simulations were conducted for each design of W and for each of the three estimators described above. The sample size, n, varied over the set [100; 250; 500; 750; 1000], and the spatial autoregressive parameter, $\rho$, varied over the set [0; 0.2; 0.4; 0.6; 0.8]. The parameters associated with variables $X_1$ and $X_2$, ${\beta }_1$ and ${\beta }_2$, respectively, were held fixed at 0.5. The intercept was set to zero.

For each experiment, 1000 replications were used. This is the usual number of replications used in Monte Carlo studies with spatial data (see [12,21,22,23,24,25], among others). The bias was calculated as the average in the 1000 replications of the difference between the estimated value of the coefficient in each simulated sample and the respective true value. The RMSE was also calculated for each estimated coefficient as the square root of the sum between the square of the bias and the empirical variance in the 1000 replications.

3.2. Monte Carlo Results

It should be noted that the results obtained referring to W1 are quite similar to those obtained with W3. This suggests that estimators should not be sensitive to the density of the spatial weighting matrix when using a contiguity criterion. For this reason, the analysis of the remaining results will focus only on experiments related to the use of W1 and W2 matrices. The results for W3 can be found in

Table A5. Bias: SAR-Poisson. SAR-LogLinear, and Aspatial ML Poisson with W3 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W3						β1-SAR-Poisson 1stStep-OLS—W3					β1-Aspatial Poisson ML—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0013	−0.0009	−0.0006	0.0000	0.0000	−0.0015	−0.0004	−0.0002	−0.0007	0.0002	−0.0008	−0.0004	−0.0004	0.0000	0.0001
0.2	−0.0002	−0.0003	0.0001	−0.0003	0.0002	−0.0068	−0.0057	−0.0049	−0.0051	−0.0049	0.0334	0.0320	0.0328	0.0327	0.0327
0.4	0.0018	0.0031	0.0030	0.0033	0.0031	−0.0074	−0.0067	−0.0061	−0.0066	−0.0059	0.0945	0.0955	0.0961	0.0950	0.0954
0.6	0.0053	0.0055	0.0049	0.0045	0.0076	−0.0043	−0.0036	−0.0040	−0.0047	−0.0041	0.2162	0.2243	0.2249	0.2249	0.2227
0.8	0.0333	0.0256	0.0271	0.0269	0.0258	0.0024	0.0008	0.0008	−0.0021	−0.0002	0.5565	0.5744	0.5710	0.5781	0.5701
β2-SAR-Poisson 1stStep-ML—W3						β2-SAR-Poisson 1stStep-OLS—W3					β2-Aspatial Poisson ML—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0016	−0.0020	0.0028	−0.0006	0.0001	−0.0069	−0.0032	−0.0038	−0.0011	−0.0013	−0.0014	−0.0011	0.0034	−0.0005	0.0000
0.2	0.0020	0.0014	0.0053	0.0043	0.0041	−0.0173	−0.0166	−0.0191	−0.0190	−0.0184	0.1017	0.1032	0.1052	0.1047	0.1046
0.4	0.0175	0.0206	0.0198	0.0221	0.0212	−0.0227	−0.0236	−0.0241	−0.0259	−0.0267	0.2822	0.2859	0.2864	0.2869	0.2874
0.6	0.0382	0.0308	0.0296	0.0263	0.0248	−0.0144	−0.0147	−0.0169	−0.0191	−0.0196	0.6469	0.6572	0.6612	0.6640	0.6653
0.8	0.0978	0.0919	0.0982	0.0490	0.0311	0.0067	0.0086	−0.0031	−0.0136	−0.0153	1.8195	1.9117	1.9474	1.9597	1.9623
Rho-SAR-Poisson 1stStep-ML—W3						Rho-SAR-Poisson 1stStep-OLS—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0006	0.0000	0.0002	−0.0013	−0.0002	0.0043	0.0006	0.0018	0.0015	−0.0005
0.2	−0.0008	0.0015	−0.0006	0.0007	0.0002	0.0193	0.0199	0.0215	0.0227	0.0219
0.4	−0.0050	−0.0086	−0.0081	−0.0088	−0.0090	0.0303	0.0329	0.0336	0.0358	0.0346
0.6	0.0002	−0.0041	−0.0048	−0.0059	−0.0114	0.0236	0.0247	0.0275	0.0296	0.0288
0.8	0.0118	0.0060	0.0003	−0.0034	−0.0015	0.0114	0.0112	0.0158	0.0239	0.0231

and

Table A6. RMSE: SAR-Poisson. SAR-LogLinear and Aspatial ML Poisson with W3 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W3						β1-SAR-Poisson 1stStep-OLS—W3					β1-Aspatial Poisson ML—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0315	0.0198	0.0134	0.0109	0.0090	0.0358	0.0212	0.0148	0.0119	0.0100	0.0276	0.0178	0.0119	0.0095	0.0082
0.2	0.0326	0.0198	0.0133	0.0117	0.0099	0.0356	0.0212	0.0147	0.0118	0.0109	0.0427	0.0367	0.0345	0.0339	0.0336
0.4	0.0316	0.0197	0.0139	0.0118	0.0107	0.0324	0.0199	0.0144	0.0119	0.0106	0.1015	0.0983	0.0976	0.0959	0.0962
0.6	0.0271	0.0174	0.0125	0.0103	0.0134	0.0266	0.0160	0.0119	0.0104	0.0087	0.2288	0.2300	0.2284	0.2275	0.2246
0.8	0.0655	0.0502	0.0472	0.0455	0.0439	0.0494	0.0361	0.0327	0.0344	0.0324	0.6065	0.6061	0.5963	0.5972	0.5870
β2-SAR-Poisson 1stStep-ML—W3						β2-SAR-Poisson 1stStep-OLS—W3					β2-Aspatial Poisson ML—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1356	0.0855	0.0592	0.0474	0.0413	0.1408	0.0865	0.0648	0.0509	0.0433	0.1224	0.0797	0.0540	0.0436	0.0394
0.2	0.1226	0.0759	0.0525	0.0437	0.0398	0.1295	0.0801	0.0573	0.0458	0.0423	0.1544	0.1242	0.1146	0.1126	0.1109
0.4	0.1143	0.0788	0.0645	0.0619	0.0558	0.1152	0.0718	0.0521	0.0456	0.0425	0.3083	0.2969	0.2915	0.2902	0.2899
0.6	0.1184	0.0851	0.0721	0.0582	0.0311	0.0898	0.0571	0.0394	0.0360	0.0341	0.6891	0.6764	0.6717	0.6715	0.6705
0.8	0.2134	0.1867	0.1864	0.1999	0.1949	0.1701	0.1297	0.1155	0.1217	0.1164	2.0317	2.0655	2.1007	2.0681	2.0419
Rho-SAR-Poisson 1stStep-ML—W3						Rho-SAR-Poisson 1stStep-OLS—W3
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1272	0.0784	0.0547	0.0436	0.0368	0.1414	0.0885	0.0621	0.0478	0.0441
0.2	0.0983	0.0554	0.0380	0.0320	0.0274	0.1116	0.0647	0.0482	0.0400	0.0368
0.4	0.0641	0.0377	0.0294	0.0268	0.0244	0.0717	0.0501	0.0419	0.0404	0.0382
0.6	0.0440	0.0308	0.0254	0.0203	0.0248	0.0422	0.0324	0.0308	0.0322	0.0307
0.8	0.0774	0.0542	0.0431	0.0417	0.0416	0.0429	0.036563	0.0342	0.0400	0.0387

Notes: (1) Bias is estimated as the average of 1000 simulation replicates of the difference between the parameter estimate and its true value. (2) RMSE is estimated as the square root of the sum between the square of the bias and the empirical variance of the estimated coefficient calculated after 1000 replicates. (3) SAR-Poisson 1stStep-ML is estimated using a two-step process. In the first step, the unobservable conditional mean, μ, is estimated using a PPML regression, and in the second step, the coefficients β1, β2, and ρ are also estimated using a PPML regression. (4) SAR-Poisson 1stStep-OLS is estimated using a two-step process. In the first step, the unobservable variable Wlog(μ) is estimated using an OLS regression of Wlog(y), adding an ad hoc constant (c = 1) when y = 0 and, in the second step, the coefficients β1, β2, and ρ are estimated using a PPML regression. (5) W1 and W3 are contiguity matrices created using the nearest neighbor criterion, where it is computationally defined that each unit will have seven units as neighbors for W1 and four units as neighbors for W2, which are the closest. W2 is created based on an inverse distance criterion, using the Euclidean distance between each unit.

of Appendix A.

Table A1. Bias: SAR-Poisson, SAR-LogLinear, and Aspatial ML Poisson with W1 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W1						β1-SAR-Poisson 1stStep-OLS—W1					β1-Aspatial Poisson ML—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0028	−0.0011	−0.0007	0.0000	0.0000	−0.0011	−0.0010	−0.0006	−0.0003	−0.0003	0.0003	−0.0004	−0.0005	0.0000	0.0000
0.2	−0.0013	0.0000	−0.0007	0.0001	−0.0004	−0.0047	−0.0041	−0.0044	0.0040	−0.0036	0.0311	0.0312	0.0312	0.0310	0.0308
0.4	−0.0006	0.0000	−0.0006	0.0002	0.0002	−0.0032	−0.0047	−0.0046	−0.0044	−0.0046	0.0891	0.0868	0.0872	0.0872	0.0870
0.6	0.0003	0.0007	0.0015	0.0012	0.0015	−0.0013	−0.0019	−0.0020	−0.0022	−0.0020	0.2021	0.1979	0.1988	0.1969	0.1971
0.8	0.0049	0.0046	0.0040	0.0032	0.0043	0.0021	0.0014	0.0017	0.0005	0.0004	0.4957	0.4846	0.4886	0.4771	0.4861
β2-SAR-Poisson 1stStep-ML—W1						β2-SAR-Poisson 1stStep-OLS—W1					β2-Aspatial Poisson ML—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0062	−0.0028	0.0026	0.0004	0.0003	−0.0105	−0.0021	−0.0010	−0.0015	−0.0015	−0.0064	−0.0004	−0.0003	−0.0004	−0.0007
0.2	−0.0016	0.0012	0.0014	−0.0004	0.0003	−0.0132	−0.0118	−0.0156	−0.0138	−0.0130	0.1018	0.1035	0.1028	0.1027	0.1036
0.4	0.0024	0.0009	0.0031	−0.0036	0.0043	−0.0143	−0.0182	−0.0188	−0.0191	−0.0179	0.2737	0.2768	0.2760	0.2789	0.2799
0.6	0.0017	0.0023	0.0040	−0.0001	0.0027	−0.0052	−0.0051	−0.0097	−0.0094	−0.0098	0.6201	0.6327	0.6391	0.6363	0.6383
0.8	−0.0015	−0.0022	−0.0058	0.0077	−0.0078	0.0072	0.0066	0.0065	0.0017	0.0017	1.7530	1.7635	1.7964	1.8182	1.8423
Rho-SAR-Poisson 1stStep-ML—W1						Rho-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0042	0.0015	0.0005	−0.0012	−0.0008	0.0037	0.0018	−0.0003	0.0018	0.0014
0.2	0.0012	−0.0025	0.0006	−0.0003	0.0003	0.0079	0.0105	0.0149	0.0131	0.0119
0.4	−0.0036	−0.0014	−0.0012	−0.0023	−0.0027	0.0169	0.0234	0.0246	0.0247	0.0249
0.6	−0.0002	−0.0004	−0.0018	−0.0016	−0.0010	0.0160	0.0181	0.0204	0.0208	0.0209
0.8	0.0047	0.0039	0.0048	0.0050	0.0050	0.0018	0.0034	0.0042	0.0061	0.0082

and

Table A3. Bias: SAR-Poisson, SAR-LogLinear, and Aspatial ML Poisson with W2 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W2						β1-SAR-Poisson 1stStep-OLS—W2					β1-Aspatial Poisson ML—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0006	−0.0010	0.0001	−0.0001	0.0003	−0.0001	−0.0010	0.0001	−0.0002	0.0002	−0.0008	−0.0004	0.0001	0.0000	0.0001
0.2	−0.0017	−0.0002	−0.0001	−0.0005	−0.0004	−0.0019	−0.0007	−0.0004	−0.0006	−0.0005	0.0290	0.0291	0.0285	0.0285	0.0286
0.4	0.0003	0.0002	−0.0006	0.0003	−0.0002	−0.0007	−0.0002	−0.0008	0.0001	−0.0003	0.0769	0.0761	0.0754	0.0750	0.0740
0.6	0.0012	−0.0005	−0.0003	0.0001	−0.0003	0.0008	−0.0004	−0.0002	0.0002	−0.0001	0.1711	0.1662	0.1609	0.1609	0.1592
0.8	0.0025	−0.0015	−0.0016	−0.0017	−0.0015	0.0003	0.0002	0.0000	0.0000	0.0001	0.4010	0.3873	0.3743	0.3630	0.3577
β2-SAR-Poisson 1stStep-ML—W2						β2-SAR-Poisson 1stStep-OLS—W2					β2- Aspatial Poisson ML—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0012	0.0010	−0.0027	0.0005	−0.0004	−0.0005	0.0010	−0.0029	−0.0009	−0.0005	−0.0785	−0.0014	−0.0011	0.0004	−0.0005
0.2	−0.0038	−0.0021	−0.0011	−0.0003	0.0001	−0.0058	−0.0035	−0.0019	−0.0003	−0.0003	−0.0520	0.0991	0.0998	0.1022	0.1009
0.4	0.0031	0.0020	−0.0003	−0.0003	0.0008	−0.0006	−0.0002	−0.0012	−0.0003	0.0002	−0.0219	0.2671	0.2663	0.2682	0.2670
0.6	0.0045	0.0018	0.0013	−0.0007	0.0005	−0.0007	0.0010	0.0015	0.0009	0.0008	0.0143	0.5905	0.6000	0.6034	0.6027
0.8	0.0397	0.0172	0.0109	−0.0069	0.0081	0.0015	0.0006	−0.0001	0.0000	0.0004	0.0364	1.6308	1.6558	1.6638	1.6922
Rho-SAR-Poisson 1stStep-ML—W2						Rho-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0152	−0.0020	0.0001	−0.0010	−0.0017	−0.0082	0.0007	0.0018	0.0010	−0.0006
0.2	−0.0033	0.0006	0.0007	0.0009	0.0015	−0.0122	−0.0108	−0.0115	−0.0121	−0.0119
0.4	−0.0038	−0.0009	0.0024	0.0002	0.0007	−0.0064	−0.0030	0.0000	−0.0026	−0.0016
0.6	0.0124	0.0072	0.0040	0.0021	0.0026	0.0126	0.0144	0.0152	0.0152	0.0154
0.8	0.0501	0.0341	0.0236	0.0166	0.0143	0.0047	0.0044	0.0045	0.0045	0.0043

in Appendix A show the results for the bias of the estimated coefficients for each estimation method, considering the spatial weighting matrix based, respectively, on the contiguity criterion, W1, and the inverse distance criterion, W2. Both SAR-PPML estimators show similar and quite satisfactory results, with the SAR-PPML 1stStep-ML presenting a lower absolute value for bias for lower and median levels of spatial dependence, while the SAR-PPML 1stStep-OLS appears to behave better for values of $\rho$ closer to 1. It is worth noting that both estimators have lower bias, as an absolute value, to estimate the coefficient of the continuous variable than that of the dummy variable. Note, also, that when $\rho$ increases, both estimators present a smaller absolute value of bias when using matrix W2 compared to matrix W1. Nevertheless, this difference is negligible, especially for a large n. Finally, the aspatial PPML estimator shows progressively worse results as $\rho$ increases, as expected, being slightly better than the SAR methods when there is no spatial dependence $(\rho =0)$.

Concerning the bias of the spatial autoregressive coefficient, $\rho$, globally, the SAR-PPML 1stStep-ML shows better performance than the remaining estimators, especially when n is large. However, for $\rho =0.8$ it shows a higher bias, in absolute value, particularly for the W2 matrix. Although slightly worse than the first, the SAR-PPML 1stStep-OLS presents satisfactory results, namely for high levels of spatial dependence. On the other hand, in general, the use of a special weighting matrix based on the inverse distance between locations produces higher bias when estimating the spatial autoregressive parameter.

Table A2. RMSE: SAR-Poisson, SAR-LogLinear, and Aspatial ML Poisson with W1 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W1						β1-SAR-Poisson 1stStep-OLS—W1					β1-Aspatial Poisson ML—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0242	0.0154	0.0102	0.0083	0.0069	0.0270	0.0163	0.0117	0.0090	0.0077	0.0171	0.0174	0.0122	0.0099	0.0084
0.2	0.0251	0.0152	0.0107	0.0087	0.0071	0.0271	0.0165	0.0119	0.0092	0.0083	0.0424	0.0353	0.0330	0.0321	0.0316
0.4	0.0223	0.0137	0.0102	0.0077	0.0070	0.0228	0.0141	0.0106	0.0089	0.0080	0.0946	0.0893	0.0885	0.0880	0.0877
0.6	0.0176	0.0106	0.0080	0.0065	0.0061	0.0181	0.0113	0.0079	0.0065	0.0056	0.2122	0.2027	0.2017	0.1989	0.1987
0.8	0.0161	0.0134	0.0124	0.0113	0.0106	0.0174	0.0155	0.0123	0.0120	0.0118	0.5380	0.5063	0.5029	0.4887	0.4954
β2-SAR-Poisson 1stStep-ML—W1						β2-SAR-Poisson 1stStep-OLS—W1					β2-Aspatial Poisson ML—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1034	0.0649	0.0454	0.0357	0.0314	0.1085	0.0648	0.0495	0.0382	0.0330	0.0716	0.0752	0.0564	0.0449	0.0382
0.2	0.0912	0.0567	0.0401	0.0333	0.0287	0.0999	0.0631	0.0437	0.0357	0.0316	0.1588	0.1263	0.1132	0.1089	0.1089
0.4	0.0811	0.0486	0.0353	0.0300	0.0261	0.0849	0.0531	0.0411	0.0354	0.0309	0.2967	0.2854	0.2800	0.2819	0.2821
0.6	0.0606	0.0395	0.0305	0.0277	0.0238	0.0643	0.0403	0.0292	0.0242	0.0215	0.6483	0.6454	0.6465	0.6423	0.6432
0.8	0.0602	0.0548	0.0558	0.0502	0.0476	0.0589	0.0487	0.0427	0.0404	0.0428	1.9053	1.8426	1.8515	1.8604	1.8784
Rho-SAR-Poisson 1stStep-ML—W1						Rho-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1112	0.0676	0.0474	0.0383	0.0320	0.1191	0.0734	0.0532	0.0398	0.0357
0.2	0.0792	0.0483	0.0330	0.0276	0.0239	0.0873	0.0542	0.0395	0.0318	0.0279
0.4	0.0512	0.0289	0.0214	0.0173	0.0149	0.0560	0.0390	0.0332	0.0307	0.0294
0.6	0.0237	0.0148	0.0122	0.0110	0.0097	0.0311	0.0244	0.0237	0.0230	0.0224
0.8	0.0142	0.0118	0.0106	0.0097	0.0094	0.0152	0.0146	0.0131	0.0137	0.0162

and

Table A4. RMSE: SAR-Poisson, SAR-LogLinear, and Aspatial ML Poisson with W2 for 1000 replicates.

β1-SAR-Poisson 1stStep-ML—W2						β1-SAR-Poisson 1stStep-OLS—W2					β1-Aspatial Poisson ML—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0387	0.0232	0.0165	0.0142	0.0114	0.0402	0.0236	0.0167	0.0143	0.0114	0.0276	0.0178	0.0120	0.0095	0.0082
0.2	0.0363	0.0207	0.0151	0.0121	0.0110	0.0360	0.0207	0.0151	0.0121	0.0110	0.0386	0.0327	0.0300	0.0296	0.0294
0.4	0.0304	0.0187	0.0136	0.0103	0.0088	0.0300	0.0186	0.0136	0.0103	0.0088	0.0820	0.0781	0.0763	0.0757	0.0745
0.6	0.0241	0.0134	0.0093	0.0076	0.0066	0.0233	0.0134	0.0093	0.0076	0.0065	0.1784	0.1694	0.1625	0.1621	0.1603
0.8	0.0305	0.0156	0.0123	0.0102	0.0092	0.0100	0.0056	0.0038	0.0031	0.0026	0.4256	0.4019	0.3816	0.3692	0.3626
β2-SAR-Poisson 1stStep-ML—W2						β2-SAR-Poisson 1stStep-OLS—W2					β2-Aspatial Poisson ML—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1487	0.0914	0.0669	0.0541	0.0488	0.1520	0.0924	0.0674	0.0543	0.0489	0.1224	0.0797	0.0551	0.0436	0.0394
0.2	0.1368	0.0843	0.0624	0.0495	0.0421	0.1376	0.0845	0.0625	0.0496	0.0422	0.1512	0.1238	0.1109	0.1095	0.1075
0.4	0.1182	0.0758	0.0515	0.0412	0.0362	0.1186	0.0761	0.0515	0.0411	0.0361	0.2866	0.2732	0.2714	0.2691	0.2709
0.6	0.0902	0.0527	0.0387	0.0302	0.0274	0.0894	0.0529	0.0385	0.0300	0.0273	0.6054	0.6064	0.6061	0.6050	0.6043
0.8	0.1310	0.0613	0.0402	0.0270	0.0277	0.0376	0.0216	0.0156	0.0126	0.0105	1.6870	1.6898	1.6781	1.7039	1.7125
Rho-SAR-Poisson 1stStep-ML—W2						Rho-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.2196	0.1283	0.0934	0.0803	0.0687	−0.0082	0.0007	0.0018	0.0010	−0.0006
0.2	0.1419	0.0860	0.0636	0.0484	0.0420	−0.0122	−0.0108	−0.0115	−0.0121	−0.0119
0.4	0.0814	0.0513	0.0360	0.0287	0.0245	−0.0064	−0.0030	0.0000	−0.0026	−0.0016
0.6	0.0571	0.0377	0.0264	0.0187	0.0157	0.0126	0.0144	0.0152	0.0152	0.0154
0.8	0.0840	0.0493	0.0290	0.0276	0.0245	0.0047	0.0044	0.0045	0.0045	0.0043

in Appendix A show the results referring to RMSE. From a general point of view and regarding ${\beta }_1$, the SAR-PPML 1stStep-ML presents the best results, particularly for W1. However, the SAR-PPML 1stStep-OLS produces a more desirable set of results for higher values of $\rho .$ In both estimators, it is noted that, as the spatial dependence and the sample size increase, the RMSE decreases. This result is only slightly altered when $\rho =0.8$. As expected, the aspatial ML estimator only shows satisfactory results when $\rho =0$. As for ${\beta }_2$, the conclusions are quite similar to ${\beta }_1$, with the disclaimer that the RMSEs for this coefficient are much higher, especially when the sample size is small. Estimations involving the W1 matrix have slightly better results. Lastly, the aspatial estimator is, again, quite far from the results of the other estimators.

Both SAR-PPML estimators present quite similar results regarding the RMSE values for the estimation of the coefficient of spatial dependence, $\rho$, with the SAR-PPML 1stStep-ML showing better results as the sample size increases. It is also important to note that the SAR-PPML 1stStep-ML exhibits a higher RMSE for matrix W2 for high levels of spatial dependence when compared to SAR-PPML 1stStep-OLS. However, in general, the use of W1 seems to trigger better results.

In summary, these results are in line with those obtained in other simulation studies such as [4,12,22,24,25], suggesting the following conclusions. First, the estimator SAR-PPML 1stStep-ML presents the best performance, except for high spatial dependence, when $\rho =0.8$. Keep in mind, however, that most empirical applications give low and median values for the spatial dependence parameter. Since this estimator does not rely on logarithmic transformation of the dependent variable and uses PPML regression instead of a loglinear estimation in the first step, this result seems to be in agreement with that found by [4]. Another interesting result is that there is a higher distortion for the estimated coefficient of the dummy variable compared to the estimated coefficient of the continuous variable, suggesting that the distribution of the explanatory variables can condition the performance of the estimators, a conclusion that is also made by [12]. Other similar conclusions between studies are the fact that the RMSE decreases as the spatial dependence and sample size increase, and that the spatial weighting matrix criteria influence the results. Several studies have already addressed this issue, such as [24], where the authors found that the RMSE of coefficient estimators appears to be generally higher for the spatial weighting matrix based on inverse distance, suggesting that the variance of the estimated coefficients may, somehow, be related to the density of the spatial weights matrix chosen. Another expected conclusion was the poor performance of the Aspatial PPML estimator in the presence of spatial dependence, which presented an accentuated upward bias for the coefficients of X₁ and X₂. This result is in agreement with [26], who found biased and inconsistent estimators when spatial dependence was not taken into consideration. In addition, it is interesting to note that the distortion of results is more significant for values of $\rho$ near 1, which is in line with the results of [22].

To assess the performance of both estimators under misspecification, a new design was considered where X₁ shows spatial dependence instead of being i.i.d. Therefore, X₁ was simulated according to the following spatial autoregressive process:

$$X_1={\left(I-0.5W_1\right)}^{-1}\mathsf{\epsilon } \mathrm{with} \epsilon \sim N\left(0,I\right).$$

The other variables were generated as before, with the coefficients retaining the same values. Estimation was implemented as if X₁ was i.i.d. (ignoring that it is spatially autocorrelated). The results obtained for 1000 replications, considering the spatial weighting matrix W1, are included in

Table A7. Bias: SAR-Poisson. SAR-LogLinear with W1 for 1000 replicates; X₁ with spatial dependence.

β1-SAR-Poisson 1stStep-ML—W1						β1-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0123	−0.0068	−0.0035	−0.0022	−0.0031	−0.0077	−0.0056	−0.0062	−0.0018	−0.0028
0.2	−0.0024	−0.0021	0.0001	0.0003	0.0008	0.0044	0.0032	0.0056	0.0066	0.0068
0.4	−0.0024	−0.0016	−0.0009	−0.0004	0.0004	0.0077	0.0086	0.0102	0.0104	0.0113
0.6	−0.0029	−0.0009	−0.0033	−0.0019	−0.0008	0.0086	0.0121	0.0107	0.0126	0.0136
0.8	−0.0093	−0.0164	−0.0191	−0.0220	−0.0226	0.0097	0.0089	0.0100	0.0107	0.0097
β2-SAR-Poisson 1stStep-ML—W1						β2-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0205	−0.0086	−0.0054	−0.0020	−0.0030	−0.0187	−0.0094	−0.0031	−0.0021	−0.0038
0.2	−0.0123	−0.0024	−0.0039	−0.0028	−0.0011	−0.0364	−0.0302	−0.0326	−0.0303	−0.0289
0.4	−0.0073	0.0012	−0.0009	0.0018	0.0017	−0.0583	−0.0553	−0.0573	−0.0554	−0.0564
0.6	0.0036	0.0033	−0.0014	0.0022	0.0021	−0.0694	−0.0757	−0.0789	−0.0795	−0.0788
0.8	0.0013	0.0045	0.0016	−0.0018	0.0001	−0.0532	−0.0631	−0.0686	−0.0789	−0.0778
Rho-SAR-Poisson 1stStep-ML—W1						Rho-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0264	0.0079	0.0035	0.0021	0.0045	0.0206	0.0124	0.0071	0.0029	0.0063
0.2	−0.0012	−0.0032	0.0034	0.0028	0.0001	−0.0056	0.0035	0.0094	0.0066	0.0051
0.4	0.0073	0.0052	0.0075	0.0062	0.0050	0.0204	0.0365	0.0419	0.0428	0.0435
0.6	0.0173	0.0088	0.0106	0.0081	0.0060	0.0586	0.0807	0.0868	0.0893	0.0895
0.8	0.0400	0.0306	0.0323	0.0329	0.0314	0.0555	0.0766	0.0837	0.0926	0.0920

for bias and in

Table A8. RMSE: SAR-Poisson. SAR-LogLinear with W1 for 1000 replicates; X₁ with spatial dependence.

β1-SAR-Poisson 1stStep-ML—W1						β1-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1142	0.0703	0.0477	0.0399	0.0359	0.1140	0.0702	0.0472	0.0395	0.0353
0.2	0.1068	0.0700	0.0477	0.0404	0.0348	0.1043	0.0663	0.0462	0.0391	0.0340
0.4	0.1103	0.0666	0.0462	0.0394	0.0325	0.1051	0.0636	0.0996	0.0824	0.0324
0.6	0.0950	0.0604	0.0408	0.0326	0.0302	0.0910	0.0586	0.0396	0.0333	0.0320
0.8	0.0763	0.0487	0.0408	0.0385	0.0373	0.0707	0.0462	0.0332	0.0291	0.0263
β2-SAR-Poisson 1stStep-ML—W1						β2-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1740	0.1077	0.0767	0.0637	0.0536	0.2156	0.1323	0.0948	0.0771	0.0659
0.2	0.1635	0.1042	0.0727	0.0617	0.0520	0.2057	0.1274	0.0920	0.0794	0.0681
0.4	0.1638	0.1003	0.0720	0.0573	0.0498	0.2003	0.1274	0.0449	0.0385	0.0769
0.6	0.1548	0.0956	0.0695	0.0586	0.0476	0.1822	0.1257	0.1042	0.0958	0.0892
0.8	0.1548	0.1073	0.1017	0.0969	0.0900	0.1582	0.1117	0.0959	0.0962	0.0900
Rho-SAR-Poisson 1stStep-ML—W1						Rho-SAR-Poisson 1stStep-OLS—W1
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.3733	0.2327	0.1647	0.1259	0.1109	0.4042	0.2554	0.1768	0.1393	0.1236
0.2	0.2934	0.1744	0.1175	0.1017	0.0818	0.3358	0.1986	0.1342	0.1173	0.0986
0.4	0.2205	0.1202	0.0830	0.0675	0.0575	0.2615	0.1477	0.1102	0.0904	0.0795
0.6	0.1368	0.0733	0.0510	0.0423	0.0357	0.1610	0.1142	0.0999	0.0968	0.0948
0.8	0.0873	0.0517	0.0451	0.0443	0.0417	0.1109	0.0927	0.0926	0.0998	0.0974

for RMSE, while

Table A9. Bias: SAR-Poisson. SAR-LogLinear with W2 for 1000 replicates; X₁ with spatial dependence.

β1-SAR-Poisson 1stStep-ML—W2						β1-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0123	−0.0024	−0.0026	−0.0002	−0.0015	−0.0064	−0.0010	−0.0022	0.0016	−0.0013
0.2	−0.0040	−0.0032	0.0005	−0.0001	0.0007	−0.0016	−0.0023	−0.0043	0.0002	0.0008
0.4	−0.0042	0.0027	−0.0006	−0.0008	0.0006	−0.0033	0.0034	−0.0002	−0.0006	0.0008
0.6	−0.0016	0.0016	0.0002	0.0009	0.0005	−0.0014	0.0020	0.0003	0.0009	0.0005
0.8	0.0001	0.0049	0.0029	0.0013	0.0005	−0.0007	0.0012	0.0013	−0.0005	−0.0007
β2-SAR-Poisson 1stStep-ML—W2						β2-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	−0.0205	−0.0105	−0.0039	−0.0002	−0.0039	−0.0101	−0.0040	−0.0026	0.0002	−0.0034
0.2	−0.0122	−0.0033	−0.0042	0.0030	−0.0031	−0.0088	−0.0045	0.0009	0.0025	−0.0038
0.4	−0.0128	−0.0010	−0.0044	0.0011	−0.0026	−0.0116	−0.0055	−0.0058	−0.0005	−0.0036
0.6	−0.0010	−0.0026	−0.0007	−0.0009	0.0013	−0.0064	−0.0063	−0.0018	−0.0018	0.0004
0.8	−0.0030	−0.0034	−0.0002	−0.0011	−0.0009	−0.0024	−0.0001	0.0024	0.0014	0.0006
Rho-SAR-Poisson 1stStep-ML—W2						Rho-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.0264	−0.0183	−0.0201	−0.0258	−0.0070	−0.0024	−0.0045	0.0007	−0.0084	0.0027
0.2	0.0155	−0.0226	−0.0082	−0.0167	−0.0048	−0.0681	−0.0597	−0.0580	−0.0656	−0.0571
0.4	0.0100	−0.0100	0.0002	−0.0058	0.0005	−0.0904	−0.0728	−0.0654	−0.0701	−0.0653
0.6	0.0071	0.0018	0.0012	0.0023	−0.0018	−0.0497	−0.0225	−0.0170	−0.0142	−0.0152
0.8	0.0150	0.0136	0.0056	0.0066	0.0044	0.0135	0.0472	0.0570	0.0597	0.0620

and

Table A10. RMSE: SAR-Poisson. SAR-LogLinear with W2 for 1000 replicates; X₁ with spatial dependence.

β1-SAR-Poisson 1stStep-ML—W2						β1-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1142	0.0709	0.0503	0.0418	0.0353	0.1117	0.0717	0.0504	0.0421	0.0354
0.2	0.1113	0.0656	0.0472	0.0395	0.0337	0.1116	0.0657	0.0473	0.0397	0.0337
0.4	0.1072	0.0668	0.0445	0.0393	0.0330	0.1063	0.0666	0.0444	0.0390	0.0330
0.6	0.0970	0.0594	0.0407	0.0345	0.0309	0.0964	0.0591	0.0406	0.0343	0.0308
0.8	0.0774	0.0457	0.0308	0.0242	0.0214	0.0770	0.0441	0.0306	0.0241	0.0213
β2-SAR-Poisson 1stStep-ML—W2						β2-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.1740	0.1337	0.0986	0.0834	0.0710	0.2279	0.1466	0.1007	0.0854	0.0723
0.2	0.2144	0.1343	0.0949	0.0797	0.0690	0.2264	0.1397	0.0962	0.0811	0.0698
0.4	0.2003	0.1293	0.0908	0.0753	0.0660	0.2143	0.1322	0.0915	0.0758	0.0661
0.6	0.1901	0.1190	0.0861	0.0670	0.0612	0.1950	0.1185	0.0854	0.0671	0.0612
0.8	0.1522	0.0970	0.0624	0.0518	0.0443	0.1542	0.0901	0.0621	0.0518	0.0442
Rho-SAR-Poisson 1stStep-ML—W2						Rho-SAR-Poisson 1stStep-OLS—W2
Rho/n	100	250	500	750	1000	100	250	500	750	1000
0.0	0.3733	0.4956	0.3579	0.3020	0.2538	0.4872	0.3062	0.2027	0.1741	0.1505
0.2	0.3981	0.3621	0.2393	0.2000	0.1718	0.4180	0.2643	0.1808	0.1573	0.1370
0.4	0.4164	0.2356	0.1561	0.1264	0.1106	0.3929	0.2239	0.1568	0.1347	0.1222
0.6	0.3435	0.1368	0.0897	0.0667	0.0619	0.2933	0.1541	0.1102	0.0852	0.0786
0.8	0.1594	0.0571	0.0338	0.0294	0.0250	0.2220	0.0781	0.0662	0.0643	0.0650

show, respectively, the bias and RMSE when the spatial weighted matrix is W2. Results show that ignoring spatial autocorrelation in the explanatory variable leads to noticeably higher bias and RMSE in the estimation of all parameters, especially in the estimation of the spatial autocorrelation coefficient. Both estimators show similar performance in estimating the coefficient of X₁, whether the spatial matrix is based on the nearest neighbor criterion (W1) or the inverse distance (W2). The new estimator introduced, SAR-PPML 1stStep-ML, shows better performance than SAR-PPML 1stStep-OLS for the coefficient of X₂ when the spatial matrix is W1. The improvement in performance of SAR-PPML 1stStep-ML over the SAR-PPML 1stStep-OLS is especially visible in the estimation of $\rho$ for both spatial weighting matrices.

4. Empirical Application

This section illustrates the usefulness of the SAR-E regression introduced in Section 2 by an empirical example that estimates a knowledge production function to explain the creation of innovation in European regions. For the sake of comparison, the estimator of [12] is also calculated.

Following the arguments of [27], the number of patents in a given region per million inhabitants (Pat) is used as a proxy for knowledge creation. See also [28] for a discussion on measuring innovation. The equation to be estimated is

$$E\left(Pa{t}_i|x_i\right)={\mu }_i=\exp \left(\rho {\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{w}_{ij}}\log ({\mu }_j)+x_i\mathsf{\beta }\right),$$

(7)

where $x_i$ is a vector with explanatory variables that will be introduced in Section 4.1, $\mathsf{\beta }$ is a vector of unknown coefficients to be estimated, $\rho$ is the unknown spatial autocorrelation coefficient, and $w_{ij}$ are the elements of a spatial weighting matrix. In this empirical application, the spatial weighting matrix was calculated based on a queen contiguity criterion and is row normalized. All estimations were conducted using R [20]. The exploratory data analysis was performed according to QGIS [29] and GeoDa [30].

4.1. Data and Variables

The data were collected from Eurostat regional statistics. They contain data on 234 NUTS II regions from 24 European countries, of which 22 belong to the European Union, with the addition of the United Kingdom and Norway. NUTS is a nomenclature of territorial units for statistics developed and regulated by the European Union, defining a hierarchical system of regions with three different levels. At the top of the hierarchy are the NUTS 0 regions, referring to countries. The next level is NUTS 1, representing major socioeconomic regions within countries, followed by NUTS 2 regions, which are subdivisions of NUTS 1, and NUTS 3 regions, which are subdivisions of NUTS 2. All data refer to 2012. Regions with no neighbors were excluded (like Portuguese and French islands). Finally, NUTS II London (UK) and Centre (France), were discarded for incongruity of data. The list of countries in the database is in Appendix B.

The description of variables used in this study, together with the expected outcome of the estimated associated coefficient, can be found in

Table 1. Variable definitions and expected sign.

Variables	Abbrev.	Unit	Expected Outcome
The number of patents registered (dependent variable)	Pat	Unit per million inhabitants	−
Intramural Expenditure on R&D by private business	R&D_B	Euros per inhabitant	+
Intramural Expenditure on R&D by the government	R&D_G	Euros per inhabitant	Ambiguous
Intramural Expenditure on R&D by universities	R&D_U	Euros per inhabitant	Ambiguous
Total R&D personnel and researchers in private business (no. of full-time workers)	Pers_B	-	+
Total R&D personnel and researchers in the government (no. of full-time workers)	Pers_G	-	Ambiguous
Total R&D personnel and researchers in universities (no. of full-time workers)	Pers_U	-	Ambiguous
% Population aged 25–64 with Bachelor’s degree	Educ	Percentage	+
Population	Pop	Number of inhabitants	−
GDP per capita	GDP	Thousand euros per capita	+
Tuberculosis mortality	Mort	Rate per 100,000 inhabitants	−

.

Since [31] introduced the knowledge production function, the use of variables related to R&D has become normal when modeling the creation of innovation. Following [27,28,32,33,34] different impacts on the creation of knowledge from expenditure and human resources in R&D were considered according to its source (from the private and business sector, from government, or from universities). It is expected that more R&D expenditure, as well as more full-time R&D employees, will trigger an increase in knowledge creation. Therefore, the expected outcome of the estimated coefficients related to these variables should be positive. However, the literature suggests that this happens only for the R&D resources of the private sector. For both the public sector and universities, the effect of those variables often appears to be negative or statistically negligible (see, e.g., [27,33,35]). This behavior may be explained in the case of universities by the fact that their main contribution to knowledge creation arises in the form of scientific articles and not patents, while for the public sector it may be due to a certain inefficiency of public institutions in the production of knowledge (see [33,34]).

Three variables aiming to capture the effect of the “innovative environment” are considered. The first is the percentage of graduates in the population between 25 and 65 years old, proxying the level of education of the population in the region. The second is GDP per capita, which proxies the technological sophistication and the size of the economy. Finally, the third is the tuberculosis mortality rate, considered as a proxy for the level of poverty of the inhabitants, as several studies relate tuberculosis with poverty (see, for example, [33]). It is expected that a better socioeconomic environment will boost innovation (as in [33]). In addition, the number of inhabitants was defined as the control variable.

Table 2. Descriptive statistics.

Variables	N	Mean	Std Dev	Min	Max
Pat	234	89.171	106.045	0.000	590
R&D_B	234	318.248	382.444	0.000	2441.700
R&D_G	234	59.819	87.684	0.000	480.600
R&D_U	234	135.917	152.388	0.000	891.700
Pers_B	234	5744.342	9291.554	0.000	97,982.000
Pers_G	234	1467.979	2628.740	0.000	17,934.000
Pers_U	234	3294.923	3558.347	0.000	34,836.000
Educ	234	27.334	8.700	11.200	50.100
Pop	234	1,982,780.5	1,563,839.6	126,620.0	11,898,502.0
GDP	234	26.922	13.874	3.561	84.047
Mort	234	1.009	1.375	0.100	8.800

includes the descriptive statistics of the variables used in this study. Additionally, we note that 6% of the regions in the sample registered no patents.

The correlation matrix of these variables is shown in

Table A11. Correlation matrix of the variables.

	Pat	R&D_B	R&D_G	R&D_U	Pers_B	Pers_G	Pers_U	Educ	Pop	GDP	Mort
Pat	1.000
R&D_B	0.717	1.000
R&D_G	0.301	0.438	1.000
R&D_U	0.390	0.533	0.510	1.000
Pers_B	0.474	0.601	0.332	0.211	1.000
Pers_G	0.153	0.215	0.591	0.117	0.622	1.000
Pers_U	0.164	0.286	0.329	0.261	0.746	0.674	1.000
Educ	0.263	0.450	0.408	0.440	0.296	0.253	0.355	1.000
Pop	0.056	0.082	0.119	−0.088	0.663	0.651	0.775	0.010	1.000
GDP	0.573	0.639	0.519	0.642	0.368	0.148	0.226	0.559	−0.043	1.000
Mort	−0.285	−0.248	−0.176	−0.266	−0.141	−0.017	−0.099	−0.275	0.082	−0.466	1.000

of Appendix B. Pairwise correlations between explanatory variables do not exceed the common threshold of 0.8, as recommended in [36], which leads us to not anticipate collinearity problems in estimation.

4.2. Exploratory Spatial Analysis

Analyzing the spatial distribution map of the variable Pat per quartile in Figure A1 in Appendix B, we see the existence of a cluster effect, with a concentration of patenting taking place in Central Europe, Southern England, and Scandinavia, with the number of new patents in Southern and Eastern Europe being modest. On the other hand, Moran’s I test for spatial autocorrelation, applied to patents, shows a value for the test statistic equal to 0.6045, with a p-value equal to 0.000, denoting evidence of positive spatial dependence. This conclusion is supported by the Moran diagram (Figure A2 in Appendix B). Analyzing the latter, it is worth noting that most of the observations are in the 1st and 3rd quadrant, and therefore, the majority of regions with a higher (smaller) number of new patents have neighboring regions where this number is also higher (smaller). Analyzing the LISA3 indicators in Figure A3 of Appendix B, we see two high patent clusters in Central Europe and Scandinavia, together with low patent clusters in the Iberian Peninsula and Eastern Europe. It is also possible to identify two other clusters where patenting tends to be low, northern Britain and southern Italy. Figure A4 of Appendix B shows the LISA significance map, inferring that the results are more significant concerning the Central European, Iberian Peninsula, and Eastern Europe clusters.

4.3. Estimation Results

Equation (7) is estimated with the introduced SAR-PPML 1stStep-ML estimator. For comparative purposes, the results obtained with the alternative estimator from [12], the SAR-PPML 1stStep-OLS, are presented as well.

In SAR-PPML 1stStep-OLS an ad hoc constant (c = 1) was added when patents were 0.

Table 3. SAR PPML coefficient estimates of the knowledge production function.

	SAR-PPML 1stStep-ML		SAR-PPML 1stStep-OLS
Variable	Coefficients	Bootstrap SE	Coefficients	Bootstrap SE
$\rho$	6.81 × 10−1 ***	0.06838	6.19 × 10−1 ***	0.07821
R&D_B	8.91 × 10−4 ***	0.00034	9.18 × 10−4 ***	0.00035
R&D_G	−2.15 × 10−3 *	0.00130	−2.01 × 10−3	0.00157
R&D_U	−3.21 × 10−4	0.00079	−5.85 × 10−4	0.00083
Pers_B	−1.33 × 10−5	0.00003	−1.59 × 10−5	0.00003
Pers_G	2.75 × 10−5	0.00006	3.43 × 10−5	0.00006
Pers_U	5.07 × 10−5	0.00005	1.86 × 10−5	0.00005
Educ	2.58 × 10−4	0.01165	7.55 × 10−5	0.01279
Pop	−3.21 × 10−9	9.62 × 10−8	6.10 × 10−8	1.14 × 10−7
GDP	3.81 × 10−2 ***	0.01003	4.42 × 10−2 ***	0.01127
Mort	−1.95 × 10−1 **	0.09624	−6.72 × 10−2	0.09071
Log Likelihood	−6557.154		−7704.048
N	234		234

Notes: Standard errors were computed using the bootstrap method. Significance levels: * 10%, ** 5%, *** 1%.

includes the estimates for coefficients of the knowledge production function, together with their bootstrap standard errors. Bear in mind that, because both estimators use an explanatory variable that is the result of a fit obtained in the first step, the usual standard errors are not valid. The introduced PPML 1stStep-ML behaves better in terms of goodness of fit, with the value for the loglikelihood being noticeably higher.

For both estimations, the coefficient related to the spatially lagged variable is positive and significant (p-value < 0.01), thus inferring that there is a clear positive spatial dependence between regions regarding innovation creation, which matches the results of [34,37].

As for the remaining explanatory variables, the variable R&D_B appears to be significant at 1% in all estimates. In contrast, R&D_U is not significant, which can be explained by the fact that university contributions are mostly in the form of scientific articles and not patents, as mentioned before. On the other hand, R&D_G is significant at 10% in SAR-PPML 1stSep-ML; however, it presents a negative sign. These results converge with those of [32,33,34], which also disclose evidence of inefficiency in the use of R&D resources of the public sector. In addition, these authors also conclude that R&D expenditures in the private sector are more important to trigger knowledge creation than those from the public sector or universities.

The variables related to the “Innovative Environment,” Educ and Pop, are not statistically significant, while GDP is statistically significant at 1% in both estimations, with a positive sign. Finally, the mortality rate appears significant at 5% only in the SAR-PPML 1stStep-ML, showing a negative sign. These results are in line with expectations, as greater technological sophistication is generally associated with lower levels of poverty and higher quality of life, which fosters the growth of innovation in a region. These results corroborate studies such as [28,33], the authors of which conclude that an “innovative environment” is important for increasing knowledge creation.

Given the nonlinearity of the model, it is through the average partial effects (APE) that it is possible to quantify the impact of the variation of the explanatory variables on the dependent variable, ceteris paribus. Given the autoregressive structure of the model, it is possible to measure the indirect partial effects, that is, spatial externalities, together with the direct ones. These are included in

Table 4. Average direct partial effects and spatial spillovers.

	SAR-PPML 1stStep-ML			SAR-PPML 1stStep-OLS
Variable	Direct	ASpill-in	ASpill-out	Direct	ASpill-in	ASpill-out
R&D_B	0.0934	0.1743	0.1683	0.0878	0.1126	0.1093
R&D_G	−0.2250	−0.4200	−0.4054	−0.1921	−0.2464	−0.2391
R&D_U	−0.0336	−0.0628	−0.0606	−0.0560	−0.0718	−0.0697
Pers_B	−0.0014	−0.0026	−0.0025	−0.0015	−0.0020	−0.0019
Pers_G	0.0029	0.0054	0.0052	0.0033	0.0042	0.0041
Pers_U	0.0053	0.0099	0.0096	0.0018	0.0023	0.0022
Educ	0.0270	0.0504	0.0487	0.7221	0.9258	0.8984
Pop	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
GDP	3.9933	7.4536	7.1953	4.2288	5.4222	5.2617
Mort	−20.4060	−38.0886	−36.7686	−6.4270	−8.2409	−7.9969

.

Concerning the average direct effects, and the SAR-PPML 1stStep-ML, an increase of 1 percentage point in the tuberculosis mortality rate in the region results, on average, in a drop of 20.4 patents per million inhabitants in the same region. On the other hand, an increase in GDP per capita of just 100 euros in the region may trigger an increase, on average, of 0.4 patents per million inhabitants in the same region.

Regarding the variables of expenditure on R&D, these can present the most interesting results for economic decision makers. An increase of 10 euros per capita in a region in public R&D entities means, on average, a decrease of 2.25 patents in that region per million inhabitants. Keep in mind that the respective coefficient estimate is statistically significant only at 10%. Now, given the inefficiency inferred there, a policy maker may transfer the financial resources of these institutions to private R&D companies, since these, for each increase of 10 euros per capita in R&D expenses, trigger an increase of approximately one patent per million inhabitants in the same region. The spatial distribution map of the SAR-Poisson 1stStep-ML direct partial effect (DPE) per quartile related to the variable R&D_B is represented in Figure A5 in Appendix B. It is clear that the regions with the most efficient companies for transforming R&D expenses into patents are located in Central Europe, southern Great Britain, and Scandinavia. Therefore, regions in Eastern and Southern Europe should initiate reforms in the private R&D creation system, seeking an increase in efficiency. These reforms require the recruitment of more qualified personnel and investment in more sophisticated technology.

As for the indirect effects, the variables related to the “innovative environment” show higher indirect effects in absolute value than direct, showing that not only the socioeconomic situation of the region is central to the creation of knowledge, but also the interregional environment.

Concerning the R&D expenditure variables, investment in government R&D institutions also does not benefit neighboring regions in the knowledge creation process, since both the spill-out effect and spill-in effect are negative. On the other hand, investment in private R&D in one region will have a positive impact in neighboring regions: a variation of 10 euros per inhabitant in private R&D expenditure in all neighboring regions of i results in an increase of 1.74 new patents in region i. Conversely, an increase of 10 euros per inhabitant in expenditure on private R&D in region i results in an increase, on average, of 1.68 in all neighboring regions. These facts highlight the presence of knowledge spillovers between regions. Figure A6 and Figure A7 in Appendix B represent the spatial distribution map of spill-in and spill-out effects per quartile obtained with the SAR-PPML 1stStep-ML, respectively, of the variable R&D_B. It can be concluded that, in addition to the Central European cluster, which shows a strong relationship in the creation of knowledge, regions in Southern and Eastern Europe, as well as some regions in Southern England, have a remarkable capacity for absorbing innovation. Regarding the spill-out effects, the Central European and Scandinavian clusters are the biggest “exporters” of knowledge spillovers. Interestingly, some regions that present a lower DPE with the investment in private R&D, as is the case in Eastern Europe and the North of the United Kingdom, present higher values of spill-out and spill-in. Therefore, one may conclude that, despite having a lower capacity for innovation, these regions have a strong interconnection, which leads to high levels of knowledge spillover. This can be explained by a possible commitment of companies to strong interregional cooperation links, so that investment in one company is positively reflected in the others. These links can be a strategy to overcome the difficulty of competing solo against regions with high numbers of patents. Therefore, political and economic decision-makers in regions with lower patent capabilities should create incentives for the creation of knowledge-sharing networks, thus enabling increased competitiveness.

5. Conclusions

Many applications in spatial econometrics concern the modeling of count outcomes or other nonnegative variables. This work proposes modeling such variables by a spatial autoregressive exponential (SAR-E) regression instead of using SAR loglinear models, in line with the reasoning of [4] in the context of cross-sectional data. A two-step PPML procedure for the SAR-E model is suggested that circumvents the problem of dealing with zeros. A simulation study verifies that the introduced estimator shows better performance than the previous estimation procedures, independently of the sample size, especially when the autoregressive coefficient is not close to 1, which is the case for many applications with economic data.

The usefulness of the proposed approach is illustrated in an empirical application to analyze the main determinants of knowledge creation and to quantify the spatial knowledge spillovers across different European NUTS II regions. There, evidence of the spatial dependence on the creation of innovation in Europe is found. In addition, it is inferred that social and economic factors determine the creation of knowledge, as is the case with quality-of-life standards and technological sophistication. It also appears that public R&D institutions are inefficient, unlike private institutions, with the latter being the major promoters of innovation creation in the analyzed regions. It is also inferred that an increase in R&D expenditure by private institutions positively influences the creation of innovation in neighboring regions. Given these results, it is possible to conclude that regions with low levels of knowledge creation try to overcome this limitation by strengthening relationships with neighboring regions, thereby increasing the absorptive capacity for innovation and creating strong clusters of knowledge sharing.

In the empirical study, there were some noteworthy differences in the results obtained with the method introduced in this paper to estimate the SAR-E, two-step PPML, and the existing method of [12] concerning both the statistical significance and magnitude of some coefficient estimates, namely the autoregressive parameter. Differences in the latter explain the visible differences in the indirect effects of variables, the spill-in and spill-out spillovers, obtained with the two methods, where spillovers obtained with the proposed method are higher in absolute value. These differences are not unexpected because the data for the response variable show a clear percentage of zeros. Given the results of the simulation study and the fact that the method introduced here is better able to handle zeros in the dependent variable, it is expected that it will deliver better estimates.