A two-step approach to account for unobserved spatial heterogeneity

Empirical analysis in economics often faces the difficulty that the data is correlated and heterogeneous in some unknown form. Spatial parametric approaches have been widely used to account for dependence structures, but the problem of directly deal with spatially varying parameters has been largely unexplored. The problem can be serious in all those cases in which we have no prior information justified by the economic theory. In this paper we propose an algorithm-based procedure which is able to endogenously identify structural breaks in space. The proposed algorithm is illustrated by using two well known house price data sets.


Introduction and literature review
It is a well-established fact that if spatial models are correctly specified then they can also be consistently and efficiently estimated by the commonly used estimators. However, incorrect functional forms, correlated omitted variables, models with near unit roots and row-normalized weighting matrices, and so on, typically produce spurious spatial autocorrelations (Fingleton, 1999;McMillen, 2003, Lauridsen andKosfeld, 2006;Lee and Yu, 2009;Lauridsen and Kosfeld, 2011), which can lead to inconsistency of the usual estimators. Spatial heterogeneity is a particular form of heterogeneity, usually unobserved, that is related to geo-referred data sets and would lead to misspecification of the model if not account for.
Empirical analysis in economics often faces the difficulty that the data are correlated and heterogeneous in some unknown form. A first attempt to explicitly model discontinuities in space is for example the work by McDonald and Owen (1986), which procedure has been then used by McMillen (1994) to study potential discontinuities in the population density of Chicago in 1980. As Anselin (1988a, p. 119) stressed, there are two distinct aspects that pertain to spatial heterogeneity: the former is the structural instability as expressed by changing functional forms or varying parameters, the latter is the heteroscedasticity which follows from missing variables or other forms of misspecification that lead to error terms with non-constant variance. In this paper, we are going to deal with the idea that coefficient estimates can vary over space leading to a spatial structural instability, i.e. when its parameters take on distinct values in subsets of the spatial sample.
Moreover, if spatial heterogeneity can be categorized into a small number of regimes, each represented by different values for the regression coefficients, the phenomenon is also known as spatial regimes. In this case, if spatial heterogeneity is present the functional form of the model will be misspecified because of the wrongly assumed constant relationships between dependent variables and regressors. The following Anselin's (2010, p. 5) statement is useful to understand "Spatial heterogeneity becomes particularly challenging since it is often difficult to separate from spatial dependence. This is known in the literature as the inverse problem. It is also related to the impossible distinction between true and apparent contagion. The essence of the problem is that cross-sectional data, while allowing the identification of clusters and patterns, do not provide sufficient information to identify the processes that led to the patterns." The problem of spatial heterogeneity in terms of spatially varying parameters has been largely unexplored by spatial econometricians, typically because of the main purpose of controlling only for spatial spillover effects. As Postiglione et al. (2013, p. 171) stressed "the problem of spatial heterogeneity is often neglected in empirical analysis of geographic data and this negligence can affect sensibly model estimates". Some authors have attempted to detect the presence of spatial heterogeneity by constructing statistical tests that are typically based on the LM statistic (Anselin, 1988b;de Graaff et al., 2001;Lauridsen and Kosfeld, 2006;Lauridsen and Kosfeld, 2011;Pede et al., 2014;among others). Unfortunately, once detected the presence of spatial heterogeneity no test is able to suggest how to correctly model our spatial data set and in which direction we have to proceed for further analyses. Recently, Ibragimov and Müller (2010) have derived the small and large sample properties of the t statistic, also in the context of spatially correlated data, by assuming a reasonable partition in q groups of the data. As Ibragimov and Müller (2010, p. 454) emphasized "some a priori knowledge about the correlation structure is required …". However, in practical cases, there is usually no reason to accept one partition instead of another, which is in some way justified by the economic theory.
Following a parametric approach, the typical starting point to estimate a spatial econometric model is usually based on the choice of a row-standardized spatial weighting matrix, , which is able to specify the relationship between neighboring observations. In some cases the significance of the spatial spillover effects through the autoregressive coefficient might be simply due to an omitted spatially-correlated regressors problem, which can easily justifies the use of the well-known more flexible spatial Durbin models (see e.g. Corrado andFingleton, 2011, LeSage, 2014, for comprehensive discussions). However, neighborhood influence is not calibrated in terms of the data but is prescribed by the specification of W 2 . Imposing a predefined spatial structure of the data can be sometimes too restrictive in practical cases and it can bias results when inappropriate 3 , so that McMillen (2012), among others, has 2 See the papers of LeSage and Pace (2014), Getis (2009Getis ( ,2007 for considerations on the spatial weight matrix and the autocorrelation coefficient. 3 Recently, both from a theoretical and a computational perspective, some excellent works on the definition of the W matrix has been proposed (see LeSage and Pace, 2007;Seya et al., 2013;Bhattacharjee and Jensen-Butler, 2013;Qu and Lee, 2015). In particular, Qu and Lee (2015) defined a particular endogenous W matrix (where the usual exogenous W matrix can be considered a particular case) and showed the consequences on the estimates by considering commonly used estimators in SAR cross-sectional models when the true W is endogenous. criticized this approach. Although our purpose is not to criticize the parametric approach, it is reasonable to assume that for some economic phenomena there is no reason, justified by the economic theory, to choose a priori a particular spatial structure of the underlined spatial process.
The main purpose of the present paper is then to propose a possible partition of the spatial data, i.e. a classification of the data due to unobserved heterogeneity, with no a priori information of the true dependence structure.
For instance, the advanced recent literature in hedonic house price models accounts for spatial spillover effects but still ignore the possibility of a spatial heterogeneity effect (Holly et al., 2010;Holly et al., 2011). Some researchers are then recognizing that the spatial structures can be sufficiently different that the data should not be pooled and estimated together and global spatial regression models usually fail in taking into account any potential variations over space, with the consequences of biased resulting estimates. Along the same line, one of the purposes of the present paper is to show how the presence of spatial heterogeneity might modify (and generally reduces) the significance of the spatial effects in a global spatial autoregressive model by simply estimating spatial autoregressive models with spatial regimes (which introduce a sort of flexibility in spatial autoregressive models). This can be interpreted in the fact that the spatial heterogeneity might generates part of the spatial autocorrelation effect, or in other words that the autoregressive coefficient is sometimes overestimated. In line with this view is the work of Basile et al. (2014). In order to simultaneously account for spatial dependence, unknown functional form and unobserved heterogeneity, they proposed the use of the so-called Spatial Autoregressive Semiparametric Geoadditive Models, which are based on a combination between spatial parametric autoregressive models and unknown smooth functions.
Semiparametric and nonparametric estimation methods, such as geographically weighted regressions (GWRs) (Brunsdon et al., 1996), are proving to be valid alternatives to parametric approaches and should be used as diagnostic tools to detect the presence of spatial heterogeneity 4 .
The above methods allow coefficient estimates to vary over space by calibrating the global model separately for each spatial units in order to produce n sets of parameter estimates, with the behind basic idea that simple econometric models represent the data best in small geographic areas. They indeed have the attracting advantage to control for misspecified spatial effects while using highly flexible functional forms, with the only condition that nearby observations need more weight when constructing an estimate for a target point. McMillen and Redfearn (2010) showed that the GWR specification can be viewed as a special case of the already known locally weighted regression (LWR) method 5 . GWRs or Geoadditive Models usually interpret the spatially varying parameter problem as a smooth changing over space. As a matter of fact, pairs of beta coefficient GWR estimates that are proximal in space could not exhibit statistically significant differences. As Anselin (2010, p. 17) underlined "…in models of spatial heterogeneity, the spatial regimes or spatially varying coefficients show evidence of the heterogeneity, but do not explain it. Ideally, one would want to make the structure of dependence and/or the structure of heterogeneity endogenous". This is the aim of the present paper. Our purpose is to iteratively find, in a way that spatial parameter variations can be described by breaks in continuity over space (or spatial regimes).
In the spatial statistics literature an adaptive weights smoothing (AWS) algorithm have been recently proposed by Spokoiny (2000, 2006) to describe, in a data-driven iterative way, a maximal possible local neighborhood of every point in space in which the local parametric assumption is justified by the data. The basic assumption of the proposed approach is that for every point there exists a vicinity of in which the underlying model can be well approximated by a parametric model with the constant set of parameters. Their method has been applied in the field of image analysis and it is based on both a successive increase of the local neighborhoods around every point and on a description of the local models by assigning weights to every spatial unit that 4 A valid alternative to the GWR approach is the NCSTAR model as pointed out by Lebreton (2005). 5 The main difference pertains the way of thinking the distance measure. In GWRs, distance is thought as a mere geographic distance (Euclidean, great circle, etc.), whereas in LWRs an economic interpretation can be also assigned. depends on the result of the previous step of the procedure. The potential of this method in an econometric environment concerns the possibility of endogenously obtaining cluster of observations due to unobserved spatial heterogeneity (i.e. spatial regimes) exhibiting similar coefficient estimates. Andreano et al. (2016) defined a first algorithm based on the work by Spokoiny (2000, 2006) for the identification of economic convergence clubs. In this paper we substantially modify the previous contribution proposing a two-step procedure (see Section 3 for further details and comparisons), which is based on the conjunction between the LWR approach and the AWS procedure. A two-step procedure to deal with both spatial dependence and spatial heterogeneity within the estimation of hedonic house price functions has been proposed for instance by Beron et al. (2004). However, their method focused on the estimation of a first set of parameters related to environmental characteristics in the first step, whereas it proposed a spatial econometric model that accounts for both spatial dependence and spatial heterogeneity in the second step, assuming a spatial trend as a quadratic form of latitudes and longitudes for the heterogeneity effects. In our two step procedure, instead, we propose a first step that focuses on the estimation of unobserved discrete spatial heterogeneity (i.e. spatial regimes) and a second step that estimates the effects of both spatial dependence and the identified spatial regimes. Our purpose is based on the idea that we can combine the potential of local estimation with the usefulness of a modified AWS procedure, which is able to identify spatial regimes, i.e. subsamples over space with an estimated set of beta coefficients for each of them.
The paper is structured in the following way. In section 2 we introduce the LWR and the GWR as a special case. In section 3 we explain the first step of our procedure, i.e. the algorithm, whereas section 4 explain the second step to estimate both spatial dependence and spatial regimes. Section 5 illustrates the data set used and their main estimation results in terms of the marginal effects obtained by using different spatial econometric models. Finally, section 6 concludes.

Locally Weighted Regressions
Spatial econometric models may be not appropriate in the presence of (unobserved) spatial heterogeneity. Locally weighted regressions (LWRs) (Cleveland and Devlin, 1988) or Geographically weighted regressions (GWRs) (Brunsdon et al., 1996;Fotheringham et al., 1998, Fotheringham et al., 2002, which are recognized to be natural evolutions of the expansion method (Casetti, 1972), allow us to estimate local rather than global parameters.
Residual terms usually exhibit a different from zero spatial autocorrelation parameter that, actually, might be not statistically different from zero if the true reason of error autocorrelation is different from a true contagion process (i.e. spurious autocorrelation). The first law of geography (Tobler, 1970, p. 236) states: "everything is related to everything else, but near things are more related than distant things". The goal should be to detect the reason why the closest things are related. The basic idea behind local geographical estimation is that simple linear functions may fit well for observations close to a site, say i, but probably they will be inappropriate when more distant observations are included. Limiting the estimation to a neighborhood of observation i eliminates much of the heteroscedasticity and autocorrelation that is endemic to spatial data set. Therefore, simple linear functions can be written to account for local parameter estimates in the following way where is an -dimensional column vector, is an by + 1 matrix with i-th row = ( 0 , 1 , 2 , … , )′, ⊙ is the Hadamard product operator in which each element of is multiplied by the corresponding element of the i-th row in , = ( 0 , 1 , 2 , … , ), is a ( + 1) column vector of ones, and is an -dimensional column vector of innovations. The idea is that, in order to estimate a set of parameters for a given point i, one can approximate the model in [2.1] in the neighborhood of i by using a simple linear model, and perform least squares with a subset of points that are closed to i. The objective function, for each local regression i, takes the following where is an n-dimensional diagonal square matrix whose diagonal elements, ( 1 , 2 , … , , … , ), denote the weights (in our case geographic distances) of each of the n observed data for regression point i. We will therefore have n diagonal spatial weighting matrices and n sets of local parameter estimates that correspond to the local marginal effects. Then, a locally weighted least squares (LWLS) estimator is simply obtained by repeated weighted least squares Initial Kernel weighting functions have to be specified to define the matrices. Such functions are typically the exponential, the Gaussian, the bi-square, the tri-cube, and they seem to be similar in placing higher weight on nearby observations and lower weight on the distant ones, but they differ according to the type of bandwidth (b).
The choice of the bandwidth value is then crucial: the neighborhood of observation i is defined by the value of b which then determines the number of observations (i.e. a subsample) that receive a weight in constructing the estimate for i, and how rapidly these weights decline as distance increases. So, an important question could be how close to i should points be considered, or in other words in which way we can define a proper neighborhood for i. The observations are used several times with a maximum equals to n 2 . Every local ̂ can be based on the same set of n observations, but in this case the local estimates would be the same. The minimum number of observations to be used for each local estimation is instead + 1, but then the ̂' s would be not significant. The point is then to choose an ideal subsample of observations for each local estimation, i.e. b should be somewhere in between n and + 1 such that we have sufficient observations for each ̂ and part of them will be significant. This problem is referred as a bias-variance trade-off in choosing the value of b (Brunsdon et al., 1996). Since the LWLS estimate for point i, ̂, is only an approximation of the true value in the same point, one can consider observations very closed to i with the justification that it is reasonable to assume a lower magnitude of difference in the beta values. However, because of the sample size reduction, the standard error of ̂ will increase (i.e. ̂ are not significant). The opposite will occur if we increase the sample size. To solve the problem of choosing an optimal bandwidth value, , a cross-validation (CV) method has been widely used, which minimizes the overall residual sum of squares obtained when observation i is deleted (i.e. the target point). A different and more general approach based on a trade-off between goodness of fit and degrees of freedom is the AIC 6 (Fotheringham et al., 2002, p. 61) = min = min 2 ln(̂) + ln (2 )  Another important issue is that the bandwidth value can be defined as fixed or variable/adaptive.
Since that adaptive bandwidths suit highly irregular sample configurations and they ensure sufficient (and constant) local information for each local calibration, the adoption of an adaptive bandwidth seems to be preferred to a fixed bandwidth for analyzing spatial data. The choice of an adaptive bandwidth is then useful for two reasons: (i) it guarantees the same number of observations (i.e. the same amount of information) for each local estimation, (ii) it can potentially reduce the edge effects in spatial data analysis. Therefore, in this paper we calculate an adaptive bandwidth following the nearest neighbor approach and we define, for each used kernel weighting function, the optimal variable bandwidth or window size, , by using the minimizing the AIC criterion in [2.4]. The bandwidth value (b) selects the b-smallest distances for each initial local estimation. Each of the elements of the geographical distances are then normalized through the maximum distance, so that it does not matter how far is the most distant point within the optimal neighborhood defined by b for each local estimation. In this way, through an appropriate kernel function, the weights that contribute to each local estimation do not depend on the effective position of the points in space. In the following two paragraphs we intend to explain our two-step procedure. The first step (i.e. Section 3) is an algorithm that is able to endogenously identify spatial regimes, whereas the second step (i.e. Section 4) is based on the specification and estimation of two spatial autoregressive models with the identified spatial regimes.

The first step: an algorithm to deal with unobserved spatial heterogeneity
In this section we intend to explain our iterative procedure, which enters as a first step of the entire estimation approach. The idea comes from the adaptive weights smoothing (AWS) procedure first proposed by Polzehl and Spokoiny (2000) is the variance-covariance matrix of ̂ and Σ j is the same matrix for ̂. This variance-covariance matrix is different from that in Andreano et al. (2016) in order to preserve degrees of freedom.
Each vector ̂ for = 1, … , , at iteration , is obtained from the LWLS problem in [2.3] by only using a subsample of observations whose size is defined by a mix of one of the possible kernel functions and a bandwidth value. If the number of observations is reasonably large in all defined subsamples (which contribute to different beta vector estimations), the approximate multivariate normality is a standard result ̂~( , Σ i ) (see Ibragimov andMüller, 2010 andBester et al., 2011 is the − of the chi-square distribution (for details see Polzehl and Spokoiny, 2006 where ∈ (0,1) is a control parameter (also called memory parameter in Polzehl and Spokoiny, 2006) which is used to stabilize the convergence procedure. We set = 0.5, so that we consider the exact arithmetic mean between the weights at the ( − 1)-th and the -th iterations. In this way, the new weights are smoothly or gradually introduced to avoid potential convergence problems due to discontinuities in the iterative process. During the convergence procedure, it is important to provide enough observations for each formed group. Because of the identification problems, groups characterized by a limited number of observations (i.e., less than the number of parameters to be estimated) will be automatically excluded from the analysis and treat as outliers. A final useful definition for the description of the convergence procedure is the weights variation function which tends to zero as the convergence procedure stabilizes (see Figure 3). This procedure does not impose any restrictions and it is fully adaptive in the sense that no prior information about the spatial structure is required. To summarize, the convergence procedure is based on the following steps  Define initial weights with a kernel that is a function of distances between pairs of units in space and an optimal bandwidth value (the econometric model must be assumed to be the true one) and use these weights in equation [

The second step: spatial autoregressive models with spatial regimes
It is already well-known that (first-order) autoregressive spatial model with (first-order) autoregressive disturbances (SARARs) models and spatial Durbin models (SDMs) provide general starting points for spatial analysis, since they subsume, in different ways, both spatial autoregressive error (SEM) and spatial autoregressive (SAR) models (LeSage and Pace, 2009). The two spatial models are defined as where is an n by p+1 matrix of regressors representing own-region characteristics with a p+1 by 1 vector of parameters , ̆ is the same matrix without the intercept and is the corresponding p vector of parameters, is the autoregressive parameter, is an n by 1 vector of autoregressive error terms with spatial coefficient and spatial weight matrix 2 , is an n by 1 vector of i.i.d. error terms. 1 is the commonly used n by n spatial weighting matrix which can be different or equal to 2 . In this paper we set 1 = 2 = and we define as a k-nearest neighbor weight matrix with = 10 and geographical distances. The squared weighting matrix has the diagonal elements wii equal to zero, the rows and columns uniformly bounded as n goes to infinity, and it is then normalized, so that the admissible parameter spaces of and are known and less than unity in absolute value.
It is worth noting that, although both the weighting matrix used in this section, i.e. , and the one neighbor criterion (i.e. based on the optimal adaptive bandwidth value) is only used to define the initial set of weights in , so that this is an issue only for the first set of local estimates within the iterative procedure. Finally, there is no need of normalization to ensure admissible parameter spaces of autoregressive coefficients in the case of 8 .
Recently it has been emerged that SDMs are useful model specifications to overcome the omitted variable problem especially with hedonic house price data (Gerkman, 2012). In order to show how estimating models with spatially varying parameters can produce substantial differences in terms of the autocorrelation effects, in this paper we estimated different global models (GMs): OLS, SAR, SEM, SARAR, and SDM. Those models are then compared with the same models adding spatially 8 The weights in the diagonal of have the role, at first, of selecting the observations that contribute to local estimations trough local weighted least squares estimations, whereas the weights in permits, trough the reduced form of the model, the expansion of the infinite number of cross-sectional effects in the entire system (i.e. global spillover effects).
varying coefficients, after using the ILWR approach to identify the number of clusters in space. We called these model specifications Endogenous Spatial Regime models (ESRMs where, in both cases, ∼ (0, 2 ), = 1, … , , and is the unknown number of clusters in space. After having defined these unknown clusters through the first step procedure, we proceed with the estimation of models in [3.7] and [3.8], which can be interpreted as spatial autoregressive models with structural breaks in space (also known in the literature as spatial regimes).

Data sets and results
In order to detect the presence of unobserved spatial heterogeneity in the form of spatial regimes and to prove the usefulness of the above two-step approach, we use baltimore and house data sets 9 which summarizes the information on house sales prices in Baltimore, MD, and Lucas County, Ohio, respectively (see Table 1). Baltimore data set was first analyzed by Dubin (1992) who emphasized the problem of spatially varying parameters. Because of the computational burden of our algorithm in the first step we selected the central subsample of the cross-sectional house data set in 1993 (see Figure 1).
INSERT with subsamples 1 = 101 and 2 = 110, mostly dividing the central and the northwest areas from the south and the northeast ones (Figure 2(a)). Considering AIC values, ESRMs are always to be preferred over GMs and in both situations the SAR specification seems to fit the data best. Then, in this case positive spatial neighbor effects seem to persist despite the presence of two spatial regimes. However, if we see the results of SDM with spatial regimes (SDM-SRs) we can confirm the hypothesis that most of the spatial autocorrelation is induced by the omitted variable problem of not considering the lagged covariates. In this case the presence of spatial regimes has a "pure" informative role: the SDM-SR specification (with lowest AIC) is more useful than the others since it can be used to guide zone-specific policies. INSERT TABLE 2 HERE   INSERT TABLE 3 HERE   Tables 4 and 5 show GMs and ESRMs estimates for house data set, instead. The bandwidth value corresponds to 19 nearest neighbors, = 19, over 382 number of observations. The first step identified 4 homogeneous clusters ( = 4) with subsamples 1 = 57, 2 = 141, 3 = 58, 4 = 126 (Figure 2(b)). Conclusions are a bit different from those obtained for Baltimore: not only ESRMs provide better results in terms of AIC values, but also spatial spillover effects seem to be attenuated by the presence of the spatial regimes. This is particularly clear in the comparison between a global SDM and an SDM-SRs where the autocorrelation coefficient goes down from 0.655 to 0.248. In this case the additional information of the identified spatial regimes makes us able to split the apparent contagion (in Anselin's words) into the true contagion (the spillover effects), the part due to the lagged covariates, and the part due to the presence of spatial heterogeneity.

INSERT FIGURE 2 HERE
INSERT It is then interesting to note that estimating spatial models with spatial regimes (endogenously identified with the first step of our procedure) can be not only useful to account for both spatial dependence and spatial heterogeneity giving more flexibility to our models, but also to provide zone-specific economic results. Considering for instance Baltimore results (Tables 2-3), we can observe that some regressors in ESRMs are highly statistically significant in the first zone (i.e., central and northwest areas) whereas not in the second (i.e., the south and the northeast areas).
Therefore, for the same model a significant global parameter estimate should be mainly due to the statistical significance of the same parameter in central and northwest areas.
It is well known in the literature that house prices do not necessarily change at a uniform rate throughout an urban area. Several papers that attempt to smoothly describe a spatial price variation or to identify housing submarkets within cities have been published (see e.g. McMillen, 2013, for a recent exhaustive review). The identification of regimes in housing markets allows us to identify a set of discrete neighborhoods with constant coefficients within each subarea and discrete changes across boundaries. Therefore, this procedure could be viewed as a solution for the identification of the potential number of submarkets (see e.g. Watkins, 2001) within an urban area (e.g. a city) especially when no prior knowledge on the geographical structure is available.

INSERT TABLE 6 HERE
INSERT FIGURE 3 HERE Table 6 shows the number of iterations of the algorithm and the time required to ensure the convergence (see Figure 3). Although the computational time seems to be still reasonable as long as the number of observations are lower than 400, we can note that the time is dramatically increased by an augmentation of only around 100 of observations, and the number of iterations are twice as much those required for baltimore data set. Conversely, the number of explanatory variables seems not to change the time required by our algorithm. Finally, in Figure 3 we can see how the convergence procedure works: after an initial period of "instability", the weights tend to zero as the number of iterations increased. The most significant and rapid decrease is between the interval of iterations (5 th -15 th ) for baltimore and (15 th -30 th ) for house.

Conclusions and future developments
Spatial heterogeneity is one form of unobserved heterogeneity, which leads to model misspecification problems. Geographically weighted regressions or geoadditive models are surely promising alternatives to account for spatially varying parameters, in which the parameter surface is not assumed to be constant over space. However, economic studies usually suggest a form of spatially varying parameters, which can be interpreted as the presence of discrete spatial heterogeneity (or spatial regimes). In this paper we then define a two-step approach able to identify these regimes in a first step and to simultaneously estimate spatial dependence and heterogeneity in a second step.
By using two house price data sets, we found that our procedure was able to identify clusters and that spatial econometric models estimated on that partition were to be preferred. Most of the results showed that spatial neighbor effects are present in our data sets but in a reduced form, and probably part of the neighborhood effects is not only due to an omitted variable problem but also to the unobserved spatial heterogeneity. This firstly suggested that accounting for both spatial dependence and some form of spatial heterogeneity (e.g. the spatial regimes, with the use of some algorithmbased approach to identify them) is a fundamental step. The additional information derived from the identified spatial regimes makes us able to split the apparent contagion into three parts: the true contagion (the spillover effects), the lagged covariates (which justifies the use of the spatial Durbin model), and the presence of spatial heterogeneity. Moreover, zone-specific economic results can be provided with these more flexible spatial models.
As expected, its main drawback is the computational time required for the convergence procedure, which we suppose can become infeasible with very large data sets. The need to overcome the ILWR computational problems is straightforward. Unobserved spatial heterogeneity is still a problem in spatial econometrics. Many author are still trying to test the presence of spatial heterogeneity, but the prospect of considering algorithm-based methods is becoming a successful way to follow in order to solve practical problems with several data sets and to consider higher flexibility in our models.
The first step of our two-step procedure does not impose any restrictions and is fully adaptive in the sense that no prior information about the spatial structure is required. It can be easily extended to any other model which can be locally estimated by weighted least squares, and with different definition of the weights (e.g. economic definitions rather than mere geographic distances).
Moreover, two possible and interesting modifications of this algorithm can be the one that permits the simultaneous consideration of both spatial dependence and local marginal effect estimation inside the iterative procedure and the extension to panel data sets to also deal with time variations.
These will be accounted for in future works.    Table 3 Parameter estimates of models with spatial regimes in the baltimore data set based on the partition in Figure 2(a).  Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Table 5 Parameter estimates of models with spatial regimes in the house data set based on the partition in Figure 2(b).