Repository files navigation

The aim of this project was to apply Itô's stochastic differential equations (SDEs) theory for the particular problem from population biology and to compare a stochastic competition model model with a deterministic one. In order to gain better understanding of this theory, several problems were covered:

1. Stochastic process modelling.
2. Itô's stochastic integral modelling and approximate calculation of its expectations.
3. Obtaining numerical solutions to stochastic differential equation (SDE) using Monte-Carlo simulation.

Let F be a function F(t, X(t)) and X be a stochastic process such that

then (assuming the necessary assumptions) function F satisfies the following SDE:

Itô's lemma allows us to find closed-form solution to some SDEs and exactly calculate some stochastic integrals.

Monte Carlo simulation is a computational algorithm that relies on random sampling to obtain numerical result. Pseude-random numbers obtained from Numpy generator were used to simulate randomness.

Most of the project (except population dynamics modelling) was implemented using laptop with the following resources:Intel(R) Core(TM) i5-8300H CPU @ 2.30Ghz, 8GB RAM, Windows 10

Python and libraries:

• Python 3.8.5
• Numpy 1.19.1
• Scipy 1.5.0
• Matplotlib 3.3.1
• Seaborn 0.11.0

Population dynamics modelling was implemented inGoogle colab (with free access plan). Language and libraries versions are the following:

• Python 3.6.9
• Numpy 1.19.4
• Scipy 1.4.1
• Matplotlib 3.2.2
• Seaborn 0.11.0
• Numba 0.48.0

There's a small library for stochastic modelling and calculus: StoCalc.py. It is built atop of Scipy and Numpy and provides functions for modelling stochastic processes and integrals, calculation of stochastic integral expectations, and solving SDEs.

Obtain one Weiner process sample path on interval [0, 50] with step = 5:

importmatplotlib.pyplotaspltimportnumpyasnpfromStoCalcimportweiner_processnp.random.seed(0)time1,proc1=weiner_process(t0=0,T=T,m=1,N=11)plt.plot(time1,proc1)

Add more points to this trajectory (so that step would be 0.5):

fromStoCalcimportbrownian_bridgenew_points=np.arange(0,T,0.5)time2,proc2=brownian_bridge(np.array([time1,proc1[0]]),new_points)plt.plot(time1,proc1[0],label='step=5')plt.plot(time2,proc2,label='step=0.5')

Consider the following Itô's integral:

Calculate first 3 moments, i.e. E(I(f)), E(I²(f)), E(I³(f)), with step=0.01 and different sample size:

fromStoCalcimportito_int_expectnp.random.seed(0)deff(t,Wt):returnt*np.exp(Wt)sample_sizes= (10,100,10000)ks= (1,2,3)# moment ordersforsizeinsample_sizes:tmp= [size]forkinks:tmp.append(ito_int_expect(f,0,1,k,m=size))print(*tmp,sep='\t')# output:# 10     0.124  1.163  0.106# 100    0.116  0.729  0.859# 10000  0.009  1.291  5.864

Generate 1000 samples, plot the first 3 of them, average and variance of trajectories:

fromStoCalcimportito_int_pathsnp.random.seed(0)time,paths=ito_int_paths(f,0,1,step=0.01,m=1000)forpinpaths[:3]:plt.plot(time,p,linewidth=1)plt.plot(time,paths.mean(axis=0),linewidth=3,label='Average of 1000 trajectories')plt.plot(time,paths.var(ddof=1,axis=0),linewidth=3,label='Variance of 1000 trajectories')

Consider the following SDE:

Obtain 1000 solutions using Euler's and Milstein's method and compare average of 1000 solutions for both methods with Scipy-provided solver for a corresponing determenistic ODE:

fromStoCalcimportsdenp.random.seed(0)defdyE(y,dt,dW):f=1/4*y+1/32*y**2g=1/4*yreturn (f,g)defdyM(y,dt,dW):f=1/4*y+1/32*y**2g=1/4*ydgdy=1/4# need to provide additional derivative for Milstein's methodreturnf,g,dgdyy0=1/2step=1/5# step size is large to see the difference between methodst= (0,5)time,pathsE=sde(dyE,y0,t,step,method='euler')# Euler's methodtime,pathsM=sde(dyM,y0,t,step,method='milstein')# Milstein's method# determenistic equationfromscipy.integrateimportodeintdefdet_dydt(y,t):return1/4*y+1/32*y**2sol=odeint(det_dydt,y0,time)plt.plot(time,pathsE.mean(axis=0),color='blue',label="Euler's method")plt.plot(time,pathsM.mean(axis=0),color='red',label="Milstein's method")plt.plot(time,sol,color='green',label='Scipy determenistic solution')

Seeexamples for more.

Finally, Itô's theory was applied to the population biology system of two populations that are in competition (seeModelling.ipynb for the exact problem formulation).

The stochastic model was compared to the corresponding determenistic model, and the solutions did not exactly match. Most noticable is the fact that the mean population size estimation is much different from the determenistic result, because in the stochastic model there are 2 possible end states and both populations have a probability to go extinct. The stochastic model was used to estimate extinction probabilities, extinction time and conditional probability density for sizes of both populations. Mean conditional extincton time for the second population is also less then extinction time from the determenistic model because of the chance to survive:

SeeModelling.ipynb for conditional probability densities and modelling source code.

1. Allen, Edward. (2007). Modeling with Itô Stochastic Differential Equations. 22. 10.1007/978-1-4020-5953-7.
2. Choongbum Lee, Peter Kempthorne, Vasily Strela, Jake Xia. 18.S096 Topics in Mathematics with Applications in Finance. Fall 2013. Massachusetts Institute of Technology: MIT OpenCourseWare,https://ocw.mit.edu/. License: Creative Commons BY-NC-SA.