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)

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

# 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

pld     <- read.csv("pld_diario_completo.csv")             # weekly spot price (PLD)

Aviso em file(file, "rt") :
  não foi possível abrir o arquivo 'pld_diario_completo.csv': No such file or directory
Error in file(file, "rt") : não é possível abrir a conexão

# 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)
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")
range(ger$date, na.rm=TRUE)
range(reserv$date, na.rm=TRUE)
range(pld$date, na.rm=TRUE)

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")
tab_cov

# 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"
  )

#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

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
a %>%   ggplot(aes(val_longitude,val_latitude))+geom_point()+
ggplot(data=georegion) +geom_sf()
a$nom_subsistema.x %>% unique
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)
gtsummary::tbl_regression(m1_gam)
summary(m1_gam)
gam.check(m1_gam)          # resíduos, QQ, k-index
concurvity(m1_gam, full=TRUE)

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)

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)

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
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)
nsim <- 500
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 <- as.data.table(nd_obs)[, .(date, year = year(date), nom_subsistema.x)]

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),
  sd  = apply(agg_mat, 1, sd)

)

m1_sim_agg %>% mutate(sd_up=mean+z*sd, sd_lwr=mean-z*sd)->m1_sim_agg

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_ribbon(aes(ymin = 0, ymax = upr/1000), fill = "red3", alpha = 0.1) +

  geom_line() +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Climate penalty on stored energy (ΔEAR) — Parametric simulation quartile",
       y = "ΔEAR (GW·month)", x = "Year") +
  theme_minimal()

m1_sim_agg %>%
  ggplot(aes(x = year, y = mean/1000)) +
  geom_ribbon(aes(ymin = sd_lwr/1000, ymax = sd_up/1000), alpha = .25, fill = "steelblue") +
  geom_ribbon(aes(ymin = 0, ymax = sd_up/1000), fill = "red3", alpha = 0.1) +

  geom_line() +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Climate penalty on stored energy (ΔEAR) — Parametric simulation sd",
       y = "ΔEAR (GW·month)", x = "Year") +
  theme_minimal()

Error in `geom_ribbon()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error:
! objeto 'sd_lwr' não encontrado
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

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  = mean_penalty-z*sd_penalty / 1000,
    upr_GWh  = mean_penalty+z*sd_penalty / 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"
  )

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))
)

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)

C <- reserv %>%
  mutate(
    delta_lwr_ger = delta_Ghidro - z * se_delta2,
    delta_upr_ger = delta_Ghidro + z * se_delta2
  )
C[,c("delta_Ghidro","delta_lwr_ger","delta_upr_ger","val_geracao")] %>% summary()

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

library(ggplot2)
m2_delta_agg %>%
  ggplot(aes(x = Year, y = mean_penalty/1000)) +
  geom_ribbon(aes(ymin = lwr/1000, ymax = upr/1000), fill = "steelblue", alpha = 0.8) +
  geom_ribbon(aes(ymin = 0, ymax = upr/1000), fill = "red3", alpha = 0.1) +

  geom_line(color = "black") +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Propagated climate penalty on hydropower generation (ΔG_hidro)",
       subtitle = "Baseline climate: mean 2001–2005",
       y = "ΔG_hidro (GW·month)", x = "Year") +
  theme_minimal()

 
m2_summary  <- C %>%
  mutate(Year = year(date)) %>% ungroup() %>% 
  group_by( nom_subsistema.x) %>%
  summarise(
    mean_penalty = sum(delta_Ghidro, na.rm = TRUE),
    lwr = sum(delta_lwr_ger, na.rm = TRUE),
    upr = sum(delta_upr_ger, na.rm = TRUE),
    .groups = "drop"
  )

library(gt)
m2_summary[,c("nom_subsistema.x","mean_penalty","lwr", "upr")] %>% janitor::adorn_totals() %>% 
  gt() %>%
  fmt_number(columns = 2:4, decimals = 1) %>%
  cols_label(
    nom_subsistema.x = "Subsystem",
    mean_penalty = "Sum Δger (MWh)",
    lwr  = "Lower 95%",
    upr  = "Upper 95%",
  ) %>%
  tab_header(
    title = md("**Aggregate climate penalty on hidropower energy (2001–2018)**"),
    subtitle = "Simulated mean and 95% confidence intervals per subsystem"
  )

A Figura apresenta o resultado da propagação da penalidade climática sobre a geração hidrelétrica (ΔG_hidro), obtida pela substituição do clima observado pelo clima médio histórico (2001–2005) no modelo de armazenamento e pela aplicação desses valores contrafactuais no modelo de geração. O gráfico mostra a diferença entre a geração prevista sob clima efetivo e o cenário contrafactual neutro, com intervalos de confiança calculados por duas abordagens: o método analítico (Delta) e a simulação paramétrica.

Os resultados indicam que a penalidade climática média sobre a geração hidrelétrica segue o mesmo padrão espacial observado na energia armazenada, mas com magnitude amplificada devido à resposta elástica do despacho hídrico. O Sudeste apresenta o maior déficit médio, com reduções acumuladas superiores a –4 milhões de MW·mês em anos secos, refletindo a forte dependência da geração ao nível de reservatórios. O Sul mostra penalidades intermediárias (entre –1 e –2 milhões de MW·mês), com alta variabilidade interanual. O Nordeste e o Norte apresentam oscilações menores, mas com tendência negativa acentuada após 2012, compatível com períodos de estiagem prolongada.

Os intervalos de confiança revelam aumento expressivo da incerteza após 2010, quando as condições climáticas se tornaram mais voláteis e os reservatórios passaram a operar próximos de limites críticos. Esse padrão reforça a interpretação de que a penalidade climática propagada à geração representa não apenas a perda de estoque hídrico, mas também a diminuição da flexibilidade operativa do sistema, um dos principais canais de transmissão dos efeitos climáticos para os custos do setor elétrico.

3.2.5 Propagação conjunta de incerteza (ΔEAR → ΔG_hidro) Este bloco realiza uma simulação Monte Carlo encadeada: cada sorteio inclui simultaneamente os coeficientes do modelo climático (m1_gam) e do modelo hidrelétrico (m2_hidro). Assim, as incertezas da previsão do armazenamento (EAR) são transmitidas integralmente para as previsões de geração, produzindo uma distribuição conjunta para \(Δ𝐺_{ℎ𝑖𝑑𝑟𝑜}\).

3.2.6 Resultados — Propagação integral de incerteza (ΔG_hidro conjunto) partindo do m2——————–

library(MASS)
library(dplyr)
library(lubridate)

set.seed(20251011)
nsim <- 500  # pode aumentar depois para 1000, mas teste com 500 primeiro


# --- 1. Covariâncias e sorteios dos coeficientes ----------------------------
Vp1 <- vcov(m1_gam)
Vp2 <- vcov(m2_hidro)
beta_draws_1 <- MASS::mvrnorm(n = nsim, mu = coef(m1_gam), Sigma = Vp1)
beta_draws_2 <- MASS::mvrnorm(n = nsim, mu = coef(m2_hidro), Sigma = Vp2)

# --- 2. Matrizes do modelo 1 (EAR ~ clima) ----------------------------------
X_obs_ear <- predict(m1_gam, newdata = nd_obs, type = "lpmatrix")
X_cf_ear  <- predict(m1_gam, newdata = nd_cf,  type = "lpmatrix")

# --- 3. Propagação Monte Carlo encadeada -----------------------------------
delta_G_joint <- matrix(NA_real_, nrow = nrow(nd_obs), ncol = nsim)

for (s in seq_len(nsim)) {
  # 1. EAR factual e contrafactual simulados
  ear_f <- as.numeric(X_obs_ear %*% beta_draws_1[s, ])
  ear_c <- as.numeric(X_cf_ear  %*% beta_draws_1[s, ])
  
  # 2. Atualizar inputs do modelo de geração
  nd_f <- nd_obs %>%
    mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_f)
  nd_c <- nd_obs %>%
    mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_c)
  
  # 3. Predições completas do modelo de geração (não linear)
  g_f <- predict(m2_hidro, newdata = nd_f, type = "response")
  g_c <- predict(m2_hidro, newdata = nd_c, type = "response")
  
  # 4. Penalidade propagada
  delta_G_joint[, s] <- g_f - g_c
  print(s)
}


write_rds(delta_G_joint, "delta_G_joint.rds")

delta_G_joint<-read_rds("data/delta_G_joint.rds")

Error in readRDS(con, refhook = refhook) : não é possível abrir a conexão

# --- 4. Agregar por subsistema e ano ----------------------------------------
nd_obs$Year
dt <- nd_obs %>% select(Year, nom_subsistema.x)
dt$Year<-nd_obs$Year


library(data.table)
dt <- as.data.table(dt)
groups <- unique(dt[, .(Year, nom_subsistema.x)])

agg_joint <- lapply(seq_len(nrow(groups)), function(g){
  idx <- which(dt$Year == groups$Year[g] & dt$nom_subsistema.x == groups$nom_subsistema.x[g])
  colSums(delta_G_joint[idx, , drop = FALSE], na.rm = TRUE)
})
agg_joint <- do.call(rbind, agg_joint)
colnames(agg_joint) <- paste0("sim_", 1:nsim)
  agg_df <- cbind(groups, agg_joint)
# --- 5. Estatísticas resumo -------------------------------------------------
m2_joint_summary <- agg_df %>%
  tidyr::pivot_longer(cols = starts_with("sim_"), values_to = "delta_G") %>%
  group_by(Year, nom_subsistema.x) %>%
  summarise(
    mean = mean(delta_G, na.rm = TRUE),
    lwr  = quantile(delta_G, 0.025, na.rm = TRUE),
    upr  = quantile(delta_G, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

library(ggplot2)
m2_joint_summary %>%
  ggplot(aes(x = Year, y = mean)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "firebrick", alpha = 0.3) +
  geom_line(color = "black") +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Joint propagated climate penalty on hydropower generation (ΔG_hidro)",
       subtitle = "Uncertainty fully propagated from climate → storage → generation",
       y = "ΔG_hidro (MW·month)", x = "Year") +
  theme_minimal()

library(ggplot2)
m2_joint_summary %>%
  ggplot(aes(x = Year, y = mean/1000)) +
  geom_ribbon(aes(ymin = lwr/1000, ymax = upr/1000), fill = "steelblue", alpha = 1) +
    geom_ribbon(aes(ymin = 0, ymax = upr/1000), fill = "firebrick", alpha = 0.1) +

  geom_line(color = "black") +
  facet_wrap(~nom_subsistema.x, scales = "free_y") +
  labs(title = "Joint propagated climate penalty on hydropower generation (ΔG_hidro)",
       subtitle = "Uncertainty fully propagated from climate → storage → generation",
       y = "ΔG_hidro (GW·month)", x = "Year") +
  theme_minimal()

A Figura apresenta a penalidade climática sobre a geração hidrelétrica (ΔG_hidro) estimada com propagação conjunta de incerteza ao longo de toda a cadeia climato-hidrológica. Diferentemente das abordagens que tratam os modelos de forma independente, esta simulação incorpora simultaneamente a incerteza dos coeficientes do modelo climático (EAR ~ clima) e do modelo hidrelétrico (geração ~ EAR), gerando uma distribuição conjunta de resultados. Cada trajetória Monte Carlo representa uma combinação plausível de parâmetros e condições climáticas, permitindo capturar a variabilidade estrutural e climática do sistema.

Os resultados mostram que a incerteza propagada amplia os intervalos de confiança das estimativas anuais de ΔG_hidro, principalmente após 2010, quando a variabilidade climática e operacional aumenta significativamente. O subsistema Sudeste concentra as maiores perdas médias (até –4,5 × 10⁶ MW·mês) e a maior dispersão interanual, refletindo tanto a magnitude de sua capacidade de armazenamento quanto sua exposição a secas prolongadas. O Sul apresenta penalidades menores, mas com amplitudes de incerteza elevadas, coerentes com sua instabilidade hidrológica. Nordeste e Norte mostram menores déficits absolutos, mas com tendência negativa contínua.
Esses resultados representam a estimativa mais completa de impacto climático propagado no sistema elétrico, pois integram incertezas climáticas, hidrológicas e operacionais dentro de uma mesma estrutura probabilística.

m2_summary <- m2_joint_summary %>%
  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
m2_summary <- m2_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)
m2_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 ΔG_hidro (GWh)",
    lwr_GWh  = "Lower 95%",
    upr_GWh  = "Upper 95%",
  ) %>%
  tab_header(
    title = md("**Aggregate climate penalty on Hidric Generation (2001–2018)**"),
    subtitle = "Simulated mean and 95% confidence intervals per subsystem"
  )

Aggregate climate penalty on Hidric Generation (2001–2018)
Simulated mean and 95% confidence intervals per subsystem
Subsystem	Mean ΔG_hidro (GWh)	Lower 95%	Upper 95%
NORDESTE	−6.8	−7.0	−6.6
NORTE	−1.9	−2.1	−1.6
SUDESTE	−79.0	−82.5	−76.0
SUL	−16.7	−17.8	−15.9
Total	−104.4	−109.4	−100.1

3.3 Modelo 3 — Geração Termica como Função da Geração Hidrelétrica

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)

O Modelo 3 estima a resposta da geração térmica às variações na geração hidrelétrica, incorporando a competição com outras fontes e o ciclo sazonal de demanda. O modelo aditivo generalizado (GAM) foi ajustado para o período 2000–2019, com estrutura semiparamétrica que permite capturar relações não lineares entre as fontes de geração e seus determinantes temporais. Os resultados indicam uma relação negativa e estatisticamente significativa entre geração hidrelétrica e térmica, evidenciando o papel da térmica como mecanismo de compensação durante períodos de escassez hídrica. O efeito marginal de ger_hidroeletrica é fortemente não linear, com respostas mais intensas em faixas intermediárias de redução, o que sugere que o despacho térmico atinge limites estruturais de expansão em anos de seca severa.

library(mgcv)
library(dplyr)

# --- Preparação dos dados ---------------------------------------------------
gera_energ_day <- gera_energ_day %>%
  mutate(
    Demanda = ger_eolieletrica + ger_fotovoltaica + ger_hidroeletrica + ger_termica + ger_nuclear,
    month = as.numeric(format(date, "%m")),
    trend = as.numeric(date - min(date)) / 30  # tendência linear em meses
  )

# --- Modelo GAM 3: Geração Térmica como função da Hidro + outras fontes -----
m3_gam_termica <- bam(
  ger_termica ~ 
    s(ger_hidroeletrica, k = 4) +    # principal efeito substitutivo
    s(ger_eolieletrica, k = 4) +
    s(ger_fotovoltaica, k = 4) +
    s(ger_nuclear, k = 3) +
    s(trend, k = 3) +
    s(month, bs = "cc", k = 4) +     # padrão sazonal
    s(nom_subsistema, bs = "re"),  # efeito aleatório por subsistema
  data = gera_energ_day,
  method = "fREML",
  discrete = TRUE,
  family = gaussian(),
  na.action = na.exclude
)

summary(m3_gam_termica)

summary(m3_gam_termica)
gtsummary::tbl_regression(m3_gam_termica)
summary(m3_gam_termica)
gam.check(m3_gam_termica)          # resíduos, QQ, k-index
concurvity(m3_gam_termica, full=TRUE)
plot(m3_gam_termica, pages=4, scheme=1, shade=TRUE, se=TRUE)

\(𝑔𝑒𝑟 𝑡𝑒𝑟𝑚𝑖𝑐𝑎_𝑡=𝑓_1(𝑔𝑒_h𝑖𝑑𝑟𝑜_𝑡)+𝑓_2(𝑔𝑒𝑟_𝑒𝑜𝑙𝑖𝑐𝑎_𝑡)+𝑓_3(𝑔𝑒𝑟_𝑠𝑜𝑙𝑎𝑟_𝑡)+_𝑓_4(𝑔𝑒𝑟_𝑛𝑢𝑐𝑙𝑒𝑎𝑟𝑡)+𝑓_5(𝑡𝑟𝑒𝑛𝑑_𝑡)+𝑓_6(𝑚𝑜𝑛𝑡ℎ_𝑡)+𝑢𝑠𝑢𝑏𝑠𝑖𝑠𝑡𝑒𝑚𝑎+𝜀\) Qualidade do ajuste

O modelo apresenta R² ajustado de 0,64 e deviance explicada de 64,1%, indicando boa capacidade preditiva diante da alta variabilidade interdiária. Todos os efeitos suaves são altamente significativos (p < 0.001), com exceção do intercepto. O histograma e o QQ-plot dos resíduos mostram distribuição aproximadamente normal, com leve assimetria nas caudas — reflexo de picos de despacho térmico em períodos de seca severa.

O gráfico “Resíduos vs. Valores Ajustados” revela dispersão crescente com o aumento da geração prevista, sugerindo heterocedasticidade moderada associada a regimes distintos de operação (base × pico). Ainda assim, os resíduos são centrados e sem padrão sistemático, confirmando a adequação geral da especificação.

Efeitos parciais

Os efeitos marginais estimados exibem comportamentos coerentes com a lógica física do sistema:

• s(ger_hidroeletrica) → forte relação inversa e não linear: reduções na geração hídrica aumentam a geração térmica, especialmente em faixas intermediárias de operação (substituição estrutural). Em níveis extremos, o efeito se estabiliza, indicando limite de despacho térmico.

-s(ger_eolieletrica) e s(ger_fotovoltaica) → ambas apresentam efeitos substitutivos negativos, evidenciando o papel das renováveis na mitigação do uso térmico.

• s(ger_nuclear) → relação levemente côncava, estável, refletindo despacho praticamente fixo.

• s(trend) → tendência ascendente, indicando aumento gradual da base térmica no longo prazo, coerente com a expansão de termelétricas entre 2005–2015.

• s(month) → padrão quase plano, com fraca sazonalidade residual após o controle hidrológico.

• Efeito aleatório (nom_subsistema) → o ajuste confirma diferenças estruturais regionais: subsistemas do Sudeste e Sul exibem maior intensidade térmica relativa que Norte e Nordeste.

Síntese O modelo confirma o papel contracíclico da geração térmica no sistema elétrico brasileiro: quando o clima reduz a geração hídrica, a resposta térmica compensa parte do déficit, com sensibilidade não linear e limites operacionais definidos.

A estrutura suavizada dos efeitos permite quantificar a elasticidade marginal entre fontes, que será usada na próxima etapa para estimar a penalidade térmica propagada (ΔG_term) resultante da penalidade hídrica (ΔG_hidro) obtida no Modelo

Predições factual vs. contrafactual (propagação de incerteza)

# --- 5. Loop principal de propagação Monte Carlo -----------------------------
for (s in seq_len(nsim)) {
  # ===== (1) EAR factual e contrafactual (M1) no nível de reservatório×dia
  ear_f <- as.numeric(X1_obs %*% b1[s, ])
  ear_c <- as.numeric(X1_cf  %*% b1[s, ])

  # ===== (2) Predição de G_hidro factual e cf (M2) no nível de reservatório×dia
  nd2_f <- reserv %>% mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_f)
  nd2_c <- reserv %>% mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_c)

  X2_f <- predict(m2_hidro, newdata = nd2_f, type = "lpmatrix")
  X2_c <- predict(m2_hidro, newdata = nd2_c, type = "lpmatrix")

  g_h_f_res <- as.numeric(X2_f %*% b2[s, ])  # por reservatório×dia
  g_h_c_res <- as.numeric(X2_c %*% b2[s, ])

  # ===== (3) Agregar G_hidro por subsistema×dia (para alimentar M3)
  df_h <- data.frame(date = keys_m2$date,
                     subs = keys_m2$nom_subsistema.x,
                     g_h_f = g_h_f_res,
                     g_h_c = g_h_c_res)

  g_h_agg <- df_h |>
    group_by(date, subs) |>
    summarise(g_h_f = sum(g_h_f, na.rm = TRUE),
              g_h_c = sum(g_h_c, na.rm = TRUE), .groups = "drop")

  # ΔG_hidro_clima (subsistema×dia)
  g_h_agg <- g_h_agg |>
    mutate(delta_g_hidro_clima = g_h_f - g_h_c)

  # ===== (4) Construir input do M3
  # factual: usa a hidrelétrica observada (gera_energ_day$ger_hidroeletrica)
  # contrafactual: obs − ΔG_hidro_clima
  nd3 <- gera_energ_day |>
    left_join(g_h_agg,
              by = c("date" = "date", "nom_subsistema" = "subs"))

  # se faltar subsistema em algum dia, trata NA como 0 delta
  nd3$delta_g_hidro_clima[is.na(nd3$delta_g_hidro_clima)] <- 0

  g_h_obs <- nd3$ger_hidroeletrica
  g_h_cfM3 <- g_h_obs - nd3$delta_g_hidro_clima# clipping em 0

  nd3_f <- nd3 %>% mutate(ger_hidroeletrica = g_h_obs)   # factual (observado)
  nd3_c <- nd3 %>% mutate(ger_hidroeletrica = g_h_cfM3)  # contrafactual (obs − Δclima)

  # ===== (5) Predição térmica factual e cf (M3) no nível subsistema×dia
  X3_f <- predict(m3_gam_termica, newdata = nd3_f, type = "lpmatrix")
  X3_c <- predict(m3_gam_termica, newdata = nd3_c, type = "lpmatrix")

  g_t_f <- as.numeric(X3_f %*% b3[s, ])
  g_t_c <- as.numeric(X3_c %*% b3[s, ])

  delta_Gterm_mat[, s] <- g_t_f - g_t_c
  print(s)
}

Aggregate climate penalty on stored energy (2001–2018)
Simulated mean and 95% confidence intervals per subsystem
Subsystem	Sum Δger Term (MWh)	Lower 95%	Upper 95%
NORDESTE	−35,567.9	−37,124.3	−34,109.4
NORTE	7,159.6	−172.2	11,126.5
SUDESTE	195,864.3	171,174.0	213,723.0
SUL	567,160.8	539,040.6	587,110.0
Total	734,616.8	672,918.1	777,850.2

#3.4 Modelo 3 PLD: Custo para o sistema

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

#pld$date<-pld$DATA_INICIO %>% as_date() 
pld %>% spread(Submercado,preco)->pld

Error in `spread()`:
Caused by error:
! objeto 'Submercado' não encontrado
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

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

NA
NA

pld$Submercado %>% unique()

NULL

gera_energ_day$nom_subsistema.x %>% unique

Aviso: Unknown or uninitialised column: `nom_subsistema.x`.

NULL

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

pld_full %>% summary()

      date            nom_subsistema     val_geracao_day_subsistema ger_eolieletrica  ger_fotovoltaica  ger_hidroeletrica
 Min.   :2001-06-30   Length:25536       Min.   :    0              Min.   :   0.00   Min.   :  0.000   Min.   :    0    
 1st Qu.:2005-11-12   Class :character   1st Qu.: 5109              1st Qu.:   0.00   1st Qu.:  0.000   1st Qu.: 3736    
 Median :2010-03-27   Mode  :character   Median : 7269              Median :   0.00   Median :  0.000   Median : 6233    
 Mean   :2010-03-29                      Mean   :10564              Mean   : 257.88   Mean   :  5.367   Mean   :10019    
 3rd Qu.:2014-08-11                      3rd Qu.:11367              3rd Qu.:  69.93   3rd Qu.:  0.000   3rd Qu.:13023    
 Max.   :2018-12-31                      Max.   :40055              Max.   :7863.72   Max.   :375.152   Max.   :35712    
                                         NA's   :2078                                                                    
  ger_nuclear      ger_termica          doy         dow           month            trend           Demanda           ANO      
 Min.   :   0.0   Min.   :   0.0   Min.   :  1.0   dom:3650   Min.   : 1.000   Min.   : 18.20   Min.   :    0   Min.   :2001  
 1st Qu.:   0.0   1st Qu.: 235.0   1st Qu.: 94.0   seg:3649   1st Qu.: 4.000   1st Qu.: 71.42   1st Qu.: 5308   1st Qu.:2005  
 Median :   0.0   Median : 974.5   Median :188.0   ter:3646   Median : 7.000   Median :124.62   Median : 7653   Median :2010  
 Mean   : 403.9   Mean   :1428.6   Mean   :185.8   qua:3646   Mean   : 6.612   Mean   :124.67   Mean   :12115   Mean   :2010  
 3rd Qu.:   0.0   3rd Qu.:1933.6   3rd Qu.:277.0   qui:3646   3rd Qu.:10.000   3rd Qu.:177.87   3rd Qu.:14941   3rd Qu.:2014  
 Max.   :2026.2   Max.   :8117.4   Max.   :366.0   sex:3649   Max.   :12.000   Max.   :231.30   Max.   :40688   Max.   :2018  
                                                   sáb:3650                                                                   
      MES             SEMANA      DATA_INICIO          DATA_FIM              week            year           preco       
 Min.   : 1.000   Min.   :  1.0   Length:25536       Length:25536       Min.   : 1      Min.   :2001    Min.   :  4.00  
 1st Qu.: 4.000   1st Qu.:232.0   Class :character   Class :character   1st Qu.:14      1st Qu.:2005    1st Qu.: 18.59  
 Median : 7.000   Median :460.0   Mode  :character   Mode  :character   Median :27      Median :2010    Median : 79.16  
 Mean   : 6.612   Mean   :459.9                                         Mean   :27      Mean   :2010    Mean   :159.29  
 3rd Qu.:10.000   3rd Qu.:688.0                                         3rd Qu.:40      3rd Qu.:2014    3rd Qu.:213.52  
 Max.   :12.000   Max.   :919.0                                         Max.   :53      Max.   :2018    Max.   :822.83  
                                                                        NA's   :21866   NA's   :21866                   

Modelo 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,
                 ))

Error in eval(predvars, data, env) : 
  objeto 'nom_subsistema.x' não encontrado

summary(gam_3_pld)

Family: gaussian 
Link function: identity 

Formula:
preco ~ 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)) + nom_subsistema

Parametric coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            123.845      5.255  23.565   <2e-16 ***
nom_subsistemaNORTE     56.400      6.080   9.276   <2e-16 ***
nom_subsistemaSUDESTE   42.572     17.231   2.471   0.0135 *  
nom_subsistemaSUL       42.793      3.838  11.149   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                            edf Ref.df      F p-value    
s(log(ger_hidroeletrica)) 9.000  9.000 171.44  <2e-16 ***
s(log(trend))             8.985  9.000 680.45  <2e-16 ***
s(log(doy))               8.881  8.996  13.21  <2e-16 ***
s(log(ger_eolieletrica))  8.479  8.922  66.83  <2e-16 ***
s(log(ger_fotovoltaica))  8.627  8.956  60.52  <2e-16 ***
s(log(Demanda))           8.952  8.999 204.85  <2e-16 ***
s(log(ger_termica))       8.963  9.000 158.85  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.614   Deviance explained = 61.5%
GCV =  14797  Scale est. = 14759     n = 25536

summary(gam_3_pld)

Family: gaussian 
Link function: identity 

Formula:
preco ~ 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)) + nom_subsistema

Parametric coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            123.845      5.255  23.565   <2e-16 ***
nom_subsistemaNORTE     56.400      6.080   9.276   <2e-16 ***
nom_subsistemaSUDESTE   42.572     17.231   2.471   0.0135 *  
nom_subsistemaSUL       42.793      3.838  11.149   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                            edf Ref.df      F p-value    
s(log(ger_hidroeletrica)) 9.000  9.000 171.44  <2e-16 ***
s(log(trend))             8.985  9.000 680.45  <2e-16 ***
s(log(doy))               8.881  8.996  13.21  <2e-16 ***
s(log(ger_eolieletrica))  8.479  8.922  66.83  <2e-16 ***
s(log(ger_fotovoltaica))  8.627  8.956  60.52  <2e-16 ***
s(log(Demanda))           8.952  8.999 204.85  <2e-16 ***
s(log(ger_termica))       8.963  9.000 158.85  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.614   Deviance explained = 61.5%
GCV =  14797  Scale est. = 14759     n = 25536

gtsummary::tbl_regression(gam_3_pld)

Characteristic	Beta	95% CI1	p-value
nom_subsistema
NORDESTE	—	—
NORTE	56	44, 68	<0.001
SUDESTE	43	8.8, 76	0.013
SUL	43	35, 50	<0.001
s(log(ger_hidroeletrica))			<0.001
s(log(trend))			<0.001
s(log(doy))			<0.001
s(log(ger_eolieletrica))			<0.001
s(log(ger_fotovoltaica))			<0.001
s(log(Demanda))			<0.001
s(log(ger_termica))			<0.001
1 CI = Confidence Interval

summary(gam_3_pld)

Family: gaussian 
Link function: identity 

Formula:
preco ~ 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)) + nom_subsistema

Parametric coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            123.845      5.255  23.565   <2e-16 ***
nom_subsistemaNORTE     56.400      6.080   9.276   <2e-16 ***
nom_subsistemaSUDESTE   42.572     17.231   2.471   0.0135 *  
nom_subsistemaSUL       42.793      3.838  11.149   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                            edf Ref.df      F p-value    
s(log(ger_hidroeletrica)) 9.000  9.000 171.44  <2e-16 ***
s(log(trend))             8.985  9.000 680.45  <2e-16 ***
s(log(doy))               8.881  8.996  13.21  <2e-16 ***
s(log(ger_eolieletrica))  8.479  8.922  66.83  <2e-16 ***
s(log(ger_fotovoltaica))  8.627  8.956  60.52  <2e-16 ***
s(log(Demanda))           8.952  8.999 204.85  <2e-16 ***
s(log(ger_termica))       8.963  9.000 158.85  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.614   Deviance explained = 61.5%
GCV =  14797  Scale est. = 14759     n = 25536

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

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 23 iterations.
The RMS GCV score gradient at convergence was 0.001482978 .
The Hessian was positive definite.
Model rank =  67 / 67 

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(log(ger_hidroeletrica)) 9.00 9.00    1.00    0.34    
s(log(trend))             9.00 8.99    0.34  <2e-16 ***
s(log(doy))               9.00 8.88    0.99    0.32    
s(log(ger_eolieletrica))  9.00 8.48    0.98    0.07 .  
s(log(ger_fotovoltaica))  9.00 8.63    0.98    0.11    
s(log(Demanda))           9.00 8.95    1.01    0.66    
s(log(ger_termica))       9.00 8.96    0.98    0.12    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

concurvity(gam_3_pld, full=TRUE)

              para s(log(ger_hidroeletrica)) s(log(trend)) s(log(doy)) s(log(ger_eolieletrica)) s(log(ger_fotovoltaica))
worst    0.9792992                 0.9994989     0.8701883  0.16445905                0.9702436                0.7582099
observed 0.9792992                 0.9973993     0.4043606  0.03943132                0.9065618                0.4626105
estimate 0.9792992                 0.9742579     0.6596454  0.03892563                0.8213786                0.4134559
         s(log(Demanda)) s(log(ger_termica))
worst          0.9996432           0.9847597
observed       0.9966872           0.7688939
estimate       0.9727015           0.8624509

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

# --- 2. Matrizes de design do M1 (clima → EAR) -------------------------------
X1_obs <- predict(m1_gam, newdata = nd_obs, type = "lpmatrix")
X1_cf  <- predict(m1_gam, newdata = nd_cf,  type = "lpmatrix")

# --- 3. Preparação de bases e matrizes de saída ------------------------------
delta_PLD_mat <- matrix(NA_real_, nrow = nrow(preco_day), ncol = nsim)

Erro: objeto 'preco_day' não encontrado

# ==========================================================
# --- 4. Loop principal de propagação Monte Carlo ----------
# ==========================================================
for (s in seq_len(nsim)) {

  # ===== (1) EAR factual e contrafactual (M1) por reservatório×dia ===========
  ear_f <- as.numeric(X1_obs %*% b1[s, ])
  ear_c <- as.numeric(X1_cf  %*% b1[s, ])

  # ===== (2) Predição de G_hidro factual e cf (M2) ===========================
  nd2_f <- reserv %>% mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_f)
  nd2_c <- reserv %>% mutate(ear_reservatorio_subsistema_proprio_mwmes = ear_c)

  X2_f <- predict(m2_hidro, newdata = nd2_f, type = "lpmatrix")
  X2_c <- predict(m2_hidro, newdata = nd2_c, type = "lpmatrix")

  g_h_f_res <- as.numeric(X2_f %*% b2[s, ])  # geração hídrica por reservatório×dia
  g_h_c_res <- as.numeric(X2_c %*% b2[s, ])

  # ===== (3) Agregar G_hidro para subsistema×dia =============================
  df_h <- data.frame(date = keys_m2$date,
                     subs = keys_m2$nom_subsistema.x,
                     g_h_f = g_h_f_res,
                     g_h_c = g_h_c_res)

  g_h_agg <- df_h |>
    group_by(date, subs) |>
    summarise(g_h_f = sum(g_h_f, na.rm = TRUE),
              g_h_c = sum(g_h_c, na.rm = TRUE), .groups = "drop") |>
    mutate(delta_g_hidro_clima = g_h_f - g_h_c)

# ===== (4) Construir input do M3 (subsistema×dia) ==========================
  nd3 <- gera_energ_day |>
    left_join(g_h_agg, by = c("date" = "date", "nom_subsistema" = "subs"))

  nd3$delta_g_hidro_clima[is.na(nd3$delta_g_hidro_clima)] <- 0

  g_h_obs <- nd3$ger_hidroeletrica
  g_h_cfM3 <- g_h_obs - nd3$delta_g_hidro_clima   # sem truncagem

  nd3_f <- nd3 %>% mutate(ger_hidroeletrica = g_h_obs)%>%
                 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,
                 )
  nd3_c <- nd3 %>% mutate(ger_hidroeletrica = g_h_cfM3)%>%
                 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,
                 )

  # ===== (5) Predição térmica factual e cf (M3) ==============================
  X3_f <- predict(m3_gam_termica, newdata = nd3_f, type = "lpmatrix")
  X3_c <- predict(m3_gam_termica, newdata = nd3_c, type = "lpmatrix")

  g_t_f <- as.numeric(X3_f %*% b3[s, ])
  g_t_c <- as.numeric(X3_c %*% b3[s, ])

  delta_g_term_clima <- g_t_f - g_t_c

# ===== (6) Preparar input para M4 (preço) ==================================
  nd4 <- pld_full2 |>
    left_join(
      nd3 %>% select(date, nom_subsistema, ger_hidroeletrica) %>% 
        mutate(g_t_f=g_t_f,g_h_cfM3=g_h_cfM3,
               g_t_c=g_t_c,g_h_obs=g_h_obs,
               ),
      by = c("date", "nom_subsistema")
    ) |>
    mutate(ger_termica_f = g_t_f,
           ger_termica_c = g_t_c,
           ger_hidroeletrica_c = g_h_cfM3,
           ger_hidroeletrica_f = g_h_obs)

  # factual
  nd4_f <- nd4 %>% mutate(
    ger_termica = ger_termica_f,
    ger_hidroeletrica = ger_hidroeletrica_f
  )%>%
                 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,
                 )
  # contrafactual (com penalidade climática propagada)
  nd4_c <- nd4 %>% mutate(
    ger_termica = ger_termica_c,
    ger_hidroeletrica = ger_hidroeletrica_c
  )%>%
                 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,
                 )
# ===== (7) Predição de preço factual e cf (M4) =============================
  X4_f <- predict(gam_4_pld, newdata = nd4_f, type = "lpmatrix")
  X4_c <- predict(gam_4_pld, newdata = nd4_c, type = "lpmatrix")

  pld_f <- as.numeric(X4_f %*% b4[s, ])
  pld_c <- as.numeric(X4_c %*% b4[s, ])

  delta_PLD_mat[, s] <- pld_f - pld_c   # penalidade climática em preço (R$/MWh)
  print(s)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50

— 6. Visualização: Fan Chart do impacto no PLD ———-

Aggregate climate penalty on stored energy (2001–2018)
Simulated mean and 95% confidence intervals per subsystem
Subsystem	mean Δpld Term (MWh)	Lower 95%	Upper 95%
NORDESTE	0.7	0.6	0.9
NORTE	1.1	0.9	1.2
SUDESTE	3.0	2.5	3.6
SUL	13.5	12.4	14.7
Total	18.4	16.5	20.3