renewal model time series example for User's Guide

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

Start by reviewing the User's Guide structure and the Stan model and CmdStanPy simulator included in the issue. Add the renewal-model example with its simulation code and the cited Fraser and Steyn references, then verify that the example is placed in the appropriate guide section and the links and code are usable.