Rihab Bedoui
email: rihab.bedoui@u-paris10.fr.

Haykel Hamdi
email: haykel.hamdi@cnam.fr.

Implied Risk-Neutral probability Density functions from options prices : A comparison of estimation methods

Keywords: JEL classification : C02, C14, C65, G13 Risk-neutral density, mixture of log-normal distributions, Edgeworth expansions, Hermite polynomials, tree-based methods, kernel regression, Heston's stochastic volatility model, jump diffusion model

This paper compares the goodness-of-fit of eight option-based approaches used to extract risk-neutral probability density functions from a high-frequency CAC 40 index options during a normal and troubled period. Our findings show that the kernel estimator generates a strong volatility smile with respect to the moneyness, and the kernel smiles shape varies with the chosen time to maturity. The mixture of log-normals, Edgeworth expansion, hermite polynomials, jump diffusion and Heston models are more in line and have heavier tails than the log-normal distribution. Moreover, according to the goodness of fit criteria we compute, the jump diffusion model provides a much better fit than the other models on the period just-before the crisis for relatively short maturities. However, during this same period, the mixture of log-normal models performs better for more than three month maturity. Furthermore, in the troubled period and the period just-after the crisis, we find that semi-parametric models are the methods with the best accuracy in fitting observed option prices for all maturities with a minimal difference towards the mixture of log-normals model.

Introduction

The Risk-Neutral Density (RND) is an interesting tool for an investor who seeks to measure how market expectations evolve over time. A huge literature arose from the early 1990s on the most appropriate way to estimate this density. At the origin of all the methods, we find the famous work of [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF], who are the pioneers who determined a relationship between option prices and the RND.

Malik and Thomas (1997) have developed a nonparametric method to estimate the American option prices. This technique has been adopted by [START_REF] Bahra | Probability distributions of future asset prices implied by option prices[END_REF] and Schernick, Garcia and Tirupattur (1996) on European options. It is a representation of RND based on a linear combination of log-normal density fonctions. These authors estimated the RND of the oil price during the Gulf crisis and found that the estimated density is different from the density obtained by standard methods. Söderlind and Svensson (1997) showed how this method can be applied to different financial assets with emphasis on its use, like monetary policy. [START_REF] Abadir | Density-Embedding Functions[END_REF] determined farm formulas of option prices using Kummer functions (hypergeometric) as a basis for RND. Backus, Foresi, Li and Wu (1997) have approached the conditional distribution of underlying assets prices by Gram-Charlier series expansion.

Another approach has been proposed by [START_REF] Jarrow | Approximate valuation for arbitrary stochastic processes[END_REF]. They have developed a method of option valuation under the assumption that the underlying asset does not follow a log-normal distribution. They showed that the RND can be obtained by an Edgeworth expansion. Carrodo and Su (1996) used the Jarrow and Rudd approach to determine the skewness and kurtosis of the price of S&P500 options. A similar approach is that of [START_REF] Madan | Contingent claims valued and hedged by pricing and investing in a basis[END_REF], Abken, Madan and Ramamurti (1996) and [START_REF] Coutant | Implied risk aversion in option prices using hermite polynomials, in Estimating and interpreting probability density function[END_REF] who determined the RND from a Hermite polynomial approximation. Similarly, El Hassan and Kucera (1998) used the Fourier-Hermite development in order to evaluate European and American index options.

Among other recent researchs, [START_REF] Heston | A closed-from solution for options with stochastic volatility with applications to bond and currency options[END_REF] has provided a quasi-analytical solution for pricing options in a stochastic volatility framework. [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF] showed that the option price second derivative over the strike gives the risk-neutral density. Based on this relationship between option prices and RND, Rzepkawski (1996) adopted the Heston method to extract the RND. Similarly, Neuhauss (1995) used the direct relationship of [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF] working with the distribution function instead of the risk-neutral density. [START_REF] Rubinstein | Implied binomial trees[END_REF]Rubinstein ( , 1996) ) and Jackwerth (1996) suggested a method based on binomial trees. [START_REF] Rubinstein | Implied binomial trees[END_REF] developed a tree to estimate the implicit prices of contingent assets1 from prices options. His method is to minimize the gap between the tree implied probabilities and the probabilities determined from the tree of Cox-Ross- Rubinstein (1979). Aït-Sahalia (1996) and Aït-Sahalia and Lo (1998) estimated the S&P500 Risk-Neutral Density by performing Kernel techniques estimation. Campa, Chang and Reider (1998) compared three different RND estimating methods of the underlying asset. These methods are a smoothing of the smile through a degree three polynomial, a Rubinstein implicit tree and a mixture of log-normal distributions.

Galati and Melick (1999) estimated the moments using a mixture of log-normal laws in order to understand how the central bank interventions are perceived by foreign exchange markets traders regarding JPY / USD options between 1993 and 1996. Weinberger (2001) and Anagnar, Bedendo, Hodges and Tompkins (2002) found a typical form of the implied S&P500 risk neutral distribution after the crisis. [START_REF] Panigirtzoglou | A new approach to modelling the dynamics of implied distribution: Theory and Evidence from the S&P 500 Options[END_REF] studied implied distributions dynamics and provided algorithms that make their results applicable to the options and risk management framework. Mark, Feike and Lo (2005) showed that the RND is a function of both the underlying asset returns and volatility. Enzo, Handel and Härdle (2006) estimated the RND using the same approach developed by Aït-Sahalia and Lo (1998). Finally, Jondeau, Ser-Huang Poon and Rockinger (2007) compared the performance of various parametric and semi-parametric methods to extract RND.

Some authors argued that the RND coresponds to real world density only if investors are risk neutral. Thus, the difference between the two densities stemps from the existence of a risk premium. Consequently, it is interesting to study the ability of the estimated RND extracted from options prices to assess the real density, its forecasting ability as well as the choice of the suitable RND estimation technique.

In this paper, we estimate parametrically, semi-parametrically and nonparametrically the RND functions of CAC 40 index options between January 1 st 2007 and December 31 st 2007, and we compare the goodness-of-fit of eight option-based approaches during a normal and troubled period. To our knowledge, this is the first application dealing with European options using structural and non-structural approaches of RND functions estimation. Specially, we use nonparametric estimation methods of the RND, that is the kernel and the tree-based methods as well as six parametric and semi-parametric option-based approaches: (i) the numerical approximation of the RND based on the second derivative of option prices with respect to the strike price, as suggested by [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF], (ii) the mixture of log-normal distributions following Melick and Thomas (1997), (iii) the Edgeworth expansion around the log-normal distribution of [START_REF] Jarrow | Approximate valuation for arbitrary stochastic processes[END_REF], (iv) the Hermite polynomials, suggested by Madan and Mline (1994), (v) Heston's stochastic volatility model (1993) and (vi) the jump diffusion model following Bates (1991). To our best knowledge, our study is the first that combines the three approaches : parametric, semi-parametric and nonparametric methods of RND estimation. The paper is organized as follows. Section 2 discusses our methodological non-structural approaches of RND functions estimation. Section 3 presents the structural methods we use to extract the RND functions. Section 4 describes our data, and contains the estimation results as well as related comments regarding the comparison of the various RND estimation methods. Section 5 concludes. [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF] were the first to derive the RND using the following price of a call option formula :

Non-structural approaches of RND functions estimation

The Breeden and Litzenberger relation

C(S t , t) = e -rτ E * [max (S T -K, 0) |S t , t] (1) = e -rτ ∞ 0 max (S T -K, 0) q (S T |S t , t) dS T .( (1))
Where for date t and maturity date T, we denote C the call price, r is the risk-free interest rate, S is the underlying asset price, K is the strike price and q(.) is the undiscounted RND. Differentiating this equation with respect to the exercise price K yields the discounted cd f

∂C ∂K = -e -rτ ∞ K q(S T )dS T . (2)
and differentiating twice yields the discounted pd f 2 ∂ 2 C ∂K 2 K=S T = e -rτ q(S T ).(3) These equations show that the second derivative of the call price yields the discounted RND 3 . This suggests that a first method to extract RND is to approximate it numerically applying the finite difference approach to ( 3). Nevertheless, this method relies on the hypothesis that there exist traded option prices for many strikes. This is not likely to be the case in practice. Also, it has been shown that RNDs estimated in this manner are very unstable. In fact, differentiating twice exacerbates even tiny errors in the prices and may be difficult 4 . That is why it is necessary to extract RND using alternative methods that put more structure on the option prices.Before describing such methods, we briefly show how the parameters of these models are estimated.

RND parameters estimation

Suppose that we have to estimate a given model with respect to a set of parameters θ. Assume that for horizon τ, we have N τ c strike prices for which we have call options and N τ p strike prices for which we have put options. For date t and horizon τ, observed call and put options prices are denoted C t,τ,i i = 1, ..., N τ c and P t,τ,i i = 1, ..., N τ p . Theoretical call and put option prices, implied by the assumed model, are denoted C t (K, τ, θ) and P t (K, τ, θ), respectively, for strike price K and horizon τ. Then, the parameter vector θ ∈ Θ is typically estimated by non-linear least squares, by minimizing for each day and each maturity

min θ ∈ Θ N τ c i=1 w c i C t,τ,i -C t (K i , τ, θ) 2 + N τ p i=1 w p i P t,τ,i -P t (K i , τ, θ) 2 . ( 4
)
where w c i and w p i are weights associated with option i and Θ is the domain of θ. These weights could be given by a measure of liquidity of a given option.

Parametric methods

Black & Scholes Model

The price dynamics of an underlying asset of a Black & Scholes Model ( BSM) is a Geometric Brownian motion. Under the log-normal assumption, the volatility is constant across both the exercice price and horizons. The BSM underlying asset price is given by :

dS t = (µ + σ 2 2 )S t dt + σS t dW t (5)
2 [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF] suggest to evaluate the RND using this approximation

∂ 2 C ∂K 2 ≈ e -rτ C(K i+1 ) -2C(K i ) + C(K i-1 ) (∆K) 2 .
where ∆K = 50. See also Roncalli (1997).

where µ and σ represent respectively the instantaneous return assets expectations and volatility. W t is a Brownian motion under the real probability measure P.

When markets are complete, Harrison and Pliska (1981) show that there exists a risk-neutral transformation that leads to the following formula :

dS t = rS t dt + σS t dW Q t (6)
where W Q t is a Brownian motion under the risk-neutral probability measure Q. The call option price C BS maturing at date T = t + τ, with strike price K and dividend yield δ t,τ , is given by :

C BS (S t , K, τ, r, δ t,τ , σ) = e -rτ ∞ 0 max[S T -K, 0] f BS (S T )dS T = S t Φ(d 1 ) -Ke -rτ Φ(d 2 ). ( 7
)
Where Φ(.) is the standard cumulative normal distribution function and

d 1 = log(S t /K) + (r -δ t,τ + σ 2 2 )τ σ √ τ , d 2 = d 1 -σ √ τ. (8)
Consequently, the RND is a log-normal density with a mean m = (ln

(S t ) + (r -δ t,τ ) -σ 2 2 )τ and variance s 2 = σ 2 τ : f BS (S T ) = e rτ ∂ 2 C BS ∂K 2 K=S T = 1 S T √ 2πσ 2 τ exp - [log(S T ) -m] 2 2s 2 . ( 9
)
All Black & Scholes parameters , except the volatility, are directly observable. Nevertheless, the BSM hypothesis, according to which the returns are homoscedastic and normally distributed, does not fit with empirical evidence showing that densities are skewed and leptokurtik.

Mixture of log-normal distributions

Bahra (1996), Melick and Thomas (1997) and Södernlind (1997) are among the first to describe the RND as a mixture of distributions. For option pricing, the most well-known distribution studied in the litterature is the mixture of log-normal densities. The reason is that it appears as a trivial extension of the Black & Scholes model that involves the single log-normal density function.

l (S T , m, σ) = 1 S T √ 2Πs 2 exp       - 1 2 log (S T ) -m s 2       . (10)
A mixture of such densities yields

q (S T ; θ) = M i=1 α i l (S T , m i , s i ) . ( 11
)
where θ regroups all the unknown parameters α i , m i and s i f or i = 1, ..., M, and M denotes the number of mixtures describing the RND. Obviously, to guarantee that q is a density, we must have α i > 0 for all i = 1, ..., M, and α 1 + ... + α M = 1. In other words, q should be a convex combinaison of various log-normal densities.

The option price for such a mixture of log-normal distributions is, for a given strike K and time to maturity

τ = T -t C LN (K, θ) = e -rτ +∞ K (S T -K) q (S T ; θ) dS T = e -rτ +∞ K (S T -K) M i=1 α i l (S T , m i , s i ) dS T = e -rτ M i=1 α i +∞ K (S T -K) l (S T , m i , s i ) dS T . ( 12
)
where we define the volatility over the horizon of the option as s i = σ i √ τ, to simplify notations. The last equality is obtained by simply inverting the sum and integral operators. There are various ways to evaluate the integral. For instance, we have 5

+∞ K (S T -K) l S T , µ, s dS T = (E [S T |S T > K] -K) Pr [S T |S T > K] = exp m + 1 2 s 2 1 -Φ log(K) -m -s 2 s -K 1 -Φ log(K) -µ s (13)
Finally, the option price is given by

C LN (K; θ) = e -rτ M i=1 α i      exp µ i + 1 2 s 2 i       1 -Φ       log(K) -m i -s 2 i s i             -K 1 -Φ log(K) -m i s i      = e -rτ M i=1 α i exp µ i + 1 2 s 2 i Φ       -log(K) + m i + s 2 i s i       -e -rτ K M i=1 α i Φ -log(K) + m i s i (14)
Under the risk-neutral probability, we have to impose the martingale condition stating that the current price S t under the RND is equal to the expected discounted price of the underlying asset e -rτ E [S T ], so that 6

S t = e -rτ E [S T ] = e -rτ M j=1 α j exp m j + 1 2 s 2 j . (15)
5 This result builds on Johnson, Kotz and Balakrishman (1994) who indicate that if S follows a log-normal distribution with mean µ and variance σ 2 , then

E[S|S > K] = exp µ + 1 2 s 2 1 -Φ (U -σ) 1 -Φ (U)
where U = (log(K)-m) σ 6 In order to estimate the parameters α i , m i and s i , we use this minimization program

min m 1 , m 2 , s 1 , s 2 , α 1 , α 2 N i=1 C t,τ,i -C t (K i , τ, θ) 2 + N i=1 P t,τ,i -P t (K i , τ, θ) 2 + [e -rτ 2 j=1 α j exp m j + 1 2 s 2 j ] 2 .
As an alternative, we can directly use the BSM formula and set, in a similar spirit to what precedes,

C LN (K, θ) = M i=1 α i S t Φ d 1,i -Ke -rτ Φ d 2,i (16)
where d 1,i = log (S T /K) + m i τ + 1 2 s 2 i /s i , and d 2,i = d 1,is i , the martingale condition may be imposed as follows

S t = CL LN (0, θ) (17)
This condition just means that, for a strike price of 0, the option will always get exercised and, hence, at maturity we will always get the underlying asset. In practice, imposing K = 0 in the BSM formula is not a good idea because log (0) is not defined. It is possible to approach this limit case by setting K equal to some very small value. There are several advantages of using this technique. First,the implementation is straightforward, because it suffices to slightly modify the usual BSM formula by introducing the additional parameter µ. Second, the implementation via the BSM formula introduces the current price of the underlying asset, S t . The parameter µ, will therefore be of the magnitude of an annualized asset return. This means that its value can be easily bounded, which is a welcome feature for numerical purposes. Third, the time to maturity τ is explicitly taken into account. This means that µ i will be an annualized number which implies that parameters extracted for different maturities are directly comparable.

The mixture of distribution approach comes, however, at a cost. A first drawback in fitting a mixture of distributions is the symmetry between the densities. To illustrate this point, assume that a given RND can be correctly described by a mixture involving exactly two normal densities. Using obvious notations, we can write

q (S T , α, m 1 , m 2 , s 1 , s 2 ) = α l (S T , m 1 , s 1 ) + (1 -α) l (S T , m 2 , s 2 ) (18)
Obviously, the same density is obtained if we invert the two log-normal distributions, i.e, q (S T , 1 -α, m 1 , m 2 , s 2 , s 1 ). Clearly, the order of the parameters plays a critical role here. For an optimization program, this means, however, that several parameter vectors are associated with a same density. This in turn could yield numerically unstable programs where the optimizer cycles in an infinite loop.

Another difficulty lies in testing the number of distributions that are involved in the RND. One may be tempted to use a likelihood-ratio type test to check if the ith density should be included in the mixture by testing if α i = 0. However, if α i = 0, the parameters m i and s i associated with the ith density are unidentified. Stated differently m i and s i could take any value, because they would not play a role in the density q. Sometimes, such parameters are called down and specific tests need to be developed. This type of problems appears also in switching regressions. There is no obvious solution to it. Critical values may be obtained by simulations. In partice, one may add distributions up to the point where adding them yields to no further improvement.

There exist various solutions to help the optimizer converge in terms of MSE. In the event of two densities, M = 2, we start with a grid for α. Since it is known that 0 ≤ α ≤ 1, one may subdivide the interval [0, 1] into equally spaced points, where α i = i∆, for i = 0, ..., N. In practice, ∆ = 0.1 often yields very good results as Jondeau, Ser-Huang Poon and Rockinger (2007) find. Furthermore, to avoid the problem of symmetry mentioned above, it is possible to impose that the densities remain in a given order. One possibility that appears to give satsisfactory results is to impose that s 1 > s 2 .

The first density will then have a larger standard deviation than the second one. Similar extensions may be given for the case where M = 3. The α parameters may be taken over a simplex. 7

Semi-parametric methods

Edgeworth expansions

This method is developed by [START_REF] Jarrow | Approximate valuation for arbitrary stochastic processes[END_REF] for which a numerical application can be found in Corrado and Su (1996). The idea is to capture deviations from log-normality by an Edgeworth expansion of RND q(S T |θ) around the log-normal density. The construction of an Edgeworth expansion is conceptually similar to computing a Taylor expansion but applies to functions in general and to densities in particular. In a conventional Taylor expansion, the function is approximated by a simple polynomial around a given point. Here, the RND is approximated by an expansion around a log-normal distribution. A clear advantage of the expansion approach is that the approximation, by involving parameters that can be varied, allows generating more general functions. Flamouris and Giamouridis (2002) and Beber and Brandt (2003) studied also the Edgeworth expansions. Let Q be the cd f of a random variable S T and q its density. Let's define the characteristic function of S as φ Q (u) = e isu q(s)ds. If moments of S T exist up to order n, then there exist cumulants of the distribution Q, denoted κ Q, j and implicity defined by the expansion

log φ Q (u) = n-1 j=1 κ Q, j (iu) j j! + •(u n-1 ) (19)
Thus, if the characteristic function φ Q (.) is known, by taking an expansion of its logarithm around u = 0, it is possible to obtain the cumulants. We have these relationships between the cumulants and moments up to the fourth order

κ Q,1 = E [S T ] , κ Q,2 = V [S T ] , κ Q,3 = E (S T -E [S T ]) 3 , κ Q,4 = E (S T -E [S T ]) 4 -3V [S T ] 2 .
Jarrow and Rudd show that an Edgeworth expansion of the fourth order for the true probabiliy distribution Q around the log-normal cd f L can be written, after imposing that the first moment of the approximating density and the true density are equal, κ Q,1 = κ L,1 . The density q(s) is

q(s) = l(s)+ κ Q,2 -κ L,2 2! d 2 l(s) ds 2 - κ Q,3 -κ L,3 3! d 3 l(s) ds 3 + κ Q,4 -κ L,4 + 3 κ Q,2 -κ L,2 2 4! d 4 l(s) ds 4 + (s) (20)
where (s) captures terms neglected in the expansion. The various terms in the expansion correspond to adjustments of the variance, skewness, and kurtosis. The interpretation of this expansion is similar to a Taylor expansion. Jarrow and Rudd further show that, under the approximated density, the price of a European call option with strike K can be approximated as

C(Q) = e -rτ ∞ K (S T -K) q(S T )dS T ≈ e -rτ ∞ K (S T -K) l(S T )dS T + e -rτ κ Q,2 -κ L,2 2! ∞ K (S T -K) d 2 l(S T ) dS 2
T dS T 7 For further details, see Jondeau, Ser-Huang Poon and Rockinger (2007).

-e -rτ κ Q

,3 -κ L,3 3! ∞ K (S T -K) d 3 l(S T ) dS 3 T dS T +e -rτ κ Q,4 -κ L,4 + 3 κ Q,2 -κ L,2 2 4! ∞ K (S T -K) d 4 l(S T ) dS 4 T dS T .
(21) Notice that the first term is simply the Black-Scholes formula. In addition, the log-normal distribution has the following property

∞ K (S T -K) d j l(S T ) dS j T dS T = d j-2 l(S T ) dS j-2 T S=K , f or j ≥ 2.( 22
)
We deduce for the call option price

C (Q) ≈ C(L) + e -rτ κ Q,2 -κ L,2 2! l(K) -e -rτ κ Q,3 -κ L,3 3! dl(K) dS T + e -rτ κ Q,4 -κ L,4 + 3 κ Q,2 -κ L,2 2 4! d 2 l(K) dS 2 T . (23)
For the log-normal density, the first cumulants are given by

κ L,1 = S t e rτ , κ L,2 = κ L,1 ϑ 2 , κ L,3 = κ L,1 ϑ 3 3ϑ + ϑ 3 , κ L,4 = κ L,1 ϑ 4 16ϑ 2 + 15ϑ 4 + 6ϑ 6 + ϑ 8 . ( 24
)
where ϑ = e σ 2 τ -1

1/2
and where the first relation follows from risk-neutral valuation. Jarrow and Rudd suggest identifying the second moment by imposing κ L,2 = κ Q,2 . This argument is also justified on numerical grounds by Corrado and Su (1996) who notice that without this condition there exists a problem of multicolinearity between the second and the fourth moments. Rather than estimating the remaining cumulants, κ Q,3 and κ Q,4 , Corrado and Su (1996) estimate standarlized skewness and kurtosis ( written respectively γ Q,1 and γ Q,2 ), which are defined through the relationships :

γ Q,1 = κ Q,3 (κ Q,2 ) 3/2 = 3ϑ + ϑ 3 , γ Q,2 = κ Q,4 (κ Q,2 ) 2 = 16ϑ 2 + 15ϑ 4 + 6ϑ 6 + ϑ 8 . ( 25
)
These expressions also hold for the log-normal density. The skewness and kurtosis of the lognormal density can therefore be derived easily from the cumulants above.

With the assumption of equality of the second cumulants for the approximating and the true distributions, it follows that

C(Q) ≈ C(L) -e -rτ γ Q,1 -γ L,1 κ 3/2 L,2 3! dl(K) dS T + e -rτ γ Q,2 -γ L,2 κ 2 L,2 4! d 2 l(K) dS 2 T ( 26
)
Using this expression, it is easy to estimate with nonlinear least squares the implied volatility, σ 2 , skewness, γ Q,1 , and kurtosis, γ Q,2 . The expression of the RND can be obtained after twice differentiating [START_REF] Jondeau | Reading the smile: the message conveyed by methods which infer risk neutral densities[END_REF] with respect to K and then evaluation over

S T q(S T ) = l(S T ) -(γ Q,1 -γ L,1 ) κ 3/2 L,2 6
d 3 l(S T ) dS 3 T + (γ Q,2 -γ L,2 ) κ 2 L,2 24
d 4 l(S T ) dS 4 T . ( 27
)
where the partial derivatives can be computed iteratively using

dl(S T ) dS T = -1 + log(S T ) -m σ 2 τ l(S T ) S T . ( 28
)
d 2 l(S T ) dS 2 T = -2 + log(S T ) -m σ 2 τ 1 S T dl(S T ) dS T - 1 S 2 T σ 2 l(S T ). ( 29
)
d 3 l(S T ) dS 3 T = -3 + log(S T ) -m σ 2 τ 1 S T d 2 l(S T ) dS 2 T - 2 S 2 T σ 2 dl(S T ) dS T + 1 S 3 1 T σ 2 l(S T ), (30)
d 4 l(S T ) dS 4 T = -4 + log(S T ) -m σ 2 τ 1 S T d 3 l(S T ) dS 3 T - 3 S 2 T σ 2 d 2 l(S T ) dS 2 T + 3 S 3 T σ 2 dl(S T ) dS T - 2 S 4 T σ 2 l(S T ), (31)
where m = log(S t ) + (r -σ 2 /2)τ. The RND in the Edgeworth case will be a polynomial whose coefficients directly command the skewness and kurtosis of the RND. We also remark that the RND involves rather complicated terms with derivatives of the log-normal density.

Hermite polynomials

An alternative, yet similar, semi-parametric approach relies on an approximation of the Gaussian density based on Hermite polynomials. The theoretical foundation of this method is elaborated in Madan and Mline (1994) and applied in Abken, Madan, and Ramamurtie (1996) and Coutant, Jondeau, and Rockinger (2001).

For numerical reasons, Madan and Mline consider the map from the actual prices into the space generated by standardized log-returns z. Let's denote the volatility over the horizon of the option as s = σ √ τ,

S T = S t exp µτ - 1 2 s 2 + sz ⇒ z = log(S T /S t ) -(µτ -1 2 s 2 ) s (32)
The standarlized log return z has zero mean and unit variance.

If we focus on a call option, the payoff of such an option, as a function z, is

c(S t , K, µ, s, τ) = max S t exp µτ - 1 2 s 2 + sz -K, 0 = g(z), (33)
Hence the price of a call option may be written as

C(S t , K, µ, s, τ) = e -rτ +∞ 0 c(S t , K, µ, s, τ)q z (z)dz = e -rτ +∞ 0 g(z)q z (z)dz (34)
Here, the notation g(z) illustrates that the result holds for all sorts of payoffs. q z (.) denotes the risk-neutral density of z. We can go from the RND defined in z space to the RND defined in S T space using the following change of variable

q(S T )dS T = q z       log(S T /S t ) -(µτ -1 2 s 2 ) s       × 1 s × dS T (35)
Now, we assume that any payoff g(z) may be expressed as a function of basis elements. It turns out that a basis for the Gaussian space is given by Hermite polynomials. In other words, there exist real numbers a k such that any payoff can be written as

g(z) = ∞ k=0 a k h k (z), (36)
with a k = z g(z)h k (z)Φ(z)dz, where Φ(z) is the Gaussian distribution and h k (z) denotes Hermite polynomials normalized to unit variance, defined as

h 0 (z) = 1, h 1 (z) = z, h 2 (z) = 1 √ 2 (z 2 -1), h 3 (z) = 1 √ 6 (z 3 -3z), h 4 (z) = 1 √ 24 (z 4 -6z 2 + 3). ( 37
)
The coefficients a k can be interpreted as the covariance between the option payoff and the kth Hermite polynomial risk. Consequently, the price of the call option becomes

C(S t , K, µ, s, τ, r) = e -rτ ∞ k=0 a k z h k (z)q z dz. (38)
Madan and Milne (1994) assume that the RND q z (z) can be represented as the product of a change of measure density and a reference measure density

q z (z) = λ(z)Φ(z). (39)
Here, the reference measure density Φ(z) is simply the Gaussian one and the risk-neutral change of measure density λ(z) is approximated by a Hermite polynomial expansion

λ(z) = e rτ ∞ l=0 π l h l (z) (40)
Under the reference measure, z is by construction normally distributed with zero mean and unit variance. But more generally, the RND q z (z) may incorporate some departure from normality. The RND can be rewritten as

q z (z) = Φ(z)        ∞ l=0 π l h l (z)        (41)
The π l are interpreted as the implicit price of polynomial risk h l (z). Obviously, since the polynomial components are not traded, these risks can not be traded either. For practical purposes, the infinite sum can be truncated up to the fourth order. Since the Hermite polynomial of order l will depend on the lth moment, we will also refer to π 3 and π 4 as the price of skewness and kurtosis, respectively. It can be easily shown that for q z (z) to be a density, we need π 0 = e -rτ . Besides, since z is the standardized log-return (with zero mean and unit variance by construction), it follows that the shifts for the mean π 1 and the variance π 2 of q z (z) relative to the reference measure can be set equal to 0. Hence, we estimate the mean µ and the variance σ 2 of log-returns and set π 1 = π 2 = 0 . We obtain

q z (z) = Φ(z)        e rτ l=0 4 π l h l (z)        = Φ(z)e rτ e -rτ + π 3 h 3 (z) + π 4 h 4 (z) = Φ(z) 1 + b 3 √ 6 (z 3 -3z) + b 4 √ 24 (z 4 -6z 2 + 3) (42)
Where the b i = e rτ π i , i = 3, 4, are the future value of the ith price of risk coefficient. The parameters b 3 and b 4 correspond to the skewness and kurtosis if the reference measure density of z is chosen to be the normal distribution. It is important to emphasize that unlike the Edgeworth case, since a further change of variable from z to S T has to made, b 3 and b 4 will not correspond to the skewness and kurtosis of the underlying price process S T but to the skewness and kurtosis of log-returns z.

Finally, the skewness and kurtosis of the expansion are given by

Skewness [z] = √ 6b 3 ( 43
)
Kurtosis [z] = 3 + √ 24b 4 ( 44
)
The parameters σ , b 3 and b 4 can be obtained using a nonlinear estimation method. The general expression for the price of a call option is given by

C(S t , K, µ, s, τ, r) = ∞ k=0 ∞ l=0 a k π l z h l (z)h k (z)Φ(z)dz = ∞ k=0 a k π k (45)
where the second equality holds because Hermite polynomials form an orthogonal system . When the infinite sum is truncated up to the fouth order, we obtain

C(S t , K, µ, s, τ, r) = e -rτ a 0 + π 3 a 3 + π 4 a 4 (46)
Therefore, it remains to be shown how to obtain the a k coefficients. Abken, Madan and Ramamurtie (1998) introduce the call option generating function

G(u, S t , K, µ, s, τ) = 1 √ 2Π +∞ -∞ c(u, S t , K, µ, s, τ)exp - 1 2 (z -u) 2 dz, ( 47
)
where u is a dummy variable. The coefficients a k can then be obtained from the generating function as

a k (S t , K, µ, s, τ) = ∂ k G(u, S t , K, µ, s, τ) ∂u k | u=0 1 √ k! (48)
The evaluation of G(u, S t , K, µ, s, τ) is shown to be

G(u, S t , K, µ, s, τ) = S t exp µτ + su Φ(d 1 (u)) -KΦ(d 2 (u)) (49)
where d 1 (u) = log(S t /K) + (µτ/s + s/2) + u and d 2 (u) = d 1 (u)s. From there on, it is possible to compute the coefficients a k . Let's introduce the notation d j = d j(0) , that is we evaluate the d j(0) at zero. We obtain after tedious computations, and introducing the notations Φ , Φ andΦ , for the first, second, and third derivatives of standard normal density with respect to its argument

a 0 = S t e µτ Φ(d 1 ) -KΦ(d 2 )
a 1 = sS t e µτ Φ(d 1 ) + S t e µτ Φ(d 1 ) -KΦ(d 2 ) (50) This method has been applied in [START_REF] Coutant | Reading PIBOR Futures Options Smile : The 1997 Snap Election[END_REF] to interest-rate options. [START_REF] Coutant | Reading PIBOR Futures Options Smile : The 1997 Snap Election[END_REF] have also shown how conditions on b 3 and b 4 could be imposed so that the polynomial approximation remains positive.

a 2 = 1 √ 2 s 2 S

Non-parametric methods

Tree-based methods

The tree-based method has been presented by [START_REF] Rubinstein | Implied binomial trees[END_REF] and is further extended in Jackwerth and Rubinstein (1996) and Jackwerth (1999). Jackwerth (1997) gives a numerical illustration of the method which is based on two steps. First of all, we take existing options nearest as well as we compute the average of their Black-Scholes implied volatilities. Secondly, we compute the RNDs that would be associated to the binomial tree of Cox, Ross and Rubinstein (1979) and that would be comparable with these implied volatilities. Formally, let's consider that the tree has N steps. A terminal node, j, of the binomial tree has a probability to realize of P j . If p is the risk-neutral probability of an up movement on each node of the tree 8 , then In that formula, it is assumed that the nodes are denoted in such a way that a small j corresponds to a low value of the support S T . Let's denote these values as S j , j = 1, ...N. Given the current value of the underlying asset, it is also possible to construct the nodal values of the underlying asset in a risk-neutral environment. Even though by working backwards, we could use the binomial tree to value options, this tree would not be compatible with the observed option smile. The third possibility is to seek terminal probabilities that are close to the tree -based probabilities, yet that also correspond to the given options. Formally, let the risk-neutral probabilities, which are compatible with actual options, be denoted by P j . These risk -neutral probabilities may then be obtained by solving the following optimization : min

P j = N! j!(N -j)! p j (1 -p ) (N-j) (51)
P j N j=1 P j -P j 2 (52) subject to N j=1 P j = 1, P j ≥ 0, j = 0, ..., N (53)
S t = e -r t (τ) N j=1 S j P j (54)

C τ,i = e -rτ N j=1 P j S j -K i + , i = 1, ..., N (55)
The first constraint provides conditions so that the values P j define probabilities. The second condition corresponds to the martingale condition. Finally, the last condition states that the probabilities P j should be compatible with existing option prices. Various metrics can be used to define the distance between the probabilities P j and P j . The idea of using other metrics is expressed in [START_REF] Rubinstein | Implied binomial trees[END_REF] and tested in Jackwerth and Rubinstein (1996). We could choose the following distances : min

P j N j=1 P j -P j , (56)
min

P j N j=1           P j -P j 2 P j           (57)
or max

P j N j=1 P j log        P j P j        (58)
The first measure is a metric based on absolute deviations. The second measure is based on percentage deviations. The third measure corresponds to seeking the maximum entropy.

Kernel regression

Aït-Sahalia and Lo (1998) propose to estimate the RND nonparametrically by exploiting Breeden and Litzenberg's (1978) insight that f * t (S T ) = exp(r t,τ )∂ 2 H(.)/∂K 2 . They suggest using market prices to estimate an option pricing formula Ĥ(.) nonparametrically, which can then be differentiated twice with respect to K to obtain ∂ 2 Ĥ(.)/∂K 2 . They use kernel regression to construct Ĥ(.) and assume that the option-pricing formula H to be estimated is an arbitrary nonlinear function of a vector of option characteristics y ≡ F t,τ K τ r t,τ where F t,τ is the forward price of the asset.

In practice, they propose to reduce the dimension of the kernel regression by using a semiparametric approach. Assume that the call pricing function is given by the parametric Black-Scholes model except that the implied volatility parameter for that option is a nonparametric function σ(K/F t,τ ) :

C(S t , K, τ, r t,τ , δ t,τ ) = C BS (F t,τ , K, τ, r t,τ , σ(K/F t,τ )). ( 59
)
We assume that the function C defined by the equation ( 59) satisfies all required conditions to be a "rational" option-pricing formula in the sense of Merton (1973, 1992) 9 . In this model, we only need to estimate non-parametrically the regression of σ on a subset Ỹ of the vector of explanatory variables Y. The rest of the call pricing functions C(.) is parametric, thereby considerably reducing the sample size n required to achieve the same degree of accuracy as the full nonparametric estimator. Aït-Sahalia and Lo (2000) partition the vector of explanatory variables y = [ Ỹ F t,τ r t,τ ] where Ỹ contains d non parametric regressors. In their empirical application, they consider Ỹ ≡ [K/F t,τ τ] ( d = 2) and form the Nadaraya-Watson Kernel estimator of E[σ|K/F t,τ τ] as :

σ K/F t,τ τ = n i=1 k K/F K/F t,τ -K i /F t i ,τ i h K/F k τ τ-τ i h τ σ i n i=1 k K/F K/F t,τ -K i /F t i ,τ i h K/F k τ τ-τ i h τ (60)
where σ i is the volatility implied by the option price H i , and the univariate kernel functions k K/F and k τ and the bandwidth parameters h K/F and h τ are chosen to optimize the asymptotic properties of the second derivative of C(.), i.e., of the RND estimator. One estimates the call pricing function as :

Ĉ(S t , K, τ, r t,τ , δ t,τ ) = C BS (F t,τ , K, τ, r t,τ , σ(K/F t,τ )) (61)
The RND estimator follows by taking the second partial derivatives of C(.) with respect to K :

f * (S T ) = e r i,τ τ ∂ 2 Ĉ(S t , K, τ, r t,τ , δ t,τ ) ∂K 2 K=S τ .( 62
)
Contrary to the RND estimation approaches previously developed and which are based on the estimation of RND for each cross-section of options, the nonparametric kernel regression provides an estimator that is based on both cross-sectional and time-series option prices. Consequently, this approach is consistent over time, but may provide a poor fit for certain dates 10 .

Structural approaches of RND functions estimation 3.1 Jump diffusion model

Modeling the underlying asset with a Black & Scholes model is not consistent. Large values of returns occur too frequently to be consistent with the normality assumption. The rare events may cause brutal assets pricing variations. In order to model such a phenomenon, we should use a Poisson jump process. We assume that the underlying asset price follows a log normal jump diffusion process i.e, the addition of a geometric Brownian motion and a Poisson jump process. This process takes into account the skewness and kurtosis effects.

The underlying asset price under the risk-neutral world follows the process :

dS t = (r t -λ E(k)) S t dt + kS t dW t + kS t dq t ( 63
)
where q t is the Poisson associated probability , λ is the average rate of jump occurence and k is the size of a jump. Bates (1991) has shown that in case of a diffusion with n jumps, the call price is :

C JMP (S t , t, T, K) = e -r t τ ∞ i=0 P(n jumps)E (S T -K) + /n jumps (64)
Ball and Torous (1985) and Malz (1997) proposed a simplified version of the jump-diffusion model, where there will be at most one jump of a constant size. In this Bernoulli version of a jump diffusion model, the equation of the call price is the following :

C JMP (S t , t, T, K) = e -r t τ 1 i=0 P(n jumps)E (S T -K) + /n jumps (1 -λ τ) S t 1 + λκτ N d 1 + σ √ τ -Ke -r t τ N (d 1 )
+ λτ S t 1 + λκτ (1 + κ)N d 2 + σ √ τ -Ke -r t N(d 2 ) (65)
with

d 1 = log S t K -log(1 + λκτ) + r t -σ 2 2 τ σ √ τ ( 66
)
d 2 = log S t K -log(1 + λκτ) + log(1 + κ) + r t -σ 2 2 τ σ √ τ (67)
where (1 -λτ) is the probability that no events occur in the option life. We can easily show that the call price is simply a combinaison of two call prices calculated by [START_REF] Black | The Pricing of Commodity Contracts[END_REF] model

C JMP (S t , t, T, K|θ) = (1 -λτ) C B S t 1 + λκτ , τ, K|θ + (λτ) C B S t 1 + λκτ (1 + κ), τ, K|θ (68)
This technique is a particular case of a mixture of log-normal distributions. The RND is :

q JMP (S T , S t , t, T, K|θ) = (1 -λτ) LN S T , α, β + λτ LN S T , α + log(1 + κ), β (69)
where :

LN (x) = 1 √ 2πβ x exp - ( Ln (x) -α) 2 2β 2 (70) α = log (S t ) + (µ - 1 2 σ) 2 (T -t) (71)
the θ = (σ, λ, κ) vector parameters are estimated using the following program :

min θ n i=1 C i -C JMP i 2 (72)

Heston's stochastic volatility model

Most diffusion models of the underlying asset suppose that the volatility of the latter is a constant. Nevertheless, if the underlying asset records strong fluctuations over a short period, it would be necessary to model it with a stochastic volatility process . [START_REF] Heston | A closed-from solution for options with stochastic volatility with applications to bond and currency options[END_REF] proposed the following model :

dS t = µ t S t dt + V t S t dW 1 t ( 73
)
dV t = κ (θ -V t ) dt + σ V t dW 2 t ( 74
)
dW 1 t dW 2 t = ρ dt ( 75
)
where {S t } | t≥0 and {V t } | t≥0 are the price and volatility process, respectively, and W 1 t t≥0 , W 2 t | t≥0 are correlated Brownian motion processes (with correlation parameter ρ ). {V t } | t≥0 is a square root mean reverting process, first used by Cox, Ingersoll and Ross (1985), with long-run mean θ, and rate of reversion κ. σ is referred to as the volatility of volatility. All the parameters µ, κ, θ, σ, ρ are time and state homogenous. There are several economic, empirical and mathematical reasons for choosing a model with such a form 11 .

Under the Heston model , the value of any option Θ(S t , V t , t, T) must satisfy the following partial differential equation :

1 2 VS 2 ∂ 2 Θ ∂S 2 + ρσ VS ∂ 2 Θ ∂S∂V + 1 2 σ 2 V ∂ 2 Θ ∂V 2 + rS ∂Θ ∂S + κ [θ -V] -Λ (S, V, t) σ √ V ∂Θ ∂V -rΘ + ∂Θ ∂t = 0 (76)
where Λ(S, V, t) is called the market price of volatility risk. Heston supposes that the market price of volatility risk is propotional to volatility, i.e.

Λ(S, V, t) = k √ V (77) ⇒ Λ(S, V, t)σ √ V = kσ V t for some constant k = λ(S, V, t) (78)
Risk neutral valuation is pricing of a contingent claim in an equivalent martingale measure. The price is evaluated as the expected discounted payoff of the contingent claim, under the Risk neutral probability, the option value is

Option value = E Q e r(T-t) H(T) (79)
where H(T) is the payoff of the option at time T and r is the risk free of interest over [t, T].

Moving from a real world measure P to a risk neutral world measure Q is achieved by Girsanov's theorem. In particular, we have

d W 2 t = dW 2 t + Λ(S, V, t)dt ( 80
)
dQ dP = exp - 1 2 t 0 ν 2 s + Λ(S, V, s) 2 ds - t 0 υ s dW 1 s - t 0 Λ(S, V, s)dW 2 s ( 81
)
υ t = µ -r √ V t ( 82
)
where W1 t | t≥0 and W2 t | t≥0 are Q-Brownian motions. Under measure Q, the equations ( 73), ( 74) and ( 75) become,

dS t = r t S t dt + V t S t d W 1 t ( 83
)
dS t = κ * (θ * -V t ) dt + σ V t d W 2 t (84) d W 1 t d W 2 t = ρ dt (85)
and the equation ( 76) becomes :

1 2 VS 2 ∂ 2 Θ ∂S 2 +ρσ VS ∂ 2 Θ ∂S∂V + 1 2 σ 2 V ∂ 2 Θ ∂V 2 +rS ∂Θ ∂S + κ * [θ * -V] -Λ (S, V, t) σ √ V ∂Θ ∂V -rΘ+ ∂Θ ∂t = 0 (86)
where,

κ * = κ + λ (87) θ * = κθ λ + λ (88)
In this important result, λ has effectively been eliminated. The closed-form solution of a European call option on a non-dividend paying asset for the Heston model is :

C S t , V t , t, T) = S t P 1 -Ke -r t (T-t) P 2 (89)
where

P j (x, V t , K, T) = 1 2 + 1 π ∞ 0 Re         e -iφ log (K) f j x, V t , T, φ iφ         dφ (90) f j x, V t , T, φ = exp C T, φ + D T, φ V t + iφ x (91) C T -t, φ = rφ iτ + α σ 2       b j -ρσφ i + d j τ -2 log       1 -g j e d j τ 1 -g j             (92) D(T, φ) = b j -ρσφ i + d σ 2       1 -e d j τ 1 -g j e d j τ       (93) g j = b j -ρσφ i + d b j -ρσφ i -d ( 94
)
d j = (ρσφ -b j ) 2 -(2u j φ i -φ 2 ) (95) x = log (S t ) (96)
for j = 1, 2 where,

u 1 = 1 2 , u 2 = -1 2 , α = κθ, b 1 = k * -ρσ, b 2 = κ * 12 .
The θ = κ * , θ * , ρ, V t , σ vector parameters are then estimated using the following algorithm :

min θ n i=1 C i -C i 2 (97)
After estimating model parameters, the risk-neutral density is obtained from the numerical calculation of the integral. A method to evaluate formulas in form of ( 92) has been proposed by Carr and Madan (1999).

Application to CAC 40 index options

In this section, we estimate parametrically, semi-parametrically and nonparametrically the RND functions of CAC 40 index options between January 1 st 2007 and December 31 st 2007. We compare the goodness-of-fit of eight option-based approaches during a normal and troubled period. To our knowledge, this is the first application on the European market dealing with European options using structural and non-structural approaches of RND functions estimation.

Data

Our data is provided by the SBF-Paris Bourse 13 and includes intraday values of the CAC 40 stock index and intraday transaction prices of CAC 40 options over the period January 1 st 2007-December 31 st 2007. CAC 40 options are traded on the MONEP, the French derivatives market.

Trading covers eight open maturities : three spot months, four quarterly maturities ( March, June, September, December) and two half-yearly maturities (March, September). The maturity date is the last trading day of each month. Trading takes place on a continuous basis between 9:00 am and 5:30 pm. The future underlying asset F t,T cannot be observed exactly at time t, where t denotes a given day and time and T is the maturity . Aït-Sahlia and Lo (1998) suggest extracting an implied future index price from the put-call parity. Given our intraday data, we cannot get contemporaneous trades concerning a call and a put with similar strike price and maturity. Thus, for each maturity, we compute :

F t,T = S t (exp(r t,T -δ t,T )τ) (98)
To circumvent the unobservability of the dividend rate δ t,T , we extract an implied value of the dividend rate between the end of day t and T from the daily closing price of the index, S t end and the setlement future index price, F t end , observed at the end of each day. The obtained dividend rate is the dividend rate expected by the market between t end and T14 .

δ t end ,T = r t,T - 1 τ ln( F t end ,T S t end ) (99)
F t,T is then computed using the dividend rate and the riskfree interest rate proxied by Euribor. Since some CAC 40 options are not traded actively, we need to filter the data carefully. Five filters are applied to the initial data. We omit the quotes of the first and the last 15 minutes each day and the option quotes characterized by a price lower or equal to one tick. We only consider options with a moneyness 15 comprised between 0.85 and 1.15. This procedure eliminates far-away-from the money observations, which are unreliable due to their low volume and low sensitivity towards volatility. Besides, options quotes violating general no-arbitrage conditions, that is put-call parity, are eliminated as well as we replace the price of all illiquid options, that is in-the-money options, with the price implied by put-call parity at the relevant strike prices. Specially, we replace the price of each in-the-money call option with P(S t , K, τ, r t,T , δ t,T ) + (F t,T -K)e -r t,T τ , where, by construction, the put with price P(S t , K, τ, r t,T , δ t,T ) is out-of-the money and therefore liquid. After this procedure, all the information included in liquid put prices is extracted and resides in corresponding call prices through the put-call parity. Put prices may now be splayed without any loss of reliable information. Since our data include only call options, a single volatility function can be estimated. We provide in Figure 2 the implied volatility observed on January 1 st 2007 and Figure 3 reports the volatility smiles for four various maturities observed on July 10 th 2007, when the market trend is bear. The implied volatility is computed using the Black & Scholes formula. For a given maturity, we compute the daily average implied volatility for each strike price. For fixed S t and T, it is shown that the implied volatility function is not constant as assumed in the Black & Scholes model but decreases with respect to the strike price. The implied volatility observed on July 10 th 2007 on the French market is then more a "skew" than a "smile". In this case, contrary to what suggests Black-Scholes model, out-the-money and in-the-money options represent increased risk on potentially very large movements in the underlying compared to at-the-money. To compensate for this risk, they tend to be priced higher.

Kernel approach estimation

We focus on the nonparametric Kernel approach proposed by Aït-Sahlia and Lo (2001). This method does not require any parametric restriction on either the underlying asset price dynamics, the family of distributions that the RND belongs to or any prior distribution for the RND. We use the whole database in order to estimate nonparametrically the volatility function given in paragraph (2.5.2), with two conditional variables (d=2): the moneyness and the time to maturity τ. We only consider two regressors. We rely on the Aït-Sahlia and Lo (1998)'s procedure to select the kernel functions and the optimal bandwidths. They are chosen to optimize the asymptotic properties of the option pricing function, H. The selected kernels associated with the two regressors are Gaussian :

k K/F = 1 √ 2π e -0.5.(K/F) 2 k τ = 1 √ 2π e -0.5.(τ) 2 (101)
For each regressor, we choose the optimal bandwidth h according to the following equations :

h K/F = c K/F σ K/F n -1/(d+2(q K/F +m)) (102)
h τ = c τ σ τ n -1/(d+2q τ ) ( 103
)
where n is the sample size, q j is the number of existing continuous partial derivatives of the function to be estimated with respect to the j th regressor (j = K/F, τ), m is the order of the partial derivative with respect to the j th regressor that we want to estimate, d is the number of regressors and σ j is the unconditional standard deviation of the j th regressor. The parameter c j depends on the choice of the kernel and the function to be estimated. It is typically of the order of one. Aït-Sahlia and Lo (1998) emphasize that small deviations from the exact value have no noticeable effects. We report in Table 2 the values of the coefficients and bandwidths.

Table 2: Selected kernels and bandwidth values

Kernel n q j m d σ j h K/F Gaussian 25868 2 2 2 0.0724 0.0262 τ Gaussian 25868 2 0 2 100.5462 18.4886

Figure 4 shows σ K/F t,τ τ , the Nadaraya-Watson estimator of the volatility parameter with two regressors X/F and τ, for different maturities (20 days , 50 days, 80 days and 170 days i.e 1,2,3 and 6 months). We note that the estimator generates a strong volatility smile with respect to the moneyness. Besides, in confirmity with Aït-Sahlia and Lo (1998) and [START_REF] Bliss | Testing the Stability of Implied Probability Density Functions[END_REF], we observe that the reported smiles shape varies with the chosen time to maturity.

The estimated RND is shown in Figure 5 for 3-month maturity on January 10 th 2007. The RND is obtained by computing the second derivative of the call price with respect to the strike price. We note from Figure 5 the limits of the kernel method. Density is unclear because of the low number of observed points. This method would probably be better suited to a study of a large sample or to estimate the returns time series.

Tree-based methods

The estimated RND on a discrete probability framework, as proposed by Jackwerth and Rubinstein (1996), is shown in Figure 6 for 3 periods on January 10 th 2007 (normal period), on July 10 th 2007 on which the market trend is bear and on October 17 th 2007 for various maturities. We note from Figure 6 that in the troubled period, the RNDs for all the maturities have heavier tails than just before and after the crisis periods. The CAC 40 index options exhibit heavy tail distribution in the troubled period; so in a risk neutral framework, the probability estimated by the tree based method does not undervalue the probability of occurrence of extreme events.

Parametric and semi-parametric RND estimation methods comparison

Several papers compare the performance of various methods to extract RND. Jondeau and Rockinger (2000) compare the mixture of log-normals, the Edgeworth and Gram-Charlier approximation, the single-jump and Heston's stochastic volatility model using exchange rates. They find that up to a certain extent, these methods provide similar RNDs. They remark that during normal periods, a mixture of log-normals provides rather good results. However during a troubled period, the jump diffusion model performs better. [START_REF] Coutant | Reading PIBOR Futures Options Smile : The 1997 Snap Election[END_REF] compare the mixture of log-normals, Gram-Charlier expansion and the entropy-based method regarding speed, estimation robustness and ease of implementation using interest rates. They find that for interest rates, the method based on Gram-Charlier expansion seems to provide rather stable results. Bliss and Panigrizoglou (2002) compare the robustness of a model based on the mixture of two log-normal distributions and a smoothed implied volatility smile model, in the spirit of [START_REF] Shimko | Bounds of probability[END_REF]. They find that their smoothed spline method performs better than the mixture of log-normal distributions. However, despite its high flexibility, the smoothed spline method does not recover the tails of the RND outside the range of available strike prices.

In addition to nonparametric estimation methods of the RND, i.e. the kernel and the tree-based methods, we now compare six parametric and semi-parametric option-based approaches used to extract RND from CAC 40 index options during a normal and troubled period using the same data as in section 4.1. The first technique is the numerical approximation of the RND based on the second derivative of option prices with respect to the strike price, as suggested by [START_REF] Breeden | Prices of state-contingent claims implicit in option prices[END_REF]. The second method is obtained using a mixture of two log-normal distributions, following Melick and Thomas (1997). The third RND technique is the Edgeworth expansion around the log-normal distribution of [START_REF] Jarrow | Approximate valuation for arbitrary stochastic processes[END_REF]. The fourth method is Hermite polynomials, suggested by Madan and Mline (1994). The fifth RND is based on Heston's stochastic volatility model (1993). Finally, we consider the jump diffusion model following Bates (1991). The results of estimated parameters of different RND models are reported in Tables 3 to 8 16 . The various resulting RNDs are plotted in Figures 7 to 12 17 . The mixture of log-normals and Edgeworth expansion, Hermite polynomials and jump-diffusion models are more in line with historical data than the numerical approximation of the RND using the Breeden and Litzenberg formula. For this reason, the nonparametric approaches such as kernel and tree-based model ,as described above, were developed. The mixture of log-normals, Edgeworth expansion, hermite polynomials, jump diffusion and Heston models have heavier tails than the log-normal distribution.

The existence of several RND estimation methods from option prices raises the question of the suitable method to be chosen. To answer this question, we provide goodness-of-fit measures, allowing ton investigate how well theoretical options prices calculated with RND estimation approaches fit observed market prices Ĉ. Synthetic measures for errors calculating are used : the Mean Squared Error (MSE) and the Average Relative Error (ARE).

MSE = 10 2 m -n m i=1 (C i -Ĉi ) 2 (
104) 16 See Appendix B. 17 See Appendix A.

ARE

= 10 4 m -n m i=1 ( C i -Ĉi C i ) 2 (105)
where m is the number of observed options prices and n is the number of method parameters. The best method is the one which has the smallest errors.

In order to compare the various RND models, not only we check the statistical properties of the different RNDs but we also calculate the MSE and the ARE for different dates 01-10-2007, 07-10-2007 (troubled period) and 10-17-2007 and various maturities 1,2 and 3 months. All results are displayed in Tables 9 to 13 18 . First of all, we compare the volatilities, skewness and kurtosis obtained under different RNDs. Concerning the volatilities (Table 9), the volatility induced by the Black & Scholes model is smaller than the one we obtain using the other methods, and this is true for all periods and all maturities, 19 which explain the bias implied by the log-normality assumption. Regarding the implicit skewness (Table 10) and the implicit kurtosis (Table 11), the log-normal model is less interesting because it does not allow for asymetry or fat tails. We observe a strongly negative skewness and a larger excess of kurtosis for all methods but the Black & Scholes model. The skewness obtained by the semi-parametric methods ( Hermite polynomial and Edgeworth expansion models) are lower than that of the other models on the normal period (01-10-2007). Whereas on the troubled period (07-07-2007), Hermite polynomials and the mixture of log-normals models have the lower skewness, on the period just after the crisis, the jump-diffusion model has the lower skewness from the two month maturities, followed by the Hermite polynomials model. The graphs of the RND corroborate these findings.

Table 12 reports the Mean Squared Error (MSE) for the different methods. According to this criterion, the jump diffusion model on the period just before the crisis for relatively short maturities (1 and 2 months) provides a much better fit than the other models. However, during this same period, the mixture of log-normal models provides the best fit of the data for the 3 month maturity. We notice that in the troubled period (07-10-2007), the Edgeworth expansion model provides the best fit for all maturities. Concerning the period just after the crisis, the jump diffusion model performs well for the two month maturity while the Edgeworth expansion model provides the best fit for the three month maturity.

The Average Relative Error (ARE) results are presented in Table 13. As for the MSE criterion, the jump diffusion model outperforms the other models on the period just before the crisis. On the troubled period, Hermite polynomial method and the mixture of log-normal models have the smallest errors. The difference between the two corresponding ARE is small. Concerning the period just after the crisis, the semi-parametric methods provide the best fit for the one and three month maturities, while the mixture of log-normals approach provides a much better fit than the other models.

The RNDs functions are of paramount importance. Once estimated, it is possible to extract a lot of various information. For instance, we can perform tests and compute confidence intervals around the expected value whose evolution offers the opportunity to investors to measure how markets are thought to evolve through time. Besides, RNDs can provide a measure of expected extreme variations in the underlying asset price, which is an important tool for risk management. Moreover, since the existence of a possibly time-varying risk premium RNDs differs from real probability distributions, the RNDs incorpate information concerning the risk aversion of investors.

Conclusion

In this paper, eight option-based approaches are presented and applied to the French derivatives market. Structural and non-structural approaches are used to estimate risk-neutral probability density functions from a high-frequency CAC 40 index options during a normal and troubled period. We use Black & Scholes model, mixture of log-normals, Edgeworth expansions, Hermite polynomials, tree-based methods, kernel regression, Heston's stochastic volatility model and jump diffusion model. We determine nonparametric, semi-parametric, parametric and stochastic implied volatility functions and Risk-Neutral Densities (RND). In order to compare the various RND models, not only we check the statistical properties of the different RNDs, but also we calculate goodness of fit measures. We find that the kernel estimator generates a strong volatility smile with respect to the moneyness, and kernel smiles shape varies with the chosen time to maturity. The mixture of log-normals, Edgeworth expansion, hermite polynomials, jump diffusion and Heston models are more in line with historical data and have heavier tails than the log-normal distribution. Moreover, according to the goodness of fit criteria, the jump diffusion model provides a much better fit than the other models on the period just-before the crisis for relatively short maturities. However, during this same period the mixture of log-normals model performs better for more than three month maturity. Furthermore, in the troubled period and the period just-after the crisis, we find that semi-parametric models are the methods with the best accuracy in fitting observed option prices for all maturities with a minimal difference towards the mixture of log-normals model.

Appendix

A.Graphs

8 p

 8 = r-d u-d , where u = e σ √ τ and d = 1u

Figure 1 :

 1 Figure 1: CAC 40 stock index evolution (January 1 st 2007-December 31 st 2007)

Figure 2 :Figure 3 : 22 4. 2

 23222 Figure 2: Implied volatility surface on the French market on January 1 st 2007

Figure 4 :Figure 5 :

 45 Figure 4: Nonparametric estimators of the volatility for different maturities

Figure 6 :

 6 Figure 6: Risk Neutral Densities of the tree based method for different dates and various maturities

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 7891011 Figure 7: The Black & Scholes Risk Neutral Densities for different dates and various maturities

  t e µτ Φ(d 1 ) + 2sS t e µτ Φ(d 1 ) + S t e µτ Φ (d 1 ) -KΦ (d 2 )

	a 3 =	1 √ 6	s 3 S

t e µτ Φ(d 1 ) + 3s 2 S t e µτ Φ(d 1 ) + 3sS t e µτ Φ (d 1 ) + S t e µτ Φ " (d 1 ) -KΦ " (d 2 )

a 4 = 1 √

24

s 4 S t e µτ Φ(d 1 ) + 4s 3 S t e µτ Φ(d 1 ) + 6s 2 S t e µτ Φ (d 1 ) + 4sS t e µτ Φ " (d 1 ) + S t e µτ Φ (d 1 ) -KΦ (d 2 )

Table 3 : Estimation of Black & Scholes parameters

 3

	Dates	τ	m	s
		20 8.6153 0.0211
	01-10-2007	50 8.6185 0.0442
		80 8.6224 0.0600
		20 8.7050 0.0253
	07-10-2007	50 8.7086 0.0485
		80 8.7127 0.0678
		13 8.6701 0.0122
	10-17-2007	43 8.6739 0.0423
		73 8.6839 0.0590

Table 4 : Estimation of Mixture of log-normal parameters

 4

	Dates	τ	m 1	m 2	s 1	s 2	α
		20 8.6092 8.6152 0.0332 0.0200 0.2876
	01-10-2007	50 8.5705 8.6319 0.0671 0.0325 0.2713
		80 8.5361 8.6465 0.0835 0.0406 0.2729
		20 8.6928 8.7079 0.0407 0.0202 0.2894
	07-10-2007	50 8.6496 8.7257 0.0681 0.0328 0.2672
		80 8.6487 8.7478 0.0835 0.0405 0.4042
		13 8.6541 8.6741 0.0154 0.0111 0.3019
	10-17-2007	43 8.6043 8.6892 0.0623 0.0311 0.2515
		73 8.5620 8.7089 0.0796 0.0393 0.2729

Table 5 : Estimation of Hermite polynomials parameters

 5

	Dates	τ	σ	b 3	b 4
		20 0.0967 -0.4280 0.4845
	01-10-2007	50 0.1298 -0.4302 0.5250
		80 0.1435 -0.4283 0.4960
		20 0.1113 -0.2324 0.1289
	07-10-2007	50 0.1405 -0.4295 0.5314
		80 0.1612 -0.4320 0.4823
		13 0.0635 -0.0699 -0.0284
	10-17-2007	43 0.1284 -0.3848 0.3246
		73 0.1582 -0.4263 0.4633

Table 12 : MSE

 12 Note that τ is the maturity, B&S is the Black & Scholes model, Hermite is the Hermite model, Edgeworth is the Edgeworth model, ML is the mixture of log-normals model, Jump is the jump-diffusion model and Heston is the Heston's stochastic volatility model.

	Dates	τ	B&S	Hermite	Edgeworth ML	Jump	Heston
		20 167.902	44.33	47.1335	1419.4	17.32	2664.2
	01-10-2007	50 2392.39	1196.2	190.317	190.65	131.2532 1342.9
		80 4866.7	1143.8	445.97	352.0697 463.5055 1018.1
		20 123.1794 62.724	40.252	203.95	85.52	2854.2
	07-10-2007	50 1694.4	691.8055 77.47	2898.4	121.99	4641.8
		80 3888.4	1953.01	205.27	276.66	327.52	5561.3
		13 44.48	69.0462	70.5022	259.2161 64.2039	5362
	10-17-2007	43 1869.8	1133	564.69	478.0566 273.79	2629.4
		73 852.52	8199.8	151.83	2236.9	777.85	7400

Table 13 :

 13 ARE (10 4 ) Note that τ is the maturity, B&S is the Black & Scholes model, Hermite is the Hermite model, Edgeworth is the Edgeworth model, ML is the mixture of log-normals model, Jump is the jump-diffusion model and Heston is the Heston's stochastic volatility model.

	Dates	τ	B&S	Hermite Edgeworth ML	Jump	Heston
		20 0.1781	0.8607	5.5447	0.2046	7.12	54.61
	01-10-2007	50 0.0452	0.0356	0.1832	0.0116	0.000138 3.09
		80 0.0408	0.0406	1.9181	0.0013	0.000615 1.61
		20 0.0358	0.0194	10.67	0.0227	0.1544	44.07
	07-10-2007	50 0.0017	0.001	4.9 10 -5	0.0043	1.3 10 -4	73.1375
		80 8.6 10 -4	8.4 10 -4 2.2 10 -5	5.2 10 -4	2.5 10 -5	67.02
		13 781.1711	31.286	46.385	2416.6	3317.2	75.52
	10-17-2007	43 0.0259	0.0122	0.0219	0.0014	0.0091	162.4
		73 2.05 10 -4 0.0012	4.52 10 -5	3.37 10 -4 1.7 10 -4	62.31

A contingent asset is an asset that pays one unit of good consumed when a world state takes place and nothing else.

It should be mentioned that q(.) is the undiscounted RND whereas e -rτ q(.) represents an Arrow -Debreu state price, which is referred to as the RND.

For instance, if option prices suffer from non-synchronicity bias (that is, the underlying asset price is not observed at the same time as the option price), if the option price is fudged because of some microstructure reason (for instance, due to the bid-ask spread).

In particular , seeMerton (1992). These conditions imply that σ(K/F t,τ ) cannot be an arbitrary function but must yield an C BS (F t,τ , K, τ, r t,τ , σ(K/F t,τ )) that satisfies all the conditions of a rational option-pricing formula.

For futher details, see Aït-Sahalia and Lo (1998, 2000) and Jondeau, Ser-Huang Poon and Rockinger (2007).

see Cont (2001) for a detailed statistical analysis.

It is important to note that the interpretations of κ * and √ θ * as respectively the rate of reversion and long term volatility are still valid.

SBF-Paris Bourse provides a monthly CD-ROM including intraday values of the CAC 40 stock index and intraday transaction prices of CAC 40 options traded on the MONEP.

The constant dividend hypothesis passes over the clustering of dividends paid by most firms during specific months. However, as far as we are concerned, the bias is not relevant since we use market prices of futures to estimate an implied dividend rate. The constant dividend hypothesis could be replaced by a more realistic actual dividend assumption, leading to the following pricing formula for the futures : F t,T = S t exp(r t,T (Tt)) + n l=1 D l,t l exp(r t l ,T (Tt l )) where t l is the payment date of the l th dividend. Nevertheless, this formula requires the prediction of all dividends amounts and payment dates paid by the 40 CAC 40 index shares and the prediction of the forward interest rates.

Strike price divided by the future index level.

See Appendix B.

except on 10-17-2007 for less than one month maturity, the jump diffusion model has the smallest volatility.


C Bellégo

L Ferrara

A factor-augmented probit model for business cycle analysis

Dimension reduction of large data sets has been recently the topic of interest of many research papers dealing with macroeconomic modelling. Especially dynamic factor models have been proved to be useful for GDP nowcasting or short-term forecasting. In this paper, we put forward an innovative factor-augmented probit model in order to analyze the business cycle. Factor estimation is carried either by standard statistical methods or by allowing a richer dynamic behaviour. An application is provided on euro area data in order to point out the ability of the model to detect recessions over the period 1974-2008.

Introduction

In business cycle analysis, parametric modelling is often used by economists in central banks or governmental institutions in order to assess periodically the probability of being in a given phase of the business cycle, defined here in the NBER sense by the alternation of expansion and recession phases (see for example [START_REF] Nber | Determination of the December 2007 peak in economic activity, Business Cycle Dating Committee[END_REF]. In this respect, binary response models, such as logit and probit models, have been proved to be of great interest for practitioneers aiming at monitoring the business cycle (see [START_REF] Layton | Comparison of regime switching, probit and logit models in dating and forecasting US business cycles[END_REF], for a comparison of various parametric models). The empirical literature on this topic is huge but among the recent papers we refer for example to Estrella, Rodrigues and Schich (2003), Chauvet andPotter (2002, 2005), [START_REF] Kauppi | Predicting US recessions with dynamic binary response models[END_REF] or Nyberg (2010) for the US and to [START_REF] Moneta | Does the yield spread predict recessions in the Euro area?[END_REF], [START_REF] Duarte | Predicting real growth and the probability of recession in the Euro-area using the yield spread[END_REF] or [START_REF] Bellégo | Forecasting euro area recessions using time-varying binary response models for financial variables[END_REF] for the euro area as a whole. Binary response models establish a relationship between a low number of explanatory variables and a binary dependent variable. As an ouput, a probability of being in a given phase of the business cycle is estimated according to the values of the explanatory variables. Explanatory variables are generally chosen among a wide range of financial (interest rate term spread, asset prices ...) and macroeconomic (industrial production, employment, ...) time series. The dependent variable, denoted (r t ) t , is binary and takes the value 1 if the economy is in recession and 0 otherwise.

Recently, another strand of the literature in applied macroeconomics has focused on econometric modelling in a data-rich environment. Specifically, econometric models have been put forward to reduce the dimension of databases in order to use the most meaningful information, in a given sense. In this respect, dynamic factor models have proved their usefulness, specially in order to nowcast or forecast macroeconomic variables. Among the recent empirical literature on this topic, we refer for example to the papers of Boivin and Ng (2006), [START_REF] Schumacher | Forecasting German GDP using alternative factor models based on large datasets[END_REF][START_REF] Schumacher | Factor forecasting using international targeted predictors: The case of German GDP[END_REF], [START_REF] Giannone | Nowcasting: The real-time informational content of macroeconomic data[END_REF] or [START_REF] Barhoumi | Are disaggregate data useful for forecasting French GDP with dynamic factor models ?[END_REF]. Generally, the factor estimation step is carried out either by using standard statistical methods of dimension reduction [START_REF] Stock | Macroeconomic forecasting using diffusion indexes[END_REF] or by using estimation methods that allow a richer dynamic behaviour for the factor [START_REF] Forni | The generalized dynamic factor model: Consistency and rates[END_REF][START_REF] Moneta | Does the yield spread predict recessions in the Euro area?[END_REF][START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF]. According to the empirical literature, there is yet no consensus about the preference towards dynamic factors in modelling.

In this paper, we reconciliate both strands of the literature by putting forward a factoraugmented probit model that allows to get for each date a probability of recession conditionally to a large informational set. This conditional probability is estimated from a probit model that includes as explanatory variables the estimated dynamic or static factors. This model contributes to the recent generalisation of econometric models by incorporating synthetic factors into the standard equations, enabling thus to estimate tractable low-dimensional econometric models conditionally to a large informational set. For example, we refer to the factor-augmented VAR (FAVAR) of Bernanke, [START_REF] Bernanke | Measuring the effects of monetary policy: A factor-augmented vector autoregressive (FAVAR) approach[END_REF] or to the factor-augmented error correction model by [START_REF] Banerjee | Factor-augmented error correction model[END_REF]. As an application, we point out the ability of the factor-augmented probit model to replicate business cycles in the euro area from 1974 to 2008 by using a set of macroeconomic and financial variables that have been already proved to be reliable in-dicators of recession. Section 2 presents the model and section 3 contains the details of the application that we have carried out in order to replicate the euro area business cycle over the period January 1974 -December 2008.

Modelling

In this section we introduce the factor-augmented probit model that we put forward to handle large databases in order estimate the occurence of turning points in the business cycle through a probit model, conditionally to this large information.

Assume we observe a vector of n stationary zero mean time series x t = [x 1t , ..., x nt ] ′ , for t = 1, . . . , T . The objective of dynamic factor modelling is to decompose (x t ) t into a sum of two mutually orthogonal unobservable components: a common component of low dimension (χ t ) t , summarizing the dynamics common to all the series, and an idiosyncratic component (ξ t ) t , specific to each series. The common component (χ t ) t is supposed to linearly summarize the common behaviour of the n series. For t = 1, . . . , T , the static factor model is defined by

x t = Λf t + ξ t, (1)
where Λ is the loading matrix of dimension (n × r), the common component χ t = Λf t is driven by a small number r of factors f t common to all the variables in the model such that f t = [f 1t , ...f rt ], and ξ t = [ξ 1t , ..., ξ nt ] ′ is a vector of n idiosyncratic mutually uncorrelated components, driven by variable-specific shocks.

To take dynamics into account in modelling, an alternative specification integrates explicitly the dynamics of the factors f t . Specifically the dynamic factor representation is supposed to be given by the following equation

x t = A(L)f t + ξ t , (2)
where the common component χ t = A(L)f t integrates a linear dynamics where A(L) is a (n × r) matrix describing the autoregressive form of the r factors. If we assume that there exists a (n × q) matrix B(L) such that B(L) = A(L)N(L) with N(L) of dimension (r × q), then the dynamic factor is such that f t = N(L)U t where U t is a (q × 1) independent vector containing the dynamic shocks. It follows that the factor dynamics are described by

A(L)f t = B(L)U t , (3)
which specifies a VAR model for the factor f t .

We assume now that we observe the values of a binary variable (r t ) t that takes for value 1 when the economy is in recession at date t and 0 otherwise. Binary response models rely on the assumption that the values of the binary dependent variable, (r t ) t stem from a latent continuous variable, denoted (y t ) t , defined by the following general linear equation, for all t,:

y t = α + β ′ 0 f t + . . . + β ′ k f t-k + ε t , (4)
where for j = 0, . . . , k, f t-j = (f 1 t-j , . . . , f r t-j ) ′ is a r-vector of lagged dynamic factors used as explanatory variables, β j = (β 1 j , . . . , β r j ) ′ is a r-vector parameter and (ε t ) t is the error term supposed to be a strong white noise process with finite variance σ 2 ε .

The binary response modelling relies on the following relationship between (r t ) t and (y t ) t

r t = 1 if y t ≤ 0, 0 if y t > 0 (5)
For each date t, it can be easily proved that the conditional probability of recession is given by

P (r t = 1|I t ) = F (-α -β ′ 0 f t -. . . -β ′ k f t-k ), (6)
where I t represents the large informational set available at date t and F (.) is the cumulative density function of (ε t ) t . We assume here that the distribution of (ε t ) t is supposed to be Gaussian in order to get the probit specification of the model (the logit specification is obtained when F (.) is the logistic function).

Parameter estimation of equations ( 1) to ( 6) is carried out using a two-step procedure. First, dynamic factors ( ft ) t are estimated, then those estimated factors are put into equation (4) in order to estimate the conditional probability of being in recession given in [START_REF] Bernanke | Measuring the effects of monetary policy: A factor-augmented vector autoregressive (FAVAR) approach[END_REF].

Note that this two-step procedure implicitly incorporates the uncertainty inherent to the measure of factors.

Regarding factor estimation in the static framework of equation ( 1), we implement the procedure of Stock and Watson (2002) using static principal component analysis (PCA) to estimate the factors ft . An eigenvalue decomposition of the empirically estimated covariance matrix, T -1 T t=1 x t x ′ t , provides the (n × r) eigenvector matrix S = ( S 1 , ..., S r ) containing the eigenvectors S j corresponding to the r largest eigenvalues for j = 1, . . . , r. The factor estimates are the first r principal components of x t defined as

f sta t = S ′ x t .
In the dynamic framework of equations ( 2) and (3), factor estimation is carried out by using a 2-step Kalman filter to compute ML estimates as proposed by Doz, [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF]. This approach consists in estimating first the parameters by PCA, then, in the second step, the model is put into a state-space form and factors are estimated via Kalman smoothing. We note f dyn t the estimated factors.

The choice of the number of factors r to include in equation ( 4) is determined by using the [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] test. We implement the version of the test that determines r by minimizing the following penalized information criterion

IC(k) = ln [V (k, F )] + k × n + T nT ln nT n + T , (7)
where V (r, F ) is the sum of squared residuals such that: [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] also put forward a variety of information criteria that vary according to the penalization function.

V (k, F ) = (nT ) -1 n i=1 T t=1 (x it -Λf t ) 2 . Note that
To assess the goodness-of-fit of binary response models, Estrella (1998) proposes a P seudo-R 2 measure given by the following expression

P seudo -R 2 = 1 - L u L c -(2/n)Lc (8)
where L u is the log-likelihood of the considered model and L c is the log-likelihood of a reference nested model to which L u is compared. By construction, the nested model must have a lower log-likelihood value than the basic model. This P seudo -R 2 measure is often used in the literature for comparison between models (see for example [START_REF] Estrella | How stable is the predictive power of the yield curve ? Evidence from Germany and the United States[END_REF].

We also decide to focus on a general goodness-of-fit criterion to assess the quality of the models, often used in business cycle analysis (see, among others, [START_REF] Anas | A system for dating and detecting turning points in the euro area[END_REF], namely the quadratic probability score (QPS) given by

QP S = 1 T T t=1 (r t -P (r t = 1|I t )) 2 , (9)
where P (r t = 1|I t ) is the estimated conditional probability of being in recession, conditionally to current large information and estimated parameters.

Empirical aspects

In this section we apply the previous factor-augmented probit model to euro area data to assess the reliability of the model to replicate recessions in the zone, from January 1974 to December 2008.

Data description

We focus only on data that have been proved to be useful to detect recession and that are available since the early seventies. In this respect, we first consider industrial production and unemployment rate, for both the US and the euro area, which are generally assessed to estimate recession dates, as for example by the NBER Dating Committee. As other useful series, we include the consumer prices index. Note that, because of a lack of historical availability, we do not consider business surveys. Moreover, business surveys have not proved so far their reliability for recession detection.

Results

First we carry out a test to determine the number r of factors to include in the factoraugmented probit model given in equation ( 6). The [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] test leads us to assume that r = 2 factors, denoted (f 1t ) t and (f 2t ) t , are sufficient to provide a good description of the database. Then the two factors are estimated over the complete period using both the static model, denoted by ( f sta 1t ) t and ( f sta 2t ) t , and the dynamic model, denoted by ( f dyn 1t ) t and ( f dyn 2t ) t . Note that the two first principal factors represent around 30 percent of the total variance and that other factors possess an idiosynchratic variance lower than 10 percent.

The first estimated factor-augmented probit model is the model ( 6) without any lags (k = 0) whose conditional probability of being in recession is given by

P (r t = 1|I t ) = F (-α -β 1 0 f1t -β 2 0 f2t ), (10)
where I t is the information set available at time t. Parameter estimation results are presented in Table 2, for both static and dynamic factors. The first factor always possesses a negative coefficient, meaning thus that a decrease in this factor leads to an improvement of the recession probability, while the second one presents a positive coefficient. The dynamic version of the model presents better goodness-of-fit criteria than the static version (lower QPS and higher Pseudo-R 2 ). Estimated probabilities are presented in Figure 1. Both approaches estimate similar conditional probabilities of being in recession, although the dynamic approach provides a smoother probability. This fact is related to the more persistent behaviour of the dynamic factor that integrates autoregressive lags by definition. 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 It turns out that upward movements in the conditional probability generally correspond to benchmark recession phases. However, the double-dip movement in early eighties is not clearly reproduced by the models and in addition the conditional probabilities barely cross the critical value of 0.50 during the 1992-93 recession. This latter point reflects the fact that the 1992-93 recession belongs to the Great Moderation period characterized by a lower amplitude in macroeconomic fluctuations of industrialized countries. The maximum value of 1.0 is only reached briefly during the first oil shock and during the last 2008 recession. We also note that the year 1977 is characterized by an upward movement in the estimated probabilies that stay however largely below the threshold of 0.50.

Dynamic Factor-Augmented Probit

In a second step, we take advantage of the leading behaviour of financial variables over recession periods. Indeed it is well known that financial data possess a forward-looking property and can be thus taken as leading indicators of recession with an average lead ranging from 3 to 12 months over the business cycle (see for example [START_REF] Estrella | How stable is the predictive power of the yield curve ? Evidence from Germany and the United States[END_REF][START_REF] Moneta | Does the yield spread predict recessions in the Euro area?[END_REF][START_REF] King | Financial market perceptions of recession risk[END_REF][START_REF] Bellégo | Forecasting euro area recessions using time-varying binary response models for financial variables[END_REF]. Moreover, real data are well known by practitioners to reflect economic activity in a coincident manner. In this respect, we allow the model to reflect both coincident and leading properties in the dataset by imposing the constraint that β 1 = . . . = β k-1 = 0 in equation ( 6), for several values of k such that k ∈ {3, 6, 9, 12}. Thus, the conditional probability of being in recession at any date t is given by

P (r t = 1|I t ) = F (-α -β 1 0 f1t -β 2 0 f2t -β 1 k f1,t-k -β 2 k f2,t-k ), (11)
where k = 3, 6, 9, 12 months.

Parameter estimation results are presented in Table 2 for the various lags k. From Table 2, we note that the coincident factors always have a significant impact on the business cycle according to the t-stat statistics. For dynamic factors, the first factor lagged k months, f dyn 1,t-k , is not significant for lags k = 9 and k = 12. In opposition, the influence of the second factor tends to increase with k, for both static and dynamic approaches. This evolution may reflect the growing impact of financial variables when the horizon increases. This could mean that the second factor is related to the financial variables while the first one reflects rather the contemporaneous behaviour of hard data. It turns out that the inclusion of a lag in the factor-augmented equation enables to improve the goodness-of-fit of the models. The optimal goodness-of-fit criterion is reached for k = 12 months, which can be interpreted as the estimated lead of financial variables over recessions. Estimated probabilities of recession for k = 12 are presented in Figure 2. We observe that conditional probability values have been overall increased and reach now the maximum value of 1.0 at several consecutive dates, leading thus to a more understandable signal. The critical value of 0.50 is now largely crossed during the 1992-93 recession, but a double-dip signal now appears, reflecting perhaps the US 1990-91 recession. Again, the dynamic factor-augmented probit model provides a smoother signal, easier to interpret. For example, regarding the timing of the last recession, the dynamic version of the model with k = 12 identifies a peak in February 2008 (benchmark peak is located in March 2008), when the probability crosses durably the critical value of 0.50, while the static version is noisier and only identifies a readable peak in June 2008.

Conclusions

In this paper, we have put forward an innovative factor-augmented probit model able to replicate the business cycle. This model enables to summarize a large database into few estimated factors that are then included into a probit model. As an output we get an estimated conditional probability of being in recession each month. An application is provided on euro area data in order to point out the ability of the model to detect recession. We provide a model specification that enables to take the leading pattern of financial variables contained in the data set into account. As further research, it would be worthwile to check the robustness of the model using real-time data in order to monitor the euro area business cycle. 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004

Dynamic Factor-Augmented Probit

Figure 1 :

 1 Figure 1: Estimated probabilities with k = 0 of being in recession in the euro area stemming from the factor-augmented probit model using both dynamic (top) and a static (bottom) factors, from Feb. 1974 to Dec. 2008. Shaded areas correspond to the benchmark chronology of recessions (source Anas et al., 2007).

Figure 2 :

 2 Figure 2: Estimated probabilities with k = 12 of being in recession in the euro area stemming from the factor-augmented probit model using both dynamic (top) and a static (bottom) factors, from Feb. 1974 to Dec. 2008. Shaded areas correspond to the benchmark chronology of recessions (source Anas et al., 2007).

  Many research papers have pointed out the leading ability of financial variables (see or exampleStock and Watson, 2003, Estrella, Rodrigues and[START_REF] Estrella | How stable is the predictive power of the yield curve ? Evidence from Germany and the United States[END_REF], or King, Levin and Perli, 2007, and the references therein). Thus in this paper, we test several variables from financial markets in the euro area, such as stock maket index, oil prices, various term spread, monetary aggegates ..., with relatively long sample (we refer to[START_REF] Bellégo | Forecasting euro area recessions using time-varying binary response models for financial variables[END_REF], for a detailed description of this database).Finally, we get a medium-size database of 18 variables from January 1974 to December 2008, 5 from real economy and 13 from financial markets (see Table1for details). For each variable, we determine whether the variable is stationary or not by using standard tests. When the variable is stationary, we use it without any transformation (i.e. in level), otherwise we transformate it by taking the growth rate, or the differences, over one month.Concerning the benchmark dating chronology (r t ) t , we use the business cycle turning point chronology proposed by[START_REF] Anas | A turning point chronology for the Euro-zone classical and growth cycles[END_REF]. They provide both a quarterly chronology based on GDP and a monthly chronology based on IPI. As we need a monthly reference for our exercice, we choose the IPI chronology except for the year 2001 in which for the first time in euro area countries an industrial recession occured without a global recession. Since 1974, we retain thus four recession periods until 2006 from peak to trough: March 1974 -March 1975, January 1980 -December 1980, September 1981 -December 1982 and February 1992 -February 1993. Regarding the last recession, we assume that a peak occured in March 2008 (see also CEPR, 2009) and that the economy is still in recession in December 2008, the last point of the sample.

Table 2 :

 2 Estimated parameters (t-stat in parenthesis) and criteria for the factoraugmented probit model with 2 coincident factors f 1t and f 2t and their lagged versions for a given lag k, f 1,t-k and f 2,t-k