Simulating Paths Under Non-Markov Transition Rules

Why Simulation Needs More Than the Current State

In a Markov multi-state model, the current state contains all the information needed to describe future transitions. In a non-Markov model, future transition rates can also depend on how long the process has been in the current state, the time on the study clock, or earlier transitions. A simulation must therefore carry forward the relevant history, not just a state label. This matters when predicting complete paths, evaluating a proposed model, or generating data for a design study.

The basic construction is sequential. At the current time, use the current state and history to draw the time until the next transition. At that event time, draw which transition occurs. Record the transition and its time, update the history, and repeat. The same procedure works for duration-dependent hazards and for hazards that use other recorded features of the path.

Definition: Let $X(t)$ denote the state occupied at time $t$, and let $H_t$ denote the observed history immediately before $t$, including any information the model uses from the path so far. If the process is in state $r$, write $\lambda_{rs}(t\mid H_t)$ for the conditional transition intensity from $r$ to destination $s$. The total intensity of leaving $r$ is $$ \lambda_{r\cdot}(t\mid H_t)=\sum_{s\ne r}\lambda_{rs}(t\mid H_t). $$

An intensity has units of inverse time. Informally, conditional on the history and on still occupying $r$, the probability of a transition to $s$ during a short interval of length $dt$ is approximately $\lambda_{rs}(t\mid H_t)dt$. The history may include the current sojourn duration, the entry time into the state, previous states, or cumulative exposure, depending on the specified model.

For simulation, assume that the transition intensities are nonnegative, locally integrable, and specified conditional on the history that has been generated so far. The process starts from a specified state and initial history. Transitions occur one at a time, with no simultaneous jumps. Any external time-varying covariates must either be known over the simulated interval or be generated by an explicitly specified rule. These conditions make the next-event distribution well defined. If the model leaves some required history or covariate unspecified, the simulation is not fully defined.

Drawing the Next Waiting Time and Destination

Suppose the process has just entered state $r$ at time $t_0$, with its history at entry recorded. Consider the path that would continue to occupy $r$ without another transition. Along that no-transition path, the model specifies the total intensity at each future time $u$. The cumulative total hazard over the next $w$ time units is the integral of that intensity. Its exponential gives the probability of no transition during that interval.

\Lambda(w)=\int_{t_0}^{t_0+w}\lambda_{r\cdot}(u\mid H_u^{(0)})\,du, \qquad P(W>w\mid H_{t_0})=\exp\{-\Lambda(w)\}.

Here $H_u^{(0)}$ is the history that applies at $u$ if no transition has yet occurred after $t_0$. Draw $U_1$ uniformly on $(0,1)$, and set $E=-\log U_1$. The next waiting time $W$ is obtained by solving $\Lambda(W)=E$. If the integrated hazard never reaches $E$, there is no further transition under this draw; the path remains in its current state for the rest of the simulation horizon.

W=\inf\{w\geq 0:\Lambda(w)\geq -\log U_1\}.

If a transition occurs at time $t_*=t_0+W$, draw its destination using a second uniform random number. For a destination $s$, the conditional probability of that transition at the event time is its intensity divided by the total intensity:

P\{X(t_*)=s\mid \text{a transition at }t_*,H_{t_*}\} = \frac{\lambda_{rs}(t_*\mid H_{t_*})} {\lambda_{r\cdot}(t_*\mid H_{t_*})}.

This draw requires a positive total intensity at the event time. If the total intensity is zero throughout an interval, no transition can occur in that interval. After selecting a destination, append the transition and its time to the history. The next waiting-time calculation must use the updated history. The waiting-time draw is not a draw from a single fixed distribution for the entire path: its distribution can change at every transition.

Worked Examples

Worked Example: Constant Competing Transition Rates

Hypothetical setup. A process has just entered state A at time 0. Until it leaves A, the rates are constant: A to B has intensity $0.20$ per day and A to C has intensity $0.10$ per day. Use uniforms $U_1=0.50$ and $U_2=0.80$, with B ordered before C in the destination draw.

The total intensity is $0.20+0.10=0.30$ per day. The waiting-time threshold is $E=-\log(0.50)=0.6931$. Since the cumulative total hazard after $w$ days is $0.30w$, solve $0.30w=0.6931$, giving $w=2.310$ days.

At that event time, the destination probabilities are $P(B)=0.20/0.30=0.6667$ and $P(C)=0.10/0.30=0.3333$. With B first, $U_2=0.80$ falls in the C interval $(0.6667,1)$. The simulated path therefore leaves A for C after 2.310 days. The rates have no duration or history variation in this example, so the integrated-hazard inversion reduces to division by a constant.

Worked Example: A Piecewise-Constant Duration Hazard

Hypothetical setup. A process enters A at time 0. The A-to-B intensity is $0.10$ per day during the first 4 days in A, then $0.30$ per day. The A-to-C intensity remains $0.05$ per day. Use $U_1=0.25$, followed, if an event occurs, by $U_2=0.70$; order B before C.

The total intensity is $0.15$ per day for the first 4 days and $0.35$ per day thereafter. The cumulative hazard by day 4 is $4(0.15)=0.60$. The threshold is $-\log(0.25)=1.3863$, which exceeds 0.60, so the transition does not occur in the first interval. The remaining hazard threshold is $1.3863-0.60=0.7863$. At $0.35$ per day, this takes $0.7863/0.35=2.247$ additional days. Thus the waiting time is $4+2.247=6.247$ days.

At day 6.247, the destination probabilities are $P(B)=0.30/0.35=0.8571$ and $P(C)=0.05/0.35=0.1429$. Since $U_2=0.70<0.8571$, the destination is B. The calculation checks the change point explicitly: using the later rate over all 6.247 days would give the wrong waiting time.

Worked Example: Updating the History After a Transition

Hypothetical setup. Suppose the process enters B from A. The model uses the immediately preceding state: after A-to-B, the intensities are B-to-D $0.40$ per day and B-to-E $0.10$ per day. Use $U_1=0.60$, then $U_2=0.75$, with D ordered before E.

The total intensity in B under this history is $0.50$ per day. The threshold is $-\log(0.60)=0.5108$, so the waiting time is $0.5108/0.50=1.022$ days. At the event time, the probability of D is $0.40/0.50=0.80$, with the remaining probability $0.20$ assigned to E. Because $0.75<0.80$, the simulated next transition is B-to-D.

The prior state is part of the history used to define the B transition rates. If the process had instead entered B from C, the stipulated rates would be B-to-D $0.10$ and B-to-E $0.30$ per day. The total rate would then be $0.40$ per day, and the destination probabilities would be $0.25$ for D and $0.75$ for E. This contrast shows why a simulator must retain the relevant prior transition rather than reset its input to “currently in B” alone.

Why the Construction Works

The waiting-time formula follows from the definition of the total intensity. Conditional on no transition through time $u$, let $S(u)$ be the probability of remaining in the current state through $u$. Over a short interval $du$, the probability of leaving is approximately $\lambda_{r\cdot}(u\mid H_u^{(0)})du$. Therefore the survival probability satisfies

\frac{dS(u)}{du} =-\lambda_{r\cdot}(u\mid H_u^{(0)})S(u), \qquad S(t_0)=1.

Dividing by $S(u)$, integrating from $t_0$ to $t_0+w$, and exponentiating gives $S(t_0+w)=\exp\{-\Lambda(w)\}$. Inverse-transform sampling with $E=-\log U_1$, which has a unit exponential distribution, then gives the stated waiting-time draw.

For the destination draw, the cause-specific density of a transition to $s$ at time $t_0+w$ is the probability of surviving without a transition until that time multiplied by $\lambda_{rs}(t_0+w\mid H_{t_0+w})$. Summing these densities over destinations replaces the transition-specific intensity by $\lambda_{r\cdot}$. Their ratio at a given event time is therefore $\lambda_{rs}/\lambda_{r\cdot}$, as used in the destination draw. This argument assumes the destination intensities describe all possible ways to leave the current state.

Initialize.
Set the initial time, state, and all history features required by the model.

Draw the next event time.
Compute the total intensity along the no-transition path, integrate it, and invert against $-\log U_1$. Stop at the simulation horizon if the threshold is not reached before then.

Draw the destination.
At the event time, divide each destination intensity by the total intensity and use $U_2$ to select a destination.

Update and repeat.
Record the event time and transition, update the history, and return to Step 2 until the horizon or an absorbing state is reached.

Verification, Alternatives, and Failure Modes

A useful verification starts with a case whose probabilities can be computed exactly. Suppose A has constant intensities $0.20$ per day to B and $0.10$ per day to C, and the simulation horizon is 5 days. The total rate is $0.30$ per day, so the probability of no transition by day 5 is $\exp(-0.30\times5)=0.2231$. The probability of a transition by day 5 is $1-0.2231=0.7769$. Splitting that probability by the destination rates gives

P(B\text{ by day }5)=\frac{0.20}{0.30}(0.7769)=0.5179,\qquad P(C\text{ by day }5)=\frac{0.10}{0.30}(0.7769)=0.2590.

Worked Example: Checking Simulated Frequencies Against Exact Probabilities

Hypothetical check. Simulate 1,000 independent paths under the constant-rate setup above, stopping each path at day 5. Suppose the outcomes are 518 paths in B, 260 in C, and 222 still in A.

The exact expected outcome counts are $1000(0.5179)=517.9$ in B, $1000(0.2590)=259.0$ in C, and $1000(0.2231)=223.1$ still in A. The simulated proportions are $0.518$, $0.260$, and $0.222$. They are close to the exact probabilities, and both the simulated counts and exact expected counts sum to 1,000. Small differences are expected from random variation; agreement in one run is a check, not proof that the code is correct.

A stronger check repeats the comparison over many runs and tests edge cases, such as a zero intensity for one destination or a hazard change exactly at a known boundary. For the piecewise-constant example, compare simulated waiting-time and destination frequencies with probabilities calculated from the integrated hazards in each segment.

For constant rates, the waiting-time inversion is especially simple. For piecewise-constant rates, integrate separately over each segment, as in the second example. More general time-varying rates may require numerical integration and a numerical root solver, or a thinning algorithm with a valid upper bound on the total intensity. A time-grid method is another approximation: it checks for transitions at discrete intervals. It can miss multiple events in one interval and introduces time-step error, so reducing the step size and checking convergence are important if it is used.

A common mistake is to draw a waiting time using only the current state’s total rate, then choose the destination using rates from an outdated history. The correction is to evaluate both the total and destination-specific intensities using the same current history and event time, then update history only after the transition. Another mistake is to use a constant-rate exponential draw when the total intensity changes with duration or calendar time; the correct draw solves the integrated-hazard equation. Diagnostics include checking that destination probabilities are nonnegative and sum to one, that simulated paths obey allowed transitions, and that event times increase and do not exceed the horizon.

Simulation is only as reliable as the model specification and numerical implementation. A finite simulation provides approximate summaries with Monte Carlo error, not exact probabilities. Rare transitions may require many paths to estimate well. Numerical integration or root-finding tolerances can also affect results, especially near sharp hazard changes. The Markov special case is simpler because the current state alone determines transition rates; it can be used as a benchmark when the history dependence is removed. In contrast, replacing a history-dependent model with a Markov model is not a neutral computational shortcut: it changes the model unless the state representation truly contains the relevant predictive information.

Key takeaway: Simulate a non-Markov path by integrating the total intensity conditional on the current history, drawing the next event time, selecting its destination from the transition-specific intensities at that time, and updating the history before continuing.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In a state, two constant outgoing intensities are 0.12 and 0.08 per hour. What is the total rate, and what is the probability that the first destination is the one with rate 0.12?
A total hazard is 0.25 per day for the first 2 days, then 0.50 per day. If $U_1=0.50$, does the event occur in the first segment? If not, find the waiting time.
Why must the destination probabilities be computed at the simulated event time rather than automatically at the time of entry into the state?
In 500 simulated paths, the exact probability of no event by a fixed horizon is 0.20. What count of paths with no event is expected, and would a count of 90 by itself prove a coding error?

Solutions. (1) The total rate is $0.12+0.08=0.20$ per hour. The probability of the first destination is $0.12/0.20=0.60$. (2) The threshold is $-\log(0.50)=0.6931$. The first segment contributes $2(0.25)=0.50$, so no event occurs in that segment. The remaining threshold is $0.1931$, which takes $0.1931/0.50=0.3863$ days in the second segment. The waiting time is $2.386$ days. (3) Intensities, and therefore their relative shares, can vary with time and with the evolving history. The conditional destination distribution at an event depends on those intensities at the event time. (4) The expected count is $500(0.20)=100$. A count of 90 does not by itself prove an error: random variation can produce a difference. Check repeated-run behavior and implementation details against the exact probabilities.