Oleg Kudryavtsev

INFORMATIQUE ET EN AUTOMATIQUE Efficient pricing options under regime switching

Keywords: Lévy processes, barrier options, regime switching models, Wiener-Hopf factorization, Laplace transform, numerical methods, numerical transform inversion Procesus de Lévy, options barrière, modèles de changements de régime, factorisation de Wiener-Hopf, transformée de Laplace, méthodes numériques, transformation numérique inverse

In the paper, we propose two new efficient methods for pricing barrier option in wide classes of Lévy processes with/without regime switching. Both methods are based on the numerical Laplace transform inversion formulae and the Fast Wiener-Hopf factorization method developed in Kudryavtsev and Levendorski ǐ (Finance Stoch. 13: 531-562, 2009). The first method uses the Gaver-Stehfest algorithm, the second one -the Post-Widder formula. We prove the advantage of the new methods in terms of accuracy and convergence by using Monte-Carlo simulations.

Introduction

In recent years more and more attention has been given to stochastic models of financial markets which depart from the traditional Black-Scholes model. Nowadays, a battery of models is available. One of the well-accepted models are jump diffusions or, more generally, Lévy processes. We concentrate on the one-dimensional case. For an introduction on these models applied to finance, we refer to Cont and Tankov (2004).

Regime switching Lévy models have already enjoyed much success in interpreting the behavior of a number of economic and financial series in a concise, yet parsimonious way. A Lévy process is used as the instrument that models the financial market, where the parameters of this Lévy process are allowed to depend on the state of an unobserved Markov chain that lives in continuous time. The state space may represent general financial market trends and/or other economic factors (also called "states of the world" or "regimes").

By now, there exist several large groups of relatively universal numerical methods for pricing of American and barrier options under Lévy driven financial models. Option valuation under Lèvy processes without/with regime-switching has been dealt with by a host of researchers, therefore, an exhaustive list is virtually impossible. We describe main groups of methods and several publications for each group, where the reader can find further references.

Monte Carlo methods

Monte Carlo methods perform well for pricing of barrier options in jumpdiffusion models when activity of jumps is finite because one can control the behavior of the process between the jump times, when the log-price follows a Browning bridge process (for details see [START_REF] Metwally | Using Browning bridge for fast simulation of jump-diffusion processes and barrier options[END_REF], or Cont and Tankov (2004)). In the infinite activity case, the Monte Carlo methods are much less accurate and more time consuming. Evaluation of American option prices by Monte Carlo simulation faces additional difficulties: it involves the computation of conditional expectations -see, e.g., Longstaff and Schwartz (2001). In the case of a regime switching model, one should combine the methods with Markov chain Monte Carlo simulations. An overview of Monte Carlo based methods for option pricing can be found in [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF], [START_REF] Broadie | Option pricing: valuation models and applications[END_REF], and [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF]. Generally, Monte Carlo methods consume much more time than other numerical methods.

Semi-analytical numerical methods

Methods of the second large group deal with analytical solution to the option pricing problem. The methods start with the reduction to a boundary problem for the generalization of the Black-Scholes equation (backward Kolmogorov equation); in the case of American options, a free boundary problem arises. [START_REF] Boyarchenko | Carrs randomization for finite-lived barrier options: proof of convergence[END_REF]Levendorski ǐ (1999, 2002) derived the equation for the price of a derivative security in the sense of the theory of generalized functions. Later, [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] (see also Cont and Tankov (2004)) expressed prices of European and barrier options in terms of solutions of partial integro-differential equations (PIDEs) that involve, in addition to a (possibly degenerate) second-order differential operator, a nonlocal integral term. In a regime switching framework a system of one-dimensional PIDEs has to be solved instead (see e.g. [START_REF] Chourdakis | Switching Levy Models in Continuous Time: Finite Distributions and Option Pricing[END_REF]).

Theory of pseudo-differential operators (PDO) extends the notion of a differential operator and is widely used to solve integro-differential equations. The essential idea is that a differential operator with constant coefficients can be represented as a composition of a Fourier transform, multiplication by a polynomial function, and an inverse Fourier transform. Moreover, the PDO technique based on the Fourier transfrom and the operator form of the Wiener-Hopf method is much more powerful than the technique based on the study of the kernel of the PIDE. This was the reason the theory of PDO was invented in the first placesee, e.g., [START_REF] Eskin | Boundary Problems for Elliptic Pseudo-Differential Equations[END_REF] and Hörmander (1985).

The straightforward idea to apply the PDO theory in the context of option pricing has been systematically pursued in a series of publications summarized in two monographs [START_REF] Boyarchenko | Carrs randomization for finite-lived barrier options: proof of convergence[END_REF]Levendorski ǐ (2002, 2007), and developed further in subsequent papers. In particular, Boyarchenko and Levendorski ǐ (2009) calculated the prices of American options in regime switching Lévy models. However, the general formulas for the prices involve the double Fourier inversion (and one more integration needed to calculate the factor in the Wiener-Hopf factorization formula), and hence it is difficult to implement them in practice depart from the particular cases of explicit formulas for the factors. See also Jiang and Pistorius (2008).

If the characteristic exponent of the underlying Lévy process is rational, the basic examples being the Brownian motion, Kou's model and its generalization constructed and studied in Levendorski ǐ (2004) (later, this model was used under the name Hyper-Exponential Jump-Diffusion model (HEJD)), the Wiener-Hopf factors can be derived explicitly. For a special case of diffusions with embedded exponentially distributed jumps, or more generally, for HEJD, Levendorski ǐ (2004) provides the reduction to a series of linear algebraic systems (and the solution of an equation with a monotone function on each time step), which makes the numerical procedure very fast and efficient.

In the case of processes with rational characteristic exponents, Laplace transform methods may be applied as well. First, one finds the Laplace transform of the value function of a given option with respect to the time to maturity. In [START_REF] Lipton | Assets with jumps[END_REF], [START_REF] Kou | First passage times of a jump diffusion process[END_REF], [START_REF] Sepp | Analytical Pricing of Double-Barrier Options under a Double-Exponential Jump Diffusion Process: Applications of Laplace Transform[END_REF], and in a number of other papers, e.g., [START_REF] Avram | On the valuation of constant barrier options under spectrally one-sided exponential Lévy models and Carr's approximation for American puts[END_REF] and [START_REF] Asmussen | Russian and American put options under exponential phase-type Lévy models[END_REF], the Laplace transform is derived from the distribution of the first passage time; the distribution is calculated applying the Wiener-Hopf factorization method in the form used in probability. See also [START_REF] Kyprianou | options and Canadization through fluctuation theory[END_REF], Alili and Kyprianou (2005). Once the Laplace transform is calculated, one uses a suitable numerical Laplace inversion algorithm to recover the option price. In other cases, one can approximate the initial process by Kou's model or by an HEJD, and then use the Laplace transform method (see e.g. [START_REF] Jeannin | A transform approach to compute prices and greeks of barrier options driven by a class of Lévy processes[END_REF], [START_REF] Crosby | Approximating general classes of Lévy processes by hyperexponential jump-diffusion processes[END_REF]).

However, the problem of the inversion of the Laplace transform is non-trivial from the computational point of view. There exist many different methods of numerical Laplace inversion, but some procedures, such as popular in computational finance the Gaver-Stehfest algorithm, usually require high precision.

INRIA

Notice that the latter is based on the Post-Widder inversion formula, which involves differentiation instead of integration.

We refer the reader to [START_REF] Abate | A unified framework for numerically inverting Laplace transforms[END_REF] for a description of a general framework for numerical Laplace inversion that contains the optimized version of the one-dimensional Gaver-Stehfest method. Notice that [START_REF] Sepp | Analytical Pricing of Double-Barrier Options under a Double-Exponential Jump Diffusion Process: Applications of Laplace Transform[END_REF], [START_REF] Crosby | Approximating general classes of Lévy processes by hyperexponential jump-diffusion processes[END_REF] found that the choice of 12-14 terms in the Gaver-Stehfest formula may result in satisfactory accuracy for the case of Kou and HEJD models, respectively. In this case, the standard double precision gives reasonable results.

However, often one must use high precision arithmetic, at least 50 significant digits, better, 100, according to Alex Lipton (lecture at the 2008 Bachelier Congress). The necessity of using high precision arithmetic decreases the computational speed of Laplace transform methods in option pricing considerably.

Another feature that often slows down the calculations is the fact that the values of the Laplace transform must first be found at several (at least a dozen) different points. Apart from a few cases where transform function is given by an explicit formula, the calculation of these values is time consuming.

Finally, when one uses a Laplace transform method to calculate the value function of an option, one must perform numerical Laplace inversion (as described above) separately for each initial spot price of the underlying.

Numerical methods

The next large group deals with numerical methods for the generalized Black-Scholes equation. There are four main approaches for solving PIDE: multinomial trees, finite difference schemes, Galerkin methods and numerical Wiener-Hopf factorization methods. [START_REF] Amin | Jump-diffusion option valuation in discrete time[END_REF] constructed a family of Markov chain approximations of jumpdiffusion models. Multinomial trees can be considered as special cases of explicit finite difference schemes. The main advantage of the method is simplicity of implementation; the drawbacks are inaccurate representation of the jumps and slow convergence.

Galerkin methods are based on the variational formulation of PIDE. While implementation of finite difference methods requires only a moderate programming knowledge, Galerkin methods use specialized toolboxes. Finite difference schemes use less memory than Galerkin methods, since there is no overhead for managing grids, but a refinement of the grid is more difficult. A complicated wavelet Galerkin method for pricing American options under exponential Lévy processes is constructed in [START_REF] Matache | Wavelet Galerkin Pricing of American Options on Lévy Driven Assets[END_REF]. A general drawback of variational methods is that, for processes of finite variation, the convergence can be proved in the H s -norm only, where s < 1/2; hence, the convergence in C-norm is not guaranteed.

In a finite difference scheme, derivatives are replaced by finite differences. In the presence of jumps, one needs to discretize the integral term as well. Finite difference schemes were applied to pricing barrier options in [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF], and to pricing American options in [START_REF] Carr | Why be backward?[END_REF], [START_REF] Hirsa | Pricing American options under Variance Gamma[END_REF] and Levendorski ǐ et al. (2006). Wang et al (2007) calculate prices of American options using the penalty method and a finite difference scheme.

Construction of any finite difference scheme involves discretization in space and time, truncation of large jumps and approximation of small jumps. Trun-cation of large jumps is necessary because an infinite sum cannot be calculated; approximation of small jumps is needed when Lévy measure diverges at zero. The result is a linear system that needs to be solved at each time step, starting from payoff function. In the general case, solution of the system on each time step by a linear solver requires O(m 2 ) operations (m is a number of space points), which is too time consuming. In [START_REF] Carr | Why be backward?[END_REF], [START_REF] Hirsa | Pricing American options under Variance Gamma[END_REF], and [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] the integral part is computed using the solution from the previous time step, while the differential term is treated implicitly. This leads to the explicit-implicit scheme, with tridiagonal system which can be solved in O(m ln m) operations. Levendorski ǐ et al. (2006) use the implicit scheme and the iteration method at each time step. The methods in [START_REF] Carr | Why be backward?[END_REF], [START_REF] Hirsa | Pricing American options under Variance Gamma[END_REF] and Levendorski ǐ et al. ( 2006) are applicable to processes of infinite activity and finite variation; the part of the infinitesimal generator corresponding to small jumps is approximated by a differential operator of first order (additional drift component). [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] use an approximation by a differential operator of second order (additional diffusion component). The discretization scheme for PIDE in Albanese and Kuznetsov (2003) is applicable to models for which the spectrum of the infinitesimal generator can be computed in analytically closed form.

In Kudryavtsev and Levendorski ǐ (2009) the fast and accurate numerical method for pricing barrier option under wide classes of Lévy processes was developed. The Fast Wiener-Hopf method (FWH-method) constructed in the paper is based on an efficient approximation of the Wiener-Hopf factors in the exact formula for the solution and Fast Fourier Transform algorithm. Apart from finite difference schemes where the application of the metohd entails a detailed analysis of the underlying Lévy model, the FWH-method deals with the characteristic exponent of the process.

The method starts with time discretization, which can interpreted as Carr's randomization [START_REF] Carr | Randomization and the American put[END_REF]. A sequence of stationary boundary problems for a PDO on the line results. Problems of the sequence are solved by using Wiener-Hopf approach. At the next step, the inverse of the operator that solves the boundary problem must be approximated depart from the finite difference schemes where an approximation of the infinitesimal generator is used. Generally, an approximation of the inverse can be expected to perform better. [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF] demonstrate the advantage of the FWHmethod over finite difference schemes in terms of accuracy and convergence using Monte-Carlo simulations.

Under the main property of a regime switching model, to take various economical factors into consideration in modelling, the state space of the driving Markov chain is inevitably large. As a consequence, the computational complexity involved in option valuation becomes a serious issue. Comparing the numerical methods described above, one is tempted to conclude that the FWHmethod should be preferred as rather simple, fairly fast and accurate method. Thus we have chosen to generalize the FWH-method to the case of regime switching Lévy models. Moreover, we suggest the improvement of the method which should significantly reduce the computational complexity.

INRIA

Enhanced FWH-methods

In the present paper, we introduce two enhanced FWH-methods based on the numerical Laplace transform inversion. The methods developed in the paper can be applied to pricing barrier option under wide classes of Lévy processes with/without regime switching; in the following publications, we will apply these methods to American options.

The idea behind our approach is to transform the problem to a space where the solution is relatively easy to obtain by using the Fast Wiener-Hopf factorization method. Apart from particular cases where Laplace transform is given by an explicit expression, the methods developed in the paper can be applicable for the general case. The Laplace transform maps the generalized Black-Scholes equation with the appropriate boundary conditions into the one-dimentional problem on the half-line parametrically dependent on the transform parameter.

In our first approach, we solve the problems obtained by using the FWHmethod at real positive values of the transform parameter specified by the Gaver-Stehfest algorithm. Then option prices are computed via the numerical inversion formula.

The second new approach is based on the Post-Widder formula; we find out the nth derivative of the transformed function at the certain transform parameter value by using an iterative procedure which is nothing but Carr's randomization in the FWH-method. We repeat the procedure several times for different values of n and apply the convergence acceleration algorithm of [START_REF] Abate | Numerical inversion of Laplace transforms of probability distributions[END_REF].

After straightforward modifications the both methods are applicable to the regime switching case.

The rest of the paper is organized as follows: in Section 2 we list the necessary facts of the theory of Lévy processes and regime switching models. Section 3 reviews Fast Wiener-Hopf factorization method developed in Kudryavtsev and Levendorski ǐ (2009), and introduces two enhanced FWH-methods based on the numerical Laplace transform inversion. In Section 4 we generalize the FWHmethods to the case of regime switching Lévy models for pricing barrier options. In Section 5, we produce numerical examples, and compare the results obtained by different methods; Section 6 concludes.

2 Lévy processes and the regime structure

Lévy processes: general definitions

A Lévy process is a stochastically continuous process with stationary independent increments (for general definitions, see e.g. Sato (1999)). A Lévy process may have a Gaussian component and/or pure jump component. The latter is characterized by the density of jumps, which is called the Lévy density. We denote it by F (dy). A Lévy process X t can be completely specified by its characteristic exponent, ψ, definable from the equality E[e iξX(t) ] = e -tψ(ξ) (we confine ourselves to the one-dimensional case). The characteristic exponent is given by the Lévy-Khintchine formula:

ψ(ξ) = σ 2 2 ξ 2 -iµξ + +∞ -∞ (1 -e iξy + iξy1 |y|≤1 )F (dy), (1)
where σ 2 ≥ 0 is the variance of the Gaussian component, and the Lévy measure

F (dy) satisfies R\{0} min{1, y 2 }F (dy) < +∞. (2)
Assume that under a risk-neutral measure chosen by the market, the stock has the dynamics S t = e Xt . Then we must have E[e Xt ] < +∞, and, therefore, ψ must admit the analytic continuation into a strip Im ξ ∈ (-1, 0) and continuous continuation into the closed strip Im ξ ∈ [-1, 0]. Further, if the riskless rate, r, is constant, and the stock does not pay dividends, then the discounted price process must be a martingale. Equivalently, the following condition (the EMMrequirement) must hold

r + ψ(-i) = 0, (3)
which can be used to express µ via the other parameters of the Lévy process:

µ = r - σ 2 2 + +∞ -∞ (1 -e y + y1 |y|≤1 )F (dy). (4)
The infinitesimal generator of X, denote it L, is an integro-differential operator which acts as follows:

Lu(x) = σ 2 2 u ′′ (x) + µu ′ (x) + +∞ -∞ (u(x + y) -u(x) -y1 |y|≤1 u ′ (x))F (dy). (5)
The infinitesimal generator L also can be represented as a pseudo-differential operator (PDO) with the symbol -ψ(ξ): L = -ψ(D). Recall that a PDO A = a(D) acts as follows:

Au(x) = (2π) -1 +∞ -∞ e ixξ a(ξ)û(ξ)dξ, ( 6
)
where û is the Fourier transform of a function u:

û(ξ) = +∞ -∞ e -ixξ u(x)dx.
Note that the inverse Fourier transform in ( 6) is defined in the classical sense only if the symbol a(ξ) and function û(ξ) are sufficiently nice. In general, one defines the (inverse) Fourier transform by duality.

Regular Lévy processes of exponential type

Loosely speaking, a Lévy process X is called a Regular Lévy Process of Exponential type (RLPE) if its Lévy density has a polynomial singularity at the origin and decays exponentially at the infinity (see [START_REF] Boyarchenko | Non-Gaussian Merton-Black-Scholes theory[END_REF]). An almost equivalent definition is: the characteristic exponent is analytic in a strip Im ξ ∈ (λ -, λ + ), λ -< -1 < 0 < λ + , continuous up to the boundary of the strip, and admits the representation

ψ(ξ) = -iµξ + φ(ξ), (7)
where φ(ξ) stabilizes to a positively homogeneous function at the infinity:

φ(ξ) ∼ c ± |ξ| ν , as Re ξ → ±∞, in the strip Im ξ ∈ (λ -, λ + ), (8)
INRIA where c ± > 0. "Almost" means that the majority of classes of Lévy processes used in empirical studies of financial markets satisfy conditions of both definitions. These classes are: Brownian motion, Kou's model [START_REF] Kou | A jump-diffusion model for option pricing[END_REF]), Hyperbolic processes [START_REF] Eberlein | Hyperbolic distributions in finance[END_REF], Eberlein at all (1998)), Normal Inverse Gaussian processes and their generalization (Barndorff-Nielsen (1998) and Barndorff-Nielsen and Levendorski ǐ (2001)), and extended Koponen's family. [START_REF] Koponen | Analytic approach to the problem of convergence of truncated Lévy flights towards the Gaussian stochastic process[END_REF]

ψ(ξ) = -iµξ + cΓ(-ν)[λ ν + -(λ + + iξ) ν + (-λ -) ν -(-λ --iξ) ν ], (9)
where c > 0, µ ∈ R, and λ -< -1 < 0 < λ + . Note that Boyarchenko and Levendorski ǐ (2000, 2002) consider a more general version with c ± instead of c, as well as the case ν = 1 and cases of different exponents ν ± . If ν ≥ 1 or µ = 0, then the order of the KoBoL process equals to the order of the infinitesimal generator as PDO, but if ν < 1 and µ = 0, then the order of the process is ν, and the order of the PDO -L = ψ(D) is 1. Example 2.2. If Lévy density is given by exponential functions on negative and positive axis:

F (dy) = 1 (-∞;0) (y)c + λ + e λ+y dy + 1 (0;+∞) (y)c -(-λ -)e λ-y ,
where c ± ≥ 0 and λ -< -1 < 0 < λ + , then we obtain Kou model. The characteristic exponent of the process is of the form

ψ(ξ) = σ 2 2 ξ 2 -iµξ + ic + ξ λ + + iξ + ic -ξ λ -+ iξ . ( 10
)
The version with one-sided jumps is due to Das and Foresi (1996), the two-sided version was introduced in Duffie, Pan and Singleton (2000), see also [START_REF] Kou | A jump-diffusion model for option pricing[END_REF].

The Wiener-Hopf factorization

There are several forms of the Wiener-Hopf factorization. The Wiener-Hopf factorization formula used in probability reads:

E[e iξXT ] = E[e iξ XT ]E[e iξX T ], ∀ ξ ∈ R, (11)
where T ∼ Exp q, and Xt = sup 0≤s≤t X s and X t = inf 0≤s≤t X s are the supremum and infimum processes. Introducing the notation

φ + q (ξ) = qE ∞ 0 e -qt e iξ Xt dt = E e iξ XT , (12)
φ - q (ξ) = qE ∞ 0 e -qt e iξX t dt = E e iξX T ( 13
) RR n°7184
we can write [START_REF] Boyarchenko | Generalizations of the Black-Scholes equation for truncated Lévy processes[END_REF] as

q q + ψ(ξ) = φ + q (ξ)φ - q (ξ). (14)
Equation ( 14) is a special case of the Wiener-Hopf factorization of the symbol of a PDO. In applications to Lévy processes, the symbol is q/(q + ψ(ξ)), and the PDO is E q := q/(q -L) = q(q + ψ(D)) -1 : the normalized resolvent of the process X t or, using the terminology of Boyarchenko and Levendorski ǐ (2005,2006,2007), the expected present value operator (EPV-operator) of the process X t . The name is due to the observation that, for a stream g(X t ),

E q g(x) = E +∞ 0 qe -qt g(X t )dt | X 0 = x .
Introduce the following operators:

E ± q := φ ± q (D), (15)
which also admit interpretation as the EPV-operators under supremum and infimum processes. One of the basic observations in the theory of PDO is that the product of symbols corresponds to the product of operators. In our case, it follows from ( 14) that

E q = E + q E - q = E - q E + q (16)
as operators in appropriate function spaces. For a wide class of Lévy models E and E ± admit interpretation as expectation operators:

E q g(x) = +∞ -∞ g(x + y)P q (y)dy, E ± q g(x) = +∞ -∞
g(x + y)P ± q (y)dy, where P q (y), P ± q (y) are certain probability densities with

P ± q (y) = 0, ∀ ± y < 0.
Moreover, characteristic functions of the distributions P q (y) and P ± q (y) are q(q + ψ(ξ)) -1 and φ ± q (ξ), respectively.

The regime switching Lévy process

Let I = {1, 2, . . . , d} be the space of all financial market states. Consider a continuous-time Markov chain Z t , taking values in I. Denote the generator of Z t with the transition rate matrix Λ = (λ kj ), where k, j belong to I. Notice that the off-diagonal elements of Λ must be non-negative and the diagonal elements must satisfy λ kk =j =k λ kj . Recall, given that the process Z t starts in a state k at time t 1 , it has made the transition to some other state j at time t 2 with probability given by

P (Z t2 = j|Z t1 = k) = {exp((t 2 -t 1 )Λ)} kj .
We will assume that the underlying asset price takes the form S t = S 0 e Xt , where the log-price process X t will be constructed from a collection of Lévy processes, as follows.

INRIA

Consider a collection of independent Lévy processes X k , k ∈ I. Given that Z t = k, we assume that the joint stock price process S t follows a one-dimensional exponential Lévy process with characteristic exponent ψ k . The drift terms µ k of each state are assumed prefixed by the EMM-requirement ψ k (-i) + r = 0, where r > 0 is a riskless rate. The increments of the log-price process will switch between the d Lévy processes, depending on the state Z t . Thus, this modeling assumption can be written as

dX t = dX Zt t . (17)
The infinitesimal generator of the process X t (see e.g. Chourdakis (2005)), conditional on X 0 = x and Z 0 = j is equal to

Lf (x, j) = (λ jj + L j )f (x, j) + d k =j λ jk f (x, k), (18)
where L j is the generator of the process X j t .

The system of the generalized Black-Scholes equations

The price of any derivative contract, V (t, X t ), will satisfy the Feynman-Kac formula, that is to say

(∂ t + L -r)V (t, x) = 0, (19)
where x denotes the (normalized) log-price, t denotes the time, and L is the infinitesimal generator (under risk-neutral measure).

For the sake of brevity, consider the down-and-out put option without rebate, with strike K, maturity T and barrier H < K, on a non-dividend paying stock S t . Therefore, for the one-state Lévy process X t = ln(S t /H) with the generator (5), the derivative price, V (t, X t ), will satisfy the following partial integro-differential equation (or more general pseudo-differential equation) with the appropriate initial and boundary conditions. See details in Boyarchenko and Levendorski ǐ (2002); and Cont and Tankov (2004).

(∂ t + L -r)V (t, x) = 0, t < T, x > 0, (20)
V (T, x) = (K -He x ) + , x > 0 (21) V (t, x) = 0, t ≤ T, x ≤ 0, (22)
where a + = max{a, 0}. In addition, V must be bounded.

If the characteristic exponent ψ is sufficiently regular (e.g. X t belongs to the class of RLPE), then the general technique of the theory of PDO can be applied to show that a bounded solution, which is continuous on supp V ⊂ (-∞, T ) × (0, +∞), is unique -see, e.g., [START_REF] Kudryavtsev | Pricing of first touch digitals under normal inverse Gaussian processes[END_REF].

In a regime switching setting we will have to deal with the conditional (on the regime j) option values V (t, x, j). Under the regime switching structure, a system of PIDEs will have to be solved.

(∂ t + λ jj + L j -r)V (t, x, j) + d k =j λ jk V (t, x, k) = 0, t < T, x > 0, (23)
V (T, x, j) = (K -He x ) + , x > 0, ( 24
) V (t, x, j) = 0, t ≤ T, x ≤ 0. ( 25
)
Here, L j represents the infinitesimal generator of the jth Lévy process. It is possible to apply any of the usual finite-difference schemes to this system of PIDEs to solve the problem. However, as discussed earlier, it faces difficulties due to the non-local integral terms. Instead, we develop the enhanced versions of the Fast Wiener-Hopf factorization algorithm (see Kudryavtsev and Levendorski ǐ (2009)) which are applicable to pricing barrier options under regime switching Lévy models.

3 Laplace transform in the context of the FWHmethod

Numerical Laplace transform inversion: an overview

The Laplace transform is one of the classical methods for solving partial (integro)differential equations which maps the problem to a space where the solution is relatively easy to obtain. The corresponding solution is referred to as the solution in the Laplace domain. In our case, the original function can not be retrieved analytically via computing the Bromwich's integral. Hence, the numerical inversion is needed.

Recall that popular in computational finance the Gaver-Stehfest algorithm for inverting Laplace transforms is related to the Post-Widder inversion formula. If f (τ ) is a function of a nonnegative real variable τ and the Laplace transform f (λ) = ∞ 0 e -λτ f (τ ) dτ is known, the approximate Post-Widder formula for f (τ ) can be written as

f (τ ) = lim N →∞ f N (τ ); ( 26
)
f N (τ ) := (-1) N N ! N + 1 τ N +1 f (N ) N + 1 τ , (27)
where f (N ) (λ) -N th derivative of the Laplace transform f at λ. It is well known that the convergence f N (τ ) to f (τ ) as N → ∞ is slow (of order N -1 ), so acceleration is needed. In order to enhance the accuracy, Abate and Whitt (1995) use a linear combination of the terms, i.e.,

f N,m (τ ) = m k=1 w(k, m)f N k (τ ), (28)
w(k, m) = (-1) m-k k m k!(m -k)! . ( 29
)
In this case, convergence f N,m (τ ) to f (τ ) is of order N -m . The methods of numerical Laplace inversion that fit the framework of Abate-Whitt (2006) have the following general feature: the approximate formula for f (τ ) can be written as

f (τ ) ≈ 1 τ N k=1 ω k • f α k τ , 0 < τ < ∞, ( 30
)
where N is a positive integer and α k , ω k are certain constants that are called the nodes and the weights, respectively. They depend on N , but not on f or INRIA on τ . In particular, the inversion formula of the Gaver-Stehfest method can written in the form [START_REF] Hirsa | Pricing American options under Variance Gamma[END_REF] with

N = 2n; (31) α k = k ln(2) ( 32
)
ω k := (-1) n+k ln(2) n! min{k,n} j=[(k+1)/2)] j n+1 C j n C j 2j C k-j j , (33)
where [x] is the greatest integer less than or equal to x and C K L = L! (L-K)!K! are the binomial coefficients. Because of the binomial coefficients in the weights, the Gaver-Stehfest algorithm tends to require high system precision in order to yield good accuracy in the calculations.

From [START_REF] Valko | Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion[END_REF], we conclude that the required system precision is about 2.2n, when the parameter is n. The precision requirement is driven by the coefficients ω k in [START_REF] Jiang | On perpetual American put valuation and first-passage in a regime-switching model with jumps[END_REF]. Such a high level of precision is not required for the computation of the transform f . In particular, for n = 7 standard double precision gives reasonable results. Since constants ω k do not depend on τ they can be tabulated for the values of n that are commonly used in computational finance (e.g., 6 or 7).

The Fast Wiener-Hopf factorization method

We briefly review the framework proposed by Kudryavtsev and Levendorski ǐ (2009). The main contribution of the FWH-method is an efficient numerical realization of EPV-operators E, E + and E -.

Recall that we consider the procedure for approximations of the Wiener-Hopf factors for the symbol q/(q + ψ(ξ) with ψ being characteristic exponent of RLPE of order ν ∈ (0; 2] and exponential type [λ -; λ + ]. The first ingredient is the reduction of the factorization problems to symbols of order 0, which stabilize at infinity to some constant. Introduce functions

Λ -(ξ) = λ ν+/2 + (λ + + iξ) -ν+/2 ; (34) Λ + (ξ) = (-λ -) ν-/2 (-λ --iξ) -ν-/2 ; (35) Φ(ξ) = q (q + ψ(ξ))Λ + (ξ)Λ -(ξ) -1 . (36)
Choices of ν + and ν -depend on properties of ψ, hence on order ν (see ( 7)-( 8)) and drift µ. See details in Kudryavtsev and Levendorski ǐ (2009). First, approximate Φ by a periodic function with a large period 2π/h, which is the length of the truncated region in the frequency domain, then approximate the latter by a partial sum of the Fourier series, and, finally, use the factorization of the latter instead of the exact one.

Explicit formulas for approximations of φ ± have the following form. For small positive h and large even M , set

b h k = h 2π π/h -π/h ln Φ(ξ)e -iξkh dξ, k = 0, b + h,M (ξ) = M/2 k=1 b h k (exp(iξkh) -1), b - h,M (ξ) = -1 k=-M/2+1 b h k (exp(iξkh) -1); Φ ± (ξ) ≈ exp(b ± h,M (ξ)), φ ± q (ξ) = Λ ± (ξ)Φ ± (ξ).
Approximants for EPV-operators can be efficiently computed by using the Fast Fourier Transform (FFT) for real-valued functions. Consider the algorithm of the discrete Fourier transform (DFT) defined by

G l = DF T [g](l) = M-1 k=0 g k e 2πikl/M , l = 0, ..., M -1. ( 37
)
The formula for the inverse DFT which recovers the set of g k 's exactly from G l 's is:

g k = iDF T [G](k) = 1 M M-1 l=0 G l e -2πikl/M , k = 0, ..., M -1. (38)
In our case, the data consist of a real-valued array {g k } M k=0 . The resulting transform satisfies G M-l = Ḡl . Since this complex-valued array has real values G 0 and G M/2 , and M/2 -1 other independent complex values G 1 , ..., G M/2-1 , then it has the same "degrees of freedom" as the original real data set. In this case, it is efficient to use FFT algorithm for real-valued functions (see Press, W. et al (1992) for technical details). To distinguish DFT of real functions we will use notation RDFT.

Fix the space step h > 0 and number of the space points M = 2 m . Define the partitions of the normalized log-price domain [-Mh 2 ; Mh 2 ) by points x k = -Mh 2 + kh, k = 0, ..., M -1, and the frequency domain [-π h ; π h ] by points ξ l = 2πl hM , l = -M/2, ..., M/2. Then the Fourier transform of a function g on the real line can be approximated as follows:

ĝ(ξ l ) ≈ he iπl RDF T [g](l), l = 0, ..., M/2.
Here and below, z denotes the complex conjugate of z. Using the notation p(ξ) = q(q + ψ(ξ)) -1 , we can approximate E q :

(E q g)(x k ) ≈ iRDF T [p. * RDF T [g]](k), k = 0, ..., M -1. ( 39
)
Here and below, . * is the element-wise multiplication of arrays that represent the functions. Further, we define

b h k ≈ iRDF T [ln Φ](k); p ± (ξ l ) = Λ ± (ξ l ) exp(b ± h,M (ξ l )), l = -M/2, ..., 0. ( 40
)
The action of the EPV-operator E ± q is approximated as follows:

(E ± q g)(x k ) = iRDF T [ p ± . * RDF T [g]](k), k = 0, ..., M -1. ( 41
)

INRIA

The Gaver-Stehfest algoritm and the FWH-method

In our study we will apply the Laplace transform to solve the problems for PIDE [START_REF] Carr | The fine structure of asset returns: an empirical investigation[END_REF] and the system of PIDEs [START_REF] Cont | Financial modelling with jump processes[END_REF]. We start with a Lévy model without regime switching and solve the corresponding problem for pricing barrier options in the Laplace domain at real positive values of the transform parameter specified by the Gaver-Stehfest algorithm.

We introduce a new variable

τ = T -t. With a new function v(τ, x) = V (T -τ, x) the problem (20)-(22) turns into (∂ τ + r -L)v(τ, x) = 0, τ > 0, x > 0, ( 42
) v(0, x) = (K -He x ) + , x > 0 (43) v(τ, x) = 0, τ ≥ 0, x ≤ 0. ( 44
)
The Laplace transform of v(τ, x) with respect to the time variable is defined by ṽ

(λ, x) := ∞ 0 e -λτ v(τ, x) dτ,
where λ is a transform variable with positive real part, Re λ > 0. To be specific, in subsequent study we assume that λ ∈ R + . The standard rules yield

∂ τ v(τ, x) → λṽ(λ, x) -v(0, x), Lv(τ, x) → Lṽ(λ, x).
Applying Laplace transform to [START_REF] Levendorski Ǐ | Early exercise boundary and option pricing in Lévy driven models[END_REF], we obtain that ṽ(λ, x) satisfies the following equation:

(λ + r -L)ṽ(λ, x) = (K -He x ) + , x > 0, [START_REF] Lipton | Assets with jumps[END_REF] subject to the corresponding transformed boundary condition

ṽ(λ, x) = 0, x ≤ 0. ( 46
)
Given n, we can use the Gaver-Stehfest inversion formula for ṽ(λ, x) provided that the solutions to the problem ( 45),( 46) are found at λ = k ln 2/τ , k = 1, . . . , N (see ( 30)-( 33)). Set q = λ + r and denote by 1 [0,+∞) (x) the indicator function of [0, +∞). A general class of boundary problems that contains the problem ( 45)-( 46) was studied in [START_REF] Boyarchenko | Non-Gaussian Merton-Black-Scholes theory[END_REF] and Levendorski ǐ (2004b). It was shown that the unique bounded solution is given by ṽ

(λ, x) = 1 q E - q 1 [0,+∞) (x)E + q (K -He x ) + . (47)
Now, the Fast Wiener-Hopf factorization method [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF] can be applied. Since the approximate expressions for the Wiener-Hopf factors φ ± q (ξ) are available (see 40), one can calculate ṽ(λ, x) quite easily using formulas [START_REF] Levendorski Ǐ | Pricing of the American put under Lévy processes[END_REF].

It follows, that the computational complexity of the developed algorithm (as well as the FWH-method) is O(N M ln M ), where M is a number of points in the log-price space; in the case of the FWH-method, M denotes the number of time steps. The Gaver-Stehfest algorithm produces rapid convergence results already using N = 10 -14 depart from the FWH-method with N being of order 400 -800. Hence, the new method is computationally much faster (often, dozen of times faster) than original FWH-method constructed in [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF].

Our new method enjoys an additional appealing feature: it produces a set of option prices at different spot levels. Notice that in the case of the known Laplace transform methods, one must perform numerical Laplace inversion separately for each initial spot price of the underlying.

Our new algorithm provides increasing accuracy as n in the Gaver-Stehfest inversion formula increases. However, if n > 7 good accuracy results can be achieved only using a multi-precision computational environment.

The method based on the Post-Widder formula (see the next subsection) achieves similar performance to the method proposed here; however, the former method does not require high precision.

The Post-Widder formula or randomization

In this subsection, we propose the second new approach to pricing barrier options which involves the numerical Laplace transform inversion formulas ( 28), [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]. Recall that we are looking for the solution v(τ, x) to the problem ( 42)-( 44) at τ = T .

Applying the Laplace transform to the corresponding PIDE, we consider the problem ( 45), [START_REF] Madan | The variance Gamma process and option pricing[END_REF] in the Laplace domain, once again. As a basis for the Gaver-Stehfest algorithm, it was established a discrete analog of the Post-Widder formula [START_REF] Eberlein | Hyperbolic distributions in finance[END_REF] involving finite differences to approximate N th derivative of the transformed function. In fact, for performing numerical inversion we need to find ∂ N λ ṽ(λ, x). We have, on differentiating both sides of the equations ( 45), [START_REF] Madan | The variance Gamma process and option pricing[END_REF] with respect to λ:

(λ + r -L)∂ λ ṽ(λ, x) = -ṽ(λ, x), x > 0, ( 48
)
∂ λ ṽ(λ, x) = 0, x ≤ 0. ( 49
)
Repeating this procedure, for all k = 1, 2, ..., N , we obtain a sequence of the following problems

(λ + r -L)∂ k λ ṽ(λ, x) = -k∂ k-1 λ ṽ(λ, x), x > 0, ( 50
) ∂ λ ṽ(λ, x) = 0, x ≤ h. ( 51
)
Fix an integer N > 1, and set ∆τ = T /(N + 1), λ = 1/∆τ . Then we introduce the following functions:

v 0 (x) = (K -He x ) + ; ( 52
)
v k+1 (x) = (-1) k k! 1 ∆t k+1 ∂ k λ ṽ 1 ∆τ , x , k = 0, ..., N. (53)
It follows that

∂ k λ ṽ 1 ∆τ , x = (-1) k k!(∆τ ) k+1 v k+1 (x), k = 0, ..., N. (54)
Substituting expressions 1/∆τ for λ and (54) for ∂ k λ ṽ 1 ∆τ , x into ( 50)-( 51), simplifying and eliminating the multipliers from the final set of equations, one

INRIA for k = 1, ..., N + 1: (q -L)v k (x) = 1 ∆τ v k-1 (x), x > 0, ( 55
) v k (x) = 0, x ≤ 0, ( 56
)
where q = r + 1/∆τ . The sequence v k (x), k = 1, 2, ..., N + 1, is determined recurrently by means of the problem (55), (56) at each step k. It follows from [START_REF] Boyarchenko | Non-Gaussian Merton-Black-Scholes theory[END_REF] that the unique bounded solution to the problem (55), ( 56) is given by

v k = 1 q∆τ E - q 1 [0,+∞) E + q v k-1 . (57)
Once again, the Fast Wiener-Hopf factorization method [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF] can be applied. Moreover, approximate formulas for E + q , E - q (41) are at the first and last steps only. At all intermediate steps, the exact analytic expression q/(q + ψ(ξ) is used (see [START_REF] Kyprianou | options and Canadization through fluctuation theory[END_REF]). Indeed, for k = 1, 2, . . . , N + 1, define

w k = 1 [0;+∞) E + q v k-1 . (58)
Then

v k = (q∆τ ) -1 E - q w k (x). (59)
Using the Wiener-Hopf factorization formula ( 16), we obtain that

w k = (q∆τ ) -1 1 [0;+∞) E q w k-1 . (60)
Finally, we take into account the Post-Widder formula ( 26)- [START_REF] Eberlein | New insights into smile, mispricing and value at risk: The hyperbolic model[END_REF]. As a result, we conjecture that the solution v N +1 (x) to our problem converges to the unknown solution v(T, x) of the problem ( 42)-( 44), as N gets arbitrarily large with T held fixed.

Unfortunately, the Post-Widder formula provide a very poor approximation (of order N -1 ). See details in Subsection 3.1. For example, v 1000 (x) may yield an estimate to v(T, x) with only two or three digits of accuracy. To achieve a good approximation, a convergence acceleration algorithm is required for the sequence v N (x). A good candidate is the summation formula ( 28)-( 29) (see Abate and Witt (1995)). We start with the choice N = 10 and m = 3, and increase them if necessary.

Given parameters N and m in ( 28)-( 29), the computational complexity of the developed algorithm is O(N 0 M ln M ), where M is a number of points in the log-price space, and

N 0 = (N +1)(m+1)m 2 .
The new enhanced FWH-method based on the Post-Widder formula produces rapid convergence results already using N = 10 and m = 3. Hence, the new method is computationally much faster than the original FWH-method developed in [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF].

The second new method achieves similar performance to the first one constructed in the previous Subsection. Our new algorithm provides increasing accuracy as N and m in the formula (28) increase. At the same time, the method does not require a multi-precision arithmetic. Remark 3.1. Notice that our value v k (x) is also the approximation for the solution v(k∆τ, x) to the problem (42)-( 44) which arises when time is discretized and the derivative ∂ v(k∆τ, x) in ( 42) is replaced with the finite difference

(1/∆τ )(v(k∆τ, x) -v((k -1)∆τ, x)).
The notion of discretizing time while leaving space continuous is known in the numerical methods literature as the method of horizontal lines or Rothe's method (see [START_REF] Rothe | Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindimensinaler Randwerttaufgaben[END_REF]). Carr's randomization procedure [START_REF] Carr | Randomization and the American put[END_REF] indicates an alternative interpretation of approximation induced by our procedure. Notice that Carr's randomization was successfully applied to the valuation of (single and double) barrier options in a number of works. Remark 3.2. As a result, we conjecture that Carr's approximation to the value of a finite-lived barrier option with bounded continuous terminal payoff function always converges to the actual value for a wide class of Lévy processes. Moreover, we provide the order of the convergence.

Notice that the detailed probability-theoretical proof of convergence of the Carr's randomization procedure for barrier options in Lévy-driven models is given by M. [START_REF] Boyarchenko | Carrs randomization for finite-lived barrier options: proof of convergence[END_REF].

Pricing barrier options under regime switching Lévy models

In the present secton, we generalize the framework proposed by Kudryavtsev and Levendorski ǐ ( 2009), and extend it to regime switching Lévy models. Recall that we consider the down-and-out put option without rebate, with strike K, maturity T and barrier H < K, on a non-dividend paying stock e Xt , where X t is defined by [START_REF] Broadie | Option pricing: valuation models and applications[END_REF].

Regime states can either be visible or hidden from market participants. We assume that the states are visible, and the initial state is given a priori.

In a regime switching setting we will have to keep track of the conditional (on the regime j) option values V (t, x, j). It follows from Subsection 2.5 that under the regime switching structure, a problem of pricing down-and-out barrier options can be reduced to the problem ( 23)- [START_REF] Crosby | Approximating general classes of Lévy processes by hyperexponential jump-diffusion processes[END_REF].

Once again, we introduce a new variable τ = T -t. With new functions v(τ, x, j) = V (T -τ, x, j) the problem ( 23)- [START_REF] Crosby | Approximating general classes of Lévy processes by hyperexponential jump-diffusion processes[END_REF] turns into

(∂ τ -λ jj + r -L j )v(τ, x, j) - d k =j λ jk v(τ, x, k) = 0, τ > 0, x > 0, (61)
v(0, x, j) = (K -He x ) + , x > 0, (62)
v(τ, x, j) = 0, τ ≥ 0, x ≤ 0. ( 63
)
Set v 0 (x, j) = 1 [0;+∞) (x)(K -He x ) + , q j = r + 1/∆τ -λ jj , for j = 1, ..., d. We divide [0, T ] into N time steps of length ∆τ = T /N , and we introduce a vector-function

V s (x) =     v s (x, 1) v s (x, 2) ... v s (x, d)     , INRIA v s (x, j
) is an approximation to the price of the barrier option at state j and time τ s = s∆τ , s = 1, 2, . . ..

We discretize the time derivative, and find the vector-function V s (x), s = 1, 2, . . ., as a unique bounded solution to the boundary problem for the following system PIDEs (q j -L j )v(τ, x, j) -

d k =j λ jk v ( x, k) = 1 ∆τ v s-1 (x, j), x > 0, (64)
v s (x, j) = 0, x ≤ 0, (65)
where j ranges over {1, 2, ..., d}. One can rewrite (64)-( 65) as

(Q -1 (Q -L) -Q -1 Λ 0 )V s (x) = 1 ∆τ Q -1 V s-1 (x), x > h, (66)
V k (x) = 0, x ≤ 0, (67)
where, in matrix notation

L =     L 1 0 0 ... 0 0 L 2 0 ... 0 ... 0 0 0 ... L d     ; Q =     q 1 0 0 ... 0 0 q 2 0 ... 0 ... 0 0 0 ... q d     ; Λ 0 =     0 λ 12 λ 13 ... λ 1d λ 21 0 λ 23 ... λ 2d ... λ d1 λ d2 λ d3 ... 0     .
First, we factorize the operators E j = q j /(q j -L j ) (see Subsection 2.3):

E j = E - j E + j , j = 1, ..., d.
Then, we introduce the following operator, in matrix notation

Ẽ =     E - 1 1 [0;+∞) E + 1 0 0 ... 0 0 E - 2 1 [0;+∞) E + 2 0 ... 0 ... 0 0 0 ... E - N 1 [0;+∞) E + N     (68)
Applying the operator Ẽ to both sides of (66), and taking into account the fact that (57) solves the problem (55),(56), we have

(I -ẼQ -1 Λ 0 )V s (x) = 1 ∆τ Q -1 ẼV s-1 (x), x > 0, ( 69
) V s (x) = 0, x ≤ 0. ( 70
)
The system (69),(70) is easily solved by using the iteration method, with any accuracy. Before stating our iterative procedure, we introduce some notations. By C 0 (R + ; R d ) we denote the space of bounded, continuous functions from R + to R d that vanish at +∞. Here R + is the set of positive real numbers. The topology on C 0 (R + ; R d ) is defined by the norm

||V || 0 = sup x∈R + ||V (x)||, (71)
where || • || the Euclidean norm. Note that (C 0 (R + ; R d ), || • ||) turns out to be a Banach space.

Recall that the EPV-operators E qj and E ± qj admit interpretation as expectation operators (see details in Subsection 2.3). Therefore the mapping Ẽ defined by ( 68) is a continuous operator on the space C 0 (R + ; R d ) with the norm || Ẽ|| = 1.

Then we can rewrite (69) as V s = Φ(V s ), where

Φ(V ) = ẼQ -1 Λ 0 V + 1 ∆τ Q -1 ẼV s-1 (x).
Clearly, Φ maps

C 0 (R + ; R d ) into C 0 (R + ; R d ), and
Φ(V ) -Φ(U ) = ẼQ -1 Λ 0 (V -U ), V, U ∈ C 0 (R + ; R d ).
With the right choice of ∆τ , || ẼQ -1 Λ 0 || < 1, hence Φ is a contraction map, and we can calculate the solution to (69),(70) using the iteration procedure, with any accuracy. If ||Q -1 Λ 0 || is sufficiently small, then several iterations are enough to achieve a good accuracy.

If a model with d regimes is used, we can write the recursive relationship as

V l s (x) = 1 ∆τ Q -1 ẼV s-1 (x) + ẼQ -1 Λ 0 V l-1 s (x), x > 0,
where l is the iteration number, and V 0 s (x) = V s-1 (x). After straightforward modifications, the methods developed in Subsection 3.3 and Subsection 3.4 are applicable to barrier options under regime switching as well.

Numerical examples

In this section, we compare the performance of the new two methods and the original FWH-method. In numerical examples, we implement the algorithms of the enhanced FWH-methods described in Subsection 3.3 and in Subsection 3.4. We will refer to these algorithms as the FWH&GS-method and FWH&PWmethod, repectively.

Pricing options without regime switching

In Subsection 5.1, we will compare the prices from the FWH&GS-method and the FWH&PW-method against prices obtained by differnent numerical methods and reported in [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF].

We use Monte Carlo method (MC-method) and the accurate finite-difference scheme of Kudryavtsev and Levendorski ǐ (2006) (FDS-method) as the benchmarks. We also compare the performance of the FWH&GS, FWH&PW-methods and the finite difference scheme constructed in Cont and Voltchkova (2005) (CVmethod). We will show the advantage of the new methods in terms of speed over the original FWH-method.

We consider the down-and-out put option with strike K, barrier H and time to expiry T . The option prices were calculated on a PC with characteristics Intel Core(TM)2 Due CPU, 1.8GHz, RAM 1024Mb, under Windows Vista.

INRIA

and Levendorski ǐ (2009) used a PC with the same characteristics. The prices calculated by MC, FDS, CV-methods we will use are the same as in Table 1 of [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF]. For the Monte Carlo calculations Kudryavtsev and Levendorski ǐ (2009) used 500, 000 paths with time step = 0.00005. The examples, which we analyze in detail below, are fairly representative.

We consider KoBoL (CGMY) model of order ν ∈ (0, 1), with parameters σ = 0, ν = 0.5, λ + = 9, λ -= -8, c = 1. We choose instantaneous interest rate r = 0.072310, time to expiry T = 0.5 year, strike price K = 100 and the barrier H = 90. In this case, the drift parameter µ is approximately zero. The localization domain is (x min ; x max ) with x min = -ln 2 and x max = ln 2; we check separately that if we increase the domain two-fold, and the number of points 4-fold, the prices change by less than 0.0001.

Table 1, reports prices for down-and-out put options calculated by using Monte-Carlo simulation, FDS, FWH, CV-methods, with very fine grids, and FWH&PW and FWH&GS methods. The options are priced at five spot levels. ExtCV labels option prices obtained by linear extrapolation of prices V h,N with h = 0.000005 and h = 0.000002.

In Table 2, Panel A and Panel B, the sample mean values are compared with the prices computed by FDS, FWH, CV-methods, FWH&PW and FWH&GS methods. The results show a general agreement between the Monte Carlo simulation results and those computed by FDS-method and the different versions of the FWH methods. The prices from the FWH&PW and FWH&GS methods converge very fast and agree with MC-prices and FDS-prices very well (see also Table 2, Panel D). We see that FWH&PW and FWH&GS methods produce sufficiently good results in just several dozens of milliseconds.

Table 3, demonstrates the advantage of the new methods in terms of accuracy and convergence in time over the original FWH-method. In 5.2, we will compare the prices from the FWH&PW-method and the original FWH-method against prices obtained by MC-method, under the regime switching Lévy model assumption.

It is well known that the convergence of Monte Carlo estimators of quantities involving first passage is very slow. Hence, a large number of paths was needed to obtain a convergence. For the Monte Carlo calculations we used 500, 000 paths with time step = 0.00001. For simulating trajectories of the tempered stable (KoBoL) process we implemented the code of J. Poirot and P. Tankov (www.math.jussieu.fr/ ∼ tankov/). The program uses the algorithm in Madan and Yor (2005), see also Poirot and Tankov (2006). We combine this method with Markov chain Monte Carlo simulations.

We consider again the down-and-out put option with strike K, barrier H and time to expiry T . We choose instantaneous interest rate r = 0.04879, time to expiry T = 0.1 year, strike price K = 100 and the barrier H = 90.

Our example uses a three-state KoBoL (CGMY) model, with parameters presented in Table 4. The rate matrix of the underlying Markov chain is given by

Λ =   -0.8 0.5 0.3 0.2 -0.7 0.5 0.2 0.4 -0.6   .
Table 5, reports prices for down-and-out put options calculated by using Monte-Carlo simulation, FWH and FWH&PW methods. The options are priced for each intial state at four spot levels.

In Table 6, the sample mean values are compared with the prices FWH and FWH&PW methods. The prices from the FWH&PW and FWH methods converge very fast and agree with each other and MC-prices very well. We see that FWH&PW method is several times faster than the generalized version of the original FWH method.

Conclusion

In the paper, propose two new fast and accurate methods for pricing barrier options in wide classes of Lévy processes with/without regime switching. Both methods use the numerical Laplace transform inversion formulae and the Fast Wiener-Hopf factorization method developed in Kudryavtsev and Levendorski ǐ (2009). The first method uses the Gaver-Stehfest algorithm, the second onethe Post-Widder formula.

Using an accurate albeit relatively slow finite-difference algorithm developed in Levendorski ǐ et al (2006) and Monte Carlo simulations, we demonstrate the accuracy and fast convergence of the two new methods. Numerical examples show that the new methods are computationally much faster (often, dozen of times faster) than the original FWH-method constructed in [START_REF] Kudryavtsev | Fast and accurate pricing of barrier options under Levy processes[END_REF]. Our new methods enjoy an additional appealing feature: they produce a set of option prices at different spot levels, simultaneously.

The method based on the Post-Widder formula achieves similar performance to the method which uses the Gaver-Stehfest algorithm; however, the former method does not require high precision.

We notice that Carr's randomization procedure [START_REF] Carr | Randomization and the American put[END_REF] indicates an alternative interpretation of the approximation induced by our second method. As a result, we conjecture that Carr's approximation to the value of a finitelived barrier option with bounded continuous terminal payoff function always converges to the actual value for a wide class of Lévy processes. Moreover, we provide the order of the convergence.

We generalize the framework proposed by Kudryavtsev and Levendorski ǐ (2009), and extend it to regime switching Lévy models. We prove the advantage of the new methods in terms of accuracy and convergence by using Monte-Carlo simulations.

INRIA 5 . 2

 52 KoBoL parameters: ν = 0.5, λ + = 9, λ -= -8, c = 1, µ ≈ 0. Option parameters: K = 100, H = 90, r = 0.072310, T = 0.5. Algorithm parameters: h -space step, N -number of time steps (or the parameter of the FWH&PW and FWH&GS methods), S -spot price. Panel A: Down-and-out put prices calculated by using MC, FDS, FWH and CV methods. Panel B: Down-and-out put prices calculated by using MC, FDS, FWH&PW and FWH&GS methods. Pricing options with regime switching

  KoBoL parameters: ν = 0.5, λ + = 9, λ -= -8, c = 1, µ ≈ 0. Option parameters: K = 100, H = 90, r = 0.072310, T = 0.5. Algorithm parameters: h -space step, N -number of time steps (or the parameter of the FWH&PW and FWH&GS methods), S -spot price. Panel A: Down-and-out put prices calculated by using FWH and FWH&PW and FWH&GS methods; h = 0.001. Panel B: Relative errors for FWH and FWH&PW and FWH&GS methods with h = 0.001; the benchmark -down-and-out put prices calculated by using FWH with h = 0.001, N = 6400. Panel C: Down-and-out put prices calculated by using FWH and FWH&PW and FWH&GS methods; h = 0.0005. Panel D: Relative errors for FWH and FWH&PW and FWH&GS methods with h = 0.0005; the benchmark -down-and-out put prices calculated by using FWH with h = 0.0005, N = 6400.

  introduced a symmetric version; Boyarchenko and Levendorski ǐ (1999, 2000), gave a non-symmetric generalization; later a subclass of this model appeared under the name CGMY-model in Carr et al. (2002), and Boyarchenko and Levendorski ǐ (2002) used the name KoBoL family. The important exception is Variance Gamma Processes (VGP; see, e.g., Madan et al. (1998)). VGP satisfy the conditions of the first definition but not the second one, since the characteristic exponent behaves like const • ln |ξ|, as ξ → ∞.

Example 2.1. The characteristic exponent of a pure jump KoBoL process (a.k.a. CGMY model) of order ν ∈ (0, 2), ν = 1 is given by

Table 1 :

 1 KoBoL(CGMY) model regime switching: option prices

	A

Table 3 :

 3 KoBoL(CGMY) model without regime switching: FWH&PW and FWH&GS methods

									A
		FWH&PW		FWH&GS			FWH
	Spot price	N = 5		N = 10	N = 14	N = 400	N = 800	N = 1600	N = 3200	N = 6400
	S = 91	0.235146	0.235774	0.235850	0.237139	0.236491	0.236168	0.236006	0.235926
	S = 101	0.568914	0.567482	0.567368	0.568303	0.567766	0.567496	0.567361	0.567294
	S = 111	0.385695	0.385778	0.385705	0.385349	0.385557	0.385661	0.385713	0.385739
	S = 121	0.208113	0.208503	0.208576	0.208223	0.208406	0.208497	0.208543	0.208566
	S = 131	0.107411	0.107560	0.107642	0.107439	0.107516	0.107554	0.107573	0.107583
	CPU-time
	(sec)	0.016		0.031		0.047		0.21	0.421	0.842	1.669	3.386
									B
				FWH&PW		FWH&GS		FWH
	Spot price	N = 5	N = 10	N = 14	N = 400	N = 800	N = 1600	N = 3200
	S = 91	-0.33%	-0.06%	-0.03%		0.51%	0.24%	0.10%	0.03%
	S = 101	0.29%	0.03%	0.01%		0.18%	0.08%	0.04%	0.01%
	S = 111	-0.01%	0.01%	-0.01%		-0.10%	-0.05%	-0.02%	-0.01%
	S = 121	-0.22%	-0.03%	0.00%		-0.16%	-0.08%	-0.03%	-0.01%
	S = 131	-0.16%	-0.02%	0.05%		-0.13%	-0.06%	-0.03%	-0.01%
									C
		FWH&PW		FWH&GS			FWH
	Spot price	N = 5		N = 10	N = 14	N = 400	N = 800	N = 1600	N = 3200	N = 6400
	S = 91	0.234730	0.235328	0.235433	0.236720	0.236073	0.235750	0.235589	0.235508
	S = 101	0.568745	0.567312	0.567199	0.568134	0.567596	0.567327	0.567192	0.567125
	S = 111	0.385633	0.385715	0.385642	0.385287	0.385495	0.385599	0.385651	0.385677
	S = 121	0.208086	0.208476	0.208549	0.208196	0.208379	0.208470	0.208516	0.208539
	S = 131	0.107376	0.107525	0.107607	0.107403	0.107480	0.107519	0.107538	0.107548
	CPU-time
	(sec)	0.062		0.078		0.094		0.468	0.904	1.84	3.62	7.27
									D
				FWH&PW		FWH&GS		FWH
	Spot price	N = 5	N = 10	N = 14	N = 400	N = 800	N = 1600	N = 3200
	S = 91	-0.33%	-0.08%	-0.03%		0.51%	0.24%	0.10%	0.03%
	S = 101	0.29%	0.03%	0.01%		0.18%	0.08%	0.04%	0.01%
	S = 111	-0.01%	0.01%	-0.01%		-0.10%	-0.05%	-0.02%	-0.01%
	S = 121	-0.22%	-0.03%	0.00%		-0.16%	-0.08%	-0.03%	-0.01%
	S = 131	-0.16%	-0.02%	0.05%		-0.13%	-0.06%	-0.03%	-0.01%

RR n°7184

The drift terms µ j are prefixed by the EMM-requirement Table 5: KoBoL(CGMY) model with regime switching: option prices


Eric Colinet
email: eric.colinet@cea.fr

Member, IEEE Jérôme Juillard
email: jerome.juillard@supelec.fr

A Weighted Least-Squares Approach to Parameter Estimation Problems Based on Binary Measurements

Keywords: parameter estimation, quantized observations, binary sensors, FIR digital filters

We present a new approach to parameter estimation problems based on binary measurements, motivated by the need to add integrated low-cost self-test features to microfabricated devices. This approach is based on the use of original weighted least-squares criteria: as opposed to other existing methods, it requires no dithering signal and it does not rely on an approximation of the quantizer. In this paper, we focus on a simple choice for the weights and establish some asymptotical properties of the corresponding criterion. To achieve this, the assumption that the quantizer's input is Gaussian and centered is made. In this context, we prove that the proposed criterion is locally convex and that it is possible to use a simple gradient descent to find a consistent estimate of the unknown system parameters, regardless of the presence of measurement noise at the quantizer's input.

I. INTRODUCTION

N this paper, we present a new parameter estimation method based on binary measurements. This work was originally motivated by the need to add integrated low-cost self-test features to microfabricated devices, such as MEMS and NEMS. Even though there exists a wide range of applications where identification methods based on binary observations are necessary or desirable [START_REF] Wang | System identification using binary sensors[END_REF], the focus is brought here on the test of microelectronic devices. It is wellknown that, as characteristic dimensions become smaller, the dispersions afflicting electronic devices tend to become larger. Typical sources of uncertainty and dispersion are variations in the fabrication process, changes in the operating conditions or imperfect knowledge of physics. As a consequence, it is usually impossible to guarantee a priori that a given device will function properly. Expensive tests must then be run after fabrication to ensure that only suitable devices are commercialized. Furthermore, self-test (and self-tuning) features, such as parameter estimation routines, must often be implemented, so that devices can adapt to changing conditions. However, most parameter estimation methods [START_REF] Ljung | System identification -theory for the user[END_REF][START_REF] Walter | Identification of parametric models from experimental data[END_REF] rely on high-resolution digital measurements. Their integration requires the implementation of high-resolution analog-todigital converters (ADCs) and, thus, results in longer design times, larger silicon areas and increased costs. Our objective is then to develop a parameter estimation method that relies on very low-resolution (ideally binary) measurements, in order to keep the added cost of testing as small as possible.

In the field of micro-electronics, the issue of parameter estimation of linear systems from binary data has partially been addressed by Negreiros [START_REF] Negreiros | Ultimate low cost analog BIST[END_REF] and Juillard and Colinet [START_REF] Juillard | Initialization of the BIMBO self-test method using binary inputs and outputs[END_REF] 1 . In [START_REF] Juillard | Initialization of the BIMBO self-test method using binary inputs and outputs[END_REF], a white Bernoulli input is used to excite the unknown system. Provided the impulse response does not vanish too quickly, it is possible to establish an analytical relationship between the cross-covariance function of the binary outputs and inputs and the impulse response of the unknown system and, thus, to identify the system.

Outside the context of micro-electronics, the most significant contributions come from Wigren [START_REF] Wigren | Adaptive filtering using quantized output measurements[END_REF], Wang [START_REF] Wang | System identification using binary sensors[END_REF][START_REF] Wang | Joint identification of plant rational models and noise distribution functions using binary-valued observations[END_REF][START_REF] Zhao | Identification of Wiener systems with binary-valued output observations[END_REF], Rafajlowicz [START_REF] Rafaljowicz | Linear systems identification from random threshold binary data[END_REF][START_REF] Rafaljowicz | System identification from cheap, qualitative output observations[END_REF] and their co-workers. In [START_REF] Rafaljowicz | Linear systems identification from random threshold binary data[END_REF][START_REF] Rafaljowicz | System identification from cheap, qualitative output observations[END_REF], the author considers a linear system excited by a deterministic signal. The output of the system is measured via a quantizer whose threshold is "randomly specified"2 . The author then uses an elegant argument based on Von Neumann's theorem to estimate the Fourier transform of the quantizer's input from the quantized data and, knowing the inputs, it becomes simple to estimate the transfer function of the system. In [START_REF] Wang | System identification using binary sensors[END_REF], Wang et al. introduced a new method for estimating parameters from binary (or quantized) data. The unknown system is excited by a periodic signal and, as in Rafajlowicz's work, the threshold of the quantizer is randomly specified. However, the hypotheses on the quantizer's input and on the random threshold are much less strict: the major constraint is that the cumulative distribution function (cdf) of the threshold must be invertible. The parameter estimation problem then reduces to solving a linear system. This approach is generalized in [START_REF] Wang | Joint identification of plant rational models and noise distribution functions using binary-valued observations[END_REF], where it is shown that the cdf of the threshold can be estimated along with the parameters of the system, and in [START_REF] Zhao | Identification of Wiener systems with binary-valued output observations[END_REF] where the identification of infinite impulse response (IIR) systems and nonlinear Wiener systems is addressed. Finally, Wigren has developed a least mean squares (LMS) approach to the problem of online parameter estimation from binary observations [START_REF] Wigren | Adaptive filtering using quantized output measurements[END_REF]. This method is based on the definition of an approximate gradient of the (quantized) least-squares criterion. Under some stationarity hypotheses and in the absence of measurement noise at the quantizer's input, it is possible to guarantee the asymptotical convergence of this method to the nominal parameters.

Each of these methods has its merits and weaknesses. It is our opinion that Wigren's approach is the one that fits better in the context of test of electronics system: as opposed to Wang's or Rafajlowicz's methods, it does not require a varying threshold (dithering signal), which can only be implemented with a high-resolution DAC or a surface-consuming analog component. Furthermore, it can be used online as well as offline. However, it does rely on an approximation of the quantizer (and, thus, it is not an exact method, stricto sensu, as opposed to [START_REF] Wang | System identification using binary sensors[END_REF]). Moreover, convergence is not guaranteed in the presence of measurement noise at the quantizer's input.

The weighted least-squares (WLS) method we introduce in the present paper bears some likeness to Wigren's, in the sense that it is also based on the minimization of a least-squares criterion. Its distinguishing features are that no approximation of the quantizer is made and that convergence can be guaranteed, even in the presence of measurement noise. Furthermore, this approach is more general than [START_REF] Wang | System identification using binary sensors[END_REF] in the sense that it does not require the presence of noise at the quantizer input.

The outline of the article is the following. In section II, the framework and the notations are introduced. In section III, we show how to construct exact WLS continuous criteria for parameter estimation from binary data. The asymptotic properties (when the number of data goes to infinity) of one of these criteria are established in section IV. Using probabilistic arguments, it is shown that this approach leads to a consistent estimate of the nominal parameters of the system, even in the presence of measurement noise. Section V contains some concluding remarks.

II. FRAMEWORK AND NOTATIONS

Let us consider a discrete-time linear time-invariant (LTI) system H . We assume H has a finite impulse response of length L , i.e. the impulse response can be represented by a column vector

( ) L k k 1 = = θ θ .
Let l u be the known scalar value of the system input at time l and l y be the (scalar) value of the system output, so that:

θ φ T L l l y , = , (1)
where

( ) l L l k k L l u 1 , + - = = φ
is the L sample-long (column) regression vector at time l . Let l d be an additive dithering signal at the system output and l b an additive noise (in other words, l d is known whereas l b is not). Let x m and 2 x σ denote the first-and second-order moments of any signal x . The system output is measured via a 1-bit ADC so that only the sign ( )

l l z S s =
of the system output is known, where ( ) ( ) designates:

   - = ≥ = otherwise , 1 0 if , 1 x S x x S , (2)
• θ padded with L L - ˆ zeros if L L > ˆ, • θ truncated to length L ˆ, if L L < ˆ.

III. DEFINITION OF WLS CONTINUOUS CRITERIA

In the following sections, we concentrate on criteria of the form:

( ) ∑ ∑ = = - = N l p l N l l l p l N p z s s z J 1 2 1 2 2 ˆ4 1 , 1 ≥ p . (3)
Note that taking 0 = p in (3) yields

( ) ∑ = - = N l l l N s s N J 1 2 0 41 ,
the "classical" unweighted least-squares criterion. However, because of the quantization of l s and l s ˆ, N J 0 is piecewise- constant with respect to θ ˆ (and, thus, not easily amenable to optimization). In the numerator of (3), the term

θ l z is continuously differentiable everywhere, ( ) θ N p J
is also continuously differentiable everywhere, except when:

0 θ φ φ =           +           ⇔ ≤ ≤ = ⇔ = ∑ = N T L N T L l N l p l d d N l z z M M 1 ,, 1 1 2 1 , 0 0 ˆ. ( 4
) Provided L N >
and under some mild assumptions on l u , there usually exists no θ ˆ that verifies (4) 3 . The major exception is when 0

= l d , N l ≤ ≤ 1
, in which case (4) has

0 θ = ˆ for solution, regardless of l u 4 . Note that 1 0 ≤ ≤ N p J .
Thus, it is possible to construct continuously differentiable WLS criteria for parameter estimation problems based on binary observations without having to rely on an approximation of the quantizer. These WLS criteria are "exact" in the sense that 0 0

0 = ⇒ = N p N J J
5 and they can be minimized using any standard optimization method [START_REF] Ljung | System identification -theory for the user[END_REF]. Note that it is impossible to guarantee their convexity a priori, without further assumptions on the inputs.

In [START_REF] Bai | Towards identification of Wiener systems with the least amount of a priori information on the nonlinearity[END_REF], Bai shows that minimizing N J 0 yields a consistent estimation of θ , supposing u is white, 0

= = d b , L L = and N goes to infinity. Since 0 0 0 = ⇒ = N p N J J
, it follows that the same is true of N p J . In the next section, we establish asymptotic properties of N J 1 under different assumptions, some less restrictive than in [START_REF] Bai | Towards identification of Wiener systems with the least amount of a priori information on the nonlinearity[END_REF] (

L L ≠ , 0 ≠ b , 0 ≠ d
, coloured input), some not (Gaussianity of z ˆ and z ).

IV. ASYMPTOTIC PROPERTIES OF ( )

θ 1N J A. Asymptotic expression of ( ) θ 1N J
In order to perform the asymptotic analysis of N J 1 , some assumptions must be made. A1l u , l d and l b are stationary ergodic Gaussian processes (not necessarily white) with mean 0. N J 1 can then be interpreted as the ratio of empirical means of stationary ergodic processes:

( ) ∑ ∑ = = - = N l l N l l l l N z N s s z N J 1 2 1 2 2 1 1 41 , ( 5
)
3 In a deterministic context, it suffices to design the input so that the first 1 ˆ+ L lines of (4) are independent. In a stochastic context, one can find simple conditions under which (4) has no solution with probability 1. 4 Note that these properties also hold when one replaces

p l z 2 ˆ by ( ) l z f ˆ in (3)
, where f is any positive even continuously-differentiable function such that ( ) 0 0 = f . 5 If l z ˆ is a continuous random variable, one can also prove that

0 0 0 = ⇒ = N N p J J with probability 1. Thus: ( ) ( ) ( ) 2 2 2 1 1 Ê Ê 4 1 lim z s s z J J N N - = = ∞ → ∞ . ( 6
)
From A1, z ˆ and z are jointly Gaussian with mean 06 . One may then write:

( ) ( ) ( ) ( ) ( ) ( )         + = - = ∫ ∫ ∫ ∫ ∫ ∫ ∞ ∞ - ∞ - ∞ ∞ ∞ - ∞ ∞ - ∞ 0 0 2 0 0 2 2 ˆ2 2 2 1 , , 4 1 , 4 1 dz
z d z z p z dz z d z z p z dz z d z z p z S z S z J z z σ σ , (7) ( ) ( ) 2
ˆ2 2 2 ˆ2 2 1 2 2 1 2 1 exp , ˆr z z r z z r z z p z z z z z z -                 - + - - = σ πσ σ σ σ σ , ( 8
)
where r is the cross-correlation coefficient of z and z ˆ. (7) can be expressed as:

( )

      - - = ∞ 2 1 1 acos 1 r r r J π . ( 9
)
The same reasoning applies to the classical least-squares criterion. This yields:

( )

r J J N N acos 1 lim 0 0 π = = ∞ → ∞ . (10) ∞ 1
J is a decreasing convex function of r ,

J is decreasing and concave on the same interval (Fig. 2). Theorem -If A1, A2, A3 and A4 hold, the following properties are true:

J 1 ∞ (r) J 0 ∞ (r)
Property 1 -If L L ≥ ˆ, a) if 0 2 ≠ d σ , the minimum of ∞ 1 J is unique and it is reached for L θ θ = . b) if 0 2 = d σ , ∞ 1 J has infinitely many minima of the form L θ θ κ = , ] [ +∞ ∈ , 0 κ . Property 2 -If L L < ˆ
and u is white, Property 1 still holds.

Property 3 - ( ) θ 1∞ J
is convex in a neighbourhood of its minimum (minima).

Property 4 -Properties 1 and 2 hold for ( )

θ 0∞ J .
Proof -The cross-correlation coefficient of two centred variables is ( ) ( ) ( )

z z z z r R θ θ 1 ˆ2 2 2 2 L L u T d L L u T d LL u T d b r + + + + = σ σ σ σ θ 0 θ 2 2 L L u L L u T d L L u L L u T d r + = + ⇔ = ∂ ∂ σ σ . ( 12
)
Thus, provided A4 holds and

0 2 ≠ + θ R θ L L u T d σ , opt L L u L L u L L u T d L L d L L u T d + + = σ σ κ and θ R R θ L L u L L u opt 1 ˆ- = . ( 13
)
Replacing θ ˆ by opt θ κ in the expression of κ yields:

( )

0 1 2 = - d σ κ . ( 14
)
One can then distinguish two cases, depending on whether

0 = d σ . If 0 2 ≠ d σ , then 1 = κ . To determine whether opt θ ˆ is a maximum, we show that ( ) θ r is concave in the neighbourhood of opt θ ˆ. Since ( ) θ r is at least 2 C on L R , we
can determine its Hessian matrix H :

( ) ( ) ( ) ( ) ( )                 + + +         + + + - = L L u T L L u L L u T d L L u T d L L u L L u T d L L u T L L u L L u T L L u g θ θ 1 λ - = + , 0 > n λ , ( ) ( ) ( ) [ ] ( ) ( ) ( ) ( )
         = - = = - = - + - = ∑ ∑ ∑ ∑ = = = = N l l L l N l l l l L l N l T L l L l N l l l T L l L l N N T d s s d s s J J 1 , 1 2 , 1 , , 1 2 , , 1 1 ˆˆ4 4 1 2 1 ˆφ q φ θ p φ φ Q φ φ θ P q p θ Q P θ Q θ θ g . ( 18)
A line-search method can be used to find an optimal value of the step parameter n λ . Because of its simplicity, we have called the identification method based on the minimization of N J 1 "BIMBO", for Basic Identification Method using Binary Observations. We show in Fig. 3

V. CONCLUSION

In this paper, we have introduced general WLS criteria for parameter estimation problems based on binary observations. Some asymptotical properties of this criterion were established in the simple case when 1 = p (BIMBO method), under assumptions of Gaussianity. In particular, we have proved that this approach leads to consistent estimates of the system parameters, regardless of the variance of the measurement noise at the quantizer's input or of a dithering signal. Moreover, we have established that the corresponding criterion is locally convex, which makes it amenable to simple optimization tools. The extension of this approach to less simple systems (IIR or nonlinear) or to a framework of online identification is also being studied.

Fig. 1 -

 1 Fig. 1 -Framework and notations.

Fig. 1 .

 1 Fig. 1. Finally, the notation L θ

  function of θ ˆ and so is N

Fig. 2

 2 Fig. 2 -Comparison of

  Fig. 3 -Contour lines of 256 0 J

The approach in[START_REF] Juillard | Initialization of the BIMBO self-test method using binary inputs and outputs[END_REF] was actually developed as a "quick and dirty" way to initialize the optimization algorithms on which the method presented in section III relies. However, because of its minimal requirements, it can also be used as a standalone parameter estimation approach.

From the point of view of electronics, this randomly specified threshold can be implemented as an addition of a dithering signal at the input of the (fixed-threshold) quantizer.

In practice, the joint Gaussianity of z ˆ and z may be approximately achieved with less stringent assumptions on the distribution of the system input.

We find:

is maximal and positive when

To pursue the analysis further, one can distinguish two more cases: