I am solving a coupled set of stochastic differential equations,

dx[t]=-kappa x[t] + g dW_t

dy[t]=-gamma x[t] + x[t] dW_t

where dW_t is a Wiener increment. Notice how they both share the same Wiener increment and y[t] is coupled to x[t].

To solve this numerically, we can easily implement the Euler-Maruyama method,

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import time
from numba import jit, njit
import os

#Parameters
kappa=1.
gamma=1.
g=1.

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

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

Mmax=15000 #number of samples

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

@njit(fastmath=True)
def SDE(Nmax,Mmax,nth):
    x_list =np.zeros((Nmax+1,Mmax),dtype=np.float64)
    y_list =np.zeros((Nmax+1,Mmax),dtype=np.float64)
    for m in range(Mmax):
        x=1
        y=0
        for n in range(Nmax+1):   
         
            y+=(-gamma*y)*dt+x*w[n,m] 
            x+=(-kappa*x)*dt+g*w[n,m]              
   
             
            x_list[n,m]=x
            y_list[n,m]=y
    return (x_list,y_list) 

(x,y)=SDE(Nmax,Mmax,nth)

x2=np.array([np.mean(x[n,:]*x[n,:]) for n in range(len(t_list))]).real
y2=np.array([np.mean(y[n,:]*y[n,:]) for n in range(len(t_list))]).real

And we are interested in the expected value of Y^2. The code works perfectly and I can benchmark analytically. However, this only works IF the order of equations is as I have written.

If I had the opposing order:

x+=(-gamma*x)*dt+g*w[n,m] 
y+=(-kappa*y)*dt+x*w[n,m]            
                           

the result is completely different.

Example: For the given parameters I give, the correct answer is around <Y^2>=0.25, while with the incorrect ordering, I get <Y^2>=1.20.

While in this case, it was somewhat inoffensive because I had the analytical result to benchmark, for a more complicated system where analytics is not possible, I believe this is a serious issue, as I would trust the given result, when it seems it was not correct.

Are there any intuitions as to why this happens?

PS: If I write my equations in matrix form, instead of ‘line by line’, I get the correct equation.