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)

Read Topix Data 2021/1/4-2026/5/29

library(readr)
topix <- read_csv("topix.csv")
## Rows: 1320 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): date
## dbl (1): close
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
topix$price <- as.numeric(topix$close)
topix$Date  <- as.Date(topix$date)
n <- length(topix$price)

return <- 100 * diff(log(topix$price))
mydate <- topix$Date[2:n]
n <- n-1

Draw time series plot of TOPIX returns

topix.df <-data.frame(dates = mydate, value = return) 
library(ggplot2)
ggplot(topix.df, aes(x = mydate, y = return)) +
  geom_line(linewidth = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
  geom_hline(yintercept = c(2, -2), linetype = "dotted", linewidth = 0.7) +
  labs(title = "TOPIX Daily Log Returns", 
       subtitle = "Dotted lines indicate ±2%",
       x = "Date", y = "Log return (%)") + 
  scale_x_date(date_breaks = "1 year", date_labels = "%Y")

Markov chain Monte Carlo for SV model with leverage

  • asv_mcmc(return vecor, # of iterations, # of burn-in period)
nobs <- n; Y <- return
nsim = 20000; burn = 3000;
set.seed(111)
out  = asv_mcmc(Y, nsim, burn);
## Start MCMC estimation
## Acceptance Rate for h     =  74.69 %
## Acceptance Rate for (mu,phi,sigma_eta,rho) =  86.49 %
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.069267 0.086218 -0.099114  0.068848  0.23808
## phi         0.895494 0.024850  0.841165  0.897655  0.93778
## sigma[eta]  0.301493 0.039250  0.228218  0.300041  0.37852
## rho        -0.516160 0.085674 -0.674827 -0.518709 -0.34118
##             ESS  IF    CD  Pr(+)
## mu         1422  14 0.570 0.7969
## phi         254  79 0.483 1.0000
## sigma[eta]  170 118 0.669 1.0000
## rho         252  79 0.360 0.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.45329 0.38498 -0.27438  0.43937 1.24573
## h[500] -0.30818 0.37464 -1.01379 -0.31966 0.45014
## h[750] -0.32325 0.40933 -1.09178 -0.33508 0.50706
##         ESS  IF    CD  Pr(+)
## h[250] 3975 5.0 0.185 0.8839
## h[500] 3119 6.4 0.532 0.1987
## h[750] 5035 4.0 0.885 0.2082

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");
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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] -1931.42 -1930.22 -1931.62 -1931.36 -1930.80 -1931.35 -1930.96 -1930.76
##  [9] -1931.24 -1931.06
#
# 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 
##   -1931.078127       0.130472
  • auxiliary particle filter
for(i in 1:m){
  loglik[i] = asv_apf(mu, phi, sigma_eta, rho, Y, I)
}
print(loglik, digits = 6)
##  [1] -1930.75 -1931.39 -1930.67 -1931.36 -1930.94 -1931.18 -1931.80 -1931.11
##  [9] -1931.39 -1931.35
#
# 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 
##   -1931.193829       0.107419

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.559

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.9768        0.0782

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 
##              -1941.7298                  0.1329
# Or
result = asv_logML(h, theta, theta_star, Y, 1000, 5000, vHyper = vhyper)
## Computing the log likelihood ...
## Log likelihood ordinate (se) = -1931.5999 (0.2440)
## Log prior ordinate = -1.5592
## Computing the log posterior density ...
## Log posterior ordinate (se) = 9.2654 (0.0666)
## Log marginal likelihood (se) = -1942.4245 (0.2529)
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 -1942.424 0.25295
## Log likelihood   -1931.600 0.24402
## Log prior           -1.559 0.00000
## Log posterior        9.265 0.06659