The Vasiček short-rate model

The Vasiček model is an interest rate model which specifies the short rate $r(t)$ under the risk-neutral dynamics (or $\mathbb{Q}$-dynamics) as $$\begin{array}{}\text{(1)}& dr(t)=\kappa (\theta -r(t))dt+\sigma dW(t),\end{array}$$ with initial condition $r(0)=r_0$ and $W(t)$ denoting a standard Brownian motion driving the stochastic differential equation. An explicit expression for $r(t)$ can be derived using Itô calculus (see, e.g. Mikosch (1998); Chapter 3). To solve $\text{(1)}$, we try the solution $f(r(t),t)=r(t)e^{\kappa t}$. The Itô lemma implies $$\begin{array}{}\text{(2)}& \begin{array}{rl}df(r(t),t)& =\kappa r(t)e^{\kappa t}dt+e^{\kappa t}dr(t)\\ & =\kappa r(t)e^{\kappa t}dt+e^{\kappa t}[\kappa (\theta -r(t))dt+\sigma dW(t)]\\ & =\kappa \theta e^{\kappa t}dt+\sigma e^{\kappa t}dW(t).\end{array}\end{array}$$ The RHS of $\text{(2)}$ no longer depends on $r(t)$ and we can thus integrate to find expressions for $f(r(t),t)$. This result is important for two reasons. First, we can integrate over the range $[0,t]$ to find $$\begin{array}{rl}f(r(t),t)& -f(r(0),0)=r(t)e^{\kappa t}-r_0={\int }_0^{t}df(r(s),s)\\ & ={\int }_0^t\kappa \theta e^{\kappa s}ds+{\int }_0^t\sigma e^{\kappa s}dW(s).\end{array}$$ After some rearrangments the explicit solution for $r(t)$ is $$\begin{array}{}\text{(3)}& r(t)=r_0{e}^{-\kappa t}+e^{-\kappa t}{\int }_0^t\kappa \theta e^{\kappa s}ds+e^{-\kappa t}{\int }_0^t\sigma e^{\kappa s}dW(s).\end{array}$$ With this explicit expression for $r(t)$ we can calculate its mean and variance. Since stochastic integrals have zero mean we have $$\begin{array}{rl}\mathbb{E}[r(t)]& =r_0{e}^{-\kappa t}+\theta e^{-\kappa t}[e^{\kappa t}-1]=r_0{e}^{-\kappa t}+\theta [1-e^{-\kappa t}]\\ & =:r_0{e}^{-\kappa t}+\theta \kappa \mathrm{\Lambda }(t),\end{array}$$ where we defined $\mathrm{\Lambda }(t)={\int }_0^t{e}^{-\kappa s}ds=\frac{1}{\kappa }[1-e^{-\kappa t}]$. The long-run mean of the process is $\underset{t\to \mathrm{\infty }}{lim}\mathbb{E}[r(t)]=\theta$. Moreover, by Itô isometry the variance of the process is $$\mathbb{V}\text{ar}[r(t)]={\sigma }^2{e}^{-2\kappa t}{\int }_0^t{e}^{2\kappa s}ds=\frac{{\sigma }^2}{2\kappa }[1-e^{-2\kappa t}]=\frac{1}{2}{\sigma }^2\mathrm{\Lambda }(2t).$$ Second, $\text{(2)}$ is the starting point to derive an exact discretization of $r(t)$. Such a discretization will allow us to simulate paths of $r(t)$ that can later be used to value interest rate derivatives. Given a time step $h$, we now integrate $\text{(2)}$ over the interval $[t-h,t]$ to obtain $$\begin{array}{rl}f& (r(t),t)-f(r(t-h),t-h)={\int }_{t-h}^{t}df(r(s),s)\\ & =\theta (e^{\kappa t}-e^{\kappa (t-h)})+{\int }_{t-h}^t\sigma e^{\kappa s}dW(s)\\ & ⟺r(t)=\theta (1-e^{-\kappa h})+e^{-\kappa h}r(t-h)+e^{-\kappa t}{\int }_{t-h}^t\sigma e^{\kappa s}dW(s).\end{array}$$ The discretized process of $r(t)$ thus follows an AR(1) model with intercept $\theta (1-e^{-\kappa h})$, autoregressive coefficient $e^{-\kappa h}$, and innovation process $e^{-\kappa t}{\int }_{t-h}^t\sigma e^{\kappa s}dW(s)$. These innovations have some interesting properties. First, note how the integration intervals of subsequent steps of the discretization do not overlap. For example, ${\int }_{t-h}^t\sigma e^{\kappa s}dW(s)$ depends on $\{W(s):t-h\le s\le t\}$, whereas ${\int }_t^{t+h}\sigma e^{\kappa s}dW(s)$ depends on $\{W(s):t\le s\le t+h\}$. With increments of Brownian motions being independent, we must conclude that these stochastic integrals are independent. Second, stochastic integrals are normally distributed and mean zero. The distributions of the innovations is thus fully specified after computing its variance, or $$\begin{array}{rl}\mathbb{V}& \text{ar}[e^{-\kappa t}{\int }_{t-h}^t\sigma e^{\kappa s}dW(s)]\\ & ={\sigma }^2{e}^{-2\kappa t}{\int }_{t-h}^t{e}^{2\kappa s}ds=\frac{{\sigma }^2}{2\kappa }[1-e^{-2\kappa h}]={\sigma }^{2}h\left[\frac{e^{-2\kappa h}-1}{-2\kappa h}\right]\\ & =:{\sigma }^{2}h\alpha (-2\kappa h),\end{array}$$ where we defined $\alpha (x)=\frac{e^x-1}{x}\text{ with }\alpha (0)=0$. Overall, we can simulate data from the Vasiček model using $r(0)=r_0$ as the starting value and moving forward according to the recursion $$r(t)=\theta \kappa \mathrm{\Lambda }(h)+e^{-\kappa h}r(t-h)+u_t^{(h)},u_t^{(h)}\stackrel{i.i.d.}{\sim }\mathcal{N}(0,{\sigma }^{2}h\alpha (-2\kappa h)).$$ Some sample paths of the Vasiček model are displayed in Figure 1.

Risk-neutral valuation of a zero-coupon bond

Let $X$ denote a contingent claim with maturity date $T$. According to the risk-neutral valuation formula (cf. Proposition 8.1.2 in Bingham and Kiesel (2004)), the price at time $t$ of this claim can be computed as ${\mathrm{\Pi }}_X(t)={\mathbb{E}}_{\mathbb{Q}}\left[X{e}^{-{\int }_t^{T}r(s)ds}\right|{\mathcal{F}}_t]$. Since the zero-coupon bond, or $T$-bond, promises a cash payment of 1 at maturity, its time-$t$ price is given by

$$\begin{array}{}\text{(4)}& P(t,T)={\mathbb{E}}_{\mathbb{Q}}\left[e^{-{\int }_t^{T}r(s)ds}\right|{\mathcal{F}}_t],\end{array}$$ where ${\mathcal{F}}_t$ is the $\sigma$-algebra containing all information up to time $t$. In this section we will evaluated $\text{(4)}$ analytically. That is, starting from the explicit expression for $r(t)$ in $\text{(3)}$, we first derive the distribution of ${\int }_t^{T}r(s)ds|{\mathcal{F}}_t$ and subsequently evaluated the conditional expectation. The first step is rather tedious and explained in detail in the Appendix. It turns out that ${\int }_t^{T}r(s)ds|{\mathcal{F}}_t$ is normally distributed with mean $$r(t)\mathrm{\Lambda }(T-t)+\theta [(T-t)-\mathrm{\Lambda }(T-t)]$$ and variance $$\frac{{\sigma }^2}{{\kappa }^2}[(T-t)-2\mathrm{\Lambda }(T-t)+\frac{1}{2}\mathrm{\Lambda }(2(T-t))].$$ The expectation in $\text{(4)}$ can be calculated rather quickly using moment generating functions. For a random variable $Y$, its moment generating function is defined as $M_Y(t)=\mathbb{E}\left[e^{tY}\right]$. If $Y$ is normally distributed, say $Y\sim \mathcal{N}(\mu ,{\sigma }^2)$, then $M_Y(t)=e^{\mu t+\frac{1}{2}{\sigma }^2{t}^2}$. To evaluate $\text{(4)}$, we use this result for $t=-1$ and find that the time-$t$ price of a zero-coupon bond with maturity $T$ equals

$$\begin{array}{}\text{(5)}& \begin{array}{rl}P(t,T)& =\exp \{-r(t)\mathrm{\Lambda }(T-t)-\theta [(T-t)-\mathrm{\Lambda }(T-t)]+\frac{{\sigma }^2}{2{\kappa }^2}[(T-t)-2\mathrm{\Lambda }(T-t)+\frac{1}{2}\mathrm{\Lambda }(2(T-t))]\}\\ & =\exp \{-r(t)\mathrm{\Lambda }(T-t)+(\frac{{\sigma }^2}{2{\kappa }^2}-\theta )(T-t)+(\theta -\frac{{\sigma }^2}{{\kappa }^2})\mathrm{\Lambda }(T-t)+\frac{{\sigma }^2}{4{\kappa }^2}\mathrm{\Lambda }(2(T-t))\}.\end{array}\end{array}$$

We can translate these zero-coupon bond prices into yields using $$y(t,T)=-\frac{1}{T-t}\log P(t,T)=(\theta -\frac{{\sigma }^2}{2{\kappa }^2})+\frac{1}{T-t}[r(t)\mathrm{\Lambda }(T-t)+(\frac{{\sigma }^2}{{\kappa }^2}-\theta )\mathrm{\Lambda }(T-t)-\frac{{\sigma }^2}{4{\kappa }^2}\mathrm{\Lambda }(2(T-t))].$$ At long maturities, as $T\to \mathrm{\infty }$, the yield converges to $\theta -\frac{{\sigma }^2}{2{\kappa }^2}$. Visualisations of the complete yield curve are shown in the section entitled Verification by Monte Carlo simulation.

Feynman-Kac formula: solving the PDE

Consider a short-rate model with $\mathbb{Q}$-dynamics given by $$dr(t)=a(t,r(t))dt+b(t,r(t))dW(t)$$ and write $P(t,T)=F(t,r(t);T)$ to explicitly indicate the dependence on $r(t)$. For brevity, we will sometimes in this section omit the function arguments, e.g. write $F$ instead of $F(t,r(t);T)$. The Feynman-Kac formula (see, e.g. Bingham and Kiesel (2004); Proposition 8.2.2) stipulates that $F(t,r(t);T)$ solves the partial differential equation (PDE)

$$\begin{array}{}\text{(6)}& \frac{\partial F}{\partial t}+a\frac{\partial F}{\partial r}+\frac{b^2}{2}\frac{{\partial }^{2}F}{\partial r^2}-rF=0,\end{array}$$

with terminal condition $F(T,r;T)=1$ for all $r\in \mathbb{R}$. We make two observations. First, we have $a(t,r(t))=\kappa (\theta -r(t))$ and $b(t,r(t))=\sigma$ for the Vasiček model. Second, it is hard (or sometimes even impossible) to solve $\text{(6)}$ analytically. We are however in the lucky situation where $a(t,r(t))$ and $b(t,r(t))$ are linear in $r(t)$. It can be shown (cf. Filipović (2009); Proposition 5.2) that this leads to an affine term structure, that is the solution $F(t,r;T)$ must take the form $$F(t,r;T)=\exp [-A(t,T)-B(t,T)r],$$ for appropriate $A(t,T)$ and $B(t,T)$.

We can now solve the PDE by inserting this specific functional form into $\text{(6)}$ and see what this implies for $A$ and $B$. Since $\frac{\partial F}{\partial t}=(-\frac{\partial A}{\partial t}-\frac{\partial B}{\partial t}r)F$, $\frac{\partial F}{\partial r}=-BF$ and $\frac{{\partial }^{2}F}{\partial r^2}=B^{2}F$, we find $$(-\frac{\partial A}{\partial t}-\frac{\partial B}{\partial t}r)F-\kappa (\theta -r)BF+\frac{{\sigma }^2}{2}{B}^{2}F-rF=0$$ or equivalently after collecting terms $$F[(-\frac{\partial A}{\partial t}-\kappa \theta B+\frac{{\sigma }^2}{2}{B}^2)+(-\frac{\partial B}{\partial t}+\kappa B-1)r]=0.$$ The boundary condition is $F(T,r;T)=\exp [-A(T,T)-B(T,T)r]=1$ or $A(T,T)+B(T,T)r=0$. If these relations need to hold for all $r\in \mathbb{R}$, then intercept terms should be zero as well as the expressions proportional to $r$. The PDE for $F(T,r(t);T)$ is now seen to reduce to a set of coupled ordinary differential equations (ODEs):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& -\frac{\partial A}{\partial t}-\kappa \theta B+\frac{{\sigma }^2}{2}{B}^2=0,\\ & -\frac{\partial B}{\partial t}+\kappa B-1=0,\\ & A(T,T)=B(T,T)=0.\end{array}\end{array}$$

The second equation in $\text{(7)}$ completely specifies $B$. With some rewriting, we have $\frac{\partial B}{\partial t}-\kappa B=e^{\kappa t}\frac{\partial }{\partial t}\left[e^{-\kappa t}B\right]=-1$. Integrating over $[t,T]$ gives $$e^{-\kappa T}B(T,T)-e^{-\kappa t}B(t,T)=-{\int }_t^T{e}^{-\kappa s}ds.$$ Together with the boundary condition $B(T,T)=0$ we conclude that $B(t,T)={\int }_t^T{e}^{-\kappa (s-t)}ds=\frac{1}{\kappa }[1-e^{-\kappa (T-t)}]=\mathrm{\Lambda }(T-t)$. Having found the explicit solution for $B(t,T)$, we can complete the derivations by $$\begin{array}{rl}& A(t,T)=-[A(T,T)-A(t,T)]={\int }_t^T-\frac{\partial A}{\partial s}ds\\ & ={\int }_t^T[\kappa \theta B(s,T)-\frac{{\sigma }^2}{2}[B(s,T){]}^2]ds\\ & =\theta [(T-t)-{\int }_t^T{e}^{-\kappa (T-s)}ds]-\frac{{\sigma }^2}{2{\kappa }^2}{\int }_t^T[1-2e^{-\kappa (T-s)}+e^{-2\kappa (T-s)}]ds,\\ & =\theta [(T-t)-\mathrm{\Lambda }(T-t)]-\frac{{\sigma }^2}{2{\kappa }^2}[(T-t)-2\mathrm{\Lambda }(T-t)+\frac{1}{2}\mathrm{\Lambda }(2(T-t))]\\ & =(\theta -\frac{{\sigma }^2}{2{\kappa }^2})(T-t)+(\frac{{\sigma }^2}{{\kappa }^2}-\theta )\mathrm{\Lambda }(T-t)-\frac{{\sigma }^2}{4{\kappa }^2}\mathrm{\Lambda }(2(T-t)),\end{array}$$ where we used ${\int }_t^T{e}^{-\kappa (T-s)}ds=\frac{1}{\kappa }[1-e^{-\kappa (T-t)}]=\mathrm{\Lambda }(T-t)$ and ${\int }_t^T{e}^{-2\kappa (T-s)}ds=\frac{1}{2}\mathrm{\Lambda }(2(T-t))$. The overall expression for $P(t,T)=F(t,r(t);T)$ coincides with the result from the previous section.

Verification by Monte Carlo simulation

If we like to avoid extensive algebraic computations, then we can opt for a simulation approach. That is, we approximate the expectation in $\text{(4)}$ by Monte Carlo simulation. Our example is $P(0,T)={\mathbb{E}}_{\mathbb{Q}}[e^{-{\int }_0^{T}r(s)ds}]$. The steps are as follows:

1. Partition the interval $[0,T]$ in $n$ intervals of equal length. For $j=0,1,\dots ,n$, the implied grid points are $t_j=j\frac{T}{n}$.
2. Start from $r(0)=r(t_0)=r_0$ and use the AR(1) recursion with stepsize $h=\frac{T}{n}$ to simulate $N_{sim}$ sample paths of the Vasiček model. We use $r^{(i)}(t_j)$ to denote the realised value of the $i^{\text{th}}$ path at grid point $t_j$.
3. Approximate $P(0,T)$ by $$\begin{array}{}\text{(8)}& \frac{1}{N_{sim}}\sum _{i=1}^{N_{sim}}{e}^{-h\sum _{j=1}^n{r}^{(i)}(t_j)}.\end{array}$$

The computation of $\text{(8)}$ requires choices for $n$ and $N_{sim}$. We like both these quantities to be large. A large number of simulated paths is needed because the expectation in $\mathbb{E}{\mathbb{Q}} \big[e^{-\int_0^T r(s) ds} \big]$isbeingreplacedbyasampleaverageacrossthesimulatedpaths.Thelargevaluefor$n$shouldensurethattheintegral$\int_0^T r(s) ds$issufficientlywellapproximatedbytheRiemannsum$\sum{j=1}^n r^{(i)}(t_j)\big[ t_j - t_{j-1} \big] = h \sum_{j=1}^n r^{(i)}(t_j)$.

Visual evidence for this simulated approach is available in Figures 2–3. We have taken $r_0=6$, $\theta =0.08$, $\kappa =0.86$, and $\sigma =0.01$. The (estimated) zero-coupon bond prices from $\text{(8)}$ are converted into yields. Figure 2 shows that: (1) differences between simulated and exact yields are small, and (2) variability between simulated yields decreases when $N_{sim}$ increases from 50 to 500. The influence of $n$ is portrayed in Figure 3. In practice, we can use these kind of graphs to decide on suitable choices for $n$ and $N_{sim}$. Simply select a pair $(n,N_{sim})$ and verify whether the computed quantity is insensitive to changes therein.

Appendix

Recall $\mathrm{\Lambda }(t)={\int }_0^t{e}^{-\kappa s}ds=\frac{1}{\kappa }[1-e^{-\kappa t}]$ and note how it implies $$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\int }_x^T{e}^{-\kappa (s-x)}ds& =e^{\kappa x}[\mathrm{\Lambda }(T)-\mathrm{\Lambda }(x)]=\frac{1}{\kappa }{e}^{\kappa x}[e^{-\kappa x}-e^{-\kappa T}]\\ & =\mathrm{\Lambda }(T-x).\end{array}\end{array}$$ The integral in this last equation will be used at several occasions in the derivations below. Using the expression for $r(t)$, we have $$\begin{array}{rl}{\int }_t^{T}r(s)ds& ={\int }_t^T{e}^{-\kappa s}r(0)ds\\ & +{\int }_t^T{\int }_0^s{e}^{-\kappa (s-x)}\kappa \theta dxds+{\int }_t^T{\int }_0^s{e}^{-\kappa (s-x)}\sigma dW(x)ds\\ & =:I+II+III.\end{array}$$ We will develop these three contributions separately. Term $I$ is easiest. Using $\text{(A1)}$, we find $$I={\int }_t^T{e}^{-\kappa s}r(0)ds=r(0)e^{-\kappa t}{\int }_t^T{e}^{-\kappa (s-t)}ds=r(0)e^{-\kappa t}\mathrm{\Lambda }(T-t).$$

Term $II$ requires a change in the order of integration. Inspecting the area of integration in Figure 4, we arrive at the following integral relation $$\begin{array}{rl}II& ={\int }_t^T{\int }_0^s{e}^{-\kappa (s-x)}\kappa \theta dxds\\ & ={\int }_0^t{\int }_t^T{e}^{-\kappa (s-x)}\kappa \theta dsdx+{\int }_t^T{\int }_x^T{e}^{-\kappa (s-x)}\kappa \theta dsdx\\ & ={\int }_0^t{\int }_t^T{e}^{-\kappa (s-t+t-x)}\kappa \theta dsdx+\kappa \theta {\int }_t^T\mathrm{\Lambda }(T-x)dx\\ & ={\int }_0^t{e}^{-\kappa (t-x)}\kappa \theta dx,\mathrm{\Lambda }(T-t)+\kappa \theta {\int }_t^T\mathrm{\Lambda }(T-x)dx.\end{array}$$ We apply exactly the same change in the order of integration to term $III$. The result is $$\begin{array}{rl}III& ={\int }_t^T{\int }_0^s{e}^{-\kappa (s-x)}\sigma dW(x)ds\\ & ={\int }_0^t{\int }_t^T{e}^{-\kappa (s-x)}\sigma dsdW(x)+{\int }_t^T{\int }_x^T{e}^{-\kappa (s-x)}\sigma dsdW(x)\\ & ={\int }_0^t{\int }_t^T{e}^{-\kappa (s-t+t-x)}\sigma dsdW(x)+\sigma {\int }_t^T\mathrm{\Lambda }(T-x)dW(x)\\ & ={\int }_0^t{e}^{-\kappa (t-x)}\sigma dW(x),\mathrm{\Lambda }(T-t)+\sigma {\int }_t^T\mathrm{\Lambda }(T-x)dW(x).\end{array}$$ If terms $I$–$III$ are added together, then we see that ${\int }_x^T{e}^{-\kappa (s-x)}ds$ is an affine transformation of the prevailing short rate $r(t)$, i.e. $$\begin{array}{rl}{\int }_t^T& r(s)ds\\ & =I+II+III\\ & =[r(0)e^{-\kappa t}+{\int }_0^t{e}^{-\kappa (t-x)}\kappa \theta dx+{\int }_0^t{e}^{-\kappa (t-x)}\sigma dW(x)]\mathrm{\Lambda }(T-t)\\ & +\kappa \theta {\int }_t^T\mathrm{\Lambda }(T-x)dx+\sigma {\int }_t^T\mathrm{\Lambda }(T-x)dW(x)\\ & =r(t)\mathrm{\Lambda }(T-t)+\kappa \theta {\int }_t^T\mathrm{\Lambda }(T-x)dx+\sigma {\int }_t^T\mathrm{\Lambda }(T-x)dW(x).\end{array}$$ Conditional on all information up to time t, i.e. conditional on ${\mathcal{F}}_t$, the first two terms are deterministic. Moreover, since (1) increments of Brownian motions are independent of the current value and (2) stochastic integrals are normally distributed, we know that ${\int }_t^{T}r(s)ds|{\mathcal{F}}_t$ is normally distributed with mean $$\begin{array}{rl}r(t)\mathrm{\Lambda }(T-t)& +\kappa \theta {\int }_t^T\mathrm{\Lambda }(T-x)dx\\ & =r(t)\mathrm{\Lambda }(T-t)+\kappa \theta {\int }_0^{T-t}\mathrm{\Lambda }(x)dx\\ & =r(t)\mathrm{\Lambda }(T-t)+\theta [(T-t)-\mathrm{\Lambda }(T-t)]\end{array}$$ and by Itô isometry a variance of $$\begin{array}{rl}{\sigma }^2{\int }_t^T& [\mathrm{\Lambda }(T-x){]}^{2}dx={\sigma }^2{\int }_0^{T-t}[\mathrm{\Lambda }(x){]}^{2}dx\\ & =\frac{{\sigma }^2}{{\kappa }^2}{\int }_0^{T-t}(1-2e^{-\kappa x}+e^{-2\kappa x})dx\\ & =\frac{{\sigma }^2}{{\kappa }^2}[(T-t)-2{\int }_0^{T-t}{e}^{-\kappa x}dx+\frac{1}{2}{\int }_0^{2(T-t)}{e}^{-\kappa x}dx]\\ & =\frac{{\sigma }^2}{{\kappa }^2}[(T-t)-2\mathrm{\Lambda }(T-t)+\frac{1}{2}\mathrm{\Lambda }(2(T-t))].\end{array}$$

References

N. H. Bingham and R. Kiesel (2004), Risk-neutral Valuation, Springer Finance

D. Filipović (2009), Term-structure Models, Springer Finance

T. Mikosch (1998), Elementary Stochastic Calculus with Finance in View, World Scientific