Efficient BEM-Based Algorithm for Pricing Floating Strike Asian Barrier Options ( with MATLAB R © Code )

This paper aims to illustrate how SABO (Semi-Analytical method for Barrier Option pricing) is easily applicable for pricing floating strike Asian barrier options with a continuous geometric average. Recently, this method has been applied in the Black–Scholes framework to European vanilla barrier options with constant and time-dependent parameters or barriers and to geometric Asian barrier options with a fixed strike price. The greater efficiency of SABO with respect to classical finite difference methods is clearly evident in numerical simulations. For the first time, a user-friendly MATLAB R © code is made available here.


Introduction
The pricing of continuously-monitored Asian options is a relevant task from both a mathematical and a financial point of view.
Asian options are quite common derivatives because they provide protection against strong price fluctuations in volatile markets and reduce the possibilities of price manipulations.The payoff of an Asian option depends on the average price of the underlying asset that is less volatile than the asset price itself.In general, Asian options are hence less valuable than their vanilla European counterparts because an option on a lower volatility asset is worth less.
On the other hand, it is more difficult to deal with Asian options than vanilla options because their price depends on the average value assumed by the underlying asset during the option's life, requiring some mathematical effort in order to describe the dynamics of the average under consideration.
In this paper, Asian options are equipped with a continuously-monitored geometric average [1].Asian options evaluated with the geometric mean, although not common among practitioners, give some information also about the evaluation of Asian options with the arithmetic mean [2].From a theoretical point of view, the method illustrated in the paper is extensible to arithmetic Asian options, as well, with slight modification, but from the numerical point of view, there are several problems that we plan to investigate in the near future.Defining the stochastic process A t := t 0 log(S τ )dτ (1) then the geometric average is defined as exp A t t .When A and S are written with subscripts (A t and S t ), they are intended as stochastic processes; otherwise, they are considered independent variables in the differential analysis context.The differential problem that describes the price evolution of this option is: V(S, A, T) assigned , S ∈ R + , A ∈ R.
Wanting to provide a further protection against excessive fluctuations of the strike price, it is possible to apply barriers in the option contract; for example, knock-out barriers make the option cease to exist if the underlying asset reaches a barrier during the life of the option.The model analyzed in this paper concerns an Asian option with an up-and-out barrier at S = B and a floating strike payoff, i.e., V(S, A, T) assigned , S ∈ (0, B), A ∈ R (5) asymptotic conditions of vanilla option , {(S, A) : with the final condition: The problem (2)-(3) of pricing a floating strike Asian option with a continuous geometric average and without a barrier has a closed-form solution in the domain (S, A, t) ∈ R + × R × [0, T) that can be formulated either as the payoff expected value (also known as the Feynman-Kac formula): with the transition probability density function G associated with the differential operator defined in Equation ( 4), which is known to be: or (see [3]) through the formula: for the call option, where N [•] is the normal cumulative distribution function, and eventually, the put-call parity: Instead, when applying barriers, no closed-formulas are available.In this context, SABO is a Semi-Analytical method conceived of for the pricing of Barrier Options, and its milestones are resumed in Section 4.1.It is quite a general method, applicable also to fixed strike payoffs [4,5], put options [6], time-dependent parameters [7] and double barriers [8].
SABO is compared here with two Finite Difference (FD) methods chosen among the wide class of numerical methods at our disposal [9].Equation ( 4) is proven to be hypo-elliptic [10][11][12], a property that guarantees a smooth solution and should benefit from approximations based on Taylor expansions.Anyway, SABO appears to be certainly more efficient looking at the results below.

Results
We have performed several simulations related to the pricing of a geometric call Asian option with a floating strike price and with an up-and-out barrier as modeled by the differential problem (4)- (8).Numerical results, some of which are displayed in the following, have been obtained by the MATLAB R codes implementing the algorithms of SABO, FD1 and FD2 described in Sections 4.1, 4.2.1 and 4.2.2, respectively.In Section 5 the SABO code is provided.
Example 1.In this example, we use the finance parameters displayed in Table 1.The floating strike call option with an up-and-out barrier at S = 150 =: B is evaluated at t = 0 and A = 0, truncating the A-domain at A min = 0 and A max = 5, in accordance with (20) and either (19) or (27).
The approximation by SABO is obtained setting the parameters described in Section 4.1 N t = N A = 20, and the option value is displayed in Figure 1 as a function of S.
The convergence at S = 100, 120, 140, 148 is shown in Table 2 refining the mesh: at each level, parameters N t and N A are doubled.
The related results are reported in Tables 3 and 4.
The approximation by FD2 is obtained setting the discretization parameters ∆t = ∆A = 0.01/2 k and either ∆S = 2 or ∆S = 1.The results are displayed in Tables 5 and 6.  1.
The convergence at S = 100, 120, 140 is shown in  Example 2. In this example, we use the finance parameters displayed in Table 7 and volatility equal to σ = 0.2, 0.3, 0.4.The floating strike call option with an up-and-out barrier at S = 115 =: B is evaluated at t = 0 and A = 0, truncating the A-domain at A min = 0 and A max = 5, in accordance with (20) and (19).
The approximation by SABO is obtained setting the parameters described in Section 4.1 N t = N A = 20 and the option value is displayed in Figure 2 as a function of S for the different values of σ (continuous lines), in comparison with the corresponding prices of Asian options without barriers (dotted lines).7 and various values of σ (continuous lines) compared with the corresponding prices of Asian options without barriers (dotted lines).

Discussion
Looking at the Example 1, the values of a call option with an up-and-out barrier obtained by SABO and displayed in Figure 1 show that the solution, as expected, assumes lower values than the analogous option without barriers whose closed formula is (12) or that can be computed through the evaluation of the payoff expected value (10).The same behavior is recovered by the two proposed FD methods (FD1 and FD2).
Talking about efficiency and convergence, we have to look at stabilization of digits in Tables 2-6 where   7 and various values of σ (continuous lines) compared with the corresponding prices of Asian options without barriers (dotted lines).

Discussion
Looking at Example 1, the values of a call option with an up-and-out barrier obtained by SABO and displayed in Figure 1 show that the solution, as expected, assumes lower values than the analogous option without barriers whose closed formula is (12) or that can be computed through the evaluation of the payoff expected value (10).The same behavior is recovered by the two proposed FD methods (FD1 and FD2).
Talking about efficiency and convergence, we have to look at the stabilization of the digits in Tables 2-6 where the option values at S = 100, 120, 140, 148 are written.SABO, the results of which are written in Table 2, appears to be faster than the FD methods: doubling parameters N t and N A , the CPU time for computation quadruples, but one more digit of accuracy is achieved.The convergence is slower near the barrier because there, the barrier option value is more different from the option value without barriers: the option value without the barrier can be quasi-exactly computed by the Feynman-Kac Formula (10), and therefore, the approximation error introduced by SABO solving the boundary integral Equation ( 15) related to the barrier case is more involved in representation Formula ( 14) as the asset nears the barrier (look at (26)).
Comparing Tables 3 and 4, with an analogous computational time, SABO appears much more accurate and therefore efficient than FD1.Furthermore, note that FD1 is still sensitive to the mesh refinement in the S-domain: to halve ∆S means a significant variation in values of V together with a big increase of the computational costs.
Analyzing Tables 5 and 6, we observe that FD2 has a superior accuracy compared to FD1 due to its higher order of consistency: approximations of derivatives in the t and A variables are both of second order.Anyway, FD2 is less efficient than SABO, and the coarseness of the S-grid still significantly affects its results.Refinements in the S-grid would result in much longer computational times, no longer comparable with those of SABO.
Looking at Example 2, SABO maintains its robustness varying the volatility values.The solutions displayed in Figure 2 show the property of smoothness proven in [11].The increase in volatility causes the expected increase in the value of vanilla options, but on the contrary, it implies a diminishing of barrier option values near the barrier.

SABO
SABO is a Semi-Analytical method for the pricing of Barrier Options.The foundation on which it is based is the integral representation of the solution of the modeling differential problem based on the knowledge of the transition probability density function.
For a problem like that in ( 4)-( 7), it is proven in [4] that the integral representation formula is: where the transition probability density function G(S, A, t; B, A, t) is defined in (11) and associated with the differential operator defined in Equation ( 4).Note that in (14), both V(S, A, t) and ∂V ∂ S (B, A, t) are unknown, but, at S = B, the Boundary Condition (6) can be applied giving rise to the Volterra integral equation of the first kind: in the unknown ∂V ∂ S (B, A, t).The unknown is approximated by: having defined piecewise constant basis functions: on a uniform time grid: and piecewise constant basis functions: on a uniform A-grid over the truncated A-domain [A min , A max ]: A careful choice of A min and A max has to be performed in such a way that the double integral: is rightly approximated by: By the analysis developed in [5], this results in seeking for A max such that: and A min such that: with a suitable tolerance.Otherwise, it is possible to consider the whole unbounded A-domain ≡ R, but with two infinite basis functions, the first and the last, as investigated in [4].
In order to numerically solve the Volterra integral Equation ( 15), we convert it into a discrete linear system of equations by means of collocation BEM.Hence we collocate (15) at points: finally obtaining, for i = 1, . . ., N A , j = 1, . . ., N t : The unknown vector is: and the matrix entries are equal to the case of fixed strike Asian options deeply investigated also from a computational point of view in [4,5]: for i, h = 1, . . ., N A , j, k = 1, . . ., N t , define ξ = i − h, ξ = −N A + 1, . . ., N A − 1 and = k − j, = 0, . . ., N t − 1, thus obtaining: As we are considering constant coefficients in (4), the fundamental solution (11) depends on t, t, A, A only through the differences t − t and A − A implying that the matrix has a block-Toeplitz structure both in time and A-space: so it is possible to adopt suitable strategies to save computational costs (as done in [13], even if this feature is not implemented in the code included here) and memory requirements.
A change in the payoff function is caught by the evaluation of the rhs entries: in the case of a floating strike call payoff, for i = 1, . . ., N A , j = 1, . . ., N t : with Erfc the complementary of error function Erf.
The approximation of ∂V ∂ S (B, A, t) through the resolution of ( 21) is an intermediate step for the final evaluation of the option price by ( 14): (26)

Finite Difference Methods
SABO is compared here (as in [5]) with two Finite Difference (FD) methods chosen among the wide class of FD methods available and deeply analyzed in [14].This is because Equation ( 4) is proven to be hypoelliptic [10][11][12], a property that guarantees a smooth solution and that should benefit from approximations based on Taylor expansions.
The existence of the solution of Problem ( 4)-( 5) and its Feynman-Kac representation (10) are proven involving stochastic arguments without the need for exact boundary conditions.Boundary conditions at S = 0 and for A → ±∞ have to be empirically deduced analogously to what was done in [5,15].At S = B, Condition (6) holds.
In order to apply FD methods, the first step is to shrink to a bounded domain [0, B] × [A min , A max ] enclosing the option evaluation point (S * , A * ) and then apply boundary conditions at the borders.

•
We have chosen A min less than or equal to the minimum between the value suggested in (20) and A * .There is not a consistent condition valid at A = −∞ in finance, and as a consequence, it is not easy to conceive of a proper condition at points (S, A min , t).The use of the upwind method makes this condition unnecessary if log(S) > 0 (as in the herein proposed examples).In fact, ∂V ∂A is approximated by a forward difference so that the boundary condition at A min becomes useless; otherwise, if log(S) < 0, a backward difference can be considered together with a condition on the derivative at A = A min , as for example The upper bound is set at: the maximum value assumed by the stochastic process A t defined in (1) if S t = B during the whole time interval [0, T].
The average exp(A/T) is a non-decreasing function of time; therefore, if exp(A/t) − B > 0 at time t, then exp(A/T) − B > 0 at maturity, and so, the option is worth nothing as shown by the Payoff Function (8).Hence, if log(S) > 0, ∂V ∂A is approximated by a forward difference and the boundary condition at A max is set: • At S = 0, Equation ( 4) is degenerate, but if the stochastic process S t = 0 at any time, then the average asset price exp(A t /t) = constant and the option value V can be considered as independent of stochastic variables S and A, deducing from Equation (4) that: The ordinary differential equation: The second step, after the definition of the computational domain, is the definition of the grid.Using a (A, t)-grid as defined by ( 17) and ( 18) in [A min , A max ] × [0, T], let us further introduce the S-grid in [0, B]: and define the approximated option value:

•
The values V N t ih can be found from the final condition: • in the herein proposed examples, log(S i ) > 0, so everywhere, we apply a forward difference and, at the boundary A = A max , Condition (28): • at the boundary S = 0, we apply Condition (29): • at the boundary S = B, we apply the Boundary Condition (6), hence:

Finite Difference Method FD1
The first scheme FD1 is the very classical FD scheme where the derivatives of V in (4) are approximated by truncations of Taylor expansions: • first order backward difference for the time derivative approximation: • second order central difference for the S derivative approximation: • second order central difference for the S second order derivative approximation: • if log(S i ) > 0, first order forward difference for the A derivative approximation: We can now write down the discrete approximation of (4) at each point (S i , A h , t k ): and rearrange the scheme in a compact form: where: It is easy to see that the scheme, backward computing the new value V k−1 i,h , is explicit in time.We want to remark about two more things: first, the weights depend only on i, therefore on the S variable, and second, the values V k i,0 having A min as the coordinate give no contribution to the scheme.

Finite Difference Method FD2
The method FD2 is suggested in [16].The PDE ( 4) is collocated at the points: Then, we use the following approximations, based on suitable Taylor expansions and standard finite difference approximations: After substituting the above approximations into the PDE and discarding the error terms, we get the following equations for the approximate values of the option prices: where, using the notation λ = ∆t ∆S 2 : The procedure for solving the option pricing equation is as follows: 1.
Fill the values For each k = N t − 1 : −1 : 0: This is a time-explicit difference scheme, as well.If we need the option price only at t = 0, then it is not necessary to store the full matrix V of approximate option prices; in fact, we need only two levels, say V old corresponding to t = t k+1 and V new corresponding to the current time level t = t k .At the beginning, V old is computed using the final condition, and at the end of each time step, the values of V new are copied to V old .Anyway, it requires the resolution of a linear system: the three-diagonal matrix M = diag(a, b, c) can be assembled and factorized at the beginning, outside the cycles, depending only on the S-grid.

MATLAB R Code Implementing SABO
All the above-provided numerical results were obtained by codes developed with MATLAB R Release 2007b running on a laptop computer (CPU Intel i5, 4 Gb RAM).The code implementing the SABO algorithm is given below.

Figure 1 .
Figure1.SABO approximation of a call up-and-out Geometric Asian option with floating strike and data in Table1.

Figure 2 .
Figure 2. SABO approximation of a call up-and-out Geometric Asian option with floating strike, data in Table7and various values of σ (continuous lines) compared with the corresponding prices of Asian options without barriers (dotted lines).

¼
the option values at S = 100, 120, 140, 148 are written.½ SABO, whose results are written in Table2appears to be faster than FD methods: doubling of ¾ parameters N t and N A the CPU time of computation quadruples but one more digit of accuracy is ¿ achieved.The convergence is slower near the barrier because there the barrier option value is more different from the option value without barriers: the option value without barrier can be quasi-exactly computed by the Feynman-Kac formula(10) and therefore the approximation error introduced by

Figure 2 .
Figure 2. SABO approximation of a call up-and-out geometric Asian option with a floating strike, the data in Table7and various values of σ (continuous lines) compared with the corresponding prices of Asian options without barriers (dotted lines).
condition at A = A max to define V k i,N A , for i = 0, • • • , N S .(b)For each h = N A − 1 : −1 : 0, solve for the three-diagonal system in the unknowns the values V k i,h , for i = 1, • • • , N S − 1, using the boundary condition at S = 0 and S = B.

Table 1 .
Floating strike up-and-out call option data.

Table 7 .
Floating strike up-and-out call option data.