Bayesian multivariate Beveridge–Nelson decomposition of I(1) and I(2) series with cointegration

3.1 VAR model

Let for d = 1, 2, x t , d be an N _d -variate I(d) sequence. Let N ≔ N ₁ + N ₂. Let for all t, x t ≔ x t , 1 ′ , x t , 2 ′ ′ , y _t,1 ≔ x _t,1, y _t,2 ≔ Δ x _t,2, and y t ≔ y t , 1 ′ , y t , 2 ′ ′ , where the prime (′) denotes the transpose of a vector or a matrix, so that { y _t } is I(1). Assume also that { y _t } is CI(1,1) with cointegrating rank r.^[9] Let for d = 1, 2, μ _d ≔ E(Δ y _t,d ). Let μ ≔ μ 1 ′ , μ 2 ′ ′ . Let { y t * } be such that for all t,

where α is an intercept vector. Assume a VAR(p + 1) model for { y t * } such that for all t,^[10]

where Π (L) is a (p + 1)th-order polynomial matrix in the lag operator L and WN( Σ ) means ‘(multivariate) white noise sequence with variance–covariance matrix Σ ’.

3.2 VECM representation

Write

where Φ (L) ≔ ( Π (L) − Π (1)L)/(1 − L). Then we have a VECM of order p for { Δ y t * } such that for all t,

where Λ , Γ ∈ R N × r . Since { y t * } is CI(1,1), the roots of det( Π (z)) = 0 must lie on or outside the unit circle. This requirement gives an implicit restriction on the VECM parameters ( Φ (.), Λ , Γ ).

We can write for all t,

where β ≔ Γ ′ α and δ ≔ Γ ′ μ . Though slightly different, the last expression is essentially a steady state VECM suggested by Villani (2009, p. 633). We have for all t,

which help us to specify an informative prior on ( μ , β , δ ). Thus a steady state VECM is useful for Bayesian analysis.

3.3 State space representation

Assume that p ≥ 1. We have for all t,

or

Let s _t be a state vector such that for all t,

which is I(0) and observable given the model parameters. A state space representation of the steady state VECM is for all t,

where

with O _m×n denoting an m × n null matrix. Note that { s _t } is I(0) if and only if the eigenvalues of A lie inside the unit circle.

We have for all t, for h ≥ 0,

or

where

3.4 Multivariate B–N decomposition

Introducing cointegration changes the state space model, but the formulae for the multivariate B–N decomposition of I(1) and I(2) series given by Murasawa (2015, Theorem 1) remain almost unchanged. Let x t * and c _t be the B–N permanent and transitory components in x _t , respectively.

Let

Then for all t,

where W depends only on the VECM coefficients and { s _t } is observable given the model parameters. This observation is useful for Bayesian analysis of { c _t }.

4.1 Conditional likelihood function

Assume Gaussian innovations for Bayesian analysis, and write the VECM as for all t,

where I N N 0 N , P − 1 means ‘independent N-variate normal distributions with mean vector 0 _N (null vector) and variance–covariance matrix P ⁻¹’. Let ψ ≔ ( β ′, μ ′)′, Φ ≔ [ Φ ₁, …, Φ _p ], and Y ≔ [ y ₀, …, y _T ].^[11] By the prediction error decomposition, the joint pdf of Y is^[12]

Our Bayesian analysis relies on ∏ t = p + 1 T p ( Δ y t | s t − 1 , ψ , Φ , P , Λ , Γ ) , which is equivalent to the conditional likelihood function of ( ψ , Φ , P , Λ , Γ ) given s _p .

4.2 Identification

Identification of ( Λ , Γ ) from Π (1) requires some restrictions, though such identification is unnecessary for the multivariate B–N decomposition. Let

Then Λ _* Γ _*′ = ΛΓ ′ and Γ _*′ Γ _* = I _r . These restrictions are common, but do not identify the signs of Λ _* and Γ _*.

Alternatively, we can apply linear normalization. Write Γ ≔ [ Γ ₁′, Γ ₂′]′, where Γ ₁ is r × r and Γ ₂ is (N − r) × r. Let

Then we can identify Λ ̄ , Γ ̄ . Let β ̄ ≔ Γ ̄ ′ α = Γ 1 − 1 ′ β and ψ ̄ ≔ β ̄ ′ , μ ′ ′ correspondingly. Note that α is not identifiable from the VECM.

4.3 Prior

4.4 Posterior simulation

We simulate p( ψ , Φ , P , Λ , Γ , ν| Y ) by a Gibbs sampler consisting of five blocks:

1. Draw ψ from p( ψ | Φ , P , Λ , Γ , ν, Y ) = p( ψ | Φ , P , Λ , Γ , Y ).

2. Draw ( Φ , P ) from p*( Φ , P | ψ , Λ , Γ , ν, Y).

3. Draw Λ from p*( Λ | ψ , Φ , P , Γ , ν, Y ) = p*( Λ | ψ , Φ , P , Γ , Y ).

4. Draw Γ from p*( Γ | ψ , Φ , P , Λ , ν, Y ) = p*( Γ | ψ , Φ , P , Λ , Y ). Accept the draw if ( Φ , Λ , Γ ) ∈ S; otherwise go back to step 2 and draw another ( Φ , P , Λ , Γ ).

5. Draw ν from p(ν| ψ , Φ , P , Λ , Γ , Y ) = p(ν| Φ , P ).

The first block builds on Villani (2009); the second block is standard; the third and fourth blocks follow the parameter-augmented Gibbs sampler proposed by Koop, León-González, and Strachan (2010); the fifth block is standard. See the Appendix for the details of each block.

4.5 Bayes factor

We use the Bayes factor for Bayesian model selection. When choosing between nested models with certain priors, the Savage–Dickey (S–D) density ratio gives the Bayes factor without estimating the marginal likelihoods; see Wagenmakers et al. (2010) for a tutorial on the S–D method.

We choose the cointegrating rank r. Consider comparing the following two models (hypotheses):

Koop, León-González, and Strachan (2008, pp. 451–452) note that the problem is the same as comparing the following two nested models:

The restriction ( Φ , Λ , Γ ) ∈ S is necessary only for the B–N decomposition, and unnecessary for forecasting in general; hence we ignore this restriction for the moment, so that the priors on Λ and Γ become independent.^[13] Then the S–D density ratio for H ₀ versus H _r is

The prior gives the denominator directly. For the numerator, we have

Let { ψ l , Φ l , P l , Γ l } l = 1 L be posterior draws. Let

See Appendix A.4 for the conditional posterior p( Λ | ψ , Φ , P , Γ , Y ). An estimator of the S–D density ratio for H ₀ versus H _r is

5.1 Data

We consider joint estimation of the natural rates (or their permanent components) and gaps of the following four macroeconomic variables in the US:

1. Output Let Y _t be output. Assume that { ln Y _t } is I(2), so that {Δln Y _t } is I(1).

2. Inflation rate Let P _t be the price level and π _t ≔ ln(P _t /P _t−1) be the inflation rate. Assume that {π _t } is I(1).

3. Interest rate Let I _t be the 3-month nominal interest rate (annual rate in per cent), i _t ≔ ln(1 + I _t /400), r _t ≔ i _t − E _t (π _t+1) be the ex ante real interest rate, and r ̂ t ≔ i t − π t + 1 be the ex post real interest rate. Assume that {r _t } is I(1).^[14]

4. Unemployment rate Let L _t be the labor force, E _t be employment, and U _t ≔ − ln(E _t /L _t ) be the unemployment rate. Assume that {U _t } is I(1).

Table 1 describes the data, which are available from FRED (Federal Reserve Economic Data). When monthly series are available, i.e., except for real GDP, we take the 3-month arithmetic means of monthly series each quarter to obtain quarterly series, from which we construct the quarterly inflation, interest, and unemployment rates as defined above.

Let for all t,

The sample period of { y _t } is 1948Q1–2018Q4 (284 observations). Table 2 shows summary statistics of { y _t } and {Δ y _t } multiplied by 100. Note that {π _t }, { r ̂ t } , and {Δln Y _t } are quarterly rates (not annualized).

5.2 Preliminary analyses

As preliminary analyses, we check if { y _t } is indeed I(1), though our method is valid even if some series in { y _t } are I(0). For each I(0) series, we treat it as an I(1) series cointegrated with itself, and increase the cointegrating rank by 1 as ‘pseudo’ cointegration; see Fisher, Huh, and Pagan (2016, sec. 2.2).

Table 3 shows the results of the ADF and ADF-GLS tests for unit root, with or without constant and/or trend terms in the ADF regression.^[15] The results depend on the number of lags included in the ADF regression, since the ADF test suffers from size distortion with short lags and low power with long lags. The ADF-GLS test mitigates the problem except when there is no constant nor trend term in the ADF regression, in which case the ADF test is asymptotically point optimal. The level 0.05 ADF-GLS test rejects H 0 : y t , i ∼ I ( 1 ) in favor of H 1 : y t , i ∼ I ( 0 ) for {Δln Y _t } (and possibly for {π _t }, whose results are mixed). Hence these unit root tests do not support our assumption that { y _t } is I(1).

Table 4 shows the results of the KPSS stationarity tests, with or without a trend term. The results depend on the lag truncation parameter for the Newey–West estimator of the long-run error variance. With a trend term, the level 0.05 KPSS test rejects H ₀ : {y _t,i } ∼ I(0) in favor of H ₁ : {y _t,i } ∼ I(1) except for {Δln Y _t }. With no trend term, the test rejects H ₀ : {y _t,i } ∼ I(0) in favor of H ₁ : {y _t,i } ∼ I(1) except for { r ̂ t } . Though the results for { r ̂ t } and {Δln Y _t } are mixed, these stationarity tests support our assumption that { y _t } is I(1).

Overall, these unit root and stationarity tests confirm that {Δ y _t } is I(0), but are inconclusive if each component of { y _t } is I(1). We can still proceed with our assumption that { y _t } is I(1), allowing for ‘pseudo’ cointegration.

As preliminary analyses, we also check if { y _t } is CI(1,1). Table 5 shows the results of the Engle–Granger cointegration tests, with or without a trend term in the cointegrating regression and for each choice of the left-hand side (LHS) variable. The results depend on the number of lags included in the ADF regression for the residual series, but the level 0.05 test rejects the null of no cointegration in favor of the alternative of cointegration only when we put Δln Y _t on the LHS of the cointegrating regression.

Table 6 shows the results of Johansen’s cointegration tests, with unrestricted constant and restricted trend terms in the VECM. The results depend on the number of lags included in the VECM, but overall, the tests suggest that the cointegrating rank is 2, possibly including ‘pseudo’ cointegration.

Since the results of these cointegration tests are mixed, we estimate VECMs with different cointegrating ranks, and use the Bayes factor to choose the cointegrating rank.

5.3 Model specification

For our data, N ≔ 4. To select p, we fit VAR models to { y _t } up to VAR(8), i.e., p = 7, and check model selection criteria. The common estimation period is 1950Q1–2018Q4. Table 7 summarizes the results of lag order selection.^[16] The level 0.05 LR test fails to reject H ₀ : { y _t } ∼ VAR(7) against H ₁ : { y _t } ∼ VAR(8) and AIC selects VAR(7), but BIC and HQC select much smaller models. Since a high-order VAR model covers low-order VAR models as special cases, to be conservative, we choose a VAR(8) model for { y _t }, i.e., p = 7, and impose a shrinkage prior on the VAR coefficients in the VECM.

For the prior on α , we set α ₀ ≔ 0 _N and Q _{0, α} ≔ I _N . For the prior on μ , we set μ 0 ≔ μ ̂ , where μ ̂ is the sample mean of {Δ y _t }, and Q _{0, μ} ≔ I _N .

Following Kadiyala and Karlsson (1997) and Bańbura, Giannone, and Reichlin (2010), we set

where for i = 1, …, N, s i 2 is an estimate of var(u _t,i ) based on the univariate AR(p + 1) model with constant and trend terms for {y _t,i }.

For the prior on Λ , we set η ₀ ≔ 1 and G ₀ ≔ I _N . For the prior on Γ , we set τ ₀ ≔ 1, which gives the flat prior on the cointegrating space.

For the prior on ν, we set A ₀ ≔ 1 and B ₀ ≔ 1, i.e., ν ∼ χ ²(1); hence the tightness hyperparameter on the VAR coefficients tends to be small, implying potentially mild shrinkage toward M ₀ ≔ O _N×pN .

Overall, our priors are weakly informative in the sense of Gelman et al. (2014, p. 55).

5.4 Bayesian computation

We run our Gibbs sampler on R 4.0.5 developed by R Core Team (2021). We use the ML estimate of ( Φ , P , Λ , Γ ) for their initial values.^[17] With poor initial values, the restriction that the eigenvalues of A lie inside the unit circle does not hold, and the iteration cannot start. Hence the choice of initial values seems important when one applies the B–N decomposition. Once the iteration starts, the restriction rarely binds for our sample.

To check convergence of the Markov chain generated by our Gibbs sampler to its stationary distribution, we perform convergence diagnoses discussed in Robert and Casella (2009, ch. 8) and available in the coda package for R. Given the diagnoses, we discard the initial 1000 draws, and use the next 4000 draws for our posterior inference. At significance level 0.05, Geweke’s Z-score rejects equality of the means of the first 10% and the last 50% of the posterior sample (the default choice in the coda package) for 6 out of 141 identified parameters (applying linear normalization to Γ ). Since the rejection rate is close to the significance level, 1000 draws are enough as burn-in. The minimum effective sample size among the 141 identified parameters is 238. Gelman et al. (2014, p. 267) write, ‘…100 independent draws are enough for reasonable posterior summaries.’ Hence 4000 draws are enough for posterior inference.

To select the cointegrating rank, we set p ≔ 7, assume the above priors, and compute the S–D density ratios for r = 0 versus r = 1, 2, 3, respectively, from which we derive the posterior probabilities of r = 0, 1, 2, 3, assuming the flat prior on r. We find that the posterior probability of r = 2 is often numerically 1, which is consistent with the results of Johansen’s cointegration tests.^[18] Thus we set r ≔ 2 in the following analysis.

5.5 Empirical results

Figure 2 plots the actual rates and our point estimates (posterior medians) of the natural rates (or their permanent components) of the four variables. For ease of comparison of the actual and natural rates, we omit error bands for the natural rates, which are identical to those for the gaps. Figure 3 plots our point estimates of the gaps and their 95% error bands. Since a VECM has both the level and difference series, to avoid confusion in estimation, we do not annualize quarterly rates. To compare our results with previous works that use annualized quarterly rates, one must multiply our inflation and interest rates by 4 or 400.

Figure 2:

Actual rates (solid red) and the posterior medians of the natural rates (dashed green). The shaded areas are the NBER recessions.

Figure 3:

Posterior medians of the gaps and their 95% error bands (posterior 0.025- and 0.975-quantiles). The shaded areas are the NBER recessions.

Our estimate of the natural rate of inflation is smoother than typical univariate estimates of trend inflation in the CPI; e.g., Faust and Wright (2013, p. 22). Our estimate looks close to a recent estimate of trend inflation in the US CPI in Morley, Piger, and Rasche (2015, p. 894) based on a bivariate UC model for the inflation and unemployment rates, which assumes independent shocks to the trend and gap components and allows for structural breaks in the variances of these shocks. Our estimate is more volatile, however, than a recent estimate of trend inflation in the PCE price index by Chan, Clark, and Koop (2018, p. 21) that uses information in survey inflation expectations. Hwu and Kim (2019) estimate trend inflation in the PCE price index using a univariate UC model with correlated shocks to the trend and gap components, Markov-switching gap persistence, and Markov-switching association between inflation and inflation uncertainty. Without the last feature, their estimate looks close to ours; see Hwu and Kim (2019, p. 2314).

Our estimate of the natural rate of interest is more volatile than a recent estimate by Del Negro et al. (2017, p. 237) based on a VAR model with common trends (or a multivariate UC model with independent shocks to the trend and gap components and a factor structure for the trend components) for short- and long-term interest rates, inflation and its survey expectations, and some other variables. Their estimate is smooth partly because they impose tight priors on the variances of the shocks to the trend components. With the ‘loosest possible prior’, their estimate is as volatile as ours; see Del Negro et al. (2017, p. 272).^[19] Despite different definitions of the natural rate, their estimate based on a DSGE model looks close to ours; see Del Negro et al. (2017, p. 237).^[20] Estimates by Laubach and Williams (2016, p. 60), Holston, Laubach, and Williams (2017, p. S61), and Lewis and Vazquez-Grande (2019) are close to trend output growth by construction; hence their estimates are quite different from those by Del Negro et al. (2017) and ours, especially before 1980.^[21]

Our estimate of the natural rate of unemployment looks more volatile than a recent estimate by Morley, Piger, and Rasche (2015, p. 898), obtained as a by-product of estimation of trend inflation. A possible reason for the difference is that they assume independent shocks to the trend and gap components, which may not hold in practice. The two estimates may coincide if one allows for dependence between the shocks, as Morley, Nelson, and Zivot (2003) show for the univariate trend-cycle decomposition. Despite the difference in volatility, our estimate of the unemployment rate gap looks close to that by Morley, Piger, and Rasche (2015, p. 901) in terms of the sign and magnitude. Recent estimates of the natural rate and gap by Crump et al. (2019, pp. 176, 180) are as smooth as those by Morley, Piger, and Rasche (2015), perhaps partly because they also assume independent shocks to the trend and gap components.^[22]

Our estimate of the output gap is at most about ±5% of the output level, and looks close to a recent estimate by Morley and Wong (2020, p. 9) based on a large Bayesian VAR(4) model with 23 variables. Our estimate also looks close to their estimate based on a Bayesian VAR(4) model with four variables using output growth, the unemployment rate, the CPI inflation rate, and the growth rate of industrial production, assuming that they are all I(0); see Morley and Wong (2020, p. 8). Though output growth may or may not be I(0) in the US, it may be clearly I(1) in some countries or regions, in which case their method may give an unreasonable estimate of the output gap with a strong upward or downward trend, as Murasawa (2015) shows for the Japanese data.

Figure 4 plots the posterior probability of positive gap for the four variables. This probability index is useful if the sign of the gap is of interest. Even if the 95% error band for the gap covers 0, the posterior probability of positive gap may be close to 0.025 or 0.975. Indeed, the probability index is often below 0.25 or above 0.75; hence we are often quite sure about the sign of the gap. Moreover, Figure 4 shows the relation between the gaps more clearly than Figure 3.

Figure 4:

Posterior probability of positive gap. The shaded areas are the NBER recessions.

Figure 5 shows the relation between the gaps more directly. The left panels are the scatter plots of the posterior medians of the gaps in each quarter. The right panels are the posterior pdfs of the correlation coefficients between the gaps. We see that the Phillips curves and Okun’s law hold between the gaps, though we do not impose such relations. Thus our estimates of the gaps seem mutually consistent from a macroeconomic point of view.

Figure 5:

Phillips curves and Okun’s law. The left panels are the scatter plots of the posterior medians of the gaps in each quarter. The right panels are the posterior pdfs of the correlation coefficients between the gaps.

Figure 6 shows the relation between the trend output growth rate Δ ⁡ ln Y t + 1 * and the natural rate of interest r t * , i.e., the dynamic IS curve; see Eq. (1). The left panel is the scatter plot of the posterior medians of Δ ⁡ ln Y t + 1 * and r t * in each quarter. The right panel is the posterior pdf of the slope coefficient, i.e., 1/σ in Eq. (1), obtained by regressing Δ ⁡ ln Y t + 1 * on r t * for each posterior draw of { Δ ⁡ ln Y t * , r t * } . On average, the slope is positive as expected, but the evidence is weak. The result is similar to Hamilton et al. (2016) and Lunsford and West (2019), who study the determinants of r t * and find positive but low and insignificant long-run correlations between Δln Y _t and r _t in the postwar US.

Figure 6:

Dynamic IS curve. The left panel is the scatter plot of the posterior medians of Δ ⁡ ln Y t + 1 * and r t * in each quarter. The right panel is the posterior pdf of the slope coefficient obtained by regressing Δ ⁡ ln Y t + 1 * on r t * for each posterior draw of { Δ ⁡ ln Y t * , r t * } .

5.6 Comparison of alternative model specifications

Figure 7 compares point estimates of the gaps under three alternative assumptions, i.e., I(1) log output, I(2) log output with no cointegration, and I(2) log output with cointegration.^[23] For ease of comparison, we omit error bands here.

Figure 7:

Posterior medians of the gaps assuming I(1) log output (solid red), I(2) log output with no cointegration (dashed green), and I(2) log output with cointegration (longdash blue). The shaded areas are the NBER recessions.

The result assuming I(1) log output is similar to that in Murasawa (2014, p. 513), who assumes a VAR(12) model with no constant term for the centered stationarized series, and chooses the tightness hyperparameter by the empirical Bayes method. In contrast to the result for Japan in Figure 1, for the US, the multivariate B–N decomposition assuming I(1) log output gives a ‘reasonable’ estimate of the output gap that fluctuates around 0. However, the output gap is persistently positive since 2010, which may not be sensible.

Assuming I(2) log output changes the estimate of the output gap. However, it hardly changes the estimates of other gaps, similar to the result for Japan in Murasawa (2015), which makes sense because the B–N transitory components are ‘forecastable movements’, and whether log output is I(1) or I(2), the two VAR models use the information in log output for forecasting other variables in almost the same way.^[24] The output gap now keeps fluctuating around 0 since 2010, perhaps because I(2) log output captures a possible structural break in the mean output growth rate after the 2008 Great Recession.

Allowing for cointegration gives much bigger estimates of the gaps for all variables. The result follows partly because a VECM uses more predictors, i.e., the error correction terms, than a VAR model; see Evans and Reichlin (1994). However, the result also implies that for all variables, the coefficients on the error correction terms are large enough, if not significantly different from 0, to obtain such different estimates of the gaps. These bigger estimates of the gaps, or the associated estimates of the natural rates, are often close to other recent estimates that focus on a particular natural rate or gap, as already noted.