# Metadata ----
## written by     : Fabio Ascione & Benedikt Bruckner
## description    : plot CO2e footprints of tourism from ICIO model
##                  run after cons_td_analysis.R
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#+++++++++++++++++++++++++++++++++++++++++
# 0 - Packages ----
#+++++++++++++++++++++++++++++++++++++++++
librarian::shelf(dplyr, readr, tidyr, readxl, unpivotr, stringr, reshape2, purrr, xlsx, 
                 DescTools, decompr, data.table, janitor, here, ggplot2, forcats)

### Set path
path <- here()

### PALETTES
pal <- met.brewer("Hokusai3")
extended_pal <- c(pal, "#6BA87F")


#+++++++++++++++++++++++++++++++++++++++++
# 1 - LOAD DATA ----
#+++++++++++++++++++++++++++++++++++++++++

## ICIO
# full footprints (total final demand, no tsa, all country-sectors)
load(paste0(path,"/data/derived/TD-consumption/fpfin_icio.Rda"))
# tourism footprint (tsa demand for germany in base price, total final demand for other country-sectors)
load(paste0(path,"/data/derived/TD-consumption/fpfin_tour_icio_bp.Rda"))
# tsa shares
load(paste0(path,"/data/derived/TD-consumption/tsa_shares_icio_fin.Rda"))
# tsa data
load(paste0(path, '/data/derived/TD-consumption/tsa_icio_fin_bp.Rda'))
# icio tsa correspondence
load(file = paste0(path, '/data/derived/TD-consumption/icio_tsa_corresp.Rda'))


## EXIOBASE

# full footprints (tsa demand for germany in base price, total final demand for other country-sectors)
load(paste0(path,"/data/derived/TD-consumption/ghg_fp_exio_bp_tot.Rda"))
# tourism footprint (tsa shares for germany, total final demand for other country-sectors)
load(paste0(path,"/data/derived/TD-consumption/ghg_fp_exio_shares_tot.Rda"))
# tsa shares
load(paste0(path,"/data/derived/TD-consumption/tsa_shares_exiobase_fin.Rda"))
# tsa data
load(paste0(path, '/data/derived/TD-consumption/tsa_exio_fin_bp.Rda'))
# exio tsa correspondence
load(file = paste0(path, '/data/derived/TD-consumption/exio_tsa_corresp.Rda'))



#+++++++++++++++++++++++++++++++++++++++++
# 2 - MANIPULATE DATA ----
#+++++++++++++++++++++++++++++++++++++++++

#+++++++++++++++++++++++++++++++++++++++++
## 2.1 - calculate footprint using TSA SHARES ---- ICIO ----
#+++++++++++++++++++++++++++++++++++++++++


## calculate footprint with tsa shares
# post multiply total footprint (total final demand, without tsa shares) with tsa shares
fpfin_icio_shares <- fpfin_icio %>% 
  mutate(year=as.numeric(year)) %>% 
  full_join(tsa_shares_icio_fin, by= c('year', 'isic3'='ICIO_code')) %>% 
  # multiply footprints with tsa shares
  mutate(fp_tot_tour_share=fp_tot*tsa_share_in_letzter_verwendung,
         fp_tot_for_tour_share=fp_tot_for*tsa_share_in_letzter_verwendung,
         fp_tot_dom_tour_share=fp_tot_dom*tsa_share_in_letzter_verwendung,
         fp_dir_tour_share=fp_dir*tsa_share_in_letzter_verwendung,
         fp_dir_dom_tour_share=fp_dir_dom*tsa_share_in_letzter_verwendung,
         fp_dir_for_tour_share=fp_dir_for*tsa_share_in_letzter_verwendung,
         fp_ind_tour_share=fp_ind*tsa_share_in_letzter_verwendung,
         fp_ind_dom_tour_share=fp_ind_dom*tsa_share_in_letzter_verwendung,
         fp_ind_for_tour_share=fp_ind_for*tsa_share_in_letzter_verwendung,
         fp_ind_energy_tour_share=fp_ind_energy*tsa_share_in_letzter_verwendung,
         fp_ind_energy_dom_tour_share=fp_ind_energy_dom*tsa_share_in_letzter_verwendung,
         fp_ind_energy_for_tour_share=fp_ind_energy_for*tsa_share_in_letzter_verwendung
         ) %>% 
  select(year, iso3c, isic3, isic3_label=ICIO_label, contains('fp'))
# prepare tsa products codes to match icio sectors
icio_tsa_corresp_prep <- icio_tsa_corresp %>% mutate(tsa_industry_label_deu=
                              case_when(ICIO_code == 'I' ~ 'Beherbergungs- und Gaststättenleistungen',
                                        ICIO_code == 'H49' ~ 'Eisenbahn-, Straßen- und Nahverkehrsleistungen',
                                        ICIO_code == 'N' ~ 'Mietfahrzeuge, Reisebüros und -veranstalter',
                                        .default = tsa_industry_label_deu)) %>% 
  group_by(ICIO_code, tsa_industry_label_deu) %>%  # Group by the two columns with duplicates
  summarise(
    tsa_codes = paste(tsa_industry_code_deu, collapse = ", "),  # Collapse the codes in col3
    .groups = "drop"  # Remove grouping structure after summarise
  )
# add tsa product codes to icio footprints
fpfin_icio_shares_prep <- fpfin_icio_shares %>% full_join(icio_tsa_corresp_prep, by=c('isic3'='ICIO_code')) %>% 
  relocate(c(tsa_industry_label_deu, tsa_codes), .after = isic3_label)



#+++++++++++++++++++++++++++++++++++++++++
## 2.2 - calculate footprint using TSA SHARES ---- EXIOBASE ----
#+++++++++++++++++++++++++++++++++++++++++

# convert list into panel dataframe
ghg_fp_exio_shares_tot_panel <- map_df(ghg_fp_exio_shares_tot, ~as_tibble(.x), .id = 'year') %>% 
  separate(rowname, into = c('iso2c','exio_code'), sep = '_', extra = 'merge') %>%
  rename(fp_tot_tour_share=last_col())

# divide tourism footprint (with tsa shares) by tsa shares to get total final demand
fpfin_exio_shares <- ghg_fp_exio_shares_tot_panel %>% 
  mutate(year=as.numeric(year), exio_code=as.numeric(exio_code)) %>%  
  full_join(tsa_shares_exiobase_fin,by= c('year', 'exio_code'='EXIOBASE_pxp_code')) %>% 
  # divide tourism footprints by tsa shares to get total footprint
  mutate(fp_tot=fp_tot_tour_share/tsa_share_in_letzter_verwendung) %>% 
  select(year, iso2c, exio_code, exio_label=EXIOBASE_pxp_label, contains('fp')) %>% 
  # convert to million tons
  mutate(across(contains('fp'), ~ . / 1e+09))

# prepare tsa products codes to match exio sectors
exio_tsa_corresp_prep <- exio_tsa_corresp %>% 
  mutate(tsa_industry_label_deu=if_else(EXIOBASE_pxp_code == '156', 'Beherbergungs- und Gaststättenleistungen', tsa_industry_label_deu)) %>%
  group_by(EXIOBASE_pxp_code, EXIOBASE_pxp_label,  tsa_industry_label_deu) %>%
  summarise(
    tsa_codes = paste(tsa_industry_code_deu, collapse = ", "),
    .groups = "drop"
  )

# add tsa product codes to exio footprints
fpfin_exio_shares_prep <- fpfin_exio_shares %>% 
  full_join(exio_tsa_corresp_prep, by=c('exio_code'='EXIOBASE_pxp_code', 'exio_label'='EXIOBASE_pxp_label')) %>%
  relocate(c(tsa_industry_label_deu, tsa_codes), .after = exio_label)




#+++++++++++++++++++++++++++++++++++++++++
# 3 - PLOT DATA ----
#+++++++++++++++++++++++++++++++++++++++++


#+++++++++++++++++++++++++++++++++++++++++
## 3.0 - Symposium ----
#+++++++++++++++++++++++++++++++++++++++++

### set ggplot theme settings
theme_fin <- theme_minimal() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
    plot.title = element_text(face = "bold", size = 14),
    axis.title = element_text(size = 11),
    legend.position = "top",
    legend.title=element_blank(),
    axis.text.x = element_text(angle=45,hjust=1)
  )

### 1) TOTAL MIO TONS EMISSIONS TOURISM SCOPE 1-3
fig_totals <- fpfin_icio_shares_prep %>%
  filter(iso3c %in% c('DEU')) %>%
  group_by(year) %>%
  summarise(scope1=sum(fp_dir_dom_tour_share, na.rm = T), 
            scope2=sum(fp_ind_energy_tour_share, na.rm = T),
            scope3=sum(fp_ind_tour_share, na.rm = T)) %>%
  filter(year<2020) %>%
  pivot_longer(., c(scope1:scope3)) %>%
  ggplot(., aes(x=year, y=value, fill=name)) +
  geom_col(stat = 'identity') +
  ylim(y=0,x=100)+
  scale_fill_manual(values = pal)+
  ylab('')+xlab('')+
  theme_fin
fig_totals
ggsave(fig_totals, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_totals_scopes.png", width = 8.5, height = 5, dpi = 600, device = "png")

### 2) SHARE TOURISM Mio tons CO2e IN TOTAL MIO TONS CO2E 
fig_shares <- fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>%
  group_by(year) %>% 
  summarise(fp_tot=sum(fp_tot, na.rm = T),
            scope1=sum(fp_dir_dom_tour_share, na.rm = T), 
            scope2=sum(fp_ind_energy_tour_share, na.rm = T),
            scope3=sum(fp_ind_tour_share, na.rm = T)) %>% 
  mutate(fp_scope1_share = scope1 / fp_tot,
          fp_scope2_share = scope2 / fp_tot,
         fp_scope3_share = scope3 / fp_tot) %>% 
  pivot_longer(., c(fp_scope1_share:fp_scope3_share)) %>%
  filter(year<2020) %>%
  ggplot(., aes(x=year, y=value, fill=name)) +
  geom_col(stat = 'identity') +
  ylab('')+xlab('')+
  scale_fill_manual(values = pal)+
  theme_fin
fig_shares
ggsave(fig_shares, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_shares_scopes.png", width = 8.5, height = 5, dpi = 600, device = "png")

### 3) TOTAL MIO TONS EMISSIONS TOURISM SCOPE 1-3 WITH DOMESTIC FOREIGN SPLIT
fig_totals_scopes_origin <- fpfin_icio_shares_prep %>%
  filter(iso3c %in% c('DEU')) %>%
  group_by(year) %>%
  summarise(scope1=sum(fp_dir_dom_tour_share, na.rm = T), 
            scope2_dom=sum(fp_ind_energy_dom_tour_share, na.rm = T),
            scope2_for=sum(fp_ind_energy_for_tour_share, na.rm = T),
            scope3_dom=sum(fp_ind_dom_tour_share, na.rm = T),
            scope3_for=sum(fp_ind_for_tour_share, na.rm = T)) %>% 
  filter(year<2020) %>%
  pivot_longer(., c(scope1:scope3_for)) %>% 
  ggplot(., aes(x=year, y=value, fill=name)) +
  geom_col(stat = 'identity') +
  ylim(y=0,x=100)+
  ylab('')+xlab('')+
  scale_fill_manual(values = pal)+
  theme_fin
fig_totals_scopes_origin
ggsave(fig_totals_scopes_origin, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_totals_scopes_origin.png", width = 8.5, height = 5, dpi = 600, device = "png")

### 4) EMISSIONS BY TOURISTIC PRODUCT / SECTOR
fig_prod_scopes_origin <- fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>%
  group_by(year, tsa_industry_label_deu, tsa_codes) %>%
  summarise(fp_tot_tour_share=sum(fp_tot_tour_share, na.rm = T),
            scope1=sum(fp_dir_dom_tour_share, na.rm = T), 
            scope2_dom=sum(fp_ind_energy_dom_tour_share, na.rm = T),
            scope2_for=sum(fp_ind_energy_for_tour_share, na.rm = T),
            scope3_dom=sum(fp_ind_dom_tour_share, na.rm = T),
            scope3_for=sum(fp_ind_for_tour_share, na.rm = T)) %>%
  filter(year==2019) %>%
  pivot_longer(., c(scope1:scope3_for)) %>% 
  ggplot(., aes(x=reorder(tsa_industry_label_deu,-value), y=value, fill=name)) +
  geom_col(stat = 'identity') +
  scale_fill_manual(values = pal)+
  theme_fin
fig_prod_scopes_origin
ggsave(fig_prod_scopes_origin, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_prod_scopes_origin.png", width = 8.5, height = 5, dpi = 600, device = "png")

### 5) TYPE OF TOURISTS
# calculate shares from tsa data
tsa_origin_type_shares <- tsa_icio_fin_bp %>%
  mutate(across(verwendung_tour_business_mill_eur_bp:last_col(), ~ . / verwendung_tour_total_mill_eur_bp)) %>%
  mutate(iso3c='DEU') %>% rename(isic3=ICIO_code) %>% 
  filter(year>2015)
fpfin_icio_origin_type <- fpfin_icio_shares_prep %>%
  left_join(tsa_origin_type_shares) %>%  
  filter(iso3c %in% c('DEU')) %>%
  # filter(str_detect(isic3, paste(tsa_sectors, collapse = "|"))) %>% 
  mutate(fp_tot_tour_business=fp_tot_tour_share*verwendung_tour_business_mill_eur_bp,
         fp_tot_tour_private=fp_tot_tour_share*verwendung_tour_private_mill_eur_bp,
         fp_tot_tour_incoming=fp_tot_tour_share*verwendung_tour_incoming_mill_eur_bp,
         fp_tot_tour_domestic=fp_tot_tour_share*verwendung_tour_domestic_mill_eur_bp) %>%
  select(year, iso3c, isic3, isic3_label, tsa_codes, tsa_industry_label_deu, fp_tot_tour_share, fp_tot_tour_business, fp_tot_tour_private, fp_tot_tour_incoming, fp_tot_tour_domestic)

# private vs. business
fig_priv_bus <- fpfin_icio_origin_type %>%
  group_by(year, iso3c) %>%
  summarise(across(c(fp_tot_tour_share:last_col()), ~ sum(., na.rm = T))) %>%
  pivot_longer(., c(fp_tot_tour_business:fp_tot_tour_private)) %>% 
  ggplot(., aes(x=year,y=value, fill=name))+
  geom_col(stat = 'identity')+
  theme_fin
fig_priv_bus
ggsave(fig_priv_bus, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_priv_bus.png", width = 8.5, height = 5, dpi = 600, device = "png")

# incoming vs. domestic
fig_origin <- fpfin_icio_origin_type %>%
  group_by(year, iso3c) %>%
  summarise(across(c(fp_tot_tour_share:last_col()), ~ sum(., na.rm = T))) %>%
  pivot_longer(., c(fp_tot_tour_incoming:fp_tot_tour_domestic)) %>% 
  ggplot(., aes(x=year,y=value, fill=name))+
  geom_col(stat = 'identity')+
  theme_fin
fig_origin
ggsave(fig_origin, file = "C:/pik/dkt/fig/symposium_presentation/fp_co2e_origin.png", width = 8.5, height = 5, dpi = 600, device = "png")



#+++++++++++++++++++++++++++++++++++++++++
## 3.1 - ICIO ----
#+++++++++++++++++++++++++++++++++++++++++

## TOTAL TOURISTIC; DOMESTIC VS FOREIGN; DIRECT VS INDIRECT

# tourism vs. rest footprint on German sectors - lines in absolute mill. t co2e
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c) %>% 
  summarise(fp_tot_sum=sum(fp_tot, na.rm = T),fp_tot_tour_share_sum=sum(fp_tot_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_tot_sum, colour='total footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_tour_share_sum, colour='tourism footprint'))

# tourism vs. rest footprint on German sectors - lines in % share
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c) %>% 
  summarise(fp_tot_sum=sum(fp_tot, na.rm = T),fp_tot_tour_share_sum=sum(fp_dir_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_tot_tour_share_sum/fp_tot_sum, colour='tourism footprint/total footprint')) +
  geom_line()

# absolute tourism footrpint by sector - lines in absolute mill. t co2e
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  ggplot(., aes(x=year, y=fp_tot, colour='total footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_tour_share, colour='tourism footprint')) + 
  facet_wrap(~isic3, scales = 'free')

# absolute tourism footrpint by tsa product - lines in absolute mill. t co2e
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>%
  group_by(year, tsa_industry_label_deu, tsa_codes) %>% 
  summarise(fp_tot=sum(fp_tot, na.rm = T),
            fp_tot_tour_share=sum(fp_tot_tour_share, na.rm = T)) %>%
  ggplot(., aes(x=year, y=fp_tot, colour='total footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_tour_share, colour='tourism footprint')) + 
  facet_wrap(~tsa_industry_label_deu, scales = 'free')

# share tourism footrpint by sector - lines in %
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  ggplot(., aes(x=year, y=fp_tot_tour_share/fp_tot, colour='tourism footprint/total footprint')) +
  geom_line()+
  facet_wrap(~isic3, scales = 'free')

# share tourism footrpint by tsa product - lines in %
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>%
  group_by(year, tsa_industry_label_deu, tsa_codes) %>% 
  summarise(fp_tot=sum(fp_tot, na.rm = T),
            fp_tot_tour_share=sum(fp_tot_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_tot_tour_share/fp_tot, colour='tourism footprint/total footprint')) +
  geom_line()+
  facet_wrap(~tsa_industry_label_deu, scales = 'free')

# tourism footprint domestic vs. foreign - lines in mill. t co2e
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c) %>% 
  summarise(fp_tot_for_tour_share=sum(fp_tot_for_tour_share, na.rm = T),fp_tot_dom_tour_share=sum(fp_tot_dom_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_tot_for_tour_share, colour='foreign footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_dom_tour_share, colour='domestic footprint'))

# tourism footprint direct vs. indirect - lines in mill. t co2e
fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c) %>% 
  summarise(fp_dir_tour_share=sum(fp_dir_tour_share, na.rm = T),fp_ind_tour_share=sum(fp_ind_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_dir_tour_share, colour='direct footprint')) +
  geom_line()+
  geom_line(aes(y=fp_ind_tour_share, colour='indirect footprint'))


## ORIGIN OF TOURISTS (INCOMING VS DOMESTIC) AND TYPE OF TRAVEL (PRIVATE VS BUSINESS)

# calculate shares from tsa data
tsa_origin_type_shares <- tsa_icio_fin_bp %>% 
  mutate(across(verwendung_tour_business_mill_eur_bp:last_col(), ~ . / verwendung_tour_total_mill_eur_bp)) %>% 
  mutate(iso3c='DEU') %>% rename(isic3=ICIO_code)

# add origin type shares
fpfin_icio_origin_type <- fpfin_icio_shares_prep %>% 
  full_join(tsa_origin_type_shares) %>% 
    filter(iso3c %in% c('DEU')) %>% 
  mutate(fp_tot_tour_mill_eur_business=fp_tot_tour_share*verwendung_tour_business_mill_eur_bp,
         fp_tot_tour_mill_eur_private=fp_tot_tour_share*verwendung_tour_private_mill_eur_bp,
         fp_tot_tour_mill_eur_incoming=fp_tot_tour_share*verwendung_tour_incoming_mill_eur_bp,
         fp_tot_tour_mill_eur_domestic=fp_tot_tour_share*verwendung_tour_domestic_mill_eur_bp) %>% 
  select(year, iso3c, isic3, isic3_label, tsa_industry_label_deu, tsa_codes, fp_tot_tour_mill_eur = fp_tot_tour_share,fp_tot_tour_mill_eur_business,fp_tot_tour_mill_eur_private,fp_tot_tour_mill_eur_incoming,fp_tot_tour_mill_eur_domestic)

# private vs. business
fpfin_icio_origin_type %>% 
  group_by(year, iso3c) %>% 
  summarise(across(c(fp_tot_tour_mill_eur:last_col()), ~ sum(., na.rm = T))) %>%
  filter(year>2015,year<2021) %>% 
  ggplot(., aes(x=year))+
  geom_line(aes(y=fp_tot_tour_mill_eur_business, colour='business'))+
  geom_line(aes(y=fp_tot_tour_mill_eur_private, colour='private'))

# incoming vs. domestic
fpfin_icio_origin_type %>% 
  group_by(year, iso3c) %>% 
  summarise(across(c(fp_tot_tour_mill_eur:last_col()), ~ sum(., na.rm = T))) %>% 
  filter(year>2015,year<2021) %>% 
  ggplot(., aes(x=year))+
  geom_line(aes(y=fp_tot_tour_mill_eur_incoming, colour='incoming'))+
  geom_line(aes(y=fp_tot_tour_mill_eur_domestic, colour='domestic'))



#+++++++++++++++++++++++++++++++++++++++++
## 3.2 - EXIOBASE ----
#+++++++++++++++++++++++++++++++++++++++++


## TOTAL TOURISTIC; DOMESTIC VS FOREIGN; DIRECT VS INDIRECT

# tourism vs. rest footprint on German sectors - lines in absolute mill. t co2e
fpfin_exio_shares_prep %>% 
  filter(iso2c %in% c('DE')) %>%  
  group_by(year, iso2c) %>% 
  summarise(fp_tot_sum=sum(fp_tot, na.rm = T),fp_tot_tour_share_sum=sum(fp_tot_tour_share, na.rm = T)) %>%
  ggplot(., aes(x=year, y=fp_tot_sum, colour='total footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_tour_share_sum, colour='tourism footprint'))

# tourism vs. rest footprint on German sectors - lines in % share
fpfin_exio_shares_prep %>% 
  filter(iso2c %in% c('DE')) %>% 
  group_by(year, iso2c) %>% 
  summarise(fp_tot_sum=sum(fp_tot, na.rm = T),fp_tot_tour_share_sum=sum(fp_tot_tour_share, na.rm = T)) %>% 
  ggplot(., aes(x=year, y=fp_tot_tour_share_sum/fp_tot_sum, colour='tourism footprint/total footprint')) +
  geom_line()

# absolute tourism footrpint by tsa product - lines in absolute mill. t co2e
fpfin_exio_shares_prep %>% 
  filter(iso2c %in% c('DE')) %>% 
  group_by(year, tsa_industry_label_deu, tsa_codes) %>% 
  summarise(fp_tot=sum(fp_tot, na.rm = T),
            fp_tot_tour_share=sum(fp_tot_tour_share, na.rm = T)) %>%
  ggplot(., aes(x=year, y=fp_tot, colour='total footprint')) +
  geom_line()+
  geom_line(aes(y=fp_tot_tour_share, colour='tourism footprint')) + 
  facet_wrap(~tsa_industry_label_deu, scales = 'free')

# share tourism footrpint by tsa product - lines in %
fpfin_exio_shares_prep %>% 
  filter(iso2c %in% c('DE')) %>%
  group_by(year, tsa_industry_label_deu, tsa_codes) %>% 
  summarise(fp_tot=sum(fp_tot, na.rm = T),
            fp_tot_tour_share=sum(fp_tot_tour_share, na.rm = T)) %>%
  ggplot(., aes(x=year, y=fp_tot_tour_share/fp_tot, colour='tourism footprint/total footprint')) +
  geom_line()+
  facet_wrap(~tsa_industry_label_deu, scales = 'free')



#+++++++++++++++++++++++++++++++++++++++++
## 3.3 - COMPARISON ICIO EXIOBASE ----
#+++++++++++++++++++++++++++++++++++++++++


## overall distribution plots
# exiobase
fpfin_exio_shares_prep %>% 
  filter(iso2c == 'DE', year>2015,year<2021) %>% 
  select(year, fp_tot_tour_share) %>% 
  ggplot(., aes(log(fp_tot_tour_share), group=year, colour=year)) +
  geom_density()
# icio
fpfin_icio_shares_prep %>% 
  filter(iso3c == 'DEU') %>% 
  select(year, fp_tot_tour_share) %>%
  ggplot(., aes(log(fp_tot_tour_share), group=year, colour=year)) +
  geom_density()

## aggregate results
# exiobase
fp_exio_deu_aggregate <- fpfin_exio_shares_prep %>% 
  filter(iso2c == 'DE', year>2015,year<2021) %>% 
  mutate(iso3c = 'DEU') %>% 
  relocate(iso3c, .after = iso2c) %>% 
  select(-iso2c) %>% 
  group_by(year, iso3c) %>% 
    summarise(fp_tot_exio=sum(fp_tot, na.rm = T),
            fp_tot_tour_share_exio=sum(fp_tot_tour_share, na.rm = T))
# icio & exiobase combined 
exio_icio_deu_aggregate <- fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c) %>% 
  summarise(fp_tot_icio=sum(fp_tot, na.rm = T), fp_tot_tour_share_icio=sum(fp_tot_tour_share, na.rm = T)) %>% 
  full_join(fp_exio_deu_aggregate)

## plots
# absolute values: line plot exio vs. icio
exio_icio_deu_aggregate %>% 
  pivot_longer(fp_tot_icio:last_col()) %>% 
  mutate(model=if_else(str_detect(name, 'icio'), 'icio', 'exio')) %>%
  mutate(linetype=if_else(str_detect(name, 'tour'), 'tourism footprint', 'total footprint')) %>% 
  ggplot(., aes(x=year, y=value, colour=model, linetype=linetype, group=interaction(model,linetype)))+
  geom_line(size = 1) +
  scale_linetype_manual(values = c("solid", "dashed"))

# relative shares
exio_icio_deu_aggregate %>% 
  mutate(icio_tour_share=fp_tot_tour_share_icio/fp_tot_icio,
         exio_tour_share=fp_tot_tour_share_exio/fp_tot_exio) %>% 
  ggplot(., aes(x=year, y=icio_tour_share))+
  geom_line(aes(colour='icio share'))+
  geom_line(aes(y=exio_tour_share, colour='exio share'))

## by tsa product
# exiobase
fp_exio_deu_tsa <- fpfin_exio_shares_prep %>% 
  filter(iso2c == 'DE', year>2015,year<2021) %>% 
  mutate(iso3c = 'DEU') %>% 
  relocate(iso3c, .after = iso2c) %>% 
  select(-iso2c) %>% 
  group_by(year, iso3c, tsa_codes, tsa_industry_label_deu) %>%
    summarise(fp_tot_exio=sum(fp_tot, na.rm = T),
            fp_tot_tour_share_exio=sum(fp_tot_tour_share, na.rm = T)) %>% 
  mutate(tsa_industry_label_deu=
                              case_when(tsa_codes %in%  c(3,4) ~ 'Eisenbahn-, Straßen- und Nahverkehrsleistungen',
                                        tsa_codes %in% c(7,8) ~ 'Mietfahrzeuge, Reisebüros und -veranstalter',
                                        .default = tsa_industry_label_deu)) %>%   
  group_by(year, tsa_industry_label_deu) %>%
  summarise(
    fp_tot_exio = sum(fp_tot_exio, na.rm = T),  
    fp_tot_tour_share_exio = sum(fp_tot_tour_share_exio, na.rm = T))
# icio & exiobase combined 
exio_icio_deu_tsa <- fpfin_icio_shares_prep %>% 
  filter(iso3c %in% c('DEU')) %>% 
  group_by(year, iso3c, tsa_codes, tsa_industry_label_deu) %>%
  summarise(fp_tot_icio=sum(fp_tot, na.rm = T), fp_tot_tour_share_icio=sum(fp_tot_tour_share, na.rm = T)) %>% 
  full_join(fp_exio_deu_tsa)

## plots
# absolute values: line plot exio vs. icio
exio_icio_deu_tsa %>% 
  pivot_longer(fp_tot_icio:last_col()) %>% 
  mutate(model=if_else(str_detect(name, 'icio'), 'icio', 'exio')) %>% 
  filter(str_detect(name, 'share')) %>% 
  ggplot(., aes(x=year, y=value, colour=model))+
  geom_line(size = 1) +
  facet_wrap(~tsa_industry_label_deu)