Why Simulate Histories Forward?
Under a fixed treatment regime, each covariate can depend on earlier covariates and treatments, while later covariates and outcomes can depend on the entire history generated so far. Forward simulation respects that order: it draws the next covariate from its conditional distribution given the history already available, then continues until an outcome is drawn. The result is a complete trajectory compatible with the fitted conditional models.
This tutorial focuses on generating trajectories conditional on a specified baseline value. It does not average over a target population’s baseline distribution. Also distinguish two questions: the simulation algorithm describes what follows from fitted models under a regime; interpreting that distribution causally additionally requires the identification assumptions, including consistency, sequential exchangeability, and positivity.
Set Up the Conditional Models
Let \(L_0\) denote baseline covariates. At times \(t=0,\ldots,K-1\), treatment \(A_t\) is followed by covariates \(L_{t+1}\); the terminal outcome is \(Y\). A fixed regime is the specified sequence \(\bar a=(a_0,\ldots,a_{K-1})\). Write \(\bar L_t=(L_0,\ldots,L_t)\) and \(\bar A_t=(A_0,\ldots,A_t)\). For a given baseline value \(l_0\), the fitted models supply conditional distributions \(\hat f_t(l_{t+1}\mid \bar l_t,\bar a_t)\) for covariate transitions and \(\hat f_Y(y\mid \bar l_K,\bar a_{K-1})\) for the outcome.
Here \(\hat f\) can be a probability mass function or a density, as appropriate. The treatment values in these models are set to the regime’s prescribed values; treatment is not drawn from a treatment-assignment model. Each fitted conditional distribution must be normalized for the histories at which it is used, and its parameters must be fixed during a simulated trajectory.
For a categorical covariate with conditional probabilities \(p_1,\ldots,p_J\), one can draw a uniform random variable \(U\) on \([0,1]\) and select category \(j\) when \(U\) falls in that category’s cumulative-probability interval. For a continuous covariate, draw from the fitted conditional distribution; for example, if \(L_{t+1}\mid H_t\sim N(\hat\mu_t(H_t),\hat\sigma_t^2)\), generate \(L_{t+1}=\hat\mu_t(H_t)+\hat\sigma_t Z\), where \(Z\sim N(0,1)\) and \(H_t\) denotes the conditioning history. Use a fresh random draw at each step.
Set the baseline history to the specified \(l_0\), and record the fixed regime \(\bar a\).
At step \(t\), set \(A_t=a_t\), evaluate the fitted transition distribution at the history generated so far, and draw \(L_{t+1}\).
After generating \(L_K\), evaluate the fitted outcome distribution at the complete history and draw \(Y\).
Retain the generated covariates, the regime values, and the outcome together. Repeat with new random draws for another trajectory.
Why the Algorithm Generates the Intended Distribution
The fitted joint distribution under the fixed regime is the product of its sequential conditional distributions. The following result makes explicit why drawing in temporal order reproduces that product.
Proof. At the first step, conditional on \(L_0=l_0\), the algorithm draws \(L_1\) from \(\hat f_0(\cdot\mid l_0,a_0)\). Thus the first conditional factor in the joint distribution is \(\hat f_0(l_1\mid l_0,a_0)\). Given the generated history \((l_0,l_1)\), the next draw has conditional distribution \(\hat f_1(\cdot\mid l_0,l_1,a_0,a_1)\). The chain rule therefore gives the joint distribution of \(L_1,L_2\), conditional on \(l_0\), as the product of these two factors. Repeating this argument at each time adds the factor \(\hat f_t(l_{t+1}\mid\bar l_t,\bar a_t)\). Given the completed covariate history, the final draw contributes \(\hat f_Y(y\mid\bar l_K,\bar a_{K-1})\). Multiplying these conditional factors gives the displayed expression. Normalization ensures that each step defines a probability distribution, and fresh randomness implements the specified conditional draw rather than imposing an additional dependence between draws. \(\square\)
The result is a statement about the fitted models: the algorithm reproduces their implied distribution. It does not show that the fitted distributions equal the true counterfactual distributions. That requires adequate model specification and, for causal interpretation, the identification assumptions.
Worked Example: Two Binary Covariate Transitions
Fix \(L_0=0\) and a two-step treatment regime. Suppose the fitted transition and outcome models under that regime are \(P(L_1=1)=0.6\), \(P(L_2=1\mid L_1=1)=0.8\), \(P(L_2=1\mid L_1=0)=0.25\), and \(P(Y=1\mid L_1,L_2)=0.1+0.2L_1+0.3L_2\). Each specified outcome probability lies between 0 and 1. Use uniform draws \(U_1=0.42\), \(U_2=0.70\), and \(U_Y=0.55\).
First, \(0.42<0.6\), so the first draw gives \(L_1=1\). Given that value, the probability of \(L_2=1\) is \(0.8\), and \(0.70<0.8\), so the second draw gives \(L_2=1\). The terminal success probability at this history is \(0.1+0.2(1)+0.3(1)=0.6\). Since \(0.55<0.6\), the outcome draw gives \(Y=1\). This simulated path is \((L_1,L_2,Y)=(1,1,1)\).
The probability of that discrete path is the product \(0.6(0.8)(0.6)=0.288\). As a check on the model’s implications, \(P(L_1=1,L_2=1)=0.6(0.8)=0.48\), and \(P(L_1=0,L_2=1)=0.4(0.25)=0.10\), so \(P(L_2=1)=0.58\). Hence \(E[Y]=0.1+0.2E[L_1]+0.3E[L_2]=0.1+0.2(0.6)+0.3(0.58)=0.394\). The single simulated outcome is not this expectation; it is one draw from the fitted regime-specific distribution.
Worked Example: A Normal Covariate and a Binary Outcome
Fix \(L_0=0\) and the treatment regime at \(A_0=1\). Suppose the fitted transition model is \(L_1\mid L_0=0,A_0=1\sim N(1.2,0.25)\), so its standard deviation is \(0.5\), and the outcome model is \(P(Y=1\mid L_1,A_0=1)=\operatorname{expit}(-0.8+0.5L_1)\), where \(\operatorname{expit}(x)=1/(1+e^{-x})\).
For a reproducible draw, take \(Z=1\) in the normal simulation and \(U_Y=0.4\) for the outcome. Then \(L_1=1.2+0.5(1)=1.7\). The linear predictor for the outcome is \(-0.8+0.5(1.7)=0.05\), so \(P(Y=1\mid L_1=1.7)=1/(1+e^{-0.05})\approx0.5124974\). Because \(0.4<0.5124974\), the simulated outcome is \(Y=1\).
For a density check, the normal density at \(1.7\) is \[ \frac{1}{0.5\sqrt{2\pi}}\exp\!\left(-\frac{(1.7-1.2)^2}{2(0.25)}\right) = \frac{1}{0.5\sqrt{2\pi}}e^{-0.5} \approx 0.4839. \] The joint density in \(L_1\) at \(1.7\), together with the probability mass for \(Y=1\), is approximately \(0.4839415(0.5124974)=0.2480188\), or \(0.2480\) to four decimal places. This is a density-times-mass, not the probability of observing the exact continuous value \(L_1=1.7\), which is zero.
Worked Example: Updating a Mixed History at Each Step
Under a fixed two-step regime, suppose \(P(L_1=1)=0.4\), \(L_2\mid L_1\sim N(2+L_1,1)\), and \(P(Y=1\mid L_1,L_2)=\operatorname{expit}(-2+0.5L_1+0.4L_2)\). Take \(U_1=0.30\), \(Z_2=-0.5\), and \(U_Y=0.20\).
Since \(0.30<0.4\), the first draw gives \(L_1=1\). The conditional mean for the next covariate must now be evaluated at this generated value: it is \(2+1=3\). Thus \(L_2=3+1(-0.5)=2.5\). The outcome linear predictor is \(-2+0.5(1)+0.4(2.5)=-0.5\), giving success probability \(\operatorname{expit}(-0.5)\approx0.3775407\). Since \(0.20<0.3775407\), the outcome is \(Y=1\).
The mixed discrete-continuous density-times-mass for this path is \(0.4\,\phi(2.5;3,1)\,0.3775407\), where \(\phi(\cdot;3,1)\) is the \(N(3,1)\) density. Since \(\phi(2.5;3,1)=(2\pi)^{-1/2}e^{-0.125}\approx0.3520653\), the product is approximately \(0.0531676\). The factor \(0.4\) is the probability of \(L_1=1\); the normal density factor is not a probability at a single point.
Implementation Checks and Common Failure Modes
A useful implementation audit follows the time order of the model. At each transition, confirm that the conditioning history contains exactly the variables available by that time, that the regime value is correctly substituted for treatment, and that the draw comes from the full fitted conditional distribution. Then verify that each stored trajectory has one coherent sequence of values and that the terminal outcome model uses the history actually generated.
- Using observed treatment values instead of the fixed regime. This simulates observed-care histories, not histories under the specified intervention. Set every \(A_t\) to \(a_t\) in both the model evaluation and the stored path.
- Using a later covariate before it is generated. A transition model at time \(t\) is evaluated using the history available at time \(t\), not future covariates or outcomes.
- Replacing random draws with conditional means. A conditional mean is generally not a representative draw from a conditional distribution. For nonlinear outcome models, even applying the outcome model to mean covariates need not give the average outcome probability. The binary example showed why the entire distribution of histories matters.
- Reusing the same random number unintentionally. Each step needs the randomness required by its conditional draw. Reuse can induce dependence not specified by the fitted kernels.
- Ignoring invalid predictions or unstable histories. Check that predicted probabilities are in \([0,1]\), categorical probabilities sum to one, and continuous-model parameters are valid. Also inspect whether simulated covariate values or histories are extreme relative to the data used to fit the models.
Compared with sequential integration, forward simulation replaces integration over histories with draws from the successive conditional distributions. When all relevant distributions are discrete and the history space is small, enumeration or exact integration can provide a direct check without random simulation. Forward simulation is convenient for complex histories but reflects model misspecification, extrapolation, and simulation variability; it cannot repair unsupported histories or incorrect conditional models. It also does not itself produce a population-average target: averaging over baseline covariates is a separate operation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Under a two-time regime, should \(A_1\) be sampled from a fitted treatment model or set to the prescribed regime value when generating \(L_2\)?
- A binary covariate has fitted probability \(P(L_1=1)=0.35\). With uniform draw \(U=0.20\), what value is generated?
- If \(L_2\mid L_1=1\sim N(4,9)\), state the standard deviation and write a normal-draw expression using \(Z\sim N(0,1)\).
- Why is the density at one continuous covariate value not the probability of that exact value?
- Give one reason a simulated distribution may fail to represent the true counterfactual distribution even when the simulation code correctly follows the algorithm.
Solutions. (1) Set \(A_1\) to its prescribed regime value; the algorithm intervenes on treatment rather than drawing it. (2) Since \(0.20<0.35\), generate \(L_1=1\). (3) The standard deviation is \(\sqrt{9}=3\), so generate \(L_2=4+3Z\). (4) For a continuous distribution, a single point has probability zero; the density describes probability relative to a small interval and its width. (5) For example, a transition or outcome model may be misspecified, or the regime may require extrapolation into histories poorly represented in the fitted data.