1 Package Installation

If any packages are missing, run this chunk once in the R Console before knitting.
It is intentionally set to eval=FALSE so knitting does not try to install packages.

install.packages(
  c("fpp3", "rugarch", "knitr", "rmarkdown"),
  repos = "https://cloud.r-project.org",
  dependencies = TRUE
)

2 Part 1: Hierarchical Time Series Forecasting

2.1 1.1 Exploring the Hierarchical Structure

For this section, I use the tourism dataset from fpp3. The data contain quarterly Australian domestic tourism trips. A natural geographic hierarchy is:

Australia to State to Region

The aggregate_key() function creates all levels of the hierarchy while preserving the relationship between State and Region.

data("tourism", package = "tsibbledata")

tourism_hts <- tourism |>
  aggregate_key(
    State / Region,
    Trips = sum(Trips)
  )

tourism_hts

## # A tsibble: 6,800 x 4 [1Q]
## # Key:       State, Region [85]
##    Quarter State        Region        Trips
##      <qtr> <chr*>       <chr*>        <dbl>
##  1 1998 Q1 <aggregated> <aggregated> 23182.
##  2 1998 Q2 <aggregated> <aggregated> 20323.
##  3 1998 Q3 <aggregated> <aggregated> 19827.
##  4 1998 Q4 <aggregated> <aggregated> 20830.
##  5 1999 Q1 <aggregated> <aggregated> 22087.
##  6 1999 Q2 <aggregated> <aggregated> 21458.
##  7 1999 Q3 <aggregated> <aggregated> 19914.
##  8 1999 Q4 <aggregated> <aggregated> 20028.
##  9 2000 Q1 <aggregated> <aggregated> 22339.
## 10 2000 Q2 <aggregated> <aggregated> 19941.
## # ℹ 6,790 more rows

The following visualization focuses on the national and state levels.

tourism_hts |>
  filter(is_aggregated(Region)) |>
  autoplot(Trips) +
  labs(
    title = "Australian Tourism: National and State-Level Trips",
    subtitle = "Quarterly overnight trips",
    x = "Quarter",
    y = "Trips ('000)"
  ) +
  facet_wrap(vars(State), scales = "free_y", ncol = 3) +
  theme(legend.position = "none")

To inspect a lower level of the hierarchy, the next graph shows regional tourism series for one state.

tourism |>
  filter(State == "New South Wales") |>
  autoplot(Trips) +
  labs(
    title = "Regional Tourism within New South Wales",
    x = "Quarter",
    y = "Trips ('000)",
    colour = "Region"
  ) +
  facet_wrap(vars(Region), scales = "free_y")

2.1.1 Interpretation

The hierarchy contains a national total at the highest level, state totals at the middle level, and individual regions at the bottom level. The plots show that different states and regions have distinct levels, trends, and seasonal patterns. This is important because a single national forecast may hide meaningful local behavior.

2.2 1.2 Train/Test Split

I reserve the final eight quarters (two years) as a test set.

tourism_cutoff <- max(tourism$Quarter) - 8

tourism_train <- tourism |>
  filter(Quarter <= tourism_cutoff) |>
  aggregate_key(
    State / Region,
    Trips = sum(Trips)
  )

tourism_test <- tourism |>
  filter(Quarter > tourism_cutoff) |>
  aggregate_key(
    State / Region,
    Trips = sum(Trips)
  )

range(tourism_train$Quarter)

## <yearquarter[2]>
## [1] "1998 Q1" "2015 Q4"
## # Year starts on: January

range(tourism_test$Quarter)

## <yearquarter[2]>
## [1] "2016 Q1" "2017 Q4"
## # Year starts on: January

2.3 1.3 Base Models and Reconciliation

An ETS model is fitted to each series. I then reconcile the forecasts using:

Bottom-up - forecasts the bottom-level regional series and aggregates upward.
Top-down - forecasts the total and distributes it downward.
Middle-out - begins from the middle of the hierarchy and combines upward aggregation with downward disaggregation.

tourism_fit <- tourism_train |>
  model(
    base = ETS(Trips)
  ) |>
  reconcile(
    bottom_up = bottom_up(base),
    top_down = top_down(base, method = "forecast_proportions"),
    middle_out = middle_out(base, split = 1)
  )

tourism_fit

## # A mable: 85 x 6
## # Key:     State, Region [85]
##    State          Region             base bottom_up    top_down     middle_out  
##    <chr*>         <chr*>          <model> <model>      <model>      <model>     
##  1 ACT          … Canberra … <ETS(A,N,N)> <ETS(A,N,N)> <ETS(A,N,N)> <ETS(A,N,N)>
##  2 ACT          … <aggregat… <ETS(A,N,N)> <ETS(A,N,N)> <ETS(A,N,N)> <ETS(A,N,N)>
##  3 New South Wal… Blue Moun… <ETS(M,N,M)> <ETS(M,N,M)> <ETS(M,N,M)> <ETS(M,N,M)>
##  4 New South Wal… Capital C… <ETS(M,N,A)> <ETS(M,N,A)> <ETS(M,N,A)> <ETS(M,N,A)>
##  5 New South Wal… Central C… <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)>
##  6 New South Wal… Central N… <ETS(M,N,M)> <ETS(M,N,M)> <ETS(M,N,M)> <ETS(M,N,M)>
##  7 New South Wal… Hunter   … <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)>
##  8 New South Wal… New Engla… <ETS(M,N,N)> <ETS(M,N,N)> <ETS(M,N,N)> <ETS(M,N,N)>
##  9 New South Wal… North Coa… <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)> <ETS(A,N,A)>
## 10 New South Wal… Outback N… <ETS(M,N,A)> <ETS(M,N,A)> <ETS(M,N,A)> <ETS(M,N,A)>
## # ℹ 75 more rows

2.4 1.4 Forecasts

tourism_fc <- tourism_fit |>
  forecast(h = 8)

tourism_fc

## # A fable: 2,720 x 6 [1Q]
## # Key:     State, Region, .model [340]
##    State  Region   .model    Quarter
##    <chr*> <chr*>   <chr>       <qtr>
##  1 ACT    Canberra base      2016 Q1
##  2 ACT    Canberra base      2016 Q2
##  3 ACT    Canberra base      2016 Q3
##  4 ACT    Canberra base      2016 Q4
##  5 ACT    Canberra base      2017 Q1
##  6 ACT    Canberra base      2017 Q2
##  7 ACT    Canberra base      2017 Q3
##  8 ACT    Canberra base      2017 Q4
##  9 ACT    Canberra bottom_up 2016 Q1
## 10 ACT    Canberra bottom_up 2016 Q2
## # ℹ 2,710 more rows
## # ℹ 2 more variables: Trips <dist>, .mean <dbl>

The national-level forecasts from all methods can be compared visually.

tourism_fc |>
  filter(is_aggregated(State), is_aggregated(Region)) |>
  autoplot(
    tourism_hts |>
      filter(is_aggregated(State), is_aggregated(Region))
  ) +
  labs(
    title = "National Tourism Forecasts by Reconciliation Method",
    x = "Quarter",
    y = "Trips ('000)",
    colour = "Method"
  )

The following plot compares state-level forecasts.

tourism_fc |>
  filter(!is_aggregated(State), is_aggregated(Region)) |>
  autoplot(
    tourism_hts |>
      filter(!is_aggregated(State), is_aggregated(Region))
  ) +
  labs(
    title = "State-Level Tourism Forecasts",
    x = "Quarter",
    y = "Trips ('000)",
    colour = "Method"
  ) +
  facet_wrap(vars(State), scales = "free_y")

2.5 1.5 Accuracy: RMSE and MAPE

Forecast accuracy is evaluated against the test data.

tourism_accuracy <- tourism_fc |>
  accuracy(
    tourism_test,
    measures = list(
      RMSE = RMSE,
      MAPE = MAPE
    )
  )

tourism_accuracy |>
  select(.model, State, Region, RMSE, MAPE) |>
  arrange(.model, RMSE) |>
  head(20) |>
  kable(
    digits = 2,
    caption = "Sample of Hierarchical Forecast Accuracy Results"
  )

Sample of Hierarchical Forecast Accuracy Results
.model	State	Region	RMSE	MAPE
base	Tasmania	Wilderness West	7.69	18.83
base	South Australia	Clare Valley	7.85	13.38
base	South Australia	Eyre Peninsula	10.36	9.56
base	South Australia	Barossa	10.93	14.52
base	Northern Territory	MacDonnell	11.07	58.99
base	South Australia	Kangaroo Island	11.15	29.91
base	Northern Territory	Barkly	12.30	43.41
base	Northern Territory	Kakadu Arnhem	13.93	18.27
base	Northern Territory	Lasseter	14.31	34.00
base	Victoria	Murray East	15.87	28.61
base	South Australia	Murraylands	16.34	18.45
base	South Australia	Adelaide Hills	18.85	40.89
base	Northern Territory	Alice Springs	19.45	21.10
base	Victoria	Wimmera	20.81	35.59
base	New South Wales	Outback NSW	21.12	14.06
base	South Australia	Limestone Coast	22.25	14.15
base	Victoria	Macedon	22.53	20.90
base	Northern Territory	Darwin	23.09	10.50
base	Victoria	Western Grampians	23.17	24.16
base	South Australia	Riverland	23.64	14.48

To make the comparison easier, I classify each series into a hierarchy level.

tourism_accuracy_by_level <- tourism_accuracy |>
  mutate(
    Level = case_when(
      is_aggregated(State) & is_aggregated(Region) ~ "National",
      !is_aggregated(State) & is_aggregated(Region) ~ "State",
      !is_aggregated(State) & !is_aggregated(Region) ~ "Region",
      TRUE ~ "Other"
    )
  ) |>
  group_by(.model, Level) |>
  summarise(
    RMSE = mean(RMSE, na.rm = TRUE),
    MAPE = mean(MAPE, na.rm = TRUE),
    .groups = "drop"
  )

tourism_accuracy_by_level |>
  arrange(Level, RMSE) |>
  kable(
    digits = 2,
    caption = "Average RMSE and MAPE by Hierarchical Level"
  )

Average RMSE and MAPE by Hierarchical Level
.model	Level	RMSE	MAPE
base	National	1720.72	5.22
top_down	National	1720.72	5.22
middle_out	National	2026.98	6.48
bottom_up	National	2514.19	8.73
top_down	Region	44.87	16.66
middle_out	Region	46.97	16.91
base	Region	52.64	17.54
bottom_up	Region	52.64	17.54
top_down	State	278.88	9.06
middle_out	State	306.85	9.83
base	State	306.85	9.83
bottom_up	State	388.67	11.34

tourism_accuracy_by_level |>
  ggplot(aes(x = .model, y = RMSE, fill = .model)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(vars(Level), scales = "free_y") +
  labs(
    title = "RMSE by Reconciliation Method and Hierarchical Level",
    x = "Forecast Method",
    y = "Average RMSE"
  )

3 Part 2: Grouped Time Series Forecasting

3.1 2.1 Exploring the Grouped Structure

For grouped forecasting, I use aus_retail. Australian retail turnover can be grouped across two dimensions:

State
Industry

Unlike a strict hierarchy, State and Industry are crossed dimensions. For example, an industry can appear in many states, and each state contains many industries.

data("aus_retail", package = "tsibbledata")

aus_retail

## # A tsibble: 64,532 x 5 [1M]
## # Key:       State, Industry [152]
##    State                        Industry           `Series ID`    Month Turnover
##    <chr>                        <chr>              <chr>          <mth>    <dbl>
##  1 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Apr      4.4
##  2 Australian Capital Territory Cafes, restaurant… A3349849A   1982 May      3.4
##  3 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Jun      3.6
##  4 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Jul      4  
##  5 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Aug      3.6
##  6 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Sep      4.2
##  7 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Oct      4.8
##  8 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Nov      5.4
##  9 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Dec      6.9
## 10 Australian Capital Territory Cafes, restaurant… A3349849A   1983 Jan      3.8
## # ℹ 64,522 more rows

The following plot compares total retail turnover by state.

aus_retail |>
  group_by(State) |>
  index_by(Month) |>
  summarise(
    Turnover = sum(Turnover, na.rm = TRUE)
  ) |>
  autoplot(Turnover) +
  labs(
    title = "Australian Retail Turnover by State",
    x = "Month",
    y = "Turnover"
  )

The next plot compares total turnover by industry.

aus_retail |>
  group_by(Industry) |>
  index_by(Month) |>
  summarise(
    Turnover = sum(Turnover, na.rm = TRUE)
  ) |>
  autoplot(Turnover) +
  labs(
    title = "Australian Retail Turnover by Industry",
    x = "Month",
    y = "Turnover"
  ) +
  facet_wrap(vars(Industry), scales = "free_y")

3.2 2.2 Construct the Grouped Time Series

The expression State * Industry creates aggregate series for State, Industry, State-by-Industry combinations, and the overall total.

retail_grouped <- aus_retail |>
  aggregate_key(
    State * Industry,
    Turnover = sum(Turnover)
  )

retail_grouped

## # A tsibble: 77,249 x 4 [1M]
## # Key:       State, Industry [181]
##       Month State        Industry     Turnover
##       <mth> <chr*>       <chr*>          <dbl>
##  1 1982 Apr <aggregated> <aggregated>    6225.
##  2 1982 May <aggregated> <aggregated>    6382.
##  3 1982 Jun <aggregated> <aggregated>    6162.
##  4 1982 Jul <aggregated> <aggregated>    6399.
##  5 1982 Aug <aggregated> <aggregated>    6163.
##  6 1982 Sep <aggregated> <aggregated>    6331.
##  7 1982 Oct <aggregated> <aggregated>    6535.
##  8 1982 Nov <aggregated> <aggregated>    7022.
##  9 1982 Dec <aggregated> <aggregated>    9322.
## 10 1983 Jan <aggregated> <aggregated>    6286.
## # ℹ 77,239 more rows

3.3 2.3 Train/Test Split

The final 24 months are reserved for testing.

retail_cutoff <- max(aus_retail$Month) - 24

retail_train <- aus_retail |>
  filter(Month <= retail_cutoff) |>
  aggregate_key(
    State * Industry,
    Turnover = sum(Turnover)
  )

retail_test <- aus_retail |>
  filter(Month > retail_cutoff) |>
  aggregate_key(
    State * Industry,
    Turnover = sum(Turnover)
  )

range(retail_train$Month)

## <yearmonth[2]>
## [1] "1982 Apr" "2016 Dec"

range(retail_test$Month)

## <yearmonth[2]>
## [1] "2017 Jan" "2018 Dec"

3.4 2.4 Fit ETS Models

I fit an ETS model separately to every series in the grouped structure.

retail_fit <- retail_train |>
  model(
    ets = ETS(Turnover)
  ) |>
  reconcile(
    bottom_up = bottom_up(ets),
    ols = min_trace(ets, method = "ols")
  )

retail_fit

## # A mable: 181 x 5
## # Key:     State, Industry [181]
##    State                       Industry            ets bottom_up    ols         
##    <chr*>                      <chr*>          <model> <model>      <model>     
##  1 Australian Capital Territo… Cafes, re… <ETS(M,A,A)> <ETS(M,A,A)> <ETS(M,A,A)>
##  2 Australian Capital Territo… Cafes, re… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  3 Australian Capital Territo… Clothing … <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  4 Australian Capital Territo… Clothing,… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  5 Australian Capital Territo… Departmen… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  6 Australian Capital Territo… Electrica… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  7 Australian Capital Territo… Food reta… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  8 Australian Capital Territo… Footwear … <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
##  9 Australian Capital Territo… Furniture… <ETS(M,A,M)> <ETS(M,A,M)> <ETS(M,A,M)>
## 10 Australian Capital Territo… Hardware,… <ETS(M,N,M)> <ETS(M,N,M)> <ETS(M,N,M)>
## # ℹ 171 more rows

Here, the unreconciled ets forecasts act as the flat / independently modeled forecasts, while mint produces forecasts that are coherent across the State and Industry grouping structure.

3.5 2.5 Forecast Grouped Series

retail_fc <- retail_fit |>
  forecast(h = 24)

retail_fc

## # A fable: 13,032 x 6 [1M]
## # Key:     State, Industry, .model [543]
##    State                        Industry                         .model    Month
##    <chr*>                       <chr*>                           <chr>     <mth>
##  1 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Jan
##  2 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Feb
##  3 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Mar
##  4 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Apr
##  5 Australian Capital Territory Cafes, restaurants and catering… ets    2017 May
##  6 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Jun
##  7 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Jul
##  8 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Aug
##  9 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Sep
## 10 Australian Capital Territory Cafes, restaurants and catering… ets    2017 Oct
## # ℹ 13,022 more rows
## # ℹ 2 more variables: Turnover <dist>, .mean <dbl>

The overall Australian retail forecast can be viewed below.

retail_fc |>
  filter(is_aggregated(State), is_aggregated(Industry)) |>
  autoplot(
    retail_grouped |>
      filter(is_aggregated(State), is_aggregated(Industry))
  ) +
  labs(
    title = "Total Australian Retail Turnover Forecast",
    x = "Month",
    y = "Turnover",
    colour = "Model"
  )

State-level grouped forecasts are shown below.

retail_fc |>
  filter(!is_aggregated(State), is_aggregated(Industry)) |>
  autoplot(
    retail_grouped |>
      filter(!is_aggregated(State), is_aggregated(Industry))
  ) +
  facet_wrap(vars(State), scales = "free_y") +
  labs(
    title = "Retail Forecasts by State",
    x = "Month",
    y = "Turnover"
  )

3.6 2.6 Grouped Forecast Accuracy

retail_accuracy <- retail_fc |>
  accuracy(
    retail_test,
    measures = list(
      RMSE = RMSE,
      MAPE = MAPE
    )
  )

retail_accuracy |>
  select(.model, State, Industry, RMSE, MAPE) |>
  arrange(.model, RMSE) |>
  head(20) |>
  kable(
    digits = 2,
    caption = "Sample of Grouped Forecast Accuracy Results"
  )

Sample of Grouped Forecast Accuracy Results
.model	State	Industry	RMSE	MAPE
bottom_up	Northern Territory	Newspaper and book retailing	0.34	23.03
bottom_up	Northern Territory	Footwear and other personal accessory retailing	0.48	5.72
bottom_up	Australian Capital Territory	Liquor retailing	0.66	4.11
bottom_up	Northern Territory	Clothing retailing	0.81	8.61
bottom_up	Northern Territory	Furniture, floor coverings, houseware and textile goods retailing	0.83	8.73
bottom_up	Tasmania	Footwear and other personal accessory retailing	0.87	9.37
bottom_up	South Australia	Newspaper and book retailing	0.90	7.03
bottom_up	Northern Territory	Other recreational goods retailing	1.06	16.75
bottom_up	Northern Territory	Pharmaceutical, cosmetic and toiletry goods retailing	1.25	12.35
bottom_up	Australian Capital Territory	Newspaper and book retailing	1.32	20.08
bottom_up	Northern Territory	Clothing, footwear and personal accessory retailing	1.32	8.27
bottom_up	Northern Territory	Hardware, building and garden supplies retailing	1.33	5.65
bottom_up	Northern Territory	Electrical and electronic goods retailing	1.48	9.67
bottom_up	Australian Capital Territory	Other recreational goods retailing	1.52	13.24
bottom_up	Australian Capital Territory	Department stores	1.55	3.39
bottom_up	Australian Capital Territory	Other specialised food retailing	1.57	10.61
bottom_up	Australian Capital Territory	Furniture, floor coverings, houseware and textile goods retailing	1.61	5.03
bottom_up	Australian Capital Territory	Footwear and other personal accessory retailing	1.69	8.28
bottom_up	Tasmania	Takeaway food services	1.85	5.78
bottom_up	Australian Capital Territory	Hardware, building and garden supplies retailing	1.96	5.48

Classify each series according to its grouping level.

retail_accuracy_by_level <- retail_accuracy |>
  mutate(
    Level = case_when(
      is_aggregated(State) & is_aggregated(Industry) ~ "Total",
      !is_aggregated(State) & is_aggregated(Industry) ~ "State",
      is_aggregated(State) & !is_aggregated(Industry) ~ "Industry",
      !is_aggregated(State) & !is_aggregated(Industry) ~ "State x Industry",
      TRUE ~ "Other"
    )
  ) |>
  group_by(.model, Level) |>
  summarise(
    RMSE = mean(RMSE, na.rm = TRUE),
    MAPE = mean(MAPE, na.rm = TRUE),
    .groups = "drop"
  )

retail_accuracy_by_level |>
  arrange(Level, RMSE) |>
  kable(
    digits = 2,
    caption = "Average Accuracy by Grouping Level"
  )

Average Accuracy by Grouping Level
.model	Level	RMSE	MAPE
ets	Industry	61.08	2.99
ols	Industry	66.44	3.55
bottom_up	Industry	73.44	5.13
ets	State	140.60	1.97
ols	State	155.25	3.10
bottom_up	State	171.59	2.47
bottom_up	State x Industry	14.70	6.92
ets	State x Industry	14.70	6.92
ols	State x Industry	16.39	13.09
bottom_up	Total	693.42	1.25
ols	Total	877.40	1.44
ets	Total	913.39	1.49

retail_accuracy_by_level |>
  ggplot(aes(x = .model, y = RMSE, fill = .model)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(vars(Level), scales = "free_y") +
  labs(
    title = "Flat ETS vs Reconciled MinT Forecasts",
    subtitle = "Lower RMSE indicates better forecast accuracy",
    x = "Model",
    y = "Average RMSE"
  )

3.7 2.7 Do Grouped Forecasts Outperform Flat Models?

The following table directly compares the unreconciled ETS forecasts with the MinT-reconciled forecasts.

retail_comparison <- retail_accuracy_by_level |>
  select(.model, Level, RMSE, MAPE) |>
  pivot_wider(
    names_from = .model,
    values_from = c(RMSE, MAPE)
  )

retail_comparison |>
  kable(
    digits = 2,
    caption = "Comparison of Flat and Reconciled Grouped Forecasts"
  )

Comparison of Flat and Reconciled Grouped Forecasts
Level	RMSE_bottom_up	RMSE_ets	RMSE_ols	MAPE_bottom_up	MAPE_ets	MAPE_ols
Industry	73.44	61.08	66.44	5.13	2.99	3.55
State	171.59	140.60	155.25	2.47	1.97	3.10
State x Industry	14.70	14.70	16.39	6.92	6.92	13.09
Total	693.42	913.39	877.40	1.25	1.49	1.44

4 Part 3: Volatility Modeling with ARCH/GARCH

4.1 3.1 Volatility Clustering

Volatility clustering refers to the tendency for periods of large changes to be followed by other large changes, while calmer periods tend to be followed by additional calm periods. Financial return series often display this behavior even when the mean return itself is difficult to predict.

To keep this code reproducible, I did a GARCH-type return series with changing conditional variance.

set.seed(123)

n <- 1200

omega_true <- 0.05
alpha_true <- 0.12
beta_true  <- 0.82

returns <- numeric(n)
variance <- numeric(n)

# Start at the theoretical unconditional variance.
variance[1] <- omega_true / (1 - alpha_true - beta_true)
returns[1] <- sqrt(variance[1]) * rnorm(1)

for (t in 2:n) {
  variance[t] <- omega_true +
    alpha_true * returns[t - 1]^2 +
    beta_true * variance[t - 1]

  returns[t] <- sqrt(variance[t]) * rnorm(1)
}

vol_data <- tibble(
  Time = 1:n,
  Return = returns,
  TrueVariance = variance
)

head(vol_data)

## # A tibble: 6 × 3
##    Time  Return TrueVariance
##   <int>   <dbl>        <dbl>
## 1     1 -0.512         0.833
## 2     2 -0.201         0.765
## 3     3  1.29          0.682
## 4     4  0.0634        0.808
## 5     5  0.109         0.713
## 6     6  1.37          0.636

4.2 3.2 Plot the Return Series

ggplot(vol_data, aes(x = Time, y = Return)) +
  geom_line() +
  labs(
    title = "Simulated Return Series with Volatility Clustering",
    x = "Time",
    y = "Return"
  )

The magnitude of returns changes over time. There are visible clusters of relatively high and low volatility.

4.3 3.3 Squared Returns

Squared returns provide a simple visual indicator of changing variance.

vol_data |>
  mutate(SquaredReturn = Return^2) |>
  ggplot(aes(x = Time, y = SquaredReturn)) +
  geom_line() +
  labs(
    title = "Squared Returns",
    subtitle = "Clusters of large squared returns suggest time-varying volatility",
    x = "Time",
    y = "Squared Return"
  )

4.4 3.4 Fit a GARCH(1,1) Model

A GARCH(1,1) model specifies the conditional variance as:

\[ \sigma_t^2 = \omega + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \]

where:

\(\omega\) is the constant variance component,
\(\alpha_1\) measures the immediate response to new shocks,
\(\beta_1\) measures persistence in volatility.

garch_spec <- ugarchspec(
  variance.model = list(
    model = "sGARCH",
    garchOrder = c(1, 1)
  ),
  mean.model = list(
    armaOrder = c(0, 0),
    include.mean = TRUE
  ),
  distribution.model = "norm"
)

garch_fit <- ugarchfit(
  spec = garch_spec,
  data = returns,
  solver = "hybrid"
)

garch_fit

## 
## *---------------------------------*
## *          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|)
## mu      0.013222    0.023211  0.56963 0.568928
## omega   0.073182    0.030446  2.40367 0.016232
## alpha1  0.098617    0.025395  3.88331 0.000103
## beta1   0.804396    0.057694 13.94238 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.013222    0.020448  0.64658 0.517901
## omega   0.073182    0.027586  2.65287 0.007981
## alpha1  0.098617    0.025744  3.83065 0.000128
## beta1   0.804396    0.053498 15.03600 0.000000
## 
## LogLikelihood : -1504.86 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       2.5148
## Bayes        2.5317
## Shibata      2.5147
## Hannan-Quinn 2.5212
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.839  0.1751
## Lag[2*(p+q)+(p+q)-1][2]     2.505  0.1913
## Lag[4*(p+q)+(p+q)-1][5]     3.151  0.3801
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.6809  0.4093
## Lag[2*(p+q)+(p+q)-1][5]    4.6427  0.1841
## Lag[4*(p+q)+(p+q)-1][9]    7.4821  0.1620
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     3.469 0.500 2.000 0.06252
## ARCH Lag[5]     6.014 1.440 1.667 0.05961
## ARCH Lag[7]     7.029 2.315 1.543 0.08562
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.4475
## Individual Statistics:              
## mu     0.03021
## omega  0.32200
## alpha1 0.30225
## beta1  0.31752
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.1224 0.2619    
## Negative Sign Bias  0.4781 0.6327    
## Positive Sign Bias  1.0756 0.2823    
## Joint Effect        1.5817 0.6635    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     9.267       0.9688
## 2    30    24.050       0.7265
## 3    40    31.733       0.7891
## 4    50    38.333       0.8641
## 
## 
## Elapsed time : 0.1433499

4.5 3.5 Parameter Estimates and Significance

garch_coef <- as.data.frame(garch_fit@fit$matcoef)

garch_coef |>
  kable(
    digits = 4,
    caption = "GARCH(1,1) Parameter Estimates"
  )

GARCH(1,1) Parameter Estimates
	Estimate	Std. Error	t value	Pr(>\|t\|)
mu	0.0132	0.0232	0.5696	0.5689
omega	0.0732	0.0304	2.4037	0.0162
alpha1	0.0986	0.0254	3.8833	0.0001
beta1	0.8044	0.0577	13.9424	0.0000

4.5.1 Parameter Interpretation

omega represents the long-run baseline contribution to variance.
alpha1 is the ARCH effect. A significant positive value means recent shocks affect current volatility.
beta1 is the GARCH effect. A significant positive value means volatility is persistent over time.
When alpha1 + beta1 is close to 1, volatility shocks tend to decay slowly.
For a covariance-stationary GARCH(1,1), we generally expect alpha1 + beta1 < 1.

coef_est <- coef(garch_fit)

alpha_hat <- unname(coef_est["alpha1"])
beta_hat  <- unname(coef_est["beta1"])
omega_hat <- unname(coef_est["omega"])

persistence <- alpha_hat + beta_hat

persistence

## [1] 0.9030125

4.6 3.6 Conditional Variance vs Unconditional Variance

The fitted conditional standard deviation is extracted with sigma(). Squaring it gives the estimated conditional variance.

conditional_sd <- as.numeric(sigma(garch_fit))
conditional_variance <- conditional_sd^2

unconditional_sample_variance <- var(returns)

garch_variance_df <- tibble(
  Time = 1:length(conditional_variance),
  ConditionalVariance = conditional_variance,
  UnconditionalVariance = unconditional_sample_variance
)

ggplot(garch_variance_df, aes(x = Time)) +
  geom_line(aes(y = ConditionalVariance)) +
  geom_hline(
    yintercept = unconditional_sample_variance,
    linetype = "dashed"
  ) +
  labs(
    title = "Conditional vs Unconditional Variance",
    subtitle = "Solid line = GARCH conditional variance; dashed line = constant sample variance",
    x = "Time",
    y = "Variance"
  )

The GARCH model allows variance to rise and fall through time, whereas the unconditional sample variance assumes one constant value for the entire series.

4.7 3.7 Estimated Long-Run Variance

For a stationary GARCH(1,1), the estimated long-run variance is:

\[ \frac{\omega}{1-\alpha_1-\beta_1} \]

if (persistence < 1) {
  estimated_long_run_variance <- omega_hat / (1 - persistence)
} else {
  estimated_long_run_variance <- NA_real_
}

variance_comparison <- tibble(
  Measure = c(
    "Sample unconditional variance",
    "Estimated GARCH long-run variance",
    "Mean fitted conditional variance"
  ),
  Variance = c(
    unconditional_sample_variance,
    estimated_long_run_variance,
    mean(conditional_variance, na.rm = TRUE)
  )
)

variance_comparison |>
  kable(
    digits = 4,
    caption = "Variance Comparison"
  )

Variance Comparison
Measure	Variance
Sample unconditional variance	0.7529
Estimated GARCH long-run variance	0.7546
Mean fitted conditional variance	0.7531

4.8 3.8 ARCH Effects Diagnostic

One way to motivate ARCH/GARCH modeling is to test whether squared returns are serially correlated. The Ljung-Box test below is applied to squared returns.

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

## 
##  Box-Ljung test
## 
## data:  returns^2
## X-squared = 103.14, df = 12, p-value < 2.2e-16

5 References

Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed.), Chapter 11: Forecasting hierarchical and grouped time series.
fabletools documentation: aggregate_key(), bottom_up(), top_down(), middle_out(), and min_trace().
Econometrics with R, Section 16.4, discussion of volatility clustering and ARCH/GARCH models.
Ghalanos, A. rugarch: Univariate GARCH models in R.

Hierarchical, Grouped, and Volatility Time Series Forecasting

Meghma Banerjee

August 08, 2026