Ji-Ming Penga

Christian Kanzowb
email: kanzow@@math.uni-hamburg.de.

Masao Fukushimac

A HYBRID JOSEPHY-NEWTON METHOD FOR SOLVING BOX CONSTRAINED VARIATIONAL INEQUALITY PROBLEMS VIA THE D-GAP FUNCTION

Keywords: Variational inequality problem, box constraints, D-gap function, Newton's method, unconstrained optimization, global convergence, quadratic convergence

A box constrained variational inequality problem can be reformulated as an unconstrained minimization problem through the D-gap function. Some basic properties of the affine variational inequality subproblems in the classical Josephy-Newton method are studied. A hybrid Josephy-Newton method is then proposed for minimizing the D-gap function. Under suitable conditions, the algorithm is shown to be globally convergent and locally quadratically convergent. Some numerical results are also presented.

INTRODUCTION

Let F be a mapping fron1 �R'' into itself and X be a nonempty closed convex subset of R11• The variational inequality problem (VIP) is to fjnd a vector x* EX such that (F(xx)Sx*) > 0, 'iy E X,

( 1)

where ( •, •) denotes the inner product in �11• rn this paper, we study the case that X is a box defined by X= { x E 1R" ili < x; < u;, i = 1, .... , n}

where I; E 3t U { -x} and 11i E � U { +oo} with li < u1 representing the lower and upper bounds on thevariables, respectively. It is not difficult to see [START_REF] Bertsekas | Parallel and Distributed Computation Numerical lvfetlrod[END_REF]Proposition 5. 7] that the box structure of the set X enables us to rewrite (1) as F,.(x*)(V;-x;) > 0, i = 1, ... , n, 'iy EX.

( 2 )

If the constraint set X is the nonnegative orthant ��' then the VIP reduces to the complemelitarity problem (CP). This class of special VIPs has numerous itnportant applications in various fields such as mathe matical programming, economics, and engineering; see [7,1 3] and references therein.

A useful way to deal with VIP (1) is to reformulate it first as a system of equations or an optitnization problen1 via a n1erit function, and then solve the resulting systen1 of equations or optimization problem. Recently.tnuch attention has been paid to the reforn1ulation ofVIPs and CPs and vari ous merit functions have been proposed and studied [START_REF] Fuk Ushima | Merit fu nctions fo r variational inequality and complementarity problems[END_REF].

Well known merit functions for VIPs include the gap.ftmction g(x) = sup (F(x), x-y) r�X first presented by A us! ender [START_REF] Auslcnder | Optimisation: Mi!Jwdes Numhiques[END_REF] and then studied by M arcotte and Dussault [START_REF] Marcotte | A new algorithm f or solving variational inequalities with application to the t raffic assignment problem[END_REF][START_REF] Marcotte | A note on a globally convergent Newton method fo r solvin g monotone variational inequalities[END_REF], the regulari�ed gapjzmction

( 1 " -)

introduced by Fukushima [START_REF] Fuk Ushima | Equivalent differentiable optimization problems and descen t methods fo r asymmetric variational inequality problems[END_REF] and Auclunuty [1 ], and the D-gap function (4)

proposed by Peng [START_REF] Peng | Equivalence of variational inequality problems to unconstrained mini mization[END_REF] and Yamashita et al. [27], where a and (3 are arbitrary positive parameters such that a < f-1. Since a box constrained VIP is actually equivalent to a system of KKT rnixed con1plen1entarity conditions, several merit functions based on the KKT systen1 of VIP have also been proposed and explored extensively, see Qi [START_REF] Qi | Regular pseudo-smooth NCP and BVIP functions and globally and quad ratically convergent generalized Newton me thods for complcmentarity and vari ati o nal inequali t y problems[END_REF] and references therein.

In this paper, we focus our attention on the D-gap function for VIP.

It is not difficult to see that gu,s{x) > 0 for all .rE R'\ and gn / �(x) = 0 if and only ifx is a solution of the VIP (1 ) . Therefore the VIP can be cast as the following unconstrained optin1ization problem: min go. a (x).

xEiR" .

(

) 5
When the mapping F is differentiable, the D-gap function ga(l is also differentiable [START_REF] Fuk Ushima | Equivalent differentiable optimization problems and descen t methods fo r asymmetric variational inequality problems[END_REF]27]. However, it is not twice differentiable in general . Therefore it is not straightforward to apply conventional second-order 1nethods toproblem [START_REF] Facchinei | Regularity properties of a scmismoot.h reformulation of v�triational inequalities[END_REF]. Asarernedy for this inconvenience, Sun et al. [2 5 ] introduced the concept of a computable generalized Hessian of the D-gap function g03 and presented a Newton-type method for solving problen1 [START_REF] Facchinei | Regularity properties of a scmismoot.h reformulation of v�triational inequalities[END_REF]. Restricting themselves to the box constrained VIP, Kanzow and Fukushin1a [START_REF] Ka | Theoretical and numerical investigation of the D gap function fo r box constrained variational ineq ua lities[END_REF] discussed a generalized Hessian of the D-gap function and proposed a Gauss-Newton-type method to minimize it. Further properties of the D-gap function have been investigated in [START_REF] Fukushima | Minimizing and sta tionary sequences of merit func tions fo r complementarity problems and vari a tional inequalities[END_REF][START_REF] Ka | Theoretical and numerical investigation of the D gap function fo r box constrained variational ineq ua lities[END_REF][START_REF] Sun | A computable generalized Hessian of the D gap function and Newton-type methods fo r variational inequality problems[END_REF].

This work is n1otivatcd by the recent paper [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] in which the D-gap function is used to globalize the classical Josephy-Newton method [START_REF] Josephy | Newto n' s method for generalized equations[END_REF] for general VIPs. The n1ain purpose of this work is to further study the algorithm proposed in [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] by restricting ourselves to box constrained VIPs, and test the effectiveness of the algorithm. In [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] , the global and local convergence of the algorithm were proven under the assumptions that Fis strongly n1onotone. By restricting to box constrained VIP, we will show in this paper that the algorithm is globall y convergent and locally q uadratica lly convergent if the mapping Fis a unifonn P-function. Hence the results here can also be viewed as a refinement of the results in [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF].

The paper is organized as follows: In Section 2, we first introduce some basic definitions that will be used in the paper. Son1e sensi tivity results fo r box constrained affine VIP are given under certain condi tions. Then, by applying these sensitivity results, we study the properties of the subproblem in theJosephy-Newtonmethod [START_REF] Josephy | Newto n' s method for generalized equations[END_REF].ln Section 3, we present a hybrid algorithrn which uses the D-gap function to globalize Josephy-Newton method. Convergence properties of the algorithm are also discussed. In Section 4 we report some numerical results. Finally we conclude the paper with some rernarks in Section 5.

SOME RESULTS FOR BOX CONSTRAINED AFFINE VIP

We first review some concepts related to the mapping F. The mapping F: �n---+ ?R'1 is said to be a P0-jimction if

m�x(xi -yt)(Fi(."\:) -F i(V)) > 0 Vx, y E rrt11,
It is known [START_REF] More | On P-and S-functions and related classes of n-dimensional nonlinear mappings[END_REF]Theorem 5.8] that if F is a differentiable P0-function, then V' F(x) T is a P0-matrix for each x. Moreover, if Fis a differentiable uniform P-function with modulus I'•> 0, th e n V' F(.x) T is a uniform P-n1atrix with modulus p. > 0 in the sense that Since we are not aware of any explicit reference for the formula (7), we include a short proof for it.

LEMMA 2.1 Let F: ��____,. �n be a d?fferentiable un{form P-function with modulus Jl. > 0. Then (7) holds.

Proof Let x E �� and z E �n be arbitrary but fixed. Let { tk} C R be a sequence of positive numbers converging to 0. Since F is a uniform P-function and the index set { 1, . .. , n} is finite, there exists an index i0 E { 1, ... , n} (independent of k) and a subsequence { tk} K such that zin i-0 and Dividing this expression by t�, taking the limit k � oo (k E K) and using the assumed differentiability ofF, we obtain This implies and completes the proof.

In the classical Josephy-Newton method, we solv e the fo11owing linearized VIP subproblem at each iteration: Given a current iterate x\ find a vector x EX such that

Here we study the properties of the affine subproblem [START_REF] Ferris | Accessing realistic mixed complementarity pro blems within MATLAB[END_REF]. In [START_REF] Taji | A glo bally convergent ;\.lcwton method for solving strongly monotone variational inequalities[END_REF], Taji et al. studied the subproblem (8) for a general closed convex constraint set X under the assumption that F is strongly monotone and con tinuously differentiable. Restricting ourselves to the box constrained VIP, we will investigate properties of the subproblem under weaker conditions.

To this end. we first consider a perturbation of affine VIP with box constraints. Let X be a box, MaP-matrix, and bE ?R". For any v E �11, let x(v) EX denote the unique solution of the following affine VIP:

( h + v + Mx, y -x) > 0� or equivalently Vv EX, o/ •, (9)
(b; +vi+ [Mx]J(,r;-x1) > 0, i = 1, . . . , n, Vy EX. ( 10
)
PROPOSITION 2.2 Let x( v) be the unique solution q[ the ({fflne VIP [START_REF] Fuk Ushima | Equivalent differentiable optimization problems and descen t methods fo r asymmetric variational inequality problems[END_REF],

where X fč a box and M a P-matrix. Let 11 > 0 be a constant satisfmfČg [START_REF] Facchinei | A semismooth Newton method fo r var iational inequalities: The case of box constraints[END_REF]. Then we have llx(v)-x(v')ll < _!_llv-v'jj, '1/v, v' E R".

(hi + Vi + [ Mx (V)] i) (Xi (V 1) -X; (V)) > 0, i = 1 , ... , 11,
Adding the above two inequalities, we obtain

(v1-v/)(x;(v')-x;(v)) > (x1(v)-x1(v'))(j\4(x(v) -x (v'))Jr� i= 1, ... ,11. '11) ( . .
Since M is a P-matrix, there exists a constant fL > 0 and some index i0 such that

This inequality together with (1 1) implies that

Since (V io -V :0) (X iu (V 1) -X io (V) ) < 11 V -V 1 11 ' 11 X (V) -X (V 1)
1 l l x (v) -x (v1 ) 1 1 < -[[v -v ' l l • ' ' fl
This c01npletes the proof of the proposition.

Next we give another result on the continuity of a solution of the box constrained affine VIP with a P-rnatrix. For any P-matrix M and vectors b,p E �'\ let x (M) EX denote the unique solution of the following affine

VIP:

or equivalently

(b + Mp + Mx, y -x) > 0, \fy E X, (13)
( b ; + [M p ]; + [MxL)(y;-Xi) > 0, i = 1, ... , n, VEy E X.
( 1 4)

PROPOSITION 2.
3 Suppose that M and N are ?-matrices and X is a box. Let J1 > 0 be a constant satisfying (6). Then 1-ve have

1 l l x(M)-x(N)II < -IIM-Nl l •llx(N) + p[[. ( 15
)
fl Proof The proposition trivially holds if x(M) = x(N). Hence we only need to consider the case where x( M ) # x(N). Since X is a box, it follows from (14} and the definitions of x(M) and x(I\l) that and

Adding the above two inequalities, we get

([Mpl--[iVjJ]; + [Mx(M)L � [Nx(N)];)(x;(N)-xi(AJ)) > 0, i= l, . .. ,n,
which implies that

[(ill/-tV)(x(N) + p)]i(x;(N)-xi(�/)) > [i\J(x(Af)-x(N))]i(x;(M)-xi(iV)), i = L ... , 11. ( 16
)
Since M is a P-matrix, it follows fr01n. (6) and .x(M)-x(N) ::f. 0 that there exists an index i0 such that

Then ( 16) and ( 1 7) imply JL[[x(AJ)-x(N")[I 2 < [(Af-N)(x(N) + p)L0(Xi0(N)-x,0(AJ)) < IIM-Nil •[[x(N) + p J J• [[x(lv')-x(Af)jJ.
This c01npletes the proof. results based only on somewhat strong assun1ptions on the n1apping F. Now we return to the affine subproblem of Josephy -Newton method. Let x be a given point in �R'' and consider the 1inearized VIP of finding a point z EX such that (F(x) + \7F(x) T (z-x),y -z) > 0, Vy EX, [START_REF] Marcotte | A new algorithm f or solving variational inequalities with application to the t raffic assignment problem[END_REF] or equivalently

(Fi( x ) + [V'F( x )T(z-x)]J(J' ;-zi ) > 0 , i = 1, ... ,n, Vy EX. ( 19
)
If \7 F(x) is a P-tnatrix, then problem (18) has a unique solution, which we denote z(x). Our next result studies the properties of the solu tion of problem [START_REF] Marcotte | A new algorithm f or solving variational inequalities with application to the t raffic assignment problem[END_REF]. The proposition is a refinement of Proposition 2.2

in [START_REF] Taji | A glo bally convergent ;\.lcwton method for solving strongly monotone variational inequalities[END_REF].

PROPOSITION 2.4 Suppose that F is a continuously d�fferentiable uni form P-function and X is a box. Then the solution z (x) of the affine VIP [START_REF] Marcotte | A new algorithm f or solving variational inequalities with application to the t raffic assignment problem[END_REF] is continuous as a function r�l x. Moreover, x is a solution of the VIP (1) if and only ifx = z(x).

Proof First we note that, since F is a uniform P-function, ( 7) is satis fied with some constant 1-L > 0 independent of x. For two arbitrary points x, x' E �n, let z( x) and z(x') be the unique solutions of the linearized VIPs (18) at x and x', respectively. Also let z denote the unique solution of the affine VIP

(F(x) + \7F(x' ) T(z-x ) ,y-z) > 0, Vy EX. It then follows from Proposition 2.2 with v := F(x)-VH F(x ' ) T x, v':= F(x')-\7 F(x')Tx', b := 0 and M:=\7 F (x')T that 1 ll z -z(x') 11 < -II F (x) -F( x ' ) + \7 F( x') T (x' -x) 11•

IL

On the other hand, by Proposition 2.3 with b := F( x ), p := -x, N :=\7 F (x)T and M:=\7F( x')T, we have Consequently, for any fixed x E ���.we obtain lim ll.:(x) -.:(x') 11 = 0.

l l z-z(x) 11 < _!_ ll \7 F( x ') T -\7 F(x) T 11 • ll z(x) -x ii
x

'-x

This proves the first half of the proposition.

To prove the second half. suppose first that .:(x) = x. Then it follovv•s itnrnediately from (18) that x solves (1). Conversely suppose that xis a solution of ( 1 ). Since X is a box and z(x) EX. (

) 21
The inequalities ( 20) and (21) give (

) 22
Since VG F(x)T is a P-matrix. it fo llows fron1 (7) that there exists an index i0 . .

such that
Combining [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] with [START_REF] Qi | Regular pseudo-smooth NCP and BVIP functions and globally and quad ratically convergent generalized Newton me thods for complcmentarity and vari ati o nal inequali t y problems[END_REF]. we get .:(x) = x.

Regularity conditions have been widely used in the study of varia tional inequal1ty problems. patticu1ar1y in the analysis of IocaJ con vergence properties of iterative methods for VTPs. A solution x"' of the VIP (l) is said to be regular in the sense of Robinson [24] (see also [START_REF] Harker | Finite d imensional variational ineq uality and nonlinear complementarity problems: A survey of theory, algorithms and applications[END_REF])

if there exist a neighborhood n of x * and a neighborhood V of 0 E 3(.11 such that, for every v E V, the perturbed VIP of finding a vector x EX such that

(F(x; v),y-x} > 0, Vy EX, (24)
where F(x;v) := F(x*) + v + V F(-x:*)T (x-x*), has a unique solution

x (v) that is Lipschitz continuous as a function of v, i.e.,

for some p > 0.

ll x (v )-x (v')ll < Pllv-v'll , Vv , v' E V
In [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF], a strong but si1nple sufficient condition for a solution x* of the VIP to be regular was given. Here we give a different condition pertaining to the box constrained VIP. The next proposition follows directly fron1 Proposition 2.2. The proof is omitted. PROPOSITION 2.5 Assume that x* is a solution of the VIP (1) lvith X being a box. !f the matrix \7 F(x*) is a P-matrix, then x* is a regular solution.

Remark P roposition 2. 5 can also be derived from a characterization of strong regularity given in [6, Theorem 3.4]. The approach in [START_REF] Facchinei | A semismooth Newton method fo r var iational inequalities: The case of box constraints[END_REF] is based on the equivalent KKT system of box VIP and is different from our approach here.

The classical Josephy-Newton method generates a sequence {xk} by the updating rule xk+l = z(xk). From Proposition 2.4, the method is well defined if F is a uniform P-function. Moreover, if the initial point x0 is sufficiently close to the solution point x * and the matrix V F(x*) is a P-matrix, then by Proposition 2.5, x* is a regular solution and hence it follows fron1 the basic result of Josephy-Newton method [START_REF] Harker | Finite d imensional variational ineq uality and nonlinear complementarity problems: A survey of theory, algorithms and applications[END_REF][START_REF] Josephy | Newto n' s method for generalized equations[END_REF] that the generated sequence {x"} converges locally quadratically to x*.

A HYBRID JOSEPHY-NEWTON METHOD

In the previous section, we have discussed the affine VIP with box constraints. In this section, we consider the method proposed in [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] for solving the VIP (I) with general convex constraints, which is Jose phy-Newton method with D-gap function globalization. Our aim is to refine the convergence results obtained in [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF] by restricting ourselves to the special case where the VIP (1) is box constrained.

The algorithm is stated as follows:

ALGORITHM

Step 0 Choose x0 E �nl wE (0,1 ) , ( E (0, 1), bE (0� 1 ), u E (0,1), and sufficiently small E> 0. Let k: = 0.

Step 1 Step 2 If I J dkll < E then stop. Otherwise find the sn1allest nonnegative integer n1k satisfying (28) and let Ak := wmk.

Find :lE . X such that and let dk:= z k -.-•l. If ( .k dk) ( k) g ct/3 X + • < ( g a/3 X , ( 2 6
Step 3 Set)?' 1 :=xk+).. kcf and k:=k + I. Go to Step 1.

Using similar arguments to those in the proof of the global con vergence theorem in [START_REF] Peng | A hybrid Newton mdhod fo r solving the variational i nequality problem via the D-gap function[END_REF], one can prove the following theoretn. To study the convergence rate of the algorithm, we need the following results concerning an error bound property of the D-gap function. We denote by Ya(x) the unique maximizer on the right-hand side in the defining equation [START_REF] Bertsekas | Parallel and Distributed Computation Numerical lvfetlrod[END_REF] of the regularized gap function j�. Note that .Yu(x) = f1x( x-a-1 F(x)), where ITx denotes the projection operator on X. We define Ra(x) := x-Ya(x). Moreover Y.B(x) and R1:J(x) are defined similarly. Let B(�) denote the closed sphere centered at x* with radius .6. > 0, i.e., B ( .6.) := { x E � lll x -x *ll < .6.}.

Following the proof of Lemma 5.1 in [START_REF] Ka | Theoretical and numerical investigation of the D gap function fo r box constrained variational ineq ua lities[END_REF], we have the next lemma. LEMMA 3.3 Let x* be a solution of the VIP (1) with X being a box. Suppose that F is a uni f orm P-. function with modulus 1-l• Suppose also that F is Lipschitz continuous lvith constant K, > 0 on B(.6.) for some � > 0. Then there exists a constant rJ > 0 such that

llx-x*ll < TJ11R; 3( x )il, Vx E B(�), li"here 7] = (K + !3)/ p,. ( 29 )
The next lemma shows that the D-gap function provides a local error bound for the VIP (1) under suitable assumptions. This result will be useful in establishing the quadratic convergence of the proposed algorithm. Supp ose that F is a un�lorm P-jimctlon with modulus J.l• Suppose also that F is Lipschit:z continuous 1vith constant f£ > 0 on B(�) for some � > 0. Then there exist constants c1, c2 > 0 such that Proof By [27, Proposition 3.1]� we have It follows from Lemma 3.3 and the left part of the above inequality that where 77 = (1'£ + /3)/p. , which shows that the left inequality in (30) is true.

Next observe that ffRn(x)f l = l lx-Yo.(x)-x* + Yn(x*)ll < llx-x*ll + IT).Jx-a-1F(x))-Ilx(x* -a-1F(x*))

< llxx*ll + llx-o-1 F(x)x'i' + 0'�1F(x*)IJ < (2 +:)fi xx* ll for all x E B(A), where the equality fofšows from the de!lnition of R(l and the fact that Ro:(x*) = 0. the first inequality follows from the triangle inequality and the definition of Yt<• the second inequality follows from the nonexpansiveness of the projection operator IT x. and the last •inequality follows from the Lipschitz continuity of F. The right inequality in (30) then follows flom the right inequality of (31 ). The proof is con1plete.

'Ve arc ready to prove quadratic convergence of the proposed algorithm. ""' THEOREM 3.5 Suppose that F is continuousfv differentiable and X is a box. Let x* he an accumulation point of the sequence {xk} generated by the algorithm. If F is a uniform P-function and \7 F is locally Lipschit zian, then x* is a solution of the VIP (1) and the sequence {xk} converges quadratically to x*.

Proof By Theorem 3.1, x* is a stationary point of the D-gap func tion gct.B• Since \7 F(x*) is a P-n1atrix by Len1n1a 2.1, it follows from [15, Theorem 3.1] that x* is already a solution of the VIP (1). More over Proposition 2.5 shows that x* is a reg ular solution. Since F is differentiable, it is locally Lipschitzian; hence there exists a � 1 > 0 such that F is Lipschitz continuous on B(�t). Hence by Len1n1a 3 . 4, we have for some ch c2 > 0. Moreover, by choosing a smaller �1 if necessary, we may assume that \7 F is Lipschitz continuous on B(�1). Since x* is a regular solution, by the basic result on Josephy-Newton's n1ethod for the VIP [START_REF] Harker | Finite d imensional variational ineq uality and nonlinear complementarity problems: A survey of theory, algorithms and applications[END_REF][START_REF] Josephy | Newto n' s method for generalized equations[END_REF], there exists a �2 > 0 such that for any initial point chosen from B(�2), the Newton iteration is well-defined and holds for some constant c3 > 0. Let �3 := min(�1, �2). Then it follows frmn (32) and (33) that Let Then it follows fr01n (34) that for any x E B(�4)

where the first inequality follows from the choice of �4 and the second inequality follows frmn the left inequality in (32). This implies that, when xk E B(.6.4), we have d" =z(x")-xk and the step size ).. k =I is accepted, i.e., .r/' + 1 = z(xk). Consequently it follows from (33) that the sequence {xk} converges to x* quadratically.

Silnilarly we may prove superlinear conve rgence of the algorithm under slightly weaker assmnptions. The proof is omitted. THEOREM 3.6 Supp ose that F is continuousfţ d({ferentiable and X is a box. Let x* he an accumulation point of the sequence {xk} generated by the afşorithm. �!' F is a untforrn P-function� then x* is a solution qf the VIP (1) and the sequence {xk} converges superfŠearfŢ to x*.

NUMERICAL RESULTS

We implemented the hybrid Josephy� Newton method suggested in this paperinMATLABandrun it on a SUNSPARC 10station. We first give a brief description of the implementation: Let In addition, the iteration was stopped if k > kmttx with kmax-100.

For the D�gap function g0;:J, we used the panuneters a = 0.9 and /1 = 1.1. However, we replaced the standard (monotone) Armijo-rule }?y a nonmonotone variant, see Grippo et al. [START_REF] Grippo | A nonmonotone linesearch technique fo r Newton's method[END_REF] for details.

As a solver for the linearized variational inequality problems, we used the semismooth Newton-type method from [START_REF] Facchinei | A semismooth Newton method fo r var iational inequalities: The case of box constraints[END_REF]. In contrast to what is said in the description of our algorithm, however, we always accept the corresponding search direction dk whenever it satisfies the descent test note that this guarantees that the Arrnijo line search is well-defi ned.

In particular, we accept this search direction dk even if we were not able to solve the corresponding linearized variational inequality problem.

In this way, we try to overcome the problem that we have to take too many gradient steps in a row which is obviously not very desirable.

In order to improve the efficiency of our algorithm, however, we also used a preprocessor; more precisely, we first try to solve our test examples by using the recently proposed method from Kanzow and Fukushima [16]. This is a nonsmooth Newton-type method applied to the residual equation r (x) = 0 and globalized by the D-gap function go!;�, see [START_REF] Kanzow | Solving box constrained variational inequal ity problems by using the natural residual with D-gap function globalization[END_REF] for details. The motivation for doing this is quite simple: The method from [16] works extremely well whenever it solves a problem successfully. Unfortunately, it does not seem to be very robust unless relatively strong assumptions are satisfied.

So we first apply the nonsmooth Newton-type method fron1 [16] in order to solve a test example, but we stop this preprocessing iteration if either the termination criterion ( 35) is satisfied or if a certain test indi cates that the preprocessor runs into difficulties. In the latter case, we switch to the hybrid Newton method introduced in this paper which is not as efficient as the method from [16], but which seems to be con siderably more reliable.

Basically, our criterion for switching from the preprocessor to the hybrid Newton method is as follows: If If the preprocessor is successful and convetges to a soluti on of the box co nstrained variational inequality problem which satisfie s the standard regularity conditions used in [16] fo r the local convergence theory. then tk = l fo r all k sufficiently large and g(tixk) = O(I. IV'g63(.: l) l l2), so none of the tests in (36) will be satisfied.

We applied the method just described to al l test problem s fr om the MCPLIB and GAMSLIB libraries, see [4,[START_REF] Ferris | Accessing realistic mixed complementarity pro blems within MATLAB[END_REF], using all the diffe rent starting points which are available within the MA TLAB environn1ent.

We report the numerical results in Table I Looking at Tables I and11, we see that we have just a fe w fa ilures on some difficult test pro blen1s, whereas the overall behavior of our n1 ethod is q ui te good. Although many of the simple problems were solved by the preprocessor (i .e .. t here are no N-and no G-steps), the hybrid Josephy Newton tnethod introduced in this paper was neces sary in order to solve a number of other test examples.

In fa ct, we made the fo llowing observation during the testing phase for our algorithm: Both the preprocessor fr om [16] and the hybri d Newton-type method discussed in this paper try to . m inimize the D-gap In fa ct, this difficulty arises, e.g., fo r the billups example. In general , ho\v ever� our observation is that the tn ethod fr om [1 6] tends to converge It is therefore our feeling that the robustness of many existing solvers can be improved by using the search direction fr om our hybrid Josephy Newton 1nethod whenever the underlying solver does not seem to converge.

Finally, let us shortly mention what happens if we do not use our preprocessor. Obviously, this changes the numerical res ults. In general, however, the changes are not very dramatic, for example , the three gafni examples can be solved in just 3-4 1 osephy-Newton steps (instead of 12-13 linear system solves when using the preprocessor), the fo ur powell_mcp problems can be solved in 5-8 Josephy-Newton iterations (compared to 6-8 linear system solves), and the cVemge example is solved in 6 Josephy-Newton iterations (compared to 1 0 iterations in the preprocessor). The diffe rences are mor.e dramatic fo r the j osephy and koj shin test problems, mainly due to the fa ct that the linearized subproble1n s quite often do not have a solution fo r these examples. For instance, it takes 97 Josephy-Newton iterations (14 N-steps and 83 G-steps) in order to solve th e josephy problem when using the fourth starting point.

CONCLUDING REMARKS

The variational inequality problem is reformulated as an unconstrained minimization problem by using the D-gap function gn/3• Some properties of the affine Vl subproblem ofthe classical Josephy-Newton method are studied. A hybrid Josephy-Newton method is then proposed to mini mize the fu nction gtt d• Under mild conditions, the proposed n1ethod is shown to be globally convergent. If some additional assumptions are satisfied, then the sequence converges quadratically or superlinearly to a solution of the original variational inequality problem . A sufficient condition is given fo r a solution x* of the VIP to be regular. This condi tion is only concerned with the mapping F,. unlike the conditions in

  (x)z(x')II < 11.:-.:(x') 11 + 11.:-:-(x) 11 ] IIF(x) -F(x') +V' F(x1) (x'-x) 11 I" + llvF(x')T-vF(x)T[I•i[:-(x)-xjl).

  (2) yields F;(X) (zi(.Y) -x;) > 0, i = L ... , n.

( 20 )

 20 Similarly, fron1 x EX and ( 19), we have Fi(x)(x;-=t(x)) + [vF(x)T (.:(x)-x)]Jx;-.:i(x)) > 0. i-1, ... , n.

  (.:;(x)-x;)[vF(x)T(.:(x)-x)L < 0, i-I. ....n.

  ) then let ).. k := 1 and go to Step 3. If the linearized VIP (25) is not sol vable or if dk does not satisfy the condition k k ) then set d := -\7 g n/lx .

  THEORE!\1 3. 1 Su ppose that the mapp ing F is continuousfŞ dfff erenti ahle. Let E = 0 and supp ose that the algorithm generates an i11/lnite sequence {xk}. Then any accumulation point x* o . f the sequence {xk} is a stationwy point of the D-gapfunction ger3• If Fis a uniform ?-function, then by [1 5, Theorem 4.1], the level sets of D-gap function gn./3 are bounded. Since {g(l',G(xk)} is nonincreasing, the boundedness of level sets guarantees the boundedness of the generated sequence {xk} and hence the existence of at least one accumulation point of{.: /}. On the other hand, by [15, Theorem 3.1], any stationary point .V; of g03 such that VF F(x) is a P-matrix is a solution of the VIP (1 ).Therefore, if Fis a uniform ?-function, then it follows fro In Theorem 3.1 that any accumulation point of the generated sequence {xk} solves the VIP (I). Because the box constrained VIP with a uniform P-function has a unique solution, we obtain the next corollary to Theorem 3.1. COROLLARY 3.2 Suppose that F is a continuously differentiable uni form ?-function and X is a box. Then for any starting point x0 E fĕ'', the sequence {xk} generated by the algorithm converges to the unique solu tion of the VIP (1).

LEMMA 3 . 4

 34 Let x* be a solution of the VIP (1) with X being a box.

  r(x) := x-Proj[/,u](x-F(x)) denote the natural residual of the box constrained variational inequality problem. We terminate our method if for sotne iterate xk, where Et : = 10-6 and c2 : ::::: 10 11.

  ln the line search rule (28), we used w = 0.5 and b = 10-4.

(3 6 )

 6 the11 termi nate the preprocessing iteration and go to the hybrid Newton 1n ethod using the previous iterate xk -l as a starting point. The actual parameters used in ( 36) are -4 fm in = 10

  test problem in MCPLIB nmnber of variables nmnber of (finite) bounds on the va ria bles x1 starting point number of iterati ons used in the preprocessing phase number of Newton steps used in the hybrid Newton phase nurn ber of gradient steps used in the hybrid Newton phase number of fu nction evaluations value of g03(x) at the final iterate x = .:/ val ue of llr(x) ll at the final iterate x = xr.

  function g03. Now, the D-gap function tn ight have a local minimum which does not correspond to a solution of the box constrained varia tional inequality problem. In that case, we wo uld expect both algorithms to run into difficulties by converging to one of these local minima, in particular, since the search directions computed by both methods are based on some local information of the v�uiational inequality problem.

  minimum of gcr,a much more often than the method discussed here. This seems to indicate that, fr om a global point of view, the search direction con1puted by our hybrid Newton-type method is a much better search direction than the one computed by the nonsmooth Newton method in [16}.

[ 6 .

 6 1 7,24] .

  11 ,

	the inequality (12) yields

TABLE I

 I Numerical results for MCPLI B test problems

	Problem bertsekas bertsekas bertsekas bert oc	n I5 15 15 5000	m I 5 I5 15 2000	I 2 3 1	0 I l 4	4 4 4 0	0 0 0 0	11 1 2 12 6	. f) g",;i(x . 7.5e-1 5 7.5e -15 7.5c -15 2.7e -7 llr(: \lll 2.7e -7 2.7c -7 3.8e -28 6. 1e -1 4
	billups	I	1	1
	bra tu	5625 11250	1	13	0	0	29	1.0e -20 1.9c-10
	choi	13	26	1	4	0	0	5	4.4e -15	2. Ie -7
	colvdual	20	20	1
	colvdual colvnlp colvnlp cycle	20 15 1 5 1	20 15 15 1	2 ") ..:.. 1	() 1 3	3 3 0	0 0 0	10 20 5	2.0e -13 3 . 5e -1 2 2.3e -21 1.5e -1 0 J.4c -6 5.9c -6
	ehl k40 ehl_k60 ehl_k80 ehl_kost ehl_kost ehl_kost explcp freehcrt freebert freebert freebert fr eebert gafni gafni gafni hanskoop hanskoop	4 1 61 8 1 101 101 101 16 1 5 15 15 15 15 5 5 5 14 1 4	40 60 80 100 lOO 100 16 1 0 10 10 10 10 10 10 10 14 14	I 1 1 1 ") ... 3 1 1 2 3 4 5 1 2 3 1 2	I2 22 24 28 28 28 15 0 I 0 0 1 13 12 1 3 0 0	8 8 8 8 8 8 0 3 3 3 4 3 0 0 0 4 6	0 () 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0	115 22 1 233 273 273 273 31 10 I6 10 11 9 46 44 46 13 17	9.6e -23 3. 1 e-11 1.7e-1 7 1.3e-8 3.1e -17 1. 7e -8 6.6e -16 8. 1e -.. 8 6.6e -16 8. le-8 6.6e -1 6 8. 1 e --8 0 () 2.4e -I6 4.8e -8 6.6e -23 2.6e -11 2.4e -16 4.8e -8 4.4e -21 2.1e -10 3.2e -19 1.8e -9 7.5c-1 9 2.7e -9 2.7e -9 7 . 5c-1 9 7.5e 19 2.7e -9 2.3e -14 4.8c -7 2. 1e •-14 4.5e -7
	hanskoop	1 4	14	3	0	5	0	1 4	5.3c -14	7.3e -7
	hanskoop	14	14	4	4	4	0	22	3.2e -14 5.7e -7
	hanskoop hydroc06 hydroc20 jel joscphy josephy josephy josephy josephy	14 29 99 6 4 4 4 4 4	14 11 39 6 4 4 4 4 4	5 I 1 I 1 2 3 4 5	0 5 8 8 10 7 11 4 3	7 0 0 0 0 0 0 0 0	0 0 0 0 0 0 0 0 0	24 7 1 0 16 23 15 24 5 4	6.0e -15 1.1e-25 l . l e -12 2.4e -7 3.7e-1 4 6.0e-7 1.4e -15 1.2e -7 2.9c-1 7 1.7e -8 4.5e -22 6 . 7e-11 2.9e -17 1. 7e -8 2. le -7 4.2e -15 4.3e -15 2. 1e -7
	josephy kojshin kojshin kojshin kojshin kojshin	4 4 4 4 4 4	4 4 4 4 4 4	6 I 2 3 4 5	6 9 7 lO 1 3	0 0 0 0 0 0	0 0 0 () 0 0	12 22 14 23 ") ... 4	2.2e -15 l.2e -20 3.5e -10 1.5e -7 1.3e -23 1.2e -11 1.2e-20 3.5e -lO 0 0 4 . 3e -1 5 2. 1 e-7
	kojshin	4	4	6	5	0	0	7	6.3e -16 7.59e -8
	mathinum	3	3	1	'1') ,.._	0	0	47	4.le-1 4	6.4e -7
	mathinum mathinum	3 3	3 3	2 3	4 32	0 0	0 0	5 76	3.4e -16 4.3e -14	5.8e -8 6.5e -7

SP P-steps N-steps G-steps F-evaf.

TABLE I (

 I Continued)SP P-steps N-steps G-steps F-eval. gl>/l(. i) JJr(x)rl l

	Problem	n	m
	mathinum mathisum mathisum mathisum mathisum methan08 nash nash obstacle	3 4 4 4 4 3 1 10 10 2500	3 4 4 4 4 1 5 10 10 5000	4 1 2 3 4 1 1 ') I	5 4 5 26 5 4 6 1 2 10	0 0 0 0 0 0 0 0 0	0 0 0 0 0 0 0 0 0	6 6 6 53 6 5 7 39 11	2 . 7e -1 4 5. le -7 2.0e -22 4.4e -11 5.4e -1 5 2.3e -7 3.7e -14 6. 1e -7 8.2e -21 2.8e -10 6.3e -24 7.9e --12 3.8e -17 1.9e-8 6.le -15 2.5e -7 3.3e -31 2.3e -15
	opt_cont3 1 opt_cont127 4096 1024 opt_cont255 8192 opt_cont 51 1 16384 16384 1 1 024 1 4096 1 8192 1 pgvon l05 105 105 1 pgvon l06 106 106 I ptes 42 52 1 powell 16 16 1 powell 16 16 2 powell 16 16 3 powell 16 16 4 powell_mcp 8 0 1 powell_m cp 8 0 2 powell_mcp 8 0 3	5 8 t l 12 20 20 4 6 0 0 6 7 8	0 0 0 0 0 6 3 4 12 7 0 0 0	0 0 0 0 0 0 0 0 0 0 0 0 0	12 32 51 52 91 153 1 7 "" -..:. 35 21 7 8 9	3.6e -30 4.3e -15 1.4e-29 9.0e -15 2.9e -29 l.2e -14 6.0e -2 9 l.6e -14 2.7e -12 5.2e -6 2. le-16 4.6e -8 7.5e -20 8.6e -10 L8e -16 4.3e -8 2.3e -2I l.Se-10 l.7e -17 1.3e -8 4.2e -24 6. 5e -12 4.7e -25 2.2e -12 2.2e -16 4.7e -8
	powell .... mcp scarf anum scarfanum scarf anum scarfasum scarf as urn scarfasmn scarfbnum scarfbnum scarfbsum scarfbsum sppe sppe to bin to bin	8 13 13 I3 14 14 14 39 39 40 40 27 27 42 42	0 13 1 3 13 1 4 14 14 39 39 40 40 27 27 42 42	4 1 2 3 1 2 3 2 2 I 2 1 2	7 0 0 9 4 0 9 13 26 23 21 15 22	0 4 5 0 0 3 0 0 0 3 4 0 0	0 0 0 0 0 0 0 0 0 0 0 0 0	8 20 20 12 6 16 12 62 66 60 59 46 83	l.?e-1 5 L3e-7 8.2e -13 2.9e -6 l.3e-16 3.7e -8 5.2e -20 7. 2e-l0 3.2e -18 5.6e -9 8.8e -13 3.0e-6 5.2e -20 7.2e -10 5.3e -14 7.3e-7 9.6e -14 9.8e -7 1.4e -22 3.7e -11 3. 1e -21 1.7e-10 9.6e -23 3J e-11 4.0e -21 2.0e-10

TABLE 11

 11 Numerical results fo r GAMSLIB test problems

	Problem	n	m
	cafemgc cammcp cmmge co2mge dmcmge ers82mcp etamge han smcp hansmge harkmcp harmge kehomge kormcp mr5mcp nsmge	101 242 242 101 9 6 20� 208 170 170 232 () 114 114 43 43 43 43 32 32 !I 9 9 9 78 0 350 350 212 212	1 I 1 1 I 1 I I I I I I I I 1	10 6 6 2 147 5 15 33 5 26 I 12 3 I 5	0 () 0 0 33 0 0 ll 5 0 6 0 0 9 14	0 0 0 0 8 0 0 16 0 0 0 0 0 0 0	29 . 8 53 1 5 1658 6 42 787 62 60 18 18 5 51 59	3.3e -16 l.le-16 2.0e -32 4.4e-16 5.7e -8 3.3e -8 1.6e -15 1.3e -7 1.5e-14 3.9e -7 1.8e -24 4.2e -12 1 .5e -• 14 3.8e -7 2.4e -I3 1.6e-6 2.8e -13 1.7e -6 l.Oe -14 3.2e-7 1.4e -13 1.2e -6 4. 3e -23 2. 1 e-ll 1.8e -25 1.3e-12 7.3e -22 8.5e -11 2.3e-l 7 1.5e-8
	oligomcp scarfmcp scarfmoe c transmcp nvo 3mcp	6 18 18 1 1 6	6 18 IR 6	I I 1 1 I	6 0 0 0 8	0 4 6 1 0	0 0 0 0 0	9 14 18 3 16	l.Oe -20 3.2e -10 5.4e -12 7.3e -6 1.4e-1 5 1.2e -7 2.0e-1 7 1.4e -8 1.4e-15 1 .2e -7

SP P-steps N-steps G-steps F-eval. f g,3(X ) lfr(xf) ff

Acknowledgment

The first author would like to thank Prof. Y. Yuan fo r his constant help and encouragement.

The research of this author was supported by the DFG (Deutsche Forschungsgemeinschaft).

t The research of this author was supported by Project 19601035 of NSFC in China.

+ The work of this author was supported in part by the Scientific Research Grant-in-Aid from the Ministry of Education, Science and Culture, Japan.


A C Or

Hugenholtz

THE DYNAMICS OF A TIPPE TOP*

Keywords: classical rigid-body dynamics, gyroscopic motions, contact frictions, rolling, sliding AMS subject classifications. 70E05, 70J25, 70K20, 73T05

 [START_REF] Hugenholtz | On tops rising by friction[END_REF], followed by Parkyn [START_REF] Parkyn | The rising of tops with rounded pegs[END_REF], [START_REF]The inverting top[END_REF], also performed similar analyses, which confirmed the importance of eccentricity and sliding friction to the stability of the Tippe top.

As later noted by Kane and Levinson [START_REF] Kane | A realistic solution of the symmetric top problem[END_REF], however, the results of [START_REF] Braams | On the influence of friction on the motion of a top[END_REF], [START_REF]The inverting top[END_REF] are rather dependent on initial conditions-a property that does not seem to fit observations well. At one time, the Tippe top had attracted such attention that numerous articles were published [START_REF] Fokker | The tracks of tops' pegs on the floor[END_REF]- [START_REF] Johnson | The tippy top[END_REF]. To summarize, many of these studies recognized the importance of sliding friction, but yet most only attempted to solve the linear stability problem. O'Brien and Synge [START_REF] Synge | The Instability of the Tippe-Top Explained by Sliding Frictwn[END_REF] argued that the simple friction law cannot account for the observed instability. Instead, the authors used a viscous friction law, which has the viscous friction varied linearly with the sliding velocity at the contact point. The authors derived the stability criterion, which shows that the top becomes unstable when exceeding a certain spin speed. Cohen [START_REF] Cohen | The tippe top revisited[END_REF] revisited the Tippe top problem in the 1970s and provided the first numerical simulation of the reversal motions of the top. Cohen's paper also showed a photograph of Wolfgang Pauli and Niels Bohr fascinated by the Tippe top in action. Kane and Levinson [START_REF] Kane | A realistic solution of the symmetric top problem[END_REF] conceived the Tippe top as a particular case of the more general rigid body problem and developed full numerical simulations for bodies of different shapes using the general equations of motions. However, their paper dealt with the Tippe top problem only in passing mention, although some other examples such as that of a hemispherical shell were given. Kane and Levinson also criticized O'Brien and Synge's stability condition for being insensitive to the frictional coefficient. They also criticized Cohen's model as not including adequate provisions for transitions from sliding to rolling. Indeed, Kane and Levinson's model permits a switch between rolling and sliding motions. The sliding friction according to Coulomb is /J>K N, where /J>K is the coefficient of kinetic or dynamic friction and N is the normal reaction. A switch to rolling is assumed to occur as the absolute contact velocity vanishes. A switch to sliding occurs when the friction required to maintain rolling exceeds p,5N, where JJs is the coefficient of static friction. Typically, the authors assume /JK to be slightly less than P,s to permit the mechanism to work. Like Cohen, however, Kane and Levinson do not include the viscous term in their friction law. For the class of problems, the contact friction is extremely important, and it dictates the motions. Descriptions of friction laws are available in the literature since the last century. For references, refer to Routh [START_REF] Routh | Dynamics of a System of Rigid Bodies[END_REF] and the references in Cohen [START_REF] Cohen | The tippe top revisited[END_REF]. The aim of this paper is to provide a direct solution to the Tippe top motion, which may serve to summarize previous findings and to resolve the puzzles that still remain. Indeed, it becomes clear in the following that both viscous and Coulomb frictions can cause the transition to reversal. For the case of the Tippe top, a transition from sliding to rolling never seems to occur. For this purpose, our model uses a general frictional law that includes both Coulomb and the viscous terms, with sliding being assumed throughout. Along with the sliding model, the model based on the assumption of pure rolling is also simulated, to provide a comparison between the two cases. Section 2 contains the formulation of the problem. Section 3 provides the linear stability analysis. Section 4 contains the full numerical solution of the nonlinear problem. In §5, we conclude with some further remarks.

2. Formulation. Governing equ ations. We refer to Fig. 2 for the geometry of spherical top. We use a nonrotating Cartesian frame attached to the sphere's center 0, Fa, which is defined by a triple of unit vectors, (f, J, K)T (T denotes transpose of a column vector). A second Cartesian frame fixed rigidly to the body, Fb, is also used. This frame's axes are defined by a second triple of unit vectors, (T, ), k)T, where the axes align with the body's principal axes. The vectrix notation [START_REF] Hughes | Spacecraft Attitude Dynamics[END_REF] is used here for developing the governing equations from a vectorial form. In this notation, an abstract vector iJ is expressed as

iJ = f[;v = f[v,
where the upper (lower) case boldface variable refers to the column vector relative to (2)

tJ(Vo + qk) = F-tJFr� 1 K, d . 1 dt (low)= f.LQV�k-tJqFr� kxK-KxF,
where 10 is the moment of inertia tensor about point 0; V 0 is the absolute velocity at point 0; w is the angular velocity, k is the unit vector along the top axis; F is the force acting on the body through the contact point. The major nondimensional parameters are the scaled top mass f.L, the Froude number Fr, the inertia ratio a, and the initial center of mass offset q. The nondimensional parameters are defined as

(3) c a= Ao'
where M is the body mass; g is the gravity acceleration; A0 and C are, respectively, principal moments of inertia about point 0. The equilibrium solution of the initial state is F= tJFr-1 K, w= DK, V0 = 0, k= K .

We further introduce three direction cosine parameters l, m, n for the attitude, to facilitate the equations to be expressed in the component form. The direction cosines satisfy identities

k = (l, m,n) T ,
It can be easily checked that 10 has an expression

lo = b + (0'-l)kk r ,
where 13 is the order-3 identity matrix. The transverse axes have a degenerate ori entation due to axisymmetry. The kinematics are governed by ( 4)

Equations ( 1), ( 2), ( 4) complete the twelfth-order six degree-of-freedom dynamical problem with the specification of F. The translational displacement X0 is not coupled into the dynamics and can be determined from

Xo = Vo .
In the case of a pure rolling, the velocity at the moving contact point vanishes iden tically throughout the motion, that is,

Vc =0,

The rolling constraint permits F to be eliminated from combining the translational and rotational equations. In general, the force at contact point F has the following representation:

F = NK + f,
where N is the normal reaction and f is the tangential force. While N can be deter mined from the normal component of the dynamical equations, f must be supplied independently. In any case, it is known that f opposes motion. In a dimensional form, f is given by [START_REF] Parkyn | The rising of tops with rounded pegs[END_REF] where f.lv is the coefficient of viscous friction; Vc is given by V c = Vcev. The above expression is a generalization of the Coulomb friction law to include the viscous term. The Coulomb's dynamic law the corresponds to the following truncated expression:

f = -tJx N,
The next term in ( 6) is proportional to the sliding velocity, which is similar to viscous friction. In nondimensional form, the friction law becomes [START_REF]The inverting top[END_REF] where now the Coulomb parameter ax and the viscous parameter av are given by 3. The results. (i) Stability. The component form of our governing equations is shown in the Appendix. O'Brien and Synge's stability condition appears qualitatively correct, but the dependence of the stability on the friction coefficient is not included.

To understand the roles of various parameters, the stability problem will be briefly discussed before we turn to the numerical sliding model. First, we consider the two classical examples of a top with a fixed-pivot and a rolling peg. Stability conditions similar to those given below can be found in [START_REF] Crabtree | An Elementary Treatment of the Theory of Spinning Tops and Gyroscopic Motion[END_REF], [START_REF] Perry | Spinning Tops and Gyroscopic Motions[END_REF]. Introducing time dependence of ei>.t, we obtain the characteristic roots after simple manipulation of the governing equations, (

) 8
where & is the axial to transverse moments of inertia ratio about the contact point, and

for a rolling top, for a fixed-pivot top.

The stability condition is and a rolling top always satisfies this condition since q < 0. A rolling top also possesses two integrals of motion (a fixed-pivot top has three integrals). One is the total energy. The other is the angular momentum projected on to the line connecting the contact point C to the center of mass G, as there is no moment applied about this line. This angular momentum hcg is given by [START_REF] Braams | The symmetrical spherical top[END_REF] reg = qk +K.

Next, we consider the sliding case. As remarked by Kane and Levinson [START_REF] Kane | A realistic solution of the symmetric top problem[END_REF], the linear stability criterion of O'Brien and Synge [START_REF] Synge | The Instability of the Tippe-Top Explained by Sliding Frictwn[END_REF] is inadequate, since it does not depend on the friction coefficient. We only need to focus on the case of aK = 0 and av =J 0. The Coulomb term is nonlinear. The subsequent numerical simulations will show that this term alone can also produce instability, but with an algebraic rather than exponential growth rate. We let W3 = n = 1, w =W I + iw2, <I>= l + im, and v = vl + iV2 and we obtain the following linearized set:

(11)

w + ij.LqV = i(-(o--1) + j.LqFr� 1 )<I> + (-av + i(o--1))w + iavV, <i> = iii>-iw, v -iqw = q<I> -( q + i � nw -:v v,
The problem depends on the parameters J.L, q, Fr, a-, and av. The growth rate and frequency dependence on these parameters will be discussed only in cases relevant to the physical problem. Figures 3(a) and 3(b) show the change of growth rate and frequency, respectively, with respect to the center-of-mass offset q for different Fr� 1 (marked alongside curves) at a-= 1.0, J.L = 2.28, and av = 0.5. We see that a weak spin relative to gravity (i.e., small Fr) is stabilizing if the center of mass is below 0 and destabilizing if it is above 0. A reversal of this stability property occurs at high Fr. The inertia ratio a-also affects the stability property, as suggested by Figs. 3 (c) and 3 (d) (O" marked alongside curves). In the limit of very strong spin or very weak gravity, i.e., Fr-1 = 0, the inertia ratio stabilizes the top above the value of O" = 1. 14. The result is consistent with the well-established principle that the stability of a body depends on the shape of its moment of inertia ellipsoid, if energy dissipation is permitted (see also [START_REF] Braams | On the influence of friction on the motion of a top[END_REF]). An oblate spheroid (O" > 1) indicates stability, while a prolate spheroid (O" < 1) indicates instability. To show the dependence of the growth rate on the friction coefficient av , we now consider the case where Fr-1 = 0 and O" = 1. Figure 4 shows a family of growth-rate curves versus q. The solid curves rise with increasing av to 0. 1 to 1.0. For this branch, a stronger friction means more destabilization. The family of curves has a maxima with respect to av , however, after which further increase in friction leads to a decrease in growth rate. The second branch corresponds to av = 2.0 and 5.0 (dotted lines). Thus, while a moderate amount of friction destabilizes, a stronger friction tends to destabilize. finite av and later we examine the effect of finite aK. The simulation corresponds to av = 0.5, JL = 2.28, q = -0. 1, and CJ = 1.0. We present two solutions in Fig. 5, which correspond to Fr-1 = 0.3 and Fr-1 = 0. 1 (marked alongside curves). With the top's radius of 0.6366 inch, the parameter values correspond to a spin speed of 7. 15 and 12. 4 revolutions per second, respectively. For the low Fr case, the initial conditions are w 1 = 0.021g, w2 = 0.0047, w3 = -o.ggg7, l = 0.0676, m = 0.0676, n = o.gg54, V 1 = 0.0g35, and V2 = 0.0384. For the high Fr case, we have w 1 = 0.0021g, W2 = 0.0047, W3 = -o.gggg, [ = 0.0080, m = 0.0150, n = 0.9gg8, V 1 = 0.0110, and V2 = 0.0038. Figure 5(a) shows that there is a drop of spin speed after reversal due to a loss of energy and the work done in raising the center of mass. Figure 5(b) shows the magnitude of transverse velocity (note that all components are measured with respect to Fa)• Where the top axis passes through a polar angle of goo, the transverse angular speed only peaks at a small value of about 0. 12. Figure 5(c) shows the polar angle of the top's overturning axis, measured from the upward vertical. The crossing at the goo plane occurs at about 78 and 41 spin periods from the start, or 10.g seconds and 3.3 seconds after spun off. For the case of Fr-1 = 0.3, we further show the components of energy. Figure 6 indicates that the translational energy at G is very small. Even the potential energy is significantly smaller than the spin en ergy (about G) for this low Fr case. The dissipation of energy is quite noticeable.

The axial and transverse angular momenta are quite similar to the angular velocities, which will not be shown. For the same simulation, we show in Fig. 7 the translational displacement of the sphere's center of mass. Along the locus, the translational veloc ity vector is indicated by arrows at a regular interval of 5/(27r) spin-period. In the enlarged portion of the curve, the arrows are shown at intervals of 1/(27r). In Fig. 8, a three-dimensional view of motion of the axis is shown by tracking k on a unit sphere at uniform intervals of 1/(27r) (cross marks). The axis spirals outward from the north pole and inward to the south pole.

Finally, we examine the friction effect. For the same case with aK = 0, av = 0.5, and Fr-1 = 0.3, we show a comparison between the magnitude of friction force ac cording to the sliding model (solid) and to the rolling model (solid lines) in Fig. g(a).

Since O:K = 0, the sliding friction is thus proportional to the contact velocity Vc, and therefore the sliding curve can be interpreted as Vc apart from a multiplica tion constant. The friction required for maintaining rolling appears fairly stationary. However, during the tumbling process, the sliding friction typically exceeds the rolling friction. Toward the final equilibrium, a transition from sliding to rolling may occur, as the contact speed tends to zero. In One last debatable issue is whether transitions between rolling and sliding are possible during the course of the motion. In all the simulations we studied, the magnitude of the contact velocity hardly ever vanishes. As an example, we show the sliding speed for case ( 4) in Fig. 9( c). The result suggests that a switch between sliding and rolling is not an important ingredient for our analysis. friction can also destabilize. It produces an initial growth that is algebraic rather than exponential. (ii) A linear stability analysis based on the viscous friction indicates that the magnitude of the friction coefficient is important to both the growth rate and stability limit. (iii) The phenomenon of the Tippe top appears to be almost exclusively sliding. Once sliding is established, there is hardly ever any switching to rolling motion.

i � 100�• •••••••••• •• • •••• •••••• • +•• •••••••••••• 7�•••• (b) � Time (no.
Appendix. The governing equations in a scalar form. The angular mo mentum equations are [START_REF] Braams | The tippe top[END_REF] w1 +(a-1)lw3 + p,q(-nV2 + mV3) =(a-1)( -w2 + m)w3-p,qFr-1 m + h , w2 +(a-1) mw3 + p,q(nV1 -lV3) =(a-1)(wl -l)w3 + p,qFr-1 l-h , aw 3 +(a-1)(lw 1 + mw2) + p,q( -mV1 + ZV2) =(a-1)(( -w2 + m)w1 + (w1 -l)w2), where h and h are the traction force. The attitude equations are

(1 3)
The translational momentum equations are V1 + q(nw2-ffiW3) = q((w� + w5)l-W1W2m-W1W3n) + �fi , (14) V2-q(nw1 -lw3) = q((w� + wr)m-W1W2l-W2W3n) + �h , v3 + q( mwl -lw2) = q((wr + w?)n-wlw3l-W2W3m) + lN-Fr-1 ,

1-'

where N is the normal force. The velocity at the point of contact is given by [START_REF] Hart | Angular momentum and tippe top[END_REF] Since the top stays in contact with the horizontal surface, we can eliminate the ver tical translational momentum equation. Thus N is determined the vertical kinemtic constraint. In the rolling model, we can eliminate h and h from the two set of mo mentum equations. Then ( 15) is used to eliminate the translational velocity. The three remaining equations, together with the attitude equations form the rolling model. For the sliding model, ( 12)-( 14) are solved in conjunction with the following friction law:

FIG. 1 .

 1 FIG. 1. Toy top configuration and tumbling sequence.

FIG. 2 .

 2 FIG. 2. Spherical model configuration.

Fa(

  Fb), that is, v = (VI , v2, V3)T and v = ( vl ,v2,v3)T. The cross product between two vectors ilx vis defined as F;ru xv or, equivalently, F[ u xv. For a column vector u = ( u 1 , u2, u3)T, u x is the skew-symmetric matrix To reduce the number of parameters, we introduce the spin period n� I as the timescale and the radius R as the lengthscale. The linear and angular momentum equations referring to Fa are then (1)

( 4 FIG. 3 .

 43 FIG. 3. Stability plots showing (a) growth rate; (b) frequency for a= 1 and varying Fr-1; (c) growth rate; (d) fre quency for Fr-1 = 0 and varying a.

FIG. 3

 3 FIG. 3 (continued).

2 FIG. 4 .

 24 FIG. 4. Stability plot showing growth rate with varying fr iction parameter, for Fr-1 = 0 and 0" = 1.

  Fig. g(b), we show the tumbling effect for four different cases of friction, in terms of the polar angle and at Fr-1 = 0.3. The same parameters and initial conditions as before apply, and the friction coefficients used are (1) O:K = 0, O:v = 0.5; (2) O:K = 0.5, O:v = 0.5; (3) O:K = 0.5, O:v = 0. 1; (4) aK = 0.5, av = 0.0. In the last case, viscous friction is absent. The result shows that tumbling can still occur but appears considerably weaker. The destabilizing role of the Coulomb friction becomes clear if a Fourier expansion of the Coulomb term is performed. The expansion contains a sinusoidal forcing that has the same sign as the viscous friction. The Coulomb forcing yields an algebraic growth, in contrast to the viscous forcing that yields an exponential growth. In reality, it is anticipated that both Coulomb and viscous frictions are present. Thus cases (2), (3) are more relevant.

FIG. 5 .FIG. 6 .FIG. 7 . 5 XFIG. 8 .

 56758 FIG. 5. Nonlinear simulations showing the (a) spin speed; (b) transverse angular rate; (c) polar angle in time for the cases Fr-1 = 0.1 and 0.3.

4 .

 4 Concluding remarks. The complete sequence of tumbling motions of the Tippe top between the two vertical spin equilibria is simulated numerically using a spherical model with an inertia tensor close to the actual top's. Several important points concerning the effects of friction are noted. (i) The viscous friction was first introduced in[START_REF] Crabtree | An Elementary Treatment of the Theory of Spinning Tops and Gyroscopic Motion[END_REF] to produce instability of the spin equilibrium. In fact, the Coulomb

  . of spin period)

FIG. 9 .

 9 FIG. 9. Time plots to illustrate the effect of changing fr ictional coefficients: (a) sliding fr iction (solid) and fr iction according to the rolling model (dotted); (b) polar angle change; (c) sliding speed for the case of zero L>v.

TABLE 1

 1 Mass property of the toy top.

	Top Mass (oz.)		0.64
	Density of wood (oz-in3)		0.42
	Trans. MOl about 0 (oz-in2)	0.067
	Trans. MOl about G (oz-in2)	0.066
	Trans. MOl about C (oz-in2)	0.220
	Axial MOl (oz-in2)		0.066
	q, Center of mass offset	-0.094
	<J, Inertia ratio about G		1.011
	J.£, Dimensionless mass		2.28
	1/Froude (spin at 100 rpm)		5.53
	(spin at 200 rpm)		1.38
	(spin at 300 rpm)		0.62
	Length measures (inch) of the toy top.
	Stick diameter, a	0.25
	Cavity diamater, b	0.75
	Cavity depth, d	0.50
	C.M. height, h	0.57
	Radius of sphere, R 0.63
	Stick length, l	1.00

Acknowledgments. I thank the anonymous referees, particularly Dr. J .R. Ockendon, for numerous suggestions and for the source of references. Also, this work would not have been completed without the inspiration of the toy top.


Corinne Berzin
email: corinne.berzin@upmf-grenoble.fr

José R León

Estimating the Hurst parameter

Keywords: central limit theorem, estimation, fractional Brownian motion, Gaussian processes, Hermite polynomials

Let {b H (t), t ∈ R} be a fractional Brownian motion with parameter 0 < H < 1. We are interested in the estimation of this parameter. To achieve this goal, we consider certain functionals of the second order increments of b H (•), using variation technics. Based on an almost-sure convergence theorem for general functionals, we single out particular functionals that allows to construct certain regression models for the parameter H. We show that this regression based estimator for H is asymptotically unbiased, consistent and that it satisfies a Central Limit Theorem (CLT).

Introduction

In this work, our main interest is to provide estimators of the Hurst parameter H (0 < H < 1) for a fractional Brownian motion b H (•). Instead of the original process we assume we observe a smoothed by convolution process, defined as b ε H (u) = ϕ ε * b H (u). Here ε → 0 is the smoothing parameter and ϕ ε is a positive convolution kernel. Our method consists in obtaining some least squares estimators in certain regression models. The asymptotic behavior of such estimators will be equivalent to those of certain non-linear functionals of the Gaussian process b ε H (•), whose asymptotic behavior can be obtained by using the method of moments, via the Diagram formula.

The estimation of the Hurst parameter, or of the other parameters associated to long rang dependence, has been recently a matter of intense research. For this different approaches have been developed. Some of them consider spectral methods: as the Whittle estimator [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent Gaussian stationary time series[END_REF], [START_REF] Velazco | Whittle pseudo-maximun likelihood estimation of nonstationary time series[END_REF], and also estimators based on the log-periodogram of the data [START_REF] Ludeña | Parametric estimation for Gaussian long-range dependent processes based on the log-periodogram[END_REF], [START_REF] Moulines | Broad band log-periodogram regression of time series with long range dependence[END_REF]. Others use variation methods as in [START_REF] Benassi | Identification of Filtered White Noises[END_REF], [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF] and [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a gaussian process[END_REF], or are based on wavelet techniques [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF].

Here the problem is different because we use a mollified version of the original process. Nevertheless, our method can be related with the variation methods. It is important to mention that the recent works of Coeurjolly [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF] and Feuerverger et al. [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], constitute our main sources of inspiration.

Let us introduce the problem. Let {b H (t), t ∈ R} be the fractional Brownian motion with parameter 0 < H < 1. We study the second order increments of b H (•) and we shall establish in section 3.1.1 an almost-sure convergence in law result stated below in Corollary 3.1: For all x ∈ R λ 0 u 1 :

b

H (u + 2ε) -2b H (u + ε) + b H (u) ε H v 2H √ 4 -2 2H x a.s.
-→ Pr{N * x}, where v 2 2H = var[b H [START_REF] Azaïs | Almost sure oscillation of certain random processes[END_REF]], λ is the Lebesgue measure and N * is a standard Gaussian random variable.

This result also implies that for smooth functions f with at most polynomial growth, we have :

1 0 f b H (u + 2ε) -2b H (u + ε) + b H (u) ε H v 2H √ 4 -2 2H du a.s. -→ E[f (N * )]. (1)
Defining

ϕ = 1 [-1,0] * 1 [-1,0] and b ε H (u) = ϕ ε * b H (u), where ϕ ε (•) = 1 ε ϕ( • ε ), we have ε 2 bε H (u) = b H (u + 2ε) -2b H (u + ε) + b H (u)
. Hence (1) can be written

1 0 f ε (2-H) bε H (u) v 2H √ 4 -2 2H du a.s. -→ E[f (N * )].
By considering other regularizations:

b ε H (u) = ϕ ε * b H (u)
, ϕ being a C 2 positive kernel with L 1 norm equal to one, we can write a result similar to (1) in the form:

1 0 f (Z ε (u)) du a.s. -→ E[f (N * )], (2)
where

Z ε (u) = ε (2-H) bε H (u) σ 2H , with σ 2 2H = var[ε (2-H) bε H (u)]
. By observing b ε H (u) at several scales of the parameter ε, i.e. h i = εc i , c i > 0, i = 1, . . . , l, we want to estimate the H parameter. To construct the estimator, we consider f β (x) = |x| β . Then using (2), we get:

log(M β (ε)) = (H -2) β log(ε) + β log(σ 2H ) + log(E[|N * | β ]) + o a.s. (1), where M β (ε) = 1 0 | bε H (u)| β du.
The following regression model can be written, for each scale h i :

Y i = aX i + b + ξ i , i = 1, . . . , l
where a = (H -2) β and for i = 1, . . . , l,

Y i = log(M β (h i )), X i = log(h i ) and b = β log(σ 2H ) + log(E[|N * | β ]).
Hence, the least squares estimator Ĥβ of H is defined as

β ( Ĥβ -2) = l i=1 z i log(M β (h i )),
where,

z i = y i l i=1 y 2 i , y i = log(c i ) - 1 l l i=1 log(c i ).
We shall prove that estimator Ĥβ is strongly consistent, asymptotically unbiased and asymptotically Gaussian. For the later, we need to study the convergence rate in (2) for a function f ∈ L 2 (φ(x) dx) where φ stands for the standard Gaussian measure. We can formulate the problem in the following way:

Let g be a function in L 2 (φ(x) dx) and g(x) = ∞ n=1 ĝn H n (x) its Hermite expansion. We have g 2 2,φ = ∞ n=1 ĝ2 n n! < +∞ and ĝ0 = 0. We prove in Theorem 3.2 that the functional defined for c > 0 by

S g (εc) = 1 √ εc 1 0 g(Z εc (u)) du,
converges in distribution towards a cylindrical centered Gaussian process X(c).

The CLT result for Ĥβ , follows by considering the asymptotic behavior of the functional associated with

g β (x) = |x| β /E[|N * | β ]-1, and from the fact that |( Ĥβ -H)/ √ ε- 1 β l i=1 z i √ c i S g β (εc i )| = o p (1)
. This theorem also allows us to consider another least squares estimator of H, weakly consistent and unbiased, defined by

Ĥlog -2 = l i=1 z i M log (h i ), where M log (ε) = 1 0 log(| bε H (u)|) du.
In this case the associated function will be

g log (x) = log(|x|) -E[log(|N * |)].
It is important to mention that working with the second order increments allows us to obtain an asymptotic Gaussian result, for all the possible values of parameter H i.e. (0, 1). In previous articles, we considered only the first order increments. This leads to two different results according to whether H belongs to (0, 3/4) or to (3/4, 1). In the first case, the limit is Gaussian and in the second one, it belongs to the second Itô-Wiener Chaos. This dichotomy makes the asymptotic behavior more involved.

The paper is organized as follows. In Section 2, we state some notation and the hypothesis under which we work. Section 3 is devoted to establishing the main results. Section 4 contains the proofs.

Hypothesis and notation

Let {b H (t), t ∈ R} be a fractional Brownian motion with parameter 0 < H < 1 (see for instance [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF]), i.e. b H (•) is a centered Gaussian process with the covariance function:

E[b H (t)b H (s)] = 1 2 v 2 2H |t| 2H + |s| 2H -|t -s| 2H with v 2 2H = [Γ(2H + 1) sin(πH)] -1 .
We define, for a C 2 density ϕ with compact support included in [-1, 1], for each t 0 and ε > 0 the regularized processes:

b ε H (t) = 1 ε ∞ -∞ ϕ t -x ε b H (x) dx and Z ε (t) = ε (2-H) bε H (t) σ 2H , with σ 2 2H := var[ε (2-H) bε H (t)] = 1 2π +∞ -∞ |x| (3-2H) | φ(-x)| 2 dx.
We can also consider the function ϕ

= 1 [-1,0] * 1 [-1,0]
. Even if it does not belong to the class defined above it is important because the second derivative of the corresponding smoothed process b ε H (t) gives place to the second order increments for fBm. This assertion is easy to see by using the following representation:

b ε H (t) = -1 -2 (u + 2) b H (t -εu) du - 0 -1 u b H (t -εu) du.
Moreover, all the results that we shall obtain below could be shown, in a similar form, for this particularly important regularization.

We shall use the Hermite polynomials, which can be defined by exp(tx

-t 2 /2) = ∞ n=0 H n (x)t n /n!.
They form an orthogonal system for the standard Gaussian measure φ(x) dx and, if h

∈ L 2 (φ(x) dx), h(x) = ∞ n=0 ĥn H n (x) and ||h|| 2 2,φ = ∞ n=0 ĥ2
n n!. Mehler's formula (see [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF]) gives a simple form to compute the covariance between two L 2 functions of Gaussian random variables:

if k ∈ L 2 (φ(x) dx), k(x) = ∞ n=0 kn H n (x) and if (X, Y ) is a standard Gaussian random vector with correlation ρ then E[h(X)k(Y )] = ∞ n=0 ĥn kn n!ρ n . (3)
Throughout the paper, C shall stand for a generic constant, whose value may change during a proof, and log for the Naperian logarithm. N * will denote a standard Gaussian random variable, and the symbol ⇒ will mean weak convergence of measures.

Results

In this section we shall establish the main results of this work. For this we need to introduce some definitions.

Let g be a function in

L 2 (φ(x) dx) such that g(x) = ∞ n=1 ĝn H n (x), with ||g|| 2 2,φ = ∞ n=1 ĝ2 n n! < +∞. We define, for x ∈ R, ρ H (x) = E[Z ε (εx + u)Z ε (u)] = -v 2 2H 2σ 2 2H ∞ -∞ φ * φ(y)|x -y| 2H dy = 1 2πσ 2 2H ∞ -∞ |y| (3-2H) e ixy | ϕ(-y)| 2 dy,
and φ(y) = φ(-y). The fact that ρ H does not depend on the variable u follows from ϕ ∈ C 2 and the stationary increments of b H .

We shall write

σ 2 g = ∞ n=1 ĝ2 n n! +∞ -∞ ρ n H (x) dx.
More generally, for x ∈ R and b, c > 0, we define

ρ H (x, b, c) = E[Z εb (εx + u)Z εc (u)] = ( -v 2 2H 2σ 2 2H )( b c ) (H+1) ∞ -∞ ∞ -∞ φ(z) φ b c (z -y) | x b -y| 2H dz dy = (bc) (2-H) 2πσ 2 2H ∞ -∞ |y| (3-2H) e ixy ϕ(-by) ϕ(cy) dy = ρ H (-x, c, b),
and

ρ g (b, c) = 1 √ bc ∞ n=1 ĝ2 n n! +∞ -∞ ρ n H (x, b, c) dx. Note that ρ H (x, c, c) = ρ H ( x c ) and then ρ g (c, c) = σ 2 g . For all m ∈ N * , for all c 1 > 0, c 2 > 0, . . . , c m > 0 and for all d 1 , d 2 ,. . . , d m ∈ R, we will denote σ 2 g,m (c, d) = m i=1 m j=1 d i d j ρ g (c i , c j ).
For ε > 0, define

S g (ε) = ε -1 2 1 0 g(Z ε (u)) du. 3.1 Convergence for (Z ε (u)) u∈[0,1] and (S g (εc)) c>0 3.1.1 Almost sure convergence Theorem 3.1 For all 0 < H < 1 and k ∈ N * , 1 0 [Z ε (u)] k du a.s. -→ ε→0 E[N * ] k . Corollary 3.1 Almost surely, for all 0 < H < 1, Z ε (•) ⇒ ε→0 N * .
The above convergence is in law, the random variable Z ε (•) is considered as a variable on ([0, 1], λ) where λ is the Lebesgue measure.

Remark : For the first order increments, a similar theorem was proved by Azaïs & Wschebor in [START_REF] Azaïs | Almost sure oscillation of certain random processes[END_REF].

Convergence in law

Theorem 3.2 For all 0 < H < 1,

S g (ε•) ⇒ ε→0 X(•), where X(•) is a cylindrical centered Gaussian process with covariance ρ g (b, c) = E[X(b)X(c)], b, c > 0.
The above convergence is in the sense of finite-dimensional distributions.

Remark 1 : If c > 0 is fixed, S g (εc) ⇒ ε→0 σ g N * .
This particular result was shown in Corollary 3.2 (i) in C. Berzin and J.R. León [START_REF] Berzin | Convergence in fractional models and applications[END_REF]. In section 4.2.1 we shall prove the following lemma

Lemma 3.1 For all m ∈ N * , for all c 1 > 0, • • •, c m > 0 and d 1 , • • •, d m ∈ R, σ 2 g,m (c, d) = lim ε→0 E[ m i=1 d i S g (εc i )] 2 0.
Remark 2: The above lemma allows us to conclude that ρ g (b, c) is a covariance function, since it is a symmetrical function in b and c.

Estimation of the Hurst parameter

Let β > 0 and

M β (ε) = 1 0 | bε H (u)| β du. ( 4
)
Thanks to Theorem 3.1 and Corollary 3.1,

1 0 |Z ε (u)| β du = ε (2-H)β σ β 2H M β (ε) a.s. -→ ε→0 E[|N * | β ].
Then,

(2 -H) β log(ε) -β log(σ 2H ) + log(M β (ε)) a.s. -→ ε→0 log(E[|N * | β ]). Thus log(M β (ε)) = (H -2) β log(ε) + β log(σ 2H ) + log(E[|N * | β ]) + o a.s. (1). ( 5
)
Let h i = εc i , c i > 0, i = 1, . . . , l. The model can be written as:

Y i = aX i + b + ξ i , i = 1, . . . , l
where a = (H -2) β and for i = 1, . . . , l,

Y i = log(M β (h i )), X i = log(h i ) and b = β log(σ 2H ) + log(E[|N * | β ]).
Hence, the least squares estimator Ĥβ of H is given by

β ( Ĥβ -2) = l i=1 z i log(M β (h i )), (6)
where

for i = 1, . . . , l, z i = y i l i=1 y 2 i and y i = log(c i ) - 1 l l i=1 log(c i ).
Note the following property l i=1 y i = 0 and

l i=1 z i y i = 1. (7)
We have the following corollary of Theorem 3.2.

Corollary 3.2 1)

Ĥβ is an asymptotically unbiased strongly consistent estimator of H.

2) Furthermore, 1 √ ε ( Ĥβ -H) ⇒ ε→0 N (0, σ 2 g β ,l (c, √ c (z/β)),
where

g β (x) = |x| β E[|N * | β ] -1 = ∞ n=1 ĝ2n,β H 2n (x), (8)
with

ĝ2n,β = n p=0 (-1) (n-p) (2p)!(n -p)!2 (n-2p) Γ(p + ( β+1 2 )) Γ( β+1 2 ) = 1 (2n)! n-1 i=0 (β -2i). ( 9
)
Remark 1: As in [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF], the variance σ 2 g β ,l (c, √ c (z/β)) is minimal for β = 2. This fact will be shown in section 4.3.1 after the proof of Corollary 3.2.

Remark 2 : Note that for β = 1, M 1 (ε) = 1 0 | bε H (u)| du = +∞ -∞ N ḃε H [0,1] (x) dx,
where

N ḃε H [0,1] (x)
denotes the number of times the process ḃε H (•) crosses level x before time 1. This remark allows to establish a link between our estimator and that defined in [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], but in this work the authors use the crossings of b ε H (•) instead of those of ḃε H (•). Now, let us define

M log (ε) = 1 0 log(| bε H (u)|) du. ( 10
)
Lemma 3.1 following Theorem 3.2, also entails that

1 0 log(|Z ε (u)|) du P -→ ε→0 E[log(|N * |)],
and also in L 2 (Ω). We proceed as before and the least squares estimator Ĥlog of H is given by

Ĥlog -2 = l i=1 z i M log (h i ). ( 11
)
Theorem 3.2 gives us the following corollary.

Corollary 3.3 1)
Ĥlog is an unbiased weakly consistent estimator of H.

2) Furthermore,

1 √ ε ( Ĥlog -H) ⇒ ε→0 N (0, σ 2 g log ,l (c, √ c z))
where

g log (x) = log(|x|) -E[log(|N * |)] = ∞ n=1 ĝ2n,log H 2n (x), (12)
with

ĝ2n,log = (-1) (n-1) 2n(2n)!! . ( 13
)
We can link the two estimators Ĥβ and Ĥlog . For this, let β(ε) a sequence of positive numbers going to zero when ε goes to zero and let Ĥβ(ε) be the corresponding estimator say

Ĥβ(ε) -2 = l i=1 z i log(M β(h i ) (h i )) β(h i )
where

M β(ε) (ε) = 1 0 | bε H (u)| β(ε) du. ( 14
)
We have the following corollary.

Corollary 3.4 If β(ε) = o( √ ε) then Ĥβ(ε)
is an asymptotically unbiased weakly consistent estimator of H and the asymptotic behavior of

1 √ ε ( Ĥβ(ε) -H) is the same one as that of 1 √ ε ( Ĥlog -H). 4 Proofs 4.1 Almost sure convergence for (Z ε (u)) u∈[0,1] Proof of Theorem 3.1. For all k ∈ N * , E[ 1 0 [Z ε (u)] k du] = E[N * ] k . Furthermore, var[ 1 0 [Z ε (u)] k du] = E[ 1 0 [Z ε (u)] k -E[N * ] k du] 2 . Let g (k) (x) = x k -E[N * ] k = ∞ n=1 ĝn,(k) H n (x). Mehler's formula (3) entails that var[ 1 0 [Z ε (u)] k du] = ∞ n=1 ĝ2 n,(k) n! 1 0 1 0 ρ n H u -v ε du dv = ∞ n=1 ĝ2 n,(k) n! 1 -1 (1 -|u|)ρ n H u ε du.
If we let u = εx, we get var[

1 0 [Z ε (u)] k du] = ε ∞ n=1 ĝ2 n,(k) n! 1 ε -1 ε (1 -ε|x|)ρ n H (x) dx. Since |ρ H (x)| 1, |ρ n H (x)| |ρ H (x)| for n 1, we have var[ 1 0 [Z ε (u)] k du] ε ∞ n=1 ĝ2 n,(k) n! 1 ε -1 ε |ρ H (x)| dx ε||g (k) || 2 2,φ 1 ε -1 ε |ρ H (x)| dx.
A Taylor expansion of order four of the function (1 -x) 2H and integration by parts

show that ρ H (x) is equivalent to -v 2 2H /σ 2 2H |x| (2H-4) 2H(2H -1)(H -1)(2H -3) for large values of |x|. Thus |ρ H (x)| is bounded from above by C |x| (2H-4) . As ||g (k) || 2 2,φ < +∞, we get var[ 1 0 [Z ε (u)] k du] C ε. Now let ε ν = ν -a , a > 1.
Using the Borel-Cantelli Lemma, one has

1 0 [Z εν (u)] k du a.s. -→ ν→∞ E[N * ] k . ( 15
)
For ε ν+1 ε ε ν let us consider the modulus of the difference

1 0 [Z εν (u)] k du - 1 0 [Z ε (u)] k du .
One has

1 0 [Z εν (u)] k du - 1 0 [Z ε (u)] k du J 1 + J 2 ,
where

J 1 = 1 - ε ε ν (2-H)k 1 0 [Z εν (u)] k du ,
and

J 2 = ε ε ν (2-H)k 1 0 [Z εν (u)] k du - 1 0 [Z ε (u)] k du .
We study J 1 and J 2 separately. Using [START_REF] Velazco | Whittle pseudo-maximun likelihood estimation of nonstationary time series[END_REF], one has sup

ε ν+1 ε εν J 1 1 - ε ν+1 ε ν (2-H)k 1 0 [Z εν (u)] k du a.s.
-→ ν→∞ 0.

Moreover, for J 2 we have

J 2 1 σ 2H k ε (2-H)k 1 0 [ bεν H (u)] k -[ bε H (u)] k du. Now bε H (u) k -bεν H (u) k = bε H (u) -bεν H (u) k-1 j=0 [ bε H (u)] j [ bεν H (u)] (k-1-j) .
Using that the trajectories of b H (•) are (H -δ)-Hölder continuous, for any δ > 0 i.e.

|b H (u + ε) -b H (u)| C(ω)ε (H-δ) , we get | bε H (u)| = | 1 ε 2 ∞ -∞ φ(v)(b H (u -εv) -b H (u))dv| Cε (H-2-δ) ,
and by similar computations bε

H (u) -bεν H (u) C ε (H-δ) ν ε -2 1 - ε ν+1 ε ν (H-δ)
.

Thus

[ bε

H (u)] k -[ bεν H (u)] k C ε (H-δ) ν ε -2 1 - ε ν+1 ε ν (H-δ) ε (H-2-δ)(k-1) ,
and then sup

ε ν+1 ε εν J 2 C 1 σ 2H k ε -δk ν+1 1 - ε ν ε ν+1 (H-δ) a.s.
-→ ν→∞ 0 as soon as δ < H (ak+1) . Thus we have proved that sup

ε ν+1 ε εν 1 0 [Z εν (u)] k du - 1 0 [Z ε (u)] k du ,
tends almost surely to zero when ν goes to infinity and Theorem 3.1 follows. Corollary 3.1 follows readily from Theorem 3.1 by the method of moments.

Convergence in law for (S g (εc)) c>0

Asymptotic variance of

m i=1 d i S g (εc i ) Proof of Lemma 3.1. Fix m ∈ N * , c 1 > 0, . . . , c m > 0 and d 1 , . . . , d m ∈ R E[ m i=1 d i S g (εc i )] 2 = m i=1 m j=1 d i d j E[S g (εc i )S g (εc j )].
For fixed b > 0 and c > 0 and by Mehler's formula (3), we get as in Section 4.1 4) , for |x| large enough. The dominated convergence theorem entails that

E[S g (εb)S g (εc)] = 1 ε √ bc ∞ n=1 ĝ2 n n! 1 0 1 0 ρ H u -v ε , b, c n du dv = 1 √ bc ∞ n=1 ĝ2 n n! 1 ε -1 ε (1 -ε|x|)[ρ H (x, b, c)] n dx. Since ρ H (x, b, c) is a correlation we have |ρ H (x, b, c)| n |ρ H (x, b, c)|. Moreover, as for ρ H (x) in section 4.1 it can be seen that ρ H (x, b, c) is equivalent to -v 2 2H /σ 2 2H |x| (2H-4) (bc) (2-H) 2H(2H- 1)(H -1)(2H -3) for |x| large enough. Thus |ρ H (x, b, c)| is bounded from above by C |x| (2H-
E[S g (εb)S g (εc)] -→ ε→0 ρ g (b, c) thus E[ m i=1 d i S g (εc i )] 2 -→ ε→0 σ 2 g,m (c, d),
this yields Remark 2 of Theorem 3.2. We want to prove that S g (εc)

⇒ ε→0 N (0; σ 2 g,m (c, d)).
Let

S g M (εc) = m i=1 d i S g M (εc i ) with g M (x) = M n=1 ĝn H n (x).
First, let us prove the following lemma.

Lemma 4.1 S g M (εc) ⇒ ε→0 N (0; σ 2 g M ,m (c, d)).
Proof of Lemma 4.1. Let fixed p ∈ N * , we want to study the asymptotic behavior of E[S g M (εc)] p . For p = 1, by Mehler's formula (see [START_REF] Berzin | Convergence in fractional models and applications[END_REF]), E[S g M (εc)] = 0. Now, for p ∈ N * -{1},

E[S g M (εc)] p = 1 ε p/2 m j 1 ,...,jp=1 d j 1 √ c j 1 • • • d jp √ c jp M l 1 ,...,lp=1 ĝl 1 • • • ĝlp × 1 0 1 0 • • • 1 0 E[H l 1 (Z εc j 1 (s 1 )) • • • H lp (Z εc jp (s p ))] ds,
where ds = ds 1 . . . ds p .

To get the asymptotic behavior of the later expression we shall use the diagram formula. Therefore, we must introduce some definitions given in [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF]. Besides that our problem is in continuous time and non-ergodic (the time scale does not tend to infinity), we can say that our theorem is a variant of the main result in that work. Nevertheless, we have decided to include the proof for completeness and to make the reading more easy .

An undirected graph G with l 1 +l 2 +. . .+l p vertices is a diagram of order (l 1 , . . . , l p ) if:

(i) The set of vertices V of the graph G has the form: V = p j=1 L j where L j = {(j, l) : 1 l l j } for j = 1, . . . , p.

(ii) Each vertex is of degree 1.

(iii) Edges may pass only between different levels.

We denote Γ = Γ(l 1 , . . . , l p ) the set of diagrams having these properties, G(V ) denotes the set of edges of G; the edges w are oriented, beginning in d 1 (w) and finishing in d 2 (w).

The diagram formula [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF] allows to write:

E[H l 1 (Z εc j 1 (s 1 )) • • • H lp (Z εc jp (s p ))] = G∈Γ w∈G(V ) d 1 (w)<d 2 (w) ρ H s d 1 (w) -s d 2 (w) ε , c j d 1 (w) , c j d 2 (w)
where G is an undirected graph with l 1 + • • • + l p vertices and p levels. A diagram is said "regular" ( [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF], p.432) if its levels can be matched in such a way that no edge passes between levels in different pairs, otherwise, it is said "irregular".

We shall classify the diagrams of Γ as in [START_REF] Ho | Limiting distributions of nonlinear vector functions of stationary Gaussian processes[END_REF], p. 1166, calling R the set of the regular graphs and R c the rest.

We start by considering R. Consider a regular diagram G * and let i be the permutation such that (i(1), i(2)), . . . , (i(p -1), i(p))

defines the diagram with p = 2q. The contribution of this diagram is

I m,p (ε)(G * ) := m j 1 ,...,jp=1 p i=1 d j i √ c j i × q k=1 ĝ2 e(k) 1 ε [0,1] 2 ρ H u -v ε , c j i(2k-1) , c j i(2k) e(k)
du dv , where e(k) is the number of edges linking i(2k -1) to i(2k).

As in Section 4.2, we get

I m,p (ε)(G * ) = m j 1 ,...,jp=1 p i=1 d j i √ c j i × q k=1 ĝ2 e(k) 1 ε -1 ε (1 -ε|x|)[ρ H (x, c j i(2k-1) , c j i(2k) )] e(k) dx . Since ρ H (x, b, c) is bounded from above by C |x| 2H-4 , for large enough |x| , ||g M || 2 2,φ < +∞ and |ρ H (x, b, c)| 1 I m,p (ε)(G * ) -→ ε→0 I m,p (G * ) := m j 1 ,...,jp=1 p i=1 d j i √ c j i × q k=1 ĝ2 e(k) ∞ -∞ [ρ H (x, c j i(2k-1) , c j i(2k) )] e(k) dx .
The right hand term can be written

I m,p (G * ) = q k=1   ĝ2 e(k) m j i(2k) =1 m j i(2k-1)=1 d j i(2k) √ c j i(2k) d j i(2k-1) √ c j i(2k-1) ∞ -∞ [ρ H (x, c j i(2k-1) , c j i(2k) )] e(k) dx   .
Since we show in Lemma 4.2 that the contribution of the irregular diagrams tends to zero, we obtain that lim

ε→0 [S gM (εc)] p = I m,p (G * ),
where the summation in goes over the regular diagrams with p levels. By computing the numbers of such diagrams, as in Breuer & Major ([4], p. 434), we get

E[S g M (εc)] p -→ ε→0 (p)!! σ p g M ,m (c, d).
This achieves the proof of Lemma 4.1. 2 The following lemma studies the contribution of the irregular graphs. Lemma 4.2 For fixed j 1 , . . . , j p we have

L(ε) := [0,1] p w∈G(V ) ρ H s d 1 (w) -s d 2 (w) ε , c j d 1 (w) , c j d 2 (w) ds = o(ε p/2 ).
Proof of Lemma 4.2.

L(ε) = [0,1] p p i=1   w∈G(V ) d 1 (w)=i ρ H s i -s d 2 (w) ε , c j i , c j d 2 (w)   ds i . Let K G (i) be the number of edges such that d 1 (w) = i, L(ε)≤ [0,1] p p i=1   1 K G (i) w∈G(V ) d 1 (w)=i ρ H s i -s d 2 (w) ε , c j i , c j d 2 (w) K G (i)   ds i .
Let A = {(c j 1 , c j 2 ), j 1 = 1, . . . , m, j 2 = 1, . . . , m}.

Using the same techniques as in the proposition of [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF] p. 435,

L(ε) ≤ p i=1 sup (c,d)∈A sup v∈[0,1] 1 0 ρ H s i -v ε , c, d K G (i) ds i ≤ p i=1 sup (c,d)∈A 1 -1 ρ H u ε , c, d K G (i)
du .

Let i be fixed in {1, . . . , p} and define g(i) = K G (i) l i , with c and d held fixed in A. We have two cases to consider:

i) K G (i) ≥ 1 1 -1 ρ H u ε , c, d K G (i) du = ε 1 ε -1 ε ρ H (x, c, d) K G (i) dx ≤ ε ∞ -∞ ρ H (x, c, d) dx ≤ C ε = C ε 1-g(i) ε g(i) ≤ C ε g(i) , ii) K G (i) = 0 1 -1 ρ H u ε , c, d K G (i) du = 2 = 2 ε 0 = 2 ε g(i) . Thus L(ε) ≤ C ε p i=1 g(i)
. As in [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF] p. 436, either p i=1 g(i) > p/2 (and in this case L(ε) = o(ε p/2 )) or there exists 1 ≤ i 0 ≤ p such that 0 < K G (i 0 ) < l i 0 , i.e. 0 < g(i 0 ) < 1, then

L(ε) ≤ C ε p i=1 g(i) ε (1-g(i 0 )) ,
and since p i=1 g(i) -p/2 ≥ 0, Lemma 4.2 follows. 2

Hence, we have proved that 1) By using ( 6), ( 4) and ( 5) we get

S g M (εc) ⇒ ε→0 N (0; σ 2 g M ,m (c, d)).
β ( Ĥβ -2) = l i=1 z i (H -2) β log(εc i ) + β log(σ 2H ) + log(E[|N * | β ]) + o a.s.,i (1)
, and property [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF] gives

β ( Ĥβ -2) = (H -2) β + o a.s. (1).
We proved that Ĥβ is a strongly consistent estimator of H. Let us see now that Ĥβ is an asymptotically unbiased estimator of H. By ( 6)

β (E[ Ĥβ ] -2) = l i=1 z i E[log(M β (εc i ))],
where M β (ε) is defined in (4). Since

1 0 |Z ε (u)| β du = ε (2-H)β σ β 2H M β (ε), (16)
by property [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], one has

β (E[ Ĥβ ] -2) = l i=1 z i E[log 1 0 |Z εc i (u)| β du ] + (H -2)β. (17)
Hence, it is enough to prove that E[log

1 0 |Z ε (u)| β du ] -→ ε→0 log E[|N * |] β . ( 18
)
For this, let us notice that since log is a concave function and log(x) x when x 0, we have

1 0 log(|Z ε (u)| β ) du log 1 0 |Z ε (u)| β du 1 0 |Z ε (u)| β du.
Thus, if we denote

X ε := 1 0 |Z ε (u)| β du + β 1 0 log(|Z ε (u)|) du ,
we have shown that log

1 0 |Z ε (u)| β du X ε . (19)
Now, since |N * | β and log(|N * |) are in L 1 (Ω), the same result is true for X ε , say

X ε ∈ L 1 (Ω). ( 20
)
Furthermore by using Lemma 3.1 following Theorem 3.2, it is easy to see that

X ε L 1 -→ ε→0 E[|N * | β ] + β E[log(|N * |)] . (21)
Finally, by using Theorem 3.1, we get:

log 1 0 |Z ε (u)| β du a.s. -→ ε→0 log (E[|N * |] β ) . (22)
Hence, ( 19), ( 20), ( 21) and ( 22) yield (18).

2) Formula ( 16) entails that

E[M β (ε)] = σ β 2H ε (H-2)β E[|N * | β ].
As in [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], let us define

A β (ε) = M β (ε) -E[M β (ε)] E[M β (ε)] .
With this definition and using a Taylor expansion for the logarithm function one has log

(M β (ε)) = log(E[M β (ε)]) + log(1 + A β (ε)) = (H -2) β log(ε) + β log(σ 2H ) + log(E[|N * | β ]) + A β (ε) + A 2 β (ε) [-1 2 + ε 1 (A β (ε))]. (23)
Let us see that

A 2 β (ε) - 1 2 + ε 1 (A β (ε)) = o p ( √ ε). (24)
By the definition of g β (see [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]), one has

A β (ε) = √ εS g β (ε),
and by Lemma 3.1 following Theorem 3.2

E[ 1 √ ε A 2 β (ε)] = √ εE[S 2 g β (ε)] = O( √ ε), so 1 √ ε A 2 β (ε) = o p (1)
and then (24) is proved. By using (23) and ( 24) we obtain log

(M β (ε)) = (H -2) β log(ε) + β log(σ 2H ) + log(E[|N * | β ]) + √ ε S g β (ε) + o p ( √ ε).
Thus, by using (6) and property [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], we have

β ( Ĥβ -2) = (H -2) β + l i=1 z i √ εc i S g β (εc i ) + o p ( √ ε). Thus ( Ĥβ -H) √ ε = 1 β l i=1 z i √ c i S g β (εc i ) + o p (1
(x) = |x| β -E[|N * | β ]. Since h β is even, we have: h β (x) = ∞ n=1 ĥ2n,β H 2n (x), with for n 1, ĥ2n,β = 2 (2n)! +∞ 0 H 2n (x)x β φ(x) dx.
By the definition of the Hermite polynomials and since for a 0 Ψ(a) := +∞ 0 x a φ(x) dx = π -1/2 2 (a-2)/2 Γ((a + 1)/2),

we get

ĥ2n,β = π -1/2 n p=0 (-1) (n-p) 2 (n-p) (n -p)!(2p)! 2 (p+β/2) Γ(p + (β + 1)/2).
To obtain the first equality in [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent Gaussian stationary time series[END_REF], it suffices to see that ĝ2n,β = ĥ2n,β 2Ψ(β) and (25) gives the required result. Moreover, Coeurjoly in [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF], has shown that ĝ2n

,β = 1 (2n)! n-1 i=0 (β - 2i). 2
Proof of Remark 1 in Corollary 3.2. Let us note that ĝ2,β = β 2 . Then for b, c > 0

ρ g β (b, c) = β 2 4 ρ g 2 (b, c) + ρ g β (b, c) where g β (x) = ∞ n=2 ĝ2n,β H 2n (x) which belongs to L 2 (φ(x) dx). Then σ 2 g β ,l (c, √ c (z/β)) = σ 2 g 2 ,l (c, √ c (z/2)) + σ 2 g β ,l (c, √ c (z/β)) σ 2 g 2 ,l (c, √ c (z/2))
since the last term in the above equality is positive by Lemma 3.1 following Theorem 3.2.

Indeed, from the definition of the function g log (see [START_REF] Ludeña | Parametric estimation for Gaussian long-range dependent processes based on the log-periodogram[END_REF]), one has We have proved that Ĥlog is a weakly consistent estimator of H. Let us show now that Ĥlog is unbiased. Since

and by [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a gaussian process[END_REF] and property [START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], we get

E[ Ĥlog ] -2 = l i=1 z i E[M log (εc i )] = H -2, ( 31
)
and Ĥlog is an unbiased estimator of H. If we divide both sides, in the above inequality, by β(ε) and using (36), we get (37). 2

Final remark: We have considered the estimation of the Hurst parameter for fractional Brownian motion. The same procedure can be implemented, with little changes, when the observed process is a Gaussian process with stationary increments. In this case the incremental variance must be: E(X t+h -X t ) 2 = h H L(h) for L(•) a slowly varying function in zero. Here the estimation refers to the H parameter.

2 4. 2 . 2

 222 Convergence in law Proof of Theorem 3.2. Fix m ∈ N * , c 1 > 0, c 2 > 0, . . . , c m > 0 and d 1 , d 2 ,. . . , d m ∈ R and S g (εc) := m i=1 d i S g (εc i ).

  g M (εc) -S g (εc)] 2 = 0. Now, since N (0; σ 2 g M ,m (c, d)) ⇒ M →∞ N (0; σ 2 g,m (c, d)), applying Lemma 1.1 of[START_REF] Dynkin | Self-intersection gauge for random walks and for Brownian motion[END_REF], Theorem 3.2 is proved.

2 4

 2 (|N * |)].

1 0

 1 log(|Z ε (u)|) du = E[log(|N * |)] + √ εS g log (ε),(27)and by Lemma 3.1 following Theorem 3.2,E[ √ εS g log (ε)] 2 = O(ε).Thus (26) follows and from the definition of M log (ε) (see[START_REF] Ho | Limiting distributions of nonlinear vector functions of stationary Gaussian processes[END_REF]), one obtainsM log (ε) = (H -2) log(ε) + log(σ 2H ) + E[log |N * |] + o p(1). (28)Now, by (11), (28) and using property[START_REF] Feuerverger | Estimation of fractal index and fractal dimension of a Gaussian process by counting the number of level crossings[END_REF], we get Ĥlog -2 = l i=1 z i M log (εc i ) = (H -2) + o p (1).

1 0

 1 log(|Z ε (u)|) du = (2 -H) log(ε) -log(σ 2H ) + M log (ε),(29)one hasE[M log (ε)] = E[log(|N * |)] + (H -2) log(ε) + log(σ 2H ),

1 0 1 0[log 1 0[ 1 0

 1111 |N * | β(ε) ] β(ε) -E[log(|N * |)] = O(β(ε)),(36)and β(ε) = o( √ ε), (35) holds. Thus we follow the same lines as in the proof of Corollary 3.2 to obtain( Ĥβ(ε) -H) c i S(εc i ) + o p (1),where Ĥβ(ε) is defined in[START_REF] Moulines | Broad band log-periodogram regression of time series with long range dependence[END_REF]. Lemma 4.3 and Corollary 3.3 allow us to conclude the asymptotic normality of Ĥβ(ε) . Let us show now that Ĥβ(ε) is an asymptotically unbiased estimator of H. As in the proof of Corollary 3.2 (see (17) and (18)), it is enough to prove that|Z ε (u)| β(ε) du ] -→ ε→0 E[log(|N * |)].(37)To show this above convergence, we use the fact that log is a concave function and thenβ(ε)E[log(|N * |)] = E[ log(|Z ε (u)| β(ε) ) du] E|Z ε (u)| β(ε) du ] log E|Z ε (u)| β(ε) du] = log E[|N * | β(ε) ] .

  ). Theorem 3.2 gives the required result. Now, to conclude the proof of Corollary 3.2, we have to compute the coefficients in the Hermite expansion of function g β (x) = |x| β E[|N * | β ] -1. To simplify the notation, first let us compute the coefficients in the Hermite expansion of h β

The research of the second author was suported in part by the project "Modelaje Estocástico Aplicado" of the Agenda Petróleo of FONACIT Venezuela.

2) Now, (29), (30), ( 27) and (31) entail that

Theorem 3.2 gives the required result. Now, to conclude the proof of Corollary 3.3, we have to compute the coefficients in the Hermite expansion of function

For the last integral we used that H m (x)φ(x) = (-1) m d m dx m (φ(x)), and integration by parts. Now the definition of the Hermite polynomials and (25) yield

Consequently, by formula for Γ(p + 1/2), one has

which gives the required result. 2

where

As in [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF], we have the following lemma

Proof of Lemma 4.3. Lemma 3.1 following Theorem 3.2 gives the asymptotic behavior of E[S g log (ε)] 2 , set

Now, (9) in Corollary 3.2 gives the expression of the Hermite's coefficients, ĝ2n,β(ε) ,

and as in the proof of Lemma 3.1 following Theorem 3.2, Mehler's formula [START_REF] Berzin | Convergence in fractional models and applications[END_REF]

The expression of the Hermite coefficients, ĝ2n,log , of g log is given in ( 13) in Corollary 3.3 and we observe that ĝ2n,β(ε)

Furthermore for small enough ε ĝ2n,β(ε)

To achieve the proof of Lemma 4.3, we compute E[ S(ε)S g log (ε)] by Mehler's formula and we proceed as for E[ S(ε)] 2 to obtain