DATA624 Spring 2021 HW8 - ARIMA

8.11.1

a.

Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers. Explain the differences among these figures.

Do they all indicate that the data are white noise?

Starting with the definition of “white noise” which is “a time series that shows no autocorrelation”.

Yes. All three random time series are indeed white noise. Since they’re each comprised of random numbers, there should not be any correlation. To replicate the problem, the three time series for 36, 360, and 1000 random entries are plotted below.

set.seed(3)
random36 <- ts(rnorm(36))
random360 <- ts(rnorm(360))
random1000 <- ts(rnorm(1000))

A series is white noise if its mean equals zero. A check of the mean for each series shows this.

cat(paste0("The mean of random36 is ", round(mean(random36),3), ". \nThe mean of random360 is ", round(mean(random360),3), ". \nThe mean of random1000 is ", round(mean(random1000),3), ". "))

## The mean of random36 is -0.095. 
## The mean of random360 is 0.046. 
## The mean of random1000 is -0.029.

autoplot(random36)

autoplot(random360)

autoplot(random1000)

Plotting the Acf for each series shows that the correlation bars never exceed the individual bounds.

x1 <- ggAcf(random36) + scale_x_continuous(limits = c(5, 20)) + 
                        scale_y_continuous(limits = c(-1, 1))

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

x2 <- ggAcf(random360) + scale_y_continuous(limits = c(-1, 1))
x3 <- ggAcf(random1000) + scale_y_continuous(limits = c(-1, 1))
grid.arrange(x1,x2,x3,ncol=3)

Ljung-Box Test

Considers the first h autocorrelation values together. A significant test (a small p value < 0.05) indicates the data are probably not white noise. IOW, we fail to reject the null hypothesis if the p value from the Ljung-Box test is greater than 0.05. We see that each distribution’s p-value is above that threshold.

Box.test(random36, lag=20, fitdf = 0, type="Lj")

## 
##  Box-Ljung test
## 
## data:  random36
## X-squared = 16.438, df = 20, p-value = 0.6891

Box.test(random360, lag=20, fitdf = 0, type="Lj")

## 
##  Box-Ljung test
## 
## data:  random360
## X-squared = 25.853, df = 20, p-value = 0.1707

Box.test(random1000, lag=20, fitdf = 0, type="Lj")

## 
##  Box-Ljung test
## 
## data:  random1000
## X-squared = 19.144, df = 20, p-value = 0.5125

White noise is a purely random time series. The Box Ljung test confirms whether you have a white noise time series.

As the author states:

we expect each autocorrelation to be close to zero. Of course, they will not be exactly equal to zero as there is some random variation. For a white noise series, we expect 95% of the spikes in the ACF to lie within \(\pm2/\sqrt{T}\) where \(T\) is the length of the time series.

b.

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

The critical boundaries lines are defined by the formula: \(\pm1.96/\sqrt{N}\) where \(N\) is the size of the sample. In this problem, the size of the samples grow, therefore the critical bounds will change. The autocorrelations are different with each sample size since as the population grows, the chances that there will be correlations however remote declines.

8.11.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.

The author states: A stationary time series is one whose properties do not depend on the time at which the series is observed. Thus, time series with trends, or with seasonality, are not stationary - the trend and seasonality will affect the value of the time series at different times.

From the plot below, the closing price shows a slow rise, followed by a plateau, a brisk downward trend, and a bounce off the lows over the time period.

autoplot(ibmclose)

The ACF plot shows all previous lags to be significant and decaying slightly given the length of the lag. The PACF shows only one significant spike.

ggAcf(ibmclose)

ggPacf(ibmclose)

The ndiffs function shows that we should take one difference.

ndiffs(ibmclose)

## [1] 1

The auto arima function shows that the time series should be an ARIMA(1,1,0)

fit <- auto.arima(ibmclose, seasonal=FALSE,
  stepwise=FALSE, approximation=FALSE)
summary(fit)

## Series: ibmclose 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.0869
## s.e.  0.0519
## 
## sigma^2 estimated as 52.36:  log likelihood=-1249.97
## AIC=2503.94   AICc=2503.98   BIC=2511.76
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set -0.252425 7.216455 5.205757 -0.07643735 1.198182 0.9993315
##                       ACF1
## Training set -0.0005345879

When one difference is applied to the ibmclose time series, we see from the plot that there is no trend nor seasonality.

autoplot(diff(ibmclose))

The ACF plot shows spikes at three different lag periods which means that there’s some correlations in the data even after differencing. However, it is reasonably stationary at this point.

ggAcf(diff(ibmclose))

ggPacf(diff(ibmclose))

Finally, because the p-value for the differenced ibmclose time series from the Box test is lower than 0.05. That means that the residuals are independent

Box.test(diff(ibmclose), lag=25,  type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  diff(ibmclose)
## X-squared = 42.308, df = 25, p-value = 0.01665

8.11.3

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

Custom function - create_plots and BoxCox_Stationary

This custom function, “create_plot”, creates plots.

#A function that plots the time series, the Acf, and the PAcf

create_plots <- function(ts, title){
   grid.arrange( 
    autoplot(ts) +
      ggtitle(title) ,
    grid.arrange(
      ggAcf(ts) +
        ggtitle(paste0(title, " ACF")), 
      ggPacf(ts) + 
        ggtitle(paste0(title, " PACF")), ncol = 2), nrow = 2)
  
  
}

Custom function - BoxCox_Stationary

This function applies a Box Cox transformation on a time series, finds and reports on the number of differences, and then applies the differences.

#The BoxCox_Stationary function takes in a time series and title for the plots 
BoxCox_Stationary <- function(ts, title){
#Use the BoxCox lambda function to find the best lambda and print it out.
  lambda <- BoxCox.lambda(ts)
#Transform the series with BoxCox function and the best lambda
  ts_transform <- BoxCox(ts, lambda)
#Find the number of differences to take with 
  diffs <- ndiffs(ts_transform)
#Print out the best lambda and the number of differences
  print(paste0("The best lambda for ", title, " is: ", round(lambda,2), " and the number of difference is: ", diffs))
#For loop differentiates the time series using the diffs.
  for (i in diffs){
    ts_transform <- diff(ts_transform)
    }
#Return the Box Cox transformed time series
  return (ts_transform)
  }

a. usnetelec

The plots below show that the usnetelec time series is not stationary. From the plots below, we see a rising trend over time. Additionally, for a stationary time series, the ACF will drop to zero relatively quickly, while the ACF of non-stationary data decreases slowly.

create_plots(usnetelec, "usnetelec")

To confirm that the series is stationary, I applied the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test:

“The Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test figures out if a time series is stationary around a mean or linear trend, or is non-stationary due to a unit root. A stationary time series is one where statistical properties - like the mean and variance - are constant over time. The null hypothesis for the test is that the data is stationary. The alternate hypothesis for the test is that the data is not stationary”

  https://www.statisticshowto.com/kpss-test/

ur.kpss(usnetelec) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 1.464 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Here, we see that the test statistic is 1.464 which means that we reject the null hypothesis that usnetelec is stationary.

The best lambda for usnetelec is: 0.52 and the number of difference is: 2

usnetelec_trans <- BoxCox_Stationary(usnetelec, "usnetelec")

## [1] "The best lambda for usnetelec is: 0.52 and the number of difference is: 2"

create_plots(usnetelec_trans, "usnetelec transformed")

Given that the test statistic is 0.4315 which is lower than the 5% significance level, we fail to reject the null hypothesis.

ur.kpss(usnetelec_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.4315 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

b. usgdp

create_plots(usgdp, "usgdp")

The best lambda for usgdp is: 0.37 and the number of difference is: 1

usgdp_trans <- BoxCox_Stationary(usgdp, "usgdp")

## [1] "The best lambda for usgdp is: 0.37 and the number of difference is: 1"

create_plots(usgdp_trans, "usgdp transformed")

ur.kpss(usgdp_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.2013 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The value of the test-statistic: 0.2013 is lower than any of the critical values so it’s a failure to reject the null hypothesis which is that the differenced time series is stationary.

c. mcopper

create_plots(mcopper, "mcopper")

mcopper_trans <- BoxCox_Stationary(mcopper, "mcopper")

## [1] "The best lambda for mcopper is: 0.19 and the number of difference is: 1"

create_plots(mcopper_trans, "mcopper transformed")

The value of the test-statistic: 0.0573 is lower than any of the critical values so it’s a failure to reject the null hypothesis which is that the differenced time series is stationary.

ur.kpss(mcopper_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 6 lags. 
## 
## Value of test-statistic is: 0.0573 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

d. enplanements

create_plots(enplanements, "enplanements")

enplanements_trans <- BoxCox_Stationary(enplanements, "mcopper")

## [1] "The best lambda for mcopper is: -0.23 and the number of difference is: 1"

create_plots(enplanements_trans, "enplanements transformed")

The value of the test-statistic: 0.0151 is lower than any of the critical values so it’s a failure to reject the null hypothesis which is that the differenced time series is stationary.

ur.kpss(enplanements_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0151 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

e. visitors

create_plots(visitors, "visitors")

visitors_trans <- BoxCox_Stationary(visitors, "visitors")

## [1] "The best lambda for visitors is: 0.28 and the number of difference is: 1"

create_plots(visitors_trans, "visitors transformed")

The value of the test-statistic: 0.0519 is lower than any of the critical values so it’s a failure to reject the null hypothesis which is that the differenced time series is stationary.

ur.kpss(visitors_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.0519 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

8.11.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.

The best lambda for liquor_sales is: -0.04 and the number of difference is: 1

retaildata <- read_excel("retail.xlsx", skip=1)
liquor_sales <- ts(retaildata[,"A3349414R"],
  frequency=12, start=c(1982,4))

create_plots(liquor_sales, "liquor_sales")

liquor_sales_trans <- BoxCox_Stationary(liquor_sales, "liquor_sales")

## [1] "The best lambda for liquor_sales is: -0.04 and the number of difference is: 1"

create_plots(liquor_sales_trans, "liquor_sales_transformed")

ur.kpss(liquor_sales_trans) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0286 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

auto.arima(liquor_sales)

## Series: liquor_sales 
## ARIMA(1,1,2)(0,1,2)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sma1    sma2
##       0.8917  -1.3669  0.3773  -0.3661  0.0926
## s.e.  0.0404   0.0686  0.0646   0.0552  0.0523
## 
## sigma^2 estimated as 25.41:  log likelihood=-1116.42
## AIC=2244.85   AICc=2245.08   BIC=2268.29

8.11.6

a.

Use the following R code to generate data from an AR(1) model with \(\phi_1=0.6\) and \(\sigma^2=1\) The process starts with \(y_1=0\)

To automate the code, the custom function, generate_AR1_model, takes in a value for \(\phi_1\) and produces time series

generate_AR1_model <- function(phi){
    set.seed(22)
    y <- ts(numeric(100))
    e <- rnorm(100)
    e[1] <- 0
    constant <- 0 
    for(i in 2:100)
      y[i] <- constant + e[i] + phi*y[i-1] 
    return(y)
}

b.

Produce a time plot for the series. How does the plot change as you change \(\phi_1\)

As the size of phi moves from small to large, the variability of the time series increases.

autoplot(generate_AR1_model(0.2), series="phi=0.6") +
  autolayer(generate_AR1_model(0.4) , series = "phi=0.4") +
  autolayer(generate_AR1_model(0.6) ,series = "phi=0.2") +
  autolayer(generate_AR1_model(0.7) ,series = "phi=0.7") +
  autolayer(generate_AR1_model(0.8) ,series = "phi=0.8") +
  xlab("Time") +
  ylab("Different Levels of Phi") +
  guides(colour=guide_legend(title="Different Levels of Phi"), 
         fill=guide_legend(title="AR"))+
      scale_color_manual(values = c("deepskyblue", "deepskyblue1","deepskyblue2",  "deepskyblue3", "deepskyblue4"))

c. 

Write your own code to generate data from an MA(1) model with \(\theta_1=0.6\) and \(\sigma^2=1\)

To automate the code, the custom function, generate_MA1_model, takes in a value for \(\theta_1\) and produces time series

generate_MA1_model <- function(theta){
    set.seed(22)
    y <- ts(numeric(100))
    e <- rnorm(100)
    e[1] <- 0
    constant <- 0
    for(i in 2:100)
      y[i] <-  constant + e[i] + theta*e[i-1] 
    return(y)
}

d.

Produce a time plot for the series. How does the plot change as you change \(\theta_1\)

\(\theta_1\) does little to nothing to change the variation of the series.

autoplot(generate_MA1_model(0.2), series="theta=0.6") +
  autolayer(generate_MA1_model(0.4) , series = "theta=0.4") +
  autolayer(generate_MA1_model(0.6) ,series = "theta=0.2") +
  autolayer(generate_MA1_model(0.7) ,series = "theta=0.7") +
  autolayer(generate_MA1_model(0.8) ,series = "theta=0.8") +
  xlab("Time") +
  ylab("Different Levels of theta") +
  guides(colour=guide_legend(title="Different Levels of theta"), 
         fill=guide_legend(title="AR"))+
      scale_color_manual(values = c("gold", "gold1","gold2",  "gold3", "black"))

e.

Generate data from an ARMA(1,1) model with \(\phi_1=0.6, \theta_1=0.6,\) and $ ^2=1$

generate_ARMA1_1_model <- function(phi, theta){
    set.seed(22)
    y <- ts(numeric(100))
    e <- rnorm(100)
    e[1] <- 0
    constant <- 0
    for(i in 2:100)
      y[i] <-  constant + e[i] + phi*y[i-1] + theta*e[i-1] 
    return(y)
}

y = generate_ARMA1_1_model(.6, .6)
autoplot(y)

f. Generate data from an AR(2) model with \(\phi_1 = -0.8\) \(\phi_2 = -0.8\) \(\sigma^2=1\)

generate_AR2_model <- function(phi_1, phi_2){
    set.seed(22)
    y <- ts(numeric(100))
    e <- rnorm(100)
    constant <- 0 
    for(i in 3:100)
      y[i] <-  constant + e[i] + phi_1*y[i-1] + phi_2*y[i-2]  
    return(y)
}

y <- generate_AR2_model(-0.8, .3)
autoplot(y)

g. Graph the latter two series and compare them.

There is greater variation in the AR 2 particularly towards the end of the series. The model whereas the ARMA(1,1) model has a lot more

autoplot(generate_ARMA1_1_model(0.6, 0.6), series="ARMA(1, 1)") +
  autolayer(generate_AR2_model(-0.8, 0.3) , series = "AR(2)") +
  xlab("Time") +
  ylab("Different Levels of theta") +
  guides(colour=guide_legend(title="Different Levels of theta"), 
         fill=guide_legend(title="ARMA"))+
      scale_color_manual(values = c( "orange",  "blue"))

8.11.7

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

a.

By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.

Following the steps that the author lays Out:

1. Plot the data and identify any unusual observations.

We see from the plots below that the time series is not stationary. There is a sharp upward trend from 1960 through 1975, a plateau from 1975 to 1995, and then a decline from 1995 past 2000. The ACF shows us that there was significant auto correlation between the time periods. We can clearly see that this time series is not stationary.

create_plots(wmurders, "Murders of Women")

2. If necessary, transform the data (using a Box-Cox transformation) to stabilise the variance.

Using the Box-Cox tranformation the variance of the time series was reduced to zero.

var(wmurders)

## [1] 0.5604752

#Use the BoxCox lambda function to find the best lambda and print it out.
  lambda <- BoxCox.lambda(wmurders)
#Transform the series with BoxCox function and the best lambda
  wmurders_transformed <- BoxCox(wmurders, lambda)

var(wmurders_transformed)

## [1] 0.04259765

3. If the data are non-stationary, take first differences of the data until the data are stationary.

The time series only becomes stationary after taking 2 difference as the KPSS Unit Root Tests show.

Taking the first difference, wmurders_transformed_diff

  wmurders_transformed_diff <- diff(wmurders_transformed)

create_plots(wmurders_transformed_diff, "Murders of Women -Diff")

ur.kpss(wmurders_transformed_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.5466 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Taking the second difference, wmurders_transformed_ddiff:

wmurders_transformed_ddiff <- diff(wmurders_transformed_diff)

ur.kpss(wmurders_transformed_ddiff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.0532 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

4. Examine the ACF/PACF to determine the correct ARIMA(p,d,q) structure.

Looking at the ACF and PACF charts, it looks like the best model would be ARIMA(1,2,1).

create_plots(wmurders_transformed_ddiff, "Murders of Women 2 Diff")

5.Try your chosen model(s), and use the AICc to search for a better model.

The AICc results confirm that the best values for ARIMA(p,d,q) would be ARIMA(1,2,1).

Arima(wmurders, order=c(0,2,0))

## Series: wmurders 
## ARIMA(0,2,0) 
## 
## sigma^2 estimated as 0.1007:  log likelihood=-14.38
## AIC=30.76   AICc=30.84   BIC=32.73

Arima(wmurders, order=c(1,2,0))

## Series: wmurders 
## ARIMA(1,2,0) 
## 
## Coefficients:
##           ar1
##       -0.6719
## s.e.   0.0981
## 
## sigma^2 estimated as 0.05471:  log likelihood=2
## AIC=0   AICc=0.24   BIC=3.94

Arima(wmurders, order=c(1,2,1))

## Series: wmurders 
## ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.2434  -0.8261
## s.e.   0.1553   0.1143
## 
## sigma^2 estimated as 0.04632:  log likelihood=6.44
## AIC=-6.88   AICc=-6.39   BIC=-0.97

Arima(wmurders, order=c(2,2,0))

## Series: wmurders 
## ARIMA(2,2,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.8289  -0.2246
## s.e.   0.1346   0.1353
## 
## sigma^2 estimated as 0.05292:  log likelihood=3.34
## AIC=-0.68   AICc=-0.19   BIC=5.23

Arima(wmurders, order=c(0,2,2))

## Series: wmurders 
## ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1     ma2
##       -1.0181  0.1470
## s.e.   0.1220  0.1156
## 
## sigma^2 estimated as 0.04702:  log likelihood=6.03
## AIC=-6.06   AICc=-5.57   BIC=-0.15

Arima(wmurders, order=c(2,2,2))

## Series: wmurders 
## ARIMA(2,2,2) 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.3200  0.1947  -0.7027  -0.1776
## s.e.   0.5955  0.2161   0.6030   0.5090
## 
## sigma^2 estimated as 0.04597:  log likelihood=7.63
## AIC=-5.26   AICc=-3.98   BIC=4.59

6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model. The ACF shows that the residuals are within the confidence boundaries as such they’re not correlated with each other. Also the histogram plot shows a normal distribution, and the portmanteau test returns a large p-value which confirms that the residuals are indeed white noise.

checkresiduals(Arima(wmurders, order=c(1,2,1)))

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,1)
## Q* = 12.419, df = 8, p-value = 0.1335
## 
## Model df: 2.   Total lags used: 10

b.

Should you include a constant in the model? Explain. No. In the ARIMA model, d=2. According to the author, Hyndman, “For \(d>1\) no constant is allowed as a quadratic or higher order trend is particularly dangerous when forecasting. The parameter \(\mu\) is called the”drift" in the R output when \(d=1\). Using a constant here would introduce drift into the forecast."

c. 

Write this model in terms of the backshift operator. \((1-\phi_1B)(1-B)^{2}y_t = \epsilon_t(1+\theta_1B)\)

d. 

Fit the model using R and examine the residuals. Is the model satisfactory?

Yes. the model is satisfactory. The residuals exhibit white noise behavior which means that they’re independent of each other.

fit$coef

##        ar1 
## 0.08688777

fit <- Arima(wmurders, order=c(1,2,1))
checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,1)
## Q* = 12.419, df = 8, p-value = 0.1335
## 
## Model df: 2.   Total lags used: 10

e.

Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.

wmurders[55]

## [1] 2.589383

forecast(fit, h=3)

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.470660 2.194836 2.746484 2.048824 2.892496
## 2006       2.363106 1.986351 2.739862 1.786908 2.939304
## 2007       2.252833 1.765391 2.740276 1.507354 2.998313

f. 

Create a plot of the series with forecasts and prediction intervals for the next three periods shown.

autoplot(forecast(fit,h=3))

### g. Does auto.arima() give the same model you have chosen? If not, which model do you think is better?

auto.arima(wmurders)

## Series: wmurders 
## ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.2434  -0.8261
## s.e.   0.1553   0.1143
## 
## sigma^2 estimated as 0.04632:  log likelihood=6.44
## AIC=-6.88   AICc=-6.39   BIC=-0.97