spsur : An R Package for Dealing with Spatial Seemingly Unrelated Regression Models

Spatial seemingly unrelated regression (spatial SUR) models are a useful multiequational econometric specification to simultaneously incorporate spatial effects and correlated error terms across equations. The purpose of the spsur R package is to supply a complete set of functions to test for spatial structures in the residual of a SUR model; to estimate the most popular specifications by applying different methods and test for linear restrictions on the parameters. The package also facilitates the estimation of so-called spatial impacts, conveniently adapted to a SUR framework. The package includes functions to simulate datasets with the features decided by the user, which may be useful in teaching activities or in more general research projects. The article concludes with a real data application showing the potential that spsur has to examine the relation of individual mobility over geographic areas and the incidence of COVID-19 in Spain during the first lockdown.


Introduction
Seemingly unrelated regression models (SUR from now on) are a type of multiequational econometric formulation that gained popularity with the publication of the seminal paper by Zellner (1962). The SUR framework assures increased efficiency by estimating the system jointly rather than by processing each equation separately. These gains explain the relevant position that SUR models occupy in the applied research agenda. However, spatial econometrics is a particular exception: the applied literature in this area is limited, in spite of the extensive methodological research on spatial SUR models (Anselin 1988a,b;Wang and Kockelman 2007;Baltagi and Bresson 2011;Baltagi and Pirotte 2011;Mur, López, and Herrera 2010;Anselin 2016;López, Mínguez, and Mur 2020). Anselin (1988a) introduced the term spatial SUR for the first time in reference to a model made of "an equation for each time period, which is estimated for a cross section of spatial units" (page 141). Anselin's approach focuses on the problem of serial dependence in the errors, where each equation stands for a time period. In this work he did not explore estimation techniques, although the likelihood function, the information matrix, and some spatial autocorrelation tests were presented. From this seminal reference in Anselin's book (Anselin 1988a), a relevant number of methodological contributions have been made present, including developing spatial autocorrelation tests and implementing estimation methods.
With respect to the tests of spatial autocorrelation developed in the residual, several papers have focused on this topic. For example, Anselin (1988b) develops a test based on the Lagrange multiplier (LM) principle for spatial error SUR models. Mur et al. (2010) design robust and non-robust LM tests including an extensive Monte Carlo exercise for small sample sizes. Baltagi and Bresson (2011) propose joint and conditional LM tests for spatial autocorrelation and random effects for a spatial SUR panel model, while López, Mur, and Angulo (2014) address the problem of selecting the best specification for spatial SUR models using a battery of LM tests.
Several papers have focused on estimation methods. The maximum likelihood (ML) approach under the normality assumption is the most frequently used method, though other alternatives have been proposed. The paper by Wang and Kockelman (2007) makes a significant contribution to this literature, implementing a three-step estimation method combining feasible generalized least squares (FGLS) and ML, which allows for heterogeneity among the individuals by introducing unobserved random effects "à la KKP" (Kapoor, Kelejian, and Prucha 2007). Baltagi and Bresson (2011) extend the results by Wang and Kockelman (2007) by developing a collection of misspecification tests under an ML framework, whereas Baltagi and Pirotte (2011) introduce the generalized method of moments (GMM) in a spatial SUR with spatial autocorrelated errors, also "à la KKP". Anselin (1988a, page 146) describes the three-least-squares (3SLS) estimation method based on instrumental variables (IV) for spatial lag SUR models and Anselin (2016) points out the choice of the instruments. López et al. (2020) develop an extensive Monte Carlo exercise to compare ML and 3SLS (based on IV) methods for spatial lag and spatial Durbin SUR models. The method of generalized spatial three-stage least squares (GS3SLS), developed by Kelejian and Prucha (2004) in the more general framework of spatial simultaneous systems, has been used to estimate spatial SUR Models. The GS3SLS method is based on an iterative process that uses IV and GMM (Kelejian and Prucha 1999). Finally, Kakamu, Polasek, and Wago (2012) develop a Bayesian approach for the estimation of spatial SUR models using a Markov chain Monte Carlo (MCMC) method.
From the discussion above, it is clear that there is growing interest in Seemingly Unrelated Regressions in a spatial context, although the lack of specific software has delayed its adoption in applied research. Empirical applications are still hindered by the lack of readily available user-friendly software. Our purpose in this paper is to introduce the R package called spsur to fill this gap. The spsur R package offers accessible tools for those researchers who need to incorporate multiequational models into their research when spatial effects and correlation between residuals are present. The spsur user has a powerful tool to test for spatial autocorrelation in the residuals of a classic SUR model. This package can estimate spatial SUR models by ML and 3SLS 1 , and it uses the Wald statistic to test hypotheses about the stability of parameters, including the parameters of spatial dependence. It can even estimate a restricted model when the parameters are equal in several equations. Direct, indirect, and total impacts can be obtained for all the specifications. Finally, it can estimate uniequational spatial models if it is necessary to compare the estimation of each equation with the result of the multiequational model.
The structure of this paper is the following. Section 2 introduces the notation and presents the main methodology developed on spatial SUR models in the literature. Section 3 presents the structure of the package and some of the main functions included in spsur. It focuses on the estimation methods, the misspecification tests, and the computation of spatial impacts (direct, indirect, and total). Section 4 introduces an additional functionality that facilitates simulation experiments that could be helpful in research and/or teaching activities. Section 5 presents an illustration of the spsur package using real data. Finally, Section 6 concludes with a brief summary and prospects for the future.

Spatial SUR models
The SUR model consists of several equations, possibly with different regressors, where the error terms are correlated. We focus on Anselin's case with a single equation, G = 1, several time periods, T m, and N individuals in each cross section. Usually, N > T m. Originally, the spatial SUR specification assumes that the cross-sectional dimension N is large and that the time dimension T m is finite. For most regional studies, this assumption of short panels holds true. Following the usual notation, we consider that each equation refers to a time period although the model can be used to estimate different variables in the same cross section.
The baseline model without spatial effects is the following: where y t and u t are N × 1 vectors, X t is an N × p t matrix, with p t the number of regressors that appear in the t-th equation, and β t is the corresponding (p t × 1) vector of coefficients. This model is spatially independent. Therefore, following the terminology in López et al. (2014), we call it SUR-SIM.
The serial dependence in the errors of (1) is not parameterized but estimated in the T m × T m covariance matrix, Σ = (σ ts ; t, s = 1, . . . , T m ).
When we are dealing with spatial data, a frequent problem is spatial dependence, an unobservable feature that can appear in different ways. A general SUR model which includes most of the possible spatial effects of interest in applied research is the following (Elhorst 2014): where X * t is the matrix of regressors excluding the intercept term; W y t , W x t , and W u t are N × N weighting matrices. The model in (2) is referred to as a general nesting model (SUR-GNM, from now on). It is normal to assume that the weighting matrices are the same, Moreover, the user should supply the exogenous W matrix based on a priori knowledge (other alternatives appear in Harris, Moffit, and Kravtsova 2011;Kelejian and Piras 2017). This matrix is vital in the specification and must conform to the requisites described in Kelejian and Prucha (2004). The terms on the main diagonal are zero and the row and column sums of W are uniformly bounded in absolute value, in the same way as which must exist. The model of (2) can be expressed in matrix form: . . . , λ T m ). Ω = Σ ⊗ I N , with Σ = [σ ij ; i, j = 1, . . . , T m ], ⊗ denotes the Kronecker product and p = p t . Moreover, X * is an N T m × (p − T m) that omits the columns of the intercepts (these terms cannot be spatially lagged).
The spatial SUR model specified in (3) can be easily extended to multi-dimensional panels with G > 1 ( Wang and Kockelman 2007;López et al. 2014), which is omitted here for simplicity. This multi-dimensional panel model with N > 1, T m > 1, and G > 1 has been used in López et al. (2014) to compare model selection strategies in the context of spatial SUR models, in López et al. (2020) to compare ML and 3SLS (based on IV) methods to estimate some types of spatial SUR models, and it has been applied by López et al. (2017) to analyze spatial spillovers in the public expenditures of local governments at a municipal level.

Testing for spatial effects
The Lagrange multipliers presented below test whether the simple nonspatial model of (1) is misspecified due to omitted, unobserved spatial factors. All of them are obtained in an ML framework. The general expression of the multipliers, including a detailed description of the evaluation process, can be found in Mur et al. (2010). First, we present a global test for the presence of spatial lags both in the SUR equation and in the errors: The model of the null is the SUR-SIM of (1), whereas the alternative is the SUR-SARAR. The corresponding Lagrange multiplier is called LM SUR SARAR . The next multiplier (LM SUR SLM ) maintains the model of the null, which is the SUR-SIM, but that of the alternative is the SUR-SLM, so: Implicitly, it is assumed that the λ parameters are zero. Finally, the LM SUR SEM places a SUR-SEM model in the alternative against, once more, the SUR-SIM in the null (assuming that the ρ parameters are zero): H 0 : λ t = 0 (∀t) vs H A : ∃t with λ t ̸ = 0 Bera and Yoon (1993) show that the LM SUR SLM and LM SUR SEM multipliers are oversized in the case of misspecifying the alternative hypothesis (i.e., we want to test the SUR-SIM model in the null using the LM SUR SLM of (5) when, in fact, the data have been generated by a SUR-SEM). This lack of robustness makes the identification of the correct spatial specification more complicated. The solution, as suggested in Anselin, Bera, Florax, and Yoon (1996), is to robustify the raw multipliers using the so-called robust Lagrange multipliers, denoted as LM * SUR SLM and LM * SUR SEM , respectively. The size of the robust multipliers is more balanced at the cost of a slight loss of power.
The LM tests can be used as a guide to select the correct specification in the case of spatial SUR models (López et al. 2014) similar to those that have been used in a cross section framework (Mur and Angulo 2009), but problems can occur under certain conditions, as explained by Piras and Prucha (2014). Alternatively, the user could decide on the most suitable specification based on deep knowledge about the problem and not on the pre-test results (for example, guided by economic theory).

Estimation of spatial SUR models
The estimations of spatial SUR models are a direct extension of the methodology used to estimate single, spatial cross-section models. In the case of spatially lagged variables, the endogeneity of the spatial lag must be dealt with and, in the case of spatially correlated errors, the non-spherical nature of the covariance matrix must be accounted for. Several estimation approaches have been suggested in the literature. The first is based on the ML principle under the hypothesis of normality. The main problem with this method is that it is highly demanding in terms of computational effort when the size of the cross section, N , increases due to the Jacobian determinant in the concentrated likelihood function. The second approach is based on the methods of 3SLS that use IV and on GMM. These methodologies are alternatives to the often unrealistic assumption of normality, and they help to avoid computational problems with big datasets. A review of these methodologies in the framework of spatial panel data, including spatial SUR models, is found in Anselin, Le Gallo, and Jayet (2008), and we present a condensed version for the specific case of spatial SUR models in this section.

The maximum likelihood estimation
Inference based on the ML method has been the most commonly used method to obtain the estimation of the coefficients of a spatial SUR model (e.g., Wang and Kockelman 2007;Baltagi and Bresson 2011;Lauridsen et al. 2010;Le Gallo and Dall'erba 2006;López et al. 2017). In this case, assuming that the errors are normally distributed, the log-likelihood function of (3) can be written as: The vector of parameters is η ⊤ = β ⊤ ; θ ⊤ ; ρ 1 ; . . . ; ρ T m ; λ 1 ; . . . ; λ T m ; σ ij . In total, there are P = (2p − T m) + 2T m + T m(T m + 1)/2 unknown coefficients. Wang and Kockelman (2007) and López et al. (2014) use numerical optimization techniques to solve the maximum likelihood estimation. Note that the log of the Jacobian terms, ln|B t | and ln|A t |, constrains the spatial parameters, ρ t and λ t , to lie inside the so-called stability interval (if the weighting matrix has been row-standardized, a rough approximation is MIN being the greatest negative eigenvalue of W). In practice, estimation consists of applying a nonlinear optimization to the concentrated log-likelihood function. Under the assumption that the spatial SUR model is correctly specified, the ML estimators will be consistent, efficient, and asymptotically normally distributed (Davidson and MacKinnon 1993).
A main obstacle in the practical implementation of the ML estimation is the need to compute a Jacobian determinant for an N -dimensional matrix. The classic solution to this problem is to decompose the Jacobian in terms of the eigenvalues of the spatial weights matrix. Other methods are based on Cholesky or lower-upper (LU) decomposition, which exploits the sparsity of the spatial weights matrix W (see the references included in Bivand, Hauke, and Kossowski 2013a and Piras 2010 about this problem).

The 3SLS estimation method based on instrumental variables
Due to the highly intensive processing requirements of the ML method when N increases and to the often unrealistic assumption of normality, other methods have been proposed to estimate spatial SUR models. In the case of SUR-SLM and SUR-SDM, the 3SLS method based on instrumental variables (IV) can be used to estimate these models. In these models, the presence of the spatially lagged dependent variables introduces a type of endogeneity that calls for an instrumental variable approach. Anselin (1988a, page 147) and Anselin et al. (2008) proposed this methodology, although it was not implemented. Anselin (2016) later proposed this same estimation method and described the adequate instruments. López et al. (2020), Tupy et al. (2021), and Páez et al. (2021) applied this methodology to estimate SUR-SLM models. The list of proposed instruments for SUR-SLM and SUR-SDM, along with a full description of the method, is available in López et al. (2020) where an extensive Monte Carlo evaluates this method and compares the results with the ML estimation.
Note that this estimation procedure is almost linear, which simplifies computational effort. On the negative side, we should keep in mind that the estimation of the ρ parameters is not restricted by the Jacobian term, so we can receive estimates of ρ well outside the stability interval (although the risk seems acceptable, see López et al. 2020).

Testing for linear restrictions in spatial SUR models
In this section we present the Wald test approach to test linear restrictions on the parameters of an estimated spatial SUR model. Specifically, we present a Wald test to check linear restrictions on the β parameters or to test linear restrictions on the spatial parameters (λ or ρ). Several papers have used the Wald test in spatial SUR models. For example, this test is applied in Cotteleer and Van Kooten (2012) and in He et al. (2015) to evaluate significant differences in the β coefficients. Páez et al. (2021) apply the Wald test to evaluate the null hypothesis of identical levels of spatial dependence in the ρ t parameters. Anselin (2016) describes the Wald test for β coefficients in a spatial SUR framework. Assuming that the model is correctly specified (emphasizing the importance of the misspecification tests discussed in Section 2.1), the ML and 3SLS estimates are asymptotically normally distributed. So, for the ML case, we can write (Lee 2004;Elhorst 2014): where η is the vector of the parameters, as defined in (6), and V(η) is the corresponding covariance matrix, as described in López et al. (2014). A similar result can be stated for the 3SLS case, for which we need the orthogonality conditions defined in Section 2.2.2 and the covariance matrix, whose details can be found in Wooldridge (2010, Chapter 7).
Using the asymptotic distribution of (7), we can test for any kind of linear restriction on the parameter vectors of the model, such as exclusion restrictions or cross-equation linear restrictions. The procedure is carried out in the usual way. First, we need to write down the (r × P ) matrix R, with r = rank(R), which captures the linear restrictions: The Wald statistic for testing (8) is Note that the restrictions on (8) must involve only βs or spatial parameters. With slight variations, we obtain a similar result for the 3SLS method. Alternatively, the LR approach could be used but it is a bit more costly in terms of computing time than its Wald counterpart.

Direct, indirect, and total impacts
A capital difference between time series and spatial models lies in the interpretation of the β parameters. In a time series context, the β parameters have an unequivocal meaning: they measure the expected increase in the explained variable as a consequence of a unitary change in a given regressor. This impact is constant across time, ∂y t ∂x t = β k , and its significance can be checked by a simple t ratio.
In a spatial setting, there are feedback effects. This means that a change in a regressor in a given location will trigger a chain of reactions that will spread all over the space, according to matrix W. Consider the case of a SUR-SLM model with only one equation, G = 1, whose reduced form can be written as: A is a block-diagonal matrix whose non-zero matrices are A t = I N − ρ t W; t = 1, 2, . . . T m.
Assuming that all the ρ parameters pertain to the so-called stability interval, A −1 t admits a power series expansion in terms of the spatial coefficients and the weighting matrix, which leads us to expression (9). Clearly, if we change any regressor in any part of the space, its effects will impact the explained variable wherever it is located. This is the key starting point for the spatial multiplier analysis. Unfortunately, each type of spatial model produces its own multiplier effects, (LeSage and Pace 2009), so it is difficult to offer general results apart from basic notions and definitions.
Using the classic notation in LeSage and Pace (2009), we can write the following for the t-th period of a SUR-SLM model: (10) We assume that each equation has the same p regressors.
, and x thn refers to regressor h in period t in region n, whereas y tn is the observed value of the explained variable in period t and region n. Note that unless there are changes in the β parameters or in the weighting matrices W, the multiplier matrices S h t (W) remain constant across time. Expression (10) is crucial because it helps us to map the chain of feedback effects inherent in the specification. For example, the impact of a marginal change in regressor h (h = 1, . . . , p) in region n on the explained variable located in region s is the (n, s) element of this matrix, Applied literature distinguishes three types of impacts: The so-called direct impact, where all the changes occur in the same region and coincide with the mean of the terms in the main diagonal of the multiplier matrix of (10), indirect impact which measures spill-over effects, evaluated as the mean of the off-diagonal terms of S h t (W), and total impact, as the sum of the two effects. In an obvious notation: with τ N as an N × 1 unitary vector and tr() as the trace operator. The subindex g denotes that these effects correspond to the g-th equation of the SUR (if G = 1, the subindex g should be dropped).
Finally, for the case of G > 1, we can aggregate the multipliers in (11) to obtain overall measures for each variable in the SUR: The multipliers of (12) should be accompanied by measures of statistical significance to be useful. To simplify, let us assume that the same list of regressors appears in the G equations of the SUR. LeSage and Pace (2009) suggest two possibilities. The first is the Bayesian MCMC approach, which draws samples of parameters from their posterior distribution and computes random sequences of multipliers to calibrate the estimates. To our knowledge, there are no results available on applying MCMC in a SUR framework. The second alternative, implemented in spsur, is simulation using the maximum likelihood or 3SLS estimates, which involves the following procedure: • Estimate the spatial SUR model either by ML (preferably) or 3SLS.
• Draw B samples of the random vectors of the parameters for our spatial SUR, η, using the estimated multivariate normal distribution, N η;V(η) , whereη andV(η) are their respective estimates. η b denotes the b-th draw.
• Re-evaluate the effects using the simulated vectors η b and maintaining the data of the regressors. In this way, a sequence of values: • Finally, calculate the corresponding dispersion measure as the standard deviation of the series of estimated multipliers: with Effect = Total, Direct, Indirect.
spdep and spatialreg are two powerful packages with high functionalities for the spatial econometric analysis of cross-sectional data, whereas spml is focused on the area of spatial panels. On the other hand, the R package systemfit deals with SUR models but does not incorporate spatial effects. To our knowledge, only the embryonic package, spse (Piras 2018), downloadable from the R-Forge repository in https://rdrr.io/rforge/spse, offers some functionalities to estimate spatial SUR models (see Bivand, Millo, and Piras 2021 for a review of software for spatial econometrics in R). A comparative vignette of the common functionalities of spsur and spse is available in the GitHub repository at https://github.com/rominsal/spsur for the development version of the spsur package. In summary, none of these packages are well-equipped to deal with the overall inference problems of spatial SUR models.
The most similar antecedent to spsur is the beta version of SpaceStat (Anselin 1992) called SpaceStat (version May91), with limited functionalities. Moreover, it requires Windows 98 (see the YouTube video: https://www.youtube.com/watch?v=6pM0tDWqt3o&t=53s). Of course, nowadays it would be very difficult to do applied research using SpaceStat (version May91). There is a second alternative in Python, in the PySAL library (Rey and Anselin 2009), which is the module spreg.sur. It provides SUR estimations of models with no spatial impacts and diagnostic measures, among which we find the raw Lagrange multipliers. Moreover, several functions provide the estimation of spatial SUR models with spatial errors or spatial lags of the explained variable, respectively. The SUR-SEM model is estimated by ML, while the SUR-SLM is estimated using spatial instrumental variables (3SLS based on IV). Comparison of the results of an estimated spatial SUR model with PySAL and spsur and spse is available in a vignette of spsur. Finally, there is some code in MATLAB connected to the work by López et al. (2014), but the source is not free software.

Introducing the spsur package
spsur (López, Mínguez, and Mur 2022) is an open-source software for the R computing platform (R Core Team 2022). The package can be freely downloaded from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=spsur. Users will find complete assistance documentation and some vignettes to guide them through the functionalities of the package, as well as examples of its use.
spsur is made of 10 different functions that allow for a thorough analysis of a spatial SUR specification. Figure 1 presents a graphic description of the package. spsur includes one function to evaluate the Lagrange multipliers for the presence of spatial correlation in a SUR-SIM model, lmtestspsur(). For the general case of G > 1, spsurml() and spsur3sls() allow the estimation of the model by ML and 3SLS (based on IV). The function spsurtime() is directed specifically toward spatial SUR models, where G = 1 and T m > 1. Four functions are designed to improve specification: anova(), wald_betas(), wald_deltas(), and lr_betas(); while impactspsur() computes the direct, indirect, and total spatial effects. Finally, dgp_spsur() allows the user to solve Monte Carlo experiments using spatial SUR models.
The functions to estimate spatial SUR models spsurml() and spsur3sls() return an object of the class 'spsur'. This class has been specifically defined in the package and several methods have been implemented to work with this object. The method summary() presents a full output of the estimated model, print() presents a short output of the results of the estimation, and plot() generates a plot with the coefficients. To conclude, it is important to point out that spsur includes several vignettes to guide the user. One vignette shows the main functionalities of the package and another presents some results from the use of spsur for single equation and baseline models. Two additional vignettes, available in the GitHub version, present a sort of Monte Carlo exercise and a comparison with spse and PySAL.

Datasets and baseline model
Three datasets have been included in the package; namely, spc, NCOVR and spain.covid. The first two are well known in the spatial econometric literature and are used in this paper to show the performance spsur. The third is included in Section 5 to illustrate an empirical case where the spatial SUR approach is useful in econometric methodology. Each dataset consists of an 'sf' object (Pebesma 2018), allowing the construction of a list of neighbors and plot maps.
The spc dataset: The first dataset, spc, constitutes a classic example of a Phillips curve, taken from Anselin (1988a, pages 203-211). This dataset includes an 'sf' object of 25 counties in Southwest Ohio with information about different variables for two time periods, 1981 and 1983. The explained variable measures the changes in wage rates (WAGE), while the regressors are unemployment rate (UN), net-migration rate (NMR), and a dummy variable (SMSA) with a value of 1 for metropolitan counties.
The SUR-SIM model estimated by Anselin (1988a) is: We use the Formula package (Zeileis and Croissant 2010) to specify the multiequational SUR model of (13): Note that we first select the two dependent variables appearing in each equation, separated by a vertical bar (|). Next, after the prime, we specify the regressors for each equation, once more using a vertical bar to separate the two groups.
The NCOVR dataset: The second dataset (NCOVR) contains data downloaded from the GeoDa Data and Lab collection about homicide rates in 3,085 continental US counties for four years (1960, 1970, 1980, and 1990). The dataset includes a large number of socio-economic characteristics for these counties. This dataset comes from the paper by Baller, Anselin, Messner, Deane, and Hawkins (2001) R> data("NCOVR", package = "spsur") Following Baller et al. (2001), we choose a W matrix based on the k-nearest-neighbors, with k = 10, using the geographic centroids of US counties.
R> library("spdep") R> library("sf") R> co <-sf::st_coordinates(sf::st_centroid(NCOVR.sf)) R> ncovrlw <-nb2listw(knn2nb(knearneigh(co, k = 10, longlat = TRUE)), + style = "W") We consider the next SUR-SIM model with three equations and different numbers of regressors in each one only to show the functionalities of spsur: HR 80 is the homicide rate per 100,000 inhabitants in 1980, PS 80 measures the population structure in 1980, UE 80 is the unemployment rate, DV 80 is the divorce rate, FP 79 is the percentage of families below the poverty line in 1980, and SOUTH is a dummy variable for Southern counties. Note that the model of (14) has only one cross section.
To input this model into R, we write, old-style crs object detected; please recreate object with a recent sf::st_crs() old-style crs object detected; please recreate object with a recent sf::st_crs() old-style crs object detected; please recreate object with a recent sf::st_crs()  Figure 2 shows the spatial distribution of homicide rates, in 1980, where a clear spatial structure emerges. Therefore, we expect to find spatial autocorrelation at least in the first equation.

Misspecification tests for omitted spatial effects
The generic function lmtestspur() computes five Lagrange multipliers for omitted spatial dependence in a SUR-SIM model. The main advantage of the LM tests is that they only need information of the null hypothesis, and it is therefore easy to obtain. The full syntax of the first one is: lmtestspsur(formula, data, listw, na.action, time = NULL, Tm = 1, zero.policy = NULL, R = NULL, b = NULL, ...) The lmtestspur() function includes a 'formula' object, a data frame in data argument, a list of neighbors in a 'listw' object (it could also be a W neighbor matrix), zero.policy, and na.action for the treatment of missing values (see Bivand and Wong 2018, for details The output of the lmtestspsur() function includes a list of 'htest' class objects for each of the LM tests. The situation depicted by this battery of multipliers shows that the SUR-SIM model of (14) is clearly misspecified, but the set of LMs do not reveal the type of misspecifications that afflict the model.

Maximum Likelihood estimation of spatial SUR models
A spatial SUR model must be estimated in case the LM tests reject the SUR-SIM alternative. The function spsurml() obtains the ML estimation of the different spatial SUR models listed in Section 2. The general syntax for this function is: spsurml(formula = NULL, data = NULL, na.action, listw = NULL, type = "sim", Durbin = NULL, method = "eigen", zero.policy = NULL, interval = NULL, trs = NULL, R = NULL, b = NULL, X = NULL, Y = NULL, G = NULL, N = NULL, Tm = NULL, p = NULL, control = list() )

Example with spc
The symptoms of spatial misspecification in the Phillips curve are rather weak. The user can obtain all the LM tests with the function lmtestspsur() shown in Section 3.2, which we omit to save space. However, as an example, suppose we conclude that the SUR-SIM model is misspecified because we have omitted a spatial lag of the explained variables on the right side of the equation (that is, we interpret that the correct specification is a SUR-SLM model, which is plausible if we fix the significance level at 10%).
To estimate a SUR-SLM model for the spc dataset, the syntax is: R> spcsur.slm <-spsurml(formula = spcformula, data = spc, type = "slm", + method = "eigen", listw = Wspc) The listw argument allows the specification of neighbors either as a list or as a matrix. The method argument (the same as the argument of the estimation functions in the spatialreg package) permits the choice of method to compute the determinant of the spatial Jacobian in the log-likelihood. For relatively small samples (hundreds of observations), the simplest method is "eigen" (which is also the default). Nevertheless, the computation of the eigenvalues becomes very expensive for medium to large samples. The function spsurml() offers other options to obtain the spatial determinant without computing the eigenvalues; for example, the Cholesky or LU factorization ("Matrix" or "LU"), or the Chebyshev approximation ("Chebyshev"). See Bivand et al. (2013a) and Chapter 4 in LeSage and Pace (2009) for details. The zero.policy, interval, and trs arguments are also similar to those available in the spatialreg package (see ?spatialreg::lagsarlm for the details of each one).
Finally, the options included in the control argument are almost identical to those available in the lagsarlm() function of spatialreg (see ?spsurml for details of all of them). Some of these options include the possibility of checking for intermediate results during the optimization process using trace = TRUE, setting the threshold for the convergence of the log-likelihood in tol, the maximum number of iterations in maxit, and the possibility of computing the numerical covariances using fdHess = TRUE (the function computes analytical covariances by default). Note that the computation of analytical covariances can be quite slow for large samples (thousands of observations).
The output of the function spsurml() is an 'spsur' object and, consistent with the conventions in the R environment, the summary() method prints the output equation by equation. The summary includes standard information for a SUR estimation, disaggregated by equation, such as estimated coefficients, standard deviations, etc. The matrices of covariances and of correlations between the SUR residuals appear next. Finally, spsurml() shows some measures of goodness of fit (overall R 2 , log-likelihood), the usual Breusch-Pagan (BP) test of diagonality (Breusch and Pagan 1980), and a marginal Lagrange multiplier (LMM) to test for the omitted spatial effects arising from the model that we have estimated (that is, if the type is slm or sdm, the LMM tests for omitted spatial errors; if the type is sem or sdem, the LMM tests for omitted spatial lags).
Additional standard methods for spsur objects include anova, coef, fitted, logLik, print, plot, residuals, and vcov. All of them work in the usual way for R standard methods.

A comparison of spatial SUR specifications using likelihood ratios
As is well known (Elhorst 2014), the spatial models listed in Section 2 are related in a nesting sequence. Given that spsurml() operates in an ML framework, we can easily compare two nested models using the corresponding likelihood ratios (LRs) to select the best specification. This is the purpose of the anova() method applied to 'spsur' objects. In the case of a Phillips curve, we simply write: R> spcsur.slm <-spsurml(formula = spcformula, data = spc, type = "slm", + listw = Wspc, control = list(trace = FALSE)) R> spcsur.sdm <-spsurml(formula = spcformula, data = spc, type = "sdm", + listw = Wspc, control = list(trace = FALSE)) R> anova(spcsur.slm, spcsur.sdm) The output of the function is an 'anova' object showing the value of the corresponding LR test (if the lrtest argument is TRUE as is the default) and its associated p values, obtained from a χ 2 distribution with degrees of freedom equal to the number of restrictions in the respective null hypothesis. Akaike information criterion (AIC) and log-likelihood (logLik) values for each model are also included. The information provided by the likelihood ratios may help the user define the spatial structure needed for the case under consideration.

3SLS based on IV estimation of SUR-SLM and SUR-SDM models
The function spsur3sls() estimates spatial SUR models by spatial IVs in a 3SLS framework, as described in Section 2.2.2. The syntax of the function is similar to spsurml(). spsur3sls(formula = NULL, data = NULL, na.action, R = NULL, b = NULL, listw = NULL, zero.policy = NULL, X = NULL, Y = NULL, G = NULL, N = NULL, Tm = NULL, p = NULL, type = "slm", Durbin = NULL, maxlagW = NULL, trace = TRUE) The main differences between the two are: (i) the argument type now only supports two spatial SUR models, slm and sdm, and (ii) a new argument is needed, maxlagW, to set the maximum order of the spatial lags of the regressors and obtain the matrix of IVs. The default value is NULL, which chooses maxlagW = 2 for SUR-SLM models and maxlagW = 3 for the SUR-SDM variant.

Example with NCOVR
The 3SLS method is especially adequate in cases of large datasets or when the Jacobian matrix is very dense. The NCOVR dataset including 3,085 US counties is a good candidate for spsur3sls(), for which we simply need: R> ncovrsur.slm.3sls <-spsur3sls(formula = ncovrformula, data = NCOVR.sf, + type = "slm", listw = ncovrlw, trace = FALSE) As an alternative to summary, the user can apply the print method to get an abbreviated version of the results, If a variable is not included in an equation, the print method shows an NA value. As in the case of spsurml(), the output of the function spsur3sls() is an 'spsur' class object.

Testing the spatial SUR estimation and re-estimating models
An important feature of spsur is its ability to thoroughly test the specification. spsurml() routinely reports the BP diagonality test as well as a Lagrange marginal multiplier for omitted spatial effects. To these basic measures, we add three functions, wald_betas(), wald_deltas(), and lr_betas(), which improve specification.

Testing for linear restrictions on coefficients
The function wald_betas() is used to test linear restrictions in the β parameters. This function needs three arguments. The first is an 'spsur' class object obtained with the functions spsurml() or spsur3sls(), the second is a matrix that must reflect the restrictions on the β parameters that we want to test, and the third is a vector with the values of the restrictions under the null hypothesis.
For example, to test whether the intercepts of the two equations of the SUR-SLM model of the Phillips curve (13) are equal: It is necessary to write the code, R> R1 <-matrix(c(1, 0, 0, 0, -1, 0, 0, 0), nrow = 1) R> b1 <-matrix(0, ncol = 1) R> wald_betas(spcsur.slm, R = R1, b = b1) Wald test on beta parameters data: spc Wald test = 0.39767, df = 1, p-value = 0.5283 The output of the test is an 'htest' object. In this case, the two intercepts are not statistically different from each other, and the correct specification must therefore include the same intercept in both equations. The constrained version of the spatial Phillips curve is: Both functions, spsurml() and spsur3sls(), allow for a restricted ML or 3SLS estimation of the SUR model, using the arguments R = R1 and b = b1, which expedites the estimation problem. The output of the constrained ML estimation appears below.

Testing for linear restrictions on the spatial parameters
Similarly, the purpose of the function wald_deltas() is to obtain a Wald test to evaluate linear restrictions on the spatial parameters of the model (λ or ρ parameters). The use of the function is totally analogous to wald_betas() and extends to the ML and 3SLS methods. For example, the user may be interested in testing whether the λ parameters of the two equations of the constrained version of the previous Phillips curve are equal; that is: R> R2 <-matrix(c(1,-1), nrow = 1) R> b2 <-matrix(0, ncol = 1) R> wald_deltas(spcsur.slm.restricted, R = R2, b = b2) Wald test on spatial delta parameters data: spc Wald test = 16.932, df = 1, p-value = 3.874e-05 In this case, the Wald test clearly rejects the hypothesis of (16), confirming that the specification (15) with two different parameters of spatial dependence is correct.

Spatial impacts in the SUR framework
spsur includes the function impactspsur() to evaluate spatial effects. The function is a wrapper in a multiequational context of the impacts() function available in 3 spatialreg (see Bivand and Piras 2015, for details). As a consequence, the arguments of impactspsur() are a subset of the corresponding arguments in impacts(). Specifically, the full syntax of impactspsur() is: impactspsur(obj, ..., tr = NULL,R = NULL, listw = NULL, evalues = NULL, tol = 1e-06, empirical = FALSE, Q = NULL) The argument obj includes an 'spsur' object (the output of spsurml(), spsur3sls() or spsurtime()). The rest of the arguments can be checked using either the help of impacts() or the help of impactspsur().

Example with spc
Below, an example of how to compute the impacts using the previous Durbin model estimated for NCOVR.sf data is shown. We also compute the vector of traces of powers of the spatial weights matrix (using the trW() function of spatialreg) to speed up the computations.

Simulated spatial SUR datasets
This section describes an additional function included in spsur that offers the possibility of generating random datasets with the required dimensions and spatial structure. This function, dgp_spsur(), has great flexibility and may be of particular interest in two cases: (i) as part of a larger simulation where the users need random datasets with specific features (i.e., for their own research projects) or (ii) with the aim of showing specific properties of spatial SUR models and that inferential procedures related to them (i.e., for teaching activities).
dgp_spsur() is a data generating process (DGP) where the user decides the dimensions of the required dataset, specifies the values of the parameters that intervene in the corresponding SUR model, and selects the distribution functions from which the regressors and the random terms are to be drawn. The syntax of the function is the following: dgp_spsur(Sigma, Tm = 1, G, N, Betas, Thetas = NULL, rho = NULL, lambda = NULL, p = NULL, listw = NULL, X = NULL, type = "matrix", pdfU = "nvrnorm", pdfX = "nvrnorm") The dimensions of the dataset are defined by the arguments T m, N , and G. The argument Sigma specifies the covariance matrix among the residuals of the G equations, which is the core of the SUR model. Of course, a (N × N ) spatial weighting matrix should be uploaded in the argument listw.
A fundamental piece of information is the argument p, which defines the number of regressors (X variables) that appear in each equation. If p is a scalar, every equation has the same number of regressors, whereas a 1×G vector indicates that the G equations contain a different number of regressors. Subsequently, the user should specify the parameters that intervene in the equations of the SUR. This is the purpose of the arguments Betas, Thetas, rho, and lambda, which are defined as row vectors of the adequate dimensions. Betas is a row vector of the order 1 × K, where K is equal to pG, or to K = Σ G j=1 p j if p is a row vector. The values assigned to these parameters determine the type of spatial model required. For example, if Thetas = NULL, rho = NULL but lambda and Betas are not NULL, this means that the user is specifying a SUR-SEM model, or a SUR-SDEM model in the case of rho = NULL but Thetas, Betas, and lambda are different from NULL.
There are two possibilities to build the matrix X of regressors. In Monte Carlo experiments, it is usual to maintain the values of the regressors fixed and draw random matrices only for the error terms. If this is the case, the user should upload the required X matrix in the argument X, which must be consistent with the dimensions of the SUR model. If X = NULL, the user randomly draws this matrix using the function dgp_spsur(). This is the purpose of the argument pdfX, which uses a multivariate standard normal distribution. The alternative is a Uniform [0, 1] distribution for each regressor. In both cases, the regressors are generated independently. Finally, the argument dfU selects the multivariate probability distribution function to draw the values of the error terms.
Output types can be selected using the type argument. By defect, the output of dgp_spsur() is a vector called Y of the order (N T mG × 1), with the values generated for the explained variable in the G equations of the spatial SUR model. If the argument X is set to NULL, the user will receive another matrix called X with the values generated for the regressors of the SUR.

A toy example. How to create random datasets
In this section, we present a toy example of a spatial SUR process with two equations for 537 individuals defined on a regular hexagonal lattice. The geometry of the lattice is defined on a 1 × 1 square using the sf package.
R> W <-nb2listw(poly2nb(as(hexs.sf, "Spatial"), queen = FALSE)) For the DGP example, we consider two equations and 2 regressors, plus the intercepts. We generate a SUR-SLM model with different levels of spatial dependence in each equation, as shown in (17)).
R> sur.dgp.sim <-spsurml(Y = dgp.spatial.sur$Y, X = dgp.spatial.sur$X, + type = "sim", G = 2, Tm = 1, N = 537, p = 3, listw = W, + control = list(trace = FALSE)) R> res <-residuals(sur.dgp.sim) As expected, the results are quite satisfactory. All the estimated parameters are within the interval (theoretical value plus/minus two standard deviations). The Breusch-Pagan test (with a value of 148.57) clearly rejects the assumption of diagonality. The correlation obtained for the SUR equations are slightly inflated (0.516 on average for an expected value of 0.5), whereas the marginal multiplier (6.47) clearly confirms the adequacy of the spatial lag model.

spsur application to mobility due to COVID-19
There is broad consensus that this pandemic has been one of the most relevant events that has occurred during recent decades with impacts (social, medical, economic, etc.) that will last for a long time. However, in spite of the dramatic consequences, social distancing and the reduction of individual mobility are seen as the most effective strategies to slow down the spread of the disease.
This section illustrates the practical use of the spsur package to develop a spatial SUR model and explain the effect of COVID-19 on two types of individual mobility (controlling for several factors). In particular, we focus on the case of Spain, using data at a provincial level (NUTS3 in Eurostat terminology) during the first weeks of the pandemic 4 . The hypotheses that we explore is that the high incidence of COVID-19 dissuades individuals from leaving their homes and therefore reduces mobility.
For this illustration, we consider two different indicators of individual mobility: "Within" and "Exits". The first is defined by the ratio between the number of intra-provincial trips within province "i" in week "t" compared to a pre-COVID week (February, 14-21, 2020). The second is defined in the same way but for inter-provincial trips (see Table 1). Figure 4 shows the weekly evolution of these mobility indices; Figures 4A for intra-provincial and 4B for inter-provincial trips. Note the sharp reduction of the mobility indices, both intra and inter, occurring in week 4. This is not accidental, since week 4 is the beginning of the state of emergency (StEm) declared by the Spanish Government. Weeks 6 and 7 (in orange), correspond to the tightening of restrictive measures when no essential economic activities were allowed. The StEm remained in effect until June 21, the end of week 17.
Several variables could have impact on both mobility indices (Table 1). We consider three categories of variables: government restriction orders, socio-economic characteristics, and week . In orange, total lockdown, weeks 6 and 7 (from March 28 to April 10). The red line marks the mobility level = 1 (the reference week).
COVID-19 incidence. With respect to government restriction orders, the Spanish Government declared a partial lockdown with stay-at-home orders for citizens for 10 weeks (weeks 4-13), but essential economic activities were permitted. In the middle of this period (two weeks), the lockdown was total, and essential economic activities were not permitted. Therefore, it is expected that during this period, mobility would decrease. Therefore, the expected sign of these variables explaining mobility is negative.
Even though mobility restrictions were dictated for the entire country, the Government's orders were unevenly followed in different provinces. Several factors related with the socioeconomic characteristics can explain these regional differences.
In the first place, differences in population density have an impact on mobility (Brodeur, Grigoryeva, and Kattan 2021). In this case, high levels of population density (the provinces of Madrid or Barcelona) can change individual behavior due to the fact that the virus spreads mainly through interpersonal interactions (expected negative sign). In the second place, something similar occurs in provinces with a high proportion of older people, whose ability to move is more restricted (Brodeur et al. 2021).They are also more responsive to restriction orders (Páez et al. 2021, Engle, Stromme, andZhou (2020)). In the third place, provinces with high levels of essential services reduce mobility less than provinces with low levels of essential services (expected negative sign, Abdullah, Dias, Muley, and Shahin 2020). Lastly, the incidence of COVID-19 could influence movement in the sense that evidence of the disease dissuades individuals from leaving their homes, and therefore mobility is reduced (Engle et al. 2020, Abdullah et al. (2020). Table 1 summarizes all the variables considered in this illustration and their expected signs. Figure 5A represents the weekly incidence of the pandemic (log of cases detected per 100,000 inhabitants). The incidence has an inverted-V shape that peaks on week 6. Looking at this figure, it is apparent that the incidence of COVID-19 has a strong impact on mobility, since movement decreases considerably as the virus intensifies. Figure 5B shows the spatial distribution of COVID-19. A strong spatial pattern emerges, confirming the spatial diffusion of the pandemic and therefore, spatial effects in the model are expected. In gray, partial lockdown weeks 4-5 and 8-13 (from March 14-27 and April 11-23). In orange, total lockdown, weeks 6 and 7 (from March 28 to April 10). The blue line marks zero incidence. (B) Quartiles of provincial COVID-19 incidence for week 6 (in log). Similar spatial pattern in other weeks.

Variable Sort Description Sign Within
Number of trips in week "t" within province "i" (compared to the pre-COVID reference week) Exit Number of trips in week "t" with origin in province "i" and arrival to another (compared to the pre-COVID reference week) Emergence Dummy variable. 1 if StEm is active in week "t". Essential activities (e.g. food, health) are allowed. 0 otherwise Percentage of population aged 65 and older in province "i" (−) Essential Percentage of firms in province "i" with essential activities (+) Incidence Weekly incidence in the week "t − 1" in the province "i" (in logs)

The SUR approach
The dataset with the two mobility indices and the control variables for the 50 spatial units (provinces) and 17 temporal periods (weeks from February 14 to June 21, 2020) is available in the spsur package, R> data("spain.covid", package = "spsur") A baseline SUR-SIM model with G = 2 equations (for inter and intra movements, respec-tively), which does not include spatial effects can be specified as follows: . . . , 17; i = 1, . . . , 50 (18) To estimate this model using the spsur package, first, the two equations of model (18) can be specified using the Formula package. In this case, only a single sequence of explanatory variables is included on the right-hand side of the formula because these factors are the same for both equations.

R> formula <-Within | Exits~Emergence + EmergenceTotal + Density + + Old65 + Essential + Incidence
This specification can be estimated for maximum likelihood with the function spsurml() R> Tm <-17 R> covid.sim <-spsurml(formula = formula, data = spain.covid, type = "sim", + Tm = Tm, control = list(trace = FALSE)) R> summary(covid.sim) Call: spsurml(formula = formula, data = spain.covid, type = "sim", Tm = Tm, control = list(trace = FALSE)) Spatial SUR model type: sim  The residuals of both equations are highly correlated (0.69), and the Breush-Pagan test (407.1) confirms the adequacy of the SUR specification. All the variables, with the exception of Density in the second equation, are significant. The signs of the two dummy variables, (Emergence and EmergenceTotal) are negative, as expected. Population density (Density) is significant and negative in the first equation and not significant in the second one. The variable for older people, Old65, contrary to expectations, has a positive sign, which is difficult to justify in this first model. The percentage of essential services (Essential) has a positive impact on movement in the sense that greater economic activity promotes mobility. Finally, a strong incidence of the disease in previous weeks, Incidence, has a negative sign, showing that a high incidence of COVID-19 discourages individuals from moving.
Given the framework that supports our research and using Spanish provinces as observation units, it is advisable to test for the presence of (omitted) spatial effects, especially spatial autocorrelation. We consider a weighting matrix based on the queen criteria to test the hypothesis of spatial independence in the residuals. The spdep package can be used at this moment. Note that there are three provinces with no neighbors (in the Balearic and Canary islands), so the option zero.policy = TRUE will be necessary, as appears in the code below.
R> listw <-spdep::nb2listw(listw, style = "W", zero.policy = TRUE) The spsur function lmtestspsur() helps test the SUR-SIM model against any other specification through the corresponding Lagrange multipliers. The results, as shown below, detect a strong spatial structure that has been omitted in the SUR-SIM: all the LM tests, robust and non-robust, reject their corresponding null hypotheses.

Spatial SUR models. Selection strategy.
The next question is obvious: if the SUR-SIM model is not appropriate, which is the best specification, from a spatial perspective, to model the relation between mobility and COVID-19? Since the LM tests don't give an a priori specification preference, we can estimate all the possible spatial SUR models listed in Figure 6 and obtain several statistical measures -LogLik, AIC, BIC (Bayesian information criterion) -to select the most suitable model.
The spsur user can estimate all the models by changing the type argument. For SUR-SLX, SUR-SDM, SUR-SDEM, and SUR-GNM specifications, we consider only spatial lags for a subset of the independent variables. The Durbin option can be used for this purpose.
R> formulaD <-~Essential + Incidence R> covid.slx <-spsurml(formula = formula, data = spain.covid, type = "slx", + Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE) R> covid.slm <-spsurml(formula = formula, data = spain.covid, type = "slm", + Tm = Tm, listw = listw, zero.policy = TRUE) R> covid.sem <-spsurml(formula = formula, data = spain.covid, type = "sem", + Tm = Tm, listw = listw, zero.policy = TRUE) R> covid.sdm <-spsurml(formula = formula, data = spain.covid, type = "sdm", + Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE) R> covid.sdem <-spsurml(formula = formula, data = spain.covid, type = "sdem", + Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE) R> covid.sarar <-spsurml(formula = formula, data = spain.covid, + type = "sarar", Tm = Tm, listw = listw, zero.policy = TRUE) R> covid.gnm <-spsurml(formula = formula, data = spain.covid, type = "gnm", + Tm = Tm, Durbin = formulaD, listw = listw, zero.policy = TRUE) At this point, the anova method can be applied to objects of the 'spsur' class to get information about the LogLik of each model and the information criteria AIC and BIC. Note that the lrtest = FALSE must be selected because the eight models are not nested. The results show a clear decrease in both AIC and BIC for the most complex models. Moreover, models SUR-SDEM and SUR-GNM show the best performance in the information criteria. In cross-section, the criteria of AIC and BIC have been used to select the correct specification (Agiakloglou and Tsimpanos 2021). Given that both models are nested (see Figure 6), we can use the LR test to choose between them. Again, the anova() method can be used, but in this case with lrtest = TRUE to show the LR test of both nested models. At this point, it is necessary to note that although in terms of likelihood, the SUR-GNM is a higher model than SUR-SDEM, the SUR-GNM could lead to an overestimated model and to problems identifying spatial dependence parameters (e.g. Burridge, Elhorst, and Zigova 2016;Vega and Elhorst 2015). Therefore, the alternative SUR-SDEM is finally selected. For example, for spatial cross-section models Elhorst (2014, page 33) points out his preference for SDEM versus GNM. In contrast, Vega and Elhorst (2015, page 341) claim that "only the parameters of a spatial econometric model with all possible spatial interaction effects based on arbitrary spatial weights matrices have not been proved to be free of this type of identification problem". Finally, we decide to follow the specific-to-general strategy and consider the SUR-GNM as the final model. The parameters of spatial error dependence λ 1 , λ 2 take values close to 1 (the maximum value) and are highly significant as a symptom of the strong spatial structure in the residuals of this model, probably due to the omission of relevant variables.
On the other hand, the parameters of substantive spatial dependence ρ 1 , ρ 2 take values close to zero but are significant, with p values lower than 0.05. In the first equation, the ρ 1 parameter has an opposite sign to λ 1 . This might be a symptom of the erroneous identification of the spatial dependence parameters. However, the signs of ρ 2 and λ 2 are both positive in the second equation. Again, the possibility of selecting the SUR-SDEM specification as an alternative to SUR-GNM is on the table.

Testing linear restrictions on the parameters of the model
Several tests can be used to reduce the model specification and simplify the model. In the first place, it is possible to test whether government restrictions have had the same impact on reducing mobility when people stay inside and when they leave the province. The function wald_betas() can be used to obtain the Wald test of coefficient equality. Following the notation of (18), the null hypothesis is, The code to test this hypothesis is (see ?spspur::wald_betas) In this model, the variable Old65 in the first equation, and lag.Incidence in the second one can be removed because they are not significant. To obtain a new model excluding these variables, we can rewrite the formula with different variables on the right-hand side of the formula. The R vector must also be adapted to the new dimension of the model,

Impacts of increased COVID-19 incidence on mobility
To finish the application, it is necessary to interpret the way that changes in the explanatory variables impact both mobility types (inter-and intra-provincial). These impacts can be obtained from the restricted SUR-GNM model. The spsur package implements the function impactspsur() to compute direct, indirect, and total impacts. The output of this function is a list where each element shows the impacts of one equation, R> impacts.covid.gnm.rr <-impactspsur(covid.gnm.rr, listw = listw, R = 1000) The impacts of inter-provincial mobility are (standard errors, z statistics, and p values, are not included here to save space, they can be obtained by changing zstats = TRUE), We will focus the discussion only on the impact of incidence rates on mobility. Similar analysis could be made for the rest of the variables. In the case of the first equation, an increase of one percent in incidence has a direct impact on intra-provincial mobility, with a reduction of 1.55%. This increase in incidence also has an indirect impact (overflow effect), and this effect represents an average reduction of 1.65%. Both direct and indirect impacts involve a total impact of 3.20% on intra-provincial mobility. The summary of the impacts of the second equation for inter-provincial mobility can be obtained with, In this case, the impact of a one-percent increase in the COVID-19 incidence rate on intraprovincial mobility is lower than the previous case for inter-provincial mobility. The direct impact has a reduction on average of 0.98% and the indirect impact, an average of 0.11%. Therefore, the total impact of an increase of one unit in the incidence rate on inter-provincial mobility is 1.09%.

Conclusions and directions for future research
This article introduces a new R package, spsur, to estimate and make inferences on spatial SUR models. Our intention is to supply researchers interested in spatial econometric methodologies with free, user-friendly software to estimate these types of specifications. spsur completes the ecosystem of the R packages, spdep, sphet, spse, spatialreg, and splm, oriented to deal with spatial econometric specifications.
The package consists of 10 inter-related functions. lmtestspsur() tests for the existence of spatial effects in a specified model without spatial mechanisms; the so-called SUR-SIM model. Assuming that some of these tests are statistically significant, three functions, spsurml(), spsur3sls(), and spsurtime() allow the estimation, by ML or 3SLS based on IV, of the preferred spatial specification. The usual methods of extracting information from estimated models, such as coef(), fitted(), logLik(), plot(), print(), residuals(), and vcov(), are also included. Four additional functions, wald_betas(), wald_deltas(), lr_betas(), and anova() can be used to improve the specification of the model if necessary. The function impactspsur() is useful to obtain the spatial multipliers of the variables in the estimated model, and it comes with a collection of measures of statistical significance. Finally, the function dgp_sur() is intended to illustrate particular features of spatial SUR specifications as an aid for teaching activities or to implement the results of the Monte Carlo in a more general research project.
Of course, this version of the package should be improved in the future. Included in a list of actions already on the coauthors' table are: (i) the extension of 3SLS to treat endogenous regressors other than the spatial lag of the explained variable and the use of external instruments; (ii) new tests to improve specification; e.g., testing for the instability of the spatial dependence coefficients between equations; (iii) more efficient treatment of the unobserved effects, allowing their estimation, expanding the options to demean the data, (iv) the introduction of time dynamics in a strictly stationary framework; and (v) the consideration of more general SUR models with different types of spatial effects in each equation. Some of this work is already underway and a development version is available from the GitHub repository, including a new function named spsurgs3sls() to estimate several spatial SUR models using the GS3SLS. This function is a wrapper of spreg() from the sphet package, and it implements the GS3SLS estimation method in a SUR framework following the steps explained in Kelejian and Piras (2017, pages 304-305). This function also allows the inclusion of additional endogenous regressors other than the spatial lag of the explained variable.