1 Introduction

Seemingly unrelated regression models (SUR from now on) are a type of multivariate econometric formulation very popular since the seminal papers of Zellner (1962), Malinvaud (1970) or Theil (1971). SUR have been used in many research areas in economics and others fields (for instance, see Fiebig 2007), especially in the topics of consumption, production and environmental economics, but also on economic growth, health economics, real estate or competition between local agents. Distinct mechanisms may result in the type of dependence that lies behind SUR models such as omitted variables, unobserved effects or common factors, with a dfferent impact for each equation but uniform among the individuals. It is clear that if the errors are connected, a SUR framework assures efficiency gains by estimating the system jointly rather than processing each equation separately. These gains explain the relevant position that SUR models occupy in the applied research agenda. Spatial econometrics is, surprisingly, an exception: the applied literature is very short in spite of a quite extensive methodological research.

Anselin (1988) introduced the term spatial SUR in reference to a model made of ‘an equation for each time period, which is estimated for a cross section of spatial units’ (p. 141). The approach of Anselin focuses on the problem of serial dependence in the errors of an equation, which is perfectly natural if the same agents intervene in every period. According to usual practice at the end of the eighties, the proposal of Anselin allows for a very limited heterogeneity. In fact, the regression coefficients are assumed to be the same across individuals, and unobserved effects are excluded. These constraints have been overcomed in the more recent works of Mur, López, and Herrera (2010) and Anselin (2016), where spatial and nonspatial parameters can take on different values in different cross-section.

Our case assumes that the model may contain G equations, with spatial structure. Additionally, it is assumed that the cross-sectional dimension (N, number of individuals) is large, and that the time dimension (Tm, number of periods) may increase. Rey and Montouri (1999), Egger and Pfaffermayr (2004), Le Gallo and DallErba (2006), Fingleton (2007), Moscone, Tosetti, and Knapp (2007), Le Gallo and Chasco (2008), Lauridsen et al. (2010), Kakamu, Polasek, and Wago (2012) or Bech et al. (2012) are well-known examples in this line. If the number of individuals is small and they are observed for a long time period (finite N and large Tm), the inference should be based on the time dimension, for example, specifying an equation for each spatial unit. Arora and Brown (1977), Hordijk and Nijkamp (1977) and White and Hewings (1982) follow the last approach.

A significant contribution is Wang and Kockelman (2007), who allow for greater heterogeneity among the individuals by introducing unobserved random effects. They extend the equation of the errors by including a mechanism of error components ‘a la KKP’, from Kapoor, Kelejian, and Prucha (2007), where the unobserved effects are intermingled with the random term of the spatial errors. This scheme is different from the proposal of Baltagi, Song, and Koh (2003), BSK, where the error is the sum of a spatial error mechanism plus an unobserved individual effect. Another point of interest in the work of Wang and Kockelman is that they built a standard SUR model, with G equations, Tm cross-sections and N individuals. Baltagi and Bresson (2011) extend the results of Wang and Kockelman by developing a collection of misspecification tests, under a maximum-likelihood, ML, framework, whereas Baltagi and Pirotte (2011) introduce generalized method of moments, GMM, in a spatial SUR with spatial errors, SEM, also ‘a la KKP’. Finally López, Mur, and Angulo (2014) address the problem of selecting the best specification for spatial SUR models, using a battery of Lagrange Multipliers obtained in a ML framework. Lundberg (2006), Zhou and Kockelman (2009), Cotteleer and Kooten (2012), Izón et al. (2016) and López, Martínez-Ortiz, and Cegarra-Navarro (2017) present nice examples in this direction.

From the discussion above, it is clear that there is a growing interest in the case of seemingly unrelated regressions in a spatial context, although the lack of specific software has slowed its adoption in applied research. Our purpose is to introduce a new R package, called spsur, to fill this gap. The preferred framework is a large N and an increasing Tm. The core of the package lies in ML estimation because, at this time, we value exactitude and reliability over quickness. There is no need to recall that the hypothesis of normality plays a key role in this approach, which may be difficult to admit in some cases. For that reason, and also to allow for speed of calculus in large samples, the ML approach is completed with a module devoted to IV procedures.

Among the limitations, let us note that we rely explicitly on the assumptions of linearity and dynamic relations are not allowed. This version of the package can remove unobserved individual effects from the equations, which are not estimated. Finally, we require the user to specify the weighting matrix, which is assumed exogenous, known, time invariant and constant across equations.

Second Section introduces notation and reviews some basic results in the literature on spatial SUR models, using the approach of Anselin. Section 3 presents the structure and some of the main functions included in spsur. Section 4 focusses on misspecification tests included in spsur. Section 5 is devoted to the issue of spatial multipliers whereas Section 6 presents the general case of a SUR model with G equations, Tm time periods and N individuals. Section 7 introduces an additional functionality that allows for simulation experiments. Finally, Section 8 concludes with a brief summary and prospects for the future.

2 SUR and Spatial SUR models

The SUR model consist of several equations, possibly with different regressors, where the error terms are correlated among equations. We focus on the Anselin’s case with a single equation, G=1, several time periods, Tm, and N individuals in the cross-section; we assume that N>Tm (note that the discussion is almost similar if Tm=1 and G>1). Our baseline model, without spatial effects, is the following: \[\begin{equation} {{\bf y}_t}={{\bf X}_t} {\boldsymbol{\beta} }_t+{\boldsymbol{\varepsilon}}_t\ ; \ \ \ \ \ \ E[{\boldsymbol \varepsilon}_t]=0\ ;\ \ \ \ \ E[{\boldsymbol \varepsilon}_t {\boldsymbol \varepsilon^{T}}_s]= \sigma_{ts} {\bf I}_N \ \ \ \ \ \ t,s=1,...,Tm \tag{2.1} \end{equation}\] where \({\bf y_t}\), \({\bf u}_t\) and \({\boldsymbol \varepsilon}_t\) are N\(\times\)1 vectors, \({\bf X}_t\) is a N\(\times p_t\) matrix, with \(p_t\) the number of regressors that appear in the t-th equation and \(\boldsymbol{\beta}_t\) the \((p_t \times 1)\) vector of coefficients. We call this model as SUR-SIM, from spatially independent model following the terminology in López, Mur, and Angulo (2014).

The singular aspect of (2.1) is that the SUR structure appears because there is intra-individual serial dependence; inter-individual serial dependence is excluded. The serial dependence in the errors is not parameterized, but estimated in the (Tm\(\times\)Tm) sampling covariance matrix, \(\boldsymbol \Sigma\)=\((\sigma_{ts})\).

The equation of (2.1) may be enough in some circumstances but not in a spatial setting where a potential problem, often present in applied research, is the existence of spatial dependence. The spatial mechanisms may have different forms (Elhorst 2014). A general SUR model that include most of the possible spatial effects of interest is:

\[\begin{equation} {\bf y}_t=\lambda_t {\bf W}_t^y {\bf y}_t+ {\bf X}_t {\boldsymbol \beta}_t+ {\bf W}_t^x {\bf X}_t^* {\boldsymbol \gamma}_t + {\bf u}_t;\ {\bf u}_t=\rho_t {\bf W}_t^u {\bf u}_t +{\boldsymbol \varepsilon}_t\ ; \ \ E[{\boldsymbol \varepsilon}_t]=0\ ;\ \ E[{{\boldsymbol \varepsilon}_t} {{\boldsymbol \varepsilon^{T}}_s}]=\sigma_{ts}{\bf I}_N \tag{2.2} \end{equation}\]

where \({\bf X}_t^*\) is the matrix of regressors excluded the intercept term; \({\bf W}_t^y\), \({\bf W}_t^x\), \({\bf W}_t^u\) are (N\(\times\)N) weighting matrices. The model in (2.2) is often referred to as a General Nesting Model, SUR-GNM. Several constrained specifications emerge from this equation, such as:

1. The SUR-SIM, spatially independent model, \(\lambda_t = \rho_t = \gamma_t = 0, (\forall t)\).
2. The SUR-SLX, spatial lags in Xs model, where \(\lambda_t = \rho_t = 0, (\forall t)\).
3. The SUR-SLM, spatial lag model, where \(\rho_t=\gamma_t = 0, (\forall t)\).
4. The SUR-SEM, spatial error model, where \(\lambda_t=\gamma_t = 0, (\forall t)\).
5. The SUR-SDM, spatial Durbin model, where \(\rho_t = 0, (\forall t)\).
6. The SUR-SDEM, spatial Durbin error model, where \(\lambda_t= 0, (\forall t)\).
7. The SUR-SARAR, spatial lag model with autoregresive errors, where \(\gamma_t = 0, (\forall t)\).

In the applied literature, it is usual to assume that the weighting matrices are the same \({\bf W}_t^x\) = \({\bf W}_t^u\) = \({\bf W}_t^{\varepsilon}\) = \({\bf W}\). The user should supply this matrix, based on apriori knowledge, which, in any case, must conform to he requisites described in Kelejian and Prucha (2004): the terms on the main diagonal are zero and the row and column sums of \(\textbf W\) are uniformly bounded in absolute value, the same as \(\textbf A_{t}^{-1} =\mathop{\left({\bf I}_N -\lambda _{t} {\bf W}\right)}\nolimits^{-1}\) and \(\textbf B_{t}^{-1} =\mathop{\left({\bf I}_N -\rho _{t} {\bf W}\right)}\nolimits^{-1}\) which must exist.

The model of (2.2), for the case of G=1, can be expressed using obvious matrix notation:

\[\begin{equation} \begin{array}{ll} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\bf Ay}={\bf X}{\boldsymbol \beta} + ({\bf I}_{Tm} \otimes {\bf W}){\bf X^*} {\boldsymbol \gamma} +{\bf u}; \ \ {\bf Bu}={\boldsymbol \varepsilon}; \ \ {\boldsymbol \varepsilon} \sim N(0,{\boldsymbol \Omega} ) \\ \\ {{\bf y}=\mathop{\left[\begin{array}{c} {\mathop{\bf y}\nolimits_{1} } \\ {\mathop{\bf y}\nolimits_{2} } \\ {\cdots } \\ {\mathop{\bf y}\nolimits_{Tm} } \end{array}\right]}\limits_{TmN\times1} ; {\bf X} =\mathop{\left[\begin{array}{cccc} {\mathop{\bf X}\nolimits_{1} } & {0} & {\cdots } & {0} \\ {0} & {\mathop{\bf X}\nolimits_{2} } & {\cdots } & {0} \\ {\cdots } & {\cdots } & {\cdots } & {\cdots } \\ {0} & {0} & {\cdots } & {\mathop{\bf X}\nolimits_{Tm} } \end{array}\right]}\limits_{TmN\times p} ;\, \, {\boldsymbol \beta} =\mathop{\left[\begin{array}{c} {\mathop{{\boldsymbol \beta}}\nolimits_{1} } \\ {\mathop{{\boldsymbol \beta}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\boldsymbol \beta}}\nolimits_{Tm} } \end{array}\right]}\limits_{p\times 1} ;{\bf u}=\mathop{\left[\begin{array}{c} {\mathop{{\bf u}}\nolimits_{1} } \\ {\mathop{{\bf u}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\bf u}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmN\times1} ;\, \, {\boldsymbol \varepsilon} =\mathop{\left[\begin{array}{c} {\mathop{{\boldsymbol \varepsilon}}\nolimits_{1} } \\ {\mathop{{\boldsymbol \varepsilon}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\boldsymbol \varepsilon}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmN\times1} } \end{array} \tag{2.3} \end{equation}\]

where \(p=\sum p_t\), \({\bf A}={\bf I}_{TmN}-{\boldsymbol \Lambda} \otimes {\bf W}\) with \({\boldsymbol \Lambda}=diag(\lambda_1,\dots,\lambda_{Tm})\); \({\bf B}={\bf I}_{TmN}- {\boldsymbol \Gamma} \otimes {\bf W}\) with \({\boldsymbol \Gamma}=diag(\rho_1,\dots,\rho_{Tm})\) and \({\boldsymbol \Omega}={\boldsymbol \Sigma} \otimes {\bf I}_{TmN}\), being \(\otimes\) the Kronecker product. Moreover, \({\bf X}^*\) is a NTm \(\times\) (p-Tm) that omits the columns of the intercepts because these terms cannot be spatially lagged. Assuming that the errors are normally distributed, the log-likelihood function of (2.2) can be written as:

\[\begin{multline} l({\bf y};{\boldsymbol \eta})=-\frac{NTm}{2} \ln (2\pi)-\frac{N}{2} \ln |{\bf \Sigma}|+Tm\sum _{t=1}^{Tm}\ln|{\bf B}_t| +Tm\sum _{t=1}^{Tm}\ln|{\bf A}_t| \\ -\frac{\left({\bf Ay}-{\bf \bar{X}}[{\boldsymbol \beta} \ {\boldsymbol \gamma}]^{T} \right)^{T}{\bf B}^{T}\mathop{{\bf \Omega}}\nolimits^{-1} {\bf B}\left({\bf Ay}-{\bf \bar{X}}[{\boldsymbol \beta} \ {\boldsymbol \gamma}]^{T} \right)}{2} \tag{2.4} \end{multline}\]

where \(\bar{\textbf X}=[\textbf X,\bf W {\bf X}^*]\) is a NTm \(\times\) (2p-Tm).

The vector of parameters is \({\boldsymbol \eta}^{T}=\left[{\boldsymbol \beta}^{T};{\boldsymbol \gamma}^{T};\lambda_{1}; \dots ;\lambda_{Tm};\rho_{1};\dots ;\rho_{Tm} ;\mathop{\sigma }\nolimits_{ij} \right]\), where there are P= (2p-Tm)+Tm+Tm(Tm+1)/2 unknown coefficients Wang and Kockelman (2007) and López, Mur, and Angulo (2014) use numerical optimization techniques to solve the maximum likelihood estimates. Note that the presence of the log of the Jacobian terms, \(\ln|{\bf B}_t|\) and \(\ln|{\bf A}_t|\), forces the spatial parameters, \(\lambda_t\) and \(\rho_t\), to lie inside the so-called stability interval (if the weighting matrix is row-standardized, a rough approximation is -1,+1). Under the assumption that the spatial SUR model is correctly specified (including normality), the ML estimators are consistent, efficient and asymptotically normally distributed (Davidson and MacKinnon 1993).

3 Introducing the spsur package

spsur is an open-source software for the R computing platform (R Core Team 2018). spsur depends on the packages Formula, MASS, Matrix, methods, numDeriv, minqa, sparseMVN, spdep and stats. The package can be freely dowloaded from CRAN and a development version is available in the link Github-spsur. The users will find a complete helping documentation and a vignette to guide them through the functionalities of the package.

spsur is made of 10 different functions which allow for a thoroughful analysis of a spatial SUR specification. Figure 1 shows the main functionalities of the package. spsur includes one function to evaluate the Lagrange Multipliers for the presence of spatial correlation in a SUR-SIM model (lmtestspsur()); two functions to estimate spatial SUR models, by maximum likelihood (spsurml()) or using instrumental variables (spsur3sls()) for the general case of G>1 (including both Tm=1 and Tm>1), and a third one designed to deal with pure spatial panels (spsurtime()) for G=1 and Tm>1; four additional functions allow to improve the specification (lrtestspsur(); wald_betas(), wald_deltas() and lr_betas_spsur()); another function obtains the spatial impacts of the regressors (impacts()). Finally dgp_spsur() allows the user to solve Monte Carlo experiments using spatial SUR models.

Figure 3.1: Main functionatities of spsur package

# install.packages("spsur")

library(spsur)
data("spc", package = "spsur")

spcformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA

data(NCOVR, package = "spsur")
NCOVRformula <- HR80 | DV80 | FP79 ~ PS80 + UE80 | PS80 + UE80 + SOUTH | PS80

LMs.spc <- lmtestspsur(Form = spcformula, data = spc, W = Wspc)
#>              LM-Stat. DF p-value  
#> LM-SUR-SLM     5.2472  2  0.0725 .
#> LM-SUR-SEM     3.3050  2  0.1916  
#> LM*-SUR-SLM    2.1050  2  0.3491  
#> LM*-SUR-SEM    0.1628  2  0.9218  
#> LM-SUR-SARAR   5.7703  4  0.2170  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LMs.NCOVR <- lmtestspsur(Form = NCOVRformula, data = NCOVR, W = W)
#>              LM-Stat. DF  p-value    
#> LM-SUR-SLM    5494.02  3  < 2e-16 ***
#> LM-SUR-SEM    5620.26  3  < 2e-16 ***
#> LM*-SUR-SLM     53.16  3 1.69e-11 ***
#> LM*-SUR-SEM    179.41  3  < 2e-16 ***
#> LM-SUR-SARAR  5741.17  6  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Y <- as.matrix(c(NCOVR$HR80,NCOVR$DV80,NCOVR$FP79))
Intercep <- matrix(1,ncol = 1,nrow = 3085)
X1 <- cbind(Intercep,NCOVR$PS80,NCOVR$UE80)
X2 <- cbind(Intercep,NCOVR$PS80,NCOVR$UE80,NCOVR$SOUTH)
X3 <- cbind(Intercep,NCOVR$PS80)
X <- as.matrix(Matrix::bdiag(X1,X2,X3))

LMs.NCOVR.matrix <- lmtestspsur(Y = Y, X = X, G = 3, Tm = 1, N = 3085, W = W)

spcsur.slm <- spsurml(Form = spcformula, data = spc, type = 'slm', W = Wspc)
#> Initial point:   log_lik:  113.197  lambdas:  -0.472 -0.446 
#> Iteration:  1   log_lik:  114.085  lambdas:  -0.506 -0.482 
#> Iteration:  2   log_lik:  114.096  lambdas:  -0.506 -0.482 
#> Time to fit the model:  1.64  seconds
#> Computing marginal test... 
#> Time to compute covariances:  0.23  seconds

summary(spcsur.slm)
#> Call:
#> spsurml(Form = spcformula, data = spc, W = Wspc, type = "slm")
#> 
#>  
#> Spatial SUR model type:  slm 
#> 
#> Equation  1 
#>                 Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_1  1.4955217  0.2467240  6.0615 5.183e-07 ***
#> UN83_1         0.8070029  0.2557439  3.1555  0.003179 ** 
#> NMR83_1       -0.5194114  0.2590550 -2.0050  0.052318 .  
#> SMSA_1        -0.0073247  0.0118519 -0.6180  0.540347    
#> lambda_1      -0.5057334  0.2405734 -2.1022  0.042401 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.6224 
#>   Equation  2 
#>                 Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_2  1.7094414  0.2925620  5.8430 1.024e-06 ***
#> UN80_2        -0.6745562  0.3870737 -1.7427   0.08969 .  
#> NMR80_2        0.7502934  0.3842670  1.9525   0.05847 .  
#> SMSA_2         0.0014181  0.0241859  0.0586   0.95356    
#> lambda_2      -0.4821428  0.2557758 -1.8850   0.06730 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.4743 
#>   Variance-Covariance Matrix of inter-equation residuals:                            
#>   0.0003085954 -0.0003561928
#>  -0.0003561928  0.0015864976
#> Correlation Matrix of inter-equation residuals:                    
#>   1.000000 -0.509062
#>  -0.509062  1.000000
#> 
#>  R-sq. pooled: 0.6603 
#>  Log-Likelihood:  114.096
#>  Breusch-Pagan: 6.516  p-value: (0.0107) 
#>  LMM: 0.50489  p-value: (0.477)

spcSUR.sim <- spsurml(Form = spcformula, data = spc, type = 'sim', W = Wspc)
spcSUR.slx <- spsurml(Form = spcformula, data = spc, type = 'slx', W = Wspc)
spcSUR.sem <- spsurml(Form = spcformula, data = spc, type = 'sem', W = Wspc)
spcSUR.sarar <- spsurml(Form = spcformula, data = spc, type = 'sarar', W = Wspc)
spcSUR.sdm <- spsurml(Form = spcformula, data = spc, type = 'sdm', W = Wspc)
spcSUR.sdem <- spsurml(Form = spcformula, data = spc, type = 'sdem', W = Wspc)

mlcontrol <- list(tol = 0.1, maxit = 20, trace = FALSE)
NCOVRSUR.slm <- spsurml(Form = NCOVRformula, data = NCOVR,
                      type = 'slm', W = W, control = mlcontrol, cov = FALSE)

mlcontrol <- list(tol = 0.1, maxit = 20, trace = FALSE)
NCOVRSUR.slm.matrix <- spsurml(Y = Y, X = X, G = 3, Tm = 1, N = 3085,
                             W = W, type = 'slm', control = mlcontrol, 
                             cov = FALSE)

NCOVRSUR.slm.3sls <-spsur3sls(Form = NCOVRformula, data = NCOVR, type = "slm", W = W)
#> Time to fit the model:  0.17  seconds
summary(NCOVRSUR.slm.3sls)
#> Call:
#> spsur3sls(Form = NCOVRformula, data = NCOVR, W = W, type = "slm")
#> 
#>  
#> Spatial SUR model type:  slm 
#> 
#> Equation  1 
#>                Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_1  8.877617   1.307821  6.7881 1.206e-11 ***
#> PS80_1         0.892351   0.138898  6.4245 1.388e-10 ***
#> UE80_1        -0.060245   0.033200 -1.8146   0.06962 .  
#> lambda_1      -0.221472   0.191153 -1.1586   0.24664    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.1405 
#>   Equation  2 
#>                Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_2 2.9935206  0.3798303  7.8812  3.61e-15 ***
#> PS80_2        0.2533626  0.0262643  9.6467 < 2.2e-16 ***
#> UE80_2        0.0914432  0.0095959  9.5294 < 2.2e-16 ***
#> SOUTH_2       0.0346575  0.0471028  0.7358    0.4619    
#> lambda_2      0.2124110  0.0900048  2.3600    0.0183 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.291 
#>   Equation  3 
#>                Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_3  7.670117   1.095699  7.0002 2.734e-12 ***
#> PS80_3        -1.585681   0.160403 -9.8856 < 2.2e-16 ***
#> lambda_3       0.386346   0.087474  4.4167 1.014e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.5611 
#>   Variance-Covariance Matrix of inter-equation residuals:                              
#>  45.377766  1.867308 18.961586
#>   1.867308  1.938112 -1.427426
#>  18.961586 -1.427426 34.820491
#> Correlation Matrix of inter-equation residuals:                                
#>  1.0000000  0.1991154  0.4770190
#>  0.1991154  1.0000000 -0.1737588
#>  0.4770190 -0.1737588  1.0000000
#> 
#>  R-sq. pooled: 0.3595

data(NCOVR,package="spsur")
N <- nrow(NCOVR)
Tm <- 4
index_time <- rep(1:Tm, each = N)
index_indiv <- rep(1:N, Tm)
pHR <- c(NCOVR$HR60, NCOVR$HR70, NCOVR$HR80, NCOVR$HR90)
pPS <- c(NCOVR$PS60, NCOVR$PS70, NCOVR$PS80, NCOVR$PS90)
pUE <- c(NCOVR$UE60, NCOVR$UE70, NCOVR$UE80, NCOVR$UE90)
pNCOVR <- data.frame(indiv = index_indiv, time = index_time,
                      HR = pHR, PS = pPS, UE = pUE)
rm(NCOVR,pHR,pPS,pUE,index_time,index_indiv)
form_pHR <- HR ~ PS + UE
pHR_slm <- spsurtime(Form = form_pHR, data = pNCOVR, W = W,
                     time = pNCOVR$time, type = "slm", method = "3sls")
#> Time to fit the model:  0.29  seconds
summary(pHR_slm)
#> Call:
#> spsur3sls(W = W, X = X, Y = Y, G = G, N = N, Tm = Tm, p = p, 
#>     type = type, maxlagW = maxlagW)
#> 
#>  
#> Spatial SUR model type:  slm 
#> 
#> Equation  1 
#>                Estimate Std. Error t value Pr(>|t|)
#> (Intercept)_1  2.501017   1.645140  1.5202   0.1285
#> PS_1          -0.029887   0.105860 -0.2823   0.7777
#> UE_1           0.029898   0.046422  0.6441   0.5196
#> lambda_1       0.405674   0.392876  1.0326   0.3018
#> R-squared: 0.2843 
#>   Equation  2 
#>                Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)_2 1.8871930  1.9497433  0.9679  0.33310  
#> PS_2          0.1705918  0.1651996  1.0326  0.30179  
#> UE_2          0.0099723  0.0497241  0.2006  0.84105  
#> lambda_2      0.6971788  0.2942217  2.3696  0.01782 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.3604 
#>   Equation  3 
#>                Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_3  6.881003   1.220794  5.6365 1.774e-08 ***
#> PS_3           0.658929   0.135891  4.8489 1.256e-06 ***
#> UE_3           0.183858   0.032266  5.6982 1.239e-08 ***
#> lambda_3      -0.172457   0.178288 -0.9673    0.3334    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.06152 
#>   Equation  4 
#>               Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept)_4 3.917895   0.369609 10.6001 < 2.2e-16 ***
#> PS_4          0.909703   0.119384  7.6200 2.723e-14 ***
#> UE_4          0.319217   0.042374  7.5332 5.295e-14 ***
#> lambda_4      0.023133   0.073747  0.3137    0.7538    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> R-squared: 0.134 
#>   Variance-Covariance Matrix of inter-equation residuals:                                    
#>  31.72416 20.95984 16.31181 13.23690
#>  20.95984 53.75722 26.29607 19.95032
#>  16.31181 26.29607 45.37777 19.69094
#>  13.23690 19.95032 19.69094 38.89755
#> Correlation Matrix of inter-equation residuals:                                        
#>  1.0000000 0.5075448 0.4299180 0.3768166
#>  0.5075448 1.0000000 0.5324158 0.4362851
#>  0.4299180 0.5324158 1.0000000 0.4686880
#>  0.3768166 0.4362851 0.4686880 1.0000000
#> 
#>  R-sq. pooled: 0.1608

4 Misspecification tests and other utilities

The estimation of spatial SUR models needs of a series of assumptions that are critical to guarantee the good properties of the (ML, 3SLS) estimates. Some of the tests required to build the spatial structure of the model have been presented in Section 2.1. Below we add a series of misspecification tests, implemented in spsur, that may help to improve the specification.

## Not included the results of execution in vignette to save time.
LRtests <- lrtestspsur(Form = spcformula, data = spc, W = Wspc)

R1 <- matrix(c(1,0,0,0,-1,0,0,0), nrow = 1)
b1 <- matrix(0, ncol = 1)
Wald_beta <- wald_betas(results = spcsur.slm, R = R1, b = b1)
#> Wald stat.: 0.396 p-value: (0.529)

## Not included the results of execution in vignette to save time.
R1 <- matrix(c(1,0,0,0,-1,0,0,0),nrow=1)
b1 <- matrix(0, ncol = 1)
spcSUR.slm.restricted <- spsurml(Form = spcformula, data = spc, 
               type = "slm", W = Wspc, R = R1, b = b1, control=list(tol=0.01,maxit=200,trace=FALSE))

summary(spcSUR.slm.restricted)

## Not included the results of execution in vignette to save time.
R2 <- matrix(c(1,-1), nrow = 1)
b2 <- matrix(0, ncol = 1)
Wald_lambda <- wald_deltas(results=spcSUR.slm.restricted, R = R2, b = b2)

## Not included the results of execution in vignette to save time.
R1 <- matrix(c(1,0,0,0,-1,0,0,0), nrow = 1)
b1 <- matrix(0, ncol = 1)
LR_SMSA <- lr_betas_spsur(Form = spcformula, data = spc, W = Wspc,
                           type = "slm", R = R1, b = b1, trace = FALSE,
                           printmodels = FALSE)

5 Direct, Indirect and Total impacts

According to LeSage and Pace (2009), a major difference between times series models and spatial models consists in the interpretation of the \(\beta\) parameters. In a time series model, the \(\beta\) parameters have a 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, \(\dfrac{\partial y_{t}}{\partial x_{t}} = \beta_{k}\), and its significance can be checked by a simple tratio.

Things are not so clear in a spatial setting because of the feed-back effects. In fact, a change in a regressor, in a given location, will induce a chain of reactions that will spill over into the space, depending on the weights matrix, \(\bf W\). Consider, for example, the reduced form of a SUR-SARAR model:

\[\begin{equation} {\bf y}={\bf A}^{-1}{\bf X}{\boldsymbol \beta}+{\bf A}^{-1}{\bf B}^{-1}{\boldsymbol \varepsilon}={\bf I}_{Tm}\otimes\left[{\bf I}_{N}+{\boldsymbol \Lambda} \otimes {\bf W}+{\boldsymbol \Lambda}^{2}\otimes {\bf W}^{2}+\ldots \right] {\bf X} {\boldsymbol \beta}+{\bf A}^{-1}{\bf B}^{-1}{\boldsymbol \varepsilon} \tag{5.1} \end{equation}\]

The expression above shows that a change in the value of a regressor, in certain point in space, will impact the explained variable in the same location and also in other points of space. Note that \({\bf A}\) is a block-diagonal matrix whose non-zero matrices are \({\bf A}_{t}={\bf I}_{N}-\lambda_{t}\textbf{W}; t=1,2, \ldots Tm\). The inverse of \({\bf A}_{t}\) admits a power expansion on the weights matrix, , and the spatial parameters \(\lambda_{t}\), which lead us to the result in (5.1).

Turning to the general case of G>1 and Tm>1 and using the notation in LeSage and Pace (2009), we may write the g-th equation of (5.1), in period t, as:

\[\begin{equation} {\mathop{{\bf y}}\nolimits_{tg} }=\mathop{\left[\begin{array}{c} {\mathop{y}\nolimits_{tg1} } \\ {\mathop{y}\nolimits_{tg2} } \\ {\cdots } \\ {\mathop{y}\nolimits_{tgN} } \end{array}\right]} = \sum_{h=1}^{p_{q}} \mathop{\left[\begin{array}{cccc} {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{11}} & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{12}} & {\cdots } & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{1N}} \\ {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{21}} & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{22}} & {\cdots } & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{2N}} \\ {\cdots } & {\cdots } & {\cdots } & {\cdots } \\ {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{N1}} & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{N2}} & {\cdots } & {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{NN}} \end{array}\right]} \mathop{\left[\begin{array}{c} {\mathop{x}\nolimits_{tg1} } \\ {\mathop{x}\nolimits_{tg2} } \\ {\cdots } \\ {\mathop{x}\nolimits_{tgN} } \end{array}\right]}+ {\mathop{{\bf A}}\nolimits_{g}^{-1}} {\mathop{{\bf B}}\nolimits_{g}^{-1}} \mathop{\left[\begin{array}{c} {\mathop{\varepsilon}\nolimits_{tg1} } \\ {\mathop{\varepsilon}\nolimits_{tg2} } \\ {\cdots } \\ {\mathop{\varepsilon}\nolimits_{tgN} } \end{array}\right]} \tag{5.2} \end{equation}\]

being \({\mathop{{\bf S}}\nolimits^{h}_{tg}({\bf W})}= {\mathop{\beta}\nolimits_{h}} {\mathop{{\bf A}}\nolimits_{g}^{-1}}\) and \({\mathop{S}\nolimits^{h}_{tg}({\bf W})_{ij}}\) its \(i,j-th\) element. From above, it is clear that the value of the partial derivative of \(y\) with respect to \(x_{h}\) depends on the point in space that we are impacting; that is: \(\dfrac{\partial y_{itg}}{\partial x_{h_{jtg}}} = {\mathop{S}\nolimits^{h}_{tg}({\bf W})_{ij}}\). Moreover, it is evident, form (5.2), that a shock in a point will also spread out through the space, but the discussion about spatial multipliers excludes the case of impacting the error terms.

Specifically, for each equation and every regressor, we can evaluate three effects, such as:

\[\begin{equation} \begin{array}{c} {\mathop{{\bf M}}\nolimits_{Total_{g}^{h}}} = \dfrac{1}{N} \tau^{T}_{N} {\bf S}^{h}_{.g} \tau_{N} \ \ \ \\ {\mathop{{\bf M}}\nolimits_{Direct_{g}^{h}}} = \dfrac{1}{N} tr {\bf S}^{h}_{.g}\\ {\mathop{{\bf M}}\nolimits_{Indirect_{g}^{h}}} = {\mathop{{\bf M}}\nolimits_{Total_{g}^{h}}} -{\mathop{{\bf M}}\nolimits_{Direct_{g}^{h}}} \\ \\ h=1,\ldots p_{g}; g=1 \ldots G \end{array} \tag{5.3} \end{equation}\]

being \(\tau_{N}\) a N\(\times\)1 unitary vector and \({\bf S}^{h}_{.g}= \sum_{t=1}^{Tm} {\bf S}^{h}_{tg}({\bf W})\).

Under the assumption that the spatial parameters and the weighting matrices are constant across time, the partial derivatives in (5.2) are also constant. Futhermore, it is important to recall that the models in spsur have no time dynamics, so the impact of a modification in a regressors must be absorved in the same time period.

Finally, we can aggregate the multipliers in (5.3) to obtain overall measures for each variable:

\[\begin{equation} \begin{array}{c} {\mathop{{\bf M}}\nolimits_{Total_{\bullet}^{h}}} = \sum_{g=1}^{G} {\mathop{{\bf M}}\nolimits_{Total_{g}^{h}}} \\ {\mathop{{\bf M}}\nolimits_{Direct_{\bullet}^{h}}} = \sum_{g=1}^{G} {\mathop{{\bf M}}\nolimits_{Direct_{g}^{h}}} \\ {\mathop{{\bf M}}\nolimits_{Indirect_{\bullet}^{h}}} = \sum_{g=1}^{G} {\mathop{{\bf M}}\nolimits_{Indirect_{g}^{h}}} \\ \\ h=1,\ldots p_{g} \end{array} \tag{5.4} \end{equation}\]

The multipliers of (5.4) should be accompanied by measures of statistical significance to be useful for the user. LeSage and Pace (2009) suggest two possibilities. First is the Bayesian Markov Chain Monte Carlo method, MCMC, which draws samples of parameters from their posterior distribution and computes random sequences of multipliers, to calibrate the estimates. However, to our knowledge, there are no results available for spatial SUR models using MCMC. The second alternative is bootstrapping. In the case of spsur, we prefer a parametric bootstrap, which means that:

• We estimate the spatial SUR model by ML.
• We draw B samples of random vectors of the parameters in our spatial SUR, \({\boldsymbol \eta}\), using a multivariate normal distribution, \(N \left(\hat{{\boldsymbol \eta}}; \hat{{\bf V}}(\hat{{\boldsymbol \eta}})\right)\), where \(\hat{{\boldsymbol \eta}}\) is the ML estimation of the parameters and \(\hat{{\bf V}}(\hat{{\boldsymbol \eta}})\) their corresponding estimated covariance matrix. Let denote by \({\boldsymbol \eta}_b\) the b-th draw.
• We re-evaluate the effects using the boostrapped vectors \({\boldsymbol \eta}_b\), maintaining the data of the regressors. In this way, a sequence of values: \({\mathop{M}\nolimits_{Total_{g}^{h^{(b)}}}}\), \({\mathop{M}\nolimits_{Direct_{g}^{h^{(b)}}}}\) and \({\mathop{M}\nolimits_{Indirect_{g}^{h^{(b)}}}}\) for \(b=1, \ldots B\) is obtained.
• Finally, we calculate the corresponding measure of dispersion as the standard deviation of the series of estimated multipliers:

\[\begin{equation} {\mathop{\sigma}\nolimits_{M_{Effect_{g}^{h}}}} = \sqrt{ \dfrac{1}{B} \sum_{b=1}^B \left( \mathop{M}\nolimits_{Effect_{g}^{h^{(b)}}} - \overline{\mathop{M}\nolimits_{Effect_{g}^{b}}} \right)^{2}} \end{equation}\]

with Effect=Total, Direct, Indirect.

## Not included the results of execution in vignette to save time.
spcsur.slm.impacts <- impacts(spcsur.slm)

6 A general case: G>1, Tm>1 and N>1

The framework of the previous Sections can be generalized quite easily to a case where the user has data for more than one single equation; that is, G>1. The equivalent to the GNM model of expression (2.2) is:

\[\begin{align} {\bf y}_{tg} &= \lambda_{g} {\bf W}{\bf y}_{tg} + {\bf X}_{tg} {\boldsymbol \beta}_{g} + {\bf W} {\bf X}_{tg}^* {\boldsymbol \gamma}_{g} + {\bf u}_{tg} \\ {\bf u}_{tg} &= \rho_{g} {\bf W} {\bf u}_{tg} + {\boldsymbol \varepsilon}_{tg} \nonumber \tag{6.1} \end{align}\]

\[\begin{equation*} E \left[ {\boldsymbol \varepsilon}_{tg} \right] = 0 \quad E \left[ {\boldsymbol \varepsilon}_{tg} {\boldsymbol \varepsilon}^{T}_{sh} \right] = \left\{ \begin{array}{ll} \sigma_{gh} {\bf I}_{N} & t = s \\ {\bf 0}_{N} & t \ne s \\ \end{array} \right. \end{equation*}\]

for g=1,\(\ldots\)G; t=1,\(\ldots\)Tm. The SUR nature of the model appears in the form of contemporaneous correlation among the errors of the different equations, corresponding to the same individual. This correlation is zero for errors belonging to equations of different individuals and/or different time periods. Serial correlation is only allowed in Anselin’s approach (treated with the function spsurtime(), as shown before). If we stack the G vectors of errors corresponding to period t in a (NG\(\times\)1) vector \({\boldsymbol \varepsilon}_{t}\), its covariance matrix can be written as \(E\left[\mathop{{\boldsymbol \varepsilon}}\nolimits_{t} \mathop{{\boldsymbol \varepsilon}}\nolimits_{t}^{T} \right]= {\boldsymbol \Sigma} \otimes \mathop{{\bf I}}\nolimits_{N}\); this is a (NG\(\times\)NG) matrix. In continuation, we can stack these Tm vectors into a general (TmNG\(\times\)1) vector by ordering the observations, first, by time period, then by equation and finally by individuals. Using obvious matrix notation, we obtain:

\[\begin{equation} \begin{array}{l} \left. \begin{array}{l} {\, \begin{array}{ll} {\bf Ay}={\bf X}{\boldsymbol \beta} +({\bf I}_{Tm} \otimes {\bf I}_{G} \otimes {W}) {\bf X}^* {\boldsymbol \gamma}+{\bf u} & {} \\ {{\bf Bu}={\boldsymbol \varepsilon};\, \, {{\boldsymbol \varepsilon} \sim (0,{\boldsymbol \Omega })} } \end{array}} \end{array}\right\} \\ {} \\ {\bf y}=\mathop{\left[\begin{array}{c} {\mathop{{\bf y}}\nolimits_{1} } \\ {\mathop{{\bf y}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\bf y}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmGN\times1} ; \, \, \, {\bf X}=\mathop{\left[\begin{array}{c} {\mathop{{\bf X}}\nolimits_{1} } \\ {\mathop{{\bf X}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\bf X}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmGN\times p} ; \, \, \, {\bf X}^*=\mathop{\left[\begin{array}{c} {\mathop{\bf X}^*_{1}} \\ {\mathop{\bf X}^*_{2}} \\ {\cdots } \\ {\mathop{\bf X}^*_{Tm}} \end{array}\right]}\limits_{TmGN\times (p-G)} ; \, \, \, {\bf u}=\mathop{\left[\begin{array}{c} {\mathop{{\bf u}}\nolimits_{1} } \\ {\mathop{{\bf u}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\bf u}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmGN\times 1} ; \, \, \, {\boldsymbol \varepsilon} =\mathop{\left[\begin{array}{c} {\mathop{{\boldsymbol \varepsilon}}\nolimits_{1} } \\ {\mathop{{\boldsymbol \varepsilon}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\boldsymbol \varepsilon}}\nolimits_{Tm} } \end{array}\right]}\limits_{TmGN\times 1};\\ {} \\ \mathop{{\bf X}}\nolimits_{t} =\mathop{\left[\begin{array}{cccc} {\mathop{{\bf X}}\nolimits_{1t} } & {0} & {\cdots } & {0} \\ {0} & {\mathop{{\bf X}}\nolimits_{2t} } & {\cdots } & {0} \\ {\cdots } & {\cdots } & {\cdots } & {\cdots } \\ {0} & {0} & {\cdots } & {\mathop{{\bf X}}\nolimits_{Gt} } \end{array}\right]}\limits_{GN\times p} ; \mathop{{\bf X}}^*_{t} =\mathop{\left[\begin{array}{cccc} \mathop{{\bf X}}^*_{1t} & {0} & {\cdots } & {0} \\ {0} & \mathop{{\bf X}}^*_{2t} & {\cdots } & {0} \\ {\cdots } & {\cdots } & {\cdots } & {\cdots } \\ {0} & {0} & {\cdots } & \mathop{{\bf X}}^*_{Gt} \end{array}\right]}\limits_{GN\times (p-G)}; {\boldsymbol \beta} =\mathop{\left[\begin{array}{c} {\mathop{{\boldsymbol \beta}}\nolimits_{1} } \\ {\mathop{{\boldsymbol \beta}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\boldsymbol \beta}}\nolimits_{G} } \end{array}\right]}\limits_{p\times 1} ; {\boldsymbol \gamma} =\mathop{\left[\begin{array}{c} {\mathop{{\boldsymbol \gamma}}\nolimits_{1} } \\ {\mathop{{\boldsymbol \gamma}}\nolimits_{2} } \\ {\cdots } \\ {\mathop{{\boldsymbol \gamma}}\nolimits_{G} } \end{array}\right]}\limits_{(p-G)\times 1} ; \end{array} \tag{6.2} \end{equation}\]

where \({\bf A}={\bf I}_{Tm} \otimes \left[{\bf I}_{GN} -{\boldsymbol \Lambda} \otimes {\bf W}\right] = {\bf I}_{Tm} \otimes diag\left\{{\bf A}_{g}; g=1,2,..,G\right\}\), and \({\bf B}={\bf I}_{Tm} \otimes \left[{\bf I}_{GN} -{\boldsymbol \Upsilon } \otimes {\bf W}\right]=diag \left\{{\bf B}_{g} ;g=1,2,..,G\right\}\), \({\boldsymbol \Lambda}\) and \({\boldsymbol \Upsilon}\) are two (G\(\times\)G) diagonal matrices containing the spatial parameters \({\lambda}_{g}\) and \({\rho}_{g}\), respectively. Moreover \({\boldsymbol \Omega} ={\bf I}_{Tm} \otimes {\boldsymbol \Sigma} \otimes {\bf I}_{N}\). The same as before, we allow for a different number of regressors in the right hand side of each equation, \(p_{g}\), where the total number of \(\beta\) parameters is: \(p=\Sigma_{g=1}^{G} p_{g}\) and p-G is the total number of \(\gamma\) parameters.

Assuming that the errors are normally distributed, the log-likelihood function of the spatial GNM-SUR model of (6.2) can be written as:

\[\begin{align} l({\bf y};{\boldsymbol \eta}) =& -\frac{NTmG}{2} \ln \left(2\pi \right)-\frac{NTm}{2} \ln \left|{\boldsymbol \Sigma} \right|+Tm\left\{\sum _{g=1}^{G}\ln \left|{\bf B}_{g} \right| +\sum _{g=1}^{G}\ln \left|{\bf A}_{g} \right| \right\} \\ & -\frac{\left\{{\bf Ay}-[{\bf X};{\bf X}_{{\bf W}}^*] {\boldsymbol \Pi} \right\}^{T}{\bf B}^{T}\mathop{{\boldsymbol \Omega} }\nolimits^{-1} {\bf B}\left\{{\bf Ay}-[{\bf X};{\bf X}_{{\bf W}}^*]{\boldsymbol \Pi} \right\}}{2} \tag{6.3} \end{align}\]

where \({\boldsymbol \eta}^{T}=\left[{\boldsymbol \Pi}^{T}; \mathop{\lambda }\nolimits_{1} ;\cdots ;\mathop{\lambda }\nolimits_{G} ;\mathop{\rho }\nolimits_{1} ;\cdots ;\mathop{\rho }\nolimits_{G} ;\mathop{\sigma }\nolimits_{ij} \right]\) is the vector parameters of the model, being \({\boldsymbol \Pi}=[{\boldsymbol \beta}; {\boldsymbol \gamma}]^{T}\). The discussion that continues from this point is totally analogous to the spatial SUR panel analyzed so far, thus we leave it now to avoid repetitions.

The example below shows the case of a SUR-SLM model, for NCOVR, where G=3, Tm=4 and N=3,085, including the reshaping of a database into a classical panel:

data(NCOVR,package="spsur")
N <- nrow(NCOVR)
Tm <- 4
index_time <- rep(1:Tm, each = N)
index_indiv <- rep(1:N, Tm)
pHR <- c(NCOVR$HR60, NCOVR$HR70, NCOVR$HR80, NCOVR$HR90)
pPS <- c(NCOVR$PS60, NCOVR$PS70, NCOVR$PS80, NCOVR$PS90)
pUE <- c(NCOVR$UE60, NCOVR$UE70, NCOVR$UE80, NCOVR$UE90)
pDV <- c(NCOVR $DV60, NCOVR$DV70, NCOVR$DV80, NCOVR$DV90)
pFP <- c(NCOVR$FP59, NCOVR$FP70, NCOVR$FP80, NCOVR$FP90)
pSOUTH <- rep(NCOVR$SOUTH, Tm)
pNCOVR <- data.frame(indiv = index_indiv, time = index_time,
                      HR = pHR, PS = pPS, UE = pUE, DV = pDV,
                      FP = pFP, SOUTH = pSOUTH)
rm(NCOVR,pHR,pPS,pUE,pDV,pFP,pSOUTH,index_time,index_indiv)
pform <- HR | DV | FP ~ PS + UE | PS + UE + SOUTH | PS

The estimation of a SLM model with G=3, Tm=4 and N=3,085 can be done as usual using the spsurml() function (to prevent overflows of memory we set cov argument equals to FALSE). The last sentence of the code just print the vector of estimates of the \(\lambda\) and \(\beta\) parameters and the correlation matrix among the residuals of the three equations:

## Not included the results of execution in vignette to save time.
 psur_slm <- spsurml(Form = pform, data = pNCOVR, W = W,
                     type = "slm", cov = FALSE, 
                  control = list(tol = 0.1, maxit = 200, trace = FALSE))
psur_slm$deltas; psur_slm$betas; psur_slm$Sigma_corr

7 How to generate random SUR datasets

This section describes an additional functionality included in spsur, which offers the possibility of generate random data sets with the dimensions and spatial structure decided by the user. This is the purpose of dgp_spsur(). The function has great flexibility to accomodate the necessities of the user and may be of interest in two different situations: (i) as part of a larger simulation experiment where the user needs random data sets with specific features (SUR nature, spatial structure, etc.) or (ii) with the aim of showing specific properties of spatial SUR models and inferential procedures related to them.

dgp_spsur() is a Data Generating Process where the user decides the dimensions of the required data set, specificies the values of the parameters that intervene in the corresponding SUR model and selects the distribution function from which the regressors and the random terms are to be drawn. The syntax of the function is the following:

dgp_spsur <- function(Sigma, Tm = 1, G, N, Betas, Thetas = NULL, durbin = FALSE, rho = NULL, lambda = NULL, p = NULL, W = NULL, X = NULL, pdfU = "nvrnorm", pdfX = "nvrnorm")

The dimensions of the data set are defined by the arguments Tm, N and G. Then the user specifies the SUR structure, begining by the covariance matrix among the residuals of the G equations, which is the core of the SUR model, with the argument Sigma. Of course, a (N\(\times\)N) spatial weighting matrix should be uploaded in the argument W.

In continuation, the user should think in 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 through row vectors of the adequate dimensions. A fundamental piece of information is the argument p that defines the number of regressors (that is, x variables) that appear in each equation; if p is a scalar, every equation has the same number of regressors whereas if we need a model with different number of regressor in each equation, p must be a (1 \(\times\) G) vector. In both cases, Betas is a row vector of order (\(1 \times G^{*}\)), where \(G^{*}\) is pG in the first case and \(\Sigma_{1}^{G}p_g\) in the second case.

Note that, depending on the values given to these arguments, the type of spatial model is unequivocally determined. For example, if Thetas=NULL, rho=NULL but lambda is not NULL, the user is specifying a SUR-SLM model or a SUR-SDEM model in the case of lambda=NULL but Thetas and rho different from NULL.

There are two possibilities to build the matrix \(\textbf{X}\) of regressors. In Monte Carlo experiments is very usual to maintain fixed the values of the regressors and repeatedly draw random matrices of error terms, to simulate the explained variable. If this is the case, the user should upload the required \(\textbf{X}\) matrix in the argument X, which must be consistent with the dimensions of the SUR model. If X=NULL the user prefers to randomly draw this matrix using the function dgp_spsur(). This is the purpose of the argument pdfX. By default, dgp_spsur() uses a multivariate standard Normal distribution to draw the observations of the regressors; the alternative is a Uniform distribution, in the interval \(\left[0,1\right]\), for each regressor. In both cases, the regressors are generated independently.

Finally, the argument dfU informs which multivariate probability distribution function should be used to draw the values of the error terms. dgp_spsur() offers two possibilities, the multivariate Normal, which is the default, and the log-Normal distribution which means to exponentiate the sample drawn from the multivariate Normal distribution. In both cases, the covariance matrix matches the \({\bf \Sigma}\) matrix uploaded previously. The log-Normal should be taken as a clear departure from the assumption of normality.

The output of dgp_spsur() is a vector, called \({\bf Y}\), of order (TmNG\(\times\)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 \({\bf XX}\), of order (TmNG\(\times \Sigma_{g=1}^{G}p_g)\), with the values generated for the regressors of the SUR.

The example below shows the random generation of 3 variables, with a SUR structure, for 500 individuals and 1 cross-section, using a SUR-SDM model. Each equation contains 3 regressors, as specified in p. The weighting matrix, W, is obtained randomly applying the 5 nearest neighbors criteria, after having drawn the lattitude and longitude of each point in space

Tm <- 2; G <- 3; N <- 500
Sigma <- matrix(0.3, ncol = G, nrow = G)
diag(Sigma) <- 1
Betas <- c(1,2,3,1,-1,0.5,1,-0.5,2)
Thetas <- c(1,-1,0.5,-0.5,1,0)
lambda <- 0.5
co <- cbind(runif(N,0,1),runif(N,0,1))
W <- spdep::nb2mat(spdep::knn2nb(spdep::knearneigh(co, k = 5,longlat = FALSE)))
DGP <- dgp_spsur(Sigma = Sigma, Betas = Betas, Thetas = Thetas, 
                 lambda = lambda, Tm = Tm, G = G, N = N, p = 3, W = W)

As said, the output of this function is a list with the data (explained variable and regressors) in matrix form. Next code estimates a SDM model for the simulated data:

## Not included the results of execution in vignette to save time.
sdm_sim <- spsurml(Y = DGP$Y, X = DGP$X, type = "sdm", G = G, 
                   Tm = Tm, N = N, p = 3, W = W, 
                   control = list(tol = 0.01, maxit = 200, 
                                  trace = FALSE))
summary(sdm_sim)

8 Conclusions and Directions for Future Research

This article introduces a new R package, spsur, directed at estimating and solving inference on spatial SUR models. Our intention was to fill a gap in the applied spatial econometrics literature, in the sense that there are few papers which use spatial SUR models as a useful research technique. Our impression is that the lack of an efficient, friendly software has slowed its adoption in this field

The package consist of 10 inter-related functions. lmtestspsur() test for the existence of spatial effects in a model specified without spatial mechanisms, or SUR-SIM model. Then, assuming that some of these tests are positive, three functions spsurml(), spsur3sls() and spsurtime() allow to estimate, by ML or 3SLS, the prefered spatial specification. Four additional functions, wald_betas(), wald_deltas(), lr_betas_spsur() and lrtestspsur() can be used to improve the specification. The function impacts() is useful to obtain the spatial multipliers of the variables in the model selected, that come with a collection of measures of statistical significance. Finally, the function dgp_sur() is intented to be used in Monte Carlo experiments, involving spatial SUR models, or to illustrate specific features of this type of data.

Of course, this is only the first, basic version of the package that should be improved in the future. A list of actions already in the desk of the coauthors are (i) the extension of 3SLS to the treatment of endogenous regressors, others than the spatial lag of the explained variable, and the use of external instruments, (ii) a more efficient treatment of the unobserved effects, which allowed for their estimation and, at the same time, enlarge the options to demean the data, (iii) the introduction of time dynamics in a strictly stationary framework.