An ETD Method for Vulnerable American Options

: 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


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

Vulnerable Option Modeling
Let us consider an American option written on the asset S, which follows the riskneutral process dS S = (r 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.
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 dZ S and dZ 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 subject to the following initial conditions where (•) + = max{•, 0}; D * is the fixed default boundary such that a default occurs if the value of the option writer's assets , 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], where S f (τ) is the unknown early exercise boundary, subject to the following initial and boundary conditions.
If default occurs prior to maturity, the option price is calculated as follows 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]: 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 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 subject to the initial conditions (4) and boundary conditions (11).

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

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 LDL 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.
Since the partial differential Equation ( 14) involves only two spatial variables, S and V, the correlation matrix and its LDL T decomposition can be written in a simplified form Then, the transformation matrix C is found as 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 The inverse transformation is 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 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 (U) is a penalty term Hence, the transformed PDE finally takes the form 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: (1−α)e σ V y 1 (e σ S (ρy 1 +y 2 ) −1) , if 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.

Semi-Discretization
Under the coordinate transformation given by (18), the solution domain 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 where N 1 , N 2 is the number of computational nodes in y 1 and y 2 , respectively.Then, the spatial nodes are Let us define an approximated solution at each spatial node by u i,j (τ) ≈ U(τ, y i 1 , y j 2 ).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
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, . . ., 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 i 1 , y j 2 ), 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 where u(τ) is a vector obtained by reshaping the matrix {u i,j } 0≤i≤N 1 −1, 0≤j≤N 2 −1 columnwise such that In (32), A = {a ij } 0≤i,j<N is sparse and contains constant coefficients of the FD approximations; see [17] for details.f (u) is the penalty term where s is a vector of the following components.
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 = nk is obtained.
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 ) 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).

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 i 1 , , y j 2 ) at time τ n , the solution can be found using the following algorithm.

1.
Compute the corresponding values S and V using inverse transformation (19).
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.

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 Figures 1 and 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.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.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.

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 2 1 , 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.

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.
Error(N 1 , N 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]:

Figure 1 .
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 2 .
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).

2 Figure 4 .
Figure 4. Maximum relative errors with respect to the number of spatial steps N 1 and N 2 .

Figure 5 .
Figure 5. Maximum relative errors with respect to temporal step size k.

Table 1 .
Price of vulnerable and non-vulnerable American options for various values of correlation ρ.

Table 2 .
The approximate parabolic mesh ratio δ for stable numerical solution by the proposed algorithm with respect to the penalty parameter µ.