Detecting Common Bubbles in Multivariate Mixed Causal–Noncausal Models

Abstract

This paper proposes concepts and methods to investigate whether the bubble patterns observed in individual time series are common among them. Having established the conditions under which common bubbles are present within the class of mixed causal–noncausal vector autoregressive models, we suggest statistical tools to detect the common locally explosive dynamics in a Student t-distribution maximum likelihood framework. The performances of both likelihood ratio tests and information criteria were investigated in a Monte Carlo study. Finally, we evaluated the practical value of our approach via an empirical application on three commodity prices.

1. Introduction

Economic and financial time series may exhibit many distinctive characteristics, including the presence of serial correlation, stochastic or deterministic trends, seasonality, time varying volatility, and nonlinearities. However, when the focus of the analysis is on the relationships among various variables, it is commonly observed that one or more of these features detected in individual series are common to several variables. We talk about common features when such features are annihilated by some suitable linear combinations of variables (Engle and Kozicki 1993). The most well-known example is probably cointegration, which is the presence of common stochastic trends (Engle and Granger 1987). Other forms of comovements have also been studied, giving rise to developments around the notions of common cyclical features (Vahid and Engle 1993), common deterministic seasonality (Engle and Hylleberg 1996), common volatility (Engle and Susmel 1993), cobreaking (Hendry and Massmann 2007), etc. Recognizing these common feature structures presents numerous advantages from an economic perspective (e.g., the whole literature on the existence of long-run relationships), but there are also several benefits for statistical modeling. Indeed, imposing the commonalities in estimation reduces the number of parameters, thus potentially leading to efficiency gains in statistical inference and improvements in forecasts accuracy (Issler and Vahid 2001). Moreover, the presence of common dynamics can be used for analyzing and forecasting large dimensional systems (Cubadda and Hecq 2011, 2022a; Bernardini and Cubadda 2015).

Building on the common feature approach, we propose the detection of common bubbles in stationary time series. Intuitively, the idea is to detect bubble patterns in univariate time series and then to investigate whether those bubbles are common to a set of assets. In the affirmative case, a portfolio composed of those series would not have such a nonlinear local explosive characteristic. There are several ways to capture bubbles in the data. We rely on mixed causal–noncausal models (denoted as MAR$(r,s)$ hereafter), namely autoregressive time series that depend on both r lags and s leads. Indeed, there has been recent interest in the properties of noncausal processes associated with a blooming of applications on commodity prices, inflation or cryptocurrency series, and the developments around the notion of nonfundamental shocks, see, i.a., Hecq and Voisin (2022) and the references therein. We choose to consider mixed causal and noncausal, models as they might also be used for forecasting. This is not necessarily the case with other approaches aimed at identifying bubble phases.

A first attempt into this direction was made by Cubadda et al. (2019), who extended the canonical correlation framework of Vahid and Engle (1993) from purely causal vector autoregressive models (namely the traditional serial correlation common feature approach within a VAR) to purely noncausal VARs (a VAR with leads only). They showed that different forms of commonalities can emerge when we also look at VARs in reverse time. However, their approach, being based on either canonical correlation analysis or the general method of moments, do not work for mixed models where non-Gaussianity of the error terms is required for identification, see e.g., Lanne and Saikkonen (2013).

In this paper, we extend their work and we propose a Student t-distribution maximum likelihood (ML henceforth) framework to compare the multivariate mixed causal–noncausal model with r lags and s leads (VMAR$(r,s)$ hereafter) with a restricted version where a reduced rank structure is imposed on the lead polynomial matrix. This is our notion of common bubbles, which is equivalent to requiring that there exist linear combinations of variables exhibiting bubbles that no longer possess the bubble feature. See for instance Cubadda and Hecq (2022b) for a recent survey on reduced rank techniques for common feature analysis. We consider both likelihood ratio tests and information criteria for our purposes.

Given that explosive roots and noncausal dynamics in VARs are intimately related (see, e.g., Gourieroux and Jasiak (2017) and the references therein), our approach has a similar spirit as the one of Engsted and Nielsen (2012), who proposed a test for the hypothesis that stock prices and dividends possess a common explosive root. Possible comparative merits of our methodology are that it does not require the prior knowledge of the value of the common explosive root and that as we work with stationary VMAR models, only standard asymptotic theory applies.

The rest of this paper is organized as follows: in Section 2 we establish the notations for multivariate mixed causal and noncausal models. Contrary to the univariate case, two distinct multivariate multiplicative representations lead to the same additive form of the VMAR$(r,s)$. Consequently, such alternative representations have the same likelihood but with different lag–lead polynomial matrices. We advocate the use of the multiplicative representation where the lead polynomial matrix is the first factor since the alternative representation does not allow for the easy unraveling of the presence of common bubbles. Within a Student t-distribution ML framework, we explain how to implement both likelihood ratio tests and information criteria to detect the existence of common bubbles. Section 3 investigates, using Monte Carlo simulations, the small sample properties of our strategy for bivariate and trivariate systems both under the null of common bubbles and the alternative of no rank reductions. Section 4 illustrates the practical value of our approach with an empirical analysis of three commodity prices. Section 5 concludes.

2. Multivariate Mixed Causal–Noncausal Models

Recall that a univariate MAR($r,s$) model is constructed as follows:

$(1-{\varphi }_{1}L-\cdots -{\varphi }_r{L}^r)(1-{\psi }_1{L}^{-1}-\cdots -{\psi }_s{L}^{-s})y_t=e_t,$

where $L^r$ is the lag operator such that $L^r{y}_t=y_{t-r}$ and $L^{-s}$ is the lead operator such that $L^{-s}{y}_t=y_{t+s}$. Since all the coefficients are scalars, the polynomial product is commutative and the representation

$(1-{\psi }_1{L}^{-1}-\cdots -{\psi }_s{L}^{-s})(1-{\varphi }_{1}L-\cdots -{\varphi }_r{L}^r)y_t=e_t.$

will yield the same model parameters as the previous one. The error term $e_t$ is assumed to be i.i.d. and non-Gaussian for identification purposes.

Let us now consider the case where $Y_t$ is an N-dimensional stationary process. For the sake of simplicity, we assume that deterministic elements are absent. Analogously to the univariate case, the VMAR$(r,s)$, is defined in its multiplicative forms as follows:

(1)$\mathsf{\Psi }\left(L^{-1}\right)\mathsf{\Phi }\left(L\right)Y_t={\epsilon }_t,$

(2)$\overline{\mathsf{\Phi }}\left(L\right)\overline{\mathsf{\Psi }}\left(L^{-1}\right)Y_t={\overline{\epsilon }}_t.$

where

$\begin{array}{cc}\hfill \mathsf{\Psi }\left(L^{-1}\right)\mathsf{\Phi }\left(L\right)& =(I_N-{\mathsf{\Psi }}_1{L}^{-1}-\cdots -{\mathsf{\Psi }}_s{L}^{-s})(I_N-{\mathsf{\Phi }}_1{L}^1-\cdots -{\mathsf{\Phi }}_s{L}^r),\hfill \\ \hfill \overline{\mathsf{\Phi }}\left(L\right)\overline{\mathsf{\Psi }}\left(L^{-1}\right)& =(I_N-{\overline{\mathsf{\Phi }}}_1{L}^1-\cdots -{\overline{\mathsf{\Phi }}}_s{L}^r)(I_N-{\overline{\mathsf{\Psi }}}_1{L}^{-1}-\cdots -{\overline{\mathsf{\Psi }}}_s{L}^{-s}).\hfill \end{array}$

Both models (1) and (2) are equivalent in the sense that they generate the same time series, but due to the noncommutativity property of the matrix product, they are two distinct representations of the same process. Specifically, the lag polynomial matrices $\mathsf{\Phi }\left(L\right)$ and $\overline{\mathsf{\Phi }}\left(L\right)$, although of the same order r, have unequal values of the coefficient matrices, and the same observation applies to the $s-$order lead polynomial matrices $\mathsf{\Psi }\left(L^{-1}\right)$ and $\overline{\mathsf{\Psi }}\left(L^{-1}\right)$ as well.

We assume that ${\epsilon }_t$ and ${\overline{\epsilon }}_t$ are i.i.d. and follow multivariate Student t-distribution with location zero. We could consider different distributions as long as they are non-Gaussian. This is indeed the condition that allows for distinguishing the genuine VMAR$(r,s)$ specification from the so-called pseudo causal and noncausal representations, see Lanne and Saikkonen (2013) for details.

We further assume that the roots of the determinant of each of the polynomial matrices $\mathsf{\Psi }\left(L^{-1}\right),\mathsf{\Phi }\left(L\right)$, $\overline{\mathsf{\Phi }}\left(L\right),\overline{\mathsf{\Psi }}\left(L^{-1}\right)$ are outside the unit circle to fulfill the stationarity condition. Furthermore, we will show later that the distribution of the errors ${\epsilon }_t$ and ${\overline{\epsilon }}_t$ have identical degrees of freedom $\lambda \in {\mathbb{R}}^{+}$ but different positive definite scale matrices, which are respectively denoted by $\Sigma$ and $\overline{\Sigma }$.

Let us respectively denote with $A\left(L\right)$ and $\overline{A}\left(L\right)$ the products of the lag and lead matrix polynomials of the two models (1) and (2):¹

$\begin{array}{cc}\hfill \mathsf{\Psi }\left(L^{-1}\right)\mathsf{\Phi }\left(L\right)& \equiv A\left(L\right)=\sum _{j=-s}^r{A}_j{L}^j\to A\left(L\right)Y_t={\epsilon }_t,\hfill \\ \hfill \overline{\mathsf{\Phi }}\left(L\right)\overline{\mathsf{\Psi }}\left(L^{-1}\right)& \equiv \overline{A}\left(L\right)=\sum _{j=-s}^r{\overline{A}}_j{L}^j\to \overline{A}\left(L\right)Y_t={\overline{\epsilon }}_t.\hfill \end{array}$

The general forms of the product of the lead and lag matrix polynomials for both the representation respectively read as follows:

$A\left(L\right)=\underset{A_0}{\underbrace{I+\sum _{i=1}^{\begin{array}{c}min\\ \{r,s\}\end{array}}{\mathsf{\Psi }}_i{\mathsf{\Phi }}_i}}-\sum _{i=1}^r\left(\underset{A_i}{\underbrace{{\mathsf{\Phi }}_i-\sum _{\begin{array}{c}\forall \{l,m\}:\\ l-m=i\end{array}}{\mathsf{\Psi }}_l{\mathsf{\Phi }}_m}}\right)L^i-\sum _{j=1}^s\left(\underset{A_j}{\underbrace{{\mathsf{\Psi }}_j-\sum _{\begin{array}{c}\forall \{l,m\}:\\ m-l=j\end{array}}{\mathsf{\Psi }}_l{\mathsf{\Phi }}_m}}\right)L^{-j},$

$\overline{A}\left(L\right)\equiv \underset{{\overline{A}}_0}{\underbrace{I+\sum _{i=1}^{\begin{array}{c}min\\ \{r,s\}\end{array}}{\overline{\mathsf{\Phi }}}_i{\overline{\mathsf{\Psi }}}_i}}-\sum _{i=1}^r\left(\underset{{\overline{A}}_i}{\underbrace{{\overline{\mathsf{\Phi }}}_i-\sum _{\begin{array}{c}\forall \{l,m\}:\\ m-l=i\end{array}}{\overline{\mathsf{\Phi }}}_m{\overline{\mathsf{\Psi }}}_l}}\right)L^i-\sum _{j=1}^s\left(\underset{{\overline{A}}_j}{\underbrace{{\overline{\mathsf{\Psi }}}_j-\sum _{\begin{array}{c}\forall \{l,m\}\\ m-l=j\end{array}}{\overline{\mathsf{\Phi }}}_m{\overline{\mathsf{\Psi }}}_l}}\right)L^{-j}$

with $1\le l\le s$ and $1\le m\le r$. Hence, both the multiplicative representations yield exactly the same additive form

(3)$\underset{\begin{array}{c}{A}_0^{-1}A\left(L\right)\\ =\\ {\overline{A}}_0^{-1}\overline{A}\left(L\right)\end{array}}{\underbrace{B\left(L\right)}}{Y}_t=\underset{\begin{array}{c}{A}_0^{-1}{\epsilon }_t\\ =\\ {\overline{A}}_0^{-1}{\overline{\epsilon }}_t\end{array}}{\underbrace{{\eta }_t}},$

where ${\eta }_t$ follows a multivariate Student t-distribution with degrees of freedom $\lambda$, as ${\epsilon }_t$ and ${\overline{\epsilon }}_t$ in representations (1) and (2), and with a scale matrix $\Omega =A_0^{-1}\Sigma {\left(A_0^{-1}\right)}^{{}^{\prime }}={\overline{A}}_0^{-1}\overline{\Sigma }{\left({\overline{A}}_0^{-1}\right)}^{{}^{\prime }}$. The lag polynomial in (3) is the following:

(4)$B\left(L\right)=I-\sum _{i=1}^r{B}_i{L}^i-\sum _{j=1}^s{B}_{-j}{L}^{-j}.$

An example of derivation of the polynomial matrix $B\left(L\right)$ for VMAR$(2,2)$ is given in Section 2.1.

In summary, contrary to the univariate case, a VMAR$(r,s)$ process has two distinct multiplicative representations. ML inference can indifferently be performed with each of the two representations (1) and (2). However, both representations will correspond to the same additive form of the model in Equation (3). This makes the interpretation of the lag and lead coefficient matrices in the multiplicative forms more intricate. Lanne and Saikkonen (2013) advocated for the use of one or the other representation depending on the analysis performed; one representation might be easier to employ for certain inquiries.

2.1. Common Bubbles in VMAR(r,s)

Having discussed the main properties of the unrestricted VMAR, we consider additional restrictions coming from a reduced rank structure in the lead polynomial matrix in order to model common bubbles. Notably, the noncausal component explains the growth phase of the bubble, whereas the causal component determines the burst phase (Gouriéroux and Zakoïan 2017). Hence, it is the multiplicative structure of the VMAR, combined with heavy tailed errors, that captures the nonlinearity of the bubbles as a whole. However, without the noncausal component, heavy tailed errors in a causal AR setting would not be able to reproduce locally explosive episodes. Although the focus in this paper is on common bubbles, our approach can be easily extended to investigate commonalities in the causal part or in both the lag and the lead components.

Definition 1.

An N-dimensional VMAR$(r,s)$ process displays common bubbles (CBs hereafter) if there exists a full-rank matrix δ of dimension $N\times k$, with $0<k<N$, such that, ${\delta }^{\prime }{B}_{-j}=0$ for $j=1,\cdots ,s$, where the coefficient matrix $B_{-j}$ is defined in (4). This implies that the coefficient matrices $B_{-j}$ can be factorized as $B_{-j}={\delta }_{\perp }{\beta }_j^{\prime }$ where ${\delta }_{\perp }$ is the $N\times (N-k)$ orthogonal complement of ${\delta }^{\prime }$ such that ${\delta }^{\prime }{\delta }_{\perp }=0$ and ${\beta }_j$ constitute a matrix with dimensions $N\times (N-k)$.

For the sake of simplicity, let us start the analysis from the case $r=s=2$. The coefficient matrices of the leads in the additive representation (3) with reduced rank restrictions are as follows:

$\begin{array}{cc}\hfill B_{-1}& =A_0^{-1}({\mathsf{\Psi }}_1-{\mathsf{\Psi }}_2{\mathsf{\Phi }}_1)={\overline{A}}_0^{-1}({\overline{\mathsf{\Psi }}}_1-{\overline{\mathsf{\Phi }}}_1{\overline{\mathsf{\Psi }}}_2)={\delta }_{\perp }{\beta }_1^{\prime },\hfill \\ \hfill B_{-2}& =A_0^{-1}{\mathsf{\Psi }}_2={\overline{A}}_0^{-1}{\overline{\mathsf{\Psi }}}_2={\delta }_{\perp }{\beta }_2^{\prime },\hfill \end{array}$

where the matrices $A_0$ and ${\overline{A}}_0$ are as follows:

$\begin{array}{cc}\hfill A_0& =(I_N+{\mathsf{\Psi }}_1{\mathsf{\Phi }}_1+{\mathsf{\Psi }}_2{\mathsf{\Phi }}_2)\hfill \\ \hfill {\overline{A}}_0& =(I_N+{\overline{\mathsf{\Phi }}}_1{\overline{\mathsf{\Psi }}}_1+{\overline{\mathsf{\Phi }}}_2{\overline{\mathsf{\Psi }}}_2).\hfill \end{array}$

When k CBs exist, the matrix ${\delta }^{\prime }$ annihilates the forward looking dynamics

${\delta }^{\prime }{B}_{-1}={\delta }^{\prime }{B}_{-2}=0.$

This implies for the second lead coefficient matrices that

${\delta }^{\prime }{B}_{-2}={\delta }^{\prime }{A}_0^{-1}{\mathsf{\Psi }}_2={\delta }^{\prime }{\overline{A}}_0^{-1}{\overline{\mathsf{\Psi }}}_2=0.$

Since ${\delta }^{\prime }{A}_0^{-1}$ (resp. ${\delta }^{\prime }{\overline{A}}_0^{-1}$) cannot be equal to zero, it follows that ${\delta }^{\prime }{A}_0^{-1}={\gamma }^{\prime }$ (resp. ${\delta }^{\prime }{\overline{A}}_0^{-1}={\overline{\gamma }}^{\prime }$), with $\gamma$ (resp. $\overline{\gamma }$) being a full-rank $N\times k$ dimensional matrix, and thus ${\gamma }^{\prime }{\mathsf{\Psi }}_2=0$ (resp. ${\overline{\gamma }}^{\prime }{\overline{\mathsf{\Psi }}}_2=0$). Hence, both ${\mathsf{\Psi }}_2$ and ${\overline{\mathsf{\Psi }}}_2$ must have rank $N-k$ but potentially different left null spaces (see also Cubadda et al. 2019).

By premultiplying the first lead coefficient matrix of representation (1) by ${\delta }^{\prime }$, we have

${\delta }^{\prime }{B}_{-1}={\gamma }^{\prime }({\mathsf{\Psi }}_1-{\mathsf{\Psi }}_2{\mathsf{\Phi }}_1)={\gamma }^{\prime }{\mathsf{\Psi }}_1=0,$

which implies that ${\mathsf{\Psi }}_1$ and ${\mathsf{\Psi }}_2$ must have the same left null space. Moreover, keeping in mind that ${\delta }^{\prime }{A}_0^{-1}={\gamma }^{\prime }$, we have

${\delta }^{\prime }={\gamma }^{\prime }{A}_0={\gamma }^{\prime }(I_N+{\mathsf{\Psi }}_1{\mathsf{\Phi }}_1+{\mathsf{\Psi }}_2{\mathsf{\Phi }}_2)={\gamma }^{\prime },$

which shows that ${\mathsf{\Psi }}_1$ and ${\mathsf{\Psi }}_2$ have the same left null space as $B_{-1}$ and $B_{-2}$, respectively.

By premultiplying the first lead coefficient matrix of the alternative representation (2) by ${\delta }^{\prime }$, we instead have

${\delta }^{\prime }{B}_{-1}={\overline{\gamma }}^{\prime }({\overline{\mathsf{\Psi }}}_1-{\overline{\mathsf{\Phi }}}_1{\overline{\mathsf{\Psi }}}_2)=0,$

which implies that the matrix ${\overline{\mathsf{\Psi }}}_1$ might not even have a reduced rank.

It is easy, although tedious, to see that the same conclusion holds for any VMAR$(r,s)$. We summarize these results in the following proposition.

Proposition 1.

In the presence of k CBs in an N-dimensional VMAR$(r,s)$ process, we have ${\delta }^{\prime }\mathsf{\Psi }\left(L^{-1}\right)={\delta }^{\prime }$ in (1), whereas ${\delta }^{\prime }\overline{\mathsf{\Psi }}\left(L^{-1}\right)\ne {\delta }^{\prime }$ in (2). Hence, the same linear combinations annihilate the lead coefficient matrix both in the additive representation (3) and in the multiplicative representation (1) but not in the alternative multiplicative representation (2). This implies that the coefficient matrices ${\mathsf{\Psi }}_j$ can be factorized as ${\mathsf{\Psi }}_j={\delta }_{\perp }{\mathsf{\Gamma }}_j^{\prime }$ for $j=1,...,s$.

2.2. Testing for Common Bubbles

In view of Proposition 1, a likelihood ratio test (LRT henceforth) for the presence of k CBs requires the comparison of the likelihood value of the unrestricted model (1) with the likelihood value of the restricted model:

(5)$(I_N-{\delta }_{\perp }{\mathsf{\Gamma }}_1^{\prime }{L}^{-1}-\cdots -{\delta }_{\perp }{\mathsf{\Gamma }}_s^{\prime }{L}^{-s})(I_N-{\mathsf{\Phi }}_{1}L-\cdots -{\mathsf{\Phi }}_r{L}^r)Y_t={\epsilon }_t,$

Since ${\delta }_{\perp }$ has dimension $N\times (N-k)$ with $0<k<N$, there are $N-1$ possible reduced-rank models to consider for all possible k. Furthermore, since the matrix ${\delta }_{\perp }$ can be normalized such that

(6)${\delta }_{\perp }^{\prime }=[I_{N-k},\omega ],$

it has only $k\times (N-k)$ free parameters in $\omega$.

A sample Y of T observations drawn from an N-dimensional VMAR($r,s$) process with i.i.d. t-distributed errors has location 0, a positive definite scale matrix $\Sigma$, and degrees of freedom; $\lambda \in {\mathbb{R}}^{+}$ has the following log-likelihood function

$\begin{array}{cc}\hfill f\left(Y|\Sigma ,\lambda \right)& =\left(T-(r+s)\right)\times ln\left(\frac{\mathsf{\Gamma }\left(\frac{\lambda +N}{2}\right)}{{\left(\lambda \pi \right)}^{N/2}\mathsf{\Gamma }\left(\frac{\lambda }{2}\right)}\right)-\frac{T-(r+s)}{2}\times ln\left(|\Sigma |\right)\hfill \\ & -\frac{\lambda +N}{2}\times \sum _{t=r+1}^{T-s}ln\left[1+\frac{1}{\lambda }\left({\epsilon }_t^T{\Sigma }^{-1}{\epsilon }_t\right)\right]\hfill \end{array}$

where $\mathsf{\Gamma }\left(x\right)={\int }_0^{\infty }{u}^{x-1}{e}^{-u}du.$

Without any commonality restrictions, ${\epsilon }_t$ is given either by (1) or (2), depending on the representation chosen for the estimation. For the estimation of the likelihood function with commonality restrictions, ${\epsilon }_t$ is given by (5), where normalization (6) is imposed for identification of matrix ${\delta }_{\perp }^{\prime }$.² Hence, imposing the restrictions within the Student’s t-ML estimation framework of the MVAR model is straightforward. The LRT is then constructed as follows

(7)$L{R}_{k|0}=2\ln \left(\frac{{\widehat{L}}_0}{{\widehat{L}}_k}\right),$

where ${\widehat{L}}_k$ and ${\widehat{L}}_0$ are, respectively, the likelihood values associated with the restricted model (5) and with the unrestricted one (1). Under the null of k CBs, (7) follows an asymptotic ${\chi }_d^2$ distribution with $d=k\left(N(s-1)+k\right)$ degrees of freedom.

One can also perform an LRT for the null hypothesis such that the $k=\overline{k}$ versus the alternative $k=$k with $1<\underline{k}<\overline{k}<N$. The associated LRT statistic is as follows:

$L{R}_{\overline{k}|\underline{k}}=2\ln \left(\frac{{\widehat{L}}_{\underline{k}}}{{\widehat{L}}_{\overline{k}}}\right),$

which is asymptotically distributed as a ${\chi }_g^2$ distribution with $g=(\overline{k}-\underline{k})\left(N(s-1)+\overline{k}+\underline{k}\right)$ degrees of freedom.

An alternative is to select the best specification according to the minimization of an information criterion as follows:

(8)$BI{C}_{{}_{\kappa }}=K\ln \left(T\right)-2\ln \left({\widehat{L}}_{{}_{\kappa }}\right),$

(9)$AI{C}_{{}_{\kappa }}=2K-2\ln \left({\widehat{L}}_{{}_{\kappa }}\right),$

with $K=r{N}^2+(N-k)(sN+k)$ being the number of coefficients estimated in a model with $\kappa$ CBs for $\kappa =0,1,...,N-1$.

3. Monte Carlo Analysis

We investigated the performance of our strategies using Monte Carlo simulations to detect the common bubbles in bivariate and trivariate VMAR(1,1) models. We considered two sample sizes, $T=500$ and 1000, and two different degrees of freedom of the error term with very leptokurtic distributions, namely $\lambda =3,1.5$, to respectively consider a finite and infinite variance case. We employed lead coefficient matrices with and without reduced rank to analyze the detection of the correct model under the null of common bubbles and under the alternative of no such comovements. The coefficients employed in the bivariate settings are displayed in Table 1.

The results, based on 3000 replications for each combination of parameters, are reported in Table 2.³ All entries are the frequency of correctly detected model. That is, under the null of a CB, we report the proportion of correctly detected CB, and under the alternative of no CB, we report the proportion of correctly rejected CB. We hence performed the test $H_0:rank\left(\mathsf{\Psi }\right)=1$ against the alternative that the rank is 2. The LRTs were performed at a 95% confidence level. The information criteria detected a CB when the IC of the restricted model was lower than the one of the unrestricted model.

We can notice that the frequency of type I errors of the LRT increases when the variance of the errors becomes infinite and that it does not significantly decrease when the sample size gets larger. With finite variance ($\lambda =3$), the LRT has an appropriate size of around 5.5%, and it increases to around 8.6% when the degrees of freedom of the error distribution reach 1.5. Under the alternative, the LRT has a power of at least 99.9% across all parameter combinations, implying that it almost never detects a CB when there are none.

Regarding the model selection using information criteria, results show that BIC outperforms AIC. Under the null of a CB, BIC selects the correct model specification in 98.9% of the cases with finite variance and a sample size of 500. The frequency increases to 99.3% when the sample size increases to 1000. AIC, on the other hand, selects the correct model in only 83.8% of the cases and does not increase with the sample size. The frequency of the correctly selected model decreases for both when in the infinite variance case, but more drastically for AIC, which decreases to around 78%. BIC still selects the correct model for 96.8% of the cases with a sample size $T=500$, and the frequency increases to 97.7% for $T=1000$. Under the alternative of no CB, however, both ICs correctly select the unrestricted specification in more than 99.4% of cases across all parameter combinations.

We now turn to the trivariate case. In the presence of a CB, the rank of the lead coefficient matrix can be either 1 or 2. We thus consider the two possible CB structures. The parameters of the data generating processes are displayed in Table 3.

We evaluated our approach with 1500 replications with each of the parameter combinations. Under the null of a CB, we tested the correct CB specification against the alternative of the unrestricted full rank model. Under the alternative of no CB, we tested for each of the CB specifications.⁴ Table 4 reports the frequencies of the correctly detected models either with the LRT or with model selection using the information criteria. Analogously to the bivariate case, the LRTs were performed at a 95% confidence level, and the information criteria detect a CB when the IC of the restricted model was lower than the one of the unrestricted model.

We can notice that the size of the LRT when the true rank of the lead coefficient matrix is 2 is similar to the bivariate case. With a finite variance error distribution, the size of the LRT is around 5%, and it increases to around 9% when the variance is infinite ($\lambda =1.5$). We can see that the size of the test decreases in the more restrictive CB specification when the rank of the matrix is 1. For the finite variance cases, the size decreases to 91.9% when $T=500$ and to 93.3% when $T=1000$. The correctly detected model frequency decreases further to 86% in the infinite variance case. Under the alternative of no CB, with finite variance and a sample size of $T=500$, the LRT incorrectly detects a bubble (2 vs. 3) in 30.5% of the cases; however, this frequency decreases to 6.8% when the sample size increases to 1000. Hence, it seems that with a smaller sample size and the finite variance of the error distribution, estimating 8 coefficients in the lead matrix instead of 9 in the unrestricted model still provide a good enough fit to not be rejected by the test. The power of the test for all other model specification is above 99.7%.⁵

As it relates to the model selection using information criteria, BIC outperformed AIC in detecting common bubbles in each of the settings. BIC correctly selected a model with CBs in more than 97.2% of the cases across all model specifications, and the frequencies increased with the sample size and the amount of restricted coefficients. Indeed, it correctly selected a restricted model with a coefficient matrix of rank 1 in at least 99.8% of the cases, whereas AIC selected the correct restricted model in fewer than 88.3%. The frequency decreased with the sample size as did the variance of the errors when the rank of the restricted matrix was closer to full rank. Hence, for the infinite variance case with a sample size $T=500$, its frequency of correctly selected CB model is around 77.5% for each of the CB specifications. Under the alternative of no CB, we observe the same pattern as that for the LRT. In the finite variance case, both information criteria overselected a restricted model with a matrix of rank 2. For a sample size of 500, BIC selected the restricted model in 51.9% of the cases although it decreased to 19.8% when the sample size increased to 1000. AIC, on the other hand, only selected the restricted model in 14.5% with $T=500$ and it even decreased to 3% with $T=1000$. For all other model specifications, both IC selected the correct model in at least 99.4% of the cases.

Overall, the size of the LRT seems to converge to 5% in the finite variance cases when the sample size increases. In the infinite variance cases, the size is around 5 percentage points lower and seems to be less affected by the sample size. The power of the test is above 93% in all model specifications except with $\lambda =3$ and $T=500$, with a restricted model that has only 1 coefficient less to estimate than the unrestricted model (2 vs. 3). For the model selection using information criteria, BIC overall outperforms AIC in correctly detecting a CB but also tends to detect a CB more often than does AIC when there is none in the 2 vs. 3 case with $\lambda =3$.⁶

4. Common Bubbles in Commodity Indices?

We illustrate our strategies in testing for common bubbles in mixed causal–noncausal processes on three commodity price indices: food and beverage, industrial inputs⁷ and fuel (energy)⁸. The sample of 362 data points ranges from January 1992 to January 2022.⁹ We can see from graphs (a) of Figure 1 and Figure 2, which respectively show the series in levels and logs, that the indices seem to follow similar trends. Long-lasting increases and crashes roughly occur at the same time. This could potentially suggest the presence of common bubbles between the series.

Following the work of Hecq and Voisin (2022), we detrend all series using the Hodrick–Prescott filter (hereafter HP filter). Although this approach to obtain stationary time series has been strongly criticized, in particular for the investigation of business cycles, Hecq and Voisin (2022) show that it is a convenient strategy to preserve the bubble features. They also show in a Monte Carlo simulation that this is the filter that preserves the best identification of the MAR($r,s$) model. Giancaterini et al. (2022) arrive at the same conclusion using analytical arguments.

The HP detrended series are displayed on graphs (b) of the two Figures. It can indeed be seen that the dynamics inherent to mixed causal–noncausal processes mentioned above are preserved. The that occurred during the financial crisis of 2007 and the COVID-19 pandemic in 2020, while being of different magnitudes, happened roughly at the same time on all three series. Furthermore, long lasting increases such as the one before the financial crash and the recovery around 2009 or after 2020 are also present in all three index prices.

We first analyzed the series individually. We estimated pseudo causal autoregressive models to identify the order of autocorrelation in each of the detrended series (both in levels and logs). All the models that we identified using BIC were ultimately AR(2) processes. The normality of the errors was rejected for all series: values of the Jarque–Bera statistics ranged between 48 and 253 for the 6 series. The next step was to identify MAR($r,s$) models for all r and s subject to the constraint $p=r+s=2$, namely MAR($2,0$), MAR($1,1$) or MAR($0,2$). Based on the ML estimator with Student’s t-distributed error term, the best fitting model for all six series was the MAR($1,1$) model.

The estimated models are shown in Table 5.¹⁰ For comparison purposes with the trivariate case shown later, we display both the coefficients estimated from the multiplicative

(10)$(1-\varphi L)(1-\psi L^{-1})y_t={\epsilon }_t,\mathrm{with}{\epsilon }_t\sim t\left(\lambda \right),$

and the associated coefficients of the linear form obtained as follows:

(11)$\begin{array}{cc}\hfill y_t& =\frac{\varphi }{1+\varphi \psi }{y}_{t-1}+\frac{\psi }{1+\varphi \psi }{y}_{t+1}+{\epsilon }_t^{*}\hfill \\ \hfill & =b_1{y}_{t-1}+b_2{y}_{t+1}+{\epsilon }_t^{*}.\hfill \end{array}$

It appears that “food and beverage” and “industrial inputs” are both mostly forward looking with lead coefficients close to 0.85 and lag coefficients around 0.4. Conversely, the “fuel index” appears more backward looking with coefficients inverted. Except for the level of industrial inputs, all models have error terms with finite variance, and as expected, one obtains lower variance for the logs of the series. The similar dynamics between food and beverages and industrial inputs could indicate commonalities. The same conclusions can be drawn from the coefficients of the linear form.

For the multivariate investigations, we analyzed both bivariate and trivariate systems. Similarly to the univariate estimation, the strategy consisted of first estimating the pseudo lag order p using a standard VAR(p) for the six bivariate combinations (three in levels and three in logs) and the two trivariate models. Using BIC, all VARs were identified as VAR(2). There were starting value issues when estimating VMARs by ML, meaning that we often reached local maxima. To avoid this issue, we used a large range of starting values to estimate VMAR(1,1) with multivariate Student’s t-distributed errors, and we kept the estimated model with the highest likelihood value.¹¹

The estimated models are shown below in Table 6. We used representation (1) for the estimation, but the coefficients displayed are those of the additive form (3), which are independent of the representation used for the estimations in the following form:

$Y_t=B_1{Y}_{t-1}+B_{-1}{Y}_{t+1}+{\eta }_t,$

where ${\eta }_t$ follows a multivariate Student’s t-distribution with $\lambda$ degrees of freedom and correlation matrix $\Omega$.

In comparison to the coefficients $b_1$ and $b_2$ of the univariate linear forms in Table 5, the directions and magnitudes of the dynamics have been preserved in the multivariate model estimations. From the off-diagonal coefficients of the bivariate models, we can notice that “Food” is impacting both “Indus” and “Fuel” with the lag and the lead, with coefficient magnitudes between 0.11 and 0.47 for the levels. However, in the other direction, the magnitude of the coefficients does not exceed 0.05 for the lag of “Fuel” on “Food”. “Fuel” slightly impacts “Indus” with coefficients of magnitude around 0.1. These dynamics can also be observed in the trivariate model.

To perform the common bubble tests, we estimated VMAR models with restrictions on the lead coefficients matrix as shown in (5).¹² In the trivariate settings, the LRTs and information criteria compare the unrestricted model where the lead matrix has full rank with both CB specifications, namely imposing rank 2 or rank 1 to the lead coefficient matrix.

The results are shown in Table 7. The LRT column displays the LRT statistic, and the IC columns are the difference in the IC values of the restricted and the unrestricted models. As can be seen from the LRTs, the null hypothesis of a common bubble in the bivariate and trivariate models is rejected for all combinations of variables at a confidence level of 95%. All information criteria also indicate a better fit for the models without commonalities since all values are positive. Even for the trivariate cases 2 vs. 3, no bubble can be detected even though in the simulations exercise, the test and information criteria overdetected a CB for the same sample size and degrees of freedom. Hence, while the series seem to follow a similar pattern in the locally explosive episodes throughout the time period, we did not find there to be a significant indication of commonalities in their forward looking components.

There are two possible explanations for these findings. First, our definition requires that bubbles occur at the same time for all the series, whereas graphical evidence may suggest the presence of some degree of nonsynchronicity in the bubble patterns among variables. Second, the series apparently display uncommon explosion rates of the locally explosive episodes.

5. Conclusions

This paper proposes methods to investigate whether the bubble patterns observed in individual series are common to various series. We detected such nonlinear dynamics using recent developments in mixed causal–noncausal autoregressive models. The lead component of the model allows for the capture of locally explosive episodes in a parsimonious and strictly stationary setting. Hence, we employed multivariate mixed causal–noncausal models and applied restrictions to the lead coefficients matrices to test for the presence of commonalities in the forward-looking components of the series. Within a Student t-distribution ML framework, we propose both a LRT and information criteria to detect the presence of common bubbles. In a simulation study, we investigated the finite sample size properties of the proposed approaches, and we found that the BIC performs well when the innovation variances are both finite and infinite. Then, after implementing our approach on three commodity prices, we did not find evidence of commonalities despite the similarities between the series. Our definition of common bubbles requires that all noncausal matrices span the same left null space. A natural extension to our approach would be to relax that hypothesis to investigate nonsynchronous common bubbles, allowing for some adjustment delays along the lines of Cubadda and Hecq (2001).

Author Contributions

All authors contributed equally to the paper. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data in this paper were retrieved from the IMF database: https://www.imf.org/en/Data, accessed on 4 November 2022.

Acknowledgments

Elisa Voisin gratefully acknowledges the University of Rome Tor Vergata for organizing a 3-month research visit, during which this paper was partially written. Gianluca Cubadda gratefully acknowledges a visiting professorship grant from the Graduate School of Business and Economics, Maastricht University. Previous versions of this paper were presented at the York-Maastricht Economics Workshop, Maastricht; the Rome–Waseda Time Series Symposium, Villa Mondragone, Rome; and the CFE 2022 in London. We thank the participants, as well as two referees, for helpful comments and suggestions. The usual disclaimers apply.

Conflicts of Interest

The authors have no declare conflict of interest to declare.

Footnotes

Notes

1. This is the restricted linear form that is used in the ML estimation. Gourieroux and Jasiak (2017) have proposed an alternative approach based on roots inside and outside the unit circle of an autoregressive polynomial.

2. The “axLik” package in R offers a routine for maximizing a given likelihood function with various optimization algorithms. We used the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.

3. Optimization algorithms to maximize the Student’s t multivariate likelihood function are known to be sensitive to starting values and might easily reach local maxima. Since our focus is not on accurate estimation of the models but instead on the detection of commonalities, in order to speed up convergence, we follow previous contributions by employing either the true coefficient matrices when the estimated model correctly imposes k CBs; otherwise, we use an approximation of them with a rank different from $(n-k)$.

4. Results for other tests, such as 1 vs. 2 when the true rank is 2 for instance, are available upon request.

5. Recall from footnote 3 that we employ as starting values an approximation of the true coefficient matrices when the estimated model has a wrong number of CBs. This entails that when the true rank is 3, estimating the restricted models with rank 1 or 2 might encounter convergence issues. This could imply an overestimation of the frequencies displayed in the 2 vs. 3 and 1 vs. 3 when the true rank is 3.

6. Note that Hannan-Quin information criterion $HQC=2Kln\left(ln\left(T\right)\right)-2\ln \left(\widehat{L}\right)$ performs exactly in between BIC and AIC both under the null and under the alternative. We thus omitted this to save space, but results are available upon request.

7. Including agricultural raw materials, such as includes timber, cotton, wool, rubber, and leather.

8. Includes crude oil, natural gas, coal and propane.

9. Data are retrieved from the IMF database. They are price indices with base year 2016.

10. We used various starting values to account for the bimodality of the coefficients (see Bec et al. 2020, for more details).

11. We fixed the starting values for the correlation matrix $\Sigma$ and the degrees of freedom $\lambda$ and performed 100 MLEs based on random lead and lag coefficient matrices fulfilling stationary conditions.

12. We also used 100 combinations of starting values to make sure we obtained the best-fitting models.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Price indices in levels.

Figure 2. Price indices in logs.

Table 1

Monte Carlo parameters for bivariate VMAR(1,1).

$\mathsf{\Phi }=\left[\begin{array}{cc}0.5& 0.1\\ 0.2& 0.3\end{array}\right]$	$\mathsf{\Sigma }=\left[\begin{array}{cc}4& 0.5\\ 0.5& 1\end{array}\right]$
$T=\left\{500,1000\right\}$
$\lambda =\left\{1.5,3\right\}$
$\mathsf{\Psi }=\left\{\begin{array}{cc}\left[\begin{array}{cc}0.3& 0.25\\ 0.6& 0.5\end{array}\right]=\left[\begin{array}{c}1\\ 2\end{array}\right]\left[\begin{array}{cc}0.3& 0.25\end{array}\right]\hfill & (H_0:\mathrm{CB})\hfill \\ \\ \left[\begin{array}{cc}0.1& 0.4\\ 0.6& 0.5\end{array}\right]\hfill & (H_1:\mathrm{no}\mathrm{CB})\hfill \end{array}\right.$

Table 2

Monte Carlo results for N = 2.

	$\lambda =3$
	T = 500				T = 1000
DGP	LR test	BIC	AIC		LR test	BIC	AIC
With CB (rank 1)	0.946	0.989	0.838		0.944	0.993	0.834
Without CB (rank 2)	0.999	0.994	1.000		1.000	1.000	1.000
	$\lambda =1.5$
	T = 500				T = 1000
DGP	LR test	BIC	AIC		LR test	BIC	AIC
With CB (rank 1)	0.913	0.968	0.779		0.914	0.977	0.783
Without CB (rank 2)	0.999	0.999	0.999		1.000	1.000	1.000

Based on 3000 iterations. All results are the frequencies of the correctly detected models. The LR test was performed at a 95% confidence level. For the IC, the favoured model was the one with the lowest IC value. The ranks refer to the rank of the lead coefficient in the DGP.

Table 3

Monte Carlo parameters for trivariate VMAR(1,1).

$\mathsf{\Phi }=\left[\begin{array}{ccc}0.5& 0.1& 0.2\\ 0.2& 0.3& 0.1\\ 0.1& 0.4& 0.6\end{array}\right]$	$\Sigma =\left[\begin{array}{ccc}2& 0.5& 0.5\\ 0.5& 1& 0.5\\ 0.5& 0.5& 4\end{array}\right]$
$T=\left\{500,1000\right\}$
$\lambda =\left\{1.5,3\right\}$
$\mathsf{\Psi }=\left\{\begin{array}{cc}\left[\begin{array}{ccc}0.3& 0.1& 0.1\\ 0.2& 0.3& 0.4\\ 0.7& 0.35& 0.4\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 2& 0.5\end{array}\right]\left[\begin{array}{ccc}0.3& 0.1& 0.1\\ 0.2& 0.3& 0.4\end{array}\right]\hfill & (H_0:1\mathrm{CB}\mathrm{feature})\hfill \\ \\ \left[\begin{array}{ccc}0.15& 0.25& 0.4\\ 0.3& 0.5& 0.8\\ 0.075& 0.125& 0.2\end{array}\right]=\left[\begin{array}{c}1\\ 2\\ 0.5\end{array}\right]\left[\begin{array}{ccc}0.15& 0.25& 0.4\end{array}\right]\hfill & (H_0:2\mathrm{CB}\mathrm{features})\hfill \\ \\ \left[\begin{array}{ccc}0.3& 0.2& 0.1\\ 0.2& 0.5& 0.4\\ 0.7& 0.125& 0.2\end{array}\right]\hfill & (H_1:\mathrm{no}\mathrm{CB}\mathrm{feature})\hfill \end{array}\right.$

Table 4

Monte Carlo results for N = 3.

			$\lambda =3$
			T = 500				T = 1000
$rank\left(\mathsf{\Psi }\right)$	Rank test		LR	BIC	AIC		LR	BIC	AIC
2	2 vs 3		0.944	0.984	0.817		0.951	0.992	0.843
1	1 vs 3		0.919	1.000	0.871		0.933	1.000	0.883
3	2 vs 3		0.695	0.481	0.855		0.932	0.802	0.970
	1 vs 3		1.000	1.000	1.000		1.000	1.000	1.000
			$\lambda =1.5$
			T = 500				T = 1000
$rank\left(\mathsf{\Psi }\right)$	Rank test		LR	BIC	AIC		LR	BIC	AIC
2	2 vs. 3		0.915	0.972	0.775		0.907	0.978	0.776
1	1 vs 3		0.857	0.998	0.774		0.860	0.999	0.783
3	2 vs. 3		0.997	0.994	0.999		1.000	1.000	1.000
	1 vs. 3		1.000	1.000	1.000		1.000	1.000	1.000

Based on 1500 iterations. All results are the frequencies of correctly detected models. The LR test was performed at a 95% confidence level. For the IC, the favored model is the one with the lowest IC value. The ranks refer to the rank of the lead coefficient. $rank\left(\mathsf{\Psi }\right)$ is the rank of the lead coefficient matrix in the DGP.

Table 5

Estimated coefficients of univariate MAR(1,1) models.

Variable	Estimated Coefficients
	Multiplicative			Linear
	$\mathit{\varphi }$	$\mathit{\psi }$	$\mathit{\lambda }$	${\mathit{b}}_1$	${\mathit{b}}_2$
Food & Beverage	0.38	0.85	3.70	0.29	0.64
log(Food & Beverage)	0.34	0.86	5.47	0.26	0.67
Industrial inputs	0.43	0.87	1.66	0.31	0.63
log(Industrial inputs)	0.42	0.89	4.62	0.31	0.65
Fuel (energy)	0.87	0.44	2.20	0.63	0.32
log(Fuel)	0.83	0.48	4.95	0.59	0.34

The coefficients in the multiplicative form are the estimated coefficients from Equation (10). The linear coefficients are the ones obtained by multiplying the estimated factors of the multiplicative as in (11).

Table 6

Estimated coefficients on the multivariate VMAR(1,1) models.

$B_1$	$B_{-1}$	$\Omega$	$\lambda$
Food and Indus
$\left[\begin{array}{cc}0.28& 0.01\\ 0.26& 0.27\end{array}\right]$	$\left[\begin{array}{cc}0.65& -0.02\\ -0.11& 0.65\end{array}\right]$	$\left[\begin{array}{cc}1.32& 0.16\\ 0.16& 3.35\end{array}\right]$	2.49
Food and Fuel
$\left[\begin{array}{cc}0.35& 0.05\\ 0.47& 0.52\end{array}\right]$	$\left[\begin{array}{cc}0.55& -0.04\\ -0.40& 0.40\end{array}\right]$	$\left[\begin{array}{cc}1.42& 0.87\\ 0.87& 12.90\end{array}\right]$	3.01
Indus and Fuel
$\left[\begin{array}{cc}0.29& 0.01\\ -0.11& 0.47\end{array}\right]$	$\left[\begin{array}{cc}0.63& 0.03\\ 0.09& 0.48\end{array}\right]$	$\left[\begin{array}{cc}2.22& 1.50\\ 1.50& 7.17\end{array}\right]$	1.67
Food, Indus, and Fuel
$\left[\begin{array}{ccc}0.27& 0.01& 0.03\\ 0.27& 0.25& 0.02\\ 0.30& -0.10& 0.56\end{array}\right]$	$\left[\begin{array}{ccc}0.64& -0.02& -0.02\\ -0.16& 0.66& -0.01\\ -0.27& 0.12& 0.37\end{array}\right]$	$\left[\begin{array}{ccc}1.34& 0.12& 0.55\\ 0.12& 3.29& 2.20\\ 0.55& 2.20& 10.02\end{array}\right]$	2.28
$B_1$	$B_{-1}$	${10}^3\Omega$	$\lambda$
Food and Indus
$\left[\begin{array}{cc}0.25& 0.01\\ 0.22& 0.24\end{array}\right]$	$\left[\begin{array}{cc}0.69& -0.03\\ -0.12& 0.70\end{array}\right]$	$\left[\begin{array}{cc}0.27& 0.03\\ 0.03& 0.46\end{array}\right]$	6.30
Food and Fuel
$\left[\begin{array}{cc}0.25& 0.02\\ 0.16& 0.38\end{array}\right]$	$\left[\begin{array}{cc}0.67& -0.02\\ -0.14& 0.55\end{array}\right]$	$\left[\begin{array}{cc}0.25& 0.06\\ 0.06& 1.21\end{array}\right]$	5.23
Indus and Fuel
$\left[\begin{array}{cc}0.26& 0.04\\ -0.09& 0.56\end{array}\right]$	$\left[\begin{array}{cc}0.67& -0.01\\ 0.09& 0.37\end{array}\right]$	$\left[\begin{array}{cc}0.42& 0.24\\ 0.24& 1.19\end{array}\right]$	4.77
Food, Indus, and Fuel
$\left[\begin{array}{ccc}0.88& -0.17& -0.02\\ -0.04& 0.27& 0.07\\ -0.04& 0.06& 0.58\end{array}\right]$	$\left[\begin{array}{ccc}0.21& 0.15& 0.00\\ 0.13& 0.76& -0.08\\ 0.02& 0.05& 0.33\end{array}\right]$	$\left[\begin{array}{ccc}0.32& 0.02& 0.07\\ 0.02& 0.51& 0.26\\ 0.07& 0.26& 1.35\end{array}\right]$	6.15

Table 7

Common bubble detection on multivariate combinations of the variables.

Levels
Food	Indus	Fuel	Rank test	LRT	BIC	AIC
∎	∎		1 vs. 2	25.93	20.04	23.93
∎		∎	1 vs. 2	59.96	54.07	57.96
	∎	∎	1 vs. 2	70.49	64.59	68.49
∎	∎	∎	2 vs. 3	16.26	10.37	14.26
			1 vs. 3	88.12	64.55	80.12
Logs
Food	Indus	Fuel	Rank test	LRT	BIC	AIC
∎	∎		1 vs. 2	16.04	10.15	14.04
∎		∎	1 vs. 2	34.36	28.47	32.36
	∎	∎	1 vs. 2	46.05	40.16	44.05
∎	∎	∎	2 vs. 3	15.81	9.92	13.81
			1 vs. 3	75.01	51.44	67.01

LRT is the likelihood ratio test statistic. For the bivariate models, the critical value of the LRT at 95% confidence level is 3.41. For the trivariate models, the critical values are 3.841 and 9.488 for 2 vs. 3 and 1 vs. 3, respectively. The column BIC and AIC show the difference between the restricted and unrestricted information criteria.

References

Bec, Frédérique; Nielsen, Heino Bohn; Saïdi, Sarra. Mixed causal–noncausal autoregressions: Bimodality issues in estimation and unit root testing 1. Oxford Bulletin of Economics and Statistics; 2020; 82, pp. 1413-28. [DOI: https://dx.doi.org/10.1111/obes.12372]

Bernardini, Emmanuela; Cubadda, Gianluca. Macroeconomic forecasting and structural analysis through regularized reduced-rank regression. International Journal of Forecasting; 2015; 31, pp. 682-91.

Cubadda, Gianluca; Hecq, Alain; Telg, Sean. Detecting co-movements in non-causal time series. Oxford Bulletin of Economics and Statistics; 2019; 81, pp. 697-715. [DOI: https://dx.doi.org/10.1111/obes.12281]

Cubadda, Gianluca; Hecq, Alain. On non-contemporaneous short-run co-movements. Economics Letters; 2001; 73, pp. 389-97. [DOI: https://dx.doi.org/10.1016/S0165-1765(01)00514-6]

Cubadda, Gianluca; Hecq, Alain. Testing for common autocorrelation in data-rich environments. Journal of Forecasting; 2011; 30, pp. 325-35. [DOI: https://dx.doi.org/10.1002/for.1186]

Cubadda, Gianluca; Hecq, Alain. Dimension reduction for high dimensional vector autoregressive models. Oxford Bulletin of Economics and Statistics; 2022a; 84, pp. 1123-52. [DOI: https://dx.doi.org/10.1111/obes.12506]

Cubadda, Gianluca; Hecq, Alain. Reduced rank regression models in economics and finance. Oxford Research Encyclopedia of Economics and Finance; 2022b; [DOI: https://dx.doi.org/10.1093/acrefore/9780190625979.013.677]

Engle, Robert F.; Granger, Clive W. J. Co-integration and error correction: Representation, estimation, and testing. Econometrica: Journal of the Econometric Society; 1987; 55, pp. 251-76. [DOI: https://dx.doi.org/10.2307/1913236]

Engle, Robert F.; Susmel, Raul. Common volatility in international equity markets. Journal of Business & Economic Statistics; 1993; 11, pp. 167-76.

Engle, Robert F.; Kozicki, Sharon. Testing for common features. Journal of Business & Economic Statistics; 1993; 11, pp. 369-80.

Engle, Robert F.; Hylleberg, Svend. Common seasonal features: Global unemployment. Oxford Bulletin of Economics and Statistics; 1996; 58, pp. 615-30. [DOI: https://dx.doi.org/10.1111/j.1468-0084.1996.mp58004003.x]

Engsted, Tom; Nielsen, Bent. Testing for rational bubbles in a coexplosive vector autoregression. The Econometrics Journal; 2012; 15, pp. 226-54. [DOI: https://dx.doi.org/10.1111/j.1368-423X.2012.00369.x]

Giancaterini, Francesco; Hecq, Alain; Morana, Claudio. Is climate change time reversible?. Econometrics; 2022; 10, 36. [DOI: https://dx.doi.org/10.3390/econometrics10040036]

Gourieroux, Christian; Jasiak, Joann. Noncausal vector autoregressive process: Representation, identification and semi-parametric estimation. Journal of Econometrics; 2017; 200, pp. 118-34. [DOI: https://dx.doi.org/10.1016/j.jeconom.2017.01.011]

Gouriéroux, Christian; Zakoïan, Jean-Michel. Local explosion modelling by non-causal process. Journal of the Royal Statistical Society: Series B (Statistical Methodology); 2017; 79, pp. 737-56. [DOI: https://dx.doi.org/10.1111/rssb.12193]

Hecq, Alain; Voisin, Elisa. Predicting bubble bursts in oil prices during the COVID-19 pandemic with mixed causal-noncausal models. arXiv; 2022; arXiv: 1911.10916

Hendry, David F.; Massmann, Michael. Co-breaking: Recent advances and a synopsis of the literature. Journal of Business & Economic Statistics; 2007; 25, pp. 33-51.

Issler, Joao Victor; Vahid, Farshid. Common cycles and the importance of transitory shocks to macroeconomic aggregates. Journal of Monetary Economics; 2001; 47, pp. 449-75. [DOI: https://dx.doi.org/10.1016/S0304-3932(01)00052-6]

Lanne, Markku; Saikkonen, Pentti. Noncausal vector autoregression. Econometric Theory; 2013; 29, pp. 447-81. [DOI: https://dx.doi.org/10.1017/S0266466612000448]

Vahid, Farshid; Engle, Robert F. Common trends and common cycles. Journal of Applied Econometrics; 1993; 8, pp. 341-60.