# 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 in house prices. This was later
improved with spatio-temporal data mining techniques by Soltani et al. [2021].
Holly et al. [2010], on the other hand, developed a spatio-temporal
econometric model for analyzing house price dynamics in the United States, by
explicitly taking into account both cross-sectional dependence and
heterogeneity. We also find it prudent to refer to the early review of spatial
and spatio-temporal models within housing literature by Gelfand et al. [2004].
The authors scrutinized several hierarchical models and briefly explored their
Bayesian implementation to demonstrate the effectiveness of a Bayesian
framework in this domain. Intriguingly, except for few studies (e.g., Bayesian
implementation of STAR by Beamonte et al. [2010], Bayesian network approach by
Teye and Ahelegbey [2017]), Bayesian spatio-temporal models have not been
utilized appropriately in related statistical analysis. We attempt to bridge
this gap in this article by proposing an effective Bayesian spatio-temporal
approach in conjunction with a divide-and-conquer technique, to analyze large
house price datasets.

It is of the essence here to briefly review a few motivating literature. The
advent of the big data era has ushered in a new set of challenges for
statistical inference, particularly in the realm of current topic of
discussion. Complex Gaussian process models, though flexible, become
computationally unwieldy for large datasets that feature numerous spatial
locations and time-points, demanding substantial computational and storage
resources. Various methods have been proposed to mitigate these challenges. A
few examples are covariance tapering [Furrer et al., 2006], low-rank
approximation [banerjee2008gaussian, Wikle, 2010], and nearest-neighbor
Gaussian process models [Datta et al., 2016, Finley et al., 2019]. On the
other hand, the limitations of conventional Markov Chain Monte Carlo (MCMC)
algorithms in dealing with vast datasets prompted the development of few
scalable algorithms, including sub-sampling-based approaches [Quiroz et al.,
2019] and stochastic gradient MCMC [Nemeth and Fearnhead, 2021]. Nevertheless,
these methods often introduce additional errors and require extensive
theoretical analysis to ensure their validity. An alternative solution in a
Bayesian framework was proposed by Guhaniyogi and Banerjee [2018] and
Guhaniyogi et al. [2022] which work as our primary motivation in this study.
The key idea here is to employ a divide-and-conquer (hereafter abbreviated as
D&C) strategy for modeling Gaussian processes on massive spatial datasets.

The effectiveness of D&C algorithms, as evident both in practical applications
and theoretical foundations, has spurred our exploration of its potential
within the realm of real estate markets. In this paper, we expand the D&C
approach into the spatio-temporal domain, dubbing it as the “Divide-and-
Conquer method for Spatio-Temporal Big Datasets” (D&C-STBD), and implement it
to comprehensively analyze the house price dynamics of London. The proposed
methodology entails data segmentation, parallel MCMC inference for each
segment, and subsequent merging of the segment posteriors. To facilitate the
last step, we employ the Wasserstein barycenter-based approach detailed in
Shyamalkumar and Srivastava [2022]. It should be noted that our proposed
method allows multiple (possibly unequal number of) and missing (zero)
observations for every space-time combination, thereby obviating the need to
aggregate individual observations within every spatio-temporal unit. Prior
studies in spatio-temporal modeling have predominantly revolved around a
balanced design regarding the number of observations at specific locations and
time-points, cf. Banerjee et al. [2014], Sahu [2022]. However, there has been
limited research that incorporates unequal number of observations per spatial
and temporal unit.

As we show in our application, by embracing individual-level data (i.e.,
property-specific information as opposed to aggregated information per
region), we can sidestep the need for aggregating covariate values into a
generalized spatio-temporal structure. Subsequently, we are able to obtain
more detailed insights about the effects of property characteristics and
regional effects, as well as spatial and temporal dependence in determining
house prices. Our analysis also shows that the model is adaptable and capable
of predicting property prices, even for locations not observed in the training
data. Several interesting aspects, such as the dependence of house price on
carbon emission levels, effects of specific types of amenities, and the trend
patterns are explored in detail. It is imperative to mention that while
several studies have analysed UK house prices [see, e.g., ahlfeldt2014form,
Feng and Jones, 2016, Cook and Watson, 2016, Chi et al., 2021b, Mete and
Yomralioglu, 2022, Blatt et al., 2023], to the best of our knowledge, there
has not been an attempt to effectively address both spatial and temporal
phenomena and dependence in this regard.

Rest of this paper is organized as follows. In Section 2, we introduce the
data, along with pertinent exploratory analysis. Section 3 present our
proposed D&C-STBD method, discussing its implementation and evaluation in
Section 4. The main application is illustrated in Section 5. We conclude with
some important remarks and future scopes of research in Section 6. Technical
details are deferred to the supplementary material.

##  2 Data

###  2.1 Data preparation

In regional science and urban research, open government data is increasingly
being used for spatial analysis of urban areas, with a growing emphasis on
data accessibility [Arribas-Bel, 2014]. In the United Kingdom, HM Land
Registry provides open access to Price Paid Data111Price Paid Dataset URL:
https://www.gov.uk/government/statistical-data-sets/price-paid-data-
downloads(PPD) detailing property sales transactions in England and Wales,
including sale records and basic physical attributes. Another important
resource is the Energy Performance Certificates (EPC) dataset222The Energy
Performance Certificates dataset URL:
https://epc.opendatacommunities.org/docs/guidance, offering insights into
attributes like floor area, energy ratings, carbon emissions, and heating
costs of buildings in England and Wales. Although the PPD dataset contains
valuable information, it lacks essential property attributes such as total
floor area. Without these details, the analysis of house price data is
incomplete. To address this gap, many researchers [e.g., Chi et al., 2021a]
have combined PPD and EPC datasets to analyze property transactions from 2011
to 2019. This combined dataset includes various geospatial attributes like
postcodes, Output Areas (OAs), Lower-layer Super Output Areas (LSOAs), Middle-
layer Super Output Areas (MSOAs), and districts for each location. We use the
same idea and combine the two datasets in the current analysis.

It is critical to point out that the OAs serve as the foundational elements
for spatial census statistics, deliberately structured to ensure a certain
level of socioeconomic coherence. LSOAs, a specialized form of census
geography, typically encompass an average of five OAs each. In turn, MSOAs
combine approximately five LSOAs on average and are situated within district
boundaries. These carefully defined geographic units are meticulously crafted
by the Office for National Statistics (ONS) specifically for analytical
purposes. They are often used as the spatial units, as demanded in specific
research studies, to examine the spatial variability and socioeconomic trends
of key variables. For the house price data, frequent transactions are not
observed at many OA and LSOA levels, resulting in a significant number of
missing values in spatio-temporal settings. Consequently, interpreting data at
these levels may not yield robust insights. Conversely, the MSOA level
exhibits a substantial number of transactions, making it a more viable option
for spatial resolution in our study. Therefore, we opt to utilize the MSOA
level for our spatial analysis, which offers a more comprehensive
understanding of the data landscape. The analysis of house prices at the MSOA
level has been previously investigated in the literature (see Chi et al.
[2021b], Chi et al. [2022]).

We however faced the challenge that latitude and longitude information,
crucial for spatial analysis, are not available in these datasets. This is
resolved by utilizing Free map tools [2023] as it helps us link the spatial
units with latitude and longitude. We utilize the centroid of the geographical
coordinates of all properties within each MSOA as the geographical coordinates
for that MSOA. For the temporal resolution, we opt for a monthly time scale,
which allows for a comprehensive analysis of monthly trends across all MSOAs.
Furthermore, there are (potentially) multiple or no property transactions for
different MSOAs at various time-points. We address both aspects within our
modeling framework, and they are taken care of appropriately in the main
analysis. The methodology will be elaborated in the following section.
Overall, this study specifically centers on the London region, resulting in
983 MSOAs spanning 106 time-points (approximately 8.5 years of monthly data).
After excluding the newly built properties that represent <1% of the
observations and do not have relevant information, we get 104,198 space-time
combinations, with a total of 651,202 property transactions. For ease of
understanding, in Table 1, we summarize the variables considered in this
study, along with their descriptions.

Table 1: Variables considered in the main analysis and their descriptions. Category |  Variable |  Type |  Description
---|---|---|---
Spatial unit |  MSOA |  Coordinates |  Longitude and latitude of a given MSOA.
Temporal unit |  Monthly |  Discrete |  Month in which the sale was completed.
Response |  Price |  Numeric |  House price per square meter.
Covariates |  Area |  Numeric |  Total floor area (in square meters).
|  Rooms |  Categorical |  Number of rooms in the property, categorized in 5 levels: cat1 ($\leqslant 2$ rooms), cat2 (3 rooms), cat3 (4 rooms), cat4 (5 rooms), cat5 ($>5$ rooms).
|  Property type |  Categorical |  Type of property, categorized in 4 levels: Detached, Semi-detached, Terraced, Flats.
|  Ventilation |  Categorical |  Type of mechanical ventilation in property, categorized in 3 levels: Natural, Extract-only, Both supply and extract.
|  Fireplace |  Binary |  1 if there is fireplace in property, 0 otherwise.
|  Wind turbine |  Binary |  1 if there is a wind turbine, 0 otherwise.
|  CO2 emission |  Numeric |  CO2 emissions per year (in $kg/m^{2}$).

###  2.2 Exploratory analysis

Let us start with a visual depiction of the distribution of the key variable
in this study, the price per square meter of all transactions. In the left
panel of Figure 1, we present the histogram and density plot for this
variable. The distribution is clearly right-skewed, which is a common
observation in real estate data [Curto et al., 2015]. This non-Gaussianity
motivates us to adopt a logarithmic transformation, rendering the data normal,
as evident in the right panel of the figure. It must be noted that the
transformation log(price per square meter) is utilized in all further
explorations in this article.

Fig. 1: Density plot on histogram for the distribution of price per square
meter for all properties in the data.

Next, the focus is shifted to exploring the spatial relationship in the house
prices across the MSOAs in the London region. Figure 2 illustrates the
logarithm of the median price per square meter (median taken over the entire
observation period) for every MSOA in the dataset. It clearly reveals a
spatial trend: central areas like Westminster, Chelsea, Camden and the City of
London exhibit higher prices, contrasting with lower prices in the boundaries
of the region, e.g., in areas like Bexley, Croydon, Havering and Hillingdon.
We find that “Westminster 019” has the highest median selling price, while
“Greenwich 001” records the lowest median selling price in the 106 months of
data. Overall, this plot highlights the spatial impact of regions on property
prices in London. Empirically, we also observe that the general trend of the
median price was increasing until September-October 2017, followed by a
decline.

Fig. 2: log median price per square meter across all time intervals for 983
MSOAs in the London region.

While the above hints at spatial correlation in the data, we can further
assess the significance of the spatial dependence using statistical procedure.
To determine this, we employ the Moran’s I statistic [Moran, 1950] for each
month. It is a measure of spatial autocorrelation and can be calculated using

$I=\frac{S}{S_{W}}\frac{\sum_{i=1}^{S}\sum_{j=1}^{S}w_{ij}\left(z_{i}-\bar{z}\right)\left(z_{j}-\bar{z}\right)}{\sum_{i=1}^{S}\left(z_{i}-\bar{z}\right)^{2}},$
(2.1)

where $z_{1},\ldots,z_{S}$ are observations from $S$ number of spatial units,
$\bar{z}$ is the sample mean, and $W$ stands for a suitable spatial weight
matrix. The weight matrix features zeros on the diagonal and other elements
represented by $w_{ij}$, defined as $\exp(-d_{ij})$, where $d_{ij}$ signifies
the great circle distance on an ellipsoid between the locations $s_{i}$ and
$s_{j}$, calculated using the Vincenty method. The term $S_{W}$ represents the
sum of all the elements of $W$. We calculate the Moran’s I statistic for the
price data at all monthly intervals in London, and present them in Figure 3.
There, a black circle denotes statistically significant spatial correlation at
1% level, and we get clear evidence of a strong spatial correlation across all
time-points.

Fig. 3: Moran’s I statistic for the log price for all the monthly time
intervals in the entire London region. Black circle indicates a significant
spatial correlation at 1% level of significance.

Further, to explore the temporal relationships in the price, we compute
autocorrelation function (ACF) of the aforementioned variable for individual
MSOAs in the London area. Our analysis reveals significant ACF values at
multiple lags. As an example, we present the ACF plots for four randomly
selected MSOAs in Figure 4. Other regions also depict similar patterns, and
the plots are omitted for conciseness. The figure demonstrates a significant
temporal correlation in log prices. Considering the pronounced spatial and
temporal correlations identified through Moran’s I statistics and ACF plots,
we proceed to propose a spatio-temporal Gaussian process model for this data.

Fig. 4: ACF plots for the log price for four randomly selected MSOAs. The
dotted lines are corresponding to the critical value at 1% significance level.

Note that all computations in this paper are executed on a shared system with
128 GB RAM and 3 GHz 24-core compute nodes, using RStudio Version
2023.12.0.369 equipped with R version 4.3.2. Moran’s I statistics is
calculated using the ape package [Paradis and Schliep, 2019] and the ACF plots
are made using the forecast package [Hyndman and Khandakar, 2008].

##  3 Divide-and-conquer method for spatio-temporal analysis of house price

Our proposed divide-and-conquer (D&C) methodology is elucidated in this
section. We start by outlining the model suitable to capture the spatio-
temporal dynamics of house price, where potentially varying number of
transactions are observed at every space-time combination. Subsequently, the
entire dataset is split into multiple subsets, and the model is implemented
separately for different subsets, bringing to the fore the D&C framework.
These steps, along with the necessary prior specifications and calculations
for the Bayesian implementation, are explained in the next two subsections.
The last segment of this section provides the details of how to obtain future
predictions from the proposed methodology.

###  3.1 Main model

Throughout this section and after, the term “location” represents specific
geographical points, denoted by latitude and longitude coordinates
$s_{i}\in\mathbb{R}^{2}$ (for $i=1,\ldots,S$). As mentioned earlier, in our
analysis, these locations correspond to the centroids of the MSOAs. We use
$t_{j}\in\mathbb{N}$ to represent the monthly time-points (for $j=1,\ldots,T$)
at which data were recorded. Each combination of location and time-point,
$(s_{i},t_{j})$, is associated with $k_{ij}$ number of observations. The value
of $k_{ij}$ is more than 1 when multiple transactions occur, whereas it
remains zero when no transaction takes place for the specified location and
time point. The latter scenario represents missing data.

We shall use $G=S\times T$ to denote the number of unique space-time
combinations in the data, while the total number of observations is given by
$N=\sum_{i=1}^{S}\sum_{j=1}^{T}k_{ij}$. Also, wherever used below,
$\mathcal{N}_{k}(\cdot,\cdot)$ indicates a $k$-variate Gaussian distribution
with suitable mean and variance parameters, $\bm{I}_{k}$ is the identity
matrix of order $k$, and $\mathbb{I}(A)$ is indicator function for set $A$.

Let $y(s_{i},t_{j},r_{ij})$ denotes the logarithm of the price per square
meter for the $r_{ij}^{th}$ residential property at location $s_{i}$ and time-
point $t_{j}$. Our proposed model, for
$i=1,\ldots,S;\;j=1,\ldots,T;\;r_{ij}=1,\ldots,k_{ij}$, assumes

$y(s_{i},t_{j},r_{ij})=\mu(s_{i},t_{j},r_{ij})+v(s_{i},t_{j})+\epsilon(s_{i},t_{j},r_{ij}),$
(3.1)

where $\mu(s_{i},t_{j},r_{ij})$ denotes the mean structure encompassing an
additive combination of covariate effects and $\epsilon(s_{i},t_{j},r_{ij})$
is an independent and identically distributed white noise process with
variance $\sigma_{\epsilon}^{2}$. For the spatio-temporal dependence, we
assume that all transactions at a given location and time-point share the same
characteristics, and thereby consider the same spatio-temporal process for
such transactions. This is denoted by $v(s_{i},t_{j})$, which is assumed to be
a zero-mean spatio-temporal process.

With these assumptions, the proposed model can be written as

$y(s_{i},t_{j},r_{ij})=\beta_{0}+\sum_{p=1}^{m}\beta_{p}x_{p}(s_{i},t_{j},r_{ij})+v(s_{i},t_{j})+\epsilon(s_{i},t_{j},r_{ij}).$
(3.2)

We can express the above model in matrix form as
$\bm{Y}=\mathbb{X}\bm{\beta}+\bm{B}\bm{V}+\bm{\epsilon}$, where $\bm{Y}$ is an
$N\times 1$ vector representing all observations of the response variable,
arranged first by location and then by time. In a similar fashion,
$\mathbb{X}$ is the required $N\times(m+1)$ design matrix with first column
having all 1 and subsequent columns having the information of the regressors.
The parameter vector
$\bm{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{m})^{\top}$ captures the
covariate effects. The term $\bm{B}$ is a sparse $N\times G$ binary matrix
wherein every row contains 1 in the column that corresponds to the suitable
location-time combination. It is clear that the count of 1s in each column of
$\bm{B}$ matches $k_{ij}$ for that specific location-time combination. The
count of 1s will be zero in case of no observation for particular location and
time point. Finally, $\bm{\epsilon}$ is the $N$-dimensional white noise
vector, while the vector $\bm{V}$, with dimensions $G\times 1$, represents the
spatio-temporal process. It is formed by concatenating columns extracted from
the following matrix:

$\mathbb{V}_{T\times
S}=\begin{bmatrix}v_{s_{1},1}&v_{s_{2},1}&\dots&v_{s_{S},1}\\\
v_{s_{1},2}&v_{s_{2},2}&\dots&v_{s_{S},2}\\\ \vdots&\vdots&\ddots&\vdots\\\
v_{s_{1},T}&v_{s_{2},T}&\dots&v_{s_{S},T}\end{bmatrix}.$ (3.3)

In order to account for the spatial-temporal dependence in observed outcomes,
$v(s_{i},t_{j})$ is modeled as a zero-mean Gaussian process. We assume
$\bm{V}\sim\mathcal{N}_{G}(0,\Sigma)$, where $\Sigma$ is the covariance matrix
with 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\left\\{v(s_{i},t_{u}),v(s_{j},t_{v})\right\\}=\sigma_{v}^{2}\rho_{s}\left(\left\|s_{i}-s_{j}\right\|;\phi_{s}\right)\rho_{t}\left(\left\|t_{u}-t_{v}\right\|;\phi_{t}\right).$
(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 $\nu=0.5$, resulting in the
exponentially decaying pattern. Thus, we have:

$\rho_{s}\left(\left\|s_{i}-s_{j}\right\|;\phi_{s}\right)=\exp\left(-\phi_{s}\left\|s_{i}-s_{j}\right\|\right);\quad\rho_{t}\left(\left\|t_{u}-t_{v}\right\|;\phi_{t}\right)=\exp\left(-\phi_{t}\left\|t_{u}-t_{v}\right\|\right),$
(3.5)

where $\left\|t_{u}-t_{v}\right\|$ represents the absolute time difference,
and $\left\|s_{i}-s_{j}\right\|$ 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
$\Sigma=\sigma_{v}^{2}\left(\Sigma_{s}\otimes\Sigma_{t}\right)$, where
$\otimes$ denotes the Kronecker product, and $\Sigma_{s}$, $\Sigma_{t}$
respectively denote the spatial and temporal correlation matrices, with

$\left(\Sigma_{s}\right)_{ij}=\exp\left(-\phi_{s}\left\|s_{i}-s_{j}\right\|\right);\quad\left(\Sigma_{t}\right)_{uv}=\exp\left(-\phi_{t}\left\|t_{u}-t_{v}\right\|\right).$
(3.6)

###  3.2 Prior Specifications

Let
$\bm{\theta}=(\bm{\beta}^{\top},\sigma_{v}^{2},\sigma_{\epsilon}^{2},\phi_{s},\phi_{t})^{\top}$
be the vector of unknown parameters. We estimate them completely within a
Bayesian framework, primarily using Gibbs sampling [gelfand2000gibbs], for
which suitable prior specifications are essential. For all coefficients in
$\bm{\beta}$, we consider independent Gaussian priors with mean 0 and variance
$c$, where $c$ is relatively large, typically on the order of $10^{4}$. For
the error variances $\sigma_{\epsilon}^{2}$ and $\sigma_{v}^{2}$, we employ
independent inverse gamma priors $IG(a,\lambda)$, where $a$ and $\lambda$ 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 $\lambda=1$.

For estimating the decay parameters $\bm{\phi}=(\phi_{s},\phi_{t})$ in the
covariance functions, many prior studies [e.g., Sahu et al., 2006,
rawat2023spatio] 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 $(\phi_{s},\phi_{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 $\phi_{s}$, selecting
values from 1 to 3 in intervals of 0.2, resulting in the set
$\mathcal{A}_{s}$. For $\phi_{t}$, similarly, we select equidistant points
between 0.2 and 1 in intervals of 0.2, and denote it as $\mathcal{A}_{t}$.
Below, the product $\mathcal{A}_{s}\times\mathcal{A}_{t}$ is indicated as
$\mathcal{A}_{st}$, while $\mathcal{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

$\bm{Y}\mid\bm{\theta},\bm{V}\sim\mathcal{N}_{N}\left(\mathbb{X}\bm{\beta}+\bm{B}\bm{V},\sigma_{\epsilon}^{2}\bm{I}_{N}\right),$
(3.7)

while the prior distributions for all parameters are specified as:

$\begin{split}&\bm{\beta}\sim\mathcal{N}_{m+1}\left(0,c\bm{I}_{m+1}\right),\\\
&\sigma_{\epsilon}^{2}\sim IG(a,\lambda),\;\sigma_{v}^{2}\sim
IG(a,\lambda),\\\ &(\phi_{s},\phi_{t})\sim\mathcal{U}(\mathcal{A}_{st}),\\\
&\bm{V}\sim\mathcal{N}_{G}\left(0,\sigma_{v}^{2}\Sigma_{s}\otimes\Sigma_{t}\right).\end{split}$
(3.8)

Using $f(\cdot)$ to denote the density function in a generic way, the complete
posterior distribution is

$f\left(\bm{\beta},\bm{V},\sigma_{\epsilon}^{2},\sigma_{v}^{2},\phi_{s},\phi_{t}\mid\bm{Y}\right)\propto
f\left(\bm{Y}\mid\bm{\beta},\bm{V},\sigma_{\epsilon}^{2},\sigma_{v}^{2},\phi_{s},\phi_{t}\right)f\left(\bm{V}\mid\sigma_{v}^{2},\phi_{s},\phi_{t}\right)f(\sigma_{v}^{2})f(\phi_{s})f(\phi_{t})f(\sigma_{\epsilon}^{2})f(\bm{\beta}).$
(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(S^{2}T^{2})$ memory units for storage and $O((S\vee
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
$\bm{S}_{q}=(s_{q_{1}},s_{q_{2}},\ldots s_{q_{m_{q}}})$ denote the set of
$m_{q}$ locations in the $q^{th}$ subset, $q\in\\{1,\ldots,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
$\sum_{q=1}^{Q}m_{q}=S$ and $\cup_{q=1}^{Q}\bm{S}_{q}$ being equal to the set
of all locations. In the $q^{th}$ subset, we have the data
$\\{\bm{Y}_{(q)},\mathbb{X}_{(q)}\\}$, with the total number of observations
being $N_{(q)}=\sum_{i=1}^{m_{q}}\sum_{j=1}^{T}k_{ij}$. Clearly,
$\bm{Y}_{(q)}$ is a $N_{(q)}\times 1$ vector, and $\mathbb{X}_{(q)}$ is a
$N_{(q)}\times(m+1)$ matrix representing the covariates corresponding to the
observations for all time-points in the spatial locations within the subset
$\bm{S}_{q}$. Following earlier notations, the spatio-temporal Gaussian
process model for the $q^{th}$ subset is:

$\bm{Y}_{(q)}=\mathbb{X}_{(q)}\bm{\beta}_{(q)}+\bm{B}_{(q)}\bm{V}_{(q)}+\bm{\epsilon}_{(q)}.$
(3.10)

Throughout the discussions in this section, the parameters
$\bm{\theta}_{(q)}=(\bm{\beta}_{(q)},\sigma_{v(q)}^{2},\sigma_{\epsilon(q)}^{2},\phi_{s(q)},\phi_{t(q)})^{\top}$,
the data vector $\bm{Y}_{(q)}$, the spatio-temporal process $\bm{V}_{(q)}$,
the matrices $\mathbb{X}_{(q)}$, $\bm{B}_{(q)}$, and all the available
information $\mathcal{H}_{(q)}$ for subset $q$ are equivalent to their full
dataset counterparts
$\bm{\theta}=(\bm{\beta},\sigma_{v}^{2},\sigma_{\epsilon}^{2},\phi_{s},\phi_{t})^{\top}$,
$\bm{Y}$, $\bm{V}$, $\mathbb{X}$, $\bm{B}$, and $\mathcal{H}$. Also, we use
the notation $p_{q}=N/N_{(q)}$ to simplify some expressions.

####  3.3.2 Second step: modified sampling from subset pseudo posterior
distributions

We note that the $q^{th}$ subset described above contains a fraction
($1/p_{q}$) 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 $q^{th}$ subset. As it contains
$1/p_{q}$ fraction of the full data, the asymptotic variances of the posterior
distributions of $\bm{\theta}_{(q)}$ and $\bm{V}_{(q)}$ are inflated by a
factor of $p_{q}$ compared to that of $\bm{\theta}$ and $\bm{V}$ [Shyamalkumar
and Srivastava, 2022]. Subsequently, we employ stochastic approximation by
raising the likelihood of the parameters
$\bm{\beta}_{(q)},\sigma_{v(q)}^{2},\sigma_{\epsilon(q)}^{2},\phi_{s(q)},\phi_{t(q)},\bm{V}_{(q)}$
in the $q^{th}$ subset to the power of $p_{q}$ [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 $q^{th}$ subset, as there are $N_{(q)}$ data points, the covariance
matrices $\Sigma_{s(q)}$ and $\Sigma_{t(q)}$ have dimensions $m_{q}\times
m_{q}$ and $T\times T$, respectively. Let $G_{(q)}=m_{q}\times T$, and
$\bm{K}_{(q)}$ be a diagonal matrix of order $G_{(q)}\times G_{(q)}$, where
the diagonal elements are $k_{ij}$ for different $(s_{i},t_{j})$ combinations
in the $q^{th}$ subset. The prior distributions for all parameters in subset
$q$ align with those in (3.8), except for $\bm{V}_{(q)}$, which has the
following structure:

$\bm{V}_{(q)}\sim\mathcal{N}_{G_{(q)}}\left(0,\sigma_{v(q)}^{2}\left(\Sigma_{s(q)}\otimes\Sigma_{t(q)}\right)\right).$
(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.

$\begin{split}&\bm{\beta}_{(q)}\mid\mathcal{H}_{(q)}\sim\mathcal{N}_{m+1}\left(\left(\frac{p_{q}\mathbb{X}_{(q)}^{\top}\mathbb{X}_{(q)}}{\sigma_{\epsilon(q)}^{2}}+\frac{\bm{I}_{m+1}}{c}\right)^{-1}\left(\frac{p_{q}\mathbb{X}_{(q)}^{\top}\left(\bm{Y}_{(q)}-\bm{B}_{(q)}\bm{V}_{(q)}\right)}{\sigma_{\epsilon(q)}^{2}}\right)\left(\frac{p_{q}\mathbb{X}_{(q)}^{\top}\mathbb{X}_{(q)}}{\sigma_{\epsilon(q)}^{2}}+\frac{\bm{I}_{m+1}}{c}\right)^{-1}\right),\\\
&\sigma_{\epsilon(q)}^{2}\mid\mathcal{H}_{(q)}\sim
IG\left(a+\frac{N}{2},\frac{p_{q}\left(\left\|\bm{Y}_{(q)}-\mathbb{X}_{(q)}\bm{\beta}_{(q)}-\bm{V}_{(q)}\right\|\right)^{2}}{2}+\lambda\right),\\\
&\sigma_{v(q)}^{2}\mid\mathcal{H}_{(q)}\sim
IG\left(a+\frac{{p_{q}G_{(q)}}}{2},\frac{p_{q}\left(\bm{V}_{(q)}^{\top}\left(\widetilde{\Sigma}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1}\right)\bm{V}_{(q)}\right)}{2}+\lambda\right),\\\
&f(\phi_{s(q)},\phi_{t(q)})\mid\mathcal{H}_{(q)}\propto|\Sigma_{s(q)}|^{-p_{q}T/2}|\Sigma_{t(q)}|^{-p_{q}m_{q}/2}\exp\left(\frac{-p_{q}\bm{V}_{(q)}^{\top}\left(\widetilde{\Sigma}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1})\bm{V}_{(q)}\right)}{\sigma_{v(q)}^{2}}\right)\mathbb{I}\left((\phi_{s(q)},\phi_{t(q)})\in\mathcal{A}_{st}\right),\\\
&\bm{V}_{(q)}\mid\mathcal{H}_{(q)}\sim\mathcal{N}_{G}\left(\left(\frac{\widetilde{\Sigma}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1}}{\sigma_{v(q)}^{2}}+\frac{p_{q}\bm{K}_{(q)}}{\sigma_{\epsilon(q)}^{2}}\right)^{-1}\left(\frac{p_{q}\bm{B}_{(q)}^{\top}(\bm{Y}_{(q)}-\mathbb{X}_{(q)}\bm{\beta}_{(q)})}{\sigma_{\epsilon(q)}^{2}}\right),\left(\frac{\widetilde{\Sigma}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1}}{\sigma_{v(q)}^{2}}+\frac{p_{q}\bm{K}_{(q)}}{\sigma_{\epsilon(q)}^{2}}\right)^{-1}\right).\end{split}$
(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 $\phi_{s(q)}$ and $\phi_{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 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 $q^{th}$ 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 $\Sigma_{s(q)}$ as
${\widetilde{\Sigma}}_{s(q)}^{-1}$, 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
$\bm{V}_{(q)}$ requires inverting a $G_{(q)}\times 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 $\bm{V}_{(q)}$ into
$m_{q}$ subvectors based on the columns in the matrix $\mathbb{V}_{(q)}$.
Then, we take $\bm{Z^{{}^{\prime}}}_{1}$ as the first subvector of
$\mathbb{V}_{(q)}$, and $\bm{Z^{{}^{\prime}}}_{2}$ as a concatenated subvector
of the remaining columns of $\mathbb{V}_{(q)}$, i.e.,
$\bm{V}_{(q)}=\begin{bmatrix}\bm{Z^{{}^{\prime}}}_{1}\
\bm{Z^{{}^{\prime}}}_{2}\end{bmatrix}$. Similarly, the covariance matrix
$\Sigma_{(q)}$ can be written as

$\Sigma_{(q)}=\sigma_{v(q)}^{2}\begin{bmatrix}\Sigma_{11(q)}\otimes\Sigma_{t(q)}\hskip
5.69046pt\Sigma_{12(q)}\otimes\Sigma_{t(q)}\\\
\Sigma_{21(q)}\otimes\Sigma_{t(q)}\hskip
5.69046pt\Sigma_{22(q)}\otimes\Sigma_{t(q)}\end{bmatrix}.$ (3.13)

Then, the conditional distribution of $\bm{Z^{{}^{\prime}}}_{1}$ given
$\bm{Z^{{}^{\prime}}}_{2}=z^{\prime}_{2}$ is
$\mathcal{N}_{T}(\mu_{c(q)},\sigma_{v(q)}^{2}\Sigma_{c(q)})$, where

$\begin{split}\mu_{c(q)}&=\left(\Sigma_{12(q)}\otimes\Sigma_{t(q)}\right)\left(\Sigma_{22(q)}^{-1}\otimes\Sigma_{t(q)}^{-1}\right)z^{\prime}_{2};\\\
\Sigma_{c(q)}&=\left(\Sigma_{11(q)}\otimes\Sigma_{t(q)}\right)-\left(\Sigma_{12(q)}\otimes\Sigma_{t(q)}\right)\left(\Sigma_{22(q)}^{-1}\otimes\Sigma_{t(q)}^{-1}\right)\left(\Sigma_{21(q)}\otimes\Sigma_{t(q)}\right).\end{split}$
(3.14)

Using the properties of the Kronecker product, the expression can be further
simplified as

$\mu_{c(q)}=\left(\left(\Sigma_{12(q)}\Sigma_{22(q)}^{-1}\right)\otimes\bm{I}_{T}\right)z^{\prime}_{2};\quad\Sigma_{c(q)}=\left(\Sigma_{11(q)}-\Sigma_{12(q)}\Sigma_{22(q)}^{-1}\Sigma_{21(q)}\right)\otimes\Sigma_{t(q)}.$
(3.15)

Now, following the above, we can write

$\begin{split}\bm{V}_{1(q)}\mid\mathcal{H}_{(q)}&\sim\mathcal{N}_{T}\left(\left(\frac{\Sigma_{c(q)}^{-1}}{\sigma_{v(q)}^{2}}+\frac{p_{q}\bm{K}_{1(q)}}{\sigma_{\epsilon(q)}^{2}}\right)^{-1}\left(\frac{\Sigma_{c(q)}^{-1}\mu_{c(q)}}{\sigma_{v(q)}^{2}}+\frac{p_{q}\left(\bm{B}_{1(q)}^{T}\left(\bm{Y}_{1(q)}-\mathbb{X}_{1(q)}\bm{\beta}_{(q)}\right)\right)}{\sigma_{\epsilon(q)}^{2}}\right)\right.,\\\
&\hskip
72.26999pt\left.\left(\frac{\Sigma_{c(q)}^{-1}}{\sigma_{v(q)}^{2}}+\frac{p_{q}\bm{K}_{1(q)}}{\sigma_{\epsilon(q)}^{2}}\right)^{-1}\right),\end{split}$
(3.16)

where $\bm{Y}_{1(q)}$ be the data for the first location and all time-points
from the $q^{th}$ subset. The matrix $\mathbb{X}_{1(q)}$ is defined
accordingly. $\bm{K}_{1(q)}$ is a $T\times T$ diagonal matrix with each
diagonal element being equal to $k_{1j}$ for each $(s_{q_{1}},t_{j})$.
$\bm{B}_{1(q)}$ is a submatrix of $\bm{B}_{(q)}$, consisting of the first $T$
columns and $\sum_{j=1}^{T}k_{1j}$ rows, resulting in a sparse matrix of size
$\sum_{j=1}^{T}k_{1j}\times T$. Other notations have the same definitions as
before. We can compute the conditional posterior distributions of
$\bm{V}_{2(q)},\ldots,\bm{V}_{m_{q}(q)}$ in an identical fashion. Finally, to
compute the posterior distribution of the full vector $\bm{V}_{(q)}$, we
iteratively employ the conditional distributions in a Gibbs sampler. This
approach generates independent posterior samples after an initial burn-in
period. Consequently, the posterior distribution of $\bm{V}_{(q)}$ is computed
efficiently without incurring the computational cost of inverting the
$G_{(q)}\times 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, $\mathcal{P}(M)$ be the space of all probability measures on $M$
and $\mathcal{P}_{2}(M)$ be the set of all probability measures on $M$ with
finite second moments. For $\varrho,\xi\in\mathcal{P}_{2}(M)$,
$\Pi(\varrho,\xi)$ denotes the set of all probability measures on $M\times M$
with marginals $\varrho$ and $\xi$. Then, with $d(\cdot,\cdot)$ indicating the
euclidean metric, the Wasserstein distance of order 2 between $\varrho$ and
$\xi$ is defined as

$W_{2}(\varrho,\xi)=\left\\{\inf_{\pi\in\Pi(\varrho,\xi)}\int_{M\times
M}d(x_{1},x_{2})^{2}\Pi(x_{1},x_{2})\right\\}^{1/2}.$ (3.17)

As Bickel and Freedman [1981] discussed, the convergence in $W_{2}$ on
$\mathcal{P}_{2}(M)$ implies weak convergence of probability measures and
convergence of the second moment. Keeping that in view, if
$\bm{\theta}_{(1)},\ldots,\bm{\theta}_{(Q)}$ are the $Q$ subset posterior
distributions, then the approximated posterior $\hat{\bm{\theta}}$ is defined
as

$\hat{\bm{\theta}}=\operatorname*{arg\,min}_{\bm{\theta}\in\mathcal{P}_{2}(M)}\frac{1}{Q}\sum_{i=1}^{Q}W_{2}^{2}(\bm{\theta},\bm{\theta}_{(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
$\hat{\bm{\theta}}$, 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.

Input : Spatio temporal dataset (Y,$\mathbb{X}$), where Y is the response
vector and $\mathbb{X}$ is the complete information on regressors. Number of
subsets $Q$, the total number of post burn-in iterations on every subset is
$L$, size of each subset, prior distribution for all the parameters.

Output :

  1. 1.

$\hat{\bm{\theta}}^{(1)},\ldots,\hat{\bm{\theta}}^{(QL)}$ as approximate WASP
draws.

  2. 2.

The approximated $\gamma$ quantile for $\bm{\theta}$ is then calculated as:

$\hat{\bm{\theta}}^{\gamma}=\frac{1}{QL}\sum_{i=1}^{QL}\hat{\bm{\theta}}^{(i)\gamma}.$

Algorithm:

  1. 1.

Divide data into $Q$ subsets, with $q^{th}$ subset as
$(\bm{Y}_{(q)},\mathbb{X}_{(q)})$, each containing data from different
locations but all time-points.

  2. 2.

Compute the subset posteriors for all parameters of all subsets as discussed
in Section 3.3.2.

  3. 3.

for $q=1,\ldots Q$ do

Obtain subset posterior draws (every tenth MCMC sample after convergence) for
all parameters of all subsets at $l^{th}$ iteration using Gibbs and slice
sampling: $\bm{\theta}_{(q)}^{(l)}(l=1,\ldots,L)$.

end for

* 4.

Compute mean vectors and covariance matrices of subset posterior distributions
and the approximate WASP posterior distribution:
$\hat{\mu}_{(q)},\hat{\bar{\mu}},\hat{\Sigma}_{(q)},\hat{\bar{\Sigma}}$. 5.
Center and scale subset posterior draws for $q=1,\ldots,Q;l=1,\ldots,L$:
$\hat{c}_{(q)}^{(l)}=\hat{\Sigma}^{-1/2}_{(q)}\left(\bm{\theta}_{(q)}^{(l)}-\hat{\mu}_{(q)}\right)$
6. Define combined posterior using $q^{th}$ subset posterior draws of
$\bm{\theta}$:
$\hat{\bm{\theta}}^{(l^{{}^{\prime}})}=\hat{\bar{\mu}}+{\hat{\bar{\Sigma}}}^{1/2}\hat{c}_{(q)}^{(l)},\hskip
2.84544ptl^{{}^{\prime}}=(q-1)L+l,$ where
$\hat{\bm{\theta}}^{(l^{{}^{\prime}})}$ is the approximate WASP draw for
$\bm{\theta}$.

Algorithm 1 Proposed divide-and-conquer methodology for large spatio temporal
data (D&C-STBD)

A brief outline of the WASP procedure is also of the essence here. First,
posterior draws of $\bm{\theta}_{(1)},\ldots,\bm{\theta}_{(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 $\hat{\mu}_{(q)}$ and covariance matrices $\hat{\Sigma}_{(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
$\hat{\bar{\mu}}$ and covariance matrix $\hat{\bar{\Sigma}}$. 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 $\bm{V}$ from the subset posterior samples of $v_{(q)}$, we
calculate the mean of the samples within each subset to get an estimate for
that subset. This leads to $Q$ number of posterior vectors, each corresponding
to a different location and time-point. Finally, these vectors are stacked to
form the final estimate of $\bm{V}$.

###  3.4 Prediction

Our defined spatio-temporal Gaussian process model offers the ability to
predict the outcome variable for any spatial location and time-point, which
are not observed in the data. This is the primary advantage of this method.
Specifically, for a residential property $r^{\prime}$ in a new location
$s^{\prime}$ at a future time-point $t^{\prime}$, the predicted outcome
variable $y(s^{\prime},t^{\prime},r^{\prime})$ is conditionally independent of
the observed outcomes $\bm{Y}$ given the spatio-temporal Gaussian process
model $v(s^{\prime},t^{\prime})$ and the estimated $\bm{\beta}$. If
$\bm{X}(s^{\prime},t^{\prime},r^{\prime})$ is the covariate vector for the new
observation, then the predicted outcome variable follows a Gaussian
distribution:

$y(s^{\prime},t^{\prime},r^{\prime})\mid
v(s^{\prime},t^{\prime})\sim\mathcal{N}\left(\bm{X}^{\top}(s^{\prime},t^{\prime},r^{\prime})\hat{\bm{\beta}}+v(s^{\prime},t^{\prime}),\sigma_{\epsilon}^{2}\right),$
(3.19)

We estimate the model parameters using a Gibbs sampler on the training data,
yielding posterior estimates for $\bm{\beta}$ and $\sigma_{\epsilon}^{2}$. To
find the conditional distribution of $v(s^{\prime},t^{\prime})$ given
$\bm{V}$, we use the posterior value of $\bm{V}$ from training data estimates
and apply the multivariate Gaussian theory related to the conditional
distribution. We treat $v(s^{\prime},t^{\prime})$ as the first subvector of
the multivariate joint distribution, and $\bm{V}$ obtained from the training
dataset as the second subvector. The joint distribution of
$v(s^{\prime},t^{\prime})$ and $\bm{V}$ is given by

$\begin{pmatrix}v(s^{\prime},t^{\prime})\\\
\bm{V}\end{pmatrix}\sim\mathcal{N}_{G+1}\begin{pmatrix}\begin{pmatrix}0\\\
0\end{pmatrix},\sigma_{v}^{2}\begin{bmatrix}1&\Sigma^{\prime\top}\\\
\Sigma^{\prime}&\Sigma_{s}\otimes\Sigma_{t}\end{bmatrix}\end{pmatrix},\;\text{with}\;\Sigma^{\prime}=\Sigma_{s}(s-s^{\prime})\otimes\Sigma_{t}(t-t^{\prime}),$
(3.20)

where $\Sigma_{s(s-s^{\prime})}$ is a vector of order $S$ with $i^{th}$ entry
as $\rho_{s}\left(||s_{i}-s^{\prime}||;\phi_{s}\right)$ and
$\Sigma_{t(t-t^{\prime})}$is a vector of order $T$ with $j^{th}$ entry as
$\rho_{t}\left(||t_{j}-t^{\prime}||;\phi_{t}\right)$. The conditional
distribution is then given by

$v(s^{\prime},t^{\prime})\mid\bm{V}\sim\mathcal{N}\left(\Sigma^{\prime}(\Sigma_{s}^{-1}\otimes\Sigma_{t}^{-1})\bm{V},\sigma_{v}^{2}(1-\Sigma^{\prime}(\Sigma_{s}^{-1}\otimes\Sigma_{t}^{-1})\Sigma^{\prime\top})\right).$
(3.21)

Thus, by predicting $v(s^{\prime},t^{\prime})$ using the above conditional
distribution and the estimates of $\bm{\beta}$ from the Gibbs sampling
procedure, we can predict $y(s^{\prime},t^{\prime},r^{\prime})$ 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^{\prime},t^{\prime})$ for
a location $s^{\prime}$ and a time-point $t^{\prime}$ 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^{\prime},t^{\prime})$ using
a similar process as for the whole dataset, resulting in $Q$ realizations of
$v(s^{\prime},t^{\prime})$. These realizations are denoted as
${v}_{q}(s^{\prime},t^{\prime})$, for $1\leqslant q\leqslant Q$. The joint
distribution of ${v}_{q}(s^{\prime},t^{\prime})$ and $\bm{V_{(q)}}$ is given
by

$\begin{pmatrix}{v}_{q}(s^{\prime},t^{\prime})\\\
\bm{V_{(q)}}\end{pmatrix}\sim\mathcal{N}_{G_{(q)}+1}\begin{pmatrix}\begin{pmatrix}0\\\
0\end{pmatrix},{\sigma_{v}}_{(q)}^{2}\begin{bmatrix}1&\Sigma_{(q)}^{\prime\top}\\\
\Sigma_{(q)}^{\prime}&\Sigma_{s(q)}\otimes\Sigma_{t(q)}\end{bmatrix}\end{pmatrix},\;\text{with}\;\Sigma_{(q)}^{\prime}=\Sigma_{s(q)}(s-s^{\prime})\otimes\Sigma_{t(q)}(t-t^{\prime}),$
(3.22)

where $\Sigma_{{s(q)}(s-s^{\prime})}$ is a vector of order $S_{(q)}$ with
$i^{th}$ entry as $\rho_{s}\left(||s_{i}-s^{\prime}||;\phi_{s(q)}\right)$ and
$\Sigma_{t(t-t^{\prime})}$is a vector of order $T$ with $j^{th}$ entry as
$\rho_{t}\left(||t_{j}-t^{\prime}||;\phi_{t(q)}\right)$. The conditional
distribution is then given by

${v}_{q}(s^{\prime},t^{\prime})\mid\bm{V_{(q)}}\sim\mathcal{N}\left(\Sigma_{(q)}^{\prime}({\widetilde{\Sigma}}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1})\bm{V_{(q)}},\sigma_{v(q)}^{2}(1-\Sigma_{(q)}^{\prime}({\widetilde{\Sigma}}_{s(q)}^{-1}\otimes\Sigma_{t(q)}^{-1})\Sigma_{(q)}^{\prime\top})\right).$
(3.23)

With the above predicted values of ${v_{(q)}}(s^{\prime},t^{\prime})$ 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
$q^{th}$ subset data.

$y_{(q)}(s^{\prime},t^{\prime},r^{\prime})\mid
v_{(q)}(s^{\prime},t^{\prime})\sim\mathcal{N}\left(\bm{X}_{(q)}^{\top}(s^{\prime},t^{\prime},r^{\prime})\hat{\bm{\beta}}_{(q)}+v_{(q)}(s^{\prime},t^{\prime}),{\sigma_{\epsilon}}_{(q)}^{2}\right).$
(3.24)

Finally, to predict $v(s^{\prime},t^{\prime})$ and the response variable
$y(s^{\prime},t^{\prime},r^{\prime})$ using the entire training dataset, we
take the median of ${v}_{q}(s^{\prime},t^{\prime})$ and
$y_{(q)}(s^{\prime},t^{\prime},r^{\prime})$ respectively, for $1\leqslant
q\leqslant 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.

##  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 sommervoll2019learning. 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\times 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 $\mathcal{D}$, with an observation for location $s_{i}$ at time-
point $t_{k}$ be denoted as $y(s_{i},t_{k},r_{ik})$. The corresponding fitted
value is $\hat{y}(s_{i},t_{k},r_{ik})$. To understand the goodness of fit of
every model, we start with the ubiquitous coefficient of determination
($R^{2}$) 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, the
variance of the absolute error (VAE) can help to assess the robustness or
stability of the model’s performance, especially in scenarios where consistent
predictions are essential.

In general, a model with lower MAE, MAPE, RMSE, and VAE values is regarded as
more accurate. Mathematically, these are defined as:

$\begin{split}&\mathrm{MAE}=\frac{1}{\left\lvert\mathcal{D}\right\rvert}\sum_{\mathcal{D}}\left\lvert\hat{y}(s_{i},t_{k},r_{ik})-y(s_{i},t_{k},r_{ik})\right\rvert,\\\
&\mathrm{MAPE}=\frac{100\%}{\left\lvert\mathcal{D}\right\rvert}\sum_{\mathcal{D}}\left\lvert\frac{\hat{y}(s_{i},t_{k},r_{ik})-y(s_{i},t_{k},r_{ik})}{y(s_{i},t_{k},r_{ik})}\right\rvert,\\\
&\mathrm{RMSE}=\sqrt{\frac{1}{\left\lvert\mathcal{D}\right\rvert}\sum_{\mathcal{D}}(\hat{y}(s_{i},t_{k},r_{ik})-y(s_{i},t_{k},r_{ik}))^{2}},\\\
&\mathrm{VAE}=\frac{1}{\left\lvert\mathcal{D}\right\rvert}\sum_{\mathcal{D}}\left(\left\lvert\hat{y}(s_{i},t_{k},r_{ik})-y(s_{i},t_{k},r_{ik})\right\rvert-\mathrm{MAE}\right)^{2}.\end{split}$
(4.1)

Finally, we steer to the prediction accuracy of our proposed model. In this
case, the whole data $\mathcal{D}$ is split into a training set (call it
$\mathcal{D}_{\mathrm{tr}}$) and a test set (call it
$\mathcal{D}_{\mathrm{te}}$). We fit the proposed model to
$\mathcal{D}_{\mathrm{tr}}$ and use the estimated parameters to predict the
log house price per square meter for each entry in
$\mathcal{D}_{\mathrm{te}}$. In this case, using the true and predicted values
for $\mathcal{D}_{\mathrm{te}}$, we primarily calculate the mean absolute
percentage error to understand the forecasting performance of the model. As we
illustrate in Section 5.4, by considering different types of test sets, we can
analyze and comment on the robustness of the predictive capability of our
methodology.

##  5 Results and discussions

###  5.1 Estimated model

We begin by presenting the estimates of the parameters obtained from fitting
our proposed model to the entire data. The results are shown in Table 2, which
displays the posterior means of the parameters along with their 95% credible
intervals. It helps us carefully examine the key factors influencing property
prices in the London region, and provide valuable insights into the complex
dynamics underlying property pricing. To begin, we note that the credible
intervals for all the covariates used in our study, except for the presence of
wind turbine, are away from zero, thereby suggesting significant impact of
these factors.

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)
$\sigma_{\epsilon}^{2}$ | 0.043 | (0.043 , 0.043 )
$\sigma_{v}^{2}$ | 0.083 | (0.070 , 0.103)
$\phi_{s}$ | 2.402 | (1.486 , 3.041)
$\phi_{t}$ | 0.528 | (0.183 , 0.824)

First, we notice that both the linear and quadratic terms of the trend
coefficients are significant, with the linear coefficient being positive and
the quadratic coefficient being negative. Moreover, the magnitude of the
linear trend coefficient is quite higher than the other one, indicating that
the house price is increasing over time. To demonstrate the impact of time
trend, we utilize the intercept and time trend coefficients to generate fitted
mean log price values for all MSOAs across the 106 time-points. As shown in
Figure 5, it depicts a curvilinear relationship between the time trend and the
response variable. The average log price per square meter increases at a
diminishing rate over time and eventually flattens out. This trend is akin to
the concept of diminishing marginal returns in economics, where the later time
periods produce smaller increases in the variable of interest. These results
are consistent with a study by Office for national statistics [2018], where
the recent trends in house price growth in London was analyzed using the UK
House Price Index from January 2012 to July 2018. The study found that the
growth in house prices had declined since 2016.

Fig. 5: Average trend of the house price (log-transformed) for all locations
over 106 time.

Turning attention to the estimates of the covariates, detailed interpretations
of their impact on the response variable are discussed next.

####  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 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 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 | $R^{2}$ | 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 ($R^{2}$
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 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 $\sigma_{v}^{2}$ is almost double in magnitude compared to
$\hat{\sigma}_{\epsilon}^{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
$\hat{\phi}_{s}=2.402$ and $\hat{\phi}_{t}=0.528$ reveal practical
implications. With the equation $\exp(-\phi d)\approx 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.

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.

Fig. 7: ACF plots for the residuals for four randomly selected MSOAs. The
dotted lines are corresponding to the critical value at 1% significance level.

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
have a symmetric distribution around 0, with a hint of a heavier tail than a
normal distribution. It might be suggestive of a non-Gaussian residual
distribution, and exploring this further could be a valuable direction for
future research. We discuss this in Section 6.

###  5.4 Predictive performance

As a final exploration, we assess our model’s predictive capabilities in two
scenarios, looking at a future horizon of one year in both cases. In the first
scenario, the training set comprises 94 monthly observations for 983 MSOAs,
totaling 597,172 sample observations, while the test set includes the last 12
monthly observations (November 2018 to October 2019) for these MSOAs, totaling
54,030 samples. In the second scenario, our objective is to assess the
predictive accuracy for unobserved MSOAs. To that end, we consider 23 randomly
selected MSOAs. The training set includes the initial 94 time points for the
remaining 960 MSOAs, totaling 584,672 samples. The test set in this case
comprises 1,068 samples from the 23 selected MSOAs over the last 12 time
points.

For both scenarios, we employ a similar division strategy to the one used when
fitting the entire dataset. In the first case, the 983 MSOAs are divided into
20 subsets, with 19 subsets containing 49 MSOAs each, and the $20^{th}$ subset
keeping 52 locations. In the second scenario, we split the 960 MSOAs into 20
subsets, with each subset containing 48 MSOAs. We fit the proposed model using
the method described in Section 3.3.3, obtain the required parameter
estimates, and make predictions for each location and time in the test set
using the procedure outlined in Section 3.4.

First, to understand the predictive ability in the spatial aspect, we
calculate the Mean Absolute Percentage Error (MAPE) for all MSOAs in the first
scenario, where the average is taken over the 12 time-points in the test set.
These values are displayed in Figure 8. There, it can be observed that the
MAPE varies across MSOAs, ranging from 1.3% to 11.25%. Notably, only four
MSOAs, namely “Westminster 018”, “Kensington and Chelsea 012”, “Westminster
019” and , “Kensington and Chelsea 008” exhibit MAPE values exceeding 10%.
This indicates potential areas where our predictive models may require further
refinement.

Fig. 8: Prediction MAPE (in %) in log price per square meter for all MSOAs for
12 months in the first case.

To gain a broader perspective, we delve into local authority (region) wise
MAPE, aiming to discern trends in error rates across different regions. Our
investigation reveals that region-wise MAPE spans from 2.3% to 8.3%.
Particularly noteworthy is the highest MAPE observed in the “Kensington and
Chelsea” region, contrasted with the lowest MAPE in the “Kingston upon Thames”
region. Conversely, for the remaining regions, MAPE values remain below 7%,
suggesting relatively lower prediction errors. These discrepancy could
potentially be attributed to market instability in the UK real estate sector
during the corresponding period in 2018 and 2019. In an attempt to gain deeper
insights, a detailed examination of property-wise errors is conducted next. In
the first scenario, the test data contains a total of 54,030 house sales. We
calculate the Absolute Percentage Error (APE) in forecasting the price for
each transaction and categorize them into different APE brackets, as detailed
under Case 1 in Table 4. It shows that about 75% of properties have APE of
less than 5%, while only about 3.5% of properties exceed an APE of 10%. These
results underscore the effectiveness of our proposed approach in providing
accurate predictions at the individual property level. The findings are
identical in the second case as well.

Table 4: Distribution of houses corresponding to different levels of absolute percentage errors for both cases. | Case 1 (total sample size: 54030) | Case 2 (total sample size: 1068)
---|---|---
APE | Number of cases | Percentage of cases | Number of cases | Percentage of cases
Less than 3% | 25604 | 47.39% | 508 | 47.57%
3 to 5% | 14267 | 26.41% | 306 | 28.65%
5 to 10% | 12238 | 22.65% | 236 | 22.10%
More than 10% | 1921 | 3.55% | 18 | 1.69%

Next, focus on the second scenario with the 23 locations unobserved in the
training data. In Figure 9, we present the MAPE for these locations according
to different time-points in the future. It helps us quantify the predictive
performance over time for unobserved locations. In the figure, the right panel
corresponds to this whereas in the left panel, we include the accuracy for the
scenario where these locations are present in the training data.
Interestingly, in both scenarios, the MAPE ranges between 2.9% and 4.6%, which
confirms that the proposed method can offer good predictive certainty for
properties from new locations as well. Moreover, the results showcase a
consistently accurate level 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%.

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 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 [2014] Arribas-Bel, D., 2014.  Accidental, open and everywhere: Emerging data sources for the understanding of cities.  Applied Geography 49, 45–53.
  * Banerjee et al. [2014] Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014.  Hierarchical Modeling and Analysis for Spatial Data.  CRC Press.
  * Beamonte et al. [2010] 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 and Freedman [1981] Bickel, P.J., Freedman, D.A., 1981\.  Some asymptotic theory for the bootstrap.  The annals of statistics 9, 1196–1217.
  * Blanco and Neri [2023] Blanco, H., Neri, L., 2023.  Knocking It Down and Mixing It Up: The Impact of Public Housing Regenerations.  Technical Report. IZA Discussion Papers.
  * Blatt et al. [2023] Blatt, D., Chaudhuri, K., Manner, H., 2023.  A changepoint analysis of uk house price spillovers.  Regional Studies 57, 1223–1238.
  * Botchkarev [2019] Botchkarev, A., 2019.  A new typology design of performance metrics to measure errors in machine learning regression algorithms.  Interdisciplinary Journal of Information, Knowledge, and Management 14, 045–076.
  * Can [1990] Can, A., 1990.  The measurement of neighborhood dynamics in urban house prices.  Economic geography 66, 254–272.
  * Chegut et al. [2016] Chegut, A., Eichholtz, P., Holtermans, R., 2016.  Energy efficiency and economic value in affordable housing.  Energy Policy 97, 39–49.
  * Chi et al. [2021a] Chi, B., Dennett, A., Oléron-Evans, T., Morphet, R., 2021a.  A new attribute-linked residential property price dataset for england and wales, 2011 to 2019.  UCL Open: Environment Preprint .
  * Chi et al. [2021b] Chi, B., Dennett, A., Oléron-Evans, T., Morphet, R., 2021b.  Shedding new light on residential property price variation in england: A multi-scale exploration.  Environment and Planning B: Urban Analytics and City Science 48, 1895–1911.
  * Chi et al. [2022] Chi, B., Dennett, A., Oléron-Evans, T., Morphet, R., 2022\.  Delineating the spatio-temporal pattern of house price variation by local authority in england: 2009 to 2016.  Geographical Analysis 54, 219–238.
  * Cook and Watson [2016] Cook, S., Watson, D., 2016.  A new perspective on the ripple effect in the uk housing market: Comovement, cyclical subsamples and alternative indices.  Urban Studies 53, 3048–3062.
  * Curto et al. [2015] Curto, R., Fregonara, E., Semeraro, P., 2015.  Listing behaviour in the italian real estate market.  International Journal of Housing Markets and Analysis .
  * Datta et al. [2016] Datta, A., Banerjee, S., Finley, A.O., Gelfand, A.E., 2016\.  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association 111, 800–812.
  * Deb and Tsay [2019] Deb, S., Tsay, R.S., 2019.  Spatio-temporal models with space-time interaction and their applications to air pollution data.  Statistica Sinica 29, 1181–1207.
  * Dubin and Sung [1990] Dubin, R.A., Sung, C.H., 1990\.  Specification of hedonic regressions: non-nested tests on measures of neighborhood quality.  Journal of Urban Economics 27, 97–110.
  * Feng and Jones [2016] Feng, Y., Jones, K., 2016.  Postcode or census geography? an examination of neighbourhood classification for house price predictions, in: The 22nd Annual Pacific Rim Real Estate Society Conference, pp. 1–12.
  * Finley et al. [2019] Finley, A.O., Datta, A., Cook, B.D., Morton, D.C., Andersen, H.E., Banerjee, S., 2019\.  Efficient algorithms for bayesian nearest neighbor gaussian processes.  Journal of Computational and Graphical Statistics 28, 401–414.
  * Fotheringham et al. [2003] Fotheringham, A.S., Brunsdon, C., Charlton, M., 2003.  Geographically weighted regression: the analysis of spatially varying relationships.  John Wiley & Sons.
  * Free map tools [2023] Free map tools, 2023.  UK Postcodes with Latitude and Longitude.  URL: https://www.freemaptools.com/download-uk-postcode-lat-lng.htm.
  * Fuerst et al. [2015] Fuerst, F., McAllister, P., Nanda, A., Wyatt, P., 2015\.  Does energy efficiency matter to home-buyers? an investigation of epc ratings and transaction prices in england.  Energy Economics 48, 145–156.
  * Furrer et al. [2006] Furrer, R., Genton, M.G., Nychka, D., 2006.  Covariance tapering for interpolation of large spatial datasets.  Journal of Computational and Graphical Statistics 15, 502–523.
  * Gelfand et al. [2004] Gelfand, A.E., Ecker, M.D., Knight, J.R., Sirmans, C., 2004\.  The dynamics of location in home price.  The journal of real estate finance and economics 29, 149–166.
  * Gerassimenko et al. [2023] Gerassimenko, A., Defau, L., De Moor, L., 2023.  The impact of energy certificates on sales and rental prices: a comparative analysis.  International Journal of Housing Markets and Analysis .
  * Geweke et al. [1991] Geweke, J.F., et al., 1991.  Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments.  Technical Report. Federal Reserve Bank of Minneapolis.
  * Guhaniyogi and Banerjee [2018] Guhaniyogi, R., Banerjee, S., 2018\.  Meta-kriging: Scalable bayesian modeling and inference for massive spatial datasets.  Technometrics 60, 430–444.
  * Guhaniyogi et al. [2022] Guhaniyogi, R., Li, C., Savitsky, T., Srivastava, S., 2022\.  Distributed bayesian inference in massive spatial data.  Statist. Sci .
  * Hijmans [2021] Hijmans, R.J., 2021.  Introduction to the “geosphere” package (version 1.5-14) .
  * Holly et al. [2010] Holly, S., Pesaran, M.H., Yamagata, T., 2010.  A spatio-temporal model of house prices in the USA.  Journal of Econometrics 158, 160–173.
  * Hyndman and Khandakar [2008] Hyndman, R.J., Khandakar, Y., 2008\.  Automatic time series forecasting: the forecast package for r.  Journal of statistical software 27, 1–22.
  * Li et al. [2017] Li, C., Srivastava, S., Dunson, D.B., 2017.  Simple, scalable and accurate posterior interval estimation.  Biometrika 104, 665–680.
  * Liu [2013] Liu, X., 2013.  Spatial and temporal dependence in house price prediction.  The Journal of Real Estate Finance and Economics 47, 341–369.
  * Ma and Li [2017] Ma, H., Li, J., 2017.  The impacts of supply and demand analysis on the price of the real estate market, in: 7th International Conference on Education, Management, Information and Mechanical Engineering (EMIM 2017), Atlantis Press. pp. 1881–1885.
  * Mete and Yomralioglu [2022] Mete, M.O., Yomralioglu, T., 2022\.  A hybrid approach for mass valuation of residential properties through geographic information systems and machine learning integration.  Geographical Analysis .
  * Minsker et al. [2014] Minsker, S., Srivastava, S., Lin, L., Dunson, D., 2014\.  Scalable and robust bayesian inference via the median posterior, in: International conference on machine learning, PMLR. pp. 1656–1664.
  * Moran [1950] Moran, P.A., 1950.  Notes on continuous stochastic phenomena.  Biometrika 37, 17–23.
  * Neal [2003] Neal, R.M., 2003.  Slice sampling.  The annals of statistics 31, 705–767.
  * Nemeth and Fearnhead [2021] Nemeth, C., Fearnhead, P., 2021\.  Stochastic gradient markov chain monte carlo.  Journal of the American Statistical Association 116, 433–450.
  * Office for national statistics [2018] Office for national statistics, 2018.  Exploring recent trends in the London housing market.  URL: https://www.ons.gov.uk/economy/inflationandpriceindices/articles/exploringrecenttrendsinthelondonhousingmarket/2018-09-19.
  * Ou et al. [2021] Ou, R., Sen, D., Dunson, D., 2021.  Scalable bayesian inference for time series via divide-and-conquer.  arXiv preprint arXiv:2106.11043 .
  * Pace et al. [1998] Pace, R.K., Barry, R., Clapp, J.M., Rodriquez, M., 1998\.  Spatiotemporal autoregressive models of neighborhood effects.  The Journal of Real Estate Finance and Economics 17, 15–33.
  * Pace and Gilley [1998] Pace, R.K., Gilley, O.W., 1998\.  Generalizing the ols and grid estimators.  Real Estate Economics 26, 331–347.
  * Paradis and Schliep [2019] Paradis, E., Schliep, K., 2019\.  ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R.  Bioinformatics 35, 526–528.  doi:10.1093/bioinformatics/bty633.
  * Quiroz et al. [2019] Quiroz, M., Kohn, R., Villani, M., Tran, M.N., 2019\.  Speeding up mcmc by efficient data subsampling.  Journal of the American Statistical Association 114, 831–843.
  * Rosen [1974] Rosen, S., 1974.  Hedonic prices and implicit markets: product differentiation in pure competition.  Journal of political economy 82, 34–55.
  * Sahu [2022] Sahu, S., 2022.  Bayesian modeling of spatio-temporal data with R.  CRC Press.
  * Sahu et al. [2006] Sahu, S.K., Gelfand, A.E., Holland, D.M., 2006.  Spatio-temporal modeling of fine particulate matter.  Journal of Agricultural, Biological, and Environmental Statistics 11, 61–86.
  * Sahu et al. [2010] Sahu, S.K., Gelfand, A.E., Holland, D.M., 2010.  Fusing point and areal level space–time data with application to wet deposition.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 59, 77–103.
  * Shyamalkumar and Srivastava [2022] Shyamalkumar, N.D., Srivastava, S., 2022\.  An algorithm for distributed bayesian inference.  Stat 11, e432.
  * Soltani et al. [2021] Soltani, A., Pettit, C.J., Heydari, M., Aghaei, F., 2021\.  Housing price variations using spatio-temporal data mining techniques.  Journal of Housing and the Built Environment , 1–29.
  * Srivastava et al. [2018] Srivastava, S., Li, C., Dunson, D.B., 2018.  Scalable bayes via barycenter in wasserstein space.  The Journal of Machine Learning Research 19, 312–346.
  * Szabó et al. [2019] Szabó, B., Van Zanten, H., et al., 2019\.  An asymptotic analysis of distributed nonparametric methods.  Journal of Machine Learning Research 20, 1–30.
  * Teye and Ahelegbey [2017] Teye, A.L., Ahelegbey, D.F., 2017\.  Detecting spatial and temporal house price diffusion in the Netherlands: A Bayesian network approach.  Regional Science and Urban Economics 65, 56–64.
  * Wikle [2010] Wikle, C.K., 2010.  Low-rank representations for spatial processes.  Handbook of spatial statistics 107, 118.
  * Zhang [2006] Zhang, F., 2006.  The Schur complement and its applications. volume 4.  Springer Science & Business Media.
  * Zhang [2016] Zhang, L., 2016.  Flood hazards impact on neighborhood house prices: A spatial quantile regression analysis.  Regional Science and Urban Economics 60, 12–19.