Penerapan Model ARIMA pada Jumlah Wisatawan Mancanegara di Malang tahun 2018-2020

Rointan Margaretta Simamora

Jumat, 20 Mei 2022

1 PENDAHULUAN
2 SOURCE CODE
3 HASIL DAN PEMBAHASAN
4 DAFTAR PUSTAKA

Library:

> # install.packages("knitr")
> # install.packages("rmarkdown")
> # install.packages("prettydoc")
> # install.packages("equatiomatic")

1 PENDAHULUAN

1.1 Latar Belakang

Indonesia merupakan salah satu negara yang terkenal akan kekayaan budaya dan keindahan alamnya. Karena itu, Indonesia mempunyai potensi wisata yang menjanjikan. Pengembangan sektor pariwisata dapat meningkatkan sumber pendapatan negara maupun daerah, penerimaan devisa negara, dan penciptaan lapangan kerja. Sektor pariwisata memberi kontribusi dalam peningkatan devisa negara dari Rp175,71 triliun pada 2015 menjadi Rp229,50 triliun pada tahun 2018.

Namun, terjadi penurunan jumlah wisatawan yang signifikan pada tahun 2020 sebagai akibat dari pandemi COVID-19. Di masa pandemi, Indonesia memberlakukan beragam kebijakan untuk meminimalisasikan angka penyebaran COVID-19, seperti penutupan akses bagi turis-turis mancanegara dan penutupan objek wisata. Terjadi penurunan sebesar 75,03% dari tahun 2019.

Kota Malang merupakan salah satu kota tujuan wisata di Jawa Timur. Kota Malang memiliki berbagai objek wisata mulai dari wisata bersejarah seperti candi dan museum hingga wisata hiburan. Contoh destinasi wisata di Kota Malang adalah Jatim Park, Batu Night Spectacular, Malang Night Paradise, Kampung Warna-Warni Jodipan, Pantai Tiga Warna, Museum Angkut, dan sebagainya. Pada tahun 2017, jumlah wisatawan mancanegara di Malang sebanyak 12.546 orang dan wisatawan lokal sebanyak 4,3 juta orang. Kemudian meningkat pada tahun 2018 menjadi 15.034 wisatawan mancanegara dan 4,8 juta wisatawan nusantara. Namun, akibat pandemi COVID-19, kunjungan wisata Kota Malang menurun sampai 66,8%. Penurunan jumlah wisatawan ini berakibat pada penurunan pendapatan daerah dari sektor wisata.

1.2 Statistika Deskriptif

Statistika deskriptif adalah serangkaian teknik pengumpulan, penyajian, dan peringkasan sebagian atau seluruh data tanpa adanya pengambilan kesimpulan. Statistika deskriptif dapat digunakan untuk menggambarkan dan menganalisis data dengan menggunakan sedikitnya satu statistik contoh dengan membangun grafik dan tabel dengan membandingkan hasil data yang lain.

1.2.1 Penyajian Data

Diagram Dahan Daun (Steam and Leaf)
Apabila selisih dari nilai maksimum dengan nilai minimum data tidak terlalu besar, data dapat disajikan dalam diagram dahan daun (steam and leaf diagram). Pada prinsipnya, item data hanya dibedakan menjadi tangkai dan daun. Diagram ini dapat menunjukkan penyebaran dan pemusatan data dengan menggunakan data yang sebenarnya (bukan interval).
Histogram
Data disajikan dalam diagram yang berbentuk pesegi panjang. Sumbu-y histogram menggambarkan frekuensi data, sedangkan sumbu-x menggambarkan bilangan-bilangan data yang dinyatakan dalam interval kelas.
Diagram Kotak Garis (Box Plot)
Data juga dapat disajikan dalam bentuk diagram kotak garis dengan menentukan pagar dalam dan pagar luar data. Diagram ini digunakan untuk mengetahui penyebaran dan pemusatan data serta data pencilan.

1.2.2 Ukuran Pemusatan Data

Rata-Rata
Rata-rata merupakan pusat massa (centroid) sehingga simpangan kiri dan simpangan kanan sama besar. Rata-rata biasa dilambangkan dengan \(\mu\) dengan rumus : \[ \mu = \frac{\sum{X_i}}{n} \]
Median
Median adalah nilai tengah dari data terurut. Median disebut juga kuantil kedua. Jika banyaknya data ganjil, maka median merupakan nilai tengah dari data tersebut. Sedangkan, jika banyaknya data genap, maka median merupakan rata-rata dari 2 nilai tengah suatu data.
Modus
Modus adalah nilai yang paling sering muncul dalam data (frekuensi terbanyak).
Kuartil
Kuartil adalah nilai-nilai yang membagi data menjadi 4 bagian sama besar. Q₀ (dibaca kuartil 0) merupakan nilai minimum dari data, Q₁ (dibaca kuartil 1) merupakan nilai yang membagi data 25% data di kiri dan 75% data di kanan, Q₂ (dibaca kuartil 2) merupakan median, membagi data menjadi 2, Q₃ (dibaca kuartil 3) merupakan nilai yang membagi data 75% data di kiri dan 25% data di sebelah kanan, dan Q₄ (dibaca kuartil 4) merupakan nilai maksimum dari data.

1.2.3 Ukuran Penyebaran Data

Jangkauan (Range)
Jangkauan merupakan selisih dari nilai maksimum dengan nilai minimum.
Jangkauan Interkuartil (Interquartile Range)
IQR mrupakan selisih dari kuartil ketiga dengan kuartil pertama.
Deviasi Deviasi merupakan selisih dari data dengan rataannya. Ukuran keragaman dari deviasi adalah rataan deviasi
\[ \frac{∑(x -μ)}{n} \]
Ragam
Untuk menghilangan nilai +/-, maka deviasi dikuadratkan terlebih dahulu sebelum dirata-ratakan. Rumus :
\[ \sigma^2=\frac{\sum(x-\mu)^2}{n} \]
Simpangan Baku (Standard Deviation)
Ragam merupakan ukuran jarak kuadrat, sehingga untuk mendapatkan jarak yang sebenarnya adalah dengan mengakarkan ragam. Simpangan baku disimbolkan dengan s dan rumus : \[ \sigma=\sqrt{\sigma^2} \]

1.3 Analisis Deret Waktu

Data deret waktu adalah sekumpulan data observasi yang disajikan dalam urutan periode waktu, misalnya bulanan, tahunan, dan sebagainya. Tujuannya adalah menemukan pola data dan melakukan peramalan. Data deret waktu dengan n pengamatan dapat dinyatakan sebagai \(X_1, X_2, ..., X_n\). Dasar pemikiran deret waktu adalah pengamatan sekarang (X_t) dipengaruhi oleh pengamatan sebelumnya (X_t-k) yang dipisahkan oleh jarak waktu k (lag k).

1.3.1 Stasioneritas

Syarat utama penggunaan model ARIMA adalah stasioneritas terhadap rata-rata dan ragam. Data deret waktu dikatakan stasioner jika rata-rata dan ragamnya konstan, tidak dipengaruhi oleh waktu. Stasioneritas data deret waktu dapat dilihat secara visual dengan plot datanya. Jika plot memiliki tren (menaik atau menurun), maka data tidak stasioner terhadap rata-rata. Jika range data berbeda-beda (lebar plot berbeda), maka data tidak stasioner terhadap ragam.

Selain itu, dapat dilakukan uji unit root atau uji Augmented Dickey-Fuller (ADF) untuk menguji kestasioneran data terhadap rata-rata. Data dikatakan stasioner terhadap rata-rata jika \(p-value<\alpha\). Jika data tidak stasioner terhadap rata-rata, maka dilakukan differencing terhadap data deret waktu. Transformasi Box Cox dilakukan untuk mengecek kestasioneran data terhadap ragam. Jika data tidak stasioner terhadap ragam, maka dilakukan transformasi seperti berikut.
\[ x_{trans}=\frac{X_t^\lambda-1}{\lambda} \]

Data stasioner terhadap ragam jika nilai \(\lambda\) mendekati 1.

1.3.2 Klasifikasi Model ARIMA

Model Autoregressive (AR) Model AR berdasar pada asumsi data periode sekarang dipengaruhi oleh data periode sebelumnya. Biasanya ditulis AR(p) dengan p adalah orde model. Bentuk umum : \[ Y_t=\phi_1 Y_{t-1}+ \phi_2 Y_{t-2}+...+\phi_p Y_{t-p}+a_t \]
Model Moving Average (MA)
Model MA berdasar pada asumsi data periode sekarang dipengaruhi oleh nilai sisaan periode sebelumnya. Biasanya ditulis MA(q) dengan q adalah orde model. Bentuk umum : \[ Y_t=a_t-\theta_1 a_{t-1}+ \theta_2 a_{t-2}+...+\theta_q a_{t-q} \]
Model Autoregressive Moving Average (ARMA)
Model ARMA merupakan kombinasi model AR dan MA yang berdasar pada asumsi data periode sekarang dipengaruhi oleh data dan nilai sisaan periode sebelumnya. Biasanya ditulis ARMA(p,q) dengan bentuk umum : \[ Y_t=\phi_1 Y_{t-1}+ \phi_2 Y_{t-2}+...+\phi_p Y_{t-p}+a_t-\theta_1 a_{t-1}+ \theta_2 a_{t-2}+...+\theta_q a_{t-q} \]
Model Autoregressive Integrated Moving Average (ARIMA)
Model ARIMA mirip seperti model ARMA, tetapi model tidak stasioner sehingga harus dilakukan differencing. Biasanya ditulis ARIMA(p,d,q) dengan d adalah orde differencing. Bentuk umum : \[ Y_t=(1+\phi_1) Y_{t-1}+ (\phi_2-\phi_1) Y_{t-2}+...+(\phi_p-\phi_{p-1}) Y_{t-p}-\phi_p Y_{t-p-1}+a_t-\theta_1 a_{t-1}+ \theta_2 a_{t-2}+...+\theta_q a_{t-q} \]

1.3.3 Metodologi Analisis Deret Waktu

Identifikasi Model :

Membentuk plot data deret waktu.
Mengecek sifat stasioneritas data terhadap rata-rata dan ragam
Kestasioneran dapat diketahui melalui transformasi Box Cox dan uji ADF.
Identifikasi model dilakukan dengan melihat plot ACF dan PACF data deret waktu yang sudah stasioner. Dari plot tersebut, dapat dibentuk model tentatif.

Pendugaan Model
Setelah menentukan model tentatif, dilakukan pendugaan parameter model. Model yang akan dipilih adalah model dengan semua parameter signifikan.
Diagnostik Sisaan
Selanjutnya, dilakukan diagnostik sisaan pada model yang sudah signifikan. Diagnostik sisaan berupa pengujian asumsi normalitas dan non autokorelasi (white noise). Uji normalitas dapat dilakukan dengan uji Lilliefors, sedangkan uji non autokorelasi dilakukan dengan uji L-jung Box.
Peramalan :

Pemilihan model terbaik dilakukan dengan melihat nilai Akaike’s Information Criterion (AIC). AIC digunakan untuk menentukan model optimum dari suatu data observasi. Model terbaik akan memiliki nilai AIC terkecil.
Setelah model terbaik diperoleh, dilakukan peramalan data di masa yang akan datang.

1.4 Data

Data yang digunakan adalah data jumlah wisatawan mancanegara di Malang tahun 2018-2020 yang diperoleh dari website Badan Pusat Statistika Malang (https://malangkota.bps.go.id).

2 SOURCE CODE

2.1 Library yang Dibutuhkan

> library(readxl)
> library(MASS)
> library(FitAR)
> library(tseries)
> library(forecast)
> library(lmtest)
> library(nortest)

2.2 Memanggil Data

> data<-read_excel("D:/Intan/KOMSTAT/PRAKTIKUM/Jumlah Wisatawan.xlsx")
> data
# A tibble: 36 x 2
   Bulan   Jlh
   <dbl> <dbl>
 1     1  1098
 2     2   859
 3     3   591
 4     4   464
 5     5   639
 6     6  1137
 7     7  1346
 8     8  1745
 9     9  2176
10    10  1772
# ... with 26 more rows

Function read_excel digunakan untuk membaca file Excel “Jumlah Wisatawan”. Argumen yang diisikan dalam function adalah lokasi file yang digunakan.

2.3 Plot Deret Waktu

> plot(data, xlab="Bulan",ylab="Jumlah Wisatawan",
+      main="Plot  Jumlah Wisatawan ",type="l")

Function plot digunakan untuk menampilkan plot data jumlah wisatawan mancanegara di Malang tahun 2018-2020. Argumen yang diisikan dalam function adalah sumber data, nama sumbu x dan y, judul plot, serta tipe plot, yaitu line. Argumen ini digunakan sebagai identitas dan komponen plot data.

2.4 Statistika Deskriptif

> datakun <-data$Jlh
> summary(datakun)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      3     473    1118    1115    1696    2788 
> var(datakun)
[1] 604442.1
> median(datakun)
[1] 1117.5
> max(datakun)
[1] 2788
> min(datakun)
[1] 3
> quantile(datakun, 0.25)
25% 
473 
> quantile(datakun, 0.75)
    75% 
1696.25

Function summary digunakan untuk menampilkan deskripsi data berupa rata-rata, median, nilai maksimum, nilai minimum, kuantil bawah, dan kuantil atas. Fungsi var digunakan untuk menampilkan ragam data. Argumen yang diisikan dalam function adalah data jumlah wisatawan (data pada kolom “Jlh”). Argumen ini digunakan sebagai data yang akan dilakukan analisis statistika deskriptif.

2.5 Boxplot

> boxplot(datakun)

Function boxplot digunakan untuk menampilkan boxplot data jumlah wisatawan mancanegara di Malang tahun 2018-2020. Argumen yang diisikan dalam function adalah data jumlah wisatawan.

2.6 Uji Stasioneritas Ragam

> b <- boxcox(datakun~1)

> lambda <- b$x
> lik <- b$y
> bc <- cbind(lambda, lik)
> sorted_bc <- bc[order(-lik),]
> head(sorted_bc, n=10)
         lambda       lik
 [1,] 0.5454545 -64.83151
 [2,] 0.5858586 -64.86621
 [3,] 0.5050505 -64.87304
 [4,] 0.6262626 -64.97143
 [5,] 0.4646465 -64.99679
 [6,] 0.6666667 -65.14233
 [7,] 0.4242424 -65.21014
 [8,] 0.7070707 -65.37461
 [9,] 0.3838384 -65.52103
[10,] 0.7474747 -65.66427

Function boxcox digunakan untuk menampilkan grafik Box Cox data jumlah wisatawan. Argumen yang diisikan dalam function adalah data jumlah wisatawan dengan \(\lambda\) mendekati 1. Function head digunakan untuk menampilkan 10 teratas nilai pada tabel “sorted_bc”, di mana data telah diurutkan dari nilai likelihood terbesar.

2.6.1 Transformasi Box Cox

> transdata=((datakun^0.545)-1)/0.545
> boxcox(transdata~1)

2.7 Uji Stasioneritas Rata-Rata

> adf.test(transdata)

    Augmented Dickey-Fuller Test

data:  transdata
Dickey-Fuller = -2.3689, Lag order = 3, p-value = 0.4296
alternative hypothesis: stationary
> adf.test(diff(transdata))

    Augmented Dickey-Fuller Test

data:  diff(transdata)
Dickey-Fuller = -3.0916, Lag order = 3, p-value = 0.1492
alternative hypothesis: stationary
> adf.test(diff(diff(transdata)))

    Augmented Dickey-Fuller Test

data:  diff(diff(transdata))
Dickey-Fuller = -4.2393, Lag order = 3, p-value = 0.01249
alternative hypothesis: stationary

Function adf.test digunakan untuk melakukan uji Augmented Dickey-Fuller (ADF). Argumen yang diisikan dalam function adalah data yang akan diuji.

2.8 Identifikasi Orde

> acf(diff(diff(transdata)))

> pacf(diff(diff(transdata)))

Function acf dan pacf digunakan untuk menampilkan plot ACF dan PACF data transformasi. Argumen yang diisikan dalam function adalah sumber data, yaitu data transformasi dengan differencing 2 kali.

2.9 Model Tentatif

2.9.1 Model 1 (ARIMA(1,2,1))

> m1 <- Arima(datakun, order = c(1,2,1))
> summary(m1)
Series: datakun 
ARIMA(1,2,1) 

Coefficients:
         ar1      ma1
      0.2576  -1.0000
s.e.  0.1718   0.0867

sigma^2 estimated as 253206:  log likelihood=-260.25
AIC=526.5   AICc=527.3   BIC=531.08

Training set error measures:
                   ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
Training set -7.78123 474.4172 347.8165 -468.8127 537.007 0.9491328 0.02639449
> coeftest(m1)

z test of coefficients:

     Estimate Std. Error  z value Pr(>|z|)    
ar1  0.257556   0.171793   1.4992   0.1338    
ma1 -0.999999   0.086684 -11.5361   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.9.2 Model 2 (ARIMA(1,1,1))

> m2 <- Arima(datakun, order = c(1,1,1))
> summary(m2)
Series: datakun 
ARIMA(1,1,1) 

Coefficients:
          ar1     ma1
      -0.0418  0.3078
s.e.   0.4204  0.3916

sigma^2 estimated as 242167:  log likelihood=-265.63
AIC=537.25   AICc=538.02   BIC=541.92

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set -23.8051 471.1541 347.3823 -501.7499 547.8695 0.9479481
                    ACF1
Training set 0.003876336
> coeftest(m2)

z test of coefficients:

     Estimate Std. Error z value Pr(>|z|)
ar1 -0.041797   0.420373 -0.0994   0.9208
ma1  0.307847   0.391597  0.7861   0.4318

2.9.3 Model 3 (ARIMA(1,0,1))

> m3 <- Arima(datakun, order = c(1,0,1))
> summary(m3)
Series: datakun 
ARIMA(1,0,1) with non-zero mean 

Coefficients:
         ar1     ma1       mean
      0.6852  0.3567  1045.4638
s.e.  0.1356  0.1668   298.4582

sigma^2 estimated as 211618:  log likelihood=-270.85
AIC=549.69   AICc=550.98   BIC=556.02

Training set error measures:
                    ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
Training set -1.952849 440.4355 330.0642 -658.0401 679.786 0.9006898 0.06670911
> coeftest(m3)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ar1          0.68523    0.13557  5.0545 4.316e-07 ***
ma1          0.35669    0.16683  2.1381 0.0325109 *  
intercept 1045.46384  298.45819  3.5029 0.0004603 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.9.4 Model 4 (ARIMA(1,2,0))

> m4 <- Arima(datakun, order = c(1,2,0))
> summary(m4)
Series: datakun 
ARIMA(1,2,0) 

Coefficients:
          ar1
      -0.2623
s.e.   0.1747

sigma^2 estimated as 368639:  log likelihood=-265.67
AIC=535.34   AICc=535.73   BIC=538.39

Training set error measures:
                    ME    RMSE      MAE       MPE     MAPE     MASE        ACF1
Training set -1.093636 581.308 410.1983 -561.6149 747.1015 1.119362 -0.07122611
> coeftest(m4)

z test of coefficients:

    Estimate Std. Error z value Pr(>|z|)
ar1 -0.26234    0.17469 -1.5018   0.1332

2.9.5 Model 5 (ARIMA(1,1,0))

> m5 <- Arima(datakun, order = c(1,1,0))
> summary(m5)
Series: datakun 
ARIMA(1,1,0) 

Coefficients:
         ar1
      0.2217
s.e.  0.1655

sigma^2 estimated as 238708:  log likelihood=-265.88
AIC=535.77   AICc=536.14   BIC=538.88

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE
Training set -21.98842 474.8118 354.9839 -490.4249 545.5829 0.9686916
                   ACF1
Training set 0.04696733
> coeftest(m5)

z test of coefficients:

    Estimate Std. Error z value Pr(>|z|)
ar1  0.22167    0.16551  1.3393   0.1805

2.9.6 Model 6 (ARIMA(1,0,0))

> m6 <- Arima(datakun, order = c(1,0,0))
> summary(m6)
Series: datakun 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1       mean
      0.7880  1037.3820
s.e.  0.0974   335.2888

sigma^2 estimated as 229492:  log likelihood=-272.72
AIC=551.45   AICc=552.2   BIC=556.2

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
Training set -1.974087 465.5559 354.5935 -591.1628 614.0418 0.9676262 0.2895128
> coeftest(m6)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ar1       7.8798e-01 9.7394e-02  8.0907 5.932e-16 ***
intercept 1.0374e+03 3.3529e+02  3.0940  0.001975 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.9.7 Model 7 (ARIMA(0,2,1))

> m7 <- Arima(datakun, order = c(0,2,1))
> summary(m7)
Series: datakun 
ARIMA(0,2,1) 

Coefficients:
          ma1
      -1.0000
s.e.   0.1032

sigma^2 estimated as 258214:  log likelihood=-261.36
AIC=526.72   AICc=527.11   BIC=529.77

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
Training set -11.52385 486.5147 347.7856 -418.7952 455.7602 0.9490484 0.2266461
> coeftest(m7)

z test of coefficients:

    Estimate Std. Error z value  Pr(>|z|)    
ma1 -0.99998    0.10317  -9.693 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.9.8 Model 8 (ARIMA(0,1,1))

> m8 <- Arima(datakun, order = c(0,1,1))
> summary(m8)
Series: datakun 
ARIMA(0,1,1) 

Coefficients:
         ma1
      0.2726
s.e.  0.1672

sigma^2 estimated as 235099:  log likelihood=-265.63
AIC=535.26   AICc=535.64   BIC=538.37

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE
Training set -23.54843 471.2091 348.0462 -502.6779 550.7821 0.9497597
                      ACF1
Training set -0.0006301334
> coeftest(m8)

z test of coefficients:

    Estimate Std. Error z value Pr(>|z|)
ma1  0.27258    0.16724  1.6298   0.1031

2.9.9 Model 9 (ARIMA(0,0,1))

> m9 <- Arima(datakun, order = c(0,0,1))
> summary(m9)
Series: datakun 
ARIMA(0,0,1) with non-zero mean 

Coefficients:
         ma1       mean
      0.7226  1093.4885
s.e.  0.0927   152.0131

sigma^2 estimated as 303453:  log likelihood=-277.64
AIC=561.27   AICc=562.02   BIC=566.02

Training set error measures:
                     ME     RMSE      MAE       MPE     MAPE     MASE      ACF1
Training set -0.3424619 535.3454 429.0219 -950.8358 974.1028 1.170729 0.4097904
> coeftest(m9)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ma1          0.72257    0.09267  7.7971 6.333e-15 ***
intercept 1093.48850  152.01305  7.1934 6.320e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Function Arima digunakan untuk melakukan pemodelan ARIMA pada data. Argumen yang diisikan dalam function adalah data jumlah wisatawan.
Function summary digunakan untuk menampilkan nilai duga parameter dan nilai kebaikan model lain. Argumen yang diisikan dalam function adalah model ARIMA yang digunakan.
Function coeftest digunakan untuk menampilkan uji signifikansi parameter model. Argumen yang diisikan dalam function adalah model ARIMA yang digunakan.

2.10 Diagnostik Sisaan

2.10.1 Asumsi Non Autokorelasi

> checkresiduals(m3)


    Ljung-Box test

data:  Residuals from ARIMA(1,0,1) with non-zero mean
Q* = 1.2713, df = 4, p-value = 0.8662

Model df: 3.   Total lags used: 7
> checkresiduals(m6)


    Ljung-Box test

data:  Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 4.8936, df = 5, p-value = 0.429

Model df: 2.   Total lags used: 7
> checkresiduals(m7)


    Ljung-Box test

data:  Residuals from ARIMA(0,2,1)
Q* = 4.9139, df = 6, p-value = 0.5549

Model df: 1.   Total lags used: 7
> checkresiduals(m9)


    Ljung-Box test

data:  Residuals from ARIMA(0,0,1) with non-zero mean
Q* = 18.285, df = 5, p-value = 0.002609

Model df: 2.   Total lags used: 7

Function checkresiduals digunakan untuk melakukan uji L-jung Box pada sisaan dan menampilkan plot sisaan. Argumen yang diisikan dalam function adalah model ARIMA yang akan diuji asumsi non autokorelasi.

2.10.2 Asumsi Normalitas

> lillie.test(m3$residuals)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  m3$residuals
D = 0.12888, p-value = 0.1362
> lillie.test(m6$residuals)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  m6$residuals
D = 0.12225, p-value = 0.1895
> lillie.test(m7$residuals)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  m7$residuals
D = 0.11676, p-value = 0.2456

Function lillie.test digunakan untuk melakukan uji Lilliefors. Argumen yang diisikan dalam function adalah sisaan model ARIMA yang akan diuji asumsi normalitas.

2.11 Pemilihan Model Terbaik

> AIC(m3)
[1] 549.69
> AIC(m6)
[1] 551.4463
> AIC(m7)
[1] 526.72

Function AIC digunakan untuk menampilkan nilai AIC model. Argumen yang diisikan dalam function adalah model ARIMA yang akan digunakan.

2.12 Peramalan

> fmo <- forecast(m7, h=12)
> fmo
   Point Forecast      Lo 80     Hi 80      Lo 95    Hi 95
37      201.08581  -459.3625  861.5342  -808.9828 1211.154
38      176.17162  -770.7270 1123.0702 -1271.9848 1624.328
39      151.25742 -1024.0190 1326.5339 -1646.1728 1948.688
40      126.34323 -1248.4896 1501.1761 -1976.2822 2228.969
41      101.42904 -1455.2626 1658.1207 -2279.3254 2482.184
42       76.51485 -1649.9398 1802.9694 -2563.8696 2716.899
43       51.60065 -1835.7877 1938.9890 -2834.9108 2938.112
44       26.68646 -2014.8956 2068.2685 -3095.6438 3149.017
45        1.77227 -2188.6871 2192.2316 -3348.2463 3351.791
46      -23.14192 -2358.1795 2311.8956 -3594.2737 3547.990
47      -48.05611 -2524.1258 2428.0136 -3834.8780 3738.766
48      -72.97031 -2687.1002 2541.1596 -4070.9370 3924.996
> plot(fmo,xlab="Bulan",ylab="Jumlah Wisatawan",
+      main="Peramalan Jumlah Wisatawan")

Function forecast digunakan untuk melakukan peramalan data. Argumen yang diisikan dalam function adalah model ARIMA dan banyak periode yang akan diramal. Argumen ini digunakan untuk melakukan peramlan data model “m7” dengan periode 12 bulan ke depan.

3 HASIL DAN PEMBAHASAN

3.1 Plot Deret Waktu

Plot data jumlah wisatawan mancanegara di Malang tahun 2018-2020 menunjukkan tren sedikit menurut, data tidak stasioner terhadap rata-rata. Dari plot ini juga dapat diketahui bahwa data tidak stasioner terhadap ragam karena range yang berbeda-beda.

3.2 Statistika Deskriptif

Berdasarkan data yang digunakan, diperoleh nilai rataan sebesar 1114,889. Artinya, sebagian besar jumlah wisatawan mancanegara di Malang tahun 2018-2020 adalah 1115 orang.

\[ \mu = \frac{\sum{X_i}}{n}=1114.889 \] Diperoleh nilai ragam sebesar 60442,1. Artinya, keragaman data jumlah wisatawan mancanegara di Malang tahun 2018-2020 adalah 60442,1. \[ \sigma^2=\sqrt{\frac{\sum({X_i-\mu})^2}{n}} =604442.1 \] Diperoleh nilai median sebesar 1117,5. Artinya, terdapat 18 sampai 19 bulan jumlah wisatawan mancanegara di Malang tahun 2018-2020 yang kurang dari 1118 orang. \[ median=X_\frac{2(n+1)}{4} =1117.5 \] Diperoleh nilai maksimum dan minimum masing-masing sebesar 2788 dan 3. Artinya, jumlah wisatawan mancanegara di Malang tahun 2018-2020 paling banyak 2788 orang dan paling kecil 3 orang. \[ max=2788 \\ min=3 \] Diperoleh nilai kuartil bawah sebesar 473. Artinya, terdapat 8 sampai 9 bulan jumlah wisatawan mancanegara di Malang tahun 2018-2020 yang kurang dari 473 orang. \[ Q_1=\frac{X_\frac{n-1}{4}+X_{\frac{n+3}{4}}}{2} =473 \]

Diperoleh nilai kuartil atas sebesar 1696,25. Artinya, terdapat 27 sampai 28 bulan di mana jumlah wisatawan mancanegara di Malang tahun 2018-2020 yang kurang dari 1697 orang. \[ Q_3=\frac{X_\frac{3n+1}{4}+X_{\frac{3n+5}{4}}}{2} =1696.25 \]

3.3 Boxplot

Dari Boxplot, dapat diketahui bahwa data jumlah wisatawan mancanegara di Malang tahun 2018-2020 cenderung menceng ke kiri yang menunjukkan data tidak menyebar normal.

3.4 Analisis Deret Waktu

3.4.1 Stasioneritas Ragam

Dari uji Box Cox, diperoleh \(\lambda=0,545\). Karena nilai \(\lambda<1\), maka dilakukan transformasi Box Cox seperti berikut \[ x_{trans}=\frac{X_t^\lambda-1}{\lambda} \] Dari data hasil transformasi, dilakukan uji Box Cox kembali dan diperoleh \(\lambda\) yang sudah mendekati 1.

3.4.2 Stasioneritas Rata-Rata

Setelah, data stasioner terhadap ragam, dilakukan uji Augmented Dickey-Fuller (ADF) untuk menguji kestasioneran data terhadap rata-rata.

Hipotesis :
H₀ = data tidak stasioner
H₁ = data stasioner

Dari uji ADF, diperoleh \(p-value=0.4296>\alpha(0.05)\), berarti data belum stasioner terhadap rata-rata. Dilakukan differencing 2 kali dan diperoleh \(p-value=0.01249<\alpha(0.05)\), data sudah stasioner.

3.5 Pemodelan

3.5.1 Penentuan Orde ARIMA

Dari plot ACF dan PACF, dapat diketahui bahwa lag signifikan pada lag 1 sehingga orde maksimal untuk p dan q adalah 1.

3.5.2 Penentuan Model Tentatif

Terdapat 9 model tentatif, yaitu :
1. ARIMA(1,2,1)
2. ARIMA(1,1,1)
3. ARIMA(1,0,1)
4. ARIMA(1,2,0)
5. ARIMA(1,1,0)
6. ARIMA(1,0,0)
7. ARIMA(0,2,1)
8. ARIMA(0,1,1)
9. ARIMA(0,0,1)

3.5.3 Uji Signifikansi Parameter Model

Dari uji signifikansi (uji z), model yang semua parameternya signifikan adalah model :
3. ARIMA(1,0,1) dengan \(\hat\phi=0,685\) dan \(\hat\theta=0,357\)
6. ARIMA(1,0,0) dengan \(\hat\phi=0,788\)
7. ARIMA(0,2,1) dengan \(\hat\theta=-0,999\)
9. ARIMA(0,0,1) dengan \(\hat\theta=0,723\)

3.6 Diagnostik Sisaan

3.6.1 Uji Asumsi Non Autokorelasi

Hipotesis :
H₀ = tidak terdapat autokorelasi
H₁ = terdapat autokorelasi

Dengan uji Ljung-Box, dilakukan pengujian terhadap keempat model dan diperoleh :
* \(p-value\) untuk model ARIMA(1,0,1), ARIMA(1,0,0), dan ARIMA(0,2,1) lebih dari \(\alpha=0.05\). Artinya, tidak terdapat autokorelasi pada sisaan ketiga model.
* \(p-value\) untuk model ARIMA(0,0,1) kurang dari \(\alpha=0.05\). Artinya, terdapat autokorelasi pada sisaan model.

3.6.2 Uji Asumsi Normalitas

Hipotesis :
H₀ = sisaan berdistribusi normal
H₁ = sisaan tidak berdistribusi normal

Dilakukan uji Lilliefors terhadap model ARIMA(1,0,1), ARIMA(1,0,0), dan ARIMA(0,2,1) yang tidak memiliki autokorelasi sisaan. Diperoleh \(p-value\) untuk ketiga model lebih dari \(\alpha=0.05\). Artinya, sisaan model berdistribusi normal.

3.7 Pemilihan Model Terbaik

Pemilihan model terbaik dilakukan dengan memeriksa AIC ketiga model. Model terbaik adalah model dengan nilai AIC terkecil, yaitu model ARIMA(0,2,1) atau IMA(2,1) (AIC=526,72).

\[ X_t=2X_{t-1}-X_{t-2}+a_t-\theta a_{t-1} \\ X_t=2X_{t-1}-X_{t-2}+a_t+0.999 a_{t-1} \]

3.8 Peramalan

Dari model IMA(2,1), dilakukan peramalan 12 bulan ke depan dan diperoleh hasil seperti berikut.

> fore <- fmo$mean
> as.data.frame(fore)
           x
1  201.08581
2  176.17162
3  151.25742
4  126.34323
5  101.42904
6   76.51485
7   51.60065
8   26.68646
9    1.77227
10 -23.14192
11 -48.05611
12 -72.97031

Data ramalan menunjukkan tren menurun. Artinya, jumlah wisatawan mancanegara di Malang 12 bulan ke depan akan mengalami penurunan.

4 DAFTAR PUSTAKA

Carep-04. (2021, 1 Februari). Kunjungan Wisata Kota Malang 2020 Jeblok Sampai 66,8 Persen. Diakses pada 20 Mei 2022, dari https://kabarmalang.com/18013/kunjungan-wisata-kota-malang-2020-jeblok-sampai-668-persen

egsa.geo.ugm.ac.id. (2021, 11 Februari). Pariwisata Indonesia di Tengah Pandemi. Diakses pada 20 Mei 2022, dari https://egsa.geo.ugm.ac.id/2021/02/11/pariwisata-indonesia-di-tengah-pandemi/

Hikmah, Al. (2017). Peramalan Deret Waktu dengan Menggunakan Autoregressive (AR), Jaringan Syaraf Tiruan Radial Basic Function (RBF), dan Hybrid AR-RBF pada Inflasi Indonesia. Skripsi. Universitas Negeri Semarang, Semarang

Kustituanto, Bambang dan Rudy Badrudin, 1994. Statistika 1 (Deskriptif). Jakarta: Gunadarma

Ludin, Jimmy. (2020). Bahan Ajar : Analisis Deret Berkala. Pusat Pendidikan dan Pelatihan Badan Pusat Statistika. Tulisan pada https://pusdiklat.bps.go.id/