Numerical investigation of Black-Scholes equation

Numerical Investigation of Black-Scholes Equation A dissertation submitted to the University of Dhaka in partial fulfillment of the requirements for the degree of Masters of Science in Mathematics SUBMITTED BY Examination Roll No. :2804 M.S. Session :2017-2018 Registration No. :2013-512-745 Department of Mathematics Faculty of Science University of Dhaka January, 2020 Acknowledgement I am privileged to have the opportunity to convey my respect and sincere gratitude to my supervisor Professor Dr. Samir Kumar Bhowmik, Department of Mathematics, University of Dhaka, whose encouraging guidance, inspiration and confidence in me have been my best motivation to accomplish this paper successfully. I would like to pay my respect and profound gratefulness to Jakobin Alam Khan, MS: 2016-2017, Department of Mathematics, whose kind guidance and benevolent assistance have helped me to be capable of composing this paper successfully. I convey my keen recognition and deep indebtedness to all academic and administrative staffs of Department of Mathematics, University of Dhaka for their cordial assistance. Finally, I am grateful and thankful to my beloved parents, brother, friends and all the wellwishers for their inspiration and faith in me, which made it possible for me to accomplish this far. i Abstract Black-Scholes Equation provides a theoretical estimate for the price of options. In this thesis work, we consider both Linear and Nonlinear Black-Scholes Model for the European option price. We analytically approach the Black-Scholes Equation by the transformation of the problem into a forward convection-diffusion equation with linear and nonlinear term. We use Spectral Collocation Method and Differential Quadrature Method for Linear Black-Scholes Equation and compare the numerical results with the exact solution of the equation. Then the Nonlinear Black-Scholes Equation has been investigated. Here, different volatility models due to the presence of transaction costs has been considered. We use Finite-Difference-Method to solve the Nonlinear Black-Scholes Equation. Finally we present the numerical results of the Nonlinear Black-Scholes Equation. We use MATLAB for numerical implementation and LATEX has been used for the write up. ii Contents Introduction viii 1 1 Black Scholes Equation 1 1.1 Introduction................................ 1 1.2 FinancialDerivatives........................... 1 1.3 OptionPrice................................ 2 1.3.1 Classification of Option Price . . . . . . . . . . . . . . . . . . 3 1.3.2 EuropeanOption......................... 4 1.4 Black-Scholes Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Limitations of Black-Scholes Equation . . . . . . . . . . . . . . . . . 5 1.6 Derivation of Black Scholes Equation By Hedging argument . . . . . 6 1.7 Derivation of Black-Scholes Equation by Binomial Model . . . . . . . 8 1.8 ParityofPut-Call............................. 9 1.9 SensitivityAnalysis............................ 9 2 Linear Black-Scholes equation’s Analytic Solution 12 2.1 Introduction................................ 12 2.2 Black Scholes Equation’s Linear Model . . . . . . . . . . . . . . . . . 12 2.3 AnalyticalSolution............................ 14 2.4 Graphical Solution of the Black-Scholes Equation . . . . . . . . . . . 19 3 Linear Black-Scholes Equation with Spectral Method 22 3.1 Introduction................................ 22 iii 3.2 Spectral Collocation Method . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Chebyshev Differentiation Matrix . . . . . . . . . . . . . . . . . . . . 23 3.4 MATLAB function : ’cheb’ . . . . . . . . . . . . . . . . . . . . . . . . 26 3.5 Numerical Solution of Black-Scholes Equation with the help of Spectral Method .................................. 29 3.6 Stability and Convergence of Spectral Method . . . . . . . . . . . . . 31 3.7 Conclusion................................. 32 4 Differential Quadrature Method for Linear Black-Scholes Equation 33 4.1 Introduction................................ 33 4.2 Differential Quadrature Method (DQM) . . . . . . . . . . . . . . . . 33 4.3 Discretization of domain and Stability . . . . . . . . . . . . . . . . . 34 4.4 Derivatives of Space by using Differential Quadrature Method . . . . 35 4.5 Matlab function : Differential Quadrature Matrix . . . . . . . . . . . 37 4.6 Numerically solving Black-Scholes PDE using Differential Quadrature Method .................................. 39 4.7 Conclusion................................. 41 5 Nonlinear Black Scholes Equation and Volatility Models 42 5.1 Introduction................................ 42 5.2 Nonlinear Black Scholes Equation . . . . . . . . . . . . . . . . . . . . 42 5.3 BoundaryConditions .......................... 43 5.3.1 European Call Option . . . . . . . . . . . . . . . . . . . . . . 43 5.3.2 European Put Option . . . . . . . . . . . . . . . . . . . . . . . 43 5.4 VolatilityModels ............................. 43 5.5 Leland’s Model for Volatility . . . . . . . . . . . . . . . . . . . . . . . 45 5.6 Boyle and Vorst Volatility Model . . . . . . . . . . . . . . . . . . . . 47 5.7 Paras and Avellaneda Volatility Model . . . . . . . . . . . . . . . . . 47 5.8 Hodges and Neuberger Volatility Model . . . . . . . . . . . . . . . . . 47 5.9 Barles’ and Soner Model . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.10 Risk adjusted Pricing Volatility Model . . . . . . . . . . . . . . . . . 51 iv 5.11AnalyticalSolution............................ 52 5.12Conclusion................................. 53 6 Numerical Investigation of Nonlinear Black-Scholes Equation 54 6.1 EuropeanCalloption........................... 55 6.2 Finite-Difference Methods . . . . . . . . . . . . . . . . . . . . . . . . 55 6.3 Grid .................................... 56 6.4 Difference Quotients . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.5 VolatilityModels ............................. 59 6.6 Existence and Convergence . . . . . . . . . . . . . . . . . . . . . . . . 60 6.7 Explicit Method(Forward-Time Central Space) . . . . . . . . . . . . . 63 6.8 Implicit Method (Backward Time-Central Time) . . . . . . . . . . . . 65 6.9 Crank-Nicolson .............................. 66 6.10 Graphical Representation . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.11Conclusion................................. 73 7 Fractional Black-Scholes Equation 74 7.1 Introduction................................ 74 7.2 Transformed Black-Scholes Equation . . . . . . . . . . . . . . . . . . 74 7.3 Fractional Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 7.4 Fractional Black-Scholes Equation of Time . . . . . . . . . . . . . . . 75 7.5 Fractional Black-Scholes Equation of Space . . . . . . . . . . . . . . . 77 7.6 Conclusion................................. 79 8 Conclusion 80 Bibliography 81 v List of Figures 1.1 European Call Option Price V(S, t) at strike price K= 100 . . . . . . 4 1.2 European Put Option Price V(S, t) at strike price K= 100 . . . . . . 5 2.1 Value of the European call option V(S, t) at time, t= 0, t=.2, t=.4, t=.6, t=.8 and t=1........................... 20 2.2 Value of the European call option V(S, t) at various times in a surface plot .................................... 20 3.1 European Call Option Price V(S, t) using Spectral Collocation Method. 30 3.2 The comparison of European Call Option Price V(S, t) with the Exact solution................................... 30 4.1 European Call Option Price V(S, t) using Differential Quadrature Method. 40 4.2 European Call Option Price V(S, t) is compared with the Exact solution. 40 6.1 Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.0001 with the Crank-Nicolson Method. . . . . . . . . . . . . 68 6.2 The impact of transaction costs Nonlinear Vnonlinear(S, t)−Vlinear(S, t) with the Crank-Nicolson Method. . . . . . . . . . . . . . . . . . . . . 69 6.3 Value of European Call Option price V(S, 0) for different volatility models vs the linear model with the Crank-Nicolson Method. . . . . . 70 6.4 Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.001 with the Implicit Method. . . . . . . . . . . . . . . . . 71 vi 6.5 The Identity model ψ(x) = xwith parameter (a=.02) for the (Vnonlinear(S, t)− Vlinear(S, t)) plot for spatial step h=.1 and time step k=.001 with theImplicitMethod. ........................... 71 6.6 Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.00001 with the Explicit Method. . . . . . . . . . . . . . . . 72 6.7 The Identity model ψ(x) = xwith parameter (a=.02) for the (Vnonlinear(S, t)− Vlinear(S, t)) plot for spatial step h=.1 and time step k=.001 with theExplicitMethod. ........................... 72 vii Introduction During the last several decades the trading of options and stocks has encountered a growing interest in both scientific work and everyday life. Fischer Black and Myron Scholes published a paper titled ’The Pricing of Options and Corporate Liabilities’ in 1973. They showed how the pricing of options is uniquely determined by their formula. After their publication of this paper, the formula was named after them. The arbitrage reasoning method was used in their paper which was developed by Robert Merton. They were given the Nobel prize in 1997 for this Option Pricing model. This Black-Scholes Model (BSM) was the first model for pricing options. After this model, there were more models developed by the mathematicians to estimate stock option price. There are some limitations to the BSM. That is why the classical model sometimes differ with the realistic price of the options [3, 44]. There had been numerous attempts on improving the BSM [28, 35]. In the Linear BSM, the volatility (the standard deviation of the stock prices) is assumed to be constant. This assumption is not realistic. That is why in the Nonlinear Black-Scholes Model (NBSM), the volatility term is not constant. Therefore, NBSM is less error prone than the linear one. There have been various schemes used to solve BSM. If one wants to solve an ODE or PDE to high accuracy on a simple domain, and if the data defining the problem are smooth, then spectral methods are usually the best tool [67]. This method has the ability to achieve ten digits of accuracy where a finite difference scheme would get two or three. This method demands less computer memory than the alternatives and was developed by Steven Orszag in 1969. Lloyd N. Trefethen’s book named ”Spectral Methods in MATLAB” provides fundamental ideas and techniques of spectral methods. This book offers a remarkable range of ODE and PDE problems which can be viii 1.3.2 European Option European Call Option is a contract between two parties where at a fixed time in the future, called expiration date or maturity T, the possessor of the option may buy an asset known as the underlying asset S(t), for a fixed price called the strike or exercise price K. The other party has the agreement to sell the asset if the possessor wishes to buy it. We have already discussed about the three cases for Options. From them we can easily deduce that the value of the European Call option at expiry, called the pay-off function: V(S, T) = max(S−K, 0). Figure 1.1: European Call Option Price V(S, t) at strike price K= 100 European Put Option is the authority to sell the underlying asset S(t) at the maturity date Tfor the strike price K. The pay-off function for European Put Option is: V(S, T) = max(K−S, 0). 4 Figure 1.2: European Put Option Price V(S, t) at strike price K= 100 1.4 Black-Scholes Equation In the Inception of the 1970’s Fisher Black, Myron Scholes and Robert Merton developed the classical Black-Scholes Equation(PDE). It was a quintessential advancement in the field of quantitative finance. There are two classes of Black-Scholes PDE which are Linear and Nonlinear. Classical model always refers to linear Black-Scholes Equation.The Black-Scholes PDE is dependent on two independent variables which are time tand Stock price S(t) . The Black-Scholes PDE is of the parabolic form. It is called the backward parabolic form because of the presence of having the opposite sign on the second derivative. There are some assumptions involved in Black-Scholes PDE. As well as ,there are some limitations in this option pricing model. In this work, we reduce the Black-Scholes PDE into heat equation. Then we use numerical schemes to solve it. 1.5 Limitations of Black-Scholes Equation There are some limitations of Black-Scholes model because of some assumptions [4,7]: •Throughout the option’s lifetime, the stock pays no dividends. But in reality most of the enterprises give the payment of dividends to their share holder. Therefore, this is a significant constraint to the model. 5 •There are no commissions charged in this model. Usually the participants of the market have to pay a commision to buy or sell options. Even floor traders have to pay a small fee. The commission paid by the individual trader’s is more significant which can sometimes deform the outcome of the Black-Scholes Model. •Interest rate and volatility are known and constant. •The trading of assets is a continous process. •Transaction costs are ignored. •Normal distribution is followed for underlying asset price. •It assumes that options has the similar properties of the European option price. 1.6 Derivation of Black Scholes Equation By Hedging argument We generally presume that the asset price S(t) always follows a Geometric Browninan motion: dS =µSdt +σSdW. or dS =S(µdt +σdW). or dS S=µdt +σdW. where Wis the stochastic variable. Wand dW operate as unpredictability in stock market. And µdt works as the expected value whereas σ2dt is the variance. From Ito’s Lemma, we can write: dV =∂V ∂t +1 2σ2S2∂2V ∂S2+µS ∂V ∂S dt +σS ∂V ∂S dW. We can rewrite this as: ∆V=∂V ∂t +1 2σ2S2∂2V ∂S2+µS ∂V ∂S ∆t+σS ∂V ∂S ∆W. 6 We can rewrite dS =µSdt +σSdW this as: ∆S=µS∆t+σS∆W Then we are going to construct a portfolio Π which can be self financed. This portfolio is composed of an amount ∆ of the underlying stock, such that the portfolio is riskless and one option [62, 19, 10]. Therefore, the value of the portfolio at time tis Π(t) = V(t)+∆S(t).(1.1) Hence we can write, dΠ =dV + ∆dS (1.2) =(∂V ∂t +µS ∂V ∂S +1 2σ2S2∂2V ∂S2+ ∆µS)dt + (σS ∂V ∂S + ∆σS)dW. (1.3) There are two properties of the portfolio. The first property states that it must riskfree which indicates that the 2nd term involving the Brownian Motion dW is zero , so that σS ∂V ∂S + ∆σS = 0, which becomes ∆ = −∂V ∂S . Now we can Substitute ∆ in Equation (1.2) to get dΠ = ∂V ∂t +1 2σ2S2∂2V ∂S2. The second property is that the portfolio has to obtain the risk free rate. dΠ =rΠdt (1.4) =(∂V ∂t +1 2σ2S2∂2V ∂S2)dt =r(V−∂V ∂S S)dt (1.5) We can easily rearrange (1.5) to get ∂V ∂t +1 2σ2S2∂2V ∂S2+rS ∂V ∂S −rV = 0. This is called the Black-Scholes Partial Differential Equation. 7 1.7 Derivation of Black-Scholes Equation by Binomial Model The stock price at the time tis Stdefined as, u= expσ√dt At the time t+dt the the stock price goes up to Su t+dt =uSt, having the probability p=exprdt −d u−d The derivative has risk-less value which gives a relation [19], Vexprdt =pVu+ (1 −p)Vd=p(Vu−V d) + Vd.(1.6) where V=V(St), Vu=V(Su t+dt and Vd=V(Sd t+dt). Now we will use Taylor series expansions of Vu, Vd,exprdt, u and dup to order dt. From Taylor series we have, Vu≈V+∂V ∂S (Su t+dt −St) + 1 2 ∂2V ∂S2(Su t+dt −St)2+∂V ∂t dt (1.7) =V+∂V ∂S St(u−1) + 1 2 ∂2V ∂S2S2 t(u−1)2+∂V ∂t dt. (1.8) We similarly get, Vd≈V+∂V ∂S St(d−1) + 1 2 ∂2V ∂S2S2 t(d−1)2+∂V ∂t dt. (1.9) The other expansions are exprdt ≈1 + rdt u≈1 + σ√dt +1 2σ2dt d≈1−σ√dt +1 2σ2dt Here we have (u−1)2= (d−1)2=σ2dt. From which, we can write p(Vu−Vd) = p(u−d)∂V ∂S St(1.10) = (rdt +σ√dt −1 2σ2dt)∂V ∂S St.(1.11) 8 We will plug in (1.9) and (1.10) in (1.6) and cancel terms to produce V(1 + rdt) = rSt ∂V ∂S dt +V+1 2σ2S2 t ∂2V ∂S2dt +∂V ∂t dt or, rdt =rSt ∂V ∂S dt +1 2σ2S2 t ∂2V ∂S2dt +∂V ∂t dt. We cancel dt from this equation to obtain the Black Scholes PDE ∂V ∂t 1 2S2σ2∂2V ∂S2+rS ∂V ∂S −rV = 0. 1.8 Parity of Put-Call We will make a portfolio with one long Call option and one short Put option [65]: VPF (S, T) = VECO(S, T)−VEP O(S, T) So it becomes, VPF (S, T) = max(S−K, 0) −max(K−S, 0) = S−K. From the solution of Black-Scholes PDE with the final condition, we have, VPF (S, t) = S−Ker(T−t). Therefore we can say, VECO(S, t)−VEP O(S, t) = S−Ker(T−t). This is called the parity of Call-Put. As we can see from above: Call −Put =Asset −Forward. 1.9 Sensitivity Analysis Sensitivity analysis of the Option price with respect to some parameters known as Greeks [65]: 9 •Delta: ∆ is defined as: ∆ = ∂V ∂S . It is the measurement of the rate of change of the option price Vwith respect to the change in the underlying asset price S. Delta hedging use this Delta because the risk-less portfolio is stable in accordance to: QS QV =−∂V ∂S =−∆. where QV,QSare Options and Stocks numbers in the Portfolio. Now Delta for European Call and Put options are: ∆EC =∂V EC ∂S =N(d1). and ∆EP =∂V EP ∂S =−N(−d1). We can observe that: ∆EC ∈(0,1) and ∆EP ∈(−1,0) Evaluation of Delta for market data time series: We can figure out the implied volatility σimplied(t) from the market data series for the option price Vreal(t) and the underlying asset price Sreal(t). We have to solve: Vreal(t) = VEC (Sreal(t), t;σimplied(t)). •Gamma: It is defined as: Γ = ∂∆ ∂S =∂2V ∂S2 It is the measurement of the change of Delta of the option price Vwith respect to the change in the asset price S: ΓEC = ΓEP =∂∆EC ∂S =N0(d1)∂d1 ∂S =e−1 2d2 1 σp2π(T−t)S. •Rho: It is the sensitivity with respect to the interest rate r: ρ=∂V ∂r . 10 So it calculates the rate of change of the option price Vwith respect to the interest rate r. •Theta: It is defined as : Θ = ∂V ∂t . So it calculates the rate of change of the option price Vwith respect to the interest rate t. •Vega: Sensitivity with respect to the volatility σ: ν=∂V ∂σ . So it calculates the rate of change of the option price Vwith respect to the volatility σ. The Black-Scholes Equation is: ∂V ∂t +1 2σ2S2∂2V ∂S2+rS ∂V ∂S −rV = 0. The Greek Model for the Black-Scholes Equation is : Θ + σ2 2S2Γ + rS∆−rV = 0. Black-Scholes equation is important in the world of financial mathematics. In this chapter, we have given a description of the fundamental parts of Black-Scholes Equation. 11 Chapter 2 Linear Black-Scholes equation’s Analytic Solution 2.1 Introduction Black Scholes Equation is prominent in the quantitative finance world for the estimation of the option price. We discuss at length of the analytical solution of the Linear Equation in this chapter. Our main focus is to transform the Linear Black Scholes Equation into a heat-diffusion equation ; Then we solve the heat equation analytically. After that, we substitute the variables to get the solution of Black-Scholes PDE. 2.2 Black Scholes Equation’s Linear Model The Black Scholes partial differential equation [10] is ∂V ∂t +rS ∂V ∂S +1 2σ2S2∂2V ∂S2−rV = 0.(2.1) Here •V≡V(S, t) is the payoff function. It is the price of the option at a time t. It is a two variable function with the asset’s price Sand time t. •S(t) stock price which the underlying asset’s price at time t. It is also non negative. •tis the time . Here tis from 0 to T. ( Tis the expiration time.) 12 •ris the riskless interest rate. •σis the volatility. •Kis the strike price of an option. •µis the drift rate. In the case of Linear Black-Scholes Equation, σ, r, µ [65] are constants. We convert the linear Black-Scholes Equation to a heat equation. We need to have the initial and boundary conditions if we want to find the solution of Black-Scholes Equation. We are considering only the European Option Price. European Option Price: In this Option price, the domain of time tis, 0≤t≤T The domain of Stock price Sis, 0< S < ∞ European Call Option: From [10, 65]: V(S, T) = max(S−K, 0) for 0< S < ∞ V(0, t) = 0 for t ∈[0, T] V(S, t) = S−Ke−r(T−t)for S → ∞ European Put Option: From [10, 65]: V(S, T) = max(K−S, 0) for 0< S < ∞ V(0, t) = Ke−r(T−t)for t ∈[0, T] V(S, t)=0 for S → ∞ 13 Option at expiry time plotted over a range of stock prices are 70 ≤S≤130. We use the Black Scholes Equation’s exact solution to evaluate the value of the option price before expiration. Figure 2.1: Value of the European call option V(S, t) at time, t= 0, t=.2, t=.4, t=.6, t=.8 and t= 1. Figure 2.2: Value of the European call option V(S, t) at various times in a surface plot 20 Black-Scholes Equation is effective for various contexts of the financial field. In this chapter, we investigated the analytical solution of the Linear Black-Scholes model. And a graphical representation of the exact solution is given too. 21 Chapter 3 Linear Black-Scholes Equation with Spectral Method 3.1 Introduction We investigate Black-Scholes partial differential equation (PDE) with different Numerical schemes. In this chapter we solve it using spectral method. This method is well known for its superior efficiency and less time consumption. It takes less time than other numerical schemes such as finite-difference,finite elements and Gelarkin method. The spectral collocation method is applied for Black-Scholes PDE. We discuss briefly about the Spectral collocation method. There is a MATLAB execution of the Spectral Method [72]. We use the Chebyshev polynomial and Chebyshev points for the spectral method. Also there is a comparison of the numerical calculation with the exact solution at the end of this chapter. 3.2 Spectral Collocation Method Spectral Method has the ability to achieve the exponential accuracy to solve the differential equations and approximating the solution [67]. The accuracy of the Spectral Method is ruled by the smoothness of functions. We use polynomials to approximate a function [31]. The spectral collocation is also called pseudospectral method. Usually we use Chebyshev and Legendre polynomials. Let the weight be wiand collocation point be xi.These weights and collocation points connected to Chebyshev polynomial can be evaluated as [31]: 22 1. Chebyshev-Gauss-Lobatto: xi= cos πi N, w0=wN=π 2N, wi=π N 2. ChebyShev-Gauss-Radau: xi= cos 2πi 2N+1, w0=π 2N+1 , wi=2π 2N+1 3. ChebyShev-Gauss: xi= cos 2i+1 2N+2, wi=π N+1 3.3 Chebyshev Differentiation Matrix We use the Chebyshev-Gauss-Lobatto points. xj= cos πj N, j = 0,1,2,··· , N And Chebyshev collocation points are evaluated in the interval [−1,1] We use these points to build up the Chebyshev Differentiation Matrix. And this matrix is used as the differential operator. Now we are given a grid function vdefined on the Chebyshev points; we acquire the discrete derivative gin two stages [72]: •We let qbe the unique polynomial with a degree ≤Mwith q(xj) = vj,0≤j≤ M •gj=q0(xj) Therefore this operation has to be linear. We can represent this by the multiplication of (M+ 1) ×(M+ 1) matrix [72] , which is denoted by DM: g=DMv. (3.1) In here Mis an arbitrary positive integer. Now for M= 1. The interpolating points are: x0= 1, and x1=−1. 23 The Lagrange form for the interpolating polynomial throughout the v0and v1is, q(x) = 1 2(1 + x)v0+1 2(1 −x)v1. The derivative of q(x) is: q0(x) = 1 2v0−1 2v1.(3.2) Here M= 1, so the matrix is 2 ×2 matrix. Let’s denote aij = [aij] The entries of the matrix where i, j = 1,2,3,··· , M. For the first column of the matrix, we set x0= 1 and get: a11 =1 2, and a21 =1 2. For the second column of the matrix, we set x1=−1 and get: a12 =−1 2, and a22 =−1 2. Hence the Chebyshev Differentiation Matrix for M= 1 is: D1=1 2−1 2 1 2−1 2. Now for M= 2. The interpolating points are: x0= 1, x1= 0, and x2=−1. The Lagrange form for the interpolating polynomial throughout the v0,v1and v2is: q(x) = 1 2x(1 + x)v0+ (1 + x)(1 −x)v1+1 2x(1 −x)v1. The derivative of q(x) is: q0(x) = (x+1 2)v0−2xv1+ (x−1 2)v2.(3.3) 24 Here M= 2, so the matrix is 3 ×3 matrix. For the first column of the matrix, we set x0= 1 and get: a11 =x+1 2=3 2, For x1= 0, a21 =x+1 2=1 2. For x2=−1, a31 =x+1 2=−1 2. For the Second column of the matrix, we set x0= 1 and get: a12 =−2x=−2, For x1= 0, a22 =−2x= 0, For x2=−1, a32 =−2x= 2. For the Third column of the matrix, we set x0= 1 and get: a13 =x−1 2=1 2, For x1= 0, a23 =x−1 2=−1 2, For x2=−1, a33 =x−1 2=−3 2. Therefore the Chebyshev Differentiation Matrix for M= 2 is: D2=      3 2−21 2 1 20−1 2 −1 22−3 2       . We can give the generalization for the entries of DMfor arbitrary M[72, 30] 25 Theorem: Let’s consider for M≥1, we have the columns and rows (M+1)×(M+1) of the Chebyshev Spectral Differentiation Matrix DMfrom 0 to M. The entries of this matrix are: (DM)00 =2M2+ 1 6,(DM)MM =−2M2+ 1 6,(3.4) (DM)jj =−xj 2(1 −x2 j), j = 1,2,··· , M −1.(3.5) (DM)ij =ci cj (−1)i+j (xi−xj), , i 6=j, i, j = 1,2,··· , M −1.(3.6) where ci=(2i= 0 or M 1elsewhere. We can write this in Matrix form, DM=                        2M2+1 62(−1)j 1−xj··· ··· 2(−1)j 1−xj 1 2(−1)M −1 2 (−1)i 1−xi−xj 2(1−x2 j) (−1)i+j (xi−xj)··· (−1)i+j (xi−xj) 1 2 (−1)M+i 1+xi . . .(−1)i+j (xi−xj)−xj 2(1−x2 j) .... . .. . . . . .. . .......(−1)i+j (xi−xj) . . . −1 2 (−1)i 1−xi (−1)i+j (xi−xj)··· (−1)i+j (xi−xj)−xj 2(1−x2 j) 1 2 (−1)M+i 1+xi −1 2(−1)M−2(−1)M+j 1+xj··· ··· −2(−1)M+j 1+xj−2M2+1 6                        . The jth column of DMcarries the derivative of the degree Mpolynomial interpolant qj(x) to the delta function bolstered at xj, represented at the grid points xi If we have zero boundary conditions, then evaluation becomes easier. We can exclude the first and last column of DMmatrix because they will be multiplied by zero. 3.4 MATLAB function : ’cheb’ AMATLAB function called cheb [72] is used to evaluate the Chebyshev differentiation matrix DM. This function returns a vector xand a matrix D. The function 26 is given below: % CHEB Evaluate Diff = Chebyshev differentiation matrix, x = Chebyshev grid function [Diff,x] = cheb(M) if M==0, Diff=0; x=1; return, end x = cos(pi*(0:M)/M)’; c1 = [2; ones(M-1,1); 2].*(-1).^(0:M)’; X = repmat(x,1,M+1); dX = X-X’; Diff = (c1*(1./c1)’)./(dX+(eye(M+1))); Diff = Diff - diag(sum(Diff’)); This MATLAB function does not evaluate DMprecisely by the (3.4,3.5,3.6) formulas. Instead this makes the best use of (3.5) for the off-diagonal entries but then acquires the diagonal (3.4,3.5) from the identity: (DM)ii =− N X j=0 j6=i (DM)ij. From this, it is slightly easy to program. And it creates a matrix with better stability properties [59, 7]. We are giving a few Chebyshev Differentiaton Matrix DMmatrices evaluated by the MATLAB function cheb: >> cheb(1) ans = 0.5000 -0.5000 0.5000 -0.5000 >> cheb(2) 27 ans = 1.5000 -2.0000 0.5000 0.5000 0 -0.5000 -0.5000 2.0000 -1.5000 >> cheb(3) ans = 3.1667 -4.0000 1.3333 -0.5000 1.0000 -0.3333 -1.0000 0.3333 -0.3333 1.0000 0.3333 -1.0000 0.5000 -1.3333 4.0000 -3.1667 >> cheb(4) ans = 5.5000 -6.8284 2.0000 -1.1716 0.5000 1.7071 -0.7071 -1.4142 0.7071 -0.2929 -0.5000 1.4142 0 -1.4142 0.5000 0.2929 -0.7071 1.4142 0.7071 -1.7071 -0.5000 1.1716 -2.0000 6.8284 -5.5000 Here we can analyze that first two Matrices D1and D2matches with the entries of the matrices we evaluated analytically. Therefore this function cheb is accurate. 28 3.5 Numerical Solution of Black-Scholes Equation with the help of Spectral Method We are not going to solve the original Black-Scholes PDE. But we solve numerically the Heat Equation using Spectral Collocation Method. uxx =uτ. After the numerical solution of Heat Equation, the variables are back-substituted to get the original solution of Black-Scholes Equation. We can evaluate the second derivative by d2 dx2=D2 M. by squaring the Chebyshev Differentiation Matrix DM. From the equation (3.1) we can say: w=D2 Mv. Here D2 Mis an (M+ 1) ×(M+ 1) matrix which will map a vector (v0, v1,··· , vM)T to a vector (g0, g1,··· , gM)T. We can write this in matrix form:          g0 g1 . . . . . . gM−1 gM          =         D2 M                  v0 v1 . . . . . . vM−1 vM          . Spectral Method is used to discretize the spatial domain xin the interval x∈[−1,1]. We used the cheb function for the second derivative. The time domain τhas no extra calculation. We have a polynomial interpolant q(x). This polynomial has been calculated using the command polyval(polyfit(···)). The main algorithm of our code follows the algorithm from Program 34 of [72] where the Allen-Cahn equation (ut=uxx +u−u3) is solved using spectral method and Forward Euler. We can decrease the error exponentially, if we increase the degree of the polynomial. 29 1. The first order derivative: From [45, 53] ,we can express diagonal and off diagonal entries of the weighting coefficients of the first order. Diagonal entry: i=j, A(1) ij = N X k=1,k6=i 1 (xi−xk)=− N X k=1,k6=i A(1) ik . Off-diagonal entry: i6=j, A(1) ij = N Y k=1,k6=i (xi−xk) (xj−xi)Y k=1,k6=j (xj−xk) =1 (xj−xi) N Y k=1,k6=i (xi−xk) Y k=1,k6=j (xj−xk) =1 (xj−xi)· N Y k=1,k6=i,k6=j (xi−xk) (xj−xk). 2. Second order Derivative Diagonal entry: i=j, A(2) ij =− N X k=1,k6=i A(2) ik . Off-diagonal entry: i6=j, A(2) ij = 2A(1) ij A(1) ii −1 xi−xj. 3. Higher Order Derivatives: As we can see from the first two cases, there is a recurence formula for the higher derivatives. That is, to evaluate weighting coefficients of higher order , we need the weighting coefficients of all the lower orders relative to the one we need. Diagonal entry: i=j, A(l) ij =− N X k=1,k6=i A(l) ik . 36 Off-diagonal entry: i6=j, A(l) ij =l A(1) ij A(l−1) ii −A(l−1) ij (xi−xj)!. We have a Heat Equation to solve. We replace the second derivative in it by the Differential Quadrature Matrix. DQM = [Aij]. uxx =∂2u ∂x2=D2 QM u. 4.5 Matlab function : Differential Quadrature Matrix Mehmet Murat Altug Bicak coded the matlab function file for Differential Quadrature Matrix [9]. The function file is given below: function [D]=Diff_Quad(N) % Differential quadrature matrix for N points based on Lobatto grid points. for ii=1:N X(ii)=0.5*(1-cos((ii-1)*pi/(N-1))); end SagTaraf=1; Taraf=0; for ii=1:N for jj=1:N if ii~=jj a1=(1/(X(jj)-X(ii))); for k=1:N if k~=ii & k~=jj SagTaraf=SagTaraf*(X(ii)-X(k))/(X(jj)-X(k)); end end 37 a(ii,jj)=a1*SagTaraf; SagTaraf=1; end if ii==jj for k=1:N if k~=ii Taraf=Taraf+(1/(X(ii)-X(k))); end end a(ii,jj)=Taraf; Taraf=0; end end end D=a; tic A few outputs of Differential Quadrature Matrix produced by the function file are given below: >> Diff_Quad(1) ans = 0 >> Diff_Quad(2) ans = -1 1 -1 1 >> Diff_Quad(3) ans = -3.0000 4.0000 -1.0000 38 -1.0000 0.0000 1.0000 1.0000 -4.0000 3.0000 >> Diff_Quad(4) ans = -6.3333 8.0000 -2.6667 1.0000 -2.0000 0.6667 2.0000 -0.6667 0.6667 -2.0000 -0.6667 2.0000 -1.0000 2.6667 -8.0000 6.333 4.6 Numerically solving Black-Scholes PDE using Differential Quadrature Method Here, we use Differential Quadrature Method (DQM) to solve the transformed Equation (Heat equation) of Black-Scholes PDE. We evaluate the Differential Quadrature Matrix (DQM) with the help of [9] for Chebyshev-Gauss-Lobatto Grid points. This differential quadrature matrix converts the differential equation into a eigenvalue problem( y0=A∗y=⇒Dy =Ay). We solve the Black-Scholes Equation for the domain x∈[−1,1]. Now for our choice of Chebyshev-Gauss-Lobatto points, we get the grid points in the interval [0,1]. To convert the interval from [0,1] to [−1,1], we use 2x−1. We evaluate European Call Option V(S, t) for the following parameters: T= 1; r= 0.02, σ = 0.2; K= 100. 39 Figure 4.1: European Call Option Price V(S, t) using Differential Quadrature Method. Then we have the comparison with the exact solution of Black-Scholes Equation. And it matches with the exact solution: Figure 4.2: European Call Option Price V(S, t) is compared with the Exact solution. Now we have the table of errors compared to Differential Quadrature Method. We use infinity norm at t= 0. 40 Stock Price SDifferential Quadrature Method Exact Value Error 137.71278 47.37890 47.37287 0.00603 139.09681 48.74329 48.73885 0.00444 141.90675 50.12376 50.12078 0.00298 143.33294 52.93267 52.93226 0.00042 144.77346 54.36104 54.36175 0.00070 146.22846 55.80532 55.80704 0.00173 147.69808 57.26551 57.26816 0.00265 149.18247 58.74161 58.74510 0.00349 150.68178 60.23366 60.23791 0.00426 152.19616 61.74168 61.74663 0.00494 153.72575 63.26574 63.27130 0.00557 155.27072 64.80588 64.81200 0.00613 156.83122 66.36218 66.36881 0.00663 Table 4.1: Comparison of Differential Quadrature Method and Exact Solution for the European Call Option Price V(S, t). 4.7 Conclusion Introduction to Differential Quadrature Method is given in this chapter. We built up the differential quadrature matrix which helped extensively in our MATLAB implementation of Black-Scholes PDE. 41 Chapter 5 Nonlinear Black Scholes Equation and Volatility Models 5.1 Introduction In this chapter we discuss about Nonlinear Black Scholes Equation and several volatility models. A short and efficient analytical approach for the Nonlinear Black Scholes Equation is given too. We focus on the European Option price. 5.2 Nonlinear Black Scholes Equation The linear Black Scholes Equation is: Vt+1 2σ2S2VSS +rSVS−rV = 0.(5.1) Here S(t)>0 and t∈(0, T) There are some restrictive assumptions on the linear Black-Scholes Equation which do not match in reality. Because of transaction costs,[5, 11, 51],large investor preferences [27, 29, 64] and and incomplete markets[68]. In this thesis, we focus on couple of transaction cost models for the nonlinear Black–Scholes equations for European option price where µthe drift rate is constant and modified volatility function is not constant, ˜σ2= ˜σ2(t, S, VS, VSS). So the equation (5.1) becomes the nonlinear Black Scholes Equation: Vt+1 2˜σ2(t, S, VS, VSS)S2VSS +rSVS−rV = 0 (5.2) where S(t)>0 and t∈(0, T). 42 5.3 Boundary Conditions We have to finalize the differential problem by stating the terminal and bounadry conditions for European Call and European Put Option price. Only then we can find a unique solution for the equation (5.2). Although we only focus on European Call Option Price in this dissertation. 5.3.1 European Call Option European call option is the solution of (5.2) on 0 ≤S < ∞, 0 ≤t≤Twith the following initial condition and boundary conditions: V(S, T) = (S−K)+, for 0≤S < ∞, V(0, t) = 0, V(S, t)≈S−Ke−r(T−t), as S → ∞. 5.3.2 European Put Option The value V(S, t) of the European put option is the solution of 1.2 on 0 ≤S < ∞, 0 ≤t≤Twith the following initial condition and boundary conditions: V(S, T) = (K−S)+, for 0≤S < ∞, V(0, t) = Ke−r(T−t), V(S, t)≈0, as S → ∞. 5.4 Volatility Models In classical Black-Scholes model, the essential parameter which can not be observed but is assumed to be constant is known as the volatility σ. Many researchers tried different approaches to improve the model by using a modified volatility function ˜σ(.) to take transaction costs into account, big traders and illiquid markets which is eventually the cause of the nonlinearity of (5.2). In this part , we discuss briefly about couple of volatility models and focus will be on the transaction cost related volatility models. 43 1. We can replace the volatility σin the Black-Scholes model by the approximated volatility from the preceding values of the underlying asset. This volatility is called the historical volatility [32]. 2. We can calculate the modified volatility from the exact solution of the Classical Black-Scholes Model if the price of Option and every other parameters are known. The implied volatility is the value σ, for which call option solution and put option solution is true in comparison to the real market data. We can evaluate the implied volatility from the difference between the observed option price V(S, t) (from stock market) and the Black-Scholes exact solution where all the parameters are taken into account from the reality except for the implied volatility. 3. Another model is developed by Hull,White[38] and Heston[35], volatility follows the dynamics of a stochastic process which is called the stochastic volatility. 4. At first we have observed volatilities at each stock price and time which is termed as local volatility ˜σ= ˜σ(S, t). We can replace constant volatilities with them. Dupire [50] investigated the reliances and expressed the local volatility as a function of implicit volatilities. 5. We assume that every security is accessible at any size and any time or that any single trading will not impact the price, is not always correct. Hence,several authors have modeled illiquid markets and big trader effects.Frey and Stremme [23]and after some time Frey and Patie [28] examined these results on the price and emerge with the outcome ˜σ=σ 1−ρλ(S)SVSS .(5.3) where ρis constant,σthe historical volatility. These nonlinear models are all consistent with the linear model if the supplementary parameters for transaction cost vanish (Le, Ψ(.), M). There are couple of modified volatility function models, 44 5.5 Leland’s Model for Volatility Leland’s model is one of the infamous implied volatility model. By the trade at discrete times [51] was the plan of Leland’s to relax the hedging condtions. This also has the potential to reduce the Portfolio adjustment expenses. Leland assumes that the transaction cost is ω|∆|S/2. where ωis the whole trip transaction cost per unit money of the transaction. ∆ is the number of assets purchased (∆ >0) or sold (∆ <0) at price S, is proportional to the monetary value of the assets purchased or sold. We take a portfolio which replicates with ∆ units of the underlying asset and the bond B: Π=∆S+B. After that, a minor change in time of the size δt the change in the portfolio becomes δΠ=∆δS +rBδt −ω 2|δ∆|S. (5.4) Here δS is the change in price S. Therefore, the first term represents the change in value. Then we have the 2nd which represents the bond growth in δt time. After that we have δ∆ which represents the change in the number of assets. Hence, the last term becomes the transaction cost because of the portfolio change. We use Ito’s lemma for the value of the option price V:= V(S, t) and get: δV =VSδS + (Vt+σ2 2S2VSS)δt. (5.5) We assume that the option Vis replicated by the portfolio Π and their values have to be equivalent at all times and there can be no risk-less profit. Because of no-arbitrage argument, we get δΠ = δV Equating the terms in the last two equations ,we get ∆ = VS, 45 unprotected portfolio. In this way the portfolio is still well protected with the Risk adjusted Pricing Methodology(RAPM) and the improved volatility is, ˜σ2=σ2(1 + 3(C2M 2πSVSS)1/3). Here M≥0 is the transaction cost measure and C≥0 is the risk premium measure. 5.11 Analytical Solution We solve (5.2) subject to the terminal and boundary conditions . We change the variables [?, 24], x= ln( S K), τ=1 2σ2(T−t), u(x, τ) = e−xV(S, t) K. Differentiation gives, Vt=uττtS=−1 2σ2Suτ, VS=uxxSS+u=ux+u, VSS =uxxxS+uxxS=1 S(uxx +ux). Plugging these values into (5.2) ;we get, −1 2σ2Suτ+1 2˜σ2S(uxx +ux) + rS(ux+u)−ruS = 0. Now we are going to multiply by −2 Sσ2 uτ−˜σ2 σ2(uxx +ux)−Dux= 0. Where D=2r σ2and ˜σ2depends on the volatility model. x∈Rand 0 ≤τ≤σ2T 2Now we plug in the volatility models: ˜σ2=σ2(1 + Le.sign(uxx +ux)) or, ˜σ2=σ2(1 + Ψ(e2rτ σ2aKex(uxx +ux)) 52 or, ˜σ2=σ2(1 + (e2rτ σ2aKex(uxx +ux)) or, ˜σ2=σ2(1 + 3(C2M 2π(uxx +ux))1/3). 5.12 Conclusion In this chapter , we have discussed fairly about the different volatility models. We have given short descriptions of the famous volatility models. This chapter ended with the analytical solution of the Nonlinear Black-Scholes Model. This solution has similarities with the Linear one. 53 Chapter 6 Numerical Investigation of Nonlinear Black-Scholes Equation There are numerous numerical methods for solving Black–Scholes equation for European, American and Asian options. We can also solve linear Black-Scholes model using numerical methods. But there are major difference between the same numerical method as the volatility becomes a variable. We show the process and some different approaches to solve the volatility term. Also the volatility term is defined as different functions because of the nonlinearity. We have already showed in the previous chapters that there are exact solution for European call and put options. But if we seek more accurate answer which matches with the reality, we have to dug a bit deeper which eventually gives us the nonlinear Black-Scholes equations with various volatility models. For these complicated contracts in general settings, analytical formulae are seldom available and numerical methods have to used to solve the problem. These Methods are different from lattice methods which includes binomial and trinomial approximations [18], Monte Carlo methods using the least-square techniques [42], analytical approximations [1, 14, 52], finite-element discretizations [74, 47] to finite-difference methods [2, 13, 17]. We have another classical method which consists of developing of the free boundary problem into a linear complementary problem (LCP) and the solution by the Projected Successive Over Relaxation (PSOR) method of Cryer [20]. And there was another method developed called penalty and front–fixing [55]. The underlying model changes which is a disadvantage in these methods. 54 A complete different procedure which is based on a recursive calculation[37] of the early exercise boundary, estimating the boundary only at some points and then approximating the whole boundary by Richardson extrapolation. Explicit boundary tracking algorithms which are a finite-difference bisection scheme [49] or the front–tracking strategy of Han and Wu [33]. Our numerical investigation is solely focused on finite difference schemes. We give a detailed analysis on this scheme. 6.1 European Call option We are going to use finite-difference schemes to solve the converted partial differential equation (PDE) of Black-Scholes Equation, uτ−˜ σ2 σ2(uxx +ux)−Dux= 0.(6.1) Where x∈Rand 0 ≤τ≤˜ Twith the respective volatilities (5.2) subject to the initial and boundary conditions, u(x, 0) = (1 −e−x)+for x ∈R, u(x, τ)=0 as x → ∞, u(x, τ)≈1−e−Dτ−xas x → ∞. We are going to give a introduction to finite difference schemes and then present the numerical results from MATLAB. 6.2 Finite-Difference Methods The main goal of finite-difference schemes is to approximate the derivatives of the difference quotients and the solution of the emerging discrete schemes. At first, we start by discretizing the domain of the converted Partial differential Equation (PDE) of Black Scholes Equation with the respective volatilities. Then 55 we are going to replace the derivatives by appropriate difference quotients. Then our investigation using finite-difference continues in terms of convergence . We introduce classical finite-difference schemes for the European Call option. 6.3 Grid We are going to transform the spatial domain and the time domain. Now x∈Rand τ∈[0,˜ T] by a bounded interval x∈[−R, R], R > 0. Now we discretize the new computational domain by a uniform grid (xi, τn) with xi=ih and τn=nk, where h > 0 is the spatial step, and k > 0 is the time step. And i∈[−N, N],−R=−Nh, R =Nh, n ∈[0, M] and ˜ T=Mk. We denote the approximate solution of the partial diffential equation of BlackScholes in xiat time τnby Un i≈u(xi, τn). And discretize the initial and boundary conditions in the following way: U0 i= (1 −e−ih)+, Un −N= 0, Un N= 1 −e−Dnk−Nh. Now for a more fitting treatment of the unbounded spatial domain x∈R, we are going to introduce artificial boundary conditions [26] to confine the unbounded domain to a bounded computational domain. 6.4 Difference Quotients The spatial derivative can be approximated with forward differences: ux(x, τ) = u(x+h, τ)−u(x, τ) h+O(h). or backward differences: ux(x, τ) = u(x, τ)−u(x−h, τ) h+O(h). 56 We can sum these differences which results in central differences and we get: ux(x, τ) = u(x+h, τ)−u(x−h, τ) 2h+O(h2). For the second spatial derivative we compute using the Taylor formula: uxx(xi, τn) = u(x+h, τ)−2u(x, τ) + u(x−h, τ) h2+O(h2). Now we can callback that a function f(h) is in O(hr), if there exists a constant M > 0,such that |f(h)|≤ M|hr|as h→0 . It means that the quantity f(h) is bounded by a constant multiple of hrfor sufficiently small h[70]. We call the error between differential quotient and difference quotient the truncation error. To discretize the partial derivatives of the transformed Black-Scholes Equation, we are going to introduce the following notation for the forward difference quotient with the spatial step size h: D+ hUn i=Un i+1 −Un i h≈ux(xi, τn). and we leave out the error term O(h). Similarly the backward difference quotient with respect to the spatial variable is denoted as: D− hUn i=Un i−Un i−1 h≈ux(xi, τn). So the central difference quotient is: D0 hUn i=Un i+1 −Un i−1 2h≈ux(xi, τn). where truncation error of O(h) and O(h2) are omitted. Now for the second spatial derivative we introduce the standard difference quotient: D2 hUn i=Un i+1 −2Un i+Un i−1 h2≈uxx(xi, τn). Where the error term is O(h2).The central differences in the time variable are never used in practice because they always lead the way to bad numerical schemes, which are inherently not stable [76]. Now using the schemes results into a system of equations that can be written in a matrix form: AnUn+1 =BnUn+dn.(6.2) 57 where Un= (Un −N+1, ...., Un 0, ..., Un N−1)T∈R2N−1. An=            a0a10 0 ··· 0 a−1.........0 0............. . . 0............0 . . ..........a1 0 0 ··· 0a−1a0            ∈R(2N−1)×(2N−1). Bn=            b0b1b20··· 0 b−1.........0 b−2............. . . 0............b2 . . ..........b1 0 0 ··· b−2b−1b0            ∈R(2N−1)×(2N−1). dn=            b−2Un −N−1+b−1Un −N−a−1Un+1 −N b−2Un −N 0 . . . 0 b2Un N b1Un N+b2Un N+1 −a1Un+1 N            ∈R(2N−1)×(2N−1). Now the matrix Anis tridiagonal. Hence the resulting systems can be solved with linear effort O(N) using the Thomas algoritm [71, 33]. This is completed by first decomposing the matrix An=LnRninto a lower and an upper bidiagonal matrix and secondly solving LnRnUn+1 =BnUn+dnby forward and backward substitution. Hence , we solve LnYn=BnUn+dnfor Ynand after that we solve RnUn+1 =Yn for Un+1. We need Un −N−1and Un N+1 for the vector dn. But these values are not in the grid points that are taken into account. They are the points before and after the boundary 58 conditions. We define the auxiliary or ghost boundary conditions [70], Un −N−1= 0 and Un N+1 = 1 −e−Dnk−(N+1)h.(6.3) And after this, we are going to assume that, 1 X i=−1 ai= 2 X i=−2 bi= 1. This is satisfied by any consistent scheme after normalization of the coefficients [61]. 6.5 Volatility Models We can deal with the derivatives in the volatility in different ways. We can write the modified volatility as, ˜ σ2=σ2(1 + s(x, τ)). Here s(x, τ) denotes the volatility correction in xat time τ. And term depends on the first and second spatial derivatives of u. In the paper of During[25], he proposes a smoother approximation of uxx for the volatility correction by choosing: uxx(xi, τn)≈Un i+2 −2Un i+Un i−2 4h2=D2 2hUn i. Here we have truncation error of order O(h2). We deal with the nonlinearity explicitly in all the schemes. Then, from Leland’s Model: sn i=r2 π κ σ√δtsign(D2 2hUn i+D0 hUn i). Barles’ and Soner’s model’s volatility correction is : sn i= Ψ(eDτn+xia2K(D2 2hUn i+D0 hUn i)). Then volatility correction for the case when Ψ(.) is considered as the identity is: sn i=eDτn+xia2K(D2 2hUn i+D0 hUn i). 59 Lastly the volatility correction for the Risk Adjusted Pricing Mthedology(RAPM) is: sn i= 3 C2M 2π(D2 2hUn i+D0 hUn i)1/3 . There is a complication with the evaluation of sn iat the boundary, as in the theory we require Un∈R2N+3 to have the ability to calculate sn N−1and sn −N+1. But this calculation demands the need of Un −N−1and Un N+1, which are outside the domain. During’s paper [25] has the statement that the impact of the nonlinearity at the boundary is not that significant and it can be neglected for large R. So we can assume that our auxiliary or ghost boundary conditions are valid and say that: sn= (sn −N+1,··· , sn 0,··· , sn N−1)T∈R2N−1. 6.6 Existence and Convergence Our desire is to give a logical approximation for the sequence of the solutions of AnUn+1 =BnUn+dn. Then these needs to be satisfied: 1. A uniform solution Unhas to exist for each n∈[0, M −1]. 2. And Un ihas to converge towards the exact solution of uτ−˜ σ2 σ2(uxx +ux)−Dux= 0, as k→0 and h→0. We are going to callback the terms and conditions of the existence and convergence for the linear case. In which case the volatility correction sn iis equal to zero and Hence An=A,Bn=Band the coefficients of dnare constants. A uniform solution to the system of equations exists [41], if the matrix Ais regular(a regular matrix Ais described as a square matrix that for all positive integer n, is such that Anhas positive entries.). The scheme (2) converges, if its consistent and stable. Hence, we describe the terms consistency and stability. 60 •Consistency A scheme Lh,k of order (a, b) is consistent, if there exists a constant M > 0 such that, max i,n |Lh,kun i|≤ M(ka+hb) for sufficiently small h, k > 0. In here, un iis the exact solution to (1) in (xi, τn) and Lh,k is the finite difference scheme that excludes the truncation error of order O(ka+hb). •Stability To establish stability of (2), we are going to take into account computer evaluated vector including the rounding errors ˆ Un. So Aˆ Un+1 =Bˆ Un+dn+rn. where rndenotes the rounding errors. The error vector is, en=ˆ Un−Un, acts in according to Aen+1 =Ben+rn. To simplify this ,we are going to assume that e06= 0 which means that there is already a rounding error when evaluating the initial condition. Simultaneously, we are going to assume that the matrix vector multiplication to obtain Un+1 works precisely, So rn= 0 n∈[0, M −1]. We have the error: en+1 =A−1Ben= (A−1B)2en−1=··· = (A−1B)n+1e0. For us to have a stable system, preceding errors have to be damped and therefore we need (A−1B)n+1e0→0as n → ∞. 61 (a) RAPM model (M=.01, C = 30). (b) Barles’ and Soner’s model (a=.02). (c) ψ(x) = xselected as the Identity (a=.02). (d) Leland’s model (δt =.01, κ =.05). Figure 6.1: Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.0001 with the Crank-Nicolson Method. The impact of transaction costs which are modeled by the volatilities and evaluated with the Crank-Nicolson finite difference scheme can be seen in the next figure, we plot the difference between the European Call Option price with Vnonlinear(S, t)−Vlinear(S, t). transaction costs and the European Call Option price without the transaction costs. We get a significant price deviation between the classical (linear) Black-Scholes Model and the nonlinear model. 68 (a) ψ(x) = xselected as the Identity (a=.02). (b) Barles’ and Soner’s model (a=.02). (c) RAPM model (M=.01, C = 30). (d) Leland’s model (δt =.01, κ =.05). Figure 6.2: The impact of transaction costs Nonlinear Vnonlinear(S, t)−Vlinear(S, t) with the Crank-Nicolson Method. The difference that we have is not symmetric for all the transaction cost models. But decreases closer to the expiry date. And this is the expected result. We have plotted all the nonlinear models with the linear model in figure (6.3). Now for each volatility model and each difference scheme, we compare the deviation of accuracy of the above computation at t= 0. We use the infinity norm l∞. If we reduce the spatial step size to h=.001 or h=.0001, it improves the accuracy considerably. However, it increases the computational time monumentally. The deviation from the linear model is evaluated only for the stock price S(t) = 269 at t= 0. We can see the difference between the Models from the table. Volatility Model Deviation from the linear model sn i= 0 Identity .00068 Leland .02945 RAPM .03813 Table 6.1: l∞deviation for different volatility models for the stock price S(t) = 269 at t= 0. 69 Figure 6.3: Value of European Call Option price V(S, 0) for different volatility models vs the linear model with the Crank-Nicolson Method. The plot (6.3) shows the price of a European Call option for V(S, 0) that is at t= 0, and how all the transaction cost models converge. We have given the algorithms for explicit and implicit methods too. We use them to implement them into MATLAB too. At first we plot the European Call Option price V(S, t) for different volatility models in Implicit Method. Then we give a Nonlinear Vs Linear graph similar to the Crank-Nicolson in implicit method. We repeat this for explicit method too. They give the same resemblance with the Crank-Nicolson graphs. But CrankNicolson gives better numerical results. 70 (a) ψ(x) = xselected as the Identity (a=.02). (b) Barles’ and Soner’s model (a=.02). (c) RAPM model (M=.01, C = 30). (d) Leland’s model (δt =.01, κ =.05). Figure 6.4: Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.001 with the Implicit Method. Figure 6.5: The Identity model ψ(x) = xwith parameter (a=.02) for the (Vnonlinear(S, t)−Vlinear(S, t)) plot for spatial step h=.1 and time step k=.001 with the Implicit Method. 71 (a) ψ(x) = xselected as the Identity (a=.02). (b) Barles’ and Soner’s model (a=.02). (c) RAPM model (M=.01, C = 30). (d) Leland’s model (δt =.01, κ =.05). Figure 6.6: Value of European Call Option price V(S, t) of Nonlinear Black-Scholes models for various volatility models for spatial step h=.1 and time step k=.00001 with the Explicit Method. Figure 6.7: The Identity model ψ(x) = xwith parameter (a=.02) for the (Vnonlinear(S, t)−Vlinear(S, t)) plot for spatial step h=.1 and time step k=.001 with the Explicit Method. 72 6.11 Conclusion We have thoroughly investigated the Nonlinear Black-Scholes Equation numerically. At first we show how the Nonlinear Black-Scholes Equation can be represented as system of equations. We use matrcies instead of System of equations. Then we give graphical represenation of the solution in implicit, explicit and Crank-Nicolson Method. We give a Comparison between the different transaction cost modeled by volatilities. 73 Chapter 7 Fractional Black-Scholes Equation 7.1 Introduction In quantitative finance, option pricing is a major problem. There are models other than Black-Scholes model for pricing options which give better results; one such model is known as Finite Moment Log Stable (FMLS) which can be put in writing as fractional diffusion class equation [40]. Fractional Calculus is the calculus of the differentiation and integration where the order is a fractional number [56]. We seek to discover the analytical solution of Black-Scholes Equation’s Fractional Derivatives by the use of Adomian Decomposition Method. 7.2 Transformed Black-Scholes Equation Here, we recall from the second chapter where the analytical solution of Linear BlackScholes Equation was given. Let V(S, t) be the option price ,then ∂V ∂t +1 2σ2∂2V ∂S2+rS ∂V ∂S −rV = 0.(7.1) Having the initial and boundary conditions which are: V(0, τ)=0, V(S, τ)≈S, as S → ∞, V(S, T) = max(S−K, 0). 74 Here, we have Kis the strike price, σis the volatility, Tis the expiry time and ris the risk free interest rate. Then we apply changes to the variables: S=Kex, t=T−τ (1/2)σ2, V=Kv(x, t). Which gives us the ∂v ∂τ =∂2v ∂x2+ (k−1)∂v ∂x −kv, where k=2r σ2. Now The initial condition is: v(x, 0) = max(ex−1,0). 7.3 Fractional Derivative The most common fractional Derivatives are: Caputo: The fractional derivative of g(x) is defined as [58]: Dβg(x) = Kn−βDng(x) = 1 Γ(n−β)Zx 0 (x−t)n−β−1gn(t)dt. for n−1< β ≤n,n∈N,x > 0. Mittag-Leffler Function: This function Mβ(x) is defined as [54]: Mβ(x) = ∞ X j=0 xj Γ(βj + 1), where β > 0 and x∈C. 7.4 Fractional Black-Scholes Equation of Time Here , we asses the fractional Black-Scholes Equation for time [40]: ∂βv ∂τβ=∂2v ∂x2+ (k−1)∂v ∂x −kv, 0< β ≤1.(7.2) It has the initial condition v(x, 0) = ex−1, x > 0. Adomian Decomposition Method is used in the equation (7.2) to get: ∂v ∂τ =∂1−β ∂τ1−β∂2v ∂x2+ (k−1)∂v ∂x −kv. 75 We are going to integrate on the both sides with respect to τto get: v(x, τ) = v(x, 0) + Zτ 0∂1−β ∂s1−β∂2v(x, s) ∂x2+ (k−1)∂v(x, s) ∂x −kv(x, s)ds. We denote v(x, 0) = v0. Then v(x, τ) = Zτ 0∂1−β ∂s1−β∂2v0(x, s) ∂x2+ (k−1)∂v0(x, s) ∂x −kv0(x, s)ds. The repetitive method yields for m > 1 vm(x, τ) = Zτ 0∂1−β ∂s1−β∂2vm−1(x, s) ∂x2+ (k−1)∂vm−1(x, s) ∂x −kvm−1(x, s)ds. (7.3) Therefore from Adomian Decomposition Method, we have the solution of (7.2) by v(x, τ) = ∞ X m=0 vm(x, τ). We are going to use the initial condition: v(x, 0) = ex−1, x > 0. So, v0=ex−1. From (7.3) , we evaluate the rest of the terms of the solution. v1(x, τ) = Zτ 0∂1−β ∂s1−β∂2v0 ∂x2+ (k−1)∂v0 ∂x −kv0ds. =−ex−kτβ Γ(β+ 1) + (ex−1) −kτβ Γ((β+ 1) v2(x, τ) = Zτ 0∂1−β ∂s1−β∂2v1(x, 0) ∂x2+ (k−1)∂v1(x, 0) ∂x −kv1(x, 0)ds. =−ex(−kτβ)2 Γ(2β+ 1) + (ex−1) (−kτβ)2 Γ((2β+ 1) v3(x, τ) = Zτ 0∂1−β ∂s1−β∂2v2(x, 0) ∂x2+ (k−1)∂v2(x, 0) ∂x −kv2(x, 0)ds. =−ex(−kτβ)3 Γ(3β+ 1) + (ex−1) (−kτβ)3 Γ((3β+ 1). 76 So, the solution of the equation (7.2) by the Adomian Decomposition Method, v(x, τ) = ∞ X m=0 vm(x, τ) =ex(−kτβ Γ(β+ 1) −(−kτβ)2 Γ(2β+ 1) +(−kτβ)3 (Γ3β+ 1) −···) + (ex−1)(1 −−kτβ Γ(β+ 1) +(−kτβ)2 Γ(2β+ 1) +···) =ex(1 −Mβ(−kτβ)) + (ex−1)Mβ(−kτβ) =ex−Mβ(−kτβ). where Mβ(x) is called the Mittag-Leffler function of one parameter. v(x, τ) = ex− Mβ(−kτβ) is the exact solution of (7.2). And Mittag-Leffler function has the property [57, 46] when β→1, the solution (7.2) agrees with the exact solution of the BlackSchole PDE. 7.5 Fractional Black-Scholes Equation of Space Here, we consider the fractional Black-Scholes PDE for space [77] . In this cas, We have fractional derivative in space. This is in the sense of Riemann-Liouville. The Levy density for a Log Stable is given below: CLS(x) = (Bq |x|−1−βx < 0, Bpx−1−βx < 0. Here B > 0, p, q ∈[−1,1], p +q= 1 and 0 < β ≤2. We replace p= 0 and q= 1. Then this procedure converts into the well known Finite Moment Log Stable Method. The following fractional partial differential equation is given by (For European Option Price): ∂V (x, τ) ∂τ +r+1 2σβsec (βπ 2)∂V (x, τ) ∂x −1 2σβsec (βπ 2)∂βV(x, τ) ∂+xβ=rV (x, τ) (7.4) Here ∂βV(x,τ) ∂+xβis called the left Riemann-Liouville fractional derivative , ris risk-less interest rate and σis volatility. Now , whenever β→2, (7.4) matches with BlackScholes PDE. We solve this very equation using Adomian Decomposition Method. At first ,we set V(x, 0) = ex−1. 77 [22] Steven R. Dunbar. Stochastic processes and advanced mathematical finance. University of Nebraska-Lincoln, 2010. [23] B. Dupire. Pricing with a smile. Risk,, 7:18–20–71, 1994. [24] B. During. Black–scholes type equations: Mathematical analysis, parameter identification and numerical solution, 2005. [25] B. During, M. Fournier, and A. J ungel. High order compact finite difference schemes for a nonlinear black–scholes equation. Int. J. Appl. Theor. Finance,, 7:767–789, 2003. [26] M. Ehrhardt and R. E. Mickens. A fast, stable and accurate numerical method for the black–scholes equation of american options. International Journal of Theoretical and Applied Finance, 11:471–501, 2008. [27] R. Frey and P. Patie. Risk Management for Derivatives in Illiquid Markets:A Simulation Study. 2002. [28] R. Frey and A. Stremme. Market volatility and feedback effects from dynamic hedging. Math. Finance,, 4:351–374, 1997. [29] R. Frey and A. Stremme. Market volatility and feedback effects from dynamic hedging. Math. Finance, 4:351–374, 1997. [30] Daniel Henry Gottlieb, M. Yousuff Hussaini, and Steven A. Orszag. Theory and applications of spectral methods. 1984. [31] Philippe Grandclement. Introduction to spectral methods. Stellar fluid dynamics and numerical simulations: from the sun to neutron stars, 153, 2006. [32] M. Gunther and A. Jungel. Finanzderivate mit matlab – mathematische modellierung und numerische simulation. Vieweg, Wiesbaden, 2003. [33] H. Han and X. Wu. Fast numerical method for the black–scholes equation of american options. SIAM J. Numer. Anal.,, 41:2081–2095, 2003. 84 [34] Matthew J. Hancock. Infinite spatial domains and the fourier transform. 2006. [35] S. Heston. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Rev. Fin. Studies, 6:327–343, 1993. [36] S. Hodges and A. Neuberger. Optimal replication of contingent claims under transaction costs. Review of Futures Markets, 8:222–239, 1989. [37] J.-Z. Huang, M.G. Subrahmanyam, and G.G. Yu. Pricing and hedging american options:a recursive integration method. Rev. Fin. Studies.,, 9:277–300, 1996. [38] J. Hull and A. White. The pricing of options on assets with stochastic volatilities. J. Finance,, 42:281–300, 1987. [39] John C. Hull. Options, Futures and Other Derivatives. 2002. [40] H.Zhang, F.Liu, I.Turner, and Q.Yang. Numerical solution of the time fractional black–scholes model governing european options. Computers And Mathematics with Applications, 71:1772–1783, 2016. [41] E. Isaacson and H. B. Keller. Analysis of Numerical Methods. 1966. [42] P. Jackel. Monte Carlo methods in finance. 2002. [43] M. Jandacka and D. Sevcovic. On the risk–adjusted pricing– methodology–based valuation of vanilla options and explanation of the volatility smile. J. Appl. Math.,, 2005:235–258, 2005. [44] Zuzana Jankova. Drawbacks and limitations of black-scholes model for options pricing. Brno University of Technology., 2018. [45] Ram Jiwari, Sapna Pandit, and R C Mittal. A differential quadrature algorithm for the numerical solution of the second-order one dimensional hyperbolic telegraph equation. International Journal of Nonlinear Science, 13:259–266, 2012. [46] G. Jumarie. Laplace’s transform of fractional order via the mittag-leffler function and modified riemann–liouville derivative, 2009. 85 [47] A. Jungel. Das kleine Finite-Elemente-Skript. 2001. [48] Y. Kabanov and M. Safarian. On leland’s strategy of option pricing with transactions costs. Finance Stoch., 1, 1997. [49] T. Kim and Y. Kwon. Finite–difference bisection algorithms for free boundaries of american call and put options. J. KSIAM, pages 271–287, 2008. [50] Y. Kwok. Mathematical models of financial derivatives. 1998. [51] H. E. Leland. Option pricing and replication with transactions costs. J. Finance, 40:1283–1301, 1985. [52] L. MacMillan. Analytic approximation for the american put option. dv. Fut. Opt. Res., 1:119–139, 1986. [53] Gulnihal Meral. Differential quadrature solution of heatand mass-transfer equations. Applied Mathematical Modelling, 37,issue 6:4350–4359, 15 March 2013. [54] G.M. Mittag-Leffler. Sur lint´egrale de laplace–abel, 1902. [55] B. F. Nielsen, O. Skavhaug, and A. Tveito. Penalty and front–fixing methods for the numerical solution of american option problems. . Comp. Finance, 5:69–97, 2002. [56] N Ozdemir and M Yavuz. Numerical solution of fractional black-scholes equation by using the multivariate pade approximation. ACTA PHYSICA POLONICA A, 132:1050–1053, 2017. [57] Jigen Peng and Kexue Li. A note on property of the mittag-leffler function. Journal of Mathematical Analysis and Applications, 370:635–638, 2010. [58] I. Podlubny. Fractional differential equations. Academic Press,New York, 1999. [59] R.Baltensperger and J.-P.Berrut. The errors in calculating the pseudospectral differentiation matrices for chebyshev-gauss-lobatto points. Comp.Math.Appl, 37:41–48, 1999. 86 [60] R. Richtmyer and K. Morton. Difference Methods for Initial-Value Problems. 1967. [61] A. Rigal. High order difference schemes for unsteady one–dimensional diffusion–convection problems., volume 114. 1994. [62] Fabrice Douglas Rouah. Four derivations of the black scholes pde. 2013. [63] Jianan S. and Zhengyou Z. A mixture differential quadrature method for solving two-dimensional incompressible navier-stokes equations. Appl Math Mech, 20:1358–1366, 1999. [64] P. Schonbucher and P. Wilmott. The feedback–effect of hedging in illiquid markets. SIAM J. Appl. Math., 61:232–272, 2000. [65] Daniel Sevcovic. Analytical and numerical methods for pricing financial derivatives. Lectures at Masaryk University, 2011. [66] Chang Shu. Differential Quadrature and its Application in Engineering. 2000. [67] Chi-Wang Shu and Peter S. Wong. A note on the accuracy of spectral method applied to nonlinear conservation laws. Journal of Scientific Computing, 10:357–369, 1995. [68] H. M. Soner and N. Touz. Super replication under gamma constraints. SIAM J. Contr. Optim, 39:73–96, 2001. [69] J. Michael Steele. Stochastic calculus and financial applications. Springer, 274:74, 2001. [70] J. C. Strikwerda. Finite difference schemes and partial differential equations. 1989. [71] L. H. Thomas. Elliptic problems in linear difference equations over a network. 1949. [72] Lloyd N. Trefethen. Spectral Methods in MATLAB. 2000. 87 [73] Marek Uhliarik. Operator splitting methods and artificial boundary conditions for a nonlinear black-scholes equation, 2010. [74] M. G unther and A. Jungel. Finanzderivate mit matlab – mathematische modellierungund numerische simulation. Vieweg, Wiesbaden,, 2003. [75] W.Chen, C.Shu, W.He, and T.Zhong. The application of special matrix product to differential quadrature solution of geometrically nonlinear bending of orthotropic rectangular plates. Computers and Structures, 74,issue 1:65–76, January 2000. [76] P. Wilmott, S. Howison, and J. Dewynne. The Mathematics of Financial Derivatives. 1995. [77] Yang Xiaozhong, Wu Lifei, Sun Shuzhen, and Zhang Xue. A universal difference method for time-space fractional black-scholes equation. Advances in Difference Equations, 71, 2016. [78] T. Y., WuG., and R. Liu. A differential quadrature as a numerical method to solve differential equations. Computational Mechanics, 24:197–205, September,1999. 88