SV example:

Stochastic Volatility Model without Leverage

Install packages

# install.packages("freqdom")
# install.packages("RcppProgress")
# install.packages("ggplot2")
# install.packages("ASV")
library(freqdom)
## Loading required package: mvtnorm
## Loading required package: matrixcalc
## 
## Attaching package: 'freqdom'
## The following object is masked from 'package:base':
## 
##     %*%
library(RcppProgress)
library(ggplot2)
library(ASV)

Generate data

set.seed(1111)
nobs = 1000;
mu = 0; phi = 0.97; sigma_eta = 0.3;
h = 0; Y = c();
for(i in 1:nobs){
  eps = rnorm(1, 0, 1)
  eta = rnorm(1, 0, sigma_eta)
  y   = eps * exp(0.5*h)
  h   = mu + phi * (h-mu) + eta
  Y   = append(Y, y)
}

Markov chain Monte Carlo for SV model without leverage

  • sv_mcmc(return vecor, # of iterations, # of burn-in period)
nsim = 5000; burn = 500;
set.seed(123)
out  = sv_mcmc(Y, nsim, burn);
## Start MCMC estimation
## Acceptance Rate for h     =  94.94 % 
## Acceptance Rate for (mu,phi,sigma_eta) =  70.20 %
vmu  = out[[1]]; vphi = out[[2]]; vsigma_eta =  out[[3]]; 
mh   = out[[4]];

Summary statistics, sample paths, sample autocorrelations and posterior densities:

ReportMCMC(cbind(vmu,vphi,vsigma_eta), vname=c(expression(mu), expression(phi),expression(sigma[eta])))
##                Mean  Std Dev     95%L   Median    95%U
## mu         0.034743 0.338016 -0.62041 0.038435 0.69592
## phi        0.970053 0.010297  0.94711 0.971302 0.98625
## sigma[eta] 0.313896 0.038177  0.24979 0.311170 0.39502
##            ESS   IF    CD  Pr(+)
## mu         949  5.3 0.854 0.5468
## phi        151 33.0 0.187 1.0000
## sigma[eta]  85 59.0 0.273 1.0000

ReportMCMC(cbind(mh[,250], mh[,500],mh[,750]), vname=c(expression(h[250]), expression(h[500]),expression(h[750])))
##            Mean Std Dev     95%L   Median     95%U
## h[250]  0.16693 0.50251 -0.79928  0.15632  1.18182
## h[500]  0.85477 0.52121 -0.14855  0.85210  1.89444
## h[750] -1.57330 0.56818 -2.67989 -1.57598 -0.41358
##         ESS  IF    CD  Pr(+)
## h[250] 1374 3.6 0.404 0.6312
## h[500]  880 5.7 0.101 0.9542
## h[750]  789 6.3 0.389 0.0034

95% credible intervals for log volatilities

mh_ci = matrix(0, nrow = nobs, ncol = 3)
for(i in 1:nobs){
  mh_ci[i,1:3] = t(quantile(mh[1:nsim,i],c(0.025,0.5, 0.975)))}

mydate = seq(1,nobs);
ggplot()+geom_line(aes(x=mydate, y=mh_ci[,1], color="95%L"),size=1)+geom_line(aes(x=mydate, y=mh_ci[,2], color="Median"),size=1)+geom_line(aes(x=mydate, y=mh_ci[,3], color="95%U"),size=1)+labs(x="Date", y="Log Volatility");

Log likelihood (given mu, phi, sigma_eta)

  • simple particle filter
mu  = mean(vmu); phi = mean(vphi); sigma_eta = mean(vsigma_eta);
# Number of Particles
I = 5000
# repeat m times to compute the standard error of the estimate
m  = 10; loglik = rep(0, m)
#
for(i in 1:m){
  loglik[i] = sv_pf(mu, phi, sigma_eta, Y, I)
}
print(loglik, digits = 6)
##  [1] -1521.54 -1520.91 -1521.39 -1521.25 -1521.51 -1521.82 -1522.02 -1521.23
##  [9] -1521.33 -1521.53
#
# Estimate and its standard error
loglik_m = mean(loglik); loglik_se = sqrt(var(loglik)/m);
result = c(loglik_m, loglik_se)
names(result) = c("Log likelihood", "Std Err")
print(result, digits = 6)
## Log likelihood        Std Err 
##   -1521.451719       0.099021
  • auxiliary particle filter
for(i in 1:m){
  loglik[i] = sv_apf(mu, phi, sigma_eta, Y, I)
}
print(loglik, digits = 6)
##  [1] -1520.15 -1521.09 -1520.38 -1521.55 -1521.40 -1521.31 -1521.18 -1521.53
##  [9] -1521.54 -1521.05
#
# Estimate and its standard error
loglik_m = mean(loglik); loglik_se = sqrt(var(loglik)/m);
result = c(loglik_m, loglik_se)
names(result) = c("Log likelihood", "Std Err")
print(result, digits = 6)
## Log likelihood        Std Err 
##    -1521.11844        0.15489

Log prior density

vhyper     = c(0, 1, 20, 1.5, 5, 0.05)
theta_star = c(mu, phi, sigma_eta) 
logprior   = sv_prior(theta_star, vhyper)
#
names(logprior) = c("Log prior")
print(logprior, digits = 4)
## Log prior 
##    -1.017

Log posterior density

h            = mh[nsim,]
theta        = c(vmu[nsim],vphi[nsim],vsigma_eta[nsim])
logposterior =  sv_posterior(h, theta, theta_star, Y, 5000, vhyper)
## Computing the log posterior density ...
#
# Estimate and its standard error
logpost_m = logposterior[1]; logpost_se = logposterior[2];
result = c(logpost_m, logpost_se)
names(result) = c("Log posterior", "Std Err")
print(result, digits = 3)
## Log posterior       Std Err 
##        6.6564        0.0327

Log marginal likelihood

#
result = c( loglik_m + logprior - logpost_m, sqrt(loglik_se^2+logpost_se^2))
names(result) = c("Log marginal likelihood", "Std Err")
print(result, digits = 4)
## Log marginal likelihood                 Std Err 
##              -1528.7919                  0.1583
# Or
result = sv_logML(h, theta, theta_star, Y, 1000, 5000, vHyper = vhyper)
## Computing the log likelihood ...
## Log likelihood ordinate (se) = -1521.9235 (0.3523)
## Log prior ordinate = -1.0170
## Computing the log posterior density ...
## Log posterior ordinate (se) = 6.9342 (0.0279)
## Log marginal likelihood (se) = -1529.8747 (0.3534)
result1     = matrix(0, 4, 2)
result1[,1] =result[[1]]
result1[,2] =result[[2]]
  
colnames(result1) = c("Estimate", "Std Err")
rownames(result1) = c("Log marginal lik", "Log likelihood", "Log prior", "Log posterior")
print(result1, digit=4)
##                   Estimate Std Err
## Log marginal lik -1529.875 0.35342
## Log likelihood   -1521.924 0.35232
## Log prior           -1.017 0.00000
## Log posterior        6.934 0.02793