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.
90,332 characters · 11 sections · 65 citation commands
Spectral estimation for mixed causal-noncausal autoregressive models
\\
The traditional assumption of causality in autoregressive models is widely adopted in classical estimation methods for economic time series. The causality implies that the time series depends on a combination of lags with roots outside the unit circle, forcing the time series to depend only on past information. Despite this, causality is not always justified and is not systematically tested. Additionally, it emerges that including only past information for estimating economic models fails to explain their structural shocks, as shown by hansen2019two. In addition, lanne2012optimal,lanne2012has show that including future information may improve the forecast accuracy, i.e., through noncausal model representations. A statistically rigorous presentation of noncausal models arises in breid1991maximum, dividing the autoregressive polynomial into two components, one causal with roots outside the unit circle and another noncausal with roots inside the unit circle, or equivalently, including leads in the autoregressive process. Furthermore, allowing the autoregressive polynomial to contain roots simultaneously inside and outside the unit circle, introduces the mixed causal-noncausal models. The noncausal dynamics offer the possibility to capture nonlinearities in weakly stationary variables, such as locally explosive patterns and bubbles, replicating characteristics of anticipation and speculation. Consequently, noncausal models can produce valuable risk measures, fueling their interest in economics, with applications in inflation, commodity prices, and stock indices, see among others gourieroux2017local,lof2013noncausality,lof2017noncausality,lanne2011noncausal,hecq2016identification, hecq2021forecasting. \\
An interesting aspect of the noncausal models is the estimation and identification. The reason is that they share the same second-order dependence structure, the autocovariance function or its frequency domain analog, the spectral density, with the causal models. Consequently, the estimations based only on second-order moments such as Gaussian-likelihood, OLS, Yule-Walker, or Whittle yield coefficients corresponding to roots outside the unit circle identifying a causal model, even though the actual roots of the coefficients are inside the unit circle, corresponding to a noncausal model. The underlying reason for this lack of identification is the assumption of Gaussianity. For the Gaussian series, the probability structure is entirely determined by second-order moments. In general, for an AR($p$) model with a Gaussian errors sequence, there exist $2^{p}$ possible specifications of its roots with the same probability structure. In contrast, when the error sequence is non-Gaussian, each of the $2^{p}$ configurations of roots corresponds to a different weakly stationary process with a unique and identifiable probability structure. Evidently, knowing the non-Gaussian probability distribution of the error sequence, the natural estimator is the Maximum Likelihood (ML), identifying the true roots of the autoregressive polynomial regardless of whether it is causal or noncausal, see basawa1976asymptotic. This estimation was initially proposed by lii1992approximate for possible noninvertible MA models, proving that it is consistent for any symmetric, strictly positive, and sufficiently smooth non-Gaussian probability density function. huang2000quasi propose a quasi-Likelihood estimation of possible noninvertible MA models using the Laplacian distribution. breid1991maximum proposes an Maximum likelihood estimation (MLE) for possibly noncausal autoregressive processes driven by non-Gaussian i.i.d. noise, evidencing that the estimated parameters are asymptotically normal and providing identification when the error sequence is heavy-tailed. lii1996maximum, lehr1998maximum, implement MLE for possible noncausal and noninertible ARMA models. Similarly, lanne2011noncausal, hecq2016identification, giancaterini2022inference propose a parameter estimation in possible noncausal AR models using the Student's t distribution. For the same distribution, andrews2006maximum proposes the estimation of the all-pass ARMA models, where all roots of the autoregressive polynomial are reciprocals of roots of the moving average polynomial and vice versa. MLE for autoregressive models using the alpha-stable distribution has also been implemented in andrews2009maximum, fries2021conditional.\\
A natural criticism of MLE is the prior knowledge about the probability distribution of the error sequence, which is uncommon empirically. Indeed, arbitrarily assuming a distribution function that does not replicate the features of the data can lead to misestimating parameters and misidentifying the autoregressive model. Alternatively, in this research, we propose to estimate the parameters of the causal, noncausal, and mixed autoregressive model without prior knowledge about the probability distribution of the error sequence. Instead, it has been proposed to employ the higher-order cumulants necessary for constructing the higher-order spectral densities, with the ability to discriminate between causal and noncausal models. Higher-order cumulants, also named higher-order correlations, have been used to estimate noncausal AR models, see rosenblatt1980linear,lii1982deconvolution,giannakis1990estimating, tugnait1986identification, terdik1999bilinear,nikias1987non. We use in this paper the procedure proposed by brillinger1985fourier, leonenko1998spectral,velasco2018frequency, constructing a function of the second and third-order spectral densities (the spectrum and the bispectrum, respectively) in a minimum distance estimation with the periodogram and biperiodogram. The estimation of the parameters is performed by minimizing the function, achieving the identification of the model in its global minimum.\\
Although this estimation is consistent and asymptotically normal, the global minimum is not guaranteed to be obtained. On the contrary, the estimation function is multimodal (not convex), leading to multiple local minima that increase as the number of parameters in the model grows. Accordingly, the solution is not unique and convergent. For a discussion of the multimodality in noncausal models, see kindop2021ubiquitous, and for mixed causal-noncausal, see bec2020mixed. Consequently, an empirical prerequisite for reaching the global minimum in the estimation function is the appropriate selection of the initial values. In this sense, our first contribution to the literature is to develop a general scheme for selecting the initial values for estimating noncausal and mixed causal-noncausal based on their causal representations. This allows us to find the correct combination of parameters in the region where the global minimum of the estimation function lies. In our second contribution, we present a decision method for identifying the model based on the existence of the global minimum in the estimation function. We explore the robustness of our method in finite samples through an extensive Monte Carlo study using the alpha-stable distribution. We find unbiased estimated parameters and a correct identification as the data depart from normality. Our third contribution is to perform the third-order spectral estimation of mixed models. To the best of our knowledge, this is the first time in mixed models. In addition, we propose an empirical application of eight monthly commodity prices, evidencing noncausal and mixed causal/noncausal dynamics. \\
The paper is divided as follows: Section 2 presents the methodology, Section 3 presents the estimation method, Section 4 introduces the Monte Carlo study, Section 5 presents the empirical application, and the last Section wraps up the conclusions and discussion.
Let $y_t$ be a time series of length $T$, with zero mean and finite $k$-th order moments, for $k\geq2$, and $k$-th weakly stationary for $t=\left\{0,\pm 1, \pm2, \pm 3, ... \right\}$, with autocovariance function
Notice that the autocovariance is considering linear relations of the data in two time periods, known as second-order dependence, also that $\kappa_{2_j}=\kappa_{2_{j-1}}...=\kappa_{2_{j-n}}=\kappa_{2_{j+1}}...=\kappa_{2_{j+n}}$. In any of these representations, the autocovariance function remains unchanged. Therefore, estimations based on the second-order structure, e.g., Gaussian-Likelihood, OLS, or Yule-Walker, fail to distinguish between causal and noncausal models. The above is evident in Figure (ref). On the left is a causal model and on the right is a noncausal model with the same coefficients and the same alpha-stable error sequence. Notice that the autocorrelation functions are identical.
However, when the time series is non-Gaussian, the probability structure of causal and noncausal models is different and identifiable using additional information from the higher-order cumulants. We specify that the first three moments are equivalent to the first three cumulants. This is not valid for cumulants of orders greater than three. In particular, the third cumulant contains the relation of the time series at three points in time. Hence it is known as the triple correlation function or the bicovariance function, see bartelt1984phase. The third-order cumulant is
Note that $\kappa_3(j,l)$ is equal to zero for the Gaussian and symmetrical distributions. The definition of $\kappa_3(j,l)$ is the third coefficient of the power expansion of the natural logarithm of the characteristic function of a random variable. Given the stationarity of $y_t$, $\kappa_3(j,l)$ satisfies the following symmetry conditions
The second and third-order cumulants are necessary for constructing the spectrum and the bispectrum. The spectrum is the Fourier transform of the second cumulant. Analogously, the bispectrum is the Fourier transform of the third-order cumulant. Defining the Fourier frequencies as $\omega=2\pi/T$, the spectrum is
where $i$ is the imaginary number. The bispectrum is
Note that for the existence of the bispectrum, the third-order cumulant must be different from zero. We assume that the spectrum and the bispectrum have continuous second-derivatives and are bounded above and away from zero in modulus for the sets of frequencies $\Lambda_1$ and $\Lambda_2$. Where $\omega_1 \subset \Lambda_1=(0,1/2)$ for $S_2$, and given the symmetries $S_3(\omega_1,\omega_2)=S_3(\omega_2,\omega_1)=S_3(\omega_1,-\omega_2-\omega_1)=S_3(-\omega_1-\omega_2,\omega_2)$ and the identity $S_3(\omega_1,\omega_2)=S_3(-\omega_1,-\omega_2)$, $\Lambda_2$ is a finite union of compact sets inside an open triangle with vertices $(0,0),(1/2,0),(1/3,1/3)$, for details see rao1989estimation,rosenblatt1965estimation,rao2012introduction, anh2007minimum.
The nonparametric estimation of the spectrum and bispectrum for autoregressive processes is usually based on methods employing the discrete Fourier transform (DFT), i.e., assuming that the time series comprises a set of harmonically related sinusoids. Let us define the DFT as
The most popular estimator of the spectrum is the periodogram. For a summary of spectrum estimates, see kay1981spectrum. The periodogram is obtained as the modulus of the output values from the DFT in Equation (ref), performed directly on the time series. The origin of the periodogram can be traced back to schuster1898investigation, which conducted a Fourier series fitting to the variation of sunspot numbers trying to find hidden periodicities in the data. Thereafter, it has been widely used in spectral analysis and autoregressive modeling, see brillinger1975time, rosenblatt1965stationary,alekseev1996asymptotic,bartlett1950periodogram,brillinger1967asymptotic. The periodogram is defined as
Similar to the periodogram, however, taking into account relations between three frequencies (bifrequencies), the biperiodogram is the bispectrum estimator, defined as
where $\overline{d_T(-\omega_1-\omega_2)}$, represents the complex conjugate of the Fourier transform of the sum of two frequencies. This component is the one that accounts for the asymmetries present in non-Gaussian data. Contrary to the periodogram, which is always real, the biperiodogram is complex due to its last factor. The periodogram and biperiodogram are asymptotically unbiased estimators of the spectrum and bispectrum but inconsistent; i.e., the variance of the estimates will not tend to zero as $T$ increases without bound, see brillinger1975time, rosenblatt1965stationary,alekseev1996asymptotic,bartlett1950periodogram. Some properties of the periodogram are developed in brillinger1967asymptotic; in particular, the mean of the periodogram $k=2$ and biperiodogram $k=3$. Denoting $\boldsymbol{\omega}={\omega_1,\omega_2}$
when $\omega_s\neq 0\;mod\;2\pi$. The variance for $k={2,3}$
as $T\rightarrow\infty$. Equivalently to the periodogram, the biperiodogram at different Fourier frequencies $(\omega\neq\omega')$ is uncorrelated
In this Section, we present causal, noncausal, and mixed autoregressive models. Let $y_t$ be a causal autoregressive process generated by
The lag polynomial is $\phi(L)=1-\phi_1 L-\phi_2 L^2...\phi_r L^r$, with a set of parameters $\Phi_r=(\phi_1,\phi_2...\phi_r)$, where $L$ stands for the backward shift operator, that is, $L^k=y_{t-k}$ for $k=\left\{0,\pm 1, \pm2, \pm 3, ... \right\}$. This model is well known in the time series literature as the AR($p$). In our notation, the AR($p$) model is replaced by the AR($r,0$) model. The main difference between the AR($p$) model and the AR($r,0$) model is that in the former, it is sufficient to assume white noise in the error sequence, while in the latter, we assume that the error sequence is i.i.d and has finite moments greater than 3, that is
The same assumption applies to the AR($0,s$) and MAR($r,s$) models, discussed in Sections (ref) and (ref), respectively. Note that we differentiate the cumulants of $y_t$, $\kappa_2,\kappa_3$, with the cumulants of the error sequence $\varepsilon_t$, $\kappa_{2_e},\kappa_{3_e}$. $\phi(z)$ is provided by the roots $1-\phi_1 z-\phi_2 z^2...-\phi_r z^r$, all lying outside the unit circle. Hence, the model has a covariance-stationary representation form $\kappa_2(j)=\mathbb{E}\left[y_t y_{t-j}\right]$, and is represented as an infinite sum of past shocks, MA($\infty$) as
where $\sum_{j=0}^\infty\lvert\Psi \rvert<\infty$. Then, the autocovariance generating function of $y_t$ is $g_y(z)=\sum_{j=-\infty}^{\infty} \kappa_2(j) z^j$. Evaluating $g_y$ at $z=e^{- i\omega}$, and dividing by $2\pi$ we obtain the spectrum of $y_t$
This function is known as the second-order spectrum since it only considers second-order relations. The spectrum is continuous, non-negative, real-valued, symmetric, and periodic around $\omega$. For the AR($r,0$) model, the spectrum can be obtained explicitly using the representation of Equation (ref). Thus, we can define the transfer function of the linear filter AR($r,0$) as
Then, the spectrum of the AR($r,0$) is the multiplication of the second-order cumulant of the errors sequence and the transfer function, and its conjugate
Defining $\overline{z}=e^{2\pi i (-\omega_1-\omega_2)}$, the transfer function of the sum of the complex conjugation of the frequencies is $\overline{\psi(\Phi_r,-\omega_1-\omega_2)}=1/(1-\phi_1\overline{z}-...-\phi_2\overline{z}^2-...-\overline{z}^r)$. In this way, the bispectrum is
Let $y_t$ be a noncausal autoregressive process generated by
where the lead polynomial is $\varphi(L^{-1})=1-\varphi_1L^{-1}-\varphi_2L^{-2}...\varphi_s L^{-s}$ with the set of parameters $\Phi_s=(\varphi_1,\varphi_2...\varphi_s)$. In this notation, the AR($0,s$) model has $s$ noncausal and zero causal components. The polynomial $\varphi(z)$ is $1-\varphi_1 z-\varphi_2z^{2}...-\varphi_s z^{s}$, with all roots outside the unit circle. Note that the polynomial of the noncausal model can be represented with its roots outside the unit circle, thus depending on its leads, or with roots inside the unit circle depending on its lags. Indeed, the polynomial $\varphi(z)$ has a causal representation with its zeros inside the unit circle, (see brockwell2009time p. 81 for a stationary solution), that is $\varphi^{*}(z)=1-\varphi_1^{*}z-\varphi_2^{*}z^2...\varphi_s^{*}z^s$, with $\varphi^{*}(z)\neq0$ for $\lvert z\rvert>1$. Following Equation 6 in lanne2011noncausal, the polynomial $\varphi^{*}(z)$ can be expressed as
where $\varphi(z^{-1})$ is as in Equation (ref), implying that $\varphi^{*}_{s-j}/\varphi^{*}_s=-\varphi_j$ for $j=1,2,...s-1$, and $1/\varphi^{*}_s=\varphi_j$. Notice that the zeros of $\varphi^{*}(z)$ lie inside the unit circle. In this way, the noncausal model can be rewritten to obtain the following causal specification
where the error term $\xi_t=-1/\varphi_{s}\varepsilon_{t-s}$, and coefficients for the polynomial $\varphi^{*}_j$ equal to $\varphi^{*}_j=-\varphi_j/\varphi_s$ for $j=1,2...s-1$, and $\varphi^{*}_s=1/\varphi_s$. Note that this causal representation with roots inside the unit circle should not be confused with the explosive model, where the error term is contemporaneously uncorrelated with lags of $y_t$ (see brockwell2009time, p. 81). As an example, in the Remark 1, we rewrite the noncausal AR(0,$s$) into a purely causal AR($r,0$), for $s=r$.\\
Remark 1. The causal AR($r,0$) representation of the noncausal AR($0,s$) model, for $s=r$. The noncausal AR($0,s$) is
Reverting the process in $s$ time periods
Solving for $y_t$, we can rewrite $y_t$ as
which highlight the correlation between the error term and the regressors. Then, the purely causal AR($r,0$) representation of the AR($0,s$) is
Notice that the coefficients and the error sequence from Equation (ref) are equivalent to the ones in Equation (ref). The causal representation is key to the estimation and identification of the model, as we shall discuss in Section (ref).\\
The transfer function for the noncausal AR(0,$s$) is $\psi(\Phi_s,\omega)=1/(1-\varphi_1z-\varphi_2z^{2}...-\varphi_rz^{s})$, and its spectrum is
Additionally, the bispectrum is
In Figure (ref), we show the same models as in Figure (ref), but their frequency domain representation. At the top of the graph is the periodogram and spectrum of the causal AR(1,0) $y_t=0.7y_{t-1}+\varepsilon_t$. At the bottom is the periodogram and spectrum of a noncausal AR(0,1) model $y_t=0.7y_{t+1}+\varepsilon_t$. $\varepsilon_t\sim$ alpha-stable $f_a(t;\alpha=1.5,\beta=0.25,\gamma=1,\delta=0)$, the sample size $T=100$. As in Figure (ref), the same error sequence is used to construct the models. The periodogram, $I_2(\omega)$, as well as the spectrums $S_2(\Phi_r\omega)=S_2(\Phi_s\omega)$ are the same in both cases. Consequently, estimations based on second-order moments identify a causal model, even though the DGP is a noncausal model. In both cases, the low frequencies (below 0.2, and above 0.8) contain most of the information in the spectrum. In contrast, the high frequencies do not contain information about the process.
Let $y_t$ be a mixed causal-noncausal autoregressive process MAR($r,s$) generated by
where $\phi(L)=1-\phi_1L-\phi_2L^2...\phi_r L^r$ corresponds to the causal component, and $\varphi(L^{-1})=1-\varphi_1L^{-1}-\varphi_2L^{-2}...\varphi_s L^{-s}$ to the noncausal. The traditional AR($p$) model is decomposed to the MAR($r,s$) model by the equivalence of its order $p=r+s$. The set of parameters is $\Upsilon_{r,s}=\left\{\phi_1,...,\phi_r,\varphi_1,...,\varphi_s\right\}$. The MAR($r,s$) is a combination of the variable $y_t$ with $r$ lags and $s$ leads, and both polynomials $\varphi(z)$ and $\phi(z)$, with roots outside the unit circle. Alternatively, the polynomial $\varphi(z)$ has a causal representation with all its roots inside the unit circle, and the polynomial $\phi(z)$ outside the unit circle. For details, see lii1996maximum. Accordingly, $y_t$ has a two-sided moving average representation as in Equation (ref). The polynomial of the MAR($r,s$) model is $\upsilon_{r,s}(z)=(1-\varphi_1z-\varphi_2z^{2}...-\varphi_sz^{s})(1-\phi_1z-\phi_2z^2...-\phi_r z^r)$, this interaction between the components of the transfer function results in nonlinearities that can hinder the estimation of the parameters.
The MAR($r,s$) has a causal purely causal representation AR($r+s$,0), as follows
This model representation can be found in gourieroux2018misspecification, Appendix A. Where $\xi^{*}_t=-1/\varphi_s\varepsilon_{t-s}$, and $\upsilon^{*}_{r+s}$ is the polynomial of the causal AR($r+s$,0) representation. To obtain it, we follow several steps. First, we apply distributive law for polynomial $\upsilon_{r,s}(L)$
Reverting the process in $s$ time periods, or equivalently multiplying by $L^s$,
Solving for $y_t$, The closed form of the AR($r+s$,0) is
The polynomial $\upsilon^{*}_{r+s}(z)$ , has the set of parameters $\Upsilon^{*}_{r+s}=\left\{\upsilon^{*}_1,\upsilon^{*}_2,...\upsilon^{*}_s,...\upsilon^{*}_{r+s}\right\}$, then the AR($r+s$,0) can be written as
Note that there is a one-to-one relation between Equations (ref) and (ref). Such representation becomes important in identifying the parameters by correctly selecting the initial values in the estimation. To shed light on this relation, we explore the transformation of the MAR($1,1$) into an AR($2,0$) in the Remark 2.\\
Remark 2. The MAR($1,1$) is represented by
Multiplying by $L$ and solving for $y_t$, we obtain the AR($2,0$) representation
The causal representation is
where the set of parameters is $\Upsilon^{*}=\left\{\upsilon^{*}_1=\frac{1+\phi_1\varphi_1}{\varphi_1},\upsilon^{*}_2=-\frac{\phi_1}{\varphi_1}\right\}$, and $\xi^{*}_t=-1/\varphi_1\varepsilon_{t-1}$. Estimating the set of parameters $\upsilon^{*}$ is not trivial, given the nonlinearity in the original model. For a similar representation of Equation (ref), see gourieroux2018misspecification proposition 3.1.\\
The spectrum of the MAR($r,s$) is
The bispectrum is
where $\psi(\Upsilon_{r,s},\omega)=1/\left[(1-\phi_1z-\phi_2z^2...-\phi_r z^r)(1-\varphi_1z^{-1}-\varphi_2z^{-2}...-\varphi_sz^{-s})\right]$, and $\overline{\psi(\Upsilon_{r,s},-\omega_1-\omega_2)}=1/\left[(1-\phi_1\overline{z}-\phi_2\overline{z}^2...-\phi_r \overline{z}^r)(1-\varphi_1\overline{z}^{-1}-\varphi_2\overline{z}^{-2}...-\varphi_s\overline{z}^{-s})\right]$.
In Figure (ref), we show the real and imaginary parts of the bispectrum and the biperiodogram of simulated data from a MAR(1,1) with coefficients $\phi_1=0.7$ and $\varphi_1=0.2$ and $T=200$, the distribution of the errors is alpha-stable $f_a(t;\alpha=1.5,\beta=0.25,\gamma=1,\delta=0)$ with leptokurtic and skewness features. The set of frequencies $\omega_1$ is in the y-axis and $\omega_2$ in the x-axis. In this Figure, the white spaces indicate that the bispectrum value is relatively small and, therefore, the relation between the set of frequencies $\omega_1$ and $\omega_2$ is low. Conversely, when the Figure turns black, there is a higher relation between the set of frequencies. We can observe that the highest amount of information is stored in the low frequencies for $\omega_1$ and $\omega_2$, from 0 to 0.2 and 0.8 to 1 (see the symmetries in Equation (ref)). This accumulation of information fades as the frequency increases, from 0.2 to 0.5, and from 0.5 to 0.8. In addition, there is an accumulation of information on the $\omega_1=\omega_2$ diagonal that fades as it approaches the frequency 0.5. This same pattern can be seen in the periodogram and spectrum in Figure (ref). Actually, the main diagonal of the biperiodogram (bispectrum) is equivalent to the periodogram (spectrum). Accordingly, the asymmetric patterns that arise from relations between frequencies outside the main diagonal are neglected in the second-order spectrum estimations, thus contributing to the correct estimation of the parameters and the identification of the model in the third-order estimation. Simultaneously, we can observe that the patterns present in the biperiodogram are replicated by the bispectrum. In the real part of the biperiodogram, we can observe two triangular patterns formed by relations between the high frequencies. One below the main diagonal at the frequencies 0.3 to 0.5 and another above the main diagonal at the frequencies 0.5 and 0.7. In addition, four triangular patterns formed a relationship between the high (low) frequencies of $\omega_1$ and low (high) frequencies of $\omega_2$. Two are below the main diagonal, and two are above the main diagonal. Additionally, in the imaginary part of the biperiodogram, there are two triangular patterns at low frequencies, one above and one below the main diagonal. In addition, there are four triangular patterns of the interaction between high frequencies and low frequencies, two above and two below the main diagonal.
In this Section, we present the estimation function that was first proposed by leonenko1998spectral. However, we use the variation presented by velasco2018frequency in terms of the scaling factor in the denominators of Equation (ref). We define $\boldsymbol{\vartheta}$ as the set of parameters of the autoregressive model, for the causal $\boldsymbol{\vartheta}=\Phi_r$, for the noncausal $\boldsymbol{\vartheta}=\Phi_s$ and for the mixed $\boldsymbol{\vartheta}=\Upsilon_{r,s}$. We first perform a preliminary Gaussian likelihood estimation and obtain $\overline{\boldsymbol{\vartheta}}$, always assuming causality, and we evaluate it in the transfer functions $\psi(\overline{\boldsymbol{\vartheta}},\omega)$ and $\psi(\overline{\boldsymbol{\vartheta}},-\omega_i-\omega_j)$ and the second order cumulant, $\overline{k}_{2}=2\pi T^{-1}\sum_{j=1}^{T-1}I_2(\omega_j)/(\psi(\overline{\boldsymbol{\vartheta}},\omega_j)\psi(\overline{\boldsymbol{\vartheta}},-\omega_j))$. This normalization by preliminary estimates has no effect on the identification of the parameters, given that they are invariant to any inversion of the polynomial roots. This step significantly simplifies the estimation compared with leonenko1998spectral. The estimation function is
where $A_{2T}=m(2\pi)^2/(4\overline{\kappa_2}^2T)$ and $A_{3T}=n(2\pi)^4/(6\overline{\kappa_2}^3T^2)$. In all cases, $m=1-n$. The weights in our approach are arbitrarily set as $m=0.5$ and $n=0.5$. terdik1999bilinear estimate an AR($1,0$) with different values for $m$, highlighting the efficiency gains with the correct selection. Nevertheless, in velasco2018frequency, the weights are optimally selected, increasing the efficiency, although they turn out to be unintuitive in some cases, increasing the complexity of the estimation. Considering theorem 2 in velasco2018frequency, we include in the spectrum and the bispectrum a consistent estimator of the standardized cumulants of orders two and three, respectively, that is
For a correct identification, the parameters and cumulants must be jointly estimated. Therefore, we modify the spectrum function to add the dependence of the cumulants to $\boldsymbol{\vartheta}$. The modified spectrum and bispectrum are
When the data is Gaussian, the first term of the estimation function $R_T(\boldsymbol{\vartheta})$ (Equation (ref)) is sufficient to estimate and identify the autoregressive model correctly. The first term in $R_T(\boldsymbol{\vartheta})$ is equivalent to whittle1953estimation or Gaussian-likelihood. Nevertheless, when the data are non-Gaussian, brillinger1975time suggests including the second term, that is, the bispectrum estimation. This estimator is consistent and asymptotic normal. For details, see terdik1999bilinear,velasco2018frequency. The set of optimal parameters $\hat{\boldsymbol{\vartheta}}$ is obtained by minimizing $R_T$ by the quasi-Newton method
For the consistency of the estimator, we use Lemma 1 in brillinger1974cross, and Section 4.2, Theorem 73 in terdik1999bilinear. Let $\boldsymbol{\vartheta}$ be the true set of parameters, $\boldsymbol{\vartheta} \subset \Theta$, and
\vskip 0.5cm Theorem 1. Suppose that $y_t$ has a bispectrum $\lvert S_3(\omega_j,\omega_i)\rvert<\infty$. The first term in $R_T(\boldsymbol{\vartheta})$, $\left(\frac{ I_2(\omega_j)-S^{*}_2(\boldsymbol{\vartheta},\omega_j)}{ \psi(\overline{\boldsymbol{\vartheta}},\omega_j)\psi(\overline{\boldsymbol{\vartheta}},-\omega_j)}\right)^2$, belong to the finite union of closed intervals $\Lambda_1$, and the second term $\frac{\lvert I_3(\omega_j,\omega_i)-S^{*}_3(\boldsymbol{\vartheta},\omega_j,\omega_{i}) \rvert^2}{ \psi(\overline{\boldsymbol{\vartheta}},\omega_j)\psi(\overline{\boldsymbol{\vartheta}},\omega_i)\psi(\overline{\boldsymbol{\vartheta}},-\omega_j-\omega_i)}$ belong to the finite union of compact domains lying inside the open triangle of frequencies $\Lambda_2$. In addition, $R(\boldsymbol{\vartheta})$ has a unique minimum at $\boldsymbol{\vartheta}$. In this circumstances, $R_T(\boldsymbol{\vartheta})\overset{p}{\to}R(\boldsymbol{\vartheta})$ and $\hat{\boldsymbol{\vartheta}}\overset{p}{\to}\boldsymbol{\vartheta}$, as $T\rightarrow\infty$. $\hat{\boldsymbol{\vartheta}}$ is obtained minimizing $R_T(\boldsymbol{\vartheta})$.
\vskip 0.5cm
Under the conditions mentioned above, and assuming that $\hat{\boldsymbol{\vartheta}}$ is close enough to $\boldsymbol{\vartheta}$ when $T\rightarrow\infty$, the Taylor expansion of $R_T(\hat{\boldsymbol{\vartheta}})$ is
As $\hat{\boldsymbol{\vartheta}}$ minimizes $R_T(\boldsymbol{\vartheta})$, $ \frac{\partial R_T(\hat{\boldsymbol{\vartheta}})}{\partial \boldsymbol{\vartheta}}=0$. Then,
Then
The closed forms for the first and second derivatives of Equation (ref) are terdik1999bilinear Section 4. It is reasonable to presume that the numerator in Equation (ref) is zero since $\hat{\boldsymbol{\vartheta}}$ is close enough to $\boldsymbol{\vartheta}$, as the third-order cumulant is different from zero and the bispectrum exists.
\vskip 0.5cm Theorem 2. Under the fulfillment of Theorem 1, $\hat{\boldsymbol{\vartheta}}$ is asymptotically normal as $T\rightarrow\infty$
where $\Sigma_{0}(\boldsymbol{\vartheta})$ and $\Sigma_{1}(\boldsymbol{\vartheta})$ are Equations 4.17 and 4.18 in terdik1999bilinear. For the autoregressive process, the asymptotic variance of $\hat{\boldsymbol{\vartheta}}$ is
where the skewness ($\zeta$) and kurtosis ($\kappa$) of the error sequence are defined by the standardized third and fourth-order cumulants as
The parameter $\eta$ is the Integral respect to $\omega$ of the derivative of the modulus of the second-order transfer function with respect the set of parameters $\boldsymbol{\vartheta}$
By solving Equation (ref) we can obtain the standard errors. Notice that the asymptotic variance depends on the third and fourth-order cumulants. Since moments and cumulants do not exist for all distributions, we use their sample counterparts when they are absent. lobato2022single developed closed formulas for the asymptotic variance of the ARMA model allowing for complex roots. \vskip 0.5cm
Theorem 1 shows that the correct set of parameters can be identified and consistently estimated by finding the global minimum of $R_T(\boldsymbol{\vartheta})$. However, this estimation function is not convex. On the contrary, it is multimodal, leading to the possibility of multiple local minima that increase as the number of parameters in the model grows. As a result, the empirical prerequisite for reaching the global minimum is the correct selection of the initial values, ensuring that $\hat{\boldsymbol{\vartheta}}$ is close enough to the true value $\boldsymbol{\vartheta}$. The latter entails certain limitations for noncausal and mixed models, given their nonlinearity and that $\boldsymbol{\vartheta}$ is not directly observable. Accordingly, the solution of this type of model is not unique and convergent. For a discussion of the multimodality in noncausal models, see kindop2021ubiquitous, and for mixed causal-noncausal, see bec2020mixed. To illustrate this situation, a model that is originally causal AR($2,0$) has three possible solutions for its parameters. A noncausal AR($0,2$), a mixed MAR($1,1$) and a causal AR($2,0$). Evidently, the solutions AR($0,2$) and MAR($1,1$) correspond to local minima of $R_T(\boldsymbol{\vartheta})$ since the global minimum must be for AR($2,0$).
The problem of identification in the noncausal model has been studied by breid1991maximum,andrews2006maximum,lanne2011noncausal in the time domain. In their procedure, they estimate assuming causality and select the order of the AR($p$). In a second step, they estimate the noncausal AR($0,s$) with $p=s$, and all possible combinations of mixed MAR($r,s$), with $r+s=p$, identifying as the appropriate model the one with the highest log-likelihood. Similarly, we propose an estimation strategy that significantly increases the probability of reaching the global minimum of the estimation function by selecting a set of initial parameters close enough to the true set of parameters, helping $\hat{\boldsymbol{\vartheta}}$ to converge in probability to $\boldsymbol{\vartheta}$. Then, we address the multimodality of the estimation function by selecting the model with the lowest $R_T(\hat{\boldsymbol{\vartheta}})$. The initial parameters are chosen based on the pure causal representation of the noncausal (Equation (ref)) and mixed models (Equation (ref)) from which their original roots can be reconstructed. In this context, we developed a tool to identify purely causal, noncausal, and mixed autoregressive models in the frequency domain.
For clarity, $R_{T}(\Phi_r)$ refers to the estimation function when the model is purely causal AR($r$,0), $R_{T}(\Phi_s)$ when the model is purely noncausal AR(0,$s$), and $R_T(\Upsilon_{r,s})$ when is mixed causal-noncausal MAR($r,s$). As a common first step, we perform Gaussian likelihood or Whittle estimation and obtain $\overline{\boldsymbol{\vartheta}}$ and $\overline{k}_2$.
After the estimation, we proceed with the model identification based on the value of the estimation function.
We perform an extensive Monte Carlo study to investigate the finite sample properties of our estimation strategy using $R_T(\boldsymbol{\vartheta})$. We consider the ability to identify the three types of models, causal, noncausal and mixed. The results are based on M=1000 simulations. We use the alpha-stable distribution as the data-generating process (DGP). For implementations of the alpha-stable distribution in noncausal models, see fries2019mixed and fries2021conditional. The probability density function of the alpha-stable is
where $\alpha\in(0,2)$ is the stability parameter, $\beta \left[-1,1\right]$ the skewness parameter, $\gamma>0$ the scale parameter, and $\delta\in\Re $ the location parameter. We generate random numbers with length $T=\left\{100,200,500\right\}$, for the parameters we choose $\alpha=\left\{1.5,1.8\right\}$, $\beta=0.25$, $\gamma=1$ and $\delta=0$. To allow for skewness and leptokurticity. In Figure (ref), we show an alpha-stable MAR($1,1$), the empirical density, the ACF, and the PACF. We can observe asymmetric cycles, bubble patterns, and high persistence in the ACF. Since we are interested in the ability of our method to identify the correct model, we consider three DGPs, a purely causal AR($2,0$), with $\phi_1=0.7, \phi_2=0.2$, a purely noncausal AR($0,2$), with $\varphi_1=0.7, \varphi_2=0.2$, and a mixed MAR($1,1$) with $\phi_1=0.7, \varphi_1=0.2$.
In Table (ref), we present the percentage with which we identify the correct model (identification rate) for the estimation of $R_T(\boldsymbol{\vartheta})$. In the first column, there are the parameters of the alpha-stable distribution, notice that the only parameter that varies is $\alpha$. The second column is the sample size $T$. To compare our method of selecting the initial values and their impact on the identification rate, we compare it with the second-order method, i.e., selecting the initial values directly from the AR($2,0$) Gaussian likelihood. A similar procedure is used in the MARX package hecq2017simulation. The identification rate is the number of times that the method succeeds in detecting that the model is causal AR($2,0$) when the DGP is AR($2,0$), given the three possibilities, causal AR($2,0$), noncausal AR($0,2$), and mixed MAR($1,1$). The same logic follows when the DGP is noncausal AR($0,2$) and mixed MAR($1,1$). We report the identification rate for the causal AR($2,0$) in the third column, in the fourth column for the AR($0,2$), and in the fifth column for the MAR($1,1$).
The identification rate reveals that even with 100 observations, our strategy performs well for all three models. In all cases, the percentage increases as the number of data increases. The model with the highest level of identification is the causal AR($2,0$), followed by the mixed MAR($1,1$) and the noncausal AR($0,2$). As was expected, the more leptokurtic distribution, $\alpha=1.5$, has a higher identification rate since it deviates more from the Gaussian distribution. The results are equivalent in other simulations with different parameters (not presented). However, the identification rate decays slightly as the parameters approach the unit circle. In the case of the parameters taken from the coefficients of the AR($2,0$) model estimation by Gaussian-likelihood, the identification decreases significantly. In this case, the identification is lower in the noncausal AR($0,2$) model than in the mixed MAR($1,1$) model. It is noteworthy that identification is possible when the sample moments of the data exist, even knowing that for the alpha-stable distribution, only the moments of orders two, three, and four exist when $\alpha=2$.
The mean and standard deviation of the estimated parameters are reported in Table (ref), (we only show the results for our estimation strategy).
In all three cases, the mean of the estimated parameters is near the true parameters, showing that the estimation is unbiased and revealing the merits of our estimation strategy. The results are better for the more leptokurtic distribution, $\alpha=1.5$, in all cases, meaning that the estimator decreases in bias as the data deviates from normality. As is expected from Theorem 1, the results improve as $T$ increases. In parallel, the standard deviation decreases as $T$ increases. For the causal and noncausal models, the second parameter deviates slightly from its original value when the sample is 100 and 200. The problem is corrected when the sample is 500. In all cases, its standard deviation is lower than the other distribution $\alpha=1.8$, revealing that the variance of the estimator depends on the sample moments.
For applying the models discussed above and our estimation strategy, we considered eight monthly commodity prices, including Cotton, measured in US cents per pound. Wheat, Barley, Zinc, Copper, Aluminium, and Nickel, are measured in US dollars per metric ton. We also consider the Brent oil, measured in US dollars per barrel. The eight series are from January 1900 through July 2022 for 383 observations. The data was obtained from the FRED (\url{https://fred.stlouisfed.org/}).
An accurate understanding of the dynamics of commodity prices is critical to setting the macroeconomic strategies of exporting and importing countries. In many of them, the stability of their national accounts depends on their behavior. The factors that affect commodity prices are the climate, geopolitical events, the internal policies of each country, transport costs, demand, supply, and profitability, among others. Although commodities are a competitive market with many buyers and sellers, there is evidence that their dynamics can be explained with noncausal autoregressive models. hecq2021forecasting found evidence of noncausality in monthly Nickel prices, hecq2019predicting,gourieroux2021convolution in crude oil monthly prices, karapanagiotidis2014dynamic in 25 commodity futures price, including soft, precious metals, energy, and livestock sectors and lof2017noncausality in the exchange rates of commodity exporters.
All series do not share the same trend, but they appear to be affected by similar shocks. We can observe in Figure (ref) peaks during the 2008 financial crisis, then a drop around 2010 and a subsequent increase from late 2010 to early 2014, then generalized increases in 2022 due to the recent Russian invasion of Ukraine. For the eight series, the null hypothesis of unit root is not rejected in the Augmented Dickey-Fuller test with constant and trend alternatives.
As the data seems non-stationary, a transformation is necessary for the estimation. For instance, hencic2015noncausal use a cubic trend, fries2019mixed use a linear trend or a simple centering, hecq2020mixed, cubadda2019detecting use the returns, hecq2017seasonal use the X-11 method to filtering the seasonal effects. However, arbitrarily selecting the type of trend or the transformation of the data can affect the dynamics of the series, greatly influencing the noncausal component. Alternatively, hecq2019predicting proposes to use the Hodrick-Prescott (HP hereafter) filter with a penalty parameter $\lambda=$129.600 for monthly data, indicating through a Monte Carlo study that the method keeps the data dynamics unaltered. Notwithstanding the above, there are criticisms about the dynamics induced by the HP filter and the selection of its penalty parameter; see hamilton2018you for details. For the sake of discussion, we consider two transformations, the HP filter with a penalty parameter $\lambda= 129.600$, using the cyclical component, and the logarithmic returns. The Figures of these two transformations are respectively (ref), (ref).
Following breid1991maximum, we estimate all series by Gaussian likelihood and selected the order $p$ by second-order information criteria, AIC and BIC. After selecting the order, we fit the AR($p$) model and obtain the residuals. The descriptive statistics, the selected $p$, and the p-value of the Ljung-Box test for the residuals are in Table (ref).
The series are not non-Gaussian for both transformations. This result was corroborated by the Jarque-Bera test (not presented). In the HP transformation, the highest standard deviations are for Nickel and Copper. All the series have a positive skewness except for Copper and Aluminium. Cotton and Nickel are the series with the highest kurtosis. The selected order is $p=2$ in all cases. We construct the residuals for the AR($2$) model and calculate the Ljung-Box test; the p-value for the first two lags is reported. For all the cases, the null hypothesis is not rejected.
For the log returns, the standard deviation is relatively similar for all cases, with Brent and Nickel having the highest values. All the series have negative skewness except for Wheat. All series are leptokurtic, with the highest kurtosis for Brent and Copper. In four cases, the selected order $p=1$, and in four cases, $p=2$. The AR($p$) model residuals are calculated, and we report the p-value of the first two lags of the Ljung-Box test. For all the cases, the null hypothesis is not rejected. For the cyclical component of the HP filter, the competing models are the AR($2,0$), the MAR($1,1$), and the AR($0,2$). In Table (ref), we report the estimation using $R_T(\boldsymbol{\vartheta})$ of the three models. The standard deviation of the coefficients is in parentheses. We also present the residual sum of squares (SSE) and the value of the estimation function for each model.
The results from Table (ref) reveal that the model selected in seven cases is MAR($1,1$), except for Cotton, which is a purely noncausal AR($0,2$). A feature of the causal AR($2,0$) estimation is that the coefficient of the first lag is outside the unit circle, and the second lag is negative and inside the unit circle. This result is equivalent in the noncausal model AR($0,2$), with the coefficient of the first lead outside the unit circle and the coefficient of the second lead, negative and inside the unit circle. Notice that the coefficient outside the unit circle does not violate the stationarity assumption since, in all cases, the modulus of the characteristic roots is greater than one. Gaussian-likelihood also obtains this estimation pattern. Alternatively, for the MAR($1,1$), the noncausal coefficient is greater than the causal coefficient, both positive and significant. The causal coefficient is between 0.319 for Barley and 0.564 for Nickel. The noncausal coefficient is more homogeneous and always positive, between 0.810 for Aluminium and 0.880 for Barley. In most cases, the AR($2,0$) model has the lower SSE. In the last row, we add the estimation using the MARX package hecq2017simulation, that is, Student's t MLE. The results are equivalent in 5 cases. However, different for Wheat AR($2,0$), Aluminium AR($0,2$), and Brent oil AR($0,2$). Similar identification and coefficients were obtained following the time domain estimation of lanne2011noncausal. In Table (ref), we present the estimation using $R_T(\boldsymbol{\vartheta})$ for the log-returns.
We identify an AR(1,0) for Cotton and Wheat, AR(0,1) for Copper and Aluminium, one AR($0,2$) for Barley, and three MAR($1,1$) for Zinc, Nickel, and Brent oil. In the case of the AR(0,1) model, the coefficients are 0.415 for Cotton and 0.207 for Wheat. Both are highly significant. Similar results are obtained for the AR(1,0) model, where the coefficients are 0.383 for Copper and 0.279 for Aluminium. For MAR($1,1$), the causal coefficients are 0.333 and 0.293, and -0.146 for Zinc, Nickel, and Brent oil, respectively. For the same commodities, the noncausal coefficients are -0.090, 0.021, and 0.440, respectively. In the last row, we add the identification of the MARX package using the 's t MLE. The results are equivalent in six cases. they are different for Aluminium AR(0,1) and Brent oil AR($0,2$).
It is noted that the possible dynamics found in the data are dependent on the transformation performed. In the case of the HP filter, in seven of the eight possible cases, the model identified in its cyclical component is MAR($1,1$). On the other hand, the log-return transformation results are more heterogeneous, with three MAR($1,1$), two AR($0,1$), one AR($0,2$), and one AR($1,0$). The results are similar using Student's t MLE using the MARX package. In the case of the HP filter, the identification is equivalent in five cases. In the case of the log returns, it is equivalent in six cases.
The real and imaginary parts of the bispectrum and the biperiodogram of the MAR($1,1$) for the Brent oil cyclical component are presented in Figure (ref). We observed that the highest value relations are found between the low frequencies, 0 to 0.25 and 0.8 to 1. In this range, we observe asymmetric relations, reflecting the non-Gaussianity of the data. No clear relations are observed between the middle frequencies, from 0.25 to 0.4, and the low frequencies. Neither between the high frequencies, from 0.4 to 0.5, with the low frequencies. The same pattern can be observed in the main diagonal of the real and imaginary parts of the biperiodogram.
In Figure (ref), we present the graphical representation of the estimation function $R_T(\hat{\Upsilon}_{p,s})$ for the cyclical component of the Brent oil. On the x-axis is the noncausal coefficient ($\varphi_1$), and on the y-axis is the causal coefficient ($\phi_1$). The surface features have been truncated to accentuate the areas where the minima are located. We observe three regions where minima occur. Two of them correspond to local minima, one for the causal model AR(2,0) and one for the noncausal AR(0,2). The global minimum of the estimation function is for the MAR(1,1) with parameters $\varphi_1=0.811$ and $\phi_1=0.459$, located at the intersection of the vertical (green) and horizontal (red) lines.
This paper presents a frequency domain estimation method based on a combination of the spectrum and the bispectrum for causal, noncausal, and mixed autoregressive models. We test its properties in finite samples by conducting a Monte Carlo study, finding unbiased estimated parameters and asymptotic normality in the variance. We also performed an empirical application on eight monthly commodity prices finding evidence of noncausality and mixed causality/noncausality.
Our contribution to the current literature is divided into three parts. First, we develop a general scheme for selecting the initial values for the estimation of noncausal and mixed causal-noncausal based on their causal representations. In our second contribution, we present a decision method for identifying the model based on the existence of the global minimum in the estimation function. Third, we propose the third-order spectral estimation of mixed models.
Our estimation method is consistent with strong asymptotic properties without assuming any probability density function. Instead, non-Gaussianity is required. Nonetheless, it may have a higher computational burden than parametric methods, which assume a distribution function to estimate by ML. Although our method is helpful for estimating causal, noncausal, and mixed models, it depends on the initial values obtained in second-order estimates.
The dynamics of commodities depend on their past history but also show a dependence on future information. Nevertheless, the possible dynamics depend on the transformation performed. In the case of the HP filter, the model identified is a MAR($1,1$) in seven cases. On the other hand, in the log-return transformation, the results are more heterogeneous. It identifies an AR($1,0$) twice, AR($0,1$) twice, one AR($0,2$) and three MAR($1,1$). The dynamics are stronger in the HP filter.
Developing a tool to detect dynamics considering higher-order cumulants is interesting for future research. This is because although the estimation is of third-order, second-order criteria are used to detect dynamics, leading to a process that may not reveal all the information in the data. Also consider the existence of the noninvertible component of the moving average model.