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.
45,070 characters · 13 sections · 69 citation commands
{.26in} \thispagestyle{empty}
Since their inception by breid1991maximum, Mixed causal-noncausal Auto-Regressive models (MAR) have gained significant popularity due to their ability to capture nonlinear dynamics in time series. Specifically, unlike conventional tests that detect explosive behaviors through unit roots (phillips2011explosive, phillips2015testing), MARs treat the variable under investigation as stationary and interpret bubble patterns as an intrinsic part of its dynamics, see, among others, fries2019mixed, gourieroux2013explosive, cavaliere2020bootstrapping and hecq2021forecasting. Moreover, MAR models are also an important tool in macroeconomics due to their ability to model expectations and their impact on the variable under investigation, see lanne2011noncausal. Finally, they are also useful for testing nonfundamentalness, see alessi2011non and lanne2013noncausal.
Over time, both parametric and semiparametric estimators have been developed for MARs. An example of a semiparametric estimator is the Generalized Covariance (Gcov), which aims to minimize a portmanteau-type objective function involving the autocovariances of linear and nonlinear transformations of model errors (see gourieroux2017noncausal and gourieroux2023generalized). Another semiparametric estimator is the one proposed by hecq2022spectral, which exploits information from higher-order cumulants by combining the spectrum and the bispectrum in a minimum distance estimation. In the parametric framework, the maximum likelihood estimator for univariate models was introduced by breid1991maximum and subsequently updated by lanne2011noncausal. For multivariate models, the maximum likelihood estimator was introduced by lanne2013noncausal and davis2020noncausal.
An alternative approach to estimating MARs is the Bayesian method introduced by lanne2012bayesian for univariate models, which is based on the Markov Chain Monte Carlo (MCMC) technique. The multivariate extension, discussed in lanne2016noncausal, takes advantage of the Metropolis-within-Gibbs sampler.
This paper aims to extend and adapt the Bayesian estimation method based on the sequential Monte Carlo algorithm (SMC), as illustrated in herbst2014sequential and bognanni2018sequential, to MARs. SMC offers an attractive alternative to MCMC for different reasons. First, this method is highly adaptable and can be easily implemented under different assumptions about the error term, as long as the posterior kernel can be evaluated point-wise. Second, it allows for fast initialization with random draws from the prior, unlike MCMC, which often requires a slow mode search. Third, SMC generates the Marginal Data Density (MDD) of the model as a byproduct, simplifying the calculation of Bayes factors without additional computational effort. Finally, SMC is parallelizable (see durham2014adaptive), which not only significantly reduces the estimation time but also helps with objective functions containing local minima, a common issue in MARs. Specifically, bec2020mixed highlights local minima in the objective function of the maximum likelihood estimator; hecq2022spectral in the spectral estimation function; cubadda2024optimization in the GCov objective function; lanne2016noncausal notes non-elliptical shapes, such as skewness and multimodality, in marginal posterior distributions.
The existing literature on noncausal models has explored Bayesian estimation only under the assumption that the error term follows a Student-$t$ distribution; see, among others, lanne2012bayesian, lanne2016noncausal, nyberg2014forecasting and moussa2023identifying. However, this assumption may be too restrictive in some cases, and alternative distributions may provide a better fit. For example, as shown in gourieroux2017local, the Cauchy distribution can more accurately capture extreme events, such as market crashes, while the Skewed $t$ distribution is better suited to model asymmetric patterns in financial time series, as discussed in proietti2023peaks. Since the distribution of the error term is usually unknown to the researcher, we propose to exploit the adaptability of SMC to distributions and introduce an identification methodology, based on the MDD and the Bayesian Information Criterion (BIC), that not only determines the causal and noncausal polynomial orders but also identifies the error term distribution that best fits. To the best of our knowledge, none of the already existing methodologies allows for simultaneously performing both model selection and identification of the error distribution.
The rest of the paper is organized as follows. Section (ref) reviews the univariate and multivariate MAR. Section (ref) introduces the SMC for MARs. In Section (ref) the ability of SMC for estimation and model identification is investigated by Monte Carlo exercises. Section (ref) reports the empirical application. Finally, Section (ref) draws some conclusions. Further results and derivations are reported in the Appendices.
Univariate MAR processes, also denoted as MAR($r,s$), have the following representation:
where $\psi(\mbox{L}) = \left(1 - \psi_1 \mbox{L} - \dots - \psi_r \mbox{L}^r\right)$ is the causal polynomial of order $r$; $\phi(\mbox{L}^{-1}) = \left(1 - \phi_1 \mbox{L}^{-1} - \dots - \phi_s \mbox{L}^{-s}\right)$ is the noncausal polynomial of order \(s\); and $p$ is defined as $p = r + s$. Moreover, both polynomials have roots outside the unit circle:
Due to the noncausal polynomial, the process in (ref) is capable of capturing nonlinear dynamics, including local trends (bubbles) and conditional heteroscedasticity (see hencic2015noncausal, hecq2016identification and gourieroux2018misspecification). Finally, to successfully identify these models, it is well known that $u_t$ must be independent and identically distributed ($i.i.d.$) and non-Gaussian, see breid1991maximum and lanne2011noncausal.\\ The series $y_t$ in (ref) admits a two-sided moving average representation:
where $\theta_0$ is equal to 1. If the process is purely causal ($s=0$) or purely noncausal ($r=0$), then $\theta_j = 0$ for $j > 0$ or for $j < 0$, see gourieroux2015uniqueness.
Vector MAR processes, denoted as VMAR($r$,$s$) were introduced by lanne2013noncausal and are given by:
where: $\mathbf{y}_t = \left(y_{1,t}, \dots , y_{n,t}\right)^{\prime}$ and $\mathbf{u}_t = \left(u_{1,t}, \dots , u_{n,t}\right)^{\prime}$ are two $n \times 1$ vectors; $\Phi(\mbox{L}^{-1}) = \mbox{I}_n - \Phi_1\mbox{L}^{-1} - \dots - \Phi_s \mbox{L}^{-s} $ is a $n \times n$ noncausal polynomial of order $s$; and $\Psi(\mbox{L}) = \mbox{I}_n - \Psi_1\mbox{L} - \dots - \Psi_r\mbox{L}^{r}$ is a $n \times n$ causal polynomial of order $r$. Finally, $p$ is defined as $p = r + s$. As in the univariate framework, both the causal and noncausal polynomials have their roots outside the unit circle:
and the identification of the noncausal component requires two conditions: the error term must be $i.i.d.$ and non-Gaussian, see gourieroux2017noncausal and cubadda2023detecting.\\ Similar to Section (ref), the process in (ref) exhibits a two-sided moving average specification:
where $\Theta_0$ is equal to a $n \times n$ identity matrix. If the process is purely causal ($s=0$), or purely noncausal ($r=0$), then $\Theta_j$ is a null matrix for $j>0$ or for $j<0$; see, among others, lanne2013noncausal.
As shown in gourieroux2017noncausal and davis2020noncausal, a potential alternative mixed model specification is given by:
where the autoregressive polynomial $|\Pi(z)| = |\mbox{I}_n - \sum_{j=1}^{p} \Pi_j z^j|$ can have roots both inside and outside the unit circle. Specifically, the process defined in (ref) is considered purely causal if all roots lie outside the unit circle, purely noncausal if all roots lie inside the unit circle, and mixed causal and non-causal if the roots lie both inside and outside the unit circle. As highlighted in lanne2013noncausal, specifications (ref) and (ref) do not always overlap, and there are cases where the specification (ref) cannot be embedded in (ref), and vice versa. The question of whether there is a feasible specification that incorporates all or the majority of noncausal VAR processes is intriguing, but it is beyond the scope of this paper. Here, we focus exclusively on the VMAR($r,s$) specification. This is because, in addition to encompassing the univariate specification in (ref) as a special case, it aligns with the representation used in lanne2016noncausal. Finally, the VMAR($r,s$) representation allows for a neater interpretation of the results by explicitly separating the causal and noncausal components.
The Bayesian analysis of mixed models has been examined by lanne2012bayesian for univariate models and by lanne2016noncausal for multivariate models. We focus on the multivariate framework, as the process in (ref) is a special case of the process in (ref) when $n=1$; for the univariate example, see Appendix (ref).
Let $\mathbf{y}_t$ be generated by (ref), and define $\mathbf{y}=\left[ \mathbf{y}_{1},\dots,\mathbf{y}_{T}\right] ^{\prime}$, where $T$ denotes the sample size. We also define $\mathbf{\Theta}_1 = \left[\text{vec}(\Psi_1)^{\prime}, \dots, \text{vec}(\Psi_s)^{\prime}, \text{vec}(\Phi_1)^{\prime}, \dots, \text{vec}(\Phi_r)^{\prime}\right] ^{\prime}$) and $\mathbf{\Theta}_2$ a column vector that contains all the parameters that depend on the distribution of $\mathbf{u}_t$ (e.g., degrees of freedom, the vector half operator of the scale matrix, etc.). By defining the $k$-vector $\mathbf{\Theta} = [\mathbf{\Theta}_1^{\prime}, \mathbf{\Theta}_2^{\prime}]^{\prime}$, Bayes' theorem computes posterior probabilities by combining prior beliefs with observed data:
where $p(\mathbf{\Theta}|\mathbf{y})$ denotes the posterior probability, $p(\mathbf{y}|\mathbf{\Theta})$ represents the likelihood function, $p(\mathbf{\Theta})$ indicates the prior probability distribution of the parameters, and $p(\mathbf{y})$ represents the MDD, or marginal likelihood, which is important for model selection and prediction.
To compute $p(\mathbf{\Theta}|\mathbf{y})$, lanne2012bayesian and lanne2016noncausal employed a tailored Gibbs sampling algorithm. For MDD, p($\mathbf{y}$), lanne2012bayesian used the method introduced by geweke2005bayesian, while lanne2016noncausal adopted a technique similar to MitISEM. In this paper, we propose to compute the posterior distribution and MDD for both the univariate and multivariate frameworks using a different approach, namely the SMC algorithm, as described in Section (ref). As mentioned above, SMC enables the calculation of MDD as a simple by-product, without the need for additional processing. It also mitigates the risk of local minima, a common problem in the posterior distributions of MARs (see lanne2016noncausal and bognanni2018sequential). Finally, SMC is highly adaptable to a wide range of error term distributions, provided that the posterior kernel can be evaluated point-wise.
The SMC algorithm proposed by herbst2014sequential is based on Importance Sampling (IS). The IS approximates the target density $f(\cdot)$ using a different and easily to sample density $g(\cdot)$ called the candidate density:
where $h:\mathbf{\Theta}\rightarrow\mathbb{R}^k$, $\pi(\mathbf{\Theta})=p(\mathbf{\Theta}|\mathbf{y})$, $f(\mathbf{\Theta})=p(\mathbf{y}|\mathbf{\Theta})p(\mathbf{\Theta})$, $\mbox{Z}=p(\mathbf{y})$, and $w(\mathbf{\Theta})=f(\mathbf{\Theta})/{g(\mathbf{\Theta})}$ are the importance sample weights. If $\mathbf{\Theta}^{i^{\text{i.i.d.}}} \sim g(\theta)$, for $ i = 1, \dots, P $, then, under suitable regularity conditions (geweke1989bayesian), the Monte Carlo estimate is given by:
As $ P $ goes to infinity, it converges almost surely (a.s.) to $\mbox{E}_\pi[h(\mathbf{\Theta})]$, see bognanni2018sequential. Each $\tilde{W}_i$, referred to as a (normalized) importance weight (IW), is assigned to an associated parameter $\mathbf{\Theta}_i$. We refer to the pair $ (\mathbf{\Theta}_i, \tilde{W}_i) $ as a particle. The collection of particles $ \{ (\mathbf{\Theta}_i, \tilde{W}_i) \}_{i=1}^{P} $ forms a discrete distribution that provides an approximation of $ \pi(\mathbf{\Theta})$. The accuracy of the approximation is determined by the divergence between $ g(\cdot) $ and $ f(\cdot) $. However, constructing a suitable proposal distribution $g(\cdot)$ is challenging, especially when little is known about the shape of $f(\cdot)$. SMC offers a clever solution to this problem by iteratively constructing particle approximations to a sequence of distributions. Starting from the prior distribution, a sequence of bridge distributions that gradually incorporates more likelihood information is constructed. This procedure is repeated until the full likelihood has been assimilated in the posterior distribution, more precisely:
where $M$ are the number of stages that depend on the complexity of the problem at hand, and the values of the tempering parameter $\rho_m$ are given by an increasing sequence of values such that $\rho_1=0$ and $\rho_{M}=1$. Notice that $\rho_m$ is constructed so that more and more likelihood information is incorporated in the posterior as $m$ increases. As in herbst2014sequential and bognanni2018sequential, we define it as:
where the hyper-parameter $\lambda$, which is greater than 0, determines the rate at which the likelihood information is incorporated into the sampler. In particular, when $\lambda = 1$ the schedule is linear, and each stage contributes equally. When $\lambda > 1$, smaller increments in likelihood are prioritized in the initial stages, with larger increments occurring in the later stages. Conversely, when $\lambda < 1$ larger increments in the early stages are prioritized, with smaller increments in later stages.
A key drawback of IS is weight degeneracy, which occurs when the target distribution is sharp relative to the proposal, resulting in most IWs approaching zero. This leads to an insufficient number of representative samples and is closely related to the curse of dimensionality. Specifically, the issue arises from the extreme values the likelihood takes on in high-dimensional spaces. SMC mitigates this challenge through a tempering approach that progressively incorporates the likelihood.
Finally, following bognanni2018sequential the MDD take the following form:
where $P$ are the IS draws, and $\tilde{w}_n^i$ are the the incremental and normalized weights, obtained by: $$ \tilde{w}_m^i = [ p(\mathbf{y}|\mathbf{\Theta}_{m-1}^i)]^{\rho_m-\rho_{m-1}}. $$ See Appendix (ref) for a detailed description of the SMC algorithm and MDD calculation.
The performance of the SMC algorithm for VMAR($r,s$) models is evaluated using three Monte Carlo exercises. In each exercise, we set $r=s=1$, $n=2$ and the error term is generated from one of the following multivariate distributions: Cauchy, Student-$t$, and Skewed-$t$. For each Data Generating Process (DGP), we estimate several VMAR($r,s$) models using these three error distributions and various combinations of polynomial orders $r$ and $s$ to examine how often the algorithm correctly identifies them. Furthermore, the accuracy of the algorithm in estimating true population parameters is evaluated under the assumption that both the error distribution and the polynomial orders are known.
Sections (ref) and (ref) describe the priors and likelihood functions used in the Monte Carlo exercises, while Section (ref) outlines the Monte Carlo settings and presents the results.
As discussed in Section (ref), $\mathbf{\Theta}=[\mathbf{\Theta}_1^{\prime}, \mathbf{\Theta}_2^{\prime}] ^{\prime}$ indicates the prior probability of the VMAR($r,s$) parameters. Specifically, $\mathbf{\Theta}_1$, which collects the elements of the VMAR($r,s$) coefficient matrices, remains unchanged regardless of the error term distribution, whereas $\mathbf{\Theta}_2$ varies with the density distribution of $\mathbf{u}_t$ since it includes the vectorization of the scale matrix $\Sigma$, which is shared by the Student-$t$, Cauchy, and Skewed-$t$ distributions; the degrees of freedom ($\nu$) for the Student-$t$ and Skewed-$t$ distributions; and the vector containing the skewness parameters ($\mbox{\boldmath $\alpha$}=\left[\alpha_1, \dots , \alpha_n\right]^{\prime}$) for the Skewed-$t$ distribution.
For $\mathbf{\Psi}$ and $\mathbf{\Phi}$ we use Multivariate Normal ($ \mathcal{MN}(\cdot, \cdot)$) priors:
where $\mbox{\boldmath ${\cal I}$}(\mathbf{\Psi})$ and $\mbox{\boldmath ${\cal I}$}(\mathbf{\Phi})$ are indicator functions equal to unity in the stationary region. As in lanne2016noncausal, the parameters $\mbox{\boldmath $\gamma$}$ and $\mbox{\boldmath $\delta$}$ decrease as the number of leads and lags increases, implying the use of Minnesota-type priors. We set $\mbox{\boldmath $\gamma$} = 2/i$ for $i = 1, \dots, r$ and $\mbox{\boldmath $\delta$} = 2/q$ for $q = 1, \dots, s$. Thus, we consider more diffuse priors compared to lanne2016noncausal, who set $\mbox{\boldmath $\gamma$} = 1/i$ and $\mbox{\boldmath $\delta$} = 1/q$.
The prior distribution for the scale matrix, which is common to all the error distributions, is an Inverse Wishart ($\mbox{\boldmath ${\cal W}$}^{-1}(\cdot, \cdot)$):
where a diffuse prior is adopted by setting $\Psi_0 =5\mbox{I}_n$ for the scale matrix and $\tilde{\nu}_{\Sigma} = 3$ for the degrees of freedom
Finally, we use Exponential prior ($\mathcal{E}(\cdot)$) for the degrees of freedom ($\nu$) and $\mathcal{MN}(\cdot, \cdot)$ for the skewness parameter $\alpha$:
where, we set $\nu_0=5$ and $\kappa=3$.
The second component of the Bayesian analysis described in Section (ref) is the likelihood function $p(\mathbf{y} \mid \mathbf{\Theta})$. We use the same (approximate) likelihood functions used in lanne2013noncausal and lanne2016noncausal, given by:
where $T$ is the sample size, and $p(\mathbf{u}_t \mid \mathbf{\Theta})$ varies according to the error term distribution considered. Specifically, for the Multivariate Student-$t$ we have:
where $\nu$ indicates the degrees of freedom of the error term, see lanne2013noncausal.
The Multivariate Skewed-$t$ distribution has the following likelihood function:
where $\alpha$ is an $n$-vector of shape parameters of the error term, $t_n(\mathbf{u}_t,\nu)$ is the probability density function of an $n$-dimensional Student-$t$ random variable with $\nu$ degrees of freedom, and $\mbox{T}_1[\cdot, \nu + n]$ is the cumulative distribution function of a scalar Student-$t$ distribution with $\nu + n$ degrees of freedom.
The Cauchy distribution has the following likelihood function:
that is a special case of (ref) when $\nu$ is fixed to 1.
azzalini2003distributions underlines that when $\alpha$ is a zero vector, then equations (ref) and (ref) coincide. Moreover, lanne2016noncausal considers the multivariate Student-$t$ only for cases where $\nu > 2$, which ensures a finite variance of the error term. We impose the same restriction on (ref) and (ref). Finally, in equation (ref), (ref) and (ref) it is common to lose the first $r$ and last $s$ observations as discussed in lanne2011noncausal, lanne2013noncausal.
As previously stated, we set $r = s = 1$, $n = 2$, and we assume that the error term follows three multivariate distributions: Cauchy, Student-$t$, and Skewed-$t$. We consider a sample size of $T = 150$ observations and $B = 200$ Monte Carlo replications. Relatively small values for $T$ and $B$ are chosen because increasing them would lead to excessive computational time for the Monte Carlo exercise.
For each DGP, we estimate twenty-one models in each replication: VMAR($1,0$), VMAR($0,1$), VMAR($2,0$), VMAR($1,1$), VMAR($0,2$), VMAR($2,1$), and VMAR($1,2$) across the Cauchy, Student-$t$, and Skewed-$t$ distributions. Finally, in the SMC algorithm, we set the number of particles $P = 10000$, the rate at which the likelihood information is incorporated into the sample $\lambda = 2$, and the number of stages $M = 100$. We select the best-fitting model using MDD and BIC. The results, summarized in Table (ref), show the good ability of both MDD and BIC to select the correct error term distribution, as well as the causal and noncausal orders. However, when such combinations of $P$, $M$, and $\lambda$ are considered, BIC outperforms MDD in detecting the correct model specification, especially when the error term follows a Skewed-$t$ distribution. In particular, under the Cauchy and Student-$t$ assumptions, the results show that the second most frequently identified model is the one with the correct error term distribution, but with overidentification of the causal order. In contrast, when Skewed-$t$ DGP is considered, the second most identified model not only shows a higher percentage than in the previous two cases, but is also associated with the wrong error term distribution (Student-$t$), despite correctly identifying the causal and noncausal orders. This suggests that due to the larger number of parameters associated with the skewness of the error term, more stages may be required to more closely track the evolution of the system. However, increasing $M$ becomes computationally impractical due to the significant increase in the Monte Carlo simulation time.
Table (ref) displays the performance of the SMC estimates when polynomial orders and the error distribution are assumed to be known. The table displays the true parameter, the Monte Carlo average of the parameter estimates and their variance, the average of the standard errors obtained in each replication, the average bias (BIAS):
and the average root mean squared error (RMSE):
$\widehat{\mathbf{\Theta}}^{\mathcal{D}}_{j,b}$ represents the estimated value, in simulation $b$ of the $j$-th element of the vector $\mathbf{\Theta}$, as defined in Section (ref), with dimension $N_\mathcal{D}$ changing according to the distribution of the error term, $\mathcal{D}$. Finally, $ \mathbf{\Theta}^{\mathcal{D}_i}_{j,0} $ denotes the true population value of the corresponding element $j$. The results show that when the error term follows a Cauchy distribution, the SMC algorithm performs well, with estimates of the parameters $\mathbf{\Psi}$, $\mathbf{\Phi}$, and $\Sigma $ being very close to the true values, as indicated by the small bias and the RMSE values. The variance of the estimates is also relatively small, suggesting a stable estimation between replications. For the Student-$t$ distribution, the parameter estimates show slightly larger deviations from the true values compared to the Cauchy case. This is reflected in the higher bias and the RMSE, but the algorithm still exhibits a good performance. Finally, for the Skewed-$t$ distribution, the SMC algorithm shows a slight negative bias for the scale and skewness parameters.
\
We illustrate the use of SMC in mixed models by considering a bivariate process involving the S$\&$P Europe 350 ESG Index and Brent crude oil prices. The S$\&$P Europe 350 ESG Index is a stock market index designed to measure the performance of companies in Europe that meet certain environmental, social, and governance (ESG) criteria. It is a variant of the broader S$\&$P Europe 350 Index, which tracks the largest and most liquid stocks in the European market. The goal here is to investigate the dynamic interactions between oil prices and the European ESG sector. The data, with monthly frequency, span from July 2017 to February 2024, comprising $T=116$ observations, and available at \url{https://www.spglobal.com/spdji/en/indices/sustainability/sp-europe-350-esg-index/#overview} and \url{https://fred.stlouisfed.org/series/DCOILBRENTEU}.
The time series are shown in Figure (ref)-(a). As both series exhibit stochastic trends, we detrend them using a third-order polynomial. This choice strikes a balance between flexibility and parsimony, avoiding the risk of overfitting while capturing the underlying trend. Additionally, the third-order polynomial preserves key features in the data, such as the bubble patterns. Figure (ref)-(b) displays the detrended series.
Next, using $P = 50000$, $M = 150$, and $\lambda = 2$, we estimate the series using the SMC algorithm, considering the following models: VMAR($1,0$), VMAR($0,1$), VMAR($2,0$), VMAR($1,1$), VMAR($0,2$), VMAR($2,1$), and VMAR($1,2$) under the Cauchy, Student-$t$, and Skewed-$t$ distributions. We select the error distribution and polynomial orders $r$ and $s$ that maximize MDD and minimize BIC (Table (ref)). Both methods suggest a VMAR$(1,1)$, with the Student-$t$ distribution providing the best fit for the error term. The estimated coefficients are shown in Table (ref).
The 0.43 coefficient in matrix $\mathbf{\Phi}_1$ shows that an expected future increase of S$\boldsymbol{\&}$P Europe 350 ESG Index tends to increase the Brent price. A possible interpretation of the results may lie in the fact that as companies improve their ESG profiles, they may invest in greener technologies or renewable energy sources, which can lead to an increase in the short-term demand for oil. In other words, companies might need oil during the transition phase toward sustainability, where they are still reliant on fossil fuels for energy. Thus, a positive shift in ESG perceptions can create a paradoxical increase in demand for Brent, as firms balance immediate operational needs with long-term sustainability goals. Furthermore, Table (ref) also suggests that an expected increase in Brent prices slightly increases the S$\&$P Europe 350 ESG Index, implying that higher oil prices may positively impact the stock performance of ESG-compliant companies in Europe.
The paper extends and adapts the Bayesian SMC algorithm to handle MAR processes. Key advantages of the SMC method include faster initialization, parallelized computations, mitigation of issues related to local minima, and automatic generation of Marginal Data Density (MDD) as a byproduct. Consequently, SMC can be easily applied to a wide range of error term distributions, provided that the MDD can be estimated point-wise. Moreover, since the distribution of the error term is typically unknown to researchers, we exploit the adaptability of SMC to propose a novel identification procedure. This procedure, based on MDD and the BIC, differs from previous studies in that it not only selects the causal and noncausal polynomial orders but also identifies the error term distribution that best fits the data. As a result, we explore the Bayesian estimation of such models for the first time, not only in the Student-$t$ framework. Monte Carlo experiments demonstrate that SMC performs well in terms of bias and RMSE, successfully identifying the causal and noncausal polynomial orders as well as the error term distribution.
In the empirical application, we applied the SMC algorithm to a 2-dimensional VMAR($r,s$) model, including the S$\&$P Europe 350 ESG index and Brent crude oil prices. The causal and noncausal polynomial orders of the bivariate process are detected as $r=s=1$. Furthermore, we find that the assumption of multivariate Student-$t$ distribution is the one that better fits the data compared to the Cauchy and Skewed-$t$ distributions. Finally, our results indicate that an expected increase in the S$\&$P Europe 350 ESG Index can lead to a higher demand for Brent crude oil, while an expected rise in Brent prices positively influences the S$\&$P Europe 350 ESG Index.