Discrete Event-Time Computation Of The Product Integral

From an Event-Time Table to a Transition Matrix

A product integral can be computed as a finite sequence of matrix updates whenever the estimated transition intensities change only at observed event times. The practical task is to turn a sorted event-time table into those updates, apply them in the correct order, and retain the estimate at each requested reporting time. Doing this explicitly makes the calculation reproducible and exposes errors that can be hidden by a software command.

Use a fixed ordering of the $K$ states. Transition matrices use rows for origin states and columns for destination states. Thus, the $(h,j)$ entry of a transition-probability matrix is the probability of being in state $j$ at the later time, conditional on being in state $h$ at the earlier time. The matrix is updated on the right as event times move forward.

Definition: At an observed transition time $u$, let $Y_h(u)$ be the number at risk in state $h$ immediately before $u$, and let $dN_{hj}(u)$ be the number of observed transitions from $h$ to $j$ at $u$, for $h\ne j$. Define the estimated increment $\widehat{\Delta A}(u)$ by $$ \widehat{\Delta A}_{hj}(u)=\frac{dN_{hj}(u)}{Y_h(u)}\quad (h\ne j),\qquad \widehat{\Delta A}_{hh}(u)=-\sum_{j\ne h}\widehat{\Delta A}_{hj}(u). $$ The event-time update matrix is $B(u)=I+\widehat{\Delta A}(u)$, where $I$ is the $K\times K$ identity matrix. If $Y_h(u)=0$, there can be no observed transition out of $h$; that row of $\widehat{\Delta A}(u)$ is defined to be zero.

The event-time table contains one row for each distinct time at which at least one transition occurs. Its transition counts include all transitions observed at that time, and its risk counts are measured just before that time. Censoring does not itself create a transition increment, although it can change who remains at risk at subsequent times. If several transition types occur at the same time, their counts belong in the same row and the corresponding updates are assembled into one matrix.

Recursion: Start with $\widehat P(0,0)=I$. For sorted transition times $u_1<\cdots<u_m$, update $$ \widehat P(0,u_r)=\widehat P(0,u_{r-1})B(u_r),\qquad r=1,\ldots,m, $$ where the initial value before the first event is $I$. Between transition times, the estimate is constant. The estimate at a requested time $t$ includes the update at $t$, if there is one; equivalently, it is the product of all updates with event times $u\le t$.

Why the Recursion Works

The product-integral recursion is the discrete event-time form of accumulating transition changes over time. At a single event time, the estimated change in the transition matrix is represented by adding the intensity increment to the identity. After that update, the current transition matrix already accounts for all earlier event times. The next update therefore composes the new change with the current estimate on its right.

The order follows from matrix multiplication for row-origin transition matrices. If a row vector describes the distribution over states at time $s$, multiplying by a transition matrix advances that distribution. Advancing first to an earlier time and then through a later increment gives $(\mathbf{v}B(u_1))B(u_2)$, not $(\mathbf{v}B(u_2))B(u_1)$. In general, these two products differ.

A useful check follows directly from the construction. Each row of $\widehat{\Delta A}(u)$ sums to zero, so each row of $B(u)$ sums to one. If the total number of departures from state $h$ at $u$ does not exceed $Y_h(u)$, then $$ B_{hh}(u)=1-\sum_{j\ne h}\frac{dN_{hj}(u)}{Y_h(u)}\ge 0, $$ and its off-diagonal entries are nonnegative. Thus $B(u)$ is a row-stochastic matrix. A product of row-stochastic matrices is also row-stochastic: its entries are nonnegative, and each row sum remains one because $B(u)\mathbf{1}=\mathbf{1}$. Consequently the recursion preserves nonnegative entries and row sums of one, apart from rounding error.

Worked Example: Two-State Event-Time Updates

Suppose a hypothetical study follows people in state A or state B. At time $u_1=1$, there are 10 people at risk in A and 2 transitions from A to B. At time $u_2=3$, there are 8 at risk in A and 5 in B, with one transition from B to A and none from A to B. Compute the estimate at times 0, 2, and 3. Assume the table lists all observed transitions through time 3.

At time 1, the estimated increment and update are

\widehat{\Delta A}(1)= \begin{pmatrix}-2/10&2/10\\0&0\end{pmatrix} = \begin{pmatrix}-0.2&0.2\\0&0\end{pmatrix}, \qquad B(1)= \begin{pmatrix}0.8&0.2\\0&1\end{pmatrix}.

Initialization gives $\widehat P(0,0)=I$, so the estimate at time 0 is the identity matrix. There are no events between times 1 and 3, so $\widehat P(0,2)=B(1)$.

At time 3, the B-to-A increment is $1/5=0.2$, and the diagonal entry for state B is $-0.2$. Hence

B(3)= \begin{pmatrix}1&0\\0.2&0.8\end{pmatrix}, \qquad \widehat P(0,3)=B(1)B(3) = \begin{pmatrix} 0.8&0.2\\ 0&1 \end{pmatrix} \begin{pmatrix} 1&0\\ 0.2&0.8 \end{pmatrix} = \begin{pmatrix}0.84&0.16\\0.2&0.8\end{pmatrix}.

For example, the probability estimate in the first row and second column is $0.8(0)+0.2(0.8)=0.16$. The row sums are $0.84+0.16=1$ and $0.2+0.8=1$. The estimates are therefore $I$ at time 0, $B(1)$ at time 2, and the displayed product at time 3.

Simultaneous Events and Chronological Storage

When distinct transition types occur at the same recorded time, use one combined update matrix for that time. This is not a claim that the events have a known within-time order. It implements the event-time increment based on all the transition counts tied at that time. Applying separate updates in an arbitrary order can give a different answer and introduces an ordering that the data do not supply.

Worked Example: Combining Tied Transitions

Consider a hypothetical three-state process with states 0, 1, and 2. At time 1, there are 10 people at risk in state 0 and 2 transitions from 0 to 1. At time 2, there are 5 at risk in state 0 and 4 in state 1; two people move from 0 to 1, none from 0 to 2, and two move from 1 to 2. There are no observed transitions out of state 2 at time 2. Find the transition matrix at time 2.

The time-1 update, with rows and columns ordered $0,1,2$, is

B(1)= \begin{pmatrix} 0.8&0.2&0\\ 0&1&0\\ 0&0&1 \end{pmatrix}.

At time 2, both observed transition types enter one increment matrix. The two transitions out of state 0 give proportions $2/5=0.4$ to state 1 and $0/5=0$ to state 2. The two transitions out of state 1 give $2/4=0.5$ to state 2. Therefore

B(2)= \begin{pmatrix} 0.6&0.4&0\\ 0&0.5&0.5\\ 0&0&1 \end{pmatrix}.

Chronological updating gives

\widehat P(0,2)=B(1)B(2)= \begin{pmatrix} 0.8&0.2&0\\ 0&1&0\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 0.6&0.4&0\\ 0&0.5&0.5\\ 0&0&1 \end{pmatrix} = \begin{pmatrix} 0.48&0.42&0.10\\ 0&0.5&0.5\\ 0&0&1 \end{pmatrix}.

The first-row, third-column entry is $0.8(0)+0.2(0.5)+0(1)=0.10$. The transition from 0 to 1 at time 1 can be followed by the time-2 transition from 1 to 2, which contributes $0.2\times0.5=0.10$. Each row sums to one. The time-2 updates must not be split into an invented sequence of transitions: the combined matrix is the update specified by the tied event-time counts.

Implementing the Table and Saving Requested Times

A reliable computation separates the event-time updates from the reporting schedule. First sort the distinct transition times and compute each $B(u)$. Then move through those times once, keeping a running matrix. For a requested time between two event times, save the current running matrix without applying an update. At a requested time that is itself an event time, apply the update first and save the resulting matrix. This convention reports $\widehat P(0,t)$, including events observed at $t$; a left-limit estimate would instead exclude the update at $t$.

Initialize.
Set the running matrix to $I$. Record this value for any requested time before the first transition event.

Prepare each event row.
For every distinct transition time, use the risk counts immediately before that time and all observed transition counts at that time to construct $\widehat{\Delta A}(u)$ and $B(u)$.

Update in order.
For each time in ascending order, replace the running matrix by its product with $B(u)$ on the right. Do not reorder the states or event times partway through the calculation.

Save requested estimates.
After processing all events with time at most $t$, save the running matrix as the estimate at $t$. If no transition occurs between two requested times, their saved matrices are equal.

Worked Example: Reporting Between Event Times

In a hypothetical two-state process, the only observed transition times are 2 and 5. Their update matrices, already computed from valid risk sets and counts, are

B(2)=\begin{pmatrix}0.9&0.1\\0&1\end{pmatrix}, \qquad B(5)=\begin{pmatrix}1&0\\0.25&0.75\end{pmatrix}.

The requested times are 0, 3, 5, and 6. Initialization gives $\widehat P(0,0)=I$. At time 3, only the time-2 update has occurred, so $\widehat P(0,3)=B(2)$. At time 5, include the update at 5:

\widehat P(0,5)=B(2)B(5) = \begin{pmatrix}0.9&0.1\\0&1\end{pmatrix} \begin{pmatrix}1&0\\0.25&0.75\end{pmatrix} = \begin{pmatrix}0.925&0.075\\0.25&0.75\end{pmatrix}.

There are no further transition events by time 6, so $\widehat P(0,6)=\widehat P(0,5)$. The saved estimates are therefore $I$ at 0, $B(2)$ at 3, and the displayed product at both 5 and 6. The estimate just before time 5 would be $B(2)$, not the product including $B(5)$.

Checks, Alternatives, and Limitations

Before interpreting a computed matrix, check the input table and the numerical result. For each event time and origin state, verify that the risk count is measured just before the event and that the total number of observed departures is no greater than the number at risk. Check that transition counts are assigned to the intended origin and destination states. After every update, inspect row sums and entries; substantial negative values or row sums far from one usually indicate a coding, indexing, or arithmetic error. Tiny discrepancies can arise from rounding, so retain full precision during multiplication and round only for presentation.

A common mistake is to apply an update at time $u$ by multiplying it on the left, or to process event times in reverse order. The correction is to keep origin states in rows, destination states in columns, and update by $\widehat P_{\mathrm{current}}B(u)$ as times increase. Another mistake is to calculate a separate matrix for each transition type tied at $u$ and multiply those matrices in an arbitrary order. Use one matrix built from all counts at $u$.

This event-time computation is distinct from replacing the intensity by a constant matrix over an interval and using a matrix exponential. The exponential approach describes evolution under a constant-intensity model; the product-integral recursion uses the estimated increments at the observed event times. They need not give the same numerical result. The recursion itself is straightforward, but the population interpretation still depends on appropriate data collection, risk-set construction, and assumptions behind estimating transition intensities, including suitable handling of censoring. It does not correct a misspecified state structure or invalid risk counts.

The method stores a matrix for any requested time, but the matrix represents transitions from the chosen starting time and uses the state ordering fixed throughout the computation. Changing the starting time, state ordering, or convention for including an event at the reporting time changes what must be computed or how the result is read. Record those choices alongside the output.

Check Your Understanding

Work through each question, then compare your reasoning with the solutions.

At one event time, $Y_1=8$, with one transition from state 1 to state 2 and two from state 1 to state 3. What are the three entries in row 1 of the update matrix?
There are transition events at times 2 and 4, and the requested reporting time is 3. Which updates are included?
Why should transitions of two different types tied at the same event time be combined into one update?
If the event-time update matrices are $B(1)$ and $B(2)$, which product gives the chronological update for row-origin transition matrices?

Solutions. (1) The departure proportion is $(1+2)/8=0.375$, so the diagonal entry is $1-0.375=0.625$. The two off-diagonal entries are $1/8=0.125$ and $2/8=0.25$; the row is $(0.625,0.125,0.25)$. (2) Include only the update at time 2; the estimate is constant from time 2 up to, but not including, the update at time 4. (3) The tied counts define one increment at that time; splitting them and multiplying in an arbitrary order would impose an unsupported within-time ordering. (4) Use $B(1)B(2)$, because the earlier update acts before the later one.