How to compute Kramers-Moyal coefficients when the time step is not to small in Python?

Question

I am coding a simple one-dimensional Ornstein-Uhlenbeck simulation in python to explain some experimental data from a video source, i.e., tracking a brownian particle in a sequence of frames. The sampling frequency of the camera is fps=30 frames per seconds, which means that my time resolution is dt=0.03 seconds, approximately. On the other hand, the smallest time scale I am able to measure is about tau=0.5 seconds, one order of magnitude larger than dt (a reasonable difference to be able to compute observables such as autocorrelation functions of the process).

Now, for the simulated data I am computing the first two Kramers-Moyal coefficients using kramersmoyal. The problem arises ONLY with the second coefficient when dt is not at least two order of magnitudes smaller than tau, in which case the estimation for the noise intensity becomes convex instead of a flat function.

So the questions at hand are:

• is there a simple way to code some corrections to address this issue?
• is there any other package that deals with this problem?
• is this an implementation problem or just a consequence that these coefficients are defined in the limit of infinitely small time steps (i.e., dt->0)?

In order to be more concrete, the code bellow shows the simulation for two different time steps, dt, that differs in two and one order of magnitude from tau, respectively. Then the second coefficient is computed for both cases and the noise intensity of the processes are estimated. The noise intensity looks flat for the first case but not for the second.

# === Ornstein-Uhlenbeck simulation ===

import numpy as np
import kramersmoyal as km
import matplotlib.pyplot as plt

np.random.seed(0)

n = 100_000       # total number of time steps
dt = [.003, .03]  # time steps two/one order of magnitud smaller than tau
tau = .5          # relaxation or characteristic time
sigma = 5         # noise intensity

scale_dW = sigma * np.sqrt(dt)    # scale parameter of the Wiener process
scale_x = sigma / np.sqrt(2/tau)  # scale parameter of the process under study

dW = np.random.normal(scale=scale_dW, size=(n, 2))  # Wiener process array
x = np.empty((n, 2))                                # x array
x[0] = np.random.normal(scale=scale_x, size=2)      # initial condition

# Euler-Maruyama method to solve the SDE
for i in range(n - 1):
    x[i + 1] = x[i] - 1/tau * x[i] * dt + dW[i]

# === Second Kramers-Moyal coefficient ===

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(9,3.5))
bins = np.array([100])         # number of bins
powers = np.array([[1], [2]])  # first and second power

for i in range(len(dt)):
    # computation of the coefficients using `kramersmoyal`
    kmc, edges = km.km(x[:,i], bins, powers)

    # sigma parameter estimation
    sigma_est = np.sqrt(2 * kmc[1] / dt[i])

    # plotting
    axs[i].plot(edges[0], sigma_est, '.-', label=f'dt = {dt[i]}')
    axs[i].legend()
    axs[i].set_title(['FLAT', 'NOT FLAT'][i])
    axs[i].set_xlabel('x')
    axs[i].set_ylabel('estimated sigma')
plt.show()

Answer 1

Interesting development, using the ISO code NZ works without error - despite documentation saying it needs to be the name of the timezone.

@convertTimeZone(utcNow(), 'UTC', 'NZ', 'yyyy-MM-ddTHH:mm:ss')

Answer 2

If you cannot using convert, you can using addHours .. in my country Indonesia is gmt +7 so Im using addHours(utcNow(), 7)

If you in new zealand you can change from 7 to 13