names=c("XOM","IBM","CVX","GE","PG","T","JNJ","PFE","WFC","JPM","KO","MRK","VZ","SLB",
"WMT","MCD","PEP","COP","C","ABT","OXY","BAC","HD","DIS","CAT","UTX","MO","UNH","CVS",
"MMM","EMC","BMY","USB","AXP","UNP","HPQ","BA","DD","F","APC","CL","MDT","NKE","DOW",
"LLY","SO","EMR","HAL","DE","LOW","TGT","BAX","EOG","PNC","DHR","NEM","WAG","DVN","D",
"KMB","DUK","AGN","BK","FDX","EXC","TRV","ITW","GIS","BHI","GD","GLW","BEN","AIG",
"ADBE","ALTR","AMGN","AAPL","AMAT","ADSK","ADP","BIIB","BMC","CA","CELG","CERN",
"CSCO","CMCSA","COST","DELL","XRAY","DTV","EXPD","FAST","FISV","INTC","KLAC","LRCX",
"LLTC","MAT","MXIM","MU","MSFT","MYL","ORCL","PCAR","PAYX","QCOM","ROST","SIAL","SPLS",
"SYMC","TEVA","VRTX","VOD","XLNX")

pdf(file="GMVportfolios.pdf",width=15,height=15)

# Reading the dataset
data  = read.table("nyse73-nasdaq42.txt",header=TRUE)
p = ncol(data)-2
p1 = p*(p-1)/2


# Time series plots
par(mfrow=c(3,3))
for (i in 1:p){
  plot(data[,2+i],xlab="Months",ylab="",main=names[i],type="l",axes=FALSE)
  axis(2);box();axis(1,at=seq(12,240,by=36),lab=seq(1992,2011,by=3))
}
  
# Computing the pairwise correlations
corr = cor(data[,3:(p+2)])
correlation = NULL
for (i in 2:p)
  for (j in 1:(i-1))
    correlation = c(correlation,corr[i,j])
maxcorr = max(correlation)
for (i in 2:p)
  for (j in 1:(i-1))
    if (corr[i,j]==maxcorr){
    	  imax = i
    	  jmax = j
    }

par(mfrow=c(1,2))
hist(correlation,xlim=c(-0.25,1),breaks=seq(-0.25,1,by=0.0495),col=grey(0.9),axes=FALSE,main="")
box();axis(2);axis(1,at=c(-0.25,0,0.25,0.5,0.75,1.0))
abline(v=0,lty=3,col=2)
abline(v=0.25,lty=3,col=2)
abline(v=0.5,lty=3,col=2)
abline(v=0.75,lty=3,col=2)
abline(v=1,lty=3,col=2)
title("115 time series => 6555 correlations")
plot(data[,2+imax],data[,2+jmax],xlab=names[imax],ylab=names[jmax],col=grey(0.5),pch=16)
title(paste("Largest correlation = ",round(corr[imax,jmax],3),sep=""))


# Portfolio allocation based on k=2 assets
par(mfrow=c(1,1))
x  = data[,2+imax]
y  = data[,2+jmax]
w  = seq(0,1,by=0.02)
N  = length(w)
mP = rep(0,N+1)
sP = rep(0,N+1)
for (i in 1:N){
  P = w[i]*x+(1-w[i])*y
  mP[i] = mean(P)
  sP[i] = sqrt(var(P))
}
wx      = 1-(var(x)-cov(x,y))/(var(x)+var(y)-2*cov(x,y))
P       = wx*x+(1-wx)*y
mP[N+1] = mean(P)
sP[N+1] = sqrt(var(P))
plot(sP,mP,pch=16,xlab="Standard deviations",ylab="Means")
title(paste("w*",names[78],"+(1-w)*",ylab=names[96],sep=""))
points(sP[N+1],mP[N+1],col=2,pch=16,cex=1.5)
text(sP[N+1]+0.005,mP[N+1],"<= This is the global minimum variance portfolio: wx=0.1235",col=2)
text(0.138,0.0121,"These are portfolios for different values of w")

# Global mininum portfolio
m = matrix(apply(data[,3:(2+p)],2,mean),p,1)
S = var(data[,3:(2+p)])
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)

plot(s,m,pch=16,xlab="Standard deviations",ylab="Means",xlim=range(s,Smv),ylim=range(m,mmv))
points(Smv,mmv,pch=16,col=2,cex=1.5)
title("Global minimum variance portfolio (red dot) based on ALL 115 assets \n Standard deviations and means (black dots) for the 115 assets")
Smv.115 = Smv
mmv.115 = mmv


# Comparing the GMV portfolio based on ALL 115 assets
# to the GMV portfolio based on k assets for k<115 
set.seed(131654)
Y = data[,3:(2+p)]
M = 200
ks = c(10,25,50,75)
par(mfrow=c(2,2))
for (k in ks){
  mmv = rep(0,M)
  Smv = rep(0,M)
  for (i in 1:M){
    ind = sort(sample(1:p,size=k,replace=FALSE,prob=rep(1/p,p)))
    m = matrix(apply(Y[,ind],2,mean),k,1)
    S = var(Y[,ind])
    s = sqrt(diag(S))
    one = matrix(1,k,)
    iS  = solve(S)
    Smv[i] = 1/(t(one)%*%iS%*%one)[1,1]
    mmv[i] = Smv[i]*(t(m)%*%iS%*%one)[1,1]
    Smv[i] = sqrt(Smv[i])
  }
  plot(Smv,mmv,xlim=c(0.018,0.06),ylim=c(0.004,0.011),pch=16,
       xlab="Standard deviations",ylab="Means")
  points(Smv.115,mmv.115,pch=16,col=2,cex=1.5)
  dist = round(abs(min(Smv)-Smv.115)/Smv.115,6)
  title(paste("Global MV portfolio based on ",k," assets \n Relative ratio = ",dist,sep=""))
}
dev.off()