Deriving second-moment equations

Second moments follow from how each event changes the square of the stochastic state.

Squaring changes the jump contribution

For current state \(i\), a birth changes \(i^2\) by

\[(i+1)^2-i^2=2i+1,\]

and a removal changes it by

\[(i-1)^2-i^2=-2i+1.\]

Derive the second-moment equation

Multiply each squared jump change by its rate and take expectation:

\[ \frac{d}{dt}\mathbb E[I^2] = \mathbb E[(2I+1)\lambda I+(-2I+1)\mu I]. \]

With \(m_1=\mathbb E[I]\) and \(m_2=\mathbb E[I^2]\),

\[ \frac{dm_2}{dt} = 2(\lambda-\mu)m_2+(\lambda+\mu)m_1. \]

This equation is closed together with \(m_1'=(\lambda-\mu)m_1\).

Obtain the variance equation

\[v=m_2-m_1^2.\]

Differentiating and substituting the moment equations gives

\[ \frac{dv}{dt} = 2(\lambda-\mu)v+(\lambda+\mu)m_1. \]

The additional \((\lambda+\mu)m_1\) term creates stochastic variation even when the initial state has zero variance.

Interactive Python laboratory

Interactive PythonSecond moment and variance

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

I0 = 5
birth = 0.18
death = 0.12
growth = birth - death
times = np.linspace(0, 20, 81)

# Exact moments of the linear birth–death process
mean_exact = I0 * np.exp(growth * times)
variance_exact = (
    I0 * (birth + death) / growth
    * np.exp(growth * times)
    * (np.exp(growth * times) - 1)
)
second_exact = variance_exact + mean_exact ** 2

def path_on_grid(rng):
    I = I0
    time = 0.0
    event_times = [0.0]
    values = [I]
    while I > 0 and time < times[-1]:
        wait = rng.exponential(1 / ((birth + death) * I))
        if time + wait > times[-1]:
            break
        time += wait
        I += 1 if rng.random() < birth / (birth + death) else -1
        event_times.append(time)
        values.append(I)
    positions = np.searchsorted(
        event_times, times, side="right"
    ) - 1
    return np.array(values)[positions]

rng = np.random.default_rng(260)
paths = np.array([path_on_grid(rng) for _ in range(5000)])
mean_mc = paths.mean(axis=0)
second_mc = np.mean(paths ** 2, axis=0)
variance_mc = paths.var(axis=0, ddof=1)

selected = [0, 5, 10, 15, 20]
idx = [np.argmin(np.abs(times - t)) for t in selected]
table = pd.DataFrame({
    "time": selected,
    "exact second moment": second_exact[idx],
    "MC second moment": second_mc[idx],
    "exact variance": variance_exact[idx],
    "MC variance": variance_mc[idx]
})
print(table.round(3).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, second_exact, color="#174c34",
         label="Exact second moment")
plt.plot(times, second_mc, "--", color="#39815e",
         label="Monte Carlo second moment")
plt.plot(times, variance_exact, color="#b42318",
         label="Exact variance")
plt.plot(times, variance_mc, "--", color="#d97706",
         label="Monte Carlo variance")
plt.xlabel("Time")
plt.ylabel("Moment value")
plt.title("Second-moment equation and simulations")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Why the second moment is not the variance

\(m_2\) includes both squared mean and spread. Variance removes the squared mean. Always calculate \(v=m_2-m_1^2\) before interpreting uncertainty.

What this lesson adds

You can now derive a second-moment equation from squared jump changes, obtain the variance equation and explain how stochastic events generate variance from a deterministic initial condition.