Data 624 HW-2
Biyag Dukuray
2024-09-11
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.0 ──## ✔ tibble 3.2.1 ✔ tsibble 1.1.5
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.1 ✔ feasts 0.3.2
## ✔ lubridate 1.9.3 ✔ fable 0.3.4
## ✔ ggplot2 3.5.1 ✔ fabletools 0.4.2## Warning: package 'tsibble' 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)## ── 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 errorslibrary(tsibble)
library(tsibbledata)
library(latex2exp)
library(seasonal)## Warning: package 'seasonal' was built under R version 4.3.3## Warning in system(paste(x13.bin, file.path(tdir, "Testairline")), intern =
## TRUE, : running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x13binary__14f225377e4e__dir/Testairline
## 2>/dev/null' had status 134## The binaries provided by 'x13binary' do not work on this
## machine. To get more information, run:
## x13binary::checkX13binary()
##
## You can set 'X13_PATH' manually if you intend to use your own
## binaries. See ?seasonal for details.
##
##
## Attaching package: 'seasonal'
##
## The following object is masked from 'package:tibble':
##
## viewDATA 624 HW-2
#1 Consider the GDP information in global_economy. Plot the GDP per capita for each country over time. Which country has the highest GDP per capita? How has this changed over time?
global_economy %>% mutate(GDP_per_capita = GDP / Population) %>%
ggplot(aes(x = Year, y = GDP_per_capita, colour = Country)) +
geom_line(show.legend = FALSE) +
labs(title = "Country GDP Per Capita Over Time",
x = "Year", y = "GDP Per Capita") + theme_minimal() ## Warning: Removed 3242 rows containing missing values or values outside the scale range
## (`geom_line()`).
Monaco has the highest GPD per capita and this has been consistent over time. It has been in the top of the category for many years including Liechtenstein. Its capita per GDP has continued to follow an uptrend as far as the data displays. These two places have had the highest GPD per capita for many years in a row.
# Countries with the highest GDP per capita
global_economy %>%
mutate(GDP_per_capita = GDP / Population) %>%
group_by(Country) %>%
summarise(Max_GDP_per_capita = max(GDP_per_capita, na.rm = TRUE)) %>%
arrange(desc(Max_GDP_per_capita))## Warning: There were 3325 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `Max_GDP_per_capita = max(GDP_per_capita, na.rm = TRUE)`.
## ℹ In group 23: `Country = "Afghanistan"` and `Year = 1982`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3324 remaining warnings.## Warning: Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `Country`, `Year` first.## # A tsibble: 15,150 x 3 [1Y]
## # Key: Country [263]
## Country Year Max_GDP_per_capita
## <fct> <dbl> <dbl>
## 1 Monaco 2014 185153.
## 2 Monaco 2008 180640.
## 3 Liechtenstein 2014 179308.
## 4 Liechtenstein 2013 173528.
## 5 Monaco 2013 172589.
## 6 Monaco 2016 168011.
## 7 Liechtenstein 2015 167591.
## 8 Monaco 2007 167125.
## 9 Liechtenstein 2016 164993.
## 10 Monaco 2015 163369.
## # ℹ 15,140 more rows# Country with the highest GDP per capita
global_economy %>%
mutate(GDP_per_capita = GDP / Population) %>%
filter(GDP_per_capita == max(GDP_per_capita, na.rm = TRUE)) %>%
select(Country, Year, GDP_per_capita)## # A tsibble: 1 x 3 [1Y]
## # Key: Country [1]
## Country Year GDP_per_capita
## <fct> <dbl> <dbl>
## 1 Monaco 2014 185153.#2 For each of the following series, make a graph of the data. If transforming seems appropriate, do so and describe the effect.
United States GDP from global_economy.
# United States GDP
global_economy %>%
filter(Country == "United States") %>%
autoplot(GDP) +
labs(title = "United States GDP Over Time")
Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock.
# Victorian livestock slaughter
aus_livestock %>%
filter(State == "Victoria", Animal == "Bulls, bullocks and steers") %>%
autoplot(Count) +
labs(title = "Victorian Livestock Slaughter ('Bulls, bullocks and steers')")
Victorian Electricity Demand from vic_elec. For this plot I transformed the data to show daily points instead of half-hour intervals. This made the graph alot cleaner so it is easier to spot trends and cycles within the data. From this we can clearly see that electricity demand is its highest in the beginning and the middle of the year which are the coldest and hottest periods of the year, respectively.
# Victorian Electricity Demand
autoplot(vic_elec, Demand) +
labs(title = "Victorian Electricity Demand")
# Aggregate half-hourly data to daily data and plot
vic_elec %>%
index_by(Date = as_date(Time)) %>%
summarise(Daily_Demand = sum(Demand)) %>%
autoplot(Daily_Demand) +
labs(title = "Victorian Electricity Daily Demand",
x = "Date",
y = "Daily Electricity Demand (MW)") +
theme_minimal()
Gas production from aus_production.
# Gas production
autoplot(aus_production, Gas) +
labs(title = "Gas Production Over Time")
#3 Why is a Box-Cox transformation unhelpful for the canadian_gas data?
canadian_gas %>%
autoplot(Volume) +
labs(title= "Monthly Canadian Gas Production", y = "B Cubic meters")
lambda <- canadian_gas %>%
features(Volume, features = guerrero) %>%
pull(lambda_guerrero)
canadian_gas %>%
autoplot(box_cox(Volume, lambda)) +
labs(y = "", title = TeX(paste0("Transformed Gas Production with $\\lambda$ = ",
round(lambda,2))))
The Box-Cox transformation is unhelpful for the canadian_gas data because the series already has a stable variance. Thus, the seasonal variation is consistent throughout the majority of the series. The transformation is primarily useful for stabilizing variance in data with non-constant variance. Applying it to this series would have no real impact.
#4 What Box-Cox transformation would you select for your retail data (from Exercise 7 in Section 2.10)?
lambda <- aus_retail %>%
filter(State == "Victoria") %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)aus_retail %>% filter(State == “Victoria”) %>% autoplot(box_cox(Turnover, lambda)) + labs(title = “Retail Data with Box-Cox Transformation”, y = “Transformed Turnover”)
#5
For the following series, find an appropriate Box-Cox transformation in order to stabilise the variance. Tobacco from aus_production, Economy class passengers between Melbourne and Sydney from ansett, and Pedestrian counts at Southern Cross Station from pedestrian.
```r
# Box-Cox transformation for Tobacco production
autoplot(aus_production, Tobacco)+
labs(title = "Tobacco Production")## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
lambda_tobacco <- aus_production %>%
features(Tobacco, features = guerrero) %>%
pull(lambda_guerrero)
aus_production %>%
autoplot(box_cox(Tobacco, lambda_tobacco)) +
labs(title = "Tobacco Production with Box-Cox Transformation", x = "Year", y = "Tobacco Production (Transformed)")## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
# Box-Cox transformation for economy class passengers
autoplot(ansett %>% filter(Class == "Economy", Airports == "MEL-SYD")) + labs(title = "Economy Passengers Between Melbourne and Sydney")## Plot variable not specified, automatically selected `.vars = Passengers`
lambda_passengers <- ansett %>%
filter(Airports == "MEL-SYD" & Class == "Economy") %>%
features(Passengers, features = guerrero) %>%
pull(lambda_guerrero)
ansett %>%
filter(Airports == "MEL-SYD" & Class == "Economy") %>% autoplot(box_cox(Passengers, lambda_passengers)) + labs(title = "Economy Passengers (Melbourne-Sydney) with Box-Cox Transformation", x = "Year", y = "Passengers (Transformed)")
# Pedestrian counts at Southern Cross Station
autoplot(pedestrian %>% filter(Sensor == 'Southern Cross Station')) + labs(title = "Pedestrian counts at Southern Cross Station")## Plot variable not specified, automatically selected `.vars = Time`
lambda <- pedestrian %>%
filter(Sensor == 'Southern Cross Station') %>%
features(Count, features = guerrero) %>%
pull(lambda_guerrero)
pedestrian %>%
filter(Sensor == 'Southern Cross Station') %>%
autoplot(box_cox(Count, lambda)) +
labs(y = "", title = "Transformed pedestrian count at Southern Cross Station")
#7 Consider the last five years of the Gas data from aus_production.
- Plot the time series. Can you identify seasonal fluctuations and/or a trend-cycle?
Seasonal fluctuations and a trend cycle are both present in this data set. The seasonal fluctuations consistently displayed over several years consist of lows in Q1 and highs in Q3. The trend cycle shows an increase over the past 5 years.
gas <- tail(aus_production, 5*4) %>% select(Gas)
# Plot time series
autoplot(gas) +
labs(title = "Gas Production (Last 5 Years)", y = "Gas Production")## Plot variable not specified, automatically selected `.vars = Gas`
- Use classical_decomposition with type=multiplicative to calculate the trend-cycle and seasonal indices.
decomp <- gas %>%
model(classical_decomposition(Gas, type = "multiplicative"))
components(decomp) %>%
autoplot() +
labs(title = "Classical Decomposition of Gas Data")## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
- Do the results support the graphical interpretation from part a?
Yes the results support the graph in part a. The graph shows an uptrend through the 5 years which is displayed in the results. The graph also shows a seasonal pattern which is also displayed in the results. The randomness component of the results is also consistent with the graph in part a.
- Compute and plot the seasonally adjusted data.
components(decomp) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Gas, colour = "Data")) +
geom_line(aes(y = season_adjust, colour = "Seasonally Adjusted")) +
labs(
title = "Seasonally Adjusted Gas Production",
y = "Petajoules"
) +
scale_colour_manual(
values = c("gray", "#1f78b4"),
breaks = c("Data", "Seasonally Adjusted")
)
e. Change one observation to be an outlier (e.g., add 300 to one
observation), and recompute the seasonally adjusted data. What is the
effect of the outlier?
gas_outlier_begin <- gas
gas_outlier_begin$Gas[1] <- gas_outlier_begin$Gas[1] + 300
decomp_outlier_begin <- gas_outlier_begin %>%
model(classical_decomposition(Gas, type = "multiplicative"))
components(decomp_outlier_begin) %>%
autoplot() +
labs(title = "Decomposition with Outlier at the Beginning")## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
gas_outlier_middle <- gas
gas_outlier_middle$Gas[10] <- gas_outlier_middle$Gas[10] + 300
decomp_outlier_middle <- gas_outlier_middle %>%
model(classical_decomposition(Gas, type = "multiplicative"))
components(decomp_outlier_middle) %>%
autoplot() +
labs(title = "Decomposition with Outlier in the Middle")## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
gas_outlier_end <- gas
gas_outlier_end$Gas[20] <- gas_outlier_end$Gas[20] + 300
decomp_outlier_end <- gas_outlier_end %>%
model(classical_decomposition(Gas, type = "multiplicative"))
components(decomp_outlier_end) %>%
autoplot() +
labs(title = "Decomposition with Outlier at the End")## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
- Does it make any difference if the outlier is near the end rather than in the middle of the time series?
Yes, the position of the outlier does make a difference in the plot results. The outlier at the beginning shows an anomaly at the beginning of the plots. The outlier in the middle shows an anomaly in the middle of the plots. The outlier in the end shows an anomaly at the end of the plots. These indicate that the where the outlier is in the time series does make a difference.
#8 Recall your retail time series data (from Exercise 7 in Section 2.10). Decompose the series using X-11. Does it reveal any outliers, or unusual features that you had not noticed previously?
x11_decomposition <- aus_retail %>%
filter(State == "Victoria") %>%
model(x11_model = X_13ARIMA_SEATS(Turnover ~ x11()))## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f23a3765c5/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f2afa865f/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f2458bcc9a/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f25d25b621/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f256d8604d/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f2179a8bc6/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f2234e8045/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f276c7441e/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f214494a7e/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f257b9a099/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f25bdec9cd/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f27f92e0d9/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f23e807fe/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f27664b6b2/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f24e169cd7/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f22e930f4e/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f238d1e1c5/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f25b3d69a7/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f21ce381bd/iofile
## -n -s 2>/dev/null' had status 134## Warning in system(paste(x13.bin, m_flag, file, flags), intern = TRUE,
## ignore.stderr = TRUE): running command
## '/Users/zigcah/Library/R/x86_64/4.3/library/x13binary/bin/x13ashtml
## /var/folders/_x/bvfjszz17kn5q7qwy09v6gwc0000gn/T//RtmpNbeMwu/x1314f21c5aae1c/iofile
## -n -s 2>/dev/null' had status 134## Warning: 20 errors (1 unique) encountered for x11_model
## [20] X-13 has returned a non-zero exist status, which means that the current spec file cannot be processed for an unknown reason.{r # Plot the decomposition autoplot(x11_decomposition) + labs( title = "X-11 Decomposition of Turnover in Queensland Takeaway Food Services", y = "Turnover", x = "Time" )
#9 Figures 3.19 and 3.20 show the result of decomposing the number of persons in the civilian labour force in Australia each month from February 1978 to August 1995.
Write about 3–5 sentences describing the results of the decomposition. Pay particular attention to the scales of the graphs in making your interpretation. Is the recession of 1991/1992 visible in the estimated components?
The decomposition results show distinct components of the labor force data. The seasonal component demonstrates periodic patterns, while the trend component reveals an overall growth trajectory with occasional dips, such as during the recession of 1991/1992. The scales of the graphs highlight the relative contributions of trend, seasonal, and irregular components to the overall variation. The recession is visible as a downturn in the trend component.