Expectiles for subordinated Gaussian processes with applications

In this paper, we introduce a new class of estimators of the Hurst exponent of the fractional Brownian motion (fBm) process. These estimators are based on sample expectiles of discrete variations of a sample path of the fBm process. In order to derive the statistical properties of the proposed estimators, we establish asymptotic results for sample expectiles of subordinated stationary Gaussian processes with unit variance and correlation function satisfying $\rho(i)\sim \kappa|i|^{-\alpha}$ ($\kappa\in \RR$) with $\alpha>0$. Via a simulation study, we demonstrate the relevance of the expectile-based estimation method and show that the suggested estimators are more robust to data rounding than their sample quantile-based counterparts.


Introduction
In the statistic literature, there has been a tremendous interest in analysis, estimation and simulation issues pertaining to the fractional Brownian motion (fBm) (Mandelbrot and Ness, 1968). This is due to the fact that the fBm process offers an adequate modeling framework for nonstationary self-similar stochastic processes with stationary increments and can be used to model stochastic phenomena relating to various fields of research. A fractional Brownian motion (fBm), denoted {B H (t), t ∈ R} with Hurst exponent 0 < H < 1, is a zero-mean continuous-time Gaussian stochastic process whose correlation function satisfies E[B H (t)B H (s)] = σ 2 2 (|t| 2H + |s| 2H − |t − s| 2H ) for all pairs (t, s) ∈ R × R and σ 2 = E(B H (1) 2 ). The fBm is H-self-similar i.e., for all α > 0, B H (αt) d = α H B H (t), where d = means the equality of all its finite-dimensional probability distributions. The process corresponding to the first-order increments of the fBm is known as the fractional Gaussian noise (fGn) whose correlation function ρ H (i) is asymptotically of the order of |i| 2H−2 for large lag lengths i. In particular, for 1/2 < H < 1, the correlations are not summable, i.e. +∞ i=−∞ ρ H (i) = ∞. This property is referred to as long-range dependence or long-memory whereas the case 0 < H < 1/2 corresponds to short memory.
Several methods aimed at estimating the Hurst characteristic exponent or long-memory exponent have been developed. Among these statistical methods figure the Fourier-based methods such as the Whittle maximum likelihood estimator (see e.g. Beran (1994); Robinson (1995)) or the spectral regression based estimator (Beran, 1994). The wavelet estimators have been also extensively investigated either with an ordinary least squares (Flandrin, 1992), a weighted least squares (see e.g. Abry et al. (2000Abry et al. ( , 2003; Bardet et al. (2000); Soltani et al. (2004)) or a maximum likelihood (see e.g. (Wornell and Oppenheim, 1992;Percival and Walden, 2000)) estimation schemes. Faÿ et al. (2009) present a deep analysis of Fourier and wavelet methods. Recently, the so-called discrete variations techniques (see e.g. Kent and Wood (1997); Istas and Lang (1997); Coeurjolly (2001)) have been introduced. Within this class of estimators, Coeurjolly (2008) proposed a new method based on sample quantiles to estimate the Hurst exponent in the more general setting of locally self-similar gaussian processes. This estimator has been proven robust when dealing with outliers (Achard and Coeurjolly, 2010). The latter, often encountered in real world applications, can induce a significant estimation bias. Actually, the advantage of quantiles is that they have a bounded gross-error-sensitivity (Hampel et al., 1986;Huber, 1981) allowing them to cope efficiently with the problem of outliers. Nevertheless, this is not the only problem faced when dealing with estimation issues. Indeed, data rounding is also a serious impediment.
Data rounding is common in finance (Bijwaard and Franses, 2009;Rosenbaum, 2009), economics (Williams, 2006), computer science (Matthieu, 2006;Bois and Vignes, 1982) and computational physics (Vilmart, 2008) and can lead to several misinterpretations. Quantiles are unfortunately not robust against data rounding since their local shift sensitivity is unbounded, see Hampel et al. (1986);Huber (1981). Newey and Powell (1987) have introduced the so-called expectile which, although similar to quantile, has a bounded local shift-sensitivity and thus can handle the rounding issue.
In this paper, we derive a Bahadur-type representation for sample expectiles of a subordinated Gaussian process with unit variance and correlation function with hyperbolic decay. This allows us to investigate the statistical properties of a new discrete variations estimator of the Hurst exponent of the fBm process. In constructing this estimator, we rely mainly on the scale and location equivariance properties of expectiles (Newey and Powell, 1987). We will show via a simulation study the robustness of the proposed estimator against data rounding.
The remainder of this paper is structured as follows: Section 2 deals with asymptotic properties of sample expectiles for a class of subordinated stationary Gaussian processes with unit variance and correlation function satisfying ρ(i) ∼ κ|i| −α (κ ∈ R) with α > 0. A short simulation study is conducted to corroborate our theoretical findings. In Section 3, we discuss a sample expectile-based estimator of the Hurst exponent and derive its statistical properties. We then perform a simulation study in order to confirm the effectiveness of the suggested estimation method.
2 Expectiles for subordinated Gaussian processes 2.1 A few notation Given some random variable Z with mean µ, F Z is referred to the cumulative distribution function of Z and ξ Z (p) for p ∈ (0, 1) to its pth quantile. It is well-known that the pth quantile of a random variable Z can be obtained by minimizing asymmetrically the weighted mean absolute deviation ξ Z (p) := argmin θ E |p − 1 Z≤θ |.|Z − θ| .
In order to limit the local shift sensitivity of the pth quantile, Newey and Powell (1987) defined the notion of expectile denoted by E Z (p) for some p ∈ (0, 1). Rather than an absolute deviation (function), a quadratic loss function is considered: (1) We may note that the 50%-expectile if nothing else than the expectation of Z. Newey and Powell (1987) argued that providing E[Z] < +∞, then for every p ∈ (0, 1) the solution of (1) is unique on the set I FZ := {x ∈ R : F Z (x) ∈ (0, 1)}. The expectile can also be defined as the solution of A key property of the expectile is that it is scale and location equivariant (Newey and Powell, 1987). The scale equivariance property means that for Y = aZ where a > 0, the pth expectile of Y satisfies: The pth expectile is location equivariant in the sense that for Y = Z + b where b ∈ R, the pth expectile of Y is such that: Now, let Z = (Z 1 , . . . , Z n ) be a sample of identically distributed random variables with common distribution F Z , the sample expectile of order p is defined as:

Main result
In order to derive asymptotic results for Hurst exponent estimates based on expectiles, we have to provide asymptotic results for sample expectiles of nonlinear functions of (centered) subordinated stationary Gaussian processes with variance 1 and with correlation function decreasing hyperbolically. This will be the setting of the rest of this section. Let {Y i } +∞ i=1 be such a Gaussian process with correlation function ρ(·) satisfying ρ(i) ∼ κ|i| −α for κ ∈ R and α > 0. Let Y = (Y 1 , . . . , Y n ) a sample of n observations and h(Y) = (h(Y 1 ), . . . , h(Y n )) its subordinated version for some measurable function h. We wish to provide asymptotic results for the sample pth expectile defined by Since the criterion is differentiable in θ, the sample pth expectile also satisfies the following In the following, we need the two following additional notation for Y ∼ N (0, 1) the latter quantity corresponding to the derivative of ψ h(Y ) (·, p) if it is well-defined. Let us note that the pth expectile of h(Y ) satisfies ψ h(Y ) (E h(Y ) ; p) = 0. We now present the assumption on the function h considered in our asymptotic result.
[A(h,p)] h(·) is a measurable function such that Eh(Y ) 2 < +∞ and such that the function Such an assumption is in particular satisfied under the following one: For the purpose of this paper, our main result will be applied with h(·) = | · | β (with β > 0) or The nature of the asymptotic result will depend on the correlation structure of the Gaussian process and on the Hermite rank, τ (p, θ) of the function We recall that the Hermite rank (see e.g. Taqqu (1977)) corresponds to the smallest integer such that the coefficient in the Hermite expansion of the considered function is not zero. For the sake of simplicity, assume that the Hermite rank of this function depends neither on θ nor p and denote it simply by τ . Again, this could be weakened since we believe that the next result could be proved with the following Hermit rank: inf θ∈V(E h(Y );p ) τ (p, θ). As an example, the Hermite We now present our main result stating a Bahadur type representation for the sample pth expectile of a subordinated Gaussian process.
a (centered) stationary Gaussian process with variance 1 and correlation function satisfying ρ(i) ∼ κ|i| −α (κ ∈ R), as |i| → +∞ with α > 0 and with a function h where the sequence r n = r n (α, τ ) is defined by Proof. Let us simplify the notation for sake of conciseness: The first thing to note is that the sequence r n corresponds to the short-range or long-range dependence characteristic of the sequence ψ(h(Y 1 ); p, θ), . . . , ψ(h(Y n ); p, θ). More precisely r 2 n corresponds to the asymptotic behavior of Eψ n (E) 2 . Indeed, if (c j ) j≥0 denotes the sequence of the Hermite coefficients of the expansion of ψ(·; p, E) in Hermite polynomials (denoted by (H j (t)) j≥0 and normalized in such , we may have using standard developments on Hermite polynomials (see e.g. Taqqu (1977)) Let us define We just have to prove that V n − W n (E) converges in probability to 0 as n → +∞. The proof is based on the application of Lemma 1 of Ghosh (1971) which consists in satisfying the two following conditions: We consider only the first limit. The second one follows similar developments. We first state that the map ψ n (·) is decreasing. Indeed, let θ ≤ θ ′ and denote by Z i (θ) the variable to the decreasing of ψ n (·). And in the in between case, which leads to the same conclusion. Let y ∈ R, then also using the fact that ψ n ( E) = ψ(E) = 0 and ψ ′ (E) < 0, we derive . Under the assumption [A(h,p)], y n → y as n → +∞. Now, let U n := ψ ′ (E) (W n (E) − W n (y × r n + E)), explicitly given by Let c j,n the jth Hermite coefficient of the function r −1 n ψ(t; p, E + y × r n ) − ψ(t; p, E) , then under the assumption [A(h,p)] and from the dominated convergence theorem we can prove that Therefore for n large enough, which leads to the convergence of U n to 0 in probability. For all ε > 0, there exists n 0 (ε) such that for all n ≥ n 0 (ε), y n ≤ y + ε/2. Therefore for n ≥ n 0 (ε) which ends the proof.
In the case of short-range dependence, i.e. ατ > 1 then, using the Bahadur type representation of expectiles, we derive immediately the following asymptotic normality for the sample expectile and some generalisations. This result is based on standard central limit theorem for means of subordinated Gaussian stationary processes (Taqqu, 1977;Arcones, 1994).Therefore the proof is omitted.

Corollary 2
(i) Under the assumptions of Theorem 1 with p ∈ (0, 1) and ατ > 1, then as n → +∞ and where c k (p) is the kth Hermite coefficient of the expansion of the function ψ(h(·); E h(Y ) (p); p) in Hermite polynomials.

Simulations
To illustrate a part of the previous results, we propose a short simulation study in this section.
The latent stationary Gaussian process we consider here is the fractional Gaussian noise with variance 1, which is obtained by taking the discretized increments from a fractional Brownian motion. The correlation function of the fractional Gaussian noise with Hurst parameter (or self-similarity parameter) H ∈ (0, 1) satisfies the hyperbolic decreasing property required in Theorem 1 with α = 2 − 2H. Discretized sample paths of fractional Brownian motion can be generated exactly using the embedding circulant matrix method popularized by Wood and Chan (1994) (see also Coeurjolly (2000)) which is implemented in the R package dvfBm.
Figures 1 and 2 illustrate the convergence of the sample expectiles. Three h functions are considered: h(·) = (·), (·) 2 and log | · |. The related Hermite rank of the function ψ is respectively 1,2 and 2 for these three h functions. The sample size of the simulation is fixed to n = 500. We can claim the convergence of the sample expectile E(p; h(Y )) towards E h(Y ) (p) for all the values of α (or H), p and for the three functions h considered. If we focus on h(·) = (·), we can also remark a higher variance of the sample estimates for α = 0.6 compared to α = 1.4. This is in agreement with the theory since for α = 0.6, ατ = 0.6 < 1 and the rate of convergence is lower than n −1/2 which means an increasing of the variance. For the two other functions considered, then ατ is always greater than 1 (it equals either 2.8 or 1.2 in our simulations) and we do not observe such an increasing of the variance.
To put emphasis on this last point, Figure 3 shows in log-scale the average (over the 9 order of expectiles considered in the simulation, i.e. p = 0.1, . . . , 0.9) of the empirical variances in terms of n for the three h functions and for the two values of α = 0.6 and α = 1.4. We clearly observe that as soon as ατ > 1, the slope of the curves is close to −1 which agrees with the result presented in Corollary 2 for example. When h(·) = 1 and α = 0.6, we observe that the slope is about −0.6 which seems to agree with the convergence in n −ατ which is expected from Theorem 1. 3 Estimation of the Hurst exponent using sample expectiles and discrete variations

Estimation method and asymptotic results
Let X = (X(i)) i=1,...,n be a discretized version of a fractional Brownian motion process and let a be a filter of length ℓ + 1 and of order ν ≥ 1 with real components i.e.: ℓ q=0 q j a q = 0, for j = 0, . . . , ν − 1 and ℓ q=0 q ν a q = 0.
Define also X a to be the series obtained by filtering X with a, then:  the orders of the expectiles and we compute σ 2 n = 1/9 × andX a as the normalized vector of X a , i.e.: X a = X a E((X a (1)) 2 ) 1/2 .
It should be noticed here that the filtering operation allows to decorrelate the increments of the discretized version of the fractional Brownian motion process. Indeed, it may be proved (see e.g. Coeurjolly (2001)) that: ρ a H (i) ∼ k H |i| 2H−2ν as |i| → +∞. Consider the sequence (a m ) m≥1 defined by: which is the filter a dilated m times. It has been shown in Coeurjolly (2001Coeurjolly ( , 2008 that: The following proposition allows us to construct an ordinary least squares (OLS) estimator of the Hurst exponent H of a fBm process based on sample expectiles.
Proposition 3 Let E p; h(X a m ) and E p; h(X a m ) be the pth order sample expectiles for the filtered series h(X a m ) and h(X a m ) respectively. Here two positive functions h(·) are considered, namely: h(·) = | · | β for β > 0 and h(·) = log | · |. We have: and E p; log |X a m | = 1 2 log(σ 2 m ) + E p; log |X a m | . Proof.
We have: Setting θ ′ = θ σ β m , the proof of the first relation (9) follows easily. Using the same methodology, we can demonstrate the result given by equation (10).
Remark 1 It should be stressed here that the scaling relationship relating the theoretical pth expectiles for the series h(X a m ) and h(X a m ) can be obtained directly using the scale equivariance property (2) for h(·) = | · | β and the location equivariance property (3) for h(·) = log | · |.
Now applying the logarithmic transformation to both sides of (9), we get: On the other hand, (10) can be reformulated in the following way: It is noteworthy here that we expect that log E p; |X a m | β /E |Y | β (p) and E p; log(|X a m | β ) − E log |Y | (p) to converge towards 0 as n → ∞. Hence, based on equations (11) and (12), we opt for an OLS regression scheme. This allows to derive the two following estimators of the hurst index defined by: and where We would like to put the stress on the fact that (13) and (14) are really similar to the ones developed in Coeurjolly (2001Coeurjolly ( , 2008. Indeed, the standard procedure developed in Coeurjolly (2001) simply consists in replacing the sample expectile by the sample variance (this method will be denoted by ST in Section 3.2). To deal with outliers, the procedure developed in Coeurjolly (2008) consists in replacing the sample expectile by either the sample median of (X a m ) 2 or the trimmed-means of (X a m ) 2 . These two last methods are denoted by MED and TM in Section 3.2. Now, we state the asymptotic results for these new estimates based on expectiles.
Proposition 4 Let a a filter with order ν ≥ 2, p ∈ (0, 1), β > 0 then as n → +∞, H β and H log converge in probability to H. Moreover, the following convergences in distribution hold and where the (M, M ) matrices Σ β and Σ log are defined by (8).
Proof. We only provide a sketch of the proof. We claim that once Theorem 1 and Corollary 2 are established, the obtention of convergences stated in Proposition 4 are semi-routine. First of all, let us notice that and Since the functions | · | β and log | · | have Hermite rank 2 and since the correlation function of the stationary sequence X a m decreases hyperbolically with an exponent α = 2ν − 2H then for any m ∈ {1, . . . , M }, Theorem 1 holds with r n = n −1/2 (since ατ > 1 for all H ∈ (0, 1)). This ensures the convergence in probability of the new estimates.
The cross-correlation between X a m 1 and X a m 2 is defined by a q a r |m 1 q − m 2 r + j| 2H .
Lemma 1 in Coeurjolly (2008) states that for all m 1 , m 2 the correlation function ρ a m 1 ,a m 2 H is also decreasing hyperbolically with an exponent α = 2ν − 2H, then Corollary 2 may be applied to prove that and where according to (8), the (M, M ) matrices Σ β and Σ log are respectively defined by where (c β k ) k≥2 and (c log k ) k≥2 are respectively the Hermite coefficients of the functions | · | β and log | · |. The convergences (17) and (18) combined with (15) and (16) and the use of the deltamethod (for the convergence of H β ) end the proof.

A short simulation study
In this section, we investigate the interest of the new estimators based on expectiles. We consider three different models in our simulations.
(b) fBm with additive outliers: we contaminate 5% of the observations of the increments of the fractional Brownian motion with an independent Gaussian noise such that the SNR of the considered components equals −20Db.
(c) rounded fBm: we assume the data are given by the integer part of a discretized sample path of an original fBm.
To fix ideas, Figure 4 provides some examples of discretized sample paths of standard and contaminated fBm. The simulation results are presented in Tables 1 and 2. For these simulations, as suggested in Coeurjolly (2001), we chose the filter a = d4 corresponding to the wavelet Daubechies filter with order 4 (see Daubechies (1992)) and the maximum number of dilated filters M = 5. Also, in other simulations not presented here, we have observed that the estimates H β perform better than H log and, among all possible choices of β, the value β = 2 seems to be a good compromise. Therefore, we present only the result for this latter estimator, that is H β with β = 2.
In a first step, we had observed a quite large sensitivity to the value of the probability p defining the expectile. In order to have an efficient data-driven procedure, we propose to choose the probability parameter p via a Monte-Carlo approach as follows: 1. Estimate the parameters H and σ 2 using the standard method (the estimation of σ 2 is not described here but it may be found for example in Coeurjolly (2001)). Denote these estimates H 0 and σ 2 0 .
2. Generate B = 100 contaminated fBm with Hurst parameter H 0 and scaling coefficient σ 2 0 , define a grid of probabilities (p 1 , . . . , p P ). For each new replication, we estimate H 0 with expectiles for all the p i . The optimal p, denoted in the tables by p opt , is then defined as the one achieving the smallest mean squared error (based on the B = 100 replications).
The procedure based on expectiles, denoted E(p) in the results, is compared to the standard method (ST) and to methods which efficiently deal with outliers, that is methods MED and TM (the last one is calculated by discarding 5% of the lowest and the highest values of (X a m ) 2 at each scale m).
The standard fBm model is used as a control to show that all methods perform well. As seen in the first two columns of Tables 1 and 2, this is indeed true. All the methods seem to be asymptotically unbiased and have a variance converging to zero. We can also remark that in this situation whatever the value of H, estimates based on expectiles exhibit a performance which is very close to the one of the standard method (wich can also be viewed as the method based on expectile with p = 0.5). Several types of expectiles are investigated. When the discretized sample path of the fBm is contaminated by outliers, we recover the results already shown in Coeurjolly (2008), Achard and Coeurjolly (2010) or Kouamo et al. (2010): methods based on medians or trimmed-means are very efficient which is in agreement with the fact that quantiles have a finite gross error sensitivity. The inefficiency of expectiles for such a contamination is also coherent since expectiles have infinite gross error sensitivity. Finally, the interest of the expectile-based method can be seen with the rounded fBm corresponding to the last two columns of each table.
In this situation, expectiles are shown to be more efficient in terms of bias and its variance seems to be not too much affected by this type of strong contamination. We also put the stress on the interest and efficiency to choose the p value based on a Monte-Carlo approach. Hurst parameter H = 0.2 and sample size n = 500, 5000 are given between brackets. Methods based on expectiles, quantiles and trimmed-means as well as the standard method are considered.
The filter a correspond to the Daubechies wavelet filter with order 4 (two zero moments) and we set M 1 = 1, M 2 = 5. According to a sample size and a model, the method achieving the lowest mean squared error is printed in bold.
Standard fBm fBm with additive outliers Rounded fBm n = 500 n = 5000 n = 500 n = 5000 n = 500 n = 5000  Hurst parameter H = 0.8 and sample size n = 500, 5000 are given between brackets. Methods based on expectiles, quantiles and trimmed-means as well as the standard method are considered.
The filter a correspond to the Daubechies wavelet filter with order 4 (two zero moments) and we set M 1 = 1, M 2 = 5. According to a sample size and a model, the method achieving the lowest mean squared error is printed in bold.