MÉTODOS OBSERVACIONAIS EM CLIMATOLOGIA E METEOROLOGIA DE MESOESCALA : NOTAS DE AULA Prof. Resposável: Dra. Leila M. Véspoli de Carvalho IAG/USP ANÁLISE DE SÉRIES TEMPORAIS I) ALGORITMOS PARA A REMOÇÃO DO CICLO ANUAL, SEMI-ANUAL E TENDÊNCIAS Referências Básicas : 1Chatfield C., 1996: The Analysis of Time Series: An introduction. Chapman & Hall, fifth edition, NY. 283 pp Jenkins, G. M. and D. G. Watts, 1968: Spectral Analysis and its Applications. Holden-day, SF, 525pp. 3 Wilks, D. S., 1995: Statistical methods in the Atmospheric Sciences. Academic Press, NY, 468 pp. 2 1) A presença de ciclos nas séries temporais de dados meteorológicos e como tratá-los Motivação: Uma série temporal pode ser definida como um conjunto de observações feitas seqüencialmente no tempo. A forma de analisá-la é o objeto de estudo desse curso. Assim, para que seja percebido o exato significado do que pretendemos fazer daqui para frente, vamos observar a importância de ciclos nas variáveis meteorológicas que estudamos. Onde conseguir dados meteorológicos: http://www.cdc.noaa.gov/cdc/data.interp_OLR.html http://www.cdc.noaa.gov/Timeseries/ http://www.ncdc.noaa.gov/oa/climate/research/2003/may/global.html http://wesley.wwb.noaa.gov/ncep_data/index.html DICAS DOS PROCEDIMENTOS PARA SE OBTER SERIES TEMPORAIS A PARTIR DE REANALYSIS-2 DO NCEP ENCONTRAM-SE AQUI (ARQUIVO .pdf): Terminologias e conceitos importantes: 1. Séries temporais e processos estocásticos: 1.1 Funções determinísticas e não-determinísticas (JW, CH) Uma série temporal pode ser uma função x aleatória ou não-determinística de uma variável independente t. Na maioria das situações, a função x(t) será uma função do tempo, mas em outras situações pode ser uma função de outro parâmetro físico, como por exemplo, do espaço. Uma característica das séries temporais é que seu comportamento futuro não pode ser previsto exatamente, como seria o caso de uma função ‘determinística’ do tempo. Exercício -1: Descreva um exemplo de uma função determinística, tendo como variável independente o tempo ou o espaço ou ambos. Desafio-1. Para acompanhar esse exercício, você precisara ler as discussões de Tsonis and Elsner(1989). Considere o sistema de 3 equações diferenciais que descrevem a convecção segundo o modelo de Lorenz (Lorenz 1963). Discuta o que são atratores estranhos e sua implicação para a previsão de tempo. Contudo, muitos fenômenos naturais, embora possuam oscilações no tempo, não podem ser descritos por funções determinísticas. Por exemplo, podemos medir a temperatura do ar em um abrigo meteorológico todos os dias e verificarmos a presença de um ciclo diurno. Entretanto, não conseguimos sempre determinar uma relação determinística que possa ser ajustada a cada intervalo dessa série de dados porque diversos fatores podem estar causando variações nessa medida (exemplo, nebulosidade, entradas de frentes, alteração dos ventos por circulações locais, etc.). Se 1 compararmos uma série temporal de temperatura em um determinado sítio em dois anos distintos, podemos verificar visualmente que esses dois trechos da série não se parecem um com outro. Essa observação leva a noção de PROCESSO ESTOCÁSTICO. 1.2 Processos Estocásticos (JW, WI) Uma vez que diferentes secções de uma série temporal se parecem uma com a outra apenas nas suas propriedades médias, é necessário descrever essas séries por leis de probabilidades ou modelos. Assim, os valores possíveis das séries temporais a um dado tempo t são descritos por uma VARIÁVEL ALEATÓRIA X(t) e sua associada DISTRIBUIÇÃO DE PROBABILIDADES. O valor observado x(t) da série temporal no tempo t é então considerado como um dos infinitos valores nos quais a variável X(t) pode ter no tempo t. Em outras palavras, o comportamento da série temporal para todos os tempos t pode ser descrito por um conjunto de variáveis aleatórias {X(t)} onde t pode ter qualquer valor entre -∞ a +∞. Assim, as propriedades estatísticas das séries são descritas por distribuições de probabilidade com qualquer conjunto de tempos t1, t2, ..., tN . O conjunto ordenado de variáveis aleatórias {X(t)} em associação com sua distribuição de probabilidades é chamado de PROCESSO ESTOCÁSTICO. Exercício -2 Veja outras interpretações em: http://pespmc1.vub.ac.be/ASC/STOCHA_PROCE.html e http://www.wikipedia.org/wiki/Stochastic_process Discuta definições e exemplos meteorológicos de: variável aleatória; função distribuição de probabilidade; processos estocásticos. Usem as demais referências se necessário. Desafio-2. Considere as series temporais de temperatura fornecidas no curso. Demonstre, com argumentos baseados no que você aprendeu no curso e leu nos sites acima, se e porque a temperatura pode ser considerada uma variável aleatória e sua serie temporal um processo estocástico. Terminologias: • Séries temporais contínuas: Medidas no tempo contínuas (Ex: medidas de temperatura a partir de um termógrafo) • Séries temporais discretas: Observações tomadas em intervalos de tempo específicos, usualmente igualmente espaçados (Ex: temperatura média mensal). 2. Conceitos importantes em séries temporais: Vimos que os processos estocásticos, a partir dos quais considera-se que a série temporal observada foi gerada, podem ser descritos por distribuições de probabilidades associadas com todos os possíveis conjuntos de pontos no tempo. Para inferir a natureza dessas distribuições de probabilidade a partir de uma única ou pequeno número de séries é um exercício impossível ou de pouca significância prática. Vamos discutir algumas das mais importantes simplificações que podem ser feitas. As suposições mais importantes feitas sobre uma série temporal são: a) o correspondente processo estocástico é ESTACIONÁRIO; b) um processo estocástico estacionário pode ser adequadamente descrito pelos mais BAIXOS MOMENTOS (ou momentos de baixa ordem) de suas distribuições de probabilidade. Esses momentos de baixa ordem incluem: MÉDIA, VARIÂNCIA, COVARIÂNCIA e a TRANSFORMADA DE FOURIER DA FUNÇÃO DE COVARIÂNCIA, O ESPECTRO DE POTÊNCIA. Assim, uma aproximação alternativa é supor que o processo estocástico pode ser adequadamente descrito por meio de um modelo contendo uns poucos parâmetros os quais podem ser estimados a partir dos dados. Vamos discutir essas simplificações. 2.1. Estacionaridade (CH; WI; JW): Suponha que você esteja examinando uma série temporal por um certo tempo limitado, por exemplo, a saída obtida de um gerador de ruídos. Suponha que a comparação de diferentes trechos dessa série mostra que diferentes secções são ‘parecidas’. Em contraste, quando você observa a concentração de CO2 global nos últimos 100 anos ou a área de florestas desmatadas do planeta no mesmo período, vai notar que existe uma tendência dessas séries temporais de crescerem com o tempo. Assim, diferentes secções dessas séries possuem características distintas. A saída do gerador de ruído é considerada um processo ESTACIONÁRIO enquanto as séries temporais do CO2 e de desmatamento são ditas NÃOESTACIONÁRIAS. Qualitativamente, uma série estacionária é aquela que está em EQUILÍBRIO ESTATÍSTICO, no sentido que contém NENHUMA TENDÊNCIA, enquanto que uma série não-estacionária é aquela cujas propriedades mudam com o tempo. Na prática, as séries são usualmente de 3 tipos: aquelas que exibem propriedades de estacionaridade em longo período, como, por exemplo, as saídas de geradores de ruído. Aquelas que possuem uma razoável estacionaridade em períodos curtos, por exemplo, medidas de turbulência na atmosfera; e séries que são obviamente não estacionárias, no sentido que suas propriedades estão continuamente mudando com o tempo. Exemplos óbvios de nãoestacionaridade: temperatura em altas e médias latitudes, ventos (apresentam ciclos diurnos e anuais). Usualmente, o termo estacionaridade é interpretado como ‘fraca estacionaridade’ ou ‘ estacionaridade da covariância’. Neste sentido, estacionaridade implica que a média e a função de auto-correlação de uma série de dados não muda com o tempo. Diferentes pedaços de uma série de dados estacionária (por exemplo, os dados observados hoje e no futuro) podem ser considerados como TENDO UMA MESMA MÉDIA E VARIÂNCIA. Além disso, uma correlação entre variáveis em uma série estacionária é determinada apenas pela sua separação no tempo (ou seja, pelo seu “lag k”) e não pela sua absoluta posição no tempo. Isso significa que valores individuais em distintas porções da série podem ser diferentes embora essas duas porções da série se pareçam. A ESTACIONARIDADE DE COVARIÂNCIA é uma suposição menos restritiva que “estacionaridade restrita”, a qual implica que a distribuição total das variáveis na serie não muda com o tempo. A maior parte dos métodos que trata com não-estacionaridade de séries temporais está baseada em técnicas para remover ou filtrar a parte não-estacionária, deixando apenas a parte que pode ser tratada como estacionária. Em climatologia, utilizamos esse tipo de técnica quando 2 desejamos conhecer o comportamento das anomalias de uma determinada variável. Existem duas aproximações para tratar-se com séries nãoestacionárias. Ambas objetivam processar os dados de forma que permitam que uma subseqüente estacionaridade seja assumida. Por exemplo: subtração de uma função periódica média a partir dos dados sujeitos a um ciclo anual produziria uma nova série transformada com média constante igual a zero. A fim de produzir uma série com média e variância constante, seria necessário transformar essas anomalias em anomalias normalizadas: z= x − x x′ = (2.1) sx sx Onde z é a anomalia padronizada, calculada simplesmente pela subtração da média da amostra (que no caso seria igual a zero após remoção do ciclo anual) e dividindo pelo respectivo desvio padrão Sx, o qual varia. Por exemplo, não apenas as temperaturas tendem a ser mais frias durante o inverno, mas sua variabilidade tende a ser mais alta em regiões de latitudes médias. Uma aproximação possível para transformar séries de temperaturas mensais em uma série (aproximadamente) estacionária seria calcular as 12 médias mensais e os 12 desvios-padrão e então aplicar a Eq. (2.1) usando diferentes médias e desvios-padrão para o mês do calendário apropriado. Uma alternativa seria a estratificação dos dados. Isto é, poder-se-ia conduzir análises separadas de subconjuntos dos dados que são curtas o suficiente para serem consideradas aproximadamente estacionárias. Por exemplo, poder-se-ia analisar observações diárias para todas os dados disponíveis de janeiro para uma dada localização, assumindo-se que cada conjunto de 31 dias de dados é uma amostra que sofreu os mesmos processos físicos. Não necessariamente os processos seriam os mesmos para julho, ou fevereiro. Um exemplo sobre o uso da Eq. 2.1 para expressar dados climáticos em termos de anomalias padronizadas é o cálculo do índice do El-NiñoOscilação Sul(ENSO). Os valores do índice (veja figura abaixo) são derivados a partir de diferenças mensais nas anomalias padronizadas da pressão ao nível do mar em duas localizações: Tahiti, no Pacifico central; e Darwin, no Norte da Austrália. Assim, em termos da Eq. 2.1 o primeiro passo para calcular os pontos do gráfico é calcular a diferença Δz = z Tahiti − z Darwin para cada mês durante o período de anos considerados. Na figura abaixo, estão mostrados o IOS desde 1993 para uma ilustração mais clara. A anomalia padronizada ZTahiti para janeiro de 1997, por exemplo, é calculada subtraindo-se a pressão média para todos os janeiros em Taiti da pressão observada em janeiro de 1997. Esta diferença é então dividida pelo desvio-padrão, caracterizando uma variação ano-a-ano das pressões atmosféricas em Taiti. Na realidade, os valores do índice mostrado na figura são, em geral, anomalias da diferença das anomalias. Assim, a Eq. 2.1 é aplicada 2 vezes ao conjunto original de dados. A primeira das duas padronizações é tomada para minimizar a influência das mudanças sazonais na média mensal das pressões e a variabilidade ano-a-ano destas médias. A segunda padronização calcula a anomalia padronizada da diferença zTahiti – zDarwin e garante que o índice resultante terá unidade de desvio-padrão. Obtém-se essa padronização considerando-se uma nova variável ZΔz = série temporal de todos os Δz. Calcula-se a média dos Δz e seu desvio padrão. Existem algumas variações no jeito de se calcular o IOS. No caso do exemplo, extraído de http://www.bom.gov.au/climate/current/soi2.shtml o IOS está multiplicado por 10. A linha vermelha mostra a média móvel de 5 meses (COMO SE DEFINI UMA MÉDIA MÓVEL DE 5 MESES?) A interpretação física do IOS é que durante eventos El Niño o centro da precipitação no Pacífico tropical muda do Pacífico Oeste (próximo a Darwin) para leste ou o Pacífico Central (próximo a Taiti). Esta mudança está associada com as pressões à superfície acima da média em Darwin e mais baixas que a média em Taiti, o que junto produz um valor negativo no índice. Eventos excepcionalmente fortes (como o 82/83 e 97/98) produziram valores bem baixos no índice. 2.2 Análise de séries que contém tendência (CH) Tendência pode ser definida grosseiramente como “uma mudança de longo-termo no nível médio”. As tendências mais simples conhecida como “TENDÊNCIA LINEAR MAIS RUÍDO”, para qual a observação no tempo t é uma variável aleatória Xt dada por: 3 Xt=α+βt + εt (2.2) Onde α, β são constantes e εt representa o termo de erro aleatório com média zero (conhecido também como ruído branco). O nível médio no tempo t é dado por mt=(α+βt), o que é algumas vezes chamado de “TERMO DE TENDÊNCIA”. Alguns autores preferem descrever o coeficiente angular β como a tendência. Em outras palavras, a tendência é a MUDANÇA NO NÍVEL MÉDIO por unidade de tempo. A tendência na Eq. 2.2 é uma função determinística do tempo e é algumas vezes chamada de “TENDÊNCIA GLOBAL”. É geralmente não realística, e por isso existe agora uma maior ênfase em TENDÊNCIAS LOCAIS onde os parâmetros α e β variam no tempo. As tendências também podem ser NÃO-LINEARES. Um crescimento exponencial, ou uma tendência quadrática são alguns exemplos. Assim, tem-se que ter em mente que a análise de uma série temporal que exibe tendência depende (1) se o pesquisador quer exatamente medir essa tendência ou (2) se o pesquisador quer remover a tendência de forma a analisar flutuações locais. Isto também depende se o dado exibe SAZONALIDADE. Com dados contendo sazonalidade, é uma boa idéia começar calculando médias anuais sucessivas porque estas podem fornecer uma descrição simples das tendências implícitas. 2.2.1 Ajustando uma curva (CH) Um método comum de se tratar com dados não-sazonais que contém uma tendência, particularmente em dados diários, é ajustar uma função simples como uma curva polinomial (linear, quadrática, etc.), ou uma curva de Gompertz ou uma curva logística. Estas últimas têm sido empregadas em economia e biologia (crescimento de populações). A curva de Gompertz é dada por: log xt = a + brt (2.3) onde a, b, r são parâmetros com 0 < r <1, enquanto a curva logística é dada por xt=a/(1+be-ct) (2.4) Notem que, conforme t vai para infinito, xt na curva logística tende assintoticamente para um determinado valor. Para entender mais sobre a construção de uma curva logística e sua aplicação, você pode consultar o site: http://astro.temple.edu/~dhill001/logistic/logistic.html o qual contém um demo explicativo. Exercício -3 Testando tendências nas séries temporais e construindo uma série temporal com uma tendência a) Nos sites de dados do CPC obtenha uma série temporal mensal de dados de temperatura à superfície ou em 1000 hPa em regiões continentais (ex: SE do Brasil, NE dos Estados Unidos, Europa ocidental) e outra em regiões tropicais (continentais ou oceânicas) . Utilizando o Excel, identifique nos dados mensais a existência de tendências. Remova essa tendência e mostre os resultados. Calcule as médias anuais da temperatura e repita o procedimento. Compare o que você obteve para as diferentes séries e discuta seus resultados. b) Usando a função no Excel que gera dados aleatórios (*random), construa uma série temporal que possua uma componente aleatória mais uma tendência. Faça esse procedimento supondo tendência: linear, quadrática, Gompertz e curva logística. Plote os resultados e discuta as diferenças. *Entre no Excel e clique em Tools, Data Analysis, Random (escolha distribuição Gaussiana, com média 0 e desvio-padrão 1). Se não aparecer “Data Analysis” quando clicar em Tools, então clique em Tools, Add Ins, e adicione Analysis ToolPak. Se o seu Excel for em Português, então entre em Ferramentas, Suplemento e adicione as análises que aparecerem nessa opção. 2.3. Existência de sazonalidade nas séries temporais: A grande maioria das séries temporais de variáveis meteorológicas exibe variações com período anual. Por exemplo, a temperatura é usualmente maior no verão que no inverno, e, em algumas localidades, a precipitação possui um ciclo sazonal bem definido. O efeito sazonal pode ser aditivo ou multiplicativo. Um efeito sazonal aditivo é do tipo: Xt = mt + St + εt (2.5) Onde mt é o nível médio dessazonalizado no tempo t, St é o efeito sazonal no tempo t e εt é o erro aleatório. Evidentemente, a análise de séries temporais que exibem uma variação sazonal depende se se deseja medir esse efeito ou eliminá-lo. Por exemplo, suponha que o objetivo seja analisar as anomalias diárias na temperatura em um determinado ponto no globo durante um certo período (por exemplo, junho a agosto). A nossa premissa é que, a cada dia, a temperatura de uma localidade é determinada pelo ciclo solar aliado a variabilidades em outras escalas temporais. Se quisermos entender o papel dessas outras escalas temporais na modulação da temperatura ao longo da estação, devemos primeiro remover o 4 efeito do ciclo sazonal causado por efeito de translação da terra em torno do sol. Assim, removido esse efeito, podemos analisar a temperatura a cada dia, comparando um dia com o outro. Veja que em algumas localidades, como veremos ao longo do curso, o ciclo anual da temperatura pode ser o ciclo de maior importância em termos de espectro, porém é, em geral, o mais bem conhecido e de pouco interesse. Um procedimento muito comum em meteorologia para se remover o ciclo anual é o seguinte (Hartmann and Michelsen, 1989, J. Atmospheric Sciences, 18, 2838-2862): 1) 2) Suponha uma série temporal Xt,y , onde t=[1,365] dias e y=[1, total], total = número de anos que corresponde à série, por exemplo, 20 anos. No caso do exemplo, o fato de considerarmos t variando de 1 a 365 dias, significa que estamos desprezando anos bissextos. Essa suposição pode ser feita usando como aproximação que em anos bissextos o valor da variável no dia 28 de fevereiro é uma média entre 28 e 29 de fevereiro. Outra alternativa é simplesmente eliminar o dia 29 (depende muito do que se está estudando). Uma possibilidade de remoção do ciclo anual é primeiro calcular uma média diária obtida nos 20 anos de dados: total X t = ∑ xt , y , a barra indica a média, xt é a observação no dia t e y representa o ano. y =1 Em geral, quando fazemos esse procedimento, devido ao reduzido número de anos, a série temporal média apresenta pequenas oscilações que são o resultado de variabilidades interanuais as quais não são ‘alisadas’ quando procedemos à operação de média descrita acima. Assim, um procedimento sugerido em Hartmann e Michelsen (1989): 3) O ciclo anual resultante X t deve ser ‘alisado’ usando um filtro com pesos do tipo 1-2-1, passado 300 vezes. O Número de vezes que se passa um filtro pode ser decidido investigando-se o comportamento final da série. O que é e como construir um filtro 1-2-1? (CH) A idéia é usar um FILTRO LINEAR o qual converte uma série temporal {xt} em outra {yt} por uma operação linear: yt = +s ∑a x r =− q r t+r (2.6) onde {a r } é o conjunto de pesos. Para alisar flutuações locais e estimar a média local, devemos escolher pesos tais que ∑a r = 1 . Essa operação é freqüentemente chamada de MÉDIA MÓVEL (MOVING AVERAGE). As médias móveis são freqüentemente simétricas com s=q e aj=aj. O exemplo mais simples de um filtro simétrico é do tipo: a r = 1 /( 2 q + 1) para r=-q,...,+q. O valor alisado de xt é dado por: Sm( xt ) = 1 +q ∑ xt +r (2.7) 2q + 1 r =− q Note que nesse caso, o peso em cada elemento é igual a 1. O filtro conhecido como 1-2-1 considera uma média móvel de três elementos, porém com pesos {a r } na Eq. 2.6 iguais a 0.25, 0.5 e 0.25. Em ambas as bordas, o procedimento é calcular a média entre to e to+1 (borda inferior) e entre tf e tf-1 (borda superior). Este procedimento é bastante útil quando se deseja determinar as anomalias em relação ao ciclo sazonal e deve ser aplicado mesmo se considerarmos pêntadas. Exercício -4 Determinando o ciclo sazonal nas séries temporais e calculando anomalias. Nos sites de dados do CPC obtenha uma série temporal de dados diários de temperatura à superfície ou em 1000 hPa em regiões continentais subtropicais ou extratropicais, e em regiões tropicais sobre o continente ou oceânicas. As séries temporais devem ter pelo menos 10 anos. Desconsidere anos bissextos (ou seja, eliminem o dia 29 de fevereiro, quando existir, calculando a média entre 28 e 29). a) Utilizando o Excel ou programando na linguagem de sua escolha, calcule a média diária dos dados e plote as séries temporais médias. Discuta as eventuais diferenças observadas entre as duas séries. b) Aplique o filtro 1-2-1 conforme descrito acima nas suas séries temporais e plote novamente a série resultante. c) Determine agora uma nova série temporal de anomalias, isto é, o valor observado menos o valor médio. Plote essa série temporal para os 3 primeiros anos. Discuta os resultados. 2.4 DOMÍNIO DE TEMPO E DE FREQUÊNCIA (WI) 5 Existem duas aproximações fundamentais para a análise de séries temporais: análise no DOMÍNIO DO TEMPO e análise no DOMÍNIO DE FREQUÊNCIA. Estas duas aproximações são processadas de forma bem diferente e podem ser vistas como bastante distintas. Contudo, não são independentes! Ao contrário, são métodos complementares que são ligados matematicamente. Os métodos de domínio temporal procuram caracterizar as séries de dados nos mesmos termos em que são observados e reportados. A ferramenta primária para a caracterização de relações entre valores de dados na aproximação do domínio temporal é a FUNÇÃO DE AUTOCORRELAÇÃO. Matematicamente, as análises do domínio temporal operam no mesmo espaço dos valores dos dados. As análises no domínio de freqüência representam as séries de dados em termos de contribuições ocorrendo em diferentes escalas temporais, ou freqüências características. Cada escala temporal é representada por um par de funções seno e co-seno. A série completa é considerada como resultante de efeitos combinados de uma coleção de ondas senoidais e co-senoidais oscilando em diferentes taxas. A soma destas ondas reproduz os dados originais, mas comumente é a intensidade relativa das componentes individuais das ondas que são de interesse primário. Análises no domínio de freqüência ocorrem no espaço matemático definido por esta coleção de senos e co-senos. Isto é, as análises no domínio de freqüência envolvem transformação dos valores de n dados originais em coeficientes que multiplicam um igual número de funções periódicas (os senos e co-senos). Estes métodos são comumente aplicados em séries temporais atmosféricas e são de grande valia para vários propósitos. 2.4.1 Função de auto-correlação (JW) Em estatística clássica as observações xt (t=1,2,...,N) de alguns parâmetros físicos podem ser consideradas independentes desde que os experimentos que geraram essas observações sejam fisicamente independentes. Se a distribuição de probabilidade fx(x) associada com as observações é NORMAL ou GAUSSIANA, a mesma pode ser completamente caracterizada pela sua média: μ=E[X]= e sua variância: ∫ ∞ xf x ( x )dx (2.8) −∞ σ 2 = E [( X − μ ) 2 ] = ∫ ( x − μ ) 2 f x ( x )dx ∞ −∞ (2.9) A média mede a localização ou centro de gravidade da distribuição e a variância a sua variabilidade em torno da média. Se as observações xt formam parte da série temporal, então apenas se o processo que gerou os dados for puramente aleatório os valores vizinhos serão independentes. Em geral, os valores vizinhos de uma série temporal são CORRELACIONADOS. Assim, além de se especificar a média μ e a variância σ2, é necessário no caso de uma série Normal estacionária que se especifique a função de auto-covariância: γ ( k ) = E [( X (t ) − μ )( X (t + k ) − μ )] (2.10) Na prática, a função de auto-covariância pode ser estimada por : c( k ) = N −k ∑ (x 1 N t =1 t − x )( xt +k − x ) (2.11) onde, x= 1 N N ∑x t =1 t (2.12) é a média da série observada. O plot de c(k) versus k (conhecido como ‘lag’ ou intervalo no tempo) é chamado de FUNÇÃO DE AUTOCOVARIÂNCIA AMOSTRAL da série temporal. Algumas vezes, é conveniente quando comparamos séries com diferentes escalas de medida, normalizarmos a Eq. (2.11) dividindo pela variância c(0), de forma a obtermos a FUNÇÃO DE AUTOCORRELAÇÃO AMOSTRAL: r(k ) = c( k ) (2.13) c(0) o que é equivalente a: N −k r (k ) = ∑ (x t =1 t N − x )( xt + k − x ) ∑ ( xt − x ) 2 (2.13b) t =1 O plot de r(k) versus k é também conhecido como “Correlograma”. A função de auto-correlação é útil em algumas situações porque fornece uma visão do jeito como dependência da série cai com o ‘lag’ ou separação k entre pontos da série. Entretanto, a função de auto-correlação é as vezes muito difícil de interpretar como veremos a seguir. OBSERVAÇÕES 1. 2. NOTEM QUE EXISTE POUCO SIGNIFICADO EM SE CALCULAR rk PARA VALORES DE k MAIORES QUE N/4 QUANDO N NÃO É MUITO GRANDE, É PREFERÍVEL CALCULAR A AUTO-COVARIÂNCIA COMO: 6 ck = 1 N −k N −k ∑ (x t =1 t − x )( xt + k − x ) Exemplos de correlogramas podem ser vistos em CH. Algumas recomendações e observações de CH: a) Séries aleatórias: Se uma série é completamente aleatória, então para grande N, r(k) ≅ 0 para todos os valores diferentes de zero de k. b) Correlação de curto-termo. Séries estacionárias freqüentemente exibem correlação de curto-termo caracterizada por um valor de r(1) razoavelmente alto, seguido por uns poucos coeficientes os quais, embora maiores que zero, tendem a ficar sucessivamente menores. Valores de r(k) para “lags” maiores (intervalos de tempo maior) tendem a ser aproximadamente iguais a zero. Séries que produzem esse tipo de correlograma são aquelas que uma observação acima da média tende a ser seguida por uma ou mais observações acima da média, e analogamente para observações abaixo da média. Séries com alternâncias: Se uma série temporal tem tendência a alternar, com sucessivas observações em diferentes lados da média geral, então o correlograma também tende a alternar. O valor de r(1) será negativo. Contudo, o valor de r(2) será positivo uma vez que as observações no lag 2 tenderão a estar do mesmo lado da média. c) d) Séries não-estacionárias: Se a série contém uma tendência, então os valores de r(k) não caem para zero exceto para valores de “lag” (intervalo de tempo) muito altos. Isto ocorre porque uma observação de um lado da média geral tende a ser seguida por um grande número de observações do mesmo lado da média por causa da tendência. Note que pouco pode ser inferido por um correlograma desse tipo porque a tendência domina todas as outras características. Por essa razão, note que a função de auto-correlação só é útil para séries temporais ESTACIONÁRIAS. Por isso as tendências nas séries temporais devem ser removidas antes de proceder à análise de autocorrelação. e) Flutuações sazonais: Se a série temporal contém uma flutuação sazonal, então o correlograma também exibirá uma oscilação na mesma freqüência. Por exemplo, com observações mensais, r(6) será grande e negativo enquanto r(12) será grande e positivo. Em particular, se xt segue um padrão senoidal então r(k) também seguirá o mesmo padrão. Por exemplo, se: xt = a cos tω (2.14) onde a é uma constante e a freqüência ω é tal que 0< ω < π. Pode ser demonstrado que (Exercício 2.3 CH): r (k ) = cos kω para N grande geralmente, um correlograma desse tipo tem pouca utilidade prática. Se a variação sazonal for removida, então o correlograma pode fornecer alguma informação útil. f) Pontos Aberrantes (“outliers”): Se a série contém um ou mais pontos aberrantes, o correlograma pode ser seriamente afetado. Neste caso, é recomendável que os pontos aberrantes sejam ajustados de alguma forma antes de começar uma análise formal. Por exemplo, se existe um ponto aberrante na série temporal e este não é ajustado, então o plot de xt contra xt+k conterá dois pontos extremos os quais irão fazer com que os coeficientes de correlação amostral caiam para zero. Se existirem dois pontos aberrantes este efeito é ainda mais notável, exceto quando o “lag” iguala-se à distância entre os pontos aberrantes. Quando isso acontece, para esse lag pode ocorrer um alto coeficiente de auto-correlação. ________________________________________________________________________ Exercício-5: Determinando a função de auto-correlação para diferentes séries temporais. Aconselha-se que os dados sejam gerados dentro do Excel (por exemplo), mas que o estudante se entusiasme a fazer um programa simples, em qualquer linguagem, para o cálculo da função de auto-correlação. a) b) c) d) Utilizando a função randômica do Excel, gere uma série temporal com 365 pontos aleatórios e faça o correlograma dessa série. Adicione uma tendência linear na série randômica e recalcule o correlograma. Compare com o observado no item anterior e discuta os resultados Adicione uma oscilação do tipo cossenoidal na série do item (a) e recalcule o correlograma. Compare com o observado no item (a) e no item (b) Adicione pontos aberrantes aleatoriamente (alguns valores com +3 desvios-padrão em relação à média) e recalcule o correlograma. Mostre diferenças com respeito ao obtido nos itens anteriores. 2.5 DOMÍNIO DE FREQUÊNCIA: ANÁLISE HARMÔNICA (WI) Análises no domínio de freqüência envolvem a representacao das séries temporais em termos de contribuições feitas em diferentes escalas de tempo. Por exemplo, uma série temporal de dados de temperatura horária de uma localização em latitudes médias irá, em geral, exibir fortes variações tanto na escala de tempo diária (correspondendo ao ciclo anual de aquecimento solar) e na escala de tempo anual (refletindo a marcha das estacdoes). No domínio do tempo, estes ciclos apareceriam como valores positivos altos na função de auto-correlação para “lags” em (ou aproximadamente em) 24 h para o ciclo diurno e 24x365 = 8760h para o ciclo anual. Quando pensamos nas séries no domínio de freqüência, estamos falando de contribuições para o total da variabilidade das séries temporais em períodos de 24 e 8760h, ou em freqüências de 1/24 = 0.0417 h-1 e 1/8760 = 0.000114 h-1. A análise harmônica consiste da representação de flutuações ou variações em uma série temporal que se originou da adição de uma série de funções seno e co-seno. Estas funções trigonométricas são “harmônicos” quee são escolhidos como tendo freqüências que são múltiplas da freqüência “fundamental” determinada pelo tamanho amostral da série de dados. 7 2.5.1 Funções seno e co-seno (WI) Vamos aqui rever brevemente a natureza das funções cos(α) e sin(α). O argumento de ambas é a quantidade α, medida em unidades de ângulo. Estas unidades podem ser tanto em graus ou radianos. A figura 2.1 mostra porções de um cos (linha sólida) e seno (linha tracejada), no intervalo angular que varia de 0 a 5π/2 (0- 450º). As funções seno e co-seno estendem-se até ângulos indefinidamente grandes positivos ou negativos. O mesmo padrão de onda se repete a cada 2π radianos ou 360º, tal que: cos(2πk + α)=cos(α) (2.15) Onde k é um inteiro qualquer. Uma equação análoga vale para a função seno. Isto é, ambos seno e co-seno são periódicos. Ambas funções oscilam em torno de seu valor médio igual a zero e tem valor máximo de +1 e mínimo de -1. A função co-seno é maximizada a 0o, 360º, e assim por diante, enquanto a função seno é maximizada a 90º, 450º, e assim por diante. Estas duas funções têm exatamente a mesma forma, mas estão deslocadas uma em relação a outra por 90º. Em outras palavras, se você deslocar a função co-seno para a direita em 90º, produz a função seno e deslocando o seno para esquerda em 90º produz a função co-seno: π cos(α − ) = sin(α ) (2.16) e 2 sin(α + π 2 ) = cos(α ) (2.17) 1.5 1 0.5 0 -0.5 0 45 90 135 180 225 270 315 360 405 450 -1 -1.5 cos sin Figure 1. Trechos da função seno (linha tracejada) e co-seno (linha sólida) para intervalols entre 0o e 450o, ou 0 a 5π/2 radianos. Cada curva executa um ciclo completo a cada 2π radianos e se estende de -∞ a +∞. 2.5.2. Representação de uma série temporal simples com uma Função Harmônica. (WI) Mesmo na situação mais simples de uma série temporal que possui um caráter senoidal e executa um único ciclo sobre o curso de n observações, exsistem três pequenas dificuldades a serem superadas a fim de se usar um seno ou co-seno para representá-la: 1) 2) 3) O argumento de uma função trigonométrica é um ângulo, enquanto os dados da série são função do tempo. As funções co-seno e seno flutuam entre +1 e -1, enquanto os dados geralmente flutuam entre diferentes limites. A função co-seno tem máximo valor para α=0 e α=2π. Ambos seno e co-seno podem assim estar posicionados arbitrariamente na horizontal com respeito aos dados. A solução para o primeiro problema aparece quando consideramos o comprimento dos dados (n) como constituindo um ciclo completo, ou período fundamental. Uma vez que o período fundamental corresponde a 360º ou 2π radianos em medida angular, é fácil re-escalar proporcionalmente o tempo à medida angular usando: ⎛ 360o ⎞ ⎛ ⎞ ⎟ ⎟⎟⎜⎜ α = ⎜⎜ = 360o ⎟ ciclo ⎠⎜ n unidades de tempo / ciclo ⎟ n ⎝ ⎝ ⎠ t unidade de tempo t (2.18) ou 8 ⎞ ⎛ ⎟ ⎜ t unidade de tempo t π 2 ⎞⎜ ⎛ ⎟ = 2π α =⎜ (2.19) ⎟ n ⎟⎟ ⎝ ciclo ⎠⎜⎜ ⎝ n unidades de tempo / ciclo ⎠ Em outras palavras, se considerarmos um ciclo completo o número total de pontos n de uma série temporal, então conforme andamos em t, as equações 2.18 e 2.19 indicam em que proporção t encontra-se desse ciclo completo. Assim, a quantidade: ω1 = 2π (2.20) n é chamada de FREQUÊNCIA FUNDAMENTAL. Esta quantidade é uma freqüência angular e tem dimensões físicas de radianos por unidade de tempo. A freqüência fundamental corresponde ao período igual ao comprimento dos dados. O subscrito ‘1’ indica que ω1 pertence a onda que executa um CICLO COMPLETO sobre a série de dados inteira. O problema (2) acima é tratado fazendo-se um deslocamento da função seno ou co-seno para cima ou para baixo do nível geral dos dados, e então “esticando” ou “comprimindo” verticalmente até que seu intervalo corresponda ao dos dados. Como isso pode ser feito? Uma vez que a média de uma onda pura co-seno ou seno é zero, simplesmente adicionar o valor médio da série de dados ao co-seno assegura que o mesmo irá flutuar em torno do valor médio. O “esticamento” ou a “compressão” pode ser obtido pela multiplicação de uma função seno ou co-seno por uma constante C1 que é conhecida como AMPLITUDE. Novamente, note que o subscrito 1 significa que trata-se da amplitude do HARMÔNICO FUNDAMENTAL. Uma vez que as funções seno e co-seno possuem máximo e mínimo entre ± 1, é bastante óbvio mostrar que o máximo/mínimo dessas funções quando multiplicadas pela amplitude C1 é igual a ± C1. Assim, se combinarmos as soluções dos problemas (1) e (2) temos para a série temporal: yt = y + C1 cos( 2πt ) (2.21) n Normalmente, há situações em que é necessário transladar uma função harmônica lateralmente a fim de ter um casamento de cristas e cavados na série de dados. Este tipo de desenvolvimento é mais conveniente quando a função co-seno é utilizada, porque esta função possui máximo valor quando o ângulo na qual ela opera é zero. Mudar a função co-seno para a direita por um ângulo φ1 resulta em uma nova função que é maximizada em ω1t=2πt/n=φ1, yt = y + C1 cos( 2πt − φ1 ) (2.22) n O ângulo φ1 é chamado de ângulo de fase ou fase. Mudar a função co-seno para a direita desta quantidade requer subtrair a fase φ1, de modo que o argumento da função co-seno é zero quando (2πt/n)=φ1, o que equivale a dizer que a yt em 2.22 é maximizada para t=φ1n/2π. Para compreender a aplicação das eqs. 2.20 a 2.22, faça o exercício proposto 6. ___________________________________________________________________________ Exercício-6. Considere os dados a seguir, correspondentes à temperatura média obtida em em Ithaca (NY): Table 1- Temperatura média mensal t=mês 1 2 3 4 5 6 7 8 9 10 11 12 a) b) c) Temperatura oC -5.44 -5.17 0.11 6.89 12.67 17.94 20.44 19.5 15.67 9.72 4.06 -2.56 Plote a curva em questão no Excel. Determine ω1 e calcule a função co-seno para todos os meses da tabela plotando no mesmo gráfico. Determine a média e a função yt como em 2.21, isto é, re-escalone a função co-seno considerando a amplitude C1 como metade do intervalo entre o máximo e o mínimo da série. Plote os resultados no mesmo gráfico. Determine a função yt como em 2.22, isto é, ajuste a fase φ1. Discuta como isso foi feito. 9 Desafio-6a. Utilize uma série temporal de dados em que você esteja trabalhando (mensal, diária ou pentadal), ou use uma das séries temporais que foram fornecidas no curso. Calcule a média mensal dos dados, obtendo uma tabela semelhante à Tabela-1. Faça os itens (a) a (c), conforme indicado acima. Note que, não necessariamente, você irá observar uma série com um ciclo anual semelhante ao que foi observado no exercício. Você também pode efetuar seu trabalho utilizando dados que não são anuais e interpretar o que significa o primeiro harmônico nesse caso. ___________________________________________________________________________ 2.5.3. Estimativa de Amplitude e fase de um único harmônico (WI) No exercício anterior, estimou-se C1 e φ1 de uma forma mais grosseira. Existem formas mais adequadas de se fazer o mesmo procedimento. O método mais simples é usar a identidade trigonométrica: cos(α − φ1 ) = cos(φ1 ) cos(α ) + sin(φ1 ) sin(α ) (2.23) Substituindo α=2πt/n (Eq. 2.19) e multiplicando-se ambos os lados pelas amplitudes C1, obtém-se: C1 cos( 2πt 2πt 2πt 2πt 2πt − φ1 ) = C1 cos(φ1 ) cos( ) + C1 sin(φ1 ) sin( ) = A1 cos( ) + B1 sin( ) (2.24) n n n n n onde, A1 = C1 cos(φ1 ) (2.25) e B1 = C1 sin(φ1 ) (2.26) A Equação 2.24 indica, portanto, que um harmônico que possui amplitude C1 e fase φ1 é equivalente a soma de um seno com um coseno, ambos sem deslocamento de fase. Note que o seno aparece justamente como uma expansão de um co-seno de α menos a fase. São justamente B1 e A1 os coeficientes que irão modular a amplitude do harmônico com o tempo. Note que ambos coeficientes B1 e A1 são determinados pelo seno e co-seno da fase, respectivamente. Fazendo a transformação de variáveis da Eq. 2.24 tal que x1=cos(2πt/n) e x2=sin(2πt/n), temos uma equação resultante que é uma equação de regressão com dois preditores. Desta forma, dada a série de dados yt, podese aplicar a essa transformação qualquer software que ache uma estimativa de regressão com o método de mínimos quadrados para encontrar os coeficientes A1 e B1, tendo yt como preditando. O mesmo pacote de regressão irá produzir a média dos valores do preditando como o intercepto bo. Pode-se então encontrar C1 resolvendo o sistema das equações 2.25 e 2.26, ou seja: [ A12 + B12 = C12 (sin 2 (φ1 ) + cos 2 (φ1 )) ⇒ C1 = A12 + B12 ] 1/ 2 (2.27) Para se encontrar a função de fase φ1, são resolvidas as seguintes equações: ⎧ −1 B1 ⎪ tan A 1 ⎪ ⎪ −1 B1 ± π , ou ± 180o φ1 = ⎨tan A 1 ⎪ ⎪π o ⎪ , or 90 2 ⎩ A1 > 0 A1 < 0 (2.28) A1 = 0 Nota: uma vez que as funções trigonométricas são periódicas, o mesmo ângulo de fase é produzido se adicionarmos ou subtrairmos 180º (meio círculo) se A1<0. Decidimos se somamos ou subtraímos 180º dependendo da relação 0 < φ ≤ 2π Para encontrarmos os parâmetros A1 e B1 podemos usar o método de mínimos quadrados. Para casos especiais em que os valores dos dados estão igualmente espaçados no tempo com nenhum dado faltante (missing values) as propriedades do seno e co-seno permitem que os parâmetros de mínimos quadrados sejam obtidos mais facilmente usando: A1 = 2 n ⎛ 2πt ⎞ yt cos⎜ ⎟ ∑ n t =1 ⎝ n ⎠ B1 = 2 ⎛ 2πt ⎞ yt sin⎜ ⎟ ∑ n t =1 ⎝ n ⎠ (2.29) e n (2.30) ___________________________________________________________________________ Exercício-7: Para os dados do Exercício 6, obtenha os coeficientes do harmônico anual da temperatura em Ithaca: A1, B1, C1 e a fase φ1. Use o Excel para fazer seus cálculos e compare com a fase (em graus) e o coeficiente C1 obtidos pelo método do “olhômetro” do exercício anterior. Comente o que essas diferenças podem significar na interpretação dos resultados. Desafio-7b. Se você utilizou uma serie temporal distinta da fornecida, encontre os coeficientes para esta série. ___________________________________________________________________________ 2.5.4 Harmônicos de ordem mais elevada: (WI) Nos cálculos efetuados nos exercícios 7 e 8 produziu-se uma única função dcosseno passando próxima aos 12 meses de temperaturas médias. Este ajuste bom ocorreu porque a forma do ciclo anual de temperatura na localidade em questão é aproximadamente senoidal, com um ciclo completo executado nos 12 pontos da série temporal. Em geral, não se espera que um único harmônico irá representar todas as demais séries temporais de temperatura em diferentes localidades. Entretanto, se adicionarmos mais preditores a uma regressão múltipla poderá melhor o ajuste de um conjunto de dados, adicionando mais ondas do tipo co-seno a uma análise harmônica irá melhorar o ajuste a uma série temporal. 10 Assim, dada uma série temporal consistindo de n pontos, a mesma pode ser exatamente reproduzida. Isto significa que é possível encontrar uma função harmônica que passa através de cada um dos pontos e isso é feito adicionando uma série de n/2 funcoes harmônicas: n/2 n/2 ⎧ ⎧ ⎡ 2πkt ⎤⎫ ⎡ 2πkt ⎤ ⎡ 2πkt ⎤ ⎫ − φk ⎥ ⎬ = y + ∑ ⎨ Ak cos⎢ + Bk sin ⎢ yt = y + ∑ ⎨Ck cos⎢ ⎬ (2.31) ⎥ ⎣ n ⎦⎭ ⎣ n ⎦ ⎣ n ⎥⎦ ⎭ k =1 ⎩ k =1 ⎩ Note que o índice k na equação acima indica que a Eq. 2.24 vale para qualquer co-seno, independentemente de sua freqüência. A onda de coseno obtida para k=1 na Eq. 2.31 é, portanto, o primeiro harmônico determinado anteriormente. Os outros n/2 - 1 termos da somatória da Eq. 2.31 são harmônicos de ordem mais alta, ou ondas co-seno com freqüências: ωk = 2πk n (2.32) que são múltiplos inteiros da freqüência fundamental ω1. Note ainda que cada onda representada pelo co-seno possui sua própria fase φ e sua própria amplitude C. Portanto, k dentro da equação é de suma importância. Primeiro harmônico: k=1, α=2πkt/n, fase φ1 e amplitude C1 → um ciclo completo (0 a 2π rad) para t=0 a t=n Segundo harmônico: k=2, α=2πkt/n, fase φ2 e amplitude C2 → um ciclo completo (0 a 2π rad) para t=0 a t=n/2 → um ciclo completo para t=n/2 e t=n Terceiro harmônico: k=3, α=2πkt/n, fase φ3 e amplitude C3 → três ciclos completos entre t=0 e t=n Os coeficientes Ak e Bk podem ser encontrados, NO CASO MAIS GERAL, (Eq. 2.31) usando-se método de regressão linear múltipla, após as transfomações de dados x1=cos(2πt/n), x2=sin(2πt/n), x3= cos(2π2t/n), x4= sin(2π2t/n), x5= cos(2π3t/n), x6= sin(2π2t/n) e assim por diante. PARA SÉRIES TEMPORAIS IGUALMENTE ESPAÇADAS NO TEMPO (SEM VALORES FALTANTES (MISSING VALUES): Podemos generalizar as equações 2.29 e 2.30 para: Ak = 2 n ⎛ 2πkt ⎞ yt cos⎜ ⎟ (2.33) ∑ n ⎠ n ⎝ t =1 e 2 n ⎛ 2πkt ⎞ Bk = ∑ yt sin⎜ ⎟ (2.34) n n ⎠ ⎝ t =1 Como determinar os coeficientes Ak e Bk? Fica evidente pelas Eqs. 2.33 e 2.34 que a determinação desses coeficientes é viável usando-se um programa de computador. Nesse caso, o algoritmo utilizaria um contador para o k (isto é, fixa-se k) e faz-se a somatória do cosseno (ou seno de α) de t=1 a t=n (ou seja, t é um outro contador que deverá fazer o ‘looping’ ou somatória que será interno ao looping do contador k). Esse raciocínio pode ser também aplicado para uso dentro do Excel, por exemplo. No entanto, para séries muito grandes, usa-se um método mais eficiente que será comentado adiante. Uma vez calculados estes coeficientes, temos então que obter a fase e a amplitude do harmônico. a) Amplitude: Ck =[ Ak2 + Bk2 ]1 / 2 (2.35a) O algoritmo para encontrar Ck é, portanto, muito simples, uma vez achados Ak e Bk e deve ser feito após o looping da somatória. Quanto à fase, deve-se primeiro testar o coeficiente Ak e resolver: ⎧ B ⎪tan −1 k Ak ⎪ ⎪ B ⎪ φk = ⎨tan −1 k ± π , ou ± 180 o Ak ⎪ ⎪π o ⎪ , ou 90 2 ⎪⎩ A1 > 0 A1 < 0 (2.35b) A1 = 0 ♦ Quantos harmônicos podem ser determinados? Notem que uma função de regressão múltipla irá passar por todos os pontos de dados e irá exibir R2=100% se o número de preditores é o mesmo que número de dados. A série de cossenos da Eq. 2.31 é um exemplo. Quer dizer, os n/2 harmônicos que aparecem na somatória da Eq. 2.31 implica na existência de n variáveis preditores. Uma vez que a média amostral da Eq. 2.31 é efetivamente um dos parâmetros estimados, correspondentes ao intercepto bo, um ajuste para a Eq. 2.31 é necessário quando n é ÍMPAR. Neste 11 caso, uma soma sobre apenas (n-1)/2 harmônicos é requerida para representar completamente a função. Isto é, (n-1)/2 harmônicos é requerido para completamente representar a função. Isto é, para n ÍMPAR: (n-1)/2 amplitudes + (n-1)/2 ângulos de fase + a média amostral = n Para n PAR: A somatória da Eq. 2.31 continua sendo até n/2. Entretanto, temos (n)/2 amplitudes + [(n)/2 – 1] ângulos de fase + a média amostral = n Neste caso, o ângulo de fase para o harmônico final e mais alto φn/2 = 0 ♦ Quantos harmônicos podemos usar? Essa pergunta depende das finalidades. Podemos usar todos os n/2 harmônicos se queremos uma função que represente exatamente através de todos os pontos da série. Contudo, se o objetivo é represnetar o ciclo anual de uma quantidade climatológica, os primeiros poucos harmônicos podem fornecer uma representação adequada. É isso que se faz na prática. ___________________________________________________________________________ Exercício-8a: Para os dados do Exercício 7, obtenha os coeficientes do segundo e terceiro harmônicos da temperatura em Ithaca. Use o Excel para fazer seus cálculos compare com o exercício anterior e comente os resultados. Exercício 8b: (Tarefa de casa): Utilize dados paleoclimáticos da radiação solar (encontram-se no link 'dados') e calcule a amplitude do harmônico correspondente a 11 e 22 anos. Discuta seus resultados ___________________________________________________________________________ 12