Authors: Marcio A. Diniz, Carlos A. B. Pereira, Julio M. Stern
Categories: Article, Bayesian inference, cointegration, hypothesis testing, time series, unit root
Source: Entropy
Doi: 10.3390/e22090968
Authors: Marcio A. Diniz, Carlos A. B. Pereira, Julio M. Stern
To perform statistical inference for time series, one should be able to assess if they present deterministic or stochastic trends. For univariate analysis, one way to detect stochastic trends is to test if the series has unit roots, and for multivariate studies it is often relevant to search for stationary linear relationships between the series, or if they cointegrate. The main goal of this article is to briefly review the shortcomings of unit root and cointegration tests proposed by the Bayesian approach of statistical inference and to show how they can be overcome by the Full Bayesian Significance Test (FBST), a procedure designed to test sharp or precise hypothesis. We will compare its performance with the most used frequentist alternatives, namely, the Augmented Dickey–Fuller for unit roots and the maximum eigenvalue test for cointegration.
Several times series present deterministic or stochastic trends, which imply that the effects of these trends on the level of the series are permanent. Consequently, the mean and variance of the series will not be constant and will not revert to a long-term value. This feature reflects the fact that the stochastic processes generating these series are not (weakly) stationary, imposing problems to perform inductive inference using the most traditional estimators or predictors. This is so because the usual properties of these procedures will not be valid under such conditions.
Therefore, when modeling non-stationary time series, one should be able to properly detrend the used series, either by directly modeling the trend by deterministic functions, or by transforming the series to remove stochastic trends. To determine which strategy is the suitable solution, several statistical tests were developed since the 1970s by the frequentist school of statistical inference.
The Augmented Dickey–Fuller (ADF) test is one of the most popular tests used to assess if a time series has a stochastic trend or, for series described by auto-regressive models, if they have a unit root. When one is searching for long term relationships between multiple series under analysis, it is crucial to know if there are stationary linear combinations of these series, i.e., if the series are cointegrated. Cointegration tests were developed, also by the frequentist school, in the late 1980s [1] and early 1990s [2]. Only in the late 1980s did the Bayesian approach to test the presence of unit roots start to be developed.
Both unit root and cointegration tests may be considered tests on precise or sharp hypotheses, i.e., those in which the dimension of the parameter space under the tested hypothesis is smaller than the dimension of the unrestricted parameter space. Testing sharp hypotheses poses major difficulties for either the frequentist or Bayesian paradigms, such as the need to eliminate nuisance parameters.
The main goal of this article is to briefly review the shortcomings of the tests proposed by the Bayesian school and how they can be overcome by the Full Bayesian Significance Test (FBST). More specifically, we will compare its performance with the most used frequentist alternatives, the ADF for unit roots, and the maximum eigenvalue test for cointegration. Since this is a review article, it is important to remark that the results presented here were published elsewhere by the same authors, see [3,4].
To accomplish this objective, we will define the FBST in the next section, also showing how it can be implemented in a general context. The following section discusses the problems of testing the existence of unit roots in univariate time series and how the Bayesian tests approach the problem. Section 4 then shows how the FBST is applied to test if a time series has unit roots and illustrates this with applications on a real data set. In the sequel, we discuss the Bayesian alternatives to cointegration tests and then apply the FBST to test for cointegration using real data sets. We conclude with some remarks and possible extensions for future work.
The Full Bayesian Significance Test was proposed in [5] mainly to deal with sharp hypotheses. The procedure has several properties, see [6,7], most interestingly the fact that it is only based on posterior densities, thus avoiding the necessity of complications such as the elimination of nuisance parameters or the adoption of priors with positive probabilities attached to sets of zero Lebesgue measure.
We shall consider general statistical models in which the parameter space is denoted by Θ⊆Rm, m∈N. A sharp hypothesis H assumes that θ, the parameter vector of the chosen statistical model, belongs to a sub-manifold ΘH of smaller dimensions than Θ. This implies, for continuous parameter spaces, that the subset ΘH has null Lebesgue measure whenever H is sharp. The sample space, the set of all possible values of the observable random variables (or vectors), is here denoted by X.
Following the Bayesian paradigm, let h(·) be a probability prior density over Θ, x∈X, the observed sample (scalar or vector), and L(·∣x) the likelihood derived from data x. To evaluate the Bayesian evidence based on the FBST, the sole relevant entity is the posterior probability density for θ given x, g(θ∣x)∝h(θ)·L(θ∣x).
It is important to highlight that the procedure may be used when the parameter space is discrete. However, when the posterior probability distribution over Θ is absolutely continuous, the FBST appears as a more suitable alternative to significance hypothesis testing. For notational simplicity, we will denote ΘH by H in the sequel.
Let r(θ) be a reference density on Θ such that the function s(θ)=g(θ∣x)/r(θ) is a relative surprise, (see [8], pp. 145–146) function. The reference density is important because it guarantees that the FBST is invariant to reparametrizations, even when r(θ) is improper, see [6,9]. Thus, when considering r(θ) proportional to a constant, the surprise function will be, in practical terms, equivalent to the posterior distribution. For the applications considered in this article, we will use the improper uniform density as reference density on Θ. The authors of [10] remark that it is possible to generalize the procedure using other reference densities such as neutral, invariant, maximum-entropy or non-informative priors, if they are available and desirable.
Notice that the tangent set Tx is the highest relative surprise set, that is, the set of points of the parameter space with higher relative surprise than any point in H, being tangential to H in this sense. This approach takes into consideration the statistical model in which the hypothesis is defined, using several components of the model to define an evidential measure favoring the hypothesis.
Definition 2 sets ev¯ as the posterior probability of the tangent set that is interpreted as an evidence value against H. Hence, the evidence value supporting H is the complement of ev¯, namely, ev=1−ev¯. Notwithstanding, ev is not evidence against A:θ∉ΘH, the alternative hypothesis (which is not sharp anyway). Equivalently, ev¯ is not evidence in favor of A, although it is against H.
Thus, we are left with the problem of deciding the critical level evc for each particular application. We briefly discuss this and other practical issues in the following subsection.
Since ev (also called e-value) is a statistic, it has a sampling distribution derived from the adopted statistical model and in principle this distribution could be used to find a threshold value. If the likelihood and the posterior distribution satisfy certain regularity conditions. See [11], p. 436. [12] proved that, asymptotically, there is a relationship between ev and the p-values obtained from the frequentist likelihood ratio procedure used to test the same hypotheses. This fact provides a way to find, at least asymptotically, a critical value to ev to reject the hypothesis being tested.
In a recent review [7], the authors discuss different ways to provide a threshold for ev. Among these alternatives, we highlight the standardized e-value, which follows, asymptotically, the uniform distribution on (0,1). See also [13] for more on the standardized version of ev.
One could also try to define the FBST as a Bayes test derived from a particular loss function and the respective minimization of the posterior expected loss. Following this strategy, [10] showed that there are loss functions which result in ev as a Bayes estimator of ϕ=IH(θ), where IA(x) denotes the indicator function, being equal to one if x∈A and zero otherwise, x∉A. Hence, the FBST is in fact a Bayes procedure in the formal sense as defined by Wald in [14].
To compute the evidence value supporting H defined in the last section, we need to follow the steps showed in Table 1. Appendix A provides detailed information about the computational resources and codes used to implement the FBST in the examples presented in this work. After defining the statistical model and prior, it is simple to find the surprise function, s(θ). In step 3, one should find the point of the parameter space in H that maximizes s(θ), that is, to solve a problem of constrained numerical maximization. In several applications, this step does not present a closed form solution, requiring the use of numerical optimizers.
Step 4 involves the integration of the posterior distribution on a subset of Θ, the tangent set Tx that can be highly complex. Once more, since in many cases it is fairly difficult to find an explicit expression for Tx, one may use various numerical techniques to compute the integral. If it is possible to generate random samples from the posterior distribution, Monte Carlo integration provides an estimate of ev, as we will show in this work. Another alternative is to use approximation techniques, such as those proposed in [15], based on a Laplace approximation. We discuss how to implement such approximations for unit root and cointegration tests in [3,4].
Before presenting the Bayesian procedures used to test the presence of unit roots, let us fix notation. We will denote by yt the t-th value of a univariate time series observed in t=1,…,T+p dates, where T and p are positive integers. The usual approach is to assume that the series under analysis is described by an auto-regressive process with p lags, AR(p), meaning that the data generating process is fully described by a stochastic difference equation of order p, possibly with an intercept or drift and a deterministic linear trend, i.e., (3)yt=μ+δ·t+ϕ1yt−1+…+ϕpyt−p+εt with εt i.i.d. N(0,σ2) for t=1,…,T+p. Using the lag or backshift operator B, we denote yt−k by Bkyt, allowing us to rewrite (3) as (4)(1−ϕ1B−…ϕpBp)yt=μ+δ·t+εt where ϕ(B)=(1−ϕ1B−…ϕpBp) is the autoregressive polynomial. The difference Equation (3) will be stable, implying that the process generating {yt}t=1T+p is (weakly) stationary, whenever the roots of the characteristic polynomial ϕ(z), z∈C, lie outside the unit circle, since there may be complex roots. The set of polynomial operators, such as lag polynomials like ϕ(B), induces an algebra that is isomorphic to the algebra of polynomials in real or complex variables, see [16].
If some of the roots lie exactly on the unit circle, it is said that the process has unit roots. In order to test such a hypothesis statistically, (3) is rewritten as (5)Δyt=μ+δ·t+Γ0yt−1+Γ1Δyt−1+…+Γp−1Δyt−p+1+εt where Δyt=yt−yt−1, Γ0=ϕ1+…+ϕp−1 and Γi=−∑j=i+1pϕj, for i=1,…,p−1. If the generating process has only one unit root, one root of the complex polynomial ϕ(z), 1−ϕ1z−ϕ2z2…ϕpzp, is equal to one, meaning that 1−ϕ1−ϕ2−…−ϕp=0 i.e., ϕ(1)=0, and all the other roots are on or outside the unit circle. In this case, Γ0=0, the hypothesis that will be tested when modeling (5). Even though tests based on these assumptions verify if the process has a single unit root, there are generalizations based on the same principles that test the existence of multiple unit roots, see [17].
The search for Bayesian unit root tests began in the late 1980s. As far as we know, [18,19] were the first works to propose a Bayesian approach for unit root tests. The frequentist critics of these articles received a proper answer in [20,21], generating a fruitful debate that produced a long list of papers in the literature of Bayesian time series. A good summary of the debate and the Bayesian papers that resulted from it is presented in [22]. We will present here only the most relevant strategies proposed by the Bayesian school to test for unit roots.
Let θ=(ρ,ψ) be the parameters vector, in which ρ=∑i=1pϕi and ψ=(μ,δ,Γ1,…,Γp−1). Assuming σ2 fixed, the prior density for θ can be factorized as h(θ)=h0(ρ)·h1(ψ∣ρ).
The marginal likelihood for ρ, denoted by Lm, Lm(ρ∣y)∝∫ΨL(θ∣y)·h1(ψ∣ρ)dψ. where y={yt}t=1T+p is the observations vector, L(θ|y) the full likelihood, and Ψ the support of the random vector ψ. This marginal likelihood, associated with a prior for ρ, is the main ingredient used by standard Bayesian procedures to test the existence of unit roots. Even though the procedure varies among authors according to some specific aspects, mentioned below, basically all of them use Bayes factors and posterior probabilities.
One important issue is the specification of the null some authors, starting from [23], consider H0:ρ=1 against H1:ρ<1. Starting from [24], this is the way the frequentist school addresses the problem, but following this approach no explosive value for ρ is considered. The decision theoretic Bayesian approach solved the problem using the posterior probabilities ratio or Bayes B01=Lm(ρ=1∣y)∫01Lm(ρ∣y)·h0(ρ)dρ.
Advocates of this solution argue that one of the advantages of this approach is that the null and the alternative hypotheses are given equal weight. However, the expression above is not defined if h0(ρ) is not a proper density since the denominator of the Bayes factor is equal to the predictive density, defined just if h0(ρ) is a proper density. There are also problems if Lm(ρ=1|y) is zero or infinite.
The problem is approached by [20,25] by testing H0:ρ≥1 against H1:ρ<1, considering explicitly explosive values for ρ. The main advantage of this strategy is the possibility to compute posterior probabilities like P(ρ>1∣y)=∫1∞gm(ρ∣y)dρ defined even for improper priors on ρ, where gm is the marginal posterior for ρ.
In [26], the authors do not choose ρ as the parameter of interest, examining instead the largest absolute value of the roots of the characteristic polynomial and then verifying if it is smaller or larger than one. Usually, this value is slightly smaller than ρ, but the authors argue that this small difference may be important. When this approach is used, unit roots are found less frequently. For an AR(3) model with a constant and deterministic trend, [26] derives the posterior density for the dominant root for the 14 series used in [27] and concluded the for eleven of the series, the dominant root was smaller than one, that is to say, the series were trend-stationary. These results were based on flat priors for the autoregressive parameters and the deterministic trend coefficient.
Another controversy is about the prior over ρ: [20] argues that the difference between the results given by the frequentist and Bayesian inferences is due to the fact that the flat prior proposed in [18] overweights the stationary region of ρ. Hence, he derived a Jeffreys prior for the AR(1) this prior quickly diverges as ρ increases and becomes larger than one. The obtained posterior led to the same results of [27], which will be discussed in detail in the following section. The critics of the approach adopted by Phillips in [20] judged the Jeffreys prior as unrealistic, from a subjective point of view. See the comments on Phillips’s paper on the Journal of Applied Econometrics, volume 6, number 4, 1991. The subsequent papers of the same number support the Bayesian approach. This is a nonsensical objection if one considers that the Jeffreys prior is crucial to ensure an invariant inferential procedure, and invariance is a highly desirable property, for either objective or subjective reasons. See [28] for more on invariance in physics and statistical models.
A final controversial point concerns the modeling of initial observations. If the likelihood explicitly models the initial observed values (it is an exact likelihood), the process is implicitly considered stationary. In fact, when it is known that the process is stationary, and it is believed that the data generating process is working for a long period, it is reasonable to assume that the parameters of the model determine the marginal distribution of the initial observations. In the simplest AR(1) model, this would imply that y1∼N(0,σ2/(1−ρ2)). In this scenario, to perform the inference conditional on the first observation would discard relevant information. On the other hand, there is no marginal distribution defined for y1 if the generating process is not stationary. Then, it is valid to use a likelihood conditional on initial observations. For the models presented here, we always work with the conditional likelihood. As argued in [18], inferences for stationary models are little affected by using conditional likelihoods, especially for large samples. He compares these inferences with the ones based on exact likelihoods under explicit modeling for initial observations.
We will now describe how to use the FBST to test for the presence of unit roots referring to the general model (5). It is also possible to include q∈N moving average terms in (3) to model the process, a case that will not be covered in this article but that, in principle, shall not imply major problems for the FBST.
See Equation (5), where εt∼i.i.d.N(0,σ2) for t=1,…,T+p, recalling also that the hypothesis being tested is Γ0=0. We slightly change the notation of the last section now using ψ to denote the vector (μ,δ,Γ0,…,Γp−1) and setting θ=(ψ,σ).
Recalling the steps to implement the FBST displayed in Table 1, we have just specified the statistical model. The likelihood, conditional on the first p observations, derived from the Gaussian model is (6)L(θ∣y)=(2π)−T/2σ−Texp−12σ2·∑t=p+1T+pεt2, in which εt=Δyt−μ−δ·t−Γ0yt−1−Γ1Δyt−1−…−Γp−1Δyt−p+1. To complete step 1 of Table 1, we need a prior distribution for θ. For all the series modeled in this article, we will use the following non informative (7)h(θ)=h(ψ,σ)∝1/σ.
We are aware of the problems caused by improper priors applied to this problem when one uses alternative approaches, like those mentioned by [22]. However, one of our goals is to show how the FBST can be implemented even for a potentially problematic prior like this one. To write the posterior, we use the following ΔY=Δyp+1Δyp+2⋮ΔyT+p,X=1p+1ypΔyp…Δy21p+2yp+1Δyp+1…Δy3⋮⋮⋮⋮⋮⋮1T+pyT+p−1ΔyT+p−1…ΔyT+1,ψ=μδΓ0⋮Γp−1, being ΔY of dimension T×1, X of dimension T×(p+2) and ψ, (p+2)×1. Thanks to this notation, we can write, using primes to denote transposed ∑t=p+1T+pεt2=(ΔY−Xψ)′(ΔY−Xψ)=(ΔY−ΔY^)′(ΔY−ΔY^)+(ψ−ψ^)′X′X(ψ−ψ^), where ψ^=(X′X)−1X′·ΔY is the ordinary least squares (OLS) estimator of ψ and ΔY^=Xψ^ its prediction for ΔY. Thus, the full posterior is (8)g(θ∣y)∝σ−(T+1)exp−12σ2[(ΔY−ΔY^)′(ΔY−ΔY^)+(ψ−ψ^)′X′X(ψ−ψ^)], a Normal-Inverse Gamma density.
Step 2 demands a reference density in order to define the relative surprise function. Since we will use the improper density r(θ)∝1, the surprise function will be equivalent to the posterior distribution in our applications. Given this, to find s* (Step 3), we need to find the maximum value of the posterior under the hypothesis being tested, in our case, Γ0=0.
This maximization step is fairly simple to implement given the modeling choices made Gaussian likelihood, non informative prior and reference density proportional to a constant. The restricted (assuming H) posterior distribution is (9)gr(θr∣y)∝σ−(T+1)exp−12σ2[(ΔY−ΔY^r)′(ΔY−ΔY^r)+(ψr−ψ^r)′Xr′Xr(ψr−ψ^r)], in which θr=(ψr,σ), ψr being vector ψ without Γ0, Xr=1p+1Δyp…Δy21p+2Δyp+1…Δy3⋮⋮⋮⋮⋮1T+pΔyT+p−1…ΔyT+1,ψ^r=(Xr′Xr)−1Xr′·ΔY,andΔY^r=Xrψ^r, that is, Xr is simply matrix X above without its third column, since under H:Γ0=0 and the coefficient of the third column of X is Γ0—see Equation (5)—ψ^r is a least squares estimator of ψr and ΔY^r denotes the predicted values for ΔY given by the restricted model. From (9), it is easy to show that the maximum a posteriori (MAP) estimator of θr is given by (ψ^r,σ^r), with σ^r=(ΔY−ΔY^r)′(ΔY−ΔY^r)T+1.
Plugging the values of ψ^r and σ^r into (9), we find s*, as requested in Step 3. Step 4 will also be easy to implement thanks to structure of the models assumed in this section. Since the full posterior, (8), is a Normal-Inverse Gamma density, a simple Gibbs sampler allows us to obtain a random sample from such distribution, suggesting a Monte Carlo approach to compute ev¯. From (8), the conditional posteriors of ψ and σ are, respectively, (10)gψ(ψ∣σ,y)∝N(ψ^,σ2(X′X)−1) (11)gσ2(σ2∣ψ,y)∝IGT+12,H in which H=0.5[(ΔY−ΔY^)′(ΔY−ΔY^)+(ψ−ψ^)′X′X(ψ−ψ^)], IG denotes the Inverse-Gamma distribution and ψ^ is the OLS estimator of ψ, as above. Appendix B brings the parametrization and the probability density function of the Inverse-Gamma distribution. With a sizable random sample from the full posterior, we estimate ev¯ as the proportion of sampled vectors that generate a value for the posterior greater than s*, found in Step 3. Hence, in Step 5, we only compute one minus the estimate of ev¯ found in Step 4. The whole procedure is summarized in Table 2. For the implementations in this article we sampled 51,000 vectors from (8) and discarded the first 1,000 as a burn-in sample.
We implemented the FBST as described above to test the presence of unit roots in 14 U.S. macroeconomic time series, all with annual frequency, first mentioned in [27]. We used the extended series, analyzed in [23]. Appendix A brings more information on the data set and the computational resources and codes used to obtain the results displayed in Table 3 below.
Table 3 reports the names of the tested series, the number of available observations or sample size, the adopted value for p—as denoted in Equation (8)—if a linear (deterministic) trend was included in the model or not, the ADF test statistic and its respective p-value. We have used the computer package described in [29] to find the ADF p-values, available in the R library urca. The last two columns bring the posterior probability of non-stationarity, Γ0≥0, and the FBST e-values for the specified models. In order to obtain comparable results, we have adopted the models chosen by [22] for all the series. All the models considered the intercept or constant term, μ in (9).
The results show that the non-stationary posterior probabilities are quite distant from the ADF p-values. These results were highlighted in [18,19]. Considering the simplest AR(1) model, they argued that, once frequentist inference is based on the distribution of ρ^|ρ=1, the non-stationary posterior probabilities provide counterintuitive conclusions since the referred distribution is skewed. Their main argument is that Bayesian inference uses a distribution (marginal posterior of ρ) that is not skewed.
As mentioned before, ref. [20] claims that the difference in results between frequentist and Bayesian approaches is due to the flat prior that puts much weight on the stationary region. He proposed the use of Jeffreys priors, which restored the conclusions drawn by the frequentist test. Phillips argued that the flat prior was, actually, informative when used in time series models like those for unit root tests. Using simulations, he shows that “ [the use of a] flat prior has a tendency to bias the posterior towards stationarity. … even when [the estimate] is close to unity, there may still be a non negligible downward bias in the [flat] posterior probabilities”. Notwithstanding, the e-values reported in the last column are quite close to the ADF p-values even using the flat prior criticized by Phillips.
Before starting our brief review of the most relevant Bayesian cointegration tests, we fix notation and present the definitions to which we will refer in the sequel.
All the tests mentioned here are based on the following multivariate framework. Let Yt=[y1t…ynt]′ be a vector with n∈N time series, all of them assumed to be integrated of order d∈N, i.e., have d unit roots. The series are said to be cointegrated if there is a nontrivial linear combination of them that has b∈N unit roots, b<d. We will assume that, as in most applications, d=1 and b=0, meaning that, if the time series in Yt is cointegrated, there is a linear combination a′Yt that is stationary, where a∈Rn is the cointegrating vector. Since the linear combination a′Yt is often motivated by problems found in economics, it is called a long-run equilibrium relationship. The explanation is that non-stationary time series that are related by a long-run relationship cannot drift too far from the equilibrium because economic forces will act to restore the relationship.
Notice also (i) the cointegrating vector is not uniquely determined since, for any scalar s, (s·a) is a cointegrating vector; and (ii) if Yt has more than two series, it is possible that there is more than one cointegrating vector generating a stationary linear combination.
It is assumed that the data generating process of Yt is described by the following vector autoregression with p∈N lags, denoted VAR(p), and given (12)Yt=c+Φ0Dt+Φ1Yt−1+…+ΦpYt−p+Et, in which c is a (n×1) vector of constants, Dt a vector (n×1) with some deterministic variable, such as deterministic trends or seasonal dummies, Φi are (n×n) coefficients matrices and Et is a (n×1) stochastic vector with multivariate normal distribution with null expected value and covariance matrix Ω, denoted Nn(0,Ω). This dynamic model is assumed valid for t=1,…,T+p, the available span of observations of Yt. As in the univariate case, one may include moving average terms in (12), i.e., lags for Et, but this, in principle, would not cause major problems in the Bayesian framework. Model (12) can be rewritten using the lag or backshift operator as (13)(In−Φ1B−…−ΦpBp)Yt=c+Φ0Dt+Et, where Φ(B)=In−Φ1B−…−ΦpBp is the (multivariate) autoregressive polynomial and In denotes the n-dimensional identity matrix. The associate characteristic polynomial in this context will be the determinant of Φ(z), z∈C. If all the roots of the characteristic polynomial lie outside the unit circle, it is possible to show that Yt has a stationary representation—see [30]—such as Equation (13). In order to determine if this is the case, model (12) is rewritten as an (vectorial) error correction model (VECM):(14)ΔYt=c+Φ0Dt+Γ1ΔYt−1+…+Γp−1ΔYt−p+1+ΠYt−1+Et, where ΔYt=[Δy1t…Δynt]′, Γi=−(Φi+1+…Φp) for i=1,2,…,p−1 and Π=−Φ(1)=−(In−Φ1−…−Φp). It is possible to show that, when all the roots of det(Φ(z)) are outside the unit circle, matrix Π in (14) has full rank, i.e., all the n eigenvalues of Π are n non null. If the rank of Π is null, this matrix cannot be distinguished from a null matrix, implying that the series in Yt has at least one unit root and a valid representation is a VAR of order p−1, i.e., model (14) without the term ΠYt−1. It is possible that the series in Yt has two unit roots each, implying that the correct VECM must be written with Δ2Yt as a dependent variable.
Finally, if the (n×n) matrix Π has rank r, 0<r<n, it has n−r non null eigenvalues, implying that the series in Yt has at least one unit root and its valid representation is given by the VECM in Equation (14). In this case, Π=αβ′, where α and β are matrices (n×r) of rank r. Matrix β denotes the one with the cointegrating vectors and matrix α is called the loading matrix, since it contains the weights of the equilibrium relationships. The tests developed in [2] focus on the rank of matrix Π.
The pioneer Bayesian works to study VAR models and reduced rank regressions are [31,32,33]. However, the main concern of these papers is to estimate the model parameters and their (marginal) posterior distributions. The usual approach is to assume a given rank for the long run matrix Π, and proceed with all the computations conditional on the given rank. The Bayesian initiatives to test the rank of the referred matrix are recent, the main reference for Bayesian inference on VECM’s being [34].
To justify inferential procedures based on prespecified ranks of matrix Π, [22] argued that an empirical cointegration analysis should be based on economic theory, which proposes models obeying equilibrium relationships. According to this view, cointegration research should be “confirmatory” rather than “exploratory”. Even though the advocated conditional inference is of simple implementation and very useful for small samples, [22] recognized that tests for the rank of matrix Π should be developed. To our knowledge, few initiatives with this purpose were developed up to now.
One common approach to test sharp hypotheses in the Bayesian framework is by means of Bayes factors. Testing the rank of matrix Π by Bayes factors implies several computational complications and requires the use of proper priors, as shown in [35]. Following an informal approach, [33] obtained the posterior distribution of the ordered eigenvalues of the “squared” long run matrix, Π′·Π, obtained from a VAR model without assuming the existence of cointegration relations. As the long run matrix has a reduced rank, it has some null eigenvalues, and this should be revealed by the fact that the smallest eigenvalues should have a lot of probability mass accumulated on values close to zero. The computations can be made straightforwardly, simulating values for the long run matrix from its (marginal) posterior distribution, which is a matrix t-Student distribution under the non informative prior (16), also considered in the sequel.
Another common procedure is to estimate the rank of Π as the value r that maximizes the (marginal) posterior distribution of the rank. Conditioned on such an estimate, one proceeds to derive the full posterior and eventually estimate the cointegration space, i.e., the linear space spanned by β.
A different approach was proposed by [36], who used the Posterior Information Criterion (PIC), developed in [37], as a criterion to choose the mode of the posterior distribution of the rank of Π. However, as highlighted in [34], one of the advantages of the Bayesian approach is the possibility to incorporate the uncertainty about the parameters in the analysis, represented by the posterior distribution of the rank and, whatever the tool the scientist uses to infer the value of r, it is derived from this posterior distribution.
The authors of [38] nested the reduced rank models in an unrestricted VAR and used Metropolis–Hastings sampling with the Savage–Dickey density ratio—see [39]—to estimate the Bayes Factors of all the models with incomplete rank up to the model with full rank. The Bayes Factor derivation requires the estimation of an error correction factor for the incomplete rank. This factor, however, is not defined for improper priors due to a problem known as Bartlett paradox, which arises whenever one compares models of different dimensions. The difficulty is relevant in the present case because, after deriving the rank posterior density, one may consider that models of different dimensions are being compared. The paradox is stated informally improper priors should be avoided when one computes Bayes Factors (except for parameters common to both models) as they depend on arbitrary constants (that are integrals).
More recently, [40] developed an efficient procedure to obtain the posterior distribution of the rank using a uniform proper prior over the cointegration space linearly normalized. The author derived solutions for the posterior probabilities for the null rank and for the full rank of Π. The posterior probabilities of each intermediate rank are derived from the posterior samples of the matrices that compose the long run matrix (α and β), properly normalized, under each rank and using the method proposed by [41].
This section describes how to implement the FBST to test for cointegration. We will proceed in the same spirit of Section 3, i.e., describing the steps given in Table 1 to implement the test for cointegration.
Let us begin recalling the VECM given by Equation (14):
See Equation (14), where t=1,…,T+p, in which Et∼i.i.d.Nn(0,Σ) with 0 a null vector of dimension n×1 and Ω a symmetric positive definite real matrix. Notice that these assumptions already specify the statistical model (Gaussian) and its implied likelihood. Before giving it explicitly, let us rewrite Equation (14) using matrix (15)ΔY=Z·η+E where ΔY=ΔYp+1′ΔYp+2′⋮ΔYT+p′, Z=1Dp+1′ΔYp′…ΔY2′Yp′1Dp+2′ΔYp+1′…ΔY3′Yp+1′⋮⋮⋮⋮⋮1DT+p′ΔYT−1′…ΔYT+p−1′YT+p−1′, η=c′Φ0Γ1⋮Γp−1Π and the error vector is given by E∼MNT×n(0,IT,Ω), denoting the matrix normal distribution. See Appendix B for more information on this distribution. Now the parameter vector is given by Θ=(η,Ω).
Notice that ΔY is formed by piling up T transposed vectors ΔYt, thus resulting in a matrix with T lines and n columns (n is the number of time series in vector Yt), those being also dimensions of matrix E. Matriz Z is constructed likewise—always piling up the transposed vectors—resulting in a matrix with T lines and pn+n+1 columns. Finally, matrix η has the matrices of coefficients, all piled up properly, resulting in a matrix with pn+n+1 lines and n columns.
Given the assumptions above, ΔY∼MNT×n(Z·η,IT,Ω), implying that the likelihood is L(Θ∣y)∝|Ω|−T/2exp−12·tr[Ω−1(ΔY−Z·η)′(ΔY−Z·η)] where y denotes the set of observed values of vectors Yt for t=1,…,T+p. As in Section 3, we will consider an improper prior for Θ, given by (16)h(Θ)=h(η,Ω)∝|Ω|−(n+1)/2, and our reference density, r(Θ), will be proportional to a constant, leading to a surprise function equivalent to the (full) posterior distribution. These choices correspond to steps 1 and 2 of Table 1. These modeling choices imply the following posterior (17)g(Θ∣y)∝|Ω|−(T+n+1)/2exp−12·tr[Ω−1(ΔY−Z·η)′(ΔY−Z·η)]=|Ω|−(T+n+1)/2exp−12·tr{Ω−1[S+(η−η^)′·Z′Z·(η−η^)]} where η^=(Z′Z)−1Z′ΔY and S=ΔY′ΔY−ΔY′Z(Z′Z)−1Z′ΔY.
To implement Step 3 of Table 1, we need to find the maximum a posteriori of (17) under the constraint Θ⊂ΘH, i.e., we need to maximize the posterior in ΘH. Since we are testing the rank of matrix Π, as discussed in the beginning of Section 4, it is necessary to maximize the posterior assuming the rank of Π is r, 0≤r≤n. Thanks to the modeling choices made here—Gaussian likelihood and Equation (16) as prior—our posterior is almost identical to a Gaussian likelihood, allowing us to find this maximum using a strategy similar to that proposed by [2], who derived the maximum of the (Gaussian) likelihood function assuming a reduced rank for Π. We will summarize Johansen’s algorithm, providing in Appendix C a heuristic argument of why it indeed provides the maximum value of the posterior under the assumed hypotheses.
We begin estimating a VAR(p−1) model for ΔYt with all the explanatory variables shown in (14) except for Yt−1. Using the matrix notation established above, this corresponds to estimate ΔY=Z1·η1+U, where Z1=1Dp+1′ΔYp′…ΔY2′1Dp+2′ΔYp+1′…ΔY3′⋮⋮⋮⋮1DT+p′ΔYT−1′…ΔYT+p−1′ and η1=e′τ0υ1⋮υp−1 showing that Z1 is obtained from matrix Z extracting its last n columns, exactly those corresponding to Yt−1.
We also estimate a second set of auxiliary equations, regressing Yt−1 on a vector of constants and Dt, ΔYt−1, …, ΔYt−p+1. By piling up all the (transposed) vectors Yt−1′ for t=p+1,…,T+p, we have a (T×n) matrix, denoted by Y−1. As above, these equations can be represented by Y−1=Z1·η2+V, where Y−1=Yp′Yp+1′⋮YT+p−1′ and η2=m′ν0ζ1⋮ζp−1.
Considering the OLS estimates of these sets of equations and their respective estimated residuals, we may write (18)ΔY^=Z1·η^1+U^ (19)Y^−1=Z1·η^2+V^ where η^1=(Z1′Z1)−1Z1′·ΔY, η^2=(Z1′Z1)−1Z1′·Y−1, U^ and V^ are the respective matrices of estimated residuals. Thanks to the Frisch–Waugh–Lovell theorem—see [42] theorem 3.3 or [43] Section 2.4—it is possible to show that the estimated residuals of these auxiliary regressions are related by Π in the following (20)U^=ΠV^+W^.
One can prove that the OLS estimates of Π obtained from (15) and from (20) are numerically identical, as the estimated residuals E^ and W^.
The second stage of Johansen’s algorithm requires the computation of the following sample covariance matrices of the OLS residuals obtained Σ^VV=1T·V^′V^Σ^UU=1T·U^′U^Σ^UV=1T·U^′V^Σ^VU=Σ^UV′ and, from these, we find the n eigenvalues of matrix Σ^VV−1·Σ^VU·Σ^UU−1·Σ^UV, ordering them decreasingly λ^1>λ^2>…>λ^n. The maximum value attained by the log posterior subject to the constraint that there are r (0≤r≤n) cointegration relationships is (21)ℓ*=K−(T+n+1)2·log|Σ^UU|−T+n+12·∑i=1rlog(1−λ^i), where K is a constant that depends only on T, n and y by means of the marginal distribution of the data set, y. Since ℓ* represents the maximum of the log-posterior, to obtain s*, one should take s*=exp(ℓ*), completing step 3 of Table 1.
As in Section 3, we compute ev¯ in step 4 by means of a Monte Carlo algorithm. It is easy to factor the full posterior (17) as a product of a (matrix) normal and an Inverse-Wishart, suggesting a Gibbs sampler to generate random samples from the full posterior. See Appendix B for more on the Inverse-Wishart distribution. Thus, the conditional posteriors for η and Ω are, respectively, (22)gη(η∣Ω,y)∝MNn×k(η^,(Z′Z)−1,Ω) (23)gΩ(Ω∣η,y)∝IW(Ω|S+(η−η^)′·Z′Z·(η−η^),T) where S=ΔY′ΔY−ΔY′Z(Z′Z)−1Z′ΔY, IW denotes the Inverse-Wishart, k=pn+n+1 is the number of lines of η, and η^ its OLS estimator, as above. From a Gibbs sampler set with these conditionals, we obtain a random sample from the full posterior to estimate ev¯ as the proportion of sampled vectors that generate a value for the posterior greater than s*. Finally, we obtain ev=1−ev¯ in the final step (5). The whole implementation for cointegration tests, following the assumptions made in this section, is summarized in Table 4. See Appendix A for more information on the computational resources needed to implement the steps given by Table 4.
Before presenting the results of the procedure applied to real data sets, it is important to remark one feature of the FBST applied to cointegration tests. The estimated eigenvalues of matrix Π, λ^i, correspond to the squared canonical correlations between ΔYt and Y−1 corrected for the variable in Z1 and therefore lie between 0 and 1. Therefore, (21) shows that ℓ0*≤ℓ1*≤…ℓn*, where ℓr* denotes the maximum of the posterior (14) assuming Π has rank 0≤r≤n. Therefore, one may say that the hypotheses rank(Π)=r are nested, in the sense that the respective e-values obtained by the FBST for these hypotheses are always non-decreasing ev(0)≤ev(1)≤…≤ev(n).
This nested formulation is also present in the frequentist procedure proposed by [2], based on the likelihood ratio statistics for successive ranks of Π. Thus, the FBST should be used, like the maximum eigenvalue test, in a sequential procedure to test for the number of cointegrating relationships. We will show how this should be done in presenting the applied results in the sequel.
Now we present, by means of four examples, the application of FBST as a cointegration test. In all the examples, we have adopted a Gaussian likelihood and the improper prior (16). The Gibbs sampler was implemented as described above, providing 51,000 random vectors from the posterior (17). The first 1000 samples were discarded as a burn-in sample, the remaining 50,000 being used to estimate the integral (2). The tables show the e-value computed from the FBST and the maximum eigenvalue test statistics with their respective p-values.
The labels of the modeled series refer to the electrodes on the scalp. As seen in Figure 1 and Figure 2, the series are called FP1-F7, FP1-F3, FP2-F4, and FP2-F8, where FP refers to the frontal lobes and F refers to a row of electrodes placed behind these. Even numbered electrodes are on the right side and odd numbered electrodes are on the left side. The electrodes for these four signals mirror each other on the left and right sides of the scalp. The recordings of the studied subject, an 11-year-old female, identified a seizure in the interval (measured in seconds) [2956,2996]. Therefore, like [45], we analyze the period of 41 seconds prior to the seizure—interval [2956,2996]—and the subsequent 41 seconds—interval [2996,3036]—the seizure period. In the sequel, we will refer to these as prior to seizure and during seizure, respectively. Since the sample frequency has 256 measurements per second, there are a total of 10,496 measurements for each of the four signals. Ref. [45] used 40 seconds for each period, obtaining slightly different results.
Figure 1 and Figure 2 display the estimated phases based on the original signals. The model proposed by [45] is a VAR(1), resulting in a VECM given by (24)ΔYt=c+ΠYt−1+Et.
Table 5 and Table 6 present the results that essentially lead to the same conclusions obtained by [45], even though they have based their findings on the trace test. See Table 8 of [45].
The comparison between p-values and the FBST e-values must be made carefully, the main reason being the fact that p-values are not measures supporting the null hypothesis, while e-values provide exactly such a kind of support. That being said, a possible way to compare them is by checking the decision their use recommend regarding the hypothesis being tested, i.e., to reject or not the null hypothesis.
Frequentist tests often adopt a significance level given an observed p-value, the hypothesis is rejected if the p-value is smaller or equal to the mentioned level, usually 0.1, 0.05, or 0.01. Since the cointegration ranks generate nested likelihoods, the hypotheses are tested sequentially, starting with null rank, r=0. For Table 5, adopting a 0.01 significance level, the maximum eigenvalue test would reject r=0 and r=1, and would not reject r=2. The same conclusions follow for Table 6. Thus, the recommended action is to work, for estimation purposes for instance, assuming two cointegration relationships.
The question on which threshold value to adopt for the FBST was already mentioned on Section 1.1, but it is worthwhile to underline it once more. We highly recommend a principled approach deriving the cut-off value from a loss function, which is specific for the problem at hand and the purposes of the analysis. A naive but simpler approach would be to reject the hypothesis if the e-value is smaller than 0.05 or 0.01, emulating the frequentist strategy. Even not recommending this path, since p-values are not supporting measures for the hypothesis being tested while e-values are, the researcher may numerically compare p-values and e-values in a specific scenario. If the researcher derived the p-values from a generalized likelihood ratio test, it is possible to asymptotically compare them. The relationship ev=1−Fm[Fm−h−1(1−p)], where m is the dimension of the full parameter space, h the dimension of the parameter space under the null hypothesis, Fm the chi-square distribution function with m degrees of freedom and p the corresponding p-value. See [9,12] for the proof of the asymptotic relationship between e-values and p-values.
Since the maximum eigenvalue test is derived as a likelihood ratio test, this comparison may be done for the results of all the examples presented here, and more appropriately to this example, given its sample size of 10,496 observations. Regarding Table 5 and Table 6, one could be in doubt regarding whether to reject or not the hypothesis r=1 since the e-values are larger than 0.01. However, for this model and hypothesis, the e-value corresponding to 0.01 is 0.436. Therefore, in both tables, one could reject the hypothesis and proceed to the next rank that has plenty of evidence in its favor. In conclusion, the practical decisions of both tests (FBST and maximum eigenvalue) would be the to not reject r=2.
The original data set has monthly temperatures for each coordinate (latitude and longitude) of the covered area. The mentioned series of temperatures (MAM and NDJ) are computed as seasonal averages from this monthly data set by averaging over consecutive three months. Since we have data available from January 1949 to May 2020, the time series of monthly and seasonal average surface temperatures of length 72 for each grid point.
The authors of [47] consider Yt as a two-dimensional vector, its first component being the seasonal (average) MAM temperature of year t and the second component the seasonal NDJ temperature of the previous year. They consider a VAR(2) without deterministic terms to model the series, resulting in a VECM (25)ΔYt=Γ1ΔYt−1+ΠYt−1+Et.
We have chosen five grid points corresponding to major Brazilian cities to test the cointegration hypothesis of the mentioned seasonal series. The coordinates chosen were the closest ones 23.5505° S, 46.6333° W for São Paulo; 22.9068° S, 43.1729° W for Rio de Janeiro; 19.9167° S, 43.9345° W for Belo Horizonte; 15.8267° S, 47.9218° W for Brasília and 12.9777° S, 38.5016° W for Salvador. Figure 3 and Figure 4 show the seasonal temperatures for São Paulo and Brasília, respectively, indicating that the cointegration hypothesis is plausible for both cities.
The results are shown in Table 7. Assuming a significance level of 0.01, the maximum eigenvalue test reject the null rank and do not reject r=1 for all the five cities. If we adopt the asymptotic relationship between p-values and e-values for the model under analysis, we obtain an e-value of 0.276 corresponding to a 0.01 p-value for r=0. Therefore, the FBST would also reject the null rank for all the cities. The hypothesis r=1 is not rejected since all the e-values are close to 1, once more agreeing with the maximum eigenvalue test.
One remark about Brasília seems in order. The city was built to be the federal capital, being officially inaugurated on 21 April 1960. The construction began circa 1957 and before that the site had no human occupation. The process of moving all the administration from Rio de Janeiro, the former capital, was slow and only the 1980 census detected a population over 1 million inhabitants. The present population is almost 3.2 million inhabitants living in the Federal District that includes Brasília and minor surrounding cities. Figure 4 indicates that the seasonal temperatures began to rise exactly after 1980.
The authors used the logarithms of the series of M1 monetary aggregate, inflation rate, real income, and the primary interest rate set by the Bank of Finland to model the money demand, which, in theory, follows a long-term relationship. The sample has 106 quarterly observations of the mentioned variables, starting at the second trimester of 1958 and finishing in the third trimester of 1984. The chosen model was a VAR(2) with unrestricted constant, meaning that the series in Yt have one unit root with drift vector c and the cointegrating relations may have a non-zero mean. For more information about how to specify deterministic terms in a VAR see [48], chapter 6. Seasonal dummies for the first three quarters of the year were also considered in the model chosen by [2]. Writing the model in the error correction form, we (26)ΔYt=c+Φ0,1D1t+Φ0,2D2t+Φ0,3D3t+Γ1ΔYt−1+ΠYt−1+Et. where Π=Φ1+Φ2−I3, Γ1=−Φ2, c is a vector with constants and Dit denote the seasonal dummies for trimester i=1,2,3. The results are displayed in Table 8.
In [2], the authors concluded that there is, at least, two cointegration vectors, a conclusion that follows if one adopts a 0.01 significance level, for instance. Using the asymptotic relationship between p-values and e-values for Equation (26), we obtain, for r=0, an e-value of 0.998, and, for r=1, an e-value of 0.999, corresponding to a 0.01 p-value. These apparently discrepant values for the e-values are due to the high dimensions of the unrestricted (m=58) and under H0 (h=42 for r=0 and h=43 for r=1) parameter spaces. Therefore, under this criterion, the FBST also rejects the null rank and r=1 (since 0.132 < 0.998 and 0.994 < 0.999, respectively) and does not reject r=2, recommending the same action as the maximum eigenvalue test.
Table 9 shows that the maximum eigenvalue test rejects r=0 and does not reject r=1 at a 0.05 significance level. Once more adopting the asymptotic relationship between p-values and e-values for the chosen model, we obtain, for r=0, an e-value of 0.247 corresponding to a 0.01 p-value. Thus, under this criterion, the FBST also rejects the null rank and does not reject r=1.
In the past few decades, the econometric literature introduced statistical tests to identify unit roots and cointegration relationships in time series. The Bayesian approach applied to these topics advanced considerably after the 1990s, developing interesting alternatives, mostly for unit root testing. The (parametric) frequentist tests mentioned here may not be suitable since these procedures rely on the distribution of the test statistic—usually assuming the hypothesis being tested is true—which depend on a particular a statistical model, usually Gaussian. When the distributions of such statistics cannot be obtained, the procedure is saved by asymptotic results. If the researcher considers different statistical models and the available sample is small, the results of the tests may be quite misleading.
The present work reviewed a simple and powerful Bayesian procedure that can be applied to both unit root and cointegration testing. We have also shown that the FBST works considerably well even when one uses improper priors, a choice that may preclude the derivation of Bayes Factors, a standard Bayesian procedure in hypotheses testing.
A long series of articles provided in [7] and the references therein, has showed the versatility and properties of FBST, such a. the e-value derivation and computation are straightforward from its general definition; b. it uses absolutely no artificial restrictions like a distinct probability measure on the hypothesis set, induced by some specific parametrization; c. it is in strict compliance with the likelihood principle; d. it can conduct the test with any prior distribution; e. it does not need closed conjectures concerning error distributions, even for small samples; f. it is an exact procedure, since it does rely on asymptotic assumptions; and g. it is invariant with respect to the null hypothesis parametrization and with respect to the parameter space parametrization. See [9], p. 253 for this property.
To proceed with this research agenda, it would be interesting to perform more simulation studies with the FBST applied to unit root testing for a larger group of parametric and semi-parametric models (likelihoods). Another possibility is to include moving average terms in the data generating processes and work with Gaussian and non-Gaussian ARMA models. Notice that, given the points made above, these extensions would not impose major problems to the FBST as they would to the frequentist procedures. Regarding cointegration, the same extensions may be studied in future works, although the adoption of statistical models outside the Gaussian family would require further efforts to numerically implement the FBST. We shall also investigate the effect of the prior choice in the estimates of cointegration relations, especially for small samples.