library(readxl)
library(forecast)
library(TTR)
library(imputeTS)
library(tseries)
library(ggplot2)
library(dplyr)
library(graphics)
library(TSA)
library(tidyverse)
library(lubridate)
library(gridExtra)
library(ggfortify)
library(cowplot)
library(graphics)
library(lmtest)
library(stats)
library(MASS)
library(fpp2)
library(hrbrthemes)

## Warning: package 'hrbrthemes' was built under R version 4.1.3

## 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 CO.xlsx")

## New names:
## * `` -> `...4`
## * `` -> `...5`
## * `` -> `...6`
## * `` -> `...7`
## * `` -> `...8`
## * `` -> `...9`
## * `` -> `...10`
## * `` -> `...11`
## * `` -> `...12`
## * `` -> `...13`

mpdw <- mpdw[,2]
mpdw$`Data CO`<- as.numeric(mpdw$`Data CO`)
ts.mpdw <- ts(mpdw)
kableExtra::kable(head(mpdw) ,caption = 'Subset Data CO DKI Jakarta 2021-2022')

Subset Data CO DKI Jakarta 2021-2022
Data CO
0.38
0.41
0.30
0.35
0.76
0.40

view(mpdw)

# Train-test Split (80/20)
train <- mpdw[1:363,]
test <- mpdw[364:454,]

# Time-Series Object
data.ts <- ts(mpdw, start = 2020, frequency = 91)
train.ts <- ts(train)
test.ts <- ts(test)

Eksplorasi

Kestasioneran Data

plot(data.ts, col = "blue", type = "o", lwd = 1.5, main = "Time Series Polusi Udara DKI Jakarta 2020-2022",
     cex = 0.8)

#Fungsi Mencari Parameter Optimum
pemulusan <- function(x,min,max,method=c("SMA","DMA"),
                      metrik=c("SSE","MSE","RMSE","MAD","MAPE")){
  df.master <- data.frame()
  if(method=="SMA"){
    for(i in min:max) {
      sma <- TTR::SMA(x,i)
      ramal <- c(NA,sma)
      df <- cbind(aktual=c(x,NA),mulus=c(sma,NA),ramal)
      error <- df[,1]-df[,3]
      sse <- sum(error^2,na.rm=T)
      mse <- mean(error^2,na.rm=T)
      rmse <- sqrt(mse)
      mad <- mean(abs(error),na.rm=T)
      rerror <- error/df[,1]*100
      mape <- mean(abs(rerror),na.rm=T)
      ak <- data.frame(n=i,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
      df.master <- rbind(df.master,ak)
    }
  }
  else if(method=="DMA") {
    for(i in min:max) {
      df_ts_sma <- TTR::SMA(x,  n=i)
      df_ts_dma <- TTR::SMA(df_ts_sma, n=i)
      At <- 2*df_ts_sma-df_ts_dma
      Bt <- df_ts_sma-df_ts_dma
      pemulusan_dma <- At+Bt
      ramal_dma <- c(NA, pemulusan_dma)
      df_dma <- cbind(df_aktual=c(x,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
      error.dma <- df_dma[, 1] - df_dma[, 3]
      SSE.dma <- sum(error.dma^2, na.rm = T)
      MSE.dma <- mean(error.dma^2, na.rm = T)
      RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))
      MAD.dma <- mean(abs(error.dma), na.rm = T)
      r.error.dma <- (error.dma/df_dma[, 1])*100 
      MAPE.dma <- mean(abs(r.error.dma), na.rm = T)
      ak <- data.frame(n=i,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
                       MAD=MAD.dma,MAPE=MAPE.dma)
      df.master <- rbind(df.master,ak)
    }
  }
  opt <- df.master$n[which.min(df.master[,metrik])]
  return(list(Akurasi=df.master,n.optimum=paste("n optimum adalah",opt)))
}

SMA

pemulusan(test.ts,min=2,max=10,method = "SMA",metrik = "MSE")

## $Akurasi
##    n       SSE       MSE      RMSE       MAD     MAPE
## 1  2 11.649900 0.1308978 0.3617979 0.2637079 26.97370
## 2  3 11.112856 0.1262824 0.3553624 0.2691288 26.90669
## 3  4 10.301725 0.1184106 0.3441085 0.2599425 25.43789
## 4  5 10.179308 0.1183640 0.3440408 0.2556977 25.06402
## 5  6 10.165919 0.1195991 0.3458310 0.2583725 25.53401
## 6  7 10.351978 0.1232378 0.3510525 0.2628741 26.13185
## 7  8  9.929719 0.1196352 0.3458832 0.2585843 26.01712
## 8  9  9.737502 0.1187500 0.3446013 0.2595528 26.21553
## 9 10  9.741168 0.1202613 0.3467872 0.2624444 26.59527
## 
## $n.optimum
## [1] "n optimum adalah 5"

pemulusan(train.ts,min=2,max=10,method = "SMA",metrik = "MSE")

## $Akurasi
##    n      SSE        MSE      RMSE       MAD     MAPE
## 1  2 31.91530 0.08840803 0.2973349 0.2343213 24.58340
## 2  3 32.81743 0.09115954 0.3019264 0.2413981 25.33014
## 3  4 32.21861 0.08974542 0.2995754 0.2412326 25.16146
## 4  5 33.00971 0.09220589 0.3036542 0.2421285 25.09973
## 5  6 33.27291 0.09320142 0.3052891 0.2432960 25.22390
## 6  7 32.62904 0.09165461 0.3027451 0.2406942 24.94686
## 7  8 31.39455 0.08843537 0.2973808 0.2379683 24.74533
## 8  9 31.18505 0.08809337 0.2968053 0.2376962 24.80928
## 9 10 31.32687 0.08874467 0.2979004 0.2393144 25.10488
## 
## $n.optimum
## [1] "n optimum adalah 9"

Terkecil n = 5 DAN n = 2

DMA

pemulusan(test.ts,min=2,max=10,method = "DMA",metrik = "MSE")

## $Akurasi
##    n      SSE       MSE      RMSE       MAD     MAPE
## 1  2 13.06147 0.1484259 0.3852608 0.2897159 28.40301
## 2  3 17.85584 0.2076260 0.4556600 0.3395220 33.59611
## 3  4 16.22562 0.1931621 0.4395022 0.3342113 32.77488
## 4  5 15.76191 0.1922184 0.4384272 0.3360927 33.96116
## 5  6 16.11703 0.2014629 0.4488462 0.3493403 34.73170
## 6  7 16.57016 0.2124380 0.4609100 0.3576949 35.59723
## 7  8 15.71507 0.2067773 0.4547277 0.3491982 34.74398
## 8  9 14.97863 0.2024139 0.4499043 0.3532616 35.67812
## 9 10 14.69351 0.2040766 0.4517483 0.3569444 36.75386
## 
## $n.optimum
## [1] "n optimum adalah 2"

pemulusan(train.ts,min=2,max=10,method = "DMA",metrik = "MSE")

## $Akurasi
##    n      SSE       MSE      RMSE       MAD     MAPE
## 1  2 43.87227 0.1218674 0.3490952 0.2783472 29.22435
## 2  3 51.99775 0.1452451 0.3811104 0.2991589 31.28305
## 3  4 51.64412 0.1450677 0.3808776 0.3036833 31.47498
## 4  5 55.08291 0.1556014 0.3944635 0.3126576 32.37333
## 5  6 54.20414 0.1539890 0.3924144 0.3139804 32.18633
## 6  7 50.94212 0.1455489 0.3815087 0.3052950 31.07811
## 7  8 47.25780 0.1357983 0.3685082 0.2964592 30.15969
## 8  9 45.61789 0.1318436 0.3631028 0.2857928 29.21157
## 9 10 44.64222 0.1297739 0.3602415 0.2836250 29.32181
## 
## $n.optimum
## [1] "n optimum adalah 2"

Parameter optimum pada data training yang didapat untuk double moving average (DMA) berdasarkan nilai MSE terkecil adalah n = 2.

Berdasarkan perbandingan kedua metode, parameter optimum didapatkan nilai MSE dari DMA lebih kecil daripada nilai MSE SMA, sehingga metode pemulusan moving average akan dilanjutkan menggunakan metode DMA.

Untuk mencari parameter optimum single exponential smoothing (SES) dan double exponential smoothing (DES) dibuat sebuah fungsi optimasi dengan memanfaatkan iterasi dan evaluasi metrik tiap-tiap parameter (alpha dan beta).

#Fungsi Mencari Parameter Optimum
pemulusan2 <- function(x,alfa=NULL,beta=NULL,method=c("SES","DES"),
                       metrik=c("SSE","MSE","MAPE")){
  df.master <- data.frame()
  if(method=="SES"){
    for(i in alfa){
      df_ses <- HoltWinters(x, alpha = i, beta=F, gamma=F)
      sse <- df_ses$SSE 
      mse <- sse/length(x)
      rmse <-sqrt(mse)
      ak <- data.frame(Alpha=i,SSE=sse,MSE=mse,RMSE=rmse)
      df.master <- rbind(df.master,ak)
      datases <- data.frame(x, c(NA, df_ses$fitted[,1]))
      colnames(datases) = c("y","yhat")
    }
  }
  else if(method=="DES"){
    for (i in alfa) {
      for (j in beta) {
        df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse)
        df.master <- rbind(df.master,ak)
        datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
        colnames(datades) = c("y","yhat")
      }
    }
  }
  if(method=="SES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt)))
  }
  else if(method=="DES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    opt2 <- df.master$Beta[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt,
                                                  ",Beta optimum adalah",opt2)))
  }
}

SES

pemulusan2(test.ts,alfa=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
           method="SES",metrik="MSE")

## $Akurasi
##   Alpha       SSE       MSE      RMSE
## 1   0.1 12.829111 0.1409792 0.3754720
## 2   0.2 11.367727 0.1249201 0.3534403
## 3   0.3 10.742184 0.1180460 0.3435782
## 4   0.4 10.351219 0.1137497 0.3372679
## 5   0.5 10.078352 0.1107511 0.3327929
## 6   0.6  9.882760 0.1086018 0.3295478
## 7   0.7  9.748028 0.1071212 0.3272937
## 8   0.8  9.669146 0.1062544 0.3259668
## 9   0.9  9.649307 0.1060363 0.3256322
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.9"

pemulusan2(train.ts,alfa=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
           method="SES",metrik="MSE")

## $Akurasi
##   Alpha      SSE        MSE      RMSE
## 1   0.1 32.61816 0.08985719 0.2997619
## 2   0.2 29.57856 0.08148365 0.2854534
## 3   0.3 28.58694 0.07875190 0.2806277
## 4   0.4 28.23120 0.07777189 0.2788761
## 5   0.5 28.20486 0.07769934 0.2787460
## 6   0.6 28.43188 0.07832474 0.2798656
## 7   0.7 28.90396 0.07962524 0.2821794
## 8   0.8 29.64031 0.08165375 0.2857512
## 9   0.9 30.68155 0.08452217 0.2907270
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.5"

DMA

pemulusan2(test.ts,alfa=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
           beta=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),method="DES",metrik="MSE")

## $Akurasi
##    Alpha Beta      SSE       MSE      RMSE
## 1    0.1  0.1 20.83535 0.2289599 0.4784975
## 2    0.1  0.2 19.19865 0.2109742 0.4593193
## 3    0.1  0.3 18.95263 0.2082707 0.4563668
## 4    0.1  0.4 19.73210 0.2168362 0.4656568
## 5    0.1  0.5 21.57287 0.2370645 0.4868927
## 6    0.1  0.6 23.18390 0.2547681 0.5047456
## 7    0.1  0.7 23.67742 0.2601914 0.5100896
## 8    0.1  0.8 23.42582 0.2574266 0.5073722
## 9    0.1  0.9 23.23540 0.2553340 0.5053059
## 10   0.2  0.1 14.45127 0.1588052 0.3985037
## 11   0.2  0.2 14.89695 0.1637028 0.4046020
## 12   0.2  0.3 15.83202 0.1739783 0.4171070
## 13   0.2  0.4 16.58054 0.1822038 0.4268533
## 14   0.2  0.5 17.21588 0.1891855 0.4349546
## 15   0.2  0.6 17.95736 0.1973336 0.4442225
## 16   0.2  0.7 18.84234 0.2070587 0.4550370
## 17   0.2  0.8 19.81244 0.2177191 0.4666038
## 18   0.2  0.9 20.80615 0.2286390 0.4781621
## 19   0.3  0.1 12.78734 0.1405202 0.3748603
## 20   0.3  0.2 13.52042 0.1485760 0.3854556
## 21   0.3  0.3 14.37507 0.1579678 0.3974516
## 22   0.3  0.4 15.25421 0.1676287 0.4094249
## 23   0.3  0.5 16.20567 0.1780843 0.4220004
## 24   0.3  0.6 17.20921 0.1891122 0.4348703
## 25   0.3  0.7 18.23381 0.2003715 0.4476287
## 26   0.3  0.8 19.25846 0.2116315 0.4600342
## 27   0.3  0.9 20.25796 0.2226149 0.4718209
## 28   0.4  0.1 11.98693 0.1317245 0.3629387
## 29   0.4  0.2 12.80041 0.1406639 0.3750519
## 30   0.4  0.3 13.70517 0.1506062 0.3880802
## 31   0.4  0.4 14.64938 0.1609822 0.4012258
## 32   0.4  0.5 15.61287 0.1715700 0.4142101
## 33   0.4  0.6 16.56761 0.1820617 0.4266868
## 34   0.4  0.7 17.49409 0.1922427 0.4384549
## 35   0.4  0.8 18.38328 0.2020141 0.4494598
## 36   0.4  0.9 19.23408 0.2113635 0.4597429
## 37   0.5  0.1 11.50103 0.1263849 0.3555066
## 38   0.5  0.2 12.34556 0.1356655 0.3683280
## 39   0.5  0.3 13.25083 0.1456135 0.3815934
## 40   0.5  0.4 14.16702 0.1556815 0.3945650
## 41   0.5  0.5 15.07008 0.1656052 0.4069462
## 42   0.5  0.6 15.94315 0.1751995 0.4185684
## 43   0.5  0.7 16.77696 0.1843622 0.4293742
## 44   0.5  0.8 17.56527 0.1930250 0.4393461
## 45   0.5  0.9 18.30373 0.2011399 0.4484862
## 46   0.6  0.1 11.17218 0.1227712 0.3503873
## 47   0.6  0.2 12.02119 0.1321010 0.3634571
## 48   0.6  0.3 12.90441 0.1418067 0.3765723
## 49   0.6  0.4 13.78056 0.1514347 0.3891462
## 50   0.6  0.5 14.63236 0.1607951 0.4009927
## 51   0.6  0.6 15.45058 0.1697866 0.4120517
## 52   0.6  0.7 16.23280 0.1783824 0.4223534
## 53   0.6  0.8 16.98383 0.1866354 0.4320132
## 54   0.6  0.9 17.71651 0.1946870 0.4412335
## 55   0.7  0.1 10.94534 0.1202784 0.3468118
## 56   0.7  0.2 11.79089 0.1295702 0.3599586
## 57   0.7  0.3 12.65397 0.1390546 0.3729003
## 58   0.7  0.4 13.50333 0.1483882 0.3852119
## 59   0.7  0.5 14.32890 0.1574604 0.3968128
## 60   0.7  0.6 15.12910 0.1662538 0.4077423
## 61   0.7  0.7 15.90944 0.1748290 0.4181256
## 62   0.7  0.8 16.68016 0.1832984 0.4281337
## 63   0.7  0.9 17.45162 0.1917760 0.4379224
## 64   0.8  0.1 10.80148 0.1186976 0.3445252
## 65   0.8  0.2 11.64678 0.1279866 0.3577521
## 66   0.8  0.3 12.50149 0.1373790 0.3706467
## 67   0.8  0.4 13.34330 0.1466296 0.3829225
## 68   0.8  0.5 14.16824 0.1556949 0.3945819
## 69   0.8  0.6 14.97949 0.1646098 0.4057213
## 70   0.8  0.7 15.78417 0.1734524 0.4164762
## 71   0.8  0.8 16.58957 0.1823030 0.4269696
## 72   0.8  0.9 17.40039 0.1912131 0.4372792
## 73   0.9  0.1 10.73963 0.1180179 0.3435373
## 74   0.9  0.2 11.59603 0.1274289 0.3569718
## 75   0.9  0.3 12.46123 0.1369366 0.3700494
## 76   0.9  0.4 13.32075 0.1463818 0.3825988
## 77   0.9  0.5 14.17581 0.1557782 0.3946874
## 78   0.9  0.6 15.03327 0.1652007 0.4064489
## 79   0.9  0.7 15.90252 0.1747529 0.4180346
## 80   0.9  0.8 16.79386 0.1845479 0.4295904
## 81   0.9  0.9 17.71903 0.1947146 0.4412648
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.9 ,Beta optimum adalah 0.1"

pemulusan2(train.ts,alfa=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
           beta=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),method="DES",metrik="MSE")

## $Akurasi
##    Alpha Beta      SSE        MSE      RMSE
## 1    0.1  0.1 36.40236 0.10028198 0.3166733
## 2    0.1  0.2 37.55930 0.10346915 0.3216662
## 3    0.1  0.3 38.85652 0.10704276 0.3271739
## 4    0.1  0.4 40.39196 0.11127262 0.3335755
## 5    0.1  0.5 41.94999 0.11556470 0.3399481
## 6    0.1  0.6 43.94310 0.12105538 0.3479301
## 7    0.1  0.7 46.26070 0.12743993 0.3569873
## 8    0.1  0.8 49.11181 0.13529424 0.3678237
## 9    0.1  0.9 52.31800 0.14412671 0.3796402
## 10   0.2  0.1 32.01798 0.08820381 0.2969913
## 11   0.2  0.2 33.93299 0.09347932 0.3057439
## 12   0.2  0.3 36.11635 0.09949408 0.3154268
## 13   0.2  0.4 38.61296 0.10637178 0.3261469
## 14   0.2  0.5 41.25954 0.11366265 0.3371389
## 15   0.2  0.6 43.91831 0.12098709 0.3478320
## 16   0.2  0.7 46.59479 0.12836031 0.3582741
## 17   0.2  0.8 49.30105 0.13581558 0.3685316
## 18   0.2  0.9 51.93202 0.14306343 0.3782373
## 19   0.3  0.1 30.96611 0.08530609 0.2920721
## 20   0.3  0.2 33.23710 0.09156225 0.3025925
## 21   0.3  0.3 35.70322 0.09835599 0.3136176
## 22   0.3  0.4 38.26222 0.10540555 0.3246622
## 23   0.3  0.5 40.82679 0.11247049 0.3353662
## 24   0.3  0.6 43.34846 0.11941723 0.3455680
## 25   0.3  0.7 45.81578 0.12621428 0.3552665
## 26   0.3  0.8 48.30703 0.13307721 0.3647975
## 27   0.3  0.9 50.94911 0.14035567 0.3746407
## 28   0.4  0.1 30.63575 0.08439602 0.2905099
## 29   0.4  0.2 33.02667 0.09098256 0.3016332
## 30   0.4  0.3 35.51646 0.09784148 0.3127962
## 31   0.4  0.4 38.02892 0.10476285 0.3236709
## 32   0.4  0.5 40.53124 0.11165630 0.3341501
## 33   0.4  0.6 43.03670 0.11855841 0.3443231
## 34   0.4  0.7 45.56386 0.12552027 0.3542884
## 35   0.4  0.8 48.07794 0.13244612 0.3639315
## 36   0.4  0.9 50.48389 0.13907407 0.3729264
## 37   0.5  0.1 30.63504 0.08439405 0.2905065
## 38   0.5  0.2 33.06703 0.09109373 0.3018174
## 39   0.5  0.3 35.54402 0.09791741 0.3129176
## 40   0.5  0.4 38.01428 0.10472254 0.3236086
## 41   0.5  0.5 40.45493 0.11144608 0.3338354
## 42   0.5  0.6 42.84321 0.11802537 0.3435482
## 43   0.5  0.7 45.13590 0.12434132 0.3526207
## 44   0.5  0.8 47.28982 0.13027498 0.3609363
## 45   0.5  0.9 49.29542 0.13580005 0.3685106
## 46   0.6  0.1 30.89944 0.08512243 0.2917575
## 47   0.6  0.2 33.36709 0.09192036 0.3031837
## 48   0.6  0.3 35.84894 0.09875741 0.3142569
## 49   0.6  0.4 38.30334 0.10551884 0.3248366
## 50   0.6  0.5 40.70552 0.11213641 0.3348677
## 51   0.6  0.6 43.03671 0.11855842 0.3443231
## 52   0.6  0.7 45.29262 0.12477306 0.3532323
## 53   0.6  0.8 47.49591 0.13084272 0.3617219
## 54   0.6  0.9 49.69237 0.13689358 0.3699913
## 55   0.7  0.1 31.44295 0.08661970 0.2943123
## 56   0.7  0.2 33.98600 0.09362535 0.3059826
## 57   0.7  0.3 36.53866 0.10065746 0.3172656
## 58   0.7  0.4 39.07464 0.10764362 0.3280909
## 59   0.7  0.5 41.58836 0.11456849 0.3384797
## 60   0.7  0.6 44.09268 0.12146744 0.3485218
## 61   0.7  0.7 46.62041 0.12843088 0.3583725
## 62   0.7  0.8 49.21901 0.13558956 0.3682249
## 63   0.7  0.9 51.93944 0.14308386 0.3782643
## 64   0.8  0.1 32.30642 0.08899840 0.2983260
## 65   0.8  0.2 34.99786 0.09641283 0.3105042
## 66   0.8  0.3 37.72539 0.10392670 0.3223766
## 67   0.8  0.4 40.48456 0.11152771 0.3339577
## 68   0.8  0.5 43.29422 0.11926782 0.3453517
## 69   0.8  0.6 46.19110 0.12724820 0.3567187
## 70   0.8  0.7 49.22306 0.13560073 0.3682400
## 71   0.8  0.8 52.43984 0.14446237 0.3800821
## 72   0.8  0.9 55.88523 0.15395380 0.3923695
## 73   0.9  0.1 33.55291 0.09243227 0.3040268
## 74   0.9  0.2 36.49502 0.10053725 0.3170761
## 75   0.9  0.3 39.53262 0.10890529 0.3300080
## 76   0.9  0.4 42.68552 0.11759097 0.3429154
## 77   0.9  0.5 45.99566 0.12670980 0.3559632
## 78   0.9  0.6 49.51890 0.13641571 0.3693450
## 79   0.9  0.7 53.31782 0.14688105 0.3832506
## 80   0.9  0.8 57.45758 0.15828533 0.3978509
## 81   0.9  0.9 62.00782 0.17082044 0.4133043
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.5 ,Beta optimum adalah 0.1"

Metode SES lebih baik

Forecasting

Double Moving Average

#Single Moving  Average dengan n=2
df_ts_sma <- TTR::SMA(train.ts,  n=2)
#Double Moving  Average dengan n=2
df_ts_dma <- TTR::SMA(df_ts_sma, n=2)
At <- 2*df_ts_sma-df_ts_dma
Bt <- df_ts_sma-df_ts_dma
pemulusan_dma <- At+Bt
ramal_dma <- c(NA, pemulusan_dma)
df_dma <- cbind(df_aktual=c(train.ts,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
head(df_dma)

##      df_aktual pemulusan_dma ramal_dma
## [1,]      0.38            NA        NA
## [2,]      0.41            NA        NA
## [3,]      0.30         0.315        NA
## [4,]      0.35         0.295     0.315
## [5,]      0.76         0.785     0.295
## [6,]      0.40         0.605     0.785

Eksploratif: DMA Plot

par(bg = '#141415')
ts.plot(mpdw$`Data CO`, xlab="periode waktu", ylab="Yt", col="blue", lty=1,
      gpars = list(col.main="white",col.axis="white",col.sub="white"))
points(mpdw$`Data CO`,col="white")
par(col = "white")
lines(pemulusan_dma,col="red",lwd=2,lty=1)
par(col = "white")
lines(ramal_dma,col="green",lwd= 2,lty=1)
title("Rataan Bergerak  Berganda N=2", cex.main=1, font.main=4 ,col.main="white")
legend("topright", c("Data Aktual","Pemulusan","Ramalan"),lty=1,col=c ("blue","red","green"))
box(col="white",lwd=2)

n = 2 lebih mendekati sebaran pola dari data aktual

Single Exponential Smoothing

#SES alpha=0.5
df_des <- HoltWinters(train.ts, alpha = 0.5,gamma=F)
datades <- data.frame(train.ts, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("y","yhat")
head(datades)

##      y      yhat
## 1 0.38        NA
## 2 0.41        NA
## 3 0.30 0.4400000
## 4 0.35 0.3994217
## 5 0.76 0.4039284
## 6 0.40 0.6126526

ramal_des <- forecast::forecast(df_des,h=3)
(df_ramal_des <- data.frame(ramal_des))

##     Point.Forecast     Lo.80    Hi.80     Lo.95    Hi.95
## 364       1.055251 0.6950855 1.415416 0.5044255 1.606076
## 365       1.057279 0.6539340 1.460623 0.4404163 1.674141
## 366       1.059306 0.6163650 1.502248 0.3818860 1.736727

Eksploratif: DES Plot

Plot

par(bg = '#141415')
ts.plot(mpdw$`Data CO`, xlab="periode  waktu", ylab="Yt", col="blue", lty=1)
points(mpdw$`Data CO`,col="white")
par(col="white")
lines(datades$yhat, col="red",lwd=2) #nilai dugaan
par(col="white")
lines(ts(df_ramal_des$Point.Forecast,start = 101), col="green",lwd=2) #nilai dugaan
title("Single Exponential Smoothing Alpha=0.5 ",cex.main=1, font.main=4 ,col.main="white")
legend("bottomright", c("Data aktual", "Fitted SES","Ramalan"), lty=1,col=c ("blue","red","green"))
box(col="white",lwd=2)

# ARIMA ## CEK TAILS OFF/CUT OFF

ggAcf(ts(mpdw$`Data CO`),col="white",lwd=2)+labs( x="Lag",y = "ACF",
         title="Plot ACF Data CO DKI Jakarta ",
        subtitle = "(Polusi udara 2021-2022)")+
  theme_ft_rc()+
  theme(
    plot.title = element_text(size = 14L,
                              face = "bold",
                              hjust = 0.5),
    plot.subtitle = element_text(size = 11L,
                              face = "plain",
                              hjust = 0.5),plot.background=element_rect(fill="#141415",color="#141415"),panel.background=element_rect(fill="#141415",color="#141415"),
      axis.title = element_text(color="white"),axis.text = element_text(color="white")
  )

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

ggPacf(ts(mpdw$`Data CO`),col="white",lwd=2)+labs( x="Lag",y = "PACF",
         title="Plot PACF Kadar CO DKI Jakarta  ",
        subtitle = "(Juli 2012-2022)",col="white")+
  theme_ft_rc()+
  theme(
    plot.title = element_text(size = 14L,
                              face = "bold",
                              hjust = 0.5),
    plot.subtitle = element_text(size = 11L,
                              face = "plain",
                              hjust = 0.5),plot.background=element_rect(fill="#141415",color="#141415"),panel.background=element_rect(fill="#141415",color="#141415"),
      axis.title = element_text(color="white"),axis.text = element_text(color="white")
  )

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

Uji Statistik

adf.test(data.ts)

## Warning in adf.test(data.ts): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.ts
## Dickey-Fuller = -4.8343, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

SUdah Stasioner

EACF

eacf(data.ts)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x o  o  o  o 
## 1 o o x x x o o o o x o  o  o  o 
## 2 x o x o x o o o o o x  o  o  o 
## 3 x x x o o o o o o o x  o  o  o 
## 4 x x x x x o o o o o x  o  o  o 
## 5 x x x x x x o o o o x  o  o  o 
## 6 x x x o o x o o o o x  o  o  o 
## 7 x x x x o x o o o o o  o  o  o

ARIMA (1,0,5) ARIMA (1,0,6) ARIMA (0,0,10) ARIMA (1,0,10)

Perbandingan Kebaikan Model Tentatif

model1 <- Arima(data.ts,order = c(1,0,5),method ="ML")
model2 <- Arima(data.ts,order = c(1,0,6),method ="ML")
model3 <- Arima(data.ts,order = c(0,0,10),method ="ML")
Model <- c("ARIMA (1,0,5)","ARIMA (1,0,6)","ARIMA (0,0,10)")
AIC <- c(model1$aic,model2$aic,model3$aic)
BIC <- c(model1$bic,model2$bic,model3$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)

Model	AIC	BIC
ARIMA (1,0,5)	105.2526	138.1974
ARIMA (1,0,6)	107.1280	144.1908
ARIMA (0,0,10)	105.6776	155.0947

paste("Model yang terbaik adalah model",Akurasi$Model[which.min(Akurasi[,"AIC"])])

## [1] "Model yang terbaik adalah model ARIMA (1,0,5)"

coeftest(model1)

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.922349   0.052995 17.4045 < 2.2e-16 ***
## ma1       -0.337115   0.071384 -4.7225 2.329e-06 ***
## ma2       -0.248415   0.057858 -4.2935 1.759e-05 ***
## ma3       -0.153593   0.052832 -2.9072  0.003647 ** 
## ma4        0.108716   0.050292  2.1617  0.030642 *  
## ma5       -0.093699   0.054761 -1.7110  0.087073 .  
## intercept  1.036335   0.043863 23.6266 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Signifikan semua

Uji Diagnostik

Diagnostik Model: Eksploratif

sisaan <- model1$residuals
par(mfrow=c(2,2))
par(bg = '#141415')
qqnorm(sisaan)
box(col="white",lwd=2)
qqline(sisaan, col = "red", lwd =1, col.lab="white",
       col.axis="white",col.sub = "white")
box(col="white",lwd=2)
plot(c(1:length(sisaan)),sisaan,col="white",col.lab="white",col.axis="white")
box(col="white",lwd=2)
acf(sisaan,col="white",col.sub = "white",col.axis="white", col.lab="white")
box(col="white",lwd=2)
pacf(sisaan,col="white",col.sub = "white",col.axis="white", col.lab="white",col.main="white")
box(col="white",lwd=2)

Berdasarkan hasil eksplorasi menggunakan Q-Q plot, terlihat bahwa sisaan berdistribusi mengikuti garis normal, sehingga dapat dikatakan bahwa sisaan menyebar normal. Kemudian, plot sisaan yang diperoleh menunjukkan bahwa sisaan memiliki pola acak dan tersebar di sekitar nilai nol. Sedangkan pada plot ACF dan PACF, nilai awal amatan tidak melewati garis signifikan, atau dapat dikatakan bahwa sisaan saling bebas.

Diagnostik Model: Uji Formal

1. Sisaan Menyebar Normal

Uji formal ini dilakukan dengan Jarque Bera test dan Shapiro-Wilk test.

Hipotesis yang diuji: H0 : Sisaan Menyebar Normal H1: Sisaan Tidak Menyebar Normal

jarque.bera.test(sisaan)

## 
##  Jarque Bera Test
## 
## data:  sisaan
## X-squared = 12.305, df = 2, p-value = 0.002128

shapiro.test(sisaan)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.98896, p-value = 0.001686

2. Sisaan Saling Bebas

Uji formal ini dilakukan dengan LJung-Box test.

Hipotesis yang diuji: H0 : Sisaan antara lag saling bebas H1: Sisaan antara lag tidak saling bebas

Box.test(sisaan, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 0.013066, df = 1, p-value = 0.909

3. Nilai Tengah Sisaan Sama dengan Nol

Uji formal ini dilakukan dengan t-test.

Hipotesis yang diuji: H0 : Nilai tengah sisaan sama dengan nol H1: Nilai tengah sisaan tidak sama dengan nol

t.test(sisaan, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.23454, df = 453, p-value = 0.8147
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.02169267  0.02757237
## sample estimates:
##   mean of x 
## 0.002939849

FORECASTING

ramalan <- forecast::forecast(Arima(train.ts, order=c(1,0,5),method="ML",include.drift = TRUE), h=8) 
data.ramalan <- ramalan$mean
par(bg = '#141415')
plot(ramalan,col="white",col.sub ="white",col.axis="white",
     col.lab="white",col.main="white",lwd=2)
box(col="white",lwd=2)

par(bg = '#141415')
ts.plot(mpdw$`Data CO`,col="white",lwd=2,main="Perbandingan Peramalan dan Aktual",gpars = list(col.main="white",col.axis="white",col.sub="white"))
lines(data.ramalan, col = "red",lwd=2)
box(col="white",lwd=2)

perbandingan.temp<-matrix(data=c(test.ts[1:5], data.ramalan[1:5]), nrow = 5, ncol = 2)
colnames(perbandingan.temp)<-c("Aktual","Hasil Forecast")
head(perbandingan.temp)

##      Aktual Hasil Forecast
## [1,]   1.80      0.9611172
## [2,]   1.91      1.0259902
## [3,]   0.65      1.0004032
## [4,]   0.63      0.9623544
## [5,]   0.90      0.9736433