library(MASS)
set.seed(23)
n<-200
Y<-rgamma(n,2,5)

hist(Y)
qqnorm(Y);qqline(Y)

qqnorm(log(Y));qqline(log(Y))
qqnorm(sqrt(Y));qqline(sqrt(Y))


Y.dot<-exp(mean(log(Y)))
Y.bc<-boxcox(Y~1,lambda=seq(-2,2,by=0.01))
plot(Y.bc,ylim=range(-600,-400),type='l')
Y.bc.lam<-Y.bc$x[which.max(Y.bc$y)]

Y.trans<-(Y^Y.bc.lam-1)/(Y.bc.lam*Y.dot^(Y.bc.lam-1))
hist(Y.trans)
qqnorm(Y.trans);qqline(Y.trans)


lval<-qchisq(0.95,1)/2

plot(Y.bc,xlim=range(-0.5,1),ylim=range(-500,-470),type='l')
abline(h=max(Y.bc$y)-lval,lty=2)

########################################################
set.seed(383)
n<-50
be0<-2
be1<-0.01
X<-runif(n,0,100)
muvec<-be0+be1*X

al<-0.75
be<-1
Y<-muvec+rgamma(n,al,be)-(al/be)

plot(X,Y,pch=18)

fit1<-lm(Y~X)
abline(fit1,col='blue')
abline(be0,be1,col='red')

#########################################
(CD<-read.csv('CobbDouglas.csv'))

y<-(as.numeric(gsub(",","",as.character(CD$Y))))
x1<-log(as.numeric(gsub(",","",as.character(CD$X2))))
x2<-log(as.numeric(gsub(",","",as.character(CD$X3))))


boxcox(lm(y~x1+x2),lambda=seq(-0.5,0.5,by=0.0001))
dev.print(pdf,file='CD-BoxCoxLam.pdf',width=8,height=8)

(ZD<-read.table('ZellnerData.txt',header=TRUE))

pdf('ZD-Data.pdf',paper='USr',width=11,height=6)
par(mfrow=c(1,2))
plot(ZD$Capital,ZD$ValueAdd,pch=18,xlab='Labor',ylab='ValueAdd')
plot(ZD$Labor,ZD$ValueAdd,pch=18,xlab='Labor',ylab='ValueAdd')
dev.off()

fit1<-lm(ValueAdd~Capital+Labor,data=ZD)
summary(fit1)

pdf('ZD-Resids.pdf',paper='USr',width=11,height=6)
par(mfrow=c(1,2))
plot(rstandard(fit1)~ZD$Capital,xlab='Capital',pch=18)
plot(rstandard(fit1)~ZD$Labor,xlab='Labor',pch=18)
dev.off()


lam.fit<-boxcox(lm(ValueAdd~Capital+Labor,data=ZD),lambda=seq(0,1.5,by=0.0001))
dev.print(pdf,file='ZD-BoxCoxLam.pdf',width=8,height=8)

(lambda.hat<-lam.fit$x[which.max(lam.fit$y)])


ytilde<-exp(mean(log(ZD$ValueAdd)))
Y<-ZD$ValueAdd
ZD$ynew<-(Y^lambda.hat-1)/(lambda.hat*ytilde^(lambda.hat-1))
fit2<-lm(ynew~Capital+Labor,data=ZD)
summary(fit2)

pdf('ZD-Trans-Resids.pdf',paper='USr',width=11,height=6)
par(mfrow=c(1,2))
plot(rstandard(fit2)~ZD$Capital,xlab='Capital',pch=18)
plot(rstandard(fit2)~ZD$Labor,xlab='Labor',pch=18)
dev.off()


pdf('ZD-QQ.pdf',paper='USr',width=11,height=6)
par(mfrow=c(1,2))
qqnorm(rstandard(fit1),main='Original',pch=19);qqline(rstandard(fit1))
qqnorm(rstandard(fit2),main='Transformed',pch=19);qqline(rstandard(fit2))
dev.off()