##########################################################################################
#
# Univariate and factor stochastic volatility estimation
#
# Copyright by Hedibert Freitas Lopes, PhD
# Date: March 21th 2019.
#
##########################################################################################
install.packages("stochvol")
install.packages("factorstochvol")
library("factorstochvol")
library("stochvol")

data = read.table("pbr-vale-google.txt",header=TRUE)
price = data[,4:6]
date = data[,1:3]
dates = c("20/3/06","17/6/09","13/9/12","15/12/15","19/3/19")

pdf(file="sv-fsv.pdf",width=12,height=7)
ind = trunc(seq(1,n,length=5))
par(mfrow=c(1,1))
plot(price[,1],xlab="Days",ylab="Price",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates)
lines(price[,2],col=2)
lines(price[,3]/20,col=3)
legend("top",legend=c("Petrobras","Vale","Google (*0.05)"),col=1:3,lty=1)

ret = apply(log(price),2,diff)
n = nrow(ret)
ind = trunc(seq(1,n,length=5))
covs = array(0,c(65,3,3))
cors = matrix(0,65,3)
for (i in 1:65){
  covs[i,,]=cov(ret[(50*(i-1)+1):(50*i),])
  cors[i,1] = covs[i,1,2]/sqrt(covs[i,1,1]*covs[i,2,2])
  cors[i,2] = covs[i,1,3]/sqrt(covs[i,1,1]*covs[i,3,3])
  cors[i,3] = covs[i,2,3]/sqrt(covs[i,2,2]*covs[i,3,3])
}

par(mfrow=c(1,1))
plot(ret[,1],xlab="Days",ylab="Log return",type="l",axes=FALSE,ylim=c(-0.3,1.0))
axis(2);box();axis(1,at=ind,lab=dates)
lines(0.4+ret[,2],col=2)
lines(0.7+ret[,3],col=3)
legend("topright",legend=c("Petrobras","Vale","Google"),col=1:3,lty=1)
abline(h=0)

par(mfrow=c(1,1))
plot(ret[,1],xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE,ylim=c(0,0.1),col=0)
axis(2);box();axis(1,at=ind,lab=dates)
lines(50*(1:65),sqrt(covs[,1,1]),col=1,pch=16,type="b")
lines(50*(1:65),sqrt(covs[,2,2]),col=2,pch=16,type="b")
lines(50*(1:65),sqrt(covs[,3,3]),col=3,pch=16,type="b")
legend("topright",legend=c("Petrobras","Vale","Google"),col=1:3,lty=1)
title("Sample standard deviations\n 50 observations")

par(mfrow=c(1,1))
plot(ret[,1],xlab="Days",ylab="Correlations",
     type="l",axes=FALSE,ylim=c(-0.2,1),col=0)
axis(2);box();axis(1,at=ind,lab=dates)
lines(50*(1:65),cors[,1],col=1,pch=16,type="b")
lines(50*(1:65),cors[,2],col=2,pch=16,type="b")
lines(50*(1:65),cors[,3],col=3,pch=16,type="b")
abline(h=cor(ret)[1,2],lty=2,col=1)
abline(h=cor(ret)[1,3],lty=2,col=2)
abline(h=cor(ret)[2,3],lty=2,col=3)
legend("bottomleft",legend=c("Petrobras,Vale","Petrobras,Google","Vale,Google"),col=1:3,lty=1)
title("Sample correlations\n 50 observations")



pbr.sv  = svsample(ret[,1])
vale.sv = svsample(ret[,2])
google.sv = svsample(ret[,3])
fsv3 = fsvsample(ret,runningstore=5)

par(mfrow=c(1,3))
boxplot(cbind(pbr.sv$para[,1],vale.sv$para[,1],google.sv$para[,1]),outline=FALSE,names=c("Petrobras","Vale","Google"))
title(expression(mu))
boxplot(cbind(pbr.sv$para[,2],vale.sv$para[,2],google.sv$para[,2]),outline=FALSE,names=c("Petrobras","Vale","Google"))
title(expression(phi))
boxplot(cbind(pbr.sv$para[,3],vale.sv$para[,3],google.sv$para[,3]),outline=FALSE,names=c("Petrobras","Vale","Google"))
title(expression(sigma))


pbr.vol = apply(exp(pbr.sv$latent/2),2,median)
vale.vol = apply(exp(vale.sv$latent/2),2,median)
google.vol = apply(exp(google.sv$latent/2),2,median)

par(mfrow=c(1,1))
plot(pbr.vol,xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE,ylim=c(0,0.12))
axis(2);box();axis(1,at=ind,lab=dates)
lines(vale.vol,col=2)
lines(google.vol,col=3)
points(50*(1:65),sqrt(covs[,1,1]),col=1,pch=16,cex=0.5)
points(50*(1:65),sqrt(covs[,2,2]),col=2,pch=16,cex=0.5)
points(50*(1:65),sqrt(covs[,3,3]),col=3,pch=16,cex=0.5)
legend("topright",legend=c("Petrobras","Vale","Google"),col=1:3,lty=1)
title("Univariate SV models")


stdu = cbind(pbr.vol,vale.vol,google.vol)
std = fsv3$runningstore$sd[,,1]

names = c("pbr","vale","goog","factor")
par(mfrow=c(1,4))
boxplot.matrix(t(fsv3$para[1,,]),main=expression(mu),outline=FALSE,names=names)
boxplot.matrix(t(fsv3$para[2,,]),main=expression(phi),outline=FALSE,names=names)
boxplot.matrix(t(fsv3$para[3,,]),main=expression(sigma),outline=FALSE,names=names)
boxplot.matrix(t(fsv3$facload[,1,]),outline=FALSE,main="Loadings",names=c(expression(beta[1]),expression(beta[2]),expression(beta[3])))

par(mfrow=c(2,2))
for (i in 1:3){
  plot(std[,i],xlab="Days",ylab="Std. Deviations",
       type="l",axes=FALSE,ylim=range(std[,i],stdu[,i]))
  axis(2);box();axis(1,at=ind,lab=dates)
  title(paste("Idiosyncratic factor ",i,sep=""))
}
i=4
plot(std[,i],xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates)
title("Common factor")


par(mfrow=c(1,1))
plot(sqrt(fsv3$runningstore$cov[,1,1]),xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates)
lines(pbr.vol,col=2)
legend("topright",legend=c("Univariate SV","Factor SV"),col=1:2,lty=1)
title("Petrobras")

plot(sqrt(fsv3$runningstore$cov[,4,1]),xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates)
lines(vale.vol,col=2)
legend("topright",legend=c("Univariate SV","Factor SV"),col=1:2,lty=1)
title("Vale")

plot(sqrt(fsv3$runningstore$cov[,6,1]),xlab="Days",ylab="Std. Deviations",
     type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates)
lines(google.vol,col=2)
legend("topright",legend=c("Univariate SV","Factor SV"),col=1:2,lty=1)
title("Google")
  
corrs = fsv3$runningstore$cor[,,1]
par(mfrow=c(1,1))
plot(corrs[,1],xlab="Days",ylab="Correlation",
     type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates)
title("Petrobras - Vale")
points(50*(1:65),cors[,1],col=2,pch=16)
legend("bottomleft",legend=c("FSV","Sample std. dev."),col=1:2,lty=1)

plot(corrs[,2],xlab="Days",ylab="Correlation",
     type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates)
title("Petrobras - Google")
points(50*(1:65),cors[,2],col=2,pch=16)

plot(corrs[,3],xlab="Days",ylab="Correlation",
     type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates)
title("Vale - Google")
points(50*(1:65),cors[,3],col=2,pch=16)


par(mfrow=c(3,3))
plot(sqrt(fsv3$runningstore$cov[,1,1]),xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[1]))
plot(corrs[,1],xlab="Days",ylab="",type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates);title(expression(rho[12]))
plot(corrs[,2],xlab="Days",ylab="",type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates);title(expression(rho[13]))
plot(fsv3$runningstore$cov[,2,1],xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[21]))
plot(sqrt(fsv3$runningstore$cov[,4,1]),xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[2]))
plot(corrs[,2],xlab="Days",ylab="",type="l",axes=FALSE,ylim=c(-0.2,1))
axis(2);box();axis(1,at=ind,lab=dates);title(expression(rho[23]))
plot(fsv3$runningstore$cov[,3,1],xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[31]))
plot(fsv3$runningstore$cov[,5,1],xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[32]))
plot(sqrt(fsv3$runningstore$cov[,6,1]),xlab="Days",ylab="",type="l",axes=FALSE)
axis(2);box();axis(1,at=ind,lab=dates);title(expression(sigma[3]))


dev.off()