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

Abstract

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 efficiency and accuracy already revealed in the case of a single asset is confirmed by the presented numerical results.

1. 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,\dots ,S_n,t)$ based on n assets $S_1,\dots ,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\le {\rho }_{ij}\le 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

$$\begin{array}{c}\hfill \frac{\partial V}{\partial t}+\frac{1}{2}\sum _{i,j=1}^n{\rho }_{ij}{\sigma }_i{\sigma }_j{S}_i{S}_j\frac{{\partial }^{2}V}{\partial S_i\partial S_j}+r\sum _{i=1}^n{S}_i\frac{\partial V}{\partial S_i}-rV=0,(S_1,\dots ,S_n,t)\in {\left[{\mathbb{R}}^{+}\right]}^n\times \left[0,T\right)\\ \hfill V(S_1,\dots ,S_n,T)=V_T(S_1,\dots ,S_n),(S_1,\dots ,S_n)\in {\left[{\mathbb{R}}^{+}\right]}^n\end{array}$$

(1)

where r is the risk-free interest rate, ${\sigma }_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]:

• $max\left({\sum }_{i=1}^n{c}_i{S}_i-K,0\right)$, for Basket options with weights $c_i$;
• for different kinds of Rainbow options,
  $max(S_1,\dots ,S_n)$n-color better-of option
  $min(S_1,S_2)$two-color worse-of option
  $max(S_2-S_1,0)$outperformance option
  $max(min(S_1-K,\dots ,S_n-K),0)$min option 
• $max(S_1-S_2-K,0)$, for spread option

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.

2. 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;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})$ (also known as Green’s or fundamental solution) from the classical theory of PDEs [13]: denoting the space differential operator by

$$\mathcal{L}\left[V\right]:=\frac{{\sigma }_1^2{S}_1^2}{2}\frac{{\partial }^{2}V}{\partial S_1^2}+r{S}_1\frac{\partial V}{\partial S_1}+\frac{{\sigma }_2^2{S}_2^2}{2}\frac{{\partial }^{2}V}{\partial S_2^2}+r{S}_2\frac{\partial V}{\partial S_2}+\rho {\sigma }_1{\sigma }_2{S}_1{S}_2\frac{{\partial }^{2}V}{\partial S_1\partial S_2}-rV,$$

(2)

G is solution of (1), i.e.,

$$\frac{\partial G}{\partial t}(S_1,S_2,t)+\mathcal{L}\left[G\right](S_1,S_2,t)=0$$

(3)

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:

$$\begin{array}{c}\hfill -\frac{\partial G}{\partial \tilde{t}}({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})+{\mathcal{L}}^{\ast }\left[G\right]({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=0({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\times \left(t,+\infty \right)\end{array}$$

(4)

$$\begin{array}{c}\hfill G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=\delta (S_1,{\tilde{S}}_1)\delta (S_2,{\tilde{S}}_2)({\tilde{S}}_1,{\tilde{S}}_2)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+};\end{array}$$

(5)

moreover it is such that

$$(\mathcal{L}\left[V\right],G)=(V,{\mathcal{L}}^{\ast }\left[G\right])$$

(6)

having set $(·,·)$ the standard $L^2$ scalar product. The adjoint operator ${\mathcal{L}}^{\ast }$ is defined as

$$\begin{array}{ccc}\hfill {\mathcal{L}}^{\ast }\left[G\right]& :=& {\displaystyle \frac{{\sigma }_1^2{\tilde{S}}_1^2}{2}\frac{{\partial }^{2}G}{\partial {\tilde{S}}_1^2}+\frac{\partial G}{\partial {\tilde{S}}_1}(-r{\tilde{S}}_1+2{\sigma }_1^2{\tilde{S}}_1+\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1)}\hfill \\ \hfill & +& {\displaystyle \frac{{\sigma }_2^2{\tilde{S}}_2^2}{2}\frac{{\partial }^{2}G}{\partial {\tilde{S}}_2^2}+\frac{\partial G}{\partial {\tilde{S}}_2}(-r{\tilde{S}}_2+2{\sigma }_2^2{\tilde{S}}_2+\rho {\sigma }_1{\sigma }_2{\tilde{S}}_2)}\hfill \\ \hfill & +& {\displaystyle \rho {\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\frac{{\partial }^{2}G}{\partial {\tilde{S}}_1\partial {\tilde{S}}_2}+G(\rho {\sigma }_1{\sigma }_2-3r+{\sigma }_1^2+{\sigma }_2^2).}\hfill \end{array}$$

(7)

From [14] the fundamental solution is known to be

$$G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=\frac{e^{-r(\tilde{t}-t)}}{2\pi (\tilde{t}-t)}\frac{exp\left(-\frac{{\mathbf{\alpha }}^{\prime }{\mathbf{\Sigma }}^{-1}\mathbf{\alpha }}{2}\right)}{{\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\sqrt{det\mathbf{\Sigma }}}$$

(8)

with

$${\alpha }_i=\frac{log\frac{S_i}{{\tilde{S}}_i}+\left(r-\frac{{\sigma }_i^2}{2}\right)(\tilde{t}-t)}{{\sigma }_i\sqrt{\tilde{t}-t}},i=1,2$$

the elements of the array $\mathbf{\alpha }$ and

$$\mathbf{\Sigma }=\left(\begin{array}{cc}1& {\rho }_{12}\\ {\rho }_{12}& 1\end{array}\right).$$

the correlation matrix.

Proposition 1.

The fundamental solution (8) of (3)–(5) verifies, $\forall \tilde{t},t\in [0,T]$, the following identities

$$\begin{array}{ccc}\hfill {\displaystyle {\mathcal{I}}_1:={\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d{\tilde{S}}_{1}d{\tilde{S}}_2}& =& {\displaystyle e^{-r(\tilde{t}-t)}}\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\displaystyle {\mathcal{I}}_2:={\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d{S}_{1}d{S}_2}& =& {\displaystyle e^{-(\tilde{t}-t)(3r-{\sigma }_1{\sigma }_2{\rho }_{12}-{\sigma }_1^2-{\sigma }_2^2)}}\hfill \end{array}$$

(10)

Proof. 

From (8),

$$G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=\frac{e^{-r(\tilde{t}-t)}}{2\pi (\tilde{t}-t)}\frac{exp\left(-\frac{{\alpha }_1^2+{\alpha }_2^2-2{\alpha }_1{\alpha }_2{\rho }_{12}}{2(1-{\rho }_{12}^2)}\right)}{{\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\sqrt{1-{\rho }_{12}^2}}$$

hence

$${\mathcal{I}}_1=\frac{e^{-r(\tilde{t}-t)}}{2\pi (\tilde{t}-t){\sigma }_1{\sigma }_2\sqrt{1-{\rho }_{12}^2}}{\int }_{{\mathbb{R}}^{+}}\frac{exp\left(-\frac{{\alpha }_1^2}{2(1-{\rho }_{12}^2)}\right)}{{\tilde{S}}_1}d{\tilde{S}}_1{\int }_{{\mathbb{R}}^{+}}\frac{exp\left(-\frac{{\alpha }_2^2-2{\alpha }_1{\alpha }_2{\rho }_{12}}{2(1-{\rho }_{12}^2)}\right)}{{\tilde{S}}_2}d{\tilde{S}}_2$$

by the changes of variables ${\xi }_i:=log\left({\tilde{S}}_i\right)$ and exploiting the property ${\int }_0^{+\infty }{e}^{-{\xi }^2}d\xi =\sqrt{\pi }/2$, we obtain

$$\begin{array}{ccc}& =& \frac{e^{-r(\tilde{t}-t)}}{2\pi (\tilde{t}-t){\sigma }_1{\sigma }_2\sqrt{1-{\rho }_{12}^2}}{\int }_{{\mathbb{R}}^{+}}{e}^{-\frac{{\alpha }_1^2}{2}}\sqrt{2\pi {\sigma }_2^2(\tilde{t}-t)(1-{\rho }_{12}^2)}d{\xi }_1\hfill \\ & =& \frac{e^{-r(\tilde{t}-t)}}{\sqrt{2\pi (\tilde{t}-t)}{\sigma }_1}{\int }_{{\mathbb{R}}^{+}}{e}^{-\frac{{\alpha }_1^2}{2}}d{\xi }_1=e^{-r(\tilde{t}-t)}.\hfill \end{array}$$

Analogously

$$\begin{array}{ccc}\hfill {\mathcal{I}}_2& =& \frac{e^{-r(\tilde{t}-t)}}{2\pi (\tilde{t}-t){\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\sqrt{1-{\rho }_{12}^2}}{\int }_{{\mathbb{R}}^{+}}{e}^{-\frac{{\alpha }_1^2}{2(1-{\rho }_{12}^2)}}d{S}_1{\int }_{{\mathbb{R}}^{+}}{e}^{-\frac{{\alpha }_2^2-2{\alpha }_1{\alpha }_2{\rho }_{12}}{2(1-{\rho }_{12}^2)}}d{S}_2\hfill \\ & =& \frac{e^{(\tilde{t}-t)(-2r+{\sigma }_2^2-\frac{{\rho }_{12}^2{\sigma }_2^2}{2})}}{\sqrt{2\pi (\tilde{t}-t)}{\sigma }_1{\tilde{S}}_1}{\int }_{{\mathbb{R}}^{+}}{e}^{-\frac{{\alpha }_1^2}{2}+{\alpha }_1{\rho }_{12}(\tilde{t}-t){\sigma }_2}d{S}_1=e^{(\tilde{t}-t)(-3r+{\sigma }_2^2+{\sigma }_1^2+{\sigma }_1{\sigma }_2{\rho }_{12})}\hfill \end{array}$$

   □

Remark 1.

Note that ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$ turn out to coincide with the coefficients of the changes of variables $\tilde{G}=e^{r(\tilde{t}-t)}G$ and $\tilde{G}=e^{(3r-{\sigma }_1{\sigma }_2{\rho }_{12}-{\sigma }_1^2-{\sigma }_2^2)(\tilde{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:

Proposition 2.

The Feynman–Kac formula

$$V(S_1,S_2,t)={\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2$$

(11)

giving the option value as expected value of the payoff function, holds.

Proof. 

Multiply (1) by G and integrate in time-space domain

$$\begin{array}{ccc}\hfill 0& =& {\int }_t^T{\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}\frac{\partial V}{\partial \tilde{t}}({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}\mathcal{L}\left[V\right]({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

then, integrating by parts in time and in space and applying (6), one obtains

$$\begin{array}{ccc}& =& {\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,t)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,t)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_t^T{\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\frac{\partial G}{\partial \tilde{t}}(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t}){\mathcal{L}}^{\ast }\left[G\right](S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

$$\begin{array}{ccc}& =& {\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,t)\delta (S_1,{\tilde{S}}_1)\delta (S_2,{\tilde{S}}_2)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_t^T{\int }_{{\mathbb{R}}^{+}\times {\mathbb{R}}^{+}}V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\left(\frac{\partial G}{\partial \tilde{t}}-{\mathcal{L}}^{\ast }\left[G\right]\right)(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

and then, thanks to property (6), these steps lead to the Feynman-Kac formula.    □

3. 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

$$V(S_1,S_2,T)=max(S_1+S_2-K,0)$$

(12)

and two knock-out barrier conditions

$$\begin{array}{ccc}\hfill V(S_1,S_2,T)& =& 0S_1+S_2\le B_1\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill V(S_1,S_2,T)& =& 0S_1+S_2\ge B_2\hfill \end{array}$$

(14)

with down-and-out barrier $B_1$ and up-and-out barrier $B_2$. Then, the resulting domain $\Omega$ 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.

Proposition 3.

The solution of problem (1), (13) and (14), $\forall (S_1,S_2,t)\in \Omega \subseteq {\mathbb{R}}^2\times [0,T)$, can be expressed by the following integral representation formula:

$$\begin{array}{ccc}\hfill V(S_1,S_2,t)& =& {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{\partial \Omega }\phi ({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

(15)

with the functional φ defined in (17) below.

Proof. 

Multiply (1) by G and integrate in time and in the domain $\Omega$

$$\begin{array}{ccc}\hfill 0& =& {\int }_t^T{\int }_{\Omega }\frac{\partial V}{\partial \tilde{t}}({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{\Omega }\mathcal{L}\left[V\right]({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

then, integrating by parts in time and applying Green’s Theorem and (6), one obtains

$$\begin{array}{ccc}& =& {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,t)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,t)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_t^T{\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\frac{\partial G}{\partial \tilde{t}}(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t}){\mathcal{L}}^{\ast }\left[G\right](S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{\partial \Omega }\mathbf{p}(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})·\mathbf{n}d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

having defined $\partial \Omega$ as the boundary of the domain, n as its unit normal vector outwardly directed and $\mathbf{p}=(p_1,p_2)$ as

$$\begin{array}{ccc}\hfill p_1& =& {\displaystyle \frac{1}{2}{\sigma }_1^2{\tilde{S}}_1^2\left(G\frac{\partial V}{\partial {\tilde{S}}_1}-V\frac{\partial G}{\partial {\tilde{S}}_1}\right)+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\left(G\frac{\partial V}{\partial {\tilde{S}}_2}-V\frac{\partial G}{\partial {\tilde{S}}_2}\right)}\hfill \\ \hfill & -& {\displaystyle VG\left({\sigma }_1^2{\tilde{S}}_1+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1-r{\tilde{S}}_1\right)}\hfill \\ \hfill p_2& =& {\displaystyle \frac{1}{2}{\sigma }_2^2{\tilde{S}}_2^2\left(G\frac{\partial V}{\partial {\tilde{S}}_2}-V\frac{\partial G}{\partial {\tilde{S}}_2}\right)+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\left(G\frac{\partial V}{\partial {\tilde{S}}_1}-V\frac{\partial G}{\partial {\tilde{S}}_1}\right)}\hfill \\ \hfill & -& {\displaystyle VG\left({\sigma }_2^2{\tilde{S}}_2+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_2-r{\tilde{S}}_2\right)}\hfill \end{array}.$$

(16)

Thus we can conclude that

$$\begin{array}{ccc}\hfill 0& =& {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,t)\delta (S_1,{\tilde{S}}_1)\delta (S_2,{\tilde{S}}_2)d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & -& {\int }_t^T{\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\left(\frac{\partial G}{\partial \tilde{t}}-{\mathcal{L}}^{\ast }\left[G\right]\right)(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \\ & +& {\int }_t^T{\int }_{\partial \Omega }\phi (S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2\hfill \end{array}$$

with

$$\phi (S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=\mathbf{p}(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})·\mathbf{n}$$

(17)

   □

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 $\partial \Omega ={\Gamma }_1\cup {\Gamma }_2\cup {\Gamma }_3\cup {\Gamma }_4$ (see Figure 1), the integral of $\phi$ reduces to

$${\int }_t^T{\int }_{\partial \Omega }\phi (S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2={\int }_t^T{\int }_{{\Gamma }_1\cup {\Gamma }_3}\phi (S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2$$

because on ${\Gamma }_2:=\left\{(S_1,0):B_1\le S_1\le B_2\right\}$ the first component of the normal vector is null and every term of $p_2$ has the factor $S_2$, trivial on ${\Gamma }_2$; on ${\Gamma }_4:=\left\{(0,S_2):B_1\le S_2\le B_2\right\}$ the second component of the normal vector is null and every term of $p_1$ has the factor $S_1$, trivial on ${\Gamma }_4$.

Moreover, considering double knock-out barriers on ${\Gamma }_1:=\left\{(S_1,S_2):S_1+S_2=B_2\right\}$ and on ${\Gamma }_3:=\left\{(S_1,S_2):S_1+S_2=B_1\right\}$, the option value V is equal to 0, justifying the representation formula: $\forall (S_1,S_2,t)\in \Omega \times [0,T)$

$$\begin{array}{ccc}\hfill V(S_1,S_2,t)& =& {\displaystyle {\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2}\hfill \\ \hfill & +& {\displaystyle {\int }_t^T{\int }_{{\Gamma }_1\cup {\Gamma }_3}G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\varphi ({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2}\hfill \end{array}$$

(18)

with $\varphi$ depending only on $({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\in {\Gamma }_1\cup {\Gamma }_3\times [0,T)$, i.e.,

$${\displaystyle \varphi ({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})=\left(\frac{1}{2}{\sigma }_1^2{\tilde{S}}_1^2\frac{\partial V}{\partial {\tilde{S}}_1}+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\frac{\partial V}{\partial {\tilde{S}}_2}\right)n_1+\left(\frac{1}{2}{\sigma }_2^2{\tilde{S}}_2^2\frac{\partial V}{\partial {\tilde{S}}_2}+\frac{1}{2}\rho {\sigma }_1{\sigma }_2{\tilde{S}}_1{\tilde{S}}_2\frac{\partial V}{\partial {\tilde{S}}_1}\right)n_2.}$$

4. Boundary Integral Equation

The representation Formula (18) of V holds in the whole domain and can be used once the unknown density $\varphi$ is recovered. To this aim, the idea is to consider the application of the representation formula at the barriers, where $\forall (S_1,S_2,t)\in {\Gamma }_1\cup {\Gamma }_3\times [0,T)$ the option value is trivial, i.e., $V(S_1,S_2,t)=0$ and therefore

$$\begin{array}{c}\hfill {\int }_t^T{\int }_{{\Gamma }_1\cup {\Gamma }_3}G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})\varphi ({\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2=\\ \hfill -{\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G(S_1,S_2,t;{\tilde{S}}_1,{\tilde{S}}_2,T)d{\tilde{S}}_{1}d{\tilde{S}}_2.\end{array}$$

(19)

This Boundary Integral Equation (BIE) in the unknown $\varphi$ is a Volterra integral equation of first kind and can be numerically solved as described in the following section.

Remark 4.

This strategy can be extended to higher dimensions, denoting by Γ the hyperplane defined by the barrier conditions. The representation formula still holds: $\forall (S_1,\dots ,S_n,t)\in \Omega \times [0,T)$

$$\begin{array}{ccc}\hfill V(S_1,\dots ,S_n,t)& =& {\displaystyle {\int }_{\Omega }V({\tilde{S}}_1,\dots ,{\tilde{S}}_n,T)G(S_1,\dots ,S_n,t;{\tilde{S}}_1,\dots ,{\tilde{S}}_n,T)d{\tilde{S}}_1\dots d{\tilde{S}}_n}\hfill \\ & +& {\displaystyle {\int }_t^T{\int }_{\Gamma }G(S_1,\dots ,S_n,t;{\tilde{S}}_1,\dots ,{\tilde{S}}_n,\tilde{t})\varphi ({\tilde{S}}_1,\dots ,{\tilde{S}}_n,\tilde{t})d\tilde{t}d{\tilde{S}}_1\dots d{\tilde{S}}_n}\hfill \end{array}$$

(20)

together with the BIE with knock-out condition: $\forall (S_1,S_2,t)\in \Gamma \times [0,T)$

$$\begin{array}{c}\hfill {\displaystyle {\int }_t^T{\int }_{\Gamma }G(S_1,\dots ,S_n,t;{\tilde{S}}_1,\dots ,{\tilde{S}}_n,\tilde{t})\varphi ({\tilde{S}}_1,\dots ,{\tilde{S}}_n,\tilde{t})d\tilde{t}d{\tilde{S}}_1,\dots ,d{\tilde{S}}_n=}\\ \hfill {\displaystyle -{\int }_{\Omega }V({\tilde{S}}_1,\dots ,{\tilde{S}}_n,T)G(S_1,\dots ,S_n,t;{\tilde{S}}_1,\dots ,{\tilde{S}}_n,T)d{\tilde{S}}_{1}d{\tilde{S}}_2}\end{array}$$

(21)

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.

5. Numerical Approximation by Collocation BEM

We introduce

$$\tilde{\varphi }(S_1,S_2,t):=\sum _{n=1}^{N_{\Delta t}}\sum _{m=1}^{M_{\Delta S}}{\alpha }_{nm}{\overline{\varphi }}_m(S_1,S_2){\widehat{\varphi }}_n\left(t\right)$$

(22)

for the numerical approximation of $\varphi$ by:

• piecewise constant functions in time

  $${\widehat{\varphi }}_n\left(t\right):=H[t-t_{n-1}]-H[t-t_n]n=1,\cdots ,N_{\Delta t}$$

  defined by the Heaviside functions $H[·]$ on a uniform decomposition of the time interval $[0,T]$:

  $$\Delta t=\frac{T}{N_{\Delta t}},N_{\Delta t}\in {\mathbb{N}}^{+}{t}_n=n{N}_{\Delta t},n=0,\cdots ,N_{\Delta t},$$
• piecewise constant functions in space ${\overline{\varphi }}_m\left(S\right)={\overline{\varphi }}_m(S_1,S_2)$, $m=1,\dots ,{\overline{M}}_{\Delta S},{\overline{M}}_{\Delta S}+1,\dots ,M_{\Delta S}$ defined on a decomposition ${\mathcal{T}}_1=\{e_1,\dots ,e_{{\overline{M}}_{\Delta S}}\}$ on ${\Gamma }_1$ and ${\mathcal{T}}_2=\{e_{{\overline{M}}_{\Delta S}+1},$$\dots ,e_{M_{\Delta S}}\}$ on ${\Gamma }_2$ constituted by $M_{\Delta S}$ segments such that $\mathrm{length}\left(e_i\right)\le \Delta S$ with $e_i\cap e_j=\varnothing$ if $i\ne j$. 

The related unknown array

$$\mathbf{\alpha }={\left({\mathbf{\alpha }}^{\left(1\right)}{\mathbf{\alpha }}^{\left(1\right)}{\mathbf{\alpha }}^{\left(2\right)}\dots {\mathbf{\alpha }}^{\left(N_{\Delta t}\right)}\right)}^{\top }{\mathbf{\alpha }}^{(\ell )}={\left({\alpha }_{\ell 1},\cdots ,{\alpha }_{\ell M_{\Delta S}}\right)}^{\top },\ell =1,\dots ,N_{\Delta t}$$

can be computed by solving a linear system

$$\mathbb{M}\mathbf{\alpha }=\mathbf{\beta },$$

(23)

of $N_{\Delta t}\times M_{\Delta S}$ equations obtained by collocating the BIE (19) at the midpoints of time decomposition interval

$${\overline{t}}_n=\frac{t_{n-1}+t_n}{2},n=1,\dots ,N_{\Delta t}$$

(24)

and at the midpoints of space decomposition segments

$${\overline{S}}_m=({\overline{S}}_1^{\left(m\right)},{\overline{S}}_2^{\left(m\right)}),m=1,\dots ,M_{\Delta S}.$$

(25)

Note that the matrix $\mathbb{M}$ has a block upper triangular Toeplitz structure

$$\left(\begin{array}{ccccc}{\mathbb{M}}^{\left(1\right)}& {\mathbb{M}}^{\left(2\right)}& {\mathbb{M}}^{\left(3\right)}& \dots & {\mathbb{M}}^{\left(N_{\Delta t}\right)}\\ 0& {\mathbb{M}}^{\left(1\right)}& {\mathbb{M}}^{\left(2\right)}& \dots & {\mathbb{M}}^{(N_{\Delta t}-1)}\\ 0& 0& {\mathbb{M}}^{\left(1\right)}& \dots & {\mathbb{M}}^{(N_{\Delta t}-2)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & {\mathbb{M}}^{\left(1\right)}\end{array}\right)$$

(26)

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 difference $\tilde{t}-t$ and, by consequence, each block depends only on $t_n-{\overline{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

$$\left\{\begin{array}{cc}{\mathbf{z}}^{(\ell )}={\mathbf{\beta }}^{(\ell )}-\sum _{j=2}^{N_{\Delta t}-\ell +1}{\mathbb{M}}^{\left(j\right)}{\mathbf{\alpha }}^{(\ell +j-1)}\hfill & \\ {\mathbb{M}}^{\left(1\right)}{\mathbf{\alpha }}^{(\ell )}={\mathbf{z}}^{(\ell )}.\hfill \end{array}\right.\ell =N_{\Delta t},\cdots ,1$$

where only the non singular block ${\mathbb{M}}^{\left(1\right)}$ needs to be inverted.

We can get the expression of a general entry of the matrix block ${\mathbb{M}}^{(\ell )},\ell =1,\dots ,N_{\Delta t},$ introducing, in the lhs of (19), a space-time piecewise constant approximation (22) and collocation points (24) and (25): $\forall \overline{m},m=1\dots ,M_{\Delta S},\forall \overline{n}=1\dots ,N_{\Delta t},\forall n=\overline{n}\dots ,N_{\Delta t}$

$$\begin{array}{ccc}\hfill {\mathbb{M}}_{\overline{m}m}^{(\ell )}={\mathbb{M}}_{\overline{m}m}^{\overline{n}n}={\mathbb{M}}_{\overline{m}m}^{(n-\overline{n}+1)}& =& {\displaystyle {\int }_{{\overline{t}}_{\overline{n}}}^T{\int }_{e_m}G({\overline{S}}_1^{\left(\overline{m}\right)},{\overline{S}}_2^{\left(\overline{m}\right)},{\overline{t}}_{\overline{n}};{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t}){\overline{\varphi }}_m({\tilde{S}}_1,{\tilde{S}}_2){\widehat{\varphi }}_n\left(\tilde{t}\right)d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2}\hfill \\ & =& {\displaystyle {\int }_{max(t_{n-1},{\overline{t}}_{\overline{n}})}^{t_n}{\int }_{e_m}G({\overline{S}}_1^{\left(\overline{m}\right)},{\overline{S}}_2^{\left(\overline{m}\right)},{\overline{t}}_n;{\tilde{S}}_1,{\tilde{S}}_2,\tilde{t})d\tilde{t}d{\tilde{S}}_{1}d{\tilde{S}}_2}\hfill \end{array}$$

(27)

and, collocating the rhs of (19) in the same way, we can write

$${\beta }_{\overline{m}}^{\ell }={\beta }_{\overline{m}}^{\overline{n}}=-{\int }_{\Omega }V({\tilde{S}}_1,{\tilde{S}}_2,T)G({\overline{S}}_1^{\left(\overline{m}\right)},{\overline{S}}_2^{\left(\overline{m}\right)},{\overline{t}}_{\overline{n}};{\tilde{S}}_1,{\tilde{S}}_2,T)d\tilde{S}.$$

(28)

Once the array $\mathbf{\alpha }$ of coefficients in (22) are computed by solving the linear system (23), the unknown functional $\varphi$ in (18) can be replaced by its approximation $\tilde{\varphi }$, 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 $\Omega$ 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\to \tilde{t}$ (look at (5)) and possible non-smoothness of the payoff function like, for example, in the case of Basket options.

6. 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

$$V(S_1,S_2,T)=max(S_1+S_2-K,0)$$

with data suggested in [12] and listed in

Table 1. Double knock-out basket call option data.

K	T	r	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	$\mathit{\rho }$	${\mathit{B}}_1$	${\mathit{B}}_2$
1	1	0.05	0.25	0.25	0.7	1	2

. 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 $\Delta S=0.125\sqrt{2}$ and $\Delta 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.

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. Comparison between FEM and SABO: option value obtained at $S_1=0.5,S_2=1,t=0$ (upper tables) and CPU time of computation in seconds (lower tables). Reference value: 0.306264.

SABO						FEM
	$N_{\Delta t}$	2	4	8	16	$N_{\Delta t}$	25	50	100
$M_{\Delta S}$
2		0.304101	0.305635	0.305527	0.305501	$R_0$	0.293627	0.293663	0.293662
4		0.305242	0.306735	0.306631	0.306586	$R_1$	0.304109	0.304129	0.304134
8		0.304903	0.306424	0.306320	0.306274	$R_2$	0.305875	0.305891	0.305894
16		0.304916	0.306429	0.306329	0.306282	$R_3$	0.306171	0.306185	0.306188
						$R_4$	0.306229	0.306243	0.306246
CPU time SABO						CPU time FEM
	$N_{\Delta t}$	2	4	8	16	$N_{\Delta t}$	25	50	100
$M_{\Delta S}$
2		2.1$\times {10}^{-1}$	5.8$\times {10}^{-1}$	8.6$\times {10}^{-1}$	1.1$\times {10}^0$	$R_0$	7.2$\times {10}^{-1}$	1.4$\times {10}^0$	3.4$\times {10}^0$
4		5.1$\times {10}^{-1}$	7.8$\times {10}^{-1}$	1.2$\times {10}^0$	2.1$\times {10}^0$	$R_1$	4.3$\times {10}^0$	8.3$\times {10}^0$	1.5$\times {10}^1$
8		9.2$\times {10}^{-1}$	1.5$\times {10}^0$	2.5$\times {10}^0$	4.5$\times {10}^0$	$R_2$	8.1$\times {10}^0$	1.3$\times {10}^1$	2.2$\times {10}^1$
16		2.1$\times {10}^0$	4.1$\times {10}^0$	6.8$\times {10}^0$	1.3$\times {10}^0$	$R_3$	3.0$\times {10}^1$	3.6$\times {10}^1$	4.8$\times {10}^1$
						$R_4$	6.4$\times {10}^2$	6.7$\times {10}^2$	7.1$\times {10}^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_{\Delta S}=N_{\Delta 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.

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 ${\Gamma }_4$ (or equivalently ${\Gamma }_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:

$$\begin{array}{cc}\hfill C(S,t)& {\displaystyle =S\sum _{n=-\infty }^{\infty }\left\{{\left(\frac{B_2^n}{B_1^n}\right)}^{\frac{2r}{{\sigma }^2}+1}\left(\mathcal{N}\left[d_1\right]-\mathcal{N}\left[d_2\right]\right)-{\left(\frac{B_1^{n+1}}{B_2^{n}S}\right)}^{\frac{2r}{{\sigma }^2}+1}\left(\mathcal{N}\left[d_3\right]-\mathcal{N}\left[d_4\right]\right)\right\}}\hfill \\ & {\displaystyle -K{e}^{-r(T-t)}\sum _{n=-\infty }^{\infty }\left\{{\left(\frac{B_2^n}{B_1^n}\right)}^{\frac{2r}{{\sigma }^2}-1}\left(\mathcal{N}[d_1-\sigma \sqrt{T-t}]-\mathcal{N}[d_2-\sigma \sqrt{T-t}]\right)\right.}\hfill \\ & {\displaystyle \left.-{\left(\frac{B_1^{n+1}}{B_2^{n}S}\right)}^{\frac{2r}{{\sigma }^2}-1}\left(\mathcal{N}[d_3-\sigma \sqrt{T-t}]-\mathcal{N}[d_4-\sigma \sqrt{T-t}]\right)\right\}}\hfill \end{array}$$

(29)

$$\begin{array}{cccccc}\hfill d_1& =& {\displaystyle \frac{ln\left(\frac{S{B}_2^{2n}}{K{B}_1^{2n}}\right)+\left(r+\frac{{\sigma }_2}{2}(T-t)\right)}{\sigma \sqrt{T-t}}}\hfill & \hfill d_2& =& {\displaystyle \frac{ln\left(\frac{S{B}_2^{2n-1}}{B_1^{2n}}\right)+\left(r+\frac{{\sigma }_2}{2}(T-t)\right)}{\sigma \sqrt{T-t}}}\hfill \\ \hfill d_3& =& {\displaystyle \frac{ln\left(\frac{B_1^{2n+2}}{KS{B}_2^{2n}}\right)+\left(r+\frac{{\sigma }_2}{2}(T-t)\right)}{\sigma \sqrt{T-t}}}\hfill & \hfill d_4& =& {\displaystyle \frac{ln\left(\frac{B_1^{2n+2}}{S{B}_2^{2n+1}}\right)+\left(r+\frac{{\sigma }_2}{2}(T-t)\right)}{\sigma \sqrt{T-t}}}\hfill \end{array}$$

with $\mathcal{N}$ the standard normal cumulative distribution function.

These “boundary conditions” on ${\Gamma }_2$ and ${\Gamma }_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 $\rho =0.1$ on the top, highly correlated assets with $\rho =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

$$V(S_1,S_2,T)=max(K-S_2,0)\mathrm{if}{\mathrm{S}}_1\in [\mathrm{B},+\infty ]$$

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_{\Delta S}=N_{\Delta t}=10$. In

Table 3. Stabilization of value $V(S_1,S_2,0)$ at $(S_1,S_2)=(2,1)$ with the refinement of discretization parameters ${\Delta }_S$ and ${\Delta }_t$.

	${\mathit{M}}_{\mathbf{\Delta }\mathit{S}}$	5	10	20
${\mathit{N}}_{\mathbf{\Delta }\mathit{t}}$
5		0.89692	0.89731	0.89731
10		0.89679	0.89719	0.89719
20		0.89675	0.89715	0.89716
40		0.89673	0.89714	0.89714

, 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_{\Delta S}$ and $N_{\Delta t}$ in (22).

7. 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.