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′<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:
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];
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?
