ASV example:

Stochastic Volatility Model with 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; rho = -0.5;
h = 0; Y = c();
for(i in 1:nobs){
  eps = rnorm(1, 0, 1)
  eta = rho*sigma_eta*eps + sigma_eta*sqrt(1-rho^2)*rnorm(1, 0, 1)
  y   = eps * exp(0.5*h)
  h   = mu + phi * (h-mu) + eta
  Y   = append(Y, y)
}

Markov chain Monte Carlo for SV model with leverage

  • asv_mcmc(return vecor, # of iterations, # of burn-in period)
nsim = 10000; burn = 3000;
set.seed(111)
out  = asv_mcmc(Y, nsim, burn);
## Start MCMC estimation
## Acceptance Rate for h     =  87.24 %
## Acceptance Rate for (mu,phi,sigma_eta,rho) =  73.75 %
vmu  = out[[1]]; vphi = out[[2]]; vsigma_eta =  out[[3]]; 
vrho = out[[4]]; mh   = out[[5]];

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

ReportMCMC(cbind(vmu,vphi,vsigma_eta, vrho), vname=c(expression(mu), expression(phi),expression(sigma[eta]), expression(rho)))
##                Mean  Std Dev     95%L   Median     95%U
## mu          0.13811 0.262104 -0.39621  0.14190  0.65962
## phi         0.96255 0.011735  0.93617  0.96375  0.98198
## sigma[eta]  0.30806 0.035997  0.24299  0.30500  0.38242
## rho        -0.33306 0.092972 -0.50629 -0.33680 -0.13807
##             ESS   IF    CD  Pr(+)
## mu         1652  6.1 0.064 0.7137
## phi         294 34.0 0.260 1.0000
## sigma[eta]  134 74.0 0.167 1.0000
## rho         250 40.0 0.223 0.0001

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.51912 0.46217 -0.37937  0.52156  1.44318
## h[500]  0.86652 0.54822 -0.21812  0.86229  1.95442
## h[750] -1.47188 0.57892 -2.59901 -1.47711 -0.32669
##         ESS  IF    CD  Pr(+)
## h[250] 2636 3.8 0.744 0.8734
## h[500] 1138 8.8 0.859 0.9432
## h[750] 1384 7.2 0.895 0.0069

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 and rho)

  • simple particle filter
mu  = mean(vmu); phi = mean(vphi); sigma_eta = mean(vsigma_eta);
rho = mean(vrho);
# 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] = asv_pf(mu, phi, sigma_eta, rho, Y, I)
}
print(loglik, digits = 6)
##  [1] -1632.19 -1632.64 -1632.12 -1632.65 -1632.66 -1632.30 -1631.74 -1632.42
##  [9] -1632.19 -1632.20
#
# 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 
##   -1.63231e+03    9.24541e-02
  • auxiliary particle filter
for(i in 1:m){
  loglik[i] = asv_apf(mu, phi, sigma_eta, rho, Y, I)
}
print(loglik, digits = 6)
##  [1] -1631.82 -1632.58 -1633.10 -1632.20 -1632.09 -1632.44 -1632.54 -1633.36
##  [9] -1632.85 -1633.03
#
# 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 
##   -1632.601268       0.153681

Log prior density

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

Log posterior density

h            = mh[nsim,]
theta        = c(vmu[nsim],vphi[nsim],vsigma_eta[nsim],vrho[nsim])
logposterior =  asv_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 
##        8.3475        0.0532

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 
##              -1642.5068                  0.1626
# Or
result = asv_logML(h, theta, theta_star, Y, 1000, 5000, vHyper = vhyper)
## Computing the log likelihood ...
## Log likelihood ordinate (se) = -1632.8819 (0.3966)
## Log prior ordinate = -1.5580
## Computing the log posterior density ...
## Log posterior ordinate (se) = 8.3153 (0.0533)
## Log marginal likelihood (se) = -1642.7553 (0.4002)
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 -1642.755 0.40018
## Log likelihood   -1632.882 0.39661
## Log prior           -1.558 0.00000
## Log posterior        8.315 0.05333