A Numerical Method for Pricing Discrete Double Barrier Option by Lagrange Interpolation on Jacobi Node
aa r X i v : . [ q -f i n . C P ] F e b A Numerical Method for Pricing Discrete Double Barrier Option by LagrangeInterpolation on Jacobi Node
Amirhossein Sobhani a, ∗ , Mariyan Milev b a Department of Applied Mathematics, Amirkabir University of Technology, No. 424, Hafez Ave., Tehran 15914, Iran b UFT-PLOVDIV, Department of Mathematics and Physics
Abstract
In this paper, a rapid and high accurate numerical method for pricing discrete single and double barrier knock-out calloptions is presented. According to the well-known Black-Scholes framework, the price of option in each monitoringdate could be calculate by computing a recursive integral formula upon the heat equation solution. We have approx-imated these recursive solutions with the aim of Lagrange interpolation on Jacobi polynomials node. After that, anoperational matrix, that makes our computation significantly fast, has been driven. The most important feature of thismethod is that its CPU time dose not increase when the number of monitoring dates increases. The numerical resultsconfirm the accuracy and e ffi ciency of the presented numerical algorithm. Keywords:
Double and single barrier options, Black-Scholes model, Option pricing, Jacobi polynomials
1. Introduction
Barrier options play a key role in financial markets where the most important problem is the so called optionvaluation problem, i.e. to compute a fair value for the option, i.e. the premium. The Nobel Prize-winning Black-Scholes option valuation theory motivates using classical numerical methods for partial di ff erential equations (PDE’s)[1]. In computational Finance numerous nonstandard numerical methods are proposed and successfully applied forpricing options [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. Numerical methods are often preferred to closed-form solutionsas it they could me more easily extended or adapted to satisfy all the financial requirements of the option contractsand continuously changing conditions imposed by financial institutions and over-the-counter market for controllingtrading of derivatives.Kunitomo and Ikeda [13] obtained general pricing formulas for European double barrier options with curvedbarriers but like for a variety of path-dependent options and corporate securities most formulas are obtained forrestricted cases as continuous monitoring or single barrier [5]. The discrete monitoring is essential as the trading yearis considered to consist of 250 working days and a week of 5 days. Thus, taking for one year T =
1, the applicationof barriers occurs with a time increment of 0 .
004 daily and 0 .
02 weekly.For discrete barrier options there are some analytical solutions. For example, Fusai reduces the problem of pricingone barrier option to a Wiener-Hopf integral equation [3]. Several other di ff erent contracts with discrete time mon-itoring are characterized by updating the initial conditions, such as Parisian options and occupation time derivatives[14]. We remark that although most real contracts specify fixed times for monitoring the asset, academic researchershave focused mainly on continuous time monitoring models as the analysis of fixed barriers could be treated mathe-matically using some techniques such as the reflection principle [15]. For example, using the reflection principle inBrownian motions, Li expresses the solution in general as summation of an infinite number of normal distributionfunctions for standard double barrier options, and in many non-trivial cases the solution consists of finite terms [16]. ∗ Corresponding author
Email addresses: [email protected], [email protected] (Amirhossein Sobhani), [email protected] (Mariyan Milev)
Preprint submitted to Elsevier February 5, 2018 elsser derives a formula for continuous double barrier knock-out and knock-in options by inverting analytically theLaplace transform by a contour integration, [17]. Broadie et. al. have found an explicit correction formula for dis-cretely monitored option with one barrier [18]. However, these three well-known methods [6, 10, 11] have not beenstill applied in the presence of two barriers, i.e. a discrete double barrier option.Although it could not be claimed that it is impossible to be found an exact or closed-form solution of the Black-Scholes equation [19] for the valuation of discrete double barrier knock-out call option, it is sure that there is asubstantial di ff erences in the option prices between continuous and discrete monitoring even for 1 000 000 monitoringdates. This could be trivially tested for a single barrier knock-in and knock-out option using formulas [3], [13][6],or the correction formula [18], for double barrier knock-out options with the numerical algorithm [5] or with a high-order accurate finite di ff erence scheme [11]. It is well-known in literature the relation when comparing the price ofcontinuous and discretely monitored barrier options with the corresponding vanilla option with same parameters andabsence of rebates. The discrete monitoring considerably complicates the analysis of barrier options [18] and theirpricing often requires nonstandard method as those presented in [2, 5, 7, 11]. Di ffi culties of pricing double barrieroptions emerge even in the case of continuous monitoring where some drawbacks of close-form formulas could beclearly observed. The analytical solutions of such options is usually expressed as infinite series of reflections andpresented with Fourier series. For fixed barriers contracts the Fourier series solution gives the same answer when allthe terms have been added up but the main drawback is that the rate of convergence of the sum to the solution can bequite di ff erent, depending on the time to expiry.Initially classical quantitative methods in Finance have been explored for pricing barrier options. This includesstandard lattice techniques, i.e. the binomial and trinomial trees of Kamrad and Ritchken [20], Boyle and Lau [21],Kwok [15], Heyen and Kat [22], Tian [23], Dai and Lyuu [24] used standard lattice techniques, the binomial andtrinomial trees, for pricing barrier options. Ahn et al. [25] introduce the adaptive mesh model (AMM) that increasesthe e ffi ciency of trinomial lattices. The Monte Carlo simulation methods were implemented in [26, 27, 28, 29, 30, 31].Also numerical algorithms based on quadrature methods have been proposed in [32, 5].Recently a great variety of more sophisticated semi-analytical methods for pricing barrier options have beendeveloped which are based on integral transforms [3, 33, 34], or on the transition probability density function of theprocess used to describe the underlying asset price [32, 5, 35, 18, 36, 4, 37, 38]. Farnoosh et al. [39, 40] have proposeda numerical algorithms for pricing discrete single and double barrier options with time-dependent parameters, whilein [41] a projection methods have been explored. These techniques are very high performing for pricing discretelymonitored single and double barrier options and our computational results are in very good agreement with them. Themain objective of this paper is present a new e ffi cient computational method for valuation of discrete barrier optionsbased on a Lagrange interpolation on Jacobi nodes that have not only a simpler computer implementation but alsodi ff er with minimum memory requirements and extreme short computational times.This article is organized as follows. In Section 2 we formulate the mathematical model for valuation of barrieroptions under the classical Black-Scholes framework. In Section 3 we briefly list definitions for Jacobi Polynomials.In section 4 we propose a new e ffi cient numerical methods where an orthogonal Lagrange interpolation is utilizedand a suitable operational matrix form has been obtained for pricing discrete double barrier options. One of the mainadvantages of this algorithm is that it do not depend on the number of monitoring dates. In the next Section 5 weobserve numerical errors of order 10 − and 10 − in maximum norm for di ff erent computational experiments accordingto the number of node points. The obtained results are in good agreement with other benchmark values in literatureand this confirms the e ffi ciency and accuracy of the presented numerical algorithm.
2. The Pricing Model
We assume that the stock price process S t follows the Geometric Brownian motion: dS t S t = rdt + σ dB t where S , r and σ are initial stock price, risk-free rate and volatility respectively. We consider the problem of pricingknock-out discrete double barrier call option, i.e. a call option that becomes worthless if the stock price touches eitherlower or upper barrier at the predetermined monitoring dates:0 = t < t < · · · < t M = T .
2e assume that that monitoring dates are equally spaced, i.e; t m = m τ where τ = TM . If the barriers are not touchedin monitoring dates, the pay o ff at maturity time is max( S T − E , E is exercise price. The price of optionis defined discounted expectation of pay o ff at the maturity time. Based on the Black-Scholes framework, the optionprice P ( S , t , m −
1) as a function of stock price at time t ∈ ( t m − , t m ), satisfies in the following partial di ff erentialequations − ∂ P ∂ t + rS ∂ P ∂ S + σ S ∂ P ∂ S − r P = , (1)subject to the initial conditions: P ( S , t , = ( S − E ) (max( E , L ) ≤ S ≤ U ) P ( S , t m , = P ( S , t m , m − ( L ≤ S ≤ U ) ; m = , , ..., M − , where P ( S , t m , m −
1) : = lim t → t m P ( S , t , m − E ∗ = ln (cid:16) EL (cid:17) ; µ = r − σ ; θ = ln (cid:16) UL (cid:17) and δ = max { E ∗ , } , we define g m ( z ) as following recursiveformula: g ( z ) = Z θ k ( z − ξ, τ )g ( ξ ) d ξ (2)g m ( z ) = Z θ k ( z − ξ, τ )g m − ( ξ ) d ξ ; m = , , ..., M (3)where g ( z ) = Le − α z (cid:16) e z − e E ∗ (cid:17) ( δ ≤ z ≤ θ ) , (4) k ( z , t ) = √ π c t e − z c t . (5)It could be shown that the price of the knock-out discrete double barrier option can be obtain as follows ( see [41] ): P ( S , t M , M − ≃ e α z + β t g M ( z /θ ) (6)where z = log (cid:16) S L (cid:17) .
3. Jacobi Polynomials
Let w ( α,β ) ( x ) = (1 − x ) α (1 + x ) β , α, β > − , and L w ( α,β ) ( − ,
1) be Hilbert space with the following inner product andnorm: < f , g > w ( α,β ) = Z − f ( x ) g ( x ) w ( x ) dx , (7) k f k w ( α,β ) = p < f , f > w ( α,β ) . (8)The Jacobi polynomials, J ( α,β ) i ( x ) are orthogonal polynomials in L w ( α,β ) ( − , Z − J ( α,β ) i ( x ) J ( α,β ) j ( x ) w ( x ) dx = λ i δ i j , (9)where λ i = k J ( α,β ) i k . These polynomials, that set an orthogonal basis in L w ( α,β ) ( − , J ( α,β )0 ( x ) = , J ( α,β )1 ( x ) =
12 ( α + β + x +
12 ( α − β ) (10) J ( α,β ) i + ( x ) = (cid:16) a ( α,β ) i x − b ( α,β ) i J ( α,β ) i ( x ) (cid:17) − c ( α,β ) i J ( α,β ) i − ( x ) (11)3here: a ( α,β ) i = (2 i + α + β + i + α + β + i + n + α + β +
1) (12) b ( α,β ) i = ( β − α )(2 n + α + β + i + n + α + β + n + α + β ) (13) c ( α,β ) i = ( n + α )( n + β )(2 n + α + β + i + n + α + β + n + α + β ) . (14)
4. Pricing by orthogonal Lagrange interpolation
In this section we consider Π n as space of all polynomials with degree less or equal to n , set points { x α,β i } ni = asroots of ( n + J ( α,β ) n + that are shifted to [0 , θ ] and I α,β n : C [0 , θ ] → Π n as orthogonal polynomialinterpolation projection operator, that is defined as follows: I α,β n ( f ) = n X i = f ( x α,β i ) L i ( x ) (15)where L i ( x ) is the i -th Lagrange polynomial basis function defined on { x α,β i } ni = : L i ( x ) = n Y j = , j , i ( x − x α,β j )( x α,β i − x α,β j ) . (16)Let operator K : L ([0 , θ ]) → L ([0 , θ ]) is defined as follows: K (cid:0) g (cid:1) ( z ) : = Z θ κ ( z − ξ, τ )g( ξ ) d ξ. (17)where κ is defined in (5). According to the definition of operator K , equations (2) and (3) can be rewritten as below:g = K g (18)g m = K g m − m = , , ..., M (19)We denote ˜g , n = I α,β n K (cid:0) g (cid:1) (20)˜g m , n = I α,β n K (cid:0) ˜g m − (cid:1) = (cid:16) I α,β n K (cid:17) m (cid:0) g (cid:1) , m ≥ . (21)where I α,β n K is as follows: ( I α,β n K )(g) = I α,β n (cid:0) K (g) (cid:1) . Since, ˜g m , n ∈ Π n for m ≥
1, we can write ˜g m , n = n X i = a mi L i ( z ) = Φ ′ n ( x ) G m , where G m = [ a m , a m , · · · , a mn ] ′ and Φ n = [ L m , L m , · · · , L n ] ′ . From equation (21) we obtain˜g m , n = ( I α,β n K ) m − (cid:0) ˜g , n (cid:1) . (22)4 β − . − . − . . e −
06 8 . e −
06 2 . e −
05 1 . e −
05 9 . e − − . . e −
06 7 . e −
06 1 . e −
05 3 . e −
05 4 . e −
050 2 . e −
05 2 . e −
05 2 . e −
05 5 . e −
05 9 . e − . e −
05 9 . e −
05 9 . e −
05 8 . e −
05 1 . e − . e −
04 1 . e −
04 1 . e −
04 1 . e −
04 1 . e − Table 1: The maximum norm error for n =
25 of example(1) with L =
95 and M = Since Π n is a finite dimensional linear space, thus the linear operator I α,β n K on Π n could be considered as a n × n matrix K . Consequently equation (22) can be written as following matrix operator form˜g m , n = Φ ′ n K m − G . (23)For evaluation of the option price by (23), it is enough to calculate the matrix operator K and the vector G . It iseasy to check (see [41]) that: G = [ a , a , · · · , a n ] ′ K = (cid:16) k i j (cid:17) n × n where a i = Z θδ κ ( x α,β i − ξ, τ )g ( ξ ) d ξ , ≤ i ≤ n . k i j = Z θ κ ( x α,β i − ξ, τ ) L j − ( ξ ) d ξ . Therefore, the price of the knock-out discrete double barrier option can be estimated as follows: P ( S , t M , M − ≃ e α z + β t ˜g M , n ( z ) (24)where z = log (cid:16) S L (cid:17) and ˜g M , n from (23). The matrix form of relation (23) implies that the computational time ofpresented algorithm be nearly fixed when monitoring dates increase. Actually, the complexity of our algorithm is O ( n ) that dose not depend on number of monitoring dates.
5. Numerical Result
In the current section, the presented method in previous section for pricing knock-out call discrete double barrieroption is compared with some other methods. The numerical results are obtained from the relation (24) with n basisfunctions. The Source code has been written in M atlab Example 1.
In the first example, the pricing of knock-out call discrete double barrier option is considered with thefollowing parameters: r = . , σ = . , T = . , S = , E = , U = and L = , , , , . . Intable (2) , numerical results of presented method with Milev numerical algorithm [5], Crank-Nicholson [42], trinomial,adaptive mesh model (AMM) and quadrature method QUAD-K200 as benchmark [43] are compared for variousnumber of monitoring dates. In addition, it can be seen that CPU time of presented method is fixed against increasesof monitoring dates. l L PM ( α = − . , β = − . =
25) Milev(200) Milev(400) Trinomial AMM-8 Benchmark80 2.4499 - - 2.4439 2.4499 2.449990 2.2028 - - 2.2717 2.2027 2.20285 95 1.6831 1.6831 1.6831 1.6926 1.6830 1.683199 1.0811 1.0811 1.0811 0.3153 1.0811 1.081199.9 0.9432 0.9432 0.9432 - 0.9433 0.9432CPU 0.035 s 1 s 5 s80 1.9420 - - 1.9490 1.9419 1.942090 1.5354 - - 1.5630 1.5353 1.535425 95 0.8668 0.8668 0.8668 0.8823 0.8668 0.866899 0.2931 0.2931 0.2931 0.3153 0.2932 0.293199.9 0.2023 0.2023 0.2023 - 0.2024 0.2023CPU 0.035 s 8 s 30 s80 1.6808 - - 1.7477 1.6807 1.680890 1.2029 - - 1.2370 1.2028 1.2029125 95 0.5532 0.5528 0.5531 0.5699 0.5531 0.553299 0.1042 0.1042 0.1042 0.1201 0.1043 0.104299.9 0.0513 0.0513 0.0513 - 0.0513 0.0513CPU 0.035 s 35 s 150 sTable 2: Double barrier option pricing of Example (1): T = . r = . σ = . S = E = Number of nodal points l og ( M a x - e rr o r) -25-20-15-10-50 ( α , β )=(0,0)( α , β )=(0.5,0.5)( α , β )=(-0.5,-0.5) (a) M = Number of nodal points l og ( M a x - e rr o r) -20-15-10-505 ( α , β )=(0,0)( α , β )=(0.5,0.5)( α , β )=(-0.5,-0.5) (b) M = Max − error for example (1) with L = Example 2.
In this example, the parameters of knock-out call discrete double barrier option is considered as r = . , σ = . , T = . , E = , U = and L = . In table (3) the option price for di ff erent spot prices are evaluatedand compared with Milev numerical algorithm [5], Crank-Nicholson [42] and the Monte Carlo (MC) method with paths [44].
80 85 90 95 100 105 110 115 120 e rr o r × -6 -8-6-4-202468 ( α , β )=(-0.5,-0.5) (a) Error S
80 85 90 95 100 105 110 115 120 E s t i m a t ed P r i c e ( α , β )=(-0.5,-0.5) (b) Estimated PriceFigure 2: The error and estimated Price in example(1) with L =
80 and M = S PM ( α = − . , β = − . , )( n =
25) Crank-Nicolson(1000) Milev(400) Milev(1000) MC (st.error)with 10 paths95 0.174498 0.1656 0.174503 0.174498 -95.0001 0.174499 ≃ ≃ Table 3: Double barrier option pricing of Example (2): T = . M = r = . σ = . E = U =
110 and L = Example 3.
Due to the fact that the probability of crossing upper barrier during option’s life when U ≥ E is toosmall, the price of discrete single down-and-out call option can be estimated by double ones by setting upper barriergreater than E (for more details see[5]). Now, we consider a discrete single down-and-out call option with thefollowing parameters: r = . , σ = . , T = . , S = , E = and L = , . , . . The price is estimatedby double ones with U = . E. The numerical results are shown in table (4) and compared with Fusai’s analyticalformula [3], the Markov chain method (MCh)[2] and the Monte Carlo method (MC) with paths [29] that showsthe validity of presented method in this case.
7M ( α = − . , β = − . =
25 n =
50 (IR17) MCH MC (st.error)95 25 6.63104 6.63156 6.63156 6.6307 6.63204 (0.0009)99.5 25 3.35644 3.35558 3.35558 3.3552 3.35584 (0.00068)99.9 25 3.00897 3.00887 3.00887 3.0095 3.00918 (0.00064)95 125 6.16940 6.16863 6.16864 6.1678 6.16879 (0.00088)99.5 125 1.95811 1.96130 1.96130 1.9617 1.96142 (0.00053)99.9 125 1.50991 1.51020 1.51068 1.5138 1.5105 (0.00046)CPU 0.038 s 0.051 s
Table 4: Single barrier option pricing of Example (3): T = . r = . σ = . S = E = U = Example 4.
In this example we estimate the price of continus monitoring call barrier down and out option, P c, withdiscrete ones, P d m , using the following formula[18]: P c ( L ) = P d m (cid:18) L e βσ √ ∆ t (cid:19) , where β = ζ (1 / / √ π ≃ . with ζ the Riemann zeta function. The parameters of this problem is considered asr = . , σ = . , T = . , E = , S = . In table (5) the option price for di ff erent Lower barriers are evaluatedand compared with continuous monitoring price that is obtained in [18]. As we can see, this estimations is accurateexcept when the barrier is close to the spot price. PM( α = − . , β = − . , M =
50) PM( α = − . , β = − . , M = =
25 n =
50 n =
25 n = Table 5: Single barrier option pricing with continuous monitoring of Example (4): T = . r = . σ = . S = E = U =
6. Conclusion and remarks
In this article, we used the Lagrange interpolation on Jacobi polynomial nodes for pricing discrete single anddouble barrier options. In section 4 we obtained a matrix relation (23) for solving this problem. Numerical resultsverify that computational time is fixed when the number of monitoring dates increase.
ReferencesReferences [1] G. D. Smith, Numerical solution of partial di ff erential equations: finite di ff erence methods, Oxford university press, 1985.[2] J.-C. Duan, E. Dudley, G. Gauthier, J.-G. Simonato, Pricing discretely monitored barrier options by a markov chain, The Journal of Derivatives10 (4) (2003) 9–31.[3] G. Fusai, I. D. Abrahams, C. Sgarra, An exact analytical solution for discrete barrier options, Finance and Stochastics 10 (1) (2006) 1–26. doi:10.1007/s00780-005-0170-y .[4] G. Fusai, M. C. Recchioni, Analysis of quadrature methods for pricing discrete barrier options, Journal of Economic Dynamics and Control31 (3) (2007) 826–860.
5] M. Milev, A. Tagliani, Numerical valuation of discrete double barrier options, Journal of Computational and Applied Mathematics 233 (10)(2010) 2468–2480.[6] A. Tagliani, M. Milev, Laplace transform and finite di ff erence methods for the black–scholes equation, Applied Mathematics and Computation220 (2013) 649–658.[7] H. Gzyl, M. Milev, A. Tagliani, Discontinuous payo ff option pricing by mellin transform: A probabilistic approach, Finance Research Letters20 (2017) 281–288.[8] A. Sobhani, M. Milev, A numerical method for pricing discrete double barrier option by legendre multiwavelet, Journal of Computationaland Applied Mathematics 328 (2018) 355–364.[9] M. J. Dilloo, D. Y. Tangman, A high-order finite di ff erence method for option valuation, Computers & Mathematics with Applications.[10] M. Milev, A. Tagliani, E ffi cient implicit scheme with positivity preserving and smoothing properties, Journal of Computational and AppliedMathematics 243 (2013) 1–9.[11] J. Ndogmo, D. Ntwiga, High-order accurate implicit methods for the pricing of barrier options, arXiv preprint arXiv:0710.0069.[12] S. Kabaivanov, V. Markovska, Modelling environment changes for pricing weather derivatives, Scientific Annals of Economics and Business64 (4) (2017) 423–430.[13] N. Kunitomo, M. Ikeda, Pricing options with curved boundaries, Mathematical finance 2 (4) (1992) 275–298.[14] G. Fusai, S. Sanfelici, A. Tagliani, Practical problems in the numerical solution of pde’s in finance.[15] Y. K. Kwok, Mathematical models of financial derivatives. 1998.[16] A. Li, The pricing of double barrier options and their variations.[17] A. Pelsser, Pricing double barrier options using laplace transforms, Finance and Stochastics 4 (1) (2000) 95–104.[18] M. Broadie, P. Glasserman, S. Kou, A continuity correction for discrete barrier options, Mathematical Finance 7 (4) (1997) 325–349.[19] O. E. Barndor ff -Nielsen, Processes of normal inverse gaussian type, Finance and stochastics 2 (1) (1997) 41–68.[20] B. Kamrad, P. Ritchken, Multinomial approximating models for options with k state variables, Management science 37 (12) (1991) 1640–1652.[21] P. P. Boyle, S. H. Lau, Bumping up against the barrier with the binomial method, The Journal of Derivatives 1 (4) (1994) 6–14.[22] R. C. Heynen, H. M. Kat, Barrier options, Exotic options: the State of the Art (1997) 125–159.[23] Y. S. Tian, Pricing complex barrier options under general di ff usion processes, The Journal of Derivatives 7 (2) (1999) 11–30.[24] T.-S. Dai, Y.-D. Lyuu, The bino-trinomial tree: A simple model for e ffi cient and accurate option pricing, The Journal of Derivatives 17 (4)(2010) 7–24.[25] D.-H. Ahn, S. Figlewski, B. Gao, Pricing discrete barrier options with an adaptive mesh model, Available at SSRN 162450.[26] L. Andersen, R. Brotherton-Ratcli ff e, Exact exotics, Risk 9 (10) (1996) 85–89.[27] D. R. Beaglehole, P. H. Dybvig, G. Zhou, Going to extremes: Correcting simulation bias in exotic option valuation, Financial AnalystsJournal (1997) 62–68.[28] P. Baldi, L. Caramellino, M. G. Iovino, Pricing general barrier options: a numerical approach using sharp large deviations, MathematicalFinance 9 (4) (1999) 293–321.[29] M. Bertoldi, M. Bianchetti, Monte carlo simulation of discrete barrier options, Financial engineering-Derivatives Modelling, Caboto SIMSpa, Banca Intesa Group, Milan, Italy 25 (2003) 1.[30] G. C. Kuan, N. Webber, Pricing barrier options with one-factor interest rate models, The Journal of Derivatives 10 (4) (2003) 33–50.[31] D. Jeong, M. Yoo, C. Yoo, J. Kim, A hybrid monte carlo and finite di ff erence method for option pricing, Computational Economics (2017)1–14.[32] A. D. Andricopoulos, M. Widdicks, P. W. Duck, D. P. Newton, Universal option valuation using quadrature methods, Journal of FinancialEconomics 67 (3) (2003) 447–471.[33] M. Broadie, Y. Yamamoto, A double-exponential fast gauss transform algorithm for pricing discrete path-dependent options, OperationsResearch 53 (5) (2005) 764–779.[34] F. Fang, C. W. Oosterlee, Pricing early-exercise and discrete barrier options by fourier-cosine series expansions, Numerische Mathematik114 (1) (2009) 27.[35] A. Golbabai, L. Ballestra, D. Ahmadian, A highly accurate finite element method to price discrete double barrier options, ComputationalEconomics 44 (2) (2014) 153–173.[36] G. Dorfleitner, P. Schneider, K. Hawlitschek, A. Buch, Pricing options with green’s functions when volatility, interest rate and barriers dependon time, Quantitative finance 8 (2) (2008) 119–133.[37] C. Skaug, A. Naess, Fast and accurate pricing of discretely monitored barrier options by numerical path integration, Computational Economics30 (2) (2007) 143–151.[38] M. A. Sullivan, Pricing discretely monitored barrier options, Journal of computational finance 3 (4) (2000) 35–52.[39] R. Farnoosh, H. Rezazadeh, A. Sobhani, M. H. Beheshti, A numerical method for discrete single barrier option pricing with time-dependentparameters, Computational Economics 48 (1) (2015) 131–145. doi:10.1007/s10614-015-9506-7 .[40] R. Farnoosh, A. Sobhani, H. Rezazadeh, M. H. Beheshti, Numerical method for discrete double barrier option pricing with time-dependentparameters, Computers & Mathematics with Applications 70 (8) (2015) 2006–2013. doi:10.1016/j.camwa.2015.08.016 .[41] R. Farnoosh, A. Sobhani, M. H. Beheshti, E ffi cient and fast numerical method for pricing discrete double barrier option by projection method,Computers & Mathematics with Applications.[42] B. Wade, A. Khaliq, M. Yousuf, J. Vigo-Aguiar, R. Deininger, On smoothing of the crank–nicolson scheme and higher order schemes forpricing barrier options, Journal of Computational and Applied Mathematics 204 (1) (2007) 144–158.[43] C.-J. Shea, Numerical valuation of discrete barrier options with the adaptive mesh model and other competing techniques, Master’s Thesis,Department of Computer Science and Information Engineering, National Taiwan University.[44] P. Brandimarte, Numerical methods in finance: a MATLAB-based introduction, Vol. 489, John Wiley & Sons, 2003.cient and fast numerical method for pricing discrete double barrier option by projection method,Computers & Mathematics with Applications.[42] B. Wade, A. Khaliq, M. Yousuf, J. Vigo-Aguiar, R. Deininger, On smoothing of the crank–nicolson scheme and higher order schemes forpricing barrier options, Journal of Computational and Applied Mathematics 204 (1) (2007) 144–158.[43] C.-J. Shea, Numerical valuation of discrete barrier options with the adaptive mesh model and other competing techniques, Master’s Thesis,Department of Computer Science and Information Engineering, National Taiwan University.[44] P. Brandimarte, Numerical methods in finance: a MATLAB-based introduction, Vol. 489, John Wiley & Sons, 2003.