A regressão quantílica é uma técnica estatística destinada a estimar e conduzir inferências sobre funções quantílicas condicionais. Assim como os métodos clássicos de regressão linear baseados na minimização de somas de resíduos quadrados permitem estimar modelos para funções médias condicionais, os métodos de regressão quantílica oferecem um mecanismo para estimar modelos para a função mediana condicional e toda a gama de outras funções quantílicas condicionais. Ao complementar a estimativa de funções médias condicionais com técnicas para estimar uma família inteira de funções quantílicas condicionais, a regressão quantílica é capaz de fornecer uma análise estatística mais completa das relações estocásticas entre variáveis aleatórias.

A regressão quantílica tem sido usada em uma ampla variedade de configurações de aplicação. As curvas de crescimento de referência para altura e peso das crianças têm uma longa história na medicina pediátrica; métodos de regressão quantílica podem ser usados para estimar curvas de referência de quantil superior e inferior em função da idade, sexo e outras covariáveis, sem impor suposições paramétricas rigorosas sobre as relações entre essas curvas. Os métodos de regressão quantílica têm sido amplamente utilizados na economia para estudar os determinantes dos salários, os efeitos da discriminação e as tendências da desigualdade de rendimentos.

Vários estudos recentes modelaram o desempenho dos alunos das escolas públicas em exames padronizados em função de características socioeconómicas, como o rendimento e o nível de escolaridade dos pais, e de variáveis políticas como o tamanho das turmas, as despesas escolares e as qualificações dos professores. Parece bastante implausível que todos esses efeitos covariáveis atuem de modo a alterar toda a distribuição dos resultados dos testes num valor fixo. É de interesse óbvio saber se as intervenções políticas alteram o desempenho dos alunos mais fortes da mesma forma que os alunos mais fracos são afectados. Tais questões são naturalmente investigadas dentro da estrutura de regressão quantílica.

Em ecologia, a teoria frequentemente sugere como as covariáveis observáveis afetam a limitação do tamanho sustentável da população, e a regressão quantílica tem sido usada para estimar diretamente modelos para os quantis superiores da distribuição condicional, em vez de inferir tais relações a partir de modelos baseados na tendência central condicional. Na análise de sobrevivência, e na análise do histórico de eventos em geral, há muitas vezes também um desejo de concentrar a atenção em segmentos específicos da distribuição condicional, por exemplo, as perspectivas de sobrevivência dos mais velhos, sem a imposição de pressupostos distribucionais globais.

A linguagem de programação R oferece vários recursos para implementar a regressão quantílica, nós utilizaremos neste texto o pacote quantreg. Neste pacote estão implementados modelos paramétricos, não paramétricos lineares e não lineares (penalizados pela variação total) para quantis condicionais de uma resposta univariada e vários métodos para lidar com dados de sobrevivência censurados. Os métodos de selecção de modelos com base no risco de quebra esperado também estão agora incluídos. Ver Koenker (2005) e Koenker, R. et al. (2017).



1. Introdução


Os quantis oferecem uma maneira conveniente de resumir distribuições de probabilidade univariadas, conforme exemplificado pelos onipresentes boxplots de Tukey. Em contraste com os momentos, que caracterizam características globais da distribuição e são, consequentemente, fortemente influenciados pelo comportamento da cauda, os quantis são inerentemente locais e quase impermeáveis a pequenas perturbações da massa distributiva.

Podemos mover a massa acima e abaixo do canteiro central sem perturbá-la, desde que, é claro, a massa não seja transferida de cima para baixo do canteiro central, ou vice-versa. Esta localidade dos quantis é altamente vantajosa pelas mesmas razões que as funções de base suportadas localmente são vantajosas na regressão não paramétrica, porque assegura uma forma de robustez que falta em muitos procedimentos estatísticos convencionais, nomeadamente aqueles baseados na minimização de somas de resíduos quadrados.

Quando existem covariáveis – e quase sempre existem covariáveis quando os problemas se tornam realmente interessantes – não podemos confiar na ordenação como estratégia para calcular quantis. Em vez disso, felizmente, existe uma alternativa de otimização simples e elegante.

Os quantis univariados surgem como soluções para o problema da perda esperada linear por partes, \[ \min_\xi \mbox{E}\big(\rho_\tau (Y-\xi) \big), \] onde \[ \rho_\tau (u) = \big(\tau-\pmb{I}(u < 0)\big)u \] e \(\tau\in (0,1)\). Explicar aqui a função acima, ver figura

Quando a distribuição de \(Y\) admite um único \(\tau\)-ésimo quantil podemos diferenciar, \[ \mbox{E}\big(\rho_\tau (Y-\xi) \big)=\tau\int_\xi^\infty (y-\xi)\mbox{d}F_Y(y)+(\tau-1)\int_{-\infty}^\xi (y-\xi)\mbox{d}F_Y(y) \] para obter a condição de primeira ordem, \[ 0=\tau\int_\xi^\infty \mbox{d}F_Y(y)+(\tau-1)\int_{-\infty}^\xi \mbox{d}F_Y(y)=F_Y(\alpha)-\tau, \] de maneira que \[ \alpha=F^{-1}_Y(\tau)\cdot \]

Quando existem múltiplos valores tais que \(F_Y(y)=\tau\), é convencional escolher o menor, ou seja, \[ \alpha = \inf\{y \, : \, F_Y(y)\geq \tau \}\cdot \]

Correspondendo a esses quantis populacionais estão expressões análogas para os quantis amostrais com \(F_Y\) substituído pela função de distribuição empírica \[ \widehat{F}_n(y) = \dfrac{1}{n} \sum_{i=1}^n \pmb{I}(Y_i\leq y)\cdot \]

A admissibilidade do quantil amostral univariado sob a perda \(\rho_\tau\) foi considerada por Fox and Rubin (1964), mas a origem de tais soluções sob perda linear assimétrica remonta pelo menos a Edgeworth (1888a).

Os estimadores de regressão que minimizam somas de resíduos absolutos também têm uma longa história. Já no século XVIII, Boscovich, e um pouco mais tarde Laplace, defenderam uma forma de regressão bivariada que restringia o resíduo médio a zero e minimizava a soma dos resíduos absolutos para encontrar a estimativa do parâmetro de inclinação restante.

Um século depois, Edgeworth (1888b) propôs remover a restrição de intercepto e determinar os parâmetros de inclinação e intercepto, minimizando a soma dos resíduos absolutos. Ele forneceu um algoritmo eficaz para calcular o estimador que antecipa os métodos simplex modernos. Para mais detalhes ver Koenker (2017). A proposta de Edgeworth definhou até ser reavivada na década de 1950, quando foi reconhecida como um programa linear.

Uma aplicação inicial da regressão mediana em economia aparece no trabalho de Arrow and Hoffenberg (1959), que a consideraram conveniente para estimar coeficientes de insumo-produto sujeitos a restrições de positividade. Embora houvesse um reconhecimento geral de que a abordagem da mediana ou \(\ell_1\) tinha a vantagem de ser mais resistente a valores discrepantes do que o estimador de mínimos quadrados usual, uma desvantagem da abordagem, além da falta de familiaridade com seus métodos computacionais, era a ausência de um aparato formal de inferência.

Não houveram quaisquer reconhecimento de que poderia ser interessante considerar modelos de regressão quantílica diferentes da mediana até que Gib Bassett e Roger Koenker começaram a explorar esse território em meados da década de 1970. Eles começaram com a observação de que, desde que o modelo “contivesse um intercepto”, ou seja, que a extensão linear da covariável/matriz de projeto \(X\), incluísse um vetor constante, então soluções para o análogo de regressão do nosso problema elementar, \[ \min_{\beta\in\mathbb{R}^p} \sum_{i=1}^n \rho_\tau \big(y_i-x_i^\top\beta\big), \] tinha a propriedade de que aproximadamente \(\tau n\) dos resíduos, \(r_i = y_i-x_i^\top\widehat{\beta}\), \(i=1,\cdots,n\) seriam positivos e \((1-\tau)n\) seriam negativos. Isto decorre imediatamente da condição do subgradiente que exige que no ótimo, \(\widehat{\beta}(\tau)\), \[ 0\in \partial_\beta \sum_{i=1}^n \rho_\tau \left.\big(y_i-x_i^\top\beta\big)\right|_{\beta=\widehat{\beta}}\cdot \]

Aqui, \[ \partial_\beta \rho_\tau\big(y_i-x_i^\top\beta\big)=-\psi_\tau\big(y_i-x_i^\top\beta\big)x_i, \] com \(\psi_\tau(u)=\tau-\pmb{I}(u<0)\) para \(u\neq0\) e é de valor definido, com \(\partial_\beta \rho_\tau\big(y_i-x_i^\top\beta\big)=[-\tau,1-\tau]x_i\) quando o resíduo é zero.

Quando as observações estão em “posição geral”, de modo que não mais do que \(p\) observações estão em um hiperplano de dimensão \(p\) no espaço amostral de regressão, a condição do subgradiente implica que \(\tau\) deve estar entre \(N/n\) e \((N+p)/n\) onde \(N\) é o número de observações abaixo do hiperplano ajustado, ou seja, com resíduos estritamente negativos.

Isso pareceu justificar a conjectura de que as soluções \(\widehat{\beta}(\tau)\) de tais problemas poderiam ser consideradas análogas aos quantis amostrais para o modelo linear, estimando os parâmetros de modelos que especificavam funções quantílicas condicionais afins para \(Y|X\).

A inferência começou a ser derivada para expressões combinatórias para a densidade finita amostral de \(\widehat{\beta}(\tau)\) com base na condição de otimalidade do gradiente anterior. Como isso envolvia um somatório de todas as \(\binom{n}{p}\) soluções de subconjuntos elementares que correspondiam a ajustes exatos de \(p\) observações, a princípio não parecia muito prático. Mas finalmente conseguiu-se mostrar que esta densidade tinha uma forma limitante gaussiana simples que justificava plenamente a terminologia da regressão quantílica que começava-se a usar.

No devido tempo, esses resultados apareceram em Koenker and Bassett (1978). Desde então, muitas pessoas contribuíram para um esforço que gradualmente construiu uma extensa caixa de ferramentas para estimativa e inferência sobre modelos quantis condicionais.


2. Inferência para modelos quantílicos condicionais


Um preceito fundamental das estatísticas é que as estimativas das magnitudes dos efeitos devem ser acompanhadas de alguma avaliação da precisão dessas estimativas. Nesta seção revisaremos uma variedade de métodos que evoluíram para resolver essa tarefa para a regressão quantílica.

O comportamento assintótico do processo de regressão quantílica \(\{\widehat{\beta}(\tau) \, : \, \tau \in (0,1)\}\) é muito paralelo à teoria dos quantis amostrais ordinários no problema de uma amostra. Koenker and Bassett (1978) mostram que no modelo linear clássico, \[ Y_i = x_i\beta +u_i, \] com \(u_i \sim F\), ou seja, os \(u_i\) são variáveis aleatórias independentes identicamente distribuídas segundo a função de distribuição \(F\), com densidade \(f(u)>0\) no suporte \(\{u \, : \, 0<F(u)<1\}\), a função de distribuição conjunta de \[ \sqrt{n}\big(\widehat{\beta}(\tau_i)-\beta(\tau_i)\big)_{i=1}^m \] é asintóticamente normal com média zero e matriz de covariâncias \(\Omega\otimes D^{-1}\), este chamado de produto tensorial entre \(\Omega\) e \(D^{-1}\).

Aqui, \(\beta(\tau)=\beta+F^{-1}_u(\tau)e_1\). \(e_1=(1,0,\cdots,0)^\top\), \(x_{1i}=1\), \(n^{-1}\sum x_ix_i^\top\to D\), uma matriz definida positiva e \[ \Omega = \big( \omega_{ij} = \min\{\tau_i,\tau_j\}-\tau_i\tau_j \big)/\big( f(F^{-1}(\tau_i))f(F^{-1}(\tau_j))\big)\cdot \]

Quando a resposta é condicionalmente independente sobre \(i\), mas não distribuída de forma idêntica, a matriz de covariância assintótica de \(\xi(\tau)=\sqrt{n}\big( \widehat{\beta}(\tau)-\beta(\tau)\big)\) é um pouco mais complicada. Seja \[ \xi_i(\tau) = x_i \beta(\tau) \] denotando a distribuição condicional da função quantílica de \(Y\) dado \(x_i\), \(f_i(\cdot)\) a correspondente função de densidade condicional, \[ J_n(\tau_i,\tau_j)=\big( \min(\tau_i,\tau_j)-\tau_i\tau_j\big)\frac{1}{n}\sum_{i=1}^n x_ix_i^\top \] e \[ H_n(\tau)=\frac{1}{n}\sum_{i=1}^n x_ix_i^\top f_i(\xi_i(\tau))\cdot \]

Sob condições de regularidade moderada nos \(\{f_i\}\)’s e \(\{x_i\}\)’s, temos normalidade assintótica conjunta para os vetores \(\big(\xi(\tau_1),\cdots,\xi(\tau_m)\big)\) com média zero e matriz de covariância \[ V_n = \Big( H_n(\tau_i)^{-1}J_n(\tau_i,\tau_j)H_n(\tau_j)^{-1}\Big)_{i,j=1}^m\cdot \]


2.1 Interpretação da regressão quantílica


A estimativa de mínimos quadrados de modelos de regressão média levanta a questão: como a média condicional de \(Y\) depende das covariáveis \(X\)? A regressão quantílica faz esta pergunta em cada quantil da distribuição condicional, permitindo obter uma descrição mais completa de como a distribuição condicional de \(Y\) dado \(X=x\) depende de \(x\).

Em vez de assumir que as covariáveis mudam apenas a localização ou escala da distribuição condicional, os métodos de regressão quantílica permitem explorar também os efeitos potenciais sobre a forma da distribuição. Assim, por exemplo, o efeito de um programa de formação profissional sobre a duração do actual período de desemprego dos participantes pode ser o de prolongar os períodos mais curtos e, ao mesmo tempo, reduzir drasticamente a probabilidade de períodos muito longos.

O efeito médio do tratamento em tais circunstâncias pode ser pequeno, mas o efeito do tratamento na forma da distribuição das durações do desemprego pode, no entanto, ser bastante significativo.

O termo “condicional” na regressão quantílica vem da especificação linear \[ Q_Y(\tau \, | \, X=x) = x^\top \beta(\tau)\cdot \]

Ao contrário da análise de mínimos quadrados, a interpretação dos coeficientes da regressão quantilica é o efeito marginal do regressor no quantil definido pelas covariáveis. O estimador de mínimos quadrados convencional modela a média condicional da resposta \(Y\), dadas as covariáveis \(X\) como \[ \mbox{E}(Y \, | \, X=x) = x^\top \beta, \] isto implica que \[ \mbox{E}(Y)=\mbox{E}\big( \mbox{E}(Y \, | \, X=x)\big) = \mbox{E}(X^\top\beta)=\mbox{E}(X)^\top\beta, \] isto é, os mínimos quadrados estimam consistentemente o efeito do aumento do valor médio de \(X\) na média incondicional de \(Y\), devido à lei das esperanças iteradas.

No caso do \(\tau\)-ésimo quantil, geralmente o modelo de regressão quantil condicional não produz o efeito no quantil incondicional como \[ Q_\tau (Y \, | \, X=x) = X^\top\beta(\tau), \] então \[ \mbox{E}\big(Q_\tau (Y \, | \, X=x)\big) = \mbox{E}\big(X^\top\beta(\tau)\big), \] logo \[ q_\tau = \mbox{E}(X^\top)\beta(\tau), \] onde \(q_\tau\) é o \(\tau\)-ésimo quantil incondicional de \(Y\).


2.2. Efeitos dos tratamentos quantílicos


A formulação mais simples de regressão quantílica é a de dois modelo de controle de tratamento de amostra. No lugar do modelo clássico de desenho experimental de Fisher, no qual o tratamento induz uma simples mudança de locação da distribuição da resposta, Lehmann (1974) propôs o seguinte modelo geral de resposta ao tratamento:

“Suuponha que o tratamento adicione a quantidade \(\Delta(x)\) quando a resposta do sujeito não tratado seria \(x\). Então a distribuição \(G\) das respostas ao tratamento é a da variável aleatória \(X+\Delta(X)\) onde \(X\) é distribuído de acordo com \(F\).”

Os casos especiais obviamente incluem o modelo de mudança de locação \(\Delta(X) = \Delta_0\) e o modelo de mudança de escala \(\Delta(X) = \Delta_0 X\), mas o caso geral é natural dentro do paradigma de regressão quantílica. Doksum (1974) mostra que se \(\Delta(x)\) for definido como a “distância horizontal” entre \(F\) e \(G\) em \(x\), então \[ F(x)=G(x+\Delta(x)), \] então \(\Delta(x)\) é definido univocamente e pode ser expresso como \[ \Delta(x)=G^{-1}(F(x))-x\cdot \]

Mudando as variáveis para \(\tau=F(x)\) pode-se definir o efeito do tratamento quantílico, \[ \delta(\tau)=\Delta\big(F^{-1}(\tau)\big)=G^{-1}(\tau)-F^{-1}(\tau)\cdot \] No cenário de duas amostras, esta quantidade é naturalmente estimável por \[ \widehat{\delta}(\tau)=\widehat{G}^{-1}_n(\tau)-\widehat{F}^{-1}_m(\tau), \] onde \(\widehat{G}^{-1}_n\) e \(\widehat{F}^{-1}_m\) denotam as funções de distribuição empírica das observações de tratamento e controle, com base em \(n\) e \(m\) observações, respectivamente.

Formulando o modelo de regressão quantílica para o problema de tratamento binário como, \[ Q_{Y_i}(\tau \, | \, D_i)=\alpha(\tau)+\delta(\tau) D_i, \] onde \(D_i\) denota o indicador de tratamento, com \(D_i = 1\) indicando tratamento, \(D_i = 0\), controle, então o efeito quantil do tratamento pode ser estimado resolvendo, \[ \big(\widehat{\alpha}(\tau)+\widehat{\delta}(\tau)\big)^\top = \mbox{argmin} \sum_{i=1}^n \rho_\tau(y_i-\alpha-\delta D_i)\cdot \]

A solução \(\big(\widehat{\alpha}(\tau)+\widehat{\delta}(\tau)\big)^\top\) produz \(\widehat{\alpha}(\tau) = \widehat{F}_n^{-1}(\tau)\), correspondente à amostra de controle, e \[ \widehat{\delta}(\tau)=\widehat{G}^{-1}_n(\tau)-\widehat{F}^{-1}_m(\tau)\cdot \]

Doksum sugere que se pode interpretar os sujeitos de controle em termos de uma característica latente: por exemplo, em aplicações de análise de sobrevivência, um sujeito de controle pode ser chamado de frágil se for propenso a morrer em uma idade precoce, e robusto se for propenso a morrer em uma idade avançada.

Esta característica latente é, portanto, indexada implicitamente por \(\tau\), o quantil da distribuição de sobrevivência em que o sujeito apareceria se não fosse tratado, ou seja, \((Y_i \, | \, D_i = 0) =\alpha(\tau)\) e o tratamento, no modelo de Lehmann, é assumido como alterar a resposta de controle do sujeito, \(\alpha(\tau)\), tornando-o \(\alpha(\tau) +\delta(\tau)\) sob o tratamento.

Se a característica latente, digamos, a propensão para a longevidade, fosse observável então poderíamos ver o efeito do tratamento \(\delta(\tau)\) como uma interação explícita com esta variável observável. Na ausência de tal variável observável, no entanto, o efeito quantílico do tratamento pode ser considerado como uma medida natural da resposta ao tratamento.

Pode-se notar que o efeito do tratamento quantílico está intimamente ligado ao gráfico QQ-plot de duas amostras o qual tem uma longa história como dispositivo gráfico de diagnóstico. A função \[ \Delta(x)=\widehat{G}^{-1}_n\big(\widehat{F}_m(x)\big)-x \] é exatamente o que é plotado no gráfico QQ-plot tradicional de duas amostras. Se \(F\) e \(G\) forem idênticos então a função \(\widehat{G}^{-1}_n\big(\widehat{F}_m(x)\big)\) ficará ao longo da linha de 45 graus; se eles diferirem apenas por uma mudança de escala de locação, então \(\widehat{G}^{-1}_n\big(\widehat{F}_m(x)\big)\) ficará ao longo de outra linha com intercepto e inclinação determinadas pela locação e mudança de escala, respectivamente. A regressão quantílica pode ser vista como um meio de estender o gráfico QQ-plot de duas amostras e métodos relacionados para configurações gerais de regressão com covariáveis contínuas.

Quando a variável de tratamento assume mais de dois valores, o efeito do tratamento quantílico de Lehmann-Doksum requer apenas uma pequena reinterpretação. Se a variável de tratamento for contínua como, por exemplo, em estudos de dose-resposta, então é natural considerar a suposição de que seu efeito é linear, e escrever, \[ Q_{Y_i}(\tau \, | \, x_i)=\alpha(\tau)+\beta(\tau)x_i\cdot \] Assumimos assim que o efeito do tratamento \(\beta(\tau)\), de mudar \(x\) de \(x_0\) para \(x_0+1\) é o mesmo que o efeito do tratamento de uma alteração de \(x\) de \(x_1\) para \(x_1+1\): Observe que esta noção de tratamento quantil o efeito mede, para cada \(\tau\), a mudança na resposta necessária para permanecer na \(r\)-ésima função quantil condicional.


2.2.1. Efeitos binários do tratamento


A configuração da regressão quantílica mais simples é a resposta binária ao tratamento ou modelo de duas amostras, onde temos um indicador de tratamento \(D_i\), que assume o valor 1 para observações “tratadas” e 0 para observações “controle”. Na versão clássica do tratamento da média o interesse concentra-se exclusivamente na diferença nas médias das duas amostras, \[ \mbox{E}(Y_i \, | \, D_i) = \alpha+\beta D_i\cdot \]

Isto é comummente justificado pelo modelo de mudança de locação expresso como, \[ Y_i=\alpha+\beta D_i+u_i, \] onde os \(u_i\) são assumidos implícita ou explicitamente como independentes e distribuídos de forma idêntica. Assim, efetivamente, acredita-se que o tratamento desloque toda a distribuição da resposta em sincronia pela quantidade \(\beta\).

Em contraste, o modelo de efeito de tratamento quantílico, \[\begin{equation*} \tag{1} Q_{Y\, | \, D}(\tau \, | \, D)=\alpha(\tau)+\beta(\tau)D \end{equation*}\] reconhece que a distribuição da resposta pode ser arbitrariamente diferente no âmbito dos regimes de tratamento e controle.

Nesta formulação \(\alpha(\tau)\) denota a função quantílica da resposta para controles, \(F^{-1}_{Y\, | \, |D}(\tau)\) e \(\beta(\tau)\) denotam a diferença entre as funções quantílicas da resposta de tratamento e controle: \[ F^{−1}_{Y \, | \, D=1}(\tau) - F^{−1}_{Y \, | \, D=0}(\tau)\cdot \]

Este modelo está intimamente relacionado com a proposta de Lehmann (1974) de generalizar o modelo de efeito de tratamento médio considerando a diferença horizontal entre as funções de distribuição de tratamento e controle, que ele definiu como a função \(\Delta(y)\) tal que, \[ F_{Y \, | \, D=1}(y) - F_{Y \, | \, D=0}(y+\Delta(y))\cdot \] Assim, o efeito de tratamento da média escalar torna-se um objeto funcional capaz de descrever completamente a diferença entre as distribuições de tratamento e controle.

Até este ponto, mantivemos silêncio sobre possíveis problemas de correlações ou seleções que possam surgir em relação à atribuição do tratamento. Uma vez que somos quase inevitavelmente incapazes de observar tanto a resposta ao tratamento como a resposta ao controle para indivíduos individuais, somos consequentemente constrangidos a inferir sobre distribuições marginais.

A presunção subjacente ao modelo é que um sujeito de controle com resposta \(\alpha(\tau)\) terá, se tratado, resposta \(\alpha(\tau)+\beta(\tau)\). Esta presunção é por vezes referida como “invariância de classificação”, como por exemplo, em Heckman, Smith and Clements (1997).

Desde que o tratamento seja atribuído aleatoriamente, a estimativa do modelo (1) pode ser implementada. Dada uma amostra \(\{(y_i, d_i) \, : \, i = 1,\cdots,n\}\) podemos simplesmente resolver, \[ \big(\widehat{\alpha}(\tau),\widehat{\beta}(\tau)\big) = \mbox{argmin}_{(\alpha,\beta)} \sum_{i=1}^n \rho_\tau \big(y_i-\alpha+\beta d_i \big)\cdot \]

Para isso nem precisamos de nenhum maquinário de programação linear, pois o problema se separa em dois problemas distintos com soluções dadas pelo quantil amostral comum para as amostras de controle e tratamento. Ver Koenker (2005) para mais detalhes.

Como consequência, a inferência sobre o modelo (1) de tratamento binário também pode se basear na teoria clássica de grandes amostras para os quantis amostrais comuns. Como as duas amostras são independentes, temos que \(\widehat{\beta}(\tau)\) tem distribuições assintóticas de dimensão finita, \[ \sqrt{n}\big(\widehat{\beta}(\tau)-\beta(\tau)\big) \sim N\big(0,\lambda_0\Omega^0 +\lambda_1\Omega^1 \big), \] onde \[ \Omega_{ij}^k = (\tau_i-\tau_j\wedge \tau_j)/\Big(f_k\big(F_k^{-1}(\tau_i)\big)f_k\big(F_k^{-1}(\tau_j)\big) \Big), \] para \(k=1,2\) e \(i,j=1,\cdots,p\) e \(\lambda_k=n/n_k\) desde que os tamanhos amostrais relativos \(n_k/n\) permaneçam limitados a zero e um quando \(n\to\infty\). Claro, isto levanta a questão de como estimar as matrizes \(\Omega^k\), uma vez que envolvem as funções de densidade condicional das duas amostras. Isso gerou uma literatura bastante extensa, e uma breve visão geral é fornecido em Koenker (2005).

A teoria anterior nos permite construir faixas de confiança para as estimativas do modelo (1) usando a matriz de covariância estimada. Bandas uniformes representam um desafio um pouco maior; uma abordagem atraente seria empregar a versão assintótica da abordagem de Hotelling (1939) descrita em Koenker (2011). Várias outras abordagens de reamostragem também foram sugeridas recentemente, nomeadamente por Belloni, Chernozhukov and Kato (2016) e Hagemann (2016).

Dada a ênfase tradicional colocada em modelos de mudança de locazação de resposta ao tratamento, por ex. Cox (1984), é de algum interesse explorar testes deste modelo hipotético. Esses testes estão intimamente relacionados aos testes clássicos de qualidade de ajuste envolvendo parâmetros estimados. Uma abordagem para tais testes, seguindo Khmaladze (1981), é descrita em Koenker and Xiao (2002).


2.2.2. Múltiplos tratamentos, covariáveis e interações concomitantes


Expandir o paradigma de tratamento binário para permitir múltiplas opções de tratamento levanta algumas questões novas, especialmente do ponto de vista dos testes, no entanto, o modelo (1) ainda pode ser baseado em diferenças de quantis amostraiss univariados e, portanto, as regiões de confiança podem ser baseadas em teoria essencialmente semelhante à já descrita.

Um problema tentador e de importância crescente em muitos campos é o da atribuição de tratamento: dado um modelo estimado dos efeitos do tratamento, como deveríamos proceder para atribuir novos sujeitos a vários regimes de tratamento? Tais questões, especialmente na área médica, exigem respostas a questões espinhosas de avaliação de risco, onde uma perspectiva de distribuição dos efeitos heterogêneos do tratamento pode ser crucial. Uma nova perspectiva sobre estas questões é oferecida no trabalho recente de Wang, Zhou, Song and Sherwood (2016), baseado parcialmente em Manski (2004).

Quando há covariáveis concomitantes além das variáveis indicadoras de tratamento, surgem mais novas questões. Se a atribuição do tratamento for totalmente aleatória, é tentador simplesmente ignorar estas covariáveis; este é o ponto de vista articulado por Freedman (2008), que argumenta que o preconceito induzido pela introdução mal especificada de efeitos de covariáveis estranhos é provavelmente mais prejudicial do que quaisquer benefícios que possam advir da redução da variância.

Este argumento tem pelo menos força igual na configuração da regressão quantílica como na regressão média. Presumivelmente, a aleatorização deixa-nos com o tratamento \(D\) que é estocasticamente independente de outras covariáveis, digamos \(X\), pelo que o condicionamento adicional em \(X\) não ajudará e poderá dificultar quando a forma funcional do condicionamento \(X\) for mal escolhida. É claro que quando o tratamento é atribuído “aos observáveis” \(X\) então o argumento para a sua inclusão é muito mais convincente. Kadane and Seidenfeld (1996) fornecem uma boa discussão sobre isso à luz da infame crítica de Student (1931) ao experimento do leite em Lanarkshire. Em vez de confiar em tais argumentos de seleção em observáveis, vários autores optaram pela reponderação do escore de propensão.

Um exemplo inicial disto é o trabalho de Lipsitz, Fitzmaurice, Molenberghs and Zhao (1997) com contribuições posteriores de Firpo (2007) e outros. O extenso trabalho recente sobre os chamados métodos “duplamente robustos” que combinam estas abordagens também poderia ser empregado, conforme sugerido recentemente por Diaz (2016).

Um tanto negligenciado na literatura econométrica sobre resposta ao tratamento e avaliação de programas é o papel potencialmente importante das interações de covariáveis com variáveis de tratamento. Embora as interações apareçam com destaque na literatura clássica de análise de variância e apareçam em algumas pesquisas recentes de modelos lineares de alta dimensão, a econometria tende a concentrar a atenção nos principais efeitos do tratamento.

As interações, se presentes, devem desempenhar um papel essencial na atribuição do tratamento pós-análise. Mais trabalhos precisam ser feitos para desenvolver melhores ferramentas de diagnóstico para incorporar tais efeitos. Cox (1984) oferece uma extensa agenda de tópicos de pesquisa abertos, muitos dos quais poderiam ser estendidos de forma proveitosa ao cenário de regressão quantílica.


2.3. Equivariância por transformações da regressão quantílica


Uma propriedade importante do modelo de regressão quantílica é que, para qualquer função monótona, \(h(\cdot)\), \[ Q_{h(T)}(\tau \, | \, x) = h\big(Q_{T}(\tau \, | \, x) \big)\cdot \]

Isto decorre imediatamente da observação de que \[ P(T<t \, | \, x) = P\big( h(T)<h(t) \, | \, x\big)\cdot \] Esta equivariância às transformações monótonas da função quantílica condicional é uma característica crucial, permitindo dissociar os objetivos potencialmente conflitantes das transformações da variável resposta.

Esta propriedade de equivariância contrasta diretamente com os conflitos inerentes à estimativa de modelos de transformações para relações médias condicionais. Como, em geral, \[ \mbox{E}(h(T) \, | \, x) \neq h\big(\mbox{E}(T\, | \, x)\big), \] a transformação altera de forma fundamental o que está sendo estimado na regressão de mínimos quadrados ordinária.

Uma aplicação particularmente importante deste resultado de equivariância, e que se revelou extremamente influente na aplicação econométrica da regressão quantílica, envolve a censura da variável de resposta observada. O modelo mais simples de censura pode ser formulado da seguinte forma. Seja \(Y_i^*\) denotar uma resposta latente, ou seja, não observável que se presume ser gerada a partir do modelo linear \[ Y_i^* = x_i^\top\beta+u_i, \qquad i=1,\cdots,n \] com \(\{u_i\}\) iid da função de distribuição \(F\). Devido à censura, os \(Y_i^*\) não são observados diretamente, mas sim observados \[ Y_i=\max\{0,Y_i^*\}\cdot \]

Powell (1986) observou que a equivariância dos quantis às transformações monótonas implicava que neste modelo as funções quantílicas condicionais da resposta dependiam apenas do ponto de censura, mas eram independentes de \(F\).

Formalmente, a \(r\)-ésima função quantílica condicional da resposta observada \(Y_i\); neste modelo pode ser expresso como \[ Q_i(\tau \, | \, x_i) = \max\{0,x_i^\top\beta+F^{-1}_u(\tau)\}\cdot \]

Os parâmetros das funções quantílicas condicionais podem agora ser estimados resolvendo \[ \min_{\beta}\sum_{i=1}^n \rho_\tau \big( y_i-\max\{0,x_i^\top\beta\}\big) \] onde se assume que os vetores \(x_i\) contêm um intercepto para absorver o efeito aditivo de \(F_u^{-1}(\tau)\). Este modelo é computacionalmente um pouco mais exigente do que a regressão quantílica linear convencional porque é não linear em parâmetros.


2.4. Robustez


A robustez das suposições distribucionais é uma consideração importante em todas as estatísticas, por isso é importante enfatizar que a regressão quantílica herda certas propriedades de robustez dos quantis comuns amostrais. As estimativas e os testes inferencias associados têm um caráter inerentemente livres de distribuição ou não paramétricos, uma vez que as estimativa dos quantis são influenciadas apenas pelo comportamento local da distribuição condicional da resposta perto do quantil especificado.

Dada uma solução \(\widehat{\beta}(\tau)\), baseada nas observações \(\{y,X\}\), desde que não se altere o sinal dos resíduos, qualquer uma das \(y\) observações pode ser alterada arbitrariamente sem alterar a solução inicial. Apenas os sinais dos resíduos são importantes na determinação das estimativas de regressão quantílica e, portanto, as respostas periféricas influenciam o ajuste na medida em que estão acima ou abaixo do hiperplano estimado, mas quão acima ou abaixo é irrelevante.

Embora as estimativas da regressão quantílica sejam inerentemente robustas à contaminação das observações da resposta, elas podem ser bastante sensíveis à contaminação das observações do projeto \(\{x_i\}\). Várias propostas foram feitas para melhorar este efeito.


2.5. Aspectos computacionais


Embora tenha sido reconhecido por vários autores iniciais, incluindo Gauss, que as soluções para o problema de regressão mediana eram caracterizadas por um ajuste exato através de \(p\) observações amostrais quando \(p\) parâmetros lineares são estimados, nenhum algoritmo eficaz surgiu até o desenvolvimento da programação linear no década de 1940.

Foi então rapidamente reconhecido que o problema da regressão mediana poderia ser formulado como um programa linear e o método simplex empregado para resolvê-lo. O algoritmo de Barrodale and Roberts (1973) forneceu a primeira implementação eficiente projetada especificamente para a regressão mediana e ainda é amplamente utilizado em softwares estatísticos. O pacote de funções R quantreg implementa este algoritmo, para isso define-se o modelo a ser estimado com a função rq e a opção method = “br”, utiliza o algoritmo mencionado como o padrão de ajuste.

Pode ser descrito concisamente da seguinte forma. Em cada etapa, temos um conjunto experimental de \(p\) “observações básicas” cujo ajuste exato pode constituir uma solução. Calculamos a derivada direcional da função objetivo em cada uma das \(2p\) direções que correspondem à remoção de uma das observações básicas atuais e à execução de um passo positivo ou negativo. Se nenhuma dessas derivadas direcionais for negativa a solução foi encontrada, caso contrário escolhe-se a mais negativa, a direção de descida mais acentuada, e segue nessa direção até que a função objetivo pare de diminuir.

Esta busca unidimensional pode ser formulada como um problema de encontrar a solução para um problema de quantil escalar ponderado. Tendo escolhido o comprimento do passo, determinamos de fato uma nova observação para entrar no conjunto básico, um pivô simplex ocorre para atualizar a solução atual e a iteração continua. Esta estratégia simplex modificada é altamente eficaz em problemas com um número modesto de observações, alcançando velocidades comparáveis às soluções de mínimos quadrados correspondentes. Mas para problemas maiores com, digamos, \(n > 100.000\) observações, a abordagem simplex eventualmente se torna consideravelmente mais lenta que os mínimos quadrados.

Para problemas grandes, o desenvolvimento recente de métodos de pontos interiores para problemas de programação linear é altamente eficaz. Portnoy and Koenker (1997) descrevem uma abordagem que combina algum pré-processamento estatístico com métodos de pontos interiores e alcança resultados comparáveis ao desempenho para soluções de mínimos quadrados, mesmo em problemas muito grandes. Este algoritmo, também implementado na função rq é especificado como àquele a ser utilizado como ajuste na opção method = “fn”.

Para problemas ainda maiores pode-se usar a abordagem Frisch-Newton após o pré-processamento utilizando a opção method = “pfn”. Também descrito detalhadamente em Portnoy and Koenker (1997), este método é principalmente adequado para problemas com \(n\) grande e \(p\) pequeno, ou seja, quando a dimensão paramétrica do modelo é modesta. Para problemas grandes com grandes dimensões paramétricas muitas vezes é vantajoso usar method = “sfn”, que também usa o algoritmo de Frisch-Newton, mas explora álgebra esparsa para calcular iterações. Isto é especialmente útil quando o modelo inclui variáveis fatores que, quando expandidas, geram matrizes de projeto ou planejamento muito esparsas. Atualmente, as opções de inferência, ou seja, métodos de summary (resumo), são um tanto limitadas quando se utiliza este método. Apenas a opção se = “nid” está disponível atualmente, mas espera-se a implementação de algumas opções bootstrap. Ainda temos outras possibilidades de ajustes.

Uma característica importante da formulação de programação linear de regressão quantílica é que toda a gama de soluções para \(\tau\in \{0,1\}\) pode ser calculada de forma eficiente por programação paramétrica. Em qualquer solução \(\widehat{\beta}(\tau_0)\) há um intervalo de \(\tau\)’s sobre o qual esta solução permanece ótima, podem-se calcular os pontos finais deste intervalo e, portanto, resolver iterativamente para todo o caminho da amostra \(\widehat{\beta}(\tau)\).

Mais dados e modelos mais complexos colocaram maior pressão sobre os recursos computacionais em todas as estatísticas e levaram a muitas inovações, incluindo métodos computacionais para regressão quantílica. Os primeiros métodos de programação linear baseados em simplex deram gradualmente lugar a métodos de pontos interiores, conforme descrito em Portnoy and Koenker (1997), e a álgebra esparsa expandiu ainda mais o escopo desses métodos, permitindo que modelos com vários milhares de parâmetros sejam estimados de forma eficiente, conforme ilustrado em Koenker (2011).

Mas as novas demandas excedem as capacidades até mesmo desses métodos e o advento da computação distribuída mudou a atenção para os métodos de gradiente descendente. Koenker (2017) descreve brevemente uma abordagem para tais estratégias de regressão quantílica com base na abordagem do método de multiplicadores de direção alternada (ADMM) de Parikh and Boyd (2014).


3. Regressão quantílica linear


A forma mais simples e irrestrita de estimativas nos modelos de regressão quantílica permite que as variáveis preditoras \(X\) exerçam mudanças na tendência central, na variância e na forma da distribuição da variável resposta \(Y\) (Koenker and Machado, 1999).

Quando o único efeito estimado é uma mudança na tendência central, por exemplo, na média da distribuição de \(Y\) condicionada aos valores de \(X\), temos o familiar modelo de regressão de variância homogênea associado à regressão de mínimos quadrados ordinários, ver a linha tracejada na Figura 1.

Todas as estimativas da inclinação da regressão quantílica \(\beta_1(\tau)\) são para um parâmetro comum, e qualquer desvio entre as estimativas da regressão quantílica é simplesmente devido à variação da amostragem.

n =90; set.seed(867)
# Gerar uma distribuição lognormal
erro <- rlnorm(90, meanlog = 0, sdlog = 0.75)
# Modelo
beta0 = 6.0; beta1 = 0.05; X = seq(0,100,length.out = 90); Y = beta0 + beta1*X + erro
par(mar=c(4,4,1,1))
plot(X,Y, pch = 19, ylim = c(6,14.5))
grid()
taus = c(0.90,0.75,0.50,0.25,0.10)
ajuste = quantreg::rq(Y ~ X, tau = taus)
for (i in 1:dim(fitted.values(ajuste))[2]){ lines(X,fitted.values(ajuste)[,i], col = i, lwd = 2)}
legend("bottom", legend = ajuste$tau, title = expression(tau), 
       horiz = TRUE, col = 1:dim(fitted.values(ajuste))[2], lwd = 2)
ajuste.lm = lm(Y ~X)
lines(X,fitted.values(ajuste.lm), lty = 2, lwd = 2)

Figura 1. Uma amostra (\(n\) = 90) de um modelo de erro homogêneo lognormal com mediana = 0 e \(\sigma\) = 0.75, \(Y = \beta_0 + \beta_1X_1 + \epsilon\), \(\beta_0 = 6.0\) e \(\beta_1 = 0.05\) com estimativas de regressão quantílica nos quantis 0.90, 0.75, 0.50, 0.25 e 0.10 (linhas sólidas) e estimativa de mínimos quadrados da função média (linha tracejada).

Uma estimativa da taxa de variação nas médias da regressão de mínimos quadrados ordinários também é uma estimativa do mesmo parâmetro na regressão quantílica. As estimativas dos interceptos \(\beta_0(\tau)\) do modelo de regressão quantílica são para o quantil paramétrico de \(Y\) quando \(X_1, X_2,\cdots, X_p = 0\), que diferem entre quantis e para a média.

Aqui, a principal virtude das estimativas do intercepto da regressão quantílica é que elas não dependem de uma forma assumida da distribuição do erro, como quando a regressão de mínimos quadrados é usada, e que assume uma distribuição de erros normal.

Abaixo apresentamos um resumo das estimativas dos parâmetros da regressão quantílica, ou seja, das estimativas de \(\beta_0(\tau)\) e \(\beta_1(\tau)\) para os valores dos quantis especificados na opção tau, no comando rq mostrado acima.

(ajuste)
## Call:
## quantreg::rq(formula = Y ~ X, tau = taus)
## 
## Coefficients:
##              tau= 0.10  tau= 0.25  tau= 0.50  tau= 0.75 tau= 0.90
## (Intercept) 6.35638305 6.81498144 7.15188364 7.90678880 8.5541826
## X           0.05086929 0.04737411 0.04567329 0.04215357 0.0574585
## 
## Degrees of freedom: 90 total; 88 residual


Por convenção, para todas as rotinas de ajuste do modelo linear no R, vemos apenas os coeficientes estimados e algumas informações sobre o modelo que está sendo estimado. Para obter uma avaliação mais detalhada do modelo ajustado, podemos usar summary.

A tabela resultante, mostrada a continuação, fornece o intercepto e a inclinação estimados na primeira coluna e os intervalos de confiança para esses parâmetros na segunda e terceira colunas. Por padrão, esses intervalos de confiança são calculados pelo método de inversão de classificação, descrito em Koenker (2005).

summary(ajuste)
## 
## Call: quantreg::rq(formula = Y ~ X, tau = taus)
## 
## tau: [1] 0.1
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 6.35638      6.25931  6.54433 
## X           0.05087      0.04795  0.05313 
## 
## Call: quantreg::rq(formula = Y ~ X, tau = taus)
## 
## tau: [1] 0.25
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 6.81498      6.53916  7.01713 
## X           0.04737      0.04392  0.05139 
## 
## Call: quantreg::rq(formula = Y ~ X, tau = taus)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 7.15188      6.97829  7.38539 
## X           0.04567      0.04279  0.04884 
## 
## Call: quantreg::rq(formula = Y ~ X, tau = taus)
## 
## tau: [1] 0.75
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 7.90679      7.17099  8.41477 
## X           0.04215      0.03556  0.05807 
## 
## Call: quantreg::rq(formula = Y ~ X, tau = taus)
## 
## tau: [1] 0.9
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept)  8.55418      8.01383 10.41035
## X            0.05746      0.01958  0.06677


No código seguinte traçamos o gráfico resumo do ajuste da regressão quantílica, Figura 2. Esta é uma ferramenta gráfica para avaliar a relação entre a variável de resposta \(Y\) e a variável preditora \(X\). A relação entre a variável resposta e as peditoras será avaliada de maneira indireta. Neste código informamos os diferentes quantis nos quais vamos avaliar a distribuição estimada dos coeficientes da regressão, especificados como tau = 5:45/50 ou

tau = 5:45/50; tau
##  [1] 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38
## [16] 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68
## [31] 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 0.86 0.88 0.90


No painel direito abaixo traçamos estimativas da densidade correspondente à \(\widehat{\beta}_0(\tau)\) ou intercepto do modelo. As estimativas da densidade empregam o método de kernel adaptativo proposto por Silverman (1986) e implementado na função quantreg akj. Esta função é particularmente conveniente porque permite associar massas desiguais a observações como as produzidas pelo processo de regressão quantílica.

No paines esquerdo temos um gráfico similar mas correspondente à \(\widehat{\beta}_1(\tau)\), rotulado com \(X\), o nome da variável preditora. O eixo \(x\) em ambos gráficos representa valores de \(\tau\) e o eixo \(y\) representa a densidade marginal da dendidade dos coeficientes estimados. As linhas vermelhas sólidas e tracejadas mostram o modelo linear e seus intervalos de confiança dos limites inferior e superior, respectivamente.

No gráfico a esquerda podemos ver que apenas uns poucos valores de \(\widehat{\beta}_1(\tau)\) se afastam do intervalo de confiança da inclinação do modelo linear, isto significa não haver diferença entre a inclinação do modelo linear e a inclinação na regressão quantílica.

Diferente acontece no gráfico a direita. Percebemos que entre os quantis 0.4 e quase 0.75 os valores médios de \(\widehat{\beta}_0(\tau)\) e boa parte de sua banda de confinaça estão dentro da banda de confiança do intercepto do modelo linear. Os valores de \(\widehat{\beta}_0(\tau)\) para \(\tau\) entremos, ou seja, menores do que 0.4 e miores do que 0.75 estão fora da banda de confiança.

Destes dois gráficos podemos afirmar que o modelo de regressão quantílico mediano, obtido quando \(\tau=0.5\) e o modelo linear não apresentam diferenças estatisticamente significativas. Ainda podemos afirmar que o coeficiênte estimado pelo modelo linear para a variável \(X\), que é constante, não difere do respectivo coeficiente estimado na regressão quantílica em quase todos os quantis.

Para o intercepto observamos que quando o quantil é menor do que 0.4 existe divergência estatística entre o intercepto da regressão quantílica e o intercepto da regressão lineaar. Também existe esta diferença quando os quantis são maiores do que 0.75 e ainda observamos assimetria na distribuição marginal de ambos coficientes nesta situação, isto sugere assimetria dos dados.

par(mfrow=c(1,2))
plot(summary(quantreg::rq(Y ~ X, tau = 5:45/50)), mfrow = c(1,2), mar = c(4,4,3,1))

Figura 2. As estimativas amostrais de \(\beta_0(\tau)\) (esquerda) e \(\beta_1(\tau)\) (direita) são mostradas como funções degrau discretas em pontos pretos. Linhas vermelhas tracejadas conectam pontos finais com intervalos de confiança de 90%.

Esta é uma forma particularmente comum de heterocedasticidade. Se olharmos com mais cuidado para o ajuste, veremos desvios interessantes da simetria que provavelmente não seriam revelados pelos testes típicos de heterocedasticidade nos livros didáticos. A solução para sintomas como estes seria reformular o modelo em termos logarítmicos. É interessante comparar o que acontece após a transformação do log com o que já vimos.

Observemos a Figura 3 abaixo. Parece que o modelo em torno do ponto mediano e quantis superiores se ajustam melhor. Em geral, devemos escolher o modelo com o valores de AIC mais baixa. AIC (Akaike Information Criterion) é uma métrica usada para avaliar a adequação de um modelo estatístico para um determinado conjunto de dados. Baseia-se na complexidade do modelo, bem como na sua qualidade de ajuste. O modelo é considerado melhor quanto menor o valor do AIC.

aic_df = data.frame(AIC = AIC(ajuste), model = paste("tau =",taus))
library(ggplot2)
ggplot(aic_df, aes(x = model, y = AIC)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = expression(paste("Valores de AIC para diferentes valores de ",tau)),
       x = expression(paste("Modelos de Regressão Quantílica com diferentes valroes de ",tau)),
       y = "Akaike Information Criterion (AIC)") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

Figura 3. Valores do Akaike Information Criterion (AIC) para os diferentes \(\tau\)’s.

Porém, é fundamental lembrar que o AIC não deve ser o único fator considerado na escolha de um modelo. Ao escolher um modelo, é importante levar em conta uma série de elementos adicionais, incluindo interpretabilidade, aplicabilidade e relevância teórica. Também é possível que nenhum dos modelos tenha um AIC visivelmente superior aos demais, caso em que outros fatores podem ser mais cruciais na seleção de um modelo.

Concluimos então que os melhores ajustes são obtidos para valores de \(\tau\) iguais ou maiores do que 0.5. Isto condiz com a assimetria à direita da densidade do erro, a distrobuição log-normal.


4. Regressão quantílica não linear


De vez em quando podemos nos deparar com especificações de modelos quantílicos condicionais com parâmetros não lineares, \[ Q_{Y_i}(\tau \, | \, x) = g\big(x,\theta(\tau)\big) \] que pode ser estimado da maneira imediatamente óbvia, \[ \widehat{\theta}(\tau)=\mbox{argmin}_\theta \sum_{i=1}^n \rho_\tau \big(y_i-g(x_i\theta)\big)\cdot \]

Os principais exemplos de tais circunstâncias são o estimador de Powell (1986) do modelo de regressão censurado (Tobit) e o estimador de escore máxima de Manski (1975) do modelo de resposta binária. Em ambos os casos, temos um modelo de resposta latente linear em parâmetros que postula, \[ Q_{Y_i^*}(\tau \, | \, x) = x^\top\beta, \] mas a resposta observável \(Y_i^* = \max\{0, Y_i^∗\}\) no primeiro caso e \(Y_i = \pmb{I}(Y_i^∗ > 0)\) no último.

Dado que \(Q_{h(Y)}(\tau) = h\big((Q_Y(\tau)\big)\) para qualquer transformação monótona \(h\), veja, por exemplo. Koenker (2005), segue-se que o estimador de Powell, \[ \widehat{\beta}(\tau)=\mbox{argmin}_\beta \sum \rho_\tau \big( y_i-\max\{0,x_i^\top\beta\}\big) \] e o estimador de Manski, \[ \widehat{\beta}(\tau)=\mbox{argmin}_{||\beta||=1} \sum \rho_\tau \big( y_i-\pmb{I}(x_i^\top\beta>0)\big) \] estima consistentemente os parâmetros do modelo de variável latente, até a escala neste último caso.

Outros modelos de transformação paramétrica oferecem outros exemplos. A venerável família de transformações de poder Box-Cox, adaptada ao cenário de regressão quantílica, afirma que, \[ Q_{h(Y,\lambda)}(\tau \, | \, X)=x_i^\top \beta(\tau) \] onde \(h(y,\lambda) = (y^\lambda-1)/\lambda\). É claro que se \(\lambda\) for conhecido podemos facilmente estimar os parâmetros lineares \(\beta(\tau)\), porém a estimativa conjunta de \((\lambda(\tau),\beta(\tau))\) requer mais esforço.

Machado and Mata (2000) propõem estimar o modelo não linear, \[ Q_{Y}(\tau \, | \, X)=h^{-1}_\lambda \big(x_i^\top \beta(\tau)\big), \] onde \(h^{-1}(z)=(\lambda z+1)^{1/\lambda}\) e Fitzenberger, Wilke and Zhang (2009) sugerem uma modificação para levar em conta circunstâncias em que \(\lambda x_i^\top \beta+1<0\).

Mais recentemente, Mu and He (2007) propuseram um estimador alternativo no qual \(\lambda\) é estimado minimizando uma soma de resíduos cumsum quadrados. O desempenho desses métodos é sensível à heterogeneidade da densidade condicional da resposta, como seria de esperar com base na teoria de grandes amostras que já esboçamos.

Uma terceira opção que parece não ter sido explorada na literatura é simplesmente redimensionar a resposta dividindo por sua média geométrica \(\widetilde{y}_i = y_i/\widetilde{y}\), onde \(\widetilde{y}= \big(\prod y_i\big)^{1/n}\) e então estimar \(\lambda\) resolvendo, \[ \min_{\lambda,\beta}\sum \rho_\tau \big( h(\widetilde{y}_i,\lambda)-x_i^\top\beta \big)\cdot \]

O redimensionamento da resposta leva em conta o termo Jacobiano da transformação Box-Cox da resposta, como na configuração de regressão média convencional. Esta formulação pode fornecer uma densidade condicional mais homogênea em algumas aplicações. Para modelos não lineares em parâmetros mais complicados, é possível usar versões iterativas dos métodos de pontos interiores que fundamentam o ajuste linear de parâmetros, conforme descrito em Koenker and Park (1996).


5. Regressão quantílica não paramétrica


Existe uma extensa literatura sobre regressão quantílica não paramétrica; esta forma relaxa a linearidade estrita nas suposições das covariáveis nos métodos anteriores, preservando ao mesmo tempo a estrutura linear conveniente nos parâmetros que facilita o cálculo eficiente.

Chaudhuri (1991) considera o comportamento assintótico de estimadores de regressão quantílica polinomial local e estabelece condições sob as quais esses estimadores alcançam taxas ótimas de convergência. Posteriormente, o trabalho de Lee (2003) e Lee, Mammen and Park (2010) estendeu a abordagem localmente polinomial para modelos parcialmente lineares e aditivos, respectivamente.

À medida que os modelos de regressão quantílica não paramétricos se tornam mais complexos, o ajuste local e o retroajuste para acomodar novos componentes tornam-se mais onerosos e a literatura evoluiu em direção aos métodos de peneira. Veja, por exemplo, o influente trabalho inicial de Stone (1994) e a pesquisa de Chen (2007). As expansões de funções básicas podem ser adaptadas a aplicações específicas e incorporar mais facilmente certos componentes parcialmente lineares e aditivos.

O desafio é controlar a dimensão paramétrica dos modelos resultantes. Os métodos de penalidade, particularmente a penalidade \(\ell_1\) de Donoho, Chen and Saunders (1998) e Tibshirani (1996), emergiram como ferramentas críticas para a redução da dimensão. O lasso é especialmente conveniente na configuração de regressão quantílica, pois mantém a estrutura de programação linear do problema original.

Esta foi a principal motivação para o uso de penalidades de rugosidade de variação total em Koenker, Ng and Portnoy (1994) e Koenker and Mizera (2004), onde penalidades \(\ell_1\) foram impostas às transformadas lineares dos parâmetros do modelo, controlando efetivamente a variação total das derivadas das funções ajustadas. Um aspecto crucial da estratégia computacional subjacente a esses métodos é a álgebra linear esparsa empregada para representar matrizes de pranejamento de alta dimensão e para resolver sistemas de equações lineares necessários em cada iteração dos algoritmos de pontos interiores usados para ajuste.

Isto é particularmente evidente em aplicações como a de Koenker (2011), onde múltiplos componentes aditivos resultam em vários milhares de parâmetros de modelo. Nesses casos, pode haver vários parâmetros Lagrangianos controlando os componentes não paramétricos aditivos, bem como um lasso mais convencional \(\lambda\) que controla o número efetivo de efeitos covariáveis lineares ativos. Uma variedade de propostas foi feita para como escolher esses parâmetros de penalidade, mas acho que é justo dizer que nenhum consenso foi alcançado.

Esta é uma variante da medida da dimensão Meyer e Woodroofe (2000) \[ \mbox{div}(\widehat{y})=\sum_{i=1}^n \partial \widehat{y}_i/\partial y_i, \] tem somandos que assumem o valor um quando \(y_i\) é interpolado e zero de outra forma na configuração de regressão quantil. No entanto, a proposta de Belloni and Chernozhukov (2011) que constrói uma distribuição de referência para \(\lambda\)’s baseada em uma representação central da condição de gradiente parece ser uma abordagem alternativa muito atraente.

Ainda mais formidável que a \(\lambda\)-seleção é a tarefa de inferência pós-seleção nessas configurações não paramétricas de alta dimensão. Consequentemente, esse tópico gerou considerável pesquisa e controvérsias recentes. A maior parte destes trabalhos se concentrou nos métodos de reamostragem, como exemplificado no contexto de regressão quantil de Belloni, Chernozhukov and Kato (2015) e Belloni, Chernozhukov and Kato (2016). Métodos do tubo de hotelling descritos em Koenker (2011) oferecem uma alternativa atraente para algumas aplicações.

Para ilustrar a abordagem podemos considerar a análise de um modelo autoregressivo simples de primeira ordem para a temperatura máxima diária em Melbourne, Austrália. Os dados são retirados de Hyndman, Bashtannyk and Grunwald (1996). Na Figura 3 fornecemos um gráfico de dispersão de 10 anos de dados de temperatura diária: a temperatura máxima diária de hoje é plotada em relação à máxima de ontem.

library(quantreg); library(splines)
data(MelTemp)
head(MelTemp)
## [1] 38.1 32.4 34.5 20.7 21.5 23.1
par(mar=c(4,4,1,1))
Temp.ontem = MelTemp[1:(length(MelTemp)-1)]
Temp.hoje = MelTemp[2:length(MelTemp)]
plot(Temp.ontem, Temp.hoje, xlab="Temperatura máxima de ontem", 
     asp = 0, xlim = c(5,45), ylim = c(5,45), ylab="Temperatura máxima de hoje",type="n")
points(Temp.ontem,Temp.hoje, cex=0.6,pch=19)
abline(a=0,b=1,lty= 3, lwd = 2)
for(i in seq(0.05,0.95,by=0.05)){
  fig2.ajuste =  rq(Temp.hoje ~ bs(Temp.ontem, df =3), tau = i)
  points(Temp.ontem,fig2.ajuste$fitted.values, col= "red", pch = 19, cex = 0.1, type = "p")
}
grid()

Figura 3: Temperatura máxima diária de Melbourne: o gráfico ilustra 10 anos de dados diários de temperatura máxima (em graus centígrados) para Melbourne, Austrália, como um gráfico de dispersão \(AR(1)\). Os dados estão espalhados pela linha (tracejada) de 45 graus, sugerindo que hoje é aproximadamente semelhante a ontem. Sobrepostas ao gráfico de dispersão estão funções quantílicas condicionais estimadas para os quantis \(\tau\in\{0.05,0.10,\cdots,0.95\}\). Observe que quando a temperatura de ontem está alta, o espaçamento entre as curvas de quantis adjacentes é mais estreito em torno da linha de 45 graus e em cerca de 20 graus centígrados do que na região intermediária. Isto sugere bimodalidade da densidade condicional no verão.

A nossa primeira observação do gráfico é que existe uma forte tendência para os dados se agruparem ao longo da linha (tracejada) de 45 graus, o que implica que com alta probabilidade o máximo de hoje está próximo do máximo de ontem. Mas um exame mais detalhado do gráfico revela que esta impressão se baseia principalmente no lado esquerdo do gráfico, onde a tendência central da dispersão segue muito de perto a linha de 45 graus.

No lado direito, porém, correspondendo às condições de verão, o padrão é mais complicado. Lá, parece que ou há outro dia quente, caindo novamente ao longo da linha de 45 graus, ou há um resfriamento dramático. Mas um leve resfriamento parece ser mais raro. Na linguagem das densidades condicionais, se hoje está quente, a temperatura de amanhã parece ser bimodal, com uma moda aproximadamente centrada no máximo de hoje e a outra moda centrada em cerca de 20°C.

Várias curvas de regressão quantílica estimadas foram sobrepostas ao gráfico de dispersão. Cada curva é especificada como uma B-spline linear. Em condições de inverno, estas curvas estão agrupadas em torno da linha de 45 graus, no entanto, no verão, parece que as curvas do quantil superior estão agrupadas em torno da linha de 45 graus e em torno de 20.

Nas temperaturas intermediárias o espaçamento das curvas quantílicas é um pouco maior indicando menor probabilidade desta faixa de temperatura. Esta impressão é reforçada ao considerar uma sequência de gráficos de densidade baseados nas estimativas de regressão quantílica. Dada uma família de funções quantílicas condicionais estimadas razoavelmente densamente espaçadas, é simples estimar a densidade condicional da resposta em vários valores da covariável de condicionamento. Na Figura 4 ilustramos esta abordagem com várias estimativas de densidade baseadas nos dados de Melbourne.

par(mar=c(4,4,1,1), mfrow=c(2,3))
plot(density(Temp.hoje[Temp.ontem == 11]), bw = "SJ",
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 11"); grid()
plot(density(Temp.hoje[Temp.ontem == 16]), bw = "SJ", 
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 16"); grid()
plot(density(Temp.hoje[Temp.ontem == 21]),  bw = "SJ",
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 21"); grid()
plot(density(Temp.hoje[Temp.ontem == 25]),  bw = "SJ",
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 25"); grid()
plot(density(Temp.hoje[Temp.ontem == 30]),  bw = "SJ",
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 30"); grid()
plot(density(Temp.hoje[Temp.ontem == 35]),  bw = "SJ",
     xlab = "Temperatura máxima hoje", main = "Temperatura ontem 35"); grid()

Figura 4. Estimativas da função de densidade condicional da temperatura máxima de hoje para vários valores da temperatura máxima de ontem com base nos dados de Melbourne. Observe que a temperatura é trimodal quando ontem estava quente.

Condicionado em uma temperatura baixa (11°C) do dia anterior, vemos uma boa densidade condicional bimodal para a temperatura máxima do dia seguinte. Condicionado em uma temperatura moderada do dia anterior, vemos uma boa densidade condicional unimodal para a temperatura máxima do dia seguinte., mas à medida que a temperatura do dia anterior aumenta, vemos uma tendência para a cauda inferior se alongar e, eventualmente, vemos uma densidade tendendo à bimodal e, no caso extremo, trimodal. Neste exemplo, a suposição clássica da regressão de que as covariáveis afetam apenas a locação da distribuição resposta, mas não a sua escala ou forma, é violada.


6. Modelos de séries temporais


Os modelos de séries temporais lineares de coeficientes constantes têm desempenhado um papel de enorme sucesso nas estatísticas e, gradualmente, várias formas de modelos de séries temporais de coeficientes aleatórios também emergiram como concorrentes viáveis em campos específicos de aplicação. Uma variante desta última classe de modelos, embora talvez não seja imediatamente reconhecível como tal, é o modelo de regressão linear quantílica.

Este modelo tem recebido considerável atenção na literatura teórica e pode ser estimado com os métodos de regressão quantílica propostos em Koenker and Bassett (1978). Aqui procuramos relaxar restrições consideradas no artigo mencionado e considerar modelos lineares de autorregressão quantílica cujos parâmetros autorregressivos, ou seja, de inclinação podem variar com os quantis \(\tau\in(0,1)\). Esperamos que estes modelos possam expandir as opções de modelagem para séries temporais que apresentam dinâmica assimétrica ou persistência local.

Um esforço considerável de pesquisa recente tem sido dedicado a modificações dos tradicionais modelos dinâmicos de coeficientes constantes para incorporar uma variedade de efeitos de inovação heterogêneos. Uma motivação importante para tais modificações é a introdução de assimetrias na dinâmica dos modelos. É amplamente reconhecido que muitas variáveis económicas importantes podem apresentar trajetórias de ajuste assimétricas, por exemplo, Neftci (1984), Enders and Granger (1998).

A observação de que as empresas estão mais propensas a aumentar do que a reduzir os preços é uma característica fundamental de muitos modelos macroeconómicos. Beaudry and Koop (1993) argumentaram que os choques positivos no PIB dos EUA são mais persistentes do que os choques negativos, indicando uma dinâmica assimétrica do ciclo económico em diferentes quantis do processo de inovação. Além disso, embora seja geralmente reconhecido que as flutuações do produto são persistentes, resultados menos persistentes também são encontrados em horizontes mais longos (Beaudry and Koop, 1993), sugerindo alguma forma de “persistência local”.

Existe uma literatura teórica substancial tratando do modelo de autorregressão quantílica linear. Neste modelo, a \(\tau\)-ésima função quantílica condicional da resposta \(Y_t\) é expressa como uma função linear dos valores defasados da resposta. Aqui estudamos a estimação e inferência em uma classe mais geral de modelos autoregressivos quantílicos (QAR) nos quais todos os coeficientes autorregressivos podem ser dependentes de \(\tau\) e, portanto, são capazes de alterar a locação, escala e forma das densidades condicionais.


6.1 Modelo de autoregressão quantílica


Seja \(\{U_t\}\) uma sequência de variáveis aleatórias padrão independentes identicamente distribuidas, considere o processo autoregressivo de order \(p\) como sendo \[\begin{equation} \tag{6.1} y_t = \theta_0(U_t)+\theta_1(U_t)y_{t-1}+\cdots+\theta_p(U_t)y_{t-p}, \end{equation}\] onde os \(\theta_j\) são funções desconhecidas tais que \(\theta_j\, : \, [0,1]\to\mathbb{R}\) que desejaremos estimar.

Desde que o lado direito de (6.1) seja monótono crescente em \(U_t\), segue-se que a \(\tau\)-ésima função quantílica condicional de \(y_t\) pode ser escrita como, \[\begin{equation} \tag{6.2} Q_{y_t}(\tau \, | \, y_{t-1},\cdots,y_{t-p})=\theta_0(\tau)+\theta_1(\tau)y_{t-1}+\cdots+\theta_p(\tau)y_{t-p}, \end{equation}\] ou um pouco mais compacto como, \[\begin{equation} \tag{6.3} Q_{y_t}(\tau \, | \, \mathcal{F}_{t-1})=\pmb{x}_t^\top \theta(\tau), \end{equation}\] onde \(\pmb{x}_t = (1,y_{t−1},\cdots, y_{t−p})^\top\) e \(\mathcal{F}_{t}\) é a \(\sigma\)-álgebra gerada por \(\{y_s, s\leq t\}\).

A transição de (6.1) para (6.2) é uma consequência imediata do fato de que para qualquer função monótona crescente \(g\) e uma variável aleatória uniforme padrão \(U\), temos \[ Q_{g(U)}(\tau)=g\big(Q_U(\tau)\big)=g(\tau), \] onde \(Q_U(\tau) =\tau\) é a função quantílica de \(U\). No modelo acima, os coeficientes autorregressivos podem ser dependentes de \(\tau\) e, portanto, podem variar ao longo dos quantis.

As variáveis condicionantes não apenas mudam a locação da distribuição de \(y_t\), mas também podem alterar a escala e a forma da distribuição condicional. Iremos nos referir a este modelo como modelo QAR(p).

Argumentaremos que os modelos QAR podem desempenhar um papel útil na expansão do território de modelagem entre modelos clássicos de séries temporais lineares estacionárias e suas alternativas de raiz unitária. Para ilustrar isso no caso QAR(1), considere o modelo \[\begin{equation} \tag{6.4} Q_{y_t}(\tau \, | \, \mathcal{F}_{t-1}) =\theta_0(\tau)+\theta_1(\tau)y_{t-1}, \end{equation}\] com \(\theta_0=\sigma \Phi^{-1}(\tau)\) e \(\theta_1(\tau)=\min\{\gamma_0+\gamma_1\tau, 1\}\) para \(\gamma_0\in(0,1)\) e \(\gamma_1>0\).

Neste modelo se \(U_t> (1−\gamma_0)/\gamma_1\) o modelo gera \(Y_t\) de acordo com o modelo gaussiano de raiz unitária padrão, mas para realizações menores de \(U_t\) temos uma tendência de reversão à média. Assim, o modelo exibe uma forma de persistência assimétrica no sentido de que sequências de inovações fortemente positivas tendem a reforçar o seu comportamento semelhante à raiz unitária, enquanto realizações negativas induzem reversão à média e, assim, prejudicam a persistência do processo. O modelo AR(1) gaussiano clássico é obtido definindo \(\theta_1(\tau)\) como uma constante.

set.seed(646)
# Gerar uma distribuição lognormal
taus = 3:48/50
n = 90 # comprimento da série temporal para cada quantil
U = rlnorm(n*length(taus), meanlog = 0, sdlog = 1)
sigma = 2.0
gama0 = 0.7
gama1 = 0.5
teta0 = sigma*qnorm(taus)
teta1 = pmin(gama0+gama1*taus,1)
# Modelo
Y = array(U, dim = c(n,length(taus)))
for(i in 2:n)
  {
  for(k in 1:length(taus))
    {
      Y[i,k] = teta0[k]+teta1[k]*lag(Y[i])
    }
}
par(mfrow=c(1,2), mar = c(4,4,1,1))
plot(ts(Y[,5]), ylim = c(-5,5), ylab = expression(paste(tau," = 0.14")), lty = 5, xlab = "Tempo")
points(Y[,5], pch =19,cex = 0.4);grid()
plot(ts(Y[,40]), ylim = c(-5,5), ylab = expression(paste(tau," = 0.84")), lty = 5, xlab = "Tempo")
points(Y[,40], pch =19, cex = 0.4);grid()

Figura 6.1. Esquerda: série temporal AR(1) simulada para \(\tau=0.14\). Esquerda: série temporal AR(1) simulada para \(\tau=0.84\).

A formulação em (6.4) revela que o modelo pode ser interpretado como uma forma especial de modelo autorregressivo de coeficiente aleatório ou modelos RCAR. Tais modelos surgem naturalmente em muitas aplicações de séries temporais. Em contraste com a maior parte da literatura sobre modelos RCAR, nos quais os coeficientes são normalmente assumidos como estocasticamente independentes uns dos outros, o modelo QAR possui coeficientes que são funcionalmente dependentes.

A monotonicidade das funções quantílicas condicionais impõe alguma disciplina nas formas assumidas pelas funções \(\theta\). Esta disciplina requer essencialmente que a função \(Q_{y_t}(\tau \, | \, y_{t−1},\cdots,y_{t−p})\) seja monótona em \(\tau\) em alguma região relevante \(\Upsilon\) do espaço \((y_{t−1},\cdots,y_{t−p})\).

A correspondência entre a formulação do coeficiente aleatório do modelo QAR (6.1) e a formulação da função quantílica condicional (6.2) pressupõe a monotonicidade desta última em \(\tau\). Na região \(\Upsilon\) onde esta monotonicidade se mantém (6.1) pode ser considerado um mecanismo válido para simulação a partir do modelo QAR (6.2). É claro que o modelo (6.1) pode, mesmo na ausência desta monotonicidade, ser tomado como um mecanismo válido de geração de dados, porém a ligação ao modelo quantílico condicional estritamente linear não é mais válida.

Nos pontos onde a monotonicidade é violada, as funções do quantil condicional correspondentes ao modelo descrito por (6.1) apresentam “torções” lineares. Tentar ajustar tais modelos lineares por partes com especificações lineares pode ser perigoso.


6.2 Propriedades do processo QAR


O modelo QSAR(p) em (6.1) pode ser reformulado em notação de coeficientes aleatórios mais convencional como, \[\begin{equation} \tag{6.5} y_t = \mu_0 +\alpha_{1,t} y_{t-1}+\cdots+\alpha_{p,t}y_{t-p}+u_t \end{equation}\] onde \(\mu_0 = \mbox{E}_{\theta_0}(U_t)\), \(u_t = \theta_0(U_t)−\mu_0\), e \(\alpha_{j,t} = \theta_j(U_t)\), para \(j = 1,\cdots,p\). Assim, \(\{u_t\}\) é uma sequência de variáveis aleatórias iid com função de distribuição \(F(\cdot) = \theta_0^{−1}(\cdot + \mu_0)\) e os coeficientes \(\alpha_{j,t}\) são funções da variável aleatória de inovação \(u_t\).

O processo QAR(p) (6.5) pode ser expresso como um processo de autorregressão vetorial \(p\)-dimensional de ordem 1: \[ Y_t = \Gamma + A_t Y_{t-1}+V_t, \] onde \[ \Gamma = \begin{pmatrix} \mu_0 \\ 0_{p-1} \end{pmatrix}, \; A_t = \begin{pmatrix} A_{p-1,t} & \alpha_{p,t} \\ \pmb{I}_{p-1} & 0_{p-1} \end{pmatrix}, \; V_t =\begin{pmatrix} u_t \\ 0_{p-1} \end{pmatrix}, \] sendo \(A_{p-1,t}=(\alpha_{1,t},\dots,\alpha_{p-1,t})\), \(Y_t=(y_t,\cdots,y_{t-(p-1)})^\top\) e \(0_{p-1}\) é o vetor \((p -1)\)-dimensional de zeros. Mostraremos que, sob condições de regularidade dadas no Teorema a seguir, uma solução \(\mathcal{F}_t\)-mensurável para (6.5) pode ser encontrada, esta solução é encontrada durante a demonstração do Teorema 6.1.

Para formalizar a discussão anterior e facilitar a análise assintótica posterior, introduzimos as seguintes condições.



Teorema 6.1. Sob as suposições A.1 e A.2, a série temporal \(y_t\) dada por (6.5) é de covariância estacionária e satisfaz o Teorema do Limite Central \[ \dfrac{1}{\sqrt{n}} \sum_{t=1}^n (y_t-\mu_y) \underset{n\to\infty}{\longrightarrow} N(0,\omega_y^2), \] onde \(\displaystyle\mu_y=\mu_0/\big(1-\sum_{j=1}^p \mu_j\big)\), \(\displaystyle\omega_y^2=\lim_{n\to\infty} \mbox{E}\big( \sum_{t=1}^n (y_t-\mu_t)^2\big)\) e \(\mu_j=\mbox{E}(\alpha_{j,t})\), \(j=1,\cdots,p\).


Demonstração Dado um processo de autorregressão de ordem \(p\) (6.5), denotamos \(\mbox{E}(\alpha_{j,t})=\mu_j\) e assumimos que \(1-\sum_j \mu_j\neq 0\). Seja \(\displaystyle\mu_y=\mu_0/\big(1-\sum_{j=1}^p \mu_j\big)\) e denotamos \[ \underline{y}_t=y_t-\mu_y \] e \[\begin{equation} \tag{6.6} \underline{y}_t=\alpha_{1,t}\underline{y}_{t-1}+\cdots+\alpha_{p,t}\underline{y}_{t-p}+\nu_t \end{equation}\] onde \[ \nu_t=u_t+\mu\sum_{\ell=1}^p (\alpha_{\ell,t}-\mu_\ell)\cdot \] Pode-se ver que \(\mbox{E}(\nu_t) = 0\) e \(\mbox{E}(\nu_t\nu_s) = 0\) para qualquer \(t\neq s\) já que \(\mbox{E}(\alpha_{\ell,t}) = \mu_\ell\) e \(u_t\) são independentes. Para derivar condições de estacionariedade para o processo \(\underline{y}_t\), primeiro encontramos uma solução \(\mathcal{F}_t\)-mensurável para (6.6). Definimos os vetores \(p\times 1\) aleatórios \[ \underline{Y}_t=\big(\underline{y}_t,\cdots,\underline{y}_{t-p+1}\big)^\top=(\nu_t,0,\cdots,0)^\top \] e a matriz \(p\times p\) \[ A_t=\begin{pmatrix} A_{p-1,t} & \alpha_{p,t} \\ \pmb{I}_{p-1} & 0_{p-1} \end{pmatrix}, \] onde \(A_{p-1,t}=(\alpha_{1,t},\cdots,\alpha_{p-1,t})\) e \(0_{p-1}\) é um vetor de dimensão \(p-1\) de zeros, então \[ \mbox{E}\big(V_tV_t^\top\big)=\begin{pmatrix} \sigma_\nu^2 & 0_{1\times (p-1)} \\ 0_{(p-1)\times 1} & 0_{(p-1)\times (p-1)}\end{pmatrix} = \Sigma \] e o processo original pode ser escrito como \[ \underline{Y}_t=A_t\underline{Y}_{t-1}+V_t\cdot \]

\(\blacksquare\)

Para ilustrar algumas características importantes do processo QAR, consideramos o caso mais simples do processo QAR(1), \[\begin{equation} \tag{6.7} y_t=\alpha_ty_{t-1}+u_t, \end{equation}\] onde \(\alpha_t=\theta_1(U_t)\) e \(u_t=\theta_0(U_t)\) correspondente a (6.4), cujas propriedades estão resumidas no seguinte corolário.


Corolário 6.1. Se \(y_t\) for determinado por (6.6) e \(\displaystyle\omega_\alpha^2 = \mbox{E}(\alpha_t)^2 < 1\), sob a suposição A.1, \(y_t\) é de covariância estacionária e satisfaz o Teorema do Limite Central \[ \dfrac{1}{\sqrt{n}}\sum_{t=1}^n y_t \underset{n\to\infty}{\longrightarrow} N(0,\omega_y^2), \] onde \(\omega_y^2=\sigma^2(1+\mu_\alpha)/\big((1-\mu_\alpha)(1-\omega_\alpha^2)\big)\) com \(\mu_\alpha=\mbox{E}(\alpha_t)<1\).


No exemplo dado na Seção 2.1, \(\alpha_t = \theta_1(U_t) = \min\{\gamma_0 +\gamma_1 U_t,1\}\leq 1\), e \(P(|\alpha_t| < 1) > 0\), a condição do Corolário 6.1 é válida e o processo \(y_t\) é globalmente estacionário, mas ainda pode apresentar persistência local e assimétrica na presença de certos tipos de choques, choques positivos, no exemplo.

O Corolário 6.1 acima também indica que, mesmo com \(\alpha_t > 1\) em algum intervalo de quantis, \(y_t\) ainda pode ser de covariância estacionária no longo prazo, desde que \(\omega_\alpha^2 = \mbox{E}(\alpha_t)^2 < 1\). Assim, um processo autoregressivo quantílico pode permitir algumas formas transitórias de comportamento explosivo, mantendo a estacionariedade no longo prazo.

Sob as suposições do Corolário 6.1, substituindo recursivamente em (6.6), podemos ver que \[\begin{equation} \tag{6.8} y_t=\sum_{j=0}^\infty \beta_{t,j}u_{t-j}, \; \mbox{onde} \; \beta_{t,0}=1, \; \beta_{t,j}=\prod_{i=0}^{j-1}\alpha_{t-i}, \; \mbox{para} \; j\leq 1, \end{equation}\] é uma solução estacionária \(\mathcal{F}_t\)-mensurável para (6.7). Além disso, se \(\sum_{j=0}^\infty \beta_{t,j} \nu_{t-j}\) converge em \(L^p\), então \(y_t\) tem um momento finito de \(p\)-ésima ordem.

A solução \(\mathcal{F}_t\)-mensurável de (6.7) fornece uma representação MA(\(\infty\)) duplamente estocástica de \(y_t\). Em particular, a resposta ao impulso de \(y_t\) a um choque \(u_{t−j}\) é estocástica e é dada por \(\beta_{t,j}\). Por outro lado, embora a resposta ao impulso do processo autoregressivo quantílico seja estocástica, ela converge em quadrado médio para zero e, portanto, em probabilidade quando \(j\to\infty\), corroborando a estacionariedade de \(y_t\).

Se denotarmos a função de autocovariância de \(y_t\) por \(\gamma_y(h)\), pode-se verificar que \[ \gamma_y(h) = \mu_\alpha^{|h|}\sigma_y^2 \] onde \(\sigma_y^2 = \sigma^2/(1 − \omega_\alpha^2)\).

Comparando com o processo \(y_t\) QAR(1), se considerarmos um processo convencional AR(1) com coeficiente autoregressivo \(\mu_\alpha\) e denotando o processo correspondente por \(\underline{y}_t\), a variância de longo prazo de \(y_t\), dada por \(\omega_y^2\) é maior que a de \(\underline{y}_t\). A variância adicional do processo QAR \(y_t\) vem da variação de \(\alpha_t\). De fato, \(\omega_y^2\) pode ser decomposto na soma da variância de longo prazo de \(y_t\) e um termo adicional que é determinado pela variância de \(\alpha_t\), \[ \omega_y^2=\underline{\omega}_y^2+\dfrac{\sigma^2}{(1-\mu_\alpha)^2(1-\omega_\alpha^2)}\mbox{Var}(\alpha_t), \] onde \(\underline{\omega}_y^2=\sigma^2/(1-\mu_\alpha)^2\) é a variancia de longo prazo de \(y_t\).

Quando \(\alpha_i(\tau)\) não depende de \(\tau\) para \(i=1,\cdots,q\), temos o familiar modelo de autorregressão do erro iid, com erros tendo função quantílica \(\alpha_0\).

De forma mais geral, temos um modelo quantílico de coeficientes aleatórios QAR(q) com \(Y_t\) gerado como, \[ Y_t=\alpha_0(U_t)+\sum_{i=1}^q \alpha_i(U_t)Y_{t-i}, \] com \(U_t\sim U(0,1)\).

No entanto, este é um modelo de coeficientes aleatórios bastante especial, uma vez que todos os seus coeficientes são conduzidos pelas mesmas variáveis aleatórias uniformes iid. Na terminologia de Schmeidler (1986) os coeficientes são comonotônicos.

O caso mais simples do modelo QAR(1),


6.3 Estimação


A estimação do modelo autoregressivo quantílico (6.3) envolve a resolução do problema \[\begin{equation} \tag{6.9} \min_{\theta\in\mathbb{R}^{p+1}} \sum_{t=1}^n \rho_\tau \big(y_t-x^\top_t \theta\big), \end{equation}\] onde \(\rho_\tau(u)=u\big(\tau-\pmb{I}(u<0)\big)\) como em Koenker and Bassett (1978). As soluções \(\widehat{\theta}(\tau)\) são chamados quantis de autorregressão.

Dado \(\widehat{\theta}(\tau)\), a \(\tau\)-ésima função quantílica condicional de \(y_t\), condicionada em \(x_t\), pode ser estimada por \[ \widehat{Q}_{y_t}(\tau \, | \, x_t) = x_t^\top \widehat{\theta}(\tau), \] e a densidade condicional de \(y_t\) pode ser estimada pelos quocientes de diferença, \[ \widehat{f}_{y_t}(\tau \, | \, x_{t-1}) = (\tau_i-\tau_{i-1})\Bigm/ \big(\widehat{Q}_{y_t}(\tau_i \, | \, x_{t-1}) - \widehat{Q}_{y_t}(\tau_{i-1} \, | \, x_{t-1}) \big), \] para alguma sequência apropriadamente escolhida de \(\tau\)’s.

Denotando \(\mbox{E}(y_t)=\mu_t\), \(\mbox{E}(y_ty_{t-j})=\gamma_j\) e \[ \Omega_0 =\mbox{E}(x_tx_t^\top)=\lim_{n\to\infty} \dfrac{1}{n}\sum_{t=1}^n x_tx_t^\top, \] então \[ \Omega_0 = \begin{pmatrix} 1 & \mu_y^\top \\ \mu_y^\top & \Omega_y \end{pmatrix}, \] onde \(\mu_y=\mu_y\times \pmb{1}_{p\times 1}\) e \[ \Omega_y = \begin{pmatrix} \gamma_0 & \cdots & \gamma_{p-1} \\ \vdots & \ddots & \vdots \\ \gamma_{p-1} & \cdots & \gamma_0 \end{pmatrix}\cdot \]

No caso especial do modelo QAR(1) (6.7), \(\Omega_0=\mbox{E}(x_tx_t^\top)=\mbox{diag}(1,\gamma_0)\), sendo \(\gamma_0=\mbox{E}(y_t^2)\). Seja \[ \Omega_1=\lim_{n\to\infty} \dfrac{1}{n}\sum_{t=1}^n f_{t-1}\Big(F^{-1}_{t-1}(\tau)\Big) x_tx_t^\top \] e definimos \(\Sigma=\Omega_1^{-1}\Omega_0\Omega_1^{-1}\). A distribuição assintótica de \(\widehat{\theta}(\tau)\) está resumida no seguinte Teorema.


Teorema 6.2. Sob as suposições A.1, A.2 e A.3, \[ \Sigma^{-1/2}{\sqrt{n}} \big( \widehat{\theta}(\tau)-\theta(\tau)\big) \underset{n\to\infty}{\longrightarrow} B_k(\tau), \] onde \(B_k(\tau)\) ewpresenta uma ponte Browniana padrão \(k\)-dimensional, \(k=p+1\).


Demonstração Dado um processo de autorregressão de ordem \(p\) (6.5), denotamos

\(\blacksquare\)


6.3.1 Exemplo


Aqui você encontrará uma breve demonstração de coisas que você pode fazer com a autorregressão quantílica em R. Os dados para este tutorial são o Índice de Miséria da Zona Euro, Índice de Miséria é um indicador econômico, criado pelo economista Arthur Okun. O índice ajuda a determinar como o cidadão médio está se saindo economicamente e é calculado adicionando a taxa de desemprego ajustada sazonalmente à taxa de inflação anual. A Zona Euro, oficialmente designada por Área do Euro e também coloquialmente referida como Eurozona, refere-se a uma união monetária constituída por 20 Estados-membros da União Europeia que adotaram oficialmente o euro como moeda comum.

Os dados considrados formam uma série temporal de frequência mensal cuja soma: (taxa de desemprego + taxa de inflação) compõe o chamado “Índice de Miséria”. Começamos carregando os dados e verificando o que os diferentes critérios de informação dizem sobre o número de defasagens no modelo:

options(digits = 4)
library(quantreg) # Para a regressão quantílica
dados = read.table("mindex-08-12.txt",sep = "\t",header = FALSE)
y = ts(rev(dados$V2), frequency = 12,  start = c(1995,1))
par(mar=c(4,4,3,1), cex = 0.9, pch = 19)
plot(y, main = "Índice de Miséria", ylab = "Índice", lwd = 2, ty = "b", cex.main = 2, xlab = "Tempo")
grid()

lagorder = function(x,maxorder){ # maxorder é o número máximo de atraso (lag) que você permite
  ordmax = maxorder ; AK = NULL ; SC = NULL ; lagmat = NULL
  l1 = length(x)
  for (i in 1:ordmax){
    lagmat = cbind(lagmat[-i,],x[(1):(l1-i)]) # lagged matrix
    armod <- lm(x[(i+1):l1]~lagmat)
    AK[i] = AIC(armod)
    SC[i] = BIC(armod)
  }
  return(c(which.min(AK),which.min(SC) ))
}
lagorder(y, maxorder = 12) 
## [1] 1 1


Todos concordam que um atraso é suficiente. Agora estimamos a autorregressão quantílica, uma para cada quantil em incrementos de 0.05. Esta estimação é de ordem 1; segundo obtido acima o atraso ou lag de ordem 1 é suficiente para a estimação. Então eliminamos o primeiro e o último elemento da série e estimamos o modelo médio, ou seja, o modelo de regressão linear em lm0. Depois disto ajustamos os modelos de autoregressão quantílicos na sequência de quantis selecionados.

TT = length(y)
lm0 = lm(y[-1]~y[-TT]) ; summary(lm0)
## 
## Call:
## lm(formula = y[-1] ~ y[-TT])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8691 -0.1427 -0.0039  0.1240  0.8092 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.1686     0.1667    1.01     0.31    
## y[-TT]        0.9856     0.0146   67.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.232 on 209 degrees of freedom
## Multiple R-squared:  0.956,  Adjusted R-squared:  0.956 
## F-statistic: 4.58e+03 on 1 and 209 DF,  p-value: <2e-16
tauseq = seq(.05,.95,.05)
qslope = NULL ; qr0 = list()
for (i in 1:length(tauseq)){
  qr0[[i]] = rq(y[-1]~y[-TT], tau = tauseq[i])
  qslope[i] = qr0[[i]]$coef[2]
}


Na figura a continuação mostramos, na parte superior, as linhas ajustadas estão sobrepostas em cinza. No caso de coeficientes AR constantes (constantes entre quantis), devemos obter linhas paralelas entre si, pois a única mudança é o quantil que você deseja ajustar. Neste caso, podemos ver no painel inferior direito que os coeficientes AR não são constantes. Para ajuste de quantis baixos, o processo se comporta como um passeio aleatório, enquanto uma forte reversão à média é observada para quantis altos. Esta assimetria sugere que o processo é heterocedástico, a variância na extremidade inferior é maior do que na extremidade superior, então obtemos uma figura semelhante a um “leque” em vez de linhas paralelas. O significado económico disto é que quando este índice é elevado, são tomadas medidas para o diminuir. Medidas como a redução do custo dos empréstimos que, por sua vez, permitem que as empresas em dificuldades se mantenham à tona e que as empresas bem-sucedidas mantenham os níveis de investimento. A questão é que quando o índice está alto, tentamos deprimi-lo, enquanto quando está na faixa intermediária ele pode ir para os dois lados, daí o formato de “leque”.

par(mar=c(4,4,1,1), pch = 19, cex = 0.9)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(y[-1]~y[-TT], col = 1, ylab = "Índice", xlab = "(Lagged) Ínide de Miséria")
for (i in 1:length(tauseq)){
  lines(y[-TT],qr0[[i]]$fit, lwd = .02, col  = rgb(0, 0, 0, 0.2) )
}
grid()
plot(y, ylab = "Índice", main  = "Índice de Miséria - Zona do Euro", 
     lwd = 1, ty = "b", xlab = "Tempo"); grid()
plot(qslope~tauseq, ty = "b", xlab = "Quantis", ylab = "Coeficiente AR"); grid()


Notas (1) Há menos dispersão quando o índice de defasagem é maior e mais dispersão quando o índice de defasagem é menor. Esta é uma autorregressão quantílica condicional, ela não responde à pergunta “o que acontecerá com a dispersão do índice quando eu alterar a defasagem em uma pequena quantidade?”, para responder a esse tipo de pergunta você precisará ler sobre regressão quantílica incondicional, o que não é discutido aqui. (incondicional ao nível do índice..)

  1. As linhas ajustadas não devem se cruzar como fazem no gráfico acima. Quando isso acontece a interpretação é muito problemática, existem maneiras de evitar esses cruzamentos para melhorar a interpretação.

Num outro contexto, pode experimentar esta ferramenta para estimar o Valor em Risco, com quantil igual a 0,05 para o valor de VaR de 5%. Tenha em mente que neste caso você precisa de uma amostra grande para precisão, pois apenas 5% das observações terão informações relevantes para determinar o valor ajustado. Talvez seja interessante ver como a estimativa do VaR da regressão quantílica se compara com o garch(1,1) comum etc.

O que eu gosto especialmente na ferramenta de regressão quantílica é que a única suposição é muito “leve”, apenas que a forma funcional é linear, sem gaussianidade ou algo assim, então é muito geral. Plus está pronto para usuários R graças a Roger Koenker.


10. Referências


  1. Abadie, A., J. Angrist, and G. Imbens (2002): “Instrumental Variables Estimates of Subsidized Training on the Quantile of Trainee Earnings,” Econometrica, 70, 91–117.
  2. Arellano, M., and S. Bonhomme (2015): “Quantile Selection Models: With an Application to Understanding Changes in Wage Inequality,” preprint.
  3. ———– (2016): “Nonlinear Panel Data Estimation via Quantile Regression,” Econometrics Journal.
  4. Arrow, K., and M. Hoffenberg (1959): A Time Series Analysis of Interindustry Demands. North-Holland, Amsterdam.
  5. Bassett, Gilbert W., J., M.-Y. S. Tam, and K. Knight (2002): “Quantile models and estimators for data analysis,” Metrika, 55, 17–26.
  6. Beaudry, P. and Koop, G. (1993). Do recessions permanently change output?, Journal of Monetary Economics, 31, 149-163.
  7. Belloni, A., and V. Chernozhukov (2011): “`1 -penalized quantile regression in high-dimensional sparse models,” The Annals of Statistics, 39, 82–130.
  8. Belloni, A., V. Chernozhukov, and K. Kato (2015): “Uniform post-selection inference for least absolute deviation regression and other Z-estimation problems,” Biometrika, 102, 77–94.
  9. ———– (2016): “Valid post-selection inference in High-Dimensional Approximately Sparse Quantile Regression Models,” arXiv:1312.7186v4.
  10. Buchinsky, M. (2001): “Quantile regression with sample selection: Estimating women’s return to education in the U.S.,” Empirical Economics, 26, 87–113.
  11. Carlier, G., V. Chernozhukov, and A. Galichon (2016): “Vector quantile regression: An optimal transport approach,” Ann. Statist., 44, 1165–1192.
  12. Chamberlain, G. (1984): “Panel data,” in Handbook of Econometrics, vol. 2, pp. 1247–1318. Elsevier.
  13. Chamberlain, G. (1994): “Quantile regression, censoring and the structure of wages,” in Advances in Econometrics, ed. by C. Sims. Elsevier: New York.
  14. Chaudhuri, P. (1991): “Nonparametric Estimates of Regression Quantiles and Their Local Bahadur Representation,” The Annals of Statistics, 19, 760–777.
  15. Chen, X. (2007): “Large Sample Sieve Estimation of Semi-Nonparametric Models,” in Handbook of Econometrics, ed. by J. J. Heckman, and E. E. Leamer, vol. 6, Part B, pp. 5549 – 5632. Elsevier.
  16. Chen, X., R. Koenker, and Z. Xiao (2009): “Copula-based nonlinear quantile autoregression,” Econometrics Journal, 12, S50–S67.
  17. Chernozhukov, V., I. Fernndez-Val, and A. Galichon (2010): “Quantile and Probability Curves Without Crossing,” Econometrica, 78, 1093–1125.
  18. Chernozhukov, V., and C. Hansen (2005a): “The Effects of 401(k) Participation on the wealth distribution: An Instrumental Quantile Regression Analysis,” Review of Economics and Statistics, 73, 735–751.
  19. ———– (2005b): “An IV Model of Quantile Treatment Effects,” Econometrica, 73, 245–262.
  20. ———– (2006): “Instrumental Quantile Regression Inference for structural and treatment effect models,” J. of Econometrics, 132, 491–525.
  21. ———– (2008): “Instrumental Variable Quantile Regression: A Robust Inference Approach,” J. of Econometrics, 142, 379–398.
  22. ———– (2013): “Quantile Models with Endogeneity,” Annual Reviews in Economics, 5, 57–81.
  23. Chesher, A. (2003): “Identification in Nonseparable Models,” Econometrica, 71, 1405–1441.
  24. ———– (2005): “Nonparametric Identification under Discrete Variation,” Econometrica, 73, 1525–1550.
  25. Choudhury, J., and P. Chaudhuri (2017): “Quantile Regression with Measurement Error and Missing Data,” in Handbook of Quantile Regression, ed. by V. Chernozhukov, X. He, R. Koenker, and L. Peng. Chapman-Hall.
  26. Cole, T. J., and P. Green (1992): “Smoothing Reference Centile Curves: The LMS Method and Penalized Likelihood,” Statistics in Medicine, 11, 1305–1319.
  27. Cox, D. R. (1984): “Interaction,” International Statistical Review, 52, 1–24.
  28. Diaz, I. (2016): “Efficient estimation of quantiles in missing data models,” arXiv:1512.08110.
  29. Donoho, D., S. Chen, and M. Saunders (1998): “Atomic decomposition by basis pursuit,” SIAM J. of Scientific Computing, 20, 33–61.
  30. Durbin, J. (1954): “Errors in Variables,” Review of the International Statistical Institute, 22, 23–32.
  31. Edgeworth, F. (1888a): “A Mathematical Theory of Banking,” J. Royal Stat Soc, 51, 113–127.
  32. ———– (1888b): “On a new method of reducing observations relating to several quantities,” Philosophical Magazine, 25, 184–191.
  33. Efron, B. (1967): “The Two Sample Problem with Censored Data,” in Proc. 5th Berkeley Sympt. Math. Statist. Prob. Prentice-Hall: New York.
  34. Enders, W. and Granger, C. (1998). Unit Root tests and asymetric adjustment with an example sing the term structure of interest rates, Journal of Business and Economic Statistics, 16, 304-311.
  35. Firpo, S. (2007): “Efficient Semiparametric Estimation of Quantile Treatment Effects,” Econometrica, 75, 259–276.
  36. Fitzenberger, B., R. Wilke, and X. Zhang (2009): “Implementing BoxCox Quantile Regression,” Econometric Reviews, 29, 158–181.
  37. Fox, M., and H. Rubin (1964): “Admissibility of Quantile Estimates of a Single Location Parameter,” Annals of Mathematical Statistics, 35, 1019–1030.
  38. Freedman, D. A. (2008): “On regression adjustments in experiments with several treatments,” Annals Applied Stat., 2, 176–196.
  39. Galichon, A. (2016): Optimal Transport Methods in Economics. Princeton U. Press.
  40. Galvao, A. F. (2011): “Quantile regression for dynamic panel data with fixed effects,” Journal of Econometrics, 164, 142–157.
  41. Gronau, R. (1974): “Wage Comparisons–A Selectivity Bias,” Journal of Political Economy, 82, 1119–1143.
  42. Hagemann, A. (2011): “Robust Spectral Analysis,” arXiv:1111.1965.
  43. ———– (2016): “Cluster-robust bootstrap inference in quantile regression models,” Journal of the American Statistical Association.
  44. Hallin, M., Z. Lu, D. Paindaveine, and M. Šiman (2015): “Local bilinear multiple-output quantile/depth regression,” Bernoulli, 21, 1435–1466.
  45. Hallin, M., D. Paindaveine, and M. Šiman (2010): “Multivariate quantiles and multiple-output regression quantiles: From L1 optimization to halfspace depth,” Ann. Statist., 38, 635–669.
  46. He, X., and H. Liang (2000): “Quantile regression estimates for a class of linear and partially linear errors-in-variables models,” Statistica Sinica, 10(1), 129–140.
  47. Heckman, J. (1974): “Shadow Prices, Market Wages, and Labor Supply,” Econometrica, 42, 679–694.
  48. Heckman, J. J. (1979): “Sample Selection Bias As a Specification Error,” Econometrica, 47, 153–161.
  49. Heckman, J. J., J. Smith, and N. Clements (1997): “Making the most out of programme evaluations and social experiments: Accounting for heterogeneity in programme impacts,” The Review of Economic Studies, 64, 487–535.
  50. Hotelling, H. (1939): “Tubes and Spheres in n-spaces, and a class of statistical problems,” Am J. Math, 61, 440–460.
  51. Imbens, G., and J. Angrist (1994): “Identification and Estimation of Local Average Treatment Effects,” Econometrica, 62, 467–475.
  52. J. Neyman, E. L. S. (1948): “Consistent Estimates Based on Partially Consistent Observations,” Econometrica, 16, 1–32.
  53. Kadane, J. B., and T. Seidenfeld (1996): “Statistical Issues in the Analysis of Data Gathered in the New Designs,” in Bayesian Methods and Ethics in a Clinical Trial Design, ed. by J. B. Kadane, pp. 115–125. Wiley.
  54. Kato, K. (2012): “Estimation in functional linear quantile regression,” Ann. Statist., 40, 3108–3136.
  55. Kato, K., A. F. Galvao, and G. V. Montes-Rojas (2012): “Asymptotics for panel quantile regression models with individual effects,” Journal of Econometrics, 170, 76–91.
  56. Khmaladze, E. V. (1981): “Martingale Approach in the Theory of Goodness-of-fit Tests,” Theory of Probability and its Applications, 26, 240–257.
  57. Kley, T., S. Volgushev, H. Dette, and M. Hallin (2016): “Quantile spectral processes: Asymptotic analysis and inference,” Bernoulli, 22, 1770–1807.
  58. Knight, K., and G. Bassett (2007): “Second order improvements of sample quantiles using subsamples,” Technical Report: Department of Statistics, U. of Toronto.
  59. Koenker, R. (2004): “Quantile Regression for Longitudinal Data,” Journal of Multivariate Anal- ysis, 91, 74–89.
  60. Koenker, R. (2005): Quantile Regression. Cambridge University Press.
  61. (2011): “Additive models for quantile regression: Model selection and confidence bandaids,” Brazilian J. Probab. Stat., 25, 239–262.
  62. (2016): quantreg: Quantile RegressionR package version 5.27.
  63. Koenker, R. (2017): “Computational Methods for Quantile Regression,” in Handbook of Quantile Regression. Chapman-Hall.
  64. Koenker, R., and G. Bassett (1978): “Regression Quantiles,” Econometrica, 46, 33–50.
  65. Koenker, R., and V. d’Orey (1987): “Computing Regression Quantiles,” Applied Statistics, 36, 383–393.
  66. Koenker, R., and O. Geling (2001): “Reappraising Medfly Longevity: A quantile regression survival analysis,” J. of Am. Stat. Assoc., 96, 458–468.
  67. Koenker, R and J.A.F. Machado (1999). Goodness of fit and related inference processes for quantile regression. J Am Stat Assoc 94: 1296–1310.
  68. Koenker, R., and I. Mizera (2004): “Penalized Triograms: Total Variation Regularization for Bivariate Smoothing,” J. Royal Stat. Society (B), 66, 145–163.
  69. Koenker, R., P. Ng, and S. Portnoy (1994): “Quantile Smoothing Splines,” Biometrika, 81, 673–80.
  70. Koenker, R., and B. Park (1996): “An interior point algorithm for nonlinear quantile regression,” Journal of Econometrics, 71, 265–283.
  71. Koenker, R., and Z. Xiao (2002): “Inference on the quantile regression process,” Econometrica, 70, 1583–1612.
  72. ———– (2006): “Quantile Autoregression, with Discussion and Rejoinder,” Journal of the American Statistical Association, 101, 980–1006.
  73. Kong, L., and I. Mizera (2012): “Quantile tomography: using quantiles with multivariate data,” Statistica Sinica, 22, 1589–1610.
  74. Lamarche, C. (2010): “Robust penalized quantile regression estimation for panel data,” Journal of Econometrics, 157, 396–408.
  75. Lee, S. (2003): “Efficient semiparametric estimation of a partially linear quantile regression model,” ET, pp. 1–31.
  76. Lee, Y. K., E. Mammen, and B. U. Park (2010): “Backfitting and smooth backfitting for additive quantile models,” Ann. Statist., 38, 2857–2883.
  77. Lehmann, E. (1974): Nonparametrics: Statistical Methods based on Ranks. Holden-Day.
  78. Li, R., and L. Peng (2017): “Survival Analysis with Computing Risks,” in Handbook of Quantile Regression. Chapman-Hall.
  79. Li, T.-H. (2008): “Laplace Periodogram for Time Series Analysis,” Journal of the American Statistical Association, 103, 757–768.
  80. ———– (2012): “Quantile Periodograms,” Journal of the American Statistical Association, 107, 765–776.
  81. Lipsitz, S., G. Fitzmaurice, G. Molenberghs, and L. Zhao (1997): “Quantile regression methods for longitudinal data with Drop-outs: Application CD4 cell counts of patients infected with HIV,” Applied Statistics, pp. 463–476.
  82. Ma, L., and R. Koenker (2006): “Quantile regression methods for recursive structural equation models,” Journal of Econometrics, 134, 471–506.
  83. Machado, J., and J. Mata (2000): “Box-Cox Quantile regression and the distribution of firm sizes,” Journal of Applied Econometrics, pp. 253–274.
  84. Manski, C. (1975): “Maximum score estimation of the stochastic utility model of choice,” J. of Eonometrics, 3, 205–228.
  85. Manski, C. (1993): “The Selection Problem in Eocnometrics and Statistics,” in Handbook of Statistics, ed. by C. R. G.S. Maddala, and H. Vinod. Elsevier: New York.
  86. Manski, C. F. (2004): “Statistical Treatment Rules for Heterogeneous Populations,” Econometrica, 72, 1221–1246.
  87. Matzkin, R. L. (2015): “Estimation of Nonparametric Models With Simultaneity,” Econometrica, 83, 1–66.
  88. Meyer, M., and M. Woodroofe (2000): “On the degrees of freedom in shape-restricted regression,” Annals of Statistics, 28, 1083–1104.
  89. Mu, Y., and X. He (2007): “Power Transformation toward a Linear Regression Quantile,” Journal of the American Statistical Association, 102, 269–279.
  90. Neftci, S. (1984). Are economic time series asymmetric over the business cycle?, Journal of Political Economy, 92, 307-328.
  91. Parikh, N., and S. Boyd (2014): “Proximal Algorithms,” Foundations and Trends in Optimization, 1, 123–231.
  92. Peng, L. (2017): “Quantile Regression for Survival Analysis,” in Handbook of Quantile Regression. Chapman-Hall.
  93. Peng, L., and Y. Huang (2008): “Survival Analysis With Quantile Regression Models,” Journal of the American Statistical Association, 103, 637–649.
  94. Portnoy, S. (2003): “Censored Quantile Regression,” Journal of the American Statistical Association, 98, 1001–1012.
  95. Portnoy, S., and R. Koenker (1997): “The Gaussian Hare and the Laplacian Tortoise: Computability of squared-error versus absolute-error estimators, with discusssion,” Statistical Science, 12, 279–300.
  96. Powell, J. L. (1986): “Censored Regression Quantiles,” J. of Econometrics, 32, 143–155.
  97. R Core Team (2016): R: A Language and Environment for Statistical ComputingR Foundation for Statistical Computing.
  98. Ramsay, J., and B. Silverman (1997): Functional Data Analysis. Springer.
  99. Rosenblatt, M. (1952): “Remarks on a Multivariate Transformation,” Ann. Math. Statist., 23, 470–472.
  100. Schennach, S. M. (2008): “Quantile Regression with Mismeasured Covariates,” Econometric Theory, 24, 1010–1043.
  101. Schmeidler, D. (1986): “Integral representation without additivity,” Proceedings of Am Math Society, 97, 255–261.
  102. Stone, C. J. (1994): “The Use of Polynomial Splines and Their Tensor Products in Multivariate Function Estimation,” Annals of Statistics, 22, 118–171.
  103. Strotz, R., and H. Wold (1960): “A Triptych on Causal Systems,” Econometrica, 28, 417–463.
  104. Student (1931): “The Lanarkshire Milk Experiment,” Biometrika, 23, 398–406.
  105. Tibshirani, R. (1996): “Regression Shrinkage and Selection Via the Lasso,” J. Royal Stat. Soc. (B), 58, 267–288.
  106. Wald, A. (1940): “The fitting of straight lines if both variables are subject to error,” Annals of Mathematical Statistics.
  107. Wang, H. J., D. Li, and X. He (2012): “Estimation of High Conditional Quantiles for Heavy-Tailed Distributions,” Journal of the American Statistical Association, 107, 1453–1464.
  108. Wang, H. J., L. A. Stefanski, and Z. Zhu (2012): “Corrected-loss estimation for quantile regression with covariate measurement errors,” Biometrika, 99, 405–421.
  109. Wang, L., Y. Zhou, R. Song, and B. Sherwood (2016): “Quantile Optimal Treatment Regimes,” Technical Report: Department of Statistics, U. of Minnesota.
  110. Wei, Y. (2008): “An Approach to Multivariate Covariate-Dependent Quantile Contours With Application to Bivariate Conditional Growth Charts,” Journal of the American Statistical Association, 103, 397–409.
  111. ———– (2017): “Quantile Regression with Measurement Error and Missing Data,” in Handbook of Quantile Regression, ed. by V. Chernozhukov, X. He, R. Koenker, and L. Peng. Chapman-Hall.
  112. Wei, Y., and R. J. Carroll (2009): “Quantile Regression With Measurement Error,” Journal of the American Statistical Association, 104, 1129–1143.
  113. Wei, Y., A. Pere, R. Koenker, and X. He (2005): “Quantile Regression for Reference Growth Charts,” Statistics in Medicine, 25, 1369–1382.
  114. Yang, X., N. Narisetty, and X. He (2016): “A New Approach to Censored Quantile Regression Estimation,” University of Michigan, Department of Statistics Technical Report.
  115. Ying, Z. (2017): “Analysis of Survival Data and Quantile Regression,” in Handbook of Quantile Regression. Chapman-Hall.