A Fast Preconditioned Penalty Method for American Options Pricing Under Regime-Switching Tempered Fractional Diffusion Models

Abstract

A fast preconditioned penalty method is developed for a system of parabolic linear complementarity problems (LCPs) involving tempered fractional order partial derivatives governing the price of American options whose underlying asset follows a geometry Lévy process with multi-state regime switching. By means of the penalty method, the system of LCPs is approximated with a penalty term by a system of nonlinear tempered fractional partial differential equations (TFPDEs) coupled by a finite-state Markov chain. The system of nonlinear TFPDEs is discretized with the shifted Grünwald approximation by an upwind finite difference scheme which is shown to be unconditionally stable. Semi-smooth Newton’s method is utilized to solve the finite difference scheme as an outer iterative method in which the Jacobi matrix is found to possess Toeplitz-plus-diagonal structure. Consequently, the resulting linear system can be fast solved by the Krylov subspace method as an inner iterative method via fast Fourier transform (FFT). Furthermore, a novel preconditioner is proposed to speed up the convergence rate of the inner Krylov subspace iteration with theoretical analysis. With the above-mentioned preconditioning technique via FFT, under some mild conditions, the operation cost in each Newton’s step can be expected to be $\mathcal{O}(N\mathrm{log}N)$, where N is the size of the coefficient matrix. Numerical examples are given to demonstrate the accuracy and efficiency of our proposed fast preconditioned penalty method.

Appendices

Appendix A: Derivation of the Eq. (4)

According to the no arbitrage condition and Itô’s formula for switching Lévy process, European claim price must satisfy the following system of partial integro-differential equations [13]:

$$\begin{aligned} \frac{\partial V_s(x,t)}{\partial t}+\mathcal {L} V_s(x,t)-r V_s(x,t)=0,\quad \forall s \in \mathcal{S}, \end{aligned}$$

(22)

where $\mathcal {L} V_s(x,t)$ is the generator of the switching Lévy process:

$$\begin{aligned} \mathcal {L} V_s(x,t)=m_s\frac{\partial V_s(x,t)}{\partial x}+\frac{1}{2}\sigma _s^2 \frac{\partial V_s^2(x,t)}{\partial x^2}+\sum _{\tilde{s}\ne s}q_{s \tilde{s}}(V_{\tilde{s}}(x,t)-V_s(x,t))\nonumber \\ +\int _{{\mathbb {R}}\backslash \{0\}}(V_s(x+z,t)-V_s(x,t)-z h(z)\frac{\partial V_s(x,t)}{\partial x})\;W_s(dz). \end{aligned}$$

Applying the Fourier transform to (22), by means of the condition (ii) of the intensity matrix $\mathcal{Q}$ and (3), leads to the following equations $\forall s \in \mathcal{S}$,

$$\begin{aligned} \frac{\partial \mathcal {F}[V_s(x,t)](\zeta ,t)}{\partial t}+(\varPsi _s(-\zeta )+q_{s s}-r-\mathrm{i}(r-v_s)\zeta )\mathcal {F}[V_s(x,t)](\zeta ,t)\nonumber \\ +\sum _{\tilde{s}\ne s}q_{s \tilde{s}}\mathcal {F}[V_{\tilde{s}}(x,t)](\zeta ,t)=0. \end{aligned}$$

(23)

where the exponential Fourier transform is defined by

$$\begin{aligned} \mathcal {F}[V_s(x,t)](\zeta ,t)\equiv \hat{V}_s(\zeta ,t)=\int ^{\infty }_{-\infty } e^{\mathrm{i} \zeta x} V_s(x,t) dx. \end{aligned}$$

From [8], we have

$$\begin{aligned} \varPsi ^\mathrm{KoBoL}_s(\zeta )=\frac{1}{2}{\sigma _s}^{\alpha _s}\big [p_s(\lambda _s-\mathrm{i} \zeta )^{\alpha _s}+\varsigma _s(\lambda _s+\mathrm{i} \zeta )^{\alpha _s}-{\lambda _s}^{\alpha _s}-\mathrm{i} \zeta {\alpha _s} {\lambda _s}^{\alpha _s-1}(\varsigma _s-p_s)\big ].\nonumber \\ \end{aligned}$$

(24)

By inserting specific characteristic exponent (24) into (23), the system of TFPDEs (4) can be derived from the inverse Fourier transform. To be noted that the Fourier transforms of the left-sided and right-sided tempered fractional operators are given by [10, 32, 33]

$$\begin{aligned}&\mathcal {F}\Big [e^{-\lambda _{s} x}{{_{-\infty }}D^{\alpha _s}_x} [e^{\lambda _{s} x} V_s(x,t)]\Big ]=(\lambda _s-\mathrm{i} \zeta )^{\alpha _s}\hat{V}_s(\zeta ,t),\quad \\&\mathcal {F}\Big [e^{\lambda _{s} x}{_{x}}D^{\alpha _s}_\infty [e^{-\lambda _{s} x} V_s(x,t)]\Big ]=(\lambda _s+\mathrm{i} \zeta )^{\alpha _s}\hat{V}_s(\zeta ,t). \end{aligned}$$

Appendix B: Monotonicity of the Newton’s Iteration

From the algorithm (18),

$$\begin{aligned} \left( \mathcal{M}{_s}-\tau \rho \mathcal{\tilde{I}}_s^k\right) \left( (\bar{u}_s^{(m)})^{k+1}-(\bar{u}_s^{(m)})^{k}\right) +\bar{g}((\bar{u}_s^{(m)})^{k})=0. \end{aligned}$$

Then replacing $\bar{g}((\bar{u}_s^{(m)})^{k})$ by its definition (17), we have

$$\begin{aligned} \left( \mathcal{M}{_s}-\tau \rho \mathcal{\tilde{I}}_s^k\right) \left( (\bar{u}_s^{(m)})^{k+1}-(\bar{u}_s^{(m)})^{k}\right) -(\bar{u}_s^{(m)})^{0}+\mathcal{M}{_s}(\bar{u}_s^{(m)})^{k}-\tau \rho {\mathcal{N}}^{(m)}_{s}=0, \end{aligned}$$

where $(\bar{u}_s^{(m)})^{0}=\tau f_{s}+u_{s}^{(m-1)}-\tau \sum _{\tilde{s}\in \mathcal{S}}q_{s \tilde{s}}u^{(m-1)}_{\tilde{s}}$.

Noting that $\mathcal{\tilde{I}}_s^k (\bar{u}_s^{(m)})^{k}-{\mathcal{N}}^{(m)}_{s}=\mathcal{\tilde{I}}_s^k u_s^*$. Finally, the Newton’s method (18) can be simplified to the following algorithm, which gives for all $s \in \mathcal{S}$,

$$\begin{aligned} \left( \mathcal{M}{_s}-\tau \rho \mathcal{\tilde{I}}_s^k \right) (\bar{u}_s^{(m)})^{k+1}= (\bar{u}_s^{(m)})^{0}-\tau \rho \mathcal{\tilde{I}}_s^k u_s^*{.} \end{aligned}$$

(25)

Given that $\mathcal{M}_s$ is an M-matrix, by following the idea from [18], we can prove that the Newton’s iteration converges monotonically. For the monotone property, writing equation (25) for iteration $k\ge 1$ as:

$$\begin{aligned} \left( \mathcal{M}{_s}-\tau \rho \mathcal{\widetilde{I}}_s^{k} \right) (\bar{u}_s^{(m)})^{k}-\tau \rho \left( \mathcal{\widetilde{I}}_s^{k-1}-\mathcal{\widetilde{I}}_s^{k} \right) (\bar{u}_s^{(m)})^{k}= (\bar{u}_s^{(m)})^{0}-\tau \rho \mathcal{\tilde{I}}_s^{k-1} u_s^*. \end{aligned}$$

(26)

Subtracting equation (26) from (25) gives

$$\begin{aligned} \left( \mathcal{M}_s-\tau \rho \mathcal{\tilde{I}}_s^{k} \right) ((\bar{u}_s^{(m)})^{k+1}-(\bar{u}_s^{(m)})^{k})=\tau \rho \left( \mathcal{\widetilde{I}}_s^{k-1}-\mathcal{\widetilde{I}}_s^{k} \right) (u_s^*-(\bar{u}_s^{(m)})^{k}). \end{aligned}$$

Under the both cases: $u_{n,s}^*-(\bar{u}_{n,s}^{(m)})^{k}\ge 0$ and $u_{n,s}^*-(\bar{u}_{n,s}^{(m)})^{k}< 0$ in the pointwise, we have

$$\begin{aligned} \left( \mathcal{M}_s-\tau \rho \mathcal{\tilde{I}}_s^{k} \right) ((\bar{u}_s^{(m)})^{k+1}-(\bar{u}_s^{(m)})^{k})\le 0. \end{aligned}$$

Note that all of the elements of the inverse of an M-matrix are nonnegative. Since the matrix $\mathcal{M}_s-\tau \rho \mathcal{\tilde{I}}_s^{k}$ is an M-matrix, we can see that the Newton’s iteration is monotonically decreasing for the iteration step $k\ge 1$. Similar with the technique used in [18], we can also prove the uniqueness of the penalty iteration given that the matrix $\mathcal{M}_s$ is an M-matrix.