#https://rpubs.com/neeraj1990/HM_model_R

rm(list = ls())

#reading Data
price <- read.csv(file = "https://raw.githubusercontent.com/Neeraj2308/HM-Model/master/Data.csv")
names = names(price)

#removing NA
price <- price[complete.cases(price), ]

#remove date
price <- price[, -1]

#number of securities
ns <- ncol(price)

#Getting return (we have taken log return)
ret <- apply(log(price), 2, diff)

#number of time periods
n=nrow(ret)

# Basic summaries
sapply(as.data.frame(ret), summary)

par(mfrow=c(2,4))
for (i in 1:4)
  ts.plot(price[,i])
for (i in 1:4)
  ts.plot(ret[,i])
  
  
  
# 1) Global minimum variance portfolio

#expected return and sd
er <- apply(ret, 2, mean) #daily expected return
std <- apply(ret, 2, sd) #daily sd

#covariance matrix
covm <- cov(ret)
covm.static = covm

# Scatter plots
pairs(ret)

#correlation matrix
corm <- cor(ret)

#Global minimum variance (gmv) portfolio
gmv.weights = function(covm,ns){
  Am <- rbind(2*covm, rep(1, ns))
  Am <- cbind(Am, c(rep(1, ns), 0))
  b  <- c( rep(0, ns), 1)
  w.gmv <- solve(Am) %*% b 
  w.gmv <- w.gmv[-(ns+1),]  #last value is lambda constraint. Hence not relevant
  return(w.gmv)
}
w.gmv = gmv.weights(covm,ns)

sum(w.gmv) #sum of weights are 1
w.gmv

port = ret%*%w.gmv
var.gmv <- t(w.gmv) %*% covm %*% w.gmv
ret.gmv <- t(w.gmv) %*% er

tab = rbind(cbind(er,std),
c(ret.gmv,var.gmv))
rownames(tab) = c(names,"P.gmv")
tab

par(mfrow=c(1,3))
ts.plot(ret,col=2:5)
lines(port,col=6)
legend("topright",legend=c(names,"P.gmv"),col=2:6,lty=1,lwd=2)
boxplot(cbind(ret,port),names=c(names,"P.gmv"),col=2:6)
plot(density(ret[,1]),col=2,xlim=range(ret),xlab="",main="",ylim=c(0,45))
lines(density(ret[,2]),col=3)
lines(density(ret[,3]),col=4)
lines(density(ret[,4]),col=5)
lines(density(port),col=6)


# Factor stochastic volatility

library(factorstochvol)

fsv = fsvsample(ret)

covmt    = array(0,c(n,ns,ns))
cormt    = array(0,c(n,ns,ns))
ws       = matrix(0,n,ns)
var.gmvt = rep(0,n)
for (t in 1:n){
  l = 0
  for (j in 1:ns)
  for (i in j:ns){
    l = l + 1
    covmt[t,i,j] = fsv$runningstore$cov[t,l,1]
    covmt[t,j,i] = covmt[t,i,j]
  }
  w.gmvt = gmv.weights(covmt[t,,],ns)
  var.gmvt[t] = t(w.gmvt) %*% covmt[t,,] %*% w.gmvt
  ws[t,] = w.gmvt 
  iD = diag(1/sqrt(diag(covmt[t,,])))
  cormt[t,,] = iD%*%covmt[t,,]%*%iD
}

par(mfrow=c(2,2))
for (i in 1:ns){
  ts.plot(sqrt(covmt[,i,i]),ylab="Standard deviation",main=names[i],ylim=c(0,0.06))
  abline(h=sqrt(covm.static[i,i]),col=2)
  lines(sqrt(var.gmvt),col=3)
}


par(mfrow=c(2,3))
ts.plot(cormt[,1,2],ylab="Correlation",main=paste(names[1],"-",names[2],sep=""),ylim=range(cormt))
abline(h=corm[1,2],col=2)
ts.plot(cormt[,1,3],ylab="Correlation",main=paste(names[1],"-",names[3],sep=""),ylim=range(cormt))
abline(h=corm[1,3],col=2)
ts.plot(cormt[,2,3],ylab="Correlation",main=paste(names[2],"-",names[3],sep=""),ylim=range(cormt))
abline(h=corm[2,3],col=2)
ts.plot(cormt[,1,4],ylab="Correlation",main=paste(names[1],"-",names[4],sep=""),ylim=range(cormt))
abline(h=corm[1,4],col=2)
ts.plot(cormt[,2,4],ylab="Correlation",main=paste(names[2],"-",names[4],sep=""),ylim=range(cormt))
abline(h=corm[2,4],col=2)
ts.plot(cormt[,3,4],ylab="Correlation",main=paste(names[3],"-",names[4],sep=""),ylim=range(cormt))
abline(h=corm[3,4],col=2)


par(mfrow=c(2,2))
for (i in 1:ns){
  ts.plot(ws[,i],ylim=range(ws),ylab="Weight",main=names[i])
  abline(h=w.gmv[i],col=2)
}

par(mfrow=c(1,1))
ts.plot(ws[,4],ylim=c(0,1),main="Metal",ylab="Weights")
lines(abs(ret[,4]),col=4,type="h")
abline(h=w.gmv[4],col=2)