lmq2.R

#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.



# This file calculates a regression using maximum likelihood.



# Read in the Ostrom data



Ostrom <- read.table("C:/users/Wood/Documents/My Teaching/Maximum Likelihood/Data/ostrom.dat", header=TRUE)

attach(Ostrom)

Ostrom

summary(Ostrom)



# Generate time series plots of the data.

par(mfrow=c(2,1))

plot(year,USSR)

plot(year,US)



# Plot US versus USSR

plot(USSR,US)



# First, calculate a regression the easy way, plot the

# residuals, and do an F test for whether the coefficient on USSR=0

OLS.model <- lm(US ~ USSR , data=Ostrom)

summary(OLS.model)

plot(residuals(OLS.model))



library(car)

# Set up the F test

r <- c(0,1)

q <- c(0)

linearHypothesis(OLS.model, r, q)



# We can also get various diagnostics from the preceding model as follows.

# These are only a few. We explore diagnostics for heteroskedasticity and autocorrelation in the next lesson.

anova(OLS.model)

predict(OLS.model)

logLik(OLS.model)

AIC(OLS.model)

influence.measures(OLS.model)

library(lmtest)

reset(US ~ USSR, power=3, type='fitted')



# Now let's fit the Normal Maximum Likelihood model.



# First, put the data into matrices for the MLE procedure

x <- cbind(1,as.matrix(USSR))

y <- as.matrix(US)

ones <- x[,1]



# Calculate K and n for later use

K <- ncol(x)

K

K1 <- K + 1

n <- nrow(x)

n



# Define the function to be optimized



llik.regress <- function(par,X,Y) {
  
  Y <- as.vector(y)
  
  X <- as.matrix(x)
  
  xbeta <- X%*%par[1:K]
  
  Sig <- par[K1:K1]
  
  sum(-(1/2)*log(2*pi)-(1/2)*log(Sig^2)-(1/(2*Sig^2))*(y-xbeta)^2)
  
}

llik.regress



# Now let's use the above function to estimate the model.



model <- optim(c(15,0.9,13),llik.regress, method = "BFGS", control = list(trace=1,maxit=100,fnscale = -1),
               
               hessian = TRUE)

model



# Notice that the coefficients from this model are the same as those produced through OLS estimation

#Now let's get the variances, standard errors, Z statistics, and log likelihood from our model.



# Calculate the variance matrix from the Hessian matrix. 



v <- -solve(model$hessian)

v



# Calculate the standard errors from the variance matrix.



se <- sqrt( diag(v))

se



# Calculate the z statistics from the coefficients and standard errors



b <- model$par

b

zstat <-b/se

zstat



# Calculate p-values for the z statistics

pzstat <- 2* (1 - pnorm(abs(zstat)))

pzstat



# Put the results together in a table.

table <- cbind(b,se,zstat,pzstat)

table





# Note that the estimate of the error variance is biased by n/n-K

# we can correct the preceding estimate of the standard error of #estimates as follows

Sig <- model$par[3]

Sig

Sigunb <- sqrt(n/(n-K)*Sig^2)

Sigunb



# The square of the z-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. 

LU <- model$value

LU



# Now let's do a likelihood ratio test that the coefficient on USSR=0 using "LU" from above.

# First, restimate the model with the restriction that USSR=0.  Save the log likelihood from this model



llik.regressR <- function(par,X,Y) {
  
  Y <- as.vector(y)
  
  X <- as.matrix(ones)
  
  xbeta <- X%*%par[1:K-1]
  
  Sig <- par[K:K]
  
  sum(-(1/2)*log(2*pi)-(1/2)*log(Sig^2)-(1/(2*Sig^2))*(y-xbeta)^2)
  
}

llik.regressR



# Now let's use the above function to estimate the model.



modelR <- optim(c(140,74),llik.regressR, method = "CG", control = list(trace=1,maxit=100,fnscale = -1),
                
                hessian = TRUE)

modelR

LR <- modelR$value

LR



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



llratio <- -2*(LR-LU)

llratio

llratioPvalue <- 1-pchisq(llratio,K-1)

llratioPvalue





# In some instances it will be useful to perform matrix calculations on MLE outputs. For example,

# in computing quadratic forms. As an example of using matrices, let's illustrate the preceding

# regression using some of the matrix capabilities of R



# We defined the appropriate matrices for this above.



# Compute coefficient estimates using matrices.



bols <- solve(t(x) %*% x) %*% t(x) %*% y

bols



# Compute and list predicted values and residuals



yhat <- x %*% bols

yhat

errors <- y - yhat   

errors



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



SSE <- t(errors) %*% errors

SSE

S2 <- (1/n) * SSE

S2



# Compute the maximum likelihood covariance matrix of coefficients. The diagonal of this

# matrix contains the standard errors.



covb <-  as.numeric(S2) * (solve(t(x) %*% x))

covb



# Extract the standard errors



se <- sqrt(diag(covb))

se



# Calculate t statistics



tstat <- bols/se

tstat