Rainfall Time-Series and GARCH Analysis
Lawrence V. Fulton
August 08, 2026
1 Rainfall Analysis
1.1 Load Packages and Data
library(fpp3)
library(kableExtra)
library(psych)
library(ggplot2)
library(fGarch)
set.seed(1234)
DATA_PATH <- paste0(
"C:/Users/lfult/OneDrive - bc.edu/",
"publications/living sustainably/rain2.csv"
)
mydata <- read.csv(DATA_PATH)
str(mydata)## 'data.frame': 856 obs. of 2 variables:
## $ YearMonth: chr "1946 Sep" "1946 Oct" "1946 Nov" "1946 Dec" ...
## $ Amt : num 15.78 1.31 1.86 2.43 2.14 ...head(mydata)## YearMonth Amt
## 1 1946 Sep 15.78
## 2 1946 Oct 1.31
## 3 1946 Nov 1.86
## 4 1946 Dec 2.43
## 5 1947 Jan 2.14
## 6 1947 Feb 0.29summary(mydata$Amt)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.7375 1.8750 2.5668 3.4000 18.07002 Descriptive Statistics
2.1 Monthly Distribution
myt <- ts(
mydata$Amt,
frequency = 12,
start = c(1946, 9)
)
mymonth <- factor(
cycle(myt),
levels = 1:12,
labels = month.abb
)
boxplot(
mydata$Amt ~ mymonth,
horizontal = FALSE,
notch = TRUE,
col = rainbow(12),
xlab = "Month",
ylab = "Rainfall (inches)",
main = "Monthly Distribution of Rainfall"
)
2.2 Descriptive Statistics Table
desc_table <- psych::describe(mydata$Amt)
desc_table %>%
kbl(
digits = 3,
caption = "Rainfall Descriptive Statistics"
) %>%
kable_classic(html_font = "Cambria")| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| X1 | 1 | 856 | 2.567 | 2.59 | 1.875 | 2.117 | 1.824 | 0 | 18.07 | 18.07 | 1.958 | 5.101 | 0.089 |
3 Build Time Series
3.1 Construct the Tsibble
mydata2 <- mydata[, -3] %>%
as_tibble() %>%
mutate(
Month = yearmonth(YearMonth)
) %>%
as_tsibble(index = Month)
print(mydata2)## # A tsibble: 856 x 3 [1M]
## YearMonth Amt Month
## <chr> <dbl> <mth>
## 1 1946 Sep 15.8 1946 Sep
## 2 1946 Oct 1.31 1946 Oct
## 3 1946 Nov 1.86 1946 Nov
## 4 1946 Dec 2.43 1946 Dec
## 5 1947 Jan 2.14 1947 Jan
## 6 1947 Feb 0.29 1947 Feb
## 7 1947 Mar 1.46 1947 Mar
## 8 1947 Apr 0.3 1947 Apr
## 9 1947 May 3.32 1947 May
## 10 1947 Jun 0.31 1947 Jun
## # ℹ 846 more rows3.2 Baseline Time-Series Plot
autoplot(mydata2, Amt) +
labs(
title = "Monthly Rainfall",
x = "Year",
y = "Rainfall (inches)"
) +
theme_minimal()
3.3 Seasonal Plot
mydata2 %>%
gg_season(Amt, labels = "both") +
labs(
title = "Seasonal Rainfall Pattern",
y = "Rainfall (inches)",
x = "Month"
) +
theme_minimal()
3.4 Seasonal Subseries
mydata2 %>%
gg_subseries(Amt) +
labs(
title = "Monthly Rainfall Subseries",
y = "Rainfall (inches)"
) +
theme_minimal()
3.5 Autocorrelation Function
mydata2 %>%
ACF(Amt) %>%
autoplot() +
labs(title = "Rainfall Autocorrelation Function") +
theme_minimal()
3.6 Partial Autocorrelation Function
mydata2 %>%
PACF(Amt) %>%
autoplot() +
labs(title = "Rainfall Partial Autocorrelation Function") +
theme_minimal()
4 Decomposition
4.1 STL Decomposition
dcmp_stl <- mydata2 %>%
model(
STL = STL(Amt)
)4.2 Classical Decomposition
dcmp_classical <- mydata2 %>%
model(
Classical = classical_decomposition(Amt)
)4.3 STL Trend
components(dcmp_stl) %>%
as_tsibble() %>%
autoplot(Amt, color = "gray") +
geom_line(
aes(y = trend),
colour = "#0072B2",
linewidth = 1
) +
labs(
title = "STL Trend",
x = "Year",
y = "Rainfall (inches)"
) +
theme_minimal()
4.4 Full STL Decomposition
components(dcmp_stl) %>%
autoplot() +
labs(title = "STL Decomposition")
4.5 Classical Trend
components(dcmp_classical) %>%
as_tsibble() %>%
autoplot(Amt, color = "gray") +
geom_line(
aes(y = trend),
colour = "#0072B2",
linewidth = 1
) +
labs(
title = "Classical Decomposition Trend",
x = "Year",
y = "Rainfall (inches)"
) +
theme_minimal()
4.6 Full Classical Decomposition
components(dcmp_classical) %>%
autoplot() +
labs(title = "Classical Time-Series Decomposition")
5 Forecasting Models
5.1 Training and Test Sets
mytrain <- mydata2 %>%
filter_index("1946 Sep" ~ "2004 Dec")
mytest <- mydata2 %>%
filter_index("2005 Jan" ~ "2017 Dec")
cat("Training observations:", nrow(mytrain), "\n")## Training observations: 700cat("Testing observations:", nrow(mytest), "\n")## Testing observations: 156cat(
"Training period:",
as.character(min(mytrain$Month)),
"to",
as.character(max(mytrain$Month)),
"\n"
)## Training period: 1946 Sep to 2004 Deccat(
"Testing period:",
as.character(min(mytest$Month)),
"to",
as.character(max(mytest$Month)),
"\n"
)## Testing period: 2005 Jan to 2017 Dec5.2 Fit Baseline Forecasting Models
myfit <- mytrain %>%
model(
SNAIVE = SNAIVE(Amt),
ETS = ETS(Amt),
ARIMA = ARIMA(Amt)
)
report(myfit)## # A tibble: 3 × 11
## .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE ar_roots ma_roots
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
## 1 SNAIVE 12.0 NA NA NA NA NA NA NA <NULL> <NULL>
## 2 ETS 6.13 -2920. 5871. 5871. 5939. 6.01 5.87 1.73 <NULL> <NULL>
## 3 ARIMA 8.51 -1739. 3486. 3486. 3504. NA NA NA <cpl [2]> <cpl [12]>5.3 Forecast Holdout Period
myforecast <- myfit %>%
forecast(new_data = mytest)5.4 Plot Forecasts
myforecast %>%
autoplot(
mydata2 %>%
filter_index("2000 Jan" ~ "2017 Dec"),
level = NULL
) +
guides(
color = guide_legend(title = "Model")
) +
labs(
title = "Rainfall Forecast Comparison",
subtitle = "Holdout period: January 2005 to December 2017",
x = "Year",
y = "Rainfall (inches)"
) +
theme_minimal()
5.5 Forecast Accuracy
forecast_accuracy <- accuracy(
myforecast,
mytest
)
forecast_accuracy## # A tibble: 3 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA Test -2.42 3.58 3.22 -2252. 2261. NaN NaN 0.212
## 2 ETS Test -0.261 2.56 2.03 -1196. 1225. NaN NaN 0.191
## 3 SNAIVE Test -1.22 4.63 3.40 -2322. 2371. NaN NaN -0.0572forecast_accuracy %>%
select(.model, RMSE, MAE, MAPE, MASE) %>%
kbl(
digits = 3,
caption = "Out-of-Sample Forecast Accuracy"
) %>%
kable_classic(html_font = "Cambria")| .model | RMSE | MAE | MAPE | MASE |
|---|---|---|---|---|
| ARIMA | 3.580 | 3.225 | 2261.043 | NaN |
| ETS | 2.558 | 2.031 | 1224.936 | NaN |
| SNAIVE | 4.631 | 3.402 | 2370.635 | NaN |
6 Volatility Models
6.1 Convert Training and
Test Samples to ts
train <- ts(
mytrain$Amt,
start = c(1946, 9),
frequency = 12
)
test <- ts(
mytest$Amt,
start = c(2005, 1),
frequency = 12
)6.2 Fit ARMA-GARCH(1,1)
The ARMA component models the conditional mean, while the GARCH component models the conditional variance.
mygarch <- garchFit(
formula = ~ arma(1, 1) + garch(1, 1),
data = train,
cond.dist = "QMLE",
trace = FALSE
)
coef(mygarch)## mu ar1 ma1 omega alpha1 beta1
## 2.563305599 -0.007431047 0.083211308 2.896684955 0.019820577 0.525622880summary(mygarch)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(1, 1) + garch(1, 1), data = train, cond.dist = "QMLE",
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x000001d25166f510>
## [data = train]
##
## Conditional Distribution:
## QMLE
##
## Coefficient(s):
## mu ar1 ma1 omega alpha1 beta1
## 2.563306 -0.007431 0.083211 2.896685 0.019821 0.525623
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 2.563306 0.271919 9.427 < 2e-16 ***
## ar1 -0.007431 0.097667 -0.076 0.93935
## ma1 0.083211 0.088874 0.936 0.34913
## omega 2.896685 0.882293 3.283 0.00103 **
## alpha1 0.019821 0.030757 0.644 0.51929
## beta1 0.525623 0.133405 3.940 8.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -1640.195 normalized: -2.343136
##
## Description:
## Sat Aug 8 11:09:03 2026 by user: lfult
##
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 1431.7292319 0.0000000
## Shapiro-Wilk Test R W 0.8140743 0.0000000
## Ljung-Box Test R Q(10) 2.6934726 0.9877466
## Ljung-Box Test R Q(15) 14.8189127 0.4645383
## Ljung-Box Test R Q(20) 21.9244980 0.3446216
## Ljung-Box Test R^2 Q(10) 4.5980260 0.9163646
## Ljung-Box Test R^2 Q(15) 12.1786974 0.6654557
## Ljung-Box Test R^2 Q(20) 20.2123488 0.4447187
## LM Arch Test R TR^2 9.6146081 0.6497299
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## 4.703414 4.742424 4.703269 4.7184947 GARCH Diagnostics
7.1 Conditional Standard Deviation
plot(
mygarch,
which = 2
)
7.2 Series with Conditional Standard-Deviation Bands
plot(
mygarch,
which = 3
)
7.3 Diagnostic Panel
old_par <- par(mfrow = c(2, 2))
plot(mygarch, which = 2)
plot(mygarch, which = 9)
plot(mygarch, which = 10)
plot(mygarch, which = 11)
par(old_par)8 GARCH Forecasting
8.1 Forecast the Holdout Period
mypred <- predict(
mygarch,
n.ahead = length(test)
)
head(mypred)## meanForecast meanError standardDeviation
## 1 2.317694 2.645452 2.645452
## 2 2.546083 2.598868 2.591124
## 3 2.544386 2.568523 2.561006
## 4 2.544398 2.551819 2.544428
## 5 2.544398 2.542662 2.535340
## 6 2.544398 2.537653 2.530369tail(mypred)## meanForecast meanError standardDeviation
## 151 2.544398 2.53163 2.524391
## 152 2.544398 2.53163 2.524391
## 153 2.544398 2.53163 2.524391
## 154 2.544398 2.53163 2.524391
## 155 2.544398 2.53163 2.524391
## 156 2.544398 2.53163 2.5243918.2 Convert GARCH Forecasts to Time Series
garch_mean <- ts(
mypred$meanForecast,
start = start(test),
frequency = frequency(test)
)
garch_sd <- ts(
mypred$standardDeviation,
start = start(test),
frequency = frequency(test)
)
garch_se <- ts(
mypred$meanError,
start = start(test),
frequency = frequency(test)
)8.3 GARCH Forecast Accuracy
garch_error <- as.numeric(test) - as.numeric(garch_mean)
GARCH_ME <- mean(
garch_error,
na.rm = TRUE
)
GARCH_MAE <- mean(
abs(garch_error),
na.rm = TRUE
)
GARCH_MSE <- mean(
garch_error^2,
na.rm = TRUE
)
GARCH_RMSE <- sqrt(GARCH_MSE)
GARCH_MAPE <- mean(
abs(
garch_error /
ifelse(
as.numeric(test) == 0,
NA,
as.numeric(test)
)
),
na.rm = TRUE
) * 100
garch_accuracy <- tibble(
Model = "ARMA(1,1)-GARCH(1,1)",
ME = GARCH_ME,
MAE = GARCH_MAE,
MSE = GARCH_MSE,
RMSE = GARCH_RMSE,
MAPE = GARCH_MAPE
)
garch_accuracy## # A tibble: 1 × 6
## Model ME MAE MSE RMSE MAPE
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARMA(1,1)-GARCH(1,1) 0.0187 1.99 6.94 2.63 1123.garch_accuracy %>%
kbl(
digits = 3,
caption = "ARMA-GARCH Out-of-Sample Accuracy"
) %>%
kable_classic(html_font = "Cambria")| Model | ME | MAE | MSE | RMSE | MAPE |
|---|---|---|---|---|---|
| ARMA(1,1)-GARCH(1,1) | 0.019 | 1.992 | 6.939 | 2.634 | 1123.166 |
9 GARCH Volatility Plots
9.1 Forecasted Conditional Volatility
plot(
garch_sd,
type = "l",
lwd = 2,
xlab = "Year",
ylab = "Conditional Standard Deviation",
main = "Forecasted Rainfall Volatility: ARMA-GARCH(1,1)"
)
grid()
9.2 Observed Rainfall and Forecast Volatility
old_par <- par(
mfrow = c(2, 1),
mar = c(4, 4, 3, 1)
)
plot(
test,
type = "l",
lwd = 2,
xlab = "",
ylab = "Rainfall (inches)",
main = "Observed Rainfall: Holdout Sample"
)
grid()
plot(
garch_sd,
type = "l",
lwd = 2,
xlab = "Year",
ylab = "Conditional SD",
main = "Forecasted GARCH Volatility"
)
grid()
par(old_par)10 GARCH Prediction Intervals
10.1 Construct 95 Percent Prediction Intervals
garch_upper95 <- garch_mean + 1.96 * garch_se
garch_lower95 <- garch_mean - 1.96 * garch_se
garch_lower95 <- pmax(
garch_lower95,
0
)10.2 Plot GARCH Forecast and Prediction Intervals
plot(
test,
type = "l",
lwd = 2,
ylim = range(
c(
test,
garch_mean,
garch_lower95,
garch_upper95
),
na.rm = TRUE
),
xlab = "Year",
ylab = "Rainfall (inches)",
main = "ARMA-GARCH Rainfall Forecast"
)
lines(
garch_mean,
lwd = 2
)
lines(
garch_upper95,
lty = 2,
lwd = 1.5
)
lines(
garch_lower95,
lty = 2,
lwd = 1.5
)
legend(
"topright",
legend = c(
"Observed",
"GARCH Mean Forecast",
"95% Upper",
"95% Lower"
),
lty = c(1, 1, 2, 2),
lwd = c(2, 2, 1.5, 1.5),
bty = "n"
)
grid()
11 Direct
fGarch Forecast Plot
predict(
mygarch,
n.ahead = length(test),
plot = TRUE,
conf = 0.95,
nx = 60
)
## meanForecast meanError standardDeviation lowerInterval upperInterval
## 1 2.317694 2.645452 2.645452 -0.333972136 9.228321
## 2 2.546083 2.598868 2.591124 -0.058889843 9.335020
## 3 2.544386 2.568523 2.561006 -0.030170604 9.254053
## 4 2.544398 2.551819 2.544428 -0.013414864 9.210431
## 5 2.544398 2.542662 2.535340 -0.004236169 9.186509
## 6 2.544398 2.537653 2.530369 0.000784302 9.173425
## 7 2.544398 2.534917 2.527654 0.003526866 9.166278
## 8 2.544398 2.533423 2.526171 0.005024028 9.162376
## 9 2.544398 2.532608 2.525362 0.005841017 9.160247
## 10 2.544398 2.532163 2.524921 0.006286749 9.159085
## 11 2.544398 2.531921 2.524680 0.006529904 9.158451
## 12 2.544398 2.531788 2.524549 0.006662541 9.158106
## 13 2.544398 2.531716 2.524477 0.006734890 9.157917
## 14 2.544398 2.531677 2.524438 0.006774353 9.157814
## 15 2.544398 2.531655 2.524417 0.006795878 9.157758
## 16 2.544398 2.531644 2.524405 0.006807619 9.157728
## 17 2.544398 2.531637 2.524399 0.006814023 9.157711
## 18 2.544398 2.531634 2.524395 0.006817516 9.157702
## 19 2.544398 2.531632 2.524394 0.006819421 9.157697
## 20 2.544398 2.531631 2.524393 0.006820461 9.157694
## 21 2.544398 2.531630 2.524392 0.006821027 9.157693
## 22 2.544398 2.531630 2.524392 0.006821337 9.157692
## 23 2.544398 2.531630 2.524392 0.006821505 9.157691
## 24 2.544398 2.531630 2.524391 0.006821597 9.157691
## 25 2.544398 2.531630 2.524391 0.006821647 9.157691
## 26 2.544398 2.531630 2.524391 0.006821675 9.157691
## 27 2.544398 2.531630 2.524391 0.006821690 9.157691
## 28 2.544398 2.531630 2.524391 0.006821698 9.157691
## 29 2.544398 2.531630 2.524391 0.006821702 9.157691
## 30 2.544398 2.531630 2.524391 0.006821705 9.157691
## 31 2.544398 2.531630 2.524391 0.006821706 9.157691
## 32 2.544398 2.531630 2.524391 0.006821707 9.157691
## 33 2.544398 2.531630 2.524391 0.006821707 9.157691
## 34 2.544398 2.531630 2.524391 0.006821707 9.157691
## 35 2.544398 2.531630 2.524391 0.006821707 9.157691
## 36 2.544398 2.531630 2.524391 0.006821707 9.157691
## 37 2.544398 2.531630 2.524391 0.006821707 9.157691
## 38 2.544398 2.531630 2.524391 0.006821707 9.157691
## 39 2.544398 2.531630 2.524391 0.006821707 9.157691
## 40 2.544398 2.531630 2.524391 0.006821707 9.157691
## 41 2.544398 2.531630 2.524391 0.006821708 9.157691
## 42 2.544398 2.531630 2.524391 0.006821708 9.157691
## 43 2.544398 2.531630 2.524391 0.006821708 9.157691
## 44 2.544398 2.531630 2.524391 0.006821708 9.157691
## 45 2.544398 2.531630 2.524391 0.006821708 9.157691
## 46 2.544398 2.531630 2.524391 0.006821708 9.157691
## 47 2.544398 2.531630 2.524391 0.006821708 9.157691
## 48 2.544398 2.531630 2.524391 0.006821708 9.157691
## 49 2.544398 2.531630 2.524391 0.006821708 9.157691
## 50 2.544398 2.531630 2.524391 0.006821708 9.157691
## 51 2.544398 2.531630 2.524391 0.006821708 9.157691
## 52 2.544398 2.531630 2.524391 0.006821708 9.157691
## 53 2.544398 2.531630 2.524391 0.006821708 9.157691
## 54 2.544398 2.531630 2.524391 0.006821708 9.157691
## 55 2.544398 2.531630 2.524391 0.006821708 9.157691
## 56 2.544398 2.531630 2.524391 0.006821708 9.157691
## 57 2.544398 2.531630 2.524391 0.006821708 9.157691
## 58 2.544398 2.531630 2.524391 0.006821708 9.157691
## 59 2.544398 2.531630 2.524391 0.006821708 9.157691
## 60 2.544398 2.531630 2.524391 0.006821708 9.157691
## 61 2.544398 2.531630 2.524391 0.006821708 9.157691
## 62 2.544398 2.531630 2.524391 0.006821708 9.157691
## 63 2.544398 2.531630 2.524391 0.006821708 9.157691
## 64 2.544398 2.531630 2.524391 0.006821708 9.157691
## 65 2.544398 2.531630 2.524391 0.006821708 9.157691
## 66 2.544398 2.531630 2.524391 0.006821708 9.157691
## 67 2.544398 2.531630 2.524391 0.006821708 9.157691
## 68 2.544398 2.531630 2.524391 0.006821708 9.157691
## 69 2.544398 2.531630 2.524391 0.006821708 9.157691
## 70 2.544398 2.531630 2.524391 0.006821708 9.157691
## 71 2.544398 2.531630 2.524391 0.006821708 9.157691
## 72 2.544398 2.531630 2.524391 0.006821708 9.157691
## 73 2.544398 2.531630 2.524391 0.006821708 9.157691
## 74 2.544398 2.531630 2.524391 0.006821708 9.157691
## 75 2.544398 2.531630 2.524391 0.006821708 9.157691
## 76 2.544398 2.531630 2.524391 0.006821708 9.157691
## 77 2.544398 2.531630 2.524391 0.006821708 9.157691
## 78 2.544398 2.531630 2.524391 0.006821708 9.157691
## 79 2.544398 2.531630 2.524391 0.006821708 9.157691
## 80 2.544398 2.531630 2.524391 0.006821708 9.157691
## 81 2.544398 2.531630 2.524391 0.006821708 9.157691
## 82 2.544398 2.531630 2.524391 0.006821708 9.157691
## 83 2.544398 2.531630 2.524391 0.006821708 9.157691
## 84 2.544398 2.531630 2.524391 0.006821708 9.157691
## 85 2.544398 2.531630 2.524391 0.006821708 9.157691
## 86 2.544398 2.531630 2.524391 0.006821708 9.157691
## 87 2.544398 2.531630 2.524391 0.006821708 9.157691
## 88 2.544398 2.531630 2.524391 0.006821708 9.157691
## 89 2.544398 2.531630 2.524391 0.006821708 9.157691
## 90 2.544398 2.531630 2.524391 0.006821708 9.157691
## 91 2.544398 2.531630 2.524391 0.006821708 9.157691
## 92 2.544398 2.531630 2.524391 0.006821708 9.157691
## 93 2.544398 2.531630 2.524391 0.006821708 9.157691
## 94 2.544398 2.531630 2.524391 0.006821708 9.157691
## 95 2.544398 2.531630 2.524391 0.006821708 9.157691
## 96 2.544398 2.531630 2.524391 0.006821708 9.157691
## 97 2.544398 2.531630 2.524391 0.006821708 9.157691
## 98 2.544398 2.531630 2.524391 0.006821708 9.157691
## 99 2.544398 2.531630 2.524391 0.006821708 9.157691
## 100 2.544398 2.531630 2.524391 0.006821708 9.157691
## 101 2.544398 2.531630 2.524391 0.006821708 9.157691
## 102 2.544398 2.531630 2.524391 0.006821708 9.157691
## 103 2.544398 2.531630 2.524391 0.006821708 9.157691
## 104 2.544398 2.531630 2.524391 0.006821708 9.157691
## 105 2.544398 2.531630 2.524391 0.006821708 9.157691
## 106 2.544398 2.531630 2.524391 0.006821708 9.157691
## 107 2.544398 2.531630 2.524391 0.006821708 9.157691
## 108 2.544398 2.531630 2.524391 0.006821708 9.157691
## 109 2.544398 2.531630 2.524391 0.006821708 9.157691
## 110 2.544398 2.531630 2.524391 0.006821708 9.157691
## 111 2.544398 2.531630 2.524391 0.006821708 9.157691
## 112 2.544398 2.531630 2.524391 0.006821708 9.157691
## 113 2.544398 2.531630 2.524391 0.006821708 9.157691
## 114 2.544398 2.531630 2.524391 0.006821708 9.157691
## 115 2.544398 2.531630 2.524391 0.006821708 9.157691
## 116 2.544398 2.531630 2.524391 0.006821708 9.157691
## 117 2.544398 2.531630 2.524391 0.006821708 9.157691
## 118 2.544398 2.531630 2.524391 0.006821708 9.157691
## 119 2.544398 2.531630 2.524391 0.006821708 9.157691
## 120 2.544398 2.531630 2.524391 0.006821708 9.157691
## 121 2.544398 2.531630 2.524391 0.006821708 9.157691
## 122 2.544398 2.531630 2.524391 0.006821708 9.157691
## 123 2.544398 2.531630 2.524391 0.006821708 9.157691
## 124 2.544398 2.531630 2.524391 0.006821708 9.157691
## 125 2.544398 2.531630 2.524391 0.006821708 9.157691
## 126 2.544398 2.531630 2.524391 0.006821708 9.157691
## 127 2.544398 2.531630 2.524391 0.006821708 9.157691
## 128 2.544398 2.531630 2.524391 0.006821708 9.157691
## 129 2.544398 2.531630 2.524391 0.006821708 9.157691
## 130 2.544398 2.531630 2.524391 0.006821708 9.157691
## 131 2.544398 2.531630 2.524391 0.006821708 9.157691
## 132 2.544398 2.531630 2.524391 0.006821708 9.157691
## 133 2.544398 2.531630 2.524391 0.006821708 9.157691
## 134 2.544398 2.531630 2.524391 0.006821708 9.157691
## 135 2.544398 2.531630 2.524391 0.006821708 9.157691
## 136 2.544398 2.531630 2.524391 0.006821708 9.157691
## 137 2.544398 2.531630 2.524391 0.006821708 9.157691
## 138 2.544398 2.531630 2.524391 0.006821708 9.157691
## 139 2.544398 2.531630 2.524391 0.006821708 9.157691
## 140 2.544398 2.531630 2.524391 0.006821708 9.157691
## 141 2.544398 2.531630 2.524391 0.006821708 9.157691
## 142 2.544398 2.531630 2.524391 0.006821708 9.157691
## 143 2.544398 2.531630 2.524391 0.006821708 9.157691
## 144 2.544398 2.531630 2.524391 0.006821708 9.157691
## 145 2.544398 2.531630 2.524391 0.006821708 9.157691
## 146 2.544398 2.531630 2.524391 0.006821708 9.157691
## 147 2.544398 2.531630 2.524391 0.006821708 9.157691
## 148 2.544398 2.531630 2.524391 0.006821708 9.157691
## 149 2.544398 2.531630 2.524391 0.006821708 9.157691
## 150 2.544398 2.531630 2.524391 0.006821708 9.157691
## 151 2.544398 2.531630 2.524391 0.006821708 9.157691
## 152 2.544398 2.531630 2.524391 0.006821708 9.157691
## 153 2.544398 2.531630 2.524391 0.006821708 9.157691
## 154 2.544398 2.531630 2.524391 0.006821708 9.157691
## 155 2.544398 2.531630 2.524391 0.006821708 9.157691
## 156 2.544398 2.531630 2.524391 0.006821708 9.15769112 Model Comparison
12.1 Compare Baseline Models with GARCH
baseline_accuracy <- forecast_accuracy %>%
transmute(
Model = .model,
RMSE = RMSE,
MAE = MAE
)
garch_comparison <- garch_accuracy %>%
transmute(
Model = Model,
RMSE = RMSE,
MAE = MAE
)
model_comparison <- bind_rows(
baseline_accuracy,
garch_comparison
)
model_comparison## # A tibble: 4 × 3
## Model RMSE MAE
## <chr> <dbl> <dbl>
## 1 ARIMA 3.58 3.22
## 2 ETS 2.56 2.03
## 3 SNAIVE 4.63 3.40
## 4 ARMA(1,1)-GARCH(1,1) 2.63 1.99model_comparison %>%
arrange(RMSE) %>%
kbl(
digits = 3,
caption = "Holdout Forecast Model Comparison"
) %>%
kable_classic(html_font = "Cambria")| Model | RMSE | MAE |
|---|---|---|
| ETS | 2.558 | 2.031 |
| ARMA(1,1)-GARCH(1,1) | 2.634 | 1.992 |
| ARIMA | 3.580 | 3.225 |
| SNAIVE | 4.631 | 3.402 |
13 GARCH Persistence
For a standard GARCH(1,1) model,
\[ \sigma_t^2 = \omega + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2. \]
The quantity \(\alpha_1 + \beta_1\) measures volatility persistence.
garch_coef <- coef(mygarch)
if (
"alpha1" %in% names(garch_coef) &&
"beta1" %in% names(garch_coef)
) {
persistence <- (
garch_coef["alpha1"] +
garch_coef["beta1"]
)
cat(
"GARCH volatility persistence (alpha1 + beta1):",
round(persistence, 4),
"\n"
)
if (persistence < 1) {
cat(
"Interpretation: volatility is mean reverting.\n"
)
} else {
cat(
"Interpretation: volatility is extremely persistent ",
"or potentially nonstationary.\n"
)
}
}## GARCH volatility persistence (alpha1 + beta1): 0.5454
## Interpretation: volatility is mean reverting.14 Session Information
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] fGarch_4052.93 psych_2.5.3 kableExtra_1.4.0 fable_0.4.1
## [5] feasts_0.4.1 fabletools_0.5.0 tsibbledata_0.4.1 tsibble_1.1.6
## [9] ggplot2_4.0.3 lubridate_1.9.4 tidyr_1.3.2 dplyr_1.2.1
## [13] tibble_3.2.1 fpp3_1.0.1
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 anytime_0.3.11 xfun_0.52
## [4] bslib_0.9.0 lattice_0.22-6 vctrs_0.7.3
## [7] tools_4.5.0 Rdpack_2.6.4 generics_0.1.4
## [10] parallel_4.5.0 gbutils_0.5.1 pkgconfig_2.0.3
## [13] Matrix_1.7-3 RColorBrewer_1.1-3 S7_0.2.2
## [16] distributional_0.5.0 lifecycle_1.0.5 compiler_4.5.0
## [19] farver_2.1.2 stringr_1.5.1 mnormt_2.1.1
## [22] htmltools_0.5.8.1 sass_0.4.10 yaml_2.3.10
## [25] pillar_1.10.2 crayon_1.5.3 jquerylib_0.1.4
## [28] ellipsis_0.3.2 cachem_1.1.0 nlme_3.1-168
## [31] fBasics_4052.98 tidyselect_1.2.1 digest_0.6.37
## [34] stringi_1.8.7 purrr_1.0.4 labeling_0.4.3
## [37] fastmap_1.2.0 grid_4.5.0 cli_3.6.4
## [40] magrittr_2.0.3 utf8_1.2.5 withr_3.0.2
## [43] scales_1.4.0 rappdirs_0.3.3 timechange_0.3.0
## [46] rmarkdown_2.29 timeDate_4041.110 progressr_0.15.1
## [49] timeSeries_4052.112 urca_1.3-4 evaluate_1.0.3
## [52] knitr_1.50 rbibutils_2.3 viridisLite_0.4.2
## [55] rlang_1.2.0 spatial_7.3-18 Rcpp_1.0.14
## [58] glue_1.8.0 xml2_1.3.8 svglite_2.1.3
## [61] rstudioapi_0.17.1 jsonlite_2.0.0 R6_2.6.1
## [64] cvar_0.6 systemfonts_1.2.2