 # Numerical valuation of Bermudan basket options via partial differential equations

We study the principal component analysis (PCA) based approach introduced by Reisinger & Wittum (2007) for the approximation of Bermudan basket option values via partial differential equations (PDEs). This highly efficient approximation approach requires the solution of only a limited number of low-dimensional PDEs complemented with optimal exercise conditions. It is demonstrated by ample numerical experiments that a common discretization of the pertinent PDE problems yields a second-order convergence behaviour in space and time, which is as desired. It is also found that this behaviour can be somewhat irregular, and insight into this phenomenon is obtained.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

This paper deals with the valuation of Bermudan basket options. Basket options have a payoff depending on a weighted average of different assets. Semi-closed analytic valuation formulas are generally lacking in the literature for these options. Consequently, research into efficient and stable methods for approximating their fair values is of much interest. The valuation of basket options gives rise to multidimensional time-dependent partial differential equations. Here the spatial dimension equals the number of different assets in the basket. When is large, it is well-known that this leads to a highly challenging numerical task. In the present paper we shall consider Bermudan-style basket options and investigate a principal component analysis based approach introduced by Reisinger & Wittum  and subsequently studied in e.g. [3, 4, 5] that renders this task feasible.

A European-style basket option is a financial contract that provides the holder the right, but not the obligation, to buy or sell a given weighted average of assets at a specified future date for a specified price . Parameter is called the maturity time and the strike price of the option. In this paper we assume the well-known Black–Scholes model. Then the asset prices () follow a multidimensional geometric Brownian motion, under the risk-neutral measure, given by the system of stochastic differential equations

 dSiτ=rSiτdτ+σiSiτdWiτ(0<τ≤T, 1≤i≤d).

Here is time, with being the time of inception of the option, is the given risk-free interest rate, () are the given volatilities and () is a multidimensional standard Brownian motion with given correlation matrix . Further, initial asset prices () are given. Let be the fair value of a European basket option if at time the -th asset price equals (), where is the time remaining till maturity of the option. From financial mathematics theory it follows that satisfies the -dimensional time-dependent partial differential equation (PDE)

 ∂u∂t(s,t)=12d∑i=1d∑j=1σiσjρijsisj∂2u∂si∂sj(s,t)+d∑i=1rsi∂u∂si(s,t)−ru(s,t) (1.1)

for . PDE (1.1) is also fulfilled whenever for any given , defining a natural boundary condition. In almost all financial applications, the correlation matrix has nonzero off-diagonal entries, and hence, (1.1) contains mixed spatial derivative terms. At maturity time of the option its fair value is known and determined by the particular option contract. This yields the initial condition

 u(s,0)=ϕ(s) (1.2)

for . Here function is the given payoff of the option.

A Bermudan-style basket option is a financial contract that provides the holder the right to buy or sell a given weighted average of assets for a specified price at one from a specified finite set of exercise times with . Let for and . Then the fair value function of a Bermudan basket option satisfies the PDE (1.1), with the natural boundary condition, on each time interval for . Next, the initial condition (1.2) holds and for one has

 u(s,αe)=max(ϕ(s),limt↑αeu(s,t)) (1.3)

whenever . Condition (1.3) stems from the early exercise feature of Bermudan options and represents the optimal exercise condition. Notice that it is nonlinear. In the present paper we shall consider the class of Bermudan basket put options. These have a payoff function of the form

 ϕ(s)=max(K−d∑i=1ωisi,0) (1.4)

with given fixed weights such that .

An outline of the rest of this paper is as follows.

Following Reisinger & Wittum , we first apply in Section 2.1 a useful coordinate transformation to (1.1) by using a spectral decomposition of the pertinent covariance matrix. This leads to a -dimensional time-dependent PDE for a transformed option value function

in which each coefficient is directly proportional to one of the eigenvalues. In Section

2.2 this feature is employed to define a principal component analysis (PCA) based approximation to . The key property of is that it is defined by only a limited number of one- and two-dimensional PDEs. In Section 2.3 a note on the optimal exercise condition is given. Section 2.4 describes a common discretization of the one- and two-dimensional PDE problems by means of finite differences on a suitable nonuniform spatial grid followed by the Brian and Douglas ADI scheme on a uniform temporal grid. In view of the nonsmoothness of the payoff function, cell averaging and backward Euler damping are applied.

The main contribution of our paper is given in Section 3. Extensive numerical experiments are presented where we study in detail the error of the discretization described in Section 2.4 of the PCA-based approximation defined in Section 2.2 for Bermudan basket options. Three financial parameter sets from the literature are considered, with number of assets . A second-order convergence behaviour is observed, which is as desired. It is also found that this behaviour can be somewhat irregular. Additional numerical experiments are performed that yield insight into this phenomenon.

Section 4 contains our conclusions and outlook.

## 2 Approximation approach

### 2.1 Coordinate transformation

In this section the PDE (1.1) for a Bermudan basket option is converted into a form that is the starting point for the solution approach discussed in the subsequent sections. In the following, the elementary functions , , and

are to be taken componentwise whenever they are applied to vectors.

Consider the covariance matrix given by for . Let

be an orthogonal matrix of eigenvectors of

and a diagonal matrix of eigenvalues of such that . As in , consider the coordinate transformation

 x(s,t)=QT(ln(s/K)−b(t)), (2.1)

with and for . Define the function by

 u(s,t)=v(x(s,t),t).

A straightforward calculation shows that satisfies

 ∂v∂t(x,t)=12d∑k=1λk∂2v∂x2k(x,t)−rv(x,t) (2.2)

whenever , , . The PDE (2.2) is a pure diffusion equation, without mixed derivatives, and with a simple reaction term.

It is convenient to perform a second coordinate transformation , which maps the spatial domain onto the -dimensional open unit cube,

 y(x)=1πarctan(x)+12. (2.3)

Define the function by

 v(x,t)=w(y(x),t).

Then it is readily verified that

 ∂w∂t(y,t)=d∑k=1λk[p(yk)∂2w∂y2k(y,t)+q(yk)∂w∂yk(y,t)]−rw(y,t) (2.4)

whenever , , with

 p(η)=12π2sin4(πη),q(η)=1πsin3(πη)cos(πη)for η∈R.

Clearly, the PDE (2.4) is a convection-diffusion-reaction equation without mixed derivative terms. Define the function by

 ψ(y,t)=ϕ(Kexp[Qx+b(t)])  with  x=tan[π(y−12)] (2.5)

whenever , . Then for (2.4) one has the initial condition

 w(y,0)=ψ(y,0) (2.6)

together with the optimal exercise condition

 w(y,αe)=max(ψ(y,αe),limt↑αew(y,t)) (2.7)

for and . At the boundary of the spatial domain we shall consider a Dirichlet condition. In Appendix A the details of its derivation are provided, where the minor Assumption A.1 on the matrix is made. For any given such that the entries of the -th column of are all strictly positive there holds

 w(y,t)=Ke−r(t−αe−1) (2.8)

whenever with and , . On the complementary part of a homogeneous Dirichlet condition is valid.

### 2.2 Principal component analysis based approximation

Let the eigenvalues of be ordered such that . In financial applications it often holds that is significantly larger than the other eigenvalues. Motivated by this observation, Reisinger & Wittum  proposed a principal component analysis (PCA) based approximation of the exact solution to the multidimensional PDE (2.4). To this purpose, consider also as a function of the eigenvalues and write with . Set

 ˆλ=(λ1,0,…,0)T  and  δλ=λ−ˆλ=(0,λ2,…,λd)T.

Assuming sufficient smoothness of , a first-order Taylor expansion at yields

 w(y,t;λ)≈w(y,t;ˆλ)+d∑l=2δλl∂w∂λl(y,t;ˆλ). (2.9)

The partial derivative (for ) can be approximated by a forward finite difference,

 ∂w∂λl(y,t;ˆλ)≈w(y,t;ˆλ+δλlel)−w(y,t;ˆλ)δλl, (2.10)

where denotes the -th standard basis vector in . Combining (2.9) and (2.10), gives

 w(y,t;λ)≈w(y,t;ˆλ)+d∑l=2[w(y,t;ˆλ+δλlel)−w(y,t;ˆλ)].

Write

 w(1)(y,t)=w(y,t;ˆλ)  and  w(1,l)(y,t)=w(y,t;ˆλ+δλlel).

 w(y,t)≈˜w(y,t)=w(1)(y,t)+d∑l=2[w(1,l)(y,t)−w(1)(y,t)] (2.11)

whenever , , . By definition, satisfies the PDE (2.4) with being set to zero for all and satisfies (2.4) with being set to zero for all , which is completed by the same initial condition, optimal exercise condition and boundary condition as for , discussed above.

In financial applications one is often interested in the option value at inception in the single point , where is the vector of known asset prices. Let

 Y0=y(x(S0,T))∈(0,1)d

denote the corresponding point in the -domain. Then can be obtained by solving a one-dimensional PDE on the line segment in the -domain that is parallel to the -axis and passes through . In other words, can be fixed at the value whenever . Next, (for ) can be obtained by solving a two-dimensional PDE on the plane segment in the -domain that is parallel to the -plane and passes through . Thus, in this case, can be fixed at the value whenever .

In view of the above key observation, computing the PCA-based approximation (2.11) for requires solving just 1 one-dimensional PDE and two-dimensional PDEs. This clearly yields a main computational advantage compared to solving the full -dimensional PDE whenever is large. Moreover, the different terms in (2.11) can be computed in parallel.

Reisinger & Wissmann  have given a rigorous analysis of the error in the PCA-based approximation relevant to European basket options. Under a mild assumption on the payoff function , they proved that in the maximum norm.

### 2.3 A note regarding the optimal exercise condition

Let and write . Let , which forms the intersection of and . By the optimal exercise condition (2.7), the natural approximation to at based on is

 w(y,αe) ≈ max(ψe(y),limt↑αe˜w(y,t)) = limt↑αemax(ψe(y),˜w(y,t)) = limt↑αemax(ψe(y),w(1)(y,t)+d∑l=2[w(1,l)(y,t)−w(1)(y,t)]).

On the other hand, by construction of and (), we have

 w(y,αe) ≈ ˜w(y,αe) = w(1)(y,αe)+d∑l=2[w(1,l)(y,αe)−w(1)(y,αe)] = limt↑αe(max(ψe(y),w(1)(y,t))+d∑l=2[max(ψe(y),w(1,l)(y,t))−max(ψe(y),w(1)(y,t))]).

It may hold that

 ˜w(y,αe)≠max(ψe(y),limt↑αe˜w(y,t)),

and hence, the PCA-based approximation does not satisfy the optimal exercise condition. A further investigation into this will be the subject of future research.

### 2.4 Discretization

To arrive at the values and (for ) in the approximation of we perform finite difference discretization of the pertinent one- and two-dimensional PDEs on a (Cartesian) nonuniform spatial grid, followed by a suitable implicit time discretization. Let and . Note that the point in the -domain corresponds to the point in the -domain if . For any given a nonuniform mesh in the -th spatial direction is defined by (see e.g. )

 yk,i=φ(ξi)  with  ξi=ξmin+iΔξ, Δξ=ξmax−ξminm+1(i=0,1,…,m+1),

with

 φ(ξ)=κ0+κ1sinh(ξ)(ξmin≤ξ≤ξmax)

and

 ξmin=−sinh−1(κ0/κ1)  and  ξmax=sinh−1((1−κ0)/κ1).

Remark that since . The parameter controls the fraction of mesh points that lie in the neighborhood of

. We make the heuristic choice

. The above mesh is smooth in the sense that there exist constants (independent of , ) such that the mesh widths satisfy

 C0Δξ≤Δyk,i≤C1Δξ  and  |Δyk,i+1−Δyk,i|≤C2(Δξ)2.

The spatial derivatives in (2.4) are discretized using central finite difference schemes. Let be any given smooth function, let be any given smooth mesh and denote the mesh widths by . Then second-order approximations to the first and second derivatives are given by

 f′(ηi)≈βi,−1f(ηi−1)+βi,0f(ηi)+βi,1f(ηi+1),f′′(ηi)≈γi,−1f(ηi−1)+γi,0f(ηi)+γi,1f(ηi+1),

with

 βi,−1=−hi+1hi(hi+hi+1) , βi,0=hi+1−hihihi+1 , βi,1=hihi+1(hi+hi+1) ,

and

 γi,−1=2hi(hi+hi+1) , γi,0=−2hihi+1 , γi,1=2hi+1(hi+hi+1) .

The above two finite difference formulas are applied with for and .

Semidiscretization of the PDE for on the plane segment

leads to a system of ordinary differential equations (ODEs)

 W′(t)=(λ1A1+λlAl)W(t) (2.12)

for , . Here is a vector of dimension and , are given matrices that are tridiagonal (possibly up to permutation) and correspond to, respectively, the first and the -th spatial direction. The ODE system (2.12) is completed by an initial condition

 W(0)=Ψ0

and, for , an optimal exercise condition

 W(αe)=max(Ψe,limt↑αeW(t)).

Here the vector is determined by the function on for . The maximum of any given two vectors is to be taken componentwise.

The payoff function given by (1.4) is continuous but not everywhere differentiable, and hence, this also holds for the function given by (2.5). It is well-known that the nonsmoothness of the payoff function can have an adverse impact on the convergence behaviour of the spatial discretization. To alleviate this, we employ cell averaging near the points of nonsmoothness in defining the initial vector , see e.g. .

For the temporal discretization of the ODE system (2.12) a standard Alternating Direction Implicit (ADI) method is applied. Consider a given step size with integer and define temporal grid points for . Assume that for some integer whenever . Let and

 N={n1,n2,…,nE−1}.

Application of the familiar second-order Brian and Douglas ADI scheme for two-dimensional PDEs leads to an approximation that is successively defined for by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Z0=Wn−1+Δt(λ1A1+λlAl)Wn−1,Z1=Z0+12Δtλ1A1(Z1−Wn−1),Z2=Z1+12ΔtλlAl(Z2−Wn−1),Wn=Z2  (if n∉N)andWn=max(Ψe,Z2)  (if n=ne∈N). (2.13)

In the scheme (2.13) a forward Euler predictor stage is followed by two implicit but unidirectional corrector stages, which serve to stabilize the predictor stage. The two linear systems in each time step can be solved very efficiently by using a priori factorizations of the pertinent matrices. As for the spatial discretization, also the convergence behaviour of the temporal discretization can be adversely effected by the nonsmooth payoff function. To remedy this, backward Euler damping (or Rannacher time stepping) is applied at initial time as well as at each exercise date, that is, with , the time step from to , is replaced by two half steps of the backward Euler method for .

Finally, discretization of the PDE for on the line segment is performed completely analogously to the above. Then a semidiscrete system is obtained with a vector of dimension and an tridiagonal matrix. Temporal discretization is done using the Crank–Nicolson scheme with backward Euler damping.

## 3 Numerical experiments

In this section we investigate by ample numerical experiments the error of the discretization described in Section 2.4 of the PCA-based approximation defined in Section 2.2. We consider three parameter sets for the basket option and the underlying asset price model.

Set A is given by Reisinger & Wittum  and has , , , and

 (ρij)di,j=1=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1.000.790.820.910.840.791.000.730.800.760.820.731.000.770.720.910.800.771.000.900.840.760.720.901.00⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

The eigenvalues of the corresponding covariance matrix are

Hence, is clearly dominant.

Sets B and C are taken from Jain & Oosterlee  and have dimensions and , respectively. Here , , and , and whenever . Sets B and C have and , respectively, and . Thus is also dominant for these parameter sets. It can be shown that for all three Sets A, B, C the matrix of eigenvectors of satisfies Assumption A.1.

We consider a Bermudan basket option with exercise times () and study the absolute error in the discretization of at the point with . For comparison, also the European basket option is included in the experiments. The number of time steps is taken as for the European option and for the Bermudan option.

Table 1 provides reference values for , which have been computed by choosing . In the case of Set A, Reisinger & Wittum  obtain the approximation for the European basket option. In the case of Sets B and C, Jain & Oosterlee  obtain, using the stochastic grid bundling method, the approximations and , respectively, for the Bermudan basket option. Clearly, these three approximations from the literature agree well with our corresponding values for given in Table 1.

Figure 1 displays the absolute error in the discretization of versus for all . Here both the European and Bermudan basket options are considered for all three parameter sets A, B, C. The favourable result is observed that the discretization error is always bounded from above by with a moderate constant , which is as desired.

For the European option and Set A, the error drop in the (less important) region is somewhat surprising, but it is easily explained from a change of sign in the error. Except for this, in the case of the European basket option, the error behaviour is always found to be regular and second-order.

For the Bermudan basket option the observed error behaviour is less regular, in particular in the interesting region of large values . To gain more insight into this phenomenon, we have computed separately the discretization error for the leading term and for the correction term in , see (2.11). Reference values for the leading term are given in Table 2.

The obtained result is shown in Figure 2, where dark squares indicate the error for the leading term and light circles the error for the correction term. It is clear that in all six cases the error for the leading term behaves regularly and the error for the correction term is small compared to this (for Set A if , as above). For the Bermudan basket option, however, the behaviour of the discretization error for the correction term is rather irregular. A subsequent study shows that for any given the error is always very close to the error , which is as expected, but the difference can be both positive and negative, leading to an irregular behaviour of . This is exacerbated when summing these differences up over . Hence, the irregular behaviour of the error for the correction term can adversely affect the regular behaviour of the error for the leading term. We remark that this has been observed in many other experiments we performed for the Bermudan basket option, for example for other points , for other numbers of exercise times , for other dimensions and for other covariance matrices , having . It is our aim of future research to find a remedy for this phenomenon.

## 4 Conclusions

In this paper we have investigated the PCA-based approach by Reisinger & Wittum  for the valuation of Bermudan basket options. This approximation approach is highly effective as it requires the solution of only a limited number of low-dimensional PDEs, supplemented with optimal exercise conditions. By numerical experiments the favourable result is shown that a common discretization of these PDE problems leads to a second-order convergence behaviour in space and time. It is also observed that this convergence behaviour can be somewhat irregular. Insight into this phenomenon is obtained by regarding the total discretization error as a superposition of discretization errors for the leading term and the correction term. Our aim for future research is to determine a suitable remedy for it. Another topic for future research concerns a rigorous analysis of the error in the PCA-based approximation for Bermudan basket options. The results obtained by Reisinger & Wissmann , relevant to European basket options, will be important here.

## References

•  K.J. in ’t Hout. Numerical Partial Differential Equations in Finance Explained. Financial Engineering Explained. Palgrave Macmillan UK, 2017.
•  S. Jain and C.W. Oosterlee. The stochastic grid bundling method: efficient pricing of Bermudan options and their Greeks. Appl. Math. Comp., 269:412–431, 2015.
•  C. Reisinger and R. Wissmann. Numerical valuation of derivatives in high-dimensional settings via partial differential equation expansions. J. Comp. Finan., 18(4):95–127, 2015.
•  C. Reisinger and R. Wissmann. Error analysis of truncated expansion solutions to high-dimensional parabolic PDEs. ESAIM: M2AN, 51(6):2435–2463, 2017.
•  C. Reisinger and R. Wissmann. Finite difference methods for medium- and high-dimensional derivative pricing PDEs. In High-Performance Computing in Finance: Problems, Methods, and Solutions, pages 175–196. 2018.
•  C. Reisinger and G. Wittum. Efficient hierarchical approximation of high-dimensional option pricing problems. SIAM J. Sci. Comp., 29(1):440–458, 2007.

## Appendix A Dirichlet boundary condition for (2.4)

Consider the following minor assumption on the matrix  of eigenvectors of the covariance matrix .

###### Assumption A.1

Each column of satisfies one of the following two conditions:

• all its entries are strictly positive;

• it has both a strictly positive and a strictly negative entry.

Then we have

###### Lemma A.2

Let the function be given by (2.5) with defined by (1.4). Let , and with fixed whenever . If the -th column of satisfies (A.1.a), then as . If the -th column of satisfies (A.1.b), then as . Finally, as .

Proof  Let and , so that .

Suppose first . Then . If the -th column of satisfies condition (A.1.a), then all entries of tend to . Consequently, all entries of tend to zero and . If the -th column of satisfies condition (A.1.b), then the entries of go to either or with at least one entry that tends to . It follows that the entries of go to either zero or with at least one entry that tends to , and therefore .

Suppose next . Then and the entries of go to either or with at least one entry that tends to . Hence, .

For any given the diffusion and convection coefficients and in (2.4) vanish as or . Accordingly, (2.4) is also satisfied on each boundary part

 {y:y=(y1,y2,…,yd)T with yk=δ and yj∈(0,1) whenever j≠k}

for . Also the initial condition (2.6) and optimal exercise condition (2.7) hold on each such boundary part, upon taking the relevant limit value for given by Lemma A.2. On each part where this limit value equals , the solution (2.8) is obtained, and on each part where the limit value equals zero, the zero solution holds. This yields the Dirichlet boundary condition for the PDE (2.4) stated in Section 2.1.