Next Article in Journal
Symmetric Mass Generation
Next Article in Special Issue
Reduction in Waiting Time in an M/M/1/N Encouraged Arrival Queue with Feedback, Balking and Maintaining of Reneged Customers
Previous Article in Journal
Joint-Prior-Based Uneven Illumination Image Enhancement for Surface Defect Detection
Previous Article in Special Issue
Generating Optimal Discrete Analogue of the Generalized Pareto Distribution under Bayesian Inference with Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Estimation of Asymmetric Spatial Autoregressive Dependence on Irregular Lattices

1
Faculty of Civil Engineering and Geodetic Science, Leibniz University Hannover, 30167 Hannover, Germany
2
Chairs of Statistics and Econometrics, University of Goettingen, Humboldtallee 3, 37073 Göttingen, Germany
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Symmetry 2022, 14(7), 1474; https://doi.org/10.3390/sym14071474
Submission received: 10 June 2022 / Revised: 8 July 2022 / Accepted: 15 July 2022 / Published: 19 July 2022

Abstract

:
In spatial econometrics, we usually assume that the spatial dependence structure is known and that all information about it is contained in a spatial weights matrix W . However, in practice, the structure of W is unknown a priori and difficult to obtain, especially for asymmetric dependence. In this paper, we propose a data-driven method to obtain W , whether it is symmetric or asymmetric. This is achieved by calculating the area overlap of the adjacent regions/districts with a given shape (a pizza-like shape, in our case). With W determined in this way, we estimate the potentially asymmetric spatial autoregressive dependence on irregular lattices. We verify our method using Monte Carlo simulations for finite samples and compare it with classical approaches such as Queen’s contiguity matrices and inverse-distance weighting matrices. Finally, our method is applied to model the evolution of sales prices for building land in Brandenburg, Germany. We show that the price evolution and its spatial dependence are mainly driven by the orientation towards Berlin.

1. Introduction

Geospatial analysis is based on Tobler’s first law of geography, which points out that everything is connected to everything else, but that nearby objects dominate (Tobler [1]). Not every process can be described by applying this rule, and there is no precise and unique definition of “nearby”. Therefore, it is often assumed that the dependence structure is known through underlying physical systems (e.g., river flows) or geographical information or network structures (e.g., public transport connections), or assumptions such as symmetry are made. Generally, to model spatial data, one needs to know the n ( n 1 ) potential interactions among the system’s n objects. A major challenge is to obtain these interactions from the n data points in the sample.
There are various attempts to model such spatial dependence. For instance, the spatial covariance could directly be described by a parametric covariance function that has fewer parameters than all possible links. This approach is commonly known as the geostatistical approach (see, e.g., Cressie and Wikle [2]). Alternatively, spatial dependence could be modeled by explicitly including spatially lagged variables. These models are often referred to as spatial econometrics models. Processes in which the objects/regions influence each other can be modeled with spatial autoregressive (SAR) models (LeSage [3], LeSage and Pace [4]). SAR models describe processes in which the observed value of one region influences the observed value in other regions and vice versa. An exemplary question that such models attempt to answer is whether the average salary in one region influences the average salary in adjacent regions, and if so, to what extent. Such interplay among n regions is described by an n × n spatial weights matrix W .
These autoregressive-type approaches involving the specification of a weights matrix have been widely used in different areas. For instance, in environmental statistics, the effects of weather conditions on fertilizer application were modelled (Billé and Rogna [5]) and environmental expenditure interactions among OECD countries were investigated. In addition, the impact of COVID on financial returns was investigated (Billé and Caporin [6]), and studies of the labour market have employed spatial weight matrices (Billé [7]). Furthermore, in health economics, Donegan et al. [8] used spatial econometrics approaches for modelling community health.
There are different methods for obtaining a suitable W . One is to assume that all regions sharing a common border influence each other. This method is also called the Queen’s contiguity matrix, alluding to the chessboard (LeSage and Pace [4]). If region A and region B share a common border, the Queen’s contiguity matrix is then constructed with W A B = W B A = 1 . Other methods involve geographical distances such as those used in Lin Lawell [9] to model air pollution. They assume that the first-neighbour regions are the regions within 500 km of each other. Another method for constructing W is to assume that the regions’ influences on each other declines as the distance between them increases. The so-called ‘inverse-distance W ’ was used by Boly et al. [10] to model real-estate valuation and by Zhao et al. [11] to describe air pollution in China. The inverse-distance matrix is constructed with W A B = W B A = 1 / d i s t a n c e ( A , B ) . An advanced method also takes geographical or economical connections into account when constructing W . Krisztin et al. [12] described the worldwide spread of the coronavirus by assuming that countries with a common border influence each other, as well as countries connected by an airline. Other suggested methods for constructing W include the parametric or semi-parametric estimation method proposed by Pinkse et al. [13], the method of Stakhovych and Bijmolt [14], where W was selected from a set of possible candidates with a goodness-of-fit criterion, requiring that the true W is in the set of candidates, and that of Bhattacharjee and Jensen-Butler [15], where W was selected from an estimated spatial autocovariance matrix under the conditions of symmetry and a finite sample size of n.
Ahrens and Bhattacharjee [16] proposed to estimate W in two steps with a least absolute shrinkage and selection operator (lasso) approach. An alternative penalized estimation approach was proposed by Lam and Souza [17]. They selected W from a linear combination of different weights matrices by setting irrelevant components to zero. This linear combination could involve higher-order spatial lags (see, e.g., Cohen and Paul [18]). A method for selecting the best weights matrix from a set of candidates was also proposed by Debarsy and LeSage [19] or by Debarsy and Ertur [20]. Zhang and Yu [21] showed that the true weights matrix can consistently be selected from the set of candidates (if included) using Mallow’s C p criterion. If the true matrix is not in the set of candidates, the selection is still optimal, and they proposed a model averaging technique over different candidates. In addition, the method was illustrated by an analysis of historical rice prices in China.
Most of the studies have in common the fact that they apply symmetric weighting schemes (if A influences B, then B influences A, but with a potentially different magnitude). To the best of our knowledge, there are only a few references that deal with the construction of weights matrices accounting for asymmetry. Exceptions include, for instance, Zhou et al. [22], who constructed their asymmetric weights matrix by element-wise multiplication of a nearest-neighbour weights matrix (mostly symmetric) with a matrix that indicated whether the influence flowed from A to B or vice versa. The direction of the influence must be known in advance. A similar method for constructing an asymmetric weights matrix was proposed by Merk and Otto [23], to model PM 2.5 concentrations. In their approach, the directional component depended on the wind direction and wind speed, which strongly influenced the spatial dependence structure.
However, in practice, we often do not have prior knowledge about the asymmetrical or directional influences. Therefore, in this paper, we develop a data-driven method to estimate the interactions directly from the sample data and explicitly allow asymmetric dependence for irregular lattice data. We assume that the interactions are not dependent on the location of the region but are equal over the entire sample area. With this approach for irregular lattice data, we extend the work of Merk and Otto [24], who proposed a similar approach for obtaining the interactions for regular lattice data based on a lasso procedure. The proposed method respects the interactions from close neighbours, like the Queen’s contiguity matrix or the inverse-distance matrix, but allows asymmetric interactions (if A influences B, B does not necessarily influence A). It is worth noting that the proposed method provides a method of obtaining a flexible weights matrix. This weights matrix can be an additional candidate for the selection of the weights matrix described above, for which good candidates must be available.
This paper is structured as follows. Section 2 describes the theoretical model. Section 3 covers a Monte Carlo simulation study. In Section 4, we apply our method to real-world data and analyse the evolution of land sale prices in Brandenburg, Germany, to show that the short-distance directional influences are coming from Berlin. We conclude by discussing new research directions in Section 5.

2. Theoretical Model

Let { Y ( s ) R : s B } be a univariate process in the spatial domain B. We explicitly consider B to be an irregular lattice, as would be the case for spatial polygon data. Shapes like this can be found in countries, states, and districts. An artificial example can be seen in Figure 1.
Furthermore, suppose that the process is observed for a finite set of locations B = { s 1 , , s n } and Y = ( Y ( s 1 ) , , Y ( s n ) ) t . We consider that the process follows a spatial autoregressive model, i.e.,
Y = W Y + X β + ϵ ,
where W is an unknown spatial weights matrix and β is a vector of p exogenous influences. Furthermore, X is a ( p + 1 ) × n matrix and ϵ is a random error vector, which is assumed to have zero mean with a diagonal, homoscedastic covariance matrix σ ϵ 2 I . The identity matrix is denoted by I . Since W consists of n ( n 1 ) unknown parameters, which must be estimated, and only n values of the response are observed, W is typically replaced by a multiple of a pre-specified matrix W ˜ , i.e., W = ϕ W ˜ , where ϕ is an unknown scalar parameter to be estimated (Anselin [25]).
In the following, we aim to estimate the full spatial dependence structure W . To obtain a meaningful and easy-to-interpret representation, suppose that the spatial dependence can be separated into k directions and d distances (see Figure 2a). Furthermore, suppose that there is a unique weighting for each segment of the sectors such that W can be decomposed as
W = i = 1 k j = 1 d ϕ i j W ˜ i j ,
where each matrix W ˜ i j has positive weights for the ( i j ) th segment only (see Figure 2b). The weights are chosen to be proportional to the overlapping areas. That is,
W ˜ i j = 0 w i j , 12 w i j , 1 n w i j , 21 0 w i j , 2 n w i j , n 1 w i j , n 2 0 = ( w i j , ι η ) ι , η = 1 , , n
where w i j , ι η is the relative area of the η th location lying in the ( i j ) th segment with respect to region ι . For example, (see Figure 2b), one w i j , ι η is the normalized area of 2. Due to normalization, each W ˜ i j is bounded by 1 in the row sums. Each of these matrices contains the relative overlapping areas of one segment with the n regions. Since a segment usually overlaps with a small fraction of the regions only, all W ˜ i j are sparse.
Figure 2b depicts an artificial map of ten regions. An exemplary segment ( ϕ i j ) determines the spatial dependence of the region S 1 towards the north by the relative sizes of its overlap with the neighbouring regions S 2 and S 3 . As we exclude self-influences, the overlap with region S 1 is ignored, and the weights for the spatial dependence are normalized with respect to the two remaining intersections. Here, the weights are given by w i j , 12 = 0.6 and w i j , 13 = 0.4 for S 2 and S 3 , respectively. With this construction, the matrix ( w i j , ι η ) is row-normalized, as w i j , 12 and w i j , 13 are the only non-zero entries in row ι = 1 .
This leads to the classical autoregressive model with higher-order spatial lags, given by
Y = i = 1 k j = 1 d ϕ i j W ˜ i j Y + X β + ϵ .
In the reduced form, we obtain with
S { ϕ i j : i = 1 , , k ; j = 1 , , d } = I i = 1 k j = 1 d ϕ i j W ˜ i j
that
Y = I i = 1 k j = 1 d ϕ i j W ˜ i j 1 ( X β + ϵ ) = S { ϕ i j } 1 ( X β + ϵ ) .
In this paper, we consider that d and k are chosen to be large enough to obtain a flexible model reflecting the true underlying spatial dependence precisely, generally leading to a rich, parameterized description.
The inverse S 1 exists if the column-sum norm W 1 < a or the row-sum norm W < a is bounded by a finite number a. With the W ˜ i j row sums bounded by 1, this can be ensured under the assumption that
i = 1 k j = 1 d | ϕ i j | < 1 .
A more general condition is discussed by Elhorst et al. [26]. Since we only want to describe the direction from which we see a stronger or weaker influence, we constrain ourselves to the first condition.
Implementing the model demands a proper choice of the parameters k and d, which determine the division of the direction and distance of the spatial dependence, respectively. Moreover, the actual length of each distance step can be optimally chosen. We perform a spatial partitioning analysis to obtain a robust choice of all three parameters.
Let l q be the distance between the qth and ( q + 1 ) th distance steps. We suggest selecting k K N , d D N , and l q L R > 0 using Akaike information criterion (AIC) selection (Akaike [27]). Therefore, we calculate all W ˜ i j ( k , d , l q ) for different sets of D × K × L . Then, we estimate the unknown parameters β , ϕ i j with maximum-likelihood methods and select the best-performing model based on its AIC. Having obtained the best-fitting ( k ^ , d ^ , l ^ q ) , we estimate the final model in the next step.
We estimate the parameters { β , ϕ i j } using the maximum-likelihood approach combined with parameter selection using minimal AIC. For both symmetric and asymmetric dependence, a large share of the ϕ i j values can be expected to be zero, because k and d are large. For this reason, we repeat the estimation d × k times and, at every step, we drop the least significant ϕ i j parameter (i.e., we set this particular ϕ i j to zero). Finally, from these d × k estimations we choose the one with the lowest AIC.
In general, the unknown parameters { β , ϕ i j } can be estimated. For that reason, the joint probability function, the log-likelihood function f Y ( Y ) , is maximized with respect to parameters ϕ i j and β . Assuming that ϵ N ( 0 , σ ϵ 2 I ) , the log-likelihood is given by
L ( Y | β , ϕ i j ) = n 2 log ( 2 π σ ϵ 2 ) 1 2 σ ϵ 2 ϵ ˜ t ϵ ˜ + log det ( S )
with
ϵ ˜ = I i = 1 k j = 1 d ϕ i j W ˜ i j Y X β = S Y X β .
The estimator is obtained by maximizing Equation (8) with respect to all parameters. For the consistency of the resulting estimators, we refer to Lee [28] and Gupta and Robinson [29].

3. Monte Carlo Simulations

For the simulation study, we simulate the irregular spatial lattice as Voronoi cells (Longley et al. [30], Sen [31]), where n { 200 , 500 , 900 } centroids are sampled from a two-dimensional uniform distribution on the interval [ 0 , n ] 2 (see Figure 1). Eventually, Y can be simulated as in Equation (6), with ϵ being uncorrelated realizations of a standard normal distribution. We choose β = ( 3 , 1 , 2 , 0 , 5 ) t and X has only ones in the first column. The remaining elements are standard normally distributed random values.
We consider four different specifications of the spatial dependence with k = 8 and k = 16 , as shown in Figure 3 and Figure 4. First, an isotropic process is considered where only the first lag has an influence, as in classical contiguity settings. Second, a directional process is considered with a clear north-to-south dependence. Third, we consider a nearest-neighbour dependence in the northwest direction only. Finally, the fourth setting is extended by adding another level of dependence strength, while retaining the mainly northwest directional dependence. Blue sections represent zero influence, orange sections represent the maximum influence, and purple sections represent 50 % of the maximum influence. The maximum influence is obtained by setting i = 1 k j = 1 d | ϕ i j | = 0.95 .
We carried out two types of simulation. In the first type, we used d = 2 , k = 8 and l = 375 , employing the spatial dependence structures depicted in Figure 3. We first performed the spatial partitioning (results in Table 1 and Figure 5) to gain the best k ˜ and l ˜ for the estimation of ϕ ^ i j (we use k ˜ to denote the optimal parameter and k ^ to denote an estimator). The procedure to determine k ˜ and l ˜ is described below, and was repeated 10 3 times. Then, using the optimal ( k ˜ , l ˜ ) from each of the 1000 estimations, we estimated ϕ ^ i j and compared it to alternative estimations with the Queen’s contiguity matrix, the inverse-distance neighbour matrix by AIC, and the sum of squared residuals ( ϵ 2 ). The comparison can be found in Table 2, showing that our method works best.
In the second type of simulation, we show that we can re-estimate even higher-dimensional specifications of the spatial dependence, as depicted in Figure 4. We simulated Y with d = 2 , k = 16 and l 1 = l 2 = 375 using this high-resolution spatial dependence and estimated ϕ ^ i j , conditionally upon selecting the true d, k and l as the best parameters. The results are summarized in Table 3. One can see that we can re-estimate these higher-dimensional settings with a true positive rate of at least 76%. The simulations were carried out in R using the packages spdep and rgdal (Bivand et al. [32], Bivand et al. [33]).
The spatial partitioning algorithm to determine the optimal parameters ( k ˜ , l ˜ ) was implemented as follows. Initially, we select the distance d = 2 and split the spatial dependence into different lengths l q = l { 225 , 250 , , 525 } and various numbers of directions k { 4 , 5 , 6 , 7 , 8 , 9 , 10 } . For each spatial dependence structure, all d × k spatial weights matrices W ^ i j are calculated. We calculate 3 × 1000 sets of Y with k = 8 and l = 375 ( Y k = 8 , l = 375 ) that differ by a random white noise element ϵ . Finally, we re-estimate Y ^ k l for all combinations of k and l and calculate their AIC values to select the best set of ( k ˜ , l ˜ ) for each of the 1000 Y k = 8 , l = 375 . The results can be seen in Table 1 and Figure 5.
Figure 5 shows that for 8-setting-1, the symmetric setting, it is more difficult to estimate the correct k, because the models are observationally equivalent. The directional influence for this radially symmetric setting looks almost identical for k = 4 and k = 8 . Since settings with fewer variables benefit from AIC selection, observation of this pattern is to be expected. The length of the slices (l; here the y-axis) can be identified.
However, 8-setting-2 is axially symmetric about the y-axis. Only the northern half of the sections can provide spatial influence. In this case, the spatial partitioning algorithm is less robust in determining the length parameter l ˜ . One plausible reason is the ambiguity between the length and the number of distance steps, which occurs particularly for this spatial dependence setting. For 8-setting-3, the true values of k and l are selected consistently. For 8-setting-4, estimation is more difficult, but the results show reasonable consistency for larger values of n. In general, k ˜ was rarely selected to be too large, but the selection of l ˜ was more difficult.
For each of the previous 10 3 estimations, we obtained the best values for k ˜ and l ˜ . Now, using these k ˜ and l ˜ values, we estimate, 10 3 times, the values for the segments ϕ ^ i j , as described in Section 2. Finally, we compare the mean AIC values and the sum of squared residuals ( ϵ 2 ) of those estimates to the commonly used Queen’s contiguity matrices (Lin Lawell [9]) and inverse-distance weighting matrices (Boly et al. [10]).
The average AIC value and the sum of squared residuals are reported in Table 2. For all settings and all n, our method outperforms the common procedure, i.e., smaller values for AIC and ϵ 2 are obtained.
For the final estimation, the spatial dependence is split into k = 16 directions and d = 2 distances (see Figure 4), leading to d × k = 32 spatial weight matrices W ˜ i j with dimension n × n . We assume knowledge of the true parameters of ( k ˜ , d ˜ , l ˜ p ) , and we re-estimate the parameters of ϕ ^ i j . Table 3 shows the results of the simulation study for the different dependence structures (see Figure 4), with 500 estimations for each n and for each setting. Since we know the true values, we compare the BIAS and the RMSE of the estimations. For the BIAS and RMSE results, we distinguish two cases. Either the true values should be chosen to be zero ( = ! 0 ), or the true values should be chosen to be greater than zero ( > ! 0 ).
The true positive values represent the rate at which ϕ ^ i j was correctly identified as positive. The false non-zero value represents the rate of falsely identifying zero values as positive. From Table 3, we can see that the results improve when a larger n is chosen. The minimal true positive rate was 76%, but, in most cases, the true positive rate exceeds 90%. Furthermore, we chose up to only 7% of the values as falsely non-zero. From the true positive and false non-zero values, one can see that, for a very simple dependence structure such as 16-setting-3, the sections ϕ ^ i j are identified almost perfectly.

4. Real-World Application

In this section, we apply our approach to the evolution of land prices in Brandenburg, Germany. In general, spatial autoregressive or spatial lag models are widely applied to model housing prices or the evolution of house prices (e.g., Fingleton [34], Osland [35], Baltagi et al. [36], Baltagi et al. [37], Jin and Lee [38]). Furthermore, Taşpınar et al. [39] showed that there are spatial interactions in the volatility of house-price returns. Below, we focus on the spatial dependence structure due to the specifics of the studied area. The state of Brandenburg, with 2370 districts, completely surrounds Berlin. Berlin does not belong to Brandenburg and is therefore not included in the dataset. It is to be expected that the distance to Berlin has a strong influence on the evolution of Brandenburg’s sales prices. It is also to be expected that the influences the districts have on each other is stronger in the vicinity of Berlin.
We consider the sales prices of land that was sold as building land for the first time. The dataset contains the sales prices from 2005 to 2014 for the 2370 districts of Brandenburg. To model the evolution of the sales prices, we take the average prices per square metre for all regions between 2005 and 2009 and compute the difference from the average price for the regions from 2010 to 2014. For regions that do not have sufficient data, we replace the difference in selling price with zero, that is, we assume no price changes between these periods. The resulting map is plotted in Figure 6, on the right.
We wish to investigate whether the districts in the dataset exert a short-distance influence on each other, whether this influence is directional, and in which direction the influence flows. Therefore, we choose the setting such that our weights ϕ i j are oriented towards Berlin (similarly to a compass pointing north, hence in Figure 6 the green pieces are pointing to Berlin). Furthermore, we include the following regressors: (1) the distance to Berlin for each district, (2) the squared distance to Berlin, (3) the inverse distance to Berlin, and (4) an intercept of 1. To find the optimal number of directions k and the best length l, we calculate the K × L W ˜ i j for a subset of n 0 = 58 of the n = 2370 districts to lower the required computational effort (see Figure 6, left).
The results of the spatial partitioning are plotted in Figure 7 for distances d = 1 and d = 2 . In both cases, the best number of directions is k ˜ = 3 . In the case where d = 1 , the optimal length of the pieces is l ˜ = 7000 m, and in the case where d = 2 , it is l ˜ = 5250 m. The global best choice, according to the AIC, is the latter setting. For the selection of the optimal length we use step sizes of δ l = 250 m.
With the optimal settings, we estimated the parameters for β ^ and ϕ ^ i j as described in Section 2. We estimated the model d × k + p times and, after each step, removed the least significant parameter. Finally, we selected the best estimation according to the AIC. In both cases, the estimation with only one parameter was selected via the AIC. In both cases, the remaining parameter is the one pointing towards Berlin (see Figure 8). Furthermore, β ^ = 0 is selected by our method, which shows that the evolution of the sales prices is not dependent on the distance to Berlin. The result from the real-world study shows that the evolution of the sales prices can be described by a directional influence within the districts of Brandenburg, coming mainly from Berlin.

5. Discussion

Spatial weights matrices for SAR models are of interest in many areas, e.g., for the description of environmental processes such as air pollution or the description of the evolution of housing prices. Numerous methods have therefore been introduced to estimate the spatial weights matrix. However, methods applicable in asymmetric scenarios are lacking.
In this paper, we presented a method for calculating the n × n spatial weights matrix W for directional processes which can be symmetric as well as anti-symmetric. The directional influence is assumed to be the same in all regions of a map and is described by the value of the segments ϕ i j . The proposed method can be used to estimate the directional influence in certain use-cases, such as in the modelling of P M 2.5 in a windy area (see Zhou et al. [22]). In other areas of application, spatial models are created in different environmental areas such as regions of ocean currents (mostly from the same direction). To date, the models have only allowed a symmetric influence or it was necessary to bias the directional influence by knowingly setting certain directions of influence to zero.
In Section 3, we presented Monte Carlo simulations to show that we can re-estimate the ϕ ^ i j consistently. Additionally, we showed that weight matrices calculated with our method outperformed a Queen’s contiguity matrix and an inverse-distance weighting matrix in the selection of four directional settings.
In Section 4, we applied our method to real-world data and showed that there is a short-distance influence in the evolution of land sales prices in Brandenburg, Germany. This directional influence flows from Berlin, at the centre of Brandenburg.
One limitation of the proposed method is that the effort required to calculate the overlapping areas of all pizza-like shapes for all regions is computationally demanding. For larger numbers of regions or larger d and k, the proposed method needs to be improved.
In future work, we intend to decrease the computational effort required to calculate the weights matrices W i j , to allow a larger number of distances d for the segment shapes. With this in hand, one could investigate whether another means of selecting parameters—other than the AIC—could boost the results, since selection using the AIC always promotes fewer parameters. This could be achieved by predicting which area of which region overlap needs to be calculated. Another suggestion is to try different segment shapes. In this paper, we considered only a pizza-like slicing. It could also be worth investigating whether the estimation itself could be enhanced by using a lasso technique.

Author Contributions

F.H.H.: conceptualization, methodology, software, data curation, writing—original draft preparation, writing—review and editing, and visualization; M.S.M.: conceptualization, methodology, and writing—review and editing; P.O.: conceptualization, methodology, writing—original draft preparation, writing—review and editing, supervision, project administration, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All software code is open-source. The data are available upon request.

Acknowledgments

The results presented here were achieved by computations carried out on the cluster system at the Leibniz University Hannover, Germany.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Tobler, W.R. A computer movie simulating urban growth in the Detroit region. Econ. Geogr. 1970, 46, 234–240. [Google Scholar] [CrossRef]
  2. Cressie, N.; Wikle, C.K. Statistics for Spatio-Temporal Data; John Wiley & Sons: Hoboken, NJ, USA, 2015. [Google Scholar]
  3. LeSage, J.P. Bayesian estimation of spatial autoregressive models. Int. Reg. Sci. Rev. 1997, 20, 113–129. [Google Scholar] [CrossRef]
  4. LeSage, J.; Pace, R.K. Introduction to Spatial Econometrics; Chapman and Hall/CRC: Boca Raton, FL, USA, 2009. [Google Scholar]
  5. Billé, A.G.; Rogna, M. The effect of weather conditions on fertilizer applications: A spatial dynamic panel data analysis. J. R. Stat. Soc. Ser. A (Stat. Soc.) 2022. [Google Scholar]
  6. Billé, A.G.; Caporin, M. Impact of COVID-19 on Financial Returns: A Spatial Dynamic Panel Data Model with Random Effects. Available SSRN 3990761 2021. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3990761 (accessed on 9 June 2022).
  7. Billé, A.G. Spatial autoregressive nonlinear models in R with an empirical application in labour economics. In Handbook of Research Methods and Applications in Empirical Microeconomics; Edward Elgar Publishing: Cheltenham, UK, 2021; pp. 23–41. [Google Scholar]
  8. Donegan, C.; Chun, Y.; Griffith, D.A. Modeling community health with areal data: Bayesian inference with survey standard errors and spatial structure. Int. J. Environ. Res. Public Health 2021, 18, 6856. [Google Scholar] [CrossRef]
  9. Lin Lawell, C.Y.C. A Spatial Econometric Approach to Measuring air Pollution Externalities. Available SSRN 675501. 2005. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=675501 (accessed on 9 June 2022).
  10. Boly, A.; Coulibaly, S.; Kéré, E.N. Tax policy, foreign direct investment and spillover effects in Africa. J. Afr. Econ. 2020, 29, 306–331. [Google Scholar] [CrossRef]
  11. Zhao, H.; Cao, X.; Ma, T. A spatial econometric empirical research on the impact of industrial agglomeration on haze pollution in China. Air Qual. Atmos. Health 2020, 13, 1305–1312. [Google Scholar] [CrossRef]
  12. Krisztin, T.; Piribauer, P.; Wögerer, M. The spatial econometrics of the coronavirus pandemic. Lett. Spat. Resour. Sci. 2020, 13, 209–218. [Google Scholar] [CrossRef]
  13. Pinkse, J.; Slade, M.E.; Brett, C. Spatial price competition: A semiparametric approach. Econometrica 2002, 70, 1111–1153. [Google Scholar] [CrossRef]
  14. Stakhovych, S.; Bijmolt, T.H. Specification of spatial models: A simulation study on weights matrices. Pap. Reg. Sci. 2009, 88, 389–408. [Google Scholar] [CrossRef]
  15. Bhattacharjee, A.; Jensen-Butler, C. Estimation of the spatial weights matrix under structural constraints. Reg. Sci. Urban Econ. 2013, 43, 617–634. [Google Scholar] [CrossRef]
  16. Ahrens, A.; Bhattacharjee, A. Two-step lasso estimation of the spatial weights matrix. Econometrics 2015, 3, 128–155. [Google Scholar] [CrossRef] [Green Version]
  17. Lam, C.; Souza, P.C. Estimation and selection of spatial weight matrix in a spatial lag model. J. Bus. Econ. Stat. 2020, 38, 693–710. [Google Scholar] [CrossRef]
  18. Cohen, J.P.; Paul, C.M. The impacts of transportation infrastructure on property values: A higher-order spatial econometrics approach. J. Reg. Sci. 2007, 47, 457–478. [Google Scholar] [CrossRef]
  19. Debarsy, N.; LeSage, J.P. Bayesian model averaging for spatial autoregressive models based on convex combinations of different types of connectivity matrices. J. Bus. Econ. Stat. 2022, 40, 547–558. [Google Scholar] [CrossRef]
  20. Debarsy, N.; Ertur, C. Interaction matrix selection in spatial autoregressive models with an application to growth theory. Reg. Sci. Urban Econ. 2019, 75, 49–69. [Google Scholar] [CrossRef]
  21. Zhang, X.; Yu, J. Spatial weights matrix selection and model averaging for spatial autoregressive models. J. Econom. 2018, 203, 1–18. [Google Scholar] [CrossRef]
  22. Zhou, H.; Jiang, M.; Huang, Y.; Wang, Q. Directional spatial spillover effects and driving factors of haze pollution in North China Plain. Resour. Conserv. Recycl. 2021, 169, 105475. [Google Scholar] [CrossRef]
  23. Merk, M.S.; Otto, P. Estimation of Anisotropic, Time-Varying Spatial Spillovers of Fine Particulate Matter Due to Wind Direction. Geogr. Anal. 2020, 52, 254–277. [Google Scholar] [CrossRef]
  24. Merk, M.S.; Otto, P. Estimation of the spatial weighting matrix for regular lattice data—An adaptive lasso approach with cross-sectional resampling. Environmetrics 2022, 31, e2705. [Google Scholar] [CrossRef]
  25. Anselin, L. Spatial Econometrics: Methods and Models; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1988; Volume 4. [Google Scholar]
  26. Elhorst, J.P.; Lacombe, D.J.; Piras, G. On model specification and parameter space definitions in higher order spatial econometric models. Reg. Sci. Urban Econ. 2012, 42, 211–220. [Google Scholar] [CrossRef]
  27. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control. 1974, 19, 716–723. [Google Scholar] [CrossRef]
  28. Lee, L.F. Consistency and efficiency of least squares estimation for mixed regressive, spatial autoregressive models. Econom. Theory 2002, 18, 252–277. [Google Scholar] [CrossRef]
  29. Gupta, A.; Robinson, P.M. Inference on higher-order spatial autoregressive models with increasingly many parameters. J. Econom. 2015, 186, 19–31. [Google Scholar] [CrossRef] [Green Version]
  30. Longley, P.; Goodchild, M.; Maguire, D.; Rhind, D. Geographic Information Systems and Science; Wiley: Hoboken, NJ, USA, 2005. [Google Scholar]
  31. Sen, Z. Spatial Modeling Principles in Earth Sciences; Springer International Publishing: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
  32. Keitt, T.; Bivand, R.; Pebesma, E.; Rowlingson, B.; Package ‘Rgdal’. Bindings for the Geospatial Data Abstraction Library. 2015. Available online: https://cran.r-project.org/web/packages/rgdal/index.html (accessed on 15 October 2017).
  33. Bivand, R.; Package ‘Spdep’. The Comprehensive R Archive Network. 2015. Available online: https://www.yumpu.com/en/document/view/9283478/package-spdep-the-comprehensive-r-archive-network (accessed on 9 June 2022).
  34. Fingleton, B. A generalized method of moments estimator for a spatial panel model with an endogenous spatial lag and spatial moving average errors. Spat. Econ. Anal. 2008, 3, 27–44. [Google Scholar] [CrossRef] [Green Version]
  35. Osland, L. An application of spatial econometrics in relation to hedonic house price modeling. J. Real Estate Res. 2010, 32, 289–320. [Google Scholar] [CrossRef]
  36. Baltagi, B.H.; Fingleton, B.; Pirotte, A. Spatial lag models with nested random effects: An instrumental variable procedure with an application to English house prices. J. Urban Econ. 2014, 80, 76–86. [Google Scholar] [CrossRef]
  37. Baltagi, B.H.; Bresson, G.; Etienne, J.M. Hedonic housing prices in Paris: An unbalanced spatial lag pseudo-panel model with nested random effects. J. Appl. Econom. 2015, 30, 509–528. [Google Scholar] [CrossRef]
  38. Jin, C.; Lee, G. Exploring spatiotemporal dynamics in a housing market using the spatial vector autoregressive Lasso: A case study of Seoul, Korea. Trans. GIS 2020, 24, 27–43. [Google Scholar] [CrossRef]
  39. Taşpınar, S.; Doğan, O.; Chae, J.; Bera, A.K. Bayesian inference in spatial stochastic volatility models: An application to house price returns in Chicago. Oxf. Bull. Econ. Stat. 2021, 83, 1243–1272. [Google Scholar] [CrossRef]
Figure 1. Example of simulated values for an n = 200 Voronoi map. The simulation is based on 8-setting-4 (see Figure 3).
Figure 1. Example of simulated values for an n = 200 Voronoi map. The simulation is based on 8-setting-4 (see Figure 3).
Symmetry 14 01474 g001
Figure 2. Left: schematic representation of dividing a spatial dependence structure into k directions and d distances. Right: evaluation of an arbitrary segment ( ϕ i j ) of the spatial dependence structure referring to district S 1 . Area 1 is not considered, since it is part of the region that is evaluated. The relative sizes of the overlaps 2 and 3 give the positive weights and are row-normalized to one. Here, we obtain approximations to the weights: w i j , S 1 S 2 = 0.6 and w i j , S 1 S 3 = 0.4 .
Figure 2. Left: schematic representation of dividing a spatial dependence structure into k directions and d distances. Right: evaluation of an arbitrary segment ( ϕ i j ) of the spatial dependence structure referring to district S 1 . Area 1 is not considered, since it is part of the region that is evaluated. The relative sizes of the overlaps 2 and 3 give the positive weights and are row-normalized to one. Here, we obtain approximations to the weights: w i j , S 1 S 2 = 0.6 and w i j , S 1 S 3 = 0.4 .
Symmetry 14 01474 g002
Figure 3. Four different spatial dependence structures implied by the parameters ϕ i j . In the following, they will be denoted 8-setting-1 for the very left setting, to 8-setting-4 on the right. Blue: ϕ i j = 0 , orange: ϕ i j = max , and purple: ϕ i j = max 2 . The value max comes from i = 1 k j = 1 d | ϕ i j | = 0.95 .
Figure 3. Four different spatial dependence structures implied by the parameters ϕ i j . In the following, they will be denoted 8-setting-1 for the very left setting, to 8-setting-4 on the right. Blue: ϕ i j = 0 , orange: ϕ i j = max , and purple: ϕ i j = max 2 . The value max comes from i = 1 k j = 1 d | ϕ i j | = 0.95 .
Symmetry 14 01474 g003
Figure 4. Four different spatial dependence structures implied by the parameters ϕ i j . In the following, they will be denoted 16-setting-1 for the very left setting, to 16-setting-4 on the right. Blue: ϕ i j = 0 , orange: ϕ i j = max , and purple: ϕ i j = m a x 2 . The value max comes from i = 1 k j = 1 d | ϕ i j | = 0.95 .
Figure 4. Four different spatial dependence structures implied by the parameters ϕ i j . In the following, they will be denoted 16-setting-1 for the very left setting, to 16-setting-4 on the right. Blue: ϕ i j = 0 , orange: ϕ i j = max , and purple: ϕ i j = m a x 2 . The value max comes from i = 1 k j = 1 d | ϕ i j | = 0.95 .
Symmetry 14 01474 g004
Figure 5. Results for spatial partitioning to choose the best k ˜ and l ˜ using the four spatial dependence structures from Figure 3. The 12 heatmaps show the best values for each estimation after 1000 realizations. Each heatmap depicts values of k { 4 , 5 , 6 , 7 , 8 , 9 , 10 } on the x-axis and values of l { 225 , 250 , , 525 } on the y-axis. The maximum of the colour scale is chosen for for each heatmap individually to achieve optimal visualization. Especially in 8-setting-3, the selections of k ˜ and l ˜ are very sharp. Therefore, every plot has its own scale, and we named the maximum value max. The blue crosses indicate the true parameters k = 8 and l = 375 .
Figure 5. Results for spatial partitioning to choose the best k ˜ and l ˜ using the four spatial dependence structures from Figure 3. The 12 heatmaps show the best values for each estimation after 1000 realizations. Each heatmap depicts values of k { 4 , 5 , 6 , 7 , 8 , 9 , 10 } on the x-axis and values of l { 225 , 250 , , 525 } on the y-axis. The maximum of the colour scale is chosen for for each heatmap individually to achieve optimal visualization. Especially in 8-setting-3, the selections of k ˜ and l ˜ are very sharp. Therefore, every plot has its own scale, and we named the maximum value max. The blue crosses indicate the true parameters k = 8 and l = 375 .
Symmetry 14 01474 g005
Figure 6. Brandenburg with its n = 2370 districts. Green cross: Berlin centre. Triangles: two example settings of the ϕ i j with two different d and l. They are always chosen to be biased towards Berlin centre. The inset shows the sub-region of Brandenburg with n 0 = 58 districts which is used to determine the parameters { d ˜ , k ˜ , l ˜ } by spatial partitioning.
Figure 6. Brandenburg with its n = 2370 districts. Green cross: Berlin centre. Triangles: two example settings of the ϕ i j with two different d and l. They are always chosen to be biased towards Berlin centre. The inset shows the sub-region of Brandenburg with n 0 = 58 districts which is used to determine the parameters { d ˜ , k ˜ , l ˜ } by spatial partitioning.
Symmetry 14 01474 g006
Figure 7. Heatmap of AIC values for different k, l, d. Lighter colours represent a smaller AIC value.
Figure 7. Heatmap of AIC values for different k, l, d. Lighter colours represent a smaller AIC value.
Symmetry 14 01474 g007
Figure 8. Results of the estimation for Brandenburg. The values represent the estimated values. There is only directional influence from Berlin. This plot visualizes the estimated directional influence between the districts in Brandenburg.
Figure 8. Results of the estimation for Brandenburg. The values represent the estimated values. There is only directional influence from Berlin. This plot visualizes the estimated directional influence between the districts in Brandenburg.
Symmetry 14 01474 g008
Table 1. Means and standard deviations of the best k ˜ and l ˜ for 10 3 estimations. For each setting, the best values are printed in bold. The greater the value of n, the closer we approach to the true values of k = 8 and l = 375 .
Table 1. Means and standard deviations of the best k ˜ and l ˜ for 10 3 estimations. For each setting, the best values are printed in bold. The greater the value of n, the closer we approach to the true values of k = 8 and l = 375 .
8-Setting-18-Setting-28-Setting-38-Setting-4
n200500900200500900200500900200500900
mean( k ˜ )5.346.176.296.86.437.517.99886.266.597.36
sd( k ˜ )1.311.391.351.371.441.070.16001.311.471.22
mean( l ˜ )325336347382393376377375374379389375
sd( l ˜ )866566583730187.14.5693637
Table 2. One thousand estimations of y ^ using this method. The results are compared, using the average AIC and ϵ 2 , to estimations with Queen’s contiguity matrices and inverse-distance weighting matrices. The values of k ˜ and l ˜ are chosen from spatial partitioning.
Table 2. One thousand estimations of y ^ using this method. The results are compared, using the average AIC and ϵ 2 , to estimations with Queen’s contiguity matrices and inverse-distance weighting matrices. The values of k ˜ and l ˜ are chosen from spatial partitioning.
8-Setting-18-Setting-28-Setting-38-Setting-4
n200500900200500900200500900200500900
AICOur Method57014312573568142625625731439259256814252563
Queen67115412846654156929657852015343864416122901
Inv. Dist.65215392859652155228687752031344563615972815
ϵ 2 Our Method0.9540.9871.0000.9470.9840.9951.0161.0351.0400.9570.9860.996
Queen1.3431.2451.3411.2991.3701.4972.5612.9422.5071.3141.3701.381
Inv. Dist.1.2011.2251.3371.2561.2931.3212.3773.1112.5151.2291.3141.268
Table 3. Mean of 500 estimates for each of the 4 different settings with 3 different values of n. We compare the BIAS, RMSE, true positive selections, and false negative selections. The values for BIAS and RMSE are to the power of 10 3 and are separated between values that should be estimated as zero ( = ! 0 ) and values that should be estimated as greater than zero ( > ! 0 ).
Table 3. Mean of 500 estimates for each of the 4 different settings with 3 different values of n. We compare the BIAS, RMSE, true positive selections, and false negative selections. The values for BIAS and RMSE are to the power of 10 3 and are separated between values that should be estimated as zero ( = ! 0 ) and values that should be estimated as greater than zero ( > ! 0 ).
16-Setting-116-Setting-216-Setting-316-Setting-4
n200500900200500900200500900200500900
BIAS = ! 0 1.8350.8920.5951.9201.1030.7290.6680.3230.2432.1131.1070.757
(in 10 3 ) > ! 0 −4.721−3.934−3.677−4.961−4.194−3.824−32.80−28.74−28.06−5.091−4.171−3.855
RMSE = ! 0 1.8610.9100.6092.6361.4450.9400.8880.4070.3312.5711.4300.967
(in 10 3 ) > ! 0 5.1254.1023.8156.7215.0624.15132.8028.7528.067.4905.1014.363
True positive0.860.970.990.780.950.991.01.01.00.760.920.98
False non-zero0.070.060.050.060.050.040.030.020.020.060.050.05
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Harke, F.H.; Merk, M.S.; Otto, P. Estimation of Asymmetric Spatial Autoregressive Dependence on Irregular Lattices. Symmetry 2022, 14, 1474. https://doi.org/10.3390/sym14071474

AMA Style

Harke FH, Merk MS, Otto P. Estimation of Asymmetric Spatial Autoregressive Dependence on Irregular Lattices. Symmetry. 2022; 14(7):1474. https://doi.org/10.3390/sym14071474

Chicago/Turabian Style

Harke, Franz H., Miryam S. Merk, and Philipp Otto. 2022. "Estimation of Asymmetric Spatial Autoregressive Dependence on Irregular Lattices" Symmetry 14, no. 7: 1474. https://doi.org/10.3390/sym14071474

APA Style

Harke, F. H., Merk, M. S., & Otto, P. (2022). Estimation of Asymmetric Spatial Autoregressive Dependence on Irregular Lattices. Symmetry, 14(7), 1474. https://doi.org/10.3390/sym14071474

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop