Continuous time markov chains, is this step by step example correct

Question

I have some questions regarding CTMC... and most importantly whether the step-by-step example I provide below is correct. My main sources about CTMC are: ([1], and [2]).

Let's assume 3 possible states $S = {1, 2, 3}$. There are the following data instances $D = {d_1, d_2, d_3, d_4}$ from which I want to model as a CTMC. In parenthesis the time spent at the state before transiting to the next.

$ d_1 = 1(3min) \rightarrow 2(3min) \rightarrow 3\\ d_2 = 1(8min) \rightarrow 2(12min) \rightarrow 3\\ d_3 = 1(2min) \rightarrow 3\\ d_4 = 2(2min) \rightarrow 3 $

The main question(s) I want to ask and answer from the CTMC is:

• Problem: What is the probability of transiting from state $i$ to state $j$ given that $s$ time has elapsed since entering state $i$?

So to form the CTMC I start with the embedded transition matrix for $D$ which is:

$ p_{ij} = \begin{bmatrix} 0 & 2/3 & 1/3 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix} $

Now the next step is to calculate the rates. From [1, p3] it says that "Each time state i is visited, the chain spends, on average, $E(Hi) = 1/a_i$ units of time there before moving on". From what I understand $a_i$ is the set of rates for state $i$, also denoted in other texts as $\lambda_i$ or $v_i$, for example an operator receives 2 calls in an hour, hence $a_i = 1$ call in half an hour. However, I am not sure how to calculate this from the above data $D$. But, I can calculate the average time spent at each state and then from the above equation, solving for $a_i = 1 / E(H_i)$ I can obtain these rates...

$E(H_i) = \begin{bmatrix}4.33 & 5.66 & 0 \end{bmatrix}$, denoting with 0 the terminal (consume) state.

Hence,

$a_i = \begin{bmatrix}0.23 & 0.18 & 0 \end{bmatrix}$, again 0 is for consume state.

• Q1: Is this valid?
• Q2: Given that in this way I assume that the probability distribution of time spent is an exponential distribution, I lose the true variance of the actual distribution as the variance of the exponential distribution is the same as the mean... Is this OK? Is there any other way?
• Q3: How can terminal/consume states be represented?

Again from [1, p3] I can now compute the transition rate matrix (infinitesimal generator) $Q$:

$Q_{ij} = -a_i$, if $i==j$

$Q_{ij} = a_i * p_{ij}$, if $i != j$

$ Q = \begin{bmatrix} -0.23 & 0.15 & 0.08\\ 0 & -0.18 & 0.18\\ 0 & 0 & 0 \end{bmatrix} $

• Q4: Although [1, p3] states that $Q_{ii} = -a_i$, [2, p227] says that $Q_{ii} = 0$. Which one is correct then?

Now that we know $Q$, from [1, p2] I can calculate $P_{ij}(t)$ as follows:

$P_{ij}(t) = e^{Qt}$ (from what I read there are a number of ways to compute that and I can use a numerical software that implements them)

• Q5: I assume $t = s$, ie the elapsed time since entering state $i$
• Q6: Why do I need to calculate the complete matrix $P_{ij}(t)$ if I am interested in one particular ${ij}$. Is there any way to calculate that.
• Q7: I read in [2, p.228] that if the process can be modelled as a Poisson process then the solution is: $P_{ij}(t) = \dfrac{(\lambda_i t)^k}{k!}e^{-\lambda_i t}$, where $k$ is the number of events in an interval of length t, which in my case above is $1$. Can I assume that the above process is a Poisson process?
• Q8: And lastly is everything I have explained in the above step-by-step example correct?

Thank you the most for your time reading this and if you can provide any answers to my questions.

Answer 1

I believe the best strategy for a problem of this kind would be to proceed in two steps:

1. Fit a continuous time Markov chain model to the data by estimating the (infinitesimal) generator $Q$.

2. Using the estimated generator and the Kolmogorov backward equations, find the probability that a Markov chain following the fitted model transitions from state $i$ to state $j$ in time $s$.

The generator can be estimated directly, no need to first go via the embedded Markov chain. A summary of methods (considering the more complicated case of missing data) can for example be found in Metzner et al. (2007). While estimating the generator is possible using the observations you list in your example, you have very little data available: Your data contains 6 observed transitions (the usable data being the time until the transition and the pair of states). From this data you need to estimate the six transition rates which make up the off-diagonal elements of the generator. Since the amount of data is not much larger than the number of parameters to be estimated, the resulting estimate may not be very reliable. In particular, with so little data it will be difficult to know whether unobserved transitions are impossible, or just have not been observed by chance.

The general form for the generator of a continuous time Markov chain with three states is $$ Q = \begin{pmatrix} -a_1 & a_{12} & a_{13} \\ a_{21} & -a_2 & a_{23} \\ a_{31} & a_{32} & -a_3 \end{pmatrix} $$ where the diagonal elements satisfy $a_i = \sum_{j\neq i}a_{ij}$. The easiest way for estimating $Q$ may be using the maximum likelihood estimate (MLE): from equation (13) in Metzner et al. (2007) we know that the MLE for $a_{ij}$ is given by $$ \hat a_{ij} = \frac{N_{ij}}{R_i} $$ where $N_{ij}$ is the number of observed transitions from $i$ to $j$ and $R_i$ is the total time spent in state $i$ by the observed process. With your data the estimates are $\hat a_{12} = 2 / (3+8+2) = 0.154$, $\hat a_{13} = 1 / (3+8+2) = 0.077$, $\hat a_{21} = 0 / (2+3+12) = 0$, $\hat a_{23} = 3 / (2+3+12) = 0.176$. This seems compatible with the results you obtained. (There is a small complication here, because the cited article assumes that the Markov chain is observed until a fixed time $T$ rather then until a certain jump, but it is easy to check that this distinction does not matter for computing the MLE.)

The data you give lists no observations for the system being in state 3, so no estimates of $a_{31}$ and $a_{32}$ are possible from the given data. If the system has been observed to stay in state 3 for a time without making a transition, then $R_3$ is positive and we would get the MLEs $\hat a_{31} = \hat a_{32} = 0$. In this case, state 3 would be absorbing (represented by a row of zeros in the generator) and the generator would be the one given by you above.

Finally, using the Kolmogorov backward equations, it is possible to compute the transition probabilities for the fitted model. From equation (14) in your first reference we know that the transition matrix over a time interval of length $s$ is given by $$ P(s) = \exp(sQ) $$ where $sQ$ is obtained by multiplying every element in $Q$ by $s$ and $\exp(sQ)$ is the matrix exponential. The matrix exponential can be directly computed in many software packages (e.g. using expm in matlab), or "by hand" by diagonalising the matrix $sQ$. The required probability that a Markov chain following the fitted model transitions from state $i$ to state $j$ in time $s$ can then be found in row $i$, column $j$ of the resulting matrix.

answers to your questions:

• problem: see above

• Q1: I believe your approach is valid but unnecessarily complicated

• Q2: You don't need to "assume" that the waiting times are exponential, this is automatically true in CTMC models. Wikipedia has a section about how to estimate the parameter of an exponential distribution. (And I believe your approach coincides with the MLE, so should be valid.)

• Q3: terminal/consume states are represented as a row of zeros in the generator.

• Q4: the diagonal elements of the generator are always equal to the negative of the sum of all non-diagonal elements in the same row, i.e. they are equal to $-a_i$. This may be 0 if the whole row is 0.

• Q5: this seems not to be a question

• Q6: I don't think there is a simple way to compute just one element of a matrix exponential. The usual thing to do is to compute the whole matrix and then to pick out the elements you need.

• Q7: whether a Poisson process is a good model for your data depends mostly on the context and origin of your data. Not though that Poisson processes can take arbitrarily large values, so if you know that you have three states a Poisson process is not a valid model to describe your data.

• Q8: I don't know.