A penalty-based method from reconstructing smooth local volatility surface from American options

This paper is devoted to develop a robust penalty-based method of reconstructing 
smooth local volatility surface from the observed American option 
prices. This reconstruction problem is posed as an inverse problem: 
given a finite set of observed American option prices, find a local 
volatility function such that the theoretical option prices matches 
the observed ones optimally with respect to a prescribed performance 
criterion. The theoretical American option prices are governed by 
a set of partial differential complementarity problems (PDCP). We 
propose a penalty-based numerical method for the solution of the PDCP. 
Typically, the reconstruction problem is ill-posed and a bicubic spline 
regularization technique is thus proposed to overcome this difficulty. 
We apply a gradient-based optimization algorithm to solve this nonlinear 
optimization problem, where the Jacobian of the cost function is computed 
via finite difference approximation. Two numerical experiments: a 
synthetic American put option example and a real market American put 
option example, are performed to show the robustness and effectiveness 
of the proposed method to reconstructing the unknown volatility surface.


1.
Introduction. In the standard Black-Scholes model of option pricing [5], the volatility is assumed to be a constant. It is also well known that with this strong assumption the Black-Scholes model cannot capture features of empirically observed option prices. In fact, the so-called volatility skew or smile phenomenon [13,10] exists in all the major stock index markets today. To capture the existence of volatility smile and skew, several different extensions of the Black-Scholes model have been proposed, such as the stochastic volatility approach [10,13], the jump diffusion model [2,17], the uncertain volatility model [3,24] and the deterministic local volatility function approach [1,9]. As the stochastic volatility model, jump diffusion model and uncertain volatility model introduce some additional non-traded sources of risk besides the risk of underlying assets, the completeness of the market is no longer maintained. On the other hand, the deterministic local volatility approach does not introduce any other source of risk, and hence it can maintain the market completeness. In a complete market, arbitrage pricing and hedging are allowed. Thus the deterministic local volatility model is very popular and attractive in financial practice.
Despite the attractiveness of local volatility model, the finding of the local volatility function is a very challenging task. One popular method is to reconstruct the volatility surface from the observed market option prices. However, since the market option prices are typically limited to a relatively few different maturities and strikes, the method is usually regarded as a nonlinear least square problem [7,1,9]. In the case of European options, this nonlinear least square problem is stated as: given a finite set of observed European option prices, find a local volatility function such that the theoretical option prices matches the observed ones optimally with respect to a prescribed performance criterion. As the theoretical European option prices are governed by a set of Black-Scholes partial differential equations (PDEs), this nonlinear least square problem actually is a distributed optimal parameter selection problem. This problem has been extensively studied both in the case of continuous space and finite-dimension space, using optimal control approach ( [1,14,12]) and mathematical program approach ( [9]) respectively. The idea of solving distributed optimal parameter selection problem is usually accomplished by developing a gradient-type optimization algorithm, coupled with some numerical implementation of solving the governing PDEs, descriing ideal option prices.
The nonlinear least square approach to the reconstruction of the volatility surface is essentially an inverse problem of a forward option pricing. Typically, it is ill-posed. To overcome the difficulty, some regularization techniques are commonly used. One technique is Tikhonov regularization, see [1,8,14,6]. Unfortunately, this approach requires the solution of a very large-scale nonlinear optimization problem, where the dimension is equal to the total number of discretization points. In addition, it requires the determination of a regularization parameter. Another regularization technique is to approximate the local volatility surface directly by using some smoothing interpolation methods, such as the cubic spline smoothing [7,21]. In [11,7], this smoothing technique has been proposed to the reconstruction of the volatility surface from the European option, where its robustness and effectiveness are being verified by several numerical experiments.
Due to the linearity of the governing PDEs, the reconstruction of the volatility surface from European options has been extensively studied in the existing literature. Some efficient and effective algorithms are being developed, see [7,11,14], just to name a few. Hower, in reality, most of the option contracts traded in exchanges are American options. Thus, it is easy to observe a large amount of market American option prices, and hence the reconstruction of the volatility surface from American options is of more practical significance. However, in the case of American options the situation is more complicated. The reason is that due to the possibility of early exercise the set of Black-Scholes PDEs become a set of Black-Scholes partial differential complementarity problems (PDCPs), which obviously is nonlinear. Furthermore, the determination of the volatility function is much more complicated as well. In the framework of infinite-dimensional space, the problem is formulated as an optimal control problem with variational inequalities constraint [1,20]. In the framework of finite-dimensional space, it is formulated as a mathematical program with equilibrium constraint (MPEC) [9]. Due to the highly complicated nature of both of the formulations, they have not yet received the deserved attention in the mathematical finance literature. The difficulty involved in the determination of the volatility function from American options is that the numerical algorithm for solving the optimal control or MPEC is computationally highly demanding. For example, to obtain a decent direction in a gradient-type optimization algorithm, it requires to solve a set of governing PDCPs satisfied by the ideal American option prices. In general, the PDCPs can only be solved numerically by some iterative numerical methods, such as projective successive over-relaxation method (PSOR) [23]. It is computationally much more demanding when compared with solving the PDEs in the case of European options.
In this paper, we will apply the bicubic spline functional approach, where the volatility function σ(s, t) is explicitly represented by a spline with a fixed set of spline knots. The values of the volatility function at the spline knots uniquely determine a local volatility surface. The spline is determined by solving a nonlinear least square problem with the local volatility values at chosen spline knots as unknown variables such that they match the market option prices as closely as possible. By limiting the number of spline knots, the dimension of the resulting optimization problem is typically small. Hence, the computational cost can be reduced significantly.
To solve this nonlinear optimization problem, we will develop a gradient-based routine. During each iteration loop of the routine, a system of PDCPs is to be solved so as to get a decent direction. This is extremely time-consuming if the traditional PSOR method is applied. Penalty method has been used very successfully for solving complementarity problems (cf. [22,4,16,18,26]). The advantages are that it does not introduce any extra or auxiliary variables and the resulting algebraic equations are easily solvable by a conventional numerical method such as Newton method. More importantly, it has been shown (see [26]) that the penalty approach is much more robust and effective than PSOR. On this basis, the reconstruction of the volatility surface from American options can be obtained more effectively. We will adopt the penalty approach and develop a gradient-based algorithm to solving the nonlinear optimization problem in finite-dimensional space. In this approach, the set of PDCPs are approximated by a set of nonlinear penalty PDEs with a penalty term. To solve the set of penalized PDEs effectively, we propose to use a finite difference discretization. This discretization results in a large-scale nonlinear implicit program. We compute the Jacobian of the cost function by using the finite difference technique. With the computed Jacobian, a gradient-based method is applied to solving the nonlinear optimization problem. To verify the effectiveness and robustness of this new algorithm, two experiments are carried out. One is an artificially constructed example, where a series of American puts are computed under the assumption that the underlying asset follows a known absolute diffusion process ( [15]). The algorithm can recover accurately the local volatility surface. Another example is a realistic one. We use real American puts data on the NASDAQ 100 index to test the robustness and effectiveness of the proposed method.
The organization of the paper is as follows. In Section 2, we start with the mathematical formulation of the forward pricing problem and the inverse volatility reconstruction problem. In Section 3, the bicubic spline approximating regularization technique is presented for the inverse problem. In Section 4, the numerical computation of the Jacobian of the cost function is proposed. In Section 5, the computational procedure for solving the proposed optimization problem is developed. Two numerical experiments are presented in Section 6 to illustrate the effectiveness and applicability of the proposed method. The first example is an artificially constructed example, where it is assumed that the underlying asset follows a known absolute diffusion process. The second example is a realistic one involving real American puts data on the NASDAQ 100 index.
2. The formulation of the problem. Under the risk neutral measure, the evolution of the underlying asset price S is governed by the following stochastic differential equation (SDE): where r is the constant risk free interest rate, q is the constant dividend yield, W is a standard Wiener process, and σ (S, t) is the deterministic local volatility function that we wish to determine. It defines the local volatility surface. Let V (σ (S, t) , S * , t * , K, T ) denote the theoretical price of an American put option with striking price K and maturity date T at time t * for the underlying asset with spot price S * . Following the standard arguments (e.g. [23]), it is known that V must satisfy the following partial differential complementarity problem (PDCP): is the payoff function, and the operator L is defined by The American option price can be uniquely determined through solving Problem 2 (numerically) once the volatility function is available.
Assume that the current market prices of vanilla American put options V j m j=1 are given, whereV j =V j (σ (S, t) , S * , t * , K j , T j ) denotes the time t = t * vanilla put option prices with striking price K j and expiry date T j . As there exit corresponding bid pricesV bid j and ask pricesV ask j for any listed option (K j , T j ), we take the average of the bid prices and ask prices as its market price, i.e., With these notations, the reconstruction problem may be stated as: given a finite set of observed American option prices, find a local volatility function such that the theoretical option prices match the observed ones optimally with respect to a prescribed performance criterion. This can be formulated as the following nonlinear optimization problem:

Problem. Given a set of market put option prices
with ω j being a weight (such as trading volume) reflecting the relative importance or significance of different options.
Problem 2 is a nonlinear least square problem. The nonlinear PDCPs (1) do not appear explicitly in the formulation, but the current V j (σ (S, t)) depends on the PDCPs . The dependence of the cost function on the parameter σ (S, t) is nonlinear, implicit and complicated, and thus analytical solution is not available. We need to find effective and robust numerical methods to solve Problem 2.
3. Regularization strategy with bicubic spline approximation. Problem 2 is an inverse problem of the forward pricing Problem 2. Since only a finite number of reported market prices V j m j=1 is available for the determination of the continuous local volatility function σ (S, t), Problem 2 is ill-posed in the sense that a small perturbation in the option price may lead to a large change in the volatility surface. Moreover, due to the lack of bid-ask spreads for some illiquidity options, the observed market prices generally contain error. Furthermore, the theoretical value of an American put option need to be computed numerically, which also inevitably introduces some numerical error. As a consequence, there can be many solutions that yield sufficient agreement with market prices. This situation will affect the accurate pricing and hedging significantly. To overcome this deficiency, we need to impose regularization condition on Problem 2. In this paper, we adopt the bicubic spline regularization technique proposed by T. F. Coleman in [7] to Problem 2.
We apply a 2-dimensional bicubic spline smoothing function to directly approximate the local volatility surface σ (S, t). Let the number of the spline knots p ≤ m be given with S j ,t j p j=1 being a set of fixed spline knots in [0, ∞) × [0, T max ]. Given these p spline knots, an interpolating cubic spline c (S, t,σ) with the natural spline end condition is uniquely determined under the condition c S j ,t j =σ j , j = 1, . . . , p. Here,σ is p−vector defined byσ = (σ 1 , . . . ,σ p ) andσ j . = σ S j ,t j corresponds to the local volatility value at the knot. Then, the local volatility values σ j , j = 1, . . . , p, are determined by fitting the market observable option values as closely as possible. Define In order to incorporate additional a priori information, we can impose lower bound L and upper bound U at the knots. Now, the final regularized version of Problem 2 can be stated formally as follows. L ≤σ ≤ U Problem 3 is a box-constraint nonlinear optimization problem with the local volatilityσ at the spline knots S j ,t j p j=1 . This problem can usually be solved by gradient-based optimization methods, such as conjugate gradient method, trust region method. In the next section, we will concentrate on the computation of the gradient of the cost function (2).

4.
Computation of the gradient of the cost function. It is clear that the gradient of the cost function (2) cannot be expressed analytically, because V j (c (S, t;σ)) can only be obtained through solving Problem 2 numerically. This is computationally very expensive. We first examine the structure of the gradient of the cost function (2). Then, a finite-difference-based numerical implementation is developed for the computation of the gradient.
With this notation, it is obvious that the cost function (2) is f (σ). Hence, the gradient of the cost function is simply where J (σ) is the m × p Jacobian matrix of F with respect toσ. However, we can only get V j (c (S, t;σ)) numerically based on discretization on a predefined grid points in the region [0, S max ]×[0, T max ], where S max denotes a large enough number to ensure the accuracy of the solution. Therefore, the computation of Jacobian matrix J (σ) also depends on the grid points. To this end, we first discretize the computation region [0, S max ] × [0, T max ] into a uniformly N S × N t spaced grids with Thus, J (σ) can be computed based on these grid points, which is expressed as follows: We use forward finite differences to approximate all the partial derivatives in J (σ). More specifically, where e j is a unit vector whose jth component equals to 1, and h is a small positive increment.
5. Numerical method. To numerically solve Problem 3, a gradient-based optimization method is usually used. In this paper, we employ the conjugate gradient method. This method is an iterative method, which requires at least the evaluation of f (σ) and ∇f , and hence it requires the evaluation of V j , j = 1, . . . m, in each iteration. Unlike European options, the pricing of American options is a complementarity problem, which is computationally demanding. Furthermore, Problem 3 is a large scale optimization problem, which is also a time-consuming iterative procedure. Thus, effective numerical methods are required to solve the forward pricing problem (i.e., Problem 2) and the optimization problem (i.e., Problem 3). We first develop a penalty-based method to solve Problem 2 in Subsection 5.1. Then, Subsection 5.2 presents a numerical algorithm for solving Problem 3.

5.1.
Solving the forward pricing problem: A penalty-based approach. Problem 2 is a partial differential complementarity problem. As mentioned in Section 1, the penalty method is an effective and robust approach to solving this complementarity problem. Here, we propose a power penalty approach [22] to solve Problem 2, which results in the following penalized problem: with the terminal condition V (S, T ) = max {K − S, 0} at t = T , and the boundary conditions In (7), V (S, t) denotes the put American option value of an underlying asset with striking K and expiry date T at (S, t), λ is the penalty parameter satisfying λ > 1, k > 0 is the power of the penalty function.
The basic idea of the penalty approach is to force the positive parts of V * − V in (7) to be close to zero as the penalty parameter λ becomes sufficiently large. In this way, the complementarity conditions in (1) are approximately satisfied, which is with C > 0 being a constant. The detailed study of the solvability and convergence property of the power penalty approach can be found in [11]. For convenience, we take, in this paper, k = 1, i.e., linear penalty method. Problem 5.1 is a nonlinear and non-smooth partial differential equation. To numerically solve it, we use the finite difference discretization (4) and (5) described in Section 4. Define V n i . = V (S i , t n ) and V * i = V * (S i ). Then, by applying a Crank-Nicholson finite difference scheme to equation (7), the following discrete equation is obtained. where and ∆S = Smax Ns−1 , ∆t = T Nt−1 . Remark 1. The first row and last row of M have been modified to incorporate the corresponding boundary conditions (8). Also, the terminal condition has been incorporated to V N S −1 as the payoff function of the put option.
It is clear that equation (9) is a nonlinear system, and thus we use the generalized Newton-type iteration method to solve this system, yielding where (V n ) l denotes the lth iteration of V n .

5.2.
Algorithm. We now present a numerical algorithm for solving the nonlinear optimization problem (2). It is a conjugate gradient method, which is detailed as follows.
We comment that the above algorithm can be easily extended to other gradientbased optimization methods, such as trust region method, quasi-Newton method. All these methods are easy to implement once the gradient of the cost function is given. Details of these methods and their convergence properties can be found in [19]. In the next section, we concentrate on testing the effectiveness and robustness of Algorithm 1 through solving a synthetic example and a real market example. 6. Numerical examples. We now describe computational tests with the proposed method for the reconstruction of the local volatility surface σ (S, t) from limited observation data. All computations are performed in Matlab R2010b environment on a dual 3.10GHz and 4GB RAM Intel PC.
6.1. A synthetic example. In order to test the accuracy and robustness of the proposed method, we first investigate an artificially constructed example. We consider American put options with different strike prices and maturity dates. In this example, we choose a given target volatility function σ (S, t). Then, through solving (7) -(8) with this given σ (S, t), we compute the prices of a series of American puts with eleven different striking prices and two different expiry dates as the 'observed' data. The chosen volatility surface, which is independent of time, is defined by σ(S, t) = θ S with θ = 15. With this volatility function, we can easily see that the underlying asset S follows the following diffusion process: dS = (r − q) Sdt + θdW under the risk neutral measure. Let the initial underlying asset price S 0 = 100, the risk-free interest rate r = 0.05, and the dividend rate q = 0.02. We present the above 22 'observed' pricesV j = V j of American put options in Table 1.
The To investigate the robustness of our algorithm, we add a white noise to the 'observed' prices in Table 1. For j = 1, . . . , 22, we set where j is a random variable drawn from the standardized normal distribution. We then use the newV j to reconstructing the local volatility. The volatility function are plotted in Figure 3. It can be seen that the reconstructed volatility surface also     Table 2. Observed American option prices fits well to the exact in our interested regions, except for some large deviation in few areas. We comment that the deviation in these areas arises from occurrence of some unrealistic extreme values of the white noise j , since j is arbitrarily drawn from a normal distribution. After eliminating these extreme values, the reconstructed volatility surface is recovered very well, which verifies the robustness of the new method.
6.2. A real market example. We now test the performance of our new method with real market data. The example is an American put option from 30/10/2000, whose underlying asset is NASDAQ 100 index. It is the same data given in [20].
To reconstruct the volatility function, we use the quoted prices with maturities  Table  2.
In Table 2, the weight ω j is set to be the trading volume of option V j , reflecting its relative importance among all the observed options. The NASDAQ 100 index on 30/10/2000 is S 0 = 76.7656. The other default parameters are r = 0.05 and q = 0. We also set S max = 200, N t = 101, N S = 51.    Figure 4. To test the effectiveness of our new algorithm, we compare the the option values generated from the recovered volatility function with the real market options values in Figures 5 -7. It can be observed that the generated option values coincide with the real one very well. This shows that the new method is very reliable and effective. 7. Conclusions. In this paper, we developed a new method to reconstruct the smooth volatility surface from the market observed American option values. This volatility reconstruction problem was first formulated as a nonlinear least square problem. To overcome the ill-poseness of the problem, a bicubic spline smoothing regularized technique was applied, where the unknown variables are confined to the   Figure 7. Comparison between the option prices computed from reconstructed volatility surface and observed American option prices with expiry date T 3 spline knots. For the regularized problem, we presented a gradient-based optimization method. When implementing the optimization method, we designed a forward finite difference approximation to compute the Jacobian matrix. In each iteration of this optimization method, the option values at spline knots were obtained through numerically pricing American options. To compute the values of American options, we proposed a penalty approach. Finally, a synthetic test and a real market test were solved by using the new method. Numerical results verified that the new method is effective and robust.