ols.Rmd

# Ordinary Least Squares
Let $Y$ be our outcome random variable and 
$$
\v{X}=\begin{bmatrix}
1 & X_{[1]} & X_{[2]} & \ldots & X_{[k]}
\end{bmatrix}
$$ 
be our predictor (or explanatory) vector containing $k$ predictors and a constant. We denote the joint distribution of $(Y,\v{X})$ by $F(y,\v{x})$, i.e.,
$$
F(y,\v{x})=\p{Y\leq y, \v{X}\leq\v{x}}
=\p{Y\leq y,X_1\leq x_1,\ldots,X_k\leq x_k}.
$$

The **dataset** or **sample** is a collection of observations $\{(Y_i,\v{X}_i): i=1,2,\ldots,n\}$. We assume that each observation
$(Y_i,\v{X}_i)$ is a random (row) vector drawn from the common distribution, sometimes referred to as the **population**, $F$.

For a given vector of (unknown) coefficients 
$\v{\beta}=\begin{bmatrix}\beta_0 & \beta_1 & \ldots & \beta_k\end{bmatrix}^T\in\mathbb{R}^{k+1}$, we define the following **cost function**:
$$
\widehat{S}(\v{\beta})=\frac{1}{n}\sum\limits_{i=1}^n(Y_i-\v{X_i}\v{\beta})^2.
$$
The cost function $\widehat{S}({\v{\beta}})$ can also be thought of as the average sum of residuals. In fact, $\widehat{S}({\v{\beta}})$ is the moment (plug-in) estimator of the mean squared error, 
$$
S(\v{\beta})=\E{(Y-\v{X}\v{\beta})^2}.
$$

We now minimize $\widehat{S}({\v{\beta}})$ over all possible choices of
$\v{\beta}\in\mathbb{R}^{k+1}$. When the minimizer exists and is unique, we call it the **least squares estimator**, denoted $\widehat{\v{\beta}}$. 

:::{.definition name="(Ordinary) Least Squares Estimator"}
The least square estimator is $$
\widehat{\v{\beta}}
=\underset{\v{\beta}\in\mathbb{R}^{k+1}}{\arg\min} \ \widehat{S}(\v{\beta}),
$$
provided it exists uniquely.
:::

## Solution of OLS
We rewrite the cost function as
$$
\widehat{S}(\v{\beta})=\frac{1}{n}SSE(\v{\beta}),
$$
where 
$SSE(\v{\beta})\stackrel{\text{def}}{=}\sum\limits_{i=1}^n(Y_i-\v{X_i}\v{\beta})^2$.

We now express $SSE(\v{\beta})$ as a quadratic function of $\v{\beta}$.
$$
\begin{aligned}
SSE &=\sum\limits_{i=1}^n(Y_i-\v{X_i}\v{\beta})^2 \\
&=\sum\limits_{i=1}^n Y_i^2 
  - 2\sum\limits_{i=1}^n Y_i(\v{X_i}\v{\beta})
  + \sum\limits_{i=1}^n (\v{X_i}\v{\beta})^2 \\
&=\sum\limits_{i=1}^n Y_i^2 
  - 2\sum\limits_{i=1}^n Y_i(\v{\beta}^T\v{X_i}^T)
  + \sum\limits_{i=1}^n (\v{X_i}\v{\beta})(\v{X_i}\v{\beta}) \\
&=\sum\limits_{i=1}^n Y_i^2 
  - 2\sum\limits_{i=1}^n \v{\beta}^T(Y_i\v{X_i}^T)
  + \sum\limits_{i=1}^n (\v{\beta}^T\v{X_i}^T)(\v{X_i}\v{\beta}) \\
&=\left(\sum\limits_{i=1}^n Y_i^2\right) 
  - 2\v{\beta}^T\left(\sum\limits_{i=1}^n\v{X_i}^TY_i\right)
  + \v{\beta}^T\left(\sum\limits_{i=1}^n \v{X_i}^T\v{X_i}\right)\v{\beta}
\end{aligned}
$$
Taking partial derivative w.r.t. $\beta_j$, we get
$$
\frac{\partial}{\partial\beta_j}SSE(\v{\beta})=-2\left[\sum\limits_{i=1}^n\v{X_i}^TY_i\right]_j 
+ 2\left[\left(\sum\limits_{i=1}^n  \v{X_i}^T\v{X_i}\right)\v{\beta}\right]_j.
$$

Therefore,
$$
\frac{\partial}{\partial\v{\beta}}SSE(\v{\beta})
=-2\left(\sum\limits_{i=1}^n\v{X_i}^TY_i\right) 
+ 2\left(\sum\limits_{i=1}^n  \v{X_i}^T\v{X_i}\right)\v{\beta}.
$$

In order to miniminize $SSE(\v{\beta})$, a necessary condition for $\widehat{\v{\beta}}$ is
$$
\frac{\partial}{\partial\v{\beta}}SSE(\v{\beta})\bigg|_{\v{\beta}
=\widehat{\v{\beta}}}=\v{0},
$$
i.e.,
$$
-2\left(\sum\limits_{i=1}^n\v{X_i}^TY_i\right) 
+ 2\left(\sum\limits_{i=1}^n  \v{X_i}^T\v{X_i}\right)\widehat{\v{\beta}}
=\v{0}
$$
So,
\begin{equation}
\left(\sum\limits_{i=1}^n\v{X_i}^TY_i\right)
=\left(\sum\limits_{i=1}^n  \v{X_i}^T\v{X_i}\right)\widehat{\v{\beta}}
(\#eq:moment-0)
\end{equation}

Both the left and right hand side of the above equation are $k+1$ vectors. So, we have a system of $(k+1)$ linear equations with $(k+1)$ unknowns---the elements of $\v{\beta}$. 

Let us define

$$
\widehat{\mathbb{Q}}_{\v{XX}}
=\frac{1}{n}\left(\sum\limits_{i=1}^n\v{X_i}^T\v{X_i}\right)
\mbox{ and }
\widehat{\mathbb{Q}}_{\v{X}Y}
=\frac{1}{n}\left(\sum\limits_{i=1}^n\v{X_i}^TY_i\right).
$$

Rewriting \@ref(eq:moment-0), we get
\begin{equation}
\widehat{\mathbb{Q}}_{\v{X}Y}=\widehat{\mathbb{Q}}_{\v{XX}}
\widehat{\v{\beta}}.
(\#eq:moment)
\end{equation}

Equation \@ref(eq:moment) is sometimes referred to as the first-order moment condition. For the uniqueness of solution, we require that
$\widehat{\mathbb{Q}}_{\v{XX}}$ is non-singular. In that case, we can solve for $\widehat{\v{\beta}}$ to get,
$$
\widehat{\v{\beta}}=\left[\widehat{\mathbb{Q}}_{\v{XX}}\right]^{-1}
\widehat{\mathbb{Q}}_{\v{X}Y}.
$$
To verify that the above choice minimizes $SSE(\v{\beta})$,
one can consider the second-order moment conditions.
$$
\frac{\partial^2}{\partial\v{\beta}\partial\v{\beta}^T}SSE(\v{\beta})
=2\widehat{\mathbb{Q}}_{\v{XX}}.
$$
If $\widehat{\mathbb{Q}}_{\v{XX}}$ is non-singular, it is also positive-definite. So, we have actually proved the following theorem.

:::{.theorem}
If \ $\widehat{\mathbb{Q}}_{\v{XX}}$ is non-singular, then the
least squares estimator is unique, and is given by
$$
\widehat{\v{\beta}}=\left[\widehat{\mathbb{Q}}_{\v{XX}}\right]^{-1}
\widehat{\mathbb{Q}}_{\v{X}Y}.
$$
:::


## Errors and Residuals
Recall that $\v{\beta}$ denotes the coefficients of the best linear predictor \@ref(thm:blp). We first define the **fitted** value as
$$
\widehat{Y}_i=\v{X}_i\widehat{\v{\beta}}\mbox{ for }
i=1,2,\ldots,n.
$$
For the least squares estimators, we define the **errors** and **residuals** in the following way:
$$
\eps_i=Y_i-\v{X}_i\v{\beta}, \mbox{ and } 
e_i=Y_i-\widehat{Y}_i.
$$

:::{.theorem #olserror name="Least Squares Error"}
If \ $\widehat{\mathbb{Q}}_{\v{XX}}$ is non-singular, then  
1. $\sum\limits_{i=1}^n\v{X}_i^Te_i=\v{0}$  
2. $\sum\limits_{i=1}^ne_i=0$
:::

:::{.proof}
$$
\begin{aligned}
\sum\limits_{i=1}^n\v{X}_i^Te_i 
&=\sum\limits_{i=1}^n\v{X}_i^T(Y_i-\widehat{Y}_i) \\
&=\sum\limits_{i=1}^n\v{X}_i^TY_i-\sum\limits_{i=1}^n\v{X}_i^T\widehat{Y}_i \\
&=\sum\limits_{i=1}^n\v{X}_i^TY_i-\sum\limits_{i=1}^n\v{X}_i^T\v{X}_i\v{\widehat{\beta}} \\
&=\widehat{\m{Q}}_{\v{X}Y}-\widehat{\m{Q}}_{\v{XX}}\v{\widehat{\beta}} \\
&=\widehat{\m{Q}}_{\v{X}Y}-\widehat{\m{Q}}_{\v{XX}}
\left( \widehat{\m{Q}}_{\v{XX}}^{-1} \widehat{\m{Q}}_{\v{X}Y} \right) \\
&=\v{0}
\end{aligned}
$$
From the first row of (1) we get

$$
\sum\limits_{i=1}^n e_i=0.
$$
Hence the result.
:::


## Matrix Notation
Taking the definition of errors from the last section, we can write down a system of $n$ linear equations:
$$
\begin{aligned}
Y_1 &= \v{X_1}\v{\beta} + \eps_1 \\
Y_2 &= \v{X_2}\v{\beta} + \eps_2 \\
& \vdots \\
Y_n &= \v{X_n}\v{\beta} + \eps_n
\end{aligned}
$$

Define 
$$
\v{Y}=\begin{bmatrix}
Y_1 \\
Y_2 \\
\vdots \\
Y_n
\end{bmatrix}_{n\times1},\ 
\mathbb{X}=\begin{bmatrix}
1 & X_{[1]1} & X_{[2]1} & \ldots & X_{[k]1} \\
1 & X_{[1]2} & X_{[2]2} & \ldots & X_{[k]2} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & X_{[1]n} & X_{[2]n} & \ldots & X_{[k]n}
\end{bmatrix},\mbox{ and }
\v{\eps}=\begin{bmatrix}
\eps_1 \\
\eps_2 \\
\vdots \\
\eps_n
\end{bmatrix}_{n\times1}.
$$
We can now rewrite the system as the following:
$$
\v{Y}=\m{X}\v{\beta}+\v{\eps}.
$$
We also note that
$$
\widehat{\m{Q}}_{\v{XX}}=\sum\limits_{i=1}^n\v{X}^T_i\v{X}_i=
\m{X}^T\m{X},
$$
and
$$
\widehat{\m{Q}}_{\v{X}Y}=\sum\limits_{i=1}^n\v{X}_i^TY_i=
\m{X}^T\v{Y}.
$$
So, we have write the least squares estimator as
$$
\widehat{\v{\beta}}=\left[\m{X}^T\m{X}\right]^{-1}\m{X}^T\v{Y}.
$$
Similarly, the residual vector is
$$
\v{e}=\v{Y}-\m{X}\widehat{\v{\beta}}.
$$
As a consequence of \@ref(thm:olserror), we can write
$$
\m{X}^T\v{e}=\v{0}.
$$


## An Example
Let $X$ and $Y$ be random variables with $X\sim U(0,1)$ and $Y=10X^2+W$, 
where $W\sim N(0,1)$ and $X\ind W$. 

- Derive the CEF $E[Y|X]$.
- How to compute the BLP?
- Why do we care about the BLP?
\newpage

## Comparing the CEF, BLP, and Regression
As it turns out, the BLP is $Y=-\frac{5}{3}+10X$.
```{r}
CEF <- function(X) {  10 * X^2 }
BLP <- function(X) {  - 5 / 3 + 10 * X}
```

Let's plot both the CEF and BLP.
```{r}
plot(NA,type="n", xlab="X", ylab="Y", xlim=c(0, 1), ylim=c(-2, 10), 
     main="Plot CEF and BLP")
lines(c(0:100)/100, CEF(c(0:100)/100), col="black", lwd=2)
lines(c(0:100)/100, BLP(c(0:100)/100), col="red", lwd=2)
```

## Draw a Random Sample and Plot Them All!
```{r}
x <- runif(50, -1, 1)
x <- sort(x)
w <- rnorm(50, 0, 1)
d <- data.frame( x,  y = 10 * x^2 + w)

plot(d$x, d$y, type = "p", pch = 1, 
     main = "Scatter plot", xlab = "x", ylab = "y")
 lines(c(-100:100)/100, CEF(c(-100:100)/100), col="black", lwd=0.5)
 #lines(c(-100:100)/100, BLP(c(-100:100)/100), col="red", lwd=0.5, lty = 2)
 abline(lm(y ~ x, data = d), col = "red", lwd = 2)
 lines(x, predict(lm(y~I(x^2),data =d)), col='green', lwd=2)
```