data = read.table("DJIA-19components-jan1981-dec2012.txt",header=TRUE)
attach(data)
p = ncol(data)-3
n = nrow(data)
company = names(data)[4:(3+p)]


ret = matrix(0,n-1,p)
for (i in 1:p)
  ret[,i] = 100*(data[2:n,3+i]/data[1:(n-1),3+i]-1)
  
data = data[2:n,]
  
ind = trunc(seq(2,n,length=5))
label=c("2/1/81","22/12/88","16/12/96","23/12/04","31/12/12")


pdf(file="DJIA-19components-jan1981-dec2012.pdf",width=15,height=10)
  
par(mfrow=c(2,2))
for (i in 1:p){
  plot(ret[,i],xlab="Days",ylab="Returns",main=company[i],ylim=range(ret),axes=FALSE,type="l")
  axis(2);box();axis(1,at=ind-1,lab=label)
}

years  = c(81:99,0:12)
ny     = length(years)
p1     = p*(p-1)/2
sd     = matrix(0,ny,p)
corr   = matrix(0,ny,p1)
for (i in 1:ny){
  	  cond = data[,3]==years[i]
    for (k in 1:p){
      sd[i,k] = sqrt(var(ret[cond,k])) 
    }
    l1 = 0
    for (k in 2:p){
      for (k1 in 1:(k-1)){
        l1 = l1 + 1
        corr[i,l1] = cor(ret[cond,k],ret[cond,k1])
      }     
    }
  }  
}


par(mfrow=c(4,5))
for (i in 1:p){
  plot(sd[,i],xlab="Days",ylab="St.Dev.",ylim=c(min(sd),max(sd)),main=company[i],axes=FALSE,type="b")
  axis(2);box();axis(1,at=1:ny,lab=1981:2012)
}

par(mfrow=c(4,4))
l1 = 0
for (k in 2:p)
  for (k1 in 1:(k-1)){
  	  l1 = l1 + 1 
    plot(corr[,l1],xlab="Days",ylab="Corr",ylim=c(-1,1),main="",axes=FALSE,type="b")
    axis(2);box();axis(1,at=1:ny,lab=1981:2012)
    title(paste(company[k]," - ",company[k1],sep=""))
    abline(h=0,lty=3)
    abline(h=-1,lty=3)
    abline(h=1,lty=3)
}

# Global mininum portfolio
ws = matrix(0,ny,p)
ss = rep(0,ny)
ms = rep(0,ny)
for (i in 1:ny){
  	  cond = (data[,3]==years[i])
    m = matrix(apply(ret[cond,],2,mean),p,1) 
    S = var(ret[cond,])
    s = sqrt(diag(S))
    one = matrix(1,p,)
    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)
    ws[i,] = wmv
    ss[i] = smv
    ms[i] = mmv
}

par(mfrow=c(1,1))
plot(ss,ms,xlab="Standard deviation",ylab="Mean",col=0)
text(ss,ms,1981:2012)
title("Annual global minimum variance portfolio")
text(ss[7]-0.1,ms[7],"Black Monday-Oct 19th=>",col=2)
text(ss[18]+0.1,ms[18],"<=Russian financial crisis",col=2)
text(ss[20]+0.1,ms[20],"<=Dotcom bubble",col=2)
text(ss[22]-0.1,ms[22],"Internet bublle bursting=>",col=2)
text(ss[28]+0.1,ms[28],"<=Financial crisis",col=2)

par(mfrow=c(4,5))
for (i in 1:p){
  plot(1981:2012,ws[,i],ylim=range(ws),type="h",xlab="YEARS",ylab="Weight",lwd=3,main=company[i])
  abline(h=0,lty=3)
  abline(h=0.25,lty=3)
  abline(h=0.5,lty=3)
  abline(h=-0.2,lty=3)
}
dev.off()