Table of Contents

  1. Introduction
  2. Loading Require Packages and Data File
  3. Descriptive analysis of Time Series Data
  4. Model Selection Among Linear, Quadratic, Seasonal, Cosine, ARMA, ARIMA models 4.1 Stationary Check: ACF/ PACF plots and ADF Test
    4.2 Normality Check: QQ plot and Shapiro Walk Test
  5. Data Preparation for ARIMA model
    5.1 Box Cox Transformation : Stationary Check, Normality Check
    5.2 Log Transformation : Stationary Check, Normality Check
    5.3 Differencing
  6. Finding possible ARIMA model using model specification tools
    6.1 Tool1: ACF-PACF
    6.2 Tool2: EACF
    6.3 Tool3: BIC table
  7. Parameter Estimation
    7.1 Models Testing
    7.2 AIC
    7.3 BIC Scores
    7.4 Goodness of Fit Test
    7.4 Testing for over-fitting
  8. Diagnostic Check for Best Fitting Model
  9. Conclusion
  10. References

1. Introduction

This report examines a dataset that contains yearly minimum extent of Arctic sea ice (measured in million square km) from 1979 to 2022. The dataset was sourced from NASA and each observation was taken by satellite in September of each year, as this is when Arctic sea ice reaches its minimum extent. The aim of the report is to analyze the time series plot and identify any trends, seasonality, changing variance, behaviors, and change points. The stationarity of the time series data will be checked, and if it is found to be non-stationary, appropriate transformations and differencing will be performed to make it stationary. The report will use various model specification tools, such as ACF-PACF, EACF, and BIC table, to find possible ARIMA(p,d,q) models. The best model will be selected based on the model fit, and the parameters will be estimated accordingly.

2. Loading Require Packages and Data File

# Set working directory not to repeat the path to the folder while loading data
setwd("~/Desktop/2nd Year 1st Sem/Time Series/Assignment 2")
artic <- read.csv("assignment2Data2023.csv")
head(artic)
##   Year Arctic.Sea.Ice.Extent..million.square.km.
## 1 1979                                      6.90
## 2 1980                                      7.54
## 3 1981                                      6.90
## 4 1982                                      7.17
## 5 1983                                      7.22
## 6 1984                                      6.43
#Changing column names
colnames(artic) <- c("Day","Value")
head(artic)
##    Day Value
## 1 1979  6.90
## 2 1980  7.54
## 3 1981  6.90
## 4 1982  7.17
## 5 1983  7.22
## 6 1984  6.43

3. Descriptive analysis of Time Series Data

#Convert to time series data
articts <- ts(as.vector(artic$Value), start=1979, end = 2022)
articts
## Time Series:
## Start = 1979 
## End = 2022 
## Frequency = 1 
##  [1] 6.90 7.54 6.90 7.17 7.22 6.43 6.49 7.16 6.96 7.13 6.91 6.04 6.30 7.21 6.18
## [16] 6.96 6.03 7.19 6.63 6.35 5.76 5.98 6.60 5.64 6.01 5.79 5.32 5.77 4.16 5.59
## [31] 5.12 4.62 4.34 3.39 5.05 5.03 4.43 4.17 4.67 4.66 4.19 3.82 4.72 4.67
#Summary Statistics
summary(articts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.390   4.707   5.995   5.800   6.900   7.540
#Time Series Plot
plot(articts,
     type = "o",
     xlab = "Time",
     ylab = "Ice minimum extent (million square km)",
     main = " Time series plot of yearly Artic sea ice minimum extent(1979-2022)",
     cex.main = 0.9)

  1. Trend: The time series displays an overall downward trend with minor changes in local trends.
  2. Seasonality: There appears to be no seasonality in the data, but we will confirm this with ACF/PACF plots.
  3. Changing Variance: There are no significant changes in variance in the time series.
  4. Behaviors (Moving Average/ Autoregressive): Autoregressive behavior is evident as some succeeding time points are correlated, and moving average behavior might be present. We will confirm this with ACF/ PACF plots.
  5. Intervention/ Change Point: A sudden dip is observed around 2007, indicating the presence of an intervention or change point.
#ACF/ PACF plots to check seaonality
acf(articts, main= "ACF plot of ice minimum extent")

pacf(articts, main= "PACF plot of ice minimum extent")

4. Model Selection Among Linear, Quadratic, Seasonal, Cosine, ARMA, ARIMA models

4.1 Stationarity Check: ACF/ PACF plots and ADF test

ACF/ PACF Plots to determine stationarity
  • We will check whether the data is stationary or not by using ACF/ PACF plots.
#ACF/ PACF Plots to determine stationarity 
acf(articts, main= "ACF plot of ice minimum extent",  cex.main = 0.9)

pacf(articts, main= "PACF plot of ice minimum extent",  cex.main = 0.9)

  • The ACF plot’s slow decaying pattern and the PACF plot’s strong correlation in the first lag indicate the existence of an autoregressive (AR) term in the data, which is associated with the trend and suggesting data is non-stationary. To achieve stationarity, differencing can be utilized to eliminate the trend element, and additional changes such as Box-Cox or log transformations might be employed to stabilize the variance if needed.
ADF test to determine stationarity
  • We will perform Augmented Dickey Fuller (ADF) test to check whether the given time series is stationary or not.
    H0 : The series is non-stationary
    HA : The series is stationary 
#ADF test to determine stationarity
adf.test(articts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  articts
## Dickey-Fuller = -2.3296, Lag order = 3, p-value = 0.4433
## alternative hypothesis: stationary
  • Since the p-value = 0.4433 is greater than the significance level of 0.05, we cannot reject the null hypothesis, indicating that the time series data is non-stationary.

4.2 Normality Check : QQ Plot and Shapiro Walk Test

QQ Plot to determine normality
  • The QQ plot will be used to assess whether the data follows a normal distribution. This plot compares the theoretical quantiles of a normal distribution to the actual quantiles of the data. If the points on the QQ plot form a straight line, then the data can be assumed to be normally distributed.
##QQ Plot to determine normality
qqnorm(articts, ylab="Ice minimum extent (million square km)", xlab="Normal Scores")
qqline(articts, col = 2)

  • Based on the QQ plot, it is observed that the data deviates from the normal distribution as most of the points do not fall on the reference line. Therefore, it can be concluded that the data is not normally distributed.
Shapiro Walk Test to determine normality
  • The Shapiro-Wilk test will be conducted to assess the normality of the data. A p-value greater than the significance level of \(\alpha=0.05\) indicates that we fail to reject the null hypothesis, which states that the data is normally distributed. If we do not reject the null hypothesis, it suggests that the residuals are normally distributed.
#Shapiro Walk Test to determine normality
shapiro.test(articts) 
## 
##  Shapiro-Wilk normality test
## 
## data:  articts
## W = 0.94543, p-value = 0.0373

  • Upon examination of the result, it was found that the p-value of 0.0373 is less than the significance level of \(\alpha=0.05\). As a result, the null hypothesis, which is that the data are normally distributed, is rejected, indicating that the data is not normal.

  • From above analysis, ARMA model is also out of consideration as our data is non-stationary. therefore, we will conclude ARIMA model is the most suitable approach as explained above.

5. Data Preparation for ARIMA model

5.1 Box Cox Transformation : Stationary Check, Normality Check

  • After analyzing the data, it was determined that the time series is non-stationary and the data is not normally distributed. In addition, for ARIMA modeling, it is necessary that the variance of the data remains constant over time. Therefore, we will apply both box cox and log transformations to the data to determine which transformation is most suitable for achieving stationarity, and then proceed with differencing the transformed series.
Box Cox Transformation
# Box Cox Transformation
BC <- BoxCox.ar(articts) 

BC$ci
## [1] 0.4 2.0
lambda <- BC$lambda[which(max(BC$loglike) == BC$loglike)]
lambda
## [1] 1.4
BC.artic <- ((articts^lambda)-1)/lambda
  • With a 95% confidence interval of the above plot, the \(\lambda\) value can be found between 0.4 and 2. Hence, \(\lambda\) is approximately 1.4 which is the midpoint of the values with the CI.
#Time Series Plot of Box-Cox Transformed Data
plot(BC.artic, 
     type='o', 
     ylab ="Ice minimum extent (million square km)", 
     main="Time Series Plot of box cox transformed ice minimum extent",
     cex.main = 0.9)

  • The time series plot remains unchanged even after applying the box-cox transformation.
ADF test to determine stationarity of Box Cox Transformed Data
  • We will perform Augmented Dickey Fuller (ADF) test to check whether the box-cox transformed time series is stationary or not.
    H0 : The series is non-stationary
    HA : The series is stationary 
#ADF test of Box-Cox Transformed Data to determine stationarity
adf.test(BC.artic)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  BC.artic
## Dickey-Fuller = -2.369, Lag order = 3, p-value = 0.4277
## alternative hypothesis: stationary
  • The p-value of 0.4277 is greater than the significance level of \(\alpha=0.05\), indicating that we fail to reject the null hypothesis. Therefore, we can conclude that the time series data is non-stationary.
QQ Plot and Shapiro Walk test to determine normality of Box Cox Transformed Data
#QQ Plot of Box-Cox Transformed Data to determine normality
qqnorm(BC.artic, ylab="Ice minimum extent (million square km)", xlab="Normal Scores")
qqline(BC.artic, col = 2)

#Shapiro Walk Test of Box-Cox Transformed Data to determine normality
shapiro.test(BC.artic)
## 
##  Shapiro-Wilk normality test
## 
## data:  BC.artic
## W = 0.94795, p-value = 0.04619

Based on the QQ plot, it can be observed that the tails of the distribution deviate significantly from normality, while the p-value of the Shapiro test is less than the significance level of \(\alpha=0.05\), providing enough evidence to reject the null hypothesis of normality. Therefore, it can be concluded that the Box Cox transformation did not have a significant impact in improving the normality of the observations.

5.2 Log Transformation : Stationary Check, Normality Check

#Log Transformation
log.artic = log(articts)

#Time Series Plot of Log-Transformed Data
plot(log.artic, 
     type='o', 
     ylab ="Ice minimum extent (million square km)", 
     main="Time Series Plot of Log transformed ice minimum extent",
     cex.main = 0.9)

  • The time series plot remains unchanged even after applying the log transformation.
ADF test to determine stationarity of Log Transformed Data
  • We will perform Augmented Dickey Fuller (ADF) test to check whether the log transformed time series is stationary or not.
    H0 : The series is non-stationary
    HA : The series is stationary 
#ADF test Log Transformed Data to determine stationarity
adf.test(log.artic)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  log.artic
## Dickey-Fuller = -2.2867, Lag order = 3, p-value = 0.4603
## alternative hypothesis: stationary
  • The p-value of 0.4603 is greater than the significance level of \(\alpha=0.05\), indicating that we fail to reject the null hypothesis. Therefore, we can conclude that the time series data is non-stationary.
QQ Plot and Shapiro Walk test to determine normality of Log Transformed Data
#QQ Plot of Log Transformed Data to determine normality
qqnorm(log.artic, ylab="Ice minimum extent (million square km)", xlab="Normal Scores")
qqline(log.artic, col = 2)

#Shapiro Walk Test of Log Transformed Data to determine normality
shapiro.test(log.artic)
## 
##  Shapiro-Wilk normality test
## 
## data:  log.artic
## W = 0.93008, p-value = 0.01051

  • Based on the QQ plot, it can be observed that the tails of the distribution deviate significantly from normality, while the p-value of the Shapiro test is less than the significance level of 0.05, providing enough evidence to reject the null hypothesis of normality. Therefore, it can be concluded that the Log transformation did not have a significant impact in improving the normality of the observations.

  • Even though we tried both box-cox and log transformations, neither resulted in significant improvement compared to the original data. However, the box-cox transformed data showed slightly better results. Therefore, we will proceed with the Box-cox transformed data to perform differencing.

5.3 Differencing

  • We will now move on to the process of differencing to obtain a stationary dataset. Sometimes the first difference may not work, then, we can go for the second difference of the series.
First Differencing to Box Cox Transformed Data
diff.BC.artic = diff(BC.artic,differences = 1)
plot(diff.BC.artic, 
     type='o',
     ylab ="Ice minimum extent (million square km)", 
     main="Time Series Plot of first difference of BC ice minimum extent",
     cex.main = 0.9)
abline(h=0, lty=2, col = "red")

  • The differenced series looks much more stationary when compared with the original time series. Time Series Plot of first differencing of box-cox transformed data shows a constant mean and variance over time, then the time series is likely to be stationary.
ACF/ PACF Plots to determine stationarity of First Differencing Box Cox Transformed Data
acf(diff.BC.artic, main= "ACF plot of first differenced BC series",  cex.main = 0.9)

pacf(diff.BC.artic, main= "PACF plot of first differenced BC series",  cex.main = 0.9)

  • The ACF plot is not having decaying pattern and PACF plot doesn’t have very high first lag. So, the series may be stationary. We will proceed with unit roots testing to make sure the data is stationary.
Unit-root tests to determine stationarity of First Differencing Box Cox Transformed Data
# Unit-root tests on first differenced series
adf.test(diff.BC.artic)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff.BC.artic
## Dickey-Fuller = -6.6147, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
pp.test(diff.BC.artic)
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff.BC.artic
## Dickey-Fuller Z(alpha) = -54.357, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary
kpss.test(diff.BC.artic)
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff.BC.artic
## KPSS Level = 0.073296, Truncation lag parameter = 3, p-value = 0.1
  • Looking at p-values of ADF and PP test, the p-value = 0.01 is less than significant level of \(\alpha=0.05\). Thus, we can reject null hypothesis H0: The series is non stationary and conclude the time series data is stationary.
  • However, in KPSS test, not rejecting the null hypothesis means implies stationary. And, as our p value = 0.1 is greater than significant level of \(\alpha=0.05\), we fail to reject null hypothesis which implies the time series data is stationary.
QQ Plot and Shapiro Walk test to determine normality of First Differenced Box Cox Transformed Data
#QQ Plot of 1st Differenced Box-Cox Transformed Data to determine normality
qqnorm(diff.BC.artic, ylab="Ice minimum extent (million square km)", xlab="Normal Scores")
qqline(diff.BC.artic, col = 2)

#Shapiro Walk Test of  1st Differenced Box-Cox Transformed Data to determine normality
shapiro.test(diff.BC.artic)
## 
##  Shapiro-Wilk normality test
## 
## data:  diff.BC.artic
## W = 0.98383, p-value = 0.796
  • Based on the QQ plot, it can be observed that the points are falling on the line which represent normality, as well the p-value of the Shapiro test is greater than the significance level of 0.05, thus, we fail to reject the null hypothesis H0: data is normal. Hence, normality.
  • To conclude, the 1st differenced to box cox transformed data makes the time series stationary and make the data normal.

6. Finding possible ARIMA model using model specification tools

  1. ACF/ PACF Plot
  2. ECAF Plot
  3. BIC Table

6.1 Tool1: ACF-PACF

# Model specification with the first differenced series
acf(diff.BC.artic, main= "ACF plot of first differenced BC series", cex.main = 0.9)

pacf(diff.BC.artic, main= "PACF plot of first differenced BC series", cex.main = 0.9)

  • From the ACF plot, we get q = 1, 2. From the PACF plot, we get p = 2, 3. So, the models suggested by ACF-PACF plots are: ARIMA(2,1,1), ARIMA(2,1,2), ARIMA(3,1,1), ARIMA (3,1,2)

6.2 Tool2: EACF

  • We have to focus on EACF to determine the values of p and q.
# Model specification with the first differenced series
eacf(diff.BC.artic,ar.max = 5, ma.max = 5)
## AR/MA
##   0 1 2 3 4 5
## 0 x o o o o o
## 1 x x o o x o
## 2 x o o o o o
## 3 x o o o o o
## 4 o o o o o o
## 5 o o o o o o
  • The top left most 0 is at location (0,1). The neighbour model is (0,2) and (1,2). Hence, from the EACF we can find that the suggested ARIMA model is: ARIMA(0,1,1), ARIMA(0,1,2), ARIMA(1,1,2)

6.3 Tool3: BIC table

  • Now, we can chart the BIC table to determine a more probable models for the data.
# Model specification with the first differenced series
res = armasubsets(y=diff.BC.artic, nar=5 ,nma=5, y.name='p', ar.method='ols')
plot(res)

  • From the BIC table above, additional probable models that can be taken into consideration is {ARIMA(1,1,2), ARIMA(1,1,4), ARIMA(2,1,2), ARIMA(2,1,4)} as the shaded columns correspond to AR(1), AR(2), MA(2), MA(4). As AR(5) is going to be a big model and not included in ACF/PACF/EACF, we decide to disregard it.

  • Thus, final set of possible models are as follows:

    • {ARIMA(0,1,1)
    • ARIMA(0,1,2)
    • ARIMA (1,1,2)
    • ARIMA (1,1,4)
    • ARIMA(2,1,1)
    • ARIMA(2,1,2)
    • ARIMA (2,1,4)
    • ARIMA (3,1,1)
    • ARIMA (3,1,2)}

7. Parameter Estimation

7.1 Models Testing

ARIMA(0,1,1)
model.011 = arima(BC.artic,order=c(0,1,1), method='ML')
coeftest(model.011)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.624619   0.089687 -6.9644 3.297e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.011CSS = arima(BC.artic,order=c(0,1,1), method='CSS')
coeftest(model.011CSS)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.633040   0.088697 -7.1371 9.529e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • According to coefficient test, MA(1) is significant as the value in less than 0.05 in both test: CSS and ML. Thus, ARIMA(0,1,1) is significant.
ARIMA(0,1,2)
model.012 = arima(BC.artic,order=c(0,1,2), method='ML')
coeftest(model.012)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.79782    0.19634 -4.0635 4.835e-05 ***
## ma2  0.18396    0.17213  1.0688    0.2852    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.012CSS = arima(BC.artic,order=c(0,1,2), method='CSS')
coeftest(model.012CSS)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.80431    0.19690 -4.0848 4.412e-05 ***
## ma2  0.18115    0.17296  1.0474    0.2949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • We are keeping AR as 0 and increased MA to 2, we saw that MA(1) is still significant whereas MA(2) is insignificant as the value is greater than 0.05 in both test: CSS and ML. Thus, ARIMA(0,1,2) is insignificant.
ARIMA(1,1,2)
model.112 = arima(BC.artic,order=c(1,1,2), method='ML')
coeftest(model.112)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.73582    0.11788  6.2419 4.323e-10 ***
## ma1 -1.81253    0.23885 -7.5886 3.233e-14 ***
## ma2  0.99998    0.26092  3.8325 0.0001268 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.112CSS = arima(BC.artic,order=c(1,1,2), method='CSS')
coeftest(model.112CSS)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.49214    0.39095 -1.2589   0.2081
## ma1 -0.20062    0.45422 -0.4417   0.6587
## ma2 -0.22168    0.32054 -0.6916   0.4892
model.112CSSML = arima(BC.artic,order=c(1,1,2),method='CSS-ML') 
coeftest(model.112CSSML) #Prefer ML
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.73582    0.11788  6.2419 4.322e-10 ***
## ma1 -1.81254    0.23886 -7.5883 3.241e-14 ***
## ma2  0.99998    0.26093  3.8324 0.0001269 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Now, we increased AR to 1 and keep MA the same as before. Two coefficient tests are showing two different results. This is because of our box cox transformed data is not normally distributed as shown in Section 4.3. In this situation, we can use CSS-ML method which uses both method in one. CSS-ML shows both AR(1), AR(2) and MA(2) coefficient is significant as well. Thus, ARIMA(1,1,1) is significant.
ARIMA(1,1,4)
model.114 = arima(BC.artic,order=c(1,1,4), method='ML')
coeftest(model.114)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.28937    0.31064  0.9315   0.35158    
## ma1 -1.29269    0.28898 -4.4732 7.705e-06 ***
## ma2  0.52742    0.29141  1.8099   0.07031 .  
## ma3 -0.24846    0.40710 -0.6103   0.54165    
## ma4  0.41647    0.32422  1.2845   0.19896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.114CSS = arima(BC.artic,order=c(1,1,4), method='CSS')
coeftest(model.114CSS)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)    
## ar1 -0.180249   0.260028 -0.6932  0.48819    
## ma1 -0.783891   0.234749 -3.3393  0.00084 ***
## ma2  0.044283   0.340892  0.1299  0.89664    
## ma3 -0.127907   0.274208 -0.4665  0.64089    
## ma4  0.411130   0.164678  2.4966  0.01254 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.114CSSML = arima(BC.artic,order=c(1,1,4),method='CSS-ML') 
coeftest(model.114CSSML) #Prefer ML
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.28937    0.31064  0.9315   0.35158    
## ma1 -1.29269    0.28898 -4.4732 7.705e-06 ***
## ma2  0.52742    0.29141  1.8099   0.07032 .  
## ma3 -0.24845    0.40710 -0.6103   0.54166    
## ma4  0.41647    0.32422  1.2845   0.19896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Now, we are keeping AR the same and increase MA to 4. ML shows MA(1) and MA(2) is significant and CSS shows MA(1) and MA(4) is significant. We tested with CSS-ML and the test favors ML result: MA(1) and MA(2) is significant. Thus, ARIMA(1,1,4) is insignificant.
ARIMA(2,1,1)
model.211 = arima(BC.artic,order=c(2,1,1), method='ML')
coeftest(model.211)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1 -0.26765    0.24156  -1.108  0.26785  
## ar2 -0.11904    0.19482  -0.611  0.54119  
## ma1 -0.47506    0.19761  -2.404  0.01622 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.211CSS = arima(BC.artic,order=c(2,1,1), method='CSS')
coeftest(model.211CSS)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1 -0.25440    0.23763 -1.0705  0.28437  
## ar2 -0.11833    0.19407 -0.6097  0.54203  
## ma1 -0.48960    0.19240 -2.5447  0.01094 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • The results show MA(1) is significant and both AR insignificant. Thus, ARIMA(2,1,1) is insignificant. Thus, ARIMA(2,1,1) is insignificant.
ARIMA(2,1,2)
model.212 = arima(BC.artic,order=c(2,1,2), method='ML')
coeftest(model.212)
## 
## z test of coefficients:
## 
##     Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.80111    0.15582   5.1414 2.727e-07 ***
## ar2 -0.10819    0.16268  -0.6651     0.506    
## ma1 -1.80370    0.14114 -12.7793 < 2.2e-16 ***
## ma2  1.00000    0.15058   6.6412 3.112e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.212CSS = arima(BC.artic,order=c(2,1,2), method='CSS')
coeftest(model.212CSS)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.38248    0.16886  2.2650   0.02351 *  
## ar2 -0.10864    0.22575 -0.4812   0.63035    
## ma1 -1.29920    0.19812 -6.5578 5.462e-11 ***
## ma2  0.60132    0.14913  4.0322 5.526e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • We increased to MA to 2, and now AR(1) and MA(2) are siginificant and as expeced MA(1) is still significant. However, AR(2) is insignificant. Thus, ARIMA(2,1,2) is insignificant.
ARIMA(2,1,4)
model.214 = arima(BC.artic,order=c(2,1,4), method='ML')
coeftest(model.214)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.25619    0.23712  1.0804   0.27996    
## ar2 -0.26754    0.25223 -1.0607   0.28884    
## ma1 -1.28149    0.24837 -5.1596 2.475e-07 ***
## ma2  0.85397    0.34559  2.4711   0.01347 *  
## ma3 -0.78053    0.45679 -1.7087   0.08750 .  
## ma4  0.71684    0.32072  2.2351   0.02541 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.214CSS = arima(BC.artic,order=c(2,1,4), method='CSS')
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
coeftest(model.214CSS)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.0823315  0.0014987  54.9346 < 2.2e-16 ***
## ar2  0.4676898  0.0089673  52.1552 < 2.2e-16 ***
## ma1 -1.5147617  0.0362826 -41.7490 < 2.2e-16 ***
## ma2  0.0584274  0.1347703   0.4335    0.6646    
## ma3  0.9061293  0.2114652   4.2850 1.827e-05 ***
## ma4 -0.0151237  0.1207263  -0.1253    0.9003    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.214CSSML = arima(BC.artic,order=c(2,1,4),method='CSS-ML') 
coeftest(model.214CSSML) #Prefer ML
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.25619    0.23712  1.0804   0.27996    
## ar2 -0.26754    0.25223 -1.0607   0.28884    
## ma1 -1.28149    0.24837 -5.1596 2.475e-07 ***
## ma2  0.85397    0.34559  2.4711   0.01347 *  
## ma3 -0.78053    0.45679 -1.7087   0.08750 .  
## ma4  0.71684    0.32072  2.2351   0.02541 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • CSS and ML are showing two different result. When checked with CSS-ML, the results shows all coefficients of MA being significant and AR being insignificant. Thus, RIMA(2,1,4) is insignificant.
ARIMA(3,1,1)
model.311 = arima(BC.artic,order=c(3,1,1), method='ML')
coeftest(model.311)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.45129    0.29735 -1.5177   0.1291
## ar2 -0.27370    0.24055 -1.1378   0.2552
## ar3 -0.17910    0.17667 -1.0138   0.3107
## ma1 -0.30602    0.27014 -1.1328   0.2573
model.311CSS = arima(BC.artic,order=c(3,1,1), method='CSS')
coeftest(model.311CSS)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.44972    0.28509 -1.5775   0.1147
## ar2 -0.29013    0.23739 -1.2222   0.2216
## ar3 -0.19020    0.17446 -1.0903   0.2756
## ma1 -0.30295    0.25925 -1.1686   0.2426
  • We dialed MA back to 1 and increased AR to 3 and now all the coefficients are insignificant. Thus, ARIMA(3,1,1) is insignificant.
ARIMA(3,1,2)
model.312 = arima(BC.artic,order=c(3,1,2), method='ML')
coeftest(model.312)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.796484   0.157318   5.0629 4.130e-07 ***
## ar2 -0.074527   0.196206  -0.3798    0.7041    
## ar3 -0.049002   0.156986  -0.3121    0.7549    
## ma1 -1.801055   0.126791 -14.2049 < 2.2e-16 ***
## ma2  0.999999   0.132851   7.5272 5.182e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.312CSS = arima(BC.artic,order=c(3,1,2), method='CSS')
coeftest(model.312CSS)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.409206   0.099848  -4.0983 4.162e-05 ***
## ar2 -0.920922   0.061145 -15.0613 < 2.2e-16 ***
## ar3 -0.505729   0.109665  -4.6116 3.996e-06 ***
## ma1 -0.192102   0.058683  -3.2736  0.001062 ** 
## ma2  1.218802   0.058466  20.8464 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.312CSSML = arima(BC.artic,order=c(3,1,2),method='CSS-ML') 
coeftest(model.312CSSML) #Prefer CSS
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.368698   0.138525  -2.6616  0.007777 ** 
## ar2 -0.848873   0.043352 -19.5807 < 2.2e-16 ***
## ar3 -0.607568   0.129032  -4.7087 2.494e-06 ***
## ma1 -0.299872   0.103530  -2.8965  0.003774 ** 
## ma2  0.999932   0.155715   6.4216 1.349e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  • As CSS and ML tests are showing different result, we performed another test CSS-ML and the results shows all the coefficients are significant. Thus, ARIMA(3,1,2) is insignificant.

  • For now, we will consider the AIC and BIC values of the models to decide the best one within the subset of possible models.

7.2 AIC Score

sort.score(AIC(model.011,model.012,model.112,model.114,model.211,model.212,model.214,model.311,model.312), score = "aic")
##           df      AIC
## model.112  4 131.4230
## model.212  5 132.9743
## model.114  6 134.1905
## model.312  6 134.8760
## model.214  7 135.7616
## model.011  2 136.0854
## model.012  3 136.8267
## model.211  4 138.9046
## model.311  5 140.0541

7.3 BIC Scores

sort.score(BIC(model.011,model.012,model.112,model.114,model.211,model.212,model.214,model.311,model.312), score = "bic" ) 
##           df      BIC
## model.112  4 138.4678
## model.011  2 139.6078
## model.212  5 141.7803
## model.012  3 142.1103
## model.114  6 144.7577
## model.312  6 145.4432
## model.211  4 145.9494
## model.214  7 148.0900
## model.311  5 148.8601
  • Based on AIC and BIC Scores, the minimum is ARIMA (1,1,2) which is 131.4230 and 138.4678 respectively. Hence, ARIMA(1,1,2) is the best model. Sometimes, AIC and BIC might not be consistent. Thus, we will go for other goodness of fit measures such as mean squared error and mean absolute error, etc.

7.4 Goodness of Fit Test

#Goodness of Fit Test
model.011A = Arima(BC.artic,order=c(0,1,1),method='ML') 
model.012A = Arima(BC.artic,order=c(0,1,2),method='ML') 
model.112A = Arima(BC.artic,order=c(1,1,2),method='ML') 
model.114A = Arima(BC.artic,order=c(1,1,4),method='ML') 
model.211A = Arima(BC.artic,order=c(2,1,1),method='ML') 
model.212A = Arima(BC.artic,order=c(2,1,2),method='ML') 
model.214A = Arima(BC.artic,order=c(2,1,4),method='ML') 
model.311A = Arima(BC.artic,order=c(3,1,1),method='ML') 
model.312A = Arima(BC.artic,order=c(3,1,2),method='ML')

Smodel.011A <- accuracy(model.011A)[1:7] 
Smodel.012A <- accuracy(model.012A)[1:7] 
Smodel.112A <- accuracy(model.112A)[1:7] 
Smodel.114A <- accuracy(model.114A)[1:7] 
Smodel.211A <- accuracy(model.211A)[1:7] 
Smodel.212A <- accuracy(model.212A)[1:7] 
Smodel.214A <- accuracy(model.214A)[1:7] 
Smodel.311A <- accuracy(model.311A)[1:7] 
Smodel.312A <- accuracy(model.312A)[1:7] 
df.Smodels <- data.frame(rbind(Smodel.011A,
                               Smodel.012A,
                               Smodel.112A, 
                               Smodel.114A,
                               Smodel.211A,
                               Smodel.212A, 
                               Smodel.214A,
                               Smodel.311A,
                               Smodel.312A))
colnames(df.Smodels) <- c("ME", "RMSE", "MAE", "MPE", "MAPE","MASE", "ACF1")
rownames(df.Smodels) <- c("ARIMA(0,1,1)", "ARIMA(0,1,2)", "ARIMA(1,1,2)", "ARIMA(1,1,4)", "ARIMA(2,1,1)", "ARIMA(2,1,2)", "ARIMA(2,1,4)", "ARIMA(3,1,1)", "ARIMA(3,1,2)")
df_sorted <- df.Smodels[order(df.Smodels$RMSE), ]
df_sorted
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## ARIMA(2,1,4) -0.2305750 0.9162520 0.7408581 -4.615092 11.37548 0.6514338
## ARIMA(1,1,4) -0.2007748 0.9301629 0.7640564 -4.394948 11.62921 0.6718320
## ARIMA(2,1,2) -0.1782252 0.9387239 0.7636391 -4.150710 11.68170 0.6714651
## ARIMA(3,1,2) -0.1875978 0.9393466 0.7632726 -4.307207 11.67088 0.6711428
## ARIMA(1,1,2) -0.1587656 0.9409922 0.7742585 -3.832943 11.85723 0.6808026
## ARIMA(3,1,1) -0.3215838 1.0762527 0.8557404 -6.790272 13.30575 0.7524494
## ARIMA(0,1,2) -0.3008403 1.0865081 0.8789888 -6.493769 13.70052 0.7728916
## ARIMA(2,1,1) -0.3068757 1.0880547 0.8814228 -6.608009 13.70375 0.7750319
## ARIMA(0,1,1) -0.3062147 1.1048327 0.9068893 -6.637663 14.04689 0.7974244
##                     ACF1
## ARIMA(2,1,4) -0.02180762
## ARIMA(1,1,4) -0.05384995
## ARIMA(2,1,2) -0.05667410
## ARIMA(3,1,2) -0.05180574
## ARIMA(1,1,2)  0.02964036
## ARIMA(3,1,1) -0.13557784
## ARIMA(0,1,2) -0.05823374
## ARIMA(2,1,1) -0.11553479
## ARIMA(0,1,1) -0.20874144
  • Based on the goodness of fit test, ARIMA(2,1,4) is the best model as it’s Root Mean Squared Error (RMSE), Mean Absolute Error (MAE) and Mean Absolute Scaled Error (MASE) value are the lowest. However, when we looked back the coefficient test of ARIMA(2,1,4), it’s not significant. The same applied to ARIMA(1,1,4), ARIMA(2,1,2). Those models are not significant. For ARIMA(3,1,2), we can see from coefficient tests that, with ML it is not significant, however, with CSS and CSS-ML, the model is significant. Here, we choose to assume this model as insignificant as we are focusing on ML coefficient test rather than CSS.
  • Thus, we conclude the model with fifth minimum RMSE, MAE, MASE value as the best fitting model which is ARIMA(1,1,2). We can support this decision by looking at AIC and BIC tables. We observe the smallest AIC and BIC with model is ARIMA(1,1,2).

7.5 Testing for over-fitting

  • Now we will test the overfitting of our best model which is ARIMA(1,1,2). To do so we will check whether ARIMA(1,1,1), ARIMA(1,1,3), ARIMA(0,1,2), ARIMA(0,1,3), ARIMA(2,1,2) is significant or not.
  • We have already checked ARIMA(0,1,2) and ARIMA(2,1,2) is significant or not and we found both of them is insignificant.
  • So, now we will only test whether ARIMA(1,1,1), ARIMA(1,1,3) and ARIMA(0,1,3) is significant or not. 
# Testing for over-fitting 
model.111 = arima(BC.artic,order=c(1,1,1),method='ML')
coeftest(model.111) 
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)    
## ar1 -0.17921    0.19099 -0.9383   0.3481    
## ma1 -0.54925    0.13656 -4.0220 5.77e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.113 = arima(BC.artic,order=c(1,1,3),method='ML')
coeftest(model.113)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.66920    0.17866  3.7457 0.0001799 ***
## ma1 -1.68351    0.25955 -6.4864 8.792e-11 ***
## ma2  0.77883    0.38440  2.0261 0.0427535 *  
## ma3  0.12246    0.19853  0.6168 0.5373467    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.013 = arima(BC.artic,order=c(0,1,3),method='ML')
coeftest(model.013)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.906733   0.167560 -5.4114 6.253e-08 ***
## ma2  0.011805   0.192911  0.0612   0.95121    
## ma3  0.300257   0.181968  1.6501   0.09893 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • We can see that all the models are insignificant. Thus, we can assume ARIMA(1,1,2) is the best model and continue with residual analysis and forecasting.

8. Diagnostic Check for Best Fitting Model

#Residual Analysis
x11()
res.112CSSML = residuals(model.112CSSML)
par(mfrow=c(3,2))
plot(y = res.112CSSML, 
     type = 'o', 
     x = as.vector(time(BC.artic)),
     xlab = 'Time', 
     ylab='Standardized Residuals',
     main = "Standardised residuals.")
abline(h=0)
hist(res.112CSSML,xlab='Standardized Residuals', main = "Histogram of standardised residuals.")
qqnorm(y=res.112CSSML, main = "QQ plot of standardised residuals.")
qqline(y=res.112CSSML, col = 2, lwd = 1, lty = 2)
shapiro.test(res.112CSSML)
## 
##  Shapiro-Wilk normality test
## 
## data:  res.112CSSML
## W = 0.97602, p-value = 0.4828
acf(res.112CSSML, main = "ACF of standardized residuals.")
pacf(res.112CSSML, main = "PACF of standardized residuals.")

9. Conclusion

10. References