{'fraction_non_alphanumeric': 0.04960312947758648, 'fraction_numerical': 0.022447633560008187, 'mean_word_length': 4.5566788828510045, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral estimation for mixed causal-noncausal autoregressive models

This paper investigates new ways of estimating and identifying causal, noncausal, and mixed causal-noncausal autoregressive models driven by a non-Gaussian error sequence. We do not assume any parametric distribution function for the innovations. Instead, we use the information of higher-order cumulants, combining the spectrum and the bispectrum in a minimum distance estimation. We show how to circumvent the nonlinearity of the parameters and the multimodality in the noncausal and mixed models by selecting the appropriate initial values in the estimation. In addition, we propose a method of identification using a simple comparison criterion based on the global minimum of the estimation function. By means of a Monte Carlo study, we find unbiased estimated parameters and a correct identification as the data depart from normality. We propose an empirical application on eight monthly commodity prices, finding noncausal and mixed causal-noncausal dynamics.


Introduction
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 Hansen and Sargent (1991). In addition, ; Lanne, Luoto, and Saikkonen (2012) 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 Breid, Davis, Lh, and Rosenblatt (1991), 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 Gouriéroux and Zakoïan (2017) ;Hecq, Lieb, and Telg (2016); Hecq and Voisin (2021); Lanne and Saikkonen (2011); Lof (2013); Lof and Nyberg (2017).
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 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 Kindop (2021), and for mixed causal-noncausal, see Bec, Nielsen, and Saïdi (2020). 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.

Methodology
Let yt be a time series of length T , with zero mean and finite k-th order moments, for k ≥ 2, and k-th weakly stationary for t = {0, ±1, ±2, ±3, ...}, with autocovariance function κ2(j) = E [ytyt−j] f or j = {0, ±1, ±2, ±3, ...} . (1) Notice that the autocovariance is considering linear relations of the data in two time periods, known as second-order dependence, also that κ2 j = κ2 j−1 ... = κ2 j−n = κ2 j+1 ... = κ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 1. 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.

Figure 1:
In the left is the causal AR(1,0) model yt = 0.7y t−1 + εt, and in the right, the noncausal AR(0,1) model yt = 0.7y t+1 + t, plot and ACF. εt ∼ alpha-stable fa(t; α = 1.5, β = 0.25, γ = 1, δ = 0). The definition of the alpha-stable distribution is in Section 4, the sample size T = 100. The advantage of using the alpha-stable distribution is the ability to have different shapes depending on the value of its parameters. It can easily replicate sample characteristics of leptokurticity and skewness. Although theoretically, the existence of its moments is linked to the values of its stability parameter (α), i.e., the variance, the skewness, and the kurtosis only exist when it is equal to two. An interesting feature of the alpha-stable distribution is its infinite divisibility, i.e., it can be expressed as a sum of i.i.d variables while preserving its distribution. See Section 4 for details.
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 Bartelt, Lohmann, and Wirnitzer (1984). The third-order cumulant is Note that κ3(j, l) is equal to zero for the Gaussian and symmetrical distributions. The definition of κ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 yt, κ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 ω = 2π/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 Λ1 and Λ2. Where ω1 ⊂ Λ1 = (0, 1/2) for S2, and given the symmetries S3(ω1, ω2) = S3(ω2, ω1) = S3(ω1, −ω2 − ω1) = S3(−ω1 − ω2, ω2) and the identity S3(ω1, ω2) = S3(−ω1, −ω2), Λ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 Anh, Leonenko, and Sakhno (2007); Gabr (1989, 2012); Rosenblatt and Van Ness (1965).
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 Kay and Marple (1981). The periodogram is obtained as the modulus of the output values from the DFT in Equation 6, performed directly on the time series. The origin of the periodogram can be traced back to Schuster (1898), 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 Alekseev (1996); Bartlett (1950);Brillinger (1975); Brillinger and Rosenblatt (1967); Rosenblatt (1965). 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 dT (−ω1 − ω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 Alekseev (1996); Bartlett (1950); Brillinger (1975); Rosenblatt (1965). Some properties of the periodogram are developed in Brillinger and Rosenblatt (1967); in particular, the mean of the periodogram k = 2 and biperiodogram k = 3. Denoting ω = ω1, ω2 when ωs = 0 mod 2π. The variance for k = 2, 3 as T → ∞. Equivalently to the periodogram, the biperiodogram at different Fourier frequencies (ω = ω ) is uncorrelated 2.1 The causal AR(r, 0) In this Section, we present causal, noncausal, and mixed autoregressive models. Let yt be a causal autoregressive process generated by φ(L)yt = εt.
The lag polynomial is φ(L) = 1 − φ1L − φ2L 2 ...φrL r , with a set of parameters Φr = (φ1, φ2...φr), where L stands for the backward shift operator, that is, L k = y t−k for k = {0, ±1, ±2, ±3, ...}. 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 2.2 and 2.3, respectively. Note that we differentiate the cumulants of yt, κ2, κ3, with the cumulants of the error sequence εt, κ2 e , κ3 e . φ(z) is provided by the roots 1 − φ1z − φ2z 2 ... − φrz r , all lying outside the unit circle. Hence, the model has a covariance-stationary representation form κ2(j) = E [ytyt−j], and is represented as an infinite sum of past shocks, MA(∞) as where ∞ j=0 |Ψ| < ∞. Then, the autocovariance generating function of yt is gy(z) = ∞ j=−∞ κ2(j)z j . Evaluating gy at z = e −iω , and dividing by 2π we obtain the spectrum of yt 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 ω. For the AR(r, 0) model, the spectrum can be obtained explicitly using the representation of Equation 14. 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 z = e 2πi(−ω 1 −ω 2 ) , the transfer function of the sum of the complex conjugation of the frequencies is In this way, the bispectrum is S3(Φr, ω1, ω2) = k3 e (2π) 2 ψ(Φr, ω1)ψ(Φr, ω2)ψ(Φr, −ω1 − ω2)

The noncausal AR(0, s)
Let yt be a noncausal autoregressive process generated by where the lead polynomial is ϕ(L −1 ) = 1 − ϕ1L −1 − ϕ2L −2 ...ϕsL −s with the set of parameters Φs = (ϕ1, ϕ2...ϕs). In this notation, the AR(0, s) model has s noncausal and zero causal components. The polynomial ϕ(z) is 1 − ϕ1z − ϕ2z 2 ... − ϕsz 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 ϕ(z) has a causal representation with its zeros inside the unit circle, (see Brockwell and Davis (1987) p. 81 for a stationary solution), that is Lanne and Saikkonen (2011), the polynomial ϕ * (z) can be expressed as where ϕ(z −1 ) is as in Equation 19, implying that ϕ * s−j /ϕ * s = −ϕj for j = 1, 2, ...s − 1, and 1/ϕ * s = ϕj. Notice that the zeros of ϕ * (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 ξt = −1/ϕsεt−s, and coefficients for the polynomial ϕ * j equal to ϕ * j = −ϕj/ϕs for j = 1, 2...s − 1, and ϕ * s = 1/ϕ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 yt (see Brockwell and Davis (1987), 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.
Reverting the process in s time periods Solving for yt, we can rewrite yt 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 23 are equivalent to the ones in Equation 21.
The causal representation is key to the estimation and identification of the model, as we shall discuss in Section 3.2.
The transfer function for the noncausal AR(0,s) is ψ(Φs, ω) = 1/(1 − ϕ1z − ϕ2z 2 ... − ϕrz s ), and its spectrum is Additionally, the bispectrum is In Figure 2, we show the same models as in Figure 1, but their frequency domain representation. At the top of the graph is the periodogram and spectrum of the causal AR(1,0) yt = 0.7yt−1 + εt. At the bottom is the periodogram and spectrum of a noncausal AR(0,1) model yt = 0.7yt+1 + εt. εt ∼ alpha-stable fa(t; α = 1.5, β = 0.25, γ = 1, δ = 0), the sample size T = 100. As in Figure 1, the same error sequence is used to construct the models. The periodogram, I2(ω), as well as the spectrums S2(Φrω) = S2(Φsω) 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.

Mixed Cuasal/noncausal MAR(r, s)
Let yt be a mixed causal-noncausal autoregressive process MAR(r, s) generated by where φ(L) = 1 − φ1L − φ2L 2 ...φrL r corresponds to the causal component, and ϕ(L −1 ) = 1 − ϕ1L −1 − ϕ2L −2 ...ϕsL −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 Υr,s = {φ1, ..., φr, ϕ1, ..., ϕs}. The MAR(r, s) is a combination of the variable yt with r lags and s leads, and both polynomials ϕ(z) and φ(z), with roots outside the unit circle. Alternatively, the polynomial ϕ(z) has a causal representation with all its roots inside the unit circle, and the polynomial φ(z) outside the unit circle. For details, see K.-S. Lii and Rosenblatt (1996). Accordingly, yt has a two-sided moving average representation as in Equation 14. The polynomial of the MAR(r, s) model is The MAR(r, s) has a causal purely causal representation AR(r + s,0), as follows This model representation can be found in Gourieroux and Jasiak (2018), Appendix A. Where ξ * t = −1/ϕsεt−s, and υ * 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 υr,s(L) Reverting the process in s time periods, or equivalently multiplying by L s , Solving for yt, The closed form of the AR(r + s,0) is The polynomial υ * r+s (z) , has the set of parameters Υ * r+s = {υ * 1 , υ * 2 , ...υ * s , ...υ * r+s }, then the AR(r + s,0) can be written as Note that there is a one-to-one relation between Equations 26 and 29. 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 yt, we obtain the AR(2, 0) representation The causal representation is where the set of parameters is Υ Estimating the set of parameters υ * is not trivial, given the nonlinearity in the original model. For a similar representation of Equation 32, see Gourieroux and Jasiak (2018) proposition 3.1.
where ψ(Υr,s, ω) In Figure 3, we show the real and imaginary parts of the bispectrum and the biperiodogram of simulated data from a MAR(1,1) with coefficients φ1 = 0.7 and ϕ1 = 0.2 and T = 200, the distribution of the errors is alpha-stable fa(t; α = 1.5, β = 0.25, γ = 1, δ = 0) with leptokurtic and skewness features. The set of frequencies ω1 is in the y-axis and ω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 ω1 and ω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 ω1 and ω2, from 0 to 0.2 and 0.8 to 1 (see the symmetries in Equation 3). 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 ω1 = ω2 diagonal that fades as it approaches the frequency 0.5. This same pattern can be seen in the periodogram and spectrum in Figure 2. 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 ω1 and low (high) frequencies of ω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.

Estimation
In this Section, we present the estimation function that was first proposed by Leonenko et al. (1998). However, we use the variation presented by Velasco and Lobato (2018) in terms of the scaling factor in the denominators of Equation 35. We define ϑ as the set of parameters of the autoregressive model, for the causal ϑ = Φr, for the noncausal ϑ = Φs and for the mixed ϑ = Υr,s. We first perform a preliminary Gaussian likelihood estimation and obtain ϑ, always assuming causality, and we evaluate it in the transfer functions ψ(ϑ, ω) and ψ(ϑ, −ωi − ωj) and the second order cumulant, k2 = 2πT −1 T −1 j=1 I2(ωj)/(ψ(ϑ, ωj)ψ(ϑ, −ω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 Leonenko et al. (1998). The estimation function is where A2T = m(2π) 2 /(4κ2 2 T ) and A3T = n(2π) 4 /(6κ2 3 T 2 ). In all cases, m = 1 − n. The weights in our approach are arbitrarily set as m = 0.5 and n = 0.5. Terdik (1999) estimate an AR(1, 0) with different values for m, highlighting the efficiency gains with the correct selection. Nevertheless, in Velasco and Lobato (2018), 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 Velasco and Lobato (2018), 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 ϑ. The modified spectrum and bispectrum are When the data is Gaussian, the first term of the estimation function RT (ϑ) ( Equation 35) is sufficient to estimate and identify the autoregressive model correctly. The first term in RT (ϑ) is equivalent to Whittle (1953) or Gaussian-likelihood. Nevertheless, when the data are non-Gaussian, Brillinger (1975) suggests including the second term, that is, the bispectrum estimation. This estimator is consistent and asymptotic normal. For details, see Terdik (1999); Velasco and Lobato (2018). The set of optimal parametersθ is obtained by minimizing RT by the quasi-Newton methodθ = arg min ϑ∈Θ RT (ϑ).
Under the conditions mentioned above, and assuming thatθ is close enough to ϑ when T → ∞, the Taylor expansion of RT (θ) is Asθ minimizes RT (ϑ), ∂R T (θ) The closed forms for the first and second derivatives of Equation 41 are Terdik (1999) Section 4. It is reasonable to presume that the numerator in Equation 41 is zero sinceθ is close enough to ϑ, as the third-order cumulant is different from zero and the bispectrum exists.
The parameter η is the Integral respect to ω of the derivative of the modulus of the second-order transfer function with respect the set of parameters ϑ By solving Equation 44 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. Lobato and Velasco (2022) developed closed formulas for the asymptotic variance of the ARMA model allowing for complex roots.

Estimation strategy
Theorem 1 shows that the correct set of parameters can be identified and consistently estimated by finding the global minimum of RT (ϑ). 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θ is close enough to the true value ϑ. The latter entails certain limitations for noncausal and mixed models, given their nonlinearity and that ϑ 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 Kindop (2021), and for mixed causal-noncausal, see Bec et al. (2020). 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 RT (ϑ) since the global minimum must be for AR(2, 0). The problem of identification in the noncausal model has been studied by Andrews et al. (2006); Breid et al. (1991); Lanne and Saikkonen (2011) 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θ to converge in probability to ϑ. Then, we address the multimodality of the estimation function by selecting the model with the lowest RT (θ). The initial parameters are chosen based on the pure causal representation of the noncausal (Equation 21) and mixed models (Equation 27) 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, RT (Φr) refers to the estimation function when the model is purely causal AR(r,0), RT (Φs) when the model is purely noncausal AR(0,s), and RT (Υr,s) when is mixed causal-noncausal MAR(r, s). As a common first step, we perform Gaussian likelihood or Whittle estimation and obtain ϑ and k2.
1. For the estimation of the Causal AR(r, 0), the selection of the initial values to minimize RT (Φr) is ϑ, that is The initial values are always obtained in the preliminary estimation, assuming causality. The minimum obtained by the estimation function is RT (Φr).
We must first conduct a Gaussian-likelihood to obtain the initial values and obtain ϑ. After, the set of parameters Υr,s is reconstructed using the relation between Equations 28 and 29. Example 2 We want to estimate a mixed causal-noncausal MAR(1, 1), the steps are (a) Estimate the causal AR(2, 0) by Gaussian likelihood, and obtainυ * 1 andυ * 2 .
Notice that this nonlinear system has two possible solutions. Both with possible complex roots.
After the estimation, we proceed with the model identification based on the value of the estimation function.

Simulation study
We perform an extensive Monte Carlo study to investigate the finite sample properties of our estimation strategy using RT (ϑ). 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 Fries and Zakoian (2019) and Fries (2021). The probability density function of the alpha-stable is where α ∈ (0, 2) is the stability parameter, β [−1, 1] the skewness parameter, γ > 0 the scale parameter, and δ ∈ the location parameter. We generate random numbers with length T = {100, 200, 500}, for the parameters we choose α = {1.5, 1.8}, β = 0.25, γ = 1 and δ = 0. To allow for skewness and leptokurticity. In Figure 4, 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 φ1 = 0.7, φ2 = 0.2, a purely noncausal AR(0, 2), with ϕ1 = 0.7, ϕ2 = 0.2, and a mixed MAR(1, 1) with φ1 = 0.7, ϕ1 = 0.2.
In Table 1, we present the percentage with which we identify the correct model (identification rate) for the estimation of RT (ϑ). In the first column, there are the parameters of the alpha-stable distribution, notice that the only parameter that varies is α. 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 . 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, α = 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 α = 2.
The mean and standard deviation of the estimated parameters are reported in Table 2, (we only show the results for our estimation strategy).  (2) noncausal AR(2) 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, α = 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 α = 1.8, revealing that the variance of the estimator depends on the sample moments.

Empirical application
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 (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. Hecq and Voisin (2021) found evidence of noncausality in monthly Nickel prices, Gourieroux, Jasiak, and Tong (2021); Hecq and Voisin (2019) in crude oil monthly prices, Karapanagiotidis (2014) in 25 commodity futures price, including soft, precious metals, energy, and livestock sectors and Lof and Nyberg (2017) 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 7 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, Hencic and Gouriéroux (2015) use a cubic trend, Fries and Zakoian (2019) use a linear trend or a simple centering, Cubadda, Hecq, and Telg (2019);Hecq, Issler, and Telg (2020) use the returns,  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, Hecq and Voisin (2019) proposes to use the Hodrick-Prescott (HP hereafter) filter with a penalty parameter λ =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 Hamilton (2018) for details. For the sake of discussion, we consider two transformations, the HP filter with a penalty parameter λ = 129.600, using the cyclical component, and the logarithmic returns. The Figures of these two transformations are respectively 8, 9.
Following Breid et al. (1991), we estimate all series by Gaussian likelihood and selected the order p by secondorder 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 3.
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 4, we report the estimation using RT (ϑ) 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 4 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 , 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 Lanne and Saikkonen (2011). In Table 5, we present the estimation using RT (ϑ) 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 5. 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 6, we present the graphical representation of the estimation function RT (Υp,s) for the cyclical component of the Brent oil. On the x-axis is the noncausal coefficient (ϕ1), and on the y-axis is the causal coefficient (φ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 ϕ1 = 0.811 and φ1 = 0.459, located at the intersection of the vertical (green) and horizontal (red) lines.

Conclusions and Discussion
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 secondorder 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.