Welcome to THETAWIKI. If you like to create or edit a page please make sure to login or register an account. All registered users please make sure to provide a valid email address.

Schoebel Zhu Hull White

From ThetaWiki
Jump to: navigation, search

This article documents and implements the simulation discretization scheme for the Schöbel-Zhu Hull-White Model.

The Schöbel-Zhu Hull-White Model

The Schöbel-Zhu Hull-White (SZHW) Model as introduced in Haastrecht et al. (2009) [1] is set up under the risk-neutral measure $ \mathbb{Q} $ as follows.

The stock price dynamics $ S(t) $ is the following equation:

$ d S(t) = r(t) S(t) dt \ + \ v(t) S(t) \ dW_{S} (t), \ \ S(0) = S_0, \qquad (SZHW.1) $

The Schöbel-Zhu (1999) [2] volatility dynamics is given as:

$ d v(t) = \kappa \left [ \varphi - v(t) \right ] \ dt \ + \ \tau \ dW_{v} (t), \ \ v(0) = v_0, \qquad (SZHW.2) $

Hull-White(1990) [2] One Factor interest rate dynamics is the following equation:

$ d r(t) = \left [ \theta (t) - a \ r(t) \right ] \ dt \ + \ \sigma \ dW_{r} (t), \ \ r(0) = r_0, \qquad (SZHW.3) $

where in equation (SZHW.1), the interest rate $ r(t) $ is stochastic and follows the process defined in equation (SZHW.3). The stochastic volatility term $ v(t) $ in equation (SZHW.1) follows another diffusion process which is known as Schöbel-Zhu (1999) ) [2] stochastic volatility model as defined in equation (SZHW.2). For reference on the Schöbel-Zhu (1999) stochastic volatility model, see [1] and [2].

In equation (SZHW.2), the constant parameters $ \kappa $, $ \varphi $, $ \tau $ denote respectively the mean reversion speed, long-run volatility, and volatility of stochastic volatility. In equation (SZHW.3), the time-varying drift $ \theta (t) $ is a term used to fit exactly to the initial yield curve; the constant parameters $ a $ and $ \sigma $ are respectively the mean reversion speed and volatility for the short rate $ r(t) $. The function form for $ \theta (t) $ is derived in Appendix A of Hull-White model. Initial values for the stock price, stochastic volatility and interest rate are respectively $ S_0 $, $ v_0 $, and $ r_0 $.

The instantaneous correlations among the diffusion processes in equations (SZHW.1) - (SZHW.3) are specified as follows:

$ \begin{align} \langle dW_{S } (t), \ dW_{v }(t) \rangle = \rho_{SV} dt, \\ \langle dW_{S } (t), \ dW_{r }(t) \rangle = \rho_{Sr} dt, \\ \langle dW_{r } (t), \ dW_{v }(t) \rangle = \rho_{rv} dt, \\ \end{align} \qquad (SZHW.4) $

The Hull White One Factor zero coupon bond $ P(t, \ T) $ has the following process under the risk-neutral measure $ \mathbb{Q} $:

$ \frac{d P(t, \ T)}{ P(t, \ T)} = r(t) \ dt \ - \ \sigma B(t, \ T) \ dW_{r} (t), \qquad (SZHW.5) $

where $ B(t, \ T) = \frac{ 1 \ - \ e^{-a(T-t)} }{a} $.

To reduce the dimension of the above processes for Monte Carlo simulations, we choose to work under the forward measure $ \mathbb{Q}^T $ where the numeraire is the zero coupon bond price $ P(t, \ T) $. To do this, define the forward stock price as

$ F(t) = \frac{S(t)}{ P(t, \ T)}, \qquad (SZHW.6) $

Ito differentiation of equation (SZHW.6) yields

$ d F(t) = d \left (  \frac{ S(t) }{ P(t, \ T) }  \right ) $
$ = d \left (  \frac{1}{ P(t, \ T) } \right ) S(t) \ + \  \frac{1}{ P(t, \ T)} d S(t) \ + \   \left \langle  d \left (  \frac{1}{ P(t, \ T)}  \right ), \  d S(t)  \right \rangle $
$  = - \frac{ S(t) }{ P(t, \ T)^2} d P(t, \ T) \ + \  \frac{S(t)}{P(t, \ T)^3}  \left \langle  d P(t, \ T), \  d P(t, \ T)  \right \rangle $
$ \  +  \  \frac{1}{ P(t, \ T) } \left [ r(t) S(t) dt  \ + \ v(t) S(t) \ d W_S (t) \right ] \ + \ \frac{ S(t) }{ P(t, \ T) } \rho_{Sr} v(t) \sigma B(t, \ T) dt $
$ = \frac{S(t)}{P(t, \ T)} \left \{ -r(t) dt \ + \  \sigma  B(t, \ T)  \  d W_r (t)  \ + \   \sigma^2 B(t, \ T)^2 dt \ + \ \left [  r(t) dt \ + \ v(t) \  d W_S (t)  \right ]  \ + \  \rho_{Sr} v(t) \sigma B(t, \ T) dt \right \} $
$ = F(t)  \left \{ \left [  \sigma^2 B(t, \ T)^2  \ + \  \rho_{Sr} v(t) \sigma B(t, \ T)  \right ] \ dt  \ + \  \sigma  B(t, \ T)  \ d W_r (t)  \ + \ v(t) \ d W_S (t) \right \} $
$ = F(t)  \left \{   \sigma  B(t, \ T)  \ d W_r^T (t)  \ + \ v(t) \ d W_S^T (t) \right \} $
$ \Rightarrow d \ln F(t) = - \frac{1}{2} v_S^T (t)^2 dt \ + \  \sigma B(t, \ T) \ dW_r^T (t) \ + \ v(t) \ dW_S^T (t),  \qquad (SZHW.7) $

where

$ v_S^T (t)^2 = \sigma^2 B(t, \ T)^2 \ + \ v(t)^2 \ + \ 2 \ \rho_{Sr} v(t) \sigma B(t, \ T), \qquad (SZHW.7a) $
$ dW_r^T (t) = d W_r (t) \ + \ \sigma B(t, \ T) dt, \qquad (SZHW.7b) $
$ dW_S^T (t) = d W_S (t) \ + \ \rho_{Sr} \sigma B(t, \ T) dt, \qquad (SZHW.7c) $

and the terms $ dW_S^T (t) $, $ dW_r^T (t) $ and $ dW_v^T (t) $ are standard Brownian motions under the forward measure $ \mathbb{Q}^T $. For derivations of the representations in (SZHW.7a) - (SZHW.7c), please see Appendix A. Define the log forward stock price $ x(t) = \ln F(t) $, then from equation (SZHW.7) we have:

$ dx(t) = - \frac{1}{2} v_S^T (t)^2 dt \ + \ \sigma B(t, \ T) dW_r^T (t) \ + \ v(t) dW_S^T (t). \qquad (SZHW.8) $

The stochastic volatility process under the forward measure $ \mathbb{Q}^T $ is the following equation:

$ dv(t) = \kappa \left [ \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(t, \ T) \ - \ v(t) \right ] dt \ + \ \tau \ dW_v^T (t), \qquad (SZHW.9) $

where $ dW_v^T (t) = dW_v (t) \ + \ \rho_{rv} \sigma B(t, \ T) dt $.

Define the SZHW model stochastic variance $ V(t) = v(t)^2 $, Ito differentiation yields

$ \begin{align} & dV(t) = d [v(t)^2] \\ & = 2 v(t) dv(t) \ + \ \langle dv(t), \ dv(t) \rangle \\ & = 2 v(t) \left \{ \kappa \left [ \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(t, \ T) \ - \ v(t) \right ] dt \ + \ \tau \ dW_v^T (t) \right \} \ + \ \tau^2 dt \\ & = 2 \kappa \left \{ \frac{\tau^2}{2 \kappa} \ + \ \left ( \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(t, \ T ) \right ) v(t) \ - \ v(t)^2 \right \} dt \ + \ 2 \tau v(t) \ dW_v^T (t) \\ & = 2 \kappa \left \{ \frac{\tau^2}{2 \kappa} \ + \ \left ( \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(t, \ T ) \right ) \sqrt{ V(t) } \ - \ V(t) \right \} dt \ + \ 2 \tau \sqrt{ V(t) } \ dW_v^T (t). \qquad (SZHW.10) \end{align} $

where $ \sqrt{ V(t) } = v(t) $. Hench the stochastic variance process in the SZHW model is similar in form to the Heston (1993) stochastic variance process. When $ \varphi = \frac{\rho_{rv} \sigma \tau }{\kappa} B(t, \ T ) $ under the forward measure $ \mathbb{Q}^T $, the two process are virtually the same, with $ \kappa_H = 2 \kappa $, $ \theta = \frac{\tau^2}{2 \kappa} $, $ \epsilon = 2 \tau $, where $ \kappa_H $, $ \theta $, and $ \epsilon $ are respectively the mean reversion speed, the long-term variance and the volatility of variance for the Heston stochastic variance. For details of the Heston model, please see the reference Heston Volatility.


Discretization Scheme for the Schöbel-Zhu Stochastic Volatility Process under the Forward Measure $ \mathbb{Q}^T $

To simulate under the forward measure $ \mathbb{Q}^T $ the log forward stock price in equation (SZHW.8) and the stochastic volatility in (SZHW.9), we first discretize the Schöbel-Zhu (1999) volatility process as follows.

Integrate equation (SZHW.9) from $ t $ to $ t + \delta $, we have

$ v(t + \delta) = v(t) e^{- \kappa \delta t} \ + \ \kappa \int_t ^{ t + \delta } \left [ \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(u, \ T ) \right ] e^{- \kappa (t + \delta - u)} du $
$ \ + \ \int_t ^{ t + \delta } \tau e^{- \kappa (t + \delta - u)} dW_v^T(u), \qquad(SZHW.11) $

where

$ \kappa \int_t ^{ t + \delta } \left [ \varphi - \frac{\rho_{rv} \sigma \tau }{\kappa} B(u, \ T ) \right ] e^{- \kappa (t + \delta - u)} du $
$ = \left ( 1 \ - \ e^{- \kappa \delta} \right ) \left (\varphi - \frac{\rho_{rv} \sigma \tau }{a \kappa} \right ) \ + \ \frac{\rho_{rv} \sigma \tau }{a ( \kappa + a)} \left [ e^{-a (T - (t + \delta))} \ - \ e^{-a(T - t)} e^{- \kappa \delta} \right ] := K2 $,

and

$ \mathbb{E} \left [ \int_t ^{ t + \delta } \tau e^{- \kappa (t+\delta -u)} dW_v^T (u) \mid \mathcal{F}_t \right ]^2 $
$ = \int_t ^{ t + \delta } \tau^2 e^{-2 \kappa (t + \delta - u)} du = \frac{\tau^2}{2 \kappa} \left ( 1 \ - \ e^{- 2 \kappa \delta} \right ) := K3^2 $,

further define $ K1 = e^{- \kappa \delta t} $, we have the following simplified representation for equation (SZHW.11):

$ v(t + \delta) = K1 * v(t) \ + \ K2 \ + \ K3 * Z_v, \qquad (SZHW.12) $

where $ Z_v \sim \mathbb{N} (0, \ 1) $, i.e. $ Z_v $ is an i.i.d standard normal variable.

The integral representation for the SZHW model stochastic variance $ V(t) = v(t)^2 $ process in equation (SZHW.10) from time $ t $ to time $ t + \delta $ is the following:

$ \int _t^{t + \delta} d V(u) = \int _t^{t + \delta} d [ v(u)^2 ] $
$ = \int _t^{t + \delta} \left \{ 2 \kappa \left [ \frac{\tau^2}{2 \kappa} \ + \ \left ( \varphi \ - \ \frac{\rho_{rv} \sigma \tau }{\kappa} B(u, \ T) \right ) v(u) \ - \ v(u)^2 \right ] \right \} du \ + \ \int _t^{t+\delta} 2 \tau v(u) d W_v^T (u) $.

Rearranging and simplifying gives

$ \int _t^{t+\delta} v(u) d W_v^T (u) $
$ = \frac{1}{2 \tau} \left [ v(t + \delta)^2 \ - \ v(t)^2 \right ] \ - \ \frac{\tau}{2} \delta \ - \ \frac{\kappa \varphi }{\tau} \int _t^{t + \delta} v(u) \ du $
$ + \ \rho_{rv} \sigma \int _t^{t + \delta} B(u, \ T) v(u) \ du \ + \ \frac{\kappa}{\tau} \int _t^{t+\delta} v(u)^2 du. \qquad (SZHW.13) $


Discretization Scheme For the Log Forward Stock Price Process under the Forward Measure $ \mathbb{Q}^T $

Cholesky decomposition of the instantaneous correlation matrix (defined in equation (SZHW.4)) yields the following lower triangular matrix:

$ L = \begin{bmatrix} 1 & & \\ \rho_{Sv} & \sqrt{1 \ - \ \rho_{Sv}^2} & \\ \rho_{rv} & \omega_{Sr} & \sqrt{1 \ - \ \rho_{rv}^2 \ - \ \omega_{Sr}^2 } \\ \end{bmatrix}, \qquad (SZHW.14) $

with $ \omega_{Sr} = \frac{\rho_{Sr} - \rho_{Sv} \rho_{rv}}{\sqrt{1 \ - \ \rho_{Sv}^2}}. $

The log forward stock price dynamics in equation (SZHW.8) can then be re-expressed in terms of uncorrelated Brownian motions as

$ dx(t) = -\frac{1}{2} v_S^T (t)^2 dt \ + \ v(t) \left ( \rho_{Sv} d W_v^T (t) \ + \ \sqrt{1 \ - \rho_{Sv}^2 } \ d \hat{W}_S^T (t) \right ) $
$ + \ \sigma B(t, \ T) \left ( \rho_{rv} d \hat{W}_v^T (t) \ + \ \omega_{Sr} d \hat{W}_S^T (t) \ + \sqrt{1 \ - \ \rho_{rv}^2 \ - \ \omega_{Sr}^2 } \ d \hat{W}_r^T (t) \right ), \qquad (SZHW.15) $

where $ d \hat{W}_S^T $, $ d W_v^T $ and $ d \hat{W}_r^T $ are independent standard Brownian Motions under the forward measure $ \mathbb{Q}^T $.

The integral representation of equation (SZHW.15) over the time step $ \delta $ is the following:

$ x(t + \delta) - x(t) $
$ = - \frac{1}{2} \int _t^{t+\delta} v_S^T (u)^2 du \ + \rho_{Sv} \int _t^{t+\delta} v(u) d W_v^T (u) \ + \ \rho_{rv} \int _t^{t+\delta} \sigma B(u, \ T) d W_v^T (u) $
$ \ + \ \int _t^{t+\delta} \left [ \sqrt{1 \ - \ \rho_{Sv}^2 } \ v(u) \ + \ \omega_{Sr} \sigma B(u, \ T) \right ] d \hat{W}_S^T (u) $
$ \ + \ \sqrt{1 \ - \ \rho_{rv}^2 \ -\ \omega_{Sr}^2 } \ \sigma \int _t^{t+\delta} B(u, \ T) d \hat{W}_r^T (u). \qquad(SZHW.16) $

Substituting in the expression for $ v_S^T (t)^2 $ defined in equation (SZHW.7a) and the equation for $ \int _t^{t+\delta} v(u) d W_v^T (u) $ in (SZHW.13), we have

$ x(t + \delta) - x(t) $
$ = - \frac{1}{2} \int _t^{t+\delta} \left [ \sigma^2 B(u, \ T)^2 \ + \ v(u)^2 \ + \ 2 \rho_{Sr} v(u) \sigma B(u, \ T) \right ] du $
$ \ + \ \left \{ \frac{ \rho_{Sv}}{2 \tau} \left [ v(t + \delta)^2 \ - \ v(t)^2 \ - \ \tau^2 \delta \right ] \ - \ \frac{ \rho_{Sv} \kappa \varphi}{\tau} \int_t^{t+\delta} v(u) du \ + \ \rho_{Sv} \rho_{rv} \sigma \int_t^{t+\delta} B(u, \ T) v(u) du \ + \ \frac{ \rho_{Sv} \kappa }{\tau} \int_t^{t+\delta} v(u)^2 du \right \} $
$ \ + \ \rho_{rv} \int_t^{t+\delta} \sigma B(u, \ T) d W_v^T (u) \ + \ \int_t^{t+\delta} \left [ \sqrt{1 \ - \ \rho_{Sv}^2} v(u) \ + \ \omega_{Sr} \ \sigma B(u, \ T) \right ] d \hat{W}_S^T (u) $
$ \ + \ \sqrt{1 \ - \ \rho_{rv}^2 \ - \ \omega_{Sr}^2 } \ \sigma \int_t^{t+\delta} B(u, \ T) d \hat{W}_r^T (u). \qquad(SZHW.17) $

We use the following approximations for the integrals $ \int_t^{t+\delta} v(u)^p du $:

$ \left ( \gamma_1 \ v(t)^p \ + \ \gamma_2 \ v(t + \delta)^p \right ) \delta, \qquad(SZHW.18) $

where $ v(t)^p $ denotes the $ p $th power of $ v(t) $, and $ \gamma_1, \ \gamma_2 $ are discretization constants, when $ \gamma_1 = \gamma_2 = 0.5 $, the central discretization scheme is used.

Furthermore, define

$ \sigma \int_t^{t+\delta} B(u, \ T) du $
$ = \sigma \int_t^{t+\delta} \frac{ 1 - e^{-a(T-u)} }{a} du $
$ = \frac{\sigma}{a} \left ( \delta \ - \ \frac{ e^{ -a(T-(t+\delta)) } \ - \ e^{-a(T-t)} }{a} \right ) := H(t, \ t + \delta), \qquad(SZHW.19) $


$ \sigma^2 \int_t^{t+\delta} B(u, \ T)^2 du $
$ = \sigma^2 \int_t^{t+\delta} \left ( \frac{ 1 - e^{-a(T-u)} }{a} \right )^2 du $
$ = \frac{\sigma^2}{a^2} \left ( \delta \ + \ \frac{ e^{ -2a(T-(t+\delta)) } \ - \ e^{ -2a(T-t) } }{2a} \ - \ 2 \frac{ e^{ -a(T-(t+\delta)) } \ - \ e^{-a(T-t)} }{a} \right ) := G(t, \ t + \delta), \qquad(SZHW.20) $


$ \int_t^{t+\delta} \left [ \sqrt{1- \rho_{Sv}^2} \ v(u) \ + \ \omega_{Sr} \sigma B(u, \ T) \right ] d \hat{W}_S^T (u) \thicksim \mathbb{N} (0, \ \sigma_{\hat{W}_S^T}^2), \qquad(SZHW.21) $

where

$ \sigma_{\hat{W}_S^T}^2 $
$ = \int_t^{t+\delta} \left [ (1 - \rho_{Sv}^2) v(u)^2 \ + \ \omega_{Sr}^2 \sigma^2 B(u, \ T)^2 \ + \ 2 \omega_{Sr} \sqrt{1 - \rho_{Sv}^2} \ \sigma v(u) B(u, \ T) \right ] du $
$ \approxeq (1 - \rho_{Sv}^2) \left ( \gamma_1 \ v(t)^2 \ + \ \gamma_2 \ v(t + \delta)^2 \right ) \delta \ + \ \omega_{Sr}^2 G(t, \ t + \delta) $
$ \ + \ 2 \omega_{Sr} \sqrt{1 - \rho_{Sv}^2} \left ( \gamma_1 \ v(t) \ + \ \gamma_2 \ v(t + \delta) \right ) H(t, \ t + \delta). $

Conditional on $ d W_v^T (u) $ - the independent driving Brownian motion for Schöbel-Zhu stochastic volatility process under the forward measure $ \mathbb{Q}^T $, the covariance of the integral term $ \int_t^{t + \delta} \sigma B(u, \ T) d W_v^T (u) $ is given as:

$ \int_t^{t + \delta} \sigma B(u, \ T) d W_v^T (u) \mid \int_t^{t + \delta} d W_v^T (u) $
$ \thickapprox \sqrt{G(t, \ t + \delta)} \left ( \rho_{vv2}(t, \ t + \delta) \ Z_v \ + \ \sqrt{1 - \rho_{vv2}(t, \ t + \delta)^2} \ Z_{v2} \right ), \qquad(SZHW.22) $

where $ \rho_{vv2} = \frac{\int_t^{t + \delta} \sigma B(u, \ T) du}{\sqrt{\delta \ G(t, \ t+ \delta)} } $, and $ Z_v, \ Z_{v2} $ are independent Gaussian random variables. Plugging equations (SZHW.18) - (SZHW.22) into equation (SZHW.17), simplifying and collecting terms, we have the following discretized process for the log forward stock price:

$ x(t + \delta) $
$ = x(t) \ + \ \left ( -\frac{1}{2} G(t, \ t + \delta) \ - \ \frac{\rho_{Sv} \tau \delta}{2} \right ) $
$ \ - \ \gamma_1 \left \{ \rho_{Sr} H(t, \ t + \delta) \ + \ \rho_{Sv} \left [ \frac{\kappa \varphi}{\tau} \delta \ - \ \rho_{rv} H(t, \ t + \delta) \right ] \right \} v(t) $
$ \ - \ \gamma_2 \left \{ \rho_{Sr} H(t, \ t + \delta) \ + \ \rho_{Sv} \left [ \frac{\kappa \varphi}{\tau} \delta \ - \ \rho_{rv} H(t, \ t + \delta) \right ] \right \} v(t + \delta) $
$ \ + \ \left \{ \left ( \frac{\rho_{Sv} \kappa}{\tau} \ - \ \frac{1}{2} \right ) \gamma_1 \delta \ - \ \frac{\rho_{Sv} }{2 \tau} \right \} v(t)^2 $
$ \ + \ \left \{ \left ( \frac{\rho_{Sv} \kappa}{\tau} \ - \ \frac{1}{2} \right ) \gamma_2 \delta \ + \ \frac{\rho_{Sv} }{2 \tau} \right \} v(t + \delta)^2 $
$ \ + \ \sqrt{\sigma_{\hat{W}_S^T}^2} Z_S \ + \ \rho_{rv} \sqrt{G(t, \ t + \delta)} \left ( \rho_{vv2}(t, \ t + \delta) \ Z_v \ + \ \sqrt{1 - \rho_{vv2}(t, \ t + \delta)^2} \ Z_{v2} \right ) $
$ \ + \ \sqrt{G(t, \ t + \delta)} \sqrt{1 - \rho_{rv}^2 - \omega_{Sr}^2} \ Z_r, \qquad(SZHW.23) $

where $ Z_S, \ Z_{v2} $ and $ Z_r $ are independent standard normal variables, and $ Z_v = \frac{v(t + \delta) - K1 \ v(t) \ - \ K2}{K3} $, by rearranging equation (SZHW.12).

Define

$ C0 = \left ( - \frac{1}{2} G(t, \ t + \delta) \ - \ \frac{\rho_{Sv} \tau \delta}{2} \right ) $,
$ C1 = - \gamma_1 \left \{ \rho_{Sr} H(t, \ t + \delta) \ + \ \rho_{Sv} \left [ \frac{\kappa \varphi }{\tau} \delta \ - \ \rho_{rv} H(t, \ t + \delta) \right ] \right \} $,
$ C2 = - \gamma_2 \left \{ \rho_{Sr} H(t, \ t + \delta) \ + \ \rho_{Sv} \left [ \frac{\kappa \varphi }{\tau} \delta \ - \ \rho_{rv} H(t, \ t + \delta) \right ] \right \} $,
$ C3 = \left ( \frac{\rho_{Sv} \kappa}{\tau} \ - \ \frac{1}{2} \right ) \gamma_1 \delta \ - \ \frac{\rho_{Sv}}{2 \tau} $,
$ C4 = \left ( \frac{\rho_{Sv} \kappa}{\tau} \ - \ \frac{1}{2} \right ) \gamma_2 \delta \ + \ \frac{\rho_{Sv}}{2 \tau} $,
$ C5 = \sqrt{\sigma_{\hat{W}_S^T}^2} $,
$ C50 = \omega_{Sr}^2 G(t, \ t + \delta) $,
$ C51 = 2 \omega_{Sr} \sqrt{1 - \rho_{Sv}^2} H(t, \ t + \delta) \gamma_1 $,
$ C52 = 2 \omega_{Sr} \sqrt{1 - \rho_{Sv}^2} H(t, \ t + \delta) \gamma_2 $,
$ C53 = (1 \ - \ \rho_{Sv}^2) \gamma_1 \ \delta $,
$ C54 = (1 \ - \ \rho_{Sv}^2) \gamma_2 \ \delta $,
$ C6 = \rho_{rv} \rho_{vv2}(t, \ t + \delta) \sqrt{G(t, \ t + \delta)} $,
$ C7 = \rho_{rv} \sqrt{1 - \rho_{vv2}(t, \ t + \delta)^2} \sqrt{G(t, \ t + \delta)} $,
$ C8 = \sqrt{1 - \rho_{rv}^2 - \omega_{Sr}^2} \sqrt{G(t, \ t + \delta)} $.

The simplified representation for equation (SZHW.23) in terms of Schöbel-Zhu (1993) stochastic volatility $ v(t) $ is the following:

$ x(t + \delta) = x(t) \ + \ C0 \ + \ C1 * v(t) \ + \ C2 * v(t + \delta) \ + \ C3 * v(t)^2 \ + \ C4 * v(t + \delta)^2 $
$ \ + \ C5 * Z_S \ + \ C6 * Z_v \ + \ C7 * Z_{v2} \ + \ C8 * Z_r, \qquad(SZHW.24) $

where

$ C5 = \sqrt{C50 \ + \ C51 * v(t) \ + \ C52 * v(t + \delta) \ + \ C53 * v(t)^2 \ + \ C54 * v(t + \delta)^2} $.

Taking exponential of both sides of equation (SZHW.24), the discretized forward stock price $ F(t) $ under the forward measure $ \mathbb{Q}^T $ is

$ F(t + \delta) $
$ = F(t) e^{ C0 \ + \ C1 * v(t) \ + \ C2 * v(t + \delta) \ + \ C3 * v(t)^2 \ + \ C4 * v(t + \delta)^2 \ + \ C5 * Z_S \ + \ C6 * Z_v \ + \ C7 * Z_{v2} \ + \ C8 * Z_r }. \qquad(SZHW.25) $


Martingale Correction for the Discretized Stock Price Process under the Forward Measure $ \mathbb{Q}^T $

The discretized forward stock price $ F(t) $ is not a Martingale due to possible discretization errors. To make it a Martingale, we do the Martingale correction for the drift in the forward stock prices $ F(t) $ as in the following. Taking the $ \mathbb{Q}^T $ expectation of $ F(t + \delta) $ conditional on $ F(t) $, and using the iterated conditional expectations property, we have:

$ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t) \right ] = \mathbb{E}^{\mathbb{Q}^T} \left \{ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t), v(t + \delta) \right ] \mid F(t) \right \}. \qquad(SZHW.26) $

Substituting equation (SZHW.25) into equation (SZHW.26) gives

$ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t) \right ] $
$ = \mathbb{E}^{\mathbb{Q}^T} \left \{ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t), v(t + \delta) \right ] \mid F(t) \right \} $
$ = \mathbb{E}^{\mathbb{Q}^T} \left \{ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t) e^{ C0 \ + \ C1 * v(t) \ + \ C2 * v(t + \delta) \ + \ C3 * v(t)^2 \ + \ C4 * v(t + \delta)^2 \ + \ C5 * Z_S \ + \ C6 * Z_v \ + \ C7 * Z_{v2} \ + \ C8 * Z_r } \mid F(t), v(t + \delta) \right ] \mid F(t) \right \} $
$ = \mathbb{E}^{\mathbb{Q}^T} \left \{ \left [ F(t) e^{ C0 \ + \ C1 * v(t) \ + \ C2 * v(t + \delta) \ + \ C3 * v(t)^2 \ + \ C4 * v(t + \delta)^2 \ + \ \frac{1}{2} C5^2 \ + \ \frac{1}{2} C7^2 \ + \ \frac{1}{2} C8^2 } e^{C6 \frac{v(t + \delta) - K1 * v(t) - K2}{K3}} \right ] \mid F(t) \right \} $
$ = \mathbb{E}^{\mathbb{Q}^T} \left \{ \left [ F(t) e^{ C0 \ + \ C1 * v(t) \ + \ C2 * v(t + \delta) \ + \ C3 * v(t)^2 \ + \ C4 * v(t + \delta)^2 \ + \ \frac{1}{2} C7^2 \ + \ \frac{1}{2} C8^2 } e^{ \frac{1}{2} \left ( C50 \ + \ C51 * v(t) \ + \ C52 * v(t + \delta) \ + \ C53 * v(t)^2 \ + \ C54 * v(t + \delta)^2 \right ) } e^{C6 \frac{v(t + \delta) - K1 * v(t) - K2}{K3} } \right ] \mid F(t) \right \} $
$ = \mathbb{E}^{\mathbb{Q}^T} \left \{ \left [ F(t) e^{ C0 \ + \ \frac{1}{2} C50 \ + \ C6 \frac{-K2}{K3} \ + \ \frac{1}{2} C7^2 \ + \ \frac{1}{2} C8^2 } e^{ \left ( C1 \ + \ \frac{1}{2} C51 \ - \ \frac{C6*K1}{K3} \right ) v(t)} e^{ \left ( C2 \ + \ \frac{1}{2} C52 \ + \ \frac{C6}{K3} \right ) v(t + \delta)} e^{ \left ( C3 \ + \ \frac{1}{2} C53 \right ) v(t)^2} e^{ \left ( C4 \ + \ \frac{1}{2} C54 \right ) v(t + \delta)^2 } \right ] \mid F(t) \right \}, \qquad(SZHW.27) $

Denote

$ D0 = C0 \ + \ \frac{1}{2} C50 \ + \ C6 \frac{-K2}{K3} \ + \ \frac{1}{2} C7^2 \ + \ \frac{1}{2} C8^2 $,
$ D1 = C1 \ + \ \frac{1}{2} C51 \ - \ \frac{C6*K1}{K3} $,
$ D2 = C2 \ + \ \frac{1}{2} C52 \ + \ \frac{C6}{K3} $,
$ D3 = C3 \ + \ \frac{1}{2} C53 $,
$ D4 = C4 \ + \ \frac{1}{2} C54 $.

Equation (SZHW.27) can be re-expressed in a simplified manner as

$ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t) \right ] $
$ = F(t) e^{D0 \ + \ D1 * v(t) \ + \ D3 * v(t)^2} \mathbb{E}^{\mathbb{Q}^T} \left \{ \left [ e^{D2 * v(t + \delta) \ + \ D4 * v(t + \delta)^2} \right ] \mid F(t) \right \}. \qquad(SZHW.28) $

For $ F(t + \delta) $ to be a Martingale, the drift $ C0 $ in $ D0 $ needs to be corrected as $ C0^{*} $, with

$ D0^{*} = C0^{*} \ + \ \frac{1}{2} C50 \ + \ C6 \frac{-K2}{K3} \ + \ \frac{1}{2} C7^2 \ + \ \frac{1}{2} C8^2 $.

To have a representation for $ C0^{*} $, first denote

$ \boldsymbol{\Psi}_Y = \mathbb{E}^{\mathbb{Q}^T} \left \{ \left [ e^{D2 * v(t + \delta) \ + \ D4 * v(t + \delta)^2} \right ] \mid F(t) \right \}. \qquad(SZHW.29) $

Equation (SZHW.28) expressed as a Martingale expectation is then

$ \mathbb{E}^{\mathbb{Q}^T} \left [ F(t + \delta) \mid F(t) \right ] = F(t) $
$ = F(t) e^{ D0^{*} \ + \ D1 * v(t) \ + \ D3 * v(t)^2 } \boldsymbol{\Psi}_Y. \qquad(SZHW.30) $

Taking log of both sides of equation (SZHW.30) and rearranging gives

$ D0^{*} = - D1 * v(t) \ - \ D3 * v(t)^2 \ - \ \ln \boldsymbol{\Psi}_Y $
$ \Rightarrow C0^{*} = - D1 * v(t) \ - \ D3 * v(t)^2 \ - \ \ln \boldsymbol{\Psi}_Y \ - \ \frac{1}{2} C50 \ - \ C6 \frac{-K2}{K3} \ - \ \frac{1}{2} C7^2 \ - \ \frac{1}{2} C8^2. \qquad(SZHW.31) $

To compute $ \boldsymbol{\Psi}_Y $ in equation (SZHW.29), define

$ Y : = D2 * v(t + \delta) \ + \ D4 * v(t + \delta)^2 $,

and note that $ v(t + \delta) $ is normally distributed, the moment generating function of $ Y $ is

$ \boldsymbol{\Psi}_Y = \mathbb{E}^{\mathbb{Q}^T} \left [ e^{Y} \ \mid F(t) \right ] = e^{ - \frac{D2^2}{4D4}} \frac{ e^{ -\frac{\lambda D4 * \sigma_v^2}{1 - 2 D4 * \sigma_v^2}}}{\sqrt{1 - 2 D4 * \sigma_v^2}}, \qquad(SZHW.32) $

with the regularity condition $ 1 – 2 D4 * \sigma_v^2 > 0 $ satisfied, and

$ \lambda = \left ( \frac{ \mu_v \ + \ \frac{D2}{2D4} }{\sigma_v} \right )^2, \qquad(SZHW.33) $
$ \mu_v = \mathbb{E}^{\mathbb{Q}^T} \left [ v(t + \delta) \mid v(t) \right ] = K1 * v(t) \ + \ K2, \qquad(SZHW.34) $
$ \sigma_v^2 = Var^{\mathbb{Q}^T} \left [ v(t + \delta) \mid v(t) \right ] = K3^2. \qquad(SZHW.35) $

Substituting equations (SZHW.33), (SZHW.34) and (SZHW.35) into equation (SZHW.32), taking logs gives

$ \ln \boldsymbol{\Psi}_Y = - \frac{D2^2}{4 D4} \ - \ \frac{1}{2} \ln (1 - 2 D4 * K3^2) \ + \ \frac{ D4 \left ( K2 + \frac{D2}{2 D4} \right )^2 }{1 - 2 D4 * K3^2} $
$ \ \ + \frac{ 2 \left ( K2 + \frac{D2}{2 D4} \right ) D4 * K1 }{1 - 2 D4 * K3^2} v(t) \ + \ \frac{ D4 * K1^2}{1 - 2 D4 * K3^2} v(t)^2, \qquad(SZHW.36) $
$ \Rightarrow C0^{*} = - \frac{1}{2} C50 \ + \ C6 \frac{K2}{K3} \ - \ \frac{1}{2} C7^2 \ - \ \frac{1}{2} C8^2 \ + \ \frac{D2^2}{4 D4} \ + \ \frac{1}{2} \ln (1 - 2 D4 * K3^2) $
$ \ \ - \frac{ D4 \left ( K2 + \frac{D2}{2 D4} \right )^2 }{1 - 2 D4 * K3^2} \ + \ \left ( -D1 \ - \ \frac{2 \left ( K2 + \frac{D2}{2 D4} \right ) D4 * K1}{1 - 2 D4 * K3^2} \right ) * v(t) $
$ \ \ + \left ( -D3 \ - \ \frac{D4 * K1^2}{1 - 2 D4 * K3^2} \right ) * v(t)^2. \qquad(SZHW.37) $

Replacing the term $ C0 $ in equation (SZHW.25) with the term $ C0^{*} $ in equation (SZHW.37), the discretized forward stock price $ F(t + \delta) $ is drift Martingale corrected for the simulation scheme. 

The Schöbel-Zhu Hull-White model implemented in ThetaML code

This section gives the code implementation of the Schöbel-Zhu Hull-White model. The simulated Hull-White discount bond prices PtT are implemented in Hull-white model and is called in the model SchoebelZhuHullWhite.

model SchoebelZhuHullWhite
%this model simulates a stock price process with stochastic volatility and stochastic interest rate;
%the volatility dynamics is the Schöbel-Zhu model and the interest rate has the Hull-White
%one factor dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Test case parameters:
%Case I:   x0=100, v0=0.2, k=0.0063, tht=0.6446, eps=0.0124, 
%          rho_xr=-0.3864, rho_xv=-0.6714,rho_rv=0.177,T=10, a=0.01613, sigma=0.012812
%Case II:  x0=100, v0=0.2, k=0.4, tht=0.2, eps=0.4, 
%          rho_xr=0.2, rho_xv=-0.7, rho_rv=0.15, T=10, a=0.01613, sigma=0.012812
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    import a            "Hull-White mean reversion parameter"
    import sigma        "Hull-White volatility parameter" 
    import k            "Volatility mean reversion parameter"
    import tht          "Long term volatility"
    import eps          "Volatility of volatility"
    import v0           "Initial volatility"
    import S0           "Initial vtock price"
    import rho_sr       "Correlation between stock price and interest rate"
    import rho_rv       "Correlation between stock volatility and interest rate"
    import rho_sv       "Correlation between stock price and volatility" 
    import PtT          "Hull-White discount bond prices" 
    import T            "Discount bond maturity"
    import t            "Simulation time horizon"
    import dt           "Discretization time step"
    export StSZHW       "Schöbel-zhu Hull-White stock prices"
 
    %cholesky decomposition Lower correlation matrix components    
    rho_sv2 = sqrt(1-rho_sv^2);
    w_sr = (rho_sr-rho_sv*rho_rv)/rho_sv2;
    rho_rw = sqrt(1-rho_rv^2-w_sr^2);
 
    gam1 = 0.5;        %discretization parameter
    gam2 = 0.5;        %discretization parameter
 
    y0 = log(S0/PtT);  %initial log forward stock price
    yt = y0;           %initialize log forward stock price
    StSZHW = S0;       %initial SZHW stock price
    vt0 = v0;          %initial volatility     
    loop t/dt
        theta dt;
        %discretized Schöbel-Zhu stochastic volatility
        K1 = exp(-k*dt);
        K2 = (tht-rho_rv*sigma*eps/(a*k))*(1-exp(-k*dt))
        -rho_rv*sigma*eps/(a*(k+a))*(exp(-a*(T-@time)-k*dt)-exp(-a*(T-@time-dt)));
        K3 = eps*sqrt((1-exp(-2*k*dt))/(2*k));
        %independent standard normal variate for Schöbel-zhu volatility
        Zv = randn(); 
        %simulated Schöbel-Zhu volatility process            
        vt1 = K1*vt0+K2+K3*Zv;
 
        %%%%%%%%%%%discretizations for the log forward stock price process%%%%%%%%%%%%
        Htdt = sigma/a*(dt-exp(-a*(T-@time-dt))/a+exp(-a*(T-@time))/a);
        Gtdt = sigma^2/a^2*(dt+exp(-2*a*(T-@time-dt))/(2*a)-2*exp(-a*(T-@time-dt))/a
        -exp(-2*a*(T-@time))/(2*a)+2*exp(-a*(T-@time))/a);
 
        rho_vv2 = Htdt/sqrt(dt*Gtdt);
 
        C0 = -0.5*(Gtdt+rho_sv*eps*dt);
        C1 = -gam1*(rho_sr*Htdt+rho_sv*(tht*k*dt/eps-rho_rv*Htdt));
        C2 = -gam2*(rho_sr*Htdt+rho_sv*(tht*k*dt/eps-rho_rv*Htdt));
        C3 = -0.5*gam1*dt+rho_sv/eps*(gam1*k*dt-0.5);
        C4 = -0.5*gam2*dt+rho_sv/eps*(gam2*k*dt+0.5);
        C50 = w_sr^2*Gtdt;
        C51 = 2*gam1*w_sr*rho_sv2*Htdt;
        C52 = 2*gam2*w_sr*rho_sv2*Htdt;
        C53 = gam1*dt*rho_sv2^2;
        C54 = gam2*dt*rho_sv2^2;
        C5 = sqrt(C50+C51*vt0+C52*vt1+C53*vt0^2+C54*vt1^2);
        C6 = rho_rv*sqrt(Gtdt)*rho_vv2;
        C7 = rho_rv*sqrt(Gtdt)*sqrt(1-rho_vv2^2);
        C8 = rho_rw*sqrt(Gtdt);
 
        %do the Martingale correction for the discretized log asset price
        D1 = C1+0.5*C51-K1/K3*C6;
        D2 = C2+0.5*C52+C6/K3;
        D3 = C3+0.5*C53;
        D4 = C4+0.5*C54;
        %check whether the regularity condition is satisfied, if so, do Martingale correction 
        %for the drift
        RegularityCondition = 1-2*D4*K3^2;
 
        if RegularityCondition > 0
            E0 = 0.5*log(1-2*D4*K3^2)-D4*(K2+D2/(2*D4))^2/(1-2*D4*K3^2)
            -0.5*C50+C6*K2/K3-0.5*C7^2-0.5*C8^2+D2^2/(4*D4);
            E1 = -D1-2*D4*K1*(K2+D2/(2*D4))/(1-2*D4*K3^2);
            E2 = -D3-D4*K1^2/(1-2*D4*K3^2);
 
            %Martingale corrected drift
            C00 = E0+E1*vt0+E2*vt0^2;
        else
            C00 = C0;
        end
 
        %independent standard normal variates
        Zs = randn();        
        Zv2 = randn();       
        Zr = randn();        
        %simulated log forward stock price       
        yt = y0 + C00 + C1*vt0 + C2*vt1 + C3*vt0^2 + C4*vt1^2 + C5*Zs + C6*Zv + C7*Zv2 + C8*Zr;
 
        %update the log forward stock price 
        y0 = yt; 
        %update the Schöbel-zhu volatilities           
        vt0 = vt1;
 
        %simulated Schöbel-Zhu Hull-White One Factor stock prices
        StSZHW = exp(yt)*PtT; 
    end      
end


References

  1. 1.0 1.1 Haastrecht, van. A., Pelsser, A. A. J., 2009, Efficient, almost Exact Simulation of the Heston Stochastic Volatility Model, International Journal of Theoretical and Applied Finance, forthcoming.
  2. 2.0 2.1 2.2 2.3 Schöbel, R., Zhu, J., 1999, Stochastic Volatility with an Ornstein-Uhlenbeck Process: an Extension, European Finance Review 4, 23-46.


Appendix A. Derivations of equations (SZHW.7a) - (SZHW.7c)

To have the stock price and Schöbel-Zhu volatility dynamics under the forward measure $ \mathbb{Q}^T $, we use the Girsanov theorem for the measure change from $ \mathbb{Q} $ to $ \mathbb{Q}^T $.

Under the risk neutral measure $ \mathbb{Q} $, the numeraire asset bank account is defined as:

$ B(t) = e^{\int_0^{t} r(u) du}, \qquad(SZHW.A.1) $

with initial bank account value $ B(0) = 1 $. The term $ r(t) $ is a positive function of $ t $, which in this article is the Hull-White One Factor short rate.

Differentiation of equation (SZHW.A.1) gives the following dynamics for bank account process under the risk neutral measure $ \mathbb{Q} $:

$ d B(t) = r(t) B(t) dt. \qquad(SZHW.A.2) $

The numeraire asset for forward measure $ \mathbb{Q}^T $ is zero coupon bond of maturity $ T $. In this case it is zero coupon bond price $ P(t, \ T) $ computed from the Hull White One Factor Model. The dynamics for $ d P(t, \ T) $ is defined in equation (SZHW.5).

The Radon-Nikodym derivative $ d \Lambda_{\mathbb{Q}}^T (t) $ for the measure change from $ \mathbb{Q} $ to $ \mathbb{Q}^T $ is:

$ d \Lambda_{\mathbb{Q}}^T (t) = \frac{d \mathbb{Q}^T}{d \mathbb{Q}} \mid \mathcal{F}_0 $
$ = d \left ( \frac{\frac{P(t, \ T)}{P(0, \ T)}}{\frac{B(t)}{B(0)}} \right ) = d \left ( \frac{P(t, \ T)}{P(0, \ T)} \frac{1}{B(t)} \right ) $
$ = \frac{1}{P(0, \ T) B(t)} d P(t, \ T) \ + \ \frac{P(t, \ T)}{P(0, \ T)} d \left ( \frac{1}{B(t)} \right ). \qquad(SZHW.A.3) $

Substituting in the dynamics for $ d P(t, \ T) $ as in equation (SZHW.5), and noting that

$ d \left ( \frac{1}{B(t)} \right ) = - \frac{1}{B(t)^2} d B(t) = - \frac{1}{B(t)} r(t) dt. \qquad(SZHW.A.4) $

We have

$ d \Lambda_{\mathbb{Q}}^T (t) $
$ = \frac{P(t, \ T)}{P(0, \ T) B(t)} \left [ r(t) dt \ - \ \sigma B(t) d W_r (t) \right ] \ + \ \frac{P(t, \ T)}{P(0, \ T)} \left [ - \frac{1}{B(t)} r(t) dt \right ] $
$ = \frac{P(t, \ T)}{P(0, \ T) B(t)} \left [ - \sigma B(t) d W_r (t) \right ] $
$ = \Lambda_{\mathbb{Q}}^T (t) \left [ - \sigma B(t) d W_r (t) \right ] $
$ \Rightarrow \frac{d \Lambda_{\mathbb{Q}}^T (t)}{ \Lambda_{\mathbb{Q}}^T (t)} = - \sigma B(t) d W_r (t). \qquad(SZHW.A.5) $

In equation (SZHW.A.5), the third equality follows by replacing $ \Lambda_{\mathbb{Q}}^T (t) $ with its definition $ \frac{P(t, \ T)}{P(0, \ T) B(t)} $. Equation (SZHW.A.5) implies that the Girsanov kernel $ \eta_t $ for measure change from $ \mathbb{Q} $ to $ \mathbb{Q}^T $ is of the form $ - \sigma B(t) $, which is the zero coupon bond volatility for the process $ \frac{d P(t, \ T)}{P(t, \ T)} $ as shown in equation (SZHW.5).

Itó integration of equation (SZHW.A.5) from $ 0 $ to $ t $ gives $ \Lambda_{\mathbb{Q}}^T (t) $ as an exponential Martingale:

$ \Lambda_{\mathbb{Q}}^T (t) = e^{-\frac{1}{2} \int_0^t (\sigma B(u) )^2 du \ - \ \int_0^t \sigma B(u) d W_r(u)}. $

Define the corresponding standard Brownian motion for $ W_r (t) $ under the forward measure $ \mathbb{Q}^T $ as $ W_r^T (t) $, we have the following relation:

$ d W_r^T (t) = \sigma B(t) dt \ + \ d W_r (t). \qquad(SZHW.A.6) $

Thus the change of measure from $ \mathbb{Q} $ to $ \mathbb{Q}^T $ results in a drift change by the Girsanov kernel $ \eta_t $.

To have the representations of $ W_S^T (t) $ and $ W_v^T (t) $ in terms of $ W_S (t) $ and $ W_v (t) $, we note that the correlations among the driving Brownian motions defined in equation (SZHW.4), rearranged in the terms of correlations with $ d W_r^T (t) $ are as follows:

$ A = \begin{bmatrix} 1 \\ \rho_{rv} \\ \rho_{Sr} \\ \end{bmatrix}. \qquad(SZHW.A.7) $

Due to correlations among the Brownian motions, change of measure from $ \mathbb{Q} $ to $ \mathbb{Q}^T $ induces changes of drift in all the correlated processes. The drift change is the Girsanov kernel $ \eta_t $ scaled by the respective correlation with the standard Brownian motion driving the numeraire asset.

Specifically, under the forward measure $ \mathbb{Q}^T $ the Brownian motion driving the SchÖbel-Zhu volatility process as defined in (SZHW.2) is now

$ d W_v^T (t) = \rho_{rv} \sigma B(t) dt \ + \ d W_v (t), \qquad(SZHW.A.8) $

and the corresponding Brownian motion driving the stock price process as defined in (SZHW.1) is

$ d W_S^T (t) = \rho_{Sr} \sigma B(t) dt \ + \ d W_S (t). \qquad(SZHW.A.9) $