A simple and efficient numerical method for pricing discretely monitored early-exercise options
aa r X i v : . [ q -f i n . C P ] J un A SIMPLE AND EFFICIENT NUMERICAL METHOD FORPRICING DISCRETELY MONITORED EARLY-EXERCISEOPTIONS
MIN HUANG AND GUO LUO
Abstract.
We present a simple, fast, and accurate method for pricing a va-riety of discretely monitored options in the Black-Scholes framework, includ-ing autocallable structured products, single and double barrier options, andBermudan options. The method is based on a quadrature technique, and it em-ploys only elementary calculations and a fixed one-dimensional uniform grid.The convergence rate is O (1 /N ) and the complexity is O ( MN log N ), where N is the number of grid points and M is the number of observation dates. Introduction
Exotic options are commonly traded throughout the world. Many popular ex-otic options are path-dependent and have early-exercise features. These optionscan often be priced using analytical formulas if they are continuously monitored(e.g. barrier options). In practice, however, most path-dependent exotic optionsare discretely monitored [4], in which case they need to be priced using numericaltechniques. Due to the complicated structures of these options, traditional pricingmodels based on Monte-Carlo simulations and finite difference methods are oftentoo time-consuming to be useful in practical situations. More recent pricing meth-ods based on advanced mathematical techniques, on the other hand, tend to be moreefficient (e.g. [14, 12, 10, 11]), but for many financial institutions, these methods areoften too difficult to understand and to properly implement. To strike a balance be-tween model performance and practical utility, we propose a new quadrature-basedmethod that is much faster and more accurate than the traditional Monte-Carloand PDE methods, yet at the same time is easy to understand and to implement.We will first give a brief review of the types of products considered, as well as thequadrature-based pricing model which is the foundation of our work. Then we willexplain our method and provide numerical examples.1.1.
Autocallable Structured Products.
Autocallable structured products be-long to a class of exotic options with early-exercise features. Many different typesof autocallable products have been created and traded in financial markets, andthey have become increasingly popular in recent years. We refer to the Appendixof [9] for a description of the main features of various autocallable products.We will consider a very common autocallable product with discrete observationdates. At each observation date there is a pre-specified barrier level. If the price of
Key words and phrases.
Discrete option pricing, quadrature method, autocallable structuredproduct, single and double barrier option, Bermudan option.Disclaimer: The statements made and opinions expressed here are solely our own, and do notreflect the views of China Mechants Bank or its employees and affiliates. the underlying asset is greater (less) than or equal to the barrier level (dependingon the terms of the product), the option is exercised and a pre-specified fixed-ratereturn is paid. If the asset price is below (above) the barriers at all observationdates, the option is never exercised and the investor receives a negative returnat maturity. In addition, autocallable products may have a knock-in feature. Inthis case, if the option is never exercised, the negative return the investor receivesdepends on whether the asset price at maturity reaches a pre-specified knock-inlevel. While the value of a continuously monitored autocallable product has a sim-ple closed-form solution, the value of a discretely monitored autocallable productcannot be calculated easily. In the discrete case, there exist analytical solutions interms of multiple integrals, cf. [25, 26, 9]. The numerical calculation of these inte-grals, however, can become prohibitive if the number of observation dates exceedsfive. In practice, discretely monitored autocallable products are commonly pricedusing Monte-Carlo simulations. This method is straightforward, but convergenceis usually slow and acceleration techniques such as variance reduction are oftenneeded (cf. [2, 13]). Another popular method for pricing discretely monitoredautocallable products is to solve the governing Black-Scholes partial differentialequation (PDE) using finite difference method (cf. [9]). Assuming a second-ordercentral difference approximation in space, the overall convergence rate of a typicalfinite difference based pricing method is O (( δx ) + δt ) if the explicit forward Eulermethod is used in time, and O (( δx ) + ( δt ) ) if the implicit Crank-Nicolson methodis used. Since two-dimensional grids are needed for finite difference methods, com-putational complexities are at least of order 1 / ( δxδt ). In addition, since the payoffsof autocallable products are discontinuous (in asset price), additional care (such assmoothing of payoff functions) must be taken to ensure the accuracy of any finitedifference approximations. Remark 1.1.
Some autocallable products may have payoffs at maturity that are ofthe same type as that of European vanilla options. These products can be effectivelyviewed as combinations of autocallable products and barrier options (see below).
Discrete Barrier Options.
Barrier options are among the most populartypes of exotic options. A barrier option may be activated (knock-in option) ordeactivated (knock-out option) when the price of the underlying asset crosses cer-tain barrier levels. Barrier options may be discretely or continuously monitored.A single barrier option has one barrier at each observation date, while a doublebarrier option has two barriers at each observation date. The final payoff of abarrier option (if it is active at maturity) may be of the same type as that of avanilla option or that of a digital option. A special type of barrier option has aconstant amount of cash as the final payoff, and such an option is called a touchoption. Most barrier options have time-independent barrier levels, but options withtime-dependent barrier levels have also been studied [7].Similar to the case of autocallable structured products, there are closed-formsolutions for discrete barrier options in terms of multiple integrals [23], but suchsolutions are often difficult to evaluate directly. In practice, discrete barrier optionscan be priced using Monte-Carlo simulations or standard binomial tree methods,but these methods are usually slow [5]. Other methods that have been proposed toprice discrete barrier options include continuity correction approximations [4, 27],Wiener-Hopf methods [14], adaptive mesh methods [1], Hilbert transform meth-ods [12], finite element methods [16], Fourier-cosine series expansion methods [10],
SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 3 and quadrature methods [3, 6]. These methods, while useful in certain contexts,have not been as widely used as the traditional Monte-Carlo and finite differencemethods, usually due to their complexity.1.3.
Bermudan Options.
Bermudan options are discrete versions of Americanoptions. A Bermudan option can be exercised at any of the prescribed observationdates, and the payoff is of the same type as that of a vanilla option. Similarto discrete barrier options, Bermudan options can be priced using Monte-Carlosimulations [17], Hilbert transform methods [11], Fourier-cosine series expansionmethods [10], and quadrature methods [22, 21, 6].1.4.
Overview of the Quadrature Method.
Among the various methods pro-posed to price discretely monitored options, the quadrature method is particularlyappealing because of its high efficiency and accuracy. The method has been appliedto discrete barrier options and Bermudan options [3, 22, 21, 6]. The main idea isto solve for option values at each observation date via backward induction in time.The risk-neutral valuation formula is expressed as a single integral, which is thenevaluated numerically to produce the option price. Specifically, let V denote thevalue of the option, S the value of the underlying asset, r the risk-neutral interestrate, t , . . . , t M the observation dates, and E the risk-neutral expectation. If theunderlying asset S does not trigger an early exercise at t m , we have V ( t m , S ) = e − r ( t m +1 − t m ) E (cid:2) V ( t m +1 , · ) | S (cid:3) = e − r ( t m +1 − t m ) Z ∞ V ( t m +1 , y ) f ( y | S ) dy, where f ( y | S ) is a probability density function whose form depends on the modelof the underlying asset. If S triggers an early exercise at t m , on the other hand,the option price V ( t m , S ) would be equal to a prescribed value. The integral abovecan be calculated using FFT [21, 22] or Fast Gauss Transform [6]. Since, however, V ( t m +1 , y ) is discontinuous (in y ) for autocallable products and barrier options,and non-differentiable (in y ) for Bermudan options, care must be taken to ensurethe accuracy of the numerical evaluation of the integral. While several shifted ornonuniform grids have been designed in previous studies to address this difficulty[21, 6], the problem becomes particularly challenging when multiple discontinuitiesare present at each observation point, for instance in the case of double barrieroptions with time-dependent barrier levels.1.5. Motivation of Our Work.
Although a good number of advanced techniqueshave been proposed to improve pricing models’ accuracy and efficiency, in mostpractical situations, simple methods that are more cost effective are usually pre-ferred over their more sophisticated counterparts. The reason for this is twofold.First of all, when data quality is not high enough, sophisticated models are notnecessarily beneficial. For instance, interest rates and volatilities are crucial com-ponents of nearly every pricing model, but they need to be estimated from availablemarket data. If the estimated parameters contain large errors, which is not uncom-mon in products sold in emerging markets, any advantages gained from the useof sophisticated models may be (more than) offset by these errors, making thesimpler models more attractive. Secondly, implementations of pricing models usu-ally involve staff members from multiple business departments, and the resultingproducts often need active maintenance and updates. As a result, models that are
MIN HUANG AND GUO LUO too complicated in nature may hinder effective business communications, whichincreases maintenance costs and operational risks. In view of these considerations,it is not difficult to see why traditional Monte-Carlo and PDE methods are stillamong the most popular methods in the valuation of discretely monitored options,even though their computational costs are already high enough to adversely impacttheir applicability in business.In view of the practical concerns mentioned above, we propose a new quad-rature method to price the aforementioned discretely monitored options in theBlack-Scholes framework. The convergence rate of our method is O (1 /N ) andthe complexity is O ( M N log N ), where N is the number of grid points and M isthe number of observation dates. The performance of our method is on par withprevious quadrature-based methods such as the CONV method [21], but it is morestraightforward, and is better suited for products with multiple discontinuities. Ourmethod differs from other quadrature methods mainly in three aspects. First, wework with probability density functions directly instead of using characteristic func-tions or Toeplitz matrices. Secondly, we use only a fixed one-dimensional uniformgrid to compute all integrals. Thirdly, we utilize explicit Black-Scholes formulasto improve the accuracy of the calculations. Due to these novel modifications, ourmethod is very easy to implement, and is capable of handling sophisticated productssuch as double-barrier options with time-dependent barrier levels.1.6. Organization of the Paper.
The rest of the paper is organized as follows.Section 2 specifies the class of (discrete) option pricing problems that our quadra-ture method is intended to solve, and Section 3 presents the main recursion formulafor our method. After detailing the implementation of our method in Section 4–5,we summarize the algorithm in Section 6 and then present numerical examples inSection 7. The Appendix collects a few useful theoretical results, which lay thefoundation of a class of (discrete) option pricing algorithms (including the one de-scribed here) but which, to our best knowledge, do not seem to have been rigorouslyproved or even properly formulated in the literature. We supply the proofs here inthe hope that they would be useful to interested readers.2.
Basic Assumptions
We assume that the price of the underlying asset S ( t ) satisfies the followingstochastic differential equation in the risk-neutral measure:(2.1) dS ( t ) = (cid:2) r ( t ) − q ( t ) (cid:3) S ( t ) dt + σ ( t ) S ( t ) dW ( t ) , where r ( t ) is the risk-neutral interest rate, q ( t ) is the yield rate, σ ( t ) is the volatility,and W ( t ) is the Wiener process.In practice, interest rates are always time dependent. Yield rates for FX productsare simply foreign interest rates, and for other types of products they may be impliedfrom futures prices. Thus yield rates are usually time dependent as well. Impliedvolatilities are time dependent, whereas historical volatilities can often be taken asconstant.The solution to (2.1) is(2.2) S ( t ) = S ( t ) exp (cid:26)Z tt (cid:2) r ( s ) − q ( s ) − σ ( s ) (cid:3) ds + Z tt σ ( s ) dW ( s ) (cid:27) , SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 5 where t is the present date. We consider a discretely monitored option with obser-vation dates t , . . . , t M , where the last observation date t M is the maturity date. Itfollows from (2.2) that each S ( t m ) for 1 ≤ m ≤ M has the lognormal distribution(2.3) S ( t m ) ∼ S ( t ) exp (cid:26)Z t m t (cid:2) r ( s ) − q ( s ) − σ ( s ) (cid:3) ds + (cid:16)Z t m t σ ( s ) ds (cid:17) / Z (cid:27) , where Z denotes the standard normal distribution. Now define r m = Z t m t m − r ( s )∆ t m ds, q m = Z t m t m − q ( s )∆ t m ds, σ m = Z t m t m − σ ( s )∆ t m ds, for 1 ≤ m ≤ M where ∆ t m = t m − t m − , and define piecewise constant functions˜ r ( t ) = r m , ˜ q ( t ) = q m , and ˜ σ ( t ) = σ m , for t m − < t ≤ t m . For the process d ˜ S ( t ) = (cid:2) ˜ r ( t ) − ˜ q ( t ) (cid:3) ˜ S ( t ) dt + ˜ σ ( t ) ˜ S ( t ) dW ( t ) , since Z t m t r ( s ) ds = m X n =1 Z t n t n − r ( s ) ds = m X n =1 r n ∆ t n = Z t m t ˜ r ( s ) ds, Z t m t q ( s ) ds = m X n =1 Z t n t n − q ( s ) ds = m X n =1 q n ∆ t n = Z t m t ˜ q ( s ) ds, and Z t m t σ ( s ) ds = m X n =1 Z t n t n − σ ( s ) ds = m X n =1 σ n ∆ t n = Z t m t ˜ σ ( s ) ds, it follows that ˜ S ( t m ) has the same distribution as S ( t m ) in (2.3) for each m . Sincethe value of the option depends only on probability distributions of the asset priceat observation dates, the option value remains the same if we replace the process S by the process ˜ S . In other words, we may safely assume that r ( t ) , q ( t ), and σ ( t )are piecewise constant functions. Thus in what follows, we shall assume r ( t ) = r m , q ( t ) = q m , and σ ( t ) = σ m , for t m − < t ≤ t m . Consider now a general class of discretely monitored options with barriers. Sincethe sum of a knock-in barrier option and a knock-out barrier option with the sameobservation dates and barrier levels is a vanilla option (or a digital option if thebarrier options are digital), to study the pricing of these discretely monitored op-tions, it suffices to consider a knock-out barrier option which ceases to exist whenbarrier levels are crossed. To this end, assume that(A) The option has two strike prices K − m , K + m ∈ [0 , ∞ ], with K − m ≤ K + m , at eachobservation date t m , m = 1 , , . . . , M .(B) The option is exercised if S ≤ K − m or S ≥ K + m at some t m , and the payoffsare given by a − m S + b − m (if S ≤ K − m ) and a + m S + b + m (if S ≥ K + m ), respectively,for some a ± m , b ± m ∈ R .(C) The final payoff at maturity is V ( t M , S ) = a M S + b M , for K − M < S < K + M . MIN HUANG AND GUO LUO
These assumptions are general enough to cover a wide class of discretely moni-tored options, such as the ones mentioned in the introduction. For instance, com-mon up-and-out autocallable products would have1 ≤ m ≤ M : 0 < K + m < ∞ , K − m = 0 , a + m = 0 , b + m > m = M : a M = 0 , b M < . Down-and-out put barrier options would have1 ≤ m ≤ M − K + m = ∞ , < K − m < ∞ , a − m = 0 , b − m = 0; m = M : K − M = 0 , < K + M < ∞ , a M = − , b M = K + M ,a + M = b + M = 0 . Double barrier knock-out call options would have1 ≤ m ≤ M − < K ± m < ∞ , a ± m = 0 , b ± m = 0; m = M : K + M = ∞ , < K − M < ∞ , a M = 1 , b M = − K − M ,a − M = b − M = 0 . Bermudan put options with strike K would have1 ≤ m ≤ M : K − m : the unique solution of K − K − m = V ( t m , K − m ) ,K + m = ∞ , a − m = − , b − m = K ; m = M : a M = 0 , b M = 0 . We will give a proof of the uniqueness of K ± m for Bermudan options in the Appendix.To summarize, our basic assumptions are(1) The underlying asset price S follows a geometric Brownian motion withpiecewise constant interest rates, yield rates, and volatilities.(2) There are finitely many observation points, and two exercise levels (possibly ∞ ) at each observation point. If S is above the upper exercise level or belowthe lower exercise level at any observation point, the option is exercised andthe payoff is a linear function in S .(3) At maturity, if S is between the two exercise levels, a payoff is incurredwhich is also a linear function in S .3. Outline of the method
Let V ( t, S ) denote the value of the option (as a function of asset price S ) at anytime t , and let V m ( S ) = V ( t m , S ) , m = 0 , , . . . , M, denote the value of the option at the observation dates. Our goal is to find V ( S ( t )),and our strategy is to use backward induction in time. Since V M ( S ) is piecewiselinear in S , V M − ( S ) has a simple explicit expression. For each m = M − , . . . , V m − ( S ) as the risk-neutral expectation of V m ( S ) for K − m − < S < K + m − ,and as a ± m − S + b ± m − otherwise. The expectation is given by an explicit integraland is calculated numerically. The core of the quadrature method is the calculationof M − τ m = 12 σ m ∆ t m = 12 σ m ( t m − t m − ) . SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 7
For each 1 ≤ m ≤ M −
1, note that S ( t m ) has a lognormal distribution as in (2.3).The relevant probability density functions are known to be [19](3.1) ρ m ( y, S ) = 12 √ πτ m y exp (cid:26) − τ m (cid:16) log yS − σ m (cid:2) r m − q m − σ m (cid:3) τ m (cid:17) (cid:27) . For simplicity of notations we define K +0 = ∞ and K − = 0. By the fundamentaltheorems of asset pricing, we have the risk-neutral pricing formula [24]: V m − ( S ) = e − r m τ m /σ m E (cid:2) V m ( · ) | S (cid:3) = e − r m τ m /σ m Z ∞ V m ( y ) ρ m ( y, S ) dy (3.2)= e − r m τ m /σ m √ πτ m Z ∞ y V m ( y ) exp (cid:26) − τ m (cid:16) log yS − σ m (cid:2) r m − q m − σ m (cid:3) τ m (cid:17) (cid:27) dy, for K − m − < S < K + m − and 1 ≤ m ≤ M −
1. By Assumption (B) from Section 2,we also have(3.3) V m − ( S ) = ( a − m − S + b − m − , S ≤ K − m − a + m − S + b + m − , S ≥ K + m − . To further study the formulas (3.2) and (3.3), we first recall some classical resultson the pricing of binary options.
Lemma 3.1.
Let
K > , and let χ A denote the characteristic function of a set A .Consider an option with no early-exercise features.(1) If the option has payoff ˆ V m ( y ) = χ [ K, ∞ ) y , then ˆ V m − ( S ) = V am ( S, K, .(2) If ˆ V m ( y ) = χ (0 ,K ] y , then ˆ V m − ( S ) = V am ( S, K, − .(3) If ˆ V m ( y ) = χ [ K, ∞ ) , then ˆ V m − ( S ) = V bm ( S, K, .(4) If ˆ V m ( y ) = χ (0 ,K ] , then ˆ V m − ( S ) = V bm ( S, K, − .The functions V am and V bm are defined as V am ( S, K, ǫ ) = e − q m τ m /σ m SN ( ǫd ) , V bm ( S, K, ǫ ) = e − r m τ m /σ m N ( ǫd ) , where N is the cumulative normal distribution function, and d = 1 √ τ m (cid:16) log SK + 2 σ m (cid:2) r m − q m + σ m (cid:3) τ m (cid:17) , d = d − √ τ m . Proof.
By definition V am is the value of an asset-or-nothing option, and V bm is thevalue of a cash-or-nothing option. The valuation formulas are just standard resultsfor binary options [18]. (cid:3) Remark 3.2.
The standard Black-Scholes formulas in Lemma 3.1 ignore the possi-ble effects of volatility smiles. If such effects need to be taken into account, one mayamend the definitions of V am and V bm (as given in the Lemma) by incorporating suit-able vega-induced correction terms [15] . For instance, the value of a cash-or-nothingcall option in the presence of volatility smiles would become V smile = V no smile − ∂V vanilla ∂σ ∂σ∂K . MIN HUANG AND GUO LUO
We can use Lemma 3.1 to obtain an explicit formula for the value of V M − ( S ).By Assumption (B)–(C) from Section 2, we have(3.4) V M ( S ) = a − M S + b − M , S ≤ K − M a M S + b M , K − M < S < K + M a + M S + b + M , S ≥ K + M . Without loss of generality we may assume 0 < K ± M < ∞ , since otherwise we maychoose some arbitrary 0 < K ± M < ∞ and set a ± M = a M , b ± M = b M . Proposition 3.3.
The value of the option at t M − is given by ˜ V M − ( S ) = a − M V aM ( S, K − M , −
1) + b − M V bM ( S, K − M , − a M (cid:2) V aM ( S, K − M , − V aM ( S, K + M , (cid:3) + b M (cid:2) V bm ( S, K − M , − V bM ( S, K + M , (cid:3) + a + M V aM ( S, K + M ,
1) + b + M V bM ( S, K + M , , for K − M − < S < K + M − .Proof. Clearly, the option from t M − to t M is equivalent to a linear combinationof binary options consisting of two put options with strike K − M , two call optionswith strike K − M and four call options with strike K + M . The result then follows fromLemma 3.1. (cid:3) With the aid of (3.2) and Proposition 3.3, we may write the main recursion ofour quadrature method as follows.
Proposition 3.4.
Let ˜ V M = V M be defined in (3.4) , and ˜ V m − be given by thefollowing recursion formula: ˜ V m − ( S ) = e − r m τ m /σ m Z K + m K − m ˜ V m ( y ) ρ m ( y, S ) dy + a + m V am ( S, K + m , b + m V bm ( S, K + m ,
1) + a − m V am ( S, K − m , −
1) + b − m V bm ( S, K − m , − , for ≤ m ≤ M . Then we have ˜ V m ( S ) = V m ( S ) for K − m < S < K + m and ≤ m ≤ M . In particular, ˜ V ( S ( t )) = V ( S ( t )) .Proof. This is merely the classical recursion formula for quadrature methods spe-cialized to the Black-Scholes model. To prove the formula, we only need to show(3.6) V m ( S ) = ˜ V m ( S ) χ ( K − m ,K + m ) + ( a + m S + b + m ) χ [ K + m , ∞ ) + ( a − m S + b − m ) χ (0 ,K − m ] , for all 0 ≤ m ≤ M . By assumption (3.6) is true for m = M . Assume now (3.6)holds for some 1 ≤ m ≤ M . Substituting the equation into (3.2), applying Lemma3.1, comparing the result with (3.5) and using (3.3), we observe that (3.6) holds for m −
1. The result then follows from induction. (cid:3)
Remark 3.5.
Recursion formula (3.5) lies at the heart of our quadrature methodand distinguishes our method from other quadrature methods, which are primarilybased on (3.2) or one of its many variants. The significance of the formula (3.5)
SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 9 lies in the fact that it makes explicit use of Black-Scholes formulas to separate theexpectation integral E (cid:2) V m ( · ) | S (cid:3) into a “quadrature part” F m − ( S ) = e − r m τ m /σ m Z K + m K − m ˜ V m ( y ) ρ m ( y, S ) dy, and an “early-exercise part” E m − ( S ) = a + m V am ( S, K + m ,
1) + b + m V bm ( S, K + m , a − m V am ( S, K − m , −
1) + b − m V bm ( S, K − m , − . Since the function ˜ V m ( S ) is smooth for S ∈ ( K − m , K + m ) (in fact for all S ∈ (0 , ∞ ) , aswe will show below), the integral F m − ( S ) can be evaluated accurately and efficientlyusing a high-order quadrature method such as Simpson’s rule. In contrast, theintegrand V m ( S ) in the original recursion formula (3.2) is discontinuous on (0 , ∞ ) (in either V m itself or in its first derivative); this makes the accurate evaluationof the expectation integral a difficult and challenging task. Although (3.5) appliesspecifically to the Black-Scholes model, the same idea can be used for other assetprice models, as long as a suitable analytical formula (exact or approximate) can befound for the probability density function ρ m ( y, S ) and early-exercise part E m − ( S ) . With ˜ V M − ( S ) given in Proposition 3.3, Proposition 3.4 implies that we mayapply (3.5) successively to obtain ˜ V M − ( S ) , ˜ V M − ( S ) , . . . , ˜ V ( S ). The value of theoption is equal to ˜ V ( S ( t )).4. Details of Implementation
In (3.5), ˜ V m − ( S ) is written as a sum of explicit functions and an integral, namely˜ V m − ( S ) = a + m V am ( S, K + m ,
1) + b + m V bm ( S, K + m , a − m V am ( S, K − m , −
1) + b − m V bm ( S, K − m , −
1) + F m − ( S ) , where F m − ( S ) = e − r m τ m /σ m Z K + m K − m ˜ V m ( y ) ρ m ( y, S ) dy. We may truncate the integral by replacing its upper and lower bounds by L + m = min (cid:8) K + m , S ( t ) C (cid:9) , and L − m = max (cid:8) K − m , S ( t ) /C (cid:9) , respectively, where C > C = 10 σ √ t M − t + (cid:0) σ (cid:1) ( t M − t ) , where σ = max ≤ m ≤ M σ m , is sufficient to reduce the truncation errors to round-off level. Heuristically this is clear from (2.3), which suggests that the chance that S ( t m ) move outside the range ( S ( t ) /C, S ( t ) C ) is negligibly small. The rigorousderivation of the error bounds can be obtained using a recursive argument, as willbe explained in the Appendix.Now we consider the truncated integral˜ F m − ( S ) = e − r m τ m /σ m Z L + m L − m ˜ V m ( y ) ρ m ( y, S ) dy. If K − m ≥ S ( t ) C or K + m ≤ S ( t ) /C , then by convention the integral is zero. Thusin what follows we shall assume K − m < S ( t ) C and K + m > S ( t ) /C . Let B ± m = log L ± m S ( t ) , and denote S = S ( t ) e x , y = S ( t ) e z ,α m = 1 σ m (cid:2) r m − q m − σ m (cid:3) , β m = 1 σ m (cid:2) r m − q m − σ m (cid:3) + 2 r m σ m ,u m ( x ) = ˜ V m ( S ( t ) e x ) , w m ( x ) = exp n − x τ m − α m x o . The truncated integral can be rewritten as(4.2) ˜ F m − ( S ( t ) e x ) = e − β m τ m √ πτ m Z B + m B − m w m ( x − z ) u m ( z ) dz. One can show by differentiating (4.2) that ˜ F m , and thus ˜ V m and u m , are smoothfunctions in x . This means we can compute the integrals efficiently using a high-order quadrature such as Simpson’s rule.In general, B ± m are different for different values of m , so they cannot all beplaced on one grid. Now we choose a uniform grid x = { x , x , . . . , x N } , where x = − log C and x N = log C . Let h = x N − x N − CN − . For each m , let p − m = min (cid:8) i : x i ≥ B − m (cid:9) , p + m = max (cid:8) i : x i < B + m (cid:9) , where by definition p − m ≥ p + m < N . Since we will use Simpson’s rule whichrequires an odd number of grid points, we define p = ( p + m − p − m ) mod 2 , and rewrite (4.2) as˜ F m − ( S ( t ) e x ) = e − β m τ m √ πτ m (cid:18)Z x p + m + p x p − m w m ( x − z ) u m ( z ) dz (4.3) + Z x p − m B − m w m ( x − z ) u m ( z ) dz + Z B + m x p + m + p w m ( x − z ) u m ( z ) dz (cid:19) . For each 2 ≤ m ≤ M −
1, we will compute ˜ F m − ( S ( t ) e x ) for all x ∈ (cid:8) x , x , . . . , x N , B − m − , B + m − , ξ − m − , ξ + m − (cid:9) , where ξ − m − = 12 ( x p − m − + B − m − ) , ξ + m − = 12 ( x p + m − + p + B + m − ) . For m = 1 we only need to compute ˜ F m − ( S ( t ) e x ) for x = 0, since the value ofthe option is given by ˜ V ( S ( t )). SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 11
Computation of the first integral in (4.3) . To compute the first integralin (4.3) using Simpson’s rule, we let U m ( i ) = u m ( x i ) , i = p − m , p + m + p u m ( x i ) , i = p − m + 1 , p − m + 3 , . . . , p + m + p − u m ( x i ) , i = p − m + 2 , p − m + 4 , . . . , p + m + p − . The integral is discretized as(4.4) Z x p + m + p x p − m w m ( x − z ) u m ( z ) dz = h p + m + p X i = p − m w m ( x − x i ) U m ( i ) + O ( h ) , since Simpson’s rule is of order 4 [8]. Note that U m ( i ) is known from the pre-vious step (or by Proposition 3.3 for m = M −
1) for all i = 1 , , . . . , N . For x ∈ { B ± m − , ξ ± m − , } , the sum (4.4) can be computed directly with complexity O ( N ). For all grid points x ∈ { x , x , . . . , x N } , on the other hand, the sum (4.4)can be computed altogether using FFT with complexity O ( N log N ). This latterfact is crucial to the efficient implementation of our quadrature method and is aconsequence of the following simple observation. Proposition 4.1.
Define (2 N − -periodic grid functions ˆ z, ˆ U m , and ˆ F m by ˆ z ( i ) = z i = − C + ( i − h, ≤ i ≤ N − , ˆ U m ( i ) = , ≤ i < p − m U m ( i ) , p − m ≤ i ≤ p + m + p , p + m + p < i ≤ N − , ˆ F m = F − n F (cid:0) w m (ˆ z ) (cid:1) F ( ˆ U m ) o , where F and F − denote the discrete Fourier transform and the inverse discreteFourier transform of size N − , respectively. Then Z x p + m + p x p − m w m ( x j − z ) u m ( z ) dz = h F m ( j + N ) + O ( h ) , for all ≤ j ≤ N . The above discrete Fourier transforms and inverse discreteFourier transform can be calculated using FFT, and the total computational com-plexity is O ( N log N ) .Proof. We consider the discrete convolution G m ( j ) = N − X i =1 w m ( z j − i ) ˆ U m ( i ) , for j ∈ Z . Note that by definition, z j + N − i = − C + ( j − i + N − h = ( j − i ) h = x j − x i , for all 1 − N ≤ j − i ≤ N −
1. Thus for 1 ≤ j ≤ N , we have G m ( j + N ) = N − X i =1 w m ( z j + N − i ) ˆ U m ( i )= p + m + p X i = p − m w m ( z j + N − i ) U m ( i ) = p + m + p X i = p − m w m ( x j − x i ) U m ( i ) . Therefore (4.4) with x = x j can be written as(4.5) Z x p + m + p x p − m w m ( x j − z ) u m ( z ) dz = h G m ( j + N ) + O ( h ) . The discrete convolution G m can be calculated using FFT as G m = F − n F (cid:0) w m (ˆ z ) (cid:1) F ( ˆ U m ) o = ˆ F m , with a complexity of O ( N log N ) [8]. (cid:3) Remark 4.2.
Our method differs from other well-known FFT-based methods (suchas [21, 22] ) in that we express the discrete quadrature rule (4.4) directly in termsof discrete Fourier transforms, instead of applying continuous Fourier transformto the integral and then discretizing the Fourier integrals (in other words, we haveexchanged the order of Fourier transform and discretization). The direct applica-tion of the discrete Fourier transform (to the discrete quadrature rule) not onlyeliminates the need for artificially-introduced damping factors, which are requiredfor the existence of the continuous Fourier transforms, but also eliminates the needfor additional specially-designed computational grids which are required to satisfyNyquist relations. This enables us to carry out the main recursion (3.5) on a fixeduniform grid, without any additional artificial parameters.
Computation of the last two integrals in (4.3) . The last two integrals in(4.3) are calculated in similar ways using Simpson’s rule. First, note that we mayuse Proposition 3.3 to calculate ˜ V M − ( S ( t ) e x ) for x ∈ { B ± M − , ξ ± M − } . Generallyall four points are needed if the option has two barriers, and only two are needed ifthe option has one barrier. For each 2 ≤ m ≤ M −
1, assume u m ( x ) = ˜ V m ( S ( t ) e x )has been calculated for x ∈ { B ± m , ξ ± m } . The last two integrals in (4.3) are calculatedusing Simpson’s rule as follows: Z x p − m B − m w m ( x − z ) u m ( z ) dz = 16 ( x p − m − B − m ) (cid:2) w m ( x − B − m ) u m ( B − m )(4.6) + 4 w m ( x − ξ − m ) u m ( ξ − m ) + w m ( x − x p − m ) u m ( x p − m ) (cid:3) + O ( h ) , Z B + m x p + m + p w m ( x − z ) u m ( z ) dz = 16 ( B + m − x p + m + p ) (cid:2) w m ( x − B + m ) u m ( B + m )(4.7) + 4 w m ( x − ξ + m ) u m ( ξ + m ) + w m ( x − x p + m + p ) u m ( x p + m + p ) (cid:3) + O ( h ) . SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 13 Finding Optimal Exercise Prices for Bermudan options
Unlike autocallable products and barrier options, Bermudan options do not havepre-specified exercise levels. Instead, one needs to solve for K ± m from the equations˜ V m ( K + m ) = K + m − K, for call options and ˜ V m ( K − m ) = K − K − m , for put options, where ˜ V m is determined by (3.5). For simplicity we assume theyield rates q m ≥
0, which is almost always the case in practice. We will demonstratehow to find K − m , as the same procedure applies to K + m . Let p = min (cid:8) i : ˜ V m ( S ( t ) e x i ) > K − S ( t ) e x i (cid:9) . If p = 1 there is no early exercise, so K − m = 0. Otherwise we have S ( t ) e x p − ≤ K − m < S ( t ) e x p by Corollary 8.3. The value of K − m can be found using classicalroot-finding methods such as the bisecting method or the secant method. Notethat the bisecting method is guaranteed to converge by Corollary 8.3, and it takes O (log N ) steps to reduce the error of the approximate root to an order of O ( h ).Since the cost for calculating ˜ V m at one point using (4.1) is O ( N ), the total cost forfinding the optimal exercise price is O ( N log N ). The secant method is superlinearand converges faster than the bisecting method, though its error estimates are notas straightforward. 6. Summary of the Algorithm
We summarize our algorithm as follows: Define the functions V am , V bm as in Lemma 3.1 Define the function ˜ V M − as in Proposition 3.3 if option style is Bermudan then Calculate K ± M − as in Section 5 end if Calculate p ± M − and p to find the bounds of integration in (4.3) Use Proposition 3.3 to compute ˜ V M − ( S ( t ) e x ) for x ∈ { B ± M − , ξ ± M − } , andassign their values to v ± , v ± respectively Define a vector S as S ( i ) ← S ( t ) e x i for i = 1 , , . . . , N Define a vector y as y ( i ) ← ˜ V M − ( S ( i )) for i = 1 , , . . . , N for m = M − do Let ˆ z and ˆ U m be as defined in Proposition 4.1 ˆ F m ← F − n F (cid:0) w m (ˆ z ) (cid:1) F ( ˆ U m ) o Define (or redefine) the vector Y as Y ( j ) ← h ˆ F m ( j + N ) for j = 1 , , . . . , N Use (4.6), (4.7), and the values of v ± , v ± , y ( p − m ) , y ( p + m + p ) to compute thelast two integrals in (4.3) at x i for i = 1 , , . . . , N , and assign their values to Y and Y if option style is Bermudan then Calculate K ± m − as in Section 5 end if Calculate p ± m − and p to find the bounds of integration in (4.3) Compute ˜ V m − ( S ( t ) e x ) for x ∈ { B ± m − , ξ ± m − } using (4.1) (where the firstintegral in (4.3) is computed using (4.4), and the last two integrals using(4.6), (4.7), and the existing values of v ± , v ± , y ( p − m ) , y ( p + m + p )), and assigntheir values to v ± , v ± respectively y ← Y + Y + Y + a + m V am ( S , K + m ,
1) + b + m V bm ( S , K + m , a − m V am ( S , K − m , −
1) + b − m V bm ( S , K − m , − , as in (4.1) (note that y now stores ˜ V m − ( S ( i )) for i = 1 , , . . . , N ) end for Compute ˜ V ( S ( t )) using (4.1), where the first integral in (4.3) is computedusing (4.4), and the last two integrals using (4.6), (4.7), and the existing valuesof v ± , v ± , y ( p − ) , y ( p +1 + p )Since the computational complexity of each step of the loop is O ( N log N ), thetotal complexity is O ( M N log N ). Remark 6.1.
While other quadrature methods typically employ multiple uniformgrids or specially-designed (nonuniform or shifted) grids, our method utilizes onlya fixed one-dimensional uniform grid, which not only eliminates the need for com-plicated inter-grid data transfer procedures, but also eliminates the need for specialsubroutines that are often required to interpolate data across discontinuities. Thismakes our method particularly easy to implement. Numerical Examples
We will demonstrate the accuracy and efficiency of the proposed method usingtwo examples, in which the value of an autocallable structured product and that ofa double barrier option with time-dependent barriers are found.7.1.
Example 1: Autocallable Structured Product.
We consider a knock-outautocallable structured product maturing in one year. The price of the underlyingasset is 3000, the nominal amount is 1, and the volatility is 20%. The observationdates (in years from now), barrier levels, and risk-free rates (in %) are given belowin Table 1.
Table 1.
An autocallable structured product.Observation date Barrier level Risk-free rate0.2 3050 20.4 3100 2.10.6 3150 2.20.8 3200 2.31 3250 2.4If the asset price reaches or goes above the barrier level at some observationdate t , the investor receives a payment of 4% × t . If the asset price is below thebarrier at every observation date, the investor will have to pay a premium of 1%.The relative errors of the computed option values with varying grid sizes are shown SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 15
500 1000 1500 2000 250000.20.40.60.811.21.4 x 10 −6 Grid Size R e l a t i v e E rr o r o f O p t i on V a l ue (a) −3 Number of Paths (Millions) R e l a t i v e E rr o r o f O p t i on V a l ue (b) Figure 1.
Errors for the autocallable structured product, com-puted using (a) the proposed method and (b) Monte-Carlo simu-lations.below in Figure 1(a), where the exact option value is taken to be the one computedon the grid of size 70,001.As a comparison, the relative errors of the option values computed using Monte-Carlo simulations with antithetic variates technique are shown in Figure 1(b). Asis clear from the figures, the error of the option value computed using the proposedmethod is well within 10 − with just 501 points in the grid, and drops very quicklyas the grid size increases. In contrast, it takes more than ten million paths forMonte-Carlo simulations to reduce the error of the computed option value to within10 − , and the error decays very slowly as the number of paths increases.Figure 2 below shows the CPU time used by the proposed method to price theautocallable structured product, where the code is developed in Matlab and is runon a personal computer. As is clear from the figure, the CPU time required by
500 1000 1500 2000 25002.533.544.55 x 10 −3 Grid Size C P U T i m e ( S e c ond s ) Figure 2.
CPU time for autocallable structured product valua-tion, using the proposed method.the proposed method is well within 0.01 seconds, and it increases approximatelylinearly as grid size increases. It is difficult to compare the speed of the proposed method with that of Monte-Carlo simulations directly, since the CPU time requiredby the latter depends largely on specific implementations. Nevertheless, the typicalCPU time consumed by a Monte-Carlo simulation with tens of millions of pathsranges from tens of seconds to a few minutes.7.2.
Example 2: Double Barrier Option.
As another example, consider aknock-out double barrier put option with time-dependent barrier levels. The priceof the underlying asset is 2500, the strike price is 2600, the nominal amount is 1,and the volatility is 25%. The option matures in two years. The observation dates(in years from now), barrier levels, and risk-free rates (in %) are given below inTable 2.
Table 2.
A double barrier option.Observation date Barrier level 1 Barrier level 2 Risk-free rate0.25 2200 2800 10.50 2100 2900 1.10.75 2000 3000 1.21 1900 3100 1.31.25 1800 3200 1.21.50 1700 3300 1.31.75 1600 3400 1.42 − −
700 900 1100 1300 1500 1700 1900 2100 2300 2500012345678 x 10 −6 Grid Size R e l a t i v e E rr o r o f O p t i on V a l ue (a) −4 Number of Paths (Millions) R e l a t i v e E rr o r o f O p t i on V a l ue (b) Figure 3.
Errors for the double barrier option, computed using(a) the proposed method and (b) Monte-Carlo simulations.
SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 17
As a comparison, the relative errors of the option values computed using Monte-Carlo simulations with antithetic variates technique are shown in Figure 3(b). Asis clear from the figures, the error of the option value computed using the proposedmethod is within 10 − with just 701 points in the grid, and drops very quickly asthe grid size increases. In contrast, it takes more than ten million paths for Monte-Carlo simulations to reduce the error of the computed option value to within 10 − ,and the error decays very slowly as the number of paths increases.Figure 4 below shows the CPU time used by the proposed method to pricethe double barrier option, where the code is developed in Matlab and is run ona personal computer. As is clear from the figure, the CPU time required by the
700 900 1100 1300 1500 1700 1900 2100 2300 25005.566.577.588.599.5 x 10 −3 Grid Size C P U T i m e ( S e c ond s ) Figure 4.
CPU time for double barrier option valuation, usingthe proposed method.proposed method is within 0.01 seconds, and it increases approximately linearlyas grid size increases. In contrast, the typical CPU time consumed by a Monte-Carlo simulation with tens of millions of paths ranges from tens of seconds to a fewminutes.
Remark 7.1.
Although the designed order of the proposed method is 4, in the abovenumerical examples, a lower order of convergence (close to 3) is actually observedfor the grid sizes considered. It is interesting to note that this apparent “loss” oforder of accuracy is not a defect of our method; rather, it is a manifestation of thesubtle influences that barrier levels can have on option pricing algorithms. Theseinfluences can be understood from two perspectives. First, in the above examples,the barrier levels K ± m are close to the spot price S ( t ) . This gives rise to a relativelysmall integration domain [ B − m , B + m ] compared with the entire computational domain [ − log C, log C ] (recall that B + m − B − m = log L + m S ( t ) − log L − m S ( t ) = log L + m L − m ≤ min n log K + m K − m , log K + m C − o ), which means that the set of grid points that are available for the discrete quad-rature rule (4.4) represents only a relatively small fraction of the set of grid pointsintroduced on the entire computational domain. Secondly, in the above examples theoption prices V m +1 ( S ) contain discontinuities at each observation date t m +1 . Thesediscontinuities necessarily show up in the form of large gradients in the (smooth)functions ˜ V m ( S ) (via the expectation integrals E (cid:2) V m +1 ( · ) | S (cid:3) ), which means that the discrete quadrature rule (4.4) is being applied to fast-varying functions with onlya relatively small number of grid points, leaving the integrands only marginally re-solved and hence explaining the degeneracy observed in the convergence rate. If thebarrier levels K ± m are pushed farther away from the spot price S ( t ) , so that theoption becomes increasingly like a vanilla option, then the discrete quadrature rule (4.4) effectively applies to slow-varying functions with a relatively large number ofgrid points, which improves the resolution of the integrands and hence the conver-gence rate (to close to 4). Despite these caveats on convergence rate, we emphasizethat our method is capable of pricing a sophisticated discretely monitored optionand obtaining five to six significant digits within a fraction of a second, while at thesame time being very easy to understand and to implement. Thus the (relativelytechnical) issue of convergence rate should pose no real concerns in practice. Remark 7.2.
Although a relative error of the order − or − may not alwaysseem necessary for option pricing problems considered in the real financial world,this extra accuracy is actually needed in the calculation of the Greeks, which are typ-ically approximated by finite difference formulas and which are much more sensitiveto numerical errors incurred in the calculation of option prices. Appendix
Estimate of Truncation Errors.
To estimate the truncation error for theintegral in (4.1), we first introduce
Lemma 8.1.
Let A = max (cid:8) | a ± | , . . . , | a ± M | , | a M | (cid:9) , B = max (cid:8) | b ± | , . . . , | b ± M | , | b M | (cid:9) ,R = min (cid:8) r , . . . , r M , (cid:9) , Q = min (cid:8) q , . . . , q M , (cid:9) . Then | ˜ V m ( S ) | ≤ e Q ( t m − t M ) AS + e R ( t m − t M ) B, ∀ ≤ m ≤ M, ∀ S ∈ ( K − m , K + m ) , and | V m ( S ) | ≤ e Q ( t m − t M ) AS + e R ( t m − t M ) B, ∀ ≤ m ≤ M, ∀ S ∈ (0 , ∞ ) . Proof.
Clearly, by assumption (cf. (3.4)), | ˜ V M ( S ) | = | V M ( S ) | ≤ AS + B, ∀ S ∈ (0 , ∞ ) . Now assume | ˜ V m ( S ) | ≤ e Q ( t m − t M ) AS + e R ( t m − t M ) B, ∀ S ∈ ( K − m , K + m ) , and | V m ( S ) | ≤ e Q ( t m − t M ) AS + e R ( t m − t M ) B, ∀ S ∈ (0 , ∞ ) , for some 1 ≤ m ≤ M . By (3.2) and Proposition 3.4, we have˜ V m − ( S ) = e − r m τ m /σ m Z ∞ V m ( y ) ρ m ( y, S ) dy, SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 19 for K − m − < S < K + m − . By definition (3.1), it is clear that ρ m ( y, S ) = ρ m ( y/S, /S .Thus a simple change of variable z = y/S yields | ˜ V m − ( S ) | = e − r m τ m /σ m (cid:12)(cid:12)(cid:12)Z ∞ V m ( Sz ) ρ m ( z, dz (cid:12)(cid:12)(cid:12) ≤ e − r m τ m /σ m Z ∞ (cid:2) e Q ( t m − t M ) ASz + e R ( t m − t M ) B (cid:3) ρ m ( z, dz = e − r m τ m /σ m (cid:18) e Q ( t m − t M ) AS Z ∞ zρ m ( z, dz + e R ( t m − t M ) B (cid:19) . The integral R ∞ zρ m ( z, dz is the expectation of the lognormal distribution, whichis simply exp { r m − q m ) τ m /σ m } [19]. Thus we have (observe R ≤ Q ≤ | ˜ V m − ( S ) | ≤ e − q m τ m /σ m + Q ( t m − t M ) AS + e − r m τ m /σ m + R ( t m − t M ) B ≤ e Q ( t m − − t M ) AS + e R ( t m − − t M ) B, for all S ∈ ( K − m − , K + m − ). By (3.3) and Proposition 3.4, we then deduce | V m − ( S ) | ≤ e Q ( t m − − t M ) AS + e R ( t m − − t M ) B, for all S ∈ (0 , ∞ ). The result then follows from induction. (cid:3) The possibly infinite integral in (3.5) is approximated by a finite integral. Tobe specific, let
C > , S = S ( t ), and ˜ G M − ( S ) = ˜ V M − ( S ). We consider ˜ G m − (1 ≤ m ≤ M −
1) defined recursively by˜ G m − ( S ) = e − r m τ m /σ m E (cid:2) ˜ G m ( · ) | S (cid:3) = e − r m τ m /σ m Z L + m L − m ˜ G m ( y ) ρ m ( y, S ) dy + a + m V am ( S, K + m , b + m V bm ( S, K + m ,
1) + a − m V am ( S, K − m , −
1) + b − m V bm ( S, K − m , − . A direct calculation shows that the errors ˜ R m = ˜ V m − ˜ G m satisfy the recursion˜ R m − ( S ) = e − r m τ m /σ m (cid:18)Z L + m L − m ˜ R m ( y ) ρ m ( y, S ) dy + Z L − m K − m ˜ V m ( y ) ρ m ( y, S ) dy + Z K + m L + m ˜ V m ( y ) ρ m ( y, S ) dy (cid:19) . We use the operator notation T m ( f )( S ) = e − r m τ m /σ m Z ∞ f ( y ) ρ m ( y, S ) dy, to write the recursion of ˜ R m as˜ R m − = T m ( ˜ R m χ ( L − m ,L + m ) ) + T m ( ˜ V m χ ( K − m ,L − m ] ∪ [ L + m ,K + m ) ) , for 1 ≤ m ≤ M −
1. Note also that ˜ R M − ( S ) = 0.Now we consider ˜ Q m defined by the recursion(8.1) ˜ Q m − = T m ( ˜ Q m ) + T m (cid:0) ( e Q ( t − t M ) Ay + e R ( t − t M ) B ) χ (0 ,S /C ] ∪ [ S C, ∞ ) (cid:1) , and ˜ Q M − = 0. It is easy to show using Lemma 8.1 and induction that | ˜ R m | ≤ ˜ Q m for all 0 ≤ m ≤ M −
1. We can also apply the recursion formula (8.1) to get theexpansion˜ Q = M − X m =1 T ◦ · · · ◦ T m (cid:0) ( e Q ( t − t M ) Ay + e R ( t − t M ) B ) χ (0 ,S /C ] ∪ [ S C, ∞ ) (cid:1) . Since T m is the risk-neutral expectation operator discounted from t m to t m − , T ◦ · · · ◦ T m ( f ) is simply the value of a European option whose payoff at t m is givenby f . Thus ˜ Q is equal to the value of a sum of 2( M −
1) binary call options withstrike S C and 2( M −
1) binary put options with strike S /C . As a result,˜ Q ( S ) = e Q ( t − t M ) A M − X m =1 (cid:2) C a ( S , S C, t m ) + P a ( S , S /C, t m ) (cid:3) (8.2) + e R ( t − t M ) B M − X m =1 (cid:2) C b ( S , S C, t m ) + P b ( S , S /C, t m ) (cid:3) , where C a , P a , C b , P b denote values of asset-or-nothing call, asset-or-nothing put,cash-or-nothing call, and cash-or-nothing put respectively. According to Lemma3.1, the values of these binary options are given by C a ( S , S C, t m ) = S e − ¯ q m ( t m − t ) N ( d ) , P a ( S , S /C, t m ) = S e − ¯ q m ( t m − t ) N ( d ) ,C b ( S , S C, t m ) = e − ¯ r m ( t m − t ) N ( d ) , P b ( S , S /C, t m ) = e − ¯ r m ( t m − t ) N ( d ) , where d = 1¯ σ m √ t m − t (cid:2) − log C + (¯ r m − ¯ q m + ¯ σ m )( t m − t ) (cid:3) ,d = 1¯ σ m √ t m − t (cid:2) − log C − (¯ r m − ¯ q m + ¯ σ m )( t m − t ) (cid:3) ,d = 1¯ σ m √ t m − t (cid:2) − log C + (¯ r m − ¯ q m − ¯ σ m )( t m − t ) (cid:3) ,d = 1¯ σ m √ t m − t (cid:2) − log C − (¯ r m − ¯ q m − ¯ σ m )( t m − t ) (cid:3) , and ¯ r m , ¯ q m , ¯ σ m represent the time-weighted averages of { r n , q n , σ n } mn =1 respectively.Generally, in practice, the absolute values of annual interest and yield rates will notexceed 50%, t M − t will not exceed 10 years, and A will not exceed 1. For generalautocallable structured products, B will not exceed t M − t , and for Bermudanoptions B will not exceed K , which is not much larger than S . We may also assume M ≤ σ = max ≤ m ≤ M σ m . If we chooselog C = 10 σ √ t M − t + (cid:0) σ (cid:1) ( t M − t ) , we can make sure that d , , , ≤ − , and thus N ( d , , , ) < − . SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 21
A crude estimate using (8.2) then shows that the error bound ˜ Q ( S ) does notexceed 10 − ( S + 1). This means the relative truncation error is negligible for allpractical purposes.8.2. Analysis of K ± m for Bermudan Options. The proper application of Propo-sition 3.4 requires the uniqueness of the exercise prices K ± m , which we now establishfor Bermudan options.To begin with, observe that the risk-neutral pricing formulas (3.2)–(3.3) appliedto Bermudan options can be written in an alternative form as(8.3) V m − ( S ) = max (cid:8) ˜ V m − ( S ) , ǫ ( S − K ) (cid:9) , ∀ ≤ m ≤ M, ∀ S ∈ (0 , ∞ ) , where ǫ = ( , if the option is a Bermudan call − , if the option is a Bermudan put , and ˜ V m − ( S ) = e − r m τ m /σ m Z ∞ V m ( y ) ρ m ( y, S ) dy (8.4) = e − r m τ m /σ m Z ∞ V m ( Sz ) ρ m ( z, dz. Since, by definition, V M ( S ) = max (cid:8) , ǫ ( S − K ) (cid:9) , (8.3) and (8.4) define V m and ˜ V m recursively for all 0 ≤ m ≤ M −
1. It is easy tosee that V m ( S ) ≥ ˜ V m ( S ) > ≤ m ≤ M − S ∈ (0 , ∞ ). Proposition 8.2.
Assume q m ≥ for all ≤ m ≤ M . For a Bermudan call optionwith strike K , the equation ˜ V m ( K + m ) = K + m − K has at most one finite solution,and < ˜ V m ( S ) − ˜ V m ( S ) < S − S , ∀ ≤ m ≤ M − , ∀ S > S > . For a Bermudan put option with strike K , the equation ˜ V m ( K − m ) = K − K − m has atmost one finite solution, and < ˜ V m ( S ) − ˜ V m ( S ) < S − S , ∀ ≤ m ≤ M − , ∀ S > S > . Proof.
Clearly Q = 0 since q m ≥
0. The proofs for call and put options are verysimilar, and we will present the argument for put options only which proceeds byinduction. Since ˜ V M − is the value of a European vanilla put option and q m ≥ − < ˜ V M − ( S ) − ˜ V M − ( S ) = − Z S S ˜ V ′ M − ( S ) dS < S − S , for all S > S >
0. Now suppose ˜ V M − ( S ) = K − S and ˜ V M − ( S ) = K − S for some S > S >
0. This means S − S = ˜ V M − ( S ) − ˜ V M − ( S ) < S − S , which is a contradiction. So the proposition is true for m = M − ≤ m ≤ M −
1, which implies that thefunction S − K + ˜ V m ( S ) is strictly increasing. Since ˜ V m ( S ) >
0, for sufficientlylarge S we have S − K + ˜ V m ( S ) >
0, or ˜ V m ( S ) > K − S . This means that K − m = inf (cid:8) S > V m ( S ) > K − S (cid:9) is well-defined and satisfies K − m < ∞ . Consider now (cf. (8.3))(8.5) V m ( S ) = ( ˜ V m ( S ) , S ≥ K − m K − S, S < K − m , and any S > S >
0. If S < K − m (and hence S < K − m ), we have V m ( S ) − V m ( S ) = ( K − S ) − ( K − S ) = S − S , by (8.5). If S ≥ K − m (and hence S ≥ K − m ), we have V m ( S ) − V m ( S ) = ˜ V m ( S ) − ˜ V m ( S ) ∈ (0 , S − S ) , by (8.5) and inductive hypothesis. If S < K − m ≤ S , we have V m ( S ) − V m ( S ) = ( K − S ) − ˜ V m ( S ) < ( K − S ) − ( K − S ) = S − S ,V m ( S ) − V m ( S ) > ( K − S ) − ˜ V m ( S ) ≥ , by (8.5), inductive hypothesis, and the definition of K − m . In conclusion, we haveshown that 0 < V m ( S ) − V m ( S ) = S − S , ∀ K − m > S > S > , (8.6a) 0 < V m ( S ) − V m ( S ) < S − S , otherwise . (8.6b)With the aid of (8.4) and (8.6), we write˜ V m − ( S ) − ˜ V m − ( S ) = e − r m τ m /σ m ( I + I ) , where 0 < I = Z K − m /S (cid:2) V m ( S z ) − V m ( S z ) (cid:3) ρ m ( z, dz = ( S − S ) Z K − m /S zρ m ( z, dz, < I = Z ∞ K − m /S (cid:2) V m ( S z ) − V m ( S z ) (cid:3) ρ m ( z, dz< ( S − S ) Z ∞ K − m /S zρ m ( z, dz. As a result, 0 < ˜ V m − ( S ) − ˜ V m − ( S )(8.7) < e − r m τ m /σ m ( S − S ) Z ∞ zρ m ( z, dz ≤ S − S , since by elementary properties of lognormal distributions [19], e − r m τ m /σ m Z ∞ zρ m ( z, dz = e − q m τ m /σ m ≤ . Now suppose ˜ V m − ( S ) = K − S and ˜ V m − ( S ) = K − S . This means˜ V m − ( S ) − ˜ V m − ( S ) = S − S , so by (8.7) we must have S = S . The proposition then follows from induction. (cid:3) SIMPLE AND EFFICIENT METHOD FOR OPTION PRICING 23
Corollary 8.3.
Assume a Bermudan put option with strike K has an optimal early-exercise level K − m > for some ≤ m ≤ M − . Then we have ˜ V m ( S ) > K − S for S > K − m and ˜ V m ( S ) < K − S for S < K − m . Similarly, for a Bermudan call optionwe have ˜ V m ( S ) < S − K for S > K + m and ˜ V m ( S ) > S − K for S < K + m .Proof. We will present the proof for put options only as the argument for calloptions is similar. It follows from Proposition 8.2 that the function ˜ V m ( S ) + S − K is increasing in S . Since ˜ V m ( K − m ) + K − m − K = 0, we have ˜ V m ( S ) > K − S for S > K − m and ˜ V m ( S ) < K − S for S < K − m . (cid:3) Acknowledgements
We are very grateful to Qingshuo Song for his valuable insights and helpfulsuggestions, as well as to Zhenan Sui for her careful reading and editing of themanuscript.
References [1] D.-H. Ahn, S. Figlewski, and B. Gao, Pricing discrete barrier options with an adaptive meshmodel, J. Deriv. 6(4), 33–43 (1999).[2] T. Alm, B. Harrach, D. Harrach, and M. Keller, A Monte Carlo pricing algorithm for auto-callables that allows for stable differentiation, J. Comput. Financ. 17(1), 43–70 (2013).[3] A. D. Andricopoulos, M. Widdicks, P. W. Duck, and D. P. Newton, Universal option valuationusing quadrature methods, J. Financ. Econ. 67(3), 447–471 (2003).[4] M. Broadie, P. Glasserman, and S. Kou, A continuity correction for discrete barrier options,Math. Financ. 7(4), 325–348 (1997).[5] M. Broadie, P. Glasserman, and S. G. Kou, Connecting discrete and continuous path-dependent options, Financ. Stoch. 3(1), 55–82 (1999).[6] M. Broadie and Y. Yamamoto, A double-exponential fast Gauss transform algorithm for pricingdiscrete path-dependent options, Oper. Res. 53(5), 764–779 (2005).[7] P. Buchen and O. Konstandatos, A new approach to pricing double-barrier options with arbi-trary payoffs and exponential boundaries, Appl. Math. Financ. 16(6), 497–515 (2009).[8] R. L. Burden, J. D. Faires, and A. M. Burden, Numerical Analysis (10th Edition), BrooksCole (2015).[9] G. Deng, J. Mallett, and C. McCann, Modeling autocallable structured products, J. Deriv.Hedge Funds 17(4), 326–340 (2011).[10] F. Fang and C. W. Oosterlee, Pricing early-exercise and discrete barrier options by Fourier-cosine series expansions, Numer. Math. 114, 27–62 (2009).[11] L. Feng and X. Lin, Pricing Bermudan options in L´evy process models, SIAM J. Finan. Math.4(1), 474–493 (2013).[12] L. Feng and V. Linetsky, Pricing discretely monitored barrier options and defaultable bondsin L´evy process models: a fast Hilbert transform approach, Math. Financ. 18(3), 337–384(2008).[13] C. P. Fries and M. S. Joshi, Perturbation stable conditional analytic Monte-Carlo pricingscheme for auto-callable products, Int. J. Theor. Appl. Finance 14(2), 197–219 (2011).[14] G. Fusai, I. D. Abrahams, and C. Sgarra, An exact analytical solution for discrete barrieroptions, Finance Stochast. 10(1), 1–26 (2006).[15] J. Gatheral, The Volatility Surface: A Practitioner’s Guide, John Wiley & Sons (2006).[16] A. Golbabai, L. V. Ballestra, and D. Ahmadian, A highly accurate finite element method toprice discrete double barrier options, Comput. Econ. 44(2), 153–173 (2014).[17] A. Ib´a˜nez and C. Velasco, The optimal method for pricing Bermudan options by simulation,Math. Financ. 28(4), 1143–1180 (2018).[18] J. C. Hull, Options, Futures, and Other Derivatives (9th Edition), Pearson (2014).[19] N. L. Johnson, S. Kotz, and N. Balakrishnan, “14: Lognormal Distributions”, ContinuousUnivariate Distributions (Volume 1, 2nd Edition), Wiley Series in Probability and Mathemat-ical Statistics, John Wiley & Sons (1994). [20] C. F. Lo, H. C. Lee, and C. H. Hui, A simple approach for pricing barrier options withtime-dependent parameters, Quant. Financ. 3(2), 98–107 (2003).[21] R. Lord, F. Fang, F. Bervoets, and C. W. Oosterlee, A fast and accurate FFT-based methodfor pricing early-exercise options under L´evy processes, SIAM J. Sci. Comput. 30(4), 1678–1705(2008).[22] C. O’Sullivan, Path dependant option pricing under L´evy processes, EFA 2005 Moscow Meet-ings Paper (2005).[23] E. Reiner, Convolution methods for path-dependent options, Financial Mathematics: RiskManagement, Modeling and Numerical Methods, IPAM UCLA (Jan. 3–12, 2001).[24] S. E. Shreve, Stochastic Calculus for Finance II: Continuous-Time Models, Springer-Verlag(2004).[25] T. Guillaume, Autocallable structured products, J. Deriv. 22(3), 73–94 (2015).[26] T. Guillaume, Analytical valuation of autocallable notes, Int. J. Financ. Eng. 2(2), 1–23(2015).[27] J. Z. Wei, Valuation of discrete barrier options by interpolations, J. Deriv. 6(1), 51–73 (1998).
Min Huang
China Merchants Bank7088 Shennan Boulevard, Shenzhen, Guangdong, ChinaEmail: [email protected]