Part 1 R Script:
library(fpp3)
tourism_data <- tourism |> mutate( State = recode(
State, Australian Capital Territory = “ACT”,
New South Wales = “NSW”,
Northern Territory = “NT”, Queensland =
“QLD”, South Australia = “SA”,
Tasmania = “TAS”, Victoria = “VIC”,
Western Australia = “WA” ) )
tourism_hts <- tourism_data |> aggregate_key( State / Region, Trips = sum(Trips) )
print(tourism_hts)
hierarchy_structure <- tourism_hts |> distinct(State, Region) |> mutate( Level = case_when( is_aggregated(State) & is_aggregated(Region) ~ “Total”, !is_aggregated(State) & is_aggregated(Region) ~ “State”, !is_aggregated(State) & !is_aggregated(Region) ~ “Region”, TRUE ~ “Other” ) ) |> count(Level, name = “Number_of_series”)
print(hierarchy_structure)
hierarchy_structure |> ggplot(aes(x = Level, y = Number_of_series)) + geom_col() + geom_text( aes(label = Number_of_series), vjust = -0.3 ) + labs( title = “Australian Tourism Hierarchy”, subtitle = “Australia → States → Regions”, x = “Hierarchy level”, y = “Number of time series” ) + theme_minimal()
total_state_data <- tourism_hts |> filter(is_aggregated(Region)) |> mutate( Series = if_else( is_aggregated(State), “Australia: Total”, as.character(State) ) )
total_state_data |> ggplot(aes(x = Quarter, y = Trips)) + geom_line() + facet_wrap( vars(Series), scales = “free_y”, ncol = 3 ) + labs( title = “Australian Overnight Tourism by Hierarchy Level”, subtitle = “National total and individual state series”, x = “Quarter”, y = “Trips, thousands” ) + theme_minimal()
victoria_regions <- tourism_hts |> filter( State == “VIC”, !is_aggregated(Region) ) |> mutate(Region = as.character(Region))
victoria_regions |> ggplot(aes(x = Quarter, y = Trips)) + geom_line() + facet_wrap( vars(Region), scales = “free_y”, ncol = 4 ) + labs( title = “Bottom-Level Tourism Series for Victoria”, subtitle = “Each panel represents one regional series”, x = “Quarter”, y = “Trips, thousands” ) + theme_minimal()
tourism_train <- tourism_hts |> filter(year(Quarter) <= 2015)
tourism_test <- tourism_hts |> filter(year(Quarter) >= 2016)
cat(“period:”) print(range(tourism_train$Quarter))
cat(“period:”) print(range(tourism_test$Quarter))
tourism_fit <- tourism_train |> model( base = ETS(Trips) ) |> reconcile( bu = bottom_up(base),
td = top_down( base, method = “forecast_proportions” ),
mo = middle_out( base, split = “State” ),
ols = min_trace( base, method = “ols” ),
mint = min_trace( base, method = “mint_shrink” ) )
print(tourism_fit)
tourism_fc <- tourism_fit |> forecast(h = “2 years”)
print(tourism_fc)
comparison_actual <- tourism_hts |> filter( year(Quarter) >= 2012, is_aggregated(Region), is_aggregated(State) | State == “VIC” ) |> mutate( Series = if_else( is_aggregated(State), “Australia: Total”, “Victoria: State total” ) )
comparison_forecasts <- tourism_fc |> filter( is_aggregated(Region), is_aggregated(State) | State == “VIC” ) |> mutate( Series = if_else( is_aggregated(State), “Australia: Total”, “Victoria: State total” ) )
ggplot() + geom_line( data = comparison_actual, aes(x = Quarter, y = Trips), linewidth = 0.8 ) + geom_line( data = comparison_forecasts, aes( x = Quarter, y = .mean, linetype = .model ), linewidth = 0.8 ) + facet_wrap( vars(Series), scales = “free_y”, ncol = 1 ) + labs( title = “Comparison of Hierarchical Forecasting Approaches”, subtitle = “Black line represents actual tourism; line types represent forecasts”, x = “Quarter”, y = “Trips, thousands”, linetype = “Method” ) + theme_minimal()
accuracy_individual <- tourism_fc |> accuracy( data = tourism_hts, measures = list( RMSE = RMSE, MAPE = MAPE ) ) |> mutate( Level = case_when( is_aggregated(State) & is_aggregated(Region) ~ “Total”, !is_aggregated(State) & is_aggregated(Region) ~ “State”, !is_aggregated(State) & !is_aggregated(Region) ~ “Region”, TRUE ~ “Other” ) )
print(accuracy_individual)
accuracy_by_level <- accuracy_individual |> group_by(Level, .model) |> summarise( RMSE = mean(RMSE, na.rm = TRUE), MAPE = mean(MAPE, na.rm = TRUE), Number_of_series = n(), .groups = “drop” ) |> mutate( Level = factor( Level, levels = c(“Total”, “State”, “Region”) ) ) |> arrange(Level, RMSE)
print(accuracy_by_level)
accuracy_table <- accuracy_by_level |> select(Level, .model, RMSE, MAPE) |> pivot_wider( names_from = .model, values_from = c(RMSE, MAPE) )
print(accuracy_table)
best_rmse <- accuracy_by_level |> group_by(Level) |> slice_min( order_by = RMSE, n = 1, with_ties = FALSE ) |> ungroup() |> select(Level, Best_RMSE_Method = .model, RMSE)
best_mape <- accuracy_by_level |> group_by(Level) |> slice_min( order_by = MAPE, n = 1, with_ties = FALSE ) |> ungroup() |> select(Level, Best_MAPE_Method = .model, MAPE)
best_methods <- best_rmse |> left_join(best_mape, by = “Level”)
print(best_methods)
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 Forecasting Method and Hierarchy Level”, x = “Forecasting method”, y = “Average RMSE” ) + theme_minimal()
accuracy_by_level |> ggplot( aes( x = .model, y = MAPE, fill = .model ) ) + geom_col(show.legend = FALSE) + facet_wrap( vars(Level), scales = “free_y” ) + labs( title = “MAPE by Forecasting Method and Hierarchy Level”, x = “Forecasting method”, y = “Average MAPE (%)” ) + theme_minimal()
coherence_test <- tourism_fc |> as_tibble() |> group_by(.model, Quarter) |> summarise( Total_forecast = .mean[ is_aggregated(State) & is_aggregated(Region) ][1],
Sum_of_regions = sum( .mean[ !is_aggregated(State) & !is_aggregated(Region) ], na.rm = TRUE ),
Difference = Total_forecast - Sum_of_regions, .groups = “drop” )
coherence_summary <- coherence_test |> group_by(.model) |> summarise( Maximum_absolute_difference = max(abs(Difference), na.rm = TRUE), .groups = “drop” )
print(coherence_summary)
write.csv( accuracy_by_level, “hierarchical_accuracy_by_level.csv”, row.names = FALSE )
write.csv( best_methods, “best_hierarchical_methods.csv”, row.names = FALSE )
Part 2 R Script:
library(fpp3)
broad_industries <- c( + “Food retailing”, + “Household goods retailing”, + “Clothing, footwear and personal accessory retailing”, + “Department stores”, + “Other retailing”, + “Cafes, restaurants and takeaway food services” + )
retail_data <- aus_retail |> + filter( + Industry %in% broad_industries, + year(Month) >= 2000 + ) |> + select( + State, + Industry, + Month, + Turnover + )
print(retail_data) # A tsibble: 10,032 x 4 [1M] # Key: State, Industry [44] State Industry Month Turnover
1 Australian Capital Territory Cafes, restaurants an… 2000 Jan 20.3 2 Australian Capital Territory Cafes, restaurants an… 2000 Feb 20.5 3 Australian Capital Territory Cafes, restaurants an… 2000 Mar 22.5 4 Australian Capital Territory Cafes, restaurants an… 2000 Apr 23.3 5 Australian Capital Territory Cafes, restaurants an… 2000 May 23.4 6 Australian Capital Territory Cafes, restaurants an… 2000 Jun 25.7 7 Australian Capital Territory Cafes, restaurants an… 2000 Jul 25.7 8 Australian Capital Territory Cafes, restaurants an… 2000 Aug 25.8 9 Australian Capital Territory Cafes, restaurants an… 2000 Sep 27.5 10 Australian Capital Territory Cafes, restaurants an… 2000 Oct 29.3 # ℹ 10,022 more rows # ℹ Use print(n = ...)to see more rows
State * Industry creates crossed groups:
1. Australian total
2. Totals by state
3. Totals by industry
4. State-by-industry series
retail_grouped <- retail_data |> + aggregate_key( + State * Industry, + Turnover = sum(Turnover) + )
print(retail_grouped) # A tsibble: 13,452 x 4 [1M] # Key: State, Industry [59] Month State Industry Turnover
<chr> <chr> 1 2000 Jan 10720. 2 2000 Feb 10042. 3 2000 Mar 10854. 4 2000 Apr 10533. 5 2000 May 10972. 6 2000 Jun 11592. 7 2000 Jul 10528. 8 2000 Aug 11106. 9 2000 Sep 11364. 10 2000 Oct 11580. # ℹ 13,442 more rows # ℹ Use print(n = ...)to see more rows
group_structure <- retail_grouped |> + distinct(State, Industry) |> + 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 × Industry”, +
+ TRUE ~ “Other” + ) + ) |> + count(Level, name = “Number_of_series”)print(group_structure) # A tibble: 4 × 2 Level Number_of_series
1 Industry 6 2 State 8 3 State × Industry 44 4 Total 1 group_structure |> + ggplot( + aes( + x = reorder(Level, Number_of_series), + y = Number_of_series + ) + ) + + geom_col() + + geom_text( + aes(label = Number_of_series), + vjust = -0.3 + ) + + labs( + title = “Structure of the Grouped Retail Dataset”, + subtitle = “State and industry are crossed grouping dimensions”, + x = “Grouping level”, + y = “Number of series” + ) + + theme_minimal()
state_totals <- retail_grouped |> + filter( + !is_aggregated(State), + is_aggregated(Industry) + ) |> + mutate(State = as.character(State))
state_totals |> + ggplot( + aes( + x = Month, + y = Turnover, + group = State, + linetype = State + ) + ) + + geom_line() + + labs( + title = “Retail Turnover by State”, + subtitle = “The six selected industries are aggregated within each state”, + x = “Month”, + y = “Turnover, AU$ millions”, + linetype = “State” + ) + + theme_minimal() + + theme(legend.position = “bottom”)
industry_totals <- retail_grouped |> + filter( + is_aggregated(State), + !is_aggregated(Industry) + ) |> + mutate(Industry = as.character(Industry))
industry_totals |> + ggplot( + aes( + x = Month, + y = Turnover, + group = Industry, + linetype = Industry + ) + ) + + geom_line() + + labs( + title = “Australian Retail Turnover by Industry”, + subtitle = “Each series is aggregated across all states and territories”, + x = “Month”, + y = “Turnover, AU$ millions”, + linetype = “Industry” + ) + + theme_minimal() + + theme(legend.position = “bottom”)
bottom_level <- retail_grouped |> + filter( + !is_aggregated(State), + !is_aggregated(Industry) + )
latest_month <- max(bottom_level$Month)
latest_group_values <- bottom_level |> + filter(Month == latest_month) |> + mutate( + State = as.character(State), + Industry = as.character(Industry) + )
latest_group_values |> + ggplot( + aes( + x = Industry, + y = State, + fill = Turnover + ) + ) + + geom_tile() + + geom_text( + aes(label = round(Turnover, 0)), + size = 3 + ) + + scale_x_discrete( + labels = function(x) stringr::str_wrap(x, width = 20) + ) + + labs( + title = “State-by-Industry Retail Turnover”, + subtitle = paste(“Turnover during”, latest_month), + x = “Industry”, + y = “State”, + fill = “AU$ millions” + ) + + theme_minimal() + + theme( + axis.text.x = element_text( + angle = 45, + hjust = 1 + ) + )
Dynamically select the final 24 months as the test set.
test_months <- retail_grouped |> + as_tibble() |> + distinct(Month) |> + arrange(Month) |> + slice_tail(n = 24) |> + pull(Month)
retail_train <- retail_grouped |> + filter(!Month %in% test_months)
retail_test <- retail_grouped |> + filter(Month %in% test_months)
cat(“period:”)
Training period: > print(range(retail_train$Month)) <yearmonth[2]> [1] “2000 Jan” “2016 Dec” > > cat(“period:”)
Test period: > print(range(retail_test\(Month))
<yearmonth[2]>
[1] "2017 Jan" "2018 Dec"
>
>
>
> # These are called "flat" base models because every
series is
> # modeled independently. No grouping constraints are imposed.
>
> retail_fit <- retail_train |>
+ model(
+ ets_flat = ETS(Turnover),
+ arima_flat = ARIMA(Turnover)
+ ) |>
+
+ # Reconcile the independently generated forecasts.
+ reconcile(
+ # Bottom-up uses the state-by-industry forecasts.
+ ets_bottom_up = bottom_up(ets_flat),
+ arima_bottom_up = bottom_up(arima_flat),
+
+ # OLS uses forecasts from every grouping level but gives
+ # them equal error-variance treatment.
+ ets_ols = min_trace(
+ ets_flat,
+ method = "ols"
+ ),
+
+ arima_ols = min_trace(
+ arima_flat,
+ method = "ols"
+ ),
+
+ # MinT shrinkage accounts for estimated forecast-error
+ # variances and correlations.
+ ets_mint = min_trace(
+ ets_flat,
+ method = "mint_shrink"
+ ),
+
+ arima_mint = min_trace(
+ arima_flat,
+ method = "mint_shrink"
+ )
+ )
>
> print(retail_fit)
# A mable: 59 x 10
# Key: State, Industry [59]
State Industry ets_flat
<chr*> <chr*> <model>
1 Australian Capital Territory Cafes, restaurants and ta…
<ETS(M,Ad,M)>
2 Australian Capital Territory Clothing, footwear and pe…
<ETS(M,Ad,M)>
3 Australian Capital Territory Department stores …
<ETS(A,Ad,A)>
4 Australian Capital Territory Food retailing …
<ETS(M,Ad,M)>
5 Australian Capital Territory Household goods retailing…
<ETS(M,Ad,M)>
6 Australian Capital Territory Other
retailing … <ETS(M,A,M)>
7 Australian Capital Territory
<aggregated> <ETS(M,Ad,M)>
8 New South Wales Cafes, restaurants and ta…
<ETS(M,Ad,M)>
9 New South Wales Clothing, footwear and
pe… <ETS(M,A,M)>
10 New South Wales Department stores …
<ETS(A,Ad,A)>
# ℹ 49 more rows
# ℹ 7 more variables: arima_flat <model>, ets_bottom_up
<model>,
# arima_bottom_up <model>, ets_ols <model>, arima_ols
<model>,
# ets_mint <model>, arima_mint <model>
# ℹ Use `print(n = ...)` to see more rows
>
>
>
> retail_forecasts <- retail_fit |>
+ forecast(h = 24)
>
> print(retail_forecasts)
# A fable: 11,328 x 6 [1M]
# Key: State, Industry, .model [472]
State Industry .model Month
<chr*> <chr*> <chr> <mth>
1 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Jan
2 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Feb
3 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Mar
4 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Apr
5 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 May
6 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Jun
7 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Jul
8 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Aug
9 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Sep
10 Australian Capital Territory Cafes, restaurants and … ets_f… 2017 Oct
# ℹ 11,318 more rows
# ℹ 2 more variables: Turnover <dist>, .mean <dbl>
# ℹ Use `print(n = ...)` to see more rows
>
>
>
> history_months <- retail_grouped |>
+ as_tibble() |>
+ distinct(Month) |>
+ arrange(Month) |>
+ slice_tail(n = 60) |>
+ pull(Month)
>
> national_actual <- retail_grouped |>
+ filter(
+ is_aggregated(State),
+ is_aggregated(Industry),
+ Month %in% history_months
+ )
>
> national_forecasts <- retail_forecasts |>
+ filter(
+ is_aggregated(State),
+ is_aggregated(Industry),
+ .model %in% c(
+ "ets_flat",
+ "ets_mint",
+ "arima_flat",
+ "arima_mint"
+ )
+ )
>
> ggplot() +
+ geom_line(
+ data = national_actual,
+ aes(
+ x = Month,
+ y = Turnover
+ ),
+ linewidth = 0.9
+ ) +
+ geom_line(
+ data = national_forecasts,
+ aes(
+ x = Month,
+ y = .mean,
+ linetype = .model
+ ),
+ linewidth = 0.8
+ ) +
+ labs(
+ title = "National Retail Turnover Forecasts",
+ subtitle = "Comparison of flat and MinT-reconciled
forecasts",
+ x = "Month",
+ y = "Turnover, AU\) millions”, + linetype = “Forecast
method” + ) + + theme_minimal() > > > > accuracy_individual
<- retail_forecasts |> + accuracy( + data = retail_test, +
measures = list( + RMSE = RMSE, + MAPE = MAPE + ) + ) |> + 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 ×
Industry”, +
+ TRUE ~ “Other” + ), +
+ Model_family = case_when( + stringr::str_detect(.model, “^ets”) ~
“ETS”, + stringr::str_detect(.model, “^arima”) ~ “ARIMA”, + TRUE ~
“Other” + ), +
+ Approach = case_when( + stringr::str_detect(.model, “flat”) ~ “Flat”,
+ stringr::str_detect(.model, “bottom_up”) ~ “Bottom-up”, +
stringr::str_detect(.model, “ols”) ~ “OLS”, +
stringr::str_detect(.model, “mint”) ~ “MinT”, + TRUE ~ “Other” + ) + )
> > print(accuracy_individual) # A tibble: 472 × 9 .model State
Industry .type RMSE MAPE Level Model_family
1 arima_bot… Australia… Cafes, re… Test 3.24 3.87 Stat… ARIMA
2 arima_bot… Australia… Clothing,… Test 2.13 5.58 Stat… ARIMA
3 arima_bot… Australia… Departmen… Test 1.50 3.64 Stat… ARIMA
4 arima_bot… Australia… Food reta… Test 11.6 5.27 Stat… ARIMA
5 arima_bot… Australia… Household… Test 5.52 4.82 Stat… ARIMA
6 arima_bot… Australia… Other ret… Test 4.74 7.37 Stat… ARIMA
7 arima_bot… Australia… <aggregat… Test 9.86 1.52 State ARIMA
8 arima_bot… New South… Cafes, re… Test 28.8 1.76 Stat… ARIMA
9 arima_bot… New South… Clothing,… Test 91.4 11.1 Stat… ARIMA
10 arima_bot… New South… Departmen… Test 17.0 2.78 Stat… ARIMA
# ℹ 462 more rows # ℹ 1 more variable: Approach print(n = ...) to see more rows > > > >
accuracy_by_level <- accuracy_individual |> + group_by( + Level, +
Model_family, + Approach, + .model + ) |> + summarise( + RMSE =
mean(RMSE, na.rm = TRUE), + MAPE = mean(MAPE, na.rm = TRUE), +
Number_of_series = n(), + .groups = “drop” + ) |> + mutate( + Level =
factor( + Level, + levels = c( + “Total”, + “State”, + “Industry”, +
“State × Industry” + ) + ) + ) |> + arrange(Level, Model_family,
RMSE) > > print(accuracy_by_level) # A tibble: 32 × 7 Level
Model_family Approach .model RMSE MAPE Number_of_series print(n = ...) to see more rows > > > >
accuracy_by_level |> + ggplot( + aes( + x = Approach, + y = RMSE, +
fill = Model_family + ) + ) + + geom_col( + position = “dodge” + ) + +
facet_wrap( + vars(Level), + scales = “free_y” + ) + + labs( + title =
“RMSE by Grouping Level”, + subtitle = “Flat models compared with
reconciled forecasts”, + x = “Forecasting approach”, + y = “Average
RMSE”, + fill = “Base model” + ) + + theme_minimal() > > > >
accuracy_by_level |> + ggplot( + aes( + x = Approach, + y = MAPE, +
fill = Model_family + ) + ) + + geom_col( + position = “dodge” + ) + +
facet_wrap( + vars(Level), + scales = “free_y” + ) + + labs( + title =
“MAPE by Grouping Level”, + subtitle = “Lower values indicate more
accurate forecasts”, + x = “Forecasting approach”, + y = “Average MAPE
(%)”, + fill = “Base model” + ) + + theme_minimal() > > > >
flat_benchmarks <- accuracy_by_level |> + filter(Approach ==
“Flat”) |> + select( + Level, + Model_family, + Flat_RMSE = RMSE, +
Flat_MAPE = MAPE + ) > > grouped_vs_flat <- accuracy_by_level
|> + filter(Approach != “Flat”) |> + left_join( + flat_benchmarks,
+ by = c( + “Level”, + “Model_family” + ) + ) |> + mutate( +
RMSE_improvement_percent = + 100 * (Flat_RMSE - RMSE) / Flat_RMSE,
+
+ MAPE_improvement_percent = + 100 * (Flat_MAPE - MAPE) / Flat_MAPE,
+
+ Lower_RMSE_than_flat = + RMSE < Flat_RMSE, +
+ Lower_MAPE_than_flat = + MAPE < Flat_MAPE, +
+ Outperforms_flat_on_both = + Lower_RMSE_than_flat & +
Lower_MAPE_than_flat + ) |> + arrange( + Level, + Model_family, +
desc(RMSE_improvement_percent) + ) > > print(grouped_vs_flat) # A
tibble: 24 × 14 Level Model_family Approach .model RMSE MAPE
Number_of_series print(n = ...) to
see more rows > > > > best_rmse_by_level <-
accuracy_by_level |> + group_by(Level) |> + slice_min( + order_by
= RMSE, + n = 1, + with_ties = FALSE + ) |> + ungroup() |> +
select( + Level, + Best_model = .model, + Model_family, + Approach, +
RMSE, + MAPE + ) > > print(best_rmse_by_level) # A tibble: 4 × 6
Level Best_model Model_family Approach RMSE MAPE
+ Sum_of_bottom_level = sum( + .mean[ + !is_aggregated(State) & +
!is_aggregated(Industry) + ], + na.rm = TRUE + ), +
+ Coherence_gap = + Total_forecast - Sum_of_bottom_level, +
+ .groups = “drop” + ) > > coherence_summary <- coherence_check
|> + group_by(.model) |> + summarise( + Maximum_absolute_gap = +
max(abs(Coherence_gap), na.rm = TRUE), +
+ Mean_absolute_gap = + mean(abs(Coherence_gap), na.rm = TRUE), +
+ .groups = “drop” + ) |> + arrange(Maximum_absolute_gap) > >
print(coherence_summary) # A tibble: 8 × 3 .model Maximum_absolute_gap
Mean_absolute_gap
Part 3 R Script:
library(rugarch) library(tidyverse)
data(“EuStockMarkets”)
dax_prices <- as.numeric(EuStockMarkets[, “DAX”])
dax_returns <- diff(log(dax_prices)) * 100
returns_data <- tibble( Observation = seq_along(dax_returns), Return = dax_returns, Squared_Return = dax_returns^2 )
print(head(returns_data))
returns_data |> ggplot(aes(x = Observation, y = Return)) + geom_line() + labs( title = “Daily DAX Percentage Log Returns”, subtitle = “Periods of large and small movements appear in clusters”, x = “Trading-day observation”, y = “Percentage log return” ) + theme_minimal()
returns_data |> ggplot(aes(x = Observation, y = Squared_Return)) + geom_line() + labs( title = “Squared DAX Returns”, subtitle = “Clusters of large squared returns indicate changing volatility”, x = “Trading-day observation”, y = “Squared return” ) + theme_minimal()
acf( dax_returns^2, lag.max = 30, main = “Autocorrelation of Squared DAX Returns” )
squared_return_test <- Box.test( dax_returns^2, lag = 20, type = “Ljung-Box” )
print(squared_return_test)
constant_model <- lm(dax_returns ~ 1)
unconditional_variance <- var( residuals(constant_model) )
unconditional_sd <- sqrt(unconditional_variance)
cat( “unconditional variance:”, unconditional_variance, “” )
cat( “Sample unconditional standard deviation:”, unconditional_sd, “” )
arch_specification <- ugarchspec( variance.model = list( model = “sGARCH”, garchOrder = c(1, 0) ), mean.model = list( armaOrder = c(0, 0), include.mean = TRUE ), distribution.model = “std” )
arch_fit <- ugarchfit( spec = arch_specification, data = dax_returns, solver = “hybrid” )
show(arch_fit)
garch_specification <- ugarchspec( variance.model = list( model = “sGARCH”, garchOrder = c(1, 1) ), mean.model = list( armaOrder = c(0, 0), include.mean = TRUE ), distribution.model = “std” )
garch_fit <- ugarchfit( spec = garch_specification, data = dax_returns, solver = “hybrid” )
show(garch_fit)
convergence_results <- tibble( Model = c(“ARCH(1)”, “GARCH(1,1)”), Convergence_Code = c( convergence(arch_fit), convergence(garch_fit) ) )
print(convergence_results)
arch_parameters <- as.data.frame( arch_fit@fit$matcoef ) |> rownames_to_column(“Parameter”)
names(arch_parameters) <- c( “Parameter”, “Estimate”, “Standard_Error”, “t_Value”, “p_Value” )
garch_parameters <- as.data.frame( garch_fit@fit$matcoef ) |> rownames_to_column(“Parameter”)
names(garch_parameters) <- c( “Parameter”, “Estimate”, “Standard_Error”, “t_Value”, “p_Value” )
cat(“parameter estimates:”) print(arch_parameters)
cat(“parameter estimates:”) print(garch_parameters)
arch_conditional_variance <- as.numeric(sigma(arch_fit))^2
garch_conditional_variance <- as.numeric(sigma(garch_fit))^2
variance_data <- tibble( Observation = seq_along(dax_returns), ARCH_Conditional_Variance = arch_conditional_variance, GARCH_Conditional_Variance = garch_conditional_variance, Constant_Unconditional_Variance = unconditional_variance )
print(head(variance_data))
variance_long <- variance_data |> pivot_longer( cols = c( ARCH_Conditional_Variance, GARCH_Conditional_Variance, Constant_Unconditional_Variance ), names_to = “Variance_Model”, values_to = “Variance” ) |> mutate( Variance_Model = recode( Variance_Model, ARCH_Conditional_Variance = “ARCH(1) conditional variance”, GARCH_Conditional_Variance = “GARCH(1,1) conditional variance”, Constant_Unconditional_Variance = “Constant unconditional variance” ) )
variance_long |> ggplot( aes( x = Observation, y = Variance, linetype = Variance_Model ) ) + geom_line() + labs( title = “Conditional Versus Unconditional Variance”, subtitle = paste( “The horizontal series assumes constant variance;”, “ARCH and GARCH allow variance to change” ), x = “Trading-day observation”, y = “Estimated variance”, linetype = “Variance estimate” ) + theme_minimal()
garch_mean <- as.numeric( fitted(garch_fit) )
garch_sigma <- as.numeric( sigma(garch_fit) )
volatility_bands <- tibble( Observation = seq_along(dax_returns), Return = dax_returns, Conditional_Mean = garch_mean, Upper_Band = garch_mean + garch_sigma, Lower_Band = garch_mean - garch_sigma )
volatility_bands |> ggplot(aes(x = Observation)) + geom_line( aes(y = Return), linewidth = 0.4 ) + geom_line( aes( y = Upper_Band, linetype = “Upper conditional SD” ) ) + geom_line( aes( y = Lower_Band, linetype = “Lower conditional SD” ) ) + labs( title = “DAX Returns and GARCH Conditional Volatility Bands”, subtitle = “The bands widen during periods of elevated volatility”, x = “Trading-day observation”, y = “Percentage return”, linetype = NULL ) + theme_minimal()
garch_coefficients <- coef(garch_fit)
alpha1 <- as.numeric( garch_coefficients[“alpha1”] )
beta1 <- as.numeric( garch_coefficients[“beta1”] )
manual_persistence <- alpha1 + beta1
garch_persistence <- as.numeric( persistence(garch_fit) )
garch_half_life <- as.numeric( halflife(garch_fit) )
garch_long_run_variance <- as.numeric( uncvariance(garch_fit) )
persistence_summary <- tibble( Measure = c( “ARCH shock coefficient, alpha1”, “Lagged variance coefficient, beta1”, “Total persistence, alpha1 + beta1”, “Package persistence estimate”, “Volatility shock half-life”, “GARCH long-run variance”, “Sample unconditional variance” ), Value = c( alpha1, beta1, manual_persistence, garch_persistence, garch_half_life, garch_long_run_variance, unconditional_variance ) )
print(persistence_summary)
arch_information <- infocriteria(arch_fit) garch_information <- infocriteria(garch_fit)
model_comparison <- tibble( Model = c(“ARCH(1)”, “GARCH(1,1)”), AIC = c( arch_information[1], garch_information[1] ), BIC = c( arch_information[2], garch_information[2] ), Shibata = c( arch_information[3], garch_information[3] ), Hannan_Quinn = c( arch_information[4], garch_information[4] ) )
print(model_comparison)
garch_standardized_residuals <- as.numeric( residuals( garch_fit, standardize = TRUE ) )
residual_test <- Box.test( garch_standardized_residuals, lag = 20, type = “Ljung-Box” )
squared_residual_test <- Box.test(
garch_standardized_residuals^2, lag = 20, type = “Ljung-Box” )
cat(“-Box test for standardized residuals:”) print(residual_test)
cat(“-Box test for squared standardized residuals:”) print(squared_residual_test)