task2.Rmd

---
title: "Project Task 2"
author: "Aleksej Hoffaerber, Egor Zmaznev"
date: "3/4/2020"
output:
  pdf_document:
    toc: yes
    toc_depth: 3 
    df_print: paged
header-includes:
 \usepackage{booktabs}
 \usepackage{longtable}
 \usepackage{array}
 \usepackage{multirow}
 \usepackage{wrapfig}
 \usepackage{float}
 \floatplacement{figure}{H}
---

\newpage
\listoffigures
\newpage

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
# Libraries and theme----
library(fpp2)
library(ggplot2)
library(PerformanceAnalytics)
library(corrplot)
library(ggpubr)
library(seasonal)
library(splines)
library(dplyr)
library(magrittr)
library(readr)
library(skimr)
library(urca)
library(tsDyn)
library(ForecastComb)

theme_set(theme_minimal())
```


```{r loading data, echo=FALSE, message=FALSE, warning=FALSE}
# 1. Loading data ----
# Loading .csv file containing the important data points on "Net national disposable income" (di) and "Final consumption expenditure" (ce)
# for all the countries available. Filtering out unnecessary columns, and changing data format to a wide format
df <- read_csv("data.csv")

df %>% 
  filter(LOCATION == "AUS") %>% 
  select(SUBJECT, TIME, Value) %>%
  tidyr::spread(SUBJECT, Value) %>% # only calling single function to reduce conflicts between tidyr and dplyr
  rename(Date = TIME,
         di = B6NS1,
         ce = P3S1)

# creating time-series object
tse <- ts(df[,2:3], start = c(1959,3), end = c(2019,4), frequency = 4)

# binding date information from ts object to original df (easier than transforming the original format to a R date format)
df <- cbind(as.double(time(tse)), df) 

df %<>% select(-Date) %>% 
  rename(Date = "as.double(time(tse))")

# 2. Split into test and train set (10 years of prediction, 39 quarters) ----
ts.train <- window (tse,
               start = c(1959,3),
               end = c(2009,4),
               frequency = 4)

ts.test <- window (tse, 
              start = c(2010,1),
              frequency = 4)
```
%>% 

## Time-series Analysis and Summaries
### Analysis of Trend, Seasons, and Spurious Regression Issues
As subject of study serves the Australian consumption and income data. Having a look at the time series, we see a high amount of correlation between both series that needs to be analysed using the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) Test, the Augmented Dickey-Fuller (ADF) Test, and corrected using differencing and other transformations. It can be observed that the seasonality and trend components are similar for both series, as seen in _Figure 1_. The distance between consumption and income increases with the later years, meaning that people in Australia earn more than they spend (see _Figure 1_). Also, spurious regression issues can be discarded as a potential issue, as disposable available income is impacts the consumption of the population because consumption depends on earnings and debt available.

```{r Figure 1 Disposable income Analysis, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap = "Correlation and Trend of Australian Consumption and Income Data"}
# Analysis of the disposable income ----
full.plot <- autoplot(tse[,1], series = "Income") +
  autolayer(tse[,2], series = "Consumption") +
  ggtitle("Net National Disposable income and Final Consumption Expenditure", 
          subtitle = "Two time series show a high degree of correlation") +
  scale_x_continuous(name = "Years",
                     limits = c(1960, 2020),
                     breaks = seq(1960, 2020, by = 10)) +
  scale_y_continuous(name = "AUD mn.",
                     limits = c(0, 500000),
                     breaks = seq(0, 500000, by = 125000)) +
  theme(legend.position = "right")

zoom.plot <- autoplot(tse[,1], series = "Income") +
  autolayer(tse[,2], series = "Consumption") +
  ggtitle("Net National Disposable income and Final Consumption Expenditure", 
          subtitle = "High degree of correlation even clearer for trend and cycles after zoom-in") +
  scale_x_continuous(name = "Years",
                     limits = c(1990, 2020),
                     breaks = seq(1990, 2020, by = 10)) +
  scale_y_continuous(name = "AUD mn.",
                     breaks = seq(0, 500000, by = 125000))

ggarrange(full.plot, zoom.plot, nrow = 2) # Figure 1


```


### Stationarity Analysis 
The KPSS tests null hypothessis is that the data is stationary (Hyndman & Athanasopoulos, 2018). Our test statistic is with `10.8` much bigger than the 1% critical level, meaning that the data is stationary, as p>0.01. The ADF test tests also for the presence of a unit root process in the time series, with the null hypothesis indicating the stationarity of the underlying time series. We use the `drift` characteristic because of the observed characteristics of the time series and in order to test for a variety of test scenarios (Enders, 2014). ADF indicates that we fail to reject our null hypothesis for `tau2` as p>0.05, meaning that gamma is unequal zero and therefore our data is non-stationary. But, we keep our the null hypothesis for `phi1` meaning that `alpha0` and `gamma` are equal to each other and to zero (Enders, 2014), meaning that we have a stationary series and a drift component. Same applies for an ADF test with trend. Despite rejecting the null hypothesis for `tau3`, saying that gamma is unequal to zero and therefore our data stataionary, `phi2` and `phi3` indicate to a significance level of 5% that we have a `drift` and `trend` component (see _Appendix_). Based on these characteristics adn the similarities in the time series, both time series are not stationary:

```{r ADF with drift test, echo=FALSE, message=FALSE, warning=FALSE}
summary(ur.kpss(tse[,1],use.lag = 4)) 
summary(ur.df(tse[,1], type = "drift", lag = 4))
```

## Autocorrelation and Residual Analysis
Checking the residuals indicates that lots of autocorrelation remains in the residuals. In other words: valuable information that is currently not used to predict the data. This autocorrelation pattern is also typical for this type of economic data and time series. The significance of autocorrelation issue is very high, with the p-value being zero, as indicated by the Ljung-Box test in the _Appendix_. Additionally, strong trend cycles and seasonality, visible in the ACF and residual plots indicate that both series are non-stationary. Based on the distribution of the residuals we also see, that stabilisation of variance and mean are particularly important in order to design a suitable forecasting model based on the AR(I)MA process.

```{r Figure 2 residuals analysis, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap="Residual Analysis of Disposable Income"}
checkresiduals(tse[,1]) # Figure 2
```

### Disposable Income Analysis
To deal with the economic characteristic we apply a Box-Cox transformation to the data. But, this is not enough, as after the transformation not only the autocorrelation and heteroscedasticity issues remain, but also because the distibution of the residuals does not have a fitting Gauss distribution, as can be seen in _Figure 2_. To fix this, we apply differencing methods to stabilize the mean and maintain the interpretability of the model and its results. As we saw with the previous KPSS test result that rejected the null hypothesis for all p-levels from 10 to 1 percent, we need to apply a differencing method.

To calculate a fitting number of differences we must scrutinise the ACF plots: Based on the previous analyses, we know that seasons and trends play a role. A way to account for both issues and to keep the interpretability of the results is to apply first-order and seasonal differencing. In this way, we account for quarterly and seasonal (so yearly) difference in the series and can interpret the results as _quarterly changes_ in the respective variables.

As can be seen in _Appendix_, in _Figures 21 and 22_ and after applying Box-Cox transformation, seasonal, and first-order differencing, the ACF plots look less impacted by issues, such as autocorrelation, heteroscedasticity, and in total: non-stationarity. The main and variance are now stabilised. But, there are still many spikes in the ACF and residuals, but those do not follow a specific pattern. Still, the Ljung-Box test indicates that a significant amount of autocorrelation remains in the residuals, which we can not get rid off based on our pre-processing toolset. Specific assumptions and calculations in the ARIMA setting will account for the last issues. Applying the KPSS test, we see that the data is now stationary. The difference in both test results is based on autocorrelation or serial correlation, which is a much stronger indication than stationarity.

```{r Ljung-Box test and KPSS summary, echo=FALSE, message=FALSE, warning=FALSE}
Box.test(diff(diff(log(tse[,1]),4),1), lag = 10, type = "Ljung-Box")
diff(diff(log(tse[,1]),4),1) %>% ur.kpss() %>% summary()
```

To be completely sure this transformation is correct we apply KPSS functions in order to determine lag values for the differencing. Unsuprisingly, KPSS indicates to use first-order and seasonal differencing.

```{r DI ndiffs, echo=FALSE, message=FALSE, warning=FALSE, results = FALSE}
tse[,1] %>% log() %>% nsdiffs() # seasonal differencing
tse[,1] %>% log() %>% diff(lag = 4) %>% ndiffs() # first-order differencing
```

### Consumption Expenditure Analysis
Exactly the same test and procedure will be applied to the Final consumption Expenditure and seperate data will not be displayed. Because both time series look very similar (see _Figure 1_), it can be inferred that a very similar transformation must be applied. This holds true, as the plots _Appendix_, _Figure 23_ show. We incorporate both Box-Cox, seasonal and first-order differencing transformed series into the data frame and ts object.

```{r CE KPSS, echo=FALSE, message=FALSE, warning=FALSE, results=FALSE}
tse[,2] %>% ur.kpss() %>% summary()
```

```{r Log diff and seasonal, echo=FALSE, message=FALSE, warning=FALSE, results=FALSE}
Box.test(diff(diff(log(tse[,2]),4),1), lag = 10, type = "Ljung-Box")
diff(diff(log(tse[,2]),4),1) %>% ur.kpss() %>% summary()

tse[,2] %>% log() %>% nsdiffs()
tse[,2] %>% log() %>% diff(lag = 4) %>% ndiffs()

df$di.adj <- c(c(rep(NA,4), diff(diff(tse[,1]),4),1))
df$ce.adj <- c(c(rep(NA,4), diff(diff(tse[,2]),4),1))

tse <- ts(df[,2:5], start = c(1959,3), end = c(2019,4), frequency = 4)

# Split again into test and train set, with slightly adjusted start date to account for differencing reduction of the data ----
ts.train <- window(tse,
                   start = c(1960,4),
                   end = c(2009,4),
                   frequency = 4)

ts.test <- window(tse, 
                  start = c(2010,1),
                  end = c(2019,4),
                  frequency = 4)
```

##  Long-Term Relationship Analysis
As already indicated above, just regressing both variables on each other might lead to spurious regressions, identifiable by a high Adj. R-sqaured and high residual autocorrelation (Hyndman & Athanasopoulos, 2018). This occurence impacts the reliance of our forecast in the long-term horizon. Based on the previously conducted KPSS test, we already know that non-stationarity and possible cointegration of our series are an issue. These unit root tests already indicated that in order to guarantuee that the characteristic equation lies within the unit circle (Hyndman & Athanasopoulos, 2018), we must take the differences in order to assure stationarity. 

As we see, from the analyses, incorporating our transformed series for consumption into the regression does resolve our stationarity issues in the regression setting. This indicates, that both variables have a significant, cointegrated long-term relationship, which can be used to design our forecasts. This also holds true if we apply the alternative values for `tau2`, being -3.43, -2.86, and -2.57 respectively.

```{r TSLM with original data, echo=FALSE, message=FALSE, warning=FALSE, results=FALSE}
fit.tse <- tslm(formula = ts.train[,1] ~ ts.train[,2]) # not white-noise adjusted
summary(fit.tse)
summary(ur.df(residuals(fit.tse), type = "none"))
```

```{r TSLM with transformed data, echo=FALSE, message=FALSE, warning=FALSE}
fit.tse.adj <- tslm(formula = ts.train[,"di.adj"] ~ ts.train[,"ce.adj"]) # white noise adjusted
summary(fit.tse.adj)
summary(ur.df(residuals(fit.tse.adj), type = "drift", lag = 1))
```

## Identify ARIMA Model for Consumption Expenditure
### Manual ARIMA Selection Approach
Firstly, based on the previously determined time series characteristics, we know that we must employ a first-order lag differences for our autoregression (AR) part and a first-order seasonal component for our quarterly data, meaining that `d` and `D` are equal to one in order to reflect the observations and KPSS test from before. 

```{r Figure 3 log-diff, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap= "Residual Analysis for Log-Transformed and Differenced Consumption Data"}
log(ts.train[,2]) %>% 
  diff(lag = 4) %>% 
  diff() %>% 
  ggtsdisplay() # Figure 3
```

Based on the resulting assumptions, we conlude that an ARIMA(0,1,3)(0,1,2) might be suitable, as indicated in the _Appendix_ in _Figure 24 and 25_. Because of the significant spikes seen in the PACF for lag 3 and 6 we determine determine `q` = 3 and `Q` = 2. The subsequent check shows some autocorrelation left, visible in the ACF plot at spike 6. Because we already adjusted the MA(q) part, this must be a detail to be adjusted in the AR(p) part. We therefore set `p` = 1 and 2 and compare their AICc, autocorrelation, stationarity, and white noise residuals.

```{r Fit of the first models, echo=FALSE, message=FALSE, warning=FALSE}
fit.1 <- Arima(ts.train[,2], order = c(0,1,3), seasonal = c(0,1,2), lambda = BoxCox.lambda(ts.train[,2]))
fit.2 <- Arima(ts.train[,2], order = c(1,1,3), seasonal = c(0,1,2), lambda = BoxCox.lambda(ts.train[,2]))
```

```{r Figure 4 ARIMA (2,1,3)(0,1,2), echo=FALSE, message=FALSE, warning=FALSE, results = FALSE, fig.height=2, fig.width=8, fig.cap = "ACF and PACF Analysis for ARIMA (2,1,3)(0,1,2) Model"}
(fit.3 <- Arima(ts.train[,2], order = c(2,1,3), seasonal = c(0,1,2), lambda = BoxCox.lambda(ts.train[,2])))
ce.acf.3 <- ggAcf(fit.3$residuals) + ylab("") + ggtitle("ACF for ARIMA(2,1,3)(0,1,2)")
ce.pacf.3 <- ggPacf(fit.3$residuals) + ylab("") + ggtitle("PACF for ARIMA(2,1,3)(0,1,2)")
ggarrange(ce.acf.3, ce.pacf.3, ncol = 2) # Figure 4

```

```{r Figure 5 residual analysis, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual Analysis for ARIMA (2,1,3)(0,1,2) Model"}
checkresiduals(fit.3) # Figure 5
```

Comparing our models for ARIMA(1,1,3)(0,1,2) (see _Appendix_, _Figure 26 and 27_) and ARIMA(2,1,3)(0,1,2), we see that the latter has better ACF and PACF spike conditions, illustrated in _Figure 5_. Also, the Ljung-Box test is supporting this assumptions as the autocorrelation in the ARIMA(1,1,3)(0,1,2) and ARIMA(0,1,3)(0,1,2) is still highly significant as opposed by the latter model. The resulting series is now a white-noise series. The distribution of the residuals also fits the assumed distribution pattern. Looking at the coefficients, we observe that `ar1`, `ma1`, and `sma2` are not statistically significant, and `ma3` being almost not statistically significant. Because `ar1` should not be altered because of high significance of `ar2`, we try new model variants:

```{r Modified fit 4, echo=FALSE, message=FALSE, warning=FALSE}
(fit.4 <- Arima(ts.train[,2], order = c(2,1,2), seasonal = c(0,1,1), lambda = BoxCox.lambda(ts.train[,2])))
ce.acf.4 <- ggAcf(fit.4$residuals) + ylab("") + ggtitle("ACF for ARIMA(2,1,2)(0,1,1)") + theme_minimal()
ce.pacf.4 <- ggPacf(fit.4$residuals) + ylab("") + ggtitle("PACF for ARIMA(2,1,2)(0,1,1)") + theme_minimal()
```

Because in our `fit4` model (see _Appendix_, _Figure 28_) we still have some issues with autocorrelation we set `P` = 1 because of the trend that can still be observed in the residuals. We design and end up with an ARIMA(2,1,2)(1,1,1) model that also surpasses all previous models in terms of coefficient significance, and AICc, AIC, and BIC. Also the Ljung-Box p-value is maximized, the residual distibution feed the ARIMA process requirements, and the unit root theorem is also satisfied, as seen in _Figures 6 to 8_:

```{r Figure 6 model, ACF and PACF, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap = "ACF and PACF Comparison of ARIMA(2,1,2)(0,1,1) and ARIMA(2,1,2)(1,1,1)"}
(fit.5 <- Arima(ts.train[,2], order = c(2,1,2), seasonal = c(1,1,1), lambda = BoxCox.lambda(ts.train[,2])))
ce.acf.5 <- ggAcf(fit.5$residuals) + ylab("") + ggtitle("ACF for ARIMA(2,1,2)(1,1,1)") + theme_minimal()
ce.pacf.5 <- ggPacf(fit.5$residuals) + ylab("") + ggtitle("PACF for ARIMA(2,1,2)(1,1,1)") + theme_minimal()
ggarrange(ce.acf.4, ce.pacf.4,
          ce.acf.5, ce.pacf.5, ncol = 2, nrow = 2) # Figure 6
```

```{r Figure 7 residual analysis, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual and ACF Analysis for ARIMA(2,1,1)(1,1,1)"}
checkresiduals(fit.5) # Figure 8
```

```{r Figure 8 characterisitc roots, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Unit Root Theorem Analysis for ARIMA(2,1,1)(1,1,1)"}
autoplot(fit.5) # Figure 8 # no cap
```

For detailed figures and graphs showing the other scenarios, please see the _Appendix_. A comprehensive evaluation of ARIMA models can also be taken from this table:

```{r ARIMA model comparison, message=FALSE, warning=FALSE, include=FALSE, results=FALSE}
arima.comp.1 <- data.frame(model=c("ARIMA(1,1,1)(0,1,2)", 
                                    "ARIMA(1,1,3)(0,1,2)", 
                                    "ARIMA(2,1,3)(0,1,2)", 
                                    "ARIMA(2,1,2)(0,1,1)", 
                                    "ARIMA(2,1,2)(1,1,1)"),
                            LB.p.value = c(checkresiduals(fit.1)$p.value, checkresiduals(fit.2)$p.value,
                                           checkresiduals(fit.3)$p.value, checkresiduals(fit.4)$p.value,
                                           checkresiduals(fit.5)$p.value),
                           aicc = c(fit.1$aicc, fit.2$aicc, fit.3$aicc, fit.4$aicc, fit.5$aicc),
                           bic = c(fit.1$bic, fit.2$bic, fit.3$bic, fit.4$bic, fit.5$bic)) %>% 
  mutate_if(is.numeric, round, digit = 3) %>% 
  arrange(desc(LB.p.value), aicc) # ordering by p.value and AICc

```

```{r ARIMA model comparison data frame, echo=FALSE, message=FALSE, warning=FALSE}
knitr::kable(arima.comp.1,
             booktabs = T,
             caption = "Evaluation of ARIMA models without Regression Component")
```

### Automated ARIMA Selection Approach
Because KPSS can only be used to determine `d` and `D`, we need to employ Information Criteria, such as AICc, to pick the correct `p`,`q`,`P`,`Q` values. This is already incorporated in the automated ARIMA model selection that calculates different ARIMA models and picks the best models based on those Informaiton Criteria. In fact, the same ARIMA(2,1,2)(1,1,1) is picked, based on the unit root space optimization to guarantee stationarity. 

```{r AutoARIMA, message=FALSE, warning=FALSE, include=FALSE, results = FALSE}
(fit.arima <- auto.arima(ts.train[,2], lambda = BoxCox.lambda(ts.train[,2])))
summary(fit.arima) # no print

```

## Comparison of ACF and PACF with ARIMA and Raw Consumption Expenditure
Comparing the ACF and PACF plots of the raw Final consumption Expenditure and ARIMA data, we can observe the following: 1) The autocorrelation in the residuals is resolved. This was important to resolve, as ARIMA assumes that historical patterns will not change during the forecast. 2) The issue of a high PACF spike at lag 1, indicating correlation between the error terms of consumption between different lags was resolved. This is important to resolve, because ARIMA assumes uncorrelated future errors.

```{r Figure 9 ACF and PACF comparison, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap = "ACF and PACF Comparison between Original Data and ARIMA TRANSFORMED Results"}
ce.acf <- ggAcf(ts.train[,2]) + ylab("") + ggtitle("ACF for consumption Exp.")
ce.pacf <- ggPacf(ts.train[,2]) + ylab("") + ggtitle("PACF for consumption Exp.")
ggarrange(ce.acf, ce.pacf,
          ce.acf.5, ce.pacf.5, 
          ncol = 2, nrow = 2) # Figure 9
```

## ARIMA Forecast on Test Data Set
As can bee seen in _Figure 10_, the manual ARIMA model fits the data quite well, despite some minor overestimation. The prediction interval increases in size throughout time because of the included differences method.

```{r ARIMA forecast, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap = "Predicted Australian Consumption Expenditure Based on ARIMA(2,1,2)(1,1,1) Model"}
f10 <- autoplot(forecast(fit.5, h =40), series = "Forecast") +
  autolayer(ts.test[,2], series = "Actual") +
  ggtitle("Consumption Exp. Prediction, based on ARIMA without dynamic component",
          subtitle = "Forecast fits actual data quite well, despite slight overestimation") +
  scale_x_continuous(name = "Years",
                     limits = c(1990,2020))+
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,630000),
                     breaks = seq(0, 625000, by = 125000)) # Fig. 10 
f10
```

## Dynamic Regression with Explanatory Variable
### Manual Dynamic Regression Model Selection
The inclusion of a new explanatory variable in the ARIMA model requires us to check the errors terms of the regression model (eta) and our ARIMA model (epsilon). In our case, our two variables for consumption and income are cointegrated. That's why we can rely on non-stationary time series (Hyndman & Athanasopoulos, 2018). In our first model, that is already adjusted with `d` = `D` = 1, as we observed with the KPSS test before in order to guarantee non-stationarity of the data, we still observe significant ACF spikes for lag 1,3, and 4 (see _Appendix, Figure 29_), suggesting a `Q`-value of 1. PACF spikes for lag 1,3, and 4 also indicate that `P` = 1. Coefficients and AICc, and BIC values will be showed at the end.

```{r First manual dynamic regression, echo=FALSE, message=FALSE, warning=FALSE}
fit.arima.adv.1 <- Arima(ts.train[,2], 
                          order = c(0,1,0), 
                          seasonal = c(0,1,0), 
                          xreg = ts.train[,3], 
                          lambda = BoxCox.lambda(ts.train[,2]))

ce.acf.adv.1 <- ggAcf(fit.arima.adv.1$residuals) + 
  ylab("") + 
  ggtitle("ACF for DR + ARIMA(0,1,0)(0,1,0)")
ce.pacf.adv.1 <- ggPacf(fit.arima.adv.1$residuals) + 
  ylab("") + 
  ggtitle("PACF for DR + ARIMA(0,1,0)(0,1,0)")

```

```{r Comparison of ARIMA DR, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap= "ACF and PACF Analysis Comparison for ARIMA(0,1,0)(0,1,0) and ARIMA(0,1,0)(1,1,1)"}
# Comparison of the ARIMA models ----
fit.arima.adv.2 <- Arima(ts.train[,2], 
                          order = c(0,1,0), 
                          seasonal = c(1,1,1), 
                          xreg = ts.train[,3], 
                          lambda = BoxCox.lambda(ts.train[,2]))


ce.acf.adv.2 <- ggAcf(fit.arima.adv.2$residuals) + 
  ylab("") + 
  ggtitle("ACF for DR + ARIMA(0,1,0)(1,1,1)")
ce.pacf.adv.2 <- ggPacf(fit.arima.adv.2$residuals) + 
  ylab("") + 
  ggtitle("PACF for DR + ARIMA(0,1,0)(1,1,1)") 

ggarrange(ce.acf.adv.1,ce.pacf.adv.1,
          ce.acf.adv.2, ce.pacf.adv.2, 
          ncol = 2, nrow = 2) # fig 11
```

This setting shows us significant ACF and PACF spikes for lags 2 and 3 as well as 6 (see _Figure 11_) and potential for improvement for the distribution of residuals (see _Appendix, Figure 30_). We set `p` = `q` = 2, as in the previous model to balance the ACF and PACF values against each other. This results in optimal models considering ACF/PACF and white-noise behaviour, residual distribution, heteroscedasticity, stationarity, and coefficient significance (see _Figure 12_ and _Figure 13_). When looking at the coefficients, we observe that `ma1` and `xreg`, are not significant, but we include the latter because of the task. We do not delete `ma1`, as this would impact the significant `ma2` coefficient and because of the needed transformation towards autocorrelation decrease.

```{r Figure 12 ACF and PACF comparison, echo=FALSE, message=FALSE, warning=FALSE, fig.height=2, fig.width=8, fig.cap = "ACF and PACF Analysis for ARIMA(2,1,2)(1,1,1)"}
(fit.arima.adv.3 <- Arima(ts.train[,2], 
                          order = c(2,1,2), 
                          seasonal = c(1,1,1), 
                          xreg = ts.train[,3], 
                          lambda = BoxCox.lambda(ts.train[,2])))

ce.acf.adv.3 <- ggAcf(fit.arima.adv.3$residuals) + 
  ylab("") + 
  ggtitle("ACF for DR + ARIMA(2,1,2)(1,1,1)")

ce.pacf.adv.3 <- ggPacf(fit.arima.adv.3$residuals) + 
  ylab("") + 
  ggtitle("PACF for DR + ARIMA(2,1,2)(1,1,1)")

ggarrange(ce.acf.adv.3,ce.pacf.adv.3, ncol = 2) # fig 12
```

```{r Figure 13 residual analysis DR, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual Analysis for ARIMA(2,1,2)(1,1,1)"}
checkresiduals(fit.arima.adv.3) # fig 13
```

```{r Figure 14 errors comparison, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Error Comparison"}
cbind("Reg Err." = residuals(fit.arima.adv.3, type = "regression"),
      "ARIMA Err." = residuals(fit.arima.adv.3, type = "innovation")) %>% 
  autoplot(facets = T) +
  ggtitle("Comparison of DR + ARIMA(2,1,2)(1,1,1) Errors",
          subtitle = "Regression Errors capute the overall time series trend, ARIMA errors resemlbe a white noise series")
# fig 14 #no cap
```

This reestimation yields in a suitable model, considering the white noise type of ARIMA residuals, ACF and PACF specifics, as well as a fitting residual distribution that is only slightly skewed because of the observation outliers during the financial crisis in 2007/08. Additionally, we could include a constant in order to mimick the trend that is displayed in our regression residuals (see _Figure 14_). For this, and because a drift cannot be included if the order of difference > 2, we must set d = 0 and also q = 0, because this drift should explain the information conveyed in the regression residuals. But because this change yields in more autocorrelation, we refrain from doing so. 

### Automated ARIMA Model Selection
On the other side, the automated approach yields in a different model variation, that was already discussed above but discarded because of its negative impact on ACF and PACF plots and white-noise properties. It yields a lower AICc and does not yield in autocorrelation reduction, as seen in _Figure 15_. In sum, the automated dynamic regression model is a worse forecast model than our ARIMA(2,1,2)(1,1,1) model, as can be seen in _Figure 16_.

```{r Figure 15 errors comparison, echo=FALSE, fig.cap="Residual Analysis for Automated Selected ARIMA Model", fig.height=3, fig.width=8, message=FALSE, warning=FALSE}
fit.arima.adv <- auto.arima(ts.train[,2], xreg = ts.train[,3], lambda = BoxCox.lambda(ts.train[,2]))
checkresiduals(fit.arima.adv) # fig 15
```

```{r DR model comparison, message=FALSE, warning=FALSE, include=FALSE, results=FALSE}
arima.comp.2 <- data.frame(model=c("ARIMA(0,1,0)(0,1,0)", 
                                   "ARIMA(0,1,0)(1,1,1)", 
                                   "ARIMA(2,1,2)(1,1,1)", 
                                   "ARIMA(1,1,1)(0,0,2)"),
                           LB.p.value = c(checkresiduals(fit.arima.adv.1)$p.value, 
                                          checkresiduals(fit.arima.adv.2)$p.value,
                                          checkresiduals(fit.arima.adv.3)$p.value, 
                                          checkresiduals(fit.arima.adv)$p.value),
                           aicc = c(fit.arima.adv.1$aicc, 
                                    fit.arima.adv.2$aicc, 
                                    fit.arima.adv.3$aicc, 
                                    fit.arima.adv$aicc),
                           bic = c(fit.arima.adv.1$bic, 
                                   fit.arima.adv.2$bic, 
                                   fit.arima.adv.3$bic, 
                                   fit.arima.adv$bic)) %>%
  mutate_if(is.numeric, round, digit = 3) %>% 
  arrange(desc(LB.p.value), aicc) # ordering by p.value and AICc

```

```{r DR plot comparison, echo=FALSE, message=FALSE, warning=FALSE}
knitr::kable(arima.comp.2,
             booktabs = T,
             caption = "Evaluation of ARIMA models with regression")
```

```{r Figure 17 forecasts comparison, echo=FALSE, message=FALSE, warning=FALSE, fig.height=6, fig.width=8, fig.cap = "Forecast Performance Comparison Between Manual and Automated Approach"}
# Forecasting Comparison
fcs.adv.3 <- forecast(fit.arima.adv.3, xreg = ts.test[,3])
fcs.adv <- forecast(fit.arima.adv, xreg = ts.test[,3])

fcs.adv.3.plot <- autoplot(fcs.adv.3, series = "Fitted consumption") +
  autolayer(ts.test[,2], series = "Actual consumption") +
  scale_x_continuous(name = "Years",
                     limits = c(1990,2020))+
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,630000),
                     breaks = seq(0, 625000, by = 125000))

fcs.adv.plot <- autoplot(fcs.adv, series = "Fitted consumption") +
  autolayer(ts.test[,2], series = "Actual consumption") +
  scale_x_continuous(name = "Years",
                     limits = c(1990,2020))+
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,625000),
                     breaks = seq(0, 625000, by = 125000))

ggarrange(fcs.adv.3.plot, fcs.adv.plot, nrow = 2, common.legend = T, legend = "bottom") # fig 17
```

## Forecast Combination: Comparison, Plotting and Measurment
In the following part, we will combine different forecasting models and assess its performance. Specifically, we will test 3 different cases in forecasting combination:

1. ARIMA model combined with the dynamic regression(DR).
In this scenario, we will combine our best performing ARIMA model (in our case, it is `ARIMA(2,1,2)(1,1,1)[4]`) and same ARIMA model with the specified `xreg` parameter.

2. AutoARIMA with regression model and manually picked ARIMA with regression.
Previously, the automated model has been compared with the manually adjusted one, and we have seen the advantage of the manual approach. Yet, it is interesting to see the performance of the combined model and to compare it to other variations of the forecast combination.

3. Manual ARIMA with regression and TBATS/ETS
In this case, we want to combine the best performing ARIMA model with other forecasting approaches. In this case, for the comparison, we decided to compare the combination with ETS and TBATS models.

Every forecast combination has two variations: averaged forecast combination and using optimal weights. Optimal weights have been calculated by the following formulas: $w=\frac{\sigma_{2}^{2}-\sigma_{12}}{\sigma_{1}^{2}+\sigma_{2}^{2}-2\sigma_{12}}$ and $1-w=\frac{\sigma_{1}^{2}-\sigma_{12}}{\sigma_{1}^{2}+\sigma_{2}^{2}-2\sigma_{12}}$

```{r Forecast combination, message=FALSE, warning=FALSE, include=FALSE}
# Dynamic and ARIMA ----
# Averaging: dynamic regression and ARIMA
fcs.5 <- forecast(fit.5, h = 40)
comb_dr <- (fcs.5[["mean"]]+fcs.adv.3[["mean"]])/2

accuracy(fcs.5, x = ts.test[,2])
accuracy(fcs.adv.3, x = ts.test[,2])
accuracy(comb_dr, x = ts.test[,2])

# Optimal weights: dynamic regression and ARIMA
w <- (var(fcs.adv.3[["mean"]]) - sd(fcs.adv.3[["mean"]]+fcs.5[["mean"]]))/(var(fcs.5[["mean"]])+var(fcs.adv.3[["mean"]])-2*sd(fcs.adv.3[["mean"]]+fcs.5[["mean"]]))
rw <- 1-w

comb_dr_w <- w*fcs.5[["mean"]]+rw*fcs.adv.3[["mean"]]
accuracy(comb_dr_w, x = ts.test[,2])

# ARIMA and autoARIMA comparison ----
# Averaging: Auto and Manual
comb_am <- (fcs.adv[["mean"]]+fcs.adv.3[["mean"]])/2

accuracy(fcs.adv, x = ts.test[,2])
accuracy(fcs.adv.3, x = ts.test[,2])
accuracy(comb_am, x = ts.test[,2])

# Optimal weights: Auto and Manual
w <- (var(fcs.adv.3[["mean"]]) - sd(fcs.adv.3[["mean"]]+fcs.adv[["mean"]]))/(var(fcs.adv[["mean"]])+var(fcs.adv.3[["mean"]])-2*sd(fcs.adv.3[["mean"]]+fcs.adv[["mean"]]))
rw <- 1-w

comb_am_w <- w*fcs.adv[["mean"]]+rw*fcs.adv.3[["mean"]]
accuracy(comb_am_w, x = ts.test[,2])

# With ETS and TBATS ----
# Averaging: ETS, TBATS
ETS <- forecast(ets(ts.train[,2]), h = 40)
TBATS <- forecast(tbats(ts.train[,2], biasadj = T), h = 40)
comb_tbats <- (fcs.adv.3[["mean"]]+TBATS[["mean"]])/2
comb_ets <- (fcs.adv.3[["mean"]]+ETS[["mean"]])/2

accuracy(comb_ets, x = ts.test[,2])
accuracy(comb_tbats, x = ts.test[,2])

# Optimal weights: ETS, TBATS
# For TBATS
w <- (var(fcs.adv.3[["mean"]]) - sd(fcs.adv.3[["mean"]]+TBATS[["mean"]]))/(var(TBATS[["mean"]])+var(fcs.adv.3[["mean"]])-2*sd(fcs.adv.3[["mean"]]+TBATS[["mean"]]))
rw <- 1-w

comb_tbats_w <- w*TBATS[["mean"]]+rw*fcs.adv.3[["mean"]]
accuracy(comb_tbats_w, x = ts.test[,2])

# For ETS
w <- (var(fcs.adv.3[["mean"]]) - sd(fcs.adv.3[["mean"]]+ETS[["mean"]]))/(var(ETS[["mean"]])+var(fcs.adv.3[["mean"]])-2*sd(fcs.adv.3[["mean"]]+ETS[["mean"]]))
rw <- 1-w

comb_ets_w <- w*ETS[["mean"]]+rw*fcs.adv.3[["mean"]]
accuracy(comb_ets_w, x = ts.test[,2])
```

_Table 3_ below provides the comparison of both sole and combined models and measures of ME, RMSE and MAE for the test set. In the performance assessment and comparison, we are primarily looking at the RMSE, therefore, the ordered comparison of the model performance can be seen in the _Figure 20_ below.

```{r Table with comparisons, echo=FALSE, message=FALSE, warning=FALSE}
all_comp <- as.data.frame(matrix(data = c(accuracy(fcs.5, x = ts.test[,2])["Test set", c("ME", "RMSE", "MAE")],
                            accuracy(fcs.adv, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(fcs.adv.3, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(ETS, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(TBATS, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_dr, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_dr_w, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_am, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_am_w, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_ets, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_ets_w, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_tbats, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")],
                            accuracy(comb_tbats_w, x = ts.test[,2])["Test set",c("ME", "RMSE", "MAE")]),
                   nrow = 13,
                   ncol = 3,
                   byrow = T,
                   dimnames = list(c("ARIMA(2,1,2)(1,1,1)", 
                                     "ARIMA(1,1,1)(0,0,2) with reg",
                                     "ARIMA(2,1,2)(1,1,1) with reg",
                                     "ETS", 
                                     "TBATS",
                                     "Averaged ARIMA+DR",
                                     "Optimal weights: ARIMA+DR",
                                     "Averaged: DRs(auto+manual)",
                                     "Optimal: DRs(auto+manual)",
                                     "Averaged: DR + ETS",
                                     "Optimal weights: DR + ETS",
                                     "Averaged: DR + TBATS",
                                     "Optimal weights: DR + TBATS"), 
                                   c("ME", "RMSE", "MAE"))))

knitr::kable(all_comp,
             booktabs = T,
             caption = "Evaluation of forecasting of sole and combined models")
```

As for the first combination case, where the ARIMA model is combined with the same model with regression, we can observe that the RMSE for the sole models is comparably close. So, in fact, there is no direct need to combine the forecasts, as the performance of the combined models and sole models will differ just slightly. 

```{r ARIMA and DR plot, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap="Forecast Combination for Dynamic Regressions Using ARIMA"}
# ARIMA + Dynamic ARIMA ----
c_1 <- autoplot(tse[,2]) +
  autolayer(comb_dr, series = "Averaged forecast")+
  autolayer(fcs.5, series = "ARIMA(2,1,2)(1,1,1)", PI = F) +
  autolayer(fcs.adv.3, series = "ARIMA(2,1,2)(1,1,1) with xreg", PI = F) +
  autolayer(comb_dr_w, series = "Optimal weights combination")+
  ggtitle("Forecasts combination: ARIMA and dynamic regression", 
          subtitle = "Forecast are overlapping because of predominance of the DR + ARIMA(2,1,2)(1,1,1)") +
  scale_x_continuous(name = "Year",
                     limits = c(2000,2020),
                     breaks = seq(2000, 2020, by = 5)) +
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,400000),
                     breaks = seq(0, 400000, by = 75000))+
  scale_color_discrete(name = "Forecasting Models")

c_1 # fig 17
```

From the _Figure 17_, we can see that all of the forecasting models perform almost equally as good, with just a slight advantage of the ARIMA(2,1,2)(1,1,1) with regression. This can be explained with the fact that both models are providing forecasts with estimations "above" the actual forecast and with a slight difference comparing to the models we will observe in the following two cases. Also, the 4 forecast are incredibly similar because they are all based on ARIMA(2,1,2)(1,1,1) model.

For the second combination, we decided to combine `auto.arima` with regression and manually adjusted model with regression, used in the previous case. As we previously observed, the difference between both cases is observable with the measured RMSE of 19351.334 and 52235.143, for the manual and auto models correspondingly.

```{r Manual and autoARIMA plot, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap="Forecast Combination for Dynamic Regressions and Automated Regression"}
# ARIMA Manual + ARIMA Auto ----
c_2 <- autoplot(tse[,2]) +
  autolayer(comb_am, series = "Averaged forecast")+
  autolayer(fcs.adv, series = "ARIMA(1,1,1)(0,0,2) with xreg", PI = F)+
  autolayer(fcs.adv.3, series = "ARIMA(2,1,2)(1,1,1) with xreg", PI = F)+
  autolayer(comb_am_w, series = "Optimal weights combination")+
  ggtitle("Forecasts combination: autoARIMA and manual approach", 
          subtitle = "Comparison of the sole models and forecast combinations") +
  scale_x_continuous(name = "Year",
                     limits = c(2000,2020),
                     breaks = seq(2000, 2020, by = 5)) +
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,400000),
                     breaks = seq(0, 400000, by = 75000))+
  scale_color_discrete(name = "Forecasting Models")

c_2 # fig 18
```

From the _Figure 18_, presented above, we can observe the similar pattern as in the previous case: combination model with optimal weights performs slightly better than the averaged model, yet the sole manually adjusted model is still seen as the best performer due to both forecasts laying above the actual data of the test set.

In the previous cases, we have looked into the combination of the different variations of the ARIMA models that we have estimated in the for the previous tasks. All of them share similar characteristics being above the test set in terms of forecast and primarily varying in the "closeness" to the actual data and capturing of the seasonality (especially visible in the last case). However, it is interesting to see the difference between the ARIMA models and other estimation approaches. For this case, we have made forecasting using TBATS and ETS models both separately combining with the best performing ARIMA model with regression.

```{r fig.height=6, message=FALSE, warning=FALSE, echo=FALSE, fig.height=6, fig.width=9.2, fig.cap="Forecast Comparison for ETS, TBATS Approaches"}
# ARIMA+ETS and ARIMA+TBATS ----
c_3 <- autoplot(tse[,2]) +
  autolayer(fcs.adv.3, series = "ARIMA with regression", PI = F) +
  autolayer(ETS, series = "ETS", PI = F)+
  autolayer(TBATS, series = "TBATS", PI = F)+
  ggtitle("Performance of the models without combination", 
          subtitle = "Comparison of ARIMA with regression, TBATS and ETS") +
  scale_x_continuous(name = "Year",
                     limits = c(2000,2020),
                     breaks = seq(2000, 2020, by = 5)) +
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,400000),
                     breaks = seq(0, 400000, by = 75000))+
  scale_color_discrete(name = "Models:")+ 
  theme(legend.position="bottom")

c_4 <- autoplot(tse[,2]) +
  autolayer(comb_tbats, series = "Averaged: TBATS+ARIMA")+
  autolayer(comb_ets, series = "Averaged: ETS+ARIMA")+
  autolayer(comb_ets_w, series = "Optimal weights: ETS+ARIMA") +
  autolayer(comb_tbats_w, series = "Optimal weights: TBATS+ARIMA")+
  ggtitle("Forecasts combination: TBATS and ETS with ARIMA", 
          subtitle = "Combining ARIMA with regression with the other approaches") +
  scale_x_continuous(name = "Year",
                     limits = c(2000,2020),
                     breaks = seq(2000, 2020, by = 5)) +
  scale_y_continuous(name = "consumption Exp. (in mn AUD)",
                     limits = c(0,400000),
                     breaks = seq(0, 400000, by = 75000))+
  scale_color_discrete(name = "Models:")+ 
  theme(legend.position="bottom")

ggarrange(c_3, c_4, nrow = 2) # fig 19
```

From _Figure 19_ above, we can see that TBATS model and ARIMA with regression are very similar in its estimations, which is logical as far as TBATS is similar to the dynamic harmonic regression in many aspects, apart from seasonality. Different estimation characteristics are observed for the ETS (M,A,M): the model itself performs very well for based on the training set, and we can see that it captures data trend and seasonality quite well. Apart from that, ETS estimations are mapped slightly below the actual data, which should be very beneficial in the combination of two models.

The results of the combined models can be found in the lower part of the _Figure 19_. High performance of the forecast combination can be observed: a combination of the ARIMA with ETS has the best results among all models. Applying the use of the optimal weights, allows us to enhance the model even further, producing the overall RMSE of 4098.904. TBATS model forecast combination with ARIMA is weaker, due to the aforementioned estimation characteristics of TBATS. In the _Figure 20_ below, we can see that overall TBATS model, as with many automated modelling frameworks, performs rather poorly with just a slighter better estimation, than autoARIMA.

```{r Figure Measurements comparison, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=8, fig.cap="Forecast Accuracy Comparison"}
all_comp <- all_comp[order(all_comp$RMSE),]
all_comp$Model <- factor(row.names(all_comp), levels = row.names(all_comp)) 

ggplot(all_comp, aes(x=Model, y= RMSE, label=RMSE)) +
  scale_x_discrete(limits = rev(levels(all_comp$Model)))+
  geom_bar(stat='identity', aes(fill=row.names(all_comp)), width=.5)+
  theme(legend.position = "none")+
  labs(subtitle="ETS and combinations with ETS perform the best followed by sole MDR", 
       title= "Comparison of forecasts accuracy, based on RMSE") + 
  coord_flip()+
  labs(caption = "DR - Dynamic Regression (in our case: ARIMA(2,1,2)(1,1,1) with xreg") # fig 20
```

From _Figure 20_, we can see the overall comparison of all the models (combined and sole) performance based on the RMSE. Noticeably, the combination of two different best performing sole models (ARIMA with regression and ETS), provides us with the best results. In particularly this case, the poor performance of the automated modelling frameworks is visible with TBATS and autoARIMA having the biggest RMSE. In all the cases of forecasts combination, the advantage of the application of the optimal weights over averaging can be mentioned.

\newpage

## Appendix

```{r Appendix 1 Stationarity Analysis, echo=FALSE, message=FALSE, warning=FALSE}
summary(ur.df(tse[,1], type = "trend", lag = 1)) 
```

```{r Appendix 2 Ljung-Box test, echo=FALSE, message=FALSE, warning=FALSE}
Box.test(tse[,1], lag = 10, type = "Ljung-Box") # Appendix 
Box.test(tse[,2], lag = 10, type = "Ljung-Box") # Appendix 
```

```{r Appendix 3 residual analysis CE, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual Analysis for Income"}
checkresiduals(tse[,2]) # fig 21
```

```{r Appendix 4 second order diff-log analysis, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual Analysis for First-Order Seasonal Differencing with Box-Cox for Income"}
checkresiduals(diff(diff(log(tse[,1]),4),1)) # fig 22
```

```{r Appendix 5 second order diff-log CE, echo=FALSE, message=FALSE, warning=FALSE, fig.height=3, fig.width=8, fig.cap = "Residual Analysis for First-Order Seasonal Differencing for Consumption"}
checkresiduals(diff(diff(log(tse[,2]),4),1)) # fig 23
```

```{r Appendix 6 first ARIMA model, fig.height=2, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "Residual Analysis for First ARIMA(0,1,3)(0,1,2) Model"}
(fit.1 <- Arima(ts.train[,2], order = c(0,1,3), seasonal = c(0,1,2), lambda = BoxCox.lambda(ts.train[,2])))
ce.acf.1 <- ggAcf(fit.1$residuals) + ylab("") + ggtitle("ACF for ARIMA(0,1,3)(0,1,2)")
ce.pacf.1 <- ggPacf(fit.1$residuals) + ylab("") + ggtitle("PACF for ARIMA(0,1,3)(0,1,2)")
ggarrange(ce.acf.1, ce.pacf.1, ncol = 2) # fig 24
```

```{r Appendix 7 residual analysis, fig.height=3, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "Residual Analysis for First ARIMA(1,1,3)(0,1,2) Model"}
checkresiduals(fit.1)$p.value # fig 25
```

```{r Appendix 8 second ARIMA model, fig.height=2, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "PACF & ACF Analysis for Second ARIMA(1,1,3)(0,1,2) Model"}
(fit.2 <- Arima(ts.train[,2], order = c(1,1,3), seasonal = c(0,1,2), lambda = BoxCox.lambda(ts.train[,2])))
ce.acf.2 <- ggAcf(fit.2$residuals) + ylab("") + ggtitle("ACF for ARIMA(1,1,3)(0,1,2)")
ce.pacf.2 <- ggPacf(fit.2$residuals) + ylab("") + ggtitle("PACF for ARIMA(1,1,3)(0,1,2)")
ggarrange(ce.acf.2, ce.pacf.2, ncol = 2) # fig 26
```

```{r Appendix 9 residual analysis, fig.height=3, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "Residual  Analysis for Second ARIMA(1,1,3)(0,1,2) Model"}
checkresiduals(fit.2) # fig 27
```

```{r Appendix 10 ACF and PACF, fig.height=4, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "PACF & ACF  Analysis for Third ARIMA(2,1,3)(0,1,2) and Fourth ARIMA(2,1,2)(0,1,1) Model"}
ggarrange(ce.acf.3, ce.pacf.3, 
          ce.acf.4, ce.pacf.4,
          ncol = 2, nrow = 2) # fig 28
```

```{r Appendix 11 DR residual analysis, fig.height=3, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "Residual Analysis for First DR + ARIMA(0,1,0)(0,1,0) Model"}
checkresiduals(fit.arima.adv.1) # fig 29 
```

```{r Appendix 12, fig.height=3, fig.width=8, echo=FALSE, message=FALSE, warning=FALSE, fig.cap = "Residual Analysis for Second DR + ARIMA(0,1,0)(1,1,1) Model"}
checkresiduals(fit.arima.adv.2) # fig 30
```

## References

* Enders, Walter (2014) Applied Econometric Time Series, 4th Edition, ISBN: 978-1-118-80856-6, Accessed on 01.04.2020
* Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. OTexts.com/fpp2. Accessed on 29.03.2020