A Finite Element Approach to the Numerical Solutions of Leland's Mode
AA Finite Element Approach to the Numerical Solutions ofLeland’s Model
Dongming Wei , Yogi Ahmad Erlangga , and Gulzat Zhumakhanova Nazarbayev University, Department of Mathematics, School of Sciences andHumanities, 53 Kabanbay Batyr Ave, Nur-Sultan 010000, Kazakhstan Zayed University, Department of Mathematics, College of Natural and Health Sciences,Abu Dhabi Campus, P.O. Box 144534, United Arab Emirates [email protected], [email protected], [email protected] * Corresponding author
Abstract
In this paper, finite element method is applied to Leland’s model for numerical simulation ofoption pricing with transaction costs. Spatial finite element models based on P1 and/or P2 elementsare formulated in combination with a Crank-Nicolson-type temporal scheme. The temporal schemeis implemented using the Rannacher approach. Examples with several sets of parameter values arepresented and compared with finite difference results in the literature. Spatial-temporal mesh-sizeratios are observed for controlling the stability of our method. Our results compare favorably withthe finite difference results in the literature for the model.
Keywords—
Option pricing, Leland’s model, finite element
A fair option price in a complete financial market with no transaction costs can be modelled by the Black-Scholesequation [7, 23]. The underlying assumption requires, however, that portfolio hedging takes place continuously.In the market with transaction costs, this assumption becomes unrealistically expensive. Modifications to theBlack-Scholes model have been proposed to count for the transaction costs, which lead to various nonlinearmodels [19, 6, 18, 11, 16].In [19], Leland proposes a modification to the Black-Scholes equation by allowing portfolio rebalancing at adiscrete time δt with the transaction costs proportional to the value of the underlying asset. For an Europeancall option with the strike price K and expiration time T , its price V at any time t can then be modelled byusing the nonlinear partial differential equation: V t + 12 σ S (1 + Le sgn( V SS )) V SS + rSV S − rV = 0 , in (0 , T ) × R + (1)where V = V ( S, t ), S is the value of the underlying asset, Le = (cid:114) π cσ √ δt is the Leland number, r is the risk-freeinterest rate, c is the round trip of the transaction cost per currency unit, and σ is the volatility. For (1), inaddition to the boundary conditions V (0 , t ) = 0 , (2) V ( S, t ) = S − K, as S → ∞ , (3)for all t ∈ [0 , T ], the following terminal condition at the expiration time T is also required: V ( S, T ) = max( S − K, . (4)The condition (4) is referred to as the pay-off function. a r X i v : . [ q -f i n . C P ] O c t few remarks are in order. Firstly, without transaction cost, V SS > V SS also holdsin the nonlinear case, then we have˜ σ := 12 σ (1 + Le sgn( V SS )) = 12 σ (1 + Le );Thus (1) is reduced to the Black-Scholes equation with an adjusted constant volatility ˜ σ > σ . Secondly, theLeland number depends highly on the transaction cost, which is typically a small percentage of the value ofthe assets and the time between rehedging, which is much smaller than the time to expiration. A small Lelandnumber, usually when Le <
1, corresponds to a small transaction cost or large time interval between rehedging,which is considered high risk. In this case, the nonlinear terminal-boundary value problem (1) is known tobe well-posed and has a unique solution V ( S, T ) for any terminal condition [5]. See also, e.g., [21, 27] formathematical analysis on the solution of (1).Practical option pricing is typically done by solving the underlying terminal-boundary value problem nu-merically. Popular numerical methods for this purpose are based on finite difference methods (FDM) and finiteelement (FEM) methods [1]. FDM are particularly popular for both linear and nonlinear cases due to the sim-plicity of the methods, especially when the computation is performed on a uniform mesh [4, 9, 28]. Developmentof high-order methods as well as mesh adaptivity used to control numerical errors may however not be triviallydone with FDM [10, 17, 20, 14]. These are not an issue with FEM, even though the implementation is morecomplex than FDM [24]. While FEM have been demonstrated to be a viable alternative to FDM in the linearcases [22, 3, 13], only limited work is presently done on the nonlinear cases, especially involving transaction costsunder Leland’s model (see [2]).Our aim with this paper is to present some effective finite element methods for the Leland model and todemonstrate that the finite element method is also a practical choice for numerical simulations of problems incomputational finance with Leland’s model.We note here that since the initial work of Leland’s, several authors have proposed modifications to theoriginal Leland model (1), to better capture the hedging strategy under transaction costs. These include, tomention a few, the model of Boyle and Vorst [8], Hoggard [15], and Zhao and Ziemba [29]. Such a modificationis reflected in the adjusted volatility, which, however, shares a nonlinear term in common: the signum functionsgn( V SS ). Along this line, [5] notes that, for Le ≥
1, Leland’s model (1) is ill-posed under nonconvex terminalconditions and proposes a hedging strategy that fixes this issue. Despite all of these variants, in this paper, weshall focus only on the application of the finite element method for the original Leland model (1). Applicationof the finite-element methods to the aforementioned variants can be easily done.The remainder of the paper is organized as follows. After introducing transformation to Leland’s model (1) inSection 2, we discuss a finite-element method and treatment for the signum term in Section 3. Section 4 discussesthe time-integration method. Numerical results are presented in Section 5, followed by concluding remarks inSection 6.
Following the standard strategy for solving the Black-Scholes-type problems, we first transform the terminal-boundary value problem (1) and (4) to an initial-boundary value problem of a simplified differential equation byusing the following change of variables: • τ = 12 σ ( T − t ) and hence t = T − τσ ; • x = ln( S ) + kτ and hence S = e x − kτ ; and • u ( x, τ ) = e kτ V ( S, t ), and thus V ( S, t ) = e − kτ u ( x, τ ).One can then show that • V t = σ e − kτ ( − ku x + ku − u τ ), • V S = 1 S e − kτ u x , and • V SS = 1 S e − kτ ( u xx − u x ).Substitution of the above derivatives of the option price V into (1) yields the transformed Leland model u τ = u xx − u x + Le | u xx − u x | . (5)Applying the above change of variables to the terminal and boundary conditions leads to1. the initial condition: u ( x,
0) = max( e x − K, . the boundary conditions • u ( x, τ ) = 0 as x → −∞ • u ( x, τ ) = e x − K as x → ∞ .As the problem is now defined in the unbounded spatial domain ( −∞ , + ∞ ), for computational purposes, wetruncate the solution domain to Ω = [ − R, R ], where 0 < R < ∞ and R is taken to be a large number. We enforcethe condition at −∞ to be satisfied at x = − R , and similarly for the other boundary condition. For the construction of the finite element approximation of the problem, we consider the following mixed formu-lation of (5): u τ = v + Le | v | , (6) v = u xx − u x . (7)The Galerkin’s finite element method for this formulation is based on the following weak form of (6) and (7).We begin by multiplying the equations (6) and (7) by the test functions w and z , respectively, and integrateeach over the domain Ω to get (cid:90) Ω w ( u τ − ( v + Le | v | ) dx = 0 , (cid:90) Ω z ( v − ( u xx − u x )) dx = 0 , which, after integration by parts, can be written as ∂∂τ (cid:90) Ω wudx − (cid:90) Ω ( wv + Lew | v | ) dx = 0 , (8) (cid:90) Ω zvdx + (cid:90) Ω z x u x dx + (cid:90) Ω zu x dx = 0 . (9)Let u = n (cid:88) i =1 u i ψ i + (cid:88) i ∈I ∂ Ω u i ψ i , I ∂ Ω = { , n + 1 } , be the finite element approximation of the solution u , where thesecond sum is the extension of the solution to the boundary ∂ Ω = {− R, R } and ψ i is the global finite elementshape function for the i th node in a spatial division − R = x < ... < x i < ... < x n +1 = R . Similarly, we have v = n +1 (cid:88) i =0 v i φ i , in which no boundary conditions are set for v . Then (8) can be written as0 = ∂∂τ (cid:90) Ω w n (cid:88) i =1 u i ψ i + (cid:88) i ∈I ∂ Ω u i ψ i dx − (cid:90) Ω ( w n (cid:88) i =1 v i φ i + Lew | n +1 (cid:88) i =0 v i φ i | ) dx = ∂∂τ n (cid:88) i =1 u i (cid:90) Ω wψ i dx + ∂∂τ (cid:88) i ∈I ∂ Ω u i (cid:90) Ω wψ i dx − n (cid:88) i =1 v i (cid:90) Ω wφ i dx − Le (cid:90) Ω w | n +1 (cid:88) i =0 v i φ i | dx. Enforcing this condition to be satisfied by n functions w j , j = 1 , . . . , n yields a system of n equations0 = ∂∂τ n (cid:88) i =1 u i (cid:90) Ω w j ψ i dx + ∂∂τ (cid:88) i ∈I ∂ Ω u i (cid:90) Ω w j ψ i dx − n (cid:88) i =1 v i (cid:90) Ω w j φ i dx − Le (cid:90) Ω w j | n (cid:88) i =1 v i φ i | dx. (10) or (9), with u x = n (cid:88) i =1 u i ψ i,x + (cid:88) i ∈I ∂ Ω u i ψ i,x , we have0 = (cid:90) Ω z n +1 (cid:88) i =0 v i φ i dx + (cid:90) Ω z x n (cid:88) i =1 u i ψ i,x + (cid:88) i ∈I ∂ Ω u i ψ i,x dx + (cid:90) Ω z n (cid:88) i =1 u i ψ i,x + (cid:88) i ∈I ∂ Ω u i ψ i,x dx = n (cid:88) i =1 v i (cid:90) Ω zφ i dx + n (cid:88) i =1 u i (cid:90) Ω z x ψ i,x dx + (cid:90) Ω zψ i,x dx + (cid:88) i ∈I ∂ Ω u i (cid:90) Ω z x ψ i,x dx + (cid:90) Ω zψ i,x dx . By enforcing the above equation to be satisfied by z j , j = 0 , . . . , n + 1 results in the system of n equations0 = n +1 (cid:88) i =0 v i (cid:90) Ω z j φ i dx + n (cid:88) i =1 u i (cid:90) Ω z j,x ψ i,x dx + (cid:90) Ω z j ψ i,x dx + (cid:88) i ∈I ∂ Ω u i (cid:90) Ω z j,x ψ i,x dx + (cid:90) Ω z j ψ i,x dx . (11)For our finite element models, we consider the Galerkin approach, where we set w i = z i = φ i = ψ i . Equa-tions (10) and (11) then become, for j = 1 , . . . , n∂∂τ n (cid:88) i =1 u i (cid:90) Ω ψ j ψ i dx + ∂∂τ (cid:88) i ∈I ∂ Ω u i (cid:90) Ω ψ j ψ i dx = n (cid:88) i =1 v i (cid:90) Ω ψ j ψ i dx + Le (cid:90) Ω ψ j | n +1 (cid:88) i =0 v i ψ i | dx. (12) n +1 (cid:88) i =0 v i (cid:90) Ω ψ j ψ i dx = − n (cid:88) i =1 u i (cid:90) Ω ψ j,x ψ i,x dx + (cid:90) Ω ψ j ψ i,x dx − (cid:88) i ∈I ∂ Ω u i (cid:90) Ω ψ j,x ψ i,x dx + (cid:90) Ω ψ j ψ i,x dx . (13)Let the domain Ω be subdivided into n E nonoverlapping elements such that Ω = n +1 (cid:83) i =1 Ω i , where Ω i = [ x i − , x i ],the i -th element with boundary nodes x i − and x i . In this way, each integral above can be written as the sumof integral over each element. For instance (cid:90) Ω ψ j ψ i dx = n (cid:88) (cid:96) =1 (cid:90) Ω (cid:96) ψ j ψ i dx, and so on. Thus, the integral over the domain Ω can be evaluated by first evaluating integral over elements andthen summing up, a process referred to as “assembly”. In the implementation, the assembly process is based onelement matrices that represents integral terms in (12) and (13) over each element Ω j . Structures of the elementmatrices depend on the specific choice of the functions ψ i . Specifically, the interpolation functions ψ i are chosensuch that, at the nodal points x j , ψ i ( x j ) = (cid:40) , i = j, , otherwise . (14)In this way, at the left boundary point x = − R , u ( x ) = (cid:88) i = n u i ψ ( x ) + (cid:88) i ∈I ∂ Ω u i ψ ( x ) = u = 0 . Similarly at the right boundary point x n = R , u ( x n ) = (cid:88) i = n u i ψ ( x n ) + (cid:88) i ∈I ∂ Ω u i ψ ( x n ) = u n = e R − K. In the sequel, we discuss two interpolation functions used in our finite element methods. .1 P1 finite element In the element Ω j = [ x j − , x j ], with the meshsize h j = x j − x j − , we define two interpolation basis function: ψ j − ( x ) = ( x − x j ) / ( x j − − x j ) = − ( x − x j ) /h j , (15) ψ j ( x ) = ( x − x j − ) / ( x j − x j − ) = ( x − x j − ) /h j . (16)This is a linear (Lagrange) interpolation polynomial, which leads to the P1 (linear) finite element.For the (cid:90) ψ j ψ i dx term, the element matrix reads M j = (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j dx (cid:90) Ω j ψ j ψ j − dx (cid:90) Ω j ψ j ψ j dx = h j (cid:20) (cid:21) . For the − (cid:90) ψ j,x ψ i,x dx term, the element matrix reads K j = − (cid:90) Ω j ψ j − ,x ψ j − ,x dx (cid:90) Ω j ψ j − ,x ψ j,x dx (cid:90) Ω j ψ j,x ψ j − ,x dx (cid:90) Ω j ψ j,x ψ j,x dx = − h j (cid:20) − − (cid:21) . For the (cid:90) ψ j ψ i,x dx term, the element matrix reads P j = (cid:90) Ω j ψ j − ψ j − ,x dx (cid:90) Ω j ψ j − ψ j,x dx (cid:90) Ω j ψ j ψ j − ,x dx (cid:90) Ω j ψ j ψ j,x dx = 12 (cid:20) − − (cid:21) . For the nonlinear term (cid:90) ψ j | (cid:88) v i ψ i | dx term, since ψ j ≥ ψ j ψ i ≥
0. Thus, (cid:90) Ω ψ j | n +1 (cid:88) i =0 v i ψ i | dx = (cid:90) Ω | n +1 (cid:88) i =0 v i ψ j ψ i | dx = (cid:90) Ω | v j − ψ j ψ j − + v j ψ j ψ j + v j +1 ψ j ψ j +1 | dx = (cid:90) Ω j | v j − ψ j ψ j − + v j ψ j ψ j + v j +1 ψ j ψ j +1 | dx + (cid:90) Ω j +1 | v j − ψ j ψ j − + v j ψ j ψ j + v j +1 ψ j ψ j +1 | dx = (cid:90) Ω j | v j − ψ j ψ j − + v j ψ j ψ j | dx + (cid:90) Ω j +1 | v j ψ j ψ j + v j +1 ψ j ψ j +1 | dx ≈ | v j − | (cid:90) Ω j | ψ j ψ j − | dx + | v j | (cid:90) Ω j | ψ j ψ j | dx + (cid:90) Ω j +1 | ψ j ψ j | dx + | v j +1 | (cid:90) Ω j +1 | ψ j ψ j +1 | dx. = | v j − | (cid:90) Ω j ψ j ψ j − dx + | v j | (cid:90) Ω j ψ j ψ j dx + (cid:90) Ω j +1 ψ j ψ j dx + | v j +1 | (cid:90) Ω j +1 ψ j ψ j +1 dx. he corresponding element matrix for the element Ω j with nodal solution values | v j − | and | v j | is given by¯ M j = (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j dx (cid:90) Ω j ψ j ψ j − dx (cid:90) Ω j ψ j ψ j dx = h j (cid:20) (cid:21) . Notice that, for the P1 finite element, ¯ M j = M j . In the basic element Ω j , we add a midpoint x j − , hence x j − x j − = h j /
2, and define three quadratic interpolationpolynomials ψ j − ( x ) = ( x − x j − )( x − x j )( x j − − x j − )( x j − − x j ) = 2( x − x j − )( x − x j ) /h j , (17) ψ j − ( x ) = ( x − x j − )( x − x j )( x j − − x j − )( x j − − x j ) = − x − x j − )( x − x j ) /h j , (18) ψ j ( x ) = ( x − x j − )( x − x j − )( x j − x j − )( x j − x j − ) = 2( x − x j − )( x − x j − ) /h j . (19)The resulting finite element method is referred to as the P2 finite element.For the (cid:90) ψ j ψ i dx term, the element matrix reads M j = (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j dx (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j − dx (cid:90) Ω j ψ j − ψ j dx (cid:90) Ω j ψ j ψ j − dx (cid:90) Ω j ψ j ψ j − dx (cid:90) Ω j ψ j ψ j dx = h j −
12 16 2 − . For the − (cid:90) ψ j,x ψ i,x dx term, the element matrix reads K j = − (cid:90) Ω j ψ j − ,x ψ j − ,x dx (cid:90) Ω j ψ j − ,x ψ j − ,x dx (cid:90) Ω j ψ j − ,x ψ j,x dx (cid:90) Ω j ψ j − ,x ψ j − ,x dx (cid:90) Ω j ψ j − ,x ψ j − ,x dx (cid:90) Ω j ψ j − ,x ψ j,x dx (cid:90) Ω j ψ j,x ψ j − ,x dx (cid:90) Ω j ψ j,x ψ j − ,x dx (cid:90) Ω j ψ j,x ψ j,x dx = − h j − − − − . For the (cid:90) ψ j ψ i,x dx term, the element matrix reads P j = (cid:90) Ω j ψ j − ψ j − ,x dx (cid:90) Ω j ψ j − ψ j − ,x dx (cid:90) Ω j ψ j − ψ j,x dx (cid:90) Ω j ψ j − ψ j − ,x dx (cid:90) Ω j ψ j − ψ j − ,x dx (cid:90) Ω j ψ j − ψ j,x dx (cid:90) Ω j ψ j ψ j − ,x dx (cid:90) Ω j ψ j ψ j − ,x dx (cid:90) Ω j ψ j ψ j,x dx = 16 − − − − . To construct the element matrix corresponding to the (cid:90) ψ j | (cid:88) v i ψ i | dx term, we begin with the splitting ψ j ( x ) = ψ + j ( x ) − | ψ − j ( x ) | here ψ + j is the nonnegative part of ψ j and ψ − j is the nonpositive part of ψ j . Then (cid:90) Ω ψ j | n (cid:88) i =1 v i ψ i | dx = (cid:90) Ω ( ψ + j − | ψ − j | ) | n (cid:88) i =1 v i ψ i | dx = (cid:90) Ω | n (cid:88) i =1 v i ψ + j ψ i | dx − (cid:90) Ω | n (cid:88) i =1 v i | ψ − j | ψ i | dx ≈ (cid:90) Ω n (cid:88) i =1 | v i || ψ + j ψ i | dx − (cid:90) Ω n (cid:88) i =1 | v i || ψ − j ψ i | dx = (cid:90) Ω j (cid:16) | v j − || ψ + j ψ j − | + | v j − || ψ + j ψ j − | + | v j || ψ + j ψ j | (cid:17) dx + (cid:90) Ω j +1 (cid:16) | v j || ψ + j ψ j | + | v j + || ψ + j ψ j + | + | v j +1 || ψ + j ψ j +1 | (cid:17) dx − (cid:90) Ω j (cid:16) | v j − || ψ − j ψ j − | + | v j − || ψ − j ψ j − | + | v j || ψ − j ψ j | (cid:17) dx − (cid:90) Ω j +1 (cid:16) | v j || ψ − j ψ j | + | v j + || ψ − j ψ j + | + | v j +1 || ψ − j ψ j +1 | (cid:17) dx = | v j − | (cid:90) Ω j (cid:0) | ψ + j ψ j − | − | ψ − j ψ j − | (cid:1) dx + | v j − | (cid:90) Ω j (cid:16) | ψ + j ψ j − | − | ψ − j ψ j − | (cid:17) dx + | v j | (cid:90) Ω j (cid:0) | ψ + j ψ j | − | ψ − j ψ j | (cid:1) dx + (cid:90) Ω j +1 (cid:0) | ψ + j ψ j | − | ψ − j ψ j | (cid:1) dx + | v j + | (cid:90) Ω j (cid:16) | ψ + j ψ j + | − | ψ − j ψ j + | (cid:17) dx + | v j +1 | (cid:90) Ω j (cid:0) | ψ + j ψ j +1 | − | ψ − j ψ j +1 | (cid:1) dx. The above leads to the element matrix¯ M j = (cid:90) Ω j (cid:0) | ψ + j − ψ j − | − | ψ − j − ψ j − | (cid:1) dx (cid:90) Ω j (cid:16) | ψ + j − ψ j − | − | ψ − j − ψ j − | (cid:17) dx (cid:90) Ω j (cid:0) | ψ + j − ψ j | − | ψ − j − ψ j | (cid:1) dx (cid:90) Ω j | ψ + j − ψ j − | dx (cid:90) Ω j | ψ + j − ψ j − | dx (cid:90) Ω j | ψ + j − ψ j | dx (cid:90) Ω j (cid:0) | ψ + j ψ j − | − | ψ − j ψ j − | (cid:1) dx (cid:90) Ω j (cid:16) | ψ + j ψ j − | − | ψ − j ψ j − | (cid:17) dx (cid:90) Ω j (cid:0) | ψ + j ψ j | − | ψ − j ψ j | (cid:1) dx . Note that for φ − j − = 0. Evaluating the integrals results in the element matrix¯ M j = h j
15 8 08 64 80 8 15 , associated with the nodal solution | v j − | , | v j − / | , and | v j | . In this case, ¯ M (cid:54) = M . The global finite element system can be assembled into the following first order ODE in matrix form ∂∂τ ( M u + b M ) = M v + Le ¯ M | v | := F , (20)where v = M − ( K u − P u + b K − b P ), and all boldface letters are used to denote the solution nodal value vectorsor the finite element matrices. For an implementation of the assembly process, see, e.g., [25]. ntegration over time is approximated using the Crank-Nicolson-type scheme:1∆ τ ( M u n +1 + b n +1 M − ( M u n + b nM )) = θ F n +1 + (1 − θ ) F n , with θ ∈ [0 , θ = 0 and 1 correspond to the explicit forward and implicit backward Euler method, respectively.Rearranging the above equation leads to M u n +1 − θ ∆ τ F n +1 = M u n + (1 − θ )∆ τ F n − b n +1 M + b nM or M u n +1 − θ ∆ τ ( M v n +1 + Le ¯ M | v n +1 | ) = M u n + (1 − θ )∆ τ F n − b n +1 M + b nM . With v n +1 = M − ( K u n +1 − P u n +1 + b n +1 K − b n +1 P n ), we get M u n +1 − θ ∆ τ ( K u n +1 − P u n +1 + Le ¯ M | v n +1 | ) = M u n + (1 − θ )∆ τ F n − b n +1 M + b nM + θ ∆ τ ( b n +1 K − b n +1 P ) . This equation is nonlinear in the ( n + 1) variables. Setting | v n +1 | = | v n | results in the linearized form: M u n +1 − θ ∆ τ ( K u n +1 − P u n +1 ) = M u n + (1 − θ )∆ τ F n + θ ∆ τ Le ¯ M | v n | − b n +1 M + b nM + θ ∆ τ ( b n +1 K − b n +1 P )or A u n +1 = M u n + (1 − θ )∆ τ F n + θ ∆ τ Le ¯ M | v n | − b n +1 M + b nM + θ ∆ τ ( b n +1 K − b n +1 P )with A = M − θ ∆ τ ( K − P ) , F n = M v n + Le ¯ M | v n | , v n = M − ( K u n − P u n + b nK − b nP )Recall that for the P2 FEM, ¯ M (cid:54) = M . As another level of approximation, we may set also ¯ M = M . Werefer to the numerical method described above as Version 1, and with the setting ¯ M = M is Version 2. Forimproved stability, the Crank-Nicolson method is implemented using the Rannacher approach [26, 12], in whichthe first Crank-Nicolson step is replaced by a few backward implicit Euler steps with smaller time steps (e.g.,∆ τ R = ∆ τ /n R ), where n R is the number of backward Euler time steps in from τ = 0 to τ = ∆ τ ). In this section, we present some numerical solutions computed using the FEM and the time-integration methoddiscussed in Sections 3 and 4. All results are computed on a uniform finite-element mesh, even though the methodcan be implemented on a nonuniform mesh.
In the absence of closed-form exact solution, to verify the accuracy of the methods, we first conducted a numericaltest using a linear case, which corresponds to the situation where no rebalancing involves. This is equivalent tohaving δt → ∞ and therefore Le = 0. For this test, we set r = 0 . σ = 0 .
2, the expiration time T = 1, and thestrike price K = 100.The numerically computed option price V is shown in Figure 1 as a function of asset price S . Comparing withthe exact solution, we see in the figure the highly accurate solution obtained by the linear (P1) and quadratic(P2) FEM. For the nonlinear Leland model, we consider cases, where (i)
Le <
1, and (ii)
Le >
1. We discuss the convergenceand numerical stability of our schemes for numerical solutions of the two these cases in the following:
The
Le < case. Our first test in this case is based on the parameters r = 0 . σ = 0 . T = 1, K = 100, c = 0 .
01, and δt = 0 .
01, yielding Le = 0 .
4. For the FEM, we set discretization parameters such that ∆ τ /h < τ /h <
1, to control the stability of the method. Otherwise mentioned, we keep ∆ τ /h = 0 .
01. Thenumerical solutions over time t ∈ [0 ,
1] are shown in Figure 2, for the P1 and P2 FEM. A look at the figure givesno indication of numerical instability.To have a clearer picture of the computed solutions, we show in Figure 3 the FEM solutions at expirationand compare them with a finite-difference (FDM) solution and the exact solution of the linear case. In this case, t ∈ [0 , r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 .
01, giving Le ≈ .
4; Left figure: h = 0 .
1, ∆ τ = 0 .
001 (thus, ∆ τ /h = 0 . ≤ h = 0 . τ = 0 . τ /h = 0 . ≤ τ /h = 0 . the option price based on Leland’s model is higher than that based on the linear Black-Scholes model. While thefinite-difference method used here is insensitive to the choice of the discretization parameters, a large differencein the solutions is observed for the finite-element method. This dependency on the discretization parameters isstronger in the P2 FEM than in the P1 FEM, whose solution remains close to the FDM solution. We also observethat setting ¯ M = M in the P2 FEM (Version 2) does not result in a solution which is significantly different fromVersion 1 (with M (cid:54) = M ).A second test in this class of problems uses the same parameters, except for c , which is now set to 0 . Le = 0 .
8. Numerical solutions overthe time t ∈ [0 ,
1] (not shown in this paper) do not indicate any instability, as in the previous case. A closerobservation of the solution at expiration also shows a higher price obtained by Leland’s model than that by thelinear Black-Scholes model (see Figure 4). The FEM results also show the strong-dependency on the numericalparameters, leading to a higher computed option price. Furthermore, Version 1 and Version 2 do not result insignificantly different solutions.
The
Le > case. For this class of problem, we set c = 0 .
03, while keeping the same values for the otherparameters as in the
Le < Le = 1 .
2. Since Version 1 and 2 do not lead to a significantdifference in the solutions, we therefore will only implement and show results using Version 1 in the next numericaltests.Figure 5 presents numerical results using the P1 FEM, with ∆ t/h = 0 . . t/h = 0 .
8, as t gets closer to r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 .
01, giving Le ≈ . Left : h = 0 .
1, ∆ τ = 0 .
001 (∆ τ /h = 0 .
01, ∆ τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 .
01, ∆ τ /h = 0 . r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 .
01, giving Le ≈ . Left : h = 0 .
1, ∆ τ = 0 .
001 (∆ τ /h = 0 .
01, ∆ τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 .
01, ∆ τ /h = 0 . t/h = 0 . τ /h , as demonstrated by the numerical results in Figure 6.Using the above-mentioned test setting, we compute the solution with the P2 FEM. Numerical solutions with∆ τ /h = 0 . τ /h = 0 .
8, with spurious oscillation on the coarsespatial-temporal mesh. Such an instability can however still be controlled by decreasing the spatial and temporalstep, even though in the current test, it is not fully eliminated.
In this work, several finite element-based models for approximations of the solutions of the Leland’s modelare built in combination with the Crank-Nicolson-type temporal scheme. It is demonstrated through severalnumerical examples that stable and accurate solutions can be obtained by these models by controlling the spatialfinite element size and the temporal steps. These numerical results compare favorably with those computed byfinite difference schemes. Our finite element models can be used as an effective alternative for numerical solutionsof the Leland model including the standard Black-Scholes model. r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 . Le ≈ . Left : h = 0 .
1, ∆ τ = 0 .
001 (∆ τ /h = 0 .
01, ∆ τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 .
01, ∆ τ /h = 0 . r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 . Le ≈ . Left : h = 0 .
05, ∆ τ = 0 . τ /h = 0 . τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 . τ /h = 0 . References [1] Y. Achdou and O. Pironneau.
Computational Methods for Option Pricing . SIAM, 2005.[2] R.M.P. Almeida, T.D.C. Chihaluca, and J.C.M. Duque. Hermite finite element method for nonlinear Black-Scholes equation governing European options. In J. Vigo-Aguiar et al., editor,
Proceedings of the 17thInternational Conference on Computational and Mathematical Methods in Science and Engineering , July2017.[3] A. Andalaft-Chacur, M.M Ali, and J.G. Salazar. Real options pricing by the finite element method.
Com-puters and Mathematics with Applications , 61:2863–2873, 2011.[4] J. Ankudinova and M. Ehrhardt. On the numerical solution of nonlinear Black-Scholes equations.
Computersand Mathematics with Applications , 56:799–812, 2008.[5] M. Avellaneda and A. Paras. Dynamic hedging portfolios for derivative securities in the presence of largetransaction costs.
Applied Mathematical Finance , 1:165–193, 1994.[6] G. Barles and H. Soner. Option pricing with transaction costs and a nonlinear Black-Scholes equation.
Finance and Stochastics , 2(4):369–397, 1998.[7] F. Black and M. Scholes. The pricing of options and corporate liabilities.
The Journal of Political Economy ,81:637–654, 1973. r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 . Le ≈ . Left : h = 0 .
1, ∆ τ = 0 .
001 (∆ τ /h = 0 .
01, ∆ τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 .
01, ∆ τ /h = 0 . r = 0 . σ = 0 . K = 100, c = 0 . δt = 0 . Le ≈ . Left : h = 0 .
05, ∆ τ = 0 . τ /h = 0 . τ /h = 0 . Right : h = 0 . τ = 0 . τ /h = 0 . τ /h = 0 . [8] P.P. Boyle and T. Vorst. Option replication in discrete time with transaction costs. The Journal of Finance ,47:271–293, 1992.[9] R. Company, L. J´odar, and J.-R. Pintos. A numerical method for European option pricing with transactioncosts nonlinear equation.
Mathematical and Computer Modelling , 50(5–6):910–920, 2009.[10] B. D¨uring, M. Fournier, and A. J¨ungel. High-order compact finite difference schemes for a nonlinear Black-Scholes equation.
International Journal of Theoretical and Applied Finance , 6(7):767–789, 2003.[11] R. Frey and A. Stremme. Market volatility and feedback effects from dynamic hedging.
MathematicalFinance , 7:351–374, 1997.[12] M.B. Giles and R. Carter. Convergence analysis of Crank–Nicolson and Rannacher time-marching.
Journalof Computational Finance , 9(4), 2006.[13] A. Golbabai, L.V. Ballestra, and D. Ahmadian. Superconvergence of the finite element solutions of theBlack–Scholes equation.
Finance Research Letters , 10:17–26, 2013.[14] S. Gulen, C. Popescu, and M. Sari. A new approach for the black–scholes model with linear and nonlinearvolatilities.
Mathematics , 7(8):760, 2019.[15] T. Hoggard, A. E. Whalley, and P. Wilmott. Hedging option portfolios in the presence of transaction costs.
Advances in Futures and Options Research , 7:21–35, 1994.[16] M. Jandaˇcka and D. ˇSevˇcoviˇc. On the risk-adjusted pricing-methodology-based valuation of vanilla optionsand explanation of the volatility smile.
Journal of Applied Mathematics , 3:235–258, 2005.
17] A.Q.M. Khaliq and W. Liao. High-order compact scheme for solving nonlinear Black-Scholes equation withtransaction costs.
International Journal of Computer Mathematics , 86:1009–1023, 2009.[18] M. Kratka. No mystery behind the smile.
Risk , 9:67–71, 1998.[19] H. Leland. Option pricing and replication with transactions costs.
The Journal of Finance , 40(5):1283–1301,1985.[20] G. Linde, J. Persson, and L. von Sydow. A highly accurate adaptive finite difference solver for theBlack–Scholes equation.
International Journal of Computer Mathematics , 86(12):2104–2121, 2009.[21] M.C. Mariani, E.K. Ncheuguim, and I. Sengupta. Solution to a nonlinear Black-Scholes equation.
ElectronicJournal of Differential Equations , 158:1–10, 2011.[22] S. Markolefas. Standard Galerkin formulation with high order Lagrange finite elements for option marketspricing.
Applied Mathematics and Computation , 195:707–720, 2008.[23] R. Merton. Theory of rational option pricing.
Bell Journal of Economics and Management Science , 4(1):141–183, 1973.[24] O. Pironneau and F. Hecht. Mesh adaption for the Black and Scholes equations.
East-West Journal ofNumerical Mathematics , 8(1), 1999.[25] Kythe K. Prem and D. Wei.
An Introduction to Linear and Nonlinear Finite Element Analysis: A Compu-tational Approach . Springer Science+Business Media, New York, 2004.[26] R. Rannacher. Finite element solution of diffusion problems with irregular data.
Numerische Mathematik ,43:309–327, 1984.[27] M. ˇZitˇnansk´a and D. ˇSevˇcoviˇc. Analysis of the nonlinear option pricing model under variable transactioncosts.
Asia-Pacific Financial Markets , 23(2):153–174, 2016.[28] W. Zhao, X. Yang, and L. Wu. Alternating segment explicit-implicit and implicit-explicit parallel differencemethod for the nonlinear Leland equation.
Advances in Difference Equations , 103, 2016. 18 pp, DOI:10.1186/s13662-016-0823-5.[29] Y. Zhao and W. T. Ziemba. On Leland’s option hedging strategy with transaction costs. Sauder School ofBusiness Working Paper, Available at SSRN: http://dx.doi.org/10.2139/ssrn.591661 , Aug 2004. 33 pp., Aug 2004. 33 pp.