Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Milstein method
open-in-new
<h2 id="description">Description</h2> <p>Consider the <a href="/facts/Autonomous_equation/ahlr3A6o">autonomous</a> <a href="/facts/It%C5%8D_calculus/2uSR02fY">Itō</a> stochastic differential equation: d X t = a ( X t ) d t + b ( X t ) d W t {\displaystyle \mathrm {d} X_{t}=a(X_{t})\,\mathrm {d} t+b(X_{t})\,\mathrm {d} W_{t}} with <a href="/facts/Initial_condition/sA4ohFP4">initial condition</a> X 0 = x 0 {\displaystyle X_{0}=x_{0}} , where W t {\displaystyle W_{t}} denotes the <a href="/facts/Wiener_process/KPh7vYd7">Wiener process</a>, and suppose that we wish to solve this SDE on some interval of time [ 0 , T ] {\displaystyle [0,T]} . Then the Milstein approximation to the true solution X {\displaystyle X} is the <a href="/facts/Markov_chain/5oP5hntk">Markov chain</a> Y {\displaystyle Y} defined as follows: </p> <ul><li>Partition the interval [ 0 , T ] {\displaystyle [0,T]} into N {\displaystyle N} equal subintervals of width Δ t > 0 {\displaystyle \Delta t>0} : 0 = τ 0 < τ 1 < ⋯ < τ N = T with τ n := n Δ t and Δ t = T N {\displaystyle 0=\tau _{0}<\tau _{1}<\dots <\tau _{N}=T{\text{ with }}\tau _{n}:=n\Delta t{\text{ and }}\Delta t={\frac {T}{N}}} </li> <li>Set Y 0 = x 0 ; {\displaystyle Y_{0}=x_{0};} </li> <li>Recursively define Y n {\displaystyle Y_{n}} for 1 ≤ n ≤ N {\displaystyle 1\leq n\leq N} by: Y n + 1 = Y n + a ( Y n ) Δ t + b ( Y n ) Δ W n + 1 2 b ( Y n ) b ′ ( Y n ) ( ( Δ W n ) 2 − Δ t ) {\displaystyle Y_{n+1}=Y_{n}+a(Y_{n})\Delta t+b(Y_{n})\Delta W_{n}+{\frac {1}{2}}b(Y_{n})b'(Y_{n})\left((\Delta W_{n})^{2}-\Delta t\right)} where b ′ {\displaystyle b'} denotes the <a href="/facts/Derivative/7Ito70Dh">derivative</a> of b ( x ) {\displaystyle b(x)} with respect to x {\displaystyle x} and: Δ W n = W τ n + 1 − W τ n {\displaystyle \Delta W_{n}=W_{\tau _{n+1}}-W_{\tau _{n}}} are <a href="/facts/Independent_and_identically_distributed/othIRaWt">independent and identically distributed</a> <a href="/facts/Normal_distribution/UapjjPyQ">normal random variables</a> with <a href="/facts/Expected_value/1XV0JKL8">expected value</a> zero and <a href="/facts/Variance/ULBJKXD1">variance</a> Δ t {\displaystyle \Delta t} . Then Y n {\displaystyle Y_{n}} will approximate X τ n {\displaystyle X_{\tau _{n}}} for 0 ≤ n ≤ N {\displaystyle 0\leq n\leq N} , and increasing N {\displaystyle N} will yield a better approximation.</li></ul> <p>Note that when b ′ ( Y n ) = 0 {\displaystyle b'(Y_{n})=0} (i.e. the diffusion term does not depend on X t {\displaystyle X_{t}} ) this method is equivalent to the <a href="/facts/Euler%E2%80%93Maruyama_method/HL4K30c4">Euler–Maruyama method</a>. </p><p>The Milstein scheme has both weak and strong order of convergence Δ t {\displaystyle \Delta t} which is superior to the <a href="/facts/Euler%E2%80%93Maruyama_method/HL4K30c4">Euler–Maruyama method</a>, which in turn has the same weak order of convergence Δ t {\displaystyle \Delta t} but inferior strong order of convergence Δ t {\displaystyle {\sqrt {\Delta t}}} .<a class="footnote-ref" id="fnref:3" href="#fn:3"><sup>3</sup></a> </p> <h2 id="intuitive-derivation">Intuitive derivation</h2> <p>For this derivation, we will only look at <a href="/facts/Geometric_Brownian_motion/0Ve9dW6c">geometric Brownian motion</a> (GBM), the stochastic differential equation of which is given by: d X t = μ X d t + σ X d W t {\displaystyle \mathrm {d} X_{t}=\mu X\mathrm {d} t+\sigma XdW_{t}} with real constants μ {\displaystyle \mu } and σ {\displaystyle \sigma } . Using <a href="/facts/It%C5%8D%27s_lemma/NF7tTcgd">Itō's lemma</a> we get: d ln X t = ( μ − 1 2 σ 2 ) d t + σ d W t {\displaystyle \mathrm {d} \ln X_{t}=\left(\mu -{\frac {1}{2}}\sigma ^{2}\right)\mathrm {d} t+\sigma \mathrm {d} W_{t}} </p><p>Thus, the solution to the GBM SDE is: X t + Δ t = X t exp { ∫ t t + Δ t ( μ − 1 2 σ 2 ) d t + ∫ t t + Δ t σ d W u } ≈ X t ( 1 + μ Δ t − 1 2 σ 2 Δ t + σ Δ W t + 1 2 σ 2 ( Δ W t ) 2 ) = X t + a ( X t ) Δ t + b ( X t ) Δ W t + 1 2 b ( X t ) b ′ ( X t ) ( ( Δ W t ) 2 − Δ t ) {\displaystyle {\begin{aligned}X_{t+\Delta t}&=X_{t}\exp \left\{\int _{t}^{t+\Delta t}\left(\mu -{\frac {1}{2}}\sigma ^{2}\right)\mathrm {d} t+\int _{t}^{t+\Delta t}\sigma \mathrm {d} W_{u}\right\}\\&\approx X_{t}\left(1+\mu \Delta t-{\frac {1}{2}}\sigma ^{2}\Delta t+\sigma \Delta W_{t}+{\frac {1}{2}}\sigma ^{2}(\Delta W_{t})^{2}\right)\\&=X_{t}+a(X_{t})\Delta t+b(X_{t})\Delta W_{t}+{\frac {1}{2}}b(X_{t})b'(X_{t})((\Delta W_{t})^{2}-\Delta t)\end{aligned}}} where a ( x ) = μ x , b ( x ) = σ x {\displaystyle a(x)=\mu x,~b(x)=\sigma x} </p><p>The numerical solution is presented in the graphic for three different trajectories.<a class="footnote-ref" id="fnref:4" href="#fn:4"><sup>4</sup></a> </p> <h3>Computer implementation</h3> <p>The following <a href="/facts/Python_(programming_language)/YbuGqofa">Python</a> code implements the Milstein method and uses it to solve the SDE describing geometric Brownian motion defined by { d Y t = μ Y d t + σ Y d W t Y 0 = Y init {\displaystyle {\begin{cases}dY_{t}=\mu Y\,{\mathrm {d} }t+\sigma Y\,{\mathrm {d} }W_{t}\\Y_{0}=Y_{\text{init}}\end{cases}}} </p> # -*- coding: utf-8 -*- # Milstein Method import numpy as np import matplotlib.pyplot as plt class Model: """Stochastic model constants.""" mu = 3 sigma = 1 def dW(dt): """Random sample normal distribution.""" return np.random.normal(loc=0.0, scale=np.sqrt(dt)) def run_simulation(): """ Return the result of one full simulation.""" # One second and thousand grid points T_INIT = 0 T_END = 1 N = 1000 # Compute 1000 grid points DT = float(T_END - T_INIT) / N TS = np.arange(T_INIT, T_END + DT, DT) Y_INIT = 1 # Vectors to fill ys = np.zeros(N + 1) ys[0] = Y_INIT for i in range(1, TS.size): t = (i - 1) * DT y = ys[i - 1] dw = dW(DT) # Sum up terms as in the Milstein method ys[i] = y + \ Model.mu * y * DT + \ Model.sigma * y * dw + \ (Model.sigma**2 / 2) * y * (dw**2 - DT) return TS, ys def plot_simulations(num_sims: int): """Plot several simulations in one image.""" for _ in range(num_sims): plt.plot(*run_simulation()) plt.xlabel("time (s)") plt.ylabel("y") plt.grid() plt.show() if __name__ == "__main__": NUM_SIMS = 2 plot_simulations(NUM_SIMS) <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Euler%E2%80%93Maruyama_method/HL4K30c4">Euler–Maruyama method</a></li></ul> <h2 id="further-reading">Further reading</h2> <ul><li>Kloeden, P.E., & Platen, E. (1999). <i>Numerical Solution of Stochastic Differential Equations</i>. Springer, Berlin. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 3-540-54062-8.{{cite book}}: CS1 maint: multiple names: authors list (link)</li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Mil'shtein, G. N. (1974). "Approximate integration of stochastic differential equations". Teoriya Veroyatnostei i ee Primeneniya (in Russian). 19 (3): 583–588. <a href="https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=tvp&paperid=2929&option_lang=eng" target="_blank">https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=tvp&paperid=2929&option_lang=eng</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> <li id="fn:2"><p>Mil’shtein, G. N. (1975). "Approximate Integration of Stochastic Differential Equations". Theory of Probability & Its Applications. 19 (3): 557–000. doi:10.1137/1119062. <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li> <li id="fn:3"><p>V. Mackevičius, Introduction to Stochastic Analysis, Wiley 2011 <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li> <li id="fn:4"><p>Umberto Picchini, SDE Toolbox: simulation and estimation of stochastic differential equations with Matlab. http://sdetoolbox.sourceforge.net/ <a href="http://sdetoolbox.sourceforge.net/" target="_blank">http://sdetoolbox.sourceforge.net/</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li> </ol>