Análise de Séries Temporais - Trabalho Final

Modelagem da Temperatura Máxima de Porto Alegre (RS)

Introdução

A análise das séries temporais de temperatura máxima diária é fundamental na meteorologia e em várias áreas decisivas, da agricultura ao planejamento urbano. O Instituto Nacional de Meteorologia (INMET) disponibiliza uma vasta quantidade de informações coletadas automaticamente por estações meteorológicas em todo o país, proporcionando uma valiosa fonte de dados para estudos e análises.

Neste trabalho, focamos na análise da estação meteorológica de Porto Alegre (RS), abrangendo os útlimos 11 anos, mas excluindo os dados de 2023 para validação do modelo, ou seja, os modelos foram feitos utilizando os dados de temperatura máxima de Porto Alegre entre 2012 e 2022. Essa análise é de grande importância para compreender e prever variações climáticas, sendo relevante na agricultura, gestão de recursos hídricos, segurança pública, turismo, etc. Usaremos técnicas de análise de séries temporais e modelos de previsão para ajustar dados até o final de 2022, e avaliaremos sua precisão prevendo os dados até 31/07/2023 e comparando com os dados já observados. Dessa forma, buscamos aplicar as metodologias estudadas na disciplina de séries temporais contribuindo para o entendimento das tendências climáticas locais e aprimorar a capacidade de prever as variações nas temperaturas máximas diárias, com potenciais aplicações em diversas áreas.

Análise Descritiva

No gráfico abaixo, temos a série temporal da temperatura máxima diária (°C) na cidade de Porto Alegre (RS) de \(01/01/2012\) até \(30/07/2023\). Em uma primeira análise, a série aparenta ter um comportamento estacionário e uma sazonalidade anual com 11 “ciclos” bem definidos ao longo dos 11 anos analisados. Para reforçar essas constatações, iremos analisar a função de autocorrelação (ACF), a função de autocorrelação parcial (ACF) e a série decomposta.

Pelo gráfico da função de autocorrelação gerado, podemos visualizar os padrões sazonais na série, com picos significativos na ACF em intervalos regulares. Além disso, temos fortes indícios da presença de componentes autoregressivos (AR).

Já com o gráfico da função de autocorrelação parcial (PACF), podemos visualizar que não há lags significativos na PACF após os 100 primeiros lags aproximadamente, indicando que a série não possui dependência de longo prazo. No entanto, a presença de lags significativos no início da série e em direções diferentes nos dão indícios da presença de componente de média móvel (MA) e a necessidade de modelos AR de diferentes ordens em partes diferentes da série.

Abaixo, o gráfico mostra os seguintes componentes:

• Série Temporal Original: A série temporal bruta, que inclui tendência, sazonalidade e erro.

• Tendência: A componente de tendência, que representa a direção geral do comportamento da série ao longo do tempo.

• Sazonalidade: A componente sazonal, que representa os padrões que se repetem em intervalos regulares, como sazonalidade anual, mensal, etc.

• Erro (Resíduos): A componente de erro, que representa o que resta na série após a remoção da tendência e da sazonalidade.

Metodologia

Teste de Raíz Unitária

Antes de partirmos para a modelagem da série, iremos realizar o teste Teste de Dickey-Fuller Aumentado (ADF), utilizado para avaliar se a série temporal é realmente estacionária ou não estacionária.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  treino
## Dickey-Fuller = -5.7667, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary

Pelo resultado acima podemos ver que o \(p.valor < 0.05\), ou seja, a um nível de significância de 5%, rejeitamos a hipótese nula. Portando, concluímos que há evidências de que a série seja estacionária, estando também de acordo com os resultados das análises gráficas realizadas anteriormente.

Modelo 1 - ARIMA(1,0,3)

Para modelar séries temporais com componentes sazonais, utilizamos modelos ARIMA sazonais, conhecidos como SARIMA (Seasonal AutoRegressive Integrated Moving Average). Um modelo SARIMA é uma extensão do ARIMA que incorpora termos sazonais para capturar a variação sazonal nos dados.

O SARIMA é denotado como \(SARIMA(p, d, q)(P, D, Q)\), onde:

• “p”: é a ordem do componente AR não sazonal.
• “d”: é a ordem de integração não sazonal.
• “q”: é a ordem do componente MA não sazonal.
• “P”: é a ordem do componente AR sazonal.
• “D”: é a ordem de integração sazonal.
• “Q”: é a ordem do componente MA sazonal.

No software R, utilizamos a função \(auto.arima()\) que automatiza o processo de seleção do modelo ARIMA. Ela explora várias combinações de parâmetros p (ordem do componente AR), d (ordem de integração) e q (ordem do componente MA) e fornece o modelo ARIMA com o menor valor de AIC ou BIC como o melhor modelo. Após identificar possíveis modelos ARIMA não sazonais, utilizando o parâmetro \(seazonal=TRUE\), o algoritmo então estende a busca para modelos SARIMA, que incluem componentes sazonais (P, D, Q, s) além dos componentes não sazonais. Ele ajusta modelos SARIMA com diferentes valores de P (AR sazonal), D (integração sazonal), Q (MA sazonal) e s (período da sazonalidade) e calcula os critérios de informação para cada modelo.

Para ajustar o modelo, utilizamos apenas os dados até 2022 (no qual chamamos de “treino”) para então fazer a previsão até os dias atuais (a qual chamamos de “teste”).

model1 <- auto.arima(treino, seasonal = TRUE)
summary(model1)

## Series: treino 
## ARIMA(1,0,3) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1      ma2     ma3     mean
##       0.9899  -0.2715  -0.3643  -0.175  25.7374
## s.e.  0.0027   0.0155   0.0164   0.015   0.9371
## 
## sigma^2 = 10.6:  log likelihood = -10437.78
## AIC=20887.56   AICc=20887.58   BIC=20925.34
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.002761661 3.254056 2.566752 -1.981664 10.87387 0.5666404
##                    ACF1
## Training set 0.01244031

O modelo pode ser definido como,

\[X_t = c + \phi X_{t-1} + \theta \epsilon_{t-1} + \theta \epsilon_{t-2} + \theta \epsilon_{t-3} + \epsilon\]

em que os parâmetros encontrados foram,

\[X_t = 25.74 + 0.99 X_{t-1} -0.27 \epsilon_{t-1} -0.36 \epsilon_{t-2} -0.17 \epsilon_{t-3} + \epsilon\]

Fazendo o forecast (previsão) do modelo para os dias atuais temos:

No entanto, para validarmos o modelo, precisamos antes verificar se os resíduos não são autocorrelacionados a partir da análise da função de autocorrelação dos resíduos e do teste Ljung-Box:

## 
##  Box-Ljung test
## 
## data:  forecasted$residuals
## X-squared = 34.597, df = 20, p-value = 0.02236

O correlograma mostra que as autocorrelações para os erros de previsão praticamente não excedem os limites de significância para defasagens 1-20. Além disso, o p-valor para o teste Ljung-Box é \(0.02236\) indicando que não há muita evidência de autocorrelações diferentes de zero nas defasagens 1-20. Portanto, podemos concluir que os resíduos não são autocorrelacionados.

Ademais, outro passo para verificar se o modelo preditivo pode ser melhorado é investigar se os erros de previsão são normalmente distribuídos com média zero e variância constante.

O teste Ljung-Box não revela evidências suficientes de que exista autocorrelações diferentes de zero nos erros de previsão na amostra, as variações nos erros de previsão parecem aproximadamente constantes e os erros de previsão parecem distribuídos normalmente. Portanto, podemos concluir que o modelo ARIMA(1,0,3) é um modelo válido.

Modelo 2 - SARIMA(3,0,0)(0,1,0)[365]

Apesar do modelo anterior ter sido considerado válido e com uma acurácia significativa ao ajuste dos dados de treino, podemos perceber que mesmo com o parâmetro \(seazonal=TRUE\) a função auto.arima() não indicou componentes sazonais que havíamos verificado que estão presentes na série nas etapas anteriores. Portanto, neste modelo iremos utilizar o parâmetro \(D=1\), indicando que se deseja aplicar uma diferenciação sazonal de primeira ordem à série temporal. Isso significa que, em cada período sazonal, os valores da série serão subtraídos dos valores do período sazonal anterior para remover a sazonalidade.

model2 <- auto.arima(treino, D = 1)
summary(model2)

## Series: treino 
## ARIMA(3,0,0)(0,1,0)[365] 
## 
## Coefficients:
##          ar1      ar2     ar3
##       0.7051  -0.1854  0.0507
## s.e.  0.0165   0.0200  0.0165
## 
## sigma^2 = 19.94:  log likelihood = -10642.58
## AIC=21293.16   AICc=21293.18   BIC=21317.98
## 
## Training set error measures:
##                       ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.03650731 4.255992 3.212511 -1.85509 13.57526 0.7091991
##                      ACF1
## Training set -0.001597896

O modelo pode ser definido como,

\[\nabla_{365}X_t = \phi X_{t-1} + \phi X_{t-2} + \phi X_{t-3} + \epsilon\]

em que os parâmetros encontrados foram

\[\nabla_{365}X_t = 0.71 X_{t-1} -0.19 X_{t-2} + 0.05 X_{t-3} + \epsilon\]

Fazendo o forecast (previsão) do modelo para os dias atuais temos:

Assim como feito anteriormento, iremos verificar se o novo modelo é válido,

## 
##  Box-Ljung test
## 
## data:  forecasted2$residuals
## X-squared = 33.538, df = 20, p-value = 0.02943

O correlograma mostra que as autocorrelações para os erros de previsão praticamente não excedem os limites de significância para defasagens 1-20. Além disso, o p-valor para o teste Ljung-Box é \(0.02943\) indicando que não há muita evidência de autocorrelações diferentes de zero nas defasagens 1-20. Portanto, podemos concluir que os resíduos não são autocorrelacionados.

Além disso, os erros de previsão têm variância constante ao longo do tempo e são normalmente distribuídos com média zero

Portanto, podemos concluir que o modelo SARIMA(3,0,0)(0,1,0)[365] também é válido.

Conclusão

	Modelo.1	Modelo.2
AIC	20887.56	21293.16

Portanto, encontramos que os dois modelos gerados são válidos. Pelo AIC, podemos ver que o menor deles é o Modelo 1, indicando que este teria um melhor ajuste aos dados. No entanto, o gráfico evidencia que o Modelo 2 teve um desempenho muito mais acurado ao compararmos as médias dos valores preditos.