#########################################################################
#
# Data on 177 chief executive officers
# Goal: Examine the effects of firm performance on CEO salary.
# Source: Wooldridge (2012) page 111
#
#  1. salary                   1990 compensation, $1000s
#  2. age                      in years
#  3. college                  =1 if attended college
#  4. grad                     =1 if attended graduate school
#  5. comten                   years with company
#  6. ceoten                   years as ceo with company
#  7. sales                    1990 firm sales, millions
#  8. profits                  1990 profits, millions
#  9. mktval                   market value, end 1990, mills.
# 10. lsalary                  log(salary)
# 11. lsales                   log(sales)
# 12. lmktval                  log(mktval)
# 13. comtensq                 comten^2
# 14. ceotensq                 ceoten^2
# 15. profmarg                 profits as % of sales
#
#########################################################################
rm(list=ls())

pdf(file="ceosal2.pdf",width=12,height=10)
names = c("years with company",
"years as ceo with company",
"firm sales, millions (LOG)",
"market value, millions (LOG)")

data = matrix(scan("ceosal2.txt"),ncol=15,byrow=TRUE)

n = nrow(data)

lsalary  = data[,10]
comten   = data[,5]
ceoten   = data[,6]
lsales   = log(data[,7])
lmrkval  = log(data[,9])

SST = sum((lsalary-mean(lsalary))^2)

par(mfrow=c(2,2))
reg = lm(lsalary~comten)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-2))
U   = max(abs(reg$res))
L   = -U
plot(comten,lsalary,ylab="Compensation, $1000s (LOG)",xlab=names[1])
abline(reg$coef,col=2,lwd=2)
title(paste("R2=",R2,"%",sep=""))
plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))

reg = lm(lsalary~ceoten)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-2))
U   = max(abs(reg$res))
L   = -U
plot(ceoten,lsalary,ylab="Compensation, $1000s (LOG)",xlab=names[2])
abline(reg$coef,col=2,lwd=2)
title(paste("R2=",R2,"%",sep=""))
plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))

reg = lm(lsalary~lsales)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-2))
U   = max(abs(reg$res))
L   = -U
plot(lsales,lsalary,ylab="Compensation, $1000s (LOG)",xlab=names[3])
abline(reg$coef,col=2,lwd=2)
title(paste("R2=",R2,"%",sep=""))
plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))

reg = lm(lsalary~lmrkval)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-2))
U   = max(abs(reg$res))
L   = -U
plot(lmrkval,lsalary,ylab="Compensation, $1000s (LOG)",xlab=names[4])
abline(reg$coef,col=2,lwd=2)
title(paste("R2=",R2,"%",sep=""))
plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))


# Modelo com 2 covariaveis
# ------------------------
X = cbind(lsales,lmrkval)

pairs(X,labels=names[3:4])


reg = lm(lsalary~lsales+lmrkval)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-3))
U   = max(abs(reg$res))
L   = -U
summary(reg)

par(mfrow=c(1,2))
plot(reg$fit,lsalary,ylab="Compensation, $1000s (LOG)",xlab="Fitted")
title(paste("R2=",R2,"%",sep=""))
abline(0,1)

plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))


# Modelo com 4 covariaveis
# ------------------------
X = cbind(comten,ceoten,lsales,lmrkval)

pairs(X,labels=names[1:4])

reg = lm(lsalary~comten+ceoten+lsales+lmrkval)
R2  = round(100*(1-sum(reg$res^2)/SST),2)
s   = sqrt(sum(reg$res^2)/(n-5))
U   = max(abs(reg$res))
L   = -U
summary(reg)

par(mfrow=c(1,2))
plot(reg$fit,lsalary,ylab="Compensation, $1000s (LOG)",xlab="Fitted")
title(paste("R2=",R2,"%",sep=""))
abline(0,1)

plot(reg$fit,reg$res,ylim=c(L,U),xlab="Fitted",ylab="Residuals")
abline(h=0,col=2,lty=2)
abline(h=-2*s,col=2,lty=2)
abline(h=2*s,col=2,lty=2)
title(paste("Coverage = ",round(100*mean(abs(reg$res)<(2*s)),1),"%",sep=""))


dev.off()