# Programme to simulate  the dynamics of debt/GDP  with the simplest model
# Traditional framework 
# Carlo A. Favero November 2024

rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
listofpackages <- c("ggplot2")

for (j in listofpackages){
  if (sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = TRUE)
}

# Load necessary libraries
library(ggplot2)

# Parameters
mu_1 <- 0.05  # Example value for R^{av}_t (mu_1)
mu_2 <- 0.03  # Example value for g_t (mu_2)

# Initialization
n_periods <- 100
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)


# Initial values for b_t and d_t
b[1] <- 1.2  # Arbitrary starting value for b_1
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

# Simulate the series
for (t in 2:n_periods) {
  R_av[t] <- mu_1
  g[t] <- mu_2
  sfa[t] <- rnorm(1, mean = 0.00, sd = 0.01)
  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
}


# 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 = b,
  d_t = d,
  disc_b_t=disc_b
)

# Plot the series using ggplot2
ggplot(data, aes(x = Date)) +
  geom_line(aes(y = b_t, color = "b_t"), linewidth = 1) +
  labs(title = "Simulation of b_t over the horizon",
       y = "Value", x = "Date") +
  scale_color_manual(name = "Series", values = c("b_t" = "blue")) +
  theme_minimal()

ggplot(data, aes(x = Date)) +
  geom_line(aes(y = d_t, color = "d_t"), linewidth = 1) +
  labs(title = "Simulation of d_t over the horizon",
       y = "Value", x = "Date") +
  scale_color_manual(name = "Series", values = c("d_t" = "red")) +
  theme_minimal()

ggplot(data, aes(x = Date)) +
  geom_line(aes(y = disc_b_t, color = "disc_b_t"), linewidth = 1) +
  labs(title = "discounted value of future d_t over the horizon",
       y = "Value", x = "Date") +
  scale_color_manual(name = "Series", values = c("disc_b_t" = "green")) +
  theme_minimal()