Monte-Carlo for SDE with square root diffusion term I've recently got a question from a Master student about a numerical simulation/integration of the SDE of the following shape
$$
  \mathrm dX_i(t) = \left(\sum_{j=1}^M \nu_{ji} a_jX_j(t)\right)\mathrm dt + \sum_{j=1}^M\left(\nu_{ji}\sqrt{a_j X_j(t)}\right)\mathrm dW_j(t)
$$
where $\nu_{ji}$ are entries of the integer-valued matrix (stoichiometry matrix), and $a_j$ are positive reals (rates of chemical reactions). The student tried the Euler-Maruyama scheme
$$
  X_i(t+\Delta t) = X_i(t) + \left(\sum_{j=1}^M \nu_{ji} a_jX_j(t)\right)\Delta t + \sum_{j=1}^M\left(\nu_{ji}\sqrt{a_j X_j(t)\Delta t}\right)\xi_j
$$
where $\xi_j$ are iid standard normal random variables. However, for $\Delta t$ small enough the scheme does not converge and blows up drastically. This is of course expected: for small $\Delta t$ the term in front of $\xi_j$ significantly dominates the deterministic drift term. 
We also looked into higher-order schemes for the integration, such as Runge-Kutta method, but they all blow up due to similar reasons. I would be happy to hear of any hint how to approach this problem.
 A: You can still use the Euler-Maruyama scheme, but it will be slower than usual.
There is nothing stopping you from using the Euler-Maruyama scheme for the SDE, but there are two factors to be taken into account. The first is that the SDE process is very likely positive. This is very similar to a CIR process, and so you can likely perform a bit of analysis on the SDE. However, to ensure your approximation is always positive (and you never take the square root of a negative number) you will likely need a truncation scheme. On a similar note it will require a bit of care to see if the boundary is absorbing or strongly reflective. For a CIR process you can assess this using the Feller condition, and your process will also require similar care. For some related material on similar ideas see A comparison of biased simulation schemes for
stochastic volatility models by 
Lord, Koekkoek, and van Dijk.
Moving away from some of the technical details, it's also worth mentioning that the Euler-Maruyama scheme will likely have a much slower rate of convergence. Had the volatility term been Lipschitz, standard theory would be able to give to weak and strong convergence of order 1 and $\tfrac{1}{2}$ respectively. However, the square root process is not so nice, and while the scheme will still converge, it will be slower. Gyongy and Rasonyi provide a rate of convergence for similar processes in A note on Euler approximations for SDEs with Holder continuous diffusion coefficients, where the rate of MC error decay can drop from $O(n^{1/2})$ down to $O(\tfrac{1}{\log(n)})$.
Are you sure you don't have a bug?

However, for $\Delta t$ small enough the scheme does not converge and blows up drastically. This is of course expected: for small $\Delta t$ the term in front of $\xi_j$ significantly dominates the deterministic drift term.

The dominance of the stochastic term is to be expected. However, that your scheme is "blowing up" is not. In honesty it sounds like a bug in your simulation. However, unbounded paths (i.e. very big numbers in your simulation) might be perfectly feasible. I think a question of note here is "how are you quantifying if the  error?". Ultimately, each approximate paths is a random variable. 
Picard iterations?
If you really want a sample path where you want increasing levels of confidence about its solution to the SDE (really best interpreted using integrals) then you could use Picard iterations. Just another idea. 
