#################################################################################
# Regression Analysis                                                           #
#################################################################################

# Figure 10 Main Paper
# Figure 11 Main Paper

#################################################################################
# Head                                                                          #
#################################################################################

# remove the old Stuff
rm(list=ls())

# load packages
library(spdep)
library(pspline)
library(stargazer)
library(ggplot2)
library(car)

# Set working Directory
setwd("")

# load the data
load("Data.RData")

# Open containers for residuals and run the regression
res_series_u<-c()
res_series_epsilon<-c()
res_series_test<-list()
res_errorsar_TC3<-list()
for (t in 1952:2016){
  print(t)
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  X_use<-cbind(1,
               X[which(X[,2]==t),which(colnames(X)=="lGDP_i")],
               X[which(X[,2]==t),which(colnames(X)=="lGDP_j")],
               X[which(X[,2]==t),which(colnames(X)=="DA_ij")],
               X[which(X[,2]==t),which(colnames(X)=="POL_ij")],
               X[which(X[,2]==t),which(colnames(X)=="D_ij")])
  
  W_all<-TC3[[t-1951]]
  diag(W_all)<-0
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  res_errorsar_TC3[[t-1951]]<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  
  insert<-(diag(length(Y_use))-res_errorsar_TC3[[t-1951]]$lambda*W_all)%*%(Y_use-X_use%*%res_errorsar_TC3[[t-1951]]$coefficients)
  res_series_epsilon<-c(res_series_epsilon,insert) 
  res_series_u<-c(res_series_u,Y_use-X_use%*%res_errorsar_TC3[[t-1951]]$coefficients)
  res_series_test[[t-1951]]<-shapiro.test(insert)
}

#################################################################################
# Figure 10                                                                     #
#################################################################################

# Figure 10: QQ-Plot (left) and Histogram (right) of standardized Residuals.

pdf(paste("2) Estimating the Models/figure10.pdf",sep=""),width = 10,height = 5)  
par(mfrow=c(1,2))
lmres_series_epsilon<-res_series_epsilon/sd(res_series_epsilon)
qqPlot(lmres_series_epsilon, main="QQ Plot",ylab="Standardized Residuals")

hist(lmres_series_epsilon, freq=FALSE, 
     main="Distribution of standardized Residuals",xlim=c(-5,5),xlab="")
sresid<-lmres_series_epsilon
xfit<-seq(-5,5,length=1000) 
yfit<-dnorm(xfit) 
lines(xfit, yfit)
dev.off()


Y<-X[,3]
Tradecorr<-res_series_u-res_series_epsilon
X_beta<-Y-res_series_u

#################################################################################
# Figure 11                                                                     #
#################################################################################

# Figure 11: Kernel Density Estimates of the Tradecorrelation Residuals
# for selected Countries in the Time Periods 1952-1991 (left) and 1992-2016
# (right).
pdf("2) Estimating the Models/figure11.pdf",width=10,height = 5)
par(mfrow=c(1,2))
t=c(1952:1991)
raneff_spat<-cbind(Tradecorr[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
exports<-cbind(Y[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
fitted_values<-cbind(X_beta[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
spat<-cbind(raneff_spat,R[which(X[,2]%in%t),2])
residuals_here<-cbind(res_series_epsilon[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
# Order them by nation
spat<-spat[order(spat[,2]),]
us<-spat[which(spat[,2]==202),1]
us_res<-residuals_here[which(fitted_values[,2]==202),1]
plot(density(us),col="black",main="Tradecorrelation 1952-1991",xlim=c(-1,2),ylim=c(0,4.6),lwd=3,xlab="")
su<-spat[which(spat[,2]==176),1]
su_res<-residuals_here[which(fitted_values[,2]==176),1]
lines(density(su),col="gray",lwd=3)
uk_res<-residuals_here[which(fitted_values[,2]==201),1]
uk<-spat[which(spat[,2]==201),1]
lines(density(uk),lty=2)
de<-spat[which(spat[,2]==72),1]
de_res<-residuals_here[which(fitted_values[,2]==72),1]
lines(density(de),lty=3)
fr<-spat[which(spat[,2]==67),1]
fr_res<-residuals_here[which(fitted_values[,2]==67),1]
lines(density(fr),lty=4)
chn<-spat[which(spat[,2]==39),1]
chn_res<-residuals_here[which(fitted_values[,2]==39),1]
lines(density(chn),lty=5)
can<-spat[which(spat[,2]==34),1]
can_res<-residuals_here[which(fitted_values[,2]==34),1]
lines(density(can),lty=6)
legend("topright",c("USA","SUN","GBR","DEU","FRA","CHN","CAN"),col=c("black","gray","black","black","black","black","black"),lty=c(1,1,2,3,4,5,6),lwd=c(3,3,1,1,1,1,1),bty = "n")

t=c(1992:2016)
raneff_spat<-cbind(Tradecorr[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
exports<-cbind(Y[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
fitted_values<-cbind(X_beta[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
spat<-cbind(raneff_spat,R[which(X[,2]%in%t),2])
residuals_here<-cbind(res_series_epsilon[which(X[,2]%in%t)],R[which(X[,2]%in%t),2])
# Order them by nation
spat<-spat[order(spat[,2]),]
us<-spat[which(spat[,2]==202),1]
us_res<-residuals_here[which(fitted_values[,2]==202),1]
plot(density(us),col="black",main="Tradecorrelation 1992-2016",xlim=c(-1,2),ylim=c(0,4.6),lwd=3,xlab="")
su<-spat[which(spat[,2]==154),1]
su_res<-residuals_here[which(fitted_values[,2]==154),1]
lines(density(su),col="gray",lwd=2)
uk_res<-residuals_here[which(fitted_values[,2]==201),1]
uk<-spat[which(spat[,2]==201),1]
lines(density(uk),lty=2)
de<-spat[which(spat[,2]==72),1]
de_res<-residuals_here[which(fitted_values[,2]==72),1]
lines(density(de),lty=3)
fr<-spat[which(spat[,2]==67),1]
fr_res<-residuals_here[which(fitted_values[,2]==67),1]
lines(density(fr),lty=4)
chn<-spat[which(spat[,2]==39),1]
chn_res<-residuals_here[which(fitted_values[,2]==39),1]
lines(density(chn),lty=5)
can<-spat[which(spat[,2]==34),1]
can_res<-residuals_here[which(fitted_values[,2]==34),1]
lines(density(can),lty=6)
legend("topright",c("USA","RUS","GBR","DEU","FRA","CHN","CAN"),col=c("black","gray","black","black","black","black","black"),lty=c(1,1,2,3,4,5,6),lwd=c(3,3,1,1,1,1,1),bty = "n")

dev.off()

rm(list=ls())






