Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,964 characters · 21 sections · 47 citation commands
Bayesian inference for dynamic spatial quantile models with interactive effects
{\bf Keywords}: Dynamic panel, endogeneity, factor models, heterogenous spatial effects, high dimensional data.
{\bf JEL classification}: C31, C33, E44
\baselineskip=23pt
Spatial panel data models explore spatial interactions across individuals or regions over time. These models have found wide-ranging applications in fields such as housing economics (e.g., BF2015), marketing (e.g., hunneman2022), and urban economics (e.g., Glaser2022). In this paper, we introduce a novel dynamic spatial panel quantile model with interactive effects. This model simultaneously accounts for spatial dependence, temporal dynamics, unobserved heterogeneity, and heterogeneous effects across quantiles—features that are often observed in economic and social data but rarely modeled together.
There is a rich literature on linear panel models with interactive effects (e.g., bai2009, baili2014; ando2017; hallin2007; moon2015; pesaran2006; stock2002; lu2016) and linear spatial panel models (e.g., aquaro2021; baltagi2011; bai2015; kelejian2004; lee2004; lin2010; li2017; lu2017; qu2015; reich2011; shi2017; yu2008). However, research on quantile regression (see KB1978) in spatial and high-dimensional panel settings remains limited, especially for models that incorporate both interactive effects and dynamic structures. Recent work by ando2020 extended quantile methods to panel data with interactive effects, capturing quantile co-movements while allowing for heterogeneous slope coefficients. However, their model does not accommodate dynamic and spatial dependencies in quantiles.
Spatial interactions and temporal dynamics are fundamental features of economics and the social sciences. Spatial dependence may arise from peer effects, spillover effects, or imitative behavior (see manski1993, anselin1988, among others), while temporal dependence often reflects persistence in individual or group behavior over time. To capture these complex features, we propose a dynamic spatial quantile panel model with interactive effects. This framework facilitates quantile-specific analysis of both spatial and temporal dependence, accommodates heterogeneous slope coefficients, and incorporates latent factor structures. Due to the model’s high-dimensional parameter space—driven by the dynamic structure and large number of cross-sectional units—we develop a new Bayesian Markov Chain Monte Carlo (MCMC) estimation method tailored to this setting.
{\bf Methodological challenges}: Several methodological challenges arise in the estimation process. First, because the quantile function is defined implicitly in terms of its coefficients, a conventional Gibbs sampler is not directly applicable. To overcome this, we introduce randomised quantile dynamics that create conditional conjugacy, allowing for efficient sampling. Second, the spatial autoregressive parameter is structural in nature and does not permit a standard Gibbs step. Inspired by methods in structural vector autoregression, we construct a bimodal mixture distribution to incorporate this parameter within our MCMC framework. Finally, to alleviate the computational burden from high-dimensional matrix operations, we apply a breadth-first search algorithm from graph theory to block-diagonalise the spatial matrix, thereby significantly reducing computational complexity. These methodological innovations enable our MCMC algorithm to efficiently estimate large-scale dynamic spatial quantile panel models with latent factor structures.
{\bf Theoretical Challenges:} To establish the validity of our Bayesian approach, we prove Bayesian consistency—a theoretical result not previously achieved for panel data models with interactive effects. This requires overcoming several substantial hurdles: the dynamic structure of the model, high-dimensional incidental parameters from latent factors and their loadings, the nonsmooth nature of the quantile objective function, the nonlinearity introduced by the spatial term, and the rotational indeterminacy of the factor structure. Although recent work (e.g., ando2020; chen2020) addresses some related issues, their results are not applicable to our setting. Consequently, we first derive new asymptotic properties within a frequentist framework (see Remark 7), which we then build upon to develop our Bayesian theory. Even with these foundations, our model demands additional theoretical innovation to rigorously establish Bayesian consistency.
Our contributions are summarised as follows. First, we introduce a dynamic spatial panel quantile model with interactive effects, designed to accommodate key features such as spatial dependence and temporal dynamics. This specification enhances the model’s flexibility in capturing complex data patterns. Second, we develop a novel Bayesian estimation procedure tailored to high-dimensional parameter spaces. Third, we establish several asymptotic properties, including Bayesian consistency. Finally, we demonstrate the utility of our framework through an empirical analysis of the Australian gasoline market.
The paper is organised as follows. Section 2 introduces a new dynamic spatial panel quantile model with interactive fixed effects. Section 3 outlines the proposed Bayesian MCMC estimation procedure. Section 4 establishes the Bayesian consistency of the method. In Section 5, we apply the proposed approach to the Australian gasoline market. Section 6 concludes the paper. To conserve space, all technical proofs are provided in the online supplementary material. The supplementary document also includes Monte Carlo simulation results, which demonstrate the effectiveness of the proposed estimation procedure.
{\bf Notations:} Let $\| A\|$ = $[\mathrm{tr}(A'A)]^{1/2}$ be the Frobenius norm of matrix $A$, where “tr” denotes the trace of a square matrix, and let $\|A\|_2$ be its spectrum norm (the largest singular value of $A$). In addition, for any $N\times N$ matrix $\|A\|_1$ is defined as $\|A\|_1=\max_{1\le j \le N}$ $\sum_{i=1}^N |a_{ij}|$ where $a_{ij}$ is the $(i, j)$-th element of $A$. Similarly, $\|A\|_\infty=\max_{1\le i \le N}\sum_{j=1}^N |a_{ij}|$. For sequences $a_n$ and $b_n$, the notation $a_n \lesssim b_n$ means $a_n = O(b_n)$, that is, there exists $C>0$ and for all $n$ large enough, $a_n \le C b_n$. We write $c_n = O_p(d_n)$ if $c_n/d_n$ is stochastically bounded, and $c_n = o_p(d_n)$ if $c_n/d_n$ converges to zero in probability.
Suppose that, for the $i$-th unit $(i=1,...,N)$ at time $t$ $(t=1,...,T)$, its response $y_{it}$ is observed together with a set of $p$ explanatory variables $\{x_{it,1},...,x_{it,p}\}$. We consider the $\tau$-th quantile function of $y_{it}$ by jointly modeling spatial effects, time effects and common shocks. To capture these effects simultaneously, we define the $\tau$-th quantile function of $y_{it}$ as
for $i=1, \dots ,N$ and $t=1, \dots ,T$. Here $w_{ij}$ $(i=1,2,\cdots,N; j=1,2,\cdots,N)$ are pre-specified spatial weights with $w_{ii}=0$, $\rho_{i,\tau}$ and $\delta_{i,\tau}$ are the heterogeneous spatial parameters capturing the strength of the spillover effects, the coefficients $\gamma_{i,\tau}$ are the heterogeneous temporal parameters, $\mbox{\boldmath $x$}_{it}=(1,x_{it,1},...,x_{it,p})'$ is $(p+1)$ -dimensional vector of explanatory variables; $B_{\tau}=(\mbox{\boldmath $b$}_{1,\tau}, \mbox{\boldmath $b$}_{2,\tau}, \dots, \mbox{\boldmath $b$}_{N,\tau})'$, $\mbox{\boldmath $b$}_{i,\tau}=(b_{i,0,\tau},b_{i,1,\tau},...,b_{i,p,\tau})' $ is a $(p+1)$-dimensional vector of regression coefficients; $\mbox{\boldmath $f$}_{t,\tau}=(f_{t1,\tau},...,f_{t r_\tau,\tau})'$ is $r_{\tau}$-dimensional unobservable common factors; $\mbox{\boldmath $\lambda$}_{i,\tau}=(\lambda_{i1,\tau},...,\lambda_{i r_\tau,\tau})'$ is $r_{\tau}$-dimensional vector of factor loadings; $X_t$ and $F_{t,\tau}$ are information on the explanatory variables and the common factors up to time $t$. The exact expression of $\mbox{\boldmath $x$}_{it}'\mbox{\boldmath $b$}_{i,\tau}$ is
where $e_{it}$ is the idiosyncratic error term and $G_{i,e_{it}}^{-1}(\tau)$ is the $\tau$-th quantile point of $e_{it}$ with $G_{i,e_{it}}(\cdot)$ being the cumulative distribution function of $e_{it}$. Thus, the $\tau$-th quantile of the idiosyncratic error $G_{i,e_{it}}^{-1}(\tau)$, which depends only on $i$ and $\tau$, is absorbed by the term $\mbox{\boldmath $x$}_{it}'\mbox{\boldmath $b$}_{i,\tau}$ since the first element of $\mbox{\boldmath $x$}_{it}$ is 1. We assume that $e_{it}$ is identically distributed over $t$ while its distribution may vary over $i$.
The preceding quantile function ((ref)) is associated with the following data-generating process, provided that the right-hand side of the equation is an increasing function of $u_{it}$,
where $y_{it, u_{it}}$ is $\tau={u}_{it}$-th quantile and ${u}_{it}$ are i.i.d. $U(0,1)$. The coefficient of the constant regressor absorbs the error term. As seen, we integrate spatial interactions and temporal dynamics into a model, which allows for both contemporaneous and dynamic spatial effects, thereby enabling the capture of temporal spillover effects and peer influences in the spatial domain.
{\remark koenker2006 considered autoregressive quantile model in the context of univariate time series. Their quantile function is expressed as a weighted sum of past observed values of the response variable. While their model is regarded as autoregressive in this sense, the quantile function itself is not autoregressive. In contrast, our quantile function in ((ref)) includes a weighted sum of past quantile functions, making the quantile itself autoregressive, or dynamic. }
Define the $N\times N$ matrix $S(\mbox{\boldmath $\rho$}_{\tau})\equiv (I-\mbox{\boldmath $\rho$}_{\tau} W)^{-1}$, where $\mbox{\boldmath $\rho$}_{\tau}=\mathrm{diag}(\rho_{1,\tau}, \dots, \rho_{N,\tau})$, and $W=[w_{ij}]$ is the $N\times N$ spatial weights matrix. Also, we define the $N\times N$ matrix $A(\mbox{\boldmath $\rho$}_{\tau},\mbox{\boldmath $\delta$}_{\tau},\mbox{\boldmath $\gamma$}_{\tau})\equiv (I-\mbox{\boldmath $\rho$}_{\tau} W)^{-1}(\mbox{\boldmath $\gamma$}_{\tau}+\mbox{\boldmath $\delta$}_{\tau} W)$ where $\mbox{\boldmath $\delta$}_{\tau}=\mathrm{diag}(\delta_{1,\tau}, \dots, \delta_{N,\tau})$, $\mbox{\boldmath $\gamma$}_{\tau}=\mathrm{diag}(\gamma_{1,\tau}, \dots, \gamma_{N,\tau})$. Stack the quantile functions over cross sections by defining \[\bm{Q}_t\Big(X_t, F_{t,\tau},B_\tau,\Lambda_\tau,\boldmath $\rho$_\tau,\boldmath $\gamma$_\tau,\boldmath $\delta$_\tau\Big)=
\] By recursive substitution, model ((ref)) can be rewritten as
where the $N\times N$ matrix $P_h(\mbox{\boldmath $\rho$}_{\tau},\mbox{\boldmath $\delta$}_{\tau},\mbox{\boldmath $\gamma$}_{\tau})$ is given as
We will impose restrictions on ((ref)) to ensure the sum to be well defined (referred to as stationarity). In addition, we assume $\bm{Q}_t(\cdot)=0$ when $t=0$. The effect of the initial condition is generally negligible if $T$ is large. Thus, an alternative expression of ((ref)) is
where $p_{ij,h}(\mbox{\boldmath $\rho$}_{\tau},\mbox{\boldmath $\delta$}_{\tau},\mbox{\boldmath $\gamma$}_{\tau})$ is the $(i,j)$th element of $P_h(\mbox{\boldmath $\rho$}_{\tau},\mbox{\boldmath $\delta$}_{\tau},\mbox{\boldmath $\gamma$}_{\tau})$ in ((ref)).
To eliminate the rotational indeterminacy of the common factor structure, we need to impose a restriction on $F_\tau=(\mbox{\boldmath $f$}_{1,\tau},...,\mbox{\boldmath $f$}_{T,\tau})'$ and $\Lambda_\tau=(\mbox{\boldmath $\lambda$}_{1,\tau},...,\mbox{\boldmath $\lambda$}_{N,\tau})'$. For example, bai2013 imposed the followings
where $I_{r_{\tau}}$ is an $r_{\tau}\times r_{\tau}$ identity matrix, and $D_{r_{\tau}}$ is a diagonal matrix whose diagonal elements are distinct and are arranged in a descending order. We refer to baing2013 for alternative restrictions on the common factor structure.
We have to estimate the unknown parameters $\mbox{\boldmath $\rho$}_{\tau}$, $\mbox{\boldmath $\gamma$}_\tau$, $\mbox{\boldmath $\delta$}_\tau$, $B_\tau$, $\Lambda_\tau$, and $F_\tau$ simultaneously. Let $\vartheta_\tau=\{\mbox{\boldmath $\rho$}_\tau, \mbox{\boldmath $\delta$}_\tau, \mbox{\boldmath $\gamma$}_\tau, B_\tau, \Lambda_\tau, F_\tau\}$. The frequentist estimator can be obtained as the minimiser of the following objective function
where $Q_{y_{it}}\left(\tau|X_t, \vartheta_\tau\right)\equiv Q_{y_{it}}(\tau|X_t, F_{t,\tau},B_\tau,\Lambda_\tau,\mbox{\boldmath $\rho$}_\tau,\mbox{\boldmath $\gamma$}_\tau,\mbox{\boldmath $\delta$}_\tau)$ is defined in ((ref)), $q_\tau (u)=u(\tau-I(u\le 0))$ is the quantile loss function, $Y \equiv \{y_{it}|i=1,...,N,t=1,...,T\}$ and $X \equiv \{\mbox{\boldmath $x$}_{it}|i=1,...,N,t=1,...,T\}$.
In this paper, we employ Bayesian MCMC sampling procedure to produce posterior samples from the posterior distribution
where $\pi(\vartheta_\tau)$ is a prior density of the parameter $\vartheta_\tau$. Details of $\pi(\vartheta_\tau)$ are given below.
The quantile structure in ((ref)) is a highly nonlinear function of parameters $\mbox{\boldmath $\rho$}_\tau, \mbox{\boldmath $\gamma$}_\tau$, $\mbox{\boldmath $\delta$}_\tau$ and $B_\tau$. Such nonlinearity costs the conditional conjugacy for inference, rendering the ideal Gibbs sampling infeasible. Simple methods such as the Metropolis-Hastings might be very inefficient due to the large number of parameters. For example, in our application, the dimension of the parameter space is greater than 44,000. In this paper, we propose to randomise the quantile equation to adapt to the Bayesian inference framework. In particular, we set
In the Bayesian framework, it is commonly assumed that $\varepsilon_{it,\tau}$ in (ref) has an Asymmetric Laplace distribution (ALD): $p(\varepsilon_{it,\tau}\mid\tau,\sigma ) =\frac{\tau(1-\tau)}{\sigma}\exp(-\frac{q_\tau(\varepsilon_{it,\tau})}{\sigma}) $. The Gibbs sampling method of kozumi2011gibbsuses a scaled mixture representation of the ALD. Under this representation, the conditional density of $y_{it}$, given the auxiliary variable $V_{it}$, is expressed as normal:
where $V_{it}\sim Exp(1)$ and $\sigma$ is the temperature parameter to facilitate Bayesian inference.
We denote $Q_{y_{it}}(\tau|X_t, \vartheta_\tau)$ as $Q_{it,\tau}$ in ((ref)), and $e^q_{it,\tau}\sim N(0, \sigma_{q,\tau}^2)$. Note that $e_{it,\tau}$ in ((ref)) and $e_{it,\tau}^q$ in ((ref)) are different. The introduction of $e^q_{it,\tau}$ is the randomisation that allows to treat the quantile values as a latent variable to augment the parameter space. Without $e^q_{it,\tau}$, the $Q_{it, \tau}$ is a deterministic function of $\vartheta_\tau$ in (ref) and hence cannot be treated as conditioning parameters. Consequently, the pseudo posterior kernel (ref) will be a complicated function of $\vartheta_\tau$, making posterior computation burdensome, particularly in high-dimensional settings. By introducing the randomised component $e^q_{it,\tau}$, the process $Q_{it,\tau}$ gains variability and can be considered part of the augmented parameter space. This randomisation enables efficient Gibbs sampling by allowing the conditional distributions to take tractable forms. In this way, equation (ref) plays a key role in making Gibbs sampling computationally feasible. See also our discussion in Remark 3.
We apply the same prior to the parameters for each quantile $\tau$. So ignore $\tau$ for notational simplicity. The parameter space includes:
{\remark
In addition to $\vartheta_\tau$ in the model ((ref)), the Bayesian parameter space also has auxiliary variables $\vartheta_\tau'=\{Q_\tau, V, \mbox{\boldmath $\phi$}, \sigma\}$, where $V$ and $Q_\tau$ are the collections of $V_{it}$ in ((ref)) and $Q_{it,\tau}$, respectively. Denote the corresponding posterior as
where $p(Y\mid X, \vartheta_\tau, V,\sigma)=\prod_{i=1}^N\prod_{t=1}^T p(y_{it} \mid X_t, \vartheta_\tau,V_{it} , \sigma) $ with $p(y_{it} \mid X_t, \vartheta_\tau,V_{it} , \sigma) $ given in ((ref)), and $\pi(\vartheta_\tau')$ is a prior density of $\vartheta_\tau'$. The concentrated posterior for comparison with the frequentist parameter set is defined as $p(\vartheta\mid X, Y) = \int p(\vartheta, \vartheta'\mid X, Y)d\vartheta'$. Although $p(\vartheta\mid X, Y)$ here and $\pi(\vartheta\mid X, Y)$ in ((ref)) are different, Proposition 1 in Section 3.2 below ensures that our MCMC posterior samples generated from $p(\vartheta\mid X, Y)$ can be regarded as those generated from $\pi(\vartheta\mid X, Y)$ through the importance sampling procedure.
}
We briefly describe the procedure in this section, with detailed techniques available in the Appendix. The computational challenge of our model arises from three main factors. First, the high dimensionality of the quantile values and their associated dynamics require a significant number of large matrix inversions. Second, due to the substantial heterogeneity in the model, generic methods, such as Metropolis-Hastings or particle filters, face high computational cost. Finally, the large dimension \( N \) and time period \( T \) result in a vast number of observations, further intensifying the computational burden. Consequently, our methods are designed to rely on the Gibbs sampler whenever possible.
The Markov Chain Monte Carlo (MCMC) algorithm for posterior inference proceeds as follows. Each step is performed conditionally on all other parameters. For clarity, subscripts are omitted, and we present the key elements below; full details are provided in the supplementary material.
This section obtains the “pseudo” posterior $p(\vartheta_\tau | X, Y)$. The exact posterior $\pi(\vartheta_\tau | X, Y)$ in ((ref)) can be built through importance sampling. Our simulation study demonstrates that this computationally efficient MCMC approach yields satisfactory results. For example, consider to compute the posterior expectation of $q(\vartheta_\tau)$ with $q(\cdot)$ being a known function of $\vartheta$: \[ E\left[q(\vartheta_\tau)\mid X, Y\right]= \int q(\vartheta_\tau) \pi(\vartheta_\tau\mid X, Y)d\vartheta_\tau. \] We can use the importance sampling: \[ E\left[q(\vartheta_\tau)\mid X, Y\right]= \frac{ \int q(\vartheta_\tau) \frac{ \pi_K(\vartheta_\tau\mid X, Y)}{p_K(\vartheta_\tau\mid X, Y)} p(\vartheta_\tau\mid X, Y)d\vartheta_\tau } { \int \frac{ \pi_K(\vartheta_\tau\mid X, Y)}{p_K(\vartheta_\tau\mid X, Y)} p(\vartheta_\tau\mid X, Y)d\vartheta_\tau, } \] where $\pi_K(\vartheta_\tau\mid X, Y)$ and $p_K(\vartheta_\tau\mid X, Y)$ are the kernels of $\pi(\vartheta_\tau\mid X, Y)$ and $p(\vartheta_\tau\mid X, Y)$, respectively. Because we have a set of posterior samples from $p(\vartheta_\tau\mid X, Y)$ and the kernels are easy to compute, the posterior expectation of $q(\vartheta_\tau)$ can be obtained without having $\pi(\vartheta_\tau\mid X, Y)$ directly. Appendix E.3 in supplementary document demonstrates Proposition 1 through a small Monte Carlo simulation.
Note that computing the posterior kernel \( p(\vartheta_\tau \mid X, Y) \) can be computationally expensive. However, after obtaining a sample using our method, it remains feasible, as the importance sampling process is parallelisable. A key aspect of successful importance sampling is having a good proposal distribution, and our method plays a critical role in ensuring that the posterior distribution is precisely targeted.
{\remark This section provides a sufficient approach, particularly when \( \sigma^2_q \) is tuned for robustness. An alternative way to address the computational burden in the importance sampling is to let $\sigma^2_q \rightarrow 0$ as the number of MCMC sampling increases. Then, we can have the posterior converging to the true posterior density \( \pi(\vartheta_\tau \mid X, Y) \). As a result, after a certain number of MCMC iterations, the generated sample from \( p(\vartheta_\tau \mid X, Y) \) is regarded as the sample from the true posterior density \( \pi(\vartheta_\tau \mid X, Y) \). }
In Bayesian theory, the number of factors can be viewed as a random variable. Its posterior can be inferred from by exploring the marginal likelihood of models with different value of $r$. However, because of the high computational cost, we adopt the idea of the sparse mixture approach if malsiner2016model by setting a large but finite number of factors. Each factor has a “switch” variables taking value of 0 and 1 indicating whether the corresponding factor is selected. This construction follows the literature of stochastic search variable selection dated back to mitchell1988bayesian.
In particular, we revise (ref) to have
where $\mbox{\boldmath $s$}_{\tau}$ is a $r_{\max} \times 1$ vector of $0$'s and $1$'s and the symbol $\circ$ means the Hadamard product. The maximum number of factors is $r_{\max}$. Each element $s_{j, \tau}$ in vector $\mbox{\boldmath $s$}_{\tau}$ controls for whether factor $f_{j,\tau}$ is selected.
Define the prior of $\mbox{\boldmath $s$}_{\tau}=(s_{1,\tau},...,s_{r_{\max},\tau})$ as \[P(s_{j,\tau}=1)=\pi,\quad P(s_{j,\tau}=0)=1-\pi\] for $j=1,...,r_{\max}$. The posterior distribution of the number of factors is the distribution of the sum of $\mbox{\boldmath $s$}_{\tau}$.
Because the only difference between (ref) and (ref) is $\mbox{\boldmath $s$}$, we only need to add one more step for $\mbox{\boldmath $s$}_{u_{it}}$ and revise the step of $F$ slightly in the original MCMC algorithm to infer the number of factors as follows.
{ \remark In practice, we tune $\sigma^2_q$ to control an $R^2$ level from (ref) or (ref). Alternatively, if we do not tune but want estimate $\sigma^2_q$ we can monitor the $R^2$ instead. For a high $R^2$ (in our application more than $98\%$), the parameters can be used as if they were drawn from the true posterior of the model. Alternatively, one can use the posterior sample from this model setting as a proposal distribution for an importance sampling scheme applied to the original model. Because the importance sampling is parallelisable, such second stage computation is much more affordable than any generic methods. }
To provide a theoretical justification for the Bayesian method, this section presents results on posterior consistency. A sequence of posterior distributions is considered consistent if, as the length of the time series and the number of cross-sectional units increase, the posterior converges to the degenerate measure at the true parameter value of the population density. Intuitively, posterior consistency ensures that the information from the quantile objective function outweighs the prior information. Before analysing the asymptotic behavior of the posterior distribution, we first need to examine the average consistency of the frequentist estimator. This is necessary to establish a set of proper conditions (conditions that have not yet been fully explored) that guarantee the convergence of the estimated model to the true population density.
Below, we denote the true spatial parameters, the lag coefficients and the true regression coefficient as $\rho_{i,0,\tau}$, $\delta_{i,0,\tau}$, $\gamma_{i,0,\tau}$ and $\mbox{\boldmath $b$}_{i,0,\tau}$, respectively. Similarly, we denote $F_{0,\tau}=(\mbox{\boldmath $f$}_{1,0,\tau},...,\mbox{\boldmath $f$}_{T,0,\tau})'$ and $\Lambda_{0,\tau}=(\mbox{\boldmath $\lambda$}_{1,0,\tau},...,\mbox{\boldmath $\lambda$}_{N,0,\tau})'$ as the true factors and loadings. A set of regularity conditions that are needed for theoretical analysis are given as follows.
Let $\mathcal{F}$ be a compact subset of $R^{r_\tau}$. The common factors $\mbox{\boldmath $f$}_{t,0,\tau}\in \mathcal{F}$ satisfy $T^{-1}\sum_{t=1}^T\mbox{\boldmath $f$}_{t,0,\tau}\mbox{\boldmath $f$}_{t,0,\tau}'=I_{r_\tau}$.
(B1) Let $\mathcal P, \mathcal D, \mathcal G$ be compact subsets of $\mathbb R$, and let $\mathcal{B}$ and $\mathcal{L}$ be compact subsets of $\mathbb R^{p+1}$ and $\mathbb R^{r_\tau}$, respectively. The spatial parameters $\rho_{i,0,\tau}$ and $\delta_{i,0,\tau}$, the lag coefficients $\gamma_{i,0,\tau}$, the regression coefficient $\mbox{\boldmath $b$}_{i,0,\tau}$, and the factor-loading $\mbox{\boldmath $\lambda$}_{i,0,\tau}$ satisfy that $\rho_{i,0,\tau}\in\mathcal P$, $\delta_{i,0,\tau}\in\mathcal D$, $\gamma_{i,0,\tau}\in\mathcal G$, $\mbox{\boldmath $b$}_{i,0,\tau}\in \mathcal{B}$ and $\mbox{\boldmath $\lambda$}_{i,0,\tau}\in \mathcal{L}$ for each $i$. \hangindent=30pt
(B2) The factor-loading matrix $\Lambda_{0,\tau}= (\mbox{\boldmath $\lambda$}_{1,0,\tau},\ldots,\mbox{\boldmath $\lambda$}_{N,0,\tau})'$ satisfies $N^{-1}\sum_{i=1}^N\mbox{\boldmath $\lambda$}_{i,0,\tau}\mbox{\boldmath $\lambda$}_{i,0,\tau}'\xrightarrow{p} \Sigma_{\Lambda_\tau}$, where $ \Sigma_{\Lambda_\tau}$ is an $r_\tau\times r_\tau$ positive definite diagonal matrix with diagonal elements distinct and arranged in the descending order. In addition, the eigenvalues of $\Sigma_{\Lambda_\tau}$ are distinct. \hangindent=30pt
(C1): The random variable \[\varepsilon_{it,\tau}=y_{it}-Q_{y_{it}}\Big(\tau| X_t, B_{\tau}, \mbox{\boldmath $f$}_{t,\tau}, \Lambda_{\tau},\mbox{\boldmath $\rho$}_{\tau}, \mbox{\boldmath $\delta$}_{\tau}, \mbox{\boldmath $\gamma$}_{\tau}\Big)\] satisfies $P(\varepsilon_{it,\tau}\le 0)=\tau$, and is independently distributed over $i$ and $t$, conditional on $X_t$, $B_{0,\tau}$, $F_{0,\tau}$, $\Lambda_{0,\tau}$, $\mbox{\boldmath $\rho$}_{0,\tau}$, $\mbox{\boldmath $\delta$}_{0,\tau}$ and $\mbox{\boldmath $\gamma$}_{0,\tau}$. \hangindent=30pt
(C2): The conditional density function of $\varepsilon_{it,\tau}$ given $\{X_t, B_{0,\tau}, F_{0,\tau},\Lambda_{0,\tau}, \mbox{\boldmath $\rho$}_{0,\tau}, \mbox{\boldmath $\delta$}_{0,\tau}, \mbox{\boldmath $\gamma$}_{0,\tau}\}$, denoted as $g_{it}(\varepsilon_{it,\tau})$, is continuous. In addition, for any compact set $\mathcal{C}$, there exists a positive constant $\underline g>0$ (depending on $\mathcal{C}$) such that $\inf_{c\in \mathcal{C}} g_{it}(c) \ge \underline g$ for all $i$ and $t$. \hangindent=30pt
\hangindent=30pt (D1): $W$ is an exogenous spatial weights matrix whose diagonal elements of $W$ are all zeros. In addition, $W$ is bounded by some constant $C$ for all $N$ under $\|\cdot\|_1$ and $\|\cdot\|_\infty$. \hangindent=30pt
(D2): The matrix $(I_N-\mbox{\boldmath $\rho$}_\tau W)^{-1}$ satisfies $\sup_{\bm{\rho}_\tau \in\mathcal P}\|I_N-\mbox{\boldmath $\rho$}_\tau W\|_2 <C$, and \[\sup_{\bm{\rho}_\tau \in\mathcal P}\Big(\Big\|(I_N-\mbox{\boldmath $\rho$}_\tau W)^{-1}\Big\|_1\vee\Big\|(I_N-\mbox{\boldmath $\rho$}_\tau W)^{-1}\Big\|_\infty\Big)<C\] where $C$ is some positive constant. \hangindent=30pt
(E1): For a positive constant $C$, explanatory variables satisfy $\sup_{it}\|\mbox{\boldmath $x$}_{it}\|\le C$ almost surely. \hangindent=30pt
(E2): Denote $\bm{Q}_{t,\tau}^0=\bm{Q}_t(X_t, F_{t,0,\tau},B_{0,\tau},\Lambda_{0,\tau},\mbox{\boldmath $\rho$}_{0,\tau},\mbox{\boldmath $\gamma$}_{0,\tau},\mbox{\boldmath $\delta$}_{0,\tau})$. Let $\mathcal X(B_{0,\tau})$ be an $N\times T$ matrix with its $(i,t)$-th entry $\mbox{\boldmath $x$}_{it}'\mbox{\boldmath $b$}_{0,\tau}$. Define $u_{it,0,\tau}$ to be the $(i,t)$th element of $\mathcal U_{0,\tau}$ with \hangindent=30pt \[\mathcal U_{0,\tau}=W(I_N-\mbox{\boldmath $\rho$}_{0,\tau}W)^{-1} \Big[\mathcal X(B_{0,\tau})+\Lambda_{0,\tau}F_{0,\tau}'+(\mbox{\boldmath $\gamma$}_{0,\tau}+\mbox{\boldmath $\delta$}_{0,\tau}W)\mathcal Q_{-1,\tau}^0\Big].\] where $\mathcal Q_{-1,\tau}^0=[\bm{Q}_{0,\tau}^0, \bm{Q}_{1,\tau}^0, \dots, \bm{Q}_{T-1,\tau}^0]$. Let $\mbox{\boldmath $z$}_{it, \tau}=(u_{it,0,\tau}, Q_{i,t-1,\tau}^0, \sum_{j=1}^N w_{ij}Q_{j,t-1,\tau}^0, \mbox{\boldmath $x$}_{it}')'$ with $Q_{i,t-1}^0$ being the $i$-th element of $\bm{Q}_{t-1}^0$, and $\mbox{\boldmath $Z$}_{i,\tau}=(\mbox{\boldmath $z$}_{i1,\tau}, \mbox{\boldmath $z$}_{i2,\tau}, \dots, \mbox{\boldmath $z$}_{iT,\tau})'$. Further define $A_{i,\tau}=\frac{1}{T}\mbox{\boldmath $Z$}_{i,\tau}'M_{F_\tau} \mbox{\boldmath $Z$}_{i,\tau}$, $B_{i,\tau}=(\mbox{\boldmath $\lambda$}_{i,0,\tau}\mbox{\boldmath $\lambda$}_{i,0,\tau}') \otimes I_T$, $C_{i,\tau}=\frac{1}{\sqrt{T}}[\mbox{\boldmath $\lambda$}_{i,0,\tau} \otimes (M_{F_\tau}\mbox{\boldmath $Z$}_{i,\tau})]'$ with $M_{F_\tau}=I-F_\tau(F_\tau'F_\tau)^{-1}F_\tau'$. Let ${\mathcal{F}}_\tau$ be the collection of $F_\tau$ such that ${\mathcal{F}}_\tau=\{F_\tau:F_\tau'F_\tau/T=I_{r_\tau}\}$. We assume that with probability approaching one, \hangindent=30pt
where $\lambda_{\min}(A)$ denotes the smallest eigenvalue of matrix $A$, and $E_{i,\tau}(F_\tau)=B_{i,\tau}-C_{i,\tau}'A_{i,\tau}^{-1}C_{i,\tau}$. \hangindent=30pt
(E3): For each $i$, we assume that there exists a constant $c>0$ such that for each $i$, with probability approaching one, \hangindent=30pt \[\liminf_{T\to\infty}\lambda_{\min}\Big(\frac1T\mbox{\boldmath $Z$}_{i,\tau}'M_{F_{0,\tau}}\mbox{\boldmath $Z$}_{i,\tau}\Big)\ge c.\]
The data generating process from ((ref)) is assumed to be stationary. To ensure the stationarity, it is assumed that $\bar c <1$ where \[\bar c =\sup_{\bm{\rho}_\tau\in\mathcal P, \bm{\delta}_\tau\in\mathcal D, \bm{\gamma}_\tau\in\mathcal G}\Big|\lambda_{\max}\Big(A(\mbox{\boldmath $\rho$}_\tau, \mbox{\boldmath $\delta$}_\tau, \mbox{\boldmath $\gamma$}_\tau)\Big)\Big|<1.\] where $\lambda_{\max}(A)$ denotes the eigenvalue of $A$ with the largest modulus.
{ \remark Assumptions A and B on the factors and factor loadings are standard in factor models. Similar to ando2020, the factors and factor loadings are treated as parameters. Assumptions C and D on the idiosyncratic errors and the spatial weighting matrix are standard assumptions in the literature. Assumption E is necessary for deriving the consistency of the frequentist estimator (See ando2020 for similar assumptions). Assumption F is a stationary condition similar to yu2008. Similar to the investigation in yu2008, a sufficient condition for Assumption F is $\sup\|A(\mbox{\boldmath $\rho$}_{\tau},\mbox{\boldmath $\delta$}_{\tau},\mbox{\boldmath $\gamma$}_{\tau})\|_2 < 1$. The stationary condition is verified accordingly in both simulation and empirical analysis. }
Now, we investigate the consistency of the frequentist estimator, defined as the minimiser of ((ref)) subject to a normalisation condition. Recall $\vartheta_\tau=\{\mbox{\boldmath $\rho$}_\tau, \mbox{\boldmath $\delta$}_\tau, \mbox{\boldmath $\gamma$}_\tau, B_\tau, \Lambda_\tau, F_\tau\}$, and $\hat{\vartheta}_\tau=\{\hat{\mbox{\boldmath $\rho$}}_\tau, \hat{\mbox{\boldmath $\delta$}}_\tau, \hat{\mbox{\boldmath $\gamma$}}_\tau, \hat{B}_\tau, \hat{\Lambda}_\tau, \hat{F}_\tau\}$ denotes the frequentist estimator. A set of assumption A--F leads to the following result.
{\rm \remark The last claim $\frac1{NT}\|\hat\Lambda_\tau\hat F_\tau'-\Lambda_{\tau,0}F_{\tau,0}'\|^2=O_p(\delta_{NT}^2)$ further implies \[\frac1N\sum_{i=1}^N \|\hat{\mbox{\boldmath $\lambda$}}_{i,\tau}-\mbox{\boldmath $\lambda$}_{i,\tau,0}\|^2=O_p(\delta_{NT}^2), \quad \frac1T\sum_{t=1}^T \|\hat{\mbox{\boldmath $f$}}_{t,\tau}-\mbox{\boldmath $f$}_{t,\tau,0}\|^2=O_p(\delta_{NT}^2). \] }
{\rm \remark The structure ((ref)) indicates that the quintile function at time $t$ is expressed as the sum of the previous quintile functions up to time $t-1$. This created a difficult technical challenge to establish the claims in Theorem (ref). More specifically, as an intermediate result, the technical proof of our Lemma 2 (in the supplementary document) establishes $\frac1N\sum_{i=1}^N \|\hat{\phi}_{i,\tau}-\phi_{i,\tau,0}\|^2=o_p(1)$, where $\mbox{\boldmath $\phi$}_{i,\tau,0}=(\rho_{i,\tau,0}, \gamma_{i,\tau,0}, \delta_{i,\tau,0}, \mbox{\boldmath $b$}_{i,\tau,0}')'$, $\hat\mbox{\boldmath $\phi$}_{i,\tau}=(\hat\rho_{i,\tau}, \hat\gamma_{i,\tau}, \hat\delta_{i,\tau}, \hat\mbox{\boldmath $b$}_{i,\tau}')'$. This result was not available from previous studies (ando2020, chen2020) due to the static nature of their model. This illustrates one of the key challenges when we establish Theorem (ref).
}
The next theorem also plays an important role when we investigate the consistency of our Bayesian MCMC procedure. Theorem (ref) implies that it is ideal to set the number of common factors equal to or greater than the true number of common factors when one's focus is the consistent estimation of parameters $\mbox{\boldmath $\rho$}_\tau$, $\mbox{\boldmath $\delta$}_\tau$, $\mbox{\boldmath $\gamma$}_\tau$ and $B_\tau$. To obtain the claim, we need additional assumption.
Let $F_\tau(k)$ be the common factor matrix with $k>r$ and $r$ being the true number of common factors, $\mathcal Z_\tau(\mbox{\boldmath $\phi$}_\tau)$ be the $N\times T$ matrix with its $(i,t)$th entry equal to $\mbox{\boldmath $z$}_{it,\tau}'\mbox{\boldmath $\phi$}_{i,\tau}$, where $\mbox{\boldmath $\phi$}_\tau=(\mbox{\boldmath $\phi$}_{1,\tau}, \mbox{\boldmath $\phi$}_{2,\tau}, \dots, \mbox{\boldmath $\phi$}_{N,\tau})'$ and $\mbox{\boldmath $\phi$}_{i,\tau}=(\rho_{i,\tau}, \gamma_{i,\tau}, \delta_{i, \tau}, \mbox{\boldmath $b$}_{i,\tau}')'$. Here $\mbox{\boldmath $z$}_{it,\tau}$ is defined in Assumption E.2. There exists a positive constant $\breve c>0$ such that with probability approaching one, \[ \inf_{F_{\tau}(k),F_{\tau}(k)'F_{\tau}(k)/T=I_k} \frac1{NT}\|M_{\Lambda_{0,\tau}}\mathcal Z_{\tau}(\mbox{\boldmath $\phi$}_\tau)M_{F_{\tau}(k)}\|^2\ge \breve c \frac1N\sum_{i=1}^N\|\mbox{\boldmath $\phi$}_{i,\tau}\|^2,\] where $M_{\Lambda_{0,\tau}}=I-\Lambda_{0,\tau}(\Lambda_{0,\tau}'\Lambda_{0,\tau})^{-1}\Lambda_{0,\tau}'$.
Now, our concern is the sequence of posterior distributions $\pi(\vartheta_\tau|Y,X)$ constructed by the size of $T\times N$ panel data, generated from the true density $f(Y|X,\vartheta_{\tau,0})$. In this paper, we show that the constructed posterior forms a Hellinger-consistent sequence. In regards to the posterior consistency based on Hellinger distance, we refer to barron1999, Ghosal, walker2001.
Recall the pseudo-likelihood based density function: \[ f(Y|X,\vartheta_\tau)\propto \exp \left[ -\frac{1}{NT} \sum_{i=1}^N\sum_{t=1}^T q_{it,\tau}(\vartheta_\tau) \right], \] where $q_{it,\tau}(\vartheta_\tau)\equiv q_\tau(y_{it}-Q_{y_{it}}(\tau|X_t, F_{t,\tau},B_\tau,\Lambda_\tau,\mbox{\boldmath $\rho$}_\tau,\mbox{\boldmath $\gamma$}_\tau,\mbox{\boldmath $\delta$}_\tau))$. We establish Bayesian consistency under the pseudo-likelihood $f(Y|X,\vartheta_\tau)$, in the sense that, for any $\mu>0$,
where $\vartheta_{\tau,0}$ is true value of $\vartheta_{\tau}$, $\pi_{N,T,\tau}(\cdot)$ is defined as
where $A\subset \Theta$ is the subset of parameter space $\Theta$, and for two density functions $h(y)$ and $g(y)$, the Hellinger distance is defined as $ H(h,g)=\left\{\int (g^{1/2}(y)-h^{1/2}(y)) d(y)\right\}^2 $. To obtain the result in ((ref)), we need an additional condition.
Let $K_\varepsilon (\vartheta_{\tau,0})$ be a Kullback--Leibler neighborhood of $\vartheta_{\tau,0}$ such that $\vartheta_\tau$ satisfies
Then, the prior density $\pi(\vartheta_\tau)$ assigns positive mass on all Kullback--Leibler neighborhoods of the pseudo-likelihood based density under the true value $f(Y|X,\vartheta_{\tau,0})$, $ \pi (K_v (\vartheta_{\tau,0}))>0$.
When Bayesian consistency is considered important, the above theorem offers guidance for designing an appropriate prior distribution. The prior distribution discussed in Section 3 is specifically constructed to satisfy Assumption H. As a result, we expect that the posterior mean under our prior will converge to the true parameter values as both $N$ and $T$ tend to infinity. This expectation is confirmed through our simulation study.
{\remark A common question is about the asymptotic properties of the posterior distribution regarding the number of common factors in our MCMC procedure. When the number of common factors is strictly smaller than the true number, there exists a positive constant such that the expected quantile loss is larger than that under the true number of common factors. As a result, our MCMC procedure asymptotically eliminates posterior samples with $r$ smaller than the true number of common factors as both $N$ and $T$ tends to infinity. }
We apply our method to the fuel prices reported by retailers in Queensland, a state located in the northeast of Australia. In Queensland, an aggregation system for fuel price reporting has been established under Section 4 of the Fair Trading (Fuel Price Reporting) Regulation 2018.\footnote{See https://www.epw.qld.gov.au/about/initiatives/fuel-price-reporting} As a result, all fuel retailers in Queensland (including all fuel stations) are required to report their fuel prices as part of the Queensland fuel price reporting scheme, which helps motorists find the cheapest fuel prices. This requirement has been in effect since 3 December 2018. The data is publicly available at \url{https://www.data.qld.gov.au/dataset/fuel-price-reporting}.
Figure (ref) (a) shows the locations of these stations across Queensland, as well as their brands.\footnote{The brand small aggregates the smaller brands with fewer than 5 stations.} In total, there are 1011 registered stations in Queensland. Figure (ref) (b) displays the locations of fuel stations in Brisbane, the capital city of Queensland, along with their respective brands. It is evident that many stations are situated close to one another, particularly in the municipal area.
We analyse unleaded gasoline, as in pinkse2002spatial, because it has the largest market share and is nearly a homogeneous product. Our data spans from February 1, 2019, to September 30, 2021, with a total of \( T = 973 \) days and no missing values. This sample includes five lockdown periods in Queensland during the Covid-19 pandemic.
Figure (ref) shows the average fuel prices for each brand over time, revealing clear seasonality. The significant drop in fuel prices in 2020 corresponds to the longest lockdown period, from March 26 to the end of April. Additionally, when aggregating over time, Figure (ref) highlights the price heterogeneity across brands.
In the following analysis, we exclude small brands with fewer than 5 stations and stations located on islands. After this cleaning process, we are left with \( N = 946 \) stations. Throughout the sample period, no new stations were built, nor were any stations decommissioned. While some stations changed brands, we account for the brand effect in the estimation.
For the term $\mbox{\boldmath $x$}_{it}'\mbox{\boldmath $b$}_{i,\tau}$ in ((ref)), we apply the following empirical specification:
In the above equation, the subscript \(i\) refers to fuel station \(i\), and \(t\) represents time (day). The independent variables in this model include several time effects: \({\rm Year}_t\), \({\rm Month}_t\), and \({\rm DayofWeek}_t\). The year effect captures any overall trend, while the month and day-of-the-week effects capture explicit seasonal variations. Each station \(i\) has its own parameter for these time variables. For instance, the effect of the year 2021 differs for station 1 compared to station 2.
The brand \({\rm Brand}_{it}\) dummy variables capture the brand effect. Note that some stations changed brands during the sample period, which is why a double subscript is used for this variable. We include all brand dummy variables and exclude the intercept for identification purposes, so the station-specific effect is naturally incorporated.
The variable \({\rm Covid}_t\) is a dummy variable that takes the value 1 if time \(t\) falls within a lockdown period and 0 otherwise. Details of the lockdown periods are provided in the appendix. Briefly, there was only one long lockdown period in 2020, from March 26 to the end of April, with all other periods being less than one week in duration.
For the weight matrix \( W \), we consider driving distance rather than geographic distance, taking into account traffic conditions and speed limits in different areas. To compute the average driving time between two stations, we use the Open Source Routing Machine (OSRM). In this application, if two stations are within a 5-minute driving distance of each other, they are classified as neighbors. We normalise each row of the matrix \( W \) such that the sum of the elements in each row equals 1. Specifically, if station 1 has three neighbors (stations 2, 7, and 9), then \( W_{1,2} = W_{1,7} = W_{1,9} = \frac{1}{3} \), and all other \( W_{1,j} = 0 \).
The prior is set to be informative but covers a broad range of the parameter space, as in the simulation. The detailed settings can be found in the appendix.
Due to the large number of parameters in the model, such as the number of \( \rho \), \( \delta \), \( \gamma \), and \( b \), \( \Lambda \), \( F \), which totals \( 3N + N(p+1) + r_f(T + N) = 44,597 \) in our application, we do not store all simulated values during the MCMC process. Instead, we focus on the posterior means, which allows us to accumulate values with minimal memory usage. This approach makes it straightforward to evaluate uncertainties. For example, if posterior variance is required, we can save the sum of the squared values, then use the sample mean of the squared values and the sample mean to compute the sample variance. Any moment-based posterior statistics can be derived in this way.
Figure (ref) presents the histogram of the posterior means of \( \rho_{i,\tau} \) for all \( i \) at the quantiles \( \tau = (0.01, 0.05, 0.5, 0.95, 0.99) \). Without imposing any restrictions, the distribution of \( \rho_{t,\tau} \) reveals two key characteristics. First, the values are positive, indicating that positive spillover effects exist between fuel prices at nearby stations. Second, the heterogeneity patterns are consistent across different quantiles. Some stations exhibit greater sensitivity to their neighbors (larger \( \rho \) values), while others show minimal sensitivity, with \( \rho \) values close to zero. Figure (ref) underscores the need for a heterogeneous coefficient model.
Figure (ref) shows the distribution of the posterior means of the Covid lockdown coefficients at different quantiles. It is clear that the distributions of these coefficients differ significantly across quantiles. For instance, the histogram for \( \tau = 0.05 \) is right-skewed, while the histogram for \( \tau = 0.99 \) is left-skewed. The mode is around zero. This is not surprising, as the factors are intended to capture any systematic price changes. Figure (ref) demonstrates that the lockdowns have caused changes in prices, not only in terms of dispersion but also in how these price responses vary across different quantiles.
Figure (ref) displays the average brand premiums across different quantiles. It is important to note that each station has its own distinct brand premium, which can be interpreted as the average station effect within the same brand. Similar patterns emerge across the brand premiums. For instance, Pacific Petroleum consistently has the lowest values, while United Petroleum consistently has the highest values across all quantiles.
Similar to the spatial coefficient $\rho_{i,\tau}$, the lag coefficient $\gamma_{i,\tau}$ and the lag-spatial coefficient $\delta_{i,\tau}$ demonstrate a strong pattern of heterogeneity, while exhibiting similar patterns across quantiles, respectively. Due to page limitations, these are not shown here.
The systemic factors and their effects on each time series are plotted in Figure (ref). Each subplot represents $N=946$ time series, depicting $\mbox{\boldmath $\lambda$}_{i,\tau}'\mbox{\boldmath $f$}_{t,\tau}$. It is evident that the factor structure captures the seasonality present in the data. Since the simulation study does not indicate the correct number of factors, we refrain from analysing individual factors in this application.
We compute the average variation from the posterior mean of the contemporaneous effect. Specifically, we consider the term $\rho_{i,\tau} \sum_{i \neq j, j = 1}^N w_{ij} Q_{y_{jt}} \left( \tau | X_t, F_{t,\tau}, B_\tau, \Lambda_\tau, \mbox{\boldmath $\rho$}_\tau, \mbox{\boldmath $\gamma$}_\tau, \mbox{\boldmath $\delta$}_\tau \right)$ in (ref). The ratio of the average variation of this term to the average variation of the posterior values of the quantiles is approximately $30\%$ for all quantiles in $(0.01, 0.05, 0.5, 0.95, 0.99)$. This suggests that the contemporaneous spatial effect plays a significant role in explaining the quantiles.
In this paper, we proposed a novel dynamic spatial panel quantile model with interactive effects. The model effectively captures several complex features simultaneously, including spatial spillover effects, heterogeneous regression coefficients, and unobserved heterogeneity that vary across quantiles. To estimate the model, we developed a Bayesian MCMC procedure capable of handling all these aspects, and we established Bayesian consistency to support its theoretical validity. The practical utility of the method was demonstrated through an application to gasoline price data in Australia.