Multi-Asset Barrier Options Pricing by Collocation BEM (with Matlab ® Code)

: In this paper, we extend the SABO technique (Semi-Analytical method for Barrier Options), based on collocation Boundary Element Method (BEM), to the pricing of Barrier Options with payoff dependent on more than one asset. The efﬁciency and accuracy already revealed in the case of a single asset is conﬁrmed by the presented numerical


Introduction to the Differential Model Problem
Partial Differential Equations (PDEs) of Mathematical Physics model a huge variety of real-life problems, from science to engineering. Since the work [1] of F. Black and M. Scholes, equations that model physical phenomena have been reconsidered to interpret financial phenomena. PDEs in space-time variables, modeling the price of the most evolved financial products, need efficient techniques to be numerically solved.
This work investigates the extensions of the so called SABO technique (Semi-Analytical method for Barrier Options), based on collocation Boundary Element Method (BEM) and applied so far for the numerical pricing of barrier options in a one dimensional asset framework [2].
The consideration of two dimensional partial differential problems can be suggested for several reasons: the desire to complicate the model to get closer to reality (such as, for example, introducing the dependence of option value on a stochastic volatility variable [3]) or the evaluation of options that depend not only on the asset value. In this last direction we have already approached Asian options whose payoff depends on the average of asset values [4,5] giving rise to a degenerate PDE based on two independent variables: the asset value and the average.
In this article, the extension is devoted to options whose payoff depends on more than one asset. The consequent differential model, described in the following, is set by a parabolic equation that, with suitable transformations, can be traced back to the heat equation [6].
From the computational point of view, the pricing of multi-asset options is recognized in the literature as a quite difficult issue and difficulties increase with the application of barriers [7]. The problem can be tackled starting from stochastic differential equations by Monte Carlo methods [8,9] or considering the problem in its partial differential formulation [10,11] with some limitations on the number of assets involved.
In principle our strategy can be applied to an undefined number of underlying assets, so theoretically the extension is straightforward, but the choice of approximation technique is crucial to face the curse of dimensionality: the application of collocation method ensures more efficiency than classical domain methods at low dimensions. This will be detailed in Section 4.
We assume the Black-Scholes-Merton scenario for the evaluation of an option V(S 1 , . . . , S n , t) based on n assets S 1 , . . . , S n during the time interval [0, T], with no dividends: the behavior of the underlying assets is described by a geometric Brownian motion and −1 ≤ ρ ij ≤ 1 is the correlation between the two assets i and j. The application of n-dimensional Itô Lemma and no-arbitrage arguments leads to the related backward parabolic differential model problem ∂V ∂t where r is the risk-free interest rate, σ i is the volatility of underlying asset S i and V T is the option contract payoff at the expiry T that may depend on a strike price K and assume several expressions, for example, looking at call options in [12]: for different kinds of Rainbow options, max(S 1 , . . . , S n ) n-color better-of option min(S 1 , The discussion in the following sections can be developed in the general dimension n [13]; however, for the sake of clarity and to simplify numerics, we will detail it considering two underlying assets S 1 , S 2 only. Hence, the convenience and reliability of SABO will be highlighted with proofs and numerical examples for a two assets framework, and, in continuity with the previous article in this journal [4], we have inserted in Appendix A a ready to use Matlab ® code.

Integral Representation Formula of the Solution for Options without Barriers
By the Green's theorem, we prove the integral representation formula for the solution of the differential model problem (1) describing the value of an option V(S 1 , S 2 , t) based on two assets S 1 , S 2 during the time interval [0, T].
We recall the notion of joint transition probability density function G(S 1 , S 2 , t;S 1 ,S 2 ,t) (also known as Green's or fundamental solution) from the classical theory of PDEs [13]: denoting the space differential operator by w.r.t. the first tern of variables and of the adjoint model problem (4) and (5) w.r.t. the second tern of variables: G(S 1 , S 2 , t;S 1 ,S 2 ,t) = δ(S 1 ,S 1 )δ(S 2 ,S 2 ) (S 1 ,S 2 ) ∈ R + × R + ; moreover it is such that having set (·, ·) the standard L 2 scalar product. The adjoint operator L * is defined as From [14] the fundamental solution is known to be with the elements of the array α and Σ = 1 ρ 12 ρ 12 1 .
the correlation matrix.
Proof. From (8), by the changes of variables ξ i := log(S i ) and exploiting the property Remark 1. Note that I 1 and I 2 turn out to coincide with the coefficients of the changes of variables G = e r(t−t) G andG = e (3r−σ 1 σ 2 ρ 12 −σ 2 2 )(t−t) G allowing to obtain (3) and (4), respectively, without the last term.
With the knowledge of the fundamental solution, we will derive, by mathematical analysis tools, the following integral representation formula: giving the option value as expected value of the payoff function, holds.
and then, thanks to property (6), these steps lead to the Feynman-Kac formula.

Integral Representation Formula of the Solution for Barrier Options
If we consider barrier options, sometimes analytical solutions are known (some examples are collected in [15]) but sometimes not. For example, in [12], the author considers the case of a European-style two-assets-basket double-barrier call option with payoff and two knock-out barrier conditions with down-and-out barrier B 1 and up-and-out barrier B 2 . Then, the resulting domain Ω is partitioned as shown in Figure 1 in order to approximate the option value by a Finite Element Method. However, if we are interested in the computation of the option value only at a desired point (S 1 , S 2 ) of the domain, it may certainly be convenient to apply the strategy already tested for various kind of options (based on one asset only) under several dynamics [2][3][4][5]16,17] and that we called SABO.
The method is based on a new analytical representation formula whose proof relies on Green's Theorem and retraces the proof of Feynman-Kac formula. (1), (13) and (14), ∀(S 1 , S 2 , t) ∈ Ω ⊆ R 2 × [0, T), can be expressed by the following integral representation formula:

Proposition 3. The solution of problem
with the functional ϕ defined in (17) below.

Proof. Multiply (1) by G and integrate in time and in the domain
then, integrating by parts in time and applying Green's Theorem and (6), one obtains having defined ∂Ω as the boundary of the domain, n as its unit normal vector outwardly directed and p = (p 1 , p 2 ) as Thus we can conclude that Remark 2. The representation formula, and SABO in general, can be naturally adapted to barrier options configurations different from the proposed example. Note also that it is not really important to know exactly the expression of the functional ϕ because, as it depends on the unknown solution V, it is unknown in its turn and therefore it will be numerically determined.

Remark 3.
The representation formula can be analytically derived in time or in space to obtain Greeks functions that can be straightforwardly evaluated by SABO without the need of evaluating option values (look at [2]).
In the example (12)- (14), due to the configuration of the boundary ∂Ω = Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 (see Figure 1), the integral of ϕ reduces to the first component of the normal vector is null and every term of p 2 has the factor S 2 , trivial on Γ 2 ; on Γ 4 := {(0, S 2 ) : B 1 ≤ S 2 ≤ B 2 } the second component of the normal vector is null and every term of p 1 has the factor S 1 , trivial on Γ 4 . Moreover, considering double knock-out barriers on Γ 1 :

Boundary Integral Equation
The representation formula (18) of V holds in the whole domain and can be used once the unknown density φ is recovered. To this aim, the idea is to consider the application of the representation formula at the barriers, where ∀(S 1 , S 2 , t) ∈ Γ 1 ∪ Γ 3 × [0, T) the option value is trivial, i.e., V(S 1 , S 2 , t) = 0 and therefore and with the general expression of fundamental solution G written in [14].

Remark 5.
If the asset domain is either a 2 or 3 or 4 dimensional domain, to solve (21) means to reduce the dimensionality of the problem by one. In this case, the collocation method can be implemented, as proposed in Section 5 (introducing suitable mesh triangulations), since finite element discretization can be easily applied up to three dimensions with presumably greater accuracy w.r.t. Monte Carlo simulations. However, of course, the collocation method suffers the curse of dimensionality. Therefore, first of all, SABO in this framework, is suggested as an alternative to deterministic methods (finite difference methods and finite element methods) for its higher efficiency.
With the increasing of the dimension n > 4, (21) is still valid and very general, not requiring special conditions as for example in [18], so the strategy of solving (21) can be anyway taken into account. An idea could be to involve Monte Carlo method, not as a method for path simulation [9], but considered as a quadrature method [19] directly applied to the integral terms in (21). Then unknown φ might be represented by linear combination of radial basis functions around some nodes spread out in the n − 1-dimensional hyperplane Γ. This will be the subject of future investigation.

Numerical Approximation by Collocation BEM
We introduceφ for the numerical approximation of φ by: • piecewise constant functions in timê  The related unknown array α = α (1) α (1) α (2) . . . α (N ∆t ) α ( ) = α 1 , · · · , α M ∆S , = 1, . . . , N ∆t can be computed by solving a linear system and this is due to the form of the fundamental solution (8): the starting PDE Equation (1) has coefficients independent of time implying that the fundamental solution depends on time only through the differencet − t and, by consequence, each block depends only on t n −t n . From the computational point of view, this means notable computational and memory savings because only the last column of blocks needs to be computed and the linear system (23) can be solved by block-backward substitution where only the non singular block M (1) needs to be inverted. We can get the expression of a general entry of the matrix block M ( ) , = 1, . . . , N ∆t , introducing, in the lhs of (19), a space-time piecewise constant approximation (22) and collocation points (24) and, collocating the rhs of (19) in the same way, we can write Once the array α of coefficients in (22) are computed by solving the linear system (23), the unknown functional φ in (18) can be replaced by its approximationφ, in order to get the approximation of the option price, at the desired time instant and assets values, without the need of introducing grids or triangulation over the domain Ω as in Figure 1.
Note that the evaluation of system entries can be simply performed by Matlab ® "integral2" function as it presents only little troubles: the degeneration of the fundamental solution in time towards the Dirac delta distribution when t →t (look at (5)) and possible non-smoothness of the payoff function like, for example, in the case of Basket options.

Numerical Results and Discussion
• Consider the problem (1), (12), (13) and (14) of evaluating a double knock-out basket call option based on two assets with data suggested in [12] and listed in Table 1. The domain is as depicted in Figure 1. Numerical results have been obtained by the code inserted in the Appendix A. The approximated solution, represented in Figure 2, is evaluated over a rectangular grid of points, vertices of 16 2 simplices that subdivide the square [0, 2] 2 . It is obtained by SABO with ∆S = 0.125 √ 2 and ∆t = 0.1. The behavior of the solution is in good agreement with that found in [12] ("Courtesy of A. Kvetnaia"), there approximated by Finite Element Method (FEM) after having transformed Equation (1) in a divergence-free form, not necessary or useful to apply SABO.  Table 1.
The greater accuracy, implying greater velocity, of SABO w.r.t. Finite Differences and Monte Carlo Methods has been already highlighted in our papers about barrier options, see for example [3,20]. Here, in Table 2, we compare the CPU time of computation (Laptop computer: CPU Intel i5, 4 Gb RAM.) of SABO w.r.t. FEM as applied in [12], i.e., linear finite elements in space coupled with Crank-Nicolson scheme in time: the chosen reference value is V(0.5, 1, 0) = 0.306264, obtained by stressing the discretization parameters in our SABO code (M ∆S = N ∆t = 30, tolerance on integrals computation equal to 10 −12 ) and both methods converge towards it. The starting space mesh for FEM, that will be denoted by R 0 from now on, is made by 48 uniform triangles as depicted in Figure 3, then any successive refinement, R i i = 1, 2, 3 , is obtained by dividing each triangle in four uniform triangles (by connecting the midpoints of the edges). Despite not having implemented a parallel code for the computation of independent blocks in (26), SABO allows to achieve the same accuracy of FEM with a computation time of at least one lower order of magnitude.   Table 1, mesh R 0 and ∆t = 0.02.
In order to quantitatively check the reliability of numerical results, we show in Figure 4 that the approximated solution evaluated inside the domain but moving towards boundary Γ 4 (or equivalently Γ 2 ) converges to the exact solution of a double knock-out call option with the same financial parameters but based on one asset, whose explicit expression C(S 1 , t) can be found in [21] and is here reported for reader's convenience: with N the standard normal cumulative distribution function. These "boundary conditions" on Γ 2 and Γ 4 are necessary to apply Finite Element Method; on the contrary, they are not set during SABO implementation, but they are naturally matched by the approximated solution.
In order to show the stability of results we have changed also some financial parameters starting from Table 1: in Figure 5, the value of the strike price (on the top) and the volatility of one asset (on the bottom) and, in Figure 6, the correlation (poorly correlated assets with ρ = 0.1 on the top, highly correlated assets with ρ = 0.9 on the bottom).    • Consider a put down and out option based on two assets whose value is the solution of (1) provided with the payoff function meaning that there is a lower out barrier on the first asset only. The boundary and the domain are shown in Figure 7 top. In Figure 7 bottom, there is the representation of the approximated solution at time t = 0, obtained with M ∆S = N ∆t = 10. In Table 3, we can observe the stabilization of digits in the numerical approximation of V(S 1 , S 2 , 0) at (S 1 , S 2 ) = (2, 1) doubling the number of discretization nodes M ∆S and N ∆t in (22). Table 3. Stabilization of value V(S 1 , S 2 , 0) at (S 1 , S 2 ) = (2, 1) with the refinement of discretization parameters ∆ S and ∆ t .

Conclusions
In this paper, we have extended the SABO method, based on collocation BEM, for pricing barrier options depending on two assets. The good performance of this technique, largely employed in the past for options on a single asset, has been shown by means of numerical results. The implemented Matlab ® code is enclosed in the Appendix A, in such a way that it can be directly used by any interested reader.
Author Contributions: All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.
Funding: This work has been partially supported by INdAM-GNCS Research Projects.