Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,261 characters · 20 sections · 0 citation commands
{.26in} \thispagestyle{empty}
\doublespacing
The causal-noncausal (mixed) Vector Autoregressive processes have nonlinear dynamics with locally explosive patterns such as local trends, spikes and bubbles [Gourieroux and Jasiak (2017), Davis and Song (2020)]. In applied research, these strictly stationary models can replicate the behavior of commodity prices, including oil prices [Lof and Nyberg (2017), Cubadda et al. (2023), Blasques et al. (2025)], cryptocurrency rates [Hencic and Gourieroux (2019), Hall and Jasiak (2024), Cavaliere et al. (2020)], financial indexes, including S&P500 and NASDAQ [Gourieroux, Zakoian (2017), Freis (2022)], climate risk on El Nino and La Nina [De Truchis, Fries and Thomas (2024)], and green stocks prices [Hecq et al (2024)]. Specifically, a mixed VAR(1) process ($Y_t$) of dimension $m$ is the strictly stationary solution of $Y_t = \Phi Y_{t-1} + \varepsilon_t$, where the errors $\varepsilon_t$ are assumed to be non-Gaussian, independent and identically distributed (i.i.d.), and the eigenvalues of of autoregressive matrix $\Phi$ are of modulus either smaller or greater than 1, and are referred to as the causal and non-causal eigenvalues, respectively. In the presence of non-causal eigenvalues, the errors $\varepsilon_t$ are correlated with the lagged values of $Y_t$ and cannot be interpreted as the (causal) innovations of the process. We show that the mixed VAR model has a causal, i.e. past-dependent representation, which is a multivariate non-linear autoregression $Y_t = a(Y_{t-1}; v_t)$, where $a$ is a non-linear function, and the errors $v_t$ are i.i.d. and independent of $Y_{t-1}$, satisfying the definition of {\it causal} innovations.
The initial papers on mixed VARs were focused on the identification of noncausal dynamics and parameter estimation. Davis and Song (2020) introduced the maximum likelihood estimator under a parametric assumption on the distribution of error $\epsilon_t$. Gourieroux and Jasiak (2017), (2023) proposed a semi-parametric one-step Generalized Covariance (GCov) estimator which is consistent and semi-parametrically efficient. This paper is focused on the Markov property of mixed VAR processes, which leads to nonlinear predictive density formulas and nonlinear impulse response functions (IRF) for these processes. We derive the closed-form expressions of forward and backward predictive densities, and build prediction intervals (sets) that account for estimation risk and are adjustable to the desired conditional coverage. We also discuss the identification of the aforementioned causal innovations and define the nonlinear IRF for tracing out the effects of shocks to the locally explosive component of the process. The aim of this paper is to provide a toolbox of inference methods for mixed VAR models, equivalent to those available for the traditional causal VARs, to facilitate the use of these models in applied research in finance and macroeconomics.
The closed-form predictive density for forecasting the mixed VAR processes is new and derived in a semi-parametric framework. This closed-form representation was deemed infeasible until recently\footnote{"The predictive density is generally not available under closed form", Fries and Zakoian (2019)}. It extends to a general framework earlier results on a mixed VAR model with a multiplicative form of autoregressive matrix, also called the multivariate Mixed Autoregressive (MAR) model [see Nyberg and Saikkonen (2014), Gourieroux and Jasiak (2016), Lanne and Luoto (2016)]. The two types of models are not equivalent, except for pure causal and noncausal processes, because the multiplicative form assumption is strongly constraining and may not exist for all mixed VAR processes. In particular, it does not exist for a bivariate, mixed VAR(1) model with at least one causal and one non-causal eigenvalues [see Appendix A.2]. The availability of a closed-form predictive density for the general VAR model considered in this paper eliminates the need for assuming a restrictive multiplicative form of autoregressive matrix and using computationally intense forecasting methods based on numerical approximations.
From the quantiles of the predictive densities, we infer the point forecasts and prediction intervals. In addition, we derive a new backcasting method for noncausal processes. The backcasting algorithms are readily available for time-reversible linear Gaussian time series, but need to be developed specifically for the time-irreversible non-Gaussian mixed VAR processes for the treatment of missing data, or implementation of backpropagation algorithms in applied research [Twumasi and Twumasi (2022)]. In our paper, we use the backcasting algorithm to build a new bootstrapped prediction interval adjustable to a desired coverage level, and we show how to account for the estimation risk by introducing the confidence set for prediction intervals. This approach is applicable to any nonlinear dynamic model where the forecast error is not a sum of the unknown future innovation and estimation error, and the process is back-castable. The adjustment of the empirical coverage of a prediction interval to the desired one is beneficial when the multivariate predictive density is estimated from data with a limited number of tail observations, or evaluated over insufficiently support points because of a computational burden.
Another contribution of this paper is in defining the nonlinear causal innovations for mixed VAR models and examining their identification issues. We describe the nonlinear IRF analysis of mixed VARs, with a shock applied to a locally explosive component of the process, and performed on- and off-bubble. The nonlinear causal innovation for the mixed VAR model is based on Gourieroux and Jasiak (2005), which has recently received attention in the context of nonlinear impulse response functions in macroeconomics [see, e.g. Gonzalves et al.(2021), Gourieroux and Lee (2025)].
The paper is organized as follows. Section 2 describes the causal-noncausal VAR model and its state-space representation. Sections 3, 4 and 5 contain the new results. Section 3 proves the Markov property and derives the closed-form formulas of multivariate predictive density for forecasting and backcasting. Section 4 introduces the inference on the random set of prediction intervals. Section 5 defines the nonlinear causal innovations and discusses their identification and filtering. A simulation study and an empirical application to a bivariate process of oil prices and real US GDP rates are presented in Section 6. Section 7 concludes. The technical results are given in Appendices A.1-A.2. Appendix A.1 contains the proof of the forward and backward predictive density formula. Appendix A.2 explains the constraints induced by the multiplicative representation of a causal-noncausal VAR model. Online Appendices B, C, D and E present the closed-form expression of a kernel-based semi-nonparametric estimator of predictive density, give additional results on simulations and the empirical application of the mixed VAR model to oil prices and real GDP rates.
\setcounter{equation}{0}
\setcounter{equation}{0} This Section reviews the state space representation of the causal-noncausal (mixed) VAR(p) model studied in Gourieroux and Jasiak (2016),(2017), and Davis and Song (2020).
The multivariate causal-noncausal VAR(p), referred to as the mixed VAR process henceforth, is defined by:
where $Y_t$ is a vector of size $m$, $\Phi_j, \, j=1,...,p$, are matrices of autoregressive coefficients of dimension $m \times m$ and $(\varepsilon_t)$ is a sequence of errors, which are serially i.i.d. random vectors of dimension $m$ with mean zero, and common joint density $g$. The errors $(\varepsilon_t)$ are assumed to have a non-Gaussian distribution and are not assumed independent of past $Y$'s. This implies that $(\varepsilon_t)$ cannot be interpreted as the (causal) innovations.
For the existence of a unique strictly stationary solution to the VAR model (2.1), we assume that the roots of the characteristic equation of the autoregressive polynomial matrix $det(Id - \Phi_1 \lambda - \cdots \Phi_p \lambda^p) = 0$ are of modulus either strictly greater, or strictly smaller than one, i.e. are either outside, or inside the unit circle, and that the joint density $g$ has a uniform tail index $\alpha$. The strictly stationary solution $(Y_t)$ to model (2.1) can be written as an infinite two-sided moving average (MA($\infty$)) in errors $\varepsilon_t$:
This is a linear time series, according to the terminology of Rosenblatt (2012). The autoregressive matrices $\Phi_1,...,\Phi_p$ and the matrices of coefficients $C_j$ on the past and future terms of this MA representation are uniquely defined when $\varepsilon_t$ is non-Gaussian, which is an identifying assumption. Process $Y_t$ is said to be causal in $\varepsilon_t$, if $Y_t = \sum_{j = 0}^{+ \infty} C_j \varepsilon_{t-j}$, noncausal in $\varepsilon_t$, if $Y_t = \sum_{j = - \infty}^{-1} C_j \varepsilon_{t-j} = \sum_{j = 1}^{+ \infty} C_{-j} \varepsilon_{t+j}$, or mixed, otherwise.
In the presence of a noncausal component, the assumption of strict stationarity of ($Y_t$) implies that, in calendar time, $(Y_t)$ has a nonlinear dynamic and past-dependent conditional heteroscedasticity. Then, the process $(Y_t)$ can be characterized by a rather complicated conditional distribution of $Y_{t+h}$ given $\underline{Y_t} = (Y_t, Y_{t-1},...)$ for $h=1,2,..$, that leads us to the nonlinear out-of-sample (oos) predictive distribution described in Section 3. The nonlinearity in the path of $(Y_t)$ manifests itself through local trends, spikes and bubbles similar to those observed in the time series of commodity (oil) prices, exchange rates, or cryptocurrency prices [Gourieroux and Zakoian (2017), Gourieroux and Jasiak (2017), Gourieroux, Jasiak and Tong (2021)].
Let us recall the state-space representation of the mixed VAR process with the latent causal and noncausal components of the observed process as state variables.
a) {\it The mixed VAR(1) representation of a mixed VAR(p) process}
As it is commonly done in the literature on multivariate causal autoregressive processes, model (2.1) can be rewritten as a $n=mp$ multivariate mixed VAR(1) model, by stacking the present and lagged values of process $(Y_t)$:
The eigenvalues of the autoregressive matrix $\Psi$ are reciprocals of the roots of the characteristic equation (2.4).
b) {\it The Change of Basis}
Matrix $\Psi$ has a real Jordan representation: $ \Psi = A \left(
\right) A^{-1}, $ where $J_1$ (resp. $J_2$) are real $(n_1 \times n_1)$ (resp. $(n_2 \times n_2)$) matrices where $n_2 = n-n_1$, with all eigenvalues of modulus strictly less than 1 [resp. strictly larger than 1], and $A$ is a ($n \times n$) invertible matrix [see, Perko (2001), Gourieroux and Jasiak (2017), Section 5.2, for real Jordan representations]. Then, the above equation can be rewritten after the change of basis $A^{-1}$ as:
Let us introduce a block decomposition of $A^{-1} $: $ A^{-1} \equiv \left(
\right), $ where $A^{1}$ is of dimension $(n_1 \times n)$, and define the transformed variables:
This leads us to two sets of state components of process ($Y_t$) such that:
Hence, the state-space representation of process ($Y_t$) consists of: {\bf state equations} (2.5) representing the causal and non-causal dynamics of state variables $Z_t$ and {\bf measurement equations} for $Y_t$ obtained by solving: $$\left(
\right) = A Z_t.$$ \noindent These measurement equations are deterministic. Therefore, the filtrations generated by processes $(Y_t)$ and $(Z_t)$ are identical.
Since the Jordan representation of matrix $\Psi$ is not unique, the state-space representation is not unique either. It depends on the choice of state variables, i.e. the latent factors $Z_1$ and $Z_2$, up to linear invertible transformations.
c) {\it State-Specific Linear Errors}
The first set of state equations in (2.5) defines the causal VAR(1) process ($Z_{1,t}$), which has the causal MA$(\infty)$ representation:
where $\eta_{1,t}$ is a function of $\varepsilon_{t}, \varepsilon_{t-1}...$. The second set of state equations in (2.5) needs to be inverted to obtain a MA representation in matrices with eigenvalues of modulus strictly less than 1. We get:
We observe that $(Z_{2,t})$ is a noncausal process with a one-sided moving average representation in future values $\varepsilon_{t+1}, \varepsilon_{t+2},...,$.
Let us now discuss the state-specific linear errors $\eta_{1,t}, \eta_{2,t}$. From equations (2.8), (2.9) and the stationarity conditions, it follows that:
where $\underline{Z}_{1,t-1} = (Z_{1,t-1}, Z_{1,t-2},...$) and $\bar{Z}_{2,t+1} = (Z_{2,t+1}, Z_{2,t+2},...$). Therefore, $\eta_{1,t}$ (resp. $\eta_{2,t}$) can be interpreted as the causal linear innovation of $Z_{1,t}$ based on the information $\underline{Z}_{1,t-1}$ (resp. noncausal linear innovation of $Z_{2,t}$ based on the information $\bar{Z}_{2,t+1}$). It is important to note that the information set $\underline{Z}_{1,t-1}$ (resp. $\bar{Z}_{2,t+1}$) differs in general from the global causal information set $\underline{Y}_{t-1} = \underline{Z}_{t-1}$ (resp $\bar{Y}_{t+1} = \bar{Z}_{t+1}$). In Section 5, we clarify this point by introducing the notion of a nonlinear innovation that takes into account all available information. Before doing that, we need to derive the expressions of forward and backward predictive densities. \setcounter{equation}{0}
This Section presents the forecasting and backcasting methods based on closed-form expressions of forward and backward predictive densities.
The expression of the predictive density of $Y_{T+1}$ given $\underline{Y}_T = (Y_T, Y_{T-1},...)$ is given below and derived in Online Appendix A.1., assuming that the matrices of autoregressive coefficients $\Phi_1, \Phi_2,...$ and joint error density $g$ are known.
{\bf Proposition 1:} The conditional probability density function (pdf) of $Y_{T+1}$ given $\underline{Y}_T $ is:
where $l_2(z_2)$ is the stationary pdf of $Z_{2,t} = A^2 \left(
\right)$, if $n_2 \geq 1$. In the pure noncausal process, $n_2=0$, we have: $ l(y| Y_T ) = g(y- \Phi_1 Y_T - \cdots - \Phi_p Y_{T-p+1}).$
In the special case of the mixed VAR(1) process with $p=1$, the predictive density becomes:
Proof: See Appendix A.1.
This predictive density is defined for the VAR representation (2.1). It is a semi-parametric function of the parameters $\Phi_1,...,\Phi_p$, determining $A^2$ and $J_2$, and of the functional parameters $g, l_2$. This predictive density is complicated, except if it is multivariate Gaussian. In general, the conditional mean is not linear in $\underline{Y}_{T}$.
From equation (3.1), the deterministic relation (2.4) between $Y_t$ and $Z_t$, and the symmetry of the calendar and reverse time scales, it follows that:
{\bf Corollary 1: Markov Property} The mixed VAR(p) process $(Y_t)$ of dimension $m$ [resp. the state process $(Z_t)$ of dimension $n=mp$] is a Markov processes of order $p$ [resp. of order 1] in calendar time for $n_2 \geq 1$. The processes $(Y_t)$ and $(Z_t)$ are both Markov of orders $p$ and 1, respectively, in reverse time too.
This corollary extends to any autoregressive order $p$ the result of Cambanis and Fakhre-Zakeri (1994), who show that a linear pure noncausal autoregressive process of order 1 is a causal Markov process of order 1. It also extends the Proposition 3.1 in Freis and Zakoian (2019) derived in the univariate case.
The predictive density summarizes the nonlinear causal dynamics of a mixed VAR process. From the predictive densities, we infer the point forecasts and prediction intervals based on its quantiles. These lead us to the oos point predictions and prediction intervals at various horizons for this nonlinear and non-Gaussian process analogous to the linear pointwise predictions and prediction intervals of the causal VAR models.
Since the mixed VAR process of order $p$ is Markov of order $p$ both in calendar and reverse times, we can derive from Proposition 1 the closed-form expression of backward predictive density in reverse time for backcasting. For ease of exposition, we present it below for $p=1$.
{\bf Corollary 2: Backcasting} Let us consider a mixed VAR(1) model. The backward predictive density of $Y_{T-1}$ given $Y_T$ is: $$l_B (y|Y_T) = \frac{l_1(A^1 y)}{ l_1(A^1 y_T)} |det \; J_2| \; g(Y_T - \Phi y),$$ where $l_1$ is the stationary density of $Z_{1t}$ and $g$ is the density of $\varepsilon$.
Proof: See Online Appendix A.1.
By considering jointly Proposition 1 and Corollary 2, we see that the mixed VAR models have a nonlinear dynamic structure extending to nonlinear dynamic framework the standard Kalman filter available for linear Gaussian processes.
In Section 2.3, we mentioned that there exists a multiplicity of real Jordan representations of matrix $\Psi$. It implies a multiplicity of matrices $A$ built from its extended real eigenspaces (in the presence of complex conjugate eigenvalues). However, $\det J_2 = \prod_{j=1}^{n_2} \lambda_j$, where $|\lambda_j| > 1, \, j=1,...,n_2$, is independent of the real Jordan representation. Similarly, the noncausal component $Z_2$ is defined up to a linear invertible transformation. Since the Jacobian is the same for the numerator and denominator of the ratio $l_2 \left[A^2 \left(
\right)\right] \; /\;l_2 \left[ A^2 \left(
\right) \right]$, it has no effect on the ratio. Thus, the expression of $l(y|Y_T)$ does not depend on the selected real Jordan representation, that is on the selected state-space representation. The same remark applies to the backward predictive density.
The closed-form expressions of the backward and forward predictive densities at horizon 1 can be used sequentially to forecast or back-cast out-of-sample at any horizon $h>1$. Alternatively, the predictive densities at horizons $h>1$ can be obtained by using the Sampling Importance Resampling (SIR) method. Then, the forward predictive density at horizon $h$ can be written as a multivariate integral over $h$ future values of the process and approximated for given $\Phi_1,..,\Phi_p$ and $g$ by drawing the future values of the process from the closed-form predictive density at horizon 1 replicated $S$ times by the SIR method [Gelfand and Smith (1992), Tanner (1993)] (see Section 6 for an illustration). Specifically, to approximate the predictive density at horizon $h>1$, for given $\Phi_1, \Phi_2,...$ and $g$, we need S independent future paths of the process. Each path $s=1,...,S$ is obtained by forecasting sequentially $Y_{T+1}^s| Y_T$, followed by $Y_{T+2}^s| Y_{T+1}^s,...,Y_{T+h}^s|Y_{T+h-1}^s$, from the predictive densities. This becomes a drawing $Y_{T+1}^s,...,Y_{T+h}^s$ of a future path $s$. By replicating it independently for $s=1,...,S$, we get $Y_{T+h}^s, \; s=1,...,S$. Then, for a large $S$, we can use the sample distribution of $Y_{T+h}^s, \; s=1,...,S$ as an estimator of the predictive distribution at horizon $h$.
\setcounter{equation}{0}
The parameters of a mixed VAR model need to be estimated before the forecasts are computed. Below, we review the existing estimation methods in the time domain and describe the predictive algorithm providing the estimated predictive densities in a semi-parametric setup. As pointed out in Section 3, from the quantiles of the estimated predictive density we obtain the point and interval forecasts, which can be impacted by the preliminary estimation step involving estimators converging at different rates. Therefore, in the second part of this section we study the forecast interval uncertainty, using the theory of random sets [see Molchanov and Molinari (2018)].
The mixed VAR model can be estimated by the maximum likelihood method based on an assumed parametric distribution of $\varepsilon_t$ [see Breidt et al. (1991), Lanne and Saikkonen (2013), Davis and Song (2020)]. This approach yields consistent estimators provided that the parametric distributional assumption is correct (and non-Gaussian for identification).
Alternatively, the mixed VAR model can be consistently estimated without any parametric assumptions on the error distribution \footnote{Except for the non-Gaussianity assumption and the uniform tail parameter.} by using the semi-parametric (Generalized) Covariance (GCov) estimator [Gourieroux and Jasiak (2023)]. The GCov estimator is consistent, asymptotically normally distributed and semi-parametrically efficient. As an alternative in the frequency domain, minimum distance estimators based on the cumulant spectral density of order 3 and 4 have been proposed in Velasco and Lobato (2018) and Velasco (2022).
The prediction methods introduced in Section 3 for given $\Phi_1,...,\Phi_p$ and $g$ can be applied in a parametric or semi-parametric framework by replacing these unknown parameters by their consistent estimates. In the semi-parametric framework, this can be done along the following lines:
step 1. Apply the GCov estimator based on zero auto-covariance conditions of nonlinear error functions to obtain the estimators of matrices of autoregressive coefficients $\hat{\Phi}_1, ..., \hat{\Phi}_p$.
step 2. Use the $\hat{\Phi}_i, i=1,...,p,$ estimates to compute the roots of the lag-polynomial and more generally an estimated real Jordan representation: $\hat{A}, \hat{J}_1, \hat{J}_2$.
step 3. Compute the approximated model errors using the estimates obtained in Step 1: $\hat{\varepsilon}_t = Y_t - \hat{\Phi}_1 Y_{t-1} -... - \hat{\Phi}_p Y_{t-p}$.
step 4. Compute $\hat{Z}_t = \hat{A}^{-1} \left(
\right), \;\;\; \hat{\eta}_t = \hat{A}^{-1} \left(
\right)$.
step 5. The following densities can be estimated by kernel estimators applied to the approximated series:
- the density $g$ of $\varepsilon_t$ can be estimated from $\hat{\varepsilon}_t, \; t=1,...,T$;
- the density $l_2$ of $Z_{2,t}$ can be estimated from $\hat{Z}_{2,t}, \; t=1,...,T$ (see Online Appendix B);
step 6. The predictive density can be estimated from the formula (3.1) by replacing $l_2$ by $\hat{l}_2$, $A^2$ by $\hat{A}^2$, and also $J_2$ by $\hat{J}_2$, $g$ by $\hat{g}$, and $\Phi_1, ..., \Phi_p$ by $\hat{\Phi}_1, ..., \hat{\Phi}_p$. The mode (median) of the predictive density provides the point forecasts and the quantiles of the estimated predictive density can be used to obtain estimated prediction intervals at horizon 1 (see, Online Appendix B).
The estimated model parameters and residuals $\hat{\varepsilon}_t$ can be used to build oos predictions and prediction intervals conditional on given values of the last observations in the sample, called the conditional prediction interval. In finite sample, there is uncertainty on the estimated prediction intervals resulting from semi-parametric and non-parametric estimators with different convergence rates. The estimation errors of the scalar and functional parameters have a nonlinear effect on this uncertainty, which should not be disregarded.
To highlight the specificities of the nonlinear dynamic model analysis, we compare our approach with the standard practice of predicting from a linear AR(1) model defined below:
{\it Example 1: Linear AR(1) model:} The model is given by: $y_t = \phi y_{t-1} + \epsilon_t, \; |\phi|<1$ where the errors $\epsilon_t$ are i.i.d. with the unknown true density $f_0$.
Below, we compare the true and estimated prediction intervals, and introduce the bootstrap-adjusted prediction intervals.
{\bf 4.2.1 True and Estimated Prediction Intervals}
i) {\bf True Prediction Interval:} For ease of exposition, let us consider the VAR(1) model, a forecast horizon $h=1$ and a future value of the first component series $Y_{1,T+1}$ to be forecast at date $T$ out of sample (oos) given $Y_T = (y_{1,T}, y_{2,T})' \equiv y$, where $(Y_{2,t})$ contains the remaining components of the series. Then, the true prediction interval at level $1-\alpha_1$ for $Y_{1,T+1}$ is:
where $P_{0}$ is the true density function of process ($Y_{t}$) and $Q(y, \alpha, P_0)$ denotes the $\alpha$-quantile of $Y_{1,T+1}$ conditional on $Y_T=y$, derived from the joint multivariate predictive density $l(y | \underline{Y}_T)$ (see Proposition 1 for the closed-form expression of the predictive density). Let $Q_l$ and $Q_u$ denote the true $\alpha_1/2$ and $1-\alpha_1/2$ conditional quantiles, respectively. Recall that under the semi-parametric approach, the true density $P_{0}$ is characterized by the parameter matrix $\Phi_{1,0}$ and functional parameter $g_0$. Then, conditional on the information at time $T$, the prediction interval is $PI(Y_T, \alpha_1)$. This is a random interval (set), a function of $Y_T$, and its conditional coverage level is $P_0 (Y_{1, T+1} \in PI(Y_T, \alpha_1)|Y_T) = 1- \alpha_1, \forall Y_T. $
{\it Example 1: Linear AR(1) model, cont.:} The theoretical prediction interval is: $PI(y, \alpha_1) = ( \phi y - Q_0 (\alpha_1/2), \phi y + Q_0 (1-\alpha_1/2)),$ where $Q_0$ is the quantile function corresponding to the true distribution function of $\epsilon_t$.
In the mixed VAR model $\epsilon_t$ is no longer a causal innovation and the marginal quantiles $Q_0$ have to be replaced by the conditional quantiles $Q(y, \alpha_1, P_0)$ depending on the current value $Y_T=y$. Moreover, this conditional quantile is not an affine function of $y$.
By using the standard expression of a Gaussian prediction interval with $\Phi$ denoting the Gaussian cumulative distribution function (c.d.f), the asymptotically valid prediction interval (4.1) for $Y_{1,T+1}$ can be equivalently written as:
with\footnote{Note that $\Phi^{-1}(\alpha_1/2)$ is negative.} $ m(y, \alpha_1; P_{0}) = 0.5[Q_l (y,\alpha_1 P_{0})+ Q_u (y,\alpha_1; P_{0})] $ and $\sigma(y,\alpha_1; P_{0}) = - [1/(2\Phi^{-1}(\alpha_1/2)] $
$[Q_u (y, \alpha_1; P_{0})- Q_l (y, \alpha_1; P_{0})]$. This normalized representation of prediction interval resembling the traditional Gaussian approach can be used even if the conditional density function of $(Y_t)$ is not Gaussian. In particular, the functions $ m(y, \alpha_1; P_{0})$ and $ \sigma(y, \alpha_1; P_{0})$ are nonlinear in $y$, in general.
{\bf ii) Estimated Prediction Interval:} The unknown marginal predictive density function $P_0$ can be consistently estimated from eq. (3.1) and denoted by $\hat{P}$, given the estimated matrix of autoregressive parameters $\hat{\Phi}_1$ and the residuals $\hat{\varepsilon}_t$ obtained from the semi-parametric GCov estimator, along with the nonparametric estimator $\hat{g}$. Then, the estimated prediction interval for $Y_{1,T+1}$ is:
where $Q_l(y, \alpha_1; \hat{P})$ and $Q_u(y, \alpha_1; \hat{P})$ are the $\alpha_1/2$ and $1-\alpha_1/2$ conditional quantiles of the estimated predictive density. This estimated prediction interval (4.3) is consistent of the true prediction interval (4.2) when the number of observations tends to infinity. We know that the estimated prediction interval $\widehat{PI}(Y_T,\alpha_1)$ does not always satisfy the conditional coverage condition in finite sample: $P_0 [Y_{1, T+1} \in \widehat{PI}(Y_T, \alpha_1) | Y_T] = 1-\alpha_1, \;\; \forall Y_T$. In general, we have: $P_0(Y_{1, T+1} \in \widehat{PI}(Y_T, \alpha_1)|Y_T) = 1 - \alpha_1 (Y_T),$ i.e. the conditional coverage depends on $Y_T$, and the unconditional coverage: $P_0(Y_{1, T+1} \in \widehat{PI}(Y_T, \alpha_1)) = E_0[1 - \alpha_1 (Y_T)],$ differs from $1-\alpha_1$\footnote{Because the process is Markov of order 1, we condition on $Y_T$ only instead of $\underline{Y}_T = (Y_T, Y_{t-1},...)$, even though the estimators depend on the entire past.}.
{\it Example 1: Linear AR(1) model, cont.:} The estimated prediction interval becomes: $\widehat{PI}(Y_T, \alpha_1) = [\hat{\phi} Y_T - \hat{Q}(\alpha_1/2), \hat{\phi} Y_T + \hat{Q}(1-\alpha_1/2)].$ In the linear dynamic model, the estimated quantiles are easily derived from the residuals $\hat{\epsilon}_t = y_t - \hat{\phi}y_{t-1}$ ranked in increasing order.
Such a simple non-parametric estimation method of the quantile function is not available in the nonlinear dynamic framework, including the mixed VAR(1) model, where $\epsilon_t$ is no longer a causal innovation. Instead, the conditional quantiles are estimated from the kernel functional estimator of the predictive density. Note also that the estimates of $Q$ are expected to be less accurate than the estimate of the autoregressive coefficient $\phi$, because they approximate non-parametrically the tails of a multivariate distribution.
iii) {\bf 4.2.2 Backward Bootstrap Adjusted Prediction Interval}
Since the estimated prediction interval is random, its finite sample distribution can be approximated by a "backward" bootstrap, i.e. by replicating the trajectory of the process by backcasting, conditional on $Y_T$ to adjust for the bias in coverage conditional on $Y_T$. More precisely, given $\hat{\Phi}_1$, $\hat{g}$ considered fixed and the residuals, we can generate by backcasting the artificial paths $Y_t^s, \; t=1,...,T$, with the same terminal condition $Y_T^s=Y_T=y$ for all the bootstrapped samples (see, Corollary 2 for the closed-form expression of the backward predictive distribution). We propose the following algorithm:
step 1: Starting from $Y_T$, we backcast $Y_{T-1}^s$, conditional on $Y_T$, next we backcast $Y_{T-2}^s$ conditional on $Y_{T-1}^s$, and so on.
step 2: By replicating the backcasted path $S$ times, we end up generating $S$ bootstrapped series $Y_t^s, t=1,...,T$, of length $T$ equal to the length of the initial series and with the same terminal value $Y_T$.
step 3: From each replicated path $(Y_t^s, t=1,...,T)$, we estimate the model parameters $\Phi^s$ and $\hat{g}^s$, $s=1,...,S$. This allows us for computing at $Y_{T}$ the predictive density estimator $\hat{P}^s, s=1,...,S$ of $Y_{T+1}$ given $Y_T$ and $S$ new prediction intervals for $Y_{1,T+1}$ from each of the replicated paths.
step 4: The bootstrapped prediction interval obtained from a replicated path is:
where $\hat{P}^s$ is the semi-parametric estimate of $P_{0}$ from the generated path $(Y_t^s, t=1,...,T)$. This bootstrap PI of $Y_{1,T+1}$ given $Y_T$ can be replicated independently $S$ times. The components of the prediction interval (4.4) are denoted by:
step 5: For large $S$, the joint sample distribution of $[\hat{m}^s (y,\alpha_1), \hat{\sigma}^s (y,\alpha_1)]$ provides an approximation of the conditional distribution of $[m(y, \hat{P}), \sigma(y,\hat{P})]$ given $Y_T$, when $T$ is finite and sufficiently large.
{\it Example 1: Linear AR(1) model, cont.:} In the linear AR(1) model, the bootstrap is usually applied by drawing independently in the sample distribution of residuals. This is justified by the additive decomposition of the conditional quantile.
In a nonlinear dynamic model with nonlinear dependence between $y$ and $Q$, the bootstrap has to be performed conditional on $Y_T$ to adjust for the bias in coverage conditional on $Y_T$.
The estimated prediction interval $\widehat{PI}(y,\alpha_1)$ is a pointwise estimator of interval $PI(y,\alpha_1)$ and, as any estimator, is random itself. This randomness is difficult to assess since $\widehat{PI}(y, \alpha_1)$ is a random interval [see Molchanov and Molinari (2018) for random sets in Econometrics]. Let us now extend the method of pointwise estimation to build confidence sets for $PI(y, \alpha_1)$. Since there does not exist a total ordering on intervals, we constrain the confidence set to be also of the form:
and search for an estimator of $q$ providing the correct asymptotic coverage. This confidence set has the following conditional coverage probability of the true prediction interval: $$
$$
\setcounter{equation}{7} because $\Phi^{-1}(\alpha_1/2) <0$. For $\alpha_2 \in (0,1)$, possibly different from and less than $\alpha_1$, there exists a value $q_0(y, \alpha_1; \alpha_2)$ such that:
Asymptotically, we get the conditional $1-\alpha_2$ coverage probability of the true conditional prediction interval, although the true predictive density $P_0$ and the true distribution of $\hat{P}$ remain unknown. In practice, this procedure can provide us a prediction interval at level 95%, for example, when the available sample does not contain sufficiently many tail observations and a reliable prediction interval at level 90% or lower is only estimable.
Then, equations (4.7) and ((ref)) can be replaced by their bootstrapped counterparts obtained from $S$ replicated paths of the series, which are backcast prior to T, given $Y_T=y$. Next, a forecast at T+1 is computed from each of the replicated paths. More precisely, the bootstrapped conditional coverage probability is defined as: $$ \hat{\Pi}^s(y,\alpha_1,q) = \frac{1}{S} \sum_{s=1}^S \delta^s, $$ where $\delta^s = \left\{
\right. $
Then, we consider a solution $\hat{q}^S(y,\alpha_1, \alpha_2)$ of:
ensuring a $1-\alpha_2$ conditional coverage probability. The bootstrap confidence set for the prediction interval is:
and $ \lim_{T \rightarrow \infty} \lim_{S \rightarrow \infty} P_0 [ \widehat{CSPI}(y, \alpha_1, \alpha_2) \supset PI (y,\alpha_1)$ $|Y_T=y] = 1-\alpha_2, \; \forall y, \forall P_0. $ Hence, the length of the estimated $\widehat{PI} (y,\alpha_1)$ is modified by a factor $\hat{q}^S (y,\alpha_1, \alpha_2)/|\Phi^{-1}(\alpha_1/2)| $, that depends on the observed value $Y_T=y$, in general. In practice, we can choose $\alpha_1=\alpha_2 = 0.05$ corresponding to the standard levels for prediction and confidence intervals, respectively. As mentioned earlier, we can choose $\alpha_1$ different from $\alpha_2$, including $\alpha_1 = 1.0$, which would correspond to the confidence set for a point prediction equal to the median of the predictive density.
The analysis of confidence set for the prediction interval is related to set identification, where a confidence set for the identified set is determined for models under partial identification [Beresteanu et al., (2017)]. The partial identification literature considers a parametric model with a partly identifiable parameter\footnote{See Imbens and Manski (2004) for confidence intervals of identified intervals in the framework of partial identification and confidence intervals that asymptotically cover the true interval with a probability larger or equal to $1-\alpha_2$.} to determine either the confidence set under the classical approach, or the credible set under the Bayesian approach. In our framework, we have two types of "parameters": $P$ and $Y_{1, T+1}$. The second one is not identifiable, although its conditional distribution is estimated, and plays the role of a conditional prior.
\setcounter{equation}{0}
The inference on the traditional causal VAR model commonly includes the IRF analysis, which is an important tool used by economists and financial policy makers for the so-called causal analysis. This terminology differs from the "causal-noncausal" terminology introduced by Rosenblatt for time series. More precisely, the "causal analysis" involves a) making inference on the future given the past, i.e. forecasting, and also b) the counterfactual analysis of the effects of current transitory shocks on the future. The standard linear causal IRF analysis cannot be applied to the mixed VAR(p) processes, because neither errors $\varepsilon_t$ in model (2.1), nor the state-specific linear errors $\eta_t$ in (2.5) are causal and independent of the lagged values of $Y_t$. Hence, it is difficult to interpret the shock on $\varepsilon_t$ that implicitly changes a past that has already been realized. This section introduces an alternative concept of innovation that eliminates such a difficulty in nonlinear IRF analysis. This new innovation accounts for nonlinear dynamics of $(Y_t)$ with bubbles and spikes. In addition, it satisfies both the serial and cross-sectional independence conditions [Gourieroux and Jasiak (2005), Gourieroux, Monfort and Renne (2017) on structural SVAR models in Macroeconomics]. We discuss below the filtering and identification of the nonlinear causal innovations in mixed VAR models.
The nonlinear causal innovations are defined below for any Markov process of order $p$, including any mixed VAR(p) model.
{\bf Definition 1:} Let us consider a Markov process of order $p$ with a positive continuous transition density. A nonlinear causal innovation of process $(Y_t)$ is a process $(v_t)$ of dimension $m$ such that: i) the vectors $v_t$ are serially i.i.d.; ii) the strictly stationary process $(Y_t)$ can be written in a nonlinear causal autoregressive form:
with $\underline{Y_{t-1}} = (Y_{t-1},...,Y_{t-p})$; iii) the future values $\bar{v}_t$ of the innovation process are independent of the lagged values $\underline{Y_{t-1}}$ of the observed process $(Y_t)$.
It is easy to see that conditions i) and ii) are equivalent to the Markov of order $p$ property of $(Y_t)$ with a continuous distribution, and can be applied in particular to the mixed model (2.1) by Corollary 1. Such a nonlinear causal autoregressive form always exists [see Rosenblatt (1952) and the discussions in Sections 5.2, 5.3]. It is not unique in a nonparametric framework since $v_t$ can always be replaced by $w_t = d(v_t)$, where $d$ is a diffeomorphism, and $a$ is replaced by $a \circ d^{-1}$. If function $a$ is invertible with respect to $v$, then $v_t$ is a nonlinear function of $Y_t$ given its past and condition iii) is satisfied. Under conditions i) and ii) the autoregressive equation (5.1) can be applied recursively.
The sequence of nonlinear causal innovations provides the basis of nonlinear impulse response functions [Gourieroux and Jasiak (2005), Gonzalves et al. (2021)]. Let us consider a transitory shock $\delta$ at date $T$ on $v_T$, and then apply recursively the autoregressive equation (5.1) to get the IRF, defined as shocked future values $Y_T^{\delta} = a(Y_{T-1}, Y_{T-2},\cdots, v_T + \delta)$, $Y_{T+1}^{\delta} = a( Y_{T}^{\delta}, Y_{T-1},..., v_{T+1})$, and so on. The conditions i) and iii) in Definition 1 are crucial for the interpretation of these shocks in terms of "causal analysis". Condition iii) means that $v_T$ can be shocked at any time T without an effect on the realized past $\underline{Y}_{T-1}$. Condition i) means that this shock has no effect on the future values $\bar{v}_{T+1} = (v_{T+1}, v_{T+2},...)$.
The nonlinear IRF defined above corresponds to a multivariate shock $\delta$. The literature on "causal analysis" is also interested in univariate shocks to specific variables. It is possible to define shocks to a specific component $v_{1,t}$ of $v_t$, if additionally the component $v_{1,t}$ and ($v_{2,t},...,v_{n,t}$) are cross-sectionally independent.
In the linear causal VAR model, there are identification problems concerning the parameters of the model and the innovation-based IRFs, because the innovations to be shocked are defined up to an orthogonal linear transform. Both these identification problems are solved by identifying the mixing matrix $D$ in the $\varepsilon_t=D v_t$ error representation, using the Independent Component Analysis (ICA) method for example [see Gourieroux, Monfort and Renne (2017)].
The nonlinear autoregression (5.1) depends on function $a$ and the distribution of the causal nonlinear innovation $v_t$. In the mixed VAR(p) model, both depend on $\Phi_1, \Phi_2,...$ and $g$, without being one-to-one functions of these parameters, because the nonlinear function $a$ depends on more arguments than the density $g$. This new nonparametric identification issue concerning the causal innovation can be examined in a functional framework by using the properties of harmonic functions. In particular, we get the following result:
{\bf Proposition 2 }[Gourieroux and Lee (2025)]: In the nonlinear causal autoregressive representation of a mixed VAR(1) process, the dimension of under-identification in the functional space of nonlinear autoregressive models is finite and equal to $2 m$.
This result reveals the challenges and complexity of the identification of the nonlinear causal autoregression. The above under-identification is very large and exceeds the under-identification of linear causal VAR models. As pointed out in Section 5.1, it is possible to restrict the analysis to nonlinear autoregressions with Gaussian innovations, such that the components $v_{1,t},...,v_{n,t}$ are independent N(0,1). However, this additional restriction is insufficient to solve the identification issue. Indeed, there exist nonlinear transformations of $v_t$, with respect to which the multivariate N(0, Id) distribution remains invariant [see Gourieroux and Lee (2025)], corresponding to "local matrix rotations".
{\bf Corollary 3: Functional multiplicity} There exists a functional multiplicity of innovation processes ($v_t$), which are cross-sectionally independent standard normal variables.
It follows that the functional identification of the nonlinear autoregressive representation of a causal nonlinear innovation and nonlinear IRFs is a complicated issue in nonlinear autoregressive models, and in particular in the mixed VAR models. It resembles the identification issue of causal VAR and Structural (SVAR) models, which has remained unsolved for a long time. Until a better solution becomes available, we can use one of the two following methods:
i) We can introduce identifying restrictions on function $a$ through parametric assumptions on the joint distribution of errors $\varepsilon_t$. For example, one could assume that this distribution is a non-Gaussian elliptical or t-Student. However, such a partial identifying assumption is difficult to justify due to the lack of a structural interpretation of errors ($\varepsilon_t$).
ii) A recursive structure of the model can be assumed, by analogy to the recursive structure of errors in causal VAR models that has been traditionally used for (S)VAR identification in the early macroeconomic literature [see Sims (1980) based on Cholesky decomposition].
In practice one may prefer to adopt method ii). Note that since the observed process ($Y_t$) and the state process ($Z_t$) satisfy a one-to-one relationship, the state variables also satisfy a nonlinear autoregressive scheme:
where $\tilde{v}_t$ can be serially conditionally i.i.d. (and possibly restricted to be standard Gaussian). A recursive structure can be imposed on model (5.2) with the noncausal state variable to be shocked first under this approach (also referred to as "shock ordering"). This choice can be further motivated by the fact that in applied research, the noncausal order is often equal to 1, i.e. $n_2 =1$ [see e.g. Hecq, Lieb and Telg (2016), Gourieroux and Jasiak (2017) and Section 6.3]. This implies that there is a common noncausal component generating the bubbles and local trends in $(Y_t)$. For example, in macroeconomic models, that noncausal component can be related to speculative bubbles in oil prices impacting jointly the price index, GDP and other macroeconomic variables. In financial applications, the single factor capturing the nonlinear dynamics can be interpreted as a systemic component, or a common bubble [Gourieroux and Zakoian (2017), Cubadda et al. (2023), Hall and Jasiak (2024)]. For a bivariate VAR(1) with $n_2=1$, we get an identifiable recursive form of the nonlinear autoregressive model (5.2):
where function $G_2$ in $Z_{2,t} = G_2 ( Z_{2,t-1}; v_{2,t})$ is the inverse of the function defining the nonlinear Gaussian causal innovation of $Z_2$, which is uniquely defined as:
where $F_2$ is the conditional cumulative distribution function (c.d.f) of $Z_{2,t}$ given $Z_{t-1}$ and $\Phi$ is the c.d.f. of the standard Gaussian distribution. From equation (a.5) given in Appendix A.1, it follows that the conditional density of $Z_{2,t}$ given $Z_{t-1}$ has a closed form given by:
where $g_{\eta_2}$ is the marginal density of $\eta_2$. The c.d.f. $F_2$ is the integral of $l(z_{2,t}| z_{t-1})$ over the set of admissible values of $z_{2,t}$. By the Markov property of $(Z_{2,t})$ in Corollary 4 below, $Z_{1,t-1}$ does not appear in equation (5.6) and we have $v_{2,t} = \Phi^{-1} [ F_2 (Z_{2,t}|Z_{2,t-1})]$.
Next, we append $v_{2,t}$ by the (Gaussian) ($v_{1,t}$) associated with the causal state variable(s) in this recursive model. In our illustration, we consider the case $n_1=1$ and define:
where $F_{1|2}$ is the conditional cumulative distribution function of $Z_{1,t}$ given $Z_{2,t}, Z_{t-1}$ \footnote{The system of equations (5.6)-(5.8) is the inverse of system (5.5)-(5.7) and it can be used for simulating the path of $Z_t$ from independent standard Gaussian drawing of $v_{1,t}, v_{2,t}$.}. From equation (a.5) given in Appendix A.1, we get the closed-form expression of the conditional density:
The c.d.f. $F_{1|2}$ is the integral of $l(z_{1,t}| z_{2,t}, z_{t-1})$ over the set of admissible values of $z_{1,t}$.
When the noncausal state variable $Z_2$ is the first one to be shocked, it is easy to check that $v_{1,t}(Z)$ differs in general from the causal innovation $\eta_{1,t}$ associated with the latent causal component $Z_{1,t}$. By applying the shock ordering, we are implementing the autoregressive model (5.2) in a recursive form (5.3)-(5.4). Moreover, we deduce from (5.5):
{\bf Corollary 4:} The noncausal state variable $(Z_{2,t})$ is a Markov process of order 1 with respect to the filtration associated with $(Y_t)$, [and also to the filtrations associated with $(Z_t)= (Z_{1,t}, Z_{2,t})$, and with $(Z_{2,t})]$.
By construction, the structural noncausal shock $v_{2,t} (Z)$ is defined in a unique way (since $n_2=1$) and is a Gaussian white noise\footnote{See Gourieroux and Jasiak (2005) for the uniqueness of a nonlinear Gaussian innovation in univariate models.}. It is also independent of the multivariate shock $v_{1,t} (Z)$. Then, we can trace out nonlinear IRFs based on the effect of a change $\delta_2$ in $v_{2,t}$, with $v_{1,t}$ held constant. It is called henceforth the Common Bubble Shock (CBS) [See Online Appendix C for quantile estimation of $F_2 ( Z_{2,t} | Z_{2,t-1})$]. In the general case $n_2 = 1, n_1 \geq 1$, the associated multivariate impulse responses of $Z_{2,t}$ are identifiable, i.e. independent of other components ($v_{1,t}, v_{3,t},...,v_{n,t}$).
In practice, when there is a single noncausal state variable, the filtering algorithm of Section 4 can be completed by the two following steps:
step 7. The approximations of nonlinear causal innovations $\hat{v}_{2,t}(Z)$ can be computed from the estimated distributions as $\hat{v}_{2,t}(Z) = \Phi^{-1}(\hat{F}_{2,T} (\hat{Z}_{2,t}| \hat{Z}_{2,t-1}))$ by applying the formula of predictive density (5.6) with $l_2$ replaced by $\hat{l}_2$ and $g_{\eta_2}$ replaced by $\hat{g}_{\eta_2}$, i.e. the kernel-smoothed empirical density of $\hat{\eta}_{2,t}, \, t=1,...,T$.
step 8. Next, $\hat{v}_{2,t}(Z)$ can be appended by the approximated nonlinear causal innovations $\hat{v}_{1,t}$, independent of $\hat{v}_{2,t}(Z)$, which are defined as: $\hat{v}_{1,t} = \Phi^{-1}(\hat{F}_{1|2,T} (\hat{Z}_{1,t}| \hat{Z}_{2,t}, \hat{Z}_{t-1})), \; \mbox{for} \; n_1=1$.
These causal innovations can be computed from the predictive density formula (5.8) with $g_{\eta_2}$ and $g_{\eta}$ replaced by their empirical counterparts.
\setcounter{equation}{0}
This Section illustrates the nonlinear forecasts from the mixed VAR(1) model and the nonlinear IRF analysis. Subsection 6.1 presents a simulation study that examines the oos forecasts from the model. In subsection 6.2, the semi-parametric GCov estimators, nonlinear innovations and IRFs are illustrated in an application to a bivariate series of US GDP rates and oil prices. Additional simulations and empirical results are provided in Online Appendices D and E.
We consider a simulated bivariate mixed VAR(1) process with the following matrix of autoregressive coefficient: $\Phi = \left(
\right), $ with eigenvalues 0.7 and 2, located inside and outside the unit circle, respectively . The errors follow a bivariate noise with independent components both t-student distributed with $\nu$ =4 degrees of freedom, mean zero, and variance equal to $\nu/(\nu-2)=2$. The matrix A is as follows: $A = \left(
\right). $ The simulated paths of the series of length 500 is displayed in Figure 1. The solid (black) line represents process $(Y_{1t})$ and the dashed (red) line represents process $(Y_{2t})$. \noindent The sequence of spikes in the noncausal component $Y_{2t}=Z_{2t}$ impacts the component $Y_{1t}$ through the recursive form of matrix $\Phi$.
Let us now consider the forecasts based on parameter estimates. The Generalized Covariance (GCov) estimate of matrix $\Phi$ is obtained by minimizing the portmanteau statistic computed from the auto- and cross-correlations up to and including lag $H=10$ of the errors $\varepsilon_t = Y_t - \Phi_1 Y_{t-1}$ and their squared values [see Gourieroux and Jasiak (2023)]\footnote{The finite sample properties of the GCov estimator are illustrated in Gourieroux and Jasiak (2023).}. The estimated autoregressive matrix is $\hat{\Phi} = \left(
\right),$ with eigenvalues $\hat{\lambda}_1=0.690$, $\hat{\lambda}_2=2.027$, which are close to the true values $\lambda_1 = J_1 = 0.7$ and $\lambda_2 = J_2 = 2.0$. The standard errors of $\hat{\Phi}$ obtained by bootstrap are 0.023, 0.308, for the elements of the first row, and 0.009, 0.120 for the elements of the second row.
After estimating $\Phi$, the GCov estimated errors $\hat{\varepsilon}_t=Y_t - \hat{\Phi} Y_{t-1}$ are computed. Matrices $A$ and $A^{-1}$ are identified from the real Jordan representation of matrix $\Phi$, up to scale factors. The estimated matrix $\hat{A}^{-1}$ computed from the normalized Jordan decomposition of $\hat{\Phi}$ is:
$ \hat{A}^{-1} = \left(
\right)$. It corresponds to the true matrix $A^{-1} = \left(
\right)$ up to scale factors of about 0.02 and 0.97 for each column. We use these estimates to approximate the causal and noncausal components displayed in Figure 2.
Let us now consider the oos nonlinear forecast one step ahead performed at date $T=500$ when the process takes values $Y_{1, 500} =-3.367$ and $Y_{2, 500}=-0.239$. The true values of $Y_{1,501}$ and $Y_{2, 501}$ are -2.260 and -0.331, respectively. The nonlinear forecasts are summarized by the predictive density, whose mode provides pointwise predictions of $Y_{1, T+1}$ and $Y_{2, T+1}$. It is estimated from formula (3.2) one-step ahead oos by using a kernel estimator over a grid of 100 values below and above $Y_1$ and $Y_2$, with Gaussian kernels and bandwidths $h_2=1$ and $h_{11} = s.d.(\varepsilon_1), h_{12} = s.d.(\varepsilon_2)$ (see Online Appendix B for kernel density estimators). The estimated point forecasts are $\hat{Y}_{1,501}= -2.80$ and $\hat{Y}_{2,501}= -0.30$. The estimated prediction intervals at level 0.80 determined from the predictive density are [-4.80, -0.80] for $Y_{1,501}$ and [-2.60, 2.10] for $Y_{2,501}$. The rationale for choosing level 80% is to ensure a sufficiently large number of observations in the tails for reliable estimation of the quantiles of predictive density. Both prediction intervals contain the true future values of the process. Additional results on the coverage of the estimated prediction interval and on the estimated prediction set uncertainty are provided in Online Appendix D.
In this Section, we apply the mixed VAR(1) model to analyse jointly the real oil prices and the real GDP growth rates.
{\bf a) The data}
We examine the quarterly series of oil prices and US GDP growth rates over the period: Q1 1986 - Q2 2019. The real US GDP growth rate series is calculated from Real Gross Domestic Product, Quarterly, Seasonally Adjusted Annual Rate available at at https://fred.stlouisfed.org from the Federal Reserve Economic Data.
The oil prices are provided online by the US Energy Information Administration under the Short-Term Energy Outlook Real and Nominal Prices, March 2023 \footnote{ called Quarterly Average Imported Crude Oil Price/barrel, Real Price, deflated by the US Consumer Price CPI index.}. This series approximates "the price of oil paid by US refiners for crude oil purchased from abroad" examined by Kilian and Vigfusson (2017). The series of GDP rates and oil prices (divided by 10) of length T=134 are displayed in Figure 3. We observe, for example, in year 2008, a pronounced bubble in oil prices with a strong negative impact on the GDP growth rate.
{\bf b) Estimation}
We estimate the causal-noncausal VAR(1) from the demeaned series of GDP growth rates (series 1) and demeaned oil prices divided by 10 (series 2). The sample mean of growth rates is 0.640 and the mean of rescaled oil prices is 6.348. The estimated autoregressive matrix is $\hat{\Phi} = \left[
\right]$ with standard deviations of autoregressive coefficients of 0.014, 0.001, 0.032, 0.004, respectively. The eigenvalues are 0.294 and 1.073. The densities of estimated errors provided in Online Appendix E, Figure a.5 are non-Gaussian.
The presence of a noncausal root reflects the nonlinear dynamic features such as the spikes and a bubble observed in Figure 3. The spikes and bubble in oil prices are accommodated by the strictly stationary mixed VAR(1) model of GDP rates and oil price levels.
{\bf c) Prediction}
We perform out-of-sample (oos) predictions from the mixed VAR(1) at the points indicated in Figure 3 below, which mark the bubble episode at the end of the sample to ensure a sufficiently large number of prior observations for estimation.
During the selected period, both series displayed several small sudden changes, before the oil price dropped. The predictive densities are evaluated at each point over a grid of 200 points, equidistant by 0.1 below and above the last observed value of each variable. We use again Gaussian kernels and bandwidths $h_2=s.d.(Z_{2})$, $h_{11} = s.d.(\varepsilon_1)$, and $h_{12} = s.d.(\varepsilon_2)$ to estimate one-step ahead oos predictive density (3.2) [See Online Appendix E, Figures a.7-a.11 for predictive density plots which vary in $T$\footnote{because of the increasing information set and the dependence on the current environment of $Y_T$. } and become bi-modal on-bubble when the probability of crash increases].
Table 1 reports the point forecasts obtained from the main modes of predictive densities, and the forecast intervals obtained from the sample quantiles of estimated predictive densities. The forecast interval is at level 80% to ensure a sufficient number of observations to estimate the quantiles of the predictive density and to alleviate the potential effect of bimodality. We also consider the predictions of the last 2 points in the trajectory. The results in Table 2 are computed according to the procedure described earlier, starting from the conditioning time T=132 when growth rate $Y_{1, 132}=-0.458$ and oil $Y_{2, 132}=0.268$ are closer to their averages.
These prediction intervals are conditional and depend on the date and distance in time from a bubble. We observe that the prediction intervals on-bubble are longer than off-bubble. In addition, we use the approach outlined in Section 4 to estimate the confidence set of the prediction interval for $Y_{1,134}$. It is based on 100 backcast replications of the series conditional on $Y_T$. The confidence set of prediction interval for the demeaned growth rate at level $1-\alpha_2 =0.95$ is: $\widehat{CSPI}(y, \alpha_1, \alpha_2) = [-1.546, 1.405]$. This interval is at a level higher than the prediction intervals in Table 2, accounts also for the estimation risk and is longer.
{\bf d) Bubble shock}
Let us now examine the effect of a set of the following values of the CBS: $\delta_2= -2, -1, 0, 1, 2$ to $z_{2,T}$ performed at date $T=110$ during the bubble episode and at time $T=133$ after the bubble, under the identifying assumption of a recursive form (5.3)-(5.4). Since $Z_t= A^{-1} Y_t$, that shock can originate from either the oil prices, or growth rates, or both.
The impulse responses in Figure 4 are computed conditionally on the path of the process up to and including date $T$. Because the dynamic model is nonlinear, the IRF is nonlinear in $\delta_2$, and the shock effects need to be compared with a baseline. The baseline path of $Z_{2, T+1}^b,...,Z_{2, T+10}^b$ is computed from 10 random values of standard Normal used as the future $v_{2, T+1},...,v_{2, T+10}$. Next, a shock $\delta$ is added to $v_{2, T+1}$ and the consecutive values of shocked $Z_{2, T+1}^{\delta},...,Z_{2, T+10}^{\delta}$ are calculated recursively by inverting the conditional c.d.f. [See Online Appendix C].
To interpret the shocks, we compare the relative size of shock effects to the baseline, and we conclude that a shock has dissipated when the shocked path overlaps with the baseline. We observe that, on-bubble, shocks to $Z_2$ are more long-lasting and more symmetric around the baseline. Off-bubble, we observe that a shock of $\delta=-2$ dissipates much slower than other shocks. Next, we perform shocks of the same size to the causal component $Z_1$ using a similar approach (Figure 5). We find that shocks to $Z_1$ have weaker effects and dissipate quickly, which is consistent with the fact that the bubble is driven mainly by the noncausal component $Z_2$.
The shocked component $Z_2$ can be combined with the values of causal component $Z_1$ conditional on its own past and the shocked values of $Z_2$. This is done by predicting $Z_1$ from the mode of the conditional density $l(z_{1,t}|z_{2,t}, z_{1, t-1})$ given in (5.8) and kernel estimated. Then, we can approximate the shocks to $y_1$ and $y_2$ as $Y^{\delta} = A Z^{\delta}$. For illustration, Figure 6 displays the shocked values of $Y_{2, T+1}^b,...,Y_{2, T+10}^b$.
We find that $Y_2$ is more responsive to on-bubble shocks. Off-bubble, the shocks dissipate faster except for the large positive shock, which again takes longer to disappear.
This paper considers oos nonlinear forecasting and backcasting in a mixed VAR model. It introduces a closed-form expression of forward (resp. backward) predictive density for forecasting (resp. backcasting) from mixed VAR models. As a post-estimation inference method, we adjust the estimated prediction interval by a conditional "backward" bootstrap and introduce the confidence set of the estimated prediction interval. A definition of causal (past-dependent) nonlinear innovations for mixed VAR models is also given. Since the causal nonlinear innovations are not uniquely defined, their identification is examined and the IRF analysis following a shock to the noncausal component is discussed.
For illustration, the proposed approach is applied to the analysis of the joint dynamics of a bivariate series of oil prices and real GDP growth rates. We find that the noncausal state variable captures the explosive patterns, including bubbles and spikes. We examine the IRFs and observe different effects of a shock to the noncausal component on- and off-bubble. The mixed VAR models offers a parsimonious representation of nonlinear dynamic processes and allow for advanced analysis of the dynamics during a bubble episode.
\singlespacing {\bf Funding:} This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).
Beresteanu, A., Molchanov, I., and F. Molinari (2017): "Partial Identification Using Random Set Theory", Journal of Econometrics, 16, 17-32.
Blasques, F., Koopman, S., Mingoli, G. and S., Telg (2025): "A Novel Test for the Presence of Local Explosive Dynamics", Journal of Time Series Analysis, forthcoming.
Breidt, F., Davis, R., Lii, K., and M., Rosenblatt (1991) : "Maximum Likelihood Estimation for Noncausal Autoregressive Processes", Journal of Multivariate Analysis, 36, 175-198.
Cambanis, S., and I., Fakhre-Zakeri (1994): "On Prediction of Heavy-Tailed Autoregressive Sequences Forward versus Reversed Time", Theory of Probability and its Applications, 39, 217-233.
Cavaliere, G., Nielsen, H., and A., Rahbek (2020): "Bootstrapping Noncausal Autoregressions with Applications to Explosive Bubbles Modelling", Journal of Business and Economic Statistics, 38, 55-67.
Cubadda, G., Hecq, A., and S., Telg (2019): "Detecting Co-Movements in Non-Causal Time Series, Oxford Bulletin on Economics and Statistics, 697-715.
Cubadda, G., Hecq, A., and E., Voisin (2023): "Detecting Common Bubbles in Multivariate Mixed Causal-Noncausal Models", Econometrics, 11, 9.
Cubadda, G., Giancaterini, F., Hecq, A. and J. Jasiak (2024): "Optimization of the Generalized Covariance Estimator in Noncausal Processes", Statistics and Computing, 34, 127.
Davis, R., and L., Song (2020) : "Noncausal Vector AR Processes with Application to Economic Time Series", Journal of Econometrics, 216, 246-267.
De Truchis, G., Freis, S. and A., Thomas (2025): "Forecasting Extreme Trajectories Using Semi-Norm Representations", working paper Paris Dauphine University.
Fries, S., and J.M., Zakoian (2019): "Mixed Causal-Noncausal Autoregressive Processes", Econometric Theory, 35, 1234-1270.
Gelfand, A., and A., Smith (1992): "Bayesian Statistics without Tears: A Sampling-Resampling Perspective", Annals of Statistics, 46, 84-88.
Giancaterini, F., Hecq, A., Jasiak, J. and A. Manafi-Neyazi (2025): "Regularized Generalized Covariance Estimator", ArXiv 6390395.
Gonzalves, S., Herrera, A., Kilian, L., and E., Pesavento (2021): 'Impulse Response Analysis for Structural Dynamic Models with Nonlinear Regressors", Journal of Econometrics, 225, 107-130.
Gourieroux, C., and A., Hencic (2015) : "Noncausal Autoregressive Model in Application to Bitcoin/USD Exchange Rates", Econometrics of Risk, Studies in Computational Intelligence, 583:17-40
Gourieroux, C., and J., Jasiak (2005) : "Nonlinear Innovations and Impulse Responses with Application to VaR Sensitivity", Annals of Economics and Statistics, 78, 1-31.
Gourieroux, C., and J., Jasiak (2016): "Filtering, Prediction, and Simulation Methods for Noncausal Processes", Journal of Time Series Analysis, 37, 405-430.
Gourieroux, C., and J., Jasiak (2017): "Noncausal Vector Autoregressive Process: Representation, Identification and Semi-Parametric Estimation", Journal of Econometrics, 200, 118-134.
Gourieroux, C., and J., Jasiak (2023): "Generalized Covariance Estimator", Journal of Business and Economic Statistics, 41, 1315 -1357.
Gourieroux, C., Jasiak, J., and M., Tong (2021): "Convolution-Based Filtering and Forecasting: An Application to WTI Crude Oil Prices", Journal of Forecasting, 40, 1230-1244.
Gourieroux, C., and Q., Lee (2025): "Identification and Impulse Response Functions for Nonlinear Dynamic Models", ArXiv 2506.13531.
Gourieroux, C., Monfort, A., and J.P., Renne (2017) : "Statistical Inference for Independent Component Analysis", Application to Structural VAR Models", Journal of Econometrics, 196, 111-126.
Gourieroux, C., and J.M., Zakoian (2017) : "Local Explosion Modelling by Noncausal Process", Journal of the Royal Statistical Society, B, 79, 737-756.
Hall, M., and J., Jasiak (2024): "Modelling Common Bubbles in Cryptocurrency Prices", Economic Modelling, 139, 106782.
Hecq, A. , Lieb, L. and S., Telg (2016): "Identification of Mixed Causal-Noncausal Models in Finite Samples", Annals of Economics and Statistics, 123/124, 307-331.
Herrera, A., Lagalo, L., and T., Wada (2015) : "Asymmetries in the Response of Economic Activity to Oil Price Increases and Decreases ?", Journal of International Money and Finance, 50, 108-133.
Imbens, G., and C., Manski (2004): "Confidence Intervals for Partially Identified Parameters", Econometrica, 72, 1845-1857.
Kilian, L., and R., Vigfusson (2017) : "The Role of US Oil Price Shocks in Causing US Recessions", Journal of Money, Credit and Banking, 40, 1747-1776.
Lanne, M., and J., Luoto (2016): "Noncausal Bayesian Vector Autoregression", Journal of Applied Econometrics, 31, 1392-1406.
Lanne, M., and P., Saikkonen (2011) : "Noncausal Autoregressions for Economic Time Series", Journal of Time Series Econometrics, 3, 1-39.
Lanne, M., and P., Saikkonen (2013) : "Noncausal Vector Autoregression", Econometric Theory, 29, 447-481.
Lof, M., and H., Nyberg (2017): "Noncausality and the Commodity Currency Hypothesis", Energy Economics, 65, 424-433.
Molchanov, I., and F., Molinari (2018): "Random Sets in Econometrics", Econometric Society Monographs, Cambridge University Press.
Nyberg, H., and P., Saikkonen (2014): 'Forecasting with a Noncausal VAR Model", Computational Statistics and Data Analysis, 76, 536-555.
Perko, L. (2001): "Differential Equations and Dynamical Systems", Springer, New York.
Rosenblatt, M. (1952) : "Remarks on Multivariate Transformations", Annals of Mathematical Statistics, 23, 470-472.
Rosenblatt, M. (2012) : "Gaussian and Non-Gaussian Linear Time Series and Random Fields", Springer Verlag.
Sims, C. (1980): "Macroeconomics and Reality", Econometrica, 48, 1-48.
Swensen, A. (2022): "On Causal and Non-Causal Cointegrated Vector Autoregressive Time Series", Journal of Time Series Analysis, 42, 178-196.
Tanner, M. (1993): "Tools for Statistical Inference", Springer Series in Statistics, 2nd edition, Springer, New York.
Twumasi, C., and J., Twumasi (2022): "Machine Learning Algorithms for Forecasting and Backcasting Blood Demand Data with Missing Values and Outliers: A Study of Tema General Hospital of Ghana", International Journal of Forecasting, 38, 1258-1277.
Velasco, C. (2023): "Identification and Estimation of Structural VARMA Models Using Higher Order Dynamics", Journal of Business and Economic Statistics, 41, 819-832.
Velasco, C., and I., Lobato (2018): "Frequency Domain Minimum Distance Inference for Possibly Noninvertible and Noncausal ARMA Models", Annals of Statistics, 46, 555-579.
\doublespacing
\setcounter{equation}{0}