3.3 Estimation by maximum likelihood method ML
As in the simple linear regression model, we assume that the error are normally and identically distributed so that:
\[\begin{align} L&=f(\varepsilon_1,\varepsilon_2,..,\varepsilon_n)\notag\\ &=f(\varepsilon_1)\times f(\varepsilon_2)\times .. \times f(\varepsilon_n)\notag\\ &=\frac{1}{\sqrt{2\pi\sigma^2}}exp\bigg(-\frac{\varepsilon_1^2}{2\sigma^2}\bigg)\times\frac{1}{\sqrt{2\pi\sigma^2}}exp\bigg(-\frac{\varepsilon_2^2}{2\sigma^2}\bigg)\times .. \times \frac{1}{\sqrt{2\pi\sigma^2}}exp\bigg(-\frac{\varepsilon_n^2}{2\sigma^2}\bigg)\notag\\ &= \bigg(\frac{1}{\sqrt{2\pi\sigma^2}}\bigg)^n exp\bigg(-\frac{\sum\varepsilon_t^2}{2\sigma^2}\bigg)\notag\\ &=\bigg(\frac{1}{\sqrt{2\pi\sigma^2}}\bigg)^n exp\bigg(-\frac{\varepsilon^t\varepsilon}{2\sigma^2}\bigg) \tag{3.25} \end{align}\]
Since \(y\) is a function of \(\varepsilon\), the expression \(g(y)=\frac{d\varepsilon_t}{d y_t}f(\varepsilon_t)\) will be used to derive the distribution of \(y\) from that of \(\varepsilon\). so using (3.1), we get:
\[\begin{equation} y_t=\beta_1+\beta_2x_{2t}+..+\beta_kx_{kt}+\varepsilon_t\implies \frac{d\varepsilon_t}{d y_t}=1 \tag{3.26} \end{equation}\]
Substituting then this result into the above likelihood function (3.25):
\[\begin{equation} L=\bigg(\frac{1}{\sqrt{2\pi\sigma^2}}\bigg)^n exp\bigg(-\frac{(y-X\beta)^t(y-X\beta)}{2\sigma^2}\bigg) \tag{3.27} \end{equation}\]
As we see the maximization of this function, with respect to \(\beta\), reduces to minimizing the term \((y-X\beta)^t(y-X\beta)\), which is The OLS criterion. That is, the ML method leads to the OLS estimators. The advantage of this method is it can estimate, in addition to the regression coefficients, the variance of the errors \(\sigma^2\) at the same time.
For simplification, Instead of maximizing the likelihood function, we maximize its log since it is a monotonic function as follows:
\[\begin{equation*} l=Log(L)=-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^t(y-X\beta) \end{equation*}\]
For convex function, the extremum (maximum or minimum) will be located when the derivatives are equal to zero:
\[\begin{equation*} \begin{cases} \frac{\partial l}{\partial \beta}=\frac{1}{\sigma^2}X^t(y-X\beta_{ML})=0\\ \frac{\partial l}{\partial \sigma^2}=-\frac{n}{\sigma_{ML}^2}+\frac{1}{(\sigma_{ML}^2)^2}(y-X\beta_{ML})^t(y-X\beta_{ML})=0 \end{cases} \end{equation*}\]
Solving the first equation for \(\beta_{ML}\), we get the maximum likelihood estimator:
\[\begin{equation} \beta_{ML}=\big(X^tX\big)^{-1}X^ty \tag{3.28} \end{equation}\]
And solving for \(\sigma_{ML}^2\):
\[\begin{equation} \sigma_{ML}^2=\frac{(y-X\beta)^t(y-X\beta)}{n}=\frac{e^te}{n} \tag{3.29} \end{equation}\]
We see that ML estimator \(\sigma_{ML}^2\) of the error variance \(\sigma^2\) is smaller than that of OLS (3.20) (because of the ML’s denominator \(n\) is larger than the OLS’ denominator \(n-k\)) and hence not unbiased. In large samples, however, there will be no much difference between these two quantities, and The ML estimators are consistent and converge to that of the OLS method.
Instead of using the distribution of each single \(\varepsilon_t\), we can use the multivariate version for \(n\) dimensions so that the likelihood will be:
\(L=\frac{1}{(2\pi)^{\frac{n}{2}}det(\Omega_{\varepsilon})^{\frac{1}{2}}}exp\bigg(-\frac{1}{2}\varepsilon^t\underbrace{\Omega_{\varepsilon}^{-1}}_{=(\sigma^2I_n)^{-1}}\varepsilon\bigg)\).
\(L=\bigg(\frac{1}{\sqrt{2\pi\sigma^2}}\bigg)^nexp\bigg(-\frac{\varepsilon^t\varepsilon}{2\sigma^2}\bigg)\)
Besides its use when the errors are normally distributed, this method is used, often, when the errors follow any other distribution such as poisson or binomial, etc.. Using the above example again 3.1, we can apply the maximum likelihood method in R with the function glm. Since the the OLS and the maximum likelihood methods are equivalents under the classical assumptions , we expect, thus, to get the same estimates.
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 7.913295 | 3.2101540 | 2.465082 | 0.0206196 |
| x2 | 1.395749 | 0.1413768 | 9.872542 | 0.0000000 |
| x3 | -0.818892 | 0.1443155 | -5.674317 | 0.0000057 |
| x4 | 2.008222 | 0.1636282 | 12.273079 | 0.0000000 |
import statsmodels.api as sm
Y = df_p_sample['Y']
X = df_p_sample[['X2','X3','X4']]
X = sm.add_constant(X)
mod_mlm_p = sm.GLM(Y, X, family=sm.families.Gaussian()).fit()
res_p = mod_mlm_p.summary2().tables[1]| Coef. | Std.Err. | z | P>|z| | [0.025 ] | ||
|---|---|---|---|---|---|---|
| const | 3.5782321 | 2.8764907 | 1.243957 | 0.2135152 | -2.0595862 | 9.2160503 |
| X2 | 1.1923048 | 0.1709840 | 6.973194 | 0.0000000 | 0.8571822 | 1.5274273 |
| X3 | -0.4177008 | 0.1537923 | -2.716006 | 0.0066075 | -0.7191282 | -0.1162735 |
| X4 | 2.0368269 | 0.1653634 | 12.317274 | 0.0000000 | 1.7127205 | 2.3609332 |