# 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)
}


### Blanchard Replication---------------------------------------------------
num_sim <- 10000
n_periods <- 10
c_list <- c(0,0.2,0.5)
b_mat_list <- list()
x_mat_list <- list()
for (c in c_list){
  b_mat <- matrix(nrow = num_sim, ncol = n_periods)
  x_mat <- matrix(nrow = num_sim, ncol = n_periods)
  for (i in 1:num_sim){
    s_x = 0.003
    s_u = 0.01
    s_s = 0.01 #Not given
    a_u = -0.02
    a_s = -0.02 #a_u = a_s, implies that under certainty, the debt ratio would remain constant.
    
    e_x = numeric(n_periods)
    e_u = numeric(n_periods)
    x = numeric(n_periods)
    u = numeric(n_periods)
    rg_diff = numeric(n_periods)
    s <- numeric(n_periods)
    b <- numeric(n_periods)
    d <- numeric(n_periods)
    
    b[1] = 1
    x[1] = 0
    rg_diff[1] = -0.02
    for (t in 2:n_periods){
      e_x[t] = rnorm(1,mean = 0, sd = s_x)
      x[t] = x[t-1] + e_x[t]
      e_u[t] = rnorm(1, mean = 0, sd = s_u)
      u[t] = a_u + e_u[t]
      rg_diff[t] <- x[t] + u[t]
      e_s <- rnorm(1, mean = 0, sd = s_s)
      #s[t] <- a_s + e_s + c*rg_diff[t]*b[t-1]
      s[t] <- (1-c)*(a_s + e_s) + c*rg_diff[t]*b[t-1]
      b[t] <- b[t-1] + rg_diff[t]*b[t-1] - s[t]
    }
    b_mat[i,] <- b
    x_mat[i,] <- x
  }
  b_mat_list <- list.append(b_mat_list, b_mat)
  x_mat_list <- list.append(x_mat_list, x_mat)
  
}


# 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)

plot_list <- list()
for (b_mat in b_mat_list){
  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 data frame for plotting
  data <- data.frame(
    Date = dates,
    b_t_mean = b_mean,
    b_t_lb = b_lb,
    b_t_ub = b_ub
  )
  
  p <- 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")
  
  plot_list <- list.append(plot_list,p)
}
# Save the plot as a PDF
pdf("debt_ratio_plot.pdf", width = 7, height = 5)  # Adjust width and height as needed
print(p)  # Print the plot inside the PDF
dev.off()



plot_list1 <- list()
for (x_mat in x_mat_list){
  x_mean <- apply(x_mat, 2, mean)
  x_sd <- apply(x_mat,2, sd)
  x_ub <- x_mean + 2*x_sd
  x_lb <- x_mean - 2*x_sd
  
  # Create a data frame for plotting
  data <- data.frame(
    Date = dates,
    x_t_mean = x_mean,
    x_t_lb = x_lb,
    x_t_ub = x_ub
  )
  
  p <- ggplot(data, aes(x = Date)) +
    geom_ribbon(aes(ymin = x_t_lb, ymax = x_t_ub),  fill = "lightblue", alpha = 0.5) +
    geom_line(aes(y = x_t_mean), size = 0.5, col = "blue") +
    labs(x = "Date", y = "r-g ", title = "Mean r-g with Upper and Lower Bounds Over Time")
  
  plot_list1 <- list.append(plot_list1,p)
}

#-------------------
#density plots
#-------------------

b_mat <- b_mat_list[[1]]
c_0 <- data.frame(deviation = b_mat[,n_periods] - b_mat[,1], cat = rep("c_0",num_sim))
b_mat <- b_mat_list[[2]]
c_02 <-  data.frame(deviation = b_mat[,n_periods] - b_mat[,1], cat = rep("c_02",num_sim))
b_mat <- b_mat_list[[3]]
c_05 <- data.frame(deviation = b_mat[,n_periods] - b_mat[,1], cat = rep("c_05",num_sim))
debt_exp <-  do.call('rbind', list(c_0,c_02,c_05))
p <- ggplot(debt_exp, aes(deviation, fill = cat)) +
  geom_density(alpha = 0.7)

# Save the plot as a PDF
pdf("SDSA_fig_1.pdf", width = 7, height = 5)  # Adjust width and height as needed
print(p)  # Print the plot inside the PDF
dev.off()
# check the quantile value that it iexceeded with probability 2 per cent 
quantile_value_c_0 <- quantile(c_0$deviation, probs = 0.98, na.rm = TRUE)
quantile_value_c_02 <- quantile(c_02$deviation, probs = 0.98, na.rm = TRUE)
quantile_value_c_05 <- quantile(c_05$deviation, probs = 0.98, na.rm = TRUE)
print(quantile_value_c_02)