import numpy as np import matplotlib.pyplot as plt

1. Parameters

N = 10000 # Number of simulations mu = 0.1 threshold = 1.0 sigma = 0.2 dt = 0.01 max_wait = 40

1. Plotting

fig, axes = plt.subplots(2, 1, figsize=(12, 12))

1. Subplot 1: Trajectories

ax = axes[1] ax.plot([0, max_wait], [1, 1], 'r--') waiting_times = [] for _ in range(N):

   X = 0
   t = 0
   X_values = [X]
   t_values = [t]
   while X < threshold and t < max_wait:
       X += mu * dt + sigma * np.random.normal(0, np.sqrt(dt))
       t += dt
       X_values.append(X)
       t_values.append(t)
   if t < max_wait:
     waiting_times.append(t)
     ax.plot(t_values, X_values, alpha=20/N, c='b')

ax.set_xlim(0, max_wait) ax.set_ylim(-1, threshold + 0.02) ax.set_xlabel('Time') ax.set_ylabel('X') ax.set_title(rf'Brownian motion with drift. (N={N}, $\mu$={mu}, $\sigma$={sigma})')

1. Subplot 2: Histogram

ax = axes[0] ax.hist(waiting_times, bins=100, density=True) ax.set_xlabel('Waiting Time') ax.set_ylabel('Frequency') ax.set_title('Waiting time distribution')

def f(x, mu, lambda_):

 return np.sqrt(lambda_ / (2 * np.pi * x**3)) * np.exp(-(lambda_ * (x - mu)**2) / (2 * mu**2 * x))

x = np.linspace(0.1, 40, 400) y = f(x, threshold / mu, (threshold / sigma)**2) ax.plot(x, y, 'r--')

plt.tight_layout() plt.savefig("inverse_gaussian.png") plt.show()