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


# SIMULATING SV DATA
set.seed(1234)
sim = svsim(500,mu=-10,phi=0.99,sigma=0.2)

pdf(file="graph1.pdf",width=12,height=10)
par(mfrow=c(2,1))
ts.plot(sim$vol,ylab="",main="Standard deviations")
ts.plot(sim$y,ylab="",main="Log-returns")
dev.off()


# PERFORMING STATISTICAL INFERENCE
draws = svsample(sim$y,draws=200000,burnin=1000,thinpara=100,
priormu=c(-10,1),priorphi=c(20,1.2),priorsigma=0.2)

names(draws)
[1] "para"     "latent"   "latent0"  "y"        "runtime" 
[6] "priors"   "thinning" "summary" 

pdf(file="graph2.pdf",width=12,height=4)
par(mfrow=c(1,3))
plot(density(draws$para[,1]),main=expression(mu),xlab="")
plot(density(draws$para[,2]),main=expression(phi),xlab="")
plot(density(draws$para[,3]),main=expression(sigma),xlab="")
dev.off()

stdevs = t(apply(exp(draws$latent/2),2,quantile,c(0.025,0.5,0.975)))

pdf(file="graph3.pdf",width=15,height=10)
par(mfrow=c(1,1))
ts.plot(stdevs,col=c(3,2,3))
lines(sim$vol)
dev.off()




# REAL DATA - PETROBRAS
data = read.csv("petrobras.csv",header=TRUE)

attach(data)

# Log returns
y = diff(log(as.numeric(petrobras)))
n = nrow(data)
ind = seq(1,n-1,length=6)

pdf(file="graph4.pdf",width=12,height=10)
plot(y,xlab="days",ylab="Log-returns",main="",type="l",axes=FALSE)
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
dev.off()

draws = svsample(y)

pdf(file="graph5.pdf",width=12,height=4)
par(mfrow=c(1,3))
plot(density(draws$para[,1]),main=expression(mu),xlab="")
plot(density(draws$para[,2]),main=expression(phi),xlab="")
plot(density(draws$para[,3]),main=expression(sigma),xlab="")
dev.off()


stdevs = t(apply(exp(draws$latent/2),2,quantile,c(0.025,0.5,0.975)))
stdevm= apply(exp(draws$latent/2),2,mean)


pdf(file="graph6.pdf",width=15,height=10)
par(mfrow=c(1,1))
plot(stdevs[,2],ylim=c(0,max(stdevs)),
       xlab="days",ylab="",main="",type="l",axes=FALSE)
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
lines(stdevs[,1],col=2)
lines(stdevs[,3],col=2)
lines(0.15*abs(y),type="h",col=grey(0.75))
dev.off()


fit.garch = garchFit(~garch(1,1),data=y,trace=F,include.mean=FALSE)    


fit.garch@fit$matcoef

pdf(file="graph7.pdf",width=15,height=10)
par(mfrow=c(1,1))
plot(fit.garch@sigma.t,ylab="Standard deviation",xlab="Days",type="l",axes=FALSE,
ylim=c(0,max(
fit.garch@sigma.t,stdevm,stdevs)))
axis(2);axis(1,at=ind,lab=data[ind,1]);box()
lines(stdevm,col=2,lwd=2)
lines(0.15*abs(y),type="h",col=grey(0.75))
lines(stdevs[,1],col=2,lty=2)
lines(stdevs[,3],col=2,lty=2)
legend("topleft",legend=c("SV-AR(1)","GARCH(1,1)","abs(y)"),col=c(2,1,grey(0.75)),bty="n",lty=1,lwd=2)
dev.off()