# Programming the Milstein method and computing the increments

In the wikipedia article on the Milstein method, the following python code to simulate a geometric Brownian motion is presented:

import numpy as np
import matplotlib.pyplot as plt

num_sims = 1  # One Example

# 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

## Initial Conditions
y_init = 1
mu    = 3
sigma = 1

# dw Random process
def dW(delta_t):
"""Random sample normal distribution"""
return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))

# vectors to fill
ts = np.arange(t_init, t_end + dt, dt)
ys = np.zeros(N + 1)
ys[0] = y_init

# Loop
for _ in range(num_sims):
for i in range(1, ts.size):
t = (i - 1) * dt
y = ys[i - 1]
# Milstein method
ys[i] = y + mu * dt * y + sigma* y* dW(dt) + 0.5* sigma**2 * y* (dW(dt)**2 - dt)
plt.plot(ts, ys)

# Plot
plt.xlabel("time (s)")
plt.grid()
h = plt.ylabel("y")
h.set_rotation(0)
plt.show()


I am wondering, in the line:

ys[i] = y + mu * dt * y + sigma* y* dW(dt) + 0.5* sigma**2 * y* (dW(dt)**2 - dt)


should this not be

rd = dW(dt)
ys[i] = y + mu * dt * y + sigma* y* rd + 0.5* sigma**2 * y* (rd**2 - dt)


i.e. using the same increment twice. Or is it really possible to compute two (in general different) increments in one step?

• Your observation seems pertinent and intelligent to me. Jul 13, 2021 at 0:17