# CEV model effective simulation

I want to simulate the following CEV process : $$dM_t = \sigma_t M_t^{\eta} dW_t$$

Using Euler discretization to $$M_t$$, if at a given time $$t$$, $$M_t$$ takes a negative value then $$M_{t+1} = M_t + \sigma_t M_t^{\eta} \sqrt{dt} \epsilon$$, where $$\epsilon$$ follows a standard Gaussian distribution, gives 'NaN's when $$\eta < 1$$. A first solution might be to take the positive value such as : $$M_{t+1} = \max(0,M_t + \sigma_t M_t^{\eta} \sqrt{dt} \epsilon)$$. The problem with this method is that I get to set a floor (here at 0) and that affects the impact results. I still get to preserve my martingale property with this method : $$E(M_t) = M_0$$, but it's not quite good because the Euler tends to give many negative values which are then floored to zero and there is no rebound, the zero becomes "absorbant".

Therefore, I went for a variable change in order to have values justifiable without setting a floor : $$dM_t = \sigma_t M_t^{\eta} dW_t$$ Let the variable $$Q_t$$ as : $$Q_t = \frac{1}{1 - \eta} M_t^{1 - \eta}$$ Applying Itô formula to $$Q_t$$ gives : $$dQ_t = \sigma_t dW_t - \frac{\eta}{2(1 - \eta)} \sigma_t^2 Q_t^{-1}dt$$

I can therefore simulate $$Q_t$$ without having to deal with the negative part and deduce $$M_t$$ as : $$M_t = [(1 - \eta) Q_t]^{\frac{1}{1 - \eta}}.$$

My problem with this approach is that I no longer verify the martingale property $$E(M_t) = M_0$$. Does anybody know why ?

Here is my Python code for this simulation :

import numpy as np
np.random.seed(1)

nb_sim = int(1e4)
nb_steps = int(1e2)

z = np.random.normal(size = [nb_sim, nb_steps])

eta = 0.5
vol = 0.2 * np.ones(nb_steps)

m = np.ones(z.shape)
q = np.power(m, 1 - eta)/ (1 - eta)

# $$dQ_t = \sigma_t dW_t - \frac{\eta}{2(1 - \eta)} \sigma_t^2 Q_t^{-1} dt$$

dt = 1

for t in range(nb_steps-1):

q[:,t+1] = q[:,t] + vol[t] * np.sqrt(dt) * z[:,t] - dt * eta * np.power(vol[t], 2) / (2 * (1 - eta) * q[:,t])

# Deduce $$M_t$$ from $$Q_t$$ :
# $$M_t = [(1 - \eta) Q_t]^{\frac{1}{1 - \eta}}$$

m = np.power((1 - eta) * q, 1 / (1 - eta))

# Martingale test : $$E(M_t) = M_0$$
E = m.mean(axis=0)

# The expectation E should be values around M_0 = 1 in this specific example, for all t.
print(E)


print(np.round(E))

[1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 4.000e+00 4.000e+00 4.000e+00 8.000e+00 8.000e+00 8.000e+00 8.000e+00 8.000e+00 8.000e+00 1.600e+01 1.600e+01 2.200e+01 2.200e+01 2.400e+01 2.400e+01 2.400e+01 2.700e+01 2.800e+01 2.800e+01 2.800e+01 2.800e+01 2.800e+01 2.900e+01 2.900e+01 3.000e+01 3.400e+01 3.600e+01 3.600e+01 5.290e+02 5.320e+02 5.320e+02 5.320e+02 5.320e+02 5.320e+02 5.330e+02 5.360e+02 5.360e+02 5.370e+02 5.370e+02 5.370e+02 5.590e+02 5.590e+02 5.600e+02 5.600e+02 5.660e+02 5.660e+02 5.660e+02 5.660e+02 5.670e+02 5.670e+02 5.670e+02 5.720e+02 5.800e+02 5.800e+02 5.810e+02 5.850e+02 5.850e+02 5.890e+02 5.890e+02 5.890e+02 5.890e+02 5.910e+02 5.910e+02 5.910e+02 5.910e+02 5.930e+02 5.950e+02 5.950e+02 5.970e+02 5.970e+02 5.970e+02 5.970e+02 5.970e+02 5.970e+02 6.000e+02 6.010e+02 6.010e+02 6.020e+02 6.030e+02 1.889e+03 1.894e+03 1.894e+03 1.895e+03 1.896e+03]

• The problem of your model (and your initial model, too) is that when $Q$ approaches zero then at the next step it jumps to high values due to the term $Q^{-1}$. Commented Jul 26, 2022 at 11:44
• @G.Gare thank you for your answer. I still can't figure out how I can improve my simulation process and I am open to any suggestion, thanks. Commented Aug 4, 2022 at 12:48
• The problem is that your diffusion function $\sigma M_t^\eta$ is not globally Lipschitz for $\eta < 1$, as in zero you have a vertical asymptote (think of the graph of $\sqrt{x}$). I think you may look into some adaptive Euler Maruyama code, for example reducing the time step whenever $M_t$ approaches zero could be resolutive. Commented Aug 4, 2022 at 13:42
• @G.Gare thank you very much for your insight it was very clear. I shall implement that and come back to you asap. Commented Aug 5, 2022 at 21:00