{"word_count": 440, "character_count": 5647}
A few months ago, I posted a note with some home made codes for quantile regression... there was something odd on the output, but it was because there was a (small) mathematical problem in my equation. So since I should teach those tomorrow, let me fix them.MedianConsider a sample [latex]\{y_1,\cdots,y_n\}[/latex]. To compute the median, solve[latex display="true"]\min_\mu \left\lbrace\sum_{i=1}^n|y_i-\mu|\right\rbrace[/latex]which can be solved using linear programming techniques. More precisely, this problem is equivalent to[latex display="true"]\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^na_i+b_i\right\rbrace[/latex]with [latex]a_i,b_i\geq 0[/latex] and [latex]y_i-\mu=a_i-b_i[/latex], [latex]\forall i=1,\cdots,n[/latex]. Heuristically, the idea is to write [latex]y_i=\mu+\varepsilon_i[/latex], and then define [latex]a_i[/latex]'s and [latex]b_i[/latex]'s so that [latex]\varepsilon_i=a_i-b_i[/latex] and [latex]|\varepsilon_i|=a_i+b_i[/latex], i.e. [latex display="true"]a_i=(\varepsilon_i)_+=\max\lbrace0,\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i&#62;0}[/latex]and[latex display="true"]b_i=(-\varepsilon_i)_+=\max\lbrace0,-\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i&#60;0}[/latex]denote respectively the positive and the negative parts.Unfortunately (that was the error in my previous post), the expression of linear programs is[latex display="true"]\min_{\mathbf{z}}\left\lbrace\boldsymbol{c}^\top\mathbf{z}\right\rbrace\text{ s.t. }\boldsymbol{A}\mathbf{z}=\boldsymbol{b},\mathbf{z}\geq\boldsymbol{0}[/latex]In the equation above, with the [latex]a_i[/latex]'s and [latex]b_i[/latex]'s, we're not far away. Except that we have [latex]\mu\in\mathbb{R}[/latex], while it should be positive. So similarly, set [latex]\mu=\mu^+-\mu^-[/latex] where [latex]\mu^+=(\mu)_+[/latex] and [latex]\mu^-=(-\mu)_+[/latex].Thus, let[latex display="true"]\mathbf{z}=\big(\mu^+;\mu^-;\boldsymbol{a},\boldsymbol{b}\big)^\top\in\mathbb{R}_+^{2n+2}[/latex]and then write the constraint as [latex]\boldsymbol{A}\mathbf{z}=\boldsymbol{b}[/latex] with [latex]\boldsymbol{b}=\boldsymbol{y}[/latex] and [latex display="true"]\boldsymbol{A}=\big[\boldsymbol{1}_n;-\boldsymbol{1}_n;\mathbb{I}_n;-\mathbb{I}_n\big][/latex]And for the objective function[latex display="true"]\boldsymbol{c}=\big(\boldsymbol{0},\boldsymbol{1}_n,-\boldsymbol{1}_n\big)^\top\in\mathbb{R}_+^{2n+2}[/latex]To illustrate, consider a sample from a lognormal distribution,n = 101 set.seed(1)y = rlnorm(n)median(y)[1] 1.077415For the optimization problem, use the matrix form, with [latex]3n[/latex] constraints, and [latex]2n+1[/latex] parameters,library(lpSolve) X = rep(1,n) A = cbind(X, -X, diag(n), -diag(n))b = yc = c(rep(0,2), rep(1,n),rep(1,n))equal_type = rep("=", n) r = lp("min", c,A,equal_type,b)head(r$solution,1)[1] 1.077415It looks like it's working well...QuantileOf course, we can adapt our previous code for quantilestau = .3quantile(y,tau)      30% 0.6741586The linear program is now[latex display="true"]\min_{q^+,q^-,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbrace[/latex]with [latex]a_i,b_i,q^+,q^-\geq 0[/latex] and [latex]y_i=q^+-q^-+a_i-b_i[/latex], [latex]\forall i=1,\cdots,n[/latex]. The R code is nowc = c(rep(0,2), tau*rep(1,n),(1-tau)*rep(1,n))r = lp("min", c,A,equal_type,b)head(r$solution,1)[1] 0.6741586So far so good...Quantile RegressionConsider the following dataset, with rents of flat, in a major German city, as function of the surface, the year of construction, etc.base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)The linear program for the quantile regression is now[latex display="true"]\min_{\boldsymbol{\beta}^+,\boldsymbol{\beta}^-,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbrace[/latex]with [latex]a_i,b_i\geq 0[/latex] and [latex display="true"]y_i=\boldsymbol{x}^\top[\boldsymbol{\beta}^+-\boldsymbol{\beta}^-]+a_i-b_i[/latex][latex]\forall i=1,\cdots,n[/latex] and [latex]\beta_j^+,\beta_j^-\geq 0[/latex] [latex]\forall j=0,\cdots,k[/latex]. So use hererequire(lpSolve) tau = .3n=nrow(base)X = cbind( 1, base$area)y = base$rent_euroK = ncol(X)N = nrow(X)A = cbind(X,-X,diag(N),-diag(N))c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))b = base$rent_euroconst_type = rep("=",N)r = lp("min",c,A,const_type,b)beta = r$sol[1:K] -  r$sol[(1:K+K)]beta[1] 148.946864   3.289674Of course, we can use R function to fit that modellibrary(quantreg)rq(rent_euro~area, tau=tau, data=base)Coefficients:(Intercept)        area  148.946864    3.289674Here again, it seems to work quite well. We can use a different probability level, of course, and get a plotplot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")),     ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5)sf=0:250yr=r$solution[2*n+1]+r$solution[2*n+2]*sflines(sf,yr,lwd=2,col="blue")tau = .9r = lp("min",c,A,const_type,b)tail(r$solution,2)[1] 121.815505   7.865536yr=r$solution[2*n+1]+r$solution[2*n+2]*sflines(sf,yr,lwd=2,col="blue")And we can adapt the later to multiple regressions, of course,X = cbind(1,base$area,base$yearc)K = ncol(X)N = nrow(X)A = cbind(X,-X,diag(N),-diag(N))c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))b = base$rent_euroconst_type = rep("=",N)r = lp("min",c,A,const_type,b)beta = r$sol[1:K] -  r$sol[(1:K+K)]beta[1] -5542.503252     3.978135     2.887234to be compared withlibrary(quantreg)rq(rent_euro~ area + yearc, tau=tau, data=base)Coefficients: (Intercept)         area        yearc -5542.503252     3.978135     2.887234 Degrees of freedom: 4571 total; 4568 residual