Programming one SDE trajectory

A reliable stochastic program should make its inputs, random generator, numerical loop, stored path and validity diagnostics easy to inspect.

Programming goal

The model equations and Euler–Maruyama method are already known. The new task is to package one complete simulation so that inputs, random-number generation, state storage, validation and outputs are clearly separated.

Scenario and purpose

We reuse the stochastic SIS model with \(N=500\), \(I(0)=10\), \(\beta=0.25\) and \(\gamma=0.15\). The purpose is no longer to derive its coefficients. It is to build a dependable function that produces one possible epidemic history.

One trajectory is one sequence of numerical states generated from one sequence of Brownian increments. It is a possible outcome under the model, not the expected path and not a probability distribution.

Inputs and outputs

Inputs	Returned outputs
Initial state and parameters	Time grid and state path
Final time and time step	Diagnostic information
A random-number generator	Values ready for tables and plots

Passing the generator into the function avoids hidden randomness and makes later repeated simulations straightforward.

The function contract

A function contract states what a function expects and what it returns.

Requirement	Reason
\(N>0\), \(T>0\) and \(\Delta t>0\)	Population, horizon and step length must be meaningful.
\(\beta,\gamma\ge0\)	Biological event rates cannot be negative.
\(0\le I_0\le N\)	The initial state must lie in the SIS state space.
\(T/\Delta t\) is an integer	The stored grid reaches the requested final time exactly.

If a requirement fails, the function raises a clear error instead of returning a misleading trajectory.

Interactive Python laboratory

Interactive PythonOne reusable SDE trajectory

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def simulate_sis_sde(
    rng, N, I0, beta, gamma, T, dt
):
    if N <= 0:
        raise ValueError("N must be positive.")

if T <= 0 or dt <= 0:
        raise ValueError("T and dt must be positive.")

if beta < 0 or gamma < 0:
        raise ValueError("Rates cannot be negative.")

if not 0 <= I0 <= N:
        raise ValueError("I0 must be between 0 and N.")

steps = int(round(T / dt))

if not np.isclose(steps * dt, T):
        raise ValueError(
            "Choose dt so that T/dt is an integer."
        )
    times = np.linspace(0, T, steps + 1)
    I = np.zeros(steps + 1)
    I[0] = I0
    corrections = 0

for n in range(steps):
        current = I[n]
        infection_rate = beta * (N - current) * current / N
        recovery_rate = gamma * current

drift = infection_rate - recovery_rate
        diffusion = np.sqrt(max(
            infection_rate + recovery_rate, 0
        ))
        dW = np.sqrt(dt) * rng.normal()

proposal = (
            current + drift * dt + diffusion * dW
        )
        I[n + 1] = np.clip(proposal, 0, N)

if I[n + 1] != proposal:
            corrections += 1

diagnostics = {
        "steps": steps,
        "boundary corrections": corrections,
        "minimum I": I.min(),
        "maximum I": I.max()
    }

return times, I, diagnostics

# Create the generator once, outside the function
seed = 71
rng = np.random.default_rng(seed)

times, I, diagnostics = simulate_sis_sde(
    rng=rng,
    N=500,
    I0=10,
    beta=0.25,
    gamma=0.15,
    T=80,
    dt=0.01
)

trajectory = pd.DataFrame({
    "time": times,
    "I": I
})

print("Diagnostics")
print(pd.Series(diagnostics).round(4).to_string())
print()
print("First 10 stored states")
print(trajectory.head(10).round(4).to_string(index=False))
print()
print("Last 10 stored states")
print(trajectory.tail(10).round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.plot(trajectory["time"], trajectory["I"],
         color="#b42318", linewidth=1.1)
plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title(f"One SIS SDE trajectory, seed = {seed}")
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Why each design choice matters

Choice	Benefit
Function parameters instead of global variables	The simulation can be reused safely with new scenarios.
Generator passed as `rng`	The source of randomness is visible and controllable.
Preallocated NumPy array	Storage is efficient and has a predictable length.
Returned diagnostics	Boundary corrections and validity can be inspected rather than hidden.

Seed, generator and reproducibility

The seed initialises the random-number generator. Recreating a generator with the same seed and using the same code reproduces the same trajectory. Continuing to use an existing generator produces new draws because its internal state has advanced.

Action	Result
Create a new generator with seed 71	Reproduces this stored path.
Call the function again with the same existing generator	Produces another path from later random draws.
Create a new generator with a different seed	Produces a different reproducible path.

A fixed seed supports checking and teaching. It does not make the stochastic model deterministic.

What one path can establish

One path demonstrates a possible history and helps debug the program. It cannot estimate a probability, mean, variance or typical range. Those require independent repetitions.

Continuous model, discrete numerical storage

The SDE defines a continuous-time process. Euler–Maruyama stores approximations only at grid times. A line plot connects those stored values; it is not evidence that the computer calculated every intermediate time.

What this lesson adds

You can now package one SDE path into a validated reusable function, pass randomness explicitly, return diagnostics, and distinguish a possible trajectory from an ensemble result.