################################################################
# Panel Quantile Regression 
# Lamarche, Journal of Econometrics (2010)
##########################################################################
rq.fit.panel <- function(X,y,s,w=c(.25,.5,.25),taus=(1:3)/4,lambda = 1){
require(SparseM)
require(quantreg)
K <- length(w)
if(K != length(taus))
stop("length of w and taus must match")
X <- as.matrix(X)
p <- ncol(X)
n <- length(levels(as.factor(s)))
N <- length(y)
if(N != length(s) || N != nrow(X))
 stop("dimensions of y,X,s must match")
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
Fidelity <- cbind(as(w,"matrix.diag.csr") %x% X,as.matrix.csr(w) %x% Z)
Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
D <- rbind(Fidelity,Penalty)
y <- c(w %x% y,rep(0,n))
a <- c((w*(1-taus)) %x% (t(X)%*%rep(1,N)),sum(w*(1-taus)) * (t(Z) %*% rep(1,N)) + lambda/2 * rep(1,n))
rq.fit.sfn(D,y,rhs=a)
}
#
################################################################
# Quantile Regression IV Panel Data Functions
# Harding and Lamarche, Economics Letters 104 (2009), 133-135 
##########################################################################
# Functions
##########################################################################
rq.fit.fe <- function(X,y,w=1,taus=tau)
	{
	require(SparseM)
	require(quantreg)
	K <- length(w)
	if(K != length(taus))
	stop("length of w and taus must match")
	X <- as.matrix(X)
	p <- ncol(X)
	n <- length(levels(as.factor(s)))
	N <- length(y)
	if( N != nrow(X))
	 stop("dimensions of y,X must match")
	D <- cbind(as(w,"matrix.diag.csr") %x% X)
	y <- c(w %x% y)
	a <- c((w*(1-taus)) %x% (t(X)%*%rep(1,N)))
	rq.fit.sfn(D,y,rhs=a)
}
#############################################
rq.fit.ivpanel <- function(d,exo,iv,y,s,tau)
	{
	# IV Quantile Regression with Fixed Effects
	# s is an strata indicator
	# d is the endogenous variable
	# exo are the exogenous variables, no intercept.
	# iv is the intrument
	Z <- model.matrix(s~as.factor(s)-1)[,-1]
	exo <- cbind(exo,Z)
	X <- cbind(exo,1)
	x <- cbind(d,X)
	w <- cbind(iv,X)
	ww <- t(w) %*% w
	ww.inv <- ginv(as.matrix(ww))
	wd <- t(w)%*%d
	dhat <- w%*%ww.inv%*%wd
	PSI <- cbind(dhat,X)
	PSI1 <- cbind(d,X)
	coef <- rq.fit.fe(PSI,y,tau=tau)$coef
	resid <- y - PSI1%*%coef
	mu1 <- mean(resid)
	sigma1 <- var(resid)
	c <- ((1-tau)*tau)/(dnorm(qnorm(tau,mu1,sqrt(sigma1)),mu1,sqrt(sigma1))^2)
	PSIinv <- diag(length(coef))
	PSIinv <- backsolve(qr(x)$qr[1:length(coef), 1:length(coef)], PSIinv)
	PSIinv <- PSIinv %*% t(PSIinv)
	vc1 <- c*PSIinv
	std <- sqrt((length(y)/100)*diag(vc1))
	alpha <- seq(coef[1]-2*std[1],coef[1]+2*std[1],by=std[1]/20)
	z <- cbind(dhat,X)
	betas <- matrix(NA,dim(z)[2],length(alpha))
	g <- matrix(NA,length(alpha),1)
	for (i in 1:length(alpha)){
		ya <- y - alpha[i]*d
		betas[,i] <- rq.fit.fe(z,ya,tau=tau)$coef
		g[i,] <- max(svd(betas[1:dim(dhat)[2],i])$d)
		}
	I <- which.min(g[,1])
	param1 <- alpha[I]
	est1 <- betas[(dim(dhat)[2]+1):dim(z)[2],I]
	c(param1,est1)
}

PQ <- function(h, id){
if(is.vector(h))
h <- matrix(h, ncol = 1)
Ph <- unique(id)
Ph <- cbind(Ph, table(id))
for(i in 1:ncol(h))
Ph <- cbind(Ph, tapply(h[, i], id, mean))
is <- tapply(id, id)
Ph <- Ph[is, - (1:2)]
Qh <- h - Ph
list(Ph=as.matrix(Ph), Qh=as.matrix(Qh), is=is)
}