Cadeias de Markov
Fernando Lucambio Pérez
2025-09-10
6. Inferência
Uma Cadeia de Markov é por vezes um modelo probabilístico adequado para determinada série temporal em que a observação em um determinado momento é uma categoria à qual um indivíduo corresponde. A mais simples Cadeia de Markov é aquela na qual existe um número finito de estados ou categorias, um número finito de pontos no tempo equidistantes em que são feitas as observações, a cadeia é de primeira ordem e as probabilidades de transição são as mesmas para cada intervalo de tempo. Estas foram as cadeias que estudamos até agora.
Vamos considerar aqui vários métodos para obter estimadores da matriz de probabilidades transição em três situações:
quando o comprimento do intervalo entre os pontos de tempo do modelo coincide com o intervalo de observação,
quando a duração deste comprimento do intervalo entre os pontos de tempo não coincide com o intervalo de observação e
quando os intervalos de observação são desiguais em comprimento.
Além disso, discutimos o uso da técnica de bootstraps como um método para avaliar a incerteza nas estimaçães e para a construção de intervalos de confiança da matriz de transição. Também estudaremos como verificar a ordem de uma Cadeia de Markov. Depois apresentamos testes para verificar a ordem da cadeia. Diferentes livros e artigos foram consultados para escrever este texto, alguns importantes são Jacobsen (1982), Billingsley (1961a), Billingsley (1961b) e Basawa and Prakasa Rao (1980). Outras referências serão também mencionadas no texto.
6.1. Matriz de transição
Seja \(\{C_1,C_2,\cdots\}\) um processo estocástico ou uma sequência de variáveis aleatórias assumindo valores em algum conjunto finito, chamado aqui de espaço de estados \(\mathcal{S}\). A variável \(C_n\) deve ser considerada como estado no tempo \(n\) de um sistema cuja evolução é regida por um conjunto de leis de probabilidade.
A típica Cadeia de Markov a tempo discreto limita a descrição do histórico de cada sujeito a pontos de tempo com igualmente espaçados. Em outras palavras, em vez de modelar a possibilidade de progressão a cada instante no tempo, ou seja, dia, mês ou ano. O intervalo entre esses pontos de tempo é conhecido como o comprimento do ciclo.
Nesta seção, presume-se que a matriz de probabilidades de transição vai ser estimada a partir de dados de coorte longitudinais, com intervalos de observação comuns a todos os sujeitos. A atenção é restrita a obtenção da estimativa de máxima verossimilhança da matriz de transição em três situações específicas crescentes de complexidade:
O primeiro caso é quando os intervalos de observação são constantes e coincidem com a duração do ciclo.
O segundo caso acontece quando os intervalos de observação são constantes, mas não coincidem com a duração do ciclo. Os métodos discutidos na presente seção somente podem ser utilizados em certas situações. Quando não puder, o método discutido para o terceiro caso é possível.
O terceiro caso, representa a situação mais comum, quando os intervalos de observação não são iguais em comprimento. A duração do ciclo pode ou não coincidir com um destes intervalos.
Vamos considerar \(\{c_0,c_1,\cdots,c_{n} \}\) uma amostra de uma Cadeia de Markov com probabilidades de transição \(p_{x,y}\) e distribuição inicial \(\pi_0\). Observe que \(\{c_0,c_1,\cdots,c_n \}\) deve ser uma sequência de \(n+1\) estados.
Então, a probabilidade de que \(c_0,c_1,\cdots,c_n\) seja essa sequência é justamente \[ \pi(x_0)p_{c_0,c_1}\cdots p_{c_{n-1},c_n}\cdot \]
Para \(x,y=1,2,\cdots,d\), seja \(n_{x,y}\) o número de transições, assim a matriz \((n_{x,y})\) vai ser chamada de matriz de contagens de transições da sequência. Dado que \[ \pi(c_0)p_{c_0,c_1}\cdots p_{c_{n-1},c_n} \, = \, \pi(c_0) \prod_{x,y} p_{x,y}^{n_{x,y}}, \] a contagem das transições junto com o estado inicial formam uma estatística suficiente. A distribuição dessa estatística vai ser nosso objetivo.
Sabemos que a probabilidade de obtiver uma sequência em particular, que começe com com \(x_0\) e tenha matriz de contagens de transição \((n_{x,y})\) foi dada acima e, com o objetivo de encontrarmos a distribuição da estatística suficiente é necessário somente contar o número de tais sequências.
Se \(n_{x,\cdot}=\sum_y n_{x,y}\) e \(n_{\cdot, y}=\sum_x n_{x,y}\) então \(\{n_{x,\cdot}\}\) e \(\{n_{\cdot, y}\}\) são as frequências das contagens de \(\{c_0,c_1,\cdots,c_{n-1}\}\) e \(\{c_1,c_2,\cdots,c_n\}\), respectivamente. Disso seque que \[ \begin{array}{rcl} n_{x,\cdot} \, - \, n_{\cdot, y} & = & \pmb{1}_{x}(x_0)-\pmb{1}_{x}(x_n) \\[0.6em] \displaystyle \sum_{x,y} n_{x,y} & = & \displaystyle \sum_x n_{x,\cdot} \, = \, \sum_y n_{\cdot, y} \, = \, n \cdot \end{array} \]
É claro, a partir da primeira dessas relações, que \((n_{x,y})\) e o estado inicial determinam completamente o estado final. Da mesma forma, \((n_{x,y})\) e o estado do final determinam o estado inicial. No entanto, \((n_{x,y})\) sozinho não determina os estados inicial e final: por exemplo, as sequências \(\{1,2,1 \}\) e \(\{2,1,2 \}\) têm contagens de transição idênticas. A resposta a este problema combinatório foi fornecida por Whittle (1961).
Teorema 6.1. Fórmula de Whittle
Seja \((n_{x,y})\) uma matriz \(d\times d\) de inteiros não negativos satisfazendo que \(\sum_{xy} n_{xy}=n\) e tais que \[ n_{x,\cdot} - n_{\cdot, y} \, = \, \pmb{1}_{x}(u) \, - \, \pmb{1}_{x}(v), \] \(x,y=1,\cdots,d\) para algum par \(u,v\). Se \(N_{u,v}^{(n)}(n_{x,y})\) representa o número de sequências \(\{x_0,x_1,\cdots,x_n\}\) tendo contagens de transição \((n_{x,y})\) e satisfazendo \(x_0=u\) e \(x_n=v\), então \[ N_{u,v}^{(n)}(n_{x,y}) \, = \, \dfrac{\displaystyle \prod_x n_{x\cdot}!}{\displaystyle \prod_{x,y} n_{x,y}!}C_{v,u}, \] onde \(C_{v,u}\) é o cofator \((v,u)\) da matriz \((n_{x,y})^*\) de componentes \[ n_{x,y}^* \, = \, \left\{ \begin{array}{lcl} \pmb{1}_{x}(y) - \dfrac{n_{x,y}}{n_{x,\cdot}} & \mbox{se} & n_{x,\cdot}>0 \\[0.8em] \pmb{1}_{x}(y) & \mbox{se} & n_{x,\cdot}=0 \end{array}\right.\cdot \]
Demonstração.
Esta forma de demonstrar a Fórmula de Whittle deve-se à Billingsley (1961b). A demonstração é por indução. O resultado é fácil de estabelecer se \(n=1\), caso em que ambos os lados de são 1. Se \((n_{u,v})\) é \((n_{x,y})\) com a \((u,v)\) entrada diminuída em 1, temos que \[ N_{u,v}^{(n)}(n_{x,y}) \, = \, \sum_w N_{w,v}^{(n-1)}(n_{u,w}), \] onde a soma se estende sobre aqueles \(w\) para so quais \(n_{u,w}>0\). Por isso, basta mostrar que o lado direito de \(N_{u,v}^{(n)}(n_{x,y})\) satisfaz esta mesma relação ou que \[ (n_{x,y})^* \, = \, \sum_w n_{u,w}n_{u,\cdot}^{-1}(n_{v,w})^*(u,w)\cdot \] Desde que \((n_{v,w})^*\) e \((n_{x,y})^*\) concordem fora da \(w\)-ésima coluna, \((n_{v,w})^*(u,w)=(n_{v,w})^*\). Com este fato, juntamente com a definição, segue que \((n_{x,y})^*\) é equivalente a \(\sum_w n_{u,w}^*(n_{v,w})^*=0\) onde a soma se estende sobre todos os \(w\). Dado que \(\sum_w n_{u,w}^*(n_{v,w})^*=\pmb{1}_{u}(v)|(n_{x,y})^*|\), a expressão recorrente acima vale para o caso no qual \(u\neq v\) e é necessário somente mostrar que \(|(n_{x,y})^*|=0\) caso \(u=v\). Suponhamos convenientemente que \(n_{x,\cdot}=n_{\cdot,x}\) seja positivo para \(x\leq r\) e zero para \(x>r\). Então \((n_{x,y})\) tem a forma \[ (n_{x,y})=\begin{pmatrix} A & 0 \\ 0 & 0 \end{pmatrix}, \] onde \(A\) é uma matriz \(r\times r\). Pela definição, \[ (n_{x,y})^*=\begin{pmatrix} A^* & 0 \\ 0 & \mbox{I} \end{pmatrix}, \] onde as linhas de \(A^*\) somam zero. Por isso, \(|(n_{x,y})^*|=|A^*|=0\).
Exemplo 6.1: Cadeia com dois estados.
Seja, por exemplo, a sequência de 12 valores observados \(\{0,1,1,0,1,0,1,1,1,0,0,1\}\). Esta sequência tem \(u=0\), \(v=1\) e matriz de contagens de transição \[ (n_{x,y}) \, = \, \begin{pmatrix} 1 & 4 \\ 3 & 3 \end{pmatrix}\cdot \]
Podemos utilizar a seguinte função R para encontrarmos a matriz de contagens de transição:
x = c(0,1,1,0,1,0,1,1,1,0,0,1)
library(markovchain); library(matlab); library(matlib)
Matriz = createSequenceMatrix(x, sanitize=FALSE)
Matriz## 0 1
## 0 1 4
## 1 3 3Vemos, do Teorema 5.1, que \[ n_{x,y}^* \, = \, \begin{pmatrix} \dfrac{4}{5} & -\dfrac{4}{5} \\[0.8em] -\dfrac{1}{2} & \dfrac{1}{2} \end{pmatrix} \] e \(C_{0,1}=4/5\). Temos que \[ N_{0,1}^{(12)}(n_{x,y}) \, = \, \dfrac{5!\cdot 6!}{1!\cdot 3!\cdot 3!\cdot 4!}\cdot \dfrac{4}{5} \, = \, 80\cdot \] Logo, 80 é o número de sequências \(\{0,x_1,\cdots,x_{10},1\}\) tendo contagens de transição \((n_{x,y})\).
Desenvolvemos uma função R para encontrarmos o número de sequências \(\{c_0,c_1,\cdots,c_n\}\) tendo contagens de transição \((n_{x,y})\) e satisfazendo \(c_0=u\) e \(c_n=v\):
Whittle = function(M, u, v){
n = length(rowSums(M))
Prod1 = 1;
for(i in 1:n) Prod1 = Prod1*gamma(rowSums(M)[[i]]+1)
Prod2 = 1;
for(i in 1:dim(M)[1]){
for(k in 1:dim(M)[2]) Prod2 = Prod2*gamma(M[i,k]+1)
}
u = which(row.names(M) == u)
v = which(row.names(M) == v)
C = cofactor(eye(n)-M/rowSums(M), v, u)
return((Prod1/Prod2)*C)
}
Whittle(Matriz, 0, 1)## [1] 806.1.1. Intervalos coincidentes
Suponhamos que nos seja dada a realização de uma Cadeia de Markov e que se deseja estimar a matriz de probabilidades de transição. Uma abordagem é encontrar as contagens de transição e estimar as probabilidades de transição de uma forma óbvia.
Exemplo 6.2: Cadeia com três estados
Esta seria uma situação hipotética. Consideremos uma Cadeia de Markov com três estados da qual é observada a sequência:
seq = c(2,3,3,2,1,1,1,1,1,2,2,1,3,1,3,2,3,3,2,1,2,2,2,2,3,2,3,2,3,
3,2,2,2,2,2,1,3,1,3,2,3,3,2,2,1,2,2,1,3,2,3,2,1,3,2,2,3,2,
3,1,3,2,3,3,2,2,2,3,2,1,3,2,3,2,3,3,1,2,3,2,2,2,3,2,3,2,3,
3,1,2,2,2,1,2,3,2,3,2,1,3,2,1,2,3,2,3,3,1,3,2,3,3,2,1,2,1)Por simples contagem, segue-se que, a matriz do número de transições entre os estados é: \[ (n_{x,y})=\begin{pmatrix} 4 & 8 & 10 \\ 13 & 17 & 22 \\ 6 & 26 & 9 \end{pmatrix}, \] onde \(n_{x,y}\) denota o número de transições observadas desde o estado \(x\) ao estado \(y\).
Uma vez que o número de transições do estado 2 para o estado 3 é 22 e o número total de transições do estado 2 é \(13+17+22\), uma estimativa empírica de \(\widehat{p}_{2,3}\) é 22/52. Uma estimativa empírica para a matriz de transições seria então \[ \widehat{\mathcal{P}}=\begin{pmatrix} \dfrac{4}{22} & \dfrac{8}{22} & \dfrac{10}{22} \\[0.8em] \dfrac{13}{52} & \dfrac{17}{52} & \dfrac{22}{52} \\[0.8em] \dfrac{6}{41} & \dfrac{26}{41} & \dfrac{9}{41} \end{pmatrix} \cdot \]
Vamos agora mostrar que este é, de fato, a estimativa de máxima verossimilhança condicional de \(\mathcal{P}\), condicionada à primeira observação. Suponhamos, então, que nós queremos estimar os \(d^2-d\) parámetros de uma Cadeia de Markov \(\{C_n\}\) com \(d\) estados a partir uma realização \(c_0,c_1,\cdots,c_T\).
A função de verossimilhança condicional à primeira observação é \[ L=\prod_{x=1}^d \prod_{y=1}^d p_{x,y}^{n_{x,y}} \cdot \]
Desta expressão obtemos que o logaritmo da função de verossimilhança é \[ \ell=\sum_{x=1}^d \left( \sum_{y=1}^d n_{x,y}\log\big( p_{x,y}\big) \right)= \sum_{x=1}^d \ell_x, \] a qual podemos maximizar maximizando cada somando separadamente.
Substituindo \(\displaystyle 1-\sum_{z\neq x} p_{x,z}\) por \(p_{x,x}\), diferenciando cada \(\ell_x\) com relação à todas as probabilidades de transição não diagonais \(p_{x,y}\) e igualando as derivadas a zero obtemos \[ 0=\dfrac{-n_{x,x}}{\displaystyle 1-\sum_{z\neq x} p_{x,z}}+\dfrac{n_{x,y}}{ p_{x,y}}=-\dfrac{n_{x,x}}{ p_{x,x}}+\dfrac{n_{x,y}}{ p_{x,y}} \cdot \] Assim, a menos que um denominador seja zero na equação acima \[ n_{x,y} p_{x,x}=n_{x,x} p_{x,y}, \] e por isso \(\displaystyle p_{x,x}\sum_{y=1}^d n_{x,y}=n_{xx}\).
Isto implica que o ponto de máximo local da função de verossimilhança das probabilidades de transição é \[ \widehat{p}_{x,x}=\dfrac{n_{x,x}}{\displaystyle \sum_{y=1}^d n_{x,y}} \quad \mbox{e} \quad \widehat{p}_{x,y}=\dfrac{n_{x,y}}{\displaystyle \sum_{y=1}^d n_{x,y}} \cdot \]
Também poderíamos utilizar multiplicadores de Lagrange para expressar as restrições \(\sum_{y=1}^d p_{x,y}=1\), sob as quais buscamos maximizar os termos de \(\ell_x\) e, portanto, a função de verossimilhança. Mas isso, não é necessário em geral.
Exemplo 6.3.
Fazendo uso do pacote de funções markovchain, construímos a matriz de contagens de transições utilizando o comando createSequenceMatrix, como mostrado a seguir:
Matriz = createSequenceMatrix(seq, sanitize=FALSE)
Matriz## 1 2 3
## 1 4 8 10
## 2 13 17 22
## 3 6 26 9Às vezes o número de Whittle é impraticável, como nesta situação. Utilizando a função anterior obtemos que:
Whittle(Matriz, u = 2, v = 1)## [1] 8.462769e+44A questão agora é transformar as frequências observadas em probabilidades, para isso utilizamos o comando markovchainFit do qual temos por resposta uma lista com diversas informações. A primeira resposta é a matriz de probabilidades de transição estimada, a qual pode ser obtida também digitando mcFitMLE[[1]].
mcFitMLE = markovchainFit(data=seq)
mcFitMLE$estimate## MLE Fit
## A 3 - dimensional discrete Markov Chain defined by the following states:
## 1, 2, 3
## The transition matrix (by rows) is defined as follows:
## 1 2 3
## 1 0.1818182 0.3636364 0.4545455
## 2 0.2500000 0.3269231 0.4230769
## 3 0.1463415 0.6341463 0.2195122mcFitMLE[[1]]## MLE Fit
## A 3 - dimensional discrete Markov Chain defined by the following states:
## 1, 2, 3
## The transition matrix (by rows) is defined as follows:
## 1 2 3
## 1 0.1818182 0.3636364 0.4545455
## 2 0.2500000 0.3269231 0.4230769
## 3 0.1463415 0.6341463 0.2195122O valor da fórmula de Whittle as vezes é tão grande que não pode ser calculado com precisão. Portanto, em vez disso, podemos aproximar o número de sequências \(\{c_0,c_1,\cdots,c_n\}\); tendo contagens de transição \((n_{x,y})\), satisfazendo \(c_0=u\) e \(c_n=\nu\), aplicando uma série de Stirling para os termos fatoriais \[ \ln(x!) \approx x\ln(x)-x+\frac{1}{2}\ln(2\pi x)+\frac{1}{12x}-\frac{1}{360x^3}+\frac{1}{1260x^5}-\frac{1}{1680x^7}, \] válido sempre que \(x>16\).
Exemplo 6.4.
Desenvolvemos as seguintes funções R para aproximar o número de Whittle, aplicando para isso a série de Stirling acima aos termos fatoriais:
# Aproximação de Stirling para o logaritmo do fatorial
fact = function(x) {
fact = x*log(x)-x+0.5*log(2*pi*x)+1/(12*x)-1/(360*x^3)+1/(1260*x^5)-1/(1680*x^7)
return(exp(fact))
}
Whittle1 = function(M, u, v){
n = length(rowSums(M))
Prod1 = 1;
for(i in 1:n) Prod1 = Prod1*fact(rowSums(M)[[i]])
Prod2 = 1;
for(i in 1:dim(M)[1]){
for(k in 1:dim(M)[2]) Prod2 = Prod2*fact(M[i,k])
}
u = which(row.names(M) == u)
v = which(row.names(M) == v)
C = cofactor(eye(n)-M/rowSums(M), v, u)
return((Prod1/Prod2)*C)
} Fazendo uso da função desenvolvida, obtemos que:
library(markovchain); library(matlab); library(matlib)
x = c(0,1,1,0,1,0,1,1,1,0,0,1)
Matriz1 = createSequenceMatrix(x, sanitize=FALSE)
Matriz1## 0 1
## 0 1 4
## 1 3 3Whittle1(Matriz1, 0, 1)## [1] 80.02461mas, caso a sequência seja muito longa, o valor do número de possíveis sequências continua explossivo.
Whittle1(Matriz, 2, 1)## [1] 8.462769e+44Teorema 6.2.
Seja \(\{C_n\}\) uma Cadeia de Markov ergódica. Então, independentemente da distribuição inicial \[ \sqrt{n_{x,y}}(\widehat{p}_{x,y}-p_{x,y}) \overset{D}{\longrightarrow} Z, \] onde \(Z\sim N(0,\Sigma)\), \(\Sigma=(\sigma_{x,y})\), \(x,y\in\mathcal{S}\) e \[ \sigma_{x,y}=\left\{\begin{array}{lcl} p_{x,y}(1-p_{x,y}) & \quad \mbox{caso} \quad \{x,y\}\ne \{z,w\} \\ -p_{x,y}p_{x,z} & \quad \mbox{caso} \quad x=z, \; y\ne w \\ 0 & \quad \mbox{caso contrário} \end{array} \right. \]
Demonstração.
Consequência do Teorema Central do Limite (ver Anderson and Goodman 1957).O resultado deste teorema implica que a covariância assintótica tem uma estrutura multinomial dentro das linhas e independência entre as linhas. Como resposta temos também o desvio padrão da estimativa, calculado segundo a expressão acima, como \(\widehat{p}_{x,y}/\sqrt{n_{x,y}}\), assim como um intervalo de confiança de 95%, os limites inferior e superior deste intervalo e o valor da função de log-verossimilhança.
Exemplo 6.5.
Mostramos agora a forma de obtermos os resultados apresentados no Teorema 6.2. Todos os resultados estão guardados no objeto mcFitMLE e para sabermos o valor do desvio padrão, por exemplo, digitamos:
mcFitMLE$standardError## 1 2 3
## 1 0.09090909 0.12856487 0.14373989
## 2 0.06933752 0.07929049 0.09020030
## 3 0.05974365 0.12436633 0.07317073Os intervalos confidenciais são de 95%, o qual verificamos digitando:
mcFitMLE$confidenceLevel## [1] 0.95Então, por final, obtemos os limites inferior e superior do intervalo de confiança como mostrado abaixo. Ainda mostramos o valor da função \(\ell\), o logaritmo da função de verosimilhança.
mcFitMLE$lowerEndpointMatrix## 1 2 3
## 1 0.003639601 0.1116538 0.17282038
## 2 0.114100921 0.1715165 0.24628754
## 3 0.029246032 0.3903928 0.07610017mcFitMLE$upperEndpointMatrix## 1 2 3
## 1 0.3599968 0.6156189 0.7362705
## 2 0.3858991 0.4823296 0.5998663
## 3 0.2634369 0.8778999 0.3629242mcFitMLE$logLikelihood## [1] -115.7695O nosso próximo exemplo é uma aplicação de Cadeias de Markov em engenharia; tem a ver com pontes e baseia-se no trabalho de Skuriat-Olechnowska (2005).
Exemplo 6.6: Inspeção de pontes
Exemplo desenvolvido em Skuriat-Olechnowska (2005). A maioria das pontes na Holanda é construída em concreto e mais de metade delas tem mais de 30 anos. À medida que as pontes se deterioram a uma velocidade acelerada devido à corrosão, à degradação do concreto e ao dano do veículo a Divisão de Engenharia Civil de Rijkswaterstaat, deve repará-las e, sempre que possível, assim impedem uma maior deterioração.
A Divisão de Engenharia Civil de Rijkswaterstaat é a organização responsável pela concepsão, construção, gestão e manutenção das principais infra-estruturas da Holanda. Isso inclui a rede rodoviária principal, a rede de hidrovia principal e os sistemas de águas. Faz parte do Ministério dos Transportes, das Obras Públicas e da Gestão da Água.
O referido ministério é o principal encarregado por de cerca de 3500 pontes na Holanda por isso procura-se conhecer a vida útil restante de suas estruturas, já que é sabido que durante a vida útil, as estruturas precisarão ser reparadas. Neste momento, uma estratégia de reparo baseada em inspeções é usada para determinar quando o reparo será feito. Esta estratégia de reparo para pontes de concreto na Holanda resulta em reparação de pontes a cada 25 até 35 anos. O reparo será normalmente de 0.5% até 15.5% da área da estrutura (van Beek et al., 2003).
O banco de dados inclui um total de 5986 eventos de inspeções registradas para 2473 superestruturas individuais. Ignorando o tempo entre a construção da ponte e uma primeira inspeção, há 3513 transições registradas entre estados de condição. Estes dados contém informações sobre o estado em que a estrutura de pontes encontraram-se durante as inspeções, ou seja, contém um histórico de inspeções e o ano de construção.
\[
\begin{array}{c|c}\hline
\mbox{Estado} & \mbox{Classificação} \\[0.8em]\hline
\mbox{Perfeito} & 0 \\
\mbox{Muito bom} & 1 \\
\mbox{Bom} & 2 \\
\mbox{Razoável} & 3 \\
\mbox{Mediocre} & 4 \\
\mbox{Ruim} & 5 \\
\mbox{Muito ruim} & 6 \\\hline
\end{array}
\] Tabela 6.1. Esquema de classificação da
condição das pontes.
Para o gerenciamento das informações utilizam-se diversos sistemas, dois deles: PONTIS e BRIDGIT são dois dos Sistemas de Gerenciamento de Ponte mais comuns atualmente disponíveis. Ambos têm suas origens no Arizona Pavement Management System desenvolvido no final da década de 1970 e são quase exclusivamente utilizados nos Estados Unidos. Todos esses modelos usam Cadeias de Markov para modelar a deterioração incerta das pontes ao longo do tempo.
Nos Países Baixos, os resultados das inspeções de ponte são registrados em um banco de dados, que é usado principalmente para manutenção de registros. Esta base de dados é uma fonte muito rica de informações, contém dados coletados ao longo de quase 20 anos, e a finalidade da pesquisa atual é utilizar esses dados para estimar a taxa de deterioração.Definida a codificação dos estados que vamos utilizar, classificado o estado de conservação das pontes inspecionadas procedemos à estimação da matriz de probabilidades de transição desta cadeia. A tabela a seguir mostra a contagem das transições de cada estado para qualquer outro estado.
\[
\begin{array}{c|ccccccc}\hline
\mbox{De} \backslash \mbox{Para} & 0 & 1 & 2 & 3 & 4
& 5 & 6 \\[0.8em]\hline
0 & 520 & 134 & 327 & 111 & 36 & 7 & 0 \\
1 & 270 & 128 & 222 & 97 & 36 & 7 & 0 \\
2 & 284 & 101 & 368 & 193 & 61 & 9 & 5 \\
3 & 94 & 33 & 119 & 131 & 42 & 3 & 1 \\
4 & 16 & 14 & 42 & 50 & 17 & 7 & 0 \\
5 & 7 & 3 & 4 & 4 & 3 & 0 & 1 \\
6 & 1 & 1 & 0 & 3 & 1 & 0 & 0
\end{array}
\] Tabela 6.2. Contagem original de transições
das classificações das condições das pontes.
Nós vemos que a informação, que vem de dados de deterioração, é bastante subjetiva. Vemos também que a classificação da condição varia de um perfeito (estado 0) para um muito ruim (estado 6), através da definição de estados. Interessante observar que pontes em estados mau e muito mal raramente acontecem, por causa disso procedeu-se à junção destes estados e codificou-se como 5.
Como resultado se obtem como resposta a matriz de probabilidades de transição estimada a seguir: \[ \widehat{\mathcal{P}}=\begin{pmatrix} 0.4581 & 0.1181 & 0.2881 & 0.0978 & 0.0317 & 0.0062 \\ 0.3553 & 0.1684 & 0.2921 & 0.1276 & 0.0474 & 0.0092 \\ 0.2782 & 0.0989 & 0.3604 & 0.1890 & 0.0597 & 0.0138 \\ 0.2222 & 0.0780 & 0.2813 & 0.3097 & 0.0993 & 0.0095 \\ 0.1096 & 0.0959 & 0.2877 & 0.3425 & 0.1164 & 0.0479 \\ 0.2857 & 0.1429 & 0.1429 & 0.2500 & 0.1429 & 0.0356 \end{pmatrix} \]
Foi realizada uma análise sobre modelagem de deterioração de pontes mas os dados foram coletados, no banco de dados denominado DISK e fornecidos pelo Divisão de Engenharia do Rijkswaterstaat. Continham um número limitado de estados de condições das pontes, foi decidido usar uma Cadeia de Markov para modelar a deterioração. O modelo de deterioração de Markov é baseado em condições. Por isso, é flexível na adaptação aos dados de inspeção visual.
Infelizmente, não foi possível observar o tempo exato das transições. Assim, foi adaptada uma Cadeia de Markov com censura de intervalo. Censura de intervalo significa que não sabemos a hora exata de um evento. Em nosso contexto, isso significa que não sabemos a hora em que a ponte se move de um estado para outro. Uma probabilidade de transição foi definida como a probabilidade de uma ponte passar de um estado para outro, igual ou pior. Assumimos que nenhuma manutenção foi realizada entre as inspeções.
Com o auxilio dos seguintes comandos investigamos o comportamento futuro desta cadeia:
Estados = c("0", "1", "2", "3", "4", "5")
Pontes = matrix(c(520, 134, 327, 111, 36, 7,
270, 128, 222, 97, 36, 7,
284, 101, 368, 193, 61, 14,
94, 33, 119, 131, 42, 4,
16, 14, 42, 50, 17, 7,
8, 4, 4, 7, 4, 1),
nrow = 6, ncol = 6, byrow = TRUE, dimnames = list(Estados, Estados))
Pontes = as(as.table(Pontes), "markovchain")
Pontes## Unnamed Markov chain
## A 6 - dimensional discrete Markov Chain defined by the following states:
## 0, 1, 2, 3, 4, 5
## The transition matrix (by rows) is defined as follows:
## 0 1 2 3 4 5
## 0 0.4581498 0.11806167 0.2881057 0.09779736 0.03171806 0.006167401
## 1 0.3552632 0.16842105 0.2921053 0.12763158 0.04736842 0.009210526
## 2 0.2781587 0.09892262 0.3604310 0.18903036 0.05974535 0.013712047
## 3 0.2222222 0.07801418 0.2813239 0.30969267 0.09929078 0.009456265
## 4 0.1095890 0.09589041 0.2876712 0.34246575 0.11643836 0.047945205
## 5 0.2857143 0.14285714 0.1428571 0.25000000 0.14285714 0.035714286steadyStates(Pontes)## 0 1 2 3 4 5
## [1,] 0.3243988 0.1092058 0.3077224 0.1852073 0.06111703 0.01234871Temos por resposta que, aproximadamente, 32.4% das pontes inspecionadas corresponderão a pontes em perfeito estado, 10.9% corresponderão a pontes em estado muito bom, 30.8% permanecerão em estado bom, 18.5% das pontes corresponderão a pontes em estado de conservação razoável, 6.1% estarão em estado medíocre na próxima avaliação enquanto 1.2% delas estarão em estado ruim ou muito ruim.
6.1.2. Intervalos não coincidentes
Consideremos a situação na qual \(L_0\) seja o intervalo de observação, porém o desejado seria \(L_d\) a duração do ciclo desejado. O estimador de máxima verossimilhança da matriz de transições de probabilidades é \(\widehat{\mathcal{P}}_0\), associada com o comprimento do ciclo de observação \(L_0\).
Pela propriedade de invariância, o estimador de máxima verossimilhança da matriz de transição associada com a duração do ciclo \(L_d\) é obtida como \[ \widehat{\mathcal{P}}_d=\widehat{\mathcal{P}}^t, \] onde \(t=L_d/L_0\).
No exemplo anterior, se em vez de observarmos por um período de um ano tivesse sido o período de observação de dois anos: \(L_0=2\) e \(L_d=1\). Então, pode-se encontrar a raiz quadrada da matriz de transição estimada, devido a que \(t=1/2\).
O cálculo da matriz \(\widehat{\mathcal{P}}_d\) é simples a partir da decomposição de \(\widehat{\mathcal{P}}_0\) em seus valores e vectores próprios, chamada de decomposição espectral. Com base nesta decomposição, esta matriz pode ser escrita como \[ \widehat{\mathcal{P}}_0=V\Lambda V^{-1}, \] onde \[ \Lambda=\begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & \lambda_N \end{pmatrix} \] é a matriz de auto-valores e \(V\) a matriz de auto-vetores correspondentes.
Segue então que \[ \widehat{\mathcal{P}}_0^t=V\Lambda^t V^{-1}, \] onde \[ \Lambda^t=\begin{pmatrix} \lambda_1^t & 0 & \cdots & 0 \\ 0 & \lambda_2^t & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & \lambda_N^t \end{pmatrix}\cdot \]
Os autovalores são transformados segundo o valor da potência \(t\), mas os autovetores não mudam. Temos diversas opções disponíveis de funções de decomposição de matrizes, de forma que estes cálculos podem ser feitos muito rapidamente, por exemplo, no R a função básica eigen permite realizar estes cálculos.
Este método é muito semelhante a um método utilizado para se obter o estimador de máxima verossimilhança da matriz de transição a tempo contínuo em Kalbfleisch and Lawless (1985). No entanto, enquanto este método funciona sempre no caso a tempo contínuo, nem sempre funciona no caso a tempo discreto.
Um modelo a tempo discreto não é necessariamente Markov em todos os ciclos. Isto é comparável a dizer que alguns dos valores próprios da matriz de transição podem ser negativos. Desde que a matriz de transição estimada \(\widehat{\mathcal{P}}\) seja semidefinida positiva, todos os valores próprios serão não-negativos, este método permitirá calcular o estimador de máxima verossimilhança diretamente.
O procedimento apresentado não é único na literatura especializada, porém o consideramos de fácil implementação. Um outro procedimento pode ser consultado em Miller and Homan (1994).
Exemplo 6.7: estudo de coorte sobre o HIV
Vejamos o exemplo a seguir o qual é um estudo de coorte sobre o HIV presentado em Sendi et al. (1999). Os pesquisadores construíram uma Cadeia de Markov estacionária para descrever a progressão mensal de indivíduos infectados por HIV em maior risco de desenvolver infecção por Complexo Mycobacterium avium.
O Complexo Mycobacterium avium é um grupo de bactérias que pode ser encontrado normalmente na rede hidráulica das cidades e em pessoas com imunossupressão, como portadores do HIV/AIDS. Esta progressão incluía a possibilidade de movimento entre três faixas de contagem de células CD4 distintas, com e sem AIDS.
Neste modelo, definiram-se quatro categorias amplas de estados de saúde para descrever a história de indivíduos infectados pelo HIV com alta depleção de CD4 e maior risco de desenvolver doença de MAC. Essas categorias são:
não AUXILIA. Indivíduos infectados pelo HIV assintomáticos sem história de doença do indicador de AIDS. Estratificada a condição “Sem AIDS” em três níveis de CD4: 0-49 células/mm\(^3\), 50-74 células/mm\(^3\), 75+ células/mm\(^3\). Condiciona-se o risco de AIDS e MAC sintomáticos na contagem de células CD4 do paciente;
AIDS (não-MAC). Os indivíduos infectados pelo HIV com AIDS de acordo com os Centros de Controle e Prevenção de Doenças dos Estados Unidos (1992) expandiram a definição de vigilância (excluindo uma contagem de células CD4 <200 células/mm\(^3\)). Estratificada a condição de AIDS em três níveis de CD4: 0-49 células/mm\(^3\), 50-74 células/mm\(^3\), 75+ células/mm\(^3\). Condiciona-se o risco de morrer de AIDS sintomático ou de desenvolver MAC após a contagem de células CD4 do paciente;
MAC. Pacientes HIV-positivos com infecção pelo complexo Mycobacterium avium e;
morte. Os pacientes que morrem por qualquer causa mudam para absorver a Cadeia de Markov.
\[ \begin{array}{c|ccc} \mbox{De} \setminus \mbox{Para} & 0-49 & 50-74 & 75+ \\[0.8em]\hline 0-49 & 682 & 33 & 25 \\ 50-74 & 154 & 64 & 47 \\ 75+ & 19 & 19 & 43 \end{array} \] Tabela 6.3. Transições observadas em seis meses na contagem de células CD4 (1993-1995). Estudo realizado na Suiça com pacientes infectados pelo HIV.
Dados coletados num estudo multicêntrico onde os pacientes infectados pelo HIV têm visitas de acompanhamento bastante regulares, a cada seis meses. Os dados estão disponíveis no arquivo de dados craigsendi, pacote markovchain e mostrados na tabela acima.
data(craigsendi)
csMc = as(craigsendi, "markovchain")
csMc## Unnamed Markov chain
## A 3 - dimensional discrete Markov Chain defined by the following states:
## 0-49, 50-74, 75-UP
## The transition matrix (by rows) is defined as follows:
## 0-49 50-74 75-UP
## 0-49 0.9216216 0.04459459 0.03378378
## 50-74 0.5811321 0.24150943 0.17735849
## 75-UP 0.2345679 0.23456790 0.53086420Devemos mencionar que os dados apresentados na tabela acima constituem a contagem das transições observadas mas, uma vez convertidos em Cadeia de Markov temos as probabilidades estimadas de transição. A apresentação do exemplo acerca da estudo de coorte HIV permite-nos a leitura dos dados na referência mencionada antes e guardados no arquivo craigsendi. Esclarecemos novamente que os dados originais foram observados num período de seis meses, o qual não é o desejado no estudo, queremos o comportamento mensal.
Mostramos a continuação nossa implementação no R do procedimento para encontrarmos a matriz de probabilidades de transição no ciclo desejado. Observe que, em nossa implementação no seguinte exemplo, todos oa auto-valores foram positivos.
Para esta análise, a duração do ciclo desejado é de um mês. Para estimar a matriz de probabilidades de transição para esse intervalo, vamos decompor \(\widehat{\mathcal{P}}^6\). Utilizando a função R eigen obtemos:
L = eigen(csMc@transitionMatrix)
L## eigen() decomposition
## $values
## [1] 1.0000000 0.5701572 0.1238380
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.5773503 -0.1276431 0.02818224
## [2,] -0.5773503 0.2866930 -0.87301666
## [3,] -0.5773503 0.9494811 0.48687542Agora vamos transformar esta matriz à situação procurada, ou seja, numa transição mensal.
csMc1 = L$vectors%*%diag((L$values)^(1/6))%*%solve(L$vectors)
csMc1 = new("markovchain", byrow=T, transitionMatrix=csMc1)Traduzindo estes resultados: a matriz de auto-valores é \[ \Lambda=\begin{pmatrix} 1.000000 & 0 & 0 \\ 0 & 0.5701572 & 0 \\ 0 & 0 & 0.123838 \\ \end{pmatrix} \] e a matriz de auto-vetores correspondente é \[ V=\begin{pmatrix} -0.5773503 & -0.1276876 & 0.02817213 \\ -0.5773503 & 0.2867671 & -0.87297660 \\ -0.5773503 & 0.9494527 & 0.48694783 \\ \end{pmatrix}\cdot \]
Como mencionado o ciclo observado foi de 6 meses, mas o ciclo desejado é de um mês. Para isso tomamos a raiz sexta de \(\Lambda\) e fazendo os cálculos sugeridos obtemos que a matriz de transição estimada de um mês é \[ \widehat{\mathcal{P}}=\begin{pmatrix} 0.9819 & 0.0122 & 0.0059 \\ 0.1766 & 0.7517 & 0.0717 \\ 0.0177 & 0.0933 & 0.8830 \\ \end{pmatrix}\cdot \]
Se esta matriz fosse multiplicada seis vezes o resultado será a matriz \(\widehat{\mathcal{P}}^6\) como esperado, esta foi a matriz apresentada obtida Tabela 6.3. Observe que este processo é muito rápido e simples. Neste exemplo, a matriz sugere que haverão demasiados pacientes no estado 0-49 após seis ciclos.
Podemos agora, inclusive, identificarmos a distribuição estacionária, para isso fazemos:
steadyStates(csMc1)## 1 2 3
## [1,] 0.8343668 0.07659214 0.08904103Significa que, a longo prazo, 83.4% dos indivíduos infectados por HIV mantêm-se na faixa 0-49 de contagem de células CD4, 7.7% apresentam contagem na faixa 50-74 em 8.9% dos casos a contagem é 75 ou mais.
6.1.3. Intervalos desiguais
Em muitas situações, os intervalos de observação podem ser desiguais em comprimento. Como exemplo, suponha que uma matriz de transição de um ano seja o desejado, mas o grupo foi observado somente no ano segundo e no terceiro ano. Nesta situação, a matriz de probabilidades de transição de um ano poderia ser estimada utilizando apenas as informações do ano dois a três, mas isso elimina metade dos dados observados. Idealmente, gostariamos de usar toda as informações.
Isto pode ser feito usando o algoritmo EM (Dempster at all., 1977). O passo E imputa os dados em falta (ausentes) calculando o número esperado de transições de ciclo único. O passo M trata o número esperado de transições de ciclo único como o verdadeiro conjunto de dados e maximiza a verosimilhança. Este processo é repetido até que as probabilidades de transição se estabilizem.
Lembremos da situação em que os intervalos de observação e a duração do ciclo coincidem, ou seja, o processo de estimação considerado na Seção V.1.1. Se \(n_{rc}\) representa o nuacute;mero de indivíduos que se movem do estado \(x\) para o estado \(y\) em um ciclo, a função de verossimilhança é \[\begin{equation*} L \, = \, \prod_{x=1}^d \prod_{y=1}^d p_{x,y}^{n_{x,y}}, \end{equation*}\] e o método descrito na Seção V.1.1 fornece o estimador de máxima verosimilhança de \(\mathcal{P}=(p_{x,y})\).
Considere que há \(T\) intervalos de observação que são múltiplos inteiros \(k_1, k_2,\cdots, k_T\) da duração do ciclo. Os dados ausentes são os estados de saúde de cada indivíduo nos ciclos não observados. Assim, o algoritmo EM envolve imputar esses estados aos ciclos não observados, calcular o número esperado de transições e, em seguida, usar o método descrito na Seção V.1.1 de intervalos de observação coincidentes para obter uma nova estimativa da matriz de transição.
Isso é repetido até que a matriz de transição se estabilize. Uma matriz de transição inicial é necessária para iniciar o algoritmo. A convergência para o estimador de máxima verossimilhança não é garantida, pode convergir para máximos locais, portanto, várias matrizes de transição iniciais são recomendadas.
Para a etapa E, a matriz de transição estimada é utilizada para calcular a probabilidade de cada caminho que um sujeito poderia ter seguido para terminar onde parou após os \(k_t\) ciclos. Por exemplo, dada uma matriz de transição \(\mathcal{P}_d\) de um ano, as probabilidades de transição de dois anos são dadas calculando \(\mathcal{P}_d\times \mathcal{P}_d\).
Rotulando a matriz de transição de um ano \[\begin{equation*} \mathcal{P}_d = \begin{pmatrix} p_{1,1} & p_{1,2} & \cdots & p_{1,d} \\ p_{2,1} & p_{2,2} & \cdots & p_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ p_{d,1} & p_{d,2} & \cdots & p_{d,d} \end{pmatrix}, \end{equation*}\] este produto pode ser expresso em termos de probabilidades de transição de um ano como \[\begin{equation*} \mathcal{P}_0 \, = \, \mathcal{P}_d\times \mathcal{P}_d = \begin{pmatrix} \displaystyle \sum_{j=1}^d p_{1,j}p_{j,1} & \displaystyle \sum_{j=1}^d p_{1,j}p_{j,2} & \cdots & \displaystyle \sum_{j=1}^d p_{1,j}p_{j,d} \\ \displaystyle \sum_{j=1}^d p_{2,j}p_{j,1} & \displaystyle \sum_{j=1}^d p_{2,j}p_{j,2} & \cdots & \displaystyle \sum_{j=1}^d p_{2,j}p_{j,d} \\ \vdots & \vdots & \ddots & \vdots \\ \displaystyle \sum_{j=1}^d p_{d,j}p_{j,1} & \displaystyle \sum_{j=1}^d p_{h,j}p_{j,2} & \cdots & \displaystyle \sum_{j=1}^d p_{d,j}p_{j,d} \end{pmatrix}, \end{equation*}\] onde cada produto de probabilidade \((p_{x,j}p_{j,y})\) representa um caminho possível do estado inicial \(x\) para o estado \(y\) após dois ciclos (anos).
Denotemos o número de indivíduos observados movendo-se entre o estado \(x\) e o estado \(y\) após \(k_t\) ciclos como \(n_{x,y}^{k_t}\). Dada a probabilidade de cada caminho possível, o restante desta etapa envolve estimar o número de sujeitos que seguem cada um desses caminhos e então calcular o número de transições de ciclo único. No modelo acima com \(d\) estados, existem \(d\) caminhos em cada célula da matriz de transição de dois anos.
Cada um dos indivíduos daquela célula deve ter seguido um dos \(d\) caminhos. O número esperado de indivíduos para seguir cada caminho é baseado na probabilidade relativa de cada caminho, ou seja, uma distribuição multinomial. Por exemplo, na célula superior esquerda, a probabilidade de um indivíduo seguir o caminho \((1\to 1\to 1)\) é \[\begin{equation*} P(C_n=1, C_{n+1}=1,C_{n+2}=1 | C_0=1, \cdots, C_{n-1}=1) \end{equation*}\] então, o número esperado de assuntos para ter seguiu este caminho é
Algoritmo EM
Neste problema, a chave para o algoritmo EM é um E-passo eficiente. Lembremos que o E-passo envolve (1) calcular a probabilidade de cada caminho possível, (2) obter o número esperado de indivíduos para seguir cada caminho e (3) calcular o número de transições de ciclo único. Nesta seção descrevemos uma abordagem orientada matricialmente para manter o controle de todos os caminhos potenciais.
Mostramos, com um exemplo, que tal matriz que contém todos os caminhos potenciais para \(k=3\) ciclos quando há apenas \(d=2\) estados: \[\begin{equation*} \mathcal{P}_3 = \begin{pmatrix} p_{1,1}p_{1,1}p_{1,1} & p_{1,1}p_{1,1}p_{1,2} \\ p_{1,1}p_{1,2}p_{2,1} & p_{1,1}p_{1,2}p_{2,2} \\ p_{1,2}p_{2,1}p_{1,1} & p_{1,2}p_{2,1}p_{1,2} \\ p_{1,2}p_{2,2}p_{2,1} & p_{1,2}p_{2,2}p_{2,2} \\ p_{2,1}p_{1,1}p_{1,1} & p_{2,1}p_{1,1}p_{1,2} \\ p_{2,1}p_{1,2}p_{2,1} & p_{2,1}p_{1,2}p_{2,2} \\ p_{2,2}p_{2,1}p_{1,1} & p_{2,2}p_{2,1}p_{1,2} \\ p_{2,2}p_{2,2}p_{2,1} & p_{2,2}p_{2,2}p_{2,2} \end{pmatrix} \end{equation*}\]
Observe a organização desses caminhos. A primeira coluna contém os caminhos que terminam no estado 1 e a segunda coluna contém os caminhos que terminam no estado 2. As primeiras quatro linhas contêm caminhos que começam no estado 1 e as últimas quatro linhas contêm caminhos que começam no estado 2. Finalmente, as três transições de ciclo único que constituem um caminho têm um padrão específico à medida que você percorre os caminhos dentro de uma coluna. Este tipo de organização torna a contabilização do E-passo fácil.
Para construir tal matriz, considere uma matriz \(d\times d\) de transição de ciclo único \(\mathcal{P}\) e os dados observados em \(T\) comprimentos de intervalo únicos iguais a \(k_t\), \(t =1,2,\cdots,T\) ciclos. Uma vez que a construção da matriz é semelhante para cada comprimento de ciclo, assuma um único intervalo igual a \(k\) ciclos. A matriz, \(\mathcal{P}_k\), uma matriz \(d^k\times d\), é construída usando a seguinte equação de multiplicação de matriz iterativa, \(\mathcal{P}_1=\mathcal{P}\) e \[\begin{equation*} \mathcal{P}_k(d(r-1)+1,j) = \mathcal{P}_{k-1}(r,c)\times \mathcal{P}(c,j)\left\{ \begin{array}{l} r=1,2,\cdots,d^{k-1} \\ c=1,2,\cdots,d \\ j=1,2,\cdots,d \end{array}\right.\cdot \end{equation*}\]
Em outras palavras, a primeira linha de \(\mathcal{P}_k\) é o produto Kronecker do primeiro elemento em \(\mathcal{P_{k-1}}\) e a primeira linha de \(\mathcal{P}\). A segunda linha é o produto Kronecker do segundo elemento \(\mathcal{P}_{k-1}(1,2)\) e a segunda linha de \(\mathcal{P}\) e assim por diante. A matriz \(\mathcal{P}_k\) tem todos os caminhos potenciais organizados de forma que cada coluna \(c\) contenha todos os caminhos que terminam no estado \(c\) com as primeiras \(d^{k-1}\) linhas contendo os caminhos que começam no estado 1, as próximas \(d^{k-1}\) linhas contendo os caminhos que começam no estado 2 e assim por diante. Isso permite o cálculo do número esperado de sujeitos para seguir cada caminho, uma vez que organiza todos os caminhos possíveis em linhas adjacentes e uma única coluna.
Na construção de cada uma das probabilidades em \(\mathcal{P}_k\), \(k\) elementos únicos de \(\mathcal{P}\) foram multiplicados juntos. Usamos o padrão de multiplicação para registrar as transições de ciclo único. Seja \(\widehat{n}(i,j)\) o número esperado de sujeitos para seguir o caminho descrito na linha \(i\) e coluna \(j\) de \(\mathcal{P}_k\). O número esperado de transi¸ões de ciclo único de \(r\) para \(c\) é \[\begin{equation*} \end{equation*}\]
Exemplo 6.8
Como o número de caminhos potenciais que um sujeito pode seguir cresce rapidamente à medida que \(d\) e \(k_t\) aumentam, este exemplo consiste em um conjunto de dados simulados de um modelo de três estados \(\mathcal{S}=\{1,2,3\}\), sendo o estado 3 um estado absorvente. Temos observações entre o segundo e terceiro ciclos, mostradas na matriz do número de transições abaixo: \[\begin{equation*} (n_{x,y})=\begin{pmatrix} 214 & 45 & 41 \\ 56 & 62 & 82 \\ 0 & 0 & 0 \end{pmatrix}\cdot \end{equation*}\]
Estados = c("1", "2", "3")
Matriz2 = matrix(c(227, 22, 21, 20, 70, 17, 0, 0, 138),
nrow = 3, ncol = 3, byrow = TRUE, dimnames = list(Estados, Estados))
Matriz2 = as(as.table(Matriz2), "markovchain")
Matriz2## Unnamed Markov chain
## A 3 - dimensional discrete Markov Chain defined by the following states:
## 1, 2, 3
## The transition matrix (by rows) is defined as follows:
## 1 2 3
## 1 0.8407407 0.08148148 0.07777778
## 2 0.1869159 0.65420561 0.15887850
## 3 0.0000000 0.00000000 1.00000000Desejamos a matriz de probabilidades de transição de um mês. Este é um problema pequeno o suficiente para que o algoritmo EM possa ser feito usando um pacote de planilha.
As etapas E e M são detalhadas a seguir. As equações da etapa E combinam as transições de um ciclo observadas (primeiro total em cada equação) com o número imputado de transições de um ciclo com base nas transições de dois ciclos observadas. Nesse caso, cada célula rc (r; c42) na matriz de transição de dois ciclos é uma soma de probabilidades de dois caminhos, yr1y1c þ yr2y2c, enquanto a última célula nas duas primeiras linhas é uma soma de probabilidades de três caminhos.6.2. Cadeias de ordem superior
Um processo de Markov é caracterizado por um conjunto (possivelmente infinito) de estados. Presume-se que o sistema que está sendo modelado por esse processo ocupe um e apenas um desses estados em um determinado momento. A evolução do sistema é representada por transições do processo de Markov de um estado para outro. Supõe-se que essas transições ocorram instantaneamente; em outras palavras, a ação de passar de um estado para outro consome tempo zero. A propriedade fundamental de um sistema Markoviano, conhecida como propriedade Markoviana, é que a evolução futura do sistema depende apenas do estado atual do sistema e não da história passada.
Muitas vezes acontece ser útil descrever um processo estocástico como um conjunto de estados discretos com transições probabilísticas e exemplos abundam em diversos campos, como o estudo de processos químicos, sequências de DNA, finanças dentre outros. Se a probabilidade de transição para o próximo estado é condicionada apenas no estado atual, chamamos este modelo de uma Cadeia de Markov e quando as probabilidades condicionais não são dadas de outra forma, elas serão estimadas a partir de uma série temporal de observações.
Mas, e caso a probabilidade de transição para o próximo estado seja condicionada não somente no estado atual? em tais situações surgem novos questionamentos. Se a ordem da Cadeia de Markov estiver em questão, o primeiro é respondermos: o que é ordem de uma Cadeia de Markov?
Definição 6.1: Ordem de uma Cadeia de Markov.
Uma sequência de observações discretas a tempo discreto \(\{C_n\}_{n\geq 1}\) formam uma Cadeia de Markov de ordem \(k\) se a probabilidade condicional satisfaz \[ P(C_{n+1} \, | \, C_n,C_{n-1},\cdots) = P(C_{n+1} \, | \, C_n,\cdots,C_{n-k+1}), \qquad \forall \, k\leq n\cdot \]
As Cadeias de Markov consideradas até o momento são cadeias de ordem um, ou seja, \(k=1\). Isso significa, como sabemos, que as probabilidades de transição para um estado futuro dependem apenas do estado atual e não de estados anteriores. Um processo de ordem \(k\) pode sempre ser lançado como de primeira ordem agrupando estados. Um processo que não tenha dependência do passado ou presente, como variáveis aleatórias independentes, é dito ser uma Cadeia de Markov de ordem zero.
Por outro lado, como deve ser facilmente percebido, cadeias de ordens superiores, ou seja, cadeias de segunda ordem ou superiores implicam numa representação mais complicada.
Considere primeiro uma Cadeia de Markov de segunda ordem. Dado que um indivíduo está no estado \(z\) no instante \(n-2\) e em \(y\) no instante \(n-1\), seja \(\gamma_{zyx}\) a probabilidade de estar o indivíduo no estado \(x\) no instante \(n\). Uma cadeia estacionária de primeira ordem é uma cadeia especial de segunda ordem, na qual \(\gamma_{zyz}\) não depende de \(z\).
Para vermos isso, considere o par de estados sucessivos \(z\) e \(y\) definidos como um estado composto \((z,y)\). A probabilidade do estado composto \((y,x)\) no instante \(n\) dado o estado composto \((z,y)\) no instante \(n-1\) é \(\gamma_{zyx}\). Vejamos isso.
Sabemos que \[ P(X_n=x \, | \, X_{n-1}=y,X_{n-2}=z) = \gamma_{zyx} \] e queremos verificar se \(P\big(X_n=(y,x) \, | \, X_{n-1}=(z,y)\big)= \gamma_{zyx}\). Logo \[ \begin{array}{rcl} P\big(X_n=(y,x) \, | \, X_{n-1}=(z,y)\big) & = & \gamma_{zyx} \\[0.8em] P\big(X_n=(y,x) \, | \, X_{n-1}=y,X_{n-2}=z\big) & = & P\big(X_n=x, X_{n-1}=y \, | \, X_{n-1}=y,X_{n-2}=z\big) \, = \, \gamma_{zyx}\cdot \end{array} \]
Claro que, a probabilidade do estado \((w,x)\), \(w\neq y\), dado \((x,y)\), é zero. Os estados compostos podem ser encontrados para formar uma cadeia com \(d^2\) estados, onde \(d\) é o número de estados e com certas probabilidades de transição 0. Esta repressentação nos ajudará na descrição dos testes de verificação da ordem de uma cadeia a serem descritos aqui.
6.2.1. Testes
Verificarmos a ordem de uma Cadeia de Markov poderá ser realizado de diversas maneiras, mas aqui consideraremos somente uma delas. Conhecido como teste aproximado, será descrito nesta seção e é baseado na estatística \(\chi^2\). Devemos mencionar que muitos dos resultados apresentados aqui foram resumidos no artigo de Anderson and Goodman (1957).
Modelos de cadeia de Markov são usados em vários campos aplicados, como análise de séries temporais, estudos longitudinais, dados de tempo de vida real e problemas ambientais. O comportamento de uma cadeia de Markov depende da matriz de transição, que contém probabilidades de transição. Na maioria dos estudos práticos, a matriz de transição é desconhecida e precisa ser estimada. Existem vários métodos para estimativa e procedimento de teste de probabilidades de transição. No entanto, a maioria dos pesquisadores trabalhou na estimativa de parâmetros. No entanto, relatórios sobre procedimentos de teste são dificilmente encontrados. Um dos testes mais importantes em modelos de cadeia de Markov é a estacionariedade das probabilidades de transição, que é interessante trabalhar. Nesta seção, um breve resumo dos testes de procedimento em cadeia de Markov é apresentado.
Anderson and Goodman (1957) obtiveram estimativas de máxima verossimilhança e sua distribuição assintótica em uma cadeia de Markov de ordem arbitrária quando há observações repetidas da cadeia. Testes de razão de verossimilhança e testes \(\chi^2\) são considerados para testar a estacionariedade e a ordem de cadeias de Markov de ordem superior.
Billingsley (1961b) utilizou a fórmula de Whittle, os métodos do qui-quadrado e da máxima verossimilhança para estimar e testar os parâmetros. Uma amostra \(\{x_1,x_2,\cdots,x_n\}\) de um processo de Markov de primeira ordem com probabilidades de transição \(p_{ij}\) e probabilidades iniciais \(p_i\) foi considerada. Se a matriz \(d\times d\) \((n_{x,y})\) for definida como a contagem de transições da sequência, então pode-se demonstrar que \[ \sum_{x,y} \dfrac{(n_{x,y}-n_x p_{x,y})^2}{n_x p_{x,y}} \] é assintoticamente qui-quadrado em distribuição com \(d(d-1)\) graus de liberdade.
Esta estatística qui-quadrado é útil para testar se as probabilidades de transição do processo têm valores especificados \(p_{x,y}\). Surge então o problema natural de testar se essas probabilidades de transição têm uma forma especificada \(p_{x,y}(t)\), onde \(t\) é uma incógnita da amostra.
Bartlett (1951) construiu um teste de razão de verossimilhança para a qualidade do ajuste, comprovando que a distribuição assintótica em cadeias de Markov era normal. Para testar se uma sequência de observações é, no máximo, dependente de r, assume-se que as probabilidades de transição são conhecidas ou, pelo menos, dependem de um número limitado de parâmetros que podem ser estimados.
Se as probabilidades de transição forem completamente desconhecidas, um teste diferente é necessário, e este teste é apresentado (Hoel, 1954), onde a derivação dependeu dos resultados e métodos de Bartlett. No entanto, os métodos anteriores de teste de parâmetros baseavam-se em probabilidades de transição e o teste estatístico dependia de probabilidades de transição.
Nas últimas décadas, mais pesquisas sobre procedimentos de estimação e teste de parâmetros do modelo de cadeia de Markov foram estendidas a novos métodos nos quais covariáveis e funções de ligação são utilizadas e medidas repetidas são consideradas. Por exemplo, Muenz e Rubinstein (1985) propuseram um modelo para cadeia de Markov baseado em covariáveis e mostraram como as covariáveis se relacionam com mudanças de estado.
Uma extensa análise dependente de covariáveis para modelos de Markov de ordem superior foi aprimorada (Islam e Chowdhury, 2006). A influência de covariáveis dependentes do tempo na distribuição marginal da resposta binária foi estudada (Azzalini, 1994). Foi demonstrado que as covariáveis se relacionam apenas com o valor médio do processo, independentemente do parâmetro de associação. Uma aplicação de modelos de Markov baseados em distribuição marginal é fornecida (Shafiqur Rahman e Islam, 2007).
Um teste de qualidade de ajuste para o modelo de regressão logística baseado em dados binários foi empregado (Tsiatis, 1980). Ele modificou o modelo relacionado à probabilidade de respostas com um conjunto de covariáveis.
Exemplificaremos a teoria a ser apresentada neste seção com o seguinte exemplo, inspirado no trabalho de Doubleday & Esunge (2011). A ideia é usar Cadeias de Markov para prever o comportamento dos preços das ações utilizando o índice Dow Jones Industrial Average (DJIA); este foi um índice criado em 1896 por Charles Dow, editor do The Wall Street Journal e fundador do Dow Jones & Company.
O DJIA é ao lado do Nasdaq Composite e do Standard & Poor’s 500 um dos principais indicadores dos movimentos do mercado americano. Dos três indicadores, o DJIA é o mais largamente publicado e discutido. O cálculo deste índice é bastante simples e baseia-se na cotação das ações de 30 das maiores e mais importantes empresas dos Estados Unidos. Como o índice não é calculado pela Bolsa de Valores de Nova Iorque, seus componentes são escolhidos pelos editores do jornal financeiro norte-americano The Wall Street Journal. Não existindo nenhum critério pré-determinado a não ser que os componentes sejam companhias norte-americanas líderes em seus segmentos de mercado.Exemplo 6.9: Tendência de mercado financeiro.
A modelagem do Índice Dow Jones Industrial Average ou DJIA é frequentemente utilizada para determinar estratégias de negociação com o máximo de recompensa. As mudanças no comportamento do DJIA são importantes, pois os movimentos podem afetar profundamente as escolhas dos investidores, sejam estes indivíduos ou corporações.
O objetivo neste exemplo é mostrar como analisar o DJIA usando um modelo estocástico de tempo discreto, ou seja, uma Cadeia de Markov. Dois modelos foram destacados, onde o DJIA foi considerado como sendo em
ganho ou perda e
ganho ou perda pequeno, moderado ou grande.
Esses modelos serão usados para obter probabilidades de transição e a distribuição estacionária.
Os preços de fechamento do mercado são considerados para que a análise possa ser feita de forma discreta e as probabilidades de transição são utilizadas como partes de Cadeias de Markov para modelar o mercado. Dada esta formulação de uma matriz de transição e seu estado estacionário, podemos configurar um sistema de classificação do Dow Jones Industrial Average (DJIA).
A ideia de usar Cadeias de Markov para prever o comportamento dos preços das ações é popular, pois os investidores potenciais estão interessados nas tendências do mercado, o que pode levar a uma estratégia de investimento ideal.
Os valores de fechamento do DJIA foram reunidos para os 252 dias de negociação entre 27 de dezembro de 2016 e 26 de dezembro de 2017. Os dados, apresentados abaixo, foram obtidos do Yahoo! Finance, código “^DJI”. Para a aplicação da estratégia (i), cada dia foi classificado como tendo fechado maior ou menor que o dia anterior, assim permitindo a classificação de dois estados, a saber:
(Estado -1) O valor de fechamento é inferior ao valor de fechamento do dia anterior.
(Estado 1) O valor de fechamento é maior ou igual ao valor de fechamento do dia anterior.
quantmod::getSymbols("^DJI", src = "yahoo", from = "2016-12-27", to ="2017-12-26")## [1] "DJI"head(DJI)## DJI.Open DJI.High DJI.Low DJI.Close DJI.Volume DJI.Adjusted
## 2016-12-27 19943.46 19980.24 19939.80 19945.04 158540000 19945.04
## 2016-12-28 19964.31 19981.11 19827.31 19833.68 188350000 19833.68
## 2016-12-29 19835.46 19878.44 19788.94 19819.78 172040000 19819.78
## 2016-12-30 19833.17 19852.55 19718.67 19762.60 271910000 19762.60
## 2017-01-03 19872.86 19938.53 19775.93 19881.76 339180000 19881.76
## 2017-01-04 19890.94 19956.14 19878.83 19942.16 280010000 19942.16plot(DJI$DJI.Close, main = "Valor de fechamento ")
grid()
Questões em aberto: como vamos construir uma Cadeia de Markov a partir dos dados relatados? Qual a ordem desta cadeia?
Com os comandos a continuação construímos uma Cadeia de Markov com os estados -1 e 1 como descritos acima.
diffY = as.numeric(sign(diff(DJI$DJI.Close, lag = 1)))[-1]
head(diffY)## [1] -1 -1 -1 1 1 -1Agora procedemos à estimação da matriz de transição e obemos por resposta uma lista contendo a estimativa, o logaritmo da verossimilhança e, quando o método “bootstrap” é usado, uma matriz de desvios-padrão e as amostras bootstrap. Quando os métodos “mle”, “bootstrap” ou “map” são usados, os limites de confiança inferior e superior são retornados juntamente com o erro-padrão. O método “map” também retorna o valor esperado dos parâmetros em relação à distribuição a posteriori.
library(markovchain)
createSequenceMatrix(diffY, sanitize=FALSE)## -1 1
## -1 50 57
## 1 57 85mcFitMLE = markovchainFit(data = diffY, method = "map")
mcFitMLE## $estimate
## Bayesian Fit
## A 2 - dimensional discrete Markov Chain defined by the following states:
## -1, 1
## The transition matrix (by rows) is defined as follows:
## -1 1
## -1 0.4672897 0.5327103
## 1 0.4014085 0.5985915
##
##
## $expectedValue
## -1 1
## -1 0.4678899 0.5321101
## 1 0.4027778 0.5972222
##
## $standardError
## [,1] [,2]
## [1,] 0.04757472 0.04757472
## [2,] 0.04073022 0.04073022
##
## $confidenceInterval
## $confidenceInterval$confidenceLevel
## [1] 0.95
##
## $confidenceInterval$lowerEndpointMatrix
## [,1] [,2]
## [1,] 0.3725936 0.4410750
## [2,] 0.3173518 0.5218625
##
## $confidenceInterval$upperEndpointMatrix
## [,1] [,2]
## [1,] 0.5589250 0.6274064
## [2,] 0.4781375 0.6826482
##
##
## $logLikelihood
## [1] -169.5858mcFitMLE$estimate^100## Bayesian Fit^100
## A 2 - dimensional discrete Markov Chain defined by the following states:
## -1, 1
## The transition matrix (by rows) is defined as follows:
## -1 1
## -1 0.4297189 0.5702811
## 1 0.4297189 0.5702811steadyStates(mcFitMLE$estimate)## -1 1
## [1,] 0.4297189 0.5702811Ainda temos também que a distribuição estacionária é \(\pi=(0.432,0.568)\). Significa que temos 43% de probabilidade de perda na nossa carteira de ações com o índice DJAI e 57% (aproximadamente) de probabilidade de ganho. Aplicamos agora os conhecimentos desenvolvidos nesta seção para verificar a ordem da Cadeia de Markov.
Para isso devemos verificar se a cadeia em questão é de ordem um ou não e, como vimos, podemos utilizar a estatística de teste \(\chi^2\). Para isto, ou seja, para aplicarmos o teste aproximado descrito fazemos:
assessOrder(diffY)## The assessOrder test statistic is: 3.909447
## The Chi-Square d.f. are: 2
## The p-value is: 0.1416036com o qual concluímos que aceitamos a hipóteses nula da cadeia com dois estados ser de ordem um.
Para a aplicação da estratégia (ii), cada dia foi classificado como tendo fechado maior, igual ou menor que o dia anterior, assim permitindo a classificação de dois estados, a saber:
(Estado -1) O valor de fechamento é inferior ao valor de fechamento do dia anterior.
(Estado 0) O valor de fechamento é igual estatisticamente ao valor de fechamento do dia anterior.
(Estado 1) O valor de fechamento é maior ou igual ao valor de fechamento do dia anterior.
Para construir estes estes estados vamos considerar os percentis das diferenças:
difer = diff(as.numeric(DJI$DJI.Close, lag = 1))[-1]
summary(difer)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -372.82 -30.95 13.01 19.76 69.17 331.67q1 = summary(difer)[2]
q3 = summary(difer)[4]
Estados = rep(0,length(difer))
Estados <- ifelse(difer > q3, 1,
ifelse(difer < q1, -1, 0))
createSequenceMatrix(Estados, sanitize=FALSE)## -1 0 1
## -1 13 18 31
## 0 22 20 25
## 1 27 29 63mcFitMLE = markovchainFit(data=Estados, method = "mle")
mcFitMLE## $estimate
## MLE Fit
## A 3 - dimensional discrete Markov Chain defined by the following states:
## -1, 0, 1
## The transition matrix (by rows) is defined as follows:
## -1 0 1
## -1 0.2096774 0.2903226 0.5000000
## 0 0.3283582 0.2985075 0.3731343
## 1 0.2268908 0.2436975 0.5294118
##
##
## $standardError
## -1 0 1
## -1 0.05815405 0.06842969 0.08980265
## 0 0.07000621 0.06674830 0.07462687
## 1 0.04366515 0.04525349 0.06669961
##
## $confidenceLevel
## [1] 0.95
##
## $lowerEndpointMatrix
## -1 0 1
## -1 0.09569755 0.1562028 0.3239900
## 0 0.19114854 0.1676832 0.2268683
## 1 0.14130862 0.1550023 0.3986829
##
## $upperEndpointMatrix
## -1 0 1
## -1 0.3236573 0.4244423 0.6760100
## 0 0.4655679 0.4293317 0.5194003
## 1 0.3124729 0.3323927 0.6601406
##
## $logLikelihood
## [1] -258.4417mcFitMLE$estimate^100## MLE Fit^100
## A 3 - dimensional discrete Markov Chain defined by the following states:
## -1, 0, 1
## The transition matrix (by rows) is defined as follows:
## -1 0 1
## -1 0.25 0.2701613 0.4798387
## 0 0.25 0.2701613 0.4798387
## 1 0.25 0.2701613 0.4798387steadyStates(mcFitMLE$estimate)## -1 0 1
## [1,] 0.25 0.2701613 0.4798387Ainda temos também que a distribuição estacionária é \(\pi=(0.25,0.27,0.48)\). Significa que temos 48% de probabilidade de ganho na nossa carteira de ações com o índice DJAI. Aplicamos novamente os conhecimentos desenvolvidos nesta seção para verificar a ordem da Cadeia de Markov.
Para isso devemos verificar se a cadeia em questão é de ordem um ou não e, como vimos, podemos utilizar a estatística de teste \(\chi^2\). Para isto, ou seja, para aplicarmos o teste aproximado descrito fazemos:
assessOrder(Estados)## The assessOrder test statistic is: 7.790757
## The Chi-Square d.f. are: 12
## The p-value is: 0.8012608com o qual concluímos que aceitamos a hipóteses nula da cadeia com três estados ser de ordem um.
Por conveniência, rotularemos os estados que cada medida pode tomar por inteiros positivos até \(d\). Uma seqüência de medidas discretas pode vir de um processo que seja naturalmente discreto, como uma seqüência de DNA ou de um processo contínuo discretizado por um dispositivo de medição analógico-digital.
A abordagem clássica consiste em formular a questão como um teste de hipótese de que uma cadeia tem uma determinada ordem e especificar uma estatística de teste discriminante. Muitas vezes, a estatística de teste escolhida tem uma distribuição limitante conhecida, como \(\chi^2\), da qual o \(p\)-valor é calculado. O uso de uma distribuição limitante é inexato porque somente será alcançado para conjuntos de dados infinitos. Para pequenos conjuntos de dados, pode ser uma aproximação muito fraca.
Existem abordagens alternativas que fornecem uma resposta positiva à seleção de modelos. O Critério de Informação de Akaike (AIC) e o Critério de Informação Bayesiano (BIC) e outros produzem rankings em várias ordens com base em suas verossimilitudes relativas e têm penalidades internas por excesso de ajuste. Infelizmente, essas abordagens dependem de aproximações que somente são válidas no limite de grandes amostras. Em pequenas situações de amostra, não se pode ter certeza de sua eficácia.
Apesar das deficiências dos testes de hipóteses, existe uma vantagem em que é possível, em princípio, realizar um teste de significância de ordem Markov exato que seja válido para qualquer tamanho de amostra. Em vez de depender de propriedades assintóticas, a distribuição da estatística de teste é descoberta por amostragem a partir do conjunto de sequências, aqui referidas como substitutos, que coincidem exatamente com as propriedades da ordem \(n\) da séerie temporal observada. O desafio é gerar eficientemente um grande número de amostras, especialmente para pedidos superiores. Para o nosso conhecimento, nenhuma solução para este problema foi relatada na literatura conhecidas por nós.
6.2.2. Estimação
Sequências de dados ou séries temporais ocorrem com frequência em muitas aplicações do mundo real. Uma das etapas mais importantes na análise de uma sequência de dados (ou série temporal) é a seleção de um modelo matemático apropriado para os dados. Isso ocorre porque ele auxilia em previsões, testes de hipóteses e descoberta de regras.
Uma sequência de dados \(\{C_n\}\), com espaço de estados \(S=\{1,2,\cdots,d\}\), pode ser logicamente representada como um vetor onde \(d\) é finito e cada elemento seja \(C_n=(1,0,\cdots,0)^\top\), se a cadeia o instante \(n\) assumir valores de ordem 1, \(C_n=(0,\cdots,1,\cdots,d)^\top\) se assumir valores de ordem i, ou seja, uma sequência possui \(d\) categorias ou estados possíveis. O modelo convencional para uma cadeia de Markov de ordem \(k\) possui \((d-1)d^k\) parâmetros do modelo.
O principal problema em usar esse tipo de modelo é que o número de parâmetros (as probabilidades de transição) aumenta exponencialmente em relação à ordem do modelo. Esse grande número de parâmetros desencoraja o uso direto de uma cadeia de Markov de ordem superior.
Em Raftery (1985) propôs um modelo de Cadeia de Markov de ordem superior que envolve apenas um parâmetro adicional para cada defasagem extra. O modelo pode ser escrito da seguinte forma: \[ \tag{6.1} P(C_n = j_0 \, | \, C_{n-1}=j_1,\cdots,C_{n-k}=j_k)=\sum_{i=1}^k \lambda_i q_{j_{_0},j_{_i}}, \] onde \[ \sum_{i=1}^k \lambda_i=1 \] e \(Q=(q_{j_{_0},j_{_i}})\) é uma matriz de transição com somas de colunas iguais a um, tal que \[ \tag{6.2} 0\leq \sum_{i=1}^k \lambda_i q_{j_{_0},j_{_i}} \leq 1, \qquad q_{j_{_0},j_{_i}}\in S\cdot \]
A restrição em (6.2) é garantir que o lado direito de (6.1) seja uma distribuição de probabilidade. O número total de parâmetros independentes neste modelo é de tamanho \(k+d^2\). Raftery (1985) provou que (6.1) é análogo ao modelo \(AR(n)\) padrão no sentido de que cada defasagem adicional, após a primeira, é especificada por um único parâmetro e que as autocorrelações satisfazem um sistema de equações lineares semelhante às equações de Yule-Walker.
Além disso, os parâmetros \(q_{j_{_0},j_{_i}}\) e \(\lambda_i\) podem ser estimados numericamente maximizando a log-verossimilhança de (6.1) sujeita às restrições (6.2). No entanto, essa abordagem envolve a resolução de um problema de otimização altamente não linear. O método numérico proposto não garante nem convergência nem um máximo global, mesmo que convirja.
6.2.3. Novo modelo
Nesta subseção mostramos a extensão do modelo de Raftery (1985) porposta por Ching, Ng, and Fung (2008) para um modelo de Cadeia de Markov de ordem superior mais geral, permitindo que \(Q\) varie com diferentes defasagens. Aqui, assumimos que o peso \(\lambda_i\) é não negativo, de modo que \[ \tag{6.3} \sum_{i=1}^k \lambda_i=1\cdot \]
Deve-se notar que (6.1) pode ser reescrito como \[ \tag{6.4} C_{n+k+1}=\sum_{i=1}^k \lambda_i Q C_{n+k+1-i} \] onde \(Q C_{n+k+1-i}\) é a distribuição de probabilidade dos estados no instante \(n+k+1-i\).
Usando (6.3) e o fato de que \(Q\) é uma matriz de probabilidade de transição, notamos que cada entrada de \(C_{n+k+1}\) está entre 0 e 1, e que a soma de todas as entradas também é igual a 1.
O modelo de Raftery em (6.4) pode ser generalizado da seguinte forma: \[ \tag{6.5} C_{n+k+1}=\sum_{i=1}^k \lambda_i Q_i C_{n+k+1-i}\cdot \]
O número total de parâmetros independentes no novo modelo é \((k+kd^2\). Notamos que se \[ Q_1=Q_2 =\dots = Q_k, \] então (6.5) é apenas o modelo de Raftery em (6.4).
Para encontrar as estimativas de \(Q_1,Q_2,\cdots_k\) e \(\lambda_1,\lambda_2,\cdots,\lambda_k\) resolvemos o seguinte problema de minimização: \[ \tag{6.6} \min_\lambda \left\{ || \sum_{i=1}^k \lambda_i \widehat{Q}_i \widehat{C}-\widehat{C}||\right\} \] restrito à \[ \sum_{i=1}^k \lambda_i=1 \qquad \mbox{e} \qquad \lambda_i\geq 0, \; \forall i\cdot \] A solução do sistema (6.6) está programada na função fitMarkovChain no pacote de funções clickstream.
Exemplo 6.10.
Cadeias de Markov são geralmente utilizadas na modelagem de muitos problemas práticos. Elas também são eficazes na modelagem de séries temporais. Neste exemplo, aplicamos o modelo de Cadeias de Markov para analisar e prever séries temporais.
Algumas séries podem ser expressas por uma Cadeia de Markov de tempo discreto de primeira ordem e outras devem ser expressas por um modelo de Cadeia de Markov de ordem superior. Os resultados mostram que o desempenho e a eficácia do modelo de cadeia de Markov para prever séries temporais são muito bons.
Liu (2010) aplicou o modelo de Cadeias de Markov para analisar e prever preço e volume de vendas de carne bovina em um supermercado na cidade de AnKang. Ankang é uma cidade da China, localizada na província de Shaanxi.
price.series = c(5,5,5,5,4,5,3,5,3,3,4,2,5,5,3,1,1,1,3,3,4,1,5,1,1,
3,3,2,5,1,5,1,5,5,5,5,2,1,4,1,1,1,2,4,5,5,1,4,2,4,1,3,4,2,2,5,2,2,5,5,
2,5,4,4,4,2,2,5,2,2,5,5,5,5,3,2,2,5,4,5,2,4,5,5,4,1,1,1,2,
2,3,2,4,5,5,5,2,5,2,5,5,2,4,2,5,5,2,5,5,1,2,3,4,3,3,1,3,1,4,3,
5,4,5,5,4,5,5,2,5,2,5,2,2,3,5,5,3,5,2,5,4,2,1,5,2,5,2,2,2,2,5,
5,4,5,5,2,2,5,2,2,2,3,4,4,4,5,4,5,1,5,5,1,3,5,5,5,1,5,2,2,
2,5,5,5,5,2,4,5,2,2,5,2,5,2,2,2,4,2,2,2,4,5,5,5,3,2,2,5,2,
5,4,4,4,5,3,3,5,3,1,1,4,2,2,5,5,2,5,5,5,2,5,5,3,5,5,4,1,5,5,1,
5,5,1,5,1,5,5,5,5,1,5,5,5,2,1,2,5,2,5,5,2,3,5,5,5,5,5,2,5)O preço semanal da carne bovina em price.series pode ser classificado em cinco estados possíveis \(\{1, 2, 3, 4, 5\}\). A séries de preço é expressa como 1 = muito baixo (≤ RMB 29/kg), 2 = baixo (RMB 29∼32/kg), 3 = médio (RMB 32∼35/kg), 4 = alto (RMB 35∼38/kg), 5 = muito alto (≥RMB 38/kg).
RMB significa Renminbi ou “Moeda do Povo” em mandarim. É a moeda oficial da República Popular da China, e sua unidade principal é o yuan. O termo “Renminbi” refere-se a todo o sistema monetário, enquanto “yuan” é a unidade específica usada para precificar bens e serviços. O RMB é emitido pelo Banco Popular da China e é amplamente utilizado no comércio e nas finanças internacionais.
library(markovchain)
mcAjusteMLE = markovchainFit(price.series)
mcAjusteMLE$estimate## MLE Fit
## A 5 - dimensional discrete Markov Chain defined by the following states:
## 1, 2, 3, 4, 5
## The transition matrix (by rows) is defined as follows:
## 1 2 3 4 5
## 1 0.25000000 0.1250000 0.15625000 0.1250000 0.3437500
## 2 0.04477612 0.3283582 0.07462687 0.1194030 0.4328358
## 3 0.15384615 0.1538462 0.19230769 0.1923077 0.3076923
## 4 0.14705882 0.2352941 0.05882353 0.1764706 0.3823529
## 5 0.11009174 0.2660550 0.08256881 0.1009174 0.4403670assessOrder(price.series)## The assessOrder test statistic is: 82.41408
## The Chi-Square d.f. are: 80
## The p-value is: 0.4045573Com a série de preços price.series, o preço de hoje depende principalmente do preço de ontem. Escolhemos o modelo de cadeia de Markov de primeira ordem. Primeiro, estimamos a matriz de probabilidade de transição de uma etapa \(\mathcal{P}\) e depois encontramos a distribuição estacionária.
steadyStates(mcAjusteMLE$estimate)## 1 2 3 4 5
## [1,] 0.119403 0.25 0.09701493 0.1268657 0.40671641/steadyStates(mcAjusteMLE$estimate)## 1 2 3 4 5
## [1,] 8.375 4 10.30769 7.882353 2.458716Isto significa que em 40% dos casos, o preço é muito alto, estado 5; 25% do casos o preço é baixo, estado 2 e as outras situações ou estados da cadeia acontecem em aproximadamente 10% das situações (estados 1, 3 e 4). Ainda calculamos o tempo de retorno médio a cada estado a partir de cada estado, percebendo que todos os estados são recorrentes positivos.
O inverso da distribuição estacionária observamos que se o preço numa semana estiver muito alto (estado 5), retornará a esse estado em aproximadamente 2.5 semanas, ou seja, se a cadeia estiver no estado 5 demorará 2.5 semanas para retornar ao mesmo estado 5. Caso observe-se o estado 2, demorará 4 semanas em retornar e nas outras situações a demora em retornar a esses estados é de entre 8 a 10 semanas.
As vendas semanais da carne bovina em sales.series pode ser classificado em cinco estados possíveis \(\{1, 2, 3, 4, 5\}\). A série de volume de vendas são expressas como 1 = muito baixo (≤ 50/kg), 2 = baixo (50∼55/kg), 3 = médio (55∼60/kg), 4 = alto (60~65/kg), 5 = muito alto (≥ 65/kg).
sales.series = c(5,1,1,1,1,2,2,3,3,3,4,3,2,2,5,1,5,5,1,2,3,3,2,3,3,2,2,
1,1,5,2,3,3,3,3,2,3,1,2,1,1,5,2,2,5,5,2,3,4,3,4,2,2,1,5,5,1,5,1,5,
5,1,1,5,3,3,3,3,3,3,4,3,3,5,1,5,5,1,1,5,1,5,5,1,5,5,1,5,5,
1,5,1,5,2,3,3,3,3,3,3,3,5,1,5,5,1,5,1,5,5,5,5,1,5,1,1,5,3,3,3,
3,3,3,5,2,2,5,1,1,1,5,1,5,1,1,1,1,1,1,1,1,1,1,1,1,2,5,3,4,4,4,
4,1,3,5,5,1,5,5,1,1,5,1,1,1,1,5,1,2,1,1,2,5,2,1,1,2,3,3,3,
3,4,4,3,2,2,5,1,5,1,1,1,5,2,1,1,1,1,4,4,3,3,3,3,2,5,1,5,5,1,
5,1,5,1,1,2,2,3,3,4,3,3,1,1,1,2,1,1,5,1,1,1,5,1,1,1,1,1,1,1,
2,3,3,1,1,4,3,1,3,2,1,1,1,1,1,5,5,1,5,1,5,1,1,1,1,1,1,1)Essas expressões codificadas, tanto do preço quanto do volume de vendas, são úteis tanto do ponto de vista de marketing quanto do planejamento de produção.
markovchain::assessOrder(sales.series)## The assessOrder test statistic is: NaN
## The Chi-Square d.f. are: 80
## The p-value is: NaNMas as séries de volume de vendas são muito mais complexas. Escolhemos arbitrariamente a ordem cinco, ou seja, \(k = 5\). Primeiramente, estimamos todas as matrizes de probabilidade de transição \(\mathcal{P}_i\) e também temos as estimativas das distribuições de probabilidade estacionárias do produto.
library(clickstream)
sl <- list(sales.series)
class(sl) <- "clikcstream"
mmcAjusteMLE <- fitMarkovChain(clickstreamList = sl, order = 5)
mmcAjusteMLE@lambda## [1] 0.8874696 0.0000000 0.0000000 0.0000000 0.1097973mmcAjusteMLE@transitions$`1`## 1 2 3 4 5
## 1 0.53398058 0.2424242 0.07272727 0.07142857 0.5806452
## 2 0.09708738 0.2424242 0.12727273 0.07142857 0.1129032
## 3 0.01941748 0.3030303 0.60000000 0.50000000 0.0483871
## 4 0.01941748 0.0000000 0.12727273 0.35714286 0.0000000
## 5 0.33009709 0.2121212 0.07272727 0.00000000 0.2580645mmcAjusteMLE@transitions$`5`## 1 2 3 4 5
## 1 0.45454545 0.30303030 0.29090909 0.42857143 0.37096774
## 2 0.07070707 0.18181818 0.20000000 0.07142857 0.12903226
## 3 0.24242424 0.24242424 0.16363636 0.14285714 0.19354839
## 4 0.04040404 0.18181818 0.03636364 0.00000000 0.03225806
## 5 0.19191919 0.09090909 0.30909091 0.35714286 0.27419355matriz = (mmcAjusteMLE@lambda[1]*as.matrix(t(mmcAjusteMLE@transitions$`1`))+
mmcAjusteMLE@lambda[5]*as.matrix(t(mmcAjusteMLE@transitions$`5`)))
matriz[,5] = matriz[,5]+(1-sum(matriz[1,]))
matriz.mmc <- new("markovchain", transitionMatrix = matriz)
markovchain::absorbingStates(matriz.mmc)## character(0)markovchain::recurrentStates(matriz.mmc)## [1] "1" "2" "3" "4" "5"markovchain::steadyStates(matriz.mmc)## 1 2 3 4 5
## [1,] 0.3879658 0.1232851 0.2053124 0.05226081 0.2311758Significando que quase 39% das semanas, as vendas estão num volume muito baixo (estado 1), 23% das semanas o volume de vendas está no estado 5, volume muito alto de vendas, 20.5% o volume de vendas é médio (estado 3), 12% as vendas estão no volume baixo e 5% das semanas o volume de vendas é alto (estado 4).
6.3 Fluxo de cliques
Varejistas online analisam seus visitantes para aprimorar suas lojas, estratégias de marketing e ofertas de produtos. Muitas informações, como avaliações de clientes, históricos de compras, características demográficas ou a sequência de cliques (clickstream) em uma loja online, estão disponíveis para esse tipo de análise.
Sequências de cliques, os chamados fluxos de cliques, podem ser analisados pela mineração de padrões sequenciais com diversos algoritmos ou com modelos probabilísticos, como Cadeias de Markov (Montgomery et al. 2004). Em contraste com outros conjuntos de dados sequenciais, sequências de cliques são coleções de sequências de dados com tamanhos diferentes. Dois cliques subsequentes podem, além disso, representar o mesmo estado.
Considere, por exemplo, o seguinte extrato de sessões de usuário de uma loja online:
Sessão 1: P1 P2 P1 P3 P4 Defer
Sessão 2: P3 P4 P1 P3 Defer
Sessão 3: P5 P1 P6 P7 P6 P7 P8 P7 Buy
Sessão 4: P9 P2 P11 P12 P11 P13 P11 Buy
Sessão 5: P4 P6 P11 P6 P1 P3 Defer
Sessão 6: P3 P13 P12 P4 P12 P1 P4 P1 P3 Defer
Sessão 7: P10 P5 P10 P8 P8 P5 P1 P7 Buy
Sessão 8: P9 P2 P1 P9 P3 P1 Defer
Sessão 9: P5 P8 P5 P7 P4 P1 P6 P4 Defer
Cada clique em uma das 13 páginas de produtos possíveis na loja online é representado pela letra P e pelo identificador da página do produto (1 a 13), enquanto o clique final representa a decisão do usuário de adiar a compra (Defer) ou comprar um produto (Buy).
Defer e Buy são estados absorventes. O método clickstream é adequado para lidar com clickstreams com e sem estados absorventes. Analisar coleções de clickstreams com R é desafiador, pois
R não suporta diretamente a importação de conjuntos de dados com comprimento de linha variável,
pacotes como o markovchain (Spedicato et al. 2024) suportam apenas a análise de uma única sequência de dados, não coleções de sequências, e
ainda não há um pacote disponível para R que forneça funções para leitura, escrita, agrupamento e análise de clickstreams.
O pacote clickstream fornece uma infraestrutura básica flexível para importar, exportar e analisar conjuntos de clickstreams registrados pela maioria das lojas online. Mais precisamente, o pacote fornece funcionalidades para agrupar clickstreams, visualizá-los e prever cliques futuros para uma determinada sessão.
6.3.1 Clickstreams
Um fluxo de cliques, como vimos, é uma sequência de eventos de clique para exatamente uma sessão com um usuário de loja online. Os fluxos de cliques de diferentes sessões geralmente diferem em tipo e número de eventos de clique. Cada evento de clique é do tipo caractere. Os fluxos de cliques para uma sessão específica podem ser modelados como um vetor, enquanto uma coleção de fluxos de cliques pode ser modelada como uma lista em R.
Sabemos que uma Cadeia de Markov é um processo estocástico \(C_n\) que assume o estado \(m\) de um conjunto finito \(S\) em cada instante de tempo \(n\). Se o estado em \(n\) depende apenas dos \(k\) estados recentes, chamamos \(C_n\) de uma Cadeia de Markov de ordem \(k\). A probabilidade de estar em qualquer um dos \(d\) estados na próxima etapa é, portanto, independente do estado atual em uma Cadeia de Markov de ordem zero.
Cadeias de Markov homogêneas no tempo, onde a probabilidade de transição é independente do tempo \(n\), podem ser descritas por matrizes de transição, \(\mathcal{P}\) onde cada elemento \(p_{ij}\) descreve a probabilidade de obter uma transição do estado \(i\) no tempo \(n-k\) para o estado \(j\) no tempo \(n\). Cada probabilidade \(p_{ij}\) corresponde a um parâmetro a ser estimado.
Cadeias de Markov de ordem superior são, portanto, caracterizadas por \((d-1)d^k\) parâmetros do modelo. O principal desafio de usar Cadeias de Markov de ordem superior para analisar dados de fluxo de cliques é o número de parâmetros que aumenta exponencialmente com a ordem da cadeia. Usuários de uma loja online normalmente decidem visitar um dos vários sites possíveis não apenas com base no site que estão visitando no momento.
Normalmente, um usuário que considera uma página de produto pode adicioná-lo ao carrinho de compras, ver avaliações de produtos, seguir uma recomendação de produto ou pesquisar outro produto. Moe (2003) propõe que a probabilidade de uma transição para qualquer um dos próximos estados possíveis depende do modo (navegação, pesquisa ou compra) em que o usuário está atualmente.
Esse modo pode ser identificado ao considerar os \(k\) estados recentes (sites) de um usuário, em vez de apenas o último estado (site). Cadeias de Markov de ordem superior são, portanto, mais promissoras ao analisar dados de fluxo de cliques.
Exemplo: Clickstream Data for Online Shopping
O Kaggle é uma plataforma comunitária online para cientistas de dados a qual permite que os usuários colaborem com outros usuários, encontrem e publiquem conjuntos de dados. O objetivo desta plataforma online (fundada em 2010 por Anthony Goldbloom e Jeremy Howard e adquirida pelo Google em 2017) é ajudar profissionais e estudantes a alcançar seus objetivos em sua jornada na ciência de dados com as poderosas ferramentas e recursos que ela oferece. Atualmente (2021), existem mais de 8 milhões de usuários registrados no Kaggle.
Um dos conjuntos de dados foi publicado por Lapczynski and Białowąs (2013). Este conjunto de dados contém informações sobre o fluxo de cliques de uma loja online que oferece roupas para gestantes. Os dados são de cinco meses de 2008 entre abril (month = 4) e agosto (month = 8) e incluem, entre outros, a categoria do produto, a localização da foto na página, o país de origem do endereço IP e o preço do produto em dólares americanos.
Informações dos Atributos:
O conjunto de dados contém 14 variáveis descritas em um arquivo separado (consulte ‘Descrição do Conjunto de Dados’). Bojan Tunguz
kaggle: Clickstream Data for Online Shopping, clickstream data for online shopping Data Set
dados = read.csv("https://estatistica.c3sl.ufpr.br/~lucambio/CM/e-shop%20clothing%202008.csv")
names(dados)## [1] "year" "month"
## [3] "day" "order"
## [5] "country" "Session.ID"
## [7] "page.1..main.category." "page.2..clothing.model."
## [9] "colour" "location"
## [11] "model.photography" "price"
## [13] "price.2" "page"head(dados)## year month day order country Session.ID page.1..main.category.
## 1 2008 4 1 1 29 1 1
## 2 2008 4 1 2 29 1 1
## 3 2008 4 1 3 29 1 2
## 4 2008 4 1 4 29 1 2
## 5 2008 4 1 5 29 1 2
## 6 2008 4 1 6 29 1 3
## page.2..clothing.model. colour location model.photography price price.2 page
## 1 A13 1 5 1 28 2 1
## 2 A16 1 6 1 33 2 1
## 3 B4 10 2 1 52 1 1
## 4 B17 6 6 2 38 2 1
## 5 B8 4 3 2 52 1 1
## 6 C56 6 1 2 57 1 4A variável DAY indica o número do dia em cada mês, ORDER indica a sequência de cliques durante cada seção, COUNTRY é uma variável indicadora do país de origem do endereço IP com as seguintes categorias: 1-Australia, 2-Austria, 3-Belgium, 4-British Virgin Islands, 5-Cayman Islands, 6-Christmas Island, 7-Croatia, 8-Cyprus, 9-Czech Republic, 10-Denmark, 11-Estonia, 12-unidentified, 13-Faroe Islands, 14-Finland, 15-France, 16-Germany, 17-Greece, 18-Hungary, 19-Iceland, 20-India, 21-Ireland, 22-Italy, 23-Latvia, 24-Lithuania, 25-Luxembourg, 26-Mexico, 27-Netherlands, 28-Norway, 29-Poland, 30-Portugal, 31-Romania, 32-Russia, 33-San Marino, 34-Slovakia, 35-Slovenia, 36-Spain, 37-Sweden, 38-Switzerland, 39-Ukraine, 40-United Arab Emirates, 41-United Kingdom, 42-USA, 43-biz (.biz), 44-com (.com), 45-int (.int), 46-net (.net), 47-org (*.org).
======================================================== 6. SESSION ID -> variable indicating session id (short record)
======================================================== 7. PAGE 1 (MAIN CATEGORY) -> concerns the main product category:
1-trousers 2-skirts 3-blouses 4-sale
======================================================== 8. PAGE 2 (CLOTHING MODEL) -> contains information about the code for each product
(217 products)
Dispomos de informações gerais acerca do produo. Na variável COLOUR indica-se o color do produto de valores: 1-beige, 2-black, 3-blue, 4-brown, 5-burgundy, 6-gray, 7-green, 8-navy blue, 9-of many colors, 10-olive, 11-pink, 12-red, 13-violet, 14-white. Em LOCATION, informa-se a localização da foto na página, a tela foi dividida em seis partes: 1-top left, 2-top in the middle, 3-top right, 4-bottom left, 5-bottom in the middle, 6-bottom right. Ainda, em MODEL PHOTOGRAPHY, temos duas categorias: 1-en face e 2-profile.
Por último, temos PRICE, o preço do produto em dólares americanos; PRICE 2, uma variável que informa se o preço de um determinado produto é maior que o preço médio de toda a categoria de produtos de valores yes ou no e PAGE, o número da página no site da loja virtual (de 1 a 5).
Nestes dados temos
N = max(dados$Session.ID)
N## [1] 24026seções diferentes, ou seja, temos esse total de usuários diferentes que acessaram a loja online.
Vamos considerar unicamente usuários polacos.
cls <- list()
for(i in 1:N) cls[[i]] <- dados[dados$Session.ID == i & dados$model.photography == 1,8]library(clickstream)class(cls) <- "Clickstreams"summary(cls)## Observations: 24026
##
## Click Frequencies:
## A1 A10 A11 A12 A13 A14 A16 A17 A18 A2 A21 A22 A24 A25 A26 A27
## 2265 2280 2789 2010 1577 926 1388 1531 938 3013 1388 622 288 361 478 420
## A29 A3 A30 A32 A33 A34 A36 A37 A38 A39 A41 A42 A43 A5 A6 A7
## 715 1932 695 643 839 558 356 756 464 473 750 666 246 2354 1802 1622
## A8 A9 B1 B10 B11 B12 B14 B15 B16 B19 B20 B21 B22 B23 B24 B25
## 1585 1923 1552 2566 1543 1369 1118 1499 1408 1110 399 843 402 1391 1879 714
## B26 B27 B28 B29 B3 B30 B31 B32 B34 B4 B6 B7 B9 C11 C12 C13
## 1271 976 504 540 1288 949 1406 1361 631 3579 672 483 840 884 866 935
## C14 C15 C16 C17 C18 C19 C2 C21 C23 C24 C25 C26 C27 C31 C32 C33
## 1178 516 595 1507 651 530 1111 771 310 463 441 521 352 335 359 673
## C34 C35 C36 C37 C38 C39 C4 C41 C42 C44 C45 C46 C48 C5 C51 C52
## 714 764 552 258 343 184 856 758 494 477 302 473 476 1834 605 187
## C53 C54 C55 C6 C7 C9 P1 P12 P13 P16 P18 P19 P2 P20 P21 P22
## 651 680 598 715 1112 928 2681 999 525 1508 570 393 1005 662 297 2
## P24 P26 P3 P30 P33 P35 P36 P37 P38 P39 P4 P40 P45 P46 P47 P49
## 135 590 1149 254 819 219 451 387 218 745 1106 241 150 374 283 671
## P5 P51 P56 P57 P6 P60 P67 P69 P72 P75 P77 P79 P8 P80
## 783 431 335 251 1536 492 218 167 388 169 445 2 661 222mc<- fitMarkovChain(clickstreamList=cls,order=1,control=list(optimizer="quadratic"))mc@transientStates## [1] "A1" "A10" "A11" "A12" "A13" "A14" "A16" "A17" "A18" "A2" "A21" "A22"
## [13] "A24" "A25" "A26" "A27" "A29" "A3" "A30" "A32" "A33" "A34" "A36" "A37"
## [25] "A38" "A39" "A41" "A42" "A43" "A5" "A6" "A7" "A8" "A9" "B1" "B10"
## [37] "B11" "B12" "B14" "B15" "B16" "B19" "B20" "B21" "B22" "B23" "B24" "B25"
## [49] "B26" "B27" "B28" "B29" "B3" "B30" "B31" "B32" "B34" "B4" "B6" "B7"
## [61] "B9" "C11" "C12" "C13" "C14" "C15" "C16" "C17" "C18" "C19" "C2" "C21"
## [73] "C23" "C24" "C25" "C26" "C27" "C31" "C32" "C33" "C34" "C35" "C36" "C37"
## [85] "C38" "C39" "C4" "C41" "C42" "C44" "C45" "C46" "C48" "C5" "C51" "C52"
## [97] "C53" "C54" "C55" "C6" "C7" "C9" "P1" "P12" "P13" "P16" "P18" "P19"
## [109] "P2" "P20" "P21" "P24" "P26" "P3" "P30" "P33" "P35" "P36" "P37" "P38"
## [121] "P39" "P4" "P40" "P45" "P46" "P47" "P49" "P5" "P51" "P56" "P57" "P6"
## [133] "P60" "P67" "P69" "P72" "P75" "P77" "P79" "P8" "P80"mc@absorbingStates## character(0)Nicolau (2014)
6.4 Exercícios
Considere a seguinte sequência categórica de três estados: \[ 1,1,2,3,2,1,1,1,2,2,3,3,2,1,2,3,2,1,2,2,1,2,3,1 \] Estime um modelo de Cadeia de Markov de terceira ordem. Encontre a distribuição estacionária.
Suponha que a sequência de demanda de um produto seja dada pela seguinte expressão: \[ 1,2,5,2,3,4,4,3,4,5,1,2,2,2,4,4,4,5,4,5,3,2,2,2,3,4,4,5,5,5,4,4,3,2 \] Construir e estimar um modelo de Cadeia de Markov de segunda ordem.
Atualizar os dados do Índice Dow Jones Industrial Average (DJIA) e verificar se ainda devemos considerar esse índice numa carteria de investimentos.
Demanda de vendas. Uma aplicação do modelo de Cadeia de Markov multivariada de ordem superior. Uma empresa de refrigerantes em Hong Kong enfrenta um problema interno de planejamento de produção e controle de estoque. Uma questão urgente é o espaço de armazenamento de seu armazém central, que frequentemente se encontra em estado de transbordamento ou próximo da capacidade máxima. A empresa, portanto, precisa urgentemente estudar a interação entre a necessidade de espaço de armazenamento e o crescimento geral da demanda de vendas. O produto pode ser classificado em seis estados possíveis 1, 2, 3, 4, 5 e 6 de acordo com seus volumes de vendas. Todos os produtos são rotulados como 1 = sem volume de vendas, 2 = giro muito lento (volume de vendas muito baixo), 3 = giro lento, 4 = padrão, 5 = giro rápido ou 6 = giro muito rápido (volume de vendas muito alto). Esses rótulos são úteis tanto do ponto de vista do marketing quanto do planejamento da produção. A empresa também gostaria de prever a demanda de vendas de um cliente importante, a fim de minimizar o acúmulo de estoque. Mais importante ainda, a empresa pode compreender o padrão de vendas desse cliente e, então, desenvolver uma estratégia de marketing para lidar com ele. Abaixo é mostrada a demanda de vendas desse cliente para cinco produtos importantes da empresa durante um ano. Esperamos que as sequências de demanda de vendas geradas pelo mesmo cliente sejam correlacionadas entre si. Portanto, ao explorar essas relações, pode-se obter um modelo de Markov multivariado de ordem superior mais adequado para tais sequências de demanda e, consequentemente, obter melhores regras de previsão.
Produto.A=c(6,6,6,6,2,6,2,6,2,2,6,2,6,6,2,6,2,4,4,4,5,6,6,1,2,2,6,
6,6,2,6,2,6,6,2,6,2,2,6,2,1,2,2,6,6,6,2,1,2,6,2,6,6,2,
2,6,2,2,2,6,2,6,2,2,2,2,2,6,2,2,6,6,6,6,1,2,2,6,2,2,2,
2,6,2,2,2,2,3,3,2,3,2,6,6,6,6,2,6,2,6,6,2,6,2,6,6,2,6,
6,2,2,3,4,3,3,1,3,1,2,1,6,1,6,6,1,6,6,2,6,2,6,2,2,2,6,
6,1,6,2,6,1,2,1,6,2,6,2,2,2,2,6,6,1,6,6,2,2,6,2,2,2,3,
4,4,4,6,4,6,1,6,6,1,6,6,6,6,1,6,2,2,2,6,6,6,6,2,6,6,2,
2,6,2,6,2,2,2,6,2,2,2,6,6,6,6,3,2,2,6,2,2,2,2,2,2,6,2,
6,2,2,2,6,2,2,6,6,2,6,6,6,2,2,2,3,3,3,4,1,6,6,1,6,6,1,
6,1,6,6,6,6,1,6,6,6,2,1,2,2,2,2,2,2,3,6,6,6,6,6,2,6)
Produto.B=c(1,6,6,1,6,1,1,1,1,1,1,6,6,6,1,2,1,6,6,1,1,1,6,6,2,1,6,
6,1,1,1,6,1,2,1,6,2,2,2,2,2,6,1,6,6,1,2,1,6,6,6,1,1,1,
6,6,1,1,1,1,6,1,1,2,1,6,1,6,1,1,6,2,6,2,6,6,6,3,6,6,1,
6,6,2,2,2,3,2,2,6,6,6,1,1,6,2,6,6,2,6,2,6,6,1,3,6,6,1,
1,1,2,2,3,2,2,6,2,2,2,1,6,1,6,1,1,6,2,1,1,1,2,2,1,6,1,
1,1,1,2,6,1,1,1,1,6,1,6,1,2,1,6,1,6,6,1,6,1,2,2,2,2,3,
3,2,2,2,6,6,6,6,2,1,1,6,1,1,1,6,1,6,1,6,1,6,1,1,6,6,2,
1,1,6,6,1,1,2,6,2,6,6,6,1,2,6,1,6,1,1,1,1,6,1,6,1,1,6,
6,1,6,6,1,6,1,6,6,1,1,6,6,2,2,2,2,2,2,2,2,2,6,6,6,6,1,
6,6,6,1,6,6,1,6,6,1,1,6,1,3,3,3,5,1,6,6,6,6,6,6,6,6)
Produto.C=c(6,6,6,6,6,6,6,2,6,6,6,6,6,6,6,2,6,6,6,6,2,6,6,6,2,2,6,
6,6,6,6,6,6,1,6,2,6,6,6,6,6,6,6,6,2,6,6,1,2,6,1,6,6,1,
6,2,6,6,6,6,6,6,6,2,6,6,6,2,6,6,1,6,6,6,6,6,6,6,3,3,6,
3,2,1,2,2,1,6,6,1,6,1,6,6,6,6,6,6,1,6,6,6,1,6,6,6,6,6,
6,6,6,6,6,6,2,6,6,6,6,6,6,6,6,2,2,6,6,2,6,1,2,6,6,6,2,
6,6,2,6,6,2,6,1,6,2,6,2,1,2,6,6,2,2,6,2,6,2,2,6,2,6,6,
6,2,2,2,6,6,2,6,6,2,2,6,1,2,1,2,6,6,2,2,6,6,1,2,2,1,6,
2,6,2,2,1,1,5,6,3,6,1,6,6,1,2,2,6,1,6,2,6,6,1,6,2,6,2,
6,6,6,1,6,1,6,6,2,2,2,1,2,3,6,1,6,1,6,1,6,1,6,6,6,1,1,
6,6,6,6,6,1,6,6,6,1,6,1,1,6,6,6,6,6,6,6,6,1,6,6,1,6)
Produto.D=c(6,2,2,2,2,3,3,4,4,4,5,4,3,3,6,2,6,6,6,3,4,4,3,3,3,3,3,
2,6,6,3,4,4,4,4,3,4,2,6,2,2,6,2,2,6,6,3,4,5,4,4,6,3,6,
6,6,2,6,2,6,6,2,2,6,4,4,5,4,3,4,3,4,4,6,2,6,6,2,2,6,2,
6,6,2,6,6,2,6,6,2,6,2,6,3,5,5,5,4,4,4,3,6,2,6,6,2,6,2,
6,2,2,6,2,6,6,2,6,4,4,4,4,4,4,6,3,6,6,2,6,2,6,2,6,2,6,
6,2,2,2,2,2,2,2,2,2,3,3,3,5,5,4,5,3,3,3,6,2,6,6,2,2,6,
2,2,2,2,6,2,3,2,2,3,6,3,2,2,3,4,4,4,4,5,5,4,4,6,6,2,6,
2,6,2,2,2,2,2,2,2,5,5,4,4,5,5,2,6,2,6,6,2,6,2,6,2,2,3,
3,4,4,5,4,4,4,3,4,3,6,2,6,2,2,2,2,2,2,2,2,2,2,2,3,4,4,
4,4,5,4,4,4,3,2,2,2,6,2,2,2,6,2,6,2,6,2,2,2,2,2,3,2)
Produto.E=c(6,2,2,2,2,3,3,4,4,4,5,4,3,3,6,2,6,6,2,3,4,4,3,4,4,3,3,
2,2,6,3,4,4,4,4,3,4,2,3,2,2,6,3,3,6,6,3,4,5,4,5,3,3,2,
6,6,2,6,2,6,6,2,2,6,4,4,4,4,4,4,5,4,4,6,2,6,6,2,2,6,2,
6,6,2,6,6,2,6,6,2,6,2,6,3,4,4,4,4,4,4,4,6,2,6,6,2,6,2,
6,6,6,6,2,6,2,2,6,4,4,4,4,4,4,6,3,3,6,2,2,2,6,2,6,2,2,
2,2,2,2,2,2,2,2,2,2,3,6,4,5,5,5,5,2,4,6,6,2,6,6,2,2,6,
2,2,2,2,6,2,3,2,2,3,6,3,2,2,3,4,4,4,4,5,5,4,3,3,6,2,6,
2,2,2,6,3,2,2,2,2,5,5,4,4,4,4,3,6,2,6,6,2,6,2,6,2,2,3,
3,4,4,5,4,4,4,4,4,3,6,2,6,2,2,2,6,2,2,2,2,2,2,2,3,4,4,
4,4,5,4,4,4,3,2,2,2,6,6,6,2,6,2,6,2,6,2,2,2,2,2,2,2)