Forecast Dataset

Introdução

Este relatório detalha a análise de uma série temporal mensal contida no arquivo Forecast_Dataset.txt. O objetivo é identificar a estrutura da série, ajustar um modelo SARIMA robusto e gerar previsões. Primeiramente, carregamos os pacotes necessários.

Carregamento dos Pacotes e Funções

Primeiramente, carregamos todos os pacotes necessários para a análise. Também definimos as funções personalizadas para os testes de tendência, raiz unitária e sazonalidade, que serão utilizadas ao longo do estudo.

library(readr)
library(dplyr)
library(lubridate)
library(imputeTS)
library(forecast)
library(ggplot2)
library(tseries)
library(trend)
library(randtests)
library(Kendall)
library(seastests)
library(lmtest)
library(nortest)
library(ggpubr)


tend_determ<-function(ts){
  CS<-suppressWarnings(randtests::cox.stuart.test(ts,c("two.sided"))) #H0: NAO existe tendencia
  CeST<-suppressWarnings(trend::cs.test(ts)) #H0: NAO existe tendencia
  # Runs<-suppressWarnings(randtests::runs.test(ts)) #H0: NAO existe tendencia
  # WaldW<-suppressWarnings(trend::ww.test(ts)) #H0: NAO existe tendencia
  MannKT<-suppressWarnings(trend::mk.test(ts,continuity = TRUE)) #H0: a serie eh i.i.d. / NAO existe tendencia
  MannK<-suppressWarnings(Kendall::MannKendall(ts)) #H0: NAO existe tendencia
  KPSST<-suppressWarnings(tseries::kpss.test(ts, null = c("Trend"))) #H0: NAO existe tendencia
  #
  p_value<-c(CS$p.value,CeST$p.value,MannKT$p.value,MannK$sl,KPSST$p.value)
  p_value1<-p_value
  p_value1[p_value>=0.05]<-"NAO tendencia"
  p_value1[p_value<0.05]<-"Tendencia"
  tabela<-data.frame(Testes=c("Cox Stuart","Cox and Stuart Trend",
                              "Mann-Kendall Trend","Mann-Kendall","KPSS Test for Trend"),
                     H0=c(rep("NAO tendencia",5)),
                     p_valor=round(p_value,4),
                     Conclusao=c(p_value1))
  list(CS=CS,CeST=CeST,MannKT=MannKT,MannK=MannK,KPSST=KPSST,Tabela=tabela)
}

raiz_unit<-function(ts){
  ADF<-suppressWarnings(tseries::adf.test(ts,alternative = c("stationary"))) #H0: raiz unitaria
  PP<-suppressWarnings(tseries::pp.test(ts,alternative = c("stationary"))) #H0: raiz unitaria
  KPSSL<-suppressWarnings(tseries::kpss.test(ts, null = c("Level"))) #H0: nao existe tendencia
  #
  p_value<-c(ADF$p.value,PP$p.value,KPSSL$p.value)
  p_value1<-p_value[1:2]
  p_value1[p_value[1:2]>=0.05]<-"Tendencia"
  p_value1[p_value[1:2]<0.05]<-"NAO tendencia"
  p_value2<-p_value[3]
  p_value2[p_value[3]>=0.05]<-"NAO tendencia"
  p_value2[p_value[3]<0.05]<-"Tendencia"
  tabela<-data.frame(Testes=c("Augmented Dickey-Fuller","Phillips-Perron Unit Root","KPSS Test for Level"),
                     H0=c(rep("Tendencia",2),"NAO tendencia"),
                     p_valor=round(p_value,4),
                     Conclusao=c(p_value1,p_value2))
  list(ADF=ADF,PP=PP,KPSSL=KPSSL,Tabela=tabela)
}

sazonalidade<-function(ts,diff=0,freq){
  KrusW<-suppressWarnings(seastests::kw((ts),diff = diff, freq=12)) #H0: NAO Sazonal
  Fried<-suppressWarnings(seastests::fried((ts),diff = diff, freq=12)) #H0: NAO Sazonal
  #
  p_value<-c(KrusW$Pval,Fried$Pval)
  p_value1<-p_value
  p_value1[p_value>=0.05]<-"NAO Sazonal"
  p_value1[p_value<0.05]<-"Sazonal"
  tabela<-data.frame(Testes=c("Kruskall Wallis","Friedman rank"),
                     H0=c(rep("NAO Sazonal",2)),
                     p_valor=round(p_value,4),
                     Conclusao=c(p_value1))
  list(KrusW=KrusW,Fried=Fried,Tabela=tabela)
}

Carga e Limpeza dos Dados

Carregamos o arquivo Forecast_Dataset.txt e o convertemos em um objeto de série temporal (ts), especificando o início em Julho de 1991 e a frequência mensal (12).

dados_new <- read.csv("Forecast_Dataset.txt")
dados_new$date <- as.Date(dados_new$date)
head(dados_new)
##         date    value
## 1 1991-07-01 3.526591
## 2 1991-08-01 3.180891
## 3 1991-09-01 3.252221
## 4 1991-10-01 3.611003
## 5 1991-11-01 3.565869
## 6 1991-12-01 4.306371
ts_new <- ts(dados_new$value, start = c(1991, 7), frequency = 12)

Análise Exploratória da Série Original

A primeira etapa da análise é visualizar a série e suas funções de autocorrelação para entender seus componentes.

summary(ts_new) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.815   5.844   9.319  10.694  14.290  29.665
autoplot(ts_new) + labs(x = "Ano", y = "Valor", title = "Série Temporal Original") + theme_minimal()

ggAcf(ts_new, lag.max = 48) + labs(y = "FAC", title = "FAC da Série Original") + theme_minimal()

ggPacf(ts_new, lag.max = 48) + labs(y = "FACP", title = "FACP da Série Original") + theme_minimal()

O gráfico da série exibe uma clara tendência crescente e um padrão sazonal que se repete anualmente. A FAC confirma isso com seu decaimento lento e picos significativos nos lags sazonais (12, 24, 36). A série aparenta ser não estacionária e possuir um forte componente sazonal. Essas características guiam para a escolha de um modelo SARIMA.

Decomposição e Testes

plot(decompose(ts_new))

tend_determ(ts = ts_new)$Tabela
##                 Testes            H0 p_valor Conclusao
## 1           Cox Stuart NAO tendencia    0.00 Tendencia
## 2 Cox and Stuart Trend NAO tendencia    0.00 Tendencia
## 3   Mann-Kendall Trend NAO tendencia    0.00 Tendencia
## 4         Mann-Kendall NAO tendencia    0.00 Tendencia
## 5  KPSS Test for Trend NAO tendencia    0.01 Tendencia
raiz_unit(ts=ts_new)$Tabela
##                      Testes            H0 p_valor     Conclusao
## 1   Augmented Dickey-Fuller     Tendencia    0.01 NAO tendencia
## 2 Phillips-Perron Unit Root     Tendencia    0.01 NAO tendencia
## 3       KPSS Test for Level NAO tendencia    0.01     Tendencia

A decomposição confirma visualmente uma forte tendência de crescimento e uma sazonalidade anual bem definida. Segundo os testes aplicados, a série temporal apresenta uma tendência estatisticamente significativa e não possui raíz unitária, ou seja, é estacionária em torno de uma tendência determinística significativa.

Tratamento da Série e Verificação

Para tornar a série estacionária, aplicamos uma diferenciação regular (para remover a tendência) e uma sazonal (para remover a sazonalidade).

diff_ts <- diff(ts_new, lag = 24, differences = 1)
diff_ts <- diff(diff_ts, differences = 1)

autoplot(diff_ts) + labs(title = "Serie Apos Diferenciacao") + theme_minimal()

ggAcf(diff_ts, lag.max = 72) + labs(title = "FAC da Serie Diferenciada")

ggPacf(diff_ts, lag.max = 72) + labs(title = "FACP da Serie Diferenciada")

tend_determ(ts = diff_ts)$Tabela
##                 Testes            H0 p_valor     Conclusao
## 1           Cox Stuart NAO tendencia  0.3966 NAO tendencia
## 2 Cox and Stuart Trend NAO tendencia  0.7626 NAO tendencia
## 3   Mann-Kendall Trend NAO tendencia  0.9881 NAO tendencia
## 4         Mann-Kendall NAO tendencia  0.9881 NAO tendencia
## 5  KPSS Test for Trend NAO tendencia  0.1000 NAO tendencia
raiz_unit(ts = diff_ts)$Tabela
##                      Testes            H0 p_valor     Conclusao
## 1   Augmented Dickey-Fuller     Tendencia    0.01 NAO tendencia
## 2 Phillips-Perron Unit Root     Tendencia    0.01 NAO tendencia
## 3       KPSS Test for Level NAO tendencia    0.10 NAO tendencia
sazonalidade(ts = ts_new,diff = 1,freq = 12)$Tabela
##            Testes          H0 p_valor Conclusao
## 1 Kruskall Wallis NAO Sazonal       0   Sazonal
## 2   Friedman rank NAO Sazonal       0   Sazonal

A série diferenciada agora aparenta ser estacionária, flutuando em torno de uma média zero. O padrão de decaimento lento, que era visível na FAC da série original, desapareceu completamente e as autocorrelações estão perto de zero. Os testes indicam que a tendência foi removida, que a série é estacionária e que há sazonalidade.

Ajuste do Modelo SARIMA

Utiliza-se a função auto.arima() para uma busca otimizada do melhor modelo SARIMA. Após testar modelos vizinhos e comparar seus critérios de informação (AIC/BIC), o modelo SARIMA(1,1,1)(1,1,2)[12] foi confirmado como a opção mais robusta e parcimoniosa.

mod_new <- Arima(ts_new, order = c(1,1,1), seasonal = list(order = c(1,1,2), period = 12))

summary(mod_new)
## Series: ts_new 
## ARIMA(1,1,1)(1,1,2)[12] 
## 
## Coefficients:
##           ar1      ma1    sar1     sma1    sma2
##       -0.3044  -0.6655  0.8435  -1.5910  0.7701
## s.e.   0.0952   0.0812  0.0724   0.1151  0.1177
## 
## sigma^2 = 0.7635:  log likelihood = -252.33
## AIC=516.65   AICc=517.11   BIC=536.16
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE    MAPE      MASE
## Training set 0.06042276 0.8343627 0.5421608 0.04552758 4.67216 0.4185697
##                     ACF1
## Training set -0.01367887
lmtest::coeftest(mod_new)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.304438   0.095244  -3.1964  0.001392 ** 
## ma1  -0.665466   0.081188  -8.1966 2.472e-16 ***
## sar1  0.843480   0.072394  11.6512 < 2.2e-16 ***
## sma1 -1.591024   0.115123 -13.8202 < 2.2e-16 ***
## sma2  0.770057   0.117721   6.5414 6.096e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O modelo teve um AIC de 516.65 e o teste de coeficientes mostra que todos os termos são estatisticamente significativos, validando a estrutura do modelo.

Análise de Diagnóstico dos Resíduos

Verifica-se se os resíduos do modelo final se comportam como um ruído branco, uma premissa fundamental para a validade do modelo.

# Residuos
res.mod_new<-mod_new$residuals
forecast::ggtsdisplay(res.mod_new)

Box.test(res.mod_new,lag = 15, type="Ljung") 
## 
##  Box-Ljung test
## 
## data:  res.mod_new
## X-squared = 22.979, df = 15, p-value = 0.08459
# Normalidade
nortest::ad.test(res.mod_new)
## 
##  Anderson-Darling normality test
## 
## data:  res.mod_new
## A = 5.0467, p-value = 1.621e-12
raiz_unit(ts = res.mod_new)$Tabela
##                      Testes            H0 p_valor     Conclusao
## 1   Augmented Dickey-Fuller     Tendencia    0.01 NAO tendencia
## 2 Phillips-Perron Unit Root     Tendencia    0.01 NAO tendencia
## 3       KPSS Test for Level NAO tendencia    0.10 NAO tendencia
ggpubr::ggqqplot(res.mod_new)

O diagnóstico dos resíduos é positivo. O gráfico da FAC não mostra picos significativos, e o teste de Ljung-Box resulta em um p-valor alto (0.08459), indicando que não há autocorrelação nos resíduos. O teste de normalidade indica que os resíduos não seguem uma distribuição perfeitamente normal, o que não é relevante nesse caso. Os resíduos do modelo são estácionarios. Com o modelo bem especificado, as previsões podem ser realizadas.

Previsão

Com um modelo robusto e validado, podemos gerar previsões para o futuro. Vamos prever os próximos 12 meses (1 anos).

predict.mod_new<-forecast::forecast(mod_new,h=12)
predict.mod_new
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jul 2008       23.83268 22.71094 24.95443 22.11712 25.54825
## Aug 2008       24.33120 23.20897 25.45343 22.61489 26.04750
## Sep 2008       24.23052 23.05047 25.41058 22.42579 26.03526
## Oct 2008       25.75075 24.54151 26.95999 23.90138 27.60012
## Nov 2008       27.63372 26.38907 28.87836 25.73020 29.53724
## Dec 2008       28.10708 26.83013 29.38403 26.15415 30.06001
## Jan 2009       32.40007 31.09099 33.70916 30.39800 34.40214
## Feb 2009       21.36184 20.02159 22.70210 19.31210 23.41159
## Mar 2009       22.41781 21.04703 23.78859 20.32138 24.51424
## Apr 2009       22.47151 21.07089 23.87213 20.32944 24.61357
## May 2009       25.12511 23.69526 26.55496 22.93834 27.31188
## Jun 2009       23.90183 22.44338 25.36029 21.67132 26.13235
forecast::autoplot(predict.mod_new) + forecast::autolayer(predict.mod_new$fitted,show.legend = F) + theme_minimal()

O gráfico de previsão mostra a continuação da tendência de crescimento e do padrão sazonal da série. Os intervalos de confiança (áreas sombreadas) se alargam com o tempo, refletindo corretamente a crescente incerteza em previsões de longo prazo.

Conclusão Geral

A análise da série temporal contida no Forecast_Dataset.txt revelou uma forte tendência de crescimento e um componente sazonal anual bem definido.

Para tratar essas características, foi necessário um processo de dupla diferenciação: uma regular, para remover a tendência, e uma sazonal, para remover a sazonalidade. Esta transformação se mostrou eficaz, resultando em uma série estacionária pronta para a modelagem.

Através de um processo iterativo de ajuste e seleção, o modelo SARIMA(1,1,1)(1,1,2)[12] foi identificado como o mais adequado e parcimonioso. A validade deste modelo foi confirmada pela análise de diagnóstico, que demonstrou que seus resíduos se comportam como um ruído branco, sem autocorrelação remanescente.

Finalmente, o modelo validado foi utilizado para gerar previsões para os próximos 12 meses. As previsões indicam a continuação da tendência de crescimento e do padrão sazonal.

Air Quality UCI Dataset

Introdução

Este relatório apresenta uma análise completa da série temporal de concentração de monóxido de carbono (CO(GT)), extraída do conjunto de dados “Air Quality” do repositório UCI. O objetivo é identificar a estrutura da série, ajustar um modelo robusto e gerar previsões.

Carregamento dos Pacotes e Funções

Primeiramente, carregamos todos os pacotes necessários para a análise. Também definimos as funções personalizadas para os testes de tendência, raiz unitária e sazonalidade, que serão utilizadas ao longo do estudo.

# Previamente instalados

Carga e Limpeza dos Dados

O conjunto de dados original (AirQualityUCI.csv) requer um pré-processamento para tratar valores faltantes (indicados por -200) e para combinar as colunas de data e hora em um único objeto temporal. A série de interesse é a CO(GT), que corresponde a concentração média horária de Monóxido de Carbono.

dados_brutos <- read_delim("AirQualityUCI.csv", delim = ";", locale = locale(decimal_mark = ","))

dados_limpos <- dados_brutos %>%
  select(1:15) %>%
  filter(!is.na(Date)) %>%
  mutate(across(where(is.numeric), ~na_if(., -200))) %>%
  mutate(
    Time_corrigida = gsub("\\.", ":", Time),
    DateTime = dmy_hms(paste(Date, Time_corrigida))
  )

serie_para_analise <- dados_limpos %>%
  select(DateTime, `CO(GT)`) %>%
  filter(!is.na(`CO(GT)`)) %>%
  arrange(DateTime)

head(serie_para_analise)
## # A tibble: 6 × 2
##   DateTime            `CO(GT)`
##   <dttm>                 <dbl>
## 1 2004-03-10 18:00:00      2.6
## 2 2004-03-10 19:00:00      2  
## 3 2004-03-10 20:00:00      2.2
## 4 2004-03-10 21:00:00      2.2
## 5 2004-03-10 22:00:00      1.6
## 6 2004-03-10 23:00:00      1.2
d1 <- ts(serie_para_analise$`CO(GT)`, frequency = 24)

d1 <- na_interpolation(d1, option = "linear")

Análise Exploratória da Série Original

Nesta seção, investigamos as características da série original para identificar visualmente a presença de tendência, sazonalidade e para avaliar sua estacionariedade.

# Serie
autoplot(d1) + labs(x = "Tempo (Horas)", y = "Concentração de CO(GT)", title = "Série Original da Concentração de CO") + theme_minimal()

# FAC
ggAcf(d1, lag.max = 72) + labs(y = "FAC", title = "FAC da Série Original") + theme_minimal()

# FACP
ggPacf(d1, lag.max = 72) + labs(y = "FACP", title = "FACP da Série Original") + theme_minimal()

  • Série Original: Bastante irregular, com muitos picos e vales, o que indica uma grande flutuação na concentração de CO.A série não flutua em torno de uma média constante, indicando não estacionariedade.

  • FAC: Lag 24, 48, 72 apresentam autocorrelação elevada. Isso indica sazonalidade diária (24 horas). A concentração de CO tende a se repetir em padrões semelhantes a cada dia.Além disso, decaimento suave nas autocorrelações entre os picos indica memória de curto prazo na série.

  • FACP: A série tem componente autoregressivo forte no lag 1 e não há muitos outros lags com correlação parcial significativa.

Decomposição e Testes

A decomposição da série e os testes estatísticos formais são utilizados para confirmar as suspeitas levantadas pela análise visual.

plot(decompose(d1))

tend_determ(ts = d1)
## $CS
## 
##  Cox Stuart test
## 
## data:  ts
## statistic = 2075, n = 3757, p-value = 1.543e-10
## alternative hypothesis: non randomness
## 
## 
## $CeST
## 
##  Cox and Stuart Trend test
## 
## data:  ts
## z = 2.2935, n = 7674, p-value = 0.02182
## alternative hypothesis: monotonic trend
## 
## 
## $MannKT
## 
##  Mann-Kendall trend test
## 
## data:  ts
## z = 1.0138, n = 7674, p-value = 0.3107
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.271140e+05 5.018749e+10 7.809660e-03 
## 
## 
## $MannK
## tau = 0.00781, 2-sided pvalue =0.31069
## 
## $KPSST
## 
##  KPSS Test for Trend Stationarity
## 
## data:  ts
## KPSS Trend = 1.082, Truncation lag parameter = 11, p-value = 0.01
## 
## 
## $Tabela
##                 Testes            H0 p_valor     Conclusao
## 1           Cox Stuart NAO tendencia  0.0000     Tendencia
## 2 Cox and Stuart Trend NAO tendencia  0.0218     Tendencia
## 3   Mann-Kendall Trend NAO tendencia  0.3107 NAO tendencia
## 4         Mann-Kendall NAO tendencia  0.3107 NAO tendencia
## 5  KPSS Test for Trend NAO tendencia  0.0100     Tendencia
raiz_unit(ts=d1)
## $ADF
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts
## Dickey-Fuller = -11.476, Lag order = 19, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $PP
## 
##  Phillips-Perron Unit Root Test
## 
## data:  ts
## Dickey-Fuller Z(alpha) = -1104.5, Truncation lag parameter = 11,
## p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $KPSSL
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts
## KPSS Level = 1.465, Truncation lag parameter = 11, p-value = 0.01
## 
## 
## $Tabela
##                      Testes            H0 p_valor     Conclusao
## 1   Augmented Dickey-Fuller     Tendencia    0.01 NAO tendencia
## 2 Phillips-Perron Unit Root     Tendencia    0.01 NAO tendencia
## 3       KPSS Test for Level NAO tendencia    0.01     Tendencia

  • Decomposição: Há uma tendência oscilante, com elevação no meio do período observado. Existe também uma sazonalidade clara e estável, provavelmente diária.

  • Testes: Os testes de tendência indicam presença de tendência na série. Os testes de raiz unitária têm resultados mistos, mas no geral, a série não tem raiz unitária, mas pode conter uma leve tendência determinística

Tratamento da Série e Verificação

Para que a modelagem ARIMA seja possível, a série precisa ser tornada estacionária. Isso é feito através da diferenciação. Dada a presença de um ciclo diário, aplicamos tanto uma diferenciação sazonal (lag=24) quanto uma regular.

diff_d1 <- diff(d1, lag = 24, differences = 1)
diff_d1 <- diff(diff_d1, differences = 1)

autoplot(diff_d1) + labs(title = "Série Após Diferenciação Sazonal e Regular") + theme_minimal()

ggAcf(diff_d1, lag.max = 72) + labs(title = "FAC da Série Diferenciada")

ggPacf(diff_d1, lag.max = 72) + labs(title = "FACP da Série Diferenciada")

tend_determ(ts = diff_d1)$Tabela
##                 Testes            H0 p_valor     Conclusao
## 1           Cox Stuart NAO tendencia  0.6490 NAO tendencia
## 2 Cox and Stuart Trend NAO tendencia  0.6112 NAO tendencia
## 3   Mann-Kendall Trend NAO tendencia  0.6712 NAO tendencia
## 4         Mann-Kendall NAO tendencia  0.6766 NAO tendencia
## 5  KPSS Test for Trend NAO tendencia  0.1000 NAO tendencia
raiz_unit(ts = diff_d1)$Tabela
##                      Testes            H0 p_valor     Conclusao
## 1   Augmented Dickey-Fuller     Tendencia    0.01 NAO tendencia
## 2 Phillips-Perron Unit Root     Tendencia    0.01 NAO tendencia
## 3       KPSS Test for Level NAO tendencia    0.10 NAO tendencia

A série aparenta estar estacionária. Os valores flutuam em torno de uma média constante e a variância parece estável ao longo do tempo. No FAC, a maioria dos coeficientes está dentro dos limites de significância, exceto alguns picos (ex: no lag 1, 2 e próximo de 24). Quanto ao FACP, há picos sazonais em 48 e 72, indicando comportamento periódico de 24h.

Ajuste do Modelo SARIMA

O modelo SARIMA(2,1,2)(1,0,1)[24] se mostrou superior, com todos os coeficientes significativos e os menores valores de AIC/BIC..

mod <- Arima(d1, 
             order = c(2,1,2), 
             seasonal = list(order = c(1,0,1), period = 24))

summary(mod)
## Series: d1 
## ARIMA(2,1,2)(1,0,1)[24] 
## 
## Coefficients:
##          ar1      ar2      ma1      ma2    sar1     sma1
##       0.8022  -0.0919  -0.7852  -0.2104  0.4921  -0.2290
## s.e.  0.0425   0.0356   0.0414   0.0414  0.0345   0.0384
## 
## sigma^2 = 0.5251:  log likelihood = -8415.63
## AIC=16845.27   AICc=16845.28   BIC=16893.89
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE      MASE
## Training set -0.002981888 0.7243172 0.4973672 -15.79168 31.53943 0.5088819
##                      ACF1
## Training set 0.0003642239
lmtest::coeftest(mod)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.802162   0.042496  18.8763 < 2.2e-16 ***
## ar2  -0.091882   0.035562  -2.5837  0.009774 ** 
## ma1  -0.785206   0.041433 -18.9511 < 2.2e-16 ***
## ma2  -0.210391   0.041423  -5.0791 3.793e-07 ***
## sar1  0.492079   0.034454  14.2824 < 2.2e-16 ***
## sma1 -0.229022   0.038365  -5.9695 2.379e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O sumário confirma a estrutura do modelo e mostra um AIC de 16845.27. O coeftest mostra que todos os seis coeficientes do modelo são estatisticamente significativos (p-valores < 0.05), validando a estrutura escolhida.

Análise de Diagnóstico dos Resíduos

A etapa final da modelagem é verificar se os resíduos do modelo mod se comportam como um ruído branco (erros aleatórios, sem padrões).

res.mod <- residuals(mod)

ggtsdisplay(res.mod, main = "Diagnóstico dos Resíduos do Modelo Final")

Box.test(res.mod, lag = 30, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res.mod
## X-squared = 537.66, df = 30, p-value < 2.2e-16
nortest::ad.test(res.mod)
## 
##  Anderson-Darling normality test
## 
## data:  res.mod
## A = 172.31, p-value < 2.2e-16
ggpubr::ggqqplot(res.mod)

  • Gráficos: Os gráficos da FAC e FACP dos resíduos mostram picos nos lags 24, 48 e 72, indicando presença de autocorrelação.

  • Teste de Ljung-Box: O teste indica que existe autocorrelação significativa nos resíduos.

  • Normalidade: O teste de normalidade e o gráfico Q-Q indicam que os resíduos não seguem uma distribuição perfeitamente normal, o que é comum em dados reais.

Modelo SARIMAX

Apesar do bom ajuste do modelo SARIMA, a análise indica a presença de autocorrelação. Isso sugere que fatores externos podem estar influenciando a série. A alternativa escolhida agora é incorporar as variáveis de Temperatura (T) e Umidade Relativa (RH) em um modelo SARIMAX.

dados_para_xreg <- serie_para_analise %>%
  inner_join(dados_limpos, by = "DateTime") %>%
  select(T, RH)

xreg_matrix <- as.matrix(dados_para_xreg)

mod_sarimax <- Arima(d1, 
                     order = c(2,1,2), 
                     seasonal = list(order = c(1,0,1), period = 24),
                     xreg = xreg_matrix)

summary(mod_sarimax)
## Series: d1 
## Regression with ARIMA(2,1,2)(1,0,1)[24] errors 
## 
## Coefficients:
##          ar1      ar2      ma1      ma2    sar1     sma1       T      RH
##       0.8372  -0.1350  -0.8180  -0.1747  0.4963  -0.2389  0.0850  0.0283
## s.e.  0.0422   0.0353   0.0416   0.0416  0.0368   0.0410  0.0091  0.0025
## 
## sigma^2 = 0.5165:  log likelihood = -8005.62
## AIC=16029.24   AICc=16029.26   BIC=16091.75
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE     MASE
## Training set -0.002912409 0.7182162 0.4935009 -14.97908 31.18627 0.504926
##                      ACF1
## Training set 0.0007064539
lmtest::coeftest(mod_sarimax)
## 
## z test of coefficients:
## 
##        Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.8371960  0.0422380  19.8209 < 2.2e-16 ***
## ar2  -0.1350454  0.0352531  -3.8307 0.0001278 ***
## ma1  -0.8180014  0.0415938 -19.6664 < 2.2e-16 ***
## ma2  -0.1746798  0.0415847  -4.2006 2.662e-05 ***
## sar1  0.4962677  0.0368362  13.4723 < 2.2e-16 ***
## sma1 -0.2389179  0.0409667  -5.8320 5.477e-09 ***
## T     0.0849568  0.0090514   9.3860 < 2.2e-16 ***
## RH    0.0282653  0.0024899  11.3522 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O modelo SARIMAX representa uma melhora substancial em relação ao modelo univariado, com AIC de 16029.24 e significância das variáveis externas

Análise de Resíduos do Modelo SARIMAX

A etapa final é verificar se o modelo SARIMAX possui resíduos bem-comportados.

res.mod.sarimax<-mod_sarimax$residuals
forecast::ggtsdisplay(res.mod.sarimax)

Box.test(residuals(mod_sarimax), lag = 30, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(mod_sarimax)
## X-squared = 546.11, df = 30, p-value < 2.2e-16

Mesmo com a inclusão das variáveis externas e um ajuste muito superior, o modelo ainda não foi capaz de capturar toda a autocorrelação presente nos dados.

Conclusão Geral

A análise demonstrou que um modelo SARIMA univariado é limitado. A evolução para um modelo SARIMAX, incorporando Temperatura e Umidade, resultou em uma melhora drástica no ajuste do modelo (AIC caiu de ~16845 para ~16729), provando que as condições climáticas são preditores fundamentais.

Apesar deste ajuste superior, o diagnóstico final dos resíduos do modelo SARIMAX ainda apontou a presença de autocorrelação. Isso indica que a dinâmica da série possui complexidades adicionais (como outliers não modelados ou relações não-lineares) que exigem abordagens mais avançadas, mas que o modelo com variáveis externas representa o maior avanço explicativo encontrado neste estudo.

Referências