Estimation with method of maximum likelihood Can anybody help me to generate the estimator of equation:
$$Y_i = \beta_0 + \beta_1X_{i1} + \beta_2X_{i2}+\cdots+\beta_4X_{i4}+\varepsilon_i$$
using method of maximum likelihood, where $\varepsilon_i$ are independent variables which have normal distribution $N(0,\sigma^2)$

This is given by least squares estimation. To see this, write
$$
L(\beta, \sigma^2 | Y) = \prod_i (2\pi\sigma^2)^{1/2} \exp\left(\frac {-1} {2\sigma^2} (Y_i - \beta_0 - \sum_j \beta_j X_{ij})^2\right) = (2\pi\sigma^2)^{n/2} \exp\left(\frac {-1} {2\sigma^2} \sum_i(Y_i - \beta_0 - \sum_j \beta_j X_{ij})^2\right)
$$
Maximization of this is equivalent to minimization of $\sum_i (Y_i - \beta_0 - \sum_ j\beta_j X_{ij})^2$, which is the definition of least squares estimation. If $\mathbf X$ is such that $\langle \mathbf X\rangle_{ij} = X_{ij}$ and $\mathbf Y$ is such that $\langle\mathbf Y\rangle_{ij} = Y_{ij}$ then we can rewrite this as the minimization of $Q(\beta) = \|\mathbf Y - \mathbf X\beta\|^2$ with respect to $\beta$, and there are several ways to see that $\hat \beta = (\mathbf X^T \mathbf X)^{-1} \mathbf X^T \mathbf Y$ is the minimizer. The most expedient is probably to calculate $\partial Q/ \partial \beta$ and set equal to $0$ which gives
$$
-2 \mathbf X^T \mathbf Y + 2 (\mathbf X^T \mathbf X) \beta = 0
$$
and note that $Q$ is a strictly convex function of $\beta$ which ensures that the solution to this equation is a minimum of $Q$.
Note however that we are relying on $\mathbf X^T \mathbf X$ being invertible, i.e. $\mathbf X$ has linearly independent columns; if this fails, $Q$ is not strictly convex and any generalized inverse $(\mathbf X^T\mathbf X)^-$ can be used to replace $(\mathbf X^T \mathbf X)^{-1}$. 
