PDeseas <- function(series,Data){ Country=Data$Country Date=Data$Date Countries=unique(Country) PDeseas=rep(NA,length(series)) for (i in 1:length(Countries)) { sa_series=series[which(Country==Countries[i])] if (sum(!is.na(sa_series))>0) { minyear=as.numeric(format(min(Date[which(Country==Countries[i])]),"%Y")) minquarter=as.numeric(substr(quarters(min(Date[which(Country==Countries[i])])), 2, 2)) maxyear=as.numeric(format(max(Date[which(Country==Countries[i])]),"%Y")) maxquarter=as.numeric(substr(quarters(max(Date[which(Country==Countries[i])])), 2, 2)) sa_series=ts(sa_series,c(minyear,minquarter),c(maxyear,maxquarter),4) PDeseas[which(Country==Countries[i])]=as.numeric(sa_series) try(sa_series <- seas(sa_series,x11 = "",na.action=na.exclude),silent=TRUE) if (class(sa_series)=="seas") PDeseas[which(Country==Countries[i])]=as.numeric(final(sa_series)) } } return(PDeseas) } PDiff <-function(vecIn,vecID){ IDs=unique(vecID) vecOut=vecIn vecOut[]=NA for (i in IDs) { Ind=which(vecID==i) vecOut[Ind]=c(NA,vecIn[Ind[2]:Ind[length(Ind)]]-vecIn[Ind[1]:Ind[length(Ind)-1]]) } PDiff<-vecOut } PDetrend <-function(vecIn,vecID){ IDs=unique(vecID) vecOut=vecIn vecOut[]=NA for (i in IDs) { Ind=which(vecID==i & ! is.na(vecIn)) if (length(Ind)>3) { fit = lm(vecIn[Ind]~seq(1,length(Ind))) vecOut[Ind]=vecIn[Ind]-fitted(fit) } } PDiff<-vecOut } PLag <-function(vecIn,vecID,Lag=1){ IDs=unique(vecID) vecOut=vecIn vecOut[]=NA for (i in IDs) { Ind=which(vecID==i) vecOut[Ind]=c(rep(NA,Lag),vecIn[Ind[1]:Ind[length(Ind)-Lag]]) } PDiff<-vecOut } CollectResults<- function(ListResults){ LR=character() SR=character() Rest=character() for (i in 1:length(ListResults)) { RowNames=rownames(ListResults[[i]]$output) LR=c(LR,RowNames[1:(which(RowNames=="Adjustment")-1)]) SR=c(SR,RowNames[(which(RowNames=="Adjustment")):which(RowNames=="Const")]) Rest=c(Rest,RowNames[(which(RowNames=="Const")+1):length(RowNames)]) } RowNames=unique(c(LR,SR,Rest)) CollectResults=data.frame(matrix(NA,length(RowNames),length(ListResults))) rownames(CollectResults)=RowNames for (i in 1:length(ListResults)) { CollectResults[,i]=ListResults[[i]]$output[RowNames,1] CollectResults[which(!rownames(CollectResults) %in% rownames(ListResults[[i]]$output)),i]="" } Order=1:length(RowNames) Order[which(RowNames=="Trend")]=which(RowNames=="Const")+0.1 CollectResults[order(Order),] } CollectResultsMG<- function(ListResults){ LR=character() SR=character() Rest=character() for (i in 1:length(ListResults)) { RowNames=rownames(ListResults[[i]]$output) LR=c(LR,RowNames[1:(which(RowNames=="Adjustment")-1)]) SR=c(SR,RowNames[(which(RowNames=="Adjustment")):which(RowNames=="Const")]) Rest=c(Rest,RowNames[(which(RowNames=="Const")+1):length(RowNames)]) } RowNames=unique(c(LR,SR,Rest)) CollectResultsMG=data.frame(matrix(NA,length(RowNames),length(ListResults))) rownames(CollectResultsMG)=RowNames for (i in 1:length(ListResults)) { CollectResultsMG[,i]=ListResults[[i]]$output[RowNames,2] CollectResultsMG[which(!rownames(CollectResultsMG) %in% rownames(ListResults[[i]]$output)),i]="" } Order=1:length(RowNames) Order[which(RowNames=="Trend")]=which(RowNames=="Const")+0.1 CollectResultsMG[order(Order),] } CollectSRCoeff <- function(ListResults){ Names=rownames(ListResults[[1]]) for (i in 2:length(ListResults)) { Names=c(Names,rownames(ListResults[[i]])) } Names=unique(Names) CollectResults=data.frame(matrix(NA,length(Names),length(ListResults))) rownames(CollectResults)=Names for (i in 1:length(ListResults)) { a=ListResults[[i]] a["Trend","PMG"]=-a["Trend","PMG"]/a["Adjustment","PMG"] b=a[(which(rownames(a)=="Adjustment_pval")+1):(which(rownames(a)=="Hausman p-val")-1),-2] names(b)=rownames(a)[(which(rownames(a)=="Adjustment_pval")+1):(which(rownames(a)=="Hausman p-val")-1)] b=round(b,3) for (j in seq(1,length(b),2)) { if (b[j+1]<.01) { b[j]=paste(b[j],"***",sep="") } else if (b[j+1]<.05) {b[j]=paste(b[j],"**",sep="") } else if (b[j+1]<.1) {b[j]=paste(b[j],"*",sep="") } } b=b[-grep("pval",names(b))] CollectResults[match(names(b),Names),i]=b } CollectSRCoeff=CollectResults[rowSums(is.na(CollectResults)) != ncol(CollectResults),] } OptimBIC <-function(Data,Names,NamePDF=NULL,MinObs,MinDate,MaxDate,Trend){ Combs=expand.grid(rep(list(2:5),2)) BIC=rep(NA,3) Index=NULL for (i in nrow(Combs):1) { Reg=NULL try(Reg<-OptimPMG(Data,Names,Combs[i,1],Combs[i,2],NamePDF,MinObs = MinObs, MinDate=MinDate,Index=Index,Trend=Trend),silent=TRUE) if(is.null(Reg)){BIC[i]=Inf}else{BIC[i]=Reg$detailed_output["BIC","PMG"]} if (i==nrow(Combs)) {Index=Reg$Index} } OptimBIC=OptimPMG(Data,Names,Combs[which(BIC==min(BIC,na.rm=TRUE)),1],Combs[which(BIC==min(BIC,na.rm=TRUE)),2],NamePDF,MinObs = MinObs, MinDate=MinDate,MaxDate=MaxDate,Trend=Trend) } OptimPMG <- function(Data,Names,p,q,NamePDF=NULL,Trend=FALSE,MinObs=0,MinDate='1900-01-01',MaxDate='2100-01-01',Index=NULL){ Country=Data$Country Date=as.Date(Data$Date) dy=PDiff(Data[,Names[1]],Country) ly=PLag(Data[,Names[1]],Country) X=matrix(NA,nrow(Data),length(Names)-1) colnames(X)=Names[2:(length(Names))] W=matrix(NA,nrow(Data),(p-1)+(q)*ncol(X)) if (p>1) { for (i in 1:(p-1)) { W[,i]=PLag(dy,Country,i) } ColnamesW=paste("d",Names[1],1:(p-1),sep="_") }else{ColnamesW=character(0)} j=max(p-1,0) for (i in 1:ncol(X)) { X[,i]=Data[,colnames(X)[i]] for (k in 0:(q-1)) { j=j+1 W[,j]=PLag(PDiff(X[,i],Country),Country,k) } ColnamesW=c(ColnamesW,paste("d",Names[1+i],0:(q-1),sep="_")) } colnames(W)=gsub("_0","",ColnamesW) if (is.null(Index)){ Index=which(complete.cases(cbind(1:length(ly),ly,dy,X,W,Country)) & Data$Date>=as.Date(MinDate) & Data$Date<=as.Date(MaxDate)) # Exclude 2007-2009 #Index=which(complete.cases(cbind(1:length(ly),ly,dy,X,W,Country)) & Data$Date>=as.Date(MinDate) & Data$Date<=as.Date(MaxDate) & !(Data$Date>="2007-01-01" & Data$Date<="2009-12-30")) Index=Index[which(Index %in% which(Country %in% names(which(table(Country[Index])>MinObs))))] } ly=ly[Index];dy=dy[Index];X=X[Index,];W=W[Index,];Country=Country[Index];Date=Date[Index] W=cbind(W,rep(1,nrow(W)));colnames(W)[ncol(W)]="Const" if (Trend) {W=cbind(W,(1:nrow(W))/400);colnames(W)[ncol(W)]="Trend"} # Exclude 2007-2009 # if (Trend) {W[which(Date>="2010-01-01"),"Trend"]=W[which(Date>="2010-01-01"),"Trend"]+3} ID=unclass(factor(Country)) vars=as.data.frame(cbind(ly,X,dy,W[,paste("d_",Names[2:length(Names)],sep="")])) vars$resid=rep(NA,length(dy));vars$coint=rep(NA,length(dy));vars$fit=rep(NA,length(dy)) maxIter=1000 dLLMin=10^-10 N=length(unique(ID)) #--------------------------------- # MG #--------------------------------- PhiMG=0 coeffsMGsingle=matrix(NA,ncol(cbind(ly,X,W)),N) R2MG=0;AdjR2MG=0 for (i in 1:N) { obsi=which(ID==i) fit=lm(dy[obsi]~ly[obsi]+X[obsi,]+W[obsi,]-1) coeffsMGsingle[,i]=fit$coefficients coeffsMGsingle[2:(ncol(X)+1),i]=-coeffsMGsingle[2:(ncol(X)+1),i]/coeffsMGsingle[1,i] R2MG=R2MG+summary(fit)$r.squared AdjR2MG=AdjR2MG+summary(fit)$adj.r.squared } R2MG=R2MG/N AdjR2MG=AdjR2MG/N coeffsMG=rowMeans(coeffsMGsingle) CovarThetaMG=matrix(NA,ncol(X),ncol(X)) for (i in 1:ncol(X)) { for (j in 1:ncol(X)) { CovarThetaMG[i,j]=sum((coeffsMGsingle[i+1,]-coeffsMG[i+1])*(coeffsMGsingle[j+1,]-coeffsMG[j+1]))/(N*(N-1)) } } #--------------------------------- # PMG #--------------------------------- thetaPMG=rep(0,ncol(X)) Tvec=as.data.frame(table(ID)) Tvec=Tvec[,2] logL1=-Inf logL2=-Inf H=list() phi=rep(NA,N) sigma2=rep(NA,N) for (i in 1:N) { obsi=which(ID==i) H[[i]]=diag(Tvec[i])-W[obsi,]%*%solve(t(W[obsi,])%*%W[obsi,])%*%t(W[obsi,]) } # Optimisation a=0 dLL=Inf while (a < maxIter & dLL>=dLLMin) { a=a+1 A=matrix(0,length(thetaPMG),length(thetaPMG)) B=matrix(0,length(thetaPMG),1) Ksi=ly-X%*%thetaPMG for (i in 1:N) { obsi=which(ID==i) phi[i]=as.numeric(solve(t(Ksi[obsi])%*%H[[i]]%*%Ksi[obsi])%*%t(Ksi[obsi])%*%H[[i]]%*%dy[obsi]) sigma2[i]=as.numeric(t(dy[obsi]-phi[i]*Ksi[obsi])%*%H[[i]]%*%(dy[obsi]-phi[i]*Ksi[obsi]))/Tvec[i] A=A+phi[i]^2/sigma2[i]*t(X[obsi,])%*%H[[i]]%*%X[obsi,] B=B+phi[i]/sigma2[i]*t(X[obsi,])%*%H[[i]]%*%t(dy[obsi]-phi[i]%*%ly[obsi]) } thetaPMG=-solve(A)%*%B logL1=logL2 logL2=0 for (i in 1:N) { obsi=which(ID==i) logL2=logL2-Tvec[i]/2*log(2*pi*sigma2[i])-.5*1/sigma2[i]*t(dy[obsi]-phi[i]*Ksi[obsi])%*%H[[i]]%*%(dy[obsi]-phi[i]*Ksi[obsi]) } dLL=logL2-logL1 } CoeffsSRPMG=matrix(NA,ncol(W)+1,N) CovarEstimPMG=matrix(0,ncol(X)+N*(1+ncol(W)),ncol(X)+N*(1+ncol(W))) for (i in 1:N) { obsi=which(ID==i) fit=lm(dy[obsi]~Ksi[obsi]+W[obsi,]-1) CoeffsSRPMG[,i]=fit$coefficients vars$resid[obsi]=fit$residuals vars$fit[obsi]=fitted(fit) vars$coint[obsi]=X[obsi,]%*%thetaPMG-mean(X[obsi,]%*%thetaPMG)+mean(ly[obsi]) if (Trend) {vars$coint[obsi]=vars$coint[obsi]-fit$coefficients[1]*(W[obsi,"Trend"]-mean(W[obsi,"Trend"])*CoeffsSRPMG[nrow(CoeffsSRPMG),i])} CovarEstimPMG[1:ncol(X),1:ncol(X)]=CovarEstimPMG[1:ncol(X),1:ncol(X)]+phi[i]^2/sigma2[i]*t(X[obsi,])%*%X[obsi,] CovarEstimPMG[1:ncol(X),ncol(X)+i]=-phi[i]/sigma2[i]*t(X[obsi,])%*%Ksi[obsi,] CovarEstimPMG[1:ncol(X),(ncol(X)+N+1+(i-1)*ncol(W)):(ncol(X)+N+i*ncol(W))]=-phi[i]/sigma2[i]*(t(X[obsi,]) %*% W[obsi,]) CovarEstimPMG[ncol(X)+i,ncol(X)+i]=(t(Ksi[obsi]) %*% Ksi[obsi]) /sigma2[i] CovarEstimPMG[ncol(X)+i,(ncol(X)+N+1+(i-1)*ncol(W)):(ncol(X)+N+i*ncol(W))]=t(Ksi[obsi]) %*% W [obsi,]/sigma2[i] CovarEstimPMG[(ncol(X)+N+(i-1)*ncol(W)+1):(ncol(X)+N+i*ncol(W)),(ncol(X)+N+1+(i-1)*ncol(W)):(ncol(X)+N+i*ncol(W))]=t(W[obsi,])%*%W[obsi,]/sigma2[i] } CovarEstimPMG=as.matrix(Matrix::forceSymmetric(CovarEstimPMG,uplo="U")) CovarEstimPMG=solve(CovarEstimPMG) StdErr=sqrt(diag(CovarEstimPMG)) LR=data.frame(matrix(NA,length(thetaPMG),4),row.names=rownames(thetaPMG)) colnames(LR)=c("Coeff","StdErr","tStat","pval") LR[,1]=thetaPMG LR[,2]=StdErr[1:length(thetaPMG)] LR[,3]=LR[,1]/LR[,2] LR[,4]=as.matrix((1 - pnorm(abs(LR[,3]), mean = 0, sd = 1, lower.tail = TRUE)) * 2) CovarThetaPMG=CovarEstimPMG[1:length(thetaPMG),1:length(thetaPMG)] #--------------------------------- # Output table #--------------------------------- Output = data.frame(matrix(NA,ncol(X)+ncol(W)+1,2)) colnames(Output)=c("PMG","MG") rownames(Output)=c(colnames(X),"Adjustment",colnames(W)) Output$PMG[1:length(thetaPMG)]=thetaPMG Output$PMG[(length(thetaPMG)+1):nrow(Output)]=rowMeans(CoeffsSRPMG) Output$MG=c(coeffsMG[2:(ncol(X)+1)],coeffsMG[1],coeffsMG[(ncol(X)+2):length(coeffsMG)]) # Pvals Output <- Output[rep(1:nrow(Output), each = 2), ] Output[1:nrow(Output) %% 2 == 0, ] <- NA rownames(Output) = gsub("\\.1","_pval",rownames(Output)) Output$PMG[seq(2,length(thetaPMG)*2,2)]=as.matrix((1 - pnorm(abs(thetaPMG/sqrt(diag(CovarThetaPMG))), mean = 0, sd = 1, lower.tail = TRUE)) * 2) sd_err=sqrt(diag(CovarThetaPMG)) names(sd_err)=Names[2:length(Names)] Output$PMG[seq((ncol(X)*2+2),nrow(Output),2)]=(1 - pnorm(abs(rowMeans(CoeffsSRPMG)/sqrt(rowSums((CoeffsSRPMG-rowMeans(CoeffsSRPMG))^2)/(N*(N-1)))),0,1, lower.tail = TRUE)) * 2 pvalsMG=(1 - pnorm(abs(rowMeans(coeffsMGsingle)/sqrt(rowSums((coeffsMGsingle-rowMeans(coeffsMGsingle))^2)/(N*(N-1)))),0,1, lower.tail = TRUE)) * 2 Output$MG[seq(2,nrow(Output),2)]=c(pvalsMG[2:(ncol(X)+1)],pvalsMG[1],pvalsMG[(ncol(X)+2):length(pvalsMG)]) # Hausman test Hausman=t(coeffsMG[2:(ncol(X)+1)]-thetaPMG)%*%solve(CovarThetaMG-CovarThetaPMG)%*%(coeffsMG[2:(ncol(X)+1)]-thetaPMG) if (Hausman>0) { HausmanPval=pchisq(Hausman, df=length(thetaPMG), lower.tail=FALSE) }else{HausmanPval=0} rownames(Output[nrow(Output)+1,])="Hausman p-val" Output[nrow(Output),1]=HausmanPval;Output[nrow(Output),1]=HausmanPval # LogL etc rownames(Output[nrow(Output)+1,])="LogL" Output[nrow(Output),1]=logL2 rownames(Output[nrow(Output)+1,])="BIC" Output[nrow(Output),1]=(log(length(dy))*(length(thetaPMG)+1+ncol(W)*N)-2*logL2)/length(dy) R2=crossprod(vars$fit-mean(dy))/crossprod(dy-mean(dy)) rownames(Output[nrow(Output)+1,])="R2" Output[nrow(Output),1]=R2 Output[nrow(Output),2]=R2MG rownames(Output[nrow(Output)+1,])="Adj. R2" Output[nrow(Output),1]=1-(1-R2)*(length(dy)-1)/(length(dy)-(length(thetaPMG)+(1+ncol(W))*N)) Output[nrow(Output),2]=AdjR2MG vars=vars[,-ncol(vars)] rownames(Output[nrow(Output)+1,])="N. obs" Output[nrow(Output),1]=length(dy);Output[nrow(Output),2]=length(dy) rownames(Output[nrow(Output)+1,])="N. countries" Output[nrow(Output),1]=N;Output[nrow(Output),2]=N # Significance stars Output2=round(Output,3) for (i in nrow(Output):2) { if (grepl("pval",rownames(Output)[i])) { for (j in 1:2) { if (!is.na(Output2[i,j])) { if (Output2[i,j]<.01) { Output2[i-1,j]=paste(Output2[i-1,j],"***",sep="") } else if (Output2[i,j]<.05) {Output2[i-1,j]=paste(Output2[i-1,j],"**",sep="") } else if (Output2[i,j]<.1) {Output2[i-1,j]=paste(Output2[i-1,j],"*",sep="") } } } Output2=Output2[-i,] } } rownames(Output2[nrow(Output2)+1,])="p/q" Output2[nrow(Output2),1]=paste(p,q,sep="/") show(Output2) colnames(vars)[1]=paste("l_",Names[1],sep="");colnames(vars)[length(Names)+1]=paste("d_",Names[1],sep="") OutputList=list(Output2,vars,Date,Country,Output,Index,coeffsMGsingle,sd_err) names(OutputList)=c("output","vars","Date","Country","detailed_output","Index","coeffsMGsingle","sd_err") #--------------------------------- # Charts #--------------------------------- if (!is.null(NamePDF)) { pdf(paste("Charts/",NamePDF,".pdf",sep=""), width=5, height=6) rows=6 cols=3 Countries=sort(unique(Country)) par(mfrow = c(rows,cols), mar=c(2,2,2,0.5), oma=c(2,2,2,2), cex=0.5, cex.main = 1, cex.axis = 1, mgp = c(2,0.5,0)) for (i in 1:length(Countries)) { Date_i=Date[which(Country==Countries[i])] vars_i=vars[which(Country==Countries[i]),] plot(Date_i,vars_i[,1],main=paste(colnames(vars)[1],"and coint (red)"),type="l",xlab="",ylab="", ylim=c(min(vars_i[,1],vars_i[,ncol(vars)]),max(vars_i[,1],vars_i[,ncol(vars)]))) lines(Date_i,vars_i[,ncol(vars)],col="red") k=1 for (j in 2:(ncol(vars)-1)) { plot(Date_i,vars_i[,j],main=colnames(vars)[j],type="l",xlab="",ylab="") if ((k+rows*cols-1)%%(rows*cols) == 0) title(Countries[i],outer=TRUE,col.main="blue",adj=0.03) k=k+1 } for (m in (j+1):(ceiling(j/(rows*cols))*rows*cols)) {plot.new()} } Resid=vars$resid/sd(vars_i$resid) for (i in 1:length(Countries)) { plot(Date[which(Country==Countries[i])], Resid[which(Country==Countries[i])] ,type="l",xlab="",ylab="",main=Countries[i], ylim=c(min(Resid),max(Resid))) if (i==floor(i/(rows*cols))*rows*cols+1) {title("STANDARDISED RESIDUALS",outer=TRUE,col.main="red",adj=0.03)} } # Data availability chart par(mfrow = c(1,1),mar=c(5,2,4,2)) Availability=data.frame(C=Country,D=Date,ones=rep(1,length(Country))) Availability=reshape(Availability, idvar = "D", timevar = "C", direction = "wide") Availability[is.na(Availability)]=0 colnames(Availability)=gsub("ones.","",colnames(Availability)) rownames(Availability)=Availability[,1];Availability=Availability[,-1] Availability=Availability[rev(order(rownames(Availability))),] Availability=Availability[,order(colnames(Availability))] image(as.matrix(t(Availability)), col=0:1,axes=FALSE, main="Data availability") axis(1, at = (1:ncol(Availability)-1)/(ncol(Availability)-1), labels=colnames(Availability), las=2,cex=.3) axis(2, at = (1:nrow(Availability))/nrow(Availability), labels=format(as.Date(rownames(Availability)), format = "%m/%Y"), las=2,cex=.3) # Coint relation par(mfrow = c(rows,cols), mar=c(2,2,2,0.5), oma=c(2,2,2,2), cex=0.5, cex.main = 1, cex.axis = 1, mgp = c(2,0.5,0)) for (i in 1:length(Countries)) { plot(Date[which(Country==Countries[i])],vars[which(Country==Countries[i]),1],main=Countries[i],type="l",xlab="",ylab="", ylim=c(min(vars[which(Country==Countries[i]),1],vars[which(Country==Countries[i]),ncol(vars)]), max(vars[which(Country==Countries[i]),1],vars[which(Country==Countries[i]),ncol(vars)]))) lines(Date[which(Country==Countries[i])],vars[which(Country==Countries[i]),ncol(vars)],col="red") } dev.off() } OptimPMG =OutputList }