R软件程序脚本：

rq.fit.panel<-function(X,y,s,w=c(.1,.25,.5,.25,.1),taus=c(0.1,0.25,0.5,0.75,0.9),lambda = 1){

# prototype function for panel data fitting of QR models

# the matrix X is assumed to contain an intercept

# 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 homogneous settings it should be the ratio of the scale parameters

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,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 * rep(1/2,n))

rq.fit.sfn(D,y,rhs=a)

}

of the fixed effects and the idiocyncratic 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 colums 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.

R统计软件操作程序：

library(“quantreg”)

source(“rq.fit.panel”)

m <- 12

n <- 8

s <- rep(1:n,rep(m,n))

data<-read.csv("zhongbu.csv")

x<-data$LNVEDU

X<-cbind(1,x)

y<-data$LNVGDP

fit<- rq.fit.panel(X,y,s)

fit

n1=length(y)

b=2000 #Number of bootstrap samples

ncoef=5*ncol(X)+n

boot_bhat=matrix(NA,b, ncoef)

block_length =8 #代表块状长度

num_blocks = n1/block_length #n/block_length

Indices = seq(1:n1) # All of the indices from 1 to n

Indices = matrix(Indices,block_length,num_blocks)

for (i in 1:b){ #Number of bootstrap samples

randblock =sample(seq(1:num_blocks),num_blocks,replace = TRUE) # Choose which blocks to use

Ind_sim = Indices[,randblock] #Find which data are in each block

Ind_sim = c(Ind_sim)

p <- length(bhat)

rdf <- n1 - p

vnames<- dimnames(x)[[2]]

coef <- array(bhat, c(p, 4))

dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value","Pr(>|t|)"))

coef[, 2] <- serr

coef[, 3] <- coef[, 1]/coef[, 2]

coef[, 4] <- if (rdf > 0) 2 * (1 - pt(abs(coef[, 3]), rdf))

coef