# Bond-pricing under the Vasicek short rate model

I'm currently studying the Vasicek model of the short interest rate

$$dr_t=a(\mu-r_t)dt+\sigma dW_t$$

I know how to solve this stochastic differential equation (SDE) and how to find expectation and variance of $$r_t$$. Then I wanted to find the function to describe the evolution of the price $$B(r_t,t)$$ of a zero-coupon bond. I've seen you can use Ito's formula to obtain this differential equation:

$$\frac{\partial B}{\partial t}+\frac{\sigma^2}{2}\frac{\partial^2 B}{\partial r^2}+(a(\mu-r)-\lambda\sigma)\frac{\partial B}{\partial r}-rB=0 \tag{1}$$

where $$\lambda$$ is the market price of risk (for reference, check pages 391-392 of Yue-Kuen Kwok's Mathematical Models of Financial Derivatives [PDF]). Other articles give this equation (sometimes considering $$\lambda=0$$) and some give the solution in a closed form. Up to here I'm okay.

Then I need the Green's function for this equation so I saw that the Vasicek model is a particular Ornstein–Uhlenbeck process with an additional drift term: the classic Ornstein–Uhlenbeck process $$dr_t=-ar_tdt+\sigma dW_t$$ can also be described in terms of a probability density function, $$P(r,t)$$, which specifies the probability of finding the process in the state $$r$$ at time $$t$$. This function satisfies the Fokker–Planck equation

$$\frac{\partial P}{\partial t}=\frac{\sigma^2}{2}\frac{\partial^2 P}{\partial r^2}+a\frac{\partial (rP)}{\partial r} \tag{2}$$

The transition probability, also known as the Green's function, $$P(r,t\mid r',t')$$ is a Gaussian with mean $$r'e^{-a(t-t')}$$ and variance $$\frac {\sigma^2}{2a}\left(1-e^{-2a(t-t')}\right)$$:

$$P(r,t\mid r',t')={\sqrt {\frac {a }{\pi \sigma^2(1-e^{-2a (t-t')})}}}\exp \left[-{\frac {a}{\sigma^2}}{\frac {(r-r'e^{-a (t-t')})^{2}}{1-e^{-2a (t-t')}}}\right] \tag{3}$$

This gives the probability of the state $$r$$ occurring at time $$t$$ given initial state $$r'$$ at time $$t′. Equivalently, $$P(r,t\mid r',t')$$ is the solution of the Fokker–Planck equation with initial condition $$P(r,t')=\delta(r-r')$$.

My aim is to test some numerical methods on this model in order to extend them on the CIR model later so I need the Green's function of this Vasicek model and the corrisponding differential equation (if equation (1) is not correct).

My try

I tried to correlate equation (1) and (2) by adding the missing drift term to the O-U process and considering $$\lambda=0$$ in (1) but I get $$-aP$$ in (2) and not $$-rP$$ as it is in (1) (also the signs are misplaced). Then I thought that maybe I should try to correlate not the forward equation (2) to (1) but the backward Kolmogorov equation (which in this case is exactly equation (1) but without the term $$-rB$$). However that would require to get rid of the term $$-rB$$ in (1) but I don't think this is possible since $$B$$ is a function of $$r$$. This is why I think correlating equation (1) with the equation (2) or his bacward Kolmogorov version is not possible.

Second attempt was then changing the Green's function according to the new O-U process, the one that matched the Vasicek model (the term $$r'e^{-a(t-t')}$$ is changed to the expected value of the Vasicek model $$\mu+[r'-\mu]e^{-a(t-t')}$$), and since this solves backward Kolmogorov, which is (1) without the term $$-rB$$, maybe I can just adjust this by a multiplying factor so that this solves (1) too.

The reasons I have are that:

1. I checked on MATLAB the new Green's function and it seems to solve Fokker-Planck forward and Kolmogorov backward; also gives $$1$$ when integrated in $$(r,t)\in\mathbb{R}\times[0,1]$$ (with $$r'=r_0$$ and $$t'=0$$) and in $$r\in\mathbb{R}$$ (with $$t=1$$, $$r'=r_0$$ and $$t'=0$$) [so it seems to be correct];

2. I plotted the surface solution given in close form on the articles and it matches perfectly with the integral solution $$V(r,t) = \int_{r_{\min}}^{r_{\max}}e^{irr'}P(r,t\mid r',0)dr'$$ where $$r_{\min}$$ and $$r_{\max}$$ are chosen to be and interval around the expect value of $$r_t$$ of radius five times the variance of $$r_t$$ [so it seems the multiplying factor is $$e^{irr'}$$ but I still don't know why...].

NOTE_1: the reason why I tryed the term $$e^{irr'}$$ is that it is the new initial condition you get for $$P(r,t\mid r',t')$$ if you use the Fourier transform on equation (2), in the place of $$\delta(r-r')$$ (for reference check page 34 of this).

NOTE_2: I also tryed the term \$e^{1r'(t-0)} and it seems to work too... now I'm getting a bad feeling, maybe I've messed up with the coding part?