Plot.Rmd

---
title: "Plot"
author: "Arpita"
date: "2023-07-15"
output: html_document
---

```{r}
library(data.table)
library(lubridate)
library(dplyr)
library(tidyr)
library(readxl)
library(sf)
library(ggplot2)
library(scico)
library(dlmtree)
library(plotrix)
library(latex2exp)
library(viridis)
```

```{r}
countries <- c("CAL", "CAR", "CENT", "FlA", "MIDA", "MIDW", "NE", "NW",
             "SE", "SW", "TEN", "TEX")

if (!file.exists("output")){
    dir.create("output")
}

# ALERT: Run DistributedModelGeneration.Rmd to create all the models before proceeding....

```
  
```{r}
library(dlmtree)
library(readxl)
library(data.table)
library(lubridate)
library(ggplot2)

solar_cap <- fread("data/solar_capacity_factor.csv")
# Change hour to 0-23 for wraparound
solar_cap[, Hour := Hour %% 24]
# Rename Florida region in capacity file
solar_cap[Region == "FLA", Region := "FlA"]


# Policy changes
perc_increase <- c(.05, .1, .15, .2)


# Region to plot (for Fig 2a and 2b)
r <- "CAL"


# DLM at x=0.15 increase for Fig 2a and 2b
x = 0.15 #solar increase %

#Specify year
y <- "2022"


# Get solar capacity and hourly median for year 2022
hr_cap <- solar_cap[Region == r & year == y]
# Load DLNM model
#load(paste0("tdlm_models/mlist_", r, ".Rda"))
#d <- combine.models(dm_mlist)
lags=12
models_filename = paste0("tdlm_models_lag",lags)

load(paste0(models_filename,"/dm_mod_", r, ".Rda"))
d<- dm
# Loop over hours and calculate CO2 offset for increased solar by x%
perc_increase <- c(x)
co2_offset <- rbindlist(lapply(0:23, function(h) {
  m <- hr_cap[Hour == h, solar_median_mwh]
  s <- summary(d, cenval = m,
               pred.at = m * (1 + perc_increase), 
               mcmc = TRUE, verbose = F)
  # Loop over percent increase, record lagged CO2 offset
  rbindlist(lapply(1:length(s$pred.at), function(p) {
    data.table(Region = r,
               Start = h,
               Med_Solar = m,
               Hour = h:(h+lags) %% 24, # Loop back to 0 if Hour > 23
               Wraparound = c(rep(0, sum(h:(h+lags) < 24)),
                              rep(1, sum(h:(h+lags) >= 24))), # Needed for line break in hourly lag plot
               Percent = perc_increase[p],
               Inc_Solar = m*perc_increase[p],
               CO2_offset = -sapply(1:(lags+1), function(l) median(s$dlm_mcmc[l,p,])), #changed to median
               C_lower = -s$ciupper[p,],
               C_upper = -s$cilower[p,])
  }))
}))
```

```{r}
#Figure 2a (for California)
p2a<- ggplot(co2_offset[Percent == x & Hour %in% c(12, 20)]) +
  geom_line(aes(x = Start, y = CO2_offset, 
                color = factor(Hour),
                group = paste0(Hour, Wraparound)), size = 0.5) +
  geom_ribbon(data = co2_offset[Percent == x & Hour %in% c(12, 20)], 
              aes(x = Start, ymin = C_lower, ymax = C_upper, 
              fill = factor(Hour)), alpha=0.3, show.legend = FALSE) +
  geom_point(data = co2_offset[Percent == x & Hour %in% c(12, 20), last(CO2_offset), by = Hour],
             aes(x = Hour, y = V1, color = factor(Hour)), size = 0.5) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "right",
        legend.key.height = unit(1.5, "cm"),
        legend.text = element_text(size = 6),  
        legend.title = element_text(size = 6)) +
  scale_x_continuous(limits=c(0,23), expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
  scale_y_continuous(limits = c(-450, 0.1))+
  scale_color_manual(values = c(2,3),
    labels = c("12:00", "20:00")) +
  labs(x = "Hour of the day", y = TeX("CO$_2$ change (metric tons)"), 
       color = "Observation\nhour",
       title = "")

  png(paste("output/", r, "_hourly_together_median.png", sep=""), width= 3900, height=1800, res=1200)
  print(p2a)
  dev.off() 
```

```{r}
#Figure 2b
co2_offset_ci <- list()
for (r in "CAL") {
  # Get solar capacity and hourly median for year 2022
  hr_cap <- solar_cap[Region == r & year == 2022]
  # Load DLNM model
  #load(paste0("tdlm_models/mlist_", r, ".Rda"))
  #d <- combine.models(dm_mlist)
  load(paste0(models_filename,"/dm_mod_", r, ".Rda"))
  d<- dm
  # Loop over hours and calculate CO2 offset for increased solar by x%
  co2_offset_ci[[r]] <- rbindlist(lapply(0:23, function(h) {
    m <- hr_cap[Hour == h, solar_median_mwh]
    s <- summary(d, cenval = m,
                 pred.at = m * (1 + perc_increase), 
                 mcmc = TRUE, verbose = F)
    # Loop over percent increase, record lagged CO2 offset
    rbindlist(lapply(1:length(s$pred.at), function(p) {
      rbindlist(lapply(1:(lags+1), function(l) {
        data.table(Region = r,
                   Start = h,
                   Hour = (h+l-1) %% 24, # Loop back to 0 if Hour > 23
                   Percent = perc_increase[p],
                   Inc_Solar = m * (1 + perc_increase[p]) - m,
                   Samp = 1:d$mcmcIter,
                   CO2_offset = -s$dlm_mcmc[l, p,])
      }))
    }))
  }))
}


mwh_scale <- round(co2_offset_ci[[r]][Percent == x, mean(Inc_Solar), by = Start][, max(V1)] / 200, 1)
co2_scale <- round(-co2_offset_ci[[r]][Percent == x, sum(CO2_offset), by = .(Hour, Samp)
][, .(lower = quantile(V1, 0.025)), by = .(Hour)][, min(lower)] / 200, 1)

p2b<-  ggplot(co2_offset_ci[[r]][Percent == 0.15, sum(CO2_offset), by = .(Hour, Samp)
][, .(mean = median(V1), lower = quantile(V1, 0.025), #changed to median
      upper = quantile(V1, 0.975)), by = .(Hour)]) +
  geom_line(data = co2_offset_ci[["CAL"]][Percent == x, mean(Inc_Solar), by = Start],
            aes(x = Start, y = V1/mwh_scale),
            size = 0.5, linetype = 2) +
  geom_ribbon(aes(x = Hour, ymin = lower/co2_scale, ymax = upper/co2_scale), fill = "grey") +
  geom_line(aes(x = Hour, y = mean/co2_scale), size = 0.5) +
  geom_hline(yintercept = 0, size = 1) +
  theme_minimal(base_size = 7) +
  theme(legend.position = "top",
        legend.key.width = unit(2, "cm"),
        axis.line = element_line(),
        axis.ticks = element_line()) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
  scale_y_continuous(breaks = c(-200, -100, 0, 100, 200), limits = c(-220, 220),
                     labels = c(-200*co2_scale, -100*co2_scale, 0, 100*mwh_scale, 200*mwh_scale)) +
  scale_color_discrete(labels = paste0(c(5, 10, 15, 20), "%")) +
  labs(x = "Hour of the day", 
       y = TeX("$\\leftarrow$ CO$_2$ change\t Solar change $\\rightarrow$"), #"<- Total CO2 Offset (metric tons) \t|\t Increase in Median Solar (MWh) ->",
       color = "Solar capacity increase", 
       title = "")#paste0(r, ": CO2 offset for alternate policy changes in solar capacity"))


#png(paste(r, "_cumulative.png", sep=""), width=1000, height=800)
png(paste("output/", r, "_cumulative_median.png", sep=""), width= 3900, height=1900, res=1200)
print(p2b)
dev.off()
```

```{r}
#For Figure 3 (hourly total)
perc_increase <- c(0.05, 0.1, 0.15, 0.2)
countries <- c('CAL', 'CAR','FlA', 'MIDA', 'MIDW', 
               'NW', 'SE', 'SW', 'TEX', 'TEN', 'CENT', 'NE')
names <- c("California", "Carolinas", "Florida", "Mid-Atlantic", "Midwest",  "Northwest", "Southeast", "Southwest", "Texas", "Tennessee","Central",   "New England")
names(names) <- countries
co2_offset_comb <- list()
for (r in countries) {
  print(r)
  # Get solar capacity and hourly median for year 2022
  hr_cap <- solar_cap[Region == r & year == y]
  # Load DLNM model
  #load(paste0("tdlm_models/mlist_", r, ".Rda"))
  #d <- combine.models(dm_mlist)
  load(paste0(models_filename, "/dm_mod_", r, ".Rda"))
  d<- dm
  # Loop over hours and calculate CO2 offset for increased solar by x%
  co2_offset_comb[[r]] <- rbindlist(lapply(0:23, function(h) {
    m <- hr_cap[Hour == h, solar_median_mwh]
    s <- summary(d, cenval = m,
                 pred.at = m * (1 + perc_increase), 
                 mcmc = TRUE, verbose = F)
    # Loop over percent increase, record lagged CO2 offset
    rbindlist(lapply(1:length(s$pred.at), function(p) {
      data.table(Region = r,
                 Start = h,
                 Med_Solar = m,
                 Hour = h:(h+lags) %% 24, # Loop back to 0 if Hour > 23
                 Wraparound = c(rep(0, sum(h:(h+lags) < 24)),
                                rep(1, sum(h:(h+lags) >= 24))), # Needed for line break in hourly lag plot
                 Percent = perc_increase[p],
                 Inc_Solar = m * (1 + perc_increase[p]) - m,
                 CO2_offset = -sapply(1:(lags+1), function(l) median(s$dlm_mcmc[l,p,]))) #changed to median
    }))
  }))
}
#Checking upper bound on CO2 emissions
solar_increase = list()
for (r in countries)
{
  solar_row=list()
  for(percent in c(0.05, 0.1, 0.15,0.2))
  {
    solar_row = cbind( solar_row, (sum(unique(co2_offset_comb[[r]][Percent==percent]$Inc_Solar))*365/1000))
  }
  solar_increase = rbind( solar_increase, solar_row)
}
solar_increase = as.data.frame(solar_increase)
rownames(solar_increase) = countries
colnames(solar_increase) = c(0.05, 0.1, 0.15,0.2)


solar_increase[] <- lapply(solar_increase, function(x) as.numeric(as.character(x)))
#           0.05       0.1      0.15       0.2
#CAL  2126.78200 4253.5640 6380.3460 8507.1280
#CAR   486.54500  973.0900 1459.6350 1946.1800
#FlA   561.88100 1123.7620 1685.6430 2247.5240
#MIDA  361.84275  723.6855 1085.5283 1447.3710
#MIDW  223.16100  446.3220  669.4830  892.6440
#NW    731.97100 1463.9420 2195.9130 2927.8840
#SE    409.47525  818.9505 1228.4258 1637.9010
#SW    256.90525  513.8105  770.7157 1027.6210
#TEX  1214.84775 2429.6955 3644.5433 4859.3910
#CENT   24.85650   49.7130   74.5695   99.4260
#TEN    68.92113  137.8423  206.7634  275.6845
#NE     46.46450   92.9290  139.3935  185.8580
column_sums = colSums(solar_increase)

# Print the column sums
print(column_sums)
#     0.05       0.1      0.15       0.2 
# 6513.653 13027.306 19540.959 26054.612 


# Combined plots
co2_plots <- list()
for (r in countries) {
  ylim=-700 #prev:-1000
  if(r %in% c('TEN', 'CENT', 'NE')){
    ylim=-40 #prev:-50
  }
  mwh_scale <- 1#floor(co2_offset_comb[[r]][, first(Inc_Solar), 
                                          # by = .(Start, Percent)][, max(V1)] / 20) / 10
  co2_scale <- 1#floor(-co2_offset_comb[[r]][, .(CO2_offset = sum(CO2_offset)), 
                                 # by = .(Hour, Percent)][, min(CO2_offset)] / 20) / 10
  
  co2_plots[[r]] <- ggplot() +
    # geom_line(data = co2_offset_comb[[r]][, first(Inc_Solar), by = .(Start, Percent)], 
    #           aes(x = Start, y = V1/mwh_scale, color = as.factor(Percent),
    #               linetype = as.factor(Percent)), size = 1) +
    geom_line(data = co2_offset_comb[[r]][, .(CO2_offset = sum(CO2_offset)), 
                                by = .(Hour, Percent)],
              aes(x = Hour, y = CO2_offset/co2_scale, 
                  color = as.factor(Percent)),# linetype = as.factor(Percent)), 
              size = 0.3) +
    geom_hline(yintercept = 0, size = 0.3) +
    theme_minimal(base_size = 10) +
    theme(legend.position = "bottom",
          legend.key.width = unit(0.5, "cm"),
          axis.line = element_line(),
          axis.ticks = element_line(),
          plot.subtitle = element_text(size = 8)) +
    scale_x_continuous(expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                       labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
    scale_y_continuous(limits = c(ylim, 0), expand=c(0,0)) +
    scale_color_discrete(labels = paste0(c(5, 10, 15, 20), "%")) +
    #scale_linetype_discrete(labels = paste0(c(5, 10, 15, 20), "%")) +
    labs(x = "", y = "", subtitle = names[r], 
         color = "Solar capacity increase")#,
         #linetype = "Solar capacity increase")
}

library(ggpubr)

# Arrange the plots into a grid with a common legend
p3 <- ggarrange(
  plotlist = co2_plots,
  nrow = 4,
  ncol = 3,
  common.legend = TRUE,
  legend = "bottom",
  align = "hv"  # Align both horizontally and vertically
)

# Annotate the plot with axis titles
annotated_plot <- annotate_figure(
  p3,
  left = text_grob(
    label = TeX("CO$_2$ change (metric tons)"),
    rot = 90,
    size = 10
  ),
  bottom = text_grob(
    label = "Hour of the day",
    size = 10
  )
)

png(paste("output/cumulative_combined_median.png", sep=""), width= 8000, height=8000, res=1200)
print(annotated_plot)
dev.off()
```

```{r}
#For Figures 4 (per day total), 5 (interregional), 6 (next two decades) and table 1
# Collect all posterior samples for 95% CI calculation
perc_increase <- c(0.05, 0.1, 0.15, 0.2)
countries <- c('CAL', 'CAR','CENT','FlA', 'MIDA', 'MIDW', 'NE',
               'NW', 'SE', 'SW', 'TEN', 'TEX')
names <- c("California", "Carolinas", "Central", "Florida", "Mid-Atlantic", "Midwest","New England", "Northwest", "Southeast", "Southwest",  "Tennessee", "Texas")
names(names) <- countries
co2_offset_comb <- list()
for (r in countries) {
  print(r)
  # Get solar capacity and hourly median for year 2022
  hr_cap <- solar_cap[Region == r & year == y]
  # Load DLNM model
  #load(paste0("tdlm_models/mlist_", r, ".Rda"))
  #d <- combine.models(dm_mlist)
  load(paste0("tdlm_models/dm_mod_", r, ".Rda"))
  d<- dm
  # Loop over hours and calculate CO2 offset for increased solar by x%
  co2_offset_comb[[r]] <- rbindlist(lapply(0:23, function(h) {
    m <- hr_cap[Hour == h, solar_median_mwh]
    s <- summary(d, cenval = m,
                 pred.at = m * (1 + perc_increase), 
                 mcmc = TRUE, verbose = F)
    # Loop over percent increase, record lagged CO2 offset
    rbindlist(lapply(1:length(s$pred.at), function(p) {
      data.table(Region = r,
                 Hour = h,
                 Percent = perc_increase[p],
                 Samp = 1:d$mcmcIter,
                 CO2_offset = -colSums(s$dlm_mcmc[, p,]))
    }))
  }))
}

co2_offset_comb <- rbindlist(co2_offset_comb)
co2_offset_comb$Region <- factor(co2_offset_comb$Region, countries, names)

# For Figure 4

p_daily<-  ggplot(data=co2_offset_comb[, sum(CO2_offset), by = .(Region, Percent, Samp)
][, .(mean = median(V1), lower = quantile(V1, 0.025), #changed to median
      upper = quantile(V1, 0.975)),
  by = .(Region, Percent)],
  aes(x = Percent, y = mean, ymin = lower, ymax = upper, color = factor(Percent))) +
  geom_point(size = 0.5) +
  geom_errorbar(width = 0.025, size = 0.3) +
  facet_wrap(~ Region) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "none", 
        axis.line = element_line(color = "black"),
        axis.ticks = element_line(color = "black")) +
  scale_x_continuous(labels = function(x) paste0(round(x * 100), "%")) +
  labs(x = "Solar increase (in %)",
       y = TeX("Change in CO$_2$ (in metric tons)"))

png(paste("output/daily_median.png", sep=""), width= 5000, height=4000, res=1200)
print(p_daily)
dev.off()
```
```{r}
# For Figure 5(cumulative annual plot for )
countries <- c('CAL', 'CAR','CENT', 'FlA', 'MIDA', 'MIDW', 'NE',
               'NW', 'SE', 'SW', 'TEN', 'TEX')
names <- c("California", "Carolinas", "Central", "Florida", "Mid-Atlantic", "Midwest",  "New England", "Northwest", "Southeast", "Southwest", "Tennessee", "Texas")
# Convert list to dataframe
w=co2_offset_comb[, sum(CO2_offset), by = .(Region, Percent, Samp)
][, .(mean = median(V1), lower = quantile(V1, 0.025), #changed to median
      upper = quantile(V1, 0.975)),
  by = .(Region, Percent)]
w$Percent=w$Percent*100
w$mean = -w$mean*365/1000000

# Stacked plot
p_total<- ggplot(data=w, aes(fill=Region, y=mean, x=Percent)) + 
    geom_bar(position="stack", stat="identity", width=2) +
    #geom_hline(yintercept=30.85) +
    theme_minimal(base_size = 7) +
    labs(x = "Increase in Solar (in percentage)", y = TeX("Change in CO$_2$ (in metric tons)"), color = "Regions")+ coord_flip()+
    theme(legend.key.size = unit(0.3, "cm"))
png(paste("output/total_median.png", sep=""), width= 3900, height=2200, res=1200)
print(p_total)
dev.off()

```

```{r}
# For table S1 in Appendix (total change in CO2 per day)
table1=cat(apply(dcast(co2_offset_comb[, sum(CO2_offset), by = .(Region, Percent, Samp)
][, .(paste0(round(median(V1), 1), " (", round(quantile(V1, 0.025), 1), ", ", #changed to median
             round(quantile(V1, 0.975), 1), ")")),
  by = .(Region, Percent)], Region ~ Percent), 1, paste, collapse = " & "), sep = "\\\\\n")
```


```{r}
#Combine within region and neighboring effect on r1 from each country r2
all_results<- c("r1", "r2", "h", "effect_5", "effect_10", "effect_15", "effect_20")
values<- list()
for(r1 in countries){
  print(r1)
  # Get solar capacity and hourly median for year 2022
  hr_cap <- solar_cap[Region == r1 & year == y]
  # Load DLNM model
  #load(paste0("tdlm_models/mlist_", r1, ".Rda"))
  #d <- combine.models(dm_mlist)
  load(paste0(models_filename, "/dm_mod_", r1, ".Rda"))
  d<- dm
  effect =0
  for(h in 0:23){
          m <- hr_cap[Hour == h, solar_median_mwh]
          s <- summary(d, cenval = m,
                   pred.at = m * (1 + perc_increase), verbose = F, mcmc=TRUE)
          mean <-  apply(s$cumulative_mcmc, 2, median)  #s$cumulative.effect$mean
          effect = effect-mean
          all_results = rbind(all_results, c(r1, r1, h, as.numeric(mean[1]), as.numeric(mean[2]), as.numeric(mean[3]), as.numeric(mean[4])))
         }
  for(r2 in countries){
      #f <- paste0("tdlm_models/mlist_",r2,"-",r1,".Rda")
      f<- (paste0(models_filename, "/dm_mod_", r2, "-", r1, ".Rda"))
      if (file.exists(f)){
        print(f)
        load(f)
        hr_cap <- solar_cap[Region == r2 & year == y]
        for(h in 0:23){
          m <- hr_cap[Hour == h, solar_median_mwh]
          #d <- combine.models(dm_mlist)
          d<- dm
          s <- summary(d, cenval = m,
                   pred.at = m * (1 + perc_increase), verbose = F,mcmc=TRUE)
          mean=apply(s$cumulative_mcmc, 2, median) #s$cumulative.effect$mean changed to median
          effect = effect-mean
          all_results = rbind(all_results, c(r2, r1, h, as.numeric(mean[1]), as.numeric(mean[2]), as.numeric(mean[3]), as.numeric(mean[4])))
         }
      }
  }
  values[[r1]]=effect
}
all_results = as.data.frame(all_results)
colnames(all_results) = all_results[1,]
all_results = all_results[-1,]
write.csv(all_results, "output/all_results_median.csv")


```

```{r}
#For Figure 6 (interregional effect)
wide_df <- as.data.frame(all_results)
colnames(wide_df) <- c("Gen_Region", "Import_Region", "Hour", "Percent_5", "Percent_10", "Percent_15", "Percent_20")
wide_df$Percent_5 = as.numeric(wide_df$Percent_5)
wide_df$Percent_10 = as.numeric(wide_df$Percent_10)
wide_df$Percent_15 = as.numeric(wide_df$Percent_15)
wide_df$Percent_20 = as.numeric(wide_df$Percent_20)
long_df <- pivot_longer(wide_df, cols = starts_with("Percent"), names_to = "Percent", values_to = "CO2_offset")
long_df$CO2_offset = as.numeric(long_df$CO2_offset)
daily_co2 <- long_df %>% group_by(Gen_Region, Import_Region, Percent) %>% summarise(CO2 = sum(CO2_offset, na.rm = TRUE))

# Define the correct levels and labels
levels <- c("CAL", "SW", "NW", "TEX", "CENT", "FlA", "MIDW", "SE", "MIDA", "TEN", "CAR", "NE")
labels <- c("California", "Southwest", "Northwest", "Texas", "Central", "Florida", "Midwest", "Southeast", "Mid-Atlantic", "Tennessee", "Carolinas", "New England")

# Update the factor levels
daily_co2$Gen_Region <- factor(daily_co2$Gen_Region, levels = levels, labels = labels)
daily_co2$Import_Region <- factor(daily_co2$Import_Region, levels = rev(levels), labels = rev(labels))
daily_co2$Percent <- factor(daily_co2$Percent, c("Percent_5", "Percent_10", "Percent_15", "Percent_20"), c("5%", "10%", "15%", "20%"))
daily_co2$CO2 <- -daily_co2$CO2

p_interregional<- ggplot(daily_co2) +
  geom_tile(aes(x = Gen_Region, y = Import_Region, fill = CO2)) +
  facet_wrap(~Percent) +
  scale_fill_viridis_c(limits = c(-10000, 0)) +
  theme_minimal(base_size = 8) +
  theme(axis.line = element_line(color = "black"),
        axis.ticks = element_line(color = "black"),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        legend.key.height = unit(0.5, "cm")) +
  coord_equal() +
  labs(x = "Generating Region", y = "Import Region",
       fill = TeX("Import region CO$_{2}$ change"))
png(paste("output/interregional_median.png", sep=""), width=5200, height=5200, res=1200)
print(p_interregional)
dev.off()

```


```{r}
#Figure in Supplementary 
#Sensitivity Analysis with varying lag hours (for California)
countries <- c('CAL', 'CAR', 'CENT', 'FlA', 'MIDA', 'MIDW', 
               'NE', 'NW', 'SE', 'SW', 'TEN', 'TEX')
solar_cap <- fread("data/solar_capacity_factor.csv")
# Change hour to 0-23 for wraparound
solar_cap[, Hour := Hour %% 24]
# Rename Florida region in capacity file
solar_cap[Region == "FLA", Region := "FlA"]


# Policy changes
perc_increase <- c(.05, .1, .15, .2)


# Region to plot (for Fig 2a and 2b)
r <- "CAL"


# DLM at x=0.15 increase for Fig 2a and 2b
x = 0.15 #solar increase %

#Specify year
y <- "2022"

# Get solar capacity and hourly median for year 2022
hr_cap <- solar_cap[Region == r & year == y]
# Load DLNM model
#load(paste0("tdlm_models/mlist_", r, ".Rda"))
#d <- combine.models(dm_mlist)

lag_vector = c(4,8,12,16)

for (lags in lag_vector){
  models_filename = paste0("tdlm_models_lag",lags)
  load(paste0(models_filename,"/dm_mod_", r, ".Rda"))
  d<- dm
  # Loop over hours and calculate CO2 offset for increased solar by x%
  perc_increase <- c(x)
  co2_offset <- rbindlist(lapply(0:23, function(h) {
    m <- hr_cap[Hour == h, solar_median_mwh]
    s <- summary(d, cenval = m,
                 pred.at = m * (1 + perc_increase), 
                 mcmc = TRUE, verbose = F)
    # Loop over percent increase, record lagged CO2 offset
    rbindlist(lapply(1:length(s$pred.at), function(p) {
      data.table(Region = r,
                 Start = h,
                 Med_Solar = m,
                 Hour = h:(h+lags) %% 24, # Loop back to 0 if Hour > 23
                 Wraparound = c(rep(0, sum(h:(h+lags) < 24)),
                                rep(1, sum(h:(h+lags) >= 24))), # Needed for line break in hourly lag plot
                 Percent = perc_increase[p],
                 Inc_Solar = m*perc_increase[p],
                 CO2_offset = -s$matfit[p,],#-sapply(1:(lags+1), function(l) mean(s$dlm_mcmc[l,p,])), #-s$matfit[p,]
                 C_lower = -s$ciupper[p,],
                 C_upper = -s$cilower[p,])
    }))
  }))
  
  p2a_lags<- ggplot(co2_offset[Percent == x & Hour %in% c(12, 20)]) +
  geom_line(aes(x = Start, y = CO2_offset, 
                color = factor(Hour),
                group = paste0(Hour, Wraparound)), size = 1) +
  geom_ribbon(data = co2_offset[Percent == x & Hour %in% c(12, 20)], 
              aes(x = Start, ymin = C_lower, ymax = C_upper, 
              fill = factor(Hour)), alpha=0.3, show.legend = FALSE) +
  geom_point(data = co2_offset[Percent == x & Hour %in% c(12, 20) & Start==Hour, last(CO2_offset), by = Hour],
             aes(x = Hour, y = V1, color = factor(Hour)), size = 2) +
  theme_minimal(base_size = 24) +
  theme(legend.position = "right",
        legend.key.height = unit(1.5, "cm"),
        legend.text = element_text(size = 18),  
        legend.title = element_text(size = 18)) +
  scale_x_continuous(limits=c(0,23), expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
  scale_y_continuous(limits = c(-700, 0.1))+
  scale_color_manual(values = c(2,3),
    labels = c("12:00", "20:00")) +
  labs(x = "Hour of the day", y = TeX("CO$_2$ change (metric tons)"), 
       color = "Observation\nhour",
       title = paste0(lags, " lag-hours considered"))

  png(paste("output/", r, "_lag", lags, "_hourly_together_median.png", sep=""), width=800, height=400)
  print(p2a_lags)
  dev.off() 
  
  #CO2effect_for_an_hour = co2_offset[Percent == x & Hour %in% c(4, 8, 12,16,20), sum(CO2_offset), by="Hour"]
  #print(paste0("CO2 offset for increasing solar at 4 pm = ", CO2effect_for_an_hour))
  
  #Figure cumulative
  co2_offset_ci <- list()
  for (r in "CAL") {
    # Get solar capacity and hourly median for year 2022
    hr_cap <- solar_cap[Region == r & year == 2022]
    # Load DLNM model
    #load(paste0("tdlm_models/mlist_", r, ".Rda"))
    #d <- combine.models(dm_mlist)
    load(paste0(models_filename,"/dm_mod_", r, ".Rda"))
    d<- dm
    # Loop over hours and calculate CO2 offset for increased solar by x%
    co2_offset_ci[[r]] <- rbindlist(lapply(0:23, function(h) {
      m <- hr_cap[Hour == h, solar_median_mwh]
      s <- summary(d, cenval = m,
                   pred.at = m * (1 + perc_increase), 
                   mcmc = TRUE, verbose = F)
      # Loop over percent increase, record lagged CO2 offset
      rbindlist(lapply(1:length(s$pred.at), function(p) {
        rbindlist(lapply(1:(lags+1), function(l) {
          data.table(Region = r,
                     Start = h,
                     Hour = (h+l-1) %% 24, # Loop back to 0 if Hour > 23
                     Percent = perc_increase[p],
                     Inc_Solar = m * (1 + perc_increase[p]) - m,
                     Samp = 1:d$mcmcIter,
                     CO2_offset = -s$dlm_mcmc[l,p,]
                     )
        }))
      }))
    }))
    
    #mwh_scale <- round(co2_offset_ci[[r]][Percent == x, mean(Inc_Solar), by = Start][, max(V1)] / 200, 1)
    #co2_scale <- round(-co2_offset_ci[[r]][Percent == x, sum(CO2_offset), by = .(Hour, Samp)
    #][, .(lower = quantile(V1, 0.025)), by = .(Hour)][, min(lower)] / 200, 1)
    
    p2b_lags<-  ggplot(co2_offset_ci[[r]][Percent == 0.15, sum(CO2_offset), by = .(Hour, Samp)
    ][, .(mean = mean(V1), lower = quantile(V1, 0.025), #changed to median
          upper = quantile(V1, 0.975)), by = .(Hour)]) +
      #geom_line(data = co2_offset_ci[["CAL"]][Percent == x, mean(Inc_Solar), by = Start],
      #          aes(x = Start, y = V1/mwh_scale),
      #          size = 1, linetype = 2) +
      geom_ribbon(aes(x = Hour, ymin = lower, ymax = upper), fill = "grey") +
      geom_line(aes(x = Hour, y = mean), size = 1) +
      #geom_hline(yintercept = 0, size = 1) +
      theme_minimal(base_size = 24) +
      theme(legend.position = "top",
            legend.key.width = unit(2, "cm"),
            axis.line = element_line(),
            axis.ticks = element_line()) +
      scale_x_continuous(expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                         labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
      scale_y_continuous(limits = c(-650, 0), breaks = c(-600, -500, -400, -300, -200, -100, 0)) +
      scale_color_discrete(labels = paste0(c(5, 10, 15, 20), "%")) +
      labs(x = "Hour of the day", 
           y = TeX("$\\leftarrow$ CO$_2$ change (metric tons)"), #"<- Total CO2 Offset (metric tons) \t|\t Increase in Median Solar (MWh) ->",
           color = "Solar capacity increase", 
           title = paste0(lags, " lag-hours considered"))#paste0(r, ": CO2 offset for alternate policy changes in solar capacity"))
    
    
    #png(paste(r, "_cumulative.png", sep=""), width=1000, height=800)
    png(paste("output/", r, "_lag", lags, "_cumulative_median.png", sep=""), width=800, height=400)
    print(p2b_lags)
    dev.off()
    
    #daily mean
    #daily = sum(co2_offset_ci[[r]][Percent == 0.15, sum(CO2_offset), by = .(Hour, Samp)][, .(mean = mean(V1), lower = quantile(V1, 0.025), upper = quantile(V1, 0.975)), by = .(Hour)]$mean)
    #print(paste(daily, "CO2 offset estimated considering ", lags, " lag-hours", sep=""))
  }
  
}
```


```{r}
#Sensitivity Analysis with varying lag hours (for California->Southwest)
countries <- c('CAL', 'CAR', 'CENT', 'FlA', 'MIDA', 'MIDW', 
               'NE', 'NW', 'SE', 'SW', 'TEN', 'TEX')
solar_cap <- fread("data/solar_capacity_factor.csv")
# Change hour to 0-23 for wraparound
solar_cap[, Hour := Hour %% 24]
# Rename Florida region in capacity file
solar_cap[Region == "FLA", Region := "FlA"]


# Policy changes
perc_increase <- c(.05, .1, .15, .2)


# Region to plot (for Fig 2a and 2b)
r <- "CAL"
r2 <- "SW"


# DLM at x=0.15 increase for Fig 2a and 2b
x = 0.15 #solar increase %

#Specify year
y <- "2022"

# Get solar capacity and hourly median for year 2022
hr_cap <- solar_cap[Region == r & year == y]
# Load DLNM model
#load(paste0("tdlm_models/mlist_", r, ".Rda"))
#d <- combine.models(dm_mlist)

lag_vector = c(4,8,12,16)

for (lags in lag_vector){
  models_filename = paste0("tdlm_models_lag",lags)
  load(paste0(models_filename,"/dm_mod_", r,"-",r2, ".Rda"))
  d<- dm
  # Loop over hours and calculate CO2 offset for increased solar by x%
  perc_increase <- c(x)
  co2_offset <- rbindlist(lapply(0:23, function(h) {
    m <- hr_cap[Hour == h, solar_median_mwh]
    s <- summary(d, cenval = m,
                 pred.at = m * (1 + perc_increase), 
                 mcmc = TRUE, verbose = F)
    # Loop over percent increase, record lagged CO2 offset
    rbindlist(lapply(1:length(s$pred.at), function(p) {
      data.table(Region = r,
                 ImportingRegion = r2, 
                 Start = h,
                 Med_Solar = m,
                 Hour = h:(h+lags) %% 24, # Loop back to 0 if Hour > 23
                 Wraparound = c(rep(0, sum(h:(h+lags) < 24)),
                                rep(1, sum(h:(h+lags) >= 24))), # Needed for line break in hourly lag plot
                 Percent = perc_increase[p],
                 Inc_Solar = m*perc_increase[p],
                 CO2_offset = -sapply(1:(lags+1), function(l) median(s$dlm_mcmc[l,p,])), #-s$matfit[p,] #-s$matfit[p,],
                 C_lower = -s$ciupper[p,],
                 C_upper = -s$cilower[p,])
    }))
  }))
  
  p2a_lags<- ggplot(co2_offset[Percent == x & Hour %in% c(12, 20)]) +
  geom_line(aes(x = Start, y = CO2_offset, 
                color = factor(Hour),
                group = paste0(Hour, Wraparound)), size = 1) +
  geom_ribbon(data = co2_offset[Percent == x & Hour %in% c(12, 20)], 
              aes(x = Start, ymin = C_lower, ymax = C_upper, 
              fill = factor(Hour)), alpha=0.3, show.legend = FALSE) +
  geom_point(data = co2_offset[Percent == x & Hour %in% c(12, 20) & Start==Hour, last(CO2_offset), by = Hour],
             aes(x = Hour, y = V1, color = factor(Hour)), size = 2) +
  theme_minimal(base_size = 24) +
  theme(legend.position = "right",
        legend.key.height = unit(1.5, "cm"),
        legend.text = element_text(size = 18),  
        legend.title = element_text(size = 18)) +
  scale_x_continuous(limits=c(0,23), expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
  scale_y_continuous(limits = c(-700, 0.1))+
  scale_color_manual(values = c(2,3),
    labels = c("12:00", "20:00")) +
  labs(x = "Hour of the day", y = TeX("CO$_2$ change (metric tons)"), 
       color = "Observation\nhour",
       title = paste0(lags, " lag-hours considered"))

  png(paste("output/", r,"-",r2, "_lag", lags, "_hourly_together_median.png", sep=""), width=800, height=400)
  print(p2a_lags)
  dev.off() 
  
  #CO2effect_for_an_hour = co2_offset[Percent == x & Hour %in% c(4, 8, 12,16,20), sum(CO2_offset), by="Hour"]
  #print(paste0("CO2 offset for increasing solar at 4 pm = ", CO2effect_for_an_hour))
  
  #Figure cumulative
  co2_offset_ci <- list()
  for (r in "CAL") {
    # Get solar capacity and hourly median for year 2022
    hr_cap <- solar_cap[Region == r & year == 2022]
    # Load DLNM model
    #load(paste0("tdlm_models/mlist_", r, ".Rda"))
    #d <- combine.models(dm_mlist)
    load(paste0(models_filename,"/dm_mod_", r,"-",r2, ".Rda"))
    d<- dm
    # Loop over hours and calculate CO2 offset for increased solar by x%
    co2_offset_ci[[r]] <- rbindlist(lapply(0:23, function(h) {
      m <- hr_cap[Hour == h, solar_median_mwh]
      s <- summary(d, cenval = m,
                   pred.at = m * (1 + perc_increase), 
                   mcmc = TRUE, verbose = F)
      # Loop over percent increase, record lagged CO2 offset
      rbindlist(lapply(1:length(s$pred.at), function(p) {
        rbindlist(lapply(1:(lags+1), function(l) {
          data.table(Region = r,
                     Start = h,
                     Hour = (h+l-1) %% 24, # Loop back to 0 if Hour > 23
                     Percent = perc_increase[p],
                     Inc_Solar = m * (1 + perc_increase[p]) - m,
                     Samp = 1:d$mcmcIter,
                     CO2_offset = -s$dlm_mcmc[l, p,])
        }))
      }))
    }))
    
    #mwh_scale <- round(co2_offset_ci[[r]][Percent == x, mean(Inc_Solar), by = Start][, max(V1)] / 200, 1)
    #co2_scale <- round(-co2_offset_ci[[r]][Percent == x, sum(CO2_offset), by = .(Hour, Samp)
    #][, .(lower = quantile(V1, 0.025)), by = .(Hour)][, min(lower)] / 200, 1)
    
    p2b_lags<-  ggplot(co2_offset_ci[[r]][Percent == 0.15, sum(CO2_offset), by = .(Hour, Samp)
    ][, .(mean = median(V1), lower = quantile(V1, 0.025), #changed to median
          upper = quantile(V1, 0.975)), by = .(Hour)]) +
      #geom_line(data = co2_offset_ci[["CAL"]][Percent == x, mean(Inc_Solar), by = Start],
      #          aes(x = Start, y = V1/mwh_scale),
      #          size = 1, linetype = 2) +
      geom_ribbon(aes(x = Hour, ymin = lower, ymax = upper), fill = "grey") +
      geom_line(aes(x = Hour, y = mean), size = 1) +
      #geom_hline(yintercept = 0, size = 1) +
      theme_minimal(base_size = 24) +
      theme(legend.position = "top",
            legend.key.width = unit(2, "cm"),
            axis.line = element_line(),
            axis.ticks = element_line()) +
      scale_x_continuous(expand = c(0, 0), breaks = c(0, 4, 8, 12, 16, 20),
                         labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00")) +
      scale_y_continuous(limits = c(-800, 0), breaks = c(-800, -700, -600, -500, -400, -300, -200, -100, 0)) +
      scale_color_discrete(labels = paste0(c(5, 10, 15, 20), "%")) +
      labs(x = "Hour of the day", 
           y = TeX("$\\leftarrow$ CO$_2$ change (metric tons)"), #"<- Total CO2 Offset (metric tons) \t|\t Increase in Median Solar (MWh) ->",
           color = "Solar capacity increase", 
           title = paste0(lags, " lag-hours considered"))#paste0(r, ": CO2 offset for alternate policy changes in solar capacity"))
    
    
    #png(paste(r, "_cumulative.png", sep=""), width=1000, height=800)
    png(paste("output/", r, "-",r2, "_lag", lags, "_cumulative_median.png", sep=""), width=800, height=400)
    print(p2b_lags)
    dev.off()
    
    #daily mean
    #daily = sum(co2_offset_ci[[r]][Percent == 0.15, sum(CO2_offset), by = .(Hour, Samp)][, .(mean = mean(V1), lower = quantile(V1, 0.025), upper = quantile(V1, 0.975)), by = .(Hour)]$mean)
    #print(paste(daily, "CO2 offset estimated considering ", lags, " lag-hours", sep=""))
  }
  
}


```
```{r}
#For Table S2 in the Supplementary Materials
countries <- c('CAL', 'CAR', 'CENT', 'FlA', 'MIDA', 'MIDW', 'NE', 'NW', 'SE', 'SW', 'TEN', 'TEX')

lags <- 12 # Set the number of lag hours (update this for analyse sensitivity to lag hours)

data_list <- list()
lag_data_list <- list()
lead_data_list <- list()
complete_list <- list()
for (i in (1:length(countries))){
  print(countries[i])
  data <- read.csv(paste("data/", countries[i], ".csv", sep=""))
  setDT(data)
  
  # Create lagged data
  if (!file.exists(paste0("tdlm_models_lag",lags))){
    dir.create(paste0("tdlm_models_lag",lags))
  }
  models_filename = paste0("tdlm_models_lag",lags)
  lag_dat <- as.matrix(data[, shift(NG..SUN, 0:lags, type = "lag")])
  #lead_dat <- as.matrix(data[, shift(NG..SUN, 10:1, type = "lead")])
  #lag_dat <- cbind(lead_dat, lag_dat)
  # data[, demand_ma10 := rollmean(demand, k = 10)] -- CHECK! 
  # Complete data rows
  
  data_list[[i]]<-data #Note: this data may contain NA values
  lag_data_list[[i]]<- lag_dat #Note: this data may contain NA values
  #lead_data_list[[i]]<-lead_dat
  complete_list[[i]] <- which(complete.cases(lag_dat) & 
                                complete.cases(data[, .(CO2.Emissions.Generated, demand)]))
  
}

# Determine final datasets (remove missing)
final_data <- lapply(1:length(countries), function(i) {
  data_list[[i]][complete_list[[i]]]
})
names(final_data) <- countries

# Calculate hourly medians
solar_medians <- lapply(final_data, function(f) {
  f[year == 2022 & Hour %in% 1:24, median(NG..SUN), by = Hour] })

#solar increase
perc_inc <- lapply(c(0.05, 0.1, 0.15, 0.2), function(p) {
  sapply(solar_medians, function(s) sum(s$V1 * p)) })

#15%increase
inc_15 = perc_inc[[3]] * 365 / 1000000

#CO2 for  COL with 15% increase
factors_COL=unlist(lapply(countries, function(c){(final_data)[[c]][year == 2022 & NG..COL != 0, 
             mean(CO2.Emissions..COL / NG..COL, na.rm = T)]}))
CO2_15_COL=factors_COL*inc_15

#CO2 for  NG with 15% increase
factors_NG=unlist(lapply(countries, function(c){(final_data)[[c]][year == 2022 & NG..NG != 0, 
             mean(CO2.Emissions..NG / NG..NG, na.rm = T)]}))
CO2_15_NG=factors_NG*inc_15

#CO2 for  OIL with 15% increase
factors_OIL=unlist(lapply(countries, function(c){(final_data)[[c]][year == 2022 & NG..OIL != 0, 
             mean(CO2.Emissions..OIL / NG..OIL, na.rm = T)]}))
CO2_15_OIL=factors_OIL*inc_15

maxtab = cbind(inc_15, factors_COL,CO2_15_COL, factors_NG, CO2_15_NG, factors_OIL, CO2_15_OIL)
maxtab = rbind(maxtab, "total"=colSums(maxtab, na.rm=T))
write.csv(maxtab, "output/theoreticalmaxCO2.csv")

```