Introduction

We have been provided by a data that represents the yearly changes in the thickness of Ozone layer for the years (1927-2016) in Dobson units. The main aim of this experiment is to:

Task 1 - Implementing model building strategies to find the best fitting model among linear, quadratic, cosine, cyclical or seasonal trend models. Using the best model, provide the prediction of yearly changes for the next 5 years.

Task 2 - Using all suitable model specification tools such as ACF-PACF, EACF, BIC table, propose a set of possible ARIMA(p, d, q) models.

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# loading dataset and adding a header.
dataset <- read.csv("data1.csv", header = FALSE)
names(dataset)[1] <- "Ozone_Level"
# creating Time Series of dataset. 
ozoneLevelTS <- ts(dataset$Ozone_Level, start =1927, end =2016, frequency =1)
ozoneLevelTS
## Time Series:
## Start = 1927 
## End = 2016 
## Frequency = 1 
##  [1]   1.35118436   0.76053242  -1.26855735  -1.46368717  -0.97920302
##  [6]   1.50856746   1.62991681   1.83242333  -0.83968364  -1.09611566
## [11]  -2.67457473  -2.78716606  -2.97317944  -0.23495548   0.97067000
## [16]   1.23652307   2.23062331   0.35671637  -2.12028099  -2.66812477
## [21]  -0.64702795   0.99608564   1.83817742   1.89922697  -0.55488121
## [26]  -1.40387419  -3.39178762   0.05777194   3.40717183   1.31488379
## [31]  -0.12882457  -2.51580137  -3.06205664  -3.33637179  -2.66332198
## [36]  -0.67958655  -2.11660422  -2.36318997  -5.36156537  -3.03103458
## [41]  -2.28838624  -1.06438684  -1.68813570  -3.16974819  -3.65647649
## [46]  -1.25151090  -1.08431732  -0.44863234  -0.17636387  -2.64954530
## [51]  -1.28317654  -4.29289634  -3.24282341  -3.60135297  -2.57288652
## [56]  -5.00586059  -6.68548244  -4.58870210  -4.32654629  -2.29370761
## [61]  -2.26456266  -2.27184846  -2.66440082  -3.79556478  -7.65843185
## [66] -10.17433972  -4.21230497  -2.82287161  -1.36776491  -4.43997062
## [71]  -3.78323838  -5.85304107  -8.55125744  -8.06501289  -7.75975806
## [76]  -6.65633206  -6.62708203  -7.83548356  -8.84424264  -7.67352209
## [81]  -7.05582939  -4.69497353  -7.12712128  -9.58954985 -10.19222042
## [86]  -9.33224686 -10.80567444 -10.86096923 -11.57941376  -9.30284452
### creating Time Series Plot
plot(ozoneLevelTS, ylab='Change in Ozone Layer', xlab='Year', type='o', 
     main ="Time Series plot - Changes in the thickness of Ozone layer")

Observations from the plot

  1. Trend - a downward trend observed.

  2. Seasonality - no repeating patterns observed.

  3. Changing Variance - change in variance is present.

  4. Behavior - mixed of both Moving Average and Auto-Regressive.

  5. Change Point - yes, there is a change point at 1992-1993.

Task-1 : Finding best fitting model

Model Building Strategies

To determine the appropriate model for time series data, we need to follow few model building strategies:

1. Model Specification

2. Model Fitting

3. Model Diagnostics

Linear Model

# extracting time from time series.
t <- time(ozoneLevelTS)

# Model 1 - Linear Regression Model

plot(ozoneLevelTS, ylab= 'Change in Ozone Layer', xlab= 'Year', type= 'o', 
     main = "Model 1 - Fitted Linear Model ")
model1 <- lm(ozoneLevelTS ~ t)
summary(model1)
## 
## Call:
## lm(formula = ozoneLevelTS ~ t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7165 -1.6687  0.0275  1.4726  4.7940 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 213.720155  16.257158   13.15   <2e-16 ***
## t            -0.110029   0.008245  -13.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.032 on 88 degrees of freedom
## Multiple R-squared:  0.6693, Adjusted R-squared:  0.6655 
## F-statistic: 178.1 on 1 and 88 DF,  p-value: < 2.2e-16
abline(model1)

par(mfrow= c(2,2))

# Residual Analysis.
plot(y= rstudent(model1), x= as.vector(time(ozoneLevelTS)), xlab= 'Year',   
     ylab= 'Standardized Residuals', type='o', main = "Time Series Plot")

# Residual Distribution.
hist(rstudent(model1), xlab= 'Standardized Residuals', main = "Histogram")

# Q-Q plots
y = rstudent(model1)
qqnorm(y, main= "Q-Q Plot")
qqline(y, col= 2, lwd= 1, lty= 2)

# ACF of linear model.
acf(rstudent(model1), main = "ACF")

# Shapiro-Wilk test of linear model.
y = rstudent(model1)
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.98733, p-value = 0.5372

Linear Model analysis

The p-value of slope < 0.05, therefore, we conclude that the slope is significant. The Adjusted R squared value is near 67% that indicates a good model. Time-series plot of the residuals represents randomness. The histogram plot slightly captures a symmetrical distribution of data that lies between -4 and +4. The Q-Q plot appears to be nearly like a straight line with slight deviations at the ends which indicates that the data points are normally distributed.

To decide the normality of residuals, Shapiro- Wilk test is used. This test uses hypothesis testing with null hypothesis H0 : data is normally distributed. If the p-value > 0.05, we fail to reject H0. Therefore, decide the normality of the residuals. After implementing this test, we obtain a p-value = 0.5372 which is greater than 0.05. Therefore, we fail to reject H0. Moreover, we can conclude that the residuals are normally distributed.

For 95% confidence limit, all the bars in ACF plot must lie between the dashed lines. In the plot there is one significant lag at the beginning and two insignificant lags after 6 lags. This indicates that all information is not been captured by the model and few information is left in the residuals. Hence, Linear model can be considered an ideal model for the given data.

Quadratic Model

# Model 2 - Quadratic Model.

t2 = t**2
model2 = lm(ozoneLevelTS ~ t + t2)
summary(model2)
## 
## Call:
## lm(formula = ozoneLevelTS ~ t + t2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1062 -1.2846 -0.0055  1.3379  4.2325 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.733e+03  1.232e+03  -4.654 1.16e-05 ***
## t            5.924e+00  1.250e+00   4.739 8.30e-06 ***
## t2          -1.530e-03  3.170e-04  -4.827 5.87e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.815 on 87 degrees of freedom
## Multiple R-squared:  0.7391, Adjusted R-squared:  0.7331 
## F-statistic: 123.3 on 2 and 87 DF,  p-value: < 2.2e-16
plot(ts(fitted(model2)), ylab='Change in Ozone Layer', main ="Model 2 - Fitted Quadratic Curve",  
  ylim = c(min(c(fitted(model2), as.vector(ozoneLevelTS))),             
  max(c(fitted(model2), as.vector(ozoneLevelTS)))))
lines(as.vector(ozoneLevelTS),type="o")

par(mfrow= c(2,2))
# Residual Analysis.
plot(y= rstudent(model2), x= as.vector(time(ozoneLevelTS)), xlab= 'Year',    
     ylab= 'Standardized Residuals', type= 'o', main = "Time Series Plot")

# Residual Distribution.
hist(rstudent(model2), xlab= 'Standardized Residuals', main = "Histogram")

# Q-Q plots.
y = rstudent(model2)
qqnorm(y, main = "QQ Plot")
qqline(y, col= 2, lwd= 1, lty= 2)

# Shapiro-Wilk test
y = rstudent(model2)
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.98889, p-value = 0.6493
# ACF of quadratic model
acf(rstudent(model2), main = "ACF")

#### Quadratic Model analysis

The p-value of slope < 0.05, therefore, we conclude that the slope is significant. The Adjusted R squared value is near 73% that indicates a better value observed than the linear model. Time-series plot of the residuals represents randomness. The histogram plot slightly captures a symmetrical distribution of data that lies between -4 and +4. The Q-Q plot appears to be nearly like a straight line with slight deviations at the ends which indicates that the data points are normally distributed.

To decide the normality of residuals, Shapiro- Wilk test is used. This test uses hypothesis testing with null hypothesis H0 : data is normally distributed. If the p-value > 0.05, we fail to reject H0. Therefore, decide the normality of the residuals. After implementing this test, we obtain a p-value = 0.6493 which is greater than 0.05. Therefore, we fail to reject H0. Moreover, we can conclude that the residuals are normally distributed.

For 95% confidence limit, all the bars in ACF plot must lie between the dashed lines. In the plot there is one significant lag at the beginning and two insignificant lags after 6 lags. This indicates that all information is not been captured by the model and few information is left in the residuals. Hence, Linear model can be considered an ideal model for the given data.

Seasonal Model

# Creating ACF plot for the original time series to determine its frequency for the Seasonal Model.
acf(ozoneLevelTS, main = "ACF of Original Time Series")

By looking at the above ACF plot, we can observe a presence of seasonality in the given data. Using the plot we can determine the frequency as 5.

newTS<- ts(dataset$Ozone_Level, start = 1927, end = 2016, frequency = 5)
plot(newTS, ylab='Change in Ozone Layer', xlab='Year', type='o', main ="TS Plot for Change in Ozone Layer")

# Seasonal Trend
trend = season(newTS)
t1 <- time(newTS)

# Model 3 - Seasonal Model
model3= lm(newTS ~ t1)
summary(model3)
## 
## Call:
## lm(formula = newTS ~ t1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9301 -1.7836  0.5112  2.3642  7.1411 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 32.501851  12.373926   2.627  0.00892 **
## t1          -0.018076   0.006276  -2.880  0.00417 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.413 on 444 degrees of freedom
## Multiple R-squared:  0.01834,    Adjusted R-squared:  0.01613 
## F-statistic: 8.296 on 1 and 444 DF,  p-value: 0.004165
plot(ts(fitted(model3)), ylab='Change in Ozone Layer', main ="Model 3 - Fitted Seaonal Model.",    
     ylim = c(min(c(fitted(model3), as.vector(newTS))),      
              max(c(fitted(model3), as.vector(newTS)))), col ="red")
lines(as.vector(newTS),type="o")

par(mfrow= c(2,2))
# Residual Analysis.
plot(y= rstudent(model3), x= as.vector(time(newTS)), xlab= 'Year',  
     ylab= 'Standardized Residuals', type= 'o', main ="Time Series Plot.")
points(y= rstudent(model3), x= as.vector(time(newTS)), pch= as.vector(season(newTS)))

# Residual Distribution
hist(rstudent(model3),xlab='Standardized Residuals', main ="Histogram")

# Q-Q Plot
y = rstudent(model3)
qqnorm(y, main ="QQ Plot")
qqline(y, col =2, lwd =1, lty =2)


# Shapiro-Wilk test of seasonal model ####
y = rstudent(model3)
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.96351, p-value = 4.438e-09
# ACF
acf(rstudent(model3), main ="ACF")

Seasonal Model Analysis

The p-value of slope < 0.05, therefore, we conclude that the slope is significant. The Adjusted R squared value is near 67% that indicates a good model. Time-series plot of the residuals represents randomness. The histogram plot slightly captures a symmetrical distribution of data that lies between -4 and +4. The Q-Q plot appears to be nearly like a straight line with slight deviations at the ends which indicates that the data points are normally distributed.

To decide the normality of residuals, Shapiro- Wilk test is used. This test uses hypothesis testing with null hypothesis H0 : data is normally distributed. If the p-value > 0.05, we fail to reject H0. Therefore, decide the normality of the residuals. After implementing this test, we obtain a p-value = 0.5372 which is greater than 0.05. Therefore, we fail to reject H0. Moreover, we can conclude that the residuals are normally distributed.

For 95% confidence limit, all the bars in ACF plot must lie between the dashed lines. In the plot there is one significant lag at the beginning and two insignificant lags after 6 lags. This indicates that all information is not been captured by the model and few information is left in the residuals. Hence, Linear model can be considered an ideal model for the given data.

Cosine Model

# Model 4 - Cosine Model
harm <- harmonic(newTS, 1)
modelData <- data.frame(newTS, harm)

model4 <- lm(newTS ~ cos.2.pi.t. + sin.2.pi.t. , data = modelData)
summary(model4)
## 
## Call:
## lm(formula = newTS ~ cos.2.pi.t. + sin.2.pi.t., data = modelData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0964 -1.8785  0.4434  2.3783  6.8902 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.1359     0.1627  -19.28   <2e-16 ***
## cos.2.pi.t.   0.1379     0.2298    0.60   0.5487    
## sin.2.pi.t.   0.4008     0.2303    1.74   0.0825 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.435 on 443 degrees of freedom
## Multiple R-squared:  0.007592,   Adjusted R-squared:  0.003111 
## F-statistic: 1.694 on 2 and 443 DF,  p-value: 0.1849
plot(ts(fitted(model4)), ylab='y', main ="Fitted cosine wave ",   
    ylim = c(min(c(fitted(model4), as.vector(newTS))),          
             max(c(fitted(model4), as.vector(newTS)))), 
    col ="red")
lines(as.vector(newTS),type="o")

par(mfrow= c(2,2))
# Residual Analysis.
plot(y= rstudent(model4), x= as.vector(time(newTS)), xlab= 'Time',   
     ylab='Standardized Residuals', type='o', main ="Time Series Plot")
points(y=rstudent(model4),x=as.vector(time(newTS)), pch=as.vector(season(newTS)))

# Residual distribution
hist(rstudent(model4),xlab='Standardized Residuals', main ="Histogram")

# Q-Q plots
y = rstudent(model4)
qqnorm(y, main ="QQ Plot")
qqline(y, col =2, lwd =1, lty =2)

# Shapiro-Wilk test
y = rstudent(model4)
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.95788, p-value = 5.518e-10
# ACF
acf(rstudent(model4), main ="ACF")

Cosine Model Analysis

The p-value of slope < 0.05, therefore, we conclude that the slope is significant. The Adjusted R squared value is near 67% that indicates a good model. Time-series plot of the residuals represents randomness. The histogram plot slightly captures a symmetrical distribution of data that lies between -4 and +4. The Q-Q plot appears to be nearly like a straight line with slight deviations at the ends which indicates that the data points are normally distributed.

To decide the normality of residuals, Shapiro- Wilk test is used. This test uses hypothesis testing with null hypothesis H0 : data is normally distributed. If the p-value > 0.05, we fail to reject H0. Therefore, decide the normality of the residuals. After implementing this test, we obtain a p-value = 0.5372 which is greater than 0.05. Therefore, we fail to reject H0. Moreover, we can conclude that the residuals are normally distributed.

For 95% confidence limit, all the bars in ACF plot must lie between the dashed lines. In the plot there is one significant lag at the beginning and two insignificant lags after 6 lags. This indicates that all information is not been captured by the model and few information is left in the residuals. Hence, Linear model can be considered an ideal model for the given data.

Prediction of Yearly Changes in Ozone Layer for next 5 years.

# Predicting future observations for the linear model
h <-5
forecastTS <- time(ozoneLevelTS)
futureTimes <- data.frame(t = seq(2017,2021,1))
frcModel1 <- predict(model1, newdata = futureTimes, interval ="prediction")
frcModel1
##         fit       lwr       upr
## 1 -8.208590 -12.33732 -4.079864
## 2 -8.318619 -12.45033 -4.186902
## 3 -8.428648 -12.56342 -4.293878
## 4 -8.538677 -12.67656 -4.400792
## 5 -8.648706 -12.78977 -4.507643
### plotting the forecast ###
plot(ozoneLevelTS, xlim= c(1927,2021), ylim = c(-20.00, 5), ylab ="Ozone change series",     
     main ="Forecasts from the linear model fitted to the ozone change series.")

lines(ts(as.vector(frcModel1[,3]), start =2017), col="blue", type="l")
lines(ts(as.vector(frcModel1[,1]), start =2017), col="red", type="l")
lines(ts(as.vector(frcModel1[,2]), start =2017), col="blue", type="l")
legend("bottomleft", lty=1, pch=1, col=c("black","blue","red"), text.width =18,      
       c("Data","5% forecast limits","Forecasts"))

Task 2

To generate the ARIMA model, we would require to remove the patterns created by ACF and PACF plots and this process of removing trend from the series is known as differencing. The process converts a non stationary series to a stationary series.

An apparent change in the variation of the data has been observed. For the building of ARIMA model, we need to get rid of the trends and variance.

Box-Cox Transformation To stabilize the variance within the data, we need to apply transformation. For this process, we need to make the entire data positive as the transformation only deals with positive data values. To encounter this issue, we will add a positive constant to the data values.

plot(ozoneLevelTS,ylab=' Change', xlab='Year',type='o',  main ="Time Series plot for change in ozone")

par(mfrow=c(1,2))
acf(ozoneLevelTS,main="ACF plot for the change in ozone series") 
pacf(ozoneLevelTS,main="PACF plot for the change in ozone series")

# Applying  Box-cox transformation  to stabilize the observed variance 

boxCox = BoxCox.ar(ozoneLevelTS + abs(min(ozoneLevelTS))+1)

boxCox$ci
## [1] 0.9 1.5
lambda <- boxCox$lambda[which(max(boxCox$loglike) == boxCox$loglike)]
lambda 
## [1] 1.2
ozoneShiftTS = ozoneLevelTS + abs(min(ozoneLevelTS))+1

boxCox.ozoneShiftTS = (ozoneShiftTS^lambda-1)/ozoneShiftTS

plot(boxCox.ozoneShiftTS, type='o', ylab= 'Simulated series', main= "Time series plot of the BC-transformed simulated series.")

# Checking Q-Q plots for the original data
qqnorm(ozoneLevelTS)
qqline(ozoneLevelTS, col= 2)


# Checking Q-Q plots for the transformed data
qqnorm(boxCox.ozoneShiftTS)
qqline(boxCox.ozoneShiftTS, col= 2)

A slight improvement in the variance was observed after performing Box-Cox transformation but the normality got worse. The comparison of the Q-Q Plots of both transformed data and original data shows that the normality of the original data was aligned better than that of the transformed data. As there was no significant change observed due transformation, we choose to proceed with the original data for differencing.

We will use Dickey Fuller Unit-Root test to test the null hypothesis for the process to be non-stationary at first but then changes to stationary after first differencing. If p-value > 0.05, then we can say that the series is stationary.

After getting a stationary series, we plot the ACF and PACF of differenced series and propose set of ARIMA Models.

# Dickey-Fuller unit-root Test

adf.test(ozoneLevelTS, alternative= c("stationary"))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ozoneLevelTS
## Dickey-Fuller = -3.2376, Lag order = 4, p-value = 0.0867
## alternative hypothesis: stationary
# Implementing Differencing

diff.ozoneLevelTS = diff(ozoneLevelTS)
plot(diff.ozoneLevelTS, type='o', ylab='ozonechange', main= "Time series plot of the first transformed simulated series.")                        

# implementing Dickey - Fuller test  and Phillip-Perron test to differenced series

adf.test(diff.ozoneLevelTS, alternative= c("stationary"))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff.ozoneLevelTS
## Dickey-Fuller = -7.1568, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
pp.test(diff.ozoneLevelTS, alternative= c("stationary"))
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff.ozoneLevelTS
## Dickey-Fuller Z(alpha) = -70.244, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

On application of Dickey Fuller unit-root test, we observed that the (p-value = 0.0867) > 0.05, therefore, we can conclude that the series is non-stationary. After the first differencing, no trend was observed in the series. To determine if the series is stationary or not, we applied Dicker Fuller test on it and observed the (p-value = 0.01) < 0.05, therefore, we can conclude that the series is stationary. We also performed Phillip-Peppron unit-root test to double check if the series is stationary. After the test we observed (p = 0.01) < 0.05, which confirms that the series is stationary.

Auto Regressive Integrated Moving Average, ARIMA(p, q, d) is a combination of AR(p) and MA(q) and differencing(d). The value of p and q is obtained by PACF and ACF respectively. Differencing d will be equal to 1 as we did differencing only once.

We required ACF and PACF to effectively identify AR(p) and MA(q) models. We use EACF(Extended-ACF) to identify ARMA(p,q) models and a table with elements in the ith row and jth column is displayed. In the EACF table there will be a vertex of zeros also known as the vertex of top left zeroes. This vertex helps in identifying the p and q values of the ARMA(p,q) model.

Bayesian Information Criterion (BIC) is a model specification which is used to identify ARMA(p,q) models. The BIC table consists of rows that exhibit p and q values of the ARMA(p,q) model where the cells of these variables are shaded. Lower BIC models are situated in higher rows and are darker in shade.

# Plotting the ACF and PACF for first difference.

par(mfrow=c(1,2)) 

acf(diff.ozoneLevelTS, main = "ACF - First Difference")
pacf(diff.ozoneLevelTS, main = "PACF - First Difference")

par(mfrow=c(1,1)) 

eacf(diff.ozoneLevelTS)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o x o o o x o o x o  o  o  o 
## 1 x o x o o o o o o x o  o  o  o 
## 2 x o x o o o x o o x o  o  o  o 
## 3 x o x o o x o o o o o  o  o  o 
## 4 x o o x o x o o o o o  o  o  o 
## 5 x x x x o x o o o o o  o  o  o 
## 6 o o o x x o o o o o o  o  o  o 
## 7 o o o x o o o o o o o  o  o  o
ozoneRes=armasubsets(y= diff.ozoneLevelTS, nar=5, nma=5, y.name='test', ar.method='ols')
plot(ozoneRes)

The observations for p and q from ACF and PACF are :

p = 1, 2, 3

q = 1, 2

The proposed models are ARIMA(1,1,1), ARIMA(1,1,2), ARIMA(2,1,1), ARIMA(2,1,2), ARIMA(3,1,1), ARIMA(3,1,2).

The observations for p and q from EACF are :

p = 1, 2

q = 3, 4

The proposed models are ARIMA(1,1,3), ARIMA(1,1,4), ARIMA(2,1,3), ARIMA(2,1,4).

The observations from BIC table are :

p = 3

q = 2

Therefore, the proposed model is ARIMA(3,1,2).

Conclusion

Task 1

As per the experiment and observations, we can conclude that the most suitable model should be linear model in comparison to the quadratic, seasonal and cyclic models as it captures most of the data and just a small amount of data is left in the residuals.

Task 2

The proposed set of ARIMA(p,q,d) models are ARIMA(1,1,1), ARIMA(1,1,2), ARIMA(2,1,1), ARIMA(2,1,2), ARIMA(3,1,1), ARIMA(3,1,2), ARIMA(1,1,3), ARIMA(1,1,4), ARIMA(2,1,3), ARIMA(2,1,4). There is a strong evidence obtained from ACF-PACF and BIC table that confirmed the possibility ARIMA(3,1,2) model.