\documentclass[a4paper,10pt]{article} 
\usepackage{amsmath} 
\usepackage{amsfonts} 
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{setspace}
\usepackage{harvard}
\usepackage{aer}
\usepackage{fullpage}
\usepackage{hyperref}
\usepackage{graphicx}

\newcommand{\indep}{\perp\!\!\!\perp}
\newcommand{\argmax}{\operatornamewithlimits{arg\,max}}
\newcommand{\argmin}{\operatornamewithlimits{arg\,min}}
\newcommand{\plim}{\operatornamewithlimits{plim}}
\newcommand{\citefull}[1]{\citename{#1} \citeyear{#1}}
\newcommand{\citeparagraph}[1]{\medskip\noindent\textbf{{\citename{#1} \citeyear{#1}}}}

\newcommand{\cov}{\text{Cov}}
\newcommand{\var}{\text{Var}}
\newcommand{\rank}{\text{rank}}
% \newcommand{\det}{\text{det}}
\def\inprobLOW{\rightarrow_p}
\def\inprobHIGH{\,{\buildrel p \over \rightarrow}\,} 
\def\as{\,{\buildrel a.s. \over \rightarrow}\,} 
\def\asu{\,{\buildrel a.s.u. \over \rightarrow}\,} 
\def\inprob{\,{\inprobHIGH}\,} 
\def\indist{\,{\buildrel d \over \rightarrow}\,} 

% defined environments
\newtheorem{thm}{Theorem} %[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}

\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
\newtheorem{ex}[thm]{Example}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}

\title{14.385 Recitation 10}
\author{Paul Schrimpf}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Local Linear Quantile Regression}

All the nonparametric methods that we have talked about for mean
regression also work for quantile regression.  The steps for showing
consistency, calculating MSE, and finding the asymptotic distribution
for nonparametric quantile regression are similar to the steps for
showing the same things for mean regression.  Here, we will analyze
local linear quantile regression.  

\subsection{Setup}
Let $q(x) = Q_Y(\tau|x)$ be the $\tau$th quantile of $y$ given $x$.
Define the local linear quantile regression estimator as:
\begin{align}
  \hat{q}(x) \argmin_\alpha \min_\beta \frac{1}{n} \sum_{i=1}^n
  \rho_\tau(y_i -\alpha - (x_i - x)\beta) K\left(\frac{x_i -
      x}{h}\right)
\end{align}
The first order condition for this problem is 
\begin{align}
  0 = & \frac{1}{n} \sum_{i=1}^n \left(\tau - \mathbf{1}(y_i -
    \alpha - (x_i - x) \beta) \right) K\left(\frac{x_i - x}{h}\right) 
\end{align}
Under appropriate conditions (the kernel has certain properties and a
LLN applies)
\begin{align}
  \frac{1}{n} \sum_{i=1}^n \left(\tau - \mathbf{1}(y_i -
    \alpha - (x_i - x) \beta) \right) K\left(\frac{x_i - x}{h}\right)  
  = & E\left[ \mathbf{1}(y_i - \alpha - (x_i - x) \beta)
    K\left(\frac{x_i - x}{h}\right) \right] + O(1/(nh)) 
\end{align}

\subsection{Consistency}
Let's verify Fischer consistency by showing that this expection is
zero at $\alpha = q(x)$.  This calculation will also give us the bias
as a function of $h$.  As in previous kernel situations, the key will
be write the expectation as an integral, make a change of variables,
and then take a Taylor expansion.  
\begin{align}
  E\left[ \mathbf{1}(y_i - \alpha - (x_i - x) \beta)
    K\left(\frac{x_i - x}{h}\right) \right] = & \int_X
  F_{y|x}\left(\alpha + (x_i - x) \beta\right)  K\left(\frac{x_i
      - x}{h}\right) f_x(x_i) dx_i \\
  & \text{change of variables $x_i = x+hu_i$} \notag \\
  = & h \int_U F_{y|x}\left(\alpha + h u_i \beta\right)  K(u_i)
  f_x(x+hu_i) du_i \\
  & \text{let $y_i = q(x_i) + \epsilon_i$} \notag \\
  = & h \int_U F_{\epsilon}\left(\alpha + h u_i \beta - q(x+hu_i)\right)  K(u_i)
  f_x(x+hu_i) du_i \\ 
  & \text{expand around $hu=0$} \\
  = & h \int_U [\tau - F_{\epsilon}(\alpha - q(x)) - hu_i
    f_\epsilon(\alpha - q(x))(\beta - q'(x)) - \notag \\
    & \,\,\,\,\, - h^2 \frac{1}{2} u_i\left(
      f_\epsilon'(\alpha - q(x))(\beta - q'(x)) q'(x)' + \nabla^2 q(x)
      f_\epsilon(\alpha-q(x))\right)u_i' + o(h^2u_iu_i') ] \times
  \notag \\
  & \,\,\,\,\,\, \times \left[f_x(x) + hu_i f_x'(x) + \frac{1}{2} h^2 u_i f_x"(x)
    u_i' + o(h^2u_iu_i') \right] K(u) du \\
  & \text{kernel has mean 0} \\
  = & h (\tau  - F_\epsilon(\alpha - q(x))f_x(x)) + \notag \\ 
  & + h^3 \int u'u K(u) du [ (\tau  - F_\epsilon(\alpha -
    q(x)))\frac{1}{2} tr(f_x"(x))
    - \notag \\ 
    & \,\,\,\,\,\, - f_x(x) \frac{1}{2}
    tr\left(f_\epsilon'(\alpha-q(x))(\beta - q'(x)) q'(x)' 
      + \nabla^2 q(x) f_\epsilon(\alpha - q(x)) \right) + (\beta -
    q'(x)) f_x'(x) f_\epsilon(\alpha - q(x)) ] + o(h^3) \label{focBias}
\end{align}
If we evaluate this at $\alpha_0 = q(x)$ and $\beta_0 = q'(x)$, we
have:
\begin{align}
  h^3 \frac{1}{2} \left(\int u'u K(u) du \right) f_x(x) tr(\nabla^2 q(x))
  f_\epsilon(0) \rightarrow 0 \text{ as } h \rightarrow 0
\end{align}
Importantly, this implies both $\hat{\beta}$ and $\hat{\alpha}$ are
consistent.  

\subsection{Bias}
To find the bias of $\hat{\alpha}$ we simply expand the first order
condition around $\alpha_0$.  Since the
first order condition is not differentiable, we expand the
approximation to the FOC.  Let $\widehat{FOC}$ denote the true first
order condition, and $FOC$ denote the asymptotic approximation.  Then
\begin{align}
  0 = & \widehat{FOC}(\hat{\alpha},\hat{\beta}) \notag \\ 
  = & FOC(\hat{\alpha},\hat{\beta}) + O_p(1/(nh)) \notag \\
  O_p(1/(nh)) = & FOC(\alpha_0,\beta_0) + (\hat{\alpha} - \alpha_0) \frac{\partial
    FOC}{\partial \alpha} + (\hat{\beta} - \beta_0) \frac{\partial
    FOC}{\partial \beta} 
\end{align}
From the above we know that $FOC(\alpha_0,\beta_0) = h^3 \frac{1}{2} \left(\int
  u'u K(u) du \right) f_x(x) tr(\nabla^2 q(x)) f_\epsilon(0) + o(h^3)$, and $\hat{\beta} - \beta_0 = O_p(1/(nh))$.  Also, by
inspecting, (ref{focBias}), $\frac{\partial
  FOC}{\partial \beta} = o(h^3)$.  From the same equation, we have:
\begin{align}
  \frac{\partial
  FOC}{\partial \alpha} = & h \left(f_\epsilon(0)f_x(x) + h^2
  \text{stuff} \right)
\end{align}
Thus,
\begin{align}
  (\hat{\alpha} - \alpha_0) = & \frac{h^2 \frac{1}{2} \left(\int
  u'u K(u) du \right) f_x(x) tr(\nabla^2 q(x))
f_\epsilon(0)}{f_\epsilon(0) f_x(x) + o(h^2)} + O_p(1/(nh))
\frac{o(h^3)}{o(h^3)} + O_p(1/(nh)) \\
= & h^2 \frac{1}{2} tr(\nabla^2 q(x)) \left(\int
  u'u K(u) du \right) + O_p(1/(nh)) 
\end{align}

\subsection{Variance}

To find the asymptotic variance, we will need to assume that
\begin{align}
  \widehat{FOC}(\hat{\alpha},\hat{\beta}) -
  \widehat{FOC}(\alpha_0,\beta_0) = FOC(\hat{\alpha},\hat{\beta}) -
  FOC(\alpha_0,\beta_0) + o_p(1/\sqrt{nh^d}) 
\end{align}
This would follow from stochastic equicontinuity of the first order
condition.  The theory section of problem set 3 showed that it holds,
but the argument is long, so I won't repeat it here.  With this
assumption, we have:
\begin{align} 
  0 = & \widehat{FOC}(\hat{\alpha},\hat{\beta}) \\
  = & \widehat{FOC}(\alpha_0,\beta_0) + FOC(\hat{\alpha},\hat{\beta}) -
  FOC(\alpha_0,\beta_0) o_p(1/\sqrt{nh^d}) \\
  = & \widehat{FOC}(\alpha_0,\beta_0) + (\hat{\alpha} - \alpha_0) \frac{\partial
    FOC}{\partial \alpha} + (\hat{\beta} - \beta_0) \frac{\partial
    FOC}{\partial \beta} + o(h^3) + o_p(1/\sqrt{nh^d}) \\
  = & \widehat{FOC}(\alpha_0,\beta_0) + (\hat{\alpha} - \alpha_0)
  (h f_x(x) f_\epsilon(0) + o(h^3) + O_p(1/(nh))o(h^3) + o(h^3) +
  o_p(1/\sqrt{nh^d}) \\
  \sqrt{nh^d} (\hat{\alpha} - \alpha_0) = & \frac{\sqrt{nh^d}\widehat{FOC}(\alpha_0,\beta_0)}{h f_x(x)
    f_\epsilon(0) } + \sqrt{nh^d} o(h)^3 + o_p(1)
\end{align}
Now, we just need to find the distribution of
$\sqrt{nh^d}\widehat{FOC}(\alpha_0,\beta_0)$.  That is,
\begin{align}
  \sqrt{nh^d}\widehat{FOC}(\alpha_0,\beta_0) = &
  \sqrt{\frac{h^d}{n}} \sum_{i=1}^n \left(\tau - \mathbf{1}(\epsilon_i -
    (x_i - x) q'(x) < 0)\right) K\left(\frac{x_i - x}{h} \right) 
\end{align}
Assume a CLT applies so that this is normal with mean given by its
exectation (which we computed above) and variance equal to 
\begin{align}
  var\left[\left(\tau - \mathbf{1}(\epsilon_i -
      (x_i - x) q'(x) < 0)\right) K\left(\frac{x_i - x}{h}
    \right)\right] = & var(f(x_i,\epsilon_i))  \notag \\
  = & E[var(f(x_i,\epsilon_i)|x_i)] + var[E[f_x,\epsilon_i]|x_i]
  \notag \\
  \\ & \text{steps like those taken to show consistency} \notag \\
  = & h^2 \tau(1-\tau) f_x(x) \int K(u)^2 du  + 0
\end{align}
Thus, we have:
\begin{align}
  \sqrt{nh^d} (\hat{\alpha} - \alpha_0) \indist N\left(h^2 \frac{1}{2} tr(\nabla^2 q(x)) \left(\int
  u'u K(u) du \right),\frac{\tau(1-\tau) \int K(u)^2}{f_x(x)
  f_\epsilon(0)^2} \right)
\end{align}
and the MSE is
\begin{align}
  MSE(\hat{\alpha} - \alpha_0) = & h^4 \frac{1}{4} tr(\nabla^2 q(x))^2 \left(\int
  u'u K(u) du \right)^2  + \frac{\tau(1-\tau) \int K(u)^2}{nh^d f_x(x)
  f_\epsilon(0)^2}
\end{align}
With $d=1$, The MSE minimizing bandwidth will then be $h = O(n^{-1/5}) $, which
makes the bias squared and variance disappear at the same
rate.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Semi-Parametric Estimators}

\subsection{Transformation Model}

A transformation model takes the form:
\begin{align}
  \Lambda(y)  = x\beta + u
\end{align}
where the transformation, $\Lambda$, and the distribution of $u$ are
unknown.  Recall from problem set 3 and recitation 6 that the
transformation model is generalization of the mixed proportional
hazard model and the accelerated failure time model.  In particular,
the mixed proportional hazard model can be written as
\begin{align}
  \log z(t) = x\beta - \log(\alpha) + u 
\end{align}
where $z(t)$ is the integrated hazard function, $\alpha$ is frailty,
and $u$ has a Gumbel distribution.  

\subsubsection{Estimation of $\beta$}
There are multiple $\sqrt{n}$ consistent estimators of $\beta$.  Han
(1987) describes a maximum rank score estimator of $\beta$.  $\beta$
can also be recovered from an estimator of average derivative of
$E[y|x]$.  We discussed average derivative estimation in class, so
let's take that approach.  We can write the transformation model as:
\begin{align}
  y = \Lambda^{-1} (x\beta + u) 
\end{align}
so
\begin{align}
  E[y|x] = E[\Lambda^{-1} (x\beta + u)|x] = g(x)
\end{align}
The derivative is
\begin{align}
  \frac{\partial g}{\partial x_k} = E[ {\Lambda^{-1}}'(x\beta+u)|x] \beta_k
\end{align}
so we can recover $\beta$ up to scale by looking at
\begin{align}
  \frac{E[w(x)(\partial g)/(\partial x_k)]}{E[w(x)(\partial g)/(\partial x_j)]} =
  \frac{\beta_k}{\beta_j} 
\end{align}
where $w(x)$ is any function that depends on $x$.  We could set $w(x)
= 1$, and then use the estimator of $E[(\partial g)/(\partial
  x_k)]$ discussed in lecture.  Alternatively, we can follow Powell,
Stock, and Stoker (1989) and take $w(x) = f(x)$.  Then, integration by
parts shows that
\begin{align}
  E[f(x) (\partial g)/(\partial x_k) ] = -2 E[y f'(x) ] 
\end{align}
The derivative of the density of $x$ can be estimated by kernel
methods,
\begin{align}
  \hat{f'}_j (x) = \frac{1}{n} \sum_{i \neq j} \frac{1}{h}
  K'\left(\frac{x_i-x}{h}\right) 
\end{align}
where we leave out observation $j$ to simplify analysis of  the
weighted average derivative estimator, which is
\begin{align}
\  \widehat{E[f(x) (\partial g)/(\partial x_k) ]} = \frac{-2}{n} \sum y_i
    \hat{f'}_i(x_i) 
\end{align}
Powell, Stock, and Stoker show that if $f$ is $\frac{d}{2} + 3$ times
differentiable, where $d$ is the dimension of $x$, and $Nh^{d+2}
\rightarrow \infty$, while $Nh^{d + 4}  \rightarrow 0$ then 
\begin{align}
    \sqrt{n} (\widehat{E[f(x)(\partial g)/(\partial x_k) ]} -
    E[f(x)(\partial g)/(\partial x_k) ]) \indist N(0,V) 
\end{align}
This implies that our estimate of $\beta$ is also asymptotically
normal and converges at a $\sqrt{n}$ rate
\begin{align}
 \sqrt{n}(  \widehat{\frac{\beta_k}{\beta_j}} -
 \frac{\beta_k}{\beta_j}) \indist N(0,V_\beta)
\end{align}

\subsubsection{Estimation of $\Lambda(\cdot)$}
Horowitz (1996) developed an estimator of $\Lambda(\cdot)$ and
$F_u(\cdot)$.  Somewhat surprisingly, both of these functions can be
estimated at a $\sqrt{n}$ rate.  Consider the cdf of $y$ conditional
on $x\beta$:
\begin{align}
  G(y|x\beta) = & F_u(\Lambda(y) - x\beta) 
\end{align}
Differentiating $G(y|z)$ with respect to $y$ and $z$ gives:
\begin{align}
  G_{y}(y|z) = \Lambda'(y) f_u(\Lambda(y) - x\beta) \\
  G_z(y|z) = f_u(\Lambda(y) - x\beta) 
\end{align}
So,
\begin{align}
  \Lambda(y) = & \Lambda(y_0) + \int_{y_0}^y \int_Z w(z) dz \Lambda'(y) dy \\
  = & \Lambda(y_0) + \int_{y_0}^y \int_Z w(z)
  \frac{G_{y}(y|z)}{G_z(y|z)} dz dy \label{2.5} \\
\end{align}
where $w(z)$ is any function that integrates to $1$.  $G_{y}$ and
$G_z$.  An estimate of $\Lambda$ can be constructed as
\begin{align}
  \hat{\Lambda}(y) = \int_{y_0}^y \int_Z w(z)
  \frac{\hat{G}_{y}(y|z)}{\hat{G}_z(y|z)} dz dy \label{lest}
\end{align}
Horowitz proposes estimating $\hat{G}_{y}$ by taking the derivative
with respect to $z$ of the kernel regression of $\mathbf{1}(y_i < y)$
on $z_i = x_i\beta$.  Local linear or polynomial regression could
probably also
be used to estimate the derivative directly.  $\hat{G}_{y}(y|z)$ can
be estimated using kernel methods as well.  Specifically,
\begin{align}
  \hat{G}_{y}(y|z) = \frac{ \sum K_{y}((y_i - y)/h) K_z((z_i
    - z)/h)} {\sum K_z((z_i -z)/h) }
\end{align}
Although the estimates of both $\hat{G}_{y}$ and $\hat{G}_z$ converge
slower than $\sqrt{n}$, the estimate of $\hat{\Lambda}$ averages over
many values of $z$ to achieve a $\sqrt{n}$ rate of convergence.  See
Horowitz's paper for details.  Important conditions are that $K_{z}$
be a higher-order kernel (ie $\int u^k K_{z}(u) du = 0$ for some
$k>1$) and that $h_y$ and $h_z$ go to zero faster than the usual
optimal rate.  Both of these extra conditions are to make the bias of
$\hat{\Lambda}$ be $o_p(1/\sqrt{n})$.  

\subsection{Estimation of $F_u(\cdot)$}
The obvious estimator of $F_u$ is the empirical distribution function
of $\hat{\Lambda}(y) - x\hat{\beta}$.  However, this estimator will
not work because $\hat{\Lambda}(y)$ is only consistently estimable in
a compact interval, $[y_1, y_2]$.  Nonetheless, a simply modification
of the EDF works.  This estimator is based on $F_u(u) P(y \in
[y_1,y_2]|u) = F_u(u|y \in [y_1,y_2]) P(y \in [y_1,y_2])$, which
rearranged suggests: 
\begin{align}
  \hat{F}_u(u) = \frac{\sum \mathbf{1}(\hat{\Lambda}(y_i) - x_i
    \hat{\beta} < u) \mathbf{1}(\hat{\Lambda}(y_1) - u \leq x_i \beta
    \leq \hat{\Lambda}(y_2) - u)}{\sum \mathbf{1}(\hat{\Lambda}(y_1) -
    u \leq x_i \beta \leq \hat{\Lambda}(y_2) - u)} 
\end{align}

\end{document}
