# Clear environment and set working directory
rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# Load necessary package
listofpackages <- c("ggplot2","dplyr","readxl","rlist")
for (j in listofpackages) {
  if (sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = TRUE)
}



### SDSA with bootstrap applied to Italian Data-------------------------
IT_DB <- read_xlsx("data/DB.xlsx") %>% filter(MS == "IT (Italy)")
IT_DB$Debt_lag <- lag(IT_DB$`Gvt Debt`)
IT_DB$Avg_cost <- (IT_DB$`Primary Surplus` - IT_DB$Surplus)/IT_DB$Debt_lag
IT_DB$sfa <- IT_DB$SFA/IT_DB$GDP
names(IT_DB)[4] <- "Debt"
names(IT_DB)[8] <- "Growth"
mod_growth <- lm(Growth ~ 1, data = IT_DB)
mu_1 <- coef(mod_growth)
resid_g <- resid(mod_growth)[-1]
mod_r <- lm(Avg_cost ~ lag(Avg_cost), data = IT_DB)
resid_r <- resid(mod_r)
mu_2 <- coef(mod_r)[1]
rho <- coef(mod_r)[2]
mod_sfa <- lm(sfa ~ 1, data = IT_DB)
mu_3 <- coef(mod_sfa)
resid_sfa <- resid(mod_sfa)[c(-1,-2)]
e_t <- cbind(resid_g,resid_r,resid_sfa)


num_sim <- 10000
n_periods <- 10
b_mat <- matrix(nrow = num_sim, ncol = n_periods)
for (i in 1:num_sim){
  # Initialization
  R_av<- numeric(n_periods)
  g<- numeric(n_periods)
  sfa <- numeric(n_periods)
  b <- numeric(n_periods)
  d <- numeric(n_periods)
  sfa <- numeric(n_periods)
  disc_b<- numeric(n_periods)
  bstar <- numeric(n_periods)
  bstar <- rep(0.6, n_periods)
  
  # Initial values for b_t and d_t
  b[1] <- IT_DB$Debt[nrow(IT_DB)]/IT_DB$GDP[nrow(IT_DB)]
  disc_b[1]<-b[1] # discounted b set equal to b in the first period 
  d[1] <- -0.03    # Initial value for d_1 (could be different)
  sfa[1] <- 0
  R_av[1] <- IT_DB$Avg_cost[nrow(IT_DB)]
  for (t in 2:n_periods) {
    e_bs <- e_t[sample(1:nrow(e_t),1),]
    g[t] <- mu_1 + e_bs[1]
    R_av[t] <- mu_2 + rho*R_av[t-1] + e_bs[2]
    sfa[t] <- mu_3 + e_bs[3]
    d[t] <- -((R_av[t] - g[t]) / (1 + g[t])) * b[t-1] #-((R_av[t] - g[t]) / (1 + R_av[t])) * b[t-1] 
    b[t] <- ((1 + R_av[t]) / (1 + g[t])) * b[t-1] + d[t] + sfa[t]
    disc_b[t]<- b[t]*((1 +  g[t]) / (1 + R_av[t]))^t
  }
  b_mat[i,] <- b
}

b_mean <- apply(b_mat, 2, mean)
b_sd <- apply(b_mat,2, sd)
b_ub <- b_mean + 2*b_sd
b_lb <- b_mean - 2*b_sd

# Create a date series starting from 2024
start_date <- as.Date("2024-01-01")
dates <- seq.Date(from = start_date, by = "year", length.out = n_periods)

# Create a data frame for plotting
data <- data.frame(
  Date = dates,
  b_t_mean = b_mean,
  b_t_lb = b_lb,
  b_t_ub = b_ub
)

ggplot(data, aes(x = Date)) +
  geom_ribbon(aes(ymin = b_t_lb, ymax = b_t_ub),  fill = "lightblue", alpha = 0.5) +
  geom_line(aes(y = b_t_mean), size = 0.5, col = "blue") +
  labs(x = "Date", y = "Debt Ratio", title = "Mean Debt Ratio with Upper and Lower Bounds Over Time")

debt_exp <- b_mat[,n_periods] - b_mat[,1]
hist(debt_exp, freq = FALSE, breaks = 50)