GARCH-t(1,1) model

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 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, nu = 10)
}
fit  <- stan(file = 'garch-t.stan', data = mydata, init = myinit,
             chains = 4, warmup = 1000, iter = 11000, thin=10,
             cores = parallel::detectCores())

param <-c("omega","alpha","beta", "mu", "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.246   0.001 0.060 0.147 0.205 0.240  0.279  0.376  4027 1.000
## alpha 0.189   0.001 0.035 0.127 0.165 0.186  0.211  0.264  3773 1.001
## beta  0.616   0.001 0.063 0.477 0.578 0.621  0.659  0.725  3972 1.001
## mu    0.097   0.000 0.027 0.045 0.079 0.096  0.115  0.149  4105 1.000
## nu    9.074   0.039 2.317 5.795 7.434 8.632 10.195 14.893  3562 1.000
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun  2 13:57:37 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= 3882.444 BIC= 3908.368

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