Análise e Modelagem do Consumo de Energia Elétrica no Brasil

1 Introdução

A energia elétrica é um recurso fundamental para o progresso do Brasil. Desde sua chegada no final do século XIX, com a construção da primeira usina hidrelétrica em Niterói, o setor elétrico brasileiro se desenvolveu rapidamente. Hoje, o Brasil é o maior produtor de energia elétrica renovável na América Latina.

Sua importância para a economia do país é indiscutível, sendo utilizada em diversos setores como indústria, comércio, agricultura e serviços, gerando emprego e renda.

Atualmente, o Brasil é o 7º maior consumidor de energia do mundo e o maior da América do Sul

O trabalho desenvolvido nesse documento tem como objetivo, analisar e modelar a série temporal de consumo de energia elétrica no Brasil, entre o período de Janeiro de 2004 a Setembro de 2023

Os dados utilizados foram obtidos da organização EPE - Empresa de Pesquisa Energética

Os dados ja limpos e organizados em um csv podem ser acessados via github: Consumo_Brasil.csv

2 Análise Exploratória

consumo_ts_br = readr::read_csv('consumo_eletrico_brasil.csv') |>
  select(-"...1") |>
  mutate(ano_mes = ano_mes_dia |> yearmonth(), .keep = 'unused') |>
  relocate(ano_mes, .before = consumo) |>
  tsibble::as_tsibble(index = ano_mes)

Os sgeuintes dados foram utilizados

consumo_ts_br |>
  rmarkdown::paged_table()

Onde eles apresentaram o seguinte comportamento

consumo_ts_br |>
  autoplot() +
  theme_minimal() +
  labs(title = 'Consumo de Energia Elétrica no Brasil', x = 'Data', y = 'Consumo em MWh')

Vemos que se trata de uma série com uma forte tendência de crescimento. Além disso, é possível observar a presença de ciclos, tais ciclos foram detalhadamente estudados ao longo do trabalho, onde um padrão temporal foi observado, se tratando portanto de uma sazonalidade

Também vemos grandes quedas no uso de energia nos anos de 2009 e 2020

A queda em 2009 pode ser explicada pela forte crise econômica mundial de 2008, como noticiado em G1 - Queda do uso de energia 2009: “Como ressaltado ao longo do ano, o mercado brasileiro de energia elétrica sofreu forte impacto da crise financeira internacional, porém seus efeitos se concentraram na classe industrial, como consequência da imediata e profunda retração da atividade deste segmento”

Já a queda em 2020 está fortemente ligada a pandemia de Covid-19, onde os setores industriais apresentaram forte queda nos primeiros meses do ‘lockdown’, como noticiado em Energe - Queda do uso de energia 2020 “Ao analisar o desempenho do consumo de energia por ramo de atividade, é possível notar que a indústria automotiva e o segmento têxtil ainda lideram as maiores quedas no mercado livre. O setor de veículos apresentou, no início da pandemia (mês de abril) recuo de 65,6% no consumo de energia. Já a área têxtil registrou retração de 46,3% também no mesmo período.”

Buscando melhor entender os dados trabalhados, principalmente sua sazonalidade, o método STL de decomposição de séries temporais foi utilizado

O Seasonal and Trend decomposition using Loess(STL) busca decompor a série em 3 componentes: tendência, sazonalidade e parte aleatória

consumo_ts_br |>
  model(
    STL(consumo ~ trend() +
          season(window = "periodic"),
        robust = TRUE)) |>
  components() |>
  autoplot() +
  theme_minimal() +
  labs(title = 'Decomposição via STL',
       subtitle = 'Consumo = Tendência + Sazonalidade + Ruído')

A partir da decomposição da série vemos a forte presença de sazonalidade com ciclos anuais, indicando que os modelos com parcelas sazonais anuais serão aqueles com melhor ajuste

Além disso, a parcela aleatória discriminou fortemente a queda no uso de energia de 2020, fato esse que não foi visto para o ano de 2009. Isso indica que os dados do ano de 2020 podem ser mostrar como valores atípicos no etapa de modelagem

3 Modelagem

Os seguintes passos foram seguidos:

• Análise da FAC e FACP
  • Utilizar os gráficos de Autocorrelação e Autocorrelação Parcial para investigar o melhor modelo para os dados, assim como os parâmetros desse modelo
• Análise de Variância da Série
  • Visualizar se a série possui variância heterocedástica, onde se necessário, uma transformação deve ser aplicada
• Estimação dos Modelos Não Sazonais e Análise de Resíduos
  • Procura dos melhores parâmetros dos modelos que não possuem uma parcela sazonal
• Estimação dos Modelos Sazonais e Análise de Resíduos
  • Procura dos melhores parâmetros dos modelos que possuem uma parcela sazonal
• Análise do Modelo Final
  • Interpretação, desempempnho e previsão com o melhor modelo estimado

3.1 FAC e FACP

Na etapa de análise exploratória vimos que a série trabalhada possui uma forte tendência de crescimento, além de apresentar uma sazonalidade anual

Assim para a aplicação de modelos que suponhem estacionariedade, esperamos que seja necessário pelo menos 1 diferenciação

A série apresentou o seguinte comportamento da FAC e FACP:

consumo_ts_br |>
  gg_tsdisplay((consumo), plot_type = "partial", lag_max = 50)

Obervando a FAC, vemos um decaimento extremamente lento dos valores de autocorrelações, indicando que a série trabalhada é não estacionária

O teste de Dickler Fuller Aumentado com as seguintes hipóteses foi aplicado:

\[H_0: \text{Série é não estacionária}\]

\[H_1: \text{Série é estacionária}\]

consumo_ts_br$consumo |>
  tseries::adf.test(k=12)

    Augmented Dickey-Fuller Test

data:  consumo_ts_br$consumo
Dickey-Fuller = -1.7759, Lag order = 12, p-value = 0.6701
alternative hypothesis: stationary

Assim como era esperado, o teste não rejeitou a hipótese de não estacionariedade da série e portanto uma diferenciação foi aplicada nos dados

A série diferenciada apresentou o seguinte comportamento

consumo_ts_br |>
  gg_tsdisplay(difference(consumo), plot_type = "partial", lag_max = 50)

A partir de uma diferenciação vemos que a tendência da série foi estabilizada

Aplicando novamente o teste de Dickler Fuller Aumentado, o seguinte resultado foi obtido

consumo_ts_br$consumo |>
  diff() |>
  tseries::adf.test(k=12)

    Augmented Dickey-Fuller Test

data:  diff(consumo_ts_br$consumo)
Dickey-Fuller = -5.0955, Lag order = 12, p-value = 0.01
alternative hypothesis: stationary

Onde \(H_0\) foi rejeitada e portanto é razoável assumir a sua estacionariedade

3.2 Análise da Variância dos Dados

A série diferenciada apresentou o seguinte comportamento

consumo_ts_br |>
  autoplot(consumo |>
             difference()) +
  theme_minimal() +
  labs(title = 'Consumo de Energia Elétrica no Brasil', x = 'Data', y = 'Consumo em MWh')

Observando a série, vemos que há um aumento da variância dos dados ao passar dos anos, sendo uma forte indicação de uma série heterocedástica

Para ajustar tal problema, as seguintes ferramentas podem ser utilizadas:

• Aplicação de uma transformação
  • Uma transformação Log ou Raiz Quadrada
• Modelagem da Variância
  • Utilização de modelos de média e variância, sendo o ARCH(Autoregressive Conditional Heteroskedasticity) o mais conhecido.

A partir do gráfico acima, vemos que o aumento da variância segue um padrão, onde esse aumento ocorre ao longo dos anos, e portanto uma simples transformação Log ou Raiz Quadrada deve reverter tal quadro.

Modelos ARCH são robustos a comportamentos irregulares da variância ao passar do tempo e portanto possuem uma maior complexidade de implementação e interpretação.

Assim, a aplicação de uma transformação Log nos dados foi utilizada para a reverter o quadro de heterocedasticidade

consumo_ts_br |>
  autoplot(consumo |>
             log() |>
             difference()) +
  theme_minimal() +
  labs(title = 'Consumo de Energia Elétrica no Brasil', x = 'Data', y = 'Log do Consumo em MWh')

A FAC e FACP tiveram o seguinte desempenho após a aplicação da transformação

consumo_ts_br |>
  gg_tsdisplay(consumo |>
                 log() |>
                 difference(), plot_type = "partial", lag_max = 50)

Após aplicar a tranformação Log e diferenciar o dados, os comportamentos da FAC e FACP retornaram uma forte autocorrelação com o lag 12 meses, ou seja, a série possui uma sazonalidade anual. Tal fato deve ser levando em conta na etapa de modelagem.

3.3 Modelos Não Sazonais

O seguinte tópico teve como objetivo estimar modelos ARIMA, onde ao final dessa etapa, o melhor modelo não sazonal foi estimado

O modelo deve ter o parametro de diferenciação igual a 1, dado que os dados trabalhados são não estacionários, onde aplicar 1 diferenciação foi suficiente para alterar tal cenário.

Para a seleção do modelo, o parâmetro de diferenciação foi fixado igual a 1, e definiu-se um vetor de 1 a 3 para a procura dos melhores parâmetros autorregressivos e de médias móveis

O método utilizado foi de forward stepwise, sendo o AICc o critério de informação utilizado. Tal critério de informação foi utilizado pois apresenta um menor viés que o AIC (BREWER, Mark J. 2016)(SEN, Liew Khim. 2002)

consumo_ts_br  = consumo_ts_br |>
  mutate(log_consumo = log(consumo))


models_arima_br = consumo_ts_br |> 
  model(arima_test = ARIMA(log_consumo ~ pdq(0:3, 1, 0:3) + PDQ(0,0,0), 
                           trace = F,
                           approximation = F,
                           ic = 'aicc'))

models_arima_br |> report()

Series: log_consumo 
Model: ARIMA(2,1,2) w/ drift 

Coefficients:
         ar1      ar2      ma1     ma2  constant
      1.5742  -0.7916  -1.7028  0.7997     5e-04
s.e.  0.0705   0.0656   0.0690  0.0722     1e-04

sigma^2 estimated as 0.0004637:  log likelihood=573.28
AIC=-1134.56   AICc=-1134.19   BIC=-1113.78

Vemos que o modelo ARIMA(2,1,2) foi o modelo a apresentar o melhor AICc

O modelo estimado apresentou o seguinte comportamento dos resíduos

models_arima_br |>
  gg_tsresiduals()

Os resíduos de um modelo ARMA devem respeitar as seguintes suposições:

• Serem não correlacionados

• Devem possuir média 0 e variância constante

A partir do gráfico de FAC, podemos ver um pico de autocorrelação no lag 12, indicando a presença de forte sazonalidade na série assim como explicado na análise exploratória.

Alguns outros lags após o lag 12 também se mostraram significativamente diferente de 0

Os testes de Box-Pierce e de Ljung-Box para testar a independência entre os resíduos foram utilizados. Os testes possuem as seguintes hipóteses:

\[H_0: \text{Observações são independentes}\] \[H_1: \text{Observações não são independentes}\]

models_arima_br |>
  augment() |>
  features(.resid, box_pierce, lag = 12)

# A tibble: 1 x 3
  .model     bp_stat   bp_pvalue
  <chr>        <dbl>       <dbl>
1 arima_test    56.1 0.000000116

models_arima_br |>
  augment() |>
  features(.innov, ljung_box, lag = 12)

# A tibble: 1 x 3
  .model     lb_stat    lb_pvalue
  <chr>        <dbl>        <dbl>
1 arima_test    59.2 0.0000000318

Vemos que ambos os testes rejeitaram fortemente a hipótese de independência dos resíduos

Assim, como esperado, o modelo sem parcela sazonal não foi capaz de modelar com sucesso o consumo de energia no pais, onde seus resíduos não respeitaram as pressuposições do modelo teórico.

O próximo tópico buscou estimar modelos que possuem ajustes sazonais

3.4 Modelos Sazonais

O seguinte tópico teve como objetivo estimar modelos com parcela sazonal, onde ao final dessa etapa, o melhor modelo sazonal foi estimado

Ao longo da análise exploratória e da modelagem via modelos não sazonais, vimos que o lag 12 possui uma alta autocorrelação, indicando que tal ponto deve ser modelado a partir de uma parcela sazonal

2 modelos foram utilizados:

• Modelo de Sazonalidade Determinística
  • O modelo calcula uma média para o período sazonal fixado e retira esse valor dos dados. Tal modelo é de simples estimação e interpretação, porém não se mostra robusto para modelar a magnitudade da sazonalidade de cada período
• ARIMA Sazonal
  • O modelo tem a capacidade de estimar termos autorregressivos e de médias móveis para a parcela sazonal. Tal modelo possui uma estimação e interpretação mais complexa, porém se mostra robusto para modelar a magnitudade da sazonalidade de cada período

3.4.1 Modelo de Sazonalidade Determinística

O modelo de sazonalidade determinística utiliza o seguinte algoritmo:

1. Dado um série \(Z_t\) estacionária com presença de sazonalidade, identifique o período de sazonalidade

2. Calcule um valor médio sazonal chamado de \(\mu_k\), dado por:

\[\hat\mu_k =\frac{Z_k + Z_{k+s} + Z_{k+2s} + ... + Z_{N-s+k}}{N/s}\]

3. Faça \(N_t = Z_t - \hat\mu_T\), onde assim \(N_t\) será um série estacionária sem sazonalidade

O seguinte bloco de código construi a série estacionária e sem sazonalidade \(N_t\)

medias = consumo_ts_br |>
  mutate(mes = ano_mes |> 
           lubridate::month()) |>
  as_tibble() |>
  select(-ano_mes) |>
  mutate(
    log_consumo_diff = 
           c(0, diff(log_consumo))) |>
  group_by(mes) |>
  summarise(media_mes = log_consumo_diff |>
              mean())




consumo_ts_br_deter_seas = consumo_ts_br |>
  mutate(mes = ano_mes |> lubridate::month()) |>
  left_join(medias, by = 'mes') |>
  mutate(log_consumo_diff_std = (log_consumo |>
           difference()) - media_mes
           )

A série \(N_t\) tem a seguinte estrutura

consumo_ts_br_deter_seas |>
  gg_tsdisplay(log_consumo_diff_std, plot_type = 'partial', lag_max = 50)

Novamente o método de forward stepwise foi utilizado para estimação do modelo de menor AICc. Para isso, definiu-se o vetor de 0 a 5 para a procura dos melhores valores de p e q do modelo \(ARMA(p,q)\) para \(N_t\)

arima_deterministic_season = consumo_ts_br_deter_seas |>
  model(teste = ARIMA(log_consumo_diff_std ~ pdq(0:5,0,0:5) + PDQ(0,0,0))) 

arima_deterministic_season |>
  report()

Series: log_consumo_diff_std 
Model: ARIMA(5,0,0) 

Coefficients:
          ar1      ar2     ar3      ar4      ar5
      -0.1692  -0.1098  0.0422  -0.1808  -0.1464
s.e.   0.0643   0.0642  0.0651   0.0646   0.0647

sigma^2 estimated as 0.0002926:  log likelihood=627.12
AIC=-1242.24   AICc=-1241.88   BIC=-1221.43

Vemos que o modelo com menor AICc estimado foi um \(ARMA(5,0)\)

O modelo apresentou os seguintes resíduos

arima_deterministic_season |>
  gg_tsresiduals()

Via FAC, vemos que os resíduos possuem uma correlação significativa no lag 7 e 12

Portanto, a suposição de independência dos resíduos foi quebrada

Utilizando os novamente os testes de Box-Pierce e de Ljung-Box de indepência, temos que:

arima_deterministic_season |>
  augment() |>
  features(.resid, box_pierce, lag = 12)

# A tibble: 1 x 3
  .model bp_stat bp_pvalue
  <chr>    <dbl>     <dbl>
1 teste     27.4   0.00674

arima_deterministic_season |>
  augment() |>
  features(.innov, ljung_box, lag = 12)

# A tibble: 1 x 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 teste     28.8   0.00429

Ambos os teste apresentaram um p-valor abaixo de 5%.

Assim, dado a análise dos resíduos e aplicação dos testes, é razoável concluir que os resíduos são correlacionados

E portanto vemos que mesmo retirando a média sazonal de 12 meses dados, um modelo ARMA não foi capaz de modelar a série

3.4.2 ARIMA Sazonal

Um modelo ARIMA Sazonal, também chamado de SARIMA, é uma extensão do modelo ARIMA, podendo possuir termos autorregressivos e de médias móveis para uma parcela sazonal

Como escrito anteriormente, a parcela sazonal a ser modelada deve ser o lag 12 (1 ano)

Para a estimação do modelo, o método stepwise baseado no critério AICc foi novamente utilizado, onde os parâmetros de diferenciação foram fixadas em 1 a variável resposta modelada foi o log de consumo

models_sarima_br = consumo_ts_br |> 
  model(arima_test = ARIMA(log_consumo ~ pdq(0:2, 1, 0:2) + PDQ(0:2, 1, 0:2, 12),
                           approximation = T,
                           stepwise = F,
                           ic = 'aicc'))

models_sarima_br |> report()

Series: log_consumo 
Model: ARIMA(2,1,2)(1,1,1)[12] 

Coefficients:
          ar1      ar2     ma1     ma2    sar1     sma1
      -1.1558  -0.9948  1.1539  0.9653  0.1093  -0.8104
s.e.   0.0097   0.0134  0.0196  0.0463  0.1030   0.0707

sigma^2 estimated as 0.0003322:  log likelihood=581.98
AIC=-1149.96   AICc=-1149.44   BIC=-1126.08

O modelo de menor AICc foi um \(SARIMA(2,1,2)(1,1,1)_{[12]}\), apresentando os seguintes resíduos

models_sarima_br |> 
  gg_tsresiduals()

Vemos que tal modelo foi o primeiro a apresentar a autocorrelação do lag 12 significativamente não diferente de 0

Aplicando novamente os teste de Box-Pierce e de Ljung-Box para testar a indepência entre os resíduos

models_sarima_br |>
  augment()  |>
  features(.resid, ljung_box, lag = 12)

# A tibble: 1 x 3
  .model     lb_stat lb_pvalue
  <chr>        <dbl>     <dbl>
1 arima_test    16.2     0.180

models_sarima_br |>
  augment()  |>
  features(.resid, box_pierce, lag = 12)

# A tibble: 1 x 3
  .model     bp_stat bp_pvalue
  <chr>        <dbl>     <dbl>
1 arima_test    15.7     0.207

Vemos que ambos os teste apresentaram p-valor acima de 0.05 e portanto a hipótese de independencia não foi rejeitada

Via histograma, vemos que a média dos resíduos foi 0

Assim vemos que o modelo estimado não quebrou nenhuma pressuposição e portanto tal modelo foi selecionado para a etapa de interpretação, análise de desempenho e previsão

4 Interpretação, Desempenho e Previsão

4.1 Termos Estimados

O modelo final selecionado foi um \(SARIMA(2,1,2)(1,1,1)_{[12]}\), o modelo teve as seguintes estimativas

models_sarima_br |> 
  coef() |>
  select(-.model, -statistic) |>
  rename(Termo = 1,
         Estimativa = 2,
         ErroPadrao = 3,
         pValor = 4) |>
  rmarkdown::paged_table()

Vemos que apenas o termo autorregressivo sazonal não foi significativamente diferente de 0 dado uma significância de 5%

4.2 Desempenho na Série Observada

models_sarima_br |>
  augment() |>
  ggplot(aes(x = ano_mes)) +
  geom_line(aes(y = log_consumo, colour = "Dados")) +
  geom_line(aes(y = .fitted, colour = "Valor Ajustado")) +
  scale_colour_manual(
  values = c(Dados = "black", "Valor Ajustado" = "#D55E00")
  ) +
  labs(y = "Log de Conusmo de Energia",
  title = "Modelo ARIMA Sazonal") +
  guides(colour = guide_legend(title = "Series")) +
  theme_minimal() +
  theme(legend.position = "bottom")

Vemos que o modelo estimado obteve um ótimo ajuste sobre os dados, sendo capaz de se ajustar com os picos sazonais e com os valores discrepantes dos anos de 2009 e 2020

4.3 Previsão

O pacote Fable utilizado para a modelagem é robusto na questão de previsão de valores a partir de dados tranformados, ou seja, previsões são feitas baseadas na série transformada e funções de transformações inversas são aplicadas para a previsão da série original, no caso desse trabalho: a série exponenciada

O modelo apresentou a seguinte previsão para 1 ano

 models_sarima_br|>
  forecast(h =12) |>
  autoplot(consumo_ts_br) +
  theme_minimal()

Observando a parte estimada

 models_sarima_br|>
  forecast(h =12) |>
  autoplot(consumo_ts_br |>
             filter(ano_mes > lubridate::ymd('2019-12-31') |>
                      tsibble::yearmonth() 
                    )
           ) +
  theme_minimal() 

Vemos que o primeiro semestre de 2024 terá um pico, chegando a \(4.5 \times10^{7}\)MWh de consumo. A partir de Agosto, o consumo sefrerá uma pequena queda, chegando aos níveis de Janeiro de 2023

5 Conclusão

Esse trabalho teve objetivo estudar o consumo de energia elétrica no Brasil.

Durante as análises, vimos que o consumo apresentou forte crescimento entre os anos de 2004 até 2023. Além disso certos momentos atípicos também foram estudados, como as quedas nos anos de 2009 e 2020.

Utilizando técnicas da teoria de sérias temporais, vimos que o consumo possui forte sazonalidade anual, significando que os valores de um mês X possui forte correlação com os valores de um mês X+12.

Para a modelagem dos dados, 3 tipos de modelos foram testados: ARIMA, ARIMA com Sazonalidade Determinística e ARIMA Sazonal, onde apenas o último desses citados foi capaz de modelar a parcela sazonal com sucesso. Tal fato foi visto a partir de uma análise de resíduos

Ao utilizar tal modelo para previsão, o primeiro semestre de 2024 apresentou um pico de consumo, onde o nível de consumo foi normalizado a partir de Agosto de 2024, chegando a níveis semelhantes a Janeiro de 2023

6 Referências

• HYNDMAN, R. J.; ATHANASOPOULOS, G. Forecasting: principles and practice. [S.l.]: OTexts, 2018. https://otexts.com/fpp3/

• MORETTIN, P. A.; TOLOI, C. M. Análise de śeries temporais: modelos lineares univariados. [S.l.]: Editora Blucher, 2018

• BOX, G. E. et al. Time series analysis: forecasting and control. [S.l.]: John Wiley & Sons, 2015.

• BREWER, Mark J.; BUTLER, Adam; COOKSLEY, Susan L. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. Methods in Ecology and Evolution. [S.l.]: Biometrika, 2016

• SEN, Liew Khim et al. The performance of AICC as an order selection criterion in ARMA time series models. [S.l.]: Pertanika Journal of Science and Technology, 2002