A computational weighted finite difference method for American and barrier options in subdiffusive Black-Scholes model
aa r X i v : . [ q -f i n . C P ] D ec A computational weighted finite di ff erence method for American and barrieroptions in subdi ff usive Black-Scholes model Grzegorz Krzy˙zanowski , Marcin Magdziarz Hugo Steinhaus Center, Faculty of Pure and Applied Mathematics, Wroclaw University of Science and Technology 50-370 Wroclaw, Poland
Abstract
Subdi ff usion is a well established phenomenon in physics. In this paper we apply the subdi ff usive dynamics toanalyze financial markets. We focus on the financial aspect of time fractional di ff usion model with moving boundaryi.e. American and barrier option pricing in the subdi ff usive Black-Scholes (B-S) model. Two computational methodsfor valuing American options in the considered model are proposed - the weighted finite di ff erence (FD) and theLongsta ff -Schwartz method. In the article it is also shown how to valuate numerically wide range of barrier optionsusing the FD approach. Keywords:
Weighted finite di ff erence method, subdi ff usion, time fractional Black-Scholes model, American optionnumerical evaluation.
1. Introduction
Subdi ff usion is the standard tool in the analysis of complex systems. In physics it is manifested by the celebratedFractional Fokker-Planck equation. This equation was derived from the continuous-time random walk scheme withheavy-tailed waiting times. Equivalently it can be described in terms of subordination, where the standard di ff usionprocess is time-scaled by the so-called inverse subordinator. Nowadays the subdi ff usion is a very well establishedphenomenon with many real-life examples - e.g. di ff usion in percolative and porous systems, charge carrier transportin amorphous semiconductors, transport on fractal geometries. (see [31, 40] and references in them)In this paper we apply the subdi ff usive dynamics to analyze financial markets. Option pricing is the core content ofmodern finance and has fundamental meaning for global economy. By the recent announcement of Futures IndustryAssociation trading activity in the global exchange-traded derivatives markets, in 2019 reached a record of 34 , ,
23 billion of them were options contracts [13]. The value of the global derivatives marketsis estimated 700 trillion dollars to upwards of 1 , ff ective tool of the optionsvaluation. The model was of such great importance that the authors were awarded the Nobel Prize for Economicsin 1997. The classical model was generalized in order to weaken its strict assumptions, allowing such features as Email addresses: [email protected] (Grzegorz Krzy˙zanowski), [email protected] (MarcinMagdziarz)
Preprint submitted to Elsevier January 1, 2021 tochastic interest [34], regime switching model [24], stochastic volatility [17], market regulations [1], and transac-tions costs [26]. Analysis of empirical financial records indicates that the data can exhibit fat tails. This feature hasbeen observed in many di ff erent markets (see e.g. [5] and the references therein). Such dynamics can be observed inemerging markets where the number of sellers and buyers is low. Also an interest rate often exhibits the feature ofconstant periods appearing - e.g. in US between 2002 and 2017 [23]. In response to the empirical evidence of fat tails, α -stable distribution as an alternative to the Gaussian law was proposed [12, 32]. The stable distribution has foundmany important applications, for example in finance [41], physics [15, 19, 35] and electrical engineering [47]. Withincreasing interest of fractional calculus and non-local di ff erential operators the family of fractional B-S equations hasemerged in the recent literature (see e.g. [22] and references therein).The subdi ff usive B-S is the generalization of the classical B-S model to the cases, where the underlying assetsdisplay characteristic periods in which they stay motionless. The standard B-S model does not take this phenomenainto account because it assumes the asset is described by continuous Gaussian random walk. As a result of an optionpricing for such underlying asset, the fair price provided by the B-S model is misestimated. In order to describe thisdynamics properly, the subdi ff usive B-S model assumes that the underlying instrument is driven by α -stable inversesubordinator [45]. Then the frequency of the constant periods appearing is dependent of subdi ff usion parameter α ∈ (0 , α →
1, the subdifussive B-S is reducing to the classical model. Due to its practicality and simplicity, thestandard B-S model is one of the most widely used in option pricing. Although in contrast to the subdi ff usive case itdoes not take into account the empirical property of the constant price periods in the underlying instrument dynamics.In Figure 1 we compare sample simulation of underlying asset in classical and subdi ff usive market model. Even shortstagnation of a market can not be simulated by standard B-S model. As a generalization of the classical B-S model,its subdi ff usive equivalent can be used in wide range of markets - including all cases where B-S can be applied.In the case of standard B-S many numerical methods for American and barrier options have been proposed. Theyinclude such methods and techniques as: linear complementarity [6, 21, 42, 50], Bermudan approximation [4, 42],penalty method [14, 33, 39], repeated Richardson extrapolation [3, 10], finite di ff erence [21, 38, 52] / element [29, 33]and binomial / trinomial methods [43, 44]. Also numerous analytical and semi analytical methods have been developedwhich are based mainly on asymptotic approximation formulas, method of lines, fast Fourier transform and integralequations. Despite their high numerical e ffi ciency they can not be adapted to many applications e.g. the standard B-Swith variable volatility (see e.g. [3] and references therein).Since the subdi ff usive B-S model was proposed [30] many open problems still have remained unsolved. Oneof them is the way of valuation American and Exotic options. In this paper we derive the Linear ComplementarityProblem (LCP) system describing the fair price of an American option in subdi ff usive B-S model. We apply theweighted scheme of the FD method and the Longsta ff -Schwartz method to solve the system numerically. We compareboth methods, moreover we show how to valuate numerially wide range of barrier options in subdi ff usive B-S modelusing the FD approach. The paper is an extension of [22] where the governing fractional di ff erential equation for theEuropean option in the subdi ff usive B-S is derived and the weighted finite di ff erence method is proposed. Moreoverin [22] the optimal value of weighting parameter θ for the wide class of fractional di ff usion-type problems is found.
2. Subdi ff usive B-S model The evolution of the market is taking place up to time horizon T and is contained in the probability space ( Ω , F , P ).Here, Ω is the sample space, F is a filtration interpreted as the information about history of asset price and P is the“objective” probability measure. The assumptions are the same as in the classical B-S model [21] with the exceptionthat we do not have to assume the market liquidity and that the underlying instrument, instead of Geometric BrownianMotion (GBM), follows subdi ff usive GBM [30]: ( Z α ( t ) = Z ( S α ( t )) , Z (0) = Z , where Z α ( t ) - the price of the underlying instrument, Z ( t ) = Z (0) exp ( µ t + σ B t ), µ - drift (constant), σ - volatility(constant), B t - Brownian motion, S α ( t ) - the inverse α -stable subordinator defined as S α ( t ) = inf( τ > U α ( τ ) > t )[30], where U α ( t ) is the α -stable subordinator [45], 0 < α <
1. Here S α ( t ) is independent of B t for each t ∈ [0 , T ]. InFigure 1 we compare the samples trajectories of GBM and subdifussive GBM for given α -stable subordinator.2 Figure 1: The sample trajectory of GBM (upper panel) with its subdi ff usive analogue (middle panel) and the corresponding inversesubordinator (lower panel). In the subdi ff usive GBM the constant periods characteristic for emerging markets can be observed.The parameters are Z = σ = µ = α = . Let us introduce the probability measure Q ( A ) = Z A exp − γ B ( S α ( T )) − γ S α ( T ) ! dP , (1)where γ = ( µ + σ ) /σ , A ∈ F . As it is shown in [30], Z α ( t ) is a Q -martingale. The subdi ff usive B-S model isarbitrage-free and incomplete [30]. Despite Q defined in (1) is not unique, but it is the “best” martingale measure inthe sense of criterion of minimal relative entropy. It means that the measure Q minimizes the distance to the measure P [31]. Between European put and call options the put-call parity holds [30].
3. Selected options
In Tables 1 and 2 we recall the payo ff functions for options considered in this article. Recall that the payo ff function f ( Z t ) is the gain of the option holder at the time t for given underlying instrument Z . European AmericanPlain max ( Z T − K ,
0) max ( Z t − K , Z T − K , { M T > H + } max ( Z t − K , { M t > H + } Knock up-and-out max ( Z T − K , { M T < H + } max ( Z t − K , { M t < H + } Knock down-and-in max ( Z T − K , { m T < H − } max ( Z t − K , { m t < H − } Knock down-and-out max ( Z T − K , { m T > H − } max ( Z t − K , { m t > H − } Knock double-out max ( Z T − K , { H + > M T , m T > H − } max ( Z t − K , { H + > M t , m t > H − } Knock double-in max ( Z T − K , − max ( Z T − K , { H + > M T , m T > H − } max ( Z t − K , − max ( Z t − K , { H + > M t , m t > H − } Table 1: Payo ff functions for selected call options. uropean AmericanPlain max ( K − Z T ,
0) max ( K − Z t , K − Z T , { M T > H + } max ( K − Z t , { M t > H + } Knock up-and-out max ( K − Z T , { M T < H + } max ( K − Z t , { M t < H + } Knock down-and-in max ( K − Z T , { m T < H − } max ( K − Z t , { m t < H − } Knock down-and-out max ( K − Z T , { m T > H − } max ( K − Z t , { m t > H − } Knock double-out max ( K − Z T , { H + > M T , m T > H − } max ( K − Z t , { H + > M t , m t > H − } Knock double-in max ( K − Z T , − max ( K − Z T , { H + > M T , m T > H − } max ( K − Z t , − max ( K − Z t , { H + > M t , m t > H − } Table 2: Payo ff functions for selected put options. Here and in the rest of the paper K - strike, Z t = Z α ( t ) - value of underlying instrument at time t , t ∈ [0 , T ], M t = max τ ∈ [0 , t ] ( Z τ ), m t = min τ ∈ [0 , t ] ( Z τ ), H + , H − - barriers.
4. Valuation of American option as free boundary problem
In whole paper we assume that the dividend rate δ =
0. The next proposition explains why in context of Americanoptions we will proceed only with the put options.
Proposition 4.1.
If the dividend rate δ = r ≥
0, then the value of American call option is equal to its Europeananalogue. Similarly it can be shown that if r =
0, then it is not worth to realize American put before T , so in this casevalue of American put is equal to his European equivalent. Proof:. [21] Let us assume δ = r ≥ t < T .It will be done only if it would be profitable, it means if Z t > K . Then, the expected profit (i.e. the future expectedprofit computed at time t =
0) is ( Z t − K ) e − rt . Option holder can also proceed following strategy: at time t short sellan underlying instrument for Z t , next buy it (using the American or European option) for min ( Z T , K ) at time T . Thenhis expected profit is Z t e − rt − min ( Z T , K ) e − rT . Let us observe that the following sequence of inequalities holds: Z t e − rt − min ( Z T , K ) e − rT ≥ ( Z t − min ( Z T , K )) e − rt ≥ ( Z t − K ) e − rt . Hence, earlier realization of the American option is at most as profitable as use of the option at t = T . If the gainof American option holder in the most profitable strategy is the same as the profit of European option holder, weconclude that both instruments have the same value. Indeed, the fair price of an option is equal to the expected gainof the most profitable strategy provided by the option. Analogously it can be shown the second statement. (cid:4) We proceed with the following main result of this section
Theorem 4.1.
The fair price of an American put option in the subdi ff usive B-S model is equal to v ( z , t ), where v ( z , t )satisfies: x = ln z , u ( x , t ) = v ( e x , T − t ) . (2)and u ( x , t ) is the solution of the system u ( x , = max (cid:0) K − exp ( x ) , (cid:1) , u ( x , t ) ≥ max (cid:0) K − exp ( x ) , (cid:1) , t ∈ [0 , T ] c D α u ( x , t ) − σ ∂ u ( x , t ) ∂ x − r − σ ! ∂ u ( x , t ) ∂ x + ru ( x , t ) ≥ t ∈ (0 , T ) (cid:0) u ( x , t ) − max (cid:0) K − exp ( x ) , (cid:1)(cid:1) · c D α u ( x , t ) − σ ∂ u ( x , t ) ∂ x − r − σ ! ∂ u ( x , t ) ∂ x + ru ( x , t ) ! = t ∈ (0 , T )lim x →∞ u ( x , t ) = , lim x →−∞ u ( x , t ) = K , (3)4here x ∈ ( −∞ , ∞ ) and c D α t is Caputo fractional derivative defined as [48]: c D α t g ( t ) = Γ (1 − α ) Z t dg ( s ) ds ( t − s ) − α ds . Proof:.
We consider the subdi ff usive B-S Equation [22], c D α t v ( z , t ) = − σ z ∂ v ( z , t ) ∂ z − rz ∂ v ( z , t ) ∂ z + rv ( z , t ) , ( z , t ) ∈ (0 , ∞ ) × (0 , T ) , (4)with the terminal condition determining put option v ( z , T ) = max ( K − z ,
0) , z ∈ (0 , ∞ ) . (5)At the time t ∈ [0 , T ] we can gain at least max ( K − z ,
0) (by exercising the option) and maybe even more. It leads usto the inequality: v ( z , t ) ≥ max ( K − z , . (6)After the optimal exercise moment v ( z , t ) can not describe the value process of the optimal strategy (because we”missed” the best moment to use the option). Hence, we can gain at most the value provided by keeping the option -i.e. by the solution of (4). In other words, the di ff erential equation (4) can underestimate the value of the option. Sofor t ∈ (0 , T ) true is the following inequality: c D α t v ( z , t ) + σ z ∂ v ( z , t ) ∂ z ≤ rv ( z , t ) − rz ∂ v ( z , t ) ∂ z . (7)Indeed, the di ff erential inequality (7) cannot formally be interpreted that over any infinitesimal time interval the sumof the losses and profits from holding the option is lower or equal to the return at the riskless rate. Note that thetime derivative from (7) can be interpreted as loss in value due to having less time for exercising the option and thederivative of second order over the price of the underlying instrument can be interpreted as the profit in holding theoption.At each moment t ∈ (0 , T ) we decide if it is worth to use the option, or keep it. Mathematically we can describe it as: v ( z , t ) = max ( K − z , , if we realize option or c D α t v ( z , t ) + σ z ∂ v ( z , t ) ∂ z + rz ∂ v ( z , t ) ∂ z − rv ( z , t ) = , if we keep it on. This can be written as follows:( v ( z , t ) − max ( K − z , · c D α t v ( z , t ) + σ z ∂ v ( z , t ) ∂ z + rz ∂ v ( z , t ) ∂ z − rv ( z , t ) ! = . (8)For the su ffi ciently high price of the underlying instrument we will not use the option. So:lim z →∞ v ( z , t ) = . (9)In analogy if the price of the underlying instrument will be low we will use the option, selling the underlying instru-ment for K : lim z → v ( z , t ) = K . (10)After the change of variables (2), (5)-(10) have the form of (3). (cid:4) . Numerical scheme for American put option In this section we derive the numerical scheme for American put option. To do so, we will approximate limits byfinite numbers and derivatives by finite di ff erences. We will proceed for a θ -convex combination of explicit ( θ = θ =
0) discrete scheme, similarly as it was done for European options in [22]. In the same workthe corresponding stability / convergence analysis can be found. We introduce parameter θ ∈ [0 ,
1] by optimizationpurposes - similarly as for the case α = θ = has the best properties in terms of the error and unconditionalstability / convergence [21]. Instead of the continuous space ( −∞ , ∞ ) × [0 , T ], we take its discrete and finite equivalent { x = x min , x . . . , x n = x max } × { t = , . . . , t N = T } , where x min , x max - lower and upper boundary of the grid. Weconsider a uniform grid, so t j = j ∆ t and x i = x min + i ∆ x , where i = , . . . , n , j = , . . . , N , ∆ t = T / N , ∆ x = ( x max − x min ) / n . After obtaining the discrete analogue of (3) we will solve it recursively. As a result we will find itsnumerical solution ˆ u ji = ˆ u ( x i , t j ), i = , . . . , n , j = , . . . , N .We begin discretizing initial condition for put option (first line of (3)), we getˆ u i = max (cid:0) K − exp ( x min + i ∆ x ) , (cid:1) , (11)for i = , . . . , n .Similarly the discrete version of boundary conditions (the last 2 lines of (3)) has the form ˆ u ln = , ˆ u l = K , (12)for l = , . . . , N .Discretizing third, second and fourth lines of (3) we get C ˆ u = ˆ u + (1 − θ ) G + θ G + θ B ˆ u , ˆ u = max( ˆ u , ˆ u ) , C ˆ u k + = k − X j = (cid:16) b j − b j + (cid:17) ˆ u k − j + b k ˆ u + (1 − θ ) G k + + θ G k + θ B ˆ u k , ˆ u k + = max (cid:16) ˆ u k + , ˆ u (cid:17) , (13)for k ≥ b j = ( j + − α − j − α , C = ( θ I + (1 − θ ) A ) , A = (cid:16) a i j (cid:17) ( n − × ( n − , such that: a i j = + ad ∆ x + cd , for j = i , i = , ..., n − , − ad ∆ x − bd ∆ x ! , for j = i − , i = ..., n − , − ad ∆ x + bd ∆ x ! , for j = i + , i = ..., n − , , in other cases, B = (cid:16) b i j (cid:17) ( n − × ( n − , such that: b i j = − ad ∆ x + cd ! , for j = i , i = ..., n − , ad ∆ x + bd ∆ x , for j = i + , i = , ..., n − , ad ∆ x − bd ∆ x , for j = i − , i = , ..., n − , , in other cases, G k = ad ∆ x − bd ∆ x ! ˆ u k , , ..., , ad ∆ x + bd ∆ x ! ˆ u kn − ! T , k = (cid:16) u k , u k , ..., u kn − (cid:17) T , a = σ , b = r − σ ! , c = r , d = Γ (2 − α ) ∆ t α , ∆ t = T / N , k = , , ... N . Note that the analogical scheme for the European option [22] is C ˆ u = ˆ u + (1 − θ ) G + θ G + θ B ˆ u , C ˆ u k + = k − X j = (cid:16) b j − b j + (cid:17) ˆ u k − j + b k ˆ u + (1 − θ ) G k + + θ G k + θ B ˆ u k , (14)with corresponding boundary conditions ˆ u ln = exp( x max ) − K exp( − r ( T − t l )) , ˆ u l = , (15)and initial condition for a call option ˆ u i = max (cid:0) exp ( x min + i ∆ x ) − K , (cid:1) , (16)where l = , . . . , N , i = , . . . , n . Figure 2: The price of American put ( P A ) in dependence of K . The parameters are n = x min = − x max = σ = r = , N = T = Z = θ = ˇ θ α .
6. Numerical schemes for barrier options
The systems (13) and (14) can be used to price di ff erent types of options, only if the initial-boundary conditionswill be properly modified. Let us treat x min and x max not as approximations of infinite values, but as logarithm of lower7nd supreme barriers H − and H + defined in double barrier option. We take a logarithm because of change of variables x = ln z made in (2) and in [22]. The initialˆ u i = max (cid:0) exp (cid:0) ln H − + i ∆ x (cid:1) − K , (cid:1) , (17)and boundary conditions ˆ u ln = , ˆ u l = , (18)where l = , . . . , N , ∆ x = (ln H + − ln H − ) / n , i = , . . . , n , together with (14) is the scheme for the European double knock-out call option. The same boundary-initial conditionswith (13) create the scheme for the American double knock-out call option. Analogously, prices of one-side barrierknock-out options can be obtained. Hence we have initial and boundary conditions for knock-down-and-out calloption ˆ u i = max (cid:0) exp (cid:0) ln H − + i ∆ x (cid:1) − K , (cid:1) , (19) ˆ u ln = exp( x max ) − K exp( − r ( T − t l )) , ˆ u l = , (20) l = , . . . , N , ∆ x = ( x max − ln H − ) / n , i = , . . . , n ,and for knock-up-and-out call option ˆ u i = max (cid:0) exp ( x min + i ∆ x ) − K , (cid:1) , (21) ˆ u l = , ˆ u ln = , (22) l = , . . . , N , ∆ x = (ln H + − x min ) / n , i = , . . . , n .If we want to price the knock-in options, it is helpful to use the fact that for fixed parameters there holds the socalled in-out parity Knock in = Van − Knock out , where Van - the price of Vanilla (plain) option,
Knock in , Knock out - option prices of knock-in and knock-out of thesame type and style.Please note, that the value of a double knock-out option for Z outside of the interval (ln H − , ln H + ) (but being apositive number) is equal 0. Analogous remark applies for one-sided barrier options.To summarize we present the way to price the considered options in Tables 3, 4 and 5.8 Figure 3: The price of European down-and-out call option ( C d − oE ) in dependence of Z . The parameters are n = σ = , r = , N = x max = H − = T = K = θ = ˇ θ α . Style of option Numerical scheme Boundary conditions Initial condition Apply in-out parity?Plain (14) (15) (16) NoKnock up-and-in (14) (22) (21) YesKnock up-and-out (14) (22) (21) NoKnock down-and-in (14) (20) (19) YesKnock down-and-out (14) (20) (19) NoKnock double-out (14) (18) (17) NoKnock double-in (14) (18) (17) Yes
Table 3: European call options.
To price the European put options we can firstly compute their call equivalents and then apply the Put-Call parity.We can also use other initial conditions than in Table 3,
Style of option Numerical scheme Boundary conditions Initial condition Apply in-out parity?Plain (14) (12) (11) NoKnock up-and-in (14) (27) (26) YesKnock up-and-out (14) (27) (26) NoKnock down-and-in (14) (25) (24) YesKnock down-and-out (14) (25) (24) NoKnock double-out (14) (18) (23) NoKnock double-in (14) (18) (23) Yes
Table 4: European put options. where ˆ u i = max (cid:0) K − exp (cid:0) ln H − + i ∆ x (cid:1) , (cid:1) , (23)for ∆ x = (ln H + − ln H − ) / n (double options),ˆ u i = max (cid:0) K − exp (cid:0) ln H − + i ∆ x (cid:1) , (cid:1) , (24)9
02 4 6 81 3 5 7 91.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5024130.51.52.53.5
Figure 4: The price of European double knock-out call option ( C − oE ) in dependence of Z . The parameters are n = σ = , r = , N = H + = H − = T = K = θ =
0. In this case, there is no clear relation between the prices of option fordi ff erent α (like e.g. in vanilla equivalent for T > α the price is higher for all Z ). We can conclude that in thisfigure there is a critical point where all plots intersect. It is unknown under which conditions (if there is any) such point exists andwhat is its value. ˆ u ln = , ˆ u l = , (25)for l = , . . . , N , ∆ x = ( x max − ln H − ) / n , i = , . . . , n , (knock-down options),ˆ u i = max (cid:0) K − exp ( x min + i ∆ x ) , (cid:1) , (26) ˆ u l = K , ˆ u ln = , (27)for l = , . . . , N , ∆ x = (ln H + − x min ) / n , i = , . . . , n (knock-up options).If there is no dividend (the case we consider in this paper), the American call is equal its European equivalent soTable 3 holds also for American call options. For American put options we have, Style of option Numerical scheme Boundary conditions Initial condition Apply in-out parity?Plain (13) (12) (11) NoKnock up-and-in (13) (27) (26) YesKnock up-and-out (13) (27) (26) NoKnock down-and-in (13) (25) (24) YesKnock down-and-out (13) (25) (24) NoKnock double-out (13) (18) (23) NoKnock double-in (13) (18) (23) Yes
Table 5: American put options. ∆ x is di ff erent.In Figures 2, 3,4 we compare the fair prices given for di ff erent values of α for American put, European down-and-out call and European double knock-out call option respectively.We recall the observation from [22] that the optimal choice of θ for given α is such that ˇ θ α = − − α − − α . Thenthe lowest boundary for an error is achieved without losing the unconditional stability / convergence. For θ > ˇ θ α thestability / convergence is not provided. In Figure 5 the relation between the fair price of American put P A and θ ispresented. The real price of the option is close to 0 ,
25. The jump presented in the figure is the result of the increasingerror. It is the consequence of lack of the stability [22].
Figure 5: The price of American put ( P A ) in dependence of θ . The explosion of the numerical error is the consequence of the lackof the unconditional stability outside of the interval [0 , ˇ θ α ] [22]. The parameters are n = σ = , r = , N = T = Z = K = α = , x min = − x max = ff -Schwartz method The Longsta ff -Schwartz (LS) method is one of the most popular approaches for valuing American / Bermudanoptions and their Asian equivalents [20, 28]. Moreover it has an important applications in solving dynamic investmentportfolio problems and in American / Bermuda style swaptions valuation (see e.g. [27] and references therein). Allthese applications have an important meaning in finance. Only for notional amount of interest rate swaps outstandingat the end of 1999 the losses caused by wrong exercise strategies were estimated on billions of dollars [27].The main idea of this method is the use of least squares to estimate the conditional expected payo ff to the optionholder from continuation. This strategy allows to find the value of an option with the optimal exercise time (i.e.the moment where exercising option is the most profitable, if there are more than one such moments we choose thelowest of them). The method was introduced for the classical B-S model but it can be extended for many other cases[28]. In this paper we focus on American option, but note that the same method can be used to price Bermudan andAmerican-style Asian options. The LS algorithm can be found in Appendix.Note that the inverse- α stable process is not Markovian so the expected value in A . A . T the method is unstable. The reason is the ill-condition of the underlying regression problem [36]. The analyticalformulas indicating regime where the method is stable and where is not are still unknown. This fact limitates the11ange of possible applications of the algorithm. Note that the error in LS method could be of di ff erent origin. Thefirst is produced by discretization the continuous stochastic processes into m nodes and assumption that the Americanoption can be exercise only at these points. The second is related to Monte Carlo method i.e. that we estimate theexpected value by the mean of size M . The next possible origin of the error is coming from approximation of theconditional expected value A . l basis functions. The last possible type of the error is producedfor non-Markovian underlying processes. Since these processes have memory, the expected value A . x min and x max ), therefore its stability / convergence is easier to analyse. Example 1.
We compare both methods presented in this paper in pricing American put option. Simulations are madefor r = , Z = K = N = x max = x min = − n = σ = m = M = ff erentvalues of T and α . Figure 6: The price of American put option computed by FD and LS for di ff erent T and α (upper panel) with corresponding runningtime of the algorithms (lower panel). For small values of T , LS does not match the real solution properly. FD is precise and fastmethod for all T and α . We take θ = ˇ θ α . FDLS
Figure 7: The price of American put option computed by FD and LS for di ff erent α (upper panels) and corresponding running timeof the algorithms (lower panels). For small values of T (left panels), LS is working visibly slower and is producing higher error as α is decreasing. For bigger T (right panels), LS is matching the FD output. For this case, increasing m and M follows approachingLS to the FD output, but also increases time of computation. FD is precise and fast method for all α and T . We take θ = Figure 6 presents a comparison between methods FD and LS for 3 particular choices of α . We see, that as T increases, the results of both methods are closer (because LS algorithm is based on the Monte Carlo, oscillations arevisible). As α gets higher, the LS result is closer to the FD output. Also the time of computation is visibly higher forlower values of T and α . Note that for T close to 0 option prices manifested by di ff erent parameters α can significantlydi ff er - like in the left end of the time interval T = ,
01 (see the values computed by FD method). For each α ∈ (0 , T → P A for “small” and “big” parameter T . Also the relation between running time of bothmethods t run and α is provided. In both figures the FD is the reference method. The LS method requires to generate(and save) M paths of Z α and that is memory / speed expensive. As α decreases, more time for generating Z α ( t ) isrequired. For t < ES α ( t ) > t , so the dynamics of Z α ( t ) is ”faster” than its classical equivalent ( t ∈ [0 , T ]). Thus,the errors caused by approximating the non-Markovian process by Markovian approach can accumulate with errorscaused by ill-condition of the regression problem. Even for ”big” T the LS method is running visibly longer and haslower precision than the FD. The advantage of LS is finding sample of optimal exercise times, which is not providedby the FD. We can conclude that the presented LS algorithm could not work correctly for ”small” T (in particular if T < α the method is ine ffi cient. The use of the method should be considered only if there is a need ofoptimal exercise strategy. For computing the fair price of American option the FD approach is recommended. Example 2.
Let us take parameters T = Z = σ = , r = , K = H − = x max = θ = ˇ θ α [22]). Note that the implicit FD method considered in [51] is theparticular case θ = α = n and N . Note that for α =
1, the optimal scheme for the parameter θ ( ˇ θ α = ,
5) iscalled the Crank-Nicolson (C-N) method. For calculating the exact price of the option we use the analytical formula[37]: C d − o = Z Φ ( y ) − Ke − rT Φ ( y − σ √ T ) − Z ( H − / Z ) λσ − + Φ ( y ) + Ke − rT ( H − / Z ) λσ − Φ ( y − σ √ T ), where Φ denotes cumulative distribution function of the standard normal distribution, λ = r − σ / y = ln( Z / K ) + ( r + σ / T σ √ T , y = ln(( H − ) / ( KZ )) + ( r + σ / T σ √ T . The C-N method is more e ffi cient than its implicit equivalent. The same conclusion is true13or American options (see e.g. [21]). ( n , N ) (20 ,
20) (40 ,
40) (100 , , , , θ = ,
98% 1 ,
03% 0 ,
39% 0 ,
18% 0 ,
06% 0 , θ = ˇ θ α ,
55% 0 ,
28% 0 ,
07% 0 ,
02% 0 0
Table 6: Relative error for knock down-and-out call option. The real value is 0 , ( n , N ) (20 ,
20) (40 ,
40) (100 , , , , θ = , × ,
02 0 ,
16 0 ,
56 4 ,
35 43 , θ = ˇ θ α , × ,
02 0 ,
15 0 ,
56 4 ,
34 43 , Table 7: The average running time [ s ] of the implicit and optimal FD method related to Table 6. In Table 8 there are presented the relative errors (columns 6-9) and the di ff erences between the option pricescalculated for implicit and optimal numerical schemes (columns 3-5) for di ff erent n , N and α . The most significantdi ff erence is observed for high α , since ˇ θ α (column 2) is increasing function of α . As α approaches 0, the optimalscheme is reducing to the implicit method. The di ff erence for n = N = α is lowerthan 10 − . Since the di ff erences between both prices for n = N = θ = θ = ˇ θ α respectively. In Table 7 and Table 9 the corresponding (mean) running times to Table 6 andTable 8 respectively are presented. The corresponding values to ( n , N ) = (1500 , n and N ) the algorithm works longer. In contrast to Table 6 and Table8 the significant di ff erences between both numerical schemes were not observed. C d − o ( ˇ θ α ) − C d − o (0) relative error θ = θ = ˇ θ α α ˇ θ α ( n , N ) = (20 ,
20) (200 , , ,
40) (100 , ,
40) (100 , , ,
48 0 ,
65% 0 ,
07% 0 ,
01% 1 ,
01% 0 ,
39% 0 ,
36% 0 , , ,
46 0 ,
51% 0 ,
05% 0 0 ,
91% 0 ,
35% 0 ,
36% 0 , , ,
44 0 ,
39% 0 ,
04% 0 0 ,
78% 0 ,
31% 0 ,
33% 0 , , ,
41 0 ,
29% 0 ,
03% 0 0 ,
64% 0 ,
26% 0 ,
28% 0 , , ,
37 0 ,
20% 0 ,
02% 0 0 ,
5% 0 ,
22% 0 ,
23% 0 , , ,
33 0 ,
13% 0 ,
01% 0 0 ,
36% 0 ,
18% 0 ,
17% 0 , , ,
27 0 ,
07% 0 0 0 ,
22% 0 ,
15% 0 ,
11% 0 , Table 8: The di ff erences and relative errors for the option prices calculated by optimal and implicit numerical schemes for di ff erent n , N and α . α θ = θ = ˇ θ α ( n , N ) = (20 ,
20) (40 ,
40) (100 , , , ,
20) (40 ,
40) (100 , , , , ,
01 0 ,
04 0 ,
24 0 ,
51 45 ,
64 0 ,
01 0 ,
04 0 ,
24 0 ,
51 45 , , ,
02 0 ,
05 0 ,
28 0 ,
51 47 ,
53 0 ,
02 0 ,
05 0 ,
28 0 ,
51 46 , , ,
02 0 ,
03 0 ,
15 0 ,
54 43 ,
34 0 ,
02 0 ,
03 0 ,
15 0 ,
53 42 , , ,
02 0 ,
04 0 ,
14 0 ,
54 43 ,
61 0 ,
02 0 ,
04 0 ,
14 0 ,
54 43 , , ,
02 0 ,
04 0 ,
14 0 ,
74 43 ,
49 0 ,
02 0 ,
04 0 ,
14 0 ,
74 43 , , ,
02 0 ,
03 0 ,
14 0 ,
72 45 , ,
02 0 ,
03 0 ,
14 0 ,
72 45 , , ,
02 0 ,
04 0 ,
14 0 ,
52 45 ,
42 0 ,
02 0 ,
04 0 ,
14 0 ,
52 46 , Table 9: The average running time [ s ] of the implicit and optimal FD method related to Table 8. . Summary In this paper:– We have derived the system describing the fair price of American put option in subdi ff usive B-S model.– We have introduced the weighted numerical scheme for this system.– We have shown how to modify previous results for valuing wide range of barrier options in frame of the samemodel.– We have given the formula for the optimal choice of weighting parameter θ in dependence of subdi ff usionparameter α .– We have applied the Longsta ff -Schwartz numerical method for the subdi ff usive B-S model. This method isworse than the FD method in terms of speed and precision of computation. Moreover for small values of T thismethod is unstable so can not be used in many di ff erent cases.– We have presented some numerical examples to illustrate introduced theory.The numerical techniques presented in this paper can successfully be repeated for other time fractional di ff usionmodels with moving boundaries problems. Acknowledgements
This research was partially supported by NCN Sonata Bis 9 grant nr 2019 / / E / ST1 / Appendix A. Longsta ff -Schwartz Algorithm The method uses a dynamic programming to find the optimal stopping time and Monte Carlo to approximate thefair price of an option. We start assuming the exercise time τ is equal T . Going backwards to 0, we replace τ by eachmoment we find where is better to exercise. Let us denote V ( Z ( t ) , t ) as the fair price of an American option with thepayo ff function f ( Z ( t )) and underlying asset Z ( t ). It is easy to conclude that V ( Z , = E (cid:0) e − r τ f ( Z ( τ )) (cid:1) . Let us divide the interval [0 , T ] into m subintervals (of the same length) using the grid [ t = , t , . . . , t m = T ],moreover we introduce H ( Z ( t i )) = E (cid:16) e − r ( τ i − t ) f ( Z ( τ i )) | Z ( t i ) (cid:17) , (A.1)where τ i is the optimal exercise moment in { t i + , . . . , t m − , t m } , i = , . . . , m −
2. The interpretation of the function H ( Z ( t i )) is the expected profit from keeping the option up to time t i . For each trajectory we will proceed using thefollowing algorithm:In other words at each grid point from { t , . . . , t m − } we compare profit from keeping and exercise the option. Thenwe decide is it more profitable to exercise the option or keep it on. In the algorithm above V = V ( Z (0) , H ( Z ( t i )) for i = , . . . , m . To do so, Longsta ff and Schwartz proposed to useleast squares regression. This can be done since the conditional expectation is an element of L space, so it can berepresented using its infinite countable orthonormal basis. To proceed with the computations, the finite set of l suchbasis elements should be chosen. For the simulations we choose 3 first elements of Laguerre polynomials L ( x ) = L ( x ) = − x , L ( x ) = / (cid:16) − x − x (cid:17) . Note that the early exercise at t can be profitable only if f ( Z ( t )) >
0, i.e. ifoption is in the money. The whole LS algorithm for the subdi ff usive case will look as follows:15 lgorithm 1 τ = t m , V = f ( Z ( τ )) for t from t m − to t do V ← e − r ∆ t V if H ( Z ( t )) < f ( Z ( t )) then τ ← t V ← f ( Z ( τ )) end if end for V ← e − r ∆ t V Return V . Algorithm 2 Generate Z j ( t i ) [30] for j = . . . , M , i = . . . , m τ = [ t m , t m , . . . , t m ], V = f ( Z ( τ )) for t from t m − to t do Find in the money trajectories i.e. w = { j , . . . , j R } s.t. f ( Z k ( t )) > k ∈ w Put Z w ← [ Z j , . . . , Z j R ], V w ← [ V j , . . . , V j R ] Find regression coe ffi cients β , . . . , β l such that l X i = β i L i ( Z w ) = e − r ∆ t V w , For k ∈ w if l X i = β i L i ( Z k ) < f ( Z k ( t )) then τ k ← t V k ← f ( Z ( τ k )) end if for i ∈ { . . . , M } \ w do V i ← e − r ∆ t V i end for end for Price ← M X i = e − r ∆ t V i M Return
Price . 16 eferences [1] Hazhir Aliahmadi, Mahsan Tavakoli-Kakhki, and Hamid Khaloozadeh. Option pricing under finite moment log stable process in a regulatedmarket: a generalized fractional path integral formulation and monte carlo based simulation.
Communications in Nonlinear Science andNumerical Simulation , page 105345, 2020.[2] Ghada Alobaidi, Roland Mallier, and A Stanley Deakin. Laplace transforms and installment options.
Mathematical Models and Methods inApplied Sciences , 14(08):1167–1189, 2004.[3] Luca Vincenzo Ballestra. Fast and accurate calculation of American option prices.
Decisions in Economics and Finance , 41(2):399–426,2018.[4] Luca Vincenzo Ballestra and Liliana Cecere. A fast numerical method to price American options under the bates model.
Computers & Mathematics with Applications , 72(5):1305–1319, 2016.[5] Szymon Borak, Adam Misiorek, and Rafał Weron. Models for heavy-tailed asset returns. In
Statistical tools for finance and insurance ,pages 21–55. Springer, 2011.[6] Michael J Brennan and Eduardo S Schwartz. The valuation of American put options.
The Journal of Finance , 32(2):449–462, 1977.[7] Paul Brockman and Harry J Turtle. A barrier option framework for corporate security valuation.
Journal of Financial Economics , 67(3):511–529, 2003.[8] Zhongdi Cen and Wenting Chen. A HODIE finite di ff erence scheme for pricing American options. Advances in Di ff erence Equations ff erence scheme for pricing American put options with Singularity-Separating method. Numerical Algorithms , 53.4: 497-510, 2010.[10] Chuang-Chang Chang, San-Lin Chung, and Richard C Stapleton. Richardson extrapolation techniques for the pricing of American-styleoptions.
Journal of Futures Markets: Futures, Options, and Other Derivative Products , 27(8):791–817, 2007.[11] Bernd Engelmann, Matthias R Fengler, Morten Nalholm, and Peter Schwendner. Static versus dynamic hedges: an empirical comparisonfor barrier options.
Review of Derivatives Research , 9(3):239–264, 2006.[12] Eugene F Fama. Risk, return and equilibrium: some clarifying comments.
The Journal of Finance , 23(1):29–40, 1968.[13] FIA. Global Futures and Options Trading Reaches Record Level in 2019, 2019.[14] Peter A Forsyth and Kenneth R Vetzal. Quadratic convergence for valuing American options using a penalty method.
SIAM Journal onScientific Computing , 23(6):2095–2122, 2002.[15] V Yu Gonchar, AV Chechkin, EL Sorokovoi, VV Chechkin, LI Grigor’eva, and ED Volkov. Stable L´evy distributions of the density andpotential fluctuations in the edge plasma of the u-3m torsatron.
Plasma Physics Reports , 29(5):380–390, 2003.[16] Steve Heston and Zhou Guofu. On the rate of convergence of discrete-time contingent claims.
Mathematical Finance , 10.1: 53-75, 2000.[17] John Hull and Alan White. The pricing of options on assets with stochastic volatilities.
The journal of finance , 42(2):281–300, 1987.[18] Patrick Jaillet, Damien Lamberton, and Bernard Lapeyre. Variational inequalities and the pricing of American options.
Acta ApplicandaeMathematica , 21.3 (1990): 263-289.[19] Aleksander Janicki and Aleksander Weron. Can one see α -stable variables and processes? Statistical Science , pages 109–126, 1994.[20] Ralf Korn, Elke Korn, and Gerald Kroisandt.
Monte Carlo methods and models in finance and insurance . CRC press, 2010.[21] Grzegorz Krzy˙zanowski. Selected applications of di ff erential equations in Vanilla Options valuation. Mathematica Applicanda , 46(2), 2018.[22] Grzegorz Krzy˙zanowski, Marcin Magdziarz, and Łukasz Płociniczak. A weighted finite di ff erence method for subdi ff usive Black–Scholesmodel. Computers & Mathematics with Applications , 80(5):653 – 670, 2020.[23] Grzegorz Krzy˙zanowski, Ernesto Mordecki, and Andr´es Sosa. A zero interest rate Black-Derman-Toy model. arXiv preprintarXiv:1908.04401 , 2019.[24] Sha Lin and Xin-Jiang He. A regime switching fractional Black–Scholes model and european option pricing.
Communications in NonlinearScience and Numerical Simulation , 85:105222, 2020.[25] Yumin Lin and Chuanju Xu. Finite di ff erence / spectral approximations for the time-fractional di ff usion equation. Journal of computationalphysics
Computers & Mathematics with Applications , 65(11):1719–1726, 2013.[27] Francis A Longsta ff , Pedro Santa-Clara, and Eduardo S Schwartz. Throwing away a billion dollars: The cost of suboptimal exercisestrategies in the swaptions market. Journal of Financial Economics , 62(1):39–66, 2001.[28] Francis A Longsta ff and Eduardo S Schwartz. Valuing American options by simulation: a simple least-squares approach. The review offinancial studies , 14(1):113–147, 2001.[29] Sofiane Madi, Mohamed Cherif Bouras, Mohamed Haiour, and Andreas Stahel. Pricing of American options, using the Brennan–Schwartzalgorithm based on finite elements.
Applied Mathematics and Computation , 339:846–852, 2018.[30] Marcin Magdziarz. Black-Scholes formula in subdi ff usive regime. Journal of Statistical Physics , 136(3):553–564, 2009.[31] Marcin Magdziarz and Janusz Gajda. Anomalous dynamics of Black-Scholes model time-changed by inverse subordinators.
Acta PhysicaPolonica B , 43(5), 2012.[32] Benoit B Mandelbrot. The variation of certain speculative prices. In
Fractals and scaling in finance , pages 371–418. Springer, 1997.[33] SAJID Memon. Finite element method for American option pricing: a penalty approach.
International Journal of Numerical Analysis andModelling, Series B: Computing and Information , 3(3):345–370, 2012.[34] Robert C Merton. On the pricing of corporate debt: The risk structure of interest rates.
The Journal of finance , 29(2):449–470, 1974.[35] T Mizuuchi, VV Chechkin, K Ohashi, EL Sorokovoy, AV Chechkin, V Yu Gonchar, K Takahashi, S Kobayashi, K Nagasaki, H Okada, et al.Edge fluctuation studies in heliotron j.
Journal of nuclear materials , 337:332–336, 2005.[36] Oleksii Mostovyi. On the stability the least squares monte carlo.
Optimization Letters , 7(2):259–265, 2013.[37] Marek Musiela and Marek Rutkowski. Martingale methods in financial modelling, 2005, 2005.
38] JC Ndogmo and DB Ntwiga. High Order Accurate Implicit Methods for Barrier Option Pricing.
Applied Mathematics and Computation ,2011.[39] Bjørn Fredrik Nielsen, Ola Skavhaug, and Aslak Tveito. Penalty and front-fixing methods for the numerical solution of American optionproblems.
Journal of Computational Finance , 5(4):69–98, 2002.[40] Sebastian Orzeł and Aleksander Weron. Calibration of the subdi ff usive Black-Scholes model. Acta Phys. Pol. B , 41(5):1051–1059, 2010.[41] Svetlozar Rachev and Stefan Mittnik. Stable Paretian models in finance.
Vol. 7. Wiley
Applied Mathematics and Computation , 251:363–377, 2015.[43] LCG Rogers and EJ Stapleton. Fast accurate binomial pricing.
Finance and Stochastics , 2(1):3–17, 1997.[44] Mark Rubinstein. On the relation between binomial and trinomial option pricing models.
The Journal of Derivatives , 8(2):47–50, 2000.[45] Ken-iti Sato, Sato Ken-Iti, and A Katok.
L´evy processes and infinitely divisible distributions . Cambridge university press, 1999.[46] Aleksandra Stankovska. Global derivatives market.
SEEU Review , 12, 01 2016.[47] Bart W Stuck and B Kleiner. A statistical analysis of telephone noise.
Bell System Technical Journal , 53(7):1263–1320, 1974.[48] Dina Tavares, Ricardo Almeida, and Delfim FM Torres. Caputo derivatives of fractional variable order: numerical approximations.
Com-munications in Nonlinear Science and Numerical Simulation , 35:69–87, 2016.[49] Aleksander Weron and Rafał Weron. In˙zynieria finansowa (financial engineering).
Wydawnictwo Naukowo-Techniczne, Warszawa , 1998.[50] Paul Wilmott.
Derivatives: The theory and practice of financial engineering . John Wiley & Son Limited, 1998.[51] Hongmei Zhang, et al. Numerical solution of the time fractional Black–Scholes model governing European options.
Computers & Mathe-matics with Applications , 71.9 (2016): 1772-1783.[52] Jichao Zhao, Matt Davison, and Robert M Corless. Compact finite di ff erence method for American option pricing. Journal of Computationaland Applied Mathematics , 206(1):306–321, 2007., 206(1):306–321, 2007.