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")

Fit GJR-GARCH-t(1,1) model using rstan

library(rstan)
## Loading required package: StanHeaders
## 
## rstan version 2.32.7 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
library(rstudioapi)
## Warning: package 'rstudioapi' was built under R version 4.5.3
mydata <- list( n = n, y = return);

myinit <- function() { 
  list( omega = 0.1, alpha = 0.1, beta = 0.1, mu = 0.0, gamma = 0.1, nu = 10)
}
fit  <- stan(file = 'gjr-garch-t.stan', data = mydata, init = myinit,
             chains = 4, warmup = 1000, iter = 11000, thin=10, cores = parallel::detectCores())

param <-c("omega","alpha","beta", "mu", "gamma", "nu")
print(fit, pars = param, digits=3)
## Inference for Stan model: anon_model.
## 4 chains, each with iter=11000; warmup=1000; thin=10; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean    sd  2.5%   25%    50%    75%  97.5% n_eff  Rhat
## omega  0.221   0.001 0.052 0.133 0.183  0.215  0.252  0.335  3975 1.000
## alpha  0.035   0.000 0.025 0.002 0.015  0.030  0.050  0.097  3872 1.000
## beta   0.664   0.001 0.058 0.542 0.629  0.668  0.705  0.766  3776 1.000
## mu     0.066   0.000 0.027 0.012 0.048  0.066  0.085  0.120  4188 1.000
## gamma  0.232   0.001 0.050 0.143 0.197  0.230  0.263  0.342  4213 1.000
## nu    11.663   0.059 3.671 6.801 9.093 10.915 13.370 20.748  3835 0.999
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun  2 14:08:17 2026.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
stan_dens(fit,  pars = param, nrow = 3, ncol = 2, separate_chains = TRUE)

stan_trace(fit, pars = param, nrow = 3, ncol = 2)

stan_ac(fit, pars = param, nrow = 3, ncol = 2, separate_chains = TRUE)

mcmcout <- rstan::extract(fit)
nparam  <- length(param)
cat("AIC=",-2*max(mcmcout$loglik)+2*nparam,"BIC=",-2*max(mcmcout$loglik)+nparam*log(n))
## AIC= 3855.85 BIC= 3886.957

Plot the median of conditional variances

v_med = rep(0, n);
for(t in 1:n){
  v_med[t] = quantile(mcmcout$sigma2[,t], 0.5);
}
variance.df <-data.frame(dates = mydate, med = v_med)  
#
ggplot(variance.df) + geom_line(aes(x=dates, y=v_med)) + 
  xlab("Date") + 
  ylab("Estimated Conditional Variance") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y")