Souleymane Laminou Abdou
email: souleymane.laminou-abdou@univ-rennes1.fr

Pr Carl Chiarella

Carole Bernard

Grégoire Leblon

Rivo Randrianarivony

Jean-Laurent Viviani

Pricing and Hedging American and Hybrid Strangles with Finite Maturity *

Keywords: JEL classification: G12, G13, C61, C65 American strangles, early exercise boundaries, Kim representation, numerical integration, Greek parameters

This paper introduces variants of strangles, called Euro-American or hybrid strangles, and it promotes a new numerical pricing technique. We highlight and compare the properties of European, American, and hybrid strangles with pricing and hedging in mind. The new quadrature approach we propose can account for systems of coupled integral equations that locate the early exercise boundaries of finite-lived contracts. We show that this method is efficient, accurate, and fast for pricing all types of early exercisable strangles. Other advantages of this technique are that it avoids the nonmonotonic gradient problem faced by others and it allows users to control for errors. We then investigate the hedging of all strangles, we derive analytical expressions for some Greek parameters, and we stress how these parameters can differ (or not) from each other.

Introduction

According to [START_REF] Chaput | Volatility trade design[END_REF], strangles and straddles represent about 80% of option strategies. These strategies are used for risk management, volatility trading, and volatility speculation (see §11.4 of [START_REF] Hull | Options, Futures and Other Derivatives[END_REF]). The classical strangle is a European-style strangle and it is comprises a long position in a European put option and a long position in a European call option, with both options being written on the same underlying asset and maturing at the same time. The call option strike is typically greater than the put option one, but in the case where they are equal, the position is termed a straddle. Strangle positions have recently been studied in terms of the American style, meaning that the holder can decide to exercise earlier than maturity (see [START_REF] Gerber | Martingale approach to pricing perpetual American options[END_REF] and [START_REF] Chiarella | Evaluation of American strangles[END_REF], among others). An American strangle may be roughly viewed as a long position in an American put and an American call option equipped with a non-standard early exercise boundary (EEB) and a self-closing mechanism. This latter mechanism ensures that the right to sell (or to buy) disappears as soon as the holder decides to exercise the right to buy (or to sell)1 . As such, American strangles are not simple portfolios of standard American puts and calls2 .

In this study, we extend the family of strangles by introducing a new variant called Euro-American or hybrid strangles. These contracts can be exercised earlier than maturity given that only one side, the call side or the put side, is exercisable before expiration. Due to this possible early extinction, Euro-American contracts greatly differ from simple portfolios made up of independent American and European options. Euro-American strangles have several clear goals. From a practical viewpoint, these contracts can fit investors' and speculators' needs better than other strangles. Hybrids of American and European strangles can indeed offer the best of each of these variants. From a theoretical perspective, they can be used to understand and model real financial decisions. With real options in mind, strangle positions can model a firm decision to expand (exercise the call side) or to transfer (exercise the put side). In many cases, not all of these decisions can occur earlier than maturity, and if one side is delayed to maturity we have a hybrid strangle contract. From an option theory point of view, they contribute to an understanding of how American strangles work and how the opportunity for early exercise makes them different from standard European strangles. Hybrid contracts involve only a "one side" early exercise so we can assess the relative importance of it. Because a Euro-American strangle lies contractually between an American strangle and a European strangle, we can expect its price to be bounded by these two strangles. However, it is not a priori straightforward to predict how much cheaper or more expensive Euro-American strangles are compared to other strangles.

Pricing American-style contracts is known to be a challenge because the holder can exercise at any time before maturity. Two interconnected questions must be answered simultaneously: what is the best time to exercise and what is the resulting payoff? Mathematically, this free boundary problem can be addressed by solving a partial differential equation subject to some (boundary) conditions (see [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF] for an early treatment) 3 . Pricing American strangles can appear even more challenging because they depend on two interdependent, selfclosing, and time-varying early exercise boundaries. For their part, Euro-American strangles depend on a single time-varying early exercise boundary only, but this is subtly influenced by the existing one-sided European feature of the contract, and furthermore, the opportunity to exercise at maturity disappears if an early exercise is decided beforehand. Such contracts have really special features to consider.

Past studies on the pricing of American strangles with finite maturity advocated the use of various advanced numerical techniques. [START_REF] Alobaidi | Laplace transforms and the American straddle[END_REF] used Laplace transforms and derived analytical formulas for pricing American straddles as well as integral equations to locate early exercise boundaries. Unfortunately, the expressions they provide cannot be inverted analytically and they give no recommendation for numerical inversion. [START_REF] Chiarella | Evaluation of American strangles[END_REF](CZ) used Fourier transforms and derived analytical formulas for finite-lived American strangles as well as a system of coupled integral equations to locate early exercise boundaries. Finally, they employed a two-step algorithm for pricing, mixing a quadrature approach and an interpolation technique. This two-step approach seemed necessary to deal 3 Various approaches have been developed for pricing American-style options. [START_REF] Macmillan | An analytical approximation for the American put prices[END_REF] provided an analytic approximation for valuing the early exercise premium, which is the price difference between an American option and its European equivalent. [START_REF] Kim | The analytic valuation of American options[END_REF] later offered an intuitive representation of the early exercise premium (see also [START_REF] Jacka | Optimal stopping and the American put[END_REF] and [START_REF] Carr | Alternative characterizations of American put options[END_REF]). However, this representation requires the knowledge of the early exercise boundary at any point in time and consequently relies on a computational technique to solve the integral equation locating the boundary. Another way to proceed is to use numerical methods for solving the free boundary problem. These methods are now well known in finance. Finite difference-based methods are the most common (see, for instance, §28 of [START_REF] Wilmott | Paul Wilmott introduces quantitative finance[END_REF]) and Monte Carlo simulations are also used (see [START_REF] Broadie | Enhanced monte carlo estimates for American option prices[END_REF] and [START_REF] Longstaff | Valuing American options by simulation: A simple least-squares approach[END_REF] with the non-monotonic gradient problem they faced4,5 . To avoid the above advanced methods, one may be tempted to consider simpler portfolios made up of two individual American call and put options. In some cases the pricing bias is indeed rather limited6 , but in general, this way to proceed is hazardous because it can lead to suboptimal early exercise decisions (see [START_REF] Moraux | On perpetual American strangles[END_REF] for examples in the perpetual case).

In this study, we derive analytical formulas for all finite-lived strangles (among which are the new variant contracts) and promote a new numerical approach able to deal with the various systems of integral equations to locate early exercise boundaries. The quadrature method we introduce combines some Newton-Cotes weights and some fourth-order Gregory weights, as presented in Linz (1985, p.98) and Press, Teukolsky, Veterling, & Flannery (2007, pp.159-160). This numerical approach in turn has several interesting features. First, it is a rather simple (one-step) numerical approach. Second, it is accurate, efficient, and fast. Third, the approach allows users to control for errors, and we show that numerical estimates tend to the true price from above as the number of discrete points increases7 . Finally, the quadrature scheme faces no non-monotonic gradient problem and consequently requires neither repeated computations of early exercise boundaries nor extrapolation techniques.

We expect the prices of Euro-American strangles to lie between those of European and American comparable strangles. Indeed, our simulations reveal that Euro-American strangles can effectively be more expensive than comparable European strangles and cheaper than comparable American strangles, but we also provide scenarios where Euro-American strangles are as expensive as their European or American counterparts. Hence, the opportunities to buy or to sell the underlying asset earlier than the expiration date may be, in some contexts, effectively worthless. This information should be useful for potential users of existing early exercisable strangles because they should not pay for such opportunities.

Beyond introducing new variants of early exercisable strangles, our (financial) analysis differs from that of CZ. We split the price of strangles following the early exercise premium (EEP) representation advocated by [START_REF] Kim | The analytic valuation of American options[END_REF] 8 . Hence, in the present study, the American strangle is essentially the sum of a European strangle plus a premium to exercise it earlier than maturity. By contrast, CZ price the American strangle as a whole and then represent the contract as a portfolio made up of an American call and an American put with adjusted early exercise boundaries9 (see their Proposition 7). Logically, most of their simulations compare early exercise boundaries of strangles to those of standard American options. Instead, we think it is important to consider European strangles as relevant benchmarks because we know their price analytically. We can then explore the hedging parameters of every strangle and highlight how they differ from each other.

The rest of the paper is organized as follows. Section 2 presents the framework and problem statement. Section 3 presents the early exercise premium representation of Americanstyle strangles. Section 4 presents and discusses the new numerical approach. Section 5 investigates the prices of American-style strangles with finite maturity. Section 6 considers the hedging issues and Greek parameters. The last section (Section 7) presents our conclusions.

The framework and problem statement

This section introduces contract specifications, notations, and hypotheses. The financial markets in our setting are perfect, efficient, and complete; trading takes place continuously and information is free. There are neither taxes nor transaction costs. There is a risk-free asset paying a known and constant interest rate denoted by r. There is also a risky asset, say a stock, paying a continuous dividend rate δ that underlies the different contracts. When the stock pays dividends, we know from standard option theory that it may be optimal to exercise standard American call options earlier than at expiration. The risk neutral price process of the underlying asset is denoted by S = (S t ) t≥0 and it is described by

dS t = (r -δ) S t dt + σS t dW t , S 0 = x
where W = (W t ) t≥0 is a standard Brownian motion, σ is the volatility, and x ≥ 0. It follows that S t = x exp [(r -δ -σ 2 /2) t + σW t ]. Throughout the paper, the normal cumulative density function is denoted by N and the normal probability density function is denoted by ϕ.

A strangle (K 1 , K 2 , T ) is a contract expiring at time T that gives the holder the right to sell the underlying asset S at price K 1 and the right to buy at price K 2 (we omit S in previous parentheses because there is no ambiguity). Hereafter, we use the terms "call side" and "upper side" interchangeably for describing the right to buy embedded in the strangle.

Similarly, the terms "put side" and "lower side" are used interchangeably for describing the right to sell. K 1 may therefore be termed the "put side" strike of the strangle and K 2 the "call side" strike of the strangle. If ever K 1 = K 2 , the contract is a straddle. It should be noted that in this paper we do not discuss "American guts," for which K 1 > K 2 . European strangles allow the holder(s) to sell the underlying asset at K 1 or to buy it at K 2 at the expiration of the contract only. American strangles allow the holder(s) to act at or before the expiration; that is, at any time between inception and termination10 . The Euro-American or hybrid strangles that we introduce in this paper allow the holder(s) to sell or buy at the expiration of the contract and provide an opportunity to act earlier than at expiration on one side only. Such contracts can, for instance, allow the holder(s) to buy at any time before T and to buy or sell at time T , or the converse (i.e., to sell at any time before T and to buy or sell at time T ), but not both.

It is well known (see for instance [START_REF] Merton | Theory of rational option pricing[END_REF]) that the price V (S, t) of every contract written on S satisfies the following fundamental partial differential equation (PDE)

∂V ∂t + (r -δ) S ∂V ∂S + 1 2 σ 2 S 2 ∂ 2 V ∂S 2 -rV = 0, t ∈ [0, T ). (1)
Of course, the contract considered here is assumed to be alive. For early exercisable contracts, this means that the time t-price of the underlying asset lies in the continuation region. Figure 1 gives a typical representation of the price of a classical American strangle that highlights the intrinsic value (which is max (K 1 -S; 0) + max (S -K 2 ; 0)) and corresponding continuation and stopping regions. These regions are characterized by the early exercise boundaries that relate to the threshold values of the underlying asset where the decision to exercise early should intervene (because it is optimal for holders to do so). In Figure 1, we emphasize these threshold values with dotted lines. For Euro-American strangles with an early exercisable call side feature (or, respectively, an early exercisable put side feature), the stopping region on the left (on the right) does not exist.

"INSERT FIGURE 1 ABOUT HERE"

Let us now denote by τ = T -t the time to expiration, by A c 1 c 2 (S, τ ) the price of the American strangle, by c 1 (τ ) the lower early exercise boundary associated with the right to sell, and by c 2 (τ ) the upper early exercise boundary associated with the right to buy. Due to the exercise policy, finite-lived American-style contract boundaries are time dependent (see Merton (1973, pp.170-171) and depend on the remaining time to expiration. As long as the underlying stock price lies between these two thresholds, it is not optimal to exercise the American strangle (cf. the continuation region in Figure 1). Consequently, the price of the American strangle satisfies the fundamental equation (1) with a final condition

A c 1 c 2 (S T , 0) = max [0, (K 1 -S T ) ψ(S T -c 1 (0))] + max [0, (S T -K 2 ) ψ (c 2 (0) -S T )] , S T ≥ 0, (2a)
and some specific boundary conditions

lim S↓c 1 (τ ) A c 1 c 2 (S, τ ) = K 1 -c 1 (τ ) , 0 ≤ τ ≤ T, ( 2b
) lim S↑c 2 (τ ) A c 1 c 2 (S, τ ) = c 2 (τ ) -K 2 , 0 ≤ τ ≤ T, (2c) lim S↓c 1 (τ ) ∂A c 1 c 2 (S, τ ) ∂S = -1, lim S↑c 2 (τ ) ∂A c 1 c 2 (S, τ ) ∂S = 1, 0 ≤ τ ≤ T. ( 2d
)
where ψ (x) is the Heaviside step function defined by ψ (x) = 0 if x ≤ 0 and ψ (x) = 1 if

x > 0. Equation (2a) depicts the value of American strangles at expiration in the case where no early exercise occurs. At expiration, a rational investor exercises the in-the-money side of the American strangle. If ever the price of the underlying asset first reaches the early exercise boundary c 1 (τ ) at a given time τ before expiration, then the holder of the American strangle immediately exercises their right to sell the stock and the received payoff is Equation (2b). If ever the underlying asset first reaches the early exercise boundary c 2 (τ ) at a given time τ before expiration, then the holder of the American strangle immediately exercises their right to buy the stock and the received payoff is as shown by Equation (2c). Equations ( 2b)-(2c) represent value-matching conditions and Equations (2d) represent smooth-pasting conditions. They can be deduced from Merton (1973, pp.170-171) or Chiarella & Ziogas (2005, p.35). Notice that including step functions in the terminal payoff above is important to limit cash flows at maturity to K 1 -c 1 (0) or c 2 (0)-K 2 , where c 1 (0) and c 2 (0) stand for the limit of early exercise boundaries as τ tends to zero. It is impossible for the underlying asset price to end beyond the boundaries without crossing them beforehand and hence provoking an early exercise.

Euro-American strangles are new variants introduced in this study. Let us first consider Euro-American strangles with an early exercisable call side. These contracts offer the right to sell or buy the underlying asset at maturity and the right to buy it earlier. Similar to standard American call options, these contracts imply, in the presence of dividends, an upper early exercise boundary (placed on the price of the underlying asset) whose value at time t can be denoted by h (τ ) with τ = T -t. Holders of these contracts should exercise early if and when the price of the underlying asset reaches the boundary h. As long as the underlying stock price lies below the threshold h, it is not optimal for holders to exercise the right (to buy) earlier than at expiration. Consequently, the price of Euro-American strangles with an early exercisable call side, denoted by EA h (S, τ ), satisfies the fundamental Equation (1) with the following final time and boundary conditions:

EA h (S T , 0) = max [0, K 1 -S T ] + max [0, (S T -K 2 ) ψ (h (0) -S T )] , S T ≥ 0, (3a) lim S↑c 2 (τ ) EA h (S, τ ) = h(τ ) -K 2 , 0 ≤ τ ≤ T, ( 3b
) lim S↓0 EA h (S, τ ) = 0, 0 ≤ τ ≤ T, (3c) lim S↑c 2 (τ ) ∂EA h (S, τ ) ∂S = 1, 0 ≤ τ ≤ T. ( 3d
)
Let us now consider Euro-American strangles with an early exercisable put side. These contracts offer the right to sell or buy the underlying asset at maturity and the right to sell it earlier. This hybrid contract implies a lower early exercise boundary (placed on the price of the underlying asset) whose value at time t is denoted by l (τ ) with τ = T -t. Holders of these contracts should exercise early if and when the price of the underlying asset reaches the boundary l. As long as the underlying stock price lies above the threshold l it is not optimal for holders to exercise the right (to sell) earlier than expiration. In view of this, the price of Euro-American strangles with an early exercisable put side, denoted by EA l (S, τ ), satisfies the fundamental Equation (1) with the following final time and boundary conditions:

EA l (S T , 0) = max [0, (K 1 -S T ) ψ(S T -l (0))] + max [0, (S T -K 2 )] , S T ≥ 0, (4a) lim S↓l(τ ) EA l (S, τ ) = K 1 -l (τ ) , 0 ≤ τ ≤ T, ( 4b
)
lim S↑∞ EA l (S, τ ) = 0, 0 ≤ τ ≤ T, (4c) lim S↓c 1 (τ ) ∂EA l (S, τ ) ∂S = -1, 0 ≤ τ ≤ T. ( 4d
)
3. Early exercise premium representation [START_REF] Kim | The analytic valuation of American options[END_REF], [START_REF] Jacka | Optimal stopping and the American put[END_REF][START_REF] Carr | Alternative characterizations of American put options[END_REF] showed that the price of standard Amer-ican options can be decomposed into the price of a European option plus an early exercise premium having a special but intuitive integral representation. These authors demonstrated the EEP representation by different means. Hereafter, we provide EEP representations for all American and hybrid strangles. Such a representation strategy appears especially useful for early exercisable strangles for at least two reasons. First, we can price European strangles analytically. Second, European strangles can then play the role of common benchmarks for all contracts11 . Proposition 1 splits the price of an American strangle (A c 1 c 2 (S, τ )) into the price of a European strangle (E(S, τ )) plus an early exercise premium (EEP c 1 c 2 (S, τ )) with an integral representation.

Proposition 1. The price of an American strangle is given by:

A c 1 c 2 (S, τ ) = E (S, τ ) + EEP c 1 c 2 (S, τ ) (5)
where

E(S, τ ) = K 1 e -rτ N (-d 2 (S, τ ; K 1 )) -Se -δτ N (-d 1 (S, τ ; K 1 )) + Se -δτ N (d 1 (S, τ ; K 2 )) -K 2 e -rτ N (d 2 (S, τ ; K 2 ))
and

EEP c 1 c 2 (S, τ ) = ∫ τ 0 [ K 1 re -r(τ -η) N (-d 2 (S, τ -η; c 1 (η))) -Sδe -δ(τ -η) N (-d 1 (S, τ -η; c 1 (η))) + Sδe -δ(τ -η) N (d 1 (S, τ -η; c 2 (η))) -K 2 re -r(τ -η) N (d 2 (S, τ -η; c 2 (η))) ] dη
with early exercise boundaries c 1 (.) and c 2 (.) defined by

c 2 (τ ) -K 2 = A c 1 c 2 (c 2 (τ ) , τ ) (6) K 1 -c 1 (τ ) = A c 1 c 2 (c 1 (τ ) , τ ) (7)
and

d 1 (S, τ ; β) = ln( S β )+ r-δ+ σ 2 2 τ σ √ τ , d 2 (S, τ ; β) = d 1 (S, τ ; β) -σ √ τ .
Proof. To prove Proposition 1, we follow [START_REF] Kim | The analytic valuation of American options[END_REF] and divide the pricing problem into two parts. We can then solve these two sub-problems independently using the results of [START_REF] Kolodner | Free boundary problem for the heat equation with applications to problems of change of phase. i. general method of solution[END_REF] and [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF] directly.

Hence, the main difference with CZ is that we have recourse to the techniques introduced by [START_REF] Kim | The analytic valuation of American options[END_REF] to deal with the standard American options. We split the pricing of strangles into two problems that we solve by using the results of [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF] and [START_REF] Kolodner | Free boundary problem for the heat equation with applications to problems of change of phase. i. general method of solution[END_REF]. We end up with the European component of the considered strangle and its EEP. By contrast, CZ consider the problem as a whole and apply a Fourier transform approach on the whole. They transform the PDE and apply the incomplete Fourier transform approach used by McKean to solve the free boundary problem. They decide to reorganize, to highlight the sum of an American-style put option and an American-style call option (with non-standard and adjusted EEPs). Of course our Proposition 1 conforms to their results.

The expression for the early exercise premium is semi-analytical only because it depends on a couple of early exercise boundaries to be determined. To this end, Equations ( 6) and ( 7) form a coupled integral equation system that we need to solve numerically by some techniques (see next section for details).

Proposition 2 details our core result on Euro-American strangles. The lengthy proof (relegated to the Appendix) uses the results of [START_REF] Kolodner | Free boundary problem for the heat equation with applications to problems of change of phase. i. general method of solution[END_REF] and [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF] on free boundary problems.

Proposition 2. i) The price of a Euro-American strangle with an American call side is given by:

EA h (S, τ ) = E (S, τ ) + EEP h (S, τ ) (8)
where

EEP h (S, τ ) = ∫ τ 0 [Sδe -δ(τ -η) N (d 1 (S, τ -η; h(η))) -K 2 re -r(τ -η) N (d 2 (S, τ -η; h(η)))
] dη the early exercise boundary h is given by

h (τ ) -K 2 = EA h (h (τ ) , τ ). (9)
ii) The Price of a Euro-American strangle with an American put side is given by:

EA l (S, τ ) = E (S, τ ) + EEP l (S, τ ) (10)
where

EEP l (S, τ ) = ∫ τ 0 [ K 1 re -r(τ -η) N (-d 2 (S, τ -η; l(η))) -Sδe -δ(τ -η) N (-d 1 (S, τ -η; l(η)))]dη
and the early exercise boundary is given by

K 1 -l (τ ) = EA l (l (τ ) , τ ). ( 11
)
Proof. See Appendix A.

The next step is to compute the Volterra integrals involved in Equations ( 5), ( 8), and (10).

Numerical implementation and algorithmic issues

We will use a quadrature approach to evaluate the integrals of Equation ( 5) and to transform Equations ( 6) and ( 7). A quadrature is a way to approximate an integral in general. If f is a real function of, say, one variable defined on [a, b], then the integral of f over [a, b] may be approximated by

∫ b a f (x)dx ≈ h n ∑ k=0 w k f (x k )
where (x k ) k=1,...,n are equally spaced sample points such that x 0 = a, x k = x 0 + kh, and

x n = b. The interval [a, b] has been divided here into n subintervals of equal width (with h = (b -a)/n). Quadrature methods differ essentially in terms of weights (w k ); for instance, CZ use Cavalieri-Simpson weights (see Table 1).

The Cavalieri-Simpson approach has significant drawbacks when applied to American strangles. It requires an even number of intervals, as highlighted by Table 1, leads to nonmonotonic gradients, requires a double computation of each early boundary, and then needs a technique to average both estimates. CZ suggest the use of Richardson's extrapolation technique to rectify this and to combine the two boundaries (found for h and h/2). Nevertheless, there is no possibility here to control for errors. To address this issue, we follow another way and propose a new quadrature method that is able to compute in the one-step integrals involved in the American and Euro-American strangles. We will see that this alternative way is accurate, efficient, and faster.

The new one-step quadrature method that we propose is as follows. Consider first a number of points n to approximate the EEB along the time line and i an index position taking values in {1, 2, ..., n}. The points of the EEB(s) are then obtained by computing integral(s) iteratively at i for i equal to 1 up to n. Each i is associated with a point of the EEB. i = 1 points to the first position before expiration (where we want to compute the EEB(s)). At each position i, integral equations are computed using i intervals and (i + 1) weighted endpoints. For weights, we use the closed-form formulas of Newton-Cotes when i is smaller than or equal to six (listed in Table 2). For i greater than six, we use the fourth-order Gregory quadrature formulas12 (derived by fitting cubic polynomials through successive groups of four points; see [START_REF] Linz | Analytical and numerical methods for Volterra equations[END_REF] and [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]). Unlike Simpson's quadrature, this quadrature method can deal with odd and even values of intervals. The general form of Gregory weights is

∫ b a f (x)dx ≈ h [ 3 8 f (x 0 ) + 7 6 f (x 1 ) + 23 24 f (x 2 ) + f (x 3 ) + f (x 4 )... +f (x n-4 ) + f (x n-3 ) + 23 24 f (x n-2 ) + 7 6 f (x n-1 ) + 3 8 f (x n ) ]
.

"INSERT TABLE 1 AND 2 ABOUT HERE"

In terms of convergence, a key advantage of this quadrature is that we have analytical expressions for computing errors for each value of i. Consequently, we can control for errors. For the Newton-Cotes part, the analytical expressions are as follows. For even values of i and any function f ∈ C (i+2) ([a, b]), the quadrature error is

E i (f ) = M i (i + 2)! h i+3 f (i+2) (ξ)
where ξ ∈ [a, b] and

M i = ∫ i 0 t ( i ∏ k=0 (t -k) ) dt < 0.
For odd values of i and any function

f ∈ C (i+1) ([a, b]), the quadrature error is E i (f ) = K i (i + 1)! h i+2 f (i+1) (η)
where η ∈ [a, b] and

K i = ∫ i 0 i ∏ k=0 (t -k) dt < 0.
The degree of exactness (DE hereafter) is equal to i + 1 for even values of i and i for other values. For the Gregory part, the analytical expressions are quite lengthy and we refer to Linz (1985, pp.100-107), where a detailed and complete error analysis is provided. Table 3 summarizes the numerical values of M i /(i + 2)! and K i /(i + 1)! and the DEs for the Newton-Cotes part13 . It should be noted that from i = 2, the integration order is higher or equal to that of the extended Simpson's quadrature14 . The DE of the fourthorder Gregory formulas is three. As emphasized earlier, our quadrature approach is suitable for even and odd values of i and only one step is necessary whatever the value of i. By comparison, the usual interpolation techniques implement the same algorithm twice (a first step for odd values of i and a second step for even values) and they then compute a weighted average of both. Consequently, our approach leads to less calculus and so should be less time consuming (this is verified below in Table 4.). A pseudo-code for the hybrid quadrature is presented in Appendix C (see Algorithm 1).

Discretizing the time interval is the next step in a numerical quadrature strategy. Following [START_REF] Linz | Analytical and numerical methods for Volterra equations[END_REF] and [START_REF] Kim | The analytic valuation of American options[END_REF], we divide the time interval into n parts of length h and define τ i = ih for i = 1, 2, 3, . . . , n and h = T /n.

Before computing the price of an American or Euro-American strangle we must determine early exercise boundaries by iteratively solving integral equations for each value τ i . The initial values for c 1 , c 2 , h and l are c

1 (0) = l(0) = min ( K 1 , r δ K 1 ) and c 2 (0) = h(0) = max ( K 2 , r δ K 2 ) 15
. Then, for each step i, starting from i = 1, the only unknown variables in the system of integral equations are c 1 (τ i ), c 2 (τ i ), h(τ i ), and l(τ i ). Given that the values of c 1 (τ j ), c 2 (τ j ), h(τ j ), and l(τ j ) for j < i are known from previous steps, we have:

-For American strangles

K 1 -c 1 (ih) = E(c 1 (ih), ih) + EEP c 1 c 2 (c 1 (ih), ih) c 2 (ih) -K 2 = E(c 2 (ih), ih) + EEP c 1 c 2 (c 2 (ih), ih)
where

E(c k (ih), ih) = K 1 e -rih N (-d 2 (c k (ih), ih; K 1 )) -c k (ih)e -δih N (-d 1 (c k (ih), ih; K 1 )) + c k (ih)e -δih N (d 1 (c k (ih), ih; K 2 )) -K 2 e -rih N (d 2 (c k (ih), ih; K 2 ))
and

EEP c 1 c 2 (c k (ih), ih) = h i ∑ j=0 w j       K 1 re -rh(i-j) N (-d 2 (c k (ih), (i -j) h; c 1 (jh))) -c k (ih)δe -δh(i-j) N (-d 1 (c k (ih), (i -j) h; c 1 (jh))) +c k (ih)δe -δh(i-j) N (d 1 (c k (ih), (i -j) h; c 2 (jh))) -K 2 re -rh(i-j) N (d 2 (c k (ih), (i -j) h; c 2 (jh)))      
for k = 1, 2, given that the weights (w j ) j depend on the above hybrid scheme.

-For Euro-American strangles with an American call side

h(ih) -K 2 = E(h(ih), ih) + EEP h (h(ih), ih)
where

E(h(ih), ih) = K 1 e -rih N (-d 2 (h(ih), ih; K 1 )) -h(ih)e -δih N (-d 1 (h(ih), ih; K 1 )) + h(ih)e -δih N (d 1 (h(ih), ih; K 2 )) -K 2 e -rih N (d 2 (h(ih), ih; K 2 ))
and

EEP h (h(ih), ih) = h i ∑ j=0 w j [ h(ih)δe -δh(i-j) N (d 1 (h(ih), (i -j) h; h(jh))) -K 2 re -rh(i-j) N (d 2 (h(ih), (i -j) h; h(jh)))
]

-For Euro-American strangles with an American put side

K 1 -l(ih) = E(l(ih), ih) + EEP l (l(ih), ih)
where

E(l(ih), ih) = K 1 e -rih N (-d 2 (l(ih), ih; K 1 )) -l(ih)e -δih N (-d 1 (l(ih), ih; K 1 )) + l(ih)e -δih N (d 1 (l(ih), ih; K 2 )) -K 2 e -rih N (d 2 (l(ih), ih; K 2 ))
and

EEP l (l(ih), ih) = h i ∑ j=0 w j [ K 1 re -rh(i-j) N (-d 2 (l(ih), (i -j) h; l(jh))) -l(ih)δe -δh(i-j) N (-d 1 (l(ih), (i -j) h; l(jh))) ]
For d 1 we have the following limits:

lim ξ→τ N (d 1 (c 1 (τ ), τ -ξ; c 2 (ξ))) = lim ξ→τ N (-d 1 (c 2 (τ ), τ -ξ; c 1 (ξ))) = 0 lim ξ→τ N (d 1 (c 2 (τ ), τ -ξ; c 2 (ξ))) = lim ξ→τ N (-d 1 (c 1 (τ ), τ -ξ; c 1 (ξ))) = 0.5 lim ξ→τ N (d 1 (h(τ ), τ -ξ; h(ξ))) = lim ξ→τ N (-d 1 (l(τ ), τ -ξ; l(ξ))) = 0.5
These expressions are derived using the limits of the function d 1 (c u (τ ), τ -ξ; c v (ξ)) given by:

lim ξ→τ ln ( cu(τ ) cv(ξ) ) + ( r -δ + σ 2 2 ) (τ -ξ) σ √ τ -ξ =          0 if u = v ∞ if u > v -∞ if u < v for u, v ∈ {1, 2}.
From a financial point of view, it is notable that lim ξ→τ N (d 1 (c 1 (τ ), τ -ξ; c 2 (ξ))) and lim ξ→τ N (-d 1 (c 2 (τ ), τ -ξ; c 1 (ξ))) represent the probabilities that c 1 (τ ) > c 2 (ξ) and c 2 (τ ) > c 1 (ξ), respectively, within a short period of time τ -ξ. By definition, at any time until maturity, we have c 1 (.) < c 2 (.). So both probabilities are zero. Along similar lines,

lim ξ→τ N (d 1 (c 2 (τ ), τ -ξ; c 2 (ξ))), lim ξ→τ N (-d 1 (c 1 (τ ), τ -ξ; c 1 (ξ))), lim ξ→τ N (d 1 (h(τ ), τ -ξ; h(ξ))),
and lim ξ→τ N (-d 1 (l(τ ), τ -ξ; l(ξ))) represent the respective probabilities that c 2 (τ ) > c 2 (ξ),

c 1 (τ ) < c 1 (ξ), h(τ ) > h(ξ)
, and l(τ ) < l(ξ) within a short period of time τ -ξ. When ξ tends to τ , due to the property of the cumulative normal distribution function, the probabilities are equal to 0.5. In discrete time, the equations yield to:

lim j→i N (d 1 (c 1 (ih), (i -j) h; c 2 (jh))) = lim j→i N (-d 1 (c 2 (ih), (i -j) h; c 1 (jh))) = 0 lim j→i N (d 1 (c 2 (ih), (i -j) h; c 2 (jh))) = lim j→i N (-d 1 (c 1 (ih), (i -j) h; c 1 (jh))) = 0.5 lim j→i N (d 1 (h(ih), (i -j) h; h(jh))) = lim j→i N (-d 1 (l(ih), (i -j) h; l(jh))) = 0.5
Similar limits hold also for d 2 . As a result, when j reaches the value of i, the numerical implementation does not involve the boundaries values at jh and ensures the existence of d 1 and d 2 . To determine the unknown value of the boundaries at step i, we introduce two different functions:

f 1 (x) = -K 1 + x + A(x, ih) f 2 (x) = -x + K 2 + A(x, ih)
where A is the price of the considered strangle contract, either a hybrid American strangle or a standard one, at time to maturity ih. The roots of f 1 (or f 2 ) computed at i correspond exactly to the lower (or upper) early exercise boundaries at that time 16 .

To analyze the complexity of our approach, assume that there are m operations in the bisection method (m depends on the targeted level of accuracy). To find each boundary at step i, we have i + 1 evaluations of the integral equation. For n evaluations we have ((n + 1)(n + 2))/2 operations. The total number of operations (considering an evaluation of the integrand as one operation) is m[((n + 1)(n + 2))/2]; its order is mn 2 . This number highlights how the runtime of the code is influenced by the choice of n. This is the same as in [START_REF] Chiarella | Evaluation of American strangles[END_REF] only when the Richardson's extrapolation step is ignored.

Appendix C shows the explicit pseudo-codes for the hybrid quadrature, the determination of the EEB, and the price of a Euro-American strangle with an early exercisable call side 17 . Algorithms 2 and 3 are iteratively and recursively linked. Algorithm 2 involves the index i, whereas Algorithm 3 refers to n. We can see that at each step i the function HybridQuad(i) gives the quadrature weights for i + 1 endpoints.

Simulations

First, it is worth schematically comparing our method to that of CZ. Figure 2 shows how the numerical scheme of CZ and ours differ in their principles. n and 2n are the numbers of the discretization length of the integrals in the price formula.

Next, Figure 3 plots the call side early exercise boundary computed with the CZ numerical scheme before (checkered gray line) and after (black line) Richardson's extrapolation. The checkered gray line in this figure clearly shows the non-monotonic gradient problem. 16 The bisection method is used with a six-digit precision for the root searching result. The Newton-Raphson method, the bisection method, and other methods are available for root finding. By the nature of the contract, the holder of the option will choose the best of a put with strike K 1 and a call with strike K 2 (with a self-closing mechanism). Accordingly, the function A(x, ih) is a decreasing function of x within [0, K 1 ] and an increasing function of x within [K 2 , ∞) (see Figure 1). This monotonic property is extended to functions f 1 and f 2 within the root searching intervals and ensures that the bisection method succeeds at any arbitrary level of confidence.

17 Parameter sigma, the volatility, is implicitly used in functions d 1 and d 2 .

Interestingly, we can also notice a) that there is a need for an extrapolation technique and b) that the Richardson's extrapolation technique is far from being a simple averaging of external points. Figure 4 compares the call side early exercise boundary of CZ after the Richardson's extrapolation (gray line) with the one we obtain from our own procedure (black line). Both EEBs are associated with the same American strangle but both boundaries differ from each other. Our boundary appears to be lower than the one provided by the CZ algorithm. Consequently, our method suggests to exercise the call side of the strangle earlier and this increases the probability of exercising the call side. Due to the smoothness of the EEB pricing the accuracy should be improved in our context. It becomes clear that the choice of the quadrature method deserves particular attention and significantly affects the pricing and the management of American-Style strangles. Table 3, introduced below, will confirm these statements.

"INSERT FIGURES 2 AND 3 ABOUT HERE" "INSERT TABLE 3 ABOUT HERE"

Using the one-step procedure saves computational time because it is run with an interval length of h/2 only, whereas the two-step algorithm uses both h and h/2. So, by construction, CZ use 50% more subintervals than we do. To assess the gain in time, we proceed as follows. We implement both of the codes for the different values of n that are the number of h/2 long subintervals. The computational times are then displayed on the bottom line of Table 4. The table reveals that for a given n, the convergence towards the common long-run value is faster with our method. The computational time for the two-step procedure is at least about 25% longer. The maximum difference is obtained for n = 400, with a difference of 36%.

"INSERT TABLE 4 ABOUT HERE"

Table 4 compares the prices computed with the new quadrature method to those provided by CZ for n = 100, 200, 400, and 800. We consider different values of the stock price S to represent five different situations: in-the-money (on the call and put sides), at-the-money (on the call and put sides), and out-of-the-money (i.e., in between K 2 and K 2 ). The prices obtained by both methods share four decimal digits. When n increases, the fifth decimal converges to the same value. We can see that our algorithm converges faster than that of CZ in the sense that for an equivalent computational time our prices are closer to the exact solution than CZ's. Table 4 also shows that, as predicted, the price limit is reached from above when using our approach. By contrast, the CZ numerical approximation can converge from above or below (see for instance S = 1.5 and S = 1.75). The bottom line is that we can avoid the non-monotonic gradient problem, save computational time, improve the convergence rate, and control for errors.

"INSERT TABLE 5 ABOUT HERE"

Table 5 compares the prices of American and Euro-American strangles for three different environments (Panels A, B, and C) and for a range of asset prices chosen to highlight the price differences in all cases. We can check that the prices of the European and American strangles bound all other values. Panel A displays a situation where the risk-free rate is far smaller than the dividend rate. Here, the Euro-American strangles with American put sides are cheaper than their American strangle counterparts, whereas the Euro-American strangles with American call sides are as expensive. This panel highlights a situation where the lower early exercise boundary does not matter. In other words, the American strangle contract offers an opportunity to sell earlier than maturity that is worth nothing. Panel B provides simulations in an environment where the risk-free rate is far greater than the dividend rate. In that situation, the upper early exercise boundary does not matter, reflecting that an American strangle contract offers, in such a context, a worthless opportunity to buy earlier than maturity. Panel C highlights a situation where the risk-free rate is equal to the dividend rate. Here, the prices of both the Euro-American strangles are different from those of the standard American strangles.

Hedging management issues

It is important for traders and writers to hedge option positions. Dynamic hedging strategies suggest the need to consider price sensibilities. Delta (∆), gamma (Γ), and theta (Θ) are in these respects key parameters. Kim-style representations of American and Euro-American strangles within the Black-Scholes model are useful to find analytical formulas. Proofs for the following hedging parameters are in Appendix B. Delta, which is the first derivative of the position with respect to the stock price, measures the sensitivity of the contract price to the underlying stock price. This hedging ratio indicates the number of stocks to hold in a dynamic hedging strategy. From our previous representations of strangles, we have

∆ V = ∂V (S, τ ) ∂S = ∆ E + ∆ EEP (12)
where ∆ E is the delta of a European strangle (∆ E = e -δτ (N (d

1 (S, τ ; K 2 )) -N (-d 1 (S, τ ; K 1 ))))
and ∆ EEP is the delta of the early exercise premium. For the American strangle, we have

∆ AS EEP = ∫ τ 0 δe -δ(τ -η) {N (d 1 (S, τ -η; c 2 (η))) -N (-d 1 (S, τ -η; c 1 (η))) + 1 σ √ τ -η [ ϕ (d 1 (S, τ -η; c 2 (η))) ( 1 - K 2 c 2 (η) r δ ) +ϕ (-d 1 (S, τ -η; c 1 (η))) ( 1 - K 1 c 1 (η) r δ )]} dη.
For a Euro-American strangle with an American upper side we have

∆ EcAS EEP = ∫ τ 0 δe -δ(τ -η) {N (d 1 (S, τ -η; h(η))) + 1 σ √ τ -η ϕ (d 1 (S, τ -η; h(η))) ( 1 - K 2 h(η) r δ )} dη.
For a Euro-American strangle with an American lower side, we have

∆ EpAS EEP = ∫ τ 0 δe -δ(τ -η) {-N (-d 1 (S, τ -η; c 1 (η))) + 1 σ √ τ -η ϕ (-d 1 (S, τ -η; l(η))) ( 1 - K 1 l(η) r δ )} dη.
Gamma is the second derivative of the position with respect to the stock price. Viewed as the first derivative of the delta, this Greek parameter provides information on the rebalancing frequency of the self-financing replicating portfolio. It is a very important parameter for risk management because it informs option sellers about their exposure to an abrupt movement in the underlying asset price. Gamma is

Γ V = ∂ 2 V (S, τ ) ∂S 2 = Γ E + Γ EEP (13)
where Γ E is the gamma of a European strangle given by

Γ E = e -δτ Sσ √ τ (ϕ (d 1 (S, τ ; K 1 )) + ϕ (d 1 (S, τ ; K 2 ))) .
For an American strangle, we have

Γ AS EEP = ∫ τ 0 e -δ(τ -η) Sσ √ τ -η [ ϕ (d 1 (S, τ -η; c 2 (η))) { 1 - d 1 (S, τ -η; c 2 (η)) σ √ τ -η ( 1 - K 2 c 2 (η) r δ )} + ϕ (d 1 (S, τ -η; c 1 (η))) { 1 - d 1 (S, τ -η; c 1 (η)) σ √ τ -η ( 1 - K 1 c 1 (η) r δ )}]
dη.

For a Euro-American strangle with an American upper side, we have

Γ EcAS EEP = ∫ τ 0 e -δ(τ -η) Sσ √ τ -η ϕ (d 1 (S, τ -η; h(η))) { 1 - d 1 (S, τ -η; h(η)) σ √ τ -η ( 1 - K 2 h(η) r δ )} dη.
For a Euro-American strangle with an American lower side, we have

Γ EpAS EEP = ∫ τ 0 e -δ(τ -η) Sσ √ τ -η ϕ (d 1 (S, τ -η; l(η))) { 1 - d 1 (S, τ -η; l(η)) σ √ τ -η ( 1 - K 1 l(η) r δ )} dη.
The time decay theta is an important Greek parameter to consider (see Taleb (1997, pp.167-170) for a discussion). Because strangles satisfy the fundamental Black-Scholes Merton partial differential equation (Equation ( 1)), this Greek parameter can be deduced from the delta and gamma. We have

Θ V = ∂V ∂τ = -(r -δ) S∆ V - 1 2 σ 2 S 2 Γ V + rV. ( 14
)
Theta mixes delta and gamma in a quite subtle way (see Equation ( 14)). For American and Euro-American strangles, it is (of course) zero in stopping regions, because the time effect disappears once the contracts are exercised. Theta reaches its minimum when the option is approximately at the money (see Figure 5). Actually, the value of an option decreases when it nears expiry without being in the money because the probability that the option expires out of the money increases.

From a computational perspective, the hedging parameters can be easily evaluated once we know the early exercise boundaries. These boundaries are computed with the recursive method introduced earlier. Clearly, this method is quite likely to be more efficient than the alternative methods relaying on perturbation schemes. In perturbation analysis, repeating the price computations is necessary and this compounds numerical errors. By contrast, the direct computation of Greek parameters, as permitted by our formulas, restricts the errors to those associated with the determination of the boundaries (see [START_REF] Huang | Pricing and hedging American options: a recursive investigation method[END_REF] for an investigation of the recursive method).

"INSERT FIGURE 5 ABOUT HERE".

Figure 5 portrays the prices and Greek parameters of all the strangles on the same graph as a function of the underlying asset (contemporaneous) price. The upper-left graph illustrates the price differences between the strangles (for r < δ). 18 . The upper-right graph deals with the deltas. As expected, the deltas lie between -1 and 1. When the stock price is larger than K 2 , the call side of the American strangle is in the money and the right to buy becomes influential. When the stock price is lower than K 1 , the put side is in the money and the right to sell takes the lead. The gammas highlight the extinction of the early exercise strangles very clearly. The values simply drop to zero in such a case. All the gammas reach their maximum near strikes. The large gammas indicate that the delta hedging may be more demanding for these contracts compared to the European equivalent contracts.

Conclusions

This study reconsiders the pricing and hedging of strangles with finite maturity. It also introduces a new variant of these contracts called Euro-American strangles, and it provides a new (efficient and accurate) one-step numerical algorithm. The quadrature approach we propose combines Newton-Cotes and Gregory weighting schemes. Our algorithm avoids the non-monotonic gradient problem faced by the Chiarella and Ziogas' approach and has no recourse to any extrapolation scheme. We show that the computational time is significantly lowered and that the convergence rate is also improved. In addition, the new scheme allows us to control for errors. Our approach relies extensively on the Kim-style representation of American-style strangles. As a result of this, we can derive analytical expressions for prices and for hedging parameters. Finally, we make simulations to compare the parameters with each other. Our results suggest that the hedging of strangles deserves specific attention.

Appendix A. Solving PDE(1) with Kim's decomposition

The goal here is to solve the PDE (1) equipped with boundary conditions (3a)-(3d). Along the lines of [START_REF] Kim | The analytic valuation of American options[END_REF], we transform the problem so we can solve it using the results of [START_REF] Kolodner | Free boundary problem for the heat equation with applications to problems of change of phase. i. general method of solution[END_REF] and [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF]. Due to the linear property of differentiation, we can divide the Euro-American strangle with an American upper side into two components

EA h (S, τ ) = F (S, τ ) + G(S, τ )
where both F and G satisfy PDE (1). The function F is precisely characterized by

F (S, 0) = 0, S ≥ 0, (15)
F (S, τ ) = h(τ ) -K 2 , 0 ≤ τ ≤ T, ( 16
) lim S→h(τ ) ∂F ∂S = 1, 0 ≤ τ ≤ T. ( 17
) while G(S, 0) = max [0, K 1 -S] + max [0, (S -K 2 )ψ(h(0) -S)] , S ≥ 0, (18)
The solution of the second problem (i.e., the partial differential Equation ( 1) with boundary condition (18)) is equivalent to the price of a portfolio made up of two independent standard options (a call and a put), given that the price of the underlying asset at maturity is at most equal to h(0). In fact, in the case where there is no early exercise (of the American side), it is impossible for the underlying asset to be above h(0) (due to continuity of the h(.)). The linear property of differentiation allows us to divide the problem into two subproblems and solve them separately. The first challenge is to solve the PDE(1) under the boundary condition

G 1 (S, 0) = max [0, K 1 -S] , S ≥ 0.
The solution is available from the Black-Scholes framework and we have:

G 1 (S, τ ) = E [ e -rτ (K 1 -S τ )1 Sτ <K 1 |S 0 = S ] = ∫ ∞ -∞ e -rτ ( K 1 -Se (r-δ-σ 2 2 )τ +σ √ τ u ) 1 Se (r-δ-σ 2 2 )τ +σ √ τ u <K 1 1 √ 2π e -u 2 2 du = ∫ -ln(S/K 1 )+(r-δ-σ 2 2 )τ σ √ τ -∞ e -rτ ( K 1 -Se (r-δ-σ 2 2 )τ +σ √ τ u ) 1 √ 2π e -u 2 2 du
where E [.] is the mathematical expectation under the risk-neutral probability.

Separating the integral above yields:

G 1 (S, τ ) = K 1 e rτ N ( - ln(S/K 1 ) + (r -δ -σ 2 2 )τ σ √ τ ) -Se -δτ N ( - ln(S/K 1 ) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) .
The second challenge is to solve the PDE(1) under boundary condition

max [0, (S -K 2 )ψ(h(0) -S)] , S ≥ 0.
Similarly, its solution is given by:

G 2 (S, τ ) = E [ e rτ (S τ -K 2 )1 K 2 <Sτ <h(0) |S 0 = S ] = ∫ ∞ -∞ e rτ ( Se (r-δ-σ 2 2 )τ +σ √ τ u -K 2 ) 1 K 2 <Se (r-δ-σ 2 2 )τ +σ √ τ u <h(0) 1 √ 2π e -u 2 2 du = ∫ -ln(S/h(0))+(r-δ-σ 2 2 )τ σ √ τ - ln(S/K 2 )+(r-δ-σ 2 2 )τ σ √ τ e rτ ( Se (r-δ-σ 2 2 )τ +σ √ τ u -K 2 ) 1 √ 2π e -u 2 2 du = ∫ ln(S/K 2 )+(r-δ-σ 2 2 )τ σ √ τ ln(S/h(0))+(r-δ-σ 2 2 )τ σ √ τ e rτ ( Se (r-δ-σ 2 2 )τ -σ √ τ u -K 2 ) 1 √ 2π e -u 2 2 du = ∫ ln(S/K 2 )+(r-δ-σ 2 2 )τ σ √ τ -∞ e rτ ( Se (r-δ-σ 2 2 )τ -σ √ τ u -K 2 ) 1 √ 2π e -u 2 2 du - ∫ ln(S/h(0))+(r-δ-σ 2 2 )τ σ √ τ -∞ e rτ ( Se (r-δ-σ 2 2 )τ -σ √ τ u -K 2 ) 1 √ 2π e -u 2 2 du = Se -δτ N ( ln(S/K 2 ) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) -K 2 e -rτ N ( ln(S/K 2 ) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) -Se -δτ N ( ln(S/h(0)) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) + K 2 e -rτ N ( ln(S/h(0)) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) .
The overall solution of the first problem is therefore given by: [START_REF] Kolodner | Free boundary problem for the heat equation with applications to problems of change of phase. i. general method of solution[END_REF] and [START_REF] Mckean | Appendix: A free boundary problem for the heat equation arising from a problem in mathematical economics[END_REF] analyzed the heat equation that is a free boundary problem similar to the system-formed PDE (1) under conditions ( 15)-( 17). They solved this by means of Fourier transforms. Because it is long and tedious, we directly apply their results to our problem. We obtain:

G(S, τ ) = G 1 (S, τ ) + G 2 (S, τ ) = K 1 e -rτ N ( - ln(S/K 1 ) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) -Se -δτ N ( - ln(S/K 1 ) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) + Se -δτ N ( ln(S/K 2 ) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) -K 2 e -rτ N ( ln(S/K 2 ) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) -Se -δτ N ( ln(S/h(0)) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) + K 2 e -rτ N ( ln(S/h(0)) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) . ( 19
)
F (S, τ ) = ∫ τ 0 e -r(τ -ξ) e -[α-β(ξ)] 2 /2(τ -ξ) √ 2π(τ -ξ) × [ h(ξ)σ 2 + ( β ′ (ξ) - α -β(ξ) 2(τ -ξ) ) (h(ξ) -K 2 )
] dξ

where

β(ξ) = [ ln(h(ξ)) + ( r -δ - 1 2 σ 2 ) ξ ] /σ, α = [ ln(S) + ( r -δ - 1 2 σ 2 ) τ ] /σ.
By rearranging the terms of F (S, τ ) we have

F (S, τ ) = ∫ τ 0 e -r(τ -ξ) h(ξ) e -[α-β(ξ)] 2 /2(τ -ξ) √ 2π(τ -ξ) × [ σ 2 + β ′ (ξ) - α -β(ξ) 2(τ -ξ) ] dξ -K 2 ∫ τ 0 e -r(τ -ξ) e -[α-β(ξ)] 2 /2(τ -ξ) √ 2π(τ -ξ) × [ β ′ (ξ) - α -β(ξ) 2(τ -ξ) ] dξ. ( 20
)
And, by noticing that

e -r(τ -ξ) h(ξ) e -[α-β(ξ)] 2 /2(τ -ξ) √ 2π(τ -ξ) × [ σ 2 + β ′ (ξ) - α -β(ξ) 2(τ -ξ) ] = -e -δ(τ -ξ) S ∂ ∂ξ N ( α -β(ξ) + σ(τ -ξ) √ τ -ξ ) and e -r(τ -ξ) e -[α-β(ξ)] 2 /2(τ -ξ) √ 2π(τ -ξ) × [ β ′ (ξ) - α -β(ξ) 2(τ -ξ) ] = e -r(τ -ξ) ∂ ∂ξ N ( α -β(ξ) √ τ -ξ )
as well as using integration by parts, Equation(20) then gives

F (S, τ ) = δS ∫ τ 0 e -δ(τ -ξ) N ( α -β(ξ) + σ(τ -ξ) √ τ -ξ ) dξ -rK 2 ∫ τ 0 e -r(τ -ξ) N ( α -β(ξ) √ τ -ξ ) dξ + Se -δτ N ( ln(S/h(0)) + ( r -δ + 1 2 σ 2 ) τ σ √ τ ) -K 2 e -rτ N ( ln(S/h(0)) + ( r -δ -1 2 σ 2 ) τ σ √ τ ) . ( 21
)
Finally, by adding Equations ( 21) and ( 19) we obtain Equation (8). The proof of Proposition (3) is similarly obtained.

Appendix B. Derivation of Greek parameters.

In this appendix we derive the deltas of American strangles. The deltas for European and Euro-American strangles, as well as gammas, may be found similarly. Once we get deltas and gammas, thetas can be deduced by using the fundamental Black-Scholes PDE (Equation ( 14)). Let us consider an American strangle option

A c 1 c 2 (S, τ ) = E(S, τ ) + EEP c 1 c 2 (S, τ )
According to the linearity of differentiation, we have

∆ A = ∂A c 1 c 2 (S, τ ) ∂S = ∆ E + ∆ EEP
where ∆ E = ∂E(S,τ )

∂S

and ∆ EEP = ∂EEPc 1 c 2 (S,τ ) ∂S

. To simplify the proof let us divide ∆ E and

∆ EEP into two components. For ∆ E , we have ∆ E = ∆ E 1 + ∆ E 2 where ∆ E 1 = ∂ ∂S [ K 1 e -rτ N (-d 2 (S, τ ; K 1 )) -Se -δτ N (-d 1 (S, τ ; K 1 )) ]
and

∆ E 2 = ∂ ∂S [ Se -δτ N (d 1 (S, τ ; K 2 )) -K 2 e -rτ N (d 2 (S, τ ; K 2 ))
] .

For ∆ EEP , we have ∆ EEP = ∆ EEP 1 + ∆ EEP 2 where

∆ EEP 1 = ∂ ∂S [∫ τ 0 [ K 1 re -r(τ -η) N (-d 2 (S, τ -η; c 1 (η))) -Sδe -δ(τ -η) N (-d 1 (S, τ -η; c 1 (η))) ] dη ] and ∆ EEP 2 = ∂ ∂S [∫ τ 0 [ Sδe -δ(τ -η) N (d 1 (S, τ -η; c 2 (η))) -K 2 re -r(τ -η) N (d 2 (S, τ -η; c 2 (η))) ] dη ] .
For ∆ E , things are straightforward and well known. We have:

∆ E 2 = e -δτ N (d 1 (S, τ ; K 2 )) + Se -δτ 1 Sσ √ τ ϕ (d 1 (S, τ ; K 2 )) -K 2 e -rτ 1 Sσ √ τ ϕ (d 2 (S, τ ; K 2 ))
where the two last terms collapse. ∆ E 1 is obtained with a similar calculation and we have ∆

E 1 = -e -δτ N (d 1 (S, τ ; K 1 )). The next step is to investigate ∆ EEP 1 and ∆ EEP 2 . ∆ EEP 2 = ∫ τ 0 ∂ ∂S [ Sδe -δ(τ -η) N (d 1 (S, τ -η; c 2 (η))) -K 2 re -r(τ -η) N (d 2 (S, τ -η; c 2 (η)))
] dη

where the derivatives are

δe -δ(τ -η) N (d 1 (S, τ -η; c 2 (η)))+Sδe -δ(τ -η) ϕ (d 1 (S, τ -η; c 2 (η))) Sσ √ τ -η -K 2 re -r(τ -η) ϕ (d 2 (S, τ -η; c 2 (η))) Sσ √ τ -η
Given that Se -δτ ϕ (d 1 (S, τ ; K)) = Ke -rτ ϕ (d 2 (S, τ ; K)) is verified for all K, we obtain after simplification

δe -δ(τ -η) N (d 1 (S, τ -η; c 2 (η))) + δe -δ(τ -η) ϕ (d 1 (S, τ -η; c 2 (η))) σ √ τ -η ( 1 - K 2 c 2 (η) r δ
) With all in one we get

∆ EEP 2 = ∫ τ 0 δe -δ(τ -η) [ N (d 1 (S, τ -η; c 2 (η))) + ϕ (d 1 (S, τ -η; c 2 (η))) σ √ τ -η ( 1 - K 2 c 2 (η) r δ )]
dη.

We can derive the expression for ∆ EEP 1 with the same procedure. We get

∆ EEP 1 = ∫ τ 0 δe -δ(τ -η) [ -N (-d 1 (S, τ -η; c 1 (η))) + ϕ (-d 1 (S, τ -η; c 1 (η))) σ √ τ -η ( 1 - K 1 c 1 (η) r δ )]
dη.

Algorithm 2: Euro-American strangle with early exercisable call side . Function: EAS Input: S, K 1 , K 2 , r, q, sigma, i, h, c 2 /* t = ih stands for time to maturity */ Output: A /* price of the strangle */ begin t ←-ih; E ←- The graphs of these figures represent the prices and Greeks of the following strangles: European (gray line), American (black line), Euro-American with American call side (marker "+"), and Euro-American with American put side (marker "o"). The parameter values are K 1 = 1, K 2 = 1.001, r = 10% (r = 0.5 for prices graph), δ = 10%, T = 1, σ = 20%. "Total Subint." represents the total number of subintervals needed for implementing the approach. "time" is the computational time needed to compute American strangle prices. The parameter values are K 1 = 1, K 2 = 1.5, r = 5%, δ = 10%, σ = 20%, and τ = 1. Notice that for n = 100, (100 + 2 × 100) = 300 subintervals have been used by CZ whereas LAM requires 2 × 100 subintervals only. The parameter values are K 1 = 1, K 2 = 1.001, T = 1, and σ = 20%. Both r, δ, and S vary.

( K 1 e -rt N (-d 2 (S, t; K 1 )) -Se -qt N (-d 1 (S, t; K 1 )) + Se -qt N (d 1 (S, t; K 2 )) -K 2 e -rt N (d 2 (S, t; K 2 )) ) ; EEP (v) ←- ( Sqe -q(i-v)h N (d 1 (S, (i -v)h; c 2 (vh))) -K 2 re -r(i-v)h N (d 2 (S, (i -v)h; c 2 (vh))) ) ; /*
; c 2 (0) ←max(K 2 , K 2 × (r/q)) ; /* recursive computation of c 2 starting with c 2 (0) */ for i ← 1 to n do f 2 (x) ← x -K 2 -EAS(x, K 1 , K 2 , r, q, sigma, i, h, c 2 ) /* new function */ c 2 (i) ←Bisection(f 2 (x)); /* Bisection(f ) gives the root of f */ return c2

Figure 1 .

 1 Figure 1. American strangles continuation & stopping regions.

Figure 2 .

 2 Figure 2. Comparing numerical strategies.

Figure 3 .

 3 Figure 3. Call side early exercise boundaries: CZ's algorithm.

Figure 4 .

 4 Figure 4. Call side early exercise boundaries: a comparison.

Figure 5

 5 Figure 5. Simulations.

Table 2 .

 2 Hybrid quadrature weights.

	n	w 0	w 1	w 2	w 3	w 4	w 5	w 6	w 7	.	w i

Table 4 .

 4 Values of American strangles: convergence and speed.

	Method	CZ

Table 5 .

 5 Prices of American strangles vs. Euro-American strangles.

	S	American	Euro-American	Euro-American	European strangle
		strangle	strangle with	strangle with
			American call side	American put side
			Panel A: r < δ (r = 5%, δ = 10%)
	1.0	0.156485	0.156485	0.152081	0.152081
	1.1	0.169608	0.169608	0.158479	0.158479
	1.2	0.219351	0.219351	0.195680	0.195680
	1.3	0.299809	0.299809	0.255461	0.255461
			Panel B: r > δ (r = 10%, δ = 5%)
	0.7	0.299959	0.247329	0.299959	0.247329
	0.8	0.205527	0.179077	0.205527	0.179077
	0.9	0.156331	0.145270	0.156331	0.145270
	1.0	0.156364	0.151920	0.156364	0.151920
			Panel C: r = δ (r = δ = 5%)
	0.9	0.165262	0.163381	0.165033	0.163152
	1.0	0.152594	0.151846	0.151853	0.151104
	1.1	0.178339	0.178052	0.176449	0.176162
	1.2	0.234506	0.234400	0.230451	0.230345

The self-closing mechanism has a number of interesting features for managers. Unlike sellers of portfolios made up of individual American options, sellers of American strangles avoid the risk of successive exercises because American strangles disappear once the early exercise is decided upon. For buyers, these contracts are cheaper than their rough synthetic portfolios.

Studying contracts with no maturity,[START_REF] Moraux | On perpetual American strangles[END_REF] found that holders of individual American options can exercise earlier than holders of American strangles, can receive a greater value at the exercise time, and if they manage their position through the synthetic portfolio characteristics may exercise their right to sell or to buy rather suboptimally.

CZ wrote, "the numerical scheme is firstly carried out using a time-step size of hand is then repeated using h/2. In each case, since it is necessary to alternate between two different numerical integration schemes (for odd and even values), it turns out that the free boundaries have non-monotonic gradients. This is rectified by combining the two estimates using Richardson's extrapolation. Pricing the American strangle is then achieved via numerical integration using Simpson's rule, combined with the estimates".

Perpetual (American) strangles have also been studied in the literature.[START_REF] Gerber | Martingale approach to pricing perpetual American options[END_REF] provided the first pricing formula, but they did not solve the system of equations locating early exercise boundaries, nor did they provide any simulations.[START_REF] Moraux | On perpetual American strangles[END_REF] analysed these contracts as asymmetric rebates of double knock-out barrier options with special payoffs and found early exercise boundaries numerically. He then provided simulations and insights on price properties and on the optimal exercise policy as well as discussing hedging issues.

E.g., CZ report a moderate 6% for the largest bias in their simulations.

By way of comparison, the existing two-step quadrature approach over/under-estimates the true price unexpectedly.

see[START_REF] Macmillan | An analytical approximation for the American put prices[END_REF] and[START_REF] Barone-Adesi | Efficient analytic approximation of American option values[END_REF] for other approaches using the early exercise premium.

Standard early exercise boundaries are not suitable anymore due to the self-closing mechanism.

We know that holding an American-style strangle is not equivalent to being long in a portfolio made up of an American call option (K 2 , T ) and an American put option (K 1 , T ). Both rights (to sell/to buy) disappear when the owner decides to exercise one.

This can be compared with CZ, who essentially favor another decomposition; that is, A c1c2 (S, τ ) = C(S, τ, c 2 (.)) + P(S, τ, c 1 (.)), where C(S, τ, c 2 (.)), and P(S, τ, c 1 (.)) stand for the prices of standard American calls and puts that use c 2 (.) and c 1 (.) for the respective early exercise boundaries.

Note that the Newton-Cotes weights may become negative for i greater than six and this causes numerical instability.

For more details on the error analysis, readers may consult[START_REF] Quarteroni | Numerical mathematics[END_REF].

For more efficiency,Linz (1985, p.98) recommends the use of the high-order quadrature when computing a Volterra integral equation type with lower values of i.

Following[START_REF] Kim | The analytic valuation of American options[END_REF], CZ give the proof for c 1 (0) and c 2 (0). We can obtain the results of h(0) and l(0) similarly.

In this case, the Euro-American strangle price with an American call side is similar to that of an American strangle.

* Both authors acknowledge the financial support received from IAE Rennes and CREM (CNRS benefited from the comments of Carol Alexander (the Co-editorin-chief), two anonymous referees, some participants at the 2013 AFFI International Conference (Van Son Lai and the discussion participant Victor Reutenauer), the 8th Bachelier Society Congress (special thanks to

Appendix C. Pseudo-codes

Algorithm 1: Hybrid quadrature Function: HybridQuad Input: n /* number of subinterval */ Output: w /* vector of weights of n + 1 endpoints */ begin Create vector w of length n + 1 and full with default value 1; switch n do case n = 1 /* Closed Newton-Cotes weights */ w(0) ←-0.5; w(1) ←-0.5;

-

With k ∈ {1, 2, 3, ...}. For odd values of n, the quadrature does not exist theoretically. Consequently, inserting an arbitrarily additional coefficient (4/3 or 2/3) for odd values of n, as CZ did, may lead to an additional bias in the integral approximation.

Table 3. Degree of Exactness (DE) and error associated with Newton-Cotes.

i DE


Daniel Marian Merezeanu
email: daniel@aii.pub.ro

Daniela Andone

Tandem Queueing Systems Maximum Throughput Problem

Keywords: tandem queues, production line, maximum throughput, QBD process

In this paper we consider the problem of maximum throughput for tandem queueing system. We modeled this system as a Quasi-Birth-Death process. In order to do this we named level the number of customers waiting in the first buffer (including the customer in service) and we called phase the state of the remining servers. Using this model we studied the problem of maximum throughput of the system: the maximum arrival rate that a given system could support before becoming saturated, or unstable. We considered different particular cases of such systems, which were obtained by modifying the capacity of the intermediary buffers, the arrival rate and the service rates. The results of the simulations are presented in our paper and can be summed up in the following conclusions: 1. The analytic formula for the maximum throughput of the system tends to become rather complicated when the number of servers increase 2. The maximum throughput of the system converges as the number of servers increases 3. The homogeneous case reveals an interesting characteristic: if we reverse the order of the servers, maximum thruoughput of the system remains unchanged The QBD process used for the case of Poisson arrivals can be extended to model more general arrival processes.

INTRODUCTION

In this paper we analyze the problem of maximum throughput of a tandem queueing system. We consider a system with Poisson arrivals and independent exponential service times. This kind of system is appropriate for representing the production lines -the basis of flexible manufacturing systems. Basically a production line consists of a one-way serial configuration of machines and buffers. To model this as a tandem queues system we will assimilate the machines to servers and the buffers to queues placed between these servers. The tandem queues system obtained is presented in figure 1, This system consists of K+1 servers (notated i S , where i=0,1,…,K) and the same number of buffers (the buffers capacity is notated i B , where i=0,1,…,K). In order to analyze the problem of maximum throughput for this system we made the following assumptions:

The first buffer ( 0 The goal of this paperis to study the maximum arrival rate ( max λ ) that the system can support before becoming saturated or unstable, the so-called maximum throughput. In order to do this we model as a quasi-birth-and-dead (QBD) process.We simulated the behavior of the QBD process in order to obtain the formula and values for the maximum throughput for different configurations. We present the results of these simulations in the paper.

GENERATING THE QBD MODEL

In order to model the given tandem queueing system as a QBD process we proceeded as follows. We named level the number of customers in the buffer before server 0 S (including the customer in service) and we called phase the state of the remaining servers. Considering these notations, the only events that can change the state of the system are the arrival of the new customer or a service completion by one of the servers. Since the probability of two of these events occurring simultaneously is zero and the level can icrease or decrease by at most one on the occurrence of a single event, this means we have infinite, continous-time QBD process.

To specify the phase we need only to consider the case in which server 0 S is occupied, so we encode the phase as an ordered K-tuple (

+ = + i i B m
is not a valid one, since it is not possible for the server i S to be both empty and blocked simultaneously. This allows us to determine the number of valid phases M as:

( ) ( )( ) [ ] ( ) ( )( ) [ ] ( )( ) 5 1 2 5 1 3 5 1 3 1 1 1 + + + + + + - + + + + + + + + B B B B B B B B M K K K
(1) where B represents the capacity of the intermediary buffers. We used this model to study the maximum tgroughput of the considered system.

MAXIMUM THROUGHPUT PROBLEM

In order to study the maximum throughput problem, we should generate the infinitesimal generator matrix Q of the defined QBD process. The Q matrix is a tridiagonal block matrix having the structure presented below:

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = L L L L L L L L L 1 2 0 1 2 0 1 2 0 1 0 0 0 0 0 0 A A A A A A A A A B Q (2)
where block 1 A accounts for those transitions where the level remains unchanged, 0 A accounts for those transitions where the level increases, and 2 A accounts for those transitions where the level desreases. Block 1 B corresponds to the transitions from level 1 to level 0 (there are no customers in the input buffer 0 B ) and it is of no interest for our purposes.

Having the Q matrix, the next step is to determine tge condition of stability of the system. Let where e is a column vector of ones. Then the system is stable if and only if:

e A e A 0 2 π π > (4)
In computing i A , i=0,1,2 it is enough to consider the case where the level is greater than 0. Then any arrival will cause the level to increase, but the phase will remain the same. In consequence

I A λ = 0 .
Therefore, when we compute A we will always have an arrival added to the diagonal of A from 0 A but by the same token we will always have an arrival subtracted from the diagonal of A because of the λ on the diagonal of 1

A that we need to make the sum of the rows of Q equal to 0. That is, in computing the A matrix, arrivals make no difference and can be ignored.

To get the maximum value of λ (maximum throughput) we need to solve: . In order to obtain the maximum throughput ( max λ ) for the steady state of the QBD process we only need to solve the resulting system. The only remaining problem is to determine the structure of the blocks for matrix Q. The 0 A block (corresponding to the transitions to the next level) is easy to obtain. As we already saw, the level can only increase with the arrival of a new customer into the system. This doesn't change the phase, so 0 A will be an Mdiagonal matrix of λ . Yet the process of generating the other two blocks 1 A and 2 A is not so straightforward. In order to generate 1

A and 2 A we should analyze each phase and decide which are the possible transitions and determine the resulting phase and level. We made the observation that the level decreases only at the moment when server 0 S completes the service of a customer and may pass it to server 1 S . There are two possible situations: if server 0 S was no the server that completed the service, then the only way th elevel can decrease is if the server 1 S is blocking 1 S . The complexity of the model made this analysis possible only by using numerical methods, so we develop an algorithm to do that. We use this algorithm to simulate different configurations of the system. The following sections describe the simulation algorithm and the results we obtained.

SIMULATION ALGORITHM

The algorithm's steps are: 1. Initialize the process; read the number of servers into the tandem queueing system, their rates and the capacities of the intermediary buffers. 2. Determine the number of valid phases and which are these valid phases. 3. Analyze each valid phase in order to find the events that can occur next and generate the blocks of the Q matrix, according to the following rules:

-for the arrival of a customer, the level increase and the phase remain unchanged; in consequence the block A are determined by the equalizing with zero the sum of the rows of the Q matrix. 4. Determine the maximum throughput of the system.

We will illustrate this process for a particular configuration, consisting of two servers (K-1) and an intermediary buffer with capacity B=2. This particular structure is presented in figure 2. The associate QBD model for the above tandem queueing is presented in figure 3.

Figure 1 -

 1 Figure 1 -Tandem queueing system's structure

Figure 2 -

 2 Figure 2 -Example of a tandem queueing system with two servers and a finite intermediary buffer

Figure

  Figure 3 -The associated QBD model of the tandem queueing system

and the formula of the maximum throughput obtained for this particular configuration is:

( )

EXPERIMENTAL RESULTS

We use this simulation algorithm for different configurations of the tandem queueing system having a different number of servers, with different service rates, with or without intermediary buffers. Table 1 summarizes some of the results of the simulation for configuration without intermediary buffers. We change only the service rate of the first server 0 In figure 4 we present a graphic representation of the maximum throughput of the system. The results for the same configurations of the tandem queueing system, using intermediary buffers with capacity B=1 are presented in table 2. The graphic representation of the evolution of maximum throughput for a tandem queueing system with intermediary buffers is shown in figure 5.

Figure 5 -Maximum throughput for the system with intermediary buffers

CONCLUSIONS

The results obtained allow us to formulate some interesting conclusions. First we should notice that the complexity of the system significantly increases when the number of servers and the capacity of the analytic formula of the system's maximum throughput tends to become rather complicated.

We also noticed that the evolution of the maximum throughput of a tandem queueing system converges as the number of servers increases. Also the results prove that the intermediary buffers increase the value of the maximum throughput of the system. Another observation is that the homogeneous structure (all the intermediary buffers having the same capacity) reveals an interesting characteristic, that if the order of the servers is reversed, the maximum throughput of the system remains unchanged. In further research our efforts will focus on extending the application of the QBD process for more general arrival processes. Our immediate plans are to adapt the method presented in the paper in studying the maximum throughput problem for a tandem queueing system with n MMPP arrivals and to find a way to improve the simulation algorithm and to reduce its complexity.


Andrei Hossu
email: hossu@aii.pub.ro

Daniela Andone

A NEW APPROACH OF GRAY IMAGES BINARIZATION WITH THRESHOLD METHODS

Keywords: Vision systems, Gray level image binarization, gray level histogram, global optimum thresholding, dynamic optimum threshold, temporal histogram, temporal thresholding and moving scene in robotic automation

The paper presents some aspects of the (gray level) image binarization methods used in artificial vision systems. It is introduced a new approach of gray level image binarization for artificial vision systems dedicated to industrial automationtemporal thresholding. In the first part of the paper are extracted some limitations of using the global optimum thresholding in gray level image binarization. In the second part of this paper are presented some aspects of the dynamic optimum thresholding method for gray level image binarization. Starting from classic methods of global and dynamic optimal thresholding of the gray level images in the next section are introduced the concepts of temporal histogram and temporal thresholding. In the final section are presented some practical aspects of the temporal thresholding method in artificial vision applications form the moving scene in robotic automation class; pointing out the influence of the acquisition frequency on the methods results.

GLOBAL THRESHOLD FOR IMAGE BINARIZATION

Threshold methods are defined as starting from the analyze of the values of a function T of the type: , y, p(x, y), f(x, y)]

T = T [x
(1)

Where: f(x, y) -represents the intensity value of the image element located on the coordinates (x, y); p(x,y) -represents the local properties of the specific point (like the average intensity of a region centered in the coordinates (x, y)). T -is the binarization threshold The goal is to obtain from an original gray level image, a binary image g(x, y) defined by:

⎩ ⎨ ⎧ ≤ > = T y x f T y x f y x g ) , ( 0
) , ( 1 ) , ( for for (2)
For T a function only of f(x, y), the obtained threshold is called global threshold.

In the case of T a function of both f(x, y) and p(x, y), the obtained threshold is named local threshold.

In the case of T a function of all f(x, y), p(x, y), x and y, the threshold is a dynamic threshold.

Intensity level normal distribution assumption

Gray level histogram represents the probability density function of the intensity values of the image.

In order to simplify the explanations, we suppose the image histogram of the gray levels is composed from two values combined with additive Gaussian noise:

-The first segment of the image histogram corresponds to the background points -the intensity levels are closer to the lower limit of the range (the background is dark) -The second segment of the image histogram corresponds to the object points -the intensity levels are closer to the upper limit of the intensity range (the objects are bright).

The problem is to estimate a value of the threshold T for which the image elements with an intensity value lower than T will contain background points and the pixels with the intensity value greater than T will contain object points, with a minimum error.

For a real image, the partitioning between the two brightness levels (the background and the object) is not so simple and also not so accurate.

The partitioning is fully accurate only if the two modes of the bimodal histogram are not overlapped. In Fig. 1 is presented a possible shape of a bimodal histogram with overlapped modes.

Fig.1 Bimodal histogram with overlapped modes

The classification is defined as the process of the distribution of the pixels in classes (the class of the object pixels and the class of background pixels).

The goal of the binarization process is the minimization of the error of classification. The optimum binarization threshold is located in the intersection position of the two normal distributions.

The estimation of the error of classification is obtained from the area of the overlapped segments:

image size B A E + = (3)
Suppose the image contains two intensity level values affected with additive Gaussian noise. The mixture probability density function is:

) ( ) ( ) ( 2 2 1 1 x p P x p P x p + = (4)
Where:

x -the random value representing the intensity level, p 1 (x), p 2 (x) -are the probability density functions, P 1 , P 2 -are the a priori probabilities of the two intensity levels (P 1 + P 2 = 1). For the normal distribution case on the two brightness levels:

2 2 2 2 2 2 2 1 2 1 1 1 2 ) ( exp 2 2 ) ( exp 2 ) ( σ µ - πσ + σ µ - πσ = x P x P x p (5)
Where: -are the standard deviations of the two statistical populations. Suppose the background is darker than the object. In this case 2 1 µ < µ and defining a threshold T, so that all pixels with intensity level below T are considered belonging to the background and all pixels with level above T are considered object points. The probability of misclassification an object point (classifying an object point as a background point) is: Similarly, E 2 :

dx x p T E dx x p T E T T ) ( ) ( ) ( ) ( 1 2 2 1 ∫ = ∫ = ∞ + ∞ -
The probability of error is given by:

) ( ) ( ) ( 1 2 2 1 T E P T E P T E + = (7)
To find the threshold value for which the error is minimum, E(T) is differentiate with respect to T:

) ( ) ( 2 2 1 1 t p P t p P = (8)
Applying the result to the Gaussian density we obtain:

AT 2 + BT + C = 0, (9)
Where:

2 2 1 1 2 2 2 1 2 1 2 2 2 2 2 1 2 1 2 2 2 1 2 2 2 1 ln ) ( 2 P P C B A σ σ σ σ + µ σ - µ σ = σ µ - σ µ = σ - σ =
If the standard deviations are equal, a single threshold is sufficient:

1 2 2 1 2 2 1 ln 2 P P T µ - µ σ + µ + µ = (11)
If the probabilities are equal P 1 = P 2, the threshold value is equal with the average of the means. A way of checking the validity of the assumption of bimodal histogram is to estimate the mean-square error between the mixture density, p(x) and the experimental histogram h(x i ).

2 1 )] ( ) ( [ 1 i N i i x h x p N M - ∑ = = (12)
Where: N -number of possible levels of the image (usually N = 256)

The image binarization is obtained changing the color attribute of each pixel according to its intensity level relative to the binarization threshold.

(6) (10) Characteristics of the global thresholding methods: -The assumption that both classes have the same standard deviation is acceptable, but the assumption the classes (two levels) have the same a priori probabilities in many applications is not acceptable. In the case of the artificial vision systems dedicated to object recognition for industrial applications there is a large amount of a priori information about the image that has to be processed. Can be obtained better results of estimation of the distribution (standard deviation and the mean value) of the image elements of the scene (background image -without the objects). Usually, in robotic applications, the illumination environment is known and controlled and also the object classes with a probability of apparition in the image are not known. In many robotic application an estimation of the ratio between the area of the objects to be analyzed and the total area of the image scene, can be made with good results (a batter estimation than the assumption of P 1 = P 2 = 0.5).

DYNAMIC THRESHOLD FOR IMAGE BINARIZATION

There are some classes of scenes of artificial vision systems where using the global threshold methods is not acceptable:

• The case of the applications where the lighting system does not supply a uniform intensity all over the analyzed surface. • Segments of the image (or some times, image elements) do not have the same behavior in the same lighting conditions.

For these types of images, for binarization of the image, the most often used are dynamic threshold methods.

The algorithm of the estimation of the dynamic threshold consist of:

The methods are based on the local analyze of the image.

• The original image is divided in regions of a prescribed size. -Large size of the region makes the method to loose the dynamic threshold characteristics and to fail into a global threshold method -Small size of the region makes to loose the statistical characteristic of the population of the image elements contained by the analyzed region (and the accuracy of the results is lost). -The last comment on the method is the fact that this method does not solve the problem of the non-uniformity of the illumination system or of the acquisition sensor.

TEMPORAL HISTOGRAM

For the class of artificial vision systems dedicated for moving scene (used very often in inspection and robotic applications) three types of image intensity level distortions can be identified:

• Illumination nonuniformity (obtaining a uniform intensity of the light on the whole area of the scene where the image is analyzed -usually 2 m -it is practical impossible). The main problem of the methods presented before represents the assumption that the image is a statistical population obtained from the addition of two ore more distributions (in the general accepted case normal distribution).

In the general case (an array image) an image represents a data set of:

{f(x, y) | x ∈ [0,N], y ∈ [0,M]}, (13)
where:

N represents the number of image elements per row (number of image columns), M represents the number of image elements per column (number of image rows),

In the linear image case, this data set become:

{f(x) | x ∈ [0,N]}, (14)
where:

N represents the number of image elements per row (number of image columns),

This assumption on the distribution of the intensity levels has the starting point the assumption that the insertion point of the noise is located on the transmission level of the information. In other words, the assumptions is that:

-The acquired image is an ideal image (with only two gray levels: the gray level of the scene pixels and the gray level of the pixels corresponding to the object) -Then a global noise is applied, transforming the two levels in two normal distributions.

The assumption is false and using it we are analyzing a histogram, which is far away of two normal distributions, and from here the results are distorted.

In reality the noise on the intensity level has its insertion point on the acquisition level (noise on the "intensity source", sampling, holding and conversion, etc) and not on image transmission level. Intensity source has the meaning of intensity signal on the acquisition element and not only the lighting system. This implies the fact that the noise on the intensity source represents the whole chain of: lighting source noise, reflective characteristics of the object surface and reflective characteristics of the scene surface and the sensitivity characteristics of the sensor.

Moving the insertion point of the noise we obtain:

In the general case (an array image) an image represents a data set of:

{f i (x, y) | i ∈ [0,L]}, (15)
where:

L represents the number of the image frames (the size of the statistic population analyzed),

x ∈ [0,N], N representing the number of image elements per row (number of image columns), y ∈ [0,M], M representing the number of image elements per column (number of image rows).

In the linear image case, this data set become:

{f i (x) | i ∈ [0,L]}, ( 16
)
where:

L represents the number of the image frames (the size of the statistic population analyzed),

x ∈ [0,N], N representing the number of image elements per row (number of image columns).

In this way several temporal built statistical populations (from intensity levels of the same image element on a set of image frames acquired on different moments) replace the spatial built statistical population (made from image elements of the same image).

The method of temporal histogram has the result the fact that each element of this set of histograms represents a bimodal histogram with two not overlapped modes (in case of a correct acquisition environment). It can be also introduce an estimation of the quality of the acquisition and binarization process using the estimation of the misclassification error analyzing the parameters of the two normal distributions (one corresponding to the object and the other one corresponding to the scene).

The method offers also the capacity of identification of the areas where some modifications should be done (on the lighting system) in order to improve the quality of the acquisition and binarization process.

The lack of the proposed method is the memory consumption (it has to be built N x M different histograms in array acquisition, or N -in linear acquisition case). This problem is not so restrictive because at the end only the threshold values have to be stored and not the whole histograms.

Another restriction is the fact that the method requires a large number of image frames acquired for construction of the statistical populations (in application setup time). In the case of the systems dedicated to industrial applications usually this does not represent a real problem. This type of applications does not require a system response in condition of a small number of image frames a priori acquired (but other type of applications has these requirements).

The vision systems dedicated to industrial applications can take the advantage on the fact that the image environment does not change a lot in time.

In this way it can be initially reserved a certain time (before the runtime of the system) for acquiring a large enough number of image frames in order to be able to identify the permanent characteristics of the environment. All the intensity level distortions (described in this section) presents permanent characteristics.

Using this method is a necessity for the artificial vision systems dedicated to applications where the errors on binarization are not acceptable (inspection and measurement systems -error free systems). In the applications dedicated exclusively to shape recognition the errors are accepted in a predefined range.

THE INFLUENCE OF THE ACQUISITION FREQUENCY ON THE BINARIZATION THRESHOLD VALUE

In moving scene applications, in order to maintain a constant resolution of the vision system along the direction of the scene movement, it is necessary the ratio between the acquisition frequency (the image lines rate -in the case of a line scan camera) and the scene speed to be constant.

The acquisition frequency determines the exposure time of the CCD sensor cells. It can be notice an important influence of the speed (of the conveyor) on the intensity level of the same image element in the same lighting environment.

In Fig. 3 are presented the experimental results obtained analyzing the influence on the intensity levels (for both: bright object and dark background) of the speed of the conveyor (acquisition frequency). The results were obtained on a statistical population from an image element on each measured speed. The second column represents the measured speed of the scene (conveyor) -V [m/min]. The 3 rd to 8 th columns represents image intensity levels (clarified in Fig. ) estimated from the analyzed statistical population (temporal histogram). The values from the Threshold column are the binarization threshold values obtained from a global optimum temporal thresholding method applied on the histogram built for each analyzed level of the speed.

In Fig. 4 are presented graphical the explanations on the meanings of the data involved in the analysis of the influence of the speed (acquisition frequency) on the intensity levels. The artificial vision system benefits from these results using a relation between the value of the binarization threshold and the speed V of the scene (acquisition frequency)

T = T (x, V) (17)
Because of the response time restrictions imposed to the artificial vision system, instead of using an explicit expression of the estimated function T(x, V), a search method in an a priori filled Fig. 3 The influence on the intensity levels of the speed of the scene (acquisition frequency) Fig. 4 The influence of the speed on the intensity levels

CONCLUSIONS

For the class of artificial vision systems dedicated for inspection and measurement industrial applications the error on binarization process is not acceptable. In this, case classic methods like global, local and dynamic threshold are not applicable. The paper introduces a new approach of gray level image binarization -temporal thresholding.

For the class of artificial vision systems dedicated for moving scene the acquisition frequency is dependent on the speed of the transmission support (usually a conveyor). To solve this problem, the artificial vision system has to estimate the influence of the acquisition frequency on the histogram and on the binarization threshold values.

The paper proposes a processing time efficient method to estimate the binarization threshold for the case of an error free vision system in the case of variation of the acquisition frequency.

  In the final stage, a second interpolation process is applied: for each image element is assigned a threshold value T(x, y) from the interpolation of the values of the neighbor image elements. The method is called dynamic thresholding because the value of the resulted threshold for each image element is dependent of the position of the element in the image -T(x, y).

	threshold values (for which the assumption
	of a bimodal histogram is considered
	acceptable). • Characteristics of the dynamic thresholding methods:
	-Lack of processing time consumption -each
	element of the image is used at least two times
	(the method requires multiple-pass of the image)
	in different steps of the algorithm (and the
	number of the elements is very large).
	-Estimation of the acceptable error value (or the
	validation of the bimodal histogram assumption)
	is a complex process.
	-To choose the size of the image regions we have
	to take into account:
	• For each region it is estimated the histogram • For each histogram it is estimated the error
	induced from the assumption of bimodal
	histogram (a histogram built from two normal
	distributions) • If the value of the error is less than an
	acceptable value, the global threshold for
	the region is estimated. • If the value of the error is too big (the
	histogram is too far from a bimodal
	histogram) the threshold value for
	binarization is estimated from the
	interpolation of the neighbors region

  table (at setup time) is more appropriate.

		250			Intensity levels relative to the scene speed
		200
									Object Min+
									Object Max
	Intensity levels	100 150							Object Min-Scene Min+ Scene Max Scene Min-Threshold
		50
		0
		10	20	30	40	50	60	70	80
					Speed V [m/min]
							The size of the table is 256 (the number of the
							possible values of the binarization thresholds),
							containing floating-point values of the speed of the
							conveyor (acquisition frequency) for which the value
							of the binarization threshold has to be changed.

  Croicu, A., Hossu, A., Dothan, E., Ellenbogen, D., Livne, Y. (1998)-ISCAN-Virtual Class based Architecture for Float Glass Lines, IsoCE'98, Sinaia, June Hossu, A., Croicu, A., Dothan, E., Ellenbogen, D., Livne, Y.(1998) -ISCAN Cold-Side Glass Inspection System for Continuous Float Lines, User Manual, Rosh-Haayn,.

	Number of	Scene Max	Object Max
	Scene Min-pixels	Object Min-Scene Min+	Object Min+
		Threshold
	Number of pixels			Intensity level
				6. REFERENCES
				Increased speed
				influence on intensity levels	Haralick, R., Shapiro, L. (1992) Computer and Robot Vision, Addison-Wesley Publishing
				Intensity level	Company, Borangiu, Th, Hossu A., Croicu, A. (1994) -
				ROBOTVISIONPro, Users Manual, ESHED
				ROBOTEC, Tel -Aviv,.