# ==============================
# 0. Setup and Packages
# ==============================

#set working directory

setwd("C:/Users/maria/Documents/R files")


#install packages

install.packages("data.table")
install.packages("lubridate")
install.packages("dplyr")
install.packages("tidyr")
install.packages("ggplot2")
install.packages("stringr")
install.packages("stringdist")
install.packages("readr")
install.packages("scales")
install.packages("plm")
install.packages("tseries")
install.packages("corrplot")
install.packages("mice")
install.packages("lmtest")
install.packages("sandwich")
install.packages("miceadds")
install.packages("tibble")
install.packages("broom")
install.packages("lmtest")
install.packages("fixest")
install.packages("urca")
install.packages("car")
install.packages("lfe")
install.packages("tidyverse")
#load packages

library(data.table)
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(stringdist)
library(readr)
library(scales)
library(plm)
library(tseries)
library(corrplot)
library(mice)
library(lmtest)
library(sandwich)
library(miceadds)
library(tibble)
library(broom)
library(lmtest)
library(fixest)
library(urca)
library(car)
library(lfe)
library(tidyverse)

# ==============================
# 1. Descriptive statistics and initial data cleaning
# ==============================


#load the data and rename columns

data <- fread("CRV_SVK_2024.csv", encoding="UTF-8")
mapping <- fread("variable_renaming_filled.csv",
                 stringsAsFactors = FALSE,
                 encoding = "UTF-8")
setnames(data, mapping$new_name)
print(names(data))
summary(data)

#correct the date format

class(data$date_first_registration)
head(data$date_first_registration, 10)
data <- data %>%
  mutate(
    date_first_registration    = as.Date(date_first_registration,
                                         format = "%d.%m.%Y"),
    date_first_registration_sk = as.Date(date_first_registration_sk,
                                         format = "%d.%m.%Y")
  )

# Find the earliest and latest registration date in SK

earliest_date_sr <- min(data$date_first_registration_sk, na.rm = TRUE)
print(earliest_date_sr)
last_date_sr <- max(data$date_first_registration_sk, na.rm = TRUE)
print(last_date_sr)

#Data inspection
# Count missing values per column

missing_counts <- sapply(data, function(x) {
  if (is.character(x)) {
    sum(is.na(x) | x == "")
  } else {
    sum(is.na(x))
  }
})
print(missing_counts)

# Calculate the percentage of missing values per column

missing_percent <- sapply(data, function(x) {
  if (is.character(x)) {
    (sum(is.na(x) | x == "") / length(x)) * 100
  } else {
    (sum(is.na(x)) / length(x)) * 100
  }
})
print(missing_percent)

#Filter passenger cars

data_cars <- data %>%
  filter(vehicle_type == "OSOBNÉ VOZIDLO")
passenger_count <- nrow(data_cars)
print(passenger_count)

#Handle missing or invalid dates

missing_only_first <- data_cars %>%
  filter(
    (is.na(date_first_registration) | date_first_registration == ""),
    !(is.na(date_first_registration_sk) | date_first_registration_sk == "")
  )
n_missing_only_first <- nrow(missing_only_first)
print(n_missing_only_first)
missing_dates <- data_cars %>%
  filter(
    is.na(date_first_registration) | date_first_registration == "" |
    is.na(date_first_registration_sk) | date_first_registration_sk == ""
  )
nrow(missing_dates)
invalid_dates <- data_cars %>%
  filter(date_first_registration_sk < date_first_registration)
n_invalid <- nrow(invalid_dates)
print(n_invalid)

data_cars <- data_cars %>%
  filter(
    !is.na(date_first_registration),                     
    !is.na(date_first_registration_sk),                  
    date_first_registration_sk >= date_first_registration
  )

#Summarize yearly registrations, distinguishing between new vs. used vehicles

yearly_passenger <- data_cars %>%
  mutate(
    year_reg = year(date_first_registration_sk),
    condition = ifelse(
      date_first_registration == date_first_registration_sk,
      "New",
      "Used"
    )
  ) %>%
  group_by(year_reg, condition) %>%
  summarise(count = n(), .groups = "drop")

#Create a wide table showing New/Used counts per year

yearly_table <- yearly_passenger %>%
  pivot_wider(
    names_from = condition,
    values_from = count,
    values_fill = 0
  )

# Print the table in the console
print(yearly_table)

# 4. Plot a bar chart showing yearly registrations, differentiated by New vs. Used
ggplot(yearly_passenger, aes(x = factor(year_reg), y = count, fill = condition)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Yearly Passenger Car Registrations",
    x = "Year",
    y = "Number of Registrations",
    fill = "Vehicle Type"
  ) +
  theme_minimal()

is_after_nov5 <- function(reg_date) {
  this_year <- year(reg_date)
  cutoff <- ymd(paste0(this_year, "-11-05"))
  reg_date > cutoff
}

historical_5years <- data_cars %>%
  mutate(
    Year = year(date_first_registration_sk),
    condition = ifelse(
      date_first_registration == date_first_registration_sk,
      "New",
      "Used"
    )
  ) %>%
  filter(Year >= 2019 & Year <= 2023)

hist_yearly <- historical_5years %>%
  group_by(Year, condition) %>%
  summarise(Total = n(), .groups = "drop")

hist_after_nov5 <- historical_5years %>%
  filter(is_after_nov5(date_first_registration_sk)) %>%
  group_by(Year, condition) %>%
  summarise(AfterNov5 = n(), .groups = "drop")

hist_joined <- hist_yearly %>%
  left_join(hist_after_nov5, by = c("Year", "condition")) %>%
  mutate(
    AfterNov5 = if_else(is.na(AfterNov5), 0L, AfterNov5),
    AfterNov5Share = AfterNov5 / Total
  )
average_share_after_nov5 <- hist_joined %>%
  group_by(condition) %>%
  summarise(AverageShare = mean(AfterNov5Share, na.rm = TRUE), .groups = "drop")

cutoff_2024 <- ymd("2024-11-05")

current_2024 <- data_cars %>%
  mutate(
    condition = ifelse(
      date_first_registration == date_first_registration_sk,
      "New",
      "Used"
    )
  ) %>%
  filter(
    year(date_first_registration_sk) == 2024,
    date_first_registration_sk <= cutoff_2024
  ) %>%
  group_by(condition) %>%
  summarise(SoFar = n(), .groups = "drop")

estimates_2024 <- current_2024 %>%
  left_join(average_share_after_nov5, by = "condition") %>%
  mutate(
    EstimatedTotal_2024 = SoFar / (1 - AverageShare),
    EstimatedAfterNov5_2024 = EstimatedTotal_2024 * AverageShare,
    Additional_AfterNov5 = EstimatedAfterNov5_2024 - SoFar
  )
print(estimates_2024)

fuel_counts <- data_cars %>%
  count(fuel_type_raw, sort = TRUE)
write.csv(fuel_counts, "fuel_counts_export.csv", row.names = FALSE)

data_cars <- data_cars %>%
  mutate(
    fuel_group = case_when(
      fuel_type_raw == "Nafta" ~ "Diesel",
      fuel_type_raw == "Benzín" ~ "Petrol",
      fuel_type_raw == "Benzín  + elektrický pohon nie plug-in" ~ "Petrol hybrid",
      fuel_type_raw == "Benzín + LPG" ~ "Other",
      fuel_type_raw == "Nafta  + elektrický pohon nie plug-in" ~ "Diesel hybrid",
      fuel_type_raw == "ELEKTRINA" ~ "Electric",
      fuel_type_raw == "Benzín  + elektrický pohon plug-in" ~ "Plug in hybrid petrol",
      fuel_type_raw == "BA 95 B + ELEKTRINA" ~ "Petrol hybrid",
      fuel_type_raw == "Benzín + CNG" ~ "Other",
      fuel_type_raw == "CNG" ~ "Other",
      fuel_type_raw == "Nafta + elektrický pohon plug-in" ~ "Plug in hybrid diesel",
      fuel_type_raw == "NM + ELEKTRINA" ~ "Diesel hybrid",
      fuel_type_raw == "Benzín  + LPG + elektrický pohon nie plug-in" ~ "Petrol hybrid",
      fuel_type_raw == "BA 95 B + LPG + ELEKTRINA" ~ "Petrol hybrid",
      fuel_type_raw == "BA 98 B + ELEKTRINA" ~ "Petrol hybrid",
      fuel_type_raw == "BA 95" ~ "Petrol",
      fuel_type_raw == "ELEKTRINA+Benzín" ~ "Petrol hybrid",
      fuel_type_raw == "Benzín prímesou oleja 1:40 (dvojtaktný motor)" ~ "Petrol",
      fuel_type_raw == "BENZÍN+ETANOL" ~ "Other",
      fuel_type_raw == "Benzín prímesou oleja 1:25 (dvojtaktný motor)" ~ "Petrol",
      fuel_type_raw == "Benzín  + LPG + elektrický pohon plug-in" ~ "Plug in hybrid petrol",
      fuel_type_raw == "BA 95 B" ~ "Petrol",
      fuel_type_raw == "Benzín + LNG" ~ "Other",
      fuel_type_raw == "Vodík" ~ "Other",
      fuel_type_raw == "BA 91" ~ "Petrol",
      fuel_type_raw == "BA 100 B + ELEKTRINA" ~ "Petrol hybrid",
      fuel_type_raw == "BA 96" ~ "Petrol",
      fuel_type_raw == "Benzín prímesou oleja 1:40 (dvojtaktný motor) + CNG" ~ "Other",
      fuel_type_raw == "ELEKTRINA+Benzín+LPG" ~ "Petrol hybrid",
      fuel_type_raw == "ETANOL + elektrický pohon plug-in" ~ "Other",
      fuel_type_raw == "LPG" ~ "Other",
      fuel_type_raw == "BA 90" ~ "Petrol",
      fuel_type_raw == "BA 98" ~ "Petrol",
      fuel_type_raw == "BIONAFTA" ~ "Diesel",
      fuel_type_raw == "Benzín prímesou oleja 1:40 (dvojtaktný motor) + LPG" ~ "Other",
      fuel_type_raw == "ELEKTRINA+Nafta" ~ "Diesel hybrid",
      fuel_type_raw == "ETANOL  + elektrický pohon nie plug-in" ~ "Other",
      fuel_type_raw == "ETANOL E85" ~ "Other",
      fuel_type_raw == "NM + NM BIO 48" ~ "Diesel",
      fuel_type_raw == "Nafta na báze parafínu" ~ "Diesel",
      TRUE ~ "Other"
    )
  )

data_cars <- data_cars %>%
  mutate(
    condition = ifelse(
      date_first_registration == date_first_registration_sk,
      "New",
      "Used"
    )
  )

yearly_fuel <- data_cars %>%
  mutate(
    year_reg = lubridate::year(date_first_registration_sk)
  ) %>%
  group_by(year_reg, condition, fuel_group) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(year_reg, condition) %>%
  mutate(
    total_year_type = sum(count),
    share = count / total_year_type * 100 
  ) %>%
  ungroup()

plot_fuel_trends <- function(data_cars, vehicle_type) {
  ggplot(
    data = data_cars %>% filter(condition == vehicle_type),
    aes(x = year_reg, y = share, color = fuel_group)
  ) +
    geom_line(size = 1) +
    scale_x_continuous(
      breaks = seq(2009, 2024, 1) 
    ) +
    scale_y_continuous(labels = scales::percent_format(scale = 1)) +
    labs(
      title = paste("Fuel Type Trends for", vehicle_type, "Cars"),
      x = "Year",
      y = "Share (%)",
      color = "Fuel Group"
    ) +
    theme_minimal(base_size=16) +
    theme(
      panel.grid.minor = element_blank(),
      legend.position = "bottom"
    )
}

plot_new  <- plot_fuel_trends(yearly_fuel, "New")
plot_used <- plot_fuel_trends(yearly_fuel, "Used")
print(plot_new)
print(plot_used)



data_cars <- data_cars %>%
  mutate(
    year_reg = year(date_first_registration_sk)
  )

monthly_data <- data_cars %>%
  filter(
    !is.na(date_first_registration),
    !is.na(date_first_registration_sk)
  ) %>%
  mutate(
    Month = floor_date(date_first_registration_sk, "month")
  ) %>%
  filter(!(year(Month) == 2024 & month(Month) == 11)) %>%
  group_by(Month, condition) %>%
  summarise(count = n(), .groups = "drop")

ggplot(monthly_data, aes(x = Month, y = count, color = condition)) +
  geom_line(size = 0.8) +
  
  geom_vline(xintercept = as.Date("2012-10-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2017-02-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2023-07-01"), color = "red", size = 1) +
  
  scale_color_manual(values = c("New" = "#f8766d", "Used" = "#00bfc4")) +
  
  scale_x_date(
    limits     = c(as.Date("2009-01-01"), as.Date("2025-01-01")),
    date_breaks = "1 year",
    date_labels = "%b %Y"
  ) +
  
  labs(
    title = "Monthly Passenger Car Registrations (New vs. Used)",
    x     = "Month",
    y     = "Number of Registrations",
    color = ""
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor        = element_blank(),
    legend.position         = "right",
    axis.text.x             = element_text(angle = 90, vjust = 0.5)
  )
monthly_power <- data_cars %>%
  mutate(
    Month = floor_date(date_first_registration_sk, "month")
  ) %>%
  group_by(Month) %>%
  summarise(
    AvgPower = mean(power_kw, na.rm = TRUE),
    .groups = "drop"
  )

monthly_co2 <- data_cars %>%
  mutate(
    Month = floor_date(date_first_registration_sk, "month")
  ) %>%
  group_by(Month) %>%
  summarise(
    AvgCO2 = mean(emission_co2, na.rm = TRUE),
    .groups = "drop"
  )

plot_power <- ggplot(monthly_power, aes(x = Month, y = AvgPower)) +
  geom_line(color = "blue", size = 1) +
  
 
  geom_vline(xintercept = as.Date("2012-10-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2017-02-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2023-07-01"), color = "red", size = 1) +
  
  
  scale_x_date(
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  
  labs(
    title = "Monthly Average Engine Power (kW)",
    x = "Year",
    y = "Average Engine Power (kW)"
  ) +
  theme_minimal(base_size=16)

plot_co2 <- ggplot(monthly_co2, aes(x = Month, y = AvgCO2)) +
  geom_line(color = "darkgreen", size = 1) +
  
  geom_vline(xintercept = as.Date("2012-10-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2017-02-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2023-07-01"), color = "red", size = 1) +
  
  scale_x_date(
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  
  labs(
    title = "Monthly Average CO2 Emissions (g CO2/km)",
    x = "Year",
    y = "Average CO2 (g CO2/km"
  ) +
  theme_minimal(base_size=16)

print(plot_power)
print(plot_co2)

monthly_stats <- data_cars %>%
  mutate(Month = floor_date(date_first_registration_sk, "month")) %>%
  group_by(Month) %>%
  summarise(
    AvgPower = mean(power_kw, na.rm = TRUE),
    AvgCO2   = mean(emission_co2, na.rm = TRUE),
    .groups  = "drop"
  )

monthly_stats_long <- monthly_stats %>%
  pivot_longer(
    cols = c("AvgPower", "AvgCO2"),
    names_to = "Measure", 
    values_to = "Value"
  )

ggplot(monthly_stats_long, aes(x = Month, y = Value, color = Measure)) +
  geom_line(size = 1) +
  
  geom_vline(xintercept = as.Date("2012-10-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2017-02-01"), color = "red", size = 1) +
  geom_vline(xintercept = as.Date("2023-07-01"), color = "red", size = 1) +
  
  scale_x_date(
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  
  scale_color_manual(
    values = c("AvgPower" = "blue", "AvgCO2" = "green"),
    labels = c("AvgPower" = "Power", "AvgCO2" = "CO2")
  ) +
  
  labs(
    title = "Monthly Averages: Power & CO2 Over Time",
    x = "Year",
    y = "Value (Different Units)",
    color = "Measure"
  ) +
  theme_minimal() +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

manufacturer_table <- data_cars %>%
  group_by(brand) %>%
  summarise(
    count = n(),        
    .groups = "drop"
  ) %>%
  mutate(
    percentage = 100 * count / sum(count)   
  ) %>%
  arrange(desc(count))

print(manufacturer_table)

write.csv(manufacturer_table, "manufacturer_table.csv", row.names = FALSE)
 

 numeric_cols <- names(data_cars)[sapply(data_cars, is.numeric)]
 
 outlier_table <- data.frame(
   Variable = character(),
   NonMissing = numeric(),
   OutlierCount = numeric(),
   OutlierPercent = numeric(),
   stringsAsFactors = FALSE
 )
 
 
 for (col_name in numeric_cols) {
   
   vec <- data_cars[[col_name]]
   vec_no_na <- vec[!is.na(vec)]
   
   if (length(vec_no_na) == 0) {
     next
   }
   
   m <- mean(vec_no_na)
   s <- sd(vec_no_na)
   
   lower_bound <- m - 2 * s
   upper_bound <- m + 2 * s
   
   outlier_count <- sum(vec_no_na < lower_bound | vec_no_na > upper_bound)
   
   outlier_table <- rbind(
     outlier_table,
     data.frame(
       Variable = col_name,
       NonMissing = length(vec_no_na),
       OutlierCount = outlier_count,
       OutlierPercent = round((outlier_count / length(vec_no_na)) * 100, 2),
       stringsAsFactors = FALSE
     )
   )
 }
 
 print(outlier_table)
 
 ggplot(yearly_passenger, aes(x = year_reg, y = count, fill = condition)) +
   geom_col(position = "stack") +
   
   scale_x_continuous(
     breaks = 2009:2024,
     limits = c(2009, 2024)
   ) +
   
   labs(
     title = "Passenger Car Registrations by Year",
     x = "Year",
     y = "Number of Registrations",
     fill = ""
   ) +
   theme_minimal() +
   theme(
     panel.grid.minor = element_blank(),
     legend.position = "right"
   )
 
 range(yearly_passenger$year_reg, na.rm = TRUE)
 
 unique(yearly_passenger$year_reg)
 
 yearly_passenger %>%
   filter(is.na(year_reg) | year_reg < 2009 | year_reg > 2024)
 
 yearly_passenger %>%
   filter(is.na(count) | count < 0)
 
 yearly_passenger <- data_cars %>%
   mutate(
     year_reg = year(date_first_registration_sk),
     condition = ifelse(
       date_first_registration == date_first_registration_sk,
       "New",
       "Used"
     )
   ) %>%
   group_by(year_reg, condition) %>%
   summarise(count = n(), .groups = "drop")

 yearly_table <- yearly_passenger %>%
   pivot_wider(
     names_from = condition,
     values_from = count,
     values_fill = 0
   )
 
 print(yearly_table)
 
 {
 year_totals <- yearly_passenger %>%
   group_by(year_reg) %>%
   summarise(total = sum(count), .groups = "drop")
 
 yearly_passenger <- yearly_passenger %>%
   mutate(opacity = ifelse(year_reg == 2024, 0.4, 1))
 
ggplot() +
  geom_bar(
    data = yearly_passenger,
    aes(
      x = year_reg,
      y = count,
      fill = condition,
      alpha = opacity
    ),
    stat = "identity",
    position = "dodge",
    show.legend = TRUE
  ) +
  geom_line(
    data = year_totals,
    aes(
      x = year_reg,
      y = total,
      group = 1,
      color = "Total"
    ),
    size = 1
  ) +
  geom_point(
    data = year_totals,
    aes(
      x = year_reg,
      y = total,
      group = 1,
      color = "Total"
    ),
    size = 2
  ) +
  scale_fill_discrete(name = "") +
  scale_alpha(range = c(0.4, 1), guide = "none") +
  scale_color_manual(
    name = "",
    values = c("Total" = "black"),
    breaks = c("Total")
  ) +
  scale_x_continuous(
    breaks = seq(min(year_totals$year_reg), max(year_totals$year_reg), by = 1)
  ) +
  labs(
    title = "Yearly Passenger Car Registrations (New vs. Used) + Total",
    subtitle = "Data for 2024 up to 4th November",
    x = "Year",
    y = "Number of Registrations"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  ) +
  guides(
    color = guide_legend(
      override.aes = list(
        linetype = 1,   
        shape = 16,     
        fill = NA,      
        size = 2
      )
    )
  ) +
  geom_vline(
    xintercept = 2023.5,
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "text",
    x = 2024,
    y = max(year_totals$total, na.rm = TRUE) * 1.05,
    label = "",
    color = "red",
    hjust = 1
  )

ggsave("yearly_reg.png", width = 16, height = 9, units = "in", dpi = 300)
 }
 year_totals <- yearly_passenger %>%
   group_by(year_reg) %>%
   summarise(total = sum(count), .groups = "drop")
 
 ggplot() +
   geom_bar(
     data = yearly_passenger,
     aes(
       x = factor(year_reg),
       y = count,
       fill = condition
     ),
     stat = "identity",
     position = "dodge"
   ) +
   
   geom_line(
     data = year_totals,
     aes(
       x = factor(year_reg),
       y = total,
       color = "Total",  
       group = 1
     ),
     size = 1
   ) +
   
   geom_point(
     data = year_totals,
     aes(
       x = factor(year_reg),
       y = total,
       color = "Total",  
       group = 1
     ),
     size = 2
   ) +
   
   scale_fill_discrete(name = "") +
   scale_color_manual(
     name = "",              
     values = c("Total" = "black"),
     breaks = c("Total")    
   ) +
   
   labs(
     title = "Yearly Passenger Car Registrations (New vs. Used) + Total",
     x = "Year",
     y = "Number of Registrations"
   ) +
   theme_minimal() +
   theme(
     panel.grid.minor = element_blank(),
     legend.position = "right"
   )
 
# ==============================
# 2. Fuel Data
# ==============================
 
 
 fuel_data <- read.csv("Fuel Data.csv")
 head(fuel_data)
 
 fuel_data <- fuel_data %>%
   arrange(Year, Month)
 
 cpi_data <- read.csv("CPI Data.csv")

 fuel_data <- fuel_data %>%
   left_join(cpi_data, by = c("Year", "Month"))
 
 head(fuel_data)
 
 fuel_data <- fuel_data %>%
   filter(
     Year > 2008 | (Year == 2009 & Month >= 1)
   )
 head(fuel_data)
 colnames(fuel_data)
 cols_to_num <- c(
   "Gasoline.95.octane..eur.1l.",
   "Gasoline.98.octane..eur.1l.",
   "LPG..eur.1l.",
   "Diesel.Oil..eur.1l.",
   "CNG..eur.1kg.",
   "LNG..eur.1kg.",
   "AdBlue..eur.10l.",
   "bioLNG..eur.1kg.",
   "Hydrogen..eur.1kg.",
   "Electrical.energy.AC.charging..eur.1kWh.",
   "Electrical.energy.DC.charging..eur.1kWh.",
   "Electrical.energy.ultra.fast.charging..eur.1kWh.",
   "X",
   "CPI"
 )
 
 fuel_data <- fuel_data %>%
   mutate(across(
     all_of(cols_to_num),
     ~ as.numeric(gsub(",", ".", .x))
   ))
 
 str(fuel_data[ , cols_to_num])
 
 fuel_data <- fuel_data %>%
   rename(
     year                        = Year,
     month                       = Month,
     petrol_95                   = Gasoline.95.octane..eur.1l.,
     petrol_98                   = Gasoline.98.octane..eur.1l.,
     lpg                         = LPG..eur.1l.,
     diesel                      = Diesel.Oil..eur.1l.,
     cng                         = CNG..eur.1kg.,
     lng                         = LNG..eur.1kg.,
     adbluel                     = AdBlue..eur.10l.,
     biolng                      = bioLNG..eur.1kg.,
     hydrogen                    = Hydrogen..eur.1kg.,
     ac_charging                 = Electrical.energy.AC.charging..eur.1kWh.,
     dc_charging                 = Electrical.energy.DC.charging..eur.1kWh.,
     ultrafast_charging          = Electrical.energy.ultra.fast.charging..eur.1kWh.,
     x                           = X,
     cpi                         = CPI
   )
 
summary(fuel_data)

fuel_data <- fuel_data %>%
  mutate(
    date = as.Date(paste0(year, "-", month, "-01")),
    real_petrol_95 = petrol_95 * (100 / cpi),
    real_diesel    = diesel    * (100 / cpi)
  )

ggplot(fuel_data, aes(x = date)) +
  geom_line(aes(y = real_petrol_95, color = "Petrol (Real)"), size = 1.2) +
  geom_line(aes(y = real_diesel,    color = "Diesel (Real)"), size = 1.2) +

  scale_color_manual(
    name   = NULL,
    values = c("Petrol (Real)" = "blue", "Diesel (Real)" = "red")
  ) +
  
  scale_x_date(
    breaks = seq(as.Date("2009-01-01"), as.Date("2025-01-01"), by = "1 years"),
    labels = date_format("%Y"),
    limits = c(as.Date("2009-01-01"), as.Date("2025-01-01"))
  ) +
  
  labs(
    title = "Real Petrol vs. Diesel Prices",
    x     = "Year",
    y     = "Price (€)"
  ) +
  
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor        = element_blank(),
    panel.grid.major.x      = element_line(color = "grey80"),
    panel.grid.major.y      = element_line(color = "grey90"),
    legend.position         = "bottom",
    legend.text             = element_text(size = 12),
    axis.text.x             = element_text(angle = 0, vjust = 0.5),
    axis.title              = element_text(face = "bold")
  )

# ==============================
# 3. Registration Fee + Model String Cleaning
# ==============================

unique(data_cars$emission_norm)
unique(data_cars$raw_euro)
unique(data_cars$emission_norm_raw)

data_cars <- data_cars %>%
  mutate(
    emission_norm = case_when(
      str_detect(emission_norm_raw, regex("EURO\\s?1", ignore_case = TRUE)) ~ "Euro 1",
      str_detect(emission_norm_raw, regex("EURO\\s?2", ignore_case = TRUE)) ~ "Euro 2",
      str_detect(emission_norm_raw, regex("EURO\\s?3", ignore_case = TRUE)) ~ "Euro 3",
      str_detect(emission_norm_raw, regex("EURO\\s?(IV|4)", ignore_case = TRUE)) ~ "Euro 4",
      str_detect(emission_norm_raw, regex("EURO\\s?5a", ignore_case = TRUE)) ~ "Euro 5a",
      str_detect(emission_norm_raw, regex("EURO\\s?5b", ignore_case = TRUE)) ~ "Euro 5b",
      str_detect(emission_norm_raw, regex("EURO\\s?5", ignore_case = TRUE)) ~ "Euro 5",
      str_detect(emission_norm_raw, regex("EURO\\s?6a", ignore_case = TRUE)) ~ "Euro 6a",
      str_detect(emission_norm_raw, regex("EURO\\s?6b", ignore_case = TRUE)) ~ "Euro 6b",
      str_detect(emission_norm_raw, regex("EURO\\s?6c", ignore_case = TRUE)) ~ "Euro 6c",
      str_detect(emission_norm_raw, regex("EURO\\s?6d", ignore_case = TRUE)) ~ "Euro 6d",
      str_detect(emission_norm_raw, regex("EURO\\s?6e", ignore_case = TRUE)) ~ "Euro 6e",
      str_detect(emission_norm_raw, regex("EURO\\s?6", ignore_case = TRUE)) ~ "Euro 6",
      str_detect(emission_norm_raw, regex("EEV", ignore_case = TRUE)) ~ "EEV",
      TRUE ~ NA_character_
    )
  )

unique(data_cars$emission_norm)
unique(data_cars$fuel_group)

data_cars <- data_cars %>%
  filter(date_first_registration_sk >= date_first_registration)

data_fee_period1 <- data_cars %>%
  filter(date_first_registration_sk >= as.Date("2009-01-01") &
           date_first_registration_sk <= as.Date("2012-09-30")) %>%
  mutate(fee_period = "2009_2012")

data_fee_period2 <- data_cars %>%
  filter(date_first_registration_sk >= as.Date("2012-10-01") &
           date_first_registration_sk <= as.Date("2017-01-31")) %>%
  mutate(fee_period = "2012_2017")

data_fee_period3 <- data_cars %>%
  filter(date_first_registration_sk >= as.Date("2017-02-01") &
           date_first_registration_sk <= as.Date("2023-06-30")) %>%
  mutate(fee_period = "2017_2023")

data_fee_period4 <- data_cars %>%
  filter(date_first_registration_sk >= as.Date("2023-07-01")) %>%
  mutate(fee_period = "2023_")

data_fee_period1 %>% summarise(min = min(date_first_registration_sk), max = max(date_first_registration_sk))
data_fee_period2 %>% summarise(min = min(date_first_registration_sk), max = max(date_first_registration_sk))
data_fee_period3 %>% summarise(min = min(date_first_registration_sk), max = max(date_first_registration_sk))
data_fee_period4 %>% summarise(min = min(date_first_registration_sk), max = max(date_first_registration_sk))

data_fee_period1 <- data_fee_period1 %>%
  mutate(reg_fee = 33)

data_fee_period2 <- data_fee_period2 %>%
  mutate(
    reg_fee = case_when(
      fuel_group == "Electric" ~ 33,
      
      power_kw <= 80 ~ 33,
      power_kw <= 86 ~ 167,
      power_kw <= 92 ~ 217,
      power_kw <= 98 ~ 267,
      power_kw <= 104 ~ 327,
      power_kw <= 110 ~ 397,
      power_kw <= 121 ~ 477,
      power_kw <= 132 ~ 657,
      power_kw <= 143 ~ 787,
      power_kw <= 154 ~ 957,
      power_kw <= 165 ~ 1157,
      power_kw <= 176 ~ 1397,
      power_kw <= 202 ~ 1697,
      power_kw <= 228 ~ 2047,
      power_kw <= 254 ~ 2467,
      power_kw > 254 ~ 2997,
      TRUE ~ NA_real_
    )
  )

summary(data_fee_period2$reg_fee)

unique(data_cars$fuel_group)
unique(data_cars$fuel_type_raw)

data_fee_period3 <- data_fee_period3 %>%
  mutate(
    years = as.integer(difftime(date_first_registration_sk, date_first_registration, units = "days") / 365.25),
    
    base_tariff = case_when(
      power_kw <= 80 ~ 33,
      power_kw <= 86 ~ 90,
      power_kw <= 92 ~ 110,
      power_kw <= 98 ~ 150,
      power_kw <= 104 ~ 210,
      power_kw <= 110 ~ 260,
      power_kw <= 121 ~ 360,
      power_kw <= 132 ~ 530,
      power_kw <= 143 ~ 700,
      power_kw <= 154 ~ 870,
      power_kw <= 165 ~ 1100,
      power_kw <= 176 ~ 1250,
      power_kw <= 202 ~ 1900,
      power_kw <= 228 ~ 2300,
      power_kw <= 254 ~ 2700,
      power_kw > 254 ~ 3900,
      TRUE ~ NA_real_
    ),
    
    age_coef = case_when(
      years == 0 ~ 1,
      years == 1 ~ 0.82,
      years == 2 ~ 0.68,
      years == 3 ~ 0.56,
      years == 4 ~ 0.46,
      years == 5 ~ 0.38,
      years == 6 ~ 0.32,
      years == 7 ~ 0.26,
      years == 8 ~ 0.23,
      years == 9 ~ 0.19,
      years == 10 ~ 0.16,
      years == 11 ~ 0.14,
      years == 12 ~ 0.12,
      years == 13 ~ 0.10,
      years == 14 ~ 0.09,
      years == 15 ~ 0.08,
      years == 16 ~ 0.07,
      years > 16 ~ 0.06,
      TRUE ~ NA_real_
    ),
    
    raw_fee = base_tariff * age_coef,
    
    reg_fee = case_when(
      fuel_group == "Electric" ~ 33,
      fuel_group %in% c("Petrol hybrid", "Diesel hybrid", "Plug in hybrid petrol", "Plug in hybrid diesel") |
        grepl("LPG|CNG|Vodík", fuel_type_raw, ignore.case = TRUE) ~ pmax(raw_fee * 0.5, 33),
      TRUE ~ pmax(raw_fee, 33)
    )
  )

summary(data_fee_period3$reg_fee)  
colnames(data_cars)
unique(data_cars$emission_norm)

data_fee_period4 <- data_fee_period4 %>%
  mutate(
    vehicle_age = as.integer(difftime(date_first_registration_sk, date_first_registration, units = "days") / 365.25),
    
    base_tariff = case_when(
      power_kw <= 80 ~ 33,
      power_kw <= 90 ~ 60,
      power_kw <= 100 ~ 90,
      power_kw <= 110 ~ 120,
      power_kw <= 125 ~ 200,
      power_kw <= 140 ~ 300,
      power_kw <= 155 ~ 500,
      power_kw <= 170 ~ 700,
      power_kw <= 210 ~ 900,
      power_kw > 210 ~ 1000,
      TRUE ~ NA_real_
    ),
    
    ekv_coef = case_when(
      vehicle_age >= 40 ~ 0.1,
      
      fuel_group %in% c("Plug in hybrid petrol", "Plug in hybrid diesel") |
        grepl("Vodík", fuel_type_raw, ignore.case = TRUE) ~ 0.2,
      
      emission_norm %in% c("Euro 1", "Euro I") ~ 1,
      emission_norm %in% c("Euro 2", "Euro II") ~ 0.8,
      emission_norm %in% c("Euro 3", "Euro III") ~ 0.7,
      emission_norm %in% c("Euro 4", "Euro IV") ~ 0.6,
      emission_norm %in% c("Euro 5", "Euro V", "EEV") ~ 0.5,
      emission_norm %in% c("Euro 6a", "Euro 6b", "Euro 6c", "Euro VIA", "Euro VIB", "Euro VIC") ~ 0.45,
      emission_norm %in% c("Euro 6d", "Euro VID", "Euro 6e", "Euro VIE") ~ 0.4,
      TRUE ~ NA_real_
    ),
    
    ekv_coef = case_when(
      !is.na(ekv_coef) ~ ekv_coef,  # keep known values
      date_first_registration <= as.Date("1996-12-31") ~ 1,
      date_first_registration <= as.Date("2001-12-31") ~ 0.8,
      date_first_registration <= as.Date("2006-12-31") ~ 0.7,
      date_first_registration <= as.Date("2011-12-31") ~ 0.6,
      date_first_registration <= as.Date("2016-08-31") ~ 0.5,
      date_first_registration <= as.Date("2020-08-31") ~ 0.45,
      date_first_registration >  as.Date("2020-08-31") ~ 0.4,
      TRUE ~ NA_real_
    ),
    
    raw_fee = base_tariff * ekv_coef,
    
    reg_fee = case_when(
      fuel_group == "Electric" ~ 33,
      TRUE ~ pmax(raw_fee, 33)
    )
  )

summary(data_fee_period4$reg_fee)

data_fee_period1 <- data_fee_period1 %>%
  mutate(fee_period = 1)

data_fee_period2 <- data_fee_period2 %>%
  mutate(fee_period = 2)

data_fee_period3 <- data_fee_period3 %>%
  mutate(fee_period = 3)

data_fee_period4 <- data_fee_period4 %>%
  mutate(fee_period = 4)

data_cars <- bind_rows(
  data_fee_period1,
  data_fee_period2,
  data_fee_period3,
  data_fee_period4
)

summary(data_cars$reg_fee)

data_cars %>%
  filter(is.na(reg_fee)) %>%
  select(date_first_registration, power_kw, fuel_group, fuel_type_raw, emission_norm, fee_period) %>%
  View()
data_cars <- data_cars %>%
  filter(!is.na(reg_fee))

data_cars_trimmed <- data_cars %>%
  filter(date_first_registration_sk <= as.Date("2024-10-31"))

monthly_avg_fee <- data_cars_trimmed %>%
  mutate(month = floor_date(date_first_registration_sk, "month")) %>%
  group_by(month) %>%
  summarise(avg_fee = mean(reg_fee, na.rm = TRUE)) %>%
  ungroup()

policy_change_dates <- as.Date(c("2012-10-01", "2017-02-01", "2023-07-01"))

ggplot(monthly_avg_fee, aes(x = month, y = avg_fee)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_vline(xintercept = policy_change_dates, color = "red", linewidth = 1) +
  scale_x_date(
    breaks = seq(as.Date("2009-01-01"), as.Date("2024-10-31"), by = "1 year"),
    labels = scales::date_format("%Y"),
    limits = as.Date(c("2009-01-01", "2024-10-31"))
  ) +
  labs(
    title = "Monthly Average Registration Fee (€)",
    x = "Year",
    y = "Average Registration Fee (€)"
  ) +
  theme_minimal(base_size = 16)

data_cars_trimmed <- data_cars %>%
  filter(date_first_registration_sk <= as.Date("2024-10-31"))

monthly_revenue <- data_cars_trimmed %>%
  mutate(month = floor_date(date_first_registration_sk, "month")) %>%
  group_by(month) %>%
  summarise(total_revenue = sum(reg_fee, na.rm = TRUE)) %>%
  ungroup()

policy_change_dates <- as.Date(c("2012-10-01", "2017-02-01", "2023-07-01"))

ggplot(monthly_revenue, aes(x = month, y = total_revenue)) +
  geom_line(color = "darkgreen", linewidth = 1) +
  geom_vline(xintercept = policy_change_dates, color = "red", linewidth = 1) +
  scale_x_date(
    breaks = seq(as.Date("2009-01-01"), as.Date("2024-10-31"), by = "1 year"),
    labels = date_format("%Y"),
    limits = as.Date(c("2009-01-01", "2024-10-31"))
  ) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Monthly Revenue from Registration Fees (€)",
    x = "Year",
    y = "Total Monthly Revenue (€)"
  ) +
  theme_minimal(base_size = 16)

# ==============================
# 4. Cleaning Model Strings
# ==============================


length(unique(data_cars$model))

data_cars %>%
  count(model, sort = TRUE) %>%
  head(100) %>%
  View()

grouped_models <- data_cars %>%
  mutate(
    brand_clean = str_to_upper(brand),
    model_clean = str_to_upper(model),
    model_clean = str_replace_all(model_clean, "[^A-Z0-9 ]", " "),      
    model_clean = str_squish(model_clean),                              
    model_clean = str_remove_all(model_clean, fixed(brand_clean)),     
    model_clean = str_squish(model_clean),                              
    grouped_model = word(model_clean, 1)                                
  ) %>%
  select(brand, model, grouped_model) %>%
  distinct() %>%
  arrange(brand, grouped_model)

write_csv(grouped_models, "grouped_models_for_review.csv")

data_cars_test <- data_cars

model_proxy_test <- data_cars_test %>%
  group_by(power_kw, length, width, height, fuel_group, body_type, seat_count) %>%
  summarise(
    model_count = n_distinct(model),
    example_models = paste(unique(model)[1:min(3, length(unique(model)))], collapse = ", "),
    n = n()
  ) %>%
  arrange(model_count, desc(n))

View(model_proxy_test)

mean(model_proxy_test$model_count == 1)

data_cars_test <- data_cars %>%
  mutate(
    model_clean = str_to_upper(model),                                      
    model_clean = str_replace_all(model_clean, "[^A-Z0-9 ]", " "),           
    model_clean = str_squish(model_clean)                                    
  )

model_proxy_test <- data_cars_test %>%
  group_by(power_kw, length, width, height, fuel_group, body_type, seat_count) %>%
  summarise(model_count = n_distinct(model_clean), .groups = "drop")

mean(model_proxy_test$model_count == 1)

problematic_specs <- model_proxy_test %>%
  filter(model_count > 1)

problematic_rows <- data_cars_test %>%
  semi_join(problematic_specs, by = c("power_kw", "length", "width", "height", "fuel_group", "body_type", "seat_count"))

model_wide <- problematic_rows %>%
  select(power_kw, length, width, height, fuel_group, body_type, seat_count, model_clean) %>%
  distinct() %>%
  group_by(power_kw, length, width, height, fuel_group, body_type, seat_count) %>%
  mutate(model_id = paste0("model_", row_number())) %>%
  pivot_wider(names_from = model_id, values_from = model_clean)

write.csv(model_wide, "model_conflict_groups.csv", row.names = FALSE)

model_grouped <- read_csv("model_conflict_groups_with_grouped.csv") %>%
  select(power_kw, length, width, height, fuel_group, body_type, seat_count, grouped_model)

data_cars_test <- data_cars_test %>%
  left_join(model_grouped,
            by = c("power_kw", "length", "width", "height", "fuel_group", "body_type", "seat_count"))

data_cars_test <- data_cars_test %>%
  mutate(
    final_model = case_when(
      !is.na(grouped_model) ~ grouped_model,
      TRUE ~ model_clean
    )
  )

head(data_cars_test$grouped_model)
head(data_cars_test$final_model,100)

colnames(data_cars_test)

model_export <- data_cars_test %>%
  distinct(brand, final_model) %>%
  arrange(brand, final_model) %>%
  mutate(manual_check = "") 

write_csv(model_export, "final_model_inspect.csv")

data_cars_test <- data_cars_test %>%
  rename(model_raw = model)

model_manual <- read_csv("final_model_inspection.csv")


data_cars_test <- data_cars_test %>%
  left_join(
    model_manual %>% select(brand, final_model, manual_check),
    by = c("brand", "final_model")
  )

data_cars_test <- data_cars_test %>%
  mutate(model = manual_check)

data_cars_test <- data_cars_test %>%
  select(-manual_check)

head(data_cars_test$model,100)

vin_4times <- data_cars %>%
  group_by(vin) %>%
  filter(n() >= 5) %>%
  ungroup()

vin_4times_selected <- vin_4times %>%
  arrange(vin) %>%
  select(vin, date_first_registration, date_first_registration_sk, everything())

write_csv(vin_4times_selected, "vin_4times_selected.csv")


vin_counts <- data_cars_test %>%
  group_by(vin) %>%
  summarise(count = n(), .groups = "drop") %>%
  filter(count >= 5)

vin_multiple_entries <- data_cars_test %>%
  filter(vin %in% vin_counts$vin) %>%
  select(vin, date_first_registration, date_first_registration_sk, everything()) %>%
  arrange(desc(vin))  

write.csv(vin_multiple_entries, "vin_multiple_entries.csv")


earliest_vin_records <- data_cars_test %>%
  arrange(vin, date_first_registration_sk) %>%
  group_by(vin) %>%
  slice(1) %>%
  ungroup()

n_distinct(earliest_vin_records$vin) == nrow(earliest_vin_records)

data_cars_test <- earliest_vin_records


# ==============================
# 5. Panel Structure Setup + Models
# ==============================

#TRANSFORM TO PANEL STRUCTURE

data_cars_test <- data_cars_test %>%
  select(-final_model)

data_panel <- data_cars_test %>%
  mutate(
    power_bin = cut(power_kw,
                    breaks = seq(0, max(power_kw, na.rm = TRUE) + 10, by = 10),
                    right = FALSE, include.lowest = TRUE)
  )

data_cars_test <- data_cars_test %>%
  mutate(
    time = floor_date(date_first_registration_sk, "month")
  )

head(data_cars_test$time)

#Create a panel

panel_data <- data_cars_test %>%
  mutate(
    time = floor_date(date_first_registration_sk, "month")
  ) %>%
  group_by(time, brand, model, fuel_group, power_kw, condition) %>%
  summarise(
    regs_no      = n(),
    reg_fee      = mean(reg_fee, na.rm = TRUE),
    co2          = mean(emission_co2, na.rm = TRUE),
    consumption  = mean(fuel_consumption_combined, na.rm = TRUE),
    seats        = median(seat_count, na.rm = TRUE),
    weight       = mean(weight_gross_max, na.rm = TRUE),
    displacement = mean(engine_stroke, na.rm = TRUE),
    .groups      = "drop"
  ) %>%
  mutate(
    distinct_car_config = paste(brand, model, fuel_group, power_kw, condition, sep = "_")
  )

duplicate_check <- panel_data %>%
  group_by(time, distinct_car_config) %>%
  summarise(duplicates = n() > 1) %>%
  filter(duplicates == TRUE)

if (nrow(duplicate_check) > 0) {
  print("Duplicates found!")
} else {
  print("No duplicates!")
}

missing_check <- panel_data %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  gather(key = "variable", value = "missing_count") %>%
  filter(missing_count > 0)

print(missing_check)

data_cars_test %>%
  summarise(
    missing_power_kw     = sum(is.na(power_kw)),
    missing_co2          = sum(is.na(emission_co2)),
    missing_consumption  = sum(is.na(fuel_consumption_combined)),
    missing_displacement = sum(is.na(engine_stroke))
  )


group_time_dist <- panel_data %>%
  count(distinct_car_config) %>%
  count(n)

print(group_time_dist, n=190)

panel_data <- panel_data %>%
  mutate(
    year = year(time),
    month = month(time)
  )
panel_data <- panel_data %>%
  left_join(
    fuel_data %>% select(year, month, cpi),
    by = c("year", "month")
  )

panel_data <- panel_data %>%
  mutate(reg_fee_real = reg_fee * (100 / cpi))

unique(panel_data$fuel_group)

colnames(fuel_data)

panel_data <- panel_data %>%
  left_join(
    fuel_data %>%
      select(year, month, real_petrol_95, real_diesel, ac_charging),
    by = c("year", "month")
  ) %>%
  mutate(
    real_fuel_price = case_when(
      fuel_group %in% c("Petrol", "Petrol hybrid", "Plug in hybrid petrol") ~ real_petrol_95,
      fuel_group %in% c("Diesel", "Diesel hybrid", "Plug in hybrid diesel") ~ real_diesel,
      fuel_group == "Electric" ~ ac_charging,
      TRUE ~ NA_real_
    )
  ) %>%
  select(-real_petrol_95, -real_diesel, -ac_charging, -year, -month)

panel_data <- panel_data %>%
  mutate(
    cost_per_100km = real_fuel_price * consumption
  )

panel_plm <- pdata.frame(
  panel_data,
  index = c("distinct_car_config", "time"),
  drop.index = FALSE,
  row.names = TRUE
)


panel_plm$panel_rowid <- seq_len(nrow(panel_plm)) 

index(panel_plm)
summary(index(panel_plm))
pdim(panel_plm)


fuel_data <- fuel_data %>%
  mutate(mean_petrol_diesel = (real_petrol_95 + real_diesel) / 2)

panel_plm <- panel_plm %>%
  mutate(time = as.Date(as.character(time)))

panel_plm <- panel_plm %>%
  left_join(fuel_data %>% select(date, mean_petrol_diesel),
            by = c("time" = "date"))

panel_plm <- panel_plm %>%
  mutate(
    consumption = ifelse(fuel_group %in% c("Electric"), 0, consumption),
    real_fuel_price = ifelse(fuel_group %in% c("Electric", "Other"), mean_petrol_diesel, real_fuel_price)
  )

panel_plm %>%
  summarise(
    regs_no_NA = sum(is.na(regs_no)),
    regs_no_0  = sum(regs_no == 0, na.rm = TRUE),
    
    reg_fee_real_NA = sum(is.na(reg_fee_real)),
    reg_fee_real_0  = sum(reg_fee_real == 0, na.rm = TRUE),
    
    real_fuel_price_NA = sum(is.na(real_fuel_price)),
    real_fuel_price_0  = sum(real_fuel_price == 0, na.rm = TRUE),
    
    consumption_NA = sum(is.na(consumption)),
    consumption_0  = sum(consumption == 0, na.rm = TRUE),
    
    co2_NA = sum(is.na(co2)),
    co2_0  = sum(co2 == 0, na.rm = TRUE),
    
    power_kw_NA = sum(is.na(power_kw)),
    power_kw_0  = sum(power_kw == 0, na.rm = TRUE),
    
    seats_NA = sum(is.na(seats)),
    
    weight_NA = sum(is.na(weight)),
    weight_0  = sum(weight == 0, na.rm = TRUE)
  )

panel_plm <- panel_plm %>%
  mutate(
    consumption = ifelse(fuel_group %in% c("Electric"), 0.0001, consumption),
    co2         = ifelse(fuel_group %in% c("Electric"), 0.0001, co2)
  )
# === A) Overall Summary (Missing and Zeros) ===

vars_to_check <- c("co2", "consumption", "power_kw", "weight", "displacement")

overall_summary <- lapply(vars_to_check, function(var) {
  x <- panel_plm[[var]]
  data.frame(
    Variable = var,
    Total_N = length(x),
    Missing_Count = sum(is.na(x)),
    Missing_Percent = round(100 * sum(is.na(x)) / length(x), 2),
    Zero_Count = sum(x == 0, na.rm = TRUE),
    Zero_Percent = round(100 * sum(x == 0, na.rm = TRUE) / length(x), 2)
  )
}) %>% bind_rows()

print(overall_summary)
write.csv(overall_summary,"missingness.csv")


# List of variables to check
vars_to_check <- c("co2", "consumption", "displacement", "power_kw", "weight")

# Create a breakdown of missing and zero shares *by fuel group as a share of total missing/zeros*
fuel_issue_breakdown <- lapply(vars_to_check, function(var) {
  
  # Total missing and zero counts for the variable
  total_missing <- sum(is.na(panel_plm[[var]]))
  total_zeros <- sum(panel_plm[[var]] == 0, na.rm = TRUE)
  
  # Group by fuel and count missing/zero
  panel_plm %>%
    group_by(fuel_group) %>%
    summarise(
      Missing_Count = sum(is.na(.data[[var]])),
      Zero_Count = sum(.data[[var]] == 0, na.rm = TRUE)
    ) %>%
    mutate(
      Variable = var,
      Missing_Percent = round(100 * Missing_Count / total_missing, 2),
      Zero_Percent = round(100 * Zero_Count / total_zeros, 2)
    )
}) %>%
  bind_rows() %>%
  select(Variable, fuel_group, Missing_Count, Missing_Percent, Zero_Count, Zero_Percent) %>%
  arrange(Variable, desc(Missing_Percent))

print(fuel_issue_breakdown, n = 100)

# List of variables to check
vars_to_check <- c("co2", "consumption", "displacement", "power_kw", "weight")

# Build breakdown of missing, zero, and combined problematic values
fuel_issue_breakdown <- lapply(vars_to_check, function(var) {
  
  total_missing <- sum(is.na(panel_plm[[var]]))
  total_zeros <- sum(panel_plm[[var]] == 0, na.rm = TRUE)
  total_problematic <- sum(is.na(panel_plm[[var]]) | panel_plm[[var]] == 0, na.rm = TRUE)
  
  panel_plm %>%
    group_by(fuel_group) %>%
    summarise(
      Missing_Count = sum(is.na(.data[[var]])),
      Zero_Count    = sum(.data[[var]] == 0, na.rm = TRUE),
      Problem_Count = sum(is.na(.data[[var]]) | .data[[var]] == 0, na.rm = TRUE)
    ) %>%
    mutate(
      Variable           = var,
      Missing_Percent    = round(100 * Missing_Count / total_missing, 2),
      Zero_Percent       = round(100 * Zero_Count / total_zeros, 2),
      Problem_Percent    = round(100 * Problem_Count / total_problematic, 2)
    )
}) %>%
  bind_rows() %>%
  select(Variable, fuel_group,
         Problem_Count, Problem_Percent) %>%
  arrange(Variable, desc(Problem_Percent))

print(fuel_issue_breakdown, n = 100)

# Filter only for the relevant variables
heatmap_data <- fuel_issue_breakdown %>%
  filter(Variable %in% c("co2", "consumption"))

heatmap_data$fuel_group <- factor(
  heatmap_data$fuel_group,
  levels = heatmap_data %>%
    group_by(fuel_group) %>%
    summarise(total = sum(Problem_Count)) %>%
    arrange(desc(total)) %>%
    pull(fuel_group)
)

filtered_data <- fuel_issue_breakdown %>%
  filter(Variable %in% c("co2", "consumption"))

ggplot(filtered_data, aes(x = Problem_Percent, 
                          y = reorder(fuel_group, Problem_Percent),
                          fill = Problem_Percent)) +
  geom_col(width = 0.7) +
  geom_text(aes(label = paste0(round(Problem_Percent, 1), "%")),
            hjust = -0.1, size = 3.5) +
  scale_fill_gradient(low = "#fddbc7", high = "#b2182b") +
  facet_wrap(~ Variable, scales = "free_y", ncol = 1) +
  labs(
    title = "Problematic Values (NA and 0) Share by Fuel Group in CO2 and consumption",
    x = "Problematic Percentage",
    y = "Fuel Group"
  ) +
  coord_cartesian(xlim = c(0, max(filtered_data$Problem_Percent, na.rm = TRUE) + 5)) +
  theme_minimal(base_size = 16) +
  theme(
    legend.position = "none",
    strip.text = element_text(face = "bold", size = 12),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    plot.title = element_text(hjust = 0.5, size = 13)
  )


total_obs <- panel_plm %>%
  filter(!is.na(fuel_group)) %>%
  nrow()

fuel_counts <- panel_plm %>%
  filter(!is.na(fuel_group)) %>%
  group_by(fuel_group) %>%
  summarise(
    Count = n(),
    Percent = round(100 * Count / total_obs, 2)
  ) %>%
  arrange(desc(Count))

print(fuel_counts)

panel_plm %>%
  filter(fuel_group %in% c("Diesel", "Petrol")) %>%
  group_by(fuel_group) %>%
  summarise(
    missing_co2 = sum(is.na(co2)),
    missing_consumption = sum(is.na(consumption))
  )

panel_plm <- panel_plm %>% filter(!is.na(power_kw))


# Common vars used in both imputations
vars_for_mice <- c("co2", "consumption", "reg_fee_real", "real_fuel_price",
                   "power_kw", "seats", "weight")


# Subset Diesel data
diesel_data <- panel_plm %>%
  filter(fuel_group == "Diesel") %>%
  select(co2, consumption, reg_fee_real, real_fuel_price, power_kw, seats, weight)

# Summaries
summary_stats <- diesel_data %>%
  summarise(across(everything(), list(
    min = ~min(. , na.rm = TRUE),
    q1 = ~quantile(. , 0.25, na.rm = TRUE),
    median = ~median(. , na.rm = TRUE),
    mean = ~mean(. , na.rm = TRUE),
    q3 = ~quantile(. , 0.75, na.rm = TRUE),
    max = ~max(. , na.rm = TRUE)
  ), .names = "{.col}_{.fn}"))

print(summary_stats)

# Boxplots
par(mfrow = c(2, 3))  # Layout
boxplot(diesel_data$co2, main = "CO2")
boxplot(diesel_data$consumption, main = "Consumption")
boxplot(diesel_data$reg_fee_real, main = "Reg Fee")
boxplot(diesel_data$real_fuel_price, main = "Fuel Price")
boxplot(diesel_data$power_kw, main = "Power (kW)")
boxplot(diesel_data$weight, main = "Weight")

# ========== Diesel ==========
# List of variables to clean
vars_to_clean <- c("co2", "reg_fee_real", "real_fuel_price", "power_kw", "seats", "weight")

# ==============================
# Diesel: Prepare Data for MICE
# ==============================

# 1. Filter only Diesel vehicles and select relevant variables
diesel_data <- panel_plm %>%
  filter(fuel_group == "Diesel") %>%
  select(all_of(vars_for_mice))

# 2. Remove outliers
diesel_data <- diesel_data %>%
  filter(
    is.na(co2) | co2 <= 500,
    is.na(consumption) | consumption <= 30,
    is.na(power_kw) | power_kw <= 400,
    is.na(weight) | weight <= 4000
  )

# Run MICE imputation
diesel_mice <- mice(diesel_data, method = "pmm", m = 5, maxit = 5, seed = 123)

# ==============================
# Diagnostics: diesel_mice
# ==============================

#
# Density plots: Imputed vs observed
densityplot(diesel_mice, ~co2,
            main = "MICE imputation - Density Plot of CO2 – Diesel Vehicles",
            xlab = "CO2 (g/km)",
            ylab = "Density")

densityplot(diesel_mice, ~consumption,
            main = "MICE imputation - Density Plot of consumption – Diesel Vehicles",
            xlab = "consumption (l/100km)",
            ylab = "Density")

# Convergence check
plot(diesel_mice)

# Inspect imputation structure
print(diesel_mice$method)
print(diesel_mice$where)
summary(diesel_mice)


# ==============================
# Summary + Boxplots for Petrol
# ==============================

# Subset Petrol data
petrol_data <- panel_plm %>%
  filter(fuel_group == "Petrol") %>%
  select(co2, consumption, reg_fee_real, real_fuel_price, power_kw, seats, weight)

# Summaries
summary_stats_petrol <- petrol_data %>%
  summarise(across(everything(), list(
    min = ~min(. , na.rm = TRUE),
    q1 = ~quantile(. , 0.25, na.rm = TRUE),
    median = ~median(. , na.rm = TRUE),
    mean = ~mean(. , na.rm = TRUE),
    q3 = ~quantile(. , 0.75, na.rm = TRUE),
    max = ~max(. , na.rm = TRUE)
  ), .names = "{.col}_{.fn}"))

print(summary_stats_petrol)

# Boxplots
par(mfrow = c(2, 3))  # Layout for 6 plots
boxplot(petrol_data$co2, main = "Petrol: CO2")
boxplot(petrol_data$consumption, main = "Petrol: Consumption")
boxplot(petrol_data$reg_fee_real, main = "Petrol: Reg Fee")
boxplot(petrol_data$real_fuel_price, main = "Petrol: Fuel Price")
boxplot(petrol_data$power_kw, main = "Petrol: Power (kW)")
boxplot(petrol_data$weight, main = "Petrol: Weight")

# ==============================
# Petrol: Prepare Data for MICE
# ==============================

# 1. Filter only Petrol vehicles and select relevant variables
petrol_data <- panel_plm %>%
  filter(fuel_group == "Petrol") %>%
  select(all_of(vars_for_mice))

# 2. Remove outliers
petrol_data <- petrol_data %>%
  filter(
    is.na(co2) | co2 <= 400,
    is.na(consumption) | consumption <= 30,
    is.na(power_kw) | power_kw <= 400,
    is.na(weight) | weight <= 4000
  )

# Run MICE imputation
petrol_mice <- mice(petrol_data, method = "pmm", m = 5, maxit = 5, seed = 123)

# ==============================
# Diagnostics: petrol_mice
# ==============================

# Density plots: Imputed vs observed
densityplot(petrol_mice, ~co2,
            main = "MICE imputation - Density Plot of CO2 – Petrol Vehicles",
            xlab = "CO2 (g/km)",
            ylab = "Density")

densityplot(petrol_mice, ~consumption,
            main = "MICE imputation - Density Plot of consumption – Petrol Vehicles",
            xlab = "consumption (l/100km)",
            ylab = "Density")


# Convergence check
plot(petrol_mice)


plot(petrol_mice, main = "MICE Convergence Diagnostics: Petrol Vehicles")

plot(petrol_mice)
title("MICE Convergence Diagnostics: Petrol Vehicles", line = -1)

# Inspect imputation structure
print(petrol_mice$method)
print(petrol_mice$where)
summary(petrol_mice)
# === Prepare clean subsets for Diesel and Petrol (with panel_rowid) ===
diesel_subset_clean <- panel_plm %>%
  filter(fuel_group == "Diesel") %>%
  select(panel_rowid, all_of(vars_for_mice)) %>%
  filter(
    is.na(co2) | co2 <= 500,
    is.na(consumption) | consumption <= 30,
    is.na(power_kw) | power_kw <= 400,
    is.na(weight) | weight <= 4000
  )

petrol_subset_clean <- panel_plm %>%
  filter(fuel_group == "Petrol") %>%
  select(panel_rowid, all_of(vars_for_mice)) %>%
  filter(
    is.na(co2) | co2 <= 400,
    is.na(consumption) | consumption <= 30,
    is.na(power_kw) | power_kw <= 400,
    is.na(weight) | weight <= 4000
  )

# === Get all 5 imputed versions from MICE ===
diesel_imputations <- complete(diesel_mice, action = "all")
petrol_imputations <- complete(petrol_mice, action = "all")

# === Define injection function (no dummy, no *_star) ===
inject_imputations <- function(base_df, diesel_imp, petrol_imp) {
  base_df <- as.data.frame(base_df)
  base_df$panel_rowid <- as.integer(base_df$panel_rowid)
  diesel_imp$panel_rowid <- as.integer(diesel_imp$panel_rowid)
  petrol_imp$panel_rowid <- as.integer(petrol_imp$panel_rowid)
  
  base_df %>%
    left_join(diesel_imp %>%
                select(panel_rowid, co2_d = co2, consumption_d = consumption),
              by = "panel_rowid") %>%
    left_join(petrol_imp %>%
                select(panel_rowid, co2_p = co2, consumption_p = consumption),
              by = "panel_rowid") %>%
    mutate(
      co2 = case_when(
        fuel_group == "Diesel" & (is.na(co2) | co2 == 0) ~ co2_d,
        fuel_group == "Petrol" & (is.na(co2) | co2 == 0) ~ co2_p,
        TRUE ~ co2
      ),
      consumption = case_when(
        fuel_group == "Diesel" & (is.na(consumption) | consumption == 0) ~ consumption_d,
        fuel_group == "Petrol" & (is.na(consumption) | consumption == 0) ~ consumption_p,
        TRUE ~ consumption
      )
    ) %>%
    select(-co2_d, -co2_p, -consumption_d, -consumption_p)
}

# === Inject imputations to create 5 completed datasets ===
imputed_panels <- lapply(1:5, function(i) {
  diesel_completed <- bind_cols(diesel_subset_clean["panel_rowid"], diesel_imputations[[i]])
  petrol_completed <- bind_cols(petrol_subset_clean["panel_rowid"], petrol_imputations[[i]])
  inject_imputations(panel_plm, diesel_completed, petrol_completed)
})


# Define function to compute summary stats
get_summary_stats <- function(df, label) {
  df %>%
    summarise(
      mean_co2 = mean(co2, na.rm = TRUE),
      q1_co2 = quantile(co2, 0.25, na.rm = TRUE),
      median_co2 = median(co2, na.rm = TRUE),
      q3_co2 = quantile(co2, 0.75, na.rm = TRUE),
      mean_cons = mean(consumption, na.rm = TRUE),
      q1_cons = quantile(consumption, 0.25, na.rm = TRUE),
      median_cons = median(consumption, na.rm = TRUE),
      q3_cons = quantile(consumption, 0.75, na.rm = TRUE)
    ) %>%
    mutate(Source = label)
}

# Get stats for original (incomplete) data
original_stats <- get_summary_stats(panel_plm, "Original")

# Get stats for each of the 5 imputations
imputed_stats <- lapply(1:5, function(i) {
  get_summary_stats(imputed_panels[[i]], paste0("Imputed_", i))
})

# Combine into one data frame
summary_comparison <- bind_rows(original_stats, do.call(rbind, imputed_stats)) %>%
  select(Source, everything())

print(summary_comparison)

write.csv(summary_comparison,"mice_comparison.csv")

imputed_panels <- lapply(imputed_panels, function(df) {
  pdata.frame(df, index = c("distinct_car_config", "time"))
})

str(imputed_panels[[1]])



imputed_panels <- lapply(imputed_panels, function(df) {
  names(df)[names(df) == "time"] <- "evid_svk"
  
  df$lnreg         <- log(df$regs_no)
  df$lnrfee        <- log(df$reg_fee_real)
  df$lnfuelprice   <- log(df$real_fuel_price)
  df$lnconsumption <- log(df$consumption)
  df$lnco2         <- log(df$co2)
  df$lnkw          <- log(df$power_kw)
  df$lnweight      <- log(df$weight)
  df$lndisplacement<- log(df$displacement)
  df
})

imputed_panels <- lapply(imputed_panels, function(df) {
  df <- df[complete.cases(df[, c("lnreg", "lnrfee", "lnfuelprice", "lnconsumption", "lnco2", "lnkw", "seats", "lnweight")]), ]
  df
})
df <- imputed_panels[[1]]
idx <- complete.cases(df[, c("lnco2", "lnconsumption")])
finite_idx <- idx & is.finite(df$lnco2) & is.finite(df$lnconsumption)
cor(df$lnco2[finite_idx], df$lnconsumption[finite_idx])

vars <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "seats", "lnweight")

df <- imputed_panels[[1]]
cor_df <- df[, vars]
cor_df <- cor_df[complete.cases(cor_df) & 
                   apply(cor_df, 1, function(x) all(is.finite(x))), ]

# Compute correlation matrix
cor_mat <- cor(cor_df)

print(round(cor_mat, 3))


dev.new(width = 6, height = 6)  
corrplot(cor_mat, method = "color", addCoef.col = "black", number.cex = 1.1)
imputed_panels <- lapply(imputed_panels, function(df) {
  df$evid_svk <- as.Date(df$evid_svk)  
  df <- df[df$evid_svk <= as.Date("2024-10-31"), ]
  df
})


####FE 1 WAY####

model_list <- lapply(imputed_panels, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config,
    data = df
  )
})

library(miceadds)
qhat <- lapply(model_list, coef)
uhat <- lapply(model_list, vcov)
pooled <- miceadds::pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)
print(pooled_summary)

pooled_summary$stars <- cut(
  pooled_summary$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "FE_1WAY_feols.csv")


r2_values <- sapply(model_list, function(m) summary(m)$sq.cor)


within_r2_values <- sapply(model_list, function(m) summary(m)$within)


n_values <- sapply(model_list, function(m) summary(m)$nobs)

cat("Average R²: ", mean(r2_values), "\n")
cat("Average N: ", mean(n_values), "\n")

summary(model_list[[1]])

model_list <- lapply(imputed_panels, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config,
    cluster = ~distinct_car_config,
    data = df
  )
})

qhat <- lapply(model_list, coef)
uhat <- lapply(model_list, vcov)
pooled <- miceadds::pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)
print(pooled_summary)

pooled_summary$stars <- cut(
  pooled_summary$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "FE_1WAY_feols.csv")


r2_values <- sapply(model_list, function(m) summary(m)$sq.cor)


within_r2_values <- sapply(model_list, function(m) summary(m)$within)


n_values <- sapply(model_list, function(m) summary(m)$nobs)

cat("Average R²: ", mean(r2_values), "\n")
cat("Average N: ", mean(n_values), "\n")

summary(model_list[[1]])
####FE 2 WAY####

model_list <- lapply(imputed_panels, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk,
    data = df
  )
})

library(miceadds)
qhat <- lapply(model_list, coef)
uhat <- lapply(model_list, vcov)
pooled <- miceadds::pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)
print(pooled_summary)

pooled_summary$stars <- cut(
  pooled_summary$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "FE_2WAY_feols.csv")


r2_values <- sapply(model_list, function(m) summary(m)$sq.cor)


within_r2_values <- sapply(model_list, function(m) summary(m)$within)


n_values <- sapply(model_list, function(m) summary(m)$nobs)

cat("Average R²: ", mean(r2_values), "\n")
cat("Average N: ", mean(n_values), "\n")

summary(model_list[[1]])


####POLS####


model_list_pols <- lapply(imputed_panels, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
    data = df
  )
})

library(miceadds)


qhat_pols <- lapply(model_list_pols, coef)
uhat_pols <- lapply(model_list_pols, vcov)
pooled_pols <- miceadds::pool_mi(qhat = qhat_pols, u = uhat_pols)
pooled_summary_pols <- summary(pooled_pols)
print(pooled_summary_pols)


pooled_summary_pols$stars <- cut(
  pooled_summary_pols$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary_pols)
write.csv(pooled_summary_pols, "POLS_pooled_summary.csv")


r2_pols <- sapply(model_list_pols, function(m) summary(m)$sq.cor)
n_pols  <- sapply(model_list_pols, function(m) summary(m)$nobs)

cat("Average R²: ", mean(r2_pols), "\n")
cat("Average N: ", mean(n_pols), "\n")

####RE####

vars_used <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "seats", "lnweight", "distinct_car_config")
ln_vars   <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "lnweight") # ln vars only!

model_list_re <- lapply(imputed_panels, function(df) {
  df <- df[
    complete.cases(df[, vars_used]) &
      apply(df[, ln_vars], 1, function(x) all(is.finite(x))),
  ]
  plm(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
    data  = df,
    index = "distinct_car_config",
    model = "random"
  )
})

library(miceadds)
qhat_re <- lapply(model_list_re, coef)
uhat_re <- lapply(model_list_re, vcov)
pooled_re <- miceadds::pool_mi(qhat = qhat_re, u = uhat_re)
pooled_summary_re <- summary(pooled_re)
print(pooled_summary_re)

pooled_summary_re$stars <- cut(
  pooled_summary_re$p,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary_re)

r2_re <- sapply(model_list_re, function(m) summary(m)$r.squared[1])
n_re  <- sapply(model_list_re, function(m) nobs(m))

cat("Average R² (RE): ", mean(r2_re, na.rm=TRUE), "\n")
cat("Average N (RE): ", mean(n_re, na.rm=TRUE), "\n")




####POLS plm####

library(plm)


df <- imputed_panels[[1]]

vars_used <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "seats", "lnweight", "distinct_car_config")
ln_vars   <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "lnweight")


df <- df[
  complete.cases(df[, vars_used]) & 
    apply(df[, ln_vars], 1, function(x) all(is.finite(x))),
]


model_pols <- plm(
  lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
  data  = df,
  index = "distinct_car_config",
  model = "pooling"
)


model_fe <- plm(
  lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
  data  = df,
  index = "distinct_car_config",
  model = "within"
)


model_re <- plm(
  lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
  data  = df,
  index = "distinct_car_config",
  model = "random"
)

# F-test for fixed effects
f_test <- pFtest(model_fe, model_pols)
print(f_test)

# Hausman test: FE vs RE
hausman_test <- phtest(model_fe, model_re)
print(hausman_test)


cat("POLS R²:   ", summary(model_pols)$r.squared[1], "\n")
cat("FE R²:     ", summary(model_fe)$r.squared[1], "\n")
cat("RE R²:     ", summary(model_re)$r.squared[1], "\n")
cat("Obs (FE):  ", nobs(model_fe), "\n")
cat("Obs (POLS):", nobs(model_pols), "\n")
cat("Obs (RE):  ", nobs(model_re), "\n")

summary(model_fe)
summary(model_re)
summary(model_pols)

coefs_pols <- summary(model_pols)$coefficients
coefs_pols_df <- as.data.frame(coefs_pols)


coefs_pols_df$stars <- cut(
  coefs_pols_df[,"Pr(>|t|)"],
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)


write.csv(coefs_pols_df, "POLS_summary_stars.csv")

coefs_re <- summary(model_re)$coefficients
coefs_re_df <- as.data.frame(coefs_re)

coefs_re_df$stars <- cut(
  coefs_re_df[,"Pr(>|z|)"],
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)

write.csv(coefs_re_df, "RE_summary.csv")



vars_used <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "seats", "lnweight", "distinct_car_config")
ln_vars   <- c("lnreg", "lnrfee", "lnfuelprice", "lnco2", "lnweight")


filter_good <- function(df) {
  df[complete.cases(df[, vars_used]) &
       apply(df[, ln_vars], 1, function(x) all(is.finite(x))),
  ]
}


model_list_pols <- lapply(imputed_panels, function(df) {
  df <- filter_good(df)
  plm(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
    data  = df,
    index = "distinct_car_config",
    model = "pooling"
  )
})
model_list_fe <- lapply(imputed_panels, function(df) {
  df <- filter_good(df)
  plm(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
    data  = df,
    index = "distinct_car_config",
    model = "within"
  )
})
model_list_re <- lapply(imputed_panels, function(df) {
  df <- filter_good(df)
  plm(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight,
    data  = df,
    index = "distinct_car_config",
    model = "random"
  )
})

pooled_pols <- miceadds::pool_mi(
  qhat = lapply(model_list_pols, coef),
  u    = lapply(model_list_pols, vcov)
)
pooled_fe <- miceadds::pool_mi(
  qhat = lapply(model_list_fe, coef),
  u    = lapply(model_list_fe, vcov)
)
pooled_re <- miceadds::pool_mi(
  qhat = lapply(model_list_re, coef),
  u    = lapply(model_list_re, vcov)
)

pooled_summary_pols <- summary(pooled_pols)
pooled_summary_fe   <- summary(pooled_fe)
pooled_summary_re   <- summary(pooled_re)

add_stars <- function(df, p_col="p") {
  df$stars <- cut(
    df[[p_col]],
    breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
    labels = c("***", "**", "*", ".", "")
  )
  df
}
pooled_summary_pols <- add_stars(pooled_summary_pols)
pooled_summary_fe   <- add_stars(pooled_summary_fe)
pooled_summary_re   <- add_stars(pooled_summary_re)

write.csv(pooled_summary_pols, "POLS_pooled_summary.csv")
write.csv(pooled_summary_fe,   "FE_pooled_summary.csv")
write.csv(pooled_summary_re,   "RE_pooled_summary.csv")

r2_pols <- sapply(model_list_pols, function(m) summary(m)$r.squared[1])
r2_fe   <- sapply(model_list_fe,   function(m) summary(m)$r.squared[1])
r2_re   <- sapply(model_list_re,   function(m) summary(m)$r.squared[1])
n_pols  <- sapply(model_list_pols, nobs)
n_fe    <- sapply(model_list_fe,   nobs)
n_re    <- sapply(model_list_re,   nobs)

cat("Average R² (POLS):", mean(r2_pols, na.rm=TRUE), "\n")
cat("Average R² (FE):  ", mean(r2_fe, na.rm=TRUE), "\n")
cat("Average R² (RE):  ", mean(r2_re, na.rm=TRUE), "\n")
cat("Average N (POLS):", mean(n_pols, na.rm=TRUE), "\n")
cat("Average N (FE):  ", mean(n_fe, na.rm=TRUE), "\n")
cat("Average N (RE):  ", mean(n_re, na.rm=TRUE), "\n")

adj_r2_re <- sapply(model_list_re, function(m) summary(m)$r.squared["adjrsq"])
cat("Average Adjusted R² (RE):", mean(adj_r2_re, na.rm=TRUE), "\n")

adj_r2_pols <- sapply(model_list_pols, function(m) summary(m)$r.squared["adjrsq"])
cat("Average Adjusted R² (POLS):", mean(adj_r2_pols, na.rm=TRUE), "\n")

adj_r2_fe <- sapply(model_list_fe, function(m) summary(m)$r.squared["adjrsq"])
cat("Average Adjusted R² (FE):", mean(adj_r2_fe, na.rm=TRUE), "\n")




imputed_panels <- lapply(imputed_panels, function(df) {
  df$segment <- with(df, ifelse(
    condition == "New" & power_kw <= 80,      "new_low",
    ifelse(condition == "New" & power_kw <= 130, "new_mid",
           ifelse(condition == "New" & power_kw > 130,  "new_high",
                  ifelse(condition == "Used" & power_kw <= 80,      "used_low",
                         ifelse(condition == "Used" & power_kw <= 130,     "used_mid",
                                ifelse(condition == "Used" & power_kw > 130,      "used_high", NA_character_)
                         ))))))
  df
})


library(fixest)

library(fixest)


fe_formula <- lnreg ~ lnrfee + lnrfee:relevel(factor(segment), ref = "new_mid") + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk


model_list_2wayFE <- lapply(imputed_panels, function(df) {
  feols(
    fml = fe_formula,
    data = df,
    cluster = ~distinct_car_config
  )
})


library(miceadds)
qhat <- lapply(model_list_2wayFE, coef)
uhat <- lapply(model_list_2wayFE, vcov)
pooled <- pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)

pooled_summary$stars <- cut(
  pooled_summary$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "2wayFE_segment_pooled.csv")

summary(model_list_2wayFE[[5]])


library(dplyr)
library(ggplot2)

df <- imputed_panels[[1]]


df$evid_svk <- as.Date(df$evid_svk)


plot_data <- df %>%
  group_by(evid_svk, segment) %>%
  summarise(mean_lnreg = mean(lnreg, na.rm = TRUE), .groups = "drop")


print(head(plot_data))


ggplot(plot_data, aes(x = evid_svk, y = mean_lnreg, color = segment, group = segment)) +
  geom_line(size = 0.7) +
  labs(
    title = "Mean ln(registrations) by Segment Over Time",
    x = "Time (evid_svk)",
    y = "Mean ln(reg_no)",
    color = "Segment"
  ) +
  theme_minimal() +
  theme(legend.position = "right")

####ROBUSTNESS CHECKS####

library(dplyr)

library(dplyr)

panels_bm <- lapply(imputed_panels, function(df) {
  # Create brand_model column
  df <- df %>%
    mutate(brand_model = paste(brand, model, sep = "_"))
  
  # Aggregate
  df %>%
    group_by(brand_model, evid_svk) %>%
    summarize(
      lnreg = sum(lnreg, na.rm = TRUE),         
      lnrfee = mean(lnrfee, na.rm = TRUE),
      lnfuelprice = mean(lnfuelprice, na.rm = TRUE),
      lnkw = mean(lnkw, na.rm = TRUE),
      lnco2 = mean(lnco2, na.rm = TRUE),
      seats = mean(seats),                      
      lnweight = mean(lnweight, na.rm = TRUE),
      .groups = "drop"
    )
})

library(fixest)
library(miceadds)

model_list <- lapply(panels_bm, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnkw + lnco2 + seats + lnweight | brand_model + evid_svk,
    data = df
  )
})

qhat <- lapply(model_list, coef)
uhat <- lapply(model_list, vcov)
pooled <- miceadds::pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)
print(pooled_summary)

pooled_summary$stars <- cut(
  pooled_summary$p, 
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "FE_2WAY_feols_brandmodel.csv")

r2_values <- sapply(model_list, function(m) summary(m)$sq.cor)
within_r2_values <- sapply(model_list, function(m) summary(m)$within)
n_values <- sapply(model_list, function(m) summary(m)$nobs)

cat("Average R²: ", mean(r2_values, na.rm = TRUE), "\n")
cat("Average N: ", mean(n_values, na.rm = TRUE), "\n")

summary(model_list[[2]])


imputed_panels_new <- lapply(imputed_panels, function(df) {
  df %>% filter(condition == "New")
})

imputed_panels_used <- lapply(imputed_panels, function(df) {
  df %>% filter(condition == "Used")
})

# For NEW cars
model_list_new <- lapply(imputed_panels_new, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk,
    data = df
  )
})
qhat_new <- lapply(model_list_new, coef)
uhat_new <- lapply(model_list_new, vcov)
pooled_new <- miceadds::pool_mi(qhat = qhat_new, u = uhat_new)
pooled_summary_new <- summary(pooled_new)
pooled_summary_new$stars <- cut(
  pooled_summary_new$p,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
write.csv(pooled_summary_new, "FE_2WAY_newcars.csv")
print(pooled_summary_new)

# For USED cars
model_list_used <- lapply(imputed_panels_used, function(df) {
  feols(
    lnreg ~ lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk,
    data = df
  )
})
qhat_used <- lapply(model_list_used, coef)
uhat_used <- lapply(model_list_used, vcov)
pooled_used <- miceadds::pool_mi(qhat = qhat_used, u = uhat_used)
pooled_summary_used <- summary(pooled_used)
pooled_summary_used$stars <- cut(
  pooled_summary_used$p,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
write.csv(pooled_summary_used, "FE_2WAY_usedcars.csv")
print(pooled_summary_used)

summary(model_list_used[[5]])


reform_date <- as.Date("2012-10-01")

imputed_panels_2012 <- lapply(imputed_panels, function(df) {
  df <- as.data.frame(df)
  evid_svk_chr <- as.character(df$evid_svk)
  evid_dates <- suppressWarnings(as.Date(evid_svk_chr))
  if (all(is.na(evid_dates))) {
    evid_dates <- as.Date(paste0(evid_svk_chr, "-01"))
  }
  df$month_date <- as.Date(evid_dates)
  df$power_kw <- as.numeric(df$power_kw)
  df <- df[!is.na(df$month_date), ]
  
  df$months_from_reform <- as.integer((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12)
  
  # Restrict to 24 months before and after
  df <- df[df$months_from_reform >= -24 & df$months_from_reform <= 24, ]
  
  df$treated <- as.integer(df$power_kw >= 81)
  df$post_reform <- as.integer(df$month_date >= reform_date)
  return(df)
})


library(ggplot2)
plot_df <- imputed_panels_2012[[1]]

plot_df_summarized <- plot_df %>%
  group_by(month_date, treated) %>%
  summarize(total_regs = sum(regs_no, na.rm = TRUE), .groups = "drop") %>%
  mutate(
    log_total_regs = log1p(total_regs),
    group_label = factor(
      treated,
      levels = c(0, 1),
      labels = c("Control (≤80 kW)", "Treated (>80 kW)")
    )
  )

ggplot(plot_df_summarized, aes(x = month_date, y = log_total_regs, color = group_label)) +
  geom_line(size = 0.8) +
  labs(
    x = "Month",
    y = "log(Total registrations + 1)",
    color = ""
  ) +
  geom_vline(xintercept = as.numeric(reform_date), linetype = "dashed", color = "black") +
  ggtitle("Log Total Registrations by Group (Pre/Post 2012 Fee Reform)") +
  theme_minimal(base_size = 16)


did_model_list <- lapply(imputed_panels_2012, function(df) {
  feols(
    lnreg ~ treated + post_reform + treated * post_reform + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk,
    data = df
  )
})


qhat <- lapply(did_model_list, coef)
uhat <- lapply(did_model_list, vcov)
pooled <- miceadds::pool_mi(qhat = qhat, u = uhat)
pooled_summary <- summary(pooled)


pooled_summary$stars <- cut(
  pooled_summary$p,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
  labels = c("***", "**", "*", ".", "")
)
print(pooled_summary)
write.csv(pooled_summary, "DID_FEOLS_summary.csv")

summary(did_model_list[[5]])

library(dplyr)
library(lubridate)

reform_date <- as.Date("2012-10-01")


imputed_panels_2012 <- lapply(imputed_panels_2012, function(df) {
  df$month_date <- as.Date(df$month_date)
  df$event_time <- interval(reform_date, df$month_date) %/% months(1)
  df$event_time <- as.integer(round((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12))
  df
})

library(fixest)

event_study_models <- lapply(imputed_panels_2012, function(df) {
  df <- df[!is.na(df$event_time), ]
  
  feols(
    lnreg ~ i(event_time, treated, ref = -4) + lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk,
    data = df
  )
})

library(ggplot2)
es_mod <- event_study_models[[1]]

summary(event_study_models[[1]])

coefs <- broom::tidy(es_mod, conf.int = TRUE)
coefs_es <- coefs[grepl("^event_time::", coefs$term), ]
coefs_es$event_time <- as.integer(gsub("event_time::(-?\\d+).*", "\\1", coefs_es$term))

ggplot(coefs_es, aes(x = event_time, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = 0, color = "grey") +
  labs(
    x = "Months from reform",
    y = "Event study coefficient (event_time x treated)",
    title = "Event Study: Effect of 2012 Fee Reform by Months Relative 4 Months Before the Reform"
  ) +
  theme_minimal(base_size = 16)

library(lubridate)
library(dplyr)

library(lubridate)
library(dplyr)

reform_date <- as.Date("2017-02-01")
window_months <- 24

make_event_groups <- function(df) {
  if (!inherits(df$evid_svk, "Date")) {
    df$month_date <- as.Date(as.character(df$evid_svk))
  } else {
    df$month_date <- df$evid_svk
  }
  
  df$event_time <- interval(reform_date, df$month_date) %/% months(1)
  
  df <- df %>%
    filter(
      month_date >= (reform_date %m-% months(window_months)) &
        month_date <= (reform_date %m+% months(window_months))
    )
  
  # Assign treatment group
  df$treatment_group <- NA_character_
  df$treatment_group[df$power_kw <= 80] <- "control"
  df$treatment_group[
    df$condition == "New" & df$power_kw > 80 & df$power_kw <= 176
  ] <- "new_dec"
  df$treatment_group[
    df$condition == "New" & df$power_kw > 176
  ] <- "new_inc"
  df$treatment_group[
    df$condition == "Used" & df$power_kw > 80 & df$power_kw <= 176
  ] <- "used_dec"
  df$treatment_group[
    df$condition == "Used" & df$power_kw > 176
  ] <- "used_inc"
  
  n_na <- sum(is.na(df$treatment_group))
  cat("Number of NAs in treatment_group:", n_na, "\n")
  if (n_na > 0) {
    print("Rows with NA treatment_group:")
    print(head(df[is.na(df$treatment_group), ]))
  }
  
  return(df)
}

imputed_panels_2017 <- lapply(imputed_panels, make_event_groups)



plot_df <- imputed_panels_2017[[1]]


plot_df_summarized <- plot_df %>%
  group_by(month_date, treatment_group) %>%
  summarize(total_regs = sum(regs_no, na.rm = TRUE), .groups = "drop")


plot_df_summarized <- plot_df_summarized %>%
  mutate(log_total_regs = log1p(total_regs))


ggplot(plot_df_summarized, aes(x = month_date, y = log_total_regs, color = treatment_group)) +
  geom_line(size = 0.9) +
  labs(
    title = "Log Total Registrations by Treatment Group (per month)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  theme_minimal() +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black")

df_sub1 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "used_inc")) %>%
  mutate(
    plot_group = case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "used_inc" ~ "Used, Fees Increased (177+ kW)"
    )
  )

ggplot(df_sub1, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.8) +
  labs(
    title = "Log Total Registrations: Control vs. Used, kW Tariff Increased (177+ kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  theme_minimal(base_size = 16) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  scale_color_manual(
    values = c("Control (<80kW)" = "#8f7868", "Used, Fees Increased (177+ kW)" = "red")
  )


df_sub2 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "used_dec")) %>%
  mutate(
    plot_group = case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "used_dec" ~ "Used, Fees Decreased (81-176 kW)"
    )
  )

ggplot(df_sub2, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.8) +
  labs(
    title = "Log Total Registrations: Control vs. Used, kW Tariff Decreased (81–176 kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  theme_minimal(base_size=16) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  scale_color_manual(
    values = c(
      "Control (<80kW)" = "#8f7868",
      "Used, Fees Decreased (81-176 kW)" = "blue"
    )
  )



df_sub3 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "new_inc")) %>%
  mutate(
    plot_group = case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "new_inc" ~ "New, Fees Increased (177+ kW)"
    )
  )

ggplot(df_sub3, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 1) +
  labs(
    title = "Log Total Registrations: Control vs. New, kW Tariff Increased (177+ kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  theme_minimal(base_size = 16) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black")+
  scale_color_manual(
    values = c(
      "Control (<80kW)" = "#8f7868",
      "New, Fees Increased (177+ kW)" = "#c730ae"
    )
  )


df_sub4 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "new_dec")) %>%
  mutate(
    plot_group = case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "new_dec" ~ "New, Fees Decreased (81–176 kW)"
    )
  )

ggplot(df_sub4, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 1) +
  labs(
    title = "Log Total Registrations: Control vs. New, kW Tariff Decreased (81–176 kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  theme_minimal(base_size = 16) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black")+
  scale_color_manual(
    values = c(
      "Control (<80kW)" = "#8f7868",
      "New, Fees Decreased (81–176 kW)" = "#30c788"
    )
  )



reform_date <- as.Date("2017-02-01")


imputed_panels_event <- lapply(imputed_panels_2017, function(df) {

  df$month_date <- as.Date(df$month_date)
  df$event_time <- as.integer(round((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12))
  df
})

imputed_panels_event <- lapply(imputed_panels_event, function(df) {
  df <- df[df$event_time >= -24 & df$event_time <= 24, ]
  df
})

df <- imputed_panels_event[[1]]


df$treatment_group <- relevel(factor(df$treatment_group), ref = "control")
df$event_time <- factor(df$event_time)


df$new_inc  <- as.integer(df$treatment_group == "new_inc")
df$used_inc <- as.integer(df$treatment_group == "used_inc")
df$new_dec  <- as.integer(df$treatment_group == "new_dec")
df$used_dec <- as.integer(df$treatment_group == "used_dec")



event_formula <- lnreg ~ 
  i(event_time, new_inc, ref = -4) +
  i(event_time, used_inc, ref = -4) +
  i(event_time, new_dec, ref = -4) +
  i(event_time, used_dec, ref = -4) +
  lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk

model_event_study <- feols(event_formula, data = df, cluster = ~distinct_car_config)

summary(model_event_study)

library(broom)
tidy_out <- broom::tidy(model_event_study, conf.int = TRUE)



event_tidy <- tidy_out %>%
  tidyr::extract(
    term,
    into = c("event_time", "group"),
    regex = "event_time::(-?\\d+):([a-z_]+)",
    remove = FALSE
  ) %>%
  mutate(
    event_time = as.integer(event_time),
    group = factor(group)
  )

plot_group <- "new_inc"

event_plot_df <- event_tidy %>% filter(group == plot_group)

library(ggplot2)

event_plot_df <- event_tidy %>% filter(group == "new_inc")

event_plot_df <- event_tidy %>% filter(group == "new_inc")
ggplot(event_plot_df, aes(x = as.integer(event_time), y = estimate)) +
  geom_point(size = 2, color = "blue") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = "blue") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +   # Reform
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +  # Zero effect
  labs(
    title = "Event Study: 2017 Reform New Vehicles, kW Tariff Fees Increased (177+ kW)",
    x = "Months from reform",
    y = "Event study coefficient (event_time x group)"
  ) +
  theme_minimal(base_size = 16)

event_plot_df <- event_tidy %>% filter(group == "used_inc")
ggplot(event_plot_df, aes(x = as.integer(event_time), y = estimate)) +
  geom_point(size = 2, color = "red") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = "red") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +   # Reform line
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +  # Zero line
  labs(
    title = "Event Study: Used, Fees Increased (177+ kW)",
    x = "Months from Reform (Feb 2017)",
    y = "Effect on ln(registrations)"
  ) +
  theme_minimal()

event_plot_df <- event_tidy %>% filter(group == "new_dec")
ggplot(event_plot_df, aes(x = as.integer(event_time), y = estimate)) +
  geom_point(size = 2, color = "forestgreen") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = "forestgreen") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  labs(
    title = "Event Study: New, Fees Decreased (81–176 kW)",
    x = "Months from Reform (Feb 2017)",
    y = "Effect on ln(registrations)"
  ) +
  theme_minimal()

event_plot_df <- event_tidy %>% filter(group == "used_dec")
ggplot(event_plot_df, aes(x = as.integer(event_time), y = estimate)) +
  geom_point(size = 2, color = "purple") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = "purple") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  labs(
    title = "Event Study: Used, Fees Decreased (81–176 kW)",
    x = "Months from Reform (Feb 2017)",
    y = "Effect on ln(registrations)"
  ) +
  theme_minimal()


reform_date <- as.Date("2017-02-01")
window_months <- 24

make_event_groups <- function(df) {
  if (!"month_date" %in% names(df)) {
    if ("evid_svk" %in% names(df)) {
      df$month_date <- as.Date(as.character(df$evid_svk))
    } else {
      stop("No date variable found!")
    }
  }
  
  df$event_time <- as.integer(round((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12))
  df <- df %>%
    filter(
      month_date >= (reform_date %m-% months(window_months)) &
        month_date <= (reform_date %m+% months(window_months))
    )
  # Assign new group names (focus on kW Bracket, not effect direction)
  df$treatment_group <- NA_character_
  df$treatment_group[df$power_kw <= 80] <- "control"
  df$treatment_group[df$condition == "New"  & df$power_kw > 80  & df$power_kw <= 176] <- "new_kwdec"
  df$treatment_group[df$condition == "New"  & df$power_kw > 176]                   <- "new_kwinc"
  df$treatment_group[df$condition == "Used" & df$power_kw > 80  & df$power_kw <= 176] <- "used_kwdec"
  df$treatment_group[df$condition == "Used" & df$power_kw > 176]                   <- "used_kwinc"
  df$treatment_group <- factor(df$treatment_group, levels = c("control", "new_kwdec", "new_kwinc", "used_kwdec", "used_kwinc"))
  return(df)
}
imputed_panels_2017 <- lapply(imputed_panels, make_event_groups)
library(dplyr)
library(ggplot2)


table(plot_df_summarized$treatment_group)



df_sub1 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "used_kwinc")) %>%
  mutate(
    plot_group = dplyr::case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "used_kwinc" ~ "Used, kW Bracket: Fee Increased (177+ kW)"
    )
  )
ggplot(df_sub1, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.7) +
  labs(
    title = "Log Total Registrations: Control vs. Used, kW Bracket Fee Increased (177+ kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control (<80kW)" = "#8f7868", "Used, kW Bracket: Fee Increased (177+ kW)" = "red")) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  theme_minimal()

# ---- Control vs Used, kW Bracket: Fee Decreased (81–176 kW) ----
df_sub2 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "used_kwdec")) %>%
  mutate(
    plot_group = dplyr::case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "used_kwdec" ~ "Used, kW Bracket: Fee Decreased (81–176 kW)"
    )
  )
ggplot(df_sub2, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.7) +
  labs(
    title = "Log Total Registrations: Control vs. Used, kW Bracket Fee Decreased (81–176 kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control (<80kW)" = "#8f7868", "Used, kW Bracket: Fee Decreased (81–176 kW)" = "blue")) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  theme_minimal()

# ---- Control vs New, kW Bracket: Fee Increased (177+ kW) ----
df_sub3 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "new_kwinc")) %>%
  mutate(
    plot_group = dplyr::case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "new_kwinc" ~ "New, kW Bracket: Fee Increased (177+ kW)"
    )
  )
ggplot(df_sub3, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.7) +
  labs(
    title = "Log Total Registrations: Control vs. New, kW Bracket Fee Increased (177+ kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control (<80kW)" = "#8f7868", "New, kW Bracket: Fee Increased (177+ kW)" = "#c730ae")) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  theme_minimal()

# ---- Control vs New, kW Bracket: Fee Decreased (81–176 kW) ----
df_sub4 <- plot_df_summarized %>%
  filter(treatment_group %in% c("control", "new_kwdec")) %>%
  mutate(
    plot_group = dplyr::case_when(
      treatment_group == "control" ~ "Control (<80kW)",
      treatment_group == "new_kwdec" ~ "New, kW Bracket: Fee Decreased (81–176 kW)"
    )
  )
ggplot(df_sub4, aes(x = month_date, y = log_total_regs, color = plot_group)) +
  geom_line(size = 0.7) +
  labs(
    title = "Log Total Registrations: Control vs. New, kW Bracket Fee Decreased (81–176 kW)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control (<80kW)" = "#8f7868", "New, kW Bracket: Fee Decreased (81–176 kW)" = "#30c788")) +
  geom_vline(xintercept = as.Date("2017-02-01"), linetype = "dashed", color = "black") +
  theme_minimal()

# ---- EVENT STUDY: Prepare variables ----
df <- imputed_panels_2017[[1]]
df$treatment_group <- relevel(factor(df$treatment_group), ref = "control")
df$event_time <- as.integer(round((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12))
df$event_time <- factor(df$event_time)

# Make dummies for each group
df$new_kwinc  <- as.integer(df$treatment_group == "new_kwinc")
df$used_kwinc <- as.integer(df$treatment_group == "used_kwinc")
df$new_kwdec  <- as.integer(df$treatment_group == "new_kwdec")
df$used_kwdec <- as.integer(df$treatment_group == "used_kwdec")

# ---- EVENT STUDY MODEL ----
event_formula <- lnreg ~ 
  i(event_time, new_kwinc, ref = -4) +
  i(event_time, used_kwinc, ref = -4) +
  i(event_time, new_kwdec, ref = -4) +
  i(event_time, used_kwdec, ref = -4) +
  lnrfee + lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk

model_event_study <- feols(event_formula, data = df, cluster = ~distinct_car_config)
tidy_out <- broom::tidy(model_event_study, conf.int = TRUE)

# Parse the term column to get event_time and group for plotting
event_tidy <- tidy_out %>%
  tidyr::extract(
    term,
    into = c("event_time", "group"),
    regex = "event_time::(-?\\d+):([a-z_]+)",
    remove = FALSE
  ) %>%
  mutate(
    event_time = as.integer(event_time),
    group = factor(group, levels = c("new_kwinc", "used_kwinc", "new_kwdec", "used_kwdec"))
  ) %>%
  filter(!is.na(group))

# ---- EVENT STUDY PLOTS ----

plot_event_study_group <- function(group_code, group_label, color){
  event_plot_df <- event_tidy %>% filter(group == group_code)
  ggplot(event_plot_df, aes(x = event_time, y = estimate)) +
    geom_point(size = 2, color = color) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = color) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "black") +   # Reform line
    geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +  # Zero effect
    labs(
      title = paste0("Event Study 2017: ", group_label),
      x = "Months from reform",
      y = "Event study coefficient (event_time x group)"
    ) +
    theme_minimal(base_size = 16)
}

# Used, kW Bracket Fee Increased (177+ kW)
plot_event_study_group("used_kwinc", "Used, kW Bracket: Fee Increased (177+ kW)", "red")
# Used, kW Bracket Fee Decreased (81–176 kW)
plot_event_study_group("used_kwdec", "Used, kW Bracket: Fee Decreased (81–176 kW)", "blue")
# New, kW Bracket Fee Increased (177+ kW)
plot_event_study_group("new_kwinc", "New, kW Bracket: Fee Increased (177+ kW)", "#c730ae")
# New, kW Bracket Fee Decreased (81–176 kW)
plot_event_study_group("new_kwdec", "New, kW Bracket: Fee Decreased (81–176 kW)", "#30c788")



reform_date <- as.Date("2023-07-01")
window_months <- 12

df <- imputed_panels[[1]]

if (!"month_date" %in% names(df)) {
  df$month_date <- as.Date(as.character(df$evid_svk))
}
df$month_date <- as.Date(df$month_date)

# Restrict to ±12 months around July 2023
reform_date <- as.Date("2023-07-01")
window_months <- 12
df <- df %>%
  filter(
    month_date >= (reform_date %m-% months(window_months)) &
      month_date <= (reform_date %m+% months(window_months))
  )

# Assign treated dummy: 1 if power_kw > 80, 0 otherwise
df$treated <- ifelse(df$power_kw > 80, 1, 0)
df_summarized <- df %>%
  group_by(month_date, treated) %>%
  summarize(total_regs = sum(regs_no, na.rm = TRUE), .groups = "drop") %>%
  mutate(
    log_total_regs = log(total_regs + 1),  # log(0) protection
    group_label = factor(
      treated,
      levels = c(0, 1),
      labels = c("Control (≤80kW)", "Treated (>80kW)")
    )
  )

ggplot(df_summarized, aes(x = month_date, y = log_total_regs, color = group_label)) +
  geom_line(size = 1.1) +
  geom_vline(xintercept = reform_date, linetype = "dashed", color = "black") +
  labs(
    title = "Log of Total Registrations by Group (±12 months around July 2023 Reform)",
    x = "Month",
    y = "log(Total Registrations + 1)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control (≤80kW)" = "#8f7868", "Treated (>80kW)" = "red")) +
  theme_minimal(base_size = 16)

library(lubridate)
library(dplyr)
library(fixest)


df <- imputed_panels[[1]]


if (!"month_date" %in% names(df)) {
  df$month_date <- as.Date(as.character(df$evid_svk))
}
df$month_date <- as.Date(df$month_date)


reform_date <- as.Date("2023-07-01")
window_months <- 12


df <- df %>%
  filter(
    month_date >= (reform_date %m-% months(window_months)) &
      month_date <= (reform_date %m+% months(window_months))
  )

# Treated dummy
df$treated <- ifelse(df$power_kw > 80, 1, 0)

# Event time in months (distance from reform, 0 = July 2023)
df$event_time <- as.integer(round((as.yearmon(df$month_date) - as.yearmon(reform_date)) * 12))



# Reference period
ref_period <- -6

event_formula <- lnreg ~ i(event_time, treated, ref = ref_period) + lnrfee +
  lnfuelprice + lnco2 + seats + lnweight | distinct_car_config + evid_svk

model_event_study <- feols(event_formula, data = df, cluster = ~distinct_car_config)

summary(model_event_study)

tidy_out <- broom::tidy(model_event_study, conf.int = TRUE)

event_tidy <- tidy_out %>%
  tidyr::extract(
    term,
    into = "event_time",
    regex = "event_time::(-?\\d+):treated",
    remove = FALSE
  ) %>%
  mutate(event_time = as.integer(event_time))

ggplot(event_tidy, aes(x = event_time, y = estimate)) +
  geom_point(size = 2, color = "darkblue") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.5, color = "darkblue") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  labs(
    title = "Event Study: 2023 Treated (81+ kW) vs Control (<80 kW)",
    x = "Months from reform",
    y = "Event study coefficients (event_time x treated)"
  ) +
  theme_minimal(base_size = 16)