1 Introduction

The class of linear quantile autoregression (QAR) models of Koenker and Xiao (2006) has proven to be particularly useful for studying asymmetric dynamics and local persistence in many economic and financial time series. Examples include real estate prices (Galvao, Montes-Rojas & Park, 2013; Lee, Lee & Lee, 2014), stock prices and exchange rates (Ferreira, 2011; Baur, Dimpfl & Jung, 2012; Yang, Tu & Zeng, 2014), and more general business cycle and inflation indicators (Manzan, 2015). Our objective in this paper is to extend the class of QAR models by allowing for the possibility of Markov-switching regimes that would influence the conditional quantiles. Indeed, the presence of regime changes in the conditional distribution of the response variable would obviously affect its quantiles. The importance of allowing for such effects can also be recognized from the work of Qu (2008) who proposes testing procedures for structural change in conditional quantiles.

Since the seminal contribution of Hamilton (1989), Markov-switching specifications have become an immensely popular approach to model structural changes in the behaviour of economic and financial time series; see Piger (2009) for an overview. A distinctive feature of the Markov-switching approach is that the regime changes are endogenous to the model, as opposed to being treated exogenously like in the classic approach to structural changes by Chow (1960) which assumes a priori knowledge about how to classify the data between regimes. Building on the earlier work of Goldfeld and Quandt (1973), the Hamilton (1989) approach specifies the regimes as being determined by a discrete-time, discrete-state Markov chain with unknown transition probabilities. This process is assumed to be recurrent, meaning that it can move from one state to any other state at any time. As in the well-known Kalman filter, the unobserved state can be inferred from the observable time series and the sample likelihood function can then be recursively computed. Our specification also complements the related extension of the QAR model proposed by Galvao, Montes-Rojas, and Olmo (2011) in which regime changes are triggered when the time series passes (quantile-specific) threshold values, like in the self-exciting threshold autoregression models of Tong (1983).

Following Albert and Chib (1993), McCulloch and Tsay (1994), Chib (1996), and Frühwirth-Schnatter (2001), we adopt a Bayesian approach to inference based on data augmentation such that the latent states can be analyzed along with the other unknown model parameters through Gibbs sampling. The advantage of the Gibbs sampling approach to the analysis of Markov-switching models has long been recognized. For example, Albert and Chib (1993) remark that such an approach avoids the direct calculation of the likelihood function needed for maximum likelihood estimation. Moreover, the posterior distributions of all unknown parameters (and functions thereof) are fully tractable and easy to simulate from. By treating the unobserved states as missing data, this approach also provides posterior distributions of the states, integrated (by marginalization) over all the other model parameters.

The Bayesian analysis of quantile regression models rests on the connection between the quantile estimation problem and the likelihood function under an asymmetric Laplace distribution, established by Yu and Moyeed (2001) and Tsionas (2003). It is important to note, however, that as in Koenker and Machado (1999) the asymmetric Laplace distribution is not assumed for the data generating process, but is merely used to obtain a quasi-likelihood function. Other examples of such a parametrization for the purpose Bayesian inference include Geraci and Bottai (2007), Kottas and Krnjajić (2009), Yue and Rue (2011), Gerlach, Chen, and Chan (2011), and Chen et al. (2012).

The asymmetric Laplace distribution can be expressed as a location-scale mixture of exponential and normal distributions (Kotz, Kozubowski & Podgórski, 2001). Following Kozumi and Kobayashi (2011), we exploit this representation to develop a Gibbs sampling approach wherein all complete conditional densities have closed-form expressions. An estimate of the marginal likelihood can then be calculated from the Gibbs output using the method of Chib (1995). The marginal likelihood is the key ingredient for model selection and discrimination in Bayesian statistics; see the discussion in Kass and Raftery (1995).

We use the Gibbs sampler to solve the well-known quantile crossing problem that may arise when quantile models are fitted separately for each considered quantile probability level, τ. This common practice of treating the quantile levels independently of one another can yield fitted quantile curves that cross one another, thereby leading to a nonsensical response distribution. Indeed, quantile monotonicity requires the quantiles to be increasing as a function of τ, meaning that any well-defined distribution must necessarily have non-crossing quantiles. Koenker and Xiao (2006), §4 remark that the crossing problem appears more acute in QAR models than in ordinary quantile regressions with exogenous covariates, since the support of the regressors is determined within the autoregressive model.

This quantile crossing problem is potentially worse for non-linear quantile autoregression models, like the threshold specification of Galvao, Montes-Rojas, and Olmo (2011) and the Markov-switching specification developed here. Fortunately, the Gibbs sampler provides the basis for a stepwise re-estimation procedure that ensures non-crossing quantiles. As in Gelfand, Smith, and Lee (1992), draws from the constrained posterior distribution are obtained straightforwardly by retaining the Gibbs draws that satisfy the non-crossing condition when sampling the unconstrained posterior distribution. As far as we know, this is the only way to carry out full Bayesian calculations while avoiding well-nigh impossible numerical integrations over high-dimensional sets defined by complex restrictions involving the model parameters and the data.

It is important to note that the proposal in a related paper by Ye et al. (2016) is to allow for Markov-switching parameters in an ordinary quantile regression with exogenous (independent) regressors, not a quantile autoregression. This is a fundamental difference with what we propose here. Indeed, as we already mentioned, QAR models are more likely to suffer from the quantile crossing problem than ordinary quantile regressions since the regressors are endogenous to the model. Ye et al. (2016) also exploit the asymmetric Laplace connection to devise a (quasi) maximum likelihood estimation (MLE) approach for point estimation, but they do not establish any distributional theory to guide inference. Moreover, their MLE-based approach treats each quantile level separately, which means that it may yield crossing quantiles even though their model is less prone to this problem. In sharp contrast to Ye et al. (2016), our Gibbs sampling approach offers a complete Bayesian methodology not only for estimation, but also inference, model selection, and ensuring non-crossing quantiles in quantile regresison models. It is also important to note that our complete closed-form Gibbs sampling approach can be applied in any quantile regression model with endogenous or exogenous covariates, and whether Markov-switching effects are allowed for or not. In our empirical application, for instance, we estimate non-crossing quantiles with the linear QAR model as well as the Markov-switching QAR model.

The current paper is organized as follows. Section 2 begins by introducing the QAR models of Koenker and Xiao (2006), then shows the asymmetric Laplace connection, and describes the proposed Markov-switching quantile autoregression models. Section 3 develops the Gibbs sampling algorithm based on a location-scale mixture representation of the asymmetric Laplace distribution. There are two variants of the approach, depending on how the state variables are sampled (single- or multi-move Gibbs sampling). Section 4 explains the computation of the marginal likelihood used for model comparisons and the stepwise re-estimation procedure to ensure non-crossing quantiles. Section 5 presents some simulation evidence about the relative performance of the Gibbs samplers in the quantile regression context. Section 6 presents the empirical application to the U.S. real interest rate, which illustrates the gains obtained by enforcing the quantile monotonicity restriction. Section 7 concludes and the computational details are given in the appendices.

2 Markov-switching quantile autoregression

Suppose we have a response variable of interest y_t whose time-t conditional quantiles we wish to model as a function of past information. The linear quantile autoregression (QAR) model of Koenker and Xiao (2006) specifies the τth conditional quantile of y_t as

where yt1:t2 refers to observations y_t1, …, y_t2, for t₁ ≤ t₂, with the convention that y=y1:T. Given a specified value of the quantile level τ ∈ (0, 1), the parameters of the QAR(p) model in (1) can be estimated by solving the following minimization problem:

by choice of c(τ), ϕ₁(τ), …, ϕ_p(τ), and where ρ_τ(⋅) is the asymmetric absolute deviation loss function defined as ρτ(u)=u(τ−I[u<0]) (Koenker & Bassett, 1978). Here 𝕀 [A] is the indicator function which equals one when event A is true, and zero otherwise. For the median τ = 0.5, the loss function in (2) becomes ρ_τ(u) = 0.5 ∣ u ∣.

Koenker and Machado (1999) explain that the parameters of linear quantile models can also be estimated by (quasi) maximum likelihood. To see how, consider the QAR(p) model specified in parametric distributional form as

where δ > 0 is a scale parameter and {ε_t} are i.i.d. according to the asymmetric Laplace distribution with probability density function

The conditional density function of y_t for a given τ then becomes

and, conditional on y1:p, the sample likelihood can be written as

Since the negative of the loss function appears in the exponent of this expression, maximization of (5) is equivalent to solving the minimization problem in (2). As Yu and Moyeed (2001) and Tsionas (2003) explain, the asymmetric Laplace distribution provides a natural pathway for the Bayesian analysis of quantile regression models. Note also that the value of δ does not matter for the estimation of the correct quantiles by maximizing (5). But rather than fixing it to a constant value, using δ as a free scale parameter clearly makes the assumed asymmetric Laplace distribution more flexible to capture the true (unknown) error distribution. Other examples of such a parametrization for Bayesian inference include Geraci and Bottai (2007), Kottas and Krnjajić (2009), Yue and Rue (2011), Gerlach, Chen, and Chan (2011), and Chen et al. (2012).

The starting point for the developed model is the observation that the QAR(p) model in form (3) is equivalent to

where {ε_t} are i.i.d. according to the asymmetric Laplace distribution with density (4) so that μ(τ) is the location of the τth conditional quantile in this model with autocorrelated errors. Indeed, since η_t = y_t − μ(τ), for all t, (6) can be rewritten as (3) with c(τ)=μ(τ)(1−∑j=1pϕj(τ)). The presence of regime changes in the true conditional distribution of y_t would obviously affect the location of its quantiles. Following Hamilton (1989), the quantile model in (6) can be generalized to account for the possible presence of such regimes as follows:

where s_t is a latent variable taking values in the set {1, …, K} according to a Markov chain with one-step transition probability matrix P whose elements are defined as

The term μ(τ, s_t) in (7) is thus the location of the τth conditional quantile of y_t given the past of y_t itself and the current regime s_t. The proposed Markov-switching quantile autoregression (MSQAR) model can then be rewritten as

since η_t = y_t − μ(τ, s_t), for all t, in (7). For a specified τ, the conditional density function of y_t given y_t−1, …, y_t−p and s_t, …, s_t−p becomes

where the τth conditional quantile function is given by

Here st1:t2 refers to s_t1, …, s_t2, for t₁ ≤ t₂ , and we let s=s1:T. Observe that it would not make sense to allow for switching in δ, since the quantile regimes cannot be identified through this reduced-form parameter. Note also that if μ_i(τ) = μ_j(τ) for all i, j, then the K-regime MSQAR(K, p) model in (11) collapses to the linear (K = 1 regime) QAR(p) model in (1). In this case we have μ(τ)=c(τ)/(1−∑j=1pϕj(τ)), assuming of course that there are no unit roots so the denominator is different from zero.

The MSQAR model in form (9) can be expressed more compactly as

where ϕ(τ,L)=(1−ϕ1(τ)L−⋯−ϕp(τ)Lp) is a pth order polynomial in the lag operator L, defined such that L^k z_t = z_t−k for k ≥ 0. The following assumptions are made:

1. μ₁(τ) < μ₂(τ) < ⋯ < μ_K(τ).

2. All the roots of ϕ(τ, L) = 0 lie outside the unit circle.

3. p_ij > 0 for all i, j ∈ {1, …, K}.

4. Pr(s₁ = i) = 1/K for all i ∈ {1, …, K}.

The first three assumptions are standard for Markov-switching time-series models. Assumption 1 ensures the identification of the regimes, while Assumption 2 imposes stationarity given the sequence of state variables. Assumption 3 guarantees that the Markov chain is irreducible, i.e. starting s_t from an arbitrary state i ∈ {1, …, K}, any state j ∈ {1, …, K} is reachable in finite time. This assumption is also needed for identification purposes because if a state is never visited then the associated parameters cannot be identified. Assumption 4 treats the initial state as a random draw from a uniform distribution, independently of the transition probabilities; see Frühwirth-Schnatter (2006), Ch. 10 for a discussion of the computational advantages of this assumption.

Under these assumptions, the joint density of yp+1:T and s, conditional on y1:p, is

which does not constitute the likelihood function. Indeed, the likelihood function for yp+1:T is obtained by integrating out the state variables. Appendix A presents an algorithm to compute the MSQAR likelihood function.^[1]

3 Gibbs sampling

The asymmetric Laplace distribution admits various mixture representations (Kotz, Kozubowski & Podgórski, 2001). Following Kozumi and Kobayashi (2011), the Gibbs sampling algorithm developed here uses a representation based on a mixture of exponential and normal distributions. Specifically, if ε_t follows the asymmetric Laplace distribution with density (4), then ε_t can be represented as

with

where et∼E(1), a standard exponential distribution, and z_t ∼ N(0, 1), independently of e_t. With this mixture representation, model (9) can be equivalently rewritten as

However, as Kozumi and Kobayashi (2011) explain, such an expression is not convenient for Gibbs sampling because the scale parameter δ appears in the conditional location of y_t. This issue is resolved by working instead with the reparameterization:

where v_t = δe_t and we let v=vp+1:T=(vp+1,…,vT).

The model parameters are collected in θ(τ)=(μ(τ)′,ϕ(τ)′,δ,p′)′, where μ(τ)=(μ1(τ),…,μK(τ))′, ϕ(τ)=(ϕ1(τ),…,ϕp(τ))′, and p=vec(P) stacks the elements of the transition probability matrix into a column vector. The Gibbs sampler is an iterative procedure that can be started using any set of values in the support of the posterior distribution π(θ(τ) | y) and will produce a Markov chain {θ(τ)n}n=1N with equilibrium distribution equal to this target posterior distribution; see Casella and George (1992) for an introduction and Tierney (1994) for a more detailed treatment. Note that v and s are latent, so data augmentation will be used when sampling from the posterior distribution. Specifically, given a previous Gibbs draw (sn−1,pn−1,μ(τ)n−1,ϕ(τ)n−1,vn−1,δn−1), the next one is obtained by sampling iteratively from the following full conditional distributions:

where Cμ(τ)={μ(τ):μ1(τ)<μ2(τ)<⋯<μK(τ)} and Cϕ(τ)={ϕ(τ):all the roots of ϕ(τ,L)=0 lie outside the unit circle} are the constraint sets for μ(τ) and ϕ(τ) under Assumption 1 and Assumption 2, respectively. These steps are repeated a large number of times until the Markov chain has achieved convergence. After the burn-in iterations, each complete pass through the Gibbs steps yields a draw from the joint posterior density π(s,p,μ(τ),ϕ(τ),v,δ | y). These draws from the sampling algorithm are denoted by {sn,pn,μ(τ)n,ϕ(τ)n,vn,δn | y}n=1N. Appendix B explains in detail how to generate draws from each of the conditional distributions making up the steps of the Gibbs sampling procedure.

4 Specification issues

The proposed MSQAR model applies to a specified quantile probability level τ. A typical quantile regression analysis, however, might involve several probability levels τ₁ < ⋯ < τ_q. Two issues arise when several quantile probability levels τ are considered. The first is that the models at the considered quantile levels each yield an inference about the latent Markov chain. This begs the question: which one should be retained?

The second issue is that the fitted quantiles can cross one another, since the models (defined for each τ) are fitted separately. In other words, if the models are estimated separately for each of the q desired probability levels, then the resulting conditional quantile functions may not be monotonically increasing in τ. This is the well-known quantile crossing problem, which obviously leads to a nonsensical response distribution since any distribution must necessarily have non-crossing quantiles in order to be well defined. We solve this problem by proposing a stepwise procedure similar to Wu and Liu (2009), whereby the quantiles are refitted sequentially while constraining the current curve not to cross the previous one.

5 Monte Carlo experiments

In this section, we examine the performance of the Gibbs sampler by means of Monte Carlo experiments. We use as data-generating process (DGP) the Garcia and Perron (1996) Markov-switching AR(2) model, given as

where μ(st)=∑i=13μiI[st=i], σ(st)=∑i=13σiI[st=i], and s_t takes values in {1, 2, 3} according to the outcome of a first-order Markov chain. This three-state Markov-switching model was used by Garcia and Perron (1996) to investigate the presence of regime changes in the conditional mean and variance of the quarterly U.S. real interest rate from 1960Q1 to 1990Q4 (124 observations).

We set the true values of the parameters appearing in (14) as μ=(−1.5,1.3,4), ϕ=(0.05,0.05), σ2=(5.5,1.5,6.5), and the transition probabilities of the Markov chain are set as p_ij = 0.95 when i = j, and p_ij = 0.025 when i ≠ j. These values correspond closely to the estimates reported by Garcia and Perron (1996) and Kim and Nelson (1999), §9.3. We examine sample sizes T = 120, 240 and three different distributions for ϵ_t: (i) a standard normal distribution; (ii) a Student-t distribution with three degrees of freedom, standardized to have unit variance; and (iii) a gamma distribution with shape 4 and scale 1, standardized to have mean zero and unit variance.

We consider the MSQAR(K, p) model in (9), correctly specified with K = 3, p = 2, and nine different quantile levels τ = 0.1, …, 0.9. We also include a misspecified QAR(2) to assess the effect of ignoring the Markov-switching structure. The prior for μ(τ) is a trivariate normal with mean μ+Qzt(τ) and covariance matrix equal to 0.12 times a diagonal matrix, where Q_zt(τ) is the standard normal quantile function. A bivariate normal distribution with mean zero and covariance matrix equal to 0.08 times a diagonal matrix is used as the prior for ϕ(τ), at all τ. Finally the prior for the scale parameter δ appearing in the asymmetric Laplace distribution is an IG(0.1/2, 0.1/2) distribution, and the transition probabilities are sampled according to Appendix B.2 using a Dirichlet distribution with hyperparameters set to 0.1. After running the Gibbs (single- and multi-move) samplers for 5000 burn-in iterations, we use the next 20,000 draws with a thinning of 2 for inference, and the simulation results are based on 400 replications of each DGP configuration.

From the Gibbs output we obtain s^1,…,s^T, the states identified a posteriori. As explained in Section 4.2 these are found from the posterior state probabilities Pr(st=i | y), which are computed by averaging the Gibbs draws of the state variables. The Gibbs sampler can then be evaluated by computing the proportion of correctly classified (PCC) states as

where s_t is the true state at time t.

Table 1 reports the median along with the corresponding lower 5% and upper 95% quantiles of the PCC statistics across the 400 replications of each DGP configuration. We see that 84 to 94% of states are correctly identified when the quantile level τ is between 0.4 and 0.6, and the PCCs deteriorate as the quantile level becomes more extreme towards either tail. Comparing the distributions, we observe that the best performance is achieved under Student-t errors. For instance, when T = 120 and τ is in the 0.4–0.6 range, the median PCCs under normal and gamma errors are around 80%, while those under the heavier-tailed Student-t errors exceed 90%. It is also interesting to note that increasing the sample size does not affect much the median PCC, but rather has a greater effect of the range of the estimated PCCs. As T doubles from 120 to 240 under each distribution, we can see in general a narrowing of the range between the lower and upper PCC quantiles, while the median remains almost unchanged. This is the same effect also seen when comparing the single- and multi-move samplers. Indeed, the relative computational efficiency of the multi-move relative to the single-move samplers appears most importantly as a reduction in the variance of correctly identified states.

The accuracy of estimation is further assessed by examining the deviations between the true conditional quantile Qyt(τ | yt−1,yt−2,st,st−1,st−2;θ) implied by the DGP in (14) and Qyt(τ | yt−1,yt−2,s^t,s^t−1,s^t−2;θ^(τ)), the estimated conditional quantile under the MSQAR(3, 2) model. More precisely, we compute the mean absolute deviation error (MADE)^[4] across observations as

where θ^(τ) are the posterior means of the parameters, and the MADE statistic is computed for each DGP replication, as we did before with the PCC statistic. The results are shown in Table 2, where for each DGP configuration the entries are the median of the 400 MADEs reported with the corresponding lower 5% and upper 95% quantiles.

We see that using the asymmetric Laplace distribution as the likelihood achieves the greatest estimation precision as defined by (16) for the central quantiles under Student-t errors. In this case, the MADEs are around 0.30, while they exceed 0.50 under normal and gamma distributed errors, and sometimes by far. As expected, the MADEs appear roughly symmetric under normal and Student-t errors, i.e. the estimated MADEs are about the same whether τ equals 0.1 or 0.9, 0.2 or 0.8, etc. On the contrary, the MADEs for the gamma distribution are consistently higher in the right tail, e.g. when τ = 0.9. This happens because the distribution is skewed to the right, meaning that there are relatively fewer observations in the right tail compared to the left tail.

Of course, increasing the sample size improves the estimation precision at all quantile levels τ. This is readily seen from the upper 95% MADE quantiles which tend to decrease as T doubles from 120 to 240. Observe further that the multi-move sampler appears preferable, since these upper 95% limits are systematically lower than with the single-move sampler. We therefore leave aside the single-move sampler and proceed to the empirical application with the computationally more efficient multi-move sampler.

Table 3 provides a comparison of the in-sample MADEs of the misspecified QAR(2) model, obtained by Gibbs sampling and MLE.^[5] The pattern is clear: even with the misspecified QAR(2) model, the Bayesian approach yields smaller in-sample median MADEs relative to the MLE approach. Focusing on the averages, we see that the pattern holds for each T and error distribution.

Table 4–Table 6 provide an assessment of the out-of-sample accuracy of the model forecasts. For this purpose we use the DGP in (14) to simulate trajectories of length T + 1, and then, using only the data up to time T, we forecast the quantiles of y_{T + 1}. Table 4 reports the average of I[s^T+1=sT+1] across DGP replications. The multi-move sampler achieves the highest out-of-sample proportion of correctly identified states across sample sizes and error distributions. For instance, under normal errors, the multi-move sampler is able to identify 74.9% of states at time T + 1 = 241 (on average across quantiles), which is higher than the single-move sampler (69.4%) and the MLE approach (71.3%).