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

source("loadData.R")

##List of countries
library(countrycode)
##UnBalanced Panels
mf_ub$country <- countrycode(mf_ub$ID, "wb", "country.name")
mf_ub$income <- NA
mf_ub$income[which(mf_ub$ID %in% HIC_i)] <- "High and Middle Income"
mf_ub$income[which(mf_ub$ID %in% MIC_i)] <- "High and Middle Income"
mf_ub$income[which(mf_ub$ID %in% LICs)] <- "Low Income"
mf_ub <- mf_ub %>% arrange(income, country)
mf_ub <- mf_ub[,c("country", "n")]
#mf_ub <- mf_ub %>% tidyr::unite(Observations, country, n, sep="")
mf_ub <- mf_ub %>% tidyr::unite(Observations, country, n, sep=": ")
stargazer(mf_ub,  type="html", summary = FALSE, out="results/FINAL_countryList_ub.html",
          title="Countries Included in the Unbalanced Panel Data Set")


###Generate summary stats.  
df_describe <- psych::describe(df_ub, skew=T)  %>% as.data.frame()
df_describe <- df_describe[-c(1:2,7),-c(1,6,7,10,13)]
df_describe <- as.data.frame(apply(df_describe, 2, round, 3))

varLabels <- c( "Territorial emissions in tCO2 per person",
                "GDP per capita", 
              "Urban population as % of total population", 
              "Military expenditures as % GDP")


df_describe$Variables <- varLabels
df_describe <- df_describe[,c(9,1:8)]
row.names(df_describe) <- NULL

stargazer(df_describe,  type="html", summary = FALSE, out="results/summaryStats.html",
          title="Descriptive Statistics of Included Variables")

##EXPORTING TO .DOC
library(rtf)
Diag <- RTF(file="results/Data_Summary.doc")
addParagraph(Diag, "\n")
addTable(Diag, df_describe, col.justify = 'L')

corrP2 <- cor(df_ub[,c(3:6)])
corrP2 <- apply(corrP2, 2, round, 3)
row.names(corrP2) <- varLabels
colnames(corrP2) <- c("1.", "2.", "3.", "4.")
stargazer(corrP2,  type="html", summary = FALSE, out="results/Correlations_unbalanced.html",
          title="Bivariate Correlations for Unbalanced Panel Dataset")

addParagraph(Diag, "\n")
addTable(Diag, corrP2, col.justify = 'L')

#################SET MODEL FORMULAs#####################################

##  - 
df$ID <- factor(df$ID)
df_ub$ID <- factor(df_ub$ID)

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

##PW REGRESSIONS WITH AR(1) CORRECTION AND PANEL-CORRECTED STANDARD ERRORS
p1 <- pw.function(f1, data=df_ub)
p2 <- pw.function(f2, data=df_ub)
p3 <- pw.function(f3, data=df_ub)
p4 <- pw.function(f4, data=df_ub)

df_ub80 <- subset(df_ub, YR >= 1980)

p5 <- pw.function(f1, data=df_ub80)
p6 <- pw.function(f2, data=df_ub80)
p7 <- pw.function(f3, data=df_ub80)
p8 <- pw.function(f4, data=df_ub80)

###EXPORTING RESULTS  
title1 <- "Coefficients for the Regression of Per Capita Carbon Dioxide Emissions, 1960-2014: Two-Way Fixed Effects Regression with Panel-Corrected Standard Errors and Prais-Winston AR(1) Correction."

models <- list(p1, p2, p3, p4, p5, p6, p7, p8)

xVarLabels <- c(
  "GDP per capita", 
  "Urban population as % of total population", 
  "Military expenditures as % GDP", 
  "Military expenditures as % GDP x Low Income", 
  "Year",
  "Military expenditures as % GDP x Year"

)
Results <- outputReg_pw(regModels=models[1:8], addVIFs = FALSE)
nC <- ncol(Results)-1
names(Results)[2:(nC+1)] <- as.character(1:nC)
Results <- Results[-1,]
vn <- length(xVarLabels)
#Results$vars[1:vn] <- xVarLabels


##EXPORTING TO .DOC
addParagraph(Diag, "\n")
addTable(Diag, Results, col.justify = 'L')
addParagraph(Diag, "\n")
done(Diag) 

stargazer(Results,  type="html", summary = FALSE, out="results/FINAL_results.html",
          title=title1)
