An ETD Method for Vulnerable American Options

1. Introduction

The Global Financial Crisis of 2007–2008 has raised concerns regarding the exposure to credit default risk in financial derivatives traded in over-the-counter (OTC) markets. In the absence of an organized exchange, option holders face the risk of default, leading to the emergence of "vulnerable options" that carry lower values compared to their non-vulnerable counterparts due to the inherent default risk.

Early studies by Johnson and Stulz [1] proposed pricing formulas for vulnerable European options under the assumption that the option represents the only liability of the involved party. Subsequent research by Klein [2] extended this work by considering the recovery of nominal claims in default and the correlation between the issuer’s asset and the underlying option asset. Further advancements [3] incorporated the Vasicek model and introduced a default boundary dependent on the option’s value. In [4], the model has been extended to accommodate default before option maturity, while Ref. [5] explores vulnerable options in incomplete markets with fluctuating interest rates.

Recently, the focus has shifted towards the pricing of American options, which offer the holder the right to exercise at any point before the expiration date. Klein and Yang [6] ventured into the realm of American vulnerable options, but a closed-form formula remained elusive. Later, Ref. [7] provided a semi-analytical solution specifically for standalone vulnerable American put options, with a focus on the early exercise boundary.

The difficulty arises in pricing vulnerable American options due to the significant discontinuity in their payoff function. Current models for these options are often based on a binomial tree method, yielding option prices and critical stock prices in tabulated form [7]. Ref. [8] utilized the two-point Geske and Johnson method to assess vulnerable American options, delivering analytical pricing formulas under risk-neutral measures. The study also delved into the price sensitivity of these options to correlations between the underlying and option writer’s assets. Further, Refs. [9,10] offered models incorporating jump-diffusion processes and correlated credit risk, respectively.

Pricing vulnerable American options poses substantial challenges due to the presence of early-exercise opportunities and the mixed derivative term in the governing nonlinear partial differential equation (PDE), which complicates numerical schemes and may result in solution instabilities [11]. To overcome these difficulties, this paper presents a novel numerical methodology that integrates multiple techniques. The penalty method, see [12] and references therein, to our knowledge, previously unexplored in the context of vulnerable American options, is employed to handle early-exercise opportunities. The resulting nonlinear PDE is then transformed into a canonical form through mixed derivative elimination [13].

The transformed PDE is solved using the Exponential Time Differencing (ETD) method [14], designed for solving ordinary differential equations (ODEs) or systems of ODEs. This method utilizes exponential integrators, known for their effectiveness in handling stiff systems characterized by solution modes with divergent decay rates. ETD techniques excel at capturing both rapid transient dynamics and slower dynamics accurately. In the context of vulnerable American options, a suitable approach is required for applying the ETD method to PDEs. The method of lines is employed to discretize space variables, such as the underlying asset price and the firm’s asset value, thereby converting the pricing PDE into a system of ODEs.

The ETD technique is then applied to this system of ODEs, efficiently yielding a stable and accurate solution for the option price. This combined approach enables the handling of the complex mixed-derivative term, which often poses numerical challenges in traditional methods. To validate the proposed method, extensive numerical experiments are conducted, comparing the results with existing approaches. Furthermore, numerical stability and convergence rate analyses reinforce the efficiency and precision of the ETD method. The findings demonstrate that the proposed method offers not only computational efficiency but also remarkable accuracy in pricing vulnerable American options.

The remainder of this paper is structured as follows. Section 2 outlines the formulation of the free-boundary PDE problem for vulnerable American options. In Section 3, the proposed methodology is detailed, encompassing the mixed-derivative removal transformation, the semi-discretization of the PDE, and the utilization of the ETD method to solve the resulting system of ODEs. Noteworthy focus is given to the solution of the default case. Section 4 presents the numerical results, including a numerical analysis of stability and convergence, as well as a comprehensive comparison with existing methods. Finally, Section 5 provides concluding remarks.

2. Vulnerable option modelling

Let us consider an American option written on the asset S which follows the risk-neutral process

$$\frac{dS}{S}=(r-q)dt+{\sigma }_{\mathrm{S}}d{Z}_{\mathrm{S}},$$

(1)

where r is the risk-free rate, q is the dividend yield, ${\sigma }_{\mathrm{S}}$ is the constant volatility of the underlying asset S, $Z_{\mathrm{S}}$ is a standard Wiener process.

If the option is traded in the over-the-counter market and there is no guaranty that participants will fulfil their obligations, the option writers become vulnerable to a counter-party credit risk, and corresponding options are known as vulnerable.

Johnson and Stulz [1] were pioneers in studying vulnerable options. They assumed that default may happen when the option price is greater than the value of the option writer’s assets. This model was extended by Klein and Inglis in [3] by allowing the option writer to hold other liabilities, but only changes in the value of the writer’s assets, including underlying asset of the option, can lead to financial distress. Later on, in [6] Klein and Yang derived the PDE formulation for vulnerable American option price considered in present study.

Let V denote the market value of the writer’s assets, which is correlated with the underlying asset price S. The risk-neutral process for V is given as follows

$$\frac{dV}{V}=rdt+{\sigma }_{\mathrm{V}}d{Z}_{\mathrm{V}},$$

(2)

where ${\sigma }_{\mathrm{V}}$ is the instantaneous standard deviation of the return on the writer’s assets, $Z_{\mathrm{V}}$ is a standard Wiener process. The instantaneous correlation coefficient between $d{Z}_{\mathrm{S}}$ and $d{Z}_{\mathrm{V}}$ is $\rho$.

Let $\tau =T-t$ be the time to maturity, then the vulnerable American put option price $P(\tau ,S;V)$ is the solution of the following nonlinear PDE

$$\frac{\partial P}{\partial \tau }=\frac{{\sigma }_{\mathrm{S}}^2}{2}{S}^2\frac{{\partial }^{2}P}{\partial S^2}+\frac{{\sigma }_{\mathrm{V}}^2}{2}{V}^2\frac{{\partial }^{2}P}{\partial V^2}+\rho {\sigma }_{\mathrm{S}}{\sigma }_{\mathrm{V}}\frac{{\partial }^{2}P}{\partial S\partial V}+(r-q)S\frac{\partial P}{\partial S}+rV\frac{\partial P}{\partial V}-rP,$$

(3)

subject to the following initial conditions

$$P(0,S;V)={(S-K)}^{+}\left({\mathbb{1}}_{(V\ge D^{*}+K-S)}+\frac{(1-\alpha )V}{D^{*}+K-S}{\mathbb{1}}_{(V<D^{*}+K-S)}\right),$$

(4)

where ${(·)}^{+}=max\{·,0\}$; $D^{*}$ is the fixed default boundary such that a default occurs if the value of the option writer’s assets $V<D^{*}+{(K-S)}^{+}$. If $V\ge D^{*}+{(K-S)}^{+}$, the entire claim is paid out. Parameter $\alpha$ is the dead-weight cost related with the bankruptcy of the writer, expressed as a percentage of the value of the writer’s assets. Expression (4) is a payoff function of a vulnerable put option with strike price K, and ${\mathbb{1}}_A$ is an indicator function.

The boundary condition when $S\to \infty$ for PDE put problem (3) is established as follows

$$\underset{S\to \infty }{lim}P(\tau ;S,V)=0,\tau >0.$$

(5)

Let us denote the value of corresponding non-vulnerable American option by $P_v(\tau ,S)$, which can be calculated as a solution of the Black-Scholes problem [15],

$$\frac{\partial P_v}{\partial \tau }=\frac{{\sigma }_{\mathrm{S}}^2}{2}{S}^2\frac{{\partial }^2{P}_v}{\partial S^2}+(r-q)S\frac{\partial P_v}{\partial S}-r{P}_v,S>S_f\left(\tau \right),$$

(6)

where $S_f\left(\tau \right)$ is the unknown early exercise boundary, subject to the following initial and boundary conditions

$$\begin{array}{ccc}\hfill P_v(0,S)& =& {(K-S)}^{+},S_f\left(0\right)=K,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill P_v(\tau ,S_f\left(\tau \right))& =& K-S_f\left(\tau \right),\frac{\partial P_v}{\partial S}(\tau ,S_f\left(\tau \right))=-1,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill P_v(\tau ,S)& =& K-S,0\le S<S_f\left(\tau \right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \underset{S\to \infty }{lim}{P}_v(\tau ,S)& =& 0.\hfill \end{array}$$

(10)

If default happens prior to maturity, the option price is calculated as follows

$$P(\tau ,S;V)=(1-\alpha )\frac{V}{D^{*}+P_v(\tau ,S)}{P}_v(\tau ,S),S\le K,V\le D^{*}+P_v(\tau ,S),$$

(11)

where $P_v(\tau ,S)$ is the value of vanilla (non-vulnerable) American option defined by PDE problem (6). Expression (11) defines the amount which the holder will receive if the value of option writer’s assets does not exceed the variable default boundary prior to maturity.

If default does not happen prior to maturity ($V>D^{*}+P_v(\tau ,S)$), then the option price is a solution of PDE (3) with the conditions imposed on the optimal stopping boundary $S_f\left(\tau \right)$ [6]:

$$P(\tau ,S_f\left(\tau \right),V)=K-S_f\left(\tau \right),\frac{\partial P}{\partial S}(\tau ,S_f\left(\tau \right),V)=-1.$$

(12)

Hence, (3) with initial conditions (4) and boundary conditions (5), is a free boundary PDE problem that can be reduced to the fixed domain problem by introducing some penalty term. In the simplest case, the penalty term is defined as follows

$$f\left(P\right)=\mu {({(K-S)}^{+}-P(\tau ,S;V))}^{+},$$

(13)

where $\mu$ is some positive sufficiently large constant. This penalty term guarantees that the solution at any moment before the maturity will not be less than the payoff of the corresponding non-vulnerable option. This penalty term allows a fixed solution domain, removing the free and moving boundary imposed by the early exercise feature of the contract. Note that if $\mu =0$, PDE problem (14) reduces to the vulnerable European option.

Finally, PDE problem to be solved is

$$\frac{\partial P}{\partial \tau }=\frac{{\sigma }_{\mathrm{S}}^2}{2}{S}^2\frac{{\partial }^{2}P}{\partial S^2}+\frac{{\sigma }_{\mathrm{V}}^2}{2}{V}^2\frac{{\partial }^{2}P}{\partial V^2}+\rho {\sigma }_{\mathrm{S}}{\sigma }_{\mathrm{V}}\frac{{\partial }^{2}P}{\partial S\partial V}+(r-q)S\frac{\partial P}{\partial S}+rV\frac{\partial P}{\partial V}-rP+f\left(P\right),$$

(14)

subject to the initial conditions (4), and boundary conditions (11).

3. Numerical solution

The PDE problem (14) with conditions (4) and (11) does not have an analytic solution and has to be solved numerically. In [6] the three-dimensional tree was used. This method is easy to implement, but computationally expensive and time-consuming for obtaining the solution not in one fixed point, but in some domain. Therefore, in present study the method of exponential time differencing proposed in [14] is employed after applying the mixed derivative terms removing transformation [13] and semi-discretization technique [16].

3.1. Mixed derivative terms removing

The presence of mixed derivative terms in a partial differential equation can lead to instabilities and inaccuracies, posing significant challenges for numerical methods. To overcome these issues, various transformation techniques have been developed to remove these mixed derivative terms. One such method, proposed in [13], is based on an $LD{L}^T$ transformation of the correlation matrix, which avoids the need for iterative methods for orthogonal diagonalization of the matrix. In this paper, we adopt a similar approach to simplify the partial differential equation (14).

To begin, we apply a logarithmic transformation to the variables in the PDE, which enables us to eliminate the mixed derivative terms and reduce the problem to a more manageable form. This transformation technique has been successfully used in previous studies and provides a straightforward way to address the computational challenges posed by mixed derivative terms.

$$x_1=\frac{log(V/K)}{{\sigma }_{\mathrm{V}}},x_2=\frac{log(S/K)}{{\sigma }_{\mathrm{S}}}.$$

(15)

Since the partial differential Equation (14) involves only two spatial variables, S and V, the correlation matrix and its $LD{L}^T$ decomposition can be written in a simplified form

$$R=\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)=\left(\begin{array}{cc}1& 0\\ \rho & 1\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& 1-{\rho }^2\end{array}\right)\left(\begin{array}{cc}1& \rho \\ 0& 1\end{array}\right)=LD{L}^T.$$

(16)

Then, the transformation matrix C is found as

$$C=L^{-1}={\left(\begin{array}{cc}1& 0\\ \rho & 1\end{array}\right)}^{-1}=\left(\begin{array}{cc}1& 0\\ -\rho & 1\end{array}\right).$$

(17)

The transformation $Y=CX$ is used to remove the mixed derivative term in the PDE. Specifically, if we define $Y={(y_1,y_2)}^T$ and $X={(x_1,x_2)}^T$, then we have

$$y_1=\frac{1}{{\sigma }_{\mathrm{V}}}log(V/K),y_2=\frac{1}{{\sigma }_{\mathrm{S}}}log(S/K)-\frac{\rho }{{\sigma }_{\mathrm{V}}}log(V/K).$$

(18)

The inverse transformation is

$$V=K{e}^{{\sigma }_{\mathrm{V}}{y}_1},S=K{e}^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}.$$

(19)

This transformation simplifies the PDE by eliminating the mixed derivative term, resulting in a new PDE for Y. Specifically, if we substitute Y into the original PDE (14), and by denoting $U(\tau ,y_1,y_2)=\frac{1}{K}P(\tau ,S;V)$, we obtain

$$\frac{\partial U}{\partial \tau }=\frac{1}{2}\sum _{i=1}^2\left(d_i\right)\frac{{\partial }^{2}U}{\partial y_i^2}+\sum _{i=1}^2\left(\sum _{j=1}^2\frac{(r-q_j-{\sigma }_j^2/2)c_{ij}}{{\sigma }_j}\right)\frac{\partial U}{\partial y_i}-rU+f\left(U\right),$$

(20)

where $d_i$ are the diagonal elements of matrix D in (16), $c_{ij}$ are the components of matrix C defined in (17), and $f\left(U\right)$ is a penalty term

$$f\left(U\right)=\mu {\left(U(0,y_1,y_2)-U(\tau ,y_1,y_2)\right)}^{+},\mu \gg 1.$$

(21)

Hence, the transformed PDE finally takes the form

$$\begin{array}{cc}\hfill \frac{\partial U}{\partial \tau }=& \frac{1}{2}\frac{{\partial }^{2}U}{\partial y_1^2}+\frac{(1-{\rho }^2)}{2}\frac{{\partial }^{2}U}{\partial y_1^2}+\frac{r-{\sigma }_{\mathrm{V}}^2/2}{{\sigma }_{\mathrm{V}}}\frac{\partial U}{\partial y_1}\hfill \\ \hfill & +\left(\frac{r-q-{\sigma }_{\mathrm{S}}^2/2}{{\sigma }_{\mathrm{S}}}-\rho \frac{r-{\sigma }_{\mathrm{V}}^2/2}{{\sigma }_{\mathrm{V}}}\right)\frac{\partial U}{\partial y_2}-rU+f\left(U\right).\hfill \end{array}$$

(22)

The initial conditions (4) for V can be transformed to the new variables Y using the same transformation, resulting in the transformed initial condition for U:

$$U(0,y_1,y_2)=\left\{\begin{array}{cc}{(e^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}-1)}^{+},\hfill & \mathrm{if}{e}^{{\sigma }_{\mathrm{V}}{y}_1}\ge \frac{D^{*}}{K}+1-e^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}\hfill \\ \frac{(1-\alpha )e^{{\sigma }_{\mathrm{V}}{y}_1}{(e^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}-1)}^{+}}{\frac{D^{*}}{K}+1-e^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}},\hfill & \mathrm{if}{e}^{{\sigma }_{\mathrm{V}}{y}_1}<\frac{D^{*}}{K}+1-e^{{\sigma }_{\mathrm{S}}(\rho y_1+y_2)}\hfill \end{array}\right..$$

(23)

To solve the transformed PDE problem, we use the ETD technique [14], which consists of two steps:

First, we apply a semi-discretization scheme to obtain a system of ordinary differential equations (ODEs) in time. In our case, we use the second-order central difference formula for space discretization, which results in a system of nonlinear ODEs.

Second, we use exponential time integration to solve the system of ODEs. The exponential time integrator is based on a splitting of the nonlinear part of the ODEs and the linear part, which can be solved exactly. This allows for an efficient and accurate numerical solution of the transformed PDE.

3.2. Semi-Discretization

Under the coordinate transformation given by (18), the solution domain $V\times S=[0,\infty )\times [0,\infty )$ is transformed to the infinite domain $[y_1,y_2]={\mathbb{R}}^2$. Therefore, the transformed PDE (22) needs to be satisfied for any fixed pair $(y_1,y_2)\in {\mathbb{R}}^2$. To compute a numerical solution, we consider a truncated finite domain $\mathsf{\Omega }=[y_{1_{min}},y_{1_{max}}]\times [y_{2_{min}},y_{2_{max}}]$ and assume that (22) is fulfilled at the boundaries of $\mathsf{\Omega }$, which can then be used as the boundary conditions. Note that the rectangular domain $\mathsf{\Omega }$ after applying the inverse transformation (19) becomes a non-rectangular one. However, it is always possible to choose $\mathsf{\Omega }$ to represent the zone of interest.

Our goal is to construct a numerical solution for the transformed PDE problem (22) with initial conditions (23) on the truncated domain $\mathsf{\Omega }$. The ETD method is based on matrix exponentials for an exact solution of a system of ordinary differential equations (ODE). Therefore, as a preliminary step, we need to perform spatial semi-discretization.

For this purpose, we introduce a uniform mesh with step sizes in each spatial dimension given by

$$h_1=\frac{y_{1_{max}}-y_{1_{min}}}{N_1-1},h_2=\frac{y_{2_{max}}-y_{2_{min}}}{N_2-1},$$

(24)

where $N_1$, $N_2$ is the number of computational nodes in $y_1$ and $y_2$, respectively. Then, the spatial nodes are

$$y_i^j=y_{i_{min}}+j{h}_i,j=0,\dots ,N_i-1,i=1,2.$$

(25)

Let us define an approximated solution at each spatial node by $u_{i,j}\left(\tau \right)\approx U(\tau ,y_1^i,y_2^j)$. Then, for interior nodes ($j\ne 0$, $j\ne N_i-1$, $i=1,2$), the differential operator in $y_1$-dimension are discretized by using the central finite difference (FD) schemes of the second order

$$\begin{array}{ccc}\hfill \frac{\partial U}{\partial y_1}(\tau ,y_1^i,y_2^j)& =& \frac{u_{i+1,j}\left(\tau \right)-u_{i-1,j}\left(\tau \right)}{2h_1}+O\left(h_1^2\right),\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill \frac{{\partial }^{2}U}{\partial y_1^2}(\tau ,y_1^i,y_2^j)& =& \frac{u_{i+1,j}\left(\tau \right)-2u_{i,j}\left(\tau \right)+u_{i-1,j}\left(\tau \right)}{h_1^2}+O\left(h_1^2\right).\hfill \end{array}$$

(27)

At the boundary nodes, equation (22) holds, thus, the discretization of the spatial derivatives is established by using one-sided FD of the second order, $j=1,\dots ,N_1-2$,

$$\begin{array}{ccc}\hfill \frac{\partial U}{\partial y_1}(\tau ,y_1^0,y_2^j)& =& \frac{1}{2h_1}\left(-3u_{0,j}\left(\tau \right)+4u_{1,j}\left(\tau \right)-u_{2,j}\left(\tau \right)\right)+O\left(h_1^2\right),\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \frac{\partial U}{\partial y_1}(\tau ,y_1^{N_1-1},y_2^j)& =& \frac{1}{2h_1}\left(3u_{N_1-1,j}\left(\tau \right)-4u_{N_1-2,j}\left(\tau \right)+u_{N_1-3,j}\left(\tau \right)\right)+O\left(h_1^2\right),\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \frac{{\partial }^{2}U}{\partial y_1^2}(\tau ,y_1^0,y_2^j)& =& \frac{1}{h_1^2}\left(2u_{0,j}\left(\tau \right)-5u_{1,j}\left(\tau \right)+4u_{2,j}\left(\tau \right)-u_{3,j}\left(\tau \right)\right)+O\left(h_1^2\right),\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \frac{{\partial }^{2}U}{\partial y_1^2}(\tau ,y_1^{N_1-1},y_2^j)& =& \frac{1}{h_1^2}\left(2u_{N_1-1,j}\left(\tau \right)-5u_{N_1-2,j}\left(\tau \right)+4u_{N_1-3,j}\left(\tau \right)-u_{N_1-4,j}\left(\tau \right)\right)\hfill \\ & & +O\left(h_1^2\right),\hfill \end{array}$$

(31)

Analogously, FD approximations can be obtained for derivatives with respect to $y_2$.

Substituting the spatial derivatives in (22) by the finite-difference approximations (26)–(28) at each spatial node $(y_1^i,y_2^j)$, $i=0,\dots ,N_1-1$, $j=0,\dots ,N_2-1$, one gets a system of $N=N_1\times N_2$ nonlinear ODEs which can be written in the following matrix form

$$\frac{d\mathbf{u}}{d\tau }\left(\tau \right)=A\mathbf{u}\left(\tau \right)+f\left(\mathbf{u}\left(\tau \right)\right),$$

(32)

where $\mathbf{u}\left(\tau \right)$ is a vector obtained by reshaping the matrix ${\left\{u_{i,j}\right\}}_{0\le i\le N_1-1,0\le j\le N_2-1}$ column-wise such that

$$u_{i,j}\left(\tau \right)={\mathbf{u}}_m\left(\tau \right),m=i+j{N}_1.$$

(33)

In (32), $A={\left\{a_{ij}\right\}}_{0\le i,j<N}$ is sparse and contains constant coefficients of the FD approximations, see [17] for details, $f\left(\mathbf{u}\right)$ is the penalty term

$$f\left(\mathbf{u}\left(\tau \right)\right)=\mu {\left({(1-\mathbf{s})}^{+}-\mathbf{u}\left(\tau \right)\right)}^{+},\mu \gg 1,$$

(34)

where $\mathbf{s}$ is a vector of the following components

$${\mathbf{s}}_m=exp\left({\sigma }_{\mathrm{S}}({\rho }_1^i+y_2^j)\right),m=i{N}_2+j,0\le i\le N_1-1,0\le j\le N_2-1.$$

(35)

Let us consider the ETD method [14] for the system (32). To achieve temporal discretization, the time-step $k=\frac{T}{N_{\tau }}$ is chosen, and the solution at each moment ${\tau }^n=nk$ is obtained.

$$\mathbf{u}\left({\tau }^{n+1}\right)=e^{Ak}\mathbf{u}\left({\tau }^n\right)+{\int }_0^k{e}^{As}f\left(\mathbf{u}({\tau }^{n+1}-s)\right)ds.$$

(36)

According to [14], Section 2.1, the integral equation (36) is equivalent to the system (32) in some given interval $[{\tau }^n,{\tau }^{n+1}]$.

Approximating unknown values $\mathbf{u}({\tau }^{n+1}-s)$ by known $\mathbf{u}\left({\tau }^n\right)$ for $s\in [0,k]$, which has the local truncation error $O\left(k^2\right)$[14,17], and applying the Simpsons’s quadrature rule for the integral term, one gets the explicit formula for $\mathbf{u}\left({\tau }^{n+1}\right)$

$$\mathbf{u}\left({\tau }^{n+1}\right)=e^{Ak}\mathbf{u}\left({\tau }^n\right)+\frac{k}{6}\left(e^{Ak}+4e^{Ak/2}+I\right)f\left(\mathbf{u}\left({\tau }^n\right)\right)+O\left(k^2\right),$$

(37)

where I is the identical $N\times N$ matrix. Note that approximation of the integral term by the the Simpson’s quadrature rule avoids the computation of matrix inverse, that is analytically impossible for singular matrix A and very challenging numerical task for its approximation.

Explicit formula (37) finds the solution $\mathbf{u}\left({\tau }^{n+1}\right)$, $n=0,\dots ,N_{\tau }-1$, level-by-level starting from the initial value $\mathbf{u}\left(0\right)$ given by (23).

3.3. Default case solution

The problem of pricing vulnerable American options is highly challenging due to its multidimensionality, the presence of mixed derivative terms, early exercise options, and the possibility of default. In the event of a default prior to maturity, the option price is calculated using equation (11), while the non-vulnerable American option price is determined by solving the free-boundary PDE problem (6)–(10).

Several numerical methods have been proposed in the literature for American option pricing, including finite difference methods (FDM) [18,19], analytical approximations [20,21], and mesh-free methods [22,23], among others. In this paper, we adopt the numerical algorithm proposed in [24], which is based on the front-fixing transformation combined with explicit FDM. The algorithm computes the solution at all time-levels using a pre-determined time-step k and an appropriate spatial step-size h to ensure the stability of the numerical solution. Since the explicit FDM is used and no matrix computations are required, the algorithm is both robust and fast.

After computing the value of the non-vulnerable American option, it can be incorporated into the proposed algorithm through interpolation. In the case of default, when $S\le K$, the vulnerable option price is computed using (11), which also requires the values of the corresponding non-vulnerable option. Therefore, for every node $(y_1^i,,y_2^j)$ at time ${\tau }^n$, the solution can be found by using the following algorithm:

• Compute corresponding values S and V by using inverse transformation (19).
• If $V\le D^{*}+P_v({\tau }^n,S)$ (default occurs prior maturity) and $S\le K$, then

  $$u_{i,j}\left({\tau }^n\right)=(1-\alpha )\frac{V}{D^{*}+P_v(\tau ,S)}\frac{1}{K}{P}_v(\tau ,S),$$
• If $V>D^{*}+P_v({\tau }^n,S)$ (no default), then

  $$u_{i,j}\left({\tau }^n\right)={\mathbf{u}}_m\left({\tau }^n\right),m=i{N}_2+j,$$

  where ${\mathbf{u}}_m\left({\tau }^n\right)$ is calculated by (37).

To evaluate the performance of the proposed numerical algorithm for pricing vulnerable American options subject to default risk, we implemented the algorithm in MATLAB R2022b. In the following section, we present the numerical results obtained for different parameter values and compare them with benchmark results from the literature.

4. Numerical results

In this section, we present the results of our simulations and compare them with existing methods from the literature. We also provide an analysis of the numerical errors and convergence rates of the proposed algorithm.

The example is based on the vulnerable American option pricing problem given in [6]. We consider an option with maturity $T=2$ years, strike price $K=200$, and volatility of the underlying asset ${\sigma }_{\mathrm{S}}=0.2$. The dividend yield is assumed to be zero ($q=0$) and the interest-free rate $r=0.05$. The option writer is assumed to be highly levered, $D^{*}=900$, $V=1000$, the volatility of the writer’s assets is ${\sigma }_{\mathrm{V}}=0.2$, the deadweight cost is $\alpha =0.25$. We assume that underlying asset S and the return on the writer’s assets V are correlated with some $\rho$.

Detailed numerical solutions for scenarios when $\rho =0$ and $\rho =0.4$ have been plotted and compared with non-vulnerable options in Figure 1 and Figure 2, respectively. It’s worth noting that in cases where $\rho =0$, the equation is devoid of a mixed derivative term and only the transformation as per Equation (15) is implemented, given that the correlation matrix R is identical. As for scenarios where $\rho =0.4$, both the payoff and the numerical solution corresponding to $\tau =T$ are graphically represented in Figure 3.

Table 1 provides a detailed comparison of option prices for varying values of $\rho$, in accordance with the solution proposed by [6]. The computations cover a range denoted by $\mathsf{\Omega }=[-2,,10]\times [-10,,1]$, strategically selected to represent the most significant area of interest. For the discretization parameters, $N_1$ and $N_2$ are both set at 100, and the temporal step size is determined as $k=0.01·min{(h_1,h_2)}^2$. Comparing our findings with those of [6], a slight discrepancy becomes noticeable. This deviation may be attributable to the linear interpolation error that arises when determining the price at the specifically designated point S.

In all considered scenarios, it’s evident that the price of a vulnerable option is invariably lower than that of its non-vulnerable counterpart. This observation is largely ascribed to the default risk associated with the option writer, as discussed in [6].

Numerical methods require two fundamental qualities: stability and convergence. Stability, in the context of numerical methods, refers to the method’s capacity to limit errors during computation. On the other hand, convergence refers to the ability of the numerical method to approach the exact solution as the computation progresses or the step size decreases.

However, due to the intricate nature of numerical algorithms, it is often tedious and even not feasible to analytically study these properties. The complexity inherent in these algorithms often precludes a straightforward mathematical analysis of their stability and convergence characteristics.

To gain a deeper understanding of these characteristics within the context of the proposed method, we conduct a numerical investigation of both stability and convergence. This approach allows us to assess the performance of the algorithm under different conditions and verify whether it maintains the desired qualities of stability and convergence, providing valuable insights into the method’s effectiveness and reliability.

4.1. Numerical stability

First, let’s examine the aspect of numerical stability. For this purpose, we revisit the previously discussed example, fixing the spatial discretization steps at $h_1$ and $h_2$. We aim to identify the parabolic mesh ratio $\delta$, such that $k=\delta min\{h_1^2,h_2^2\}$ ensures a stable solution. It is evident that this ratio is influenced by the specific parameters of the problem. However, it was also observed that it depends on the penalty parameter $\mu$ that is related to the nonlinearity of the PDE. Therefore, in our study, we will also vary this parameter to explore its influence on stability.

For the previously discussed example, wherein $[N_1,N_2]=[50,50]$ are held constant and the values of $\mu$ are varied, the corresponding parabolic mesh ratios, $\delta$, are tabulated in Table 2. Notably, up to $\mu ={10}^3$, the parabolic mesh ratio remains unchanged. This implies that for such values of $\mu$, the stability is primarily determined by the problem parameters, and k should be selected accordingly to ensure this stability condition. However, as $\mu$ becomes significantly large, the influence of the penalty term starts to dominate over the problem parameters.

4.2. Numerical convergence

Drawing upon the ideas presented in [24], it is straightforward to demonstrate that the proposed scheme is consistent with the PDE problem, with a local truncation error of $\mathcal{O}(h^2,k)$. In this subsection, our aim is to numerically study the convergence of the proposed method. We revisit the previous example with all but one step size held constant, which allows us to analyze the convergence with respect to the chosen variable.

We calculate the approximated option value for integer values of S in the interval $[157,170]$, replicating the process used in the previous example. The option price computed using the binomial tree approach, as described in [6], is used as our reference value. For each simulation - that is, for each set of fixed step sizes - we compute the maximum relative error as follows

$$Error(N_1,N_2,k)={∥\frac{P_{\mathrm{ref}}(S,0)-P_{\mathrm{num}}(S,0)}{P_{\mathrm{ref}}(S,0)}∥}_{\infty }=\underset{S=157,\dots ,170}{max}\frac{|P_{\mathrm{ref}}(S,0)-P_{\mathrm{num}}(S,0)|}{|P_{\mathrm{ref}}(S,0)|},$$

(38)

where $P_{\mathrm{ref}}$ is the reference value and $P_{\mathrm{num}}$ is the computed numerical solution.

The errors for various $N_1$ and fixed $N_2=100$ and $k={10}^{-4}$, as well as corresponding errors for various $N_2$ (fixed $N_1=100$ and $k={10}^{-4}$) are plotted in Figure 4. As expected, increasing number of steps leads to more precise solution.

Similar plot for errors with respect to time-step k for fixed $N_1,N_2=100$ is presented in Figure 5. Note that apart from $k=1.6·{10}^{-3}$ the behaviour of error grows exponentially due to the broken stability condition. Obviously, if k is not small enough to guarantee the stability of numerical solution, the relative error increases.

Let us define the log-ratios of errors as in [26]:

$${\gamma }_{N_1}={log}_2\frac{Error(N_1,N_2,k)}{Error(2N_1,N_2,k)},{\gamma }_{N_2}={log}_2\frac{Error(N_1,N_2,k)}{Error(N_1,2N_2,k)}$$

(39)

$${\gamma }_k={log}_2\frac{Error(N_1,N_2,2k)}{Error(N_1,N_2,k)}.$$

(40)

These ratios are calculated for all possible pairs, and the mean values are reported in Table 3.

5. Conclusions

The present study addresses the challenging problem of pricing vulnerable American options through the development of a novel numerical algorithm. We focus on the two-dimensional Black-Scholes equation, which presents several computational difficulties that must be overcome to obtain an accurate and stable solution.

One of the main challenges is the presence of a cross-derivative term in the PDE, which is caused by the correlation between the asset price and the value of the writer’s assets. This term requires a larger computational stencil and may lead to instabilities. To address this issue, we apply a cross-derivative removing transformation, which simplifies the numerical scheme and improves its stability.

Another difficulty arises from the early-exercise opportunity of American-style options, which requires the consideration of a linear complementarity problem. We propose a penalty method to convert this problem into a non-linear PDE, which can be more easily solved numerically.

The choice of appropriate boundary conditions is also crucial for obtaining an accurate solution. We take advantage of the cross-derivative removing transformation to assume that the PDE holds at the boundaries of the computational domain. As a result, we substitute fixed boundary conditions with one-sided finite difference schemes, which provide greater flexibility and versatility in the algorithm.

To solve the resulting numerical problem, we utilize the ETD method, which is a fast and easy-to-implement approach that allows us to compute the option value for a range of reasonable values of the writer’s assets. We conduct numerical experiments to compare our proposed method with other approaches in the literature and to establish its convergence rate and study the stability. Our results demonstrate the effectiveness of the proposed algorithm in accurately and efficiently pricing vulnerable American options.

In conclusion, the proposed numerical algorithm provides a valuable tool for risk management in financial markets, offering a reliable and efficient solution for pricing vulnerable American options.