All libraries needed for the Homework

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 --
## v tibble      3.2.1     v tsibble     1.1.5
## v dplyr       1.1.2     v tsibbledata 0.4.1
## v tidyr       1.3.0     v feasts      0.3.2
## v lubridate   1.9.2     v fable       0.3.4
## v ggplot2     3.5.1     v fabletools  0.4.2
## Warning: package 'ggplot2' was built under R version 4.3.3
## Warning: package 'tsibble' was built under R version 4.3.3
## Warning: package 'tsibbledata' was built under R version 4.3.3
## Warning: package 'feasts' was built under R version 4.3.3
## Warning: package 'fabletools' was built under R version 4.3.3
## Warning: package 'fable' was built under R version 4.3.3
## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date()    masks base::date()
## x dplyr::filter()      masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval()  masks lubridate::interval()
## x dplyr::lag()         masks stats::lag()
## x tsibble::setdiff()   masks base::setdiff()
## x tsibble::union()     masks base::union()
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tidyverse)
## -- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
## v forcats 1.0.0     v readr   2.1.4
## v purrr   1.0.1     v stringr 1.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter()     masks stats::filter()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag()        masks stats::lag()
## i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(lubridate)
library(tsibble)
library(pracma)
## Warning: package 'pracma' was built under R version 4.3.3
## 
## Attaching package: 'pracma'
## 
## The following object is masked from 'package:purrr':
## 
##     cross
library(mlbench)
## Warning: package 'mlbench' was built under R version 4.3.3
library(corrplot)
## corrplot 0.92 loaded
library(e1071)
## Warning: package 'e1071' was built under R version 4.3.1
## 
## Attaching package: 'e1071'
## 
## The following object is masked from 'package:pracma':
## 
##     sigmoid
## 
## The following object is masked from 'package:fabletools':
## 
##     interpolate
library(psych)
## Warning: package 'psych' was built under R version 4.3.1
## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:pracma':
## 
##     logit, polar
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(imputeTS)
## Warning: package 'imputeTS' was built under R version 4.3.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.3.1
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

9.1 - Figure 9.32 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

  1. Explain the differences among these figures. Do they all indicate that the data are white noise?

The differences among the figures are the height of the spike levels and the area between the dashed lines. It looks like as the time series increases, the spike levels seem to decrease as well as the areas between the dashed lines seems to decrease. The data between the two dashed lines is white noise.

  1. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

The critical values are at different distances from the mean to zero as the are give by the formula: ± 1.96/√T,T denoting the length of the time series. As T increases, the critical value decrease and get closer to zero. The autocorrelations are different in each figure when they refer to white noise due to larger series sizes with random numbers, they decrease the chance of autocorrelation.

9.2- A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

gafa_stock %>%
  filter(Symbol == "AMZN") %>%
  gg_tsdisplay(Close, plot_type='partial') +
  labs(title = "Amazon Closing Stock Price")
## Warning: Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.
## Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.

The time series has an increasing trend,making it not stationary. The ACF plot shows a slow decrease and has a r1 which is positive,however, a stationary time series will basically end up being zero. Since, the PACF plot has a lag of 1, this indicates that the series is not stationary. In consequence, the mean needs to be stationarized and should be differenced. It is evident that the closing stock price is related to the closing price (y{t-1}) but is not really affected by prices of previous days.

9.3- For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

(a).Turkish GDP from global_economy (b).Accommodation takings in the state of Tasmania from aus_accommodation. (c).Monthly sales from souvenirs

(a).Turkish GDP from global_economy

#Turkey GDP
turkey <- global_economy %>%
  filter(Country == "Turkey")
#plot Turkey GDP
turkey_plt <- turkey %>% 
  autoplot(GDP) +
  labs(title="Turkey GDP")

#ACF plot of Turkey GDP
turkey_ACF <- turkey %>%
  ACF(GDP) %>%
  autoplot() + labs(title="Correlation of Turkey GDP ")

#PCF plot of Turkey GDP
turkey_PCF <- turkey %>%
  PACF(GDP) %>%
  autoplot() + labs(title="Partial Autocorrelation of Turkey GDP ")

grid.arrange(turkey_plt, turkey_ACF,turkey_PCF , nrow = 2)

#utilize guerrero feature to observe transformation of 0.16
lambda_Turkey <- turkey %>%
  features(GDP,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_Turkey
## [1] 0.1572187
################## Turkey GDP after Transformation ##################

#Turkey GDP after transformation
turkey_plot_box_cox <- turkey %>% 
  autoplot(difference(box_cox(GDP, lambda_Turkey), 1)) +
  labs(title="Turkey GDP After Box Cox Transformation")

turkey_ACF_box_cox<- turkey %>%
  ACF(difference(box_cox(GDP, lambda_Turkey), 1)) %>%
  autoplot() + labs(title="Correlation of Turkey GDP After Box Cox Transformation ")

turkey_PCF_box_cox <- turkey %>%
  PACF(difference(box_cox(GDP, lambda_Turkey), 1)) %>%
  autoplot() + labs(title="Partial Autocorrelation of Turkey GDP After Box Cox Transformation ")

grid.arrange(turkey_plot_box_cox, turkey_ACF_box_cox, turkey_PCF_box_cox, nrow = 2)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

#funciton which shows whether additional differencing is needed -  not needed in this case
turkey %>%
  features(box_cox(GDP, lambda_Turkey), unitroot_ndiffs)
## # A tibble: 1 x 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1

It is clear from the original plots that the Turkish GDP data is not stationary. A Box-Cox transformation with lambda 0.16 and order differencing of 1 is applied.

(b). Accommodation takings in the state of Tasmania from aus_accommodation

################## Tasmania Accommodation  before Transformation ##################
tasmania <- aus_accommodation %>%
  filter(State == "Tasmania")

tasmania_plot <- tasmania %>% 
  autoplot(Takings) +
  labs(title="Accomodation Takings in Tasmania")

tasmania_ACF <- tasmania %>%
  ACF(Takings) %>%
  autoplot() + labs(title="Correlation of Takings in Tasmania ")

tasmania_PCF <- tasmania %>%
  PACF(Takings) %>%
  autoplot() + labs(title="Partial Autocorrelation of Takings in Tasmania ")

grid.arrange(tasmania_plot, tasmania_ACF, tasmania_PCF, nrow = 2)

#utilize guerrero feature to observe transformation of -0.05
lambda_Tasmania <- tasmania %>%
  features(Takings,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_Tasmania
## [1] 0.001819643
################## Tasmania Accommodation after Transformation ##################
tasmania_plot_Box_Cox<- tasmania %>%
  autoplot(difference(box_cox(Takings, lambda_Tasmania), 4)) +
  labs(title="Accomodation Takings in Tasmania After Transformation")

tasmania_ACF_Box_Cox <- tasmania %>%
  ACF(difference(box_cox(Takings, lambda_Tasmania), 4)) %>%
  autoplot() + labs(title="Correlation of Takings After Tranformation ")

tasmania_PCF_Box_Cox <- tasmania %>%
  PACF(difference(box_cox(Takings, lambda_Tasmania), 4)) %>%
  autoplot() + labs(title="Partial Autocorrelation of Takings After Transformation ")

grid.arrange(tasmania_plot_Box_Cox, tasmania_ACF_Box_Cox, tasmania_PCF_Box_Cox, nrow = 2)
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

#funciton which shows whether additional differencing is needed -  not needed in this case
tasmania %>%
  features(box_cox(Takings, lambda_Tasmania), unitroot_ndiffs)
## # A tibble: 1 x 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      1

We can see that the Tasmania Accommodation is not stationary. A Box-Cox transformation with lambda -0.05 and order differencing of 1 is applied.

(c).Monthly sales from souvenirs

################## Souvenirs  before Transformation ##################
souvenirs_plot <- souvenirs %>%
  autoplot(Sales) +
  labs(title="Monthly Sales")

souvenirs_ACF <- souvenirs %>%
  ACF(Sales) %>%
  autoplot() + labs(title="Correlation of Monthly Sales ")

souvenirs_PCF <- souvenirs %>%
  PACF(Sales) %>%
  autoplot() + labs(title="Partial Autocorrelation of Monthly Sales ")

grid.arrange(souvenirs_plot, souvenirs_ACF, souvenirs_PCF, nrow = 2)

#utilize guerrero feature to observe transformation of 0.002
lambda_souvenirs <- souvenirs %>%
  features(Sales,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_souvenirs
## [1] 0.002118221
################## Souvenirs after Transformation ##################
souvenirs_plot_Box_Cox <- souvenirs %>%
  autoplot(difference(box_cox(Sales,lambda_souvenirs), 12)) +
  labs(title="Monthly Sales After Transformation")

souvenirs_ACF_Box_Cox <- souvenirs %>%
  ACF(difference(box_cox(Sales, lambda_souvenirs), 12)) %>%
  autoplot() + labs(title="Correlation of Monthly Sales After Tranformation ")

souvenirs_PCF_Box_Cox <- souvenirs %>%
  PACF(difference(box_cox(Sales, lambda_souvenirs), 12)) %>%
  autoplot() + labs(title="Partial Autocorrelation of Monthly Sales After Transformation ")

grid.arrange(souvenirs_plot_Box_Cox, souvenirs_ACF_Box_Cox, souvenirs_PCF_Box_Cox , nrow = 2)
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

#funciton which shows whether additional differencing is needed -  not needed in this case
souvenirs %>%
  features(box_cox(Sales, lambda_Tasmania), unitroot_ndiffs)
## # A tibble: 1 x 1
##   ndiffs
##    <int>
## 1      1

We can see that the Tasmania Accommodation is not stationary. A Box-Cox transformation with lambda 0.002 and order differencing of 1 is applied.

9.5-For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

#retrieve the data from Excercise 7 in Section 2.10
set.seed(2356)
myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))


################## myseries  before Transformation ##################

#Apply the unitroot_nsdiffs() function to the data to determine if seasonal differencing is needed
myseries_plot <- myseries %>%
  autoplot(Turnover) +
  labs(title="Monthly Turnover")

myseries_ACF <- myseries %>%
  ACF(Turnover) %>%
  autoplot() + labs(title="Correlation of Monthly Turnover ")

myseries_PCF <- myseries %>%
  PACF(Turnover) %>%
  autoplot() + labs(title="Partial Autocorrelation of Monthly Turnover ")

grid.arrange(myseries_plot, myseries_ACF,myseries_PCF, nrow = 2)

#function which determines the lambda, which is approximately 0.08
lambda_myseries <- myseries %>%
  features(Turnover,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_myseries
## [1] 0.07875777
#function which determines the differencing order number
lambda_myseries <- myseries %>%
  features(Turnover,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_myseries
## [1] 0.07875777
################## myseries after Transformation ##################
myseries_plot_Box_Cox <- myseries %>%
  autoplot(difference(box_cox(Turnover, lambda_myseries), 12)) +
  labs(title="Monthly Turnover After Transformation")

myseries_ACF_Box_Cox<- myseries %>%
  ACF(difference(box_cox(Turnover, lambda_myseries), 12)) %>%
  autoplot() + labs(title="Correlation of Monthly Turnover After Tranformation ")

myseries_PCF_Box_Cox <- myseries %>%
  PACF(difference(box_cox(Turnover, lambda_myseries), 12)) %>%
  autoplot() + labs(title="Partial Autocorrelation of Monthly Turnover After Transformation ")

grid.arrange(myseries_plot_Box_Cox, myseries_ACF_Box_Cox, myseries_PCF_Box_Cox , nrow = 2)
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

It is clear that the myseries is not stationary. A Box-Cox transformation with lambda 0.08 and differencing to the first order is applied. Looking at the plots of the transformation, the data does not look fully stationary.

9.6- Simulate and plot some data from simple ARIMA models.

  1. Use the following R code to generate data from an AR(1) model with ϕ1=0.6 and σ2=1. The process starts with y1=0.

  2. Produce a time plot for the series. How does the plot change as you change Ï•1?

  3. Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1

  4. Produce a time plot for the series. How does the plot change as you change θ1?

  5. Generate data from an ARMA(1,1) model with ϕ1=0.6, θ1=0.6 and σ2=1.

  6. Generate data from an AR(2) model with ϕ1=−0.8, ϕ2=0.3 and σ2=1. (Note that these parameters will give a non-stationary series.)

  7. Graph the latter two series and compare them.

  8. Use the following R code to generate data from an AR(1) model with ϕ1=0.6 and σ2=1. The process starts with y1=0

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim_data <- tsibble(indx = seq_len(100), y = y, index = indx)
autoplot(sim_data)
## Plot variable not specified, automatically selected `.vars = y`

  1. Produce a time plot for the series. How does the plot change as you change Ï•1?
#creating y2 and y3 variables to clearly see the data behavior
set.seed(2345)
y2 <- numeric(100)
y3 <- numeric(100)
e2 <- rnorm(100)
e3 <- rnorm(100)

#utilizing formula for 0.7 and 0.8
for(i in 2:100){
    y2[i] <- 0.7*y2[i-1] + e2[i]
    y3[i] <- 0.8*y3[i-1] + e3[i]
    
}

sim <- tsibble(indx = seq_len(100), y = y, y2 = y2, y3 = y3, index = indx)

sim_plot_1 <- sim %>% autoplot(y) + labs( title = "Phi = 0.6")
sim_plot_2 <- sim %>% autoplot(y2) + labs( title = "Phi = 0.7")
sim_plot_3 <- sim %>% autoplot(y3) + labs( title = "Phi = 0.8")

grid.arrange(sim_plot_1 ,sim_plot_2,sim_plot_3, nrow = 2)

The plots shows that as there is an increase in PHi, yt becomes white noise.

  1. Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1
set.seed(3576)
y4 <- numeric(100)
e4 <- rnorm(100)

for(i in 2:100)
  y4[i] <- 0.6*e4[i-1] + e4[i]
  
sim_2 <- tsibble(indx = seq_len(100), y4 = y4, index = indx)
  1. Produce a time plot for the series. How does the plot change as you change θ1?
set.seed(6566)
set.seed(3454)
y5 <- numeric(100)
y6 <- numeric(100)
e5 <- rnorm(100)
e6 <- rnorm(100)

for(i in 2:100){
  y5[i] <- 0.7*e5[i-1] + e5[i]
  y6[i] <- 0.8*e6[i-1] + e6[i]
}
sim_2 <- tsibble(indx = seq_len(100), y4 = y4, y5 = y5, y6 = y6, index = indx)

sim_plot_4 <- sim_2 %>% autoplot(y4) + labs( title = "Theta = 0.6")
sim_plot_5 <- sim_2 %>% autoplot(y5) + labs( title = "Theta = 0.7")
sim_plot_6 <- sim_2 %>% autoplot(y6) + labs( title = "Theta = 0.8")

grid.arrange(sim_plot_4, sim_plot_5, sim_plot_6, nrow = 2)

Judging by the plots it looks like the model places greater emphasis on the most recent data points when making predictions.

  1. Generate data from an ARMA(1,1) model with ϕ1=0.6, θ1=0.6 and σ2=1.
set.seed(6785)
y7 <- (numeric(100))
e7 <- rnorm(100)

for(i in 2:100)
  y7[i] <- 0.6*y7[i-1] + 0.6*e7[i-1] + e7[i] 

sim_3 <- tsibble(indx = seq_len(100), y7 = y7, index = indx)
  1. Generate data from an AR(2) model with ϕ1=−0.8, ϕ2=0.3 and σ2=1. (Note that these parameters will give a non-stationary series.)
set.seed(2345)
y8 <- (numeric(100))
e8 <- rnorm(100)

for(i in 3:100)
  y8[i] <- -0.8*y8[i-1] + 0.3*y8[i-2] + e8[i]

sim_4 <- tsibble(indx = seq_len(100), y8 = y8, index = indx)
  1. Graph the latter two series and compare them.
sim_plot_7 <- sim_3 %>% autoplot(y7) 
sim_plot_8 <- sim_4 %>% autoplot(y8) 

grid.arrange(sim_plot_7, sim_plot_8, nrow = 2)

9.7- Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.

  1. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

  2. Write the model in terms of the backshift operator.

  3. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.

  4. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.

  5. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

  6. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

#selecting the appropriate ARIMA model
arima_aus <- aus_airpassengers %>%
  model(ARIMA(Passengers))

#checking the residuals
arima_aus%>% gg_tsresiduals()

#plotting the forecasts for the next 10 periods
arima_aus %>% forecast(h=10) %>%
  autoplot(aus_airpassengers)

The forecasts look to be on an upward trend for the next 10 years.

  1. Write the model in terms of the backshift operator.

ARIMA(0,2,1): ((1−B)^2)Yt=(1+(Θ1)B)ϵt

  1. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
#selecting the appropriate ARIMA model
arima_aus_2 <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))

#plotting the forecasts for the next 10 periods
arima_aus_2 %>% forecast(h=10) %>%
  autoplot(aus_airpassengers)

The forecasts with ARIMA(0,1,0) look to be on an upward trend for the next 10 years.

  1. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.
#selecting the appropriate ARIMA model
arima_aus_3 <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ 1 + pdq(2,1,2)))

#plotting the forecasts for the next 10 periods
arima_aus_3 %>% forecast(h=10) %>%
  autoplot(aus_airpassengers)

The forecasts with ARIMA(2,1,2) look to be on an upward trend for the next 10 years.

  1. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
#selecting the appropriate ARIMA model
arima_aus_4 <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))
## Warning: Model specification induces a quadratic or higher order polynomial trend. 
## This is generally discouraged, consider removing the constant or reducing the number of differences.
#plotting the forecasts for the next 10 periods
arima_aus_4 %>% forecast(h=10) %>%
  autoplot(aus_airpassengers)

The forecasts with ARIMA(0,2,1) look to be on an upward trend for the next 10 years. There is a warning that using a polynomial trend is discouraged. The reason is we get higher forecasts, meaning, the previous model for (d) went to 100 passengers and this one goes to 120.

9.8- For the United States GDP series (from global_economy):

  1. if necessary, find a suitable Box-Cox transformation for the data;

  2. fit a suitable ARIMA model to the transformed data using ARIMA();

  3. try some other plausible models by experimenting with the orders chosen;

  4. choose what you think is the best model and check the residual diagnostics;

  5. produce forecasts of your fitted model. Do the forecasts look reasonable?

  6. compare the results with what you would obtain using ETS() (with no transformation).

  7. if necessary, find a suitable Box-Cox transformation for the data;

#filter for the United States from global_economy series
us_economy <- global_economy %>%
  filter(Country == "United States")

#plot the Unites States from global_economy series
autoplot(us_economy)
## Plot variable not specified, automatically selected `.vars = GDP`

#checking the lambda---- it is 0.28
lambda_united_states <- us_economy %>%
  features(GDP,features=guerrero) %>%
  pull(lambda_guerrero)
lambda_united_states 
## [1] 0.2819443
#using Box-Cox and plotting the transformation with lambda 0.28
us_economy %>% autoplot(box_cox(GDP, lambda_united_states))

A Box-Cox transformation was necessary and the suitable lambda was 0.28. It is evident that the second plot had a more straightened line than the first plot.

  1. fit a suitable ARIMA model to the transformed data using ARIMA();
us_arima <- us_economy %>%
  model(ARIMA(box_cox(GDP, lambda_united_states )))

report(us_arima)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_united_states) 
## 
## Coefficients:
##          ar1  constant
##       0.4586  118.1822
## s.e.  0.1198    9.5047
## 
## sigma^2 estimated as 5479:  log likelihood=-325.32
## AIC=656.65   AICc=657.1   BIC=662.78

The ARIMA() suggests a model of: ARIMA(1,1,0) W/ drift

  1. try some other plausible models by experimenting with the orders chosen;
#some plausible models
arima_plausible <- us_economy %>%
  model(
    "ARIMA(2,1,2)" = ARIMA(box_cox(GDP, lambda_united_states) ~ 1 + pdq(2,1,2)),
    "ARIMA(0,1,0)" = ARIMA(box_cox(GDP, lambda_united_states) ~ 1 + pdq(0,1,0)),
    "ARIMA(1,1,0)" = ARIMA(box_cox(GDP, lambda_united_states) ~ 1 + pdq(1,1,0))       
      )

#comparing the AIC and BIC for each model
glance(arima_plausible) %>%
  arrange (AIC) %>%
  select(.model:BIC)
## # A tibble: 3 x 6
##   .model       sigma2 log_lik   AIC  AICc   BIC
##   <chr>         <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(1,1,0)  5479.   -325.  657.  657.  663.
## 2 ARIMA(2,1,2)  5734.   -325.  662.  664.  674.
## 3 ARIMA(0,1,0)  6774.   -332.  668.  668.  672.

We can see that ARIMA(1,1,0) has both the lowest AIC and BIC which means it is the best fit model out of the three.

  1. choose what you think is the best model and check the residual diagnostics;
us_arima %>% gg_tsresiduals()

The residuals look like white noise on the AC Plot. According to the time plot, the residuals seem to have an even distribution.

  1. produce forecasts of your fitted model. Do the forecasts look reasonable?
us_arima %>% 
  forecast(h=8) %>%
  autoplot(us_economy)

I believe the forecasts look reasonable because a trend (upward trend) can be identified.

  1. compare the results with what you would obtain using ETS() (with no transformation).
#fitting data into ETS()
us_ets<-  us_economy %>%
  model(ETS(GDP))

#plotting ETS()
us_ets %>%
  forecast(h=8) %>%
  autoplot(us_economy)

#EtS results
report(us_ets)
## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9990876 
##     beta  = 0.5011949 
## 
##   Initial states:
##          l[0]        b[0]
##  448093333334 64917355687
## 
##   sigma^2:  7e-04
## 
##      AIC     AICc      BIC 
## 3190.787 3191.941 3201.089

Comparing the ARIMA() to the ETS() model, it looks like the ARIMA() model is narrower than the ETS() model. More importantly, The ARIMA() has a lower AIC and lower BIC () than the ETS(), model making the ARIMA() model a better fit than the ETS() model in forecasting accuracy.