Formular regressão quantílica como problema de programação linear?

Como formular regressão quantílica como um problema de programação linear? Ao olhar para o problema mediano dos quantis, sei que é

$$\begin{align} \text{minimize } & \sum_{i=1}^n |\beta_0 + X_i \beta_1-Y_i|\\ \text{transforms into } & \\ \text{minimize } & \sum_{i=1}^n e_i\\ \text{s.t.} & \\ & e_i\geq \beta_0 + X_i\beta_{1}-Y_i\\ & e_i\geq -(\beta_0 + X_i\beta_{1}-Y_i) \end{align}$$ mas como transformar a minimização para outros quantis?

Você usa o estimador de regressão quantil

$$\hat \beta(\tau) := \arg \min_{\theta \in \mathbb R^K} \sum_{i=1}^N \rho_\tau(y_i - \mathbf x_i^\top \theta).$$

onde $\tau \in (0,1)$ é constante escolhido de acordo com o qual o quantil precisa ser estimado e a função $\rho_\tau(.)$ é definida como

$$\rho_\tau(r) = r(\tau - I(r<0)).$$

Para ver o objetivo de Considere primeiro que ele leva os resíduos como argumentos, quando são definidos como . A soma do problema de minimização pode, portanto, ser reescrita como$\rho_\tau(.)$$\epsilon_i =y_i - \mathbf x_i^\top \theta$

$$\sum_{i=1}^N \rho_\tau(\epsilon_i) =\sum_{i=1}^N \tau \lvert \epsilon_i \lvert I[\epsilon_i \geq 0] + (1-\tau) \lvert \epsilon_i \lvert I[\epsilon_i < 0]$$

de modo que resíduos positivos associados à observação acima da linha de regressão quantílica sugerida recebem o peso de enquanto resíduos negativos associados à observação abaixo da linha de regressão quantílica sugerida são ponderados com .$y_i$$\mathbf x_i^\top \theta$$\tau$$y_i$$\mathbf x_i^\top \theta$$(1-\tau)$

Intuitivamente:

Com resíduos positivos e negativos são "punidos" com o mesmo peso e um número igual de observação está acima e abaixo da "linha" no ideal, de modo que a linha seja a regressão mediana "linha".$\tau=0.5$$\mathbf x_i^\top \hat \beta$

Quando cada resíduo positivo é ponderado 9 vezes o de um resíduo negativo com peso e, portanto, é ideal para todas as observações acima da "linha" aproximadamente 9 ser colocado abaixo da linha. Portanto, a "linha" representa o quantil 0,9. (Para uma declaração exata disso, consulte THM. 2.2 e Corollary 2.1 em Koenker (2005) "Regressão quantílica")$\tau=0.9$$1-\tau= 0.1$$\mathbf x_i^\top \hat \beta$

Os dois casos são ilustrados nessas parcelas. Painel esquerdo e painel direito .$\tau=0.5$$\tau=0.9$

Os programas lineares são predominantemente analisados ​​e resolvidos usando o formulário padrão

$$(1) \ \ \min_z \ \ c^\top z \ \ \mbox{subject to } A z = b , z \geq 0$$

Para chegar a um programa linear na forma padrão, o primeiro problema é que em um programa (1) todas as variáveis ​​sobre as quais a minimização é realizada devem ser positivas. Para isso, os resíduos são decompostos em partes positivas e negativas usando variáveis ​​de folga:$z$

$$\epsilon_i = u_i - v_i$$

onde parte positiva e é a parte negativa. A soma dos pesos residuais atribuídos pela função de verificação é então vista como$u_i = \max(0,\epsilon_i) = \lvert \epsilon_i \lvert I[\epsilon_i \geq 0]$$v_i = \max(0,-\epsilon_i) =\lvert \epsilon_i \lvert I[\epsilon_i < 0]$

$$\sum_{i=1}^N \rho_\tau(\epsilon_i) = \sum_{i=1}^N \tau u_i + (1-\tau) v_i = \tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v,$$

onde e e é vetor todas as coordenadas igual a .$u = (u_1,...,u_N)^\top$$v=(v_1,...,v_N)^\top$$\mathbf 1_N$$N \times 1$$1$

Os resíduos devem satisfazer as restrições que$N$

$$y_i - \mathbf x_i^\top\theta = \epsilon_i = u_i - v_i$$

Isso resulta na formulação como um programa linear

$$\min_{\theta \in \mathbb R^K,u\in \mathbb R_+^N,v\in \mathbb R_+^N}\{ \tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v\lvert y_i= \mathbf x_i\theta + u_i - v_i, i=1,...,N\},$$

como indicado na equação de Koenker (2005) "Regressão quantílica", página 10 (1.20).

No entanto, é perceptível que ainda não está restrito a ser positivo, conforme exigido no programa linear no formulário padrão (1). Portanto, novamente se decompõe em parte positiva e negativa.$\theta\in \mathbb R$

$$\theta = \theta^+ - \theta^-$$

onde novamente é parte positiva e é parte negativa. As restrições podem então ser escritas como$\theta^+=max(0,\theta)$$\theta^- = \max(0,-\theta)$$N$

$$\mathbf y:= \begin{bmatrix} y_1 \\ \vdots \\ y_N\end{bmatrix} = \begin{bmatrix} \mathbf x_1^\top \\ \vdots \\ \mathbf x_N^\top \end{bmatrix}(\theta^+ - \theta^-) + \mathbf I_Nu - \mathbf I_Nv ,$$

onde .$\mathbf I_N = diag\{\mathbf 1_N\}$

Em seguida, defina e a matriz de design armazenando dados em variáveis ​​independentes como$b:=\mathbf y$$\mathbf X$

$$\mathbf X := \begin{bmatrix} \mathbf x_1^\top \\ \vdots \\ \mathbf x_N^\top \end{bmatrix}$$

Para reescrever restrição:

$$b= \mathbf X(\theta^+ - \theta^-) + \mathbf I_N u- \mathbf I_N v= [\mathbf X , -\mathbf X , \mathbf I_N ,\mathbf - \mathbf I_N] \begin{bmatrix} \theta^+ \\ \theta^- \\ u \\ v\end{bmatrix}$$

Definir a matriz$(N \times 2K + 2N )$

$$A := [\mathbf X , -\mathbf X , \mathbf I_N ,\mathbf - \mathbf I_N]$$ e introduza e como variáveis ​​sobre as quais minimizar, para que façam parte de para obter$\theta^+$$\theta^-$$z$

$$b = A \begin{bmatrix} \theta^+ \\ \theta^- \\ u \\ v\end{bmatrix} = Az$$

Como e afetam apenas o problema de minimização através da restrição a da dimensão deve ser introduzida como parte do vetor coeficiente que pode ser definido adequadamente como$\theta^+$$\theta^-$$\mathbf 0$$2K\times 1$$c$

$$c = \begin{bmatrix}\mathbf 0 \\ \tau \mathbf 1_N \\ (1-\tau) \mathbf 1_N \end{bmatrix},$$

garantindo assim que$c^\top z = \underbrace{\mathbf 0^\top(\theta^+ - \theta^-)}_{=0}+\tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v = \sum_{i=1}^N \rho_\tau(\epsilon_i).$

Portanto e são então definidos e o programa, conforme indicado em completamente especificado.$c,A$$b$$(1)$

Provavelmente é melhor digerido usando um exemplo. Para resolver isso em R, use o pacote quantreg de Roger Koenker. Aqui está também uma ilustração de como configurar o programa linear e resolver com um solucionador de programas lineares:

base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)
attach(base)
library(quantreg)
library(lpSolve)
tau <- 0.3


# Problem (1) only one covariate
X <- cbind(1,base$area)
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_euro
const_type <- rep("=",N)

linprog <- lp("min",c,A,const_type,b)
beta <- linprog$sol[1:K] -  linprog$sol[(1:K+K)]
beta
rq(rent_euro~area, tau=tau, data=base)


# Problem (2) with 2 covariates
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_euro
const_type <- rep("=",N)

linprog <- lp("min",c,A,const_type,b)
beta <- linprog$sol[1:K] -  linprog$sol[(1:K+K)]
beta
rq(rent_euro~ area + yearc, tau=tau, data=base)

Reescrevi o código de Jesper Hybel em Python usando cvxopt. Estou postando aqui, caso alguém também precise disso em Python.

import pandas as pd
import io
import requests
import numpy as np
url="http://freakonometrics.free.fr/rent98_00.txt"
s=requests.get(url).content
base=pd.read_csv(io.StringIO(s.decode('utf-8')), sep='\t')


tau = 0.3

from cvxopt import matrix, solvers

X = pd.DataFrame(columns=[0,1])
X[1] = base["area"] #data points for independent variable area
X[2] = base["yearc"] #data points for independent variable year
X[0] = 1 #intercept

K = X.shape[1]
N = X.shape[0]

# equality constraints - left hand side

A1 = X.to_numpy() # intercepts & data points - positive weights
A2 = X.to_numpy() * - 1 # intercept & data points - negative weights
A3 = np.identity(N) # error - positive
A4 = np.identity(N)*-1 # error - negative

A = np.concatenate((A1,A2,A3,A4 ), axis= 1) #all the equality constraints 

# equality constraints - right hand side
b = base["rent_euro"].to_numpy()

#goal function - intercept & data points have 0 weights
#positive error has tau weight, negative error has 1-tau weight
c = np.concatenate((np.repeat(0,2*K), tau*np.repeat(1,N), (1-tau)*np.repeat(1,N) ))

#converting from numpy types to cvxopt matrix

Am = matrix(A)
bm = matrix(b)
cm = matrix(c)

# all variables must be greater than zero
# adding inequality constraints - left hand side
n = Am.size[1]
G = matrix(0.0, (n,n))
G[::n+1] = -1.0

# adding inequality constraints - right hand side (all zeros)
h = matrix(0.0, (n,1))

#solving the model
sol = solvers.lp(cm,G,h,Am,bm, solver='glpk')

x = sol['x']

#both negative and positive components get values above zero, this gets fixed here
beta = x[0:K] - x[K :2*K]

print(beta)
```