rm(list = ls())
#   
date()
sessionInfo()
#########################################################################
# MONTE CARLO EXPERIMENT 4
##########################################################################
##########################################################################
# Libraries
##########################################################################
library(quantreg)
library(splines)
library(SparseM)
#library(plm)
library(MASS)
##########################################################################
# Functions
##########################################################################
source("rq.fit.panel.all.revised.R")
#########################################################################
# Data
##########################################################################
n <- c(100,100,100,200,200,200)
m <- c(50,100,200,50,100,200)

#n <- c(100,100,100)
#m <- c(50,100,200)

S = 200

R <- 400

lambda = 0.5
beta1 = 1
beta2 = 0.5
rho_x = 0.8
rho_f = 0.9
rho_x_scale = sqrt(1-rho_x^2)
rho_f_scale = sqrt(1-rho_f^2)
z0 = 0
f0 = 0

delta = 0.1

Egamma = 0.5
Vgamma = 1
Evartheta = 0.5
Vvartheta = 1
Emu = 0.5
Vmu = 1

GeoLambda <- array(NA,c(length(S)))

MG.F00 <- array(0,c(3,R,length(n)))
MG.G00 <- array(0,c(3,R,length(n)))
DYN.F00 <- array(0,c(3,R,length(n)))
DYN.G00 <- array(0,c(3,R,length(n)))

set.seed(14)

for (ii in 1:length(n))
	{
#ii=i=1
	s <- rep(1:n[ii],each=m[ii])
	st <- rep(1:m[ii],n[ii])

for(i in 1:R){
	z1it <- z2it <- y0 <- Y <- YL1 <- YL2 <- u <- NULL

	deltai <- rep(runif(n[ii],0,0.2),each=(m[ii]+S))
	sigmai <- rep(runif(n[ii],0.9,1.1),each=(m[ii]+S))

	for(t in 1:(m[ii]+S))
	{
	# f
			if(t==1){f1 <- f0}
			if(t!=1){f1 <- rbind(f1,rho_f * f1[t-1] + rho_f_scale * rnorm(1))}

			if(t==1){f2 <- f0}
			if(t!=1){f2 <- rbind(f2,rho_f * f2[t-1] + rho_f_scale * rnorm(1))}


	}

	f1 <- rep(f1,n[ii])
	f2 <- rep(f2,n[ii])

	for(j in 1:n[ii])
	{
		for(t in 1:(m[ii]+S))
		{
		if(t==1){z1 <- z0}
		if(t!=1){z1 <- rbind(z1,rho_x * z1[t-1] + rho_x_scale * rnorm(1))}

		if(t==1){z2 <- z0}
		if(t!=1){z2 <- rbind(z2,rho_x * z2[t-1] + rho_x_scale * rnorm(1))}

		}
		z1it <- c(z1it,z1)
		z2it <- c(z2it,z2)

	}

	mu <- sqrt(Vmu)*rnorm(n[ii],Emu)
	mu <- rep(mu,each=(m[ii]+S))

	vartheta <- sqrt(Vvartheta)*rnorm(n[ii],Evartheta)
	vartheta <- rep(vartheta,each=(m[ii]+S))

	x1 <- mu + vartheta * f1 + z1it
	x2 <- mu + vartheta * f2 + z2it

	X1 <- matrix(x1,m[ii]+S,n[ii])
	X1S <- X1[c(1:S),]
	X1T <- X1[-c(1:S),]

	X2 <- matrix(x2,m[ii]+S,n[ii])
	X2S <- X2[c(1:S),]
	X2T <- X2[-c(1:S),]

	x1bar <- rep(apply(X1T,2,mean),each=m[ii])
	x2bar <- rep(apply(X2T,2,mean),each=m[ii])

	gamma1 <- sqrt(Vgamma)*rnorm(n[ii],Egamma)
	gamma1 <- rep(gamma1,each=(m[ii]+S))
	gamma2 <- sqrt(Vgamma)*rnorm(n[ii],Egamma)
	gamma2 <- rep(gamma2,each=(m[ii]+S))

	Gamma1 <- matrix(gamma1,m[ii]+S,n[ii])
	gamma1T <- c(Gamma1[-c(1:S),])
	Gamma2 <- matrix(gamma2,m[ii]+S,n[ii])
	gamma2T <- c(Gamma2[-c(1:S),])

	U <- sigmai * (1 + deltai * x1) * rnorm( (m[ii]+S) * n[ii] )
 	
	epsilon <- gamma1 * f1 + gamma2 * f2 + U
	Epsilon <- matrix(epsilon,m[ii]+S,n[ii])
	epsilonS <- Epsilon[c(1:S),]
	epsilonT <- Epsilon[-c(1:S),]

	F1 <- matrix(f1,m[ii]+S,n[ii])
	F1S <- F1[c(1:S),]
	F1T <- F1[-c(1:S),]

	F2 <- matrix(f2,m[ii]+S,n[ii])
	F2S <- F2[c(1:S),]
	F2T <- F2[-c(1:S),]

	f1bar <- rep(apply(F1T,2,mean),each=m[ii])
	f2bar <- rep(apply(F2T,2,mean),each=m[ii])

	U <- matrix(U,m[ii]+S,n[ii])
	US <- U[c(1:S),]
	UT <- U[-c(1:S),]

	ubar <- rep(apply(UT,2,mean),each=m[ii])

	alpha <- x1bar + gamma1T * f1bar + ubar + rep(rnorm(n[ii]),each=m[ii])
	Alpha <- matrix(alpha,m[ii],n[ii])

	beta1ix <- rep(runif(n[ii],-0.25,0.25),each=m[ii])
	beta1i <- beta1 + beta1ix
	betauniq <- beta1i[st==1]

	for (oi in 1:S){ GeoLambda[oi] <- lambda^(oi-1)}

	for (sj in 1:n[ii])
	{
		y0[sj] <- (1-lambda)^(-1) * Alpha[1,sj] +  betauniq[sj] * sum(GeoLambda * rev(X1S[,sj])) + beta2 * sum(GeoLambda * rev(X2S[,sj])) + sum(GeoLambda * rev(epsilonS[,sj]))
	}

	for (j in 1:n[ii])
	{
		for(t in 1:m[ii])
		{
			if(t==1){y <-  Alpha[t,j] + lambda * y0[j] + betauniq[j] * X1T[t,j] + beta2 * X2T[t,j] + epsilonT[t,j]}
			if(t!=1){y <- rbind(y, Alpha[t,j] + lambda * y[t-1] + betauniq[j] * X1T[t,j] + beta2 * X2T[t,j] + epsilonT[t,j])}
		}
		Y <- c(Y,y)
	}

	for (j in 1:n[ii])
	{
		Ysj <- cbind(Y,s)[s==j,]
		YL1 <- c(YL1,c(NA,Ysj[-dim(Ysj)[1],1]))
	}	

	for (j in 1:n[ii])
	{
		YL1sj <- cbind(YL1,s)[s==j,]
		YL2 <- c(YL2,c(NA,YL1sj[-dim(YL1sj)[1],1]))
	}

	YY <- cbind(Y,YL1)
	X <- cbind(c(X1T),c(X2T))
	mYX <- PQ(cbind(YY,X),st)$Ph
	YX <- cbind(s,YY,X,mYX)
	YX <- na.omit(YX)
	
	BetasHP <- matrix(0,6,n[ii])

	#

	for (isi in 1:n[ii]){
		fit <- rq(YX[YX[,1]==isi,2]~YX[YX[,1]==isi,3]+YX[YX[,1]==isi,4]+YX[YX[,1]==isi,5]+YX[YX[,1]==isi,6]+YX[YX[,1]==isi,7]+YX[YX[,1]==isi,8]+YX[YX[,1]==isi,9],tau=0.5)
		BetasHP[1,isi] <- fit$coef[2]
		BetasHP[2,isi] <- fit$coef[3]
		BetasHP[3,isi] <- fit$coef[3]/(1-fit$coef[2])

		fit <- rq(YX[YX[,1]==isi,2]~YX[YX[,1]==isi,3]+YX[YX[,1]==isi,4]+YX[YX[,1]==isi,5]+YX[YX[,1]==isi,6]+YX[YX[,1]==isi,7]+YX[YX[,1]==isi,8]+YX[YX[,1]==isi,9],tau=0.25)
		BetasHP[4,isi] <- fit$coef[2]
		BetasHP[5,isi] <- fit$coef[3]
		BetasHP[6,isi] <- fit$coef[3]/(1-fit$coef[2])

		}

	#
	# Hashem
	#
	
	MG.F00[1,i,ii] <- mean(BetasHP[1,])
	MG.F00[2,i,ii] <- mean(BetasHP[2,])
	MG.F00[3,i,ii] <- mean(BetasHP[3,])

	MG.G00[1,i,ii] <- mean(BetasHP[4,])
	MG.G00[2,i,ii] <- mean(BetasHP[5,])
	MG.G00[3,i,ii] <- mean(BetasHP[6,])

	#
	# Galvao
	#
	
	YYY <- cbind(Y,YL1,YL2)
	YYX <- cbind(s,YYY,X)
	YYX <- na.omit(YYX)

	fit <- rq.fit.ivpanel(YYX[,3],cbind(YYX[,5],YYX[,6]),YYX[,4],YYX[,2],YYX[,1],tau=0.5)

	DYN.F00[1,i,ii] <- fit[1]
	DYN.F00[2,i,ii] <- fit[2]
	DYN.F00[3,i,ii] <- fit[2]/(1-fit[1])

	fit <- rq.fit.ivpanel(YYX[,3],cbind(YYX[,5],YYX[,6]),YYX[,4],YYX[,2],YYX[,1],tau=0.25)

	DYN.G00[1,i,ii] <- fit[1]
	DYN.G00[2,i,ii] <- fit[2]
	DYN.G00[3,i,ii] <- fit[2]/(1-fit[1])

	}

}

#
# Functions to obtain bias and rmse
#

bias <- function(x,beta,delta,tau){(mean(x) - (beta+delta*qnorm(tau)))}
se <- function(x){sqrt(var(x))}
rmse <- function(x,beta,delta,tau){ ( (mean(x) - (beta+delta*qnorm(tau)))^2 + var(x) )^.5}

bias.theta <- function(x,beta,delta,lambda,tau){
	theta <- (beta+delta*qnorm(tau))/(1-lambda)
	(mean(x) - theta)
	}
se <- function(x){sqrt(var(x))}
rmse.theta <- function(x,beta,delta,lambda,tau){ 
	theta <- (beta+delta*qnorm(tau))/(1-lambda)
	( (mean(x) - theta)^2 + var(x) )^.5
}

#
# Generate a Table
#
B <- array(NA,c(12,6))
colnames(B) <- c("DYN", "QMG", "DYN", "QMG", "DYN", "QMG")

j=1

for (i in 1:length(n))
{

	j = (2 * i) - 1
	B[j,3] <- bias(DYN.F00[1,,i],lambda,delta=0,tau=0.5)
	B[j+1,3] <- rmse(DYN.F00[1,,i],lambda,delta=0,tau=0.5)
	B[j,4] <- bias(MG.F00[1,,i],lambda,delta=0,tau=0.5)
	B[j+1,4] <- rmse(MG.F00[1,,i],lambda,delta=0,tau=0.5)

	B[j,1] <- bias(DYN.F00[2,,i],beta1,delta=delta,tau=0.5)
	B[j+1,1] <- rmse(DYN.F00[2,,i],beta1,delta=delta,tau=0.5)
	B[j,2] <- bias(MG.F00[2,,i],beta1,delta=delta,tau=0.5)
	B[j+1,2] <- rmse(MG.F00[2,,i],beta1,delta=delta,tau=0.5)

	B[j,5] <- bias.theta(DYN.F00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.5)
	B[j+1,5] <- rmse.theta(DYN.F00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.5)
	B[j,6] <- bias.theta(MG.F00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.5)
	B[j+1,6] <- rmse.theta(MG.F00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.5)

}

#
round(B,4)
#

#
# Generate a Table
#

Bt <- array(NA,c(12,6))
colnames(Bt) <- c("DYN", "QMG", "DYN", "QMG", "DYN", "QMG")

j=1

for (i in 1:length(n))
{

	j = (2 * i) - 1
	Bt[j,3] <- bias(DYN.G00[1,,i],lambda,delta=0,tau=0.25)
	Bt[j+1,3] <- rmse(DYN.G00[1,,i],lambda,delta=0,tau=0.25)
	Bt[j,4] <- bias(MG.G00[1,,i],lambda,delta=0,tau=0.25)
	Bt[j+1,4] <- rmse(MG.G00[1,,i],lambda,delta=0,tau=0.25)

	Bt[j,1] <- bias(DYN.G00[2,,i],beta1,delta=delta,tau=0.25)
	Bt[j+1,1] <- rmse(DYN.G00[2,,i],beta1,delta=delta,tau=0.25)
	Bt[j,2] <- bias(MG.G00[2,,i],beta1,delta=delta,tau=0.25)
	Bt[j+1,2] <- rmse(MG.G00[2,,i],beta1,delta=delta,tau=0.25)

	Bt[j,5] <- bias.theta(DYN.G00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.25)
	Bt[j+1,5] <- rmse.theta(DYN.G00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.25)
	Bt[j,6] <- bias.theta(MG.G00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.25)
	Bt[j+1,6] <- rmse.theta(MG.G00[3,,i],beta1,delta=delta,lambda=lambda,tau=0.25)

}

#
round(Bt,4)

ns <- c(rep(n,each=2))
ms <- c(rep(m,each=2))
round(cbind(ns,ms,B),4)
round(cbind(ns,ms,Bt),4)

#
# Final Tables for the Normal Case:
#

Table1.beta <- round(cbind(ns,ms,B[,1:2],Bt[,1:2]),4)
Table1.beta 

Table1.lambda.theta <- round(cbind(ns,ms,B[,-c(1:2)],Bt[,-c(1:2)]),4)
Table1.lambda.theta