ESTIMATING A LINEAR REGRESSION USING MLE

/*The purpose of this session is to introduce you to the MLE of the normal general linear model. This approach to linear regression” forms the statistical basis for hypothesis testing found in most econometrics textbooks. Thus, the normal general linear model is important from the standpoint of justifying least squares. More important, this model serves as a tool for understanding maximum likelihood estimation of many time series models, models with heteroskedastic disturbances, and models with non-normal disturbances.*/

clear

use "c:\users\B. Dan Wood\My Documents\My Teaching\Maximum Likelihood\Data\ostrom.dta", clear

summarize us ussr

/* Now let's fit a regression the easy way, plot and save the residuals and do an F test for whether the coefficient on USSR=0.*/

regress us ussr

predict res, res

test ussr=0

graph twoway (scatter res year), yline(0)

/* Now let's fit the Normal Maximum Likelihood model. Recall from last week the procedure that was used to define the function is:*/

program define regressest

args lnf theta1 theta2

quietly replace `lnf'=-ln(sqrt(2*_pi*`theta2'^2))-1/(2*`theta2'^2)*($ML_y1-`theta1')^2

end

*Now lets use the above function to estimate the model

ml model lf regressest (eq1:us=ussr) (eq2:us=)

ml maximize

ml graph

*Now let's estimate the standard errors with the outer product of gradient method.

* Note that this would be given normally if the technique is bhhh.

ml model lf regressest (eq1:us=ussr) (eq2:us=), technique(bhhh) vce(opg)

ml maximize

/*Notice that the coefficients from this model are the same as those produced through OLS estimation. Also recall that the estimate of the disturbance variance is biased. The following code will correct for this.*/

scalar sigma=[eq2]_cons

display sigma

scalar sigunb=sqrt(sigma^2*(_N/(_N-2)))

display sigunb

test ussr=0

/* The square of the t-statistic on USSR is a Wald test that USSR=0. We might want to test this hypothesis with other methods, including either a likelihood ratio or LaGrange multiplier. Therefore, we save the log-likelihood for later testing below. We save it using two methods, the first in the standard way, and the second to implement STATA's canned likelihood ratio test procedure.*/

scalar lu=e(ll)

lrtest, saving(0)

/* Now let's do a likelihood ratio test that the coefficient on USSR=0 using “LU” from above and also using STATA's canned procedure for likelihood ratio tests. First, restimate the model with the restriction that USSR=0. */

ml model lf regressest (eq1:us=) (eq2:us=)

ml maximize

/*Now calculate the likelihood ratio test using the saved value of the log likelihood function.*/

scalar lr=e(ll)

display lr

display lu

scalar llratio=-2*(lr-lu)

display llratio

/* Now calculate the test using STATA's canned procedure */

lrtest, using(0)

*Now clear the program "regresset"

program drop regressest

/* In some instances it will be useful to perform matrix calculations on MLE outputs. For example, in computing quadratic forms. You may also want to use the “mkmat” command in simplifying your commands. As an example of using matrices, let's illustrate the preceding regression using some of the matrix capabilities of Stata. */

*Recall from previous statistics courses that the constant in an OLS model is just a column of 1s

gen one=1

* Now, define the dependent and independent variables in terms of named matrices.

mkmat us, matrix(y)

mkmat one ussr, matrix(X)

/* Compute coefficient estimates using matrices. Note that "inv" before a matrix is STATA's language for matrix inverse. */

matrix bols=inv(X'*X)*X'*y

matrix list bols

* Compute and list predicted values and residuals

matrix yhat=X*bols

matrix list yhat

matrix errors=y-yhat

matrix list errors

* Compute and display the maximum likelihood estimate of the error variance

matrix SSE=errors'*errors

matrix list SSE

matrix S2=(1/_N)*SSE

matrix list S2

/* Compute the maximum likelihood covariance matrix of coefficients. The diagonal of this matrix contains the standard errors. */

matrix covb=S2*inv(X'*X)

matrix list covb

/* Extract the mle standard errors by first getting the diagonal of covb and then taking the sqrt. */

matrix varb=vecdiag(covb)

matrix list varb

matrix seconstant = sqrt(varb[1,1])

matrix list seconstant

matrix seussr = sqrt(varb[1,2])

matrix list seussr

/*compute the mle t statistics */

matrix t_constant = bols[1,1]/sqrt(varb[1,1])

matrix list t_constant

matrix t_ussr = bols[2,1]/sqrt(varb[1,2])

matrix list t_ussr