Peramalan Jumlah Kunjungan Wisatawan Mancanegara ke Indonesia (Januari 2013 - Agustus 2022) dengan Metode Pemulusan dan Pemodelan ARIMA

P1 - Kelompok 7 STA1341 2022

Anggota Kelompok:

Aprilia Permata Putri (G1401201002)
Nana Oktaviana (G1401201006)
Ferista Wahyu Saputri (G1401201008)
Ainaini Salsabila (G1401201055)
Akmal Riza Wibisono (G1401201086)

Latar Belakang

Pada tahun 2020, Kementerian Pariwisata dan Ekonomi Kreatif menetapkan sektor pariwisata sebagai leading sector penyumbang devisa Indonesia setelah sektor industri kelapa sawit. Penetapan pariwisata sebagai leading sector dilakukan usai sektor ini berhasil mencatatkan capaian mentereng pada periode 2013-2019 dengan jumlah kunjungan wisatawan mancanegara yang terus meningkat. Peningkatan jumlah kunjungan wisatawan selayaknya diiringi dengan partisipasi pemerintah dalam pembangunan sektor pariwisata berupa sarana dan prasarana publik, infrastruktur, transportasi, dan layanan lainnya agar terjadinya pembangunan ekonomi yang efektif dan efisien (Herawati 2016).

Di tengah meningkatnya semangat pembangunan pariwisata nasional, berita buruk datang ketika pandemi virus korona melanda dunia sejak tahun 2020. Kebijakan lockdown berimbas besar pada sektor-sektor utama nasional, tak terkecuali sektor pariwisata. Kini, pada akhir tahun 2022, pandemi mulai terkendali diikuti dengan memulihnya kegiatan ekonomi dunia. Dengan predikat leading sector penyumbang devisa negara, tentunya pemerintah dan industri pariwisata membutuhkan peramalan kunjungan wisatawan mancanegara yang akurat agar mampu menyusun dan mengambil kebijakan yang tepat sehingga peningkatan geliat ekonomi kreatif yang sempat terhambat saat pandemi dapat kembali dalam waktu yang sesegera mungkin.

Tinjauan Pustaka

Data Deret Waktu (Time Series)

Data time series (dalam bahasa Indonesia disebut sebagai data deret waktu) adalah kumpulan amatan yang terurut dalam waktu (Hanke dan Wichern 2005). Proses analisis melalui metode time series dapat dimaknai sebagai suatu proses analisis pola hubungan antara suatu peubah tertentu dengan peubah waktu. Tujuan akhir dalam proses analisis data time series adalah menghasilkan ramalan nilai amatan di masa depan yang optimum dengan galat minimum (Wei 2006).

Metode Pemulusan

Double Moving Average

DMA merupakan metode yang menghitung rata-rata bergerak yang kedua atau menghitung rata-rata bergerak dari rata-rata bergerak atau SMA. DMA menghilangkan lebih banyak keacakan daripada rata-rata bergerak sederhana yang sama panjangnya (DMA dapat mengakumulasi dan mengakomodir trend yang ada) (Makridakis et al. 1997).

Holt Winters

Metode Holt-Winters diperkenalkan oleh Holt (1957) dan Winters (1960) merupakan perkembangan dari metode exponential smoothing yang menggunakan tiga parameter pemulusan, yaitu α parameter untuk pemulusan level, β parameter pemulusan kecenderungan (trend), dan γ parameter pemulusan musiman (seasonal) (Makridakis et al. 1997).

Analisis Data

Load Package

library(googlesheets4)
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(FinTS)

Import Data

Data historis yang digunakan merupakan data sekunder sejumlah 116 amatan yang bersumber dari situs resmi Badan Pusat Statistik (bps.go.id) berupa total kunjungan bulanan wisatawan mancanegara ke Indonesia dari periode Januari 2013 hingga Agustus 2022.

gs4_deauth()
wisata <- read_sheet("https://docs.google.com/spreadsheets/d/1oWrxKrK3FVRTbVCGb-ijLp0OY8_WH5dChhS4CvZnk6k/edit?usp=sharing")
wisata <- wisata[,-c(1,4)]
colnames(wisata) <- c("Bulan","Wisatawan")
wisata$Bulan <- as.Date(wisata$Bulan,"%y/%m/%d")
kableExtra::kable(head(wisata) ,caption = 'Subset Data Jumlah Kunjungan Wisatawan Mancanegara ke Indonesia 2013-2022',  align = rep('c',7))

Subset Data Jumlah Kunjungan Wisatawan Mancanegara ke Indonesia 2013-2022
Bulan	Wisatawan
2013-01-31	614328
2013-02-28	678415
2013-03-31	725316
2013-04-30	646117
2013-05-31	700708
2013-06-30	789594

Plot Time Series

ggplot(wisata, aes(x=Bulan, y=Wisatawan)) +
  geom_line(lwd=0.7) + 
  scale_x_date(date_breaks = "12 month", date_labels = "%b %Y") +
  labs(title = "Plot Deret Waktu Data Jumlah Kunjungan Wisatawan Mancanegara Bulanan",
       subtitle = "Periode Januari 2013 - Agustus 2022") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA,color="black"))

Berdasarkan plot deret waktu, dapat diamati bahwa total kunjungan bulanan wisatawan mancanegara ke Indonesia terus mengalami peningkatan sejak tahun 2013 hingga awal tahun 2020. Hal ini menjadi salah satu indikator keberhasilan pemerintah Indonesia, khususnya Kementerian Pariwisata, dalam mempromosikan destinasi wisata nasional di lingkup internasional.

Terhitung sejak pandemi COVID-19 melanda dunia sejak Maret 2020, total kunjungan wisatawan mancanegara ke Indonesia ikut menurun drastis 88% hingga ke titik 150.000 kunjungan wisatawan per bulan. Sejak April 2020 hingga saat ini, total kunjungan wisatawan masih stabil dengan rata-rata 150.000 wisatawan per bulannya meskipun sudah mulai menunjukkan tren peningkatan di pertengahan tahun 2022 seiring dengan mulai meredanya penyebaran virus korona di tengah masyarakat internasional.

Splitting Data

Sebelum dilakukan pemulusan, data aktual dibagi menjadi data training dan data testing. Dengan memperhatikan pola data, ditetapkan data training merupakan data pertama sampai data ke-102, sedangkan data testing adalah data ke-103 sampai data ke-116. Kemudian, data diubah menjadi format data time series.

train <- wisata[1:102,2]
test <- wisata[103:116,2]
wisata.ts = ts(wisata$Wisatawan, frequency = 12, start= 2013)
training.ts <- ts(train, start=2013, frequency = 12)
testing.ts <- ts(test, start=2021, frequency = 12)

plot(training.ts, main = "Kunjungan Wisatawan Mancanegara Periode 2013-2022", 
     xlab = "Tahun", ylab="Wisatawan", xlim=c(2013,2022), ylim=c(0,1600000), lwd=2)
points(training.ts)
points(testing.ts)
lines(training.ts, col ="black", lwd=2)
lines(testing.ts, col ="red", lwd=2)
legend("bottomleft", c("Data Training", "Data Testing"), lty=1, col=c("black", "red"))

Berdasarkan plot deret waktu data training dan testing jumlah kunjungan wisatawan mancanegara, terlihat bahwa data tidak stasioner karena plot data tidak menyebar di sekitar rata-rata serta ragam yang konstan. Jadi data tersebut tidak stasioner dalam rata-rata dan ragam.

Metode Smoothing

Double Moving Average

#Fungsi Mencari Parameter Optimum
n = 10
df_acc_dma <- data.frame()

for(p in 2:n){
  mulussma <- TTR::SMA(training.ts, n=p)
  ts_dma <- TTR::SMA(mulussma, n=p)
  At <- 2*mulussma - ts_dma
  Bt <- 2/(2-1)*mulussma - ts_dma
  mulusdma <- At+Bt
  ramal_dma <- c(NA, mulusdma)
  df_dma <- cbind(Aktual=c(training.ts,NA), pemulusan=c(mulusdma,NA),
                  ramal_dma)
  head(df_dma)
  tail(df_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)
  
  vektor <- c(p, SSE.dma,MSE.dma,MAPE.dma,RMSE.dma,MAD.dma)
  df_acc_dma <- rbind(df_acc_dma, vektor)
}
colnames(df_acc_dma) <- c("n", "SSE","MSE","MAPE", "RMSE","MAD")
kableExtra::kable(df_acc_dma, align = rep('c',7), caption = "Perbandingan Error Pemulusan Metode DMA")

Perbandingan Error Pemulusan Metode DMA
n	SSE	MSE	MAPE	RMSE	MAD
2	9.804104e+13	9.903135e+11	105.9268	995145	903283.3
3	9.827669e+13	1.013162e+12	109.2182	1006559	913753.5
4	9.877374e+13	1.039724e+12	115.4625	1019668	929014.2
5	9.929874e+13	1.067728e+12	122.4162	1033309	945516.4
6	9.994287e+13	1.098273e+12	129.0438	1047985	960642.6
7	1.007192e+14	1.131676e+12	136.6127	1063803	977078.9
8	1.012844e+14	1.164189e+12	144.9530	1078976	991761.0
9	1.018066e+14	1.197725e+12	157.2560	1094406	1010581.8
10	1.023814e+14	1.233511e+12	172.2050	1110635	1033334.1

Parameter optimum yang didapat untuk double moving average (DMA) berdasarkan metrik akurasi di atas adalah n = 2.

# Forecasting
 mulussma <- TTR::SMA(training.ts, n=2)
  ts_dma <- TTR::SMA(mulussma, n=2)
  At <- 2*mulussma - ts_dma
  Bt <- 2/(2-1)*(mulussma-ts_dma)
  mulusdma <- At+Bt
  ramal_dma <- c(NA, mulusdma) 
  df_dma <- cbind(Aktual=wisata.ts, pemulusan = c(mulusdma,rep(NA,14)), ramal_dma= c(ramal_dma, 
  rep(NA,13)))
  df_dma

##           Aktual  pemulusan  ramal_dma
## Jan 2013  614328         NA         NA
## Feb 2013  678415         NA         NA
## Mar 2013  725316  785106.50         NA
## Apr 2013  646117  661493.00  785106.50
## May 2013  700708  654956.50  661493.00
## Jun 2013  789594  852758.75  654956.50
## Jul 2013  717784  766496.00  852758.75
## Aug 2013  771009  730457.75  766496.00
## Sep 2013  770878  810764.00  730457.75
## Oct 2013  719903  707061.00  810764.00
## Nov 2013  807422  791070.50  707061.00
## Dec 2013  860655  939602.50  791070.50
## Jan 2014  753079  766109.75  939602.50
## Feb 2014  702666  609380.75  766109.75
## Mar 2014  765607  743532.50  609380.75
## Apr 2014  726332  763719.00  743532.50
## May 2014  752363  729414.50  763719.00
## Jun 2014  851475  895776.25  729414.50
## Jul 2014  777210  832977.75  895776.25
## Aug 2014  826821  783525.00  832977.75
## Sep 2014  791296  819623.00  783525.00
## Oct 2014  808767  786491.00  819623.00
## Nov 2014  764461  766487.75  786491.00
## Dec 2014  915334  919822.75  766487.75
## Jan 2015  724698  790193.75  919822.75
## Feb 2015  794302  668726.00  790193.75
## Mar 2015  792804  844632.50  668726.00
## Apr 2015  750999  739424.25  844632.50
## May 2015  794294  773764.00  739424.25
## Jun 2015  815307  853031.50  773764.00
## Jul 2015  815351  831121.75  853031.50
## Aug 2015  853244  862750.25  831121.75
## Sep 2015  870351  903047.50  862750.25
## Oct 2015  826196  827987.50  903047.50
## Nov 2015  777976  732804.75  827987.50
## Dec 2015  913828  911626.00  732804.75
## Jan 2016  814303  891310.75  911626.00
## Feb 2016  888309  832166.75  891310.75
## Mar 2016  915019  977201.00  832166.75
## Apr 2016  901095  917646.50  977201.00
## May 2016  915206  908290.75  917646.50
## Jun 2016  857651  853845.50  908290.75
## Jul 2016 1032741 1033347.25  853845.50
## Aug 2016 1031986 1163114.75 1033347.25
## Sep 2016 1006653  999753.50 1163114.75
## Oct 2016 1040651 1030150.75  999753.50
## Nov 2016 1002333 1018252.00 1030150.75
## Dec 2016 1113328 1112338.25 1018252.00
## Jan 2017 1107968 1189874.25 1112338.25
## Feb 2017 1023388  998223.00 1189874.25
## Mar 2017 1059777 1005439.25  998223.00
## Apr 2017 1171386 1226580.00 1005439.25
## May 2017 1148588 1226595.25 1226580.00
## Jun 2017 1144001 1125755.75 1226595.25
## Jul 2017 1370591 1423798.25 1125755.75
## Aug 2017 1393243 1568848.50 1423798.25
## Sep 2017 1250231 1231467.00 1568848.50
## Oct 2017 1161565 1032139.50 1231467.00
## Nov 2017 1062030  970646.75 1032139.50
## Dec 2017 1147031 1093630.00  970646.75
## Jan 2018 1097839 1149291.75 1093630.00
## Feb 2018 1197503 1185525.00 1149291.75
## Mar 2018 1363426 1479654.75 1185525.00
## Apr 2018 1302321 1411487.00 1479654.75
## May 2018 1242705 1181972.25 1411487.00
## Jun 2018 1322674 1297954.25 1181972.25
## Jul 2018 1547231 1663347.00 1297954.25
## Aug 2018 1511021 1670386.25 1663347.00
## Sep 2018 1370943 1308766.00 1670386.25
## Oct 2018 1291605 1166712.00 1308766.00
## Nov 2018 1157483 1064449.00 1166712.00
## Dec 2018 1405554 1366980.25 1064449.00
## Jan 2019 1201735 1336833.50 1366980.25
## Feb 2019 1243996 1101697.00 1336833.50
## Mar 2019 1311911 1360585.50 1101697.00
## Apr 2019 1274231 1315747.25 1360585.50
## May 2019 1249536 1215102.25 1315747.25
## Jun 2019 1434103 1461723.50 1215102.25
## Jul 2019 1468173 1615115.75 1461723.50
## Aug 2019 1530268 1571344.25 1615115.75
## Sep 2019 1388719 1399903.00 1571344.25
## Oct 2019 1346434 1229701.00 1399903.00
## Nov 2019 1280781 1232654.00 1229701.00
## Dec 2019 1377067 1351898.75 1232654.00
## Jan 2020 1290411 1340961.50 1351898.75
## Feb 2020  872765  703361.50 1340961.50
## Mar 2020  486155   76268.00  703361.50
## Apr 2020  158066 -213913.75   76268.00
## May 2020  161842  -83280.75 -213913.75
## Jun 2020  156561  158072.75  -83280.75
## Jul 2020  155742  151576.50  158072.75
## Aug 2020  161549  162386.50  151576.50
## Sep 2020  148984  150198.00  162386.50
## Oct 2020  152293  143696.50  150198.00
## Nov 2020  144476  145003.50  143696.50
## Dec 2020  164079  163117.00  145003.50
## Jan 2021  126515  131826.25  163117.00
## Feb 2021  105788   72433.25  131826.25
## Mar 2021  119979  107981.50   72433.25
## Apr 2021  112756  121593.50  107981.50
## May 2021  139433  140685.00  121593.50
## Jun 2021  126844  143704.50  140685.00
## Jul 2021  127249         NA  143704.50
## Aug 2021  118533         NA         NA
## Sep 2021  120100         NA         NA
## Oct 2021  146137         NA         NA
## Nov 2021  150577         NA         NA
## Dec 2021  163619         NA         NA
## Jan 2022  143578         NA         NA
## Feb 2022   18455         NA         NA
## Mar 2022   40790         NA         NA
## Apr 2022  111057         NA         NA
## May 2022  212332         NA         NA
## Jun 2022  345438         NA         NA
## Jul 2022  476970         NA         NA
## Aug 2022  510246         NA         NA

# DMA Plot
ts.plot(df_dma[,1], xlab="periode waktu", ylab="Kujungan Wisatawan", col="blue", lty=1, ylim=c(-100000, 1600000))
points(df_dma[,1],col="black")
lines(df_dma[,2],col="#66FF66",lwd=2,lty=1)
lines(df_dma[,3],col="#FF3300",lwd= 2,lty=1)
title("Peramalan Kunjungan Wisatawan Mancanegara dengan Metode DMA n=2", cex.main=1, font.main=4 ,col.main="black")
legend("topleft", c("Data Aktual","Pemulusan","Ramalan"),lty=1,col=c ("blue","#66FF66","#FF3300"))
box(col="black",lwd=2)

#Akurasi data testing DMA
      df.master <- data.frame()
      df_ts_sma <- TTR::SMA(testing.ts,  n=2)
      df_ts_dma <- TTR::SMA(df_ts_sma, n=2)
      At <- 2*df_ts_sma-df_ts_dma
      Bt <- 2/(2-1)*(df_ts_sma-df_ts_dma)
      pemulusan_dma <- At+Bt
      ramal_dma <- c(NA, pemulusan_dma)
      df_dma <- cbind(df_aktual=c(testing.ts,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=2,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
                       MAD=MAD.dma,MAPE=MAPE.dma)
      df.master <- rbind(df.master,ak)
      df.master

##   n         SSE        MSE    RMSE      MAD     MAPE
## 1 2 66431040002 6039185455 77712.2 60954.55 103.0809

Double Exponential Smoothing

# Iterasi Mencari Nilai Alpha Beta Optimum
a = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
b = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
output1 = data.frame()
for (i in a) {
  for (j in b) {
    des1 <- HoltWinters(training.ts, alpha = i, beta=j, gamma=F)
    sse1 <- des1$SSE
    mse1 <- sse1/length(training.ts)
    rmse1 <-sqrt(mse1)
    akurasi1 <- cbind("SSE"=sse1, "MSE"=mse1, "RMSE"=rmse1)
    output1 <- rbind(output1, akurasi1)
  }
  
}
output_des1 <- cbind("Alpha" = rep(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9), each=9), "Beta"=b, output1)
kableExtra::kable(output_des1 ,caption = 'Mencari Parameter Optimum dalam Metode Double Exponential Smoothing',  align = rep('c',7))

Mencari Parameter Optimum dalam Metode Double Exponential Smoothing
Alpha	Beta	SSE	MSE	RMSE
0.1	0.1	8.300512e+12	81377569954	285267.5
0.1	0.2	6.319724e+12	61958074637	248913.8
0.1	0.3	6.305676e+12	61820353376	248637.0
0.1	0.4	6.536981e+12	64088046458	253156.2
0.1	0.5	6.719991e+12	65882268895	256675.4
0.1	0.6	6.717820e+12	65860981222	256633.9
0.1	0.7	6.578284e+12	64492981631	253954.7
0.1	0.8	6.510730e+12	63830685653	252647.4
0.1	0.9	6.580969e+12	64519304985	254006.5
0.2	0.1	4.092829e+12	40125771374	200314.2
0.2	0.2	4.017021e+12	39382559860	198450.4
0.2	0.3	4.063170e+12	39834998769	199587.1
0.2	0.4	4.077612e+12	39976589941	199941.5
0.2	0.5	4.149743e+12	40683759708	201702.2
0.2	0.6	4.278384e+12	41944944123	204804.6
0.2	0.7	4.394697e+12	43085267604	207569.9
0.2	0.8	4.425173e+12	43384049494	208288.4
0.2	0.9	4.341777e+12	42566442488	206316.4
0.3	0.1	2.851895e+12	27959753984	167211.7
0.3	0.2	2.835231e+12	27796382308	166722.5
0.3	0.3	2.844583e+12	27888065671	166997.2
0.3	0.4	2.865362e+12	28091779773	167606.0
0.3	0.5	2.865492e+12	28093062176	167609.9
0.3	0.6	2.814049e+12	27588719575	166098.5
0.3	0.7	2.721521e+12	26681580114	163345.0
0.3	0.8	2.619594e+12	25682297884	160257.0
0.3	0.9	2.533766e+12	24840846754	157609.8
0.4	0.1	2.193549e+12	21505384775	146647.1
0.4	0.2	2.175004e+12	21323572742	146025.9
0.4	0.3	2.167874e+12	21253670792	145786.4
0.4	0.4	2.149419e+12	21072734192	145164.5
0.4	0.5	2.110741e+12	20693538220	143852.5
0.4	0.6	2.064370e+12	20238921741	142263.6
0.4	0.7	2.026774e+12	19870336950	140962.2
0.4	0.8	2.005558e+12	19662328547	140222.4
0.4	0.9	2.001049e+12	19618125003	140064.7
0.5	0.1	1.801492e+12	17661684735	132897.3
0.5	0.2	1.785219e+12	17502149347	132295.7
0.5	0.3	1.774745e+12	17399459880	131907.0
0.5	0.4	1.758997e+12	17245070049	131320.5
0.5	0.5	1.742804e+12	17086316433	130714.6
0.5	0.6	1.735423e+12	17013951193	130437.5
0.5	0.7	1.741618e+12	17074681975	130670.1
0.5	0.8	1.761950e+12	17274023421	131430.7
0.5	0.9	1.795307e+12	17601046704	132668.9
0.6	0.1	1.554628e+12	15241449877	123456.3
0.6	0.2	1.546069e+12	15157544040	123116.0
0.6	0.3	1.544217e+12	15139381174	123042.2
0.6	0.4	1.544828e+12	15145373039	123066.5
0.6	0.5	1.552557e+12	15221143610	123374.0
0.6	0.6	1.571147e+12	15403405034	124110.5
0.6	0.7	1.601014e+12	15696211208	125284.5
0.6	0.8	1.640429e+12	16082636915	126817.3
0.6	0.9	1.686596e+12	16535253631	128589.5
0.7	0.1	1.393924e+12	13665919193	116901.3
0.7	0.2	1.395519e+12	13681558512	116968.2
0.7	0.3	1.405462e+12	13779036681	117384.1
0.7	0.4	1.421523e+12	13936497306	118052.9
0.7	0.5	1.446015e+12	14176621583	119065.6
0.7	0.6	1.479508e+12	14504983451	120436.6
0.7	0.7	1.520456e+12	14906434806	122091.9
0.7	0.8	1.566320e+12	15356082103	123919.7
0.7	0.9	1.614529e+12	15828717332	125812.2
0.8	0.1	1.288592e+12	12633252108	112397.7
0.8	0.2	1.300991e+12	12754811242	112937.2
0.8	0.3	1.322694e+12	12967584331	113875.3
0.8	0.4	1.351755e+12	13252504434	115119.5
0.8	0.5	1.388789e+12	13615580462	116685.8
0.8	0.6	1.433095e+12	14049952640	118532.5
0.8	0.7	1.483071e+12	14539916369	120581.6
0.8	0.8	1.537245e+12	15071028426	122764.1
0.8	0.9	1.594963e+12	15636893657	125047.6
0.9	0.1	1.222513e+12	11985424883	109478.0
0.9	0.2	1.246380e+12	12219414836	110541.5
0.9	0.3	1.280312e+12	12552074479	112036.0
0.9	0.4	1.322493e+12	12965613461	113866.6
0.9	0.5	1.372996e+12	13460748488	116020.5
0.9	0.6	1.431368e+12	14033022590	118461.1
0.9	0.7	1.497216e+12	14678591984	121155.2
0.9	0.8	1.570910e+12	15401079268	124101.1
0.9	0.9	1.653861e+12	16214325839	127335.5

Parameter optimum yang didapat untuk double exponential smoothing (DES) berdasarkan metrik ukuran kebaikan model di atas adalah alpha = 0.9 dan beta=0.1, sehingga metode pemulusan DES akan dilanjutkan dengan parameter tersebut.

df_des1 <- HoltWinters(training.ts, alpha = 0.9, beta=0.1, gamma=F)
head(df_des1$fitted)

##              xhat    level    trend
## Mar 2013 742502.0 678415.0 64087.00
## Apr 2013 789574.9 727034.6 62540.26
## May 2013 710091.8 660462.8 49629.05
## Jun 2013 750430.9 701646.4 48784.51
## Jul 2013 837986.9 785677.7 52309.19
## Aug 2013 771295.2 729804.3 41490.93

datades1 <- data.frame(training.ts, c(NA, NA, df_des1$fitted[,1]))
colnames(datades1) = c("y","yhat")
head(datades1)

##        y     yhat
## 1 614328       NA
## 2 678415       NA
## 3 725316 742502.0
## 4 646117 789574.9
## 5 700708 710091.8
## 6 789594 750430.9

tail(datades1)

##          y      yhat
## 97  126515 111381.03
## 98  105788  81384.18
## 99  119979  61926.54
## 100 112756  77977.40
## 101 139433  76211.86
## 102 126844 105734.51

# Hasil Peramalan
ramal_des1 <- forecast::forecast(df_des1,h=14) #periode ramalan sebanyak 14
(df_ramal_des1 <- data.frame(ramal_des1))

##          Point.Forecast       Lo.80    Hi.80      Lo.95     Hi.95
## Jul 2021       99256.53   -42577.06 241090.1  -117659.2  316172.3
## Aug 2021       73780.00  -125802.58 273362.6  -231455.2  379015.2
## Sep 2021       48303.48  -203286.27 299893.2  -336469.8  433076.8
## Oct 2021       22826.95  -278561.90 324215.8  -438107.5  483761.4
## Nov 2021       -2649.57  -353038.81 347739.7  -538523.7  533224.5
## Dec 2021      -28126.09  -427419.75 371167.6  -638793.0  582540.9
## Jan 2022      -53602.62  -502099.61 394894.4  -739519.6  632314.3
## Feb 2022      -79079.14  -577316.93 419158.6  -841068.1  682909.8
## Mar 2022     -104555.67  -653222.53 444111.2  -943669.2  734557.9
## Apr 2022     -130032.19  -729914.48 469850.1 -1047473.0  787408.6
## May 2022     -155508.72  -807457.40 496440.0 -1152578.2  841560.8
## Jun 2022     -180985.24  -885893.87 523923.4 -1259050.0  897079.5
## Jul 2022     -206461.76  -965251.56 552328.0 -1366930.6  954007.1
## Aug 2022     -231938.29 -1045547.70 581671.1 -1476246.5 1012369.9

# Plot Peramalan
plot(ramal_des1, xlab="Tahun", ylab="Jumlah Wisatawan",main="Double Exponential Smoothing (Alpha=0.9, Beta=0.1)")

# Gabungan Data Aktual, Pemulusan, dan Ramalan
data.des1 <- cbind(aktual=c(training.ts, rep(NA,14)),
                   pemulusan=c(NA, df_des1$fitted[,2], as.numeric(df_des1$coefficients[1]+df_des1$coefficients[2]), rep(NA,14)),
                   ramalan = c(NA, NA, df_des1$fitted[,1], df_ramal_des1$Point.Forecast))
data.des1 <- ts(data.des1, start=2013, freq=12)
data.des1

##           aktual  pemulusan    ramalan
## Jan 2013  614328         NA         NA
## Feb 2013  678415  678415.00         NA
## Mar 2013  725316  727034.60  742502.00
## Apr 2013  646117  660462.79  789574.86
## May 2013  700708  701646.38  710091.84
## Jun 2013  789594  785677.69  750430.89
## Jul 2013  717784  729804.29  837986.88
## Aug 2013  771009  771037.62  771295.22
## Sep 2013  770878  775040.48  812502.79
## Oct 2013  719903  729188.64  812759.42
## Nov 2013  807422  802534.85  758550.50
## Dec 2013  860655  858219.01  836295.15
## Jan 2014  753079  767188.27  894171.70
## Feb 2014  702666  711443.66  790442.61
## Mar 2014  765607  761726.11  726798.11
## Apr 2014  726332  731756.14  780573.36
## May 2014  752363  751698.87  745721.66
## Jun 2014  851475  842953.71  766262.11
## Jul 2014  777210  786007.61  865186.11
## Aug 2014  826821  824171.12  800322.16
## Sep 2014  791296  796253.46  840870.57
## Oct 2014  808767  808739.42  808491.19
## Nov 2014  764461  770115.10  821001.98
## Dec 2014  915334  901529.50  777288.97
## Jan 2015  724698  744340.94  921127.42
## Feb 2015  794302  789497.82  746260.22
## Mar 2015  792804  793097.69  795740.86
## Apr 2015  750999  755806.74  799076.41
## May 2015  794294  790610.45  757458.49
## Jun 2015  815307  813334.04  795577.40
## Jul 2015  815351  815823.57  820076.65
## Aug 2015  853244  850133.69  822140.87
## Sep 2015  870351  869240.93  859250.27
## Oct 2015  826196  831512.06  879356.58
## Nov 2015  777976  783862.73  836843.26
## Dec 2015  913828  900834.79  783895.87
## Jan 2016  814303  824128.88  912561.82
## Feb 2016  888309  882179.36  827012.63
## Mar 2016  915019  912575.08  890579.78
## Apr 2016  901095  903303.00  923175.02
## May 2016  915206  914876.97  911915.75
## Jun 2016  857651  864264.48  923785.84
## Jul 2016 1032741 1016189.02  867221.22
## Aug 2016 1031986 1032191.65 1034042.53
## Sep 2016 1006653 1010973.71 1049860.08
## Oct 2016 1040651 1039061.25 1024753.49
## Nov 2016 1002333 1007526.88 1054271.81
## Dec 2016 1113328 1103801.50 1018062.95
## Jan 2017 1107968 1109462.34 1122911.42
## Feb 2017 1023388 1033771.94 1127227.36
## Mar 2017 1059777 1058018.44 1042191.41
## Apr 2017 1171386 1161049.46 1068020.62
## May 2017 1148588 1151764.65 1180354.52
## Jun 2017 1144001 1146421.97 1168210.73
## Jul 2017 1370591 1349600.82 1160689.17
## Aug 2017 1393243 1392194.62 1382759.18
## Sep 2017 1250231 1267837.55 1426296.53
## Oct 2017 1161565 1174017.86 1286093.56
## Nov 2017 1062030 1073933.63 1181066.30
## Dec 2017 1147031 1139354.78 1070268.80
## Jan 2018 1097839 1102314.96 1142598.55
## Feb 2018 1197503 1187905.74 1101530.37
## Mar 2018 1363426 1346659.27 1195758.68
## Apr 2018 1302321 1309049.13 1369602.27
## May 2018 1242705 1251028.18 1325936.82
## Jun 2018 1322674 1316449.10 1260425.01
## Jul 2018 1547231 1525652.73 1331448.34
## Aug 2018 1511021 1515926.14 1560072.41
## Sep 2018 1370943 1388441.82 1545931.19
## Oct 2018 1291605 1302714.29 1402697.93
## Nov 2018 1157483 1172431.90 1306972.04
## Dec 2018 1405554 1381322.16 1163235.64
## Jan 2019 1201735 1220954.96 1393934.55
## Feb 2019 1243996 1241223.34 1216269.38
## Mar 2019 1311911 1304623.22 1239033.16
## Apr 2019 1274231 1277707.10 1308992.04
## May 2019 1249536 1252477.14 1278947.44
## Jun 2019 1434103 1415799.74 1251070.45
## Jul 2019 1468173 1464442.30 1430865.98
## Aug 2019 1530268 1525527.82 1482866.16
## Sep 2019 1388719 1404668.88 1548217.85
## Oct 2019 1346434 1353091.00 1413004.02
## Nov 2019 1280781 1288246.38 1355434.84
## Dec 2019 1377067 1367747.44 1283871.37
## Jan 2020 1290411 1298545.90 1371760.03
## Feb 2020  872765  915012.21 1295237.09
## Mar 2020  486155  524907.59  873680.90
## Apr 2020  158066  187129.30  448698.95
## May 2020  161842  154134.17   84763.69
## Jun 2020  156561  146775.46   58705.61
## Jul 2020  155742  146183.19   60153.89
## Aug 2020  161549  152210.55   68164.55
## Sep 2020  148984  142345.25   82596.51
## Oct 2020  152293  144934.31   78706.09
## Nov 2020  144476  138820.20   87917.96
## Dec 2020  164079  156360.51   86894.08
## Jan 2021  126515  125001.60  111381.03
## Feb 2021  105788  103347.62   81384.18
## Mar 2021  119979  114173.75   61926.54
## Apr 2021  112756  109278.14   77977.40
## May 2021  139433  133110.89   76211.86
## Jun 2021  126844   99256.53  105734.51
## Jul 2021      NA         NA   99256.53
## Aug 2021      NA         NA   73780.00
## Sep 2021      NA         NA   48303.48
## Oct 2021      NA         NA   22826.95
## Nov 2021      NA         NA   -2649.57
## Dec 2021      NA         NA  -28126.09
## Jan 2022      NA         NA  -53602.62
## Feb 2022      NA         NA  -79079.14
## Mar 2022      NA         NA -104555.67
## Apr 2022      NA         NA -130032.19
## May 2022      NA         NA -155508.72
## Jun 2022      NA         NA -180985.24
## Jul 2022      NA         NA -206461.76
## Aug 2022      NA         NA -231938.29

# Plot Hasil Pemulusan DES
ts.plot(data.des1[,1], xlab="Tahun", ylab="Jumlah Wisatawan", col="blue", lty=1)
points(data.des1[,1],col="black")
lines(data.des1[,2], col="#66FF66",lwd=2) #nilai dugaan
lines(data.des1[,3], col="#FF3300",lwd=2) #nilai dugaan
title("Double Exponential Smoothing (Alpha=0.9, Beta=0.1)",cex.main=1, font.main=4, col.main="black")
box(col="black",lwd=2)
legend("topleft", c("Data aktual", "Fitted DES","Ramalan"), lty=1,col=c ("blue","#66FF66","#FF3300"))

Berdasarkan eksplorasi di atas, plot data training hasil pemulusan menggunakan metode pemulusan DES dengan α = 0.9 dan β = 0.1 memiliki pola hampir sama dengan plot data aktual dibandingkan dengan data peramalannya.

# Akurasi data testing DES
selisihdes<-ramal_des1$mean-testing.ts
SSEtestingdes<-sum(selisihdes^2)
MSEtestingdes<-mean(selisihdes^2)
RMSEtestingdes<-sqrt(mean(selisihdes^2))
MADtestingdes<-mean(abs(selisihdes^2))
r.error.des <- (selisihdes/testing.ts)*100 
      MAPEtestingdes <- mean(abs(r.error.des), na.rm = T)
akDES <- data.frame(SSE=SSEtestingdes,MSE=MSEtestingdes,RMSE=RMSEtestingdes,
                       MAD=MADtestingdes, MAPE=MAPEtestingdes)
akDES

##            SSE          MSE     RMSE          MAD     MAPE
## 1 827444879247 103430609906 321606.3 103430609906 108.0806

Robust Exponential Smoothing

library(forecast)
library(robets)
model <- robets(training.ts, model = "ZZZ")
summary(model)

## ROBETS(A,N,A) 
## 
## Call:
##  robets(y = training.ts, model = "ZZZ") 
## 
##   Smoothing parameters:
##     alpha = 0.4512 
##     gamma = 0.1968 
## 
##   Initial states:
##     sigma = 128890.5723 
##     l = 818737.5833 
##     s=96596.42 -11315.58 7458.417 51613.42 34506.42 -3386.583
##            32737.42 -24443.58 -67738.58 -25933.58 -24435.58 -65658.58
## 
##   sigma:  182826.4
## 
##   robAIC  robAICc   robBIC 
## 2732.796 2732.918 2738.046 
## Training set: non robust error measures: 
##            ME          RMSE           MAE           MPE          MAPE 
## -31624.926504 182826.414613  90593.890294     -2.454522     12.290675 
## Training set: robust error measures: 
##      MedianE         RTSE        RTSPE 
## 4.175765e+03 4.074174e+09 7.239773e+01

##            ME          RMSE           MAE           MPE          MAPE 
## -3.162493e+04  1.828264e+05  9.059389e+04 -2.454522e+00  1.229068e+01 
##       MedianE          RTSE         RTSPE 
##  4.175765e+03  4.074174e+09  7.239773e+01

# Hasil Peramalan
ramalrobets <- forecast(model, h=14)
ramalrobets

##          Point Forecast       Lo 80    Hi 80      Lo 95    Hi 95
## Jul 2021       193812.5   45066.672 342558.4  -33674.63 421299.7
## Aug 2021       231052.3   67867.316 394237.3  -18517.60 480622.2
## Sep 2021       180203.5    3757.079 356650.0  -89648.04 450055.1
## Oct 2021       162188.2  -26590.448 350966.8 -126523.81 450900.1
## Nov 2021       124515.6  -75837.576 324868.7 -181898.12 430929.2
## Dec 2021       244647.5   33352.926 455942.0  -78499.65 567794.6
## Jan 2022       154198.6  -67498.009 375895.3 -184857.12 493254.4
## Feb 2022       129925.4 -101706.632 361557.5 -224325.23 484176.1
## Mar 2022       162055.4  -79103.141 403213.8 -206764.73 530875.4
## Apr 2022       119932.7 -130389.981 370255.3 -262902.79 502768.1
## May 2022       100497.9 -167741.269 368737.0 -309738.50 510734.2
## Jun 2022       137120.6 -139386.656 413627.8 -285760.74 560001.9
## Jul 2022       193812.5  -90722.625 478347.7 -241346.44 628971.5
## Aug 2022       231052.3  -61290.413 523395.0 -216047.31 678151.9

# Plot Peramalan
plot(ramalrobets)

# Akurasi data testing RES
selisihres<-ramalrobets$mean-testing.ts
SSEtestingres<-sum(selisihres^2)
MSEtestingres<-mean(selisihres^2)
RMSEtestingres<-sqrt(mean(selisihres^2))
MADtestingres<-mean(abs(selisihres^2))
r.error.res <- (selisihres/testing.ts)*100 
      MAPEtestingres <- mean(abs(r.error.res), na.rm = T)
akRES <- data.frame(SSE=SSEtestingres,MSE=MSEtestingres,RMSE=RMSEtestingres,
                       MAD=MADtestingres, MAPE=MAPEtestingres)
akRES

##            SSE         MSE     RMSE         MAD     MAPE
## 1 336467223278 42058402910 205081.5 42058402910 223.4415

Winter’s Aditif

# Pemulusan 
aditif <- HoltWinters(training.ts, seasonal = "additive")
aditif

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0
##  gamma: 0
## 
## Coefficients:
##           [,1]
## a    60940.979
## b     4278.619
## s1  -23074.063
## s2   23359.188
## s3   20538.937
## s4  -35457.146
## s5   46567.271
## s6   95069.604
## s7  -17560.854
## s8  -72775.437
## s9  -13010.687
## s10 -56839.104
## s11 -32720.729
## s12  65903.021

# Forecasting
ramalan1 <- forecast(aditif, h=14)
ramalan1

##          Point Forecast     Lo 80    Hi 80     Lo 95     Hi 95
## Jul 2021       42145.54 -109977.0 194268.1 -190505.8  274796.8
## Aug 2021       92857.40 -122276.3 307991.1 -236161.2  421876.0
## Sep 2021       94315.77 -169168.2 357799.7 -308648.1  497279.7
## Oct 2021       42598.31 -261646.7 346843.3 -422704.3  507900.9
## Nov 2021      128901.34 -211254.9 469057.6 -391322.8  649125.5
## Dec 2021      181682.29 -190940.3 554304.8 -388194.7  751559.3
## Jan 2022       73330.46 -329147.9 475808.8 -542207.1  688868.0
## Feb 2022       22394.49 -407873.0 452661.9 -635642.8  680431.8
## Mar 2022       86437.86 -369929.7 542805.4 -611516.1  784391.8
## Apr 2022       46888.06 -434165.6 527941.7 -688820.0  782596.1
## May 2022       75285.05 -429248.3 579818.4 -696332.1  846902.2
## Jun 2022      178187.42 -348780.4 705155.3 -627740.4  984115.2
## Jul 2022       93488.96 -454996.6 641974.5 -745347.3  932325.2
## Aug 2022      144200.83 -424989.5 713391.2 -726300.7 1014702.3

# Akurasi data training
sse1.train <- aditif$SSE
sse1.train

## [1] 1.266537e+12

# Plot Pemulusan
plot(ramalan1, xlab="Tahun", ylab="Jumlah Kunjungan")

Winter’s Damped Multiplikatif

# Pemulusan 
multi <- HoltWinters(training.ts, seasonal = "multiplicative")
multi

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, seasonal = "multiplicative")
## 
## Smoothing parameters:
##  alpha: 0.8255676
##  beta : 0.05675182
##  gamma: 1
## 
## Coefficients:
##              [,1]
## a    1.190900e+05
## b   -2.813304e+04
## s1   1.132100e+00
## s2   1.128662e+00
## s3   1.043335e+00
## s4   1.023979e+00
## s5   9.923898e-01
## s6   1.161302e+00
## s7   1.039876e+00
## s8   9.551069e-01
## s9   9.323945e-01
## s10  7.416652e-01
## s11  9.646528e-01
## s12  1.065110e+00

# Forecasting
ramalan2 <- holt(training.ts, damped = TRUE, seasonal = "multiplicative", h=14)
ramalan2

##          Point Forecast      Lo 80    Hi 80     Lo 95     Hi 95
## Jul 2021       125880.1  -12491.63 264251.8  -85741.2  337501.4
## Aug 2021       125108.0  -77609.28 327825.3 -184921.3  435137.3
## Sep 2021       124490.3 -131290.02 380270.6 -266691.9  515672.5
## Oct 2021       123996.2 -178900.76 426893.1 -339244.7  587237.0
## Nov 2021       123600.8 -222406.38 469608.0 -405571.5  652773.2
## Dec 2021       123284.6 -262788.98 509358.1 -467164.0  713733.1
## Jan 2022       123031.6 -300635.31 546698.4 -524911.0  770974.1
## Feb 2022       122829.2 -336340.20 581998.5 -579409.8  825068.1
## Mar 2022       122667.2 -370192.28 615526.7 -631096.3  876430.8
## Apr 2022       122537.7 -402415.13 647490.5 -680308.4  925383.8
## May 2022       122434.1 -433189.09 678057.2 -727318.2  972186.3
## Jun 2022       122351.2 -462663.81 707366.1 -772352.0 1017054.4
## Jul 2022       122284.8 -490966.00 735535.7 -815601.4 1060171.1
## Aug 2022       122231.8 -518204.47 762668.0 -857230.9 1101694.5

# Akurasi data training
sse2.train <- multi$SSE
sse2.train

## [1] 1.043219e+12

# Plot Pemulusan
plot(ramalan2, xlab="Tahun", ylab="Jumlah Kunjungan")

selisih1<-as.numeric(ramalan1$mean)-as.numeric(testing.ts)
SSEtesting1<-sum(selisih1^2)
MSEtesting1<-mean(selisih1^2)
RMSEtesting1<-sqrt(mean(selisih1^2))
MADtesting1<-mean(abs(selisih1^2))
r.error.testing1 <- (selisih1/testing.ts)*100 
      MAPEtesting1 <- mean(abs(r.error.testing1), na.rm = T)

selisih2<-as.numeric(ramalan2$mean)-as.numeric(testing.ts)
SSEtesting2<-sum(selisih2^2)
MSEtesting2<-mean(selisih2^2)
RMSEtesting2<-sqrt(mean(selisih2^2))
MADtesting2<-mean(abs(selisih2^2))
r.error.testing2 <- (selisih1/testing.ts)*100 
      MAPEtesting2 <- mean(abs(r.error.testing2), na.rm = T)


akurasi <- matrix(c(SSEtesting1, MSEtesting1, RMSEtesting1, MADtesting1, MAPEtesting1, SSEtesting2,  MSEtesting2,  RMSEtesting2, MADtesting2, MAPEtesting2), nrow=5, ncol=2)
row.names(akurasi)<- c("SSE", "MSE", "RMSE", "MAD", "MAPE")
colnames(akurasi) <- c("Aditif", "Damped Multiplikatif")
akurasi

##            Aditif Damped Multiplikatif
## SSE  3.590361e+11         3.552670e+11
## MSE  2.564543e+10         2.537621e+10
## RMSE 1.601419e+05         1.592991e+05
## MAD  2.564543e+10         2.537621e+10
## MAPE 5.081117e+01         5.081117e+01

#metode Winter yang dipilih : Damped Multiplikatif

Perbandingan Nilai Akurasi

akurasi <- matrix(c(SSE.dma, MSE.dma, RMSE.dma, MAD.dma, MAPE.dma, SSEtestingdes,  MSEtestingdes,  RMSEtestingdes,  MADtestingdes, MAPEtestingdes, SSEtestingres,  MSEtestingres,  RMSEtestingres,  MADtestingres, MAPEtestingres, SSEtesting1, MSEtesting1, RMSEtesting1, MADtesting1, MAPEtesting1, SSEtesting2,  MSEtesting2,  RMSEtesting2, MADtesting2, MAPEtesting2), nrow=5)
row.names(akurasi)<- c("SSE", "MSE", "RMSE", "MAD", "MAPE")
colnames(akurasi) <- c("DMA", "DES", "RES", "Winter Aditif", "Winter Damped Multiplikatif")
kableExtra::kable(akurasi,caption = 'Penentuan Model Smoothing Terbaik',  align = rep('c',7))

Penentuan Model Smoothing Terbaik
	DMA	DES	RES	Winter Aditif	Winter Damped Multiplikatif
SSE	6.643104e+10	8.274449e+11	3.364672e+11	3.590361e+11	3.552670e+11
MSE	6.039185e+09	1.034306e+11	4.205840e+10	2.564543e+10	2.537621e+10
RMSE	7.771220e+04	3.216063e+05	2.050815e+05	1.601419e+05	1.592991e+05
MAD	6.095455e+04	1.034306e+11	4.205840e+10	2.564543e+10	2.537621e+10
MAPE	1.030809e+02	1.080806e+02	2.234415e+02	5.081117e+01	5.081117e+01

Berdasarkan data testing, model Holt Winter Damped Multiplicative merupakan model terbaik dilihat dari perbandingan ukuran kebaikan model di atas.

Pemodelan

Pengecekan Kestasioneran Data

Stasioner Terhadap Rataan

par(mfrow=c(1,1))
ggAcf(training.ts, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()

ggPacf(training.ts, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()

Berdasarkan plot time series, terlihat data memiliki pola gabungan. Pada periode pra-covid terdapat pola musiman (fluktuatif) pada rata-rata dan lebar pita setiap waktu cukup berbeda, sehingga mengindikasikan data tidak stasioner pada rata-rata dan ragam. Kemudian, terjadi penurunan yang ekstrem akibat adanya pandemi Covid-19. Setelah pandemi, kunjungan wisatawan terlihat memiliki pola yang cenderung stasioner.

Sedangkan pada plot ACF, terlihat bahwa plot membentuk pola tails off (slowly), sehingga mengindikasikan bahwa data tersebut tidak stasioner.

Uji Formal

Uji Augmented Dickey-Fuller (ADF) merupakan salah satu uji yang dapat dilakukan untuk mengetahui apakah suatu data deret waktu stasioner atau tidak. Uji ADF telah mempertimbangkan kemungkinan adanya autokorelasi pada error term jika series yang digunakan nonstasioner (Nkoro dan Uko, 2016).

\[ \begin{aligned} H_0&: \rho = 0 \text{ (Data Tidak Stasioner)}\\ H_1&: \rho \neq 0 \text{ (Data Stasioner)} \end{aligned} \]

adf.test(training.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.ts
## Dickey-Fuller = -1.0461, Lag order = 4, p-value = 0.9272
## alternative hypothesis: stationary

Hasil dari uji formal ADF di atas diperoleh nilai \(P-value>0.05\) sehingga tak tolak \(H_0\) pada taraf nyata 5%. Hal tersebut berarti cukup bukti untuk mengatakan bahwa data tidak stasioner. Kemudian perlu dilakukan differencing atau pembedaan untuk mengatasi ketidakstasioneran pada data.

Stasioner Terhadap Ragam

Untuk menguji kestasioneran terhadap ragam, dilakukan uji Langrage multiplier test atau Uji Arch dengan hipotesis sebagai berikut:

\[ \begin{aligned} H_0&: \text{ Data Stasioner dalam Ragam (Homoskedastisitas)}\\ H_1&: \text{ Data Tidak Stasioner dalam Ragam (Heteroskedastisitas)} \end{aligned} \]

ArchTest(training.ts)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  training.ts
## Chi-squared = 79.476, df = 12, p-value = 5.196e-12

Hasil dari uji formal Arch di atas diperoleh nilai \(P-value<0.05\) sehingga tolak \(H_0\) pada taraf nyata 5%. Hal tersebut berarti cukup bukti untuk mengatakan bahwa data tidak stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

Differencing Pertama

training.dif <- diff(training.ts, differences = 1)
training.dif

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 2013           64087   46901  -79199   54591   88886  -71810   53225    -131
## 2014 -107576  -50413   62941  -39275   26031   99112  -74265   49611  -35525
## 2015 -190636   69604   -1498  -41805   43295   21013      44   37893   17107
## 2016  -99525   74006   26710  -13924   14111  -57555  175090    -755  -25333
## 2017   -5360  -84580   36389  111609  -22798   -4587  226590   22652 -143012
## 2018  -49192   99664  165923  -61105  -59616   79969  224557  -36210 -140078
## 2019 -203819   42261   67915  -37680  -24695  184567   34070   62095 -141549
## 2020  -86656 -417646 -386610 -328089    3776   -5281    -819    5807  -12565
## 2021  -37564  -20727   14191   -7223   26677  -12589                        
##          Oct     Nov     Dec
## 2013  -50975   87519   53233
## 2014   17471  -44306  150873
## 2015  -44155  -48220  135852
## 2016   33998  -38318  110995
## 2017  -88666  -99535   85001
## 2018  -79338 -134122  248071
## 2019  -42285  -65653   96286
## 2020    3309   -7817   19603
## 2021

par(mfrow=c(1,1))
plot.ts(training.dif, lty=1, ylab=expression(Y[t]))

Berdasarkan plot, secara eksploratif dapat diamati bahwa data hasil differencing pertama menjadi lebih stasioner pada rataan, tetapi masih kurang stasioner pada ragam. Untuk menghasilkan kesimpulan yang valid, perlu dilakukan pengecekan ulang kestasioneran data menggunakan uji formal sebagai berikut,

adf.test(training.dif)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.dif
## Dickey-Fuller = -4.0294, Lag order = 4, p-value = 0.01051
## alternative hypothesis: stationary

ArchTest(training.dif)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  training.dif
## Chi-squared = 31.583, df = 12, p-value = 0.001604

Berdasarkan hasil uji formal di atas, dapat diamati bahwa data hasil differencing dengan paramater d = 1 sukses menstasionerkan data terhadap rataan, tetapi belum terhadap ragam sehingga perlu dilakukan proses diferensiasi tahap kedua.

Differencing Kedua

training.diff <- diff(training.ts, differences = 2)
training.diff

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 2013                  -17186 -126100  133790   34295 -160696  125035  -53356
## 2014 -160809   57163  113354 -102216   65306   73081 -173377  123876  -85136
## 2015 -341509  260240  -71102  -40307   85100  -22282  -20969   37849  -20786
## 2016 -235377  173531  -47296  -40634   28035  -71666  232645 -175845  -24578
## 2017 -116355  -79220  120969   75220 -134407   18211  231177 -203938 -165664
## 2018 -134193  148856   66259 -227028    1489  139585  144588 -260767 -103868
## 2019 -451890  246080   25654 -105595   12985  209262 -150497   28025 -203644
## 2020 -182942 -330990   31036   58521  331865   -9057    4462    6626  -18372
## 2021  -57167   16837   34918  -21414   33900  -39266                        
##          Oct     Nov     Dec
## 2013  -50844  138494  -34286
## 2014   52996  -61777  195179
## 2015  -61262   -4065  184072
## 2016   59331  -72316  149313
## 2017   54346  -10869  184536
## 2018   60740  -54784  382193
## 2019   99264  -23368  161939
## 2020   15874  -11126   27420
## 2021

Setelah dilakukan differencing satu kali masih ditemukan ketidakstasioneran pada data sehingga dilakukan sekali lagi differencing.

Pengecekan Ulang Kestasioneran Data

par(mfrow=c(1,1))
plot.ts(training.diff, lty=1, ylab=expression(Y[t]))

Berdasarkan plot, secara eksploratif dapat diamati bahwa data hasil differencing lebih stasioner, baik itu terhadap rataan maupun ragam. Untuk menghasilkan kesimpulan yang valid, perlu dilakukan pengecekan ulang kestasioneran data menggunakan uji formal sebagai berikut,

adf.test(training.diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.diff
## Dickey-Fuller = -7.3044, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

ArchTest(training.diff)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  training.diff
## Chi-squared = 20.829, df = 12, p-value = 0.05295

Berdasarkan hasil uji formal di atas, dapat diamati bahwa data hasil differencing dengan paramater d = 2 sukses menstasionerkan data terhadap rataan dan ragamnya.

Identifikasi Model ARIMA

Plot ACF dan PACF

Menurut Tinungki (2018), Autocorrelation Function (ACF) adalah fungsi yang menunjukkan korelasi antara suatu amatan ke-\(t\) [\(y_t\)] dengan amatan sebelumnya [\(y_{t-1}\)], sekaligus mengunjukkan koefisien autokorelasi antar amatan.

Sementara itu, Partial Autocorrelation Function (PACF) adalah fungsi yang menunjukkan korelasi antara suatu amatan dengan amatan sebelumnya setelah amatan-amatan yang mengintervensi dihapuskan.

par(mfrow=c(1,1))
ggAcf(training.diff, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()

ggPacf(training.diff, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()

Dari plot ACF dan PACF di atas, terlihat bahwa terdapat cuts off pada lag ke=1 pada kedua plot. Maka dari itu, dapat dideteksi dua potensi model data deret waktu, yaitu adalah ARIMA(0,2,1) dan ARIMA(1,2,0).

Plot EACF

Extended Autocorrelation Function (EACF) merupakan metode grafis untuk mengidentifikasi model ARMA. Dalam tabel EACF, proses ARMA (p, q) memiliki pola teoritis segitiga nol dengan bagian atas simpul kiri yang sesuai ordo ARMA.

eacf(training.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  x  x  o 
## 1 x o o o o o o o o o o  x  x  x 
## 2 x o o o o o o o o o o  x  x  x 
## 3 x x x x o o o o o o o  x  o  o 
## 4 x o o x o o o o o o o  o  o  o 
## 5 x o o x o o o o o o o  o  o  o 
## 6 o o o x o o o o o o o  o  o  o 
## 7 o x o o o o o o o o o  x  o  o

Berdasarkan tabel EACF di atas, dapat teramati adanya beberapa kemungkinan terbentuknya model ARIMA berdasarkan pola triangle of zero, antara lain model ARIMA(4,2,4) dan ARIMA(3,2,4).

Pemilihan Model Terbaik

Berdasarkan hasil sebelumnya diperoleh beberapa model tentatif yaitu ARIMA(0,2,1), ARIMA(1,2,0), ARIMA(4,2,4), ARIMA(3,2,4). Selanjutnya akan dilakukan pemilihan model terbaik berdasarkan nilai AIC.

# Perbandingan Kebaikan Model Tentatif
#Melihat summary dan signifikansi parameter dari setiap model
model1 <- Arima(training.ts, order=c(0,2,1), method="ML")
coeftest(model1)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.956875   0.039771  -24.06 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ramalan1 <- forecast::forecast(training.ts,model=model1,h=14)
data.ramalan1 <- ramalan1$mean
MAPE.1 <- mean(abs((testing.ts-data.ramalan1)/testing.ts),na.rm = T)*100

model2 <- Arima(training.ts, order=c(1,2,0), method="ML")
coeftest(model2)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.445376   0.088747 -5.0185 5.208e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ramalan2 <- forecast::forecast(training.ts,model=model2,h=14)
data.ramalan2 <- ramalan2$mean

model3 <- Arima(training.ts, order=c(4,2,4), method="ML")
coeftest(model3)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.691620   0.807077 -0.8569 0.3914756    
## ar2 -0.960806   0.287103 -3.3466 0.0008182 ***
## ar3 -0.229613   0.658170 -0.3489 0.7271903    
## ar4 -0.046831   0.134129 -0.3492 0.7269745    
## ma1 -0.139147   0.803174 -0.1732 0.8624576    
## ma2  0.368987   0.402787  0.9161 0.3596217    
## ma3 -0.773820   0.453221 -1.7074 0.0877515 .  
## ma4 -0.341228   0.798492 -0.4273 0.6691311    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ramalan3 <- forecast::forecast(training.ts,model=model3,h=14)
data.ramalan3 <- ramalan3$mean

model4 <- Arima(training.ts, order=c(3,2,4), method="ML")
coeftest(model4)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.883566   0.542231 -1.6295   0.10321    
## ar2 -1.000914   0.231966 -4.3149 1.597e-05 ***
## ar3 -0.366850   0.481899 -0.7613   0.44650    
## ma1  0.058928   0.515228  0.1144   0.90894    
## ma2  0.266394   0.257034  1.0364   0.30001    
## ma3 -0.666180   0.306503 -2.1735   0.02974 *  
## ma4 -0.539430   0.507065 -1.0638   0.28741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ramalan4 <- forecast::forecast(training.ts,model=model4,h=14)
data.ramalan4 <- ramalan4$mean

Model <- c("ARIMA (0,2,1)","ARIMA (1,2,0)","ARIMA (4,2,4)","ARIMA (3,2,4)")
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, align = rep('c',7))

Model	AIC	BIC
ARIMA (0,2,1)	2607.030	2612.240
ARIMA (1,2,0)	2636.722	2641.932
ARIMA (4,2,4)	2613.148	2636.594
ARIMA (3,2,4)	2611.245	2632.086

Dari output di atas, diketahui bahwa pada model ARIMA(0,2,1), ARIMA(1,2,0) seluruh parameter signifikan, sementara pada model lainnya mayoritas parameter tidak signifikan. Dalam penetuan model terbaik, selain dilihat dari signifikasi pendugaan parameter model juga perlu memperhatikan nilai AIC nya. Dari kedua model tentatif tersebut, model ARIMA(0,2,1) memiliki nilai AIC terkecil. Oleh karena itu, dugaan model terbaik sementara adalah ARIMA (0,2,1).

Overfitting Model

Selanjutnya dilakukan overfitting model ARIMA(0,2,1) yaitu ARIMA(1,2,1) dan ARIMA(0,2,2). Kedua model tentatif tersebut dibandingkan dengan model ARIMA (0,2,1).

# Perbandingan Kebaikan Model Tentatif
#Melihat summary dan signifikansi parameter dari setiap model
model5 <- Arima(training.ts, order=c(1,2,1),method="ML") 
coeftest(model5)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.116698   0.103690   1.1255   0.2604    
## ma1 -0.968868   0.036772 -26.3482   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model6 <- Arima(training.ts, order=c(0,2,2),method="ML") 
coeftest(model6)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)    
## ma1 -0.856423   0.097855 -8.7520   <2e-16 ***
## ma2 -0.107150   0.096930 -1.1054    0.269    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model <- c("ARIMA (0,2,1)", "ARIMA (1,2,1)","ARIMA (0,2,2)")
AIC <- c(model1$aic,model5$aic,model6$aic)
BIC <- c(model1$bic,model5$bic,model6$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi, align = rep('c',7))

Model	AIC	BIC
ARIMA (0,2,1)	2607.030	2612.240
ARIMA (1,2,1)	2607.766	2615.581
ARIMA (0,2,2)	2607.820	2615.636

Berdasarkan output uji signifikansi di atas, kedua model overfitting ARIMA (0,2,1), yaitu ARIMA (1,2,1) dan ARIMA (0,2,2)ada beberapa parameter yang tidak signifikan. Berdasarkan hasil perbandingan nilai AIC, model ARIMA(0,2,1) merupakan model terbaik sementara karena memiliki parameter yang seluruhnya signifikan dengan nilai AIC paling minimum.

Identifikasi Model Seasonal ARIMA

Model ARIMA(0,2,1)(1,0,0)[12]

model.s3 <- training.ts %>%
  Arima(order = c(0,2,1),
        seasonal = c(1,0,0))
fc.3 <- model.s3 %>%
          forecast(h=14)
MAPE.s3 <- mean(abs((testing.ts-fc.3$mean)/testing.ts),na.rm = T)*100

AIC_model.s3 <- model.s3$aic
BIC_model.s3 <- model.s3$bic
coeftest(model.s3)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.975713   0.043884 -22.2341 < 2.2e-16 ***
## sar1  0.298706   0.093239   3.2037  0.001357 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(0,2,1)(0,0,1)[12]

model.s4 <- training.ts %>%
  Arima(order = c(0,2,1),
        seasonal = c(0,0,1))
fc.4 <- model.s4 %>%
          forecast(h=14)
MAPE.s4 <- mean(abs((testing.ts-fc.4$mean)/testing.ts),na.rm = T)*100

AIC_model.s4 <- model.s4$aic
BIC_model.s4 <- model.s4$bic
coeftest(model.s4)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.975786   0.045925 -21.2476 < 2.2e-16 ***
## sma1  0.311258   0.103813   2.9983  0.002715 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(4,2,0)(1,0,0)[12]

model.s5 <- training.ts %>%
  Arima(order = c(4,2,0),
        seasonal = c(1,0,0))
fc.5 <- model.s5 %>%
          forecast(h=14)
MAPE.s5 <- mean(abs((testing.ts-fc.5$mean)/testing.ts),na.rm = T)*100

AIC_model.s5 <- model.s5$aic
BIC_model.s5 <- model.s5$bic
coeftest(model.s5)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.649367   0.099317 -6.5383  6.22e-11 ***
## ar2  -0.442252   0.115591 -3.8260 0.0001302 ***
## ar3  -0.341392   0.114111 -2.9918 0.0027738 ** 
## ar4  -0.179345   0.098469 -1.8213 0.0685568 .  
## sar1  0.356613   0.092870  3.8399 0.0001231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(4,2,0)(0,0,1)[12]

model.s7 <- training.ts %>%
  Arima(order = c(4,2,0),
        seasonal = c(0,0,1))
fc.7 <- model.s7 %>%
          forecast(h=14)
MAPE.s7 <- mean(abs((testing.ts-fc.7$mean)/testing.ts),na.rm = T)*100

AIC_model.s7 <- model.s7$aic
BIC_model.s7 <- model.s7$bic
coeftest(model.s7)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.660756   0.098307 -6.7213 1.801e-11 ***
## ar2  -0.443423   0.115307 -3.8456 0.0001203 ***
## ar3  -0.362024   0.114402 -3.1645 0.0015536 ** 
## ar4  -0.190627   0.097734 -1.9505 0.0511204 .  
## sma1  0.405998   0.111401  3.6445 0.0002679 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perbandingan Model Tentatif Terbaik

Model.s1 <- c("ARIMA(0,2,1)(1,0,0)[12]","ARIMA(0,2,1)(0,0,1)[12]")
AIC.s1 <- round(c(AIC_model.s3, AIC_model.s4),2)
BIC.s1 <- round(c(BIC_model.s3, BIC_model.s4),2)
perbandingan.aic <- data.frame(cbind(Model=Model.s1, AIC=AIC.s1, BIC=BIC.s1))
kableExtra::kable(perbandingan.aic, align=rep("c",7), caption = "Perbandingan Model Tentatif Terbaik")

Perbandingan Model Tentatif Terbaik
Model	AIC	BIC
ARIMA(0,2,1)(1,0,0)[12]	2599.35	2607.17
ARIMA(0,2,1)(0,0,1)[12]	2599.61	2607.43

Overfitting Model Seasonal ARIMA

Model ARIMA(1,2,1)(1,0,0)[12]

model.s8 <- training.ts %>%
              Arima(order = c(1,2,1),
                    seasonal = c(1,0,0))
fc.8 <- model.s8 %>%
          forecast(h=14)
MAPE.s8 <- mean(abs((testing.ts-fc.8$mean)/testing.ts),na.rm = T)*100

AIC_model.s8 <- model.s8$aic
BIC_model.s8 <- model.s8$bic
coeftest(model.s8)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.213492   0.100659   2.1209  0.033927 *  
## ma1  -0.999995   0.056534 -17.6885 < 2.2e-16 ***
## sar1  0.348313   0.091756   3.7961  0.000147 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(0,2,2)(1,0,0)[12]

model.s9 <- training.ts %>%
              Arima(order = c(0,2,2),
                    seasonal = c(1,0,0))
fc.9 <- model.s9 %>%
          forecast(h=14)
MAPE.s9 <- mean(abs((testing.ts-fc.9$mean)/testing.ts),na.rm = T)*100

AIC_model.s9 <- model.s9$aic
BIC_model.s9 <- model.s9$bic
coeftest(model.s9)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.822857   0.223343 -3.6843 0.0002294 ***
## ma2  -0.175339   0.101449 -1.7283 0.0839269 .  
## sar1  0.335853   0.091849  3.6566 0.0002556 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perbandingan Model Overfitting Terbaik

Model.s2 <- c("ARIMA(0,2,1)(1,0,0)[12]","ARIMA(1,2,1)(1,0,0)[12]")
AIC.s2 <- round(c(AIC_model.s3, AIC_model.s8),2)
BIC.s2 <- round(c(BIC_model.s3, BIC_model.s8),2)
MAPE.s2 <- round(c(MAPE.s3, MAPE.s8),2)
perbandingan.aic2 <- data.frame(cbind(Model=Model.s2, AIC=AIC.s2, BIC=BIC.s2, MAPE=MAPE.s2))
kableExtra::kable(perbandingan.aic2, align = rep("c",7), caption = "Perbandingan Model Overfitting Terbaik")

Perbandingan Model Overfitting Terbaik
Model	AIC	BIC	MAPE
ARIMA(0,2,1)(1,0,0)[12]	2599.35	2607.17	123.95
ARIMA(1,2,1)(1,0,0)[12]	2597.1	2607.52	129.69

Perbandingan Model ARIMA dan Seasonal ARIMA

Model.terbaik <- c("ARIMA (0,2,1)","ARIMA(0,2,1)(1,0,0)[12]")
AIC.s2 <- round(c(model1$aic, AIC_model.s3),2)
BIC.s2 <- round(c(model1$bic, BIC_model.s3),2)
MAPE.terbaik <- round(c(MAPE.1, MAPE.s3),2)
perbandingan.aic2 <- data.frame(cbind(Model=Model.terbaik, AIC=AIC.s2, BIC=BIC.s2,
                                      MAPE=MAPE.terbaik))
kableExtra::kable(perbandingan.aic2, align = rep("c",7), caption = "Perbandingan Model Terbaik")

Perbandingan Model Terbaik
Model	AIC	BIC	MAPE
ARIMA (0,2,1)	2607.03	2612.24	110.25
ARIMA(0,2,1)(1,0,0)[12]	2599.35	2607.17	123.95

Diagnostik Model

Eksploratif

Secara eksploratif analisis sisaan dapat dilihat menggunakan QQ plot, residuals plot, ACF dan PACF plot.

sisaan <- model.s8$residuals
par(mfrow = c(1,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(sisaan, type="o", 
     ylab = "Sisaan", xlab = "Order", main = "Sisaan M1 vs Order")
abline(h = 0, col='red')

Acf <- ggAcf(sisaan, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()
Pacf <- ggPacf(sisaan, lag.max = 100) +
  scale_x_continuous(breaks = seq(0,100,10)) +
  theme_minimal()
grid.arrange(Acf, Pacf, ncol=2)

Normal Q-Q Plot

Berdasarkan output di atas, terlihat bahwa beberapa amatan cenderung tidak mengikuti garis qq-plot distribusi normal. Sehingga secara eksploratif, dapat disimpulkan bahwa sisaan tidak menyebar normal.

Residual vs Order

Berdasarkan hasil eksplorasi di atas, titik pada plot kebebasan sisaan mayoritas bergerak di sekitar titik nol. Namun, terdapat beberapa titik amatan yang terletak cukup jauh dari titik nol. Sehingga, belum dapat disimpulkan apakah terdapat autokorelasi atau tidak.

Plot ACF dan PACF

Berdasarkan hasil eksplorasi di atas, berdasarkan plot ACF maupun plot PACF terlihat terdapat garis vertikal di lag tertentu yang melebihi tinggi garis horizontal berwarna biru. Hal tersebut mengindikasikan adanya autokorelasi pada model.

Uji Formal

1) Sisaan Menyebar Normal

\[ \begin{aligned} H_0&: \text{ Sisaan mengikuti sebaran normal}\\ H_1&: \text{ Sisaan tidak mengikuti sebaran normal} \end{aligned} \]

ks.test(sisaan,"pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.5098, p-value < 2.2e-16
## alternative hypothesis: two-sided

Berdasarkan hasil uji, didapatkan nilai P−Value<0.05 yang berarti tolak H0. Jadi pada taraf nyata 5%, cukup bukti untuk menyatakan bahwa sisaan tidak mengikuti sebaran normal.

2) Sisaan saling bebas/tidak ada autokorelasi

Hipotesis yang digunakan adalah sebagai berikut: \[ \begin{aligned} H_0&: \text{ Tidak ada autokorelasi}\\ H_1&: \text{ Ada autokorelasi} \end{aligned} \]

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

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

Berdasarkan hasil uji, didapatkan nilai P−Value>0.05 yang berarti tak tolak H0. Jadi pada taraf nyata 5%, tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada data.

3) Nilai tengah sisaan sama dengan nol

Hipotesis yang digunakan adalah sebagai berikut: \[ \begin{aligned} H_0&: \text{ Nilai tengah sisaan bernilai nol}\\ H_1&: \text{ Nilai tengah sisaan tak bernilai nol} \end{aligned} \]

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

## 
##  One Sample t-test
## 
## data:  sisaan
## t = -0.85163, df = 101, p-value = 0.3964
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -27519.91  10988.09
## sample estimates:
## mean of x 
## -8265.913

Berdasarkan hasil uji, didapatkan nilai P−Value>0.05 yang berarti tak tolak H0. Jadi pada taraf nyata 5%, belum cukup bukti untuk menyatakan bahwa nilai tengah sisaan tidak bernilai nol.

4. Sisaan Homogen

Secara uji formal dengan menggunakan ARCH test. \[ \begin{aligned} H_0&: \text{ Sisaan homogen}\\ H_1&: \text{Sisaan tidak homogen} \end{aligned} \]

ArchTest(sisaan)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan
## Chi-squared = 19.221, df = 12, p-value = 0.08332

Berdasarkan ARCH test diperoleh p-value < α (0.05) maka tolak \(H_0\). Artinya, cukup bukti untuk menyatakan bahwa sisaan tidak homogen pada taraf nyata 5%.

Jadi berdasarkan uji formal yang telah dilakukan, diperoleh hasil bahwa sisaan tidak menyebar normal, tidak ada autokorelasi, nilai tengah sisaan bernilai nol, dan sisaan tidak homogen.

Forecasting

ramalan <- forecast::forecast(training.ts,model=model.s8,h=14)
ramalan

##          Point Forecast       Lo 80    Hi 80      Lo 95     Hi 95
## Jul 2021      122472.63   -6736.188 251681.5  -75135.21  320080.5
## Aug 2021      121831.78  -82203.561 325867.1 -190213.35  433876.9
## Sep 2021      115095.43 -147397.906 377588.8 -286353.49  516544.3
## Oct 2021      113953.04 -197582.376 425488.5 -362499.27  590405.4
## Nov 2021      108949.17 -245681.248 463579.6 -433411.25  651309.6
## Dec 2021      113498.99 -280148.288 507146.3 -488532.57  715530.6
## Jan 2022       98137.43 -331543.406 527818.3 -559002.70  755277.6
## Feb 2022       88640.55 -374783.451 552064.5 -620105.30  797386.4
## Mar 2022       91306.10 -404038.263 586650.5 -666257.73  848869.9
## Apr 2022       86512.88 -439261.000 612286.8 -717588.88  890614.6
## May 2022       93527.47 -461430.334 648485.3 -755207.25  942262.2
## Jun 2022       86865.20 -496218.002 669948.4 -804883.60  978614.0
## Jul 2022       83065.24 -542056.727 708187.2 -872976.31 1039106.8
## Aug 2022       80564.67 -587864.901 748994.2 -941710.14 1102839.5

data.ramalan <- ramalan$mean
plot(ramalan)
lines(wisata.ts,lty=1,col="black")
points(wisata.ts,cex=0.5)

perbandingan<-matrix(data=c(testing.ts, data.ramalan), nrow = 14, ncol = 2)
colnames(perbandingan)<-c("Aktual","Hasil Forecast")
perbandingan

##       Aktual Hasil Forecast
##  [1,] 127249      122472.63
##  [2,] 118533      121831.78
##  [3,] 120100      115095.43
##  [4,] 146137      113953.04
##  [5,] 150577      108949.17
##  [6,] 163619      113498.99
##  [7,] 143578       98137.43
##  [8,]  18455       88640.55
##  [9,]  40790       91306.10
## [10,] 111057       86512.88
## [11,] 212332       93527.47
## [12,] 345438       86865.20
## [13,] 476970       83065.24
## [14,] 510246       80564.67

Referensi

Hanke JE, Wichern DW. 2005. Business Forecasting Eight Edition. New Jersey (US): Pearson Prentice Hall.

Herawati S. 2016. Peramalan kunjungan wisatawan mancanegara menggunakan generalized regression neural networks. Jurnal Infotel. 8(1).

Kementerian Pariwisata dan Ekonomi Kreatif [Kemenparekraf]. 2020. Outlook Pariwisata & Ekonomi Kreatif Indonesia. Jakarta: Kemenparekraf.

Kementerian Pariwisata dan Ekonomi Kreatif/ Badan Pariwisata dan Ekonomi Kreatif [Kemenparekraf/Baparekraf]. 2020. Rencana Strategis Kemenparekraf/Baparekraf 2020-2024. Jakarta: Kemenparekraf/Baparekraf.

Makridakis, S., S. Wheelwright, R. Hyndman, and Y. Chang. 1998. Forecasting Methods and Applications. 3rd ed. New York: John Wiley & Sons

Nkoro Emeka, Uko Aham Kelvin . 2016. Autoregressive Distributed Lag (ARDL) cointegration technique:application and interpretation. Journal of Statistical and Econometric Methods. 5:63-91

Tinungki GM. 2018. The analysis of partial autocorrelation function in predicting maximum wind speed. The Electrochemical Society. 1-12. doi: 10.1088/1755-1315/235/1/012097

Wei WWS. 2006. Time Series Analysis. Boston (US): Department of Statistics Temple University.