library(fpp3)## 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.4 ✔ fable 0.4.1
## ✔ ggplot2 3.5.1## Warning: package 'ggplot2' was built under R version 4.4.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)## ── 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 errorsSys.time()## [1] "2025-04-07 08:56:33 EDT"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:
library(ggplot2)
library(dplyr)
global_economy %>%
filter(Country == "Algeria") %>%
mutate(GDP_per_Capita = GDP / Population) %>%
ggplot(aes(x = Year, y = GDP_per_Capita)) +
geom_line(color = "steelblue", linewidth = 1.2) + # Updated from size to linewidth
geom_point(color = "darkred", size = 2) +
labs(
title = "Algeria's GDP per Capita Over Time",
x = "Year",
y = "GDP per Capita (USD)",
caption = "Source: Global Economy Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
# 1.Answer: This is GDP per capita for Australia. It's GDP annually since 1960. The GDP per capita has grown over the time. 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:
global_economy %>%
filter(Country == "Algeria") %>%
mutate(GDP_per_Capita = GDP / Population) %>%
ggplot(aes(x = Year, y = GDP_per_Capita)) +
geom_line(color = "darkgreen", linewidth = 1.2) + # Smooth line
geom_point(color = "orange", size = 2) + # Highlight points
labs(
title = "Algeria's GDP per Capita Over Time",
x = "Year",
y = "GDP per Capita (USD)",
caption = "Source: Global Economy Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
) 
2b.
United States GDP from global_economy.
# 2b.Answer:
global_economy %>%
filter(Country == "United States") %>%
ggplot(aes(x = Year, y = GDP)) +
geom_line(color = "blue", linewidth = 1.5) + # Bold line for better visibility
labs(
title = "United States GDP Over Time",
x = "Year",
y = "Gross Domestic Product (USD Billions)",
caption = "Source: Global Economy Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
# 2b.Answer: There has been a steady increase in GDP over the years.2c.
Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock
# 2c.Answer:
aus_livestock## # A tsibble: 29,364 x 4 [1M]
## # Key: Animal, State [54]
## Month Animal State Count
## <mth> <fct> <fct> <dbl>
## 1 1976 Jul Bulls, bullocks and steers Australian Capital Territory 2300
## 2 1976 Aug Bulls, bullocks and steers Australian Capital Territory 2100
## 3 1976 Sep Bulls, bullocks and steers Australian Capital Territory 2100
## 4 1976 Oct Bulls, bullocks and steers Australian Capital Territory 1900
## 5 1976 Nov Bulls, bullocks and steers Australian Capital Territory 2100
## 6 1976 Dec Bulls, bullocks and steers Australian Capital Territory 1800
## 7 1977 Jan Bulls, bullocks and steers Australian Capital Territory 1800
## 8 1977 Feb Bulls, bullocks and steers Australian Capital Territory 1900
## 9 1977 Mar Bulls, bullocks and steers Australian Capital Territory 2700
## 10 1977 Apr Bulls, bullocks and steers Australian Capital Territory 2300
## # ℹ 29,354 more rowsvictorian_cattle_slaughter <- aus_livestock %>%
filter(State == "Victoria", Animal == "Bulls, bullocks and steers")
ggplot(victorian_cattle_slaughter, aes(x = Month, y = Count)) +
geom_line() +
labs(title = "Slaughter of Victorian Bulls, Bullocks, and Steers",
x = "Month",
y = "Count") +
scale_x_yearmonth()
# 2c.Answer: We can observe an upward and downward trend, also seasonality between them.2d.
Victorian Electricity Demand from vic_elec.
# 2d.Answer:
vic_elec## # A tsibble: 52,608 x 5 [30m] <Australia/Melbourne>
## Time Demand Temperature Date Holiday
## <dttm> <dbl> <dbl> <date> <lgl>
## 1 2012-01-01 00:00:00 4383. 21.4 2012-01-01 TRUE
## 2 2012-01-01 00:30:00 4263. 21.0 2012-01-01 TRUE
## 3 2012-01-01 01:00:00 4049. 20.7 2012-01-01 TRUE
## 4 2012-01-01 01:30:00 3878. 20.6 2012-01-01 TRUE
## 5 2012-01-01 02:00:00 4036. 20.4 2012-01-01 TRUE
## 6 2012-01-01 02:30:00 3866. 20.2 2012-01-01 TRUE
## 7 2012-01-01 03:00:00 3694. 20.1 2012-01-01 TRUE
## 8 2012-01-01 03:30:00 3562. 19.6 2012-01-01 TRUE
## 9 2012-01-01 04:00:00 3433. 19.1 2012-01-01 TRUE
## 10 2012-01-01 04:30:00 3359. 19.0 2012-01-01 TRUE
## # ℹ 52,598 more rowsggplot(vic_elec, aes(x = Time, y = Demand)) +
geom_line() +
labs(title = "Victorian Electricity Demand Over Time",
x = "Time",
y = "Demand (GWh)")
# 2d.Answer: We can observe seasonality.2e.
Gas production from aus_production.
aus_production## # A tsibble: 218 x 7 [1Q]
## Quarter Beer Tobacco Bricks Cement Electricity Gas
## <qtr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1956 Q1 284 5225 189 465 3923 5
## 2 1956 Q2 213 5178 204 532 4436 6
## 3 1956 Q3 227 5297 208 561 4806 7
## 4 1956 Q4 308 5681 197 570 4418 6
## 5 1957 Q1 262 5577 187 529 4339 5
## 6 1957 Q2 228 5651 214 604 4811 7
## 7 1957 Q3 236 5317 227 603 5259 7
## 8 1957 Q4 320 6152 222 582 4735 6
## 9 1958 Q1 272 5758 199 554 4608 5
## 10 1958 Q2 233 5641 229 620 5196 7
## # ℹ 208 more rowsggplot(aus_production, aes(x = Quarter, y = Gas)) +
geom_line() +
labs(title = "Gas Production in Australia Over Time",
x = "Quarter",
y = "Gas Production (units)")
# 2e.Answer: There is an trend and seasonality over the years in the quarters.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:
canadian_gas## # A tsibble: 542 x 2 [1M]
## Month Volume
## <mth> <dbl>
## 1 1960 Jan 1.43
## 2 1960 Feb 1.31
## 3 1960 Mar 1.40
## 4 1960 Apr 1.17
## 5 1960 May 1.12
## 6 1960 Jun 1.01
## 7 1960 Jul 0.966
## 8 1960 Aug 0.977
## 9 1960 Sep 1.03
## 10 1960 Oct 1.25
## # ℹ 532 more rows# Time series plot of Canadian Gas Production
autoplot(canadian_gas, Volume) +
labs(
title = "Canadian Gas Production Over Time",
x = "Year",
y = "Gas Volume (Million Cubic Meters)",
caption = "Source: Canadian Gas Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
) 
# Subseries plot of Canadian Gas Production (Seasonality)
gg_subseries(canadian_gas, Volume) +
labs(
title = "Seasonal Subseries Plot of Canadian Gas Production",
x = "Month",
y = "Gas Volume (Million Cubic Meters)",
caption = "Source: Canadian Gas Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
gg_season(canadian_gas, Volume) +
labs(
title = "Seasonal Plot of Canadian Gas Production",
x = "Month",
y = "Gas Volume",
caption = "Source: Canadian Gas Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
# 3a.Answer: From the gg season plot we can see the monthly volume trends from 1969 to 1999, showing an overall upward trend. Volume generally peaks in January and December and dips mid-year, indicating a seasonal pattern. Earlier years (1969 to around 1979) show volume levels between 5 and 10, while later years (around 1989 and onward) reach volumes above 15, suggesting that whatever is being measured has grown significantly over time. During the year after 1989 there is much variation in gas production over the years with lowest production during Febraury and max production during March and December. Later years also show more pronounced fluctuations, suggesting increased volatility."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:
# Perform STL decomposition
gas_stl <- canadian_gas %>%
model(STL(Volume ~ season(window = "periodic"))) %>%
components()
# Plot the decomposition components
autoplot(gas_stl) +
labs(
title = "STL Decomposition of Canadian Gas Production",
x = "Year",
y = "Gas Volume (Million Cubic Meters)",
caption = "Source: Canadian Gas Dataset"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
3c.
How does the seasonal shape change over time? [Hint: Try plotting the seasonal component using gg_season().]
# 3c.Answer:
dcmp<-canadian_gas%>%
model(stl=STL(Volume))
components(dcmp)%>%
gg_season(Volume)
# 3c.Answer: During the years above 1999 the seasonality shape is constant for the years and is fluctuationg in the same manner for all with lowest in Febrarury and highest in December. Over the decades, volume levels increase, with recent years (1990s) reaching above 15, while earlier years remain between 5 and 10. 3d.
produce a plausible seasonally adjusted series? What are these numbers, plot the series.
# 3d.Answer:
dcmp<-canadian_gas%>%
model(stl=STL(Volume))
components(dcmp)## # A dable: 542 x 7 [1M]
## # Key: .model [1]
## # : Volume = trend + season_year + remainder
## .model Month Volume trend season_year remainder season_adjust
## <chr> <mth> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 stl 1960 Jan 1.43 1.08 0.520 -0.172 0.911
## 2 stl 1960 Feb 1.31 1.11 0.215 -0.0178 1.09
## 3 stl 1960 Mar 1.40 1.13 0.307 -0.0395 1.09
## 4 stl 1960 Apr 1.17 1.16 0.0161 -0.00627 1.15
## 5 stl 1960 May 1.12 1.18 -0.116 0.0476 1.23
## 6 stl 1960 Jun 1.01 1.21 -0.356 0.159 1.37
## 7 stl 1960 Jul 0.966 1.23 -0.403 0.136 1.37
## 8 stl 1960 Aug 0.977 1.26 -0.349 0.0677 1.33
## 9 stl 1960 Sep 1.03 1.28 -0.340 0.0870 1.37
## 10 stl 1960 Oct 1.25 1.31 -0.0899 0.0329 1.34
## # ℹ 532 more rowscanadian_gas%>%
autoplot(Volume, color='grey')+
autolayer(components(dcmp), season_adjust, color='red') +
labs(y="Volume in Billions Cubic Metres",
title="STL Decomposition of Canadian Gas Production with Seasonal Adjusted")
# 3d.Answer: By removing season year from the STL Model using Season Adjust in the Autolayer function, we can produce seasonal adjusted series."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?
# Answer:
#residuals plot
gg_tsresiduals(fit) ## 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()`).
# Answer: The ACF plot shows significant autocorrelation at several lags, indicating that residuals are not completely uncorrelated. The histogram of residuals suggests that not perfectly normal distribution, it shows slight skewness and a few outliers.")4e.
Produce forecasts for the test data with given code below:
# 4e Answer: Forecast for test data
# 4e Answer (Fixed)
#fc <- fit %>%
# forecast(new_data = anti_join(myseries, myseries_train, by = c("State", "Industry", "Series ID", "Month")))
#fc %>% autoplot(myseries)
library(fpp3)
set.seed(12345678)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
# Split into training and test sets
myseries_train <- myseries %>%
filter(year(Month) < 2011)
# Plot to check the split
autoplot(myseries, Turnover) +
autolayer(myseries_train, Turnover, colour = "red")
# Fit a seasonal naïve model
fit <- myseries_train %>%
model(SNAIVE(Turnover))
# Verify model
report(fit)## Series: Turnover
## Model: SNAIVE
##
## sigma^2: 1.2856# Prepare test data (exclude Turnover to avoid floating-point issues)
new_data <- anti_join(myseries, myseries_train,
by = c("State", "Industry", "Series ID", "Month"))
# Generate forecasts
fc <- fit %>% forecast(new_data)
# Plot forecast
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.768fc %>% 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 Answer: The error values are almost similar for both the sets. However, it works well for the training set than the testing set.4g.
How sensitive are the accuracy measures to the amount of training data used?
# 4g Answer:The slight differences in accuracy measures between the training and test datasets suggest that the model may exhibit some sensitivity to the amount of training data used.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:
myseries_train <- aus_retail %>%
filter(year(Month) < 2018)
autoplot(myseries, Turnover) + autolayer(myseries_train, Turnover, colour = "blue")
5b.
Fit all the appropriate benchmark methods to the training set and forecast the periods covered by the test set.
# 5b.Answer:
library(forecast)## Warning: package 'forecast' was built under R version 4.4.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoofit_auto_arima <- auto.arima(myseries_train$Turnover)
fit_ets <- ets(myseries_train$Turnover)
forecast_horizon <- 48 # Assuming forecasting for 4 years (48 months)
fc_auto_arima <- forecast(fit_auto_arima, h = forecast_horizon)
fc_ets <- forecast(fit_ets, h = forecast_horizon)5c.
Compute the accuracy of your forecasts. Which method does best?
accuracy_auto_arima <- accuracy(fc_auto_arima, actual = actual_values)
print(accuracy_auto_arima)## ME RMSE MAE MPE MAPE MASE
## Training set 0.005148139 50.67927 14.98015 -9.695522 18.73685 0.9403715
## ACF1
## Training set 0.002090398accuracy_ets <- accuracy(fc_ets, actual = actual_values)
print(accuracy_ets)## ME RMSE MAE MPE MAPE MASE
## Training set 0.002981578 53.32009 15.92701 -6.486425 16.30695 0.9998102
## ACF1
## Training set -0.2331982# 5c.Answer:The best method would be the Auto Arima method as it has a lower RMSE when compared to ETS method.5d.
Do the residuals from the best method resemble white noise?
residuals_auto_arima <- residuals(fit_auto_arima)
# Create a data frame for ggplot
residuals_df_auto_arima <- data.frame(Time = time(residuals_auto_arima),
Residuals = residuals_auto_arima)
# Plot residuals using ggplot2
ggplot(residuals_df_auto_arima, aes(x = Time, y = Residuals)) +
geom_line(color = "blue") +
labs(title = "Residuals of ARIMA Model",
x = "Time",
y = "Residuals") +
theme_minimal()## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# 5d.Answer: Yes, we can observe the white noise.6.
Using the code below, get a series (it gets a series randomly by using sample() function):
library(fpp3)
if (!exists("aus_retail")) {
stop("Error: 'aus_retail' dataset not found. Make sure fpp3 is installed and loaded.")
}
# Select a random time series
set.seed(12345678)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
# Ensuring no missing values
myseries <- myseries %>%
filter(!is.na(`Series ID`), !is.na(Turnover))
# Checkingdata
glimpse(myseries)## Rows: 369
## Columns: 5
## Key: State, Industry [1]
## $ State <chr> "Northern Territory", "Northern Territory", "Northern Terr…
## $ Industry <chr> "Clothing, footwear and personal accessory retailing", "Cl…
## $ `Series ID` <chr> "A3349767W", "A3349767W", "A3349767W", "A3349767W", "A3349…
## $ Month <mth> 1988 Apr, 1988 May, 1988 Jun, 1988 Jul, 1988 Aug, 1988 Sep…
## $ Turnover <dbl> 2.3, 2.9, 2.6, 2.8, 2.9, 3.0, 3.1, 3.0, 4.2, 2.7, 2.5, 2.4…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 3myseries = myseries %>% filter(!is.na(`Series ID`))6a.
What is the name of the series you randomly choose? Write it.
# 6a.Answer:
#The series randomly chosen is "Clothing, footwear, and personal accessory retailing" in the Northern Territory having series ID : A3349767W.6b.
Run a linear regression of Turnover on trend.(Hint: use TSLM() and trend() functions)
# 6b.Answer:
model <- myseries %>%
model(TSLM(Turnover ~ trend()))6c.
See the regression result by report() command.
# 6c.Answer:
report(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-166d.
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_values <- forecast(model, h = 12 * 3)
forecast_values## # 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) +
autolayer(forecast_values, series = "Forecasted Turnover") +
labs(title = "Forecasted Turnover for the Next 3 Years",
y = "Turnover") +
theme_minimal()## Plot variable not specified, automatically selected `.vars = Turnover`## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
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)
augment(fit)## # A tsibble: 273 x 8 [1M]
## # Key: State, Industry, .model [1]
## State Industry .model Month Turnover .fitted .resid .innov
## <chr> <chr> <chr> <mth> <dbl> <dbl> <dbl> <dbl>
## 1 Northern Territory Clothing, … SNAIV… 1988 Apr 2.3 NA NA NA
## 2 Northern Territory Clothing, … SNAIV… 1988 May 2.9 NA NA NA
## 3 Northern Territory Clothing, … SNAIV… 1988 Jun 2.6 NA NA NA
## 4 Northern Territory Clothing, … SNAIV… 1988 Jul 2.8 NA NA NA
## 5 Northern Territory Clothing, … SNAIV… 1988 Aug 2.9 NA NA NA
## 6 Northern Territory Clothing, … SNAIV… 1988 Sep 3 NA NA NA
## 7 Northern Territory Clothing, … SNAIV… 1988 Oct 3.1 NA NA NA
## 8 Northern Territory Clothing, … SNAIV… 1988 Nov 3 NA NA NA
## 9 Northern Territory Clothing, … SNAIV… 1988 Dec 4.2 NA NA NA
## 10 Northern Territory Clothing, … SNAIV… 1989 Jan 2.7 NA NA NA
## # ℹ 263 more rowsgg_tsresiduals(fit)## 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()`).
# 6e.Answer: The histogram looks somewhat normal but the time series plot shows some patterns and clusters, indicating non-randomness. Also the ACF plot reveals significant autocorrelation, which implies dependency in the residuals. So the residuals do not satisfy the white noise requirements.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?
library(dplyr)
ggplot(jan_vic_elec, aes(x = Temperature, y = Demand)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Electricity Demand vs. Temperature (January 2014)",
x = "Temperature (°C)",
y = "Electricity Demand")## `geom_smooth()` using formula = 'y ~ x'
# 7a.Answer: The positive relationship between temperature and electricity demand, can be observed by the statistically significant positive coefficient of the temperature variable in the regression model.7b.
Produce a residual plot. Is the model adequate? Are there any outliers or influential observations?
lm_model <- lm(Demand ~ Temperature, data = jan_vic_elec)
summary(lm_model)##
## 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-11residuals <- resid(lm_model)
plot(fitted(lm_model), residuals,
main = "Residual Plot",
xlab = "Fitted Values",
ylab = "Residuals")
abline(h = 0, col = "red") # Add horizontal line at y = 0
# 7b.Answer:There are a few outliers and so it can be concluded that the model is not adequate.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?
fc<-jan_vic_elec%>%
model(reg_model=TSLM(Demand~Temperature))
report(fc)## 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-11fc %>%
forecast(new_data(jan_vic_elec, 1) %>%
mutate(Temperature = 15)) %>%
autoplot(jan_vic_elec) + labs(title = "Next day temperature if max temp is 15∘C")
# 7c.Answer:Yes, it is believable as the temperature increases, demand also increased.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: Yes, it is believable as the temperature increases, the demand is also increasing. However, it might not seem well in certain scenarios. # 7d.Answer:
jan_vic_elec %>%
model(TSLM(Demand ~ Temperature)) %>%
forecast(
new_data(jan_vic_elec, 1) %>%
mutate(Temperature = 15)
) %>%
autoplot(jan_vic_elec)
7e.
Give prediction intervals for your forecasts.
# 7e.Answer:
forecast_15 <- predict(lm_model, newdata = data.frame(Temperature = 15), interval = "prediction", level = 0.95)
forecast_35 <- predict(lm_model, newdata = data.frame(Temperature = 35), interval = "prediction", level = 0.95)
forecast_15_point <- forecast_15[, 1]
forecast_15_interval <- forecast_15[, 2:3]
forecast_15_interval## lwr upr
## 97951.22 204845.488.
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.4.3shampoo_data <- read_excel("C:/Users/harsh/Downloads/shampoo-2.xlsx")
head(shampoo_data)## # A tibble: 6 × 2
## Month sales
## <dttm> <dbl>
## 1 1995-01-01 00:00:00 266
## 2 1995-02-01 00:00:00 146.
## 3 1995-03-01 00:00:00 183.
## 4 1995-04-01 00:00:00 119.
## 5 1995-05-01 00:00:00 180.
## 6 1995-06-01 00:00:00 168.num_rows <- nrow(shampoo_data)
num_rows## [1] 36num_variables <- ncol(shampoo_data)
num_variables## [1] 2#8a.Answer: It has a total of 36 Observations and 2 varaibles. It consists of 36 rows and 2 columns.
#b. Is the data annual, monthly, quarterly?
#8b.Answer: From the data, it appears to be monthly data as each row represents sales for a specific month.
#c. Convert the data into tibble , then tsibble
# Convert to tibble
library(tibble)
shampoo_tibble <- as_tibble(shampoo_data)
shampoo_tibble$Month <- as.Date(shampoo_tibble$Month)
shampoo_tsibble <- as_tsibble(shampoo_tibble, index = Month)
# Convert to tsibble
library(tsibble)
shampoo_tsibble <- shampoo_data %>%
mutate(Month = yearmonth(Month)) %>%
as_tsibble(index = Month)
#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.
plot(shampoo_data$sales,ylab = "Sales",main = "Shampoo Sales Over Time") 
#8d.Answer: From the data pattern, we can see some fluctuations in shampoo sales over time. We can observe seasonality with steady increase in shampoo sales over the months. The x-axis represents the months.
#e. What is the average, and median of shampoo sales. Put it on a histogram.
shampoo_tsibble %>%
autoplot(sales) +labs(y = "Sales",main = "Shampoo Sales Over Time") 
#f. Get seasonal plot. What do you see/ is there any pattern, is tehre any seasonality.
shampoo_tsibble%>%gg_season(sales)
#g. Get a linear regression line with trend and dummy for each month (Hint: use trend and season in regression equation).
shampoo_lm <- shampoo_tsibble %>%
model(TSLM(sales ~ trend() + season()))
report(shampoo_lm) ## 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#h. Comment on each estimated coefficient of the model.Are they statistically significant at 5 % significance level?
augment(shampoo_lm) %>%
ggplot(aes(x = Month)) +
geom_line(aes(y = sales, colour = "Data")) +
geom_line(aes(y = .fitted, colour = "Fitted")) +
labs(y = " Sales", title = "Monthly Sales") +
scale_colour_manual(values = c(Data = "black", Fitted = "red"))
#i. Which month has the highest sales?
shampoo_data%>%
group_by(Month = month(Month)) %>%
summarise(AvgSales = mean(sales)) %>%
arrange(desc(AvgSales))## # A tibble: 12 × 2
## Month AvgSales
## <dbl> <dbl>
## 1 11 394.
## 2 12 392.
## 3 9 388.
## 4 7 344.
## 5 10 340.
## 6 8 312.
## 7 6 298.
## 8 4 277.
## 9 1 267.
## 10 5 258.
## 11 2 245.
## 12 3 236.#j. Forecast it for the next year. What are the values
fcast_next_yr<-shampoo_lm%>%
forecast(h=12)
fcast_next_yr ## # A fable: 12 x 4 [1M]
## # Key: .model [1]
## .model Month
## <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>#k. Plot the forecast with original data.
fcast_next_yr%>%
autoplot(shampoo_tsibble) +
autolayer(fcast_next_yr, series = "Forecast") +
labs(title = "Shampoo Sales Forecast for Next Year", x = "Month", y = "Shampoo Sales")## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.
#l. Check if the residuals of the model is white noise.
gg_tsresiduals(shampoo_lm)
#8l.Answer: The residuals do not satisfy all requirements to be considered white noise, as they exhibit autocorrelation and some structure.
#m. By using the regression model, forecast the 1 year ahead, and then check the accuracy of the forecast. What is MSE, RMSE values?
if (!exists("shampoo_lm")) {
stop("Error: Model 'shampoo_lm' not found. Ensure you have trained the model first.")
}
# Generate a 1-year forecast
forecast_trends <- shampoo_lm %>% forecast(h = "1 year")
# Print the forecast
print(forecast_trends)## # A fable: 12 x 4 [1M]
## # Key: .model [1]
## .model Month
## <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>MSE_fc <- mean((shampoo_tsibble$sales - forecast_trends$.mean)^2)
MSE_fc ## [1] 98031.53RMSE_fc <- sqrt(mean((shampoo_tsibble$sales - forecast_trends$.mean)^2))
RMSE_fc## [1] 313.0999