Time Series Final Assignment

Setup

Setting up the working directory

setwd("C:/Users/Saikat/Desktop/Time Series/Final Assignment")

Packages

Loading the necessary packages

library(TSA)

## Warning: package 'TSA' was built under R version 3.5.3

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(fUnitRoots)

## Warning: package 'fUnitRoots' was built under R version 3.5.3

## Loading required package: timeDate

## 
## Attaching package: 'timeDate'

## The following objects are masked from 'package:TSA':
## 
##     kurtosis, skewness

## Loading required package: timeSeries

## Warning: package 'timeSeries' was built under R version 3.5.3

## Loading required package: fBasics

## Warning: package 'fBasics' was built under R version 3.5.3

library(forecast)
library(CombMSC)

## Warning: package 'CombMSC' was built under R version 3.5.3

## 
## Attaching package: 'CombMSC'

## The following object is masked from 'package:stats':
## 
##     BIC

library(lmtest)

## Warning: package 'lmtest' was built under R version 3.5.3

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following object is masked from 'package:timeSeries':
## 
##     time<-

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(fGarch)

## Warning: package 'fGarch' was built under R version 3.5.3

library(rugarch)

## Warning: package 'rugarch' was built under R version 3.5.3

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:stats':
## 
##     sigma

library(dLagM)

## Warning: package 'dLagM' was built under R version 3.5.3

## Loading required package: nardl

## Warning: package 'nardl' was built under R version 3.5.3

## Loading required package: dynlm

## Warning: package 'dynlm' was built under R version 3.5.3

## 
## Attaching package: 'dLagM'

## The following object is masked from 'package:forecast':
## 
##     forecast

library(FitAR)

## Warning: package 'FitAR' was built under R version 3.5.3

## Loading required package: lattice

## Loading required package: leaps

## Warning: package 'leaps' was built under R version 3.5.3

## Loading required package: ltsa

## Warning: package 'ltsa' was built under R version 3.5.2

## Loading required package: bestglm

## Warning: package 'bestglm' was built under R version 3.5.3

## 
## Attaching package: 'FitAR'

## The following object is masked from 'package:forecast':
## 
##     BoxCox

Reading the data

Here the read.csv() function has been used to read the data of the number of New Zealand Air Passengers departed the country from January,2000 to October,2012.We checked the column names of the dataset, then considered the only column on which we have worked.The datatype of the dataset is currently in ‘integer’ format.

Departure<-read.csv("C:/Users/Saikat/Desktop/Time Series/Final Assignment/NZAirPassenger.csv")

View(Departure)

Dept <- Departure[,3]
View(Dept)

class(Dept)

## [1] "integer"

Objective

The objective of this analysis is to analyze the number of New Zealand air passengers leaving their country during Jnauary,2000 to October, 2012 by using various time series analysis methods and finally to choose the best model among a set of possible models for this dataset and give forecasts of departure of New Zealand air passengers for the next 10 time units(in this case, months).

Converting the dataset into Time Series and plotting the graph

Here the dataset has been converted into a time series, to be ensured the class() function has been used, And finally the time series values has been put into a graph.

Dept <- ts(as.vector(Dept), start=c(2000,1), end=c(2012,10),frequency = 12) # Convert to the TS object!
plot(Dept,type='o',ylab='Departures',main = "Fig 1:NZ Departures from 2000 to 2012")

class(Dept)

## [1] "ts"

From the plot we are trying to find out the 5 main characteristics of the data:

Trends: The overall series seems to have an upward trend.

Seasonality: Seasonality is very much obvious.

Intervention/Changing point:The data has upward and downward movements throughout the series, but that is solely because of the presence of the seasonality.The data does not look to have any changing point as there is no drastic shift in the flow which is followed by constant data.

Changing Variance: Changing variance is not clear. We have to check for ARCH and GARCH.

Behaviour: The series looks like non-stationary and has a non constant mean which changes with respect to time. The series seems to have Auto Regressive as well as Moving Average characteristics, but due to the presence of seasonality, we cannot fully determine the auto-correlation structure.

Confirmation of the presence of Seasonality and Trend

# Plotting ACF and PACF of base plots
par(mfrow=c(1,2))
acf(Dept, lag.max = 36)
pacf(Dept, lag.max = 36)

par(mfrow=c(1,1))

From the ACF, it is clearly visible that the seasonal lags at 1S, 2S, and 3S(where S=12) have a slowly decaying pattern, which confirms the presence of seasonality in the timeseries.

On the other hand, the ordinary bits also have a slowly decaying pattern which confirms the existence of the ordinary trend in the dataset.

Checking the presence of Changing Variance

To be sure of the presence of the of changing variance first we need to observe the volatility clustering of the returned series. Also we need to plot the ACF and PACF of the returned series as well as need to do the McLeod-Li test and the QQ test to check whether ARCH-GARCH effect present in the dataset or not.

#Plotting the returned series

par(mfrow=c(1,1))
r.dept=diff(log(Dept))*100
plot(r.dept,main="Fig 2:Monthly Departures: January 2000 to October 2012")
abline(h=0)

#Plotting the ACF,PACF and EACF of the returned series

par(mfrow=c(1,2))
acf(r.dept, lag.max = 36)
pacf(r.dept, lag.max = 36)

par(mfrow=c(1,1))

eacf(r.dept)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o x o x o x x x  x  x  x 
## 1 x x x o x o x o o x o  x  x  x 
## 2 o o x o o o x o o o o  x  o  o 
## 3 x o x o o o x x o o o  x  x  o 
## 4 x x o o x o x o o o o  x  o  o 
## 5 x x o o x o x o o o o  x  x  x 
## 6 o x o o x x x o o o o  x  o  x 
## 7 x o x x x o x o o o o  x  x  o

#Doing the McLeod-Li test and QQ test on the returned series

McLeod.Li.test(y=r.dept,main="Fig 3:McLeod-Li Test Statistics for Departure series")

qqnorm(r.dept,main="Fig 4:Q-Q Normal Plot of Departure series")
qqline(r.dept)

From the returned series plot we do not observe any volatility clustering. The variation present in the plot is only due to the effect of seasonality. Furthermore, we have also plotted the ACF and PACF of the returned series so that we can tally that with the ACF and PACF obtained from the non-linear transformations(in this case absolute value and square root).If the number of significant lags obtained from ACF and PACF of the returned series increase in the ACF and PACF obtained from the transformed series, that means there is no ARCG-GARCH effect present in the series. Also we have done the McLeod-Li test and the QQ test to be sure of the presence of changing variance in the series. From the McLeod-Li test we can observe that not all the p-values are statistically significant,some are above the significant level, which suggests that there is no strong ARCH-GARCH effect present in the series. From the QQ plot we can observe that it does not have fat tails which also denotes the absence of the ARCH-GARCH effect in the series.

Doing the non-linear transformations

#Absolute value transformation

abs.dept = abs(r.dept)


#Plotting the ACF,PACF and EACF of the absolute value transformed series

par(mfrow=c(1,2))
acf(abs.dept,lag.max = 36)
pacf(abs.dept, lag.max = 36)

par(mfrow=c(1,1))
eacf(abs.dept)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o x o x o x o x o x o  x  o  o 
## 1 x x o x o x o x o o o  x  x  x 
## 2 x x o x o o o x o o o  x  x  x 
## 3 x x o x o o o x o o o  x  o  x 
## 4 o x x o o o o o o o o  x  x  o 
## 5 o x o o o o o o o o o  x  x  x 
## 6 o x o o o x o o o o o  x  o  x 
## 7 o x x o o x o o o o o  x  o  x

#Doing the McLeod-Li test and QQ test on the absolute value transformed series

McLeod.Li.test(y=abs.dept,main="Fig 5:McLeod-Li Test Statistics for Departure series")

qqnorm(abs.dept,main="Fig 6:Q-Q Normal Plot of Departure series")
qqline(abs.dept)

The ACF and PACF obtained from the absolute value transformed series shows that number of significant lags have been increased from the number of significant lags present in the ACF and PACF of the returned series. This suggests that there is no changing variance present in the series. Moreover, the McLeod-Li and the QQ test of the absolute value transformed series also confirms the same.

#Square root transformation

sq.dept = r.dept^2


#Plotting the ACF,PACF and EACF of the square root transformed series

par(mfrow=c(1,2))
acf(sq.dept,lag.max = 36)
pacf(sq.dept, lag.max = 36)

par(mfrow=c(1,1))
eacf(sq.dept)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o x o x o x o x o  x  o  o 
## 1 x o o x o x o x o o o  x  x  o 
## 2 x x o o o o o x o o o  x  x  o 
## 3 x x o x o o o x o o o  x  o  x 
## 4 o x x o o o o o o o o  x  x  x 
## 5 o x x o o o o o o o o  x  x  x 
## 6 x x x o o x o o o o o  x  o  x 
## 7 x x x o o x o o o o o  x  o  x

#Doing the McLeod-Li test and QQ test on the  series

McLeod.Li.test(y=sq.dept,main="Fig 7:McLeod-Li Test Statistics for Departure series")

par(mfrow=c(1,1))

qqnorm(sq.dept,main="Fig 8:Q-Q Normal Plot of Departure series")
qqline(sq.dept)

The ACF and PACF obtained from the square root transformed series shows that number of significant lags have been increased from the number of significant lags present in the ACF and PACF of the returned series. This suggests that there is no changing variance present in the series. Moreover, the McLeod-Li and the QQ test of the sqaure root transformed series also confirms the same.

Removing Seasonal Trend

As we have become sure that there is no ARCH-GARCH effect present in the series, now we should remove the seasonality trend present in the series. To do so, we have applied the residual approach of the SARIMA modelling. At first, we have opted for a plain model with only the first seasonal difference with order D = 1 and checked whether we have successfully got rid of the seasonal trend effect by inspecting the autocorrelation structure of the residuals.

#D=1

m1.dept = arima(Dept,order=c(0,0,0),seasonal=list(order=c(0,1,0), period=12))
res.m1 = residuals(m1.dept);  
par(mfrow=c(1,1))
plot(res.m1,xlab='Time',ylab='Residuals',main="Fig 9:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m1, lag.max = 36, main = "Fig 10:The sample ACF of the residuals")
pacf(res.m1, lag.max = 36, main = "Fig 11:The sample PACF of the residuals")

From the ACF we can observe there is one significant lag at the second seasonal lag(2S) which means Q=1.

From the PACF, we can observe a slowly decaying pattern at the seasonal lags, which implies P=0.

Both ACF and PACF suggests that there are still seasonal components present in the series. So, we have added Q=1 and checked whether we have successfully got rid of the seasonal trend by inspecting the autocorrelation structure of the residuals.

#D=1,Q=1

m2.dept = arima(Dept,order=c(0,0,0),seasonal=list(order=c(0,1,1), period=12))
res.m2 = residuals(m2.dept);  
par(mfrow=c(1,1))
plot(res.m2,xlab='Time',ylab='Residuals',main="Fig 12:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m2, lag.max = 50, main = "Fig 13:The sample ACF of the residuals")
pacf(res.m2, lag.max = 50, main = "Fig 14:The sample PACF of the residuals")

Although we have got one significant lag at the first seasonal lag of the ACF, we have already captured that one before. But the trend in the ordinary bit is still present which implies that the local trend in the series is stil present.

The PACF again comes up with a slowly decaying pattern at the seasonal lags, which implies P=0

So we can say that we have successfully detrended the seasonal trend but the local trend is still present.

Removing LOcal Trend

To remove the local trend, we have put d=1 (first differencing) and checked whether we have successfully got rid of the local trend by inspecting the autocorrelation structure of the residuals.

m3.dept = arima(Dept,order=c(0,1,0),seasonal=list(order=c(0,1,1), period=12))
res.m3 = residuals(m3.dept);  
par(mfrow=c(1,1))
plot(res.m3,xlab='Time',ylab='Residuals',main="Fig 15:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m3, lag.max = 36, main = "Fig 16:The sample ACF of the residuals")
pacf(res.m3, lag.max = 36, main = "Fig 17:The sample PACF of the residuals")

From the ACF and PACF we can clearly see the oridnary trend is gone, that means we have successfully detrended the series.

Moreover,a slowly decaying pattern of the significant lags is clearly visible in the ACF which denotes q=0

From the PACF, we can observe two significant lags which means p=1,2.

But this ACF and PACF is not conclusive enough to make the final model. So we will plot an EACF of the residuals for better understanding.

eacf(res.m3)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x x o x x o x x o  o  o  o 
## 1 x o x x o o x o x o o  o  o  o 
## 2 x o x o o o o o o o o  o  o  o 
## 3 x x x x o o o o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x x x o o o o o o o  o  o  o 
## 7 x x x x o o o o o o o  o  o  o

From the EACF, we have got the value of q= 3,4.

So our probable models for prediction are:-

SARIMA{(1,1,0)X(0,1,1)_12} SARIMA{(2,1,0)X(0,1,1)_12} SARIMA{(2,1,3)X(0,1,1)_12} SARIMA{(2,1,4)X(0,1,1)_12}

Model Selection

To check which model is best, we have done the coefficient test, plotted the residuals and also plotted the ACF and PACF of the residuals of every probable models.

SARIMA{(1,1,0)X(0,1,1)_12}

model2_110_011.dept = arima(Dept,order=c(1,1,0),seasonal=list(order=c(0,1,1), period=12))
coeftest(model2_110_011.dept)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.501737   0.073360 -6.8394 7.953e-12 ***
## sma1 -0.561061   0.085163 -6.5881 4.455e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.m2 = residuals(model2_110_011.dept);  
par(mfrow=c(1,1))
plot(res.m2,xlab='Time',ylab='Residuals',main="Fig 18:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m2, lag.max = 36, main = "Fig 19:The sample ACF of the residuals")
pacf(res.m2, lag.max = 36, main = "Fig 20:The sample PACF of the residuals")

Although the coefficient test is suggesting that the SARIMA{(1,1,0)(0,1,1)_12} is significant, but we have got significant lags in both the ACF and PACF of this model. Thus, we are not considering this model as our best model.

SARIMA{(2,1,0)X(0,1,1)_12}

model3_210_011.dept = arima(Dept,order=c(2,1,0),seasonal=list(order=c(0,1,1), period=12))
coeftest(model3_210_011.dept)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.705727   0.077870 -9.0629 < 2.2e-16 ***
## ar2  -0.404252   0.079215 -5.1032 3.339e-07 ***
## sma1 -0.650975   0.083086 -7.8349 4.691e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.m3 = residuals(model3_210_011.dept);  
par(mfrow=c(1,1))
plot(res.m3,xlab='Time',ylab='Residuals',main="Fig 21:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m3, lag.max = 36, main = "Fig 22:The sample ACF of the residuals")
pacf(res.m3, lag.max = 36, main = "Fig 23:The sample PACF of the residuals")

Although the coefficient test is suggesting that the SARIMA{(2,1,0)(0,1,1)_12} is significant, but we have got significant lags in both the ACF and PACF of this model. Thus, we are not considering this model as our best model.

SARIMA{(2,1,3)X(0,1,1)_12}

model4_213_011.dept = arima(Dept,order=c(2,1,3),seasonal=list(order=c(0,1,1), period=12))
coeftest(model4_213_011.dept)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.955053   0.057797 -16.5244 < 2.2e-16 ***
## ar2  -0.929906   0.094399  -9.8508 < 2.2e-16 ***
## ma1   0.274306   0.102012   2.6890  0.007167 ** 
## ma2   0.455904   0.115656   3.9419 8.084e-05 ***
## ma3  -0.431960   0.155682  -2.7746  0.005526 ** 
## sma1 -0.722717   0.108330  -6.6714 2.533e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.m4 = residuals(model4_213_011.dept);  
par(mfrow=c(1,1))
plot(res.m4,xlab='Time',ylab='Residuals',main="Fig 24:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m4, lag.max = 36, main = "Fig 25:The sample ACF of the residuals")
pacf(res.m4, lag.max = 36, main = "Fig 26:The sample PACF of the residuals")

The coefficient test of the SARIMA{(2,1,3)(0,1,1)_12} model is suggesting that this model is significant. Moreover from the ACF and PACF we are getting white noise (no significant lags). So it is one of our probable best models for prediction.

SARIMA{(2,1,4)X(0,1,1)_12}

model5_214_011.dept = arima(Dept,order=c(2,1,4),seasonal=list(order=c(0,1,1), period=12))
coeftest(model5_214_011.dept)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.92321    0.12850 -7.1844 6.749e-13 ***
## ar2  -0.83078    0.12876 -6.4520 1.104e-10 ***
## ma1   0.29064    0.14381  2.0210  0.043275 *  
## ma2   0.35287    0.11467  3.0771  0.002090 ** 
## ma3  -0.35272    0.11091 -3.1803  0.001471 ** 
## ma4  -0.14497    0.13805 -1.0502  0.293647    
## sma1 -0.69364    0.08948 -7.7519 9.055e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.m5 = residuals(model5_214_011.dept);  
par(mfrow=c(1,1))
plot(res.m5,xlab='Time',ylab='Residuals',main="Fig 27:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m5, lag.max = 36, main = "Fig 28:The sample ACF of the residuals")
pacf(res.m5, lag.max = 36, main = "Fig 29:The sample PACF of the residuals")

The coefficient test of the SARIMA{(2,1,4)(0,1,1)_12} model is suggesting that this model is significant. Moreover from the ACF and PACF we are getting white noise (no significant lags). So it is one of our probable best models for prediction.

Sorting AIC and BIC values

Now we will sort the AIC and BIC value of all the probable models. The model which will have the least value among all the models will be our best fit.

source("C:\\Users\\Saikat\\Desktop\\Time Series\\Final Assignment\\sort.score.R")

sort.score(AIC(model4_213_011.dept, model5_214_011.dept), score = "aic")

##                     df      AIC
## model4_213_011.dept  7 3099.567
## model5_214_011.dept  8 3100.317

#sort.score(BIC(model4_213_011.dept, model5_214_011.dept), score = "bic")

According to the AIC and BIC table we can conclude that SARIMA{(2,1,3)(0,1,1)} is the best fit to predict on this particular data.

Model Overfitting

Now we will test the overfitting of our best model SARIMA{(2,1,3)(0,1,1)_12}. To do so we have to check whether SARIMA{(2,1,4)(0,1,1)_12} and SARIMA{(3,1,3){0,1,1}_12} is significant or not.Among them, we have already checked SARIMA{(2,1,4)(0,1,1)_12} and choose SARIMA{(2,1,3)(0,1,1)_12} over it for final model.So now we will only test whether SARIMA{(3,1,3){0,1,1}_12} is significant or not. If SARIMA{(3,1,3){0,1,1}_12} comes insignificant, we can proceed with the SARIMA{(2,1,3)(0,1,1)_12} model to do the forecast.

model6_313_011.dept = arima(Dept,order=c(3,1,3),seasonal=list(order=c(0,1,1), period=12))
coeftest(model6_313_011.dept)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.760190   0.256941 -2.9586 0.0030902 ** 
## ar2  -0.694947   0.288517 -2.4087 0.0160101 *  
## ar3   0.173006   0.207506  0.8337 0.4044274    
## ma1   0.106191   0.240668  0.4412 0.6590420    
## ma2   0.332735   0.172496  1.9289 0.0537375 .  
## ma3  -0.466900   0.132509 -3.5235 0.0004258 ***
## sma1 -0.699087   0.096862 -7.2173 5.302e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.m6 = residuals(model6_313_011.dept);  
par(mfrow=c(1,1))
plot(res.m6,xlab='Time',ylab='Residuals',main="Fig 30:Time series plot of the residuals")

par(mfrow=c(1,2))
acf(res.m6, lag.max = 36, main = "Fig 31:The sample ACF of the residuals")
pacf(res.m6, lag.max = 36, main = "Fig 32:The sample PACF of the residuals")

From the coefficient test it is visible that SARIMA{(3,1,3)(0,1,1)_12} is insignificant, so we can proceed with the SARIMA{(2,1,3)(0,1,1)_12} modle for our final prediction.

Residual Analysis

Here we will do the residual analysis of our predictive model which is SARIMA{(2,1,3)(0,1,1)_12}. This residual analysis includes time series plot of the standardised residuals,histogram of the standardised residuals, QQ Plot of the standardised residuals,ACF of the standardised residuals,Shapiro-Wilk test of the standardised residuals and the Ljung-Box test.

residual.analysis <- function(model, std = TRUE){
  library(TSA)
  library(FitAR)
  if (std == TRUE){
    res.model = rstandard(model)
  }else{
    res.model = residuals(model)
  }
  par(mfrow=c(3,2))
  plot(res.model,type='o',ylab='Standardised residuals', main="Time series plot of standardised residuals")
  abline(h=0)
  hist(res.model,main="Histogram of standardised residuals")
  qqnorm(res.model,main="QQ plot of standardised residuals")
  qqline(res.model, col = 2)
  acf(res.model,main="ACF of standardised residuals")
  print(shapiro.test(res.model))
  k=0
  LBQPlot(res.model, lag.max = length(model$residuals)-1 , StartLag = k + 1, k = 0, SquaredQ = FALSE)
}

residual.analysis(model =model4_213_011.dept )

## 
##  Shapiro-Wilk normality test
## 
## data:  res.model
## W = 0.96778, p-value = 0.001141

## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple
## of shorter object length

The residual analysis indicates

-standardised normals are randomly distributed over time. -histogram looks normally distributed and form a bell shaped curve. -QQ plot looks almost alligned. -ACF plot does not have any significant lags -Although Shapiro-Wilk normality test indicates that residuals are not normally distributed with the p-value of 0.001 -All p-values in Ljung-Box test plot are above the significant level.

Forecasting

We will now do the final forecastong based on the SARIMA{(2,1,3)(0,1,1)_12} model for the upcoming 10 months.

library(forecast)

m4_213.dept = Arima(Dept,order=c(2,1,3),seasonal=list(order=c(0,1,1), period=12))
#future = forecast(m4_213.dept, h = 10)
#plot(future, main ="Fig 33:Forecast of NZ Departures for Next 10 Months[SARIMA{(2,1,3)(0,1,1)_12}]", xlab="Years", ylab="Departures")

Conclusion

The values for the next 10 time points have been predicted using the SARIMA{(2,1,3)(0,1,1)_12} model.The model has predicted with an accuracy of 80% that it will lie in the dark blue region and with an accuracy of 95% that it will lie in the dark grey area.