A fast Fourier transform method for Mellin-type option pricing
aa r X i v : . [ q -f i n . P R ] M a r A fast Fourier transform method for Mellin-type option pricing
D.J. Manuge a, , P.T. Kim b a Global Risk Management, Scotiabank, 4 King Street West, Toronto, ON, Canada, M5H 1B6 b Department of Mathematics and Statistics, University of Guelph, 50 Stone Road, Guelph, ON, Canada, N1G 2W1
Abstract
Analytical pricing formulas and Greeks are obtained for European and American basket put options usingMellin transforms. We assume assets are driven by geometric Brownian motion which exhibit correlationand pay a continuous dividend rate. A novel approach to numerical Mellin inversion is achieved via thefast Fourier transform, enabling the computation of option values at equidistant log asset prices. Numericalaccuracy is verified among existing methods for American call options.
Keywords: basket option, American option, Mellin transform, fast Fourier transform, Black-Scholesformula, geometric Brownian motion
1. Introduction An option is a financial contract that presents its holder with the right, but not the obligation, to buy( call ) or sell ( put ) a given amount of asset at some future date. In practice, the underlying asset is oftenthe price of a stock, commodity, foreign exchange rate, financial index or futures contract. Although many styles of options exist, we are concerned with the valuation of European and
American varieties. Americanoptions may be exercised at any time t < T , while European options can only be exercised at time T . Inboth cases, their definitions can be extended to basket options, which differ by their dependence on n ∈ N underlying assets.Since the seminal paper of [2], much of the literature assumes assets are driven by geometric Brown-ian motion (GBM). Under this assumption, European option valuation relies on solving the Black-Scholespartial differential equation (PDE). With American options, the early-exercise condition gives rise to a freeboundary , in which no closed-form solution exists. The corresponding PDE is given by the inhomogeneousBlack-Scholes equation as in [6, 23, 27], where integral-based solutions are obtained. However, using theMellin transform to solve the PDE has only recently been considered. The novelty of the Mellin transformis threefold; one, the technique requires no change of variables or reduction to a diffusion equation; two,it enables option formulas to be expressed in terms of market asset prices, rather than logarithmic assetprices; and three, there exists a numerically fast scheme to compute multi-asset option prices. For pric-ing financial derivatives, the Mellin technique was first introduced in [12], where the authors consider theEuropean call option. Thereafter it was implemented in [34, 33], where the authors provide solutions forEuropean, American, and basket options on n = 2 underlying assets. For European options, weak pay-off functions [9], discrete dividends [8], transaction costs [30], and the Black-Scholes matrix equation [10]have since been considered in detail. However, in all of these cases continuous dividends are omitted. Thesingle-asset case with a continuous dividend is solved in [17] via an approach analogous to [34], in [16] via the discounted expectation formula for options, and [13, 35] via Mellin convolution. For American options,the single-asset case is solved in [17, 18] via an approach analogous to [34]. To the authors’ knowledge, thegeneral multi-asset formula for European and American basket options is not known. Disclaimer: The statements and conclusions of this paper do not necessarily reflect those of Scotiabank.
Email addresses: [email protected] (D.J. Manuge), [email protected] (P.T. Kim)
Preprint submitted to Elsevier October 4, 2018 or multi-asset options, scaling numerical procedures to higher dimensions can pose a challenge due tothe curse of dimensionality. Since assets are modelled by GBM, prices are log-normally distributed, but thesum of log-normal variables is not. In fact, the sum has no closed-form distribution function, making basketoption pricing a non-trivial task. By using the Mellin transform we are able to circumvent this by replacingthe distribution function with the characteristic function of the log price process. The majority of numericalpricing approaches for basket options rely on estimating analytical approximations via
Monte-Carlo methods[28]. However, this can be computationally expensive. Since basket options are ( n + 1)-dimensional in spaceand time, it is important to consider the complexity of the algorithm prior to computation and is ourmotivation for employing the fast Fourier transform (FFT).This paper makes three non-trivial contributions to the literature. First, we extend the existing Mellin-type pricing formulas for European, American, and more generally basket options to include n assets withcontinuous dividend rates. Second, we obtain new expressions for the Greeks of multi-asset European andAmerican options. Third, the American put option expression is discretized, yielding a new solution to thenumerical pricing problem. Computation of the solution relies on numerical Mellin inversion. In this paper,two methods to treat Mellin inversion are considered: a sine-cosine series expansion and a novel approachutilizing the FFT. Our FFT solution extends the European pricing method of [22].This work is organized as follows. In section 2, we derive the European put option formula on n assets.In section 3, we derive the American put option formula on n assets. In section 4, the American put optionformula for n assets is recast into a numerical procedure involving the FFT. This proposed solution alsoenables the explicit pricing of European put options. In section 5, we demonstrate the accuracy of theproposed method by comparing American call option prices computed using existing methods.
2. European Options
In this section, Mellin transforms are used to derive the formula for the price of a European basket putoption where assets have a continuous dividend rate and correlation. We begin by re-deriving the single-assetcase in [17], followed by the general multi-asset case.
For an option issued on a single asset, the value V = V ( S, t ; K ; T ; σ ; r ; q ) is dependent on an underlyingasset price 0 ≤ S ( t ) < ∞ , the exercise price K >
0, the maturity time 0 ≤ t ≤ T , the volatility (or standarddeviation) σ ≥ r ≥
0, and continuous dividend rate q ≥
0. TheBlack-Scholes equation for the price of a European option with dividend assets driven by geometric Brownianmotion is ∂V∂t + ( r − q ) S ∂V∂S + σ S ∂ V∂S − rV = 0 . (2.1)Equation (2.1) must satisfy the boundary conditions V ( S, T ) = θ ( S ) = max( K − S ) = ( K − S ) + and V ( S, t ) → S → ∞ . (2.2)Let M{ f ( x ); w } denote the Mellin transform of a function f ( x ) ∈ R + given by,ˆ f ( w ) := M{ f ( x ); w } = Z ∞ f ( x ) x w − dx (2.3)where complex variable w exists on an appropriate strip of convergence in C . Conversely, the inverse Mellintransform of a function ˆ f ( x ) ∈ C is defined by f ( x ) = M − { ˆ f ( w ); x } = 12 πi Z a + i ∞ a − i ∞ ˆ f ( w ) x − w dw (2.4)2here a ∈ ℜ ( w ), the real part of a ∈ C . Thus, to find the Mellin transform of the Black-Scholes equationapply (2.3) to equation (2.1): ∂ ˆ V ( w, t ) ∂t + (cid:18) σ w + w ) − ( r − q ) w − r (cid:19) ˆ V ( w, t ) = 0 . (2.5)By the final time condition (2.2), the general solution becomesˆ V ( w, t ) = ˆ θ ( w ) exp (cid:18) − σ (cid:0) w + (1 − k ) w − k (cid:1) ( T − t ) (cid:19) (2.6)where k = 2 r/σ , k = 2( r − q ) /σ , and ˆ θ ( w ) is the Mellin transform of the payoff function. Hence, byMellin inversion we obtain an expression for the price of a European put option on one asset, V PE ( S, t ) = 12 πi Z a + i ∞ a − i ∞ ˆ θ ( w ) e σ α ( w )( T − t ) S − w dw (2.7)where α ( w ) = w + (1 − k ) w − k and the Mellin transform of the put payoff function isˆ θ ( w ) = K w +1 w ( w + 1) (2.8)for ℜ ( w ) >
0. By setting q = 0, (2.7) reduces to (2 . .
11) in [34].
For an option issued on n assets, let S = ( S , ..., S n ) ′ , σ = ( σ , ..., σ n ) ′ and q = ( q , ..., q n ) ′ . The value V = V ( S , t ; K ; T ; σ ; r ; q ) is dependent on the underlying asset prices 0 ≤ S i ( t ) < ∞ , the exercise price K >
0, the maturity time 0 ≤ t ≤ T , the asset volatilities (or standard deviations) σ i ≥
0, the risk-freeinterest rate r ≥
0, and continuous dividend rates q i ≥ ∀ i . The assets are assumed to be driven bygeometric Brownian motion, dS i = µ i S i dt + σ i S i dW i (2.9)where the Wiener processes satisfy dW i ∼ Normal(0 , dt ) and corr( dW i , dW j ) = ρ ij for ρ ij ∈ [ − , µ i = r − q i − σ i X t , the charac-teristic function Φ( u ; t ) := exp[ − t Ψ( u )] = E [exp( i u ′ X t )] is given by the exponentΨ( u ) = 12 u ′ Σ u − i µ ′ u . (2.11)It is known under these conditions that the corresponding PDE for the price of a European basket optionis the generalized Black-Scholes equation: ∂V∂t + 12 n X i,j =1 ρ ij σ i σ j S i S j ∂ V∂S i ∂S j + n X i =1 ( r − q i ) S i ∂V∂S i − rV = 0 . (2.12)We note (2.12) must satisfy the boundary conditions V ( S , T ) = θ ( S ) = (cid:0) K − n X i =1 S i (cid:1) + and V ( S , t ) → S → ∞ . (2.13)3et M{ f ( x ); w } denote the multidimensional Mellin transform of a function f ( x ) ∈ R n + given by,ˆ f ( w ) := M{ f ( x ); w } = Z R n + f ( x ) x w − d x (2.14)where complex variable w = ( w , ..., w n ) ′ exists in an appropriate domain of convergence in C n . Conversely,the inverse multidimensional Mellin transform of a function ˆ f ( w ) ∈ C n is defined by f ( x ) = M − { ˆ f ( w ); x } = (2 πi ) − n Z γ ˆ f ( w ) x − w d w (2.15)where γ = n × j =1 γ j are strips in C n defined by γ j = { a j + ib j : a j ∈ R , b j = ±∞} with a j ∈ ℜ ( w j ). Thus, tofind the multidimensional Mellin transform of the generalized Black-Scholes equation apply (2.14) to (2.12): ∂ ˆ V∂t + 12 n X i,j =1 ρ ij σ i σ j w i w j ˆ V + 12 n X i =1 σ i w i ˆ V + ( r − q i ) n X i =1 w i ˆ V − r ˆ V = 0 . (2.16)By use of (2.10) and (2.11) we may rearrange the expression to obtain the ordinary differential equation d ˆ V ( w , t ) dt = (Ψ( w i ) + r ) ˆ V ( w , t ) . (2.17)Solving via the final time condition (2.13) yieldsˆ V ( w , t ) = ˆ θ ( w ) e − (Ψ( w i )+ r )( T − t ) . (2.18)Hence, by Mellin inversion we obtain our result. Theorem 1.
The Mellin-type formula for a European basket put option on n assets is given by V PE ( S , t ) = M − (cid:8) ˆ θ ( w )Φ( w i, T − t ) (cid:9) e − r ( T − t ) . (2.19) where Φ( ∗ ) is the characteristic function of a multivariate Brownian motion with drift and the Mellin trans-form of the payoff function is given by ˆ θ ( w ) = β n ( w ) K P w ( P w )( P w + 1) (2.20) for multinomial beta function β n ( w ) = Q nj =1 Γ( w j ) / Γ( P ni =1 w i ) , w ∈ C n , and ℜ ( w ) > J -dimensional Mellintransform of the put payoff function on J assets: Z R J + ( K − J X j =1 S i ) + J Y j =1 S w j − j dS j = Q Jj =1 Γ( w j )Γ(2 + P Jj =1 w j ) K P Jj =1 w j . (2.21)When J = 1 the expression equals (2.8) and thus holds. Assume J = n , then for J = n + 1 LHS = Z R ( n +1)+ ( K − n +1 X j =1 S i ) + n +1 Y j =1 S w j − j dS j = Q nj =1 Γ( w j )Γ(2 + P nj =1 w j ) Z K ( K − S n +1 ) P nj =1 w j S w n +1 − n +1 dS n +1 = Q n +1 j =1 Γ( w j )Γ(2 + P n +1 j =1 w j ) K P n +1 j =1 from Fubini’s theorem and (3.191.1) in [20]. The result follows from the definition of the multinomial betafunction and properties of gamma functions. 4 emark 1. An application of generalized put-call parity computes the price of a European call from a put(see [29]).
3. American Options
In this section, Mellin transforms are used to derive the formula for the price of an American basket putoption where assets have a continuous dividend rate and correlation. We begin by re-deriving the single-assetcase in [17], followed by the general multi-asset case.
As mentioned, the early exercise condition of American options produces a free boundary, which wedenote by the critical asset price S ∗ ( t ). For a put, when S > S ∗ (known as the continuation region ) itis optimal to hold the option, while when S < S ∗ (known as the exercise region ) it is optimal to exercisethe option. In order for the transition at the boundary to be smooth, the option and its gradient must becontinuous. The smooth pasting conditions supply this: ∂V ( S ∗ , t ) ∂S = − θ ( S ) = K − S ∗ . (3.1)The value V = V ( S, t ; K ; T ; σ ; r ; q ) of an American option on one asset is known to satisfy the inhomogeneousBlack-Scholes equation: ∂V∂t + ( r − q ) S ∂V∂S + σ S ∂ V∂S − rV = f (3.2)where the early exercise function is f ( S, t ) = ( − rK + qS, ≤ S ≤ S ∗ , S ∗ ( t ) ≤ S ≤ ∞ (3.3)and the final time condition is inherited from the European case: V ( S, T ) = θ ( S ) = ( K − S ) + . Furthermore,the boundary conditions imposed on (3.2) are V ( S, T ) = θ ( S ) = ( K − S ) + and V ( S, t ) → S → ∞ . (3.4)Similar to the European put case, the Mellin transform of (3.2) is given by ∂ ˆ V ( w, t ) ∂t + (cid:18) σ w + w ) − ( r − q ) w − r (cid:19) ˆ V ( w, t ) = ˆ f ( w, t ) . (3.5)The Mellin transform of the early exercise function isˆ f ( w, t ) = Z ∞ ( − rK + qS ) S w − dS = − rK S ∗ ( t ) w w + q S ∗ ( t ) w +1 w + 1 . (3.6)Solving (3.5) according to (3.11) and (3.17) yieldsˆ V ( w, t ) = ˆ θ ( w ) e σ α ( w )( T − t ) + Z Tt rKw S ∗ ( s ) w e σ α ( w )( s − t ) ds − Z Tt qw + 1 S ∗ ( s ) w +1 e σ α ( w )( s − t ) ds (3.7)where α ( w ) and ˆ θ ( w ) are defined in section 2.1. By Mellin inversion we obtain the price of an Americanoption on asset driven by geometric Brownian motion: V PA ( S, t ) = 12 πi Z c + i ∞ c − i ∞ ˆ θ ( w ) e σ α ( w )( T − t ) S − w dw + 12 πi Z a + i ∞ a − i ∞ Z Tt rKw (cid:18) S ( t ) S ∗ ( s ) (cid:19) − w e σ α ( w )( s − t ) dsdw − πi Z a + i ∞ a − i ∞ Z Tt qS ∗ ( s ) w + 1 (cid:18) S ( t ) S ∗ ( s ) (cid:19) − w e σ α ( w )( s − t ) dsdw. (3.8)5he first term is the European option formula (2.7), while the second and third term represent the contri-bution of the early exercise premium. Note that when q = 0, we obtain (3.1.9) in [34]. For multiple assets, the continuation region exists for P ni =1 S i > S ∗ , while the exercise region exists for P ni =1 S i < S ∗ . The value V = V ( S, t ; K ; T ; σ ; r ; q ) of an American option on one asset is known to satisfythe inhomogeneous generalized Black-Scholes equation: ∂V∂t + 12 n X i,j =1 ρ ij σ i σ j S i S j ∂ V∂S i ∂S j + n X i =1 ( r − q i ) S i ∂V∂S i − rV = f (3.9)where the early exercise function is f ( S , t ) = ( − rK + P ni =1 q i S i , < P ni =1 S i ≤ S ∗ ( t )0 , S ∗ ( t ) < P ni =1 S i < ∞ . (3.10)Similar to the European case, the boundary conditions imposed on (3.9) are V ( S , T ) = θ ( S ) = (cid:0) K − n X i =1 S i (cid:1) + and V ( S , t ) → S → ∞ . (3.11)The smooth pasting conditions along the boundary are ∂V ( S , t ) ∂S i (cid:12)(cid:12)(cid:12)(cid:12) P ni =1 S i = S ∗ = − θ ( S ) = K − S ∗ . (3.12)The multidimensional Mellin transform of (3.9) is given by the expression ∂ ˆ V∂t + 12 n X i,j =1 ρ ij σ i σ j w i w j ˆ V + 12 n X i =1 σ i w i ˆ V + ( r − q i ) n X i =1 w i ˆ V − r ˆ V = ˆ f . (3.13)By use of (2.10) and (2.11) we may rearrange (3.13) to obtain the ordinary differential equation d ˆ V ( w , t ) dt − (Ψ( w i ) + r ) ˆ V ( w , t ) = ˆ f ( w , t ) . (3.14)Solving via the final time condition (3.12) and applying Duhamel’s principle yieldsˆ V ( w , t ) = ˆ θ ( w ) e − (Ψ( w i )+ r )( T − t ) − T Z t ˆ f ( w , s ) e − (Ψ( w i )+ r )( s − t ) ds. (3.15)Hence, by Mellin inversion we obtain our result. Theorem 2.
The Mellin-type formula for an American basket put option on n assets is given by V PA ( S , t ) = e − r ( T − t ) M − n ˆ θ ( w )Φ( w i, T − t ) o − M − n Z Tt ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) ds o (3.16) where Φ( ∗ ) is the characteristic function of a multivariate Brownian motion with drift, ˆ θ ( ∗ ) is the Mellintransform of the payoff function given by (2.20) , and the Mellin transform of the early exercise function isgiven by ˆ f ( w , t ) = β n ( w )( S ∗ ) P w P w (cid:20) q ′ w S ∗ P w + 1 − rK (cid:21) (3.17) for critical asset price S ∗ ( t ) , multinomial beta function β n ( w ) = Q nj =1 Γ( w j ) / Γ( P ni =1 w i ) , w ∈ C n , and ℜ ( w ) > . J -dimensional Mellintransform of the early exercise function on J assets: Z R J + (cid:18) − rK + J X i =1 q i S i (cid:19) J Y j =1 S w j − j dS j = Q Jj =1 Γ( w j )( S ∗ ) P Jj =1 w j Γ(1 + P Jj =1 w j ) (cid:20) S ∗ P Jj =1 q j w j P Jj =1 w j + 1 − rK (cid:21) . When J = 1 the expression equals (3.17) and thus holds. Assume J = n , then for J = n + 1 LHS = Z R ( n +1)+ (cid:18) − rK + n +1 X i =1 q i S i (cid:19) n +1 Y j =1 S w j − j dS j = − rK Q nj =1 Γ( w j )Γ(1 + P nj =1 w j ) Z S ∗ ( S ∗ − S n +1 ) P nj =1 w j S w n +1 − n +1 dS n +1 + P n +1 j =1 q j w j Q nj =1 Γ( w j )Γ(2 + P nj =1 w j ) Z S ∗ ( S ∗ − S n +1 ) P nj =1 w j S w n +1 − n +1 dS n +1 = − rK Q n +1 j =1 Γ( w j )Γ(1 + P n +1 j =1 w j ) ( S ∗ ) P n +1 j =1 + P n +1 j =1 q j w j Q n +1 j =1 Γ( w j )Γ(2 + P n +1 j =1 w j ) ( S ∗ ) P n +1 j =1 from Fubini’s theorem and equation (3.191.1) in [20]. The result follows from the definition of the multinomialbeta function and properties of gamma functions. Remark 2.
An application of generalized put-call symmetry gives the price of an American call option froma put (see [32]).
Note that the early exercise premium only contributes to the price of the option when P ni =1 S i ( s ) ≤ S ∗ ( s ).Otherwise the second term of (3.16) is zero. By imposing the smooth pasting conditions (3.12) on (3.16),we obtain an implicit equation describing the free boundary. Corollary 1.
The critical asset price S ∗ ( t ) is given by the solution of the expression K − S ∗ ( t ) = e − r ( T − t ) πi Z γ ˆ θ ( w )Φ( w i, T − t ) S ∗ ( t ) − w d w − Z γ Z Tt ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) S ∗ ( t ) − w dsd w . (3.18)The critical asset price can be obtained by solving for S ∗ ( t ) where S ∗ ( t ) = ( S ∗ , ...S ∗ n ) over the space ofpossible prices in R n + such that S ∗ = P ni =1 S ∗ i .
4. Option Sensitivities
Option sensitivities or Greeks describe the relationship between the value of an option and changes inone of its underlying parameters. They play a vital role for risk management and portfolio optimization,since they have the ability to describe how vulnerable an option is to a particular risk factor. They areeasily obtained for European and American options by passing the appropriate derivative operator under thecomplex integral in (3.16). For succinctness, the variable change τ = T − t is used in some of the followingexpressions. The first partial derivative with respect to a given asset, Delta, is given by∆ := ∂V∂S i = − e − rτ M − n w i S i ˆ θ ( w )Φ( w i, τ ) o + M − n w i S i Z τ ˆ f ( w , τ − s )Φ( w i, s ) e − rs ds o . (4.1)The cross partial derivative with respect to two independent assets is given by∆ := ∂ V∂S i ∂S j = − e − rτ M − n w i S i w j S j ˆ θ ( w )Φ( w i, τ ) o + M − n w i S i w j S j Z τ ˆ f ( w , τ − s )Φ( w i, s ) e − rs ds o . (4.2)7amma, the second derivative with respect to the asset price is given byΓ := ∂ V∂S i = − e − rτ M − n w i (1 − w i )ˆ θ ( w )Φ( w i, τ ) S − i o − M − n S − Z Tt w i (1 − w i ) ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) ds o . (4.3)Theta, the first partial derivative with respect to time isΘ := − ∂V∂t = − e − r ( T − t ) M − n (Ψ( w i ) + r )ˆ θ ( w )Φ( w i, T − t ) o + M − n Z Tt (Ψ( w i ) + r −
1) ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) ds o . (4.4)Rho, the first partial derivative with respect to the risk-free rate of return is given by ρ := ∂V∂r = − τ e − rτ M − n ( n X j =1 w i − T − t )ˆ θ ( w )Φ( w i, τ ) o − M − n Z Tt ( n X j =1 w i − s − t ) ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) ds o . (4.5)Nu, the first partial derivative with respect to volatility is given by ν := ∂V∂σ i = τ e − rτ M − n(cid:2) n X i,j =1 i = j ρ ij σ j w i w j + n X i =1 σ i w i ( w i − (cid:3) ˆ θ ( w )Φ( w i, τ ) o − M − n(cid:2) n X i,j =1 i = j ρ ij σ j w i w j + n X i =1 σ i w i ( w i − (cid:3) Z τ s ˆ f ( w , τ − s )Φ( w i, s ) e − rs ds o . (4.6)Finally, the first partial derivative with respect to the dividend rate is given byΞ := ∂V∂q i = − τ e − rτ M − n w i ˆ θ ( w )Φ( w i, τ ) o + M − n Z Tt w i ( s − t ) ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) ds o . (4.7)By eliminating the second term for each Greek we obtain the corresponding European option sensitivities.Since most payoff functions are independent of the derivative operator, these expressions also hold for manypath-independent multi-asset options. The American case differs because the exercise region varies with timeand depends on the payoff function. Even in the simplest case of the basket option, the Mellin transformof the early exercise function is dependent on the derivative operator and must be considered to obtainexpressions for other multi-asset Greeks. Remark 3.
By direct substitution of the above expressions, we may prove that ( i ) formula (2.19) is a classicalsolution to the European pricing problem (2.12) - (2.13) and ( ii ) formula (3.16) is a classical solution to theAmerican pricing problem (3.9) - (3.12) . . Numerical Solution using the Fast Fourier Transform Valuing options on n underlying assets is a difficult problem due to the curse of dimensionality. Theissues stem from multiple integration, where the order of complexity does not scale linearly as n increases.The fast Fourier transform (FFT), a numerically efficient discrete Fourier transform, is able to circumventthis problem by reducing the number of floating point operations from O ( N n ) to O ( N n log N n ) (when thenumber of transformed points N are equal across dimensions). In this section we present a new FFT-basedmethod that enables the pricing of both European and American basket options. Recall (3.16), for strip ofconvergence γ = a + i b where b → ±∞ , V ( S , t ) = e − r ( T − t ) (2 πi ) n lim b →∞ Z a + i ba − i b ˆ θ ( w )Φ( w i, T − t ) S − w d w − (2 πi ) − n lim b →∞ Z a + i ba − i b Z Tt ˆ f ( w , s )Φ( w i, s − t ) e − r ( s − t ) S − w dsd w . (5.1)Make a change of variables by setting w = a + i b so that d w = id b . Then, V ( S , t ) = e − r ( T − t ) (2 π ) n lim b →∞ Z b − b ˆ θ ( a + i b )Φ( a i − b , T − t ) S − ( a + i b ) d b − (2 π ) − n lim b →∞ Z b − b Z Tt ˆ f ( a + i b , s )Φ( a i − b , s − t ) e − r ( s − t ) S − ( a + i b ) dsd b . (5.2)Induce the time change τ = T − t and discretize the integrals over b and s by invoking the Trapezoid rule. V ( S , τ ) ≃ ∆ b e − rτ (2 π ) n N − X j ,...,j n =0 ˆ θ ( a + i b j )Φ( a i − b j , τ ) e − ( a + i b j ) ′ ln( S ) − ∆ b ∆ τ (2 π ) n N − X j ,...,j n =0 M − X l =0 ˆ f ( a + i b j , τ − t l )Φ( a i − b j , t l ) e − rt l − ( a + i b j ) ′ ln( S ) . (5.3)Time is parameterized by t l := lτ / ( M −
1) for stepsize ∆ τ = τ /M and vector l = 0 , ..., M −
1. Similarly,the Mellin integrals are defined by b j = ( b j , ..., b j n ) where b j i := ( j i − N )∆ i for j i = 0 , ..., N −
1, and∆ b = Q ni =1 ∆ i . Hence, the multiple integral in b is approximated by a multiple sum over the lattice, B = { b j = ( b j , ..., b j n ) | j = ( j , ..., j n ) ∈ { , ..., N − } n } . To evaluate the price inputs, define the initial log asset prices by the reciprocal lattice S = { s k = ( s k , ..., s k n ) | k = ( k , ..., k n ) ∈ { , ..., N − } n } where s k := ( k i − N ) λ i for k i = 0 , ..., N −
1. The well-known European FFT procedures of [7, 14] mostnoticeably differ from our approach by using log exercise prices rather than log asset prices for the FFTgrid. The idea of using log-asset prices comes from [22]. Although the integral extension from Europeanto American options is quite natural, existing FFT-based algorithms do not rely on approximating integralsolutions. Rather, there are two main approaches. One approach is the FFT convolution method forBermudan options in [31]. Bermudan options are able to provide an approximation to American options whenthe number of early exercise points reach infinity. The other approach is based on the linear complementarityformulation of American options. Coined the Fourier time-stepping method, the method relies on enforcingthe condition V ( S , t ) ≥ V ( S , T ) at each timestep over the lifetime of the option [24]. While both of thesemethods sufficiently price American options, the decomposition derived in (5.3) allow us to price bothEuropean and American options by considering a single formula. Further simplification can be made by9ecognizing that each sum in j i is truncated to N evaluation points. By setting ∆ i λ i = 2 π/N one obtains V ( S , τ ) ≃ ( − P k ∆ b e − rτ (2 π ) n N − X j ,...,j n =0 ζ E e − a ′ s e − πiN j ′ k − ( − P k ∆ b ∆ τ (2 π ) n N − X j ,...,j n =0 M − X l =0 ζ EEP e − rt l − a ′ s e − πiN j ′ k (5.4)where ζ E ( b j ) = ( − P j ˆ θ ( a + i b j )Φ( a i − b j , τ ) e − rτ (5.5)and ζ EEP ( b j , t l ) = ( − P j ˆ f ( a + i b j , τ − t l )Φ( a i − b j , t l ) e − rt l . (5.6)Under a change of variables u = wi , (5.4) is equivalent to the method of [22] when solving for Europeanoptions of unit exercise price. To obtain American options, two FFT procedures must be computed withinput arrays ζ E ( b j ) and ζ EEP ( b j , t l ). Alternatively, by combining the integrands we need only compute oneFFT, thus reducing the speed of the algorithm. An improvement in accuracy can be made by introducingthe composite Simpson’s rule over k and j . This allows the integrand to be approximated using quadraticpolynomials rather than line segments. By defining α = (3 + ( − P j − δ P j ) /
3, this weighted smoothingimplies V PA ( S , τ ) ≃ ( − P k ∆ b (2 π ) n FFT (cid:8) αζ E ( b j ) (cid:9) e − a ′ s − ( − P k ∆ b ∆ τ (2 π ) n FFT (cid:26) M − X l =0 αζ EEP ( b j , t l ) (cid:27) e − a ′ s (5.7)= ( − P k ∆ b (2 π ) n FFT (cid:26) αζ E ( b j ) − α ∆ τ M − X l =0 ζ EEP ( b j , t l ) (cid:27) e − a ′ s (5.8)where the Kronecker delta function δ P j = 1 for P j = 0 and zero otherwise. The first term of (5.7)computes the price of a European put option, while the second corresponds to the early exercise premium.Evaluating both terms, or equivalently (5.8), retrieves the value of an American put option. The error ofthe numerical procedure will depend highly on the choice of N , M , ∆ i (or λ i ), and a . Careful selectionmust be made with ∆ i (or λ i ) for the reciprical FFT grid to land on the initial log asset price specified atinput. One way this can be achieved is by solving for the root of f ( λ i ) = ln( S i ) − ( k i − N/ λ i , satisfying0 ≤ λ i ≤ ≤ k i ≤ N −
1. Since the size of the grid step shrinks as k i increases, λ i must be largeenough so that the log price is contained by the range of the FFT grid, yet small enough to obtain a finegrid between log prices. A fine grid may also be achieved by increasing the number of evaluation points N .Further details on parameter selection and computational error with FFT-based pricing are given in [22].
6. Application: Pricing American Call Options
In this section we explicitly compute American call options using (5.7) and applying put-call symmetry: V CA ( S, K, r, q, t ) = V PA ( K, S, q, r, t ). Numerical methods are coded in R . Experiments are run on a Windows7 OS machine in R Studio with Intel Core i3 CPU @ 2.53 GHz and 4 GB RAM. Prior to computing (5.7), thecritical asset price S ∗ must be determined. A typical procedure is to recursively solve (3.18) for S ∗ at eachtimestep t l ∈ [0 , T ]. If the parameters { K, r, q, σ, T } of an American option are known prior to pricing, thecritical asset price can be calculated ex-ante and stored to reduce runtime. In practice, parameters such asvolatility and time to maturity continuously change. Hence, one may wish to reduce runtime by computingan analytical approximation; often posed as an implicit function of S ∗ . We consider proposition 5.3.3. of[17, 18] for a dividend-paying asset: S ∗ ( t ) = K rσ √ δ N ( p δ ( T − t )) − e q ( T − t ) [ N ( κ ) − N ( p ( δ − q )( T − t ))] + ω + (6.1)10here ω = 2 q + σ √ δ − q σ √ δ [2 N ( p δ ( T − t )) −
1] (6.2) δ = σ q − rσ + 2 r (6.3) κ = ln( S ∗ /K ) + ( r − q + σ / T − t ) σ √ T − t (6.4)and N ( · ) is the cumulative Normal distribution. For our proposed experiments, prices are computed usingequation (6.1) for the critical asset price. The implicit function is numerically solved by Brent’s method.In table 1 existing methods for computing 6-month American call option prices are compared againstthe benchmark binomial options formula with 10000 timesteps (True) in [11]. Including the proposed FFTmethod of section 5 (FFT), the following numerical procedures are considered: the method of Barone-Adesi and Whaley (BAW) [1], the four-point method of Geske and Johnson (GJ4) [19], the modified two-point Geske-Johnson approach of Bunch and Johnson (BJ2) [4], the four-point schemes of Huang et al.(HSY4) [21], the lower and upper bound approximation of Broadie and Detemple (LUBA) [3], the four-pointrandomization method of Carr (RAN4) [5], the three-point multi-piece exponential boundary approximationof Ju (EXP3) [25], an approximation of Ju and Zhong (JZ) [26], the Gauss-Laguerre quadrature method ofFrontczak and Schobel (GL) [18], and the Mellin-inversion scheme of Dishon and Weiss (DW) (see AppendixA). S True FFT DW BAW GJ4 BJ2 HSY4 LUBA RAN4 EXP3 JZ GL
80 0.2194 0.2198 0.2198 0.2300 0.2191 0.2186 0.2199 0.2195 0.2188 0.2196 0.2216 0.218590 1.3864 1.3894 1.3895 1.4050 1.3849 1.3818 1.3898 1.3862 1.3802 1.3872 1.3857 1.3851100 4.7825 4.7942 4.7943 4.7821 4.7851 4.7862 4.8044 4.7821 4.7728 4.7837 4.7682 4.7835110 11.0978 11.1269 11.1270 11.0409 11.0889 11.2553 11.0686 11.0976 11.0893 11.0993 11.0794 11.1120120 20.0004 20.0594 20.0591 20.0000 20.0073 20.0000 20.0531 20.0000 20.0000 20.0005 20.0000 20.000080 2.6889 2.6921 2.6921 2.7108 2.6864 2.6827 2.6897 2.6893 2.6787 2.6899 2.6871 2.678890 5.7223 5.7298 5.7297 5.7416 5.7212 5.7163 5.7361 5.7231 5.7113 5.7237 5.7110 5.7195100 10.2385 10.2539 10.2538 10.2417 10.2451 10.2351 10.2752 10.2402 10.2205 10.2404 10.2143 10.2265110 16.1812 16.2076 16.2074 16.1520 16.1831 16.2107 16.2012 16.1817 16.1629 16.1831 16.1456 16.1756120 23.3598 23.4013 23.4010 23.2883 23.3419 23.4771 23.3288 23.3574 23.3389 23.3622 23.3211 23.382880 1.6644 1.6643 1.6644 1.6645 1.6644 1.6644 1.6644 1.6644 1.6604 1.6644 1.6644 1.664490 4.4947 4.4946 4.4947 4.4950 4.4946 4.4947 4.4947 4.4947 4.4959 4.4947 4.4947 4.4947100 9.2504 9.2505 9.2506 9.2513 9.25091 9.2506 9.2506 9.2506 9.2513 9.2506 9.2507 9.2506110 15.7977 15.7974 15.7975 15.7988 15.7973 15.7975 15.7975 15.7975 15.7994 15.7975 15.7977 15.7980120 23.7061 23.706 23.7062 23.07086 23.7082 23.7062 23.7062 23.7062 23.7027 23.7062 23.7066 23.7060
Table 1: American call option prices calculated using twelve different pricing methods at varying risk-free rate r , dividendrate q , and volatility σ . All options have a 6-month expiry and are calculated with exercise price K = 100 for asset prices S = { , , , , } . The first grouping is calculated with r = 0 . q = 0 .
07, and σ = 0 .
2. The second grouping iscalculated with r = 0 . q = 0 .
07, and σ = 0 .
4. The third grouping is calculated with r = 0 . q = 0 .
03, and σ = 0 . Parameter selection for each method coincides with the original references, excluding the Mellin-basedFFT and DW methods which were introduced here. As previously mentioned, the spacing for the FFT grid∆ b must be chosen a priori for the panel of option prices to land on the appropriate asset price. Since we arepricing call options by pull-call symmetry, our grid spacing will depend on K instead of S . By fixing N = 2 evaluation points and K = 100 across all experiments, solving for the root of f ( λ ) = ln( K ) − ( k − N/ λ yields a grid spacing of ∆ b = 0 . ℜ ( w ) = a >
0, itmust evaluated at a given point. The integrand of (3.8) tends to oscillate as a approaches the endpointson (0 , ∞ ). For this reason, the arbitrary selection of a = 1 is made. In addition, M = 250 timesteps ischosen to evaluate the trapezoid rule in the Mellin transform of the early exercise function. For the DWexperiments we adopt N = 250 evaluation points, a = 1 for the strip of convergence, M = 250 timesteps forthe Mellin transform of the early exercise function, and L = 10 for the bounds on the log-price range.Note that we may alternatively obtain the American option price by directly computing equation (5.3).In this case, pricing error is comparable to the FFT method; the benefit of the FFT method stems from its11peed, not necessarily an improvement in accuracy. Computationally, the FFT method most notably differsby generating 2 N option prices, while the trapezoid rule generates one. It is often the case that one wishes todetermine a single option price, however computing a panel of option prices may be viable when the initialasset price is unknown. For example, suppose a stock option is issued at some future date. By forecastingan expected price range for the stock on the date of issuance, one can determine the corresponding pricerange for the option. This eliminates having to compute multiple valuations at different forecasted assetprices. Even if one option price is required, the FFT algorithm along with a simple index search returns therequired option price in less runtime than computing the equivalent trapezoid rule. For example, using thesame parameters as in table 1 the FFT method takes ∼ ∼ { K, r, q, σ, τ } prior to pricing. Or as mentioned,if the tradeoff in error is warranted, one may compute an exact form for the critical asset price. Otheranalytical approximations may also be explored. Although we concern ourselves with American optionpricing, European options are easily obtained by computing the first term in (5.7). Since no free boundaryexists in this case, there will be less error in the option price. For example, using the same parameterchoices as row 1 of table 1, the absolute pricing error is on the order of 10 − when compared against theBlack-Scholes formula. We should note that the small error may be the result of precision in R and is wellwithin tolerances required by practitioners. Although omitted from this manuscript, numerically pricinghigher dimensional European, American, and exotic options is feasible by the proposed method.
7. Conclusion
In the context of Mellin transforms, we obtain analytic solutions for the fair value of basket put optionsand Greeks on n assets with continuous dividend rates and correlation. Solutions are obtained for bothEuropean and American option styles. By expanding on the European framework of [22], we obtain anumerical solution to the American basket put option via the fast Fourier transform. The decompositionof the solution enables the direct computation of either European or American basket option prices. Bysolving for the Mellin transform of alternate payoff functions, the results presented here may be used to pricemore complicated multi-asset options. Numerical results are compared against twelve methods for pricingAmerican call options, including two additional approaches to treat the Mellin inversion in our main result.The results verify the efficiency and accuracy of the proposed numerical solution. Acknowledgments
This research was supported in part by the Natural Sciences and Engineering Research Council of Canada,Grant DG 46204. The first author would like to thank participants of the 4th New York Conference onApplied Mathematics at Cornell University where this research was presented.
References [1] Barone-Adesi, G., Whaley, R. E., 1987. Efficient analytic approximation of American option values. Journal of Finance42 (2), 301–20.[2] Black, F., Scholes, M., 1973. The pricing of options and corporate liabilities. Journal of Political Economy 81 (3), 637.[3] Broadie, M., Detemple, J. B., 1994. American option valuation: New bounds, approximations, and a comparison of existingmethods. Cirano working papers, CIRANO.[4] Bunch, D. S., Johnson, H., 1992. A simple and numerically efficient valuation method for American puts using a modifiedGeske-Johnson approach. Journal of Finance 47 (2), 809–16.[5] Carr, P., Jul. 1998. Randomization and the American put. Review of Financial Studies 11 (3), 597–626.[6] Carr, P., Jarrow, R., Myneni, R., 1992. Alternative characterizations of American put options. Mathematical Finance2 (2), 87–106.
7] Carr, P. P., Madan, D. B., 1999. Option valuation using the fast Fourier transform. Journal of Computational Finance 2,61–73.[8] Company, R., Gonz´alez, A. L., J´odar, L., 2006. Numerical solution of modified Black-Scholes equation pricing stock optionswith discrete dividend. Mathematical and Computer Modelling 44 (11-12), 1058–1068, eng.[9] Company, R., J´odar, L., Rubio, G., Villanueva, R.-J., 2007. Explicit solution of Black-Scholes option pricing mathematicalmodels with an impulsive payoff function. Mathematical and Computer Modelling 45 (1-2), 80–92.[10] Cort´es, J. C., J´odar, L., Sala, R., Sevilla-Peris, P., 2005. Exact and numerical solution of Black-Scholes matrix equation.Applied Mathematics and Computation 160 (3), 607–613.[11] Cox, J. C., Ross, S. A., Rubinstein, M., 1979. Option pricing: A simplified approach. Journal of Financial Economics7 (3), 229–263.[12] Cruz-B´aez, D., Gonz´alez-Rodr´ıguez, J., 2002. Semigroup theory applied to options. J. Appl. Math. 2 (3), 131–139.[13] Cruz-B´aez, D. I., Gonz´alez-Rodr´ıguez, J. M., 2005. A different approach for pricing European options. In: Proceedings ofthe 8th WSEAS International Conference on Applied Mathematics. MATH’05. World Scientific and Engineering Academyand Society (WSEAS), Stevens Point, Wisconsin, USA, pp. 373–378.[14] Dempster, M. A. H., Hong, S. G., 2002. Spread option valuation and the fast Fourier transform. In: Mathematical Finance- Bachelier Congress 2000. Springer Berlin Heidelberg, pp. 203–220.[15] Dishon, M., Weiss, G. H., 1978. Numerical inversion of Mellin and two-sided Laplace transforms. Journal of ComputationalPhysics 28 (1), 129 – 132.[16] Dufresne, D., Garrido, J., Morales, M., 2009. Fourier inversion formulas in option pricing and insurance. Methodology andComputing in Applied Probability 11 (3), 359–383.[17] Frontczak, R., 2010. On the application of Mellin transforms in the theory of option pricing. Ph.D. thesis, T¨ubingenUniversity.[18] Frontczak, R., Sch¨obel, R., 2010. On modified Mellin transforms, Gauss-Laguerre quadrature, and the valuation of Amer-ican call options. J. Computational Applied Mathematics 234 (5), 1559–1571.[19] Geske, R., Johnson, H. E., December 1984. The American put option valued analytically. Journal of Finance 39 (5),1511–24.[20] Gradshteyn, I. S., Ryzhik, I. M., 2000. Table of Integrals, Series and Products, 6th Edition. Academic Press, San Diego.[21] Huang, J., Subrahmanyam, M. G., Stern, L. N., Yu, G. G., Feinstein, S., Figlewski, S., Gao, B., Jensen, B., Johnson, H.,We, N. W., 1996. Pricing and hedging American options: A recursive integration method. Review of Financial Studies 9,277–300.[22] Hurd, T. R., Zhou, Z., 2010. A Fourier transform method for spread option pricing. SIAM J. Financial Math. 1 (1),142–157.[23] Jacka, S. D., 1991. Optimal stopping and the American put. Mathematical Finance 1 (2), 1–14.[24] Jackson, K. R., Jaimungal, S., Surkov, V., Building, S. F., 2008. Fourier space time-stepping for option pricing with L´evymodels. Journal of Computational Finance 12 (2), 1–29.[25] Ju, N., 1998. Pricing an american option by approximating its early exercise boundary as a multipiece exponential function.Review of Financial Studies 11 (3), 627–46.[26] Ju, N., Zhong, R., 1999. An approximate formula for pricing American options. Journal of Derivatives 7 (2), 31–40.[27] Kim, I. J., 1990. The analytic valuation of American options. Review of Financial Studies 3 (4), 547–72.[28] Krekel, Martin. de Kock, J. K. R. M. T.-K., 2004. An analysis of pricing methods for basket options. Wilmott Magazine (3),82–89.[29] Laurence, P., Wang, T.-H., 2005. Sharp upper and lower bounds for basket options. Applied Mathematical Finance 12 (3),253–282.[30] Lin Cheng, F., 2011. Mellin transform solution for the model of European option. In: EMEIT. IEEE, pp. 329–331.[31] Lord, R., Fang, F., Bervoets, F., Oosterlee, C. W., February 2008. A Fast And Accurate FFT-Based Method For PricingEarly-Exercise Options Under Levy Processes. SIAM Journal on Scientific Computing 30, 1678 – 1705.[32] Molchanov, I. S., Schmutz, M., 2010. Multivariate extension of put-call symmetry. SIAM J. Financial Math. 1 (1), 396–426.[33] Panini, R., 2004. Option pricing with Mellin transforms. Ph.D. thesis, State University of New York at Stony Brook.[34] Panini, R., Srivastav, R. P., 2004. Option pricing with Mellin transforms. Mathematical and Computer Modelling 40 (1-2),43–56.[35] Rodrigo, M. R., Mamon, R. S., 2007. An application of Mellin transform techniques to a Black-Scholes equation problem.Analysis and Applications 05 (01), 51–66.
Appendix A. Dishon and Weiss Method
The Mellin transform is equivalent to the two-sided Laplace transform under a negative logarithmicchange of variables. By exploiting this relationship, numerical Mellin inversion is possible via a seriesexpansion of sine and cosine functions as in [15]. By letting S = e − x for − L ≤ x ≤ L , equation (3.8) can beadapted to use this scheme. The value of the American option is the sum of the European option and early13xercise premium given by V EP ( S, t ) = e ax L ˆ g ( a ) + e ax L N X j =1 (cid:26) ℜ h ˆ g (cid:16) a + πijL (cid:17)i cos (cid:16) πjxL (cid:17) − ℑ h ˆ g (cid:16) a + πijL (cid:17)i sin (cid:16) πjxL (cid:17)(cid:27) (A.1)and V EEPP ( S, t ) = e ax L ˆ h ( a ) + e ax L N X j =1 (cid:26) ℜ h ˆ h (cid:16) a + πijL (cid:17)i cos (cid:16) πjxL (cid:17) − ℑ h ˆ h (cid:16) a + πijL (cid:17)i sin (cid:16) πjxL (cid:17)(cid:27) (A.2)respectively. By defining the time change τ = T − t and imposing the Trapezoid rule we obtain ˆ g ( w ; τ ) =exp( − rτ )ˆ θ ( w )Φ( wi ; τ ) and ˆ h ( w ; τ ) = ∆ t P M − l =0 exp( − rt l ) ˆ f ( w, τ − t l )Φ( wi ; t l ). The time integral in (3.16) isapproximated with truncation M and stepsize ∆ t = τ /M . To achieve a faster rate of convergence, L shouldbe chosen so that | x/L | ≤ / aa