Complementaria 7: Cadenas Absorbentes en Python

En este tutorial se abordará el tema de análisis de cadenas de Markov absorbentes a través de dos ejemplos, uno continuo y discreto. Para desarrollar estos ejercicio utilizaremos las librerías jmarkov y numpy

import numpy as np
from jmarkov.dtmc import dtmc
from jmarkov.ctmc import ctmc

1. Cadena de Markov de Tiempo Continuo

Literal a

Vamos a realizar el primer problema del archivo “Complementaria 7(Q).pdf” que se encuentra en Bloque Neon. En primer lugar, se modelará la situación como una cadena de Markov de tiempo continuo:

\[ X(t) = \text{Estado del paciente en la EPS en el tiempo } t \]

\[ S_X = \{D, TNI, TI, CEPS, M\} \]

\[\begin{split} \mathbb{Q}_{i \rightarrow j} = \begin{bmatrix} -(\phi + \omega + \mu + \alpha ) & \phi & \omega & \mu & \alpha \\ 0 & -(\gamma + 0.5\mu + \kappa) & \gamma & 0.5\mu & \kappa \\ 0 & 0 & -(0.2\mu + \tau) & 0.2\mu & \tau \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix} \end{split}\]

Es importante notar que los estados absorbentes Cambio de EPS ($CEPS$) y Muerte ($M$) se representan en la matriz de tasas $\mathbb{Q}$ por tener filas llenas de 0 únicamente, ya que no existe ninguna tasa de salida de un estado absorbente. Con la formulación, podemos definir la lista de estados y la matriz generadora.

estados = ['D',     # Diagnosticado
           'TNI',   # Tratamieno No Invasivo
           'TI',    # Tratamiento Invasivo
           'CEPS',  # Cambio de EPS
           'M']     # Muerte

phi = 3
omega = 4
gamma = 3.5
mu = 1
alpha = 2
kappa = 2.5
tao = 1.5

matriz = np.array([[-(phi + omega + mu + alpha), phi, omega, mu, alpha], 
                    [0, -(gamma + 0.5*mu + kappa), gamma, 0.5*mu, kappa],
                    [0, 0, -(0.2*mu + tao), 0.2*mu, tao],
                    [0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0]])

Q = ctmc(matriz, states = np.array(estados))

Literal b

El literal b solicita determinar cuánto tiempo permanece un paciente diagnosticado en el mismo estado antes de morir ($M$) o cambiar de EPS ($CEPS$), siendo estos dos últimos los estados absorbentes.

Este análisis puede resolverse utilizando los tiempos antes de absorción. Recordando la estructura de la matriz de tasas de transición para una cadena de Markov absorbente en tiempo continuo, esta tiene la forma:

\[\begin{split} \mathbb{Q} = \begin{bmatrix} \mathbb{U} & \mathbb{V} \\ 0 & 0 \end{bmatrix} \end{split}\]

Donde $\mathbb{U}$ contiene las tasas de transición entre los estados transitorios. Así, la matriz de tiempos antes de absorción en una cadena continua está dada por:

\[ (-\mathbb{U})^{-1} \]

En esta matriz, el elemento $[i, j]$ representa el tiempo promedio que la cadena pasa en el estado transitorio $j$, dado que comenzó en el estado transitorio $i$, antes de alcanzar un estado absorbente.

Por lo tanto, para responder el literal, basta con calcular esta matriz y seleccionar el componente $[i, i]$, donde $i$ corresponde al estado Diagnosticado ($D$), ya que se desea conocer el tiempo que el paciente permanece en ese mismo estado antes de ser absorbido. Con la libreria jmarkov este tiempo se puede calcular usando el atributo absorbtion_times

Q.absorbtion_times(start = 'D', target = 'D')

array([[0.1]])

Literal c

Para resolver este literal, podemos hacer uso de las probabilidades de absorción, las cuales están dadas por la siguiente matriz:

\[ -\mathbb{U}^{-1} \cdot \mathbb{V} \]

Donde $\mathbb{U}$ contiene las tasas de transición entre los estados transitorios, y $\mathbb{V}$ aquellas entre estados transitorios y absorbentes. La componente $[i, j]$ de esta matriz indica la probabilidad de que la cadena sea absorbida por el estado $j$, dado que inició en el estado $i$.

Teniendo esto en cuenta, para responder el literal, se debe seleccionar la componente en la que $i$ corresponde al estado Tratamiento No Invasivo ($TNI$) y $j$ al estado de ($CEPS$). Esto se puede hacer con el atributo absorbtion_probabilities de la libreria jmarkov

Q.absorbtion_probabilities(start = 2, target = 3)

array([[0.11764706]])

2. Cadena de Markov de Tiempo Discreto

Literal a

Continuando con el segundo problema del archivo “Complementaria 7(Q).pdf” que se encuentra en Bloque Neon. En primer lugar, se modelará la situación como una cadena de Markov de tiempo discreto:

\[\begin{split} X_n = \text{Clasificación del proyecto al final del año } n\\ \end{split}\]

\[\begin{split} S_X = \{\text{Poco Rentable (PR)},\text{Rentable (R)},\text{Muy Rentable (MR)},\text{No Rentable (NR)},\text{Liquidado (L)}\}\\ \end{split}\]

\[\begin{split} \mathbb{P}_{i \rightarrow j} = \begin{bmatrix} 0.3 & 0.2 & 0 & 0.3 & 0.2 \\ 0.1 & 0.5 & 0.2 & 0.1 & 0.1 \\ 0 & 0.2 & 0.8 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{split}\]

Es importante notar que los estados absorbentes No Rentable ($NR$) y Liquidado ($L$) se representan en la matriz de transición $\mathbb{P}$ mediante la condición $p_{ij} = 0$ y $p_{ii} = 1$, ya que $i$ es absorbente, lo cual indica que, una vez que se alcanza un estado absorbente, no existe posibilidad de ir a otro estado.

Con la formulación, podemos definir la lista de estados y la matriz de probabilidades.

estados = ['PR','R','MR','NR','L']

matriz = np.array([[0.3, 0.2, 0.0, 0.3, 0.2],
                [0.1, 0.5, 0.2, 0.1, 0.1],
                [0.0, 0.2, 0.8, 0.0, 0.0],
                [0.0, 0.0, 0.0, 1.0, 0.0],
                [0.0, 0.0, 0.0, 0.0, 1.0]])

P = dtmc(matriz, states =np.array(estados))

Literal b

Para verificar si la afirmación del analista es correcta, se debe calcular el tiempo que permanece en el portafolio un proyecto que entra con la clasificación Muy Rentable ($MR$), y compararlo con el tiempo que permanece un proyecto que ingresa con la clasificación Poco Rentable ($PR$). Estos tiempos pueden estimarse mediante el tiempo esperado antes de absorción, el cual está dado por la matriz:

\[ (\mathbb{I} - \mathbb{Q})^{-1} \]

Cada elemento $[i,j]$ de esta matriz representa el número esperado de veces que la cadena estará en el estado transitorio $j$, dado que comenzó en el estado transitorio $i$, antes de alcanzar un estado absorbente. La librería jmarkov permite calcular estos tiempos mediante el atributo absorbtion_times, el cual recibe como parámetros el estado inicial y el estado final.

A modo de ejemplo, se calculará el tiempo que permanece un proyecto en la clasificación Rentable ($R$), dado que comenzó en la clasificación Poco Rentable ($PR$).

P.absorbtion_times(start = 0, target = 1) # Tiempo antes de absorción desde el estado 'Poco Rentable' (índice 0) hasta el estado 'Rentable' (índice 1)

1.0526315789473688

El resultado obtenido fue 1.0526, lo que indica que, en promedio, un proyecto que comienza en la clasificación Poco Rentable tarda aproximadamente 1.05 años en alcanzar la clasificación Rentable, antes de ser absorbido por un estado como No Rentable o Liquidado.

Retomando el ejercicio, para calcular el tiempo total que permanece un proyecto en el portafolio, se deben sumar los tiempos antes de absorción en cada uno de los tres estados transitorios. La suma de estos tiempos representa el número esperado de años que el proyecto permanece en el portafolio antes de ser absorbido por uno de los dos estados absorbentes (No Rentable o Liquidado).

Estos tiempos se presentan a continuación y se denotarán como TAA_MR y TAA_PR, que significan Tiempos Antes de Absorción para los proyectos que ingresan con clasificación Muy Rentable ($MR$) y Poco Rentable ($PR$), respectivamente.

TAA_MR = sum(P.absorbtion_times(start=estados.index('MR'), target = i) for i in range(0, 3))
TAA_PR = sum(P.absorbtion_times(start=estados.index('PR'), target = i) for i in range(0, 3))

Finalmente, la división entre los dos tiempos (TAA_MR / TAA_PR) permite evaluar si la afirmación del analista es correcta.

Si el resultado es mayor a 4, se concluye que el analista tenía razón al afirmar que un proyecto clasificado inicialmente como Muy Rentable ($MR$) permanece al menos cuatro veces más tiempo en el portafolio que uno clasificado como Poco Rentable ($PR$).

resultado = TAA_MR / TAA_PR 
resultado

3.500000000000001

Como el proyecto que entra al portafolio con clasificación muy rentable (MR) permanece dentro del portafolio en promedio 3.5 veces mas que un proyecto que entra al portafolio con clasificación poco rentable (PR) se puede concluir que la afirmación del analista es incorrecta

Literal c

EL literal c pide calcular la ganancia neta promedio que generara un proyecto que inicialmente es clasificado como muy rentable. Para esto se definira un vector de costos que corresponden al ingreso o costo de acuerdo a la lista de estados que se definio previamente.

c = [38, 76, 114, 152, 190]

Como se sabe, las ganancias se calculan como la diferencia entre los ingresos y los costos. Para estimar los ingresos, se multiplicará el ingreso asociado a cada estado transitorio por el tiempo esperado que el proyecto permanece en dicho estado antes de ser absorbido. Esto se ve de la siguiente manera

\[ \text{Ingresos} = \$38 \cdot (\mathbb{I} - \mathbb{Q})^{-1}_{[MR,PR]} + \$76 \cdot (\mathbb{I} - \mathbb{Q})^{-1}_{[MR,R]} + \$114 \cdot (\mathbb{I} - \mathbb{Q})^{-1}_{[MR,MR]} \]

Por otro lado, para calcular los costos se debe tener en cuenta la probabilidad de que el proyecto sea liquidado o deje de ser rentable. Dichas probabilidades se pueden calcular mediante las probabilidades antes de absorción dadas por la siguiente matriz:

\[ (\mathbb{I} - \mathbb{Q})^{-1} \cdot \mathbb{R} \]

Teniendo en cuenta esto se pueden estimar los costos como:

\[ \text{Costos} = \$152 \cdot (\mathbb{I} - \mathbb{Q})^{-1} \cdot \mathbb{R}_{[MR,NR]} + \$190 \cdot (\mathbb{I} - \mathbb{Q})^{-1} \cdot \mathbb{R}_{[MR,L]} \]

Dicho procedimiento se presenta a continuación:

# Calculando los ingresos
ingresos = sum(c[i] * P.absorbtion_times(start = estados.index('MR'), target = i) for i in range(0, 3))

# Calculando los costos
costos = sum(c[i] * P.absorbtion_probabilities(start = estados.index('MR'), target = i)[0][0] for i in range(3, 5))

ganancias = ingresos - costos
ganancias

1120.0000000000005