Stationary distribution of Cox-Ingersoll-Ross process I am uncertain how to go about the following problem from the lecture notes on a course in SDE's. We are given the following SDE.
$dX_t=\lambda\left(\xi-X_t\right) dt+\gamma\sqrt{|X_t|}dB_t$
Where $\lambda,\xi,\gamma$ are positive constants. Show that $X_t$ has Gamma distribution, with a rate parameter $\omega=2\lambda/\gamma^2$ and shape parameter $\nu=2\lambda\xi/\gamma^2$. Finally, what is the mean and variance in stationarity? 
I know I have to find the stationary distribution by isolating $\phi$ in the following.
$-\nabla\cdot\left( u\phi-D\nabla \phi\right)=0$
Inserting $D=\frac{1}{2}\gamma^2|X_t|, $ $u=f-\nabla D=\lambda\left(\xi-X_t\right)-\frac{1}{2}\gamma^2\frac{X_t}{|X_t|}$, I get an ugly expression, involving second derivative of absolute value of the stochastic process. Is there a mistake so far? (I take the gradient and divergence w.r.t $x=X_t$)
 A: From the Kolmogorov forward equation,
$$
\frac{d^2}{dy^2}\left(\frac{1}{2}\gamma^2y\phi(y)\right)=\frac{d}{dy}\left(\lambda(\xi-y)\phi(y)\right).
\tag{1}\label{eq:asdf}
$$
The left-hand side is equal to
$$
\frac{d}{dy}\left(\frac{1}{2}\gamma^2\phi(y)+\frac{1}{2}\gamma^2y\phi'(y)\right).
$$
So integrating \eqref{eq:asdf} from $0$ to $y$ we get
$$
\frac{1}{2}\gamma^2\phi(y)+\frac{1}{2}\gamma^2y\phi'(y)=\lambda(\xi-y)\phi(y).
$$
Here, I've used the fact that $\phi(0)=0$ (if $X$ hits zero, it'll instantly get reflected). [$\phi(0)$ is not necessarily zero. See below.] Rearranging the last equation,
$$
\frac{\phi'(y)}{\phi(y)}
=
\frac{\lambda(\xi-y)-\frac{1}{2}\gamma^2}{\frac{1}{2}\gamma^2y}
=
\left(\frac{2\lambda\xi}{\gamma^2}-1\right)\frac{1}{y}-\frac{2\lambda}{\gamma^2}.
$$
Integrating this from $\xi$ to $y$,
$$
\log\frac{\phi(y)}{\phi(\xi)}
= \left(\frac{2\lambda\xi}{\gamma^2}-1\right)\log\frac{y}{\xi} - \frac{2\lambda}{\gamma^2}(y-\xi).
$$
So, finally,
$$
\phi(y)\propto y^{\frac{2\lambda\xi}{\gamma^2}-1}e^{- \frac{2\lambda}{\gamma^2}y},
$$
which is the pdf of a Gamma distribution.
The mean and variance are accordingly given by
$$
\frac{2\lambda\xi}{\gamma^2}\frac{\gamma^2}{2\lambda}=\xi
\quad\text{and}\quad
\frac{2\lambda\xi}{\gamma^2}\left(\frac{\gamma^2}{2\lambda}\right)^2
=\frac{\gamma^2\xi}{2\lambda}.
$$

Addendum
If $\frac{2\lambda\xi}{\gamma^2}\le 1$, then $\phi(0)\ne 0$; and if $\frac{2\lambda\xi}{\gamma^2}< 1$, then $\phi'(0)=-\infty$, contrary to what I assumed above. In any case, these assumptions are not necessary. Here's the correction.
We still have
$$
\frac{1}{2}\gamma^2\phi(y)+\frac{1}{2}\gamma^2y\phi'(y)
-\lambda(\xi-y)\phi(y)=\text{constant}.
$$
Since $\phi$ is a pdf defined on $[0,\infty)$, $\phi(\infty)=0$ and $\phi'(\infty)=0$. Also, if $X_t$ is integrable with respect to the stationary distribution (which I will assume), then $y\phi(y)\rightarrow 0$ and $y\phi'(y)\rightarrow 0$ as $y\rightarrow \infty$. It follows that the constant in the last equation is zero, which brings us back to the derivation above.
