**ESTIMATING THE HETEROSKEDASTIC/AUTOCORRELATED LINEAR
REGRESSION USING MLE**

The purpose of
this session is to show you how to estimate and test the heteroskedastic
and/or autocorrelated normal general linear models
using MLE.

HETEROSKEDASTICITY:
Heteroskedasticity can be treated directly in the
context of the normal MLE simply by specifying an equation to reflect the form
of the heteroskedasticity in place of the variance
term in the log likelihood function. This equation can take many different forms
to correspond with the type of heteroskedasticity.
There could simply be a weight that reflects the variance of each disturbance
(White's approach), or the heteroskedasticity could
be some linear combination of independent variables that may (or may not) be in
the equation for the conditional mean. Other forms include dependent variable heteroskedasticity, where the weighting term is some
function of the dependent variable (logs, absolute values, etc.), ARCH, and
GARCH. All of these are easily treated in the context of the normal MLE. The
example below shows one form, but this program could easily be modified to
include all of the other forms.

/* This file
demonstrates maximum likelihood estimation of normal models with heteroskedasticity

of various forms. */

clear

/* The next line
will read a data file from Greene, Chapter 11. Change the path to find data */

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

/* Set the sample
and generate a required variable to replicate Greene. */

gen id = _n

keep if (1<=id<=100)

drop if avgexp==0

summarize

gen income2=income^2

/* First, let's
estimate the regression using least squares and look at the residuals.

This makes it
clear that income is the source of the heteroskedasticity.
*/

regress avgexp age ownrent income income2

predict res, res

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

hettest

hettest income

hettest income income2

hettest , mtest(b)

whitetst

imtest, white

/* Other
available tests in STATA */

/* Influence Statistics
*/

predict dfits, dfits

list dfits

predict cooksd, cooksd

list cooksd

predict dfits, dfits

list dfits

predict wd, welsch

list wd

predict covr, covratio

list covr

predict dincome, dfbeta(income)

list dincome

dfbeta

/* Ramsey's reset
test for omitted variables */

ovtest

/* Now, let's
estimate the heteroskedastic regression the easy way.
The following

perform's White's procedure for robust standard
errors. */

regress avgexp age ownrent income income2, robust

/* You can also
specify a weighted least squares procedure. */

regress avgexp age ownrent income income2 [aweight=income]

/*You can test
linear hypotheses using a Wald procedure following STATA's canned

estimation commands. Here are some examples.*/

test income=0

test income=200

test income income2

test income-age=1

test income=age

test .5*income+2*age=1

/* Now let's introduce STATA's margins command to enhance the interpretation of
the model. Will use this later.

i. identifies a
variable as a factor variable. */

regress avgexp age i.ownrent income income2

margins ownrent

marginsplot

/* Using margins
to interpret an interaction effect. c. identifies a
variable as a continuous variable.

## creates the interaction. */

regress avgexp age i.ownrent i.ownrent##c.age income income2

margins ownrent,
at(age=(20(1)55))

marginsplot, xdimension(at(age))

/* Now, let's
develop our own procedure using the maximize command. Fit the Normal

Maximum
Likelihood heteroskedastic regression model. We use the same procedure as last time.*/

program define regressest

args lnf theta1 theta2

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

end

/*Now lets use
the above function to estimate the model.*/

ml model lf regressest
(eq1:avgexp=age ownrent income income2)
(eq2:avgexp=income), technique(bhhh)

ml maximize

/* Not close to
the wls estimates */

/* Let's try
putting in our own starting values */

ml model lf regressest
(eq1:avgexp=age ownrent income income2) (eq2:avgexp=income),
technique(bfgs)

ml init -3.57 -3.81 254.82 -16.41 -261 1 8 ,
copy

ml maximize

/* We can use a
Wald statistic to test for heteroskedasticity. The square of the t-statistic

on income in the second equation is a Wald test
that the coefficient on income=0.

Note that the test command does not allow
testing other hypotheses than zero when using

the ml
procedure. */

test [eq2]income=0

/*We might also
want to test the joint hypothesis that both income coefficients in the main

equation are zero. */

test [eq1]income income2

/*
Below save the log likelihood from the ml model for additional hypothesis
testing. */

est store A

/* Now let's do a
likelihood ratio test for no heteroskedasticity using
STATA's canned

procedure for likelihood ratio tests. First, restimate
the model with the restriction

of homoskedasticity. */

ml model lf regressest
(eq1:avgexp=age ownrent income income2) (eq2:avgexp=)

ml maximize

/* Now calculate
the likelihood ratio test for no heteroskedasticity using
STATA's canned

procedure */

lrtest A

/*Now clear the
program "regresset"*/

program drop regressest

exit

/* This file
demonstrates various methods for dealing with autocorrelation within STATA. */

clear

/* The next line
will read a data file. Change the path to find data */

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

/* Set the Sample
*/

gen id = _n

keep if (1<=id<=22)

summarize us ussr

/* Now let's fit a
regression the easy way, save and plot the residuals and do a Wald

test for whether the coefficient on USSR=0.*/

regress us ussr

predict res, res

test ussr=0

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

tsset year

/* Now let's
compute robust (autocorrelation and heteroskedasticity
consistent) standard

errors using the Newey-West
procedure. */

newey us ussr, lag(1)

/* Now let's do Prais-Winston regression.
This is an iterative procedure that jointly

estimates the coefficients and autocorrelation coefficient. */

prais us ussr

/*Now let's do a
maximum likelihood procedure for Cochrane-Orcutt
estimation. Note that

unlike the prais
procedure above, we are not fully efficient since we are dropping the

first observation.
It would take a bit more programming to retrieve the first observation

using the Prais-Winston
transformation. */

gen uslag=us[_n-1]

gen ussrlag=ussr[_n-1]

drop if(year==67)

program define arregressest

args lnf theta1 theta2
theta3 theta4

quietly replace `lnf'= /*

*/ -(1/2)*ln(2*_pi)-(1/2)*ln(`theta4'^2)-(1/(2*`theta4'^2))*(($ML_y1-`theta1'-`theta2'*ussr)
-`theta3'*(uslag-`theta1'-`theta2'*ussrlag))^2+(1/(2*_N))*ln(1-`theta3'^2)

end

/*Now lets use
the above function to estimate the model.
*/

ml model lf arregressest
(b0:us=) (b1:) (rho:) (sigma:)
, technique(nr)

ml init 7 .9 .9 7, copy

ml maximize

/* Now clear the
program "arregresset." */

program drop arregressest

exit