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

# packages needed
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()
### Either the numbers put forth for the variable itself or the constructed based on fiscal flows
df_debt = read_excel("github_def_pdf_debt_0110.xlsx", sheet = "debt_pdf")%>%na.exclude()


###### DROPPING accute crises obs

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=c("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") %>%  
  rename(date = date.x) %>%  # Ensure consistent column name
  group_by(date) %>%
  slice_min(abs(date_diff)) %>%  
  ungroup() %>%  # Ungroup after slicing
  dplyr::select(-date_diff) # Remove the distance column if not needed

# add population growth
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")
write.csv(df_debt, 'debtregressionsZIO.csv')
write.csv(df_deficit, 'deficitregressionsZIO.csv')

#################################################
## create sub-samples ########################### 
#################################################

df_debt_precovid = subset(df_debt, date<as.Date("2020-02-01"))
df_deficit_precovid = subset(df_deficit, date<as.Date("2020-02-01"))
df_debt_preGFC = subset(df_debt, date<as.Date("2008-02-01"))
df_deficit_preGFC = subset(df_deficit, date<as.Date("2008-02-01"))

###############################################
######## PRE GFC REGS #########################
###############################################

m1a=lm(data = df_deficit_preGFC, 
t10 ~ projected_primary+tbill3m+
  growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m1b=lm(data = df_deficit_preGFC, 
svenf510 ~ projected_primary+tbill3m+
  growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m1c=lm(data = df_deficit_preGFC, 
svenf1015 ~ projected_primary+tbill3m+growth+
  foreignholdings+USREC+growth_pop+
  poly(projected_fiscal_year,2))

m2a=lm(data = df_deficit_preGFC, 
t10 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m2b=lm(data = df_deficit_preGFC, 
svenf510 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m2c=lm(data = df_deficit_preGFC, 
svenf1015 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))

m3a=lm(data = df_debt_preGFC, 
t10 ~ projected_debt+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3b=lm(data = df_debt_preGFC, 
svenf510 ~ projected_debt+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3c=lm(data = df_debt_preGFC, 
svenf1015 ~ projected_debt+tbill3m+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m3a)

### SE corrections
nw_se_m1a <- sqrt(diag(NeweyWest(m1a)))
nw_se_m1b <- sqrt(diag(NeweyWest(m1b)))
nw_se_m1c <- sqrt(diag(NeweyWest(m1c)))
nw_se_m2a <- sqrt(diag(NeweyWest(m2a)))
nw_se_m2b <- sqrt(diag(NeweyWest(m2b)))
nw_se_m2c <- sqrt(diag(NeweyWest(m2c)))
nw_se_m3a <- sqrt(diag(NeweyWest(m3a)))
nw_se_m3b <- sqrt(diag(NeweyWest(m3b)))
nw_se_m3c <- sqrt(diag(NeweyWest(m3c)))

nw_se <- lapply(list(m1a, m1b, m1c, 
  m2a, m2b, m2c), function(model) {
  sqrt(diag(NeweyWest(model)))
})

# Generate LaTeX table with Newey-West standard errors
stargazer(m1a, m1b, m1c, m2a, m2b, m2c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Balances: Pre-GFC period",
          font.size = "small",
          column.sep.width = "1pt",
          omit=c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

nw_se <- lapply(list(m3a, m3b, m3c), 
                function(model) {
      sqrt(diag(NeweyWest(model)))
                     })
stargazer(m3a, m3b, m3c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Debt: Pre-GFC period",
          font.size = "small",
          omit = c("poly", "Constant"),
          column.sep.width = "1pt",
        se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

######### Regressions using term premia ###########
############### ACM ################################

m1a=lm(data = df_deficit_preGFC, 
       acm5 ~ projected_primary+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m1b=lm(data = df_deficit_preGFC, 
       acm10 ~ projected_primary+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))

m2a=lm(data = df_deficit_preGFC, 
       acm5 ~ projected_deficit+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m2b=lm(data = df_deficit_preGFC, 
       acm10 ~ projected_deficit+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3a=lm(data = df_debt_preGFC, 
       acm5 ~ projected_debt+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3b=lm(data = df_debt_preGFC, 
       acm10 ~ projected_debt+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
### SE corrections
nw_se_m1a <- sqrt(diag(NeweyWest(m1a)))
nw_se_m1b <- sqrt(diag(NeweyWest(m1b)))
nw_se_m2a <- sqrt(diag(NeweyWest(m2a)))
nw_se_m2b <- sqrt(diag(NeweyWest(m2b)))
nw_se_m3a <- sqrt(diag(NeweyWest(m3a)))
nw_se_m3b <- sqrt(diag(NeweyWest(m3b)))

nw_se <- lapply(list(m1a, m1b,
                     m2a, m2b, m3a, m3b), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })

# Generate LaTeX table with Newey-West standard errors
stargazer(m1a, m1b, m2a, m2b,m3a, m3b,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Term premia: Full Sample",
          font.size = "small",
          column.sep.width = "1pt",
          omit=c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)


#####################################################
######## PRE COVID ##################################
#####################################################

m1a=lm(data = df_deficit_precovid, 
       t10 ~ projected_primary+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))
m1b=lm(data = df_deficit_precovid, 
       svenf510 ~ projected_primary+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))
m1c=lm(data = df_deficit_precovid, 
       svenf1015 ~ projected_primary+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))

m2a=lm(data = df_deficit_precovid, 
       t10 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC
       +growth_pop+poly(projected_fiscal_year,2))
m2b=lm(data = df_deficit_precovid, 
       svenf510 ~ projected_deficit+tbill3m+growth+foreignholdings
         +growth_pop+USREC+poly(projected_fiscal_year,2))
m2c=lm(data = df_deficit_precovid, 
       svenf1015 ~ projected_deficit+tbill3m+growth+foreignholdings
       +USREC+growth_pop+poly(projected_fiscal_year,2))

m3a=lm(data = df_debt_precovid, 
       t10 ~ projected_debt+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))
m3b=lm(data = df_debt_precovid, 
       svenf510 ~ projected_debt+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))
m3c=lm(data = df_debt_precovid, 
       svenf1015 ~ projected_debt+tbill3m+growth+foreignholdings+
         USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m3c)

### SE corrections
nw_se_m1a <- sqrt(diag(NeweyWest(m1a)))
nw_se_m1b <- sqrt(diag(NeweyWest(m1b)))
nw_se_m1c <- sqrt(diag(NeweyWest(m1c)))
nw_se_m2a <- sqrt(diag(NeweyWest(m2a)))
nw_se_m2b <- sqrt(diag(NeweyWest(m2b)))
nw_se_m2c <- sqrt(diag(NeweyWest(m2c)))
nw_se_m3a <- sqrt(diag(NeweyWest(m3a)))
nw_se_m3b <- sqrt(diag(NeweyWest(m3b)))
nw_se_m3c <- sqrt(diag(NeweyWest(m3c)))

nw_se <- lapply(list(m1b, m1c,m1a, 
                     m2b, m2c, 
                     m2a), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })

# Generate LaTeX table with Newey-West standard errors
stargazer(m1b, m1c, m1a, m2b, m2c, m2a,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/precovid.tex",
          title = "Balances: Pre-covid period",
          font.size = "small",
          omit = c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

nw_se <- lapply(list(m3b, m3c, m3a), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })
stargazer(m3b, m3c,m3a,
          type = "latex",
          omit = c("poly", "Constant"),
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/precovid.tex",
          title = "Debt: Pre-covid period",
          font.size = "small",
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

############### PRECOVID ACM ################################
m1a=lm(data = df_deficit_precovid, 
       acm5 ~ projected_primary+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m1b=lm(data = df_deficit_precovid, 
       acm10 ~ projected_primary+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))

m2a=lm(data = df_deficit_precovid, 
       acm5 ~ projected_deficit+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m2b=lm(data = df_deficit_precovid, 
       acm10 ~ projected_deficit+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3a=lm(data = df_debt_precovid, 
       acm5 ~ projected_debt+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
m3b=lm(data = df_debt_precovid, 
       acm10 ~ projected_debt+growth+foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
### SE corrections
nw_se_m1a <- sqrt(diag(NeweyWest(m1a)))
nw_se_m1b <- sqrt(diag(NeweyWest(m1b)))
nw_se_m2a <- sqrt(diag(NeweyWest(m2a)))
nw_se_m2b <- sqrt(diag(NeweyWest(m2b)))
nw_se_m3a <- sqrt(diag(NeweyWest(m3a)))
nw_se_m3b <- sqrt(diag(NeweyWest(m3b)))

nw_se <- lapply(list(m1a, m1b,
                     m2a, m2b, m3a, m3b), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })

# Generate LaTeX table with Newey-West standard errors
stargazer(m1a, m1b, m2a, m2b,m3a, m3b,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Term premia: Full Sample",
          font.size = "small",
          column.sep.width = "1pt",
          omit=c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)


############################################
### Full Sample ############################
############################################

m1a=lm(data = df_deficit, 
       t10 ~ projected_primary+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m1b=lm(data = df_deficit, 
       svenf510 ~ projected_primary+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m1c=lm(data = df_deficit, 
       svenf1015 ~ projected_primary+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m1d = lm(data = df_deficit, 
         acm10 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))

m2a=lm(data = df_deficit, 
       t10 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m2b=lm(data = df_deficit, 
       svenf510 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m2c=lm(data = df_deficit, 
       svenf1015 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m2d = lm(data = df_deficit, 
         acm10 ~ projected_deficit+tbill3m+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))

m3a=lm(data = df_debt, 
t10 ~ projected_debt+tbill3m+growth_pop+
growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
summary(m3a)
m3b=lm(data = df_debt, 
       svenf510 ~ projected_debt+tbill3m+growth_pop+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m3c=lm(data = df_debt, 
       svenf1015 ~ projected_debt+tbill3m+growth_pop+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))
m3d = lm(data = df_debt, 
         acm10 ~ projected_debt+tbill3m+growth_pop+growth+foreignholdings+USREC+poly(projected_fiscal_year,2))

### SE corrections
nw_se_m1a <- sqrt(diag(NeweyWest(m1a)))
nw_se_m1b <- sqrt(diag(NeweyWest(m1b)))
nw_se_m1c <- sqrt(diag(NeweyWest(m1c)))
nw_se_m2a <- sqrt(diag(NeweyWest(m2a)))
nw_se_m2b <- sqrt(diag(NeweyWest(m2b)))
nw_se_m2c <- sqrt(diag(NeweyWest(m2c)))
nw_se_m3a <- sqrt(diag(NeweyWest(m3a)))
nw_se_m3b <- sqrt(diag(NeweyWest(m3b)))
nw_se_m3c <- sqrt(diag(NeweyWest(m3c)))

nw_se <- lapply(list(m1a, m1b, m1c, 
                     m2a, m2b, m2c), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })

# Generate LaTeX table with Newey-West standard errors
stargazer(m1a, m1b, m1c, m2a, m2b, m2c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Balances: Full sample",
          font.size = "small",
          omit = c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

nw_se <- lapply(list(m3a, m3b, m3c), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })
stargazer(m3a, m3b, m3c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Debt: Full sample",
          font.size = "small",
          omit = c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

nw_se <- lapply(list(m1d, m2d, m3d), function(model) {
  sqrt(diag(NeweyWest(model)))
})
stargazer(m1d, m2d, m3d,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "ACM",
          font.size = "small",
          omit = c("poly", "Constant"),
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

### ADD CONTROLS ####################################
#####################################################
####  DEBT  REGS ####################################

m5b=lm(data = df_debt_precovid, 
       svenf1015 ~ projected_debt+tbill3m+growth+
       foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5b)
m5a=lm(data = df_debt_precovid, 
       svenf510 ~ projected_debt+tbill3m+growth+
         foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5a)
m5c=lm(data = df_debt_precovid, 
       t10 ~ projected_debt+tbill3m+growth+
         foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5c)

model=c(m5a, m5b, m5c)
nw_se <- lapply(list(m5a, m5b, m5c), function(model) {
  sqrt(diag(NeweyWest(model)))
})
stargazer(m5a, m5b, m5c,
          font.size = "small",
          title = "Extra Controls",
          type = "latex",  # Change to "latex" or "html" as needed
          keep.stat = c("rsq", "n"),  # Keep only R-squared and number of observations
          column.sep.width = "0.5pt")


####  TOTAL BALANCE  #########################################

m5b=lm(data = df_deficit, 
       svenf1015 ~ projected_deficit+tbill3m+growth+
         foreignholdings+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5b)
m5a=lm(data = df_deficit, 
       svenf510 ~ projected_deficit+tbill3m+growth+
         foreignholdings+ USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5a)
m5c=lm(data = df_deficit, 
       t10 ~ projected_deficit+tbill3m+growth+
         foreignholdings+ USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5c)

model=c(m5a, m5b, m5c)
nw_se <- lapply(list(m5a, m5b, m5c), function(model) {
  sqrt(diag(NeweyWest(model)))
})
stargazer(m5a, m5b, m5c,
          font.size = "small",
          title = "Extra Controls",
          type = "latex",  # Change to "latex" or "html" as needed
          keep.stat = c("rsq", "n"),  # Keep only R-squared and number of observations
          column.sep.width = "0.5pt")

###########################################################
####  primary BALANCE  ####################################

m5b=lm(data = df_deficit, 
       svenf1015 ~ projected_primary+tbill3m+growth+
         foreignholdings+ USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5b)
m5a=lm(data = df_deficit, 
       svenf510 ~ projected_primary+tbill3m+growth+
         foreignholdings+ ebp+USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5a)
m5c=lm(data = df_deficit, 
       t10 ~ projected_primary+tbill3m+growth+
         foreignholdings+ USREC+growth_pop+poly(projected_fiscal_year,2))
summary(m5c)

model=c(m5a, m5b, m5c)
nw_se <- lapply(list(m5a, m5b, m5c), function(model) {
  sqrt(diag(NeweyWest(model)))
})
stargazer(m5a, m5b, m5c,
          font.size = "small",
          title = "Extra Controls",
          type = "latex",  # Change to "latex" or "html" as needed
          keep.stat = c("rsq", "n"),  # Keep only R-squared and number of observations
          column.sep.width = "0.5pt")


#### TERM PREMIA ###

acm10=lm(data = df_deficit, 
      acm10 ~ projected_deficit+foreignholdings+ growth+
        USREC+growth_pop+
      poly(projected_fiscal_year,2))
acm510=lm(data = df_deficit, 
        acm510 ~ projected_deficit+foreignholdings+ growth+
         growth+USREC+growth_pop+
           poly(projected_fiscal_year,2))
acm10b=lm(data = df_deficit, 
        acm10 ~ 
  projected_primary+foreignholdings+ growth+USREC+growth_pop+
  poly(projected_fiscal_year,2))
acm510b=lm(data = df_deficit, 
        acm510 ~ projected_primary+foreignholdings+ growth+USREC+growth_pop+
        poly(projected_fiscal_year,2))
acm10c=lm(data = df_debt, 
          acm10 ~ projected_debt+foreignholdings+ growth+
            USREC+growth_pop+
            poly(projected_fiscal_year,2))
acm510c=lm(data = df_debt, 
           acm510 ~ projected_debt+foreignholdings+growth+USREC+growth_pop+
            poly(projected_fiscal_year,2))
### SE corrections
nw_se_acm10 <- sqrt(diag(NeweyWest(acm10)))
nw_se_acm510 <- sqrt(diag(NeweyWest(acm510)))
nw_se_acm10b <- sqrt(diag(NeweyWest(acm10b)))
nw_se_acm510b <- sqrt(diag(NeweyWest(acm510b)))
nw_se_acm10c <- sqrt(diag(NeweyWest(acm10c)))
nw_se_acm510c <- sqrt(diag(NeweyWest(acm510c)))

nw_se <- lapply(list(acm10,acm510, acm10b, acm510b, 
acm10c, acm510c), function(model) {
 sqrt(diag(NeweyWest(model)))
})

# Generate LaTeX table with Newey-West standard errors
stargazer(acm10, acm510, acm10b, acm510b, acm10c, acm510c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Term Premia: Full sample",
          font.size = "small",
           # Keep only R-squared and number of observations
          column.sep.width = "0.5pt",
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)

####################################################################
#### EXCLUDING COVID ###############################################

df_debt_precovid = subset(df_debt, date.x<as.Date("2020-02-01"))
df_deficit_precovid = subset(df_deficit, date.x<as.Date("2020-02-01"))
df_debt_preGFC = subset(df_debt, date.x <as.Date("2008-02-01"))
df_deficit_preGFC = subset(df_deficit, date.x <as.Date("2008-02-01"))

acm10=lm(data = df_deficit_preGFC, 
         acm10 ~ projected_deficit+foreignholdings+ ebp+USREC+growth_pop+
           poly(projected_fiscal_year,2))
acm510=lm(data = df_deficit_preGFC, 
          acm510 ~ projected_deficit+foreignholdings+ ebp+USREC+growth_pop+
            poly(projected_fiscal_year,2))
acm10b=lm(data = df_deficit_preGFC, 
          acm10 ~ projected_primary+foreignholdings+ ebp+USREC+growth_pop+
            poly(projected_fiscal_year,2))
acm510b=lm(data = df_deficit_preGFC, 
           acm510 ~ projected_primary+foreignholdings+ ebp+USREC+growth_pop+
             poly(projected_fiscal_year,2))
acm10c=lm(data = df_debt_preGFC, 
          acm10 ~ projected_debt+foreignholdings+ ebp+USREC+growth_pop+
            poly(projected_fiscal_year,2))
acm510c=lm(data = df_debt_preGFC, 
           acm510 ~ projected_debt+foreignholdings+ ebp+USREC+growth_pop+
             poly(projected_fiscal_year,2))
### SE corrections
nw_se_acm10 <- sqrt(diag(NeweyWest(acm10)))
nw_se_acm510 <- sqrt(diag(NeweyWest(acm510)))
nw_se_acm10b <- sqrt(diag(NeweyWest(acm10b)))
nw_se_acm510b <- sqrt(diag(NeweyWest(acm510b)))
nw_se_acm10c <- sqrt(diag(NeweyWest(acm10c)))
nw_se_acm510c <- sqrt(diag(NeweyWest(acm510c)))

nw_se <- lapply(list(acm10,acm510, acm10b, acm510b, 
                     acm10c, acm510c), function(model) {
                       sqrt(diag(NeweyWest(model)))
                     })

# Generate LaTeX table with Newey-West standard errors
stargazer(acm10, acm510, acm10b, acm510b, acm10c, acm510c,
          type = "latex",
          out="C:/Users/cgoncalves/OneDrive - International Monetary Fund (PRD)/Debt and Interest Rates (with Hongchi and Furceri)/preGFC.tex",
          title = "Term Premia: Full sample",
          font.size = "small",
          # Keep only R-squared and number of observations
          column.sep.width = "0.5pt",
          se = nw_se,  # Use Newey-West standard errors
          keep.stat = c("n", "rsq"),  # Keep only number of observations and R²
          single.row = FALSE)
