How to find the closed form formula for $\hat{\beta}$ while using ordinary least squares estimation? According to Wikipedia's article on Linear Regression:

Given a data set $\{y_i,x_{i1},\ldots,x_{ip}\}_{i=1}^{i=n}$ of $n$
  statistical units, a regression model assumes that the relationship
  between the dependent variable $y_i$ and the $p-\text{vector}$ or
  regressors $x_i$ is linear. This relationship is modelled through a
  disturbance term or error variable $\varepsilon_i$-an unobserved random
  variable that adds random noise to the to the linear relationship
  between the dependent variables and regressors. The model takes the
  form
  $y_i=\beta_0(1)+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\epsilon_i=x_i^T \beta + \varepsilon_i$

These equations can be written in vector form as $$y=\mathbf{X\beta+\epsilon}$$
For the Ordinary Least Square estimation they say that the closed form expression for the estimated value of the unknown parameter $\beta$ is 
$$\hat{\mathbf{\beta}}=(\mathbf{X^{T}X})^{-1}\mathbf{X}^{T}y$$ 
I'm not sure how they get this formula for $\hat{\beta}$. It would be very nice if someone can explain me the derivation.
 A: I'm going to show this using partial differentiation.
Consider the assumed linear model
$$y_i = \mathbf{x}_i^{T}\boldsymbol\beta + \epsilon_i$$
where $y_i, \epsilon_i \in \mathbb{R}$ and $\mathbf{x}_i=\begin{bmatrix}
x_{i0} \\
x_{i1} \\
\vdots \\
x_{ip}
\end{bmatrix}, \boldsymbol\beta = \begin{bmatrix}
\beta_0 \\
\beta_1 \\
\vdots \\
\beta_p
\end{bmatrix} \in \mathbb{R}^{p+1}$ for $i = 1, \dots, n$, with $x_{i0} = 1$.
Our aim is to solve for $\hat{\boldsymbol\beta}$ by minimizing the residual sum of squares, or minimizing $$\text{RSS}(\boldsymbol\beta) = \sum_{i=1}^{n}(y_i-\mathbf{x}_i^{T}\boldsymbol\beta)^2\text{.}$$
To compute this sum, consider the vector of residuals
$$\mathbf{e}=\begin{bmatrix}
y_1 - \mathbf{x}_1^{T}\boldsymbol\beta \\
y_2 - \mathbf{x}_2^{T}\boldsymbol\beta \\
\vdots \\
y_n - \mathbf{x}_n^{T}\boldsymbol\beta
\end{bmatrix}$$
Then $\text{RSS}(\boldsymbol\beta) = \mathbf{e}^{T}\mathbf{e}$. Our next step is to find the "partial derivatives" of $\text{RSS}(\boldsymbol\beta)$. 
To do this, note that for $k = 1, \dots, p$,
$$\dfrac{\partial \text{RSS}}{\partial \beta_k}=\dfrac{\partial}{\partial\beta_k}\left\{\sum_{i=1}^{n}\left[y_i- \sum_{j=0}^{p}\beta_jx_{ij}\right]^2 \right\}=-2\sum_{i=1}^{n}x_{ik}\left(y_i - \sum_{j=0}^{p}\beta_jx_{ij}\right)\text{.}$$
"Stacking" these, we obtain
$$\begin{align}
\dfrac{\partial \text{RSS}}{\partial \boldsymbol\beta}&=\begin{bmatrix}
\dfrac{\partial \text{RSS}}{\partial \beta_0} \\
\dfrac{\partial \text{RSS}}{\partial \beta_1} \\
\vdots \\
\dfrac{\partial \text{RSS}}{\partial \beta_p}
\end{bmatrix} \\
&= \begin{bmatrix}
-2\sum_{i=1}^{n}x_{i0}\left(y_i - \sum_{j=0}^{p}\beta_jx_{ij}\right) \\
-2\sum_{i=1}^{n}x_{i1}\left(y_i - \sum_{j=0}^{p}\beta_jx_{ij}\right) \\
\vdots \\
-2\sum_{i=1}^{n}x_{ip}\left(y_i - \sum_{j=0}^{p}\beta_jx_{ij}\right)
\end{bmatrix} \\
&= -2\begin{bmatrix}
\sum_{i=0}^{n}x_{i0}(\mathbf{y}-\mathbf{x}_1^{T}\boldsymbol\beta)\\
\sum_{i=0}^{n}x_{i1}(\mathbf{y}-\mathbf{x}_1^{T}\boldsymbol\beta) \\
\vdots \\
\sum_{i=0}^{n}x_{ip}(\mathbf{y}-\mathbf{x}_1^{T}\boldsymbol\beta)
\end{bmatrix} \\
&= -2\left(\begin{bmatrix}
\sum_{i=0}^{n}x_{i0}\mathbf{y}\\
\sum_{i=0}^{n}x_{i1}\mathbf{y} \\
\vdots \\
\sum_{i=0}^{n}x_{ip}\mathbf{y}
\end{bmatrix} - \begin{bmatrix}
\sum_{i=0}^{n}x_{i0}\mathbf{x}_1^{T}\boldsymbol\beta)\\
\sum_{i=0}^{n}x_{i1}\mathbf{x}_1^{T}\boldsymbol\beta) \\
\vdots \\
\sum_{i=0}^{n}x_{ip}\mathbf{x}_1^{T}\boldsymbol\beta)
\end{bmatrix}\right)\\
&= -2(\mathbf{X}^{T}\mathbf{y}-\mathbf{X}^{T}\mathbf{X}\boldsymbol\beta)\text{.}
\end{align}$$
where $$\mathbf{X} = \begin{bmatrix}
\mathbf{x}_1^{T} \\
\mathbf{x}_2^{T} \\
\vdots \\
\mathbf{x}_n^{T}
\end{bmatrix}\text{.}$$
Setting $\dfrac{\partial \text{RSS}}{\partial \boldsymbol\beta} = \mathbf{0}$, we obtain $$\mathbf{X}^{T}\mathbf{X}\boldsymbol\beta = \mathbf{X}^{T}\mathbf{y}$$
and assuming $\mathbf{X}^{T}\mathbf{X}$ is invertible, 
$$\hat{\boldsymbol\beta} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}$$
A: Proposition: The value of $\beta$ that minimizes $\|Y-X\beta\|^2$ is the one that makes $X\beta$ equal to the orthogonal projection onto the column space of $X.$
Proof: The desired value is called $\widehat\beta.$ Let $X\beta_P$ be the orthogonal projection of $Y$ onto the column space of $X.$ Then for any value of $\beta$ we have
\begin{align}
\|Y - X\widehat\beta_P\|^2 & = \|(Y-X\beta_P) + (X\beta_P - X\widehat\beta)\|^2 \\[10pt]
& = \|Y-X\beta_P\|^2 + (Y-X\beta_P)^T(X\beta_P-X\widehat\beta) + \| X(\beta_P-\widehat\beta) \|^2.
\end{align}
The middle term is $0$ because it's the dot product of something orthogonal to the column space and something in the column space. We therefore have
$$
\|Y-X\beta_P\|^2 + \|X(\beta_P-\widehat\beta)\|^2.
$$
The value of $\widehat\beta$ that minimizes this is the one that makes the makes the second term $0,$ and that is $\widehat\beta=\beta_P. \quad \blacksquare$
Proposition: Suppose $X$ has linearly independent columns (but not rows, since $X$ has many more rows than columns). Then the orthogonal projection onto the column space of $X$ is
$$
y \mapsto X(X^TX)^{-1}X^T y.
$$
Proof: Suppose first that $y$ is in the column space of $X.$ Then $y = X\alpha$ (for some $\alpha$). So
$$
X(X^TX)^{-1}X^Ty = X(X^TX)^{-1}X^T (X\alpha) = X\alpha = y.
$$
Next suppose $y$ is orthogonal to the column space of $X.$ Then $X^Ty=0$ so $X(X^TX)^{-1}X^Ty = 0. \qquad \blacksquare$
Therefore, for the desired value $\widehat\beta,$ we have
$$
X\widehat\beta = X(X^TX)^{-1}X^T y. \tag 1
$$
Since $X$ has more rows than columns, it does not have a inverse matrix, but since its columns are linearly independent, it has a left inverse (more than one of them, in fact) and one such left inverse is $(X^TX)^{-1}X^T,$ as you can see by multiplying. So multiply both sides of $(1)$ by that left inverse, and you've got it.
A: Essentially we're trying to minimize $\epsilon ^T \epsilon$ with respect to $\beta$. 
So we're trying to get 
\begin{align} 
\frac{ d \epsilon ^ T \epsilon }{ d \beta } &= 0\\
\\
\frac{ d ( y - X\beta ) ^ T ( y - X\beta) }{ d \beta } &= 0\\
\\
\frac{d}{d\beta}(-2y^T(X\beta) + \beta^TX^TX\beta) &= 0 \\
\\
2(X^TX)\beta &= 2 X^Ty \\
\\
\hat \beta &= (X^TX)^{-1} X^Ty
\end{align}
As required, assuming $X^TX$ is invertible.
