data = read.table("http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF6-1.txt",header=TRUE)

attach(data)

n = nrow(data)

lc = log(C)
lq = log(Q)
lf = log(PF)
ll = log(LF)

pdf(file="graph1.pdf",width=10,height=7)
par(mfrow=c(2,2))
plot(I,Q,col=I,pch=16,xlab="Airline company",ylab="Output",main="Output\nrevenue passenger miles")
plot(I,C/1000,col=I,pch=16,xlab="Airline company",ylab="Total cost",main="Total cost\n in $millions")
plot(I,PF/1000,col=I,pch=16,xlab="Airline company",ylab="Fuel price",main="Fuel price\nin $1000s")
plot(I,LF,col=I,pch=16,xlab="Airline company",ylab="Load factor",main="Load factor\naverage capacity utilization of the fleet")
dev.off()

pdf(file="graph2.pdf",width=10,height=7)
par(mfrow=c(2,2))
plot(1969+T,Q,col=I,pch=16,xlab="Year",ylab="Output",main="Output\nrevenue passenger miles")
for (i in 1:6) lines(1969+T[I==i],Q[I==i],col=i)
plot(1969+T,C/1000,col=I,pch=16,xlab="Year",ylab="Total cost",main="Total cost\n in $millions")
for (i in 1:6) lines(1969+T[I==i],C[I==i]/1000,col=i)
plot(1969+T,PF/1000,col=I,pch=16,xlab="Year",ylab="Fuel price",main="Fuel price\nin $1000s")
for (i in 1:6) lines(1969+T[I==i],PF[I==i]/1000,col=i)
plot(1969+T,LF,col=I,pch=16,xlab="Year",ylab="Load factor",main="Load factor\naverage capacity utilization of the fleet")
for (i in 1:6) lines(1969+T[I==i],LF[I==i],col=i)
dev.off()

pdf(file="graph3.pdf",width=10,height=7)
par(mfrow=c(2,2))
plot(lq,lc,col=I,xlab="log output",ylab="log total cost",pch=16)
plot(lf,lc,col=I,xlab="log fuel price",ylab="log total cost",pch=16)
plot(ll,lc,col=I,xlab="log load factor",ylab="log total cost",pch=16)
plot(lf,ll,col=I,xlab="log fuel price",ylab="log load factor",pch=16)
dev.off()

reg = lm(lc~lq+lf+ll)
plot(ll,reg$res,pch=16,xlab="Load factor",ylab="residuals",main="Model 2 (logs)")
abline(h=0,lty=2)
coef.ols = reg$coef
se.ols = sqrt(diag(summary(reg)$cov.unscaled))*summary(reg)$sigma

pdf(file="graph4.pdf",width=10,height=7)
plot(ll,reg$res,xlab="log load factor",ylab="Residuals")
abline(h=0,lty=2)
dev.off()


# Computing hetereskedasticity-robust standard errors
k = 3
p = k + 1
x = cbind(lq,lf,ll)
X = cbind(1,x)
se.hr = rep(0,ncol(X))
for (i in 1:p){
  reg1 = lm(X[,i]~X[,-i]-1)
  se.hr[i] = sqrt(sum((reg1$res^2)*(reg$res^2))/(sum(reg1$res^2))^2)
}
round(cbind(coef.ols,se.ols,se.hr),3)

# Hetereskedasticity-robust F test
esq     = reg$res^2
R2.e   = summary(lm(esq~X))$r.sq
Ftest   = R2.e/(1-R2.e)*(n-p)/k
pvalue = 1-pf(Ftest,k,n-p)
pvalue

# Generalized least squares
lesq = log(esq)
reg1 = lm(lesq~x)
g    = reg1$fit
hhat = exp(g)
reg.gls = lm(lc~x,weights=1/hhat)
se.gls = sqrt(diag(summary(reg.gls)$cov.unscaled))*summary(reg.gls)$sigma
coef.gls = reg.gls$coef


round(cbind(coef.ols,se.ols,se.hr,coef.gls,se.gls),3)

# Obtaining the maximum likelihood estimator
like = function(beta){
  -sum(dnorm(lc,beta[1]+beta[2]*lq+beta[3]*lp+beta[4]*ll,exp((beta[5]+beta[6]*ll)/2),log=TRUE))	
}
beta = c(lm(lc~lq+lp+ll)$coef,lm(lesq~ll)$coef)
nlm(like,beta)




# Changing the mean equation
comp2 = rep(0,n)
comp3 = rep(0,n)
comp4 = rep(0,n)
comp5 = rep(0,n)
comp6 = rep(0,n)
comp2[I==2]=1
comp3[I==3]=1
comp4[I==4]=1
comp5[I==5]=1
comp6[I==6]=1

reg = lm(lc~comp2+comp3+comp4+comp5+comp6+lq+lf+ll)
plot(ll,reg$res)

# Hetereskedasticity-robust F test
k = 8
p = k + 1
x = cbind(comp2,comp3,comp4,comp5,comp6,lq,lf,ll)
X = cbind(1,x)
esq     = reg$res^2
R2.e   = summary(lm(esq~x))$r.sq
Ftest   = R2.e/(1-R2.e)*(n-p)/k
pvalue = 1-pf(Ftest,k,n-p)
pvalue