Simulación e Inferencia de procesos de Markov con

aplicaciones actuariales

F. Baltazar-Larios

20 de febrero de 2023
ii
Índice general

1. Conceptos básicos de Estadı́stica 1

1.1. Estimación máximo verosı́mil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2. Propiedades asintóticas de los estimadores . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. Procesos estocásticos 5

2.1. Definiciones elementales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.1. Clasificación general. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2. Cadenas de Markov a tiempo discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.1. Conceptos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.2. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2.3. Resultados importantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.4. Clases de comunicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.5. Periodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.6. Propiedad Fuerte de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.7. Tiempos de llegada y tiempos de absorción . . . . . . . . . . . . . . . . . . . . . . 10

2.2.8. Clasificación de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2.9. Distribuciones invariantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.10. Cadenas regulares, absorbentes y convergentes. . . . . . . . . . . . . . . . . . . . . 14

2.2.11. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3. Procesos de Saltos de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.1. Definición . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.2. Resultados elementales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

iii
iv ÍNDICE GENERAL

2.3.3. Generador infinitesimal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3.4. Distribuciones estacionarias y distribuciones lı́mite . . . . . . . . . . . . . . . . . . 24

2.3.5. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3. Simulación de procesos de Markov 31

3.1. Cadenas de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.1.1. Algoritmo para cadenas de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.1.2. Cálculo numérico de propiedades de cadenas de Markov . . . . . . . . . . . . . . . 33

3.2. Procesos de saltos de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.2.1. Algoritmo para PSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.2.2. Cálculo numérico de propiedades de PSM . . . . . . . . . . . . . . . . . . . . . . . 37

3.2.3. Ejemplo. Modelo SIR-PSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4. Inferencia sobre procesos de Markov 41

4.1. Cadenas de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.1.1. Estimación máximo verosı́mil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.1.2. Propiedades asintóticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2. Procesos de saltos de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2.1. Con información a tiempo continuo. . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2.2. Con información a tiempo discreto. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5. Ecuaciones diferenciales estocásticas (EDE) 65

5.1. Movimiento Browniano o proceso de Wiener . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.1.1. Simulación de trayectorias del movimiento Browninano . . . . . . . . . . . . . . . 66

5.1.2. Aproximando al movimiento Browniano por una caminata aleatoria . . . . . . . . 66

5.1.3. Puente Browniano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2. Ecuaciones diferenciales estocásticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.3. Procesos de difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5.3.1. Proceso de difusión ergódicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.3.2. Ecuaciones de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.3.3. Puentes de difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

ÍNDICE GENERAL v

5.4. Fórmula de Itô . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.4.1. Transformación de Lamperti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.4.2. Cociente de verosimilitud para procesos de difusión . . . . . . . . . . . . . . . . . . 77

5.5. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6. Inferencia estadı́stica para ecuaciones diferenciales estocásticas 79

6.1. Simulación de procesos de difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.1.1. Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.1.2. Método de Milstein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.1.3. Simulación de puentes de difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6.1.4. Estimadores Monte-Carloy técnicas de reducción de varianza . . . . . . . . . . . . 82

6.2. Estimación por máxima verosimı́litud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6.2.1. Modelo Continuo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.3. Caso discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

6.3.1. Inferencia por verosimilitud exacta . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

6.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

A. Simulación de variables aleatorias 87

B. Convergencia de variables aleatorias y el Teorema central del lı́mite. 89

C. Códigos en R 93

D. Códigos en Julia 95

Bibliografı́a 111
vi ÍNDICE GENERAL
Capı́tulo 1

Conceptos básicos de inferencia

estadı́stica por máxima verosimilitud

En este capı́tulo se estudian los conceptos básicos para realizar inferencia estadı́stica mediante el método
de máxima verosimilitud. La función de verosimilitud, el estimador de máxima verosimilitud y sus pro-
piedades como: insesgamiento, consistencia y eficiencia asintótica son estudiados en este capı́tulo. Estos
conceptos se utilizan en el Capı́tulo 4 para encontrar el estimador máximo verosı́mil de los procesos de
Markov con espacio de estados finito y probar sus propiedades estadı́sticas. Para este capı́tulo se utilizó
principalmente [? ].

1.1. Estimación máximo verosı́mil

Sea (X , A , Q) un modelo estadı́stico donde X es el espacio muestral, A una σ-álgebra de X y Q = {Qθ |θ ∈

Θ} una clase parametrizada de medidas de probabilidad sobre (X , A ). Se asume que el espacio parametral
es un subconjunto del espacio Euclidiano d-dimensional Θ ⊆ Rn . Adicionalmente consideremos (Ω, F )
un espacio medible y X : Ω −→ X una función medible y P = {Pθ |θ ∈ Θ} medidas de probabilidad sobre
(Ω, F ) tales que Qθ es la medida de probabilidad inducida por X.

Definición 1.1. Un estimador es una función θ̃ : X −→ Θ y representa la forma de estimar el parámetro

θ de X, apartir de una muestra observada de X.

Definición 1.2. La función de verosimilitud L(θ) = L(θ; x) es la función de θ que para cada x ∈ X está
dada por
L(θ) = L(θ; x) = fθ (x), θ ∈ Θ.
La función de log-verosimilitud está dada por

l(θ) = l(θ; x) = log(L(θ; x)).

Donde fθ (x) es la función de densidad de x inducida por Pθ .

Nuestro objetivo es encontrar un estimador de θ a partir de una muestra observada, la función de verosi-
militud es una herramienta muy importante para determinar buenos candidatos. Por ejemplo si tenemos
que la razón de verosimilitud L(θ2 )/L(θ1 ) = 10 para θ1 , θ2 ∈ Θ, esto significa que bajo Pθ2 la observación
x es 10 veces más probable que bajo la medida Pθ1 . En esta situación nosotros estamos más inclinados a

1
2 CAPÍTULO 1. CONCEPTOS BÁSICOS DE ESTADÍSTICA

creer que θ2 es mejor candidato que θ1 para representar al parámetro desconocido θ .

El ejemplo antes mencionado nos da la motivación para dar la siguiente definición.
Definición 1.3. Si θ̂ = θ̂(x) es tal que
l(θ) ≤ l(θ̂), para todo θ ∈ Θ,
entonces θ̂(x) es llamado el estimador máximo verosı́mil .
Mientras que la ecuación
∂l(θ)
=0
∂θ
∂l(θ)
es llamada la ecuación de verosimilitud. Suponiendo claramente que la derivada de arriba exista ∂θ
denota el vector de derivadas parciales de l(θ).

Una observación que vale la pena mencionar es que cuando el estimador máximo verosı́mil θ̂(x) se en-
cuentra en el interior de Θ, éste resuelve la ecuación de verosimilitud.
Dado que ahora tenemos un posible candidato para estimar el verdadero parámetro θ, una situación ópti-
ma serı́a cuando la distribución del estimador bajo la medida Pθ puede ser calculada. Esto permite hacer
declaraciones de la cercanı́a con respecto al verdadero parámetro.
Definición 1.4. Se define el error cuadrático medio (ECM (θ̃)) de un estimador θ̃ como
ECM (θ̃) = Eθ [(θ̃(X) − θ)2 ] = V arθ [θ̃(X)] + Eθ [(E[θ̃(X)] − θ)2 ]
donde Eθ y V arθ son la esperanza y varianza respecto a la medida de probabilidad Pθ respectivamente.

El error cuadrático medio es una medida de la cercanı́a que tiene el estimador con respecto al verdadero
parámetro.
Definición 1.5. Un estimador θ̃ : X −→ Θ se dice insesgado si
Eθ [θ̃(X)] = θ , para cada θ ∈ Θ.
Por otro lado si un estimador no es insesgado, diremos que es sesgado.

Los estimadores insesgados son preferidos respecto a los sesgados, pues su error cuadrático medio es
menor.
Ejemplo 1.1. Sea X1 , . . . , Xn variables aleatorias i.i.d. (independientes e idénticamente distribuidas) con
densidad
β exp(−βx), x > 0,
Donde β > 0 es el parámetro. La función de log-verosimilitud es
l(β) = n(log(β) − β x̄),
aquı́ el estimador máximo verosı́mil es β̂ = X̄ −1 .
Ya que X̄ tiene una distribución gamma, la media del estimador es
n
Eβ [β̂] = β 6= β,
n−1
De este modo β̂ no es un estimador insesgado. Si en lugar del parámetro β nosotros consideramos el
parámetro 1/β , encontramos que el estimador máximo verosı́mil del nuevo parámetro es
1/β̂ = X̄.
Claramente E1/β [1/β̂] = 1/β, entonces X̄ es un estimador insesgado respecto al nuevo parámetro.
1.2. PROPIEDADES ASINTÓTICAS DE LOS ESTIMADORES 3

El ejemplo anterior muestra que el estimador máximo verosı́mil no es en general insesgado. También este
ejemplo muestra que el insesgamiento se relaciona tanto con el estimador y la parametrización en cuestión,
por tanto algunas parametrizaciones son más útiles que otras.
Definición 1.6. La variable aleatoria
∂l(θ, X)
U (θ) = ,
∂θ
es llamada la función score. La matriz
∂ 2 l(θ, X)
j(θ) = − ,
∂θ∂θT
es llamada la información observada, mientras que la esperanza de j(θ)

i(θ) = Eθ [j(θ)],

es llamada la información de Fisher o información esperada, esta matriz es importante porque bajo ciertas
condiciones los estimadores cumplen
d
(θ̃(X) − θ)i(θ)1/2 −
→ N (0̄, In )

donde In denota la matriz identidad de dimensiones n × n.

1.2. Propiedades asintóticas de los estimadores

Intuitivamente mientras más grande sea la muestra, es posible tener una noción más precisa del valor
del verdadero parámetro. Por esta razón es de vital importancia estudiar propiedades de los estimadores
cuando el tamaño de la muestra crece.
Una propiedad deseable es que la esperanza del estimador converja al verdadero parámetro cuando el
tamaño de la muestra tienda a infinito.
Definición 1.7. Un estimador θ̃ : X −→ Θ se dice asintóticamente insesgado si

Eθ [θ̃(X)] → θ, para cada θ ∈ Θ.

La propiedad anterior es más débil que la propiedad de insesgamiento, esto nos permite estudiar una gran
cantidad de estimadores, ya que algunos estimadores pueden ser asintóticamente insesgados pero no son
insesgados (ver Ejemplo 1.1).
La siguiente propiedad es muy importante para el estudio de propiedades asintóticas de los estimadore.
Definición 1.8. Un estimador θ̃ : X −→ Θ se dice consistente si para cada θ ∈ Θ
c.s.
θ̃(X) −−→ θ.

Hasta ahora sólo hemos analizado propiedades de la media de los estimadores la cual nos permite medir el
promedio del valor del estimador, sin embargo debemos analizar la varianza de este para poder determinar
qué tan alejados en promedio nos encontramos del verdadero parámetro.
Definición 1.9. Decimos que un estimador θ̃ : X −→ Θ es asintóticamente eficiente si
d
(θ̃(X) − θ)i(θ)1/2 −
→ N (0̂, In ).

Donde In denota la matriz identidad de dimensiones n × n

4 CAPÍTULO 1. CONCEPTOS BÁSICOS DE ESTADÍSTICA

Para una matriz semi-definida positiva B, nosotros denotamos por B 1/2 su correspondiente raı́z cuadrada.
La propiedad anterior muestra que nosotros podemos aproximar la distribución del máximo verosı́mil
mediante una normal multivariada con media θ y matriz de covarianzas i(θ)−1 cuando el tamaño de la
muestra es grande. Un ejemplo de estimador que cumpla ser asintóticamente eficiente es el estimador de
1/β del Ejemplo 1.1. La importancia de esta propiedad radica en que si la muestra es grande podemos
calcular la probabilidad de que nuestro verdadero parámetro θ se encuentre en cierta región de Rd .
Capı́tulo 2

Procesos Estocásticos con espacio de

estados a lo más numerable

En este capı́tulo se presenta una introducción breve al tema de los procesos estocásticos con espacio
de estados a lo más numerable. Se explican algunos conceptos y propiedades generales de este tipo de
procesos y se estudian brevemente algunos ejemplos particulares de procesos estocásticos. El contenido
de este capı́tulo será usado en los capı́tulos tres y cuatro, en donde consideraremos procesos estocásticos
con espacio de estados finito. La bibligrafı́a básica para este capı́tulo es principalmente [? ], [? ] y [? ].

2.1. Definiciones elementales

En esta sección se presentan los resultados fundamentales de procesos estocásticos. Se describen sus
caracterı́sticas principales, se realiza una clasificación básica y se dan algunos ejemplos.
Comenzaremos dando una definición formal de un proceso estocástico.

Definición 2.1. Un proceso estocástico X = {Xt }t∈T es una colección de variables aleatorias (definidas
en el mismo espacio de probabilidad) que toman valores en un conjunto E con parámetro T , es decir,
para cada t ∈ T , Xt es una variable aleatoria.

En esta definición a T se le conoce como el espacio parametral y es interpretado como el tiempo, de

esta forma a Xt se le entenderá como el estado del proceso al tiempo t. Además se considera que las
variables aleatorias Xt para toda t ∈ T , están definidas sobre el mismo espacio de probabilidad (Ω, F , P).
Un proceso estocástico puede ser interpretado como la evolución en el tiempo de algún fenomeno cuya
dinámica se rige por el azar.

Siguiendo la definición de proceso estocástico, se le puede ver como una función de dos variables

X : T × Ω → E,

tal que a cada pareja (t, ω) ∈ (T , Ω) se le asocia X(t, ω) = Xt (ω) y de esta manera para cada t ∈ T la
función de
ω → Xt (ω)

5
6 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

es una variable aleatoria y para cada ω ∈ Ω, la función

t → Xt (ω)

es una trayectoria o realización del proceso.

2.1.1. Clasificación general.

Según las caracteristicas del espacio parametral se puede hacer la primera distinción entre procesos es-
tocásticos.

Definición 2.2. Proceso estocástico a tiempo discreto Si T es un conjunto a lo más numerable, se

dice que X es un proceso estocástico a tiempo discreto, para nosotros será T = M, donde M ⊂ N y en
general se denota por X = {Xn }n≥0 .

Definición 2.3. Proceso estocástico a tiempo continuo Cuando T sea un conjunto infinito no
numerable, se dice que X es un proceso estocástico a tiempo continuo, para este caso T = [0, s] ó [0, ∞).
y en será denotado por X = {Xt }t≥0 .

En cuanto al conjunto donde las variables aleatorias toman valores, se le conoce como espacio de estados
del proceso y de manera análoga este también puede ser discreto o continuo.

2.2. Cadenas de Markov a tiempo discreto

Para esta sección de cadenas de Markov se utilizaron los libros [? ] y [? ], en donde se pueden consultar
las demostraciones de los resultados aquı́ enunciados.

2.2.1. Conceptos básicos

Definición 2.4. Consideremos una colección de variables aleatorias X0 , X1 , . . . , interpretando a Xn

como el estado de un sistema al tiempo n y supongamos que el conjunto de posibles valores de Xn es a lo
más numerable, enumerándolos, tenemos que E = {1, 2, . . . , N } para el caso finito y E = {1, 2, . . .} para
el caso numerable, E es llamado espacio de estados. Si este proceso satisface la propiedad de Markov

P[Xn+1 = xn+1 |Xn = xn , . . . , X0 = x0 ] = P[Xn+1 = xn+1 |Xn = xn ]

con x0 , . . . , xn+1 ∈ E, decimos que {Xn }n≥0 es una cadena de Markov.

A la distribución de X0 se le llama distribución inicial y la denotaremos por π = (π1 , . . . , πN ) , es decir,

si πi = P[X0 = i] para toda i ∈ E.

Algunas observaciones importantes:

1. Si Xn = xn , . . . , X0 = x0 entonces decimos que se ha presentado la trayectoria x0 , . . . , xn .

2.2. CADENAS DE MARKOV A TIEMPO DISCRETO 7

2. A la familia P[Xn+1 = j|Xn = i] se le conoce como familia de probabilidades de transición de la

cadena de Markov y describe la evolución de la cadena en el tiempo.

Definición 2.5. Se dice que la cadena de Markov {Xn }n≥0 con espacio de estados E y distribución inicial
π es homogénea si
P[Xn+1 = j|Xn = i] = P[X1 = j|X0 = i]
para cada i, j ∈ E y para cada n ≥ 0. En este caso definimos para cada i, j ∈ E

pij = P[X1 = j|X0 = i]

y para m ≥ 1
(m)
pij = P[Xm+n = j|Xn = i],

que representa la probabilidad de pasar del estado i al estado j en m unidades de tiempo.

Definición 2.6. Matriz de transición:

Es la matriz cuadrada de orden N formada por las probabilidades de transición

 
p11 p12 . . . p1N
 p21 p22 . . . p2N 
P = .
 
.. .. .. 
 .. . . . 
pN 1 pN 2 . . . pN N

con entradas no negativas y la suma de cada renglón es uno (matriz estocástica).

2.2.2. Ejemplos

Veamos algunos ejemplos de cadenas de Markov

Ejemplo 2.1. Cadena de rachas.

Considere la observación de un experimento con propabilidad p de ocurrencia, definiendo como transición

un ensayo y a Xn = como los estados de estar en racha de éxitos o no, entonces X = {Xn }n≥0 define una
cadena de Markov.

Ejemplo 2.2. Caminata aleatoria simple.

Un individuo está parado e inica una caminata en el cual sólo se mueve hacia adelante con probabilidad p
o hacia atrás con probabilidad 1 − p. Definiendo a Xn como la posición del individuo después de n−pasos,
se tiene entonces que E = Z y las probabilidades de transición para cualquier n ≥ 0 e i, j ∈ E están dadas
por:


p
 j =i+1
pij = 1 − p j = i − 1

0 e.o.c.

8 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Ejemplo 2.3. Problema del jugador.

Suponga que un jugador A apuesta $1 sucesivamente contra otro jugador B. A inicia con k unidades y
B con N − k. En cada apuesta, A gana con probabilidad p y pierde con probabilidad 1 − p. Si definimos
Xn como la fortuna de A al tiempo n, entonces X = {Xn }n≥0 define una cadena de Markov.
Ejemplo 2.4. Cadena de Ehrenfest.

Se tienen dos urnas, entre ambas hay un total de k bolas. En cada paso se elige una de las bolas al azar y
se cambia de urna. Si Xn = número de bolas en la urna A después de n ensayos, entonces X = {Xn }n≥0
define una cadena de Markov.

2.2.3. Resultados importantes

A continuación presentaremos algunos resultados importantes.

Proposición 2.1. Ecuación de Chapman-Kolmogorov.

Sea X = {Xn }n≥0 una cadena de Markov con matriz de transición

P = {pij }ij∈E entonces

(n+m) (n) (m)

X
pij = pik pkj .
k∈E

(n)
Este resultado nos dice que pij es la entrada ij de la n−ésima potencia de P. En general, no es posible
calcular explı́citamente las potecias de la matriz de transición. Sin embargo, nos podemos auxiliar de
paquetes computacionales para esta tarea.
Proposición 2.2. Propiedad de corrimiento.

Sea X = {Xn }n≥0 una cadena de Markov con espacio de estados E y con matriz de transición P . Sean
n, k ∈ N fijos y x0 , x1 , . . . , xn , . . . , xn+k ∈ E entonces

P[Xn+k = xn+k , . . . , Xn+1 = xn+1 |Xn = xn , . . . , X0 = x0 ] = P[Xk = xn+k , . . . , X1 = xn+1 |X0 = xn ].

Este resultado es una generalización de la propiedad de Markov.

2.2.4. Clases de comunicación

Definición 2.7. Sea X = {Xn }n≥0 una cadena de Markov con espacio de estados E y con matriz de
transición P . Sean i e j dos estados de E, se dice que el estado j es accesible desde el estado i si existe
(n)
n ∈ N, tal que pij > 0 (i → j). Si i → j e j → i entonces, i y j se dice que se comunican, y cuando
esto ocurre escribimos i ↔ j.

El siguiente resultado nos muestra que la relación ↔ es de equivalencia.

Proposición 2.3. La relación x ↔ y (comunicación) es una relación de equivalencia, por lo tanto, induce
a una partición en el espacio de estados.
2.2. CADENAS DE MARKOV A TIEMPO DISCRETO 9

A las clases de equivalencia inducidas por esta relación les llamaremos clases de comunicación.
Definición 2.8. Sea X = {Xn }n≥0 una cadena de Markov con espacio de estados E y con matriz de
transición P , decimos que la cadena es irreducible si para cada i, j ∈ E tenemos que i ↔ j.

Sea i ∈ E, a su clase de cominicación la dentaremos por Ci , es decir

Ci = {j ∈ E : i ↔ j}.

En cuanto a la partición inducida por la comunicación de estados decimos que, tomando C ⊂ E , es un

conjunto cerrado si para toda j ∈ E − C, i ∈ C no puede acceder a j, en caso contrario decimos que es
una clase abierta.

2.2.5. Periodo

Definición 2.9. El periodo de un estado i ∈ E es un número entero no negativo denotado por d(i), y
definido como sigue:
(n)
d(i) = m.c.d{n ≥ 1 : pii > 0},
(n)
en donde m.c.d significa máximo común divisor. Cuando pii = 0 para toda n ≥ 1, se define d(i) = 0. En
particular, se dice que un estado i ∈ E es aperiódico si d(i) = 1. Cuando d(i) = k ≥ 2 se dice que i es
periódico de periodo k.
Proposición 2.4. Si los estados i, j ∈ E pertenecen a la misma clase de comunicación, entonces tienen
el mismo periodo.

2.2.6. Propiedad Fuerte de Markov

Una extensión de la propiedad de Markov es la presentada a tiempos aleatorios conocida como propiedad
fuerte de Markov. Comenzaremos esta sección dando la definición de tiempos aleatorios.
Definición 2.10. Un tiempo aleatorio es una variable aleatoria

τ : Ω → N ∩ {∞}.

Definición 2.11. Sea X = {Xn }n≥0 una cadena de Markov con espacio de estados E y con matriz de
transición P . Un tiempo aleatorio τ es un tiempo de paro respecto a la cadena X = {Xn }n≥0 si para
n ∈ N existe An ⊂ E n+1 tal que

{τ = n} = {(X0 , . . . , Xn ) ∈ An }.

A la variable aleatoria tiempo de paro la podemos interpretar como el tiempo que se obtiene de observar
una trayectoria hasta que se cumpla cierta condición.
Teorema 2.1. Propiedad Fuerte de Markov. Sea {Xn }n≥0 una cadena de Markov con matriz de
transición P y τ un tiempo de paro. Entonces condicionado sobre {τ < ∞} ∩ {Xτ = i} el proceso
{Xτ +m }m≥0 es una cadena de Markov que comienza en i con matriz de transición P y es independiente
de las variables X0 , X1 , . . . , Xτ .
10 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

2.2.7. Tiempos de llegada y tiempos de absorción

En el estudio de las cadenas de Markov resulta de particular interés el estudio de los tiempos aleatorios
de la ocurrencia de eventos relacionados con la llegada a un estado o un conjunto de estados del espacio
de la cadena. En está sección estudiaremos este tipo de variables aleatorias.
A la m-ésima vez que la cadena accede a A un subconjunto del espacio de estados,
Definición 2.12. Sea τm
es decir (
A ∞ Xn ∈ E − A para toda n
τm = A
mı́n{n > τm−1 : Xn ∈ A} e.o.c
donde se define τ0A = 0
A < ∞ para cada m y por el Teorema 2.1 tenemos que {τ A }
Nota 2.1. Bajo la hipótesis de que τm m m≥1 son
i.i.d.

Veamos un ejemplo donde podemos ilustrar esta definición.

Ejemplo 2.5. Racha de éxitos. Sea E1 , E2 , . . . , una sucesión de ensayos Bernoulli con probabilidad de
éxito p y probabilidad de fracaso q = 1 − p. Si Xn es el número de éxitos consecutivos previos al tiempo n
(incluido el tiempo n). Se dice que una cadena de r éxitos ocurre al tiempo n si en el ensayo n − r ocurre
un fracaso y los resultados de los ensayos n − r + 1 al n son éxitos.
(n)
Definición 2.13. Para cada n ≥ 1, denotamos por fij a la probabilidad de que una cadena que inicia
en el estado i, llegue al estado j por primera vez en exactamente n pasos, es decir,
(n)
fij = P[Xn = j, Xn−1 6= j . . . , X1 6= j, X0 = i).

Además definimos (
(0) 1 si i = j
fij =
0 e.o.c

Para la cadena de Markov de racha de éxitos tenemos (por ejemplo)

(n)
1. f01 = q n−1 p.
(n)
2. f00 = pn−1 q.

Un resultado importante es el que permite descomponer la probabilidad de visitar el estado j iniciando

en i en términos de todas las posibilidades de tiempos que pueden ser la primera visita.
Proposición 2.5. Para cada n ≥ 1,
n
(n) (k) (n−k)
X
pij = fij pij .
k=1

2.2.8. Clasificación de estados

De acuerdo a las propiedades que tiene cada estado en una cadena de Markov en cuanto a la probabilidad
de regresar a él, los estados pueden ser clasificados de la siguiente manera.
2.2. CADENAS DE MARKOV A TIEMPO DISCRETO 11

Definición 2.14. Sea X = {Xn }n≥0 una cadena de Markov con espacio de estados E y con matriz de
transición P . Se dice que un estado i es recurrente si la probabilidad de (eventualmente) regresar a i,
es uno, es decir, si
P[τ1i < ∞|X0 = i] = 1.
Un estado que no es recurrente se llama transitorio.

Tratar de clasificar los estado en recurrentes y transitorios partiendo de la definición no resulta, en general,
tarea fácil. Una manera distinta de enunciar estos conceptos es la siguiente forma.

Definición 2.15. Definimos a

∞
(n)
X
fi = fii ,
n=1

que es la probabilidad de un eventual regreso al estado i.

Nota 2.2. Decimos que:

(1)
El estado i es recurrente si fi = 1, en particular si fxx = 1 el estado se le llama absorbente.

El estado x es transitorio si fx < 1.

Definición 2.16. Definiendo al número de visitas al estado j hasta el tiempo n como la variable aleatoria

n
X
Vij (n) = I{Xk =j} , cuando X0 = i
k=1

Proposición 2.6.

P (n)
1. Un estado i es recurrente si y sólo si n pii = ∞.
P (n)
2. Un estado i es transitorio si y sólo si n pii < ∞.

Este resultado nos proporciona una equivalencia, en términos de las matrices de transición, para que un
estado sea recurrente.

A continuación presentamos una serie de resultados respecto a la clasificación de estados y la comunicación.

Proposición 2.7. Sean i, j ∈ E tales que se comunican, si i es transitorio entonces j también es transi-
torio.

Este resultado tiene como consecuencia que la transitoriedad o recurrencia son propiedades de clase.

Un criterio fácil para saber si una clase es transitoria es el siguiente.

Proposición 2.8. Sea C ⊂ E una clase abierta. Entonces C es transitoria.

Usando estos resultados tenemos la siguiente conclusión.

12 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Proposición 2.9. Si el espacio de estados es finito, una clase es recurrente si y sólo si es cerrada.

Corolario 2.1. En una cadena de Markov irreducible con espacio de estados finito todos sus estados son
recurrentes.

Corolario 2.2. Toda cadena de Markov con espacio de estados finito tiene por lo menos un estado
recurrentes.

Nota 2.3. En general, el análisis de recurrencia y transitoriedad, en cadenas con espacio de estados
infinito, es más complicado.

Siguiendo la clasificación de estados, resulta importante conocer el número esperado de pasos para regresar
a i, es decir,
∞
(n)
X
µi = nfii ,
n=1

a la cual llamaremos tiempo medio de recurrencia.

Nota 2.4. Puede darse el caso de que siendo i un estado recurrente el tiempo medio de recurrencia, µi ,
sea infinito; siendo éste el caso en el que aunque el regreso a i es cierto se necesita un tiempo infinito. Ası́,
realizamos la siguiente distinción entre los estados recurrentes obteniendo una subclasificación de estos.

Definición 2.17. Decimos que

1. El estado i es recurrente positivo si µi < ∞.

2. El estado i es recurrente nulo si µi = ∞.

Se tiene directamente el siguiente resultado.

Proposición 2.10. En una cadena irreducible con espacio de estados finito, todos los estados son recu-
rrentes positivos.

Ejemplo 2.6. Veamos que la cadena de Markov de racha de éxitos es una cadena recurrente positiva.

∞ ∞ ∞
X (n)
X X 1
µ0 = nf00 = n(1 − p)pn−1 = (1 − p) npn−1 = < ∞.
1−p
n=1 n=1 n=1

Otro resultado de gran utilidad es el siguiente.

Proposición 2.11. No existen estados recurrentes nulos en cadenas de Markov finitas.

El siguiente resultado se sigue de la teorı́a que los promedios temporales son iguales a los promedios
espaciales.

Teorema 2.2. Teorema ergódico para cadenas de Markov Sea X = {Xn }n≥0 una cadena de Markov
con espacio de estados E y con matriz de transición P . Sean i, j cualesquiera estados de una cadena de
Markov irreducible se cumple que
Vij (n) 1
lı́m =
n→∞ n µj
casi seguramente.
2.2. CADENAS DE MARKOV A TIEMPO DISCRETO 13

Ejemplo 2.7. Considere una cadena de Markov con espacio de estados {0, . . . , 5} y matriz de transición

 1 1 
2 2 0 0 0 0
1 2

3 3 0 0 0 0 
1 7
 
 0 0 0 0 
P = 8 8 
 1 1
0 0 1 1 
 4 4 4 4 
 0 3 1
0 4 0 4 0 
1 1 1 2
0 5 0 5 5 5
Determine que estados son transitorios y cuales recurrentes. Calcular y estimar tiempos medios de recu-
rrencia y número esperado de visitas a cada estado.

2.2.9. Distribuciones invariantes

Definición 2.18. Sea ν ∈ R|E| , decimos que ν es una distribución invariante o estacionaria de la
cadena de Markov con matriz de transición P , si satisface

1. Tiene todas sus entradas no negativas.

P
2. i∈E νi = 1.
P
3. j∈E νi pij = νj .

En términos matriciales, la distribución de probabilidad ν es invariante si

νP = ν.

Nota 2.5. La interpretación de una distribución invariante es que una variable aleatoria con dicha
distribución mantendrá su distribución en la evolución de la cadena de Markov en el tiempo.

Desde que encontrar una distribución estacionaria en una cadena de Markov consiste en resolver un sis-
tema de |E| + 1 ecuaciones y |E| incógnitas podemos no tener solución, una infinidad de soluciones o una
única solución.

Los siguientes resultados son de gran utilidad en el cálculo de la distribución estacionaria de una cadena
de Markov.

Proposición 2.12. Sea ν una distribución estacionaria para una cadena de Markov. Si i es un estado
transitorio o recurrente nulo, entonces νi = 0.

Una propiedad importante de una distribución invariante es que si X0 tiene distribución ν, entonces Xn
tendrá la misma distribución para toda n.

Proposición 2.13. Toda cadena de Markov que es irreducible y recurrente positiva tiene una única
distribución estacionaria dada por
1
νi = > 0.
µi
Teorema 2.3. Para una cadena irreducible las siguientes condiciones son equivalentes.
14 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

1. Todos los estados son recurrentes positivos.

2. Algún estado es recurrente positivo.

3. La cadena tiene una única distribución invariante

Ejemplo 2.8.
 
1/2 1/2
P =
1/4 3/4

ν = (1/3, 2/3).

Ejemplo 2.9. Cadena de Ehrenfest. Como esta cadena de Markov es irreducible y con espacio de
estados finito se tiene que es recurrente positiva y por los resultados anteriores tiene una única distribución
invariante ν.

2.2.10. Cadenas regulares, absorbentes y convergentes.

En esta sección presentaremos resultados del comportamiento de una cadena de Markov recurrente cuando
la trayectoria es lo suficientemente grande. Otro concepto importante en el estudio de las cadenas de
Markov es el relacionado con la posible periodicidad con la que se visita a cada estado.

Definición 2.19. Decimos que una cadena de Markov {Xn } con espacio de estados finito es regular si
existe una potencia de la matriz de transición con todas las entradas positivas.

Ahora presentamos resultados fundamentales en el estudio de la convergencia de las cadenas de Markov.

Definición 2.20. Si los lı́mites

(n)
νj = limn→∞ pij

existen para cada estado j en el espacio de estados entonces π, es la distribución lı́mite.

Nota 2.6. Las siguientes observaciones son consecuencia de la definición anterior.

1. La distribución lı́mite es inducida por la matriz de transición P .

2. La distribución lı́mite no depende de la distribución inicial.

3. El lı́mite de las potencias de P es una matriz estocástica para la cual todos sus renglones son la
distribución lı́mite.

4. La distribución lı́mite es una distribución estacionaria.

Teorema 2.4. (Existencia de la distribución lı́mite.) Si la cadena de Markov es irreducible, aperiódica y

con distribución invariante ν, entonces
(n)
νj = lı́m pij
n→∞

Teorema 2.5. Si la cadena de Markov es irreducible, aperiódica y recurrente positiva existen las proba-
bilidades lı́mite y son la única distribución invariante.
2.2. CADENAS DE MARKOV A TIEMPO DISCRETO 15

2.2.11. Ejemplos

Ejemplo 2.10. Consideremos la cadena de Markov del Ejemplo 2.4, el espacio de estados es E =
{0, 1, 2, . . . , k} y la matriz de transición asociada es
 
0 1 0 0 ... 0 0
 1/k 0 (k − 1)/k 0 ... 0 0 
 
 0 2/k 0 (k − 2)/k . . . 0 0 
P = . ..  .
 
.. .. .. .. ..
 .. . . . . . . 
 
 0 0 0 0 . . . 0 1/k 
0 0 0 0 ... 1 0

Es fácil mostrar que esta cadena es irreducible, por lo tanto existe una única clase de comunicación

{0, 1, 2, . . . , k}.

Por la Proposición 2.4 sabemos que el periodo para todos los estados es el mismo, además sabemos
(n)
{n ≥ 1 : p00 > 0} = {n ≥ 1 : n es par},

es decir, tiene periodo 2.

Debido a que esta cadena es irreducible, por la Proposición 2.10 es recurrente positiva. Por lo tanto por
la Proposición 2.13 tiene una única distribución estacionaria.
Desarrollando el sistema de ecuaciones νP = ν, tenemos
1
ν0 = ν1
k
2
ν1 = ν0 + ν2
k
k−1 3
ν2 = ν1 + ν3
k k
..
.
2
νk−1 = νk−2 + νk
k
1
νk = νk−1 .
k

Reescribiendo cada ecuación en términos de ν0 obtenemos

 
νj = kj ν0 , para todo j ∈ E.

Sumando todas estas ecuaciones

k  
X
k
1 = ν0 j = ν0 2k ,
j=0

por lo que concluimos que   1

k
νj =j , para todo j ∈ E.
2k
Por la Proposición 2.13 podemos calcular los tiempos medios de recurrencia mediante la distribución
estacionaria
1 2k
µj = =   , para todo j ∈ E.
νj k
j
16 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

2.3. Procesos de Saltos de Markov

En esta sección consideraremos procesos de Markov a tiempo continuo, homogéneos y con espacio de
estados numerable, describiremos caracterı́sticas relevantes de dichos procesos que nos serán de utilidad
para simular trayectorias de ellos.

2.3.1. Definición

Definición 2.21. Un proceso de saltos de Markov (que abreviaremos como PSM) es un proceso de Markov
en tiempo continuo {Xt }t≥0 , con espacio de estados a lo más numerable E y que cumple las siguientes
condiciones:

1. Cada trayectoria t −→ Xt es continua por la derecha con lı́mite por la izquierda.

2. Para cualesquiera estados x0 , x1 , . . . , xn+1 y para cualesquiera ti < ti+1 = si+1 +ti con i ∈ {0, . . . , n},
con t0 = 0 se cumple:

P[Xtn+1 = xn+1 |Xtn = xn , . . . , Xt0 = x0 ] = P[Xtn+1 = xn+1 |Xtn = xn ]

= P[Xsn+1 = xn+1 |X0 = xn ].

Nota 2.7. En lo sucesivo supondremos que el espacio de estados E = {1, 2, . . . , N }.

Proposición 2.14. Un PSM {Xt }t≥0 es constante por pedazos.

Demostración. Sea s ≥ 0, consideremos una trayectoria fija del proceso, supongamos que Xs = x0 , como
el espacio de estado es finito, existe u > 0 tal que ningún estado en E salvo x0 se encuentra en el intervalo
(x0 − u, x0 + u). Como el proceso es continuo por la derecha entonces existe δ > 0 tal que el proceso en
[s, s + δ) se encuentra en el intervalo (x0 − u, x0 + u), pero como en este intervalo solo vive x0 , entonces
se concluye que el proceso es constante en [s, s + δ).

Definición 2.22. La función definida por

pxy (t) = P[Xt = y|X0 = x]

con t ≥ 0 y x, y ∈ E es llamada la función de transición.

Definición 2.23. Considere un PSM {Xt }t≥0 con espacio de estados E, se define la matriz de transición
al tiempo t ≥ 0 como

P (t) = {pxy (t)}x,y∈E .

Definición 2.24. A la familia {P (t)}t≥0 se le llama semigrupo de transición del proceso de Markov.

Proposición 2.15. El semigrupo de transición {P (t)}t≥0 de un PSM {Xt }t≥0 es estándar, es decir

lı́m P (t) = I,
t→0+

donde I es la matriz identidad de dimensión N × N.

2.3. PROCESOS DE SALTOS DE MARKOV 17

Demostración. Por la Proposición 2.14

c.s.
I{X0 =x,Xt =x} −−−→ I{X0 =x}
t→0+

y por el Teorema de la convergencia dominada (ver [? ] página 62)

lı́m E[I{X0 =x,Xt =x} ] = E[I{X0 =x} ].

t→0+

Concluimos que
(
1 si x = y,
lı́m pxy (t) =
t→0+ 0 si x 6= y.

2.3.2. Resultados elementales

Proposición 2.16. (Ecuación de Chapman-Kolmogorov) Considere un PSM {Xt }t≥0 con espacio
de estados E, entonces para cualesquiera par de numeros t y s positivos y para todo x, y ∈ E tenemos
X
pxy (t + s) = pxz (t)pzy (s).
z∈E

Debido a que un PSM es constante por pedazos, ahora nos interesa con qué frecuencia este proceso cambia
de estado.

Definición 2.25. Se define el tiempo de estancia en el estado Xt de un PSM {Xt }t≥0 como

∆(t) := ı́nf{s > 0 : Xt+s 6= Xt }

para todo t ≥ 0.

Nota 2.8. Esta variable aleatoria está bien definida porque el proceso es constante por pedazos.

Definición 2.26. Se define el tiempo de primer salto de un PSM {Xt }t≥0 como

T1 := ı́nf{t > 0 : Xt 6= X0 }

y el n-ésimo tiempo de salto como

Tn = ı́nf{t > Tn−1 ; Xt 6= XTn−1 } para todo n ≥ 1.

Por convención definimos T0 = 0.

Teorema 2.6. Considere un PSM {Xt }t≥0 con espacio de estados E, si en el tiempo t ≥ 0 ∆(t) < ∞
entonces para todo x ∈ E existe λx > 0, independiente de t, tal que

P[∆(t) > s|Xt = x] = exp(−λx s)

para toda s > 0.

18 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Demostración. Debido a que un PSM es homogéneo en el tiempo

P[∆(t) > s|Xt = x] = P[∆(0) > s|X0 = x],

definimos h(s) = P[∆(0) > s|X0 = x], entonces se tiene que

h(t + s) = P[∆(0) > t + s|X0 = x] = P[∆(0) > t, ∆(t) > s|X0 = x]

= P[∆(0) > t|X0 = x]P[∆(t) > s|∆(0) > t, X0 = x]
= h(t)P[∆(t) > s|Xt = x]
= h(t)h(s).

La función h es continua por la derecha, ya que es función de sobrevivencia, además 0 ≤ h(s) ≤ 1 para
toda s ≥ 0, lo cual implica que la única solución debe ser

h(s) = exp(−λx s),

donde λx > 0.

Ahora demostraremos un importante teorema que nos muestra la manera de como simular procesos de
tipo PSM.

Teorema 2.7. Considere un PSM {Xt }t≥0 con T1 < ∞, entonces para cualesquiera estados x0 , x1 , . . . , xk
y para cualesquiera tiempos ti < ti+1 = si+1 + ti con i ∈ {0, . . . , k − 1} y t0 = 0, x0 = y 6= x se tienen las
siguientes dos afirmaciones

1.
P[XT1 = y|T1 ≥ t, X0 = x]
no depende de t, es decir la distribución de salto es independiente del tiempo de primer salto.

2.
k−1
Y
P[XT1 +ti = xi , ∀i ≤ k|T1 ≥ t, X0 = x] = P[XT1 = y|T1 ≥ t, X0 = x] P[Xsj+1 = xj+1 |X0 = xj ],
j=0

es decir, el proceso después de T1 es independiente de T1 y se comporta en ley de la misma forma

que el proceso original, con distribución inicial P[XT1 = y|T1 ≥ t, X0 = x].

Demostración. Considere el proceso estocástico definido por Ynm = Xn/2m , con n ∈ {0, 1, . . .}, es fácil
darse cuenta que es de Markov homogéneo en el tiempo.
Se define

σm = (1/2m )mı́n{n > 0 : Ynm 6= X0 }.

Evidentemente se tiene que σm ≥ σm+1 ≥ T1 para toda m ∈ N, ya que

Y2m+1
m+1 σ
m
= X2m+1 σm /2m+1 = Xσm = Y2m
mσ
m
6= X0 ,
1
por lo que tenemos 2m+1 σm ≥ 2m+1 σm+1 para toda m ∈ N, además σm − T1 < 2m para toda m ∈ N.
Demostraremos la segunda afirmación del teorema.
2.3. PROCESOS DE SALTOS DE MARKOV 19

n
Denotemos a nm como el menor entero n tal que 2m ≥ t, debido a que {Xt }t≥0 es continuo por la derecha

P[XT1 +ti = xi , ∀i ≤ k, T1 ≥ t, X0 = x]

= lı́m P[Xσm +ti = xi , ∀i ≤ k, σm ≥ t, X0 = x]

m→∞
nm
= lı́m P[Xσm +ti = xi , ∀i ≤ k, σm ≥ m , X0 = x]
m→∞ 2
∞
X n
= lı́m P[Xσm +ti = xi , ∀i ≤ k, σm = m , X0 = x]
m→∞
n=n
2
m
∞
X
= lı́m P[Xn/2m +ti = xi , ∀i ≤ k, X0/2m = x, . . . , Xn−1/2m = x]
m→∞
n=nm

y por la propiedad de Markov tenemos

∞
X k−1
Y
= lı́m P[X0 = x]P[X1/2m = x|X0 = x]n−1 P[X1/2m = y|X0 = x] P[Xsj+1 = xj+1 |X0 = xj ]
m→∞
n=nm j=0
k−1
Y ∞
X
= lı́m P[X1/2m = y|X0 = x] P[Xsj+1 = xj+1 |X0 = xj ] P[X0 = x]P[X1/2m = x|X0 = x]n−1
m→∞
j=0 n=nm

por lo anterior
∞
X
P[T1 ≥ t|X0 = x] = lı́m (1 − P[X1/2m = x|X0 = x]) P[X1/2m = x|X0 = x]n−1
m→∞
n=nm

por lo tanto
P[XT1 +ti = xi , ∀i ≤ k|T1 ≥ t, X0 = x]

k−1
Y P[X1/2m = y|X0 = x]
= P[Xsj+1 = xj+1 |X0 = xj ] lı́m
m→∞ 1 − P[X1/2m = x|X0 = x]
j=0

tomando k = 0, se concluye

P[X1/2m = y|X0 = x]
P[XT1 = y|T1 ≥ t, X0 = x] = lı́m
m→∞ 1 − P[X1/2m = x|X0 = x]

esto demuestra la primer y segunda afirmación.

Definición 2.27. Dado un PSM {Xt }t≥0 con espacio de estados E, se define la matriz de transición
asociada a la distribución de salto como:

P = {pxy }x,y∈E = {P[XT1 = y|X0 = x]}x,y∈E .

Teorema 2.8. Considere un PSM {Xt }t≥0 con espacio de estados E, supongamos que T1 < . . . <
Tn < ∞, entonces para cualesquiera tiempos t1 , . . . , tn−1 , tn , 0 = s0 < . . . < sm y cualesquiera estados
x0 , x1 , . . . , xn = y0 , y1 , . . . , ym ∈ E.
Bajo estas condiciones se tienen las siguientes dos afirmaciones:
20 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

1. P[Tk − Tk−1 ≤ tk , XTk = xk , 1 ≤ k ≤ n, XTn +sj = yj , 0 ≤ j ≤ m|X0 = x0 ]

= P[Tk − Tk−1 ≤ tk , XTk = xk , 1 ≤ k ≤ n|X0 = x0 ]P[Xsk = yk , 0 ≤ k ≤ m|X0 = y0 ].

2. P[XTn +sj = yj , 0 ≤ j ≤ m|X0 = x0 , Tk − Tk−1 ≤ tk , XTk = xk , 1 ≤ k ≤ n]

= P[Xsk = yk , 0 ≤ k ≤ m|X0 = y0 ].

Demostración. El teorema se sigue de aplicar el Teorema 2.7 multiples veces.

Nota 2.9. Por el Teorema 2.6 y el Teorema 2.8 podemos concluir que para todo n ∈ N

1. El proceso después del tiempo Tn es condicionalmente independiente a la trayectoria antes de Tn

dado XTn y se distribuye igual que el proceso original empezando en el estado XTn .

2. Tn+1 − Tn condicionado a XTn = x se distribuye exponencial de parámetro λx , además es indepen-

diente a la distribución del siguiente salto.
Corolario 2.3. Considere un PSM {Xt }t≥0 con espacio de estados E, supongamos que T1 < . . . < Tn <
∞, entonces para cualesquiera tiempos t1 , . . . , tn y cualesquiera estados x0 , x1 , . . . , xn−1 ∈ E
n−1
Y
P[∆(T0 ) ≤ t1 , . . . , ∆(Tn−1 ) ≤ tn |XTk = xk , 0 ≤ k ≤ n − 1] = P[∆(0) ≤ tk+1 |X0 = xk ].
k=0

Definición 2.28. Dado un PSM {Xt }t≥0 , el Teorema 2.8 nos asegura que el proceso estocástico definido
como:

{Zn = XTn }n≥0 .

Conforma una cadena de Markov con matriz de transición dada por la Definición 2.27.

2.3.3. Generador infinitesimal

En esta parte del trabajo, nosotros caracterizaremos a un PSM mediante los cambios infinitesimales que
ocurren en intervalos de tiempo muy pequeños.
Considere un PSM {Xt }t≥0 con espacio de estados E, con Xt = x. Ahora, estamos interesados en analizar
que pasa en intervalos de tiempo infinitesimales (t, t + h), tenemos los siguientes casos.

1. El proceso puede estar en el mismo estado al tiempo t+h y esto ocurre con probabilidad pxx (h)+o(h),
donde o(h) representa la posibilidad que el proceso se salga y regrese al estado i en el intervalo de
tiempo.

2. El proceso se mueva al estado y en t + h y esto ocurre con probabilidad pxy (h) + o(h), donde o(h)
representa la posibilidad que el proceso realize dos o más transiciones en el intervalo de tiempo.

Ahora estamos interesados en el comportamiento de pxy (h) para h pequeña, para esto consideremos la
siguiente proposición.
2.3. PROCESOS DE SALTOS DE MARKOV 21

Proposición 2.17. Considere un PSM {Xt }t≥0 con espacio de estados E, entonces para cualesquiera
x, y ∈ E

pxy (h) − δxy

lı́m = qxy .
h→0+ h
Donde
(
1 para x = y
δxy =
0 para x 6= y.

y qxy = λx pxy si x 6= y y qxx = −λx .

Para probar esta proposición necesitamos el siguiente lema.

Lema 2.1. Considere un PSM {Xt }t≥0 con espacio de estados E, supongamos T1 < . . . < Tn < ∞.
Bajo estas condiciones

βt
sn−1 exp(−s)
Z
P[Tn ≤ t] ≤ ds.
0 (n − 1)!
Además si sumamos sobre n obtenemos
∞
X
P[Tn ≤ t] ≤ βt.
n=1

Donde β = máx{λx : x ∈ E}.

Demostración. Sabemos que

n
X λxk (Tj − Tj−1 )
P[Tn ≤ t|XTk = xk , 0 ≤ k ≤ n − 1] = P[ ≤ t|XTk = xk , 0 ≤ k ≤ n − 1]
λ xk
j=1
n
X λxk (Tj − Tj−1 )
≤ P[ ≤ t|XTk = xk , 0 ≤ k ≤ n − 1]
β
j=1
Xn
= P[ λxk (Tj − Tj−1 ) ≤ βt|XTk = xk , 0 ≤ k ≤ n − 1].
j=1

Por el Corolario 2.3

n
X
λxk (Tj − Tj−1 )|XTk = xk , 0 ≤ k ≤ n − 1
j=1

se distribuye Gamma de parámetros (n, 1), por lo que se tiene

Z βt n−1
s exp(−s)
P[Tn ≤ t] ≤ ds.
0 (n − 1)!
Sumando sobre n y aplicando Teorema de la convergencia dominada (ver [? ] página 62)
∞ ∞ n−1
Z βt X Z βt
X s exp(−s)
P[Tn ≤ t] ≤ ds = ds = βt.
0 (n − 1)! 0
n=1 n=1
22 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Definición 2.29. Dado un PSM {Xt }t≥0 se define la variable que cuanta los saltos hasta el instante t
como
Nt = máx{n ≥ 1 : Tn ≤ t}.

Prueba de la Proposición 2.17

Demostración. La demostración se divide en dos casos.

Caso 1. x 6= y.

Por el Lema 2.1

P[Xh = y|X0 = x] = P[Xh = y, Nh = 1|X0 = x] + P[Xh = y, Nh ≥ 2|X0 = x]

exp(λy −λx )h −1
= λx pxy exp−λy h [hI{λx =λy } + ( )I{λx 6=λy } ] + o(h).
λy − λx

Por lo anterior
lı́m h−1 pxy (h) = qxy .
h→0+
Caso 2. x = y.
Por el Caso 1 se tiene

lı́m h−1 (pxx (h) − 1) = −λx = qxx .

h→0+

Definición 2.30. El generador infinitesimal de un PSM, es la matriz Q = qxy x,y∈E con entradas.
qxy = λx pxy si x 6= y y qxx = −λx .

Nota 2.10. Observemos que el generador infinitesimal contiene toda la información del proceso, tanto
las intensidades de los tiempos de salto λx como la distribución de los saltos exp(−λx t).

El siguiente Teorema nos muestra que el generador infinitesimal de un PSM es suficiente para construir
las probabilidades de transición {pxy (t)}x,y∈E .

Teorema 2.9. Dado un PSM {Xt }t≥0 con espacio de estados E y Q su correspondiente generador infi-
nitesimal, entonces se satisface el siguiente sistema de ecuaciones diferenciales.

∂P (t)
= QP (t).
∂t

Demostración. Sea ε > 0

Por la igualdad de Chapman-Kolmogorov(2.16) tenemos

pxy (t + h) − pxy (t) X (pxz (h) − δxz )pzy (t)

= .
h h
z∈E

Ahora basta mostrar que

pxz (h) − δxz
lı́m = qxz .
h→0+ h
Pero esto es resultado de la Proposición 2.17.
2.3. PROCESOS DE SALTOS DE MARKOV 23

Nota 2.11. El teorema 2.9 muestra que P (t) cumple la ecuación diferencial
∂P (t)
= QP (t),
∂t
P∞ tn Qn
si la matriz Q es diagonalizable, entonces P (t) = exp(tQ) = n=0 n! , es la única solución al sistema.

Ahora bien, como la tasa de saltos λx es independiente del tiempo t ≥ 0 y sólo depende del estado x ∈ E,
clasificaremos a los estados.

Definición 2.31. Sea x ∈ E con tasa de salto asociada λx , entonces x se define como:

1. Permanente, si λx = 0.
2. Estable, si 0 < λx < ∞.
3. Instantáneo, si λx = ∞.

Sea {Xt }t≥0 un PSM con Xt = x a tiempo t ≥ 0, estudiaremos las consecuencias del proceso dependiendo
cada tipo de estado.

1. Si x es permanente, P (∆(t) = ∞ | Xt = x) = 1 para toda s > t. Es decir, el PSM permanece

siempre en el estado i ∈ E.
2. Si x es estable, P (0 < ∆(t) < ∞ | Xt = x) = 1 , es decir el proceso saldrá del estado i en un tiempo
de espera positivo y finito.
3. Si x es instantáneo, P (∆(t) = 0 | Xt = x) = 1 , entonces el PSM existe en el estado tan pronto
entra en él.

Proposición 2.18. (Regularidad). Considere un PSM {Xt }t≥0 y T1 < . . . < Tn sus respectivos n
tiempos de salto, asumamos que su correspondiente generador infinitesimal tiene entradas estrictamente
negativas y finitas en la diagonal, bajo estas condiciones se tiene que

1. T1 < . . . < Tn < ∞ casi seguramente.

c.s.
2. Tn −−−→ ∞.
n→∞

Demostración. Por el Lema 2.1 podemos argumentar que para toda t ≥ 0

n−1
X
P[Tn ≤ t] ≤ 1 − exp(−βt) (βt)j /j!.
j=0

Donde β=máx{λx : x ∈ E}, ası́ concluimos que

T1 < . . . < Tn < ∞ casi seguramente.
La segunda es resultado de hacer tender n a infinito y notar que P[ lı́m Tn ≤ t] ≤ lı́m P[Tn ≤ t] = 0 para
n→∞ n→∞
cada t ≥ 0.
24 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Nota 2.12. A partir de ahora solo consideraremos procesos cuyos generadores infinitesimales tengan
entradas en la diagonal estrictamente negativas. La existencia de estos procesos se probará más adelante,
especı́ficamente en el Teorema 3.2.

2.3.4. Distribuciones estacionarias y distribuciones lı́mite

Definición 2.32. Sea un PSM {Xt }t≥0 con espacio de estados E.

Se dice que {Xt }t≥0 es irreducible si su cadena de Markov asociada a los saltos {Zn }n≥0 (ver Definición
2.28) es irreducible.
Por otro lado decimos que el proceso PSM es recurrente si {Zn }n≥0 es recurrente.

Definición 2.33. Dado un PSM {Xt }t≥0 se define el tiempo de llegada hacia el estado x como:

τx = ı́nf{t > T1 : Xt = x}.

Un PSM se dice ser recurrente positivo si para todo x ∈ E se tiene que E[τx |X0 = x] < ∞.

Definición 2.34. Un PSM {Xt }t≥0 se dice ergódico si es irreducible y recurrente positivo.

En lo que resta de este capı́tulo enunciaremos resultados que son análosgos al caso de cadenas de Markov
en tiempo discreto y sus demostraciones serán omitidas y podran ser consultadas en [? ].

Proposición 2.19. Considere un PSM {Xt }t≥0 que es recurrente e irreducible, entonces {Xt }t≥0 es
ergódico.

Definición 2.35. Considere un PSM {Xt }t≥0 , una medida de probabilidad ν = {νy }y∈E sobre E es una
distribución estacionaria si para todo y ∈ E y t ≥ 0
X
νy = νx pxy (t).
x∈E

Teorema 2.10. Considere un PSM {Xt }t≥0 ergódico, entonces las siguientes afirmaciones son válidas.

1. Existe una única distribución estacionaria ν dada por

Rτ
E[ 0 x I{Xt =y} dt|X0 = x]
νy = , para todo y ∈ E.
E[τx |X0 = x]
La cual es independiente del estado x.

2. La distribución ν satisface la ecuación

X
νx qxy = 0 , para todo y ∈ E.
x∈E

Demostración. Para su demostración ver ([? ] página 261).

Teorema 2.11. (Ergódico). Dado un PSM ergódico {Xt }t≥0 con distribución estacionaria ν y una
función f : E → R, Z t X
−1
lı́m t f (Xs )ds = f (x)νx casi seguramente.
t→∞ 0 x∈E
2.3. PROCESOS DE SALTOS DE MARKOV 25

Demostración. Para su demostración ver ([? ] página 265).

Definición 2.36. Dado un PSM {Xt }t≥0 se define la distribución lı́mite como

lı́m pxy (t) para todo y ∈ E.

t→∞

Siempre que está exista y no dependa de x.

Teorema 2.12. (Distribución lı́mite) Dado un PSM ergódico {Xt }t≥0 , existe la distribución lı́mite y
está dada por

lı́m pxy (t) = νy para todo y ∈ E.

t→∞

Donde ν = {νy }y∈E está definida en el Teorema 2.10

Notemos que este lı́mite no depende de x.
Teorema 2.13. Dado un PSM ergódico {Xt }t≥0 y A, B ⊆ E se tiene
Nt
X X X
−1
lı́m t I{Zn ∈A,Zn+1 ∈B} = νx qxy .
t→∞
n=0 x∈A y∈B

Este lı́mite puede entenderse como el promedio de transiciónes de B desde A por unidad de tiempo.

Demostración. La demostración se encuentra en ([? ] página 269).

2.3.5. Ejemplos

Ejemplo 2.11. (Proceso de nacimiento y muerte) Un proceso de nacimiento y muerte es un PSM

con espacio de estados E = {0, 1, . . .} o E = {0, 1, . . . , k}, si E = {0, 1, . . .} el generador infinitesimal es
 
−λ0 λ0 0 0 ...
 α1 −(λ1 + α1 ) λ1 0 ... 
Q= 0 .
 
 α2 −(λ2 + α2 ) λ2 . . .  
.. .. .. .. ..
. . . . .

Mientras que si E = {0, 1, . . . , k} el generador infinitesimal es

 
−λ0 λ0 0 0 ... 0
 α1 −(λ1 + α1 ) λ 1 0 ... 0 
 
Q=
 0 α 2 −(λ 2 + α 2 ) λ 2 ... 0 .

 .. .. .. .. .. .. 
 . . . . . . 
0 0 0 0 αk −αk

Donde λ0 , λ1 , . . . y α1 , α2 , . . . son constantes positivas conocidas como las tasas instantáneas de nacimiento
y muerte, respectivamente. De acuerdo a lo desarrollado antes, a partir de esta matriz podemos concluir
que el tiempo de estancia en cualquier estado x ∈ E tiene distribución exp(λx + αx ), en donde se define
α0 = 0. Las probabilidades de saltos de un estado a otro son

λx
 λx +αx si y = x + 1,

αx
pxy = λx +αx si y = x − 1, (2.1)

0 otro caso.

26 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

De este modo, un proceso de nacimiento y muerte permanece en cada uno de sus estados un tiempo
exponencial, al cabo del cual salta una unidad hacia arriba o una unidad hacia abajo de acuerdo a las
probabilidades de (2.1). Un salto hacia arriba se interpreta como un nacimiento, mientras que un salto
hacia abajo representa una muerte. La variable Xt puede interpretarse como el número de individuos en
una población al tiempo t, en donde pueden presentarse nacimientos y muertes de individuos, uno a la
vez. Como λ0 > 0 y α0 = 0, la población puede crecer cuando se encuentre en el estado cero, pero nunca
decrecer por debajo de ese nivel.
Nos concentraremos en el caso donde E = {0, 1, . . . , k} y las tasas sean estrictamente positivas, podemos
notar que para este caso λk = 0, esto quiere decir que la población no puede ser mayor que k.
Por el Teorema 2.9 sabemos que P (t) cumple el siguiente sistema de ecuaciones diferenciales.

∂P (t)
= QP (t).
∂t
Donde Q es el generador infinitesimal y P (t) es la matriz de transición al tiempo t. Para encontrar
P (t) debemos probar primero que Q es diagonalizable, y por la teorı́a de ecuaciones diferenciales podemos
concluir que
P (t) = exp(Qt).
Sea D ∈ Mk+1×k+1 una matriz diagonal con entradas positivas en la diagonal
 
d0 0 0 0 . . . 0
 0 d1 0 0 . . . 0 
 
D =  0 0 d2 0 . . . 0  ,
 
 .. .. .. .. .. .. 
 . . . . . . 
0 0 0 0 0 dk

la matriz resultante del producto D−1 QD esta dada por

λ0 dd01
 
−λ0 0 0 ... 0
 α d0 −(λ + α ) d2
λ1 d1 0 ... 0 
 1 d1 1 1 
d1
−1 0 α2 d2 −(λ2 + α2 ) λ2 dd23 ... 0
 
D QD =   .
.
 . .. .. .. .. .. 
 . . . . . .


dk−1
0 0 0 0 αk dk −αk

Nuestro objetivo es encontrar una matriz D diagonal tal que D−1 QD sea una matriz simétrica, de este
modo se forma el siguiente sistema de ecuaciones
d0 d1
α1 = λ0
d1 d0
d1 d2
α2 = λ1
d2 d1
..
.
dk−1 dk
αk = λk−1 .
dk dk−1

Por lo anterior obtenemos r

αn
dn = dn−1 para toda 1 ≤ n ≤ k,
λn−1
haciendo d0 = 1 obtenemos s
Mn
dn = para toda 1 ≤ n ≤ k,
Ln−1
2.3. PROCESOS DE SALTOS DE MARKOV 27

donde Mn = ni=1 αi y Ln−1 = n−1

Q Q
i=0 λi .
−1
Sea U = D QD, por lo anterior
 √ 
√−λ 0 λ 0 α1 √0 0 ... 0
 λ0 α1 −(λ1 + α1 ) λ1 α2
√ √0 ... 0 
 
U =
 0 λ α
1 2 −(λ 2 + α 2 ) λ α
2 3 . . . 0 ,

 .. .. .. .. .. .. 
 . . . . p . . 
0 0 0 0 λk−1 αk −αk

por el Teorema espectral ([? ] página 481) U puede descomponerse en U = V −1 W V , con W una matriz
diagonal con valores reales.
En conclusión Q es diagonalizable y se puede descomponer en

Q = DV −1 W V D−1 = D̃−1 W D̃,

con D̃ = V D−1 . Por lo tanto

P (t) = D̃−1 exp(W t)D̃.
El proceso es irreducible y recurrente, pues las probabilidades de (2.1) producen una cadena de Markov
irreducible y recurrente. Por la Proposición 2.19 sabemos que también es ergódico, y por el Teorema 2.10
hay una única distribución estacionaria que cumple
X
νx qxy = 0 , para todo y ∈ E,
x∈E

que se traduce en

−λ0 ν0 + α1 ν1 = 0
λ0 ν0 − (λ1 + α1 )ν1 + α2 ν2 = 0
λ1 ν1 − (λ2 + α2 )ν2 + α3 ν3 = 0
..
.
λk−2 νk−2 − (λk−1 + αk−1 )νk−1 + αk νk = 0
−αk νk + λk−1 νk−1 = 0

acomodando los términos obtenemos

α1 ν1 − λ0 ν0 = 0
α2 ν2 − λ1 ν1 = α1 ν1 − λ0 ν0
α3 ν3 − λ2 ν2 = α2 ν2 − λ1 ν1
..
.
αk νk − λk−1 νk−1 = αk−1 νk−1 − λk−2 νk−2
αk νk − λk−1 νk−1 = 0

por lo anterior αn νn − λn−1 νn−1 = 0 para toda 1 ≤ n ≤ k, es decir

λn−1 λ0 . . . λn−1
νn = νn−1 = . . . = ν0 .
αn α1 . . . αn
Sumando todo el espacio de estados E tenemos
1
ν0 = Pk λ0 ...λi−1
.
1+ i=1 α1 ...αi
28 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS

Por el Teorema 2.10 sabemos que

1
νy = , para todo y ∈ E,
λy E[τy |X0 = y]
1
por tanto E[τy |X0 = y] = νy λy , para todo y ∈ E.

Ejemplo 2.12. (Proceso de Poisson) El proceso de Poisson homogéneo es un PSM que empieza en 0
es decir, la distribución inicial tiene el valor 1 en el estado 0. Los tiempos de estancia son exponenciales
de parámetro λ = qx y las probabilidades de salto de un estado a otro son
(
1 si y = x + 1,
pxy =
0 si y 6= x + 1.

Ası́ con nuestro generador infinitesimal podremos ir evolucionando el proceso de Poisson hasta un tiempo
T , el cual es definido como: (
λ si y = x + 1,
qxy =
0 e.o.c.

Considerando Q ∈ MN ×N : (
qxy si y 6= x,
Q(x, y) =
−λ si y = x.

Ası́ tenemos:  
−λ λ 0 0 0 ...
 0 −λ λ 0 0 ...
 
Q= 0
 0 −λ λ 0 ...

 0
 0 0 −λ λ ...

.. .. .. .. . .
. . . . . ...
Como el proceso siempre va creciendo, el proceso se degenera y no tiene solución. Ası́ sugeriremos condi-
ciones de frontera con espacio de estados E = {0, 1, 2, 3, ..., N }, tomando el sistema de N + 1 ecuaciones,
ası́ obtendrı́amos el siguiente generador infinitesimal:
 
−λ λ 0 0 0 ··· 0 0 0
 0 −λ λ 0 0 · · · 0 0 0 
 
 0 0 −λ λ 0 · · · 0 0 0 
Q= .
 
. .. .. .. . . .. .. .. 
 .
 . . . . ··· . . .
 0 0 0 0 0 · · · 0 −λ λ
0 0 0 0 0 ··· 0 0 0
Del cual obtenemos el siguiente sistema de ecuaciones:


 −λν0 = 0

λν0 − λν1 = 0





λν1 − λν2 = 0

..


 .

λνN −1 = 0





ν1 + ν2 + ν3 + ... + νN = 1

2.3. PROCESOS DE SALTOS DE MARKOV 29

Ası́ νN = 1 y la distribución estacionaria del proceso de poisson con condiciones de frontera serı́a 1 en la
última entrada: ν = (0, 0, 0, · · · , 1)
30 CAPÍTULO 2. PROCESOS ESTOCÁSTICOS
Capı́tulo 3

Simulación de procesos de Markov

En este capı́tulo se pretende mostrar las principales técnicas de simulación de algunos procesos estocásti-
cos, en particular nos concentraremos en los procesos estocásticos de Markov cuyas aplicaciones van desde
las finanzas computacionales e investigación de operaciones hasta la genética y genómica humana. Este
capı́tulo se centra en dos tipos de procesos de Markov, primero estudiaremos las técnicas de simulación
sobre cadenas de Markov en tiempo discreto ası́ como su respectiva justificación matemática, para poste-
riormente resolver algunos problemas analı́ticos mediante la simulación de cadenas de Markov. Posterior-
mente trabajaremos con procesos de saltos de Markov, mostraremos algunos resultados matemáticos que
nos permitirán justificar la simulación de dichos procesos y finalmente resolveremos problemas analı́ticos
mediante la simulación de dichos procesos. La bibligrafı́a básica para este capı́tulo es [? ].

3.1. Cadenas de Markov

3.1.1. Algoritmo para cadenas de Markov

Sabemos que dada una cadena de Markov homogénea podemos obtener su matriz de transición y su
distribución inicial, ahora nos preguntamos si dada una matriz estocástica y una distribución inicial, es
posible generar una cadena de Markov. El siguiente resultado nos asegura que sı́ se puede y además
proporciona herramientas para generar trayectorias de cadenas de Markov.

Teorema 3.1. Dada una matriz estocástica P = {pij }1≤i,j≤N ∈ MN ×N , una distribución inicial π sobre
el espacio de estados E = {1, 2, . . . , N }, existe un espacio de probabilidad en el que están definidas una
sucesión de variables aleatorias {Xn }n∈N en E, que conforman una cadena de Markov homogénea con
espacio de estados E, matriz de transición P y distribución inicial π.

Demostración. Consideremos un espacio de probabilidad (Ω, F , P) en el cual está definida una sucesión
{Un }n∈N de variables aleatorias independientes uniformes en el intervalo (0,1). Definimos φi (y) como la
función de cuantiles (ver A.1) sobre E asociada a la distribución {pij }1≤j≤N y a φ0 (y) como la función
de cuantiles sobre E asociada a la distribución π, con y ∈ (0, 1).
Sea X0 = φ0 (U0 ) y Xn+1 = φXn (Un+1 ) para n ≥ 0.

Analicemos las probabilidades de transición

31
32 CAPÍTULO 3. SIMULACIÓN DE PROCESOS DE MARKOV

P[Xn = xn |Xn−1 = xn−1 ] = P[φXn−1 (Un ) = xn |Xn−1 = xn−1 ]

Como Xn−1 sólo depende de {Uk }k≤n−1 y además las variables aleatorias {Uk }k∈N son independientes se
tiene que
P[φXn−1 (Un ) = xn |Xn−1 = xn−1 ] = P[φxn−1 (Un ) = xn ].
Por último P[φxn−1 (Un ) = xn ] = pxn−1 xn por definición de función de cuantiles.

Mostraremos que la sucesión de variables aleatorias {Xn }n∈N , cumple la propiedad de Markov.

P[X0 = x0 , X1 = x1 , . . . , Xn = xn ] = P[φ0 (U0 ) = x0 , φX0 (U1 ) = x1 , . . . , φXn−1 (Un ) = xn ]

= P[φ0 (U0 ) = x0 , φx0 (U1 ) = x1 , . . . , φxn−1 (Un ) = xn ]
por la independencia de {Uk }k∈N
= P[φ0 (U0 ) = x0 ]P[φx0 (U1 ) = x1 ] . . . P[φxn−1 (Un ) = xn ]
= P[X0 = x0 ]P[X1 = x1 |X0 = x0 ] . . . P[Xn = xn |Xn−1 = xn−1 ]
n−1
Y
= P[X0 = x0 ] pxi xi+1 .
i=0

Por lo tanto {Xn }n∈N es una cadena de Markov, con distribución inicial π y matriz de transición P .

El teorema anterior nos proporcionó la forma de generar trayectorias de una cadena de Markov homogénea
sobre E con distribución inicial y matriz de transición dadas. A continuación daremos una construcción
alternativa que nos servirá más adelante, sin embargo no se probará que esta construcción produce una
cadena de Markov ya que es similar a la propuesta en el Teorema 3.1

Definición 3.1. (Construcción 2 Cadenas de Markov). Supongamos las mismas hipótesis que en
el Teorema 3.1. Consideremos un espacio de probabilidad donde están definidas {wi,k }i∈E,k∈N variables
aleatorias independientes tales que

P[wi,k = j] = pij para todo i, j ∈ E y k ∈ N. (3.1)

Definimos a X0 = φ0 (U ) y Xk+1 = wYk ,nY (k)+1 para cada k > 0.

k
Donde U es uniforme en (0,1) independiente de {wi,k }i∈E,k∈N , φ0 la función de cuantiles correspondiente
a la distribución inicial π = (π1 , π2 , . . . , πN ) y ni (k) es el número de visitas al estado i hasta el tiempo
k − 1, definiendo ni (0) = 0 para todo i ∈ E y k ∈ N.
Análogamente como se procedió en el Teorema 3.1 se demuestra que {Xn }n∈N es una cadena de Markov
sobre E con matriz de transición P y distribución inicial π.

El resultado anterior nos proporciona una justificación de la existencia de cadenas de Markov, ası́ como
también nos señala una técnica para su simulación misma que se describe a continuación.
3.1. CADENAS DE MARKOV 33

Algoritmo 3.1. Cadenas de Markov a tiempo discreto. Este algoritmo genera trayectorias de una
cadena de Markov homogénea con distribución inicial π y matriz de transición P hasta el horizonte de
tiempo m.

1. Generar X0 de la distribución inicial π.

2. Inicializamos n=1.

3. Generamos Xn de la distribución correspondiente al renglón de la matriz de transición asociado al

estado Xn−1 .

4. Si n es igual al horizonte de tiempo m, entonces detenemos el proceso, si no es ası́ actualizamos

n = n + 1 y volvemos al paso 3.
Ejemplo 3.1. Se simuló una trayectoria de una cadena de Markov usando el Algoritmo 3.1 con la función
del Código D.1 en el software Julia con las siguientes caracterı́sticas:
Matriz de transición

 
0.2 0.1 0.5 0.2
 0.5 0.2 0.2 0.1 
P =
 0.2

0.5 0.1 0.2 
0.3 0.5 0.1 0.1
vector de distribución inicial π:

0.2 0.4 0.2 0.2
vector de espacio de estados E:

1 2 3 4 .
Con horizonte de tiempo m=40.
En la figura ?? podemos observar la gráfica de una trayectoria simulada.

3.1.2. Cálculo numérico de propiedades de cadenas de Markov

A continuación aplicaremos técnicas de simulación desarrolladas anteriormente para realizar cálculos

numéricos sobre propiedades de cadenas de Markov.
Ejemplo 3.2. Considere la cadena de Markov del Ejemplo 2.10 con k=3, por lo anterior el espacio de
estados es E = {0, 1, 2, 3} y la matriz de transición es
 
0 1 0 0
 1/3 0 2/3 0 
P =  0 2/3 0 1/3  .


0 0 1 0

En el Ejemplo 2.10 mostramos analı́ticamente que la cadena es irreducible. La función del Código D.7
determina cuales son las clases de comunicación de una cadena de Markov, como resultado nos muestra
que para la cadena de Ehrenfest existe una única clase de comunicación

{0, 1, 2, 3},

con esto podemos concluir numéricamente que esta cadena es irreducible.

Con la función del Código D.8 es posible calcular el periodo de una cadena de Markov irreducible, como
34 CAPÍTULO 3. SIMULACIÓN DE PROCESOS DE MARKOV

resultado nos muestra que la cadena de Ehrenfest tiene periodo 2, este resultado coincide con las conclu-
siones del Ejemplo 2.10.
Siguiendo con nuestro análisis la función del Código D.9 nos muestra cuales son los estados recurrentes
de una cadena de Markov dada, en el caso de la cadena de Ehrenfest nos muestra que todos los estados
son recurrentes.
En el Ejemplo 2.10 calculamos analı́ticamente la única distribución estacionaria de la cadena de Ehrenfest,
para el caso k = 3 obtenemos
ν = (1/8, 3/8, 3/8, 1/8),
mientras que con la función del Código D.10 con un horizonte de tiempo m = 10000 produce una distri-
bución estacionaria aproximada de
(0.1247, 0.3787, 0.3753, 0.1213).

Propiedad Resultado Código utilizado

Clases de comunicación {0,1,2,3} D.7
Irreducible si D.7
Periodo 2 D.8
Estados recurrentes {0,1,2,3} D.9
Estados transitorios ∅ D.9
Ergódica si D.7
Distribución estacionaria (0.1247, 0.3787, 0.3753, 0.1213) D.10

Cuadro 3.1: Tabla resumen de resultados .

En la figura ?? podemos observar el comportamiento de la convergencia a la distribución estacionaria

hasta el horizonte de tiempo m = 100000, para este análisis se implementó la función del Código D.11

3.2. Procesos de saltos de Markov

En esta sección daremos una justificación formal para la simulación de procesos de Markov en tiempo
continuo con espacio de estados finito E = {1, . . . , N }.
Anteriormente estudiamos los procesos de saltos de Markov, ahora nos interesa simular trayectorias de
dichos procesos, para ello debemos primero mostrar su existencia, debido a la Nota 2.9 podemos plan-
tearnos como hacer la construcción de un PSM.

3.2.1. Algoritmo para PSM

Definición 3.2. Definamos S0 = 0 y {Sn }n≥1 variables aleatorias i.i.d. con distribución exponencial de
parámetro 1 y {λx }x∈E números estrictamente positivos.
Tomemos también {Zn }n∈N una cadena de Markov con espacio de estados E independiente de la sucesión
de variables {Sn }n≥0 , definimos el proceso en tiempo continuo dado por :
 
∞ n n+1
X X Sj X Sj 
Xt = Zn I ≤t< . (3.2)
 λZj−1 λZj−1 
n=0 j=0 j=0
3.2. PROCESOS DE SALTOS DE MARKOV 35

Nos preguntamos si este nuevo proceso es un proceso PSM, la respuesta es afirmativa y es parte del
siguiente Teorema.
Teorema 3.2. (Existencia de PSM).
Sea Q = (qxy )x,y∈E una matriz con entradas qxy = λx pxy si x 6= y y qxx = −λx , en donde {λx }x∈E son
números mayores que cero y para todo x ∈ E {pxy }y∈E conforma una distribución sobre E tal que pxx = 0,
permitir además ser a π ∈ RN una distribucion inicial en E.
Bajo estas condiciones existe un espacio de probabilidad, donde está definido un PSM {Xt }t≥0 con espacio
de estados E, con distribución inicial π y cuyo generador infinitesimal está dado por Q.

Demostración. Definamos el proceso estocástico {Xt }t≥0 como en la Definición 3.2 en donde la cadena
de Markov {Zn }n∈N tiene matriz de transición P = {pxy }x,y∈E y distribucion inicial π.
Para demostrar la afirmación basta con demostrar que para x0 , x1 , . . . , xn+1 ∈ E y para cualesquiera
ti < ti+1 = si+1 + ti con i ∈ {0, . . . , n}, con t0 = 0

P[Xtn+1 = xn+1 |Xtn = xn , . . . , Xt0 = x0 ] = P[Xtn+1 = xn+1 |Xtn = xn ].

P∞
Definamos Nt = n=1 I{Tn ≤t} el número de saltos que el proceso ha realizado hasta el instante t, con
Tn = nj=0 Sj /λZj−1 si n ≥ 0, entonces
P

∞
X
P[Xtn+1 = xn+1 |Xtn = xn , . . . , Xt0 = x0 ] = [am + bm ]P[Ntn = m|Xtn = xn , . . . , Xt0 = x0 ]. (3.3)
m=0

En donde

am = P[Xtn +sn+1 = xn+1 , Tm+1 > tn + sn+1 |Fm ]

bm = P[Xtn +sn+1 = xn+1 , Tm+1 ≤ tn + sn+1 |Fm ]
Fm = {Ntn = m, Xtn = xn , . . . , Xt0 = x0 }.

Debido a que Sm+1 es independiente a {(Zk , Tk )}k≤m y por la propiedad de pérdida de memoria de la
exponencial se tienen las siguientes igualdades.

am = I{xn+1 =xn } P[Tm+1 > tn + sn+1 |Ntn = m, Xtn = xn , . . . , Xt0 = x0 ]

= I{xn+1 =xn } P[Tm+1 > tn + sn+1 |Tm+1 > tn ≥ Tm , Zm = xn , Xtn−1 = xn−1 , . . . , Xt0 = x0 ]
= I{xn+1 =xn } P[Sm+1 /λxn > (tn − Tm ) + sn+1 |Sm+1 /λxn > (tn − Tm ), (tn − Tm ) ≥ 0]
= I{xn+1 =xn } exp(−λxn sn+1 ).

Por otro lado condicionando respecto a Tm+1 y respecto a XTm+1 = Zm+1 y haciendo un par de desarrollos
se tiene:
X Z sn+1
bm = fXtn +sn+1 |Tm+1 ,Zm+1 ,Ntn ,Xtn ,...,Xt0 (xn+1 |tn + v, k, m, xn , . . . , x0 )
k6=xn 0

×fZm+1 ,Tm+1 |Ntn ,Xtn ,...,Xt0 (k, tn + v|m, xn , . . . , x0 )dv.

Por construcción del proceso {Xt }t≥0 , se sabe que

fXtn +sn+1 |Tm+1 ,Zm+1 ,Ntn ,Xtn ,...,Xt0 (xn+1 |tn +v, k, m, xn , . . . , x0 ) = P[Xsn+1 −v = xn+1 |X0 = k] = pkxn (sn+1 −v).

Por un argumento análogo al que se utilizó para calcular am y por la independencia que existe entre la
cadena de Markov {Zn }n∈N y la sucesión {Sn }n∈N obtenemos

fZm+1 ,Tm+1 |Ntn ,Xtn ,...,Xt0 (k, tn + v|m, xn , . . . , x0 ) = λxn exp(−λxn v)pxn k .
36 CAPÍTULO 3. SIMULACIÓN DE PROCESOS DE MARKOV

Uniendo estas igualdades podemos al fin calcular fácilmente bm .

X Z sn+1
bm = pkxn (sn+1 − v) ∗ λxn exp(−λxn v)pxn k dv.
k6=xn 0

Uniendo los cálculos anteriormente hechos y notando que am y bm son independientes de m y de tn se

tiene que (3.3) es igual a:
X Z sn+1
I{xn+1 =xn } exp(−λxn sn+1 ) + pkxn (sn+1 − v) ∗ λxn exp(−λxn v)pxn k dv.
k6=xn 0

Cuya cantidad no depende del tiempo anterior a tn , ni tampoco depende de tn y haciendo un cálculo
sencillo se tiene además que
P[Xsn+1 = xn+1 |X0 = xn ]
es igual a
X Z sn+1
I{xn+1 =xn } exp(−λxn sn+1 ) + pkxn (sn+1 − v) ∗ λxn exp(−λxn v)pxn k dv.
k6=xn 0

Por lo que concluimos

P[Xtn+1 = xn+1 |Xtn = xn , . . . , Xt0 = x0 ] = P[Xsn+1 = xn+1 |X0 = xn ].

Por lo anterior afirmamos que es un proceso de Markov continuo por la derecha y por los Teoremas 2.6,
2.7 y 2.8 se sigue que el proceso anteriormente construido posee la distribución π y generador infinitesimal
Q.

Definición 3.3. (Construcción 2 de un PSM). Supongamos las mismas hipótesis que en el Teorema
3.2. Sean {Si,0 = 0}i∈E , {Si,n }i∈E,n≥1 variables aleatorias i.i.d. con distribución exponencial de parámetro
1 y {λx }x∈E números estrictamente positivos, consideremos también {Zn }n∈N una cadena de Markov con
espacio de estados E, matriz de transición P y distribución inicial π, dicha cadena como en la Definición
3.1, e independiente de la sucesión de variables {Si,n }i∈E,n≥1 , definimos el proceso en tiempo continuo
dado por
 
∞ n S n+1
X SZj−1 ,nZj−1 (j)   
X X Zj−1 ,nZj−1 (j) SZ0 ,1
Xt = Zn I ≤ t < + Z0 I 0 ≤ t < . (3.4)
 λZj−1 λZj−1  λZ0
n=1 j=1 j=1

Donde ni (m) es la variable asociada al número de transiciónes de i a cualquier estado hasta el tiempo m
asociada a la cadena {Zn }n∈N .

Análogamente como en el Teorema 3.2 se demuestra que {Xt }t≥0 es un PSM con los mismos parámetros
dados en este teorema.
El Teorema 3.2 da lugar al siguiente algoritmo para simular un PSM.
Algoritmo 3.2. PSM. Este algoritmo genera trayectorias de un PSM con distribución inicial π y gene-
rador infinitesimal Q, hasta el horizonte de tiempo T.

1. Generar Z0 de la distribución inicial π y hacer i = 0, T0 = 0.

2. Dado que se está en el estado Zi , simular Ti+1 v.a. exponencial de parámetro λZi .
3.2. PROCESOS DE SALTOS DE MARKOV 37

3. Hacer Xt = Zi si Ti ≤ t < Ti+1 .

4. Simular una v.a. Zi+1 de la distribución de salto {fXT1 |X0 (y|Zi )}y∈E independientemente del tiempo
en el que ocurrió.

5. Hacer i = i + 1.

6. Si Ti < T volver al paso 2, de lo contrario ir al paso 7.

7. Tomar la trayectoria del proceso solamente en el intervalo [0,T], en donde T es el horizonte de

tiempo.

Ejemplo 3.3. Se simuló una trayectoria de un PSM usando el Algoritmo 3.2 con la función del Código
D.4 en el software Julia con las siguientes caracterı́sticas:
Con generador infinitesimal:  
−.13 .04 .01 .08
 .02 −.10 .04 .04 
Q= 
 .2 .1 −.6 .3 
.05 .15 .1 −.3
vector de espacio de estados E:

1 2 3 4
vector de distribución inicial π:


.2 .2 .1 .5 .

Con horizonte de tiempo T=50.

En la figura ?? podemos observar la gráfica de una trayectoria simulada.

3.2.2. Cálculo numérico de propiedades de PSM

A continuación aplicaremos técnicas de simulación desarrolladas anteriormente para realizar cálculos

numéricos sobre propiedades de los procesos de saltos de Markov

Ejemplo 3.4. Considere la cadena de Markov del Ejemplo 2.11 con k=3 y λ0 = 1, λ1 = 2, λ2 = 1,
µ1 = 2, µ2 = 1, µ3 = 2 por lo anterior el espacio de estados es E = {0, 1, 2, 3} y el generador infinitesimal
es  
−1 1 0 0
 2 −4 2 0 
Q= ,
 0 1 −2 1 
0 0 2 −2
mientras que la matriz de transición asociada a la cadena de saltos es
 
0 1 0 0
 0.5 0 0.5 0 
P =  0 0.5 0 0.5  .


0 0 1 0
38 CAPÍTULO 3. SIMULACIÓN DE PROCESOS DE MARKOV

Por otro lado implementando la función del Código D.12

P (t) = D̃−1 exp(W t)D̃,

donde  
0.210527 0.2 0.5 −0.640487
 −0.891808 −0.4 0.5 −0.151198 
D̃−1 = ,
 0.34064 −0.4 0.5 0.395843 
−0.210527 0.8 0.5 0.640487

 
exp(−5.23607t) 0.0 0.0 0.0
 0.0 exp(−3t) 0.0 0.0 
exp(W t) =  ,
 0.0 0.0 exp((1.53006e − 19)t) 0.0 
0.0 0.0 0.0 exp(−0.763932t)

 
−0.362866 −0.768563 0.58713 −0.181433
 0.333333 −0.333333 −0.666667 0.666667 
D̃ = 
 0.666667
.
0.333333 0.666667 0.333333 
−0.817513 −0.0964944 0.505251 0.408757
La distribución estacionaria esta dada por

(1/3, 1/6, 1/3, 1/6),

mientras que la aproximación implementando la función del Código D.15 con horizonte de tiempo T =
10000 es
(0.332264, 0.16706, 0.332738, 0.16797).
3.2. PROCESOS DE SALTOS DE MARKOV 39

Propiedad Resultado Código utilizado

Clases de comunicación {0,1,2,3} D.13
Irreducible si D.13
Estados recurrentes {0,1,2,3} D.14
Estados transitorios ∅ D.14
Ergódica si D.13
Distribución estacionaria (0.332264, 0.16706, 0.332738, 0.16797) D.15

Cuadro 3.2: Tabla resumen de resultados .

En la figura ?? podemos observar el comportamiento de la convergencia a la distribución estacionaria

hasta el horizonte de tiempo T = 100000, para este análisis se implementó la función del código D.16.

3.2.3. Ejemplo. Modelo SIR-PSM

El modelo SIR (Susceptibles-Infectados-Recuperados) es un modelo de compartimentos donde la población

bajo estudio se divide en clases epidemiológicas y se describe un flujo entre ellas.

Para poder modelar estocásticamente el fenómeno consideremos las variables aleatorias discretas:

S(t) = número de suceptibles a tiempo t

I(t) = número de infectados

Como el modelo no considera muertes consideraremos una población con N habitantes, es decir S(t), I(t) ∈
{0, 1, 2, ..., N } con t ∈ [0, ∞).

Para entender el modelo como uno de PSM, necesitamos encontrar las probabilidades de transición.
Denotamos a las as probabilidades de transición asociadas al proceso estocástico como:

p(s,i),(s+k,i+j) (M t) = P (S(t+ M t), I(t+ M t) = (s + k, i + j) | S(t), I(t) = (s, i))

Las cuales dependen del tiempo entre los eventos, es decir de M t y no especı́ficamente del tiempo t ası́
tenemos un proceso de tiempo homogéneo. Además dado el estado actual al tiempo t > 0, el estado
futuro al tiempo t+ M t no depende de tiempos anteriores a t, por lo que cumple la propiedad de Markov;
entonces consideraremos las probabilidades de transición definidas de la siguiente manera:

 S

 βi N M t + o(M t) si (k, j) = (−1, +1),




 γi M t + o(M t) si (k, j) = (0, −1),


p(s,i),(s+k,i+j) (M t) =
S


1 − βi N + γi M t + o(M t) si (k, j) = (0, 0),









o(M t) e.o.c.
Donde β es la tasa de transmisión y γ la tasa de recuperación.

En resumen mostraremo los cambios asociados con los eventos de infección y recuperación en la siguiente
tabla.
40 CAPÍTULO 3. SIMULACIÓN DE PROCESOS DE MARKOV

Supuestos del modelo SIR

Evento Cambio (M S, M T ) Probabilidad
S
Infección (−1, +1) βi N M t + o(M t)
Recuperación (0, −1) γi M t + o(M t)

Figura 3.1: Gráfica de la norma entre P y P̂ (n) respecto al tiempo

Tamaño de la población N=200 Tiempo de observación t=50 tasa de infección beta=0.30 tasa de recu-
peración gamma=0.15 condiciones iniciales S0 = N-2 I0 = 2 R0 = 0
Capı́tulo 4

Inferencia sobre procesos de Markov

A continuación estudiaremos la estimación de los parámetros de trayectorias de cadenas de Markov, ası́

como también demostraremos algunas caracterı́sticas estadı́sticas importantes sobre los estimadores e
implementaremos ejemplos numéricos donde se ilustren los principales resultados. La bibliografı́a básica
para este capı́tulo es principalmente [? ].

4.1. Cadenas de Markov

4.1.1. Estimación máximo verosı́mil

Supóngase que se tiene observada una trayectoria de una cadena de Markov con espacio de estados E,
nuestro problema ahora es tratar de estimar las probabilidades de transición de la cadena mediante
el método de máxima verosimilitud. La siguiente proposición nos proporciona las condiciones para la
existencia y unicidad del estimador máximo verosı́mil para la matriz de transición de una cadena de
Markov homogénea con espacio de estados finito.

Proposición 4.1. Dada una trayectoria observada {xn }n∈{0,1,...,m} de una cadena de Markov homogénea e
irreducible X = {Xn }n∈N con espacio de estados E = {1, 2, . . . , N } y matriz de transición P, el estimador
máximo verosı́mil de la matriz de transición está dado por

nij (m)
P̂ (m) = {p̂ij (m) = }i,j∈E .
ni (m)

Donde m denota el horizonte de tiempo, nij (m) es el número de veces en el que la trayectoria pasó del
estado
PN i al j en un paso y ni (m) es el número de visitas al estado i hasta el tiempo m-1, y está dado por
j=1 nij (m) para todo i ∈ E.

Nota 4.1. Una observación que vale la pena mencionar es que la existencia de este estimador está
condicionada a que ni (m) > 0 para todo i ∈ E, que es equivalente a pedir que hasta el tiempo m − 1 el
proceso haya visitado todos los estados, lo cual se logra debido a la irreducibilidad.

Demostración. Supongamos que tenemos observada la trayectoria {X0 = x0 , X1 = x1 , . . . , Xm = xm }

41
42 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

entonces la función de verosimilitud está dada por:

N Y
N
n (m)
Y
L(P ; {xn }m
n=0 ) = P[X0 = x0 , X1 = x1 , . . . , Xm = xm ] = P[X0 = x0 ] pijij .
i=1 j=1

Entonces la correspondiente log-verosimilitud es:

N X
X N
L(P ; {xn }m
n=0 ) = log(P[X0 = x0 ]) + nij (m) log(pij ).
i=1 j=1

El espacio parametral está dado por

N
X
P = {P ∈ MN ×N : pij ≥ 0 y pij = 1 para todo i ∈ E}.
j=1

Ahora maximizaremos la log-verosimilitud en función de las probabilidades de transición mediante mul-

tiplicadores de Lagrange, con las restricciones : N
P
j=1 i,j =1 para todo i ∈ E.
p
Ası́ planteamos el problema de maximización
N X
X N N
X N
X
H(P ) = log(P[X0 = y0 ]) + nij (m) log(pij ) + λi (1 − pij ).
i=1 j=1 i=1 j=1

Derivando con respecto a pij e igualando a cero tenemos las siguientes ecuaciones:
nij (m)
− λi = 0 para todo i, j ∈ E. (4.1)
pij
Entonces
nij (m)
pij = para todo i, j ∈ E.
λi
Sumando sobre j para todo i ∈ E se tiene
N
X ni (m)
1= pij = para todo i ∈ E,
λi
j=1

por tanto
nij (m)
λi = ni (m) y p̂ij (m) = para todo i, j ∈ E.
ni (m)
Hasta aquı́ se vió que {p̂ij (m) = nij (m)/ni (m)}i,j∈E es un punto crı́tico de H en P. Ahora queremos
ver que en este punto hay un máximo. Sabemos que L(P ; {xn }m n=0 ) es una función continua y el espacio
parametral anteriormente definido es un conjunto compacto en el espacio euclidiano, por esto existe P̂ (m)
tal que L(P ; {xn }m m m
n=0 ) ≤ L(P̂ (m); {xn }n=0 ) y L(P̂ (m); {xn }n=0 ) > 0, entonces se tiene para cada P en el
espacio parametral y donde esté bien definida la log-verosimilitud
l(P̂ (m); {xn }m m
n=0 ) ≤ l(P̂ (m); {xn }n=0 ).

Por el Teorema del multiplicador de Lagrange (ver [? ] página 414) P̂ (m) satisface las ecuaciones (4.1) y
por lo anteriormente mostrado se tiene
nij (m)
P̂ (m) = {p̂ij (m) = }i,j∈E ,
ni (m)
es el estimador máximo verosı́mil.
4.1. CADENAS DE MARKOV 43

Ejemplo 4.1. Estimación por máxima verosimilitud de la matriz de transición correspondiente a la tra-
yectoria simulada en el Ejemplo 3.1, la función que se utilizó se encuentra en el Código D.2.
 
0.181818 0.0909091 0.636364 0.0909091
 0.454545 0.181818 0.181818 0.181818 
P̂ (m) = 
 0.166667 0.416667 0.166667
.
0.25 
0.333333 0.5 0.166667 0.0

Ahora mostraremos que el estimador máximo verosı́mil de la matriz de transición cumple propiedades
estadı́sticas muy importantes como son consistencia y eficiencia asintótica.
Denotemos por ν = (ν1 , . . . , νN ) a la distribución estacionaria (si existe) de la cadena de Markov {Xn }n≥0 .

Nota 4.2. Una observación que vale la pena mencionar es que si la cadena {Xn }n∈N es irreducible y
recurrente positiva entonces ν existe y νi es estrictamente positiva para todo i ∈ E.

4.1.2. Propiedades asintóticas

Teorema 4.1. (consistencia). Dada una Cadena de Markov irreducible y recurrente positiva con espacio
n (m)
de estados E y matriz de transición P, los estimadores p̂ij (m) = niji (m) convergen casi seguramente al
verdadero parámetro pij cuando m tiende a infinito para todo i, j ∈ E. Más aún
c.s.
||P − P̂(m)|| −−−−→ 0.
m→∞

Donde || ∗ || denota la norma usual en el espacio euclidiano. En otras palabras el estimador P̂(m) es
consistente.

Demostración. Podemos escribir a nuestro estimador de la siguiente forma:

Pni (m)
nij (m) k=1 I{X1+τ i =j} + I{X0 =i,X1 =j}
k
p̂ij (m) = = .
ni (m) ni (m)
Por el Teorema 2.1 sabemos que {I{X =j} }k≥0 son i.i.d. Sabemos que la cadena {Xn }n≥0 es recurrente
1+τ i
k
positiva, lo que implica que lı́m ni (m) = ∞ y aplicando la ley fuerte de los grandes números obtenemos
m→∞

lı́m p̂ (m) = P[X1+τ i = j].

m→∞ ij 1

Por otro lado sabemos que {Xn }n≥0 es recurrente positiva, lo que implica P[τ1i < ∞] = 1. Por la propiedad
de Markov y debido a que τ1i es un tiempo de paro tenemos:
∞
X
P[X1+τ i = j] = P[X1+τ i = j, τ1i = k]
1 1
k=1
X∞
= P[X1+k = j, τ1i = k]
k=1
∞
X
= pij P[τ1i = k]
k=1
= pij .
44 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Por lo tanto
c.s.
||P − P̂(m)|| −−−−→ 0.
m→∞

Ejemplo 4.2. Se simuló una trayectoria de la cadena de Markov del Ejemplo 3.1 con horizonte m=10000.
En la figura 4.1 se muestra ||P̂ (n) − P || para cada n ≤ 10000, para este análisis se utilizó el Código D.3 .

Figura 4.1: Gráfica de la norma entre P y P̂ (n) respecto al tiempo

Ahora nos concentraremos en demostrar que el estimador máximo verosı́mil de la matriz de transición
P̂ (m) tiene la propiedad de eficiencia asintótica, esta propiedad nos ayudará a hacer estimación por
intervalo de la matriz de transición de una cadena de Markov.
Para lograr este objetivo primero debemos demostrar el siguiente lema.

Lema 4.1. Sean {(gi1 (m), . . . , giN (m))}1≤i≤N una familia independiente de vectores aleatorios sobre RN
tal que para cada i ∈ E

(gi1 (m), . . . , giN (m)) ∼ multinomial(m, pi1 , . . . , piN ).

Entonces
d
G(m) −−−−→ N (0̄, Γ).
m→∞

En donde
gij (m) − mpij
G(m) = { √ }i,j∈E
m
4.1. CADENAS DE MARKOV 45

y 0̄ denota la matriz de tamaño N × N con todas las entradas iguales a cero, mientras que la matriz de
varianzas y covarianzas tiene dimensión N 2 × N 2 y está dada por

Γ = {δik (δjl pij − pij pil )}i,j,l,k∈E .

En donde
(
1 si i = k,
δik =
0 si i 6= k.

g (m)−mp
Demostración. Denotemos por Gi (m) el vector cuyas entradas están dadas por Gi (m)j = ij √m ij para
todo i, j ∈ E.
A continuación probaremos que la función generadora de momentos de Gi (m) converge a la generadora
de momentos de un vector aleatorio que se distribuye normal, y ası́ concluir lo que se desea demostrar.
Sea t̄ = (t1 , . . . , tN ) ∈ RN y h∗, ∗i que denota el producto interno usual en el espacio euclidiano de dos
matrices de igual dimensión.
Si Z = (z1 , . . . , zN ) ∼ multinomial(m, p1 , p2 , . . . , pN ) su función generadora de momentos está dada por
"N #m
X
MZ (t̄) = pk exp(ti ) .
k=1

Por tanto  !#m

N N
"
X 1 X
E [exp (ht̄, Gi (m)i)] =  pij exp √ tj − tk pik  .
m
j=1 k=1

Esta última igualdad se deduce fácilmente de la función generadora de momentos de una multinomial de
parámetros (m, pi1 , . . . , piN ).

Desarrollando la serie de Taylor de la función generadora de momentos E [exp (ht̄, Gi (m)i)] en el punto
(0, . . . , 0) ∈ RN tenemos que

N N N X
N N
" !#
X 1 X X X pij
pij exp √ tj − tk pik =1+ tk1 tk2 (δk j − pik1 )(δk2 j − pik2 )
m 2m 1
j=1 k=1 k1 =1 k2 =1 j=1

N X
N X
N N N
" !#
X X pij 1 X
+ tk1 tk2 tk3 3/2
exp √ tm,j − tm,k pik (δk1 j − pik1 )(δk2 j − pik2 )(δk3 j − pik3 ).
6m m
k1 =1 k2 =1 k3 =1 j=1 k=1

Donde (tm,1 , . . . , tm,N ) = αm t̄ para alguna αm ∈ (0,1) y las derivadas parciales de E [exp (ht̄, Gi (m)i)]
están dadas por
N N N N
" !# " !#
∂ X 1 X X 1 X 1
pij exp √ tj − tk pik = pij exp √ tj − tk pik √ (δjk1 − pik1 ) = 0.
∂tk1 m m m
j=1 k=1 j=1 k=1

N N
" !#
∂2 X 1 X
pij exp √ tj − tk pik
∂tk2 ∂tk1 m
j=1 k=1
46 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

N N
" !#
X 1 X 1 1
= pij exp √ tj − tk pik √ (δjk1 − pik1 ) √ (δjk2 − pik2 ).
m m m
j=1 k=1
y por último
N N
" !#
∂3 X 1 X
pij exp √ tj − tk pik
∂tk3 ∂tk2 ∂tk1 m
j=1 k=1
N N
" !#
X 1 X 1 1 1
= pij exp √ tj − tk pik √ (δjk1 − pik1 ) √ (δjk2 − pik2 ) √ (δjk3 − pik3 ).
m m m m
j=1 k=1

Es fácil darse cuenta que el último sumando de la serie de Taylor de la función E [exp (ht̄, Gi (m)i)] en
valor absoluto puede ser acotado por mM 3/2 para alguna M > 0, por lo anterior concluimos que para toda
m > 0 suficientemente grande se tiene
 m
N X N N
X X p ij M
1 + tk1 tk2 (δk j − pik1 )(δk2 j − pik2 ) − 3/2  ≤ E [exp(ht̄, Gi (m)i)]
2m 1 m
k1 =1 k2 =1 j=1
 m
N X N N
X X pij M
≤ 1 + tk1 tk2 (δk1 j − pik1 )(δk2 j − pik2 ) + 3/2  .
2m m
k =1 k =1
1 2 j=1

Haciendo m −→ ∞ se tiene que

 
N X
N N
X X pij
E [exp(ht̄, Gi (m)i)] −→ exp  tk1 tk2 (δk1 j − pik1 )(δk2 j − pik2 ) .
2
k1 =1 k2 =1 j=1

Este último término se simplifica en

 
N X
N
X pik1
exp  tk1 tk2 (δk2 k1 − pik2 ) .
2
k1 =1 k2 =1

De donde obtenemos
 
1 T i
E [exp(ht̄, Gi (m)i)] −−−−→ exp t̄ Σ t̄ .
m→∞ 2
Donde Σi = {pik1 (δk2 k1 − pik2 )}k1 ,k2 ∈E y t̄T es el vector transpuesto asociado a t̄, por lo anterior y por la
independencia de los vectores {Gi (m)}i∈E podemos afirmar que la matriz G(m) cuyo i-ésimo renglón es
Gi (m) tiende en distribución a una normal con media cero y matriz de varianzas y covarianzas
Γ = {δik (δjl pij − pij pil )}i,j,l,k∈E .

Teorema 4.2. Sea X = {Xn }n∈N una cadena de Markov con matriz de transición P = {pij }i,j∈E y
con espacio de estados E, suponemos además que es irreducible y recurrente positiva, con distribución
estacionaria ν = (ν1 , . . . , νN ), entonces
√2 d
m(P̂(m) − P ) −−−−→ N (0̄, Σ).
m→∞
1
Donde Σ es una matriz de dimensión N 2 × N 2 con entradas Σij,kl = νi δik (δjl pij − pij pil ) para todo
i, j, l, k ∈ E.
4.1. CADENAS DE MARKOV 47

Demostración. Sin perdida de generalidad podemos suponer que la cadena {Xn }n∈N está definida como
en la Definición 3.1, pues como se vio anteriormente esta construcción conserva la misma ley.
Se sabe que (ni1 (m), . . . , niN (m)) es el vector de frecuencias de (wi,1 , . . . , wi,ni (m) ). Ahora denotemos por
(hi1 (m), . . . , hiN (m)) el correspondiente vector de frecuencias de (wi,1 , . . . , wi,dmνi e ), es facil darse cuenta
que sigue una distribución multinomial de parámetros (dmνi e , pi1 , . . . , piN ).
Además (hi1 (m), . . . , hiN (m)) es independiente al vector (hj1 (m), . . . , hjN (m)) si i 6= j.
Por el Lema 4.1 podemos afirmar que
d
H(m) −−−−→ N (0, Γ). (4.2)
m→∞

En donde
hij (m) − dmνi e pij
H(m) = { p }i,j∈E
dmνi e
y la matriz de varianzas y covarianzas está dada por

Γ = {δik (δjl pij − pijpil)}i,j,l,k∈E .

Sabemos que
( ! √ !)
√   nij (m) − pij ni (m) m
m P̂ (m) − P = p p
ni (m) ni (m) i,j∈E
  
nij (m) − pij ni (m) m
= √
m ni (m) i,j∈E
  
nij (m) − pij ni (m) m
= √ − H̃(m)i,j + H̃(m)i,j .
m ni (m) i,j∈E

Donde
hij (m) − dmνi e pij
H̃(m) = { √ }1≤i,j≤N .
m
Supongamos por el momento que

  
nij (m) − pij ni (m) m p
√ − H̃(m)i,j −−−−→ 0 para todo i, j ∈ E.
m ni (m) m→∞

Por lo anterior y por la Proposición B.9 basta con mostrar la convergencia para
  
m
H̃(m)i,j .
ni (m) i,j∈E

Por otro lado

   ( p √ !)
m dmνi e m
H̃(m)i,j = H(m)i,j p p ,
ni (m) i,j∈E ni (m) ni (m) i,j∈E

por el Teorema ergódico (Teorema 2.2)

p √
dmνi e c.s. m c.s. 1
p −−−−→ 1 y p −−−−→ √ ,
ni (m) m→∞ ni (m) m→∞ νi
48 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

por lo anterior, por (4.2) y por la Proposición B.7 concluimos que basta con mostrar la convergencia para
el vector   
1
H(m)i,j √ .
νi i,j∈E

Debido a que H(m) tiende en distribución a un vector aleatorio normal, entonces se deduce que

√
2 d
m(P̂(m) − P ) −−−−→ N (0̄, Σ).
m→∞
1
Donde Σ es una matriz de dimensiones N 2 ×N 2 , con Σ la matriz con entradas Σij,kl = νi δik (δjl pij −pijpil)
para todo i, j, l, k ∈ E.
Ahora mostraremos que
  
nij (m) − pij ni (m) m p
√ − H̃(m)i,j −−−−→ 0.
m ni (m) m→∞

Por la Proposición B.7 y el Teorema ergódico (2.2) basta con probar la convergencia
 
nij (m) − pij ni (m) p
√ − H̃(m)i,j −−−−→ 0.
m m→∞

Sea ε > 0.
Por el Teorema ergódico (2.2) podemos afirmar que existe k tal que para todo m > k se cumple

P[| ni (m) − dmνi e |> mε3 ] < ε.

Denotemos a
n
X
Snij = I{wi,r =j} − pij n.
r=1

Por lo anterior
√ i
 
nij (m) − pij ni (m) h
ij
P | √ − H̃(m)i,j |> ε ≤ ε + P max{| Snij − Sdmνie
|:| n − dmνi e |≤ mε 3
} > ε m
m
 √ 
 ij 3 ε m
≤ ε + 2P max | Sn |: 1 ≤ n ≤ mε > .
2

Esta última desigualdad se debe a que {wi,r }i≥1 son i.i.d.

Además por la desigualdad de Kolmogorov (ver [? ] página 350)
 
nij (m) − pij ni (m)
P | √ − H̃(m)i,j |> ε ≤ ε(1 + 8pij (1 − pij )).
m

Lo anterior se traduce en
  
nij (m) − pij ni (m) m p
√ − H̃(m)i,j −−−−→ 0.
m ni (m) m→∞

Sea X = {Xn }n∈N una cadena de Markov con espacio de estados E y con matriz de transición P =
{pij }i,j∈E . Supongamos que todas las probabilidades de transición son estrictamente positivas, esto con
4.1. CADENAS DE MARKOV 49

el objeto de tener un espacio parametral abierto y ası́ poder definir la matriz de información de Fisher.
Definimos a nuestro espacio parametral como
N
X −1
Θ = {P̃ ∈ MN ×N −1 : pij > 0 para todo i ∈ E , 1 ≤ j ≤ N − 1 y pij < 1 para todo i ∈ E}.
j=1

El espacio parametral anterior ahora es un conjunto abierto con la norma usual, además nos permitió
quitar las dependencias entre las probabilidades de transición correspondientes a cada renglón.

Nota 4.3. Por simplicidad adoptamos la notación L(P ; {xn }mn=0 ), para referirnos a la función de vero-
similitud de la muestra {xn }m
n=0 de la cadena de Markov {Xt } t≥0 con parámetro P̃ .

En lo que resta de esta sección vamos a trabajar especı́ficamente con el espacio parametral anteriormente
definido.

Proposición 4.2. Sea {Xn }n∈N una Cadena de Markov homogénea con espacio de estados E = {1, 2, . . . , N }
y P su correspondiente matriz de transición con todas sus entradas estrictamente positivas.
Entonces la matriz de información de Fisher está dada por
 
1 1
(im (P̃ ))i1 j1 ,i2 j2 = δi1 i2 E[ni1 (m)] δj1 j2 + .
p i1 j1 p i1 N

Donde im (P̃ ) ∈ MN 2 −N ×N 2 −N y los indices corren en i1 , i2 ∈ E y 1 ≤ j1 , j2 ≤ N − 1.

Además si suponemos que la cadena es irreducible y recurrente positiva, la matriz de información es
invertible y su inversa está dada por

δi1 i2 (δj1 j2 pi1 j1 − pi1 j1 pi1 j2 )

(im (P̃ ))−1
i1 j1 ,i2 j2 =
E[ni1 (m)]

con i1 , i2 ∈ E y 1 ≤ j1 , j2 ≤ N − 1.

Demostración. Se sabe que el espacio Θ es abierto y no vacı́o, además se tiene que por la definición de Θ
El conjunto

A = {(x0 , x1 , . . . , xm )|L(P ; {xn }m

n=0 ) > 0}

es independiente del parámetro en el espacio P̃ ∈ Θ, más aún resulta ser E m .

Debido a que L(P ; {xn }mn=0 ) es en realidad un polinomio que depende solamente de potencias del conjunto
{pij }i∈E,1≤j≤N −1 se tiene que
∂ 2 L(P ; {xn }m n=0 )
∂pi1 j1 ∂pi2 j2
existe y es finita para todo i1 , i2 ∈ E y 1 ≤ j1 , j2 ≤ N − 1.
Por estas últimas observaciones se sabe que la matriz de información de Fisher también puede ser calcu-
lada como
 2
∂ log(L(P ; {Xn }m

n=0 ))
im (P̃ )i1 j1 ,i2 j2 = E −
∂pi1 j1 ∂pi2 j2
que resulta ser igual a " #
ni j (m) ni1 N (m)
δi1 i2 E δj1 j2 1 21 + P −1
p i1 j1 (1 − N j=1 pi1 j )
2
50 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

.
Por otra parte ya que X
L(P ; {xn }m
n=0 ) = 1
x0 ,x1 ,...,xm ∈E

la derivada puede introducirse a la suma y en consecuencia

∂log(L(P ; {Xn }m
 
n=0 ))
E =0
∂pi1 j1

para todo i1 ∈ E y 1 ≤ j1 ≤ N − 1 Lo que significa que im (P̃ )i1 j1 ,i2 j2 es igual a la covarianza entre los
siguientes términos

∂log(L(P ; {Xn }m
n=0 ))
∂pi1 j1
y
∂log(L(P ; {Xn }m
n=0 ))
.
∂pi2 j2
Por lo que se concluye que la matriz de información de Fisher es una matriz semi-definida positiva.
Además
" #
ni1 j1 (m) ni1 N (m)
im (P̃ )i1 j1 ,i2 j2 = δi1 i2 E δj1 j2 + P −1
p2i1 j1 (1 − N j=1 pi1 j )
2
" #
1 1
= δi1 i2 E[ni1 (m)] δj1 j2 + P −1 .
pi1 j1 (1 − N j=1 pi1 j )

h i
Sean {Ai }i∈E matrices de dimensión N − 1 × N − 1 tales que Aij1 j2 = δj1 j2 pij1 + piN
1
para todo i ∈ E y
1
0 < j1 , j2 < N , entonces se tiene

E[n1 (m)]A1
 
0 ··· 0
 0 E[n2 (m)]A2 · · · 0 
im (P̃ ) =  .
 
.. .. .. ..
 . . . . 
0 0 · · · E[nN (m)]A N

Ahora queremos probar que esta matriz es invertible, basta con probar que las matrices {Ai }i∈E lo son.
Realizando operaciones elementales se muestra que la matrices {Ai }i∈E son invertibles y además

(Ai )−1
j1 j2 = (δj1 j2 pij1 − pij1 pij2 )

para todo i ∈ E y 0 < j1 , j2 < N

Por lo anterior la matriz de información de Fisher es definida positiva ya que es semi-definida positiva e
invertible. Por último concluimos que
δi1 i2 (δj1 j2 pi1 j1 − pi1 j1 pi1 j2 )
(im (P̃ ))−1
i1 j1 ,i2 j2 =
E[ni1 (m)]
para todo i1 , i2 ∈ E y 1 ≤ j1 , j2 ≤ N − 1.
Nótese que para que la inversa este bien definida necesitamos que

E[ni1 (m)] > 0 para todo i1 ∈ E.

Esto se obtiene por las hipótesis de recurrencia positiva e irreducibilidad cuando m es grande.
4.2. PROCESOS DE SALTOS DE MARKOV 51

Teorema 4.3. (Eficiencia Asintótica) Sea {Xn }n∈N una Cadena de Markov homogénea con espacio
de estados E = {1, 2, . . . , N } suponemos además que es irreducible y recurrente positiva entonces

√
m(P̃ˆ (m) − P̃ ) −−−−→ N (0̃, Σ̃)
2 d
m→∞

Donde
Σ̃ = lı́m m(im (P̃ ))−1 ,
m→∞

P̃ˆ (m) es la matriz P̂ (m) quitándole la última columna y 0̃ es la matriz de dimensión N × (N − 1) con
todas sus entradas iguales a cero.

Demostración. La prueba de este Teorema se sigue de los Teoremas 2.2, 4.2 y la Proposición 4.2.

Ejemplo 4.3. Se simuló una cadena de Markov con espacio de estados E = {1, 2}, distribución inicial

π = (0.5, 0.5)

y matriz de transición
 
0.2 0.8
P = ,
0.8 0.2

con horizonte m = 100000.

Es fácil mostrar que (0.5, 0.5) es la única distribución estacionaria, por lo tanto
 
0.32 0.0
Σ̃ = .
0.0 0.32

Mientras que la inversa de la matriz de información de Fisher estimada esta dada por

−1
   
0.32075 0.0
m i\
m (P̃ ) = .
0.0 0.323248

Para la estimación anterior se utilizó la función del Código D.17

4.2. Procesos de saltos de Markov

4.2.1. Con información a tiempo continuo.

Nuestro objetivo en esta sección es estimar el generador infinitesimal dada una trayectoria de un PSM,
ası́ como también demostrar algunas propiedades estadı́sticas del estimador de máxima verosimilitud.
En secciones pasadas estudiamos la importancia del generador infinitesimal, por lo cual nuestro objetivo
será estimar sus parámetros vı́a máxima verosimilitud.
A partir de ahora sólo consideraremos procesos cuyos generadores infinitesimales tengan entradas en la
diagonal estrictamente negativas, esto implica que los tiempos de saltos sean finitos casi seguramente. La
existencia de estos procesos se probó en el Teorema 3.2.
52 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Estimación máximo verosı́mil

La siguiente proposición demuestra la existencia del estimador máximo verosı́mil correspondiente a una
observación en tiempo continuo {xt }0≤t≤T de un PSM.

Proposición 4.3. Supongamos que tenemos una trayectoria observada en tiempo continuo {xt }0≤t≤T de
un PSM {Xt }t≥0 con espacio de estados E y generador infinitesimal Q, entonces el estimador máximo
verosı́mil del generador infinitesimal dado que el proceso inicio en el estado x0 está dado por :

Q̂(T ) = {q̂xy (T )}x,y∈E

donde
(
Nxy (T )/Rx (T ) para x 6= y
q̂xy (T ) =
−Nx (T )/Rx (T ) para x = y.

Donde T es el horizonte de tiempo, Q̂(T ) denota el estimador máximo verosı́mil del generador infinitesi-
mal, Nxy (T ) es el número de veces observadas que el proceso saltó del estado x al estado y en un tiempo
RT
Pdt, es el tiempo total observado en el cual el proceso se mantuvo
de salto, mientras que Rx (T ) = 0 I{Xt =x}
en el estado x, denotemos por Nx (T ) = y6=x Nxy (T ).

Nota 4.4. Una observación que vale la pena mencionar es que la existencia de este estimador está
condicionada a que Nx (T ) > 0 para todo x ∈ E.

Demostración. Asumimos que el proceso tuvo n saltos, sean 0 = t0 < t1 < . . . < tn los tiempos observados
de salto hasta el tiempo T . Por los Teoremas 2.6, 2.7 y 2.8 la función de verosimilitud es igual a

L(Q; {xt }0≤t≤T ) = fT1 ,XT1 ,T2 ,XT2 ,...,Tn ,XTn ,NT |X0 (t1 , xt1 , t2 , xt2 , . . . , tn , xtn , n|x0 )
= fT1 ,XT1 ,T2 ,XT2 ,...,Tn ,XTn |X0 (t1 , xt1 , t2 , xt2 , . . . , tn , xtn |x0 )P[Tn+1 − Tn > T − tn |XTn = xtn ]
= (λx0 e−λx0 t1 px0 xt1 )fT1 ,XT1 ,T2 ,XT ,...,Tn−1 ,XTn−1 |X0 (t2 , xt2 , . . . , tn , xtn |xt1 ) exp(−λxtn (T − tn )).

Por lo que obtenemos

n
Y
L(Q; {xt }0≤t≤T ) = [ (λxti−1 exp (−λxti−1 (ti − ti−1 ))pxti−1 xti )] exp(−λxtn (T − tn ))
i=1
N (T )
Y Y
= λN
x
x(T )
exp(−λx Rx (T )) pxyxy
x∈E y6=x
YY
= (λx pxy )Nxy (T ) exp(−λx pxy Rx (T ))
x∈E y6=x
N (T )
YY
= qxyxy exp(−qxy Rx (T )).
x∈E y6=x

El espacio parametral está dado por

X
Q = {Q ∈ MN ×N : qxy ≥ 0 para todo y 6= x y qxy = −qxx > 0 para todo x ∈ E}.
y6=x
4.2. PROCESOS DE SALTOS DE MARKOV 53

Por otro lado la función de log-verosimilitud es

XX
l(Q; {xt }0≤t≤T ) = Nxy (T ) log(qxy ) − qxy Rx (T ).
x∈E y6=x

La última expresión depende sólo de las entradas fuera de la diagonal de la matriz Q, esto quiere decir
que no tenemos ninguna restricción adicional a la no negatividad sobre las {qxy }x,y∈E,x6=y , por lo cual
derivando la log-verosimilitud e igualando a cero tenemos como punto crı́tico

(
Nxy (T )/Rx (T ) para x 6= y
q̂xy (T ) = (4.3)
−Nx (T )/Rx (T ) para x = y.

Hasta aquı́ se mostró que

Q̂(T ) = {q̂xy (T )}x,y∈E (4.4)

es un punto crı́tico de l(Q; {xt }0≤t≤T ) sobre Q, ahora sólo falta probar que es un punto máximo. Nótese
que
X X X
l(Q; {xt }0≤t≤T ) = ( Nxy (T ) log(qxy )) + ( qxx Rx (T )), (4.5)
x∈E y∈Ex x∈E

donde Ex = {y 6= x : Nxy (T ) > 0}. P

Por (4.5) sabemos que si Q es tal que qxy > 0 para algún y ∈ E −Ex , el valor de x∈E qxx Rx (T ) disminuye
ya que qxx se hace más pequeño.
Por lo anterior qxy = 0 para todo x, y ∈ E tal que y ∈ E − Ex , en resumen maximizar l(Q; {xt }0≤t≤T ) se
reduce a maximizar en el conjunto {qxy : x, y ∈ E, y ∈ Ex }.
Por otro lado para mostrar que (4.4) es un punto máximo de l(Q; {xt }0≤t≤T ), por el Teorema 6.9.4 (página
364 de [? ]) es suficiente mostrar que

∂ 2 l(Q̂(T ); {xt }0≤t≤T )

{ }x1 ,x2 ∈E,y1 ∈Ex1 ,y2 ∈Ex2 ,
∂qx1 y1 ∂qx2 y2
es una matriz diagonal y que en la diagonal sus elementos son estrictamente negativos, desarrollando las
derivadas de segundo orden obtenemos

∂ 2 l(Q̂(T ); {xt }0≤t≤T ) −Nx1 y1 (T )

= δx1 x2 δy1 y2 , para todo x1 , x2 ∈ E, y1 ∈ Ex1 , y2 ∈ Ex2 .
∂qx1 y1 ∂qx2 y2 (q̂x1 y1 (T ))2

Lo anterior prueba que Q̂(T ) de la expresión (4.4) es el estimador máximo verosı́mil.

Ejemplo 4.4. Estimación por máxima verosimilitud del generador infinitesimal correspondiente a la tra-
yectoria simulada en el Ejemplo 3.3, la función que se utilizó se encuentra en el Código D.5.
 
−0.144964 0.0362409 0.0181205 0.0906023
 0.0304693 −0.182816 0.0761733 0.0761733 
Q̂(T ) =  .
 0.187188 0.140391 −0.655159 0.32758 
0.0172995 0.121096 0.138396 −0.276792

A continuación mostraremos algunas propiedades del estimador de máxima verosimilitud del generador
infinitesimal de un PSM dado, dichas caracterı́sticas prueban la importancia de utilizar este estimador.
El siguiente teorema prueba que nuestro estimador converge al verdadero parámetro conforme la muestra
crece.
54 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Propiedades asintóticas

Teorema 4.4. (Consistencia) Sea un PSM {Xt }t≥0 ergódico con generador infinitesimal Q y con
espacio de estados E, el estimador máximo verosı́mil del generador infinitesimal Q̂(T ) dado por Q̂(T ) =
{q̂xy (T )}x,y∈E converge casi seguramente al verdadero parámetro Q cuando T(horizonte) tiende a infinito.
En otras palabras el estimador es consistente.

Demostración. Por los Teoremas 2.11 y 2.13 concluimos

Nxy (T ) T Nxy (T ) c.s.
= −−−−→ qxy .
Rx (T ) Rx (T ) T T →∞

Ya que el espacio de estados es finito se puede afirmar lo siguiente

c.s.
||Q − Q̂(T )|| −−−−→ 0 cuando T −→ ∞.
T →∞

Ejemplo 4.5. Se simuló una trayectoria del PSM del Ejemplo 3.3 con horizonte T=10000.
En la figura 4.2 se muestra ||Q̂(t) − Q|| para cada t ≤ 10000, para este análisis se utilizó el Código D.6.

Figura 4.2: Gráfica de la norma entre Q y Q̂(T ) respecto al tiempo

Al igual que en el caso de las cadenas de Markov a tiempo discreto, es posible probar en el caso de PSM
que su estimador máximo verosı́mil asociado al generador infinitesimal, es asintóticamente eficiente. Para
mostrar esta propiedad, requerimos el siguiente lema.
4.2. PROCESOS DE SALTOS DE MARKOV 55

Lema 4.2. Sean 0 = T0 < T1 < . . . < Tn < ∞ los tiempos de salto de un PSM {Xt }t≥0 con espacio de
estados E y generador infinitesimal Q, entonces

E[T − TNT ]
√ −−−−→ 0.
T T →∞

En particular por la Proposición B.3 tenemos convergencia en probabilidad.

Demostración.
∞
X
E[T − TNT ] = T P[NT = 0] + E[(T − TNT )I{NT =n} ]
n=1
X∞
= T P[NT = 0] + E[(T − Tn )I{Tn ≤T <Tn+1 } ]
n=1
X∞ Z T Z ∞
= T P[NT = 0] + (T − u)fTn ,Tn+1 (u, v)dvdu
n=1 0 T
X∞ X Z T Z ∞
= T P[NT = 0] + (T − u)fTn ,XTn (u, x)fT0 |X0 =x (v − u)dvdu.
n=1 x∈E 0 T

La igualdad anterior se debe al Teorema 2.8, dicho lo anterior integrando y sumando sobre x ∈ E
obtenemos
XZ T Z ∞ XZ T
(T − u)fTn ,XTn (u, x)fT0 |X0 =x (v − u)dvdu = (T − u)fTn ,XTn (u, x) exp−λx (T −u) du
x∈E 0 T x∈E 0
Z T
≤ (T − u)fTn (u) exp−α(T −u) du.
0

Donde α=mı́n{λx : x ∈ E}, ahora integramos por partes el término anterior y obtenemos

Z T
(T − u)fTn (u) exp−α(T −u) du =
0

1
Z T−α Z T
−FTn (u) exp−α(T −u) [α(T − u) − 1]du + −FTn (u) exp−α(T −u) [α(T − u) − 1]du. (4.6)
1
0 T−α

El primer término de (4.6) es negativo mientras que el segundo es positivo, por esta razón se sigue del
Lema 2.1 la siguiente desigualdad, donde β =máx{λx : x ∈ E}
∞ Z
X T Z T
(T − u)fTn (u) exp−α(T −u) du ≤ βu exp−α(T −u) [1 − α(T − u)]du.
1
n=1 0 T−α

Por esta desigualdad

Z T
E[T − TNT ] ≤ T P[NT = 0] + βu exp−α(T −u) [1 − α(T − u)]du. (4.7)
1
T−α
56 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Integrando por partes el segundo miembro de la desigualdad (4.7) tenemos

Z T Z T
−α(T −u) exp(−1)
βu exp [1 − α(T − u)]du = + β(T − u) exp−α(T −u) du
T−α1 α T−α 1

Z 1
exp(−1) α
= + βv exp−αv dv.
α 0

El segundo sumando de la última igualdad es menor o igual que la esperanza de una exponencial de
parámetro α multiplicada por αβ , por esta razón tenemos
X exp(−1) +β/α
E[T − TNT ] ≤ T P[X0 = 0] exp−λx T + .
α
x∈E

Cocluimos la convergencia
E[T − TNT ]
√ −−−−→ 0.
T T →∞

Por último debido a la Proposición B.2 tenemos

T − TNT p
√ −−−−→ 0.
T T →∞

Teorema 4.5. Considere un PSM {Xt }t≥0 ergódico con espacio de estados E y generador infinitesimal
Q, entonces
√
2 ˆ ) − Q̃) −−− d
T (Q̃(T −→ N (0̃, Λ).
T →∞
ˆ ) y Q̃ son el estimador y el generador infinitesimal sin la diagonal, es decir matrices de
Donde Q̃(T
dimensión N × (N − 1) y Λ ∈ M(N 2 −N )×(N 2 −N ) es la matriz con entradas Λx1 y1 ,x2 y2 = νx1 δx1 x2 δy1 y2 qx1 y1
1
para todo x1 , y1 , x2 , y2 ∈ E con x1 6= y1 y x2 6= y2 .

Demostración. Sin pérdida de generalidad suponemos que el proceso {Xt }t≥0 está construido como en la
Definición 3.3.

√ √
  
ˆ ) − Q̃) = Nxy (T )
T (Q̃(T T − qxy
Rx (T ) x∈E,y6=x
√
  
Nxy (T ) Nxy (T ) Nxy (T )
= T − λx + λx − qxy
Rx (T ) Nx (T ) Nx (T ) x∈E,y6=x
√
     
Nxy (T ) Nx (T ) Nxy (T )
= T − λx + λx − pxy
Nx (T ) Rx (T ) Nx (T ) x∈E,y6=x
    
T Nxy (T ) Nx (T ) − λx Rx (T ) λx T Nxy (T ) − pxy Nx (T )
= √ + √ .
Nx (T )Rx (T ) T Nx (T ) T x∈E,y6=x

Supongamos por el momento las siguientes afirmaciones:

PdT νx λx e Sx,r
Nx (T ) − λx Rx (T ) dT νx λx e − λx r=0 λx p
| √ − √ |−−−−→ 0. (4.8)
T T T →∞
4.2. PROCESOS DE SALTOS DE MARKOV 57

PdT νx λx e
Nxy (T ) − pxy Nx (T ) r=1 I{wx,r =y} − pxy dT νx λx e p
| √ − √ |−−−−→ 0. (4.9)
T T T →∞

PdT ν λ e Sx,r PdT νx λx e

dT νx λx e−λx r=0x x I{wx,r =y} −pxy dT νx λx e
Definimos ξ x (T ) = √
T
λx
, ρxy (T ) = r=1
√
T
, también
    
T T Nxy (T ) Nx (T ) − λx Rx (T ) x λx T Nxy (T ) − pxy Nx (T ) xy
Y = √ − ξ (T ) + √ − ρ (T )
Nx (T )Rx (T ) T Nx (T ) T x∈E,y6=x

y  
T Nxy (T ) λx T
ZT = (ξ x (T )) + (ρxy (T )) .
Nx (T )Rx (T ) Nx (T ) x∈E,y6=x
√ ˆ
Bajo la notación anterior tenemos que T (Q̃(T ) − Q̃) = Y T + Z T .
Por otro lado por los Teoremas 2.11 y 2.13 sabemos que

Nxy (T ) c.s. T c.s. 1 T c.s. 1

−−−−→ pxy , −−−−→ , −−−−→ ,
Nx (T ) T →∞ Rx (T ) T →∞ νx Nx (T ) T →∞ νx λx
p
y debido al Lema B.1 y las afirmaciones (4.8) y (4.9) tenemos que Y T −−−−→ 0. y por la Proposición B.9
T →∞
sabemos que es suficiente verificar la convergencia para Z T .

Establescamos la siguiente notación

    
T T Nxy (T ) pxy x λx T 1 xy
Z1 = − ξ (T ) + − ρ (T )
Nx (T )Rx (T ) νx Nx (T ) νx x∈E,y6=x

y     
pxy 1
Z2T = x
ξ (T ) + xy
ρ (T ) .
νx νx x∈E,y6=x
p
Bajo la notación anterior tenemos que Z T = Z1T + Z2T , a continuación mostraremos que Z1T −−−−→ 0 y
T →∞
por la Proposición B.9 bastarı́a con mostrar la convergencia para Z2T .

Por la Definición 3.3 sabemos que {Sx,r }r≥1 y {wi,r }i≥1 son independientes y por el Teorema central
del lı́mite tenemos que
d
ξ x (T ) −−−−→ N (0, νx λx )
T →∞
y
d
ρxy (T ) −−−−→ N (0, νx λx pxy ).
T →∞

Por otro lado por los Teoremas 2.11 y 2.13 y por la Proposición B.7 tenemos que
p
Z1 (T ) −−−−→ 0.
T →∞

Por lo anterior basta ver la convergencia de Z2 (T ) y por la Definición 3.3 sabemos que {Sx,r }r≥1 y
{wi,r }i≥1 son independientes, por tanto ξ x (T ) es independiente a ρxy (T ). Por otro lado por el Teorema
de Curtiss (B.3) basta analizar la convergencia para los siguientes dos términos
  !
N X N PdT νx λx e
X pxy dT νx λx e − r=0 Sx,r 
E exp  uxy √
νx T
x∈E y6=x
58 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

y  
N X
N PdT νx λx e !
X 1 r=1 I{wx,r =y} − pxy dT νx λx e
E exp  uxy √  .
νx T
x∈E y6=x

Donde U = {uxy }x,y∈E,x6=y es cualquier matriz de dimensión N × (N − 1). El primer término se puede
escribir como PdT ν λ e p
Y dT νx λx e − r=0x x Sx,r X pxy dT νx λx e
E[exp(( p ) uxy √ )]
x∈E
dT νx λx e y6=x
νx T

el cual por la Proposición B.9 y el Teorema B.2 converge a

λx P
!
2
νx ( y6=x uxy pxy )
Y
exp . (4.10)
2
x∈E

Mientras que el segundo término puede ser escrito como

  PdT νx λx e !
N X N r
X λx r=1 I{wx,r =y} − pxy dT νx λx e
E exp  uxy p  .
νx dT νx λx e
x∈E y6=x

Por el Lema 4.1 converge al número

 
Y X X λx pxk1
exp  uxk1 uxk2 (δk2 k1 − pxk2 ) . (4.11)
νx 2
x∈E k1 6=x k2 6=x

Uniendo (4.10) y (4.11) se concluye que

√ ˆ ) − Q̃) −−−
d
T (Q̃(T −→ N (0̃, Λ).
T →∞

Donde Λx1 y1 ,x2 y2 = νx1 δx1 x2 δy1 y2 qx1 y1 para todo x1 , y1 , x2 , y2 ∈ E con x1 6= y1 y x2 6= y2 .
1
Ahora sólo resta probar

PdT ν λ e Sx,r
Nx (T )−λx Rx (T ) dT νx λx e−λx r=0x x λx p
1. | √
T
− √
T
|−−−−→ 0.
T →∞
PdT νx λx e
Nxy (T )−pxy Nx (T ) r=1 I{wx,r =y} −pxy dT νx λx e p
2. | √
T
− √
T
|−−−−→ 0.
T →∞

Mostraremos la primer afirmación. Sabemos que

Nx (T )
X Sx,r
Rx (T ) = + (T − TNT )I{XT =x} .
λx NT
r=0

Por el Lema 4.2 tenemos PNx (T ) Sx,r

Rx (T ) − r=0 λx T − TNT p
| √ |≤ √ −−−−→ 0.
T T T →∞

Por esto último basta probar

PNx (T ) Sx,r PdT νx λx e Sx,r
Nx (T ) − λx r=0 λx dT νx λx e − λx r=0 λx p
| √ − √ |−−−−→ 0.
T T T →∞
4.2. PROCESOS DE SALTOS DE MARKOV 59

Sea ε > 0.
Por el Teorema 2.13 podemos afirmar que existe k tal que para todo T > k se cumple

P[| Nx (T ) − dT νx λx e |> T ε3 ] < ε.

Denotemos a
n
X
Vnxy = n − Sx,r .
r=0
Por lo anterior

" PNx (T ) PdT ν λ e #

Nx (T ) − r=0 Sx,r dT νx λx e − r=0x x Sx,r
P | √ − √ |> ε ≤
T T

xy
√
ε + P[max{| Vnxy − VdT 3
νx λx e |:| n − dT νx λx e |≤ T ε } > ε T ]

que a su vez debido a que {Sx,r }r≥1 son i.i.d.

" √ #
ε T
≤ ε + 2P max{| Vnxy |: 1 ≤ n ≤ T ε3 } > .
2

Además por la desigualdad de Kolmogorov (ver [? ] página 350)

" PNx (T ) PdT ν λ e #
Nx (T ) − r=0 Sx,r dT νx λx e − r=0x x Sx,r
P | √ − √ |> ε ≤ 9ε.
T T

Lo anterior se traduce en
PdT νx λx e Sx,r
Nx (T ) − λx Rx (T ) dT νx λx e − λx r=0 λx p
| √ − √ |−−−−→ 0.
T T T →∞

Ahora sólo resta probar la segunda afirmación.

Sea ε > 0. Por el Teorema 2.13 podemos afirmar que existe k tal que para todo T > k se cumple

P[| Nx (T ) − dT νx λx e |> T ε3 ] < ε.

Denotemos a
n
X
Snxy = I{wx,r =y} − pxy n.
r=1
Por lo anterior

" PdT νx λx e #
Nxy (T ) − pxy Nx (T ) r=1 I{wx,r =y} − pxy dT νx λx e
P | √ − √ |> ε
T T
es menor o igual que
xy
√
ε + P[max{| Snxy − SdT νx λx e |:| n − dmνi e |≤ T ε3
} > ε T]

a su vez dedido a que {wx,r }r≥1 es i.i.d. es menor o igual que

√
ε T 3
ε + 2P[max{| Sn |: 1 ≤ n ≤ T ε } > ].
2
60 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Esta última desigualdad se debe a que {wi,r }i≥1 son i.i.d.

Además por la desigualdad de Kolmogorov (ver [? ] página 350)
" PdT νx λx e #
Nxy (T ) − pxy Nx (T ) r=1 I{wx,r =y} − p xy dT νx λ x e
P | √ − √ |> ε ≤ ε(1 + 8pij (1 − pij )).
T T

Lo anterior se traduce en
PdT νx λx e
Nxy (T ) − pxy Nx (T ) r=1 I{wx,r =y} − pxy dT νx λx e p
| √ − √ |−−−−→ 0.
T T T →∞

Sea un PSM X = {Xt }t≥0 con espacio de estados E y con generador infinitesimal Q = {qxy }x,y∈E .
Supongamos que todas las probabilidades de transición {pxy }x∈E,y∈E−{x} son estrictamente positivas,
esto con el objeto de tener un espacio parametral abierto y ası́ poder definir la matriz de información de
Fisher.
Definimos a nuestro espacio parametral como

Θ = {Q̃ ∈ MN ×N −1 : qxy > 0 para todo x ∈ E , y ∈ E − {x}}.

El espacio parametral anterior ahora es un conjunto abierto con la norma usual, además nos permitió
quitar las dependencias entre la diagonal del generador infinitesimal y las otras entradas de la matriz.

Nota 4.5. Por simplicidad adoptamos la notación L(Q; {xt }0≤t≤T ), para referirnos a la función de ve-
rosimilitud de la muestra {xt }0≤t≤T del PSM {Xt }t≥0 con parámetro Q̃.

En lo que resta de esta sección vamos a trabajar especı́ficamente con el espacio parametral anteriormente
definido.
Proposición 4.4. Sea un PSM {Xt }t≥0 con espacio de estados E = {1, 2, . . . , N } y Q su correspondiente
generador infinitesimal, entonces la matriz de información de Fisher está dada por
E[Nx1 y1 (T )]
(iT (Q̃))x1 y1 ,x2 y2 = δx1 x2 δy1 y2 .
qx21 ,y1

Donde iT (Q̃) ∈ MN 2 −N ×N 2 −N y los ı́ndices corren en x1 , x2 ∈ E y y1 ∈ E − {x1 }, y2 ∈ E − {x2 }.

Además si suponemos que el proceso es ergódico, la matriz de información es invertible y su inversa está
dada por
qx21 y1
(iT (Q̃))−1
x1 y1 ,x2 y2 = δ δ
x1 x2 y1 y2 .
E[Nx1 y1 (T )]
Donde x1 , x2 ∈ E y y1 ∈ E − {x1 }, y2 ∈ E − {x2 }.

Demostración. Se sabe que el espacio Θ es abierto y no vacı́o, además se tiene que por la definición de Θ
el conjunto

A = {{xt }0≤t≤T |L(Q; {xt }0≤t≤T ) > 0}

es independiente del parámetro Q̃ ∈ Θ.
Debido a que Y Y Nxy (T )
L(Q; {xt }0≤t≤T ) = qxy exp(−qxy Rx (T ))
x∈E y6=x
4.2. PROCESOS DE SALTOS DE MARKOV 61

se tiene que
∂ 2 L(Q; {xt }0≤t≤T )
∂qx1 y1 ∂qx2 y2
existe y es finita para todo x1 , x2 ∈ E, y1 ∈ E − {x1 }, y2 ∈ E − {x2 }.
Por estas últimas observaciones se sabe que la matriz de información de Fisher también puede ser calcu-
lada como
 2 
∂ log(L(Q; {Xt }0≤t≤T ))
(iT (Q̃))x1 y1 ,x2 y2 = E −
∂qx1 y1 ∂qx2 y2
por la demostración de la Proposición 4.3
 
Nx y (T ) E[Nx1 y1 (T )]
(iT (Q̃))x1 y1 ,x2 y2 = E δx1 x2 δy1 y2 12 1 = δx1 x2 δy1 y2 .
qx1 y1 qx21 ,y1
Por lo que se concluye que la matriz de información de Fisher es semi-definida positiva, pues es una matriz
diagonal y en la diagonal sus entradas son no negativas.
Gracias a la forma de la matriz de información de Fisher tenemos que
qx21 y1
(iT (Q̃))−1
x1 y1 ,x2 y2 = δx1 x2 δy1 y2 .
E[Nx1 y1 (T )]
Nótese que para que la inversa este bien definida necesitamos que
E[Nx1 y1 (T )] > 0 para todo x1 ∈ E, y1 ∈ E − {x1 }.
Esto se obtiene por las hipótesis de ergodicidad cuando T es grande.

Teorema 4.6. (Eficiencia Asintótica) Sea un PSM {Xt }t≥0 con espacio de estados E = {1, 2, . . . , N }
y Q su correspondiente generador infinitesimal, suponemos además que es ergódico entonces
√2 ˆ ) − Q̃) −−−
d
T (Q̃(T −→ N (0̃, Λ).
T →∞

Donde
Λ = lı́m T (iT (Q̃))−1 .
T →∞

Demostración. La prueba de este teorema se obtiene de aplicar el Teorema 4.5, la Proposición 4.4 y el
Teorema 2.13.

Ejemplo 4.6. Se simuló un PSM con espacio de estados E = {1, 2, 3}, distribución inicial
π = (1/3, 1/3, 1/3)
y generador infinitesimal  
−2 1 1
Q =  1 −2 1  ,
1 1 −2
con horizonte T = 10000.
Es fácil mostrar que (1/3, 1/3, 1/3) es la única distribución estacionaria, por lo tanto
 
3 0 0 0 0 0
 0 3 0 0 0 0 
 
 0 0 3 0 0 0 
Λ=  0
.
 0 0 3 0 0  
 0 0 0 0 3 0 
0 0 0 0 0 3
62 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV

Mientras que la inversa de la matriz de información de Fisher estimada esta dada por
 
3.13155 0 0 0 0 0
 0 2.96953 0 0 0 0 
−1
   
0 0 3.00593 0 0 0
T i\
 
T ( Q̃) = .

 0 0 0 2.92098 0 0 

 0 0 0 0 3.09605 0 
0 0 0 0 0 3.04055

Para la estimación anterior se utilizó la función del Código D.18

4.2.2. Con información a tiempo discreto.

En esta sección estudiaremos un método para generar trayectorias de un PSM condicionado a puntos
finales. Además analizaremos la rapidez y precisión del método. Usaremos el método para estimar el
generador infinitesimal de un PSM a través del método Monte Carlo de una Cadena de Markov y la
Estimación de Máxima Verosimilitud.

Método Rejection Sampling REJ

Empezaremos con la definción de Puente de Markov nedesaria para la explicación de esta sección.
Definición 4.1. Un puente de Markov con parámetros s, a, b es un proceso estocástico con parámetro
t ∈ [0, s] y con la misma distribución de un Proceso de Saltos de Markov {Xt }st=0 condicionado a X0 = a
y Xs = b, con a, b ∈ E y s > 0.

Especı́ficamente trayectorias que no cumplieron la condición este algoritmo simula una trayectoria del
PSM condicionando la posición inicial X0 = a y termina hasta generar una trayectoria con posición final
Xs = b, rechazando todas aquellas.

Puentes de Markov con parámetros s, a, b.

Algoritmo 4.1. 1. Hacemos X0 = a.

2. Simulamos una trayectoria Xs ∼ P SM (s) (PSM con tiempo de horizonte s, generador infinitesimal
Q y distribución inicial α).

3. Si Xs = b se termina, de lo contrario pasamos a 2.

Ya que el algoritmo rechaza trayectorias hasta encontrar una trayectoria que cumpla las condiciones men-
cionadas anteriormente, entonces es importante mencionar la siguiente proposición:

Proposición 4.5. En el método REJ, la probabilidad de alcanzar el estado final observado es pab (T ).

1. Si T es grande entonces pab (T ) ≈ πb .

2. Si T es pequeño entonces pab (T ) ≈ qab T.

Tendremos en consideración los tiempos de ejecución excesivos, ya que muchas rutas se rechazan si:
4.2. PROCESOS DE SALTOS DE MARKOV 63

T es grande y πb es pequeño, o,

T es pequeño y a 6= b.

Ası́, junto al código D.1 visualizaremos estos casos con los siguientes ejemplos. En el primer caso conside-
ramos una trayectoria con T = 1000 y para que π1 sea pequeño bastara hacer el generador infinitesimal

−56.00 5 3 2 3 8 7 9 10 9.00
 
 0.00
 −43 9 7 2 7 4 4 2 8.00 

 0.00
 10 −56 8 5 7 9 8 8 1.00 

 0.00
 4 10 −55 10 6 6 9 1 9.00 

 0.00 1 8 8 −41 1 6 7 8 2.00 
Q= 
 0.00
 1 4 9 4 −42 2 7 6 9.00 

 0.00
 2 1 1 7 5 −40 9 7 8.00 

 0.00
 3 7 3 7 5 5 −41 2 9.00 

 0.00 3 7 5 6 10 1 2 −42 8.00 
0.05 3 3 3 7 8 9 7 5 −45.05

y π = (1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10). Para el segundo caso donde el tiempo es
pequeño y a distinto de b, consideremos un puente de Markov con T =0.0005, a = 1 y b = 2.
64 CAPÍTULO 4. INFERENCIA SOBRE PROCESOS DE MARKOV
Capı́tulo 5

Ecuaciones diferenciales estocásticas

(EDE)

5.1. Movimiento Browniano o proceso de Wiener

En 1828 el botánico Robert Brown observó que los granos de polen suspendidos en una cierta substancia
presentan movimientos irregulares (véase [5]). Con la contribución del trabajo presentado por Albert Eins-
tein en 1905 (véase [8]), se explicó que este movimiento es producto de las múltiples colisiones aleatorias
de las moléculas de la substancia con los granos de polen. Se puede decir que un movimiento de este tipo
tiene las siguientes caracterı́sticas:

1. Es un movimiento continuo.

2. Tiene desplazamientos independientes en intervalos de tiempo disjuntos.

3. Los desplazamientos pueden modelarse por variables aleatorias gausianas (teorema central del lı́mite,
véase [9]).

En 1923 el matemático Norbert Wiener demostró que existe un proceso estocástico con estas propiedades.

Definición 5.1. Un proceso de Wiener (Movimiento Browniano) estándar unidimensional es un

proceso estocástico W = {Wt }t≥0 tal que

1. W0 = 0 casi seguramente.

2. Tiene trayectorias continuas.

3. Tiene incrementos independientes.

4. La variable aleatoria Wt − Ws tiene distribución N (0, t − s) para 0 ≤ s < t.

Pruebas a detalle de la existencia un proceso con estas caracterı́sticas se pueden encontrar en [4].

65
 66 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

5.1.1. Simulación de trayectorias del movimiento Browninano

Desde que los incrementos de un movimento Browniamo siguen una distribución Normal, son indepen-
dientes y las propiedades de la distribución Nomal, podemos simular una trayectia de un movimiento
Browniano en un intervalo de tiempo [0, T ] con T > 0 usando el siguiente algoritmo
Algoritmo 5.1. Dividimos al intervalo [0, T ] en n subintervalos de la misma longitud ∆ = ti − ti−1 para
i = 1, . . . , n donde 0 = t0 < t1 < . . . < tn = T . Hacemos i = 0, Wi = 0 y luego

1. Generamos x ∼ N(0, 1).

2. i = i + 1.
√
3. Hacemos Wi = Wi−1 + x ∆.
4. Si i = n paramos, en otro caso regresamos al paso 1.

El Código 5.1 genera trayectorias de un movimento Brownoiano.

Código 5.1.
1 # ##################################
2 # Drawn Paths of Brownian Motion
3 # ##################################
4 BM = function (n , Delta )
5
6 trayMB = zeros ( Float64 ,( n +1) )
7 trayMB [1]=0
8 times =[0: Delta : T ]
9
10 for i =1: n
11 trayMB [ i +1]= trayMB [ i ]+ rand ( Normal (0 ,1) ) * Delta
12 end
13
14 plot ( times , trayMB )
15 ylabel ( " W " ) ;
16 xlabel ( " Tiempo " )
17
18 return ( trayMB )
19
20 end

En la Figura ?? podemos ver una trayectoria de un movimiento Browniano con T = 1 y n = 1, 000.

5.1.2. Aproximando al movimiento Browniano por una caminata aleatoria

En la Figura ?? podemos ver una trayectoria de un caminata aleatoria con T = 1, y cuando n = 10,
n = 100 y n = 1, 000. Podemos ver la convergencia conforme n crece.

Consideremos una sucesión {Xi }ni=1 de variables aleatorias independientes tales que

Pr[Xi = 1] = Pr[Xi = −1] = 0.5.

 5.1. MOVIMIENTO BROWNIANO O PROCESO DE WIENER 67

Si definimos
n
X
Sn = Xi ,
i=1

entonces si n crece a ∞ tenemos por el Teorema del Lı́mite Central tenemos que

S[nt]
Pr[ √ < x] −→ Pr[Wt < x],
n

es decir
S d
√n −
→ N(0, 1).
n

La aproximación de la Figura ?? fue generada usando el Código 5.2 para diferentes valores de n.

Código 5.2.
1
2 # ##################################
3 # Approx MOVIMIENTO BROWNNIANO CAMINATA ALEATORIA
4 # ##################################
5 AMBCA = function (n , T )
6
7 Delta = T / n
8

9 E =[ -1 ,1]
10 p =[0.5 ,0.5]
11 p = WeightVec ( p )
12 S = cumsum ( sample (E , n ) )
13 times =[0: Delta : T ]
14

15
16
17 labels = n * times
18 tl = length ( labels ) - count ( c -> ( 0 < c ) , labels )
19
20

21 W = zeros ( length ( labels ) )

22 W [ tl +1: end ]= S [ labels [ tl +1]: labels [ end ]]
23 W = W / sqrt ( n )
24 plot ( times , W )
25 ylabel ( " CA " ) ;
26 xlabel ( " Tiempo " )
27
28 return ( W )
29
30
31 end

Nota 5.1. Es importante notar que es remondable usar el Algoritmo 5.1 cuando sólo se necesita simular
una posición especı́fica de un movimiento Browniano, por ejemplo, en la evaluación del pauoff de un
subyacente en una opción Europea. Por otro lado, si se requiere conocer información sobre la evolución
de la trayectoria del proceso, se pueden utilizar cualquiera de las formas de simulación presentadas, por
ejemplo, en la valuación de opciones Americanas, Asiáticas y exóticas.
 68 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

5.1.3. Puente Browniano

Definición 5.2. Se define como puente Browniano al movimiento Browniano {Wtp , t ∈ [t0 , T ] : Wt0 =
x, WT = y} que empieza en x al tiempo t0 y termina en y al tiempo T y está dado por
t − t0
Wtp = x − Wt−t0 − (WT −t0 − y + x).
T − t0

El Código 5.3 es una implementación para simular trayectorias de un puente Browniano usando la Defi-
nición 5.2.

Código 5.3.
1 # ###############################
2 # Puente Browniano de t0 a T
3 # ###############################
4 BBridge = function (x ,y ,n , t0 , T )
5
6 Delta =( T - t0 ) / n
7 t0 = t0 / n
8

9 bridge = zeros ( n +1)

10
11
12 W = BM (n , Delta )
13
14 for i =1: n +1
15 bridge [ i ]= x + W [ i ] -(( i -1) / n ) *( W [ n +1] - y + x )
16 end
17
18 times =[0: Delta : T ]
19 plot ( times , bridge )
20 ylabel ( " Puente Browniano " ) ;
21 xlabel ( " Tiempo " )
22
23 return ( bridge )
24
25

26 end

En la Figura 5.1 presentamos algunas trayectorias de un puente Browniano {Wtp , t ∈ [0, 1] : Wt0 = 0, WT =
0}.

5.2. Ecuaciones diferenciales estocásticas

Es común modelar la evolución en el tiempo de un sistema por ecuaciones diferenciales, las cuales usual-
mente describen la razón de cambio del sistema. Si denotamos por Xt el estado del sistema al tiempo t
para t ≥ 0, podemos caracterizar la evolución del sistema en el tiempo por la siguiente ecuación,
d(Xt )
= a(t)Xt , (5.1)
dt
5.2. ECUACIONES DIFERENCIALES ESTOCÁSTICAS 69

Figura 5.1: Trayectorias de puente Browniano de 0 a 0

donde X0 = x0 es el estado inicial del sistema y a(t) representa la función de cambio al tiempo t. Este
modelo es adecuado cuando la función de cambio es completamente conocida. En muchos casos esta
función no es completamente conocida, pero es razonable asumir que está sujeta a algún tipo de efecto
aleatorio, entonces tenemos que

a(t) = b(t) + ruido, (5.2)

donde no es totalmente conocido el comportamiento del término ruido, pero sabemos que está regido por
una distribución de probabilidad. Entonces una forma (heurı́stica) más general de escribir a la ecuación
(5.1) es:
d(Xt )
= b(t, Xt ) + σ(t, Xt ) × ruido (5.3)
dt

donde b, σ son funciones reales conocidas. Para modelar el término de ruido se utilizan los incrementos
de un proceso de Wiener. Si la evolución del sistema es observado en el intervalo de tiempo [0, T ] y
consideramos una partición de este intervalo, 0 = t0 < t1 < . . . < tn = T , podemos escribir a la ecuación
(5.3) en forma de incrementos como sigue

Xi+1 − Xi = b(ti , Xi )∆ti + σ(ti , Xi )∆Wi , (5.4)

70 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

donde Xi = Xti , ∆Wi = Wti+1 −Wti y ∆ti = ti+1 −ti . De esta forma, utilizando la ecuación (5.3) podemos
describir la evolución del proceso X = {Xt }Tt=0 por:

n−1
X n−1
X
Xt = X0 + b(ti , Xi )(∆ti ) + σ(ti , Xi ) × (∆Wi ). (5.5)
i=1 i=1

Dado que esta ecuación depende de la partición del intervalo, el modelo resulta ser más adecuado en el
lı́mite ∆ti → 0. Esto nos lleva a la siguiente ecuación
Z t Z t
Xt = X0 + b(ti , Xi )ds + σ(ti , Xi )dWs , (5.6)
0 0

donde el último término es conocido como integral estocástica. Podemos escribir la ecuación anterior en
forma diferencial
dXt = b(t, Xt )dt + σ(t, Xt )dWt (5.7)
para toda t ≥ 0 y X0 = x0 , a ecuaciones de este tipo se les conoce como ecuaciones diferenciales estocásti-
cas.

La ecuación (5.6) involucra dos tipos de integrales:

Z t
f (ti , Xi )ds,
0

que es una integral de Riemann y

Z t
f (ti , Xi )dWs
0

llamada integral de Itô o integral estocástica respecto al proceso de Wiener y puede ser interpretada
de la siguiente forma.

Definición 5.3. Sea f (t, x) una función continua en t y x, tal que

Z T
E(f 2 (t, Xt ))dt < ∞,
0

entonces la integral de Itô de f (t, Xt ) con respecto a Wt está definida como

Z T n−1
X
f (ti , Xi )dWt = lı́m f (ti , Xi )∆Wi , (5.8)
0 n→∞
i=0

donde 0 = t0 < t1 < . . . < tn = T entendiendo el lı́mite en el sentido que la distancia máxima entre dos
puntos de la partición tiende a cero.

Para que la definición anterior tenga sentido se necesita que el proceso X sea adaptado a la filtración
natural del proceso de Wiener.

Proposición 5.1. Propiedades de la integral de Itô Si X es Itô integrable entonces:

RT
1. E[ 0 Xs dWs ] = 0.
RT RT
2. V ar[ 0 Xs dWs ] = 0 E[Xs2 ]dWs , isometrı́a de Itô.
5.3. PROCESOS DE DIFUSIÓN 71

3. Si X y Y son Itô integrables y a, b son dos constantes, entonces

Z T Z T Z T
(aXs + bYs )dWs = a Xs dWs + b Ys dWs ,
0 0 0

linealidad.
RT
4. 0 aWs dWs = aW (T ).
RT WT2 T
5. 0 Ws dWs = 2 − 2.
RT
6. El proceso Mt = M0 + 0 Xs dWs es una martingala.

Una prueba de estas propiedades se puede consultar en [15].

5.3. Procesos de difusión

Definición 5.4. Un proceso de difusión X = {Xt }t≥0 es la solución a una ecuación de la forma
(5.6) o (5.7), donde los coeficientes b(t, x) y σ(t, x) son funciones reales de t y x, y se les conoce como
coeficiente de deriva y difusión respectivamente. Denotemos por I = [r, s] al espacio de estados.

En la definición anterior b(t, x) y σ(t, x) pueden ser interpretados como la media y varianza infinitesimal
de los incrementos del proceso X, es decir, si ∆hXt = Xt+h − Xt y x ∈ I, entonces
1
limh↓0 E[∆hXt |Xt = x] = b(t, x)
h
y
1
limh↓0 P[(∆hXt )2 |Xt = x] = σ(t, x)
h
Al proceso de difusión se le interpreta como el estado del sistema que evoluciona de manera determinista
gobernado por la parte no aleatoria de la ecuación pero perturbado por un ruido aditivo dado por la
integral estocástica.

Para que una ecuación de la forma (5.7) tenga alguna solución se deben imponer condiciones a sus
coeficientes. Existen resultados de existencia y unicidad para ecuaciones diferenciales estocásticas que
establecen condiciones de regularidad para sus coeficientes. El siguiente resultado es básico para dichas
condiciones.
Teorema 5.1. Si los coeficientes b(t, x) y σ(t, x) de la ecuación (5.7) satisfacen la condición de Lipschitz
para toda x ∈ R,
|b(t, x) − b(t, y)|2 + |σ(t, x) − σ(t, y)|2 ≤ K|x − y|2 ,
y la condición de crecimiento en x
|b(t, x)|2 + |σ(t, x)|2 ≤ K(1 + |x|2 ),
con K > 0 una constante, entonces existe un proceso de difusión X = {Xt }t≥0 solución de (5.7) que es
adaptado, continuo, uniformemente acotado en L2 y es único en el sentido de indistinguibilidad 1 .

A esta solución se le conoce como solución fuerte. Una prueba de este resultado se puede encontrar en
[15].
1
{Xt } y {Yt } son indistinguibles si P (Xt = Yt ) = 1 para cada t ≥ 0
72 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

5.3.1. Proceso de difusión ergódicos

Si un proceso de difusión es ergódico entonces, para cualquier función medible f , se tiene que con proba-
bilidad 1:
1 T
Z Z
f (Xs )ds → f (x)π(x)dx = E[h(U )],
T 0 R

donde π es la densidad estacionaria del proceso de difusión y U es una varibale aleatoria con densidad π.

Nota 5.2. Si un proceso de difusión X tiene distribución estacionaria, ésta está dada por

m(x)
π(x) = , (5.9)
M

donde
1
m(x) = , (5.10)
σ 2 (x)s(x)
es la medida de velocidad,  Z x 
b(y)
s(x) = exp −2 dy , (5.11)
x0 σ 2 (y)
es la medida de escala y M la constante de normailzación.

5.3.2. Ecuaciones de Kolmogorov

Dado que los procesos de difusión tienen la propiedad de Markov, tiene sentido definir las probabilidades
de transición del estado x al tiempo s al estado y en un tiempo posterior t (0 ≤ s < t). Denotamos esas
propabilidades por
p(t − s, y|x) := P[Xt = y|Xs = x].
Las probabilidades de transición satisfacen las ecuaciones de Kolmogorov (ver [15]), forward,

∂p(t − s, y|x) ∂p(t − s, y|x) 1 ∂ 2 p(t − s, y|x)

= −b(x) + σ(x) (5.12)
∂t ∂x 2 ∂x2
y backward
∂p(t − s, y|x) ∂p(t − s, y|x) 1 ∂ 2 p(t − s, y|x)
= −b(x) − σ(x) . (5.13)
∂s ∂x 2 ∂x2
Una relación muy importante para la simulación de trayectorias de procesos de difusión y a su vez para
hacer inferencia, es obtenida de la ecuación 5.12 cuando t → −∞, entonces tenemos que

d2 2 d
2
(σ (x)π(x)) = 2 (b(x)π(x)).
dx dx
De esta ecuación es fácil obtener
1 d
b(x) = 2
(2π(x)) (σ 2 (x)π(x)),
dx dx
e integrando tenemos que Z x
2 2
σ (x) = b(y)π(y)dy.
π(x) 0
5.4. FÓRMULA DE ITÔ 73

Generador infinitesimal

Si el proceso X es solución a la ecuación (5.5) con X0 = x, al operador diferencial L definido por

σ 2 (x) 00 0
(L f )(x) := f (x) + b(x)f (x) (5.14)
2
para f dos veces diferenciable, es llamado el generador infinitesimal del proceso de difusión X.

5.3.3. Puentes de difusión

Definición 5.5. A una solución de la ecuación 6.2 en el intervalo [t1 , t2 ] (0 ≤ t1 < t2 ≤ ∞) tal que
Xt1 = a y Xt2 = b, con a, b ∈ I, se le llama puente de difusión puente (t1 , a, t2 , b).
Ejemplo 5.1. Consideremos el proceso Ornstein-Uhlenbeck que es solución a la ecuación

dXt = −αXt dt + σdWt ,

condicionado a X0 = a y X1 = b para a, b ∈ R. Desde que las densidades de transición de este proceso

siguen una ley Gaussiana es posible calcular las densidades de transión del puente de difusión (0, a, 1, b)
(ver [1]).

Por otro lado, en [3] proponen un algoritmo para generar el puente de difusión (0, a, 1, b) de forma exacta
pero que resulta numericamente más estable. El algortimo se presenta en el siguiente lema.
Lema 5.1. Consideremos los momentos en el tiempo 0 = t0 < t1 ≤< tn < tn + 1 = 1 y sea X0 = a. Si

Xti = e−α(ti −ti−1 ) Xti−1 + Wi

para i = 1, . . . , n y Wi ∼ N (0, σ 2 (1 − e−2α(ti −ti−1 ) /(2α)) y definimos

eαti − e−αti
Zti = Xti + (b − Xtn+1 )
eαtn+1 − e−αtn+1
para i = 0, . . . , n + 1. Entonces {Zti }n+1
i=0 tiene la misma dsitribució que un puente de difusión Ornstein-
Uhlenbeck (0, a, 1, b).

5.4. Fórmula de Itô

Otro resultado importante en teorı́a de procesos de difusión es la fórmula de Itô, que ayuda a encontrar
soluciones a ecuaciones diferenciales estocásticas.
Teorema 5.2. (Fórmula de Itô) Si X = {Xt }t≥0 es un proceso de difusión dado por la ecuación (5.7)
y f (t, x) una función de clase C 1 como función de t y de clase C 2 como función de x, entonces el proceso
Yt = f (t, Xt ) es también un proceso de difusión y es solución a la ecuación

1
dYt = ft (t, Xt )dt + fx (t, Xt )dXt + fxx (t, Xt )(dXt )2 . (5.15)
2

En la ecuación (5.15) los subı́ndices denotan las derivadas, los detalles de las demostraciones se pueden
ver en [15].
 74 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

Ejemplo 5.2. Movimiento Browniano Geométrico o modelo Black-Scholes-Merton. Sea

St = S0 exp (r − σ 2 /2)t + σWt

(5.16)

y si hacemos
f (t, x) = S0 exp (r − σ 2 /2)t + σx ,

entonces f (t, Wt ) = St y por lo tanto

dSt = rSt + σSt dWt .

Este modelo fue utilizado por Black-Scholes-Merton en el contexto financiero para modelar precios de
acciones (ver [13]) donde r es intepretado como la tasa de interés y σ como la volatilidad de los activos
riesgosos. Si S0 = s0 y hacemos
η = log(s0 ) + (r − σ 2 /2)t
y
ν = σ 2 t,
de 5.16 tenemos que para este proceso las densidades de tarnsición tienen distribución log-normal de
parámetros η y ν, es decir,

(log(y) − η)2
 
1
p(t − s, y|s0 ) = √ exp .
yν 2π 2ν 2

El Código 5.4 simula trayectorias de una difusión determinada por la ecuación 5.2 basado en las proba-
bilidades de transición obtenidas.
Código 5.4.
1 # ####################################
2 # BG - BSM
3 # ####################################
4 BGBSM = function (x ,r , sigma ,n , T )
5
6 Delta = T / n
7

8 path = zeros ( Float32 , n +1)

9 path [1]= x
10
11 for i =2: n +1
12 etha = log ( path [i -1]) +( r -( sigma ^2) /2) * Delta
13 nu =( sigma ^2) * Delta
14 path [ i ]= rand ( LogNormal ( etha , nu ) )
15 end
16
17 times =[0: Delta : T ]
18 plot ( times , path )
19 ylabel ( " BG - BSM " ) ;
20 xlabel ( " Tiempo " )
21
22 return ( path )
23
24

25 end
 5.4. FÓRMULA DE ITÔ 75

La Figura ?? se presentan 3 trayectorias de un movimiento Browniano Geométrico con parámetros r = 0.4

(azul), r = 0.5 (verde) y r = 0.6 (rojo), todos con σ = 1.
Ejemplo 5.3. Proceso Ornstein-Uhlenbeck o proceso de Vasicek. Un importante proceso de difu-
sión X = {Xt }t≥0 que será de gran utilidad en el desarrollo de este trabajo es la solución a la ecuación
diferencial estocástica
dXt = (λ − αXt )dt + σdWt , (5.17)
donde α, σ > 0 son constantes y X0 = x. Utilizando la fórmula de Itô se llega a la siguiente solución
explicita de esta ecuación:
  Z t
λ −αt λ
Xt = + e x− +σ eα(s−t) dWs . (5.18)
α α 0

De la solución tenemos que las probabildides de transición son Gausianas de parámetros

 
λ −αt λ
η = +e x−
α α
y
σ 2 (1 − e−2αt )
ν= .
2α

Con base a las densidades de transición podemos simular trayectorias de este proceso de difusió. Con el
Código 5.5 de podemos generar trayactorias.
Código 5.5.
1 # ######################################
2 # OU - Vasicek
3 # ######################################
4 Vasicek = function (x , alpha , lambda , sigma ,n , T )
5

6 Delta = T / n
7
8 vas = zeros ( Float32 , n +1)
9 vas [1]= x
10
11

12 for i =2: n +1
13 vas [ i ]= lambda / alpha + exp ( - alpha * Delta ) *( vas [i -1] - lambda / alpha ) + sigma ^2*(1 -
exp ( -2* alpha * Delta ) ) /(2* alpha ) * rand ( Normal () )
14
15 end
16

17 times =[0: Delta : T ]

18 plot ( times , vas )
19 ylabel ( " OU - Vasicek " ) ;
20 xlabel ( " Tiempo " )
21
22 return ( vas )
23
24 end

En la Figura ?? se presenta la simulación de la solución a la ecuación 5.17 con λ = α = 1 y tres diferentes

valores para σ = 0.5 (rojo), σ = 1.0 (azul) y σ = 2.0 (verde).
76 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)

λ σ2
Nota 5.3. Este proceso es ergódico y tiene distribución invariante Gausiana de parámetros α y 2α .

Ejemplo 5.4. Cox-Ingersoll-Ross. Este proceso de difusión es usado en el ambiente financiero para
modelar tasas de interés a corto plazo (ver [6].) y es la solución a la ecuación diferencias estocástica
p
dXt = (λ − αXt )dt + σ Xt dWt , (5.19)

donde X0 = x0 > 0 y α, λ, σ > 0. Si 2λ > sigma2 el proceso es estacionario.

La solución a la ecuación 5.19 es:

  Z t
λ
e−αt + σe−αt
p
Xt = X0 − eαs Xs dWs .
α 0

Difusiones en sistemas dinámicos.

Una de las principales caracteristicas estudiadas en sistemas dinámicos es la estacionalidad de los procesos
que describen el equilibrio del sistema. Aquı́ enunciaremos los procesos de disfusión más populares es esta
área.

Ejemplo 5.5. Modelo tipo N . Este proceso es solución a la siguiente ecuación diferencial estocástica
√
dXt = r(θ − Xt )dt + dWt ,

con r,  > 0.

Ejemplo 5.6. Modelo tipo N . Este proceso es solución a la siguiente ecuación diferencial estocástica
p
dXt = r(θ − Xt )dt + Xt (1 − Xt )dWt ,

con r,  > 0.

Ejemplo 5.7. Modelo tipo N . Este proceso es solución a la siguiente ecuación diferencial estocástica
√
dXt = r(θ − Xt )dt + dWt ,

con r,  > 0.

Difusiones en epidemiologı́a

Para tener un panorama general de la teorı́a matemática en modelos de epidemiologı́a se puede consultar
[2]. Supongamos que Xt representa la fracción de población infectada al tiempo t, en general, la tasa de
cambio es modelada con el proceso de difusión que es la solución a:
p
dXt = (aXt (1 − Xt ) − bXt + c(1 − Xt ))dt + Xt (1 − Xt )dWt , (5.20)

donde a > 0 es la tasa de transmisión (persona a persona), b > 0 representa la tasa de recuperación y
c > 0 es la tasa de transmisión desde una fuente externa y  > 0.

Los modelos con Arno

5.5. EJERCICIOS 77

5.4.1. Transformación de Lamperti

Una aplicación de la fórmula de Itô muy importante en las áreas de simulación y estimación de parámetros
para escuaciones diferenciales estocásticas es transformación estándar conocida como transformada de
Lamperti dada por transformation Z x
1
h(x) = dy, (5.21)
x ∗ σ(y)
donde x∗ es un elemento arbitrario en el espacio de estados de X. Entonces, por la fórmula de Itô

Yt = hβ (Xt )

satisface la ecuación diferencial estocástica

dYt = µ(Yt )dt + dWt , (5.22)

donde
b(h−1 (y)) 1
µ(y) = − σ 0 (h−1 (y)),
σ(h−1 (y)) 2
donde σ 0 (x) denota la derivada de σ con respecto a x.

5.4.2. Cociente de verosimilitud para procesos de difusión

En el contexto de inferencia estadı́stica para procesos de difusión el Teorema de Girsanov es usado

para obtener el cociente de verosimilitud para hacer inferencia por máxima verosı́militud.
Teorema 5.3. Teorema de Girsanov. Sea P la medida de probabilidad inducida en ([0, T ], BR ([0, T ]))
por Y , la solución de la ecuación (5.22) donde µ satisaface las condicones de la fórmula de Itô y sea Q la
medida de Wiener, entonces las medidas son equivalentes y su correspondiente derivada de Radon-Nikodı́n
es Z T
1 T 2
Z 
dP
= exp µ(Ys )dYs − µ (Ys )ds
dQ 0 2 0

Para una prueba de este resultado ver [11].

Sea S

5.5. Ejercicios
1. Encontrar las probabilidades de transición del CIR e implementar un algoritmo para simular tra-
yectorias de este proceso de difusión.

2. Encontrar solución, la distribución estacionaria, las densidades de transición, la media, la moda y

la varianza de la difusión del Ejemplo 5.5.

3. Encontrar solución, la distribución estacionaria, las densidades de transición, la media, la moda y

la varianza de la difusión del Ejemplo 5.6.

4. Encontrar solución, la distribución estacionaria, las densidades de transición, la media, la moda y

la varianza de la difusión del Ejemplo 5.7.
78 CAPÍTULO 5. ECUACIONES DIFERENCIALES ESTOCÁSTICAS (EDE)
Capı́tulo 6

Inferencia estadı́stica para ecuaciones

diferenciales estocásticas

6.1. Simulación de procesos de difusión

En esta sección presentaremos los dos métodos más importantes para la simulación de trayectorias a
soluciones de ecuaciones diferenciales estocásticas. Estos métodos son utilizados cuando no se conoce
explicitamente las probabilidades de transición o aún cuando se conocen resulta difı́cil de simular. Debido
a que estos métodos están basados en aproximaciones a tiempo discreto a los procesos en tiempo continuo
es importante empezar definiendo el tipo de convergencia de los procesos discretos a los continuos.

Definición 6.1. Convergencia Fuerte. Se dice que una aproximación a tiempo discreto X δ de un
proceso a tiempo continuo X, donde δ denota el incremento de tiempo máximo en la discretización, se
dice que es de nivel fuerte de convergencia τ si

E[|Xtδ − Xt |] ≤ Cδ τ

para todo momento t ≥ 0 , δ < δ0 con δ0 > 0 y C una constante independiente de δ.

Definición 6.2. Convergencia Débil. Considerando otra vez la aproximación a tiempo discreto X δ de
un proceso a tiempo continuo X, se dice que es de nivel débil de convergencia τ si

|E[g(Xtδ )] − E[g(Xt )]| ≤ Cδ τ

para todo t ≥ 0, δ < δ0 con δ0 > 0, C una constante independiente de δ, cualquier función polinómica
creciente continua y diferenciable de orden 2(τ + 1).

6.1.1. Método de Euler

Consideremos un proceso de difusión X = {Xt }t≥0 que es solución a la ecuación diferencial estocástica

dXt = b(Xt , t)dt + σ(Xt , t)dWt

79
 80CAPÍTULO 6. INFERENCIA ESTADÍSTICA PARA ECUACIONES DIFERENCIALES ESTOCÁSTICAS

con valor incial determinista Xt0 = x0 y un nivel de discretización

0 = t0 < t1 < ... < tn = T

en el intervalo [0, T ]. La aproximación de Euler de X es un proceso estocástico continuo Y que satisface

el siguiente sistema iterativo:

Yti+1 = Yti + b(Yti , ti )∆i + σ(Yti , ti )∆Wi

para i = 0, 1, . . . , n − 1, con Y0 = x0 y
∆i = ti+1 − ti ,
∆Wi = Wti+1 − Wti

entonces ∆Wi ∼ N (0, ∆i ).

Por otro lado como entre cualesquiera dos momentos de observación ti y ti+1 , el proceso puede ser definido
diferenciable. Una propuesta natural es considererar una interpolación lineal para Yt definida por

t − ti
Yt = Yti + (Yti+1 − Yti )
ti+1

para t ∈ [ti , ti+1 ).

Código 6.1.
1 #
######################################################################################

2 # Drawn Paths of diffusion processes by Euler Scheme

3 #
######################################################################################

4
5 function Euler ( theta , init , delta , n )
6
7 brow_inc = rand ( Normal (0 , sqrt ( delta ) ) ,n )
8 tray = zeros ( Float64 ,( n +1) )
9 tray [1]= init
10
11 for i =1: n
12 tray [ i +1]= tray [ i ]+ b ( tray [ i ] , theta ) * delta + sigma ( tray [ i ] , theta ) * brow_inc [
i]
13 end
14

15 tray
16
17 end

Ejemplo 6.1. Aquı́, se presenta una simulación de un proceso de difusión con el esquema de Euler,
consideramos el proceso de Ornstein-Uhlenbeck, que es una solución de la ecuación diferencial estocástica

dXt = −αXt dt + σdWt ,

 6.1. SIMULACIÓN DE PROCESOS DE DIFUSIÓN 81

donde α > 0 y σ > 0 son los parámetros del proceso W es un proceso de Wiener estándar. Si suponemos
que X0 = 0 y consideremos un intervalo de realización [0, 10] con un nivel de discretización 0 = t0 < t1 <
... < tn = 10, donde n = 100 y ∆ = ∆i = 0.1, para todo i = 1, 2, ..., n − 1, entonces el método de Euler
para este proceso es:

Yti+1 = Yti − αYti ∆ + σ∆W.

En la Figura ?? podemos ver trayectorias caracterı́sticas de un Ornstein-Uhlenbeck con diferentes paráme-

tros α = 1, σ = 1 (azul), α = 1.5, σ = 0.3 (rojo) y α = 3.5, σ = 0.3 (verde). Simular de las transiciones y
comparar los métodos qqplots

6.1.2. Método de Milstein

En la aproximación de Milstein se utiliza el lema de Itô para aumentar la presición en la aproximación

mediante la incorporación del término de segundo orden, podemos escribir a esta aproximación por:

1
Yti+1 = Yti + b(Yti , ti )∆i + σ(Yti , ti )∆Wi + σ(Yti , ti )σ 0 (Yti , ti )ti [(∆Wi )2 − ∆i ],
2
Código 6.2.
1 #
######################################################################################

2 # Drawn Paths of diffusion processes by Milstein Scheme

3 #
######################################################################################

4
5 function Milstein ( theta , init , delta , n )
6
7 brow_inc = rand ( Normal (0 , sqrt ( delta ) ) ,n )
8 tray = zeros ( Float64 ,( n +1) )
9 tray [1]= init
10
11 for i =1: n
12 tray [ i +1]= tray [ i ]+ b ( tray [ i ] , theta ) * delta + sigma ( tray [ i ] , theta ) * brow_inc [
i ]+(1/2) * sigma ( tray [ i ] , theta ) * sigma_x ( tray [ i ] , theta ) *( brow_inc [ i ]^2 -
delta )
13 end
14
15 tray
16
17 end

Ejemplo 6.2. Para ejemplificar el uso de este método de simulación cosideremos el movimiento Brow-
niano geométrico, el cual es solución a la ecuación diferencial estocástica

dXt = θ1 Xt dt + θ2 Xt dWt ,

para este proceso, tenemos que b(x, t) = θ1 x, σ(x, t) = θ2 x y σx (x, t) = θ2 entonces el esquema de Milstein
para este proceso es
82CAPÍTULO 6. INFERENCIA ESTADÍSTICA PARA ECUACIONES DIFERENCIALES ESTOCÁSTICAS

1
Yti+1 = Yti + θ1 Yti ∆i + θ2 ∆Yti Wi + θ2 Yti [(∆Wi )2 − ∆i ]
2
θ 2 θ2
= Yti (1 + ∆θ1 − 2 ) + θ2 ∆Wi + 2 Yti (∆Wi )2
2 2

En la Figura ?? podemos ver trayectorias caracterı́sticas de un Movimiento Browniano Geométricocon

diferentes parámetros θ1 = 1, θ2 = 0.1 (azul), θ1 = 1, θ2 = 1 (rojo) y θ1 = 0.1, θ2 = 0.1 (verde).

Simular de las transiciones y comparar los métodos qqplots

6.1.3. Simulación de puentes de difusión

En esta sección presentaremos un algoritmo para simular puentes de difusión presentado en [3].

Sea X un proceso de difusión solución de 6.2, con espacio de estados I = (r, s) y que los coeficientes
satisfacen las condiciones del teorema 5.1. Sea a, b ∈ I y W 1 y W 2 dos procesos de Wiener independientes
y X 1 y X 2 soluciones de
dXti = b(Xti )dt + σ(Xti )dWti (6.1)
i = 1, 2, X01 = a y X02 = b. En el siguiente teorema utilizaremos estos dos procesos para aproximar el
puente (0, a, ∆, b)
Teorema 6.1. Sea τ = inf{0 ≤ t ≤ ∆|Xt1 = X∆−t
2 }y

Xt1 si 0 ≤ t ≤ τ


Zt =
 2
X∆−t si τ < t ≤ ∆
Entonces la distribución de Z = {Zt }0≤t≤∆ , condicionado a {τ ≤ ∆}, es igual a la distribución de un
puente (0, a, ∆, b), condicionado al evento que el puente es alcanzado por una difusión (independiente)
con EDE (6.2) y distribución inicial con densidad p∆ (b, ∗).

La prueba de este resultado se puede consultar a detalle en [3].

OU con este algoritmo y comparar con el exacto

6.1.4. Estimadores Monte-Carloy técnicas de reducción de varianza

Usando los métodos de Monte-Carlo se puede estimar algunas cosas, poner ejemplos financieros (precios
de opciones). Revisar las técnicas de Malliavin en [14] y el método presentado en [12].

6.2. Estimación por máxima verosimı́litud

Abordaremos el problema de estimación de parámetros en procesos de difusión cuando éstos han sido
observados a tiempo discreto. Consideremos que las difusiones son observadas en el intervalo de tiempo
[0, T ] en n momentos 0 = t1 < t2 < . . . < tn = T y el tiempo entre cada observación es igual entre
cualesquiera dos observaciones consecutivas, es decir, ∆ = ti − ti−1 para i = 2, . . . , n con los siguientes
tipos de información:
6.2. ESTIMACIÓN POR MÁXIMA VEROSIMÍLITUD 83

1. Muestra grande: Se mantiene fijo a ∆ = ∆n y se aumenta el número de observaciones, es decir,

n crece. En este caso, aumenta el intervalo de observación [0, T ]. Para este tipo de información se
asumirá que los procesos de difusión son estacionarios y/o ergódicos.

2. Alta frecuencia: Se hace tender a cero a ∆ conforme aumenta n y el intervalo de observación [0, T ]
se mantiene fijo.

3. Crecimiento rápido: ∆ = ∆n decrece a cero cuando n crece y además el intervalo de observación

[0, T ] crece como función de n, es decir, n∆n → ∞. Para este caso tembién se asumira que los
procesos de difusión son estacionarios y/o ergódicos.

6.2.1. Modelo Continuo

Consideremos un proceso de difusión X = {Xt }t≥0 que es solución a la ecuación diferencial estocástica

dXt = b(Xt , θ)dt + σ(Xt , θ)dWt , (6.2)

donde θ ∈ Θ ⊂ R es un parámetro multidimensional y θ0 es el verdadero valor del parámetro. Las funciones

b:R×Θ→R

y
σ : R × Θ → (0, ∞)

son conocidas y cumplen con condiciones para la existencia de una solución de (6.2) (ver Teorema 5.1).
Denotamos por I = (l, r) con −∞ ≤ l < r < ∞ al espacio de estados y supongamos que es el mismo para
toda θ. Supongamos la solución X de (6.2) tiene valor inicial determinista Xt0 = x0 y que se cumplen las
condiciones para la exitencia de la distribución estacionaria πθ .

Cuando se tiene el proceso observado a tiempo continuo, la estimación del parámetro resulta eficiente. En
particular, desde que la variación cuadrática de un proceso de difusión que es solución de (6.2) está dada
por (ver [15])
Z T
< X, X >T = σ 2 (Xs , θ)ds, (6.3)
0

entonces,
2n
Z T X
2
< X, X >T = σ (Xs , θ)ds = lı́m (XT ∧k/2n − XT ∧(k−1)/2n )2
0 n→∞
k=1

en probabilidad bajo la medida Pθ (La distribución de X). Si asumimos que σ(x, θ) = σ(x) y que el resto
de los parámetros están sólo en la deriva, entonces, éstos pueden ser estimados por máxima verosimilitud.

Para el intervalo de tiempo [0, T ] la función de verosimilitud está dada por

T T
µ2 (Xs , θ)
Z Z 
µ(Xs , θ) 1
LT (θ) = exp dXs − ds . (6.4)
0 σ 2 (Xs ) 2 0 σ 2 (Xs )
escribir como llevar a esta

Nota 6.1. Podemos afirmar que esta situación ocurre con probabilidad cero.
84CAPÍTULO 6. INFERENCIA ESTADÍSTICA PARA ECUACIONES DIFERENCIALES ESTOCÁSTICAS

6.3. Caso discreto

Suponemos que tenemos el proceso observado a tiempo discreto en el intervalo [0, T ] a los tiempos {0 =
t1 < t2 < . . . < tn = T }, es decir, la información que se tiene es Xobs = {Xt1 = x1 , . . . , Xtn = xn }.
Denotamos por Fn = σ{xi ; 1 ≤ i ≤ n} a la σ-álgebra generada por las n observaciones y por F0 a la
σ-álgebra trivial. Si denotamos a las densidades de transición (entre las observaciones) por pθ (∆, xi |xi−1 ),
por la propiedad de Markov de los procesos de difusión podemos escribir la función de verosimilitud
basada en Xobs como
n
Y
Ln (θ) = pθ (∆, xi |xi−1 )pθ (x0 ) (6.5)
i=1

y la log-verosimilitud es
n
X
`n (θ) = log pθ (∆, xi |xi−1 ) + log pθ (x0 ).
i=1

Nota 6.2. Bajo condiciones de regularidad las probabildides de transición pθ (∆, xi |xi−1 ) satisfacen las
ecuación de Kolmogorov 5.12 y 5.13 pero unicamente en casos muy sencillos se puede resolver estas
ecuaciones diferenciales parciales.

6.3.1. Inferencia por verosimilitud exacta

Consideremos el siguiente conjunto de hipótesis necesarias para asegurar suficiencia y consistencia asintóti-
ca de los estimadores máximo verosimiles que presentaremos, para ver los detalles ver [7] y [10].

1. Creciento lineal. Existe una constante C (independiente de θ) tal que, para todo x ∈ I

|µ(x, θ)| + |σ(x, θ)| ≤ (1 + |x|).

2. Lipschitz. Existe una constante C (independiente de θ) tal que, para todo x, y ∈ I

|µ(x, θ) − µ(y, θ)| + |σ(x, θ) − σ(y, θ)| ≤ (|x − y|).

3. Coeficiente de difusión positivo.

ı́nf σ 2 (x, θ) > 0.
x∈I

4. Momentos acotados Para todo k > 0, todos los momentos de orden k existen y

supt∈T Exp|Xt |k < ∞.

Nota 6.3. Si los parámetros de los coeficientes de deriva y difusión son separables, la tasa de convergencia
es más rápida en el parámetro de difusión y es posibe utilizar la propuesta presentanda en [16] para
√
encontrar√los estimadores. Generalmente la tasa de convergencia del coefiente de difusión es n y la de
la deriva n∆n (condicionado a n∆3n → 0).

Ejemplo 6.3. Ornstein-Uhlenbeck o Vasicek

6.4. EJERCICIOS 85

Large sample, ∆ = 1
n α=1 2.5 % - 97.5 % β=3 2.5 % - 97.5 % σ=2 2.5 % - 97.5 %
100 1.2016 (0.71783,2.1760) 3.5340 (2.05939,6.4328) 2.0797 (1.68729,2.8068)
200 0.99640 (0.69757,1.4247) 3.19439 (2.19813,4.6007) 2.05284 (1.779,2.4380)
500 1.0033 (0.80208,1.2556) 3.057849 (2.42023,3.848) 2.038970 (1.85789,2.2646)
1000 1.0111 (0.86386,1.1839) 3.1094 (2.63914,3.6573) 2.0143 (1.88388,2.1667)

Cuadro 6.1:
High-frequency, T = 100
∆ α=1 2.5 % - 97.5 % β=3 2.5 % - 97.5 % σ=2 2.5 % - 97.5 %
1 0.9326945 (0.55658,1.5430) 2.91223 (1.68594,4.8569) 1.99326 (1.4897,2.548636)
1/2 0.8664383 (0.56497,1.2216) 2.74615 (1.71730,3.9431) 2.01957 (1.7898,2.2935)
1/5 0.9653574 (0.67449,1.2742) 2.99912 (2.02426,4.0291) 1.97410 (1.8987,2.1169)
1/10 1.0138 (0.72397,1.3124) 2.7642 (1.88794,3.6641) 1.94155 (1.85554,2.0347)

Cuadro 6.2:

Consideremos la difusión X solucioón de 5.17 wcon X0 = x0 , λ, α ∈ R, σ ∈ R+ y es ergódico si α > 0.

En este caso, la densidad condicional pθ (∆, xi |xi−1 ) tiene distribución Gaussiana de parámetros
β −α∆ 2 1 − e−2α∆
 
β
Xt+∆ |Xt ∼ N + (xt − )e ,σ .
α α 2α
Entonces, tenmos una expresión explicita para (6.5) y podemos encontrar los estimadores máximo ve-
rosı́miles.

El proceso de Ornstein-Uhlenbeck, en la versión de dos parámetros, definido por la ecuación difencial

estocástica
dXt = −αXt dt + σdWt ,
tiene una forma explı́cita para los estimadores máximo verosimil dado por :
 Pn 
1 i=1 Xi−1 Xi
α̂ = − log P n 2
∆ i=1 Xi−1

definido siempre que ocurra ni=1 Xi−1 Xi > 0, y

P

n
2α̂ X
σ̂ 2 = −2∆α̂
(Xi − Xi−1 e−∆α̂ )2 .
n(1 − e )
i=1

6.4. Ejercicios
1. Mostrar que el método de Euler tiene orden fuerte de convergencia con τ = 1/2.
2. Mostrar que esta aproximación tiene convergencia débil y fuerte de orden uno.

High-frequency
∆ n α=1 2.5 % - 97.5 % β = 3 2.5 % - 97.5 % σ=2 2.5 % - 97.5 %
1/1000 1000 3.9832 (-1.2793, 9.2737) 8.0231 (-1.9707,18.0656) 2.0301 (1.9442, 2.1227)

Cuadro 6.3:
86CAPÍTULO 6. INFERENCIA ESTADÍSTICA PARA ECUACIONES DIFERENCIALES ESTOCÁSTICAS
Apéndice A

Simulación de variables aleatorias

Definición A.1. Sea X variable aleatoria con soporte en E = {x1 , x2 , . . . , xn } con x1 < x2 < . . . < xn y
f su función de densidad asociada
Se define

φ(y) = xmin{i;Pi f (xk )≥y} .

k=1

Se dice que φ es la función de cuantiles asociada a la variable aleatoria X.

Lema A.1. Considérese U una v.a. uniforme en (0, 1) y X una v.a. con soporte en E = {x1 , x2 , . . . , xn }
con x1 < x2 < . . . < xn , con función de distribución f y con función de cuantiles φ, entonces se tiene que
φ(U ) tiene la misma distribución que X.

Demostración.
i
X
P[φ(U ) ≤ xj ] = P[min{i; f (xk ) ≥ U } ≤ j].
k=1
j
X
= P[U ≤ f (xk )].
k=1
j
X
= f (xk ).
k=1

por lo tanto φ(U ) tiene distribución f .

87
88 APÉNDICE A. SIMULACIÓN DE VARIABLES ALEATORIAS
Apéndice B

Convergencia de variables aleatorias y el

Teorema central del lı́mite.

En este apéndice se presentan los modos de convergencia de vectores aleatorios, ası́ como también sus
respectivas relaciones y los principales resultados que se utilizan a lo largo del trabajo. Este apéndice está
basado principalmente en [? ], [? ] y [? ].
Convergencia casi segura
En algunas situaciones la convergencia puntual resulta ser una condición muy fuerte pues se pide la con-
vergencia de la sucesión evaluada en todos y cada uno de los elementos del espacio muestral. Se puede ser
menos estricto y pedir, por ejemplo, que la convergencia se verifique en todo el espacio Ω excepto en un
subconjunto de probabilidad cero.

Definición B.1. (Convergencia casi segura). La sucesión de variables aleatorias ξ1 , ξ2 , . . .converge casi
seguramente a la variable ξ, si P {ω ∈ Ω : lı́m ξn (ω) = ξ(ω)} = 1.
n−→∞

Es decir, en la convergencia casi segura se permite que para algunos valores de ω, la sucesión numérica
ξ1 (ω), ξ2 (ω), . . . pueda no converger, sin embargo el subconjunto de Ω en donde esto suceda debe tener
c.s.
probabilidad cero. Para indicar la convergencia casi segura se escribe ξn −−→ ξ , o bien lı́m ξn = ξ c.s. A
n−→∞
menudo se utiliza el término convergencia casi donde quiera, o bien convergencia casi siempre para de-
notar este tipo de convergencia. Observe que omitiendo el argumento ω, la condición para la convergencia
casi segura se escribe en la forma más corta: P[lı́mn−→∞ ξn = ξ] = 1 , o simplemente P[ξn −→ ξ] = 1. Es
posible demostrar que el conjunto {ω ∈ Ω : ξn (ω) −→ ξ(ω)} es medible de modo que tiene sentido aplicar
la probabilidad. Puede también demostrarse que bajo este tipo de convergencia, el lı́mite es único casi
seguramente, es decir, si ξn converge a ξ c.s. y también converge a η c.s., entonces ξ = η casi seguramente.

Convergencia en probabilidad
Un tipo de convergencia aún menos restrictiva que la convergencia casi segura es la convergencia en pro-
babilidad la cual se define a continuación.

Definición B.2. (Convergencia en probabilidad). La sucesión de variables aleatorias ξ1 , ξ2 , . . . converge

en probabilidad a ξ, si para cada  > 0,
lı́m P {ω ∈ Ω : |ξn (ω) − ξ(ω)| > } = 0.
n−→∞

89
90APÉNDICE B. CONVERGENCIA DE VARIABLES ALEATORIAS Y EL TEOREMA CENTRAL DEL LÍMITE.

Cuando hablamos de vectores solo cambiamos el valos absoluto por la norma euclidiana.

p
Para denotar la convergencia en probabilidad se escribe ξn →
− ξ, y omitiendo el argumento ω la condición
se escribe lı́mn−→∞ P[| ξn − ξ |> ] = 0. Nuevamente puede comprobarse que el lı́mite es único casi segu-
ramente.

Convergencia en media
En este tipo de convergencia se usa la esperanza para determinar la cercanı́a entre dos variables aleatorias.

Definición B.3. (CONVERGENCIA EN MEDIA). La sucesión de variables aleatorias integrables ξ1 , ξ2 , . . .

converge en media a la variable aleatoria integrable ξ si
lı́m E[| ξn − ξ |] = 0.
n−→∞

m L1
A este tipo de convergencia también se le llama convergencia en L1 y se le denota por ξn −→ ξ, o ξn −→ ξ.
A partir de la definición de convergencia en media es inmediato preguntarse si de allı́ se sigue la conver-
gencia de la sucesión de medias. La respuesta es afirmativa.

Convergencia en media cuadrática

Nuevamente usando el concepto de esperanza pero ahora aplicado al segundo momento se tiene la con-
vergencia en media cuadrática. Más adelante demostraremos que la convergencia en media cuadrática
implica la convergencia en media.

Definición B.4. (CONVERGENCIA EN MEDIA CUADRÁTICA ). La sucesión de variables aleatorias

ξ1 , ξ2 , . . . converge en media cuadrática a ξ, si

lı́m E[| ξn − ξ |2 ] = 0.
n−→∞

En este tipo de convergencia se presupone que tanto los elementos de la sucesión como el lı́mite mismo
son variables aleatorias con segundo momento finito. A este tipo de convergencia también se le llama
m.c. L2
convergencia en L2 , y se le denota por ξn −−→ ξ, o ξn −→ ξ.

En general puede definirse la convergencia en Lk , para cada entero k ≥ 1, cuando se cumple la condición
E[| ξn − ξ |k ] −→ 0. Resulta que mientras mayor es el valor de k, más restrictiva es la condición de
convergencia.

Convergencia en distribución
Este es el tipo de convergencia menos restrictiva de todas las mencionadas. En contextos más generales
se le llama también convergencia débil.

Definición B.5. (Convergencia en distribución). La sucesión de variables aleatorias ξ1 , ξ2 , . . . converge

en distribución a ξ, si para todo punto x en donde la función Fξ (x) es continua, se cumple que

lı́m Fξn (x) = Fξ (x).

n→∞
d ley
Cuando esto ocurre escribimos ξn −
→ ξ o bien ξn −−→ ξ.
 91

Resultados más importantes de los tipos de convergencia

Proposición B.1. Convergencia c.s. =⇒ convergencia en probabilidad.

Proposición B.2. Convergencia en Lp =⇒ convergencia en media.

Proposición B.3. Convergencia en media =⇒ convergencia en probabilidad.

Proposición B.4. Convergencia en probabilidad. =⇒ convergencia en distribución.

Todos los resultados anteriores de esta sección pueden ser consultados en [? ] en la parte de Convergencia
Capı́tulo 7.

Definición B.6. Sean ξ, ξ1 , ξ2 , . . . vectores aleatorios, entonces la sucesión {ξn }n∈N converge casi segura-
mente (en probabilidad) a ξ si la sucesión || ξ − ξn || converge c.s. (en probabilidad), en tal caso escribimos
c.s. p
ξn −−→ ξ (ξn →
− ξ).
p
Lema B.1. (Criterio de subsucesión)Sean ξ, ξ1 , ξ2 , . . . vectores aleatorios. Entonces ξn →
− ξ si y sólo si
cada subsuceción N 0 ⊂ N tiene una subsucesión N 00 ⊂ N 0 de tal manera que ξn → ξ c.s. através de la
p
subsucesión N 00 .En particular, ξn → ξ c.s. implica ξn →
− ξ.

Para su demostración ver [? ] página 40.

p
Proposición B.5. (aplicaciones continuas) Sean ξ, ξ1 , ξ2 , . . . variables aleatorias en Ω con ξn →
− ξ y
p
considera una función f : Ω → R de tal manera que f es c.s. continua. Entonces f (ξn ) → − f (ξ).

Demostración. Se deduce fácilmente del Lema B.1.

Proposición B.6. (Operaciones elementales) Sean ξ, ξ1 , ξ2 , . . . y η, η1 , η2 , . . . variables aleatorias con

p p p p p
ξn →
− ξ y ηn →
− η. Entonces aξn + bηn → − aξ + bη, para todo a, b ∈ R, y ξn ηn → − ξη. Además ξn /ηn → − ξ/η
cuando c.s. η 6= 0 y ηn 6= 0 para toda n.

Demostración. Se deduce del Lema B.1.

d c.s.
Proposición B.7. Sean ξ, ξ1 , ξ2 , . . . y η1 , η2 , . . . variables aleatorias con ξn −
→ ξ y ηn −−→ 0. Entonces
p
ξn ηn →
− 0.

Demostración. Sea  > 0 y δ > 0, sabemos que existe N ∈ N tal que para toda n ≥ N se tiene que
P(| ηn |> δ) < δ y por tanto

P(| ξn ηn |> ) ≤ P(| ξn ηn |> , | ηn |≤ δ) + δ (B.1)

1
= δ + P(| ξn |> ). (B.2)
δ
Sabemos que P(| ξn |> 1δ ) −−−→ P(| ξ |> 1δ ), entonces
n→∞

1
lı́m sup P(| ξn ηn |> ) ≤ δ + P(| ξ |> ).
n→∞ δ

Por lo tanto P(| ξn ηn |> ) −−−→ 0.

n→∞
92APÉNDICE B. CONVERGENCIA DE VARIABLES ALEATORIAS Y EL TEOREMA CENTRAL DEL LÍMITE.

Proposición B.8. Sean X̂, X̂1 , X̂2 , . . . vectores aleatorios, las siguientes afirmaciones son equivalentes

1. F (t̂) → F (t̂) para todo punto t̂ de continuidad.

X̂n X̂

2. Para toda función de valores reales f uniformemente continua y acotada E[f (X̂n )] → E[f (X̂)].

3. Para toda función de valores reales f continua y acotada E[f (X̂n )] → E[f (X̂)].

Si alguna de estas condiciones se cumple se dice que X̂n converge en distribución o converge débilmente
d
a X̂ que se escribe X̂n −
→ X̂

La demostración se encuentra en [? ] página 16.

Proposición B.9. Sean X̂, X̂1 , X̂2 , . . . y Ŷ, Ŷ1 , Ŷ2 , . . . vectores aleatorios sobre Rk y Rp respectivamente,
d p
tales que X̂n −
→ X̂ y Ŷn →
− Ŷ entonces
d
(X̂n , Ŷn ) − → (X̂, Ŷ)
Donde Ŷ es un vector constante

Demostración. Sea H : Rk+p → R una función uniformemente continua y acotada, entonces

| E[H(X̂n , Ŷn )] − E[H(X̂, Ŷ)] |≤ E[| H(X̂n , Ŷn ) − H(X̂n , Ŷ) |]+ | E[H(X̂n , Ŷ)] − E[H(X̂, Ŷ)] |

El primer sumando tiende a cero por que H es acotada y uniformemente continua, mientras que el segundo
tiende a cero por que Ŷ es constante y por la Proposición anterior

Los siguientes resultados pueden ser consultados en [? ]

Teorema B.1. (Ley fuerte de los grandes números). Sea ξ1 , ξ2 , . . . una sucesión de variables alea-
torias independientes e idénticamente distribuidas con media µ. Entonces
n
1 X c.s.
ξi −−→ µ.
n
i=1

Teorema B.2. (Teorema Central del Lı́mite). Sea ξ1 , ξ2 , . . . una sucesión de variables aleatorias
independientes e idénticamente distribuidas con media µ y varianza finita σ. Entonces
n
X ξi − µ d
√ −
→ N (0, 1).
nσ
i=1

Para la demostración del siguiente teorema ver [? ].

Teorema B.3. (Teorema de Curtiss). Sea ξ, ξ1 , ξ2 , . . . una sucesión de vectores aleatorios con funcio-
nes generadoras de momentos M, M1 , M2 , . . . finitas en algún abierto D que contenga al cero.
d
Si Mn (t) −−−→ M (t) para t ∈ D se tiene que ξn −−−→ ξ.
n→∞ n→∞
 Apéndice C

Códigos en R

En el código C.1 se realizan trayectorias que reciben como parámetros:

1. N ; tamaño de la población

2. t; tiempo de observación / tiempo de horizonte

3. beta; tasa de infección

4. gamma; tasa de recuperación

5. S0, I0, R0; condiciones iniciales.

Código C.1.
1 sir _ gillespie <- function ( S0 , I0 , R0 , beta , gamma , t ) {
2 S = S0
3 I = I0
4 R = R0
5 N=S+I+R
6 t _ ev =0
7 times <- c (0)
8 St <- c ( S )
9 It <- c ( I )
10 proba <- numeric (2)
11
12 while ( t _ ev < t & & I > 0) {
13
14 lambda = (( beta * S * I ) / N ) + ( gamma * I )
15 t _ ev <- t _ ev + rexp (1 , rate = lambda )
16 times <- append ( times , t _ ev )
17 pi = c (( beta * S * I / N ) / lambda , ( gamma * I ) / lambda )
18 eventos <- sample ( c (1 , -1) ,1 , prob = pi )
19
20 if ( eventos ==1) {
21 I <- I + eventos
22 S <- S - 1

93
 94 APÉNDICE C. CÓDIGOS EN R

23 R <- R
24 } else {
25 I <- I + eventos
26 S <- S
27 R <- R + 1
28 }
29 It <- append ( It , I )
30 St <- append ( St , S )
31 }
32 return ( data . frame ( It = It , St = St , times = times ) )
33 }
 Apéndice D

Códigos en Julia

El código D.1 crea un método llamado trayectoria que recibe como parámetros:

1. matrixtrans; es la matriz de transición de la cadena de Markov a simular, es de tamaño N × N

2. distinicial; es el vector de transición inicial de la cadena de Markov, es de tamaño N

3. tiempo; es el tiempo hasta el cual se quiere simular la cadena de Markov

4. E; es el espacio de estados de la cadena de Markov, es de tamaño N.

El método nos proporciona una trayectoria de la cadena de Markov, es decir un vector y nos muestra la
gráfica de la trayectoria simulada.
Código D.1.
1 Pkg . add ( " QuantEcon " )
2 Pkg . add ( " StatsBase " )
3 Pkg . add ( " Distributions " )
4 Pkg . add ( " PyPlot " )
5 using PyPlot
6 using StatsBase
7 using Distributions
8 using QuantEcon
9
10
11 function trayectoria ( matrixtrans , distinicial , tiempo , E )
12 dim = length ( E )
13 posicion = collect (1:1: dim )
14 estado =0
15 simulaciones = zeros ( tiempo +1)
16 estado = wsample ( posicion , distinicial ,1) [1]
17 simulaciones [1]= E [ estado ]
18

19 k =2
20 sim = zeros ( dim )
21 while k <= tiempo +1
22 h =1

95
 96 APÉNDICE D. CÓDIGOS EN JULIA

23 while h <= dim

24 sim [ h ]= matrixtrans [ estado , h ]
25 h = h +1
26 end
27 estado = wsample ( posicion , sim ,1) [1]
28 simulaciones [ k ]= E [ estado ]
29 k = k +1
30 end
31 x = collect (0:1: tiempo )
32

33
34 return hcat (x , simulaciones )
35 end
36
37 matrixtrans =[.2 .1 .5 .2;.5 .2 .2 .1;.2 .5 .1 .2;.3 .5 .1 .1]
38 distinicial =[.2 ,.4 ,.2 ,.2]
39 E =[1 ,2 ,3 ,4]
40 tiempo =30
41 trayecto = trayectoria ( matrixtrans , distinicial , tiempo , E )
42 plot ( trayecto [1: end ,1] , trayecto [1: end ,2] , " . " )
43 axis ([0 , tiempo , E [1] -1 , E [ end ]+1])
44 title ( " Trayectoria de una Cadena de Markov " )

El código D.2 crea un método llamado inf erencia que calcula los estimadores máximo verosı́miles de
las probabilidades de transición teniendo como información una trayectoria observada de una cadena de
Markov con espacio de estados E, recibe como parámetros:

1. trayecto; es una trayectoria de la cadena de Markov, es un vector

2. E; es el espacio de estados en el que vive la cadena de Markov, es de tamaño N.

El método nos proporciona como salida una matriz de tamaño N × N.

Código D.2.
1 Pkg . add (" QuantEcon ")
2 Pkg . add (" StatsBase ")
3 Pkg . add (" Distributions ")
4 Pkg . add (" PyPlot ")
5 using PyPlot
6 using StatsBase
7 using Distributions
8 using QuantEcon
9

10 function inferencia ( trayecto , E )

11 N = length ( E )
12 P = zeros (N , N )
13 n = zeros ( N )
14 t = length ( trayecto )
15 k =1
16 l =2
17 u =1
18 while ( trayecto [1]!= E [ k ])
19 k = k +1
 97

20 end
21 while (l <= t )
22 u =1
23 while ( trayecto [ l ]!= E [ u ])
24 u = u +1
25 end
26 P [k , u ]=1+ P [k , u ]
27 n [ k ]=1+ n [ k ]
28 k=u
29 l = l +1
30 end
31 i =1
32 j =1
33 while (i <= N )
34 j =1
35 while (j <= N )
36 P [i , j ]= P [i , j ]/ n [ i ]
37 j = j +1
38 end
39 i = i +1
40 end
41 return P
42
43 end
44 estimador = inferencia ( trayecto [1: end ,2] , E )

El código D.3 crea un método para analizar la convergencia de la distancia entre la matriz de transición
real y la estimada conforme el horizonte de tiempo converge a infinito, recibe los siguientes parámetros:

1. matrixtrans; es la matriz de transición de la cadena de Markov a simular, es de tamaño N × N

2. distinicial; es el vector de distribución inicial de la cadena de Markov, es de tamaño N

3. tiempo; es el tiempo hasta el cual se quiere analizar la distancia

4. E; es el espacio de estados de la cadena de Markov, es de tamaño N.

El método nos proporciona la gráfica de la convergencia de la distancia.

Código D.3.
1 function convergencia ( matrixtrans , distinicial , tiempo , E )
2 distancia = zeros (1)
3 vectortiempo = zeros (1)
4 t =100
5 simulaciones = trayectoria ( matrixtrans , distinicial , max ( tiempo ,100) ,E )
6 k =1
7 i =1
8 while (t <= max ( tiempo ,100) )
9 i =1
 98 APÉNDICE D. CÓDIGOS EN JULIA

10 while (t > simulaciones [i ,1])

11 i = i +1
12 end
13 estimador = inferencia ( simulaciones [1: i ,2] , E )
14 distancia [ k ]= norm ( estimador - matrixtrans ,2)
15 vectortiempo [ k ]= t
16 distancia = vcat ( distancia ,0)
17 vectortiempo = vcat ( vectortiempo ,0)
18 k = k +1
19 t = t +100
20 end
21
22 plot ( vectortiempo [1: end -1] , distancia [1: end -1])
23 hlines ( xmin =0 , y =0 , xmax = vectortiempo [ end -1])
24 axis ([0 , vectortiempo [ end -1] , -.01 , maximum ( distancia [1: end -1]) ])
25 xlabel (" tiempo ")
26 ylabel (" norma ")
27 title (" Convergencia del estimador ")
28 return hcat ( vectortiempo [1: end -1] , distancia [1: end -1])
29 end

El código D.4 crea un método llamado P SM que simula una cadena de Markov en tiempo continuo
(PSM), recibe como parámetros :

1. generador; es la matriz de intensidades (generador infinitesimal), es de tamaño N × N

2. distinicial; es el vector de distribución inicial del proceso, es de tamaño N

3. tiempo; es el tiempo hasta el cual se quiere simular el proceso

4. E; es el espacio de estados del proceso, es de tamaño N.

El método nos proporciona la trayectoria simulada.

Código D.4.
1
2
3 function PSM ( generador , distinicial , tiempo , E )
4 dim = length ( E )
5 posicion = collect (1:1: dim )
6 estado =0
7 saltos = zeros (1)
8 tiempos = zeros (1)
9 estado = wsample ( posicion , distinicial ,1) [1]
10 saltos [1]= E [ estado ]
11 tiempos [1]= log ( rand () ) / generador [ estado , estado ]
12 k = tiempos [1]
13 sim = zeros ( dim )
14 cont =2
15 while k <= tiempo
16 h =1
 99

17 while h <= dim

18 sim [ h ]= - generador [ estado , h ]/ generador [ estado , estado ]
19 h = h +1
20 sim [ estado ]=0
21 end
22 estado = wsample ( posicion , sim ,1) [1]
23 saltos = vcat ( saltos ,0)
24 saltos [ cont ]= E [ estado ]
25 tiempos = vcat ( tiempos ,0)
26 tiempos [ cont ]=( log ( rand () ) / generador [ estado , estado ]) + k
27 k = tiempos [ cont ]
28 cont = cont +1
29 end
30 ## A q u A~ definimos 2 nuevos arreglos llamados tiempoinicial y tiempofinal ,
31 ## estos representan en su i - esima entrada el tiempo inicial y tiempo final
32 ## en el que el proceso se mantuvo constante en el (i -1) - esimo salto .
33 tiempoinicial = zeros ( length ( saltos ) )
34 tiempofinal = zeros ( length ( saltos ) )
35 tiempofinal = tiempos [1:( length ( saltos ) ) ]
36 tiempoinicial [2:( length ( saltos ) ) ]= tiempos [1:( length ( saltos ) -1) ]
37 tiempofinal [ length ( saltos ) ]= tiempo
38 tiempoinicial [1]=0
39 return hcat ( saltos , tiempoinicial , tiempofinal )
40 end
41 generador =[ -.13 .04 .01 .08;.02 -.10 .04 .04;.2 .1 -.6 .3;.05 .15 .1 -.3]
42 distinicial =[.2 ,.2 ,.1 ,.5]
43 E =[1 ,2 ,3 ,4]
44 tiempo =50
45 simulaciones = PSM ( generador , distinicial , tiempo , E )
46 n = div ( length ( simulaciones ) ,3)
47
48 hlines ( y = simulaciones [1: n ,1] , xmin = simulaciones [1: n ,2] , xmax = simulaciones [1: n ,3] ,
color =" blue " , linewidth =3)
49 plot ( simulaciones [1: n ,2] , simulaciones [1: n ,1] ," bo " , color =" blue ")
50 plot ( simulaciones [1:( n -1) ,3] , simulaciones [1:( n -1) ,1] ," o " , color =" white ")
51 axis ([0 , tiempo , E [1] -1 , E [ end ]+1])
52 title (" Trayectoria del PSM ")

El código D.5 crea un método llamado inf erenciaP SM , calcula el estimador máximo verosı́mil del ge-
nerador infinitesimal de una trayectoria dada de un PSM, los parámetros del método son:

1. tiempos; es un vector que almacena los tiempos entre saltos, es decir los tiempos de inter arribo del
proceso.

2. estados; es un vector que almacena los estados visitados del proceso en cada transición

3. E; es el espacio de estados del proceso, es de tamaño N.

El método nos proporciona una matriz de N × N, la cual representa el estimador máximo verosı́mil del
generador infinitesimal de la trayectoria dada.
Código D.5.
1
2
 100 APÉNDICE D. CÓDIGOS EN JULIA

3 function inferenciaPSM ( tiempos , estados , E )

4 N = length ( E )
5 Q = zeros (N , N )
6 R = zeros ( N )
7 t = length ( estados )
8 k =1
9 l =2
10 u =1
11 while ( estados [1]!= E [ k ])
12 k = k +1
13 end
14
15 while (l <= t )
16 u =1
17 while ( estados [ l ]!= E [ u ])
18 u = u +1
19 end
20 Q [k , u ]=1+ Q [k , u ]
21 Q [k , k ]=1+ Q [k , k ]
22 R [ k ]= R [ k ]+ tiempos [l -1]
23 k=u
24 l = l +1
25 end
26 R [ u ]= R [ u ]+ tiempos [ t ]
27 i =1
28 while (i <= N )
29 j =1
30 while (j <= N )
31 Q [i , j ]= Q [i , j ]/ R [ i ]
32
33 j = j +1
34 end
35 Q [i , i ]= - Q [i , i ]
36 i = i +1
37 end
38 return Q
39
40 end
41
42 gene rador estima do = inferenciaPSM ( simulaciones [1: end ,3] - simulaciones [1: end ,2] ,
simulaciones [1: end ,1] , E )

El código D.6 crea un método llamado convergencia que sirve para analizar la convergencia de la distancia
entre el generador infinitesimal real y el estimado conforme el horizonte de tiempo tiende a infinito, recibe
los siguientes parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso, es de tamaño N × N

2. distinicial; es el vector de distribución inicial del proceso, es de tamaño N

3. tiempo; es el tiempo hasta el cual se quiere analizar la convergencia, al menos debe ser 100
 101

4. E; es el espacio de estados del proceso, es de tamaño N.

El método nos proporciona la gráfica de la convergencia de la distancia.

Código D.6.
1
2 function convergencia ( generador , distinicial , tiempo , E )
3 distancia = zeros (1)
4 vectortiempo = zeros (1)
5 t =100
6 d =0
7 simulaciones = PSM ( generador , distinicial , max ( tiempo ,100) ,E )
8 k =1
9 i =1
10 while (t <= max ( tiempo ,100) )
11 i =1
12 while (t > simulaciones [i ,3])
13 i = i +1
14 end
15 d = simulaciones [i ,3] - t
16 simulaciones [i ,3]= t
17 distancia [ k ]= norm ( inferenciaPSM ( simulaciones [1: i ,3] - simulaciones [1: i ,2] ,
simulaciones [1: i ,1] , E ) - generador ,2)
18 vectortiempo [ k ]= t
19 simulaciones [i ,3]= t + d
20 distancia = vcat ( distancia ,0)
21 vectortiempo = vcat ( vectortiempo ,0)
22 k = k +1
23 t = t +100
24 end
25

26 plot ( vectortiempo [1: end -1] , distancia [1: end -1])

27 hlines ( xmin =0 , y =0 , xmax = vectortiempo [ end -1])
28 axis ([0 , vectortiempo [ end -1] , -.01 , maximum ( distancia [1: end -1]) ])
29 xlabel (" tiempo ")
30 ylabel (" norma ")
31 title (" Convergencia del estimador ")
32 return hcat ( vectortiempo [1: end -1] , distancia [1: end -1])
33 end

El código D.7 crea un método llamado clasescomunicacion recibe como parámetros:

1. matrixtrans; es la matriz de transición de una cadena de Markov

2. E; es el espacio de estados de la cadena de Markov.

Como resultado nos proporciona todas las clases de comunicación de la cadena de Markov
Código D.7.
1 function c la se sco mu ni ca ci on ( matrixtrans , E )
2 mc = MarkovChain ( matrixtrans )
3 clases = c o m m u n i c a t i o n _ c l a s s e s ( mc )
 102 APÉNDICE D. CÓDIGOS EN JULIA

4 tam1 = length ( c o m m u n i c a t i o n _ c l a s s e s ( mc ) )
5 for i in 1: tam1
6 tam2 = length ( c o m m u n i c a t i o n _ c l a s s e s ( mc ) [ i ])
7 for j in 1: tam2
8 clases [ i ][ j ]= E [ clases [ i ][ j ]]
9 end
10 end
11
12 return clases
13 end
14
15 matrixtrans =[0 1 0 0;1/3 0 2/3 0;0 2/3 0 1/3;0 0 1 0]
16 E =[0 ,1 ,2 ,3]
17 cl as es co mu ni ca ci on ( matrixtrans , E )

El código D.8 crea un método llamado periodo recibe un sólo parámetro:

1. matrixtrans; es la matriz de transición de una cadena de Markov.

Como resultado nos muestra el periodo asociado a la cadena de Markov, suponiendo que dicha cadena es
irreducible.
Código D.8.
1 function periodo ( matrixtrans )
2 mc = MarkovChain ( matrixtrans )
3 return period ( mc )
4 end
5
6 matrixtrans =[0 1 0 0;1/3 0 2/3 0;0 2/3 0 1/3;0 0 1 0]
7 periodo ( matrixtrans )

El código D.9 crea un método llamado estadosrecurrentes recibe como parámetros:

1. matrixtrans; es la matriz de transición de una cadena de Markov

2. E; es el espacio de estados de la cadena de Markov.

Como resultado nos muestra el conjunto de estados de la cadena que son recurrentes.
Código D.9.
1 function e st ad osr ec ur re nt es ( matrixtrans , E )
2 mc = MarkovChain ( matrixtrans )
3 recurrentes = recu rrent_ classe s ( mc ) [1][1: end ]
4 for i in 1: length ( recurrentes )
5 recurrentes [ i ]= E [ recurrentes [ i ]]
6 end
7
8 return recurrentes
9 end
10
11 matrixtrans =[0 1 0 0;1/3 0 2/3 0;0 2/3 0 1/3;0 0 1 0]
12 E =[0 ,1 ,2 ,3]
13 es ta do sr ec ur re nt es ( matrixtrans , E )
 103

El código D.10 crea un método llamado estacionaria que sirve para dar una aproximación de la distribu-
ción estacionaria, recibe como parámetros:

1. matrixtrans; es la matriz de transición de una cadena de Markov

2. E; es el espacio de estados de la cadena de Markov

3. tiempo; es el tiempo hasta el cual se quiere analizar la convergencia a la distribución estacionaria.

Como resultado nos muestra la aproximación del vector estacionario.

Código D.10.
1 function estacionaria ( matrixtrans ,E , tiempo )
2 dim = length ( E )
3 nu = zeros ( dim )
4 distinicial = zeros ( dim )
5 distinicial [1]=1
6 trayecto = trayectoria ( matrixtrans , distinicial , tiempo , E ) [1: end ,2]
7 for i in 2: tiempo +1
8 for j in 1: dim
9 if trayecto [ i ]== E [ j ]
10 nu [ j ]= nu [ j ]+1
11 end
12

13 end
14
15 end
16 return nu / tiempo
17 end
18

19 tiempo =10000
20 matrixtrans =[0 1 0 0;1/3 0 2/3 0;0 2/3 0 1/3;0 0 1 0]
21 E =[0 ,1 ,2 ,3]
22 estacionaria ( matrixtrans ,E , tiempo )

El código D.11 crea un método llamado convergenciaestacionaria recibe como parámetros:

1. matrixtrans; es la matriz de transición de una cadena de Markov

2. E; es el espacio de estados de la cadena de Markov

3. tiempo; es el tiempo hasta el cual se quiere analizar la convergencia de la distribución estacionaria

4. vectorinvariante; es la distribución estacionaria de la cadena de Markov.

Como resultado nos muestra la gráfica de la convergencia a la distribución estacionaria.

 104 APÉNDICE D. CÓDIGOS EN JULIA

Código D.11.
1

2 function c o n v e r g e n c i a e s t a c i o n a r i a ( matrixtrans ,E , tiempo , vectorinvariante )

3 dim = length ( E )
4 nu = zeros ( dim )
5 d = div ( tiempo ,1000)
6 convergencia = zeros ( d )
7

8 distinicial = zeros ( dim )

9 distinicial [1]=1
10
11 trayecto = trayectoria ( matrixtrans , distinicial , tiempo , E ) [1: end ,2]
12 for k in 1: d
13 for i in 1000*( k -1) +2:1000* k +1
14 for j in 1: dim
15 if trayecto [ i ]== E [ j ]
16 nu [ j ]= nu [ j ]+1
17 end
18 end
19

20 end
21 convergencia [ k ]= norm (( nu /(1000* k ) ) - vectorinvariante ,2)
22 end
23
24 plot ( collect (1000:1000: tiempo ) , convergencia [1: end ])
25 hlines ( xmin =0 , y =0 , xmax = d *1000)
26 axis ([0 , d *1000 , -.01 , maximum ( convergencia [1: end ]) ])
27 xlabel (" tiempo ")
28 ylabel (" norma ")
29 title (" Convergencia al vector estacionario ")
30 return convergencia
31 end
32
33 tiempo =100000
34 matrixtrans =[0 1 0 0;1/3 0 2/3 0;0 2/3 0 1/3;0 0 1 0]
35 E =[0 ,1 ,2 ,3]
36 vectorinvariante =[1/8 ,3/8 ,3/8 ,1/8]
37 c o n v e r g e n c i a e s t a c i o n a r i a ( matrixtrans ,E , tiempo , vectorinvariante )

El código D.12 llamado probatransicion sirve para encontrar la matriz de probabilidades de transición al
tiempo t del PSM, recibe un sólo parámetro:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso, es de tamaño N × N.

La función muestra la diagonalización del generador infinitesimal, para después expresar a la matriz de
probabilidades de transición al tiempo t del PSM.
Código D.12.
1 function probatransicion ( generador )
2 dim = length ( generador [1: end ,1])
3 I = zeros ( dim , dim )
4 for i in 1: dim
5 I [i , i ]= eig ( generador ) [1][ i ]
6 end
 105

7 return Array [( eig ( generador ) [2]) ,I , inv ( eig ( generador ) [2]) ]

8 end
9 generador =[ -1 1 0 0;2 -4 2 0;0 1 -2 1;0 0 2 -2]
10 probatransicion ( generador )

El código D.13 crea un método llamado clasescomunicacionP SM recibe como parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso

2. E; es el espacio de estados del PSM.

Como resultado nos proporciona todas las clases de comunicación de la cadena de saltos asociada al PSM.
Código D.13.
1 function c l a s e s c o m u n i c a c i o n P S M ( generador , E )
2
3 d = length ( E )
4 matrixtrans = zeros (d , d )
5 k =0
6 for i in 1: d
7 k = generador [i , i ]
8 for j in 1: d
9 matrixtrans [i , j ]= - generador [i , j ]/ k
10 end
11 matrixtrans [i , i ]=0
12 end
13
14 mc = MarkovChain ( matrixtrans )
15 clases = c o m m u n i c a t i o n _ c l a s s e s ( mc )
16 tam1 = length ( c o m m u n i c a t i o n _ c l a s s e s ( mc ) )
17 for i in 1: tam1
18 tam2 = length ( c o m m u n i c a t i o n _ c l a s s e s ( mc ) [ i ])
19 for j in 1: tam2
20 clases [ i ][ j ]= E [ clases [ i ][ j ]]
21 end
22 end
23

24 return clases
25 end
26
27 generador =[ -1 1 0 0;2 -4 2 0;0 1 -2 1;0 0 2 -2]
28 E =[0 ,1 ,2 ,3]
29 cl a s e s c o m u n i c a c i o n P S M ( generador , E )

El código D.14 crea un método llamado estadosrecurrentesP SM recibe como parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del PSM

2. E; es el espacio de estados del PSM.

Como resultado nos muestra el conjunto de estados del PSM que son recurrentes.
 106 APÉNDICE D. CÓDIGOS EN JULIA

Código D.14.
1 function e s t a d o s r e c u r r e n t e s P S M ( generador , E )
2 d = length ( E )
3 matrixtrans = zeros (d , d )
4 k =0
5 for i in 1: d
6 k = generador [i , i ]
7 for j in 1: d
8 matrixtrans [i , j ]= - generador [i , j ]/ k
9 end
10 matrixtrans [i , i ]=0
11 end
12 mc = MarkovChain ( matrixtrans )
13 recurrentes = recu rrent_ classe s ( mc ) [1][1: end ]
14 for i in 1: length ( recurrentes )
15 recurrentes [ i ]= E [ recurrentes [ i ]]
16 end
17
18 return recurrentes
19 end
20
21 generador =[ -1 1 0 0;2 -4 2 0;0 1 -2 1;0 0 2 -2]
22 E =[0 ,1 ,2 ,3]
23 es t a d o s r e c u r r e n t e s P S M ( generador , E )

El código D.15 crea un método llamado estacionariaP SM recibe como parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso, es de tamaño N×N

2. tiempo; es el tiempo hasta el cual se desea calcular la aproximación de la distribución estacionaria

3. E; es el espacio de estados del PSM, es un vector de tamaño N.

Como resultado nos muestra la aproximación de la distribución estacionaria que es un vector de la misma
longitud que E.
Código D.15.
1

2 function estacionariaPSM ( generador , tiempo , E )

3 dim = length ( E )
4 distinicial = zeros ( dim )
5 distinicial [1]=1
6 posicion = collect (1:1: dim )
7 tiempoestancia = zeros ( dim )
8 estado =0
9 estado = wsample ( posicion , distinicial ,1) [1]
10 tiempoestancia [ estado ]= log ( rand () ) / generador [ estado , estado ]
11 k = tiempoestancia [ estado ]
12 s=k
13 sim = zeros ( dim )
14
15 while k <= tiempo
16 h =1
17 while h <= dim
 107

18 sim [ h ]= - generador [ estado , h ]/ generador [ estado , estado ]

19 h = h +1
20 sim [ estado ]=0
21 end
22 estado = wsample ( posicion , sim ,1) [1]
23 s =( log ( rand () ) / generador [ estado , estado ])
24 k=s+k
25 tiempoestancia [ estado ]= s + tiempoestancia [ estado ]
26 end
27 tiempoestancia [ estado ]= tiempoestancia [ estado ] - s +k - tiempo
28 return tiempoestancia / tiempo
29 end
30 generador =[ -1 1 0 0;2 -4 2 0;0 1 -2 1;0 0 2 -2]
31 E =[0 ,1 ,2 ,3]
32 tiempo =10000
33 estacionariaPSM ( generador , tiempo , E )

El código D.16 llamado convergenciaestacionariaP SM recibe como parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso, es de tamaño N×N

2. tiempo; es el tiempo hasta el cual se desea analizar la convergencia a la distribución estacionaria

3. E; es el espacio de estados del PSM, es un vector de tamaño N

4. vectorinvariante; es la distribución estacionaria del PSM, es un vector de tamaño N.

Como resultado nos muestra la gráfica de la convergencia hacia la distribución estacionaria.

Código D.16.
1 function c o n v e r g e n c i a e s t a c i o n a r i a P S M ( generador , tiempo ,E , vectorinvariante )
2 dim = length ( E )
3 d = div ( tiempo ,1000)
4 convergencia = zeros ( d )
5 distinicial = zeros ( dim )
6 distinicial [1]=1
7 posicion = collect (1:1: dim )
8 tiempoestancia = zeros ( dim )
9 estado =0
10 estado = wsample ( posicion , distinicial ,1) [1]
11 tiempoestancia [ estado ]= log ( rand () ) / generador [ estado , estado ]
12 k = tiempoestancia [ estado ]
13 s=k
14 sim = zeros ( dim )
15 for l in 1: d
16 while k <= l *1000
17 h =1
18 while h <= dim
19 sim [ h ]= - generador [ estado , h ]/ generador [ estado , estado ]
20 h = h +1
21 sim [ estado ]=0
 108 APÉNDICE D. CÓDIGOS EN JULIA

22 end
23 estado = wsample ( posicion , sim ,1) [1]
24 s =( log ( rand () ) / generador [ estado , estado ])
25 k=s+k
26 tiempoestancia [ estado ]= s + tiempoestancia [ estado ]
27 end
28 tiempoestancia [ estado ]= tiempoestancia [ estado ] - s +k - l *1000
29 convergencia [ l ]= norm ( vectorinvariante -( tiempoestancia /( l *1000) ) ,2)
30 tiempoestancia [ estado ]= tiempoestancia [ estado ]+ s - k + l *1000
31 end
32
33 plot ( collect (1000:1000: tiempo ) , convergencia [1: end ])
34 hlines ( xmin =0 , y =0 , xmax = d *1000)
35 axis ([0 , d *1000 , -.01 , maximum ( convergencia [1: end ]) ])
36 xlabel (" tiempo ")
37 ylabel (" norma ")
38 title (" Convergencia al vector estacionario ")
39 return convergencia
40
41 return convergencia
42 end
43

44
45 generador =[ -1 1 0 0;2 -4 2 0;0 1 -2 1;0 0 2 -2]
46 E =[0 ,1 ,2 ,3]
47 tiempo =100000
48 vectorinvariante =[1/3 ,1/6 ,1/3 ,1/6]
49 c o n v e r g e n c i a e s t a c i o n a r i a P S M ( generador , tiempo ,E , vectorinvariante )

El código D.17 crea un método llamado asintotica que recibe los siguientes parámetros:

1. matrixtrans; es la matriz de transición de la Cadena de Markov, es de tamaño N × N

2. distinicial; es el vector de la distribución inicial de la cadena de Markov, es de tamaño N

3. tiempo; es el tiempo hasta el cual se quiere aproximar a la matriz Σ̃

4. E; es el espacio de estados de la cadena de Markov, es de tamaño N.

El método nos proporciona la estimación de Σ̃ (ver Teorema 4.3).

Código D.17.
1 function asintotica ( matrixtrans , distinicial , tiempo , E )
2 trayecto = trayectoria ( matrixtrans , distinicial , tiempo , E ) [1: end ,2]
3 N = length ( E )
4 matrixFisher = zeros ( N ^2 -N , N ^2 - N )
5 P = zeros (N , N )
6 n = zeros ( N )
7 t = length ( trayecto )
8 k =1
9 l =2
 109

10 u =1
11 while ( trayecto [1]!= E [ k ])
12 k = k +1
13 end
14 while (l <= t )
15 u =1
16 while ( trayecto [ l ]!= E [ u ])
17 u = u +1
18 end
19 P [k , u ]=1+ P [k , u ]
20 n [ k ]=1+ n [ k ]
21 k=u
22 l = l +1
23 end
24
25 for i in 1: N
26 for j1 in 1: N -1
27 for j2 in 1: N -1
28 if j1 == j2
29 matrixFisher [(( i -1) *( N -1) ) + j1 ,(( i -1) *( N -1) ) + j2 ]=(( P [i , j1 ]/ n [ i ]) -( P [i , j1 ]/ n [ i ])
^2) /( n [ i ]/ tiempo )
30 else
31 matrixFisher [(( i -1) *( N -1) ) + j1 ,(( i -1) *( N -1) ) + j2 ]=( -( P [i , j1 ]/ n [ i ]) ^2) /( n [ i ]/
tiempo )
32 end
33 end
34 end
35 end
36 return matrixFisher
37
38 end
39 E =[1 ,2]
40 matrixtrans =[.2 .8 ; .8 .2]
41 distinicial =[.5 ,.5]
42 tiempo =100000
43 asintotica ( matrixtrans , distinicial , tiempo , E )

El código D.18 crea un método llamado asintoticaP SM recibe como parámetros:

1. generador; es la matriz de intensidades (generador infinitesimal) del proceso, es de tamaño N× N

2. distinicial; es la distribución inicial del PSM, es de tamaño N

3. E; es el espacio de estados del PSM, es un vector de tamaño N

4. tiempo; es el tiempo hasta el cual se quiere aproximar a Λ.

Como resultado nos muestra la estimación de Λ (ver Teorema 4.5).

Código D.18.
1
 110 APÉNDICE D. CÓDIGOS EN JULIA

2
3 function asintoticaPSM ( generador , distinicial , tiempo , E )
4 simulaciones = PSM ( generador , distinicial , tiempo , E )
5 tiempos = simulaciones [1: end ,3] - simulaciones [1: end ,2]
6 estados = simulaciones [1: end ,1]
7 N = length ( E )
8 matrixFisher = zeros ( N ^2 -N , N ^2 - N )
9 Q = zeros (N , N )
10 R = zeros ( N )
11 t = length ( estados )
12 k =1
13 l =2
14 u =1
15 while ( estados [1]!= E [ k ])
16 k = k +1
17 end
18
19 while (l <= t )
20 u =1
21 while ( estados [ l ]!= E [ u ])
22 u = u +1
23 end
24 Q [k , u ]=1+ Q [k , u ]
25 Q [k , k ]=1+ Q [k , k ]
26 R [ k ]= R [ k ]+ tiempos [l -1]
27 k=u
28 l = l +1
29 end
30 R [ u ]= R [ u ]+ tiempos [ t ]
31
32
33 for i in 1: N
34
35 j1 =1
36 l =0
37 while j1 <= N
38 if i == j1
39 l =1
40 j1 = j1 +1
41 else
42 matrixFisher [(( i -1) *( N -1) ) + j1 -l ,(( i -1) *( N -1) ) + j1 - l ]=( Q [i , j1 ]/ R [ i ]) ^2/( Q [i , j1 ]/
tiempo )
43 j1 = j1 +1
44 end
45
46 end
47 end
48 return matrixFisher
49 end
50
51
52 E =[1 ,2 ,3]
53 distinicial =[1/3 ,1/3 ,1/3]
54 generador =[ -2 1 1 ;1 -2 1;1 1 -2 ]
55 tiempo =10000
56 asintoticaPSM ( generador , distinicial , tiempo , E )
 111

Listing D.1: Puente Markoviano

1 puente _ rej <- function (t ,a ,b , Q ) {
2 n <- nrow ( Q )
3 mu <- rep (0 , n )
4 mu [ a ] <- 1
5 k <- 0
6 while ( T ) {
7 k <- k +1
8 trayectoria <- ctmc (Q , mu , t )
9 if ( tail ( trayectoria $ Xn ,1) == b ) {
10 break
11 }
12 }
13 print ( k )
14 return ( trayectoria )
15 }
112 APÉNDICE D. CÓDIGOS EN JULIA
Bibliografı́a

[1] S. Asmussen and P. W. Glynn. Stochastic Simulation: Algorithms and Analysis. Springer, Stochastic
Modelling and Applied Probability, 2007.

[2] N. Bailey. The Mathematical Theory of Epidemics. Grin, London, 1957.

[3] M. Bladt and M. Sõrensen. Simple simulation of diffusion bridges with application to likelihood
inference for diffusions. Bernoulli, 2(20):645–675, 2014.

[4] G. U. Bravo. Procesos Estocásticos.

[5] R. Brown. A brief account of microscopical observations made in the months of june, july, and
august, 1827, on the particles contained in the pollen of plants; and on the general existence of active
molecules in organic and inorganic bodies. Philosophical Magazine N. S., (4):161–173, 1828.

[6] I. J. R. S. Cox, J.C. A theory of the term structure of interest rates,. Econometrica, (53):385–408,
1985.

[7] D. Dacunha-Castelle and D. Florens-Zmirou. Estimation of the coefficients of a diffusion from discrete
observations. Stochastics,, 19:263–284, 1986.

[8] A. Einstein. Investigations on the theory of the Brownian movement. Dover, 1956.

[9] W. Feller. An Introduction to Probability Theory and Its Applications, volume 1. Wiley Series in
Probability and Mathematical Statistics, 3 edition, 196.

[10] S. Iacus. Simulation and Inference for Stochastic Differential Equations with R Examples. Springer
Series in Statistics, 2008.

[11] Y. Kutoyants. Statistical Inference for Ergodic Diffusion Processes. Springer-Verlag, London, 2004.

[12] P. É. Lapayre, B. and R. Sentis. Introduction to Monte-Carlo Methods for Transportation and
Diffusion Equations. Oxford University Press, New York, 2003.

[13] R. Merton. Theory of rational option pricing. Bell J. Econ. Manage. Sci., 4(1):141–183, 1973.

[14] N. Newton. Variance reduction for simulated diffusions. SIAMJ.Appl. Math, 54(6):1780–1805, 1994.

[15] B. Ø ksendan. Stochastic Differential Equations. 1998.

[16] N. Yoshida. Estimation for diffusion processes from discrete observation. J. Multivar. Anal.,
41(2):220–242, 1992.

113