################################################################################
# Replication Script for:
# "Unveiling the Environmental Cost of Bitcoin: A Quantile-Based Investigation
#  of Energy Use and Carbon Emissions"
# Journal: Humanities and Social Sciences Communications (Springer Nature)
################################################################################

# --- 0. Package Installation and Loading ---
required_packages <- c(
  "readxl", "zoo", "xts", "quantmod", "vars", "urca", "fUnitRoots", 
  "aod", "tseries", "dplyr", "tidyr", "ggplot2", "gridExtra", 
  "reshape2", "quantreg", "KernSmooth", "np", "ConnectednessApproach", 
  "PerformanceAnalytics", "plotly", "htmlwidgets"
)

new_packages <- required_packages[!(required_packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages, dependencies = TRUE)

invisible(lapply(required_packages, library, character.only = TRUE))

# --- 1. Data Ingestion & Preprocessing ---

# 1.1 Bitcoin Electricity Consumption (CBECI)
load_excel_to_zoo <- function(file_path, sheet = 1, date_col = "Date and Time", 
                              value_col = "annualised consumption GUESS, TWh", 
                              date_format = "%Y-%m-%dT%H:%M:%S") {
  data <- read_excel(file_path, sheet = sheet)
  data[[date_col]] <- as.POSIXct(data[[date_col]], format = date_format)
  data <- data[!is.na(data[[date_col]]), ]
  zoo_obj <- zoo(data[[value_col]], order.by = as.Date(data[[date_col]]))
  return(zoo_obj)
}

electricity_file <- "Cleaned_Historical_Annualised_Electricity_Consumption_With_Headers.xlsx"
electricity_zoo  <- load_excel_to_zoo(electricity_file)

# 1.2 Bitcoin Price (Daily Close)
btc_raw   <- read.csv("bitcoin_price_2014_09_17_to_2024_12_17.csv", header = TRUE)
btc_zoo   <- zoo(btc_raw[, "Close"], order.by = as.Date(strptime(as.character(btc_raw[, "Date"]), "%d/%m/%Y")))
BTCClose  <- btc_zoo

# 1.3 Bitcoin Carbon Emissions (CBECI GHG)
btcco_raw <- read.csv("Export_Historical_Bitcoin_Emissions.csv", header = TRUE)
BTCCO.zoo <- zoo(btcco_raw[, 2], order.by = as.Date(strptime(as.character(btcco_raw[, 1]), "%d/%m/%Y")))

# 1.4 Optional / Auxiliary Data
if (file.exists("data_gpr_daily.csv")) {
  gprd_raw <- read.csv("data_gpr_daily.csv", header = TRUE)
  GPRD.zoo <- zoo(gprd_raw[, -1], order.by = as.Date(strptime(as.character(gprd_raw[, 1]), "%d/%m/%Y")))
}
if (file.exists("fear_and_greed_index_all_data.csv")) {
  fgi_raw <- read.csv("fear_and_greed_index_all_data.csv", header = TRUE)
  FGI.zoo <- zoo(fgi_raw[, -1], order.by = as.Date(strptime(as.character(fgi_raw[, 1]), "%d/%m/%Y")))
}

# --- 2. Alignment and Return Transformation ---

# 2.1 BTC Price and Electricity Consumption
start_common <- as.Date("2014-09-17")
end_common   <- as.Date("2024-12-14")

BTCPEC    <- merge(window(BTCClose, start = start_common, end = end_common),
                   window(electricity_zoo, start = start_common, end = end_common))
BTCPECRet <- diff(log(BTCPEC))

# 2.2 BTC Price and Carbon Emissions
BTCCO2    <- merge(window(BTCClose, start = start_common, end = end_common),
                   window(BTCCO.zoo, start = start_common, end = end_common))
BTCCO2Ret <- diff(log(BTCCO2))

# 2.3 Master Dataset
dataall <- merge(BTCPECRet[, 1], BTCPECRet[, 2], BTCCO2Ret[, 2])
dataall <- na.omit(dataall)
colnames(dataall) <- c("BTC_Price", "BTC_El_Cons", "BTC_CO2_Cons")

# Summary Statistics Export
summary_stats <- data.frame(
  Mean     = apply(dataall, 2, mean),
  SD       = apply(dataall, 2, sd),
  Median   = apply(dataall, 2, median),
  Min      = apply(dataall, 2, min),
  Max      = apply(dataall, 2, max),
  Skewness = apply(dataall, 2, function(x) PerformanceAnalytics::skewness(x)),
  Kurtosis = apply(dataall, 2, function(x) PerformanceAnalytics::kurtosis(x))
)
write.csv(summary_stats, "Returns_Summary_Statistics.csv")

# Plot Returns Series
jpeg("Returns_TimeSeries.jpg", width = 12, height = 8, units = "in", res = 300)
par(mfrow = c(3, 1), mar = c(3, 4, 2, 2))
for (k in 1:3) {
  chart.TimeSeries(dataall[, k], colorset = "darkblue", lwd = 1.2, 
                   ylab = colnames(dataall)[k], main = colnames(dataall)[k])
}
dev.off()

# --- 3. Quantile-on-Quantile (Q-Q) Connectedness Approach ---

run_qq_connectedness <- function(data_sub, prefix_name) {
  NAMES <- colnames(data_sub)
  m <- 5
  quantiles <- seq(0.05, 0.95, length.out = m)
  nlag <- 1
  t_len <- nrow(data_sub)
  window.size <- 200 + nlag
  t0 <- t_len - window.size + 1 + nlag
  
  NET_arr <- TCI_arr <- array(NA, c(t0, m, m), dimnames = list(1:t0, quantiles, quantiles))
  
  for (i in 1:m) {
    for (j in 1:m) {
      dca <- suppressMessages(ConnectednessApproach(
        data_sub, nlag = nlag, nfore = 20, model = "QVAR", connectedness = "Time",
        window.size = window.size, corrected = TRUE,
        VAR_config = list(QVAR = list(tau = c(quantiles[i], quantiles[j]), method = "fn"))
      ))
      TCI_arr[, i, j] <- dca$TCI
      NET_arr[, i, j] <- dca$NET[, 1]
    }
  }
  
  TCI_agg <- round(apply(TCI_arr, 2:3, mean), 1)
  NET_agg <- round(apply(NET_arr, 2:3, mean), 1)
  
  melt_TCI <- reshape2::melt(TCI_agg)
  melt_NET <- reshape2::melt(NET_agg)
  
  # Heatmap of Average Total Connectedness Index
  p_tci <- ggplot(data = melt_TCI, aes(x = Var2, y = Var1, fill = value)) + 
    geom_tile(color = "black") + 
    geom_text(aes(label = value), color = "black", size = 4) + 
    labs(x = NAMES[2], y = NAMES[1], title = paste("Averaged Q-Q TCI:", prefix_name)) +
    scale_fill_gradientn(colours = c("gold4", "white", "steelblue4")) +
    theme_minimal()
  
  ggsave(paste0(prefix_name, "_QQ_TCI_Heatmap.png"), p_tci, width = 8, height = 6, dpi = 300)
  
  # Dynamic Direct vs Reverse TCI
  net_mat  <- FROM_mat <- TO_mat <- array(NA, c(t0, m))
  nTCI_vec <- pTCI_vec <- array(NA, c(t0, 1))
  
  for (i in 1:t0) {
    pTCI_vec[i, ] <- mean(diag(TCI_arr[i, , ]))
    nTCI_vec[i, ] <- mean(diag(TCI_arr[i, m:1, ]))
    FROM_mat[i, ] <- rowMeans(TCI_arr[i, , ])
    TO_mat[i, ]   <- colMeans(TCI_arr[i, , ])
    net_mat[i, ]  <- TO_mat[i, ] - FROM_mat[i, ]
  }
  
  dates_sub <- as.Date(tail(index(data_sub), t0))
  tci_df <- data.frame(
    Date = dates_sub,
    Reverse_TCI = as.numeric(nTCI_vec),
    Direct_TCI  = as.numeric(pTCI_vec),
    Delta_TCI   = as.numeric(nTCI_vec - pTCI_vec)
  )
  
  write.csv(zoo(net_mat[, c(1, 3, 5)], order.by = dates_sub), paste0("net_", prefix_name, ".csv"))
  return(list(tci_df = tci_df, dates = dates_sub))
}

# 3.1 BTC Price vs Electricity Consumption
res_el  <- run_qq_connectedness(dataall[, c("BTC_Price", "BTC_El_Cons")], "BTC_El_Cons")

# 3.2 BTC Price vs Carbon Emissions
res_co2 <- run_qq_connectedness(dataall[, c("BTC_Price", "BTC_CO2_Cons")], "BTC_CO2_Cons")

# --- 4. Quantile Regression Models (QQR) ---
# Note: For non-parametric QQR bivariate surfaces, ensure local custom functions 
# (such as QQR / plot.QQR from your working directory) or standard 'quantreg::rq' are sourced:
taus <- seq(0.05, 0.95, by = 0.10)
qr_el  <- rq(BTC_El_Cons ~ BTC_Price, tau = taus, data = as.data.frame(dataall))
qr_co2 <- rq(BTC_CO2_Cons ~ BTC_Price, tau = taus, data = as.data.frame(dataall))

save.image(file = "Replication_Workspace.RData")
cat("Replication pipeline completed successfully.\n")