The CIR model for spot rate $r_t$ is:
$$dr_t=(\eta-\gamma r_t)dt+\sqrt{\alpha r_t} dW_t$$
where $\eta, \gamma, \alpha$ are constants.
How to express this SDE in discrete form using Milstein scheme?
The one I derived is:
$$r_{t+1}=r_t+(\eta-\gamma r_t)\delta t+\sqrt{\alpha r_t}\cdot\sqrt{\delta t}\phi +\frac{1}{2}\sqrt{\alpha r_t}\cdot\left(\frac{1}{2}\frac{\alpha}{\alpha r_t}\right)[\delta t(\phi^2-1)]$$
where $\phi$ is normal RV.
Can anyone help me to identify my error? Or is it correct?
The Milstein scheme for the following CIR model
$$dr_t=(\eta-\gamma r_t)dt+\sqrt{\alpha r_t} dW_t$$
should be
$$r_{t+1}=r_t+(\eta-\gamma r_t)\delta t+\sqrt{\alpha r_t}\cdot\sqrt{\delta t}\phi +\frac{1}{2}\sqrt{\alpha r_t}\cdot\left(\frac{1}{2}\sqrt{\frac{\alpha}{r_t}}\right)[\delta t(\phi^2-1)]$$
$$r_{t+1}=r_t+(\eta-\gamma r_t)\delta t+\sqrt{\alpha r_t}\cdot\sqrt{\delta t}\phi +\frac{1}{4}\alpha(\phi^2-1)\delta t$$ where $\phi$ is normal RV.
I think that you wrongly derived $\frac{\partial\left(\sqrt{\alpha r_t}\right)}{\partial r_t}$ in the last term.
${{r}_{t+\Delta t}}={{r}_{t}}+(\eta-\gamma r_t)\delta t+\sigma \sqrt{{{r}_{t}}}\sqrt{\delta t}\,{{Z}}+\frac{1}{4}{{\sigma }^{2}}\delta t({{Z}}^{2}-1)$
Euler is just bad. Milstein in my experience is not much better. A better scheme in general is Predictor Corrector.
Another approach (in Glasserman) is to simulate the SDE using the non-central Chi^2 distribution which resolves Euler problem.
