# RECORD 549587ac03dcb9d550c7fe68da027277dd876b12
# TITLE: A divide-and-conquer approach for spatio-temporal analysis of large house price data from Greater London
# YEAR: 2024.0  VENUE: Annals of Operations Research  DOI: 10.1007/s10479-025-06889-x
# AUTHORS: Kapil Gupta; Soudeep Deb
# PAGES: 27  CANDIDATE TABLE PAGES (1-based): [8, 9, 11, 14, 17, 19, 23]
# 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.


===== PAGE 1 (score 8, numeric-paren 0, table-words 1, var-words 6) =====

A divide-and-conquer approach for spatio-temporal analysis of large house price
data from Greater London
Kapil Gupta, Soudeep Deb
Indian Institute of Management Bangalore
Bannerghatta Main Road, Bangalore, Karnataka 560076, India.
Abstract
Statistical research in real estate markets, particularly in understanding the spatio-temporal dynamics
of house prices, has garnered significant attention in recent times. Although Bayesian methods are
common in spatio-temporal modeling, standard Markov chain Monte Carlo (MCMC) techniques are
usually slow for large datasets such as house price data. To tackle this problem, we propose a divide-
and-conquer spatio-temporal modeling approach.
This method involves partitioning the data into
multiple subsets and applying an appropriate Gaussian process model to each subset in parallel. The
results from each subset are then combined using the Wasserstein barycenter technique to obtain the
global parameters for the original problem. The proposed methodology allows for multiple observations
per spatial and time unit, thereby offering added benefits for practitioners. 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. The methodology provides insightful findings about the
effects of various amenities, trend patterns, and the relationship between prices and carbon emissions.
Furthermore, as demonstrated through a cross-validation study, it shows good predictive accuracy
while balancing computational efficiency.
Keywords:
Bayesian analysis, Real estate, Space-time prediction, Wasserstein barycenter
1. Introduction
The valuation of residential properties is immensely important for various economic stakeholders,
including property dealers, sellers, and buyers. Traditionally, house price modeling has revolved around
the hedonic framework (Rosen, 1974).
Specifically, log-linear regression techniques have been at
the forefront in the analysis of house price transactions (Liu, 2013). However, the incorporation of
location attributes into the mean structure of hedonic regression, as recommended by Dubin and Sung
(1990), has shown limitations in capturing all significant neighborhood attributes, potentially giving
rise to spatial autocorrelation concerns. To overcome this challenge, researchers have advocated for the
inclusion of spatial effects into the error structure (Can, 1990; Pace and Gilley, 1998). Simultaneously,
the dynamics of house prices are profoundly influenced by temporal autocorrelation as well.
Pace
et al. (1998) leveraged the concepts of spatio-temporal modeling to comprehensively capture both
spatial and temporal effects within house price data.
This modeling approach has since become
a cornerstone in academic research. Building on this foundation, Liu (2013) harnessed the spatio-
temporal autoregressive (STAR) model to address correlated errors in traditional hedonic regression,
delving into spatial and temporal dependence considerations and further developing a house price
index as a pivotal part of their analysis. Interestingly, their methodology assumed the spatial weight
matrix to be lower triangular, which implies that two spatial units cannot mutually influence each
other. A different approach was taken by Fotheringham et al. (2003) who integrated geographically
weighted regression with time series forecasting techniques to capture the spatio-temporal variations
arXiv:2407.15905v1  [stat.AP]  22 Jul 2024


===== PAGE 8 (score 7, numeric-paren 4, table-words 0, var-words 3) =====

a separable structure for the spatial dependence and the temporal dependence. In other words, it is
written as a product of purely spatial and purely temporal covariance functions:
Cov {v(si, tu), v(sj, tv)} = σ2
vρs (∥si −sj∥; ϕs) ρt (∥tu −tv∥; ϕt) .
(3.4)
The separable covariance structure is favored as it significantly reduces the parameters in the
covariance matrix, and is a popular choice in similar spatio-temporal studies (see, e.g., Sahu et al.,
2010; Deb and Tsay, 2019). For both spatial and temporal covariance functions, we consider the Matérn
structure with the specific choice of ν = 0.5, resulting in the exponentially decaying pattern. Thus, we
have:
ρs (∥si −sj∥; ϕs) = exp (−ϕs ∥si −sj∥) ;
ρt (∥tu −tv∥; ϕt) = exp (−ϕt ∥tu −tv∥) ,
(3.5)
where ∥tu −tv∥represents the absolute time difference, and ∥si −sj∥is calculated via the Vincenty
Ellipsoid method, which computes the great circle distance on an ellipsoid. The latter is chosen for its
superior accuracy compared to other approaches (Hijmans, 2021). It is critical to note that, because
of the separability assumption imposed on the joint covariance function, we can write the covariance
matrix Σ = σ2
v (Σs ⊗Σt), where ⊗denotes the Kronecker product, and Σs, Σt respectively denote the
spatial and temporal correlation matrices, with
(Σs)ij = exp (−ϕs ∥si −sj∥) ;
(Σt)uv = exp (−ϕt ∥tu −tv∥) .
(3.6)
3.2. Prior Specifications
Let θ = (β⊤, σ2
v, σ2
ϵ , ϕs, ϕt)⊤be the vector of unknown parameters. We estimate them completely
within a Bayesian framework, primarily using Gibbs sampling (?), for which suitable prior specifications
are essential.
For all coefficients in β, we consider independent Gaussian priors with mean 0 and
variance c, where c is relatively large, typically on the order of 104. For the error variances σ2
ϵ and
σ2
v, we employ independent inverse gamma priors IG(a, λ), where a and λ are the shape and scale
parameters, respectively. We choose a = 2, making these priors non-informative. We also noticed that
the parameter estimates and model selection criteria are insensitive to the initial values of the scale
parameter, and throughout this study, we take λ = 1.
For estimating the decay parameters ϕ = (ϕs, ϕt) in the covariance functions, many prior studies
(e.g., Sahu et al., 2006; ?) have used cross-validation to find the optimal choices. Albeit this approach
is computationally more tractable, it has some limitations, as it explores a small set of values and
can miss viable solutions. An alternative is the single-variable slice sampling algorithm (Neal, 2003),
which provides direct estimation of the parameters within a Bayesian framework. While this method
thoroughly explores the parameter space, it demands extensive computation as the covariance matrix
must be recalculated and inverted in every iteration of the Bayesian computation.
Keeping the
advantages of both approaches in view, we aim to achieve a middle ground – somewhat like a discrete
slice sampling – wherein a discretized sample space is considered for the parameters (ϕs, ϕt), and the
parameters are estimated directly within the Bayesian setting. This approach allows us to provide
credible intervals for these parameters. In terms of the prior specifications, we adopt a discretized
uniform prior for ϕs, selecting values from 1 to 3 in intervals of 0.2, resulting in the set As. For ϕt,
similarly, we select equidistant points between 0.2 and 1 in intervals of 0.2, and denote it as At. Below,
the product As × At is indicated as Ast, while U indicates the mentioned uniform distribution.
As mentioned before, we use concepts of Gibbs sampling and discrete slice sampling to implement
the proposed model for estimation and prediction purposes. As the first step, following Section 3.1, we
can write the full model as
Y | θ, V ∼NN
Xβ + BV , σ2
ϵ IN


,
(3.7)
8


===== PAGE 9 (score 17, numeric-paren 3, table-words 4, var-words 6) =====

while the prior distributions for all parameters are specified as:
β ∼Nm+1 (0, cIm+1) ,
σ2
ϵ ∼IG(a, λ), σ2
v ∼IG(a, λ),
(ϕs, ϕt) ∼U(Ast),
V ∼NG
0, σ2
vΣs ⊗Σt


.
(3.8)
Using f(·) to denote the density function in a generic way, the complete posterior distribution is
f
β, V , σ2
ϵ , σ2
v, ϕs, ϕt | Y


∝f
Y | β, V , σ2
ϵ , σ2
v, ϕs, ϕt


f
V | σ2
v, ϕs, ϕt


f(σ2
v)f(ϕs)f(ϕt)f(σ2
ϵ )f(β).
(3.9)
For brevity, detailed calculations for the conditional posteriors in case of the complete data are
deferred to Appendix A.1 of the supplementary material.
In the next subsection, we present the
posterior distributions and calculations for the D&C setup.
3.3. Divide-and-conquer approach for large dataset
The above model and its Bayesian implementation are helpful in analyzing house price dynamics.
However, even with the assumption of separability and the use of Kronecker product to reduce the
cost of inverting dispersion matrices, this approach becomes computationally challenging for a dataset
with large number of locations or time-points. Specifically, each iteration in the above implementation
requires O(S2T 2) memory units for storage and O((S ∨T)3) floating-point operations for computation.
Naturally, this approach is infeasible for large-scale data, which motivates us to develop a divide-and-
conquer approach in this section building upon the above methodology. Although the use of divide-and-
conquer is not a new concept in computer science, its implementation in scalable Bayesian inference
is a recent advancement.
We present the technique in a general three-step distributed framework
applicable for scaling posterior computations in our spatio-temporal model for house prices presented
in Section 3.1.
3.3.1. First step: partitioning the data
Consider a large dataset with S spatial locations, T time-points, and multiple observations as set
out before. Our approach involves splitting these data points into Q non-overlapping subsets for model
feasibility. We use a random allocation scheme for the S locations into these Q subsets. This process
ensures that each subset contains an adequate number of locations and their corresponding data for
all time-points. To avoid bias from locations with significantly more or fewer observations, one should
take precautions during the random assignment, striving for balanced representation across the subsets.
Our specific approach for the main analysis is elaborated in Section 4.
To set the notations for this section, let Sq = (sq1, sq2, . . . sqmq ) denote the set of mq locations in
the qth subset, q ∈{1, . . . , Q}. Although in principle a spatial location can be included in multiple
subsets, in this work we only consider disjoint partition of the spatial index set, thereby implying
PQ
q=1 mq = S and ∪Q
q=1Sq being equal to the set of all locations. In the qth subset, we have the
data {Y(q), X(q)}, with the total number of observations being N(q) = Pmq
i=1
PT
j=1 kij. Clearly, Y(q)
is a N(q) × 1 vector, and X(q) is a N(q) × (m + 1) matrix representing the covariates corresponding
to the observations for all time-points in the spatial locations within the subset Sq. Following earlier
notations, the spatio-temporal Gaussian process model for the qth subset is:
Y(q) = X(q)β(q) + B(q)V(q) + ϵ(q).
(3.10)
Throughout the discussions in this section, the parameters θ(q) = (β(q), σ2
v(q), σ2
ϵ(q), ϕs(q), ϕt(q))⊤,
the data vector Y(q), the spatio-temporal process V(q), the matrices X(q), B(q), and all the available
9


===== PAGE 10 (score 5, numeric-paren 3, table-words 0, var-words 2) =====

information H(q) for subset q are equivalent to their full dataset counterparts θ = (β, σ2
v, σ2
ϵ , ϕs, ϕt)⊤,
Y , V , X, B, and H. Also, we use the notation pq = N/N(q) to simplify some expressions.
3.3.2. Second step: modified sampling from subset pseudo posterior distributions
We note that the qth subset described above contains a fraction (1/pq) of the full data. Thus,
using its posterior will result in an overestimation of the posterior uncertainty compared to using the
full dataset. To address this mismatch, we consider a modification to the Bayesian algorithm on the
subsets without sacrificing its efficiency. Particularly, a powered likelihood is applied to modify the
likelihood of all the parameters before employing the sampling algorithm on the qth subset. As it
contains 1/pq fraction of the full data, the asymptotic variances of the posterior distributions of θ(q)
and V(q) are inflated by a factor of pq compared to that of θ and V (Shyamalkumar and Srivastava,
2022). Subsequently, we employ stochastic approximation by raising the likelihood of the parameters
β(q), σ2
v(q), σ2
ϵ(q), ϕs(q), ϕt(q), V(q) in the qth subset to the power of pq (Minsker et al., 2014).
These
adjustments account for the data present in the other subsets, and ensure accurate estimation of
posterior uncertainty while maintaining computational efficiency.
For the qth subset, as there are N(q) data points, the covariance matrices Σs(q) and Σt(q) have
dimensions mq × mq and T × T, respectively. Let G(q) = mq × T, and K(q) be a diagonal matrix
of order G(q) × G(q), where the diagonal elements are kij for different (si, tj) combinations in the qth
subset. The prior distributions for all parameters in subset q align with those in (3.8), except for V(q),
which has the following structure:
V(q) ∼NG(q)

0, σ2
v(q)
Σs(q) ⊗Σt(q)


.
(3.11)
Utilizing the abovementioned stochastic approximation, we now obtain the following subset pseudo
posterior distribution. For the sake of brevity, we omit the technical derivations in the main text, and
refer to the supplementary material for detailed information.
β(q) | H(q) ∼Nm+1


 pqX⊤
(q)X(q)
σ2
ϵ(q)
+ Im+1
c
!−1  pqX⊤
(q)
Y(q) −B(q)V(q)


σ2
ϵ(q)
!  pqX⊤
(q)X(q)
σ2
ϵ(q)
+ Im+1
c
!−1
,
σ2
ϵ(q) | H(q) ∼IG
 
a + N
2 , pq


Y(q) −X(q)β(q) −V(q)



2
2
+ λ
!
,
σ2
v(q) | H(q) ∼IG

a + pqG(q)
2
,
pq

V ⊤
(q)

eΣ−1
s(q) ⊗Σ−1
t(q)

V(q)

2
+ λ

,
f(ϕs(q), ϕt(q)) | H(q) ∝|Σs(q)|−pqT/2|Σt(q)|−pqmq/2 exp


−pqV ⊤
(q)

eΣ−1
s(q) ⊗Σ−1
t(q))V(q)

σ2
v(q)

I
(ϕs(q), ϕt(q)) ∈Ast


,
V(q) | H(q) ∼NG


 eΣ−1
s(q) ⊗Σ−1
t(q)
σ2
v(q)
+ pqK(q)
σ2
ϵ(q)
!−1  pqB⊤
(q)(Y(q) −X(q)β(q))
σ2
ϵ(q)
!
,
 eΣ−1
s(q) ⊗Σ−1
t(q)
σ2
v(q)
+ pqK(q)
σ2
ϵ(q)
!−1
.
(3.12)
To maintain consistency in our analysis, we adopt identical potential values for decay parameters
across all subsets.
Therefore, we utilize the same combination of ϕs(q) and ϕt(q).
An essential
modification in our posterior calculation involves capturing spatial dependence between subsets. This is
achieved by adjusting the inverse of the subset covariance matrix to account for spatial interdependencies
10


===== PAGE 11 (score 9, numeric-paren 4, table-words 0, var-words 5) =====

among the subsets. To properly follow the procedure, one typically needs to invert the entire spatiotemporal
matrix (or the spatial dependence matrix under the separability assumption) before subsetting it for
the qth subset. However, inverting such a large matrix is very time-consuming. To address this, we
use the concept of Schur’s complement (Zhang, 2006) and find the inverses of the diagonal matrices
efficiently. To elaborate this approach, note that we first obtain the inverse of the last diagonal part
using Schur’s complement, which involves the inverse of a small partitioned matrix. The other diagonal
part is also the inverse of a matrix comprising only the inverse of a small matrix. In the next iteration,
we treat this inverse part as the inverse of the updated main matrix and repeat the procedure until
only one diagonal matrix remains. By employing this iterative method, we can find the inverses of
all the block diagonal matrices without needing to invert the entire matrix.
Our extension of the
Schur’s complement approach significantly speeds up the computation of the inverse. Let us denote
the inverse of Σs(q) as eΣ−1
s(q), reflecting the modified version after incorporating spatial dependence.
Detailed calculations of this procedure are shown in Appendix B of the supplementary material.
In our Gibbs sampling procedure, computing the posterior of the vector V(q) requires inverting a
G(q) × G(q) matrix in each iteration, which can be computationally demanding. To mitigate this, we
adopt a more efficient approach based on the theory of multivariate normal distribution. Specifically,
we partition the vector V(q) into mq subvectors based on the columns in the matrix V(q). Then, we
take Z
′
1 as the first subvector of V(q), and Z
′
2 as a concatenated subvector of the remaining columns
of V(q), i.e., V(q) =

Z
′
1 Z
′
2

. Similarly, the covariance matrix Σ(q) can be written as
Σ(q) = σ2
v(q)
Σ11(q) ⊗Σt(q) Σ12(q) ⊗Σt(q)
Σ21(q) ⊗Σt(q) Σ22(q) ⊗Σt(q)

.
(3.13)
Then, the conditional distribution of Z
′
1 given Z
′
2 = z′
2 is NT (µc(q), σ2
v(q)Σc(q)), where
µc(q) =
Σ12(q) ⊗Σt(q)

 
Σ−1
22(q) ⊗Σ−1
t(q)

z′
2;
Σc(q) =
Σ11(q) ⊗Σt(q)


−
Σ12(q) ⊗Σt(q)

 
Σ−1
22(q) ⊗Σ−1
t(q)
 Σ21(q) ⊗Σt(q)


.
(3.14)
Using the properties of the Kronecker product, the expression can be further simplified as
µc(q) =

Σ12(q)Σ−1
22(q)

⊗IT

z′
2;
Σc(q) =

Σ11(q) −Σ12(q)Σ−1
22(q)Σ21(q)

⊗Σt(q).
(3.15)
Now, following the above, we can write
V1(q) | H(q) ∼NT


 Σ−1
c(q)
σ2
v(q)
+ pqK1(q)
σ2
ϵ(q)
!−1 
Σ−1
c(q)µc(q)
σ2
v(q)
+
pq

BT
1(q)
Y1(q) −X1(q)β(q)


σ2
ϵ(q)

,
 Σ−1
c(q)
σ2
v(q)
+ pqK1(q)
σ2
ϵ(q)
!−1
,
(3.16)
where Y1(q) be the data for the first location and all time-points from the qth subset. The matrix X1(q)
is defined accordingly. K1(q) is a T × T diagonal matrix with each diagonal element being equal to
k1j for each (sq1, tj). B1(q) is a submatrix of B(q), consisting of the first T columns and PT
j=1 k1j
rows, resulting in a sparse matrix of size PT
j=1 k1j × T. Other notations have the same definitions
as before. We can compute the conditional posterior distributions of V2(q), . . . , Vmq(q) in an identical
fashion. Finally, to compute the posterior distribution of the full vector V(q), we iteratively employ the
11


===== PAGE 12 (score 6, numeric-paren 2, table-words 0, var-words 4) =====

conditional distributions in a Gibbs sampler. This approach generates independent posterior samples
after an initial burn-in period. Consequently, the posterior distribution of V(q) is computed efficiently
without incurring the computational cost of inverting the G(q) × G(q) matrix in each Gibbs sampler
iteration. It is imperative to point out that these subset pseudo posteriors significantly speed up the
computations and can be run in parallel for all partitions.
3.3.3. Third step: combination of subset posteriors
The final and most crucial step of the D&C algorithm is to combine the subset posterior draws.
As discussed in Section 1, various existing methods can be employed in this regard. We follow the
Wasserstein Barycenter (WB) technique to approximate the complete posterior distribution.
This
approach has been successfully utilized by Li et al. (2017) and Srivastava et al. (2018) in the setting
of independent data, whereas Ou et al. (2021) demonstrated the accuracy of using this for combining
subset posteriors in the case of time series data. Motivated by these studies, we apply a similar rationale
for spatio-temporal data and utilize the WB to combine the subset posteriors. It is worth noting that
Wasserstein-based posteriors exhibit appealing asymptotic properties, as shown by Szabó et al. (2019).
Before going into the detail, it is critical to discuss the mathematical underpinning of the Wasserstein
distance.
Let (M, d) be a complete separable metric space, P(M) be the space of all probability
measures on M and P2(M) be the set of all probability measures on M with finite second moments.
For ϱ, ξ ∈P2(M), Π(ϱ, ξ) denotes the set of all probability measures on M × M with marginals ϱ and
ξ. Then, with d(·, ·) indicating the euclidean metric, the Wasserstein distance of order 2 between ϱ and
ξ is defined as
W2(ϱ, ξ) =

inf
π∈Π(ϱ,ξ)
Z
M×M
d(x1, x2)2Π(x1, x2)
1/2
.
(3.17)
As Bickel and Freedman (1981) discussed, the convergence in W2 on P2(M) implies weak convergence
of probability measures and convergence of the second moment. Keeping that in view, if θ(1), . . . , θ(Q)
are the Q subset posterior distributions, then the approximated posterior ˆθ is defined as
ˆθ = arg min
θ∈P2(M)
1
Q
Q
X
i=1
W 2
2 (θ, θ(i)).
(3.18)
This approach is referred to as the Wasserstein Posterior (WASP). The exact computation of the
WASP is known to be computationally intensive and remains a subject of ongoing research. We adopt
an approximation for ˆθ, following Shyamalkumar and Srivastava (2022). This completes the last step
in our proposed approach (D&C-STBD), and the pseudocode of the entire procedure is shown in
Algorithm 1.
A brief outline of the WASP procedure is also of the essence here.
First, posterior draws of
θ(1), . . . , θ(Q) are taken from the Q subsets using the Bayesian algorithm developed earlier.
In all
instances of Gibbs and discrete slice samplers, the convergence of the chains are confirmed through the
diagnostic method of Geweke et al. (1991). After convergence, we apply thinning to collect posterior
samples for each subset. These sample points (say, L for each subset) are then consolidated into a
unified set of QL samples. Using these posterior samples, we compute mean vectors ˆµ(q) and covariance
matrices ˆΣ(q) for each subset. This results in a total of Q mean vectors and covariance matrices. Next,
we take the mean of these and get the overall mean ˆ¯µ and covariance matrix ˆ¯Σ. Now, using the subset
mean vectors and covariance matrices, we standardize and center the subset posterior samples. Then,
we use the overall mean vector and covariance matrix to get the approximate WASP draws for the
actual posterior draws, which leads to QL number of WASP posterior samples for each parameter. We
estimate both the posterior mean of these parameters and their corresponding credible intervals using
the consolidated sample. On the other hand, for the final estimate of V from the subset posterior
12


===== PAGE 14 (score 11, numeric-paren 8, table-words 0, var-words 3) =====

follows a Gaussian distribution:
y(s′, t′, r′) | v(s′, t′) ∼N

X⊤(s′, t′, r′) ˆβ + v(s′, t′), σ2
ϵ

,
(3.19)
We estimate the model parameters using a Gibbs sampler on the training data, yielding posterior
estimates for β and σ2
ϵ . To find the conditional distribution of v(s′, t′) given V , we use the posterior
value of V from training data estimates and apply the multivariate Gaussian theory related to the
conditional distribution. We treat v(s′, t′) as the first subvector of the multivariate joint distribution,
and V obtained from the training dataset as the second subvector. The joint distribution of v(s′, t′)
and V is given by
v(s′, t′)
V

∼NG+1
0
0

, σ2
v
 1
Σ′⊤
Σ′
Σs ⊗Σt

, with Σ′ = Σs(s −s′) ⊗Σt(t −t′),
(3.20)
where Σ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 order
T with jth entry as ρt (||tj −t′||; ϕt). The conditional distribution is then given by
v(s′, t′) | V ∼N

Σ′(Σ−1
s
⊗Σ−1
t )V , σ2
v(1 −Σ′(Σ−1
s
⊗Σ−1
t )Σ′⊤)

.
(3.21)
Thus, by predicting v(s′, t′) using the above conditional distribution and the estimates of β from
the Gibbs sampling procedure, we can predict y(s′, t′, r′) using (3.19).
Now, we shift our focus to the prediction stage in D&C case.
Keeping parity with the earlier
notations, we need to forecast v(s′, t′) for a location s′ and a time-point t′ in the test dataset. Our
training dataset is divided into Q subsets, each treated as a separate training dataset. For each subset,
we predict v(s′, t′) using a similar process as for the whole dataset, resulting in Q realizations of v(s′, t′).
These realizations are denoted as vq(s′, t′), for 1 ⩽q ⩽Q. The joint distribution of vq(s′, t′) and V(q)
is given by
vq(s′, t′)
V(q)

∼NG(q)+1
 0
0

, σv2
(q)
"
1
Σ′⊤
(q)
Σ′
(q)
Σs(q) ⊗Σt(q)
#!
, with Σ′
(q) = Σs(q)(s −s′) ⊗Σt(q)(t −t′),
(3.22)
where Σs(q)(s−s′) is a vector of order S(q) with ith entry as ρs
||si −s′||; ϕs(q)


and Σt(t−t′)is a vector
of order T with jth entry as ρt
||tj −t′||; ϕt(q)


. The conditional distribution is then given by
vq(s′, t′) | V(q) ∼N

Σ′
(q)(eΣ−1
s(q) ⊗Σ−1
t(q))V(q), σ2
v(q)(1 −Σ′
(q)(eΣ−1
s(q) ⊗Σ−1
t(q))Σ′⊤
(q))

.
(3.23)
With the above predicted values of v(q)(s′, t′) and the estimated parameters within each subset q,
we can predict the response variable, i.e. the price of a house in logarithmic scale, for any MSOA at a
future time-point in the test data using (3.24), which is based on the qth subset data.
y(q)(s′, t′, r′) | v(q)(s′, t′) ∼N

X⊤
(q)(s′, t′, r′) ˆβ(q) + v(q)(s′, t′), σϵ2
(q)

.
(3.24)
Finally, to predict v(s′, t′) and the response variable y(s′, t′, r′) using the entire training dataset, we
take the median of vq(s′, t′) and y(q)(s′, t′, r′) respectively, for 1 ⩽q ⩽Q. This approach leverages the
D&C framework to efficiently handle large datasets by breaking them into smaller, manageable subsets
while maintaining prediction accuracy. It is noteworthy that we incorporate the original inverse of the
spatial covariance matrix for each subset to address spatial dependence.
14


===== PAGE 15 (score 8, numeric-paren 1, table-words 3, var-words 1) =====

4. Model implementation and evaluation
In this section, we explore the key implementation steps for our model and discuss other competing
models along with the evaluation criteria used for benchmarking.
Recall that our dataset comprises data on N = 651202 transactions from S = 983 locations and
T = 106 time-points. As covariates, we include the variables listed in Table 1, along with linear and
quadratic time trends, as well as the interaction of time with carbon emissions. We take logarithmic
transformation for all the numeric variables.
In the first step of our methodology, the dataset is
partitioned into Q = 20 disjoint subsets. The initial 19 subsets are comprised of randomly selected
49 locations, while the last subset contains the remaining 52 locations. Each subset encompasses data
for all the time-points and their corresponding multiple records. Next, for each subset, we employ
the model specified in (3.10) and fit the models separately using the posterior calculations shown
there.
After convergence of the chains, 2000 MCMC samples are taken in each of the Q subsets.
Finally, we adopt the subset posterior combination method, as described in Section 3.3.3, to obtain
final estimates for all the parameters. Additionally, we determine the values of v(s, t) corresponding
to each combination of (s, t) in the dataset.
To assess the fit of the proposed approach, we conduct a comparative analysis with four alternative
models. The first model is a standard hedonic regression that employs the same covariates as our
approach. As discussed in Section 1, it is one of the most popular approaches in related research.
The second model incorporates an additional fixed effect for the geographic region where the MSOA
is situated. We refer to this model as the spatial hedonic model, and it closely parallels the approach
described in ?. While these two approaches are conventional techniques in the extant literature of house
price analysis, we include two other models in the comparison study for understanding the effectiveness
of the proposed spatial and temporal dependence structures. The third model incorporates temporal
dependence in a manner similar to the main model proposed in this paper, but without the spatial
component. Given that there are 106 time points, the inversion of a 106 × 106 matrix during iteration
is computationally efficient, making this model’s computation manageable without the D&C step. The
fourth model, on the other hand, does not consider temporal dependence but introduces only a spatial
error process. With 983 locations, matrix inversion in this model could be time-consuming in each
iteration. Therefore, this spatial error process model is integrated into the proposed D&C framework
and implemented in the proposed fashion. This model will be referred to as D&C-SBD. By comparing
our method with these models, we aim to assess their respective data-fitting capabilities. Specifically,
we evaluate the accuracy of each model’s fit to the data and examine whether the resulting estimates are
consistent with existing literature. This comparison will help determine the robustness and reliability
of our approach in relation to established models.
We evaluate the efficacy of these models by fitting the data. Label the entire dataset as D, with
an observation for location si at time-point tk be denoted as y(si, tk, rik). The corresponding fitted
value is ˆy(si, tk, rik). To understand the goodness of fit of every model, we start with the ubiquitous
coefficient of determination (R2) which quantifies the amount of variability in the response variable
that can be accounted for by the predictor variables, thereby indicating the degree of adequacy of
the model.
Further, we consider the mean absolute error (MAE), mean absolute percentage error
(MAPE) and root mean square error (RMSE) of the fitted values, as recommended in statistical
literature (Botchkarev, 2019). The MAE, i.e. the average absolute difference between actual and fitted
values, provides a sense of the extent of the errors without considering their direction (overestimation
or underestimation).
MAPE, on the other hand, computes the percentage deviation and offers a
straightforward interpretation of the relative error. RMSE is one of the most popular metrics and
is advantageous because it penalises large errors more severely than small errors, thereby providing a
more impartial evaluation of the models’ accuracy. Finally, in some cases, a model with a lower average
error (as measured by MAE, MAPE, or RMSE) may exhibit higher variability in errors. To that end,
15


===== PAGE 17 (score 21, numeric-paren 0, table-words 3, var-words 9) =====

Table 2: Parameter estimates (posterior mean) and the corresponding credible intervals (CI) when the proposed model
is fitted to the whole data with 983 MSOAs and 106 time-points (651202 observations).
Variable
Estimate
95% CI
Intercept
9.675
(9.640 , 9.711)
Linear time trend
0.982
(0.845 , 1.103)
Quadratic time trend
−0.646
(−0.757, −0.534)
log (1 + Area)
−0.319
(−0.321, −0.316)
Rooms (3 Rooms)
0.066
(0.065 , 0.068)
Rooms (4 Rooms)
0.059
(0.057 , 0.062)
Rooms (5 Rooms)
0.064
(0.061 , 0.067)
Rooms (> 5 Rooms)
0.109
(0.106 , 0.113)
Property type (Terraced)
−0.052
(−0.054, −0.050)
Property type (Detached)
0.141
(0.138, 0.144)
Property type (Flats)
−0.250
(−0.252, −0.247)
Fireplace
0.064
(0.062, 0.065)
Ventilation (Extract-only)
0.029
(0.022 , 0.036)
Ventilation (Both supply and extract)
0.055
(0.043 , 0.065)
Wind turbine
0.006
(−0.019 , 0.031)
log (1 + CO2 emission)
−0.031
(−0.035, −0.028)
Interaction of trend and CO2 emission
0.076
(0.070 , 0.081)
σ2
ϵ
0.043
(0.043 , 0.043 )
σ2
v
0.083
(0.070 , 0.103)
ϕs
2.402
(1.486 , 3.041)
ϕt
0.528
(0.183 , 0.824)
9.0
9.5
10.0
10.5
11.0
Jan 2011
Jul 2011
Feb 2012
Sep 2012
Apr 2013
Nov 2013
Jun 2014
Jan 2015
Aug 2015
Mar 2016
Oct 2016
May 2017
Dec 2017
Jul 2018
Feb 2019
Oct 2019
Month
log price
Fig. 5: Average trend of the house price (log-transformed) for all locations over 106 time.
5.1.1. Effect of area and number of rooms
Our analysis unveils a straightforward relationship between property size and price in London. The
negative coefficient associated with log area indicates that as a property’s size increases, its value per
17


===== PAGE 18 (score 17, numeric-paren 0, table-words 2, var-words 13) =====

square meter tends to decrease. This observation is in line with common expectations and matches
the results reported in extant literature, e.g., Chegut et al. (2016).
Furthermore, the number of rooms within a property emerges as a significant factor. Properties
with more than two rooms command notably higher prices compared to those with 1 or 2 rooms. This
effect becomes more pronounced for properties with over five rooms, where the impact on price per
square meter is most substantial. It highlights the premium placed on larger, more spacious homes in
the London region, reflecting the diverse housing preferences of potential buyers.
In summary, larger properties with same number of rooms tend to have lower prices per square
meter than smaller properties. This finding aligns with the results reported in Fuerst et al. (2015).
5.1.2. Effect of property characteristics
The type of property stands out as a pivotal factor influencing prices. We find that terraced and flat
properties tend to have lower prices per square meter compared to semi-detached houses. Strikingly,
flats are at the lower end of the price spectrum, while detached houses command the highest prices.
This finding aligns with the study conducted by Fuerst et al. (2015), where terraced properties were
considered the base category, and both detached and semi-detached properties were found to be on the
premium side, while flats were on the lower end. This result underscores the importance of property
type as a driver of price volatility in the London real estate market and reflects the diverse architectural
landscape of the region.
The presence of fireplaces within properties is strongly associated with higher prices. It aligns
with the perception of fireplaces as desirable amenities for homeowners, contributing positively to
property values, and is consistent with the observations of Zhang (2016). Meanwhile, we find that
the ventilation systems exert a discernible influence on house prices. Properties equipped with extract
mechanical ventilation and both supply and extract mechanical ventilation systems tend to exhibit
higher prices compared to those reliant on natural ventilation. One may connect this to the growing
emphasis on indoor air quality and environmental factors among property buyers. It reinforces the
desirability of properties that offer superior air quality, facilitated by mechanical ventilation systems.
5.1.3. Effect of carbon emission
Another intriguing aspect of our analysis is the investigation of the relationship between carbon
emissions and house prices.
The increasing need of a cleaner environment is likely to impact the
economy of real estate and we attempt to capture the evolving dynamics of sustainability in the
London housing market through this analysis. By considering appropriate regressors, we delve into
the multifaceted aspects of this relationship, including its impact, interaction with time, and the
implications for supply and demand.
Impact of carbon emission on house price: Our findings highlight the significant impact of
carbon emissions on house prices. We observe a negative relationship, indicating that properties with
lower carbon emission values tend to command higher prices. This emphasizes the growing importance
of sustainability and environmental consciousness among real estate buyers. It suggests that buyers are
increasingly willing to invest more in homes with a reduced carbon footprint, reflecting the evolving
trend toward eco-friendly living spaces. The findings align with those of Gerassimenko et al. (2023),
who noted that energy-efficient properties, characterized by low carbon emissions, tend to command a
price premium.
The interaction with time: A more nuanced observation emerges when we consider the interaction
between carbon emissions and time. Over time, we observe a decline in willingness to pay premium
for lower carbon emissions, reflecting a shift in buyer behavior. This observation is supported by the
positive coefficient associated with the interaction term. In practice, while carbon emissions continue
to be an important factor for buyers, their impact on house prices has gradually diminished over the
18


===== PAGE 19 (score 16, numeric-paren 0, table-words 5, var-words 3) =====

years. This trend raises interesting questions about the evolution of buyer preferences and the broader
real estate market.
Supply and demand implications: Finally, the interplay between carbon emissions and house
prices also has implications for supply and demand dynamics. As sustainability becomes a central
consideration for buyers, properties with lower carbon emissions are likely to see increased demand.
This heightened demand may prompt the builders to focus on constructing more environmentally
friendly properties, further influencing the supply side of the equation. However, our observation that
the price differential between eco-friendly and conventional properties is diminishing suggests that the
market is finding a balance. As supply increases and the premium on sustainability lessens, a broader
range of buyers can access these eco-friendly properties (see Ma and Li, 2017; Blanco and Neri, 2023,
for pertinent discussions). This equilibrium benefits both buyers and the environment.
5.2. Comparison between the fits of different models
As mentioned in Section 4, we conduct a comprehensive comparative study to assess the fitting
and accuracy of different models. Table 3 displays the accuracy metrics on overall data to evaluate
the goodness of fit. The results unequivocally demonstrate that our proposed model outperforms the
other models, indicating its superior ability to capture the spatial and temporal dependencies in the
data and achieve a better fit. Notably, D&C-STBD explains the highest variation, while the temporal
model falls short in explaining the variance effectively. We also observe that the D&C-SBD method
requires significantly less time compared to the spatio-temporal model or the spatial hedonic model,
yet it provides a great fit. This highlights the efficiency and accuracy of the spatial component in the
method, especially in the context of the divide-and-conquer approach.
Table 3: Comparison between different models for fitting the entire data.
Model
R2
MAE
MAPE
RMSE
VAE
Computation time (in mins)
Simple hedonic
0.166
0.327
3.821%
0.426
0.075
2.19
Spatial hedonic
0.605
0.218
2.568%
0.293
0.038
4.92
Temporal
0.173
0.325
3.80%
0.424
0.074
25.23
SBD
0.744
0.169
1.984%
0.236
0.027
60
D&C-SBD
0.739
0.171
2.013%
0.238
0.028
1.50
D&C-STBD
0.796
0.150
1.767%
0.211
0.022
540
It is worth highlighting that, in an attempt to understand the efficacy of the divide-and-conquer
step, we tried to directly implement the model used in D&C-SBD. When comparing the performance
of the two approaches, we find that the method without the D&C step achieved only a marginally
better fit (R2 of 0.744, as compared to 0.739 in D&C-SBD). However, this improvement came at a
significant computational cost, with the spatial model taking approximately 40 times longer.
Due
to memory constraints, we could not perform the complete spatio-temporal model on the entire
dataset. Nevertheless, our findings indicate that the proposed model provides approximate results
closely aligned with what the actual model would produce.
We acknowledge that our proposed
model requires more time due to its incorporation of both spatial and temporal dimensions. However,
this additional computational effort is justified by the model’s ability to capture relationships more
effectively compared to other models.
To further support the above, it is important to recognize that the hedonic or spatial hedonic
model offer quite different results in terms of the coefficient estimates for some parameters.
For
instance, contrary to the common knowledge and existing literature, the hedonic model identifies a
significantly positive impact on the price if the property type is flat and these models fail to account
19


===== PAGE 20 (score 5, numeric-paren 1, table-words 0, var-words 4) =====

for the significant premium associated with more spacious homes. These inconsistencies underscores
the intricate nature of housing dynamics and emphasizes the necessity for sophisticated analysis beyond
conventional linear models. As we explained in the previous subsection, suitable specification of spatial
and temporal dependence in our model helps in capturing the effect sizes in an appropriate way, and
renders valuable insights about the real estate markets. Furthermore, the proposed methodology offers
a superior balance between computational efficiency and accuracy, and can be a valuable option for
large datasets with spatial and temporal relationships.
5.3. Error parameters and residual diagnostics
We now focus on the estimated error structures in the proposed model and aim to diagnose the
residuals to assess whether the model has appropriately captured the dependence patterns in the
dataset. The estimated parameters of the Gaussian processes governing the residuals show that the
spatio-temporal error process explains more variation in the data than the white noise, as the estimated
σ2
v is almost double in magnitude compared to ˆσ2
ϵ . This results in a significant contribution of the
v(s, t) process to the overall variance and endorses the use of a spatio-temporal model to enhance
our understanding of London’s house price dynamics. Furthermore, we note that the estimated decay
parameters ˆϕs = 2.402 and ˆϕt = 0.528 reveal practical implications. With the equation exp(−ϕd) ≈
0.05 in mind, we can argue that the estimated values of these decay parameters imply that the spatial
correlation remains significant for up to a distance of 1.25 kilometers, while the temporal dependence
remain significant over a period of 5-6 months.
We next look at a few diagnostic tests to assess the proposed model’s adequacy in capturing spatial
and temporal dependence patterns. In the same spirit as in Section 2, we compute the Moran’s I
statistics for the residuals at all time-points following (2.1). The plot is presented in Figure 6, and it
evidently indicates that spatial correlations are no longer significant at 1% level. Contrasting it with
Figure 3, we can ascertain that the spatial correlation has been successfully captured by our method.
−0.01
0.00
0.01
0
25
50
75
100
Month
Moran’s I Statistic
not significant
Fig. 6: Moran’s I statistic for the residuals for all the monthly time intervals in the entire London region. Black circle
indicates a significant spatial correlation at 1% level of significance.
For identifying any remaining temporal correlation, we analyze ACF for the residuals across different
locations. We find that the residual series for majority of the MSOAs do not have significant temporal
autocorrelation.
For instance, the ACF plots of residuals for the same four MSOAs as before are
illustrated in Figure 7, which signifies that the temporal correlation does not exist anymore.
In summary, our proposed model effectively captures both spatial and temporal dependence patterns,
providing a good fit for the data. In addition, we also looked at the histogram and the quantile-quantile
(QQ) plot of the residuals, to assess the normality assumption. The analysis indicates that the residuals
20


===== PAGE 23 (score 34, numeric-paren 0, table-words 3, var-words 2) =====

of forecasting for the response variable over all time-points in future. Finally, we present the individual
MAPE values for the 23 MSOAs across the 12 time points in Table 5. Our analysis reveals that, barring
three specific instances involving MSOAs and time points, we have accurately replicated the original
log house prices, achieving a MAPE of less than 10%.
3
3.3
3.6
3.9
4.2
4.5
Nov 2018
Dec 2018
Jan 2019
Feb 2019
Mar 2019
Apr 2019
May 2019
Jun 2019
Jul 2019
Aug 2019
Sep 2019
Oct 2019
Future time−point
MAPE
3
3.3
3.6
3.9
4.2
4.5
Nov 2018
Dec 2018
Jan 2019
Feb 2019
Mar 2019
Apr 2019
May 2019
Jun 2019
Jul 2019
Aug 2019
Sep 2019
Oct 2019
Future time−point
MAPE
Fig. 9: Prediction MAPE (in %) for 23 locations: left panel corresponds to the case where the locations are included in
the training set, right panel corresponds to the case where they are not part of the training set.
Table 5: MAPE values calculated for the 23 out-of-sample MSOAs across all 12 future time points, which were not part
of the training dataset, employing the proposed model. Missing cells indicate that there was no sale in those months.
More than 10% errors are marked in bold.
MSOA
Nov
Dec
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
2018
2018
2019
2019
2019
2019
2019
2019
2019
2019
2019
2019
Barnet 027
2.88
1.01
1.93
3.25
1.38
2.11
7.14
1.69
1.13
2.02
5.30
2.79
Barnet 035
1.11
2.45
1.02
1.98
2.76
3.14
3.81
4.62
2.74
2.19
0.60
1.88
Brent 007
1.91
3.03
0.67
2.62
2.21
1.88
4.25
4.54
2.72
2.49
3.59
5.01
Brent 009
1.65
2.82
1.96
1.84
9.08
2.00
2.06
2.92
2.18
6.71
Brent 013
4.60
6.51
4.84
8.21
1.43
3.30
1.68
4.23
5.50
3.16
3.28
Brent 026
2.84
0.36
8.11
0.38
1.54
5.28
2.14
2.45
2.01
0.81
5.15
Croydon 021
7.33
5.87
5.21
4.54
6.23
3.51
5.91
3.70
4.73
4.89
4.53
6.58
Ealing 030
1.49
4.02
0.72
2.60
0.11
2.03
0.83
1.55
1.09
1.64
2.31
2.05
Hackney 024
3.64
2.22
4.70
4.92
3.69
1.99
1.21
2.59
4.67
3.85
1.63
3.53
Hammersmith 024
4.83
5.02
3.01
4.23
3.99
4.66
3.37
7.08
5.99
4.79
3.37
4.01
Harrow 019
5.56
3.61
1.02
1.40
2.00
1.34
2.59
20.47
2.35
2.08
2.03
Harrow 024
2.37
2.50
2.26
2.92
4.21
4.30
3.87
1.83
4.13
1.77
2.79
5.95
Islington 008
4.60
6.73
4.55
2.89
2.12
1.86
5.75
3.94
3.04
3.76
0.03
4.74
Islington 019
2.58
3.23
4.09
1.57
2.12
1.78
2.45
3.47
3.01
1.81
4.77
Lambeth 022
3.10
0.85
2.35
3.92
2.65
1.60
0.97
0.93
2.27
2.53
2.97
1.43
Merton 010
2.69
2.73
1.83
2.10
3.57
1.68
0.76
2.88
2.24
2.96
2.28
1.93
Merton 011
1.54
1.91
2.31
2.20
2.00
1.51
2.11
1.62
1.62
2.25
1.61
1.87
Merton 013
3.45
3.83
3.43
5.90
4.68
3.86
4.33
3.33
1.47
2.72
4.94
3.88
Newham 033
2.15
2.96
4.57
2.60
3.03
2.43
3.20
0.82
Southwark 002
4.73
4.70
2.44
2.89
4.57
2.84
4.61
3.40
4.12
3.60
2.57
10.93
Sutton 011
3.30
4.08
3.92
3.92
4.54
3.87
5.52
5.88
4.38
4.89
4.12
3.91
Sutton 019
4.51
5.26
1.42
4.41
5.21
4.68
4.52
5.03
5.28
4.96
11.87
Westminster 016
7.10
7.21
8.21
5.02
7.33
4.90
7.00
6.45
5.82
5.23
6.27
6. Conclusion
In this study, we have introduced a novel divide-and-conquer spatio-temporal modeling approach
specially designed for large datasets. One distinctive aspect of our proposed methodology is its ability
to accommodate multiple observations for the same location at the same time-points, as well as handling
23


===== PAGE 24 (score 8, numeric-paren 0, table-words 1, var-words 6) =====

instances of zero transactions or missing data within the modeling framework. The model is applied to
the house price data from London. The dataset comprises 106 monthly observations from 983 MSOAs.
Through exploratory analysis, we identify significant spatial and temporal correlations, leading us
to propose a spatio-temporal model. Built upon a Bayesian framework, our modeling process utilized
Gibbs and discrete slice samplers for parameter estimation, allowing us to gain several key insights into
the factors influencing house prices in London. Carbon emissions significantly influence property prices,
reflecting the growing importance of sustainability. However, this influence is gradually diminishing,
prompting a market shift. Supply and demand dynamics are adapting, with eco-friendly construction
becoming more viable. Property characteristics, such as type and amenities, also affect prices. Our
proposed model elucidates market complexity, emphasizing the need for adaptability to changing buyer
preferences and sustainability trends in London’s real estate sector.
Before concluding the paper, we want to suggest a few directions for future research. Exploratory
data analysis revealed variations in the time trend effect on house prices among different MSOAs.
An exciting avenue for future investigation involves exploring a more generalized modeling framework
capable of accounting for regionally varying coefficients for covariates alongside the spatio-temporal
error process. This extension could help identify specific regions experiencing significant slowdowns in
growth rates or areas where people are willing to pay more for lower carbon emissions. In the same
spirit, imposing a region-wise heteroskedastic structure on the white noise can improve the fit as well
as predictive capability of the model. One may also relax the Gaussianity assumption and extend the
model to heavy-tailed distributions.
From an application standpoint, our method’s computational efficiency enhances its utility, making
it an advantageous tool for a wide range of practical situations. The proposed approach extends beyond
the confines of London’s housing market, becoming a versatile tool applicable to various other problems,
thus providing value for stakeholders and policymakers across different domains. For instance, it can
be effectively adapted to analyze space-time characteristics in datasets where high-level granularity is
available with minimal missingness. A few examples can be found in environmental research (e.g., air
pollution, rainfall analysis), epidemiology (e.g., COVID-19 data) or other economic applications (e.g.,
agricultural production).
Declaration of interest
The authors declare no conflict of interest.
References
Arribas-Bel, D., 2014. Accidental, open and everywhere: Emerging data sources for the understanding
of cities. Applied Geography 49, 45–53.
Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014. Hierarchical Modeling and Analysis for Spatial Data.
CRC Press.
Beamonte, A., Gargallo, P., Salvador, M., 2010. Analysis of housing price by means of STAR models
with neighbourhood effects: a Bayesian approach. Journal of Geographical Systems 12, 227–240.
Bickel, P.J., Freedman, D.A., 1981. Some asymptotic theory for the bootstrap. The annals of statistics
9, 1196–1217.
Blanco, H., Neri, L., 2023. Knocking It Down and Mixing It Up: The Impact of Public Housing
Regenerations. Technical Report. IZA Discussion Papers.
24