GJR-GARCH(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 GJR-GARCH(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)
mydata <- list( n = n, y = return);

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

param <-c("omega","alpha","beta", "mu", "gamma")
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.216   0.001 0.051  0.127 0.181 0.214 0.249 0.322  3918 1.001
## alpha 0.030   0.000 0.021  0.002 0.014 0.026 0.042 0.080  3970 1.000
## beta  0.669   0.001 0.057  0.555 0.631 0.669 0.709 0.777  3926 1.001
## mu    0.045   0.000 0.027 -0.008 0.027 0.045 0.064 0.099  4032 1.000
## gamma 0.244   0.001 0.045  0.163 0.212 0.242 0.272 0.339  3954 1.000
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun  2 12:36:48 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= 3870.377 BIC= 3896.3

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