GARCH Volatility Modelling

Libraries

library(tidyverse)
library(forecast)
library(gridExtra)
library(Metrics)
library(lmtest)
library(rugarch)

Dataset

df_raw <- read_csv("data_stock.csv")

df <- df_raw %>%
  mutate(Date = mdy(Date)) %>%        # "06/15/2026" -> Date
  arrange(Date) %>%                   # oldest -> newest
  mutate(
    return_log = (log(Price) - log(lag(Price))) * 100   # log return (%)
  )

# ts object for returns (drop the first NA row)
return_ts <- ts(na.omit(df$return_log))

n_total <- length(return_ts)
n_total

[1] 2517

Exploratory Analysis

p_price <- ggplot(df, aes(x = Date, y = Price)) +
  geom_line(color = "steelblue", linewidth = 0.7) +
  labs(title = "Stock Price", x = "Date", y = "Price") +
  theme_minimal()

p_ret <- ggplot(df %>% filter(!is.na(return_log)),
                aes(x = Date, y = return_log)) +
  geom_line(color = "darkred", linewidth = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  labs(title = "Log Return (%)", x = "Date", y = "Return (%)") +
  theme_minimal()

grid.arrange(p_price, p_ret, ncol = 1)

Train / Test Split

h <- 30
n <- length(return_ts)

train_ts <- subset(return_ts, end   = n - h)        # training
test_ts  <- subset(return_ts, start = n - h + 1)    # last 30 days

length(train_ts)

[1] 2487

length(test_ts)

[1] 30

ACF / PACF

grid.arrange(
  ggAcf(train_ts)  + labs(title = "ACF of Return"),
  ggPacf(train_ts) + labs(title = "PACF of Return"),
  ncol = 2
)

ARCH Effect Check

# Ljung-Box on squared (mean-removed) returns -> detect ARCH effect
res0 <- train_ts - mean(train_ts)
Box.test(res0^2, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  res0^2
X-squared = 484.62, df = 12, p-value < 2.2e-16

A small p-value indicates volatility clustering (ARCH effect), motivating a GARCH-type model.

ARIMA Modelling

Before modelling the variance, we model the conditional mean of the returns with an ARIMA process.

Order identification

auto.arima(train_ts)

Series: train_ts 
ARIMA(2,0,1) with zero mean 

Coefficients:
          ar1     ar2     ma1
      -0.9596  0.0067  0.9822
s.e.   0.0236  0.0204  0.0124

sigma^2 = 8.574:  log likelihood = -6199.4
AIC=12406.8   AICc=12406.82   BIC=12430.08

auto.arima selects ARIMA(2,0,1) with zero mean, consistent with the near-white-noise ACF/PACF of returns.

Fit ARIMA(2,0,1)

fit_arima <- Arima(train_ts, order = c(2, 0, 1), include.mean = FALSE)
fit_arima

Series: train_ts 
ARIMA(2,0,1) with zero mean 

Coefficients:
          ar1     ar2     ma1
      -0.9596  0.0067  0.9822
s.e.   0.0236  0.0204  0.0124

sigma^2 = 8.574:  log likelihood = -6199.4
AIC=12406.8   AICc=12406.82   BIC=12430.08

Parameter significance test

coeftest(fit_arima)   # from lmtest: estimate, std.error, z, p-value

z test of coefficients:

      Estimate Std. Error  z value Pr(>|z|)    
ar1 -0.9596042  0.0236139 -40.6372   <2e-16 ***
ar2  0.0067159  0.0204429   0.3285   0.7425    
ma1  0.9821589  0.0123988  79.2143   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual diagnostics

checkresiduals(fit_arima)

    Ljung-Box test

data:  Residuals from ARIMA(2,0,1) with zero mean
Q* = 6.7378, df = 7, p-value = 0.4567

Model df: 3.   Total lags used: 10

# explicit Ljung-Box on residuals -> no remaining autocorrelation expected
Box.test(residuals(fit_arima), lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  residuals(fit_arima)
X-squared = 10.998, df = 12, p-value = 0.5291

# Ljung-Box on SQUARED residuals -> tests for remaining ARCH effect
Box.test(residuals(fit_arima)^2, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  residuals(fit_arima)^2
X-squared = 481.1, df = 12, p-value < 2.2e-16

Mean forecast

fc_arima <- forecast(fit_arima, h = h, level = c(80, 95))

arima_df <- tibble(
  day      = seq_len(h),
  forecast = as.numeric(fc_arima$mean),
  actual   = as.numeric(test_ts),
  lo80 = as.numeric(fc_arima$lower[, 1]), hi80 = as.numeric(fc_arima$upper[, 1]),
  lo95 = as.numeric(fc_arima$lower[, 2]), hi95 = as.numeric(fc_arima$upper[, 2])
)

ggplot(arima_df, aes(x = day)) +
  geom_ribbon(aes(ymin = lo95, ymax = hi95, fill = "95% PI"), alpha = 0.15) +
  geom_ribbon(aes(ymin = lo80, ymax = hi80, fill = "80% PI"), alpha = 0.25) +
  geom_line(aes(y = forecast, color = "Forecast"), linewidth = 0.9) +
  geom_line(aes(y = actual, color = "Actual"), linewidth = 0.7) +
  geom_point(aes(y = actual, color = "Actual"), size = 1.3) +
  scale_color_manual(values = c("Forecast" = "red", "Actual" = "steelblue")) +
  scale_fill_manual(values  = c("95% PI" = "grey70", "80% PI" = "grey50")) +
  labs(title = "ARIMA(2,0,1) Forecast vs Actual (Test Set)",
       x = "Days Ahead", y = "Return (%)", color = NULL, fill = NULL) +
  theme_minimal() +
  theme(legend.position = "top")

# point-forecast accuracy of the MEAN
tibble(
  Metric = c("RMSE", "MAE"),
  Value  = c(Metrics::rmse(as.numeric(test_ts), as.numeric(fc_arima$mean)),
             Metrics::mae(as.numeric(test_ts),  as.numeric(fc_arima$mean)))
)

Modelling Function

# --- build the 30 candidate orders: ARCH(1..5) + GARCH(1..5, 1..5) ---
build_orders <- function() {
  arch_orders  <- lapply(1:5, function(q) c(q, 0))            # 5 ARCH models
  garch_orders <- list()
  for (p in 1:5) for (q in 1:5)                               # 25 GARCH models
    garch_orders[[length(garch_orders) + 1]] <- c(p, q)
  c(arch_orders, garch_orders)                                # 30 total
}

# --- build a single spec ---
make_spec <- function(order, arma, incmean) {
  ugarchspec(
    variance.model     = list(model = "sGARCH", garchOrder = order),
    mean.model         = list(armaOrder = arma, include.mean = incmean),
    distribution.model = "norm"
  )
}

# --- fit all 30 candidates, return ranked-by-BIC table ---
grid_select <- function(arma, incmean, data) {
  orders <- build_orders()
  map_dfr(orders, function(o) {
    spec <- make_spec(o, arma, incmean)
    fit  <- tryCatch(ugarchfit(spec, data, solver = "hybrid"),
                     error = function(e) NULL)
    if (is.null(fit) || convergence(fit) != 0) return(NULL)
    ic <- infocriteria(fit)
    label <- if (o[2] == 0) paste0("ARCH(", o[1], ")")
             else           paste0("GARCH(", o[1], ",", o[2], ")")
    tibble(model = label, arch = o[1], garch = o[2],
           AIC = ic[1], BIC = ic[2])
  }) %>% arrange(BIC)
}

# --- evaluate static k-step AND rolling 1-step forecasts ---
eval_forecasts <- function(spec, train, ret_full, test, h = 30) {
  fit <- ugarchfit(spec, train, solver = "hybrid")

  # (1) static k-step: one fit, forecast h steps ahead
  fc_static    <- ugarchforecast(fit, n.ahead = h)
  sigma_static <- as.numeric(sigma(fc_static))
  mean_static  <- as.numeric(fitted(fc_static))

  # (2) rolling 1-step: refit each step, forecast 1 step ahead
  roll <- ugarchroll(spec, data = ret_full, n.ahead = 1,
                     forecast.length = h, refit.every = 1,   # set 5 for speed
                     refit.window = "recursive", solver = "hybrid")
  rd <- as.data.frame(roll)
  sigma_roll <- rd$Sigma
  mean_roll  <- rd$Mu

  actual <- as.numeric(test)
  proxy  <- abs(actual - mean(train))   # fixed, method-independent vol proxy

  df_out <- tibble(
    day          = seq_len(h),
    actual       = actual,
    proxy        = proxy,
    sigma_static = sigma_static,
    sigma_roll   = sigma_roll,
    mean_static  = mean_static,
    mean_roll    = mean_roll
  )
  list(fit = fit, roll = roll, df = df_out)
}

# --- tidy metric table for both methods ---
metric_table <- function(df) {
  tibble(
    Method = c("Static k-step", "Rolling 1-step"),
    RMSE   = c(Metrics::rmse(df$proxy, df$sigma_static),
               Metrics::rmse(df$proxy, df$sigma_roll)),
    MAE    = c(Metrics::mae(df$proxy, df$sigma_static),
               Metrics::mae(df$proxy, df$sigma_roll))
  )
}

# --- side-by-side forecast plot ---
plot_forecasts <- function(df, unc_sigma, ttl) {
  ggplot(df, aes(x = day)) +
    geom_line(aes(y = proxy,        color = "Proxy |return|"),  linewidth = 0.6, alpha = 0.6) +
    geom_line(aes(y = sigma_static, color = "Static k-step"),   linewidth = 0.9) +
    geom_line(aes(y = sigma_roll,   color = "Rolling 1-step"),  linewidth = 0.9) +
    geom_hline(yintercept = unc_sigma, linetype = "dashed", color = "grey50") +
    scale_color_manual(values = c("Proxy |return|" = "grey60",
                                  "Static k-step"   = "red",
                                  "Rolling 1-step"  = "steelblue")) +
    labs(title = ttl,
         subtitle = paste0("Dashed grey = unconditional sigma = ", round(unc_sigma, 3)),
         x = "Days Ahead", y = "Volatility (%)", color = NULL) +
    theme_minimal() +
    theme(legend.position = "top")
}

---

GARCH Models

Mean model: ARFIMA(0,0,0) with no mean term; distribution "norm".

Model selection

sel_A <- grid_select(arma = c(0, 0), incmean = FALSE, data = train_ts)
sel_A

best_A       <- sel_A %>% slice(1)
best_order_A <- c(best_A$arch, best_A$garch)
best_A

spec_A <- make_spec(best_order_A, arma = c(0, 0), incmean = FALSE)
fit_A  <- ugarchfit(spec_A, train_ts, solver = "hybrid")
fit_A

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(0,0,0)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
omega   1.133767    0.255466   4.4380    9e-06
alpha1  0.096335    0.017114   5.6291    0e+00
beta1   0.767663    0.042404  18.1034    0e+00

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
omega   1.133767    0.375845   3.0166 0.002556
alpha1  0.096335    0.021800   4.4191 0.000010
beta1   0.767663    0.053602  14.3214 0.000000

LogLikelihood : -6085.917 

Information Criteria
------------------------------------
                   
Akaike       4.8966
Bayes        4.9036
Shibata      4.8966
Hannan-Quinn 4.8991

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                   0.005762  0.9395
Lag[2*(p+q)+(p+q)-1][2]  0.362349  0.7620
Lag[4*(p+q)+(p+q)-1][5]  1.462390  0.7492
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      2.865 0.09052
Lag[2*(p+q)+(p+q)-1][5]     3.243 0.36445
Lag[4*(p+q)+(p+q)-1][9]     5.939 0.30523
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.1998 0.500 2.000  0.6549
ARCH Lag[5]    0.8160 1.440 1.667  0.7882
ARCH Lag[7]    3.7645 2.315 1.543  0.3818

Nyblom stability test
------------------------------------
Joint Statistic:  0.5078
Individual Statistics:              
omega  0.04529
alpha1 0.11127
beta1  0.04574

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         0.846 1.01 1.35
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias            0.257 0.7972    
Negative Sign Bias   1.141 0.2541    
Positive Sign Bias   0.744 0.4569    
Joint Effect         2.807 0.4223    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     337.6    3.797e-60
2    30     461.9    1.891e-79
3    40     551.3    1.059e-91
4    50     771.8   1.275e-130


Elapsed time : 0.02546692 

Residual diagnostics

z_A <- as.numeric(residuals(fit_A, standardize = TRUE))

# (a) no remaining autocorrelation in standardized residuals
Box.test(z_A, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  z_A
X-squared = 9.9407, df = 12, p-value = 0.6212

# (b) no remaining ARCH effect in squared standardized residuals
Box.test(z_A^2, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  z_A^2
X-squared = 11.422, df = 12, p-value = 0.4931

grid.arrange(
  ggAcf(z_A)   + labs(title = "ACF of Std. Residuals"),
  ggAcf(z_A^2) + labs(title = "ACF of Squared Std. Residuals"),
  ncol = 2
)

z_df_A <- tibble(z = z_A)
p_hist_A <- ggplot(z_df_A, aes(x = z)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "steelblue", color = "white", alpha = 0.7) +
  geom_density(color = "darkred", linewidth = 0.9) +
  stat_function(fun = dnorm, args = list(mean = mean(z_A), sd = sd(z_A)),
                color = "black", linetype = "dashed", linewidth = 0.9) +
  labs(title = "Standardized Residuals vs Normal",
       x = "Standardized Residuals", y = "Density") +
  theme_minimal()

p_qq_A <- ggplot(z_df_A, aes(sample = z)) +
  stat_qq(color = "steelblue") +
  stat_qq_line(color = "darkred", linewidth = 0.8) +
  labs(title = "Normal Q-Q Plot", x = "Theoretical Quantiles",
       y = "Sample Quantiles") +
  theme_minimal()

grid.arrange(p_hist_A, p_qq_A, ncol = 2)

Volatility forecasting

res_A <- eval_forecasts(spec_A, train_ts, return_ts, test_ts, h = h)
fc_A  <- res_A$df

fc_A %>%
  transmute(day, sigma_static, sigma_roll, diff = sigma_roll - sigma_static)

plot_forecasts(fc_A, unc_sigma_A,
               "Static k-step vs Rolling 1-step Volatility")

# return with rolling +/-2 sigma band
ggplot(fc_A, aes(x = day)) +
  geom_ribbon(aes(ymin = mean_roll - 2*sigma_roll,
                  ymax = mean_roll + 2*sigma_roll, fill = "Rolling +/-2 sigma"),
              alpha = 0.25) +
  geom_line(aes(y = actual, color = "Actual return"), linewidth = 0.7) +
  geom_point(aes(y = actual, color = "Actual return"), size = 1.2) +
  scale_color_manual(values = c("Actual return" = "steelblue")) +
  scale_fill_manual(values = c("Rolling +/-2 sigma" = "grey50")) +
  labs(title = "Actual Return with Rolling Forecast Band",
       x = "Days Ahead", y = "Return (%)", color = NULL, fill = NULL) +
  theme_minimal() + theme(legend.position = "top")

metric_table(fc_A)

---

ARIMA + GARCH Models

Mean model: ARFIMA(2,0,1) with zero mean (armaOrder = c(2, 1), include.mean = FALSE); distribution "norm". The same 30 variance candidates are searched.

Model selection

sel_B <- grid_select(arma = c(2, 1), incmean = FALSE, data = train_ts)
sel_B

best_B       <- sel_B %>% slice(1)
best_order_B <- c(best_B$arch, best_B$garch)
best_B

spec_B <- make_spec(best_order_B, arma = c(2, 1), incmean = FALSE)
fit_B  <- ugarchfit(spec_B, train_ts, solver = "hybrid")
fit_B

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(2,0,1)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
ar1    -0.998549    0.019626  -50.8801 0.000000
ar2    -0.031564    0.018066   -1.7472 0.080610
ma1     0.980547    0.000203 4820.1011 0.000000
omega   1.151061    0.253890    4.5337 0.000006
alpha1  0.096574    0.017367    5.5606 0.000000
beta1   0.764926    0.042772   17.8840 0.000000

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
ar1    -0.998549    0.018557  -53.8110 0.000000
ar2    -0.031564    0.016615   -1.8998 0.057459
ma1     0.980547    0.000214 4579.6872 0.000000
omega   1.151061    0.387656    2.9693 0.002985
alpha1  0.096574    0.022127    4.3645 0.000013
beta1   0.764926    0.055833   13.7001 0.000000

LogLikelihood : -6081.903 

Information Criteria
------------------------------------
                   
Akaike       4.8958
Bayes        4.9098
Shibata      4.8958
Hannan-Quinn 4.9009

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.6588  0.4170
Lag[2*(p+q)+(p+q)-1][8]     2.1438  1.0000
Lag[4*(p+q)+(p+q)-1][14]    4.4432  0.9488
d.o.f=3
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      3.064 0.08005
Lag[2*(p+q)+(p+q)-1][5]     3.393 0.33990
Lag[4*(p+q)+(p+q)-1][9]     6.002 0.29800
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.1861 0.500 2.000  0.6662
ARCH Lag[5]    0.7353 1.440 1.667  0.8127
ARCH Lag[7]    3.6017 2.315 1.543  0.4077

Nyblom stability test
------------------------------------
Joint Statistic:  1.1641
Individual Statistics:              
ar1    0.16954
ar2    0.15249
ma1    0.28269
omega  0.04544
alpha1 0.11246
beta1  0.04497

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.49 1.68 2.12
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias          0.01845 0.9853    
Negative Sign Bias 0.98993 0.3223    
Positive Sign Bias 0.56995 0.5688    
Joint Effect       2.40082 0.4935    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     157.4    7.943e-24
2    30     166.1    3.647e-21
3    40     176.9    1.777e-19
4    50     187.1    5.212e-18


Elapsed time : 0.08000898 

Residual diagnostics

z_B <- as.numeric(residuals(fit_B, standardize = TRUE))
Box.test(z_B, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  z_B
X-squared = 8.1826, df = 12, p-value = 0.7707

Box.test(z_B^2, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  z_B^2
X-squared = 11.316, df = 12, p-value = 0.502

grid.arrange(
  ggAcf(z_B)   + labs(title = "ACF of Std. Residuals"),
  ggAcf(z_B^2) + labs(title = "ACF of Squared Std. Residuals"),
  ncol = 2
)

z_df_B <- tibble(z = z_B)
p_hist_B <- ggplot(z_df_B, aes(x = z)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "steelblue", color = "white", alpha = 0.7) +
  geom_density(color = "darkred", linewidth = 0.9) +
  stat_function(fun = dnorm, args = list(mean = mean(z_B), sd = sd(z_B)),
                color = "black", linetype = "dashed", linewidth = 0.9) +
  labs(title = "Standardized Residuals vs Normal",
       x = "Standardized Residuals", y = "Density") +
  theme_minimal()

p_qq_B <- ggplot(z_df_B, aes(sample = z)) +
  stat_qq(color = "steelblue") +
  stat_qq_line(color = "darkred", linewidth = 0.8) +
  labs(title = "Normal Q-Q Plot", x = "Theoretical Quantiles",
       y = "Sample Quantiles") +
  theme_minimal()

grid.arrange(p_hist_B, p_qq_B, ncol = 2)

Volatility forecasting

res_B <- eval_forecasts(spec_B, train_ts, return_ts, test_ts, h = h)
fc_B  <- res_B$df

fc_B %>%
  transmute(day, sigma_static, sigma_roll, diff = sigma_roll - sigma_static)

plot_forecasts(fc_B, unc_sigma_B,
               "Static k-step vs Rolling 1-step Volatility")

ggplot(fc_B, aes(x = day)) +
  geom_ribbon(aes(ymin = mean_roll - 2*sigma_roll,
                  ymax = mean_roll + 2*sigma_roll, fill = "Rolling +/-2 sigma"),
              alpha = 0.25) +
  geom_line(aes(y = actual, color = "Actual return"), linewidth = 0.7) +
  geom_point(aes(y = actual, color = "Actual return"), size = 1.2) +
  scale_color_manual(values = c("Actual return" = "steelblue")) +
  scale_fill_manual(values = c("Rolling +/-2 sigma" = "grey50")) +
  labs(title = "Actual Return with Rolling Forecast Band",
       x = "Days Ahead", y = "Return (%)", color = NULL, fill = NULL) +
  theme_minimal() + theme(legend.position = "top")

metric_table(fc_B)

---

Overall Comparison

bind_rows(
  metric_table(fc_A) %>% mutate(Mean_Model = "Zero-mean",      Best = best_A$model),
  metric_table(fc_B) %>% mutate(Mean_Model = "ARIMA(2,0,1)",   Best = best_B$model)
) %>%
  select(Mean_Model, Best, Method, RMSE, MAE)