Next Article in Journal
Reliability of Partitioning Metric Space Data
Previous Article in Journal
Local Second Order Sobolev Regularity for p-Laplacian Equation in Semi-Simple Lie Group
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An ETD Method for Vulnerable American Options

1
Instituto de Matemática Multidisciplinar, Universitat Politècnica de València, Camino de Vera, s/n, 46022 Valencia, Spain
2
Departamento de Matemática Aplicada y Ciencias de la Computación, Universidad de Cantabria, Avenida de los Castros, s/n, 39005 Santander, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(4), 602; https://doi.org/10.3390/math12040602
Submission received: 1 February 2024 / Revised: 14 February 2024 / Accepted: 15 February 2024 / Published: 17 February 2024
(This article belongs to the Section Difference and Differential Equations)

Abstract

:
This paper introduces the exponential time differencing (ETD) technique as a numerical method to efficiently solve vulnerable American options pricing. We address several challenges, including removing cross-derivative terms through appropriate transformations, treating early-exercise opportunities using the penalty method, and substituting fixed boundary conditions with corresponding one-sided finite differences. The proposed method is shown to be both accurate and efficient through numerical experiments, which also compare the results with existing methods and analyze the numerical stability and convergence rate.
MSC:
65M06; 65M20; 65M22; 91G20

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 was extended to accommodate default before option maturity, while Ref. [5] explored 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 vulnerable American 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 Modeling

Let us consider an American option written on the asset S, which follows the risk-neutral process
d S S = ( r q ) d t + σ S d Z S ,
where r is the risk-free rate, q is the dividend yield, σ S is the constant volatility of the underlying asset S, and Z S is a standard Wiener process.
If the option is traded in the over-the-counter market, and there is no guarantee that participants will fulfill their obligations, the option writers become vulnerable to a counter-party credit risk, and the corresponding options are known as vulnerable.
Johnson and Stulz [1] were pioneers in studying vulnerable options. They assumed that default may occur 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 the underlying asset of the option, can lead to financial distress. Later on, in [6], Klein and Yang derived the PDE formulation for the 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.
d V V = r d t + σ V d Z V ,
where σ V is the instantaneous standard deviation of the return on the writer’s assets, and Z V is a standard Wiener process. The instantaneous correlation coefficient between d Z S and d Z V is ρ .
Let τ = T t be the time to maturity; then, the vulnerable American put option price P ( τ , S ; V ) is the solution of the following nonlinear PDE
P τ = σ S 2 2 S 2 2 P S 2 + σ V 2 2 V 2 2 P V 2 + ρ σ S σ V 2 P S V + ( r q ) S P S + r V P V r P ,
subject to the following initial conditions
P ( 0 , S ; V ) = ( S K ) + 1 ( V D * + K S ) + ( 1 α ) V D * + K S 1 ( V < D * + K S ) ,
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 D * + ( K S ) + , the entire claim is paid out. Parameter α 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 1 A is an indicator function.
The boundary condition when S for the PDE put problem (3) is established as follows.
lim S P ( τ ; S , V ) = 0 , τ > 0 .
Let us denote the value of the corresponding non-vulnerable American option by P v ( τ , S ) , which can be calculated as a solution of the Black–Scholes problem [15],
P v τ = σ S 2 2 S 2 2 P v S 2 + ( r q ) S P v S r P v , S > S f ( τ ) ,
where S f ( τ ) is the unknown early exercise boundary, subject to the following initial and boundary conditions.
P v ( 0 , S ) = ( K S ) + , S f ( 0 ) = K ,
P v ( τ , S f ( τ ) ) = K S f ( τ ) , P v S ( τ , S f ( τ ) ) = 1 ,
P v ( τ , S ) = K S , 0 S < S f ( τ ) ,
lim S P v ( τ , S ) = 0 .
If default occurs prior to maturity, the option price is calculated as follows
P ( τ , S ; V ) = ( 1 α ) V D * + P v ( τ , S ) P v ( τ , S ) , S K , V D * + P v ( τ , S ) ,
where P v ( τ , S ) is the value of a vanilla (non-vulnerable) American option defined by the PDE problem (6). Expression (11) defines the amount the holder will receive if the value of the option writer’s assets does not exceed the variable default boundary prior to maturity.
If default does not occur prior to maturity ( V > D * + P v ( τ , S ) ), then the option price is a solution of PDE (3) with the conditions imposed on the optimal stopping boundary S f ( τ ) [6]:
P ( τ , S f ( τ ) , V ) = K S f ( τ ) , P S ( τ , S f ( τ ) , V ) = 1 .
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 ( P ) = μ ( ( K S ) + P ( τ , S ; V ) ) + ,
where μ 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 μ = 0 , the PDE problem (14) reduces to the vulnerable European option.
Finally, the PDE problem to be solved is
P τ = σ S 2 2 S 2 2 P S 2 + σ V 2 2 V 2 2 P V 2 + ρ σ S σ V 2 P S V + ( r q ) S P S + r V P V r P + f ( P ) ,
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 to obtain the solution not in one fixed point but in some domain. Therefore, in the present study, the method of exponential time differencing proposed in [14] is employed after applying the mixed derivative terms removing transformation [13] and the 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 L D 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 = log ( V / K ) σ V , x 2 = log ( S / K ) σ S .
Since the partial differential Equation (14) involves only two spatial variables, S and V, the correlation matrix and its L D L T decomposition can be written in a simplified form
R = 1 ρ ρ 1 = 1 0 ρ 1 1 0 0 1 ρ 2 1 ρ 0 1 = L D L T .
Then, the transformation matrix C is found as
C = L 1 = 1 0 ρ 1 1 = 1 0 ρ 1 .
The transformation Y = C X 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 = 1 σ V log ( V / K ) , y 2 = 1 σ S log ( S / K ) ρ σ V log ( V / K ) .
The inverse transformation is
V = K e σ V y 1 , S = K e σ S ( ρ y 1 + y 2 ) .
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), denoting U ( τ , y 1 , y 2 ) = 1 K P ( τ , S ; V ) , we obtain
U τ = 1 2 i = 1 2 ( d i ) 2 U y i 2 + i = 1 2 j = 1 2 ( r q j σ j 2 / 2 ) c i j σ j U y i r U + f ( U ) ,
where d i are the diagonal elements of matrix D in (16), c i j are the components of matrix C defined in (17), and f ( U ) is a penalty term
f ( U ) = μ U ( 0 , y 1 , y 2 ) U ( τ , y 1 , y 2 ) + , μ 1 .
Hence, the transformed PDE finally takes the form
U τ = 1 2 2 U y 1 2 + ( 1 ρ 2 ) 2 2 U y 1 2 + r σ V 2 / 2 σ V U y 1 + r q σ S 2 / 2 σ S ρ r σ V 2 / 2 σ V U y 2 r U + f ( U ) .
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 ) = ( e σ S ( ρ y 1 + y 2 ) 1 ) + , if e σ V y 1 D * K + 1 e σ S ( ρ y 1 + y 2 ) ( 1 α ) e σ V y 1 ( e σ S ( ρ y 1 + y 2 ) 1 ) + D * K + 1 e σ S ( ρ y 1 + y 2 ) , if e σ V y 1 < D * K + 1 e σ S ( ρ y 1 + y 2 ) .
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 × S = [ 0 , ) × [ 0 , ) is transformed to the infinite domain [ y 1 , y 2 ] = R 2 . Therefore, the transformed PDE (22) needs to be satisfied for any fixed pair ( y 1 , y 2 ) R 2 . To compute a numerical solution, we consider a truncated finite domain Ω = [ y 1 min , y 1 max ] × [ y 2 min , y 2 max ] and assume that (22) is fulfilled at the boundaries of Ω , which can then be used as the boundary conditions. Note that the rectangular domain Ω after applying the inverse transformation (19) becomes a non-rectangular one. However, it is always possible to choose Ω 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 Ω . 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 = y 1 max y 1 min N 1 1 , h 2 = y 2 max y 2 min N 2 1 ,
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 , , N i 1 , i = 1 , 2 .
Let us define an approximated solution at each spatial node by u i , j ( τ ) U ( τ , y 1 i , y 2 j ) . Then, for interior nodes ( j 0 , j N i 1 , i = 1 , 2 ), the differential operators in the y 1 -dimension are discretized using the central finite difference (FD) schemes of the second order.
U y 1 ( τ , y 1 i , y 2 j ) = u i + 1 , j ( τ ) u i 1 , j ( τ ) 2 h 1 + O ( h 1 2 ) ,
2 U y 1 2 ( τ , y 1 i , y 2 j ) = u i + 1 , j ( τ ) 2 u i , j ( τ ) + u i 1 , j ( τ ) h 1 2 + O ( h 1 2 ) .
At the boundary nodes, Equation (22) holds; thus, the discretization of the spatial derivatives is established by using a one-sided FD of the second order, j = 1 , , N 1 2 ,
U y 1 ( τ , y 1 0 , y 2 j ) = 1 2 h 1 3 u 0 , j ( τ ) + 4 u 1 , j ( τ ) u 2 , j ( τ ) + O ( h 1 2 ) ,
U y 1 ( τ , y 1 N 1 1 , y 2 j ) = 1 2 h 1 3 u N 1 1 , j ( τ ) 4 u N 1 2 , j ( τ ) + u N 1 3 , j ( τ ) + O ( h 1 2 ) ,
2 U y 1 2 ( τ , y 1 0 , y 2 j ) = 1 h 1 2 2 u 0 , j ( τ ) 5 u 1 , j ( τ ) + 4 u 2 , j ( τ ) u 3 , j ( τ ) + O ( h 1 2 ) , 2 U y 1 2 ( τ , y 1 N 1 1 , y 2 j ) = 1 h 1 2 2 u N 1 1 , j ( τ ) 5 u N 1 2 , j ( τ ) + 4 u N 1 3 , j ( τ ) u N 1 4 , j ( τ )
+ O ( h 1 2 ) ,
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)–(31) at each spatial node ( y 1 i , y 2 j ) , i = 0 , , N 1 1 , j = 0 , , N 2 1 , one obtains a system of N = N 1 × N 2 nonlinear ODEs, which can be written in the following matrix form
d u d τ ( τ ) = A u ( τ ) + f u ( τ ) ,
where u ( τ ) is a vector obtained by reshaping the matrix { u i , j } 0 i N 1 1 , 0 j N 2 1 column-wise such that
u i , j ( τ ) = u m ( τ ) , m = i + j N 1 .
In (32), A = { a i j } 0 i , j < N is sparse and contains constant coefficients of the FD approximations; see [17] for details. f ( u ) is the penalty term
f ( u ( τ ) ) = μ ( 1 s ) + u ( τ ) + , μ 1 ,
where s is a vector of the following components.
s m = exp σ S ( ρ 1 i + y 2 j ) , m = i N 2 + j , 0 i N 1 1 , 0 j N 2 1 .
Let us consider the ETD method [14] for the system (32). To achieve temporal discretization, the time step k = T N τ is chosen, and the solution at each moment τ n = n k is obtained.
u ( τ n + 1 ) = e A k u ( τ n ) + 0 k e A s f ( u ( τ n + 1 s ) ) d s .
According to [14], Section 2.1, the integral Equation (36) is equivalent to the system (32) in some given interval [ τ n , τ n + 1 ] .
Approximating unknown values u ( τ n + 1 s ) by known u ( τ n ) for s [ 0 , k ] , which has the local truncation error O ( k 2 ) [14,17], and applying the Simpsons’s quadrature rule for the integral term, one obtains the explicit formula for u ( τ n + 1 )
u ( τ n + 1 ) = e A k u ( τ n ) + k 6 e A k + 4 e A k / 2 + I f ( u ( τ n ) ) + O ( k 2 ) ,
where I is the identical N × N matrix. Note that the approximation of the integral term using the the Simpson’s quadrature rule avoids the computation of the matrix inverse, which is analytically impossible for a singular matrix A and a very challenging numerical task for its approximation. Moreover, the second-order discretization scheme provides an optimal balance between computational efficiency and accuracy, combining simplicity in implementation with robust performance. Although higher-order schemes can achieve greater precision, this advantage often comes at the cost of increased computational resources.
The explicit formula (37) finds the solution u ( τ n + 1 ) , n = 0 , , N τ 1 , level-by-level, starting from the initial value u ( 0 ) 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 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 τ n , the solution can be found using the following algorithm.
  • Compute the corresponding values S and V using inverse transformation (19).
  • If V D * + P v ( τ n , S ) (default occurs prior maturity) and S K , then
    u i , j ( τ n ) = ( 1 α ) V D * + P v ( τ , S ) 1 K P v ( τ , S ) ,
  • If V > D * + P v ( τ n , S ) (no default), then
    u i , j ( τ n ) = u m ( τ n ) , m = i N 2 + j ,
    where u m ( τ n ) 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 σ 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 σ V = 0.2 , and the deadweight cost is α = 0.25 . We assume that underlying asset S and the return on the writer’s assets V are correlated with some ρ .
Detailed numerical solutions for scenarios when ρ = 0 and ρ = 0.4 have been plotted and compared with non-vulnerable options in Figure 1 and Figure 2, respectively. It is worth noting that in cases where ρ = 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 ρ = 0.4 , both the payoff and the numerical solution corresponding to τ = T are graphically represented in Figure 3.
Table 1 provides a detailed comparison of option prices for varying values of ρ , in accordance with the solution proposed by [6]. The computations cover a range denoted by Ω = [ 2 , 10 ] × [ 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 the considered scenarios, it is 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. When it comes to ETD schemes, there are not many studies that look into their theoretical stability and convergence. A few studies, like [16,26,27], have made some progress in this area. However, our main interest is in examining these algorithms through numerical studies and experiments. We focus on using numerical methods and experiments to understand the qualitative properties of these algorithms. This approach helps us learn about the stability and convergence of these algorithms in a practical way, filling the gap left by the lack of theoretical studies. Moreover, this methodology allows us to assess the algorithm’s performance under various conditions, scrutinizing its consistency in preserving the essential characteristics of stability and convergence.

4.1. Numerical Stability

First, let us 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 δ , such that k = δ 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 μ 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 μ are varied, the corresponding parabolic mesh ratios, δ , are tabulated in Table 2. Notably, up to μ = 10 3 , the parabolic mesh ratio remains unchanged. This implies that for such values of μ , the stability is primarily determined by the problem parameters, and k should be selected accordingly to ensure this stability condition. However, as μ 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 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.
E r r o r ( N 1 , N 2 , k ) = P ref ( S , 0 ) P num ( S , 0 ) P ref ( S , 0 ) = max S = 157 , , 170 | P ref ( S , 0 ) P num ( S , 0 ) | | P ref ( S , 0 ) | ,
where P ref is the reference value, and P 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, an increasing number of steps leads to a more precise solution.
A 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 behavior of error grows exponentially due to the broken stability condition. Obviously, if k is not small enough to guarantee the stability of the numerical solution, the relative error increases.
Let us define the log-ratios of errors as in [28]:
γ N 1 = log 2 E r r o r ( N 1 , N 2 , k ) E r r o r ( 2 N 1 , N 2 , k ) , γ N 2 = log 2 E r r o r ( N 1 , N 2 , k ) E r r o r ( N 1 , 2 N 2 , k )
γ k = log 2 E r r o r ( N 1 , N 2 , 2 k ) E r r o r ( N 1 , N 2 , k ) .
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 nonlinear 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.

Author Contributions

Conceptualization: R.C., V.N.E. and L.J.; Methodology, R.C., V.N.E. and L.J.; Validation: R.C., V.N.E. and L.J.; Formal analysis: R.C., V.N.E. and L.J.; Investigation: R.C., V.N.E. and L.J.; Writing—review and editing: R.C., V.N.E. and L.J. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially supported by the Spanish Ministry of Economy and Competitiveness MINECO through the project PID2019-107685RB-I00 and by the Spanish State Research Agency (AEI) through the project PDC2022-133115-I00.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
ETDExponential time differencing
FDFinite difference
FDMFinite difference method
PDEPartial differential equation
ODEOrdinary differential equation

References

  1. Johnson, H.; Stulz, R. The Pricing of Options with Default Risk. J. Financ. 1987, 42, 267–280. [Google Scholar] [CrossRef]
  2. Klein, P.G. Pricing Black-Scholes options with correlated credit risk. J. Bank. Financ. 1996, 20, 1211–1229. [Google Scholar] [CrossRef]
  3. Klein, P.; Inglis, M. Pricing vulnerable European options when the option’s payoff can increase the risk of financial distress. J. Bank. Financ. 2001, 25, 993–1012. [Google Scholar] [CrossRef]
  4. Liao, S.L.; Huang, H.H. Pricing Black–Scholes options with correlated interest rate risk and credit risk: An extension. Quant. Financ. 2005, 5, 443–457. [Google Scholar] [CrossRef]
  5. Hung, M.W.; Liu, Y.H. Pricing vulnerable options in incomplete markets. J. Futur. Mark. 2005, 25, 135–170. [Google Scholar] [CrossRef]
  6. Klein, P.; Yang, J. Vulnerable American options. Manag. Financ. 2010, 36, 414–430. [Google Scholar] [CrossRef]
  7. Zhou, R. Back to CVA: The Case of American Option. 2019. Available online: https://ssrn.com/abstract=3189805 (accessed on 14 February 2024).
  8. Chang, L.F.; Hung, M.W. Valuation of vulnerable American options with correlated credit risk. Rev. Deriv. Res. 2006, 9, 137–165. [Google Scholar] [CrossRef]
  9. Wang, G.; Wang, X.; Liu, Z. Pricing vulnerable American put options under jump-diffusion processes. Probab. Eng. Inform. Sci. 2016, 31, 121–138. [Google Scholar] [CrossRef]
  10. Wang, S.; Zhou, Q.; Xiao, W. Pricing vulnerable American put options under jump-diffusion processes when corporate liabilities are random. Commun. Stat. Simul. Comput. 2021, 52, 5462–5482. [Google Scholar] [CrossRef]
  11. Zvan, R.; Forsyth, P.; Vetzal, K. Negative coefficients in two-factor option pricing models. J. Comput. Financ. 2003, 7, 37–73. [Google Scholar] [CrossRef]
  12. Gyulov, T.B.; Koleva, M.N. Penalty method for indifference pricing of American option in a liquidity switching market. Appl. Numer. Math. 2022, 172, 525–545. [Google Scholar] [CrossRef]
  13. Company, R.; Egorova, V.N.; Jódar, L.; Soleymanim, F. A mixed derivative terms removing method in multi-asset option pricing problems. Appl. Math. Lett. 2016, 60, 108–114. [Google Scholar] [CrossRef]
  14. Cox, S.M.; Matthews, P.C. Exponential Time Differencing for Stiff Systems. J. Comput. Phys. 2002, 176, 430–455. [Google Scholar] [CrossRef]
  15. Black, F.; Scholes, M. The Pricing of Options and Corporate Liabilities. J. Political Econ. 1973, 81, 637–654. [Google Scholar] [CrossRef]
  16. Company, R.; Egorova, V.N.; Jódar, L. Conditional full stability of positivity-preserving finite difference scheme for diffusion-advection-reaction models. J. Comput. Appl. Math. 2018, 341, 157–168. [Google Scholar] [CrossRef]
  17. Company, R.; Egorova, V.N.; Jódar, L. A front-fixing ETD numerical method for solving jump-diffusion American option pricing problems. Math. Comput. Simul. 2021, 189, 69–84. [Google Scholar] [CrossRef]
  18. Nielsen, B.F.; Skavhaug, O.; Tvelto, A. Penalty and front-fixing methods for the numerical solution of American option problems. J. Comput. Financ. 2002, 5, 69–97. [Google Scholar] [CrossRef]
  19. Tangman, D.Y.; Gopaul, A.; Bhuruth, M. A fast high-order finite difference algorithm for pricing American options. J. Comput. Appl. Math. 2008, 222, 17–29. [Google Scholar] [CrossRef]
  20. Geske, R.; Johnson, H.E. The American Put Option Valued Analytically. J. Financ. 1984, 39, 1511–1524. [Google Scholar] [CrossRef]
  21. Ju, N. Pricing an American option by approximating its early exercise boundary as a multipiece exponential function. Rev. Financ. Stud. 1998, 11, 627–646. [Google Scholar] [CrossRef]
  22. Hon, Y.C. A Quasi-Radial Basis Functions Method for American Options Pricing. Comput. Math. Appl. 2002, 43, 513–524. [Google Scholar] [CrossRef]
  23. Giribone, P.G.; Ligato, S. Option pricing via radial basis functions: Performance comparison with traditional numerical integration scheme and parameters choice for a reliable pricing. Int. J. Financ. Eng. (IJFE) 2015, 2, 1–30. [Google Scholar] [CrossRef]
  24. Company, R.; Egorova, V.N.; Jódar, L. Solving American option pricing models by the front fixing method: Numerical analysis and computing. Abstr. Appl. Anal. 2014, 2014, 146745. [Google Scholar] [CrossRef]
  25. Hull, J.; White, A. The impact of default risk on the prices of options and other derivative securities. J. Bank. Financ. 1995, 19, 299–322. [Google Scholar] [CrossRef]
  26. Chen, W.; Li, W.; Luo, Z.; Wang, C.; Wang, X. A stabilized second order exponential time differencing multistep method for thin film growth model without slope selection. ESAIM Math. Model. Numer. Anal. 2020, 54, 727–750. [Google Scholar] [CrossRef]
  27. Cheng, K.; Qiao, Z.; Wang, C. A Third Order Exponential Time Differencing Numerical Scheme for No-Slope-Selection Epitaxial Thin Film Model with Energy Stability. J. Sci. Comput. 2019, 81, 154–185. [Google Scholar] [CrossRef]
  28. Company, R.; Egorova, V.N.; Jódar, L.; Peris, J. A front-fixing method for American option pricing on zero-coupon bond under the Hull and White model. Math. Meth. Appl. Sci. 2022, 45, 3334–3344. [Google Scholar] [CrossRef]
Figure 1. Payoff of the vulnerable option (blue line), vulnerable American option price at τ = 0 (orange line) for V = 1000 and ρ = 0.0 , non-vulnerable American option price (dashed black line).
Figure 1. Payoff of the vulnerable option (blue line), vulnerable American option price at τ = 0 (orange line) for V = 1000 and ρ = 0.0 , non-vulnerable American option price (dashed black line).
Mathematics 12 00602 g001
Figure 2. Payoff of the vulnerable option (blue line), vulnerable American option price at τ = 0 (orange line) for V = 1000 and ρ = 0.4 , non-vulnerable American option price (dashed black line).
Figure 2. Payoff of the vulnerable option (blue line), vulnerable American option price at τ = 0 (orange line) for V = 1000 and ρ = 0.4 , non-vulnerable American option price (dashed black line).
Mathematics 12 00602 g002
Figure 3. Vulnerable American option price ( ρ = 0.4 ) at maturity (left) and at τ = T (right).
Figure 3. Vulnerable American option price ( ρ = 0.4 ) at maturity (left) and at τ = T (right).
Mathematics 12 00602 g003
Figure 4. Maximum relative errors with respect to the number of spatial steps N 1 and N 2 .
Figure 4. Maximum relative errors with respect to the number of spatial steps N 1 and N 2 .
Mathematics 12 00602 g004
Figure 5. Maximum relative errors with respect to temporal step size k.
Figure 5. Maximum relative errors with respect to temporal step size k.
Mathematics 12 00602 g005
Table 1. Price of vulnerable and non-vulnerable American options for various values of correlation ρ . Prices in the corresponding left columns are obtained by following the three-dimensional binomial tree approach [25] and published in [6]. Prices in the corresponding right columns for each ρ are calculated by the proposed method.
Table 1. Price of vulnerable and non-vulnerable American options for various values of correlation ρ . Prices in the corresponding left columns are obtained by following the three-dimensional binomial tree approach [25] and published in [6]. Prices in the corresponding right columns for each ρ are calculated by the proposed method.
S ρ = 0.0 ρ = 0.4 ρ = 0.8 Non-Vulnerable
1574342.994342.994343.0043
1584241.994241.994242.0042.00
1594140.994140.9941.0141.0341.08
1604039.994039.9940.0740.0940.17
1613938.993939.0239.1439.1339.27
1623837.993838.0538.2238.2238.37
1633737.0137.0237.0937.3037.3237.49
1643636.0236.1136.1436.3836.4236.65
1653535.0835.1835.2235.4735.5235.83
16634.0334.1434.2634.3134.5634.6335.01
16733.1033.1833.3433.4033.7233.7434.19
16832.1632.3032.4332.5232.8432.8533.40
16931.2331.4231.5231.6631.9631.9732.63
17030.3930.5530.6330.7931.1131.1231.89
Table 2. The approximate parabolic mesh ratio δ for stable numerical solution by the proposed algorithm with respect to the penalty parameter μ .
Table 2. The approximate parabolic mesh ratio δ for stable numerical solution by the proposed algorithm with respect to the penalty parameter μ .
μ 010 10 2 10 3 10 4
[ N 1 , N 2 ] = [ 50 , 50 ] 0.094 0.094 0.094 0.038 0.003
Table 3. Numerical convergence rates in directions y 1 , y 2 and time τ .
Table 3. Numerical convergence rates in directions y 1 , y 2 and time τ .
γ N 1 ¯ 1.9246
γ N 2 ¯ 1.7748
γ k ¯ 1.0765
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Company, R.; Egorova, V.N.; Jódar, L. An ETD Method for Vulnerable American Options. Mathematics 2024, 12, 602. https://doi.org/10.3390/math12040602

AMA Style

Company R, Egorova VN, Jódar L. An ETD Method for Vulnerable American Options. Mathematics. 2024; 12(4):602. https://doi.org/10.3390/math12040602

Chicago/Turabian Style

Company, Rafael, Vera N. Egorova, and Lucas Jódar. 2024. "An ETD Method for Vulnerable American Options" Mathematics 12, no. 4: 602. https://doi.org/10.3390/math12040602

APA Style

Company, R., Egorova, V. N., & Jódar, L. (2024). An ETD Method for Vulnerable American Options. Mathematics, 12(4), 602. https://doi.org/10.3390/math12040602

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop