Monte Carlo Estimation of Renewal Functions

From Renewal Sequences to an Estimate

For a renewal model, the renewal function \(m(t)=\mathbb{E}N(t)\) is an expected count. In many models its exact value is not available in a convenient form, but one can simulate the interarrival times, count the renewals by time \(t\), and average the counts across independent repetitions. The simulation must implement the endpoint convention correctly: an epoch at exactly \(t\) is counted, while the first epoch after \(t\) is not.

Let \(X_1,X_2,\ldots\) be iid, strictly positive interarrival times on a probability space \((\Omega,\mathcal{F},\mathbb{P})\), with common distribution \(F\). Define the renewal epochs \(S_n=X_1+\cdots+X_n\) and the renewal count \(N(t)=\sum_{n\ge1}\mathbf{1}_{\{S_n\le t\}}\), for \(t\ge0\). The renewal function is \(m(t)=\mathbb{E}N(t)\). The sum uses \(\le t\), so an epoch at the endpoint is included. A simulation stops when it first observes an epoch strictly greater than \(t\); that overshooting interarrival is needed to know that the process has stopped, but its epoch does not contribute to \(N(t)\).

To simulate an interarrival time from \(F\), one may use an independent uniform random variable \(U\) on \((0,1)\) and the generalized inverse \(F^{-1}(u)=\inf\{x:F(x)\ge u\}\). In the examples below, the distribution is simple enough to state the transformation directly. A fixed list of uniform inputs makes the arithmetic reproducible; when inputs are instead generated randomly and independently, the sampling results below apply.

Why the Stopping Rule Gives the Renewal Count

Proof. Since \(X_1>0\) almost surely, there is an \(\varepsilon>0\) such that \(p=\mathbb{P}(X_1\ge\varepsilon)>0\). Indeed, the union of the events \(\{X_1\ge1/k\}\), over positive integers \(k\), has probability one, so at least one such event has positive probability. Let \(I_j=\mathbf{1}_{\{X_j\ge\varepsilon\}}\). The \(I_j\) are iid Bernoulli variables with success probability \(p\). With probability one there are infinitely many successes: for any fixed integer \(r\), the probability of at most \(r\) successes among the first \(n\) trials is \(\sum_{k=0}^r {n\choose k}p^k(1-p)^{n-k}\), which tends to zero as \(n\) tends to infinity when \(0<p<1\); if \(p=1\), every trial is a success. Thus \(S_n\ge\varepsilon\sum_{j=1}^n I_j\) tends to infinity almost surely, and some finite \(S_n\) exceeds \(t\). Because every \(X_j\) is positive, the sequence \(S_n\) is strictly increasing. If \(S_q\) is its first epoch exceeding \(t\), then exactly \(S_1,\ldots,S_{q-1}\) are at most \(t\), so the algorithm reports \(q-1=\sum_{n\ge1}\mathbf{1}_{\{S_n\le t\}}=N(t)\). \(\square\)

The count is also suitable for estimating an expectation and a variance under these assumptions. In particular, its second moment is finite for each fixed \(t\). To see why, choose \(\varepsilon\) and \(p\) as in the proof, and put \(r=\lfloor t/\varepsilon\rfloor\). If \(N(t)\ge n\), then \(S_n\le t\), so at most \(r\) of the first \(n\) interarrivals can be at least \(\varepsilon\). For \(0<p<1\), this gives \(\mathbb{P}(N(t)\ge n)\le\sum_{k=0}^r{n\choose k}p^k(1-p)^{n-k}\). The right side is bounded by a polynomial in \(n\) times a geometrically decreasing factor, so its sum weighted by \(2n-1\) is finite. If \(p=1\), \(N(t)\le r\) almost surely. Finally, for any nonnegative integer-valued \(N\), \(N^2=\sum_{n=1}^N(2n-1)\), and hence \(\mathbb{E}[N(t)^2]=\sum_{n\ge1}(2n-1)\mathbb{P}(N(t)\ge n)<\infty\).

Proof. Linearity of expectation gives \(\mathbb{E}[\widehat m_M(t)]=M^{-1}\sum_{r=1}^M\mathbb{E}[Y_r]=m(t)\). Independence gives \(\operatorname{Var}(\widehat m_M(t))=M^{-2}\sum_{r=1}^M\operatorname{Var}(Y_r)=\sigma^2/M\), because the covariance between distinct \(Y_r\) is zero. To verify the sample-variance claim, use the identity \(\sum_{r=1}^M(Y_r-\widehat m_M)^2=\sum_{r=1}^M(Y_r-m(t))^2-M(\widehat m_M-m(t))^2\). Taking expectations, the first sum on the right has expectation \(M\sigma^2\), and the second term has expectation \(M(\sigma^2/M)=\sigma^2\). The expected left side is therefore \((M-1)\sigma^2\), which proves \(\mathbb{E}[s_M^2]=\sigma^2\). The estimated standard error follows by taking the square root of the estimated variance \({s_M^2}/{M}\). \(\square\)

The standard error describes the estimated spread of the average across repeated simulations; it is not the actual error \(\lvert\widehat m_M(t)-m(t)\rvert\), which is unknown unless the true value is known. A small estimated standard error is useful evidence of numerical precision, but it is not a guarantee that one particular run is close to the truth.

Worked Simulations

Common Pitfalls and Practical Use

Common mistake: count only epochs strictly before \(t\). In replication 2 at \(t=5\), the epochs are \(2,3,5,6\). Counting only those with \(S_n<5\) gives \(2\), but the renewal-count convention is \(S_n\le5\), so the correct count is \(3\). The equality case matters for distributions with atoms, such as the two-point law used here.

Common mistake: keep drawing after the stopping epoch. If the first epoch above \(5\) is the fourth epoch, the simulation has already stopped after the fourth interarrival. Any later values are not used and must not appear in a column labelled “draws used.” The overshooting draw is included because it identifies the stopping point; it is not counted as a renewal by time \(t\).

In practice, keep the tasks distinct: generate one independent sequence per replication, stop each sequence at its first epoch strictly greater than \(t\), and record the number of preceding epochs. Then average those counts and report the sample standard deviation divided by \(\sqrt M\) as the estimated standard error. Increasing \(M\) generally reduces the standard deviation of the average at the rate \(1/\sqrt M\), as the variance formula shows, though it does not change the variance of an individual renewal count.