Negative values in CIR model

Question

I'm having difficulty understanding the well known property of the CIR model that it can't go below zero. Wikipedia says that this is because the random shock on the rate will grow very small as r moves closer to zero, but won't the drift term also do so? Especially if the volatility term is high, isn't it possible for the random shock to be of a greater negative value than the drift is positive even as r goes to zero?

I'm trying to implement the method in matlab currently but it happens to me that r becomes negative if I increase the volatility. Could it be a problem with the discretization as well perhaps? The code snippet is below if it's of interest.

theta=0.5; %Long run mean
sigma=14; %Volatility of drawdowns
k = 7.326; %Mean reversion constant
n = 100; %number of time steps, t.
dt = T/n; %time step
M=10^3; %Number of realizations
d0 = theta;
d=ones(M,1).*d0;%Starting value for d
for i = (j-1:n)
    dW = sqrt(dt)*randn(M,1); % Wiener increments
    d(:,i+1) = d(:,i) + k.*(theta-d(:,i)).*dt + sigma.*sqrt(d(:,i)).*dW; %drawdown rate
end

```

Answer 1

Analytically the Feller condition ($ 2 \kappa \theta > \sigma ^2)$ guarantees that the process doesn't become negative but this is not enough when you are simulating. Even if you choose parameters that satisfy the Feller condition, you still may have the problem of getting negative values inside the square root giving bad results. This is a consequence of the discretization. The easiest way to solve the problem is to substitute in your code sqrt(d(:,i)) by sqrt(max(d(:,i),0)). There are other more sophisticated ways to deal with this, you can see here https://www.deriscope.com/docs/Andersen_Jaeckel_Kahl_2010.pdf

Answer 2

If you don't want that your CIR process goes below zero, this condition should be satisfied: $$2k\theta>\sigma^2$$ Thus, you can't choose a very big value for volatility. It is limited by this condition.