Léonard Gallardo
email: leonard.gallardo@lmpt.univ-tours.fr

Chaabane Rejeb
email: chaabane.rejeb@gmail.com

Support properties of the intertwining and the mean value operators in Dunkl's analysis

Keywords: primary: 31B05, 33C52, 47B39; secondary: 43A32, 51F15 Dunkl-Laplacian operator, Dunkl's intertwining operator, Generalized volume mean value operator and harmonic kernel, Rösler's measure, Dunkl harmonic functions

In this paper we show that the Dunkl intertwining operator has a compact support which is invariant by the associated Coxeter-Weyl group. This property enables us to determine explicitely the support of the volume mean value operator, a fundamental tool for the study of harmonic functions relative to the Dunkl-Laplacian operator.

Introduction and statement of the results

Let R be a (finite) root system in R d with associated Coxeter-Weyl group W (see [START_REF] Humphreys | Reflection groups and Coxeter groups[END_REF] or [START_REF] Kane | Reflection Groups and Invariant Theory[END_REF] for details on root systems) and for ξ ∈ R d , let D ξ be the Dunkl operator defined by

D ξ f (x) = ∂ ξ f (x) + ∑ α∈R + k(α) ⟨α, ξ⟩ f (x) -f (σ α x) ⟨α, x⟩ , f ∈ C 1 (R d ),
where R + is a subsystem of positive roots, σ α is the reflection directed by the root α ∈ R + , k is a nonnegative multiplicity function (defined on R) and ∂ ξ f is the usual ξ-directional derivative of f . These operators, introduced by C. F. Dunkl in the nineties (see [START_REF] Dunkl | Differential-difference operators associated to reflection groups[END_REF]), are related to partial derivatives by means of an intertwining operator V k (see [START_REF] Dunkl | Hankel transforms associated to finite reflection groups[END_REF] or [START_REF] Dunkl | Orthogonal Polynomials of Several variables[END_REF]) as follows

∀ ξ ∈ R d , D ξ V k = V k ∂ ξ . (1.1)
We know that V k is a topological isomorphism from the space C ∞ (R d ) (carrying its usual Fréchet topology) onto itself satisfying (1.1) and V k (1) = 1 (see [START_REF] Trimèche | The Dunkl intertwining operator on spaces of functions and distributions and integral representation of its dual[END_REF]) and V k commutes with the W -action (see [START_REF] Rösler | Dunkl Operators: Theory and Applications[END_REF]) i.e.

∀ f ∈ C ∞ (R d ), ∀ g ∈ W, g -1 .V k (g.f ) = V k (f ), (1.2)
where g.f (x) = f (g -1 x).

A fundamental fact due to M. Rösler (see [START_REF] Rösler | Positivity of Dunkl's intertwining operator[END_REF] or [START_REF] Rösler | Dunkl Operators: Theory and Applications[END_REF]) is that for every x ∈ R d , there exists a unique compactly supported probability measure µ k

x on R d with

supp µ k x ⊂ C(x) := co{gx, g ∈ W } (1.3)
(the convex hull of the orbit of x under the group W ) such that

∀ f ∈ C ∞ (R d ), V k (f )(x) = ∫ R d f (y)dµ k x (y). (1.4)
Note that the property (1.3) has been proved in [START_REF] De Jeu | The Dunkl transform[END_REF]. Throughout this paper, the notation k > 0 means that k(α) > 0 for all α ∈ R.

Concerning the measure µ k x (which we call Rösler's measure at point x), the first result of our paper is the following

Theorem A. For every x ∈ R d , we have 1) x ∈ supp µ k x , 2) If k > 0, the support of µ k x is a W -invariant set, 3) If k > 0 then W.x (the W -orbit of x) is contained in supp µ k x .
A question strongly related to the support of Rösler's measures concerns the Dunklmean value operator introduced by the authors in [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF] in the study of harmonic functions for the Dunkl-Laplacian operator ∆ k = ∑ d i=1 D 2 i where D i = D e i with (e i ) 1≤i≤d an orthonormal basis of R d . Precisely for x ∈ R d and r > 0, the mean value of a continuous function f at (x, r) is defined by

M r B (f )(x) := 1 m k (B(0, r)) ∫ R d f (y)h k (r, x, y)ω k (y)dy,
where y → h k (r, x, y) is a compactly supported measurable function (see Section 2) given by

h k (r, x, y) := ∫ R d 1 [0,r] ( √ ∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩)dµ k y (z), (1.5)
m k is the measure dm k (x) := ω k (x)dx and ω k is the weight function

ω k (x) := ∏ α∈R + ⟨α, x⟩ 2k(α) . (1.6)
In particular we have shown that a

C 2 (R d )-function u is ∆ k -harmonic in R d if and only if for all (x, r) ∈ R d ×R + , u(x) = M r B (u)(x)
. For a further thorough study of ∆ k -harmonicity on a general W -invariant open set, it would be crucial to get information on the support of the mean value operators. We already know that the measures

dη k x,r = 1 m k (B(0, r)) h k (r, x, y)ω k (y)dy (x ∈ R d , r > 0), (1.7)
are probability measures with compact support equal to supp h k (r, x, .) and satisfying the following inclusion ( [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF]):

supp h k (r, x, . ) ⊂ B W (x, r) := ∪ g∈W B(gx, r), (1.8)
where B(x, r) denotes the usual closed ball of radius r centered at x. In fact, the second main result of this paper, intimately related to Theorem A, is a precise description of the support of h k (r, x, .). It states that Theorem B: Let x ∈ R d and r > 0.

1)

We have B(x, r) ⊂ supp h k (r, x, .).

2)

If k > 0, then we have supp h k (r, x, .) = B W (x, r) := ∪ g∈W B(g.x, r).
We will call B W (x, r) the closed Dunkl ball centered at x and with radius r > 0 associated to the Coxeter-Weyl group W .

The harmonic kernel and the mean value operator

In this section we recall some results of [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF].

Let (r, x, y) → h k (r, x, y) be the harmonic kernel defined by (1.5). We note that in the classical case (i.e. k = 0), we have µ k y = δ y and h 0 (r, x, y) = 1 [0,r] (∥x -y∥) = 1 B(x,r) (y). The harmonic kernel satisfies the following properties (see [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF]):

1. For all r > 0 and x, y ∈ R d , 0 ≤ h k (r, x, y) ≤ 1.

2. For all fixed x, y ∈ R d , the function r -→ h k (r, x, y) is right-continuous and non decreasing on ]0, +∞[.

3. Let r > 0 and x ∈ R d . For any sequence (φ ε ) ⊂ D(R d ) of radial functions such that for every ε > 0, 0 ≤ φ ε ≤ 1, φ ε = 1 on B(0, r) and ∀ y ∈ R d , lim ε→0 φ ε (y) = 1 B(0,r) (y),
we have

∀ y ∈ R d , h k (r, x, y) = lim ε→0 ∫ R d φ ε ( √ ∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩)dµ k y (z),
where φ ε is the profile function of

φ ε i.e. φ ε (x) = φ ε (∥x∥).
4. For all r > 0, x, y ∈ R d and g ∈ W , we have

h k (r, x, y) = h k (r, y, x) and h k (r, gx, y) = h k (r, x, g -1 y). (2.1)
5. For all r > 0 and x ∈ R d , we have

∥h k (r, x, .)∥ k,1 := ∫ R d h k (r, x, y)ω k (y)dy = m k (B(0, r)) = d k r d+2γ d + 2γ , ( 2.2)
where d k is the constant

d k := ∫ S d-1 ω k (ξ)dσ(ξ) = c k 2 d/2+γ-1 Γ(d/2+γ) .
Here dσ(ξ) is the surface measure of the unit sphere S d-1 of R d and c k is the Macdonald-Mehta constant (see [START_REF] Opdam | Dunkl operators, Bessel functions and the discriminant of a finite Coxeter group[END_REF], [START_REF] Etingof | A uniform proof of the Macdonald-Mehta-Opdam identity for finite Coxeter groups[END_REF]) given by

c k := ∫ R d e -∥x∥ 2 2 ω k (x)dx. 6. Let r > 0 and x ∈ R d . Then the function h k (r, x, .) is upper semi-continuous on R d .
7. The harmonic kernel satisfies the following geometric inequality: if ∥a -b∥ ≤ 2r with r > 0, then

∀ ξ ∈ R d , h k (r, a, ξ) ≤ h k (4r, b, ξ)
(see [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF], Lemma 4.1). Note that in the classical case (i.e. k = 0), this inequality says that if ∥a -b∥ ≤ 2r, then B(a, r) ⊂ B(b, 4r).

8. Let x ∈ R d . Then the family of probability measures dη k x,r (y) defined by (1.7) is an approximation of the Dirac measure δ x as r -→ 0. That is

∀ α > 0, lim r→0 ∫ ∥x-y∥>α dη k x,r (y) = 0
and if f is a locally bounded measurable function on a W -invariant open neighborhood of x and if f is continuous at x, then (see [START_REF] Gallardo | A new mean value property for harmonic functions relative to the Dunkl-Laplacian operator and applications[END_REF], Proposition 3.2):

lim r→0 ∫ R d f (y)dη k x,r (y) = lim r→0 M r B (f )(x) = f (x). (2.3)

Proof of the results

For convenience we group together the first points of Theorem A and Theorem B in the following Proposition.

Proposition 3.1 Let x ∈ R d . Then i) for every r > 0, x ∈ supp h k (r, x, .), ii) x ∈ supp µ k x ,
iii) for every r > 0, B(x, r) ⊂ supp h k (r, x, .).

Proof: i) Suppose that there exists r > 0 such that x / ∈ supp h k (r, x, .). Then we can find ε > 0 such that h k (r, x, y) = 0, for all y ∈ B(x, ε). Let f be a nonnegative continuous functions on R d such that supp f ⊂ B(x, ε) and f = 1 on B(x, ε/2). Since t → h k (t, x, y) is increasing on ]0, +∞[, we deduce that

∀ t ∈]0, r], 0 ≤ M t B (f )(x) ≤ 1 m k [B(0, t)] ∫ R d f (y)h k (r, x, y)ω k (y)dy = 0.
Hence, we obtain M t B (f )(x) = 0, for all t ∈]0, r]. Letting t → 0 and using the relation (2.3), we get a contradiction.

ii) Let x ∈ R d be fixed. At first, we claim that

∀ r > 0, ∀ y ∈ R d , h k (r, x, y) ≤ µ k x [B(y, r)]. (3.1)
Indeed, from the inclusion supp µ k x ⊂ B(0, ∥x∥), we see that

∀ y ∈ R d , ∀ z ∈ supp µ k x , ∥y -z∥ 2 ≤ ∥y∥ 2 + ∥x∥ 2 -2 ⟨y, z⟩ .
This implies for any y ∈ R d and r > 0 that

∀ z ∈ supp µ k x , 1 [0,r] ( √ ∥y∥ 2 + ∥x∥ 2 -2 ⟨y, z⟩ ) ≤ 1 [0,r] (∥y -z∥) = 1 B(y,r) (z).
If we integrate the two terms of the previous inequality with respect to the measure µ k x , we obtain h k (r, y, x) ≤ µ k x ( B(y, r)

) and then (3.1) follows from (2.1). Now, if x / ∈ supp µ k x , there exists ϵ > 0 such that µ k x ( B(x, ϵ) ) = 0. Thus, we have µ k x ( B(y, ϵ/2) ) = 0 whenever y ∈ B(x, ϵ/2). Using (3.1), we deduce that h k (ϵ/2, x, .) = 0 on B(x, ϵ/2), a contradiction with the result of i).

iii) Let y ∈ R d such that ∥x -y∥ < r. As lim z→y (∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩) = ∥x -y∥ 2 , there exists η > 0 such that √ ∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩ ≤ r for every z ∈ B(y, η).

Therefore, by using the fact that y ∈ supp µ k y we obtain

h k (r, x, y) ≥ µ k y [B(y, η)] > 0. Remark 3.1 For α ∈ R, let H α := {x ∈ R d , ⟨x, α⟩ = 0}
be the hyperplane directed by α. Note that in [START_REF] Rösler | Short-time estimates for heat kernels associated with root systems[END_REF] (Corollary 3.6) and under the condition x / ∈ ∪ α∈R H α , Rösler has proved that x ∈ supp µ k x by using the asymptotic behavior of the Dunkl kernel E k (x, y) which is defined by

E k (x, y) := V k ( e ⟨.,y⟩ ) (x) = ∫ R d e ⟨z,y⟩ dµ k x (z).
We turn now to the second statement of Theorem A that we recall below: Theorem 3.1 Let x ∈ R d and assume that k > 0. Then the set supp µ k x is W -invariant.

Proof: In order to simplify the formulas, we will assume here that the root system R is normalized i.e. ∥α∥ 2 = 2 for all α ∈ R. In particular, for reflections we have σ α x = x -⟨α, x⟩ α. We will prove that if y ∈ supp µ k x , then σ α y ∈ supp µ k x for every α ∈ R. Let then y ∈ supp µ k x and suppose that there is a root α ∈ R such that σ α y / ∈ supp µ k x . Write y ′ := σ α y to simplify notations. There is a ball B(y

′ , ϵ) (ϵ > 0) such that for all f ∈ C ∞ (R d ) with compact support included in B(y ′ , ϵ), we have ∫ R d f (z)µ x (dz) = V k f (x) = 0. Let us denote by C ∞ y ′ ,ϵ (resp. C y ′ ,ϵ ) the set of all functions f ∈ C ∞ (R d ) (resp. f ∈ C(R d )) with compact support in B(y ′ , ϵ). For all ξ ∈ R d and all f ∈ C ∞ y ′ ,ϵ , we also have ∂ ξ f ∈ C ∞ y ′ ,ϵ . By the intertwining relation (1.1) we obtain ∀ ξ ∈ R d , ∀ f ∈ C ∞ y ′ ,ϵ , D ξ V k f (x) = 0. Suppose f ∈ C ∞ y ′ ,ϵ and f ≥ 0 and let g := V k f . We have g ≥ 0 on R d (because V k preserves positivity) and ∀ ξ ∈ R d , D ξ g(x) = ∂ ξ g(x) + ∑ α∈R + k(α)⟨α, ξ⟩ g(x) -g(σ α x) ⟨x, α⟩ = 0. (3.2)
But as g(x) = 0, x is a minimum of g so ∂ ξ g(x) = 0 and relation (3.2) implies

∀ ξ ∈ R d , ∑ α∈R + k(α)⟨α, ξ⟩ g(x) -g(σ α x) ⟨x, α⟩ = 0. (3.3) Now, consider the set R x := {α ∈ R + ; x ∈ H α }.
There are two possible locations for x:

• First case: Suppose that R x = ∅ i.e x / ∈ ∪ α∈R H α (i.e. for all root α ∈ R, ⟨x, α⟩ ̸ = 0). Applying (3.3) with ξ = x and using the fact that g(x) = 0, we get

∑ α∈R + k(α)g(σ α x) = 0.
As g ≥ 0 and k > 0, we obtain that g(σ

α x) = V k f (σ α x) = 0 for all α ∈ R + and all f ∈ C ∞ y ′ ,ϵ
and f ≥ 0. By uniform approximation, we deduce that for all f ∈ C y ′ ,ϵ and f ≥ 0, we also have V k f (σ α x) = 0. Finally for every f ∈ C y ′ ,ϵ , by decomposing f = f + -f -with f + = max(f, 0) and f -= -min(f, 0) and using the linearity and W -equivariance of V k (relation (1.2)), we obtain that

∀ f ∈ C y ′ ,ϵ , ∀ α ∈ R + , V k f (σ α x) = V k (σ α .f )(x) = 0, where σ α .f is the function z → f (σ α z). As it is easy to see that σ α .C y ′ ,ϵ = C σαy ′ ,ϵ , we deduce that ∀ α ∈ R + , ∀f ∈ C σαy ′ ,ϵ , V k f (x) = 0.
But this implies in particular that V k f (x) = 0 for all f ∈ C y,ϵ in contradiction with the hypothesis y ∈ supp µ k x . The result of the theorem follows in the first case. • Second case: Suppose that R x ̸ = ∅. For every β ∈ R x , clearly we have x = σ β x. Therefore, since g(x) = 0, we get g(σ β x) = 0, for all β ∈ R x . But, as x is a minimum of g, we have

∀ β ∈ R x , g(x) -g(σ β x) ⟨x, β⟩ = ∫ 1 0 ∂ β g(x -t ⟨x, β⟩ β)dt = ∂ β g(x) = 0.
Hence, the relation

(3.3) with ξ = x implies ∑ α∈R + \Rx k(α)g(σ α x) = 0.
Consequently, we obtain g(σ α x) = 0 for all α ∈ R. The end of the proof of the first case applies and gives also the result in this case. This completes the proof of the theorem.

From the W -invariance property of the support of µ k x and the fact that x ∈ supp µ k x , we obtain immediately the last assertion of Theorem A: Proof: Let g ∈ W and y ∈ R d such that ∥gx -y∥ < r. We will proceed as in the proof of Proposition 3.1, iii). We have lim

z→g -1 y √ ∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩ = ∥x -g -1 y∥.
Hence, there exists η > 0 such that for all z ∈ B(g -1 y, η), √ ∥x∥ 2 + ∥y∥ 2 -2 ⟨x, z⟩ ≤ r and thus h k (r, x, y) ≥ µ k y [B(g -1 y, η)]. But, from the fact that g -1 y ∈ supp µ k y we deduce that y ∈ supp h k (r, x, .). This completes the proof. Remark 3.2 When k ≥ 0, we will say that a root α ∈ R is active if k(α) > 0. Let us denote by R A = {α ∈ R; k(α) > 0} the set of active roots and F the vector subspace of R d generated by {α, α ∈ R A }. Then we can generalize the results of Theorems A and B in the following form a) The set R A is a root system. Indeed, using the fact that k is W -invariant, we can see that for every α

, β ∈ R A , k ( σ α β ) = k(β) > 0. Thus ∀ α ∈ R A , R A ∩ Rα = {±α} and σ α (R A ) = R A . b) Let W A be
= ξ ′ + ξ ′′ ∈ F + F ⊥ = R d (where F ⊥ is the orthogonal complement of F in R d ). • Let x ∈ R d . Rösler's measure µ k
x is of the form (see [START_REF] Rösler | Asymptotic analysis for the Dunkl kernel[END_REF])

µ k x = µ k A x ′ ⊗ δ x ′′ , ( 3.5
)

where µ k A x ′ is Rösler's measure associated to (R A , k A ) and δ x ′′ is the Dirac measure at x ′′ . We have supp µ k x = x ′′ + supp µ k A x ′ . From (1.3), the support of µ k A x ′ is contained in the convex hull of W A .x ′ (the W A -orbit of x ′
). Furthermore, by Theorem A, it is invariant under the action of the group W A and contains the whole orbite W A .x ′ .

• Let x ∈ R d and r > 0. According to (1.5)

and (3.5) the harmonic kernel is given by

h k (r, x, y) = ∫ R d 1 [0,r] ( √ ∥x ′′ -y ′′ ∥ 2 + ∥x ′ ∥ 2 + ∥y ′ ∥ 2 -2 ⟨x ′ , z ′ ⟩ ) dµ k A y ′ (z ′ ), y ∈ R d .
The support of h k (r, x, .) takes the following form supp h k (r, x, .) = x ′′ + B W A (x ′ , r) = x ′′ + ∪ g∈W A B(gx ′ , r) = ∪ g∈W A B(gx, r).

Example 3.1 Let (e 1 , e 2 ) be the canonical basis of R 2 . Then, the set R := {±e 1 , ±e 2 } is a root system in R 2 , its Coxeter-Weyl group is Z 2 2 and the multiplicity function can be identified to a pair k = (k 1 , k 2 ), with k i = k(e i ) ≥ 0, i = 1, 2. Take x = (x 1 , x 2 ) ∈ R 2 with x 1 , x 2 > 0. In this case, according to [START_REF] Xu | Orthogonal polynomials for a family of product weight functions on the spheres[END_REF], Rösler's measure is given by µ k x = µ k 1 x 1 ⊗ µ k 2 x 2 , where µ k i x i = δ x i if k i = 0 and

⟨µ k i x i , f ⟩ = Γ(k i + 1/2) √ πΓ(k i ) ∫ 1 -1 f (tx i )(1 -t) k i -1 (1 + t) k i dt.
if k i > 0 (see [START_REF] Dunkl | Integral kernels with reflection group invariance[END_REF]).

• If k = (0, 0), µ k x = δ x and h k (r, x, y) = 1 B(x,r) (y).

• If k = (k 1 , 0) with k 1 > 0, then supp µ k
x is the line segment between x and σ e 1 x = (-x 1 , x 2 ) and supp h k (r, x, .) = B(x, r) ∪ B(σ e 1 x, r).

Corollary 3 . 1 Corollary 3 . 2

 3132 Let x ∈ R d and assume that k > 0 . Then, for all g ∈ W , gx ∈ supp µ k x . Now, we can turn to the proof of the second statement of Theorem B. Let x ∈ R d and r > 0. If k > 0, then supp h k (r, x, .) = B W (x, r) := ∪ g∈W B(gx, r).(3.4)

  the Coxeter-Weyl group associated to the root system R A . Then the restriction k A of k to R A is clearly invariant under the W A -action. In other words, it is a multiplicity function.

c) For any ξ ∈ R d , we will use the notation ξ

•

  If k 1 , k 2 > 0, the support of µ k x is the convex hull of Z 2 2 .x and the closed Dunkl ball is given by

	B Z 2 2 (x, r) = supp h k (r, x, .) = B	(	(x 1 , x 2 ), r	)	∪B	(	(-x 1 , x 2 ), r	)	∪B	(	(x 1 , -x 2 ), r	)	∪B	(	(-x 1 , -x 2 ), r	)	.


Madalina Olteanu
email: madalina.olteanu@univ-paris1.fr

A Descriptive Method to Evaluate the Number of Regimes in a Switching Autoregressive Model

Keywords: switching autoregressive models, hierarchical clustering, Ward distance, SOM

This paper proposes a descriptive method for an open problem in time series analysis : determining the number of regimes in a switching autoregressive model. We will translate this problem into a classification one and define a criterion for clustering hierarchically different model fittings. Finally, the method will be tested on simulated examples and real-life data.

Introduction

In the past few years, several nonlinear autoregressive models were proposed for time series analysis. Some of these models are based on the idea that the process is characterized not by a unique linear autoregression, but by the fact that two or more regimes are driving the series behaviour. In each regime, an autoregressive function is fitted. We are interested in the case where the autoregressive functions are linear in every regime. The most classical examples are TAR (Treshold Autoregressive) models introduced by [START_REF] Tong | On a treshold model[END_REF] with regime switching according to the magnitude of a treshold variable, the smoothed version of TAR models (STAR), or the more recent Markov switching autoregressive models, first used by [START_REF] Hamilton | New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle[END_REF] to model the U.S. Gross National Product. Estimating the parameters of these models is usually done by maximizing the likelihood function, but under a very strong hypothesis, a fixed number of regimes. Choosing the "true" number of regimes is still an open problem, as this is equivalent to testing with lack of identifiability under the null hypothesis. This leads to a degenerated Fisher information matrix and thus the chi-square theory and the likelihood ratio tests fail to apply. An empirical method to detect this kind of non-linearity using Kohonen maps and hierarchical clustering of linear regressions is given below. The second section describes the method, while the third provides examples on simulated and real-life data. Finally, a conclusion will follow.

The Method

The problem of finding the "true" number of regimes can be rewritten as a classification problem by using a sliding window as follows . Suppose that we have observed the values of a time series {y t } t=1,T and we decide to fit an autoregressive model. Once the order of the model has been determined (with an AIC criterion, for example), we can consider the data set of dimension (Tp)x(p + 1), {y t , y t-1 , ...y t-p } t=p+1,T . Looking for the number of regimes is actually equivalent to looking for the number of regression lines (or hyperplanes) which will best fit the data. The idea is simple and is based on the possibility of finding patterns in data which will identify the regression hyperplanes. Given the data set in Figure 1, fitting one regression line to the data is clearly not the good choice. If we now suppose that we managed to cluster the data into two groups and we perform a regression within each of these groups, we get two lines which seem to describe better the sample. This is confirmed by the "within squared error", which is equal to the sum of squared residuals if there is only one regime and, if there are several, to the total sum of squared residuals within each group. In the general case, we would like to start with some "good" initial clusters which will be then classified hierarchically using some squared error criterion and we would expect to have an important break in the increasing values of this criterion, once we pass from the true number of regimes to a smaller one.A "good" cluster should contain observations belonging to the same regime and, at the same time, have enough points to estimate a regression line.

Initial clustering

For the initial clustering, self-organizing maps were employed. This choice was motivated by several properties of the Kohonen maps such as the vector quantization property (the last steps of the algorithm use "0-neighbours"and thus are equivalent to classical VQ techniques), the homogeneity of the clusters as well as the topology conservation used to fasten the hierarchical classification algorithm.

One problem arises once we get the clusters and try to fit a regression within each one : are there enough points in every cluster? No cluster will be allowed to have less points than the number of lags or regressor variables. For this, either we eliminate from the analysis those which do not verify this condition, either we force the points to move to a different cluster, by assigning them to the closest sufficiently large cluster.Since very few observations are concerned with this problem and in order to shorten the computing time, the first approach was preferred.

The assumption we make here is that the property of homogeneity of self-organizing maps manages to create clusters in which observations belong to one regime. This could be justified by the fact that the variables used for the classification contain information concerning the regime of the observation and similar profiles will belong to the same cluster. Indeed, in the simulated examples where the different regimes are known, we'll see that the map clusters are generally homogeneous from this point of view. This property will no longer apply if the regression hyperplanes are too close and the noise is important.

Hierarchical classification

As we actually need to compare different data fits, which is also equivalent to different numbers of hyperplanes, we need to adapt a hierarchical classification to our case (let us first make the convention to call "clusters" the result of the Kohonen map and "classes" two or more "clusters" joined together by the hierarchical method). We will choose a new "distance" between classes by developping a squared error criterion.

A very popular method used in classification is to minimize a within-class variation criterion, the variation within a class being defined as the sum of squared distances from the individuals to the barrycenter. By considering a data set x i ∈ R n , i = 1, ..., N and a fixed number of clusters k, the total within-class variation can be written as I w = k i=1 w i I i , where w i is the weight of class i, I i = N i j=1 d 2 (x i,j , g i ) being the variation within the class i and g i is the corresponding barrycenter. We denoted by N i the number of individuals in class i and by x i,j the data point indexed by j in class i. In hierarchical classification, this principle was adapted by Ward and the algorithm consists in clustering together the individuals who minimize the increase of the within-class variation. Our idea was to build an algorithm similar to Ward's, but, as our interest is to estimate the number of hyperplanes characterizing the data, we will replace the barrycenters by regression lines and the within-class variation becomes the within sum of squared errors. Obviously, in this case, the inter-class variation cannot be defined. For a fixed number of classes, k, the within sum of squared errors is defined as 2 is the sum of squared residuals and ŷt is the predicted value of y t by the linear regression of order p fitted in class C l , l = 1, k. Now, in the frame of hierarchical clustering, when passing from k to k -1 classes, if classes i and j were clustered, the within sum of squared errors becomes :

SSE w,k = k l=1 SSE C l , where SSE C l = t∈C l (y t -ŷt )
SSE i,j w,k-1 = k l=1,l =i,l =j SSE C l + SSE C i ∪C j
Following the same principle as Ward's, we want to minimize the increase in the within inertia, which in this case is defined by the within sum of squared errors. This is equivalent to finding i and j which minimize

∆S i,j w,k,k-1 = SSE i,j w,k-1 -SSE w,k = SSE C i ∪C j -SSE C i -SSE C j 2.2.

The within sum of squared errors criterion

Now, let us take a closer look at this difference and give an explicit expression as a function of the data and the residuals. The following notations are necessary :

-

Y i = {Y t } t∈C i ∈ R n i , Y j = {Y t } t∈C j ∈ R n j
, where n i and n i are the cardinalities of C i and, respectively, C j

-X i = {1, Y t-1 , ..., Y t-p } t∈C i ∈ R n i ×(p+1) , X j = {1, Y t-1 , ..., Y t-p } t∈C j ∈ R n j ×(p+1) .
The linear regressions fitted in classes C i and C j can be written as :

Y i = X i • β i + u i = X i • βi + e i Y j = X j • β j + u j = X j • βj + e j ,
where β i , β j , βi , βj ∈ R p+1 , βi and βj are the least squares estimates of β i and β j , u i , e i ∈ R n i and u j , e j ∈ R n j are the error and, respectively, the residuals vectors, u i ∼ N (0, σ 2 i I n i ) and u j ∼ N (0, σ 2 j I n j ). The linear regression fitted in the joint class C i ∪ C j is written as :

Y = X • β + u = X • β + e, where Y = Y i Y j ∈ R n i +n j , X = X i X j ∈ R (n i +n j )×(p+1) , β, β ∈ R p+1
, β is the least squares estimate of β, u, e ∈ R n i +n j are the error and the residuals vector and u ∼ N (0, Ω), Ω =

σ 2 i I n i 0 0 σ 2 j I n j .
Then, we can compute ∆S i,j w,k,k-1 as :

• ∆S i,j w,k,k-1 = SSE C i ∪C j -SSE C i -SSE C j = e e -(e i e i + e j e j )
Remark 1 : Using this form, [START_REF] Toyoda | Use of the Chow Test under Heteroscedasticity[END_REF] proves that ∆S i,j w,k,k-1 is approximately distributed as σ 2 χ 2 (p + 1), where σ 2 is any well-chosen weighted average of σ 2 i and σ 2 j and p is the number of lags considered.

• ∆S i,j w,k,k-1 = SSE C i ∪C j -SSE C i -SSE C j = Y -X • β 2 -Y i -X i • βi 2 -Y j -X j • βj 2 (1) But Y -X • β = Y i -X i • β Y j -X j • β = Y i -X i • βi Y j -X j • βj + X i • β i -X i • β X j • βj -X j • β (2) and thus Y -X • β 2 = Y i -X i • β Y j -X j • β 2 = Y i -X i • βi Y j -X j • βj 2 + X i • β i -X i • β X j • βj -X j • β 2
since the cross-product on the right term in (2) can easily be seen to be zero. We get that :

Y -X • β 2 = Y i -X i • βi 2 + Y j -X j • βj 2 + X i • β i -X i • β 2 + X j • βj -X j • β 2 (3)
and by replacing (3) in (1) :

∆S i,j w,k,k-1 = X i • β i -X i • β 2 + X j • βj -X j • β 2 = = ( βi -β) X i X i ( βi -β) + ( βj -β) X j X j ( βj -β) (4)
Using the form of the least squares estimators βi , βj and

β βi -β = β i -β + (X i X i ) -1 X i , 0 -(X i X i + X j X j ) -1 [X i , X j ] • e i e j βj -β = β j -β + 0, (X j X j ) -1 X j -(X i X i + X j X j ) -1 [X i , X j ] • e i e j
and (4) becomes

∆S i,j w,k,k-1 = (β i -β) X i X i (β i -β) + (β j -β) X j X j (β j -β)+ +(β i -β) X i e i -X i X i (X i X i + X j X j ) -1 (X i e i + X j e j ) + + e i X i -(e i X i + e j X j ) (X i X i + X j X j ) -1 X i X i (β i -β)+ +(β j -β) X j e j -X j X j (X i X i + X j X j ) -1 (X i e i + X j e j ) + + e j X j -(e i X i + e j X j ) (X i X i + X j X j ) -1 X j X j (β j -β)+ +e i X i (X i X i ) -1 X i e i + e j X j (X j X j ) -1 X j e j - -(e i X i + e j X j ) (X i X i + X j X j ) -1 (X i e i + X j e j )
If classes i and j come form the same regime, that is β i = β j = β, the increase in the within sum of sqared errors is only

∆S i,j w,k,k-1 = e i X i (X i X i ) -1 X i e i + e j X j (X j X j ) -1 X j e j - -(e i X i + e j X j ) (X i X i + X j X j ) -1 (X i e i + X j e j )
This quantity is very close to zero if the classes contain enough points. Thus, together with remark 1,we get that if the joint classes are from the same regime, the within sum of squared errors should be close to zero and if the classes are from different regimes, the larger the difference between the parameters of the two regimes, the increase in the within sum of squared errors should be larger.

The algorithm

Now, we can write the steps of the algorithm which, at the same time, clusters the data and models the dependencies within each class: We consider k going from M to 2, where M is the number of clusters resulting from the Kohonen map. For each k, we have the following steps :

• Compute the k regressions lines corresponding to the k classes • Find (i 0 , j 0 ) which minimize ∆ i,j SSE w

• Join these classes together, put k = k -1 and go to the first step.

At the last step all points are clustered together and there is a single regression line. The next thing to do is draw the dendrogram and look how the within sum of squared errors increases over the classification. As it was mentioned earlier, one might expect an important break when the number of classes is smaller than the real number of regimes.

Examples and Results

The method was tested on several nonlinear autoregressive regime switching models. Historically, the first introduced were the treshold models (we won't speak here about the other variants of these models, smoothed etc), followed by the switching Markov and next we will consider both examples.

TAR Models

The example is a TAR of order two and the coefficients were taken from the paper of Gon-zalo&Pitarkis.

y t = -3 + 0.5y t-1 -0.9y t-2 + ε t 2 + 0.3y t-1 + 0.2y t-2 + ε t , y t-2 ≤ 1.5 , y t-2 > 1.5
, where {y t } is the observed series and ε t is i.i.d. standard gaussian. Three samples containing 200, 400 and 800 points, respectively, were simulated. Let us examine the 200 points sample. The self-organizing map dimension was fixed equal to 5x5 and {y t , y t-1 , y t-2 } were the variables used for the classification. Crossing the map with a boolean variable which distinguishes whether y t-2 is above or below the treshold value allows to see that the clusters are homogeneous and each of them takes values in one regime.

If the hierarchical clustering algorithm is run on the twenty-five clusters, the squared error criterion increases as in Figure 3 and suggests that a good choice would be a two-regimes model. The hierarchical classification provides also good estimators for the parameters of the model when choosing two regimes as shown in Table 1.

The threshold value was not estimated and we'll see that the decision on the type of the model (TAR, Markov switching) cannot be made on this basis only. When increasing the number of points in the sample, we have chosen to increase also the map size (we considered a 7x7 map for the 400 sample and a 9x9 for the 800). This was done in the purpose of conserving clusters homogeneity, but a good criterion for choosing the size of the map is still to be defined. The results are very similar and in both cases the hierarchical algorithm suggests a two-regimes model.

A Two Regime Markov Switching Model

For this example, we will first define an autoregressive Markov switching process. If {y t } t∈N is the observed time series, let us suppose that it follows a linear autoregressive process of order p and that we have two regimes, the passage from one regime to the other being driven by an unobserved Markov chain, {x t } t∈N , which has a transition probability matrix A.

For two regimes and two lags of time, we can write the model as follows :

y t = f xt (y t-1 , y t-2 )+ σ xt ε t f xt (y t-1 , y t-2 ) ∈ {f 1 , f 2 }, σ xt ∈ {σ 1 , σ 2 }, ε t i.i.d.

noise (usually a standard gaussian) and

A = p 1 -p 1 -q q
is the transition probability matrix of x t .

The data used here were simulated with the parameters below (a globally stationary process was chosen):

f 1 (y t-1 , y t-2 ) = 0.2 + 0.5y t-1 + 0.1y t-2 f 2 (y t-1 , y t-2 ) = 0.3 + 0.9y t-1 -0.1y t-2
, σ 1 = 0.03 σ 2 = 0.02 and A = 0.9 0.1 0.2 0.8 As for the previous example, three samples (200, 400, 800 points) were considered. We will only list the results for the 400 sample and remark that the outputs for the other two cases were very similar. The initial clustering was performed using a 6x6 Kohonen map and {y t , y t-1 , y t-2 }as variables. Figure 4 shows that the map is well organized, the clusters are homogeneous and when crossing with the variable giving the regime, there is a good separation of them in the initial clusters. Afterwards, from the hierarchical classification of these clusters, we get, again, a huge jump when passing from two classes to one, as shown in Figure 5.The estimated coefficients in each of the two classes are shown in Table 2 (we will also note that the two clusters are homogeneous, the pourcentage of explained variance being larger than 92% in each of them). The following results are from a 400 points sample, with an initial 8x8 map. By crossing the map with the regime variable, there is a relatively good separation, although we may notice that the first two regimes seem to come closer together with respect to the third one.

y t = f xt (y t-1 , y t-2 ) + σ xt ε t , f xt (y t-1 , y t-2 ) ∈ {f 1 , f 2 , f 3 }, σ xt ∈ {σ 1 , σ 2 , σ 3 }, ε t is i.i.d. standard gaussian and    f 1 (y t-1 , y t-2 ) = 0.2 + 0.5y t-1 + 0.1y t-2 f 2 (y t-1 , y t-2 ) = 0.3 + 0.9y t-1 -0.1y t-2 f 3 (y t-1 , y t-2 ) = 0.5 + 1.2y t-1 + 0.
Here, a first conclusion would be that there are four regimes. But let us take a closer look at the hierarchical classification. One of the four final classes contains only one cluster, one cell of the map. Moreover, this cell (the 8th) is isolated from the rest of the map and it contains only four observations with very high values. If we project the data on a two or three-dimensional space, the same four observations are far from the rest. The algorithm has identified a small class of outliers which are considered as a separate regime and from this point of view the method is close to Ward's which is also sensitive to this kind of observations.

We can't continue with the examples before making an important remark. We've seen that this method of identifying the number of regression lines works quite good in the examples above. Concerning the parameter estimation and making a decision about the model (TAR or Markov switching?, for instance), the hierarchical classification provides only the estimators for the regression lines, a likelihood approach should be used instead, once we fixed the number of regimes, to estimate the rest of the parameters : threshold value, transition matrix etc. As for the second question, no theoretical result is available yet, the econometricians preferring to decide on other criteria (economic, social etc).

What about Real Life Data?

The results on the simulated examples being encouraging, we decided to run the algorithm on real data sets. Three examples were chosen, the first two are the benchmarks Old Faithful Geyser Data and Santa Fe Competition Laser Data and the third is the U.S. GNP (Gross National Product) series, used by Hamilton to introduce the switching Markov models.

Old Faithful Data

The first set of data is the classical Old Faithful Geyser in Yellowstone National Park, consisting of 299 pairs of measurements referring to the waiting time between two successive eruptions, w t , and the duration of the subsequent eruption, d t . The data was collected between August 1st and August 15th, 1985 and the two variables are recorded in minutes. Several studies of this sample are available, most authors trying to assess either the clustering of the data, either the dependency between successive events. A literature overview, as well as an analysis using time series while assuming "a priori" the existence of two patterns of dependency, can be found in Azzalini and Bowman (1990) paper. Although there are no autoregressors in this case, the problem is the same : find the number of clusters and fit a regression within each one. Our approach would be to detect the clustering of the data and, at the same time, model the dependency within each class. The idea of addressing clustering and dependency at the same time was firstly introduced by Hennig (2000) who uses regression fixed point clusters.

Here, the duration of the eruption was modeled as a linear function of the waiting time before the eruption. While plotting the duration against the waiting time, one can see that there are at least two clusters of points, depending on the waiting time. Let us make one last remark on the data, which is that due to inexact observations during the night, there are 53 points with duration=4 (long eruption) and 20 with duration=2 (short eruption). The medium eruptions (duration=3) appear only once.

For the Kohonen classification the data set {w t , d t } t=1,299 was considered (no lags of time were introduced). The map was chosen to be a 6x6 grid and 1000 iterations were performed. Afterwards, an hierarchical clustering minimizing the within sum of squared errors criterion was applied to the 33 "valid" clusters (one cluster was void and two others contained only one point).

In Figure 8, the within sum of squared errors as well as the difference ∆SSE k,k-1 are plotted. The heterogeneity of the data is obvious and one can see that at least two clusters should be considered. On the contrary, passing from three clusters to two is less obvious and we shall see next that there's overlapping of the clusters and that the regression coefficients are very close. Once we get the hierarchical classification, we get back to the self-organizing map to see how the classes spread over the grid. On the first graph in Figure 9, the tree-cut at two classes in presented : the first class in the right upper corner in light-grey circles and the second beyond the diagonal of the grid in dark-grey circles. Clusters 7 and 28 do not contain enough points and were not considered in the classification. While the 2-classes model is meaningful, in the 3-classes case things seem to be more complicated : the first class is the same, well isolated in the right upper part of the grid, while classes 2 and 3 are mixed on the grid. This is even more obvious if we plot the duration d t against the waiting time w t and identify the classes as shown in Figure 10. The same situation as for the Kohonen map occurs. In both cases, 2-classes and 3-classes model, a first class (the right upper part of the self-organizing map) is well isolated from the rest and is concentrated around the waiting=2 points. This class corresponds to short time eruptions preceeded by rather long waiting times. The second class on the first graph containts all the long duration points including all duration=4 points. This class is splitted into two when considering the 3-classes case : one sub-class is formed by the duration=4 points and the long time erruptions preceeded by long waiting times, while the second contains the data with a moderately decreasing tendency in the duration for increasing waiting times. Although there is an interpretation for the 3-classes model, the within squared error criterion used for the classification suggests a 2-classes model, because it corresponds to an important break in the increase of the squared error. Besides, in the discussion of Azzalini and Bowman (1990) geological evidence for the existence of two distinct patterns of eruptions is given and thus our conclusion is enhanced.

Santa Fe Competition Laser Series

For the laser series, a highly nonlinear data set, the algorithm selects three classes as shown in Figure 11. In order to compare the results with the existing literature, ten lags of time were used. The 10000 observations were initially clustered on a 9x9 map. Once the number of regimes was fixed, the three regimes were supposed to be the states of a discrete Markov chain and the parameters were estimated using the EM algorithm as described in [START_REF] Rynkiewicz | Hybrid HMM/MLP Models for Time Series Prediction[END_REF]. The prediction results are weaker than those obtained by [START_REF] Weigend | Nonlinear gated experts for time series : discovering regimes and avoiding overfitting[END_REF] or [START_REF] Rynkiewicz | Hybrid HMM/MLP Models for Time Series Prediction[END_REF], but the number of estimated parameters is much smaller. But, as the nonlinearities of the series are not entirely explained by a mixture of linear regressions, an adaptation of the method by replacing the linear functions with nonlinear ones could be more interesting for this kind of data.

GNP Series

Concerning the GNP series, Hamilton's approach was based on the assumption that the mean growth rate is subject to occasional, discrete shifts. We dispose of 136 trimestrial observed values of the series, from 1952 until 1984. The maximal lag to be considered was determined with the AIC criterion and was fixed at p = 3. The data was initially clustered with a 4x4 map and considering {y t , y t-1 , y t-2 , y t-3 } as variables. The sixteen clusters were grouped hierarchically by the squared error criterion and the results are displayed in Figure 12. The first graph contains the within sum of squared errors plotted against the number of classes and the sencond its percentage increase. A first break appears when considering six classes instead of seven, but this can be interpreted as being due to possible stronlgy homogeneous clusters from the same regime which get mixed. There is a second break(13%) when passing from two classes to one but this is less obvious than in previous cases and the decision to model this series by a two-regime model is questionable from our point of view. We've introduced a descriptive method to assess the presence of regime changes in nonlinear time series analysis. As there is no theoretical answer and no statistical test to solve this problem for the moment, this method may be used, but with precaution. Indeed, self organizing maps could mix the regimes if the regression hyperplanes are too close and the square error criterion seems to be sensitive to outliers. Thus, several improvements should be made in the future, like considering a different distance that would take into account the temporal dependency of the data for the initial clustering or like looking for a smarter initial clustering that would avoid mixing regimes in the same cluster. Replacing the linear regressors by nonlinear functions could also be a possibility, but then the number of parameters becomes larger and the fiability of the estimators decreases.

Figure 1 :

 1 Figure 1: Fitting clustered data

Figure 2 :

 2 Figure 2: Initial clustering with a Kohonen map for TAR model

Figure 3 :

 3 Figure 3: Squared error criterion for TAR model

Figure 4 :Figure 5 :

 45 Figure 4: Initial clustering with a Kohonen map for two regimes Markov

5y t- 2 ,

 2

Figure 6 :Figure 7 :

 67 Figure 6: Initial clustering with a Kohonen map

Figure 8 :

 8 Figure 8: Squared error criterion for the old faithful data

Figure 9 :

 9 Figure 9: The self-organizing map for the old faithful data

Figure 10 : 2 -

 102 Figure 10: 2-clusters and 3-clusters of the old faithful data

Figure 11 :

 11 Figure 11: Squared error criterion for the laser series

Figure 12 :

 12 Figure 12: Squared error criterion for GNP data

Table 1 :

 1 Coefficients for the TAR model

			Regime1			Regime2
		Intercept	y t-1	y t-2 Intercept	y t-1	y t-2
	value	-2.82	0.59	-0.89	1.4	0.37	0.32
	t-value	-19.15	13.09	15.89	8.24	5.96	4.79

Table 2 :

 2 Coefficients for the two regime Markov model3.3 A Three Regime Markov Switching ModelNow let us see what happens if we add a new regime to the model, which will moreover be explosive and drive the process into a nonstationary one. The following example was considered :

			Regime1			Regime2
		Intercept	y t-1	y t-2 Intercept	y t-1	y t-2
	value	0.3	0.88	-0.09	0.18	0.39	0.24
	t-value	41.56	32.89	-3.64	20.91	12.07	8.87