Mistake in `jd_process.py`

[JChonpca]

Thank you for your great effort in the jumpdiff, it is the first package in both simulation and estimation for the research in SDEwJ with Python language. But I find some mistakes in your simulation parts.

In the code file jd_process.py, you directly treat the jump as the dJ[i] * np.random.normal(loc = 0, scale = np.sqrt(xi)) for both Euler and Milstein schemes, which is a wrong simulation method. The correct method is treat the diffusion as the sum of d[j] iid samples that from normal distribution. The data generated with wrong data will leads to the estimation bias of your program. The example you used in the introduction in github and publication has very small lambda and the dJ[i] is mostly equal to 1, so the data generated with wrong simulation method is almost the same with the right one. If you change the lambda to 50, the estimation will crash with big bias of about 30, which is showed the code file attached. After I change the simulation data using the correct simulation method, the estimation task is worked for both small and big lambda.

I think your estimation method is good without any problem since I have tested other sample with the right simulation code, but you should correct the mistakes in jd_process.py.

Wrong Code in jd_process.py:

# Integration, either Euler
if solver == 'Euler':
    for i in range(1, length):
        X[i] = X[i-1] + a(X[i-1]) * delta_t + b(X[i-1]) * dw[i]
        if dJ[i] > 0.:
            X[i] += dJ[i] * np.random.normal(loc = 0, scale = np.sqrt(xi))

if solver == 'Milstein':
    # Generate corrective terms of the Milstein integration method
    dw_2 = (dw**2 - delta_t) * 0.5

    for i in range(1, length):
        X[i] = X[i-1] + a(X[i-1]) * delta_t + b(X[i-1]) * dw[i] \
                + b(X[i-1]) * b_prime(X[i-1]) * dw_2[i]
        if dJ[i] > 0.:
            X[i] += dJ[i] * np.random.normal(loc = 0, scale = np.sqrt(xi))

Correct Code for jd_process.py:

# Integration, either Euler
if solver == 'Euler':
    for i in range(1, length):
        X[i] = X[i-1] + a(X[i-1]) * delta_t + b(X[i-1]) * dw[i]
        if dJ[i] > 0.:
            for j in range(int(dJ[i])):
                X[i] += np.random.normal(loc = 0, scale = np.sqrt(xi))

if solver == 'Milstein':
    # Generate corrective terms of the Milstein integration method
    dw_2 = (dw**2 - delta_t) * 0.5

    for i in range(1, length):
        X[i] = X[i-1] + a(X[i-1]) * delta_t + b(X[i-1]) * dw[i] \
                + b(X[i-1]) * b_prime(X[i-1]) * dw_2[i]
        if dJ[i] > 0.:
            for j in range(int(dJ[i])):
                X[i] += np.random.normal(loc = 0, scale = np.sqrt(xi))

Simulation and Estimation wth big lambda:

# -*- coding: utf-8 -*-
"""
Created on Fri Feb 10 15:19:29 2023

@author: JChonpca_Huang
"""
import jumpdiff as jd
import matplotlib.pyplot as plt

# integration time and time sampling
t_final = 1000
delta_t = 0.001

# A drift function
def a(x):
    return -0.5*x

# and a (constant) diffusion term
def b(x):
    return 0.75

# Now define a jump amplitude and rate

xi = 2.5

lamb = 50

# and simply call the integration function
X = jd.jd_process(t_final, delta_t, a=a, b=b, xi=xi, lamb=lamb)

plt.plot(X)
plt.show()

edges, moments = jd.moments(timeseries = X)

# we want the first power, so we need 'moments[1,...]'
plt.plot(edges, moments[1,...])
plt.show()

plt.plot(edges, moments[2,...])
plt.show()

# first estimate the jump amplitude
xi_est = jd.jump_amplitude(moments = moments)

print(xi_est)

# and now estimated the jump rate
lamb_est = jd.jump_rate(moments = moments)

print(lamb_est/delta_t)

[LRydin]

Dear @JChonpca you are absolutely correct! Thank you very much. Created PR #6, altering line 109 and line 119.

[JChonpca]

Thank you for your immediate reply and bug fixing. Jumpdiff is a very nice python package and greatly helps with my research.

[LRydin]

Thanks again @JChonpca, I have merged PR #6. I'll push a new version, v0.4.3 with this correction.