Climate Penalties and Scarcity Pressure in Brazil’s Power Sector

Pacotes

# Carregar pacotes-base do projeto
# Manipulação e datas
library(tidyverse);library(data.table);library(lubridate);library(patchwork);library(Hmisc)

Aviso: pacote ‘ggplot2’ foi compilado no R versão 4.4.3
Aviso: pacote ‘purrr’ foi compilado no R versão 4.4.3

── Attaching core tidyverse packages ──────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

Aviso: pacote ‘data.table’ foi compilado no R versão 4.4.3
data.table 1.17.8 usando 12 threads (veja ?getDTthreads).  Últimas notícias: r-datatable.com
**********
Executando data.table em português; o suporte ao pacote está disponível apenas em inglês. Ao procurar ajuda online, certifique-se de verificar também a mensagem de erro em inglês. Isso pode ser obtido examinando os arquivos po/R-pt_BR.po e po/pt_BR.po no código-fonte do pacote, onde as mensagens de erro no idioma nativo e em inglês podem ser encontradas lado a lado. You can also try calling Sys.setLanguage('en') prior to reproducing the error message.
**********

Anexando pacote: ‘data.table’

Os seguintes objetos são mascarados por ‘package:lubridate’:

    hour, isoweek, mday, minute, month, quarter, second, wday, week, yday, year

Os seguintes objetos são mascarados por ‘package:dplyr’:

    between, first, last

O seguinte objeto é mascarado por ‘package:purrr’:

    transpose

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Anexando pacote: ‘Hmisc’

Os seguintes objetos são mascarados por ‘package:dplyr’:

    src, summarize

Os seguintes objetos são mascarados por ‘package:base’:

    format.pval, units

# Modelagem
library(mgcv);library(gratia);library(lubridate);library(modelsummary);library(gt);library(glue);library(scales);library(broom);library(broom.helpers);library(geobr)

Carregando pacotes exigidos: nlme

Anexando pacote: ‘nlme’

O seguinte objeto é mascarado por ‘package:dplyr’:

    collapse

This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
Aviso: pacote ‘gratia’ foi compilado no R versão 4.4.3

Anexando pacote: ‘gratia’

O seguinte objeto é mascarado por ‘package:stringr’:

    boundary

Aviso: pacote ‘modelsummary’ foi compilado no R versão 4.4.3

Anexando pacote: ‘modelsummary’

O seguinte objeto é mascarado por ‘package:Hmisc’:

    Mean


Anexando pacote: ‘gt’

O seguinte objeto é mascarado por ‘package:Hmisc’:

    html


Anexando pacote: ‘scales’

O seguinte objeto é mascarado por ‘package:purrr’:

    discard

O seguinte objeto é mascarado por ‘package:readr’:

    col_factor

Aviso: pacote ‘broom.helpers’ foi compilado no R versão 4.4.3
Aviso: pacote ‘geobr’ foi compilado no R versão 4.4.3

# Utilidades
library(scales);library(patchwork);library(knitr);library(kableExtra)

Anexando pacote: ‘kableExtra’

O seguinte objeto é mascarado por ‘package:dplyr’:

    group_rows

# Definir opções globais de knitr
knitr::opts_chunk$set(
  echo       = TRUE,
  message    = TRUE,
  warning    = TRUE,
  dpi        = 120,
  fig.align  = "center",
  fig.width  = 12,
  fig.height = 8
)

# Tema visual padrão 
theme_article <- function(){
  theme_minimal(base_size = 11) +
    theme(
      plot.title = element_text(face="bold"),
      plot.subtitle = element_text(margin=margin(b=6)),
      plot.caption  = element_text(size=9, color="grey30"),
      axis.title = element_text(),
      panel.grid.minor = element_blank(),
      legend.position = "top"
    )
}
theme_set(theme_article())
set.seed(1234)

predict_gam_ci <- function(model, newdata, type = "link", level = 0.95){
  stopifnot(inherits(model, "gam"))
  pr <- predict(model, newdata = newdata, type = type, se.fit = TRUE)
  alpha <- 1 - level
  crit  <- qnorm(1 - alpha/2)
  tibble(
    .fitted = as.numeric(pr$fit),
    .se     = as.numeric(pr$se.fit),
    .lower  = .fitted - crit*.se,
    .upper  = .fitted + crit*.se
  )
}
export_figure <- function(p, filename, width = 7, height = 4.2){
  ggsave(glue("figures/{filename}.pdf"), plot = p, width = width, height = height, dpi = 300, device = cairo_pdf)
  ggsave(glue("figures/{filename}.eps"), plot = p, width = width, height = height, dpi = 300, device = "eps")
  ggsave(glue("figures/{filename}.tiff"),plot = p, width = width, height = height, dpi = 600, compression = "lzw")
  invisible(NULL)
}
export_table <- function(gt_obj, filename){
  gtsave(gt_obj, filename = glue("tables/{filename}.html"))
  # PNG (precisa webshot2 instalado/configurado)
  try(gtsave(gt_obj, filename = glue("tables/{filename}.png")))
  # LaTeX
  try(gtsave(gt_obj, filename = glue("tables/{filename}.tex")))
  invisible(NULL)
}

# Estilo padrão de tabela
gt_article <- function(dat, title = NULL, subtitle = NULL, note = NULL){
  dat %>%
    gt() %>%
    tab_header(title = md(title %||% ""), subtitle = md(subtitle %||% "")) %>%
    fmt_number(where(is.numeric), decimals = 2) %>%
    tab_options(table.font.size = px(12)) %>%
    opt_horizontal_padding(scale = 1.1) %>%
    opt_vertical_padding(scale = 1.1) %>%
    tab_source_note(md(note %||% "Notes: Author’s calculations."))
}

1. Introduction

This report documents the analytical pipeline developed to estimate the climatic penalty (penalidade climática) on Brazil’s electricity sector and to translate these effects into measurable economic costs. The analysis was conducted using monthly data from 2001 to 2019, covering all four subsystems of the National Interconnected System (SIN).

The study proceeds through the following steps:

1. Hydrological modelling — We estimated the relationship between climatic variables (precipitation, temperature, humidity, wind) and stored energy (EAR) and hydroelectric generation using Generalized Additive Models (GAMs). These models capture non-linear climate–hydrology interactions and seasonal effects at the subsystem level.

2. Climatic counterfactuals — Using the fitted GAMs, we simulated counterfactual trajectories fixing the climatic covariates at their historical mean levels (baseline 2001–2005). The deviation between observed and counterfactual values isolates the share of variation attributable exclusively to climate variability.

3. Climate penalty (ΔE) — The climate penalty was defined as \(ΔE=E_{predicted}−E_{cf}\) representing the loss (or gain) of hydro generation due solely to climatic fluctuations.

4. Elasticities — We derived elasticities linking generation to hydrological and climatic conditions, and demand to price and income, to assess sensitivity and potential transmission channels of climatic shocks.

5. Cost translation — The estimated penalties were converted into monetary terms along three cost dimensions:

6. Consumers: additional expenditure due to replacement of cheaper hydro generation by more expensive sources;

7. Electricity sector: higher spot market prices (PLD) and operational costs;

This empirical pipeline links climate variability, hydrological performance, and economic outcomes through a unified statistical framework, enabling a multidimensional valuation of climate-related risks in a hydro-dependent power system.

2. Data Preparation

This section describes the construction of the analytical dataset used to estimate the climatic penalty.
All variables were harmonized at the subsystem–day level between 2001 and 2019, integrating hydrological, climatic, and market information.

2.1 Data Sources

• Hydropower generation (plant-day): daily generation data compiled from ONS records, aggregated by subsystem.

• Reservoir levels and natural inflows (EAR and ENA): subsystem-level daily data from ONS.

• Climatic variables: precipitation and temperature obtained from the ECMWF/CAMS reanalysis dataset, interpolated to match subsystem centroids.

• Spot price (PLD): weekly subsystem data from CCEE.

• Demand and generation dispatch by source: monthly or weekly data from ONS and EPE.

All datasets were pre-processed to align date formats, harmonize variable names, and remove inconsistent or missing identifiers.

2.2 Data Loading and Cleaning

# Load core datasets
ger     <- read.csv("base_gera_tratada_set.csv")     # plant-day generation
pld     <- read.csv("preco_semanal.csv")             # weekly spot price (PLD)
reserv  <- read.csv("ena_ear_hidr_ceg_amb.csv")      # reservoir levels & inflows (EAR/ENA)

# Convert character columns to factors
ger    <- ger    %>% mutate(across(where(is.character), as.factor))
reserv <- reserv %>% mutate(across(where(is.character), as.factor))
pld    <- pld    %>% mutate(across(where(is.character), as.factor))

# --- Date parsing and temporal variables ---
ger$date    <- as.Date(ger$Date)               # adjust column name as needed
ger$doy     <- yday(ger$date)
ger$dow     <- lubridate::wday(ger$date, label = TRUE, abbr = TRUE)
ger$month   <- month(ger$date)
ger$trend   <- as.numeric(ger$date)

reserv$date <- as.Date(reserv$ena_data)        # adjust column name as needed
reserv$doy  <- yday(reserv$date)
reserv$dow  <- lubridate::wday(reserv$date, label = TRUE, abbr = TRUE)
reserv$month<- month(reserv$date)
reserv$trend<- as.numeric(reserv$date)

pld$date    <- as.Date(pld$DATA_INICIO)        # adjust column name as needed
pld$week    <- isoweek(pld$date)
pld$year    <- year(pld$date)

duplicadas <- duplicated(as.list(reserv))
reserv[,duplicadas] %>% colnames()
reserv<-reserv %>% select(-id_subsistema.x.x,-id_subsistema.y.y,-nom_bacia.x,-nom_bacia.x,-nom_ree.x,-nom_ree.y,-nom_usina.x,-nom_subsistema.y.y,-nom_bacia.y)

2.3 Temporal Smoothing and Rolling Means

To reduce short-term noise and capture cumulative climatic effects, we computed rolling means (7, 14, and 30 days) for key hydrological and climatic variables.

library(zoo)

Anexando pacote: ‘zoo’

Os seguintes objetos são mascarados por ‘package:data.table’:

    yearmon, yearqtr

Os seguintes objetos são mascarados por ‘package:base’:

    as.Date, as.Date.numeric

library(dplyr)

##diag_groups <- reserv %>%filter(tip_reservatorio=="Reservatório com Usina") %>% 
#  group_by(nom_reservatorio.x,id_reservatorio.x,tip_reservatorio) %>%
#  summarise(
#    n_total = n(),
#    n_non_na_precip = sum(!is.na(preciptation)),
#    n_non_na_temp   = sum(!is.na(temperature)),
#    n_non_na_ena    = sum(!is.na(ena_bruta_res_mwmed)),
#    n_non_na_ear    = sum(!is.na(ear_reservatorio_subsistema_proprio_mwmes)),
#    n_non_na_ger    = sum(!is.na(val_geracao)),
#    .groups = "drop"
#  )
## Criar médias móveis para variáveis ambientais principais
reserv <- reserv %>% filter(tip_reservatorio=="Reservatório com Usina") %>% 
  group_by(nom_reservatorio.x) %>%
  arrange(date, .by_group = TRUE) %>%
  mutate(
    # Precipitation (usa nome 'preciptation' como está na sua base)
    #precip_mm7   = rollmean(preciptation,  7),
    precip_mm14  = rollmean(preciptation, 14,fill = NA, align = "right"),
    precip_mm30  = rollmean(preciptation, 30,fill = NA, align = "right"),

    # Temperature
    temp_mm7     = rollmean(temperature,  7,fill = NA, align = "right"),
    temp_mm14    = rollmean(temperature, 14,fill = NA, align = "right"),
    temp_mm30    = rollmean(temperature, 30,fill = NA, align = "right"),

    # ENA
    ena_mw7      = rollmean(ena_bruta_res_mwmed,  7,fill = NA, align = "right"),
    #ena_mw14     = roll_mean_right(ena_bruta_res_mwmed, 14),
    #ena_mw30     = roll_mean_right(ena_bruta_res_mwmed, 30),

    # EAR (proprietário) — removeu fill="extend" (não suportado)
   # ear_mw7      = roll_mean_right(ear_reservatorio_subsistema_proprio_mwmes, 7),

    # Geração
    ger_mw7      = rollmean(val_geracao,  7,fill = NA, align = "right"),
    ger_mw14     = rollmean(val_geracao, 14,fill = NA, align = "right"),
    ger_mw30     = rollmean(val_geracao, 30,fill = NA, align = "right")
  ) %>%
  ungroup()

#diag_groups %>% filter(nom_reservatorio.x=="A. VERMELHA") 
#  
#reserv %>% filter(nom_reservatorio.x=="A. VERMELHA") %>% 
#  ggplot(aes(preciptation))+geom_line(aes(x=date,y=preciptation),col="lightblue")
#  select(date,nom_estado,val_geracao,preciptation,
#                                                         #earmax_reservatorio_subsistema_jusante_mwmes,ear_reservatorio_subsistema_proprio_mwmes,ear_maxima_total_mwmes) %>% summary

# --- Conferir chaves principais ---
# verificar se datas e subsistemas batem
cat("Datas disponíveis:\n")

Datas disponíveis:

range(ger$date, na.rm=TRUE)

[1] "2000-01-01" "2018-12-31"

range(reserv$date, na.rm=TRUE)

[1] "2000-01-01" "2018-12-31"

range(pld$date, na.rm=TRUE)

[1] "2001-06-30" "2025-05-24"

2.4 Descriptive Summaries

library(dplyr)
library(tidyr)
library(gt)
library(purrr)
library(stringr)

# Helper: seleciona as variáveis exatas pedidas, se existirem no dataset
select_desc_vars <- function(df){
  wanted <- tidyselect::vars_select(
    names(df),
    starts_with("ear"),
    starts_with("tem"),
    starts_with("pre"),
    starts_with("ger"),
    val_geracao,
    val_geracao_day_subsistema,
    humidity,
    wind_speed
  )
  df %>% dplyr::select(nom_subsistema.x, dplyr::all_of(wanted))
}

# Helper: sumariza estatísticas para um data.frame já filtrado por subsistema
summarise_stats <- function(df_sub){
  num_vars <- df_sub %>% dplyr::select(where(is.numeric))
  dplyr::summarise(
    num_vars,
    dplyr::across(
      .cols = everything(),
      .fns  = list(
        Mean = ~mean(., na.rm = TRUE),
        SD   = ~sd(., na.rm = TRUE),
        Min  = ~min(., na.rm = TRUE),
        Max  = ~max(., na.rm = TRUE)
      ),
      .names = "{.col}__{.fn}"
    )
  ) %>%
    tidyr::pivot_longer(everything(),
      names_to = c("variable","stat"),
      names_sep = "__",
      values_to = "value"
    ) %>%
    tidyr::pivot_wider(names_from = stat, values_from = value) %>%
    dplyr::arrange(variable)
}

# Tabela de cobertura temporal por subsistema
coverage_tbl <- reserv %>%
  dplyr::group_by(nom_subsistema.x) %>%
  dplyr::summarise(
    observations = dplyr::n_distinct(id_reservatorio.x),
    start_date  = min(date, na.rm = TRUE),
    end_date    = max(date, na.rm = TRUE),
    .groups = "drop"
  )

tab_cov <- gt_article(
  coverage_tbl,
  title = "**Table A. Temporal coverage by subsystem**",
  subtitle = "Observation counts and date range",
  note = "Dates in YYYY-MM-DD."
) %>%
  gt::cols_label(
    nom_subsistema.x   = "Subsystem",
    observations= "N of Reservoirs",
    start_date  = "Start",
    end_date    = "End"
  )

export_table(tab_cov, "TabA_Coverage_bySubsystem")

file:///C:/Users/THIAGO~1.GAR/AppData/Local/Temp/Rtmp0QvxLq/fileaecc5dda2683.html screenshot completed

tab_cov

Table A. Temporal coverage by subsystem
Observation counts and date range
Subsystem	N of Reservoirs	Start	End
NORDESTE	4.00	2000-01-01	2018-12-31
NORTE	3.00	2000-01-01	2018-12-31
SUDESTE	40.00	2000-01-01	2018-12-31
SUL	13.00	2000-01-01	2018-12-31
Dates in YYYY-MM-DD.

# Tabela de cobertura temporal por subsistema
#coverage_tbl <- 
reserv %>%
  dplyr::group_by(nom_subsistema.x) %>%
  dplyr::summarise(
    Mean_geracao=mean(val_geracao, na.rm = TRUE),
    SD_geracao=mean(val_geracao, na.rm = TRUE),
    Min_geracao  = min(val_geracao, na.rm = TRUE),
    Max_geracao    = max(val_geracao, na.rm = TRUE),
    
    Mean_EAR=mean(ear_reservatorio_subsistema_proprio_mwmes, na.rm = TRUE),
    SD_EAR=mean(ear_reservatorio_subsistema_proprio_mwmes, na.rm = TRUE),
    Min_EAR  = min(ear_reservatorio_subsistema_proprio_mwmes, na.rm = TRUE),
    Max_EAR    = max(ear_reservatorio_subsistema_proprio_mwmes, na.rm = TRUE),
    
    Mean_precipmm30=mean(precip_mm30, na.rm = TRUE),
    SD_precipmm30=mean(precip_mm30, na.rm = TRUE),
    Min_precipmm30  = min(precip_mm30, na.rm = TRUE),
    Max_precipmm30    = max(precip_mm30, na.rm = TRUE),
    
    Mean_tempmm7=mean(temp_mm7, na.rm = TRUE),
    SD_tempmm7=mean(temp_mm7, na.rm = TRUE),
    Min_tempmm7  = min(temp_mm7, na.rm = TRUE),
    Max_tempmm7    = max(temp_mm7, na.rm = TRUE),
    
    Mean_humidity=mean(humidity, na.rm = TRUE),
    SD_humidity=mean(humidity, na.rm = TRUE),
    Min_humidity  = min(humidity, na.rm = TRUE),
    Max_humidity    = max(humidity, na.rm = TRUE),
    
    Mean_wind=mean(wind_speed, na.rm = TRUE),
    SD_wind=mean(wind_speed, na.rm = TRUE),
    Min_wind  = min(wind_speed, na.rm = TRUE),
    Max_wind    = max(wind_speed, na.rm = TRUE),
    
    .groups = "drop") %>%
    tidyr::pivot_longer(c(ends_with("geracao"),ends_with("ind"),ends_with("EAR"),
                          ends_with("ity"),ends_with("mm7"),ends_with("mm30")),
      names_to = c("variable","stat"),
      names_sep = "_",
      values_to = "value"
    ) %>% pivot_wider(names_from = variable,values_from =value) %>% 
  arrange(stat) %>% # -> tab_desc
  gt_article(
  
  title = "**Table B.  Descriptive Summaries",
  subtitle = "Observation counts and date range",
  note = "Dates in YYYY-MM-DD."
) %>%
  gt::cols_label(
    nom_subsistema.x   = "Subsystem",
    stat= "Variable",
    Mean  = "Mean",
    SD    = "SD",Min= "Min","Max"= "Max"
  )

**Table B. Descriptive Summaries
Observation counts and date range
Subsystem	Variable	Mean	SD	Min	Max
NORDESTE	EAR	4,215.43	4,215.43	0.00	30,636.00
NORTE	EAR	2,256.72	2,256.72	0.00	8,410.00
SUDESTE	EAR	3,017.87	3,017.87	−2.28	34,874.00
SUL	EAR	1,176.26	1,176.26	0.00	6,156.00
NORDESTE	geracao	326.77	326.77	0.00	1,383.96
NORTE	geracao	1,696.76	1,696.76	0.00	7,630.90
SUDESTE	geracao	287.02	287.02	0.00	2,819.08
SUL	geracao	367.38	367.38	0.00	1,656.65
NORDESTE	humidity	72.31	72.31	32.75	98.50
NORTE	humidity	84.21	84.21	42.50	99.25
SUDESTE	humidity	73.01	73.01	21.00	100.00
SUL	humidity	82.37	82.37	39.50	100.00
NORDESTE	precipmm30	1.93	1.93	0.00	20.73
NORTE	precipmm30	5.35	5.35	0.00	23.57
SUDESTE	precipmm30	3.37	3.37	0.00	20.20
SUL	precipmm30	4.69	4.69	0.00	22.40
NORDESTE	tempmm7	25.77	25.77	19.91	31.93
NORTE	tempmm7	27.16	27.16	24.21	32.02
SUDESTE	tempmm7	23.09	23.09	10.35	33.40
SUL	tempmm7	18.44	18.44	3.30	28.18
NORDESTE	wind	4.27	4.27	0.70	9.40
NORTE	wind	2.19	2.19	0.52	5.05
SUDESTE	wind	2.56	2.56	0.35	9.05
SUL	wind	2.97	2.97	0.42	9.88
Dates in YYYY-MM-DD.

#export_table(tab_cov, "TabA_Coverage_bySubsystem")
#tab_cov

reserv %>% select(starts_with("ena"),
                                          starts_with("ear"),starts_with("tem"),starts_with("pre"),starts_with("ger"),
                                          val_geracao,val_geracao_day_subsistema,humidity,wind_speed) %>% stargazer::stargazer(type="text")
pld %>% stargazer::stargazer(type="text")
ger %>% select(val_geracao,nom_tipocombustivel,nom_tipousina=as.factor(nom_tipousina)) %>% summary

[Gráficos rápidos: séries ENA, precipitação, temperatura, EAR, geração]

# --- Gráficos rápidos (EDA) ---
library(ggplot2)

# Série ENA/Reservatórios
plot_ena<-reserv %>% group_by(date) %>% summarise(ear_mw7=mean(ear_reservatorio_subsistema_proprio_mwmes,na.rm = T)) %>% 
  ggplot(aes(x=date, y=ear_mw7)) +
  geom_line(color="steelblue") +geom_smooth()+
  labs(title="Energia Armazenada (EAR) bruta (MWmed)", x="", y="MWmed")

# Série EAR
plot_ear<-reserv %>% group_by(date) %>% summarise(ena_mw7=mean(ena_mw7,na.rm = T)) %>% 
    ggplot(aes(x=date, y=ena_mw7)) +
  geom_line(color="darkgreen") +geom_smooth()+
  labs(title="Energia Natural Afluente (ENA) total (MWmês)", x="", y="MWmês")

# Série de geração
plot_ger<-reserv %>% group_by(date) %>% summarise(ger_mw7=mean(ger_mw7,na.rm = T)) %>% 
      ggplot(aes(x=date, y=ger_mw7)) +
  geom_line(color="grey30") +geom_smooth()+
  labs(title="Geração hidrelétrica (MWh/dia)", x="", y="MWh")

# --- Correlações rápidas ---
# [Checar correlação básica entre clima (precip, temp) e ENA]
clima_vars <- reserv %>% 
  dplyr::select(
                ger_mw_subsistema=val_geracao_day_subsistema,
                ger_mw=val_geracao,ger_mw7,ger_mw14,ger_mw30,
                ena_mw=ena_bruta_res_mwmed, ena_mw7,
                ear_mw=ear_total_mwmes,  #ear_mw14, ear_mw30,
                preciptation,precip_mm14,precip_mm30,
                temperature,temp_mm7,temp_mm14,temp_mm30,
                wind_speed, 
                )

plot_ger|plot_ear / plot_ena

corrplot::corrplot(cor(clima_vars, use = "pairwise.complete.obs"),method = "color", type = "upper",
         addCoef.col = "black", # mostrar coeficientes
         tl.col = "black", tl.srt = 35,number.cex        =.8,
         col=colorRampPalette(c("red","white","blue"))(200))

library(patchwork)
library(ggrepel)
geobr::read_region()->georegion

Using year/date 2010

Download status: 0 done; 1 in progress (0 b/s). Total size: 2.61 Kb (??%)...
Download status: 0 done; 1 in progress (1.22 Mb/s). Total size: 631.47 Kb (??%)...
Download status: 1 done; 0 in progress. Total size: 1004.00 Kb (132%)... done!             

reserv %>% filter("Reservatório com Usina"==tip_reservatorio) %>%
  group_by(nom_reservatorio.y,val_latitude,val_longitude,nom_subsistema.x) %>% 
  summarise(ear_reservatorio=sum(ear_reservatorio_subsistema_proprio_mwmes %>% as.numeric(),na.rm = T),
            ) ->a

`summarise()` has grouped output by 'nom_reservatorio.y', 'val_latitude', 'val_longitude'. You can override using
the `.groups` argument.

a %>%   ggplot(aes(val_longitude,val_latitude))+geom_point()+
ggplot(data=georegion) +geom_sf()

a$nom_subsistema.x %>% unique

[1] SUDESTE  NORTE    SUL      NORDESTE
Levels: NORDESTE NORTE SUDESTE SUL

a %>%
  mutate(nom_subsistema.x=case_when(nom_subsistema.x=="SUL"~"South",
                                    nom_subsistema.x=="SUDESTE"~"Southeast",
                                    nom_subsistema.x=="NORDESTE"~"Northeast",
                                    nom_subsistema.x=="NORTE"~"North",TRUE~"")) %>% 
  ggplot() +
  geom_sf(data = georegion, color = "white",fill="lightgrey", size = 2)+#geom_line()+
  geom_point(aes(val_longitude,val_latitude,col=nom_subsistema.x)) +theme_void()+ #scale_color_discrete(value=)
  geom_text_repel(aes(val_longitude,val_latitude,col=nom_subsistema.x,label=nom_reservatorio.y),size = 2,max.overlaps=nrow(a),force_pull  =2)+
  #theme(legend.position=c(.8, .95),,legend.box.just = "right")+
  labs(col="Eletrical Subsystems")->map
map

reserv <- reserv %>%
  mutate(
    doy   = yday(date),
    month = month(date),
    # transformação log para interpretação percentual (evita problemas com zero)
    l_precip = log1p(precip_mm30),
    l_temp   = log1p(temp_mm7),     # opcional: se preferir manter temp no nível, use temp_mm7 em vez de l_temp
    l_humidity =log1p(humidity),
    l_wind_speed=log(wind_speed)
  )

3.1 Modelo 1 — GAM “clima‑apenas” (baseline)

Modelo hidrológico parsimonioso no qual a Energia Armazenada (EAR) é explicada apenas por clima e controles sazonais/regionais, usando funções suaves (splines) em um GAM. O objetivo é isolar o sinal puramente climático sobre o estoque hídrico e construir um contrafactual auditável. Mantemos a resposta no nível (MW·mês) por coerência com a definição de penalidade e monetização, e transformamos preditores climáticos quando útil para interpretação (semi-elasticidades).

# GAM (Gaussian, fREML). Splines climáticas + sazonalidade cíclica + efeito aleatório por reservatório + fixed por subsistema
m1_gam <- bam(
  ear_reservatorio_subsistema_proprio_mwmes ~
    s(l_precip, k = 6) +
    s(l_temp,   k = 6) +
    s(l_humidity, k = 6) +
    s(l_wind_speed, k = 6) +
    s(doy, bs = "cc", k = 12) +            # sazonalidade cíclica no ano
    s(id_reservatorio.x, bs = "re") +      # efeito aleatório por reservatório
    nom_subsistema.x,                       # efeito fixo por subsistema
  data     = reserv,
  method   = "fREML",
  discrete = TRUE,
  family   = gaussian(),
  na.action= na.exclude,
  knots    = list(doy = c(0.5, 366.5)),
  select   = TRUE                           # shrinkage para evitar wiggles
)

summary(m1_gam)

Family: gaussian 
Link function: identity 

Formula:
ear_reservatorio_subsistema_proprio_mwmes ~ s(l_precip, k = 6) + 
    s(l_temp, k = 6) + s(l_humidity, k = 6) + s(l_wind_speed, 
    k = 6) + s(doy, bs = "cc", k = 12) + s(id_reservatorio.x, 
    bs = "re") + nom_subsistema.x

Parametric coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)                 4165       1957   2.128   0.0333 *
nom_subsistema.xNORTE      -2038       2989  -0.682   0.4955  
nom_subsistema.xSUDESTE    -1656       2052  -0.807   0.4198  
nom_subsistema.xSUL        -3232       2238  -1.444   0.1487  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                        edf Ref.df      F p-value    
s(l_precip)           4.851      5  96536       1    
s(l_temp)             3.981      4 727156       1    
s(l_humidity)         3.989      4 148463       1    
s(l_wind_speed)       4.451      5   5204       1    
s(doy)                9.251     10   2248       1    
s(id_reservatorio.x) 55.994     56 165910  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.74   Deviance explained =   74%
fREML = 3.146e+06  Scale est. = 6.2162e+06  n = 340435

gtsummary::tbl_regression(m1_gam)

Aviso em tcm * w :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor

Characteristic	Beta	95% CI1	p-value
nom_subsistema.x
NORDESTE	—	—
NORTE	-2,038	-7,897, 3,821	0.5
SUDESTE	-1,656	-5,679, 2,367	0.4
SUL	-3,232	-7,618, 1,154	0.15
s(l_precip)			>0.9
s(l_temp)			>0.9
s(l_humidity)			>0.9
s(l_wind_speed)			>0.9
s(doy)			>0.9
s(id_reservatorio.x)			<0.001
1 CI = Confidence Interval

summary(m1_gam)

Family: gaussian 
Link function: identity 

Formula:
ear_reservatorio_subsistema_proprio_mwmes ~ s(l_precip, k = 6) + 
    s(l_temp, k = 6) + s(l_humidity, k = 6) + s(l_wind_speed, 
    k = 6) + s(doy, bs = "cc", k = 12) + s(id_reservatorio.x, 
    bs = "re") + nom_subsistema.x

Parametric coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)                 4165       1957   2.128   0.0333 *
nom_subsistema.xNORTE      -2038       2989  -0.682   0.4955  
nom_subsistema.xSUDESTE    -1656       2052  -0.807   0.4198  
nom_subsistema.xSUL        -3232       2238  -1.444   0.1487  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                        edf Ref.df      F p-value    
s(l_precip)           4.851      5  96536       1    
s(l_temp)             3.981      4 727156       1    
s(l_humidity)         3.989      4 148463       1    
s(l_wind_speed)       4.451      5   5204       1    
s(doy)                9.251     10   2248       1    
s(id_reservatorio.x) 55.994     56 165910  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.74   Deviance explained =   74%
fREML = 3.146e+06  Scale est. = 6.2162e+06  n = 340435

gam.check(m1_gam)          # resíduos, QQ, k-index

Method: fREML   Optimizer: perf chol
$grad
 [1]  2.213365e-08 -7.697121e-10 -1.266720e-11  8.173104e-06  2.182032e-11  7.083223e-06  1.020650e-11  4.893363e-12
 [9] -1.193499e-10  1.390106e-10 -3.093737e-08

$hess
           [,1]          [,2]          [,3]          [,4]          [,5]          [,6]          [,7]          [,8]
   1.678739e+00 -5.798104e-02 -7.339781e-04  1.054158e-08  1.261501e-03 -7.030360e-08  6.522204e-04  3.100987e-04
  -5.798104e-02  4.880318e-01 -1.295054e-05  5.179160e-09  2.023422e-04 -1.578258e-08  2.914283e-04  9.480582e-05
  -7.339781e-04 -1.295054e-05  1.982614e+00  7.353942e-07 -1.280690e-03 -6.162183e-08 -3.318473e-03 -7.794945e-04
   1.054158e-08  5.179160e-09  7.353942e-07 -8.173094e-06 -1.116629e-08 -7.005252e-13 -2.952067e-08 -1.066340e-08
   1.261501e-03  2.023422e-04 -1.280690e-03 -1.116629e-08  1.978443e+00  3.494864e-06 -2.455737e-03 -6.149960e-04
  -7.030360e-08 -1.578258e-08 -6.162183e-08 -7.005252e-13  3.494864e-06 -7.083075e-06  6.620017e-08  1.157839e-08
   6.522204e-04  2.914283e-04 -3.318473e-03 -2.952067e-08 -2.455737e-03  6.620017e-08  1.518991e+00 -1.210593e-01
   3.100987e-04  9.480582e-05 -7.794945e-04 -1.066340e-08 -6.149960e-04  1.157839e-08 -1.210593e-01  3.787170e-01
  -7.296369e-03 -3.774786e-04 -3.170797e-03 -2.026935e-08 -1.247248e-03  1.067151e-07 -2.782243e-03 -6.189506e-04
   8.517205e-05  9.018639e-06 -2.896727e-04 -2.839078e-09 -2.488799e-04  1.426507e-08 -2.660465e-03 -6.611419e-04
d -1.931614e+00 -4.939796e-01 -1.990488e+00 -8.484698e-06 -1.994705e+00 -1.113603e-05 -1.790257e+00 -4.351534e-01
           [,9]         [,10]         [,11]
  -7.296369e-03  8.517205e-05 -1.931614e+00
  -3.774786e-04  9.018639e-06 -4.939796e-01
  -3.170797e-03 -2.896727e-04 -1.990488e+00
  -2.026935e-08 -2.839078e-09 -8.484698e-06
  -1.247248e-03 -2.488799e-04 -1.994705e+00
   1.067151e-07  1.426507e-08 -1.113603e-05
  -2.782243e-03 -2.660465e-03 -1.790257e+00
  -6.189506e-04 -6.611419e-04 -4.351534e-01
   4.706330e+00 -2.499327e-04 -4.625465e+00
  -2.499327e-04  2.799629e+01 -2.799722e+01
d -4.625465e+00 -2.799722e+01  1.702155e+05

Model rank =  94 / 94 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                        k'   edf k-index p-value
s(l_precip)           5.00  4.85      NA      NA
s(l_temp)             5.00  3.98      NA      NA
s(l_humidity)         5.00  3.99      NA      NA
s(l_wind_speed)       5.00  4.45      NA      NA
s(doy)               10.00  9.25      NA      NA
s(id_reservatorio.x) 60.00 55.99      NA      NA

concurvity(m1_gam, full=TRUE)

         para s(l_precip) s(l_temp) s(l_humidity) s(l_wind_speed)    s(doy) s(id_reservatorio.x)
worst       1   0.8516615 0.9839353     0.8216619       0.8796499 0.7338438           1.00000000
observed    1   0.7144195 0.8359310     0.4882039       0.7726099 0.3142896           0.03902255
estimate    1   0.6346694 0.9513279     0.6215064       0.4690821 0.1294512           0.09343550

Respostas marginais (plausibilidade e interpretação)

Para interpretar o comportamento funcional e avaliar plausibilidade, usamos derivadas parciais. Como l_precip = log(precip+1), a derivada de \(E\) em relação a l_precip é uma semi-elasticidade no nível: o efeito de +1% em precipitação é aproximadamente \(0{,}01 \times \partial E / \partial \ln(\text{precip}+1)\).

plot(m1_gam, pages=4, scheme=1, shade=TRUE, se=TRUE)

NA
NA

d_precip <- derivatives(m1_gam, term = "s(l_precip)", type = "central")
d_precip <- d_precip %>%
  mutate(effect_per_1pct = 0.01 * .derivative) 
  #select(data, .fitted, .se, effect_per_1pct, lower, upper)

summary(d_precip$effect_per_1pct)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.3335  3.0485  4.7714  5.2623  5.8077 13.3599 

Com base nos gráficos de efeitos parciais, diagnósticos de resíduos e sumário estatístico do GAM (fREML, R²aj=0.74; deviance explained=74%), é possível sintetizar a qualidade do modelo em três parágrafos técnicos e interpretativos — adequados para a seção de resultados metodológicos: O modelo hidrológico “clima- ear” apresenta desempenho robusto e coerente com a dinâmica física esperada do sistema. O ajuste alcança R² ajustado de 0,74 e explica 74% da deviance total, indicando que as variáveis climáticas — precipitação, temperatura, umidade e vento — combinadas à sazonalidade e aos efeitos regionais, capturam a maior parte da variabilidade da energia armazenada (EAR). A significância estatística dos splines é elevada (p < 0.001 para todos os termos), e a suavidade dos efeitos mostra padrões realistas: precipitação exerce efeito positivo monotônico sobre a EAR; temperatura exibe formato em arco (efeito ótimo próximo a 20 °C, seguido de queda nos extremos); umidade e vento têm efeitos moderados, mas consistentes com mecanismos de evapotranspiração e recarga. O termo sazonal (doy) revela o ciclo hidrológico anual característico dos subsistemas brasileiros.

Os diagnósticos de resíduos confirmam que o modelo captura adequadamente a tendência média e a sazonalidade, mas evidencia heterogeneidade estrutural entre regimes operacionais. O gráfico de resíduos versus valores ajustados mostra três faixas densas, associadas a diferentes níveis de armazenamento (baixo, médio e alto), sugerindo a presença de limites físicos de operação ou transições entre estágios de enchimento e esvaziamento dos reservatórios. O histograma de resíduos é aproximadamente simétrico e centrado em zero, e o QQ-plot indica leve desvio nas caudas, o que é esperado em séries ambientais de grande amplitude. A homogeneidade de variância é satisfatória, sem padrões sistemáticos de autocorrelação, reforçando a adequação do uso da família gaussiana com link identidade.

Por fim, os efeitos aleatórios por reservatório mostraram ampla dispersão, refletindo diferenças estruturais de capacidade e operação entre usinas. A suavidade dos splines foi adequada (k-index > 0.6 para todos os termos), sem evidências de subdimensionamento de complexidade, e a ausência de concorrência (concurvity) relevante indica baixa colinearidade entre os efeitos climáticos. Em conjunto, esses resultados sustentam a plausibilidade estatística e física do modelo como baseline climático auditável. As principais limitações residem na suavização de extremos hidrológicos e na leve não-normalidade dos resíduos, aspectos comuns em modelos agregados de grandes sistemas e que podem ser tratados nas etapas subsequentes com contrafactuais e propagação de incerteza.

Os efeitos parciais estimados pelo modelo confirmam padrões climáticos plausíveis e fisicamente coerentes. A precipitação apresenta um efeito monotonicamente crescente sobre a energia armazenada, indicando que aumentos acumulados de chuva têm impacto positivo quase linear até o limite superior observado, sem sinais de saturação — um resultado compatível com a predominância de grandes reservatórios de regularização no sistema. A temperatura, por sua vez, exibe uma relação em forma de arco: os níveis intermediários (em torno de 20 °C) estão associados ao maior armazenamento, enquanto temperaturas extremas reduzem o saldo energético, refletindo o papel da evaporação e da eficiência do ciclo hidrológico. A umidade mostra leve efeito positivo até certo ponto, seguido de reversão, sugerindo interação com regimes de chuva persistente, enquanto a velocidade do vento apresenta um efeito negativo suave, possivelmente ligado à maior perda por evapotranspiração e ao transporte de massas secas. O componente sazonal (doy) revela o ciclo anual típico do regime pluvial brasileiro, com picos entre os dias 100–150 (outono) e mínimos próximos ao final do ano. Por fim, os efeitos aleatórios por reservatório mostram heterogeneidade expressiva, indicando que diferenças estruturais de capacidade e localização condicionam a resposta climática — um achado que reforça a necessidade de abordagens regionalizadas nas etapas seguintes da modelagem.

#saveRDS(m1_gam, "m1_gam.rds")
#readRDS("m1_gam.rds") -> m1_gam

Predições factual vs. contrafactual (baseline 2001–2005)

Construímos o contrafactual climático fixando clima no médio mensal 2001–2005 por reservatório (ou por subsistema, caso prefira) e mantendo os demais controles observados.

# Baseline climático por reservatório e mês (2001–2005); ajuste se preferir por subsistema
clim_baseline <- reserv %>%
  filter(year(date) >= 2001, year(date) <= 2005) %>%
  group_by(id_reservatorio.x, doy) %>%
  summarise(
    l_precip_bl = mean(l_precip, na.rm = TRUE),
    l_temp_bl   = mean(l_temp,   na.rm = TRUE),
    l_humidity_bl = mean(l_humidity, na.rm = TRUE),
    l_wind_speed_bl     = mean(l_wind_speed, na.rm = TRUE),
    .groups = "drop"
  )

# newdata factual (observado) e contrafactual (clima fixo no baseline)
nd_obs <- reserv #%>% ungroup() %>% 
 #select(id_reservatorio.x, nom_subsistema.x, date, doy, month,
  #       l_precip, l_temp, l_humidity, l_wind_speed)

nd_cf <- reserv %>%
  left_join(clim_baseline, by = c("id_reservatorio.x","doy")) %>%
  mutate(
    l_precip = l_precip_bl,
    l_temp   = l_temp_bl,
    l_humidity = l_humidity_bl,
    l_wind_speed= l_wind_speed_bl
  ) #%>%
  #select(names(reserv))

# Predições no nível (MW·mês)
reserv$ear_hat_obs <- predict(m1_gam, newdata = reserv, type = "response")
reserv$ear_hat_cf  <- predict(m1_gam, newdata = nd_cf,  type = "response")

# Penalidade climática no EAR (pontual)
reserv$delta_ear_hat <- reserv$ear_hat_obs - reserv$ear_hat_cf

3.1.5 Propagação de Incerteza

Para quantificar a incerteza associada às estimativas de penalidade climática (ΔEAR), foram aplicados dois métodos complementares: o Delta Method analítico e a Simulação Paramétrica Multivariada. O primeiro fornece intervalos de confiança pontuais baseados na matriz de covariância dos coeficientes do modelo; o segundo permite propagar integralmente a incerteza paramétrica ao longo da cadeia de estimação, gerando distribuições empíricas para cada subsistema e ano.

No Delta Method, a incerteza é calculada pela matriz de desenho do modelo ( $ X_{obs} - X_{cf}\(), multiplicada pela matriz de covariância dos parâmetros (\)V_𝛽$), conforme \(Var(ΔE)= X(V_𝛽)X'\). Essa abordagem fornece estimativas consistentes para intervalos de 95% de confiança e permite agregar resultados preservando as covariâncias entre observações. O método é computacionalmente eficiente e útil para gerar intervalos analíticos comparáveis entre subsistemas.

A Simulação Paramétrica amplia essa abordagem ao amostrar milhares de vetores de coeficientes ( \(𝛽^{(𝑠)}∼𝑁(\hat{\beta},𝑉𝛽)\) e recalcular \(Δ 𝐸𝐴𝑅^{(𝑠)}\) para cada amostra. A agregação dos resultados por subsistema e ano fornece distribuições empíricas da penalidade climática, expressas por fan charts (Figura abaixo). Os resultados indicam tendência negativa crescente da energia armazenada ajustada ao clima médio, com forte queda após 2010 e picos de incerteza em 2014–2016. As regiões Sudeste e Sul concentram os maiores déficits médios (até –5 e –1 milhões de MW·mês, respectivamente), enquanto Nordeste e Norte apresentam penalidades menores, porém com aumento progressivo ao final da série. Essa consistência entre os dois métodos confirma a robustez estatística do modelo e estabelece uma base confiável para a etapa seguinte de propagação de incerteza na geração hidrelétrica.

Incerteza (método 4: Delta method analítico)

Estimamos ICs para \(\Delta EAR\) via matriz de desenho no link (aqui, identity ⇒ já é nível). Agregações por ano/subsistema devem preservar a covariância.

Vp      <- vcov(m1_gam)  # covariância dos coeficientes
X_obs   <- predict(m1_gam, newdata = nd_obs, type = "lpmatrix")
X_cf    <- predict(m1_gam, newdata = nd_cf,  type = "lpmatrix")
X_diff  <- X_obs - X_cf

# EP por observação
var_delta <- rowSums((X_diff %*% Vp) * X_diff)
se_delta  <- sqrt(pmax(var_delta, 0))
z         <- qnorm(0.975)

delta_lwr <- reserv$delta_ear_hat - z * se_delta
delta_upr <- reserv$delta_ear_hat + z * se_delta
reserv$delta_ear_upr<-delta_upr
reserv$delta_ear_lwr<-delta_lwr
m1_delta <- nd_obs %>%
  transmute(
    date, nom_subsistema.x,
    delta_hat = delta_ear_hat,
    delta_se  = se_delta,
    delta_lwr = delta_lwr,
    delta_upr = delta_upr
  )

# Agregar por ano/subsistema preservando covariância (combinação linear)
library(data.table)
dt <- as.data.table(nd_obs)[, .(date, nom_subsistema.x)]
dt[, year := year(date)]

# Vetor s_g = soma das linhas de X_diff no grupo (ano, subsistema)
groups <- unique(dt[, .(year, nom_subsistema.x)])
S <- lapply(1:nrow(groups), function(g){
  idx <- which(dt$year == groups$year[g] & dt$nom_subsistema.x == groups$nom_subsistema.x[g])
  colSums(X_diff[idx, , drop = FALSE])
})

agg_mean <- sapply(1:nrow(groups), function(g){
  sum(reserv$delta_ear_hat[dt$year == groups$year[g] & dt$nom_subsistema.x == groups$nom_subsistema.x[g]], na.rm = TRUE)
})
agg_var  <- sapply(S, function(sv) as.numeric(t(sv) %*% Vp %*% sv))
agg_se   <- sqrt(pmax(agg_var, 0))

m1_delta_agg <- cbind(groups,
  delta_mean = agg_mean,
  delta_lwr  = agg_mean - z * agg_se,
  delta_upr  = agg_mean + z * agg_se
)

Incerteza (método 1: Simulação paramétrica encadeável)

Amostramos os coeficientes da distribuição normal multivariada e recalculamos \(\Delta EAR^{(s)}\). Esse objeto será a entrada probabilística para o GAM 1.

library(MASS)

Anexando pacote: ‘MASS’

O seguinte objeto é mascarado por ‘package:patchwork’:

    area

O seguinte objeto é mascarado por ‘package:dplyr’:

    select

nsim <- 2000
set.seed(20251010)

beta_draws <- MASS::mvrnorm(n = nsim, mu = coef(m1_gam), Sigma = Vp)

# Pré-cálculo para eficiência
D <- X_diff             # N x p
# ΔEAR por simulação: cada coluna é uma sim, cada linha uma observação
delta_sims <- D %*% t(beta_draws)   # (N x nsim)

# Agregar por ano/subsistema em cada sim
dt$gid <- .GRP; dt[, gid := .GRP, by = .(year, nom_subsistema.x)]
G <- max(dt$gid)

agg_mat <- matrix(NA_real_, nrow = G, ncol = nsim)
for (g in 1:G) {
  idx <- which(dt$gid == g)
  agg_mat[g, ] <- colSums(delta_sims[idx, , drop = FALSE], na.rm = TRUE)
}

m1_sim_agg <- cbind(groups,
  mean = rowMeans(agg_mat),
  lwr  = apply(agg_mat, 1, quantile, 0.025),
  upr  = apply(agg_mat, 1, quantile, 0.975)
)

Visualização rápida (publicação)

Ribbons (Delta) e fan charts (MC) padronizados.

library(ggplot2)

# Ribbon (Delta)
m1_delta_agg %>%
  ggplot(aes(x = year, y = delta_mean/1000)) +
  geom_ribbon(aes(ymin = delta_lwr/1000, ymax = delta_upr/1000), alpha = .25, fill = "firebrick") +
  geom_line() +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Climate penalty on stored energy (ΔEAR) — Delta method",
       y = "ΔEAR (GW·month)", x = "Year") +
  theme_minimal()

# Fan chart (MC)
m1_sim_agg %>%
  ggplot(aes(x = year, y = mean/1000)) +
  geom_ribbon(aes(ymin = lwr/1000, ymax = upr/1000), alpha = .25, fill = "steelblue") +
  geom_line() +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Climate penalty on stored energy (ΔEAR) — Parametric simulation",
       y = "ΔEAR (GW·month)", x = "Year") +
  theme_minimal()

Penalidade climática agregada e intervalos de confiança

m1_summary <- m1_sim_agg %>%
  group_by(nom_subsistema.x) %>%
  summarise(
    mean_penalty = mean(mean, na.rm = TRUE),
    lwr_95 = mean(lwr, na.rm = TRUE),
    upr_95 = mean(upr, na.rm = TRUE),
    sd_penalty = sd(mean, na.rm = TRUE),
    .groups = "drop"
  )

# Convertendo para GWh/ano para interpretação
m1_summary <- m1_summary %>%
  mutate(
    mean_GWh = mean_penalty / 1000,
    lwr_GWh  = lwr_95 / 1000,
    upr_GWh  = upr_95 / 1000,
  ) #%>% select(-mean_penalty, -lwr_95, -upr_95)

# Tabela formatada (publicação)
library(gt)
m1_summary[,c("nom_subsistema.x","mean_GWh","lwr_GWh", "upr_GWh")] %>% janitor::adorn_totals() %>% 
  gt() %>%
  fmt_number(columns = 2:4, decimals = 1) %>%
  cols_label(
    nom_subsistema.x = "Subsystem",
    mean_GWh = "Mean ΔEAR (GWh)",
    lwr_GWh  = "Lower 95%",
    upr_GWh  = "Upper 95%",
  ) %>%
  tab_header(
    title = md("**Aggregate climate penalty on stored energy (2001–2018)**"),
    subtitle = "Simulated mean and 95% confidence intervals per subsystem"
  )

Aggregate climate penalty on stored energy (2001–2018)
Simulated mean and 95% confidence intervals per subsystem
Subsystem	Mean ΔEAR (GWh)	Lower 95%	Upper 95%
NORDESTE	−188.4	−200.5	−176.1
NORTE	−22.9	−32.5	−13.4
SUDESTE	−1,178.4	−1,262.4	−1,097.4
SUL	−216.5	−250.7	−183.3
Total	−1,606.3	−1,746.2	−1,470.2

A Tabela apresenta a penalidade climática média sobre a energia armazenada (ΔEAR) no período 2001–2018, agregada por subsistema. Os valores foram estimados a partir das distribuições obtidas via simulação paramétrica de 2.000 sorteios da matriz de covariância dos coeficientes do modelo hidrológico. As penalidades médias foram expressas em gigawatts-mês (GWh·month), com intervalos de confiança de 95% com um total de **–1,606.3 (IC95%: −1,746.2;−1,470.2) de GWh.

Os resultados indicam que o subsistema Sudeste concentra o maior déficit médio de armazenamento, estimado em aproximadamente −1,178.4 (IC95%:−1,262.4;−1,097.4) de GWh·mês, refletindo sua predominância no sistema interligado e sua maior exposição a anomalias de precipitação. O Sul apresenta penalidade intermediária (XXXXXXXXXXXX), mas com maior variabilidade interanual, compatível com seu regime pluviométrico mais irregular. No Nordeste, as perdas médias são menores (XXXXXXXXXXXXXX), embora com tendência negativa acentuada após 2010. Já o Norte apresenta valores próximos de zero, com intervalos de confiança que incluem a nulidade, sugerindo ausência de penalidade líquida significativa no período analisado.

Esses resultados corroboram a hipótese de que a penalidade climática sobre o armazenamento hídrico é concentrada e estruturalmente assimétrica, afetando com maior intensidade os subsistemas Sudeste e Sul — justamente aqueles com maior papel de regulação e suporte energético para o restante do país. Essa heterogeneidade espacial fornece o elo empírico para a próxima etapa da modelagem, em que as variações simuladas de \(\Delta EAR\) serão propagadas à geração hidrelétrica (GAM 1), à geração térmica (GAM 2) e, por fim, ao preço de curto prazo (GAM 3).

3.2 Modelo 2 — Geração Hidrelétrica como Função do Armazenamento O segundo modelo da cadeia estima a geração hidrelétrica (G_hidro) em função da energia armazenada (EAR) e de controles operacionais e sazonais. O objetivo é capturar a elasticidade da geração em relação ao estoque hídrico, permitindo traduzir as variações simuladas de ΔEAR (penalidade climática) em impactos sobre a geração efetiva. Assim como no modelo climático, emprega-se um GAM para capturar relações não lineares entre variáveis físicas e operacionais, mantendo parcimônia e interpretabilidade.

3.2 Modelo 2 — Geração Hidrelétrica como Função do Armazenamento

m2_hidro <- bam(
  val_geracao ~  Year+
    s(ear_reservatorio_subsistema_proprio_mwmes, k = 4)+
    s(val_geracao_day_subsistema, k = 4)+
    s(doy, bs = "cc", k = 2) +     # sazonalidade do dia do ano (cíclica)
    s(id_reservatorio.x, bs = "re")+# efeito aleatório por reservatório
   #s(val_vazaodefluente,k = 2)+

    nom_subsistema.x
  ,
  data     = reserv,
  method   = "fREML",
  discrete = TRUE,
  family   = gaussian(),            # família: normal; ajustar se necessário
  na.action= na.exclude,
  #knots    = list(doy = c(0.5, 366.5))
)

Aviso em smooth.construct.cc.smooth.spec(object, dk$data, dk$knots) :
  basis dimension, k, increased to minimum possible

summary(m2_hidro)
gam.check(m2_hidro)
concurvity(m2_hidro, full=TRUE)

plot(m2_hidro, pages=3, scheme=1, shade=TRUE, se=TRUE, scales="free")

O modelo hidrelétrico (GAM 2) apresentou forte desempenho explicativo, com R² ajustado de 0,85 e deviance explicada de 85,2%, o que confirma sua alta capacidade de reproduzir a geração observada a partir da energia armazenada e variáveis operacionais. O termo linear para o ano é negativo e altamente significativo (–13,5 MW·mês/ano; p < 0,001), refletindo a tendência estrutural de redução da participação hídrica ao longo do período, possivelmente em função da expansão de fontes térmicas e renováveis não hídricas. Os efeitos fixos por subsistema indicam diferenças modestas, com geração média ligeiramente superior no Norte e inferior no Sudeste, consistentes com o papel estrutural de cada região na oferta de energia e no despacho do sistema.

Os diagnósticos de resíduos evidenciam ajuste adequado, com distribuição aproximadamente normal e ausência de heterocedasticidade sistemática. O histograma dos resíduos é centrado em zero, e o QQ-plot mostra leve curvatura nas caudas, indicando pequena assimetria associada a regimes operacionais extremos (reservatórios cheios ou vazios). O gráfico de resíduos versus preditos exibe duas faixas bem definidas, correspondentes a diferentes regimes de operação hidráulica — um de baixa geração com EAR limitada e outro de operação plena com altos níveis de despacho. Embora esses regimes introduzam leve estratificação, não há indícios de viés direcional nem autocorrelação relevante, o que reforça a validade do modelo gaussiano com link identidade. O teste de dimensionalidade efetiva (k-index > 0.9) e o baixo nível de concorrência (concurvity < 0.3) confirmam a estabilidade numérica dos splines e a ausência de sobreajuste.

Os efeitos parciais ilustram de forma clara a dinâmica física do sistema. O spline da energia armazenada (EAR) apresenta um formato crescente com leve saturação, sugerindo rendimento decrescente: aumentos iniciais de armazenamento produzem grandes acréscimos de geração, mas o efeito marginal diminui à medida que o reservatório se aproxima da capacidade máxima. O spline da geração diária do subsistema mostra relação quase linear positiva, representando o efeito de escala e coordenação entre usinas interligadas. A sazonalidade anual (doy) é fraca, refletindo o fato de que as oscilações sazonais já são internalizadas no nível de armazenamento, enquanto o efeito aleatório por reservatório indica ampla heterogeneidade estrutural — alguns reservatórios operam como usinas de base, com efeito próximo de zero, e outros com forte sensibilidade à variação de EAR. Em conjunto, esses resultados confirmam que o modelo capta adequadamente o mecanismo de conversão hidrológica entre estoque e geração, oferecendo base consistente para a estimativa contrafactual da geração hídrica sob diferentes cenários climáticos.

##Predições factual vs. contrafactual com incerteza (baseline 2001–2005) O objetivo é comparar o armazenamento previsto pelo modelo sob condições climáticas observadas e sob condições neutras, representadas pelo clima médio do período 2001–2005. A diferença entre as duas séries (\(ΔEAR=E^{obs}−E^{cf}\)) representa a penalidade climática líquida, e sua incerteza é estimada via propagação dos erros paramétricos do GAM.

# --- Predições factual e contrafactual ---
reserv$ear_hat_obs <- predict(m1_gam, newdata = nd_obs, type = "response")
reserv$ear_hat_cf  <- predict(m1_gam, newdata = nd_cf,  type = "response")

# Penalidade média
reserv <- reserv %>%
  mutate(delta_ear = ear_hat_obs - ear_hat_cf)

reserv$ger_hat_factual <- predict(
  m2_hidro, newdata = reserv,
  type = "response"
)
reserv$ger_hat_contraf <- predict(
  m2_hidro, newdata = reserv %>%
    mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_hat_cf),
  type = "response"
)
reserv <- reserv %>%
  mutate(delta_Ghidro = ger_hat_factual - ger_hat_contraf)
reserv[,c("delta_Ghidro","ger_hat_factual","ger_hat_contraf","val_geracao")] %>% summary()

(3) Incerteza via Delta Method ———————————————

Vp2 <- vcov(m2_hidro)
X_f <- predict(m2_hidro, newdata = reserv ,
               type = "lpmatrix")
X_c <- predict(m2_hidro, newdata = reserv %>%
                 mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_hat_cf),
               type = "lpmatrix")

X_diff2 <- X_f - X_c
var_delta2 <- rowSums((X_diff2 %*% Vp2) * X_diff2)
se_delta2  <- sqrt(pmax(var_delta2, 0))
z <- qnorm(0.975)

reserv <- reserv %>%
  mutate(
    delta_lwr_ger = delta_Ghidro - z * se_delta2,
    delta_upr_ger = delta_Ghidro + z * se_delta2
  )

# (4) Agregação por ano e subsistema -----------------------------------------
m2_delta_agg <- reserv %>%
  mutate(Year = year(date)) %>%
  group_by(Year, nom_subsistema.x) %>%
  summarise(
    mean_penalty = mean(delta_Ghidro, na.rm = TRUE),
    lwr = mean(delta_lwr_ger, na.rm = TRUE),
    upr = mean(delta_upr_ger, na.rm = TRUE),
    .groups = "drop"
  )
m2_delta_agg %>% summary()
  

Predições do modelo Geração

newdata_m1 %>% 
  mutate(ger_hat = predict(m2_hidro, newdata = newdata_m1, type="response",na.action = na.pass)) -> newdata_m1

# 1) Série temporal: EAR observado vs previsto
newdata_m1 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(mes,nom_subsistema.x) %>%summarise(val_geracao =mean(val_geracao   ,na.rm=T),
                             ger_hat=mean(ger_hat,na.rm=T)) %>%
  ggplot(aes(x = mes,)) +
  geom_line(aes(y = val_geracao, color = "Observado"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ger_hat, color = "Previsto"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "EAR observado vs previsto (GAM)",
       y = "EAR (MWmês)", x = "Data") +
  scale_color_manual(values = c("Observado" = "black", "Previsto" = "red")) +
  theme_minimal()

newdata_m1 %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>%summarise(resid=mean(val_geracao,na.rm=T),
                                              ear_hat=mean(ger_hat,na.rm=T)) %>%
  
  ggplot( aes(x = mes,)) +
  #geom_point(aes(color=resid_col))+
  facet_wrap(~nom_subsistema.x,scales = "free")+
  geom_line(aes(y = resid),col="red") +  geom_line(aes(y=ear_hat)) +
  geom_ribbon(aes(ymin=resid,ymax=ear_hat),alpha=.4,fill = "firebrick",)+
  labs(title = "Resíduos do modelo (observado - previsto)",
       y = "Resíduo", x = "Data") +
  theme_minimal()

Ajuste do GAM: séries temporais e dispersão

A Figura 1 apresenta a comparação entre os níveis de EAR observados e os valores previstos pelo modelo GAM, desagregados por subsistema.

• Nordeste e Sudeste: o modelo acompanha razoavelmente bem as tendências de longo prazo, mas subestima os períodos de seca severa (pós-2012) e superestima alguns picos de recuperação.

• Norte: a previsão tende a suavizar a alta variabilidade interanual, captando melhor a média de longo prazo do que os extremos.

• Sul: há maior aderência entre as curvas, embora os choques hidrológicos mais intensos sejam suavizados pelo modelo.

De forma geral, o modelo é eficaz em capturar o nível médio e a sazonalidade estrutural, mas tende a subestimar eventos críticos de cheia ou de seca.

# 2) Scatterplot Observado vs Previsto
newdata_m1 %>% 
   group_by(date,nom_subsistema.x) %>%summarise(val_geracao=sum(val_geracao,na.rm=T),
                             ger_hat=sum(ger_hat,na.rm=T)) %>%
  ggplot( aes(x = val_geracao, y = ger_hat)) +
  geom_point(alpha = 0.1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "blue")  +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Dispersão EAR observado vs previsto",
       x = "Previsto (GAM)", y = "Observado") +
  theme_minimal()

A Figura 2 mostra os diagramas de dispersão entre os valores observados e previstos de EAR, também por subsistema, com a linha tracejada representando o ajuste perfeito (45°).

• A maior parte dos pontos se concentra próxima da linha de referência, confirmando a boa capacidade explicativa do modelo.

• O Sudeste concentra maior dispersão, coerente com a sua relevância e maior heterogeneidade hidrológica.

• O Nordeste mostra bom alinhamento, mas com alguns desvios em baixos níveis de EAR.

• O Sul apresenta a maior aderência, com nuvem de pontos fortemente colada à linha 45°.

Esses resultados corroboram que o GAM é adequado para capturar a dinâmica média e estrutural do sistema, mas que os valores estimados devem ser interpretados com cautela em condições extremas de escassez ou excesso hídrico.

Contrafactural Geração Termica gam

newdata_m1 %>% 
  mutate(ear_reservatorio_subsistema_proprio_mwmes=ear_reservatorio_subsistema_proprio_mwmes-clim_in$penalidade_ear)->newdata_m2
#newdata_m1 %>% 
#  mutate(ear_reservatorio_subsistema_proprio_mwmes=clim_in$clim2026_hat)->newdata_m2
newdata_m2$ger_hat_cf <- predict(m2_hidro, newdata = newdata_m2, type="response",na.action = na.pass)

# 1) Série temporal: EAR observado vs previsto
newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(,nom_subsistema.x) %>%summarise(ger_hat_cf =mean(ger_hat_cf   ,na.rm=T),
                             ger_hat=mean(ger_hat,na.rm=T)) %>%
  ggplot(aes(x = mes,)) +
  geom_line(aes(y = ger_hat_cf, color = "Previsto - Real"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ger_hat, color = "Previsto - Contrafactual"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Geração Predito vs previsto (Contra-factual)",
       y = "EAR (MWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto - Real" = "black", "Previsto - Contrafactual" = "red")) +
  theme_minimal()
newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(mes,) %>%summarise(ger_hat_cf =mean(ger_hat_cf   ,na.rm=T),
                             ger_hat=mean(ger_hat,na.rm=T)) %>%
  ggplot(aes(x = mes,)) +
  geom_line(aes(y = ger_hat_cf, color = "Previsto - Real"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ger_hat, color = "Previsto - Contrafactual"), size = 1, linetype = "dashed") +#facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Geração Predito vs previsto (Contra-factual)",
       y = "EAR (MWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto - Real" = "black", "Previsto - Contrafactual" = "red")) +
  theme_minimal()

A figura compara, por subsistema, a geração hidrelétrica prevista com clima observado e a geração contrafactual simulada ao substituir o armazenamento (EAR) pelo baseline climático (2001–2010).

• Nordeste e Sudeste: a série com clima observado apresenta quedas acentuadas após 2012, enquanto o contrafactual se mantém em patamares menores e mais estáveis. A divergência crescente reflete a penalidade climática recente sobre a geração.

• Norte: o contrafactual acompanha de perto a trajetória média; a série real mostra oscilações adicionais em anos específicos, sugerindo choques não reproduzidos pelo clima médio.

• Sul: há forte variabilidade interanual nas duas curvas; mesmo assim, a série real tende a ficar abaixo do contrafactual em anos secos, evidenciando substituição por outras fontes nesses períodos.

Em conjunto, os resultados indicam que o clima recente reduziu a geração hídrica em relação ao que seria esperado sob um regime climático médio, com intensidade maior no Sudeste e Nordeste.

newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(Year) %>%summarise(ena_mw7=sum(ger_hat,na.rm=T),
                             clim2026_hat=sum(ger_hat_cf,na.rm=T)) %>%
  ggplot(aes(x = Year,)) +
  geom_line(aes(y = ena_mw7/1000, color = "Previsto"), size = 0.7,alpha=.6) +
  geom_line(aes(y = clim2026_hat/1000, color = "Contrafactual 2010"), size = 1, linetype = "dashed") +#facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Contrafactual baseline (clima médio 2001-2010)",
       y = "Geração (GWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto" = "black", "Contrafactual 2010" = "blue")) +
  theme_minimal()
newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(Year,nom_subsistema.x) %>%summarise(ena_mw7=sum(ger_hat,na.rm=T),
                             clim2026_hat=sum(ger_hat_cf,na.rm=T)) %>%
  ggplot(aes(x = Year,)) +
  geom_line(aes(y = ena_mw7/1000, color = "Previsto"), size = 0.7,alpha=.6) +
  geom_line(aes(y = clim2026_hat/1000, color = "Contrafactual 2010"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Contrafactual baseline (clima médio 2001-2010)",
       y = "Geração (GWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto" = "black", "Contrafactual 2010" = "blue")) +
  theme_minimal()

A Figura 1 compara as trajetórias previstas de geração hidrelétrica sob duas condições:

• Previsto – Real (linha preta): resultados do GAM com clima observado no período;

  Previsto – Contrafactual (linha vermelha tracejada): resultados do mesmo modelo, mas substituindo as variáveis climáticas pelas médias históricas (baseline 2001–2010).

newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(Year,nom_subsistema.x) %>%summarise(ena_mw7=sum(ger_hat,na.rm=T),
                             clim2026_hat=sum(ger_hat_cf,na.rm=T)) %>%
  ggplot(aes(x = Year,)) +
  geom_line(aes(y = ena_mw7/1000, color = "Previsto"), size = 0.7,alpha=.6) +
  geom_line(aes(y = clim2026_hat/1000, color = "Contrafactual 2010 - baseline"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free",ncol = 2)+
  labs(title = "Contrafactual climático por subsistema (baseline 2001–2010)",
       y = "Geração (GWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto" = "black", "Contrafactual 2010 - baseline" = "blue")) +
  theme_minimal()

newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade_geração_hidro=ger_hat-ger_hat_cf,
          #penalidade2=ena_mw7-clim2026_hat,
          )->newdata_m2
newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade1=ger_hat-ger_hat_cf,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(Year) %>%summarise(ena_mw7=sum(penalidade1,na.rm=T),
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>%
  ggplot(aes(x = Year,)) +
  #geom_line(aes(y = ena_mw7, color = "previsto - Contrafactual"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ena_mw7/1000, color = "Observado - Contrafactual 2026"), size = 1, linetype = "dashed") +#facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Impacto líquido do clima recente sobre a Geração",subtitle = "Observado - Contrafactual 2026",
       y = "Penalidade Climática na Geração Hidroelética (MWmês)", x = "Data") +
    geom_ribbon(aes(ymin=ena_mw7/1000,ymax=0),alpha=.4,fill = "firebrick",)+

  scale_color_manual(values = c("Observado" = "black", "Contrafactual 2026" = "blue")) +
  theme_minimal()

newdata_m2 %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade1=ger_hat-ger_hat_cf,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(Year,nom_subsistema.x) %>%summarise(penalidade_geração_hidro=sum(penalidade_geração_hidro,na.rm=T),
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>%
  ggplot(aes(x = Year,)) +
  #geom_line(aes(y = ena_mw7, color = "previsto - Contrafactual"), size = 0.7,alpha=.6) +
  geom_line(aes(y = penalidade_geração_hidro/1000, color = "Observado - Contrafactual 2026"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Impacto líquido do clima recente sobre a Geração por subsistema",subtitle = "Observado - Contrafactual 2026",
       y = "Penalidade Climática na Geração Hidroelética (MWmês)", x = "Data") +
    geom_ribbon(aes(ymin=penalidade_geração_hidro/1000,ymax=0),alpha=.4,fill = "firebrick",)+

  scale_color_manual(values = c("Observado" = "black", "Contrafactual 2026" = "blue")) +
  theme_minimal()

newdata_m2 %>% 
   mutate(penalidade1=ger_hat-ger_hat_cf,
          penalidade_perc=ger_hat/ger_hat_cf,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(,nom_subsistema.x) %>%summarise(penalidade1=sum(penalidade1,na.rm=T),
                                          ger_hat=sum(ger_hat,na.rm=T),
                                          ger_hat_cf=sum(ger_hat_cf,na.rm=T),
                                          penalidade_perc=mean(penalidade_perc,na.rm=T)
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>% mutate(custo_financeiro_marginal=penalidade1*100)->penalidade_hidrico
penalidade_hidrico

A Tabela resume os resultados quantitativos. As colunas reportam a geração prevista real (ger_hat), o contrafactual (ger_hat_cf), a penalidade absoluta (MWmês) e a penalidade percentual média (%). Os valores indicam que:

• O Sudeste registrou a maior penalidade absoluta (≈ –3,5 milhões de MWmês), com perda média de 4,4%.

• O Nordeste perdeu ≈ –1 milhão de MWmês (–13,5%), confirmando vulnerabilidade estrutural.

• O Sul apresentou penalidade menor em termos relativos (–9,2%), mas com alta variabilidade interanual.

• O Norte foi o único subsistema a manter média contrafactual próxima à real, refletindo equilíbrio mais favorável entre clima e armazenamento.

newdata_m1$ger_hat_cf <- newdata_m2$ger_hat_cf 
newdata_m1$penalidade_geração_hidro <- newdata_m2$penalidade_geração_hidro 

Modelo de Propagação para outras fontes de energia

ger$nom_tipocombustivel %>% unique()
ger %>%mutate(usina_fonte=paste0("ger_",nom_tipousina)) %>% 
  filter(!nom_subsistema=="PARAGUAI") %>% #,"_",nom_tipousina)) %>% 
  mutate(usina_fonte=case_when(usina_fonte=="ger_BOMBEAMENTO"~"ger_HIDROELÉTRICA",
                   TRUE~usina_fonte)) %>% #,"_",nom_tipousina)) %>% 

  group_by(date,nom_subsistema,usina_fonte,val_geracao_day_subsistema) %>% 
  summarise(ger=sum(val_geracao,na.rm = T)) %>% ungroup() %>% 
  spread(usina_fonte,ger,fill = 0) %>% janitor::clean_names() %>% 
  mutate()->gera_energ_day
gera_energ_day$doy    <- lubridate::yday(gera_energ_day$date)
gera_energ_day$dow    <- lubridate::wday(gera_energ_day$date, label=TRUE, abbr=TRUE)
gera_energ_day$month  <- lubridate::month(gera_energ_day$date)
gera_energ_day$trend  <- lubridate::year(gera_energ_day$date)

ger %>%
  mutate(usina_fonte=paste0("ger_",nom_tipousina),
         mes = floor_date(date, unit = "month")) %>%# select(usina_fonte) %>% unique()
  mutate(usina_fonte=case_when(usina_fonte=="ger_BOMBEAMENTO"~"ger_HIDROELÉTRICA",
                   TRUE~usina_fonte)) %>% #,"_",nom_tipousina)) %>% 
  group_by(mes,usina_fonte,nom_subsistema) %>% #nom_subsistema
  summarise(ger=mean(val_geracao,na.rm = T)) %>% 
    ggplot(aes(x = mes,y=ger,col=nom_subsistema)) +
    geom_line() +facet_wrap(~usina_fonte,scales = "free",ncol = 2)+
  labs(title = "Impacto líquido do clima recente sobre a ENA por subsistema",subtitle = "Observado - Contrafactual 2026",
       y = "ger (GWmês)", x = "Data") +
  theme_minimal()

pld_full %>%
                 mutate(
                   ger_hidroeletrica = ifelse(ger_hidroeletrica == 0, 1, ger_hidroeletrica),
                   ger_fotovoltaica = ifelse(ger_fotovoltaica == 0, 1, ger_fotovoltaica),
                   ger_eolieletrica = ifelse(ger_eolieletrica == 0, 1, ger_eolieletrica),
                   ger_termica = ifelse(ger_termica == 0, 1, ger_termica),
                   Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,
                 ) %>% 
  ggplot(aes(date,))+
    geom_line(aes(y=Demanda),col="red")+
  geom_line(aes(y=val_geracao_day_subsistema,alpha = 0.1),col="black")+facet_wrap(~nom_subsistema.x,scales = "free_y")

gera_energ_day %>% str()
options(scipen = 8)
lm(ger_termica~
     ger_hidroeletrica+
     ger_eolieletrica+ger_fotovoltaica+ger_nuclear+    
     nom_subsistema.x+0+
     trend+as.factor(month),
   data =gera_energ_day) %>% summary()

lm(log(ger_termica)~
     log(ger_hidroeletrica)+
     ger_eolieletrica+ger_fotovoltaica+ger_nuclear+    
     nom_subsistema.x+0+
    trend+as.factor(month),
,
   data =gera_energ_day %>% 
     mutate(ger_hidroeletrica=case_when(ger_hidroeletrica<0.01~1,
                                             TRUE~ger_hidroeletrica),
            ger_termica=case_when(ger_termica<0.01~1,
                                             TRUE~ger_termica))
   ) %>% summary()

gera_energ_day$date %>% summary()
lm(ger_termica~
     ger_hidroeletrica+
     ger_eolieletrica+ger_fotovoltaica+ger_nuclear+    
     nom_subsistema.x+0+
     trend+as.factor(month),
,
   data =gera_energ_day %>% filter(date>2018-01-01 )) %>% summary()

lm(ger_termica~
     ger_hidroeletrica+
     ger_eolieletrica+ger_fotovoltaica+ger_nuclear +
     trend+as.factor(month)
     #nom_subsistema+0
     ,
   data =gera_energ_day) %>% summary()

lm(ger_termica~
     ger_hidroeletrica*nom_subsistema.x+
     ger_eolieletrica+ger_fotovoltaica+ger_nuclear +
     trend+as.factor(month),
     #nom_subsistema+0,
   data =gera_energ_day) %>% summary()

library(mgcv)
    # ger_eolieletrica+ger_fotovoltaica+ger_nuclear +
    # trend+as.factor(month)
m_gam_fontes <- bam(ger_termica  ~ s(ger_hidroeletrica, k=3) + s(ger_eolieletrica, k=4) +
                      s(ger_fotovoltaica, k=3)   +             + s(ger_nuclear, k=3) +
                      s(trend, k=3) + s(month, k=4) +
                      #s(Demanda, k=4)+
                      s(nom_subsistema.x, bs = "re"),
             data = gera_energ_day %>% 
               mutate(Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,
))
summary(m_gam_fontes)

gam.check(m_gam_fontes)                   # checa resíduos, QQ-plot, k-adequação

concurvity(m_gam_fontes, full=TRUE)       # multicolinearidade não linear
plot(m_gam_fontes, pages=4, scheme=1, shade=TRUE, se=TRUE,scales = "free")

Contrafactual Outras fontes de energia

#newdata_m1$ear_hat_cf<-clim_in$clim2026_hat

newdata_m1 %>% 
  group_by(nom_subsistema.x,date) %>% 
  summarise(ear=sum(ear_reservatorio_subsistema_proprio_mwmes,na.rm = T),
            val_geracao=sum(val_geracao),
            ear_hat=sum(ear_hat,na.rm = T),
            ger_hat=sum(ger_hat,na.rm = T),
            penalidade_geração_hidro=sum(penalidade_geração_hidro,na.rm = T),
            ger_hat_cf=sum(ger_hat_cf,na.rm = T),
            ear_hat_cf=sum(ear_hat_cf,na.rm = T)) %>% 
  right_join(gera_energ_day,by = c("nom_subsistema.x"="nom_subsistema","date"="date"))->gera_energ_day
gera_energ_day

gera_energ_day %>% 
  group_by(trend) %>% summarise(val_geracao=sum(ger_hat,na.rm = T),
                                ger_hidroeletrica=sum(ger_hidroeletrica,na.rm = T)) %>% 
  ggplot(aes(trend))+
  geom_line(aes(y=val_geracao),col="blue")+
  geom_line(aes(y=ger_hidroeletrica),col="red")+theme_minimal()
predict(m_gam_fontes,newdata = ger %>% 
  mutate(val_geracao=ger_hat) )

gera_energ_day$ger_termica_prevista<-predict(m_gam_fontes,
                                             newdata = gera_energ_day %>%
                                               mutate(
                                                         nom_subsistema=nom_subsistema.x,
                                                      Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear)
                                             )
gera_energ_day$ger_termica_prevista_cf<-predict(m_gam_fontes,
                                                newdata = gera_energ_day %>% 
                                                  mutate(ger_hidroeletrica=ger_hidroeletrica+penalidade_geração_hidro,
                                                         nom_subsistema=nom_subsistema.x,
                                                         Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear)
                                                )

gera_energ_day %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(trend) %>%summarise(ger_termica_prevista=sum(ger_termica_prevista,na.rm=T),
                             ger_termica_prevista_cf=sum(ger_termica_prevista_cf,na.rm=T)) %>%
  ggplot(aes(x = trend,)) +
  geom_line(aes(y = ger_termica_prevista/1000, color = "Previsto"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ger_termica_prevista_cf/1000, color = "Contrafactual 2010"), size = 1, linetype = "dashed") +#facet_wrap(~nom_subsistema.x,scales = "free")+
  geom_ribbon(aes(ymin=ger_termica_prevista/1000,ymax=ger_termica_prevista_cf/1000),alpha=.4,fill = "firebrick",)+

  labs(title = "Contrafactual baseline (clima médio 2001-2010)",
       y = "Geração Termica (GWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto" = "black", "Contrafactual 2010" = "blue")) +
  theme_minimal()


gera_energ_day %>% 
   mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(trend,nom_subsistema.x) %>%summarise(ena_mw7=sum(ger_termica_prevista,na.rm=T),
                             clim2026_hat=sum(ger_termica_prevista_cf,na.rm=T)) %>%
  ggplot(aes(x = trend,)) +
  geom_line(aes(y = ena_mw7/1000, color = "Previsto"), size = 0.7,alpha=.6) +
  geom_line(aes(y = clim2026_hat/1000, color = "Contrafactual 2010"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Contrafactual baseline (clima médio 2001-2010)",
       y = "Geração Termica (GWmês)", x = "Data") +
  scale_color_manual(values = c("Previsto" = "black", "Contrafactual 2010" = "blue")) +
  theme_minimal()

gera_energ_day %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade_termico=ger_termica_prevista_cf-ger_termica_prevista,
          #penalidade2=ena_mw7-clim2026_hat,
          )->gera_energ_day
gera_energ_day %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade1=ger_termica_prevista_cf-ger_termica_prevista,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(trend) %>%summarise(ena_mw7=sum(penalidade1,na.rm=T),
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>%
  ggplot(aes(x = trend,)) +
  #geom_line(aes(y = ena_mw7, color = "previsto - Contrafactual"), size = 0.7,alpha=.6) +
  geom_line(aes(y = ena_mw7/1000, color = "Observado - Contrafactual 2026"), size = 1, linetype = "dashed") +#facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Impacto líquido do clima recente sobre a Geração",subtitle = "Observado - Contrafactual 2026",
       y = "Penalidade Climática na Termica (MWmês)", x = "Data") +
    geom_ribbon(aes(ymin=ena_mw7/1000,ymax=0),alpha=.4,fill = "firebrick",)+

  scale_color_manual(values = c("Observado" = "black", "Contrafactual 2026" = "blue")) +
  theme_minimal()

gera_energ_day %>% 
   mutate(mes = floor_date(date, unit = "month"),
          penalidade_termica=ger_termica_prevista_cf-ger_termica_prevista,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(trend,nom_subsistema.x) %>%summarise(penalidade_termica=sum(penalidade_termica,na.rm=T),
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>%
  ggplot(aes(x = trend,)) +
  #geom_line(aes(y = ena_mw7, color = "previsto - Contrafactual"), size = 0.7,alpha=.6) +
  geom_line(aes(y = penalidade_termica/1000, color = "Observado - Contrafactual 2026"), size = 1, linetype = "dashed") +facet_wrap(~nom_subsistema.x,scales = "free")+
  labs(title = "Impacto líquido do clima recente sobre a Geração",subtitle = "Observado - Contrafactual 2026",
       y = "Penalidade Climática na Termica (MWmês)", x = "Data") +
    geom_ribbon(aes(ymin=penalidade_termica/1000,ymax=0),alpha=.4,fill = "firebrick",)+

  scale_color_manual(values = c("Observado" = "black", "Contrafactual 2026" = "blue")) +
  theme_minimal()

am maior volatilidade. Essa decomposição regional é crucial para avaliar a elasticidade da escassez e o redirecionamento da pressão para outras fontes de geração.

gera_energ_day %>% 
   mutate(penalidade1=ger_termica_prevista_cf-ger_termica_prevista,
          penalidade_perc=ger_termica_prevista_cf/ger_termica_prevista,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(,nom_subsistema.x) %>%summarise(penalidade1=sum(penalidade1,na.rm=T),
                                          ger_hat=sum(ger_hat,na.rm=T),
                                          ger_hat_cf=sum(ger_hat_cf,na.rm=T),
                                          penalidade_perc=mean(penalidade_perc,na.rm=T)
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>% mutate(custo_financeiro_usotermico=penalidade1*300)->penalidade_termico

bind_rows(penalidade_hidrico %>% mutate(estima="Hidro"),penalidade_termico %>% mutate(estima="Térmico")) %>% 
  select(-starts_with("custo")) %>% 
  mutate(custo_Reais=case_when(estima=="Hidro"~penalidade1*-100,
                             TRUE~penalidade1*300),
         custo_USD=custo_Reais/5.2
         ) 

bind_rows(penalidade_hidrico %>% mutate(estima="Hidro"),penalidade_termico %>% mutate(estima="Térmico")) %>% 
  select(-starts_with("custo")) %>% 
  mutate(custo_Reais=case_when(estima=="Hidro"~penalidade1*-100,
                             TRUE~penalidade1*300),
         custo_USD=custo_Reais/5.2
         ) %>% 
  ggplot(aes(x=nom_subsistema.x,y=custo_USD/1000000,fill=estima))+geom_col(position = "dodge")+
  geom_text(aes(y=custo_USD/1000000+2,label = round(custo_USD/1000, 1)), position=position_dodge(width=1), size = 3,check_overlap=F)+coord_flip()+
  labs(title = "Custo econômico estimado da penalidade climática (2010-2023)",
       subtitle = "Custo financeiro aproximado do impacto climático sobre a geração hidrelétrica e térmica",
       y = "Custo (US$ mil)", x = "Subsistema") +theme_minimal()

PLD: Custo para o sistema

pld     <- read.csv("preco_semanal.csv")           # PLD semanal
pld$date <- as.Date(pld$datas)   # idem para a base de PLD
pld$week  <- lubridate::isoweek(pld$datas)
pld$year  <- lubridate::year(pld$datas)

#pld$date<-pld$DATA_INICIO %>% as_date() 
pld %>% spread(Submercado,preco)->pld
pld$date %>% min()->inicio
pld %>% filter(DATA_INICIO==inicio)
datas<-gera_energ_day$date %>% unique
data.frame(datas) %>% filter(datas>inicio-1) %>% arrange(datas) %>% left_join(pld,by = c("datas"="date"))->daas_completas
daas_completas
i=2
for (i in 1:length(daas_completas$ANO)) {#length(daas_completas$ANO)
  print(i)
  if (is.na(daas_completas$ANO[i])) {
    daas_completas$ANO[i]<-daas_completas$ANO[i-1]
    daas_completas$MES[i]<-daas_completas$MES[i-1]
    daas_completas$SEMANA[i]<-daas_completas$SEMANA[i-1]
    daas_completas$DATA_FIM[i]<-daas_completas$DATA_FIM[i-1]
    daas_completas$SUDESTE[i]<-daas_completas$SUDESTE[i-1]
    daas_completas$SUL[i]<-daas_completas$SUL[i-1]
    daas_completas$NORDESTE[i]<-daas_completas$NORDESTE[i-1]
    daas_completas$NORTE[i]<-daas_completas$NORTE[i-1]

    
  }
  
}  
daas_completas %>% gather(key = "nom_subsitema",value = "preco",c("NORDESTE","SUDESTE","SUL","NORTE"))->pld

serie <- daas_completas %>%
  pull(NORDESTE) %>%
  ts(start =  c(2001, 181), frequency = 365)
stl(serie,s.window = "periodic") %>% plot

serie <- daas_completas %>%
  pull(SUL) %>%
  ts(start =  c(2001, 181), frequency = 365)
stl(serie,s.window = "periodic") %>% plot

serie <- daas_completas %>%
  pull(SUDESTE) %>%
  ts(start =  c(2001, 181), frequency = 365)
stl(serie,s.window = "periodic") %>% plot

serie <- daas_completas %>%
  pull(NORTE) %>%
  ts(start =  c(2001, 181), frequency = 365)
stl(serie,s.window = "periodic") %>% plot

pld$Submercado %>% unique()
gera_energ_day$nom_subsistema.x %>% unique

#newdata_m1 %>% write.csv("base valores previstos.csv")
#gera_energ_day %>% write.csv("valores_previstos_day.csv")

gera_energ_day %>% inner_join(pld,by=c("date","nom_subsistema.x"="nom_subsitema"))->pld_full
pld_full %>% summary()

pld_full %>% group_by(date) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(date,preco))+geom_line()
pld_full %>% 
  group_by(trend) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(trend,preco))+geom_line()
pld_full %>% 
  group_by(month) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(month,preco))+geom_line()
pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(mes) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(mes,preco))+geom_line()

pld_full %>% 
  ggplot(aes(date,preco,nom_subsistema.x))+geom_line()+facet_wrap(~nom_subsistema.x)

pld_full %>% 
  group_by(trend,nom_subsistema.x) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(trend,preco))+geom_line()+facet_wrap(~nom_subsistema.x)
pld_full %>% 
  group_by(month,nom_subsistema.x) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(month,preco))+geom_line()+facet_wrap(~nom_subsistema.x)
pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 

  group_by(mes,nom_subsistema.x) %>% summarise(preco=mean(preco,na.rm = T)) %>% 
  ggplot(aes(mes,preco))+geom_line()+facet_wrap(~nom_subsistema.x)

modelos ingenuos

pld_full %>% str()

lm(log(preco)~log(ger_hidroeletrica)+log(ger_eolieletrica)+log(ger_fotovoltaica)+log(ger_termica)+
     log(Demanda)+nom_subsistema.x,
   data=pld_full %>%
     mutate(ger_hidroeletrica=ifelse(ger_hidroeletrica==0,1,ger_hidroeletrica),
            ger_fotovoltaica=ifelse(ger_fotovoltaica==0,1,ger_fotovoltaica),
            ger_eolieletrica=ifelse(ger_eolieletrica==0,1,ger_eolieletrica),
            ger_termica=ifelse(ger_termica==0,1,ger_termica),
            #Demanda=ifelse(val_geracao_day_subsistema==0,1,val_geracao_day_subsistema),
            Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,

            )) %>% summary()
pld_full$ger_hidroeletrica %>% summary()

• Aumento de 1% na geração hídrica reduz o preço em 1,7%.

• +1% em eólica → +0,21% no preço.

• Aumento de 1% na solar reduz preço em 0,40%.

• Elasticidade forte: +1% na demanda → +1,95% no preço.

• R² = 0.346 (~35%)

  GAM Preço PLD

gam_mod <- gam(preco ~ #trend+
                 s(log(ger_hidroeletrica)) + 
                 s(log(trend)) +                
                 s(log(doy))+
                 s(log(ger_eolieletrica)) + 
                 s(log(ger_fotovoltaica)) + 
                 s(log(Demanda)) + 
                 s(log(ger_termica))+
                 #s(nom_subsistema.x, bs = "re"),# efeito aleatório por reservatório
                 nom_subsistema.x,
               data = pld_full %>%
                 mutate(
                   ger_hidroeletrica = ifelse(ger_hidroeletrica == 0, 1, ger_hidroeletrica),
                   ger_fotovoltaica = ifelse(ger_fotovoltaica == 0, 1, ger_fotovoltaica),
                   ger_eolieletrica = ifelse(ger_eolieletrica == 0, 1, ger_eolieletrica),
                   ger_termica = ifelse(ger_termica == 0, 1, ger_termica),
                   #Demanda = ifelse(val_geracao_day_subsistema == 0, 1, val_geracao_day_subsistema)
                   Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,
                 ))

summary(gam_mod)
plot(gam_mod, pages = 3, shade = TRUE)

• R² ajustado = 0.50 → o GAM explicou metade da variância dos preços, contra ~0.35 do modelo linear anterior.

  • Deviance explained = 50% → muito mais aderência.

• Hidrelétrica (s(log(ger_hidro)))

  • O efeito é forte até certo ponto e depois satura

• Eólica

  • Linha relativamente plana, levemente negativa em valores altos.

  • Mostra que a eólica tem efeito fraco ou até neutro em faixas normais de operação, mas pode ajudar a baixar preços em condições extremas.

• Solar

  • Efeito praticamente linear e negativo, mas fraco (edf ≈ 1 → quase linear)

• Demanda

-   Claramente crescente, e mais íngreme em valores altos. Isso mostra que o preço responde de forma não linear:

    -   em demandas baixas/médias, efeito suave,

-   em demandas muito altas, preço dispara

-   Exatamente o comportamento esperado em sistemas com despacho por ordem de mérito.

Previsto PLD

#pld_full$prec_previsto<-predict(gam_mod,)

predict(gam_mod,newdata=pld_full %>%
                 mutate(
                   ger_hidroeletrica = ifelse(ger_hidroeletrica == 0, 1, ger_hidroeletrica),
                   ger_fotovoltaica = ifelse(ger_fotovoltaica == 0, 1, ger_fotovoltaica),
                   ger_eolieletrica = ifelse(ger_eolieletrica == 0, 1, ger_eolieletrica),
                                      ger_termica = ifelse(ger_termica == 0, 1, ger_termica),
                   Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,

                   #Demanda = ifelse(val_geracao_day_subsistema == 0, 1, val_geracao_day_subsistema)
                 ),type="response") -> pld_full$prec_previsto

pld_full %>% 

  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>% summarise(prec_previsto=mean(prec_previsto,na.rm = T),
                               preco=mean(preco,na.rm = T),
                               ) %>% 
  ggplot(aes(mes,prec_previsto))+
  geom_line(size = 1, linetype = "dashed",col="red") +
  geom_line(aes(y=preco),col = "black")+
  labs(y="$",x="Year")+theme_minimal()+
  facet_wrap(~nom_subsistema.x,scales = "free")

pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>%summarise(preco=mean(preco,na.rm=T),
                                              prec_previsto=mean(prec_previsto,na.rm=T)) %>%
  ggplot( aes(x = mes,)) +
  #geom_point(aes(color=resid_col))+
  facet_wrap(~nom_subsistema.x,scales = "free")+
  geom_line(aes(y = prec_previsto),col="red") +  #geom_line(aes(y=exp(prec_previsto))) +
  geom_ribbon(aes(ymin=prec_previsto,ymax=preco),alpha=.4,fill = "firebrick",)+
  labs(title = "Resíduos do modelo (observado - previsto)",
       y = "Resíduo", x = "Data") +

  theme_minimal()

pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>% summarise(prec_previsto=mean(prec_previsto,na.rm = T),
                               preco=mean(preco,na.rm = T),
                               ) %>% 
  ggplot(aes(mes,prec_previsto))+
  geom_line(size = 1, linetype = "dashed",col="red") +
  geom_line(aes(y=preco),col = "black")+
  labs(y="$",x="Year")+theme_minimal()+
  facet_wrap(~nom_subsistema.x,scales = "free")
pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>%summarise(preco=mean(preco,na.rm=T),
                                              prec_previsto=mean(prec_previsto,na.rm=T)) %>%
  ggplot( aes(x = mes,)) +
  #geom_point(aes(color=resid_col))+
  facet_wrap(~nom_subsistema.x,scales = "free")+
  geom_line(aes(y = prec_previsto),col="red") +  geom_line(aes(y=prec_previsto)) +
  geom_ribbon(aes(ymin=prec_previsto,ymax=preco),alpha=.4,fill = "firebrick",)+
  labs(title = "Resíduos do modelo (observado - previsto)",
       y = "Resíduo", x = "Data") +
  theme_minimal()

Contrafactual PLD

gera_energ_day$ger_termica_prevista<-predict(m_gam_fontes,
                                             newdata = gera_energ_day %>%
                                               mutate(
                                                         nom_subsistema=nom_subsistema.x,
                                                      Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear)
                                             )
gera_energ_day$ger_termica_prevista_cf<-predict(m_gam_fontes,
                                                newdata = gera_energ_day %>% 
                                                  mutate(ger_hidroeletrica=ger_hidroeletrica+penalidade_geração_hidro,
                                                         nom_subsistema=nom_subsistema.x,
                                                         #ger_termica=ger_termica-penalidade_termico,
                                                         Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear)
                                                )

predict(gam_mod,newdata=pld_full %>%
                 mutate(
                  Demanda = ger_eolieletrica+ger_fotovoltaica+ger_hidroeletrica+ger_termica+ger_nuclear,
                   ger_hidroeletrica = ifelse(is.na(ger_hidroeletrica+penalidade_geração_hidro),
                                              1,
                                              ger_hidroeletrica+penalidade_geração_hidro),
                   ger_hidroeletrica = ifelse(ger_hidroeletrica<1 , 1, ger_hidroeletrica),
                   ger_fotovoltaica = ifelse(ger_fotovoltaica == 0, 1, ger_fotovoltaica),
                   ger_eolieletrica = ifelse(ger_eolieletrica == 0, 1, ger_eolieletrica),
                  ger_termica = ifelse(ger_termica+penalidade_termico< 0.1, 1, ger_termica+penalidade_termico),

                 ),type="response") -> pld_full$prec_previsto_cf

pld_full %>% ungroup() %>% 
  #mutate(penalidade_termico)
  select(prec_previsto,prec_previsto_cf,preco,ger_termica,penalidade_termico,penalidade_geração_hidro) %>% summary()

                 

pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>% summarise(prec_previsto=mean(prec_previsto,na.rm = T),
                               preco=mean(preco,na.rm = T),
                               prec_previsto_cf=mean(prec_previsto_cf,na.rm = T),
                               ) %>% 
  ggplot(aes(mes,(prec_previsto)))+
  geom_line(size = 1, linetype = "dashed",col="red") +
  #geom_line(aes(y=preco),col = "black")+
    geom_line(aes(y=(prec_previsto_cf)),col = "blue")+

  labs(y="$",x="Year")+theme_minimal()+
  facet_wrap(~nom_subsistema.x,scales = "free")

pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(mes,nom_subsistema.x) %>%summarise(prec_previsto_cf=mean(prec_previsto_cf,na.rm=T),
                                              prec_previsto=mean(prec_previsto,na.rm=T)) %>%
  ggplot( aes(x = mes,)) +
  #geom_point(aes(color=resid_col))+
  facet_wrap(~nom_subsistema.x,scales = "free")+
  #geom_line(aes(y = exp(prec_previsto_cf)),col="red") + 
  geom_line(aes(y=(prec_previsto_cf)-(prec_previsto))) +
  geom_ribbon(aes(ymin=0,ymax=(prec_previsto_cf)-(prec_previsto)),alpha=.4,fill = "firebrick",)+
  labs(title = "Penalidade (Contrafactual- Predito)",
       y = "Penalidade Climática no Preço", x = "Data") +
  theme_minimal()

pld_full %>% 
  mutate(mes = floor_date(date, unit = "month")) %>% 
  group_by(year,nom_subsistema.x) %>%summarise(prec_previsto_cf=mean(prec_previsto_cf,na.rm=T),
                                              prec_previsto=mean(prec_previsto,na.rm=T)) %>%
  ggplot( aes(x = year,)) +
  #geom_point(aes(color=resid_col))+
  facet_wrap(~nom_subsistema.x,scales = "free")+
  #geom_line(aes(y = exp(prec_previsto_cf)),col="red") + 
  geom_line(aes(y=(prec_previsto_cf)-(prec_previsto))) +
  geom_ribbon(aes(ymin=0,ymax=(prec_previsto_cf)-(prec_previsto)),alpha=.4,fill = "firebrick",)+
  labs(title = "Penalidade (Contrafactual- Predito)",
       y = "Penalidade Climática no Preço", x = "Data") +
  theme_minimal()

pld_full %>% 
   mutate(penalidade1=prec_previsto_cf-prec_previsto,
          penalidade_perc=prec_previsto_cf/prec_previsto,
          #penalidade2=ena_mw7-clim2026_hat,
          ) %>% 

  group_by(nom_subsistema.x) %>%summarise(penalidade1=sum(penalidade1,na.rm=T),
                                          prec_previsto=sum(prec_previsto,na.rm=T),
                                          prec_previsto_cf=sum(prec_previsto_cf*1000,na.rm=T),
                                          prec=mean(preco,na.rm = T),
                                          penalidade_perc=mean(penalidade_perc,na.rm=T)
                             #clim2026_hat=sum(penalidade2,na.rm=T)
                             ) %>% mutate(custo_financeiro_pld=penalidade1)->penalidade_pld
penalidade_pld

ddsad

dsa

bind_rows(penalidade_hidrico %>% mutate(estima="Hidro*"),

          penalidade_termico %>% mutate(estima="Térmico**"),
          penalidade_pld %>% mutate(estima="Operacional***")) %>% 
  select(-starts_with("custo")) %>% 
  mutate(custo_Reais=case_when(estima=="Hidro*"~(-penalidade1)*(296.8),
                               estima=="Térmico**"~penalidade1* 478.9,
                             TRUE~penalidade1*10000/7),
         custo_USD=custo_Reais/5.2
         ) %>% 
  ggplot(aes(x=nom_subsistema.x,y=(custo_USD/100000),fill=estima))+geom_col(position = "dodge")+
  geom_text(aes(label = round(custo_USD/100000, 2)), position=position_dodge(width=1), size = 3,check_overlap=F)+coord_flip()+
  labs(title = "Custo econômico estimado da penalidade climática (2001-2019)",
       subtitle = "Custo financeiro aproximado do impacto climático sobre a geração hidrelétrica e térmica",
       caption = "*296,8 R$/MWh
**478,9 R$/MWh
*** Estimado em 10 GWh por semana",
       y = "Custo (US$ milhões)", x = "Subsistema") +theme_minimal()

bind_rows(penalidade_hidrico %>% mutate(estima="Hidro*"),

          penalidade_termico %>% mutate(estima="Térmico**"),
          penalidade_pld %>% mutate(estima="Operacional***")) %>% 
  select(-starts_with("custo")) %>% 
  mutate(custo_Reais=case_when(estima=="Hidro*"~(-penalidade1)*(296.8),
                               estima=="Térmico**"~penalidade1* 478.9,
                             TRUE~penalidade1*10000/7),
         custo_USD=custo_Reais/5.2
         ) %>% 
  group_by(estima) %>% summarise(custo_USD=sum(custo_USD)) %>% 
  ggplot(aes(x=estima,y=(custo_USD/100000),fill=estima))+geom_col(position = "dodge")+
  geom_text(aes(y=custo_USD/200000,label = round(custo_USD/100000, 1)), position=position_dodge(width=1), size = 3,check_overlap=F)+coord_flip()+
  labs(title = "Custo econômico estimado da penalidade climática (2000-2019)",
       subtitle = "Custo financeiro aproximado do impacto climático sobre a geração hidrelétrica e térmica",
       caption = "*296,8 R$/MWh
**478,9 R$/MWh
*** Estimado em 10 GWh por semana",
       y = "Custo (US$ Milhoes)", x = "Subsistema") +theme_minimal()