Ioakeim Kyrgiafinis 1.0 2026-05-22 Rolling Garch (:2, :2) with confidence bands C58 rolling_garch_p_q Version 1.0 Estimates a rolling GARCH(p,q) model over a fixed window and tracks parameter instability over time, with 95% confidence intervals via the delta method. ARGUMENTS: y series Return series to model window int Rolling window size (minimum 250 recommended) arch_order int ARCH order p: 1 or 2 garch_order int GARCH order q: 1 or 2 Mean_regressors list Additional mean equation regressors (pass null for none) include_const bool Include constant in mean equation (1 = yes, 0 = no) RETURNS: Bundle R containing: Coefficients (sparse): roll_alpha0 Rolling variance intercept (omega) roll_alpha1/2 Individual ARCH coefficients (if arch_order >= 1/2) roll_beta1/2 Individual GARCH coefficients (if garch_order >= 1/2) roll_alpha Sum of ARCH coefficients roll_beta Sum of GARCH coefficients Variance (sparse): roll_sigma2 One-step-ahead conditional variance forecast roll_sigma One-step-ahead conditional volatility forecast Interpolated for plotting: roll_alpha0_plot Interpolated omega roll_alpha1/2_plot roll_beta1/2_plot roll_alpha_plot Interpolated sum of ARCH coefficients roll_beta_plot Interpolated sum of GARCH coefficients ab_plot Interpolated alpha + beta (persistence) roll_sigma2_plot Interpolated variance roll_sigma_plot Interpolated volatility Information criteria (sparse and interpolated): roll_aic/bic/logL roll_aic/bic/logL_plot 95% CI bundles: CI_alpha0 CI for omega CI_alpha1/2 CI for individual ARCH coefficients CI_beta1/2 CI for individual GARCH coefficients CI_alpha CI for sum of ARCH coefficients (if arch_order >= 2) CI_beta CI for sum of GARCH coefficients (if garch_order >= 2) CI_ab CI for persistence (alpha + beta) NOTES: - Windows where estimation fails are stored as NA and interpolated for plots - Persistence (ab_plot) >= 1 indicates non-stationarity in that window - GARCH(2,2) may show instability due to overparameterization - Minimum recommended window: 250 obs for daily data - CI widths exceeding 5 are trimmed to NA to suppress extreme values bundle R scalar ci_trim = 5 scalar WINDOW = window scalar T = $nobs scalar nroll = T - WINDOW # Storage series series roll_alpha0 = NA series roll_alpha1 = NA series roll_alpha2 = NA series roll_beta1 = NA series roll_beta2 = NA series roll_alpha = NA series roll_beta = NA series roll_sigma2 = NA series roll_sigma = NA series width_ab = NA series width_alpha = NA series width_beta = NA series width_alpha0 = NA series width_alpha1 = NA series width_alpha2 = NA series width_beta1 = NA series width_beta2 = NA series roll_aic = NA series roll_bic = NA series roll_logL = NA series h_series = NA series e_series = NA set verbose off if include_const Mean_regressors -= const endif # ROLLING LOOP loop i = 1 .. nroll --quiet scalar t_start = i scalar t_end = i + WINDOW - 1 scalar pos = t_end + 1 smpl t_start t_end # ESTIMATION if include_const if nelem(Mean_regressors) > 0 catch garch arch_order garch_order ; y const Mean_regressors else catch garch arch_order garch_order ; y const endif else if nelem(Mean_regressors) > 0 catch garch arch_order garch_order --nc ; y Mean_regressors else catch garch arch_order garch_order --nc ; y endif endif scalar err = $error if err == 0 h_series = $h e_series = $uhat scalar ncoeff = nelem($coeff) # dimensions scalar n_garch = 1 + arch_order + garch_order scalar n_mean_actual = ncoeff - n_garch # coefficient indices scalar idx_alpha0 = n_mean_actual + 1 scalar idx_alpha1 = n_mean_actual + 2 scalar idx_alpha2 = n_mean_actual + 3 scalar idx_beta1 = n_mean_actual + 1 + arch_order + 1 scalar idx_beta2 = n_mean_actual + 1 + arch_order + 2 # omega scalar o = NA if idx_alpha0 <= ncoeff o = $coeff[idx_alpha0] endif # one-step variance forecast scalar sigma2 = o loop j = 1 .. arch_order --quiet scalar idx_a = n_mean_actual + 1 + j if idx_a <= ncoeff sigma2 += $coeff[idx_a] * e_series[t_end + 1 - j]^2 endif endloop loop j = 1 .. garch_order --quiet scalar idx_b = n_mean_actual + 1 + arch_order + j if idx_b <= ncoeff sigma2 += $coeff[idx_b] * h_series[t_end + 1 - j] endif endloop scalar sigma = NA if !missing(sigma2) && sigma2 > 0 sigma = sqrt(sigma2) else sigma2 = NA endif # sums scalar sum_alpha = 0 loop j = 1 .. arch_order --quiet scalar idx_a = n_mean_actual + 1 + j if idx_a <= ncoeff sum_alpha += $coeff[idx_a] endif endloop scalar sum_beta = 0 loop j = 1 .. garch_order --quiet scalar idx_b = n_mean_actual + 1 + arch_order + j if idx_b <= ncoeff sum_beta += $coeff[idx_b] endif endloop # STORE INDIVIDUAL COEFFICIENTS if arch_order >= 1 && idx_alpha1 <= ncoeff roll_alpha1[pos] = $coeff[idx_alpha1] endif if arch_order >= 2 && idx_alpha2 <= ncoeff roll_alpha2[pos] = $coeff[idx_alpha2] endif if garch_order >= 1 && idx_beta1 <= ncoeff roll_beta1[pos] = $coeff[idx_beta1] endif if garch_order >= 2 && idx_beta2 <= ncoeff roll_beta2[pos] = $coeff[idx_beta2] endif # DELTA METHOD matrix V = $vcv scalar n_params_actual = rows(V) #alpha0 CI matrix grad_alpha0 = zeros(1, n_params_actual) if idx_alpha0 <= n_params_actual grad_alpha0[1, idx_alpha0] = 1 endif scalar var_alpha0 = grad_alpha0 * V * grad_alpha0' scalar se_alpha0 = (var_alpha0 > 0) ? sqrt(var_alpha0) : NA # a+b matrix ab = zeros(1, n_params_actual) loop j = 1 .. arch_order --quiet scalar idx_a = n_mean_actual + 1 + j if idx_a <= n_params_actual ab[1, idx_a] = 1 endif endloop loop j = 1 .. garch_order --quiet scalar idx_b = n_mean_actual + 1 + arch_order + j if idx_b <= n_params_actual ab[1, idx_b] = 1 endif endloop scalar var_ab = ab * V * ab' scalar se_ab = (var_ab > 0) ? sqrt(var_ab) : NA # alpha (sum) matrix grad_alpha = zeros(1, n_params_actual) loop j = 1 .. arch_order --quiet scalar idx_a = n_mean_actual + 1 + j if idx_a <= n_params_actual grad_alpha[1, idx_a] = 1 endif endloop scalar var_alpha = grad_alpha * V * grad_alpha' scalar se_alpha = (var_alpha > 0) ? sqrt(var_alpha) : NA # beta (sum) matrix grad_beta = zeros(1, n_params_actual) loop j = 1 .. garch_order --quiet scalar idx_b = n_mean_actual + 1 + arch_order + j if idx_b <= n_params_actual grad_beta[1, idx_b] = 1 endif endloop scalar var_beta = grad_beta * V * grad_beta' scalar se_beta = (var_beta > 0) ? sqrt(var_beta) : NA # alpha1 matrix grad_alpha1 = zeros(1, n_params_actual) if idx_alpha1 <= n_params_actual grad_alpha1[1, idx_alpha1] = 1 endif scalar var_alpha1 = grad_alpha1 * V * grad_alpha1' width_alpha1[pos] = (var_alpha1 > 0) ? 1.96 * sqrt(var_alpha1) : NA # alpha2 if arch_order >= 2 matrix grad_alpha2 = zeros(1, n_params_actual) if idx_alpha2 <= n_params_actual grad_alpha2[1, idx_alpha2] = 1 endif scalar var_alpha2 = grad_alpha2 * V * grad_alpha2' width_alpha2[pos] = (var_alpha2 > 0) ? 1.96 * sqrt(var_alpha2) : NA endif # beta1 matrix grad_beta1 = zeros(1, n_params_actual) if idx_beta1 <= n_params_actual grad_beta1[1, idx_beta1] = 1 endif scalar var_beta1 = grad_beta1 * V * grad_beta1' width_beta1[pos] = (var_beta1 > 0) ? 1.96 * sqrt(var_beta1) : NA # beta2 if garch_order >= 2 matrix grad_beta2 = zeros(1, n_params_actual) if idx_beta2 <= n_params_actual grad_beta2[1, idx_beta2] = 1 endif scalar var_beta2 = grad_beta2 * V * grad_beta2' width_beta2[pos] = (var_beta2 > 0) ? 1.96 * sqrt(var_beta2) : NA endif # STORE results roll_alpha0[pos] = o roll_alpha[pos] = sum_alpha roll_beta[pos] = sum_beta roll_sigma2[pos] = sigma2 roll_sigma[pos] = sigma roll_aic[pos] = $aic roll_bic[pos] = $bic roll_logL[pos] = $lnl if !missing(se_alpha0) width_alpha0[pos] = 1.96 * se_alpha0 endif if !missing(se_ab) width_ab[pos] = 1.96 * se_ab endif if !missing(se_alpha) width_alpha[pos] = 1.96 * se_alpha endif if !missing(se_beta) width_beta[pos] = 1.96 * se_beta endif endif endloop smpl full set verbose on # INTERPOLATION series roll_alpha0_plot = interpol(roll_alpha0) series roll_sigma2_plot = interpol(roll_sigma2) series roll_sigma_plot = interpol(roll_sigma) series roll_alpha_plot = interpol(roll_alpha) series roll_beta_plot = interpol(roll_beta) series ab_plot = roll_alpha_plot + roll_beta_plot series roll_aic_plot = interpol(roll_aic) series roll_bic_plot = interpol(roll_bic) series roll_logL_plot = interpol(roll_logL) series width_alpha0_plot = interpol(width_alpha0) series width_ab_plot = interpol(width_ab) series width_alpha_plot = interpol(width_alpha) series width_beta_plot = interpol(width_beta) series roll_alpha1_plot = interpol(roll_alpha1) series roll_alpha2_plot = interpol(roll_alpha2) series roll_beta1_plot = interpol(roll_beta1) series roll_beta2_plot = interpol(roll_beta2) series width_alpha1_plot = interpol(width_alpha1) series width_alpha2_plot = interpol(width_alpha2) series width_beta1_plot = interpol(width_beta1) series width_beta2_plot = interpol(width_beta2) # Clean extreme ci widths series width_alpha0_clean = (width_alpha0_plot > ci_trim) ? NA : width_alpha0_plot series width_alpha_clean = (width_alpha_plot > ci_trim) ? NA : width_alpha_plot series width_beta_clean = (width_beta_plot > ci_trim) ? NA : width_beta_plot series width_ab_clean = (width_ab_plot > ci_trim) ? NA : width_ab_plot series width_alpha1_clean = (width_alpha1_plot > ci_trim) ? NA : width_alpha1_plot series width_alpha2_clean = (width_alpha2_plot > ci_trim) ? NA : width_alpha2_plot series width_beta1_clean = (width_beta1_plot > ci_trim) ? NA : width_beta1_plot series width_beta2_clean = (width_beta2_plot > ci_trim) ? NA : width_beta2_plot # CI BUNDLES bundle CI_alpha0 = _(center="roll_alpha0_plot", width="width_alpha0_clean", style="fill") CI_alpha0.color = 0xCDCDCD bundle CI_alpha = _(center="roll_alpha_plot", width="width_alpha_clean", style="fill") CI_alpha.color = 0xCDCDCD bundle CI_beta = _(center="roll_beta_plot", width="width_beta_clean", style="fill") CI_beta.color = 0xCDCDCD bundle CI_ab = _(center="ab_plot", width="width_ab_clean", style="fill") CI_ab.color = 0xCDCDCD bundle CI_alpha1 = _(center="roll_alpha1_plot", width="width_alpha1_clean", style="fill") CI_alpha1.color = 0xCDCDCD bundle CI_alpha2 = _(center="roll_alpha2_plot", width="width_alpha2_clean", style="fill") CI_alpha2.color = 0xCDCDCD bundle CI_beta1 = _(center="roll_beta1_plot", width="width_beta1_clean", style="fill") CI_beta1.color = 0xCDCDCD bundle CI_beta2 = _(center="roll_beta2_plot", width="width_beta2_clean", style="fill") CI_beta2.color = 0xCDCDCD # PLOTS gnuplot roll_alpha0_plot --time-series --with-lines --band=CI_alpha0 {set title "Rolling Variance Intercept (a0)"; set xlabel "Date"; set ylabel "alpha(0)"; set grid; set yrange [0:*];} string ab_title = "" loop j = 1 .. arch_order --quiet if j == 1 ab_title = ab_title ~ "alpha(1)" else ab_title = ab_title ~ " + alpha(" ~ sprintf("%d", j) ~ ")" endif endloop loop j = 1 .. garch_order --quiet ab_title = ab_title ~ " + beta(" ~ sprintf("%d", j) ~ ")" endloop gnuplot ab_plot --time-series --with-lines --band=CI_ab {set title "@ab_title"; set xlabel "Date"; set ylabel "@ab_title"; set grid; } if arch_order >= 1 gnuplot roll_alpha1_plot --time-series --with-lines --band=CI_alpha1 {set title "Rolling a1 (95% CI)"; set xlabel "Date"; set ylabel "a1"; set grid;} endif if arch_order >= 2 gnuplot roll_alpha2_plot --time-series --with-lines --band=CI_alpha2 {set title "Rolling a2 (95% CI)"; set xlabel "Date"; set ylabel "a2"; set grid;} endif if garch_order >= 1 gnuplot roll_beta1_plot --time-series --with-lines --band=CI_beta1 {set title "Rolling b1 (95% CI)"; set xlabel "Date"; set ylabel "b1"; set grid;} endif if garch_order >= 2 gnuplot roll_beta2_plot --time-series --with-lines --band=CI_beta2 {set title "Rolling b2 (95% CI)"; set xlabel "Date"; set ylabel "b2"; set grid;} endif if arch_order >=2 gnuplot roll_alpha_plot --time-series --with-lines --band=CI_alpha {set title "Rolling Sum of ARCH Coefficients (95% CI)"; set xlabel "Date"; set ylabel "sum(alpha)"; set grid;} endif if garch_order >=2 gnuplot roll_beta_plot --time-series --with-lines --band=CI_beta {set title "Rolling Sum of GARCH Coefficients (95% CI)"; set xlabel "Date"; set ylabel "sum(beta)"; set grid;} endif gnuplot roll_aic_plot --time-series --with-lines {set title "Rolling AIC"; set xlabel "Date"; set ylabel "AIC"; set grid;} gnuplot roll_bic_plot --time-series --with-lines {set title "Rolling BIC"; set xlabel "Date"; set ylabel "BIC"; set grid;} gnuplot roll_logL_plot --time-series --with-lines {set title "Rolling Log-Likelihood"; set xlabel "Date"; set ylabel "log L"; set grid;} # Pack results R["window"] = WINDOW R["arch_order"] = arch_order R["garch_order"] = garch_order R["roll_alpha0"] = roll_alpha0 R["roll_alpha"] = roll_alpha R["roll_beta"] = roll_beta R["roll_sigma2"] = roll_sigma2 R["roll_sigma"] = roll_sigma R["roll_alpha_plot"] = roll_alpha_plot R["roll_beta_plot"] = roll_beta_plot R["ab_plot"] = ab_plot R["roll_alpha0_plot"] = roll_alpha0_plot R["width_ab_plot"] = width_ab_plot R["width_alpha_plot"] = width_alpha_plot R["width_beta_plot"] = width_beta_plot R["CI_alpha"] = CI_alpha R["CI_beta"] = CI_beta R["width_alpha_clean"] = width_alpha_clean R["width_beta_clean"] = width_beta_clean R["width_alpha0_plot"] = width_alpha0_plot R["width_alpha0_clean"] = width_alpha0_clean R["CI_alpha0"] = CI_alpha0 R["roll_aic"] = roll_aic R["roll_bic"] = roll_bic R["roll_logL"] = roll_logL R["roll_aic_plot"] = roll_aic_plot R["roll_bic_plot"] = roll_bic_plot R["roll_logL_plot"]= roll_logL_plot R["CI_ab"] = CI_ab if arch_order >= 1 R["roll_alpha1"] = roll_alpha1 R["roll_alpha1_plot"] = roll_alpha1_plot R["width_alpha1_plot"]= width_alpha1_plot R["width_alpha1_clean"]= width_alpha1_clean R["CI_alpha1"] = CI_alpha1 endif if arch_order >= 2 R["roll_alpha2"] = roll_alpha2 R["roll_alpha2_plot"] = roll_alpha2_plot R["width_alpha2_plot"]= width_alpha2_plot R["width_alpha2_clean"]= width_alpha2_clean R["CI_alpha2"] = CI_alpha2 endif if garch_order >= 1 R["roll_beta1"] = roll_beta1 R["roll_beta1_plot"] = roll_beta1_plot R["width_beta1_plot"]= width_beta1_plot R["width_beta1_clean"]= width_beta1_clean R["CI_beta1"] = CI_beta1 endif if garch_order >= 2 R["roll_beta2"] = roll_beta2 R["roll_beta2_plot"] = roll_beta2_plot R["width_beta2_plot"]= width_beta2_plot R["width_beta2_clean"]= width_beta2_clean R["CI_beta2"] = CI_beta2 endif return R include rolling_garch_p_q.gfn # EXAMPLE 1: GARCH(1,1) no mean regressors, no constant bundle R1 = rolling_garch_(y, 500, 1, 1, null, 0) # EXAMPLE 2: GARCH(1,1) with constant bundle R2 = rolling_garch_(y, 500, 1, 1, null, 1) # EXAMPLE 3: GARCH(1,1) with constant and mean regressors list X = x1 x2 bundle R3 = rolling_garch_(y, 500, 1, 1, X, 1) # EXAMPLE 4: GARCH(2,1) with constant bundle R4 = rolling_garch_(y, 500, 2, 1, null, 1) # EXAMPLE 5: GARCH(1,2) with constant bundle R5 = rolling_garch_(y, 500, 1, 2, null, 1) # EXAMPLE 6: GARCH(2,2) with constant (note: may be unstable) bundle R6 = rolling_garch_(y, 500, 2, 2, null, 1)