AN ITERATIVE METHOD FOR PRICING THE AMERICAN PUT OPTIONS WITH DIVIDENDS by Florica Coman B.Sc , Babes - Bolyai University, Romania, 1994 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES (Mathematics) THE UNIVERSITY OF BRITISH C O L U M B I A August 2005 © Florica Coman, 2005 Abstract In this thesis we present a fixed-front approach for pricing the American Put Op-tions with dividends, combined with an iterative method. This approach transforms the free boundary problem into a nonlinear parabolic differential equation defined on a fixed spatial domain. The resulting equation has been solved using an implicit numerical scheme. The method can be easily applied to Barrier Options with divi-dends and numerical results were shown to compare well with the optimal exercise boundary for different dividend yields and different values of the barrier. A l l the programming work has been done in M A T L A B and sample of codes can be found in Appendix B . i i Contents Abstract ii Contents iii Acknowledgements v 1 Introduction 1 1.1 Background Concepts 1 1.2 Literature Overview 2 1.3 The American Options and the Free Boundary Problem. Economic Intuition 6 2 Formulation of the Problem 9 2.1 The American Put with Dividends 9 2.2 Fixed domain transformation of the free boundary problem 10 2.3 Direct Discretization of P D E 13 2.4 Method Description 18 3 Numerical Results 21 3.1 Numerical Example 21 3.2 Convergence 28 i i i 3.2.1 The time step refinement 29 3.2.2 The time and spatial step refinements 34 3.3 The influence of the Dividend Yie ld on the Position of the free Boundary 35 4 The American Barrier Put with Dividends , 37 4.1 Preliminaries and Problem Formulation 37 4.2 Numerical solution of the Barrier 39 4.3 Errors 41 4.4 The Influence of the Barrier on the Position of the Free Boundary . 43 5 Conclusions 46 Bibliography 50 Appendix A 53 A . l Optimal exercise boundary near expiry 53 A.2 The Dimensionless Problem 54 Appendix B Matlab Codes 56 B.0.1 Optimal Exercise Boundary for American Put Option with Continuous Dividend Yie ld for D > r 56 B.0.2 The Change of Variable from x to S 60 B.0.3 Function ip and its Drivatives 63 B.0.4 Convergence for time step refinement 64 B.0.5 Convergence for the time and spatial step refinements . . . . 70 B.0.6 Optimal Exercise Boundary for American Put Barrier Option with Continuous Dividend Yie ld for D > r • • • 76 iv Acknowledgements This research was conducted under the supervision of Dr. Rachel Kuske. I would like to thank her for providing me with such an interesting problem as the topic of my thesis. I am also thankful to her for all the effort, patience, guidance, advice, and support throughout. I learned a lot from Dr. Anthony Peirce by taking two of his classes and doing project. I thank him for all the support he gave me, and for reading my thesis. I would also like to thank my friend Catherine Dupuis who helped me a lot through my early stage of confusion with Matlab. Final ly I wish to thank my husband, Catalin Alexandru, for his support and under-standing during my work over the past year. F L O R I C A C O M A N The University of British Columbia August 2005 Chapter 1 Introduction 1.1 Background Concepts A n option is a contract to buy or sell an underlying asset at a specified time, called expiry time, for a specified price called strike price. There are two types of options: c a l l options and put options. • The c a l l option is an option that gives the buyer the right, but not the obligation to purchase the underlying asset at a specified time and price from the writer of the call option. • The put option is an option that gives the buyer the right, but not the obligation to sell the underlying asset at a specified time and price to the writer of the put option. The underlying assets of an option can be commodity, stocks, stocks indices, foreign currencies, or future contracts. In terms of the exercise time there are two types of options: American and European options. European options can only be exercised at the expiry, while the American options can be exercised at any time up until expiry. American options are generally 1 more valuable than their European counterparts because of the greater flexibility given to the holder. A n exception is the American call options on assets without dividends, which wil l be discussed briefly in the next section. The valuation of American options either with or without dividends has been an interesting and continuous research topic for many authors for many years, going back to Merton (1973) up to today. One of the most popular approaches used to model the pricing of options is the Black - Scholes equation. This equation has been used to value both the European and American options with or without div-idends . For the European options this equation is solved as a boundary value problem whose analytic solution can easily be found [22]. For the American options the Black - Scholes approach leads to a free boundary value problem. There are two ways to solve this type of problem: analytical and numerical. Due to the fact that the price of American options depends on the history of the underlying asset price as well as its present value, there is no exact solution for the free boundary [22]. 1.2 Literature Overview To our knowledge, a number of papers present analytical approximations for Amer-ican options [7], [11] , [14], [15]. Some of these results provide complicated solutions which may be impractical, while some of them give approximated expressions only near expiry . On the other hand, various numerical methods have been developed to obtain numer-ical values of the optimal exercise boundary. Binomial methods are one of the most commonly used approaches to solve this problem, but in order to obtain accurate results, considerable computational effort is necessary. Cox, Ross and Rubinstein 2 (see [3]) developed a standard binomial model by improving the fluctuating behavior of the optimal exercise boundary through an appropriate interpolation technique, and Breen ( see [2]) used an "accelerated" binomial model by using Richardson's extrapolation technique. More efficient numerical methods than the binomial meth-ods are those using trinomial models [12], but all of the above procedures have some issues regarding the speed of computation or accuracy. Another type of numerical method widely used, is the finite difference approach implementing either explicit or implicit schemes. For example in ([10]) Ikonen and Toivanen used the operator splitting method to solve the complementarity problem, or [21] in which Sevcovic used the Fourier integral transformation method. Some other numerical approaches besides the ones mentioned above, are the methods of lines used by Meyer in [17], the moving index method to solve the linear complementarity problem [13] by find-ing the approximate position of the moving boundary which corresponds to the optimal exercise boundary of the option, or the path integral technique which uses Fourier-Hermite series expansions [4]. The academic literature presents a lot of papers that try to solve the problem of pricing barrier options as well. Analytical solutions have been found for the price of a continuously monitored down-and-out European call [16], standard barrier E u -ropean options [20], and also analytical formulas have been found for some other exotic options such as: various types of double barrier options, partial barrier op-tions and rainbow barrier options. The analytical methods work only for simple barrier options such as European-style options with an assumed lognormal stock price evolution. Due to their nature, the American barrier options are very difficult or almost impossible to solve analytically. Gao, Huang, and Subrahmanyam (see [6]) found a quasi-analytic expression for 3 American barrier options using the decomposition technique, by separating the E u -ropean option value from the early exercise premium. These results have been found under the assumption that the underlying asset price follows geometric Brownian motion. As we have seen, based on a literature search, in most of the cases there are no general analytical solutions for the value of the American barrier options, therefore numerical methods are required. The most popular numerical methods to price an American barrier options are those using binomial and trinomial tree. The meth-ods using binomial trees converge very slowly as the number of the binomial levels increases, and can lead to big errors even if very small time steps are used. Another problem with these methods is that the barrier has to lie on the horizontal layer of the nodes in the tree which is very difficult to control. More advantageous meth-ods known are those using trinomial tree [19], or "adaptive mesh" [5]. Some other methods used to price the American barrier options are Monte Carlo methods of integral evaluation [1], but the numerical accuracy of all these methods becomes an important issue. However, because of the importance of American options with or without barrier, and because not all methods are practical, the research into efficient numerical meth-ods for valuation of these options and optimal exercise boundary, continues. A n interesting approach to solve the problem of valuation of American options is that suggested firstly by Landau in 1950. He proposed a certain change of variables in order to transform the free boundary problem into a nonlinear parabolic differen-tial equation defined on a fixed spatial domain. Based on this idea, two interesting methods have been developed to remove the problem of the free and moving bound-ary; the penalty method [18] and the front-fixing method [18] and [23]. 4 Nielsen, Skavhaug, and Tveito in [18] proposed, for the front-fixing method, a new spatial variable x such that x = where S is the stock price and S/(t) is the optimal exercise boundary at time t. The resulting parabolic differential equation has been solved using two methods: an implicit-upwind and an explicit-upwind dif-ference schemes. It has been concluded that the results provided by both schemes are almost identical. The second method, which is proposed in this paper is the penalty method, which transforms the free boundary problem into a problem posed on a fix spatial do-main, by adding a penalty term to the differential equation. The resulting equation has been solved numerically using explicit, semi-implicit, and fully implicit schemes. Due to the fact that there is no analytical solution of the above problem, the approx-imate solution obtained from the front-fixing method has been used as a reference solution in order to study the performance of the penalty method. Wu and Kwok in [23] presented a front-fixing method to solve the same problem by considering a different spatial variable given by y = lngjjTj where 5.is the stock price and Sj(t) is the optimal exercise boundary at time t. The resulting equa-tion has been solved numerically using a semi-implicit scheme in order to solve a tridiagonal linear system rather than a nonlinear system of equations. The results obtained from the front-fixing method were compared with those obtained by the binomial method concluding that the two methods have similar levels of accuracy. When the asset price is very close to the optimal exercise boundary it seems that the front-fixing method is "slightly more accurate". Both of the above papers, [18] and [23], study the case of standard American put options, specifying that the discussed methods can be applied for any "standard American options, barrier options, Asian options, and lookback options" as long as 5 a proper change of variable exists, but none of these situations were discussed. 1.3 The American Options and the Free Boundary Prob-lem. Economic Intuition A n American style derivative is a contract whose cash flows can be influenced by the holder of the derivative. The holder affects the cash flows of the contract through an exercise strategy. The optimal exercise strategy is the one which wil l provide the holder of the option with maximum value. Of course, we don't know ahead of time what the optimal strategy is. In practice, the optimal strategy is found simultaneously with the price. It is this optimal strategy that is associated with the concept of the free boundaries. Therefore we can view the American option as a free boundary value problem. Let us consider S be the asset price, io = 0 be the current time, t = Tp be the option maturity time and K the strike price. There must be some value of S for which it is optimal from the holder's point of view to exercise the American option otherwise, we should hold the option for all possible S. Then this option is identical to a European option. But we have seen in [22] that this is not the case. In other words, there is a 5/(t), such that if S < Sj(t), one should exercise the put option, which maximizes the payoff function K — S and for S > Sf(t) , we should hold the option. In contrast, if S < Sf(t) one should hold the call option and exercise it when S > Sf(t) in order to maximize its payoff. This S/(t) is referred as the optimal exercise price. In other words, the put option value at asset price S and time t should be given by P(S,t) = max(K — 5,0) for S < Sf(t), and P(S,t) satisfies the Black-Scholes equation for S > Sf(t). Similarly 6 the call option value should be P(S,t) = max(S - K,0) for S > S/(t), and P(S,t) satisfies the Black-Scholes equation for S < Sf(t) (see [22]). However, we do not know S/(t) a priori. We should treat Sj(t) as a new unknown called a free boundary. We claim that there are some conditions on Sf(t) such that 1. P{S,t) is continuous across (Sf(t),t), 2. d g^g'^ is also continuous across (Sf(t),t). We will use this framework to analyze the put case, since the analysis of the call case is very similar. Case I: Without Dividends. When the underlying asset pays no dividends, it is never optimal to exercise the American call early (see [8]). If the call is exercised early, the holder must immediately pay K. From that time forward, if the underlying asset price S rises or remains the same, the net benefit to the holder would be no greater than if he had held the option. However if the asset price falls, the holder will now take a loss. Therefore, early exercise locks in the time value of strike price K and incurs a downside price risk while introducing no new benefits. Because exercise would only occur at expiry, the American and European values for the call should be the same in this case. In contrast, it may be optimal to exercise the American put early. If S is sufficiently low, it may be more worthwhile for the holder to exercise and invest K at the risk-free rate r than to wait for further declines in S that will increase its payoff [8]. Case II. With Dividends. Many stocks pay out dividends which are in fact "payments to shareholders out of the profit made by the company" [22]. When the underlying asset pays dividends, it may be worthwhile to exercise the American call prior to the ex-dividend date, since after the dividend payment the option becomes less attractive. In this case, the dividend income from time to onward may outweigh 7 the expected return to the holder from a rise in S from to to Tp. Similarly, early exercise of the put may be optimal in the presence of dividends. The payment of dividends generally causes a downward movement in S , all else equal, and this only enhances the incentive for early exercise [22]. Lastly, we remark that two common models of dividend structure are: 1. Discrete dividend payments 2. Continuous dividend yield In the first case, the discrete dividend payments will cause sharp downward jumps in S immediately prior to the ex-dividend dates. In the second case, the effect of the dividends is smoothed out over time, but the qualitative intuition regarding the optimality of the early exercise remains the same. In this project, we will assume a continuous dividend yield. 8 \ Chapter 2 Formulation of the Problem 2.1 The American Pu t wi th Dividends For a constant dividend yield D, the Black-Scholes analysis of the option pricing problem for the American put can be described as a free boundary value problem where the value of the put P{S, t) at asset price S and time t and the optimal exercise boundary S — Sf(t) at time t solve the following backward-parabolic system [9]: r ) P 1 r)2P f)P ^ + r2s2ht + ir-D)sh~rP=°> i n o < * < 7 > , s>sfm-i) P = K-S}, |f = - 1 - a t 5 = 5/(t) (2.2) P ~ 0, as 5 -> +oo (2.3) I K if D < r , P(S,7» = max (K - S, 0) , Sf(Tp) = { (2.4) |^ —A" if D > r . Here K > 0 is the strike price, Tp is the expiry time of the option, r is the risk-free rate, and a > 0 is the volatility of the underlying asset. The option satisfies the early exercise condition (2.2) where the first condition is used to determine the option value on the free boundary, while the second condition is used to determine the location of the free boundary. Equation (2.3) represents the boundary condition for the option: as the asset price increases without bound it becomes unlikely that the option will be exercised and so P ~ 0. The payoff at maturity (t = Tp) determines the option price at any time prior expiry. Thus, the option value known at expiry is given by the final condition in (2.4). Also, there is an intuition that, for D < r the free boundary starts at Sf(Tp) = K, whereas Sf(Tp) = %K in the case D > r [9]. 2.2 Fixed domain transformation of the free boundary problem In this section we will perform a fixed domain transformation of the free boundary problem (2.1) - (2.4) into a parabolic equation defined on a fixed spatial domain. To transform equation (2.1) that holds on a time dependent spatial domain (S/(t), +oo), we define new variables x and T in place of S and t by introducing a change of co-ordinates: t = T (2.5) x = S - 0(T)ip(S) (2.6) 1 if S < K~ , ^ 5 ) = ) \-\zrf{c-(S -(K- +K+)/2) if K~ < S < K+ , (2.7) 0 if 5" > where c is a general coefficient that can be adjusted. K~ and K+ are considered some values expressed as a function of the strike price K and they can also be adjusted. (p(S) is a continuous function whose graph is presented below. 10 4.5 5.5 6.5 S($) Figure 2.1: Graph of <p(S) vs. S. If 5 < K~ we have x — S — 8(T) which means that all the points less than K~ are shifted to the left with S(T) while for S > K+ we have x = S which means that all the points greater than K+ remain the same. For any points between K~ and K+ the transformation is given by (2.6) for a choice of c to keep <p(S) continuous and approximately smooth as shown in Figure 2.1. We note that /3(Tp) = Sf(Tp) and discuss our treatment of 0 later in (2.49) and When 5 = 5 / we conclude from equations (2.6) - (2.7) that <p{Sf) = 1 and that x = Sf — 8(TF)ip(Sf) = 0. Therefore the time dependent spatial domain (Sj(t), +oo) has been transformed into a fixed spatial domain (0,-(-oo). Also, it is clear that (2.50). Te[0,TF}. Using the change of variables (2.6), we compute the partial derivatives: (2.8) 11 32P ±d_P d_dP( 8,, as2 ds'os' ds'dx v ^ v 'as £<£>('-«^K(-< )^ a2P 9x 2 (2.9) dP 8P8T dPdx dP dP ( , m 3 0 \ -at = df-dt+^m = dr + ^ {-^df) ( 2 - 1 0 ) B y substituting these partial derivatives into (2.1) we get the modified form of the Black-Scholes partial differential equation: }2 D / a,„\ 2 This equation is equivalent to: dP_ dx dP 1 2 2 / 1 2 2 5V , D)S{,-m%)-m^rP (2 ,2 ) Hence the previous equation becomes: d P d 2 P 8P - § = ^ S ^ + b ^ - r P ( 2 - 1 3 ) a(T,5) = ^ 2 5 2 ( l - / 3 ( T ) | | ) 2 (2.14) b(T, S) = - I a 2 5 2 ^ ( T ) | ^ + (r - D) S (l - / J ( T ) | | ) - <p{S)% (2-15) 12 The free boundary condition (2.2) has to be changed using the new variable (x) as well. Thus, equation for S = Sf gives ff = Then, the free boundary condition corresponding to the new variable is dP — = - 1 , at a; = 0 (2.17) ox It also follows from the boundary condition (2.3) that P ~ 0 , a s x - > + o o (2.18) To summarize, the backward - parabolic system (2.1) - (2.4) becomes after the previous change of variables: dP d2P dP ~dT = a { T > S ) d x ^ + b { T ' S ) f a ~ r P > ™0<t<TFy x>0 (2.19) P = K - S f , ^ = -1, atx = 0 (2.20) J ox P ~ 0 , &sx->+oo (2.21) [ K i f l > < r , P(S, TF) = max(K - 5,0) , Sf(Tp) = { (2.22) [ —K HD>r. where the coefficients a(T,S) and b(T,S) are given by (2.14) and (2.15) 2.3 Direct Discretization of P D E The Black-Scholes P D E (2.19) has three terms with continuous partial derivatives of which, two are first order and one is a second order partial derivative. Approxi-mations to the first and second order spatial partial derivatives wil l be derived using 13 central differences. It is well known, that the backward difference approximation for leads to an implicit method. The Black-Scholes equation is: dP d2P dP backward diff. central diff. central diff. Consider the discretization of the spatial variable (x) into M equally spaced units of dx. Given these (x) values, we find the values of (S) that implicit ly solve (2.5) and (2.6) and represent them by (S). Also, for the time variable (T) consider the discretization into N equally spaced units of dt. The discretization of (x) and (T) can then be written as: x = (0, dx, 2dx, 3dx,(M — l)dx) = (x\, x2, £ 3 , XM) (2.24) T = (TF-dt,TF-2dt,TF-3dt,...,0) = t2,h, ...,tN) (2.25) Then, we have in general that xm = (m — l)dx and xm+i = rndx from (2.24) and that t„ = TF- ndt and tn+1 = TF- (n+ l)dt from (2.25) Also the discretization for (S) is: S = (SUS2,S3,...,SM) (2.26) Therefore we consider the following approximations for the existent partial deriva-tives: dP{x,T) P{x,T - dt) - P{x,T) dT -dt dP{x, T) _ P{x + dx, T) - P(x - dx, T) dx 2dx 14 (2.27) (2.28) d2P(x,T) P[x + dx,T) - 2P(x,T) + P(x - dx,T) dx2 " dx2 [ y j Substituting the approximations as indicated in (2.27) - (2.29) into (2.23) and letting P ( ( m - l)dx, TF - ndt) = P(xm, tn) = P™ gives: pn _ p n + 1 p n + 1 _ n p n + 1 , pn+1 pn+1 _ pn+1 r m r m _ „ n + l m+1 m ^ J m - 1 , m+1 •* m - 1 _ p n + 1 /Q o n \ dt ~ a m dx2 + ° m 2dx m [ 6 V ) Thus (2.30) equation becomes: P- = P^_\(-a^ A + e 1 ^ ) + ^ ( 1 + reft + a ^ 1 ^ ) + pn+1, n + l _ ^ _ _ i n + l J j ^ N (n S11 ^ m + l l « m ^ "rn 2 d x ) The boundary condition (2.18) is given by P ^ + \ = 0 for n = 0,1, 2 , T V - 1 which leads to the equation: PS, = PMM^T^ + b n M + 1 ^ ) + ^ M + 1 ( 1 + rdt + < > | | ) (2.32) The free boundary condition (2.17) can be written as follows by using the central difference to discretize the first order partial derivative pn+1 pn+1 2dx (2.33) Then P 0"+! = P2+1 + 2dx (2.34) A t x — 0 (i.e. m = 1) the equation (2.31) is: Substituting the free boundary condition as indicated in (2.34) into (2.35) gives: ir-t-^'^-tr'^P.*) 15 Then - 2 a " + 1 - ^ + 6" + 1 r i t (2.37) Therefore equation (2.37) becomes: pn + cn+l = p n + l i l + r d t + a n + l ^ ) + pn+l(_2an+l_^) ( 2 3 g ) Where c ? + 1 = 2 a ? + 1 f - b^+1dt. Now for n = 1 , 2 , 3 , N , we assign the following initial values to the M-length vectors A, B and P: / l - i ^ ( l - f l , | | ) 2 (2.39) B - ^ A . g + <r - D) S ( l - fl,g) - » ( 2 4 0 ) P <- m a ^ i f - 5 , 0 ) r (2.41) where K — S = (K — Si, K — S2, K — SM) a n d 0 represents the M-dimensional row vector. We now assign the following values to A4± and M.2 the M x M matrices: diag(A) = lMAT (2.42) diag(B) = IMBT (2.43) Mi = -2diag{A) + diag(A)I$ + diag(A)Ij}, Mi{l, 2) = 2 a ? + 1 (2.44) M2 = -diag{B)lj} + diag{B)I+}, M2(l, 2) = 0 (2.45) 16 Where I ^ 1 and I M X are defined as shifted identity matrices with Is shifted up by one row in I ^ 1 and down one row in IT^1: I + 1 -0 1 0 0 ... 0 0 0 0 1 0 ... 0 0 0 0 0 1 ... 0 0 0 0 0 0 ... 0 0 0 0 0 0 ... 0 1 \ 0 0 0 0 ... 0 0 / I - 1 M 0 0 0 0 . . 0 0 1 0 0 0 . . 0 0 0 1 0 0 . . 0 0 0 0 1 0 . . 0 0 0 0 0 0 \ 0 0 0 0 0 0 1 0 / We have: Mi -2a?+1 _ n + l 0 0 0 0 n + l 2a -24+1 ,51+1 , n + l -2a n + l . .n+l 0 0 0 0 ar1 - 2 < + 1 -2nn+1 n + l aM 0 0 0 0 , n + l M / -2an+1 1 0 0 0 0 . .. 0 0 \ -bn2+1 0 bn2+l 0 . . 0 0 0 0 hn+l . 0 0 M2 = 0 0 0 . . 0 0 0 0 0 0 . 0 hn+l \ 0 0 0 0 . hn+l • °M 0 J 17 If we introduce A4, the main M x M matrix: M = -dt\Mx{^) + M2{^) - rIM] + IM (2.46) and define C to be an M x 1 column vector where C ( l ) = — + ^ a " + 1 and all other elements are 0, we can rewrite the discrete approximation (2.31) to the P D E with boundary condition (2.32) and free boundary condition (2.38) in matrix form: MPn+1 =Pn + C; (2.47) Where ( pn \ pn r2 pn Pn+1 = / pn+l\ pn+1 r2 pn+1 ^ 3 pn+1 rM-l \ Pn+1 I C (2.48) J Pn being the option value at time tn = Tp — ndt and Pn+i being the option value at time £ n +i = Tp — (n + l)dt. 2.4 Method Description This numerical procedure is an iterative method that uses the near-expiry expres-sions of the optimal exercise boundary ( A . l ) - (A.3) as the init ial guess for the global values of the optimal exercise boundary. This method involves one main iteration and one subiteration. 18 M a i n Iteration. The numerical procedure involves a z - step main iteration and an N - step subiteration within each main iteration step. A t the start of the main iteration, we first define 0{T) which is the first guess of the global value of the optimal exercise boundary. For fl = 1 , 2 , z , we assign: 0(Q) <- (Sf(TF-),Sf(t1),S,{t2),...,Sf(tn)), fi = l • (2.49) P(il) - 0(j - 1), = j > 1 (2.50) where S/(Tp) is defined in (2.4) and S/(t) for any t < Tp are defined in ( A . l ) -(A.3). We wil l refer to the i — th element of 0 as 0{. Next , we discretize the derivative of 6 with respect to t, ^ = 0' as follows: for the first n — 1 components we use the central difference and for the last component we use the forward difference. Recall that we solve the partial differential equation backward in time from t — Tp to t = 0. Thus, „, , 0 1 - 0 3 02-04 Pn-l ~ 0n+l 0n ~ Pn+1, , R , ai a , s. , „ P ={-2dr'-2dT> 2dt ' dt ) = ( ^ ( 3 2 , - , 0 n ) (2-51) Subiteration. Under each main iteration step, we have an N — steps subiteration that gives the option value at each time from t = Tp to t = 0. In other words, we find the option value using an implicit method backward in time. W i t h the above definitions of 0 and 0', we begin the subiteration by first discretizing the variable x as in (2.24). Given these x values, and the values of 0 from (2.49), we find the values of S = (SI,S2,S3,...,SM) (2.52) by solving implicitly the equation (2.6). For this problem we have used the M A T L A B code called x2s.m (see Appendix B ) . 19 For each subiteration step, we also assign: C^Cdt (2.53) where C has been defined in (2.48) . After that, we solve the system given in (2.47) by assigning P M~l{P + C) (2.54) A t the end of the subiteration procedure, we obtain an 1 x N row vector P = ( P f = \ P f = 2 , . " , - P i n = A r ) w h e r e pi=k f o r k = 1,2,...,N is the first element in the vector P on the k — th step of the subiteration, and it represents the option value for x = 0 at time t = tn=k. Therefore, P = ( P ^ 1 , P f = 2 , p « = w ) is the value of the option at S — Sf from t = Tp to t = 0. The first free boundary condition in (2.2) gives us an update of the free boundary which wil l be used in the next step of the main iteration. Thus, we assign: (3(n + 1) <- K - (Pr\Pr2,-,PTN) (2-55) This process wil l be repeated for each step in the main iteration and terminates at n = z. 20 Chapter 3 Numerical Results 3.1 Numerical Example In this section we focus on numerical results. We wil l compute the free boundary profile for three different situations: D < r, D = r, and D > r. Time steps in terms of t are the units of the vertical axis and asset price S are the units of the horizontal axis. The truncation of the spatial domain has been done at x = 15 for each of the following numerical example. Case: D < r In what follows we consider the following set of parameters: the volatility a = 0.25, the risk-free rate r = 0.1, the dividend yield D = 0.05, the maturity time Tp = 1, the strike price K = 5, the number of time steps n = 20, and the number of spatial steps m — 60. In this case we insert the time values of (2.25) into ( A . l ) . In other words , we start with the initial guess for the global values of the optimal exercise boundary given by ( A . l ) . 21 For sufficiently small values of D with respect to r, the original expression of Sf(t) will produce imaginary values away from expiry. In order to avoid the errors caused by these imaginary values we wil l use the following expression for time away from expiry: Xf(j) = x0 - 2 • AQ{T) (3.1) This is an extension of the approximation for the optimal exercise boundary, derived from the dimensionless problem, to values of t away from expiry ([9]). We invert the change of variables (see Appendix A , (A.8) and (A.10)), and convert this to: Sf(t) = K e W - T f - t ^ a ^ (3.2) In this situation (see Figure 3.1), the convergence of the procedure is very quick. Notice that only after a few iterations the procedure is converging. 22 ti towards solution- progroi towirdi solution - ti Figure 3.1: On the left, we show the graph of 5/ as a function of 5 and t, for r = 0.1, D = 0.05, a = 0.25, n = 20, and m = 60. Progress towards solution of the optimal exercise boundary for 4 iterations. The stars indicate the approximation at step fl = 4 and the solid line indicates the approximation at the previous step. On the right, we show the graph of 5/ as a function of S and t for the same choice of parameters when the number of iterations is fi = 8. Case: D = r In this case (see Figure 3.2) we consider the following set of parameters : the volatil-ity a = 0.25, the risk-free rate r = 0.1, the dividend yield D = 0.1, the maturity time Tp = 1, the strike price K = 5, the number of time steps n = 20, and the number of spatial steps m = 60. The initial guess for the global values of the optimal exercise boundary is given by (A.2) 23 3.6 38 4 4,2 44 4.6 4.8 5 35 4 45 5 Si l l SII) Figure 3.2: On the left, we show the graph of Sf as a function of S and t, for r = 0.1, Z) = 0.1, (T = 0.25, n = 20, and m = 60. Progress towards solution of the optimal exercise boundary for 4 iterations. The stars indicate the approximation at step Q = 4 and the solid line indicates the approximation at the previous step. On the right, we show the graph of Sj as a function of S and t for the same choice of parameters when the number of iterations is Q = 8. Case: D > r In this case (see Figure 3.3) we consider the following set of parameters: the volatility a = 0.25, the risk-free rate r = 0.1, the dividend yield D = 0.15, the maturity time Tp = 1, the strike price K = 5, the number of time steps n = 20, and the number of spatial steps m = 60. The initial guess for the global values of the optimal exercise boundary is given by (A.3) 24 pfogtoti towards solution - oxortiB boundary pfojrsi toward* lolulion - exercte bound»ry Figure 3.3: On the left, we show the graph of Sf as a function of S and t, for r = 0.1, D = 0.15, <T = 0.25, n = 20, and m = 60. Progress towards solution of the optimal exercise boundary for 4 iterations. The stars indicate the approximation at step CI = 4 and the solid line indicates the approximation at the previous step. On the right, we show the graph of 5/ as a function of S and t for the same choice of parameters when the number of iterations is CI = 8. Next, we consider the case when D = r and zoom in the graph (Figure 3.4) very close to expiry for the first three iterations and the near-expiry expression of the optimal exercise boundary given in (A.2). 25 0.9995 0.999 J3 0.9985 0.998 0.9975 0.997 4.7 4.75 4.8 4.85 4.9 4.95 5 SIS) Figure 3.4: Zoomed graph of 5/ as a function of S and t of the first three iterations and the near-expiry analytical expression (the dotted curve) of the optimal exercise boundary when r = 0.1, D = 0.1, n = 20, m = 60. A close look at this graph suggests us that the numerical approximation of the optimal exercise boundary doesn't follow the same shape as the near-expiry analytical expression, the former being sharper than the last one. This is the cause of too large time step near expiry and the influence of the boundary condition imposed in the problem. The same situation will occur in the other two cases D < r, D > r, as well. In order to solve this problem, we have considered the same time step away from expiry as before, and we decreased linearly the time step as we get closer to the time near expiry. Hence a nonuniform time step has been used in order to obtain an improved optimal exercise boundary. Therefore, the time interval from 0 to 1 has been divided into 20 equal intervals, i.e. [0, 0.05, 0.1, 0.15,..., 0.9, 0.95, 1]. Hence the time step for this grid is dtfog= 0.05. After that, we refined the time step for the time interval near expiry, [0.95, 1], which has been divided into 10 equal intervals, i.e [0.95,0.955, 0.96, 0.965,..., 0.995, progress towards solution - exercise boundary 26 1]. Hence the time step corresponding to this grid is dtsmau— 0.005. Then we set a linear function for dt: dt = 2 • (dtsmaii - dtbig) -t+(2- dtbig - dtsmau) (3.3) with dtbig= 0.05 being the time step away from expiry and dtsmaa = 0.005 being the time step near expiry t = Tp-In Figure 3.5 and Figure 3.6 we wil l present the comparisons between the opti-mal exercise boundary previously found and the improved optimal exercise boundary obtained using the nonuniform time step given by (3.3) for each case of the divi-dends. The same parameters as before have been used to make these comparisons. sin S I S i Figure 3.5: On the left, we show the comparison of the optimal exercise boundary for an uniform time step (the curve with circles) and the improved optimal exercise boundary for nonuniform time step (the curve with stars) when r = 0.1, D = 0.05, n = 20, m = 60, and = 25 iterations. On the right, we show comparison of the optimal exercise boundary for an uniform time step (the curve with circles) and the improved optimal exercise boundary for nonuniform time step (the curve with stars) when r = 0.1, D = 0.1, n = 20, m = 60and fi = 25 iterations. 27 0.9 • 0.8 -07 • 0.6 -10.5 -0.4 -0.3 -0.2 -0.1 -oL 2.85 2.9 2.95 3 3.05 3.1 3.15 3.2 3.25 3.3 3.35 S(S) Figure 3.6: Comparison of the optimal exercise boundary for uniform time step (the curve with circles) and the improved optimal exercise boundary for nonuniform time step (the curve with stars) when r = 0.1, D = 0.15, n = 20, m = 60, and f2 = 25 iterations. 3.2 Convergence The numerical scheme is said to converge if the numerical solution tends to the exact solution as the discretization is refined i.e. as the space and time steps tend to zero. From the beginning we have mentioned that the exact analytical expression for the global free boundary profile is not known yet. Therefore, we have to look for an "exact" solution by considering very small time and spatial steps. For that, we considered the following set of parameters: m = 60 • 2 3 which gives us the spatial step dx « 0.03, and n = 20 • 2 7 which provides a time step away from expiry dt^g ~ 0.0004. In what follows we will be using the improved solution of the optimal exercise boundary. We have also considered about fl= 40 iterations in order to obtain the benchmark of the "exact" solution. 28 3.2.1 The time step refinement We want to study the convergence of our method by studying the behavior of the optimal exercise boundary as the time step tends to zero and the spatial step is kept constant at dx = = 0.25. Nonuniform time steps have been used over the entire domain in order to obtain the improved solution for each time step refinement. The time step refinements were achieved by successively halving the time step away from the expiry , starting from dtbig — 25- I n other words, we have used the following time step values away from the expiry: 1 1 1 1 b i 9 E { 2 0 ' 2 0 _ 2 ' 20~22' 20~¥} We performed a sequence of runs to evaluate the optimal exercise boundary for each of the time steps previously considered and after that we compared these results with the "exact" optimal exercise boundary obtained for dx ss 0.03, dtbig ~ 0.0004 and Cl= 40 iterations . We have treated separately each possible case for the dividends, i.e. D < r (Figure 3.7), D — r (Figure 3.8), and D > r (Figure 3.9). 29 Figure 3.7: On the left, we show the graphs of the optimal exercise boundary for different time steps vs. t, when r = 0.1, D = 0.05, n = 20, rn = 120, and fi = 25, compared with the "exact" optimal exercise boundary (the curve with stars) obtained for dx « 0.03, dtbiy « 0.0004 arid fi= 40. On the right, we show the zoomed graph near expiry of the optimal exercise boundary for different time steps vs. t for the same set of parameters as for the first figure 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 S($) Figure 3.8: Graphs of the optimal exercise boundary for different time steps vs. t, for r = 0.1, D = 0.1, a = 0.25, n = 20, rn = 120, and Q, = 25 iterations, compared with the "exact" optimal exercise boundary (the curve with stars) obtained for dx w 0.03, dt^g « 0.0004 and 0= 40. 30 dt = 0.05 dt = 0.0125 dt = 0.00078125| - « — exact Figure 3.9: Graph of the optimal exercise boundary for different time steps vs. t, for r = 0.1, D = 0.15, a = 0.25, n = 20, m = 120, and fl = 25 iterations, compared with the "exact" optimal exercise boundary (the curve with stars) obtained for dx « 0.03, dtbig ~ 0.0004 and ft= 40. Notice that the optimal exercise boundary for different time steps may not approach the "exact" solution at all values of t since dx is fixed (see Figure 3.7). In what follows we will study the errors between the " exact" solution and the optimal exercise boundary obtained as a result of the decrease in the time step. The error could be calculated using the infinite norm function, but this error will be based only on the region near t = Tp, since the maximum error is always near t = Tp (see Figure 3.11), so it does not give us a clear idea about the convergence of the method. Since we need to see the convergence over the whole curve, we have to involve all points and one option for that would be to use the average error over the entire domain: Eil'i \S) - Sfact\ error = (3.4) where Sj- is the optimal exercise boundary corresponding to the time step dt = —T^T, gexact j g "exact" solution, and 7Vt is the number of data points. This error in-31 volves all the data points which gives us a better idea about the convergence over the whole curve not only near t = Tp. In the case of D > r (see Figure 3.10) the figure shows clearly that the error is decreasing as the time step tends to zero. Therefore, the optimal exercise boundary tends to the "exact" solution as dt —> 0. In the other two cases when D < r and D = r (see Figure 3.10) the error decreases up to a certain point and after that it starts increasing. Figure 3.10: Graph of the average error for different time steps , for r = 0.1, a = 0.25, D = 0.05 (the diamond curve), D = 0.1 (the star curve), D = 0.15 (the circle curve), and Q. = 25 iterations 32 If we consider, for example, the case when D < r we observe that the error is decreasing up to the second time step refinement, i.e. dt = a n c ^ f ° r the n e x * three refinements it is increasing. We have plotted the error vector |5y — Sjxact\ vs. time for these last three refinements in order to check where this increase in the error is coming from (see Figure 3.11). • dl = 0.00625 dl = 0.003125 • dl = 0.0015625 0.91 0.92 0.93 0.95 0.96 t (years) 0.97 0.98 0.99 Figure 3.11: Zoomed plot of the error vectors vs. t for different time steps: dt = 2 0 1 2 3 (the dashed curve), ^^(the dotted curve), and 5^5-(the dashdot curve) for r = 0.1, a = 0.25, D = 0.05,and fi = 25 iterations Near t = Tp we can observe on the above figure that there is a maximum which becomes larger as the time step decreases. This explains why the error in-creases for these situations. Similarly, if we solve the problem for fixed dt while refining the spatial step, the error does not decrease monotonically, since Sf is not well resolved for t near Tp. 33 3.2.2 The time and spatial step refinements In order to achieve a reasonable convergence due to the above issues, we wil l consider the time step refinement of the form: dt = and the spatial step refinement of the form : dx — • ^ where i = 0,1,2, 3,4, 5. In Figure 3.12 we present the plot of the errors for each case ( D < r, D = r, and D > r ) when both, the spatial step, and time step have been refined. We observe a continuously decrease in the error as we refine both steps which means that the error has been improved. -3 , 1 1 1 , , 1 1 . . 1 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 log2(N) Figure 3.12: Graph of the errors for American Put Option with Dividends when the spatial step is dx = | | • ^ py and the time step is dt = where i = 0,1,2,3,4,5, for r = 0.1, a = 0.25, D = 0.05 (the diamond curve), D = 0.1 (the star curve), D = 0.15 (the circle curve), and CI = 25 iterations 34 3.3 The influence of the Dividend Yield on the Position of the free Boundary The results of our method show that the optimal exercise boundary and the option price depend on the dividends yield D, the interest rate r, and the volatility a [9]. To confirm this observation, we firstly plot in Figure 3.13 the "exact" optimal exercise boundary Sefxact vs. t for each particular case when D < r, D = r, and D > r . Figure 3.13: The changes in the "exact" Optimal Exercise Boundary for American Put Option with Dividends for r = 0.1, a = 0.25, and different values of the dividends yield: D = 0.05 (the dotted curve), D = 0.1 (the dashdot curve), D — 0.15 (the solid curve) From the above graph we can observe that the optimal exercise boundary is pushed to the left as the dividend yield D increases. The financial explanation for this observation is that, in general, the dividend payments are causing a decrease in the value of the stock and implicity an increase in the option price. For each partic-ular time we have the value of Sf(t) given by P(Sf(t), t) — max(K — Sf(t), 0) which represents the intersection point between the option price curve and the payoff. For 35 large dividend yield we wil l have a large increase in the option price, so the price curve wi l l be pushed up (see Figure 3.14) and the intersection point between this curve and the payoff wil l provide a small value of the optimal exercise boundary Sf(t) at that time. Since the payoff function wi l l have high values for any S < Sf(t) the option holder wil l immediately exercise the option. In Figure 3.14 we show the comparison of the option price for different times when the dividend yield is D < r and D > r. Figure 3.14: The graph of the option price for r = 0.1, a = 0.25, D = 0.05 (on the left) D = 0.15 (on the right) for different times. In both cases the option price has been graphed at the same times. 36 Chapter 4 The American Barrier Put with Dividends 4.1 Preliminaries and Problem Formulation In this section we study the problem of the American barrier put option with dividends. Barrier options are the most common exotic options. They are path-dependent options, with payoffs that depend on the price of the underlying asset at expiration and whether or not the asset price crosses a barrier during the life of the option. There are two categories of Barrier options: "knock-in" and "knock-out". a. Knock-in options become worthless unless the asset price hits the barrier before the expiration date. b. Knock-out options become worthless if the asset price hits the barrier before the expiration date. More explicitly, there are four types of put or call single barrier options [22]: • up-and-in: the option comes into existence if S > X before expiry. • down-and-in: the option comes into existence if S < X before expiry. 37 • up-and-out: the option ceases to exist if S > X before expiry. • down-and-out: the option ceases to exist if S < X before expiry. Depending on the location of the barrier related to the strike price the option can be in-the-money (if the strike price K is lower than the barrier X for calls and higher than the barrier X for puts) or out-of-the-money (if the strike price K is higher than the barrier X for calls and lower than the barrier X for puts). Addit ional features of some barrier options are the payment of a rebate when the option reached the barrier, or the requirement to pay an additional premium when the option established. Some other barrier options can be made increasingly com-plicated by allowing multiple barriers, e.g. double-knock-out. In what follows we wil l consider American barrier put with dividends but with zero rebate. For these options the only cases of practical interest are the up-and-out barrier options. The down-and-out cases reduce to a standard American put option if the barrier is in-the-money ( see [6]), while if it is out-of-the-money the barrier option wi l l expire worthless. The value, P(S,t) of an American put up-and-out barrier option on an asset with constant dividend yield D, is a function of time t and current asset price S and it satisfies the following conditions: d2P + (r-D)S—-rP 0 in 0 < t < TF , as2 Sf{t) <S<X (4.1) 1 &tS = Sf{t) (4.2) P 0 at 5 = X (4.3) 38 The barrier is located at 5 = X with K < X considering that the option is out-of-the-money. The Black - Scholes P D E holds in (4.1) if the asset remains above the barrier X. However as soon as the asset price hits the barrier the option wi l l expire worthless. This provides the boundary condition (4.3). A t expiry the option must have the put payoff given by (4.4). Therefore, from the above problem, we notice that if the asset price stays in some predefined region, given by the barrier, the payoff at the maturity wi l l coincide with the payoff of vanilla American put with dividends. If the asset at some point in time exist from the region the option becomes and remains worthless. Such a con-struction wil l make the barrier option cheaper than the vanilla option, although the possible gains at maturity of both options remain identical. 4.2 N u m e r i c a l so lu t ion of the B a r r i e r The approach of finding a numerical solution of the above problem, is similar with the one used for the American put options with dividends. We consider the same change of variables as in (2.5) - (2.6). There are a few modifications that have to be done for the barrier options. In this case the Black - Scholes P D E (4.1) holds on a finite domain whose right endpoint is defined by the barrier X. This means that our barrier coincides with one of the spatial mesh points. Because the spatial domain becomes smaller and smaller as the barrier gets closer to the strike price K, we have to make sure that the barrier X is bigger than K+ in order to be able to apply the same transformation as the one used for the American put option (see x 2 s . m in Appendix B ) . So, in this case we wil l want to have a different choice of K+ and also have to adjust the value of 39 constant c involved in the definition'of function </?(5). We have changed the choice for K+ such that it is slightly less than X: where AX is an integer multiple of dx. We know that the definition of <p(S) concentrates the spatial points in the region where the graph of the option price P(S, t) is curved. As X is getting close to K, then also K+ and K~ are approaching each other, so that tp{S) may be sharply sloped and its derivatives may give large coefficients in the equation for P(S,t). So, in order to have a good approximation for the optimal exercise boundary and option price P we have to chose large enough values for c in the definition of <p(S) in (2.7) when the difference K+ — K~ is not very big so that <p(S) is continuous and smooth. One more thing that has to be considered in finding the optimal exercise boundary in the barrier case is to keep the boundary value as zero always (4.3). Hence through the given time interval, we keep checking the spatial variable. When the spatial variable is equal to the barrier level, zero is assigned for the option price, while the option expires worthless whenever the asset price hits the barrier (4.3). It was found that as X approaches K (i.e. X < 5.3) we need to modify the code to allow for ip(S) which varies with t, in order to handle the changing size of the domain as barrier moves to the left. We do not treat this case here, as it adds additional complication. The optimal exercise boundary for the barrier option has been found using a uniform spatial step dx = ^ where X is the barrier and m is the number of spatial points, and a nonuniform time step given by (3.3) where dtsmau is the time step near expiry and dtbig is the time step away from expiry. We notice that the iterative method (4.5) 40 for the barrier put option converges as fast as the procedure for the put option ( see Figure 4.1). To emphasis this conclusion, we consider an example for the following set of parameters: the volatility a = 0.25, the risk-free rate r = 0.1, the dividend yield D = 0.05, the maturity time Tp = 1, the strike price K = 5, the number of time steps n = 20, the number of spatial steps m = 120, and the barrier X = 5.5. projre* toward* tolutian - txtrau boundary ptogrsit toward* tolutsn - siordte boundary Figure 4.1: Progress towards solution of the optimal exercise boundary for 4 iterations for American Put Option with Dividends (the left figure) and American Barrier Put Option with Dividends (right figure)when for a = 0.25, r = 0.1, D = 0.05 and barrier level X = 5.5 4.3 E r r o r s In what follows we wil l consider the value of the barrier X = 5.5 and as in the previ-ous chapter, we have to find the "exact" solution of the optimal exercise boundary in the barrier case. For that, we consider very small spatial and time steps. We take m — 60 • 2 3 which gives us the spatial step dx w 0.01, and n — 20 • 2 7 which provides a time step away from expiry dtbig ~ 0.0004. We have also considered about tl = 40 41 iterations to obtain this benchmark. Moreover we want to check the convergence of this procedure for the barrier option. As we have seen in the case of American Put with dividends, we have to decrease the spatial step and the time step in the same time in order to improve the error. We consider the time step refinement of the form: dt = and the spatial step refinement of the form : dx = ^ • ^ where i — 0 ,1 ,2 ,3 ,4 ,5 . After performing a sequence of runs from i — 0 to i = 5 we can calculate the error between the "exact" solution and the optimal exercise boundary obtained for each run. We used the average error as in the previous chapter in order to show the convergence of the method over the whole curve: ENt I ci qexact | i = l \ ° f I . . error = — (4.6) where Sxj is the optimal exercise boundary corresponding to the time step dt = 20^, and spatial step dx — ^ • Sjxact is the "exact" solution, and Nt is the number of data points. In Figure 4.2 we observe a continuously decrease in the error for each particular case (D < r, D = r, and D > r ). 42 '0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 log2(N) Figure 4.2: Graph of the errors for Barrier American Put Option with Dividends with the barrier at X = 5.5 when the spatial step is dx = ~ • J^J and the time step is dt = where i = 0,1,2,3,4,5, for r = 0.1, a = 0.25, D = 0.05 (the diamond curve), D = 0.1 (the star curve), D = 0.15 (the circle curve), and H = 25 iterations. 4.4 T h e Inf luence of the B a r r i e r on the P o s i t i o n of the Free B o u n d a r y In this section we want to illustrate the influence of the barrier on the position of the free boundary for the American up-and-out barrier option with dividends. For that, we plot the optimal exercise boundary for different levels of the barrier (Figure 4.3). We choose three different values for the barrier when K < X, specifically: X — 5.4, X — 5.8, and X = 6.2. For the other parameters involved in the model we wil l use the following values: Tp — 1, r = 0.1, and o — 0.25. There wil l be considered two cases when D < r and D > r. The "exact" optimal exercise boundary for each considered value of the barrier was found using a spatial step dx = • ^ and a time step away from expiry dt = 2Q\>,. 43 Figure 4.3: The graph of the "exact" Optimal exercise boundary for American Put Option (with no barrier, i.e. X —> oo) and the "exact" Optimal exercise boundary for American Barrier Put Option with Dividends for different barrier levels, for r = 0.1, a = 0.25, D = 0.05 (the left graph) and D = 0.15 (the right graph) As shown in the above graphs, we extract useful information regarding the barrier option. We know that the asset price for the barrier options has to satisfy an additional condition for the option holder to receive the payoff. For an up-and-out option this condition is given by the fact that the asset price S is less than the bar-rier X in which case the holder owns a standard option. If the barrier is touched, the option dies worthless. As the barrier level decreases the probability that the option wil l die worthless increases and this wil l cause a decrease in the option price. For small values of the barrier we wil l have small option price, so the price curve wil l be pushed down (see Figure 4.4) and the intersection point between this curve and the payoff wi l l provide a big value of the optimal exercise boundary Sj. So, for any S < Sj one should exercise the put option, which maximizes the payoff function K - S and for S > Sf one should hold the put option. The above graphs confirm 44 our economic intuition that as X is increasing the optimal exercise boundary for the barrier option is getting closer and closer to the optimal exercise boundary for the standard option (X —> oo) and therefore the former one becomes more valuable as X gets larger. Another important observation is that the optimal exercise boundary does not change much as the barrier decreases when D > r, since in the first case, when D < r, the changes are more significant. As we have mentioned in the previous chapter, dividends have the effect of reducing the value of the optimal exercise boundary Sf(t). Therefore, the increase in both the dividends yield D and the barrier level X wil l cause a downward movement in 5 / and this only enhances the incentive for early exercise. In Figure 4.4 we can notice the difference in the option price P(S,t) for different values of the barrier X. Figure 4.4: The graph of the barrier option price for r = 0.1, a = 0.25, D = 0.05, X = 5.4 (on the left) and X = 6.2 (on the right) for different times. In both cases the barrier option price has been graphed at the same times. 45 Chapter 5 Conclusions In this thesis we have studied a new method to find the numerical solution of the Black-Scholes equation in order to valuate American options with continuous divi-dends yield. This method is a front-fixing method which transforms the free bound-ary problem into a nonlinear parabolic differential equation defined on a fixed spatial domain. This has been done by using a change of variables involving a smooth and continuous function ip(S), that allows the placement of a large number of spatial points on the curved portion of the option price P(S, t). The resulting equation has been solved using an implicit numerical scheme. The big advantage of this method, besides the fact that the equation is solved on the fixed domain, is that it finds simultaneously the value of the optimal exercise boundary and the option price without other extra computation. Another impor-tant advantage of this procedure, is that it can be easily applied on the barrier case when the infinite domain for the American put option can be transformed into a finite one without any truncation. We have treated separately the American put option and the American barrier put option in the presence of continuous dividend yield. 46 There are a few important factors that have to be considered when we discuss about any numerical method: the accuracy and convergence, the computation speed , flex-ibility, and the error estimates. In terms of accuracy and speed we present below a table which compares these results for American put options with dividends with and without barrier. For that we used a spatial step dx « 0.08 and a time step away from expiry dt « 0.05. These values are relatively large compared to the values used for the "exact" solutions shown in Figure 3.12 and Figure 4.2. Table 1: The number of iterations to reach a specified average error between suc-cessive iterations for n = 20 and m = 120. A P D is the abbreviation for American Put with Dividends and A B P D is the abbreviation for American Barrier Put with dividends. D < r D = r D > r A P D A B P D A P D A B P D A P D A B P D 1.00e-3 7 7 8 8 4 4 1.00e-4 16 16 22 20 17 16 1.00e-5 27 27 40 39 53 51 The convergence and the errors of this method has been studied in Chapter 3 for the option without barrier and Chapter 4 for the option with barrier. Regarding the flexibility, we have seen that the front-fixing method developed ini-tially for the American put with dividends can be easily used and implemented for the American put option without dividends and for the American barrier option with or without dividends. As we have mentioned in the first chapter there are two papers in the academic liter-ature, that have been deeply studied in order to make a comparison of our method's 47 results and their corresponding results; [18] and [23]. Both papers present the val-uation of American options without dividends, using fixed-front methods, but none of them is focusing on the barrier case. Therefore, we have compared our numerical results for American put option when the dividend yield D = 0 with the numerical results obtained by Nielsen, Skavhaug, and Tveito in [18]. The parameters used for this comparison were: r = 0.1 , a — 0.2, K = 1, and Tp = 1 and the required average error between successive iterations was error = l.OOe — 05. The truncation of the domain has been done at x m a x = 3 Table 2: The comparison of the optimal exercise boundary (OEB) at time t = 0 from our method and the front-fixing method developed by Nielsen, Skavhaug, and Tveito in [18]. O u r O E B T h e O E B deve loped i n [18] dtbig « 0.0008 and dx « 0.003 dt = dx = 0.001 Sf = 0.8627 Sf = 0.8616 Our method converged after 17 iterations to reach the above optimal exercise boundary value at time t = 0 and the average error between successive iterations was error = l.OOe — 05. For comparisons made in the barrier case (Table 3) we used online derivative pricer from [25] based on the trinomial tree model. We have used the following set of parameters: r = 0.1, a — 0.25, K = 5, Tp = 1, D — 0.05 for the results in the first table and D = 0.15 for the result in the second table, and different values of the barrier. For the stock price we used the optimal exercise boundary Sf at time t = 0 corresponding to each "exact" solution. The results are presented in the table below: 48 Table 3: The comparison of the option price P(Sf,Q) for the barrier case at time t = 0 corresponding to the O E B from our front-fixing method and the trinomial tree method [25]. In the first table we have the results for the case D < r and in the second table we have the results for the case D > r D < r barrier O E B at t = 0 Front-Fixing Method Trinomial Tree X = 5.4 X = 5.8 X - 6.2 3.9982107 3.8654037 3.8346735 1.00210 1.13570 1.16533 1.00179 1.13460 1.16225 D > r barrier O E B at t = 0 Front-Fixing Method Trinomial Tree X - 5.4 X = 5.8 X = 6.2 2.8511479 2.8419684 2.8406980 2.15030 2.15970 2.16100 2.14885 2.15803 2.15930 49 Bibliography [1] P.P. Boyle. Options: A monte carlo approach. J. Fin. Econ, 4:323-338, 1977. [2] R. Breen. The accelerated binomial option pricing model. J. Fin. Quan. Anal., 26:153-164, 1991. [3] J .C . Cox S.A. Ross and M . Rubinstein. Option pricing: A simplified approach. J. Fin. Econ, 7:229-263, 1979. [4] C. Chiarella N.El-Hassan and A.Kucera . Evaluation of american option prices in a path integral framework using fourier-hermite series expansions. J. Econom. Dynam. Control, 23(9-10), 1999. [5] S. Figlewski and B . Gao. The adaptive mesh model: a new approach to efficient option pricing. J. Fin. Econ, 53:313-351, 1999. [6] B . Gao J-z. Huang and M . Subrahmanyam. The valuation of american options using the decomposition thechnique. J. Econom. Dynam. Control, 24:1783-1827, 2000. [7] R. Geske and H . E . Johnson. The american put options valued analytically. J. Fin, 39:1511-1524, 1984. 50 [8] J . Hul l . Options, Futures, and Other Derivatives. A SimonSchuster Company, 1997. [9] J .D. Evans R . A . Kuske and J .B . Keller. Behaviuor near expiry of american options with dividends. J. Math. Finance, 12:219-237, 2002. [10] S. Ikonen and J . Toivanen. Operator splitting methods for american option pricing. Applied Mathematics Letters, 17:809-814, 2004. [11] H . Johnson. A n analytic approximation for the american put price. J. Fin. Quan. Anal., 18:141-148, 1983. [12] b. Kamrad and P. Ritchken. Mult inomial approximating models for options with k-state variables. Management Science, 37:1640 -1652, 1991. [13] M . D . Koulisianis and T . S . A . Papatheodorou. Moving index method for the solution of the american options valuation problem. Math. Comput. Simulation, 54(4-5), 1999. [14] R . A . Kuske and J .B . Keller. Optimal exercise boundary for an american put option. Appl. Math. Finance, 5:107-116, 1998. [15] L . W . MacMil lan . A n analytic approximation for the american put price. Ad-vanced in Futures and Options Research, 1:119-139, 1986. [16] R . C . Merton. Theory of rational option pricing. Bell J. of Econ. and Manage-ment Science, 4:141-183, 1973. [17] G . Meyer and J.van der Hoek. The valuation of american options with the method of lines. Advanced in Futures and Options Research, 9:265-286, 1997. 51 \ [18] B . F . Nielsen O. Skavhaug and A . Tveito. Penalty and front - fixing methods for the numerical solution of american option problems. J. Comput. Fin., 5(4):69-97, 2002 [19] P.Ritchken. On pricing barrier options. Journal of Derivatives, 3:19-28, 1995. [20] E . Reiner and M . Rubinstein. Breaking down the barriers. Risk, 4:28-35, 1991. [21] D . Sevcovic. Analysis of the free boundary for the pricing of an american call option. J. of Applied Mathematics, 12:25-37, 2001. [22] P. Wil lmot S. Howison and J . Dewynne. The Mathematics of Financial Deriv-atives. Cambridge University Press, 1996. [23] L . W u and Y . K Kwok. A front - fixing finite difference method for the valuation of american options. J. Fin. Eng., 6(2):83-97, 1997 [24] K . R . Vetzal R. Zvan and P . A . Forsyth. P D E methods for pricing barrier op-tions, computational aspects of complex securities. J. Econom. Dynam. Con-trol, 24( 11-12): 1563-1590, 2000. [25] http:/ /www.derivativesmodels.com/FDX/Barrier.htm 52 Appendix A A . l Optimal exercise boundary near expiry We know that there is no explicit global expression of Sf(t), but for times near expiry i.e. TF — t <C 1, Keller, Kuske and Evans ([9]) have developed the following expressions: Sf{t) ~K- Ka^/2{TF - t) a [a2(TF - t)/2] 0 < D < r, ( A . l ) (A.2) Sf(t) ~K- Kay]2(TF - t) \n[l/{A^D{TF - t))\ D = r, Sf(t) ~ L K ( \ - aa0^2{TF - t)^j D > r . (A .3) Here the positive function Q(T) is the solution of the equation: „ - a 2 ( r ) r , x / rl ~ l im ' 7 T x—>Xf(r) (r1/2{p - u) + 2VTOL(T)^ ( / erfc B(x, z, r)dz - r l | 1 B T [ V l T ] e - B 2 l * ) ^ ] r f z j - [p + u - (p - u)2} j ^ z e ^ ^ d z ^ J B[x, Z,T]BT[X, z,T}e-B^x'z'Tkz J 0 / (A.4) where B[X,Z,T] is defined by B(x,z,-[xQ - x] / 2 r V 2 - zV2a{Tz) (A .5) 53 Expressions (A.4) and (A.5) are derived from an application of Green's function to the solution of the dimensionless formulation of equations (2.1) - (2.4) as a forward - parabolic system (see [9]). In equation (A.3), ao is a constant that solves the transcendental equation which follows from the two boundary conditions in (2.2). A n application of the bisection method yields ao « 0.4517. In equation (A.3) we have that for D > r the boundary is parabolic, and as t tends to Tp, Sj(t) tends to the value rK/D, which is less than the strike price K. In equation ( A . l ) we have that for D < r, as t tends to Tp, Sf(t) tends to K, but it is not parabolic. (A.6) A . 2 T h e Dimens ion less P r o b l e m C h a n g e o f V a r i a b l e s S = K • e ,x (A.7) t = T F - ^ r (A.8) a P = Ke-Tp{x:t) + K-S (A.9) Sf = K- e (A.10) 54 D i m e n s i o n l e s s fo rward p a r a b o l i c s y s t e m dp x/(0) = x 0 , ^ -log(-) if v > p P Xf(r) < x p = — = 0 , , at ir = x / ( r ) p ~ epT(ex — 1) , as cc —> + 0 0 p = mace (e x — 1 ,0) , at r = 0 0 if v < p \ 55 Appendix B Matlab Codes B.0.1 Optimal Exercise Boundary for American Put Option with Continuous Dividend Y i e l d for D > r '/, so lve the American Put Opt ion w i t h continuous D i v i d e n d y i e l d (D > r ) 7, u s ing f i x e d po in t method °/0 i npu t : parameters of model 7, output : e x e r c i s e boundary f u n c t i o n 0EB_f ix = i m p . d i v i d e n d s C s i g m a . r . D . n . m . I T ) ; %the s t r i k e p r i c e and the e x p i r y t ime K = 5; Tf = 1; ' / . spa t ia l d i s c r e t i z a t i o n xmax = 15; dx = xmax/m; x = d x * [ 0 : m - l ] ; 56 tO = 0; '/ . i n i t i a l time '/.time step away from expiry dt_big=(Tf-tO)/n; '/.define the points where the time step w i l l change t l = 0.95; t2 = 0.75; t = [ t l t 2 ] ; nn = n/2; '/.define the time step near expiry dt_small = (Tf - t l ) / n n ; '/.define the l i n e a r function f o r the nonuniform time step dt = 2*(dt_small - d t _ b i g ) * t + (2*dt_big - dt_small); '/.the time mesh from 0.95 to l(near expiry) t t _ 0 = [ T f - d t _ s m a l l : - d t _ s m a l l : t ( l ) ] ; °/0the time mesh from 0.75 to 0.95 t t _ l = [ t ( l ) - d t ( D : - d t ( l ) : t(2)] ; '/.the time mesh from 0.5 to 0.75 tt_ 2 = [t(2)-dt(2) : - dt(2):0.5]; '/, the time mesh from 0 to 0.5(away from expiry) t_away = [0.5-dt_big:-dt_big:t0]; 57 '/.the nonuniform time mesh over the whole option's l i f e T = [tt_0 t t _ l tt_2 t_away]; tplot = [Tf T]; '/.define the time step for any time from Tf to 0 dt2 = [dt_small*ones(l,length(tt_0)) dt(1)*ones(1,length(tt_l)) dt(2)*ones(l,length(tt_2)) dt_big*ones(1,length(t_away))]; '/.the total number of points in the time grid NN = length(tt_0)+ length(tt_l)+length(tt_2)+length(t_away); 7,define Ksmall and Kbig used in the change of variables Kshift = 1; Ksmall = K; Kbig = K+Kshift; '/.initial guess for boundary alfa = 0.4517; b = K*(r/D)*(l- sigma*alfa*sqrt(2*(Tf-T))); b = [r*K/D,b]; '/.plot the i n i t i a l guess plot(b,tplot) pause(0.2) '/, begin fixed point iterations for iter = 1:IT '/.exercise boundary and t-derivative 58 beta = b(2:NN+l); '/.the derivatives for j = 1: NN - 1 betap(j) = (b(j) - b(j+2))/(dt2(j + 1)+ d t 2 ( j ) ) ; end betap(NN) = (b(NN)-b(NN + l))/dt_Mg; "/.the change of variable for pl o t t i n g s = x2s(x(l:m),b(l).Ksmall,Kbig); '/.the f i n a l condition Pplot = max(K-s,0); P = max(K-s,0)'; '/.time stepping backwards to solve Black-Scholes equation for i t = 1:NN '/.this i s coordinate change from (s,t) to x s = x2s(x(l:m),beta(it),Ksmall,Kbig); [phi,phip,phipp] = cutoff(s.Ksmall,Kbig); '/.the coefficients for PDE A = 0.5*sigma"2*s.~2.*(l-beta(it)*phip) . ' 2 ; B = -0.5*sigma"2*s.~2.*beta(it).*phipp ... -betap(it)*phi+(r-D)*s.*(l-beta(it)*phip); 59 '/.differential operators Ml = -2*diag(A)+diag(A(l:m-l),+l)+diag(A(2:m),-1); Ml(l ,2) = 2*A(1) M2 = -diag(B(2:m),-l)+diag(B(l:m-l),1); M2(l,2) = 0; M = -dt2(it)*(Ml/(dx-2)+M2/(2*dx)-r*eye(m))+eye(m); "/.solve the system for each time C = zeros(m,l); C(l) = -B(l)+A(l)*2/dx; C = dt2(it)*C; P = M\(P+C); Pplot = [Pplot;P'] ; end p lot (b , tplot ,K-Pplot ( : ,1) , tp lot ) ; xlabel('S($)'); y l a b e K ' t (Years)'); t i t le( 'progress towards solution - exercise boundary'); pause (0.2); '/exercise boundary for next i teration b=K-Pplot(: ,1) ; end B . 0 . 2 T h e C h a n g e o f V a r i a b l e f r o m x t o S function s = x2s(x,beta,Ks,Kb) '/.given an array x, find array s satisfying x=s-beta*phi(s) °/,Kbig and Ksmall are parameters used to define phi(s) 60 •/.define dK dK = 0 .1*K; dKs = m a x ( . 0 1 , ( K s - b e t a ) . / 2 ) ; '/.bottom end of curvy b i t x = - b e t a + K s m a l l , top end x = Kbig + dK '/.find a l l indeces f o r which S < K s m a l l , i . e . x < K s m a l l - b e t a ; I l e f t = f i n d ( x < -beta+Ks); '/.find a l l indeces f o r which S > Kbig + dK, i . e . x < Kbig + dK; I r i g h t = f i n d ( x > Kb+dK); '/.find a l l indeces f o r which Ksmall < x < Kbig + dK; i f ( i s e m p t y ( I l e f t ) = = 1 ) Imid = 1 : m i n ( I r i g h t - l ) ; e l s e Imid = m a x ( I l e f t ) + l : m i n ( I r i g h t - l ) ; end '/.leftend of curvy b i t K s m a l l , r i g h t e n d of curvy b i t Kbig+dK '/.sample curvy b i t at 100 p o i n t s nc = 1 0 0 ; '/.define the step f o r sample p o i n t s s i d s i = (Kb + dK-(Ks - d K s ) ) . / n c ; °/,define sample p o i n t s s i 61 s i = (Ks - d K s ) : d s i : ( K b + dK) ; "/.define f u n c t i o n p h i and i t s f i r s t and second d e r i v a t i v e s '/.for the sample p o i n t s [ p h i s a m p . p h i p . p h i p p ] = c u t o f f ( s i , K s , K b ) ; 7,now use phisamp to h o l d approximat ion to x phisamp=si-beta*phisamp; ' / .find s p o i n t s f o r K s m a l l < x < K b i g + dK u s i n g Newton's method; f o r j = Imid I = f ind(x( j ) -phisamp<=0); i t = m i n ( I ) ; p h i s a m p ( i t ) ; i f a b s ( p h i s a m p ( i t ) - x ( j ) ) <= l . e - 6 ; s ( j ) = s i ( i t ) ; e l s e s ( j ) = s i ( i t ) - d s i * ( p h i s a m p ( i t ) - x ( j ) ) / ( p h i s a m p ( i t ) - p h i s a m p ( i t - l ) ) ; end end '/ .find s p o i n t s f o r x < Ksmal l and x > K b i g + dK i f ( i s e m p t y ( l l e f t ) = = 0 ) s ( I l e f t ) = b e t a + x ( I l e f t ) ; end 62 s ( I r i g h t ) = x ( I r i g h t ) ; B.0.3 Function </? and its Drivatives f u n c t i o n [phi,phip,phipp] = cutoff(s,Ks,Kb) °/,cutoff f u n c t i o n where '/.phi = 1 s < Ksmall '/.phi = 0 s > Kbig+dK '/.phi smooth and decreasing i n between '/.define the constants involved i n the d e f i n i t i o n of cc = 4.25; •/.define dK dK = 0.1*K; '/.find a l l the i n d i c e s f o r which s < Ksmall I l e f t = find(s<Ks); m a x ( I l e f t ) y °/,find a l l the i n d i c e s f o r which s > Kbig + dK I r i g h t = find(s>Kb+dK); m i n ( I r i g h t ) ; •/.define the midpoint of the i n t e r v a l [Ksmall, Kbig] K2 = (Ks+Kb)./2; 63 "/.define function phi on [Ksmall, Kbig + dK] phi = 0.5-0.5*erf(cc*(s-K2)/dK); "/.define the f i r s t and the second derivative of phi on (Ksmall, Kbig + dK) phip = -cc/(sqrt(pi)*dK)*exp(-(cc*(s-K2)/dK).~2); phipp = 2.*cc.~3/(sqrt(pi)*dK~3)*(s-K2).*exp(-(cc*(s-K2)/dK)."2); "/.define function phi and i t s f i r s t and the second derivatives for s < Ksmall phi(Ileft) = ones(size(Ileft)); phip(Ileft)=zeros(size(Ileft)); phipp(Ileft)=zeros(size(Ileft)); "/.define function phi and i t s f i r s t and the second derivatives for s > Kbig + dK phi(Iright) = zeros(size(Iright)); phip(Iright)=zeros(size(Iright)); phipp(Iright)=zeros(size(Iright)); B . 0 . 4 C o n v e r g e n c e f o r t i m e s t e p r e f i n e m e n t °/,the convergence of the method when dx - constant and the time step is refined "/.plot a l l the errors on the same graph (D = r , D < r , and D > r) "/.parameters of the model 64 sigma = 0.25; Tf = l ; t 0 = 0; r = 0.1; '/.the "exact" OEB load exact.mat '/.dividends y i e l d f o r each case D_greater = 0.15; D_equal =0.1; D_div = 0.05; '/.the i n i t i a l number of time points n= 20; '/.the i n i t i a l number of s p a t i a l points m = 60; '/.define the time refinement (dt =1/(20*H)) H = 2.'(0:1:5); Nt = length(H); '/. i n i t i a l i z e the vector error f o r each case ERRDR_equal = zeros(NN+1,Nt); ERR0R_D = zeros(NN+1,Nt); ERR0R_div = zeros(NN+l,Nt); HH = z e r o s ( l , N t ) ; DT =1./(20*H); 65 7,the number of iterations IT = 25; '/.define the time step for time away from expiry dt_big=(Tf-tO)/n; 7, define the points where the time step w i l l change t l = 0.95; t2 = 0.75; t = [tl t2]; nn =n/2; '/.define the time step for time near expiry dt_small = (Tf - tl)/nn; '/.define the linear function for the nonuniform time step dt = 2*(dt_small - dt_big)*t + (2*dt_big - dt.small); '/.the time mesh from 0.95 to 1 (near expiry) tt_0 = [Tf-dt_small:-dt_small:t(l)] ; '/.the time mesh from 0.75 to 0.95 t t _ l = [t(l)-dt ( D : - dt ( l ) : t(2)] ; '/.the time mesh from 0.5 to 0.75 tt_2 = [t(2)-dt(2) : - dt(2):0.5]; 66 7,the time mesh from 0 to 0.5(away from expiry) t_away = [0 .5 -dt_big: -dt_big : t0 ] ; '/.the nonuniform time mesh over the whole option's l i f e T = [tt_0 t t _ l tt_2 t_away]; tplot = [Tf T]; '/.the number of the data points over the whole option's l i f e NN = length(tt_0)+ length(tt_l)+length(tt_2)+length(t_away); '/ . initialize OEB vectors for each step refinement for each case Sf_equal = zeros (M+1, Nt) ; Sf_D = zeros(NN+l,Nt); Sf_div = zeros(NN+1,Nt); '/.the case when D = r for i = l:Nt '/define the time step refinement dt=l/(n*H(i)); '/.perform a sequence of runs for each time step refinement defined previously '/.and dx kept constant 0EB_fix_equal = imp_equal(sigma,r,D_equal,n*H(i),m,IT); 67 70the optimal exercise boundary vector for each spatial step refinement Sf_equal (: , i ) = 0EB_fix_equal(l:H(i):NN*H(i)+l); "/.find the vector error ERROR_equal(:,i) = abs(Sf_equal(:,i) - Sf_exact_equal); "/.find the average of the error error_dt_equal(i) = (sum(ERROR_equal(:,i)))/NN; HH(i) = H(i); end "/.the case when D > r for i = l:Nt "/.define the time step refinement dt=l/(n*H(i)); "/.perform a sequence of runs for each time step refinement defined previously "/.and dx kept constant 0EB_fix_D = imp_dividends(sigma,r,D_greater,n*H(i),m,IT); "/.the optimal exercise boundary vector for each spatial step refinement Sf_D (:, i) = OEB_fix_D(l:H(i):NN*H(i)+l); "/.find the vector error 68 ERROR_D(:,i) = abs(Sf_D(:,i) - Sf_exact_D); '/.find the average of the error error_dt_D(i) = (sum(ERROR_D(:,i)))/NN; HH(i) = H(i ) ; end '/.the case when D < r for i = l:Nt "/define the time step refinement dt=l/(n*H(i)); '/perform a sequence of runs for each time step refinement defined previously '/and dx kept constant OEB_fix_div = imp_d(sigma,r,D_div,n*H(i),m,IT); '/the optimal exercise boundary vector for each spatial step refinement Sf_div (: , i ) = 0EB_fix_div(l:H(i):NN*H(i)+l); '/.find the vector error ERR0R_div(:,i) = abs(Sf_div(:,i) - Sf_exact_div); '/.find the average of the error error_dt_div(i) = (sum(ERROR_div(:,i)))/NN; 69 HH(i) = H(i); end '/.plot the error vs. the number of time steps plot(log2(HH),error_dt_div,'d-', log2(HH),error_dt_equal,'kp-', log2(HH),error_dt_D,'o-' , ' l inewidth',2); '/.define the legend legl=['D = ',num2str(D_div)]; leg2=['D = ',num2str(D_equal)]; leg3=['D = ',num2str(D_greater)]; legend(legl,leg2,leg3) '/.define the axis' label xlabel('log2(N)') y labe lC ERROR') B . 0 . 5 C o n v e r g e n c e f o r t h e t i m e a n d s p a t i a l s t e p r e f i n e m e n t s '/.convergence of the method when both the time step and the spatial step are refined °/,plot a l l the errors on the same graph (D = r , D < r , and D > r) '/.parameters of the model sigma = 0.25; Tf = 1; tO = 0; r = 0.1; '/.the "exact" OEB 70 load exact.mat '/.dividends y ie ld for each case D_greater =0.15; D_equal = 0.1; D_div = 0.05; °/,the i n i t i a l number of time points n= 20; '/.the i n i t i a l number of spatial points m = 60; '/.define the refinement for time step (dt =1/(20*H)) p = 0:1:5; H = 2.~p; Nt = length(H); '/ . initialize the error matrices for each case ERR0R_equal = zeros(NN+1,Nt); ERR0R_D = zeros(NN+1,Nt); ERR0R_div = zeros(NN+l,Nt); '/ . initialize the vector error for each case error_dt_equal = zeros(l .Nt); error_dt_D = zeros(l ,Nt); error_dt_div = zeros(l .Nt); 71 HH = zeros(l ,Nt); DT = l./(20*H); °/,the number of iterations IT = 25; "/.define the time step for time away from expiry dt_big=(Tf-tO)/n; '/, define the points where the time step w i l l change t l = 0.95; t2 = 0.75; t = [ t l t2]; nn =n/2; '/.define the time step for time near expiry dt.small = (Tf - t l ) /nn; '/.define the l inear function for the nonuniform time step dt = 2*(dt_small - dt_big)*t + (2*dt_big - dt.small); '/.the time mesh from 0.95 to l(near expiry) tt_0 = [Tf-dt_small:-dt_small:t(l)] ; '/.the time mesh from 0.75 to 0.95 t t _ l = [ t ( l ) -dt (D : - d t ( l ) : t(2)] ; 72 '/.the time mesh from 0.5 to 0.75 tt_2 = [t(2)-dt(2) : - dt(2):0.5]; '/, the time mesh from 0 to 0.5(away from expiry) t_away = [0.5-dt_big:-dt_big:t0]; '/.the nonuniform time mesh over the whole option's l i f e T = [tt_0 t t _ l tt_2 t_away]; tplot = [Tf T]; '/.the number of the data points over the whole option's l i f e NN = length(tt_0)+ length(tt_l)+length(tt_2)+length(t_away); '/.define the OEB matrix for each case Sf_equal = zeros(NN+1,Nt); Sf_D = zeros(NN+1,Nt); Sf_div = zeros(NN+1,Nt); '/.the case when D = r for i = l:Nt '/.define the time step refinement dt=l/(n*H(i)); '/.perform a sequence of runs for each time step and spatial step refinements '/.defined previously 0EB_fix_equal = imp_equal(sigma,r,D_equal,n*H(i),m*(p(i)+l),IT); 73 '/.define the optimal exercise boundary vector for each time step and '/.spatial step refinement Sf_equal (: , i ) = 0EB_fix_equal(1:H(i):NN*H(i)+l); '/.find the vector error ERR0R_equal(:,i) = abs(Sf_equal(:,i) - Sf_exact_equal); '/.find the average of the error error_dt_equal(i) = (sum(ERROR_equal(:,i)))/NN; HH(i) = H(i); end '/.the case when D > r for i = l:Nt "/.define the time step refinement dt=l/(n*H(i)) ; "/.perform a sequence of runs for each time step and spatial step refinements '/.defined previously 0EB_fix_D = imp_dividends(sigma,r,D_greater,n*H(i),m*(p(i)+l),IT); '/.define the optimal exercise boundary vector for each time step and '/.spatial step refinement Sf_D (:, i) = OEB_fix_D(l:H(i):NN*H(i)+l); 74 "/.find the vector error ERR0R_D(:,i) = abs(Sf_D(:,i) - Sf_exact_D); '/.find the average of the error error_dt_D(i) = (sum(ERROR_D(:,i)))/NN; HH(i) = H(i); end "/.the case when D < r for i = l:Nt "/.define the time step refinement dt=l/(n*H(i)) ; "/.perform a sequence of runs for each time step and spatial step refinements "/.defined previously OEB_fix_div = imp_d(sigma,r,D_div,n*H(i),m*(p(i)+l),IT); "/.define the optimal exercise boundary vector for each time step and "/.spatial step refinement Sf_div (: , i ) = 0EB_fix_div(l:H(i):NN*H(i)+l); "/.find the vector error ERR0R_div(: , i ) = abs(Sf_div(:,i) - Sf_exact_div); "/.find the average of the error error_dt_div(i) = (sum(ERROR_div(:,i)))/NN; HH(i) = H(i); 75 end "/.plot the error vs. the number of time steps plot(log2(HH),error_dt_div,'d-', log2(HH),error_dt_equal,'kp-', log2(HH),error_dt_D , 'o-' , ' l inewidth',2); °/,define the legend legl=['D = ',num2str(D_div)] ; leg2=['D = ',num2str(D_equal)]; leg3=['D = ',num2str(D_greater)]; legend(legl,leg2,leg3) "/.define the axis' label xlabel('log2(N)') ylabel('ERROR') B.0.6 Optimal Exercise Boundary for American Put Barrier tion with Continuous Dividend Yie ld for D > r function 0EB_fix = barrier_dividends(sigma,r,D,n,m,IT,Xb); '/.the strike price and the expiry time K = 5; Tf = 1; "/.spatial discretization xmax = Xb; dx = xmax/m; x = dx*[0:m-l]'; 76 tO = 0; "/.initial time "/time step away from e x p i r y dt_big = (Tf-t0)/n; "/define the p o i n t s where the time step w i l l change t l = 0.95; t2 = 0.75; t = [ t l t2]; nn = n /2; "/define the time step near e x p i r y dt_small = (Tf - t l ) / n n ; "/define the l i n e a r f u n c t i o n f o r the nonuniform time step dt = 2*(dt_small - d t _ b i g ) * t + (2*dt_big - d t _ s m a l l ) ; "/the time g r i d from 0.95 to K n e a r e x p i r y ) tt_0 = [ T f - d t _ s m a l l : - d t _ s m a l l : t ( l ) ] ; '/the time g r i d from 0.75 to 0.95 t t _ l = [ t ( l ) - d t U ) : - d t ( l ) : t (2)]; "/the time g r i d from 0.5 to 0.75 tt_2 = [t(2)-dt(2) : - dt(2):0.5]; '/the time g r i d from 0 to 0.5(away from expiry) 77 t_away = [0.5-dt_big:-dt_big:t0]; "/the nonuniform time grid over the whole option's l i f e T = [tt_0 t t _ l tt_2 t_away]; tplot = [Tf,T]; "/.define the time step for any time from Tf to 0 dt2 = [dt_small*ones(l,length(tt_0)) dt(l)*ones(l , length(tt_l)) dt(2)*ones(l,length(tt_2)) dt_big*ones(l,length(t_away))]; "/.the total number of points in the time grid NN = length(tt_0)+ length(tt_l)+length(tt_2)+length(t_away); "/.define Ksmall = K- and Kbig = K+ esed in the change of variable Kshift = 1; Ksmall = K; ndx =3; np = ndx;*dx; Kbig = min(K+Kshift,Xb-np); "/.initial guess for boundary alfa = 0.4517; b = K*(r/D)*( l - sigma*alfa*sqrt(2*(Tf-T))); b = [r*K/D,b]; "/.plot the i n i t i a l guess plot(b,tplot) pause(0.2) "/begin fixed point iterations for i ter = 1:IT 78 '/.exercise boundary and t-derivative beta = b(2:NN+l); '/.the derivatives for j = 1: NN - 1 betap(j) = (b(j) - b(j+2))/(dt2(j + 1)+ dt2(j)); end betap(NN) = (b(NN)-b(NN + l)) /dt_big; '/.the change of variables s = x2sl(x(l:m),b(l).Ksmall,Kbig,np); '/.the f ina l condition Pplot = max(K-s.O).*(s < Xb); P = Pplot> ; '/.time stepping backwards to solve Black-Scholes equation for it=l:NN '/.this is coordinate change from (s,t) to x Ksmall = K; s = x2sl(x(l:m),beta(it).Ksmall,Kbig,np); [phi.phip.phipp] = cutoffl(s,Ksmall,Kbig,np); °/,the coefficients for PDE A = 0.5*sigma~2*s•"2.*(l-beta(it)*phip)."2; B = -0.5*sigma~2*s."2.*beta(it).*phipp . . . 79 -betap(it)*phi+(r-D)*s.*(l-beta(it)*phip); ' /differential operators Ml = -2*diag(A)+diag(A(l:m-l),+l)+diag(A(2:m),-l); Ml(l ,2) = 2*A(1) M2 = -diag(B(2:m),-l)+diag(B(l:m-l),1); M2(l,2) = 0; M = -dt2(it)*(Ml/(dx~2)+M2/(2*dx)-r*eye(m))+eye(m); '/solve the system for each time C = zeros(m,l); C(l)= -B(l)+A(l)*2/dx; C = dt2(it)*C; P = M\(P+C); P = P.*(s < Xb)'; Pplot = [Pplot;P']; end plot (b , tplot ,K-Pplot ( : ,1) , tp lot ) ; xlabel('S($)'); ylabel('t(Years)') t i t le( 'progress towards solution - exercise boundary'); pause (0.2); '/exercise boundary for next i teration b = K-Pplot(: ,1) ; end 80
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- An iterative method for pricing the American put options...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
An iterative method for pricing the American put options with dividends Coman, Florica 2005
pdf
Page Metadata
Item Metadata
Title | An iterative method for pricing the American put options with dividends |
Creator |
Coman, Florica |
Date Issued | 2005 |
Description | In this thesis we present a fixed-front approach for pricing the American Put Options with dividends, combined with an iterative method. This approach transforms the free boundary problem into a nonlinear parabolic differential equation defined on a fixed spatial domain. The resulting equation has been solved using an implicit numerical scheme. The method can be easily applied to Barrier Options with dividends and numerical results were shown to compare well with the optimal exercise boundary for different dividend yields and different values of the barrier. All the programming work has been done in MATLAB and sample of codes can be found in Appendix B. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0079512 |
URI | http://hdl.handle.net/2429/16506 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2005-0413.pdf [ 2.28MB ]
- Metadata
- JSON: 831-1.0079512.json
- JSON-LD: 831-1.0079512-ld.json
- RDF/XML (Pretty): 831-1.0079512-rdf.xml
- RDF/JSON: 831-1.0079512-rdf.json
- Turtle: 831-1.0079512-turtle.txt
- N-Triples: 831-1.0079512-rdf-ntriples.txt
- Original Record: 831-1.0079512-source.json
- Full Text
- 831-1.0079512-fulltext.txt
- Citation
- 831-1.0079512.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079512/manifest