1# RECORD 549587ac03dcb9d550c7fe68da027277dd876b122# TITLE: A divide-and-conquer approach for spatio-temporal analysis of large house price data from Greater London3# YEAR: 2024.0 VENUE: Annals of Operations Research DOI: 10.1007/s10479-025-06889-x4# AUTHORS: Kapil Gupta; Soudeep Deb5# PAGES: 27 CANDIDATE TABLE PAGES (1-based): [8, 9, 11, 14, 17, 19, 23]6# ABSTRACT: Recent statistical research in real estate markets focuses on the spatio-temporal dynamics of house prices. Bayesian methods are common but often slow for a large dataset in that regard. We propose a divide-and-conquer (D&C) approach, partitioning the data into subsets and applying Gaussian process models in parallel. The results are combined using the Wasserstein barycenter technique to obtain global parameters. This method allows for multiple observations per spatial and time unit. As a real-life application, we analyze house price data of more than 0.6 million transactions from 983 middle layer super output areas in London over a period of eight years. A couple of our major findings is that different types of property command different levels of premiums across the region; and that lower carbon emissions correlate positively with house prices, though the magnitude of this premium has diminished over recent years. Overall, our proposed model offers a better fit than competing approaches, attaining the highest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2$$\end{document} (0.796) and lowest root mean squared error (0.211). Specifically, it improves the error by approximately 50% compared to traditional hedonic and temporal models, while achieving improvements of 28% over fixed-effects model, and improvement of 12% over a spatial model. Most importantly, the D&C spatial variant achieves a 40-fold computational speed-up with minimal accuracy reduction. For unobserved properties, predictions remain highly reliable, achieving a mean absolute percentage error below 10% in over 98% of cases.789===== PAGE 1 (score 8, numeric-paren 0, table-words 1, var-words 6) =====1011A divide-and-conquer approach for spatio-temporal analysis of large house price12data from Greater London13Kapil Gupta, Soudeep Deb14Indian Institute of Management Bangalore15Bannerghatta Main Road, Bangalore, Karnataka 560076, India.16Abstract17Statistical research in real estate markets, particularly in understanding the spatio-temporal dynamics18of house prices, has garnered significant attention in recent times. Although Bayesian methods are19common in spatio-temporal modeling, standard Markov chain Monte Carlo (MCMC) techniques are20usually slow for large datasets such as house price data. To tackle this problem, we propose a divide-21and-conquer spatio-temporal modeling approach.22This method involves partitioning the data into23multiple subsets and applying an appropriate Gaussian process model to each subset in parallel. The24results from each subset are then combined using the Wasserstein barycenter technique to obtain the25global parameters for the original problem. The proposed methodology allows for multiple observations26per spatial and time unit, thereby offering added benefits for practitioners. As a real-life application,27we analyze house price data of more than 0.6 million transactions from 983 middle layer super output28areas in London over a period of eight years. The methodology provides insightful findings about the29effects of various amenities, trend patterns, and the relationship between prices and carbon emissions.30Furthermore, as demonstrated through a cross-validation study, it shows good predictive accuracy31while balancing computational efficiency.32Keywords:33Bayesian analysis, Real estate, Space-time prediction, Wasserstein barycenter341. Introduction35The valuation of residential properties is immensely important for various economic stakeholders,36including property dealers, sellers, and buyers. Traditionally, house price modeling has revolved around37the hedonic framework (Rosen, 1974).38Specifically, log-linear regression techniques have been at39the forefront in the analysis of house price transactions (Liu, 2013). However, the incorporation of40location attributes into the mean structure of hedonic regression, as recommended by Dubin and Sung41(1990), has shown limitations in capturing all significant neighborhood attributes, potentially giving42rise to spatial autocorrelation concerns. To overcome this challenge, researchers have advocated for the43inclusion of spatial effects into the error structure (Can, 1990; Pace and Gilley, 1998). Simultaneously,44the dynamics of house prices are profoundly influenced by temporal autocorrelation as well.45Pace46et al. (1998) leveraged the concepts of spatio-temporal modeling to comprehensively capture both47spatial and temporal effects within house price data.48This modeling approach has since become49a cornerstone in academic research. Building on this foundation, Liu (2013) harnessed the spatio-50temporal autoregressive (STAR) model to address correlated errors in traditional hedonic regression,51delving into spatial and temporal dependence considerations and further developing a house price52index as a pivotal part of their analysis. Interestingly, their methodology assumed the spatial weight53matrix to be lower triangular, which implies that two spatial units cannot mutually influence each54other. A different approach was taken by Fotheringham et al. (2003) who integrated geographically55weighted regression with time series forecasting techniques to capture the spatio-temporal variations56arXiv:2407.15905v1 [stat.AP] 22 Jul 2024575859===== PAGE 8 (score 7, numeric-paren 4, table-words 0, var-words 3) =====6061a separable structure for the spatial dependence and the temporal dependence. In other words, it is62written as a product of purely spatial and purely temporal covariance functions:63Cov {v(si, tu), v(sj, tv)} = σ264vρs (∥si −sj∥; ϕs) ρt (∥tu −tv∥; ϕt) .65(3.4)66The separable covariance structure is favored as it significantly reduces the parameters in the67covariance matrix, and is a popular choice in similar spatio-temporal studies (see, e.g., Sahu et al.,682010; Deb and Tsay, 2019). For both spatial and temporal covariance functions, we consider the Matérn69structure with the specific choice of ν = 0.5, resulting in the exponentially decaying pattern. Thus, we70have:71ρs (∥si −sj∥; ϕs) = exp (−ϕs ∥si −sj∥) ;72ρt (∥tu −tv∥; ϕt) = exp (−ϕt ∥tu −tv∥) ,73(3.5)74where ∥tu −tv∥represents the absolute time difference, and ∥si −sj∥is calculated via the Vincenty75Ellipsoid method, which computes the great circle distance on an ellipsoid. The latter is chosen for its76superior accuracy compared to other approaches (Hijmans, 2021). It is critical to note that, because77of the separability assumption imposed on the joint covariance function, we can write the covariance78matrix Σ = σ279v (Σs ⊗Σt), where ⊗denotes the Kronecker product, and Σs, Σt respectively denote the80spatial and temporal correlation matrices, with81(Σs)ij = exp (−ϕs ∥si −sj∥) ;82(Σt)uv = exp (−ϕt ∥tu −tv∥) .83(3.6)843.2. Prior Specifications85Let θ = (β⊤, σ286v, σ287ϵ , ϕs, ϕt)⊤be the vector of unknown parameters. We estimate them completely88within a Bayesian framework, primarily using Gibbs sampling (?), for which suitable prior specifications89are essential.90For all coefficients in β, we consider independent Gaussian priors with mean 0 and91variance c, where c is relatively large, typically on the order of 104. For the error variances σ292ϵ and93σ294v, we employ independent inverse gamma priors IG(a, λ), where a and λ are the shape and scale95parameters, respectively. We choose a = 2, making these priors non-informative. We also noticed that96the parameter estimates and model selection criteria are insensitive to the initial values of the scale97parameter, and throughout this study, we take λ = 1.98For estimating the decay parameters ϕ = (ϕs, ϕt) in the covariance functions, many prior studies99(e.g., Sahu et al., 2006; ?) have used cross-validation to find the optimal choices. Albeit this approach100is computationally more tractable, it has some limitations, as it explores a small set of values and101can miss viable solutions. An alternative is the single-variable slice sampling algorithm (Neal, 2003),102which provides direct estimation of the parameters within a Bayesian framework. While this method103thoroughly explores the parameter space, it demands extensive computation as the covariance matrix104must be recalculated and inverted in every iteration of the Bayesian computation.105Keeping the106advantages of both approaches in view, we aim to achieve a middle ground – somewhat like a discrete107slice sampling – wherein a discretized sample space is considered for the parameters (ϕs, ϕt), and the108parameters are estimated directly within the Bayesian setting. This approach allows us to provide109credible intervals for these parameters. In terms of the prior specifications, we adopt a discretized110uniform prior for ϕs, selecting values from 1 to 3 in intervals of 0.2, resulting in the set As. For ϕt,111similarly, we select equidistant points between 0.2 and 1 in intervals of 0.2, and denote it as At. Below,112the product As × At is indicated as Ast, while U indicates the mentioned uniform distribution.113As mentioned before, we use concepts of Gibbs sampling and discrete slice sampling to implement114the proposed model for estimation and prediction purposes. As the first step, following Section 3.1, we115can write the full model as116Y | θ, V ∼NN117 Xβ + BV , σ2118ϵ IN119120,121(3.7)1228123124125===== PAGE 9 (score 17, numeric-paren 3, table-words 4, var-words 6) =====126127while the prior distributions for all parameters are specified as:128β ∼Nm+1 (0, cIm+1) ,129σ2130ϵ ∼IG(a, λ), σ2131v ∼IG(a, λ),132(ϕs, ϕt) ∼U(Ast),133V ∼NG134 0, σ2135vΣs ⊗Σt136137.138(3.8)139Using f(·) to denote the density function in a generic way, the complete posterior distribution is140f141 β, V , σ2142ϵ , σ2143v, ϕs, ϕt | Y144145∝f146 Y | β, V , σ2147ϵ , σ2148v, ϕs, ϕt149150f151 V | σ2152v, ϕs, ϕt153154f(σ2155v)f(ϕs)f(ϕt)f(σ2156ϵ )f(β).157(3.9)158For brevity, detailed calculations for the conditional posteriors in case of the complete data are159deferred to Appendix A.1 of the supplementary material.160In the next subsection, we present the161posterior distributions and calculations for the D&C setup.1623.3. Divide-and-conquer approach for large dataset163The above model and its Bayesian implementation are helpful in analyzing house price dynamics.164However, even with the assumption of separability and the use of Kronecker product to reduce the165cost of inverting dispersion matrices, this approach becomes computationally challenging for a dataset166with large number of locations or time-points. Specifically, each iteration in the above implementation167requires O(S2T 2) memory units for storage and O((S ∨T)3) floating-point operations for computation.168Naturally, this approach is infeasible for large-scale data, which motivates us to develop a divide-and-169conquer approach in this section building upon the above methodology. Although the use of divide-and-170conquer is not a new concept in computer science, its implementation in scalable Bayesian inference171is a recent advancement.172We present the technique in a general three-step distributed framework173applicable for scaling posterior computations in our spatio-temporal model for house prices presented174in Section 3.1.1753.3.1. First step: partitioning the data176Consider a large dataset with S spatial locations, T time-points, and multiple observations as set177out before. Our approach involves splitting these data points into Q non-overlapping subsets for model178feasibility. We use a random allocation scheme for the S locations into these Q subsets. This process179ensures that each subset contains an adequate number of locations and their corresponding data for180all time-points. To avoid bias from locations with significantly more or fewer observations, one should181take precautions during the random assignment, striving for balanced representation across the subsets.182Our specific approach for the main analysis is elaborated in Section 4.183To set the notations for this section, let Sq = (sq1, sq2, . . . sqmq ) denote the set of mq locations in184the qth subset, q ∈{1, . . . , Q}. Although in principle a spatial location can be included in multiple185subsets, in this work we only consider disjoint partition of the spatial index set, thereby implying186PQ187q=1 mq = S and ∪Q188q=1Sq being equal to the set of all locations. In the qth subset, we have the189data {Y(q), X(q)}, with the total number of observations being N(q) = Pmq190i=1191PT192j=1 kij. Clearly, Y(q)193is a N(q) × 1 vector, and X(q) is a N(q) × (m + 1) matrix representing the covariates corresponding194to the observations for all time-points in the spatial locations within the subset Sq. Following earlier195notations, the spatio-temporal Gaussian process model for the qth subset is:196Y(q) = X(q)β(q) + B(q)V(q) + ϵ(q).197(3.10)198Throughout the discussions in this section, the parameters θ(q) = (β(q), σ2199v(q), σ2200ϵ(q), ϕs(q), ϕt(q))⊤,201the data vector Y(q), the spatio-temporal process V(q), the matrices X(q), B(q), and all the available2029203204205===== PAGE 10 (score 5, numeric-paren 3, table-words 0, var-words 2) =====206207information H(q) for subset q are equivalent to their full dataset counterparts θ = (β, σ2208v, σ2209ϵ , ϕs, ϕt)⊤,210Y , V , X, B, and H. Also, we use the notation pq = N/N(q) to simplify some expressions.2113.3.2. Second step: modified sampling from subset pseudo posterior distributions212We note that the qth subset described above contains a fraction (1/pq) of the full data. Thus,213using its posterior will result in an overestimation of the posterior uncertainty compared to using the214full dataset. To address this mismatch, we consider a modification to the Bayesian algorithm on the215subsets without sacrificing its efficiency. Particularly, a powered likelihood is applied to modify the216likelihood of all the parameters before employing the sampling algorithm on the qth subset. As it217contains 1/pq fraction of the full data, the asymptotic variances of the posterior distributions of θ(q)218and V(q) are inflated by a factor of pq compared to that of θ and V (Shyamalkumar and Srivastava,2192022). Subsequently, we employ stochastic approximation by raising the likelihood of the parameters220β(q), σ2221v(q), σ2222ϵ(q), ϕs(q), ϕt(q), V(q) in the qth subset to the power of pq (Minsker et al., 2014).223These224adjustments account for the data present in the other subsets, and ensure accurate estimation of225posterior uncertainty while maintaining computational efficiency.226For the qth subset, as there are N(q) data points, the covariance matrices Σs(q) and Σt(q) have227dimensions mq × mq and T × T, respectively. Let G(q) = mq × T, and K(q) be a diagonal matrix228of order G(q) × G(q), where the diagonal elements are kij for different (si, tj) combinations in the qth229subset. The prior distributions for all parameters in subset q align with those in (3.8), except for V(q),230which has the following structure:231V(q) ∼NG(q)2322330, σ2234v(q)235 Σs(q) ⊗Σt(q)236237.238(3.11)239Utilizing the abovementioned stochastic approximation, we now obtain the following subset pseudo240posterior distribution. For the sake of brevity, we omit the technical derivations in the main text, and241refer to the supplementary material for detailed information.242β(q) | H(q) ∼Nm+1243244245 pqX⊤246(q)X(q)247σ2248ϵ(q)249+ Im+1250c251!−1 pqX⊤252(q)253 Y(q) −B(q)V(q)254255σ2256ϵ(q)257! pqX⊤258(q)X(q)259σ2260ϵ(q)261+ Im+1262c263!−1264,265σ2266ϵ(q) | H(q) ∼IG267 268a + N2692 , pq270
Y(q) −X(q)β(q) −V(q)271
22722273+ λ274!275,276σ2277v(q) | H(q) ∼IG278279a + pqG(q)2802281,282pq283284V ⊤285(q)286287eΣ−1288s(q) ⊗Σ−1289t(q)290291V(q)2922932294+ λ295296,297f(ϕs(q), ϕt(q)) | H(q) ∝|Σs(q)|−pqT/2|Σt(q)|−pqmq/2 exp298299300−pqV ⊤301(q)302303eΣ−1304s(q) ⊗Σ−1305t(q))V(q)306307σ2308v(q)309310I311 (ϕs(q), ϕt(q)) ∈Ast312313,314V(q) | H(q) ∼NG315316317 eΣ−1318s(q) ⊗Σ−1319t(q)320σ2321v(q)322+ pqK(q)323σ2324ϵ(q)325!−1 pqB⊤326(q)(Y(q) −X(q)β(q))327σ2328ϵ(q)329!330,331 eΣ−1332s(q) ⊗Σ−1333t(q)334σ2335v(q)336+ pqK(q)337σ2338ϵ(q)339!−1340.341(3.12)342To maintain consistency in our analysis, we adopt identical potential values for decay parameters343across all subsets.344Therefore, we utilize the same combination of ϕs(q) and ϕt(q).345An essential346modification in our posterior calculation involves capturing spatial dependence between subsets. This is347achieved by adjusting the inverse of the subset covariance matrix to account for spatial interdependencies34810349350351===== PAGE 11 (score 9, numeric-paren 4, table-words 0, var-words 5) =====352353among the subsets. To properly follow the procedure, one typically needs to invert the entire spatiotemporal354matrix (or the spatial dependence matrix under the separability assumption) before subsetting it for355the qth subset. However, inverting such a large matrix is very time-consuming. To address this, we356use the concept of Schur’s complement (Zhang, 2006) and find the inverses of the diagonal matrices357efficiently. To elaborate this approach, note that we first obtain the inverse of the last diagonal part358using Schur’s complement, which involves the inverse of a small partitioned matrix. The other diagonal359part is also the inverse of a matrix comprising only the inverse of a small matrix. In the next iteration,360we treat this inverse part as the inverse of the updated main matrix and repeat the procedure until361only one diagonal matrix remains. By employing this iterative method, we can find the inverses of362all the block diagonal matrices without needing to invert the entire matrix.363Our extension of the364Schur’s complement approach significantly speeds up the computation of the inverse. Let us denote365the inverse of Σs(q) as eΣ−1366s(q), reflecting the modified version after incorporating spatial dependence.367Detailed calculations of this procedure are shown in Appendix B of the supplementary material.368In our Gibbs sampling procedure, computing the posterior of the vector V(q) requires inverting a369G(q) × G(q) matrix in each iteration, which can be computationally demanding. To mitigate this, we370adopt a more efficient approach based on the theory of multivariate normal distribution. Specifically,371we partition the vector V(q) into mq subvectors based on the columns in the matrix V(q). Then, we372take Z373′3741 as the first subvector of V(q), and Z375′3762 as a concatenated subvector of the remaining columns377of V(q), i.e., V(q) =378379Z380′3811 Z382′3832384385. Similarly, the covariance matrix Σ(q) can be written as386Σ(q) = σ2387v(q)388Σ11(q) ⊗Σt(q) Σ12(q) ⊗Σt(q)389Σ21(q) ⊗Σt(q) Σ22(q) ⊗Σt(q)390391.392(3.13)393Then, the conditional distribution of Z394′3951 given Z396′3972 = z′3982 is NT (µc(q), σ2399v(q)Σc(q)), where400µc(q) =401 Σ12(q) ⊗Σt(q)402 403Σ−140422(q) ⊗Σ−1405t(q)406407z′4082;409Σc(q) =410 Σ11(q) ⊗Σt(q)411412−413 Σ12(q) ⊗Σt(q)414 415Σ−141622(q) ⊗Σ−1417t(q)418 Σ21(q) ⊗Σt(q)419420.421(3.14)422Using the properties of the Kronecker product, the expression can be further simplified as423µc(q) =424425Σ12(q)Σ−142622(q)427428⊗IT429430z′4312;432Σc(q) =433434Σ11(q) −Σ12(q)Σ−143522(q)Σ21(q)436437⊗Σt(q).438(3.15)439Now, following the above, we can write440V1(q) | H(q) ∼NT441442443 Σ−1444c(q)445σ2446v(q)447+ pqK1(q)448σ2449ϵ(q)450!−1 451Σ−1452c(q)µc(q)453σ2454v(q)455+456pq457458BT4591(q)460 Y1(q) −X1(q)β(q)461462σ2463ϵ(q)464465,466 Σ−1467c(q)468σ2469v(q)470+ pqK1(q)471σ2472ϵ(q)473!−1474,475(3.16)476where Y1(q) be the data for the first location and all time-points from the qth subset. The matrix X1(q)477is defined accordingly. K1(q) is a T × T diagonal matrix with each diagonal element being equal to478k1j for each (sq1, tj). B1(q) is a submatrix of B(q), consisting of the first T columns and PT479j=1 k1j480rows, resulting in a sparse matrix of size PT481j=1 k1j × T. Other notations have the same definitions482as before. We can compute the conditional posterior distributions of V2(q), . . . , Vmq(q) in an identical483fashion. Finally, to compute the posterior distribution of the full vector V(q), we iteratively employ the48411485486487===== PAGE 12 (score 6, numeric-paren 2, table-words 0, var-words 4) =====488489conditional distributions in a Gibbs sampler. This approach generates independent posterior samples490after an initial burn-in period. Consequently, the posterior distribution of V(q) is computed efficiently491without incurring the computational cost of inverting the G(q) × G(q) matrix in each Gibbs sampler492iteration. It is imperative to point out that these subset pseudo posteriors significantly speed up the493computations and can be run in parallel for all partitions.4943.3.3. Third step: combination of subset posteriors495The final and most crucial step of the D&C algorithm is to combine the subset posterior draws.496As discussed in Section 1, various existing methods can be employed in this regard. We follow the497Wasserstein Barycenter (WB) technique to approximate the complete posterior distribution.498This499approach has been successfully utilized by Li et al. (2017) and Srivastava et al. (2018) in the setting500of independent data, whereas Ou et al. (2021) demonstrated the accuracy of using this for combining501subset posteriors in the case of time series data. Motivated by these studies, we apply a similar rationale502for spatio-temporal data and utilize the WB to combine the subset posteriors. It is worth noting that503Wasserstein-based posteriors exhibit appealing asymptotic properties, as shown by Szabó et al. (2019).504Before going into the detail, it is critical to discuss the mathematical underpinning of the Wasserstein505distance.506Let (M, d) be a complete separable metric space, P(M) be the space of all probability507measures on M and P2(M) be the set of all probability measures on M with finite second moments.508For ϱ, ξ ∈P2(M), Π(ϱ, ξ) denotes the set of all probability measures on M × M with marginals ϱ and509ξ. Then, with d(·, ·) indicating the euclidean metric, the Wasserstein distance of order 2 between ϱ and510ξ is defined as511W2(ϱ, ξ) =512513inf514π∈Π(ϱ,ξ)515Z516M×M517d(x1, x2)2Π(x1, x2)5181/2519.520(3.17)521As Bickel and Freedman (1981) discussed, the convergence in W2 on P2(M) implies weak convergence522of probability measures and convergence of the second moment. Keeping that in view, if θ(1), . . . , θ(Q)523are the Q subset posterior distributions, then the approximated posterior ˆθ is defined as524ˆθ = arg min525θ∈P2(M)5261527Q528Q529X530i=1531W 25322 (θ, θ(i)).533(3.18)534This approach is referred to as the Wasserstein Posterior (WASP). The exact computation of the535WASP is known to be computationally intensive and remains a subject of ongoing research. We adopt536an approximation for ˆθ, following Shyamalkumar and Srivastava (2022). This completes the last step537in our proposed approach (D&C-STBD), and the pseudocode of the entire procedure is shown in538Algorithm 1.539A brief outline of the WASP procedure is also of the essence here.540First, posterior draws of541θ(1), . . . , θ(Q) are taken from the Q subsets using the Bayesian algorithm developed earlier.542In all543instances of Gibbs and discrete slice samplers, the convergence of the chains are confirmed through the544diagnostic method of Geweke et al. (1991). After convergence, we apply thinning to collect posterior545samples for each subset. These sample points (say, L for each subset) are then consolidated into a546unified set of QL samples. Using these posterior samples, we compute mean vectors ˆµ(q) and covariance547matrices ˆΣ(q) for each subset. This results in a total of Q mean vectors and covariance matrices. Next,548we take the mean of these and get the overall mean ˆ¯µ and covariance matrix ˆ¯Σ. Now, using the subset549mean vectors and covariance matrices, we standardize and center the subset posterior samples. Then,550we use the overall mean vector and covariance matrix to get the approximate WASP draws for the551actual posterior draws, which leads to QL number of WASP posterior samples for each parameter. We552estimate both the posterior mean of these parameters and their corresponding credible intervals using553the consolidated sample. On the other hand, for the final estimate of V from the subset posterior55412555556557===== PAGE 14 (score 11, numeric-paren 8, table-words 0, var-words 3) =====558559follows a Gaussian distribution:560y(s′, t′, r′) | v(s′, t′) ∼N561562X⊤(s′, t′, r′) ˆβ + v(s′, t′), σ2563ϵ564565,566(3.19)567We estimate the model parameters using a Gibbs sampler on the training data, yielding posterior568estimates for β and σ2569ϵ . To find the conditional distribution of v(s′, t′) given V , we use the posterior570value of V from training data estimates and apply the multivariate Gaussian theory related to the571conditional distribution. We treat v(s′, t′) as the first subvector of the multivariate joint distribution,572and V obtained from the training dataset as the second subvector. The joint distribution of v(s′, t′)573and V is given by574v(s′, t′)575V576577∼NG+157805790580581, σ2582v583 1584Σ′⊤585Σ′586Σs ⊗Σt587588, with Σ′ = Σs(s −s′) ⊗Σt(t −t′),589(3.20)590where Σs(s−s′) is a vector of order S with ith entry as ρs (||si −s′||; ϕs) and Σt(t−t′)is a vector of order591T with jth entry as ρt (||tj −t′||; ϕt). The conditional distribution is then given by592v(s′, t′) | V ∼N593594Σ′(Σ−1595s596⊗Σ−1597t )V , σ2598v(1 −Σ′(Σ−1599s600⊗Σ−1601t )Σ′⊤)602603.604(3.21)605Thus, by predicting v(s′, t′) using the above conditional distribution and the estimates of β from606the Gibbs sampling procedure, we can predict y(s′, t′, r′) using (3.19).607Now, we shift our focus to the prediction stage in D&C case.608Keeping parity with the earlier609notations, we need to forecast v(s′, t′) for a location s′ and a time-point t′ in the test dataset. Our610training dataset is divided into Q subsets, each treated as a separate training dataset. For each subset,611we predict v(s′, t′) using a similar process as for the whole dataset, resulting in Q realizations of v(s′, t′).612These realizations are denoted as vq(s′, t′), for 1 ⩽q ⩽Q. The joint distribution of vq(s′, t′) and V(q)613is given by614vq(s′, t′)615V(q)616617∼NG(q)+1618 06190620621, σv2622(q)623"6241625Σ′⊤626(q)627Σ′628(q)629Σs(q) ⊗Σt(q)630#!631, with Σ′632(q) = Σs(q)(s −s′) ⊗Σt(q)(t −t′),633(3.22)634where Σs(q)(s−s′) is a vector of order S(q) with ith entry as ρs635 ||si −s′||; ϕs(q)636637and Σt(t−t′)is a vector638of order T with jth entry as ρt639 ||tj −t′||; ϕt(q)640641. The conditional distribution is then given by642vq(s′, t′) | V(q) ∼N643644Σ′645(q)(eΣ−1646s(q) ⊗Σ−1647t(q))V(q), σ2648v(q)(1 −Σ′649(q)(eΣ−1650s(q) ⊗Σ−1651t(q))Σ′⊤652(q))653654.655(3.23)656With the above predicted values of v(q)(s′, t′) and the estimated parameters within each subset q,657we can predict the response variable, i.e. the price of a house in logarithmic scale, for any MSOA at a658future time-point in the test data using (3.24), which is based on the qth subset data.659y(q)(s′, t′, r′) | v(q)(s′, t′) ∼N660661X⊤662(q)(s′, t′, r′) ˆβ(q) + v(q)(s′, t′), σϵ2663(q)664665.666(3.24)667Finally, to predict v(s′, t′) and the response variable y(s′, t′, r′) using the entire training dataset, we668take the median of vq(s′, t′) and y(q)(s′, t′, r′) respectively, for 1 ⩽q ⩽Q. This approach leverages the669D&C framework to efficiently handle large datasets by breaking them into smaller, manageable subsets670while maintaining prediction accuracy. It is noteworthy that we incorporate the original inverse of the671spatial covariance matrix for each subset to address spatial dependence.67214673674675===== PAGE 15 (score 8, numeric-paren 1, table-words 3, var-words 1) =====6766774. Model implementation and evaluation678In this section, we explore the key implementation steps for our model and discuss other competing679models along with the evaluation criteria used for benchmarking.680Recall that our dataset comprises data on N = 651202 transactions from S = 983 locations and681T = 106 time-points. As covariates, we include the variables listed in Table 1, along with linear and682quadratic time trends, as well as the interaction of time with carbon emissions. We take logarithmic683transformation for all the numeric variables.684In the first step of our methodology, the dataset is685partitioned into Q = 20 disjoint subsets. The initial 19 subsets are comprised of randomly selected68649 locations, while the last subset contains the remaining 52 locations. Each subset encompasses data687for all the time-points and their corresponding multiple records. Next, for each subset, we employ688the model specified in (3.10) and fit the models separately using the posterior calculations shown689there.690After convergence of the chains, 2000 MCMC samples are taken in each of the Q subsets.691Finally, we adopt the subset posterior combination method, as described in Section 3.3.3, to obtain692final estimates for all the parameters. Additionally, we determine the values of v(s, t) corresponding693to each combination of (s, t) in the dataset.694To assess the fit of the proposed approach, we conduct a comparative analysis with four alternative695models. The first model is a standard hedonic regression that employs the same covariates as our696approach. As discussed in Section 1, it is one of the most popular approaches in related research.697The second model incorporates an additional fixed effect for the geographic region where the MSOA698is situated. We refer to this model as the spatial hedonic model, and it closely parallels the approach699described in ?. While these two approaches are conventional techniques in the extant literature of house700price analysis, we include two other models in the comparison study for understanding the effectiveness701of the proposed spatial and temporal dependence structures. The third model incorporates temporal702dependence in a manner similar to the main model proposed in this paper, but without the spatial703component. Given that there are 106 time points, the inversion of a 106 × 106 matrix during iteration704is computationally efficient, making this model’s computation manageable without the D&C step. The705fourth model, on the other hand, does not consider temporal dependence but introduces only a spatial706error process. With 983 locations, matrix inversion in this model could be time-consuming in each707iteration. Therefore, this spatial error process model is integrated into the proposed D&C framework708and implemented in the proposed fashion. This model will be referred to as D&C-SBD. By comparing709our method with these models, we aim to assess their respective data-fitting capabilities. Specifically,710we evaluate the accuracy of each model’s fit to the data and examine whether the resulting estimates are711consistent with existing literature. This comparison will help determine the robustness and reliability712of our approach in relation to established models.713We evaluate the efficacy of these models by fitting the data. Label the entire dataset as D, with714an observation for location si at time-point tk be denoted as y(si, tk, rik). The corresponding fitted715value is ˆy(si, tk, rik). To understand the goodness of fit of every model, we start with the ubiquitous716coefficient of determination (R2) which quantifies the amount of variability in the response variable717that can be accounted for by the predictor variables, thereby indicating the degree of adequacy of718the model.719Further, we consider the mean absolute error (MAE), mean absolute percentage error720(MAPE) and root mean square error (RMSE) of the fitted values, as recommended in statistical721literature (Botchkarev, 2019). The MAE, i.e. the average absolute difference between actual and fitted722values, provides a sense of the extent of the errors without considering their direction (overestimation723or underestimation).724MAPE, on the other hand, computes the percentage deviation and offers a725straightforward interpretation of the relative error. RMSE is one of the most popular metrics and726is advantageous because it penalises large errors more severely than small errors, thereby providing a727more impartial evaluation of the models’ accuracy. Finally, in some cases, a model with a lower average728error (as measured by MAE, MAPE, or RMSE) may exhibit higher variability in errors. To that end,72915730731732===== PAGE 17 (score 21, numeric-paren 0, table-words 3, var-words 9) =====733734Table 2: Parameter estimates (posterior mean) and the corresponding credible intervals (CI) when the proposed model735is fitted to the whole data with 983 MSOAs and 106 time-points (651202 observations).736Variable737Estimate73895% CI739Intercept7409.675741(9.640 , 9.711)742Linear time trend7430.982744(0.845 , 1.103)745Quadratic time trend746−0.646747(−0.757, −0.534)748log (1 + Area)749−0.319750(−0.321, −0.316)751Rooms (3 Rooms)7520.066753(0.065 , 0.068)754Rooms (4 Rooms)7550.059756(0.057 , 0.062)757Rooms (5 Rooms)7580.064759(0.061 , 0.067)760Rooms (> 5 Rooms)7610.109762(0.106 , 0.113)763Property type (Terraced)764−0.052765(−0.054, −0.050)766Property type (Detached)7670.141768(0.138, 0.144)769Property type (Flats)770−0.250771(−0.252, −0.247)772Fireplace7730.064774(0.062, 0.065)775Ventilation (Extract-only)7760.029777(0.022 , 0.036)778Ventilation (Both supply and extract)7790.055780(0.043 , 0.065)781Wind turbine7820.006783(−0.019 , 0.031)784log (1 + CO2 emission)785−0.031786(−0.035, −0.028)787Interaction of trend and CO2 emission7880.076789(0.070 , 0.081)790σ2791ϵ7920.043793(0.043 , 0.043 )794σ2795v7960.083797(0.070 , 0.103)798ϕs7992.402800(1.486 , 3.041)801ϕt8020.528803(0.183 , 0.824)8049.08059.580610.080710.580811.0809Jan 2011810Jul 2011811Feb 2012812Sep 2012813Apr 2013814Nov 2013815Jun 2014816Jan 2015817Aug 2015818Mar 2016819Oct 2016820May 2017821Dec 2017822Jul 2018823Feb 2019824Oct 2019825Month826log price827Fig. 5: Average trend of the house price (log-transformed) for all locations over 106 time.8285.1.1. Effect of area and number of rooms829Our analysis unveils a straightforward relationship between property size and price in London. The830negative coefficient associated with log area indicates that as a property’s size increases, its value per83117832833834===== PAGE 18 (score 17, numeric-paren 0, table-words 2, var-words 13) =====835836square meter tends to decrease. This observation is in line with common expectations and matches837the results reported in extant literature, e.g., Chegut et al. (2016).838Furthermore, the number of rooms within a property emerges as a significant factor. Properties839with more than two rooms command notably higher prices compared to those with 1 or 2 rooms. This840effect becomes more pronounced for properties with over five rooms, where the impact on price per841square meter is most substantial. It highlights the premium placed on larger, more spacious homes in842the London region, reflecting the diverse housing preferences of potential buyers.843In summary, larger properties with same number of rooms tend to have lower prices per square844meter than smaller properties. This finding aligns with the results reported in Fuerst et al. (2015).8455.1.2. Effect of property characteristics846The type of property stands out as a pivotal factor influencing prices. We find that terraced and flat847properties tend to have lower prices per square meter compared to semi-detached houses. Strikingly,848flats are at the lower end of the price spectrum, while detached houses command the highest prices.849This finding aligns with the study conducted by Fuerst et al. (2015), where terraced properties were850considered the base category, and both detached and semi-detached properties were found to be on the851premium side, while flats were on the lower end. This result underscores the importance of property852type as a driver of price volatility in the London real estate market and reflects the diverse architectural853landscape of the region.854The presence of fireplaces within properties is strongly associated with higher prices. It aligns855with the perception of fireplaces as desirable amenities for homeowners, contributing positively to856property values, and is consistent with the observations of Zhang (2016). Meanwhile, we find that857the ventilation systems exert a discernible influence on house prices. Properties equipped with extract858mechanical ventilation and both supply and extract mechanical ventilation systems tend to exhibit859higher prices compared to those reliant on natural ventilation. One may connect this to the growing860emphasis on indoor air quality and environmental factors among property buyers. It reinforces the861desirability of properties that offer superior air quality, facilitated by mechanical ventilation systems.8625.1.3. Effect of carbon emission863Another intriguing aspect of our analysis is the investigation of the relationship between carbon864emissions and house prices.865The increasing need of a cleaner environment is likely to impact the866economy of real estate and we attempt to capture the evolving dynamics of sustainability in the867London housing market through this analysis. By considering appropriate regressors, we delve into868the multifaceted aspects of this relationship, including its impact, interaction with time, and the869implications for supply and demand.870Impact of carbon emission on house price: Our findings highlight the significant impact of871carbon emissions on house prices. We observe a negative relationship, indicating that properties with872lower carbon emission values tend to command higher prices. This emphasizes the growing importance873of sustainability and environmental consciousness among real estate buyers. It suggests that buyers are874increasingly willing to invest more in homes with a reduced carbon footprint, reflecting the evolving875trend toward eco-friendly living spaces. The findings align with those of Gerassimenko et al. (2023),876who noted that energy-efficient properties, characterized by low carbon emissions, tend to command a877price premium.878The interaction with time: A more nuanced observation emerges when we consider the interaction879between carbon emissions and time. Over time, we observe a decline in willingness to pay premium880for lower carbon emissions, reflecting a shift in buyer behavior. This observation is supported by the881positive coefficient associated with the interaction term. In practice, while carbon emissions continue882to be an important factor for buyers, their impact on house prices has gradually diminished over the88318884885886===== PAGE 19 (score 16, numeric-paren 0, table-words 5, var-words 3) =====887888years. This trend raises interesting questions about the evolution of buyer preferences and the broader889real estate market.890Supply and demand implications: Finally, the interplay between carbon emissions and house891prices also has implications for supply and demand dynamics. As sustainability becomes a central892consideration for buyers, properties with lower carbon emissions are likely to see increased demand.893This heightened demand may prompt the builders to focus on constructing more environmentally894friendly properties, further influencing the supply side of the equation. However, our observation that895the price differential between eco-friendly and conventional properties is diminishing suggests that the896market is finding a balance. As supply increases and the premium on sustainability lessens, a broader897range of buyers can access these eco-friendly properties (see Ma and Li, 2017; Blanco and Neri, 2023,898for pertinent discussions). This equilibrium benefits both buyers and the environment.8995.2. Comparison between the fits of different models900As mentioned in Section 4, we conduct a comprehensive comparative study to assess the fitting901and accuracy of different models. Table 3 displays the accuracy metrics on overall data to evaluate902the goodness of fit. The results unequivocally demonstrate that our proposed model outperforms the903other models, indicating its superior ability to capture the spatial and temporal dependencies in the904data and achieve a better fit. Notably, D&C-STBD explains the highest variation, while the temporal905model falls short in explaining the variance effectively. We also observe that the D&C-SBD method906requires significantly less time compared to the spatio-temporal model or the spatial hedonic model,907yet it provides a great fit. This highlights the efficiency and accuracy of the spatial component in the908method, especially in the context of the divide-and-conquer approach.909Table 3: Comparison between different models for fitting the entire data.910Model911R2912MAE913MAPE914RMSE915VAE916Computation time (in mins)917Simple hedonic9180.1669190.3279203.821%9210.4269220.0759232.19924Spatial hedonic9250.6059260.2189272.568%9280.2939290.0389304.92931Temporal9320.1739330.3259343.80%9350.4249360.07493725.23938SBD9390.7449400.1699411.984%9420.2369430.02794460945D&C-SBD9460.7399470.1719482.013%9490.2389500.0289511.50952D&C-STBD9530.7969540.1509551.767%9560.2119570.022958540959It is worth highlighting that, in an attempt to understand the efficacy of the divide-and-conquer960step, we tried to directly implement the model used in D&C-SBD. When comparing the performance961of the two approaches, we find that the method without the D&C step achieved only a marginally962better fit (R2 of 0.744, as compared to 0.739 in D&C-SBD). However, this improvement came at a963significant computational cost, with the spatial model taking approximately 40 times longer.964Due965to memory constraints, we could not perform the complete spatio-temporal model on the entire966dataset. Nevertheless, our findings indicate that the proposed model provides approximate results967closely aligned with what the actual model would produce.968We acknowledge that our proposed969model requires more time due to its incorporation of both spatial and temporal dimensions. However,970this additional computational effort is justified by the model’s ability to capture relationships more971effectively compared to other models.972To further support the above, it is important to recognize that the hedonic or spatial hedonic973model offer quite different results in terms of the coefficient estimates for some parameters.974For975instance, contrary to the common knowledge and existing literature, the hedonic model identifies a976significantly positive impact on the price if the property type is flat and these models fail to account97719978979980===== PAGE 20 (score 5, numeric-paren 1, table-words 0, var-words 4) =====981982for the significant premium associated with more spacious homes. These inconsistencies underscores983the intricate nature of housing dynamics and emphasizes the necessity for sophisticated analysis beyond984conventional linear models. As we explained in the previous subsection, suitable specification of spatial985and temporal dependence in our model helps in capturing the effect sizes in an appropriate way, and986renders valuable insights about the real estate markets. Furthermore, the proposed methodology offers987a superior balance between computational efficiency and accuracy, and can be a valuable option for988large datasets with spatial and temporal relationships.9895.3. Error parameters and residual diagnostics990We now focus on the estimated error structures in the proposed model and aim to diagnose the991residuals to assess whether the model has appropriately captured the dependence patterns in the992dataset. The estimated parameters of the Gaussian processes governing the residuals show that the993spatio-temporal error process explains more variation in the data than the white noise, as the estimated994σ2995v is almost double in magnitude compared to ˆσ2996ϵ . This results in a significant contribution of the997v(s, t) process to the overall variance and endorses the use of a spatio-temporal model to enhance998our understanding of London’s house price dynamics. Furthermore, we note that the estimated decay999parameters ˆϕs = 2.402 and ˆϕt = 0.528 reveal practical implications. With the equation exp(−ϕd) ≈10000.05 in mind, we can argue that the estimated values of these decay parameters imply that the spatial1001correlation remains significant for up to a distance of 1.25 kilometers, while the temporal dependence1002remain significant over a period of 5-6 months.1003We next look at a few diagnostic tests to assess the proposed model’s adequacy in capturing spatial1004and temporal dependence patterns. In the same spirit as in Section 2, we compute the Moran’s I1005statistics for the residuals at all time-points following (2.1). The plot is presented in Figure 6, and it1006evidently indicates that spatial correlations are no longer significant at 1% level. Contrasting it with1007Figure 3, we can ascertain that the spatial correlation has been successfully captured by our method.1008−0.0110090.0010100.011011010122510135010147510151001016Month1017Moran’s I Statistic1018not significant1019Fig. 6: Moran’s I statistic for the residuals for all the monthly time intervals in the entire London region. Black circle1020indicates a significant spatial correlation at 1% level of significance.1021For identifying any remaining temporal correlation, we analyze ACF for the residuals across different1022locations. We find that the residual series for majority of the MSOAs do not have significant temporal1023autocorrelation.1024For instance, the ACF plots of residuals for the same four MSOAs as before are1025illustrated in Figure 7, which signifies that the temporal correlation does not exist anymore.1026In summary, our proposed model effectively captures both spatial and temporal dependence patterns,1027providing a good fit for the data. In addition, we also looked at the histogram and the quantile-quantile1028(QQ) plot of the residuals, to assess the normality assumption. The analysis indicates that the residuals102920103010311032===== PAGE 23 (score 34, numeric-paren 0, table-words 3, var-words 2) =====10331034of forecasting for the response variable over all time-points in future. Finally, we present the individual1035MAPE values for the 23 MSOAs across the 12 time points in Table 5. Our analysis reveals that, barring1036three specific instances involving MSOAs and time points, we have accurately replicated the original1037log house prices, achieving a MAPE of less than 10%.1038310393.310403.610413.910424.210434.51044Nov 20181045Dec 20181046Jan 20191047Feb 20191048Mar 20191049Apr 20191050May 20191051Jun 20191052Jul 20191053Aug 20191054Sep 20191055Oct 20191056Future time−point1057MAPE1058310593.310603.610613.910624.210634.51064Nov 20181065Dec 20181066Jan 20191067Feb 20191068Mar 20191069Apr 20191070May 20191071Jun 20191072Jul 20191073Aug 20191074Sep 20191075Oct 20191076Future time−point1077MAPE1078Fig. 9: Prediction MAPE (in %) for 23 locations: left panel corresponds to the case where the locations are included in1079the training set, right panel corresponds to the case where they are not part of the training set.1080Table 5: MAPE values calculated for the 23 out-of-sample MSOAs across all 12 future time points, which were not part1081of the training dataset, employing the proposed model. Missing cells indicate that there was no sale in those months.1082More than 10% errors are marked in bold.1083MSOA1084Nov1085Dec1086Jan1087Feb1088Mar1089Apr1090May1091Jun1092Jul1093Aug1094Sep1095Oct1096201810972018109820191099201911002019110120191102201911032019110420191105201911062019110720191108Barnet 02711092.8811101.0111111.9311123.2511131.3811142.1111157.1411161.6911171.1311182.0211195.3011202.791121Barnet 03511221.1111232.4511241.0211251.9811262.7611273.1411283.8111294.6211302.7411312.1911320.6011331.881134Brent 00711351.9111363.0311370.6711382.6211392.2111401.8811414.2511424.5411432.7211442.4911453.5911465.011147Brent 00911481.6511492.8211501.9611511.8411529.0811532.0011542.0611552.9211562.1811576.711158Brent 01311594.6011606.5111614.8411628.2111631.4311643.3011651.6811664.2311675.5011683.1611693.281170Brent 02611712.8411720.3611738.1111740.3811751.5411765.2811772.1411782.4511792.0111800.8111815.151182Croydon 02111837.3311845.8711855.2111864.5411876.2311883.5111895.9111903.7011914.7311924.8911934.5311946.581195Ealing 03011961.4911974.0211980.7211992.6012000.1112012.0312020.8312031.5512041.0912051.6412062.3112072.051208Hackney 02412093.6412102.2212114.7012124.9212133.6912141.9912151.2112162.5912174.6712183.8512191.6312203.531221Hammersmith 02412224.8312235.0212243.0112254.2312263.9912274.6612283.3712297.0812305.9912314.7912323.3712334.011234Harrow 01912355.5612363.6112371.0212381.4012392.0012401.3412412.59124220.4712432.3512442.0812452.031246Harrow 02412472.3712482.5012492.2612502.9212514.2112524.3012533.8712541.8312554.1312561.7712572.7912585.951259Islington 00812604.6012616.7312624.5512632.8912642.1212651.8612665.7512673.9412683.0412693.7612700.0312714.741272Islington 01912732.5812743.2312754.0912761.5712772.1212781.7812792.4512803.4712813.0112821.8112834.771284Lambeth 02212853.1012860.8512872.3512883.9212892.6512901.6012910.9712920.9312932.2712942.5312952.9712961.431297Merton 01012982.6912992.7313001.8313012.1013023.5713031.6813040.7613052.8813062.2413072.9613082.2813091.931310Merton 01113111.5413121.9113132.3113142.2013152.0013161.5113172.1113181.6213191.6213202.2513211.6113221.871323Merton 01313243.4513253.8313263.4313275.9013284.6813293.8613304.3313313.3313321.4713332.7213344.9413353.881336Newham 03313372.1513382.9613394.5713402.6013413.0313422.4313433.2013440.821345Southwark 00213464.7313474.7013482.4413492.8913504.5713512.8413524.6113533.4013544.1213553.6013562.57135710.931358Sutton 01113593.3013604.0813613.9213623.9213634.5413643.8713655.5213665.8813674.3813684.8913694.1213703.911371Sutton 01913724.5113735.2613741.4213754.4113765.2113774.6813784.5213795.0313805.2813814.96138211.871383Westminster 01613847.1013857.2113868.2113875.0213887.3313894.9013907.0013916.4513925.8213935.2313946.2713956. Conclusion1396In this study, we have introduced a novel divide-and-conquer spatio-temporal modeling approach1397specially designed for large datasets. One distinctive aspect of our proposed methodology is its ability1398to accommodate multiple observations for the same location at the same time-points, as well as handling139923140014011402===== PAGE 24 (score 8, numeric-paren 0, table-words 1, var-words 6) =====14031404instances of zero transactions or missing data within the modeling framework. The model is applied to1405the house price data from London. The dataset comprises 106 monthly observations from 983 MSOAs.1406Through exploratory analysis, we identify significant spatial and temporal correlations, leading us1407to propose a spatio-temporal model. Built upon a Bayesian framework, our modeling process utilized1408Gibbs and discrete slice samplers for parameter estimation, allowing us to gain several key insights into1409the factors influencing house prices in London. Carbon emissions significantly influence property prices,1410reflecting the growing importance of sustainability. However, this influence is gradually diminishing,1411prompting a market shift. Supply and demand dynamics are adapting, with eco-friendly construction1412becoming more viable. Property characteristics, such as type and amenities, also affect prices. Our1413proposed model elucidates market complexity, emphasizing the need for adaptability to changing buyer1414preferences and sustainability trends in London’s real estate sector.1415Before concluding the paper, we want to suggest a few directions for future research. Exploratory1416data analysis revealed variations in the time trend effect on house prices among different MSOAs.1417An exciting avenue for future investigation involves exploring a more generalized modeling framework1418capable of accounting for regionally varying coefficients for covariates alongside the spatio-temporal1419error process. This extension could help identify specific regions experiencing significant slowdowns in1420growth rates or areas where people are willing to pay more for lower carbon emissions. In the same1421spirit, imposing a region-wise heteroskedastic structure on the white noise can improve the fit as well1422as predictive capability of the model. One may also relax the Gaussianity assumption and extend the1423model to heavy-tailed distributions.1424From an application standpoint, our method’s computational efficiency enhances its utility, making1425it an advantageous tool for a wide range of practical situations. The proposed approach extends beyond1426the confines of London’s housing market, becoming a versatile tool applicable to various other problems,1427thus providing value for stakeholders and policymakers across different domains. For instance, it can1428be effectively adapted to analyze space-time characteristics in datasets where high-level granularity is1429available with minimal missingness. A few examples can be found in environmental research (e.g., air1430pollution, rainfall analysis), epidemiology (e.g., COVID-19 data) or other economic applications (e.g.,1431agricultural production).1432Declaration of interest1433The authors declare no conflict of interest.1434References1435Arribas-Bel, D., 2014. Accidental, open and everywhere: Emerging data sources for the understanding1436of cities. Applied Geography 49, 45–53.1437Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014. Hierarchical Modeling and Analysis for Spatial Data.1438CRC Press.1439Beamonte, A., Gargallo, P., Salvador, M., 2010. Analysis of housing price by means of STAR models1440with neighbourhood effects: a Bayesian approach. Journal of Geographical Systems 12, 227–240.1441Bickel, P.J., Freedman, D.A., 1981. Some asymptotic theory for the bootstrap. The annals of statistics14429, 1196–1217.1443Blanco, H., Neri, L., 2023. Knocking It Down and Mixing It Up: The Impact of Public Housing1444Regenerations. Technical Report. IZA Discussion Papers.1445241446