We simulate a one dimensional Ornstein-Uhlenbeck process by using an Euler approximation of the solution to the corresponding stochastic differential equation.
# Euler-Maryuama Approximation for solution of SDE dX_t = -gamma X_t dt + sigma dB_t
import numpy as np
# makes numpy routines and data types available as np.[name ouf routine or data type]
import matplotlib.pyplot as plt
# makes plotting command available as plt.[name of command]
x0 = 5.
# initial value at time t=0
sigma = 2.
# volatility
gamma = .5
# friction coefficient
tmax = 5.
# simulation of process from time 0 to tmax
stepslist = [10,100,1000,10000]
# produce simulations with step numbers chosen from steplist
for steps in stepslist:
h = tmax/steps
# stepsize for time discretization
std = np.sqrt(h)
# standard deviation for the distribution of each step
k = 100
# number of samples that will be simulated
noise = np.random.randn(steps,k)*std
# create a steps times k dimensional matrix of normal random numbers
# with variance h
sde = np.ones((steps,k))
sde = sde*x0
# initialize the array sde with the initial condition x0
for n in range(steps-1):
sde[n+1] = sde[n]-gamma*h*sde[n]+sigma*noise[n]
t = np.arange(0,steps,1)*h
# creates vector of time points
plt.figure(figsize=(10,6), dpi=80)
# sets size of plot
plt.plot(t,sde,linewidth=0.3)
# produces a plot of the sample solutions in the components of sde
# versus t with thin lines
plt.show()
# output of plot