examples/kmc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A one dimensional Ornstein−Uhlenbeck process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple Ornstein--Uhlenbeck process takes the form\n",
    "$$\\mathrm{d}y(t) = -\\theta y(t)\\mathrm{d}t + \\sigma \\mathrm{d}W(t), $$\n",
    "with $\\theta$ the drift term and $\\sigma$ diffusion term. The term $\\mathrm{d}W(t)$ comprises the Brownian motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrating with an Euler−Maruyama scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the equivalent of the simplest Euler integrator for stochastic processes. You can find more specifics on the Euler−Maruyama onwikipedia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The drift and diffusion\n",
    "theta = 0.3\n",
    "sigma = 0.1\n",
    "\n",
    "# The total integration time\n",
    "t_start = 0\n",
    "t_final = 500\n",
    "\n",
    "# The desired timestep of integration\n",
    "delta_t = 0.001\n",
    "\n",
    "# time array of the process\n",
    "time = np.linspace(0, t_final, t_final * int(1 / delta_t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the integration initialise an array $y$ and a normally distributed noise $\\mathrm{d}W(t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the array y\n",
    "y = np.zeros([time.size])\n",
    "\n",
    "# Generate a Wiener process with a scale of np.sqrt(delta_t)\n",
    "dW = np.random.normal(loc = 0, scale = np.sqrt(delta_t), size = [time.size,1])\n",
    "# you can equivalently use np.sqrt(delta_t) * np.random.normal(loc=0, scale=1,size=[time.size,1])\n",
    "\n",
    "# Give some small random initial conditions\n",
    "y[0]=np.random.normal(size = 1) / 10\n",
    "\n",
    "# Integrate the process\n",
    "for i in range(1, time.size):\n",
    "    y[i] = y[i-1] - theta * y[i-1] * delta_t + sigma * dW[i] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising the process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is the stochastic trajectory over time\n",
    "plt.plot(time, y , label = r'Trajectory of O–U process')\n",
    "\n",
    "plt.xlabel(r'time $t$')\n",
    "plt.ylabel(r'$y(t)$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Kramers−Moyal coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now implement the KM package to calculate the Kramers−Moyal coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import your (favourite) package to calculate Kramers−Moyal coefficients\n",
    "# You will need only km to perform the calculation\n",
    "from kramersmoyal import km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the size of your target space \n",
    "bins = np.array([5000])\n",
    "\n",
    "# Introduce the desired orders to calculate\n",
    "# Please keep the  [0]  term. It is the normalisation. \n",
    "powers = np.array([[0], [1], [2], [3], [4], [5], [6], [7], [8]])\n",
    "\n",
    "# Choose a desired bandwidth bw\n",
    "# this is similar to the inverse of number of bins\n",
    "bw = 0.1\n",
    "\n",
    "# Calculate the Kramers−Moyal coefficients\n",
    "kmc, edges = km(y, bw = bw, bins = bins, powers = powers)\n",
    "\n",
    "# The K−M coefficients are stacked along the fisrt dim of the\n",
    "# kmc array, so kmc[1,:] is the first K−M coefficient, kmc[2,:]\n",
    "# is the second, etc. etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising the Kramers−Moyal coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now visualise the obtained Kramers−Moyal coefficients\n",
    "Let's plot the first Kramers−Moyal coefficients, i.e., the drift coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets restrict the plot to where we have good statistics\n",
    "# but include the timestep normalisation 1 / delta_t\n",
    "plt.plot(edges[0][1500:-1500], kmc[1,1500:-1500], label = r'First K−M coefficient: drift')\n",
    "\n",
    "# And to guide the eye, here is the actual theoretical value,\n",
    "plt.plot(edges[0][1500:-1500], -edges[0][1500:-1500] * theta * delta_t, ':', label = r'Theory', color=r'black')\n",
    "\n",
    "# For a physics understanding of the process, recall to always divide \n",
    "# the Kramers−Moyal coefficients by / delta_t\n",
    "\n",
    "plt.xlabel(r'$y(t)$')\n",
    "plt.ylabel(r'K−M coefficient')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now the second Kramers−Moyal coefficients, i.e., the diffusion coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets restrict the plot to where we have good statistics\n",
    "plt.plot(edges[0][1500:-1500], kmc[2,1500:-1500], label = r'Second K−M coefficient: diffusion')\n",
    "\n",
    "# And to guide the eye, here is the actual theoretical value,\n",
    "# which is a constant of value sigma**2 * delta_t /2\n",
    "plt.plot(edges[0][1500:-1500], np.ones(edges[0][1500:-1500].size) * sigma**2 * delta_t / 2, ':', label = r'Theory', color=r'black')\n",
    "\n",
    "plt.xlabel(r'$y(t)$')\n",
    "plt.ylabel(r'K−M coefficient')\n",
    "plt.ylim([sigma**2 * delta_t / 4, sigma**2 * delta_t])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A 2-dimensional diffusion process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import matplotlib's 3d plotting library\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A two-dimensionan Ornstein−Uhlenbeck  diffusion process takes the form\n",
    "$$\n",
    "    \\begin{pmatrix}\n",
    "    \\mathrm{d}y_1(t) \\\\ \\mathrm{d}y_2(t)\n",
    "    \\end{pmatrix}= \n",
    "    \\begin{pmatrix}\n",
    "     N_1(y) \\\\ N_2(y)\n",
    "    \\end{pmatrix}\n",
    "    \\mathrm{d} t +\n",
    "    \\begin{pmatrix}\n",
    "    g_{1,1}(y) & g_{1,2}(y) \\\\\n",
    "    g_{2,1}(y) & g_{2,2}(y)\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "    \\mathrm{d}W_1 \\\\ \\mathrm{d}W_2\n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "with $N$ the drift vector and $g$ the diffusion matrix, which can be state dependent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define as the previous example a process identical to the Ornstein−Uhlenbeck process, with\n",
    "$$\n",
    "    N = \\begin{pmatrix}\n",
    "    - N_1 y_1 \\\\ - N_2 y_2\n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "and we take $N_1=2.0$ and $N_2=1.0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this particular case a more involved diffusion matrix $g$ will be used. Let the matrix $g$ be state-dependent, i.e., dependent of the actual values of $y_1$ and $y_2$ via\n",
    "$$\n",
    "    g = \\begin{pmatrix}\n",
    "    \\frac{g_{1,1} }{1+e^{y_1^2}} & g_{1,2} \\\\\n",
    "    g_{2,1} & \\frac{g_{2,2} }{1+e^{y_2^2}}\n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "and we will take $g_{1,1}=g_{2,2}=0.5$ $g_{1,2} = g_{2,1} = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation of the diffusion noise above in 2-dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the figure for 3d ploting\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "x_1 = np.linspace(-2, 2, 30)\n",
    "x_2 = np.linspace(-2, 2, 30)\n",
    "X_1, X_2 = np.meshgrid(x_1, x_2)\n",
    "Z = 1/(1+np.exp(X_1**2)) * 1/(1+np.exp(X_2**2))\n",
    "ax.contour3D(X_1, X_2, Z, 50, cmap='gray')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrating with an Euler−Maruyama scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define first the drift vector N with some exemplary values\n",
    "N = np.array([2.0, 1.0])\n",
    "\n",
    "# Define the diffusion matrix g with some exemplary values, but below\n",
    "# we add a functional form to it\n",
    "g = np.array([[0.5, 0.0],[0.0, 0.5]])\n",
    "\n",
    "# The total integration time\n",
    "t_start = 0\n",
    "t_final = 2000\n",
    "\n",
    "# notice in 2-dimensional we needs more points, i.e., longer integration \n",
    "\n",
    "# The desired timestep of integration\n",
    "delta_t = 0.001\n",
    "\n",
    "# time array of the process\n",
    "time = np.linspace(0, t_final, t_final * int(1 / delta_t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialise an (n,2)-array $y$, (n,2)-normally distributed noise $\\mathrm{d}W_{1,2}(t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the array y\n",
    "y = np.zeros([time.size, 2])\n",
    "\n",
    "# Generate two Wiener processes with a scale of np.sqrt(delta_t)\n",
    "dW = np.random.normal(loc = 0, scale = np.sqrt(delta_t), size = [time.size, 2])\n",
    "\n",
    "# Give some small random initial conditions\n",
    "y[0,:]=np.random.normal(size = 2) / 10\n",
    "\n",
    "# Integrate the process (takes about 20 secs)\n",
    "for i in range(1, time.size):\n",
    "    y[i,0] = y[i-1,0]  -  N[0] * y[i-1,0] * delta_t + g[0,0]/(1 + np.exp(y[i-1,0]**2)) * dW[i,0]  +  g[0,1] * dW[i,1]\n",
    "    y[i,1] = y[i-1,1]  -  N[1] * y[i-1,1] * delta_t + g[1,0] * dW[i,0]  +  g[1,1]/(1 + np.exp(y[i-1,1]**2)) * dW[i,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising the process in 2 dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the figure for 3d ploting\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "# Ploting the trajector for the first 10 time units, which are\n",
    "# 10000 = 10 / delta_t data points\n",
    "ax.plot(y[:10000,0], time[:10000], y[:10000,1], color='black')\n",
    "ax.view_init(30, 20)\n",
    "\n",
    "ax.set_xlabel(r'$y_1$'); ax.set_ylabel(r'time'); ax.set_zlabel(r'$y_2$');\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Kramers−Moyal coefficients in 2 dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now implement the KM package to calculate the Kramers−Moyal coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import your (favourite) package to calculate Kramers−Moyal coefficients\n",
    "# You will need only km to perform the calculation\n",
    "from kramersmoyal import km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the size of your target space in two dimensions \n",
    "bins = np.array([300, 300])\n",
    "# Be careful, this is a grid with 300*300 = 90000 points. If you want something\n",
    "# large, make sure not to blow up your memory\n",
    "\n",
    "# Introduce the desired orders to calculate, but in 2 dimensions\n",
    "# Please keep the [0,0] term. It is the normalisation. \n",
    "powers = np.array([[0,0], [1,0], [0,1], [1,1], [2,0], [0,2], [2,2]])\n",
    "# insert into kmc:   0      1      2      3      4      5      6\n",
    "\n",
    "# Notice that the first entry in [,] is for the first dimension, the \n",
    "# second for the second dimension...\n",
    "\n",
    "# Choose a desired bandwidth bw\n",
    "bw = 0.1\n",
    "\n",
    "# Calculate the Kramers−Moyal coefficients\n",
    "kmc, edges = km(y, bw = bw, bins = bins, powers = powers)\n",
    "\n",
    "# The K−M coefficients are stacked along the first dim of the\n",
    "# kmc array, so kmc[1,...] is the first K−M coefficient, kmc[2,...]\n",
    "# is the second. These will be 2-dimensional matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising the 2-dimensional Kramers−Moyal coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets visualise the Kramers−Moyal coefficients [1,0], [0,1], [2,0], and [0,2], i.e., the drifts and diffusion terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the figure for 3d ploting\n",
    "fig = plt.figure(figsize = (12,8))\n",
    "\n",
    "#lets fix a range where we have good statistics and generate a meshgrid\n",
    "X_1, X_2 = np.meshgrid(edges[0][80:-80],edges[1][80:-80])\n",
    "\n",
    "\n",
    "# the Kramers−Moyal coefficients [1,0]\n",
    "ax_00 = fig.add_subplot(2, 2, 1, projection='3d')\n",
    "\n",
    "# kmc[1,...] i.e, power [1,0] (transpose, python stores arrays transposed)\n",
    "ax_00.contour3D(X_1, X_2, kmc[1,80:-80,80:-80].T / delta_t, 50, cmap='Greens')\n",
    "# The actual theoretical value\n",
    "ax_00.contour3D(X_1, X_2, - X_1 * N[0], 50, cmap='gray')\n",
    "\n",
    "ax_00.set_title(r'K−M coefficient [1,0]');\n",
    "\n",
    "\n",
    "# the Kramers−Moyal coefficients [0,1]\n",
    "ax_01 = fig.add_subplot(2, 2, 2, projection='3d')\n",
    "\n",
    "# kmc[2,...] i.e, power [0,1] (transpose, python stores arrays transposed)\n",
    "ax_01.contour3D(X_1, X_2, kmc[2,80:-80,80:-80].T / delta_t, 50, cmap='Greens')\n",
    "# The actual theoretical value\n",
    "ax_01.contour3D(X_1, X_2, - X_2 * N[1], 50, cmap='gray')\n",
    "\n",
    "ax_01.set_title(r'K−M coefficient [0,1]');\n",
    "\n",
    "\n",
    "# the Kramers−Moyal coefficients [2,0]\n",
    "ax_10 = fig.add_subplot(2, 2, 3, projection='3d')\n",
    "\n",
    "# kmc[4,...] i.e, power [2,0] (transpose, python stores arrays transposed)\n",
    "ax_10.contour3D(X_1, X_2, kmc[4,80:-80,80:-80].T / delta_t, 50, cmap='Greens')\n",
    "# The actual theoretical value\n",
    "ax_10.contour3D(X_1, X_2,  ( g[0,0] / (1 + np.exp(X_1**2)) )**2 / 2 , 50, cmap='gray')\n",
    "\n",
    "ax_10.set_title(r'K−M coefficient [2,0]');\n",
    "\n",
    "\n",
    "# the Kramers−Moyal coefficients [0,2]\n",
    "ax_11 = fig.add_subplot(2, 2, 4, projection='3d')\n",
    "\n",
    "# kmc[5,...] i.e, power [0,2] (transpose, python stores arrays transposed)\n",
    "ax_11.contour3D(X_1, X_2, kmc[5,80:-80,80:-80].T / delta_t, 50, cmap='Greens')\n",
    "# The actual theoretical value\n",
    "ax_11.contour3D(X_1, X_2,  ( g[1,1] / (1 + np.exp(X_2**2)) )**2 / 2 , 50, cmap='gray')\n",
    "\n",
    "ax_11.set_title(r'K−M coefficient [0,2]');\n",
    "\n",
    "# Rotate views and add labels\n",
    "ax_00.view_init(30, 20); ax_01.view_init(30, 20); ax_10.view_init(30, 20); ax_11.view_init(30, 20)\n",
    "ax_00.set_xlabel(r'$y_1$'); ax_01.set_xlabel(r'$y_1$'); ax_10.set_xlabel(r'$y_1$'); ax_11.set_xlabel(r'$y_1$')\n",
    "ax_00.set_ylabel(r'$y_2$'); ax_01.set_ylabel(r'$y_2$'); ax_10.set_ylabel(r'$y_2$'); ax_11.set_ylabel(r'$y_2$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A bivariate (two-dimensional) Jump-Diffusion process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more involved example is a bivariate Jump-Diffusion process, that takes the form\n",
    "$$\\begin{pmatrix}\n",
    "    \\mathrm{d}y_1(t) \\\\ \\mathrm{d}y_2(t)\n",
    "    \\end{pmatrix}=\n",
    "    \\begin{pmatrix}\n",
    "    N_1 \\\\ N_2\n",
    "    \\end{pmatrix}\n",
    "    \\mathrm{d} t +\n",
    "    \\begin{pmatrix}\n",
    "    g_{1,1} & g_{1,2} \\\\\n",
    "    g_{2,1} & g_{2,2}\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "    \\mathrm{d}W_1 \\\\ \\mathrm{d}W_2\n",
    "    \\end{pmatrix}+\n",
    "    \\begin{pmatrix}\n",
    "    \\xi_{1,1} & \\xi_{1,2} \\\\\n",
    "    \\xi_{2,1} & \\xi_{2,2}\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "    \\mathrm{d}J_1 \\\\ \\mathrm{d}J_2\n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "with $(N_1, N_2)$ the drift vector, $g$ the diffusion matrix, and $\\xi$ the jump matrix. The term $(\\mathrm{d}W_1(t),\\mathrm{d}W_2(t))$ comprises the Brownian motions and $(\\mathrm{d}J_1, \\mathrm{d}J_2)$ comprises the Poisson jumps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrating with an Euler−Maruyama scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define first the drift vector N with some exemplary values\n",
    "N = np.array([1.0, 0.5])\n",
    "\n",
    "# Define the diffusion matrix g with some exemplary values\n",
    "g = np.array([[0.5, 0.2],[0.1, 1.0]])\n",
    "\n",
    "# and the jump amplitude matrix xi\n",
    "xi = np.array([[0.2, 0.5],[0.4, 0.5]])\n",
    "\n",
    "# lastly, the jump frequency lambda (l)\n",
    "l = np.array([0.1, 0.2])\n",
    "\n",
    "# The total integration time\n",
    "t_start = 0\n",
    "t_final = 500\n",
    "\n",
    "# The desired timestep of integration\n",
    "delta_t = 0.001\n",
    "\n",
    "# time array of the process\n",
    "time = np.linspace(0, t_final, t_final * int(1 / delta_t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialise an (n,2)-array $y$, (n,2)-normally distributed noise $\\mathrm{d}W_{1,2}(t)$, and (n,2)-Poisson processes $\\mathrm{d}J_{1,2}(t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the array y\n",
    "y = np.zeros([time.size, 2])\n",
    "\n",
    "# Generate two Wiener processes with a scale of np.sqrt(delta_t)\n",
    "dW = np.random.normal(loc = 0, scale = np.sqrt(delta_t), size = [time.size, 2])\n",
    "\n",
    "# Generate two Poisson process with a delta_t scale\n",
    "dJ = np.zeros([ time.size, 2 ])\n",
    "dJ[:,0] = np.random.poisson(lam = l[0] * delta_t, size = time.size)\n",
    "dJ[:,1] = np.random.poisson(lam = l[1] * delta_t, size = time.size)\n",
    "\n",
    "# Give some small random initial conditions\n",
    "y[0,:]=np.random.normal(size = 2) / 10\n",
    "\n",
    "# Integrate the process, with a two step process\n",
    "for i in range(1, time.size):\n",
    "    y[i,0] = y[i-1,0]  -  N[0] * y[i-1,0] * delta_t  +  g[0,0] * dW[i,0] +  g[0,1] * dW[i,1]\n",
    "    y[i,1] = y[i-1,1]  -  N[1] * y[i-1,1] * delta_t  +  g[1,0] * dW[i,0] +  g[1,1] * dW[i,1]\n",
    "    \n",
    "    if dJ[i,0]>0 or dJ[i,1]:\n",
    "        # Generate here the normally distributed matrix xi of the jumps\n",
    "        s_00 = np.random.normal(loc = 0, scale = np.sqrt(xi[0,0]))\n",
    "        s_01 = np.random.normal(loc = 0, scale = np.sqrt(xi[0,1]))\n",
    "        s_10 = np.random.normal(loc = 0, scale = np.sqrt(xi[1,0]))\n",
    "        s_11 = np.random.normal(loc = 0, scale = np.sqrt(xi[1,1]))\n",
    "        \n",
    "        y[i,0] += s_00 * dJ[i,0]  +  s_01 * dJ[i,1]\n",
    "        y[i,1] += s_10 * dJ[i,0]  +  s_11 * dJ[i,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still in progress, now the Kramers−Moyal coefficients in two dimensions"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}