# Exact solution stock price with Vasicek interest rate model

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.)

$$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.$$