A New MARMA-GARCH Model with Applications to Cryptocurrency Volatility

1. Introduction

Financial returns often exhibit weak autocorrelation but strong and persistent volatility, posing challenges for traditional models to achieve independent errors. The mixed causal-noncausal invertible-noninvertible autoregressive moving average (MARMA) model emerges as a generalization of the ARMA model, allowing the autoregressive and moving average polynomials to have roots outside but also inside the unit disk. Whereas the non-invertible process is well known, a more recent literature studies the presence of a non-causal component that can generate local explosions commonly observed in financial asset prices, motivating their use in economic applications.1

A key feature, and one of the most challenging aspects, of a MARMA process is that it admits multiple ARMA representations, each sharing the same second-order probability structure. This is advantageous for estimation, as all representations possess an identical variance-covariance matrix and, thus, the same spectral density. However, this multiplicity poses a challenge for identification. Roots inside and outside the unit disk result in the same values of the estimation function under Gaussian estimation methods like least squares, Yule-Walker, or maximum likelihood. This is the case with the most common ARMA representation of the MARMA process, which assumes that the roots are outside the unit disk, known as the causal and invertible ARMA representation [12]. Identification becomes feasible when the assumption of normality is relaxed, each representation having a unique and identifiable higher-order probability structure. These differences are evident in the autocorrelation functions of the squared residuals and the higher-order spectral densities, revealing significant dynamics.

The MARMA model can generate nonlinear dynamics similar to those produced by the Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model [8]. However, the MARMA model has a distinct probability structure, incorporating second-, third-, and fourth-order relationships. [39] show that the nonlinear dynamics inherent in the MARMA model differ significantly from those of the GARCH model due to the involvement of third and fourth-order cumulants and the interactions between errors and future values of the time series. Furthermore, [12] noted that although all-pass models can produce nonlinear dynamics, their high-order probability structure is not sufficient to compete with the GARCH model, as they cannot replicate the well-known clustering effect, where small variations are followed by small variations and large variations are followed by large variations [22].

This limitation is particularly evident in financial returns that exhibit weak autocorrelation but extreme and persistent volatility, where no MARMA model can achieve independent errors. Conversely, the squared errors display significant dynamics, as seen in the log-returns of major cryptocurrencies [51], stocks, commodities, and market indexes [16]. To replicate this characteristic, it becomes a natural option to relax the assumption of i.i.d. errors and to introduce a more flexible form that allows for time-varying volatility. This approach implies adding a GARCH process to the errors. Such models could provide a better fit for financial returns, especially during periods of high volatility and extreme values, as observed during recent financial crises.

This new approach presents challenges in estimation and identification, as it departs from the traditional assumption of i.i.d. errors in the MARMA model. While the theory surrounding the estimation of the MARMA model typically assumes that the sequence of errors follows an i.i.d. process, [47] consider maximum likelihood estimation of non-invertible ARMA models with ARCH errors, providing conditions under which a strongly consistent and asymptotically normal distributed solution to the likelihood equations exists under a finite variance. We generalize [47] to the estimation and identification of the MARMA-GARCH model, which can be not only noninvertible but also noncausal. We investigate the use of the Whittle estimation in the frequency domain as well as the maximum likelihood. In both methods we first estimate the MARMA model and compute the error process on which we estimate the GARCH parameters. This thorough comparison is valuable because both methods are equivalent in estimating the MARMA process, sharing similar characteristics of consistency and asymptotic normality. For example, both require a stationary environment and finite variance, and they perform well with heavy-tailed errors. However, for the GARCH processes, maximum likelihood estimation only requires the existence of the variance of the process, while the Whittle estimation requires the existence of eight moments, representing a limitation for financial returns that often exhibit extreme tails. Therefore, we compare the performance of Whittle and maximum likelihood estimation methods, providing a comprehensive understanding of their strengths and limitations. Our Monte Carlo study reveals that the MLE produces parameter estimates with lower bias and less standard deviation.

The paper is organized as follows: In Section 2, we present the MARMA-GARCH model along with the methods for estimation and identification and the techniques for detrending the data. Section 3 provide the results of the Monte Carlo simulations, and the empirical application to cryptocurrencies. Finally, we conclude with a discussion and conclusion in Section 4.

2. Materials and Methods

2.1. MARMA($r,s,r^{\prime },s^{\prime }$)-GARCH($p,q$) Processes

Let $\left\{Y_t\right\}$ be a univariate mixed causal-noncausal invertible-noninvertible autoregressive moving-average generalized autoregressive conditional heteroskedasticity MARMA ($r,s,r^{\prime },s^{\prime }$)-GARCH($p,q$) process that is stationary for every $t=\{0,\pm 1,\pm 2,\dots \}$, and represented by:

$${\varphi }^{+}\left(L\right){\varphi }^{*}\left(L^{-1}\right)Y_t={\theta }^{+}\left(L\right){\theta }^{*}\left(L^{-1}\right){\epsilon }_t,$$

(1)

$${\epsilon }_t={\eta }_t{\sigma }_t,{\eta }_t\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathfrak{F}\left(\zeta \right),$$

(2)

$${\sigma }_t^2=\omega +\alpha \left(L\right){\epsilon }_t^2+\beta \left(L\right){\sigma }_t^2,$$

(3)

The sequence of innovations $\left\{{\eta }_t\right\}$ is assumed to be i.i.d. with zero mean and unit variance, and $\mathfrak{F}\left(\zeta \right)$ is a probability density function with $\zeta$ as its parameters.

$$\mathbb{E}\left[{\eta }_t\right]=0,\mathbb{E}\left[{\eta }_t^2\right]=1,\mathrm{and}\mathrm{finite}\mathrm{fourth}-\mathrm{order}\mathrm{moments}\text{:}\mathbb{E}\left[{\eta }_t^4\right]<\infty$$

where the sequence $\left\{{\epsilon }_t\right\}$ is a martingale difference with moments and cumulants:

$$\mathbb{E}\left[{\epsilon }_t\right]=0,\mathbb{E}\left[{\epsilon }_t^2\right]={\sigma }_t^2,\mathbb{E}\left[{\epsilon }_t^3\right]={\kappa }_3,\mathbb{E}\left[{\epsilon }_t^4\right]={\mu }_4,\mathrm{and}{\kappa }_4={\mu }_4-3{\sigma }_t^4.$$

Let $L^k$ denotes the backward shift operator for $k=\{0,\pm 1,\pm 2,\dots \}$. The lag polynomials ${\varphi }^{+}\left(L\right)$ and ${\theta }^{+}\left(L\right)$, and the lead polynomials ${\varphi }^{*}\left(L^{-1}\right)$ and ${\theta }^{*}\left(L^{-1}\right)$ are defined as:

$$\begin{array}{cc}\hfill {\varphi }^{+}\left(L\right)& =1-\sum _{j=1}^r{\varphi }_j^{+}{L}^j,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\varphi }^{*}\left(L^{-1}\right)& =1-\sum _{k=1}^s{\varphi }_k^{*}{L}^{-k},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\theta }^{+}\left(L\right)& =1+\sum _{l=1}^{r^{\prime }}{\theta }_l^{+}{L}^l,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\theta }^{*}\left(L^{-1}\right)& =1+\sum _{m=1}^{s^{\prime }}{\theta }_m^{*}{L}^{-m}.\hfill \end{array}$$

(7)

We define the MARMA($r,s,r^{\prime },s^{\prime }$) parameter vector as

$$\vartheta =({\varphi }_1^{+},\dots ,{\varphi }_r^{+},{\varphi }_1^{*},\dots ,{\varphi }_s^{*},{\theta }_1^{+},\dots ,{\theta }_{r^{\prime }}^{+},{\theta }_1^{*},\dots ,{\theta }_{s^{\prime }}^{*})\in {\mathbb{C}}^{r+s+r^{\prime }+s^{\prime }}.$$

Note that for $s=0$ and $s^{\prime }=0$, the MARMA($r,s,r^{\prime },s^{\prime }$) model becomes the causal and invertible ARMA($p,q$) model, under the assumption that the dependency structure is exclusively determined by lagged terms. On the other hand, the GARCH polynomials are defined as:

$$\begin{array}{cc}\hfill \alpha \left(L\right)& =\sum _{n=1}^p{\alpha }_n{L}^n,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \beta \left(L\right)& =\sum _{m=1}^q{\beta }_m{L}^m.\hfill \end{array}$$

(9)

For $p\ge 0$, $q>0$, $\omega >0$, also ${\alpha }_i>0$, and ${\beta }_j>0$ ensures the non-negativity of the GARCH process ${\sigma }_t^2$. The parameter vector of the GARCH($p,q$) is defined as

$$\xi =(\omega ,{\alpha }_1,\dots ,{\alpha }_p,{\beta }_1,\dots ,{\beta }_q)\in {\mathbb{R}}^{p+q}.$$

Equation 1 can be expanded as follows:

$$\begin{array}{c}\hfill \left(1-\sum _{j=1}^r{\varphi }_j^{+}{L}^j-\sum _{k=1}^s{\varphi }_k^{*}{L}^{-k}+\sum _{j=1}^r\sum _{k=1}^s{\varphi }_j^{+}{\varphi }_k^{*}{L}^{j-k}\right)Y_t=\left(1+\sum _{l=1}^{r^{\prime }}{\theta }_l^{+}{L}^l+\sum _{m=1}^{s^{\prime }}{\theta }_m^{*}{L}^{-m}+\sum _{l=1}^{r^{\prime }}\sum _{m=1}^{s^{\prime }}{\theta }_l^{+}{\theta }_m^{*}{L}^{l-m}\right){\epsilon }_t.\end{array}$$

(10)

We assume that all common roots between ${\varphi }^{+}\left(z\right){\varphi }^{*}\left(z\right)$ and ${\theta }^{+}\left(z\right){\theta }^{*}\left(z\right)$ are canceled. Moreover, the polynomials ${\varphi }^{+}\left(z\right)$, ${\varphi }^{*}\left(z\right)$, ${\theta }^{*}\left(z\right)$ and ${\theta }^{+}\left(z\right)$ have their zeros outside the unit disk. For the GARCH polynomials, we assume no common roots and roots strictly outside the unit disk. The polynomials can be expressed in terms of their roots:

$$\begin{array}{cc}\hfill {\varphi }^{+}\left(z\right)& =1-\sum _{j=1}^r{\varphi }_j^{+}{z}^j=0\mathrm{for}\left|z\right|>1,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\varphi }^{*}\left(z\right)& =1-\sum _{k=1}^s{\varphi }_j^{*}{z}^j=0\mathrm{for}\left|z\right|>1,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\theta }^{+}\left(z\right)& =1+\sum _{l=1}^{r^{\prime }}{\theta }_j^{+}{z}^j=0\mathrm{for}\left|z\right|>1,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\theta }^{*}\left(z\right)& =1+\sum _{m=1}^{s^{\prime }}{\theta }_j^{*}{z}^j=0\mathrm{for}\left|z\right|>1.\hfill \end{array}$$

(14)

Note that the MARMA inverse polynomials admit individual Laurent expansions. Consequently, the MARMA($r,s,r^{\prime },s^{\prime }$) model admits a two-sided representation, an MA(∞), and an AR(∞) as discussed in p. 88 [14]. The MA(∞) representation is given by:

$$Y_t=\sum _{j=-\infty }^{\infty }{\Psi }_j{\epsilon }_{t+j},\mathrm{where}{\displaystyle \frac{{\theta }^{+}\left(z\right){\theta }^{*}\left(z\right)}{{\varphi }^{+}\left(z\right){\varphi }^{*}\left(z\right)}}=\sum _{j=-\infty }^{\infty }{\Psi }_j{z}^j=\Psi \left(z\right),$$

(15)

where ${\Psi }_j$ are the coefficients of $z^j$ in the Laurent expansion of ${\displaystyle \frac{{\theta }^{+}\left(z\right){\theta }^{*}\left(z\right)}{{\varphi }^{+}\left(z\right){\varphi }^{*}\left(z\right)}}\stackrel{\mathrm{def}}{=}\Psi \left(z\right)$. Note that $\Psi \left(z\right)$ is absolutely convergent in the annulus $D=\left\{z:\delta <\left|z\right|<{\delta }^{-1}\right\},\delta <1$. Therefore, $m_1,\cdots ,m_r,m_{r+1},\cdots ,m_{s+r}$ are the zeros of the polynomials ${\varphi }^{+}\left(z\right)$ and ${\varphi }^{*}\left(z\right)$. Similarly, $n_1,\cdots ,n_{r^{\prime }},n_{r^{\prime }+1},\cdots ,n_{r^{\prime }+s^{\prime }}$ are the zeros of the polynomials ${\theta }^{+}\left(z\right)$ and ${\theta }^{*}\left(z\right)$. The polynomials can be factorized as follows:

$$\begin{array}{cc}\hfill {\varphi }^{+}\left(z\right)& =\prod _{j=1}^r(1-m_{j}z),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\varphi }^{*}\left(z\right)& =\prod _{j=r+1}^s(1-m_j{}_z),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\theta }^{+}\left(z\right)& =\prod _{j=1}^{r^{\prime }}(1-n_{j}z),\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\theta }^{*}\left(z\right)& =\prod _{j=r^{\prime }+1}^{s^{\prime }}(1-n_j{}_z).\hfill \end{array}$$

(19)

Note that in Equations (16-), any root z can be replaced by its inverse $z^{-1}$ without altering the second-order moments of the MARMA process. Indeed, there exist $2^{r+s+r^{\prime }+s^{\prime }}$ possible MARMA models to represent the $\left\{Y_t\right\}$ process. Each model can produce uncorrelated but not independent errors in the finite variance setting since they share the same second-order moments. Consequently, the MARMA($r,s,r^{\prime },s^{\prime }$) and the ARMA($p,q$) processes have the same ACF and second-order spectral density. As a result, the two processes are not identifiable under the assumption of Gaussianity. On the contrary, under a non-Gaussian distribution, each process has a different and identifiable probability structure that relies on the high-order moments.

According to [22], a squared GARCH($p,q$) model can be reformulated as an ARMA($max(p,q),q$) model. Thus, Equation (3) can be rewritten as:

$${\epsilon }_t^2-\sum _{i=1}^{\max (p,q)}({\alpha }_i+{\beta }_i){\epsilon }_{t-i}^2=\omega +v_t-\sum _{j=1}^q{\beta }_j{v}_{t-j},$$

(20)

where $v_t={\epsilon }_t^2-{\sigma }_t^2={\sigma }_t^2({\eta }_t^2-1)$. We assume that ${\sigma }_t^2$ is strictly stationary, with $\mathbb{E}\left[{\epsilon }_1^2\right]<\infty$, therefore, $\left(v_t\right)$ are martingale differences. The RHS of Equation 20 can be factorized in terms of the polynomials:

$$(\alpha +\beta )\left(z\right)=1-({\alpha }_1+{\beta }_1)z-\dots -({\alpha }_{\max (p,q)}+{\beta }_{\max (p,q)})z^{\max (p,q)},$$

(21)

and

$$\beta \left(z\right)=1-{\beta }_{1}z-\dots -{\beta }_q{z}^q.$$

(22)

We can observe that, as in Equation 16, the roots of the polynomials of the GARCH process can also lie outside the unit disk without altering the second-order probability structure of the model. As a consequence, the squared errors (${\epsilon }_t^2$) can be expressed as an ARMA($\max (p,q),q$) with martingale differences ($v_t$):

$$(\alpha +\beta )\left(L\right){\epsilon }_t^2=\omega +\beta \left(L\right)v_t.$$

(23)

The equivalence between Equations (3) and (20) was proposed by Bollerslev (1986). This equivalence allows us to estimate the parameters in Equation (23) using the properties of ${\epsilon }_t$, which is uncorrelated but not independent. In this context, the parameters of the GARCH model can be estimated using standard procedures such as Gaussian maximum likelihood and least squares. However, it is important to note that $\left(v_t\right)$ is uncorrelated but not independent and identically distributed (i.i.d.) since it depends on the past values of ${\epsilon }_t$. Therefore, this representation is called the squared of a GARCH process.

2.1.1. Causal and Invertible Representation

A MARMA($r,s,r^{\prime },s^{\prime }$) process can be represented in multiple ARMA($p,q$) forms, where $p=r+s$ and $q=r^{\prime }+s^{\prime }$, all maintaining the same second-order dependency structure [39]. One such form is the causal and invertible representation,2 characterized by polynomial roots outside the unit disk. For noncausal and noninvertible data-generating processes, this representation yields errors that are uncorrelated but not independent ([12]). However, estimating the invertible all-pass model allows for identifying the roots of the noncausal and noninvertible polynomials through the factorization in Equation 16. The causal and invertible representation is described by:

$$\varphi \left(L\right)Y_t=\theta \left(L\right){\epsilon }_t^{*},$$

with ${\epsilon }_t^{*}$ being uncorrelated but not independent, and where $\varphi \left(z\right)$ and $\theta \left(z\right)$ are polynomials of order p and q, respectively, and only depending on lagged terms. Therefore, Equation 10 can be expressed as:

$$\begin{array}{c}\hfill \left(1-\sum _{j=1}^r{\varphi }_j^{+}{L}^j-\sum _{i=1}^s{\varphi }_i^{*}{L}^i+\sum _{j=1}^r\sum _{i=1}^s{\varphi }_j^{+}{\varphi }_i^{*}{L}^{j+i}\right)Y_t=\left(1+\sum _{j=1}^{r^{\prime }}{\theta }_j^{+}{L}^j+\sum _{i=1}^{s^{\prime }}{\theta }_i^{*}{L}^i+\sum _{j=1}^{r^{\prime }}\sum _{i=1}^{s^{\prime }}{\theta }_j^{+}{\theta }_i^{*}{L}^{j+i}\right){\epsilon }_t^{*}.\end{array}$$

(24)

Considering that the coefficients ${\varphi }_j^{+}$ and ${\varphi }_i^{*}$ are invariant to the allocation of the roots, we establish $\varphi \left(L\right)={\varphi }^{+}\left(L\right){\varphi }^{*}\left(L\right)$, and $\theta \left(L\right)={\theta }^{+}\left(L\right){\theta }^{*}\left(L\right)$. The causal and invertible representation of the MARMA($r,s,r^{\prime },s^{\prime }$) process is:

$$\left(1-\sum _{j=1}^{r+s}{\varphi }_j{L}^j\right)Y_t=\left(1+\sum _{i=1}^{r^{\prime }+s^{\prime }}{\theta }_j{L}^{-i}\right){\epsilon }_t^{*}.$$

(25)

The set of parameters is $\beta =({\varphi }_1,\dots ,{\varphi }_{r+s},{\theta }_1,\dots ,{\theta }_{r^{\prime }+s^{\prime }})\in {\mathbb{R}}^{p+q}$. The representation is causal since there exists a sequence of constants $\left\{{\pi }_j\right\}$, such that $Y_t={\sum }_{j=0}^{\infty }{\pi }_j{\epsilon }_{t-j}^{*}$, with ${\sum }_{j=0}^{\infty }\left|{\pi }_j\right|<\infty$ and $\left\{Y_t\right\}$ expressed in terms of the present and past of the error sequence, which is an MA(∞). On the other hand, it is invertible since there exists a sequence of constants $\left\{{\pi }_j^{*}\right\}$, such that ${\epsilon }_t^{*}={\sum }_{j=0}^{\infty }{\pi }_j^{*}{Y}_{t-j}$, with ${\sum }_{j=0}^{\infty }\left|{\pi }_j^{*}\right|<\infty$. Hence, ${\epsilon }_t^{*}$ is formulated in terms of both the present and historical values of $Y_t$, characterizing an AR(∞) process.

The causal and invertible representation is typically obtained using traditional estimation methods. However, there is a risk that misspecification of the roots may lead to residuals that are uncorrelated but not independent, especially showing dynamics in the squared residuals that differ from those of the GARCH process (see [39], for details). This is why we expect that the errors arising from the MARMA process and its causal and invertible representation will differ regarding the GARCH estimated parameters and the higher-order dynamics.

2.1.2. Spectral Representations

The spectral representation of the model is divided into two parts: the MARMA part and the GARCH part. The spectral density of the MARMA model was proposed by [39], where it is shown to be equivalent to that of the MARMA($p,q$) model in second-order moments. The spectral density of the GARCH part is generated through its ARMA($max(p,q),q$) representation with conditionally heteroskedastic innovations proposed by [8].

Let $\left\{Y_t\right\}$ be a MARMA($r,s,r^{\prime },s^{\prime }$)-GARCH($p,q$) process satisfying Equation 1. The transfer function for the conditional mean is:

$$\psi (\vartheta ,\omega )={\displaystyle \frac{{\theta }^{+}\left(z\right){\theta }^{*}\left(z\right)}{{\varphi }^{+}\left(z\right){\varphi }^{*}\left(z\right)}}={\displaystyle \frac{1+{\sum }_{j=1}^r{\theta }_j^{+}{z}^j+{\sum }_{i=1}^s{\theta }_i^{*}{z}^{-i}+{\sum }_{j=1}^r{\sum }_{i=1}^s{\theta }_j^{+}{\theta }_i^{*}{z}^{j-i}}{1-{\sum }_{j=1}^{r^{\prime }}{\varphi }_j^{+}{z}^j-{\sum }_{i=1}^{s^{\prime }}{\varphi }_i^{*}{z}^{-i}+{\sum }_{j=1}^{r^{\prime }}{\sum }_{i=1}^{s^{\prime }}{\varphi }_j^{+}{\varphi }_i^{*}{z}^{j-i}}}.$$

(26)

The k-th order spectral density is given by:

$$S_k(\vartheta ,{\kappa }_k;{\omega }_1,\dots ,{\omega }_{k-1})={\displaystyle \frac{{\kappa }_k}{{\left(2\pi \right)}^{k-1}}}\left(\prod _{j=1}^{k-1}\psi (\vartheta ,{\omega }_j)\right)\psi \left(\vartheta ,\sum _{j=1}^{k-1}-{\omega }_j\right),$$

(27)

where ${\kappa }_k$ is the k-th order cumulant, ${\omega }_{\ell }$ are the Fourier frequencies ${\displaystyle \frac{2\pi \ell }{T}}$ for $\ell =1,\dots ,T-1$, where T is the length of the data, and $z=e^{-i\omega }$ is a root represented by a complex number. Following Equation 27, the second-order ($k=2$), third-order ($k=3$), and fourth-order ($k=4$) spectral densities are:

$$\begin{array}{cc}\hfill S_2(\vartheta ,{\kappa }_2,\omega )& ={\displaystyle \frac{{\kappa }_2}{2\pi }}\psi (\vartheta ,\omega )\psi (\vartheta ,-\omega ),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill S_3(\vartheta ,{\kappa }_3,{\omega }_1,{\omega }_2)& ={\displaystyle \frac{{\kappa }_3}{4{\pi }^2}}\psi (\vartheta ,{\omega }_1)\psi (\vartheta ,{\omega }_2)\psi (\vartheta ,-{\omega }_1-{\omega }_2),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill S_4(\vartheta ,{\kappa }_4,{\omega }_1,{\omega }_2,{\omega }_3)& ={\displaystyle \frac{{\kappa }_4}{8{\pi }^3}}\psi (\vartheta ,{\omega }_1)\psi (\vartheta ,{\omega }_2)\psi (\vartheta ,{\omega }_3)\psi (\vartheta ,-{\omega }_1-{\omega }_2-{\omega }_3).\hfill \end{array}$$

(30)

On the other hand, the transfer function for the conditional variance is defined as:

$$g(\xi ,\omega )={\displaystyle \frac{\beta \left(z\right)}{(\alpha +\beta )\left(z\right)}}={\displaystyle \frac{1-{\sum }_{j=1}^q{\beta }_j{z}^j}{1-{\sum }_{i=1}^{\max (p,q)}({\alpha }_i+{\beta }_i)z^i}}.$$

(31)

The k-th order spectral density is given by:

$$S_k^v(\xi ,{\kappa }_k^v;{\omega }_1,\dots ,{\omega }_{k-1})={\displaystyle \frac{{\kappa }_k^v}{{\left(2\pi \right)}^{k-1}}}\left(\prod _{j=1}^{k-1}g(\xi ,{\omega }_j)\right)g\left(\xi ,\sum _{j=1}^{k-1}-{\omega }_j\right),$$

(32)

where ${\kappa }_k^v$ is the k-th order cumulant of ($v_t$). In the same reasoning as in the MARMA part, the spectral densities for $k=2,3,4$ are:

$$\begin{array}{cc}\hfill S_2^v(\xi ,{\kappa }_2^v,\omega )& ={\displaystyle \frac{{\kappa }_2^v}{2\pi }}g(\xi ,\omega )g(\xi ,-\omega ),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill S_3^v(\xi ,{\kappa }_3^v,{\omega }_1,{\omega }_2)& ={\displaystyle \frac{{\kappa }_3^v}{4{\pi }^2}}g(\xi ,{\omega }_1)g(\xi ,{\omega }_2)g(\xi ,-{\omega }_1-{\omega }_2),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill S_4^v(\xi ,{\kappa }_4^v,{\omega }_1,{\omega }_2,{\omega }_3)& ={\displaystyle \frac{{\kappa }_4^v}{8{\pi }^3}}g(\xi ,{\omega }_1)g(\xi ,{\omega }_2)g(\xi ,{\omega }_3)g(\xi ,-{\omega }_1-{\omega }_2-{\omega }_3).\hfill \end{array}$$

(35)

The k-th order periodogram gives the sample counterpart of the k-th order spectral density:

$$I_k({\omega }_1,\dots ,{\omega }_{k-1})={\displaystyle \frac{1}{{\left(2\pi \right)}^{k-1}T}}\left(\prod _{j=1}^{k-1}{d}_T\left({\omega }_j\right)\right)d_T\left(-\sum _{j=1}^{k-1}{\omega }_j\right),$$

(36)

where $d_T\left(\omega \right)={\sum }_{t=1}^{T-1}{Y}_t{e}^{-it\omega }$ is the discrete Fourier transform. For $k=2$, Equation 36 is the periodogram, for $k=3$ is the biperiodogram, and for $k=4$ is the triperiodogram:

$$\begin{array}{cc}\hfill I_2\left(\omega \right)& ={\displaystyle \frac{1}{2\pi T}}{d}_T\left(\omega \right)d_T(-\omega ),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill I_3({\omega }_1,{\omega }_2)& ={\displaystyle \frac{1}{{\left(2\pi \right)}^{2}T}}{d}_T\left({\omega }_1\right)d_T\left({\omega }_2\right)\overline{d_T(-{\omega }_1-{\omega }_2)},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill I_4({\omega }_1,{\omega }_2,{\omega }_3)& ={\displaystyle \frac{1}{{\left(2\pi \right)}^{3}T}}{d}_T\left({\omega }_1\right)d_T\left({\omega }_2\right)d_T\left({\omega }_3\right)\overline{d_T(-{\omega }_1-{\omega }_2-{\omega }_3)},\hfill \end{array}$$

(39)

where $\overline{d_T(-{\omega }_1-{\omega }_2)}$ and $\overline{d_T(-{\omega }_1-{\omega }_2-{\omega }_3)}$ are the complex conjugates of the Fourier transform evaluated at the sum of two and three frequencies, respectively.

2.2. Parameter Estimation and Identification

We propose a two-step estimation procedure. First, we estimate the parameters of the MARMA($r,s,r^{\prime },s^{\prime }$) model using the causal and invertible ARMA($p,q$) representation. We select the order p and q based on information criteria (BIC and AIC) and the significance of coefficients, preferring the most parsimonious form. Then, we factorize the AR and MA polynomials to explore all combinations of roots that result in noncausal and noninvertible models. To select the best-fitting model, we choose the one that produces residuals with the least high-order dependence. We use a high-order spectral function proposed by [50] and a Portmanteau-type test on the residuals and squared residuals.

In the second step, we estimate the GARCH($p,q$) model, again using information criteria to select p and q. However, it is commonly known that GARCH(1,1) often adequately represents the conditional volatility dynamics in financial assets. We compare the Whittle and time-domain maximum likelihood estimations to determine which method better captures the high-order dynamics.

2.2.1. Estimation and Identification of the MARMA($r,s,r^{\prime },s^{\prime }$) Process

We consider the Whittle estimation of the causal and invertible ARMA($p,q$) model based on the observations $Y_1,\dots ,Y_T$ of the MARMA($r,s,r^{\prime },s^{\prime }$) process. Minimizing the Whittle estimator is asymptotically equivalent to the Gaussian-likelihood given that the distribution of the Fourier transforms of $\left\{Y_t\right\}$ converges in distribution to independent complex normal variables. Therefore, the periodogram tends in probability to the second-order spectral density.

The causal and invertible representation satisfies the usual Yule-Walker equations as the ACF of its residuals converges to the ACF of the MARMA($r,s,r^{\prime },s^{\prime }$) model, also as they have identical second-order spectral density. Consequently, it can be consistently estimated using the Whittle function; see [25] for a detailed overview. The Whittle estimator of ${\beta }_0=\left\{{\varphi }_{0,1},\dots ,{\varphi }_{0,r+s},{\theta }_{0,1},\dots ,{\theta }_{0,r^{\prime }+s^{\prime }}\right\}$ is:

$$W(\beta ,\omega )=\sum _{\ell =1}^{T/2}log\left(S_2(\beta ,{\kappa }_2,\omega )\right)+\sum _{\ell =1}^{T/2}{\displaystyle \frac{I_2\left(\omega \right)}{S_2(\beta ,{\kappa }_2,\omega )}}.$$

(40)

The estimated parameter set is obtained by minimizing the function:

$$\widehat{\beta }=\underset{\beta \in \mathbf{B}}{arg\; min}W(\beta ,{\kappa }_2,\omega ).$$

(41)

Once obtained the set of estimated parameters $\widehat{\beta }$ of the causal and invertible ARMA($p,q$) representation, we proceed to construct all possible $2^{r+s+r^{\prime }+s^{\prime }}$ representations, restricted to the order $p=r+s$ and $q=r^{\prime }+s^{\prime }$. For this purpose, we factorize the polynomials $\varphi \left(z\right)$ and $\theta \left(z\right)$. Second, we distribute all the roots using Equation 16, to construct polynomials ${\varphi }^{+}\left(z\right){\varphi }^{*}\left(z\right)$ and ${\theta }^{+}\left(z\right){\theta }^{*}\left(z\right)$.

The global minimum is achieved if $W(\widehat{\beta },\omega )\in \mathbf{B}\subset (-1,1)$, and $W\left(\widehat{\beta }\right)\le W\left(\beta \right),\forall \beta \in \mathbf{B}$. [35] demonstrated that the Whittle estimator $\widehat{\beta }$ is consistent: $\widehat{\beta }\stackrel{p}{⟶}{\beta }_0$, as $T\to \infty$. The asymptotic distribution of the estimator is Gaussian: $\sqrt{T}(\widehat{\beta }-{\beta }_0)\stackrel{d}{⟶}N(0,\mathbf{V}\left({\beta }_0\right))$. On the other hand, the asymptotic variance-covariance matrix of the estimator $\mathbf{V}\left({\beta }_0\right)$ can be approximated by the inverse of the Fisher information matrix: $\mathbf{V}\left({\beta }_0\right)\approx {\mathbf{I}}^{-1}\left({\beta }_0\right)$. The Fisher information criteria are computed in the frequency domain:

$$\mathbf{I}\left({\beta }_0\right)={\displaystyle \frac{1}{4\pi }}{\int }_{-\pi }^{\pi }{\displaystyle \frac{1}{S_2{({\beta }_0,{\kappa }_2,\omega )}^2}}\left({\displaystyle \frac{\partial S_2({\beta }_0,{\kappa }_2,\omega )}{\partial \beta }}\right){\left({\displaystyle \frac{\partial S_2({\beta }_0,{\kappa }_2,\omega )}{\partial \beta }}\right)}^{\prime }d\omega .$$

(42)

Notice that the Whittle estimator convergence to its true value occurs when $\mathbb{E}\left[Y_t^4\right]<\infty$, which ensures that the periodogram converges to the second-order spectral density ([49]). In a context of infinite variance, the mean of the periodogram cannot converge to the spectral density, leading to biased parameters and invalid inference.

The Hessian matrix $H_T\left(\beta \right)$ is defined as the second-order partial derivatives of the Whittle function with respect to the parameter vector $\beta =\{{\varphi }_{0,1},\dots ,{\varphi }_{0,r+s},{\theta }_{0,1},\dots ,{\theta }_{0,r^{\prime }+s^{\prime }}\}$:

$$\mathbf{H}\left(\beta \right)=\left[\begin{array}{ccc}{\displaystyle \frac{{\partial }^{2}W(\beta ,{\kappa }_2,\omega )}{\partial {\varphi }_{0,1}^2}}& \cdots & {\displaystyle \frac{{\partial }^{2}W(\beta ,{\kappa }_2,\omega )}{\partial {\varphi }_{0,1}\partial {\theta }_{0,r^{\prime }+s^{\prime }}}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{{\partial }^{2}W(\beta ,{\kappa }_2,\omega )}{\partial {\theta }_{0,r^{\prime }+s^{\prime }}\partial {\varphi }_{0,1}}}& \cdots & {\displaystyle \frac{{\partial }^{2}W(\beta ,{\kappa }_2,\omega )}{\partial {\theta }_{0,r^{\prime }+s^{\prime }}^2}}\end{array}\right]$$

2.2.2. Spectral Identification Function

[50] propose to combine the third- and fourth-order spectral densities into a single function, which captures the information about the respective order dynamics. These dynamics align with those observed in the residuals of a misspecified model, as illustrated in [39] Example 1. The combined spectral function is:

$$L_{kT}\left(\vartheta \right)={\displaystyle \frac{{\left(2\pi \right)}^{k-2}}{2k{T}^{k-1}}}\sum _{\mathbf{j}=1}^{T-1}{\displaystyle \frac{{\left|I_k({\omega }_1,\dots ,{\omega }_{k-1})-S_k(\vartheta ,{\kappa }_k\left(\vartheta \right);{\omega }_1,\dots ,{\omega }_{k-1})\right|}^2}{{\left|\left({\prod }_{j=1}^{k-1}\psi (\beta ,{\omega }_j)\right)\psi \left(\beta ,{\sum }_{j=1}^{k-1}-{\omega }_j\right)\right|}^2}},$$

(43)

where $\mathbf{j}=j_1,\dots ,j_{k-1}$. The denominator represents the k-th order transfer function, evaluated using the causal and invertible model parameters. When high-order dynamics are absent, the discrepancy between the k-th order periodogram and the spectral density follows a multivariate normal distribution:

$${\displaystyle \frac{I_k({\omega }_1,\dots ,{\omega }_{k-1})-S_k(\vartheta ,{\kappa }_k\left(\vartheta \right);{\omega }_1,\dots ,{\omega }_{k-1})}{\left({\prod }_{j=1}^{k-1}\psi (\beta ,{\omega }_j)\right)\psi \left(\beta ,{\sum }_{j=1}^{k-1}-{\omega }_j\right)}}\to N(0,I_k).$$

(44)

For $k=3$ and $k=4$, the spectral identification function is:

$$L_T\left(\vartheta \right)=L_{3T}\left(\vartheta \right)+L_{4T}\left(\vartheta \right).$$

(45)

The optimal model is identified by minimizing $L_T\left(\vartheta \right)$. [38] and [39] show that $L_T\left(\vartheta \right)$ effectively identifies the correct model across various non-Gaussian distributions. This spectral function is also useful for parameter estimation. Additionally, for $k=2$, the function $L_{2T}$ converges to the Whittle estimator.

2.2.3. Portmanteau-Type i.i.d. Test

We use the Portmanteau-type test proposed by [42]; however, we employ a simplified version that only considers the autocorrelation of the residuals and the squared residuals.3 Under these conditions, the test converges to the one proposed by [19], which combines the autocorrelation function (ACF) of the errors, ${\rho }_{\epsilon }\left(k\right)$, and the squared errors, ${\rho }_{{\epsilon }^2}\left(k\right)$, into a single function. The premise is that if ${\epsilon }_t$ is i.i.d., then ${\rho }_{\epsilon ,k}={\rho }_{{\epsilon }^2,k}=0$. The test includes both individual lag and cumulative versions:

$$J_{\epsilon ,{\epsilon }^2,k}={\displaystyle \frac{T^2}{T-k}}\left({\widehat{\rho }}_{\epsilon ,k}^2+{\widehat{\rho }}_{{\epsilon }^2,k}^2\right),C_{\epsilon ,{\epsilon }^2,m}=\sum _{k=1}^m{J}_{\epsilon ,{\epsilon }^2,k}.$$

(46)

The null hypothesis for the individual lag is $H_0:{\rho }_{\epsilon ,k}=0,{\rho }_{{\epsilon }^2,k}=0$ for $k\ge 1$, and for the cumulative test, $H_0:{\rho }_{\epsilon ,k}=0,{\rho }_{{\epsilon }^2,k}=0$ for $k=1,\dots ,m\ge 1$, where m is the number of lags. The asymptotic distributions of the tests are given by:

$$J_{\epsilon ,{\epsilon }^2,k}\stackrel{d}{\to }{\chi }_2^2,C_{\epsilon ,{\epsilon }^2,m}\stackrel{d}{\to }{\chi }_{2m}^2\mathrm{as}T\to \infty .$$

(47)

These tests effectively detect dependence, conditional heteroscedasticity, or nonstationarity. The ${\chi }_{2m}^2$ distribution is a suitable approximation for the asymptotic distribution of the test for non-Gaussian processes. The identification of the MARMA process is determined by the model for which the null hypothesis is not rejected. If the null hypothesis is rejected for two models, the model with the smallest value of the statistic $C_{\epsilon ,{\epsilon }^2,m}$ is selected.

The spectral function integrates third- and fourth-order spectral density functions, effectively identifying the correct MARMA model when the error distribution exhibits skewness or kurtosis. Conversely, the test by [19] utilizes information from the ACF of the residuals and the squared residuals in a cumulative function, enabling the identification of the correct model based on the minimum value of the function evaluated across all possible representations. The performance of both functions in identifying MARMA is proposed in [39]. Under non-Gaussianity, both methods are reliable for identification. However, the function by [50] yields better results for asymmetric distributions.

2.2.4. Whittle Estimation for the GARCH($p,q$)

Whittle estimation for the GARCH($p,q$) model was initially proposed by [27], using the ARMA(max($p,q$),q) representation as suggested in Equation 23. Although this relationship indicates that Whittle estimation can be used for parameter estimation of the GARCH model, the assumptions on the moments of $Y_t$ are more stringent; specifically, $\mathbb{E}\left[Y_t^8\right]<\infty$, which is quite restrictive. This is in contrast with empirical evidence on financial returns, which sometimes exhibit infinite second, third, or fourth moments. Later, [48] demonstrated that Whittle estimation for the squared GARCH process is feasible with only the existence of the fourth moment of $Y_t$. However, the convergence of the estimated parameters is significantly slower compared to the time-domain maximum likelihood estimation.

The Whittle estimator of ${\xi }_0=\left\{{\alpha }_{0,1},\dots ,{\alpha }_{0,\max \left(q\right)},{\beta }_{0,1},\dots ,{\beta }_{0,q}\right\}$ is:

$$W(\xi ,\omega )=\sum _{\ell =1}^{T/2}log\left(S_2^v(\xi ,{\kappa }_2,\omega )\right)+\sum _{\ell =1}^{T/2}{\displaystyle \frac{I_2\left(\omega \right)}{S_2^v(\xi ,{\kappa }_2,\omega )}}.$$

(48)

The estimated parameter set is obtained by minimizing the function:

$$\widehat{\xi }=\underset{\xi \in \Xi }{arg\; min}W(\xi ,{\kappa }_2,\omega ).$$

(49)

Unlike the maximum likelihood, the Whittle estimator relies solely on the values $Y_1^2,\dots ,Y_T^2$. This means that the unobservable values ${\sigma }_t^2$ do not need to be evaluated, and there is no need to choose initial values for ${\sigma }_t^2$. This characteristic can be seen as an advantage of the Whittle estimator. The parameter $\omega$ is estimated by the unconditional variance, $\mathrm{var}\left(Y_t\right)=\omega /(1-\alpha -\beta )$:

$$\widehat{\omega }=\mathrm{var}\left(Y_t\right)(1-\widehat{\alpha }-\widehat{\beta }).$$

(50)

The question arises regarding the potential use of Whittle estimation for the GARCH model in the context of noncausality or noninvertibility of the mean model. Although its implementation can be straightforward, there are gaps concerning its validity, given that the mean and volatility models are separated. While the mean can depend on future values, the volatility, in this case, only depends on past shocks.

2.2.5. Maximum Likelihood Estimation for the GARCH($p,q$)

The most common estimation method for the GARCH model is through the maximum likelihood function, as it is consistent and asymptotically normal for stationary GARCH processes [24]. In contrast to the Whittle estimator, which requires the existence of the fourth moments of $Y_t$ for GARCH models, the maximum likelihood function does not make assumptions about the moments. Given that we do not make assumptions about the distribution of the innovations ${\eta }_t$, but we use the likelihood function of the normal distribution, this estimator is usually referred to as quasi-maximum likelihood.

Given the parameter vector of the GARCH($p,q$): $\xi =({\alpha }_1,\dots ,{\alpha }_p,{\beta }_1,\dots ,{\beta }_q)\in {\mathbb{R}}^{p+q}$, the likelihood function is conditioned to the initial values ${\epsilon }_0,{\epsilon }_{1-q},{\tilde{\sigma }}_0^2,\dots ,{\tilde{\sigma }}_{1-p}$, and can be defined by:

$${\mathrm{Lik}}_T(\xi ,{\epsilon }_1,\dots ,{\epsilon }_T)=\prod _{t=1}^T{\displaystyle \frac{1}{\sqrt{2\pi {\tilde{\sigma }}_t^2}}}exp\left(-{\displaystyle \frac{{\epsilon }_t^2}{2{\tilde{\sigma }}_t^2}}\right),$$

(51)

where ${\tilde{\sigma }}_t^2$ is as in Equation 3, and recursively generated for $t\ge 1$, ${\tilde{\sigma }}_t^2\left(\xi \right)=\omega +{\sum }_{i=1}^q{\epsilon }_{t-i}^2+{\sum }_{j=1}^p{\epsilon }_{t-j}^2$. The log-likelihood function is defined as:

$$log\left({\mathrm{Lik}}_T\left(\xi \right)\right)=\sum _{t=1}^T\left(-{\displaystyle \frac{1}{2}}log\left(2\pi \right)-{\displaystyle \frac{1}{2}}log\left({\tilde{\sigma }}_t^2\right)-{\displaystyle \frac{{\epsilon }_t^2}{2{\tilde{\sigma }}_t^2}}\right).$$

(52)

The minimization obtains the estimation of the parameters:

$${\widehat{\xi }}_T=\underset{\xi \in \Xi }{arg\; min}log\left({\mathrm{Lik}}_T\left(\xi \right)\right).$$

(53)

The estimator is consistent, ${\widehat{\xi }}_T\stackrel{p}{\to }{\xi }_0$, given that ${\widehat{\xi }}_T$ converges in probability to the true parameter set ${\xi }_0$ as $T\to \infty$. The true parameter ${\xi }_0$ is identifiable if the likelihood function has a unique minimum at ${\xi }_0$. This identifiability condition ensures that no other parameter vector can produce the same likelihood as the true parameter vector. Furthermore, the estimation function converges uniformly to its expectation. Specifically, the log-likelihood function converges in probability to its expected value as T grows. This is expressed as ${\displaystyle \frac{1}{T}}log\left({\mathrm{Lik}}_T\left(\xi \right)\right)\stackrel{p}{\to }\mathbb{E}[log\left({\mathrm{Lik}}_T\left(\xi \right)\right)]$ for all $\xi \in \Xi$. The uniform convergence implies that for any small $ϵ>0$ and $\delta >0$, there exists a sample size $T_0$ such that for all $T>T_0$, the probability that the supremum of the absolute difference between the normalized log-likelihood function and its expectation is less than $ϵ$ exceeds $1-\delta$. Therefore, the convergence of the log-likelihood function to its expectation does not depend on the specific parameter value but holds uniformly over the entire parameter space $\Xi$.

The estimator ${\widehat{\xi }}_T$ is asymptotically normal as $T\to \infty$. This is expressed as $\sqrt{T}({\widehat{\xi }}_T-{\xi }_0)\stackrel{d}{\to }\mathcal{N}(0,I^{-1}\left({\xi }_0\right))$, where $I\left({\xi }_0\right)$ is the Fisher information matrix evaluated at the true parameter values. The Fisher information matrix is defined as:

$$I\left({\xi }_0\right)=\mathbb{E}\left[-{\displaystyle \frac{{\partial }^{2}logf\left(Y_t\right|{\mathcal{F}}_{t-1};{\xi }_0)}{\partial \xi \partial {\xi }^{\prime }}}\right],$$

where $f\left(Y_t\right|{\mathcal{F}}_{t-1};{\xi }_0)$ is the conditional density of $Y_t$ given past information ${\mathcal{F}}_{t-1}$. The score function, defined as $S_T\left(\xi \right)={\displaystyle \frac{\partial log\left({\mathrm{Lik}}_T\left(\xi \right)\right)}{\partial \xi }}$, satisfies the Central Limit Theorem; therefore, the score vector evaluated at the true parameter values has an asymptotically normal distribution. Additionally, the Hessian matrix of the log-likelihood function, normalized by T, must converge in probability to the Fisher information matrix. Formally, ${\displaystyle \frac{1}{T}}{\sum }_{t=1}^T{\displaystyle \frac{{\partial }^{2}log\left({\mathrm{Lik}}_T\left(\xi \right)\right)}{\partial \xi \partial {\xi }^{\prime }}}\stackrel{p}{\to }I\left({\xi }_0\right)$. Consequently, the distribution of the estimator is approximately normal, centered at the true parameter values, and with a covariance matrix given by the inverse of the Fisher information matrix.

The assumptions regarding the moments of the process are as follows. First, the variance of the $Y_t$ process must be finite, $\mathbb{E}\left[Y_t^2\right]<\infty$, to ensure the model parameters are estimable. Additionally, the existence of fourth-order moments, $\mathbb{E}\left[Y_t^4\right]<\infty$, is necessary to derive the asymptotic properties of the estimators, ensuring that the sample moments converge to their population counterparts. Specifically, the Fisher information matrix, $I\left(\xi \right)$, must be finite and positive definite to ensure that the parameter estimates have a well-defined asymptotic distribution and that the covariance matrix of the estimators is invertible.

These assumptions hold under the condition that $Y_t$ is strictly stationary. However, financial asset prices typically exhibit trends and are not stationary. Therefore, we propose techniques for detrending in the following sections.

2.3. Detrending Method

Removing the trend from data is a topic of considerable debate since the choice of method can directly influence the emergence of dynamics. For instance, the econometrician must make assumptions about the nature of the trend, whether it is stochastic or deterministic. To address this, the econometrician may utilize tests such as the augmented Dickey-Fuller or Phillips-Perron tests. However, the outcomes of these tests depend on the choice of lag length, and the empirical distribution of the test may not align with the data generating process, potentially leading to inference errors.

When selecting a stochastic trend, the most common method is to differentiate the data, though the optimal number of differences needs to be well-established. Another approach could be a logarithmic transformation. In both cases, the primary concern is that incorrect differencing or transformation can remove or add relevant information about the time series dynamics, leading to the emergence of spurious dynamics.

In the case of choosing a deterministic trend, the econometrician must determine the optimal number of polynomials to fit the trend. Practically, this is done through visual inspection and information criteria, although there is no clear rule for this choice, which can result in either the appearance of spurious dynamics or the abstraction of significant dynamics. Various methods have been proposed in the literature. For instance, [41] employed a cubic trend, while [25] and [1] utilized a linear trend or simple centering. [18] used returns, whereas [37] applied the X-11 method to filter seasonal effects. Another method is the Hodrick-Prescott filter, which has been used by [40] and shows that for monthly data, the filter may not induce dynamics. However, in practice, the selection of the lambda parameter influences the dynamics obtained in the time series; see [34] for details. Recently, a filtering approach was proposed by [5], based on a joint model that includes both a random-walk fundamental component and a bubble component specified as a MAR model. Both components are treated simultaneously in the filter updating process, allowing the test to be implemented in non-stationary settings.

[33] propose to use cubic splines, which preserve the characteristics of bubble formation and local explosions observed in Bitcoin and Ethereum. The detrended series is found to be stationary and may exhibit noncausal dynamics. Despite using a second-order spline with a knot every 30 observations to achieve a stationary series with possible noncausal dynamics, the correct order of splines for the given data remains to be determined, as these splines can also influence the existence of dynamics.

Lastly, a well-known technique for detrending data in the frequency domain involves setting the frequencies corresponding to trend information above a certain time threshold to zero (see [7]; [13]). For example, frequencies representing trends of one year, six months, or one month can be removed while preserving the medium and high frequencies that are associated with the dynamics and white noise component of the time series. The key question in this case is determining the appropriate level so that the detrending process neither adds nor removes significant dynamics from the time series. In [45], the trend component shorter than one year is removed, resulting in a stationary series for the trading volumes of various stocks.

In this paper, we employ two methods for detrending: cubic splines and the Fourier method for eliminating low frequencies. We aim to explore whether these methods can adjust the dynamics of MARMA models and allow for a GARCH process in the squared residuals. The analysis is empirical, with theoretical exploration reserved for future research. The choice of parameters in the cubic splines method and the choice of frequencies in the Fourier method directly influence the emergence of dynamics; the more the series is penalized, the fewer dynamics it exhibits. However, this detrending is done to apply our methodology, as determining the optimal type of detrending is not the primary focus of this paper.

2.3.1. Cubic Splines

Cubic splines are piecewise polynomial functions that provide a smooth fit to the data while maintaining continuity and differentiability. They are particularly advantageous in detrending because they offer a high degree of smoothness, which is necessary for accurately modeling and removing trends without distorting the underlying signal. For details on the nature of splines and their use, see [20,23,30,36]. Given the time series $Y_t$ for $t=1,\dots ,T-1$, a cubic spline $S_t\left(Y\right)$ is a function that satisfies: 1. On each interval $[Y_t,Y_{t+1}]$, the cubic spline function $S_t\left(Y\right)=a_t+b_t(Y-Y_t)+c_t{(Y-Y_t)}^2+d_t{(Y-Y_t)}^3$, for $Y_t\le Y\le Y_{t+1}$. Where Y represents any point within the interval $[Y_t,Y_{t+1}]$. 2. The function $S_t\left(Y\right)$, the first and second derivatives $S_t^{\prime }\left(Y\right)$, $S_t^{\prime \prime }\left(Y\right)$ are continuous. 3. The second derivatives at the beginning and end of the time series are zero, that is, $S^{\prime \prime }\left(Y_0\right)=0$, and $S^{\prime \prime }\left(Y_T\right)=0$.

We divide the time series into n knots to implement the cubic splines. In total, there are $T/n$ knots. In each interval, we obtain the parameters $a_t,b_t,c_t$, and $d_t$ by solving the tridiagonal system:

$${\displaystyle \frac{h_{t-1}}{6}}{M}_{t-1}+{\displaystyle \frac{h_{t-1}+h_t}{3}}{M}_t+{\displaystyle \frac{h_t}{6}}{M}_{t+1}={\displaystyle \frac{y_{t+1}-y_t}{h_t}}-{\displaystyle \frac{y_t-y_{t-1}}{h_{t-1}}}$$

where $M_t=S^{\prime \prime }\left(y_t\right)$, and $h_t=y_{t+1}-y_t$. Once the second derivatives $M_t$ are determined, follows the computation of the coefficients $a_t,b_t,c_t,d_t$:

$$a_t=y_t,b_t={\displaystyle \frac{y_{t+1}-y_t}{h_t}}-{\displaystyle \frac{h_t}{3}}(2M_t+M_{t+1}),c_t=M_t,d_t={\displaystyle \frac{M_{t+1}-M_t}{3h_t}}.$$

Once the parameters are determined, we calculate the estimated values at each knot. By aggregating these, we construct the trend of $Y_t$, defined as $t\left(Y\right)$. The detrended series is

$$Y_t^{\mathrm{detrended}}=Y_t-t\left(Y\right).$$

Fourier Detrending

For this method, we first perform $d_T\left(\omega \right)={\sum }_{t=1}^T{Y}_t{e}^{-it\omega }$ is the discrete Fourier transform. The transformed series $d_T\left(\omega \right)$ is defined for the frequencies ${\omega }_{\ell }$, ${\displaystyle \frac{2\pi \ell }{T}}$, $\ell =1,\dots ,T$. The low frequencies are defined in the set $\ell \in [1,D]$, where $D\in \mathbb{N}$ and $d=D/T$. [45] find that they can detrend stock trading volume by nullifying the frequencies related to cycles longer than one year, $D=10$.

The proportion d is the choice that the econometrician must make to detrend the time series. Low frequencies are associated with the long-term trend of the data, while high frequencies are related to short-term variations and the white noise process. The larger the value of d, the greater the penalization on these frequencies, which can result in the loss of relevant information about the dynamics. Typically, the choice of d is made through visual inspection, observing the prominent peaks in the periodogram.

After establishing d, we set $d_T\left(\omega \right)=0$ in the frequencies $\ell =1,\dots ,D$, and $\ell =T-D,\dots ,D$. In a subsequent step, the inverse Fourier transform converts the modified frequency domain series back into the time domain:

$$Y_t^{\mathrm{detrended}}=\mathrm{Re}\left({\displaystyle \frac{1}{T}}\sum _{j=0}^T{d}_T\left(\omega \right)e^{ij\omega }\right).$$

The new series has T observations and is trend-free since the frequencies associated with long-term components have been set to zero.

3. Results

3.1. Simulation Study

In this section, we compare Whittle and maximum likelihood estimation methods. Our focus is on the parameters of the GARCH model. Indeed, the behavior of the Whittle estimator has already been investigated for the MARMA model by [39].

To simulate the MARMA($r,s,r^{\prime },s^{\prime }$)-GARCH($p,q$) process, we follow the next six steps:

• Generate an error sequence: Simulate T realizations of an i.i.d. random error sequence: ${\left\{{\eta }_t\right\}}_{t=0}^T$, where ${\eta }_t\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathfrak{F}\left(\zeta \right)$. Here, $\mathfrak{F}$ is a chosen distribution function with parameters $\zeta$. We ensure $\eta$ has zero mean and unit variance.
• Generate the conditional variance:

  (a)

  Initialize ${\tilde{\sigma }}_1^2$ and compute ${\sigma }_t^2=\omega +(\alpha {\eta }_{t-1}+\beta ){\sigma }_{t-1}^2$.

  (b)

  Determine ${\epsilon }_t={\eta }_t{\sigma }_t$ for $t=2,\cdots ,T$.

• Compute the discrete Fourier transform (DFT): Perform the DFT of ${\epsilon }_t$: $d_{\epsilon }\left(\omega \right)={\sum }_{t=1}^T{\epsilon }_t{e}^{-it\omega }$.
• Select MARMA parameters: Choose parameters and orders for the MARMA($r,s,r^{\prime },s^{\prime }$) process and define the transfer function:

  $$\psi (\vartheta ,\omega )={\displaystyle \frac{(1+{\theta }_1^{+}z+\dots +{\theta }_{r^{\prime }}^{+}{z}^{r^{\prime }})(1+{\theta }_1^{*}{z}^{-1}+\dots +{\theta }_{s^{\prime }}^{*}{z}^{-s^{\prime }})}{(1-{\varphi }_1^{+}z-\dots -{\varphi }_r^{+}{z}^r)(1-{\varphi }_1^{*}{z}^{-1}-\dots -{\varphi }_s^{*}{z}^{-s})}}.$$
• Compute DFT of $Y_t$:$d_T\left(\omega \right)=d_{\epsilon }\left(\omega \right)\psi (\vartheta ,\omega )$.
• Recover $Y_t$: Obtain $Y_t$ as the real part of the inverse Fourier transform of $d_T\left(\omega \right)$:

  $$Y_t=Re\left({\displaystyle \frac{1}{T}}\sum _{j=0}^T{d}_T\left(\omega \right)e^{ij\omega }\right).$$

Figure 1 displays the trajectory of the MARMA(1,1,0,0)-GARCH(1,1) model along with the conditional variance. The top part of the figure shows the conditional variance of the GARCH(1,1) model, exhibiting volatility clustering and high sensitivity to innovations. The bottom part displays the time series $Y_t$, which features asymmetric cycles, local explosions, and patterns reminiscent of financial bubbles. These simulated patterns mirror those seen in real financial markets, highlighting the ability of the model to capture non-linear behaviors and extreme events.

Table 1 reports the performance of the Whittle and the maximum likelihood estimation methods for the parameters of a weak white noise GARCH(1,1) model ($\omega$, $\alpha$, $\beta$). Several density specifications as well as three large sample sizes are considered. We observe that maximum likelihood estimation consistently outperforms the Whittle approach in terms of both bias and RMSE across different distributions. The Whittle estimator is particularly sensitive to the heavy tails of the GARCH process, leading to higher bias and RMSE values. This sensitivity results in slower convergence of parameters to their theoretical values compared to the maximum likelihood method. The Whittle estimator has larger bias and RMSE values for all parameters and distributions, but it is very pronounced under heavy-tailed ones, such as Student’s t (4.5) and skew-t (4.5, 1.5). This method also gives significantly high sensitivity to the properties of the distribution, providing poorer parameter estimates in practice. The maximum likelihood estimation, on the other hand, shows smaller bias and lower RMSE values, giving higher levels of precision and accuracy. The results also reveal that the two estimation methods do get better with an increase in the sample size. However, the maximum likelihood estimation is superior even with smaller samples. These results show that maximum likelihood estimation is very robust and is generally the preferred choice when specifying models for GARCH(1,1) with its financial applications.

In Table 2, we simulate and estimate the MARMA(1,1,0,0)-GARCH(1,1) process using maximum likelihood. The purpose of this exercise is to compare the parameter estimation using the causal AR(2)-GARCH(1,1) representation with the correct model from which the data originate, the MARMA(1,1,0,0)-GARCH(1,1). In the case of conditional mean fitting, we can assert that both error sequences are uncorrelated. However, the causal and invertible representation is not independent; it indeed depends on past values of ($Y_t$), interactions of the conditional variance with lags and leads, and on the third and fourth cumulants of ($Y_t$) [39]. This complexity is reflected in the estimation of conditional variance, which does not correctly capture the higher-order dependency structures. This translates into biased parameters and significant patterns in the autocorrelation function of the squared residuals.

Clearly, as seen in Figure 2, the estimation of GARCH model parameters on the residuals of the causal AR(2) representation of the MARMA(1,1) model leads to biased parameter estimates. The parameter ($\omega$) tends to be higher than its true value, indicating an erroneous long-term volatility level. The parameter ($\alpha$) also tends to be higher than its true value, incorrectly suggesting a greater impact of innovations on conditional volatility. Lastly, the parameter ($\beta$) tends to be lower than its true value, implying a lower level of conditional volatility clustering. On the other hand, the estimation of the MARMA(1,1)-GARCH(1,1) model yields unbiased parameters, which decrease in variance as the sample size increases. In summary, the choice of estimation method for GARCH model parameters has significant practical implications. The maximum likelihood estimation yields more accurate and reliable parameter estimates compared to Whittle estimation, particularly in the presence of heavy-tailed distributions. This accuracy is crucial for applications in financial risk management, where precise modeling of volatility and tail risk is essential.

3.2. Empirical Applications in Cryptocurrencies

We consider the two main cryptocurrencies, the Bitcoin and the Ethereum, observed from January 2021 to July 2024. Cryptocurrencies are ideal to implement the MARMA-GARCH model as they exhibit extreme volatility [15,17,32] for which we will be using GARCH models as well as speculative bubbles on which we estimate MARMA models ([33]). For these two cryptocurrencies, the second-order correlation structure is very weak. However, they have strong high-order dynamics, including extreme and persistent volatility, as well as volatility clustering ([2]). The notable existence of time-varying volatility has motivated the use of a large number of models from the ARCH family, resulting in independent residuals, among others see [3,9,10,11,21,28,31,43]

In Figure 3, the prices of Bitcoin and Ethereum are displayed. In 2021, the prices of Bitcoin and Ethereum reached approximately 60,000 and 4,000 dollars respectively. This was mainly driven by the interest of institutional investors who massively bought the cryptocurrencies, including them in their investment portfolios. Despite this growing demand, the prices of both cryptocurrencies plummeted dramatically when scrutiny over the legality of their possession, mining, and trading intensified, with a direct ban in China, reaching a value in July 2021 of 30,000 dollars for Bitcoin and 1,800 for Ethereum. This episode consisted of a speculative bubble, which has been empirically confirmed in [4]. In 2022, there was a massive recovery in prices driven by uncertainty regarding the economic consequences of Russia’s invasion of Ukraine. Despite this, the market anticipated an unrestrained increase in inflation, driven by rising energy prices. Consequently, investors continued the trend of holding less risky assets. This trend persisted until the second half of 2023. During this period, major central banks indicated the possibility of interest rate cuts, reactivating investors’ risk appetite. Additionally, the creation of standardized products such as cryptocurrency futures generated a positive sentiment of confidence, resulting in a widespread increase in prices.

The prices of Bitcoin and Ethereum are nonstationarity. They exhibit frequent episodes of local explosions that tend to collapse abruptly. This behavior resembles the well-known formation of financial bubbles. In addition, conditional volatility and clustering are evident. The prices also display non-Gaussian characteristics, including heavy tails resulting from extreme events and positive skewness. These features complicate methods for detrending data, making it more challenging to accurately extract trends without altering the short-term variations that capture the underlying dynamics. Following [33], we opted to use a cubic spline every 30 observations. We have also tried various interval sizes in days and a quadratic spline; however, the cubic spline with a 30-day interval looks more robust for our investigation. For the Fourier detrended approach, we also tested various values for d, and after visual inspection and using the same tests, we settled on $d=0.025$, which corresponds to the first 32 frequencies. In other studies, [45] exclude cycles longer than a year. Both methodologies are subject to discussion, as the choice of trend can induce contradictory dynamics in the time series. Despite this, we leave this discussion for future research, as it is not the objective of this paper.

Figure 4 and Figure 5 show the detrended prices of Bitcoin and Ethereum, respectively. Both methodologies correctly extract the price trends; however, extracting the trend becomes more challenging during episodes of local explosions, such as those in 2021 and 2022. The detrended prices exhibit signs of non-Gaussianity, extreme events, heavy tails, and skewness. In this context, conditional volatility is evident. The high volatility scenarios in 2021 and 2022 were driven by extreme volatility. After these periods, volatility decreased and was followed by moments of lower volatility. From the second half of 2023, volatility increased again, though it remained lower than the levels experienced in 2021. The descriptive statistics of the prices and the detrended prices are reported in Table 3. Both detrending methods yield similar series; however, the Fourier approach results in series with higher kurtosis. In various attempts with different values for d, we concluded that as more frequencies are penalized, extreme values decrease and kurtosis tends to approach normality. The same trade-off applies to the cubic splines.

Table 4 reports the estimation results. Using both AIC and BIC and for both cubic splines and Fourier detrending approaches, the model that best fits the conditional mean is the AR(1). Once this model is established, we begin exploring its possible noncausal representations. In this case, we invert the root of the polynomial $(1-{\varphi }^{+}z)$, resulting in $(1-{\varphi }^{*}{z}^{-1})$. We proceed to calculate the spectral function $L_T\left(\vartheta \right)$ and the Portmanteau test $C_{\epsilon ,{\epsilon }^2,m}$. For both series and both detrending methods, the identified model was the purely noncausal MARMA(0,1,0,0) with an estimated parameter around 0.825. The frequency domain residuals are obtained with both specifications using the following formula:

$$d_{\widehat{\epsilon }}\left(\omega \right)=d_T\left(\omega \right)(1-{\widehat{\varphi }}_1^{+}z-...-{\widehat{\varphi }}_r^{+}{z}^r)(1-{\widehat{\varphi }}_1^{*}{z}^{-1}-...-{\widehat{\varphi }}_s^{*}{z}^{-s}).$$

(54)

The time domain residuals ${\widehat{\epsilon }}_t$ are obtained by the real part of the inverseFourier transform of the frequency domain residuals:

$${\widehat{\epsilon }}_t=Re\left({\displaystyle \frac{1}{T}}\sum _{j=0}^T{d}_{\widehat{\epsilon }}\left(\omega \right)e^{ij\omega }\right).$$

(55)

In a subsequent step, we estimate the parameters of the GARCH(1,1) model. We observe that the parameter ($\omega$), which represents the long-term variance, is always higher for the MARMA(0,1,0,0)-GARCH(1,1) model than for its causal AR(1)-GARCH(1,1) representation. The parameter ($\alpha$), which dictates the response of conditional volatility to innovations, is consistently greater for the MARMA(0,1,0,0)-GARCH(1,1) model than for its causal AR(1)-GARCH(1,1) counterpart. The greater impact of innovations reveals a stronger reaction of conditional volatility to new market information. Unlike the traditional ARMA-GARCH model, the MARMA-GARCH model allows this information to come from future price expectations. Conversely, the parameter ($\beta$) is consistently lower in all cases, indicating that the squared residual dependence structure of the MARMA(0,1,0,0)-GARCH(1,1) model has a lower tendency to form volatility clusters. We report the p-value of the Ljung-Box test for the residuals and the McLeod-Li test for the standardized residuals. In no case are there significant lags, indicating the absence of significant dynamics and providing evidence that the residuals are i.i.d.

The magnitude of the estimated parameters reveals that the causal AR(1)-GARCH(1,1) model underestimates the impact of news or new information entering the market. Consequently, this approach underestimates the potential occurrence of extreme values arising from abnormal events in the financial market. This, in turn, results in a lower impact on conditional volatility and, therefore, an underestimation of the potential risk associated with the cryptocurrency. In contrast, the MARMA(0,1,0,0)-GARCH(1,1) model more realistically captures the impact of such abnormal events on conditional volatility. In Figure 6, we present the conditional volatility of Bitcoin and Ethereum. It is evident that the red line reacts more sharply than the black line during periods of high volatility, such as in 2021 and 2022. This difference is more pronounced for Bitcoin than for Ethereum. Conversely, during periods of low volatility, the conditional volatility represented by the MARMA(0,1,0,0)-GARCH(1,1) model is lower.

The identification of causal, noncausal, invertible, or noninvertible dynamics, as well as the magnitude of the estimated coefficients, directly depends on the method used to remove the trend from the data. In this process, assumptions are made about the nature of the trend and the adjustment of the model parameters used to remove it. This can result in the initial properties of the data being removed or fictitious dynamics being introduced. In the two methods we proposed, cubic splines and Fourier detrending, we found that the properties of the detrended series are similar. However, this is not always the case.

4. Discussion and Conclusions

In this paper, we propose the new MARMA-GARCH model, which can appropriately accommodate the explosive episodes observed during financial bubbles and conditional volatility. We propose two methods for its estimation: the Whittle estimation function in the frequency domain and Gaussian likelihood in the time domain. For the estimation of the MARMA component, both methods yield similar results, which only require a finite variance environment. Conversely, for the estimation of the GARCH model, maximum likelihood estimation is superior, as it is more flexible than Whittle estimation in the presence of high kurtosis and skewness. This result is confirmed in a Monte Carlo simulation. Furthermore, the incorrect representation of the mean model, such as omitting noncausality or noninvertibility of the series and instead erroneously imposing roots outside the unit disk, leads to biased parameter estimation in the GARCH model. This misrepresentation fails to capture the effect of news correctly and underestimates the clustering effect. We apply this model to cryptocurrencies, fitting a MARMA(0,1,0,0) model to the mean. This model is simply a noncausal AR(1) model. The results indicated that incorrectly imposing a causal AR(1) model leads to erroneous results in the estimation of the GARCH(1,1) model, specifically underestimating the parameter $\alpha$ and overestimating the parameter $\beta$. Such biases in the parameters necessarily result in poorer data fitting and cast doubt on the reliability of statistics derived from the model, such as those used in financial quantitative risk assessment, like Value at Risk (VaR).