Using a Shared Brownian Increment in the Milstein Update

Summary

The document examines a coding question about simulating geometric Brownian motion with the Milstein method. Its update formula contains a Brownian increment in both the first-order stochastic term and the squared-increment correction. The author asks whether both appearances should use the same sampled increment for a given time step, or whether independently sampled increments can be used.

That is a consequential implementation detail: the Milstein correction is defined using the same step increment as the stochastic term, so the code should draw the increment once and reuse it in both places. Drawing twice would combine unrelated random values and would not implement the displayed Milstein update correctly. The document itself poses the question but supplies no answer, numerical comparison, or convergence evidence. It also gives no broader guidance on step-size selection or simulation validation, so its useful scope is limited to identifying the increment-consistency issue.

Key ideas

• The Milstein update includes a Brownian increment in both its linear stochastic term and its squared-increment correction.
• Both terms for a given time step should use the same sampled increment.
• Sampling independently for the two appearances does not implement the displayed Milstein formula.
• The document raises the implementation issue but provides no simulation results or convergence analysis.

Tags

# Programming the Milstein method and computing the increments


# 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?

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.