Laurent Gosse

Analysis and short-time extrapolation of stock market indexes through projection onto discrete wavelet subspaces

Keywords: Multi-resolution analysis, wavelet decomposition, scale-limited extrapolation, empirical finance. 1991 MSC: 91B84, 94A12, 42C40, 60G25

come

Introduction

Extrapolation is one of the fundamental problems arising in Signal Processing; during the last decades, a lot of work has been devoted to it, in the particular context of band-limited signals, see e.g. [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Youla | Image Restoration by the Method of Convex Projections: Part 1 Theory[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Marks | Alternating Projections onto Convex Sets[END_REF][START_REF] Potter | Energy concentration in band-limited extrapolation[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF]. The theoretical basis for band-limited extrapolation has to do with the classical Paley-Wiener theorem which states that a function belonging to L 2 (R) whose Fourier transform has compact support can be extended to the complex plane as an entire function of exponential type. Hence, its knowledge restricted to any connected Email address: l.gosse@area.ba.cnr.it (Laurent Gosse).

Preprint submitted to Elsevier Science

open set still allows for a reconstruction up to any arbitrarily big domain by relying onto classical analytic functions theory. Numerically, things are more involved, especially because this analytic extension turns out to be an ill-conditioned problem which is quite sensitive to noise and truncation errors; in particular, methods based on Taylor series have been completely given up for Fast Fourier Transform (FFT) routines, see again [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Marks | Alternating Projections onto Convex Sets[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF].

In this paper, we follow a similar methodology to produce reliable extrapolations of liquid stock market indexes relying on so-called scale-limited approximation. Roughly speaking, in the case of band-limited extrapolation, one assumes that the signal under consideration can be approximated by a function of compact support in the Fourier space (hence having only finite frequency oscillations) whereas in the present context, one considers that the signal can be well represented by its projection onto a subspace of L 2 (R) containing functions with rather "thick" oscillations. Such subspaces can be found by considering what is now usually called a Multi-Resolution Analysis of L 2 (R) (see Definition 1 below): they are spanned by bases of shift-invariant scaling functions with various smoothness properties [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Hence this work applies primarily to indexes which are dominated by their low-frequency approximation and for which fine scale features can be considered as detail not being of great importance. This is not a very restrictive assumption since the composition of most of modern stock indexes is regularly updated to keep their variations under control even when some share's prices may have erratic behavior. For instance, on june 18, 2007, AGF and Thomson have been removed from the French index CAC 40; they have been replaced by Air France-KLM and Unibail. Similarly, the Dow Jones Industrials 30 Index will see its composition modified as General Motors should leave because its capitalization dropped too strongly during march 2009. An estimate of the volatility of an investment can be given by a statistical measure known as the standard deviation of the return rate. Standard deviation can be thought as being synonymous with volatility in such a context. An S&P 500 index fund has a standard deviation of about 15%; a standard deviation of zero would mean a return rate that never varies, like a bank account paying compound interest at a guaranteed rate. Exchange rates on FOREX markets seem to display quite an opposite behavior and the methods presented in this paper don't seem to apply (see e.g. [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF] for some empirical differences between stock and FX markets, especially the gain/loss asymmetry).

Recently, a whole theory has been proposed in a series of papers [START_REF] Xia | On the limit of sampled signal extrapolation with a wavelet signal model[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF][START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF][START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF] (see also [START_REF] Croft | Wavelet-based signal extrapolation[END_REF]) to extend the well-known band-limited extrapolation theory to a more flexible framework of scale-limited spaces which have the nice property to be able to handle finite-duration signals. The roadmap really follows classical Fourier extrapolation and proceeds through linear integral operator theory (see e.g. [START_REF] Zeidler | Applied Functional Analysis[END_REF] for a general reference): this is recalled in §2.1. However, one big difference is that from a compact set of observations, one cannot theoretically reconstruct the signal on the whole real line (except in the special case where the multi-resolution analysis involves bases of analytic functions), see Theorem 3. So one can merely hope to reconstruct a small part of the missing signal on both sides of the measures interval; hence the "short-time extrapolation" in the title. The size of this "small part" is related to the smoothness of the analyzing wavelets and the scale at which one decides to work (and this scale itself is related to the volatility of the index).

In §2.2, we present the standard algorithms to produce these extrapolations: they originate from the now classical Gerchberg-Papoulis algorithm [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF] widely used for band-limited extrapolation and whose convergence has been studied in various papers, see [START_REF] Donoho | Uncertainty principles and signal recovery[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Marks | Alternating Projections onto Convex Sets[END_REF][START_REF] Youla | Image Restoration by the Method of Convex Projections: Part 1 Theory[END_REF] and references therein. In the present case, this routine converges toward the so-called Minimumnorm Solution of the extrapolation problem. Because of its poor practical convergence properties, an acceleration by the Conjugate Gradient (CG) algorithm is generally needed; this is explained in [START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF] in case of band-limited extrapolation, and in [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF] for scale-limited extrapolation. We don't completely follow this paper in §2.3 as we propose an easier way to produce a CG acceleration.

Our specific application is considered in §3. First, in §3.1, we explain that because of the features of the market indexes, particular wavelet bases are needed. Namely, since these indexes are strongly non-periodic, one cannot consider usual wavelet bases defined on the whole real line as they process finite duration signals either by zero-padding or by artificial periodization. Both techniques will lead to the appearance of strong discontinuities at the edge of the observation interval, which will be reflected in the wavelet coefficients at very fine scales. Hence such bases will yield poor representations at medium scales, especially in the vicinity of the borders. A solution is to use the wavelet bases on a bounded interval: these bases have been first proposed by Meyer [START_REF] Meyer | Ondelettes sur l'intervalle[END_REF], then by Cohen, Daubechies, Jawerth and Vial [START_REF] Cohen | Multiresolution analysis, wavelets and fast algorithms on an interval[END_REF][START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF]; see also e.g. [START_REF] Cheng | About wavelets on the interval[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Melkemi | Wavelet bases on the interval: a new approach[END_REF][START_REF] Siong Lee | Signal and Image Approximation using Interval Wavelet Transform[END_REF]. Then, in §3.2, we explain through an example that results can be of quite poor quality when scale-limited extrapolation is applied to individual share prices displaying locally strong variations on a fine scale. Smoothness and weak volatility are therefore requested to produce an acceptable result.

In §4, we present some numerical results on real-life data. First, in §4.1, a quiet bull market on the French CAC 40 is considered: this gives the best conditions in which the scale-limited extrapolation can be used because the market follows a stable trend and the volatility is low. Then, in §4.2, we look at the real-estate EPRA Eurozone index as a prototype of a moderately volatile stock index in a bear market; the quality of the extrapolation is lower, but there are signs (like the size of some wavelet coefficients) which give a hint about the difficulty in treating such a situation. Finally, we considered what happened on the Hang-Seng Chinese index during the first days of october 2008 when it lost around 25% of its value in less than 10 days; obviously, extrapolation routines cannot give reliable information in such a case, but the analysis of what is left behind after projecting the index onto a scalelimited subspace of L 2 (R) can somewhat reveal itself as a warning sign its high volatily (which makes violent corrections plausible in a difficult economic period). In Appendix A, we briefly present results on scale-limited extrapolation of assets log-returns (see [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF]).

Let us recall finally that other works previously addressed the wavelet techniques for stock market data: see the review paper [START_REF] Ramsey | The contribution of wavelets to the analysis of economic and financial data[END_REF] and the book [START_REF] Gençay | An Introduction to Wavelets and Other Filtering Methods in Finance and Economics[END_REF]. Concerning this book, the authors works primarily with the discontinuous Haar wavelet which in general forbids to use any of these scale-limited extrapolation techniques. Other forecasting algorithms are proposed in [START_REF] Aminghafari | Forecasting time-series using wavelets[END_REF][START_REF] Ario | Forecasting Time Series via Discrete Wavelet Transform[END_REF][START_REF] Fryzlewicz | Forecasting non-stationary time series by wavelet process modeling[END_REF][START_REF] Renaud | Wavelet-based Forecasting of Short and Long Memory Time Series[END_REF][START_REF] Yu | Time Series Forecasting using Wavelets with Predictor-Corrector Boundary Treatment[END_REF]. A book is devoted to wavelet techniques for time series analysis, see [START_REF] Percival | Wavelet Methods for Time Series Analysis[END_REF]. Some results of empirical finance are reviewed in [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF].

2 From band-limited to scale-limited extrapolation

In this section, we aim at reviewing some of the results from [START_REF] Xia | On the limit of sampled signal extrapolation with a wavelet signal model[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF][START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF][START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF] which will be useful in the context of processing stock market indexes. Let's begin by introducing the concept of Multi-Resolution Analysis (MRA): (see e.g. [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for details)

Definition 1 A sequence of nested subspaces V j is called a Multi-Resolution Analysis of L 2 (R) if: {0} ⊂ • • • ⊂ V -1 ⊂ V 0 ⊂ V 1 ⊂ • • • ⊂ L 2 (R).
Moreover, the following properties must hold:

• for all f ∈ L 2 (R), P Vj f -f L 2 → 0 as j → +∞ also, P Vj f → 0 as j → -∞. • if f (t) ∈ V j , then f (t/2) ∈ V j-1 and for all k ∈ Z, f (t -2 j k) ∈ V j .
• there exists a shift-invariant orthonormal base of V 0 given by the scaling function

φ n (t) = φ(t -n) for n ∈ Z.
In this définition, P Vj stands for the orthogonal projector onto the subspace V j . Intuitively, it asks for the V j 's to be vector subspaces of functions with increasing temporal resolution: when j decreases, functions in V j tend to become constants. Oppositely, when j increases, they are allowed to oscillate with high instantaneous frequency. The wavelet spaces W j are defined as the orthogonal complement of V j inside V j+1 , that is to say: for all j ∈ Z,

V j+1 = V j ⊕ W j .
From φ n , the base of V 0 , one can deduce a base of V j by simple dilatation,

φ j,n (t) = √ 2 j φ n (2 j t) = √ 2 j φ(2 j t -n). (1)
Thus, the orthogonal projection of f onto V j reads:

P Vj f = n∈Z < f, φ j,n > φ j,n , < f, φ j,n >= R f (t)φ j,n (t).dt, (2)
which is the best approximation of f in V j in the least-squares sense. Similarly, for any value of T ∈ R + , we denote by P T the time-truncation operator:

P T f (t) = f (t) if |t| ≤ T ,
zero in the other case. Finally, χ A stands for the indicator function of any set A.

Reproducing kernel Hilbert spaces (RKHS) and continuous extrapolation

The general problem of extrapolation for scale-limited functions can be presented as follows: we fix a MRA in the sense of Definition 1 and an index j ∈ Z. Any function f ∈ V j can be represented by means of the scaling function through the formula (2); in particular, if we choose a compact time interval [-T, T ], with T ∈ R + , we can define g, a "truncated version" of f the following way: (notice that in general,

g ∈ V j ) ∀t ∈ [-T, T ], g(t) = (P T • P Vj f )(t) = n∈Z < f, φ j,n > φ j,n (t). (3)
Implicitly, this defines a linear bounded operator H j : V j → L 2 (-T, T ). The problem of extrapolating a time-truncated, scale-limited signal g can be formulated as finding a unique solution f ∈ V j to H j f = g. In order to proceed with this inversion, we recall from linear operator theory the adjoint H * j which we can compute easily:

< Hf, ϕ > = < n∈Z < f, φ j,n > φ j,n , ϕ > = n∈Z < f, φ j,n >< φ j,n χ t∈[-T,T ] , ϕ > = < f, n∈Z φ j,n < φ j,n , ϕ χ t∈[-T,T ] >> = < f, H * ϕ > .
Thus H * is a projector L 2 (-T, T ) → V j ; in particular, the composition

H j • H * j defines an integral operator from L 2 (-T, T ) onto itself, ∀g ∈ L 2 (-T, T ), H j • H * j g (t) = T -T g(s)Q j (s, t).ds,
where Q j is called the Reproducing Kernel for the reproducing kernel Hilbert space [START_REF] Yao | Applications of reproducing kernel Hilbert spaces-Band-limited signals models[END_REF] V j :

Q j (s, t) = n∈Z φ j,n (s)φ j,n (t) = Q j (t, s), (s, t) ∈ R 2 . ( 4
)
In the sequel, we assume Q j to be bounded1 as a function of R 2 and continuous in [-T, T ] 2 .

Lemma 1 (see [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF]) The following properties hold for Q j bounded and continuous:

• if the scaling function φ is real-valued, H j • H * j is self-adjoint and positive-semidefinite; in particular, all its eigenvalues are real and nonnegative, so

+∞ > λ 1 ≥ λ 2 ≥ ... ≥ 0. • H j • H * j is compact; in particular, its set of eigenfunctions r k (t) satisfying (H j • H * j )r k = λ k r k forms an orthogonal basis of L 2 (-T, T ).
The second part of the lemma is classical: see for instance §4.4 (page 240) of [START_REF] Zeidler | Applied Functional Analysis[END_REF]. These eigenfunctions are the key to extrapolate scale-limited signals. In case some eigenvalues vanish, let K stand for the biggest integer for which there holds:

k ≤ K ⇒ λ k > 0 and λ K+1 = 0.
We then extend the domain of definition of the eigenfunctions

r k from [-T, T ] to R: recall that λ k r k (t) = T -T r k (s)Q j (s, t).ds.
Hence, by definition of the reproducing kernel Q j (4), we can build for all k ≤ K,

∀t ∈ R, rk (t) = 1 λ k n∈Z T -T r k (s)φ j,n (s) φ j,n (t) ∈ V j .
It turns out that these modified eigenfunctions have interesting properties:

Lemma 2 (see [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF]) For Q j bounded and continuous, the set of functions rk

∈ V j , k ≤ K is orthonormal in L 2 (R)
and orthogonal in L 2 (-T, T ); their corresponding eigenvalues are positive and lower than 1. Moreover, there holds:

∀f ∈ V j , T -T f (s)r k (s).ds = λ k R f (s)r k (s).ds. (5)
In this context, λ k can be seen as the energy contribution of rk inside the interval [-T, T ]; to see this, we choose f = rk inside the formula ( 5), apply Cauchy-Schwarz inequality, and use rk L 2 (R) ≥ rk L 2 (-T,T ) which leads to:

rk L 2 (R) ≤ 1 λ k r k L 2 (-T,T ) .
We can now define U j , the linear subspace of L 2 (R) spanned by the set rk , k ≤ K.

Theorem 3 (see [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF]) The following decomposition holds for any j ∈ Z:

V j = U j ⊕ V j ∩ P T L 2 (R) ⊥ , U j = P Vj • P T L 2 (R) .
As a consequence of the definition of U j ⊂ V j , any function f ∈ V j can be written as

f 1 + f 2 with f 1 ∈ U j and f 2 ∈ U ⊥ j . Obviously, f 1 (t) = k≤K a k rk (t).
The fact is that, from (5), it is readily seen that the coefficients a k depend on P T f only:

a k = R f (s)r k (s).ds = 1 λ k T -T f (s)r k (s).ds.
At this point, there is no clue to compute the missing part f 2 which belongs to the set of functions in V j vanishing on [-T, T ]. However, we can get rid of this part by a smoothness argument if we aim at extrapolating the signal on a small enough interval. Indeed, for finite j, the space V j ∩ L 2 (-T, T ) ⊥ cannot contain discontinuous functions; hence there is a small interval beyond T , say [T, T + δT ], inside which f 2 must vanish. The value of δT is not known precisely, but clearly, δT is lower than the support of φ j,0 . Hence, in case the following rough estimate holds,

|δT | < |supp(φ j,0 )|, (6)
Theorem 3 ensures that a signal g = P T f , f ∈ V j , can be extrapolated uniquely up to [-T -δT, T + δT ]. The value of δT thus decreases when j is getting bigger, hence a rule of thumb could be stated like "the more details one decides to keep inside the signal, the shorter the extrapolation length will be". It can happen that stock indexes, during periods of high volatility, ask for a big value of j in order to keep track of some of the small scale features: in such a case, the length of the extrapolation interval given by δT will be reduced accordingly through [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF] (see an illustration in §4). Remark 1 In case the scaling function is analytic, the space V j ∩ L 2 (-T, T ) ⊥ is reduced to zero because a non-zero analytic function cannot vanish on a non-trivial interval. In this case, we recover the usual framework of band-limited extrapolation (see [START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF]).

Generalized Gerchberg-Papoulis (GGP) algorithm and scale-limited extrapolation

The leading idea of the series of papers [START_REF] Xia | On the limit of sampled signal extrapolation with a wavelet signal model[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF][START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF][START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF] is to extend the classical Gerchberg-Papoulis algorithm, originally proposed for band-limited extrapolation (see for instance [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Youla | Image Restoration by the Method of Convex Projections: Part 1 Theory[END_REF][START_REF] Marks | Alternating Projections onto Convex Sets[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Donoho | Uncertainty principles and signal recovery[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF]), to the processing of scale-limited functions belonging to wavelet subspaces. To begin with, let's recall the original construction: a function f ∈ L 2 (R) is said to be W-band-limited if its Fourier transform f has compact support. In such a case, one has:

∀t ∈ R, f (t) = F -1 ( f ) (t) := W -W f (ξ) exp(2iπtξ).dξ with F : f → f .
The number W ∈ R + is called the band-width of f . Moreover, the classical Paley-Wiener theorem ensures that f is actually an entire function of exponential type which can be extended to the whole complex plane as a result of analytic continuation theory. This paves the way to various extrapolation routines as we explain now; first, we denote by P W the projector onto the subspace of W-band-limited functions. More precisely,

∀f ∈ L 2 (R), P W (f ) = F -1 f χ [-W,W ] .
By analogy with the preceding section, we define g as a time-truncated version of a Wband-limited signal f as follows:

∀t ∈ [-T, T ], g(t) = P T • P W (f ) (t).
Similarly, this defines a bounded linear operator H W from the subspace of W-band-limited functions onto its restriction to functions defined in [-T, T ]. The extrapolation problem can be formulated as finding the unique W-band-limited signal f for which H W f = g. In this direction, we introduce a reproducing kernel Q W for which the following equality holds [START_REF] Yao | Applications of reproducing kernel Hilbert spaces-Band-limited signals models[END_REF]:

∀f ∈ L 2 (R), P W (f ) = f, f (t) = R f (s)Q W (s, t).ds, t ∈ [-T, T ].
Since f is known to be W-band-limited, an explicit form can be derived for Q W :

Q W (s, t) = W -W exp(2iπ(t -s)ξ).dξ = sin(2π(t -s)W ) π(t -s) = Q W (t, s),
which turns out to be the low-pass filtering kernel. The whole procedure can be repeated for this particular choice of Q W corresponding indeed to the analytic scaling function being the cardinal sine, "sinc". The eigenfunctions of H W • H * W are called the Slepian Prolate Spheroidal wave functions (PSWF), see for instance [START_REF] Moore | Prolate Spheroidal wave functions, an introduction to the Slepian series and its properties[END_REF][START_REF] Zayed | A generalization of the prolate spheroidal wave functions[END_REF], and are known to realize the maximum energy concentration for a W-band-limited function onto the interval [-T, T ]. They constitute an orthogonal basis of both L 2 (-T, T ) and a subspace of L 2 (R), known as the Paley-Wiener space of W-band-limited functions. This means in particular that there's no need for a decomposition like the one of Theorem 3. This theoretical construction doesn't lead to a simple numerical algorithm able to generate the W-band-limited extrapolation of g. However, if we denote by P T c the orthogonal projector such that P T c (f )(t) = f (t) if |t| > T , zero in the other cases, we can observe that:

∀t ∈ R, g(t) = (P T • P W (f ))(t) = (Id -P T c • P W )(f )(t).
From [START_REF] Donoho | Uncertainty principles and signal recovery[END_REF], we know that the operator (Id -P T c • P W ) can be invertible on the space of Wband-limited functions. In this case, (Id -P T c • P W ) -1 = k≥0 (P T c • P W ) k . An iterative procedure to recover f from g originally due to Gerchberg and Papoulis [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF] reads:

f (0) = g, f (1) = g + P T c • P W (f (0) ), ... , f (k+1) = g + P T c • P W (f (k) ). ( 7
)
Several proofs of convergence are available, see [START_REF] Youla | Image Restoration by the Method of Convex Projections: Part 1 Theory[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Donoho | Uncertainty principles and signal recovery[END_REF]. The discrete algorithm is shown to construct the so-called Minimum-Norm Least-Squares (MNLS) solution in [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF], see also [START_REF] Potter | Energy concentration in band-limited extrapolation[END_REF].

At this level, it clear how to derive an iterative algorithm for the preceding problem of scalelimited extrapolation following the same stepping stones: it suffices to replace the projector P W by P Vj inside the preceding formula: (see [START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF])

f (0) = g, f (1) = g + P T c • P Vj (f (0) ), ... , f (k+1) = g + P T c • P Vj (f (k) ). ( 8
)
This procedure is called the Generalized Gerchberg-Papoulis algorithm (GGP) in [START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF]. We finally give the main convergence result from [START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF]:

Theorem 4 (see [START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF]) Let φ be an orthogonal real-valued scaling function and f (k) , k ∈ N be the sequence of functions generated by the generalized Gerchberg-Papoulis algorithm (8) for some f ∈ V j , j ∈ Z. Then, as k → +∞, f (k) converges toward the minimum norm solution f defined by:

• f ∈ V j • f L 2 (R) = inf ϕ L 2 (R) ; ϕ ∈ V j , P T ϕ = P T f
Recalling Theorem 3, it is likely that the minimum norm solution f belongs to U j alone; actually, f is the orthogonal projection of f onto U j (see Corollary 1 in [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF]).

Fig. 2. Illustration of the alternating projections algorithm. Then, since F 1 ∩ (F 2 ) ⊥ = 0 as a consequence of the fact that analytic functions have isolated zeros, the iterative procedure [START_REF] Croft | Wavelet-based signal extrapolation[END_REF] converges toward the unique point in F 1 starting from its orthogonal projection on F 2 (see Figure 2). Clearly, the transversality condition isn't met in the general case of scale-limited subspaces (especially because

F 1 ∩ (F 2 ) ⊥ = V j ∩ P T L 2 (R) ⊥ = 0).

Acceleration by Conjugate Gradient resolution

The iterative procedures from the preceding section can be readily implemented on a computer; however, it is a well-known fact that they are endowed with rather poor practical convergence properties. In order to fix this drawback, several acceleration routines have been proposed: we shall concentrate onto the ones based on the Conjugate Gradient algorithm (see [START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF] for the band-limited extrapolation and [START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF] for the scale-limited one). The idea is the same in both cases of band-limited and scale-limited extrapolation: one considers the pointwise collection of sampled values for g(t), and while seeking for an approximation of f (t), aim at solving efficiently the linear system given by g(t) = P T • P Vj (f )(t) (for band-limited extrapolation, one inserts P W in place of P Vj in the former formula).

In this direction, it is necessary to define g ∈ R M as the collection of N observed values (the samples of g(t)) completed by zero (M > N ), the discrete operators T N and P j being respectively the diagonal truncation matrix and the discrete projector onto V j ; the unknown is therefore f ∈ R M such that g = T N P j f . Following [START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF], we begin by observing that these matrices have several interesting properties: in case the base of V j given by the scaling function φ is orthonormal, the matrix P j is symmetric. Moreover, we clearly have that (T N P j ) * = P j T N and T * N T N = T N . The Conjugate Gradient method is to be applied to the normal equations yielding the least-squares solution of the discrete extrapolation problem:

x = inf x∈R M Ax -b L 2 ⇔ (A * A)x = A * b.
In the present case, this reduces to find efficiently f solving (T N P j ) * g = P j T N g = P j g = (T N P j ) * (T N P j ) f

= P * j T * N T N P j f = P j T N P j f = P j T N f ∈ U j ,
see Theorem 3. Intermediate steps are done by noting that T N g = g and P j f = f by their very definition. Since P j T N P j is symmetric positive definite (observe that x T P j T N P j x = (P j x) T T N (P j x) > 0 for a smooth scaling function φ in case (6) holds), the Conjugate Gradient algorithm can be set up to solve in a finite number of operations the linear system: b := P j g = P j T N P j f := Ax. We recall it briefly: Fix r 0 = b and p 0 = r 0 ; while r k > ε,

• q k = Ap k • If k=1 p 1 = r 0 , else β k-1 = r ⊤ k-1 r k-1 r ⊤ k-2 r k-2 p k = r k-1 + β k-1 p k-1 • α k = r ⊤ k-1 r k-1 p ⊤ k q k • x k = x k-1 + α k p k • r k = r k-1 -α k q k
This algorithm will hopefully converge in a finite number of iterations toward the vector x ≃ f which contains the sampled values of the minimum-norm solution f of the scalelimited extrapolation problem up to a given (small) tolerance 0 < ε ≪ 1. We close this section by noticing that in practice, it is not necessary to implement the matrix P j : the WaveLab package contains commands for computing the Fast Wavelet Transform and its inverse. It suffices thus to use it, to cut the coefficients corresponding to scales finer than j, and to apply the inverse transform. For big values of M , it will also result in a quicker execution.

Application to stock market indexes

One theoretical obstruction for considering stock market indexes in the context of multiresolution analysis of L 2 (R) is that these signals aren't typically functions of finite energy defined on the whole real line. Hence a way around this issue is to extend them with zero (zero-padding) outside of their definition interval (which is obviously compact); the resulting function belongs to L 2 (R) and we can apply the preceding developments. Moreover, if considering time scales not shorter than daily quotations, stock market signals can be treated as Lipschitz continuous W 1,∞ (R) functions (with possibly big Lipschitz constants depending on volatility). These are the only hypotheses really necessary to be made on the signal in order to set up the preceding algorithms.

The need for wavelet bases on a bounded interval

As a consequence of the aforementioned minor issue concerning our signals, it may result difficult to set up the usual wavelet algorithms; in particular, since stock indexes are strictly positive and strongly non-periodic, extending them either by zero of periodically will result in big discontinuities on the border of the interval of observations. This leads typically to the emergence of big and spurious wavelet coefficients at the high resolution scales which render for these discontinuities (in practice, standard algorithms for wavelet transforms 2 extend the source signal to the whole real line periodically). Hence, in order to be able to get a correct compression of such a signal by projecting it onto a scale-limited space, the use of wavelet bases on a bounded interval became absolutely necessary. The theory of multi-resolution on a bounded interval has been developed after the pioneering work by Meyer [START_REF] Meyer | Ondelettes sur l'intervalle[END_REF]. However, his construction couldn't be used to produce a robust numerical algorithm. This problem has been fixed later by Cohen, Daubechies, Jawerth and Vial [START_REF] Cohen | Multiresolution analysis, wavelets and fast algorithms on an interval[END_REF][START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF]; see also e.g. [START_REF] Cheng | About wavelets on the interval[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Melkemi | Wavelet bases on the interval: a new approach[END_REF][START_REF] Siong Lee | Signal and Image Approximation using Interval Wavelet Transform[END_REF]. This results in a local modification of the scaling functions at the edge of the interval so that all of the properties enjoyed by the wavelets on R are preserved. In particular, as used in the characterization of spaces C s (R) via the decay of wavelet coefficients, it is important to preserve both the wavelet's smoothness and its vanishing moments. This corresponds to introducing modified filters (the sequences usually denoted by {h n } and {g n } in the cited papers) for the treatment of the edges of the interval; otherwise, away from the borders, the algorithm remains essentially the same as for the case of the infinite real line. On Figure 3, one can visualize in the middle column the same scaling function and mother wavelet than in Figure 1; however, in both left and right columns, their modifications able to handle properly discontinuities possibly appearing at the edges of a finite interval are displayed. We stress the fact that, since we aim at extrapolating data on the right edge of the observations interval by projecting the whole set onto a scale-limited subspace, it is of crucial importance to prevent the apparition of high wavelet coefficients in this region because of artificial periodicity. In this last case, the scaling function coefficients kept when annihilating all the others corresponding to finer scales are likely to be modified because of this "Gibbs-like" phenomenon. Consequently, the initial data given to the GGP algorithm ( 8) would be rather inaccurate, and would thus result in a poor extrapolation since such an inverse problem is quite sensitive to noise (see [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF]).

Multi-scale decomposition and scale-limited extrapolation

Hereafter, the notation I(t) ∈ R + will refer to some nonnegative stock index and I, to its sampled values on a time-discrete grid. In practice, we shall always rely on daily quotations, which are freely available at http://fr.finance.yahoo.com. The general framework presented in this paper can be set up by first, selecting a multi-resolution analysis on a bounded interval which is compatible with I (in the sense that base functions should have vanishing moments and smoothness properties compatible with the volatility of the index under consideration), then choosing a scale j ∈ Z which realizes a good balance between accuracy and data compression. With both these elements in hand, it is possible to perform short-time extrapolation as explained in the previous sections on real-life examples.

According to the WaveLab package, only the standard compact support Daubechies wavelet with either 2 or 3 vanishing moments is available (see again Figure 3 for an illustration). For most of the indexes, the one with only 2 vanishing moments gave the best results for observations ranging from 512 to 2048 days of quotations (recall that one year gives roughly 250 quotation days). Once the multi-resolution analysis has been selected, the remaining issue is to choose a scale parameter j ∈ Z which retains enough information in the time series, but is low enough to allow for a non-trivial extrapolation. This means that P Vj I -I L 2 should be small; in other words, the wavelet coefficients associated to subspaces W ℓ with ℓ > j should be rather small (the smaller, the better). In this case, one can say that the index I is well compressed in the base φ j,n of V j . This has to do with the fact that oscillations in the index at a scale finer than j don't really carry important information; volatility (in the sense formerly mentioned of standard deviation) may exist, but what matters is that its effect onto the global trend of I is very limited in such a situation. This favorable context doesn't happen systematically: on Figure 4, we display the 1024-days evolution (ending on october 17, 2008), wavelet coefficients and conjugate gradient convergence which produce a 20-days extrapolation for the famous Buffett's holding, Berkshire-Hathaway. In this rather spectacular case, one sees clearly that during most of the observation interval, prices follow rather well the trend obtained by projecting them onto V 5 ;

however, in the last days, a very rapid growth cancels the validity of such an approximation in a scale-limited space. The consequences are easily seen on the right figure: in blue are represented the coefficients of φ j=5,n . In red, green, light blue and black are represented the wavelet coefficients in W 5 , W 6 , W 7 and W 8 , W 9 multiplied by 10 respectively (the abscissa is the scale index). Obviously, for this pathological example, we clearly have important energy concentrated in the fine scales: hence the conjugate gradient algorithm converges poorly, and the resulting extrapolation is quite poor because the scale-limited space doesn't really match all the meaningful features of the evolution of the share price. As a consequence, working on indexes is safer than to do so on individual share's prices, mainly because the former are more likely to display a smoother behavior because they are computed as convex combinations of the latter.

4 Numerical experiments on real data

A quiet bull market

We begin by considering the French leading index CAC 40 on a period of 1024 days (roughly 4 years beginning in 2003) covering the end of the bear market which resulted from the "internet krach" and the September 11 events. The very low interest rates, cut by FED Chairman Alan Greenspan, helped to sustain the prices and allowed for a bull market to take place. It was characterized by a moderate growth and a low volatility; this is easily noticeable by looking at the wavelet coefficients in the middle/right part of Figure 5 since the coefficients in W j , j > 5 look quite small even if they are again multiplied by 10. This is the best context in which to set up our scale-limited extrapolation algorithms. For a 30 days forecast (roughly one month and a half), one gets a monotone convergence of the conjugate gradient routine in less than 10 iterations. On the top graphic of Figure 5, one sees that starting from the red arrow (at the abscissa t = 1004), the MNLS extrapolated part matches quite well the evolutions of the index. Obviously, one cannot hope to obtain exactly the real oscillations, as what is sought is an approximation in V 5 only, but the idea is more to extrapolate a general trend and to define a confidence interval around it which should contain the fluctuations of the near-future.

In order to set up this confidence interval, one considers the discarded part of the signal (see bottom left of Figure 5), namely I(t) -P Vj I(t) (which has some vanishing moments because it belongs to the wavelet spaces ∪ k W k , k ≥ j), computes its distribution function together with its standard deviation. In case it would be a standard Gaussian distribution (see for instance the bottom of Figure 5), then it would suffice to trace both red curves at a distance of 2 times the standard deviation to ensure a 95% confidence interval. In realistic experiments, we found that the discarded signal approaches a Gaussian distribution when overall volatility was low thus one can only speak about an approximated confidence interval.

The stopping criterion ε for the CG algorithm was 10 -9 .

A moderately volatile bear market

During years 2003-2007, not only the shares prices have tremendously increased; real-estate market value has inflated more than stocks in certain areas as a consequence of wide and easy credit allocation. Some of the gains obtained on financial markets have also been recycled in "concrete" during this period. On Figure 6, we are interested in the FTSE EPRA Eurozone index (EPRA meaning European Real Estate Association 3 ), which contains es-sentially shares linked to european real-estate markets, like banks, insurance companies and promoters. Obviously, starting from august 2007, this index entered into a bear market as a consequence of the so-called subprime crisis4 ignited in the USA in february 2007 with the depreciations of the bank HSBC.

We work with a measures interval of 512 days (roughly 2 years, ending in february 2008) mainly because this index is too recent to allow for a 1024 days sample. A bear market being in general more volatile than a bull one (one says that shares prices decrease twice as fast as they grow [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF]), the Gaussian approximation of the "noise" defined as I(t) -P Vj I(t) with j = 4 is less convincing; however, the approximated confidence interval remains rather satisfying even in the clearly decreasing period ranging from day 200 to 500. From the repartition of the energy in the wavelet coefficients (middle in Figure 6, it may sound reasonable to work in the V 5 space rather than in V 4 as we did. However, this wouldn't lead to a sensible improvement in the quality of the results. On the top of Figure 6, one sees that the MNLS 15 days extrapolation (starting from the red arrow at the abscissa t = 497) is rather good in determining a trend around which a reliable confidence interval can be built. The convergence of the CG algorithm was satisfying; the stopping criterion ε for the CG algorithm was 10 -10 .

What happens for a krach ?

On Figure 7, we display what happened on the Chinese index Hang-Seng in the first days of october 2008: namely, because of recession fears in both the USA and the Eurozone, operators sold massively emerging countries shares regardless of the so-called decoupling theory 5 . The index under consideration lost around 25% of its value in less than 10 days and authorities decided to close the market several times and to ban short-selling in order to prevent a complete collapse. In such a case, there's no hope that the scale-limited extrapolation (starting from the red arrow at the abscissa t = 1004) can give a reliable result (see top of Figure 7); however, what can be considering as interesting is that such an analysis highlights the fact that the discarded part of the signal (see bottom of Figure 7) has a now repartition function which is different from the standard Gaussian: it is quite peaked which a thick tail on the left side. Hence this suggests a potentially volatile index which is likely to have violent corrections. Increasing the value of j may improve the result by including more fine scale features inside P Vj I(t): however, a strong reduction of the size of the extrapolation interval would also come as a direct consequence of (6).

Conclusion and outlook

We presented in this paper a quite general methodology to handle a standard problem in Empirical Finance, that is, trying to extrapolate stock market indexes during periods of mild volatility relying on the nowadays well-known theory of multi-resolution and discrete wavelet subspaces of L 2 (R). The numerical results obtained on real-life data coming from daily quotations indicate that this process can be considered as reliable as long as the basic assumption of the relevance of a scale-limited approximation P Vj I of the index I holds. Especially, during strong corrections (like e.g. in october 2008), one notices rather easily that I -P Vj I ceases to be small and is even not distributed like a standard Gaussian density of the same average and variance. When this "empirical density" is very peaked (like in Figure 7), the scale-limited extrapolation cannot be taken as a meaningful possible forecast for the index. Concerning the size of the extrapolation interval, it depends on the level of resolution j of the scale-limited approximation (a lower j gives a longer extrapolation, but at the price of a poorer accuracy) and on the smoothness of the scaling function (smoother yields a bigger interval). The extreme case being the one of analytic functions from which one can theoretically extrapolate onto the whole real line R).

A Treating the very oscillating "log-returns" function

There exists another way, maybe more natural, to extrapolate a stock index I(t): it is to work with its so-called "log-returns" function (see e.g. [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF]) defined as follows:

S(t) = ln I(t) I(t -1) = ln(I(t)) -ln(I(t -1)) 1 ≃ d dt ln(I(t)) = I ′ (t) I(t) . (A.1)
As the quotations are taken on a daily base, a straightforward piecewise linear interpolation produces an overall positive Lipschitz continuous function belonging to the Sobolev space W 1,∞ . Therefore, extracting its logarithmic derivative makes sense at least within the framework of L ∞ , or maybe even in functions of bounded variation. The advantage is that, since S is a very oscillating function around zero, the need for multi-resolution bases defined on a bounded interval may be relaxed as (in general) one shouldn't expect too big variations in the L ∞ norm. Hence the discontinuities possibly showing up at the edges of the interval of observations can be really smaller compared to the ones arising when working with the index I itself. This can open the way to using much smoother wavelet bases, like for instance the Symmlet: see After the separation of S between low and fine scales in order to produce P Vj (S), I(0) becomes useless, but one can use another value. Let C j be the support of the scaling function φ j,n=0 (t), as an initial point, we take: Ĩ0 = 1 |C j | Cj I(s).ds.

This allows to compute a scale-limited index at level j like:

∀t ≥ 0, Ĩj (t) := Ĩ0 exp t 0 P Vj (S)(s).dx .

We applied these ideas to the example of EPRA Eurozone (see §4.2); on the top of Figure A.2, we display the multi-resolution decomposition of its log-returns S(t) ranging from V 3 , W 3 up to W 5 . On the middle, the extrapolation of the index obtained through the scalelimited extrapolation of its log-returns function is shown (staring from the red arrow at abscissa t = 497). It is considerably smoother than the result from Figure 6; the first reason is that the analyzing wavelet is smoother, the second one is that moreover, S is integrated to deduce I. On the bottom, one can see that the convergence of the Conjugate Gradient algorithm is satisfying even for this type of extrapolation problem. However, when looking at the wavelet coefficients, one may think it became more difficult to justify why the space V 4 can be chosen to approximate the function S(t) as the fine scale coefficients aren't small like in the examples shown in the preceding sections. But since we work now with a derivative, big values on a fine scale correspond to both important and rapid (maybe erratic) variations of the index I(t). It is therefore perfectly normal to remove them in order to derive a smooth trend; the working subspace V j should correspond now to a set of moderately big coefficients in the discrete wavelet transform of S.

Fig. 1 .

 1 Fig. 1. Compact support of a scaling function and its mother wavelet.

Remark 2

 2 It is possible to give a direct convergence proof for the band-limited iterative extrapolation algorithm relying on the so-called Alternating Projections Theorem in Hilbert spaces. One defines two subsets of L 2 (R): F 1 is the subspace of W-band-limited functions, and F 2 is made of functions vanishing outside of [-T, T ].

Fig. 3 .

 3 Fig. 3. Scaling function (top) and its mother wavelet (bottom) in the interval.

Fig. 4 .

 4 Fig. 4. Quotations of Berkshire-Hathaway A (left) and wavelet coefficients (right).

0 S

 0 Figure A.1. However, a small issue arises because, since all the filtering and the extrapolation is done on the logarithmic derivative (A.1), the starting point I(0) is lost; indeed, by definition, there holds: ∀t ≥ 0, I(t) = I(0) exp t (s).dx .

Fig. 5 .Fig. 6 .

 56 Fig. 5. CAC 40 index on 1024 days: 30 days extrapolation (top), scale-limited decomposition and convergence (middle) and remaining noise repartition (bottom).

Fig. A. 1 .

 1 Fig. A.1. Comparison between the Daubechies 4 (top) and Symmlet 8 (bottom) basis functions.

This is the case as soon as |φ(t)| ≤ O(1 + |t| 0.5+ǫ ) -1 for some ǫ > 0, see[START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF]; especially, this is true for scaling functions with compact support.

See http://www.epra.com

See http://en.wikipedia.org/wiki/Subprime mortgage crisis

See http://en.wikipedia.org/wiki/Decoupling#Decoupling and the stock market declines of January 2