A model-based method with geometric solutions for gaze correction in eye-tracking

1.   Introduction

An FX-TARN is an exotic financial product on a currency pair that saw a re-emergence of activity starting from early 2010s, with increased volumes in Europe, Asian and the Middle East. The TARN imposes positive and negative cash flows on scheduled dates (fixing dates) and the negative cash flow is typically leveraged by a gear factor. These cash flows can take the form of call or put option payoffs and the investor accumulates both positive (in the money) and negative (out of the money) cash flows until a certain target accrual level (cap) is reached. Once the target has been breached, the product is terminated (knocked out) and all the future cash flows are canceled but the accumulated ones are kept. For the last payment on which the target is reached, there are three knock-out types: no gain, part gain or full gain, which are specified in the contract at the initiation of the trade. The TARN is traded mostly by corporations and investors and often seen as a low-enhanced hedge for currency exposures. The TARN transaction allows the customer to exchange one currency for another at a contract rate that is more attractive compared to the rate on a traditional forward contract. However, the higher rate is accompanied with a higher level of downside currency exchange risk if the exchange rate were to move in the wrong direction from that expected. If spot rate moves in the wrong direction, the holders can be forced to trade regularly at unfavourable rates for the full life of the product.

The TARN is a path dependent product as the payments on the fixing dates depend on the FX spot rate as well as the accumulated payment amount up to the knock-out date. The TARN is typically priced with Monte Carlo (MC) simulations (see e.g. Caspers, 2015) but there are also some developments in the finite difference methods. Luo and Shevchenko (2015) describe a finite difference scheme for pricing a TARN under Black-Scholes (BS) model. The key steps involve tracking multiple one dimensional finite difference solutions, applying jump conditions at each fixing date and a cubic spline interpolation of results after each jump. Luo and Shevchenko (2015) demonstrate that the finite difference method can be faster than the MC method by an order of magnitude while achieving the same accuracy in price. Bandelier (2017) enlarges the finite difference method of Luo and Shevchenko (2015) to the Lévy processes with jumps. This leads to large complexities in the implementation of the method and the computational time grows up with the complexity. Bandelier (2017) then combines the method of Luo and Shevchenko (2015) and the convolution method proposed by Lord et al. (2008) where FFT-based method is applied to early-exercise options. Arregui and Ráfales (2020) generalize the work of Luo and Shevchenko (2015) by introducing a stochastic local volatility (SLV) technique and solve a partial differential equation for TARN by a finite difference alternating directions implicit (ADI) method.

In this paper, we propose a new method to solve the TARN pricing problem under exponential Lévy processes. Our method is based on a combination of the Fourier cosine expansion and high-order quadrature and belongs to numerical integration methods (also referred to as transform methods).

The COS method is an alternative to finite difference method and works particulary well for option pricing problems under Lévy processes. The COS method was initiated by Fang and Oosterlee (2008) and the key insight is the close relation of the characteristic function with the series coefficients of the Fourier-cosine expansion of the density function. In most cases, the convergence rate of the COS method is exponential and the computational complexity is linear for European options. Fang and Oosterlee (2009) extend the COS method of Fang and Oosterlee (2008) to early-exercise and discretely-monitored barrier options under exponential Lévy asset price models. They show the error convergence is exponential for processes characterized by very smooth transitional probability density functions. They also propose the FFT method for the efficient calculation of the Fourier coefficients. The extension of one-dimensional COS method to two-dimensional has been investigated by Fang and Oosterlee (2011), Ruijter and Oosterlee (2012) and Zhang and Oosterlee (2014). In Fang and Oosterlee (2011), the pricing of Bermudan and barrier options under the Heston model is dealt with by a combination of a Fourier cosine expansion and high-order quadrature rules. The error analysis and numerical experiments confirm a fast error convergence. Zhang and Oosterlee (2014) propose a pricing method for early-exercise Asian options by a two-dimensional integration and a backward recursion of the Fourier coefficients. In Zhang and Oosterlee (2014), numerical techniques such as Fourier cosine expansion, Clenshaw-Curtis quadrature and the FFT are employed. The rapid convergence of the pricing is analysed and demonstrated by various numerical examples. Ruijter and Oosterlee (2012) propose a two-dimensional COS method different from Fang and Oosterlee (2011) and Zhang and Oosterlee (2014) in that Fourier cosine series expansions are applied in both dimensions, rather than just one dimension. As a result, the method requires the availability of the bivariate characteristic function of the two state variables.

We apply the two-dimensional COS methods of Fang and Oosterlee (2011) and Zhang and Oosterlee (2014) and extend them to the TARN pricing problem. The two state variables are the spot rate and the cumulated positive cash flow. To value the path-dependent options such as TARN, we need to know the joint distribution of two state variables for which an analytical formula does not exist but we can deduce it from the Fourier domain. We cannot derive the joint characteristic function for the two state variables and therefore, the two-dimensional COS method of Ruijter and Oosterlee (2012) is non-applicable here. Fortunately, we can employ the COS methods of Fang and Oosterlee (2011) and Zhang and Oosterlee (2014) as the joint distribution can be factored into two marginal conditional density functions that can be approximated by Fourier cosine expansions separately. To recover the Fourier coefficients recursively, we approximate the two-dimensional integration by higher-order quadratures such as Gauss-Legendre or Clenshaw-Curtis quadrature for the integration over the spot rate. We also need to extend the existing COS methods to suit the specific requirements for TARN pricing. For the quadrature to converge exponentially, we need the integrating function to be smooth. Unlike the situations in Fang and Oosterlee (2011) and Zhang and Oosterlee (2014), where the integrating function is always smooth, the integrand in our TARN pricing problem is only piecewise smooth and we solve this issue by splitting the integral according to the points of discontinuity. For the COS methods, the Fourier coefficients can be computed from matrix and vectors multiplications. In Fang and Oosterlee (2009), Fang and Oosterlee (2011), Zhang and Oosterlee (2014) and other COS applications, a key matrix appears repeatedly and it plays a critical role for the efficiency of the COS method. This matrix can be expressed as a sum of Hankel and Toeplitz matrices and it is well-known (see e.g. Gohberg and Olshevsky, 1994) that the matrix-vector multiplication for Hankel and Toeplitz matrices can be carried out efficiently via FFT. For the TARN pricing, we need to distinguish between the full/no gain case and part gain case. This distinction is made because the Fourier coefficients involve two different matrices for the two cases. In the full/no gain case, the same matrix appearing in Fang and Oosterlee (2009), Fang and Oosterlee (2011), Zhang and Oosterlee (2014) emerges again and we can therefore apply FFT method as before. In the part gain case, we encounter a new matrix and we are able to prove that it can still be decomposed into a Hankel matrix and a Toeplitz matrix and therefore, FFT is applicable in this case as well.

We also investigate the performance and accuracy of the two-dimensional COS method through a couple of numerical experiments. For comparison purpose, we also run some MC simulations. We focus on three important models based on Lévy processes, namely BS, Merton Jump Diffusion (JD) and NIG models. The experiments indicate that the COS method converges fast and in order to achieve the same level of accuracy, MC method will require a significantly large number of simulations. We also vary some parameters that define the TARN product and demonstrate that the changes in the price produced using the COS method are consistent with the expectations.

It is well-known that the convergence rate of the COS method depends on the smoothness of the underlying probability density function and the payoff function; see e.g. Fang and Oosterlee (2008), Fang et al. (2010), Ruijter et al. (2015) and Arias et al. (2022). A typical example is VG process, which is a special case of the Lévy process. Fang and Oosterlee (2008) show that for a peaked VG density, the convergence of the COS method is no longer exponential but geometric. For these cases, Ruijter et al. (2015) apply spectral filters to achieve faster convergence. In this paper, it is not necessary to apply the spectral filters of Ruijter et al. (2015) as the density functions of the three specific examples of Lévy processes in the paper have no discontinuity and after we split the integral according to the points of discontinuity, the payoff function becomes smooth.

The rest of the paper is organized as follows. In section 2, we introduce Lévy process and the pricing framework. We provide a description of TARN product in section 3. The main part of the paper is in section 4, where we discuss the COS method for TARN. We distinguish between the cases of full/no gain and part gain as the pricing formulas differ in these two cases. In section 5, we test the accuracy and performance of the COS method through various numerical experiments. Section 6 concludes. We also provide a brief introduction to Clenshaw-Curtis quadrature in the appendix.

2.   Valuation model

Let $(\Omega, \mathcal{F}, Q)$ be a probability space with an information filtration $(\mathcal{F}_t)$. We assume both domestic short rate $r_d$ and foreign short rate $r_f$ are constants. In addition, under the risk neutral probability measure $Q$, we assume the dynamics of FX spot rate $S(t)$ follow an exponential Lévy process:

where $X(t)$ is defined as:

where $L(t)$ is a Lévy process and the deterministic function $\xi(t)$ is selected so that the following no-arbitrage condition is satisfied:

The Lévy process $L(t)$ is infinitely divisible and its characteristic function can be expressed using Lévy-Khintchine formula (see e.g. Sato, 1999).

Proposition 1. The characteristic function of the Lévy process $L(t)$ has the form

where $\psi_{L}(u)$ is the characteristic exponent of the Lévy process and given by

where $\mu\in \mathbb{R}$, $\sigma\in \mathbb{R_+}$ and $\nu$ is a positive Radon measure on $\mathbb{R}\setminus \{0\}$ satisfying:

Given the characteristic function of $L(t)$, we can determine the function $\xi(t)$ from the following result:

Lemma 1. The function $\xi(t)$ in (2) is given by

We can also obtain the $n$-th cumulant $\kappa_n$ of the process $X(t)$ from the $n$-th derivative of the characteristic function of $X(t)$, that is

The class of Lévy processes encompasses Gaussian, jump-diffusion and pure jump processes. In our numerical studies in section 5, we focus on Gaussian, Merton JD and NIG processes. When the process $L(t)$ follows a Gaussian process, we obtain the BS model. Under this model, the implied volatility surface will be flat, which contradicts the observed smile shaped (or skew shaped) volatility surface. On the other hand, Merton JD process is a jump-diffusion process and NIG process is an infinite activity pure jump process and both of them can reproduce the observed volatility smile or skew. The NIG process can be considered as an Inverse Gaussian time-changed Gaussian process. For more information on Lévy processes, we refer to Schoutens (2003).

In the following examples, we list the characteristic exponent $\psi_{L}(u)$, the deterministic function $\xi(t)$ and the cumulants $\kappa_1$, $\kappa_2$ and $\kappa_4$ for Gaussian, Merton JD and NIG processes.

Example 2. Gaussian Process. This is a pure diffusion process with $\mu = 0$. Furthermore,

Example 3. Merton JD Process. This is a jump-diffusion process with $\mu = 0$. Furthermore,

where $\zeta$, $\beta$ and $\delta$ are the parameters governing the Lévy measure $\nu$ in (5) and satisfying $\zeta > 0$. In addition,

Example 4. NIG Process. This is an infinite activity pure jump process with $\mu = 0$ and $\sigma = 0$. Furthermore,

where $\alpha$, $\beta$ and $\delta$ are the parameters governing the Lévy measure $\nu$ in (5) and satisfying $\alpha > 0$, $-\alpha < \beta < \alpha$ and $\delta > 0$. In addition,

3.   Product description

Let $t_0$ be the valuation date and $t_N$ the maturity date of a TARN contract. On each fixing date $t_n$, $n = 1, \ldots, N$, if the target level $U$ is not breached by the accumulated amount $A_n$, the positive payoff per unit of notional foreign amount is defined by

and the negative payoff by

where $X_n: = X(t_n)$, $E$ is a strike and $\gamma$ is a strategy on foreign exchange ($\gamma = 1$ for a call option and $\gamma = -1$ for a put option). Furthermore, $g$ is a leverage factor with the typical value of 1 or 2.

The accumulated amount $A_n$ is determined by the positive payoffs accumulated from $t_1$ to $t_n$, defined as

with $A_0 = 0$.

Now denote $\tilde{N}$ the first fixing date when the target is breached:

If $U$ is not breached, we set $\tilde{N} = N+1$.

On each fixing date $t_n$, where $n = 1, \ldots, N$, if $n < \tilde{N}$, the total payoff is determined by $C^+(X_n)+C^-(X_n)$. However, if $t_n = \tilde{N}$, the product is knocked out and all the future cash flows are canceled but the last payment $W(X_n, A_{n-1})$ at time $t_n$ will be determined by one of three knock-out types: no gain, part gain or full gain. We can write $W(X_n, A_{n-1})$ as

In summary, the total payoff $C_n$ on each fixing date $t_n$, where $n = 1, \ldots, N$, is determined by the following:

The price of the TARN at $t_0$ in domestic currency is then

where $N_f$ is the notional amount per fixing date and $B_d(t_0, t_n)^{-1}$ is the domestic discounting factor from the fixing date $t_n$ to $t_0$.

4.   COS method for TARN

From the description of the product payoff function, we can calculate the time $t_n$ value of the TARN in a recursive way:

where $I$ is defined as

and

where $f(x_{n+1}, a_{n+1}|x_{n}, a_{n})$ is the joint conditional density of $X_{n+1}$ and $A_{n+1}$. The integration limits for $a_{n+1}$ are derived based on (19), from which we have $a_{n+1} = a_{n}+ C^+(x_{n+1})$. Conditioning on $x_{n+1}$ and employing the fact that $a_{n}\ge 0$, we have $a_{n+1}\ge C^+(x_{n+1})$. Furthermore, before the FX-TARN is knocked out or matures, $a_n < U$. Therefore, $a_{n+1} < U+C^+(x_{n+1})$.

The joint conditional density of $x_{n+1}$ and $a_{n+1}$ can also be written as

where $f_X$ and $f_A$ are the marginal conditional densities of $X_{n+1}$ and $A_{n+1}$, respectively. In (27) we have also used the facts

The key idea of the two-dimensional COS method is to approximate $c(t_n, x_n, a_n)$ by following several steps of approximations:

1. Truncate the integration range from $(-\infty, +\infty)$ to $[a, b]$ for $x_{n+1}$;

2.Approximate the conditional density by truncated Fourier cosine expansion, based on characteristic functions;

3. Approximate the integral with respect to $x_{n+1}$ by a high-order quadrature.

When restricted to $[a, b]$, $f_X(x_{n+1}|x_{n})$ can be approximated by Fourier cosine expansion as follows:

where $\sum{\vphantom{\sum}}'$ indicates the first term of the series is divided by 2 and $\Phi_{X(t_{n+1}-t_{n})}$ is the characteristic function of $X(t_{n+1}-t_{n})$, which can be calculated using Lévy-Khintchine formula.

We can also approximate $f_A(a_{n+1}|x_{n+1}, a_{n})$ using the truncated Fourier cosine expansion as

Plugging (29) and (30) into (26), we obtain the first COS approximation of $c(t_n, x_n, a_n)$:

The last part of (31) involves a two-dimensional integration and we will follow Fang and Oosterlee (2011) and Zhang and Oosterlee (2014) to numerically integrate over $x_{n+1}$ by some high-order quadrature. In Fang and Oosterlee (2011), Gauss-Legendre quadrature is proposed and their numerical experiments demonstrate that the accuracy of the lower-order equidistant rules such as composite trapezoidal rule were insufficient compared with Gauss quadrature. In Zhang and Oosterlee (2014), Clenshaw-Curtis quadrature is preferred as it appears to be computationally cheaper than Gauss quadrature. In the appendix, we provide a summary of Clenshaw-Curtis quadrature that has been implemented in the numerical experiments in section 5.

For both Gauss and Clenshaw-Curtis quadratures, if the integrand function is smooth, the quadrature will achieve an exponential convergence. This is indeed the case for the problems studied in Fang and Oosterlee (2011) and Zhang and Oosterlee (2014). In our case, however, the integrand function is only piecewise smooth and the discontinuity point can be identified at $E^* = \log(E/S(0))$ for $x_{n+1}$ by inspecting the functions $C^+(x_{n+1})$ and $C^-(x_{n+1})$, defined in (17) and (18), respectively. We then split the integration range $[a, b]$ into two intervals $[a^{(k)}, b^{(k)}]$, for $k = 1, 2$. We set $a^{(1)} = a$, $b^{(1)} = a^{(2)} = E^*$ and $b^{(2)} = b$. For each double integral, we apply ($n^{(k)}+2$)-point quadrature integration rule to the outer integral. Then the double integral in (31) can be approximated by

where $w_n^{(k)}$ and $\delta_n^{(k)}$ are the weights and quadrature nodes corresponding to the $k$-th integral with respect to $x_{n+1}$. In addition, $V_{n+1, j}(\delta_{n})$ is the Fourier coefficient defined by

where $v_{n+1}(\delta_n)$ denotes the value of $v_{n+1}$ evaluated at $x_{n+1} = \delta_n$.

By combining Fourier cosine expansion and high-order quadrature approximation, we arrive at the second COS approximation of $c(t_n, x_n, a_n)$:

Following Fang and Oosterlee (2009), a suitable integration range for $x_{n+1}$ can be taken in terms of the cumulants:

where $\lambda$ is a constant and $\kappa_k$ is the $k$-th cumulant of $X(t_N)$, which can be obtained from (8).

We note that recently, Junike and Pankrashkin (2022) derive a new formula for the truncation range based on Markov's inequality and prove that the range is large enough to ensure convergence of the COS method. Junike and Pankrashkin (2022) also show that the method of Fang and Oosterlee (2009) may lead to serious mispricing. However, the examples in Junike and Pankrashkin (2022) are based on simple European options. More tests are need to demonstrate the performance of the method for the path-dependent options, such as Bermudan and barrier options.

The valuation of TARN based on the COS method then boils down to the determination of the Fourier coefficients, which need to be solved recursively. We distinguish between the TARN with full/no gain case and part gain case since for the latter the function $W(X_n, A_{n-1})$ defined in (21) is dependent on $A_{n-1}$ whereas for the former, it is not. This distinction is made because the Fourier coefficients involve two different matrices for the two cases.

4.1.   TARN with full/no gain

In this case, for simplification purpose, we can write $W(X_n, A_{n-1})$ defined in (21) as $W(X_n)$ as it is independent of $A_{n-1}$. We can also write the function $I(X_n, A_n, A_{n-1})$ defined in (25) as $I(X_n, A_n)$. We first define the following function

The above integral can be expressed in closed-form from the following lemma.

Lemma 5.

Another integral that has appeared in the COS applications repeatedly is

Fang and Oosterlee (2009) prove the following result:

Lemma 6.

where

and

To derive the analytical formulas for the Fourier coefficients, we also define the following integral:

We can derive the formula for $H_j(\delta_n)$ as follows.

Lemma 7.

where

Proof. Using the definition of $I(x_n, a_n)$ in (25), we have

Using the definition of $\Psi_j$ in (36), we obtain the result.

The $t_0$ price of TARN computed with COS method can then be obtained from the following result.

Proposition 2. The price of the TARN with full/no gain under the two-dimensional COS method given $X(0) = x_0$ can be calculated by

where the Fourier coefficients $V_{n, j}(\delta_{n})$ can be obtained recursively as follows:

For $1\le l\le N-1$,

Proof. From (24) and (26), we know

The approximation of $v_{N-1}$ based on the COS method becomes

where

Using (42), we immediately obtain

For $0\le l\le N-1$, the COS approximation of $v_{l}$ is

where for $1\le l\le N-1$,

Finally, we can obtain the formula for $\hat{v}_0$ by realizing $a_0 = 0$. □

The Fourier coefficients $V_{l, j}(\delta_{n})$ can be computed efficiently. First, we define two $N_2$-dimension vectors $B^k_{l, n}$ for $k\in \{1, 2\}$ with elements:

Also define two $N_2\times N_2$ matrices $M^H_{l, n}$ and $M^T_{l, n}$ with elements:

and

It is clear that $M^H_{l, n}$ is a Hankel matrix and $M^T_{l, n}$ is a Toeplitz matrix.

Define the Fourier coefficient matrix $V_l$ with elements $V_l(j, n) = V_{l, j}(\delta_{n})$. We can compute $V_l$ by

where the $n$-th column of $V_l^{(1)}$ can be computed by $(M^H_{l, n}+M^T_{l, n})(B^1_{l, n}+B^2_{l, n})$ and the matrix $V_l^{(2)}$ has elements $V_l^{(2)}(j, n) = H_j(\delta_n)$.

It is well known (see e.g. Gohberg and Olshevsky, 1994) that matrix-vector product, with the special matrices such as Hankel and Toeplitz matrices, can be transformed into circular convolutions. Therefore, the Fourier coefficients can be computed efficiently using FFT. We refer to Fang and Oosterlee (2009) for the details about the implementation of FFT algorithm for the matrix-vector product.

4.2.   TARN with part gain

In this case, the function $W(X_n, A_{n-1})$ defined in (21) depends on $A_{n-1}$. For the efficient implementation of the COS method, we encounter a new integral defined by

We prove the following result:

Lemma 8.

where

where $M^H_{k, j}(x_1, x_2, a, b)$ and $M^T_{k, j}(x_1, x_2, a, b)$ are defined in (40) and (41), respectively.

Proof. We define $\bar{Q}^H_{k, j}(x_1, x_2, a, b, c)$ and $\bar{Q}^T_{k, j}(x_1, x_2, a, b, x)$ as two functions of $c$, where

and

where

and

We have

It is easy to show that

where

We finish the proof by setting

□

To obtain analytical formulas for the Fourier coefficients, we define the following integral:

We can derive the formula for $\bar{H}_j(\delta_n, a_{k-1})$ as follows.

Lemma 9.

where

and

Proof. Using the definition of $I(\delta_n, a_k, a_{k-1})$ in (25) and (21), we have

□

We provide the analytical formula for the $t_0$ price of TARN computed with COS method in the following proposition.

Proposition 3. The price of the TARN with part gain under the two-dimensional COS method given $X(0) = x_0$ can be calculated by

where the Fourier coefficients $V^{(1)}_{1, j}(\delta_{n})$ can be obtained recursively from $V^{(1)}_{l, j}(\delta_{n})$ and $V^{(2)}_{l, j}(\delta_{n})$ for $1\le l\le N$, described as follows.

and

For $1\le l\le N-1$,

and

Proof. The approximation of $v_{N-1}$ based on the COS method is

where

Using Lemma 9, we immediately obtain

where

and

For $0\le l\le N-1$, the COS approximation of $v_{l}$ is

where for $1\le l\le N-1$,

where

and

□

As in the case of full/no gain, the Fourier coefficients $V_{l, j}(\delta_{n})$ in the part gain case can be computed using matrix-vector products. However, in this case, besides Hankel matrix $M^H_{l, n}$ and Toeplitz matrix $M^T_{l, n}$ defined in (69) and (70), respectively, we require two $N_2\times N_2$ matrices $Q^H_{l, n}$ and $Q^T_{l, n}$ with elements:

and

From (74) and (75), it is clear that $Q^H_{l, n}$ is a Hankel matrix and $Q^T_{l, n}$ is a Toeplitz matrix. Therefore, the Fourier coefficients can be computed efficiently using FFT for TARN with part gain as well.

5.   Numerical results

In this section, a number of numerical experiments are performed to evaluate the accuracy and efficiency of the two-dimensional COS method for the pricing of TARN under the exponential Lévy processes. We will consider three different knock-out types: no gain, part gain and full gain. We focus on three specific models: BS, Merton JD and NIG, with the parameters listed below:

● BS: $r_d = 0$, $r_f = 0$, $\sigma = 0.2$.

● Merton JD: $r_d = 0$, $r_f = 0$, $\sigma = 0.2$, $\zeta = 3$, $\beta = -0.05$, $\delta = 0.05$.

● NIG: $r_d = 0$, $r_f = 0$, $\delta = 0.2$, $\alpha = 20$, $\beta = -5$.

In practice, to obtain the parameters for the Lévy processes, one proceeds to obtain a set of market prices of put and call options on FX and minimize the mean square error between market prices of the options and the corresponding prices predicted by the model.

For the TARN trades, we assume the time interval is a constant, i.e. $t_i-t_{i-1} = \Delta$, for $i = 1, \ldots, N$. The baseline parameters for the TARN trades are as follows.

● $S(0) = 1.05$, $E = 1$, $g = 2$, $\gamma = 1$, $\Delta = 1/12$, $t_0 = 0$, $t_N = 1$, $N_f = 1$.

For the COS method, we set $\lambda = 10$ throughout all the experiments. We implemented both Clenshaw-Curtis and Gauss-Legendre quadrature rules and found them to perform equally well. As a result, we only report the results based on Clenshaw-Curtis quadrature rule. For the quadrature, we set $n^{(1)} = n^{(2)} = 2^8$ and therefore $w_n^{(1)} = w_n^{(2)}$. Furthermore, in the numerical experiments, we set $N_1 = N_2 = 2^d$ and evaluate the performance of the COS method by varying $d$. $N_1$ and $N_2$ represent the truncation terms for the approximation of the marginal conditional densities of $X_{n+1}$ and $A_{n+1}$, repectively; see (29) and (30). Therefore, by varying $d$, we are testing the convergence of the COS method with respect to the truncation terms for marginal conditional densities of $X_{n+1}$ and $A_{n+1}$.

For comparison purpose, we also perform MC simulations for each model. The simulations are based on $200,000$ paths which are considered as reasonable in practice. We report both average and the standard error for the simulation results.

In Tables 1–3, we report the prices of TARN based on COS method by varying the target $U$ and the parameter $d$. We also report the corresponding results from MC. We set the prices obtained under $d = 11$ as the reference values and report in parentheses the absolute differences between the reference values and those from MC and COS under other values of $d$. For all of three models, it is clear that the price under the COS method converges fast. The price calculated based on $d = 7$ is already a good estimate for the reference value. In fact, the price difference computed from $d = 7$ and $d = 11$ is at most $0.0005$ for the BS model, $0.0006$ for the Merton JD model and $0.0006$ for the NIG model. When we raise $d$ to $9$, the price difference from the reference value is at most $0.0001$ for three models. We also confirm that all the prices computed under COS method with $d = 11$ lie within the $95%$ confidence interval formed based on MC. However, the minimum standard error from MC is 0.0015, which indicates that MC should increase by a large factor to achieve the same accuracy as the COS method since the MC standard error is proportional to the inverse square root of number of simulations.

The results in Tables 1–3 confirm that as expected, the TARN price under full gain is always highest, followed by part gain and then no gain. We need to emphasize that in the paper, we focus on the derivation of the analytical formula for the FX-TARN price. In order to hedge the risks in the FX-TARN, it is necessary to calculate the option Greeks, such as delta and gamma. An important advantage of the COS method over MC is that we can obtain the precise values for the options Greeks easily based on the analytical formulas instead of numerical approximation.

To further verify the accuracy of the COS method, we perform some additional tests under the three specific models. It turns out that the test results under the three models are consistent with each other and therefore, we only report the results under the selected models in order to save the space.

In the first experiment, we set $E = 0.6$, $\Delta = 1/4$ and the remaining parameters are same as the baseline ones. In this case, $S(0) = 1.05$, which is reasonably bigger than $E$, so we expect when the target $U$ is small enough, the TARN will always get knocked out on the first fixing date $t_1$. As we incrementally increase $U$, we expect the following things:

● Initially, the price under no gain will remain zero;

● Initially, the price under partial gain will be a linear function of the incremental increase in $U$;

● Initially, the price under full gain will stay at the level determined by the discounted value of payoff at $t_1$;

● When we continue to increase $U$, the price under the three knock-out types will all increase.

Figures 1–3 produce the results which are as expected. In the second experiment, we take the same parameters as the first experiment. When the target $U$ is big enough, the TARN will never be knocked out and in this case, the price of the TARN will converge to the same value determined by the underlying calls and puts, irrespective of knock-out types. The results in Figure 4 indicate it is indeed the case.

In the third experiment, we consider the TARN with part gain. We set $U$ = 0.9 and vary $E$ and $\Delta$, while fixing other parameters at the baseline ones. When $\Delta$ increases, we have more fixing dates since the maturity is fixed. In this case, the change in the price will depend on $E$. When $E$ is small, the payoff function at each fixing date tends to be positive. However, as $E$ increases, the payoff tends to be negative. Therefore, when $E$ is smaller and we increase $\Delta$, we expect the price will rise. On the other hand, when $E$ is bigger and we increase $\Delta$, the price will drop. We confirm these expectations in Figures 5.

In the final experiment, we investigate the impact of the leverage factor $g$. We consider a part gain TARN with $U = 0.9$, $\Delta = 1/4$ and the rest of parameters are same as the baseline ones. We expect that when $g$ switches from 1 to 2, the price will drop as the payoff on each fixing date becomes lower. The result in Fig. 6 is consistent with our expectation. We also observe that when $E$ is small, the prices under the cases of $g = 1$ and $g = 2$ are almost the same. This is because the negative payoff term is close to 0 when $E$ is small. Before the TARN is knocked out, its payoff at each fixing date with $g = 1$ minus that with $g = 2$ is in fact the payoff of a put option. As the diffusion parameter $\sigma$ takes same values under both BS and Merton JD models, we expect that the price difference between the cases of $g = 1$ and $g = 2$ in Merton JD model is larger than that in BS model. This is because in the Merton JD model, the spot rate has a larger volatility due to the existence of the jumps, which in turn leads to a higher put option value. The result in Figures 7 is consistent with our expectation.

6.   Conclusions

In this paper, we develop an efficient pricing method for FX-TARN under Lévy processes, based on the two-dimensional Fourier cosine expansion combined with higher-order quadrature rules. We derive the analytical formulas for TARN price under three different knock-out types. Furthermore, we demonstrate that the Fourier coefficients can be computed fast by using FFT. The performance and accuracy of the method have been confirmed through various numerical experiments.

As for future work, it would be interesting to extend the method for the model with both stochastic volatility and jumps. One such an example is the time changed Lévy model of Carr and Wu (2004). For this type of models, the characteristic function is still available and therefore, the COS method should continue to work.

Appendix: Clenshaw-Curtis quadrature

Clenshaw–Curtis quadrature is the method for numerical integration and it is a quadrature based on an expansion of the integrand in terms of Chebyshev polynomials. Besides having fast-converging accuracy comparable to Gaussian quadrature rules, Clenshaw–Curtis quadrature naturally leads to nested quadrature rules (where different accuracy orders share points).

A ($N+2$)-th order Clenshaw-Curtis rule takes the following form

where $w_n$ and $\delta_n$, for $n = 0, \ldots, N+1$ are the weights and nodal points, which are defined as

and $w_n$ is $n$-th element of a ($N+2$)-dimension vector $w = (D'd, D'd)'$, where $D$ is a $(N/2+1)\times(N/2+1)$-matrix, with elements:

and $d$ is a vector defined by $d_0 = 1$, $d_{N/2} = 1/(1-N^2)$, $d_k = \frac{2}{1-(2k)^2}$, for $k = 1, \ldots, N/2-1$ and can be written as

Use of AI tools declaration

The author declares that Artificial Intelligence (AI) tools haven't been used in the creation of this article.

Conflict of interest

The author declares no conflicts of interest in this paper.