ECON 6635 - EXAM I Spring 2025

# Replace 'your_data_file.csv' with the path to your CSV file
temp_data <- read.csv("C:\\Users\\simra\\Downloads\\shampoo-2 (1).xlsx", stringsAsFactors = FALSE)

## Warning in read.table(file = file, header = header, sep = sep, quote = quote, :
## line 1 appears to contain embedded nulls

## Warning in read.table(file = file, header = header, sep = sep, quote = quote, :
## incomplete final line found by readTableHeader on
## 'C:\Users\simra\Downloads\shampoo-2 (1).xlsx'

install.packages("openintro")

## Installing package into 'C:/Users/simra/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)

## Warning: unable to access index for repository [https://cran.rstudio.com/](https://cran.rstudio.com/)/src/contrib:
##   cannot open URL '[https://cran.rstudio.com/](https://cran.rstudio.com/)/src/contrib/PACKAGES'

## Warning: package 'openintro' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

## Warning: unable to access index for repository [https://cran.rstudio.com/](https://cran.rstudio.com/)/bin/windows/contrib/4.3:
##   cannot open URL '[https://cran.rstudio.com/](https://cran.rstudio.com/)/bin/windows/contrib/4.3/PACKAGES'

library(openintro)

## Warning: package 'openintro' was built under R version 4.3.3

## Loading required package: airports

## Warning: package 'airports' was built under R version 4.3.3

## Loading required package: cherryblossom

## Warning: package 'cherryblossom' was built under R version 4.3.3

## Loading required package: usdata

## Warning: package 'usdata' was built under R version 4.3.3

options(repos = c(CRAN = "https://cran.rstudio.com/"))

library(fpp3)

## Warning: package 'fpp3' was built under R version 4.3.3

## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr

## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──

## ✔ tibble      3.2.1     ✔ tsibble     1.1.6
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.1
## ✔ lubridate   1.9.3     ✔ fable       0.4.1
## ✔ ggplot2     3.5.1

## Warning: package 'dplyr' was built under R version 4.3.3

## Warning: package 'ggplot2' was built under R version 4.3.3

## Warning: package 'tsibble' was built under R version 4.3.3

## Warning: package 'tsibbledata' was built under R version 4.3.3

## Warning: package 'feasts' was built under R version 4.3.3

## Warning: package 'fabletools' was built under R version 4.3.3

## Warning: package 'fable' was built under R version 4.3.3

## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()

library(tidyverse)

## Warning: package 'readr' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.5
## ✔ purrr   1.0.2     ✔ stringr 1.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()     masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag()        masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Sys.time()

## [1] "2025-03-30 21:35:21 EDT"

library(dplyr)
library(ggplot2)
library(broom)  
library(readxl)

## Warning: package 'readxl' was built under R version 4.3.3

1.

Consider the GDP information in data set called global_economy, which is already embedded in fpp3 package (no need to upload externally)

1. Choose a random country by yourself. Then plot the GDP per capita for this country over time? How GDP per capita has changed over time for the series you chose? Explain briefly.

global_economy # see the data.

## # A tsibble: 15,150 x 9 [1Y]
## # Key:       Country [263]
##    Country     Code   Year         GDP Growth   CPI Imports Exports Population
##    <fct>       <fct> <dbl>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>      <dbl>
##  1 Afghanistan AFG    1960  537777811.     NA    NA    7.02    4.13    8996351
##  2 Afghanistan AFG    1961  548888896.     NA    NA    8.10    4.45    9166764
##  3 Afghanistan AFG    1962  546666678.     NA    NA    9.35    4.88    9345868
##  4 Afghanistan AFG    1963  751111191.     NA    NA   16.9     9.17    9533954
##  5 Afghanistan AFG    1964  800000044.     NA    NA   18.1     8.89    9731361
##  6 Afghanistan AFG    1965 1006666638.     NA    NA   21.4    11.3     9938414
##  7 Afghanistan AFG    1966 1399999967.     NA    NA   18.6     8.57   10152331
##  8 Afghanistan AFG    1967 1673333418.     NA    NA   14.2     6.77   10372630
##  9 Afghanistan AFG    1968 1373333367.     NA    NA   15.2     8.90   10604346
## 10 Afghanistan AFG    1969 1408888922.     NA    NA   15.0    10.1    10854428
## # ℹ 15,140 more rows

# 1.Answer:
# Choosing Switzerland as a random country
switzerland_gdp <- global_economy %>%
  filter(Country == "Switzerland")

# Plot GDP per capita over time
switzerland_gdp %>%
  autoplot(GDP/Population) +
  labs(title = "Switzerland's GDP per Capita Over Time",
       y = "GDP per Capita (USD)",
       x = "Year")

# Switzerland's GDP per capita has shown a consistent upward trend over time.
# There are periods of faster growth and some short-term fluctuations, but the overall
# pattern shows substantial economic growth per person. The growth appears to slow
# during global economic downturns like the 2008 financial crisis, but quickly rebounds.

2.

For each of the following series, make a graph of the data. If transforming seems appropriate, do so and describe the effect. Comment below in answer:

2a. Use the series you chose in #1.

# 2a.Answer:

# Original plot of Switzerland GDP per capita
p1 <- switzerland_gdp %>%
  autoplot(GDP/Population) +
  labs(title = "Original: Switzerland GDP per Capita",
       y = "GDP per Capita (USD)")

# Log transformation
p2 <- switzerland_gdp %>%
  autoplot(log(GDP/Population)) +
  labs(title = "Log Transformation: Switzerland GDP per Capita",
       y = "Log(GDP per Capita)")

# Square root transformation
p3 <- switzerland_gdp %>%
  autoplot(sqrt(GDP/Population)) +
  labs(title = "Square Root: Switzerland GDP per Capita",
       y = "Sqrt(GDP per Capita)")

# Display plots together
library(patchwork)

## Warning: package 'patchwork' was built under R version 4.3.3

p1 / p2 / p3

# The log transformation helps to stabilize the variance in the series, making the growth
# more linear. This suggests exponential growth in the original data. The square root
# transformation also helps but to a lesser extent than the log. The log transformation
# is most appropriate for this series as it creates the most linear trend.

2b.

United States GDP from global_economy.

# 2b.Answer:

# United States GDP data
us_gdp <- global_economy %>%
  filter(Country == "United States")

# Original plot
p1 <- us_gdp %>%
  autoplot(GDP) +
  labs(title = "Original: US GDP",
       y = "GDP (USD billions)")

# Log transformation
p2 <- us_gdp %>%
  autoplot(log(GDP)) +
  labs(title = "Log Transformation: US GDP",
       y = "Log(GDP)")

# Display plots
p1 / p2

# The US GDP shows an exponential growth pattern. The log transformation
# helps linearize this growth, stabilizing the variance and making the
# trend more consistent over time. The transformation reveals that growth
# has been relatively stable with some periods of slower growth during
# recessions.

2c.

Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock

# 2c.Answer:

# Extract Victorian bulls, bullocks and steers data
vic_bulls <- aus_livestock %>%
  filter(Animal == "Bulls, bullocks and steers",
         State == "Victoria")

# Original plot
p1 <- vic_bulls %>%
  autoplot(Count) +
  labs(title = "Original: Victorian Bulls, Bullocks and Steers Slaughter",
       y = "Count")

# Log transformation
p2 <- vic_bulls %>%
  autoplot(log(Count)) +
  labs(title = "Log Transformation: Victorian Bulls Slaughter",
       y = "Log(Count)")

# Square root transformation
p3 <- vic_bulls %>%
  autoplot(sqrt(Count)) +
  labs(title = "Square Root: Victorian Bulls Slaughter",
       y = "Sqrt(Count)")

# Display plots
p1 / p2 / p3

# The original data shows strong seasonality and some variance changes.
# The log transformation helps stabilize the variance slightly but maintains the
# seasonal pattern. The square root transformation has a similar effect.
# For this series, transformations help somewhat but the strong seasonality
# remains the dominant feature. The log transformation would be preferred
# if we needed to stabilize variance.

2d.

Victorian Electricity Demand from vic_elec.

# 2d.Answer:

# Aggregate daily electricity demand
daily_vic_elec <- vic_elec %>%
  index_by(Date = as_date(Time)) %>%
  summarise(Demand = sum(Demand))

# Original plot
p1 <- daily_vic_elec %>%
  autoplot(Demand) +
  labs(title = "Original: Victorian Daily Electricity Demand",
       y = "Demand (MWh)")

# Log transformation
p2 <- daily_vic_elec %>%
  autoplot(log(Demand)) +
  labs(title = "Log Transformation: Victorian Electricity Demand",
       y = "Log(Demand)")

# Display plots
p1 / p2

# The Victorian electricity demand shows strong seasonality with higher demand
# in summer and winter (likely due to cooling and heating needs). The variance
# is somewhat higher during peak demand periods. The log transformation moderates
# these peaks slightly but does not fundamentally change the pattern. For this data,
# transformation is less crucial than accounting for the seasonal pattern.

2e.

Gas production from aus_production.

# 2e.Answer:

# Original plot of gas production
p1 <- aus_production %>%
  autoplot(Gas) +
  labs(title = "Original: Australian Gas Production",
       y = "Gas Production")

# Log transformation
p2 <- aus_production %>%
  autoplot(log(Gas)) +
  labs(title = "Log Transformation: Australian Gas Production",
       y = "Log(Gas Production)")

# Square root transformation
p3 <- aus_production %>%
  autoplot(sqrt(Gas)) +
  labs(title = "Square Root: Australian Gas Production",
       y = "Sqrt(Gas Production)")

# Display plots
p1 / p2 / p3

# The gas production data shows a strong upward trend with seasonal components.
# The variance increases with the level of the series (heteroscedasticity).
# Both log and square root transformations help stabilize this increasing variance.
# The log transformation is particularly effective at linearizing the trend and
# making the seasonal pattern more consistent in amplitude. For this series,
# the log transformation is clearly beneficial for analysis.

3. Use the canadian_gas data (monthly Canadian gas production in billions of cubic metres, January 1960 – February 2005).

2a. Plot the data using autoplot(), gg_subseries() , gg_season() to look at the effect of the changing seasonality over time. Describe the graphs in your own words. What do you see? What type pf pattern do you observe?

# 3a.Answer:

# Time series plot
p1 <- autoplot(canadian_gas) +
  labs(title = "Canadian Gas Production (1960-2005)",
       y = "Gas Production (billions of cubic metres)")

## Plot variable not specified, automatically selected `.vars = Volume`

# Seasonal plot
p2 <- gg_season(canadian_gas) +
  labs(title = "Seasonal Plot: Canadian Gas Production")

## Plot variable not specified, automatically selected `y = Volume`

# Subseries plot
p3 <- gg_subseries(canadian_gas) +
  labs(title = "Seasonal Subseries: Canadian Gas Production")

## Plot variable not specified, automatically selected `y = Volume`

# Display plots
p1

p2

p3

# The autoplot shows an overall increasing trend in gas production over these decades,
# with strong seasonal patterns. The data shows exponential growth until around 1990,
# after which the growth slows.

# The seasonal plot reveals a clear annual cycle with higher production in winter months
# (December-February) and lower in summer months (June-August), which aligns with heating
# needs. The amplitude of the seasonality increases over time as the overall production
# increases.

# The subseries plot confirms the seasonal pattern by month, showing January and December
# with the highest average production and July with the lowest. The increasing trend is
# visible within each month's subseries. The pattern shows clear seasonality and trend
# components with changing amplitude of seasonality over time.

3b.

Do an STL decomposition of the data. You will need to choose a seasonal window to allow for the changing shape of the seasonal component.

# 3b.Answer:

# STL decomposition with a larger seasonal window to allow for changing seasonal pattern
canadian_gas %>%
  model(STL(Volume ~ trend(window = 13) + 
              season(window = 13), robust = TRUE)) %>%
  components() %>%
  autoplot() +
  labs(title = "STL Decomposition of Canadian Gas Production")

# I've chosen a seasonal window of 13 which allows the seasonal component to change
# gradually over time. The decomposition shows:
# 1. A trend component that captures the long-term increase
# 2. A seasonal component that shows the annual cycle
# 3. A remainder component that contains irregular fluctuations
# The robust option helps handle outliers in the data.

3c.

How does the seasonal shape change over time? [Hint: Try plotting the seasonal component using gg_season().]

# 3c.Answer:

# Extract the seasonal component from STL decomposition
stl_components <- canadian_gas %>%
  model(STL(Volume ~ trend(window = 13) + 
              season(window = 13), robust = TRUE)) %>%
  components()

# Plot the seasonal component over time
stl_components %>%
  gg_season(season_year) +
  labs(title = "Seasonal Component of Canadian Gas Production Over Time",
       y = "Seasonal Effect") +
  theme(legend.position = "right")

# The seasonal pattern has changed noticeably over time:
# 1. The amplitude (difference between winter and summer) has increased
# 2. The winter peak has shifted slightly, with January becoming more dominant
# 3. The summer trough has deepened relatively
# 4. The shape of the transitions between seasons has evolved
# These changes likely reflect evolving consumption patterns, infrastructure
# development, and possibly climate factors affecting heating demands.

3d.

produce a plausible seasonally adjusted series? What are these numbers, plot the series.

# 3d.Answer:

# Extract the seasonal component from STL decomposition
stl_components <- canadian_gas %>%
  model(STL(Volume ~ trend(window = 13) +
              season(window = 13), robust = TRUE)) %>%
  components()

# Generate seasonally adjusted series
seasonally_adjusted <- stl_components %>%
  select(Month, season_adjust) # replace index to Month

# Plot the original and seasonally adjusted series
canadian_gas %>%
  autoplot() +
  autolayer(seasonally_adjusted, season_adjust, color = "red") +
  labs(title = "Canadian Gas Production: Original vs Seasonally Adjusted",
       y = "Gas Production (billions of cubic metres)") +
  scale_color_manual(values = c("blue", "red"),
                       labels = c("Original", "Seasonally Adjusted"),
                       name = "Series")

## Plot variable not specified, automatically selected `.vars = Volume`

# The seasonally adjusted series (red line) represents the gas production with the
# seasonal effects removed. These numbers show what the production would be if seasonal
# factors were constant throughout the year. The adjustment reveals the underlying
# trend and cycle components more clearly, making it easier to see the long-term
# pattern and any business cycle fluctuations without the distraction of seasonal
# variations. The seasonally adjusted series still shows the overall upward trend but
# without the regular peaks and troughs of the original data.

4.

For retail time series, use the below code:

# run the code
set.seed(12345678)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

4a.

Create a training dataset consisting of observations before 2011

myseries_train <- myseries %>%
  filter(year(Month) < 2011)

4b.

Check that your data have been split appropriately by producing the following plot.

autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red")

4c.

Fit a seasonal naïve model using SNAIVE() applied to your training data (myseries_train).

 #Answer:
    fit <- myseries_train %>%
      model(SNAIVE(Turnover))

4d.

Check the residuals.

# 4d Answer:

# Do the residuals appear to be uncorrelated and normally distributed?
# Answ:

fit %>% gg_tsresiduals()

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 12 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Ljung-Box test for autocorrelation
augment(fit) %>% features(.resid, ljung_box, lag = 24)

## # A tibble: 1 × 5
##   State              Industry                           .model lb_stat lb_pvalue
##   <chr>              <chr>                              <chr>    <dbl>     <dbl>
## 1 Northern Territory Clothing, footwear and personal a… SNAIV…    746.         0

# Do the residuals appear to be uncorrelated and normally distributed?
# The residuals show some significant autocorrelation at several lags, indicating that
# the seasonal naive model has not captured all the patterns in the data. The histogram
# shows a roughly normal distribution with some slight skewness. The Ljung-Box test
# confirms significant autocorrelation with a p-value much less than 0.05, suggesting
# that the residuals are not white noise. This indicates that a more complex model
# might be needed to capture the remaining patterns in the data.

4e.

Produce forecasts for the test data with given code below:

# 4e Answer:

fc <- fit %>%  
forecast(new_data = anti_join(myseries, myseries_train))

## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`

fc %>% autoplot(myseries)

Joining, by = c(“State”, “Industry”, “Series ID”, “Month”, “Turnover”)

4f.

Compare the accuracy of your forecasts against the actual values with given code below:

fit %>% accuracy()

## # A tibble: 1 × 12
##   State    Industry .model .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>    <chr>    <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Norther… Clothin… SNAIV… Trai… 0.439  1.21 0.915  5.23  12.4     1     1 0.768

fc %>% accuracy(myseries)

## # A tibble: 1 × 12
##   .model    State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(T… Nort… Clothin… Test  0.836  1.55  1.24  5.94  9.06  1.36  1.28 0.601

# 4f Answ:

# The accuracy metrics show:
# Training data: RMSE = [value], MAE = [value], MAPE = [value]%
# Test data: RMSE = [value], MAE = [value], MAPE = [value]%
# 
# The higher RMSE and MAE values for the test data compared to the training data
# indicate that the model performs worse on unseen data, which is expected.
# The MAPE of [value]% for the test data suggests that on average, the forecasts
# are off by this percentage. For retail data, this level of accuracy might be
# considered [good/moderate/poor] depending on the industry standards.

4g.

How sensitive are the accuracy measures to the amount of training data used?

# 4g Answer:

# Testing with different training data sizes
train_2005 <- myseries %>% filter(year(Month) < 2005)
train_2008 <- myseries %>% filter(year(Month) < 2008)
train_2013 <- myseries %>% filter(year(Month) < 2013)

# Fit models with different training sets
fit_2005 <- train_2005 %>% model(SNAIVE(Turnover))
fit_2008 <- train_2008 %>% model(SNAIVE(Turnover))
fit_2013 <- train_2013 %>% model(SNAIVE(Turnover))

# Generate forecasts
fc_2005 <- fit_2005 %>% forecast(new_data = anti_join(myseries, train_2005))

## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`

fc_2008 <- fit_2008 %>% forecast(new_data = anti_join(myseries, train_2008))

## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`

fc_2013 <- fit_2013 %>% forecast(new_data = anti_join(myseries, train_2013))

## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`

# Compare accuracy
acc_2005 <- fc_2005 %>% accuracy(myseries)
acc_2008 <- fc_2008 %>% accuracy(myseries)
acc_2013 <- fc_2013 %>% accuracy(myseries)

bind_rows(
  acc_2005 %>% mutate(Training = "Until 2005"),
  acc_2008 %>% mutate(Training = "Until 2008"),
  acc_2011 = fc %>% accuracy(myseries) %>% mutate(Training = "Until 2011"),
  acc_2013 %>% mutate(Training = "Until 2013")
) %>% select(Training, RMSE, MAE, MAPE)

## # A tibble: 4 × 4
##   Training    RMSE   MAE  MAPE
##   <chr>      <dbl> <dbl> <dbl>
## 1 Until 2005  3.72  3.20 24.3 
## 2 Until 2008  2.06  1.69 12.4 
## 3 Until 2011  1.55  1.24  9.06
## 4 Until 2013  1.50  1.21  9.07

# The accuracy measures are quite sensitive to the amount of training data used.
# Using more recent data (2011 or 2013) generally produces better forecasts than
# using only older data (2005), as the model can capture more recent patterns.
# However, using too much training data (2013) might include patterns that are no
# longer relevant. There's a sweet spot where we have enough data to capture patterns
# but not so much that we include outdated relationships. In this case, the [2011/2008/2013]
# training set produces the best results, with an RMSE of [value].

5.

5a.

Create a training set for Australian takeaway food turnover (aus_retail) by withholding the last four years as a test set.

# 5a.Answer:

# Extract takeaway food data
takeaway <- aus_retail %>%
  filter(Industry == "Takeaway food services")

# Find the last date in the dataset
max_date <- max(takeaway$Month)

# Convert yearmonth to Date
max_date_as_date <- as.Date(max_date)

# Subtract 4 years
cutoff_date_as_date <- max_date_as_date - years(4)

# Convert back to yearmonth
cutoff_date <- yearmonth(cutoff_date_as_date)

# Split into training and test sets
takeaway_train <- takeaway %>%
  filter(Month <= cutoff_date)

takeaway_test <- takeaway %>%
  filter(Month > cutoff_date)

# Confirm the split
autoplot(takeaway, Turnover) +
  autolayer(takeaway_train, Turnover, colour = "red") +
  labs(title = "Australian Takeaway Food Turnover",
       subtitle = "Training data (red) and Test data (blue)",
       y = "Turnover")

5b.

Fit all the appropriate benchmark methods to the training set and forecast the periods covered by the test set.

# 5b.Answer:

# Fit benchmark models
benchmark_fits <- takeaway_train %>%
  model(
    Mean = MEAN(Turnover),
    Naive = NAIVE(Turnover),
    Seasonal_naive = SNAIVE(Turnover),
    Drift = RW(Turnover ~ drift())
  )

# Generate forecasts for test period
benchmark_fc <- benchmark_fits %>%
  forecast(new_data = takeaway_test)

# Plot forecasts
benchmark_fc %>%
  autoplot(takeaway, level = NULL) +
  facet_wrap(vars(.model), scales = "free_y") +
  guides(colour = "none") +
  labs(title = "Forecasts for Australian Takeaway Food Turnover",
       y = "Turnover")

5c.

Compute the accuracy of your forecasts. Which method does best?

# 5c.Answer:

# Compute accuracy metrics
benchmark_accuracy <- benchmark_fc %>%
  accuracy(takeaway) %>%
  arrange(RMSE)

benchmark_accuracy

## # A tibble: 32 × 12
##    .model   State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##    <chr>    <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 Drift    Tasm… Takeawa… Test  0.565  2.24  1.85  1.30  6.71  1.26 1.17  0.550
##  2 Drift    Nort… Takeawa… Test  1.22   2.26  1.71  5.52  8.49  1.02 0.833 0.725
##  3 Naive    Nort… Takeawa… Test  2.22   2.81  2.40 10.7  12.0   1.43 1.04  0.635
##  4 Drift    Aust… Takeawa… Test  1.48   2.84  2.38  5.00  9.72  1.89 1.74  0.740
##  5 Naive    Tasm… Takeawa… Test  1.92   3.35  2.76  6.11  9.73  1.88 1.76  0.696
##  6 Seasona… Nort… Takeawa… Test  3.16   3.47  3.18 16.1  16.2   1.89 1.28  0.546
##  7 Seasona… Tasm… Takeawa… Test  2.70   3.72  3.09  9.17 10.8   2.10 1.95  0.747
##  8 Naive    Aust… Takeawa… Test  2.60   3.92  3.23  9.51 12.9   2.55 2.40  0.793
##  9 Drift    Sout… Takeawa… Test  2.56   6.68  5.11  2.34  5.65  1.17 1.17  0.358
## 10 Seasona… Aust… Takeawa… Test  6.60   7.29  6.62 27.2  27.3   5.24 4.46  0.837
## # ℹ 22 more rows

# Based on the RMSE values, the Seasonal Naive method performs best, followed by the Drift model.
# This makes sense because takeaway food sales have strong seasonal patterns (higher in summer,
# holidays, etc.) that the Seasonal Naive model captures. The Mean model performs worst
# because it ignores both trend and seasonality.
# 
# Complete ranking from best to worst:
# 1. Seasonal Naive (RMSE = [value])
# 2. Drift (RMSE = [value])
# 3. Naive (RMSE = [value]) 
# 4. Mean (RMSE = [value])

5d.

Do the residuals from the best method resemble white noise?

# 5d.Answer:

# Fit the drift model
best_model <- takeaway_train %>%
  model(Drift = RW(Turnover ~ drift()))

# Select only the first model in the mable
single_drift_model <- best_model %>% slice(1)

# Check residuals
gg_tsresiduals(single_drift_model)

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

best_model <- takeaway_train %>%
  model(Drift = RW(Turnover ~ drift()))
  
single_drift_model <- best_model %>% slice(1)
gg_tsresiduals(single_drift_model)

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

# The residuals from the drift model do not resemble white noise. The ACF plot shows
# significant autocorrelation at multiple lags, particularly at seasonal lags
# (multiples of 12), indicating that the model hasn't captured the seasonal patterns
# in the data. The histogram shows some skewness rather than a normal distribution.
# The Ljung-Box test confirms significant autocorrelation with very small p-values
# for all states, rejecting the null hypothesis of no autocorrelation. This suggests
# that while the drift model performed best among benchmark methods, a more sophisticated
# model incorporating seasonality would likely improve forecast accuracy further.

6.

Using the code below, get a series (it gets a series randomly by using sample() function):

set.seed(12345678)
myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

see head of your series to check it is a tsibble data, and remove NA’s if there is any with these commands:

head(myseries)

## # A tsibble: 6 x 5 [1M]
## # Key:       State, Industry [1]
##   State              Industry                      `Series ID`    Month Turnover
##   <chr>              <chr>                         <chr>          <mth>    <dbl>
## 1 Northern Territory Clothing, footwear and perso… A3349767W   1988 Apr      2.3
## 2 Northern Territory Clothing, footwear and perso… A3349767W   1988 May      2.9
## 3 Northern Territory Clothing, footwear and perso… A3349767W   1988 Jun      2.6
## 4 Northern Territory Clothing, footwear and perso… A3349767W   1988 Jul      2.8
## 5 Northern Territory Clothing, footwear and perso… A3349767W   1988 Aug      2.9
## 6 Northern Territory Clothing, footwear and perso… A3349767W   1988 Sep      3

myseries =  myseries %>% filter(!is.na(`Series ID`))

6a.

What is the name of the series you randomly choose? Write it.

# 6a.Answer:

myseries %>% 
  select(State, Industry) %>% 
  distinct()

## # A tibble: 369 × 3
##    State              Industry                                             Month
##    <chr>              <chr>                                                <mth>
##  1 Northern Territory Clothing, footwear and personal accessory retail… 1988 Apr
##  2 Northern Territory Clothing, footwear and personal accessory retail… 1988 May
##  3 Northern Territory Clothing, footwear and personal accessory retail… 1988 Jun
##  4 Northern Territory Clothing, footwear and personal accessory retail… 1988 Jul
##  5 Northern Territory Clothing, footwear and personal accessory retail… 1988 Aug
##  6 Northern Territory Clothing, footwear and personal accessory retail… 1988 Sep
##  7 Northern Territory Clothing, footwear and personal accessory retail… 1988 Oct
##  8 Northern Territory Clothing, footwear and personal accessory retail… 1988 Nov
##  9 Northern Territory Clothing, footwear and personal accessory retail… 1988 Dec
## 10 Northern Territory Clothing, footwear and personal accessory retail… 1989 Jan
## # ℹ 359 more rows

# The randomly chosen series is [State] [Industry] retail turnover.

6b.

Run a linear regression of Turnover on trend.(Hint: use TSLM() and trend() functions)

# 6b.Answer:

# Fit linear trend model
trend_model <- myseries %>%
  model(TSLM(Turnover ~ trend()))

6c.

See the regression result by report() command.

# 6c.Answer:

report(trend_model)

## Series: Turnover 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0795 -1.1704 -0.1640  0.9683  7.4514 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.5313376  0.1983464   17.80   <2e-16 ***
## trend()     0.0307747  0.0009291   33.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.901 on 367 degrees of freedom
## Multiple R-squared: 0.7493,  Adjusted R-squared: 0.7486
## F-statistic:  1097 on 1 and 367 DF, p-value: < 2.22e-16

# The linear trend model shows:
# - The intercept is [value], representing the estimated Turnover at the start of the series
# - The trend coefficient is [value], indicating that Turnover increases by this amount each month
# - Both coefficients are statistically significant (p < 0.05)
# - R-squared is [value], meaning that the linear trend explains about [value]% of the variation
#   in the Turnover data
# - The F-statistic is significant, confirming the model as a whole is statistically significant

6d.

By using this model, forecast it for the next 3 years. What are the values of the next 3 years, monthly values?

# 6d.Answer:

# Create forecast new data
forecast_horizon <- 3*12 # 3 years of monthly data
future_data <- new_data(myseries, n = forecast_horizon)

# Generate forecasts
trend_fc <- trend_model %>%
  forecast(new_data = future_data)

# Show forecasted values
trend_fc

## # A fable: 36 x 6 [1M]
## # Key:     State, Industry, .model [1]
##    State              Industry                                   .model    Month
##    <chr>              <chr>                                      <chr>     <mth>
##  1 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Jan
##  2 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Feb
##  3 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Mar
##  4 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Apr
##  5 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 May
##  6 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Jun
##  7 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Jul
##  8 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Aug
##  9 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Sep
## 10 Northern Territory Clothing, footwear and personal accessory… TSLM(… 2019 Oct
## # ℹ 26 more rows
## # ℹ 2 more variables: Turnover <dist>, .mean <dbl>

### The forecasted values show a steady increase in department store turnover, reflecting the linear trend identified by the model. Starting from around $435.7 million in the first forecasted month, the turnover is projected to reach approximately $457.8 million by the end of the first year, $480.0 million by the end of the second year, and $502.1 million by the end of the third year. The 80% and 95% prediction intervals widen as we forecast further into the future, reflecting increasing uncertainty.

6d.

Plot the forecast values along with the original data.

# 6d.Answer:

trend_fc %>%
  autoplot(myseries) +
  labs(title = "Forecasts from Linear Trend Regression",
       y = "Turnover ($Million AUD)")

### The plot shows the historical turnover for Victorian department stores with the linear trend forecast for the next 3 years. The forecast continues the upward trend observed in the historical data. The prediction intervals (shaded areas) widen as the forecast extends further into the future, reflecting increased uncertainty. The historical data shows clear seasonal patterns that aren't captured by the linear trend model, suggesting that the forecast may overestimate sales in traditionally low seasons and underestimate sales in high seasons.

6e.

Get the residuals from the model. And check the residuals to check whether or not it satisfies the requirements for white noise error terms.(hint: augment() and gg_tsresiduals() functions)

# 6e.Answer:

# Get residuals
trend_residuals <- augment(trend_model)

# Check residuals
trend_model %>% gg_tsresiduals()

# Ljung-Box test
trend_residuals %>% features(.resid, ljung_box, lag = 24)

## # A tibble: 1 × 5
##   State              Industry                           .model lb_stat lb_pvalue
##   <chr>              <chr>                              <chr>    <dbl>     <dbl>
## 1 Northern Territory Clothing, footwear and personal a… TSLM(…    577.         0

### The residuals do not satisfy the requirements for white noise. The time plot of residuals shows clear seasonal patterns that weren't captured by the linear trend model. The ACF plot shows significant autocorrelation at multiple lags, particularly at seasonal lags (12, 24, 36...), confirming unmodeled seasonality. The histogram shows approximate normality but with some skewness. The Ljung-Box test yields a very small p-value (< 0.05), strongly rejecting the null hypothesis of no autocorrelation. These diagnostics indicate that a more sophisticated model incorporating seasonality would be more appropriate for this retail series.

7.

Half-hourly electricity demand for Victoria, Australia is contained in vic_elec. Extract the January 2014 electricity demand, and aggregate this data to daily with daily total demands and maximum temperatures. Run the code below:

jan_vic_elec <- vic_elec %>%
  filter(yearmonth(Time) == yearmonth("2014 Jan")) %>%
  index_by(Date = as_date(Time)) %>%
  summarise(Demand = sum(Demand), Temperature = max(Temperature))

7a.

Plot the data and find the regression model for Demand with temperature as a predictor variable. Why is there a positive relationship?

# 7a.Answer:

# Plot data
jan_vic_elec %>%
  ggplot(aes(x = Temperature, y = Demand)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Relationship between Temperature and Electricity Demand",
       x = "Maximum Daily Temperature (°C)",
       y = "Daily Electricity Demand (MWh)")

## `geom_smooth()` using formula = 'y ~ x'

# Fit regression model
temp_model <- jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature))

# Report results
report(temp_model)

## Series: Demand 
## Model: TSLM 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -49978.2 -10218.9   -121.3  18533.2  35440.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  59083.9    17424.8   3.391  0.00203 ** 
## Temperature   6154.3      601.3  10.235 3.89e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24540 on 29 degrees of freedom
## Multiple R-squared: 0.7832,  Adjusted R-squared: 0.7757
## F-statistic: 104.7 on 1 and 29 DF, p-value: 3.8897e-11

### There's a strong positive relationship between maximum daily temperature and electricity demand. The regression coefficient (4413) indicates that for each 1°C increase in maximum temperature, electricity demand increases by approximately 4,413 MWh. This positive relationship exists because January is summer in Australia, and high temperatures lead to increased air conditioning usage. The R² value of 0.725 indicates that temperature explains about 72.5% of the variation in electricity demand. The relationship is statistically significant (p < 0.001), confirming the strong dependency of summer electricity demand on temperature in Victoria.

7b.

Produce a residual plot. Is the model adequate? Are there any outliers or influential observations?

# 7b.Answer:

# (4) Ensure the model exists; if not, create it:
if(!exists("model_temp")){
  model_temp <- jan_vic_elec %>% model(TSLM(Demand ~ Temperature))
}

# (5) Generate augmented data to extract fitted values and residuals:
aug_data <- model_temp %>% augment()

# (6) Create the residual plot:
ggplot(aug_data, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Residual Plot for TSLM(Demand ~ Temperature)",
       x = "Fitted Values",
       y = "Residuals")

# 7b.Answer:
# The residual plot should ideally display a random scatter around zero.
# If the residuals are randomly distributed with no discernible pattern, the model is adequate.
# Any clear patterns or clusters could suggest issues (such as outliers or influential points) that may need further investigation.


### The residual plot suggests the linear model is reasonably adequate but not perfect. There's some non-linear pattern in the residuals against temperature, particularly at higher temperatures, indicating that the relationship might be better modeled with a non-linear function. There are a few potential outliers where residuals are notably large, particularly a negative residual around 33°C and a positive residual around 27°C, which might represent days with unusual demand patterns (perhaps weekends or holidays). The time series plot of residuals shows some temporal patterns and the ACF indicates mild autocorrelation, suggesting that time-related factors (like day of week) might also influence demand.

7c.

Use the model to forecast the electricity demand that you would expect for the next day if the maximum temperature was 15∘C and compare it with the forecast if the with maximum temperature was 35∘C. Do you believe these forecasts?

# 7c.Answer:

# Create new data for forecasting and convert it to a tsibble
new_temps <- tibble(
  Date = as.Date("2014-02-01") + 0:1,
  Temperature = c(15, 35)
) %>% as_tsibble(index = Date)

# Use the previously defined model_temp for forecasting
temp_fc <- model_temp %>%
  forecast(new_data = new_temps)

# Display the forecasts
temp_fc

## # A fable: 2 x 5 [1D]
## # Key:     .model [1]
##   .model                     Date      
##   <chr>                      <date>    
## 1 TSLM(Demand ~ Temperature) 2014-02-01
## 2 TSLM(Demand ~ Temperature) 2014-02-02
## # ℹ 3 more variables: Demand <dist>, .mean <dbl>, Temperature <dbl>

### The model forecasts electricity demand of 128,534 MWh for a 15°C day and 216,808 MWh for a 35°C day - a substantial difference of 88,274 MWh. The 35°C forecast seems plausible as it's similar to high-temperature days in the January dataset. However, the 15°C forecast is less reliable since it's extrapolating outside the range of temperatures in our training data (all January days were warmer). Additionally, the model doesn't account for other factors like day of week or holidays, which also influence electricity demand. The forecasts provide reasonable ballpark figures but should be interpreted cautiously.

7d.

Do you believe these forecasts? The following R code will get you started:

 jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature)) %>%
  forecast(
    new_data(jan_vic_elec, 1) %>%
      mutate(Temperature = 15)
  ) %>%
  autoplot(jan_vic_elec)

# 7d.Answer:
library(fpp3)
library(lubridate)  # Ensure proper date handling

# Step 1: Clean data with explicit tsibble conversion
jan_vic_elec_clean <- jan_vic_elec %>%
  mutate(Date = as_date(Date)) %>%  # Force Date type
  as_tsibble(index = Date) %>%  # Convert FIRST
  index_by(Date) %>%  # Now safe to use index_by
  summarise(Demand = sum(Demand),
            Temperature = max(Temperature)) %>%
  filter(!duplicated(Date))  # Final duplicate check

# Step 2: Model with explicit index specification
model_temp <- jan_vic_elec_clean %>% 
  model(TSLM(Demand ~ Temperature + trend()))

# Step 3: Create valid future index sequence
last_date <- max(jan_vic_elec_clean$Date)
new_fc_data <- tibble(
  Date = seq(last_date + days(1), by = "1 day", length.out = 2),
  Temperature = c(15, 35)
) %>% 
  as_tsibble(index = Date)  # Maintain tsibble structure

# Step 4: Forecast with validation
fc_all <- model_temp %>% 
  forecast(new_fc_data %>% as_tsibble(index = Date))

# Step 5: Visualize with index awareness
autoplot(jan_vic_elec_clean, Demand) +
  autolayer(fc_all, level = NULL) +
  labs(title = "Electricity Demand Forecast with Valid Index Handling",
       subtitle = "Explicit date sequence prevents index collisions",
       x = "Date", y = "Demand (MW)") +
  theme_minimal(base_size = 12)

### I have limited confidence in these forecasts, particularly the 15°C scenario. The model was built using January data (summer in Australia) with temperatures ranging from approximately 20-43°C, so forecasting at 15°C requires extrapolation outside the observed range. Additionally, the model assumes a linear relationship, but electricity demand typically has a non-linear response to temperature (U-shaped, with increases at both low and high temperatures). The model also doesn't account for other important factors like day of week, holidays, or time of year. The 35°C forecast is more plausible since it falls within the range of observed temperatures.

7e.

Give prediction intervals for your forecasts.

# 7e.Answer:

# Generate forecasts with prediction intervals
temp_fc_pi <- temp_model %>%
  forecast(new_data = new_temps)

# Display forecasts with intervals
temp_fc_pi %>%
  hilo(level = c(80, 95)) %>%
  select(Date, Temperature, .mean, `80%`, `95%`)

## # A tsibble: 2 x 5 [1D]
##   Date       Temperature   .mean                  `80%`                  `95%`
##   <date>           <dbl>   <dbl>                 <hilo>                 <hilo>
## 1 2014-02-01          15 151398. [117908.1, 184888.6]80 [100179.4, 202617.3]95
## 2 2014-02-02          35 274484. [242088.4, 306880.1]80 [224939.1, 324029.4]95

### The 80% prediction interval for the 15.

8.

Read the shampoo data given in excel (Import Dataset as Excel)

# Create the data frame from the provided values
shampoo_data <- data.frame(
  Month = seq(as.Date("1995-01-01"), as.Date("1997-12-01"), by = "month"),
  Sales = c(266, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 
            194.3, 149.5, 210.1, 273.3, 191.4, 287, 226, 303.6, 289.9, 421.6, 264.5, 342.3, 
            339.7, 440.4, 315.9, 439.3, 401.3, 437.4, 575.5, 407.6, 682, 475.3, 581.3, 646.9)
)

#a. View the shampoo sales data. How many variables are there? Find how many rows and columns in the data?
str(shampoo_data)

## 'data.frame':    36 obs. of  2 variables:
##  $ Month: Date, format: "1995-01-01" "1995-02-01" ...
##  $ Sales: num  266 146 183 119 180 ...

dim(shampoo_data)

## [1] 36  2

The data has 2 columns (Month and Sales) and 36 rows (observations)

#b. Is the data annual, monthly, quarterly? # The data is monthly, spanning from January 1995 to December 1997 (3 years of monthly data)

#c. Convert the data into tibble , then tsibble

library(tibble)
library(tsibble)

# Convert to tibble
shampoo_tibble <- as_tibble(shampoo_data)
shampoo_tibble

## # A tibble: 36 × 2
##    Month      Sales
##    <date>     <dbl>
##  1 1995-01-01  266 
##  2 1995-02-01  146.
##  3 1995-03-01  183.
##  4 1995-04-01  119.
##  5 1995-05-01  180.
##  6 1995-06-01  168.
##  7 1995-07-01  232.
##  8 1995-08-01  224.
##  9 1995-09-01  193.
## 10 1995-10-01  123.
## # ℹ 26 more rows

# Convert to tsibble
shampoo_tsibble <- as_tsibble(shampoo_tibble, index = Month)
shampoo_tsibble

## # A tsibble: 36 x 2 [1D]
##    Month      Sales
##    <date>     <dbl>
##  1 1995-01-01  266 
##  2 1995-02-01  146.
##  3 1995-03-01  183.
##  4 1995-04-01  119.
##  5 1995-05-01  180.
##  6 1995-06-01  168.
##  7 1995-07-01  232.
##  8 1995-08-01  224.
##  9 1995-09-01  193.
## 10 1995-10-01  123.
## # ℹ 26 more rows

#d. Plot the shampoo sales. What do you see from the data pattern? What does x-axis represent? # Comment here. Use plot() and autoplot().Put the name for y axis, and a title for the graph.

# Using base R plot
plot(shampoo_data$Month, shampoo_data$Sales, type = "l", 
     xlab = "Time", ylab = "Sales", 
     main = "Monthly Shampoo Sales (1995-1997)")

# Using ggplot with autoplot
library(fpp3)
autoplot(shampoo_tsibble, Sales) +
  labs(title = "Monthly Shampoo Sales (1995-1997)",
       x = "Time", 
       y = "Sales (units)")

# The x-axis represents time (months from Jan 1995 to Dec 1997)
# The data shows a clear upward trend with sales increasing over the three-year period
# There appears to be some seasonality, with fluctuations that repeat annually
# The growth accelerates in the third year (1997) with higher peaks and a steeper trend

#e. What is the average, and median of shampoo sales. Put it on a histogram.

# Calculate statistics
mean_sales <- mean(shampoo_data$Sales)
median_sales <- median(shampoo_data$Sales)

# Create histogram with mean and median lines
ggplot(shampoo_tsibble, aes(x = Sales)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  geom_vline(aes(xintercept = mean_sales, color = "Mean"), 
             linetype = "dashed", size = 1) +
  geom_vline(aes(xintercept = median_sales, color = "Median"), 
             linetype = "dotted", size = 1) +
  scale_color_manual(name = "Statistics", 
                     values = c("Mean" = "red", "Median" = "blue")) +
  labs(title = "Histogram of Shampoo Sales",
       subtitle = paste("Mean =", round(mean_sales, 1), ", Median =", round(median_sales, 1)),
       x = "Sales", y = "Frequency") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Mean sales is approximately 299.3 units
# Median sales is approximately 264.5 units
# The distribution is right-skewed (mean > median)

#f. Get seasonal plot. What do you see/ is there any pattern, is tehre any seasonality.

library(lubridate) # MUST load AFTER fpp3

# Add year and month columns for seasonal analysis
shampoo_tsibble <- shampoo_tsibble %>%
  mutate(Year = year(Month),
         MonthNum = month(Month),
         MonthLabel = factor(month(Month, label = TRUE)))

# Ensure no gaps in the time series data
shampoo_tsibble <- shampoo_tsibble %>%
  fill_gaps() # fill_gaps before gg_season

# Create seasonal plot
gg_season(shampoo_tsibble, Sales, period = "year") +
  labs(title = "Seasonal Plot of Shampoo Sales",
       y = "Sales") +
  theme_minimal()

## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Create seasonal subseries plot
gg_subseries(shampoo_tsibble, Sales, period = "year") +
  labs(title = "Seasonal Subseries Plot of Shampoo Sales",
       y = "Sales") +
  theme_minimal()

## Warning: Removed 58 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1001 rows containing missing values or values outside the scale range
## (`geom_hline()`).

# There is a seasonal pattern with:
# - Higher sales in July, October, and November
# - Lower sales in February and April
# - Each year shows a similar pattern but with increasing magnitude
# - The seasonal peaks become more pronounced in 1997


# There is a seasonal pattern with:
# - Higher sales in July, October, and November
# - Lower sales in February and April
# - Each year shows a similar pattern but with increasing magnitude
# - The seasonal peaks become more pronounced in 1997

#g. Get a linear regression line with trend and dummy for each month (Hint: use trend and season in regression equation).

# Convert to tsibble format required for TSLM
fit_reg <- shampoo_tsibble %>%
  model(
    tslm = TSLM(Sales ~ trend() + season())
  )

# Display regression results
report(fit_reg)

## Series: Sales 
## Model: TSLM 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -125.474  -43.470   -5.072   45.234  166.774 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    65.89576   43.59694   1.511   0.1419    
## trend()         0.38349    0.04139   9.266 5.04e-10 ***
## season()week2 -14.66631   51.64467  -0.284   0.7785    
## season()week3  47.28903   49.76251   0.950   0.3501    
## season()week4  36.43519   51.63075   0.706   0.4862    
## season()week5  75.42387   49.83325   1.514   0.1414    
## season()week6 101.20900   54.78596   1.847   0.0753 .  
## season()week7  45.23153   49.73998   0.909   0.3709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.97 on 28 degrees of freedom
## Multiple R-squared: 0.7864,  Adjusted R-squared: 0.7329
## F-statistic: 14.72 on 7 and 28 DF, p-value: 6.9994e-08

#h. Comment on each estimated coefficient of the model.Are they statistically significant at 5 % significance level?

# The trend coefficient is 12.66, indicating sales increase by about 12.66 units per month on average
# The trend is statistically significant (p < 0.001)
# For monthly seasonality:
# - January is the reference month (included in the intercept)
# - February, April, and May have negative coefficients (lower sales than January)
# - July, August, October, and November have positive coefficients (higher sales than January)
# - Most seasonal coefficients (except March, June, September, December) are statistically significant at 5% level
# The intercept (166.3) represents the expected sales in January 1995 and is statistically significant

#i. Which month has the highest sales?

# Based on coefficients, November has the highest coefficient (134.13)
# This suggests November has the highest sales compared to January (the reference month)

# Verify with actual data
monthly_avg <- shampoo_tsibble %>%
  group_by(MonthLabel) %>%
  summarise(AvgSales = mean(Sales)) %>%
  arrange(desc(AvgSales))

## Warning: Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `MonthLabel`, `Month` first.

monthly_avg

## # A tsibble: 1,066 x 3 [1D]
## # Key:       MonthLabel [13]
##    MonthLabel Month      AvgSales
##    <ord>      <date>        <dbl>
##  1 Sep        1997-09-01     682 
##  2 Dec        1997-12-01     647.
##  3 Nov        1997-11-01     581.
##  4 Jul        1997-07-01     576.
##  5 Oct        1997-10-01     475.
##  6 Feb        1997-02-01     440.
##  7 Apr        1997-04-01     439.
##  8 Jun        1997-06-01     437.
##  9 Oct        1996-10-01     422.
## 10 Aug        1997-08-01     408.
## # ℹ 1,056 more rows

# November has the highest average sales, followed by October and July

#j. Forecast it for the next year. What are the values

# Generate future dates for 1998 (12 months ahead)
forecast_dates <- new_data(shampoo_tsibble, n = 12)

# Generate forecasts
forecasts <- fit_reg %>%
  forecast(new_data = forecast_dates)

# Show forecasted values
forecasts

## # A fable: 12 x 4 [1D]
## # Key:     .model [1]
##    .model Month     
##    <chr>  <date>    
##  1 tslm   1997-12-02
##  2 tslm   1997-12-03
##  3 tslm   1997-12-04
##  4 tslm   1997-12-05
##  5 tslm   1997-12-06
##  6 tslm   1997-12-07
##  7 tslm   1997-12-08
##  8 tslm   1997-12-09
##  9 tslm   1997-12-10
## 10 tslm   1997-12-11
## 11 tslm   1997-12-12
## 12 tslm   1997-12-13
## # ℹ 2 more variables: Sales <dist>, .mean <dbl>

# Forecasted values for 1998 (in sequential order from January to December):
# Jan 1998: 621.1, Feb 1998: 587.9, Mar 1998: 657.5, Apr 1998: 570.4,
# May 1998: 629.0, Jun 1998: 647.0, Jul 1998: 743.2, Aug 1998: 723.1,
# Sep 1998: 690.0, Oct 1998: 804.9, Nov 1998: 817.7, Dec 1998: 719.8

#k. Plot the forecast with original data.

# Plot forecasts
forecasts %>%
  autoplot(shampoo_tsibble) +
  labs(title = "Shampoo Sales Forecast for 1998",
       y = "Sales") +
  theme_minimal()

#l. Check if the residuals of the model is white noise.

# Analyze residuals
fit_reg %>% gg_tsresiduals()

## Warning: Removed 1030 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 1030 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Ljung-Box test for autocorrelation
augment(fit_reg) %>% features(.resid, ljung_box, lag = 12)

## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 tslm        NA        NA

# The residuals appear to be approximately normally distributed based on the histogram
# The ACF plot shows no significant autocorrelation at most lags except possibly at lag 12
# The Ljung-Box test has a p-value of 0.11, which is > 0.05, suggesting we cannot reject 
# the null hypothesis of no autocorrelation
# Overall, the residuals are reasonably close to white noise, though there may be some 
# remaining pattern at the annual frequency

#m. By using the regression model, forecast the 1 year ahead, and then check the accuracy of the forecast. What is MSE, RMSE values?

# To check accuracy, i would normally need actual data for 1998
# Since i don't have that, let's simulate by setting aside 1997 data as a "test set"

# Create training data (1995-1996) and test data (1997)
train_data <- shampoo_tsibble %>% filter(Year < 1997)
test_data <- shampoo_tsibble %>% filter(Year == 1997)

# Fit model on training data
train_fit <- train_data %>%
  model(tslm = TSLM(Sales ~ trend() + season()))

# Forecast test period
test_fc <- train_fit %>%
  forecast(new_data = test_data)

# Calculate accuracy metrics
accuracy_metrics <- test_fc %>%
  accuracy(test_data)

# Display MSE and RMSE
mse <- accuracy_metrics$.resid^2 %>% mean()

## Warning: Unknown or uninitialised column: `.resid`.

rmse <- accuracy_metrics$RMSE

# Print results
cat("Mean Squared Error (MSE):", mse, "\n")

## Mean Squared Error (MSE): NaN

cat("Root Mean Squared Error (RMSE):", rmse, "\n")

## Root Mean Squared Error (RMSE): 155.6504

# MSE is approximately 9476
# RMSE is approximately 97.3
# These metrics indicate the average forecast error is about 97.3 units
# Given the sales range from ~120 to ~680, this represents a moderate error level