set.seed(314)

Priors

Truncated Normal Prior

hist(((vix$CLOSE / 100)^2)/252, col = "lightblue",breaks = 50, main = "Realized Volatility of VIX Index", xlab = "Daily ReaL Vol")

z = rtruncnorm(10000, a = 0, b = Inf, mean = 0, sd = .003)
hist(z, col = "lightblue", main = "Truncated Normal Prior Distribution w/ Mu = 0 & SD = .003", xlab = "Logarithmic Returns Modeling", breaks = 50)

# Gamma Prior

y = rgamma(10000, 2, 100)
hist(y, col = "lightblue", main = "Gamma Prior Distribution w/ Alpha = 2 & Beta = 100", xlab = "Logarithmic Returns Modeling", breaks = 50)

# Exponential Prior

x = rexp(10000,500)
hist(x, col = "lightblue", main = "Exponential Prior Distribution w/ Rate Parameter = 500", xlab = "Logarithmic Returns Modeling", breaks = 50)

Likelihood Distribution

hist(spys$LogReturns, breaks = 50, col = "lightblue", main = "SPY Daily Logarithmic Returns", xlab = "Daily Log Returns")

shapiro.test(spys$LogReturns)
## 
##  Shapiro-Wilk normality test
## 
## data:  spys$LogReturns
## W = 0.83195, p-value < 2.2e-16

Posterior Distribution

Truncated Normal Prior

modelStringN = "model {
#Priors
sigmasq ~ dnorm(0,100)T(0,)  #Truncated Normal Prior for realized vol
b ~ dunif(-1,1)
mu ~ dnorm(0,.1)
inv = 1 / sigmasq
r0 ~ dnorm(mu, inv * (1-b * b))


#Likelihood
r[1] ~ dnorm(b * r0 + mu, inv);  #Autoregressive AR(1) 
for(i in 2:n) {   
  r[i] ~ dnorm(mu + (b * r[i-1]), inv);
   }

}" 

for (i in 31:length(spys$LogReturns)){
  datas = list("r" = spys$LogReturns[i-30:i], "n" = length(spys$LogReturns[i-30:i])) #last 30 days trailing
  model = jags.model(file = textConnection(modelStringN), data = datas, n.chains = 2)
 


}
update(model, 1000, progress.bar = "none")

Nrep = 1000

posterior_sampleN <- coda.samples(model, 
                                 variable.names = c("sigmasq"), 
                                 n.iter = Nrep,
                                 progress.bar = "none")



realvolN = 100*sqrt(252 * as.numeric(unlist(posterior_sampleN)))   
plotPost(realvolN, main = "Posterior Distribution of Realized Volatility w/ Exponential Prior", xlab = "Annual Realized Volatility Levels")

diagMCMC(posterior_sampleN)
quantile(realvolN, c(.25,.75))

quantile(realvolN, c(.10,.9))
quantile(realvolN, c(.01,.99))

#Posterior Predictive Distribution for Truncated Normal Prior

#Posterior Predictive Distribution
theta_sim = as.matrix(posterior_sampleN)
sim = 100 * sqrt(252 *theta_sim)
y_sim = rtruncnorm(nrow(sim),a = 0, b = Inf,  mean(realvolN), sim)

hist(y_sim, freq = FALSE, col= "lightblue", xlab = "Realized Volatility (Annum)",
     main = "Posterior Preditive Distribution for Realized Volatility")
lines(density(y_sim))
abline(v = quantile(y_sim, c(0.025, 0.975)), col = "orange")

quantile(y_sim, c(.025,.975))

Gamma Prior

modelStringG = "model {
#Priors
sigmasq ~ dgamma(2,1000) #Gamma Prior   #Truncated Normal Prior for realized vol
b ~ dunif(-1,1)
mu ~ dnorm(0,.1)
inv = 1 / sigmasq
r0 ~ dnorm(mu, inv * (1-b * b))


#Likelihood
r[1] ~ dnorm(b * r0 + mu, inv);  #Autoregressive AR(1) 
for(i in 2:n) {   
  r[i] ~ dnorm(mu + (b * r[i-1]), inv);
   }

}" 

for (i in 31:length(spys$LogReturns)){
  datas = list("r" = spys$LogReturns[i-30:i], "n" = length(spys$LogReturns[i-30:i])) #last 30 days trailing
  model = jags.model(file = textConnection(modelStringG), data = datas, n.chains = 2)
 


}
update(model, 1000, progress.bar = "none")

Nrep = 1000

posterior_sampleG <- coda.samples(model, 
                                 variable.names = c("sigmasq"), 
                                 n.iter = Nrep,
                                 progress.bar = "none")



realvolG = 100 * sqrt(252 *  as.numeric(unlist(posterior_sampleG)))   

plotPost(realvolG, main = "Posterior Distribution of Realized Volatility w/ Exponential Prior", xlab = "Annual Realized Volatility Levels")

diagMCMC(posterior_sampleG)
quantile(realvolG, c(.25,.75))

quantile(realvolG, c(.10,.9))
quantile(realvolG, c(.01,.99))

Exponential Prior

modelStringE = "model {
#Priors
sigmasq ~ dexp(500) #Exponential Prior 
b ~ dunif(-1,1)
mu ~ dnorm(0,.1)
inv = 1 / sigmasq
r0 ~ dnorm(mu, inv * (1-b * b))


#Likelihood
r[1] ~ dnorm(b * r0 + mu, inv);  #Autoregressive AR(1) 
for(i in 2:n) {   
  r[i] ~ dnorm(mu + (b * r[i-1]), inv);
   }

}" 
for (i in 31:length(spys$LogReturns)){
  datas = list("r" = spys$LogReturns[i-30:i], "n" = length(spys$LogReturns[i-30:i])) #last 30 days trailing
  model = jags.model(file = textConnection(modelStringE), data = datas, n.chains = 2)
 


}

update(model, 1000, progress.bar = "none")
Nrep = 1000


posterior_sampleE <- coda.samples(model, 
                                 variable.names = c("sigmasq"), 
                                 n.iter = Nrep,
                                 progress.bar = "none")



realvolE = 100*sqrt(252 *  as.numeric(unlist(posterior_sampleE)))   
plotPost(realvolE, main = "Posterior Distribution of Realized Volatility w/ Exponential Prior", xlab = "Annual Realized Volatility Levels")

diagMCMC(posterior_sampleE)
quantile(realvolE, c(.25,.75))

quantile(realvolE, c(.10,.9))
quantile(realvolE, c(.01,.99))
#Compare realized vol with implied vol data, no need to make posterior of imp vol just actual data I think
imp = vix$CLOSE
x <- imp  ## data forming your empirical distribution
ll <- min(realvolN)   ## lower bound of interval of interest
ul <- max(realvolN)     ## upper bound of interval of interest

sum(x > ll & x < ul)/length(x)


hist(x, col = "lightblue", main = "Implied Volatility Distribution (VIX)", xlab = "Annual Implied Volatility of SPY")