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

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.


[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]

  • $\begingroup$ 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}$. $\endgroup$
    – G. Gare
    Jul 26, 2022 at 11:44
  • $\begingroup$ @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. $\endgroup$
    – H K Y
    Aug 4, 2022 at 12:48
  • $\begingroup$ 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. $\endgroup$
    – G. Gare
    Aug 4, 2022 at 13:42
  • $\begingroup$ @G.Gare thank you very much for your insight it was very clear. I shall implement that and come back to you asap. $\endgroup$
    – H K Y
    Aug 5, 2022 at 21:00


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.