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