---
title: "Sectoral_analysis"
output: html_document
date: "2025-10-23"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "C:/Users/ignasicortesar/OneDrive - Delft University of Technology/Paper3_Private_Adaptation_CGE")
```

```{r, message=FALSE}
library(tidyverse)
library(readxl)
library(ggplot2)
library(gamstransfer)

```

## Setup Components

```{r}
# define function to quickly read in specific variables from the gdx

get_gdx_records <- function(gdx_list, var_name) {
  if (!var_name %in% names(gdx_list)) {
    stop(paste("Variable", var_name, "not found in GDX data."))
  }
  gdx_list[[var_name]]$records
}


```


```{r}

# for mapping purposes

library(eurostat)
library(sf)
library(RColorBrewer)
library(giscoR)

# add Croatia

geodata <- get_eurostat_geospatial(output_class = "sf", resolution = "1", nuts_level = 2, year = 2010, cache = TRUE, update_cache = TRUE)
geodata_1 <- get_eurostat_geospatial(output_class = "sf", resolution = "1", nuts_level = 1, year = 2013, cache = TRUE, update_cache = TRUE)
geodata_1 <- subset(geodata_1, id == "HR0")
geodata_1$id[1] <- "HR00"
geodata_1$NUTS_ID[1] <- "HR00"
geodata_1$FID[1] <- "HR00"
geodata_1$geo[1] <- "HR00"
geodata_1$LEVL_CODE[1] <- 2
geodata_final <- rbind(geodata, geodata_1)
geodata <- geodata_final
```

# Read in all data for 3C scenarios

```{r}

# 3C with adaptation (Ad_Growth disabled)
results_3C_Ad1_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_0_AdSpike_0_RegDiff_0_Optimal_0SMALL.gdx")

# 3C with adaptation (optimal)
results_3C_AdOPT_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_0_AdSpike_0_RegDiff_0_Optimal_1SMALL.gdx")

# 3C without adaptation
results_3C_Ad0_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_0_AdGrowth_0_AdSpike_0_RegDiff_0_Optimal_0SMALL.gdx")

# no climate change
results_nodamages <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_NODAMAGES.gdx")

```

# Read in regional clusters

Cluster 1: High Damage, Service Oriented
Cluster 2: Low Damage, Service Oriented
Cluster 3: High Damage, Non-Service Oriented
Cluster 4: Low Damage, Non-Service Oriented

```{r}

cluster_classification <- readxl::read_excel("output-processor/cluster_definition.xlsx", sheet = "clustering_output", range = "A1:B272")

cluster_classification <- cluster_classification %>% 
  rename(geo = "cnt")

```

# Extract GVA and Output data (basic adaptation)

```{r}
# Counterfactual with basic adaptation (G0, R0)

GVA_3C_Ad1 <- get_gdx_records(results_3C_Ad1_gdx, "GVA_t") %>% 
  rename(GVA_Ad1 = "value")
Output_3C_Ad1 <-get_gdx_records(results_3C_Ad1_gdx, "XD_dyn") %>% 
  rename(XD_Ad1 = "value")

# Baseline without adaptation

GVA_3C_Ad0 <- get_gdx_records(results_3C_Ad0_gdx, "GVA_t") %>% 
  rename(GVA_Ad0 = "value")
Output_3C_Ad0 <-get_gdx_records(results_3C_Ad0_gdx, "XD_dyn") %>% 
  rename(XD_Ad0 = "value")

```


# Compute relative changes

```{r}

# compute relative output changes
GVA_3C_rel <- merge(GVA_3C_Ad1, GVA_3C_Ad0, by = c("cnt", "t", "sec")) %>% 
  transform(GVA_change = 100 * (GVA_Ad1-GVA_Ad0)/GVA_Ad0) %>% 
  rename(geo = "cnt") %>% 
  rename(Year = "t") %>% 
  rename(Sector = "sec") %>% 
  merge(cluster_classification, by = "geo")

# compute relative GVA changes
Output_3C_rel <- merge(Output_3C_Ad1, Output_3C_Ad0, by = c("cnt", "t", "sec")) %>% 
  transform(Output_change = 100 * (XD_Ad1-XD_Ad0)/XD_Ad0) %>% 
  rename(geo = "cnt") %>% 
  rename(Year = "t") %>% 
  rename(Sector = "sec") %>% 
  merge(cluster_classification, by = "geo")

#choose a year for analysis: 2100

Output_3C_rel_2100 <- Output_3C_rel %>% 
  filter(Year == 2100) %>% 
  select(-Year)

GVA_3C_rel_2100 <- GVA_3C_rel %>% 
  filter(Year == 2100) %>% 
  select(-Year)

```

# Violin plots

```{r}

# Create a violin plot for Output in 2100 using ggplot2
Output_change_2100_violins = ggplot(Output_3C_rel_2100, aes(x = Sector, y = Output_change, fill = Sector)) +
  geom_violin(trim = FALSE) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.7) +
  #geom_point(data = summary_stats,aes(x = Sector, y = mean), color = "black", size = 2) +
  scale_fill_viridis_d() +
  coord_cartesian(ylim = c(-2,2)) +
  labs(title = "Changes in Output by Sector in 2100",
       x = "Economic Sector",
       y = "Change in Output (%)") +
  theme_minimal() +
  facet_wrap(~ Cluster, ncol = 2, labeller = labeller(Cluster = function(x) paste("Cluster", x)))


# Create a violin plot for GVA in 2100 using ggplot2
GVA_3C_rel_2100_violins = ggplot(GVA_3C_rel_2100, aes(x = Sector, y = GVA_change, fill = Sector)) +
  geom_violin(trim = FALSE) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.7) +
  scale_fill_viridis_d() +
  coord_cartesian(ylim = c(-2,2)) +
  labs(title = "Changes in GVA by Sector in 2100",
       x = "Economic Sector",
       y = "Change in GVA (%)") +
  theme_minimal() +
  facet_wrap(~ Cluster, ncol = 2, labeller = labeller(Cluster = function(x) paste("Cluster", x)))



```

Write up summary stats and export

```{r}

# Output

summary_stats_output_Table4 <- Output_3C_rel_2100 %>%
  group_by(Cluster, Sector) %>%
  summarise(
    n = n(),
    mean = mean(Output_change, na.rm = TRUE),
    median = median(Output_change, na.rm = TRUE),
    sd = sd(Output_change, na.rm = TRUE),
    q25 = quantile(Output_change, 0.25, na.rm = TRUE),
    q75 = quantile(Output_change, 0.75, na.rm = TRUE),
    min = min(Output_change, na.rm = TRUE),
    max = max(Output_change, na.rm = TRUE)
  ) %>%
  arrange(Cluster, Sector)

writexl::write_xlsx(Output_3C_rel_2100, path = "output-analysis/output_results_A1R0.xlsx")
writexl::write_xlsx(summary_stats_output_Table4, path = "output-analysis/summary_stats_output.xlsx")


# GVA

summary_stats_GVA <- GVA_3C_rel_2100 %>%
  group_by(Cluster, Sector) %>%
  summarise(
    n = n(),
    mean = mean(GVA_change, na.rm = TRUE),
    median = median(GVA_change, na.rm = TRUE),
    sd = sd(GVA_change, na.rm = TRUE),
    q25 = quantile(GVA_change, 0.25, na.rm = TRUE),
    q75 = quantile(GVA_change, 0.75, na.rm = TRUE),
    min = min(GVA_change, na.rm = TRUE),
    max = max(GVA_change, na.rm = TRUE)
  ) %>%
  arrange(Cluster, Sector)


writexl::write_xlsx(summary_stats_GVA, path = "output-analysis/summary_stats_GVA.xlsx")

```

