# 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)
#########################
## Q1 
#########################
# 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)
bstar <- numeric(n_periods)
bstar <- rep(0.6, n_periods)

# Initial values for b_t and d_t
b[1] <- 1.4  # 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
set.seed(123)
# 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 + R_av[t])) * b[t-1] #-((R_av[t] - g[t]) / (1 + g[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()

#########################
## Q2 part 1
#########################
# 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)
bstar <- numeric(n_periods)
bstar <- rep(0.6, n_periods)

# Initial values for b_t and d_t
b[1] <- 1.4  # 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] -0.1*(b[t-1]-bstar[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()

#########################
## Q2 part 2
#########################

# Parameters
mu_1 <- 0.03  # real rate 
mu_2 <- 0.01  # real gdp growth 
mu_3 <- 0.02  # inflation 

# Initialization
n_periods <- 100
R_r<- numeric(n_periods)
R_av<- numeric(n_periods)
g_r<- numeric(n_periods)
infl<- numeric(n_periods)
g<- numeric(n_periods)
sfa <- numeric(n_periods)
b <- numeric(n_periods)
d <- numeric(n_periods)
disc_b<- numeric(n_periods)
Y<-numeric(n_periods)
B<-numeric(n_periods)
# Initial values for b_t and d_t
b[1] <- 1.4  # 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
Y[1]<- 1
B[1]<- 1.4
# Simulate the series
for (t in 2:n_periods) {
  sfa[t] <- rnorm(1, mean = 0.00, sd = 0.01)
  R_r[t] <- mu_1
  g_r[t] <- mu_2
  infl[t]<-mu_3
  Y[t]<-Y[t-1]*(1+g_r[t]+infl[t])
  R_av[t]<-R_r[t]+infl[t] 
  
  g[t]<- g_r[t]+infl[t]
  d[t] <- -((R_av[t] ) / (1 + g[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[t] <-b[t]*Y[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,
  B_t=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 = B_t, color = "B_t"), linewidth = 1) +
  labs(title = "nominal debt",
       y = "Value", x = "Date") +
  scale_color_manual(name = "Series", values = c("B_t" = "green")) +
  theme_minimal()

#########################
## Q3 
#########################

# Define constants
rr <- 0.03
beta <- 2

# Define the range for pi_t
pi_t <- seq(0, 1, length.out = 100)

# Calculate y_t as a function of pi_t
y_t <- pi_t*exp(-beta * (rr + pi_t))

# Create the plot using ggplot2
plot <- ggplot(data = data.frame(pi_t, y_t), aes(x = pi_t, y = y_t)) +
  geom_line(color = "blue", linewidth = 1.2) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.8, linetype = "dashed") +  # Add zero line
  labs(x = expression(paste("Tax Rate")),
       y = expression(paste("Tax Revenue as % of GDP")),
       title = expression(paste("The Laffer Curve for Inflation Tax"))) +
  theme_minimal()

# Save the plot as a PDF file
ggsave("Laffer_Curve_Inflation_Tax.pdf", plot = plot, width = 8, height = 6)


############################
# The Debt Dynamics 
############################

# Parameters
mu_1 <- 0.03  # real rate 
mu_2 <- 0.01  # real gdp growth 
mu_31 <- 0.03  # money growth first period 
mu_32 <- 0.06  # money growth after policy switch 

# Initialization
n_periods <- 100
R_r<- numeric(n_periods)
R_av<- numeric(n_periods)
g_r<- numeric(n_periods)
infl<- numeric(n_periods)
mg<- numeric(n_periods)
g<- numeric(n_periods)
sfa <- numeric(n_periods)
itb <- numeric(n_periods)
b <- numeric(n_periods)
d <- 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] <- 1.4  # 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
itb[1] <- 1

# Simulate the series
for (t in 2:n_periods) {
  sfa[t] <- rnorm(1, mean = 0.00, sd = 0.01)
  R_r[t] <- mu_1
  g_r[t] <- mu_2
  # Set money growth conditionally 
  if (t <= 20) {
    mg[t] <- mu_31 
  } else {
    mg[t] <- mu_32
  }
  infl[t]<-mg[t]-g_r[t]
  R_av[t]<-R_r[t]+infl[t] 
  g[t]<- g_r[t]+infl[t]
  itb[t]<-exp(-beta *  R_av[t])
  d[t] <- 0.0125 #-((R_av[t] - g[t]) / (1 + g[t])) * b[t-1] 
  b[t] <- ((1 + R_av[t]) / (1 + g[t])) * b[t-1] -mg[t]*itb[t] + 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()