The following stochastic volatility model for the stock price dynamic in an incomplete market was introduced by Heston in 1993.
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, & \langle dW^1_t, dW^2_t \rangle = \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 ( \overline{S}_T - K \right )_+\]
where \(\overline{S}_T = \frac 1T \int_0^T S_t dt\).
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)\) and \((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)). \]