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.
48,362 characters · 11 sections · 38 citation commands
Bayesian Median Autoregression for Robust Time Series Forecasting
Time series forecasting is a long-standing problem in econometrics and statistics, where the overwhelming focus has been on mean-based models prado2010time, hyndman2018. Although conditional means are the optimal forecast under the squared error loss, complex characteristics violating model assumptions that are present in real data may hamper the predictive performance. Flexible nonparametric methods ferraty2006nonparametric, fan2008nonlinear building on minimal assumptions are an appealing remedy; however, they often have disadvantages in terms of not only interpretability but also in involving large numbers of parameters that often lead to daunting computation and communication issues. We were thus motivated by an attempt to propose a simple, interpretable, and principled strategy that can improve upon mean-based models in real-world time series forecasting.
There is a rich literature on robust time series forecasting including categorizing outliers Fox1972, HR2014, adjusting autoregressive (AR) models to offset effects of outliers CL1993, CLJ1993, exponential smoothing and Holt-Winters seasonal methods to M-estimation CGF2008, weighted forecasts JW2008, and detecting structural changes Q2008, OQ2011, CP2019, just to name a few. As a nonlinear alternative, the median is more robust than the mean. While the application of median-based methods in time series at least dates back to 1974 when John Tukey introduced the running median method JWT1974, there is surprisingly little work to comprehensively investigate the modeling, fitting, and uncertainty quantification of a median-based model in the context of time series forecasting, and in particular, how it compares with state-of-the-art mean-based methods using real data and under various loss functions.
In the quantile regression literature that encompasses the median as a special case KB1978, KX2006 proposes quantile autoregression (QAR) models which depict the conditional distributions of the response more comprehensively at various quantile levels. EM2004 proposed the conditional autoregressive value at risk (CAViaR) model for risk management at extreme quantile levels, and this model has been followed by many others GK2005, CS2006, KG2007, GCC2011. However, they were developed either for general quantile levels or extreme quantile levels and have not been applied to time series forecasting. Moreover, QAR and CAViaR are semiparametric approaches, resorting to minimizing the check loss function for estimation and necessitating non-trivial modifications to likelihood-based order selection criteria when the order is unknown.
In this paper, we propose a simple strategy by extending the traditional AR model to a median AR model (MAR) for time series forecasting. The AR model is arguably one of the most popular methods in time series, serving as the building block for other models such as generalized autoregressive conditional heteroskedasticity models bollerslev1986 and time-varying vector autoregressive models primiceri2005time. The proposed method utilizes time-varying quantile regression but focuses on the median, favorably inheriting the robustness of median regression in contrast to the widely used mean-based methods. It relies on parametric assumptions, and this aids interpretation and enables convenient uncertainty quantification and propagation through a principled Bayesian framework. Numerical experiments using U.S. macroeconomic data show that this simple MAR approach leads to favorable and often superior predictive performances compared to selected state-of-the-art mean-based methods that are much more complicated in nature. It is remarkable that the comparison is made not only under the absolute but also the squared error loss for point forecasts and is extended to probablistic forecasts. The proposed methods are generic and can be used to complement any methods that build on the AR models by altering the Gaussian error assumption therein.
The MAR model has a close connection with the working asymmetric Laplace likelihood approach in Bayesian quantile regression. The working likelihood formalized by YM2001 has recently gained increasing attention GCC2011, SK2016, BD2018, LLM2018. It provides a principled and convenient framework to quantify uncertainties in the parameter estimation eliminating the challenging task of estimating the unknown conditional density functions that are required by the conventional quantile regression for inference yang2016posterior. Although the asymmetric Laplace likelihood is generally not the true data generating likelihood, a pragmatic view to support its use is that the maximum a posterior estimates resemble the usual quantile regression estimates which optimizes the check loss function. Theoretically, the posterior distribution of parameters concentrates on what minimizes the Kullback-Leibler divergence with respect to the true data-generating models kleijn2006. Unlike quantile regression where one may vary the quantile levels, the median is the primary quantile level of interest in time series forecasting. Therefore, the MAR model uses the Laplace distribution as the likelihood, alleviating the concern that the working likelihood is not a valid likelihood when considering multiple quantile levels. This fully parametric model enables routine posterior sampling. We estimate model parameters by Markov chain Monte Carlo, and propose using a Bayesian model averaging (BMA) approach hoeting1999 to propagate the model uncertainty in the unknown autoregressive order in addition to a Bayesian model selection strategy.
The rest of the paper is organized as follows. Section (ref) introduces the MAR model, estimation, and forecasting procedure. In Section (ref) we conduct simulations to compare the proposed approach with other competitive methods to assess parameter estimation. Section (ref) consists of a variety of applications using real-world economic data. Section (ref) concludes the paper.
We propose a median autoregressive (MAR) model for time series forecasting, as an alternative to the widely used mean-based models. Suppose we observe a vector of time series $\bm{y} = (y_1, \ldots, y_T)$. A MAR model with order $p$, denoted as MAR$(p)$, first assumes a time-varying quantile regression structure
where $\bm{y}_{t-1} = (1, y_{t-1}, ..., y_{t-p})'$, $\bm{\beta} = (\beta_0, \beta_1, ..., \beta_p)^\prime$ is vector of unknown coefficient, and $\epsilon_t$ is random error with median 0. Model (ref) is semiparametric as the error distribution is left unspecified other than the constraint of possessing a zero median. The classic AR model with a given order $p$, denoted as AR$(p)$, assumes a Gaussian error distribution with mean 0 and standard deviation $\sigma$, i.e., $\epsilon_t \sim N(0, \sigma^2).$ The MAR model further assumes $\epsilon_t \sim \mathrm{Laplace}(0, 2\tau)$ whose probability density function is
where $\tau>0$ is a scale parameter. The Laplace error assumption combined with the semiparametric structure in Equation (ref) yielding the following likelihood function for the MAR model
which is parametric. The parametric assumption in (ref) aids interpretation and enables convenient uncertainty quantification and propagation through a principled Bayesian framework. In addition, the autoregressive structure in the MAR model resembles the widely used AR model provides enormous flexibility and potential to complement the rich literature that builds on the AR model.
The use of Laplace distributions is common in literature the Bayesian quantile regression where the $\theta$th quantile of the error in Equation (ref) is assumed to be 0. For general $\theta \in (0, 1)$, a working likelihood method adopts the asymmetric Laplace distribution $\mathrm{AL}(\mu, \tau, \theta)$ as the error distribution, which has the probability density function
where $\mu$ is a location parameter and $\mathbbm{1}(\cdot)$ is the indicator function. The asymmetric Laplace distribution reduces to the Laplace distribution at the median by setting $\mu = 0$ and $\theta = 0.5$. In view of this intimate connection with Bayesian quantile regression as well as the subsequent Bayesian estimation and prediction, we also refer to the MAR model with Laplace errors given by (ref) and (ref) as Bayesian MAR, or BayesMAR, and use it exchangeably with the MAR model.
We use an uninformative prior for $\bm{\beta}$ and the Jeffreys prior for $\tau$, namely,
The posterior distribution of $(\bm{\beta}, \tau)$ is
which is proper CH2013. The regression coefficients $\bm{\beta}$ are instrumental for time series forecasting, and we derive their marginal posterior distributions by integrating out $\tau$ for efficient sampling:
The posterior sampling of $\bm{\beta}$ proceeds by Markov chain Monte Carlo (MCMC) via the Metropolis-Hastings (MH) algorithm MRRT1953, H1970:
We tune the parameter $a$ in Step 1 such that the final acceptance rate is between $20\%$ and $50\%$ GRG1996. In addition, we have also implemented Gaussian and heavy-tailed student-$t$ proposals, which are recommended by GCC2011 when studying extreme quantile levels under the CAViaR model. We did not observe empirical advantages of using such proposals over a uniform proposal under the MAR model, suggesting that one may choose more flexible proposals for median regression. For all experiments, we use 40,000 MCMC samples with 25,000 burn-ins, initialize $\bm{\beta}$ randomly within the unit interval, and use the posterior mean as the Bayes estimate of ${\bm{\beta}}$. Trace plots indicate the MCMC samples converge quickly, mostly within thousands of iterations.
The proposed Bayesian approach provides a convinient method for density forecasting beyond point forecasts. The $1$-step ahead predictive density conditional on $\bm{y}_{1:T} = (y_1, \ldots, y_T)$ is
For large $T$ the integrand $\left[ \sum^{T}_{t=p+1} |y_t - \bm{y}'_{t-1}\bm{\beta}| \right]^{-(T-p+1)}$ may quickly decay to zero, leading to considerable numerical errors in direct evaluation of this integral. Alternatively, the MCMC samples of $(\bm{\beta}, \tau)$ allow a convenient sampling strategy to approximate the predictive density. In particular, we draw samples of $y_{T + 1}$ from the Laplace likelihood $p\left(y_{T+1} | \bm{\beta}, \tau , \bm{y}_{1:T}\right)$ conditional on each posterior sample of $(\bm{\beta}, \tau)$. This strategy easily generalizes to $q$-step ahead predictive densities for any $q \geq 2$ by drawing samples jointly for $(y_{T + 1}, \ldots, y_{T + q})$ through iterative conditional distributions, which are all Laplace distributions.
The order $p$ in MAR($p$) is typically unknown. We address the problem of unknown $p$ in the Bayesian framework by putting a prior on $p$. In practice, we can usually specify a maximum order; otherwise, a $p$ that is too large hampers the interpretability. We endow the order $p$ with a uniform prior on $\{1, 2, \ldots, K\}$ with $K$ being the specified maximum order. Then the posterior distribution of $p$ in the prior support is
where $ \pi(\bm{y} \mid p) = \int_{\mathbb{R}^+} \int_{\mathbb{R}^{p+1}} \pi(\bm{y}, \bm{\beta}, \tau \mid p) \pi(\bm{\beta}) \pi (\tau) d \bm{\beta} d{\tau}$ is the marginal likelihood of $p$. The order $p$ can be selected by using the maximum a posteriori (MAP) estimate:
For time series forecasting, a more appealing perspective is to use Bayesian model averaging (BMA) to propagate uncertainties in the model space, i.e.,
where $\hat{y}^{(p)}_{T+q}$ and $p(y_{T + q}^{(p)} | \bm{y}_{1:T})$ are the $q$-step ahead point prediction and predictive density of $y_{T + q}$ under order $p$, respectively.
The main challenge in implementing MAP and BMA lies in prior specifications of model parameters and the evaluation of the marginal likelihood at given $p$. To date, a consensus on the default choice for prior specifications in the context of model selection is still lacking, and one needs to be cautious about using improper priors---which are typically the default choice for a given model---in view of prior sensitivity and the Jeffreys-Lindley paradox. We refer interested readers to the rich literature on BMA, e.g., cp2014, R2001, YVSG2018, LD2020, to name a few. We here resort to an approximation to $\pi(p \mid \bm{y})$ using the Bayesian information criterion or BIC RA1995,AJ2012, which is appealing as it eliminates the need to deal with prior specification and approximates Bayes factors reasonably well in certain cases KW1995. We observe that the BIC-based implementation tends to choose the oracle order with large probability in simulations.
Letting $(\hat{\bm{\beta}}_{\text{MLE}}, \hat{\tau}_{\text{MLE}} )$ be the maximum likelihood estimates (MLEs) of $(\bm{\beta}, \tau)$, then the the BIC of MAR($p$) is
where $n$ is the sample size. We approximate $\pi(p \mid \bm{y})$ by $\exp\{-\text{BIC}_p/2\}$ up to multiplicative constants, leading to aggregated predictions
It turns out that we can calculate the MLEs $(\hat{\bm{\beta}}_{\text{MLE}}, \hat{\tau}_{\text{MLE}} )$ efficiently. To see this, first note
For any $\tau > 0$, the likelihood function $L(\bm{\beta}, \bm{\tau})$ in Equation (ref) attains its maximum at
provided $T > p$. This corresponds to the estimators of minimizing absolute error in median regression, which can be efficiently solved by linear programming Koenker2005. An analytical solution of $\hat{\tau}_{\text{MLE}}$ is available through a Gamma kernel:
Before substituting $\hat{\bm{\beta}}_\text{MLE}$ and $\hat{\tau}_{\text{MLE}}$ into Equation (ref), we notice that both the likelihood function and sample size depend on the order $p$. To reconcile various sample sizes at different orders, we use the last $T-K$ samples to evaluate the likelihood function for any $p$. Consequently, $\text{BIC}_p$ is given by
The same methods to approach the unknown order and provide predictive densities apply to AR, where the likelihood function changes to Gaussian and the MLE of $\bm{\beta}$ that resembles the least square estimates has a simple closed-form expression.
In this section, we conduct simulations to assess the performances of BayesMAR with mean-based methods, focusing on parameter estimation under various model assumptions. To this end, we choose the AR model and the generalized autoregressive conditional heteroscedasticity model (GARCH), and defer predictive comparisons and more recent mean-based methods to real data application in Section (ref).
We generate data according to the model
where $\bm{\beta} = (\beta_0, \beta_1, \beta_2) = (0.3, 0.75, -0.35).$ We consider two scenarios depending on the distribution of $\varepsilon_t$: Gaussian error where $\varepsilon_t \sim N(0,1)$ and Laplace error where $\varepsilon_t \sim \mathrm{Laplace}(0,1)$, which correspond to the model assumptions of the AR and BayesMAR models, respectively.
For each error assumption, we generate 200 observations and replicate such simulation 100 times. We estimate AR models via a Bayesian procedure with priors $\pi(\bm{\beta}) \propto 1$ and $\pi(\sigma) \propto \sigma^{-1}$. We also conduct the maximum likelihood estimation GHP1980 for AR using the `arima' function in the R package stats, which leads to almost identical performance and is thus not reported here. As such, we use AR and BayesAR exchangeably throughout this paper. We use AR($p$)-GARCH(1,1) when implementing GARCH, i.e.,
We fit the model using the R package fGarch, where all parameters are estimated by quasi-maximum likelihood bollerslev1992quasi. In addition, we implement the Quantile Autoregression (QAR) method proposed by KX2006 using the R package quantreg, to compare its finite sample performance with BayesMAR.
We assess estimates of $(\beta_0, \beta_1, \beta_2)$ by each method based on mean squared error (MSE). For a generic parameter $\theta$, letting the estimate be $\hat{\theta}_i$ in the $i$th simulation and $\bar{\theta} = \frac{1}{100} \sum^{100}_{i=1} \hat{{\theta}}_i$, then the MSE and its standard error are $$ \text{MSE} = \frac{1}{100} \sum_{i = 1}^{100} (\hat{\theta}_i - \theta)^2, \enskip \text{SE}_{\text{MSE}} = \frac{1}{10}\sqrt{ \text{sample variance of } \{(\hat{\theta}_{i} - {\theta})^2\}_{i = 1}^{100}}.$$
We first provide the true order $p = 2$ to all models and compare their performances. Table (ref) reports the MSE of all methods. We can see that all methods benefit from a correctly specified error distribution: AR and GARCH have the smallest MSEs under Gaussian error, while BayesMAR and QAR have smaller MSEs when data are generated from Laplace distributions. However, BayesMAR appears to suffer less than AR from model misspecification; for example, the increase of MSE of $\beta_2$ under Gaussian error from AR to BayesMAR is 0.13, which is within two standard errors, while AR doubles the MSE of $\beta_2$, and so is beyond three standard errors of BayesMAR under Laplace error. It is reassuring that BayesMAR gives either the same or better MSEs than QAR in all cases, although all differences are within one standard error. This finite sample performance is consistent with the findings in GCC2011 when comparing sampling-based Bayesian approaches with optimization-based counterparts for extreme quantile levels.
We next investigate the selection of the unknown order $p$ in BayesMAR using the BIC approach described in Section (ref). We provide a large upper bound $K = 20$ for the order $p$. Figure (ref) plots the distribution of $p$ in both scenarios and suggest that the selected orders almost always concentrate around the oracle value $p = 2$, even when the model is misspecified under Gaussian error. The overall accuracy across all simulations to select $p = 2$ using MAP is $95\%$ for normal errors and $98\%$ for Laplace errors.
In this section, we compare the predictive performances of the proposed methods to that of mean-based methods using various real-world data. We use three economic series from Federal Reserve Economic Data: the quarterly data Producer Price Index for all commodities FRED:p, 3-Month Treasury Bill: secondary market Rate FRED:r, and Unemployment Rate FRED:u, coded as PPI, TBR, and UR, respectively. Each time, the series ranges from 1968Q3 to 2018Q2, containing 200 observations. The unemployment rate data is seasonally adjusted by the method of seasonal-trend decomposition using Loess CCMT1990. The three time series $(y_t)$ and the lagged data of order one $(y_t - y_{t-1})$ are plotted in Figure (ref). These three data sets have distinct patterns: the lagged PPI tends to be more stable before the crisis in the year of 2008, in contrast to the substantial fluctuation since then; the lagged TBR has more dramatic changes in earlier periods than in later periods; and the lagged UR appears to contain several extreme values, while a periodic pattern may still persist even after seasonal adjustment. These complex characteristics of the data enable a comparison between model-based methods when there is no guarantee for model assumptions to hold.
In addition to AR and GARCH, we implement selected methods proposed in research on the dynamic linear model. In particular, we consider the time-varying vector AR (TV-VAR) model proposed by NW2013, which links time-varying parameters to latent threshold processes and achieves state-of-the-art predictive performance in selected applications. A vector with a length of three that stacks the three time series $\bm{y}_t^{(3)}$, expands the dynamic linear model by utilizing a latent threshold vector $\bm{d} = (d_1, ..., d_k)$ with $d_i \geq 0$ for $i = 1,\ldots, k $ :
where $\bm{c}_t$ is the $3 \times 1$ vector of time-varying intercepts, $\bm{B}_{jt}$ is the $3 \times 3$ matrix of time-varying coefficients at lag $j$, $\bm{b}_t$ stacks $\bm{c}_t$ and $\bm{B}_{jt}$ by rows and by order $j$ from 1 to $p$, $k = 3(1 + 3p)$, and $\bm{\beta}_{t} = (\beta_{1t}, ..., \beta_{kt})$ is a latent time-varying parameter vector whose dynamic updates are specified by $(\bm{\mu}_{\bm{\beta}}, \bm{\Phi}_{\bm{\beta}}, \bm{V}_{\bm{\beta}})$ via a standard VAR model. The latent threshold vector $\bm{d}$ shrinks the time-varying coefficients $\bm{\beta}_t$ to zero, leading to dynamic sparsity and improved prediction. We implement two versions of TV-VAR: TV-VAR without latent threshold (NT) and TV-VAR with latent threshold (LT). Both methods are applicable to multivariate time series, and we adapt the OxMetrics code provided in NW2013 to the vector observations consisting of the three time series.
For BayesMAR, we use BayesMAR-BMA to denote the Bayesian model averaging strategy and BayesMAR-MAP when the MAP estimate of $p$ is used, and we adopt the same convention for AR, i.e., BayesAR-BMA and BayesAR-MAP, which are implemented similarly to BayesMAR but with the Gaussian likelihood. GARCH(1,1) uses the same autoregressive order $p$ as chosen in BayesAR-MAP; although orders other than (1, 1) can be used, we observe that the GARCH(1,1) consistently led to the best predictive performances of GARCH in our real data application. Since there are no immediately available model selection methods for NT or LT, we run the order from 1 to 5, and then choose the optimal results to favor NT and LT under each criterion. We observe that LT considerably outperformed NT in all scenarios, and thus only present the results for LT.
We use recursive out-of-sample forecasting to assess each method after $t_0$=2008Q3. In particular, we fit each model with data up to quarter $t$ and conduct an $h$-step-ahead prediction for $h = 1, 2, 3, 4$. Then we move one period ahead and repeat the same procedure until we reach $t = T$. All methods are applied to the lagged data of order one to remove local trends and forecast the changes that lead to predictions of $y_t$.
We calculate the root MSE (RMSE) and mean absolute error (MAE) from $t_0$ to $T$ to compare performance of each method, i.e.,
We additionally calculate the the relative change using BayesMAR-BMA as the reference to ease comparison, that is,
Table (ref) reports the RMSEs of all methods for the three data sets. The results indicate that BayesMAR outperforms alternatives in nearly all scenarios. In particular, the RMSEs of BayesMAR-BMA is uniformly smaller than its AR counterpart. Compared to BayesMAR-BMA, the methods of AR, GARCH, and LT increase the RMSE by 10% to 89.2% for TBR, and 3.0% to 35.6% for UR, with two exceptions for the UR when compared to GARCH at three- and four-steps-ahead predictions where GARCH slightly reduces the RMSE by 0.3% and 2.8%, respectively. Although the gap in the prediction for PPI is smaller (between 0.8% to 8.9%), all competing methods have a larger RMSE than BayesMAR-BMA. GARCH accounts for comprehensive variance structures and LT dynamically updates regression coefficients, which are arguably powerful methods with considerable complexity; it is remarkable that the proposed simple BayesMAR leads to favorable and often superior performance using real-world data.
Table (ref) reports the MAEs of all methods, suggesting similar observations as those made from Table (ref). The proposed BayesMAR methods give the smallest MAEs in nearly all scenarios, with one exception for the UR when compared to GARCH.
The RMSEs of BayesMAR-BMA and BayesMAR-MAP are close to each other with at most 1.1% (PPI), -0.6% (TBR), 2.3% (UR) relative differences, and similar observations hold for MAEs. This is partly caused by highly concentrated weights of model orders when using BIC for order selection. In particular, we find BayesMAR-BMA puts most weight on the order selected in BayesMAR-MAP, leading to minimal differences between the two variants of BayesMAR. Comparing BayesAR-BMA and BayesAR-MAP leads to the same conclusion for AR.
Figure (ref) compares up to four-steps-ahead absolute predictive errors at each $t$ from $t_0$ to $T$ given by the three static models, MAR, AR, and GARCH, which provides insights into understanding the performance of BayesMAR. We choose AR and GARCH and implement MAR and AR using MAP such that all methods in the figure build on similar model structures. We can see that AR yields a large prediction error at the beginning (PPI), which is substantially reduced by GARCH, which incorporates heterogeneous variance structures. It is reassuring that the proposed BayesMAR method, which uses a simple constant variance structure, achieves the same predictive gain (for PPI) or even further improvements (for TBR and UR). At later time points when the time series stabilizes without usually large deviations, BayesMAR tends to perform similarly to AR and GARCH. Since BayesMAR is a parametric model bearing the same model structure as AR, these comparisons suggest that further performance gains may be possible by following the rich literature that extends AR to more advanced models such as GARCH and dynamic models, by simply altering the error assumption from Gaussian to Laplace.
We plot the 95% credible intervals of AR and MAR in Figure (ref), both using the BMA version. For both AR and MAR, four-steps-ahead predictive intervals appear wider than one-step-ahead predictive intervals. This makes sense, as more uncertainty propagates, and this is particularly the case for the UR data. MAR is more robust than AR at the early stage of the PPI time series. Since an informative comparison between MAR and AR through this visualization seems difficult, we next turn to numerical comparisons.
\color{black}
The proposed method provides density forecasts in addition to point forecasts; see Section (ref). We evaluate the probabilistic forecasts produced by the proposed method and AR through the continuous ranked probability score (CRPS) TFA2007, which is implemented in the R package scoringRules. We can see that the proposed MAR model produces a considerably smaller CRPS than AR for the TBR and UR data, while being close to AR for PPI. This suggests that the advantage of MAR may extend to probabilistic forecasts beyond point forecasts (see Table (ref)).
\color{black}
This article proposes a Bayesian median autoregressive (BayesMAR) model for robust time series forecasting. The proposed method has close connections with time-varying quantile regression. BayesMAR adopts a parametric model bearing the same structure as AR models by altering the Gaussian error to Laplace, leading to a simple, robust, and interpretable modeling strategy for time series forecasting with principled uncertainty quantification through Bayesian model averaging and Bayesian model selection. Real data applications using U.S. macroeconomic data show that BayesMAR leads to favorable and often superior predictive performances compared to the selected state-of-the-art mean-based alternatives under various loss functions that encompass both point and probabilistic forecasts.
BayesMAR enjoys practical benefits as a technical tool to introduce robustness and improve predictions. In addition, since the autoregressive structure in BayesMAR resembles the widely used AR model, BayesMAR can be used to complement a rich class of methods that build on the AR model. The AR model is arguably one of the most popular methods in time series, serving as the building block for other models such as GARCH and TV-VAR in research on the dynamic linear model. The proposed MAR model shows the potential for further performance gains by following the rich literature that extends AR to more advanced models with Laplace error assumptions rather than Gaussian ones.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
This work was partially supported by the grant DMS-2015569 from the National Science Foundation.
The R code to implement the proposed methods is publicly available at \url{https://github.com/xylimeng/BayesMAR}.