Nos últimos anos, a Análise do Espectro Singular (SSA), uma técnica poderosa utilizada na análise de séries temporais, tem sido desenvolvida e aplicada a diversos problemas práticos. Neste artigo, o desempenho da técnica SSA foi avaliado por meio de sua aplicação a um conjunto de dados de séries temporais bem conhecido, a saber, o número mensal de mortes acidentais nos EUA. Os resultados demonstram que a técnica SSA fornece uma previsão precisa.

A técnica SSA é um método de decomposição de séries temporais, decompõe uma série temporal em um conjunto de componentes somáveis que são agrupados e interpretados como tendência, periodicidade e ruído. A SSA enfatiza a separabilidade dos componentes subjacentes e pode separar periodicidades que ocorrem em diferentes escalas de tempo, mesmo em dados de séries temporais muito ruidosos. A série temporal original é recuperada somando-se todos os seus componentes.

Assim, a SSA pode ser usada para analisar e reconstruir uma série temporal com ou sem diferentes componentes, conforme desejado. Por exemplo, você pode aplicar a SSA para:

  1. construir uma versão suavizada de uma série temporal usando um pequeno subconjunto de seus componentes;

  2. investigar os componentes periódicos de uma série temporal para entender os processos subjacentes que a geraram;

  3. reconstruir a série temporal original sem seus componentes periódicos;

  4. remover toda a tendência e os componentes periódicos da série, deixando apenas o ‘ruído’, que pode ser significativo por si só…

e assim por diante.

Ao contrário do método ARIMA (Autoregressive Integrated Moving Average), comumente utilizado, o SSA não faz suposições sobre a natureza da série temporal e possui apenas um único parâmetro ajustável.

1 Intodução


A Análise do Espectro Singular (SSA, na sigla em inglês) é uma técnica inovadora e poderosa de análise de séries temporais que incorpora elementos da análise clássica de séries temporais, estatística multivariada, geometria multivariada, sistemas dinâmicos e processamento de sinais.

As possíveis áreas de aplicação da SSA são diversas: da matemática e física à economia e matemática financeira, da meteorologia e oceanografia às ciências sociais e pesquisa de mercado. Qualquer série aparentemente complexa com uma estrutura potencial pode fornecer um exemplo de aplicação bem-sucedida da SSA (N. Golyandina, Nekrutkin, and Zhigljavsky 2001).

O objetivo da SSA é decompor a série original na soma de um pequeno número de componentes independentes e interpretáveis, como uma tendência de variação lenta, componentes oscilatórios e ruído sem estrutura. A SSA é uma ferramenta muito útil que pode ser usada para resolver os seguintes problemas:

  1. encontrar tendências de diferentes resoluções;
  2. suavização;
  3. extração de componentes de sazonalidade;
  4. extração simultânea de ciclos com períodos curtos e longos;
  5. extração de periodicidades com amplitudes variáveis;
  6. extração simultânea de tendências e periodicidades complexas;
  7. identificação de estrutura em séries temporais curtas; e
  8. detecção de pontos de mudança.

A resolução de todos esses problemas corresponde às capacidades básicas da SSA. Para alcançar as capacidades da SSA mencionadas acima, não precisamos conhecer o modelo paramétrico da série temporal considerada.

O surgimento da SSA é geralmente associado à publicação de artigos de Broomhead, por exemplo, Broomhead and King (1986), embora as ideias da SSA tenham sido desenvolvidas independentemente na Rússia (São Petersburgo, Moscou) e em diversos grupos no Reino Unido e nos EUA. Atualmente, existem centenas de artigos que abordam os aspectos metodológicos e as aplicações da SSA ver, por exemplo, Vautard, Yiou, and Ghil (1992); Ghil and Taricco (1997); Allen and Smith (1996); Danilov and Zhigljavsky (1997); Yiou, Sornette, and Gill (2000) e referências neles contidas.

Uma descrição completa dos fundamentos teóricos e práticos da técnica SSA (com diversos exemplos) pode ser encontrada em Danilov and Zhigljavsky (1997) e N. Golyandina, Nekrutkin, and Zhigljavsky (2001). Uma introdução elementar ao tema pode ser encontrada em Elsner and Tsonis (1996).

O fato de a série temporal original dever satisfazer uma fórmula linear recorrente (FLR) é uma propriedade importante da decomposição SSA. Geralmente, o método SSA deve ser aplicado a séries temporais regidas por fórmulas lineares recorrentes para prever novos pontos de dados. Existem dois métodos para construir intervalos de confiança com base na técnica SSA: o método empírico e o método bootstrap.

Os intervalos de confiança empíricos são construídos para a série completa, assumindo-se que ela manterá a mesma estrutura no futuro. Os intervalos de confiança bootstrap são construídos para a continuação do sinal, que são os principais componentes da série completa (N. Golyandina, Nekrutkin, and Zhigljavsky 2001).

Séries temporais reais frequentemente contêm dados faltantes, o que impede a análise e reduz a precisão dos resultados. Existem diferentes métodos baseados em SSA para preencher conjuntos de dados faltantes ver, por exemplo, Schoellhamer (2001); Kondrashov, Feliks, and Ghil (2005); Nina Golyandina and Osipov (2007); Kondrashov and Ghil1 (2006).

A detecção de pontos de mudança em séries temporais é um método que identifica se a estrutura da série sofreu alterações em algum ponto no tempo, causadas por algum fator. O método de detecção de pontos de mudança descrito em Moskvina and Zhigljavsky (2003) baseia-se na aplicação sequencial do SSA a sub-séries da série original e monitora a qualidade da aproximação das demais partes da série por aproximações adequadas.

Vale mencionar também que métodos automáticos para a identificação dos principais componentes de séries temporais dentro da estrutura da SSA foram desenvolvidos recentemente ver, por exemplo, Alexandrov and Golyandina (2004a); Alexandrov and Golyandina (2004b).

Neste artigo, iniciamos com uma breve descrição da metodologia da SSA e finalizamos aplicando essa técnica à série original, especificamente, o número mensal de mortes acidentais nos EUA, série USAccDeaths.

2 Metodologia


Considere a série temporal real não nula \(Y_T = (y_1,\cdots,y_T)\) de comprimento \(T\) suficiente. O principal objetivo da Análise de Espectro Singular (SSA) é decompor a série original em uma soma de séries, de modo que cada componente dessa soma possa ser identificado como uma tendência, um componente periódico ou quase periódico (possivelmente modulado em amplitude) ou ruído. Isso é seguido pela reconstrução da série original.

A técnica SSA consiste em duas etapas complementares: decomposição e reconstrução, ambas com dois passos distintos. Na primeira etapa, decompomos a série e, na segunda, reconstruímos a série original e utilizamos a série reconstruída (sem ruído) para prever novos pontos de dados. A seguir, apresentamos uma breve discussão sobre a metodologia da técnica SSA. Para isso, seguimos principalmente N. Golyandina, Nekrutkin, and Zhigljavsky (2001).

2.1 Etapa 1: Decomposição


Primeiro passo: Incorporação

A incorporação pode ser vista como um mapeamento que transforma uma série temporal unidimensional \(Y_T = (y_1,\cdots,y_T)\) em uma série multidimensional \(X_1,\cdots,X_K\) com vetores \(X_i = (y_i,\cdots,y_{i+L-1})^\top\in\mathbb{R}^L\) onde \(K=T-L+1\). Os vetores \(X_i\) são chamados vetores \(L\)-defasados ou, simplesmente, vetores defasados.

O único parâmetro da incorporação é o comprimento da janela \(L\), um inteiro tal que \(2\leq L\leq T/2\). O resultado deste passo é a matriz de trajetória \[ \pmb{X}=\big[ X_1,\cdots,X_K\big] = (x_{ij})_{i,j=1}^{L\, K}\cdot \]

Note que a matriz de trajetória \(\pmb{X}\) é uma matriz de Hankel, o que significa que todos os elementos ao longo da diagonal \(i+j = \mbox{const}\) são iguais. O embedding ou incorporação é um procedimento padrão em análise de séries temporais. Com o embedding realizado, a análise subsequente depende do objetivo da investigação.

Segundo passo: Decomposição em valores singulares (SVD)

O segundo passo, a etapa de SVD, realiza a decomposição em valores singulares da matriz de trajetória e a representa como uma soma de matrizes elementares biortogonais de posto um.

Denotemos por \(\lambda_1,\cdots,\lambda_L\) os autovalores de \(\pmb{X}\pmb{X}^\top\) em ordem decrescente de magnitude \(\lambda_1\geq \cdots \geq \lambda_L\geq 0\) e por \(U_1,\cdots,U_L\) o sistema ortonormal, isto é, \(\big<U_i,U_j\big> = 0\) para \(i = j\), propriedade de ortogonalidade e \(||U_i|| = 1\), propriedade da norma unitária dos autovetores da matriz \(\pmb{X}\pmb{X}^\top\) correspondentes a esses autovalores. \(\big<U_i,U_j\big>\) é o produto interno dos vetores \(U_i\) e \(U_j\) e \(||U_i||\) é a norma do vetor \(U_i\).

Seja \[ d=\max\{i\, : \, \lambda>0\}=\mbox{posto}(\pmb{X})\cdot \] Se denotarmos \(V_i = \pmb{X}^\top U_i/\lambda_i\), então o SVD da matriz de trajetória pode ser escrita como: \[ \tag{1} \pmb{X}=\pmb{X}_1+\cdots+\pmb{X}_d, \] onde \(\pmb{X}_i=\sqrt{\lambda}_i \,U_iV_i^\top\), \(i=1,\cdots,d\).

As matrizes \(\pmb{X}_i\) têm posto 1; portanto, são matrizes elementares, \(U_i\) na literatura de SSA, são chamadas de ‘funções ortogonais empíricas fatoriais’ ou simplesmente EOFs e \(V_i\) frequentemente chamadas de ‘componentes principais’ representam os autovetores esquerdo e direito da matriz de trajetória.

O conjunto \(\big(\lambda_i,U_i,V_i\big)\) é chamado de \(i\)-ésimo autovalor triplo da matriz \(\pmb{X}\) ou autotriplo, \(\sqrt{\lambda}_i\), \(i = 1,\cdots,d\) são os valores singulares da matriz \(\pmb{X}\) e o conjunto \(\{\lambda_i\}\) é chamado de espectro da matriz \(\pmb{X}\). Se todos os autovalores tiverem multiplicidade um, então a expansão (1) é definida de forma única.

O SVD (1) é ótimo no sentido de que, entre todas as matrizes \(\pmb{X}^{(r)}\) de posto \(r < d\), a matriz \(\sum_{i=1}^r X_i\) fornece a melhor aproximação para a matriz de trajetória \(\pmb{X}\), de modo que \(||\pmb{X}-\pmb{X}^{(r)}||\) é mínimo.

Observe que \[ ||\pmb{X}||^2 = \sum_{i=1}^d \lambda_i \qquad \mbox{e} \qquad ||X_i|| =\lambda_i \qquad \mbox{para} \qquad i = 1,\cdots,d\cdot \]

Assim, consideramos a razão \(\lambda_i/\sum_{k=1}^d \lambda_k\) como a característica da contribuição da matriz \(\pmb{X}_i\) para a expansão (1). Consequentemente, \(\sum_{i=1}^r \lambda_i/\sum_{k=1}^d \lambda_k\), a soma das primeiras \(r\) razões, é a característica da aproximação ótima da matriz de trajetória pelas matrizes de posto \(r\).

2.2 Etapa 2: Reconstrução


Primeiro passo: Agrupamento

O passo de agrupamento corresponde à divisão das matrizes elementares \(\pmb{X}_i\) em vários grupos e à soma das matrizes dentro de cada grupo. Seja \(I=\{i_1,\cdots,i_p\}\) um grupo de índices \(i_1,\cdots,i_p\). Então, a matriz \(\pmb{X}_I\) correspondente ao grupo \(I\) é definida como \[ \pmb{X}_I = \pmb{X}_{i_1} + \cdots + \pmb{X}_{i_p}\cdot \]

A decomposição do conjunto de índices \(J = 1,\cdots, d\) em subconjuntos disjuntos \(I_1,\cdots, I_m\) corresponde à representação \[ \tag{2} \pmb{X}_I = \pmb{X}_{i_1} + \cdots + \pmb{X}_{i_p}\cdot \]

O procedimento de escolha dos conjuntos \(I_1,\cdots,I_m\) é chamado de agrupamento de autovalores e autovetores. Para um dado grupo \(I\), a contribuição do componente \(\pmb{X}_I\) na expansão (1) é medida pela proporção dos autovalores correspondentes: \[ \sum_{i\in I} \lambda_i/\sum_{k=1}^d \lambda_k\cdot \]

Segundo passo: Média diagonal

A média diagonal transforma cada matriz \(I\) em uma série temporal, que é um componente aditivo da série inicial \(Y_T\). Se \(z_{ij}\) representa um elemento de uma matriz \(\pmb{Z}\), então o \(k\)-ésimo termo da série resultante é obtido pela média de \(z_{ij}\) sobre todos os \(i\), \(j\) tais que \(i + j = k + 2\). Este procedimento é chamado de média diagonal, ou hankelização da matriz \(\pmb{Z}\).

O resultado da hankelização de uma matriz \(\pmb{Z}\) é a matriz de Hankel \(\mathcal{H}\pmb{Z}\), que é a matriz de trajetória correspondente à série obtida como resultado da média diagonal. Note que a hankelização é um procedimento ótimo, no sentido de que a matriz \(\mathcal{H}\pmb{Z}\) é a mais próxima de \(\pmb{Z}\), com relação à norma da matriz, dentre todas as matrizes de Hankel de tamanho correspondente, para mais informações, veja N. Golyandina, Nekrutkin, and Zhigljavsky (2001).

Por sua vez, a matriz de Hankel \(\mathcal{H}\pmb{Z}\) define unicamente a série, relacionando os valores nas diagonais aos valores na série. Aplicando o procedimento de hankelização a todos os componentes da matriz de (2), obtemos outra expansão: \[ \tag{3} \pmb{X}=\widetilde{\pmb{X}}_{I_1}+\cdots+\widetilde{\pmb{X}}_{I_m}, \] onde \(\widetilde{\pmb{X}}_{I_1}=\mathcal{H}\pmb{X}\).

Isso é equivalente à decomposição da série inicial \(Y_T = (y_1,\cdots, y_T )\) em uma soma de \(m\) séries: \[ \tag{4} y_t=\sum_{k=1}^m \widetilde{y}_t^{(k)}, \] onde \(\widetilde{Y}_T^{(k)}=\big(\widetilde{y}_1^{(k)},\cdots,\widetilde{y}_T^{(k)}\big)\) corresponde à matriz \(\pmb{X}_{I_k}\).

3 Aplicação


3.1 Série de óbitos


A série de óbitos USAccDeaths mostra as mortes acidentais mensais nos EUA entre 1973 e 1978. Esses dados foram utilizados por muitos autores (ver, por exemplo, Brockwell and Davis (2002)) e podem ser encontrados em diversas bibliotecas de dados de séries temporais.

data(USAccDeaths)
print(USAccDeaths)
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1973  9007  8106  8928  9137 10017 10826 11317 10744  9713  9938  9161  8927
## 1974  7750  6981  8038  8422  8714  9512 10120  9823  8743  9129  8710  8680
## 1975  8162  7306  8124  7870  9387  9556 10093  9620  8285  8466  8160  8034
## 1976  7717  7461  7767  7925  8623  8945 10078  9179  8037  8488  7874  8647
## 1977  7792  6957  7726  8106  8890  9299 10625  9302  8314  8850  8265  8796
## 1978  7836  6892  7791  8192  9115  9434 10484  9827  9110  9070  8633  9240
par(mar=c(3,3,1,1))
plot(USAccDeaths, main="Mortes acidentais mensais nos Estados Unidos (1973-1978)", 
     ylab="Mortes", xlab="Year")
grid()

Figura 1: Mortes acidentais mensais nos EUA (1973–1978).

Aplicamos a técnica SSA a esse conjunto de dados para ilustrar a capacidade da técnica SSA em extrair tendências, oscilações, ruídos e realizar previsões. Todos os resultados e figuras na aplicação a seguir foram obtidos por meio do pacote Rssa. A Figura 1 mostra a série de óbitos no período de 1973 a 1978.

3.2 Decomposição: comprimento da janela e SVD


Como mencionado anteriormente, o comprimento da janela \(L\) é o único parâmetro na etapa de decomposição. A seleção do comprimento de janela adequado depende do problema em questão e de informações preliminares sobre a série temporal.

Resultados teóricos indicam que \(L\) deve ser suficientemente grande, mas não maior que \(T/2\). Além disso, se soubermos que a série temporal pode ter um componente periódico com período inteiro, por exemplo, se esse componente for sazonal, então, para obter uma melhor separabilidade desse componente periódico, é aconselhável utilizar um comprimento de janela proporcional a esse período.

Usando essas recomendações, consideramos \(L = 24\). Assim, com base nesse comprimento de janela e na SVD da matriz de trajetória (\(24\times 24\)), temos 24 autovetores triplos, ordenados por sua contribuição (participação) na decomposição.

Observe que as linhas e colunas da matriz de trajetória \(\pmb{X}\) são subséries da série temporal original. Portanto, os autovetores à esquerda \(U_i\) e os componentes principais \(V_i\) (autovetores à direita) também possuem uma estrutura temporal e, consequentemente, podem ser considerados séries temporais. Consideremos o resultado da etapa de SVD. A Figura 2 representa os componentes principais relacionados aos 12 primeiros autovetores triplos.

library(Rssa)
s <- ssa(USAccDeaths, L = 24)
plot(s, type = "vectors", idx = 1:12)

Figura 2: Componentes principais relacionados aos primeiros 12 autovalores e autovetores.

3.3 Informações Suplementares


Vamos descrever algumas informações que se mostram muito úteis na identificação dos autovalores e autovetores da SVD da matriz de trajetória da série original.

As informações suplementares nos ajudam a formar os grupos adequados para extrair a tendência, os componentes harmônicos e o ruído. Assim, as informações suplementares podem ser consideradas uma ponte entre as etapas de decomposição e reconstrução: \[ \mbox{Decomposição} \longmapsto \mbox{Informações Suplementares} \longmapsto \mbox{Reconstrução} \]

A seguir, explicamos brevemente alguns métodos úteis na separação do componente de sinal do ruído.

Informações auxiliares

A disponibilidade de informações auxiliares em muitas situações práticas aumenta a capacidade de construir o modelo adequado. Certamente, as informações auxiliares sobre a série inicial sempre tornam a situação mais clara e auxiliam na escolha dos parâmetros dos modelos.

Essas informações não apenas nos ajudam a selecionar o grupo adequado, mas também são úteis para previsão e detecção de pontos de mudança com base na técnica SSA. Por exemplo, a suposição de que existe uma periodicidade anual na série de óbitos sugere que devemos prestar atenção à frequência \(k/12\), (\(k=1,\cdots,12\)). Obviamente, também podemos usar as informações auxiliares para selecionar o comprimento de janela adequado.

Valores singulares

Normalmente, cada componente harmônica com uma frequência diferente produz dois autovalores e autotriplos com valores singulares próximos, exceto para a frequência 0.5, que fornece um autovalor e autotriplo com vetor singular em forma de dente de serra. Isso ficará mais claro se \(N\), \(L\) e \(K\) forem suficientemente grandes.

Outra informação útil é fornecida pela verificação de quebras nos espectros de autovalores. Como regra geral, uma série de ruído puro produz uma sequência de valores singulares que diminui lentamente. Portanto, platôs explícitos nos espectros de autovalores indicam os números ordinais dos autovalores e autotriplos emparelhados. A Figura 3 mostra o gráfico dos logaritmos dos 24 valores singulares para a série USAccDeath.

plot(s, type = "values")

Figura 3: Logaritmos dos 24 autovalores.

Cinco pares evidentes com valores singulares principais quase iguais correspondem a cinco componentes (quase) harmônicos da série: os pares de autovalores triplos 2-3, 4-5, 7-8, 9-10 e 11-12 estão relacionados a harmônicos com períodos específicos (mostraremos mais adiante que eles correspondem aos períodos 12, 6, 2.5, 4 e 3).

Diagramas de dispersão aos pares

Na prática, os valores singulares de dois autotriplos de uma série harmônica são frequentemente muito próximos entre si, e esse fato simplifica a identificação visual dos componentes harmônicos. Uma análise dos diagramas de dispersão aos pares dos vetores singulares permite identificar visualmente os autovetores e autovetores que correspondem aos componentes harmônicos da série, desde que esses componentes sejam separáveis do componente residual.

Considere uma harmônica pura com frequência \(\omega\), fase e amplitude definidas, e uma situação ideal onde \(P = 1/\omega\) é um divisor do comprimento da janela \(L\) e \(K\). Como \(P\) é um número inteiro, ele representa o período da harmônica. Na situação ideal, os autovetores à esquerda e os componentes principais têm a forma de sequências de seno e cosseno com o mesmo \(P\) e a mesma fase. Assim, a identificação dos componentes gerados por uma harmônica se reduz à determinação desses pares.

As funções seno e cosseno puras, com frequências, amplitudes e fases iguais, criam um diagrama de dispersão com os pontos dispostos em um círculo. Se \(P = 1/\omega\) for um número inteiro, então esses pontos são os vértices de um polígono regular de \(P\) vértices.

Para a frequência racional \(\omega = m/n < 0.5\), com \(m\) e \(n\) números inteiros primos entre si, os pontos são os vértices dos diagramas de dispersão de um polígono regular de \(n\) vértices. A Figura 4 mostra diagramas de dispersão de 6 pares de sequências seno/cosseno (sem ruído) com fase zero, mesma amplitude e períodos 12, 6, 4, 3, 2,5 e 2,4.

Figura 4: Diagramas teóricos de dispersão dos 6 pares de senos/cossenos.

A Figura 5 apresenta diagramas de dispersão dos autovetores emparelhados na série de USAccDeath, correspondentes aos harmônicos com períodos 12, 6, 4, 3 e 2.5. Eles estão ordenados por sua contribuição (participação) na etapa de SVD.

plot(s, type = "paired", idx = c(2,4,9)) 

plot(s, type = "paired", idx = c(11,7))

Figura 5: Diagramas de dispersão dos autovetores harmônicos emparelhados.

3.4 Separabilidade


O principal conceito no estudo das propriedades da SSA é a ‘separabilidade’, que caracteriza o quão bem os diferentes componentes podem ser separados uns dos outros. A decomposição SSA da série \(Y_T\) somente pode ser considerada bem-sucedida se os componentes aditivos resultantes da série forem aproximadamente separáveis entre si.

A seguinte quantidade, chamada de correlação ponderada ou \(\omega\)-correlação, é uma medida natural da dependência entre duas séries \(Y_T^{(1)}\) e \(Y_T^{(2)}\): \[ \rho_{12}^{(\omega)}=\dfrac{\big<Y_T^{(1)},Y_T^{(2)}\big>_\omega}{||Y_T^{(1)}||_\omega ||Y_T^{(2)}||_\omega}, \] onde \[ ||Y_T^{(i)}||=\sqrt{\big<Y_T^{(i)},Y_T^{(i)}\big>_\omega}, \] \[ \big<Y_T^{(i)},Y_T^{(j)}\big>_\omega=\sum_{k=1}^T \omega_k y_T^{(i)} y_T^{(j)}, \] \(i,j=1,2\), \(\omega_k =\min\{k, L, T-k\}\), aqui assumimos \(L\leq T/2\).

Uma dica natural para o agrupamento é a matriz dos valores absolutos das \(\omega\)-correlações, correspondente à decomposição completa, nesta decomposição, cada grupo corresponde a apenas um componente da matriz da SVD).

Se o valor absoluto das \(\omega\)-correlações for pequeno, as séries correspondentes são quase \(\omega\)-ortogonais, mas, se for grande, as duas séries estão longe de serem \(\omega\)-ortogonais e, portanto, são pouco separáveis. Assim, se dois componentes reconstruídos tiverem \(\omega\)-correlação zero, significa que esses dois componentes são separáveis.

Valores altos de \(\omega\)-correlações entre componentes reconstruídos indicam que os componentes possivelmente devem ser agrupados e correspondem ao mesmo componente na decomposição SSA.

w <- wcor(s, groups = 1:12)
plot(w, grid = c(2,4, 5,7))

Figura 6: Matriz de \(\omega\)-correlações para 12 componentes reconstruídos.

A Figura 6 acima mostra as \(\omega\)-correlações para 12 componentes reconstruídos em uma escala de cinza, do branco ao preto, correspondendo aos valores absolutos das correlações de 0 a 1.

4 Reconstrução: agrupamento e média diagonal


A reconstrução é a segunda etapa da técnica SSA. Como mencionado anteriormente, esta etapa inclui dois passos distintos: agrupamento, a identificação do componente de sinal e do ruído e média diagonal, utilização dos autovalores agrupados para reconstruir a nova série sem ruído.

Normalmente, o autovalor principal descreve a tendência geral da série. Como na maioria dos casos os autovalores com menor participação estão relacionados ao componente de ruído da série, precisamos identificar o conjunto de autovalores principais.

4.1 Agrupamento: tendência, harmônicos e ruído


Identificação de tendência:

A tendência é o componente de variação lenta de uma série temporal que não contém componentes oscilatórios. Suponha que a própria série temporal seja um componente desse tipo. A prática mostra que, nesse caso, um ou mais dos autovetores principais também variarão lentamente.

Sabemos que os autovetores têm, em geral, a mesma forma que os componentes correspondentes da série temporal inicial. Portanto, devemos encontrar autovetores de variação lenta. Isso pode ser feito considerando gráficos unidimensionais dos autovetores, veja a Figura 2.

Em nosso caso, o autovetor principal tem definitivamente a forma desejada. A Figura 7 mostra a tendência extraída sobre o fundo da série original, obtida a partir do primeiro autotriplo. Observe que podemos construir uma aproximação mais complexa da tendência se usarmos outros autotriplos. No entanto, a precisão que obteríamos seria muito pequena e o modelo da tendência se tornaria muito mais complicado.

par(mar=c(3,3,1,1))
plot(USAccDeaths, main="Mortes acidentais mensais nos Estados Unidos (1973-1978)", 
     ylab="Mortes", xlab="Year")
grid()
s1 <- reconstruct(s, groups = list(Trend = 1))
lines(s1$Trend, col = "red", lty = 1, lwd = 2)

Figura 7: Extração de tendência, primeiro autotriplo.

A Figura 8 mostra a tendência extraída, obtida a partir do primeiro e do sexto autotriplos. Observa-se que a análise conjunta do primeiro e do sexto autotriplos revela a tendência geral da série de forma mais precisa do que a análise isolada do primeiro autotriplo. Contudo, o sexto autotriplo não pertence completamente à componente de tendência, podendo ser considerado uma combinação da tendência com a componente harmônica.

par(mar=c(3,3,1,1))
plot(USAccDeaths, main="Mortes acidentais mensais nos Estados Unidos (1973-1978)", 
     ylab="Mortes", xlab="Year")
grid()
s2 <- reconstruct(s, groups = list(Trend = c(1,6)))
lines(s2$Trend, col = "red", lty = 1, lwd = 2)

Figura 8: Extração de tendência, primeiro e sexto autotriplos.

Identificação harmônica:

O problema geral aqui é a identificação e separação dos componentes oscilatórios da série que não fazem parte da tendência. A formulação do problema SSA é especificada principalmente pela natureza independente de modelo do método.

A escolha de \(L = 24\) permite extrair simultaneamente todos os componentes sazonais: 12, 6, 4, 3 e 2.5 meses, bem como a tendência. A Figura 9 mostra a oscilação da nossa série obtida pelos autotriplos 2 a 12. Comparando a Figura 9 com a Figura 1, fica claro que os autotriplos selecionados para identificar os componentes harmônicos foram escolhidos corretamente.

par(mar=c(3,3,3,1))
s3 <- reconstruct(s, groups = list(Trend = 1, Seasonal = 2:12))
plot(s3$Seasonal, main="Mortes acidentais mensais nos Estados Unidos (1973-1978) \n componente harmônico", 
     ylab="Mortes", xlab="Year", col = "darkred", lty = 1, lwd = 2)
grid()

Figura 9: Extração de oscilações: autotriplos 2 a 12.

A Figura 10 mostra a oscilação da nossa série obtida pelos autotriplos 2 a 5 e 7 a 12. Neste caso, consideramos o sexto autotriplo como um componente da tendência. Parece não haver grande discrepância entre selecionar o sexto autotriplo nos componentes de tendência ou oscilação, como se pode observar nas Figuras 9 e 10.

par(mar=c(3,3,3,1))
s4 <- reconstruct(s, groups = list(Trend = 1, Seasonal = c(2:5,7:12)))
plot(s4$Seasonal, main="Mortes acidentais mensais nos Estados Unidos (1973-1978) \n componente harmônico", 
     ylab="Mortes", xlab="Year", col = "cyan", lty = 1, lwd = 2)
grid()

Figura 10: Extração de oscilações: autotriplos 2–5 e 7–12.

Detecção de ruído:

O problema de encontrar uma estrutura refinada de uma série por meio do SSA equivale à identificação dos autovalores e autovetores da decomposição em valores singulares (SVD) da matriz de trajetória dessa série, que correspondem à tendência, a vários componentes oscilatórios e ao ruído.

Do ponto de vista prático, uma maneira natural de extrair ruído é agrupar os autovalores e autovetores que aparentemente não contêm elementos de tendência e oscilações. Vamos discutir o autovalor e autovetor 13. Consideramos que ele pertence ao ruído porque o período do componente reconstruído pelo autovalor e autovetor 13 é uma mistura dos períodos 3, 10, 14 e 24, e como o periodograma indica, isso não pode ser interpretado no contexto da sazonalidade para essa série. Assim, classificaremos o autovalor e autovetor 13 como parte do ruído. A Figura 11 mostra os resíduos obtidos pelos autovalores e autovetores 13–24.

par(mar=c(3,3,3,1))
s5 <- reconstruct(s, groups = list(Trend = 1, Seasonal = 2:12, Noise = 13:24))
plot(s5$Noise, main="Mortes acidentais mensais nos Estados Unidos (1973-1978) \n componente ruído", 
     ylab="Mortes", xlab="Year", col = "darkcyan", lty = 1, lwd = 2)
grid()

Figura 11: Séries residuais: autotriplos 13–24.

Média diagonal

A última etapa da técnica SSA é a média diagonal. A média diagonal aplicada a uma matriz resultante \(X_{I_k}\), obtida na etapa de agrupamento para o \(k\)-ésimo grupo de \(m\), produz \(\widetilde{Y}_T^k = \big(\widetilde{y}_1^k,\cdots,\widetilde{y}_T^k\big)\) e, portanto, a série inicial \(Y_T = (y_1,\cdots,y_T)\) é decomposta na soma de \(m\) séries \[ y_t = \sum_{k=1}^m \widetilde{y}_t^k, \qquad 1\leq t\leq T\cdot \]

Se considerarmos apenas a tendência (autotriplos 1 ou 1 e 6), o componente harmônico (autotriplos 2–12 ou (2–5,7–12)) e o ruído (autotriplos 13–24) como grupos, então temos 3 grupos (\(m = 3\)). No entanto, podemos ter 8 grupos se considerarmos cada grupo em detalhe, como por exemplo: autotriplos 1, 2–3, 4–5, 6, 7–8, 9–10, 11–12 (que correspondem ao sinal) e 13–24, ou 7 grupos se combinarmos os autotriplos 1 e 6 em um único grupo. A Figura 12 mostra o resultado da extração ou reconstrução do sinal sem ruído, obtido a partir dos autotriplos 1–12.

par(mar=c(3,3,3,1))
s5 <- reconstruct(s, groups = list(Trend = 1, Seasonal = 2:12, Noise = 13:24))
plot(USAccDeaths, main="Mortes acidentais mensais nos Estados Unidos (1973-1978)", 
     ylab="Mortes", xlab="Year", col = "black", lty = 1, lwd = 2)
lines(s5$Trend+s5$Seasonal, col = "darkred", lty = 2, lwd = 2)
grid()

Figura 12: Série reconstruída (autotriplos 1–12).

A linha tracejada e a linha contínua correspondem, respectivamente, à série reconstruída e à série original. Como indicado nesta figura, os grupos considerados para a reconstrução da série original são ótimos (lembrando que a etapa de SVD possui propriedades ótimas). Se somarmos as séries das Figuras 7 e 9 (ou 8 e 10), obteremos a série refinada na Figura 12.

5 Previsão


A previsão por meio de SSA pode ser aplicada a séries temporais que satisfazem aproximadamente fórmulas lineares recorrentes (LRF). Diremos que a série \(Y_T\) satisfaz uma LRF de ordem \(d\) se houver números \[ y_{i+d}=\sum_{k=1}^d a_k \, y_{i+d-k}, \qquad 1\leq i\leq T-d\cdot \]

A classe de séries regidas por fórmulas lineares recorrentes (FLRs) é bastante ampla; ela inclui séries harmônicas, polinomiais e exponenciais, e é fechada sob adição e multiplicação termo a termo. As FLRs também são importantes para aplicações práticas.

Para encontrar os coeficientes \(a_1,\cdots,a_d\), podemos usar os autovetores obtidos na etapa de decomposição em valores singulares (SVD) ou o polinômio característico, para mais informações, consulte N. Golyandina, Nekrutkin, and Zhigljavsky (2001).

Suponha que temos uma série \(Y_T = Y_T^{(1)} + Y_T^{(2)}\) e o problema de prever seu componente \(Y_T^{(1)}\). Se \(Y_T^{(2)}\) puder ser considerado ruído, então o problema consiste em prever o sinal \(Y_T^{(1)}\) na presença de um ruído \(Y_T^{(2)}\).

As principais hipóteses são:

  1. a série \(Y_T^{(1)}\) admite uma continuação recorrente com o auxílio de uma função de resposta linear (LRF) de dimensão relativamente pequena \(d\), e

  2. existe um número \(L\) tal que as séries \(Y_T {(1)}\) e \(Y_T^{(2)}\) são aproximadamente separáveis para o comprimento da janela \(L\).

A Figura 13 mostra a série original (linha contínua), a série reconstruída (linha pontilhada) e sua previsão após 1978, os seis pontos de dados de 1979. A linha pontilhada vertical mostra a truncagem entre o último ponto da série original e o ponto inicial da previsão. A Figura 13 mostra que a série reconstruída, obtida a partir dos autovalores 1-12, e a série original estão próximas, indicando que os valores previstos são razoavelmente precisos.

previstos <- forecast(s, groups = list(1:12), 
                      method = "recurrent", interval = "confidence", 
                      only.intervals = FALSE, len = 6, R = 100, level = 0.95)
nova.serie = ts(c(USAccDeaths,previstos$mean), start = c(1973,1), frequency = 12)
par(mar=c(3,3,3,1))
plot(nova.serie, main="Mortes acidentais mensais nos Estados Unidos (1973-1978)", 
     ylab="Mortes", xlab="Year", col = "black", lty = 1, lwd = 2)
lines(s5$Trend+s5$Seasonal, col = "darkred", lty = 2, lwd = 2)
abline(v = tail(time(USAccDeaths),1), lty = 2)
points(previstos$mean, pch = 19)
lines(previstos$lower, col = "green"); lines(previstos$upper, col = "green")
polygon(c(time(previstos$mean), rev(time(previstos$mean))), 
        c(previstos$lower, rev(previstos$upper)), col = rgb(0.1, 0.5, 0.8, 0.3), border = NA)
grid()

Figura 13: Série original (linha contínua), série reconstruída (linha pontilhada) e os 6 pontos de dados previstos para 1979.

6 Referências


Alexandrov, Th., and N. Golyandina. 2004a. “The Automatic Extraction of Time Series Trend and Periodical Components with the Help of the Caterpillar SSA Approach.” Exponenta Pro, no. 3-4: 54–61.(In Russian.).
———. 2004b. “Thresholds Setting for Automatic Extraction of Time Series Trend and Periodical Components with the Help of the Caterpillar SSA Approach.” Proc. IV International Conference SICPRO’, no. 05: 25–28.
Allen, M. R., and L. A. Smith. 1996. “Monte Carlo SSA: Detecting Irregular Oscillations in the Presence of Coloured Noise.” Journal of Climate, no. 9: 3373–3404.
Brockwell, P. L., and R. A. Davis. 2002. Introduction to Time Series and Forecasting. 2nd edition. Springer.
Broomhead, D. S., and G. P. King. 1986. “Extracting Qualitative Dynamics from Experimental Data.” Physica D, no. 20: 217–36.
Danilov, D., and A. Zhigljavsky. 1997. Principal Components of Time Series: The “Caterpillar” Method. University of St. Petersburg Press. (In Russian.).
Elsner, J. B., and A. A. Tsonis. 1996. Singular Spectral Analysis. A New Tool in Time Series Analysis. Plenum Press.
Ghil, M., and C. Taricco. 1997. “Advanced Spectral Analysis Methods. In Past and Present Variability of the Solar-Terrestrial System: Measurement.” In Data Analysis and Theoretical Model, 137–59. IOS Press.
Golyandina, Nina, and E Osipov. 2007. “The ‘Caterpillar’-SSA Method for Analysis of Time Series with Missing Values.” Journal of Statistical Planning and Inference 137 (8): 2642–53.
Golyandina, N., V. Nekrutkin, and A. Zhigljavsky. 2001. Analysis of Time Series Structure: SSA and Related Techniques. Chapman & Hall/CRC.
Kondrashov, D., Y. Feliks, and M. Ghil. 2005. “Oscillatory Modes of Extended Nile River Records (a. D. 622–922).” Geophysical Research Letters, no. 32: L10702, doi:10.1029/2004GL022156.
Kondrashov, D., and M. Ghil1. 2006. “Spatio-Temporal FIlling of Missing Points in Geophysical Data Sets.” Nonlinear Processes in Geophysics, no. 13: 151–59.
Moskvina, V. G., and A. Zhigljavsky. 2003. “An Algorithm Based on Singular Spectrum Analysis for Change-Point Detection.” Communication in Statistics - Simulation and Computation, no. 32: 319–52.
Schoellhamer, D. H. 2001. “Singular Spectrum Analysis for Time Series with Missing Data.” Geophysical Research Letters, no. 28: 3187–90.
Vautard, R., P. Yiou, and M. Ghil. 1992. “Singular-Spectrum Analysis: A Toolkit for Short, Noisy Chaotic Signal.” Physica D, no. 58: 95–126.
Yiou, P., D. Sornette, and M. Gill. 2000. “Data-Adaptive Wavelets and Multi-Scale Singular Spectrum Analysis.” Physica D, no. 142: 254–90.