
# This code reproduces empirical results for the time-varying estimates in:
# The Impact of Debt and Deficits on Long-Term Interest Rates in the US
# IMF WP, Furceri, Goncalves and Li, 2025


# packages
pacman::p_load(magrittr, lubridate, kableExtra, zoo, ggthemes, WDI, reshape, gtools,
               readxl, forcats, scales, hrbrthemes, readr, leaflet, ggplot2, ggrepel, readxl, EpiEstim, purrr, tidyr,
               sandwich, lmtest, RColorBrewer, wesanderson, ggpubr, rnaturalearth, sf,cowplot, googleway, data.table, wbstats,
               reshape, viridis, xtable, fredr, ggridges, dplyr)
pacman::p_load(tvthemes,urbnthemes,ggthemr,ggtech)
library(sandwich)
library(lmtest)
library(forecast)
library(tseries)
library(stargazer)
setwd("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/codes")
################################
# read data from directory

df_deficit = read_excel("github_def_pdf_debt_0110.xlsx")%>%na.exclude()

df_debt = read_excel("github_def_pdf_debt_0110.xlsx", sheet = "debt_pdf")%>%na.exclude()

################################################
###### Accute Crises, drop

df_deficit = df_deficit%>%filter(date!=as.Date("2008-09-01")&
                                   date!=as.Date("2009-01-01")&
                                   date!=as.Date("2020-09-01"))
df_debt = df_debt%>%filter(date!=as.Date("2008-09-01")&
                                   date!=as.Date("2009-01-01")&
                                   date!=as.Date("2020-09-01"))

####################
## add controls

controls=read_excel("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/codes/controls.xlsx")
df_deficit=merge(df_deficit, controls, by="date")
df_debt=merge(df_debt, controls, by="date")
df_deficit=df_deficit%>%
  filter(projected_deficit!="NA")
df_debt=df_debt%>%
  filter(projected_debt!="NA")

### add recession date dummies
rc=read_excel("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/Other datasources/RECESSIONS1.xlsx")
rc$date=as.Date(rc$date, "%Y-%m-%d")
df_debt$date=as.Date(df_debt$date, "%Y-%m-%d")
df_deficit$date=as.Date(df_deficit$date, "%Y-%m-%d")
df_deficit=merge(df_deficit, rc, by="date")
df_debt=merge(df_debt, rc, by="date")
acm= read_excel("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/Other datasources/ACMTermPremium.xlsx")
acm$date=as.Date(acm$date, "%Y-%m-%d")
df_deficit=merge(df_deficit, acm, by="date")
df_debt=merge(df_debt, acm, by="date")

library(fuzzyjoin)
fhold= read_excel("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/Other datasources/Foreignholdings.xlsx")
fhold$date=as.Date(fhold$date)
df_debt$date=as.Date(df_debt$date)
df_deficit=merge(df_deficit, fhold, by="date")
df_debt = 
  fuzzyjoin::difference_left_join(
  df_debt, fhold, by = "date",
  max_dist = 100, # Allow any distance
  distance_col = "date_diff")
df_debt = df_debt %>%  
 dplyr::rename(date = date.x) %>%  
  dplyr::select(-date.y)# Ensure consistent column name
df_debt = df_debt%>%
  group_by(date) %>%
  slice_min(abs(date_diff)) %>%  
  ungroup() %>%  # Ungroup after slicing
  dplyr::select(-date_diff) # Remove the distance column if not needed

pop = read_excel("C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/Other datasources/UN_populationestimates.xlsx")
pop$year=as.numeric(pop$date)
pop=pop%>%select(-date)
df_debt$year=df_debt$projected_fiscal_year-4
df_debt=merge(df_debt, pop, by="year")
df_deficit$year=df_deficit$projected_fiscal_year-4
df_deficit=merge(df_deficit, pop, by="year")

window_size <- 40 # our benchmark

##############
# DEBT as X
# Forward 5-10 as Y

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
    # Subset data for the current window
  subset_data <- df_debt[index, ]
  
  # Fit the linear model
  fit <- lm(svenf1015~projected_debt+
              tbill3m+growth+USREC+
              poly(projected_fiscal_year, 2), 
            data = subset_data)
beta <- summary(fit)$coefficients[2]
beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_debt), 
  width = window_size, 
  FUN = rolling_regression, 
  by.column = FALSE, 
  partial = FALSE ,
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_debt",
                         "se_debt")

# Combine with original data for visualization
final_debt <- 
cbind(df_debt[window_size:nrow(df_debt),], 
      result_df)
final_debt$coeff_debt_1plus=
final_debt$coeff_debt+final_debt$se_debt  
final_debt$coeff_debt_1minus=
final_debt$coeff_debt-final_debt$se_debt  

library(ggtext)
f1015debt=ggplot(subset(final_debt,year>2006), 
                 aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_debt))+
  geom_line(size=.7, aes(y=coeff_debt_1minus))+
  geom_line(size=.7, aes(y=coeff_debt_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
 # geom_hline(size=.7,yintercept = 0.026, color="red2")+
  theme_ridges()+
  labs(x="", y="Point coeff +/- 1SE",
       title="Y variable:F10-15",
subtitle="Coeff of Debt")+
  coord_cartesian(ylim = c(-0.02, 0.13))+
  theme(plot.caption = element_markdown(),
    plot.subtitle = element_markdown())
f1015debt

# DEBT as X
# T10 as Y
rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_debt[index, ]
  
  # Fit the linear model
  fit <- lm(t10~projected_debt+
              tbill3m+foreignholdings+growth+USREC+
              poly(projected_fiscal_year, 2), 
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_debt), 
  width = window_size, 
  FUN = rolling_regression, 
  by.column = FALSE, 
  partial = FALSE ,
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_debt",
                         "se_debt")

# Combine with original data for visualization
final_debt <- 
  cbind(df_debt[window_size:nrow(df_debt),], 
        result_df)
final_debt$coeff_debt_1plus=
  final_debt$coeff_debt+final_debt$se_debt  
final_debt$coeff_debt_1minus=
  final_debt$coeff_debt-final_debt$se_debt  

t10debt=ggplot(subset(final_debt,year>2006), 
                 aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_debt))+
  geom_line(size=.7, aes(y=coeff_debt_1minus))+
  geom_line(size=.7, aes(y=coeff_debt_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  # geom_hline(size=.7,yintercept = 0.026, color="red2")+
  theme_ridges()+
  labs(x="", y="Point coeff +/- 1SE",
       title="Y variable:T10",
       subtitle="Coeff of Debt")+
  coord_cartesian(ylim = c(-0.02, 0.13))+
  theme(plot.caption = element_markdown(),
        plot.subtitle = element_markdown())

# DEBT as X
# Tprem10 as Y
rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_debt[index, ]
  
  # Fit the linear model
  fit <- lm(acm10~projected_debt+
              tbill3m+foreignholdings+growth+USREC+
              poly(projected_fiscal_year, 2), 
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_debt), 
  width = window_size, 
  FUN = rolling_regression, 
  by.column = FALSE, 
  partial = FALSE ,
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_debt",
                         "se_debt")

# Combine with original data for visualization
final_debt <- 
  cbind(df_debt[window_size:nrow(df_debt),], 
        result_df)
final_debt$coeff_debt_1plus=
  final_debt$coeff_debt+final_debt$se_debt  
final_debt$coeff_debt_1minus=
  final_debt$coeff_debt-final_debt$se_debt  

ACM10debt=ggplot(subset(final_debt,year>2006), 
               aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_debt))+
  geom_line(size=.7, aes(y=coeff_debt_1minus))+
  geom_line(size=.7, aes(y=coeff_debt_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  # geom_hline(size=.7,yintercept = 0.026, color="red2")+
  theme_ridges()+
  labs(x="", y="Point coeff +/- 1SE",
       title="Y variable:ACM10",
       subtitle="Coeff of Debt")+
  coord_cartesian(ylim = c(-0.02, 0.13))+
  theme(plot.caption = element_markdown(),
        plot.subtitle = element_markdown())
ACM10debt

################################
# DEFICITS regs
#####################
# ACM10 as Y

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(acm10~projected_deficit+
              tbill3m+USREC+foreignholdings+growth+
              poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

acm10DEF=ggplot(subset(final_def, year>2006), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
#  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
    labs(x="", y="Point coeff +/- 1SE",title = "Y variable:ACM10",
       subtitle="Total Balances")
acm10DEF

# T10 as Y

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(t10~projected_deficit+
      tbill3m+USREC+foreignholdings+growth+
      poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

T10DEF=ggplot(subset(final_def, year>2006), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  #  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
  labs(x="", y="Point coeff +/- 1SE",
       title = "Y variable:T10",
       subtitle="Total Balances")

# F1015 as Y

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(svenf1015~projected_deficit+
              tbill3m+USREC+foreignholdings+growth+
              poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

F10DEF=ggplot(subset(final_def, year>2006), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  #  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
  labs(x="", y="Point coeff +/- 1SE",
       title = "Y variable:F10",
       subtitle="Total Balances")
F10DEF

ggarrange(F10DEF, T10DEF, acm10DEF, nrow = 1)
ggarrange(f1015debt, t10debt, ACM10debt, nrow = 1)

############################
# primary deficit regs ##
############################

# DEFICIT ACM10

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(acm10~projected_primary+
              tbill3m+USREC+foreignholdings+growth+
              poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

acm10PDEF=ggplot(subset(final_def, year>2006&year<2024), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  #  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
  labs(x="", y="Point coeff +/- 1SE",title = "Y variable:ACM10",
       subtitle="Primary Balances")
acm10PDEF

# T10 as Y

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(t10~projected_primary+
              tbill3m+USREC+foreignholdings+growth+
              poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

T10PDEF=ggplot(subset(final_def, year>2006&year<2024), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  #  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
  labs(x="", y="Point coeff +/- 1SE",
       title = "Y variable:T10",
       subtitle="Primary Balances")
T10PDEF

# DEFICIT F10

rolling_regression <- function(index) {
  if (length(index) < window_size) {
    return(c(NA, NA))  # Not enough data for regression
  }
  # Subset data for the current window
  subset_data <- df_deficit[index, ]
  
  # Fit the linear model
  fit <- lm(svenf1015~projected_primary+
              tbill3m+USREC+foreignholdings+growth+
              poly(projected_fiscal_year, 2),
            data = subset_data)
  beta <- summary(fit)$coefficients[2]
  beta_se <- summary(fit)$coefficients[2,"Std. Error"]
  # Return the coefficients
  return(c(beta, beta_se))
}

# Apply rolling regression
result <- rollapply(
  1:nrow(df_deficit), 
  width = window_size, 
  FUN = rolling_regression, 
  partial = FALSE ,
  by.column = FALSE, 
  align = "right"
)

# Convert the result to a data frame
result_df <- as.data.frame(result)
colnames(result_df) <- c("coeff_deficit",
                         "se_deficit")

# Combine with original data for visualization
final_def <- 
  cbind(df_deficit[window_size:nrow(df_deficit),], 
        result_df)
final_def$coeff_deficit_1plus=
  final_def$coeff_deficit+final_def$se_deficit  
final_def$coeff_deficit_1minus=
  final_def$coeff_deficit-final_def$se_deficit  

F10PDEF=ggplot(subset(final_def, year>2006&year<2024), aes(x=date))+
  geom_line(size=1.5, aes(y=coeff_deficit))+
  geom_line(size=.5, aes(y=coeff_deficit_1minus))+
  geom_line(size=.5, aes(y=coeff_deficit_1plus))+
  geom_hline(yintercept = 0, linetype="dashed")+
  #  geom_hline(size=.7,yintercept = -0.3, color="red2")+
  theme_ridges()+
  coord_cartesian(ylim = c(-0.75, 0.25))+
  labs(x="", y="Point coeff +/- 1SE",
       title = "Y variable:F10",
       subtitle="Primary Balances")
FP10DEF

ggarrange(F10PDEF, T10PDEF, acm10PDEF, nrow = 1)
