Pricing perpetual timer options under Heston Model by finite difference method: Theory and implementation

1.   Introduction

Timer option is a kind of exotic option whose payoff is the same as a vanilla option while the expiration date is not given in advance, but is determined by the cumulative variance of the underlying asset. When the cumulative variance reaches a preset threshold, the timer option automatically expires. A timer option is called a perpetual timer option if the option will not expire until the threshold is hit. As mentioned by Sawyer ^[1], timer options were first introduced for trading by Société Générale and Investment Banking since 2007. At that time, such an exotic option provided a more flexible risk management tool which allows market participants specifying the volatility of underlying asset by choosing a certain cumulative variance threshold. Li ^[2] demonstrated that a timer option can be cheaper than a vanilla option with the same strike price and the same expected investment horizon.

Long before timer options were traded in financial markets, Neuberger ^[3] proposed the "mileage" option with variant expiration date and Bick ^[4] discussed the replication of such option under a continuous-time model. As mentioned by Carr & Lee ^[5], these works can be seen as pioneer works of timer options. There are three major approaches to pricing timer options. The most commonly used approach is pricing through the Monte Carlo method. Li ^[2] first proposed a Black-Scholes-Merton type formula for pricing timer options and implemented a Monte Carlo simulation under the Heston model. They employed an Euler scheme on the Bessel process with predictor and corrector to improve accuracy. Bernard & Cui ^[6] successfully reduced the timer option pricing problem to a one-dimensional problem where only the dynamics of the volatility needs to be simulated. They implemented their method for Heston model and Hull White model. Cui ^[7] further improved the simulation method. The second approach is to find a closed-form approximation formula through asymptotic analysis technique. Li & Mercurio ^[8] developed a closed-form approximation for timer options under Heston model and the 3/2-model with small volatility of variance. They also provided an analysis framework for general stochastic volatility models. Furthermore, Li & Mercurio ^[9,10] discussed the asymptotic approximation for both perpetual and finite time timer options. In a recent work of Wang et al. ^[11], the author provided a closed-form approximation formula for timer options under a second-order multi-scale stochastic volatility model. Kwok & Zheng ^[12] introduced a numerical Fourier transform algorithm for finite maturity timer options under the 3/2-model. The third approach is through the path-integral method, which is borrowed from quantum field theory. By the path-integral method, Liang et al. ^[13] provided the multivariate integral expression of the timer option price under the Heston model and the 3/2-model. Cui et al. ^[14] used a stochastic time change technique to develop explicit formula for prices of timer options in the Heston model and the 3/2 model. Zhang et al. ^[15] developed a closed-form pricing formula for timer options under the Hull White model using this method.

The goal of this paper is to provide the fourth approach for pricing perpetual timer option under the Heston model, which is the finite difference method, theoretically and empirically. Different from the three approaches introduced before, especially the most widely used Monte Carlo method, the finite difference method is not a single-point calculation method. It can not only directly give the entire pricing surface, but its calculation results can also be more conveniently used for the calculation of the Greek value surface. In the Monte Carlo method, a huge amount of numerical simulation is required to obtain the entire pricing surface, and the simulation results cannot be directly used for the calculation of the Greek values. In addition, our method can give a priori theoretical error estimates, while the Monte Carlo method can only give confidence intervals for simulation results, which is a posteriori estimate based on relative errors. Compared with the finite difference method, although the calculation efficiency of the asymptotic approximation approach is higher, the accuracy of the approximate formula is limited by specific model parameters and cannot be flexibly adjusted according to actual needs. The closed-form formula obtained by the path-integral method can guarantee high calculation accuracy, but it is very complicated to implement and cannot be used to calculate Greeks. In fact, the finite difference method is one of the most commonly used methods for option pricing problems. It has the advantages of easy implementation, flexibility and high efficiency. However, due to the degenerate nature of the perpetual timer option pricing equation, the common finite difference method can hardly work. Li ^[2] illustrated a numerical example for pricing timer option under the Heston model using the ADI scheme. The numerical results showed that the convergence of this method is problematic. Since then, there is almost no relevant literature discussing the application of finite difference method to the timer option pricing problem.

Let's recall the pricing PDE of the perpetual timer call option under the Heston model from Li & Mercurio ^[9]:

where the variables $x, y$ represent the logarithmic price and instantaneous volatility of the underlying asset respectively, the variable $t$ represents the time to maturity of the option, and the constant $K$ represents the strike price of the option. The operator $\mathcal{L}$ is defined as

The coefficients defining $\mathcal{L}$ are all constants and satisfy following conditions:

where $r$ is the risk-free rate, $\theta$ is the long-term mean of volatility, $\kappa$ indicating the speed of mean reversion, $\sigma$ is the volatility of volatility, and $\rho$ is the correlation coefficient between volatility dynamics and underlying price dynamics. For more details about the financial meaning of these coefficients, we recommend Heston ^[16] and Li & Mercurio ^[9].

From (1.2) we can see that the equation shows degeneracy along the boundary at $\{y = 0\}$ since the coefficients of $u_x, u_y$ and $u$ explode there. We also notice that no boundary conditions are specified along $\{y = 0\}$. Such degeneracy indicates that in order to ensure the existence and uniqueness of the solution of (1.1) we must add suitable asymptotic growth conditions. Furthermore, in order to obtain a numerical solution, we also need to add appropriate boundary conditions at the degenerate boundary.

The main contribution of this paper are as follows. Firstly, we provide an asymptotic growth condition which theoretically guarantee the existence and uniqueness of the solution of the pricing problem (1.1), and we discuss the choice of artificial boundary value conditions at the degenerate boundary and infinity. We show that the error introduced by artificial boundary value conditions decays rapidly inside the solution domain under the control of the growth condition. We provide an a priori estimate of this error. Secondly, we developed a finite difference scheme with second-order convergence in spatial directions to solve the pricing problem numerically. The convergence of our scheme follows from a discrete Aleksandrov-Bakelman-Krylov maximum principle given by Kuo & Trudinger ^[17]. Finally, we conduct numerical experiments with our method.

The remainder of this paper is organized as follows. In section 2, we discuss the asymptotic condition and give a priori error estimation of the artificial boundary value error. In section 3, we provide the finite difference scheme and prove the convergence. In section 4 we conduct numerical experiments to show the convergence and accuracy of our method. Section 5 concludes the paper.

2.   Theoretical analysis

We first make some notations and state the asymptotic growth condition in the following subsection.

2.1.   Notations

Let domain $D\subset\mathbb{R}^n$, $I = (a, b]$ is an interval on $\mathbb{R}^1$, $Q = D\times I$, we thus define boundaries of $Q$ as follows:

We also define the the closure and interior of $Q$ respectively as:

For positive parameters $\delta, R, R_1, R_2$, we define domains that will be used in this article as follows.

For all subdomain $D$ of $\Omega$ and constant $\mathcal{T} > 0$, we denote

For all real number $a$, we define

For any domain $\mathcal{D}\subset\Omega\times\mathbb{R}$, we can define a function class $\mathcal{H}(\mathcal{D})$:

where $f(x, y, t):\mathcal{D}\mapsto\mathbb{R}$. We claim

as the asymptotic growth condition.

2.2.   Existence and uniqueness

The asymptotic growth condition (2.3) actually provides a special boundary condition at the degenerate boundary. With such growth condition we can prove a comparison principle in $\Omega_\mathcal{T}$, and the uniqueness follows.

Theorem 2.1 (Existence and Uniqueness). There exists unique solution $u\in\mathcal{H}(\bar{\Omega}_\mathcal{T})$ of the pricing problem (1.1).

Remark 1. The existence for the solution of (1.1) can be obtained by directly using a classical Schauder estimate (here we refer to Lieberman ^[18]) and standard diagonal sub-sequence procedures. And the uniqueness follows directly from the comparison principle in the following Lemma 2.2.

Lemma 2.2 (Comparison Principle). If $u, v\in \mathcal{H}(\bar{\Omega}_\mathcal{T})$ such that:

then $u\le v$ in $\Omega_\mathcal{T}$.

Proof. By the definition of $\mathcal{H}$ in (2.2), we obtain that $(u-v)\in\mathcal{H}(\bar{\Omega}_\mathcal{T})$, and thus there exists a constant $C > 0$ such that

We define constants $n_0, \lambda$ as:

and constants $N_1, N_2, N_3$ as:

We also denote constant $N = \max(N_1, N_2, N_3)$. To get the comparison principle, we define a barrier function $p(x, y, t)$ as follows

By directly calculation, we can check that $P(x, y, t) > 0$ and $\mathcal{L}P(x, y, t) > 0$ in $\Omega_\mathcal{T}$. Furthermore, we have:

Considering that $u, v\in \mathcal{H}(\bar{\Omega}_\mathcal{T})$, for all $\epsilon > 0$ there exists $\delta_0 > 0, R_0 > 0$ such that

as long as $0 < \delta < \delta_0$ and $R\ge R_0$. We also know that

Then by the standard comparison principle for non-degenerate equations, we have

For all $\epsilon > 0$, let $\delta\rightarrow0$ and $R\rightarrow \infty$, we can see that

Then

is achieved by letting $\epsilon\rightarrow0$. ■

Corollary 2.3. The solution $u\in\mathcal{H}(\bar{\Omega}_\mathcal{T})$ of (1.1) satisfies the following condition

Remark 2. Since $\mathcal{L}e^x = 0, \mathcal{L}0 = 0$ and $\mathcal{L}K = ry^{-1}K > 0$, the condition (2.7) follows directly from Lemma 2.2. And this is also in line with financial intuition, that is, the price of the call option should not be less than its final payoff and should not exceed the underlying asset itself.

2.3.   Boundary error estimation

In practice, finding the numerical solution for the pricing problem (1.1) forces us to make a cut-off from the infinity of $y = \infty$, $x = \pm\infty$ and from the degenerate boundary of $y = 0$. This cut-off leads to a bounded domain problem on $(\Omega_{\delta, R_1, R_2})_\mathcal{T}$. We have to provide additional artificial boundary conditions on $\partial_x(\Omega_{\delta, R_1, R_2})_\mathcal{T}$. The most straightforward way is that we directly provide the exact solution of the pricing problem (1.1) on $\partial_x(\Omega_{\delta, R_1, R_2})_\mathcal{T}$ as the boundary value condition, however, we do not know the exact solution there.

In this section we prove that as long as the artificial boundary value satisfies the asymptotic growth condition (2.3), the error brought by the boundary value is convergent within the interior domain, which means we may only need to provide a not too outrageous boundary value to ensure sufficient accuracy inside the solution area.

We first setup the boundary error equation. The bounded domain pricing problem with artificial boundary value condition $\psi_a(x, y, t)$ is shown as follows:

Let $u$ be the solution of the original pricing problem (1.1). We define the error function as

Naturally, $e_b$ solves following problem:

We provide a priori estimate for $e_b$ in the following Theorem 2.4.

Theorem 2.4 (Boundary Error Estimation). Let $e_b$ be the solution of (2.9), where the artificial boundary value $\psi_a$ satisfies the asymptotic growth condition (2.3). We assume there is a constant $C > 0$ such that $|\psi_a|\le C(e^x + 1)$ in $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$ and choose constants $\lambda, n_0, N$ be the same as in Lemma 2.6. Then we have the following estimate:

where

Remark 3. $E_0$, $E_1$ and $E_2$ control the boundary value errors from $y = \delta$, $x = \pm R_1$ and $y = R_2$ respectively. In practice, for all error tolerance level $\epsilon > 0$ at an interior point $(x, y, t)$ of $(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, we can select sufficiently large $R_2$, $R_1$ and a sufficiently small $\delta$ accordingly to ensure $|e_b(x, y, t)|\le \epsilon$.

Corollary 2.5. Let $e_b$ be the solution of (2.9), if $\psi_a$ satisfies the following condition:

then we have a better estimates as follows:

To show the proof of Theorem 2.4, we first discuss the boundary value error from different directions separately. Since both $u$ and $\psi_a$ satisfy the asymptotic growth condition (2.3), we discuss the constant-controlled boundary value error and the exponentially growing boundary value error in the following Lemma 2.6 and Lemma 2.7, respectively.

Lemma 2.6. If $e_b(x, y, t)$ is the solution of (2.9), and we choose constants $\lambda, n_0, N$ be the same as in (2.4) and (2.5), then we have following estimates:

(i) If $|u(x, y, t) - \psi_a(x, y, t)| \le {\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{y = \delta\}}$ on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have

(ii) If $|u(x, y, t) - \psi_a(x, y, t)| \le {\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{x = -R_1\mathit{\text{or}} x = R_1\}}$ on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have

(iii) If $|u(x, y, t) - \psi_a(x, y, t)| \le {\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{y = R_2\}}$ on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have

Lemma 2.7. If the error function $e_b(x, y, t)$ is the solution of (2.9) and we choose constants $\lambda, n_0, N$ be the same as in Lemma 2.6, then we have following estimates:

(i) If $|u(x, y, t) - \psi_a(x, y, t)| \le e^x{\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{y = \delta\}}$ on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have:

(ii) If $|u(x, y, t) - \psi_a(x, y, t)| \le e^x{\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{x = R_1\}}$ holds on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have:

(iii) If $|u(x, y, t) - \psi_a(x, y, t)| \le e^x{\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{x = -R_1\}}$ holds on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have:

(iv) If $|u(x, y, t) - \psi_a(x, y, t)| \le e^x{\bf 1}_{\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}\cap\{y = R_2\}}$ holds on $\partial_p(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, then we have:

Remark 4. The proof of Lemma 2.6 and Lemma 2.7 can be achieved by directly using a standard comparison principle for non-degenerate equation. We omit the details of the proof in this paper.

Proof of Theorem 2.4. From the Theorem 2.1 we know that the solution $u(x, y, t)$ of the pricing problem (1.1) lies between $(e^x - K)^+$ and $e^x$. Since $|\psi_a|\le C(e^x + 1)$, we have

which indicates that the artificial boundary error is controlled by a linear combination of $e^x$ and $1$ at each direction. Then we employ the boundary error estimates from Lemma 2.6 and Lemma 2.7 and find that:

■

Remark 5. Since $u$ lies between $(e^{x} - K)^+$ and $e^x$, it is reasonable to choose $\psi_a$ accordingly, and thus

Then Corollary 2.5 follows from a similar procedure as in the proof of Theorem 2.4.

3.   Finite difference scheme

In this section, we provide a finite difference scheme for solving the bounded domain problem (2.8). Our scheme is quite different from the ADI scheme, which perform poorly, as shown by Li ^[2]. We prove the convergence of our scheme through a discrete maximum principle given by Kuo & Trudinger ^[17].

We define meshing parameters as follows:

For all $i, j, n\in\mathbb{N}$, we set

We denote $E_{\tau, h}$ as a time-space mesh on $(\Omega_{\delta, R_1, R_2})_\mathcal{T}$, where $E_{\tau, h}\subset(\Omega_{\delta, R_1, R_2})_\mathcal{T}$ and consists of meshing points from (3.1), and we also denote $u(x_i, y_j, t_n)$ as $u_{ij}^n$. Then we define the following difference operators:

The finite difference operator $\mathcal{L}^{\Delta}$ corresponding to $\mathcal{L}$ is defined as follows:

Using (3.2), we can write the definition (3.3) of $\mathcal{L}^{\Delta}$ as follows:

where parameters $a_1-a_6$ are defined as:

It can be seen that our finite difference scheme is actually a five-point scheme. Then we consider the finite difference problem corresponding to (2.8):

where $u_{\tau, h}$ is a real functions defined on $E_{\tau, h}$ and $\psi_a$ is the artificial boundary value condition. Here we present an algorithm for solving the finite difference problem (3.6).

Considering the sparsity of the coefficient matrix, memory usage can be reduced by matrix compression techniques. Linear systems can be solved using the Jacobi iteration method.

Remark 6. If $u_{b}$ is the solution of (2.8) we can give the truncation error of $\mathcal{L}^{\Delta}$ by using Taylor expansion as:

which indicates that our scheme has second-order accuracy in the spatial directions.

Theorem 3.1. If $u_{\tau, h}$ is the solution of the finite difference problem (3.6) and $u_b$ is the solution of bounded domain pricing problem (2.8) then $u_{\tau, h}$ converge to $u_b$ uniformly on $E_{\tau, h}$.

Theorem 3.2. If $u(x, y, t)$ is the solution of (1.1), we choose two artificial boundary conditions $\psi_a^{(u)}(x, y, t) = e^x$ and $\psi_a^{(l)}(x, y, t) = (e^x - K)^+$. Denote $u_b^{(u)}$ and $u_b^{(l)}$ as the solution for bounded domain problem (2.8) under $\psi_a = \psi_a^{(u)}$ and $\psi_a = \psi_a^{(l)}$ respectively. We also denote $u_{\tau, h}^{(u)}$ and $u_{\tau, h}^{(l)}$ as solutions for the finite difference problem (3.6) similarly. Define

then we have:

For the proof of Theorem 3.1, we first introduce a discrete Aleksandrov-Bakelman-Krylov maximum principle, which mainly follows from Kuo & Trudinger ^[17], in the following lemma.

Lemma 3.3 (Discrete Maximum Principle). If $u_{\tau, h}$ is a real function defined on $E_{\tau, h}$, $\mathcal{L}^{\Delta}$ defined as (3.3), $\mathcal{L}^{\Delta}u_{\tau, h}\ge 0$ in $E^o_{\tau, h}$ and $u_{\tau, h}\le0$ on $\partial_pE_{\tau, h}$ then there exists $h^* > 0$, such that as long as $0 < h\le h^*$ we have:

Proof. According to Kuo & Trudinger ^[17], the inequality (3.9) holds as long as the operator $\mathcal{L}^{\Delta}$ is evolving, monotone, time-wise non-degenerate, spatially non-degenerate and weak-positive. We will follow their definitions and show that these conditions can be satisfied if $h$ is sufficiently small.

Recalling $\mathcal{L}^{\Delta}$ in (3.4) and parameters $a_1-a_6$ defined by (3.5). According to Kuo & Trudinger ^[17], $\mathcal{L}^{\Delta}$ is evolving since we do not use point with $t > t_n$ in (3.4). And it is time-wise non-degenerate since the sum of the parameters for all points with $t < t_n$ is positive, that is, $a_6 > 0$.

Now we check the spatially non-degenerate condition. Denote

by definition from ^[17], $\mathcal{L}^{\Delta}$ is spatially non-degenerate if there is a constant $c > 0$ such that for all $\xi\in\Omega$ the following inequality holds:

Since $-1 < \rho < 1$, this can be achieved by directly calculation:

Then we check the monotone and weak-positive conditions. According to Kuo & Trudinger ^[17], $\mathcal{L}^{\Delta}$ is monotone if $a_1, a_2, a_3, a_4$ are all strictly positive, and $\mathcal{L}^{\Delta}$ is weak-positive if it is monotone and $\sum_{i = 1}^5a_i < 0$. Noting that $r, \sigma, \kappa, \theta > 0$, $|\rho| < 1$, $\sigma^2 < 2\kappa\theta$ and $y > \delta$, we have:

It is easy to check $a_i > 0$ as long as

which indicates that $\mathcal{L}^{\Delta}$ is monotone. And is it quit directly to show that $\mathcal{L}^{\Delta}$ is weak-positive since $\sum_{i = 1}^5a_i = -\frac{1}{\tau} -\frac{r}{y} < 0$.

Finally, we define

to finish the proof. ■

Equipped with the Lemma 3.3, we can prove the convergence of our finite difference scheme.

Proof of Theorem 3.1. Define a perturbation function with parameter $\epsilon > 0$ to $u_b$:

We can check that

Considering the truncation error given by (3.7), for all $\epsilon > 0$, there exists $\tau_*, h_*$ such that for all $0 < \tau < \tau_*$ and $0 < h < h_*$, we have

and thus

Since $u_b^{(\epsilon)}(x, y, t) = u_b(x, y, t)+\epsilon t \ge u_{\tau, h}$ at the boundary $\partial_pE_{\tau, h}$, by Lemma 3.3, we have:

Symmetrically, we define $u_{b}^{(-\epsilon)}(x, y, t) = u_b(x, y, t)+\epsilon t$. One can check that

Then we can say for all $\epsilon > 0$, there exists $\tau_*, h_*$ such that as long as $0 < \tau < \tau_*$ and $0 < h < h_*$,

which implies the convergence. ■

Proof of Theorem 3.2. According to Theorem 2.1, we know that

Then by the standard comparison principle, we can check that

Together with (3.7) and the Theorem 3.1, we have:

■

4.   Numerical experiments

In this part we implement some numerical experiments to verify the results established by the previous sections. We choose the parameters of $\mathcal{L}$ from Liang et al. ^[13], shown in Table 1.

We will solve the finite difference problem (3.6), which is corresponding to the bounded domain problem (2.8) and the original pricing problem (1.1), in two cases, $K = 0$ and $K = 1$, respectively. Then we show the effect of artificial boundary value error and the accuracy and convergence speed of the numerical method.

4.1.   Case of K=0

We start from a trivial case, in which the strike price $K = 0$. With $K = 0$, one can check that $u(x, y, t) = e^x$ is the solution to (1.1). Consider the following finite difference problem (3.6) corresponding to $K = 0$:

where $E_{\tau, h}$ is the time-space mesh of $(\Omega_{\delta, R_1, R_2})_\mathcal{T}$ defined through (3.1).

Clearly, if we choose $\psi_a(x, y, t) = e^x$, that is there is no error in the artificial boundary value, then the error between $u_{\tau, h}$ and the exact solution $u$ of (1.1) only comes from the truncation error, which is controlled by the meshing parameters $\tau$ and $h$. In this part, we first set $\psi_a = e^x$ and calculate $u_{\tau, h}$ with different $\tau, h$ to verify the convergence speed of our method. Then we fix $\tau$ and $h$, but deliberately choose $\psi_a$ with error to observe the effect of boundary value errors.

Given $\tau_0 > 0, h_0 > 0$, the spatial discrete error and time discrete error are defined as:

where $\tau > \tau_0, h > h_0$.

We implement our numerical method on $E_{\tau, h}\subset(\Omega_{\delta, R_1, R_2})_\mathcal{T}$ with geometric parameters $\delta = 0.01,$ $R_1 = 2, R_2 = 1$.

To show the convergence speed of $e_\tau$ and $e_h$, we choose $\psi_a(x, y, t) = e^x$ and set $\tau_0 = 1.25\times10^{-3}$, $h_0 = 5\times10^{-2}$, and let $\tau$ varies from $1.45\times10^{-3}$ to $2\times10^{-3}$, $h$ varies from $5\times10^{-2}$ to $5\times10^{-1}$. In particular, we choose $\tau, h$ to be integer multiples of $\tau_0, h_0$ to ensure the grid points of $E_{\tau, h}$ are included in that of $E_{\tau_0, h_0}$. The numerical results are shown in Figure 1. It can be seen that $e_\tau$ and $e_h$ converges linearly with respect to $\tau$ and $h^2$, respectively, which is consistent with the truncation error of the finite difference scheme given by (3.7).

Now we discuss the effect of boundary errors. The grid parameters $\tau_0, h_0$ and the geometric constants $\delta, R_1, R_2$ are set as before. Let $u_0(x, y, t)$ be the solution of (4.1) under $\psi_a = e^x$ on $E_{\tau_0, h_0}$. Then we choose two artificial boundary value conditions $\psi_a^{(1)}$ and $\psi_a^{(2)}$ with error, where $\psi_a^{(1)} = e^x + 1$ has unit error and $\psi_a^{(2)} = 2e^x$ has exponential error. We calculate $u_{\tau_0, h_0} - u_0$ under the artificial boundary values $\psi_a^{(1)}$ and $\psi_a^{(2)}$ respectively. The numerical results are shown in Figure 2. We represent the values of $u_{\tau_0, h_0} - u_0$ with a color gradient.

As can be seen from Figure 2(a), there is a constant error equal to one on the boundary of $E_{\tau_0, h_0}$, but the error rapidly decays to a level close to zero inside the solution domain far from the boundary. Figure 2(b) shows that the exponential error have a similar decay. These results are consistent with the conclusion in subsection 2.3 that the boundary value errors converge rapidly inside the solution domain.

4.2.   Case of K=1

In this part, we consider a non-trivial case with $K = 1$. The corresponding finite difference problem is defined as follows:

where $E_{\tau_1, h_1}$ is the time-space mesh of $(\Omega_{\delta, R_1, R_2})_\mathcal{T}$. The grid parameters $\tau_1 = 2\times10^{-5}, h_1 = 2.5\times10^{-3}$ and the geometric parameters $\delta = 2.5\times10^{-4}, R_1 = 2$, $R_2 = 1$. We follow the statement of the Theorem 3.2 and choose two artificial boundary value conditions $\psi_a^{(u)} = e^x$, $\psi_a^{(l)} = (e^x - 1)^+$, then $u_{\tau_1, h_1}^{(u)}, u_{\tau_1, h_1}^{(l)}$ and $\hat{u}_{\tau_1, h_1}$ are defined accordingly. Figure 3 shows the difference between $u_{\tau_1, h_1}^{(u)}$ and $u_{\tau_1, h_1}^{(l)}$. The color gradient represent values of $log(|u_{\tau_1, h_1}^{(u)}(x, y, \mathcal{T}) - u_{\tau_1, h_1}^{(l)}(x, y, \mathcal{T})|)$.

The contour lines show the difference between the solutions under two different boundary value functions gradually diminish when leaving the boundary. Typically, such difference is less than $1\times10^{-5}$ in the area enclosed by the innermost contour, where we can find that the error comes from boundary values is less than $5\times10^{-6}$ and it has the same magnitude as the truncation error. In practical application, such small error is negligible.

To further demonstrate the accuracy of our method, we compare our results with those obtained from the Monte Carlo method, the asymptotic approximation method, and the path-integral method, respectively. Results obtained from the Monte Carlo method and the path-integral method are given by Liang et al. ^[13], and the results obtained from the asymptotic approximation method are given by Li and Mecurio ^[9]. In their work, numerical solutions for the pricing problem (1.1) are provided with the same parameters as ours at one single point: $u(x_0, y_0, \mathcal{T}) = 0.125789$, where $x_0 = 0$ and $y_0 = 0.087$. Since $y = 0.087$ is not on the mesh grid of $E_{\tau_1, h_1}$, we find the numerical solution there by a standard second order interpolation and the results are shown in Figure 4.

As can be seen from Figure 4(a), the red line representing $u_{\tau_1, h_1}^{(l)}$ and the green line representing $u_{\tau_1, h_1}^{(u)}$ are almost indistinguishable from the blue line representing $\hat{u}_{\tau_1, h_1}$ when $x\in[-1, 1]$. Figure 4(b) further shows the difference between $u_{\tau_1, h_1}^{(u)}$ and $u_{\tau_1, h_1}^{(l)}$, which according to (3.8) mainly controls the accuracy of $\hat{u}_{\tau_1, h_1}$. Typically, at the point of $(x_0, y_0, \mathcal{T})$, the results of each method are compared in Table 2 below. Since the path-integral method is considered to be a relatively accurate method by Li and Mercurio ^[9], we can use its results as a benchmark to calculate the relative errors of the results of the other three methods.

It can be seen that the relative error of our method is an order of magnitude smaller than the Monte Carlo method and the asymptotic approximation method. The gap between our result and the results of the path-integral method proposed by Liang et al. ^[13] is only $3\times10^{-6}$, which is negligible in practical applications. Therefore our method can be considered to have sufficient accuracy.

It is worth noting that, unlike the path-integral method proposed by Liang et al. or the asymptotic approximation method proposed by Li and Mercurio, our method does not give the solution at a single point, but the solution at every point on $E_{\tau_1, h_1}$. If we choose a bounded domain $\mathcal{D} = [-0.5, 0.5]\times[0.05, 0.5]$, which is enclosed by the innermost contour of Figure 3. The option prices are shown in Figure 5. The error in $\mathcal{D}$ is on the order of $10^{-5}$.

5.   Conclusions

In this paper, we provide a finite difference method for pricing perpetual timer options under Heston volatility model. Since the pricing equation is degenerate, we first prove that the solution of the pricing problem exists and is unique under the asymptotic growth condition. Then we give a priori estimate of the artificial boundary value error for the bounded domain problem (2.8) and show that the solution of (2.8) internally converges to the solution to the original pricing problem (1.1). Furthermore, we provide a finite difference scheme with second-order convergence to solve the bounded domain problem (2.8). We implement our finite difference method with parameters in Liang et al. ^[13], and visualize the convergence and accuracy of this method.

Conflict of interest

The authors declare no conflict of interest.