Given the CIR process: $$\ dX_t = (a − bX_t ) dt +
\sigma dW_t$$
The Milstein scheme is:
$$\ X_{i+1} = X_i + \mu \Delta+\sigma\epsilon_i +\frac{1}{2}\sigma\sigma_x(\epsilon_i^2-\Delta)$$
and can be according to Elerian (1998)
written as
$$\ X_{i+1} = A_i(\epsilon_{i+1} + \sqrt{C_i}^2) + B_i $$
where $A_i$, $B_i$, $C_i$ are functions of of $X_i$ and I want to give $A_i$, $B_i$, $C_i$ for the CIR.
When defining $A_i = \frac{\sigma \sigma_x \Delta}{2}, B_i = -\frac{\sigma }{2\sigma_x}+X_i + \mu\Delta - A_i, C_i = \frac{1 }{\sigma_x^2\Delta }$
So what I did to prove this is I plugged in $A_i, B_i$ and $C_i $ to get back to the Milstein scheme:
$X_{i+1} = A_i\epsilon_i^2 + A_iC_i + 2\sqrt{C_i}A_i\epsilon_i +B_i$
$X_{i+1}= \frac{\sigma \sigma_x \epsilon_i^2\Delta}{2}+\frac{\sigma \sigma_x \Delta}{2\sigma_x^2\Delta}+\frac{2\sigma \sigma_x \epsilon_i\Delta}{2\sigma_x\sqrt{\Delta}}-\frac{\sigma }{2\sigma_x}+X_i + \mu\Delta - \frac{\sigma \sigma_x\Delta}{2}$ refining:
$X_{i+1} = \frac{\sigma \sigma_x \epsilon_i^2\Delta}{2}+\frac{\sigma}{2\sigma_x}+\sigma\epsilon_i\sqrt{\Delta}-\frac{\sigma }{2\sigma_x}+X_i + \mu\Delta - \frac{\sigma \sigma_x \Delta}{2}$ into:
$X_{i+1} = \frac{\sigma \sigma_x \epsilon_i^2\Delta}{2}+\sigma\epsilon_i\sqrt{\Delta}+X_i + \mu\Delta - \frac{\sigma \sigma_x \Delta}{2}$ but now I dont know how to further refine to get back to the Milstein scheme