On the ARCH model with stationary liquidity

The classical ARCH model together with its extensions have been widely applied in the modeling of financial time series. We study a variant of the ARCH model that takes account of liquidity given by a positive stationary process. We provide minimal assumptions that ensure the existence and uniqueness of the stationary solution for this model. Moreover, we give necessary and sufficient conditions for the existence of the autocovariance function. After that, we derive an AR(1) characterization for the stationary solution yielding Yule–Walker type quadratic equations for the model parameters. In order to define a proper estimation method for the model, we first show that the autocovariance estimators of the stationary solution are consistent under relatively mild assumptions. Consequently, we prove that the natural estimators arising out of the quadratic equations inherit consistency from the autocovariance estimators. Finally, we illustrate our results with several examples and a simulation study.


Introduction
The ARCH and GARCH models have become important tools in time series analysis. The ARCH model was introduced in Engle (1982) and then it has been generalized to the GARCH model by Bollerslev (1986). Since, a large collection of variants and extensions of these models has been produced by many authors. See for example Bollerslev (2008) for a glossary of models derived from ARCH and GARCH. On a related work we also mention (Han 2013), where GARCH-X model with liquidity arising from a certain fractional ARMA process is considered.
In this work, we also focus on a generalization of the ARCH model, namely the model (1). Our contribution proposes to include in the expression of the squared volatility σ 2 t a factor L t−1 , which we will call liquidity. The motivation to consider such a model comes from mathematical finance, where the factor L t , which constitutes a proxi for the trading volume at day t, has been included in order to capture the fluctuations of the intra-day price in financial markets. A more detailed explanation can be found in Bahamonde et al. (2018) or Tudor and Tudor (2014). In the work (Bahamonde et al. 2018) we considered the particular case when L t is the squared increment of the fractional Brownian motion (fBm in the sequel), i.e. L t = (B H t+1 − B H t ) 2 , where B H is a fBm with Hurst parameter H ∈ (0, 1). In this work, our purpose is twofold. Firstly, we enlarge the ARCH with fBm liquidity in Bahamonde et al. (2018) by considering, as a proxi for the liquidity, a general positive (strictly stationary) process (L t ) t∈Z . This includes, besides the above mentioned case of the squared increment of the fBm, many other examples.
The second purpose is to provide a method to estimate the parameters of the model. As mentioned in Bahamonde et al. (2018), in the case when L is a process without independent increments, the usual approaches for the parameter estimation in ARCH models (such as least squares method and maximum likelihood method) do not work, in the sense that the estimators obtained by these classical methods are biased and not consistent. Here we adopt a different technique, based on the AR(1) characterization of the ARCH process, which has also been used in Voutilainen et al. (2017). The AR(1) characterization leads to Yule-Walker type equations for the parameters of the model. These equations are of quadratic form and then we are able to find explicit formulas for the estimators. We prove that the estimators are consistent by using extended version of the law of large numbers and by assuming enough regularity for the correlation structure of the liquidity process. We also provide a numerical analysis of the estimators.
The rest of the paper is organised as follows. In Sect. 2 we introduce our model and discuss the existence and uniqueness of the solution. We also provide necessary and sufficient conditions for the existence of the autocovariance function. We derive the AR(1) characterization and Yule-Walker type equations for the parameters of the model. Section 3 is devoted to the estimation of the model parameters. We construct estimators in a closed form and we prove their consistency via extended versions of the law of large numbers and a control of the behaviour of the covariance of the liquidity process. Several examples are discussed in details. In particular, we study squared increments of the fBm, squared increments of the compensated Poisson process, and the squared increments of the Rosenblatt process. We end the paper with a numerical analysis of our estimators. All the proofs and auxiliary lemmas are postponed to the appendix.

The model
Our variant of the ARCH model is defined for every t ∈ Z as where α 0 ≥ 0, α 1 , l 1 > 0, and ( t ) t∈Z is an i.i.d. process with E( 0 ) = 0 and E( 2 0 ) = 1. Moreover, we assume that (L t ) t∈Z is a positive process and independent of ( t ) t∈Z .
Remark 1 By setting L t−1 = σ 2 t−1 in (1) we would obtain the defining equations of the classical GARCH(1, 1) process. However, in this case, the processes L and would not be independent of each other.
In Sect. 3, in the estimation of the model parameters, we further assume that (L t ) t∈Z is strictly stationary with E(L 0 ) = 1. However, we first give sufficient conditions to ensure the existence of a solution in a general setting where L is only assumed to be bounded in L 2 . This allows one to introduce ARCH models that are not based on stationarity. Note that we have a recursion Let us denote A t = α 1 2 t and B t = α 0 + l 1 L t for every t ∈ Z.

Existence and uniqueness of the solution
In the case of strictly stationary L, the existence and uniqueness of the solution is studied in Brandt (1986) and Karlsen (1990). However, in order to allow flexibility and non-stationarity in our class of ARCH models, we present a general existence and uniqueness result. Furthermore, our result allows one to consider only weakly stationary sequences. In addition, our proof is based on a different technique, adapted from Bahamonde et al. (2018). We start with the following theorem providing the existence and uniqueness of a solution under relatively weak assumptions (we only assume the boundedness of L in L 2 and the usual condition α 1 < 1 [see e.g. Francq and Zakoian (2010)].
Theorem 1 Assume that sup t∈Z E(L 2 t ) < ∞ and α 1 < 1. Then (1) has the following solution Moreover, the solution is unique in the class of processes satisfying sup t∈Z E(σ 2 t ) < ∞. The following result provides us a strictly stationary solution provided that L is strictly stationary. While the result is a special case of Karlsen (1990), for the reader's convenience we present a different proof that can be applied to the case of weak stationarity as well (see Corollary 2).
In the sequel, we consider a strictly or weakly stationary liquidity (L t ) t∈Z and the corresponding unique solution given by (4). Therefore, we will implicitly assume that E(L 2 0 ) < ∞ and α 1 < 1.
In order to study covariance function or weak stationarity of the solution (4), we require that the moments E(σ 4 t ) exist. Necessary and sufficient conditions for this are given in the following lemma.
. Remark 2 As expected, in order to have finite moments of higher order we needed to pose a more restrictive assumption α 1 < 1 in the case of Gaussian innovations we obtain the well-known condition α 1 < 1 √ 3 [see e.g. Francq and Zakoian (2010) or Lindner (2009)]. An explicit expression of the fourth moment can be obtained when L is the squared increment of fBm [see Lemma 4 in Bahamonde et al. (2018)].
We end this section with the following result similar to Corollary 1 on the existence of weakly stationary solutions.
and L is weakly stationary. Then the unique solution (4) is weakly stationary.

Computation of the model parameters
In this section we consider the stationary solution and compute the parameters α 0 , α 1 , l 1 in (1) by using the autocovariance functions of X 2 and L. To this end, we use an AR(1) characterization of the ARCH process. From this characterization, we derive, using an idea from Voutilainen et al. (2017), a Yule-Walker equation of quadratic form for the parameters, that we can solve explicitly. This constitutes the basis of the construction of the estimators in the next section. From (1) it follows that if (σ 2 t ) t∈Z is stationary, then so is (X 2 t ) t∈Z . In addition Let us define an auxiliary process (Y t ) t∈Z by Now Y is a zero-mean stationary process satisfying By denoting we may write corresponding to the AR(1) characterization (Voutilainen et al. 2017) of Y t for 0 < α 1 < 1.
Assuming a(γ γ γ ) = 0 we obtain the following solutions for the model parameters α 1 and l 1 : and Again, α 0 is given by Remark 3 Note that here we assumed s(n 2 ) = 0 and a(γ γ γ ) = 0 which means that we choose n 1 , n 2 in a suitable way. Notice however, that these assumptions are not a restriction. Firstly, the case where s(n 2 ) = 0 for all n 2 = 0 corresponds to the more simple case where L is a sequence of uncorrelated random variables. Secondly, if s(n 2 ) = 0 and a(γ γ γ ) = 0, the second order term vanishes and we get a linear equation for α 1 . For detailed discussion on this phenomena, we refer to Voutilainen et al. (2017).
Remark 4 At first glimpse Eqs. (8) and (13) may seem useless as one needs to choose between signs. However, it usually suffices to know additional values of the covariance of the noise [see Voutilainen et al. (2017)]. In particular, it suffices that s(n) → 0 [see Voutilainen et al. (2019)].

Parameter estimation
In this section we discuss how to estimate the model parameters consistently from the observations under some knowledge of the covariance of the liquidity L. At this point we simply assume that the covariance structure of L is completely known. However, this is not necessary, as discussed in Remarks 5 and 6. As mentioned in the introduction, classical methods like maximum likelihood (MLE) or least squares method (LSE) may fail in the presence of memory. Indeed, while MLE is in many cases preferable, it requires the knowledge of the distributions so that the likelihood function can be computed. Compared to our method, we only require certain kind of asymptotic independence for the process L in terms of third and fourth order covariances (see Lemma 3). Unlike MLE, the LSE estimator does not require distributions to be known. However, in our model it may fail to be consistent. Indeed, this happens already in the case of squared increments of the fractional Brownian motion (Bahamonde et al. 2018). Based on formulas for the parameters provided in Sect. 2.2, it suffices that the covariances of X 2 can be estimated consistently. For simplicity, we assume that the liquidity (L t ) t∈Z is a strictly stationary sequence. The main reason why we prefer to keep the assumption of strict stationarity is that it simplifies the third and fourth order assumptions of Lemma 3 and also because our main examples of liquidities are strictly stationary processes (see Sect. 3.3). Nevertheless, the results either hold directly or can be modified with only a little effort to cover the case of weakly stationary sequences. We leave the details to the reader.

Consistency of autocovariance estimators
Assume that (X 2 1 , X 2 2 , . . . , X 2 N ) is an observed series from (X t ) t∈Z that is given by our model (1). We use the following estimator of the autocovariance function of X 2 whereX 2 is the sample mean of the observations. We show that the estimator above is consistent in two steps. Namely, we consider the sample mean and the term separately. If the both terms are consistent, consistency of the autocovariance estimator follows. Furthermore, if this holds for the lags involved in Theorems 2 and 3, also the corresponding model parameter estimators are consistent.
In addition, assume that for every fixed n, n 1 and n 2 it holds that Cov t+n converges in probability to E(X 2 0 X 2 n ) for every n ∈ Z. The next lemma provides us the missing piece of consistency of the covariance estimators. It can be proven similarly as Lemma 3, but as it is less involved, we leave the proof to the reader.
, then the sample meanμ converges in probability to E(X 2 0 ).

Remark 5
As one would expect, the assumptions of Lemma 4 are implied by the assumptions of Lemma 3, which on the other hand are only sufficient for our purpose. Indeed, it suffices that Y t = X 2 t X 2 t+n is mean-ergodic for the relevant lags n (see Theorems 2 and 3 , and Remark 6). This happens when the dependence within the stationary process Y vanishes sufficiently fast. The condition α 1 < E( 8 0 ) − 1 4 ensures the existence of an autocovariance function of Y . Furthermore, the assumptions made related to the asymptotic behavior of the covariances of L guarantee that the dependence structure of Y (measured by the autocovariance function) is such that the desired consistency of the autocovariance estimators follows. Finally, the assumptions related to the liquidity are very natural. Indeed, we only assume that the (linear) dependencies within the process L t vanish over time. Examples of L satisfying the required assumptions can be found in Sect. 3.3.
Remark 6 In addition to the mean-ergodicity discussed in Remark 5, it suffices that the autocovariance function s(·) of the liquidity L is known for the chosen lags 0 and n. Furthermore, if we can observe L, which is often the case, these quantities can also be estimated.
The proof of the the next theorem is basically the same as the proof of Theorem 2.
Remark 7 -Statements of Theorems 2 and 3 hold true also when g 0 (γ γ γ 0 ) = 0 and g(γ γ γ ) = 0, but in these cases the estimators do not necessarily become real valued as the sample size grows. -The estimators from Definitions 1 and 2 are of course related. In practice (see the next section) we use those from Definition 1 while those from Definition 2 are needed just in case when we need more information in order to choose the correct sign forα 1 , see Remark 4.
-Note that here we implicitly assumed that the correct sign can be chosen inα 1 . However, this is not a restriction as discussed.
The approach of Voutilainen et al. (2017) was motivated by the classical Yule-Walker equations of an AR( p) process, for which the corresponding estimators have the property that they yield a causal AR( p) process agreeing with the underlying assumption of the equations [see e.g. Brockwell and Davis (2013)]. In comparison, with finite samples, the above introduced estimators may produce invalid values, such as complex numbers or negative reals. Moreover, for α 1 we may obtain an estimatê α 1 ≥ E( 8 0 ) − 1 4 violating assumptions of Lemma 3. However, we would like to emphasize that, as discussed before, the assumptions of the lemma are not necessary for a consistent estimation procedure. It may also happen thatα 1 ≥ 1 and in this case, Theory 1 does not guarantee the existence of the unique solution for (1) together with its stationarity. However, even in this case, there might still exist a (unique) stationary solution (cf. AR(1) with |φ| > 1). A further analyze of properties of (1) could be a potential topic for future research. The above described unwanted estimates are of course more prone to occur with small samples, although the probability of producing such values depends also on the different components of the model (1), such as the true values of the parameters. In practice, the issue can be avoided e.g. by using indicator functions as in Voutilainen et al. (2017) forcing the estimators to the demanded intervals. We also refer to our simulation study in Sect. 4 and Appendix B, which show that with the largest sample size we always obtained valid estimates.

Examples
We will present several examples of stationary processes for which our main result stated in Theorem 2 apply. Our examples are constructed as where (Y t ) t∈R is a stochastic process with stationary increments. We discuss below the case when Y is a continuous Gaussian process (the fractional Brownian motion), a continuous non-Gaussian process (the Rosenblatt process), or a jump process (the compensated Poisson process).

The fractional Brownian motion
t∈R is a two-sided fractional Brownian motion with Hurst parameter H ∈ (0, 1). Recall that B H is a centered Gaussian process with covariance Let us verify that the conditions from Lemma 3 and Theorem 2 are satisfied by L t = (B H t+1 − B H t ) 2 . First, notice that [see Lemma 2 in Bahamonde et al. (2018)] that for since r H (t) behaves as t 2H −2 for t large.
Let us now turn to the third-order condition, i.e. Cov(L 0 L n , L t ) = E(L 0 L n L t ) − E(L 0 L n ) → 0 as t → ∞. We can suppose n ≥ 1 is fixed and t > n. For any three centered Gaussian random variables Z 1 , Z 2 , Z 3 with unit variance we have E(Z 2 1 Z 2 2 ) = 1 + 2(E(Z 1 Z 2 )) 2 and

By applying this formula to
where r H is given by (22). By (22), the above expression converges to zero as t → ∞.
Similarly for the fourth-order condition, the formulas are more complex but we can verify by standard calculations that, for every n 1 , n 2 ≥ 1 and for every t > max(n 1 , n 2 ), the quantity can be expressed as a polynomial (without term of degree zero) in r H (t), r H (t − n 1 ), r H (t + n 2 ), r H (t + n 2 − n 1 ) with coefficients depending on n 1 , n 2 . The conclusion is obtained by (22).

The compensated Poisson process
Let (N t ) t∈R be a Poisson process with intensity λ = 1. Recall that N is a cadlag adapted stochastic process, with independent increments, such that for every s < t, the random variable N t − N s follows a Poisson distribution with parameter t − s. Define the compensated Poisson process (Ñ t ) t∈R byÑ t = N t − t for every t ∈ R and let L t = (Ñ t+1 −Ñ t ) 2 . Clearly EL t = 1 for every t and, by the independence of the increments ofÑ , we have that for t large enough Cov(L 0 , L t ) = Cov(L 0 L n , L t ) = Cov(L 0 L n 1 , L t L t+n 2 ) = 0, so the conditions in Theorem 2 are fulfilled.

The Rosenblatt process
The (one-sided) Rosenblatt process (Z H t ) t≥0 is a self-similar stochastic process with stationary increments and long memory in the second Wiener chaos, i.e. it can be expressed as a multiple stochastic integral of order two with respect to the Wiener process. The Hurst parameter H belongs to ( 1 2 , 1) and it characterizes the main properties of the process. Its representation is where (W (y)) y∈R is Wiener process and f H is deterministic function such that R R f H (y 1 , y 2 ) 2 dy 1 dy 2 < ∞. See e.g. Tudor (2013) for a more complete exposition on the Rosenblatt process. The two-sided Rosenblatt process has been introduced in Coupek (2018). In particular, it has the same covariance as the fractional Brownian motion, so E(L t ) = E(Z H t+1 − Z H t ) 2 = 1 for every t. The use of the Rosenblatt process can be motivated by the presence of the long-memory in the emprical data for liquidity in financial markets, see Tsuji (2002).
The computation of the quantities Cov(L 0 , L t ), Cov(L 0 L n , L t ) and Cov(L 0 L n 1 , L t L t+n 2 ) requires rather technical tools from stochastic analysis including properties of multiple integrals and product formula which we prefer to avoid here. We only mention that the term Cov(L 0 , L t ) can be written as P(r H (t), r H ,1 (t)) where P is a polynomial without term of degree zero, r H is given by (22), while Since |u 1 − u 2 | H −1 |u 2 − u 3 + t| H −1 |u 3 − u 4 | H −1 |u 4 − u 1 + t| H −1 converges to zero as t → ∞ for every u i and since this integrand is bounded for t large by |u 1 − u 2 | H −1 |u 2 − u 3 | H −1 |u 3 − u 4 | H −1 |u 4 − u 1 | H −1 , which is integrable over [0, 1] 4 , we obtain, via the dominated convergence theorem, that Cov(L 0 , L t ) → t→∞ 0. Similarly, the quantities Cov(L 0 L n , L t ) and Cov(L 0 L n 1 , L t L t+n 2 ) can be also expressed as polynomials (without constant terms) of r H , r H ,k , k = 1, 2, 3, 4 where    where at least one set A i is (t, t + 1). Thus we may apply a similar argument as above.

Simulations
This section provides a visual illustration of convergence of the estimators (16), (17) and (18) when the liquidity process L is given by L t = (B H t+1 − B H t ) 2 with H = 4 5 . The simulation setting is the following. The i.i.d. process ( t ) t∈Z is assumed to be a sequence of standard normals. In this case the restriction given by Lemma 3 reads α 1 < 1 105 1 4 ≈ 0.31. The used lag is n = 1 and the true values of the model parameters are α 0 = 1, α 1 = 0.1 and l 1 = 0.5. The used sample sizes are N = 100, N = 1000 and N = 10000. The initial X 2 0 is set equal to 1.7. After the processes L t with t = 0, 1, . . . N − 2 and t with t = 1, 2, . . . N − 1 are simulated, the initial is used to generate σ 2 1 using (1). Together with 1 this gives X 2 1 , after which (1) yields the sample {X 2 0 , X 2 1 , . . . , X 2 N −1 }. We have simulated 1000 scenarios with each sample size and the corresponding histograms of the model parameter estimates are provided in Figs. 1, 2 and 3. Our simulations show that the behaviour of the limit distributions are close to Gaussian ones, as N increases. We also note that, since the estimators involve square roots, they may produce complex valued estimates. However, asymptotically the estimates become real. In the simulations, the sample sizes N = 100 and N = 1000 resulted complex valued estimates in 47.9% and 4.3% of the iterations respectively, whereas with the largest sample size all the estimates were real. For the histograms the complex valued estimates have been simply removed. Some illustrative tables are given in Appendix B.
It is straightforward to repeat the simulations with other Hurst indexes, or with completely different liquidities such as squared increments of the compensated Poisson process. In these cases, we obtain similar results.
The simulations have been carried out by using the version 1.1.456 of RStudio software on Ubuntu 16.04 LTS operating system. Fractional Brownian motion was simulated by using circFBM function from the package dvfBm.       Acknowledgements Open access funding provided by Aalto University. The authors wish to thank two anonymous referees for their insightful comments that helped to improve the manuscript. Pauliina Ilmonen and Lauri Viitasaari wish to thank Vilho, Yrjö and Kalle Väisälä Fund, Marko Voutilainen wishes to thank Magnus Ehrnrooth Foundation, and Soledad Torres wishes to thank National Fund for Scientific and Technological Development, Fondecyt 1171335, for support.

Compliance with ethical standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Auxiliary lemmas and proofs
The following lemma ensures that we are able to continue the recursion in (3) infinitely many times.
and sup t∈Z E(σ 4 t ) ≤ M 2 < ∞, then the convergence holds also almost surely.
Proof By independence of , we have proving the first part of the claim. For the second part, Chebysev's inequality implies which is summable by assumptions. Borel-Cantelli then implies almost surely proving the claim.

Proof of Theorem 1
We begin by showing that (4) is well-defined. That is, we prove that defines an almost surely finite random variable. First we observe that the summands above are non-negative and hence, the pathwise limits exist in and denote a n = By the root test it suffices to prove that −→ e E log A 0 = α 1 e E log 2 0 by the law of large numbers and continuous mapping theorem. By Jensen's inequality we obtain that where we have used Consider now the function f x (a):=x a for x ≥ 1 and a ≥ 0. Since f x (a) = x a log x we obtain by the mean value theorem that Hence On the other hand, for n ≥ 2 and L t−n ≥ 1 it holds that since for x ≥ 1, the function g(x):= (log x) 2 x −1 has the maximum g(e 2 ) = 4e −2 . Consequently, Hence Borel-Cantelli implies which by (26) shows (24). Let us next show that (4) satisfies (2): This shows that (4) is a solution. It remains to prove the uniqueness. By (3) we have for every t ∈ Z and k ∈ {0, 1, . . .} that Suppose now that there exists two solutions σ 2 t andσ 2 t satisfying sup t∈Z E(σ 2 t ) < ∞ and sup t∈Z E(σ 2 t ) < ∞. Then As both terms on the right-side converges in L 1 to zero by Lemma 5, we observe that for all t ∈ Z which implies the result.

Proof of Corollary 1 Let k be fixed and define
If L is strictly stationary, then so is (A t , B t ). Consequently, we have Since the limits of the both sides exist as k → ∞ we have Treating multidimensional distributions similarly concludes the proof.

Proof of Lemma 1
The "if" part follows from Theorem 3.2 of Karlsen (1990). For the converse, denote E( 4 0 ) = C and E(L 2 0 ) = C L . By (4) we have and since all the terms above are positive, both sides are simultaneously finite or infinite. Note also that, as the terms all positive, we may apply Tonelli's theorem to change the order of summation and integration obtaining Consider now the first term above. By independence, we obtain Consequently, E(σ 4 0 ) < ∞ implies α 1 < 1 √ C , since it is the radius of convergence of the series above.

Proof of Corollary 2 Let H k,t be defined by (27). By the proof of Lemma 1 we get
Thus H k,t converges to σ 2 t+1 in L 2 . To conclude the proof, it is straightforward to check that weak stationarity of L implies weak stationarity of H k,t for every k.

Proof of Lemma 2 First we notice that
by Lemma 1. Hence, the stationary processes Y and Z have finite second moments. Furthermore, the covariance of Y coincides with the one of X 2 . Applying Lemma 1 of Voutilainen et al. (2017) we get for every n ∈ Z, where r (·) is the autocovariance function of Z . For r (n) with n ≥ 1 we obtain since the sequences ( t ) t∈Z and (L t ) t∈Z are independent of each other, and t is independent of σ s for s ≤ t. By the same arguments, for n = 0 we have Now using (28) and γ (−1) = γ (1) completes the proof.
In the remaining we denote Lemma 6 Let t, s ∈ Z. Then Proof By (4) and Fubini-Tonelli where the series converges since α 1 < 1 and E(L 2 0 ) < ∞.
We use the following variant of the law of large numbers for consistency of the covariance estimators. The proof based on Chebyshev's inequality is rather standard and hence omitted. In the case of a weakly stationary sequence, a proof can be found from (Kreiß and Neuhaus 2006, p. 154), or from (Lindgren 2012, p. 65) concerning the continuous time setting.
Lemma 7 Let (U 1 , U 2 , . . .) be a sequence of random variables with a mutual expectation. In addition, assume that Proof of Lemma 3 By Lemma 7 it suffices to show that Cov(X 2 0 X 2 n , X 2 t X 2 t+n ) converges to zero as t tends to infinity. Hence we assume that t > n. By (4) Since the summands are non-negative, we can take the expectation inside. Furthermore, by independence of the sequences t and L t we observe Next we justify the use of the dominated convergence theorem in order to change the order of the summations and taking the limit. Consequently, it suffices to study the limits of the terms Step 1: finding summable upper bound. First note that the latter term is bounded by a constant. Indeed, by stationarity of (B t ) t∈Z we can write which is bounded by a repeated application of Cauchy-Schwarz inequality and the fact that the fourth moment of L 0 is finite. Consider now the first term in (31). First we recall the elementary fact Next note that the first term in (31) is bounded for every set of indices. Indeed, this follows from the independence of and the observation that we obtain terms up to power 8 at most. That is, terms of form 8 t and by assumption, By computing similarly for n = 0, using stationarity of A and Eq. (33), we deduce that where C is a constant. Moreover, by using similar arguments we observe Combining all the estimates above, it thus suffices to prove that and a 3 = α 1 E( 4 0 ).
Then we need to show that For this suppose first that 1 / ∈ S: ={a 1 , a 2 , a 3 , a 1 a 2 , a 2 a 3 , a 1 a 2 a 3 }. Then we are able to use geometric sums to obtain Continuing like this in the iterated sums in (34) we deduce Consequently, it suffices that the following three series converge . Finally, if 1 ∈ S it simply suffices to replace a 1 , a 2 , a 3 with such that 1 / ∈ {ã 1 ,ã 2 ,ã 3 ,ã 1ã2 ,ã 2ã3 ,ã 1ã2ã3 }.
Choosing δ < 0 small enough the claim follows from the fact that the inequality α 1 < 1 E( 8 0 ) 1 4 is strict.
Consequently, we conclude that lim t→∞ E(X 2 0 X 2 n X 2 t X 2 t+n ) = E(X 2 0 X 2 n ) 2 proving the claim.

Proof of Theorem 2
Since the assumptions of Lemma 3 are satisfied, so are the assumptions of Lemma 4 implying that the autocovariance estimators, the mean and the second moment estimator of X 2 t are consistent. The claim follows from the continuous mapping theorem.

Appendix B: Tables
In Table 1 we have presented means and standard deviations of the estimates with different sample sizes. In addition, we have provided Table 2 demonstrating how the estimates match their theoretical intervals 0 ≤ α 0 , 0 < α 1 < 105 − 1 4 and 0 < l 1 . We can see that multiplying the mean squared error (RMSE) provided by Table 1 with N H , the power H of the sample size, gives us evidence of the convergence rates of the estimators.