Em Estimação dos parâmetros mencionamos a tarefa de selecionar um HMM apropriado para uma série de observações e também a tarefa de verificar o ajuste do modelo selecionado. Agora apresentaremos uma introdução à seleção de modelos nos HMMs e ao respectivo diagnóstico de modelos usando pseudo-resíduos.



1 Seleção de modelos


Um problema que surge naturalmente ao trabalhar com HMMs é o de selecionar um modelo apropriado, isto é, de escolher o número apropriado de estados \(m\). Como os resultados assintóticos acerca da ordem \(m\) de um HMM ainda não são claros veja, por exemplo, Rydén (1995), é necessário especificar um critério para a comparação de modelos.

Suponha que as observações \(s_1,\cdots, s_T\) foram geradas pelo modelo operacional desconhecido \(f\) e que este se encaixa em modelos de duas famílias aproximadas diferentes, \(\{g_1\in G_1\}\) e \(\{g_2\in G_2\}\). A intenção da seleção de modelos é encontrar o modelo que seja melhor em um certo sentido.

Existem duas abordagens diferentes para a seleção de modelos. Na abordagem clássica, seleciona-se a família estimada como a mais próxima do modelo verdadeiro. Para esse propósito, define-se uma discrepância entre os modelos operacional e aproximado, \(\Delta(f,g_1)\) e \(\Delta(f,g_2)\) e, como o modelo operacional seja desconhecido, derivamos a discrepância esperada como critério de seleção, \(\widehat{\mbox{E}}\big( \Delta(f,g_1)\big)\) e \(\widehat{\mbox{E}}\big(\Delta(f,g_2)\big)\).

À escolha da discrepância de Kullback-Leibler leva ao chamado critério de informação de Akaike (AIC), que pode ser calculado da seguinte forma: \[\begin{equation} AIC \, = \, -2\log(L)+2p, \end{equation}\] onde \(\log(L)\) é a log-verossimilhança do modelo ajustado e \(p\) representa o número de parâmetros do modelo. O primeiro termo é uma medida de ajuste, diminui com o aumento do número de estados \(m\). O segundo termo é um termo de penalidade, aumenta com o aumento de \(m\). O AIC é a escolha canônica para comparação de modelos. Para uma introdução mais detalhada ao AIC e sua derivação, ver Zucchini (2000).

A idéia da segunda abordagem da seleção de modelos, a abordagem bayesiana, é selecionar a família estimada com maior probabilidade de ser verdadeira. Em uma primeira etapa, antes de considerar as observações, especificamos os chamados prévios, isto é, as probabilidades que \(f\) derive da família aproximada, \(P(f\in G_1)\) e \(P(f\in G_2)\). Em seguida, em um segundo passo, calcula-se e compara-se as posteriores, ou seja, as probabilidades de que \(f\) pertença à família que se aproxima, dadas as observações \(P(f\in G_1 \, | \, s_1,\cdots,s_T)\) e \(P(f\in G_2 \, | \, s_1,\cdots,s_T)\). Essa abordagem resulta no Critério de Informação Bayesiano (BIC), que tem um termo de penalidade ligeiramente modificado em comparação com o AIC: \[\begin{equation} BIC \, = \, -2\log(L)+p\log(T), \end{equation}\] onde \(\log(L)\) e \(p\) são como para o AIC e \(T\) é o número de observações. Usamos \(T\) aqui para o número de observações, pois estamos lidando com HMMs; mais geralmente, denotamos o número de observações por \(n\). Para mais detalhes sobre o BIC, nos referimos a Zucchini (2000) novamente.

Em comparação com o AIC, a penalidade do BIC tem mais peso, pois em geral detêm esse \(\log(T) > 2\). Assim, o Critério de Informação Bayesiana geralmente favorece modelos com menos parâmetros do que o Critério de Informação Akaike.

Isso também se aplica à série de vendas de sabão, como pode ser visto na Figura 1, que representa o AIC e o BIC em relação ao número de estados \(m\). Os valores exatos de ambos os critérios são fornecidos na Tabela 1. Compare esses valores com as log-verossimilhanças obtidas para os modelos de mistura independentes em Fundamentos. Embora os HMMs exijam mais parâmetros, os valores resultantes da AIC e da BIC são inferiores aos obtidos para os modelos de mistura independentes.

dados.AIC = c(1426,1245,1239,1240,1255,1268)
dados.BIC = c(1429,1259,1270,1296,1342,1393)
n.estados = seq(1,6)
par(mar=c(4,3,1,1), pch=19)
plot(dados.AIC~n.estados, main="", xlab="Número de estados", type="b", ylab="", 
      xlim=c(0.5,6.5), ylim=c(1200,1450), axes=FALSE)
lines(dados.BIC~n.estados, type="b", col="darkviolet")
axis(1)
lablist.y<-as.vector(seq(1200,1450, by=50))
axis(2, at=seq(1200, 1450, by=50), labels = FALSE)
text(y = seq(1200, 1450, by=50), par("usr")[1], labels = lablist.y, srt = 45, pos = 2, xpd = TRUE)
grid()
box()
text(6,1300, "AIC", lwd=3)
text(6,1420, "BIC", lwd=3)

Figura 1: Critérios de seleção de modelos AIC e BIC para a série de vendas de sabão.

\[\begin{equation} \begin{array}{c|cccc} \mbox{Modelos} & k & \log(L) & \mbox{AIC} & \mbox{BIC} \\\hline \mbox{Poisson} & 1 & -712 & 1426 & 1429\\ \mbox{2-estados} & 4 & -619 & 1245 & 1259 \\ \mbox{3-estados} & 9 & -611 & 1239 & 1270 \\ \mbox{4-estados} & 16 & -604 & 1240 & 1296 \\ \mbox{5-estados} & 25 & -603 & 1255 & 1342 \\ \mbox{6-estados} & 36 & -598 & 1268 & 1393 \\\hline \end{array} \end{equation}\]

Tabela 1: Critérios de seleção de modelos AIC e BIC para a série de vendas de sabão.

De acordo com a AIC, o modelo com \(m = 3\) é o mais apropriado, enquanto a BIC seleciona o modelo de dois estados. Assim, o modelo selecionado depende da abordagem de seleção de modelos que se deseja seguir.

Suponha que se considere a AIC como critério de seleção, ou seja, escolhemos o modelo de três estados. As respectivas estimativas dos parâmetros são fornecidas por: \[\begin{equation} \widehat{\Gamma} \, = \, \begin{pmatrix} 0.86 & 0.12 & 0.02 \\ 0.44 & 0.54 & 0.02 \\ 0.00 & 0.30 & 0.70 \end{pmatrix}, \qquad \begin{array}{rcl} \widehat{\delta} & = & (0.72 \quad 0.22 \quad 0.06) \\ \widehat{\lambda} & = & (3.74 \quad 8.44 \quad 14.93) \end{array}\cdot \end{equation}\] As distribuições de estado dependentes componentes desse modelo, juntamente com as distribuições marginais, são ilustradas na Figura 2.

Figura 2: Distribuições de estado dependentes componentes do modelo HMM Poisson com \(m=3\) estados nos três primeiros gráficos, de esquerda para direita. último gráfico mostrando as distribuições marginais. Dados da série de vendas de sabão.

Também é interessante considerar a função de autocorrelação do modelo HMM Poisson com três estados selecionado. Na Figura 3, a função de autocorrelação do modelo ajustado é justaposta à função de autocorrelação amostral. É óbvio que ambas as funções de autocorrelação correspondem em grande parte. Além disso, os momentos marginais do HMM, que também são apresentados na Figura 3, estão próximos dos momentos empíricos das observações. No entanto, é possível aplicar métodos mais sofisticados para o diagnóstico do modelo, conforme serão mostrados na seção a seguir.

Figura 3: Funções de autocorrelação amostral e do modelo HMM de três estados.


2 Validação de modelos


Na seção anterior, consideramos dois critérios para a seleção de modelos nos HMMs. Esses critérios fornecem uma regra de decisão para selecionar o melhor de vários modelos ajustados, no entanto, eles não garantem que o modelo selecionado seja realmente apropriado. Portanto, ainda é preciso avaliar o ajuste do modelo escolhido.

Em geral, os resíduos são uma ferramenta popular para avaliar o ajuste de um modelo. No caso ideal, eles devem distribuídos de forma independente e idêntica. Além disso, é de grande vantagem se, no caso de um modelo válido, eles forem distribuídos \(U(0;1)\) ou \(N(0;1)\), independentemente do modelo ajustado.

É bem sabido que numa regressão padrão, por exemplo, da forma \(y_i = \theta_0 + \theta_1 x_i + \epsilon_i\), os resíduos \(\epsilon_i\) são distribuídos de forma independente e idêntica; no entanto, os estimadores dos resíduos \(\widehat{\epsilon}_i = y_i -\widehat{\theta}_0-\widehat{\theta}_1 x_i\) são distribuídos apenas aproximadamente de forma independente e idêntica. No entanto, pode-se analisar os resíduos por meio de ferramentas gráficas, a fim de avaliar o ajuste do modelo.

Por exemplo, para procurar outliers ou verificar estruturas e dependência específicas, pode-se desenhar um gráfico do índice dos resíduos ou mostrá-los na variável dependente ou covariáveis. Caso certos padrões ocorram, o modelo deve ser revisado e melhorado, por exemplo, adicionando termos quadráticos ou cúbicos ou mesmo novas covariáveis. Outra possibilidade é testar as premissas distribucionais usando histogramas ou gráficos \(qq\) dos resíduos. Algumas dessas parcelas residuais são ilustradas na Figura 4.

Figura 4: Gráficos de resíduos.

Quando se trata de HMMs, a análise de resíduos torna-se um pouco mais complicada devido ao fato de que os resíduos não estão sequer aproximadamente distribuídos de forma independente e idêntica. Normalmente, cada um dos estados do HMM emite observações a partir de uma distribuição diferente. Assim, distribuições diferentes também devem ser ajustadas aos resíduos de acordo com os diferentes estados do HMM. A dificuldade reside no fato de que, em geral, a Cadeia de Markov subjacente é desconhecida e, portanto, os resíduos não podem ser ordenados de forma adequada. Uma solução para este problema é o uso dos chamados pseudo-resíduos, que iremos introduzir a seguir.

O conceito de pseudo-resíduos, que se baseia no conceito de \(p\)-valores, satisfaz as propriedades desejadas dos resíduos delineadas no início desta seção pelo menos aproximadamente, sendo assim considerados um meio adequado para avaliar o ajuste de um modelo. Os pseudo-resíduos podem ser definidos para quase todos os modelos e, caso o modelo seja válido, são distribuídos de forma independente e idêntica de forma uniforme ou normal.

Para a construção de pseudo-resíduos consideramos o seguinte teorema: Que \(X\) seja uma variável aleatória com função de distribuição \(F\). Então, \(U = F(X)\) é uniformemente distribuída no intervalo \((0;1)\), ou seja: \[\begin{equation} U \, = \, F(X) \, \sim \, U(0,1)\cdot \end{equation}\]

A prova é deixada como um exercício para o leitor. Com base neste teorema, uma primeira versão do pseudo-resíduo de uma observação \(x_i\) de uma variável aleatória contínua \(X_i\) pode ser definida como a probabilidade de obter uma observação inferior a \(x_i\) no modelo ajustado: \[\begin{equation} u_i \, = \, P\big(X_i\leq x_i \big) \, = \, F(x_i)\cdot \end{equation}\]

Dada a validade do modelo ajustado este tipo de pseudo-resíduo é distribuído \(U(0;1)\), com resíduos para observações extremas próximas de 0 ou 1 (Zucchini and MacDonald, 1999). A construção do chamado pseudo-resíduo uniforme é ilustrada na Figura 5.

Figura 5: Construção de pseudo-resíduos uniformemente distribuídos.

Assim, utilizando o conceito de pseudo-resíduos uniformes, observações de diferentes distribuições podem ser comparadas. Suponha ter observações \(x_1,\cdots,x_n\) e um modelo \(X_i \sim F_i\), \(i = 1,\cdots,n\); isto é, cada \(x_i\) tem a sua própria função de distribuição e portanto os valores \(x_i\) não podem ser comparados diretamente. No entanto, o conceito de pseudo-resíduos pode ser utilizado para gerar pseudo-resíduos \(u_i\) que, no caso do modelo estar correto, são independentes e identicamente distribuídos \(U(0;1)\) e, portanto, podem ser comparados. Os respectivos gráficos de resíduos podem parecer como os mostrados na Figura 6.

Figura 6: Gráficos de diagnóstico para pseudo-resíduos uniformes.

Se o histograma e o gráfico \(qq\)-plot dos pseudo-resíduos uniformes \(u_i\) não parecerem como deveriam, pode-se deduzir que os resíduos não são distribuídos uniformemente e, portanto, o modelo não é válido. Certamente, o conceito de pseudo-resíduos uniformes é muito útil, no entanto, pode levar a problemas na identificação externa. Por exemplo, considere o gráfico de índice fornecido na Figura 6 e observe os valores próximos a 0 ou 1. É difícil ver se um valor é muito improvável ou não. Como um valor de 0.999 é difícil de distinguir de um valor de 0.97, o gráfico de índice é quase inútil para uma análise visual rápida.

No entanto, essa deficiência dos pseudo-resíduos uniformes pode ser superada facilmente usando outro teorema. Seja \(\Phi\) a função de distribuição Normal padrão e \(X\) uma variável aleatória com a função de distribuição \(F\). Então, \[\begin{equation} Z \, = \, \Phi^{-1}\big( F(X)\big) \, \sim \, N(0,1), \end{equation}\] ou seja, \(Z\) definida acima é Normal padrão.

A prova é deixada novamente como um exercício para o leitor. Usando esse teorema, uma segunda versão de pseudo-resíduos surge de uma simples transformação dos pseudo-resíduos uniformes por meio da função de distribuição Normal padrão: \[\begin{equation} z_i \, = \, \Phi^{-1}(u_i) \, = \, \Phi^{-1}\big( F_{X_i}(x_i)\big)\cdot \end{equation}\]

Caso o modelo ajustado seja válido, esses chamados pseudo-resíduos normais são distribuídos \(N(0;1)\), com o valor do resíduo igual a 0, caso a observação atinja à mediana. Na medida em que os pseudo-resíduos normais medem o desvio da mediana e não da esperança (Zucchini, 2002). A construção de pseudo-resíduos normais é ilustrada na Figura 7.

Figura 7: Construção de pseudo-resíduos normais.

Assim, dado que as observações \(x_1,\cdots,x_n\) foram efectivamente geradas pelo modelo \(X_i\sim F_i\), os pseudo-resíduos uniformes \(z_i\) devem seguir uma distribuição normal padrão que pode ser verificada através da análise visual do histograma ou \(qq\)-plot ou através da realização de testes de normalidade. Alguns exemplos de gráficos de diagnóstico para pseudo-resíduos normais são apresentados na Figura 8.

Figura 8: Gráficos de diagnóstico para pseudo-resíduos normais.

A segunda versão dos pseudo-resíduos tem a vantagem de que o valor absoluto do resíduo aumenta com um desvio crescente da mediana e que os extremos podem ser identificadas mais facilmente numa escala normal (Zucchini and MacDonald, 1999). Isto torna-se óbvio quando se comparam os gráficos de índizes da escala uniforme e normal dos pseudo-resíduos apresentados nas Figuras 6 e 8. Enquanto no caso dos pseudo-resíduos uniformes foi difícil distinguir um valor de 0.999 de um valor de 0.97, no gráfico de índizes dos pseudo-resíduos normais é claramente visível que a observação pertencente ao resíduo \(\Phi^{-1}(0.999)\) é um outlier, enquanto que o residual \(\Phi^{-1}(0.97)\) não difere notavelmente dos outros resíduos.

Note-se que a teoria dos pseudo-resíduos delineada até agora pode ser aplicada a apenas distribuições contínuas. No caso de observações discretas, os pseudo-resíduos devem ser modificados para permitir a discrepância dos dados, para o seguinte compare os tesxtos de Zucchini (2002) ou Stadie (2003). Os pseudo-resíduos na situação de distribuições discretas não são definidos como pontos, mas como intervalos. Assim, para uma variável discreta aleatória \(X_i\) com função de distribuição \(F_{X_i}\) obtém-se os segmentos pseudo-resíduos uniformes \[\begin{equation} \big( u^-; u^+ \big) \, = \, \big( F_{X_i}(x_i^-); F_{X_i}(x_i)\big), \end{equation}\] com \(x_i\) denotando a maior realização possível, que é menor que \(x_i\) e os segmentos dos pseudo-resíduos Normais \[\begin{equation} \big( u^-; u^+ \big) \, = \, \big( F_{X_i}(x_i^-); F_{X_i}(x_i)\big) \cdot \end{equation}\]

A construção do segmento dos pseudo-resíduos normais de uma variável aleatória discreta é ilustrada na Figura 9.

Figura 9: Construção dos pseudo-resíduos normais no caso discreto.

Ambas as versões de segmentos de pseudo-resíduos contêm informações sobre o quão extremas e raras são as observações. Por exemplo, o limite inferior \(u_i^-\) do intervalo pseudo-residual uniforme especifica a probabilidade de observar um valor inferior a \(x_i\), \(1-u_i^+\) fornece o valor da probabilidade de um valor superior a \(x_i\) e a diferença \(u_i^+-u_i^-\) é igual à probabilidade da observação \(x_i\) sob o modelo ajustado.

Em analogia com o caso contínuo, os segmentos pseudo-residuais podem ser interpretados como realizações censuradas ao intervalo de uma distribuição uniforme ou Normal padrão, se o modelo ajustado for válido. Embora isto esteja correcto apenas no caso dos parâmetros do modelo ajustado serem conhecidos, ainda assim está aproximadamente correcto se o número de parâmetros estimados for pequeno comparado com o tamanho da amostra (Stadie, 2003).

Figura 10: Gráficos de diagnóstico para pseudo-resíduos discretos.

Certamente, é fácil construir um gráfico de índizes dos segmentos pseudo-residuais ou plotar contra qualquer variável independente ou dependente. No entanto, para construir um \(qq\)-plot dos segmentos pseudo-residuais é preciso encontrar um método para definir a ordem dos pseudo-resíduos. Uma possibilidade de ordenar os segmentos pseudo-residuais é considerar os chamados pseudo-resíduos medianos ou mid-pseudo-resíduos que são definidos como \[\begin{equation} z_i^m \, = \, \Phi^{-1}\Big( \dfrac{u_i^-+u_i^+}{2}\Big)\cdot \end{equation}\]

O \(qq\)-plot mostrado na Figura 10 foi construído ordenando os pseudo-resíduos de acordo com os respectivos mid-pseudo-resíduos. Além disso, os pseudo-resíduos medianos podem ser utilizados para a verificação da normalidade, por exemplo, por meio de um histograma dos pseudo-resíduos medianos.

Agora, tendo delineado a teoria básica dos pseudo-resíduos, podemos considerar o uso destes no contexto dos HMMs. A análise dos pseudo-resíduos de um HMM serve para dois propósitos, a saber, a avaliação do ajuste geral de um modelo selecionado e a detecção de observações aberrantes. Dependendo dos aspectos do modelo a ser analisado é possível distinguir entre duas abordagens de computação de pseudo-resíduos de um HMM.

A primeira técnica considera as observações uma de cada vez e procura aquelas que, em relação ao modelo e todas as outras observações da série, são suficientemente extremas para sugerir que são de natureza ou origem diferente das outras. Isto significa que se calculam os pseudo-resíduos a partir da distribuição condicional de \(S_u\) dadas todas as outras observações, \(S^{(-u)} = s^{(-u)}\) de acordo com \[\begin{equation} z_i \, = \, \Phi^{-1}\Big( P\big(S_u \leq s_u \, | \, S^{(-u)} = s^{(-u)}\big)\Big) \end{equation}\] para distribuições contínuas de componentes e \[\begin{equation} \big(z_i^-; z_i^+\big) \, = \, \left( \Phi^{-1}\Big( P\big(S_u < s_u \, | \, S^{(-u)} = s^{(-u)}\big)\Big); \Phi^{-1}\Big( P\big(S_u \leq s_u \, | \, S^{(-u)} = s^{(-u)}\big)\Big) \right) \end{equation}\] para distribuições discretas de componentes, onde em ambos os casos a distribuição condicional

Apresentamos um novo exemplo, desta vez completo, ou seja, são dados diferentes daqueles tratados em exemplos anteriores e mostraremos neles os diversos procedimentos estudados: modelagem, estimação, seleção de modelos, resíduos, previsão. Depois apresentaremos os resultados da estimação da parte oculta, a Cadeia de Markov, de duas formas: as distribuições de estado dependentes e os estados estimados para cada observação. Mostraremos também previsões.


Exemplo 1: Terremotos de grandes magnitudes.


Considere, como exemplo, a série de contagens anuais de grandes terremotos, ou seja, magnitude 7 e acima, para os anos 1900-2019, ambos inclusive. Um exame da Figura 11 sugere que pode haver alguns períodos com uma baixa taxa de terremotos e alguns com uma taxa relativamente alta. Os HMMs, que permitem que a distribuição de probabilidade de cada observação dependa do estado não observado ou oculto de uma Cadeia de Markov, podem acomodar tanto a dispersão excessiva quanto a dependência serial.

# Dados
#
no = c(10,10,12,9,13,15,17,12,9,14,11,10,9,13,11,9,17,11,12,5,
       6,6,9,12,6,6,7,7,14,11,5,14,9,5,8,5,7,10,10,7,7,9,7,14,
       9,3,13,7,10,11,15,6,6,10,8,7,10,18,8,10,6,9,6,13,8,11,9,
       6,19,11,13,11,11,10,7,11,15,10,11,6,5,6,4,9,4,10,4,8,4,
       3,10,7,12,6,9,13,8,10,5,8,8,8,11,11,11,9,7,15,6,13,16,
       13,12,12,11,12,13,6,12,9)
#
# Fonte: https://www.ngdc.noaa.gov/nndc/struts/form?t=101650&s=35&d=35
#
anos = seq(as.Date("1900-01-01"), length = 120, by = "year")
dados.HMM = data.frame(anos, no)
library(ggplot2)
ggplot(data = dados.HMM, aes(x = anos, y = no)) + 
    geom_line(color = "#00AFBB", size = 1) + labs(x = "Anos", y = "No. de grandes terremotos")

Figura 11: Número de grandes terremotos, magnitude 7 e maior, no mundo, 1900-2019.

A proposta é utilizarmos a libraria de funções HiddenMarkov, nesta utilizamos a função dthmm para estrurarmos o modelo HMM Poisson com número de estados 2 e 3. Ainda utilizamos a função BaumWelch para estimarmos os parâmetros pelo algorito EM.

x = runif(101)
library(HiddenMarkov)
#### No. de estados = 2
Pi = matrix(rep(1/2,4), byrow = TRUE, nrow = 2)
delta = cbind(x,1-x)
ajuste.2 = dthmm(no, Pi, delta[1,], "pois", list(lambda=c(2,1)))
ajuste.2.fit = BaumWelch(ajuste.2, control = bwcontrol(prt = FALSE, maxiter = 1500))
## EM converge a um máximo local, por isso, é aconselhável que se ajuste a um modelo
## repetidamente com diferentes valores iniciais
for(i in 1:100){
      ajuste.2 = dthmm(no, Pi, delta[i+1,], "pois", list(lambda=c(2,1)))
      try(tfm1 <- BaumWelch(ajuste.2, control = bwcontrol(prt = FALSE, maxiter = 1500)))
      if (logLik(tfm1) > logLik(ajuste.2.fit)) ajuste.2.fit = tfm1
}
#
#### No. de estados = 3
Pi = matrix(rep(1/3,9), byrow = TRUE, nrow = 3)
delta = cbind(x,(1-x)/3,2*(1-x)/3)
ajuste.3 = dthmm(no, Pi, delta[1,], "pois", list(lambda=c(3,2,1)))
ajuste.3.fit = BaumWelch(ajuste.3, control = bwcontrol(prt = FALSE, maxiter = 1500))
#
# Isso a seguir pode ser um pouco demorado
for(i in 1:100){
   ajuste.3 = dthmm(no, Pi, delta[i+1,], "pois", list(lambda=c(3,2,1)))
   try(tfm1 <- BaumWelch(ajuste.3, control = bwcontrol(prt = FALSE, maxiter = 1500)))
   if (logLik(tfm1) > logLik(ajuste.3.fit)) ajuste.3.fit = tfm1
}
#
#### No. de estados = 4
Pi = matrix(rep(1/4,16), byrow = TRUE, nrow = 4)
delta = cbind(x,(1-x)/6,2*(1-x)/6,3*(1-x)/6)
ajuste.4 = dthmm(no, Pi, delta[1,], "pois", list(lambda=c(4,3,2,1)))
ajuste.4.fit = BaumWelch(ajuste.4, control = bwcontrol(prt = FALSE, maxiter = 1500))
#
# Isso a seguir pode ser um pouco demorado
for(i in 1:100){
   ajuste.4 = dthmm(no, Pi, delta[i+1,], "pois", list(lambda=c(4,3,2,1)))
   try(tfm1 <- BaumWelch(ajuste.4, control = bwcontrol(prt = FALSE, maxiter = 1500)))
   if (logLik(tfm1) > logLik(ajuste.4.fit)) ajuste.4.fit = tfm1
}

Observemos agora os resíduos de cada um dos modelos. Como sempre, observar os resíduos tem como objetivo identificar observações atípicas, desvios de comportamento. Quaisquer desses detalhes e outros podem indicar afastamento do comportamento Normal padrão dos resíduos e, portanto, indicar que os modelos propostos não sejam apropriados.

ajuste.2.res = residuals(ajuste.2.fit)
ajuste.3.res = residuals(ajuste.3.fit)
ajuste.4.res = residuals(ajuste.4.fit)
library(car); library(astsa)
par(mar=c(4,4,1,1), pch=19, mfrow=c(2,2))
plot(ajuste.2.res, ylab="Resíduos", xlab="", ylim = c(-3,3))
abline(h=c(-2.58,-1.96,0,1.96,2.58), lty=2, lwd=2, col="blue")
grid()
hist(ajuste.2.res, freq = FALSE, main = NULL, ylab = NULL, xlab = NULL, xlim = c(-3,3), col="gray")
curve(dnorm(x, mean=mean(ajuste.2.res), sd = sd(ajuste.2.res)), add=TRUE, col="red", lwd = 2) #linha
grid()
box()
qqPlot(ajuste.2.res, pch = 19, ylab = "Quantis amostrais", xlab = "Quantis normais", id = FALSE)
plot(acf(ajuste.2.res,plot=F)[1:20],lwd = 2)
grid()

par(mar=c(4,4,1,1), pch=19, mfrow=c(2,2))
plot(ajuste.3.res, ylab="Resíduos", xlab="", ylim = c(-3,3))
abline(h=c(-2.58,-1.96,0,1.96,2.58), lty=2, lwd=2, col="blue")
grid()
hist(ajuste.3.res, freq = FALSE, main = NULL, ylab = NULL, xlab = NULL, xlim = c(-3,3), col="gray")
curve(dnorm(x, mean=mean(ajuste.3.res), sd = sd(ajuste.3.res)), add=TRUE, col="red", lwd = 2) #linha
grid()
box()
qqPlot(ajuste.3.res, pch = 19, ylab = "Quantis amostrais", xlab = "Quantis normais", id = FALSE)
plot(acf(ajuste.3.res,plot=F)[1:20],lwd = 2)
grid()

par(mar=c(4,4,1,1), pch=19, mfrow=c(2,2))
plot(ajuste.4.res, ylab="Resíduos", xlab="", ylim = c(-3,3))
abline(h=c(-2.58,-1.96,0,1.96,2.58), lty=2, lwd=2, col="blue")
grid()
hist(ajuste.4.res, freq = FALSE, main = NULL, ylab = NULL, xlab = NULL, xlim = c(-3,3), col="gray")
curve(dnorm(x, mean=mean(ajuste.4.res), sd = sd(ajuste.4.res)), add=TRUE, col="red", lwd = 2) #linha
grid()
box()
qqPlot(ajuste.4.res, pch = 19, ylab = "Quantis amostrais", xlab = "Quantis normais", id = FALSE)
plot(acf(ajuste.4.res,plot=F)[1:20],lwd = 2)
grid()

Figura 12: Terremotos de grandes magnitudes: pseudo-resíduos.

Percebemos claramente que os resíduos mostram-se similares conforme o número de estados. Significa que desta maneira não podemos decidir qual destes modelos escolher. Mostramos a seguir o comportamento dos critérios AIC e BIC com este objetivo. 

AIC.HMM = c(-2*ajuste.2.fit$LL+2*(2+1*2),-2*ajuste.3.fit$LL+2*(3+2*3),-2*ajuste.4.fit$LL+2*(4+3*4))
BIC.HMM = c(-2*ajuste.2.fit$LL+(2+1*2)*log(length(ajuste.2.fit$x)),
            -2*ajuste.3.fit$LL+(3+2*3)*log(length(ajuste.3.fit$x)),
            -2*ajuste.4.fit$LL+(4+3*4)*log(length(ajuste.4.fit$x)))
n.estados = seq(2,4)
par(mar=c(4,3,1,1), pch=19)
plot(AIC.HMM~n.estados, main="", xlab="Número de estados", type="b", ylab="", 
         xlim=c(1.5,4.5), ylim=c(600,700), axes=FALSE)
lines(BIC.HMM~n.estados, type="b", col="darkviolet")
axis(1)
lablist.y<-as.vector(seq(600,700, by=10))
axis(2, at=seq(600, 700, by=10), labels = FALSE)
text(y = seq(600, 700, by=10), par("usr")[1], labels = lablist.y, srt = 45, pos = 2, xpd = TRUE)
grid()
box()
text(4.2,645, "AIC", lwd=3)
text(4.2,690, "BIC", lwd=3)

Figura 13: Critérios de seleção AIC e BIC para a série de terremotos de grandes magnitudes.

Observando a Figura 13 percebemos claramente que o modelo mais adequado a esta situação é aquele com 2 estados. Vejamos agora a estimação dos estados para cada observação.

# Mostrando os resultados do modelos com dois estados
##-- Obtendo os Parâmetros --##
u = ajuste.2.fit$pm$lambda # garantindo que o estado 1 tenha lambda menor
if (u[1] <= u[2]) { lambda.mle = c(u[1], u[2]) } else { lambda.mle = c(u[2], u[1]) }
mtrans = ajuste.2.fit$Pi
lams = lambda.mle
pi1 = mtrans[2,1]/(2 - mtrans[1,1] - mtrans[2,2]); pi2 = 1-pi1
##-- Gráficos --##
layout(matrix(c(1,2,1,3), 2))
par(mar = c(3,3,1,1), mgp = c(1.6,.6,0))
# dados e estados
plot(ts(no, start=1900), main="", xlab="", ylab='No.', type='h', col=gray(.7))
text(ts(no, start=1900), col=6*Viterbi(ajuste.2.fit)-2, labels=Viterbi(ajuste.2.fit), cex=.9)
grid()
# probabilidade do estado 2
plot(ts(ajuste.2.fit$u[,2], start=1900), ylab = expression(hat(pi)[~2]*'(t|n)'), xlab="")
abline(h=.5, lty=2)
grid()
# histograma
hist(no, breaks=30, prob=TRUE, main="", ylab="Densidade", xlab="")
box()
xvals = seq(1,45)
u1 = pi1*dpois(xvals, lams[1])
u2 = pi2*dpois(xvals, lams[2])
lines(xvals, u1, col=4); lines(xvals, u2, col=2)
grid()

Figura 14: Superior: Série do número anual de terremotos de grande magnitude e estados estimados. Em baixo, à esquerda: Probabilidades de suavização. Em baixo à direita: histograma dos dados com as duas funções de probabilidade estimadas sobrepostas (linhas sólidas).


3 Exercícios


    1. Seja \(X\) uma variável aleatória de valor contínuo com a função de distribuição \(F\). Mostre que a variável aleatória \(U = F(X)\) é distribuída uniformemente no intervalo \([0,1]\), ou seja, \(U\sim U(0,1)\).

    2. Suponha que \(U\sim U(0,1)\) e cpnsidere \(F\) como a função de distribuição de uma variável aleatória com valor contínuo. Mostre que a variável aleatória \(X = F^{-1}(U)\) tem função de distribuição dada por \(F\).

      1. Forneça a expressão explícita para \(F^{-1}\) para a distribuição exponencial, ou seja, a função de densidade \(f(x)=\lambda e^{-\lambda x}\), \(x\geq 0\).

      2. Como verificação do resultado anterior, considere gerar 1000 variáveis aleatórias distribuídas uniformemente, transformando-as usando \(F^{-1}\) e examinando o histograma dos valores resultantes.

    4. Mostre que para uma variável aleatória \(X\) de valor contínuo com função de distribuição \(F\), a variável aleatória \(Z = \Phi^{-1}\big(F(X)\big)\), onde \(\Phi\) é a função de distribuição Normal padrão, isto é, \(Z\sim N(0,1)\).


  2. Considere um modelo HMM Exponencial com \(m\) estados.

    1. Estime os parâmetros do modelo utilizando a seguinte série de observações para \(m = 1,2,3\): \[\begin{equation} \begin{array}{cccccccccc} 14.8, & 35.0, & 13.2, & 4.4, & 0.1, & 0.3, & 0.8, & 65.6, & 1.0, & 1.1, \\ 2.1, & 10.0, & 0.1, & 8.2, & 17.5, & 69.4, & 15.4, & 0.3, & 0.1, & 0.7 \end{array} \end{equation}\]

    2. Calcule a previsão \(h\) passos à frente para um HMM Exponencial utilizando a série anteior e considere \(h=1,2\).

    3. Escreva uma função que calcule os pseudo-resíduos normais um passo à frente (\(h=1\)) para um HMM exponencial.

    Observe que a função de distribuição da previsão com uma etapa de antecedência no momento \(t\) é dada por \[\begin{equation} P\big( S_{t+1}\leq s \, | \, S^{(t)}=s^{(t)}\big) \, = \, \sum_{i=1}^m \xi_i F_i(s), \end{equation}\] onde \(\xi_i\) é a \(i\)-ésima entrada do vetor \(\displaystyle \dfrac{\alpha_t \Gamma}{\alpha_t\pmb{1}^\top}\). Note também que a função de distribuição no momento \(t = 0\) é \[\begin{equation} P\big(S_1\leq s \big) \, = \, \sum_{i=1}^m \xi_i F_i(s)\cdot \end{equation}\] Use sua função para computar os pseudo-resíduos normais de previsão para a séérie dada em (a). Use gráficos, por exemplo, histograma, \(qq\)-plot, etc., para examinar estes resíduos.

  3. Considere o processo \(AR(1)\) \[\begin{equation} S_t \, = \, \phi S_{t-1}+\epsilon_t, \end{equation}\] com \(\epsilon_t\sim N(0,1)\) independentes. Segue que \(|\phi| < 1\) e \(\mbox{Var}(S-t)=1/(1-\phi^2)\).

  4. Gerar um modelo HMM Poisson estacionário com dois estados e matriz de probabilidades de transição \[\begin{equation} \Gamma \, = \, \begin{pmatrix} 0.6 & 0.4 \\ 0.4 & 0.6 \end{pmatrix} \end{equation}\] e com médias das distribuições estado dependentes \(\lambda_1=2\) e \(\lambda_2=5\) para \(t=1,2,\cdots,100\) e \(\lambda_1=2\) e \(\lambda_2=7\) para \(t=101,\cdots,120\). Ajuste o modelo às primeiras 80 observações e utilize os pseudo-resíduos de previsão para monitorar as próximas 40 observações para verificar evidências de uma mudança.