MONTHLY INDUSTRIAL ELECTRICITY PRODUCTION FROM 1985 TO 2018
A Time Series Analysis and Prediction
Anh Phan (s3258110) and Neeraj Sehrawat (s3711712)
Last updated: 22 June, 2019
Read the dataset
Load the dataset
setwd("~/Time Series/Project")
electric <- read.csv("Electric_Production.csv")
Data Preprocessing
Rename the columns
names(electric) <- c("date", "production")
electric.ts <- matrix(electric$production, nrow = 397, ncol = 1)
electric.ts <- as.vector(t(electric.ts))
electric.ts <- ts(electric.ts,start=c(1985,1), end=c(2018,1), frequency=12)
Plot the time series
plot(electric.ts,type='l',ylab='Electricity Production')
points(y=electric.ts,x=time(electric.ts), pch=as.vector(season(electric.ts)))
Investigate the correlation
y1 = electric.ts
x1 = zlag(electric.ts)
index = 2: length(x1)
cor(y1[index], x1[index])
## [1] 0.8717309
plot(y = electric.ts, x = zlag(electric.ts), ylab = 'Electric production', xlab = 'Previous Year electric production')
ACF/PACF
par(mfrow=c(1,2))
acf(electric.ts,lag.max=36)
pacf(electric.ts, lag.max=36)
Residual approach - When D = 1
m0.electric = arima(electric.ts,order=c(0,0,0),seasonal=list(order=c(0,1,0), period=12))
res.m0_electric = residuals(m0.electric)
par(mfrow=c(3,1))
plot(res.m0_electric, xlab='Time',ylab='Residuals',main="Time series plot of the residuals")
acf(res.m0_electric, lag.max = 36, main = "The sample ACF of the residuals")
pacf(res.m0_electric, lag.max = 36, main = "The sample PACF of the residuals")
Residual approach - When D = 1, P = 0, Q = 2
m1.electric = arima(electric.ts,order=c(0,0,0),seasonal=list(order=c(0,1,2), period=12))
res.m1_electric = residuals(m1.electric)
par(mfrow=c(3,1))
plot(res.m1_electric,xlab='Time',ylab='Residuals',main="Time series plot of the residuals")
acf(res.m1_electric, lag.max = 36, main = "The sample ACF of the residuals")
pacf(res.m1_electric, lag.max = 36, main = "The sample PACF of the residuals")
Residual approach - When D = 1, P = 3, Q = 2
m2.electric = arima(electric.ts,order=c(0,0,0),seasonal=list(order=c(3,1,2), period=12))
res.m2_electric = residuals(m2.electric)
par(mfrow=c(3,1))
plot(res.m1_electric,xlab='Time',ylab='Residuals',main="Time series plot of the residuals")
acf(res.m2_electric, lag.max = 36, main = "The sample ACF of the residuals")
pacf(res.m2_electric, lag.max = 36, main = "The sample PACF of the residuals")
Residual approach - When d = 1, p = 4, q = 3
m5.electric = arima(log(electric.ts),order=c(4,1,3),seasonal=list(order=c(3,1,2), period=12))
par(mfrow=c(3,1))
res.m5_electric = residuals(m5.electric)
plot(res.m5_electric,xlab='Time',ylab='Residuals',main="Time series plot of the residuals")
acf(res.m5_electric, lag.max = 36, main = "The sample ACF of the residuals")
pacf(res.m5_electric, lag.max = 36, main = "The sample PACF of the residuals")
ADF test
order = ar(diff(res.m5_electric))$order
adfTest(res.m5_electric, lags = order, title = NULL,
description = NULL)@test$p.value
##
## 0.01
EACF
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o o o o o o o o o o o o o o
## 3 o x o o o o o o o o o o o o
## 4 x x x x o o o o o o o o o o
## 5 x x x x o o o o o o o o o o
## 6 x o x o x x o o o o o o o o
## 7 o x x o x x o o o o o o o o
Our candidate models
** SARIMA(4,1,3)*(3,1,2)_12**
** SARIMA(0,1,0)*(3,1,2)_12**
** SARIMA(0,1,1)*(3,1,2)_12**
** SARIMA(1,1,1)*(3,1,2)_12**
Model Fitting - SARIMA(4,1,3)x(3,1,2)_12
m5_413.electric = arima(log(electric.ts),order=c(4,1,3),seasonal=list(order=c(3,1,2), period=12))
coeftest(m5_413.electric)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.160909 NA NA NA
## ar2 0.121409 NA NA NA
## ar3 0.360098 0.011290 31.895 < 2.2e-16 ***
## ar4 0.061895 NA NA NA
## ma1 0.774346 0.039270 19.718 < 2.2e-16 ***
## ma2 -0.849614 NA NA NA
## ma3 -0.666124 NA NA NA
## sar1 -0.181698 0.013009 -13.967 < 2.2e-16 ***
## sar2 -0.260763 NA NA NA
## sar3 -0.119683 NA NA NA
## sma1 -0.632295 NA NA NA
## sma2 -0.063948 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model Fitting - SARIMA(0,1,0)x(3,1,2)_12
m5_010.electric = arima(log(electric.ts),order=c(0,1,0),seasonal=list(order=c(3,1,2), period=12))
coeftest(m5_010.electric)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## sar1 -0.237361 NA NA NA
## sar2 -0.190652 0.064105 -2.974 0.002939 **
## sar3 -0.062712 NA NA NA
## sma1 -0.553206 NA NA NA
## sma2 -0.190620 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model Fitting - SARIMA(0,1,1)x(3,1,2)_12
m5_011.electric = arima(log(electric.ts),order=c(0,1,1),seasonal=list(order=c(3,1,2), period=12))
coeftest(m5_011.electric )
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.448606 0.069355 -6.4683 9.912e-11 ***
## sar1 0.404387 0.447229 0.9042 0.365886
## sar2 -0.261358 0.064010 -4.0831 4.444e-05 ***
## sar3 -0.018447 0.124764 -0.1479 0.882456
## sma1 -1.183910 0.444580 -2.6630 0.007745 **
## sma2 0.374207 0.348297 1.0744 0.282647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model Fitting - SARIMA(1,1,1)x(3,1,2)_12
m5_111.electric = arima(log(electric.ts),order=c(1,1,1),seasonal=list(order=c(3,1,2), period=12))
coeftest(m5_111.electric )
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.502650 0.054785 9.1749 < 2.2e-16 ***
## ma1 -0.925915 0.023949 -38.6624 < 2.2e-16 ***
## sar1 0.400445 0.499027 0.8025 0.4222913
## sar2 -0.245524 0.067396 -3.6430 0.0002695 ***
## sar3 0.022966 0.135656 0.1693 0.8655642
## sma1 -1.176522 0.494813 -2.3777 0.0174205 *
## sma2 0.358720 0.376902 0.9518 0.3412191
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual Diagnostics - SARIMA(4,1,3)x(3,1,2)_12
residual.analysis(model = m5_413.electric)
##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.99425, p-value = 0.1412
##
##
## Box-Ljung test
##
## data: (res.model)
## X-squared = 1.8619, df = 6, p-value = 0.932
Residual Diagnostics - SARIMA(0,1,0)x(3,1,2)_12
residual.analysis(model = m5_010.electric)
##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.99278, p-value = 0.05282
##
##
## Box-Ljung test
##
## data: (res.model)
## X-squared = 34.261, df = 6, p-value = 5.99e-06
Residual Diagnostics - SARIMA(0,1,1)x(3,1,2)_12
residual.analysis(model = m5_011.electric)
##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.99495, p-value = 0.2211
##
##
## Box-Ljung test
##
## data: (res.model)
## X-squared = 26.028, df = 6, p-value = 0.0002199
Residual Diagnostics - SARIMA(1,1,1)x(3,1,2)_12
residual.analysis(model = m5_111.electric)
##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.99456, p-value = 0.1725
##
##
## Box-Ljung test
##
## data: (res.model)
## X-squared = 4.4203, df = 6, p-value = 0.62
AIC/BIC checking
sc.AIC=AIC(m5_111.electric, m5_413.electric)
sc.BIC=AIC(m5_111.electric, m5_413.electric, k = log(length(electric.ts)))
sort.score(sc.AIC, score = "aic")
sort.score(sc.BIC, score = "aic")
Overfitting
m5_112.electric = arima(log(electric.ts),order=c(1,1,2),seasonal=list(order=c(3,1,2), period=12))
coeftest(m5_112.electric)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.470568 0.125232 3.7576 0.0001716 ***
## ma1 -0.885216 0.142029 -6.2326 4.586e-10 ***
## ma2 -0.034943 0.119518 -0.2924 0.7700044
## sar1 0.381860 0.538944 0.7085 0.4786132
## sar2 -0.243940 0.067793 -3.5983 0.0003203 ***
## sar3 0.022348 0.141958 0.1574 0.8749072
## sma1 -1.157844 0.535168 -2.1635 0.0305014 *
## sma2 0.342781 0.408321 0.8395 0.4011947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Overfitting - Continue
m5_211.electric = arima(log(electric.ts),order=c(2,1,1),seasonal=list(order=c(3,1,2), period=12), method = "CSS")
coeftest(m5_211.electric)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.428928 0.061406 6.9851 2.847e-12 ***
## ar2 -0.091170 0.052225 -1.7457 0.080864 .
## ma1 -0.835403 0.038923 -21.4631 < 2.2e-16 ***
## sar1 -0.411551 0.226570 -1.8164 0.069303 .
## sar2 -0.324247 0.074182 -4.3710 1.237e-05 ***
## sar3 -0.189186 0.066409 -2.8488 0.004388 **
## sma1 -0.374970 0.226860 -1.6529 0.098356 .
## sma2 -0.197498 0.192407 -1.0265 0.304673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Forecasting
m5_111.electric = Arima(log(electric.ts),order=c(1,1,1),seasonal=list(order=c(3,1,2), period=12))
preds1 = forecast(m5_111.electric, h = 10)
plot(preds1)
