General Definitions

A DTMC is a sequence of random variables \(X_{1},\: X_{2}\: ,\ldots,\:X_{n},\ldots\) characterized by the Markov property (also known as memoryless property, see Equation \(\ref{eq:markovProp}\)). The Markov property states that the distribution of the forthcoming state \(X_{n+1}\) depends only on the current state \(X_{n}\) and doesn’t depend on the previous ones \(X_{n-1},\: X_{n-2},\ldots,\: X_{1}\).

\[\begin{equation} Pr\left(X_{n+1}=x_{n+1}\left|X_{1}=x_{1},X_{2}=x_{2},\ldots,X_{n}=x_{n}\right.\right)=Pr\left(X_{n+1}=x_{n+1}\left|X_{n}=x_{n}\right.\right). \label{eq:markovProp} \end{equation}\]

The set of possible states \(S=\left\{ s_{1},s_{2},...,s_{r}\right\}\) of \(X_{n}\) can be finite or countable and it is named the state space of the chain.

The chain moves from one state to another (this change is named either ‘transition’ or ‘step’) and the probability \(p_{ij}\) to move from state \(s_{i}\) to state \(s_{j}\) in one step is named transition probability:

\[\begin{equation} p_{ij}=Pr\left(X_{1}=s_{j}\left|X_{0}=s_{i}\right.\right). \label{eq:trProp} \end{equation}\]

The probability of moving from state \(i\) to \(j\) in \(n\) steps is denoted by \[p_{ij}^{(n)}=Pr\left(X_{n}=s_{j}\left|X_{0}=s_{i}\right.\right).\]

A DTMC is called time-homogeneous if the property shown in Equation \(\ref{eq:mcHom}\) holds. Time homogeneity implies no change in the underlying transition probabilities as time goes on. \[\begin{equation} Pr\left(X_{n+1}=s_{j}\left|X_{n}=s_{i}\right.\right)=Pr\left(X_{n}=s_{j}\left|X_{n-1}=s_{i}\right.\right). \label{eq:mcHom} \end{equation}\]

If the Markov chain is time-homogeneous, then \(p_{ij}=Pr\left(X_{k+1}=s_{j}\left|X_{k}=s_{i}\right.\right)\) and \(p_{ij}^{(n)}=Pr\left(X_{n+k}=s_{j}\left|X_{k}=s_{i}\right.\right)\), where \(k>0\).

The probability distribution of transitions from one state to another can be represented into a transition matrix \(P=(p_{ij})_{i,j}\), where each element of position \((i,j)\) represents the transition probability \(p_{ij}\). E.g., if \(r=3\) the transition matrix \(P\) is shown in Equation \(\ref{eq:trPropEx}\)

\[\begin{equation} P=\left[\begin{array}{ccc} p_{11} & p_{12} & p_{13}\\ p_{21} & p_{22} & p_{23}\\ p_{31} & p_{32} & p_{33} \end{array}\right]. \label{eq:trPropEx} \end{equation}\]

The distribution over the states can be written in the form of a stochastic row vector \(x\) (the term stochastic means that \(\sum_{i}x_{i}=1, x_{i} \geq 0\)): e.g., if the current state of \(x\) is \(s_{2}\), \(x=\left(0\:1\:0\right)\). As a consequence, the relation between \(x^{(1)}\) and \(x^{(0)}\) is \(x^{(1)}=x^{(0)}P\) and, recursively, we get \(x^{(2)}=x^{(0)}P^{2}\) and \(x^{(n)}=x^{(0)}P^{n},\, n>0\).

DTMC are explained in most theory books on stochastic processes, see and for example. Valuable references online available are: , and .

Properties and classification of states

A state \(s_{j}\) is said accessible from state \(s_{i}\) (written \(s_{i}\rightarrow s_{j}\)) if a system starting in state \(s_{i}\) has a positive probability to reach the state \(s_{j}\) at a certain point, i.e., \(\exists n>0:\: p_{ij}^{(n)}>0\). If both \(s_{i}\rightarrow s_{j}\) and \(s_{j}\rightarrow s_{i}\), then \(s_{i}\) and \(s_{j}\) are said to communicate.

A communicating class is defined to be a set of states that communicate. A DTMC can be composed of one or more communicating classes. If the DTMC is composed of only one communicating class (i.e., if all states in the chain communicate), then it is said to be irreducible. A communicating class is said to be closed if no states outside of the class can be reached from any state inside it.

If \(p_{ii}=1\), \(s_{i}\) is defined as absorbing state: an absorbing state corresponds to a closed communicating class composed by one state only.

The canonical form of a DTMC transition matrix is a matrix having a block form, where the closed communicating classes are shown at the beginning of the diagonal matrix.

A state \(s_{i}\) has period \(k_{i}\) if any return to state \(s_{i}\) must occur in multiples of \(k_{i}\) steps, that is \(k_{i}=gcd\left\{ n:Pr\left(X_{n}=s_{i}\left|X_{0}=s_{i}\right.\right)>0\right\}\), where \(gcd\) is the greatest common divisor. If \(k_{i}=1\) the state \(s_{i}\) is said to be aperiodic, else if \(k_{i}>1\) the state \(s_{i}\) is periodic with period \(k_{i}\). Loosely speaking, \(s_{i}\) is periodic if it can only return to itself after a fixed number of transitions \(k_{i}>1\) (or multiple of \(k_{i}\)), else it is aperiodic.

If states \(s_{i}\) and \(s_{j}\) belong to the same communicating class, then they have the same period \(k_{i}\). As a consequence, each of the states of an irreducible DTMC share the same periodicity. This periodicity is also considered the DTMC periodicity. States can also be classified according to how often the chain returns to them. Let \(T^{x\rightarrow x}\) be the number of steps needed to go back to state \(x\), knowing that the chain starts in \(x\).

• A state \(x\) is recurrent if \(P(T^{x\rightarrow x}<+\infty)=1\) (equivalently \(P(T^{x\rightarrow x}=+\infty)=0\)). In addition:
  1. A state \(x\) is null recurrent if in addition \(E(T^{x\rightarrow x})=+\infty\).
  2. A state \(x\) is positive recurrent if in addition \(E(T^{x\rightarrow x})<+\infty\).
  3. A state \(x\) is absorbing if in addition \(P(T^{x\rightarrow x}=1)=1\).
• A state \(x\) is transient if \(P(T^{x\rightarrow x}<+\infty)<1\) (equivalently \(P(T^{x\rightarrow x}=+\infty)>0\)).

It is possible to analyze the timing to reach a certain state. The first passage time (or hitting time) from state \(s_{i}\) to state \(s_{j}\) is the number \(T_{ij}\) of steps taken by the chain until it arrives for the first time to state \(s_{j}\), given that \(X_{0} = s_{i}\). The probability distribution of \(T_{ij}\) is defined by Equation \(\ref{eq:fpt1}\)

\[\begin{equation} {h_{ij}}^{\left( n \right)} = Pr\left( {T_{ij} = n} \right) = Pr\left( X_n = s_j,X_{n - 1} \ne s_{j}, \ldots ,X_1 \ne s_j |X_0 = s_i \right) \label{eq:fpt1} \end{equation}\]

and can be found recursively using Equation \(\ref{eq:ftp2}\), given that \({h_{ij}}^{\left( n \right)} = p_{ij}\).

\[\begin{equation} {h_{ij}}^{\left( n \right)} = \sum\limits_{k \in S - \left\{ s_{j} \right\}}^{} {{p_{ik}}{h_{kj}}^{\left( {n - 1} \right)}}. \label{eq:ftp2} \end{equation}\]

A commonly used quantity related to \(h\) is its average value, i.e. the (also expected hitting time), namely \(\bar h_{ij}= \sum_{n=1}^{\infty} n \,h_{ij}^{(n)}\).

If in the definition of the first passage time we let \(s_{i}=s_{j}\), we obtain the first recurrence time \(T_{i}=\inf \{ n\geq1:X_{n}=s_{i}|X_{0}=s_{i} \}\). The mean recurrence time of state \(s_{i}\) is the expected value of \(T_{i}\):

\[ r_i = \sum_{k = 1}^{\infty} k \cdot P(T_i = k). \]

Revisiting the definition of recurrence and transience: a state \(s_{i}\) is said to be recurrent if it is visited infinitely often, i.e., \(Pr(T_{i}<+\infty|X_{0}=s_{i})=1\). On the opposite, \(s_{i}\) is called transient if there is a positive probability that the chain will never return to \(s_{i}\), i.e., \(Pr(T_{i}=+\infty|X_{0}=s_{i})>0\).

Given a time homogeneous Markov chain with transition matrix , a stationary distribution is a stochastic row vector such that \(z=z\cdot P\), where \(0\leq z_{j}\leq 1 \: \forall j\) and \(\sum_{j}z_{j}=1\).

If a DTMC \(\{X_{n}\}\) is irreducible and aperiodic, then it has a limit distribution and this distribution is stationary. As a consequence, if \(P\) is the \(k\times k\) transition matrix of the chain and \(z=\left(z_{1},...,z_{k}\right)\) is the unique eigenvector of \(P\) such that \(\sum_{i=1}^{k}z_{i}=1\), then we get

\[\begin{equation} \underset{n\rightarrow\infty}{lim}P^{n}=Z, \label{eq:limMc} \end{equation}\]

where \(Z\) is the matrix having all rows equal to \(z\). The stationary distribution of \(\{X_{n}\}\) is represented by \(z\).

A matrix \(A\) is called positive if all of its entries are strictly positive, and we write \(A > 0\). If some power of the transition matrix \(P\) of a DTMC is positive, i.e. there exists \(m > 0\) such that \(P^m > 0\), then \(P\) is called primitive and the DTMC is said to be regular. For a finite state space, being regular is equivalent to being irreducible and aperiodic. All regular DTMCs are irreducible, but the converse is not true.

Given two disjoint sets of states \(A\) and \(B\), the \(q_j^{(AB)}\) is the probability that a chain starting in state \(s_j\) reaches \(A\) before \(B\) (Noé et al. 2009). For \(j \notin A \cup B\) it solves

\[\begin{equation} q_j^{(AB)} = \sum_{k \in S} P_{jk}q_k^{(AB)} \quad \mbox{with} \quad q_k^{(AB)} = 1 \mbox{ for } k \in A \quad \mbox{and} \quad q_k^{(AB)} = 0 \mbox{ for } k \in B. \end{equation}\]

The from \(i\) to \(j\) is the probability of ever reaching the state \(j\) if the initial state is \(i\):

\[\begin{equation} h_{i,j} = Pr(T_{ij} < \infty) = \sum_{n = 0}^{\infty} h_{ij}^{(n)} \label{eq:hitting-probs} \end{equation}\]

In a DTMC with a finite set of states, every transient state can reach at least one recurrent state. If the chain starts in a transient state, once it hits a recurrent state it is caught in the corresponding recurrent class and never goes back to the initial state. Given a transient state \(i\), the \(f^{*}_{ij}\) into the recurrent state \(j\) is the probability that \(j\) is the first recurrent state visited by the chain (which is then absorbed by the recurrent class of \(j\)). The \(b_i\) is the mean number of steps the chain takes, starting from the transient state \(i\), until it hits any recurrent state..

A short example

Consider the following numerical example. Suppose we have a DTMC with a set of 3 possible states \(S=\{s_{1}, s_{2}, s_{3}\}\). Let the transition matrix be: \[\begin{equation} P=\left[\begin{array}{ccc} 0.5 & 0.2 & 0.3\\ 0.15 & 0.45 & 0.4\\ 0.25 & 0.35 & 0.4 \end{array}\right]. \label{eq:trPropExEx1} \end{equation}\]

In \(P\), \(p_{11}=0.5\) is the probability that \(X_{1}=s_{1}\) given that \(X_{0}=s_{1}\), and so on. It is easy to see that the chain is irreducible since all the states communicate (it is made by one communicating class only).

Suppose that the current state of the chain is \(X_{0}=s_{2}\), i.e., \(x^{(0)}=(0\:1\:0)\), then the probability distribution of states after 1 and 2 steps can be computed as shown in Equations @ref(eq:trPropExEx2) and @ref(eq:trPropExEx3).

\[\begin{equation} x^{(1)}=\left(0\:1\:0\right)\left[\begin{array}{ccc} 0.5 & 0.2 & 0.3\\ 0.15 & 0.45 & 0.4\\ 0.25 & 0.35 & 0.4 \end{array}\right]=\left(0.15\:0.45\:0.4\right). \label{eq:trPropExEx2} \end{equation}\]

\[\begin{equation} x^{(n)}=x^{(n-1)}P \to \left(0.15\:0.45\:0.4\right)\left[\begin{array}{ccc} 0.5 & 0.2 & 0.3\\ 0.15 & 0.45 & 0.4\\ 0.25 & 0.35 & 0.4 \end{array}\right]=\left(0.2425\:0.3725\:0.385\right). \label{eq:trPropExEx3} \end{equation}\]

If we were interested in the probability of being in the state \(s_{3}\) in the second step, then \(Pr\left(X_{2}=s_{3}\left|X_{0}=s_{2}\right.\right)=0.385\).

Creating markovchain objects

The package is loaded within the command line as follows:

library("markovchain")

The markovchain and markovchainList S4 classes (Chambers 2008) are defined within the package as displayed:

## Class "markovchain" [package "markovchain"]
## 
## Slots:
##                                                                           
## Name:            states            byrow transitionMatrix             name
## Class:        character          logical           matrix        character

## Class "markovchainList" [package "markovchain"]
## 
## Slots:
##                                 
## Name:  markovchains         name
## Class:         list    character

The first class has been designed to handle homogeneous Markov chain processes, while the latter (which is itself a list of markovchain objects) has been designed to handle semi-homogeneous Markov chains processes.

Every markovchain object has the following slots:

1. states: a character vector, listing the states for which transition probabilities are defined.
2. byrow: a logical element, indicating whether transition probabilities are shown by row or by column.
3. transitionMatrix: the probabilities of the transition matrix.
4. name: optional character element to name the DTMC.

markovchainList objects have the following slots:

1. markovchains: a list of markovchain objects.
2. name: optional character element to name the DTMC.

A markovchain object can be created either in the long way, as the following code shows,

weatherStates <- c("sunny", "cloudy", "rain")
byRow <- TRUE
weatherMatrix <- matrix(data = c(0.70, 0.2, 0.1,
                       0.3, 0.4, 0.3,
                       0.2, 0.45, 0.35), byrow = byRow, nrow = 3,
                     dimnames = list(weatherStates, weatherStates))
mcWeather <- new("markovchain", states = weatherStates, byrow = byRow, 
               transitionMatrix = weatherMatrix, name = "Weather")

or in a shorter way, as displayed below

mcWeather <- new("markovchain", states = c("sunny", "cloudy", "rain"),
                 transitionMatrix = matrix(data = c(0.70, 0.2, 0.1,
                       0.3, 0.4, 0.3,
                       0.2, 0.45, 0.35), byrow = byRow, nrow = 3), 
                 name = "Weather")

When new("markovchain") is called alone, a default Markov chain is created.

defaultMc <- new("markovchain")

The shorter way is made possible by the initialize S4 method, which checks that:

• the transitionMatrix, of class matrix or Matrix, is a transition matrix, i.e., all its entries are probabilities and either all rows or all columns sum up to one;
• the row and column names of transitionMatrix, if present, coincide with the states slot; if they are missing, they are taken from states.

markovchain objects can be collected in a markovchainList S4 object, as the following example shows.

mcList <- new("markovchainList", markovchains = list(mcWeather, defaultMc), 
                  name = "A list of Markov chains")

Handling markovchain objects

Table @ref(tab:methodsToHandleMc) lists which methods handle and manipulate markovchain objects.

The examples that follow show how operations on markovchain objects can be easily performed. For instance, using the previously defined chain we can find the probability distribution of the weather states in two and in seven days, given that today is cloudy.

initialState <- c(0, 1, 0)
after2Days <- initialState * (mcWeather * mcWeather)
after7Days <- initialState * (mcWeather ^ 7)
after2Days

##      sunny cloudy  rain
## [1,]  0.39  0.355 0.255

round(after7Days, 3)

##      sunny cloudy  rain
## [1,] 0.462  0.319 0.219

The same answer can be obtained by defining the vector of probabilities as a column vector. A column-defined transition matrix can be set up either by creating a new matrix or by transposing an existing markovchain object with the t method.

initialState <- c(0, 1, 0)
after2Days <- (t(mcWeather) * t(mcWeather)) * initialState
after7Days <- (t(mcWeather) ^ 7) * initialState
after2Days

##         [,1]
## sunny  0.390
## cloudy 0.355
## rain   0.255

round(after7Days, 3)

##         [,1]
## sunny  0.462
## cloudy 0.319
## rain   0.219

The initial state vector does not need to be a probability vector, as the code that follows shows, where it holds the number of individuals in each state:

fvals<-function(mchain,initialstate,n) {
  out<-data.frame()
  names(initialstate)<-names(mchain)
  for (i in 0:n)
  {
    iteration<-initialstate*mchain^(i)
    out<-rbind(out,iteration)
  }
  out<-cbind(out, i=seq(0,n))
  out<-out[,c(4,1:3)]
  return(out)
}
fvals(mchain=mcWeather,initialstate=c(90,5,5),n=4)

##   i    sunny   cloudy     rain
## 1 0 90.00000  5.00000  5.00000
## 2 1 65.50000 22.25000 12.25000
## 3 2 54.97500 27.51250 17.51250
## 4 3 50.23875 29.88062 19.88062
## 5 4 48.10744 30.94628 20.94628

Basic methods return the states and the dimension of the transition matrix.

states(mcWeather)

## [1] "sunny"  "cloudy" "rain"

names(mcWeather)

## [1] "sunny"  "cloudy" "rain"

dim(mcWeather)

## [1] 3

The name of a markovchain object can be read and changed with name and name<-.

name(mcWeather)

## [1] "Weather"

name(mcWeather) <- "New Name"
name(mcWeather)

## [1] "New Name"

The states, and the transition matrix accordingly, can be sorted alphabetically:

markovchain:::sort(mcWeather)

## New Name 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##        cloudy rain sunny
## cloudy   0.40 0.30   0.3
## rain     0.45 0.35   0.2
## sunny    0.20 0.10   0.7

Transition probabilities can be accessed directly through the transitionProbability method and the "[" method.

transitionProbability(mcWeather, "cloudy", "rain")

## [1] 0.3

mcWeather[2,3]

## [1] 0.3

The transition matrix of a markovchain object can be displayed with the print or show methods (the latter being less verbose). The transition diagram can be plotted with the plot method (Figure @ref(fig:mcPlot)), a wrapper of plot.igraph from the package (Csardi and Nepusz 2006); additional parameters can be passed to plot to control the layout of the graph. The and packages can be used instead (Figure @ref(fig:mcPlotdiagram)). The plot method uses communicatingClasses to separate the states of different communicating classes, and all the states of a class have the same color. The graph can also be exported as text for Graphviz or Mermaid, as shown in Section @ref(sec:serialization).

print(mcWeather)

##        sunny cloudy rain
## sunny    0.7   0.20 0.10
## cloudy   0.3   0.40 0.30
## rain     0.2   0.45 0.35

show(mcWeather)

## New Name 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  sunny, cloudy, rain 
##  The transition matrix  (by rows)  is defined as follows: 
##        sunny cloudy rain
## sunny    0.7   0.20 0.10
## cloudy   0.3   0.40 0.30
## rain     0.2   0.45 0.35

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

Weather example. Markov chain plot

## Loading required package: shape

Weather example. Markov chain plot with diagram

Coercion from and to some other classes is possible, as shown in Figure @ref(fig:fromAndTo) and in the following code.

mcDf <- as(mcWeather, "data.frame")
mcNew <- as(mcDf, "markovchain")
mcDf

##       t0     t1 prob
## 1  sunny  sunny 0.70
## 2  sunny cloudy 0.20
## 3  sunny   rain 0.10
## 4 cloudy  sunny 0.30
## 5 cloudy cloudy 0.40
## 6 cloudy   rain 0.30
## 7   rain  sunny 0.20
## 8   rain cloudy 0.45
## 9   rain   rain 0.35

mcIgraph <- as(mcWeather, "igraph")

if (requireNamespace("msm", quietly = TRUE)) {
require(msm)
Q <- rbind ( c(0, 0.25, 0, 0.25),
             c(0.166, 0, 0.166, 0.166),
             c(0, 0.25, 0, 0.25),
             c(0, 0, 0, 0) )
cavmsm <- msm(state ~ years, subject = PTNUM, data = cav, qmatrix = Q, death = 4)
msmMc <- as(cavmsm, "markovchain")
msmMc
  } else {
  message("msm unavailable")
  }

## Loading required package: msm

## Unnamed Markov chain 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  State 1, State 2, State 3, State 4 
##  The transition matrix  (by rows)  is defined as follows: 
##             State 1    State 2    State 3    State 4
## State 1 0.853958721 0.08836953 0.01475543 0.04291632
## State 2 0.155576908 0.56663284 0.20599563 0.07179462
## State 3 0.009903994 0.07853691 0.65965727 0.25190183
## State 4 0.000000000 0.00000000 0.00000000 1.00000000

The following code coerces an object of the package into a markovchain ( was archived from CRAN in September 2020, so the example runs only if the package is installed):

if (requireNamespace("etm", quietly = TRUE)) {
library(etm)
data(sir.cont)
sir.cont <- sir.cont[order(sir.cont$id, sir.cont$time), ]
for (i in 2:nrow(sir.cont)) {
  if (sir.cont$id[i]==sir.cont$id[i-1]) {
    if (sir.cont$time[i]==sir.cont$time[i-1]) {
      sir.cont$time[i-1] <- sir.cont$time[i-1] - 0.5
    }
  }
}
tra <- matrix(ncol=3,nrow=3,FALSE)
tra[1, 2:3] <- TRUE
tra[2, c(1, 3)] <- TRUE
tr.prob <- etm::etm(sir.cont, c("0", "1", "2"), tra, "cens", 1)
tr.prob
etm2mc<-as(tr.prob, "markovchain")
etm2mc
  } else {
  message("etm unavailable")
}

## etm unavailable

The markovchain methods for import and export

A square probability matrix can also be coerced into a markovchain object, which is another way to create one, as the code below shows.

myMatr<-matrix(c(.1,.8,.1,.2,.6,.2,.3,.4,.3), byrow=TRUE, ncol=3)
myMc<-as(myMatr, "markovchain")
myMc

## Unnamed Markov chain 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  s1, s2, s3 
##  The transition matrix  (by rows)  is defined as follows: 
##     s1  s2  s3
## s1 0.1 0.8 0.1
## s2 0.2 0.6 0.2
## s3 0.3 0.4 0.3

Semi-homogeneous Markov chains are handled by markovchainList objects. The example that follows arises from health insurance, where the costs associated with the patients of a Continuous Care Health Community (CCHC) are modeled by a semi-homogeneous Markov chain, since the transition probabilities change every year. The methods explicitly written for markovchainList objects are print, show, dim and [.

stateNames = c("H", "I", "D")
Q0 <- new("markovchain", states = stateNames, 
        transitionMatrix =matrix(c(0.7, 0.2, 0.1,0.1, 0.6, 0.3,0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t0")
Q1 <- new("markovchain", states = stateNames, 
        transitionMatrix = matrix(c(0.5, 0.3, 0.2,0, 0.4, 0.6,0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t1")
Q2 <- new("markovchain", states = stateNames, 
        transitionMatrix = matrix(c(0.3, 0.2, 0.5,0, 0.2, 0.8,0, 0, 1), 
        byrow = TRUE,nrow = 3), name = "state t2")
Q3 <- new("markovchain", states = stateNames, 
          transitionMatrix = matrix(c(0, 0, 1, 0, 0, 1, 0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t3")
mcCCRC <- new("markovchainList",markovchains = list(Q0,Q1,Q2,Q3), 
      name = "Continuous Care Health Community")
print(mcCCRC)

## Continuous Care Health Community  list of Markov chain(s) 
## Markovchain  1 
## state t0 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  H, I, D 
##  The transition matrix  (by rows)  is defined as follows: 
##     H   I   D
## H 0.7 0.2 0.1
## I 0.1 0.6 0.3
## D 0.0 0.0 1.0
## 
## Markovchain  2 
## state t1 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  H, I, D 
##  The transition matrix  (by rows)  is defined as follows: 
##     H   I   D
## H 0.5 0.3 0.2
## I 0.0 0.4 0.6
## D 0.0 0.0 1.0
## 
## Markovchain  3 
## state t2 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  H, I, D 
##  The transition matrix  (by rows)  is defined as follows: 
##     H   I   D
## H 0.3 0.2 0.5
## I 0.0 0.2 0.8
## D 0.0 0.0 1.0
## 
## Markovchain  4 
## state t3 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  H, I, D 
##  The transition matrix  (by rows)  is defined as follows: 
##   H I D
## H 0 0 1
## I 0 0 1
## D 0 0 1

The elements of a markovchainList can be accessed directly, and dim returns the number of markovchain objects it contains.

mcCCRC[[1]]

## state t0 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  H, I, D 
##  The transition matrix  (by rows)  is defined as follows: 
##     H   I   D
## H 0.7 0.2 0.1
## I 0.1 0.6 0.3
## D 0.0 0.0 1.0

dim(mcCCRC)

## [1] 4

The package contains some data found in the literature on DTMC models (see Section @ref(sec:applications)). Table @ref(tab:datasets) lists the datasets included in the current release.

Finally, Table @ref(tab:demos) lists the demos included in the directory of the package.

Conditional distributions

The conditional distribution of the weather states, given that today is sunny, is computed by the following code.

conditionalDistribution(mcWeather, "sunny")

##  sunny cloudy   rain 
##    0.7    0.2    0.1

Stationary states

A stationary (steady state, or equilibrium) vector is a probability vector such that Equation \(\ref{eq:steadystat2}\) holds

\[\begin{equation} \begin{matrix} 0\leq \pi_j \leq 1\\ \sum_{j \in S} \pi_j = 1\\ \pi \cdot P = \pi \end{matrix} \label{eq:steadystat2} \end{equation}\]

Steady states are associated with the eigenvalues of \(P\) equal to one. One could be tempted to compute them by solving the eigenproblem of the matrix, taking the real parts of the eigenvectors (if \(u + iv\) is an eigenvector of \(P\) for the eigenvalue one, then so are \(u = Re(u + iv)\) and \(v = Im(u + iv)\)) and normalizing them by their sum. This approach, however, raises some concerns:

1. If \(u, v \in \mathbb{R}^n\) are linearly independent eigenvectors associated with the eigenvalue \(1\), \(u + iv\) and \(u - iv\) are also linearly independent eigenvectors, and their real parts coincide. Taking real parts, we would lose an eigenvector, because we cannot know in advance whether the algorithm that computes the eigenvectors will return something similar to what we described. We should be agnostic to the underlying eigenvector computation algorithm.

2. Imagine the identity \(P\) of dimensions \(2 \times 2\). Its eigenvectors associated to the \(1\) eigen value are \(u = (1, 0)\) and \(v = (0, 1)\). However, the underlying algorithm to compute eigenvectors could return \((1, -2)\) and \((-2, 1)\) instead, that are linear combinations of the aforementioned ones, and therefore eigenvectors. Normalizing by their sum, we would get: \((-1, 2)\) and \((2, -1)\), which obviously are not probability measures. Again, we should be agnostic to the underlying eigen computation algorithm.

3. Algorithms to compute eigenvalues and eigenvectors are computationally expensive: they are iterative, and we cannot predict a fixed number of iterations for them. Moreover, each iteration takes \(\mathcal{O}(m^2)\) or \(\mathcal{O}(m^3)\) algorithmic complexity, with \(m\) the number of states.

We use the fact that every irreducible DTMC has a unique steady state: if \(M\) is the transition matrix of an irreducible DTMC (all states communicate with each other), then there exists a unique \(v \in \mathbb{R}^m\) such that

\[ v \cdot M = v, \qquad \sum_{i = 1}^m v_i = 1 \]

We also use the fact that a steady state of a DTMC assigns \(0\) to the transient states. The canonical form of a row-stochastic matrix looks as follows:

\[ \left(\begin{array}{c|c|c|c|c} M_1 & 0 & 0 & \ldots & 0 \\ \hline 0 & M_2 & 0 & \ldots & 0 \\ \hline 0 & 0 & M_3 & \ldots & 0 \\ \hline \vdots & \vdots & \vdots & \ddots & \vdots \\ \hline A_1 & A_2 & A_3 & \ldots & R \end{array}\right) \]

where \(M_i\) corresponds to irreducible sub-chains, the blocks \(A_i\) correspond to the transitions from transient states to each of the recurrent classes and \(R\) are the transitions from the transient states to themselves.

Moreover, a Markov chain has exactly as many steady states as recurrent classes. Therefore, we have coded the following algorithm ¹:

1. Identify the recurrent classes \([C_1, \ldots, C_l]\) with the function.
2. For each class \(C_i\), extract the corresponding sub-matrix \(M_i\).
3. Solve the system \(v \cdot C_i = v, \, \sum_{j = 1}^{|C_i|} v_j = 1\) which has a unique solution, for each \(i = 1, \ldots, l\).
4. Map each state \(v_i\) back to the original order in \(P\) and assign \(0\) to the slots corresponding to transient states.

The result is returned in matrix form, one steady state per row.

steadyStates(mcWeather)

##          sunny    cloudy      rain
## [1,] 0.4636364 0.3181818 0.2181818

A Markov chain can have more than one stationary distribution, as the gambler’s ruin example shows.

gamblerRuinMarkovChain <- function(moneyMax, prob = 0.5) {
  m <- markovchain:::zeros(moneyMax + 1)
  m[1,1] <- m[moneyMax + 1,moneyMax + 1] <- 1
  states <- as.character(0:moneyMax)
  rownames(m) <- colnames(m) <- states
  
  for(i in 2:moneyMax){ 
    m[i,i-1] <- 1 - prob
    m[i, i + 1] <- prob   
  }
  
  new("markovchain", transitionMatrix = m, 
      name = paste("Gambler ruin", moneyMax, "dim", sep = " "))
}

mcGR4 <- gamblerRuinMarkovChain(moneyMax = 4, prob = 0.5)
steadyStates(mcGR4)

##      0 1 2 3 4
## [1,] 0 0 0 0 1
## [2,] 1 0 0 0 0

Entropy rate

For a stationary first-order Markov chain, the entropy rate is the conditional entropy of the next state given the current state:

\[ H(X) = H(X_{t+1}\mid X_t) = -\sum_i \pi_i\sum_j p_{ij}\log_b p_{ij}, \]

where \(\pi\) is the stationary distribution and \(0\log(0)\) is interpreted as zero. With \(b=2\) the result is measured in bits per transition; with \(b=e\) it is measured in nats per transition . The method entropyRate uses the unique stationary distribution of a finite irreducible chain. Aperiodicity is not required; reducible chains are rejected because the result can depend on the chosen stationary distribution .

entropyRate(mcWeather)

## [1] 1.366258

entropyRate(mcWeather, base = exp(1))

## [1] 0.9470176

A deterministic chain has entropy rate zero, whereas a chain whose rows are all uniform on \(n\) states has the maximum value \(\log_b(n)\).

topologicalEntropy(mcWeather)

## [1] 1.584963

normalizedEntropyRate(mcWeather)

## [1] 0.8620126

Classification of states

The absorbing states are returned by the absorbingStates method.

absorbingStates(mcGR4)

## [1] "0" "4"

absorbingStates(mcWeather)

## character(0)

The methods that need to know the communicating classes, the recurrent states or the transient states rely on the internal function .commclassKernel, a modification of Tarjan’s algorithm . It takes a transition matrix of dimension \(n\) and returns a list of two items:

1. classes, a matrix whose \((i, j)\) entry is TRUE if \(s_i\) and \(s_j\) are in the same communicating class;
2. closed, a vector whose \(i\)-th entry indicates whether the communicating class to which \(s_i\) belongs is closed.

The summary method for markovchain objects is built on two internal functions that use it.

The example matrix of illustrates the purpose of the function.

P <- markovchain:::zeros(10)
P[1, c(1, 3)] <- 1/2;
P[2, 2] <- 1/3; P[2,7] <- 2/3;
P[3, 1] <- 1;
P[4, 5] <- 1;
P[5, c(4, 5, 9)] <- 1/3;
P[6, 6] <- 1;
P[7, 7] <- 1/4; P[7,9] <- 3/4;
P[8, c(3, 4, 8, 10)] <- 1/4;
P[9, 2] <- 1;
P[10, c(2, 5, 10)] <- 1/3;
rownames(P) <- letters[1:10] 
colnames(P) <- letters[1:10]
probMc <- new("markovchain", transitionMatrix = P, 
              name = "Probability MC")
summary(probMc)

## Probability MC  Markov chain that is composed by: 
## Closed classes: 
## a c 
## b g i 
## f 
## Recurrent classes: 
## {a,c},{b,g,i},{f}
## Transient classes: 
## {d,e},{h},{j}
## The Markov chain is not irreducible 
## The absorbing states are: f

With details = TRUE, summary also prints, after the output above, a block of further properties in the spirit of the printout of a PyDTMC chain: size and rank of the transition matrix, classes, irreducibility, period, regularity, whether the chain is absorbing, reversible, stochastically monotone or symmetric, and, when they are defined, entropy rate, SLEM, spectral gap and Kemeny’s constant. The values are also returned, invisibly, in the details element of the result.

summary(mcWeather, details = TRUE)

## New Name  Markov chain that is composed by: 
## Closed classes: 
## sunny cloudy rain 
## Recurrent classes: 
## {sunny,cloudy,rain}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE
## Further properties: 
##   Size                    : 3 
##   Rank                    : 2 
##   Communicating classes   : 1 (1 recurrent, 0 transient) 
##   Irreducible             : yes 
##   Period                  : 1 
##   Regular (ergodic)       : yes 
##   Absorbing chain         : no 
##   Reversible              : no 
##   Stochastically monotone : yes 
##   Symmetric               : no 
##   Entropy rate (bits)     : 1.36626 
##   SLEM                    : 0.45 
##   Spectral gap            : 0.55 
##   Kemeny constant         : 2.81818

The states that belong to a transient class are called transient, and the transientStates method returns them.

transientStates(probMc)

## [1] "d" "e" "h" "j"

The canonicForm method turns a Markov chain into its canonical form, reordering the states so that the recurrent classes come first and the transient states last.

probMcCanonic <- canonicForm(probMc)
probMc

## Probability MC 
##  A  10 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c, d, e, f, g, h, i, j 
##  The transition matrix  (by rows)  is defined as follows: 
##     a         b    c         d         e f         g    h         i         j
## a 0.5 0.0000000 0.50 0.0000000 0.0000000 0 0.0000000 0.00 0.0000000 0.0000000
## b 0.0 0.3333333 0.00 0.0000000 0.0000000 0 0.6666667 0.00 0.0000000 0.0000000
## c 1.0 0.0000000 0.00 0.0000000 0.0000000 0 0.0000000 0.00 0.0000000 0.0000000
## d 0.0 0.0000000 0.00 0.0000000 1.0000000 0 0.0000000 0.00 0.0000000 0.0000000
## e 0.0 0.0000000 0.00 0.3333333 0.3333333 0 0.0000000 0.00 0.3333333 0.0000000
## f 0.0 0.0000000 0.00 0.0000000 0.0000000 1 0.0000000 0.00 0.0000000 0.0000000
## g 0.0 0.0000000 0.00 0.0000000 0.0000000 0 0.2500000 0.00 0.7500000 0.0000000
## h 0.0 0.0000000 0.25 0.2500000 0.0000000 0 0.0000000 0.25 0.0000000 0.2500000
## i 0.0 1.0000000 0.00 0.0000000 0.0000000 0 0.0000000 0.00 0.0000000 0.0000000
## j 0.0 0.3333333 0.00 0.0000000 0.3333333 0 0.0000000 0.00 0.0000000 0.3333333

probMcCanonic

## Probability MC 
##  A  10 - dimensional discrete Markov Chain defined by the following states: 
##  a, c, b, g, i, f, d, e, h, j 
##  The transition matrix  (by rows)  is defined as follows: 
##     a    c         b         g         i f         d         e    h         j
## a 0.5 0.50 0.0000000 0.0000000 0.0000000 0 0.0000000 0.0000000 0.00 0.0000000
## c 1.0 0.00 0.0000000 0.0000000 0.0000000 0 0.0000000 0.0000000 0.00 0.0000000
## b 0.0 0.00 0.3333333 0.6666667 0.0000000 0 0.0000000 0.0000000 0.00 0.0000000
## g 0.0 0.00 0.0000000 0.2500000 0.7500000 0 0.0000000 0.0000000 0.00 0.0000000
## i 0.0 0.00 1.0000000 0.0000000 0.0000000 0 0.0000000 0.0000000 0.00 0.0000000
## f 0.0 0.00 0.0000000 0.0000000 0.0000000 1 0.0000000 0.0000000 0.00 0.0000000
## d 0.0 0.00 0.0000000 0.0000000 0.0000000 0 0.0000000 1.0000000 0.00 0.0000000
## e 0.0 0.00 0.0000000 0.0000000 0.3333333 0 0.3333333 0.3333333 0.00 0.0000000
## h 0.0 0.25 0.0000000 0.0000000 0.0000000 0 0.2500000 0.0000000 0.25 0.2500000
## j 0.0 0.00 0.3333333 0.0000000 0.0000000 0 0.0000000 0.3333333 0.00 0.3333333

The function is.accessible checks whether a state \(s_{j}\) is accessible from state \(s_i\), that is whether the probability of eventually reaching \(s_j\) starting from \(s_{i}\) is greater than zero.

is.accessible(object = probMc, from = "a", to = "c")

## [1] TRUE

is.accessible(object = probMc, from = "g", to = "c")

## [1] FALSE

In Section @ref(sec:properties) we observed that, if a DTMC is irreducible, all its states share the same periodicity. The period function therefore returns the period of a DTMC, provided that it is irreducible. The example that follows checks with is.irreducible whether a DTMC is irreducible and, if so, computes its period with period.

E <- matrix(0, nrow = 4, ncol = 4)
E[1, 2] <- 1
E[2, 1] <- 1/3; E[2, 3] <- 2/3
E[3,2] <- 1/4; E[3, 4] <- 3/4
E[4, 3] <- 1

mcE <- new("markovchain", states = c("a", "b", "c", "d"), 
        transitionMatrix = E, 
        name = "E")
is.irreducible(mcE)

## [1] TRUE

period(mcE)

## [1] 2

The next example is the Markov chain from the website , plotted in Figure @ref(fig:mcMathematics).

mathematicaMatr <- markovchain:::zeros(5)
mathematicaMatr[1,] <- c(0, 1/3, 0, 2/3, 0)
mathematicaMatr[2,] <- c(1/2, 0, 0, 0, 1/2)
mathematicaMatr[3,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[4,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[5,] <- c(0, 0, 0, 0, 1)
statesNames <- letters[1:5]
mathematicaMc <- new("markovchain", transitionMatrix = mathematicaMatr,
                   name = "Mathematica MC", states = statesNames)

Mathematica 9 example. Markov chain plot.

## Mathematica MC  Markov chain that is composed by: 
## Closed classes: 
## c d 
## e 
## Recurrent classes: 
## {c,d},{e}
## Transient classes: 
## {a,b}
## The Markov chain is not irreducible 
## The absorbing states are: e

First passage time distributions and means

provides code to compute first passage time (within \(1,2,\ldots, n\) steps) given the initial state to be \(i\). The listings translated into on which the firstPassage function is based are:

.firstpassageKernel <- function(P, i, n){
  G <- P
  H <- P[i,]
  E <- 1 - diag(size(P)[2])
  for (m in 2:n) {
    G <- P %*% (G * E)
    H <- rbind(H, G[i,])
  }
  return(H)
}

We conclude that the probability for the first rainy day to be the third one, given that the current state is sunny, is given by:

firstPassagePdF <- firstPassage(object = mcWeather, state = "sunny", 
                                n = 10)
firstPassagePdF[3, 3]

## [1] 0.121

To compute the mean first passage times, i.e. the expected number of days before it rains given that today is sunny, we can use the meanFirstPassageTime function:

meanFirstPassageTime(mcWeather)

##           sunny   cloudy     rain
## sunny  0.000000 4.285714 6.666667
## cloudy 3.725490 0.000000 5.000000
## rain   4.117647 2.857143 0.000000

indicating e.g. that the average number of days of sun or cloud before rain is 6.67 if we start counting from a sunny day, and 5 if we start from a cloudy day. Note that we can also specify one or more destination states:

meanFirstPassageTime(mcWeather,"rain")

##    sunny   cloudy 
## 6.666667 5.000000

The implementation follows the matrix solutions by (Grinstead and Snell 2006). We can check the result by averaging the first passage probability density function:

firstPassagePdF.long <- firstPassage(object = mcWeather, state = "sunny",  n = 100)
sum(firstPassagePdF.long[,"rain"] * 1:100)

## [1] 6.666664

Kemeny’s constant

For a finite irreducible chain with stationary distribution \(\pi\), let \(m_{ij}\) denote the mean time to hit state \(j\) from state \(i\), using \(m_{ii}=0\). Kemeny’s constant is

\[ K=\sum_j \pi_j m_{ij}. \]

Remarkably, \(K\) does not depend on the starting state \(i\). Equivalently, if \(Z=(I-P+\mathbf{1}\pi^T)^{-1}\), then \(K=\operatorname{tr}(Z)-1\) . The function kemenyConstant uses this formula and requires irreducibility, but not aperiodicity. If mean return times \(m_{jj}=1/\pi_j\) are placed on the diagonal instead, the weighted sum is \(K+1\).

kemenyConstant(mcWeather)

## [1] 2.818182

piWeather <- as.numeric(steadyStates(mcWeather))
meanFirstPassageTime(mcWeather) %*% piWeather

##            [,1]
## sunny  2.818182
## cloudy 2.818182
## rain   2.818182

The second result has one entry for each possible starting state; apart from numerical rounding, all entries equal Kemeny’s constant.

Sensitivity of the stationary distribution

Kemeny’s constant uses the fundamental matrix \(Z=(I-P+\mathbf{1}\pi^T)^{-1}\); the same matrix also gives the first-order sensitivity of the stationary distribution \(\pi\) to a small change in the chain’s transition probabilities. If row \(k\) of \(P\) is perturbed by a zero-sum vector \(d\) (so that \(P+e_k d^T\) is still row-stochastic to first order), the resulting change in \(\pi\) is, to first order, \[ \Delta\pi \approx d^T S, \qquad S_{lj} = \pi_k\bigl(Z_{lj}-\pi_j\bigr) \] . The function sensitivity(object, state) returns the matrix \(S\) for a single perturbed state k = state; it does not itself pick a direction \(d\), since that depends on which two transition probabilities are being traded off.

S <- sensitivity(mcWeather, "sunny")
S

##             sunny     cloudy        rain
## sunny   0.4479264 -0.2471450 -0.20078137
## cloudy -0.3529001  0.3850864 -0.03218633
## rain   -0.4371976 -0.0364012  0.47359880

For example, moving probability mass out of \(P(\text{sunny}\to\text{sunny})\) and into \(P(\text{sunny}\to\text{rain})\) (so \(d=(-\varepsilon,0,\varepsilon)\)) changes \(\pi\), to first order, by \(\varepsilon\) times the difference of the "rain" and "sunny" rows of \(S\):

eps <- 1e-4
predictedChange <- eps * (S["rain", ] - S["sunny", ])
predictedChange

##         sunny        cloudy          rain 
## -8.851240e-05  2.107438e-05  6.743802e-05

perturbedMatrix <- mcWeather@transitionMatrix
perturbedMatrix["sunny", "sunny"] <- perturbedMatrix["sunny", "sunny"] - eps
perturbedMatrix["sunny", "rain"] <- perturbedMatrix["sunny", "rain"] + eps
mcWeatherPerturbed <- new("markovchain", states = states(mcWeather),
                          transitionMatrix = perturbedMatrix)

as.numeric(steadyStates(mcWeatherPerturbed)) - as.numeric(steadyStates(mcWeather))

## [1] -8.849550e-05  2.107036e-05  6.742514e-05

The two agree to \(O(\varepsilon^2)\), as the first-order theory predicts. Unlike a positional, entry-by-entry sensitivity, S is built so that this holds for any zero-sum perturbation direction of row k, not only the one shown here.

Spectral diagnostics: SLEM, spectral gap and implied timescales

For a finite, irreducible chain, the Perron-Frobenius theorem guarantees a unique eigenvalue \(\lambda_1=1\); every other eigenvalue satisfies \(|\lambda_k|\le 1\). The second largest eigenvalue modulus (SLEM) is

\[ \mathrm{SLEM} = \max_{k>1} |\lambda_k|, \]

and the spectral gap is \(1-\mathrm{SLEM}\). Both are lightweight, eigenvalue-based diagnostics of how quickly a chain approaches its stationary distribution: a spectral gap close to \(1\) indicates fast mixing, while a spectral gap close to \(0\) indicates slow mixing . Only irreducibility is required; unlike some implementations, slem and spectralGap do not require aperiodicity, and correctly report \(\mathrm{SLEM}=1\) (spectral gap \(0\)) for a periodic chain, reflecting that such a chain never contracts towards its stationary distribution.

slem(mcWeather)

## [1] 0.45

spectralGap(mcWeather)

## [1] 0.55

Each non-trivial eigenvalue \(\lambda_k\), \(k>1\), also has its own implied timescale \(\tau_k=-1/\log|\lambda_k|\), measured in number of steps: it is the time the corresponding mode of the chain takes to decay by a factor of \(1/e\). The largest such value, \(\tau_2\), corresponds to the SLEM and is the slowest, dominant relaxation mode. The term comes from the Markov state model literature in molecular kinetics , where implied timescales are additionally compared across chains fitted at increasing lag times as a check of the Markov approximation; impliedTimescales evaluates a single already-fitted chain, leaving that comparison to the user.

impliedTimescales(mcWeather)

##     tau2     tau3 
## 1.252336 0.000000

As a textbook illustration with a known closed form, consider simple random walk on the cycle with \(n\) states, which moves to either neighbour with probability \(1/2\) . Its eigenvalues are \(\lambda_k = \cos(2\pi k/n)\), \(k=0,\ldots,n-1\), so the SLEM and spectral gap can be checked in closed form. For an odd \(n\) the walk is aperiodic. Since \(\cos\) is largest in modulus near \(0\) and near \(\pi\), the SLEM is the larger of the moduli attained near \(k=1\) (angle near \(0\)) and near \(k=n/2\) (angle near \(\pi\)); for \(n=7\) the latter dominates, at \(k=3\) (and, symmetrically, \(k=4\)):

randomWalkOnCycle <- function(n) {
  cycleP <- matrix(0, n, n)
  for (i in seq_len(n)) {
    cycleP[i, ((i - 2) %% n) + 1] <- 1 / 2
    cycleP[i, (i %% n) + 1] <- 1 / 2
  }
  new("markovchain", transitionMatrix = cycleP,
      states = as.character(seq_len(n)),
      name = paste("Random walk on the", n, "-cycle"))
}

mcCycle7 <- randomWalkOnCycle(7)
slem(mcCycle7)

## [1] 0.9009689

abs(cos(2 * pi * 3 / 7)) # closed-form SLEM, attained at k = 3 and k = 4

## [1] 0.9009689

spectralGap(mcCycle7)

## [1] 0.09903113

impliedTimescales(mcCycle7)

##      tau2      tau3      tau4      tau5      tau6      tau7 
## 9.5891459 9.5891459 2.1167476 2.1167476 0.6654537 0.6654537

For an even \(n\) the cycle graph is bipartite and the walk is periodic (period \(2\)): it can only return to its starting state after an even number of steps. Correspondingly \(\lambda_k\) reaches \(-1\) at \(k=n/2\), so the SLEM is \(1\) regardless of \(n\):

mcCycle8 <- randomWalkOnCycle(8)
slem(mcCycle8) # 1, attained at k = 4, where cos(2*pi*4/8) = cos(pi) = -1

## [1] 1

round(spectralGap(mcCycle8), 12) # numerically 0, up to floating-point error

## [1] 0

impliedTimescales(mcCycle8)["tau2"] # the k = 4 mode, sorted first, never decays

## tau2 
##  Inf

slem, spectralGap and impliedTimescales handle this periodic case directly, rather than treating it as an error: a spectral gap of (numerically) \(0\) and an infinite implied timescale are the mathematically correct description of a chain that does not converge pointwise to its stationary distribution. The tiny residual before rounding is ordinary floating-point error from eigen(), of the same order as .Machine$double.eps, not a sign that the chain is “almost” aperiodic.

Reversibility

A chain is reversible with respect to its stationary distribution \(\pi\) if the detailed balance equations hold:

\[ \pi_i P_{ij} = \pi_j P_{ji} \quad \text{for every } i,j. \]

Intuitively, a long stationary run of a reversible chain looks statistically the same whether played forwards or backwards in time: the flow of probability from \(i\) to \(j\) exactly balances the flow back . is.reversible checks this directly, and only requires irreducibility, not aperiodicity – reversibility and periodicity are independent properties.

A canonical family of reversible chains is simple random walk on an undirected graph (move to a uniformly random neighbour at each step): the symmetry of “being neighbours” is exactly what makes detailed balance hold. The random walk on a cycle from the previous section, whether the aperiodic mcCycle7 or the periodic mcCycle8, is reversible either way:

is.reversible(mcCycle7)

## [1] TRUE

is.reversible(mcCycle8)

## [1] TRUE

By contrast, a directed cycle – where probability only ever flows “one way around” – is not reversible, since there is a net clockwise flow at stationarity with nothing flowing back:

states3 <- c("a", "b", "c")
directedCycle <- new("markovchain", states = states3,
  transitionMatrix = matrix(c(0, 1, 0,
                              0, 0, 1,
                              1, 0, 0), byrow = TRUE, nrow = 3,
                            dimnames = list(states3, states3)))
is.reversible(directedCycle)

## [1] FALSE

The weather chain mcWeather is a more realistic example of a non-reversible chain: nothing forces the flow from “sunny” to “rain” to balance the flow back, and indeed it does not:

is.reversible(mcWeather)

## [1] FALSE

approximation <- closestReversible(directedCycle)
approximation$chain

## Unnamed Markov chain (closest reversible) 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c 
##  The transition matrix  (by rows)  is defined as follows: 
##     a   b   c
## a 0.0 0.5 0.5
## b 0.5 0.0 0.5
## c 0.5 0.5 0.0

is.reversible(approximation$chain)

## [1] TRUE

weatherReversible <- closestReversible(mcWeather)
weatherReversible$chain

## New Name (closest reversible) 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  sunny, cloudy, rain 
##  The transition matrix  (by rows)  is defined as follows: 
##            sunny    cloudy       rain
## sunny  0.7000000 0.2029412 0.09705882
## cloudy 0.2957143 0.4000000 0.30428571
## rain   0.2062500 0.4437500 0.35000000

weatherReversible$distance

## [1] 0.01074924

steadyStates(mcWeather)

##          sunny    cloudy      rain
## [1,] 0.4636364 0.3181818 0.2181818

steadyStates(weatherReversible$chain) # unchanged

##          sunny    cloudy      rain
## [1,] 0.4636364 0.3181818 0.2181818

Mixing time

Where spectralGap summarises convergence to stationarity through the transition matrix’s spectrum, mixingTime answers a more operational question directly: after how many steps is the chain’s distribution close to \(\pi\), no matter where it started? Formally, with total variation distance \(d(t) = \max_i \tfrac12\sum_j |P^t_{ij}-\pi_j|\), the mixing time is the smallest \(t\) with \(d(t)\le\varepsilon\), for a threshold \(\varepsilon\) (conventionally \(0.25\)) .

mixingTime(mcWeather)

## [1] 2

mixingTime(mcWeather, epsilon = 0.01) # a tighter threshold needs more steps

## [1] 6

Spectral gap and mixing time are connected by a standard bound :

\[ t_{\mathrm{mix}}(\varepsilon) \le \frac{1}{\mathrm{gap}} \log\!\left(\frac{1}{\varepsilon\,\pi_{\min}}\right), \]

which we can check directly against the values already computed for mcWeather in the previous section:

piWeather <- as.numeric(steadyStates(mcWeather))
upperBound <- log(1 / (0.25 * min(piWeather))) / spectralGap(mcWeather)
upperBound

## [1] 5.288583

mixingTime(mcWeather) <= upperBound

## [1] TRUE

mixingTime requires aperiodicity in addition to irreducibility, unlike slem, spectralGap and impliedTimescales: a periodic chain’s distribution never settles down at all (it keeps cycling), so “the number of steps until it is close to \(\pi\)” is not a meaningful question to ask of it. Calling it on the periodic mcCycle8 therefore raises an informative error rather than returning a misleading number:

mixingTime(mcCycle8)

## Error in `mixingTime()`:
## ! mixingTime is defined here only for aperiodic (ergodic) Markov chains, because a periodic chain never converges to its stationary distribution. See slem() or spectralGap() for a spectral diagnostic that remains meaningful for periodic chains.

Lazy chains

The standard fix for periodicity is to make the chain lazy: stay in place with probability \(\alpha\) at every step, and otherwise take a step of the original chain . For transition matrix \(P\), the lazy chain’s transition matrix is

\[ L = \alpha I + (1-\alpha) P. \]

Adding positive probability of staying in every state rules out any period greater than \(1\), while leaving the stationary distribution and reachability structure unchanged. Applying it to the periodic mcCycle8 from before:

lazyCycle8 <- lazyChain(mcCycle8, alpha = 0.5)
period(mcCycle8)

## [1] 2

period(lazyCycle8)

## [1] 1

steadyStates(mcCycle8)

##          1     2     3     4     5     6     7     8
## [1,] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

steadyStates(lazyCycle8) # unchanged by laziness

##          1     2     3     4     5     6     7     8
## [1,] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

Now that the chain is aperiodic, mixingTime is well-defined for it:

mixingTime(lazyCycle8)

## [1] 6

Laziness is not free, however: every non-trivial eigenvalue of \(L\) equals \(\alpha+(1-\alpha)\lambda\) for the corresponding eigenvalue \(\lambda\) of \(P\), so the whole non-trivial spectrum is pulled towards \(\alpha\) and mixing generally gets slower as \(\alpha\) grows:

spectralGap(mcCycle7) # aperiodic already, no laziness needed

## [1] 0.09903113

spectralGap(lazyChain(mcCycle7, alpha = 0.5))

## [1] 0.1882551

spectralGap(lazyChain(mcCycle7, alpha = 0.9))

## [1] 0.03765102

Evolution of a distribution

redistribute propagates an initial distribution through a given number of steps and returns the whole trajectory, one row per time point (including the initial time \(0\)). The initial distribution can be omitted (uniform), be a single state name (a point mass) or a numeric vector. relaxationTime is the reciprocal of the spectral gap: the characteristic number of steps needed for the distribution to approach stationarity.

redistribute(mcWeather, steps = 5, initial = "rain")

##       sunny    cloudy      rain
## 0 0.0000000 0.0000000 1.0000000
## 1 0.2000000 0.4500000 0.3500000
## 2 0.3450000 0.3775000 0.2775000
## 3 0.4102500 0.3448750 0.2448750
## 4 0.4396125 0.3301938 0.2301938
## 5 0.4528256 0.3235872 0.2235872

relaxationTime(mcWeather)

## [1] 1.818182

An observed sequence of states defines an observable \(f\) on the states, \(f_j\) being the number of times state \(j\) occurs in it. timeRelaxations returns \(\mu P^t f\), the expected value of \(f\) after \(t\) steps from the initial distribution \(\mu\), and timeCorrelations returns \(\sum_i \pi_i f_i (P^t g)_i\), the stationary correlation between \(f\) at time \(0\) and the observable \(g\) of a second sequence (by default the same) at time \(t\) . Both converge at the speed set by the second largest eigenvalue modulus, and they correspond to the functions of the same name in PyDTMC.

days <- c("sunny", "sunny", "rain", "cloudy", "sunny")
timeCorrelations(mcWeather, days, timePoints = c(0, 1, 2, 10))

##        0        1        2       10 
## 4.709091 4.152727 3.911636 3.714712

timeRelaxations(mcWeather, days, initial = "rain",
                timePoints = c(0, 1, 2, 10))

##        0        1        2       10 
## 1.000000 1.400000 1.690000 1.926874

The trajectory and the eigenvalues of the transition matrix can also be plotted with the type argument of autoplot (Figure \(\ref{fig:autoplot-flow}\)).

ggplot2::autoplot(mcWeather, type = "flow", steps = 10, initial = "rain")

Evolution of the distribution of the weather chain, starting from rain.

With type = "comparison", autoplot places two or more chains defined on the same states side by side, with the states matched by name. By default the transition matrices are drawn as heatmaps on a common probability scale; with what = "stationary" the stationary distributions are compared instead, which requires irreducible chains. As an example, a lazy version of the weather chain built with lazyChain has exactly the same stationary distribution but visibly different, “stickier” transitions (Figure \(\ref{fig:autoplot-comparison}\)).

mcLazyWeather <- lazyChain(mcWeather, alpha = 0.5)
mcLazyWeather@name <- "Lazy weather"
ggplot2::autoplot(mcWeather, type = "comparison", other = mcLazyWeather)

Transition matrices of the weather chain and of its lazy version with alpha = 0.5.

Mean recurrence time

The meanRecurrenceTime method gives the first mean recurrence time (expected number of steps to go back to a state if it was the initial one) for each recurrent state in the transition probabilities matrix for a DTMC. Let’s see an example:

meanRecurrenceTime(mcWeather)

##    sunny   cloudy     rain 
## 2.156863 3.142857 4.583333

Another example, with not all of its states being recurrent:

recurrentStates(probMc)

## [1] "a" "b" "c" "f" "g" "i"

meanRecurrenceTime(probMc)

##        f        b        g        i        a        c 
## 1.000000 2.555556 2.875000 3.833333 1.500000 3.000000

Absorption probabilities and mean absorption time

We are going to use the Drunkard’s random walk from (Grinstead and Snell 2006). We have a drunk person walking through the street. Each move the person does, if they have not arrived to either home (corner 1) or to the bar (corner 5) could be to the left corner or to the right one, with equal probability. In case of arrival to the bar or to home, the person stays there.

drunkProbs <- markovchain:::zeros(5)
drunkProbs[1,1] <- drunkProbs[5,5] <- 1
drunkProbs[2,1] <- drunkProbs[2,3] <- 1/2
drunkProbs[3,2] <- drunkProbs[3,4] <- 1/2
drunkProbs[4,3] <- drunkProbs[4,5] <- 1/2

drunkMc <- new("markovchain", transitionMatrix = drunkProbs)
drunkMc

## Unnamed Markov chain 
##  A  5 - dimensional discrete Markov Chain defined by the following states: 
##  1, 2, 3, 4, 5 
##  The transition matrix  (by rows)  is defined as follows: 
##     1   2   3   4   5
## 1 1.0 0.0 0.0 0.0 0.0
## 2 0.5 0.0 0.5 0.0 0.0
## 3 0.0 0.5 0.0 0.5 0.0
## 4 0.0 0.0 0.5 0.0 0.5
## 5 0.0 0.0 0.0 0.0 1.0

Recurrent (in fact absorbing states) are:

recurrentStates(drunkMc)

## [1] "1" "5"

Transient states are the rest:

transientStates(drunkMc)

## [1] "2" "3" "4"

The probability of either being absorbed by the bar or by the sofa at home are:

absorptionProbabilities(drunkMc)

##      1    5
## 2 0.75 0.25
## 3 0.50 0.50
## 4 0.25 0.75

which means that the probability of arriving home / bar is inversely proportional to the distance to each one.

But we also would like to know how much time does the person take to arrive there, which can be done with meanAbsorptionTime:

meanAbsorptionTime(drunkMc)

## 2 3 4 
## 3 4 3

Thus, the expected time to absorption is three steps from states 2 or 4 and four steps from state 3.

The drunkard’s walk is a special case of the classical gambler’s ruin chain : a walk on \(\{0,1,\ldots,n\}\) that is absorbed at both ends and, at every interior state, moves up with probability \(p\) and down with probability \(1-p\). gamblersRuin(upperBound, prob, states) builds exactly this chain, so it reproduces drunkMc above without having to write out the transition matrix by hand (\(n=4\), \(p=1/2\), relabelled to states "1", …, "5"):

ruinMc <- gamblersRuin(upperBound = 4, prob = 0.5, states = as.character(1:5))
identical(unclass(ruinMc@transitionMatrix), unclass(drunkMc@transitionMatrix))

## [1] TRUE

For a general (asymmetric) gambler’s ruin with \(p\ne 1/2\), the probability of ruin (absorption at \(0\)) starting from state \(i\) has the closed form \[ \Pr(\text{ruin} \mid i) = \frac{r^i-r^n}{1-r^n}, \qquad r=\frac{1-p}{p}, \] , which absorptionProbabilities reproduces numerically:

upperBound <- 6
prob <- 0.35
ruinAsymmetric <- gamblersRuin(upperBound = upperBound, prob = prob)
ap <- absorptionProbabilities(ruinAsymmetric)

r <- (1 - prob) / prob
i <- seq_len(upperBound - 1)
theoreticalRuin <- (r^i - r^upperBound) / (1 - r^upperBound)

cbind(state = i, package = as.numeric(ap[as.character(i), "0"]),
      formula = theoreticalRuin)

##      state   package   formula
## [1,]     1 0.9785860 0.9785860
## [2,]     2 0.9388171 0.9388171
## [3,]     3 0.8649606 0.8649606
## [4,]     4 0.7277986 0.7277986
## [5,]     5 0.4730691 0.4730691

fundamentalMatrix(drunkMc)

##     2 3   4
## 2 1.5 1 0.5
## 3 1.0 2 1.0
## 4 0.5 1 1.5

Q <- subchain(drunkMc, transientStates(drunkMc), method = "submatrix")
Q

##     2   3   4
## 2 0.0 0.5 0.0
## 3 0.5 0.0 0.5
## 4 0.0 0.5 0.0

rowSums(Q)

##   2   3   4 
## 0.5 1.0 0.5

solve(diag(nrow(Q)) - Q)

##     2 3   4
## 2 1.5 1 0.5
## 3 1.0 2 1.0
## 4 0.5 1 1.5

watchedDrunk <- subchain(drunkMc, transientStates(drunkMc), method = "renormalize")
watchedDrunk

## Unnamed Markov chain (subchain) 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  2, 3, 4 
##  The transition matrix  (by rows)  is defined as follows: 
##     2 3   4
## 2 0.0 1 0.0
## 3 0.5 0 0.5
## 4 0.0 1 0.0

rowSums(watchedDrunk@transitionMatrix)

## 2 3 4 
## 1 1 1

Committor probability

The committor probability is the probability of reaching a given set of states before another one. Suppose that we start on a cloudy day: the probability of experiencing a rainy day before a sunny one is 0.5, as the output below shows (the function takes the target set \(A\) first):

committorAB(mcWeather,3,1)

##  sunny cloudy   rain 
##    0.0    0.5    1.0

Hitting probabilities

For states \(i\) and \(j\), the hitting probability \(h_{ij}\) is the probability of ever visiting \(j\) when the chain starts in \(i\). For \(i = j\), the package uses the return probability after at least one transition, rather than the trivial time-zero probability. The values are the minimal non-negative solution of the first-step equations

\[\begin{equation} h_{ij}=p_{ij}+\sum_{k\ne j}p_{ik}h_{kj}. \end{equation}\]

hittingProbabilities first analyses the directed graph of positive transitions. This identifies probabilities that are structurally zero or one, without discarding very small positive transitions. This is important for nearly absorbing chains: a small transition probability may imply a long waiting time, but the probability of eventually taking that transition can still be exactly one. The remaining probabilities solve the linear system \((I-Q)h=r\) restricted to the undecided states. The argument solver chooses how: "direct", the default, is an LU factorisation, the most accurate choice; "bicgstab" is an iterative Krylov method working on the sparse matrix, the fastest choice on large sparse chains; "doubling" is the doubled Neumann series of earlier versions of the package, kept for reproducibility. The iterative solvers stop at the relative residual tol or after maxIter iterations, with a warning if tol is not reached.

M <- markovchain:::zeros(5)
M[1,1] <- M[5,5] <- 1
M[2,1] <- M[2,3] <- 1/2
M[3,2] <- M[3,4] <- 1/2
M[4,2] <- M[4,5] <- 1/2

hittingTest <- new("markovchain", transitionMatrix = M)
hittingProbabilities(hittingTest)

##     1     2     3         4   5
## 1 1.0 0.000 0.000 0.0000000 0.0
## 2 0.8 0.375 0.500 0.3333333 0.2
## 3 0.6 0.750 0.375 0.6666667 0.4
## 4 0.4 0.500 0.250 0.1666667 0.6
## 5 0.0 0.000 0.000 0.0000000 1.0

When only some target states matter, targets restricts the computation to the corresponding columns, which coincide with those of the full matrix. On this small chain the three solvers agree to rounding:

hittingProbabilities(hittingTest, targets = c("1", "5"))

##     1   5
## 1 1.0 0.0
## 2 0.8 0.2
## 3 0.6 0.4
## 4 0.4 0.6
## 5 0.0 1.0

all.equal(hittingProbabilities(hittingTest, solver = "direct"),
          hittingProbabilities(hittingTest, solver = "bicgstab"))

## [1] TRUE

For the irreducible finite chain mcWeather, every state is eventually visited with probability one:

hittingProbabilities(mcWeather)

##        sunny cloudy rain
## sunny      1      1    1
## cloudy     1      1    1
## rain       1      1    1

Stochastic monotonicity

A transition matrix \(P\) is stochastically monotone if the probability of moving to higher states is non-decreasing with respect to the starting state. In other words, row \(i\) is stochastically dominated by row \(j\) whenever \(i < j\). This property is fundamental in many applied models, such as credit ratings or social mobility, because it guarantees that the Markov chain preserves the initial ordering of states over time.

For instance, in an income mobility model, if you start in a higher income class, your probability of moving to or staying in a high income class should ideally be greater than or equal to someone starting in a lower income class.

The is.stochasticallyMonotone method verifies this property for any markovchain object:

# Define the state space for income classes
incomeStates <- c("low", "mid", "high")

# Example 1: A stochastically monotone transition matrix
# Higher starting income leads to higher probabilities of reaching/staying in higher brackets
monotoneMatrix <- matrix(c(0.6, 0.3, 0.1,
                           0.3, 0.5, 0.2,
                           0.1, 0.4, 0.5), 
                         byrow = TRUE, nrow = 3, 
                         dimnames = list(incomeStates, incomeStates))

mcIncomeMonotone <- new("markovchain", states = incomeStates, 
                        transitionMatrix = monotoneMatrix, 
                        name = "Monotone Income MC")

is.stochasticallyMonotone(mcIncomeMonotone)

## [1] TRUE

# Example 2: A non-monotone transition matrix
# Here, a person in the "mid" class has an unusually high chance of dropping 
# to "low" (0.7), which is worse than a person already in the "low" class (0.5).
# This violates the stochastic monotonicity property.
notMonotoneMatrix <- matrix(c(0.5, 0.3, 0.2,
                              0.7, 0.2, 0.1,
                              0.1, 0.4, 0.5), 
                            byrow = TRUE, nrow = 3, 
                            dimnames = list(incomeStates, incomeStates))

mcIncomeNotMonotone <- new("markovchain", states = incomeStates, 
                           transitionMatrix = notMonotoneMatrix, 
                           name = "Non-Monotone Income MC")

is.stochasticallyMonotone(mcIncomeNotMonotone)

## [1] FALSE

Classical Markov chain constructors

The package also offers ready-made constructors for a few textbook Markov chains that are frequently the starting point of an applied model, rather than requiring every state and transition of, e.g., a birth-death or diffusion chain to be typed in by hand as in Section @ref(sec:structure).

bd <- birthDeath(p = c(0.3, 0.4, 0.5), q = c(0.2, 0.3, 0.1))
bd

## Birth-Death Chain 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  1, 2, 3, 4 
##  The transition matrix  (by rows)  is defined as follows: 
##     1   2   3   4
## 1 0.7 0.3 0.0 0.0
## 2 0.2 0.4 0.4 0.0
## 3 0.0 0.3 0.2 0.5
## 4 0.0 0.0 0.1 0.9

reflected <- toBoundedChain(bd, "reflecting")
reflected@transitionMatrix[c(1, 4), ]

##   1 2 3 4
## 1 0 1 0 0
## 4 0 0 1 0

ehrenfest <- urnModel(balls = 6)
ehrenfest@transitionMatrix

##           0         1         2         3         4         5         6
## 0 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1 0.1666667 0.0000000 0.8333333 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.3333333 0.0000000 0.6666667 0.0000000 0.0000000 0.0000000
## 3 0.0000000 0.0000000 0.5000000 0.0000000 0.5000000 0.0000000 0.0000000
## 4 0.0000000 0.0000000 0.0000000 0.6666667 0.0000000 0.3333333 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.8333333 0.0000000 0.1666667
## 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000

piEhrenfest <- as.numeric(steadyStates(ehrenfest))
cbind(package = piEhrenfest, binomial = dbinom(0:6, 6, 0.5))

##       package binomial
## [1,] 0.015625 0.015625
## [2,] 0.093750 0.093750
## [3,] 0.234375 0.234375
## [4,] 0.312500 0.312500
## [5,] 0.234375 0.234375
## [6,] 0.093750 0.093750
## [7,] 0.015625 0.015625

drySeason <- new("markovchain", states = c("sunny", "rain"),
                 transitionMatrix = matrix(c(0.9, 0.1, 0.2, 0.8), byrow = TRUE,
                   nrow = 2, dimnames = list(c("sunny", "rain"), c("sunny", "rain"))),
                 name = "Dry season")
wetSeason <- new("markovchain", states = c("sunny", "rain"),
                 transitionMatrix = matrix(c(0.5, 0.5, 0.6, 0.4), byrow = TRUE,
                   nrow = 2, dimnames = list(c("sunny", "rain"), c("sunny", "rain"))),
                 name = "Wet season")

mergeWith(drySeason, wetSeason, gamma = 0.3)

## Dry season + Wet season (merged, gamma = 0.3) 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  sunny, rain 
##  The transition matrix  (by rows)  is defined as follows: 
##       sunny rain
## sunny  0.78 0.22
## rain   0.32 0.68

identical(unclass(toNthOrder(mcWeather, 5)@transitionMatrix),
          unclass((mcWeather ^ 5)@transitionMatrix))

## [1] TRUE

n <- 8
neutralMoran <- populationGeneticsModel(model = "moran", n = n, s = 0)
neutralWF <- populationGeneticsModel(model = "wright-fisher", n = n, s = 0)

i <- 1:(n - 1)
fixation <- function(mc) {
  B <- absorptionProbabilities(mc)
  as.numeric(B[as.character(i), as.character(n)])
}
cbind(i = i, moran = fixation(neutralMoran),
      wrightFisher = fixation(neutralWF), theoretical = i / n)

##      i moran wrightFisher theoretical
## [1,] 1 0.125        0.125       0.125
## [2,] 2 0.250        0.250       0.250
## [3,] 3 0.375        0.375       0.375
## [4,] 4 0.500        0.500       0.500
## [5,] 5 0.625        0.625       0.625
## [6,] 6 0.750        0.750       0.750
## [7,] 7 0.875        0.875       0.875

favouredMoran <- populationGeneticsModel(model = "moran", n = n, s = 0.5)
favouredWF <- populationGeneticsModel(model = "wright-fisher", n = n, s = 0.5)

# fixation probability starting from 4 copies out of 8, neutral vs favoured
c(neutral = as.numeric(absorptionProbabilities(neutralMoran)["4", "8"]),
  favoured = as.numeric(absorptionProbabilities(favouredMoran)["4", "8"]))

##   neutral  favoured 
## 0.5000000 0.8350515

fixed <- matrix(NA, 3, 3)
fixed[1, 1] <- 0.5
randomMarkovChain(3, zeros = 2, mask = fixed, seed = 1)@transitionMatrix

##           1         2          3
## 1 0.5000000 0.4572463 0.04275369
## 2 0.0000000 0.0000000 1.00000000
## 3 0.5397077 0.3594227 0.10086963

dirichletChain(4, diffusion = 1, diagonalBias = 2,
               seed = 1)@transitionMatrix

##           1          2            3           4
## 1 0.8533391 0.07032067 0.0752879483 0.001052291
## 2 0.2089573 0.75664931 0.0228492982 0.011544060
## 3 0.1796169 0.28299040 0.4809133706 0.056479302
## 4 0.3989103 0.14352165 0.0002538273 0.457314239

identityChain(states = c("a", "b"))@transitionMatrix

##   a b
## a 1 0
## b 0 1

Simulation

A random sequence from a DTMC is simulated with the function rmarkovchain. The following code generates 220 days of weather states from the process underlying mcWeather.

weathersOfDays <- rmarkovchain(n = 220, object = mcWeather, t0 = "sunny")
weathersOfDays[1:30]

##  [1] "cloudy" "rain"   "rain"   "rain"   "cloudy" "rain"   "rain"   "cloudy"
##  [9] "sunny"  "rain"   "cloudy" "cloudy" "cloudy" "rain"   "rain"   "rain"  
## [17] "cloudy" "rain"   "rain"   "rain"   "rain"   "cloudy" "cloudy" "sunny" 
## [25] "sunny"  "sunny"  "sunny"  "sunny"  "cloudy" "sunny"

Similarly, one or more sequences can be simulated from a semi-homogeneous Markov chain, as the following code (applied to the CCHC example) shows.

patientStates <- rmarkovchain(n = 5, object = mcCCRC, t0 = "H", 
                              include.t0 = TRUE)
patientStates[1:10,]

##    iteration values
## 1          1      H
## 2          1      H
## 3          1      D
## 4          1      D
## 5          1      D
## 6          2      H
## 7          2      H
## 8          2      D
## 9          2      D
## 10         2      D

Two advanced parameters of rmarkovchain select the implementation, among four options: , in parallel, and in parallel. The boolean parameters useRcpp and parallel choose the implementation; the default is and , i.e. the sequential implementation, which is generally faster than the one. On multicore processors, setting parallel = TRUE takes advantage of several cores; when both useRcpp = TRUE and parallel = TRUE, the parallelization relies on the package . Every simulated sequence draws from its own random number generator, seeded from ’s generator, so that set.seed() gives the same result whatever the number of threads. That number is set by the argument num.cores; when it is NULL, the default, the options RcppParallel.numThreads and Ncpus and the environment variables RCPP_PARALLEL_NUM_THREADS and OMP_NUM_THREADS are honoured, in this order, and otherwise at most two threads are used, as required by the CRAN policy.

Lumpability and aggregation of states

In some applications the original state space is too detailed, and it is useful to aggregate several micro-states into a smaller number of macro-states. A Markov chain is said to be strongly lumpable with respect to a partition of the state space if, for every pair of macro-states, the total probability of moving from a micro-state in the source block to the destination block is the same for all micro-states in the source block . Under this condition, the aggregated process is again Markovian and its transition matrix is well defined.

The function is.lumpable checks exact strong lumpability for a user-supplied partition. The function lump returns the aggregated Markov chain. By default, lump requires exact lumpability; setting force = TRUE performs an approximate weighted aggregation.

The following example is the classical Land of Oz weather chain. We aggregate rainy and snowy into a single Bad_Weather macro-state, leaving nice as a separate macro-state.

states_oz <- c("rainy", "nice", "snowy")

P_oz <- matrix(
  c(0.5,  0.25, 0.25,
    0.5,  0.00, 0.5,
    0.25, 0.25, 0.5),
  byrow = TRUE,
  nrow = 3,
  dimnames = list(states_oz, states_oz)
)

mcOz <- new(
  "markovchain",
  states = states_oz,
  transitionMatrix = P_oz,
  name = "Land of Oz"
)

partition_oz <- list(
  Bad_Weather = c("rainy", "snowy"),
  Nice_Weather = c("nice")
)

is.lumpable(mcOz, partition_oz)

## [1] TRUE

lumped_oz <- lump(mcOz, partition_oz)
lumped_oz

## Land of Oz (Lumped) 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  Bad_Weather, Nice_Weather 
##  The transition matrix  (by rows)  is defined as follows: 
##              Bad_Weather Nice_Weather
## Bad_Weather         0.75         0.25
## Nice_Weather        1.00         0.00

The lumped chain has two states. From either rainy or snowy, the probability of moving to nice is 0.25, so the transition probability from Bad_Weather to Nice_Weather is also 0.25. The remaining probability, 0.75, stays inside the Bad_Weather macro-state.

If the exact lumpability condition is not satisfied, lump stops by default:

P_bad <- matrix(
  c(0.5, 0.3, 0.2,
    0.7, 0.2, 0.1,
    0.1, 0.4, 0.5),
  byrow = TRUE,
  nrow = 3,
  dimnames = list(c("A1", "A2", "B"), c("A1", "A2", "B"))
)

mcBad <- new(
  "markovchain",
  states = c("A1", "A2", "B"),
  transitionMatrix = P_bad
)

partition_bad <- list(
  Group_A = c("A1", "A2"),
  Group_B = c("B")
)

is.lumpable(mcBad, partition_bad)

## [1] FALSE

try(lump(mcBad, partition_bad))

## Error in lump(mcBad, partition_bad) : 
##   The Markov chain is not exactly lumpable. Use force = TRUE to perform an approximate weighted lumping.

For exploratory analyses, approximate lumping can be forced:

lump(mcBad, partition_bad, force = TRUE)

## Unnamed Markov chain (Lumped) 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  Group_A, Group_B 
##  The transition matrix  (by rows)  is defined as follows: 
##           Group_A   Group_B
## Group_A 0.8389831 0.1610169
## Group_B 0.5000000 0.5000000

When force = TRUE, the aggregated transition probabilities are computed using stationary weights when available. If a macro-state has zero stationary mass, an arithmetic average over its micro-states is used as a fallback. This produces a valid coarse transition matrix, but it should be interpreted as an approximation: the aggregated process is not guaranteed to be exactly Markovian unless the strong lumpability condition holds.

The helper autoLump provides a heuristic way to discover macro-states using spectral clustering. It should be regarded as an exploratory approximation, not as a test of exact lumpability.

rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")

creditMatrix <- matrix(
  c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01,
    0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01,
    0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06,
    0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18,
    0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06,
    0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20,
    0.21, 0.00, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79,
    0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100.00),
  nrow = 8,
  byrow = TRUE,
  dimnames = list(rc, rc)
) / 100

creditMc <- new(
  "markovchain",
  states = rc,
  transitionMatrix = creditMatrix,
  name = "Credit ratings"
)

set.seed(123)
auto_lumped <- autoLump(creditMc, k = 3)

auto_lumped$partition
auto_lumped$lumped_chain

To keep the vignette fast, the example uses a small transition matrix. In larger applications, the automatically discovered partition should be inspected and validated substantively before being used for modelling or forecasting.

agg_oz <- aggregateStates(mcOz, k = 2)
agg_oz$partition

## $Macro_1
## [1] "rainy" "snowy"
## 
## $Macro_2
## [1] "nice"

isTRUE(all.equal(unname(agg_oz$aggregatedChain@transitionMatrix),
                  unname(lumped_oz@transitionMatrix)))

## [1] TRUE

agg_oz$klDivergence

## [1] 0.0490225

states_dmm <- as.character(1:4)
P_dmm <- matrix(
  c(0.5, 0.4, 0.0, 0.1,
    0.4, 0.5, 0.1, 0.0,
    0.0, 0.1, 0.5, 0.4,
    0.1, 0.0, 0.4, 0.5),
  byrow = TRUE, nrow = 4,
  dimnames = list(states_dmm, states_dmm)
)
mcDMM <- new("markovchain", states = states_dmm, transitionMatrix = P_dmm,
             name = "Deng-Mehta-Meyn (2009) example")

agg_dmm <- aggregateStates(mcDMM, k = 2)
agg_dmm$partition

## $Macro_1
## [1] "1" "2"
## 
## $Macro_2
## [1] "3" "4"

agg_dmm$aggregatedChain@transitionMatrix

##         Macro_1 Macro_2
## Macro_1     0.9     0.1
## Macro_2     0.1     0.9

agg_dmm$klDivergence

## [1] 0.1080315

Serialization and export

A markovchain object is often built once and then reused elsewhere: shared with a colleague, stored alongside a dataset, or handed to a program written in another language. toDictionary and fromDictionary, together with toFile and fromFile, give the object a plain, language-neutral representation for exactly this purpose, mirroring PyDTMC’s own to_dictionary, from_dictionary, to_file and from_file – with one deliberate difference: PyDTMC represents a chain as a flat mapping from every (from_state, to_state) pair to its probability (\(n^2\) entries with no state grouping), whereas here the representation is nested by source state, which is both more compact to read and round-trips through R’s own list-of-lists idiom without any special tuple-key handling.

We reuse mcOz, the Land of Oz chain introduced in the previous section, a well-known textbook example : unlike a synthetic chain built only to demonstrate file input and output, it is worth saving and reloading.

oz_dict <- toDictionary(mcOz)
str(oz_dict)

## List of 4
##  $ name            : chr "Land of Oz"
##  $ states          : chr [1:3] "rainy" "nice" "snowy"
##  $ byrow           : logi TRUE
##  $ transitionMatrix:List of 3
##   ..$ rainy:List of 3
##   .. ..$ rainy: num 0.5
##   .. ..$ nice : num 0.25
##   .. ..$ snowy: num 0.25
##   ..$ nice :List of 3
##   .. ..$ rainy: num 0.5
##   .. ..$ nice : num 0
##   .. ..$ snowy: num 0.5
##   ..$ snowy:List of 3
##   .. ..$ rainy: num 0.25
##   .. ..$ nice : num 0.25
##   .. ..$ snowy: num 0.5

toDictionary always normalizes to a row-stochastic representation, so it faithfully captures mcOz regardless of the chain’s own byrow storage convention. Reversing it exactly recovers the original transition matrix:

identical(fromDictionary(oz_dict)@transitionMatrix, mcOz@transitionMatrix)

## [1] TRUE

toFile/fromFile write this same representation to disk, inferring the format (JSON, YAML or CSV) from the file extension:

oz_json <- tempfile(fileext = ".json")
toFile(mcOz, oz_json)
cat(readLines(oz_json), sep = "\n")

## {
##   "name": "Land of Oz",
##   "states": ["rainy", "nice", "snowy"],
##   "byrow": true,
##   "transitionMatrix": {
##     "rainy": {
##       "rainy": 0.5,
##       "nice": 0.25,
##       "snowy": 0.25
##     },
##     "nice": {
##       "rainy": 0.5,
##       "nice": 0,
##       "snowy": 0.5
##     },
##     "snowy": {
##       "rainy": 0.25,
##       "nice": 0.25,
##       "snowy": 0.5
##     }
##   }
## }

identical(fromFile(oz_json)@transitionMatrix, mcOz@transitionMatrix)

## [1] TRUE

JSON and YAML preserve every field, including the chain’s name, at full numeric precision. CSV, which has no natural place for a name, only stores the transition matrix itself – the same limitation PyDTMC’s own CSV format has:

oz_csv <- tempfile(fileext = ".csv")
toFile(mcOz, oz_csv)
cat(readLines(oz_csv), sep = "\n")

## "","rainy","nice","snowy"
## "rainy",0.5,0.25,0.25
## "nice",0.5,0,0.5
## "snowy",0.25,0.25,0.5

fromFile(oz_csv)@name # the package's own default, not "Land of Oz"

## [1] "Unnamed Markov chain"

The XML format is PyDTMC’s own, so a chain can travel between the two packages in both directions; the name is kept as an attribute that PyDTMC ignores. Reading XML requires the package.

oz_xml <- tempfile(fileext = ".xml")
toFile(mcOz, oz_xml)
cat(readLines(oz_xml)[1:4], sep = "\n")

## <?xml version='1.0' encoding='utf-8' standalone='yes' ?>
## <MarkovChain name="Land of Oz">
##  <Item state_from="rainy" state_to="rainy" probability="0.5"/>
##  <Item state_from="rainy" state_to="nice" probability="0.25"/>

Finally, toDot and toMermaid export the transition graph as Graphviz DOT or Mermaid text, to be rendered outside or pasted into a Markdown document; transitions below minProbability can be left out.

cat(toMermaid(mcOz, digits = 2), "\n")

## flowchart LR
##   s1(("rainy"))
##   s2(("nice"))
##   s3(("snowy"))
##   s1 -->|"0.5"| s1
##   s1 -->|"0.25"| s2
##   s1 -->|"0.25"| s3
##   s2 -->|"0.5"| s1
##   s2 -->|"0.5"| s3
##   s3 -->|"0.25"| s1
##   s3 -->|"0.25"| s2
##   s3 -->|"0.5"| s3

Estimation

A time-homogeneous Markov chain can be fitted from observed state sequences. The current version of implements four fitting methods: maximum likelihood (MLE), MLE with Laplace smoothing, bootstrap estimation and maximum a posteriori (MAP) estimation.

Equation \(\ref{eq:MLE}\) shows the maximum likelihood estimator of the transition probability \(p_{ij}\). Here \(n_{ij}\) is the number of observed transitions \(X_t=s_i, X_{t+1}=s_j\) and \(n_i = \sum_{u=1}^k n_{iu}\) is the number of observed outgoing transitions from state \(s_i\):

\[\begin{equation} {\hat p^{MLE}}_{ij} = \frac{n_{ij}}{n_i}. \label{eq:MLE} \end{equation}\]

For a fixed origin state \(s_i\), the counts \((n_{i1}, \ldots, n_{ik})\) are treated as a multinomial sample with row total \(n_i = \sum_{u=1}^k n_{iu}\). The standardError returned for each marginal transition probability is therefore the usual binomial standard error in Equation @ref(eq:SE).

\[\begin{equation} SE_{ij} = \sqrt{\frac{{\hat p^{MLE}}_{ij}(1-{\hat p^{MLE}}_{ij})}{n_i}} \label{eq:SE} \end{equation}\]

weatherFittedMLE <- markovchainFit(data = weathersOfDays, method = "mle",name = "Weather MLE")
weatherFittedMLE$estimate

## Weather MLE 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##           cloudy       rain     sunny
## cloudy 0.3947368 0.35526316 0.2500000
## rain   0.4909091 0.38181818 0.1272727
## sunny  0.2045455 0.07954545 0.7159091

weatherFittedMLE$standardError

##            cloudy       rain      sunny
## cloudy 0.05606853 0.05489838 0.04966996
## rain   0.06740884 0.06550962 0.04493923
## sunny  0.04299929 0.02884480 0.04807467

The Laplace smoothing approach is a smoothed version of the MLE in which each count \(n_{ij}\) is replaced by \(n_{ij}+\alpha\) (Equation \(\ref{eq:LAPLACE}\)). The smoothing parameter \(\alpha\) must be non-negative; \(\alpha=0\) gives the unsmoothed MLE.

\[\begin{equation} {\hat p^{LS}}_{ij} = \frac{{{n_{ij}} + \alpha }}{{\sum\limits_{u = 1}^k {\left( {{n_{iu}} + \alpha } \right)} }} \label{eq:LAPLACE} \end{equation}\]

weatherFittedLAPLACE <- markovchainFit(data = weathersOfDays, 
                                    method = "laplace", laplacian = 0.01,
                                    name = "Weather LAPLACE")
weatherFittedLAPLACE$estimate

## Weather LAPLACE 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##           cloudy       rain     sunny
## cloudy 0.3947126 0.35525450 0.2500329
## rain   0.4908232 0.38179175 0.1273851
## sunny  0.2045893 0.07963194 0.7157787

The confidence interval option is enabled by default; set confint = FALSE to reduce computation time when intervals are not needed. Both MLE and Laplace fitting are based on createSequenceMatrix, which returns the raw transition-count matrix.

createSequenceMatrix(stringchar = weathersOfDays)

##        cloudy rain sunny
## cloudy     30   27    19
## rain       27   21     7
## sunny      18    7    63

stringchar could contain NA values, and the transitions containing NA would be ignored.

If a state is observed only as the last element of a sequence, or is supplied through the argument possibleStates but has no observed outgoing transition, its row has zero total count. In this case the point estimate follows the package convention of assigning a uniform row, \(p_{\beta j}=1/k\). For MLE confidence intervals, such a row carries no empirical information, so the corresponding cell-wise intervals are uninformative.

The uniform row is an assumption about the data rather than a consequence of them. The argument sanitize controls it: TRUE (the default of markovchainFit), or equivalently "uniform", gives the uniform row; "absorbing" makes every such state absorbing instead, adding no transition that was never observed; FALSE leaves the row at zero, which gives a matrix that is not stochastic.

partial <- c("a", "b", "a", "b", "a")
markovchainFit(partial, possibleStates = c("a", "b", "c"),
               sanitize = "uniform")$estimate@transitionMatrix

##           a         b         c
## a 0.0000000 1.0000000 0.0000000
## b 1.0000000 0.0000000 0.0000000
## c 0.3333333 0.3333333 0.3333333

markovchainFit(partial, possibleStates = c("a", "b", "c"),
               sanitize = "absorbing")$estimate@transitionMatrix

##   a b c
## a 0 1 0
## b 1 0 0
## c 0 0 1

createSequenceMatrix can also be used to obtain raw transition counts from a given \(n \times 2\) matrix of (from, to) pairs, as the following example shows:

myMatr<-matrix(c("a","b","b","a","a","b","b","b","b","a","a","a","b","a"),ncol=2)
createSequenceMatrix(stringchar = myMatr,toRowProbs = TRUE)

##           a         b
## a 0.6666667 0.3333333
## b 0.5000000 0.5000000

A bootstrap estimation approach is also available to assess the variability of \(\hat p_{ij}\). The implementation in proceeds as follows:

1. generate bootstrap state sequences using the transition probabilities estimated from the original data. The number of bootstrap samples is controlled by the nboot argument of markovchainFit;
2. apply MLE estimation to each bootstrap sequence; the bootstrap estimates are saved in the bootStrapSamples component of the returned list;
3. compute the bootstrap transition estimate as the row-normalized average of the bootstrap MLE estimates. The reported standardError is the empirical standard deviation of the bootstrap estimates for each transition probability.

weatherFittedBOOT <- markovchainFit(data = weathersOfDays, 
                                    method = "bootstrap", nboot = 10)
weatherFittedBOOT$estimate

## BootStrap Estimate 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##           cloudy       rain     sunny
## cloudy 0.4039584 0.34401538 0.2520262
## rain   0.4619274 0.42383613 0.1142365
## sunny  0.2139972 0.09615059 0.6898522

weatherFittedBOOT$standardError

##            cloudy       rain      sunny
## cloudy 0.07466146 0.06249693 0.05123581
## rain   0.03943952 0.06057494 0.04185424
## sunny  0.05526397 0.02727959 0.04988424

The bootstrapping process can be done in parallel thanks to package . The parallelized implementation is recommended when the data sample size or the required number of bootstrap runs is high.

weatherFittedBOOTParallel <- markovchainFit(data = weathersOfDays, 
                                    method = "bootstrap", nboot = 200, 
                                    parallel = TRUE)
weatherFittedBOOTParallel$estimate
weatherFittedBOOTParallel$standardError

By default the parallel bootstrapping uses at most two threads, as required by the CRAN policy, unless the options or environment variables listed above for rmarkovchain ask for more. The argument num.cores sets the number of threads explicitly; for example, the following code uses four threads. With the same seed, the bootstrap samples are the same whatever the number of threads.

markovchainFit(data = weathersOfDays, method = "bootstrap", nboot = 200,
               parallel = TRUE, num.cores = 4)

For more details, please refer to web site.

For all the fitting methods, the logLikelihood denoted in Equation \(\ref{eq:LLH}\) is provided.

\[\begin{equation} LLH = \sum_{i,j} n_{ij} * log (p_{ij}) \label{eq:LLH} \end{equation}\] where \(n_{ij}\) is the entry of the frequency matrix and \(p_{ij}\) is the entry of the transition probability matrix.

weatherFittedMLE$logLikelihood

## [1] -203.3674

weatherFittedBOOT$logLikelihood

## [1] -203.7924

Confidence matrices of estimated parameters are available as well. The confidenceInterval component contains two matrices, lowerEndpointMatrix and upperEndpointMatrix. The confidence level is 0.95 by default and can be changed through the confidencelevel argument of markovchainFit.

For MLE fits, reports cell-wise Wilson score intervals for the marginal transition probabilities . Compared with the simple Wald interval \(\hat p_{ij} \pm z SE_{ij}\), Wilson intervals have better small-sample behaviour and do not collapse to a zero-width interval when a transition is never observed or always observed. With \(z = \Phi^{-1}(1 - \alpha/2)\), \(n_i = \sum_u n_{iu}\), \(\hat p_{ij}=n_{ij}/n_i\), and confidence level \(1-\alpha\), the interval is

\[\begin{align} LowerEndpoint_{ij} &= \frac{\hat p_{ij} + z^2/(2n_i) - z\sqrt{\hat p_{ij}(1-\hat p_{ij})/n_i + z^2/(4n_i^2)}}{1+z^2/n_i}, \label{eq:CIL} \\ UpperEndpoint_{ij} &= \frac{\hat p_{ij} + z^2/(2n_i) + z\sqrt{\hat p_{ij}(1-\hat p_{ij})/n_i + z^2/(4n_i^2)}}{1+z^2/n_i}. \label{eq:CIU} \end{align}\]

If no outgoing transition is available for a state, the point estimate follows the package convention of a uniform row, while the interval is left maximally uninformative. These intervals are cell-wise marginal intervals, not simultaneous multinomial confidence regions for a full row or matrix. For simultaneous multinomial intervals, use multinomialConfidenceIntervals. For bootstrap fits, intervals use the empirical standard deviation of the bootstrap transition-probability estimates and the two-sided normal quantile corresponding to confidencelevel. For MAP fits, credible intervals are equal-tailed beta posterior intervals, obtained from the inverse beta distribution .

weatherFittedMLE$confidenceInterval

## NULL

weatherFittedBOOT$confidenceInterval

## $confidenceLevel
## [1] 0.95
## 
## $lowerEndpointMatrix
##           cloudy       rain      sunny
## cloudy 0.2576246 0.22152365 0.15160591
## rain   0.3846273 0.30511143 0.03220371
## sunny  0.1056818 0.04268359 0.59208089
## 
## $upperEndpointMatrix
##           cloudy      rain     sunny
## cloudy 0.5502921 0.4665071 0.3524466
## rain   0.5392274 0.5425608 0.1962693
## sunny  0.3223126 0.1496176 0.7876235

The function multinomialConfidenceIntervals computes simultaneous multinomial confidence intervals. It is an translation of the functions of the package , which are in turn based on .

multinomialConfidenceIntervals(transitionMatrix = 
        weatherFittedMLE$estimate@transitionMatrix, 
        countsTransitionMatrix = createSequenceMatrix(weathersOfDays))

## $confidenceLevel
## [1] 0.95
## 
## $lowerEndpointMatrix
##           cloudy      rain     sunny
## cloudy 0.2894737 0.2500000 0.1447368
## rain   0.3636364 0.2545455 0.0000000
## sunny  0.1250000 0.0000000 0.6363636
## 
## $upperEndpointMatrix
##           cloudy      rain     sunny
## cloudy 0.5262267 0.4867530 0.3814898
## rain   0.6329951 0.5239042 0.2693588
## sunny  0.3049268 0.1799268 0.8162905

Most of the functions for fitting DTMC have been rewritten in using since version 0.2.

A DTMC can also be fitted from matrix or data.frame objects, as shown in the following code.

data(holson)
singleMc<-markovchainFit(data=holson[,2:12],name="holson")

The same applies for markovchainList (output length has been limited).

mcListFit<-markovchainListFit(data=holson[,2:6],name="holson")
mcListFit$estimate

## holson  list of Markov chain(s) 
## Markovchain  1 
## Unnamed Markov chain 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  1, 2, 3 
##  The transition matrix  (by rows)  is defined as follows: 
##            1          2         3
## 1 0.94609164 0.05390836 0.0000000
## 2 0.26356589 0.62790698 0.1085271
## 3 0.02325581 0.18604651 0.7906977
## 
## Markovchain  2 
## Unnamed Markov chain 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  1, 2, 3 
##  The transition matrix  (by rows)  is defined as follows: 
##           1         2         3
## 1 0.9323410 0.0676590 0.0000000
## 2 0.2551724 0.5103448 0.2344828
## 3 0.0000000 0.0862069 0.9137931
...

Finally, given a list object, it is possible to fit a markovchain object or to obtain the raw transition matrix.

c1<-c("a","b","a","a","c","c","a")
c2<-c("b")
c3<-c("c","a","a","c")
c4<-c("b","a","b","a","a","c","b")
c5<-c("a","a","c",NA)
c6<-c("b","c","b","c","a")
mylist<-list(c1,c2,c3,c4,c5,c6)
mylistMc<-markovchainFit(data=mylist)
mylistMc

## $estimate
## MLE Fit 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c 
##  The transition matrix  (by rows)  is defined as follows: 
##     a         b         c
## a 0.4 0.2000000 0.4000000
## b 0.6 0.0000000 0.4000000
## c 0.5 0.3333333 0.1666667
## 
## 
## $standardError
##           a         b         c
## a 0.1549193 0.1264911 0.1549193
## b 0.2190890 0.0000000 0.2190890
## c 0.2041241 0.1924501 0.1521452
## 
## $confidenceLevel
## [1] 0.95
## 
## $lowerEndpointMatrix
##           a          b          c
## a 0.1681803 0.05668215 0.16818033
## b 0.2307243 0.00000000 0.11762077
## c 0.1876163 0.09677141 0.03005337
## 
## $upperEndpointMatrix
##           a         b         c
## a 0.6873262 0.5098375 0.6873262
## b 0.8823792 0.4344825 0.7692757
## c 0.8123837 0.7000067 0.5635028

The same works for markovchainFitList.

markovchainListFit(data=mylist)

## $estimate
##   list of Markov chain(s) 
## Markovchain  1 
## Unnamed Markov chain 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c 
##  The transition matrix  (by rows)  is defined as follows: 
##     a   b   c
## a 0.5 0.5 0.0
## b 0.5 0.0 0.5
## c 1.0 0.0 0.0
## 
## Markovchain  2 
## Unnamed Markov chain 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
...

If any transition contains NA, it will be ignored in the results as the above example showed.

Prediction

The \(n\)-step forward predictions can be obtained using the predict methods explicitly written for markovchain and markovchainList objects. The prediction is the mode of the conditional distribution of \(X_{t+1}\) given \(X_{t}=s_{j}\), being \(s_{j}\) the last realization of the DTMC (homogeneous or semi-homogeneous).

predict(object = weatherFittedMLE$estimate, newdata = c("cloudy", "sunny"),
        n.ahead = 3)

## [1] "sunny" "sunny" "sunny"

predict(mcCCRC, newdata = c("H", "H"), n.ahead = 5)

## [1] "H" "D" "D"

predict(mcCCRC, newdata = c("H", "H"), n.ahead = 5, continue = TRUE)

## [1] "H" "D" "D" "D" "D"

Statistical tests

In this section, we describe the statistical tests: assessing the independence of consecutive states (assessIndependence), the Markov property (verifyMarkovProperty), the order (assessOrder), the time-homogeneity (assessStationarity) of a Markov chain sequence, and the divergence tests (verifyEmpiricalToTheoretical, verifyHomogeneity) for empirically estimated transition matrices. Most of these tests are based on the \(\chi ^2\) statistics; the Pearson and likelihood-ratio (G) versions of each test are provided, together with a Monte Carlo option (method = "simulation") that avoids the asymptotic \(\chi^2\) approximation altogether by simulating the null distribution of the G statistic directly from the fitted model. Test results are returned as objects of class htest, following R’s standard convention for hypothesis tests, so that they print, and can be manipulated, like the output of stats::chisq.test. Relevant references are and .

As the sample size increases, even slight departures from the Markov property are easily detected, so these tests are most informative on samples of small or moderate size. The Monte Carlo option is particularly useful in this regime, since the asymptotic \(\chi^2\) approximation can be unreliable when expected counts are small; the simulated \(p\)-value is computed as \((1+\#\{G^*\ge G\})/(B+1)\) over \(B\) parametric replicates drawn under the null hypothesis, which is a finite-sample-valid Monte Carlo test . The accuracy of the statistical inference functions for very small samples is nonetheless an area of ongoing work.

data(rain)
assessIndependence(rain$rain)

## 
##  Pearson's Chi-squared test for independence of consecutive states
## 
## data:  rain$rain
## X-squared = 190.61, df = 4, p-value < 2.2e-16

data(preproglucacon)
assessIndependence(preproglucacon$preproglucacon, method = "G")

## 
##  Likelihood-ratio test for independence of consecutive states
## 
## data:  preproglucacon$preproglucacon
## G-squared = 64.455, df = 9, p-value = 1.841e-10

chisq.test(createSequenceMatrix(rain$rain), correct = FALSE)$statistic

## X-squared 
##  190.6124

koeberg <- as.character(read.csv(system.file("extdata", "koeberg_wind.csv",
                                             package = "markovchain"))$state)
assessIndependence(koeberg[14:length(koeberg)], method = "G")

## 
##  Likelihood-ratio test for independence of consecutive states
## 
## data:  koeberg[14:length(koeberg)]
## G-squared = 1083.1, df = 9, p-value < 2.2e-16

sample_sequence <- c("a", "b", "a", "a", "a", "a", "b", "a", "b", "a", 
                     "b", "a", "a", "b", "b", "b", "a")
verifyMarkovProperty(sample_sequence, method = "G")

## 
##  Likelihood-ratio test
## 
## data:  sample_sequence
## G-squared = 0.59916, df = 2, p-value = 0.7411

burkett_counts <- matrix(
  c(127, 59,
     11,  7,
     36, 20,
      9,  4),
  nrow = 4,
  byrow = TRUE,
  dimnames = list(
    c("previous state 1", "previous state 2",
      "previous state 4", "previous state 5"),
    c("transition to state 1", "transition elsewhere")
  )
)
chisq.test(burkett_counts, correct = FALSE)

## Warning in chisq.test(burkett_counts, correct = FALSE): Chi-squared
## approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  burkett_counts
## X-squared = 0.63611, df = 3, p-value = 0.8881

data(rain)
assessOrder(rain$rain)

## 
##  Pearson's Chi-squared test for Markov order
## 
## data:  rain$rain
## X-squared = 26.096, df = 12, p-value = 0.0104

assessStationarity(rain$rain, 10)

## 
##  Pearson's Chi-squared test for time-homogeneity
## 
## data:  rain$rain
## X-squared = 110.02, df = 54, p-value = 1.05e-05

structural.zeros <- matrix(FALSE, 3, 3,
                           dimnames = list(c("a", "b", "c"),
                                           c("a", "b", "c")))
structural.zeros["a", "a"] <- TRUE

structural_sequence <- c(
  "a", "b", "a", "c", "a", "b", "c", "a",
  "b", "c", "b", "a", "c", "b", "a", "c",
  "a", "b", "c", "a", "b", "c", "b", "a"
)

assessStationarity(structural_sequence, nblocks = 4,
                   structural.zeros = structural.zeros)

## 
##  Pearson's Chi-squared test for time-homogeneity
## 
## data:  structural_sequence
## X-squared = 1.7639, df = 15, p-value = 1

kullback_counts <- matrix(
  c(51, 11,  8,
    12, 31,  9,
     6, 11, 10),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

kullback_P <- matrix(
  c(5/8, 1/4, 1/8,
    1/4, 1/2, 1/4,
    1/4, 3/8, 3/8),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

kullback_mc <- new("markovchain", states = c("a", "b", "c"),
                   transitionMatrix = kullback_P)

verifyEmpiricalToTheoretical(kullback_counts, kullback_mc, method = "G")

## 
##  Likelihood-ratio test
## 
## data:  kullback_counts
## G-squared = 6.5184, df = 6, p-value = 0.3677

homogeneity_P <- matrix(
  c(.70, .20, .10,
    .20, .50, .30,
    .10, .30, .60),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

homogeneity_mc <- new("markovchain", states = c("a", "b", "c"),
                      transitionMatrix = homogeneity_P)

set.seed(123)
x1 <- rmarkovchain(500, homogeneity_mc, t0 = "a")
x2 <- rmarkovchain(500, homogeneity_mc, t0 = "a")

verifyHomogeneity(list(x1, x2), method = "G")

## 
##  Likelihood-ratio test
## 
## data:  list(x1, x2)
## G-squared = 1.8881, df = 6, p-value = 0.9297

Pseudo-Bayesian estimation

show an empirical quasi-Bayesian method to estimate transition matrices, given an empirical \(\hat{P}\) transition matrix (estimated using the classical approach) and an a priori estimate \(Q\). In particular, each row of the matrix is estimated using the linear combination \(\alpha \cdot Q+\left(1-\alpha\right) \cdot P\), where \(\alpha\) is defined for each row as Equation \(\ref{eq:pseudobayes}\) shows

\[\begin{equation} \left\{\begin{matrix} \hat{\alpha_i}=\frac{\hat{K_i}}{v\left(i \right )+\hat{K_i}}\\ \hat{K_i}=\frac{v\left(i \right)^2 - \sum_{j}Y_{ij}^2}{\sum_{j}(Y_{ij}-v\left(i \right)*q_{ij})^2} \end{matrix}\right. \label{eq:pseudobayes} \end{equation}\]

The following code returns the pseudo Bayesian estimate of the transition matrix:

pseudoBayesEstimator <- function(raw, apriori){
  v_i <- rowSums(raw) 
  K_i <- numeric(nrow(raw))
  sumSquaredY <- rowSums(raw^2)
  #get numerator
  K_i_num <- v_i^2-sumSquaredY
  #get denominator
  VQ <- matrix(0,nrow= nrow(apriori),ncol=ncol(apriori))
  for (i in 1:nrow(VQ)) {
    VQ[i,]<-v_i[i]*apriori[i,]
  }
  
  K_i_den<-rowSums((raw - VQ)^2)
  
  K_i <- K_i_num/K_i_den
  
  #get the alpha vector
  alpha <- K_i / (v_i+K_i)
  
  #empirical transition matrix
  Emp<-raw/rowSums(raw)
  
  #get the estimate
  out<-matrix(0, nrow= nrow(raw),ncol=ncol(raw))
  for (i in 1:nrow(out)) {
    out[i,]<-alpha[i]*apriori[i,]+(1-alpha[i])*Emp[i,]
  }
  return(out)
}

We then apply it to the weather example:

trueMc<-as(matrix(c(0.1, .9,.7,.3),nrow = 2, byrow = 2),"markovchain")
aprioriMc<-as(matrix(c(0.5, .5,.5,.5),nrow = 2, byrow = 2),"markovchain")

smallSample<-rmarkovchain(n=20,object = trueMc)
smallSampleRawTransitions<-createSequenceMatrix(stringchar = smallSample)
pseudoBayesEstimator(
  raw = smallSampleRawTransitions, 
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

##            s1         s2
## s1 -0.1000000  0.1000000
## s2  0.0806874 -0.0806874

biggerSample<-rmarkovchain(n=60,object = trueMc)
biggerSampleRawTransitions<-createSequenceMatrix(stringchar = biggerSample)
pseudoBayesEstimator(
  raw = biggerSampleRawTransitions,
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

##             s1         s2
## s1 -0.10000000 0.10000000
## s2 -0.09574646 0.09574646

bigSample<-rmarkovchain(n=100,object = trueMc)
bigSampleRawTransitions<-createSequenceMatrix(stringchar = bigSample)
pseudoBayesEstimator(
  raw = bigSampleRawTransitions,
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

##              s1          s2
## s1 -0.026752827 0.026752827
## s2 -0.001470588 0.001470588

Bayesian estimation

The package provides functionality for maximum a posteriori (MAP) estimation of the chain parameters (at the time of writing this document, only first order models are supported) by Bayesian inference. It also computes the probability of observing a new data set, given a (different) data set. This section describes the mathematics of the methods employed by the package.

weatherStates <- c("sunny", "cloudy", "rain")
byRow <- TRUE
weatherMatrix <- matrix(data = c(0.7, 0.2, 0.1, 
                                 0.3, 0.4, 0.3, 
                                 0.2, 0.4, 0.4), 
                        byrow = byRow, nrow = 3, 
                        dimnames = list(weatherStates, weatherStates))
mcWeather <- new("markovchain", states = weatherStates, 
                 byrow = byRow, transitionMatrix = weatherMatrix, 
                 name = "Weather")      
weathersOfDays <- rmarkovchain(n = 365, object = mcWeather, t0 = "sunny")

hyperMatrix<-matrix(c(1, 1, 2, 
                      3, 2, 1,
                      2, 2, 3), 
                    nrow = 3, byrow = TRUE,
                    dimnames = list(weatherStates,weatherStates))
markovchainFit(weathersOfDays[1:200], method = "map", 
               confidencelevel = 0.92, hyperparam = hyperMatrix)

## $estimate
## Bayesian Fit 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##           cloudy       rain     sunny
## cloudy 0.2923077 0.35384615 0.3538462
## rain   0.5365854 0.31707317 0.1463415
## sunny  0.2277228 0.03960396 0.7326733
## 
## 
## $expectedValue
##           cloudy       rain     sunny
## cloudy 0.2941176 0.35294118 0.3529412
## rain   0.5227273 0.31818182 0.1590909
## sunny  0.2307692 0.04807692 0.7211538
## 
## $standardError
##            cloudy       rain      sunny
## cloudy 0.05485322 0.05753055 0.05753055
## rain   0.07445856 0.06943297 0.05452441
## sunny  0.04111711 0.02087731 0.04376245
## 
## $confidenceInterval
## $confidenceInterval$confidenceLevel
## [1] 0.92
## 
## $confidenceInterval$lowerEndpointMatrix
##           cloudy       rain      sunny
## cloudy 0.2024315 0.25531059 0.25531059
## rain   0.3917409 0.20261227 0.07519196
## sunny  0.1624958 0.01813559 0.64169687
## 
## $confidenceInterval$upperEndpointMatrix
##           cloudy      rain     sunny
## cloudy 0.3941775 0.4565555 0.4565555
## rain   0.6522800 0.4452100 0.2643798
## sunny  0.3061880 0.0898865 0.7947391
## 
## 
## $logLikelihood
## [1] -169.9034

predictiveDistribution(weathersOfDays[1:200], 
                       weathersOfDays[201:365],hyperparam = hyperMatrix) 

## [1] -156.7419

hyperMatrix2<- hyperMatrix[c(2,3,1), c(2,3,1)]
markovchainFit(weathersOfDays[1:200], method = "map", 
               confidencelevel = 0.92, hyperparam = hyperMatrix2)

## $estimate
## Bayesian Fit 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  cloudy, rain, sunny 
##  The transition matrix  (by rows)  is defined as follows: 
##           cloudy       rain     sunny
## cloudy 0.2923077 0.35384615 0.3538462
## rain   0.5365854 0.31707317 0.1463415
## sunny  0.2277228 0.03960396 0.7326733
## 
## 
## $expectedValue
##           cloudy       rain     sunny
## cloudy 0.2941176 0.35294118 0.3529412
## rain   0.5227273 0.31818182 0.1590909
## sunny  0.2307692 0.04807692 0.7211538
## 
## $standardError
##            cloudy       rain      sunny
## cloudy 0.05485322 0.05753055 0.05753055
## rain   0.07445856 0.06943297 0.05452441
## sunny  0.04111711 0.02087731 0.04376245
## 
## $confidenceInterval
## $confidenceInterval$confidenceLevel
## [1] 0.92
## 
## $confidenceInterval$lowerEndpointMatrix
##           cloudy       rain      sunny
## cloudy 0.2024315 0.25531059 0.25531059
## rain   0.3917409 0.20261227 0.07519196
## sunny  0.1624958 0.01813559 0.64169687
## 
## $confidenceInterval$upperEndpointMatrix
##           cloudy      rain     sunny
## cloudy 0.3941775 0.4565555 0.4565555
## rain   0.6522800 0.4452100 0.2643798
## sunny  0.3061880 0.0898865 0.7947391
## 
## 
## $logLikelihood
## [1] -169.9034

predictiveDistribution(weathersOfDays[1:200], 
                       weathersOfDays[201:365],hyperparam = hyperMatrix2)

## [1] -156.7419

inferHyperparam(transMatr = weatherMatrix, scale = c(10, 10, 10))

## $scaledInference
##        cloudy rain sunny
## cloudy      4    3     3
## rain        4    4     2
## sunny       2    1     7

inferHyperparam(data = weathersOfDays[1:15])

## $dataInference
##        cloudy rain sunny
## cloudy      4    4     2
## rain        4    3     1
## sunny       2    1     2

hyperMatrix3 <- inferHyperparam(transMatr = weatherMatrix, 
                                scale = c(10, 10, 10))
hyperMatrix3 <- hyperMatrix3$scaledInference
hyperMatrix4 <- inferHyperparam(data = weathersOfDays[1:15])
hyperMatrix4 <- hyperMatrix4$dataInference

data(preproglucacon)
preproglucacon <- preproglucacon[[2]][1:200]
MLEest <- markovchainFit(preproglucacon, method = "mle")
MAPest <- markovchainFit(preproglucacon, method = "map")
MLEest$estimate

## MLE Fit 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  A, C, G, T 
##  The transition matrix  (by rows)  is defined as follows: 
##           A          C         G         T
## A 0.4084507 0.07042254 0.2112676 0.3098592
## C 0.3809524 0.14285714 0.0952381 0.3809524
## G 0.3414634 0.07317073 0.1707317 0.4146341
## T 0.3181818 0.15151515 0.2424242 0.2878788

MAPest$estimate

## Bayesian Fit 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  A, C, G, T 
##  The transition matrix  (by rows)  is defined as follows: 
##           A          C         G         T
## A 0.4084507 0.07042254 0.2112676 0.3098592
## C 0.3809524 0.14285714 0.0952381 0.3809524
## G 0.3414634 0.07317073 0.1707317 0.4146341
## T 0.3181818 0.15151515 0.2424242 0.2878788

Definition and generator matrix

The package provides functionality for continuous time Markov chains (CTMCs). CTMCs are a generalization of discrete time Markov chains (DTMCs) in that we allow time to be continuous. We assume a finite state space \(S\), so that the generator matrix introduced below can be stored in memory. We can think of CTMCs as Markov chains in which state transitions can happen at any time.

More formally, we require our CTMCs to satisfy the following two properties:

• The Markov property - let \(F_{X(s)}\) denote the information about \(X\) up to time \(s\). Let \(j \in S\) and \(s \leq t\). Then, \(P(X(t) = j|F_{X(s)}) = P(X(t) = j|X(s))\).
• Time homogeneity - \(P(X(t) = j|X(s) = k) = P(X(t-s) = j|X(0) = k)\).

If both the above properties are satisfied, it is referred to as a time-homogeneous CTMC. If a transition occurs at time \(t\), then \(X(t)\) denotes the new state and \(X(t)\neq X(t-)\).

Now, let \(X(0)=x\) and let \(T_x\) be the time a transition occurs from this state. We are interested in the distribution of \(T_x\). For \(s,t \geq 0\), it can be shown that $ P(T_x > s+t | T_x > s) = P(T_x > t) $

This is the memoryless property that only the exponential random variable exhibits. Therefore, the exponential distribution is the relevant distribution, and each state \(s \in S\) has an exponential holding parameter \(\lambda(s)\). Since \(\mathrm{E}T_x = \frac{1}{\lambda(x)}\), the higher the rate \(\lambda(x)\), the smaller the expected time of transitioning out of the state \(x\).

However, specifying this parameter alone for each state gives an incomplete picture of our CTMC. To see why, consider a state \(x\) that may transition to either state \(y\) or \(z\). The holding parameter enables us to predict when a transition may occur if we start off in state \(x\), but tells us nothing about which state will be next.

To this end, we also need transition probabilities associated with the process, defined as follows (for \(y \neq x\)) - \(p_{xy} = P(X(T_x) = y | X(0) = x)\). Note that \(\sum_{y \neq x} p_{xy} = 1\). Let \(Q\) denote this transition matrix (\(Q_{ij} = p_{ij}\)). What is key here is that \(T_x\) and the state \(y\) are independent random variables. We define \(\lambda(x, y) = \lambda(x) p_{xy}\).

We now look at Kolmogorov’s backward equation. Let’s define \(P_{ij}(t) = P(X(t) = j | X(0) = i)\) for \(i, j \in S\). The backward equation is given by (it can be proved) \(P_{ij}(t) = \delta_{ij}e^{-\lambda(i)t} + \int_{0}^{t}\lambda(i)e^{-\lambda(i)s} \sum_{k \neq i} Q_{ik} P_{kj}(t-s) ds\). The first term is non-zero if and only if \(i=j\) and represents the probability that the first transition from state \(i\) occurs after time \(t\), so that at time \(t\) the state is still \(i\). The second term accounts for the transitions that occur before time \(t\): it is the probability of being in state \(j\) at time \(t\) after a first transition at time \(s \le t\).

In compact form, the equation reads \(P'(t) = AP(t)\), where \(A\) is the generator matrix: \[ A(i, j) = \begin{cases} \lambda(i, j) & \mbox{if } i \neq j \\ -\lambda(i) & \mbox{else.} \end{cases} \] The entries of each row sum to zero, and a CTMC is completely specified by its generator matrix.

Stationary distributions

The following theorem concerns the existence of a stationary distribution for CTMCs. Note that \(X(t)\) is irreducible and recurrent if and only if the embedded chain \(X_n\) is irreducible and recurrent.

Suppose that \(X(t)\) is irreducible and recurrent. Then \(X(t)\) has an invariant measure \(\eta\), which is unique up to multiplicative factors. Moreover, for each \(k \in S\), we have

\[\eta_k = \frac{\pi_k}{\lambda(k)}\]

where \(\pi\) is the invariant measure of the embedded discrete time Markov chain \(X_n\), unique up to multiplicative factors. Finally, \(\eta\) satisfies

\[0 < \eta_j < \infty, \forall j \in S\]

and if \(\sum_i \eta_i < \infty\) then \(\eta\) can be normalized to get a stationary distribution.

Estimation of the generator matrix

Let the data set be \(D = \{(s_0, t_0), (s_1, t_1), ..., (s_{N-1}, t_{N-1})\}\) where \(N=|D|\). Each \(s_i\) is a state from the state space \(S\) and during the time \([t_i,t_{i+1}]\) the chain is in state \(s_i\). Let the parameters be represented by \(\theta = \{\lambda, P\}\) where \(\lambda\) is the vector of holding parameters for each state and \(P\) the transition matrix of the embedded discrete time Markov chain.

Then the probability is given by

\[ {Pr(D | \theta) \propto \lambda(s_0)e^{-\lambda(s_0)(t_1-t_0)}Pr(s_1|s_0) \cdot\ldots\cdot \lambda(s_{N-2})e^{-\lambda(s_{N-2})(t_{N-1}-t_{N-2})}Pr(s_{N-1}|s_{N-2})} \]

Let \(n(j|i)\) denote the number of \(i \to j\) transitions in \(D\), and \(n(i)\) the number of times \(s_i\) occurs in \(D\). Let \(t(s_i)\) denote the total time the chain spends in state \(s_i\).

Then the MLEs are given by

\[ \hat{\lambda}(s) = \frac{n(s)}{t(s)}, \qquad \widehat{Pr}(j|i)=\frac{n(j|i)}{n(i)}. \]

Expected hitting time

The function ExpectedTime calculates the expected hitting time from one state to another. Let \(j\) be the target state and let \(q_{kl}\) be the entries of the generator matrix, with holding rate \(q_k=-q_{kk} > 0\) for every \(k \neq j\). Under this condition, the vector \(p\) of the expected hitting times of \(j\) is the minimal non-negative solution of the system of linear equations \[\begin{equation} \begin{cases} p_{k} = 0 & k = j \\ -\sum_{l \in S} q_{kl}p_{l} = 1 & k \neq j. \end{cases} \label{eq:EHT} \end{equation}\]

Transition probabilities over time

The function probabilityatT computes the probability of every state at time \(t\) for a ctmc object. It solves Kolmogorov’s backward equation, whose solution is \(P(t) = e^{tA}\) for \(t \geq 0\) (with \(P(0) = I\)), where \(A\) is the generator matrix. The entry \(P_{ij}(t)\) is the probability that the chain is in state \(j\) at time \(t\) given that it was in state \(i\) at time \(0\). The function also handles generators stored by columns. If the initial state is not provided, it returns the whole transition matrix \(P(t)\).

Examples

To create a CTMC object, one provides a valid generator matrix. A ctmc object has the slots states, generator, byrow and name (see the documentation of the class for details). Consider the following example in which we aim to model the transition of a molecule from the \(\sigma\) state to the \(\sigma^*\) state. When in the former state, if it absorbs sufficient energy, it can make the jump to the latter state and remains there for some time before transitioning back to the original state. Let us model this by a CTMC:

energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3,
                       1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
molecularCTMC <- new("ctmc", states = energyStates, 
                 byrow = byRow, generator = gen, 
                 name = "Molecular Transition Model")      

To generate random CTMC transitions, we provide an initial distribution of the states. This must be in the same order as the dimnames of the generator. The output can be returned either as a list or a data frame.

statesDist <- c(0.8, 0.2)
rctmc(n = 3, ctmc = molecularCTMC, initDist = statesDist, out.type = "df", include.T0 = FALSE)

##       states               time
## 1 sigma_star 0.0204442738157116
## 2      sigma   1.27491917071131
## 3 sigma_star   1.39663587353391

The argument n is the number of transitions to generate. The optional argument T is the time at which the simulation ends; to use it, set n to Inf (the number of transitions is not known beforehand) and set T accordingly.

statesDist <- c(0.8, 0.2)
rctmc(n = Inf, ctmc = molecularCTMC, initDist = statesDist, T = 2)

## [[1]]
## [1] "sigma"      "sigma_star" "sigma"      "sigma_star" "sigma"     
## [6] "sigma_star"
## 
## [[2]]
## [1] 0.00000000 0.07218092 0.07432856 0.78912847 1.58864151 1.78023148

The stationary distribution is returned by the steadyStates function:

steadyStates(molecularCTMC)

##      sigma sigma_star
## [1,]  0.25       0.75

The expected time needed to move from the state sigma to the state sigma_star is \(1/3\) of a time unit, the inverse of the holding rate of sigma, since the only possible transition leads to the target. The probabilities of the two states half a time unit after the chain was in sigma follow from \(P(t) = e^{tA}\):

ExpectedTime(molecularCTMC, 1, 2)

## [1] 0.3333333

probabilityatT(molecularCTMC, 0.5, 1)

## [1] 0.3515015 0.6484985

The ctmcFit function fits a CTMC from a sequence of states and the times at which they are observed. It returns the maximum likelihood estimates of the parameters along with confidence intervals.

data <- list(c("a", "b", "c", "a", "b", "a", "c", "b", "c"), 
             c(0, 0.8, 2.1, 2.4, 4, 5, 5.9, 8.2, 9))
ctmcFit(data)

## $estimate
## An object of class "ctmc"
## Slot "states":
## [1] "a" "b" "c"
## 
## Slot "byrow":
## [1] TRUE
## 
## Slot "generator":
##            a          b          c
## a -0.9090909  0.6060606  0.3030303
## b  0.3225806 -0.9677419  0.6451613
## c  0.3846154  0.3846154 -0.7692308
## 
## Slot "name":
## [1] ""
## 
## 
## $errors
## $errors$dtmcConfidenceInterval
## $errors$dtmcConfidenceInterval$confidenceLevel
## [1] 0.95
## 
## $errors$dtmcConfidenceInterval$lowerEndpointMatrix
##            a          b          c
## a 0.00000000 0.20765960 0.06149194
## b 0.06149194 0.00000000 0.20765960
## c 0.09453121 0.09453121 0.00000000
## 
## $errors$dtmcConfidenceInterval$upperEndpointMatrix
##           a         b         c
## a 0.5614970 0.9385081 0.7923404
## b 0.7923404 0.5614970 0.9385081
## c 0.9054688 0.9054688 0.6576198
## 
## 
## $errors$lambdaConfidenceInterval
## $errors$lambdaConfidenceInterval$lowerEndpointVector
## [1] 0 0 0
## 
## $errors$lambdaConfidenceInterval$upperEndpointVector
## [1] 1.937805 2.062825 1.835311

One approach to obtain the generator matrix is to apply the logm function from the package on a transition matrix. Numeric issues arise, see . For example, applying the standard method (‘Higham08’) on mcWeather raises an error, whilst the alternative method (eigenvalue decomposition) is OK. The following code estimates the generator matrix of the mcWeather transition matrix.

mcWeatherQ <- expm::logm(mcWeather@transitionMatrix,method='Eigen')
mcWeatherQ

##             sunny     cloudy        rain
## sunny  -0.4633266  0.4270614  0.03626515
## cloudy  0.6405921 -1.8129810  1.17238885
## rain    0.0725303  1.5631851 -1.63571543

The “half-day” transition probability of the mcWeather chain is therefore

mcWeatherHalfDayTM <- expm::expm(mcWeatherQ*.5)
mcWeatherHalfDay <- new("markovchain",transitionMatrix=mcWeatherHalfDayTM,name="Half Day Weather Transition Matrix")
mcWeatherHalfDay

## Half Day Weather Transition Matrix 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  sunny, cloudy, rain 
##  The transition matrix  (by rows)  is defined as follows: 
##             sunny    cloudy       rain
## sunny  0.81711847 0.1367957 0.04608582
## cloudy 0.20519357 0.5226767 0.27212969
## rain   0.09217163 0.3628396 0.54498878

The package provides various functions to estimate the generator matrix (GM) of a CTMC process using different methods. The following code provides a way to join and computations.

if(requireNamespace(package='ctmcd', quietly = TRUE)) {
require(ctmcd)
require(expm)
#defines a function to transform a GM into a TM
gm_to_markovchain<-function(object, t=1) {
  if(!inherits(object, c("gm","matrix","Matrix")))
    stop("Error! Expecting either a matrix or a gm object")
  if (inherits(object, c("matrix","Matrix"))) generator_matrix<-object else generator_matrix<-as.matrix(object[["par"]])
  #must add importClassesFrom("markovchain",markovchain) in the NAMESPACE
  #must add importFrom(expm, "expm")
  transitionMatrix<-expm(generator_matrix*t)
  out<-as(transitionMatrix,"markovchain")
  return(out)
}
#loading ctmcd dataset
data(tm_abs)
gm0=matrix(1,8,8) #initializing
diag(gm0)=0
diag(gm0)=-rowSums(gm0)
gm0[8,]=0
gmem=gm(tm_abs,te=1,method="EM",gmguess=gm0) #estimating GM
mc_at_2=gm_to_markovchain(object=gmem, t=2) #converting to TM at time 2
} else {
  warning('package ctmcd unavailable')
}

## Loading required package: ctmcd

## Loading required package: expm

## 
## Attaching package: 'expm'

## The following object is masked from 'package:Matrix':
## 
##     expm

Weather forecasting

Markov chains provide a simple model to predict the next day’s weather given the current meteorological condition. The first application herewith shown is the “Land of Oz example” from , the second is the “Alofi Island Rainfall” from .

mcWP <- new("markovchain", states = c("rainy", "nice", "snowy"),
         transitionMatrix = matrix(c(0.5, 0.25, 0.25,
                                   0.5, 0, 0.5,
                                   0.25,0.25,0.5), byrow = T, nrow = 3))

W0 <- t(as.matrix(c(0, 1, 0)))
W1 <- W0 * mcWP; W1

##      rainy nice snowy
## [1,]   0.5    0   0.5

W2 <- W0 * (mcWP ^ 2); W2

##      rainy nice snowy
## [1,] 0.375 0.25 0.375

W3 <- W0 * (mcWP ^ 3); W3

##        rainy   nice   snowy
## [1,] 0.40625 0.1875 0.40625

W7 <- W0 * (mcWP ^ 7)
W7

##          rainy      nice     snowy
## [1,] 0.4000244 0.1999512 0.4000244

q <- steadyStates(mcWP)
q

##      rainy nice snowy
## [1,]   0.4  0.2   0.4

R0 <- t(as.matrix(c(1, 0, 0)))
R7 <- R0 * (mcWP ^ 7); R7

##          rainy      nice     snowy
## [1,] 0.4000244 0.2000122 0.3999634

S0 <- t(as.matrix(c(0, 0, 1)))
S7 <- S0 * (mcWP ^ 7); S7

##          rainy      nice     snowy
## [1,] 0.3999634 0.2000122 0.4000244

data("rain", package = "markovchain")
table(rain$rain)

## 
##   0 1-5  6+ 
## 548 295 253

mcAlofi <- markovchainFit(data = rain$rain, name = "Alofi MC")$estimate
mcAlofi

## Alofi MC 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  0, 1-5, 6+ 
##  The transition matrix  (by rows)  is defined as follows: 
##             0       1-5        6+
## 0   0.6605839 0.2299270 0.1094891
## 1-5 0.4625850 0.3061224 0.2312925
## 6+  0.1976285 0.3122530 0.4901186

steadyStates(mcAlofi)

##              0       1-5        6+
## [1,] 0.5008871 0.2693656 0.2297473

Finance and Economics

Other relevant applications of DTMC can be found in Finance and Economics.

rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
creditMatrix <- matrix(
  c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01,
    0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01,
    0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06,
    0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18,
    0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06,
    0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20,
    0.21, 0, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79,
    0, 0, 0, 0, 0, 0, 0, 100
   )/100, 8, 8, dimnames = list(rc, rc), byrow = TRUE)

creditMc <- new("markovchain", transitionMatrix = creditMatrix, 
                name = "S&P Matrix")
absorbingStates(creditMc)

## [1] "D"

statesNames <- c("customer", "non customer")
P <- markovchain:::zeros(2); P[1, 1] <- .9; P[1, 2] <- .1; P[2, 2] <- .95; P[2, 1] <- .05;
rownames(P) <- statesNames; colnames(P) <- statesNames
mcP <- new("markovchain", transitionMatrix = P, name = "Telephone company")
M <- markovchain:::zeros(2); M[1, 1] <- -20; M[1, 2] <- -30; M[2, 1] <- -40; M[2, 2] <- 0

c1 <- 100 + conditionalDistribution(mcP, state = "customer") %*% M[1,]
c2 <- 0 + conditionalDistribution(mcP, state = "non customer") %*% M[2,]

as.numeric((c(1, 0)* mcP ^ 5) %*% (as.vector(c(c1, c2))))

## [1] 48.96009

rho <- 0.979
sigmaEps <- 0.0072
theoreticalVar <- sigmaEps^2 / (1 - rho^2)
theoreticalVar

## [1] 0.001247383

relativeError <- function(size) {
  outT <- tauchen(alpha = 0, sigma = sigmaEps, rho = rho, size = size)
  outR <- rouwenhorst(alpha = 0, sigma = sigmaEps, rho = rho, size = size)
  piT <- as.numeric(steadyStates(outT$chain))
  piR <- as.numeric(steadyStates(outR$chain))
  varT <- sum(piT * (outT$states - sum(piT * outT$states)) ^ 2)
  varR <- sum(piR * (outR$states - sum(piR * outR$states)) ^ 2)
  c(tauchen = 100 * abs(varT - theoreticalVar) / theoreticalVar,
    rouwenhorst = 100 * abs(varR - theoreticalVar) / theoreticalVar)
}

sapply(c(5, 10, 25), relativeError)

##                     [,1]         [,2]         [,3]
## tauchen     8.139476e+01 4.979587e+01 9.328706e+00
## rouwenhorst 6.953450e-14 1.112552e-12 1.077785e-12

Actuarial science

Markov chains are widely applied in the field of actuarial science. Two classical applications are policyholders’ distribution across Bonus Malus classes in Motor Third Party Liability (MTPL) insurance (Section @ref(sec:bm)) and health insurance pricing and reserving (Section @ref(sec:hi)).

getBonusMalusMarkovChain <- function(lambda) {
    bmMatr <- markovchain:::zeros(5)
    bmMatr[1, 1] <- dpois(x = 0, lambda)
    bmMatr[1, 3] <- dpois(x = 1, lambda)
    bmMatr[1, 5] <- 1 - ppois(q = 1, lambda)
    
    bmMatr[2, 1] <- dpois(x = 0, lambda)
    bmMatr[2, 4] <- dpois(x = 1, lambda)
    bmMatr[2, 5] <- 1 - ppois(q = 1, lambda)
    
    bmMatr[3, 2] <- dpois(x = 0, lambda)
    bmMatr[3, 5] <- 1 - dpois(x=0, lambda)
 
    bmMatr[4, 3] <- dpois(x = 0, lambda)
    bmMatr[4, 5] <- 1 - dpois(x = 0, lambda)
  
    bmMatr[5, 4] <- dpois(x = 0, lambda)
    bmMatr[5, 5] <- 1 - dpois(x = 0, lambda)
    stateNames <- as.character(1:5)
    out <- new("markovchain", transitionMatrix = bmMatr, 
             states = stateNames, name = "BM Matrix")
    return(out)
}

bmMc <- getBonusMalusMarkovChain(0.05)
as.numeric(steadyStates(bmMc))

## [1] 0.895836079 0.045930498 0.048285405 0.005969247 0.003978772

sum(as.numeric(steadyStates(bmMc)) * c(0.5, 0.7, 0.9, 1, 1.25))

## [1] 0.534469

ltcDemoPath<-system.file("extdata", "ltdItaData.txt", 
                         package = "markovchain")
ltcDemo<-read.table(file = ltcDemoPath, header=TRUE, 
                    sep = ";", dec = ".")
head(ltcDemo)

##   age          pAD        pID          pAI       pAA
## 1  20 0.0004616002 0.01083364 0.0001762467 0.9993622
## 2  21 0.0004824888 0.01079719 0.0001710577 0.9993465
## 3  22 0.0004949938 0.01177076 0.0001592333 0.9993458
## 4  23 0.0005042935 0.01159394 0.0001605731 0.9993351
## 5  24 0.0005074193 0.01260574 0.0001606504 0.9993319
## 6  25 0.0005154267 0.01526364 0.0001643603 0.9993202

ltcDemo<-transform(ltcDemo,
                   pIA=0,
                   pII=1-pID,
                   pDD=1,
                   pDA=0,
                   pDI=0)

possibleStates<-c("A","I","D")
getMc4Age<-function(age) {
  transitionsAtAge<-ltcDemo[ltcDemo$age==age,]
  
  myTransMatr<-matrix(0, nrow=3,ncol = 3,
                      dimnames = list(possibleStates, possibleStates))
  myTransMatr[1,1]<-transitionsAtAge$pAA[1]
  myTransMatr[1,2]<-transitionsAtAge$pAI[1]
  myTransMatr[1,3]<-transitionsAtAge$pAD[1]
  myTransMatr[2,2]<-transitionsAtAge$pII[1]
  myTransMatr[2,3]<-transitionsAtAge$pID[1]
  myTransMatr[3,3]<-1
  
  myMc<-new("markovchain", transitionMatrix = myTransMatr,
            states = possibleStates,
            name = paste("Age",age,"transition matrix"))
  
  return(myMc)

}

getFullTransitionTable<-function(age){
  ageSequence<-seq(from=age, to=120)
  k=1
  myList=list()
  for ( i in ageSequence) {
    mc_age_i<-getMc4Age(age = i)
    myList[[k]]<-mc_age_i
    k=k+1
  }
  myMarkovChainList<-new("markovchainList", markovchains = myList,
                         name = paste("TransitionsSinceAge", age, sep = ""))
  return(myMarkovChainList)
}
transitionsSince100<-getFullTransitionTable(age=100)

rmarkovchain(n = 10, object = transitionsSince100,
             what = "matrix", t0 = "A", include.t0 = TRUE)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,] "A"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [2,] "A"  "A"  "A"  "D"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [3,] "A"  "A"  "A"  "A"  "A"  "A"  "A"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [4,] "A"  "A"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [5,] "A"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [6,] "A"  "A"  "A"  "I"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [7,] "A"  "A"  "A"  "I"  "I"  "I"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [8,] "A"  "A"  "A"  "A"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##  [9,] "A"  "A"  "A"  "A"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
## [10,] "A"  "A"  "I"  "I"  "D"  "D"  "D"  "D"  "D"  "D"   "D"   "D"   "D"  
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
##  [1,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [2,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [3,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [4,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [5,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [6,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [7,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [8,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
##  [9,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"  
## [10,] "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"   "D"

transitionsSince80<-getFullTransitionTable(age=80)
lifeTrajectories<-rmarkovchain(n=1e2, object=transitionsSince80,
                               what="matrix",t0="A",include.t0=TRUE)
temp<-matrix(0,nrow=nrow(lifeTrajectories),ncol = ncol(lifeTrajectories))
temp[lifeTrajectories=="I"]<-1
expected_period_disabled<-mean(rowSums((temp)))
expected_period_disabled