ECON 6635 - EXAM I Spring 2025

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:

global_economy %>%
  filter(Country == "India") %>%
  autoplot(GDP / Population)

#India's GDP per capita has seen significant growth over time, with some years experiencing a decline. However, the overall trend has been one of rapid growth

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:
gdp <- global_economy %>%
  filter(Country == "India")

#GDP Per Capita Income Plot
autoplot(gdp, GDP/Population)

#Log transformation
gdp <- gdp %>%
  mutate(gdpPerCapita = GDP/Population, logGdp = log(gdpPerCapita))

#Log Transformation Plot
autoplot(gdp, logGdp)

gdp <- gdp %>%
  mutate(difference = difference(logGdp))

autoplot(gdp, difference)

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

#The data demonstrates an upward trend, signaling growth, while log transformation helps stabilize the values, and differencing eliminates the trends.

2b.

United States GDP from global_economy.

# 2b.Answer:
gdp <- global_economy %>%
  filter(Country == "United States")

#GDP Per Capita Plot
autoplot(gdp, GDP/Population)

#Log transformation
gdp <- gdp %>%
  mutate(gdpPerCapita = GDP/Population, logGdp = log(gdpPerCapita))

#Log transformation plot
autoplot(gdp, logGdp)

gdp <- gdp %>%
  mutate(difference = difference(logGdp))

autoplot(gdp, difference)

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

#By using the United States data, we can observe significant GDP growth, but applying a log transformation provides a clearer view of the trends.

2c.

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

# 2c.Answer:
slVic <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Bulls, bullocks and steers")

#Plot
autoplot(slVic, Count)

#Log transformation
logSlVic <- slVic %>%
  mutate(logC = log(Count))

#Plot 
autoplot(logSlVic, logC)

2d.

Victorian Electricity Demand from vic_elec.

# 2d.Answer:
autoplot(vic_elec, Demand)

logVic <- vic_elec %>%
  mutate(logD = log(Demand))

autoplot(logVic, logD)

2e.

Gas production from aus_production.

# 2e.Answer:
autoplot(aus_production, Gas)

logAusProd <- aus_production %>%
  mutate(logG = log(Gas))

autoplot(logAusProd, logG)

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:
autoplot(canadian_gas, Volume)

gg_subseries(canadian_gas, Volume)

gg_season(canadian_gas, Volume)

#The plots below clearly show a rapid increase overall, as well as on a monthly and seasonal basis.

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:

canadian_gas %>%
  model(STL(Volume ~ season(window=12))) %>%
  components() %>%
  autoplot(Volume)

3c.

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

# 3c.Answer:
canadian_gas %>%
  model(STL(Volume ~ season(window=12))) %>%
  components() %>%
  gg_season(Volume)

#Seasonally, we observe a decline in the summer and an increase in the winter.

3d.

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

# 3d.Answer:

canadian_gas %>%
  model(STL(Volume ~ season(window = 12))) %>%
  components() %>%
  mutate(seasonally_adjusted = Volume - season_year) %>%
  autoplot(seasonally_adjusted)

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))
fit

## # A mable: 1 x 3
## # Key:     State, Industry [1]
##   State              Industry                                 `SNAIVE(Turnover)`
##   <chr>              <chr>                                               <model>
## 1 Northern Territory Clothing, footwear and personal accesso…           <SNAIVE>

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()`).

#Yes, the residuals seem to be uncorrelated and follow a normal distribution. There is no obvious pattern, and their bell-shaped form is typical

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:
fc <- fit %>%  
  forecast(new_data = anti_join(myseries, myseries_train, by = c("State", "Industry", "Series ID", "Month", "Turnover")))

fc %>% 
  autoplot(myseries)

train <- fit %>% accuracy()

test <- fc %>% accuracy(myseries)

train

## # 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

test

## # 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

4g.

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

# 4g Answer:
#The accuracy metrics suggest that the model is sensitive to the training data, as it performs better on the training set than on the test set

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:
aus_takeaway <- aus_retail %>%
  filter(Industry == "Takeaway food services")

train_set <- aus_takeaway %>%
  filter(year(Month) <= 2015)

test_set <- aus_takeaway %>%
  filter(year(Month) > 2015)

5b.

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

# 5b.Answer:
fit_benchmarks <- train_set %>%
  model(
    naive = NAIVE(Turnover),
    snaive = SNAIVE(Turnover),
    mean = MEAN(Turnover),
    drift = RW(Turnover ~ drift())
  )

fc_benchmarks <- fit_benchmarks %>%
  forecast(data = test_set)

fc_benchmarks %>%
  autoplot(aus_takeaway) +
  autolayer(train_set) +
  autolayer(test_set)

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

5c.

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

# 5c.Answer:
fc_benchmarks %>%
  accuracy(test_set)

## # A tibble: 32 × 12
##    .model State  Industry .type      ME   RMSE    MAE      MPE  MAPE  MASE RMSSE
##    <chr>  <chr>  <chr>    <chr>   <dbl>  <dbl>  <dbl>    <dbl> <dbl> <dbl> <dbl>
##  1 drift  Austr… Takeawa… Test    2.27    2.89   2.57   8.84   10.3    NaN   NaN
##  2 drift  New S… Takeawa… Test  -11.1    38.1   29.8   -2.57    5.71   NaN   NaN
##  3 drift  North… Takeawa… Test   -1.98    2.60   2.23 -11.0    12.1    NaN   NaN
##  4 drift  Queen… Takeawa… Test  -30.7    35.3   31.3  -10.2    10.4    NaN   NaN
##  5 drift  South… Takeawa… Test  -10.3    11.9   11.0  -11.9    12.6    NaN   NaN
##  6 drift  Tasma… Takeawa… Test    0.193   2.15   1.76   0.0342  6.25   NaN   NaN
##  7 drift  Victo… Takeawa… Test  -34.8    41.8   36.7  -11.4    11.9    NaN   NaN
##  8 drift  Weste… Takeawa… Test    4.78   12.8   11.4    2.31    6.92   NaN   NaN
##  9 mean   Austr… Takeawa… Test   13.3    13.4   13.3   54.8    54.8    NaN   NaN
## 10 mean   New S… Takeawa… Test  320.    323.   320.    58.9    58.9    NaN   NaN
## # ℹ 22 more rows
## # ℹ 1 more variable: ACF1 <dbl>

fc_benchmarks

## # A fable: 768 x 6 [1M]
## # Key:     State, Industry, .model [32]
##    State                        Industry               .model    Month
##    <chr>                        <chr>                  <chr>     <mth>
##  1 Australian Capital Territory Takeaway food services naive  2016 Jan
##  2 Australian Capital Territory Takeaway food services naive  2016 Feb
##  3 Australian Capital Territory Takeaway food services naive  2016 Mar
##  4 Australian Capital Territory Takeaway food services naive  2016 Apr
##  5 Australian Capital Territory Takeaway food services naive  2016 May
##  6 Australian Capital Territory Takeaway food services naive  2016 Jun
##  7 Australian Capital Territory Takeaway food services naive  2016 Jul
##  8 Australian Capital Territory Takeaway food services naive  2016 Aug
##  9 Australian Capital Territory Takeaway food services naive  2016 Sep
## 10 Australian Capital Territory Takeaway food services naive  2016 Oct
## # ℹ 758 more rows
## # ℹ 2 more variables: Turnover <dist>, .mean <dbl>

#The Drift method performs the best, as it yields the lowest RMSE.

5d.

Do the residuals from the best method resemble white noise?

# 5d.Answer:
Drift_model <- train_set %>%
  model(Drift = RW(Turnover ~ drift()))

residuals <- residuals(Drift_model) %>%
  drop_na(.resid)

residual_values <- as.numeric(residuals$.resid)

acf_residuals <- ACF(residuals, .resid, lag_max = 20)
autoplot(acf_residuals) 

hist(residual_values, breaks = 20)

qqnorm(residual_values)
qqline(residual_values)

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))
autoplot(myseries, Turnover)

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`))
is_tsibble(myseries)

## [1] TRUE

6a.

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

# 6a.Answer:
series <- unique(myseries$`Series ID`)
series

## [1] "A3349767W"

industry <- unique(myseries$Industry)
industry

## [1] "Clothing, footwear and personal accessory retailing"

6b.

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

# 6b.Answer:
Turnover_Trend_lr <- myseries %>%
  model(TSLM(Turnover ~ trend()))

6c.

See the regression result by report() command.

# 6c.Answer:
report(Turnover_Trend_lr)

## 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

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:
forecast_th <- Turnover_Trend_lr %>%
  forecast(h = "3 years")

forecast_th %>%
  as_tibble() %>%
  select(Month, .mean)

## # A tibble: 36 × 2
##       Month .mean
##       <mth> <dbl>
##  1 2019 Jan  14.9
##  2 2019 Feb  14.9
##  3 2019 Mar  15.0
##  4 2019 Apr  15.0
##  5 2019 May  15.0
##  6 2019 Jun  15.1
##  7 2019 Jul  15.1
##  8 2019 Aug  15.1
##  9 2019 Sep  15.2
## 10 2019 Oct  15.2
## # ℹ 26 more rows

forecast_th

## # 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>

6d.

Plot the forecast values along with the original data.

# 6d.Answer:
autoplot(myseries, Turnover) +
  autolayer(forecast_th, .mean, color = "blue")

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:
residuals <- augment(Turnover_Trend_lr)

gg_tsresiduals(Turnover_Trend_lr)

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:
library(ggplot2)

ggplot(jan_vic_elec, aes(x = Temperature, y = Demand)) +
  geom_point() 

demand_ml <- lm(Demand ~ Temperature, data = jan_vic_elec)

summary(demand_ml)

## 
## Call:
## lm(formula = Demand ~ Temperature, data = jan_vic_elec)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -49978 -10219   -121  18533  35441 
## 
## 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.89e-11

7b.

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

# 7b.Answer:
library(broom)

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

demand_ml <- lm(Demand ~ Temperature, data = jan_vic_elec)

residual_pl <- augment(demand_ml, data = jan_vic_elec)

ggplot(residual_pl, aes(x = Temperature, y = .resid)) +
  geom_point() 

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:
demand_ml <- lm(Demand ~ Temperature, data = jan_vic_elec)

next_day <- data.frame(Temperature = c(15, 35))

forecast <- predict(demand_ml, newdata = next_day)

forecast <- data.frame(
  Temperature = c(15, 35),
  Forecasted_Demand = forecast
)

forecast

##   Temperature Forecasted_Demand
## 1          15          151398.4
## 2          35          274484.2

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:
jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature)) %>%
  forecast(
    new_data(jan_vic_elec, 1) %>%
      mutate(Temperature = 15)
  ) %>%
  autoplot(jan_vic_elec)

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

forecast_3 <- jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature)) %>%
  forecast(
    new_data = new_data(jan_vic_elec, 1) %>%
      mutate(Temperature = 35)
  )

autoplot(jan_vic_elec, Demand) +
  autolayer(forecast_2, .mean) +
  autolayer(forecast_3, .mean) 

## Scale for colour_ramp is already present.
## Adding another scale for colour_ramp, which will replace the existing scale.

#I believe the forecasts show that higher temperatures drive up power demand, but other factors such as the day of the week, seasonality, and humidity also play a significant role in influencing power consumption.

7e.

Give prediction intervals for your forecasts.

# 7e.Answer:
forecast_2 <- jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature)) %>%
  forecast(
    new_data = new_data(jan_vic_elec, 1) %>%
      mutate(Temperature = 15),
    PI = TRUE
  )

forecast_3 <- jan_vic_elec %>%
  model(TSLM(Demand ~ Temperature)) %>%
  forecast(
    new_data = new_data(jan_vic_elec, 1) %>%
      mutate(Temperature = 35),
    PI = TRUE
  )

autoplot(jan_vic_elec, Demand) +
  autolayer(forecast_2, .mean,PI = TRUE) +
  autolayer(forecast_3, .mean,PI = TRUE)

## Warning in ggdist::geom_interval(intvl_mapping, data =
## dist_qi_frame(object[single_row[["TRUE"]], : Ignoring unknown parameters: `PI`

## Warning in ggplot2::geom_point(mapping = without(mapping, "linetype"), data =
## unpack_data(object[single_row[["TRUE"]], : Ignoring unknown parameters: `PI`

## Warning in ggdist::geom_interval(intvl_mapping, data =
## dist_qi_frame(object[single_row[["TRUE"]], : Ignoring unknown parameters: `PI`

## Warning in ggplot2::geom_point(mapping = without(mapping, "linetype"), data =
## unpack_data(object[single_row[["TRUE"]], : Ignoring unknown parameters: `PI`

## Scale for colour_ramp is already present.
## Adding another scale for colour_ramp, which will replace the existing scale.

#The prediction intervals represent the expected range of demand, where wider intervals indicate greater uncertainty and narrower intervals suggest higher confidence in the forecasts.

8.

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

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

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

shampoo_2 <- read_excel("shampoo-2.xlsx")
View(shampoo_2)
nrow(shampoo_2)

## [1] 36

ncol(shampoo_2)

## [1] 2

#There are 36 rows and 2 columns in the data
  

#b. Is the data annual, monthly, quarterly?
#After initial analysis, we can observe that the data is monthly
  
#c. Convert the data into tibble , then tsibble 
shampoo_data <- as_tibble(shampoo_2)

shampoo_tsibble <- shampoo_data %>%
  mutate(Date = yearmonth(Month)) %>%
  as_tsibble(index = Date)

shampoo_tsibble

## # A tsibble: 36 x 3 [1M]
##    Month               sales     Date
##    <dttm>              <dbl>    <mth>
##  1 1995-01-01 00:00:00  266  1995 Jan
##  2 1995-02-01 00:00:00  146. 1995 Feb
##  3 1995-03-01 00:00:00  183. 1995 Mar
##  4 1995-04-01 00:00:00  119. 1995 Apr
##  5 1995-05-01 00:00:00  180. 1995 May
##  6 1995-06-01 00:00:00  168. 1995 Jun
##  7 1995-07-01 00:00:00  232. 1995 Jul
##  8 1995-08-01 00:00:00  224. 1995 Aug
##  9 1995-09-01 00:00:00  193. 1995 Sep
## 10 1995-10-01 00:00:00  123. 1995 Oct
## # ℹ 26 more rows

#d. Plot the shampoo sales. What do you see from the data pattern? What does x-axis represent? 
plot(shampoo_tsibble$Date, shampoo_tsibble$sales, 
     xlab = "Date", 
     ylab = "Shampoo Sales", 
     main = "Shampoo Sales")

autoplot(shampoo_tsibble, sales) +
  labs(title = "Shampoo Sales",
       y = "Shampoo Sales",
       x = "Date")

# Comment here. Use plot() and autoplot().Put the name for y axis, and a title for the graph.
  
#e. What is the average, and median of shampoo sales. Put it on a histogram.
average <- mean(shampoo_tsibble$sales)

median <- median(shampoo_tsibble$sales)

hist(shampoo_tsibble$sales)

abline(v = average, col = "red", lwd = 2, lty = 2)
abline(v = median, col = "blue", lwd = 2, lty = 2)

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

gg_season(shampoo_tsibble, sales) +
  labs(title = "Seasonal Plot of Shampoo Sales",
       x = "Date",
       y = "Shampoo Sales")

#g. Get a linear regression line with trend and dummy for each month (Hint: use trend and season in regression equation).
shampoo_ml <- shampoo_tsibble %>%
  model(TSLM(sales ~ trend() + season()))
  
#h. Comment on each estimated coefficient of the model.Are they statistically significant at 5 % significance level?
report(shampoo_ml)

## Series: sales 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -129.60  -62.32   -4.84   53.76  152.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     113.867     55.740   2.043   0.0527 .  
## trend()          11.754      1.534   7.664 8.88e-08 ***
## season()year2   -33.154     73.630  -0.450   0.6567    
## season()year3   -53.808     73.678  -0.730   0.4726    
## season()year4   -24.628     73.757  -0.334   0.7415    
## season()year5   -56.015     73.869  -0.758   0.4560    
## season()year6   -27.802     74.012  -0.376   0.7106    
## season()year7     7.244     74.187   0.098   0.9231    
## season()year8   -37.043     74.393  -0.498   0.6233    
## season()year9    27.536     74.629   0.369   0.7155    
## season()year10  -32.518     74.897  -0.434   0.6682    
## season()year11    9.895     75.194   0.132   0.8964    
## season()year12   -4.259     75.522  -0.056   0.9555    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90.16 on 23 degrees of freedom
## Multiple R-squared: 0.7592,  Adjusted R-squared: 0.6336
## F-statistic: 6.043 on 12 and 23 DF, p-value: 0.00011612

#i. Which month has the highest sales?
Sales_Highest <- shampoo_tsibble %>%
  filter(sales == max(sales))

Sales_Highest

## # A tsibble: 1 x 3 [1M]
##   Month               sales     Date
##   <dttm>              <dbl>    <mth>
## 1 1997-09-01 00:00:00   682 1997 Sep

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

forecast <- shampoo_tsibble %>%
  model(TSLM(sales ~ trend() + season())) %>%
  forecast(h = "1 year")

forecast %>%
  as_tibble() %>%
  select(Date, .mean)

## # A tibble: 12 × 2
##        Date .mean
##       <mth> <dbl>
##  1 1998 Jan  549.
##  2 1998 Feb  527.
##  3 1998 Mar  518.
##  4 1998 Apr  559.
##  5 1998 May  540.
##  6 1998 Jun  580.
##  7 1998 Jul  627.
##  8 1998 Aug  594.
##  9 1998 Sep  670.
## 10 1998 Oct  622.
## 11 1998 Nov  676.
## 12 1998 Dec  674.

#k. Plot the forecast with original data.
autoplot(shampoo_tsibble, sales) +
  autolayer(forecast, .mean, colour = "blue") +
  labs(title = "Forecast for Next Year",
       x = "Date",
       y = "Shampoo Sales")

#l. Check if the residuals of the model is white noise.
shampoo_ml <- shampoo_tsibble %>%
  model(TSLM(sales ~ trend() + season()))

residuals <- augment(shampoo_ml)

gg_tsresiduals(shampoo_ml)

#m. By using the regression model, forecast the 1 year ahead, and then check the accuracy of the forecast. What is MSE, RMSE values?
forecast <- shampoo_tsibble %>%
  model(TSLM(sales ~ trend() + season())) %>%
  forecast(h = "1 year")

forecast %>%
  as_tibble() %>%
  select(Date, .mean)

## # A tibble: 12 × 2
##        Date .mean
##       <mth> <dbl>
##  1 1998 Jan  549.
##  2 1998 Feb  527.
##  3 1998 Mar  518.
##  4 1998 Apr  559.
##  5 1998 May  540.
##  6 1998 Jun  580.
##  7 1998 Jul  627.
##  8 1998 Aug  594.
##  9 1998 Sep  670.
## 10 1998 Oct  622.
## 11 1998 Nov  676.
## 12 1998 Dec  674.

forecast

## # A fable: 12 x 4 [1M]
## # Key:     .model [1]
##    .model                               Date
##    <chr>                               <mth>
##  1 TSLM(sales ~ trend() + season()) 1998 Jan
##  2 TSLM(sales ~ trend() + season()) 1998 Feb
##  3 TSLM(sales ~ trend() + season()) 1998 Mar
##  4 TSLM(sales ~ trend() + season()) 1998 Apr
##  5 TSLM(sales ~ trend() + season()) 1998 May
##  6 TSLM(sales ~ trend() + season()) 1998 Jun
##  7 TSLM(sales ~ trend() + season()) 1998 Jul
##  8 TSLM(sales ~ trend() + season()) 1998 Aug
##  9 TSLM(sales ~ trend() + season()) 1998 Sep
## 10 TSLM(sales ~ trend() + season()) 1998 Oct
## 11 TSLM(sales ~ trend() + season()) 1998 Nov
## 12 TSLM(sales ~ trend() + season()) 1998 Dec
## # ℹ 2 more variables: sales <dist>, .mean <dbl>