Bayesian P-Splines Quantile Regression of Partially Linear Varying Coefﬁcient Spatial Autoregressive Models

: This paper deals with spatial data that can be modelled by partially linear varying coefﬁcient spatial autoregressive models with Bayesian P-splines quantile regression. We evaluate the linear and nonlinear effects of covariates on the response and use quantile regression to present comprehensive information at different quantiles. We not only propose an empirical Bayesian approach of quantile regression using the asymmetric Laplace error distribution and employ P-splines to approximate nonparametric components but also develop an efﬁcient Markov chain Monte Carlo technique to explore the joint posterior distributions of unknown parameters. Monte Carlo simulations show that our estimators not only have robustness for different spatial weight matrices but also perform better compared with quantile regression and instrumental variable quantile regression estimators in ﬁnite samples at different quantiles. Finally, a set of Sydney real estate data applications is analysed to illustrate the performance of the proposed method.


Introduction
Spatial econometric models can deal with the spatial correlation and heterogeneity of variables in cross-sectional data and panel data, which expands the application scope of traditional econometric models.Among many spatial econometric models, there is a large amount of literature focusing on the spatial autoregressive (SAR) model [1].For instance, Lee [2] studied the asymptotic properties of the quasi-maximum likelihood estimator of the spatial autoregressive model.Lee [3] proposed the generalized method of moments method and classical two-stage least-squares method to estimate mixed regressive and spatial autoregressive models.Kakamu and Wago [4] compared the Bayesian method and maximum likelihood method to study small-sample properties of the panel spatial autoregressive model.Xu and Lee [5] considered the instrumental variable and maximum likelihood estimation for a spatial autoregressive model with a nonlinear transformation of the dependent variable among others.
However, these studies mainly focused on parametric models and could not well explain complex economic phenomena.In fact, the relationship among many economic variables is nonlinear [6][7][8][9].In order to improve the flexibility and applicability of spatial econometric models, the research of semiparametric spatial econometric models has been gradually increasing.For example, Sun et al. [10] proposed semiparametric varying coefficient spatial autoregressive models.Chen et al. [11] developed a two-stage Bayesian estimation method for semiparametric spatial autoregressive models.Dai et al. [12] applied the quantile regression approach to partially linear varying coefficient spatial autoregressive models.Cai and Xu [13] constructed varying coefficient quantile regression models for time series data.
Although nonparametric spatial autoregressive models can improve the performance of the parametric spatial autoregressive model, it is unavoidable to suffer from the "curse of dimensionality" problem [14].To solve the problem, a few dimension-reduction approaches have been developed in some literature, including additive model [15], singleindex model [16] and varying coefficient model [17], to name a few.Among the different semiparametric models, the partially linear varying coefficient model is perhaps the most widely used.It not only contains the advantages of a linear model but also retains the robustness of nonparametric regression, which can overcome the "curse of dimensionality" problem well.Many scholars have enriched the estimation methods of varying coefficient models both theoretically and empirically, such as the splines method [18][19][20], the kernel method [21], local polynomial estimation [22,23], basis function approximation [24,25] and so on.
According to the regression object of the model, the semiparametric spatial econometric model usually includes mean regression model and quantile regression model.Most of the spatial econometric models involved in the existing literature belong to the former [26][27][28][29].The mean regression model can only reflect the location information of the conditional distribution of dependent variable and cannot describe its scale and shape.On the contrary, the quantile regression model [30] can find the location, scale and shape of the conditional distribution of dependent variable and, in particular, capture the tail characteristics of the distribution.While linear regression needs to assume that the random error term obeys the normal distribution or generalized Gauss-Laplace distribution [31], the quantile regression approach has strong robustness without making any distribution assumption for the random error terms.Many early studies have been summarized in Koenker and Bassett [30], Koenker and Machado [32], Zerom [33], Chernozhukov and Hansen [34], Su and Yang [35], which considered quantile regression from both a frequentist and a Bayesian point of view.Form the former perspective, the estimation method relies on the minimum asymmetric absolute loss function [36].Concerning the Bayesian approach, Yu and Moyeed [37] introduced the asymmetric Laplace distribution [38] as a working likelihood to perform the inference.Bayesian quantile regression has since been implemented in a wide range of applications, including models for longitudinal studies [39], Lasso regression [40] and spatial analysis [41], among others.As a result, we develop Bayesian quantile regression for a partially linear varying coefficient spatial autoregressive model.
The estimation methods of semiparametric spatial autoregressive models include quasi-maximum likelihood estimation method [42], instrumental variable method [43], generalized method of moments [44] and Bayesian estimation methods [45].Compared with the frequentist solutions, Bayesian estimation methods can infer the posterior distributions of parameters by utilizing prior information and allow for parameter uncertainty to be considered.To the best of our knowledge, Bayesian inference for quantile regression has been applied by few authors, such as [37,46,47].In addition, P-splines [48] have been a popular way for estimating nonlinearities in semiparametric models.This has attractive properties, including that each piecewise polynomial only forms a local basis with unit integrals and overlaps with a limited number of other polynomials [49].As they are composed of piecewise polynomials, the upper range of basis function is limited, and the differentials of splines are readily available [50].These characteristics ensure the P-splines are available both numerically and analytically, and P-splines can approximate nonparametric components in the semiparametric spatial autoregressive model.Hence, we can apply a Bayesian P-splines method to approximate nonparametric functions using penalized splines with fixed number and location of knots inferred through Markov chain Monte Carlo (MCMC).Gibbs sampler is a common sampling method employing the MCMC technique, which can generate simple random samples from various distributions, including uniform distribution [51].
In this paper, we developed the partially linear varying coefficient spatial autoregressive (PLVCSAR) models with Bayesian quantile regression using the asymmetric Laplace error distribution for spatial data.It allows for different degrees of spatial dependence at different quantile points of the response distribution.The PLVCSAR model is a good balance between flexibility and parsimony, which can simultaneously capture linearity, non-linearity and the spatial correlation relationship of exogenous variables in a response variable.We employ a Bayesian P-splines method to estimate the unknown parameters and approximate the varying coefficient functions, and we also design a Gibbs sampler to explore the joint posterior distributions using the MCMC technique.It may update the iterations to draw parameters from the full conditional posterior distributions of the unknown quantities through appropriate selection of prior information, which makes Bayesian inference efficient and useful even in complicated situations.The proposed model combines a spatial autoregressive model with a semiparametric framework in an adaptive way, and estimators of this article may be a great breakthrough and improvement in the field of related research.
The rest of the paper is organized as follows.In Section 2, we introduce Bayesian quantile regression of PLVCSAR models for spatially dependent responses and discuss the identifiability conditions, and then we obtain the likelihood function by approximating the varying coefficient functions with the P-splines method.In Section 3, we give the prior distributions and infer the full conditional posteriors of latent variables and unknown parameters.We also describe the detailed Gibbs sampler procedure.Simulation studies for assessing the finite sample performance of the proposed method are reported, and an empirical example is illustrated in Section 4. We summarize the article in Section 5.

Model
Given the following partially linear varying coefficient spatial autoregressive model with quantile regression where y i is the dependent variable, ) and u i are the associated explanatory variables, w ij is the (i, j)th element of an exogenously given spatial weight matrix with known constants, α τ (•) = (α τ1 (•), • • • , α τq (•)) consists of a q-dimensional vector of unknown smooth functions, ρ τ denotes the τth quantile spatial regression parameters and is restricted to the condition parameters, u i is the smoothing variable, and ε τ,i is the random error with the τth quantile on (x i , z i , u i ), which equals zero for τ ∈ (0, 1).

Likelihood
We assume that ε τ,i are mutually independent and identically distributed random variables from an asymmetric Laplace distribution with the density where µ is the location parameter, σ 0 is the scale parameter, and λ τ (ε) = ε(τ − I(ε < 0)) is called the check function.Then, the conditional distribution of y is in the form of Quantile regression is typically based on the check loss function to solve a minimization problem.With model (1), the specific problem is estimating ρ τ , β τ and α τ (•) by minimizing the following objective function and a likelihood-based interpretation.By introducing the location-scale mixture representation of the asymmetric Laplace distribution [52], model (1) can be equivalently written as where e i ∼ exp(1/σ 0 ) with mean σ 0 and In the following expressions, we omit τ for ease of notation.Considering the advantages of the Bayesian P-splines method, we intend to approximate varying coefficient function α j (•) in (1) with P-splines.For j = 1, • • • , q, the unknown function α j (•) is a polynomial spline of degree t j with k j order interior knots where ) is a K j × 1 vector of spline basis, which is determined by the knots, γ j = (γ j1 , • • • , γ jK j ) is a K j × 1 vector of spline coefficients, and boundary knots are It follows from (5) that model (4) can be written as where We view e i as latent variables for i = 1, • • • , n and define e = (e 1 , • • • , e n ) .The matrix form of the model ( 7) is where tor with all elements being 1, W = (w ij ) is an n × n specified constant spatial weight matrix, and , where D j (z, u) is an n × K j matrix.

Bayesian Estimation
In this section, we construct a Bayesian P-splines method with a Gibbs sampler to analyse the proposed model.First of all, we specify the prior distributions of the unknown parameters, and then we infer the full conditional posteriors and describe the detailed Gibbs sampler procedure.

Priors
According to the Bayesian P-splines method, we need to provide appropriate prior distributions for all the unknown parameters, including spatial autocorrelation coefficient ρ, regression coefficient vector β, spline coefficient vector γ and the location and scale parameters µ and σ 0 .
Firstly, we choose a hierarchical prior for β, which consists of conjugate normal prior and an inverse-gamma prior where r τ 0 and s 2 τ 0 are pre-specified hyper-parameters.Secondly, we choose the random walk prior for γ j where d is the order of the random walk, M γ j is the penalty matrix that equals for d-order random walk prior, P l is a (K j − l − 1) × (K j − l) matrix with the form For j = 0, 1, • • • , q, the prior of hyper-parameters τ j is given by where r τ j and s 2 τ j are pre-specified hyper-parameters.In addition, we give no prior information for the location parameter µ and a conjugate normal inverse-gamma prior for the scale parameter σ 0 π(µ) ∝ 1, where r 0 and s 2 0 are also pre-specified hyper-parameters.We select r 0 = s 2 0 = 1 to obtain a Cauchy distribution of σ 0 and use r τ j = 1 and s 2 τ j = 0.005 to obtain a highly dispersed inverse gamma prior for each hyper-parameter of τ j for j = 0, 1, • • • , q. Lastly, the spatial autocorrelation coefficient ρ is set a uniform prior for ρ ∼ U(λ −1 min , λ −1 max ), where λ min and λ max are the minimum and maximum eigenvalues of the standardized spatial weight matrix W π(ρ) ∝ 1.
The joint prior distribution of all the unknown quantities are presented by where τ = (τ 0 , τ 1 , • • • , τ q ) is a parameter vector that contains all the unknown hyperparameters for computational convenience.

The Full Conditional Posterior Distributions of the Latent Variables
According to the likelihood function (9) together with a standard exponential density, we can derive the full conditional posterior distributions of latent variables e i for i = 1, • • • , n under the condition of observation data (y, x, z, u) and the remaining unknown quantities, as follows where 11) is the kernel of a generalized inverse Gaussian distribution, we infer where the probability density function of and K v (•) is a modified Bessel function of the third kind [53].There exist efficient algorithms to simulate from a generalized inverse Gaussian distribution [54] so that our Gibbs sampler can be easily applied to quantile regressive estimation.

The Full Conditional Posterior Distributions of the Parameters
In this section, because the joint posterior of the parameters is complicated and it is not easy to draw samples directly, we propose a hybrid Gibbs sampler [55], also derive the full conditional posterior of all parameters and describe the detailed sampling procedure.
We can obtain the conditional posterior distribution of the spatial autocorrelation coefficient ρ from the likelihood function (9), which is proportional to p(ρ|y, x, z, u, e, β, γ, µ, σ 0 , τ) However, it is difficult to directly simulate from (12) without the form of any standard density function.We apply the Metropolis-Hastings algorithm [56,57] .

Numerical Illustration
In the section, Monte Carlo simulations are implemented to demonstrate the finite sample performance of the proposed model and estimation method.We also apply to analyse a real dataset example.In order to ensure the robustness and applicability, two kinds of matrices are chosen to investigate the spatial influence of the spatial weight matrix W on the estimation effects.One is the Rook weight matrix as [35], and the Rook weight matrix is generated according to Rook contiguity, which allocates the n spatial units on a lattice of m × m (≥ n) squares and finds the neighbours for unit with row normalizing.The other is the Case weight matrix as in [59], we consider the spatial scenario with r districts and m members in each district, and each neighbour of a member in a district is given equal weight.

Simulation
The samples are generated from the following model: where the covariate vectors x i = (x i1 , x i2 ) follows a bivariate normal distribution with mean vector 0 and covariance matrix , where F is the common cumulative distribution function of i ∼ N(0, 1).By subtracting the τth quantile, the error term is equal to zero at the τth quantile.The varying coefficient functions α(u) = (α 1 (u 1 ), α 2 (u 2 )) with α 1 (u 1 ) = 2 cos(2πu 1 ) + 1 and α 2 (u 2 ) = 0.5 exp{−2(2u 2 − 1) 2 } + 2u 2 .Furthermore, we chose three different values of spatial parameters ρ = {0.2,0.5, 0.8} at three different quantile points τ = {0.25,0.5, 0.75}, two kinds of the Rook and Case weight matrix as the spatial weight matrix W, respectively.The sample sizes are n = {100, 400} for the Rook weight matrix, districts and members are (r, m) = {(20, 5), (80, 5)} for the Case weight matrix.We conducted each simulation with 1000 replications.For j = 1, • • • , p, we use a quadratic P-splines in which the number of knots K j = 18 are placed at equally spaced interval of the predictor variables and design hyper-parameters (r 0 , s 2 0 , r τ 0 , s 2 τ 0 , r τ j0 , s 2 τ j0 ) = (1, 1, 1, 0.005, 1, 0.005) in our computation.The second-order random walk penalties are used for the Bayesian P-splines to approximate the unknown smooth functions.The unknown parameters are drawn from their respective prior distributions.The tuning parameter σ ρ is used to control the resultant acceptable rate for parameter around 25% by incrementally increasing or decreasing value.
We generated 6000 sampled values following the proposed Gibbs sampler and deleted the first 3000 values as a burn-in period for each of the replications until the Markov Chains reach steady state.According to the last 3000 values, we calculate the corresponding means across 1000 replications for the posterior mean (Mean), standard error (SE) and 2.5th and 97.5th percentiles of the parameters, namely the 95% posterior credible intervals (95% CI), which are defined by the posterior probability of the parameters falling into the intervals is 95% based on the highest posterior density.
We also computed the standard derivations (SD) of the estimated posterior means to compare them with the means of the estimated posterior SE.From the model (19), LeSage and Pace [60] suggested scalar summary measures for the marginal effects, which are given by ∂y ∂x j = (I n − ρW) −1 I n β j for j = 1, • • • , q.The direct effects are labeled as the average of the diagonal elements.The average of either the row sums or the column sums of the non-diagonal elements are used as the indirect effects, and the total effects are the sum of the direct and indirect effects.
To check the convergence of the MCMC algorithm, five different Markov Chains corresponding to different starting values have been ran through the Gibbs sampler to perform each replication.Figure 1 displays the sampled traces of parts of the unknown quantities, including model parameters and fitting functions on grid points.It is clear that the five parallel sequences mix reasonably well.We further calculate the "potential scale reduction factor"

√
R for all unknown parameters and varying coefficient functions on 10 selected grid points based on the five parallel sequences.Figure 2 shows the values of √ R after iterating 3000 times.We observe that all the values of √ R are less than 1.2 following the suggestion of Gelman and Rubin [61] after 3000 burn-in iterations, which is sufficient for convergence.In order to investigate the finite sample performance of varying coefficient functions, the variability measures of the mean absolute deviation errors (MADE) and global mean absolute deviation errors (GMADE) are used to measure the estimation performance.MADE and GMADE are defined as  We can see that the MADE and GMADE values not only decrease when the number of n increase but also become smaller under the Case weight matrix than the Rook weight matrix, meaning the varying coefficient functions become more accurate when increasing the sample size with application of the Case weight matrix.This shows that the proposed model and estimation method with both the Rook weight matrix and the Case weight matrix in the finite sample can obtain reasonable estimation and good performance.Tables 1 and 2 summarize the estimation results.The parameter estimates are quite different at three quantiles of the response distributions.Under the same spatial weight matrix, the accuracy of the results improves with the increasing of the sample sizes.We can see that the means of the unknown estimators are close to the respective true values, and the average values of the SE are close to the corresponding SD, indicating that the parameter estimates and the standard errors are more precise.For the parameter ρ under the same sample sizes, we find the SE and SD of parameter ρ with the Case weight matrix are slightly better than that with the Rook weight matrix.In addition, the general pattern from the estimates reported in Tables 1 and 2 is that all estimators impose relatively larger bias on the total effect estimates when there is strong positive spatial dependence for similar sample sizes.When we repeat the aforementioned experiences with different starting values, the estimation results are similar, all of which indicate that the proposed Gibbs sampler performs quite well.
Figure 4 compares the estimation results of varying coefficient functions at different quantiles, along with its 95% pointwise posterior credible intervals of α 1 (u) and α 2 (u) from a typical sample under (ρ, n) = (0.5, 100) and (ρ, n) = (0.5, 400), respectively.The typical sample is selected in such a way that its MADE value is equal to the median in the 1000 replications.We can see that the three fitting curves are fairly close to the solid curve, and the corresponding credible bandwidth is narrow.With the increasing of the sample sizes, the gaps between the fitting curves and the true value become short.There also exist visible differences at different quantiles of the response distributions.It illustrates that the varying coefficient function estimation procedure works well for small samples.
We compare the performance of the Bayesian quantile regression (BQR) estimator in this paper to the instrumental variable quantile regression (IVQR) estimator in Dai et al. [12] with two examples.
Example 1.The model is given as follows , F is the common cumulative distribution function of i , and the τth quantile of random error i is centred to zero.x i and u i are generated from N[0, 1] and U[0, 2], z i = (z i1 , z i2 ) are bivariate.z i1 and z i2 are generated independently from U[−2, 2] and N(1, 1).Table 3 summarizes the comparison results of QR, IVQR and BQR estimators with a homoscedastic error term.
The spatial weight matrix W = (w ij ) is generated based on mechanism that w ij = 0.3 |i−j| for i, j = 1, • • • , n, and then standardized transformation is applied to convert the matrix W to have row-sums of unit [12].After repeating the estimation procedure 1000 times for each case, we calculate the Bias and RMSE between the parameter estimates and true values, the MADE of the estimation accuracy of the varying coefficient functions.Tables 3 and 4 report the results of QR, IVQR and BQR corresponding to example 1 and example 2. It can be seen that the influence of explanatory variables on the response is quite different at different quantiles of the response distributions.When the sample sizes enhance, all the bias, RMSE and MADE of the estimators will decrease significantly.Comparing with the three methods QR, IVQR and BQR, the BQR estimator can obtain more robust results in the same condition with less bias, RMSE and MADE.We think that BQR algorithm is superior to QR and IVQR, although the later can also achieve reasonable estimations.

Application
As an application of the proposed model and methods to a real data example, we use the well-known Sydney real estate data with detailed description in [62].The data set contains 37,676 properties sold in the Sydney Statistical Division (an official geographical region including Sydney) in the calendar year of 2001, which is available from HRW package in R. We focus on the last week of February only to avoid the temporal issue including 538 properties.
In this application, the house price (Price) is explained by four variables, which are the distance from house to the nearest coastline location in kilometres (DC), distance from house to the nearest main road in kilometres (DR), inflation rate measured as a percentage (IR) and average weekly income (Income).The DC and DR have linear effects on the response Price, while the IR and Income have nonlinear effects on the response Price.Moreover, we make Price and DC logarithmic transformation to avoid the trouble caused by big gaps in the domain.In addition, Income is transformed so that the marginal distribution is approximately N(0, 1).Therefore, the following partially linear varying coefficient spatial autoregressive model will be developed: where the response variable Regarding the choice of the weight matrix, according to the practice in Sun et al. [10], we use the Euclidean distance in terms of any two houses to calculate the spatial weight matrix W. The location is represented with longitude and latitude, denoted as s i = (Lon i , Lat i ).The spatial weight w il is For this dataset, we adopt quadratic P-splines and hyper-parameters (λ, r 0 , s 2 0 , r τ α0 , s 2 τ α0 , r τ β j 0 , s 2 τ β j 0 ) = (2, 1, 1, 1, 0.005, 1, 0.005) for j = 1, • • • , p.The tuning parameter σ ρ is used to control the acceptable rate for updating ρ around 25%.
We run the proposed Gibbs sampler five times with different starting values and generate 10,000 sampled values following a burn-in of 20,000 iterations in each run.Traces of parts of the unknown quantities are plotted in Figure 5, and the five parallel sequences aggregate very well.Based on the five parallel sequences, we further calculate the "potential scale reduction factor" √ R, which is plotted in Figure 6.It is clear that all the values of √ R are less than 1.2 after 20,000 burn-in iterations.The proposed estimators can realize excellent convergence effects applying to the actual data.
Table 5 lists the estimated parameters together with their standard errors and 95% posterior credible intervals.It shows that the estimation of the spatial coefficient ρ is 0.57 with the standard deviation SE = 0.003 at τ = 0.5 quantile, which means that there exists positive and significant spatial spillover effects for the housing prices of Sydney real estate.However, the spatial coefficient decreases with the increase of quantiles, when the house prices are lower, the spatial effects and the interaction between different regions become stronger.The coefficients of the two covariates log(DC) and DR are β1 = 0.6039 and β2 = 0.4033 at τ = 0.5 quantile, they also have promotional effects on housing prices at the other two quantiles, log(DC) and DR will play an important positive role with the house prices rising because the two parameters present an increasing trend at higher quantiles.Figure 7 presents the estimated varying coefficient functions together with its 95% pointwise posterior credible intervals, which includes three quantiles τ = 0.25, τ = 0.5 and τ = 0.75 by a dotted line, star line and forked line, respectively.The curves totally show an upward trend, especially when u becomes larger, the curves rise up more.This shows that the effect of covariate Income on the response has a U-shaped nonlinear relationship.More specifically, when the quantile at τ = 0.5, the varying coefficient function α(u) is greater than the other two, meaning that Income has a significant promoting influence in areas with higher housing prices.The empirical result confirms the robustness and practicability of the Bayesian P-splines method.

Summary
This article focused on studying Bayesian estimation and inference in a quantile regression of partially linear varying coefficient spatial autoregressive models with Psplines.This can analyse the linear and nonlinear effects of the covariates on the response for spatial data, reduce the high risk of misspecification of the traditional SAR models and avoid certain serious drawbacks of fully nonparametric models.
We developed Bayesian quantile regression of PLVCSAR models using the asymmetric Laplace error distribution, which can capture comprehensive features at different quantile points without strict restrictions.Moreover, we considered a fully Bayesian P-splines approach to analyse the PLVCSAR models and designed a Gibbs sampler to explore the full conditional posterior distributions.Compared with the QR and IVQR estimators in the same condition, our methodology obtained more robust and precise results.Finally, the proposed model and method with an application were used to analyse a real dataset.
In this article, we considered spatial data with homoscedasticity or heteroscedastic error term, which does not need any specification of error distribution.Although we used a partially linear varying coefficient SAR model, the other models, such as a partially linear single-index SAR model and partially linear additive SAR model can also be considered.In addition, we also need to study variable selection and model selection in a large sample.

Figure 1 .
Figure 1.Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities, where (a-c) are trace plots of parameters (ρ, β 1 , β 2 ), (d-h) are trace plots of the varying coefficient function α 1 (u), and (i-l) are trace plots of the varying coefficient function α 2 (u) (only a replication with (r, m) = (80, 5), and ρ = 0.5 is displayed).
Figure 3a displays the boxplots of the MADE and GMADE values with sample size n = 100 and ρ = 0.5 at τ = 0.5 quantile point.Based on the Rook weight matrix on the left three panels, the medians are MADE 1 = 0.2049, MADE 2 = 0.1976 and GMADE = 0.2048.Based on the Case weight matrix on the right three panels, the medians are MADE 1 = 0.1993, MADE 2 = 0.1952 and GMADE = 0.2028.Figure 3b shows the boxplots of the MADE and GMADE values with sample size n = 400 and ρ = 0.5 at τ = 0.5 quantile point.Based on the Rook weight matrix, the medians are MADE 1 = 0.1049, MADE 2 = 0.1164 and GMADE = 0.1126 on the left three boxplots.Based on the Case weight matrix, the medians are MADE 1 = 0.1019, MADE 2 = 0.1160 and GMADE = 0.1101 on the right three boxplots.

Figure 2 .Figure 3 .
Figure 2. The "potential scale reduction factor" √ R for simulation results (the case of spatial parame-

Figure 4 .
Figure 4.The estimated varying coefficient functions (dotted line at τ = 0.25 quantile, star line at τ = 0.5 quantile and forked line at τ = 0.75 quantile) and their 95% pointwise posterior credible intervals (dot-dashed lines) for a typical sample (the left panels are based on the Rook weight matrix and the right panels are based on the Case weight matrix with ρ = 0.5).The solid lines denote the true varying coefficient functions.

Figure 5 .
Figure 5. Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities, where (a-c) are trace plots of parameters (ρ, β 1 , β 2 ), and (d-i) are trace plots of the varying coefficient function α 1 (u).

Figure 6 .
Figure 6.The "potential scale reduction factor" √ R for Sydney real estate data.

Figure 7 .
Figure 7.The estimated function (dot line) and its 95% pointwise posterior credible intervals (dotdashed lines) in the model (20) for Sydney real estate data.

Table 3 .
Simulation results of the parameter estimation.

Table 4 .
Simulation results of parameter estimation.

Table 5 .
(20)meter estimation in the model(20)for Sydney real estate data.