renewal model time series example for User's Guide

[bob-carpenter]

We should include a renewal model---they're popular for tracking infectious diseases. Here's the basic model with Stan code and CmdStanPy simulation.

functions {
  real random_walk_lpdf(vector y, real sigma, real mu0, real sigma0) {
    int N = rows(y);
    return normal_lpdf(y[1] | mu0, sigma0)
           + normal_lpdf(y[2:N] | y[1:N - 1], sigma);
  }
}
data {
  int<lower=0> T, U;
  simplex[U] omega;
  array[T]  int<lower=0> C;
  real mu0;
  real<lower=0> sigma0;
}
transformed data {
  row_vector<lower=0>[T] Cv = to_row_vector(C);
}
parameters {
  vector<lower=0>[T] R;
  real<lower=0> sigma;
}
model {
  // C[1] not modeled
  for (t in 2:U) {
    // shorten history with renormalized simplex
    C[t] ~ poisson(R[t] * (Cv[1:t-1] * omega[1:t-1]) / sum(omega[1:t-1]));
  }
  for (t in U + 1:T) {
    C[t] ~ poisson(R[t] * (Cv[t-U:t-1] * omega));
  }
  sigma ~ lognormal(0, 0.25);
  R ~ random_walk(sigma, mu0, sigma0);
}

And here's the simulator.

import numpy as np
import cmdstanpy as csp


def simulate_renewal_data(
    T: int = 150,
    U: int = 5,
    seed: int = 2020,
):
    rng = np.random.default_rng(seed)

    alpha = np.arange(U, 0, -1, dtype=float)
    omega = alpha / sum(alpha)

    mu0 = 1.0
    sigma0 = 0.05
    sigma = 0.025

    R = np.empty(T)
    R[0] = abs(rng.normal(mu0, sigma0))
    for t in range(1, T):
        R[t] = abs(rng.normal(R[t - 1], sigma))

    C = np.zeros(T, dtype=int)
    C[0] = 1000
    for t in range(1, U):
        lam = R[t] * np.dot(C[0:t], omega[0:t] / sum(omega[0:t]))
        C[t] = rng.poisson(lam)
    for t in range(U, T):
        lam = R[t] * np.dot(C[t-U:t], omega)
        C[t] = rng.poisson(lam)

    return {
        "T": T,
        "U": U,
        "omega": omega,
        "C": C,
        "R": R,
        "sigma": sigma,
        "mu0": mu0,
        "sigma0": sigma0,
        "alpha": alpha,
    }

sim = simulate_renewal_data()

print("alpha: ")
print(np.round(sim["alpha"], 3))
print("\nomega: ")
print(np.round(sim["omega"], 3))
print("\nsigma: ", sim["sigma"])
print("\nR: ")
print(np.round(sim["R"], 3))
print("\nC: ")
print(sim["C"])

model = csp.CmdStanModel(stan_file="poisson.stan")
fit = model.sample(data=sim)
print(fit.summary())

The model is taken from Steyn et al. (2025), Model 1 (section 6.1), which also introduces overdispersed alternatives (negative binomial instead of Poisson), latent infectiousness states, day-of-week effects, etc. They cite Fraser (2007) as the source of the basic model.

@article{steyn2025primer,
  title={A primer on inference and prediction with epidemic renewal models and sequential {M}onte {C}arlo},
  author={Steyn, Nicholas and Parag, Kris V and Thompson, Robin N and Donnelly, Christl A},
  journal={Statistics in Medicine},
  volume={44},
  number={18-19},
  pages={e70204},
  year={2025}
}

@article{fraser2007estimating,
  title={Estimating individual and household reproduction numbers in an emerging epidemic},
  author={Fraser, Christophe},
  journal={PloS one},
  volume={2},
  number={8},
  pages={e758},
  year={2007}
}

Here are links to the papers:

• https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.70204
• https://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0000758&type=printable