summary:This paper is devoted to barrier options and the main objective is to develop a sufficiently robust, accurate and efficient method for computation of their values driven according to the well-known Black-Scholes equation. The main idea is based on the discontinuous Galerkin method together with a spatial adaptive approach. This combination seems to be a promising technique for the solving of such problems with discontinuous solutions as well as for consequent optimization of the number of degrees of freedom and computational cost. The appended numerical experiment illustrates the potency of the proposed numerical scheme

Hozman, Jiří

English

Institute of Mathematics AS CR, v. v. i.

PANM 16Jiří HozmanValuing barrier options using the adaptive discontinuous Galerkin methodIn: Jan Chleboun and Karel Segeth and Jakub Šístek and Tomáš Vejchodský (eds.): Programs and Algorithms ofNumerical Mathematics, Proceedings of Seminar. Dolní Maxov, June 3-8, 2012. Institute of Mathematics AS CR,Prague, 2013. pp. 94--99.Persistent URL: http://dml.cz/dmlcz/702712Terms of use:© Institute of Mathematics AS CR, 2013Institute of Mathematics of the Academy of Sciences of the Czech Republic provides access to digitized documentsstrictly for personal use. Each copy of any part of this document must contain these Terms of use.This paper has been digitized, optimized for electronic delivery and stamped with digitalsignature within the project DML-CZ: The Czech Digital Mathematics Libraryhttp://project.dml.czPrograms and Algorithms of Numerical Matematics 16J. Chleboun, K. Segeth, J. Š́ıstek, T. Vejchodský (Eds.)Institute of Mathematics AS CR, Prague 2013VALUING BARRIER OPTIONS USING THE ADAPTIVEDISCONTINUOUS GALERKIN METHODJǐŕı HozmanTechnical University of Liberec, Faculty of Science, Humanities and EducationStudentská 2, 461 17, Liberec, Czech Republicjiri.hozman@tul.czAbstractThis paper is devoted to barrier options and the main objective is to developa sufficiently robust, accurate and efficient method for computation of their valuesdriven according to the well-known Black-Scholes equation. The main idea is basedon the discontinuous Galerkin method together with a spatial adaptive approach.This combination seems to be a promising technique for the solving of such problemswith discontinuous solutions as well as for consequent optimization of the numberof degrees of freedom and computational cost. The appended numerical experimentillustrates the potency of the proposed numerical scheme.1. IntroductionDuring the last decade, financial models have acquired increasing popularity inoption pricing. The valuation of different types of option contracts is very importantin modern financial theory and practice – especially exotic options such as discretebarrier options. Most of the analytical formulas for these options is limited by strongassumptions, which led to the application of numerical methods instead.Therefore, the main goal of this paper is to develop an efficient, robust andaccurate numerical method for the barrier option pricing problem, which arisesfrom the concept of the discontinuous Galerkin (DG) approach for the space semi-discretization, for more details see [5], and the backward Euler scheme for the dis-cretization of the resulting ODE systems. In order to increase the efficiency of theproposed method additionally, this approach is equipped with an h-adaptivity tech-nique based on regularity and residual indicators, cf. [1, 2]. The resulting numericalscheme is applied to a standard problem of discrete double barrier option pricing.2. Barrier option pricing modelIn what follows, we consider the double time-independent discrete barrier knock-out option, i.e. option that expires worthless if one of the two barriers has been hitat a monitoring date, see e.g. [1] and [6]. We denote by x the price of an underlying94asset (e.g. stock) and by t the time to expiry of the option and let M := {0 = tM0 <tM1 < . . . < tMl−1 < tMl = T} be the set of monitoring dates and B− be the lowerbarrier and B+ the upper barrier active only at discrete instances tMl ∈M .Let Ω≡(0, Smax), 0 < B− < B+ < Smax, be a bounded open interval and T standsfor maturity. The price u : QT = Ω × (0, T ) → IR of the discrete barrier optionsatisfies the Black-Scholes partial differential equation with initial and boundaryconditions:∂∂tu(x, t)−12σ2x2∂2∂x2u(x, t)− rx∂∂xu(x, t) + ru(x, t) = 0 in QT , (1)u(0, t) = 0 and u(Smax, t) = 0, (2)u(x, 0) ={max(x−K, 0) · χ[B−,B+], (call)max(K − x, 0) · χ[B−,B+], (put), x ∈ Ω, (3)where σ > 0 and r > 0 are constant model parameters denoting the volatility ofstock price and the risk-free interest rate, respectively.From the mathematical point of view the problem (1)–(3) represents a convection-diffusion-reaction equation equipped with a set of two homogeneous Dirichlet bounda-ry conditions (2) prescribed at the endpoints of interval (0, Smax) and with the initialcondition (3), where symbol K stands for the strike price and χ[B−,B+] denotes thecharacteristic function of the barrier interval.Moreover the discrete monitoring of the contract introduces an updating of thesolution u(x, t) at the monitoring dates tMl ∈M :u(x, tMl ) = limε→0+u(x, tMl − ε) · χ[B−,B+]. (4)3. Discontinuous Galerkin discretizationLet Th (h > 0) be a family of partitions of the closure Ω = [0, Smax] of the do-main Ω into N closed mutually disjoint subintervals Ik = [xk−1, xk] with length hk :=xk − xk−1. Then we set Th = {Ik, 1 ≤ k ≤ N} with spatial step h := max1≤k≤N hkand call interval Ik an element. We additionally assume that the following conditionsare satisfied:∃Cq ≥ 1 : hk ≤ Cqhk′ ∀ Ik, Ik′ ∈ Th sharing a node (local quasi-uniformity) (5)∃ k1, k2 ∈ IN such that xk1 = B− and xk2 = B+ (barrier consistency) (6)The DG method can handle different polynomial degrees over elements. There-fore, we assign a positive integer pk as a local polynomial degree to each Ik ∈ Th.Then we set the vector p = {pk, Ik ∈ Th}. Over the triangulation Th we define thefinite dimensional space of discontinuous piecewise polynomial functions:Shp ≡ Shp(Ω, Th) = {v; v|Ik ∈ Ppk(Ik) ∀ Ik ∈ Th}, (7)where Ppk(Ik) denotes the space of all polynomials of degree ≤ pk on Ik, Ik ∈ Th.Consequently, the approximate solution of the continuous problem (1)–(4) is soughtin the space Shp.95Let us denote v(x±k ) = limε→0+ v(xk ± ε). Then we define the jump and averageof v at inner points xk of Ω by [v(xk)]=v(x−k )−v(x+k ) and 〈v(xk)〉=12(v(x−k ) + v(x+k )),respectively. We also extend the definition of jump and mean value for endpointsof Ω, i.e. [v(x0)] = −v(x+0 ), 〈v(x0)〉 = v(x+0 ), [v(xN)] = v(x−N ) and 〈v(xN)〉 = v(x−N).Firstly, we recall the space semi-discrete DG scheme presented in [4] and [5]. Tothis end we introduce the following bilinear forms:aΘh (u, v) =N−1∑k=0∫ xk+1xk12σ2x2∂u(x, t)∂xv′(x) dx−N∑k=0〈12σ2x2k∂u(xk, t)∂x〉[v(xk)]+ΘN∑k=0〈12σ2x2k v′(xk)〉[u(xk, t)], (8)bh(u, v) =−N−1∑k=0∫ xk+1xk(σ2 − r)xu(x, t) v′(x) dx+N∑k=0H(u(x−k , t), u(x+k , t))[v(xk)], (9)Jωh (u, v) =N∑k=0ωk[u(xk, t)] [v(xk)]. (10)The crucial item of the DG formulation is the treatment of the linear convectionand diffusion terms. For the convection form bh we treat its terms with the aid ofa numerical flux H , see [3]. The diffusion form aΘh includes stabilization terms whichare added to the formulation of the problem in order to guarantee the stability of thenumerical scheme. According to the value of parameter Θ, we speak of symmetric(Θ = −1), incomplete (Θ = 0) or nonsymmetric (Θ = 1) variants. Furthermore, inorder to replace the inter-element discontinuities, the semi-discrete scheme is com-pleted with the penalty Jωh weighted by the penalty parameter function ωk definedin the spirit of [4]. Let us note that the right-hand side term vanishes due to theprescribed homogeneous Dirichlet boundary conditions in (2).In order to simplify the notation we define the bilinear form:BΘh (u, v) := aΘh (u, v) + bh(u, v) + αJωh (u, v) + (2r − σ2)(u, v), α > 0, (11)where (·, ·) denotes inner product and the forms aΘh (·, ·), bh(·, ·) and Jωh (·, ·) are givenby (8), (9) and (10), respectively. The value of multiplicative constant α before thepenalty form Jωh depends on the properties of diffusion term, see [4]. Finally, we endup with the following DG formulation for the semi-discrete solution uh(t) ∈ Shp:ddt(uh(t), vh) + BΘh (uh(t), vh) = 0 ∀ vh ∈ Shp, ∀ t ∈ (0, T ), (12)which represents an ODE system and due to bilinearity of form (11) we can easilydiscretize (12) by the implicit Euler method. Let 0 = t0 < t1 < . . . < tr = T bea partition of [0, T ] with time steps τl ≡ tl+1 − tl, l = 0, 1, . . . , r − 1. We definethe approximate solution of problem (1)–(4) as functions ulh ≈ uh(tl), tl ∈ [0, T ],l = 0, . . . , r − 1, satisfying the following numerical scheme:(ul+1h , vh)+ τlBΘh(ul+1h , vh)=(ulh, vh)∀ vh ∈ Shp, (13)ul+1h := ul+1h · χ[B−,B+] ∀ tl+1 ∈M, (14)96where u0h is Shp-approximation of u0. The discrete problem (13) is equivalent toa system of linear algebraic equations at each time level tl+1 ∈ [0, T ].4. Mesh adaptationIn this section, we introduce an h-adaptive DG technique for the solution ofproblem (1)–(4). Since we deal with nonstationary problems, it is suitable to useadaptive mesh refinement during the computation in order to improve the numericalsolution and to optimize the number of degrees of freedom and computational cost,consequently.We start from a uniform coarse grid T0,h := Th and construct at each time instancetl ∈ [0, T ] a new mesh Tl,h depending on the previous grid Tl−1,h through the followingh-adaptation operations: cutting (C) one element Ik into Ik1 and Ik2 and gluing (G)two elements Ik1 and Ik2 together into Ik. The described adaptation process has tocomply with restrictions on a minimal admissible size of mesh step hmin, a maximaladmissible size of mesh step hmax, a maximal number of elements Nmax and keepingof local quasi-uniformity (5) and barrier consistency (6), respectively.The main idea of the proposed h-adaptive strategy is based on• mesh refinement in domains with irregular solution (low regularity) or withhigh value of residual estimate,• mesh coarsening in domains with solution of high regularity and low value ofresidual estimate.The estimation of the regularity of the solution is essential for mesh refinement. Thepresented approach is based on a measure of inter-element jumps arising from theshock capturing techniques in hyperbolic problems, for a survey see [2].We have employed the following element-wise regularity indicator:gIk(uh) :=1h2pk+1kk∑i=k−1[uh(xi)]2 , k = 1, . . . N, (15)which recognizes the subdomains of Ω where the solution is smooth (gIk ≈ 0) fromthe areas with discontinuities or with a very steep gradient (gIk ≫ 1).The second key ingredient of the mesh refinement is the residual estimator whichis chosen proportionally to the strong formulation of the local residue from [1] asrIk(uh) :=∂uh∂t−12σ2x2∂2uh∂x2− rx∂uh∂x+ ruh, Ik ∈ Th. (16)Then the local and global residual estimators of the approximate solution uh aredefined by resIk(uh) := ‖rIk‖L2(Ik) and resG(uh) :=√∑Ik∈Thres2Ik , respectively.Our interest is to find a solution ũh ∈ Shp such that resG(ũh) ≤ TOL, whereTOL > 0 is a given tolerance. In order to satisfy this condition we prescribe thefollowing stopping criterion for the h-adaptivity: resIk ≤TOLN, ∀ Ik ∈ Th, whichguarantees the uniform distribution of the global residue.The whole h-adaptation DG algorithm can be schematically written as971. let TOL > 0, 0 < hmin ≤ hmax and Nmax be given,2. let B−, B+ ←→ T0h and Shp be set up, let u0 ←→ u0h be given,3. repeat time loop (until tl > T ) (l = 1, . . . , r)(a) solve problem (13)–(14) on Tl−1,h =⇒ ulh,(b) evaluate indicators gIk(ulh), resIk(ulh), ∀ Ik ∈ Tl−1,h =⇒ resG(ulh),(c) if resG(ulh) > TOL⇒ h-refinement,(C) h-refine elements with resIk >TOLN,(G) h-derefine elements with resIk < δTOLN∧ gIk(ulh) ≈ 0,(•) construct new mesh T newh −→ Tl−1,h and space Shp, go to (a),(d) if resG(ulh) ≤TOLβ⇒ h-coarsening,{(G) h-derefine elements with resIk < δTOLN∧ gIk(ulh) ≈ 0,(•) construct new mesh T newh −→ Tl−1,h and space Shp, go to (a),where β > 1 and δ ∈ (0, 1) are user-defined parameters, in our computations theyare typically chosen as β = 3.0 and δ = 0.1.5. Numerical exampleThe presented numerical example represents the case of a discrete double barriercall option with the expiration date T = 812(e.g. 8 months) and the strike priceK = 6.0. The prescribed barriers are B− = 4.0, B+ = 8.0 and computationaldomain was set as Ω = [0, 9]. The Black-Scholes model parameters were the risk-freeinterest rate r = 1.0y−1 and volatility σ2 = 0.01y−1. The initial uniform mesh withspatial step h = 0.25 was adaptively refined according to h-adaptation parametershmin = 10−3 and hmax = 0.5. The time step is τ =1120and we consider monthlymonitoring. We carried out computations by piecewise quadratic approximations,set Θ = 0 and used the restarted GMRES for the solving of linear systems (13).Table 1 illustrates the development of the global residue and the number of ele-ments during the computation in comparison with an adapted and uniform mesh.One can easily observe that for approximately the same values of the global residue,it is sufficient to use less elements in the adapted case than for the uniform one.Figure 1 shows the corresponding isolines of option price and global residue in space-time plot with well-resolved monthly monitoring.6. ConclusionWe have dealt with the numerical solution of the discrete barrier option pricingmodels, represented by the linear convection-diffusion-reaction equation. We havepresented DG approach together with simple h-adaptivity technique. Presented nu-merical example illustrated the potency of the resulting scheme.AcknowledgementsThis work was supported by the ESF Project No. CZ.1.07/2.3.00/09.0155.98time (bimonthly) resG (adapted) #Tlh∗ resG (uniform) #Tlh0.000000 19.960869 36 10.888578 1200.166667 1.498133 178 1.459563 1200.333333 0.615153 58 0.476494 1200.500000 0.572154 44 0.475957 1200.666667 0.119596 58 0.124287 120Table 1: Comparison of h-adaptive and uniform approach w.r.t. resG; Tlh∗ denotesinput meshes after monitoring without the updated h-refinement or h-coarsening.0.10.20.30.40.50.63 4 5 6 7 8 9time (t)Stock Price (S)0.10.20.30.40.50.63 4 5 6 7 8 9time (t)Stock Price (S)Figure 1: The isolines of price u (left) and corresponding global residue resG (right).References[1] Achdou Y., and Pironneau, O.: Computational methods for option pricing.Society for Industrial and Applied Mathematics, Philadelphia, 2005.[2] Doleǰśı, V., Feistauer, M., and Schwab C.: On some aspects of the discontinuousGalerkin finite element method for conservation laws. Math. Comput. Simulation61 (2003), 333–346.[3] Feistauer, M., Felcman, J., and Straškraba, I.: Mathematical and computationalmethods for compressible flow. Oxford University Press, Oxford, 2003.[4] Hozman, J.: Discontinuous Galerkin method for convection-diffusion problems.Ph.D. thesis, Charles University Prague, 2009.[5] Riviére, B.: Discontinuous Galerkin methods for solving elliptic and parabolicequations: theory and implementation. Frontiers in Applied Mathematics. Societyfor Industrial and Applied Mathematics, Philadelphia, 2008.[6] Seydel, R.: Tools for computational finance: 4th edition. Springer, Berlin, 2008.99

Valuing barrier options using the adaptive discontinuous Galerkin method

http://dml.cz/bitstream/handle/10338.dmlcz/702712/PANM_16-2012-1_17.pdf

Valuing barrier options using the adaptive discontinuous Galerkin method

Abstract

Similar works

Full text

Available Versions

Institute of Mathematics AS CR, v. v. i.