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]
