############################################################################################
#
# Bayesian linear regression with conjugate prior
#
############################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoBooth.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
############################################################################################


# Dataset: National Longitudinal Survey of Youth dataset
data = matrix(scan("logwages-yearseducation.txt"),byrow=TRUE,ncol=2)
n = nrow(data)
y = data[,1]
x = data[,2]
In = diag(1,n)

pdf(file="nlsy-graph1.pdf",width=10,height=10)
par(mfrow=c(1,1))
plot(table(x)/n*100,xlab="Years of schooling completed",ylab="Percentage",lwd=20)
dev.off()

pdf(file="nlsy-graph2.pdf",width=10,height=10)
par(mfrow=c(1,1))
hist(y,xlab="Log hourly wage",prob=TRUE,main="",nclass=20);box()
zzz = seq(min(y),max(y),length=1000)
m = mean(y)
sd = sqrt(var(y))
lines(zzz,dnorm(zzz,m,sd),col=2,lwd=3)
dev.off()

pdf(file="nlsy-graph3.pdf",width=10,height=10)
par(mfrow=c(1,1))
plot(x,y,xlab="Years of schooling completed",ylab="Log hourly wage",pch=16)
abline(lm(y~x),col=2,lwd=3)
dev.off()

##################
# COMPLETE MODEL
##################
d = 2
X = cbind(1,x)

# Prior hyperparameters
mu      = rep(0,d)
VbI    = diag(10,d)
VbII   = diag(10^100,d)
lambda = 0.1333
nu     = 6

# Marginal posterior distributions of sigma2 and beta
par1   = (nu+n)/2

# Under prior I
C1I    = solve(t(X)%*%X+solve(VbI))
mu1I   = C1I%*%(t(X)%*%y+solve(VbI)%*%mu)
par2I  = (nu*lambda+t(y)%*%y+t(mu)%*%solve(VbI)%*%mu-t(mu1I)%*%solve(C1I)%*%mu1I)/2
EI   = par2I/(par1+1)
SDI  = EI/sqrt(par1+2)
nu1lambda1I=nu*lambda+t(y-X%*%mu)%*%(In-X%*%solve(solve(VbI)+t(X)%*%X)%*%t(X))%*%(y-X%*%mu)
sdI=sqrt((nu1lambda1I/(nu+n))*diag(C1I))

# Under prior II
C1II   = solve(t(X)%*%X+solve(VbII))
mu1II  = C1I%*%(t(X)%*%y+solve(VbII)%*%mu)
par2II = (nu*lambda+t(y)%*%y+t(mu)%*%solve(VbII)%*%mu-t(mu1II)%*%solve(C1I)%*%mu1II)/2
EII  = par2II/(par1+1)
SDII = EII/sqrt(par1+2)
nu1lambda1II=nu*lambda+t(y-X%*%mu)%*%(In-X%*%solve(solve(VbII)+t(X)%*%X)%*%t(X))%*%(y-X%*%mu)
sdII=sqrt((nu1lambda1II/(nu+n))*diag(C1II))
round(rbind(cbind(mu1I,sdI,mu1II,sdII),c(EI,SDI,EII,SDII)),5)

# Computing log prior predictives (or log marginal likelihoods)
PI       = (In-X%*%solve(solve(VbI)+t(X)%*%X)%*%t(X))/lambda
PII      = (In-X%*%solve(solve(VbII)+t(X)%*%X)%*%t(X))/lambda
logdetI  = n*log(lambda)+log(det(solve(VbI)+t(X)%*%X))+log(det(VbI))
logdetII = n*log(lambda)+log(det(solve(VbII)+t(X)%*%X))+log(det(VbII))
lpredI   = -0.5*logdetI-0.5*(nu+n)*log(1+t(y-X%*%mu)%*%PI%*%(y-X%*%mu)/nu)
lpredII  = -0.5*logdetII-0.5*(nu+n)*log(1+t(y-X%*%mu)%*%PII%*%(y-X%*%mu)/nu)
c(lpredI,lpredII)

##################
# RESTRICTED MODEL
##################
d   = 1
X = matrix(1,n,1)

# Prior hyperparameters
mu     = rep(0,d)
VbI    = diag(10,d)
VbII   = diag(10^100,d)
lambda = 0.1333
nu     = 6

# Marginal posterior distributions of sigma2 and beta
par1   = (nu+n)/2

# Under prior I
C1I    = solve(t(X)%*%X+solve(VbI))
mu1I   = C1I%*%(t(X)%*%y+solve(VbI)%*%mu)
par2I  = (nu*lambda+t(y)%*%y+t(mu)%*%solve(VbI)%*%mu-t(mu1I)%*%solve(C1I)%*%mu1I)/2
EI   = par2I/(par1+1)
SDI  = EI/sqrt(par1+2)
nu1lambda1I=nu*lambda+t(y-X%*%mu)%*%(In-X%*%solve(solve(VbI)+t(X)%*%X)%*%t(X))%*%(y-X%*%mu)
sdI=sqrt((nu1lambda1I/(nu+n))*diag(C1I))

# Under prior II
C1II   = solve(t(X)%*%X+solve(VbII))
mu1II  = C1I%*%(t(X)%*%y+solve(VbII)%*%mu)
par2II = (nu*lambda+t(y)%*%y+t(mu)%*%solve(VbII)%*%mu-t(mu1II)%*%solve(C1I)%*%mu1II)/2
EII  = par2II/(par1+1)
SDII = EII/sqrt(par1+2)
nu1lambda1II=nu*lambda+t(y-X%*%mu)%*%(In-X%*%solve(solve(VbII)+t(X)%*%X)%*%t(X))%*%(y-X%*%mu)
sdII=sqrt((nu1lambda1II/(nu+n))*diag(C1II))
round(rbind(cbind(mu1I,sdI,mu1II,sdII),c(EI,SDI,EII,SDII)),5)

# Computing log prior predictives (or log marginal likelihoods)
PI       = (In-X%*%solve(solve(VbI)+t(X)%*%X)%*%t(X))/lambda
PII      = (In-X%*%solve(solve(VbII)+t(X)%*%X)%*%t(X))/lambda
logdetI  = n*log(lambda)+log(det(solve(VbI)+t(X)%*%X))+log(det(VbI))
logdetII = n*log(lambda)+log(det(solve(VbII)+t(X)%*%X))+log(det(VbII))
lpredIr  = -0.5*logdetI-0.5*(nu+n)*log(1+t(y-X%*%mu)%*%PI%*%(y-X%*%mu)/nu)
lpredIIr = -0.5*logdetII-0.5*(nu+n)*log(1+t(y-X%*%mu)%*%PII%*%(y-X%*%mu)/nu)
c(lpredIr,lpredIIr)

# Computing log Bayes factors
c(lpredI-lpredIr,lpredII-lpredIIr)
# Log Bayes factor (Prior I):   84.7294
# Log Bayes factor (Prior II): -29.8886