setwd('~/PYE1213')
library(markovchain)

## Package:  markovchain
## Version:  0.8.5-2
## Date:     2020-09-07
## BugReport: https://github.com/spedygiorgio/markovchain/issues

library(prettydoc)

Introducción a los procesos estocástico

En la teoría de la probabilidad, un proceso estocástico es un concepto matemático que sirve para usar magnitudes aleatorias que varían con el tiempo o para caracterizar una sucesión de variables aleatorias (estocásticas) que evolucionan en función de otra variable, generalmente el tiempo. Cada una de las variables aleatorias del proceso tiene su propia función de distribución de probabilidad y pueden o no estar correlacionadas entre sí.

Cada variable o conjunto de variables sometidas a influencias o efectos aleatorios constituye un proceso estocástico. Un proceso estocástico \(xt\) puede entenderse como una familia uniparamétrica de variables aleatorias indexadas mediante el tiempo t. Los procesos estocásticos permiten tratar procesos dinámicos en los que hay cierta aleatoriedad.

Ejemplos

Los siguientes son ejemplos dentro del amplio grupo de las series temporales: • señales de telecomunicación; • señales biomédicas (electrocardiograma, encefalograma, etc.); • señales sísmicas; • el número de manchas solares año tras año; • el índice de la bolsa segundo a segundo; • la evolución de la población de un municipio año tras año; • el tiempo de espera en la cola de cada uno de los usuarios que van llegando a una ventanilla; • el clima, un gigantesco conjunto de procesos estocásticos interrelacionados (velocidad del viento, humedad del aire, etcétera) que evolucionan en el espacio y en el tiempo; • los procesos estocásticos de orden mayor a uno, como el caso de una serie de tiempo de orden 2 y una correlación de cero con las demás observaciones.

Cádenas de markov

Una cadena de Markov es una serie de eventos, en la cual la probabilidad de que ocurra un evento depende del evento inmediato anterior. En efecto, las cadenas de este tipo tienen memoria, “Recuerdan” el último evento y esto condiciona las posibilidades de los eventos futuros.

Esta dependencia del evento anterior distingue a las cadenas de Markov de las series de eventos independientes, como tirar una moneda al aire o un dado. En los negocios, las cadenas de Markov se han utilizado para analizar los patrones de compra,los deudores morosos, para planear las necesidades de personal y para analizar el reemplazo de equipo.

-> Ejemplo: Suponga que la posibilidad que llueva mañana depende de las condiciones del estado del clima de hoy. No importa las condiciones de los días anteriores, solo del estado del clima de hoy.

Suponga también que si llueve hoy, entonces lloverá mañana con una probabilidad α, y si no llueve hoy, entonces lloverá mañana con una probabilidad β.

Ejercicio de cadena de markov en R

documentación del paquete markovchain: https://cran.r-project.org/web/packages/markovchain/markovchain.pdf

Esta libreria pretende proveer objetos para realizar analisis estadísticos de cadenas de markov a tiempos discretos. Asumamos que tenemos una cadena de markov X={X1,X2,…} definida en el espacio de estados S={a,b,c} y cuya matriz de transición es:

\[ P = \left( {\begin{array}{ccc} 0 & 0.5 & 0.5 \\ 0.5 & 0 & 0.5 \\ 0.5 & 0.5 & 0 \\ \end{array} } \right)\]

Dicha cadena podemos crearla en R, de la siguiente forma:

Crear la matriz de transicion P:

P = matrix(c(0,0.5,0.5,.5,0,.5,.5,.5,0),nrow = 3,byrow = TRUE) 
P

##      [,1] [,2] [,3]
## [1,]  0.0  0.5  0.5
## [2,]  0.5  0.0  0.5
## [3,]  0.5  0.5  0.0

El argumento “nrows” de la funcion matrix es para declarar el numero de filas que deseamos que nuestra matriz P posea, y el argumento “byrows” es para que almacene los elementos de la matriz almacenados en c(), fila por fila.

Crear la matriz de transición creamos el objeto “markovchain” de la siguiente forma:

mc = new("markovchain",transitionMatrix=P,states=c("a","b","c"),name="Cadena 1")

La estructura del objeto mc (cadena de markov) esta dad por str:

str(mc)

## Formal class 'markovchain' [package "markovchain"] with 4 slots
##   ..@ states          : chr [1:3] "a" "b" "c"
##   ..@ byrow           : logi TRUE
##   ..@ transitionMatrix: num [1:3, 1:3] 0 0.5 0.5 0.5 0 0.5 0.5 0.5 0
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3] "a" "b" "c"
##   .. .. ..$ : chr [1:3] "a" "b" "c"
##   ..@ name            : chr "Cadena 1"

Resumen de la cadena de markov 1:

summary(mc)

## Cadena 1  Markov chain that is composed by: 
## Closed classes: 
## a b c 
## Recurrent classes: 
## {a,b,c}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

Para visualizar la transición de la cadena, utilizamos el comando plot:

plot(mc)

Otras funciones importantes son:

absorbingStates(): Identifica los estados Absorbentes
transientStates(): Identifica los estados Transitorios
recurrentClasses(): Identifica las clases recurrentes

Para la cadena de markov definida se obtiene que:

recurrentClasses(mc)

## [[1]]
## [1] "a" "b" "c"

transientStates(mc)

## character(0)

absorbingStates(mc)

## character(0)

Análisis probabilístico

Para conocer la probabilidad de transición en 1 paso entre un estado y otro basta con utilizar la función transitionProbability(), con los argumentos:

object: la cadena de markov
t0: el estado en el tiempo 0
t1: el estado en el tiempo 1

La probabilidad de transición en un paso del estado “a” al estado “c” es:

transitionProbability(object = mc , t0="a", t1="c")

## [1] 0.5

Recuerde que dicha probabilidad es un elemento de la matriz de transición P, por lo tanto, la probabilidad de transición del estado “a” al estado “b” es simplemente P23

mc[2,3]

## [1] 0.5

Es posible computar la matriz de transición en n pasos, simplemente computando la n-ésima potencia de la matriz de transición P, como ejemplo calcularemos la matriz de transición en n = 5 pasos.

n = 5 #el numero de pasos al futuro
mc^5

## Cadena 1^5 
##  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.31250 0.34375 0.34375
## b 0.34375 0.31250 0.34375
## c 0.34375 0.34375 0.31250

También se pueden conocer la distribución de la cadena en n pasos adelante (P(Xn)) multiplicando la distribución inicial de X0 por la matriz de transición en n pasos (Pn), calcule la distribución de la cadena en el tiempo n = 6, si la distribución inicial de la cadena es “(0.5, 0.2, 0.3)”.

x0 <- c(0.5,0.2,0.3) # La distribución de X en t = 0
n = 6
Xn = x0*(mc^n)
Xn

##              a       b         c
## [1,] 0.3359375 0.33125 0.3328125

Puesto que Xn es una función de densidad, la suma de las probabilidades en todos los estados debe ser 1.

sum(Xn)

## [1] 1

Finalmente encontrar la distribución estacionaria de la cadena se obtiene mediante la función “steadyStates” de la siguiente forma:

DistEst <- steadyStates(mc)
DistEst

##              a         b         c
## [1,] 0.3333333 0.3333333 0.3333333

Recuerde que los tiempos medio de recurrencia son los inversos multiplicativos de la distribución estacionaria y pueden ser computados fácilmente.

M <- 1/DistEst
M

##      a b c
## [1,] 3 3 3

-> Asignación

Dibuje el diagrama de transición, determine la clase de comunicación de las siguientes cadenas de Markov, clasifique éstas como recurrentes o transitorias (20%), y encuentre la distribución estacionaria si existe (10%).

$$ P = 
\left( {\begin{array}{cccc}
1/2 & 1/2 & 0 & 0 \\ 
0 & 1/2 & 1/2 & 0 \\ 
0 & 1/2 & 1/2 & 0 \\ 
1/4 & 1/4 & 1/4 & 1/4 \\ 
\end{array} } \right) $$

P = matrix(c(1/2,1/2,0,0,0,1/2,1/2,0,0,1/2,1/2,0,1/4,1/4,1/4,1/4),nrow = 4,byrow = TRUE) 
P

##      [,1] [,2] [,3] [,4]
## [1,] 0.50 0.50 0.00 0.00
## [2,] 0.00 0.50 0.50 0.00
## [3,] 0.00 0.50 0.50 0.00
## [4,] 0.25 0.25 0.25 0.25

mc = new("markovchain",transitionMatrix=P,states=c("a","b","c","d"),name="Cadena 2")

La estructura del objeto mc (cadena de markov) esta dad por str:

str(mc)

## Formal class 'markovchain' [package "markovchain"] with 4 slots
##   ..@ states          : chr [1:4] "a" "b" "c" "d"
##   ..@ byrow           : logi TRUE
##   ..@ transitionMatrix: num [1:4, 1:4] 0.5 0 0 0.25 0.5 0.5 0.5 0.25 0 0.5 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:4] "a" "b" "c" "d"
##   .. .. ..$ : chr [1:4] "a" "b" "c" "d"
##   ..@ name            : chr "Cadena 2"

Resumen de la cadena de markov 1:

summary(mc)

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

Para visualizar la transición de la cadena, utilizamos el comando plot:

plot(mc)

recurrentClasses(mc)

## [[1]]
## [1] "b" "c"

transientStates(mc)

## [1] "a" "d"

absorbingStates(mc)

## character(0)

Análisis probabilístico

Para conocer la probabilidad de transición en 1 paso entre un estado y otro basta con utilizar la función transitionProbability(), con los argumentos:

object: la cadena de markov
t0: el estado en el tiempo 0
t1: el estado en el tiempo 1

La probabilidad de transición en un paso del estado “a” al estado “c” es:

transitionProbability(object = mc , t0="a", t1="c")

## [1] 0

Se puede ver que la probabilidad de que pase del estado “a” al estado “c” es nula, ahora del estado “b” al “c” (clases de transición)

transitionProbability(object = mc , t0="b", t1="c")

## [1] 0.5

Esta es la probabilidad de transición del estado “b” al estado “c”

mc[2,3]

## [1] 0.5

Es posible computar la matriz de transición en n pasos, simplemente computando la n-ésima potencia de la matriz de transición P, como ejemplo calcularemos la matriz de transición en n = 5 pasos.

n = 5 #el numero de pasos al futuro
mc^5

## Cadena 2^5 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c, d 
##  The transition matrix  (by rows)  is defined as follows: 
##            a         b         c            d
## a 0.03125000 0.5000000 0.4687500 0.0000000000
## b 0.00000000 0.5000000 0.5000000 0.0000000000
## c 0.00000000 0.5000000 0.5000000 0.0000000000
## d 0.03027344 0.4990234 0.4697266 0.0009765625

X0 <- c(0.1,0.2,0.3,0.4) # La distribución de X en t = 0
n = 6
Xn = X0*(mc^n)
Xn

##                a         b         c            d
## [1,] 0.007714844 0.4999023 0.4922852 9.765625e-05

Puesto que Xn es una función de densidad, la suma de las probabilidades en todos los estados debe ser 1.

sum(Xn)

## [1] 1

Finalmente encontrar la distribución estacionaria de la cadena se obtiene mediante la función “steadyStates” de la siguiente forma:

DistEst <- steadyStates(mc)
DistEst

##      a   b   c d
## [1,] 0 0.5 0.5 0

Recuerde que los tiempos medio de recurrencia son los inversos multiplicativos de la distribución estacionaria y pueden ser computados fácilmente.

M <- 1/DistEst
M

##        a b c   d
## [1,] Inf 2 2 Inf

Diagrama de transición.

plot(mc)

Clases recurrentes

recurrentClasses(mc)

## [[1]]
## [1] "b" "c"

Clases transitoria

transientStates(mc)

## [1] "a" "d"

Distribución estacionaria

DistEst <- steadyStates(mc)
DistEst

##      a   b   c d
## [1,] 0 0.5 0.5 0

Se puede ver que esta distribución es muy extraña, porque unicamente se distribuye en el estado “b” y “c”

-> Asignación:

Encontrar un ejemplo práctico en código implementado en R de la cadena de markov a algún problema en particular y explicarlo.

En un supuesto en el cual se analizan las clases sociales, se puede hacer una matriz que diga que si unos padres de 40 años de edad pertenencen a una determinada clase social, ¿Cuál es la probabilidad de que su hijo pertenezca a una determinada clase social a la edad de 40 años?

$$ P = 
 \left( {\begin{array}{ccc}
   0.7 & 0.21 & 0.09 \\
   0.12 & 0.7 & 0.18 \\
   0.01 & 0.5 & 0.49 \\
  \end{array} } \right)$$

Lo que se tiene que analizar son estos datos para poder sacar de vista desde un punto de vista sociológico.

P = matrix(c(.7,.21,.09,.12,.7,.18,.01,.5,.49),nrow = 3,byrow = TRUE) 
P

##      [,1] [,2] [,3]
## [1,] 0.70 0.21 0.09
## [2,] 0.12 0.70 0.18
## [3,] 0.01 0.50 0.49

mc = new("markovchain",transitionMatrix=P,states=c("Alta","Media","Baja"),name="Cadena 3")

En primera instancia se tiene 3 datos para las 3 clases, se puede ver de esta manera.

inicial <- c(.05,.6,.35)
round(inicial*mc^3,2)

##      Alta Media Baja
## [1,] 0.18  0.57 0.25

Se puede preguntar si una persona de clase baja, que probabilidad hay de que dentro de 3 cambie o se mantenga en la clase socioeconomica.

inicial <- c(0,0,1)
round(inicial * mc ^3,2)

##      Alta Media Baja
## [1,] 0.13   0.6 0.28

Esto significa que cada 100 personas de clase baja, en 3 generaciones, sólo 28 de sus descendientes seguirán en la clase baja, 13 estarán en clase alta y 6 en clase media.

inicial <- c(1,0,0)
round(inicial * mc ^3,2)

##      Alta Media Baja
## [1,]  0.4  0.42 0.18

Se puede ver que de cada 100 personas, en 3 generaciones, solo 4 de ellos seguirán en la clase alta, 42 estarán en clase media y 18 estarán en clase baja

inicial <- c(0,1,0)
round(inicial * mc ^3,2)

##      Alta Media Baja
## [1,] 0.19  0.57 0.23

Se puede ver que de cada 100 personas, en 3 generaciones, 57 seguirán en clase media, 19 estarán en clase alta y 23 estarán en clase baja.

str(mc)

## Formal class 'markovchain' [package "markovchain"] with 4 slots
##   ..@ states          : chr [1:3] "Alta" "Media" "Baja"
##   ..@ byrow           : logi TRUE
##   ..@ transitionMatrix: num [1:3, 1:3] 0.7 0.12 0.01 0.21 0.7 0.5 0.09 0.18 0.49
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3] "Alta" "Media" "Baja"
##   .. .. ..$ : chr [1:3] "Alta" "Media" "Baja"
##   ..@ name            : chr "Cadena 3"

summary(mc)

## Cadena 3  Markov chain that is composed by: 
## Closed classes: 
## Alta Media Baja 
## Recurrent classes: 
## {Alta,Media,Baja}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

plot(mc)