Movatterモバイル変換


[0]ホーム

URL:


Skip to content

Navigation Menu

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Sign up
This repository was archived by the owner on Feb 5, 2023. It is now read-only.
/SDEPublic archive

1st semester project at Bioinformatics institute, Fall 2020

NotificationsYou must be signed in to change notification settings

iam28th/SDE

Repository files navigation

Goals and objectives

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.

Methods

Itô's lemma

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

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.

System requirements

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

Results and examples

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.

1. Weiner process

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)

Brownian bridge 5

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')

Brownian bridge 0.5

2. Itô's stochastic integrals

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')

Integral

3. Stochastic differential equations

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.SDE

4. Population dynamics applications

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:mean_TE

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

References

  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.

Releases

No releases published

Packages

No packages published

[8]ページ先頭

©2009-2025 Movatter.jp