#############################################################################################
#
# Worked Problem I - Bernoulli trials
# Bayesian Econometrics
# Copyright by Hedibert F. Lopes
# January 27th 2014
#
# Model: x1,...,xn iid Bernoulli(theta)
#
# Prior: theta ~ Beta(a,a) for a greater than or equal to 1.
#
# Data: Number of Petrobras ups/downs over the first n=16 business days of 2014.
#
############################################################################################
sn = 4
n  = 16
theta = seq(0,1,length=1000)

# Plot of the four prior distributions along with the likelihood function
pdf(file="bernoullitrials-likelihood-priors.pdf",width=10,height=8)
plot(theta,dbeta(theta,sn+1,n-sn+1),xlab=expression(theta),
     ylab="Density",type="l",ylim=c(0,10),lwd=2)
lines(theta,dbeta(theta,1,1),col=2,lwd=2)
lines(theta,dbeta(theta,2,2),col=3,lwd=2)
lines(theta,dbeta(theta,10,10),col=4,lwd=2)
lines(theta,dbeta(theta,50,50),col=6,lwd=2)
legend(0.6,10,legend=c("Likelihood","Beta(1,1)=U(0,1)","Beta(2,2)","Beta(10,10)","Beta(50,50)"),
      col=c(1:4,6),lty=1,lwd=2,bty="n",cex=1.5)
dev.off()


# Individual plots of prior, likelihood and posterior
pdf(file="bernoullitrials-likelihood-priors-posteriors.pdf",width=12,height=8)
par(mfrow=c(2,2))

plot(theta,dbeta(theta,5,13),xlab=expression(theta),ylab="Density",type="l",ylim=c(0,10),lwd=2)
a  = 1
lines(theta,dbeta(theta,a,a),col=2,lwd=2,lty=3)
lines(theta,dbeta(theta,a+sn,a+n-sn),col=2,lwd=2)
title("Prior: Beta(1,1)")
legend(0.6,10,legend=c("Likelihood","Prior","Posterior"),col=c(1,2,2),lty=c(1,3,1),lwd=2,bty="n",cex=1.5)

plot(theta,dbeta(theta,5,13),xlab=expression(theta),ylab="Density",type="l",ylim=c(0,10),lwd=2)
a  = 2
lines(theta,dbeta(theta,a,a),col=3,lwd=2,lty=3)
lines(theta,dbeta(theta,a+sn,a+n-sn),col=3,lwd=2)
title("Prior: Beta(2,2)")
legend(0.6,10,legend=c("Likelihood","Prior","Posterior"),col=c(1,3,3),lty=c(1,3,1),lwd=2,bty="n",cex=1.5)

plot(theta,dbeta(theta,5,13),xlab=expression(theta),ylab="Density",type="l",ylim=c(0,10),lwd=2)
a  = 10
lines(theta,dbeta(theta,a,a),col=3,lwd=2,lty=3)
lines(theta,dbeta(theta,a+sn,a+n-sn),col=3,lwd=2)
title("Prior: Beta(10,10)")
legend(0.6,10,legend=c("Likelihood","Prior","Posterior"),col=c(1,3,3),lty=c(1,3,1),lwd=2,bty="n",cex=1.5)

plot(theta,dbeta(theta,5,13),xlab=expression(theta),ylab="Density",type="l",ylim=c(0,10),lwd=2)
a  = 50
lines(theta,dbeta(theta,a,a),col=6,lwd=2,lty=3)
lines(theta,dbeta(theta,a+sn,a+n-sn),col=6,lwd=2)
title("Prior: Beta(50,50)")
legend(0.6,10,legend=c("Likelihood","Prior","Posterior"),col=c(1,3,3),lty=c(1,3,1),lwd=2,bty="n",cex=1.5)

dev.off()



# In this second part of the worked exercise, we downloaded a total of 3379 ups/downs  
# for Petrobras over a period of 14 years.  

data = scan("petrobras-2000-2014.txt")
sn = cumsum(data)
n = 1:length(sn)

# Vague (noninformative) prior
a1 = 1
b1 = 1
qlower1  = qbeta(0.025,a1+sn,b1+n-sn)
qmedian1 = qbeta(0.5,a1+sn,b1+n-sn)
qupper1  = qbeta(0.975,a1+sn,b1+n-sn)

# More informative, symmetric prior 
a2 = 50
b2 = 50
qlower2  = qbeta(0.025,a2+sn,b2+n-sn)
qmedian2 = qbeta(0.5,a2+sn,b2+n-sn)
qupper2  = qbeta(0.975,a2+sn,b2+n-sn)

# More informative, asymmetric prior 
a3 = 40
b3 = 10
qlower3  = qbeta(0.025,a3+sn,b3+n-sn)
qmedian3 = qbeta(0.5,a3+sn,b3+n-sn)
qupper3  = qbeta(0.975,a3+sn,b3+n-sn)

pdf(file="bernoullitrials-posteriorsummary.pdf",width=12,height=8)
ts.plot(cbind(qlower1,qmedian1,qupper1),ylim=c(0,1),xlab="Business days")
lines(qlower2,col=2)
lines(qmedian2,col=2)
lines(qupper2,col=2)
lines(qlower3,col=4)
lines(qmedian3,col=4)
lines(qupper3,col=4)
abline(h=0.5,lty=3)
legend(1000,0.9,legend=c(paste("Noninformative - Beta(",a1,",",b1,")",sep=""),
                         paste("Informative, symmetric - Beta(",a2,",",b2,")",sep=""),
                         paste("Informative, asymmetric - Beta(",a3,",",b3,")",sep="")),
       col=c(1,2,4),lwd=2,bty="n",cex=1.5)
dev.off()