rm(list=ls())
library(stochvol)

data = read.csv("vix-sp500.csv",header=TRUE)
n      = nrow(data)
attach(data)

ret = log(price[2:n]/price[1:(n-1)])

names  = c("SP500 Prices","SP500 Returns","VIX")



pdf(file="sp500-vix.pdf",width=12,height=8)
par(mfrow=c(2,1))
ind  = seq(1,n,length=5)
plot(ret,xlab="days",ylab="",main="SP500 returns",type="l",axes=FALSE)    
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
plot(vix,xlab="days",ylab="",main="VIX",type="l",axes=FALSE)    
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
dev.off()


draws = svsample(ret)
apply(draws$para,2,summary)

vol = apply(exp(draws$latent/2),2,mean)

pdf(file="sp500-vix-svol.pdf",width=12,height=8)
par(mfrow=c(2,1))
plot(vix,ylab="Implied volatility",xlab="Days",type="l",axes=FALSE,main="VIX")
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
plot(vol,ylab="Standard deviation",xlab="Days",type="l",axes=FALSE,main="")
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
title("Model-based volatility")
dev.off()