Chapter8.rmd

# Chapter 8

```{r echo=FALSE, message=FALSE, warning=FALSE, Load_packages}

library(fpp2)
library(xlsx)
library(rdatamarket)
library(tseries)

```

2. A classic example of a non-stationary series is the daily closing IBM stock price series (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

```{r echo=FALSE, message=FALSE, warning=FALSE, Question2}

ggtsdisplay(ibmclose)
# ACF plot shows that the autocorrelation values are bigger than critical value and decrease slowly. Also, r1 is large(near to 1) and positive. It means that the IBM stock data are non-stationary(that is, predictable using lagged values).
# PACF plot shows that there is a strong correlation between IBM stock data and their 1 lagged values. It means that IBM stock data can be predicted by 1 lagged values and they aren't stationary.

# To get stationary data, IBM stock data need differencing. Differencing can help stabilize the mean of a time series by removing changes in the level of a time series. Therefore it will eliminate or reduce trend and seasonality. And the effect can make non-staionary data stationary.

```


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

```{r echo=FALSE, message=FALSE, warning=FALSE, Question3}

# a. usnetelec
autoplot(usnetelec)
# It is almost linearly increasing data. It looked like that the data only need first differencing.
Box.test(diff(usnetelec), type = "Ljung-Box")
# first differenced usnetelec data can be thought of as a white noise series.
kpss.test(diff(usnetelec))
# kpss test result also shows that first differencing made the data stationary.

# b. usgdp
autoplot(usgdp)
# It is almost linearly increasing data. It looked like that the data only need first differencing.
Box.test(diff(usgdp), type = "Ljung-Box")
# first differenced usnetelec data cannot be thought of as a white noise series.
autoplot(diff(usgdp))
# There is still a trend left in the differenced data. It looked like one more differencing would be enough, but use ndiffs function to check the number of differencing needed.
ndiffs(usgdp)
# One more differencing would be enough.
autoplot(diff(diff(usgdp)))
# Plot shows that the twice differenced data is like white noise series.
Box.test(diff(diff(usgdp)), type = "Ljung-Box")
# But it couldn't pass Ljung-Box test.
ggAcf(diff(diff(usgdp)))
# There are still some autocorrelations left.
kpss.test(diff(diff(usnetelec)))
# But kpss test result shows that differencing twice was enough to make the data stationary. Therefore in usgdp data case, even if twice differencing didn't make the data like white noise series, it made the data stationary.

# c. mcopper
autoplot(mcopper)
# mcopper data have increasing trend. And they have bigger variation for bigger prices. Therefore I'll use Box-Cox transformation before differencing.
lambda_mcopper <- BoxCox.lambda(mcopper)
autoplot(diff(BoxCox(mcopper, lambda_mcopper)))
Box.test(diff(BoxCox(mcopper, lambda_mcopper)),
         type = "Ljung-Box")
# Plot result looked like BoxCox transformation and first differencing made the data like white noise series. But Ljung-Box test shows that it didn't.
ggAcf(diff(BoxCox(mcopper, lambda_mcopper)))
# There are still some autocorrelations left.
kpss.test(diff(BoxCox(mcopper, lambda_mcopper)))
# But kpss test result shows that differencing with Box-Cox transformation was enough to make the data stationary. 
# Even if differencing with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.

# d. enplanements
autoplot(enplanements)
# enplanements data have seasonality and increasing trend even if the number of enplanements fell in 2001. Therefore, I think that the data need seasonal differencing, too. 
# The variations are bigger for bigger numbers. Therefore I'll use Box-Cox transformation before differencing.
lambda_enplanements <- BoxCox.lambda(enplanements)
ndiffs(enplanements)
nsdiffs(enplanements)
# the data need 1 first differencing and 1 seasonal differencing.
autoplot(
  diff(
    diff(
      BoxCox(enplanements, lambda_enplanements),
      lag = 12
    )
  )
)
Box.test(
  diff(
    diff(
      BoxCox(enplanements, lambda_enplanements),
      lag = 12
    )
  ),
  type = "Ljung-Box"
)
# Plot result looked like BoxCox transformation and multiple differencings made the data like white noise series. But Ljung-Box test shows that it didn't.
ggAcf(
  diff(
    diff(
      BoxCox(enplanements, lambda_enplanements),
      lag = 12
    )
  )
)
# There are still some autocorrelations left.
kpss.test(
  diff(
    diff(
      BoxCox(enplanements, lambda_enplanements),
      lag = 12
    )
  )
)
# But kpss test result shows that differencings with Box-Cox transformation was enough to make the data stationary. In enplanements data case, even if differencings with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.

# e. visitors
autoplot(visitors)
# visitors data are similar to enplanements data. They have seasonality and increasing trend. It looked like they also need Box-Cox transformation, first and seasonal differencing.
lambda_visitors <- BoxCox.lambda(visitors)
ndiffs(visitors)
nsdiffs(visitors)
# visitors data need 1 first and 1 seasonal differencing.
autoplot(
  diff(
    diff(
      BoxCox(visitors, lambda_visitors),
      lag = 12
    )
  )
)
Box.test(
  diff(
    diff(
      BoxCox(visitors, lambda_visitors),
      lag = 12
    )
  ),
  type = "Ljung-Box"
)
# Plot result looked like BoxCox transformation and multiple differencings made the data like white noise series. But Ljung-Box test shows that it didn't.
ggAcf(
  diff(
    diff(
      BoxCox(visitors, lambda_visitors),
      lag = 12
    )
  )
)
# There are still some autocorrelations left.
kpss.test(
  diff(
    diff(
      BoxCox(visitors, lambda_visitors),
      lag = 12
    )
  )
)
# But kpss test result shows that differencings with Box-Cox transformation was enough to make the data stationary. In visitors data case, even if differencings with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.

```


4. For the enplanements data, write down the differences you chose above using backshift operator notation.

### the data needed 1 first difference, 1 seasonal difference after Box-Cox transformation. The model of the data can be written as ARIMA(0, 1, 0)(0, 1, 0)12 with Box-Cox transformation(lambda = -0.227).
### The model expression using backshift operator notation B: 
### first equation : wt = (yt^(-0.227) - 1)/(-0.227)
### second equation : (1 - B)(1 - B^12)wt = et, where et is a white noise series.


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

```{r echo=FALSE, message=FALSE, warning=FALSE, Question5}

retail <- read.xlsx("retail.xlsx",
                    startRow = 2,
                    sheetIndex = 1)

retail.ts <- ts(retail[,"A3349873A"], 
                frequency=12, 
                start=c(1982,4))

autoplot(retail.ts)
# the data have increasing trend and strong seasonality. And there are bigger variations for bigger numbers. Therefore I think that I need to use first differencing and seasonal differencing. And it would be better to do Box-Cox transformation.
ndiffs(retail.ts)
nsdiffs(retail.ts)
# I'm going to do 1 first differencing and 1 seasonal differencing.
kpss.test(
  diff(
    diff(
      BoxCox(retail.ts, BoxCox.lambda(retail.ts)),
      lag = 12
    )
  )
)
# To make retail.ts data stationary, I did Box-Cox transformation, 1 first differencing and 1 seasonal differencing.

```

6. Use R to simulate and plot some data from simple ARIMA models.

```{r echo=FALSE, message=FALSE, warning=FALSE, Question6}

# a. Use the following R code to generate data from an AR(1) model with phi1 = 0.6 and sigma^2 = 1. The process starts with y1 = 0.

y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
   y[i] <- 0.6*y[i-1] + e[i]
}

# b. Produce a time plot for the series. How does the plot change as you change phi1?
ar1generator <- function(phi1){
  # generate 100 data points from an AR(1) model with input phi1.
  y <- ts(numeric(100))
  # error 'e's have variation sigma^2 as 1.
  e <- rnorm(100)
  for(i in 2:100){
    y[i] <- phi1*y[i-1] + e[i]
  }
  return(y)
}

# produce plots changing phi1 value.
autoplot(ar1generator(0.3), series = "0.3") +
  geom_line(size = 1, colour = "red") +
  autolayer(y, series = "0.6", size = 1) +
  autolayer(ar1generator(0.9), size = 1, series = "0.9") +
  ylab("AR(1) models") +
  guides(colour = guide_legend(title = "Phi1"))
# As phi increases, the variation of y increased.
    
# c. Write your own code to generate data from an MA(1) model with theta1 = 0.6 and sigma^2 = 1.
ma1generator <- function(theta1){
  # generate 100 data points from an MA(1) model with input theta1.
  y <- ts(numeric(100))
  # error 'e's have variation sigma^2 as 1.
  e <- rnorm(100)
  for(i in 2:100){
    y[i] <- theta1*e[i-1] + e[i]
  }
  return(y)
}
    
# d. Produce a time plot for the series. How does the plot change as you change theta1?
# produce plots changing theta1 value.
autoplot(ma1generator(0.3), series = "0.3") +
  geom_line(size = 1, colour = "red") +
  autolayer(y, series = "0.6", size = 1) +
  autolayer(ar1generator(0.9), size = 1, series = "0.9") +
  ylab("MA(1) models") +
  guides(colour = guide_legend(title = "Theta1"))
# As theta increases, the variation of y increased.
    
# e. Generate data from an ARMA(1,1) model with phi1 = 0.6, theta1 = 0.6 and sigma^2 = 1.
y_arima.1.0.1 <- ts(numeric(50))
e <- rnorm(50)
for(i in 2:50){
   y_arima.1.0.1[i] <- 0.6*y_arima.1.0.1[i-1] + 0.6*e[i-1] + e[i]
}
    
# f. Generate data from an AR(2) model with phi1 = -0.8, phi2 = 0.3 and sigma^2 = 1. (Note that these parameters will give a non-stationary series.)
y_arima.2.0.0 <- ts(numeric(50))
e <- rnorm(50)
for(i in 3:50){
   y_arima.2.0.0[i] <- -0.8*y_arima.2.0.0[i-1] + 0.3*y_arima.2.0.0[i-2] + e[i]
}
    
# g. Graph the latter two series and compare them.
autoplot(y_arima.1.0.1, series = "ARMA(1, 1)") +
  autolayer(y_arima.2.0.0, series = "AR(2)") +
  ylab("y") +
  guides(colour = guide_legend(title = "Models"))

autoplot(y_arima.1.0.1)
# data from an AR(2) model increased with oscillation. They are non-staionary data. But data from an ARMA(1, 1) model were stationary.

```


7. Consider the number of women murdered each year (per 100,000 standard population) in the United States. (Data set wmurders).

```{r echo=FALSE, message=FALSE, warning=FALSE, Question7}

# a. By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data. 
autoplot(wmurders)
# It looked like the data don't need seasonal differencing or Box-Cox transformation.
autoplot(diff(wmurders))
# It looked like 1 more differencing would be needed to make the data stationary. Differenced data slowly go to minus infinity.
ndiffs(wmurders)
# ndiffs function shows that the data need 2 differencing.
autoplot(diff(wmurders, differences = 2))
kpss.test(diff(wmurders, differences = 2))
# twice differencing made the data stationary.
diff(wmurders, differences = 2) %>% ggtsdisplay()
# PACF is decaying. And there are significant spikes at lag 1, and 2 in the ACF, but none beyond lag 2. If the data can be modelled by ARIMA(0, 2, q) or ARIMA(p, 2, 0), I'm going to model the data by ARIMA(0, 2, 2).

# b. Should you include a constant in the model? Explain.
# ARIMA model of the data includes twice differencing. If there is a constant in the model, twice integrated contant will yield quadratic trend, which is dangerous for forecasting. Therefore I won't include a constant in the model.

# c. Write this model in terms of the backshift operator.
# (1 - B)^2*yt = (1 + theta1*B + theta2*B^2)*et

# d. Fit the model using R and examine the residuals. Is the model satisfactory?
wmurders_arima.0.2.2 <- Arima(wmurders, 
                              order = c(0, 2, 2))

checkresiduals(wmurders_arima.0.2.2)
# The residuals of the model can be thought of as white noise series. A little sorry that they aren't normally distributed. But it is satisfactory to get them.

# e. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
fc_wmurders_arima.0.2.2 <- forecast(
  wmurders_arima.0.2.2, h = 3
)

# forecasts by Arima function
fc_wmurders_arima.0.2.2$mean

# get forecasts by manual calculation
fc_wmurders_arima.0.2.2$model
# formula
# (1 - B)^2*yt = (1 - 1.0181*B + 0.1470*B^2)*et
# yt = 2yt-1 - yt-2 + et - 1.0181*et-1 + 0.1470*et-2
years <- length(wmurders)
e <- fc_wmurders_arima.0.2.2$residuals

fc1 <- 2*wmurders[years] - wmurders[years - 1] - 1.0181*e[years] + 0.1470*e[years - 1]

fc2 <- 2*fc1 - wmurders[years] + 0.1470*e[years]

fc3 <- 2*fc2 - fc1

# forecasts by manual calculation
c(fc1, fc2, fc3)
# the forecasts are almost similar to the ones got by Arima function.

# f. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
autoplot(fc_wmurders_arima.0.2.2)

# g. Does auto.arima give the same model you have chosen? If not, which model do you think is better?
fc_wmurders_autoarima <- forecast(
  auto.arima(wmurders), h = 3
)

# Without RMSE, all errors show that ARIMA(0, 2, 2) is better than ARIMA(1, 2, 1).
accuracy(fc_wmurders_arima.0.2.2)
accuracy(fc_wmurders_autoarima)

# try using auto.arima function with stepwise and approximation options false.
fc_wmurders_autoarima2 <- forecast(
  auto.arima(wmurders, stepwise = FALSE, approximation = FALSE), 
  h = 3
)
# It is ARIMA(0, 2, 3) model.

accuracy(fc_wmurders_autoarima2)
# In this case, some errors were better while others were worse. I'll check residuals and ACF, PACF plots.
ggtsdisplay(diff(wmurders, differences = 2))
# It looked like that the data are similar to ARIMA(0, 2, 2) rather than ARIMA(0, 2, 3).

checkresiduals(fc_wmurders_arima.0.2.2)
checkresiduals(fc_wmurders_autoarima2)
# almost similar residuals.

# Therefore I'll choose ARIMA(0, 2, 2).

```


8. Consider the total international visitors to Australia (in millions) for the period 1980-2015. (Data set austa.)

```{r echo=FALSE, message=FALSE, warning=FALSE, Question8}

# a. Use auto.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.
autoplot(austa)

fc_austa_autoarima <- forecast(
  auto.arima(austa), h = 10
)

fc_austa_autoarima$model
# ARIMA(0, 1, 1) with drift model was chosen.

checkresiduals(fc_austa_autoarima)
# The residuals are like white noise.

autoplot(fc_austa_autoarima)

# b. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part (a). Remove the MA term and plot again.
fc_austa_arima.0.1.1 <- forecast(
  Arima(austa, order = c(0, 1, 1)), h = 10
)

autoplot(fc_austa_arima.0.1.1)

fc_austa_arima.0.1.0 <- forecast(
  Arima(austa, order = c(0, 1, 0)), h = 10
)

autoplot(fc_austa_arima.0.1.0)
# the forecasts of both models are like the result of naive forecast. Increasing trend isn't reflected in the forecasts.

fc_austa_arima.0.1.1$upper - fc_austa_arima.0.1.0$upper
fc_austa_arima.0.1.0$lower - fc_austa_arima.0.1.1$lower
# But prediction interval of ARIMA(0, 1, 1) model was generally larger than the one of ARIMA(0, 1, 0) model. I think that it is because of one more error term in ARIMA(0, 1, 1) model.

# c. Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.
fc_austa_arima.2.1.3.drift <- forecast(
  Arima(austa, order = c(2, 1, 3), include.drift = TRUE),
  h = 10
)

autoplot(fc_austa_arima.2.1.3.drift)
# The forecasts are increasing, but the speed of the increase is decreasing.

drift_austa <- fc_austa_arima.2.1.3.drift$model$coef[6]
fc_austa_arima.2.1.3.nodrift <- fc_austa_arima.2.1.3.drift$mean - drift_austa*seq_len(10)

autoplot(fc_austa_arima.2.1.3.drift) +
  autolayer(fc_austa_arima.2.1.3.nodrift)
# Without drift constant, the forecasts are unlikely.

# d. Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.
fc_austa_arima.0.0.1.const <- forecast(
  Arima(
    austa, order = c(0, 0, 1), include.constant = TRUE
    ),
  h = 10
)

autoplot(fc_austa_arima.0.0.1.const)
# the forecasts are fastly decreased to the mean of the data history.

fc_austa_arima.0.0.0.const <- forecast(
  Arima(austa, order = c(0, 0, 0), include.constant = TRUE),
  h = 10
)

autoplot(fc_austa_arima.0.0.0.const)
# All of the forecasts are the mean of the data history. It is like the result of mean method.

# e. Plot forecasts from an ARIMA(0,2,1) model with no constant.
fc_austa_arima.0.2.1 <- forecast(
  Arima(austa, order = c(0, 2, 1)),
  h = 10
)

autoplot(fc_austa_arima.0.2.1)
# the forecasts show increasing trend. PI is being larger for the farther future forecast.

```


9. For the usgdp series:

```{r echo=FALSE, message=FALSE, warning=FALSE, Question9}

# a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(usgdp)
autoplot(BoxCox(usgdp, BoxCox.lambda(usgdp)))
# When I transformed the original data, I could get more linearly increasing line. Therefore I'm going to do Box-Cox transformation.
lambda_usgdp <- BoxCox.lambda(usgdp)

# b.fit a suitable ARIMA model to the transformed data using auto.arima();
usgdp_autoarima <- auto.arima(usgdp, 
                              lambda = lambda_usgdp)

autoplot(usgdp, series = "Data") +
  autolayer(usgdp_autoarima$fitted, series = "Fitted")
# It looked like the model fits well to the data.

usgdp_autoarima
#ARIMA(2, 1, 0) with drift model after Box-Cox transformation.

# c. try some other plausible models by experimenting with the orders chosen;
ndiffs(BoxCox(usgdp, lambda_usgdp))
# the data need 1 first differencing to be stationary.
ggtsdisplay(diff(BoxCox(usgdp, lambda_usgdp)))
# ACF plot shows sinusoidal decrease while PACF plot shows significant spikes at lag 1 and 12. I think that I can ignore the spike at lag 12 because the data are aggregated quarterly, not monthly. Therefore, I'll experiment with ARIMA(1, 1, 0) model.
usgdp_arima.1.1.0 <- Arima(
  usgdp, lambda = lambda_usgdp, order = c(1, 1, 0)
)

usgdp_arima.1.1.0

autoplot(usgdp, series = "Data") +
  autolayer(usgdp_arima.1.1.0$fitted, series = "Fitted")

# I'll also try ARIMA(1, 1, 0) with drift model.
usgdp_arima.1.1.0.drift <- Arima(
  usgdp, lambda = lambda_usgdp, order = c(1, 1, 0),
  include.drift = TRUE
)

usgdp_arima.1.1.0.drift

autoplot(usgdp, series = "Data") +
  autolayer(usgdp_arima.1.1.0.drift$fitted, series = "Fitted")
# It looked like that these models also fit well to the data.

# d. choose what you think is the best model and check the residual diagnostics;
accuracy(usgdp_autoarima)
accuracy(usgdp_arima.1.1.0)
accuracy(usgdp_arima.1.1.0.drift)
# Some errors show that ARIMA(2, 1, 0) with drift is the best model while others show that ARIMA(1, 1, 0) with drift is the best. Check the residuals of both cases.
checkresiduals(usgdp_autoarima)
checkresiduals(usgdp_arima.1.1.0.drift)
# In either case, the residuals are like white noise series and are not normally distributed.
# I'll choose the best model as ARIMA(2, 1, 0) with drift model. With the model, RMSE and MASE values were lower. And there wasn't significant spike at lag 2 in ACF plot of ARIMA(2, 1, 0) with drift model, even if it exists in ARIMA(1, 1, 0) with drift model.

# e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fc_usgdp_autoarima <- forecast(
  usgdp_autoarima
)

autoplot(fc_usgdp_autoarima)
# It looked like the forecasts are reasonable.

# f. compare the results with what you would obtain using ets() (with no transformation).
fc_usgdp_ets <- forecast(
  ets(usgdp)
)

autoplot(fc_usgdp_ets)
# It looked like these forecasts are more likely than the ones with ARIMA model. When trend is obvious, is ETS better than ARIMA model? I wonder about it.

```


10. Consider austourists, the quarterly number of international tourists to Australia for the period 1999-2010. (Data set austourists.)

```{r echo=FALSE, message=FALSE, warning=FALSE, Question10}

# a. Describe the time plot.
autoplot(austourists)
# the data have strong seasonality and increasing trend. Also the size of variations increased as the number increased.

# b. What can you learn from the ACF graph?
ggAcf(austourists)
# autocorrelations are slowly decreasing. And the values at the lags of multiple of 4 were big compared to the others.

# c. What can you learn from the PACF graph?
ggPacf(austourists)
# there are 5 significant spikes, and then no significant spikes thereafter (apart from one at lag 8, which are probably related with quarterly seasonality).

# d. Produce plots of the seasonally differenced data  (1 - B^4)Yt. What model do these graphs suggest?
ggtsdisplay(diff(austourists, lag = 4))
# the seasonally differenced data are looked like to need at least one more differencing to make it stationary.
# For the values at the lags of multiple of 4, there are just significant spikes at lag 4. It is same for ACF and PACF plots. The order of seasonal ARIMA model can be (1, 1, 0)[4] or (0, 1, 1)[4]. I'll choose (1, 1, 0)[4] order. 
# Disregarding the values at lag 4, autocorrelation values are looked like decreasing sinusoidally while partial autocorrelation values have spikes at lag 1 and 5. There aren't significant spikes at lag 2 and 3. 
# I think that the spike at lag 5 doesn't mean that there are still important information unused at 5th lagged values. It should've been positively significant because of the big negatively significant spike at lag 4 after seasonal differencing. 
# I can find this by drawing PACF plot after doing one more differencing.
ggtsdisplay(diff(diff(austourists, lag = 4)))
# Therefore I suggest ARIMA(1, 1, 0)(1, 1, 0)[4] model.

# e. Does auto.arima give the same model that you chose? If not, which model do you think is better?
fc_austourists_autoarima <- forecast(
  auto.arima(austourists)
)

fc_austourists_autoarima$model
# auto.arima gave ARIMA(1, 0, 0)(1, 1, 0)[4] model.

fc_austourists_arima.1.1.0.1.1.0.4 <- forecast(
  Arima(austourists, 
        order = c(1, 1, 0), 
        seasonal = c(1, 1, 0))
)

autoplot(fc_austourists_autoarima)
autoplot(fc_austourists_arima.1.1.0.1.1.0.4)
# ARIMA(1, 1, 0)(1, 1, 0)[4] shows more fastly increasing trend.

accuracy(fc_austourists_autoarima)
accuracy(fc_austourists_arima.1.1.0.1.1.0.4)
# ARIMA(1, 0, 0)(1, 1, 0)[4] with drift model was fitted better. Therefore I think that this model is better than ARIMA(1, 1, 0)(1, 1, 0)[4] model.

checkresiduals(fc_austourists_autoarima)
# The residuals are like white noise series.

# f. Write the model in terms of the backshift operator, then without using the backshift operator.
fc_austourists_autoarima$model
# (1 - phi1*B)(1 - phis1*B)(1 - B^4)(yt - c*t) = et
# c = drift*(1 - phi1)(1 - phis1)*m^D = 1.7793

# (1 - phi1*B - phis1*B + phi1*phis1*B^2)(1 - B^4)(yt - c*t) =
# (1 - phi1*B - phis1*B + phi1*phis1*B^2 - B^4 + phi1*B^5 + phis1*B^5 - phi1*phis1*B^6)(yt - c*t) = et
# yt = c + (phi1 + phis1)*yt-1 - phi1*phis1*yt-2 + yt-4 - (phi1 + phis1)*yt-5 + phi1*phis1*yt-6 + et
# yt = 1.7793 - 0.06*yt-1 + 0.2496*yt-2 + yt-4 + 0.06*yt-5 - 0.2496*yt-6 + et

```


11. Consider the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 - June 2013). (Data set usmelec.) In general there are two peaks per year: in mid-summer and mid-winter.

```{r echo=FALSE, message=FALSE, warning=FALSE, Question11}

# a. Examine the 12-month moving average of this series to see what kind of trend is involved.
usmelec_ma2x12 <- ma(usmelec, order = 12, centre = TRUE)

autoplot(usmelec, series = "Data") +
  autolayer(usmelec_ma2x12, series = "2X12-MA") +
  ylab(expression(paste("Electricity(x", 10^{9}, "KWh)"))) + 
  ggtitle("Monthly total net generation of electricity") +
  scale_color_discrete(breaks = c("Data", "2X12-MA"))
# Total net generation amount increased first but stoped increasing from about 2008.

# b. Do the data need transforming? If so, find a suitable transformation.
# The data show bigger variation for bigger amount. Therefore I think that Box-Cox transformation would be suitable for the data.
lambda_usmelec <- BoxCox.lambda(usmelec)

# c. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
# The data are non-stationary.
ndiffs(usmelec)
nsdiffs(usmelec)
# I need to do 1 seasonal differencing to make the data stationary. If seasonal differencing isn't enough to make them stationary, I need to do first differencing, too.

# d. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
ggtsdisplay(diff(
  BoxCox(usmelec, lambda_usmelec),
  lag = 12
  ))
# Definitely, I need to use first differencing, too.
ggtsdisplay(
  diff(
    diff(
      BoxCox(usmelec, lambda_usmelec),
      lag = 12
    )
  )
)
# I think that ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model might describe the data well. I'll try ARIMA(0, 1, 3)(0, 1, 1)[12] with Box-Cox transformation model, too.
usmelec_arima.0.1.2.0.1.1.12 <- Arima(
  usmelec,
  lambda = lambda_usmelec,
  order = c(0, 1, 2),
  seasonal = c(0, 1, 1)
)

usmelec_arima.0.1.3.0.1.1.12 <- Arima(
  usmelec,
  lambda = lambda_usmelec,
  order = c(0, 1, 3),
  seasonal = c(0, 1, 1)
)

usmelec_arima.0.1.2.0.1.1.12$aic
usmelec_arima.0.1.3.0.1.1.12$aic
# ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model was the best.

# e. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
usmelec_arima.0.1.2.0.1.1.12
#theta1 = -0.4317, theta2 = -0.2552, phis1 = -0.8536

checkresiduals(usmelec_arima.0.1.2.0.1.1.12)
# Ljung-Box test result shows that the residuals can be thought of as white noise. And they are normally distributed.

# I want to know what model was selected if I used auto.arima function. I'll try it.
usmelec_autoarima <- auto.arima(
  usmelec,
  lambda = lambda_usmelec
)

usmelec_autoarima
# The result is ARIMA(2, 1, 4)(0, 0, 2)[12] with drift after Box-Cox transformation model. AIC is -4722. But I can't compare the AIC value with what I got above, because the number of differencing was different(Differencing changes the way the likelihood is computed).

checkresiduals(usmelec_autoarima)
# And the residuals aren't like white noise. Therefore I'll choose ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model.

# f. Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from https://goo.gl/WZIItv to check the accuracy of your forecasts.
fc_usmelec_arima.0.1.2.0.1.1.12 <- forecast(
  usmelec_arima.0.1.2.0.1.1.12,
  h = 12*15
)

# Get the latest figures.
usmelec.new <- read.csv("MER_T07_02A.csv", sep = ",")

# need to do data munging before using the data.
# make new columns Year, Month using YYYYMM column.
usmelec.new[, "Year"] <- as.numeric(
  substr(usmelec.new[, "YYYYMM"], 1, 4)
  )
usmelec.new[, "Month"] <- as.numeric(
  substr(usmelec.new[, "YYYYMM"], 5, 6)
  )
# make usmelec.new only have Year, Month and Value columns with net generation total data.
usmelec.new <- subset(
  usmelec.new, 
  Description == "Electricity Net Generation Total, All Sectors", 
  select = c("Year", "Month", "Value")
  )
# remove data if month is 13. They are old yearly data.
usmelec.new <- subset(usmelec.new, Month != 13)
# change the Value column data type to number. And divide the numbers by 1000 because the unit of the values in usmelec.new are Million KWh, not Billion KWh.
usmelec.new[, "Value"] <- as.numeric(
  as.character(usmelec.new[, "Value"])
  )/1000
# as.numeric(usmelec.new[, "Value"]) yields wrong data. Need to recognize the letters as character first, and then change the type as number. 
head(usmelec.new)
tail(usmelec.new)
# first observation was taken in January, 1973. Final observation was taken in October, 2017.

# make ts time series using usmelec.new Value column data.
usmelec.new.ts <- ts(
  usmelec.new[, "Value"], 
  start = c(1973, 1), 
  frequency = 12
  )
tail(usmelec.new.ts)
# final observation was taken in October, 2017 as expected. 

# get accuracy for 4 years of forecast horizon.
usmelec.new.ts_next4years <- subset(
  usmelec.new.ts, 
  start = length(usmelec) + 1,
  end = length(usmelec) + 12*4
)
accuracy(
  fc_usmelec_arima.0.1.2.0.1.1.12, 
  usmelec.new.ts_next4years
  )

# plot the results
autoplot(fc_usmelec_arima.0.1.2.0.1.1.12, series = "Forecasts") +
  autolayer(usmelec.new.ts, series = "Real data") +
  scale_x_continuous(limits = c(2010, 2030)) +
  ggtitle("Forecast from ARIMA(0,1,2)(0,1,1)[12] with real data")
# Real data are really similar to the forecasts. Even when they were different, real data didn't get out of the prediction interval.

# g. How many years of forecasts do you think are sufficiently accurate to be usable?
# In usmelec data case, even 4 years of forecasts were sufficiently accurate to be usable. I think that it happened because the pattern in the data almost didn't change.

```


12. For the mcopper data:

```{r echo=FALSE, message=FALSE, warning=FALSE, Question12}

# a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(mcopper)
# they are monthly data but there isn't seasonality in them.
autoplot(BoxCox(mcopper, BoxCox.lambda(mcopper)))
# It looked like Box-Cox transformation makes the variations in the data evenly over time. Therefore I'm going to use the transformation.
lambda_mcopper <- BoxCox.lambda(mcopper)

# b. fit a suitable ARIMA model to the transformed data using auto.arima();
mcopper_autoarima <- auto.arima(
  mcopper,
  lambda = lambda_mcopper
)

mcopper_autoarima
# auto.arima yielded ARIMA(0, 1, 1) with Box-Cox transformation model. AICc was -86.08.

# c. try some other plausible models by experimenting with the orders chosen;
ndiffs(mcopper)
nsdiffs(mcopper)
# the data need 1 first differencing.
ggtsdisplay(diff(mcopper))
# It looked like autocorrelation values are sinusoidally decreasing. I'll choose ARIMA model's order as (1, 1, 0) and (5, 1, 0).

mcopper_arima.1.1.0 <- Arima(
  mcopper, order = c(1, 1, 0), lambda = lambda_mcopper
)

mcopper_arima.1.1.0
# AICc was -75.64.

mcopper_arima.5.1.0 <- Arima(
  mcopper, order = c(5, 1, 0), lambda = lambda_mcopper
)

mcopper_arima.5.1.0
# AICc was -78.48.

# I'll try auto.arima function without approximation and stepwise options.
mcopper_autoarima2 <- auto.arima(
  mcopper, lambda = lambda_mcopper,
  approximation = FALSE, stepwise = FALSE
)

mcopper_autoarima2
# the result model is the same as when I didn't use the options.

# d. choose what you think is the best model and check the residual diagnostics;
# When I compared AICc values, I got the smallest when I used auto.arima function. (I could've used AICc in comparing because the differencing was the same for all models I chose.) Best model is ARIMA(0, 1, 1) with Box-Cox transformation.
checkresiduals(mcopper_autoarima)
# The residuals are like white noise. I'll select the model.

# e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fc_mcopper_autoarima <- forecast(
  mcopper_autoarima
)

autoplot(fc_mcopper_autoarima)
# The forecasts aren't reasonable.

# I'll try other models I made.
fc_mcopper_arima.1.1.0 <- forecast(
  mcopper_arima.1.1.0
)

autoplot(fc_mcopper_arima.1.1.0)
# got almost same result

fc_mcopper_arima.5.1.0 <- forecast(
  mcopper_arima.5.1.0
)

autoplot(fc_mcopper_arima.5.1.0)
# got almost same result, too

# f. compare the results with what you would obtain using ets() (with no transformation).
fc_mcopper_ets <- forecast(
  ets(mcopper)
)

autoplot(fc_mcopper_ets)
# These forecasts are more reasonable than what I got above.

```


13. Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.

```{r echo=FALSE, message=FALSE, warning=FALSE, Question13}

# a. Do the data need transforming? If so, find a suitable transformation.
# I'll analyze and forecast qauselec data.
autoplot(qauselec)
# The data need Box-Cox transformation to make the variations evenly over time.
lambda_qauselec <- BoxCox.lambda(qauselec)

# b. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
# The data have strong seasonality and increasing trend. It means that the data aren't stationary.
nsdiffs(qauselec)
ndiffs(qauselec)
# The data need 1 seasonal differencing.
kpss.test(diff(qauselec, lag = 4))
# The data don't need first differencing. But I'll try with first differencing, too.

# c. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
ggtsdisplay(diff(
  BoxCox(qauselec, lambda_qauselec), lag = 4
  ))
ggtsdisplay(diff(diff(
  BoxCox(qauselec, lambda_qauselec), lag = 4
  )))

# suggest models:
# ARIMA(0, 0, 1)(0, 1, 1)[4]
qauselec_arima0.0.1.0.1.1.4 <- Arima(
  qauselec, lambda = lambda_qauselec,
  order = c(0, 0, 1), seasonal = c(0, 1, 1)
)

qauselec_arima0.0.1.0.1.1.4
# AIC = -97.01

# ARIMA(0, 1, 1)(0, 1, 1)[4]
qauselec_arima0.1.1.0.1.1.4 <- Arima(
  qauselec, lambda = lambda_qauselec,
  order = c(0, 1, 1), seasonal = c(0, 1, 1)
)

qauselec_arima0.1.1.0.1.1.4
#AIC = -292.59

# ARIMA(0, 1, 1)(0, 1, 2)[4]
qauselec_arima0.1.1.0.1.2.4 <- Arima(
  qauselec, lambda = lambda_qauselec,
  order = c(0, 1, 1), seasonal = c(0, 1, 2)
)

qauselec_arima0.1.1.0.1.2.4
#AIC = -292.73

# try using auto.arima function
qauselec_autoarima <- auto.arima(
  qauselec, lambda = lambda_qauselec
)

qauselec_autoarima
# AIC = -300.84

# According to AIC values, ARIMA(1, 1, 1)(1, 1, 2)[4] with Box-Cox transformation is the best model. But above 2 models and the model made by auto.arima function used different number of differencings.

# d. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
# phi1 = 0.2523, theta1 = -0.6905, phis1 = 0.8878, thetas1 = -1.6954, thetas2 = -0.7641 
checkresiduals(qauselec_autoarima)
# The residuals aren't like white noise.

# try using other models.
checkresiduals(qauselec_arima0.0.1.0.1.1.4)
checkresiduals(qauselec_arima0.1.1.0.1.1.4)
checkresiduals(qauselec_arima0.1.1.0.1.2.4)
# The residuals don't resemble white noise regardless of the model. Therefore I'm going to use the best model.

# e. Forecast the next 24 months of data using your preferred model.
fc_qauselec_autoarima <- forecast(
  qauselec_autoarima, h = 8
)

autoplot(fc_qauselec_autoarima)
# The forecasts are reasonable.

# f. Compare the forecasts obtained using ets().
fc_qauselec_ets <- forecast(
  ets(qauselec), h = 8
)

autoplot(fc_qauselec_ets)
# These forecasts also are reasonable.

```


14. For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. The stlf() function will make the calculations easy (with method="arima"). Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?

```{r echo=FALSE, message=FALSE, warning=FALSE, Question14}

fc_qauselec_stlf <- stlf(
  qauselec, lambda = BoxCox.lambda(qauselec),
  s.window = 5, robust = TRUE, method = "arima",
  h = 8
)

autoplot(fc_qauselec_stlf) +
  scale_x_continuous(limits = c(2005, 2012)) +
  scale_y_continuous(limits = c(50, 70))

autoplot(fc_qauselec_ets) +
  scale_x_continuous(limits = c(2005, 2012)) +
  scale_y_continuous(limits = c(50, 70))

autoplot(fc_qauselec_autoarima) +
  scale_x_continuous(limits = c(2005, 2012)) +
  scale_y_continuous(limits = c(50, 70))
# I don't know which forecasts are best. Forecasts from STL + ARIMA(0, 1, 1) with drift model yielded highest forecasts. It looked like it followed latest trend a lot. ETS(M, A, M) model yielded broadest PI. ARIMA(1, 1, 1)(1, 1, 2)[4] model yielded least varianced forecasts. 
# For qauselec data, the choice of model didn't affect forecasts much because they already have strong seasonality and almost non-changing trend.

```

15. For your retail time series (Exercise 5 above):

```{r echo=FALSE, message=FALSE, warning=FALSE, Question15}

# a. develop an appropriate seasonal ARIMA model;
fc_retail_autoarima <- forecast(
  auto.arima(retail.ts),
  h = 36
)

autoplot(fc_retail_autoarima)
# ARIMA(1, 0, 2)(0, 1, 1)[12] with drift model was chosen.

# b. compare the forecasts with those you obtained in earlier chapters;
# In chapter 3, I used seasonal naive method to forecast.
fc_retail_snaive <- snaive(retail.ts, h = 36)

autoplot(fc_retail_snaive)

# In chapter 7, I thought that Holt-Winters' multiplicative method was best among ets models.
fc_retail_ets <- forecast(
  ets(retail.ts, lambda = BoxCox.lambda(retail.ts)), 
  h = 36
)

autoplot(fc_retail_ets)

# c. Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?
# Get the latest figures.
retail.new <- read.xlsx("8501011.xlsx", 
                        sheetName = "Data1", 
                        startRow = 10)

retail.new.ts <- ts(retail.new[, "A3349873A"],
                    start = c(1982, 4),
                    frequency = 12)

retail.new.test <- subset(
  retail.new.ts,
  start = length(retail.ts) + 1
  )

# draw plots for the data in forecast horizon.
autoplot(fc_retail_autoarima$mean, series = "ARIMA") +
  geom_line(size = 1, color = "red") +
  autolayer(fc_retail_ets$mean, series = "Holt-Winters'", size = 1) +
  autolayer(fc_retail_snaive$mean, series = "Seasonal Naive", size = 1) +
  autolayer(retail.new.test, series = "Real data", size = 1) +
  ggtitle("Turnover of other retailings not elsewhere classified",
          subtitle = "- From New South Wales in Australia") +
  ylab("Trade amount($ Millions)") +
  scale_color_discrete(breaks = c("Real data", "ARIMA", "Holt-Winters'", "Seasonal Naive"))
# plots shows that the forecasts from Holt-Winters' model were the most accurate.

# get accuracy for each models' forecasts.
accuracy(fc_retail_autoarima, retail.new.test)
accuracy(fc_retail_ets, retail.new.test)
accuracy(fc_retail_snaive, retail.new.test)
# Without ACF1, all errors show that the Holt-Winters' model was best. The next one was ARIMA model. Seasonal naive model showed the poorest forecasting accuracy.

```


16. 
```{r echo=FALSE, message=FALSE, warning=FALSE, Question16}

# a. Produce a time plot of the sheep population of England and Wales from 1867-1939 (data set sheep).
autoplot(sheep)
# There is decreasing trend, but no definite seasonality.

# b. Assume you decide to fit the following model:
# yt = yt-1 + phi1(yt-1 - yt-2) + phi2(yt-2 - yt-3) + phi3(yt-3 - yt-4) + et
# where et is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?

# (yt - yt-1) - phi1(yt-1 - yt-2) - phi2(yt-2 - yt-3) - phi3(yt-3 - yt-4) = et
# (1 - B)yt - phi1*B(1- B)yt - phi2*B^2(1- B)yt - phi3*B^3(1- B)yt = et
# (1 - phi1*B - phi2*B^2 - phi3*B^3)(1 - B)yt = et
# It is ARIMA(3, 1, 0) model.

# c. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
ggtsdisplay(diff(sheep))
# ACF plot shows sinusoidally decreasing autocorrelation values while PACF plot shows significant spikes at lag 1 to 3, but no beyond lag 3. Therefore ARIMA(3, 1, 0) is appropriate.

# d. The last five values of the series are given below:
# Year	            1935	1936	1937	1938	1939
# Millions of sheep	1648	1665	1627	1791	1797
# The estimated parameters are phi1 = 0.42, phi2 = -0.20, and phi3 = -0.30. Without using the forecast function, calculate forecasts for the next three years (1940-1942).

sheep.1940 = 1797 + 0.42*(1797 - 1791) -0.20*(1791 - 1627) - 0.30*(1627 - 1665)

sheep.1941 = sheep.1940 + 0.42*(sheep.1940 - 1797) -0.20*(1797 - 1791) - 0.30*(1791 - 1627)

sheep.1942 = sheep.1941 + 0.42*(sheep.1941 - sheep.1940) -0.20*(sheep.1940 - 1797) - 0.30*(1797 - 1791)

c(sheep.1940, sheep.1941, sheep.1942)

# e. Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?
fc_sheep_arima.3.1.0 <- forecast(
  Arima(sheep, order = c(3, 1, 0)),
  h = 3
)

fc_sheep_arima.3.1.0$mean
# calculated forecasts were a little bigger than the forecasts from Arima function. And the differences between the forecasts at the same time became bigger as time goes on. I think that it happened because of the differences of the coefficients. 
# Small differences in the coefficients made the difference between the first forecasts. And then the forecast values were used to calculate the next time point's forecasts. When the next time point's forecasts of Arima function were calculated, the difference became bigger. It looked like such situation repeated.
ar1 <- fc_sheep_arima.3.1.0$model$coef[1]
ar2 <- fc_sheep_arima.3.1.0$model$coef[2]
ar3 <- fc_sheep_arima.3.1.0$model$coef[3]

sheep.1940.new = 1797 + ar1*(1797 - 1791) + ar2*(1791 - 1627) + ar3*(1627 - 1665)

sheep.1941.new = sheep.1940.new + ar1*(sheep.1940.new - 1797) + ar2*(1797 - 1791) + ar3*(1791 - 1627)

sheep.1942.new = sheep.1941.new + ar1*(sheep.1941.new - sheep.1940.new) + ar2*(sheep.1940.new - 1797) + ar3*(1797 - 1791)

c(sheep.1940.new, sheep.1941.new, sheep.1942.new)
# above calculation confirms what I said about the differences.

```


17.
```{r echo=FALSE, message=FALSE, warning=FALSE, Question17}

# a. Plot the annual bituminous coal production in the United States from 1920 to 1968 (data set bicoal).
autoplot(bicoal)
# It looked like there isn't any particular trend or seasonality.

# b. You decide to fit the following model to the series:
# yt = c + phi1*yt-1 + phi2*yt-2 + phi3*yt-3 + phi4*yt-4 + et
# where yt is the coal production in year t and et is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?

# (1 - phi1*B - phi2*B^2 - phi3*B^3 - phi4*B^4)*yt = c + et
# if mu is the mean of yt,
# c = mu*(1 - phi1*B - phi2*B^2 - phi3*B^3 - phi4*B^4)
# This model is ARIMA(4, 0, 0) or AR(4).

# c. Explain why this model was chosen using the ACF and PACF.
ggAcf(bicoal, lag.max = 36)
ggPacf(bicoal, lag.max = 36)
# ACF plot shows sinusoidally decreasing autocorrelation values. PACF plot shows significant spikes at lag 1 and 4, but none beyond lag 4. Therefore AR(4) model is the appropriate choice.

# d. The last five values of the series are given below.
# Year	           1964	1965	1966	1967	1968
# Millions of tons	467	 512	 534	 552 	 545
# The estimated parameters are c = 162.00, phi1 = 0.83, phi2 = -0.34, phi3 = 0.55 and phi4 = -0.38. Without using the forecast function, calculate forecasts for the next three years (1969-1971).
c = 162.00
phi1 = 0.83 
phi2 = -0.34
phi3 = 0.55
phi4 = -0.38

bicoal.1969 <- c + phi1*545 + phi2*552 + phi3*534 + phi4*512

bicoal.1970 <- c + phi1*bicoal.1969 + phi2*545 + phi3*552 + phi4*534

bicoal.1971 <- c + phi1*bicoal.1970 + phi2*bicoal.1969 + phi3*545 + phi4*552

c(bicoal.1969, bicoal.1970, bicoal.1971)

# e. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
fc_bicoal_ar4 <- forecast(ar(bicoal, 4), h = 3)

fc_bicoal_ar4$mean
# The forecasts from ar function were a little bigger than the calculated forecasts. It also happened because of the small differences of coefficients.
phi1 <- fc_bicoal_ar4$model$ar[1]
phi2 <- fc_bicoal_ar4$model$ar[2]
phi3 <- fc_bicoal_ar4$model$ar[3]
phi4 <- fc_bicoal_ar4$model$ar[4]
c <- fc_bicoal_ar4$model$x.mean*(1 - phi1 - phi2 - phi3 - phi4)

bicoal.1969.new <- c + phi1*545 + phi2*552 + phi3*534 + phi4*512

bicoal.1970.new <- c + phi1*bicoal.1969.new + phi2*545 + phi3*552 + phi4*534

bicoal.1971.new <- c + phi1*bicoal.1970.new + phi2*bicoal.1969.new + phi3*545 + phi4*552

c(bicoal.1969.new, bicoal.1970.new, bicoal.1971.new)
# These calculations and results confirm the causation of the differences in forecasts.

```


18. 
```{r echo=FALSE, message=FALSE, warning=FALSE, Question18}

# a. Install the rdatamarket package in R using
#  install.packages("rdatamarket")

# b. Select a time series from http://datamarket.com/data/list/?q=pricing:free. Then copy its short URL and import the data using
#  x <- ts(rdatamarket::dmseries("shorturl")[,1], start=??, frequency=??)
# (Replace ?? with the appropriate values.)
oil.NA <- ts(dmseries("http://data.is/1qk5Uvf")[, "Total.North.America"], start=1965, frequency=1)

str(oil.NA)
head(oil.NA)

# c. Plot graphs of the data, and try to identify an appropriate ARIMA model.
autoplot(oil.NA)
ndiffs(oil.NA)
ggtsdisplay(diff(oil.NA))
# I think that ARIMA(1, 1, 0) model will fit well to the data.

# d. Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?
oil.NA_arima.1.1.0 <- Arima(
  oil.NA, order = c(1, 1, 0)
)

checkresiduals(oil.NA_arima.1.1.0)
# The residuals are like white noise.

# e. Use your chosen ARIMA model to forecast the next four years.
fc_oil.NA_arima.1.1.0 <- forecast(
  oil.NA_arima.1.1.0, h = 4
)

autoplot(fc_oil.NA_arima.1.1.0)
# The forecast values of oil production are increasing. But the increase speed is damping.

# f. Now try to identify an appropriate ETS model.
oil.NA_ets <- ets(oil.NA)

oil.NA_ets
# chosen model is ETS(A, A, N).

# g. Do residual diagnostic checking of your ETS model. Are the residuals white noise?
checkresiduals(oil.NA_ets)
# The residuals are like white noise.

# h. Use your chosen ETS model to forecast the next four years.
fc_oil.NA_ets <- forecast(
  oil.NA_ets, h = 4
)

autoplot(fc_oil.NA_ets)
# The forecast values of oil production are increasing. And the increase speed doesn't almost change.

# i. Which of the two models do you prefer?
# I prefer ARIMA model. Because I think that it is unlikely that the oil production amounts increase a lot for more than 10 years in a row. It is likely that the increase of the production amounts will be decreased.
# But the real production amount was about 931 Million Ton in 2015. The forecasts from ETS model were more similar to the real production amounts for about 3 years of forecast horizon.
# http://tools.bp.com/energy-charting-tool.aspx?_ga=2.211706901.1002175336.1518264329-1800202796.1518264329#/st/oil/dt/production/unit/KBD/region/NOA/SCA/EU/MIE/AFR/AP/view/area/

```


### Question 1 isn't related with coding that I didn't include in this file.