#install.packages("quantmod")
#install.packages("stochvol")
library(quantmod)
library(stochvol)

ticker <- "PBR"
start_date <- "2000-01-01"
end_date <- Sys.Date()
getSymbols(Symbols=ticker,src="yahoo",from=start_date,to=end_date)

price  = as.numeric(PBR$PBR.Adjusted)
chartSeries(PBR$PBR.Adjusted,name="Petrobras Prices",theme="white",TA = NULL)

return = diff(log(PBR$PBR.Adjusted))
chartSeries(return,name="Petrobras Returns",theme="white",TA = NULL) 

n = length(price)
return = (price[2:n]/price[1:(n-1)]-1)
n = length(return)
par(mfrow=c(2,2))
ts.plot(price)
ts.plot(return)
acf(price)
acf(return^2)



ndraw = 10000
svfit = svsample(return,draws=ndraw)
param = svfit$para[[1]][,1:3]
stdev = exp(svfit$latent[[1]]/2)
qmcmc = t(apply(stdev,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(2,2))
hist(param[,1],main="mu",prob=TRUE,xlab="")
hist(param[,2],main="phi",prob=TRUE,xlab="")
hist(param[,3],main="sigma",prob=TRUE,xlab="")
ts.plot(qmcmc[,2],ylab="Standard deviation")


par(mfrow=c(1,1))
ts.plot(qmcmc[,2],ylim=c(0,max(qmcmc[,2])),ylab="Standard deviation")
lines(0.1*abs(return),col=3,type="h")
legend("top",legend=c("Abs returns","tv-stdev"),col=c(3,1),lty=1,lwd=2,bty="n")
title("Petrobras - SV modeling")


n0 = 1500
par(mfrow=c(1,1))
plot((n-n0):n,qmcmc[(n-n0):n,2],ylim=c(-0.02,max(qmcmc[(n-n0):n,2])),
     ylab="Standard deviation",type="l",lwd=2,xlab="Day")
abline(h=sqrt(var(return[(n-n0):n])),lty=2,lwd=2,col=4)
abline(h=sqrt(var(return)),lty=4,lwd=2,col=4)
lines((n-n0):n,0.01+0.1*return[(n-n0):n],col=2)
legend("topright",legend=c("returns","tv-stdev","stdev (last 1500)","stdev (all)"),
       col=c(2,1,4,4),lty=c(1,1,2,4),lwd=2,bty="n")
title("Petrobras - SV modeling")



par.median = apply(param,2,median)
mu  = par.median[1]
phi = par.median[2]
sigma =par.median[3]

m0 = mu
C0 = sigma^2/(1-phi^2)


# Bootstrap filter
set.seed(12345)
M = 10000
h = rnorm(M,m0,sqrt(C0))
hs = matrix(0,M,n+1)
hs[,1] = h
for (t in 1:n){
  h1 = rnorm(M,mu*(1-phi)+phi*h,sigma)
  w = dnorm(return[t],0,exp(h1/2),log=TRUE)
  w = exp(w-max(w))
  h = sample(h1,size=M,replace=TRUE,prob=w)
  hs[,t+1] = h
}

qsdev.bf = t(apply(exp(hs/2),2,quantile,c(0.05,0.5,0.95)))


n0 = 1000
ts.plot(qsdev.bf[(n+1-n0):(n+1),2])
lines(qmcmc[(n-n0):n,2],col=2)