https://github.com/cran/quantreg
Tip revision: 3ad414306dc4e0be8f04abb2f3abff2cb11a429b authored by Roger Koenker on 13 December 2019, 10:10:02 UTC
version 5.54
version 5.54
Tip revision: 3ad4143
Panel.R
library(quantreg)
rq.fit.panel <- function(X,y,s,w=c(1/3,1/3,1/3),taus=c(0.25,0.5,0.75),lambda = 0){
# prototype function for fixed effect panel data fitting of QR models
# the vector s is a strata indicator assumed (so far) to be a one-way layout
# NB:
# 1. The value of the shrinkage parameter lambda is an open research problem in
# the simplest homogeneous settings it should be the ratio of the scale parameters
# of the fixed effects and the idiosyncratic errors
# 2. On return the coefficient vector has m*p + n elements where m is the number
# quantiles being estimated, p is the number of columns of X, and n is the
# number of distinct values of s. The first m*p coefficients are the
# slope estimates, and the last n are the "fixed effects"
# 3. Like all shrinkage (regularization) estimators, asymptotic inference is somewhat
# problematic... so the bootstrap is the natural first resort.
require(SparseM)
require(quantreg)
m <- length(w)
if(m != 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, cbind(w) %x% Z)
Penalty <- cbind(as.matrix.csr(0,n,m*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 * rep(1,n))
rq.fit.sfn(D,y,rhs=a)
}
n<-3
T<-50
nT<-n*T
u1<-rnorm(T)
u2<-rnorm(T)
u3<-rnorm(T)
x1<-rnorm(T,1,0.85)
x2<-rnorm(T,4,1)
x3<-rnorm(T,7,1)
beta1<-1
beta2<-0.5
y1<- 0+beta1*x1+(beta2*x1)*u1
y2<- 4+beta1*x2+(beta2*x2)*u2
y3<- 8+beta1*x3+(beta2*x3)*u3
plot(c(0,9), c(0,25), type='n', xlab=expression(x[it]), ylab=expression(y[it]))
points(x1,y1,pch=15)
points(x2,y2,pch=15,col="blue")
points(x3,y3,pch=15,col="red")
legend(1,17,paste("i = ",1:3,sep = ""),pch = 15, col = c("black","blue","red"))
ya<-c(y1,y2,y3)
xa<-c(x1,x2,x3)
# Naive cross-section QR
taus <- c(.25,0.5,.75)
xx <- seq(min(xa),max(xa),0.25)
f <- coef(rq(ya~xa,tau=taus))
yy <- cbind(1,xx)%*%f
for(i in 1:3)
lines(xx,yy[,i],col = "grey")
# Fixed effect QR
s <- rep(1:n,rep(T,n))
fp<-rq.fit.panel(xa,ya,s)$coef
bhat <- fp[1:3]
fehat <- fp[4:6]
xx1 <- seq(min(x1),max(x1),0.25)
for(i in 1:3){
yy1 <- fehat[1] + bhat[i] * xx1
lines(xx1,yy1,col = "black")
}
xx2 <- seq(min(x2),max(x2),0.25)
for(i in 1:3){
yy2 <- fehat[2] + bhat[i] * xx2
lines(xx2,yy2,col = "blue")
}
xx3 <- seq(min(x3),max(x3),0.25)
for(i in 1:3){
yy3 <- fehat[3] + bhat[i] * xx3
lines(xx3,yy3,col = "red")
}