Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(readxl)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TTR)
library(imputeTS)
library(tseries)
## 
## Attaching package: 'tseries'
## The following object is masked from 'package:imputeTS':
## 
##     na.remove
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(graphics)
library(TSA)
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ tibble  3.1.8     ✔ purrr   0.3.5
## ✔ tidyr   1.2.1     ✔ stringr 1.5.0
## ✔ readr   2.1.3     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ readr::spec()   masks TSA::spec()
library(lubridate)
## Loading required package: timechange
## 
## Attaching package: 'lubridate'
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(ggfortify)
## Registered S3 methods overwritten by 'ggfortify':
##   method                 from    
##   autoplot.Arima         forecast
##   autoplot.acf           forecast
##   autoplot.ar            forecast
##   autoplot.bats          forecast
##   autoplot.decomposed.ts forecast
##   autoplot.ets           forecast
##   autoplot.forecast      forecast
##   autoplot.stl           forecast
##   autoplot.ts            forecast
##   fitted.ar              forecast
##   fortify.ts             forecast
##   residuals.ar           forecast
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(graphics)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following object is masked from 'package:imputeTS':
## 
##     na.locf
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(stats)
library(MASS)
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(fpp2)
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.4 ──
## ✔ fma       2.4     ✔ expsmooth 2.3
library(hrbrthemes)
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
##       Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
##       if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
setwd("D:/SATRIA DATA")
mpdw <- read_xlsx("data pm.xlsx")
mpdw <- mpdw[,1:2]
mpdw$`Data PM 2.5`<- as.numeric(mpdw$`Data PM 2.5`)
ts.mpdw <- ts(mpdw)
kableExtra::kable(head(mpdw) ,caption = 'Subset Data CO DKI Jakarta 2021-2022')
Subset Data CO DKI Jakarta 2021-2022
t Data PM 2.5
2021-01-01 16.28
2021-01-02 19.66
2021-01-03 16.13
2021-01-04 15.28
2021-01-05 21.90
2021-01-06 17.71 
view(mpdw)
# Train-test Split (80/20)
train <- mpdw[1:363,]
test <- mpdw[364:454,]
datampdw <- ts(mpdw$`Data PM 2.5`)
ts.plot(datampdw, main = "Kadar PM 2.5 DKI Jakarta Tahun 2021 - 2022", ylab = "Harga", lwd = 1.5)
points(datampdw)
legend("topleft", c("Premium"), cex = 0.7,
       col = c("black", "red", "blue"), lty = 1)

plot(x = mpdw$t,
     y = mpdw$`Data PM 2.5`,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Kadar PM 2.5 DKI Jakarta",
     main = "Time Series Plot Kadar PM 2.5 DKI Jakarta")

# Time-Series Object
data.ts <- ts(mpdw[,2], start = 2021, end = 2022, frequency = 91)
train.ts <- ts(train)
test.ts <- ts(test)
acf(mpdw$`Data PM 2.5`, lag = 35, start= 1)
## Warning in plot.window(...): "start" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "start" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "start" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "start" is not a
## graphical parameter
## Warning in box(...): "start" is not a graphical parameter
## Warning in title(...): "start" is not a graphical parameter

pacf(mpdw$`Data PM 2.5`, lag =35)

adf.test(data.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.ts
## Dickey-Fuller = -3.8522, Lag order = 4, p-value = 0.01982
## alternative hypothesis: stationary
data.diff <- diff(data.ts, differences = 1)
plot(data.diff)

adf.test(data.diff)
## Warning in adf.test(data.diff): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff
## Dickey-Fuller = -6.7856, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
acf(data.diff, lag.max = 50, start = 1)
## Warning in plot.window(...): "start" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "start" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "start" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "start" is not a
## graphical parameter
## Warning in box(...): "start" is not a graphical parameter
## Warning in title(...): "start" is not a graphical parameter

pacf(data.diff, lag.max =50)

eacf(data.diff)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x 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 x o o o o o o o o  o  o  o 
## 4 x x x o o o x o o o o  o  o  o 
## 5 x x x x x o o o o o o  o  o  o 
## 6 x x x o x x o o o o o  o  o  o 
## 7 x o o x x o o o o o o  o  o  o
model1 <- Arima(data.ts,order = c(0,1,1),method ="ML")
model2 <- Arima(data.ts,order = c(1,1,1),method ="ML")
model3 <- Arima(data.ts,order = c(3,1,3),method ="ML")
model4 <- Arima(data.ts,order = c(5,1,5),method = "ML")
summary(model1)
## Series: data.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.4122
## s.e.   0.1339
## 
## sigma^2 = 105.6:  log likelihood = -340.74
## AIC=685.49   AICc=685.62   BIC=690.51
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.01150391 10.16553 7.597971 -13.14697 37.07492 28.14063
##                    ACF1
## Training set 0.03853498
summary(model2)
## Series: data.ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.4903  -0.9448
## s.e.  0.1198   0.0605
## 
## sigma^2 = 95.29:  log likelihood = -336.09
## AIC=678.17   AICc=678.45   BIC=685.7
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.7618441 9.601392 7.158062 -11.80964 35.16014 26.51134
##                     ACF1
## Training set -0.05827018
summary(model3)
## Series: data.ts 
## ARIMA(3,1,3) 
## 
## Coefficients:
##           ar1      ar2     ar3      ma1     ma2      ma3
##       -0.2763  -0.2593  0.4890  -0.2303  0.0641  -0.7492
## s.e.   0.2202   0.2019  0.1319   0.2235  0.2149   0.1674
## 
## sigma^2 = 93.58:  log likelihood = -333.38
## AIC=680.75   AICc=682.1   BIC=698.33
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
## Training set 0.7856919 9.298257 6.831124 -10.71364 33.13944 25.30046 0.01450311
summary(model4)
## Series: data.ts 
## ARIMA(5,1,5) 
## 
## Coefficients:
##           ar1     ar2     ar3      ar4      ar5      ma1      ma2      ma3
##       -0.3936  0.3809  0.3531  -0.2220  -0.5181  -0.1121  -0.6237  -0.3249
## s.e.   0.2521  0.2491  0.2474   0.2136   0.1412   0.2759   0.2883   0.3330
##          ma4     ma5
##       0.1190  0.3765
## s.e.  0.2347  0.2514
## 
## sigma^2 = 91.86:  log likelihood = -330.3
## AIC=682.59   AICc=685.93   BIC=710.21
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.1424416 8.993033 6.566148 -11.50557 32.79353 24.31907
##                     ACF1
## Training set -0.01047346
coeftest(model1)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)   
## ma1 -0.41221    0.13394 -3.0776 0.002086 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(model2)
## 
## z test of coefficients:
## 
##     Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.49030    0.11976   4.0939 4.243e-05 ***
## ma1 -0.94475    0.06052 -15.6104 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(model3)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.276263   0.220188 -1.2547 0.2095981    
## ar2 -0.259299   0.201943 -1.2840 0.1991361    
## ar3  0.489006   0.131860  3.7085 0.0002085 ***
## ma1 -0.230286   0.223539 -1.0302 0.3029238    
## ma2  0.064072   0.214943  0.2981 0.7656342    
## ma3 -0.749195   0.167424 -4.4748 7.647e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model <- c("ARIMA (0,1,1)","ARIMA (1,1,1)","ARIMA (3,1,3)","ARIMA (5,1,5)")
AIC <- c(model1$aic,model2$aic,model3$aic,model4$aic)
BIC <- c(model1$bic,model2$bic,model3$bic,model4$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)
Model AIC BIC
ARIMA (0,1,1) 685.4861 690.5078
ARIMA (1,1,1) 678.1700 685.7026
ARIMA (3,1,3) 680.7514 698.3274
ARIMA (5,1,5) 682.5923 710.2117 
model5 <- auto.arima(data.ts)
summary(model5)
## Series: data.ts 
## ARIMA(0,1,4) 
## 
## Coefficients:
##           ma1      ma2      ma3      ma4
##       -0.4993  -0.0302  -0.1514  -0.2159
## s.e.   0.1030   0.1115   0.1103   0.0916
## 
## sigma^2 = 92.4:  log likelihood = -333.73
## AIC=677.45   AICc=678.16   BIC=690
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.7591565 9.347633 6.948725 -11.18935 34.28335 25.73602
##                     ACF1
## Training set -0.01447111
paste("Model yang terbaik adalah model",Akurasi$Model[which.min(Akurasi[,"AIC"])])
## [1] "Model yang terbaik adalah model ARIMA (1,1,1)"
over1 <- Arima(train.ts[,2], order = c(0,1,2))
summary(over1)
## Series: train.ts[, 2] 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.5385  -0.2442
## s.e.   0.0492   0.0506
## 
## sigma^2 = 88.03:  log likelihood = -1323.53
## AIC=2653.05   AICc=2653.12   BIC=2664.73
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE    MAPE     MASE       ACF1
## Training set 0.1201814 9.343553 7.240042 -10.59197 28.9449 0.918639 0.01814187
coeftest(over1)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ma1 -0.538452   0.049192 -10.9459 < 2.2e-16 ***
## ma2 -0.244225   0.050583  -4.8283 1.377e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
over2 <- Arima(train.ts[,2], order= c(2,1,0))
summary(over2)
## Series: train.ts[, 2] 
## ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.3875  -0.2380
## s.e.   0.0511   0.0511
## 
## sigma^2 = 97.58:  log likelihood = -1341.87
## AIC=2689.74   AICc=2689.81   BIC=2701.41
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.02579747 9.837491 7.518211 -9.446969 29.18456 0.9539339
##                    ACF1
## Training set -0.0380757
coeftest(over2)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.387484   0.051088 -7.5847 3.332e-14 ***
## ar2 -0.238034   0.051114 -4.6570 3.209e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
over3 <- Arima(train.ts[,2], order = c(1,1,2))
summary(over3)
## Series: train.ts[, 2] 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.4322  -0.9570  0.0480
## s.e.  0.2015   0.2134  0.1648
## 
## sigma^2 = 87.42:  log likelihood = -1321.83
## AIC=2651.66   AICc=2651.77   BIC=2667.22
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.1489196 9.298027 7.163221 -10.82299 28.8532 0.9088916
##                      ACF1
## Training set -0.001188257
over4 <- Arima(train.ts[,2], order = c(2,1,1))
summary(over4)
## Series: train.ts[, 2] 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.3782  0.0167  -0.9021
## s.e.  0.0694  0.0637   0.0453
## 
## sigma^2 = 87.42:  log likelihood = -1321.84
## AIC=2651.67   AICc=2651.78   BIC=2667.24
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1483485 9.298259 7.164343 -10.81358 28.85456 0.9090339
##                      ACF1
## Training set -0.001964789
train.ts <- ts(train$`Data PM 2.5`,start = 1)
ramalan <- forecast::forecast(Arima(train.ts, order=c(1,1,2),method="ML",include.drift = TRUE), h=10) 
data.ramalan <- ramalan$mean
plot(ramalan,lwd=2)

accuracy(ramalan)
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.01869328 9.297048 7.165763 -11.3473 28.99974 0.9092141
##                       ACF1
## Training set -0.0008830543