library(readxl)
Panel_nas<- library(readxl)
Panel_nas<- read_excel("DataBase.xlsx")  # recordar renombrar base de datos

options("scipen"=100, "digits"=3)
# cargando paquetes para trabajar con data de panel
library(plm)
library(foreign)
library(gplots)
library(stargazer)
library(plm)

#ESTIMACION DE ESPECIFICACIONES

#Especificaciones estimadas a partir del estimados MCO

reg.P= lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt), data = Panel_nas)
reg.P3 = lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt), data = Panel_nas)
reg.P4 = lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt), data = Panel_nas)
reg.P5 = lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit), data = Panel_nas)
reg.P6 = lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit), data = Panel_nas)
reg.P7= lm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+(RCOLONij)+(PLINGij)+BGEOGRAFij, data = Panel_nas)

summary(reg.P)
summary(reg.P3)
summary(reg.P4)
summary(reg.P5)
summary(reg.P7)

stargazer(reg.P,reg.P3,reg.P4, reg.P5, reg.P6,reg.P7, digits = 3, dep.var.labels = c("Ln (Numero de emigrantes i en j/Población total i)*100,000"), covariate.labels=c("Razón PIB per cápita entre i y j","Distancia Física entre i y j","Tasa de desempleo i","Tasa de desempleo j","Esperanza de vida i","Esperanza de vida j","Años de escolarización prom. i","Años de escolarización prom. j","Índice de Gini i","Indice de Gini j","Tasa de homicidio i","Tasa de homicidio j","Conflicto político interno j (dummy)","Visado para entrada en j (dummy)","Inversión extranjera directa i","Relación colonial entre i y j(dummy)","Proximidad lingüística entre i y j(dummy)","Bordes geograficos compartidos entre i y j"), column.labels = c("Esp. 1","Esp.2","Esp.3","Esp.4","Esp.5","Esp.6","Esp.7"), single.row=FALSE, title = "Tabla 1: Estimación por Minímos Cuadrados Ordinarios(MCO)",type = "text", style = "aer",out = "tableMCO.htm")

#Especificaciones estimadas a partir del modelo MCO (agrupado)

reg.P1= plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt), index=c("id", "year"), model = "pooling", data = Panel_nas)
reg.P31 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt), index=c("id", "year"), model = "pooling",data = Panel_nas)
reg.P41 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt), index=c("id", "year"), model = "pooling", data = Panel_nas)
reg.P51 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit),index=c("id", "year"), model = "pooling", data = Panel_nas)
reg.P61 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit), index=c("id", "year"), model = "pooling",data = Panel_nas)
reg.P71= plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+(RCOLONij)+(PLINGij)+BGEOGRAFij,index=c("id", "year"), model = "pooling", data = Panel_nas)

summary(reg.P1)
summary(reg.P31)
summary(reg.P41)
summary(reg.P51)
summary(reg.P61)
summary(reg.P71)

stargazer(reg.P1,reg.P31,reg.P41, reg.P51, reg.P61,reg.P71, digits = 3, dep.var.labels = c("Ln (Numero de emigrantes i en j/Población total i)*100,000"), covariate.labels=c("Razón PIB per cápita entre i y j","Distancia Física entre i y j","Tasa de desempleo i","Tasa de desempleo j","Esperanza de vida i","Esperanza de vida j","Años de escolarización prom. i","Años de escolarización prom. j","Índice de Gini i","Indice de Gini j","Tasa de homicidio i","Tasa de homicidio j","Conflicto político interno j (dummy)","Visado para entrada en j (dummy)","Inversión extranjera directa i","Relación colonial entre i y j(dummy)","Proximidad lingüística entre i y j(dummy)","Bordes geograficos compartidos entre i y j"), column.labels = c("Esp. 1","Esp.2","Esp.3","Esp.4","Esp.5","Esp.6","Esp.7"), single.row=FALSE, title = "Tabla 1: Estimación por MCO Agrupado",type = "text", style = "aer",out = "tableMCO2.htm")


#Especificaciones estimadas a partir del estimador de efectos fijos
library(stargazer)
reg.W = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt), index=c("id", "year"), model = "within", data = Panel_nas)
reg.W3 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt), index=c("id", "year"), model="within", data = Panel_nas)
reg.W4 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt), index=c("id", "year"), model="within", data = Panel_nas)
reg.W5 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit), index=c("id", "year"), model="within", data = Panel_nas)
reg.W6 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit),index=c("id", "year"), model="within", data = Panel_nas)
reg.W7 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+(RCOLONij)+(PLINGij)+BGEOGRAFij,index=c("id", "year"), model="within", data = Panel_nas)

stargazer(reg.W,reg.W3,reg.W4, reg.W5, reg.W6,reg.W7, digits = 3, dep.var.labels = c("Ln (Numero de emigrantes i en j/Población total i)*100,000"), covariate.labels=c("Razón PIB per cápita entre i y j","Distancia Física entre i y j","Tasa de desempleo i","Tasa de desempleo j","Esperanza de vida i","Esperanza de vida j","Años de escolarización prom. i","Años de escolarización prom. j","Índice de Gini i","Indice de Gini j","Tasa de homicidio i","Tasa de homicidio j","Conflicto político interno j (dummy)","Visado para entrada en j (dummy)","Inversión extranjera directa i","Relación colonial entre i y j(dummy)","Proximidad lingüística entre i y j(dummy)","Bordes Geograficos compartidos entre i y j"), column.labels = c("Esp. 1","Esp.2","Esp.3","Esp.4","Esp.5","Esp.6","Esp.7"), single.row=FALSE, title = "Tabla 3: Estimación por Efectos Fijos",type = "text", style = "aer",out = "tablewithin.htm")

summary(reg.W)
summary(reg.W3)
summary(reg.W4)
summary(reg.W5)
summary(reg.W6)
summary(reg.W7)


#Especificaciones estimadas a partir del estimados de efectos aleatorios
library(plm)

reg.A= plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt), index=c("id", "year"), model = "random", data = Panel_nas)
reg.A3 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt), index=c("id", "year"), model="random", data = Panel_nas)
reg.A4 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt), index=c("id", "year"), model="random", data = Panel_nas)
reg.A5 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit), index=c("id", "year"), model="random", data = Panel_nas)
reg.A6 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit),index=c("id", "year"), model="random", data = Panel_nas)
reg.A7 = plm(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+(RCOLONij)+(PLINGij)+BGEOGRAFij,index=c("id", "year"), model="random", data = Panel_nas)

summary(reg.A)
summary(reg.A3)
summary(reg.A4)
summary(reg.A5)
summary(reg.A6)
summary(reg.A7)

stargazer(reg.A,reg.A3,reg.A4, reg.A5, reg.A6,reg.A7, digits = 3, dep.var.labels = c("Ln (Numero de emigrantes i en j/Población total i)*100,000"), covariate.labels=c("Razón PIB per cápita entre i y j","Distancia Física entre i y j","Tasa de desempleo i","Tasa de desempleo j","Esperanza de vida i","Esperanza de vida j","Años de escolarización prom. i","Años de escolarización prom. j","Índice de Gini i","Indice de Gini j","Tasa de homicidio i","Tasa de homicidio j","Conflicto político interno j (dummy)","Visado para entrada en j (dummy)","Inversión extranjera directa i","Relación colonial entre i y j(dummy)","Proximidad lingüística entre i y j(dummy)","Bordes geograficos compartidos entre i y j"), column.labels = c("Esp. 1","Esp.2","Esp.3","Esp.4","Esp.5","Esp.6","Esp.7"), single.row=FALSE, title = "Tabla 4: Estimación por Efectos Aleatorios",type = "text", style = "aer",out = "tablerandom.htm")

# SELECCIÓN DEL MODELO

# H0: modelo agrupado (MCO) vs H1: efectos fijos
pFtest(reg.W, reg.P) # preferible efectos fijos, se rechaza H0
pFtest(reg.W3, reg.P3) # preferible efectos fijos, se rechaza H0
pFtest(reg.W4, reg.P4) # preferible efectos fijos, se rechaza H0
pFtest(reg.W5, reg.P5) # preferible efectos fijos, se rechaza H0
pFtest(reg.W6, reg.P6) # preferible efectos fijos, se rechaza H0
pFtest(reg.W7, reg.P7) # preferible efectos fijos, se rechaza H0

#H0: modelo agrupado (MCO) vs H1: efectos aleatorios
plmtest(reg.P1, type=c("bp")) # preferible efectos aleatorios, se rechaza H0
plmtest(reg.P31, type=c("bp")) # preferible efectos aleatorios, se rechaza H0
plmtest(reg.P41, type=c("bp")) # preferible efectos aleatorios, se rechaza H0
plmtest(reg.P51, type=c("bp")) # preferible efectos aleatorios, se rechaza H0
plmtest(reg.P61, type=c("bp")) # preferible efectos aleatorios, se rechaza H0
plmtest(reg.P71, type=c("bp")) # preferible efectos aleatorios, se rechaza H0

#Ante la sospecha de variaciones en el tiempo y efectos individuales entre las variables que conforman la muestra y los resultados de este test; se debe proceder a las tecnicas propias utilizadas en datos de panel y descartar los resultados del estimados por efectos agrupados o pool de datos.

#H0: efectos aleatorios vs H1: efectos fijos
phtest(reg.W, reg.A) #preferible efectos fijos, se rechaza H0
phtest(reg.W3, reg.A3) #preferible efectos fijos, se rechaza H0
phtest(reg.W4, reg.A4) #preferible efectos fijos, se rechaza H0
phtest(reg.W5, reg.A5) #preferible efectos aleatorios, no se rechaza H0
phtest(reg.W6, reg.A6) #preferible efectos aleatorios, no se rechaza H0
phtest(reg.W7, reg.A7) #preferible efectos aleatorios, no se rechaza H0

#Otras pruebas al modelo seleccionado

#test de depedencia trasnversal/ correlacion contemporanea
pcdtest(reg.A, test = c("lm")) 
pcdtest(reg.A3, test = c("lm"))
pcdtest(reg.A4, test = c("lm"))
pcdtest(reg.A5, test = c("lm"))
pcdtest(reg.A6, test = c("lm"))
pcdtest(reg.A7, test = c("lm"))

pctest(reg.A, test = c("cd")) 
pcdtest(reg.A3, test = c("cd"))
pcdtest(reg.A4, test = c("cd"))
pcdtest(reg.A5, test = c("cd"))
pcdtest(reg.A6, test = c("cd"))
pcdtest(reg.A7, test = c("cd"))

#Tras los resultados del test CD (cross-sectional dependence) se detecta correlacion seccional entre las variables, al igual que al aplicar B-P/LM sí. No obtante, dada la estructura tipo micro-panel (pocos periodos (6) y una buena cantidad de observaciones) de la data, se entiende que este no debe ser considerado un problema.

#Test de correlacion serial
pbgtest(reg.A)
pbgtest(reg.A3)
pbgtest(reg.A4)
pbgtest(reg.A5)
pbgtest(reg.A6)
pbgtest(reg.A7)

#Hay evidencia de correlacion serial en los tests de Breusch-Godfrey/Wooldridge, no obtante dado que estamos trabajando con una base de datos que tiene caracteristicas de micro panel, tampoco debe ser considerada como un problema a los fines.

# Test de la raiz unitaria o estacionariedad
Panel.set <- pdata.frame(Panel_nas, index = c("id", "year"))
library(tseries)
adf.test(Panel.set$EMIGit_n, k=2) 

# El test de Dickey-Fuller indica que no se veriifca la existencia de una raiz unitaria en la variable dependiente, siendo esta estacionaria.

library(lmtest)
#Test de heterocedasticidad
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+factor(id), data = Panel_nas, studentize=F) 
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+factor(id), data = Panel_nas, studentize=F) 
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+factor(id), data = Panel_nas, studentize=F) 
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+factor(id), data = Panel_nas, studentize=F) 
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+factor(id), data = Panel_nas, studentize=F)
bptest(log(EMIGit_n)~(PIBij_n)+log(DISTFij)+(DESEMPit)+(DESEMPjt)+(ESPVIDAit)+(ESPVIDAjt)+(ESCit)+(ESCjt)+(GINIit)+(GINIjt)+(HOMit)+(HOMjt)+(CONFLICTjt)+(PMIit)+(IEDit)+(RCOLONij)+(PLINGij)+BGEOGRAFij+factor(id), data = Panel_nas, studentize=F) 

#Se evidencia hetrocedasticidad tanto en la espcificacion con el estimador de efectos fijos y efectos aleatorios. A los fines de correcion del estimador con efectos fijos se recomienda la correcccion por el metodo de "arellano" que da estimadores robustos corregidos por heterocedasticidad y correlacion serial. En lo que respecta al modelo de efecttos aleatorios se recomienda la aplicacion la estimacion de matriz de covarianza consistente de heterocedasticidad y autocorrelación (HAC),


#Coeficientes consistentes ante heterocedasticidad y correlacion serial en el estimador de efectos aleatorios (vcovHC) por metodo de arellano (1989)

library(sandwich)
a<-coeftest(reg.A, vcovHC(reg.A, method = "arellano")) 
c<-coeftest(reg.A3, vcovHC(reg.A3, method = "arellano")) 
d<-coeftest(reg.A4, vcovHC(reg.A4, method = "arellano"))  
e<-coeftest(reg.A5, vcovHC(reg.A5, method = "arellano")) 
f<-coeftest(reg.A6, vcovHC(reg.A6, method = "arellano")) 
h<-coeftest(reg.A7, vcovHC(reg.A7, method = "arellano")) 

stargazer(a,c,d,e,f,h, digits = 3, dep.var.labels = c("Ln (Numero de emigrantes i en j/Población total i)*100,000"), covariate.labels=c("Razón PIB per cápita entre i y j","Distancia Física entre i y j","Tasa de desempleo i","Tasa de desempleo j","Esperanza de vida i","Esperanza de vida j","Años de escolarización prom. i","Años de escolarización prom. j","Índice de Gini i","Indice de Gini j","Tasa de homicidio i","Tasa de homicidio j","Conflicto político interno j (dummy)","Visado para entrada en j (dummy)","Inversión extranjera directa i","Relación colonial entre i y j(dummy)","Proximidad lingüística entre i y j(dummy)","Bordes geograficos compartidos entre i y j"), column.labels = c("Esp. 1","Esp.2","Esp.3","Esp.4","Esp.5","Esp.6","Esp.7"), single.row=FALSE, title = "Tabla 3: Estimación por Efectos Aleatorios",type = "text", style = "aer", df=FALSE ,out = "tablerandomajustado.htm")