library(dplyr)
library(tidyr)
library(readr)
library(panelAR)
library(lmtest)
library(stargazer)
library(formula.tools)
source("myFunctions.R")

source("loadData.R")

f1 <- as.formula("CO2 ~ GDP + URBAN + MIL + ID + factor(YR)")
f2 <- as.formula("CO2 ~ GDP + URBAN + MIL + MIL:GDP + ID + factor(YR)")
f3 <- as.formula("CO2 ~ GDP + URBAN + ID + MIL + YR:MIL + factor(YR)")
f4 <- as.formula("CO2 ~ GDP + URBAN + ID + MIL + MIL:GDP + YR:MIL + factor(YR) ")


cNames <- sort(unique(as.character(df$ID)))
recordMatrix <- as.data.frame(matrix(ncol=6, nrow=45, NA))
names(recordMatrix) <- c("yr2", "coefficient", "pcse", "tvalue", "pvalue", "ncountries")

i <- 1
for(y in 2004:1960){
  
  newdf <- subset(df_ub, YR >= y)
  recordMatrix$ncountries[i] <- length(unique(newdf$ID))
  newdf$ID <- factor(newdf$ID)
  p <- pw.function(f1, data=newdf)
  mil_ind <- which(names(p$coefficients)=="MIL")
  
  recordMatrix$yr2[i] <- y
  recordMatrix$coefficient[i] <- p$coefficients[mil_ind]
  recordMatrix$pcse[i] <- p$pcses[mil_ind]
  recordMatrix$tvalue[i] <- p$tvals[mil_ind]
  recordMatrix$pvalue[i] <- p$pvals[mil_ind]
  
i <- i+1
}


##Now go forward in time from 1970 to 2014
# recordMatrixR <- as.data.frame(matrix(ncol=6, nrow=50, NA))
# names(recordMatrixR) <- c("yrR", "coefficientR", "pcseR", "tvalueR", "pvalueR", "ncountriesR")

i <- 1
for(y in 1970:2014){
  
  newdf <- subset(df_ub, YR <= y)
  recordMatrix$ncountriesR[i] <- length(unique(newdf$ID))
  newdf$ID <- factor(newdf$ID)
  p <- pw.function(f1, data=newdf)
  mil_ind <- which(names(p$coefficients)=="MIL")
  
  recordMatrix$yrR[i] <- y
  recordMatrix$coefficientR[i] <- p$coefficients[mil_ind]
  recordMatrix$pcseR[i] <- p$pcses[mil_ind]
  recordMatrix$tvalueR[i] <- p$tvals[mil_ind]
  recordMatrix$pvalueR[i] <- p$pvals[mil_ind]
  
  i <- i+1
}

recordMatrix$nyears <- 2014 - recordMatrix$yr2
recordMatrix$pvCode1 <- "Not Sig."
recordMatrix$pvCode1[which(recordMatrix$pvalue < .1)] <- "p<.1"
recordMatrix$pvCode1[which(recordMatrix$pvalue < .05)] <- "p<.05"
recordMatrix$pvCode1[which(recordMatrix$pvalue < .01)] <- "p<.01"
recordMatrix$pvCode1[which(recordMatrix$pvalue < .001)] <- "p<.001"

recordMatrix$pvCode1 <- factor(recordMatrix$pvCode1, 
                               levels=c("p<.001","p<.01", "p<.05", "p<.1" , "Not Sig." ))


recordMatrix$pvCodeR <- "Not Sig."
recordMatrix$pvCodeR[which(recordMatrix$pvalueR < .1)] <- "p<.1"
recordMatrix$pvCodeR[which(recordMatrix$pvalueR < .05)] <- "p<.05"
recordMatrix$pvCodeR[which(recordMatrix$pvalueR < .01)] <- "p<.01"
recordMatrix$pvCodeR[which(recordMatrix$pvalueR < .001)] <- "p<.001"

recordMatrix$pvCodeR <- factor(recordMatrix$pvCodeR, 
                               levels=c("p<.001","p<.01", "p<.05", "p<.1" , "Not Sig." ))

recordMatrix$yr2Label <- paste(recordMatrix$yr2, "-2014", sep="")
recordMatrix$yr2Label[which(recordMatrix$nyears %% 5 != 0)] <- ""

recordMatrix$yrLabelR <- paste("1960-",  recordMatrix$yrR, sep="")
recordMatrix$yrLabelR[which(recordMatrix$nyears %% 5 != 0)] <- ""

#write.csv(recordMatrix, file="recordMatrix.csv")
library(ggExtra)
library(ggrepel)
library(RColorBrewer)

myPalette <- colorRampPalette(rev(brewer.pal(5, "Spectral")))
t <-  theme(text=element_text(size=14),
            plot.title = element_text(size=18), axis.title.x= element_text(size=16))

arrow <- arrow(length = unit(0.4, "cm"), type = "closed")

ggplot(data=recordMatrix, aes(y=coefficient, x=nyears)) +
  # geom_segment(aes(y=coefficient-pcse/2, yend=coefficient+pcse/2, xend=nyears),
  #              size=.2, alpha=0.5) +
  geom_point(aes(shape=pvCode1, color=pvCode1), size=3) +
  geom_line(linetype=3) +
  geom_hline(yintercept=0, linetype=2) +
  geom_text(aes(label=yr2Label), nudge_y = .002, size=3) +
  ggthemes::theme_gdocs() +
  # geom_segment(aes(y=coefficientR-pcseR/2, yend=coefficientR+pcseR/2, xend=nyears),
  #              size=.2, alpha=0.5) +
  geom_point(aes(y=coefficientR, shape=pvCodeR, color=pvCodeR), size=3) +
  geom_line(aes(y=coefficientR),  linetype=4) +
  geom_text(aes(y=coefficientR, label=yrLabelR), nudge_y = -.002, size=3) +
  scale_x_continuous( breaks = c(10, 15, 20, 25, 30, 35, 40, 45, 50)) + t + 
  scale_colour_brewer(palette = "Set2") +
  ggtitle(label="Coefficients of Military Expenditures from Regressions Including Different Years", 
          subtitle="(Two-Way Fixed Effects Models with Panel-Corrected Standard Errors and AR(1) Correction)") +
  xlab("Number of Years Included in the Regression") +
  ylab("Coefficient of Military Expenditures as % of GDP") +
  labs(color="p-values", shape="p-values") +
  theme(
    axis.line = element_line(arrow = arrow)
  ) +
  theme(legend.position = c(.9,.875))

###Coefficients Plot of Interaction Terms Models 3 & 7
library(GGally)
dfResults <- data.frame(term=names(p3$coefficients), 
                        estimate=p3$coefficients,
                        conf.low=p3$coefficients-p3$pcses/2,
                        conf.high=p3$coefficients+p3$pcses/2)
dfResults <- dfResults[219:272,]
dfResults$term <- as.character(1961:2014)
dfResults$variable <- "Model 3 (1960-2014)"
##model 7
dfResults7 <- data.frame(term=names(p7$coefficients), 
                        estimate=p7$coefficients,
                        conf.low=p7$coefficients-p7$pcses/2,
                        conf.high=p7$coefficients+p7$pcses/2)
dfResults7 <- dfResults7[199:232,]
dfResults7$term <- as.character(1981:2014)
dfResults7$variable <- "Model 7 (1980-2014)"

dfResults <- dplyr::bind_rows(dfResults, dfResults7)

library(ggthemes)
t <-  theme(text=element_text(size=12),
            plot.title = element_text(size=18), axis.title.x= element_text(size=16))
ggcoef(dfResults, conf.int = FALSE,
       mapping = aes(x = estimate, y = term, color = variable, shape=variable),
       size=3) +
  theme_gdocs() +
  scale_x_continuous( breaks = c(0, .025, .05, .075, .1, .125)) + 
  scale_colour_brewer(palette = "Set1") + t +
  theme(legend.position = c(.9,.15),
        axis.text.y = element_text(size=7, angle=0)) +
  ggtitle(label="Coefficients of Interaction between Military Expenditures and Time", 
          subtitle="(Coefficients and panel-corrected standard errors from models 3 and 7)") +
  xlab("Coefficient of Time Interaction") +
  labs(color="", shape="") +
  ylab(" Year of Interaction") +
  geom_errorbarh(data=dfResults, aes(x=estimate, 
                                     xmin=conf.low,
                                     xmax=conf.high,
                                     y=term,
                                     color=variable))
