Exact solution stock price with Vasicek interest rate model

Question

Define two correlated stock price- and interest rate (Vasicek) processes, governed by the Wiener processes $W^{S}(t)$ and $W^{r}(t)$

$$dS(t)=r(t)S(t)dt+\sigma S(t)dW^{S}(t)$$

$$dr(t)=\kappa(\theta-r(t))dt+\gamma dW^{r}(t)$$

with constant scalars $S_{0}>0$, $T>0$, $r_{0}>0$, $\sigma>0$, $\theta>0$, $\gamma>0$, $\kappa>0$, and $t>0$, $dW^{S}(t)dW^{r}(t)=\rho dt$ with $t\in[0,T]$.

The exact solution to the stock price at time $T$ is as following

$$S(T)=S_{0}\exp(\int^{T}_{0}r(s)ds-\frac{1}{2}\sigma^{2}T+\sigma W^{S}(T))$$

By drawing $N$ times from $W(T)\sim\mathcal{N}(0,T)$ an approximation of the expected value can be made through a Monte Carlo simulation; however, the term $\int^{T}_{0}r(s)ds$ is stochastic, since the exact solution for $r(s)$ for the Vasicek model is as following

$$r(s)=r_{0}e^{-\kappa s}+\theta(1-e^{-\kappa s})+\gamma e^{-\kappa s}\int^{s}_{0}e^{\kappa\bar{s}}dW^{r}(\bar{s})$$

with $\bar{s}$ a dummy variable and $\int^{s}_{0}e^{\kappa\bar{s}}dW^{r}(\bar{s})\sim\mathcal{N}(0,\frac{1}{2\kappa}(e^{2\kappa s}-1))$ according to Ito's lemma, therefore

$$\int^{T}_{0}r(s)ds=\frac{r_{0}}{\kappa}(1-e^{\kappa T})+\theta T+\frac{\theta}{\kappa}(e^{\kappa T}-1)+\int^{T}_{0}\gamma e^{-\kappa s}\sqrt{\frac{1}{2\kappa}(e^{2\kappa s}-1)}Zds$$

with $Z\sim\mathcal{N}(0,1)$. Next, using integration by parts and substitution the integral $\int^{T}_{0}\gamma e^{-\kappa s}\sqrt{\frac{1}{2\kappa}(e^{2\kappa s}-1)}$ can be solved by firstly choosing $dU=\gamma e^{-\kappa s}$ and $V=\sqrt{\frac{1}{2\kappa}(e^{2\kappa s}-1)}$, which leads to

$$\int VdU=UV-\int UdV\implies\int^{T}_{0}\gamma e^{-\kappa s}\sqrt{\frac{1}{2\kappa}(e^{2\kappa s}-1)}ds=\frac{-\gamma e^{-\kappa T}}{\kappa}\sqrt{\frac{1}{2\kappa}(e^{2\kappa T}-1)}-\frac{\gamma}{2\kappa}\int^{T}_{0}\frac{1}{\sqrt{e^{\kappa s}-1}}ds$$

and secondly choosing $x=e^{\kappa s}-1$ with $ds=\frac{1}{\kappa x}dx$ leading to

$$\frac{-\gamma e^{-\kappa T}}{\kappa}\sqrt{\frac{1}{2\kappa}(e^{2\kappa T}-1)}-\frac{\gamma}{2\kappa}\int^{T}_{0}\frac{1}{\sqrt{e^{\kappa s}-1}}ds=\frac{-\gamma e^{-\kappa T}}{\kappa}\sqrt{\frac{1}{2\kappa}(e^{2\kappa T}-1)}+\frac{\gamma}{\kappa^{2}}(e^{\frac{-\kappa T}{2}}-1)$$

which is deterministic. Finally, a Monte Carlo simulation can be done by drawing $N$ times from $\int^{T}_{0}r(s)ds\sim\mathcal{N}(\mu=\frac{r_{0}-\theta}{\kappa}(1-e^{\kappa T})+\theta T,\sigma=\frac{-\gamma e^{-\kappa T}}{\kappa}\sqrt{\frac{1}{2\kappa}(e^{2\kappa T}-1)}+\frac{\gamma}{\kappa^{2}}(e^{\frac{-\kappa T}{2}}-1)$ to approximate the expectation of the exact solution of $S(T)$.

My questions are: i) is the above reasoning correct, where I write the integral with respect to $dW^{r}(t)$ as an integral with respect to $ds$ multiplied by a standard normal random variable, $Z$ and ii) since W^{r}(T) is now governed by $Z\sim\mathcal{N}(0,1)$ can I still compute $W^{S}(T)$ as following

$$W^{S}(T)=\sqrt{T}(\rho Z+\sqrt{1-\rho^{2}\bar{Z}})$$

with again $\bar{Z}\sim\mathcal{N}(0,1)$? (This is an uncoventional notation, but convenient in my opinion.)

Answer 1

Your approach seems to be fine, however, I think it is worth considering a similar way to obtain the price of a vanilla option. Let us consider your option dynamics:

$$dS(t)=r(t)S(t)dt+\sigma S(t)dW^{S}(t)$$ $$dr(t)=\kappa(\theta-r(t))dt+\gamma dW^{r}(t)$$

Let us now compute the process under the risk-neutral measure using two independent Brownian motions: $$dr(t)=\kappa(\theta-r(t))dt+\gamma d\bar{W}^{r}(t)$$ $$\frac{dS(t)}{S(t)}=r(t)dt+\rho\sigma \bar{W}^{r}(t)+\sigma\sqrt{1-\rho^2}\bar{W}^{S}.$$ Then, let us apply a change of measure using the ZCB $P(t,T)$ as the numéraire, which follows Vasicek's dynamics. $$d\kappa^T_\mathbb{Q}(t):=\frac{d\mathbb{Q}^T}{d\mathbb{Q}}|_{\mathcal{F}(t)}=\frac{P(t,T)}{P(0,T)}\frac{N(0)}{N(t)},$$ and by Ito's lemma: $$d\kappa^T_\mathbb{Q}(t)=\frac{1}{N(t)}dP(t,T)-\frac{P(t,T)}{N(t)^2}dN(t)\Leftrightarrow\frac{d\kappa^T_\mathbb{Q}(t)}{\kappa^T_\mathbb{Q}(t)}=\gamma B_r(t,T)d\bar{W}^r(t).$$ Thus, we obtain the $T$-forward measure $$d\bar{W}^r(t)=\gamma B_r(t,T)dt+d\bar{W}^{r,T}(t),$$ $$d\bar{W}^S(t)=d\bar{W}^{S,T}(t),$$ so that now we can price the ZCB under this new measure: $$\frac{P(t,T)}{P(t,T)}=r(t)dt+\gamma B(t,T)\left(\gamma B(t,T)dt+d\bar{W}^{r,T}(t)\right)$$ $$=\left(r(t)+\gamma^2B^2(t,T)\right)dt+\gamma B(t,T)d\bar{W}^{r,T}(t).$$ The stock dynamics are then: $$\frac{dS(t)}{S(t)}=\left(r(t)+\rho\gamma\sigma B(t,T)\right)dt+\sigma\left(\rho d\bar{W}^{r,T}(t)+\sqrt{1-\rho^2}d\bar{W}^{S,T}(t)\right),$$ with $B(t,T)=\frac{1}{\kappa}\left(e^{-\kappa(T-t)}-1\right)$ and $d\bar{W}^{S,T}(t)d\bar{W}^{r,T}(t)=\rho dt$. Finally, we use this for pricing the option by using the fact that the forward stock price $Z(t,T):=\frac{S(t)}{P(t,T)}$ is a $\mathbb{Q}^T$-martingale. Let us express the option value in terms of the expectation (in terms of the risk-neutral measure at time $T$) of its payoff function $H(\cdot)$: $$V(t_0,S)=N(t_0)\mathbb{E}^\mathbb{Q}\left[\frac{1}{N(T)}H(T,S)|\mathcal{F}(t_0\right]=P(t_0,T)\mathbb{E}^T[H(T,S)|\mathcal{F}(t_0)].$$ Again, by Ito's lemma: $$dZ(t,T)=\frac{1}{P(t,T)}dS(t)-\frac{S(t)}{P(t,T)}dP(t,T)+\frac{S(t)}{P^3(t,T)}(dP(t,T))^2-\frac{1}{P^2(t,T)}dP(t,T)dS(t)$$ $$\Leftrightarrow \frac{dZ(t,T)}{Z(t,T)}=\sigma dW^{S,T}(t)-\gamma B(t,T)dW^{r,T}(t).$$ It is clearly a martingale, moreover we can re-express the two Brownian motions: $$dZ(t,T)=\sigma_Z(t)dW^Z(t),$$ where $\sigma_Z(t)=\sqrt{\sigma^2+\gamma^2B(t,T)^2-2\rho\sigma\gamma B(t,T)}$. Note that the pricing formula does not depend on rates, but its volatility is time-varying. This means that this approach works not only for Vasicek, but also for other interest rate models and it is quite convenient because $Z(t,T)=S(t)/P(t,T)$ does not have $r(t)$ (since it is driftless) and $Z(T,T)=S(T)$. We are now only left with computing the discounted payoff: $$V(t_0,S)=P(t_0,T)\mathbb{E}^T[(Z(T,T)-K)^+|\mathcal{F}(t_0)]=Z(0,T)P(t_0,T)\Phi(d_+)-KP(t_0,T)\Phi(d_-)$$

$$d_+ =\frac{\log\left(\frac{Z(0,T)}{K}\right)+\frac{1}{2}\sigma_T^2(T-t_0)}{\sigma_T\sqrt{T-t_0}}$$ $$d_- =\frac{\log\left(\frac{Z(0,T)}{K}\right)-\frac{1}{2}\sigma_T^2(T-t_0)}{\sigma_T\sqrt{T-t_0}}$$ $$\sigma_T^2=\frac{1}{T-t_0}\int_{t_0}^T\sigma^2_Z(s)ds.$$