Pricing Asian options in the Heston model

By Jacques Printems

Introduction

The following stochastic volatility model for the stock price dynamic in an incomplete market was introduced by Heston in 1993 [1]. Under a Risk-Neutral probability $ \mathbb{P} $, it writes:

$$
\left \{
\begin{array}{rclll}
dS_t & = & S_t \, ( r dt + \sqrt{v_t} \, dW_t^1), &  S_0 = s_0, & \\
&&& &\\
dv_t & = & k(a-v_t) + \vartheta \sqrt{v_t} \, dW_t^2, & v_0>0, & d\langle W^1, W^2 \rangle_t = \rho \, dt
\end{array}
\right.
$$

where $ \rho \in [-1,1] $ and where $ \theta, a, k >0 $ are such that $ \theta^2/(4ka) < 1 $. Here $ W^1 $ and $ W^2 $ are two standard Brownian motions under the probability measure $ \mathbb{P} $. Consider the Asian Call option of maturity $ T $ and strike $ K > 0 $ for which there is no explicit formula:

$$
e^{-rT} {\mathbb{E}} \, \left ( \frac 1T \int_0^T S_s ds - K \right )_+. 
$$

As a first step, we project $ W^1 $ onto $ W^2 $, so that $ W^1 = \rho W^2 + \sqrt{1-\rho^2} \widetilde{W^1} $ where $  \widetilde{W^1} $ is a standard Brownian motion independent of $ W^2 $ under the probability measure $ \mathbb{P} $. Then $ S_t $ writes

\begin{align*} 
S_t  & = s_0 \, \exp \left ( \Big ( r- \frac 12 \bar v_t \Big ) \, t + \rho \int_0^t \sqrt{v_s} \,dW^2_s \right ) \exp\left ( \sqrt{1 - \rho^2} \int_0^t \sqrt{v_s} d \widetilde W^1_s \right ) \\
	& = s_0 \exp \left ( \left (  r - \frac{\rho a k}{\vartheta}  + \bar v_t \Big ( \frac{\rho k}{\vartheta} - \frac 12 \Big )  \right ) \, t  + \frac{\rho}{\vartheta} \Big ( v_t - v_0\Big )  \right ) \, \exp \left ( \sqrt{1 - \rho^2} \int_0^t \sqrt{v_s} d \widetilde W^1_s \right ),
 \end{align*}

where $ \bar v_t = \frac 1t \int_0^t v_s ds $.

Consider now $ X = \{ X^1, \dots, X^{N_1} \} $ and $ Y = \{ Y^1, \dots, Y^{N_2} \} $ two functional quantizers of the Brownian motion.

\begin{align*} 
X^i(t) & = \sqrt{\frac{2}{T}} \sum_{n=1}^{d(N_1)} x_n^i \sin \left ( \frac{\pi t}{T} \left ( n - \frac 12 \right ) \right ), \quad 1\leq i \leq N_1\\
Y^j(t)  & =  \sqrt{\frac{2}{T}} \sum_{n=1}^{d(N_2)} y_n^j \sin \left ( \frac{\pi t}{T} \left ( n - \frac 12 \right ) \right ), \quad 1\leq j \leq N_2,
 \end{align*}

where tabs $ (x_n^i) $ andt $ (y_n^j) $ are available here.

For $ i \in \{ 1,\cdots, N_1 \} $ and $ j \in \{1,\cdots,N_2\} $, we numerically solve the following ordinary differential equations.

\begin{align*} 
\dot{y}_i(t) & = k \left ( a - y_i(t) - \frac{\theta^2}{4k} \right ) + \theta \sqrt{y_j(t)} \dot{X^i}(t), \quad y_i(0) = v_0, \\
\dot{s}_{i,j}(t)& = s_{i,j} (t)\left ( r - \frac 12 y_i(t) - \frac{\rho \vartheta}{4} + \sqrt{y_i(t)} \left ( \rho \dot{X^i}(t) + \sqrt{1-\rho^2} \, \dot{Y^j}(t) \right ) \right ), \quad s_{i,j}(0) = x.
 \end{align*}

(Here, we used a Runge Kutta IV method.) The option price is now approximated by

$$
e^{-rT} \sum_{i,j} \left ( \frac 1T \int_0^T s_{i,j}(t) dt-K \right)_+ \mathbb{P}\left(W^2 \in C_i(X) \right) \times {\mathbb{P}}(\widetilde{W}^1\in C_j(Y)).
$$

Numerical test

Here, we propose a method to compute the Asian call option price based on the functional quantization of the Brownian motion as described in Article [2] (Section 8).

Here $ (N_1,N_2) $ et $ (N_3,N_4) $ are the sizes of optimal functional quantizers of the standard Brownian motion ($ N_1 $ “points”, i.e. paths, for $ W^2 $ et $ N_2 $ for $ \widetilde{W}^1 $). Parameter $ n $ stands for the number of time-steps used to solve the ordinary differential equation written above. (We used a fourth order Runge Kutta scheme). Each grid couple yields a price approximation. The final result is a Romberg extrapolation between both prices, based on a $ O(1/\log(N_1 N_2)) $ rate of convergence for the quadrature error.

The following program was developed with the C programming language, and interfaced with Ruby. You can contact the authors for the source code.

Spot $ (x) $
Constant interest rate $ (r) $
Constant dividend rate
Initial variance $ (v_0) $
Mean reversion $ (k) $
Asymptotic variance $ (a) $
Volatility of volatility $ (\theta) $
Correlation $ (\rho) $
Maturity $ (T) $
Strike $ (K) $
$ N_1 $
$ N_2 $
$ N_3 $
$ N_4 $
$ n $
This can take up to one minute of computation before displaying the results.

The site team would like to thank David Delavennat (CNRS ingeneer at LAMA-UMR 8050 Univ. MLV - Paris 12 from 2003 to 2007) for his advice on the Ruby interface.


References