I am trying to figure out the moment generating function (MGF) of $R_t$ from the CIR model through deriving the PDE, after applying Ito formula and taking expectation. However, it seems to be inconsistent with the MGF I found from other sources, as the MGF I found does not satisfy the PDE I obtain. I am not sure what is wrong with my method, as the solution I copied from different sources seems correct. Could someone help me to investigate the issue?
(The typesetting below is by courtesy of @TheBridge)
CIR Process :
$$dR_t=(a-bR_t)dt+ \sigma \sqrt {R_t}dW_t$$
It is known $R_t \sim \chi^2$ distribution with MGF $$m_t=\left(1-\frac{\sigma^2}{2b}(1-e^{-bt})\theta\right)^{-\frac{2a}{\sigma^2}}\exp\left(\frac{e^{-bt}\theta}{1-\frac{\sigma^2}{2b}(1-e^{-bt})\theta}R_0\right)$$ which is found in some literature.
I intend to derive this result. As $m(t,\theta)=\mathbb E (e^{\theta R_t})$, after applying Ito's formula, we have
$$de^{\theta R_t}=\theta e^{\theta R_t}dR_t+\frac{1}{2}\theta^2 e^{\theta R_t}\sigma^2 R_tdt=(a+\frac{1}{2}\theta\sigma^2 R_t-bR_t)\theta e^{\theta R_t}dt+(\cdots)dW_t$$
Taking expectation will produce $$m(t,\theta)= e^{\theta R_0}+\theta \left[ \int_0^t a\mathbb E[e^{\theta R_s}]ds+ \int_0^t(\frac{1}{2}\theta\sigma^2-b)\mathbb E[R_s e^{\theta R_s}]ds \right]$$
With $\frac{\partial m}{\partial \theta}=\mathbb E[R_s e^{\theta R_s}]$ then taking derivative wrt $t$ results in the PDE as follows:
$$\frac{\partial m}{\partial t}=\theta \left[a m(t,\theta)+(\frac{1}{2}\theta\sigma^2-b)\mathbb E[R_s e^{\theta.R_s}]\right]=\theta \left[a m(t,\theta)+(\frac{1}{2}\theta\sigma^2-b)\frac{\partial m}{\partial \theta}\right]$$ given initial condition $m(0,\theta)= e^{\theta R_0}$.
However, after using Matlab symbolic math for computation, the MGF I found (at the beginning of this post) does not satisfy my derived PDE. Could someone help to examine what is the potential issue in my derivation?