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)

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 without leverage

  • sv_mcmc(return vecor, # of iterations, # of burn-in period)
nobs <- n; Y <- return
nsim = 15000; burn = 500;
set.seed(123)
out  = sv_mcmc(Y, nsim, burn);
## Start MCMC estimation
## Acceptance Rate for h     =  92.23 % 
## Acceptance Rate for (mu,phi,sigma_eta) =  86.17 %
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.042269 0.091708 -0.22478 -0.041962 0.13475
## phi         0.869863 0.033055  0.79902  0.872624 0.92512
## sigma[eta]  0.362822 0.050808  0.27313  0.361754 0.46804
##             ESS   IF    CD  Pr(+)
## mu         1946  7.7 0.099 0.3199
## phi         211 71.0 0.054 1.0000
## sigma[eta]  151 99.0 0.336 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.36252 0.48404 -0.55227  0.34813 1.35150
## h[500] -0.24167 0.48289 -1.14594 -0.25577 0.74664
## h[750] -0.17225 0.47317 -1.04805 -0.18947 0.79728
##         ESS  IF    CD  Pr(+)
## h[250] 5449 2.8 0.763 0.7717
## h[500] 2443 6.1 0.215 0.2976
## h[750] 2784 5.4 0.085 0.3474

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)

  • 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] -1947.52 -1946.62 -1947.94 -1946.66 -1947.09 -1946.30 -1947.30 -1946.84
##  [9] -1947.84 -1947.01
#
# 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 
##   -1947.111275       0.170071
  • auxiliary particle filter
for(i in 1:m){
  loglik[i] = sv_apf(mu, phi, sigma_eta, Y, I)
}
print(loglik, digits = 6)
##  [1] -1946.69 -1946.27 -1946.77 -1946.77 -1947.14 -1946.43 -1946.95 -1946.61
##  [9] -1947.19 -1946.24
#
# 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 
##   -1946.706347       0.104439

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 
##    -2.225

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.2064        0.0481

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 
##               -1955.137                   0.115
# Or
result = sv_logML(h, theta, theta_star, Y, 1000, 5000, vHyper = vhyper)
## Computing the log likelihood ...
## Log likelihood ordinate (se) = -1946.7848 (0.4883)
## Log prior ordinate = -2.2246
## Computing the log posterior density ...
## Log posterior ordinate (se) = 6.2237 (0.0488)
## Log marginal likelihood (se) = -1955.2330 (0.4908)
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 -1955.233 0.49075
## Log likelihood   -1946.785 0.48832
## Log prior           -2.225 0.00000
## Log posterior        6.224 0.04878