# Stochastic Differential Euations

## February 10/17, 2023

https://hackmd.io/@ubcs3

The paper that Patrick found: https://epubs.siam.org/doi/10.1137/S0036144500378302

@article{doi:10.1137/S0036144500378302,
author = {Higham., Desmond J.},
title = {An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations},
journal = {SIAM Review},
volume = {43},
year = {2001},
doi = {10.1137/S0036144500378302},
URL = {https://doi.org/10.1137/S0036144500378302},
}

import numpy as np
import matplotlib.pyplot as plt

T = 1; N = 100; dt = T/N; sigma = np.sqrt(dt); mu = 0; npaths = 5;
dW = sigma*np.random.randn(N,npaths) + mu
X = np.vstack([np.zeros(npaths),dW.cumsum(axis=0)])
t = np.linspace(0,T,N+1)
plt.plot(t,X,'b-',alpha=0.5), plt.grid(True)
plt.show()

T = 1; N = 1000; dt = T/N; R = 1; dtau = R*dt; 
dW = np.sqrt(dt)*np.random.randn(N,1)
W = np.vstack([np.zeros(1),dW.cumsum(axis=0)])
t = np.linspace(0,T,N+1)
X = np.zeros(N+1)
a = 1; b = 2;
X[0] = 1
for j in range(1,N+1):
    X[j] = X[j-1] + a*X[j-1]*dtau + b*X[j-1]*dW[j-1]
plt.plot(X)
plt.show()

def sdeEM(f,g,X0=1,T=1,R=10,L=100,nsamples=20):
    N = R*L; dt = T/N; Dt = R*dt;
    dW = np.sqrt(dt)*np.random.randn(N,nsamples)
    W = np.vstack([np.zeros(nsamples),dW.cumsum(axis=0)])
    X = np.zeros((L+1,nsamples))
    X[0,:] = X0*np.ones(nsamples)
    for j in range(1,L+1):
        X[j,:] = X[j-1,:] + f(X[j-1,:])*Dt + g(X[j-1,:])*(W[j*R,:] - W[(j-1)*R,:])
    t = np.linspace(0,T,L+1)
    return (t,X,W)

alpha = 2; beta = 0.2; 
f = lambda x: alpha - x
g = lambda x: beta*np.sqrt(np.abs(x))
t,X,W = sdeEM(f,g,T=5)
plt.plot(t,X,'b',alpha=0.2)
plt.plot(t,X.mean(axis=1),'r')
plt.show()