Python - Euler-Maruyama method for stochastic differential equation

Question

I am trying to numerically calculate the expectation value of some stochastic process.

Given the following stochastic differential equation, also known as Ornstein-Uhlenbeck process:

da[t]=-k a[t]+sqrt{k n} w[t]

This can be solved numerically to get a[t] using:

kappa=1
nth=1
gamma=1
Mmax=1000 #number of samples

Tmax=5. ##max value for time
dt=0.5
Nmax=int(Tmax/dt) ##number of steps

t_list=np.arange(0,Tmax+dt/2,dt)

w = np.random.randn(Nmax+1,Mmax)

a_list =np.zeros((Nmax+1,Mmax),dtype=np.complex64)

for m in range(Mmax):
    a=0+1j*0
    for n in range(Nmax+1):
        a=a+(-kappa*a)*dt+np.sqrt(dt)*np.sqrt(kappa*nth)*w[n,m]
        a_list[n,m]=a

a_mean=np.array([np.mean(a_list[n,:]) for n in range(len(t_list))])

To numerically integrate this, I use the Euler-Maruyama method. This means that the deterministic part of the ODE is solved by the usual Euler method, while the stochastic part goes with sqrt[dt].

With this, I can numerically simulate it and compute the mean values of a and higher moments. I also have an analytical expression for it and it agrees.

I would like now to calculate the expectation value of the variable with its time derivative, i.e. a[t]*da[t]. I can it do analytically and I get it should be sqrt{k n}/2. However, I do not get the same results when I do it numerically. Following the code from before, I extend it to:

for m in range(Mmax):
    a=0+1j*0
    for n in range(Nmax+1):
        a=a+(-kappa*a)*dt+np.sqrt(dt)*np.sqrt(dt*kappa*nth)*w[n,m]
        
        da=(-kappa*a)+np.sqrt(kappa*nth)*w[n,m] ### I define the da[t] here

        ex=a*da

   
        ex_list[n,m]=ex
ex_mean=np.array([np.mean(ex_list[n,:]) for n in range(len(t_list))])

The problem is that what I am trying to calculate is somewhat odd: the expectation value of a variable and its time derivative. And while I can do it analytically, numerical methods seem to fail. Is there any issue with the Euler-Maruyama integration scheme that doesn't work here?

Note: As per the analytical result, here is the idea. We need the expectation value of <a[t] da[t]>. By inserting the expression of <a[t] da[t]>=-k<a[t]a[t]>+sqrt[k n] <a[t]w[t]>. By formal integration of the equation I define above, a[t]=sqrt[k n]int_0^t exp[-k (t-t')] w[t]dt'. (This is a well-known result for the given process). (https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process)

The first term is <a[t]a[t]>=n/2.

The second term needs more work <a[t]w[t]>=sqrt[k n]\int_0^t dt' Exp[-k (t-t')]<w[t]w[t']>

Now we use that the white noise is delta correlated <w[t]w[t']>=d(t-t') and the integral above simplifies to k nth. The final expression is then <a[t]da[t]>=- k n/2 + k n= k n/2.

For the simple case of k=n=1, I should get 1/2. However, the numerical code I set gives -0.4 (more or less, due to averaging.)