Import Data Air Passenger and Library tseries, lmtest, forecast, and ggplot2.
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(forecast)
library(ggplot2)
print(AirPassengers)
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
AP <- AirPassengers
ts.plot(AP,xlab="Time",ylab="Air Passenger")
acf(AP,main="Air Passenger")
pacf(AP,main="Air Passenger")
dugaan1 = arima(AP, order = c(1,0,0))
dugaan2 = arima(AP, order = c(0,0,1))
dugaan3 = arima(AP, order = c(1,0,1))
coeftest(dugaan1)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.964632 0.021359 45.1633 < 2.2e-16 ***
## intercept 278.464942 67.114062 4.1491 3.337e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(dugaan2)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.964211 0.021393 45.070 < 2.2e-16 ***
## intercept 280.646447 10.578799 26.529 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(dugaan3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.937312 0.030163 31.0748 < 2.2e-16 ***
## ma1 0.426400 0.091068 4.6822 2.838e-06 ***
## intercept 281.542574 53.613469 5.2513 1.510e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res1=dugaan1$residuals
res2=dugaan2$residuals
res3=dugaan3$residuals
Box.test(res1,lag=6)
##
## Box-Pierce test
##
## data: res1
## X-squared = 33.417, df = 6, p-value = 8.715e-06
Box.test(res2,lag=6)
##
## Box-Pierce test
##
## data: res2
## X-squared = 436.44, df = 6, p-value < 2.2e-16
Box.test(res3,lag=6)
##
## Box-Pierce test
##
## data: res3
## X-squared = 10.475, df = 6, p-value = 0.106
accuracy(dugaan1)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.944642 33.44577 25.7074 -0.5877058 9.116557 0.9940936 0.3076485
accuracy(dugaan2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.4138396 64.84889 52.86729 -10.70417 23.16898 2.044354 0.7748281
accuracy(dugaan3)
## ME RMSE MAE MPE MAPE MASE
## Training set 1.339968 31.12125 24.46295 -0.7897265 8.913893 0.9459713
## ACF1
## Training set -0.01710536
checkresiduals(dugaan1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 313.23, df = 23, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 24
checkresiduals(dugaan2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1) with non-zero mean
## Q* = 1282.6, df = 23, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 24
checkresiduals(dugaan3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 214.03, df = 22, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 24
fitting1 = fitted(dugaan1)
fitting2 = fitted(dugaan2)
fitting3 = fitted(dugaan3)
prediksi1 = forecast(dugaan1, h=20)
prediksi2 = forecast(dugaan2, h=20)
prediksi3 = forecast(dugaan3, h=20)
plot(AP,col='black',xlab="Time",ylab="Air Passenger")
lines(fitting1,col='blue')
lines(fitting2,col='green')
lines(fitting3,col='magenta')
legend("topleft", # Position
legend = c("Data Actual","ARIMA(1,0,0)", "ARIMA(0,0,1)","ARIMA(1,0,1)"), # Legend texts
lty = c(1, 1, 1,1,2), # Line types
col = c('black','blue','green','magenta','blue'), # Line colors
lwd = 2)
plot(prediksi1,col='blue',lty=1)
plot(prediksi2,col='green',lty=1)
plot(prediksi3,col='magenta',lty=1)
arimaAP <- auto.arima(AP)
arimaAP
## Series: AP
## ARIMA(2,1,1)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ma1
## 0.5960 0.2143 -0.9819
## s.e. 0.0888 0.0880 0.0292
##
## sigma^2 = 132.3: log likelihood = -504.92
## AIC=1017.85 AICc=1018.17 BIC=1029.35
forecastAP <- forecast(arimaAP, level = c(95), h = 36)
autoplot(forecastAP,xlab="Time",ylab="Air Passenger")
Note: According to the results using function auto.arima, then the three possible models ARIMA that we have created are nothing accepted for forecasting (ARIMA(2,1,1) is the model).