Regional compositional trajectories and structural change: A spatiotemporal multivariate autoregressive framework

Spatiotemporal Autoregressive Models for Areal Compositional Data

{\it Keywords: Spatiotemporal compositional data, economic compositions, spatial econometrics, real estate}

Introduction

Regional economic systems frequently undergo transformations driven by variations in the contributions of fundamental elements, such as the types of housing or industrial sectors prevalent within them. Formally defined, these fundamental components denote interrelated proportions that collectively form a complete entity, which must total a constant within each regional context. This includes situations like the allocation of economic activities among various sectors or the categorisation of property transactions by their types, and acknowledging this intrinsic dependence structure is crucial when analysing such datasets. Compositional data analysis, as outlined by MR1873662, offers a comprehensive methodological framework that facilitates the examination of such constrained relational data. Nonetheless, the analysis of intricate regional economic systems, where these fundamental components significantly influence, remains challenging. Aside from utilising traditional compositional data analysis techniques doi:10.1002/9781119976462 for the non-spatial components of the data, it is valuable to examine the spatial distribution characteristics of these compositions. A principal assumption in this examination is the influence of geographical proximity on local data, consistent with Tobler's First Law of Geography tobler\footnote{To quote Sir Ronald A. Fisher: “patches in close proximity are commonly more alike, as judged by the yield of crops, than those which are further apart” fisher1935design.}. In addition to the spatial dependence emanating from neighbouring areas, compositional data are frequently recorded multiple times over a period. In such environments, it is necessary to analyse both spatial interactions amidst the compositions and their temporal autocorrelation alongside dynamic progression. However, even with the increased accessibility of complex non-scalar spatio-temporal data and a pronounced need for spatio-temporal areal regression models, there exists a noticeable absence of appropriate model specifications for composition-valued information distributed across regions, i.e., areal spatial data. Although some literature documents applications of Whittle's simultaneous autoregressive whittle and Besag’s or Mardia's (multivariate) conditional autoregressive models Besag:74, mardia in spatial compositional datasets, extensions to more elaborate data structures are rare. This paper seeks to bridge this gap by presenting a multivariate simultaneous autoregressive model where the response at each site and temporal instance is a composition-valued data point. Methodological advancements, driven by regional economic dynamics data, are exemplified through case studies on the spatiotemporal composition of the Berlin real estate market and diverse economic sectors in Spain \textcolor{black}{(presented in the Appendix)}.

Concentrating on compositions of benthic species with exclusively non-zero parts, Billheimer1997 implemented an additive log ratio transformation (alr) 10.5555/17272 to position the data within Euclidean space and utilised a multivariate spatial conditional autoregressive model. Leininger2013 proposed a multivariate conditional autoregressive model to account for spatial configurations in land use and cover amidst multiple zero components. Beyond formulating the aforementioned conditional model, literature also addresses hierarchical PIRZAMANBEIN201814, simultaneous autoregressive Nguyen03042021, and spatial lag quantile regression models Zhao11122024. Furthermore, Yoshida2018 and Laurent2023 explored spatial relationships and covariate effects within compositional spatial area models, correspondingly. \textcolor{black}{Turning to compositional data dynamics, Paciorek2009 develops a Bayesian hierarchical spatio-temporal model to reconstruct historical forest composition from fossil pollen, Mastrantonio2019 models smooth temporal changes in compositions via a logistic–Gaussian process, and Kettunen2024 combines latent Gaussian processes with a Dirichlet–multinomial observation model to capture competition for space.}

Yet, to our best knowledge, no models exist to elucidate spatio-temporal dynamics in compositional areal data. Motivated by urban economic applications, this paper introduces a multivariate simultaneous autoregressive model tailored for spatial areal data, where the focus for each spatial site is a sequence of composition-valued quantities tracked over consecutive time steps. The proposed modelling framework, targeting composition-valued areal data observed over time, merges a multivariate simultaneous autoregressive structure, \textcolor{black}{building upon recent multivariate spatiotemporal autoregressive heteroscedasticity models by otto2024multivariate and the multivariate simultaneous autoregressive model by Kelejian1998,} with isometric log-ratio transformations to address compositional constraints.

The methodology is motivated by constrained regional economic dynamics, providing detailed case studies on Berlin's housing market compositions and sectoral compositions within Spanish municipalities (Section (ref)). Estimation utilises quasi-maximum likelihood, with proof of consistency and asymptotic normality as spatial and temporal domains expand (Section (ref)). A simulation study assesses the estimator's finite-sample characteristics (Section (ref)), and empirical findings (Section (ref)) illustrate the model’s capability to represent interpretable spatiotemporal trends in real-world constrained multivariate systems. Discussion on the findings and potential advancements concludes the paper (Section (ref)).

Cases Studies of Constrained Regional Economic Dynamics

This research is based on an in-depth analysis of two substantial case studies, which serve to illustrate the complexities and potentials of modelling regional economic dynamics under constraints. The first case study scrutinises the monthly composition of real estate transactions within Berlin's 24 different postcode regions over a span from 1995 to 2015. This analysis concentrates on the proportions of transactions differing among condominiums, developed parcels of land, and undeveloped land within the urban real estate market. This dataset provides invaluable insights into the temporal evolution of distinct sectors within the property market confined to a single city. Furthermore, fluctuations in one market segment, such as a rise in condominium sales, are assessed for their potential impact on the demand dynamics for developed or undeveloped properties in nearby areas. The foremost advantage of this study is the capacity to pinpoint spatial dependencies, indicating how property trends in one postal region might have implications for neighbouring ones, thereby highlighting the intrinsic connectedness of the urban real estate market. Figure (ref) presents comprehensive descriptive overviews of this data set. Analysis of the Berlin data reveals a gradual transformation in the intra-urban housing market composition, alongside indications of moderate spatial autocorrelation among proximate postcode regions.

The second case study delves into annual data collected from 2012 through 2021, focusing on the allocation of local enterprises across three primary economic sectors within 2,793 municipalities in Spain. These sectors include {services} (encompassing communications, financial and insurance services, administrative and support services, as well as educational, health and social services, and arts, recreation, and entertainment activities), {industry} (comprising extractive processes, manufacturing industries, energy and water provision, sanitation, and waste management operations), and {construction}. \textcolor{black}{The detailed results of the second case study can be found in the Appendix.} The two settings examined herein represent disparate spatiotemporal scales. The data from Berlin offers a high level of temporal granularity across a moderately large geographical area, whereas the dataset pertaining to Spain is characterised by extensive cross-sectional breadth, albeit with a relatively short temporal span.

These empirical patterns motivate the need for a flexible and interpretable model that jointly captures (i) the constrained nature of the data, (ii) temporal dependence within regions, (iii) spatial dependence across regions, and (iv) cross-component interactions between parts of the composition. In the next section, we introduce a spatiotemporal multivariate autoregressive model that meets these requirements and provides a coherent inferential framework for compositional areal panel data.

Modelling Framework for Spatiotemporal Compositional Data

This section introduces the modelling framework for spatiotemporal areal data with compositional responses. We begin by formalising the structure of composition-valued observations on spatial units over time (Section (ref)), and subsequently present a multivariate spatiotemporal autoregressive model that accounts for both spatial and temporal dependence while respecting compositional constraints (Section (ref)).

Our approach builds upon recent advances in spatial econometric modelling, in particular the multivariate simultaneous autoregressive (MSAR) model of yang2017identification, \textcolor{black}{which itself extends the multivariate simultaneous autoregressive model by Kelejian1998}. In that model, the outcome is a matrix of multivariate observations over areal units, and spatial dependence is captured through both within- and cross-component spatial autoregression:

where \(\mathbf{Y}\) is an \(n \times p\) matrix of outcomes for \(n\) spatial units and \(p\) response components, \(\mathbf{W}\) is a known spatial weight matrix, \(\mathbf{\Psi}\) is a \(p \times p\) matrix of spatial autoregressive coefficients, and \(\mathbf{\Pi}\) captures regression effects. The matrix \(\mathbf{V}\) contains spatially unstructured noise.

We extend this framework to incorporate temporal dependence, analogous to classical spatial dynamic panel data models Yu08. In the univariate case, a widely used model is given by

where \(\boldsymbol{Y}_t\) is the vector of observations at time \(t\), and the parameters \(\lambda\), \(\gamma\), and \(\rho\) control spatial contemporaneous, temporal, and spatiotemporal lag effects, respectively. Moreover, the series \(\{\boldsymbol{V}_t: t = 1,\ldots, T\}\) contains spatially and temporally independent error vectors.

Our goal is to extend these ideas to the case where \(\boldsymbol{Y}_t\) consists of composition-valued multivariate responses. This requires a formulation that captures multivariate spatial and temporal dependencies, while ensuring that model outputs respect the simplex geometry inherent to compositions. To this end, we introduce a spatiotemporal MSAR model that operates on isometrically log-ratio-transformed compositions, and derive an associated quasi-maximum likelihood estimator suitable for high-dimensional panels.

Compositional areal data in space-time

Extending compositional data to spatiotemporal areal processes let $\mathbf{Z}_t=\lbrace \boldsymbol{Z}_{t}(\boldsymbol{s}_n) \rbrace$ denote a collection of $nT$ constrained vectors which quantitatively describe the relative contribution of $D$ parts of some whole on $\mathcal{S}\times \mathcal{T}$ with $\mathcal{S} = \lbrace\boldsymbol{s}_1, \ldots, \boldsymbol{s}_n\rbrace$ denoting a countable set of $n$ (at least partially) interconnected areal units and $\mathcal{T} =\lbrace t_j\rbrace^T_{j=1}$ a set of $T$ distinct equidistant steps in time. In general, we assume that the areal units in $\mathcal{S}$ and the compositional $D$ parts remain consistent across all temporal instances $t\in \mathcal{T}\subset \mathds{R}_+$ such that at each location $\boldsymbol{s}_j$ and each $t\in \mathcal{T}$ $Z_i=Z_t(\boldsymbol{s}_i)$ is an element in the simplex $\mathds{S}^{D}$, \[ \operatorname{\mathds{S}}^{D}=\left\{Z =(Z_{1},Z_{2},\dots ,Z_{D})^{\top}\in \mathds {R}^{D}\,\left|\,Z_{l}\geq 0,l=1,2,\dots ,D;\sum_{l=1}^{D}Z_{l}=\kappa \right.\right\}, \] consisting of the same number of $D$ non-negative mutually dependent components that \textcolor{black}{sum} to a constant $\kappa\in\mathds{R}$. \textcolor{black}{We note that in the original definition of $\operatorname{\mathds{S}}^{D}$, all $D$ parts are strictly assumed to be positive to avoid problems with zeros in log-ratio transformations. More recently, however, Tsagris2016 introduced so-called $\alpha$-transformations which extend classic log-ratio transformations to work with zeros in some compositional parts}. Consequently, for the specific time instance $t\in \mathcal{T}$ $\boldsymbol{Z}_t$ constitutes an $n \times D$ dimensional matrix that encapsulates the entirety of compositions across all $n$ spatial areal units for the specific time instance $t\in \mathcal{T}$. We note that any vector $\widetilde{Z}_i$ of $D$ real positive components can always be translated into composition-valued information by applying the closure operator $\operatorname{cls}(\widetilde{Z}_i)$ with \[ Z_i=\operatorname{cls}(\widetilde{Z}_i)= \left(\frac{\tilde{Z}_{i1}}{\sum_{l=1}^D \tilde{Z}_{il}},\ldots,\frac{\tilde{Z}_{iD}}{\sum_{l=1}^D \tilde{Z}_{il}}\right)^\top. \] Making use of the so-called Aitchison geometry MR1873662, $\mathds{S}^D$ together with the perturbation and powering operations $\oplus$ and $\odot$ where \[ Z_i\oplus Z_j=\operatorname{cls}(Z_{i1}Z_{j1},Z_{i2}Z_{j2},\ldots, Z_{iD}Z_{jD}) \] and \[ \xi\odot Z_i=\operatorname{cls}(Z_{i1}^{\xi},Z_{i2}^{\xi},\ldots, Z_{iD}^{\xi}) \] with $Z_i,Z_j\in \operatorname{\mathds{S}}^D$ and $\xi\in\mathds{R}$ can be turned into a Hilbert space with the Aitchison inner product

Aitchison norm \textcolor{black}{$\|Z_i\|_A = \sqrt{\langle Z_i, Z_i\rangle_A}$} and the associated Aitchison \textcolor{black}{metric} \[ d_A(Z_i,Z_j)=\Vert Z_i\ominus Z_j \Vert_A \] where $\Vert Z_i\ominus Z_j \Vert_A=Z_i\oplus((-1)\odot Z_j)$ is the \textcolor{black}{negative perturbation} operation. Instead of working on the simplex $\operatorname{\mathds{S}}^D$, it is often more convenient to work on $\mathds{R}^{\tilde{D}}$ by applying \textcolor{black}{a} map $\psi:\mathds{S}^D\to \mathds{R}^{\tilde{D}}, Z_i\mapsto \psi(Z_i)$ and performing statistical analysis methods in $\mathds{R}^{\tilde{D}}$ where ${\tilde{D}}$ is determined by the particular choice of $\psi$ doi:https://doi.org/10.1002/9781119976462.ch3. Denoting the $k$-th coordinate of the transformed spatial composition $\psi(Z_i)=(\psi_1(Z_i),\ldots, \psi_{\tilde{D}}(Z_i))^\top$ by $\psi_k(Z_i)$, the underlying idea is to express $Z_i\in \operatorname{\mathds{S}}^D$ in the form of a canonical basis function representation \[ Z_i=\bigoplus^{\tilde{D}}_{k=1}\psi_k(Z_i)\odot\mathbf{w}_k \] where $\mathbf{w}_k=\operatorname{cls}(\exp(\mathbf{d}_k)),~k=1,\ldots,{\tilde{D}}$, with $\mathbf{d}_k\in\mathds{R}^{\tilde{D}}$ denoting the unit vector associated with the $k$-th coordinate doi:10.1002/9781119976462. While different transformations like the additive log-ratio (alr) transformation 10.1093/biomet/67.2.261 have been discussed in the literature, here we only review the centred log-ratio transformation (clr) 10.1093/biomet/70.1.57 and isometric log-ratio transformation\textcolor{black} {s} (ilr) Egozcue2003 which establish an isometric isomorphism between $\operatorname{\mathds{S}}^D$ and $\mathds{R}^{\Tilde{D}}$ BillheimerEtAl2001, PawlowskyGlahn2001. The recognised isometric isomorphism bridging Aitchison geometry with Euclidean geometry ensures that the inner product in Aitchison space, as well as distances and metrics, corresponds precisely to their Euclidean counterparts once variables have undergone transformation.

Imposing a sum-to-zero constraint, the clr transformation with \textcolor{black}{$\operatorname{clr}_l(Z_i)=\log(Z_{il}/g(Z_i))$}, where $g(Z_i)$ denotes the geometric mean of the composition $Z_i$, maps the composition from the simplex to a hyperplane $\mathds{H}\subset\mathds{R}^{D}$ that is orthogonal to the vector of ones. Unlike the clr transformation that yields degenerate distributions and singular covariance matrices, the ilr transformation provides a map between $\mathds{S}^D$ and $\mathds{R}^{D-1}$ which corresponds to a class of orthonormal coordinate representations that are derived by applying the Gram-Schmidt procedure to an orthonormal basis $(\mathbf{e}_1, \mathbf{e}_2,\ldots,\mathbf{e}_{D-1})$ on the simplex $\mathds{S}^D$ \[ Z_i=\bigoplus_{k=1}^{D-1}\langle Z_i,\mathbf{e}_k\rangle_A\odot\mathbf{e}_k, \ \ \] yielding $\operatorname{ilr}(Z_i) = \left(\langle Z_i,\mathbf{e}_1\rangle_A,\langle Z_i,\mathbf{e}_2\rangle_A,\ldots,\langle Z_i,\mathbf{e}_{D-1}\rangle_A\right)$. Both the clr and ilr transformations are related to each other through \textcolor{black}{$\operatorname{ilr}(Z_i)=\mathbf{V}_\mathrm{D}\operatorname{clr}(Z_i)$} with $\mathbf{V}_\mathrm{D}$ denotes a $((D-1)\times D)$-dimensional Helmert matrix with rows $\mathbf{v}_j=\operatorname{clr}(\mathbf{e}_j), j=1,\ldots,D-1$ satisfying $\mathbf{V}_\mathrm{D}\mathbf{V}_\mathrm{D}^{\top}=\mathbf{I}_{\mathrm{D}-1}$ and $\mathbf{V}_\mathrm{D}^{\top}\mathbf{V}_\mathrm{D}=\mathbf{I}_{\mathrm{D}}-D^{-1}\mathds{1}_\mathrm{D}\mathds{1}^{\top}_\mathrm{D}$ where $\mathbf{I}_{\mathrm{D}}$ is the identity matrix of dimension $(D\times D)$, and $\mathds{1}_\mathrm{D}$ a $(D\times 1)$ vector of ones. When dealing with a composition comprising of two parts ($D=2$), the $\operatorname{ilr}$ transformation aligns with the logit function typically applied in logistic regression. In scenarios where the number of components exceeds two ($D>2$), a multitude of orthonormal basis systems are possible, all of which substantially influence how projected data is interpreted. Specific configurations of these bases include coordinate representations through balances Egozcue2005:Balances where each balancing element can be understood as the normalised log-ratio of the geometric centres of two defined groups. This concept is rooted in the sequential binary partition technique, which entails dividing the composition into two distinct sections. In another configuration called pivot coordinates Fiserova2011, Hron2017, the very first ilr coefficient mirrors the numerator of the initial clr coefficient, adjusted by a factor of $\sqrt{D/(D-1)}$, rendering it comprehensible as the log-ratio between a specific component and the geometric mean. Conversely, the interpretation of subsequent coefficients proves to be more intricate. To manage this complexity, literature has put forward generalised pivot coordinates, which utilise permuted compositions, alongside symmetric pivot coordinates sym:pivot, Hron2021 to aid in better understanding and application.

It is worth noting that the aforementioned formulation necessitates strict positivity across all $D$ compositional elements under investigation. However, this requirement has been addressed by extending both the centred log-ratio (clr) transformation and the isometric log-ratio (ilr) transformation, leading to the development of the centred $\alpha$ ($\operatorname{c\alpha t}$) and isometric $\alpha$ ($\operatorname{i\alpha t}$) transformations. These extensions, as introduced by Tsagris2016 and further elaborated by CLAROTTO2022100570, enable the inclusion of zero values in some of the $D$ components, thereby broadening the applicability of these transformations in practical scenarios.

Dynamic spatiotemporal autoregressive panel model

Consider a compositional spatiotemporal areal process $\mathbf{Z}_t$ as described in Section (ref) carrying information on $D$ compositional part for a discrete set of spatial units comprising $n$ distinct geographical locations $\{\boldsymbol{s}_1, \ldots, \boldsymbol{s}_n\}$ that remains consistent across all temporal instances $t = 1, \ldots, T$. Instead of $\mathbf{Z}_t$ and making use of the isometric isomorphism between the simplex and the Euclidean spaces, let $\mathbf{Y}_t$ denote the $\psi$-transformed process $\psi(\mathbf{Z}_t)$, such that for any temporal instance $t$, $\boldsymbol{Y}_t$ is a $n \times \tilde{D}$-dimensional matrix. In particular, to avoid any degenerative distribution and non-singular covariance matrices, we set $\psi=\operatorname{ilr}$ such that $\mathbf{Y}_t$ is a $n \times (D-1)$-dimensional matrix. Further, we assume that $\mathbf{Y}_t$ follows a multivariate spatiotemporal autoregressive process \textcolor{black}{with higher-order temporal lags, i.e.,}

where $\mathbf{E}_t = (\boldsymbol{\varepsilon}_{1,t}, \ldots, \boldsymbol{\varepsilon}_{D-1,t})$ is the $n \times (D-1)$-dimensional matrix of disturbances with independent and identically distributed random vectors $\boldsymbol{\varepsilon}_{j,t} = (\varepsilon_{j,t}(\boldsymbol{s}_1), \ldots, \varepsilon_{j,t}(\boldsymbol{s}_n))^\top$ with $\mathds{E}(\boldsymbol{\varepsilon}_{j,t}) = \boldsymbol{0}$ and $\mathds{C}\text{ov}(\boldsymbol{\varepsilon}_{j,t}) = \hat{\sigma}^2 \mathbf{I}_n$ for all $j = 1,\ldots, D-1$ and $t = 1,\ldots, T$, $\mathbf{I}_n$ is the $n$-dimensional identity matrix and $\mathbf{X}_{t,i}$ is an $n \times (D-1)$-dimensional matrix of the $i$-th exogenous regressors at time point $t$ which enter the regression model with slope coefficients $\boldsymbol{\beta}_i = (\beta_{i,1}, \ldots, \beta_{i,D-1})^\top$. In what follows, we will make use of matrix notation to represent $\boldsymbol{\beta}_i$ in compact form by $\mathbf{B} = (\beta_{ij})_{i = 1, ..., q, j = 1,...,D-1}$. Moreover, the model includes a spatial autoregressive term with an $n \times n$-dimensional matrix $\mathbf{W}$ of spatial weights \textcolor{black}{and an unknown coefficient matrix $\mathbf{\Psi}$ to be estimated}. This matrix defines the spatial proximity structure, i.e., which locations are considered to be adjacent and can thereby influence each other. In practice, this matrix is assumed to be known and determined by the underlying geography. \textcolor{black}{Moreover, the model includes $1 \leq \tau_1 < \ldots < \tau_L$ temporal lags and $\{\mathbf{\Pi}_1, \ldots, \mathbf{\Pi}_L\}$ are the corresponding $(D-1)$-dimensional square coefficient matrices.} The cross-component spatial effects are represented by the off-diagonal elements of $\mathbf{\Psi}$, and the temporally lagged cross-component effects are given by the off-diagonal elements of \textcolor{black}{$\{\mathbf{\Pi}_1, \ldots, \mathbf{\Pi}_L\}$}. In addition, the component-wise spatial and temporal autoregressive effects are summarised by the diagonal entries of $\mathbf{\Psi}$ and \textcolor{black}{$\{\mathbf{\Pi}_1, \ldots, \mathbf{\Pi}_L\}$}, respectively.

The model will be estimated using a quasi maximum likelihood (QML) approach; combining the results of otto2024multivariate on multivariate spatiotemporal GARCH models, yang2017identification, who derived identifiability conditions the consistency and asymptotic normality of a QML estimator for multivariate \textcolor{black}{purely} spatial data, and Yu08, who derived asymptotic results for a QML estimator for spatiotemporal but univariate processes when both $n$ and $T$ are large. Using the $\operatorname{vec}$-operator for the vectorisation of a matrix, we can rewrite (ref) to get the vectorised form

The Kronecker product is denoted by $\otimes$. Interestingly, using such vec-representation, one can see that the multivariate spatiotemporal autoregressive model with $n$ spatial units is a special case of a (univariate) $n(D-1)$-dimensional spatiotemporal autoregressive model with a weight matrix $\mathbf{\Psi}^\top \otimes \mathbf{W}$ (i.e., having $n(D-1)$ artificial spatial units).

\textcolor{black}{For an easier notation, let $\ddot{\boldsymbol{Y}}_t = \operatorname{vec}(\mathbf{Y}_t)$, $\boldsymbol{u}_t = \operatorname{vec}(\mathbf{E}_t)$, $\mathbf{A}_\ell(\mathbf{\Pi}_\ell)=\mathbf{\Pi}_\ell^\top\otimes \mathbf{I}_n$, $\mathbf{S}_{n(D-1)}(\mathbf{\Psi}) = \mathbf{I}_{n(D-1)} - \mathbf{\Psi}^\top\otimes \mathbf{W}_n$, and $\boldsymbol{x}_t(\mathbf{B}) = \operatorname{vec}\!\Big(\sum_{i=1}^q \mathbf{X}_{t,i}\boldsymbol{\beta}_i\Big)$. Then, we can express the vectorised model in reduced form as

which provides the residual map for $\vartheta=(\operatorname{vec}(\mathbf{B})^\top, \operatorname{vec}(\mathbf{\Psi})^\top, \operatorname{vec}(\mathbf{\Pi}_1)^\top,\ldots,\operatorname{vec}(\mathbf{\Pi}_L)^\top)^\top$, \[ \boldsymbol{r}_t(\vartheta) = \mathbf{S}_{n(D-1)}(\mathbf{\Psi}) \ddot{\boldsymbol{Y}}_{t} -\boldsymbol{x}_t(\mathbf{B}) -\sum_{\ell=1}^{L}\mathbf{A}_\ell(\mathbf{\Pi}_\ell)\,\ddot{\boldsymbol{Y}}_{t-\tau_\ell}. \] Recall $\mathbf{S}_{n(D-1)}= \mathbf{I}_{n(D-1)}- \mathbf{\Psi} ^\top\!\otimes \mathbf{W}$. A convenient sufficient condition for invertibility is $\rho(\mathbf{\Psi}^\top\!\otimes \mathbf{W})<1$, where $\rho(\cdot)$ denotes the spectral radius. Using $\rho(\mathbf{A}\otimes \mathbf{B})=\rho(\mathbf{A})\rho(\mathbf{B})$, this is satisfied whenever $\rho(\mathbf{\Psi})\rho(\mathbf{W})<1$. For row-standardised spatial weight matrices, $\rho(\mathbf{W})\le 1$, implying that $\rho(\mathbf{\Psi})<1$ ensures existence of the spatial multiplier $\mathbf{S}_{n(D-1)}^{-1}=\sum_{k=0}^{\infty}(\mathbf{\Psi}^\top\!\otimes \mathbf{W})^k$. In the empirical applications, we verify this condition numerically by checking that the eigenvalues of $\mathbf{S}_{n(D-1)}$ are bounded away from zero.}

\textcolor{black}{Under the data-generating parameter $\vartheta_0 = (\operatorname{vec}(\mathbf{B}_0)^\top, \operatorname{vec}(\mathbf{\Psi}_0)^\top, \operatorname{vec}(\mathbf{\Pi}_{1,0})^\top, \ldots, \operatorname{vec}(\mathbf{\Pi}_{\ell,0})^\top)^\top$, we have $\boldsymbol{r}_t(\vartheta_0)=\boldsymbol{u}_t$. Let $\tau_{\max} = \max_{\ell=1,\ldots,L}\tau_\ell$ and $T^\star = T-\tau_{\max}$ and $\sigma^2 > 0$ denote the (unknown) innovation variance parameter. The Gaussian quasi log-likelihood (up to additive constants independent of $(\vartheta,\sigma^2)$) is given by {

} Let $(\hat\vartheta_{nT},\hat\sigma^2_{nT})$ maximise (ref) over the parameter space $\Theta\times[\underline\sigma^2,\overline\sigma^2]$, where $0<\underline\sigma^2<\sigma_0^2<\overline\sigma^2<\infty$.}

To establish the consistency of the quasi-maximum likelihood (QML) estimator for our spatiotemporal autoregressive compositional panel data model, we require a set of assumptions. These assumptions govern the behaviour of the error process, the parameter space, some standard assumptions on the spatial weight matrix, and the asymptotic framework under which consistency can be established.

Assumption (ref) ensures that \textcolor{black}{innovation process has well-defined, finite second and higher-order moments and is serially independent over time. The finite $(4+\eta)$-moment condition ensures that quadratic forms of the residuals satisfy a law of large numbers and a central limit theorem, which is required for consistency and asymptotic normality of the QML estimator.} From a practical perspective, \textcolor{black}{this assumption is satisfied} in many empirical applications where extreme shocks are \textcolor{black}{rare}. In case of many outliers or (extremely) heavy-tailed distributions, the estimation results \textcolor{black}{and inference} could be distorted. \textcolor{black}{If the regressors are stochastic rather than deterministic, the following additional assumption is required; in the case of deterministic and uniformly bounded regressors, this assumption can be omitted.}

\textcolor{black}{To ensure that the spatiotemporal autoregressive system in (ref) is well-defined and admits a stable causal solution (uniformly over the admissible parameters),} we impose the following stability and compactness condition; \textcolor{black}{the innovation variance is restricted to $\sigma^2\in[\underline\sigma^2,\overline\sigma^2]$ as in (ref).}

\textcolor{black}{For $r=1,\ldots,\tau_{\max}$, let \[ \mathbf{\Phi}_r(\vartheta)=

\] Further, the homogeneous part of the reduced-form dynamics can be written as \[ \ddot{\boldsymbol{Y}}_t=\sum_{r=1}^{\tau_{\max}}\mathbf{\Phi}_r(\vartheta)\ddot{\boldsymbol{Y}}_{t-r}. \] Let $\boldsymbol{Z}_t= (\ddot{\boldsymbol{Y}}_t^\top, \ddot{\boldsymbol{Y}}_{t-1}^\top, \ldots, \ddot{\boldsymbol{Y}}_{t-\tau_{\max}+1}^\top)^\top$. Then, the recursion can be written as \[ \boldsymbol{Z}_t=\mathbf{F}(\vartheta)\boldsymbol{Z}_{t-1}, \] where $\mathbf{F}(\vartheta)$ is the block companion matrix \[ \mathbf{F}(\vartheta)=

. \] }

Since our model includes a spatial autoregressive component, we must also impose regularity conditions on the spatial weight matrix \( \mathbf{W} \), as formulated in Assumption (ref).

These conditions ensure that spatial dependence remains well-behaved as $n$ increases, thereby maintaining the stability of the estimation procedure. In practical applications, this assumption is typically met with standard choices of the spatial weight matrix, such as those based on geographic adjacency or distance-based kernel functions, which naturally enforce decay in spatial interactions. This prevents scenarios where each location has an unbounded number of influential neighbours, which could otherwise result in singularities in the estimation process. \textcolor{black}{To rule out collinearity between regressors and lagged responses and to ensure identification of $(\mathbf{B},\mathbf{\Pi})$, we impose the following full-rank condition.}

Finally, we need an appropriate asymptotic framework to derive the consistency result, as provided in Assumption (ref).

Given the above assumptions, we can now establish the main theoretical result showing the consistency of the QML estimator in Theorem (ref). This result ensures that our estimates converge to the true parameter values as the length of the panel increases, providing the foundation for inference in the spatiotemporal autoregressive compositional panel model. The imposed assumptions collectively guarantee that the estimation problem is well-posed, the likelihood function behaves regularly, and the model remains stable under increasing sample sizes, leading to a consistent QML estimator.

The proof of the Theorem is provided in the Appendix. \textcolor{black}{While Theorem (ref) establishes convergence in probability, asymptotic normality of $(\hat\vartheta_{nT},\hat\sigma^2_{nT})$ requires additional smoothness and moment conditions ensuring that a Taylor expansion of the score is valid and that the (normalised) score satisfies a central limit theorem under the joint asymptotic regime $T\to\infty$ with $n=n(T)$ non-decreasing. For that reason, we need some additional growth constraints on $n$.}

\textcolor{black}{

}

\textcolor{black}{ Assumption (ref) introduces only the additional conditions required to strengthen consistency to asymptotic normality. Smoothness of the Gaussian quasi log-likelihood and local uniform boundedness of its derivatives follow from the structural properties of the model and the earlier assumptions. In particular, the residual map is affine in the parameters and the log-determinant term is smooth whenever $\mathbf{S}_{n(D-1)}(\mathbf{\Psi})$ is invertible, which holds under Assumption (ref). Combined with compactness of the parameter space and the moment bounds on the innovations, this implies that the likelihood is twice continuously differentiable in a neighbourhood of $(\vartheta_0,\sigma_0^2)$ and that its Hessian converges locally uniformly in probability to its expectation. The asymptotic theory is therefore driven by temporal aggregation. The temporal dimension satisfies $T\to\infty$, while the spatial dimension may remain fixed or increase with $T$. The mild growth restriction $n(T)/T\to0$ ensures that the cross-sectional dimension does not expand too quickly relative to the temporal domain. Since each single-period likelihood contribution is already normalised by the cross-sectional dimension $n(D-1)$, this condition guarantees that the score satisfies a central limit theorem and that the Hessian can be replaced by its population counterpart, yielding the usual $\sqrt{T}$-rate for the quasi-maximum likelihood estimator. }

\textcolor{black}{

}

\textcolor{black}{The proof of the theorem is provided in the Appendix. As a consequence of Theorem (ref), asymptotically valid} standard errors can be obtained from the \textcolor{black}{estimated asymptotic covariance matrix. In general, this is given by the sandwich form $\mathcal{J}^{-1}\mathcal{I}\mathcal{J}^{-1}$, where $\mathcal{J}$ and $\mathcal{I}$ are consistently estimated by their sample analogues. If the Gaussian likelihood is correctly specified, the information identity holds and the covariance matrix simplifies to $\mathcal{J}^{-1}$, which can be estimated} by the inverse of the observed Hessian of the log-likelihood evaluated at the QML estimator.

\textcolor{black}{In the present spatiotemporal panel setting, the quasi log-likelihood can be written as \[ \ell_{nT}(\vartheta,\sigma^2) = \frac{1}{T}\sum_{t=1}^T \ell_t(\vartheta,\sigma^2), \] where $\ell_t$ denotes the normalised contribution of time point $t$. Under correct specification of the temporal dynamics and Assumption (ref), the score contributions $\nabla_{(\vartheta,\sigma^2)} \ell_t(\vartheta_0,\sigma_0^2)$ form a martingale difference sequence and are serially uncorrelated across $t$. In this case, the matrix $\mathcal{J}$ can be consistently estimated by the negative time-average of the observed Hessian, \[ \widehat{\mathcal{J}} = - \frac{1}{T} \sum_{t=1}^T \nabla^2_{(\vartheta,\sigma^2)} \ell_t(\hat\vartheta_{nT},\hat\sigma^2_{nT}), \] and $\mathcal{I}$ by the sample covariance of the score contributions, \[ \widehat{\mathcal{I}} = \frac{1}{T} \sum_{t=1}^T \hat{s}_t \hat{s}_t^\top, \qquad \hat{s}_t = \nabla_{(\vartheta,\sigma^2)} \ell_t(\hat\vartheta_{nT},\hat\sigma^2_{nT}). \] }

\textcolor{black}{If the temporal dependence is not fully captured by the specified lag structure, the score process may exhibit residual serial correlation. In that case, $\mathcal{I}$ should be estimated using a long-run variance estimator, for example a heteroscedasticity and autocorrelation consistent (HAC) estimator with suitable kernel and bandwidth. Such estimation requires a sufficiently long time dimension $T$. When the dynamic specification is adequate and $T$ is moderate, the simpler time-averaged Hessian and score covariance estimators are typically preferred.}

\textcolor{black}{Model interpretation and effects on the simplex}

\textcolor{black}{The parameter matrices $\mathbf{\Psi}$ and $\mathbf{\Pi}$ describe linear spatiotemporal dependence in the transformed Euclidean space induced by the mapping $\psi:\mathds{S}^D \rightarrow \mathds{R}^{\tilde D}$. For $\psi=\operatorname{ilr}$, we have $\tilde D = D-1$. Due to the non-linear nature of the inverse transformation $Z=\psi^{-1}(Y)$, these parameters do not correspond to constant linear effects on the original compositional scale. Instead, the implied dynamics on the simplex are inherently non-linear and state dependent.}

\textcolor{black}{Interpretation on the compositional scale is therefore based on local marginal effects. Let $Z_t(\boldsymbol{s}_i)=\psi^{-1}(Y_t(\boldsymbol{s}_i))$ denote a reference composition at spatial unit $\boldsymbol{s}_i$ and time $t$, and define the Jacobian matrix of the inverse transformation as \[ \mathbf{J}_\psi\!\left(Y_t(\boldsymbol{s}_i)\right) = \frac{\partial \psi^{-1}(Y_t(\boldsymbol{s}_i))}{\partial Y_t(\boldsymbol{s}_i)^\top} \in \mathds{R}^{D\times (D-1)}. \] Conditional on a reference composition $Z_t(\boldsymbol{s}_i)$, local effects of past outcomes are characterised lag-wise. Specifically, the local marginal effect of $Y_{t-\tau_\ell}(\boldsymbol{s}_i)$ on $Z_t(\boldsymbol{s}_i)=\psi^{-1}(Y_t(\boldsymbol{s}_i))$ is governed by $\mathbf{J}_\psi(Y_t(\boldsymbol{s}_i))\,\mathbf{\Pi}_\ell^\top$, where the Jacobian is evaluated at the contemporaneous state $Y_t(\boldsymbol{s}_i)$ and $\mathbf{\Pi}_\ell$ governs the propagation of past innovations from $t-\tau_\ell$ to $t$.}

\textcolor{black}{Contemporaneous spatial dependence propagates through the reduced form representation \[ \small{\mathrm{vec}(\mathbf{Y}_t) = \left(\mathbf{I}_{n(D-1)}-\mathbf{\Psi}^\top\!\otimes \mathbf{W}\right)^{-1} \left[ \operatorname{vec}\!\left(\sum_{i=1}^{q}\mathbf{X}_{t,i}\boldsymbol{\beta}_i\right) + \sum_{\ell=1}^{L}(\mathbf{\Pi}_\ell^\top\otimes \mathbf{I}_n)\,\mathrm{vec}(\mathbf{Y}_{t-\tau_\ell}) + \operatorname{vec}(\mathbf{E}_t) \right],} \] where the inverse matrix captures the spatial transmission of shocks across units and components. Accordingly, the impact of a perturbation at spatial unit $\boldsymbol{s}_j$ on unit $\boldsymbol{s}_i$ is governed by the $(i,j)$ block of the reduced form matrix, and its effect on the original compositions is obtained by post-multiplication with $\mathbf{J}_\psi(Y_t(\boldsymbol{s}_i))$.}

\textcolor{black}{This construction yields natural definitions of direct effects (own-unit impacts, $i=j$) and indirect or spillover effects (cross-unit impacts, $i\neq j$) on the simplex. All effects are evaluated locally at meaningful reference compositions and are expressed in terms of changes in compositional shares. While the parameter matrices $\mathbf{\Psi}$ and $\mathbf{\Pi}$ depend on the chosen ilr basis, the induced effects on the simplex are invariant under orthonormal reparametrisations, providing an interpretable and basis-independent link between the linear spatiotemporal model in transformed space and the non-linear geometry of composition-valued data. Similar arguments for emphasising simplex-intrinsic interpretations rather than coordinate-dependent coefficients have recently been discussed in the compositional regression literature DARGEL2024107945.}

\textcolor{black}{The interpretation of the regression coefficients $\mathbf{B}$ associated with the exogenous covariates follows similar principles. Since the covariates enter the model linearly in ilr-transformed space, each coefficient $\beta_{i,k}$ represents a constant marginal effect of the corresponding regressor on the $k$-th ilr coordinate of the composition in the structural form of the model.} \textcolor{black}{However, due to the presence of contemporaneous spatial dependence, covariate effects propagate through space according to the reduced form of the model. In particular, a change in an exogenous regressor at spatial unit $\boldsymbol{s}_j$ affects $\mathbf{Y}_t$ through the spatial multiplier \[ \left(\mathbf{I}_{n(D-1)}-\mathbf{\Psi}^\top\!\otimes \mathbf{W}\right)^{-1}, \] giving rise to both direct effects (on the same unit) and indirect or spillover effects (on other units).} \textcolor{black}{On the original compositional scale, the marginal effect of a covariate is obtained by mapping the reduced-form impact back to the simplex via the inverse ilr transformation. Locally, at a reference composition $Z_t(\boldsymbol{s}_i)$, the effect of a covariate perturbation is therefore characterised by \[ \mathbf{J}_\psi(Y_t(\boldsymbol{s}_i)) \left[ \left(\mathbf{I}_{n(D-1)}-\mathbf{\Psi}^\top\!\otimes \mathbf{W}\right)^{-1} \operatorname{vec}(\mathbf{X}_{t,i}\boldsymbol{\beta}_i) \right], \] yielding changes in compositional shares that respect the unit-sum constraint. This construction allows covariate effects to be decomposed into direct and indirect components on the simplex, in close analogy to effect decomposition in spatial econometric models, while accounting for the non-linear geometry of compositional data induced by the inverse log-ratio transformation.}

Monte-Carlo Simulations

In the following section, we will present the results obtained from a series of simulations assessing the \textcolor{black}{performance of the estimators} in finite sample sizes. Therefore, we simulated a spatiotemporal model for 9 different sample sizes with an increasing size of the spatial field $n \in \{4^2, 6^2, 8^2\}$ and the length of the time series $T \in \{30, 80, 160\}$. All processes were simulated on a quadratic two-dimensional grid $\{\mathds{Z}^2: 0 < s_1, s_2 \leq \sqrt{n}\}$ and all cells were considered to be adjacent if they share a common border, i.e., $\mathbf{W}$ was considered to be a Queen's contiguity matrix. Moreover, the weight matrix was row-standardised, and the simulation study was conducted with \textcolor{black}{$1064$} replications. \textcolor{black}{In each replication, covariates (including an intercept) are generated independently and innovations are drawn either from a Gaussian distribution or from a Gaussian mixture (0.95 of a standard normal and 0.05 of a zero-mean Gaussian with variance 25) to induce deviations from normality. The regression coefficient matrix is fixed across settings as $\mathbf B=\bigl((1,2)^\top,(-2,1)^\top,(3,-2)^\top\bigr)^\top$ (i.e., $k=3$ and $p=2$). Moreover, we considered two different settings for the degree of spatiotemporal dependence.}

\textcolor{black}{ Setting A (moderate overall dependence). The data-generating process uses a single temporal lag, $\tau_{\text{dgp}}=\{1\}$. Instantaneous spatial dependence is strong and diagonal, \[ \mathbf\Psi=

, \qquad \mathbf\Pi_{1}=

. \] For the lattice sizes considered, the resulting spatiotemporal transition operator has spectral radius of about $0.6$, implying moderate persistence and a comparatively well-conditioned estimation problem. For setting A, we included two shorter time horizons of 5 and 10 to assess the performance for very short time series.}

\textcolor{black}{ Setting B (high persistence, multi-lag temporal dynamics). The DGP includes two temporal lags, $\tau_{\text{dgp}}=\{1,12\}$. Spatial dependence is slightly weaker but includes cross-coordinate effects, \[ \mathbf\Psi=

, \qquad \mathbf\Pi_{1}=

, \qquad \mathbf\Pi_{12}=0.3\,\mathbf\Pi_{1}. \] Here the implied transition operator has spectral radius close to $0.9$, producing substantially higher persistence and slower effective information accumulation in $T$. This setting is designed to stress finite-sample inference, in particular standard error estimation and coverage under both Gaussian and non-Gaussian innovations. The smallest time series length for setting B was chosen to be 30, leading to an effective length of 18 time points when removing the first 12 time points due to the multi-lag structure.}

\textcolor{black}{ Figures (ref)--(ref) reveal three broad patterns that are consistent with the dependence regimes induced by Settings A and B. First, the average bias (averaged across replications and matrix elements for each component) is small throughout and decreases with $T$ (Figure (ref)), with the most visible finite-sample bias occurring for the regression block in Setting B at small $T$. This is expected: the multi-lag temporal component in Setting B implies substantially higher persistence (spectral radius close to $0.9$), so the effective information accumulation in time is slower and the estimator behaves more like a near-unit-root panel. In contrast, Setting A (spectral radius around $0.6$) is markedly less persistent, and the estimator is close to unbiased already for moderate $T$. Second, (average) RMSE decreases monotonically in both $n$ and $T$ (Figure (ref)). The steep RMSE decline in the regression block is mainly described by improved estimation of the conditional mean as the panel grows, while the relatively flat curves for $\mathbf{\Psi}$ and the $\mathbf{\Pi}$ blocks reflect that these parameters are already well-identified at the chosen signal-to-noise ratio. The slower convergence in Setting B is again consistent with high persistence: the process spends longer in correlated “runs”, reducing the effective sample size in time.}

\textcolor{black}{ The coverage results shown in Figure (ref) reflect the different dependence regimes considered in the two simulation settings. In Setting A, coverage is below the nominal $0.95$ level for small time dimensions, particularly under mixture innovations, but increases steadily with $T$ and approaches the nominal level for larger samples. Differences between the naïve, sandwich (OPG), and HAC variance estimators are generally small in this regime for $T > 30$. In Setting B, which features stronger dynamic persistence, coverage is similar to setting A, slightly lower though. Robust variance estimators improve coverage in these cases, but the need a sufficient length of the time series. As $T$ increases, coverage stabilises across all parameter blocks and innovation types, and the differences between the three estimators become negligible. Overall, the results indicate that inference improves with the temporal dimension and that robust covariance estimators mainly matter in the more persistent regime and under non-Gaussian innovations. }

Empirical Results

In this section, we aim to thoroughly examine the practical implementation of the spatiotemporal autoregressive model as applied to compositional data \textcolor{black}{based on the real-estate market example from above, i.e.,} composition of real estate transactions in Berlin over a 20-year span from 1994 to 2014. \textcolor{black}{In addition to this case study, we discuss the results of a second empirical example in the Appendix, Section (ref). In contrast to the real-estate transaction compositions, this second case has a much higher spatial resolution, but smaller time horizon.}

\textcolor{black}{Parameter estimates in ILR space}

\textcolor{black}{Our} empirical investigation examines the spatial and temporal dynamics of real estate transactions within the Berlin metropolitan area, concentrating on three distinct property categories: condominiums, developed plots, and undeveloped plots. This investigation is supported by a dataset containing monthly transaction compositions for all postcode regions in Berlin, covering the period from 1995 to 2015, encompassing $n = 24$ 3-digit postcode areas and $T = 240$ time periods. Through this analysis, we are able to analyse the interactions over time and space, as well as the market share distribution among various property types within the city's internal real estate environment. Such analysis holds particular significance for comprehending the temporal progression of different segments within the urban property market. It illuminates how an increase in the sales of one property type, such as condominiums, might affect the demand for developed or undeveloped plots in surrounding localities. Furthermore, it facilitates the identification of spatial linkages, where patterns in one postcode region may exert influence on neighbouring areas, demonstrating the interconnected nature of the real estate market throughout the city. Figure (ref) presents descriptive plots that visualise the spatiotemporal dynamics of the composition data. In the subsequent discussion, extending from the brief introduction in Section (ref), we provide additional details for the first case study. Initially, the two left-side plots in the first row demonstrate the compositions for two different years—at the onset and conclusion of the time series—on the simplex. The observation colours correspond to their spatial position, with the map on the right functioning as a colour legend. A subtle spatial dependence is evident, as similar colours tend to cluster nearby on the simplex. Moving to the second row, two simplices illustrate the temporal progression of selected regions, marked on the map by a yellow and a green cross. Both regions exhibit a trend over time shifting towards different property types. Specifically, the green region, represented by postcode area 126xx, transitions from predominantly undeveloped land transactions to a focus on developed plots and condominiums. Finally, the bottom row showcases the time series plots of compositional evolution for these two regions, including a colour reference for the simplices \textcolor{black}{in the middle row}.

Table (ref) shows the estimated coefficients (coefficient matrices) of the spatiotemporal autoregressive process along with \textcolor{black}{robust HAC sandwich standard errors, $t$ statistics, and $p$-values, and summary statistics of the model}. \textcolor{black}{We applied an isometric log-ratio transformation to all compositional observations, using a common orthonormal basis derived from a balanced sequential binary partition that splits the components into groups of approximately equal size. For the three-part composition, this yields two coordinates: the first contrasts developed land against condominium, while the second contrasts the joint component (developed land and condominium) against undeveloped land.} The spatial weight matrix has been chosen as a row-standardised contiguity matrix, i.e., all spatial effects should be interpreted as the impact of the average observation of all directly neighbouring regions, sharing common borders. \textcolor{black}{The temporal lags were selected based on AIC and BIC over a set of candidate schemes, which clearly favours the specification $\tau={1,6,12}$ (see Table (ref)), indicating that short-run persistence, semi-annual dynamics, and annual seasonality jointly contribute to the temporal dependence structure of the Berlin real-estate market.}

\textcolor{black}{Upon inspection of the estimated parameters in Table (ref), temporal dependence dominates contemporaneous spatial interaction. All diagonal elements of the lag-1 temporal matrix $\mathbf{\Pi}_{1}$ are large and statistically significant, indicating strong short-run persistence in both ilr coordinates. In addition, the diagonal elements of the seasonal lag-12 matrix $\mathbf{\Pi}_{12}$ are also significant, providing clear evidence of annual recurrence in the compositional dynamics. The off-diagonal elements of $\mathbf{\Pi}_{1}$ are statistically significant but substantially smaller in magnitude than the diagonal terms, indicating modest cross-coordinate temporal transmission. In contrast, off-diagonal elements of the higher-order temporal lags are not statistically significant, suggesting that seasonal feedback operates primarily within, rather than across, ilr coordinates. Spatial autoregressive effects are comparatively weak. Under robust HAC inference, the diagonal elements $\psi_{1,1}=0.0535$ and $\psi_{2,2}=0.0883$ is statistically significant at an $\alpha=0.05$ level, while the off-diagonal entries are not significant. This indicates that contemporaneous spatial dependence is present but limited, and concentrated in the second ilr balance, rather than uniformly across both coordinates.}

\textcolor{black}{Residual diagnostics indicate that the model adequately captures the main spatiotemporal dependence structure. The pooled Ljung–Box test does not reject global serial independence, and average residual autocorrelations is significantly different from zero for about two third of the locations, while Moran’s I statistics are close to zero, suggesting no systematic remaining spatial autocorrelation. The negative autocorrelation of the residuals indicates a slight overcompensation of the temporal autocorrelation. Although residuals exhibit moderate excess kurtosis and deviations from Gaussianity, stability conditions are satisfied and robust HAC inference accounts for remaining weak dependence and distributional departures. The reported goodness-of-fit measures indicate that the model captures a substantial share of the variation in the first ilr coordinate and a moderate share in the second, while the root mean squared errors on the simplex suggest a reasonably accurate fit across all three market segments, with the smallest prediction error observed for developed land and slightly larger deviations for condominiums.}

\textcolor{black}{Economic interpretation in original compositional space}

\textcolor{black}{Interpreting these findings on the original compositional scale requires caution due to the non-linear nature of the inverse ilr transformation. All reported effects on the simplex are evaluated locally at a representative reference composition, chosen as the overall mean composition across space and time. This choice reflects the typical market structure in the sample and ensures that reported effects correspond to empirically relevant states of the system. Since the covariate specification is seasonal (Fourier terms), the intercept component represents the baseline level of the seasonal trajectory in ilr space around which periodic fluctuations occur. When mapped back to the simplex and evaluated at the mean reference composition, the implied effects describe systematic reallocations of sales shares across segments rather than uniform shifts. As reported in Table (ref), the baseline intercept effect primarily operates through direct (own-location) adjustments, while spatial spillovers are economically negligible on average. Quantitatively, the baseline shift increases the share of condominium transactions by about $2.2$ percentage points and reduces the shares of developed land and undeveloped land by roughly $0.9$ and $1.3$ percentage points, respectively. This pattern suggests a gradual reallocation of market activity towards condominium transactions within the urban housing market, while the relative importance of land transactions declines. Such a shift is consistent with mature metropolitan markets, where development activity increasingly concentrates on condominium and apartment segments rather than on land turnover.}

\textcolor{black}{To illustrate contemporaneous spatial propagation, Table (ref) also reports the immediate simplex-scale effects of unit innovation shocks in both ilr coordinates at a representative district ($j=24$, postcode area 126xx, green compositions in Figure (ref)). By construction of the ilr basis, the first coordinate contrasts developed land with condominium, while the second coordinate contrasts the joint component (developed land and condominium) against undeveloped land. For a unit innovation in the first coordinate ($k=1$), the direct effect at the shocked district implies a substantial reallocation from developed land towards condominium transactions. At the reference composition, the condominium share increases by roughly $26.7$ percentage points, while developed land decreases by about $23.6$ percentage points and undeveloped land declines slightly (around $3.1$ percentage points). Spatial spillovers transmitted through the spatial multiplier remain small but reinforce the same directional adjustment, yielding total effects of approximately $+29.2$ percentage points for condominiums, $-25.0$ percentage points for developed land, and $-4.2$ percentage points for undeveloped land. A unit innovation in the second coordinate ($k=2$) generates a different compositional adjustment: the shares of \emph{developed land} and \emph{condominium} transactions increase moderately (about $2.9$ and $8.3$ percentage points, respectively), while the share of \emph{undeveloped land} decreases by roughly $11.2$ percentage points. Again, spatial spillovers remain small but slightly amplify the response, resulting in total effects of approximately $+3.2$, $+9.3$, and $-12.6$ percentage points for developed land, condominiums, and undeveloped land, respectively. Economically, this balance represents a shift of market activity away from undeveloped land and towards developed real estate segments. Taken together, these results indicate that local shocks primarily affect the composition of transactions within the originating district, while spatial propagation plays a secondary role. The dominant adjustment occurs through substitution between developed market segments and reductions in the relative share of undeveloped land transactions, which is consistent with the highly urbanised structure of the Berlin housing market.}

\textcolor{black}{Temporal dependence further shapes the evolution of compositional adjustments. Figure (ref) presents impulse response functions on the simplex for unit innovations in ilr coordinates $k=1$ (left column) and $k=2$ (right column) at the same spatial unit ($j=24$, eastern part, postcode area 126xx). All responses incorporate both spatial and temporal feedback. The impulse responses remain bounded and gradually decay over time, confirming that the estimated spatiotemporal system is dynamically stable. The immediate effects differ across the two balances implied by the ilr coordinates. A shock in the first coordinate ($k=1$) primarily reallocates activity from developed land towards condominium transactions, while the share of undeveloped land declines moderately. In contrast, a shock in the second coordinate ($k=2$) generates a stronger shift away from undeveloped land, accompanied by moderate increases in both developed land and condominium transactions. Following the initial response, the impulse responses decay over the first few periods, reflecting strong short-run persistence driven by the first-order temporal lag. At longer horizons, smaller secondary peaks appear around the 6- and 12-month horizons, indicating seasonal feedback effects captured by the higher-order lag matrices $\mathbf{\Pi}_{6}$ and $\mathbf{\Pi}_{12}$. These recurring but damped responses suggest that part of the initial disturbance re-enters the system through seasonal transaction cycles. Across all horizons, direct (own-unit) effects account for the bulk of the response, while average spillover effects across neighbouring districts remain comparatively small. Consequently, spatial interactions mainly propagate the shock locally without substantially altering the overall dynamic pattern. Taken together, the impulse responses indicate that shocks in the Berlin housing market primarily trigger within-district compositional reallocation between developed market segments and undeveloped land, while seasonal persistence and modest spatial feedback govern the subsequent adjustment path.}

Conclusion

Understanding how regional economic structures evolve over space and time requires models that accommodate both spatial and temporal dependence while respecting the compositional nature of the data. In this paper, we addressed this challenge through \textcolor{black}{a} real-world applications: the monthly evolution of housing market compositions across postcode areas in Berlin. \textcolor{black}{The example highlighted} key features of compositional areal data in practice---structural constraints, spatiotemporal dependence, and component-specific dynamics, which conventional models would fail to accommodate.

To address these challenges, we proposed a spatiotemporal multivariate autoregressive framework tailored for composition-valued panel data. The model combines a quasi maximum likelihood estimator with isometric log-ratio transformations, ensuring identifiability and asymptotic consistency under increasing space and time domains. In the Berlin housing market, we found evidence of moderate spatial spillovers and strong persistence in the share of undeveloped land. Cross-component interactions were weak in both cases, but statistically significant.

The flexibility and interpretability of the model make it a valuable tool for analysing compositionally constrained processes in regional economics and related fields. Future research could extend the framework in several directions. Methodologically, non-Gaussian extensions, such as compositional count models or models for zero-inflated data, could allow applications in domains like ecology, demography, or epidemiology. Allowing spatially or temporally varying autoregressive coefficients may help capture localised dynamics or structural heterogeneity. Furthermore, models that incorporate covariates in both the mean and dependency structure, or that allow for endogenous lagged effects across compositional components, would increase applicability. Finally, change-point extensions or regime-switching versions of the model could account for structural breaks or policy-induced shifts in regional dynamics.