Pemodelan Data Deret Waktu Kunjungan Wisatawan Mancanegara ke Indonesia

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 Permasalahan

Pada tahun 2020, Kementerian Pariwisata dan Ekonomi Kreatif/Badan 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.

Pra-Processing Data

Memanggil Package di R

Untuk melakukan proses pengecekan kestasioneran data, diperlukan beberapa library untuk menunjang proses analisis data sebagai berikut.

library(googlesheets4)
library(tseries)
library(forecast)
library(TTR)
library(TSA)
library(imputeTS)
library(ggplot2)
library(FinTS)

Impor Data

Data yang digunakan hanyalah data Wisatawan dengan syntax berikut:

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")
wisata

## # A tibble: 116 × 2
##    Bulan      Wisatawan
##    <date>         <dbl>
##  1 2013-01-31    614328
##  2 2013-02-28    678415
##  3 2013-03-31    725316
##  4 2013-04-30    646117
##  5 2013-05-31    700708
##  6 2013-06-30    789594
##  7 2013-07-31    717784
##  8 2013-08-31    771009
##  9 2013-09-30    770878
## 10 2013-10-31    719903
## # … with 106 more rows

Analisis Eksplorasi Data

wisata.ts = ts(wisata$Wisatawan, frequency = 12, start= 2013)
plot(wisata.ts, main = "Kunjungan Wisatawan Mancanegara Periode 2013-2022", 
     xlab = "Tahun", ylab="Wisatawan")
points(wisata.ts)

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]

training.ts <- ts(train, start=2013, frequency = 12)
testing.ts <- ts(test, start=2021, frequency = 12)

plot(training.ts, col="blue", 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(testing.ts, col="red", lwd=2)

legend("bottomleft",c("Data Training", "Data Testing"), lty=1, col=c("blue", "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.

Pengecekan Kestasioneran Data

Stasioner Terhadap Rataan

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.

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 plot ACF dan pada lag ke-4 pada plot PACF. Maka dari itu, dapat dideteksi dua potensi model data deret waktu, yaitu adalah ARIMA(0,2,1) dan ARIMA(4,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(1,2,1), ARIMA(2,2,1), ARIMA(4,2,4), dan ARIMA(3,2,4).

Kesimpulan

Berdasarkan eksplorasi dan pengujian statistik di atas, didapati beberapa potensi model ARIMA pada amatan jumlah kunjungan wisatawan mancanegara ke Indonesia dengan rincian sebagai berikut,

ARIMA(0,2,1)
ARIMA(4,2,0)
ARIMA(1,2,1)
ARIMA(2,2,1)
ARIMA(4,2,4)
ARIMA(3,2,4)

Selanjutnya, perlu dilakukan pengujian signifikansi parameter pada tiap model potensial serta asumsi-asumsi yang perlu dipenuhi untuk mendapatkan model yang terbaik.

Daftar Pustaka

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