The CTMC–Heston model: calibration and exotic option pricing with SWIFT

Álvaro Leitao, J. Lars Kirkby and Luis Ortiz-Gracia

Abstract

This work presents an efficient computational framework for pricing a general class of exotic and vanilla options under a versatile stochastic volatility model. In particular, we propose the use of a finite state continuous time Markov chain (CTMC) model, which closely approximates the classic Heston model but enables a simplified approach for consistently pricing a wide variety of financial derivatives. Under the CTMC–Heston model, we show that the shape of the implied volatility is preserved (hence, it has an equivalent ability to calibrate market smiles), yet it may price complex derivatives such as Asian options, variance swaps/options and cliquets with great efficiency. To accomplish this, we employ Markov chain approximation techniques, which have been gaining traction in the recent literature. We expect that this framework will have significant practical appeal, given the widespread adoption of the Heston model and the present difficulties that arise in pricing complex products.

Download PDF Appendix 1

Introduction

1 Introduction

The Heston model (Heston 1993) is one of the most widely used stochastic volatility (SV) models in the option pricing literature as well as in practice. For a fixed time horizon, the characteristic function (ChF) is known in closed form, and European option pricing is accomplished with any standard Fourier method (Carr and Madan, 1999; Fang and Oosterlee, 2008; Kirkby, 2015; Ortiz-Gracia and Oosterlee, 2016), or with a two-dimensional partial differential equation (PDE). The availability of fast pricing for European options enables calibration of the Heston model parameters to match observed volatility surfaces as required in practice. However, after calibration there is still great difficulty in pricing exotic contracts based on the Heston model, as is the case with most multifactor equity models. Numerous extensions of the Heston model have also been proposed – such as randomization of the initial variance (Jacquier and Shi 2019), fractional and rough volatility (Guennoun et al 2018; Euch and Rosenbaum 2018; Jaber and Euch 2019), and time-dependent correlation (Benhamou et al 2010) – yet pricing exotic derivatives even under the standard formulation is still very challenging. To price contracts such as Asian options and variance swaps, Monte Carlo (MC) methods are the traditional surrogates in these cases, since they can handle high-dimensional problems with relative ease. Unfortunately, MC suffers from a number of well-known deficiencies, and complicated simulation schemes are often required to overcome not only the boundary effects that accompany models such as Heston (Broadie and Kaya 2006; Andersen 2008), but also the lethargic convergence rates that necessitate substantial computational resources to solve problems at industry scale. Some specialized Fourier approaches have been developed for contracts such as barrier options under the Heston model (Fang and Oosterlee 2011; Ruijter and Oosterlee 2012), but at this stage there is a lack of more general purpose techniques that can handle complex products under this model.

Markov chain approximation has emerged in the last few years as an effective approach for solving pricing problems of moderate dimensionality. In the seminal work of Mijatović and Pistorius (2013), an approach is proposed to solve the barrier option problem for Lévy models in one dimension. Recently, extensions of Markov chain approximations have been derived to handle multifactor models such as stochastic volatility (Kirkby et al 2017) and stochastic local volatility (Cui et al 2018a). Other important recent works in this area include Li and Zhang (2017), Zhang and Li (2019) and Li and Zhang (2019), which derive rigorous error analyses for designing the state space of the univariate approximating Markov chain. Markov chain approximation techniques are now available for several exotic options, including realized variance derivatives (Cui et al 2017a), Asian options (Kirkby and Nguyen 2020) and cliquets (Cui et al 2017b). These approaches use a Markov chain approximation to obtain a regime switching (RS) model, which greatly simplifies option pricing. These techniques are closely related to RS methods, which already have a vast and growing literature (see Elliott et al 2007; Ramponi 2012; Costabile et al 2014; Boyarchenko and Levendorskii 2014; Jiang et al 2016; Elliott et al 2005; Dang et al 2016; Song et al 2016). This work also has ties with option pricing under general Markov processes (Cai et al 2015; Cui et al 2018b) as well as trees (Nguyen 2018; Bayraktar et al 2018). As Markov chain approximation techniques have proven their effectiveness, researchers have started to study ways to improve the main computational burden of the method, which is a matrix exponential computation (Cui and Taylor 2019). Over time, we expect that general computational improvements will further extend the applicability of this approach.

The practical objective of this work is to formalize a model that reproduces vanilla market quotes but, at the same time, is amenable to complex derivatives pricing in a manner that is consistent with the calibrated model. Here, we propose a model and framework to accomplish these goals, which is based on the ubiquitous Heston model. We call this the continuous time Markov chain–Heston (CTMC–Heston) model, as it uses a finite state CTMC approximation to the Cox–Ingersoll–Ross (CIR) variance process. We show that, even for a small number of variance states, this model produces volatility smiles that are quite similar to those of the Heston model and converge in the limit as the state space is refined. The advantage is that for a moderate number of states, exotic options are efficiently priced with a model that calibrates well to the market by utilizing the Shannon wavelet (SWIFT) machinery of Ortiz-Gracia and Oosterlee (2016) and Leitao et al (2018). In contrast to the standard approach taken in the RS literature, we impose a structure on the resulting RS model, particularly its transition matrix, through the Heston parameterization process. This approach reduces the required number of calibration parameters significantly and produces a model that possesses the inherent characteristics of a well-known (and understood) stochastic volatility model, yet is more tractable and efficient for pricing exotic derivatives. In this sense, we retain the benefits of the stochastic volatility formulation, but we gain the tractability of an RS model. Alternative approaches for modeling and calibrating RS models can be found in He and Zhu (2017), Mitra and Date (2010), Janczura and Weron (2012), Goutte et al (2017), He and Zhu (2018) and Jourdain and Zhou (2020).

Among the exotic contracts, we consider the arithmetic Asian option. Asian options have been widely studied in the pricing literature, with several powerful approaches – including spectral expansion (Linetsky 2004), Laplace inversion (Geman and Yor 1993; Cui and Nguyen 2017) and the use of analytical bounds (Fusai and Kyriakou 2016) – that are all limited to one-dimensional Markov models such as Lévy processes and square-root diffusions. Fourier methods in particular have seen tremendous success in recent years for pricing exotic options, especially under Lévy processes. For Asian options, the works of Levendorskii (2018), Zhang and Oosterlee (2013), Kirkby (2016) and Leitao et al (2018) have provided fast numerical schemes that rely on the ChF recursion introduced in Carverhill and Clewlow (1990), which effectively reduces the problem dimensionality for Fourier methods. The use of wavelets in Fourier-based option pricing has also gained momentum, with works by Ortiz-Gracia and Oosterlee (2013), Kirkby (2015) and Ortiz-Gracia and Oosterlee (2016). While these approaches were initially introduced to solve one-dimensional problems, there is great interest in extending these results to multifactor settings, such as stochastic volatility and stochastic interest rates (Colldeforns-Papiol et al 2017; Kirkby et al 2017; Dang and Ortiz-Gracia 2018).

Also included in our framework are discretely sample realized variance swaps and options (Bernard et al 2014; Bernard and Cui 2014; Pun et al 2015; Zheng and Kwok 2014; Zhu and Lian 2011), which are a more recent entrant to the class of exotic contracts. While closed-form pricing of variance swaps under the Heston model has been solved by Bernard and Cui (2014), variance and volatility options are very difficult contracts to price and require numerical methods.

The contributions of this work are as follows.

• •

  We demonstrate the effectiveness of the CTMC–Heston model for calibrating market parameters in a manner that is consistent with the Heston model yet tractable enough for exotic option pricing. We further analyze several competing grid designs and provide extensive numerical experiments comparing the trade-offs of each.

• •

  We propose an efficient Fourier approach for pricing exotic derivatives in the CTMC–Heston model, including Asian options, variance derivatives and cliquets. This method combines Markov chain approximation with the SWIFT Fourier method, resulting in an accurate and efficient pricing framework. Robustness to model parameters is demonstrated numerically.

• •

  We provide a detailed theoretical analysis of the method, including derivations that demonstrate convergence to the Heston model as the state space is refined. These results are supported by numerical studies that demonstrate convergence.

This paper is organized as follows. Section 2 introduces the CTMC–Heston model and the basic ingredients of the proposed framework. Section 3 introduces wavelet-based Fourier inversion using the SWIFT method, applies it to European option pricing under the CTMC–Heston model and provides error bounds. Section 4 conducts several calibration studies that focus on investigating the appropriate grid design for the underlying CTMC, the choice of $m_0$ and the influence of model parameters on convergence. Section 5 extends the framework to exotic contracts by providing a unified, general approach that applies to Asian options, variance swaps/options and cliquets, to name a few. Numerical experiments are provided. Section 6 concludes the work and suggests directions for future research.

2 From Heston model to CTMC–Heston model

This section describes the modeling methodology we employ in this work and the Markov chain approximation we use to simplify the SV dynamics in terms of RS models. For conciseness, we focus on the case of the Heston model, although the framework extends naturally to cover alternative SV models, as discussed in Kirkby et al (2017). Jumps are easily included as well: see, for example, the Bates model in Bates (1996), which is a natural extension of the Heston model.

2.1 The Heston model

Consider the Heston stochastic volatility model (Heston 1993):

(2.1)

where $\mathrm{d}{W}_t^1$ and $\mathrm{d}{W}_t^2$ are correlated Brownian motions, ie, $\mathrm{d}{W}_t^1\mathrm{d}{W}_t^2=\rho \mathrm{d}t$, with $\rho \in (-1,1)$. Here, $r$ denotes the interest rate and $q$ is the dividend yield.

The stochastic volatility (or variance), $v_t$, is driven by a CIR process, having a mean-reversion component. The value $v_0$ is the initial volatility, $\eta$ controls the mean-reversion speed, $\theta$ is the long-term volatility and ${\sigma }_v$ corresponds to the volatility of the variance process $v_t$: the so-called volatility-of-volatility (vol-vol). The model parameters are therefore $\mathrm{\Theta }=\{v_0,\eta ,\theta ,{\sigma }_v,\rho \}$. Note that to ensure $v_t>0$, Feller’s condition $2\eta \theta >{\sigma }_v^2$ is imposed (see Heston 1993). We assume the existence of a risk-neutral measure on $(\mathrm{\Omega },ℱ,\mathbb{Q})$, under which all expectations are taken.

The solution of (2.1) can be re-expressed in the form

	$\log \left({\displaystyle \frac{S_t}{S_0}}\right)$	$={\displaystyle \frac{\rho }{{\sigma }_v}}(v_t-v_0)+(r-q)t-{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^t}{v}_{s}ds-{\displaystyle \frac{\rho }{{\sigma }_v}}{\displaystyle {\int }_0^t}\eta (\theta -v_s)ds$
		$\mathrm{\hspace{1em}}+\sqrt{1-{\rho }^2}{\displaystyle {\int }_0^t}\sqrt{v_s}d{W}_s^{*},$		(2.2)

from the Cholesky decomposition $W_t^1:=\rho W_t^2+\sqrt{1-{\rho }^2}{W}_t^{*}$, where $W_t^{*}$ is independent of $W_t^2$.

After some rearranging, we introduce the auxiliary process ${({\tilde{X}}_t)}_{t\ge 0}$:

	${\tilde{X}}_t$	$:=\log \left({\displaystyle \frac{S_t}{S_0}}\right)-{\displaystyle \frac{\rho }{{\sigma }_v}}(v_t-v_0)$
		$=\left(r-q-{\displaystyle \frac{\rho \eta \theta }{{\sigma }_v}}\right)t+\left({\displaystyle \frac{\rho \eta }{{\sigma }_v}}-{\displaystyle \frac{1}{2}}\right){\displaystyle {\int }_0^t}{v}_{s}ds+\sqrt{1-{\rho }^2}{\displaystyle {\int }_0^t}\sqrt{v_s}d{W}_s^{*}.$

We thus have the following uncoupled two-factor representation with independent Brownian motions:

(2.3)

where $\overline{\omega }:=(r-q-\rho \eta \theta /{\sigma }_v)$. By uncoupling the two driving factors, we are able to approximate the joint dynamics with an RS diffusion, which is described next.

2.2 CTMC–Heston model

To simplify the model, we introduce a CTMC approximation of the variance dynamics, ${(v_t)}_{t\ge 0}$, which can be rewritten as

$$\mathrm{d}{v}_t=\mu (v_t)\mathrm{d}t+\sigma (v_t)\mathrm{d}{W}_t^2,$$

(2.4)

where $\mu (v_t):=\eta (\theta -v_t)$ and $\sigma (v_t):={\sigma }_v\sqrt{v_t}$. We fix a finite state space $?:=\{v_1,v_2,\mathrm{\dots },v_{m_0}\}$ and a CTMC $\{\alpha (t),t\ge 0\}$, which transitions between the state indexes $\{1,\mathrm{\dots },m_0\}$, given a time increment $\mathrm{\Delta }t$, according to

$$\mathbb{P}\{\alpha (t+\mathrm{\Delta }t)=j\mid \alpha (t)=k,\alpha (t^{\prime }),0\le t^{\prime }\le t\}={\delta }_{jk}+q_{jk}\mathrm{\Delta }t+o(\mathrm{\Delta }t),$$

where ${\delta }_{jk}=1$ if $j=k$, and zero otherwise. The set of transition rates $q_{jk}$ between states $j$ and $k$ is collected in the generator matrix $Q$ (of size $m_0\times m_0$) and chosen so that the dynamics of ${(v_{\alpha (t)})}_{t\ge 0}$ are locally consistent with those of ${(v_t)}_{t\ge 0}$, that is, by matching the first two dynamic moments as described in Section 2.3. So far, we have assumed that the state space $?$ is given. In Section 4.1, more insight into its optimal construction is provided.

Given the CTMC process ${(v_{\alpha (t)})}_{t\ge 0}$, we can now approximate the dynamics in (2.3) by an RS diffusion. We define

	${\tilde{X}}_t^{\alpha }$	$=\overline{\omega }t+{\displaystyle {\int }_0^t}\left({\displaystyle \frac{\rho \eta }{{\sigma }_v}}-{\displaystyle \frac{1}{2}}\right)v_{\alpha (s)}ds+\sqrt{1-{\rho }^2}{\displaystyle {\int }_0^t}\sqrt{v_{\alpha (s)}}d{W}^{*}(s)$
		$={\displaystyle {\int }_0^t}{\zeta }_{\alpha (s)}ds+{\displaystyle {\int }_0^t}{\beta }_{\alpha (s)}d{W}^{*}(s),$		(2.5)

where, for $\alpha (s)\in \{1,\mathrm{\dots },m_0\}$,

(2.6)

In particular, ${\tilde{X}}_t^{\alpha }$ is governed by an RS diffusion.

The advantage of this formulation is that it yields a closed-form expression for the conditional ChF. For any duration of time between transitions of the underlying CTMC, the dynamics in (2.5) behave as a simple diffusion with constant coefficients. Given a time increment of size $\mathrm{\Delta }t>0$, we define for each state $j=1,\mathrm{\dots },m_0$,

	${\tilde{\varphi }}_{{\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }}^j(\xi )$	$:=?[{\mathrm{e}}^{\mathrm{i}\xi {\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }}|\alpha (0\le s\le \mathrm{\Delta }t)=j]$
		$=?[\exp (\mathrm{i}\xi ({\zeta }_j\mathrm{\Delta }t+{\beta }_j{W}^{*}(\mathrm{\Delta }t)))]:=\exp ({\psi }_j(\xi )\mathrm{\Delta }t),$		(2.7)

where

$${\tilde{\varphi }}_{{\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }}^j(\xi )$$

is the ChF of a process that remains in state $j$ for the duration $\mathrm{\Delta }t$, and ${\psi }_j(\xi )$ is its Lévy symbol. By the independence of $\alpha (t)$ and $W^{*}(t)$, the Lévy symbols satisfy

$${\psi }_j(\xi )=\mathrm{i}{\zeta }_j\xi -\frac{1}{2}{\xi }^2{\beta }_j^2,j=1,\mathrm{\dots },m_0,$$

(2.8)

where ${\zeta }_j$, ${\beta }_j$ are defined in (2.6). In fact, the process ${\tilde{X}}_t^{\alpha }$ is completely characterized by the set of Lévy symbols in each of its possible states, ${\{{\psi }_j(\xi )\}}_{j=1}^{m_0}$, together with the generator $Q$. From Chourdakis (2004), the characteristic function of ${\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }$, $\mathrm{\Delta }t\ge 0$, conditioned on the initial state $\alpha (0)=j_0$, satisfies

$$?[{\mathrm{e}}^{\mathrm{i}{\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }\xi }\mid \alpha (0)=j_0]={\mathrm{?}}^{\prime }ℳ(\xi ;\mathrm{\Delta }t){?}_{j_0},j_0\in \{1,\mathrm{\dots },m_0\},$$

where we define the matrix exponential

$$ℳ(\xi ;\mathrm{\Delta }t):=\exp (\mathrm{\Delta }t(Q^{\prime }+\text{diag}({\psi }_1(\xi ),\mathrm{\dots },{\psi }_{m_0}(\xi )))),$$

(2.9)

and $\mathrm{?}\in {ℝ}^{m_0}$ represents a column vector of ones, while ${?}_j\in {ℝ}^{m_0}$ is a unit column vector with a one in position $j$. The ${}^{\prime }$ symbol indicates the matrix (vector) transpose operation.

Note that from (2.2), ${\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }$ induces the following CTMC–Heston model for the underlying $S_{\mathrm{\Delta }t}$, namely

$$S_{\mathrm{\Delta }t}^{\alpha }=S_0\exp \left({\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }+\frac{\rho }{{\sigma }_v}(v_{\alpha (\mathrm{\Delta }t)}-v_{\alpha (0)})\right).$$

From (2.9), the conditional characteristic function of

$$R_{\mathrm{\Delta }t}^{\alpha }:=\log \left(\frac{S_{\mathrm{\Delta }t}^{\alpha }}{S_0}\right)={\tilde{X}}_{\mathrm{\Delta }t}^{\alpha }+\frac{\rho }{{\sigma }_v}(v_{\alpha (\mathrm{\Delta }t)}-v_{\alpha (0)}),$$

is recovered in closed form as

	$?[{\mathrm{e}}^{\mathrm{i}{R}_{\mathrm{\Delta }t}^{\alpha }\xi }\mid \alpha (0)=j,\alpha (\mathrm{\Delta }t)=k]$	$={ℳ}_{k,j}(\xi ;\mathrm{\Delta }t)\exp \left(\mathrm{i}\xi {\displaystyle \frac{\rho }{{\sigma }_v}}(v_k-v_j)\right)$
		$:={\widetilde{ℳ}}_{k,j}(\xi ;\mathrm{\Delta }t),$		(2.10)

which follows from conditional independence. We can view the CTMC–Heston model as both an approximation to the Heston model and a tractable model in its own right.

2.3 Generator construction

In order to evaluate the characteristic function above, a generator $Q$ describing the transitions of $\alpha (t)$ is required. Contrary to a standard RS model, where the generator must be specified as an input of the model (it drives the changes in regimes), we need to construct the generator $Q$ according to the CTMC approximation of the Heston volatility process.¹¹ 1 In this sense, the CTMC–Heston model can be thought of as a parsimonious RS model that mimics the dynamics of the Heston model. Rather than explicitly specifying the transition rates in each regime, we will determine them implicitly via the Heston model parameters. We follow the formula proposed in Lo and Skindilias (2014).

Given a grid of points $?=\{v_1,v_2,\mathrm{\dots },v_{m_0}\}$ with grid spacings $h_i=v_{i+1}-v_i$, and assuming that $v_{\alpha (t)}$ takes values on $?$, the elements $q_{ij}$ of the generator $Q$ for the CTMC approximation of the process $v_t$ defined in (2.4) read as follows:

(2.11)

with the notation $z^{\pm }=\max (\pm z,0)$. Further, to guarantee a well-defined probability matrix, the following condition must be satisfied:

(2.12)

So far, we have assumed a vector of grid points, $?$, as a given. In Section 4.1, more insight into its optimal construction is provided. However, now that the model has been specified, the next step is to estimate the transition densities using Fourier inversion.

3 Wavelet-based Fourier inversion

This section presents a highly effective approach to estimating probability densities in terms of Shannon wavelets, which uses knowledge of the ChF. Compared with local bases, such as the B-splines considered in Ortiz-Gracia and Oosterlee (2013) and Kirkby (2015), Shannon wavelets result in an exponentially convergent representation of the density in terms of the scale parameter defined below. They also alleviate truncation issues that are a known limitation of existing methods, such as the “state-of-the-art” COS method (Fang and Oosterlee 2008) or the recently introduced PROJ method (Kirkby 2017). Further, as we will describe in Section 5.2.1, the use of Shannon wavelets presents an additional advantage when pricing exotic options, since it enables a very fast and accurate approximation of a computationally expensive integral, which otherwise (by employing other bases) has to be solved via quadrature numerical techniques.

Here, we give a brief review of the SWIFT method (Ortiz-Gracia and Oosterlee 2016), an efficient Fourier inversion technique based on Shannon wavelets. For completeness, we devote a section to the basic theory on Shannon wavelets.

3.1 Multiresolution analysis and Shannon wavelets

Consider the space of square-integrable functions denoted by $L^2(ℝ)$. A general structure for wavelets in $L^2(ℝ)$ is called a multiresolution analysis. We start with a family of closed nested subspaces in $L^2(ℝ)$,

	$$\mathrm{\cdots }\subset {?}_{-2}\subset {?}_{-1}\subset {?}_0\subset {?}_1\subset {?}_2\subset \mathrm{\cdots },$$
	$$\bigcap _{m\in \mathbb{Z}}{?}_m=\{0\},\overline{\bigcup _{m\in \mathbb{Z}}{?}_m}=L^2(ℝ),$$

where

$$f(x)\in {?}_m⟺f(2x)\in {?}_{m+1}.$$

If these conditions are met, then there exists a function $\phi \in {?}_0$ that generates an orthonormal basis, denoted by

$${\{{\phi }_{m,l}\}}_{l\in \mathbb{Z}},$$

for each ${?}_m$ subspace, ${\phi }_{m,l}(x)=2^{m/2}\phi (2^{m}x-l)$. The function $\phi$ is usually referred to as the scaling function or father wavelet.

For any $f\in L^2(ℝ)$, a projection map of $L^2(ℝ)$ onto ${?}_m$, denoted by ${?}_m:L^2(ℝ)\to {?}_m$, is defined by means of

$${?}_{m}f(x)=\sum _{l\in \mathbb{Z}}{c}_{m,l}{\phi }_{m,l}(x),\text{with }{c}_{m,l}=⟨f,{\phi }_{m,l}⟩,$$

(3.1)

where $⟨f,g⟩:={\int }_{ℝ}f(x)\overline{g}(x)dx$ denotes the inner product in $L^2(ℝ)$, with $\overline{g}$ being the complex conjugate of $g$. When we consider higher $m$ values (ie, when more terms are used), the truncated series representation of the function $f$ improves. As opposed to Fourier series, a key fact regarding the use of wavelets is that wavelets can be moved (by means of the $l$ value), stretched or compressed (by means of the $m$ value) to accurately represent the local properties of a function. A basic reference work on wavelets is Daubechies (1992).

In this paper, we employ Shannon wavelets (Cattani 2008). A set of Shannon scaling functions ${\phi }_{m,l}$ in the subspace ${?}_m$ is defined as

$${\phi }_{m,l}(x)=2^{m/2}\frac{\sin (\pi (2^{m}x-l))}{\pi (2^{m}x-l)}=2^{m/2}\phi (2^{m}x-l),l\in \mathbb{Z},$$

(3.2)

where

is the scaling function, or cardinal sine function.

3.2 SWIFT method

Given a function $f\in L^2(ℝ)$, we consider its expansion in terms of Shannon scaling functions at the level of resolution $m$. In the problem addressed here, $f$ will always be a probability density function that will be recovered from its Fourier transform, ie, the ChF, denoted here by $\varphi$ (see, for example, the ChF in (2.10)).²² 2 Here, we consider the strict definition of the Fourier transform – without sign reversal in the complex exponential – that is usually employed in the definition of the ChF.

Following wavelets theory, a density function $f$ can be approximated at resolution level $m$ by

$$f(x)\approx {?}_{m}f(x)=\sum _{k\in \mathbb{Z}}{c}_{m,l}{\phi }_{m,l}(x),$$

(3.3)

where ${?}_{m}f$ converges to $f$ in $L^2(ℝ)$, ie, ${\parallel f-{?}_{m}f\parallel }_2\to 0$ when $m\to +\mathrm{\infty }$. Here, the coefficients $c_{m,l}$ and the scaling functions ${\phi }_{m,l}$ are defined in (3.1) and (3.2), respectively. In Maree et al (2017), it was shown that

$${\parallel f-{?}_{m}f(x)\parallel }_{\mathrm{\infty }}\le H(2^m\pi ;f),H(\gamma ;f):=\frac{1}{2\pi }{\int }_{|\xi |>\gamma }|\varphi (\xi )|d\xi ,$$

(3.4)

where, as mentioned, $\varphi$ is the ChF of $f$.

For many models in finance, the rapidly decaying tails of the ChF yield exponential uniform convergence of ${?}_{m}f(x)$ to $f$.

The infinite series in (3.3) can be well approximated (see Ortiz-Gracia and Oosterlee (2016, Lemma 1) for details) by a finite summation,

$$f_m(x)\approx \sum _{l=l_1}^{l_2}{c}_{m,l}{\phi }_{m,l}(x),$$

(3.5)

for certain accurately chosen values $l_1$ and $l_2$.

The next step is the computation of the density coefficients $c_{m,l}$ in (3.5). Recalling (3.1),

$$c_{m,l}=⟨f,{\phi }_{m,l}⟩={\int }_{ℝ}f(x){\overline{\phi }}_{m,l}(x)dx=2^{m/2}{\int }_{ℝ}f(x)\phi (2^{m}x-l)dx,$$

(3.6)

or, equivalently, by Parseval’s identity,

$$c_{m,l}=\frac{1}{2\pi }⟨\varphi ,{\widehat{\phi }}_{m,l}⟩=\frac{1}{2\pi }{\int }_{ℝ}\varphi (x){\widehat{\phi }}_{m,l}(x)dx,$$

(3.7)

where ${\widehat{\phi }}_{m,l}(\xi )$ denotes the Fourier transform of the scaling function ${\phi }_{m,l}(x)$.

Both suggested expressions can be employed to come up with efficient procedures to estimate the coefficients $c_{m,l}$. Here, we will present the final expressions only. We refer the reader to Ortiz-Gracia and Oosterlee (2016) for a more detailed description.

First, by applying the classical Vieta formula truncated with $J$ terms to (3.6), the cosine product-to-sum identity and the definition of the ChF, and after some algebraic manipulations, the coefficients $c_{m,l}$ can be accurately approximated by

$$c_{m,l}\approx {\overline{c}}_{m,l}:=\frac{2^{m/2}}{2^{J-1}}\sum _{\lambda =1}^{2^{J-1}}\mathrm{Re}\left[\varphi \left(\frac{(2\lambda -1)\pi 2^m}{2^J}\right){\mathrm{e}}^{\mathrm{i}l\pi (2\lambda -1)/2^J}\right].$$

(3.8)

On the other hand, if we start from (3.7) and note that ${\widehat{\phi }}_{m,l}(\xi )$ (ie, the Fourier transform of ${\phi }_{m,l}(x)$) reads

$${\widehat{\phi }}_{m,l}(\xi )=\frac{{\mathrm{e}}^{\mathrm{i}(l/2^m)\xi }}{2^{m/2}}\mathrm{rect}\left(\frac{\xi }{2^{m+1}\pi }\right),$$

where $\mathrm{rect}$ represents the rectangular function, and if we apply the change of variables $\omega =\xi /2^{m+1}\pi$, the $c_{m,l}$ can be computed by

$$c_{m,l}=2^{m/2}{\int }_{-1/2}^{1/2}\mathrm{Re}(\varphi (2^{m+1}\pi \omega ){\mathrm{e}}^{\mathrm{i2}\pi l\omega })d\omega .$$

(3.9)

Remark 3.1.

Equations (3.8) and (3.9) can be conveniently rearranged, after applying the trapezoidal rule in the second case, in order to compute the coefficients via the fast Fourier transform (FFT) algorithm (see Ortiz-Gracia and Oosterlee 2016), with computational complexity $?(\log (2)J{2}^{J+1})$.

3.3 Optimal scale $m$ and series truncation bounds $l_1$ and $l_2$, and $J$

The quality of the approximation provided by the SWIFT method is mainly affected by the scale $m$ and $l_1$ and $l_2$ in (3.5). In order to optimally determine these parameters, we follow the analysis presented in Colldeforns-Papiol and Ortiz-Gracia (2018) and Maree et al (2017).

Let us start with the choice of $m$. From (3.4), the error in the projection approximation of function $f$ in (3.3) is bounded. As the ChF, $\varphi$, is assumed to be known, we can compute $m$ given a prescribed tolerance ${\epsilon }_m$. The integral is typically not known in closed form, so it can be approximated by numerical techniques. Applying a very simple quadrature rule, an error bound can be obtained:

$$\frac{1}{2\pi }(|\varphi (-2^m\pi )|+|\varphi (2^m\pi )|)\le {\epsilon }_m.$$

This approximation is considered sufficient for our purposes. More involved numerical quadratures have been tested, but the observed differences are negligible.

Regarding truncation limits $l_1$ and $l_2$, these can be computed based on the predefined interval $[a,b]$, where the density mass in concentrated. Thus,

$$l_1:=\lfloor 2^{m}a\rfloor \mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{l}_2:=\lceil 2^{m}b\rceil ,$$

where $m$ is the scale of approximation. Therefore, we must choose the interval limits, $a$ and $b$, in such a way that the loss of density mass is minimized (see Appendix B.1).

Further, in (3.8), the parameter $J$ needs to be selected. It can be chosen to be constant based on previously determined quantities.³³ 3 It could be selected as a function of $l$. By doing so, we can benefit from the use of the FFT algorithm (see Colldeforns-Papiol and Ortiz-Gracia (2018) for more details). Thus,

where $M_{m,l}:=\max (|2^m?-l|,|2^m?+l|)$ and $?:=\max (|a|,|b|)$.

3.4 SWIFT for European options under CTMC–Heston

The starting point to solve an option pricing problem is the risk-neutral option valuation formula, where the discounted value of an option at future time $t$ and given the current state $x$, $V(x,t)$, is an expectation under the risk-neutral pricing measure, $\mathbb{Q}$, ie,

$$V(x,t)={\mathrm{e}}^{-rt}?[h(y)\mid x]={\mathrm{e}}^{-rt}{\int }_{ℝ}h(y)f(y\mid x)dy,$$

(3.10)

with $r$ the risk-free rate, $f(y\mid x)$ the transitional density and $h(y)$ the payoff function. The density function $f$ is typically unknown, even for the case of European option pricing, where the value of the variable $y$ depends solely on the state at maturity $T$. In order to compute the value of the option $V$ in (3.10), we can replace the transitional probability density function $f$ by its approximation in (3.5). Accordingly, the SWIFT pricing formula reads

$$V(x,t)\approx {\mathrm{e}}^{-rt}\sum _{l=l_1}^{l_2}(c_{m,l}{\int }_a^{b}h(y){\phi }_{m,l}(y\mid x)dy)={\mathrm{e}}^{-rt}\sum _{l=l_1}^{l_2}{c}_{m,l}{B}_{m,l},$$

(3.11)

with the payoff coefficients $B_{m,l}$ defined as

$$B_{m,l}:={\int }_a^{b}h(y){\phi }_{m,l}(y\mid x)dy.$$

We have reduced the computation of the option price to the dot product between the so-called density coefficients, $c_{m,l}$, and the payoff coefficients, $B_{m,l}$. The solution of the integral in the computation of $B_{m,l}$ is analytically available for many situations in finance, including those studied here.

For the CTMC–Heston model, European options are simple to value using SWIFT. If we let $f_R^{j_0}(y;T)$ denote the density of log return $R_T^{\alpha }:=S_T^{\alpha }/S_0$, conditioned on $v_{\alpha (0)}=v_{j_0}$,

$${\varphi }_R^{j_0}(\xi ;T):=?[{\mathrm{e}}^{\mathrm{i}{R}_T^{\alpha }\xi }\mid \alpha (0)=j_0]={\mathrm{?}}^{\prime }\widetilde{ℳ}(\xi ;T){?}_{j_0},j_0\in \{1,\mathrm{\dots },m_0\},$$

where $\widetilde{ℳ}(\xi ;T)$ is defined in (2.10). The coefficients $c_{m,l}$ are then found using (3.8), with ${\varphi }_R^{j_0}(\xi ;T)$ in place of $\varphi (\xi )$. To assess the accuracy of this approximation, we require the following lemma.

Lemma 3.2.

For fixed $m_{\mathrm{0}}$, $\xi \mathrm{\in }\mathrm{R}$, and $\mathrm{\Delta }\mathit{}t\mathrm{>}\mathrm{0}$,

$$\underset{1\le k,j\le m_0}{sup}|{\widetilde{ℳ}}_{k,j}(\xi ;\mathrm{\Delta }t)|=|{ℳ}_{k,j}(\xi ;\mathrm{\Delta }t)|\le \exp (-{\xi }^2\frac{1}{2}\mathrm{\Delta }t(1-{\rho }^2)v_1).$$

(3.12)

Moreover, ${\varphi }_R^j\mathrm{(}\xi \mathrm{;}\mathrm{\Delta }t\mathrm{)}\mathrm{=}\mathrm{E}\mathrm{[}{\mathrm{e}}^{\mathrm{i}\mathit{}{R}_{\mathrm{\Delta }\mathit{}t}^{\alpha }\mathit{}\xi }\mathrm{\mid }\alpha \mathrm{(}\mathrm{0}\mathrm{)}\mathrm{=}j\mathrm{]}$ satisfies the same bound for any $j\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{\dots }\mathrm{,}{m}_{\mathrm{0}}$.

Proof.

See Appendix A (available online). ∎

Let $V(x,\mathrm{\Delta }t)$ denote the true value of a claim with $\mathrm{\Delta }t$ until maturity, $V(x,\mathrm{\Delta }t):={\mathrm{e}}^{-r\mathrm{\Delta }t}?[h(R_{\mathrm{\Delta }t})]$, where $R_t:=\log (S_t/S_0)$, and $V_{m_0}^m(x,\mathrm{\Delta }t)$ denotes the approximation based on $R_{\mathrm{\Delta }t}^{\alpha }$ (defined above), with $m_0$ CTMC states and resolution level $m$.

Proposition 3.3.

Let $h\mathit{}\mathrm{(}y\mathrm{)}\mathrm{=}h\mathit{}\mathrm{(}y\mathrm{|}{S}_{\mathrm{0}}\mathrm{)}$ denote a payoff in the space of $R_{\mathrm{\Delta }\mathit{}t}$. Choose a parameterization of $Q$ from (2.11) that has an algebraic convergence order $\mathrm{O}\mathit{}\mathrm{(}{m}_{\mathrm{0}}^{\mathrm{-}\delta }\mathrm{)}$, eg, a uniform grid approximation (for which $\delta \mathrm{=}\mathrm{2}$). Then, for $\mathrm{[}a\mathrm{,}b\mathrm{]}$ sufficiently wide,

	${\mathrm{e}}^{r\mathrm{\Delta }t}|V(x,\mathrm{\Delta }t)-V_{m_0}^m(x,\mathrm{\Delta }t)|$
	$\mathrm{\hspace{1em}\hspace{1em}}\mathit{\hspace{1em}\hspace{1em}}\le \tau ([a,b])+\parallel h{\parallel }_{\mathrm{\infty }}^{[a,b]}?(2^m\tau ([a,b])+{\mathrm{e}}^{-2^{2m}})+?(m_0^{-\delta }),$

where

$$V_{m_0}^m(x,\mathrm{\Delta }t)={\mathrm{e}}^{-r\mathrm{\Delta }t}\sum _{l=l_1}^{l_2}{\overline{c}}_{m,l}^{j_0}{\int }_a^{b}h(y){\phi }_{m,l}(y\mid x)dy,$$

and $\tau \mathit{}\mathrm{(}\mathrm{[}a\mathrm{,}b\mathrm{]}\mathrm{)}\mathrm{:=}{\mathrm{\int }}_{{\mathrm{[}a\mathrm{,}b\mathrm{]}}^c}{f}_R^{j_{\mathrm{0}}}\mathit{}\mathrm{(}y\mathrm{\mid }x\mathrm{)}\mathit{}h\mathit{}\mathrm{(}x\mathrm{)}\mathit{}\mathrm{d}y$ is the truncation error.

Proof.

See Appendix A (available online). ∎

Remark 3.4.

As Shannon wavelets present an oscillating nature, the approximated density could potentially reach negative values, leading to an arbitrage opportunity. However, with a sufficient number of expansion terms, this effect is negligible, since the densities typically appearing in finance are very smooth. Further, as shown, the SWIFT method provides a mechanism to control the error in the density tails, where this undesired behavior is more probable. All in all, we have not encountered any related issue in the conducted experiments.

4 Calibration of the CTMC–Heston model with SWIFT

The proposed approach combines CTMC approximation techniques with the SWIFT valuation method, which has several important advantages that make it well suited for a calibration procedure. Here, we will summarize some of them.

Error control.

This is probably the most relevant property within an optimization problem. Thanks to the use of Shannon wavelets, SWIFT establishes a bound in the error given any scale $m$ of approximation.

Robustness.

SWIFT provides mechanisms to determine all of the free parameters in the approximation made based on the scale $m$, which, as mentioned in the previous point, determines the committed error.

Performance efficiency.

As with other Fourier inversion techniques, SWIFT is an extremely fast algorithm, allowing FFTs, vectorized operations or even parallel computing features.

Accuracy.

Although an error bound is provided, SWIFT has demonstrated a very high precision in most situations, far below the predicted error bound and (at least) comparable to the state-of-the-art methodologies.

Remark 4.1.

The properties mentioned above are particularly desirable in the calibration process, providing high-quality estimations and reducing the chances of any possible malfunctioning or divergence in the optimization procedure.

Here, we focus on the calibration aspects related to the application of the CTMC approximation, that is, we assume the pure Heston model parameters to be given. Thus, there are two major components affecting the quality of the CTMC–Heston version of the volatility process: the length and the distribution of points in the Markov chain, ie, $m_0$ and the locations $v_1,v_2,\mathrm{\dots },v_{m_0}$, respectively. We will study the “optimal” selection of these two components. Let us start with the parameter $m_0$. It has been proved in Mijatović and Pistorius (2013) that $v_{\alpha (t)}$ weakly converges to its continuous counterpart $v_t$, so⁴⁴ 4 This assumes that the domain $[v_1,v_{m_0}]$ eventually “covers” the domain of $v_t$. See Mijatović and Pistorius (2013) for some technical details.

$$v_{\alpha (t)}\to v_t\mathit{\hspace{1em}}\text{as }\max \{|h_i|:i=1,\mathrm{\dots },m_0-1\}\to 0,$$

(4.1)

recalling $h_i=v_{i+1}-v_i$. An immediate conclusion of this result is that increasing the width of the chain does not necessarily ensure convergence, but it is also a requirement that the gap between grid points approaches zero. So, as expected, larger $m_0$ values cannot be analyzed without studying the distribution of the newly added points, which should be carefully selected to ensure fast convergence in the approximations. Moreover, we must comply with the condition established in (2.12).

4.1 Grid selection

One of the key aspects in designing the CTMC–Heston model is specifying the variance state space (grid), $?$. Our goal is to form a model that parsimoniously resembles Heston, and grid design has a strong influence here. We will analyze several conceptually different approaches that are available in the literature, which can be classified into two main categories based on state domain truncation and distribution inversion methodologies. In the context of barrier option pricing, several designs have been studied in Cui and Taylor (2019). Alternative designs, which modify the grid to include barriers and strikes for pricing single options, have been considered in Li and Zhang (2017, 2019) and Zhang and Li (2019). As our goal is model calibration across a spectrum of strikes, the grid design should accommodate the full range of the market.

The first step in designing a grid (or state space) for the process $v_{\alpha (t)}$ is to determine its boundaries, $v_1$ and $v_{m_0}$, where $m_0$ is the number of variance states. Let us therefore define

$$\overline{\mu }(t):=?[v_t\mid v_0]\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}\overline{\sigma }(t):=\sqrt{\mathrm{var}[v_t\mid v_0]},$$

which are known for the Heston model:

$$\overline{\mu }(t)={\mathrm{e}}^{-\eta t}{v}_0+\theta (1-{\mathrm{e}}^{-\eta t}),{\overline{\sigma }}^2(t):=\frac{{\sigma }_v^2}{\eta }{v}_0({\mathrm{e}}^{-\eta t}-{\mathrm{e}}^{-2\eta t})+\frac{\theta {\sigma }_v^2}{2\eta }{(1-{\mathrm{e}}^{-\eta t})}^2.$$

Next, we can choose the bounds of the truncated space as $v_1=\overline{\mu }-\gamma \overline{\sigma }$ and $v_{m_0}=\overline{\mu }+\gamma \overline{\sigma }$ for some prescribed parameter $\gamma$, and we can further restrict $v_1>0$. The time frame $t$ is selected depending on the financial product at hand. For European options, it is set to the maturity time, ie, $t=T$. For path-dependent options (as in Section 5), a choice providing acceptable results is $t=\frac{1}{2}T$.

As the most simple approach, an equally spaced (uniform) grid of points is determined by

$$v_i=v_1+(i-1)\frac{v_{m_0}-v_1}{m_0-1},i=2,\mathrm{\dots },m_0-1.$$

More involved schemes can be introduced by considering nonuniform grids, with the aim of better capturing the relevant parts of the truncated state space and improving the convergence rate. Thus, the algorithm proposed by Tavella and Randall (2000) reads

$$v_i=v_0+\overline{\alpha }\sinh \left(c_2\frac{i}{m_0}+c_1\left(1-\frac{i}{m_0}\right)\right),i=2,\mathrm{\dots },m_0-1,$$

with

$$c_1=\text{arc}\sinh \left(\frac{v_1-v_0}{\overline{\alpha }}\right)\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{c}_2=\text{arc}\sinh \left(\frac{v_{m_0}-v_0}{\overline{\alpha }}\right),$$

where parameter $\overline{\alpha }$ controls the nonuniformity of the grid. Mijatović and Pistorius (2013) presented a slight modification of the Tavella and Randall algorithm based on the use of subgrids. The grid points are therefore chosen as follows:

Finally, the last approach takes advantage of all of the available information to select the grid of discrete points. Since the distribution of the volatility process in the Heston model is known ($v_t$ is driven by a CIR process), it seems natural to exploit this fact in order to optimally determine the location of the points. That is precisely the alternative suggested by Lo and Skindilias (2014), based on the minimization of the Kolmogorov–Smirnov distance between the original cumulative distribution function and its CTMC counterpart. By adapting their result to our problem, we find that the grid points $v_1,v_2,\mathrm{\dots },v_{m_0}$ can be chosen automatically as

$$v_i=F_{v_t}^{-1}\left(\frac{2i-1}{2m_0}\right),$$

where $F_{\{\cdot \}}$ is the cumulative distribution function of the given random variable, which, in our case, is

$$v_t\sim \frac{{\chi }^2(4\eta \theta /{\sigma }_v^2,v_0{c}_t{\mathrm{e}}^{-\eta t})}{c_t}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{c}_t=\frac{4\eta }{{\sigma }_v^2(1-{\mathrm{e}}^{-\eta t})},$$

where ${\chi }^2(d,c)$ is the noncentral chi-square distribution with $d$ degrees of freedom and a noncentrally parameter $c$.

In Figure 1, given a model parameter setting, the grid point locations ($m_0=30$) obtained by each of the approaches presented above are depicted. We can clearly identify the different distribution schemes among them.

Figure 1: Grid locations. Maturity: $T=\text{1}$. Heston parameters: $v_{\text{0}}=\text{0.4}$, $\eta =\text{3}$, $\theta =\text{0.4}$, ${\sigma }_v=\text{0.5}$, $\rho =-\text{0.1}$. Additional parameters: $\gamma =\text{4.5}$ and $\overline{\alpha }=(v_{m_{\text{0}}}-v_{\text{1}})/\text{5}$. (a) Uniform. (b) Mijatovic–Pistorius. (c) Tavella–Randall. (d) Lo–Skindilias.

4.2 CTMC: numerical study

Next, we perform a numerical study of the four proposed grid selection approaches in the context of the CTMC–Heston model. The experiments will consist of assessing the error in pricing European options, for which a reference price is easily obtained from the Heston characteristic function (available in closed form) by (3.11). In order to cover a complete spectrum of possible situations, we choose two Heston parameter settings that are directly related to the volatility process: one corresponding to a regular scenario (low initial volatility, low long-term volatility and moderate vol-vol), and one corresponding to a stressed scenario (high initial volatility, high long-term volatility and mid-high vol-vol). The specific parameter values are presented in Table 1.

Table 1: Data sets.

	Scenario	${?}_{\text{?}}$	$?$	$?$	${?}_{?}$	$?$
Set I	Regular market	0.03	3.0	0.04	0.25	$-$0.7
Set II	Stressed market	0.4	3.0	0.4	0.5	$-$0.1

4.2.1 Convergence in $m_0$

First, the convergence with respect to the number of grid points, $m_0$, is studied. We consider European option prices as the observed variable, with $V$ the reference value and $V_{m_0}$ the CTMC–Heston estimation, measuring the differences by means of the relative error. As mentioned, the convergence in the CTMC approximation of the variance process $v_t$ is ensured when the condition in (4.1) is fulfilled. Each grid selection methodology presented here satisfies this requirement, so a similar convergence behavior should be translated to the option valuation results. Indeed, in Figure 2 we can see the obtained convergence curves. In particular, Proposition 3.3 predicts a polynomial convergence rate of $?(m_0^{-\delta })$, where $\delta$ depends on the grid design. Here, we observe second-order convergence ($\delta =2$) for each grid design considered, with the exception of the Kolmogorov–Smirnov minimizer (Lo–Skindilias) provided in Lo and Skindilias (2014), for test set I.

Figure 2: Convergence in $m_{\text{0}}$. Market parameters: put option, $S_{\text{0}}=\text{100}$, $K=\text{100}$, $r=\text{0.05}$ and $T=\text{1}$. Additional parameters (truncated approaches): $\gamma =\text{10}$ and $\overline{\alpha }=(v_{m_{\text{0}}}-v_{\text{1}})/\text{5}$. (a) Set I. (b) Set II.

4.2.2 Influence of the model parameters

Now a more sophisticated analysis is carried out by systematically varying the model-related parameters affecting the variance process, ie, $v_0,\eta ,\theta ,{\sigma }_v,\rho$ and $T$, and studying the effectiveness of each grid selection methodology among several configurations. The length of the Markov chain is fixed and set to $m_0=100$, which, according to Figure 2, provides a relatively high accuracy (${10}^{-4}-{10}^{-7}$). The results are presented in Figures 3 and 4 for set I and set II, respectively.

Figure 3: Impact of the model parameters. Market parameters: put option, $S_{\text{0}}=\text{100}$, $K=\text{100}$, $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set I. Additional parameters (truncated approaches): $\gamma =\text{10}$ and $\overline{\alpha }=(v_{m_{\text{0}}}-v_{\text{1}})/\text{5}$.

Figure 4: Impact of the model parameters. Market parameters: put option, $S_{\text{0}}=\text{100}$, $K=\text{100}$, $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set II. Additional parameters (truncated approaches): $\gamma =\text{10}$ and $\overline{\alpha }=(v_{m_{\text{0}}}-v_{\text{1}})/\text{5}$.

4.2.3 Calibration experiment

Finally, a real calibration test is carried out on market data. We consider European call options on Microsoft and Google stocks quoted in January 2019. Because it is common in practice, we adjust the model parameters to fit the implied volatility observed in the market (calibration in volatilities). First, we calibrate the classical Heston model, employing the SWIFT method (together with the ChF of the Heston model for European pricing) to get the prices required in the definition of the objective function. The prices are inverted to obtain the corresponding implied volatilities. Then, this objective function is used in an optimization procedure, which aims to determine the Heston parameters minimizing the differences between market and model values. The mean square error is chosen as our error measure, while the interior-point algorithm is considered as our optimization method.

Once the Heston parameters are determined, we aim to select $m_0$ in such a way that the CTMC–Heston implied volatility approximates the Heston implied volatility up to some prescribed tolerance. Again, we analyze the four methodologies presented in Section 4.1 for variance process grid selection.

Figures 5 and 6 show the calibration outcomes for the options on Microsoft and Google, respectively. In the captions, the observed market data and the calibrated Heston parameters are presented in each case. The thick black curve corresponds to the calibrated Heston implied volatility, which is considered here as our reference. The bulleted curves correspond to CTMC–Heston estimations with an increasing value of $m_0$. We stop the procedure when the error between the Heston curve and the CTMC–Heston curve with a certain $m_0$ prescribes a given tolerance. The tolerance is set to ${10}^{-3}$, and the fit between curves is measured by the mean relative error along the strikes, $K$. We also employ a limit in the maximum chain length, set to $m_0=200$.

Figure 5: Microsoft calibration curves for varying $m_{\text{0}}$. Market parameters: call options, $S_{\text{0}}=\text{105.36}$, $K=\{\text{65},\text{70},\mathrm{\dots },\text{150},\text{155}\}$ (21 strikes), $r=\text{0.0246}$, $q=\text{0}$ and $T=\text{0.4986}$. Heston parameters: $v_{\text{0}}=\text{0.0906},\eta =\text{0.8549},\theta =\text{0.1379},{\sigma }_v=\text{0.9976},\rho =-\text{0.6187}$. (a) Uniform. (b) Mijatović–Pistorius. (c) Tavella–Randall. (d) Lo–Skindilias.

Figure 6: Google calibration curves for varying $m_{\text{0}}$. Market parameters: call options, $S_{\text{0}}=\text{1080.66}$, $K=\{\text{880},\text{890},\mathrm{\dots },\text{1290},\text{1295}\}$ (57 strikes), $r=\text{0.0249}$, $q=\text{0}$ and $T=\text{0.9972}$. Heston parameters: $v_{\text{0}}=\text{0.1482},\eta =\text{0.7752},\theta =\text{0.0722},{\sigma }_v=\text{0.9278},\rho =-\text{0.5444}$. (a) Uniform. (b) Mijatović–Pistorius. (c) Tavella–Randall. (d) Lo–Skindilias.

Several interesting lessons can be extracted from the previous experiments.

• •

  All of the approaches provide numerical convergence in $m_0$. The Mijatović–Pistourius and Lo–Skindilias schemes present a singularity in the error curve for set I that can be attributed to this particular parameter configuration, since this issue does not appear in the test with set II.

• •

  In all of the cases, the error decays very fast at the beginning, with smaller $m_0$, and smooths for bigger $m_0$, suggesting a damping effect.

• •

  The grid distribution proposed by Lo and Skindilias, despite being more intuitive, provides, in general, poorer estimations. In the first experiment, we can even see that increasing the number of states hardly improves the convergence. According to the parameter experiments, only when the variance of the process is relatively low are the estimations comparable to those of the other methodologies. This fact suggests that the state space is not sufficiently covered.

• •

  Although the uniform approach performs surprisingly well in the test with synthetic parameters, the real calibration experiment (using market data) shows a pretty inaccurate estimation for options far from an at-the-money strike, as the method is not able to even achieve the required tolerance, showing as well a rotation over the at-the-money value.

• •

  The two methodologies relying on moment-based domain truncation and nonuniform distribution perform similarly. It is worth noting that the approach of Mijatović and Pistorius explodes when the initial and long-term volatilities differ greatly from one other.

• •

  While the performance varies across the different parameter settings, these experiments indicate that there is some advantage to the nonuniform method of Tavella and Randall, which happens to be the most robust and precise choice in general.

• •

  The method of Mijatović and Pistorius turns out to be highly precise for the real market data studied here, requiring very few Markov states to achieve the prescribed accuracy.

• •

  By focusing on the correlation parameter, $\rho$, in the second test, we observe that the error tends to be at its minimum close to the no-correlation point ($\rho =0$), and it degrades when $\rho$ ventures far from zero.

Remark 4.2.

In addition to these numerical observations, we observe from the calibration of Google and Microsoft data in Figures 5 and 6 that the CTMC–Heston model with finitely many states is able to reliably reproduce market quotes; this is consistent with the implied volatility curves that would be generated by the Heston model. Yet, the CTMC–Heston model is a much more tractable model when it comes to pricing exotic contracts. It therefore offers a promising alternative to the Heston model that is useful in practice.

5 Application: pricing exotic options under CTMC–Heston model

Once a model has been calibrated (typically to European options), it is commonly employed to price more involved products with early exercise features, path-dependent payoffs, etc, under the calibrated model. In this work, we shall present several applications under the CTMC–Heston model which address well-known complex valuation problems that appear in finance.

The methodology presented encompasses a wide range of exotic contract types, which we define implicitly in terms of a generic recursion that they satisfy. Consider a fixed horizon $[0,T]$ on which there are defined $N+1$ monitoring dates, . We define the log return $R_n$ between monitoring dates as

$$R_n:=\log \left(\frac{S_n}{S_{n-1}}\right),S_n:=S(t_n),n=1,\mathrm{\dots },N,$$

where $t_{n+1}-t_n$ is not required to be uniform. The contracts of interest satisfy a very general sequence of equations

(5.1)

where $h,g:ℝ\to ℝ$ are continuous functions, ${\{w_n\}}_{n=1}^N$ is a set of weights and ${\{{\varrho }_n\}}_{n=1}^N$ is a set of shift parameters. This includes all contracts of the form

$$G(\sum _{n=1}^N{w}_{n}h(R_n);\mathrm{\Pi }),$$

where $\mathrm{\Pi }$ is a set of generic contract parameters. For example, when $h(R_n)=\exp (R_n)-1$ or $h(R_n)=R_n$, we have the two common varieties of realized variance derivatives as well as cliquets/ratchets/equity indexed annuities. This formulation also covers more general iterated compositions of the form

$$G(\underset{n=1}{\overset{N}{\otimes }}{h}_n(R_n);\mathrm{\Pi }),$$

where ${\otimes }_{n=1}^N$ denotes the composition $h_1(R_1)\circ h_2(R_2)\circ \mathrm{\cdots }\circ h_N(R_N)$, where summation is, of course, a special case. Prominent examples of contracts that fall within this framework include the following.

Realized variance swaps and options:

$$A_N=\frac{1}{T}\sum _{n=1}^N{(R_n)}^2\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{A}_N=\frac{1}{T}\sum _{n=1}^N{(\exp (R_n)-1)}^2,$$

(5.2)

with $G(A_N):=A_N-K$ for a swap, and $G(A_N):={(A_N-K)}^{+}$ for a call option.

Cliquets:

with local floor and cap $F$ and $G$, respectively, and global floor and cap $F_g$ and $G_g$,

$$A_N=\sum _{n=1}^N\max (F,\min (C,\exp (R_n)-1)),$$

with $G(A_N)=K\min (C_g,\max (F_g,A_N))$.

Asian options:

	$A_N$	$:={\displaystyle \frac{1}{N+1}}{\displaystyle \sum _{n=0}^N}{w}_n{S}_n$
		$={\displaystyle \frac{S_0}{N+1}}(w_0+{\mathrm{e}}^{R_1}(w_1+{\mathrm{e}}^{R_2}(\mathrm{\cdots }{\mathrm{e}}^{R_{N-1}}(w_{N-1}+w_N{\mathrm{e}}^{R_N})))),$

where $G(A_N):={(A_N-K)}^{+}$ for a call option.

In the remainder of this work, we develop a highly efficient approach for tackling these problems in a unified way. The only real difference in the procedure between the various contracts is the specification of the terminal payoff, which will also be handled generically using the formula provided in (5.7).

5.1 Conditional ChF recursion for exotics

The pricing methodology presented in this work relies on obtaining the ChF of the recursion introduced in (5.1) when the underlying dynamics follow an RS model. Let $Y_N^{\alpha }$ denote the recursion with regime transitions driven by the process $\alpha$. For ease of exposition, we will focus on the case of unit weights, $w_n\equiv 1$,

$Y_1^{\alpha }:=h(R_N^{\alpha }),Y_n^{\alpha }:=h(R_{N-(n-1)}^{\alpha })+g(Y_{n-1}^{\alpha }),n=2,\mathrm{\dots },N,$

(5.3)

and uniform monitoring, where $t_{n+1}-t_n=\mathrm{\Delta }t:=T/N$, although the method is easily applied in the general case. Therefore, the increments are identically and independently distributed, ie, $R_n^{\alpha }\stackrel{?}{=}{R}_{\mathrm{\Delta }t}^{\alpha }$. For a large class of models, this formulation can be recast as a recursion in the frequency domain, where $h(R_{N-(n-1)}^{\alpha })$ and $g(Y_{n-1}^{\alpha })$ are (conditionally) independent, from which the ChFs factorize.

Conditional on ${\alpha }_{N-1}:=\alpha (t_{N-1})=j$, we define the ChF for $\xi \in ℝ$ and $j\in ?:=\{1,\mathrm{\dots },m_0\}$ as

	${\varphi }_{Y_1^{\alpha }}^j(\xi ):={?}_{N-1}[{\mathrm{e}}^{\mathrm{i}\xi Y_1^{\alpha }}\mid {\alpha }_{N-1}=j]$	$={?}_{N-1}[{\mathrm{e}}^{\mathrm{i}\xi h(R_N^{\alpha })}\mid {\alpha }_{N-1}=j]$
		$=?[{\mathrm{e}}^{\mathrm{i}\xi h(R_{\mathrm{\Delta }t}^{\alpha })}\mid \alpha (0)=j],$

where ${?}_n[\cdot ]:=?[\cdot \mid ℱ(t_n)]$; here, $ℱ(t)$, $0\le t\le T$, is the filtration generated by $S(t)$.

Similarly, for $n=2,\mathrm{\dots },N$, we have

	${\varphi }_{Y_n^{\alpha }}^j(\xi )$	$:={?}_{N-n}[{\mathrm{e}}^{\mathrm{i}\xi Y_n^{\alpha }}\mid {\alpha }_{N-n}=j]$
		$\mathrm{\hspace{0.17em}}={?}_{N-n}[\exp (\mathrm{i}\xi h(R_{N-n+1}^{\alpha })+\mathrm{i}\xi g(Y_{n-1}^{\alpha }))\mid {\alpha }_{N-n}=j],j\in ?,$

which corresponds to the density $f_{Y_n^{\alpha }}^j$.

To simplify our notation, we introduce the random variables $Z_{n-1}:=g(Y_{n-1}^{\alpha })$, $n=2,\mathrm{\dots },N$, and the corresponding conditional ChF

$${\varphi }_{Z_{n-1}}^j(\xi ):={?}_{N-n}[{\mathrm{e}}^{\mathrm{i}\xi Z_{n-1}}\mid {\alpha }_{N-(n-1)}=j],j\in ?.$$

Proposition 5.1.

Let $Y_N^{\alpha }$ be the solution to (5.3), and define for $\xi \mathrm{\in }\mathrm{R}$

$${\mathrm{\Phi }}_{j,k}(\xi ):=P_{j,k}^{\mathrm{\Delta }t}?[\exp (\mathrm{i}\xi h(R_{\mathrm{\Delta }t}^{\alpha }))\mid \alpha (0)=j,\alpha (\mathrm{\Delta }t)=k],j,k\in ?,$$

(5.4)

where $P_{j\mathit{}k}^{\mathrm{\Delta }\mathit{}t}\mathrm{=}\mathrm{P}\mathrm{[}\alpha \mathrm{(}t\mathrm{+}\mathrm{\Delta }t\mathrm{)}\mathrm{=}k\mathrm{\mid }\alpha \mathrm{(}t\mathrm{)}\mathrm{=}j\mathrm{]}$. The ChF ${\varphi }_{Y_N^{\alpha }}^j\mathit{}\mathrm{(}\xi \mathrm{)}$ for $j\mathrm{\in }\mathrm{J}$ is the solution to the following recursion. For $n\mathrm{=}\mathrm{1}$,

${\varphi }_{Y_1^{\alpha }}^j(\xi )$

$={\displaystyle \sum _{k=1,\mathrm{\dots },m_0}}{\mathrm{\Phi }}_{j,k}(\xi ),j\in ?,$

and for $n\mathrm{=}\mathrm{2}\mathrm{,}\mathrm{\dots }\mathrm{,}N$,

	${\varphi }_{Z_{n-1}}^j(\xi )$	$={\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi g(y))f_{Y_{n-1}^{\alpha }}^j(y)dy,j\in ?,$
	${\varphi }_{Y_n^{\alpha }}^j(\xi )$	$={\displaystyle \sum _{k=1,\mathrm{\dots },m_0}}{\varphi }_{Z_{n-1}}^k(\xi ){\mathrm{\Phi }}_{j,k}(\xi ),j\in ?.$

In particular, when $h\mathit{}\mathrm{(}x\mathrm{)}\mathrm{=}x$ (eg, for Asian options), we have ${\mathrm{\Phi }}_{j\mathrm{,}k}\mathit{}\mathrm{(}\xi \mathrm{)}\mathrm{\equiv }{\stackrel{\mathrm{~}}{\mathrm{M}}}_{k\mathrm{,}j}\mathit{}\mathrm{(}\xi \mathrm{;}\mathrm{\Delta }\mathit{}t\mathrm{)}$. Similarly, when $g\mathit{}\mathrm{(}x\mathrm{)}\mathrm{=}x$ (eg, for variance derivatives and cliquets), we have ${\varphi }_{Z_{n\mathrm{-}\mathrm{1}}}^j\mathit{}\mathrm{(}\xi \mathrm{)}\mathrm{\equiv }{\varphi }_{Y_{n\mathrm{-}\mathrm{1}}^{\alpha }}^j\mathit{}\mathrm{(}\xi \mathrm{)}$.

5.2 SWIFT recursive exotic valuation

Applying this framework in practice to price exotic options is accomplished by approximating the sequence of steps in Proposition 5.1 numerically. The integral terms in Proposition 5.1 can be well approximated using a Shannon wavelets density projection. From (5.4),

	${\mathrm{\Phi }}_{j,k}(\xi )$	$=P_{j,k}^{\mathrm{\Delta }t}?[\exp (\mathrm{i}\xi h(R_{\mathrm{\Delta }t}^{\alpha }))\mid \alpha (0)=j,\alpha (\mathrm{\Delta }t)=k]$
		$={\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi h(x))f_R^{j,k}(x)dx,j,k\in ?,$		(5.5)

where $f_R^{j,k}$ is the conditional density of $R_{\mathrm{\Delta }t}^{\alpha }$, which can be accurately recovered using SWIFT by Fourier inversion from its ChF.

Thus, the SWIFT projected density described in Section 3.2 is given by

$$f_R^{j,k}(x)\approx {?}_m{f}_R^{j,k}(x)=\sum _{l\in \mathbb{Z}}{c}_{m,l}^{j,k}{\phi }_{m,l}(x),c_{m,l}^{j,k}:=⟨f_R^{j,k},{\phi }_{m,l}⟩,$$

where $f_R^{j,k}(x)$ is defined by the ChF in (2.10), ${\widetilde{ℳ}}_{k,j}(\xi ;\mathrm{\Delta }t)$. After truncating the series, we have

$$f_R^{j,k}(x)\approx {\tilde{f}}_R^{j,k}(x):=\sum _{l_1\le l\le l_2}{c}_{m,l}^{j,k}{\phi }_{m,l}(x).$$

5.2.1 Approximating the integral ${\mathrm{\Phi }}_{j,k}(\xi )$

One of the crucial benefits of the SWIFT framework is that it greatly simplifies the computation of integrals that do not have closed-form expressions, such as those arising in (5.5), thanks to the Shannon wavelet basis. This fact has already been highlighted in Leitao et al (2018) for Asian options under Lévy processes, and it applies in a similar manner to the problem addressed in this work. The computational cost reduction is, however, significantly greater in the present application since the integral ${\mathrm{\Phi }}_{j,k}(\xi )$ needs to be computed for each state of the CTMC approximation, $j=1,\mathrm{\dots },m_0$.

Theorem 5.2 (Theorem 1.3.2 of Stenger (2010)).

Let $\theta$ be defined on $\mathrm{R}$ with its Fourier transform $\widehat{\theta }$. Then,

$$\left|\frac{1}{\delta }{\int }_{ℝ}\theta (x)S_{j,\delta }(x)dx-\theta (j\delta )\right|\le \frac{1}{2\pi }{\int }_{|\omega |>\pi /\delta }|\widehat{\theta }(\omega )|d\omega ,$$

where $S_{j\mathrm{,}\delta }\mathit{}\mathrm{(}x\mathrm{)}\mathrm{:=}\mathrm{sinc}\mathit{}\mathrm{(}x\mathrm{/}\delta \mathrm{-}j\mathrm{)}$ for $j\mathrm{\in }\mathrm{Z}$. In addition, if it holds for some $d\mathrm{>}\mathrm{0}$, $\mathrm{|}\widehat{\theta }\mathit{}\mathrm{(}\omega \mathrm{)}\mathrm{|}\mathrm{=}\mathrm{O}\mathit{}\mathrm{(}{\mathrm{e}}^{\mathrm{-}d\mathit{}\mathrm{|}\omega \mathrm{|}}\mathrm{)}$ for $\omega \mathrm{\to }\mathrm{\pm }\mathrm{\infty }$, and then, as $\delta \mathrm{\to }\mathrm{0}$,

$$\frac{1}{\delta }{\int }_{ℝ}\theta (x)S_{j,\delta }(x)dx-\theta (j\delta )=?({\mathrm{e}}^{-\pi d/\delta }).$$

It follows that we can approximate ${\mathrm{\Phi }}_{j,k}(\xi )$ by

	${\mathrm{\Phi }}_{j,k}(\xi )$	$\approx {\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi h(x)){\displaystyle \sum _{l\in \mathbb{Z}}}{c}_{m,l}^{j,k}{\phi }_{m,l}(x)\mathrm{d}x$
		$\approx {\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi h(x)){\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^{j,k}{\phi }_{m,l}(x)\mathrm{d}x$
		$\approx 2^{m/2}{\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^{j,k}{\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi h(x))\mathrm{sinc}(2^{m}x-l)dx$
		$\approx 2^{-m/2}{\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^{j,k}\exp \left(\mathrm{i}\xi h\left({\displaystyle \frac{l}{2^m}}\right)\right).$

This closed-form integral approximation simplifies the method’s implementation, especially compared with methods that require complicated quadrature algorithms. In particular, it avoids the truncation error induced by the numerical integration methods, since the integral is approximated in the whole domain, $ℝ$. Most importantly for the problem at hand, we can use the same generic pricing scheme to price a large variety of contracts without special handling.

Assumption 5.3.

Suppose that $h\mathit{}\mathrm{(}x\mathrm{)}\mathrm{:}\mathrm{R}\mathrm{\to }\mathrm{R}$ is such that $\theta \mathit{}\mathrm{(}x\mathrm{)}\mathrm{:=}\exp \mathit{}\mathrm{(}\mathrm{i}\mathit{}\xi \mathit{}h\mathit{}\mathrm{(}x\mathrm{)}\mathrm{)}$ satisfies, for some $d_h\mathrm{>}\mathrm{0}$, $\mathrm{|}\widehat{\theta }\mathit{}\mathrm{(}\omega \mathrm{)}\mathrm{|}\mathrm{=}\mathrm{O}\mathit{}\mathrm{(}{\mathrm{e}}^{\mathrm{-}{d}_h\mathit{}\mathrm{|}\omega \mathrm{|}}\mathrm{)}$ for $\omega \mathrm{\to }\mathrm{\pm }\mathrm{\infty }$, uniformly in $\xi$. We assume the same for $g\mathit{}\mathrm{(}x\mathrm{)}\mathrm{:}\mathrm{R}\mathrm{\to }\mathrm{R}$ for some $d_g\mathrm{>}\mathrm{0}$.

Lemma 5.4 (Quadrature error).

Let $\xi \mathrm{\in }\mathrm{R}$ and $j\mathrm{,}k\mathrm{\in }\mathrm{J}$, and define

$${?}_{j,k}(\xi ;h):={\int }_{ℝ}\exp (\mathrm{i}\xi h(x)){\tilde{f}}_R^{j,k}(x)dx-2^{-m/2}\sum _{l_1\le l\le l_2}{c}_{m,l}^{j,k}\exp \left(\mathrm{i}\xi h\left(\frac{l}{2^m}\right)\right).$$

If $h$ satisfies Assumption 5.3 for some $d_h\mathrm{>}\mathrm{0}$, then

$$\underset{\xi \in ℝ,j,k\in ?}{sup}|{?}_{j,k}(\xi ;h)|\le C_{?}^h(l_2-l_1+1){\mathrm{e}}^{-\pi d_h{2}^m}.$$

Proof.

See Appendix A (available online).∎

Similarly, the integral term defined by $Z_{n-1}$ in Proposition (5.1) is well approximated using a SWIFT representation of $f_{Y_{n-1}^{\alpha }}^j(y)$. In this case,

	${\varphi }_{Z_{n-1}}^j(\xi )$	$={\displaystyle {\int }_{ℝ}}\exp (\mathrm{i}\xi g(y))f_{Y_{n-1}^{\alpha }}^j(y)dy$
		$\approx 2^{-m/2}{\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^j\exp \left(\mathrm{i}\xi g\left({\displaystyle \frac{l}{2^m}}\right)\right),$

where

$$c_{m,l}^j:=⟨f_{Y_{n-1}^{\alpha }}^j,{\phi }_{m,l}⟩$$

are the projection coefficients for the density $f_{Y_{n-1}^{\alpha }}^j$ in terms of the ChF ${\varphi }_{Y_{n-1}^{\alpha }}^j$. The resulting error bounds are analogous to those in Lemma 5.4.

Remark 5.5.

The computational complexity to recover the terminal ChF, ${\varphi }_{Y_N^{\alpha }}^j(\xi )$, along a $\xi$ grid of $2^J$ points is determined as follows. First, from Remark 3.1, the full set of coefficients $c_{m,l}^{j,k}$ corresponding to the projection ${\tilde{f}}_R^{j,k}$ can be precomputed at a cost of $?(m_0^{2}J{2}^{J+1})$, while the matrix exponential $ℳ(\xi ;\mathrm{\Delta }t)$ in (2.9) costs $?(m_0^3{2}^J)$ using standard techniques. Precomputing ${\mathrm{\Phi }}_{j,k}(\xi )$ along the same grid requires $?(m_0^2{2}^J)$ operations. Thus, initialization requires $?(m_0^2(J+m_0)2^J)=?(m_0^3{2}^J)$ operations, where the number of states $m_0$ generally dominates $J$. The recursive algorithm in Proposition 5.1 then requires $?(N(m_0^2{2}^J+m_0{2}^{2J}))$ operations to calculate ${\varphi }_{Y_N^{\alpha }}^j(\xi )$, for a total computational complexity of

$$?(m_0^3{2}^J+N(m_0^2{2}^J+m_0{2}^{2J})).$$

This cost dominates the final valuation step described in Section 5.2.2 and is thus the computational complexity of pricing those exotic contracts falling within the proposed framework.

5.2.2 Option valuation

In the final stage of the SWIFT recursion, the ChF of $Y_N^{\alpha }$, ${\varphi }_{Y_N^{\alpha }}^j(\xi )$, is obtained for each of the starting states $j\in ?$. The goal is to price contracts of the form $G(A_N^{\alpha })$ (recall the beginning of Section 5). In particular, the payoffs can be expressed as $\tilde{G}(Y_N^{\alpha })=G(A_N^{\alpha })$. For example, the Asian option payoff with uniform weights can be written as

$$G(A_N^{\alpha })=G\left(\frac{1}{N+1}\sum _{n=0}^N{S}_n^{\alpha }\right)=G\left(\frac{S_0(1+e^{Y_N^{\alpha }})}{N+1}\right):=\tilde{G}(Y_N^{\alpha }),$$

(5.6)

and the contract value at $t_0=0$ can be obtained from the discounted expectation over $Y_N^{\alpha }$:

	$V(x,T)$	$={\mathrm{e}}^{-rT}?[h(Y_N^{\alpha },T)\mid x,\alpha (0)=j]$
		$={\mathrm{e}}^{-rT}{\displaystyle {\int }_{ℝ}}h(y,T)f_{Y_N^{\alpha }}^j(y)dy={\mathrm{e}}^{-rT}{\displaystyle {\int }_{ℝ}}\tilde{G}(y)f_{Y_N^{\alpha }}^j(y)dy,$

where $x:=\log (S_0)$. From the recovered ChF ${\varphi }_{Y_N^{\alpha }}^j(\xi )$, we again have the representation

$$f_{Y_N^{\alpha }}^j(y)\approx \sum _{l_1\le l\le l_2}{c}_{m,l}^j{\phi }_{m,l}(y):={\tilde{f}}_{Y_N^{\alpha }}^j(y).$$

It follows that

	$V(x,T)$	$\approx {\mathrm{e}}^{-rT}{\displaystyle {\int }_{ℝ}}\tilde{G}(y){\tilde{f}}_{Y_N^{\alpha }}^j(y)dy$
		$={\mathrm{e}}^{-rT}{\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^j{\displaystyle {\int }_a^b}\tilde{G}(y){\phi }_{m,l}(y)dy:={\mathrm{e}}^{-rT}{\displaystyle \sum _{l_1\le l\le l_2}}{c}_{m,l}^j{B}_{m,l},$

where

	$B_{m,l}$	$:={\displaystyle {\int }_a^b}\tilde{G}(y){\phi }_{m,l}(y)dy$
		$\mathrm{\hspace{0.17em}}\approx {\displaystyle \frac{2^{m/2}}{2^{J-1}}}{\displaystyle \sum _{\lambda =1}^{2^{J-1}}}{\displaystyle {\int }_a^b}\tilde{G}(y)\cos \left({\displaystyle \frac{2\lambda -1}{2^J}}\pi (2^{m}y-l)\right)dy,$		(5.7)

which follows from applying Vieta’s formula again (see Ortiz-Gracia and Oosterlee (2016) for details).

While the coefficients $B_{m,l}$ can be computed numerically for generic payoffs $\tilde{G}$, closed-form expressions exist for some of the more common contracts, such as Asian options, given in Appendix B.2.

5.3 Numerical experiments

In this section, we will present some experiments that aim to numerically validate the introduced CTMC–Heston model in the context of exotic option valuation. We will consider several exotic contracts: realized variance swaps, realized variance options and Asian options. The realized variance swaps are chosen for comparison purposes, since an exact solution of the Heston model is available (Bernard et al 2014). That is not the case for the other two products, which often require the use of MC methods.

The numerical experiments were carried out on a computer system with a CPU Intel Core i7-4720HQ 2.6GHz processor and memory of 16GB RAM, under the software package MATLAB R2017b.

Based on the numerical tests conducted in Section 4.2, the CTMC volatility grid is selected by the Tavella–Randall scheme, which has been proven to be the most reliable and consistent approach among those analyzed. The two volatility scenarios in Table 1 are again considered.

The accuracy is measured by the relative error (RE), defined in terms of the provided reference value. When the reference solutions are computed by MC, we employ the quadratic exponential discretization scheme (Andersen 2008) with ${10}^6$ paths and $360$ discretization time steps. Note that the SWIFT method does not require any time discretization; it simply iterates over the number of monitoring dates specified by the financial contract, $N$, as shown in (5.1).

5.3.1 Realized variance swaps

We first analyze the performance of our methodology by pricing realized variance swaps (see the payoff definition in (5.2)). As mentioned, there is an analytic solution for variance swaps under Heston dynamics (Bernard et al 2014), which is employed here as a reference value. As a first experiment, we study the convergence in terms of the number of CTMC states, $m_0$, for several realized variance swaps definitions, namely, a varying number of monitoring dates, $N=\{5,12,50180,360\}$. Further, for each volatility scenario in Table 1, two values for the correlation, $\rho =-0.1$ and $\rho =-0.7$, are used since, as observed in Figures 3 and 4, its impact in the approximation can be relatively important, determining as it does the influence of the stochastic volatility on the underlying dynamics. In Figures 7 and 8, the obtained convergence in $m_0$ is depicted. We can observe relatively stable patterns of convergence, reaching significantly acceptable precision, ${10}^{-5}-{10}^{-7}$, with $m_0=30\text{–}40$.

Figure 7: Convergence in $m_{\text{0}}$. Variance swaps: $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set I (regular market). (a) $\rho =-\text{0.1}$. (b) $\rho =-\text{0.7}$.

Figure 8: Convergence in $m_{\text{0}}$. Variance swaps: $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set II (stressed market). (a) $\rho =-\text{0.1}$. (b) $\rho =-\text{0.7}$.

Next, in Tables 2 and 3, the values for $m_0=40$ are presented. We include the MC prices and standard errors (s.e.) to establish a baseline for the cases without reference values. As expected, SWIFT overcomes MC in all of the cases, with the latter achieving an accuracy of ${10}^{-4}$ at most.

Table 2: Pricing realized variance swaps with SWIFT and MC: estimations and errors. [Variance swaps: $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set I.]

	(a) $\rho =-\text{0.1}$
$?$	Ref. (Bernard and Cui 2014)	SWIFT	MC (s.e.)	${\text{??}}_{\text{?????}}$	${\text{??}}_{\text{??}}$
5	0.03712054	0.03712058	0.03712422 $(\text{2.79}\times {\text{10}}^{-\text{5}})$	$\text{9.32}\times {\text{10}}^{-\text{7}}$	$\text{9.91}\times {\text{10}}^{-\text{5}}$
12	0.03695709	0.03695710	0.03696272 $(\text{2.01}\times {\text{10}}^{-\text{5}})$	$\text{4.06}\times {\text{10}}^{-\text{7}}$	$\text{1.52}\times {\text{10}}^{-\text{4}}$
50	0.03686316	0.03686300	0.03687360 $(\text{1.40}\times {\text{10}}^{-\text{5}})$	$\text{4.48}\times {\text{10}}^{-\text{6}}$	$\text{2.83}\times {\text{10}}^{-\text{4}}$
180	0.03684115	0.03684144	0.03683574 $(\text{1.22}\times {\text{10}}^{-\text{5}})$	$\text{8.00}\times {\text{10}}^{-\text{6}}$	$\text{1.46}\times {\text{10}}^{-\text{4}}$
360	0.03683689	0.03683422	0.03684662 $(\text{1.18}\times {\text{10}}^{-\text{5}})$	$\text{7.23}\times {\text{10}}^{-\text{5}}$	$\text{2.64}\times {\text{10}}^{-\text{4}}$
	(b) $\rho =-\text{0.7}$
$?$	Ref. (Bernard and Cui 2014)	SWIFT	MC (s.e.)	${\text{??}}_{\text{?????}}$	${\text{??}}_{\text{??}}$
5	0.03757379	0.03757400	0.03755392 $(\text{2.95}\times {\text{10}}^{-\text{5}})$	$\text{5.56}\times {\text{10}}^{-\text{6}}$	$\text{5.28}\times {\text{10}}^{-\text{4}}$
12	0.03716852	0.03716873	0.03715321 $(\text{2.10}\times {\text{10}}^{-\text{5}})$	$\text{5.55}\times {\text{10}}^{-\text{6}}$	$\text{4.11}\times {\text{10}}^{-\text{4}}$
50	0.03691728	0.03691720	0.03691997 $(\text{1.43}\times {\text{10}}^{-\text{5}})$	$\text{2.18}\times {\text{10}}^{-\text{6}}$	$\text{7.30}\times {\text{10}}^{-\text{5}}$
180	0.03685641	0.03685499	0.03685140 $(\text{1.23}\times {\text{10}}^{-\text{5}})$	$\text{3.83}\times {\text{10}}^{-\text{5}}$	$\text{1.35}\times {\text{10}}^{-\text{4}}$
360	0.03684454	0.03684890	0.03685224 $(\text{1.19}\times {\text{10}}^{-\text{5}})$	$\text{1.18}\times {\text{10}}^{-\text{4}}$	$\text{2.09}\times {\text{10}}^{-\text{4}}$

Table 3: Pricing realized variance swaps with SWIFT and MC: estimations and errors. [Variance swaps: $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set II.]

	(a) $\rho =-\text{0.1}$
$?$	Ref. (Bernard and Cui 2014)	SWIFT	MC (s.e.)	${\text{??}}_{\text{?????}}$	${\text{??}}_{\text{??}}$
5	0.40670787	0.40670865	0.40654230 $(\text{2.79}\times {\text{10}}^{-\text{4}})$	$\text{1.91}\times {\text{10}}^{-\text{6}}$	$\text{4.07}\times {\text{10}}^{-\text{4}}$
12	0.40290560	0.40290600	0.40279574 $(\text{1.88}\times {\text{10}}^{-\text{4}})$	$\text{9.99}\times {\text{10}}^{-\text{7}}$	$\text{2.72}\times {\text{10}}^{-\text{4}}$
50	0.40071390	0.40071394	0.40084167 $(\text{1.14}\times {\text{10}}^{-\text{4}})$	$\text{1.00}\times {\text{10}}^{-\text{7}}$	$\text{3.18}\times {\text{10}}^{-\text{4}}$
180	0.40019941	0.40019937	0.40022385 $(\text{8.87}\times {\text{10}}^{-\text{5}})$	$\text{1.06}\times {\text{10}}^{-\text{7}}$	$\text{6.10}\times {\text{10}}^{-\text{5}}$
360	0.40009981	0.40009976	0.40001819 $(\text{8.31}\times {\text{10}}^{-\text{5}})$	$\text{1.46}\times {\text{10}}^{-\text{7}}$	$\text{2.04}\times {\text{10}}^{-\text{4}}$
	(b) $\rho =-\text{0.7}$
$?$	Ref. (Bernard and Cui 2014)	SWIFT	MC (s.e.)	${\text{??}}_{\text{?????}}$	${\text{??}}_{\text{??}}$
5	0.41662864	0.41663104	0.41672516 $(\text{3.07}\times {\text{10}}^{-\text{4}})$	$\text{5.74}\times {\text{10}}^{-\text{6}}$	$\text{2.31}\times {\text{10}}^{-\text{4}}$
12	0.40751372	0.40751047	0.40730029 $(\text{2.01}\times {\text{10}}^{-\text{4}})$	$\text{7.98}\times {\text{10}}^{-\text{6}}$	$\text{5.23}\times {\text{10}}^{-\text{4}}$
50	0.40189025	0.40188670	0.40197021 $(\text{1.18}\times {\text{10}}^{-\text{4}})$	$\text{8.84}\times {\text{10}}^{-\text{6}}$	$\text{1.98}\times {\text{10}}^{-\text{4}}$
180	0.40053090	0.40052728	0.40066117 $(\text{9.00}\times {\text{10}}^{-\text{5}})$	$\text{9.03}\times {\text{10}}^{-\text{6}}$	$\text{3.25}\times {\text{10}}^{-\text{4}}$
360	0.40026602	0.40026239	0.40014305 $(\text{8.38}\times {\text{10}}^{-\text{5}})$	$\text{9.07}\times {\text{10}}^{-\text{6}}$	$\text{3.07}\times {\text{10}}^{-\text{4}}$

5.3.2 Realized variance options

As a second numerical test, realized variance options are priced, the payoff of which has been defined above. Contrary to its swaps counterpart, there is no closed-form solution for the valuation of this product. Hence, MC is commonly employed (included here as a reference) at a rather high computational cost. However, our unified framework allows us to efficiently address the realized variance option pricing using the SWIFT method and the introduced CTMC-based methodology ($m_0=40$). In Tables 4 and 5, the obtained option prices are presented. In this case, we have fixed the number of monitoring dates ($N=12$) and used several strike values, $K$. We can observe that our technique is very precise (with errors measured against MC), taking into account that, as seen in the previous experiment, MC only ensures $3\text{–}4$ significant digits. We therefore expect an even better performance from our method when it is compared with a more reliable reference.

Table 4: Pricing realized variance options with SWIFT: estimations and errors. [Variance call options: $r=\text{0.05}$, $q=\text{0}$, $T=\text{1}$ and $N=\text{12}$. Heston parameters: set I.]

(a) $\rho =-\text{0.1}$
$?$	Ref. MC (s.e.)	SWIFT	RE
0.01	0.02567765 $(\text{1.90}\times {\text{10}}^{-\text{5}})$	0.02568006	$\text{9.40}\times {\text{10}}^{-\text{5}}$
0.02	0.01699106 $(\text{1.81}\times {\text{10}}^{-\text{5}})$	0.01701443	$\text{1.37}\times {\text{10}}^{-\text{3}}$
0.03	0.01045427 $(\text{1.59}\times {\text{10}}^{-\text{5}})$	0.01044283	$\text{1.09}\times {\text{10}}^{-\text{3}}$
0.04	0.00613621 $(\text{1.31}\times {\text{10}}^{-\text{5}})$	0.00613145	$\text{7.75}\times {\text{10}}^{-\text{4}}$
0.05	0.00351388 $(\text{1.03}\times {\text{10}}^{-\text{5}})$	0.00352261	$\text{2.48}\times {\text{10}}^{-\text{3}}$
(b) $\rho =-\text{0.7}$
$?$	Ref. MC (s.e.)	SWIFT	RE
0.01	0.02587552 $(\text{1.99}\times {\text{10}}^{-\text{5}})$	0.02586686	$\text{3.34}\times {\text{10}}^{-\text{4}}$
0.02	0.01712666 $(\text{1.91}\times {\text{10}}^{-\text{5}})$	0.01710650	$\text{1.17}\times {\text{10}}^{-\text{3}}$
0.03	0.01053463 $(\text{1.70}\times {\text{10}}^{-\text{5}})$	0.01054260	$\text{7.57}\times {\text{10}}^{-\text{4}}$
0.04	0.00631007 $(\text{1.43}\times {\text{10}}^{-\text{5}})$	0.00633681	$\text{4.23}\times {\text{10}}^{-\text{3}}$
0.05	0.00380057 $(\text{1.16}\times {\text{10}}^{-\text{5}})$	0.00380673	$\text{1.62}\times {\text{10}}^{-\text{3}}$

Table 5: Pricing realized variance options with SWIFT: estimations and errors. [Variance call options: $r=\text{0.05}$, $q=\text{0}$, $T=\text{1}$ and $N=\text{12}$. Heston parameters: set II.]

(a) $\rho =-\text{0.1}$
$?$	Ref. MC (s.e.)	SWIFT	RE
0.1	0.28810430 $(\text{1.79}\times {\text{10}}^{-\text{4}})$	0.28826035	$\text{5.41}\times {\text{10}}^{-\text{4}}$
0.2	0.19753250 $(\text{1.73}\times {\text{10}}^{-\text{4}})$	0.19753794	$\text{2.75}\times {\text{10}}^{-\text{5}}$
0.3	0.12269943 $(\text{1.55}\times {\text{10}}^{-\text{4}})$	0.12276166	$\text{5.07}\times {\text{10}}^{-\text{4}}$
0.4	0.07050097 $(\text{1.27}\times {\text{10}}^{-\text{4}})$	0.07054838	$\text{6.72}\times {\text{10}}^{-\text{4}}$
0.5	0.03826162 $(\text{9.77}\times {\text{10}}^{-\text{5}})$	0.03836334	$\text{2.65}\times {\text{10}}^{-\text{3}}$
(b) $\rho =-\text{0.7}$
$?$	Ref. MC (s.e.)	SWIFT	RE
0.1	0.29269761 $(\text{1.91}\times {\text{10}}^{-\text{4}})$	0.29262732	$\text{2.40}\times {\text{10}}^{-\text{4}}$
0.2	0.20203744 $(\text{1.86}\times {\text{10}}^{-\text{4}})$	0.20184267	$\text{9.64}\times {\text{10}}^{-\text{4}}$
0.3	0.12730330 $(\text{1.68}\times {\text{10}}^{-\text{4}})$	0.12737500	$\text{5.63}\times {\text{10}}^{-\text{4}}$
0.4	0.07568155 $(\text{1.41}\times {\text{10}}^{-\text{4}})$	0.07567259	$\text{1.18}\times {\text{10}}^{-\text{4}}$
0.5	0.04341057 $(\text{1.12}\times {\text{10}}^{-\text{4}})$	0.04352151	$\text{2.55}\times {\text{10}}^{-\text{3}}$

5.3.3 Arithmetic Asian options

Arithmetic Asian options are particularly interesting exotic derivatives that are highly appreciated and widely traded in the market. They are a cheaper alternative with a less volatile option price than their European counterparts. As for the realized variance-like derivatives, the arithmetic Asian payoff definition falls within the recursive form previously presented. For each of the number of monitoring dates $N\in 12,50,250$, we consider several strike values, $K$, centered around an initial spot value of $S_0=100$. Again, the prices provided by SWIFT together with the CTMC–Heston model ($m_0=40$) are highly satisfactory, as demonstrated in Tables 6 and 7.

Table 6: Pricing Asian options with SWIFT: estimations and errors. [Asian call options: $S_{\text{0}}=\text{100}$, $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set I.]

(a) $N=\text{12}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	21.52858352 $(\text{9.94}\times {\text{10}}^{-\text{3}})$	21.52703662	$\text{7.18}\times {\text{10}}^{-\text{5}}$
090%	12.58238080 $(\text{8.98}\times {\text{10}}^{-\text{3}})$	12.58965477	$\text{5.78}\times {\text{10}}^{-\text{4}}$
100%	05.40026210 $(\text{6.56}\times {\text{10}}^{-\text{3}})$	05.40005466	$\text{3.84}\times {\text{10}}^{-\text{5}}$
110%	01.38805277 $(\text{3.33}\times {\text{10}}^{-\text{3}})$	01.39065989	$\text{1.87}\times {\text{10}}^{-\text{3}}$
120%	00.17363304 $(\text{1.07}\times {\text{10}}^{-\text{3}})$	00.17310340	$\text{3.05}\times {\text{10}}^{-\text{3}}$
(b) $N=\text{50}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	21.53863923 $(\text{1.00}\times {\text{10}}^{-\text{2}})$	21.53392805	$\text{2.18}\times {\text{10}}^{-\text{4}}$
090%	12.62396585 $(\text{9.05}\times {\text{10}}^{-\text{3}})$	12.61821271	$\text{4.55}\times {\text{10}}^{-\text{4}}$
100%	05.45042203 $(\text{6.61}\times {\text{10}}^{-\text{3}})$	05.44996341	$\text{8.41}\times {\text{10}}^{-\text{5}}$
110%	01.42955791 $(\text{3.38}\times {\text{10}}^{-\text{3}})$	01.42754719	$\text{1.40}\times {\text{10}}^{-\text{3}}$
120%	00.18249250 $(\text{1.10}\times {\text{10}}^{-\text{3}})$	00.18314737	$\text{3.58}\times {\text{10}}^{-\text{3}}$
(c) $N=\text{250}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	21.52663462 $(\text{1.00}\times {\text{10}}^{-\text{2}})$	21.53593134	$\text{4.31}\times {\text{10}}^{-\text{4}}$
090%	12.62698599 $(\text{9.07}\times {\text{10}}^{-\text{3}})$	12.62613250	$\text{6.75}\times {\text{10}}^{-\text{5}}$
100%	05.45348823 $(\text{6.63}\times {\text{10}}^{-\text{3}})$	05.46369394	$\text{1.87}\times {\text{10}}^{-\text{3}}$
110%	01.44408194 $(\text{3.40}\times {\text{10}}^{-\text{3}})$	01.43781010	$\text{4.34}\times {\text{10}}^{-\text{3}}$
120%	00.18757760 $(\text{1.11}\times {\text{10}}^{-\text{3}})$	00.18602981	$\text{8.25}\times {\text{10}}^{-\text{3}}$

Table 7: Pricing Asian options with SWIFT: estimations and errors. [Asian call options: $S_{\text{0}}=\text{100}$, $r=\text{0.05}$, $q=\text{0}$ and $T=\text{1}$. Heston parameters: set II.]

(a) $N=\text{12}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	25.55856787 $(\text{3.28}\times {\text{10}}^{-\text{2}})$	25.59888606	$\text{1.57}\times {\text{10}}^{-\text{3}}$
090%	19.66707259 $(\text{3.03}\times {\text{10}}^{-\text{2}})$	19.62786895	$\text{1.99}\times {\text{10}}^{-\text{3}}$
100%	14.89623827 $(\text{2.75}\times {\text{10}}^{-\text{2}})$	14.85527167	$\text{2.75}\times {\text{10}}^{-\text{3}}$
110%	11.15178957 $(\text{2.47}\times {\text{10}}^{-\text{2}})$	11.14635032	$\text{4.87}\times {\text{10}}^{-\text{4}}$
120%	08.31652993 $(\text{2.19}\times {\text{10}}^{-\text{2}})$	08.32127121	$\text{5.70}\times {\text{10}}^{-\text{4}}$
(b) $N=\text{50}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	25.78240367 $(\text{3.30}\times {\text{10}}^{-\text{2}})$	25.77787944	$\text{1.75}\times {\text{10}}^{-\text{4}}$
090%	19.82638995 $(\text{3.05}\times {\text{10}}^{-\text{2}})$	19.82724668	$\text{4.32}\times {\text{10}}^{-\text{5}}$
100%	15.05301658 $(\text{2.77}\times {\text{10}}^{-\text{2}})$	15.05297239	$\text{2.93}\times {\text{10}}^{-\text{6}}$
110%	11.32914392 $(\text{2.49}\times {\text{10}}^{-\text{2}})$	11.32709690	$\text{1.80}\times {\text{10}}^{-\text{4}}$
120%	08.46141915 $(\text{2.21}\times {\text{10}}^{-\text{2}})$	08.47720283	$\text{1.86}\times {\text{10}}^{-\text{3}}$
(c) $N=\text{250}$
$?\mathbf{(}\mathbf{\%}\mathbf{\text{ of }}{?}_{\text{?}}\mathbf{)}$	Ref. MC (s.e.)	SWIFT	RE
080%	25.86414658 $(\text{3.32}\times {\text{10}}^{-\text{2}})$	25.82553402	$\text{1.49}\times {\text{10}}^{-\text{3}}$
090%	19.92192034 $(\text{3.07}\times {\text{10}}^{-\text{2}})$	19.88065606	$\text{2.07}\times {\text{10}}^{-\text{3}}$
100%	15.12457605 $(\text{2.79}\times {\text{10}}^{-\text{2}})$	15.10643333	$\text{1.19}\times {\text{10}}^{-\text{3}}$
110%	11.37933056 $(\text{2.50}\times {\text{10}}^{-\text{2}})$	11.37652663	$\text{2.46}\times {\text{10}}^{-\text{4}}$
120%	08.52543663 $(\text{2.22}\times {\text{10}}^{-\text{2}})$	08.52037902	$\text{5.93}\times {\text{10}}^{-\text{4}}$

6 Conclusions

This work provides a general, computationally efficient and robust framework for pricing a wide class of exotic contracts under the CTMC–Heston model. We demonstrate that this model approximation provides a parsimonious and faithful representation of the now standard Heston model, and it is able to reproduce the same volatility smile structure with a modest number of states. Based on this approximation, we can efficiently price a large variety of contracts that are exceptionally difficult to handle under the Heston model as originally formulated. We demonstrate this advantage with Asian options and derivatives based on the discretely sampled realized variance. The efficiency of the method is obtained by combining Markov chain approximation of the underlying variance process with the SWIFT method. The sources of numerical error are analyzed, along with an extensive set of numerical experiments. There are several promising future research directions, including extensions to the case of the Heston model with stochastic interest rates, as in Grzelak and Oosterlee (2011), or the well-established SABR model (Hagan et al 2002) to parameterize the CTMC transition matrix. Markov chain approximation techniques have only recently been introduced in the option pricing literature, and there are potentially many applications that could benefit from this approach.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Acknowledgements

Álvaro Leitao acknowledges financial support from the Spanish Ministry of Science and Innovation through the María de Maeztu Programme for Units of Excellence in R&D (MDM-2014-0445). Luis Ortiz-Gracia thanks the Spanish Ministry of Economy and Competitiveness for funding under grants ECO2016-76203-C2-2 and MTM2016-76420-P (MINECO/FEDER, UE).

References

• Andersen, L. B. G. (2008). Efficient simulation of the Heston stochastic volatility model. The Journal of Computational Finance 11(3), 1–22 (https://doi.org/10.21314/JCF.2008.189).
• Bates, D. S. (1996). Jumps and stochastic volatility: exchange rate processes implicit in Deutsche mark options. Review of Financial Studies 9(1), 69–107 (https://doi.org/10.1093/rfs/9.1.69).
• Bayraktar, E., Dolinsky, Y., and Guo, J. (2018). Recombining tree approximations for optimal stopping for diffusions. SIAM Journal on Financial Mathematics 9, 602–633 (https://doi.org/10.1137/17M1118865).
• Benhamou, E., Gobet, E., and Miri, M. (2010). Time dependent Heston model. SIAM Journal on Financial Mathematics 1(1), 289–325 (https://doi.org/10.1137/090753814).
• Bernard, C., and Cui, Z. (2014). Prices and asymptotics for discrete variance swaps. Applied Mathematical Finance 21(2), 140–173 (https://doi.org/10.1080/1350486X.2013.820524).
• Bernard, C., Cui, Z., and McLeish, D. (2014). Convergence of the discrete variance swap in time-homogeneous diffusion models. Quantitative Finance Letters 2(1), 1–6 (https://doi.org/10.1080/21649502.2014.920513).
• Boyarchenko, S., and Levendorskii, S. (2014). Efficient variations of the Fourier transform in applications to option pricing. The Journal of Computational Finance 18(2), 57–90 (https://doi.org/10.21314/JCF.2014.277).
• Broadie, M., and Kaya, O. (2006). Exact simulation of option Greeks under stochastic volatility and jump diffusion models. Operations Research 54(2), 217–231 (https://doi.org/10.1287/opre.1050.0247).
• Cai, N., Song, Y., and Kou, S. (2015). A general framework for pricing Asian options under Markov processes. Operations Research 63(3), 540–554 (https://doi.org/10.1287/opre.2015.1385).
• Carr, P., and Madan, D. B. (1999). Option valuation using the fast Fourier transform. The Journal of Computational Finance 2, 61–73 (https://doi.org/10.21314/JCF.1999.043).
• Carverhill, A. P., and Clewlow, L. J. (1990). Flexible convolution: valuing average rate (Asian) options. Risk 3(4), 25–29.
• Cattani, C. (2008). Shannon wavelets theory. Mathematical Problems in Engineering 2008, 164808 (https://doi.org/10.1155/2008/164808).
• Chourdakis, K. (2004). Non-affine option pricing. Journal of Derivatives 11(3), 10–25 (https://doi.org/10.3905/jod.2004.391032).
• Colldeforns-Papiol, G., and Ortiz-Gracia, L. (2018). Computation of market risk measures with stochastic liquidity horizon. Journal of Computational and Applied Mathematics 342, 431–450 (https://doi.org/10.1016/j.cam.2018.03.038).
• Colldeforns-Papiol, G., Ortiz-Gracia, L., and Oosterlee, C. W. (2017). Two-dimensional Shannon wavelet inverse Fourier technique for pricing European options. Applied Numerical Mathematics 117(Supplement C), 115–138 (https://doi.org/10.1016/j.apnum.2017.03.002).
• Costabile, M., Leccadito, A., Massabó, I., and Russo, E. (2014). Option pricing under regime-switching jump-diffusion models. Journal of Computational and Applied Mathematics 256, 152–167 (https://doi.org/10.1016/j.cam.2013.07.046).
• Cui, Z., and Nguyen, D. (2017). First hitting time of integral diffusions and applications. Stochastic Models 33, pp. 376–391 (https://doi.org/10.1080/15326349.2017.1300920).
• Cui, Z., and Taylor, S. (2019). Pricing discretely monitored barrier options under Markov processes using a Markov chain approximation. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.3382236).
• Cui, Z., Kirkby, J. L., and Nguyen, D. (2017a). A general framework for discretely sampled realized variance derivatives in stochastic volatility models with jumps. European Journal on Operational Research 262, 381–400 (https://doi.org/10.1016/j.ejor.2017.04.007).
• Cui, Z., Kirkby, J. L., and Nguyen, D. (2017b). Equity-linked annuity pricing with cliquet-style guarantees in regime-switching and stochastic volatility models with jumps. Insurance: Mathematics and Economics 74, 46–62 (https://doi.org/10.1016/j.insmatheco.2017.02.010).
• Cui, Z., Kirkby, J. L., and Nguyen, D. (2018a). A general valuation framework for SABR and stochastic local volatility models. SIAM Journal on Financial Mathematics 9, 520–563 (https://doi.org/10.1137/16M1106572).
• Cui, Z., Lee, C., and Liu, Y. (2018b). Single-transform formulas for pricing Asian options in a general approximation framework under Markov processes. European Journal of Operational Research 266(3), 1134–1139 (https://doi.org/10.1016/j.ejor.2017.10.049).
• Dang, D.-M., and Ortiz-Gracia, L. (2018). A dimension reduction Shannon-wavelet based method for option pricing. Journal Scientific Computing 75(2), 733–761 (https://doi.org/10.1007/s10915-017-0556-y).
• Dang, D.-M., Nguyen, D., and Sewell, G. (2016). Numerical schemes for pricing Asian options under state-dependent regime switching jump diffusion models. Computers and Mathematics with Applications 71, 443–458 (https://doi.org/10.1016/j.camwa.2015.12.017).
• Daubechies, I. (1992). Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics, Philadelphia, PA (https://doi.org/10.1137/1.9781611970104).
• Elliott, R. J., Chan, L., and Siu, T. K. (2005). Option pricing and Esscher transform under regime switching. Annals of Finance 1(4), 423–432 (https://doi.org/10.1007/s10436-005-0013-z).
• Elliott, R. J., Siu, T. K., Chan, L., and Lau, J. W. (2007). Pricing options under a generalized Markov modulated jump-diffusion model. Stochastic Analysis and Applications 25, 821–843 (https://doi.org/10.1080/07362990701420118).
• Euch, O. E., and Rosenbaum, M. (2018). The characteristic function of rough Heston models. Mathematical Finance 29(1), 3–38 (https://doi.org/10.1111/mafi.12173).
• Fang, F., and Oosterlee, C. W. (2008). A novel pricing method for European options based on Fourier-cosine series expansions. SIAM Journal on Scientific Computing 31, 826–848 (https://doi.org/10.1137/080718061).
• Fang, F., and Oosterlee, C. W. (2011). A Fourier-based valuation method for Bermudan and barrier options under Heston’s model. SIAM Journal on Financial Mathematics 2, 439–463 (https://doi.org/10.1137/100794158).
• Fusai, G., and Kyriakou, I. (2016). General optimized lower and upper bounds for discrete and continuous arithmetic Asian options. Mathematics of Operations Research 41(2), 1–29 (https://doi.org/10.1287/moor.2015.0739).
• Geman, H., and Yor, M. (1993). Bessel processes, Asian options, and perpetuities. Mathematical Finance 3(4), 349–375 (https://doi.org/10.1111/j.1467-9965.1993.tb00092.x).
• Goutte, S., Ismail, A., and Pham, H. (2017). Regime-switching stochastic volatility model: estimation and calibration to VIX options. Applied Mathematical Finance 24(1), 38–75 (https://doi.org/10.1080/1350486X.2017.1333015).
• Grzelak, L. A., and Oosterlee, C. W. (2011). On the Heston model with stochastic interest rates. SIAM Journal on Financial Mathematics 2(1), 255–286 (https://doi.org/10.1137/090756119).
• Guennoun, H., Jacquier, A., Roome, P., and Shi, F. (2018). Asymptotic behavior of the fractional Heston model. SIAM Journal on Financial Mathematics 9(3), 1017–1045 (https://doi.org/10.1137/17M1142892).
• Hagan, P. S., Kumar, D., Lesniewski, A. S., and Woodward, D. E. (2002). Managing smile risk. Wilmott Magazine, pp. 84–108.
• He, X., and Zhu, S. (2017). How should a local regime-switching model be calibrated? Journal Economic Dynamics and Control 78, 149–163 (https://doi.org/10.1016/j.jedc.2017.03.005).
• He, X., and Zhu, S. (2018). On full calibration of hybrid local volatility and regime-switching models. Journal of Futures Markets 38, 586–606 (https://doi.org/10.1002/fut.21901).
• Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options. Review of Financial Studies 6(2), 327–343 (https://doi.org/10.1093/rfs/6.2.327).
• Jaber, E. A., and Euch, O. E. (2019). Multifactor approximation of rough volatility models. SIAM Journal on Financial Mathematics 10(2), 309–349 (https://doi.org/10.1137/18M1170236).
• Jacquier, A., and Shi, F. (2019). The randomized Heston model. SIAM Journal on Financial Mathematics 10(1), 89–129 (https://doi.org/10.1137/18M1166420).
• Janczura, J., and Weron, R. (2012). Efficient estimation of Markov regime-switching models: an application to electricity spot prices. Advances in Statistical Analysis 96(1), 385–407 (https://doi.org/10.1007/s10182-011-0181-2).
• Jiang, J., Liu, R. H., and Nguyen, D. (2016). A recombining tree method for option pricing with state-dependent switching rates. International Journal of Theoretical and Applied Finance 19, 1–26 (https://doi.org/10.1142/S0219024916500126).
• Jourdain, B., and Zhou, A. (2020). Existence of a calibrated regime switching local volatility model. Mathematical Finance 30, 501–546 (https://doi.org/10.1111/mafi.12231).
• Kirkby, J. L. (2015). Efficient option pricing by frame duality with the fast Fourier transform. SIAM Journal on Financial Mathematics 6(1), 713–747 (https://doi.org/10.1137/140989480).
• Kirkby, J. L. (2016). An efficient transform method for Asian option pricing. SIAM Journal on Financial Mathematics 7(1), 845–892 (https://doi.org/10.1137/16M1057127).
• Kirkby, J. L. (2017). Robust option pricing with characteristic functions and the B-spline order of density projection. The Journal of Computational Finance 21, 61–100 (https://doi.org/10.21314/JCF.2017.332).
• Kirkby, J. L., and Nguyen, D. (2020). Efficient Asian option pricing under regime-switching jump diffusions and stochastic volatility models. Annals of Finance 16, 307–351 (https://doi.org/10.1007/s10436-020-00366-0).
• Kirkby, J. L., Nguyen, D., and Cui, Z. (2017). A unified approach to Bermudan and barrier options under stochastic volatility models with jumps. Journal of Economic Dynamics and Control 80, 75–100 (https://doi.org/10.1016/j.jedc.2017.05.001).
• Leitao, A., Ortiz-Gracia, L., and Wagner, E. I. (2018). SWIFT valuation of discretely monitored arithmetic Asian options. Journal of Computational Science 28, 120–139 (https://doi.org/10.1016/j.jocs.2018.07.004).
• Levendorskii, S. (2018). Pricing arithmetic Asian options under Lévy models by backward induction in the dual space. SIAM Journal on Financial Mathematics 9(1), 1–27 (https://doi.org/10.1137/16M1108133).
• Li, L., and Zhang, G. (2017). Error analysis of finite difference and Markov chain approximations for option pricing with non-smooth payoffs. Mathematical Finance 28, 877–919 (https://doi.org/10.1111/mafi.12161).
• Li, L., and Zhang, G. (2019). Analysis of Markov chain approximation for option pricing and hedging: grid design and convergence behavior. Operations Research 67(2), 407–427.
• Linetsky, V. (2004). Spectral expansions for Asian (average price) options. Operations Research 52, 856–867 (https://doi.org/10.1287/opre.1040.0113).
• Lo, C. C., and Skindilias, K. (2014). An improved Markov chain approximation methodology: derivatives pricing and model calibration. International Journal of Theoretical and Applied Finance 17(07), 1450047 (https://doi.org/10.1142/S0219024914500472).
• Maree, S. C., Ortiz-Gracia, L., and Oosterlee, C. W. (2017). Pricing early-exercise and discrete barrier options by Shannon wavelet expansions. Numerische Mathematik 136(4), 1035–1070 (https://doi.org/10.1007/s00211-016-0858-2).
• Mijatović, A., and Pistorius, M. (2013). Continuously monitored barrier options under Markov processes. Mathematical Finance 23(1), 1–38 (https://doi.org/10.1111/j.1467-9965.2011.00486.x).
• Mitra, S., and Date, P. (2010). Regime switching volatility calibration by the Baum–Welch method. Journal of Computational and Applied Mathematics 234(12), 3243–3260 (https://doi.org/10.1016/j.cam.2010.04.022).
• Nguyen, D. (2018). A hybrid Markov chain-tree valuation framework for stochastic volatility jump diffusion models. International Journal of Financial Engineering 5(4), 1850039 (https://doi.org/10.1142/S2424786318500391).
• Ortiz-Gracia, L., and Oosterlee, C. W. (2013). Robust pricing of European options with wavelets and the characteristic function. SIAM Journal on Scientific Computing 35(5), B1055–1084 (https://doi.org/10.1137/130907288).
• Ortiz-Gracia, L., and Oosterlee, C. W. (2016). A highly efficient Shannon wavelet inverse Fourier technique for pricing European options. SIAM Journal on Scientific Computing 38(1), 118–143 (https://doi.org/10.1137/15M1014164).
• Pun, C. S., Chung, S. F., and Wong, H. Y. (2015). Variance swap with mean reversion, multifactor stochastic volatility and jumps. European Journal of Operational Research 245(2), 571–580 (https://doi.org/10.1016/j.ejor.2015.03.026).
• Ramponi, A. (2012). Fourier transform methods for regime-switching jump-diffusions and the pricing of forward starting options. International Journal of Theoretical and Applied Finance 15, 1–26 (https://doi.org/10.1142/S0219024912500379).
• Ruijter, M. J., and Oosterlee, C. W. (2012). Two-dimensional Fourier cosine series expansion method for pricing financial options. SIAM Journal on Scientific Computing 34, B642–B671 (https://doi.org/10.1137/120862053).
• Song, Y., Cai, N., and Kou, S. (2016). A unified framework for options pricing under regime-switching models. In Proceedings of the First PKU-NUS Annual International Conference on Quantitative Finance and Economics. HKUST SPD.
• Stenger, F. (2010). Handbook of Sinc Numerical Methods. CRC Press, Boca Raton, FL.
• Tavella, D., and Randall, C. (2000). Pricing Financial Instruments: The Finite Difference Method. Wiley.
• Zhang, B., and Oosterlee, C. W. (2013). Efficient pricing of European-style Asian options under exponential Lévy processes based on Fourier cosine expansions. SIAM Journal on Financial Mathematics 4(1), 399–426 (https://doi.org/10.1137/110853339).
• Zhang, G., and Li, L. (2019). Analysis of Markov chain approximation for diffusion models with non-smooth coefficients. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.3387751).
• Zheng, W., and Kwok, Y. K. (2014). Closed form pricing formulas for discretely sampled generalized variance swaps. Mathematical Finance 24(4), 855–881 (https://doi.org/10.1111/mafi.12016).
• Zhu, S.-P., and Lian, G.-H. (2011). A closed-form exact solution for pricing variance swaps with stochastic volatility. Mathematical Finance 21(2), 233–256 (https://doi.org/10.1111/j.1467-9965.2010.00436.x).

Download PDF Appendix 1

• Keywords
• Asian options
• Computational finance
• Heston process
• Model calibration
• Option pricing
• Stochastic volatility
• Variance swaps
• Original research
• Wavelets