En este tutorial se utilizará la librería markovchain de R para trabajar con cadenas de Markov continuas y discretas. En primer lugar, instalaremos el paquete markovchain y lo llamaremos:

library(markovchain)
## Package:  markovchain
## Version:  0.6.9.8-1
## Date:     2017-08-15
## BugReport: http://github.com/spedygiorgio/markovchain/issues

Para crear una cadena de Markov discreta o continua, utilizaremos la función new. Esta función recibe, entre otros, los parámetros Class,states, byrow,y transitionMatrix. El parámetro Class es de tipo String y permite especificar el tipo de cadena de Markov a crear. Para crear una CMTD usamos Class=“markovchain”, y para crear una CMTC usamos Class=“ctmc”. El parámetro states es un vector de tipo Character que indica el espacio de estados de la cadena de Markov. El parámetro byrow de tipo Boolean indica si la matriz de transición se encuentra organizada por filas o por columnas. Así pues, si la matriz se especifica por filas, debemos asignar by.row=TRUE. Por otro lado, el parámetro transitionMatrix es una matriz con las probabilidades o tasas de transición entre estados. Enseguida se presentan ejemplos sobre cómo crear cadenas de Markov de tiempo discreto y continuo.

Ejemplo 1. Cadena de Markov de Tiempo Discreto Ergódica.

Sea \(X_n\) el clima en Bogotá el día \(n\), donde \(S_x=\{Soleado, Nublado, Lluvioso\}\). La matriz de probabilidades de transición a un paso se presenta enseguida:

\[ \textbf{P}= \begin{bmatrix} 0.7 & 0.2 & 0.1\\ 0.3 & 0.4 & 0.3\\ 0.2 & 0.45 & 0.35\\ \end{bmatrix}\]

En primer lugar, vamos a crear un vector con los estados, y una matriz con las probabilidades de transición:

estadosClima = c("Soleado", "Nublado", "Lluvioso")
probabilidadesClima = matrix(data = c(0.70, 0.2, 0.1, 0.3, 0.4, 0.3, 0.2, 0.45, 0.35), byrow = TRUE, nrow = 3, ncol=3)

probabilidadesClima
##      [,1] [,2] [,3]
## [1,]  0.7 0.20 0.10
## [2,]  0.3 0.40 0.30
## [3,]  0.2 0.45 0.35

Ahora bien, crearemos esta cadena de Markov en R, con el nombre cmClima. Especificamos como parámetros el espacio de estados, y la matriz de probabilidades de transición definidos antes.

cmClima = new(Class="markovchain", states = estadosClima, byrow = TRUE, transitionMatrix = probabilidadesClima)

Una vez creada la cadena de Markov, podemos conocer sus estados con el comando states, su número de estados con el comando dim, y su matriz de probabilidades de transición con el comando print:

states(cmClima)
## [1] "Soleado"  "Nublado"  "Lluvioso"
dim(cmClima)
## [1] 3
print(cmClima)
##          Soleado Nublado Lluvioso
## Soleado      0.7    0.20     0.10
## Nublado      0.3    0.40     0.30
## Lluvioso     0.2    0.45     0.35

Para acceder a la probabilidad de transición entre dos estados particulares, utilizamos el comando transitionProbability.

transitionProbability(cmClima, "Nublado", "Lluvioso")
## [1] 0.3

Supongamos que queremos hallar la probabilidad de que la cadena se encuentre en cada uno de los estados en dos días, teniendo un vector de probabilidades iniciales: \(\alpha=(0, 1, 0)\). Para esto, en primer lugar, creamos el vector de probabilidades iniciales en R:

alpha <- c(0, 1, 0)

Ahora bien, para hallar la probabilidad solicitada debemos hacer la siguiente operación:

\[ \alpha \cdot \textbf{P}^2 \] Así pues, realizamos lo siguiente:

probabilidadesDosPasos <- alpha * (cmClima^2)

probabilidadesDosPasos
##      Soleado Nublado Lluvioso
## [1,]    0.39   0.355    0.255

Podemos verificar si la cadena de Markov es irreducible utilizando el método is.irreducible. Para cadenas de Markov de tiempo discreto irreducibles, también podemos hallar el periodo utilizando el método period:

is.irreducible(cmClima)
## [1] TRUE
period(cmClima)
## [1] 1

Debido a que la cadena es aperiódica e irreducible, concluimos que es ergódica y podemos hallar las probabilidades en estado estable, utilizando el comando steadyStates:

steadyStates(cmClima)
##        Soleado   Nublado  Lluvioso
## [1,] 0.4636364 0.3181818 0.2181818

Adicionalmente, podemos verificar si un estado \(j\) es alcanzable desde un estado \(i\), utilizando el método is.accessible. Este método recibe los parámetros object, from y to. El parámetro object es la cadena de Markov, el parámetro from es el estado inicial \(i\), y el parámetro to es el estado final \(j\).

is.accessible(object=cmClima,from="Soleado", to="Lluvioso")
## [1] TRUE

Ejemplo 2. Cadena de Markov de Tiempo Continuo.

Sea \(\{X(t), t \geq 0\}\) el estado del clima en el tiempo \(t\), y sea \(S_x=\{Soleado, Nublado, Lluvioso\}\). La matriz de tasas de transición (\(\text{días}^{-1}\)) se presenta enseguida:

\[ \textbf{Q}= \begin{bmatrix} -4.00 & 2.43 & 1.57\\ 4.29 & -20.12 & 15.83\\ 4.41 & 24.18 & -28.59\\ \end{bmatrix}\]

En primer lugar, vamos a crear un vector con los estados, y una matriz con las probabilidades de transición:

estadosContinua = c("Soleado", "Nublado", "Lluvioso")

tasas = matrix(data = c(-4,2.43,1.57,4.29,-20.12,15.83,4.41,24.18,-28.59), byrow = TRUE, nrow = 3, ncol=3)

Ahora, creamos la cadena de Markov utilizando el parámetro class=“ctmc”:

cadenaContinua <- new(Class="ctmc", states = estadosContinua,
byrow = TRUE, generator = tasas)

Para verificar si una cadena de Markov de tiempo continuo es irreducible, utilizamos el método is.CMTCirreducible:

is.CTMCirreducible(cadenaContinua)
## [1] TRUE

Ya que la cadena es irreducible y el espacio de estados es finito, podemos hallar las probabilidades en estado estable:

steadyStates(cadenaContinua)
##       Soleado   Nublado  Lluvioso
## [1,] 0.520232 0.2904053 0.1893627

Ejemplo 3: Ejercicio CMTD

Realizaremos el Problema 1 del archivo Complementaria 4.pdf que se encuentra en Sicua Plus. Se modelará la situación como una cadena de Markov de tiempo discreto. Para esta situación se definen dos variables de estado:

\[\begin{eqnarray*} & Y_n = \text{ resultado de la revisión }n \\ & Z_n = \text{ resultado de la revisión }n-1\\ & X_n = (Z_n, Y_n) & \\ & S_X = \{(0,0),(0,1),(1,1),(1,2),(0,2),(2,0),(2,2)\} & \end{eqnarray*}\]

Donde \(0\) indica perfecto estado, \(1\) indica un daño de nivel amarillo y \(2\) un daño de nivel rojo.

Tenemos la siguiente matriz de probabilidades:

\[\textbf{P}= \begin{bmatrix} 0.6 & 0.25 & 0 & 0 & 0.15 & 0 & 0 \cr 0 & 0 & 0.7 & 0.3 & 0 & 0 & 0 \cr 0 & 0 & 0.4 & 0.6 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5 \cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5 \cr 1 & 0 & 0 & 0 & 0 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5\cr \end{bmatrix}\]

Vamos a crear el espacio de estados y la matriz de tasas de transición en R:

estadosEjercicio = c("(0,0)","(0,1)","(1,1)","(1,2)","(0,2)","(2,0)","(2,2)")

probabilidadesEjercicio= matrix(data = c(0.6, 0.25,0,0,0.15,0,0,0,0,0.7,0.3,0,0,0,0,0,0.4,0.6,0,0,0,0,0,0,0,0,0.5,0.5,0,0,0,0,0,0.5,0.5,1,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5), byrow = TRUE, nrow = 7, ncol=7)

probabilidadesEjercicio
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  0.6 0.25  0.0  0.0 0.15  0.0  0.0
## [2,]  0.0 0.00  0.7  0.3 0.00  0.0  0.0
## [3,]  0.0 0.00  0.4  0.6 0.00  0.0  0.0
## [4,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5
## [5,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5
## [6,]  1.0 0.00  0.0  0.0 0.00  0.0  0.0
## [7,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5

Ahora, creamos la cadena:

cmEjercicio = new(Class="markovchain", states = estadosEjercicio, byrow = TRUE,
transitionMatrix = probabilidadesEjercicio)

Vamos a verificar que la cadena sea ergódica. Para esto, hallamos el periodo y verificamos si es irreducible:

is.irreducible(cmEjercicio) 
## [1] TRUE
period(cmEjercicio)
## [1] 1

Ya que la cadena es ergódica, podemos hallar las probabilidades en estado estable:

steadyStates(cmEjercicio) 
##          (0,0)      (0,1)    (1,1)      (1,2)      (0,2)     (2,0)
## [1,] 0.3647416 0.09118541 0.106383 0.09118541 0.05471125 0.1458967
##          (2,2)
## [1,] 0.1458967