Ekonometrika

Assignment 11

Preparation

library(tsibbledata)

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

library(tsibble)

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

library(DT)
library(ggplot2)

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

library(fpp2)

## Warning: package 'fpp2' was built under R version 4.0.4

## Warning: package 'forecast' was built under R version 4.0.4

## Warning: package 'fma' was built under R version 4.0.4

## Warning: package 'expsmooth' was built under R version 4.0.4

library(fpp3)

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

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

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

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

library(autoplotly)

## Warning: package 'autoplotly' was built under R version 4.0.4

library(dplyr)
library(forecast)
library(feasts)
library(lubridate)
library(expsmooth)
library(fma)
library(fabletools)
library(latex2exp)

## Warning: package 'latex2exp' was built under R version 4.0.5

library(seasonal)

## Warning: package 'seasonal' was built under R version 4.0.5

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 %>%
  tsibble(key = Code, index = Year)%>%
  autoplot(GDP/Population) +
  labs(title= "GDP per capita",
       y = "$US")

## Warning: Removed 3242 row(s) containing missing values (geom_path).

2 For each of the following series, make a graph of the data. If transforming seems appropriate, do so and describe the effect.

2.1 United States GDP from global_economy

global_economy %>%
  filter(Country == "United States") %>%
  autoplot(GDP)+
  labs(title = "United States GDP", y = "$US")+
  theme_replace()+
  geom_line(col = "#1B89D3")

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

aus_livestock %>%
  filter(Animal == "Bulls, bullocks and steers",
         State == "Victoria")%>%
  autoplot(Count)+
  labs(title = "Slaughter of Victorian “Bulls, bullocks and steers”")+
  theme_replace()+
  geom_line(col = "#1B89D3")

2.3 Victorian Electricity Demand from vic_elec

vic_elec %>%
  autoplot(Demand)+
  labs(title = "Victorian Electricity Demand",
       y = "Demand")+
  theme_replace()+
  geom_line(col = "#1B89D3")

2.4 Gas production from aus_production

aus_production %>%
  autoplot(Gas)+
  labs(title = "Gas production")+
  theme_replace()+
  geom_line(col = "#1B89D3")

3 Why is a Box-Cox transformation unhelpful for the canadian_gas data?

canadian_gas %>%
  autoplot(Volume)+
  labs(title = "Canadian Gas Production",
       y = "Monthly Canadian Gas Production (billions of cubic meter)")+
  theme_replace()+
  geom_line(col = "#1B89D3")

lambda_cangas <- canadian_gas %>%
                  features(Volume, features = guerrero) %>%
                  pull(lambda_guerrero)
canadian_gas %>%
  autoplot(box_cox(Volume, lambda = lambda_cangas))+
  labs(title = latex2exp::TeX(paste0(
         "Box Cox Transformation of Canadian Gas Production with $\\lambda$ = ",
         round(lambda_cangas,2))))+
  theme_replace()+
  geom_line(col = "#1B89D3")

4 What Box-Cox transformation would you select for your retail data (from Exercise 8 in Section 2.10)?

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

myseries %>%
  autoplot(Turnover)+
  labs(title = "Australian Retail Trade Turnover")+
  theme_replace()+
  geom_line(col = "#1B89D3")

lambda_retail <- myseries %>%
                  features(Turnover, features = guerrero) %>%
                  pull(lambda_guerrero)

myseries %>%
  autoplot(box_cox(Turnover, lambda_retail))+
  labs(title = latex2exp::TeX(paste0(
         "Box Cox Transformation of Australian Retail Trade Turnover with $\\lambda$ = ",
         round(lambda_retail,2))))+
  theme_replace()+
  geom_line(col = "#1B89D3")

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

5.1 Tobacco from aus_production

lambda_tobacco <- aus_production %>%
                   features(Tobacco, features = guerrero) %>%
                   pull(lambda_guerrero)
aus_production %>%
  autoplot(box_cox(Tobacco, lambda_tobacco)) +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed gas production with $\\lambda$ = ",
         round(lambda_tobacco,2))))+
  theme_replace()+
  geom_line(col = "#1B89D3")

## Warning: Removed 24 row(s) containing missing values (geom_path).

## Warning: Removed 24 row(s) containing missing values (geom_path).

5.2 Economy class passengers between Melbourne and Sydney from ansett

lambda_class <- ansett %>%
                 filter(Class == "Economy",
                        Airports == "MEL-SYD")%>%
                 features(Passengers, features = guerrero) %>%
                 pull(lambda_guerrero)
ansett %>%
  filter(Class == "Economy",
         Airports == "MEL-SYD")%>%
  mutate(Passengers = Passengers/1000) %>%
  autoplot(box_cox(Passengers, lambda = lambda_class)) +
  labs(y = "Passengers ('000)",
       title = latex2exp::TeX(paste0(
         "Transformed Ansett Airlines Economy Class with $\\lambda$ = ",
         round(lambda_class,2))),
       subtitle = "Melbourne-Sydney")+
  theme_replace()+
  geom_line(col = "#1B89D3")

5.3 Pedestrian counts at Southern Cross Station from pedestrian

lambda_count <- pedestrian %>%
                filter(Sensor == "Southern Cross Station") %>%
                 features(Count, features = guerrero) %>%
                 pull(lambda_guerrero)
pedestrian %>%
  filter(Sensor == "Southern Cross Station") %>%
  autoplot(box_cox(Count,lambda_count))+
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed Pedestrian Counts at Southern Cross Station with $\\lambda$ = ",
         round(lambda_count,2))))+
  theme_replace()+
  geom_line(col = "#1B89D3")

6 Show that a 3×5 MA is equivalent to a 7-term weighted moving average with weights of 0.067, 0.133, 0.200, 0.200, 0.200, 0.133, and 0.067

\[ \begin{align} 3 × 5MA= {\ 1 \over 15}Y_1 + {\ 2 \over 15}Y_2 + {\ 3 \over 15}Y_3 + {\ 3 \over 15}Y_4 + {\ 3 \over 15}Y_5 + {\ 2 \over 15}Y_6 + {\ 1 \over 15}Y_7 \\ \end{align} \]

7 Consider the last five years of the Gas data from aus_production

gas <- tail(aus_production, 5*4) %>% select(Gas)
datatable(gas)

7.1 Plot the time series. Can you identify seasonal fluctuations and/or a trend-cycle?

gas %>%
  autoplot()+
  labs(title = "Last Five Years of The Gas Data")+
  theme_replace()+
  geom_line(col = "#1B89D3")

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

7.2 Use classical_decomposition with type=multiplicative to calculate the trend-cycle and seasonal indices

gas %>%
  model(classical_decomposition(Gas,type = "multiplicative")) %>%
  components() %>%
  autoplot() + 
  ggtitle("Last Five Years of The Gas Data")

## Warning: Removed 2 row(s) containing missing values (geom_path).

Hasil dekomposisi multiplikatif menunjukkan komponen musiman triwulanan dengan frekuensi 1 tahun. Ada tren peningkatan dari tahun 2006 hingga pertengahan 2007. Setelah tahun 2007, tidak ada tren hingga awal 2008. Setelah itu, tidak ada tren yang meningkat akhir 2009.

7.3 Do the results support the graphical interpretation from part a?

Hasilnya mendukung interpretasi grafis dari bagian 7.1, yang merupakan musiman frekuensi 1 tahun dan tren yang meningkat. Karena dekomposisi multiplikatif klasik bergantung pada rata-rata bergerak, tidak ada data di awal dan akhir siklus tren.

7.4 Compute and plot the seasonally adjusted data

gas_decom <- gas %>%
             model(classical_decomposition(Gas,type = "multiplicative")) %>%
             components()
gas_decom %>%
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "",
       title = "Last Five Years of The Gas Data") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

## Warning: Removed 4 row(s) containing missing values (geom_path).

7.5 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?

new_gas <- gas
new_gas$Gas[10] <- new_gas$Gas[10]+300

new_gas %>%
  model(classical_decomposition(Gas,type = "multiplicative")) %>%
  components() %>%
  autoplot() + 
  ggtitle("Last Five Years of The Gas Data with 300 added to the 10th observation")

## Warning: Removed 2 row(s) containing missing values (geom_path).

new_gas %>%
  model(classical_decomposition(Gas,type = "multiplicative")) %>%
  components() %>%
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "",
       title = "Last Five Years of The Gas Data with 300 added to the 10th observation") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

## Warning: Removed 4 row(s) containing missing values (geom_path).

7.6 Does it make any difference if the outlier is near the end rather than in the middle of the time series?

new_gas2 <- gas
new_gas2$Gas[20] <- new_gas2$Gas[10]+300

new_gas2 %>%
  model(classical_decomposition(Gas,type = "multiplicative")) %>%
  components() %>%
  autoplot() + 
  ggtitle("Last Five Years of The Gas Data with 300 added to the 10th observation")

## Warning: Removed 2 row(s) containing missing values (geom_path).

new_gas2 %>%
  model(classical_decomposition(Gas,type = "multiplicative")) %>%
  components() %>%
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "",
       title = "Last Five Years of The Gas Data with 300 added to the last observation") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

## Warning: Removed 4 row(s) containing missing values (geom_path).

8 Recall your retail time series data (from Exercise 8 in Section 2.10). Decompose the series using X-11. Does it reveal any outliers, or unusual features that you had not noticed previously?

myseries %>%
  model(classical_decomposition(Turnover,type = "multiplicative")) %>%
  components() %>%
  autoplot() + 
  ggtitle("Multiplicative decomposition of my retail time series data")

## Warning: Removed 6 row(s) containing missing values (geom_path).

myseries %>%
  model(x11 = X_13ARIMA_SEATS(Turnover ~ x11())) %>%
  components() %>%
  autoplot()+
  labs(title = "X-11 decomposition of my retail time series data")

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

9.1 Write about 3–5 sentences describing the results of the decomposition. Pay particular attention to the scales of the graphs in making your interpretation

Mengisolasi komponen tren dari komponen musiman menunjukkan bahwa tren telah meningkat melalui sebagian besar kerangka waktu, dengan beberapa periode stasioner terjadi pada awal 90-an. Rincian bulanan komponen musiman menunjukkan bahwa beberapa bulan menunjukkan velo kota yang lebih besar dalam variasi mereka daripada bulan-bulan lainnya.

9.2 Is the recession of 1991/1992 visible in the estimated components?

Yes, kami melihat penurunan pekerjaan selama 1991/1992 yang tidak dijelaskan oleh musim atau tren positif.

10 This exercise uses the canadian_gas data (monthly Canadian gas production in billions of cubic metres, January 1960 – February 2005)

10.1 Plot the data using autoplot(), gg_subseries() and gg_season() to look at the effect of the changing seasonality over time

canadian_gas %>%
  autoplot(Volume)+
  labs(title = "Monthly Canadian Gas Production",
       subtitle = "autoplot()",
       y = "billions of cubic meter")+
  theme_replace()+
  geom_line(col = "#1B89D3")

canadian_gas %>%
  gg_subseries(Volume)+
  labs(title = "Monthly Canadian Gas Production",
       subtitle = "gg_subseries()",
       y = "billions of cubic meter")

canadian_gas %>%
  gg_season(Volume)+
  labs(title = "Monthly Canadian Gas Production",
       subtitle = "gg_season()",
       y = "billions of cubic meter")

Data canadian_gas memiliki tren peningkatan yang jelas dan musim yang kuat.Pengamatan ini diverifikasi di plot subseries dan plot musiman. Secara umum, produksi gas menurun di musim panas dan meningkat di musim dingin. Musiman meningkat secara dramatis dari 1975 hingga 1990.

10.2 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

canadian_gas %>%
  model(
    STL(Volume ~ trend(window = 21) +
                   season(window = 13),
    robust = TRUE)) %>%
  components() %>%
  autoplot()+
  labs(title = "STL decomposition of Canadian Gas Production")

Hasil penguraian STL ditunjukkan di atas. Komponen tren cukup mewakili tren dalam data asli. Komponen musiman meningkat dari 1975 hingga 1985, dan kemudian menurun. Pengamatan ini konsisten dengan plot waktu. Selain itu, komponen sisanya sekitar nol.

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

Seperti yang ditunjukkan di atas, bentuk musiman datar dari awal dan kemudian seiring berjalannya waktu bentuk musiman meningkat. Pada tahun 1960 tidak ada siklus tren, kita dapat mengatakan produksi gas tidak benar-benar tren pada waktu itu. Setelah tahun 1975 ada siklus tren, maka produksi gas meningkat pada saat itu dan sebagainya.

10.4 Can you produce a plausible seasonally adjusted series?

canadian_gas %>%
 model(
    STL(Volume ~ trend(window = 21) +
                   season(window = 13),
    robust = TRUE)) %>%
  components() %>%
  ggplot(aes(x = Month)) +
  geom_line(aes(y = Volume, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(title = "STL decomposition of Canadian Gas Production") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

10.5 Compare the results with those obtained using SEATS and X-11. How are they different?

canadian_gas %>%
  model(x11 = X_13ARIMA_SEATS(Volume ~ x11())) %>%
  components() %>%
  autoplot()+
  labs(title = "X-11 decomposition of Canadian Gas Production")