#########################################################################################################
#
# Prices (and returns) of a few brazilian companies
#
# Period: 12/31/2002 to 12/21/2012 (more or less 2500 observations)
#
# Companies:
#
# PETROBRAS - Petroleo Brasileiro S.A. (PBR)
# VALE      - Vale S.A. (VALE)
# AMBEV     - Companhia de Bebidas das Americas (ADR)
# VIVO      - Telefonica Brasil S.A. (VIV)
# ITAU      - Itau Unibanco Banco Holding S.A. (ITUB)
# BRADESCO  - Banco Bradesco S.A. (BBD)
#
#########################################################################################################
#
# Hedibert Freitas Lopes
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# Email: hlopes@chicagobooth.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes
#
#########################################################################################################
data  = read.table("brazil.txt",header=TRUE)
names = c("PETROBRAS","VALE","AMBEV","VIVO","ITAU","BRADESCO")
n     = nrow(data)
ind   = c(1,253,505,757,1008,1259,1512,1764,2016,2268)


pdf(file="Brazil.pdf",width=20,height=15)
par(mfrow=c(1,1))
plot(data[,4],xlab="January 2nd 2003 - December 31st 2012",ylab="Prices",main="",axes=FALSE,type="l",lwd=1.5)
box();axis(2);axis(1,at=ind,lab=2000+data[ind,3]);axis(1,at=n,2013)
for (i in 2:6)
  lines(data[,3+i],col=i,lwd=1.5)
legend(4,60,legend=names,col=1:6,lty=1,lwd=2,bty="n")
  
par(mfrow=c(3,2))
for (i in 1:6){
  plot(data[,9+i],xlab="January 2nd 2003 - December 31st 2012",ylab="Returns",main="",
       axes=FALSE,type="l",ylim=range(data[,10:15]))
  box();axis(2);axis(1,at=ind,lab=2000+data[ind,3]);axis(1,at=n,2013)
  title(names[i])
}


L = min(data[,10:15])
U = max(data[,10:15])
breaks = seq(L,U,length=50)
par(mfrow=c(3,2))
for (i in 1:6)
hist(data[,9+i],breaks=breaks,prob=TRUE,xlab="",main=names[i])



years    = 3:12
months   = 1:12
corr     = array(0,c(120,6,6))
l        = 0 
for (year in years){
  for (month in months){
  	  l = l + 1
    corr[l,1,1] = sqrt(var(data[cond,10]))
    for (i in 2:6){
      for (j in 1:(i-1)){
  	      cond  = (data[,3]==year)&(data[,1]==month)
        corr[l,i,j] = cor(data[cond,9+i],data[cond,9+j])
      }  
      corr[l,i,i] = sqrt(var(data[cond,9+i]))
    }
  }
}

par(mfrow=c(6,6))
for (i in 1:6)
  for (j in 1:6){
    if (i==j){
      plot(corr[,i,j],main=names[i],ylab="St.Dev.",xlab="",type="l",axes=FALSE)
      axis(2);box();axis(1,at=c(seq(1,120,by=24),120),lab=c(seq(2003,2012,by=2),2013))
    }
    if (i>j){
      plot(corr[,i,j],main=paste(names[i],"x",names[j],sep=""),
           ylab="Correlation",xlab="",type="l",ylim=c(-0.304,1),axes=FALSE)
      axis(2);box();axis(1,at=c(seq(1,120,by=24),120),lab=c(seq(2003,2012,by=2),2013))
      abline(h=0,lty=3)
      abline(h=1,lty=3)
      abline(h=0.5,lty=3)
    }
    if (i<j){
      plot(0,axes=FALSE,col=0,xlab="",ylab="")
    } 
}

# Global mininum portfolio
p = 6
mmvs = rep(0,10)
Smvs = rep(0,10)
par(mfrow=c(2,5))
for (year in 3:12){
  Y = data[data[,3]==year,10:15]
  m = matrix(apply(Y,2,mean),p,1)
  S = var(Y)
  s = sqrt(diag(S))
  one = matrix(1,p,1)
  iS  = solve(S)
  Smv = 1/(t(one)%*%iS%*%one)[1,1]
  mmv = Smv*(t(m)%*%iS%*%one)[1,1]
  wmv = Smv*iS%*%one
  Smv = sqrt(Smv)
  plot(s,m,pch=16,xlab="Standard deviations",ylab="Means",xlim=c(0,7),ylim=c(-0.3,0.5),col=0)
  text(s,m,names,cex=0.5)
  text(Smv,mmv,"GMV",col=2)
  title(paste("YEAR = ",2000+year,sep=""))
  mmvs[year-2]=mmv
  Smvs[year-2]=Smv
}

par(mfrow=c(1,1))
plot(Smvs,mmvs,col=0,xlab="Standard deviations",ylab="Means")
text(Smvs,mmvs,2000+years,cex=1.5)
title("Global minimum variance portfolio by year")
text(1.65,0.292,"<= Good times",col=2,cex=1.5)
text(3.5,0.015,"Bad times =>",col=2,cex=1.5)

dev.off()