# === 加载库与参数 ===
library(data.table)
library(forecast)
set.seed(12345)
alpha <- 0.1
r <- 0.5

LRBR <- function(S, q, r) {
  u <- rbinom(1, 1, r)
  v <- rbinom(1, 1, 0.5)
  if (u == 1) {
    return(as.numeric(q > S)) 
  } else {
    return(v)
  }
}
# The file is located in the folder and named "electricity-normalized.csv".
# Example: file_path <- "D:/path/to/your/folder/electricity-normalized.csv"
# =====================================
file_path <- "D:/electricity-normalized.csv" 
data <- fread(file_path)
y_vals <- data$nswdemand  
n <- length(y_vals)

# Generate AR(3) predictions
ar3_model <- Arima(y_vals, order = c(3, 0, 0)) 
y_hat <- fitted(ar3_model) 


s <- numeric(n); s[1] <- 0
W <- numeric(n)
lambda <- numeric(n); lambda[1] <- 0
coverage <- numeric(n)
rolling_coverage <- numeric(n)
rolling_interval_length <- numeric(n)
longrun_coverage <- numeric(n)
longrun_interval_width <- numeric(n)


for (i in 1:n) {
  S_t <- abs(y_vals[i] - y_hat[i])
  L <- LRBR(S_t, s[i], r)
  g_t <- ifelse(L == 1, 1 - (r * (1 - alpha) + (1 - r) * 0.5), 
                -(r * (1 - alpha) + (1 - r) * 0.5))
  
  W[i] <- ifelse(i == 1, 1 - g_t * s[1], W[i-1] - g_t * s[i])
  lambda[i + 1] <- (i / (i + 1)) * lambda[i] - (1 / (i + 1)) * g_t
  s[i + 1] <- lambda[i + 1] * W[i]
  
  coverage[i] <- as.numeric(y_vals[i] >= (y_hat[i] - s[i]) & y_vals[i] <= (y_hat[i] + s[i]))
  rolling_coverage[i] <- mean(coverage[max(1, i - 199):i])
  rolling_interval_length[i] <- mean(2 * s[max(1, i - 199):i])
  longrun_coverage[i] <- mean(coverage[1:i])
  longrun_interval_width[i] <- mean(2 * s[1:i])
}





output_data <- data.frame(
  Index = 1:n,
  Coverage = coverage,
  Rolling_Coverage = rolling_coverage,
  Longrun_Coverage = longrun_coverage,
  Rolling_Interval_Width = rolling_interval_length,
  Longrun_Interval_Width = longrun_interval_width
)

# 保存
output_dir <- "results"
dir.create(output_dir, showWarnings = FALSE)
# Modify the directory where the output file will be saved.
# Example: setwd("D:/path/to/your/folder")
# =====================================
setwd("D:")

filename <- paste0("elec2", r, ".csv")
filepath <- file.path(output_dir, filename)

write.csv(output_data, filepath, row.names = FALSE)
cat("✅ 已保存至：", filepath, "\n")

