Sampling change in the driving brownian motion of a CIR process

Question

I have volatility driven by a CIR process:

$$\mathrm{d}v_t = \kappa (\bar{v}-v_t)\mathrm{d}t + \omega \sqrt{v_t}\mathrm{d}W_v\text{.}\tag{1}$$

I am working with several (complicated) approximations of this process (for example, QE from the Andersen paper). Given $v(t)$, these approximations sample a $v(t+\epsilon)$. However, in addition to having a $v(t+\epsilon)$, I would like to sample

$$\int_t^{t+\epsilon} \mathrm{d}W_v = W_v(t+\epsilon)-W_v(t)\text{.}\tag{2}$$

Question: Given $v(t)$ and $v(t+\epsilon)$, how can I sample from the conditional distribution of (2)?

Answer 1

One can use the Euler-Maruyama discretization scheme for CIR, 'fixed' for $v$ positivity, to get:

$$ v(t+\epsilon) -v(t)\approx \kappa (\bar{v} -v(t)^+)\epsilon + \omega \sqrt{v(t)^+} (W_v(t+\epsilon) - W_v(t)). $$

So, one approximation of the Brownian increment, when $v(t)$ and $v(t+\epsilon)$ are given, is:

$$ W_v(t+\epsilon) - W_v(t) \approx \frac{v(t+\epsilon) -v(t) - \kappa (\bar{v} -v(t)^+)\epsilon}{\omega \sqrt{v(t)^+} } \;\;\;\;\;({\rm when} \; v(t)\not= 0)$$

Note: In the Heston model context, one usually gets rid of the integral of $\sqrt{v(t)}dW_v(t) $ (integral against $W_v$) using the exact equality (equation (10), page 7 in Andersen's paper):

$$ \int_t^{t+\epsilon} \sqrt{v(u)}dW_v(u) = \omega^{-1} \left(v(t+\epsilon) -v(t) - \kappa \bar{v} \epsilon - \kappa\int_t^{t+\epsilon} v(u)du \right),$$

after employing Cholesky decomposition on $W_X$, leaving to compute an integral against a new Brownian motion $W$ that is independent of $v$, $\int_t^{t+\epsilon} \sqrt{v(u)}dW(u)$.

Answer 2

Edit: This is probably incorrect.

The Quadratic Exponential scheme is the best one I have seen as it converges in distribution and is pretty fast, so nice choice there!

When $\eta$ is constant you can simplify the integral $$ \int_t^{t+\varepsilon}\eta dW(u)=\eta\int_t^{t+\varepsilon}dW(u)=\eta\left(W(t+\varepsilon)-W(t)\right) $$

In the QE scheme you either use a standard uniform variable or a standard normal variable. Denote them $U_V$ and $Z_V$, respectively. We know that changes in a Wiener process is normally distributed as follows $$ \Delta W\equiv W(t+\varepsilon)-W(t)\sim \mathcal{N}(0,\varepsilon) $$ where the second argument is the variance. So in the simulation we can find $$ \Delta W=\begin{cases} \sqrt{\varepsilon}\cdot Z_V&\text{if }\psi\leq\psi_c\\ \sqrt{\varepsilon}\cdot\Phi^{-1}(U_V)&\text{if }\psi>\psi_c \end{cases} $$