Analisis deret waktu musiman adalah proses memahami, memodelkan, dan meramalkan data yang memiliki pola musiman atau siklus teratur. Data deret waktu musiman adalah data yang menunjukkan fluktuasi periodik atau pola teratur yang terjadi dalam interval waktu tertentu, seperti harian, bulanan, kuartalan, atau tahunan. Analisis deret waktu musiman penting dalam berbagai bidang, termasuk ekonomi, keuangan, meteorologi, dan produksi.

Berikut adalah beberapa langkah pendahuluan yang umum dilakukan dalam analisis deret waktu musiman:

Pengamatan Data: Langkah pertama adalah mengamati data deret waktu untuk mengidentifikasi pola musiman. Ini melibatkan visualisasi data menggunakan grafik deret waktu dan histogram, serta pengamatan terhadap fluktuasi yang berulang dalam interval waktu tertentu.

Differencing: Jika data tidak stasioner, langkah berikutnya adalah menerapkan differencing untuk membuat data menjadi stasioner. Differencing dapat dilakukan untuk menghilangkan tren dan/atau musim.

Identifikasi Model: Setelah data menjadi stasioner, langkah selanjutnya adalah mengidentifikasi model yang sesuai untuk meramalkan deret waktu musiman. Ini melibatkan pemilihan model ARIMA (Autoregressive Integrated Moving Average) atau model musiman seperti SARIMA (Seasonal ARIMA) berdasarkan analisis ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function).

Estimasi Parameter: Setelah model diidentifikasi, parameter model harus diestimasi menggunakan metode seperti Maksimum Likelihood (ML) atau Metode Kuadrat Terkecil (OLS).

Uji Model: Penting untuk menguji kecocokan model terhadap data, seperti uji kecocokan, uji heteroskedastisitas, dan uji autokorelasi untuk memastikan bahwa model yang dibangun cocok dengan data.

Validasi dan Peramalan: Setelah model dinyatakan sesuai, itu dapat digunakan untuk membuat peramalan untuk periode yang akan datang. Validasi peramalan dapat dilakukan dengan membandingkan hasil peramalan dengan data aktual dan mengukur tingkat akurasi model.

Evaluasi: Akhirnya, model dapat dievaluasi berdasarkan kinerja peramalan, kesesuaian dengan data, dan kemampuan model untuk menangkap pola musiman yang mungkin terjadi dalam data.

Analisis deret waktu musiman memungkinkan pemahaman yang lebih baik tentang pola yang terjadi dalam data dan membantu dalam pengambilan keputusan yang lebih baik dalam merencanakan, memprediksi, dan mengelola proses atau fenomena yang memiliki sifat musiman.

library(tseries)
## Warning: package 'tseries' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(NormalityAssessment)
## Warning: package 'NormalityAssessment' was built under R version 4.3.2
library(nortest)
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2

library(tseries): Ini adalah perintah untuk memuat paket “tseries” di R. Paket ini berisi berbagai fungsi untuk analisis deret waktu. library(forecast): Ini adalah perintah untuk memuat paket “forecast” di R. Paket ini berisi fungsi-fungsi yang digunakan untuk meramalkan data deret waktu. library(lmtest): Ini adalah perintah untuk memuat paket “lmtest” di R. Paket ini menyediakan berbagai tes untuk memeriksa asumsi-asumsi dalam model regresi linear. library(NormalityAssessment): Ini adalah perintah untuk memuat paket “NormalityAssessment” di R. Paket ini menyediakan fungsi-fungsi untuk mengevaluasi normalitas distribusi data. library(nortest): Ini adalah perintah untuk memuat paket “nortest” di R. Paket ini berisi berbagai tes untuk menguji normalitas distribusi data. library(readxl): Ini adalah perintah untuk memuat paket “readxl” di R. Paket ini menyediakan fungsi-fungsi untuk membaca file Excel (.xls dan .xlsx).

data <- read.csv("C:/Users/dayah/Downloads/penumpang.csv", sep=";")
data
##       t  zt
## 1     1 112
## 2     2 118
## 3     3 132
## 4     4 129
## 5     5 121
## 6     6 135
## 7     7 148
## 8     8 148
## 9     9 136
## 10   10 119
## 11   11 104
## 12   12 118
## 13   13 115
## 14   14 126
## 15   15 141
## 16   16 135
## 17   17 125
## 18   18 149
## 19   19 170
## 20   20 170
## 21   21 158
## 22   22 133
## 23   23 114
## 24   24 140
## 25   25 145
## 26   26 150
## 27   27 178
## 28   28 163
## 29   29 172
## 30   30 178
## 31   31 199
## 32   32 199
## 33   33 184
## 34   34 162
## 35   35 146
## 36   36 166
## 37   37 171
## 38   38 180
## 39   39 193
## 40   40 181
## 41   41 183
## 42   42 218
## 43   43 230
## 44   44 242
## 45   45 209
## 46   46 191
## 47   47 172
## 48   48 194
## 49   49 196
## 50   50 196
## 51   51 236
## 52   52 235
## 53   53 229
## 54   54 243
## 55   55 264
## 56   56 272
## 57   57 237
## 58   58 211
## 59   59 180
## 60   60 201
## 61   61 204
## 62   62 188
## 63   63 235
## 64   64 227
## 65   65 234
## 66   66 264
## 67   67 302
## 68   68 293
## 69   69 259
## 70   70 229
## 71   71 203
## 72   72 229
## 73   73 242
## 74   74 233
## 75   75 267
## 76   76 269
## 77   77 270
## 78   78 315
## 79   79 364
## 80   80 347
## 81   81 312
## 82   82 274
## 83   83 237
## 84   84 278
## 85   85 284
## 86   86 277
## 87   87 317
## 88   88 313
## 89   89 318
## 90   90 374
## 91   91 413
## 92   92 405
## 93   93 355
## 94   94 306
## 95   95 271
## 96   96 306
## 97   97 315
## 98   98 301
## 99   99 356
## 100 100 348
## 101 101 355
## 102 102 422
## 103 103 465
## 104 104 467
## 105 105 404
## 106 106 347
## 107 107 305
## 108 108 336
## 109 109 340
## 110 110 318
## 111 111 362
## 112 112 348
## 113 113 363
## 114 114 435
## 115 115 491
## 116 116 505
## 117 117 404
## 118 118 359
## 119 119 310
## 120 120 337
## 121 121 360
## 122 122 342
## 123 123 406
## 124 124 396
## 125 125 420
## 126 126 472
## 127 127 548
## 128 128 559
## 129 129 463
## 130 130 407
## 131 131 362
## 132 132 405
## 133 133 417
## 134 134 391
## 135 135 419
## 136 136 461
## 137 137 472
## 138 138 535
## 139 139 622
## 140 140 606
## 141 141 508
## 142 142 461
## 143 143 390
## 144 144 432
head(data)
##   t  zt
## 1 1 112
## 2 2 118
## 3 3 132
## 4 4 129
## 5 5 121
## 6 6 135
tail(data)
##       t  zt
## 139 139 622
## 140 140 606
## 141 141 508
## 142 142 461
## 143 143 390
## 144 144 432

head ini berfungsi untuk melihat data awal dan untuk fungsi tail untuk melihat akhir dari suatu data

#plot time series
#melabel data zt dengan"Penumpang"
penumpang<-data$zt
Bulan<-data$t

ztmengambil nilai−nilai dari kolom”zt”dan menyimpannya dalam objek penumpang,sementara sintaks Bulan<−data t mengambil nilai-nilai dari kolom “t” dan menyimpannya dalam objek Bulan.

#Tahap identifikasi
ts.plot(penumpang)

acf(penumpang, lag.max = 80)

membuat plot deret waktu dari data yang disimpan dalam objek penumpang. Ini akan menampilkan grafik yang menunjukkan bagaimana data penumpang berubah sepanjang waktu. Plot ini bisa membantu untuk melihat pola-pola atau tren yang mungkin ada dalam data. Fungsi acf() digunakan untuk menghitung dan memplot fungsi korelasi otokorelasi (ACF) dari deret waktu. Dalam kasus ini, deret waktu yang digunakan adalah penumpang. Argumen lag.max = 144 menentukan jumlah maksimum lag yang akan dihitung untuk ACF. Dalam hal ini, maksimum lag yang dihitung adalah 144, yang mungkin sesuai dengan asumsi bulanan jika data waktu diukur dalam bulan. #Kedua perintah ini membantu dalam menganalisis pola dan korelasi dalam data deret waktu, yang penting untuk pemodelan dan peramalan. dan plotnya belum stasioner dalam rata-rata maupun dalam variansi karena terdapat kecenderungan naik tidak bisa membentuk garis horizontal. yang perlu dilakukan untuk lebih lanjut adalah menstasionerkan terlebih dahulu. `}

#plot autokorelasi dan partial
acf(penumpang, lag.max = 75)
par(mfrow=c(1,2))
acf(penumpang, lag.max = 75)
pacf(penumpang, lag.max = 75)

Fungsi acf() digunakan untuk menghitung dan memplot fungsi korelasi otokorelasi (ACF) dari deret waktu yang disimpan dalam objek penumpang. Dalam hal ini, argumen lag.max = 75 menentukan bahwa hanya korelasi hingga lag maksimum 75 yang akan dihitung dan ditampilkan. Ini membantu dalam melihat pola korelasi otokorelasi dalam deret waktu. par(mfrow=c(1,2)): Fungsi par() digunakan untuk mengatur parameter gambar di R. Argumen mfrow=c(1,2) mengatur layout gambar menjadi satu baris dan dua kolom. Ini berarti gambar-gambar berikutnya akan diatur dalam dua kotak secara horizontal. Setelah pengaturan parameter gambar, dua plot akan ditampilkan: Plot pertama adalah ACF (Autocorrelation Function) dari deret waktu penumpang, dengan lag maksimum 75. Plot kedua adalah PACF (Partial Autocorrelation Function) dari deret waktu yang sama, juga dengan lag maksimum 75. Ini membantu untuk memahami pola korelasi otokorelasi dan parsial otokorelasi dalam deret waktu, yang berguna dalam pemodelan dan peramalan.

#differencing 1 non musiman dengan beberapa cara
differencing1 <- diff(penumpang, differences=1)
differencing1
##   [1]    6   14   -3   -8   14   13    0  -12  -17  -15   14   -3   11   15   -6
##  [16]  -10   24   21    0  -12  -25  -19   26    5    5   28  -15    9    6   21
##  [31]    0  -15  -22  -16   20    5    9   13  -12    2   35   12   12  -33  -18
##  [46]  -19   22    2    0   40   -1   -6   14   21    8  -35  -26  -31   21    3
##  [61]  -16   47   -8    7   30   38   -9  -34  -30  -26   26   13   -9   34    2
##  [76]    1   45   49  -17  -35  -38  -37   41    6   -7   40   -4    5   56   39
##  [91]   -8  -50  -49  -35   35    9  -14   55   -8    7   67   43    2  -63  -57
## [106]  -42   31    4  -22   44  -14   15   72   56   14 -101  -45  -49   27   23
## [121]  -18   64  -10   24   52   76   11  -96  -56  -45   43   12  -26   28   42
## [136]   11   63   87  -16  -98  -47  -71   42

Differencing pertama menghasilkan deret waktu baru yang nilainya adalah selisih antara nilai deret waktu pada waktu tertentu dengan nilai deret waktu pada waktu sebelumnya.Jadi, differencing1 akan berisi deret waktu yang merupakan hasil dari differencing pertama pada penumpang

differencing2 <- diff(penumpang)
differencing2
##   [1]    6   14   -3   -8   14   13    0  -12  -17  -15   14   -3   11   15   -6
##  [16]  -10   24   21    0  -12  -25  -19   26    5    5   28  -15    9    6   21
##  [31]    0  -15  -22  -16   20    5    9   13  -12    2   35   12   12  -33  -18
##  [46]  -19   22    2    0   40   -1   -6   14   21    8  -35  -26  -31   21    3
##  [61]  -16   47   -8    7   30   38   -9  -34  -30  -26   26   13   -9   34    2
##  [76]    1   45   49  -17  -35  -38  -37   41    6   -7   40   -4    5   56   39
##  [91]   -8  -50  -49  -35   35    9  -14   55   -8    7   67   43    2  -63  -57
## [106]  -42   31    4  -22   44  -14   15   72   56   14 -101  -45  -49   27   23
## [121]  -18   64  -10   24   52   76   11  -96  -56  -45   43   12  -26   28   42
## [136]   11   63   87  -16  -98  -47  -71   42

Differencing2 akan berisi deret waktu yang merupakan hasil dari differencing pertama pada penumpang Ini dilakukan dengan mengurangkan setiap titik data dari titik data sebelumnya dalam deret waktu hasil differencing pertama. Hasil differencing kedua disimpan dalam objek baru yang disebut differencing2. Anda dapat menggunakan objek ini untuk analisis stasioneritas lebih lanjut atau untuk membangun model dan melakukan peramalan.

differencing3 <- diff(penumpang, lag=1)
differencing3
##   [1]    6   14   -3   -8   14   13    0  -12  -17  -15   14   -3   11   15   -6
##  [16]  -10   24   21    0  -12  -25  -19   26    5    5   28  -15    9    6   21
##  [31]    0  -15  -22  -16   20    5    9   13  -12    2   35   12   12  -33  -18
##  [46]  -19   22    2    0   40   -1   -6   14   21    8  -35  -26  -31   21    3
##  [61]  -16   47   -8    7   30   38   -9  -34  -30  -26   26   13   -9   34    2
##  [76]    1   45   49  -17  -35  -38  -37   41    6   -7   40   -4    5   56   39
##  [91]   -8  -50  -49  -35   35    9  -14   55   -8    7   67   43    2  -63  -57
## [106]  -42   31    4  -22   44  -14   15   72   56   14 -101  -45  -49   27   23
## [121]  -18   64  -10   24   52   76   11  -96  -56  -45   43   12  -26   28   42
## [136]   11   63   87  -16  -98  -47  -71   42
#plot time series, autokorelasi dan patial
ts.plot(differencing1)

#differencing 1 musiman 12 dari data yang sudah didiferencing
diffnonmus_mus12 <- diff(differencing1, lag = 12)
diffnonmus_mus12
##   [1]   5   1  -3  -2  10   8   0   0  -8  -4  12   8  -6  13  -9  19 -18   0
##  [19]   0  -3   3   3  -6   0   4 -15   3  -7  29  -9  12 -18   4  -3   2  -3
##  [37]  -9  27  11  -8 -21   9  -4  -2  -8 -12  -1   1 -16   7  -7  13  16  17
##  [55] -17   1  -4   5   5  10   7 -13  10  -6  15  11  -8  -1  -8 -11  15  -7
##  [73]   2   6  -6   4  11 -10   9 -15 -11   2  -6   3  -7  15  -4   2  11   4
##  [91]  10 -13  -8  -7  -4  -5  -8 -11  -6   8   5  13  12 -38  12  -7  -4  19
## [109]   4  20   4   9 -20  20  -3   5 -11   4  16 -11  -8 -36  52 -13  11  11
## [127] -27  -2   9 -26  -1

differencing1, lag = 12) digunakan untuk melakukan differencing pada deret waktu yang telah mengalami differencing pertama dan disimpan dalam objek differencing1, dengan menggunakan lag atau jarak 12. Ini berarti nilai pada titik waktu ke-t akan dikurangkan dengan nilai pada titik waktu ke-(t-12).

#plot time series data hasil differencing 
ts.plot(diffnonmus_mus12)

Dari plot yang di hasilkan sudah stasioner

#cek kestasioneran data
adf.test(differencing1)
## Warning in adf.test(differencing1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  differencing1
## Dickey-Fuller = -7.0177, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
# plot Autokorelasi data hasil differencing
par(mfrow=c(1,2))
acf(diffnonmus_mus12, lag.max = 75)
pacf(diffnonmus_mus12, lag.max = 75) 

#terdapat beberapa dugaan model ARIMA sementara yaitu ARIMA(0,1,1)(0,1,1) periode=12 atau
#ARIMA(0,1,1)(1,1,1) atau ARIMA(0,1,1)(1,1,0)
#ARIMA(0,1,1)(0,1,1) 
fit1<-arima(penumpang, order=c(0,1,1), seasonal = list(order=c(0,1,1), periode=12), method = "ML")
fit1
## 
## Call:
## arima(x = penumpang, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     periode = 12), method = "ML")
## 
## Coefficients:
##           ma1     sma1
##       -0.3087  -0.1074
## s.e.   0.0890   0.0828
## 
## sigma^2 estimated as 135.4:  log likelihood = -507.5,  aic = 1021

ARIMA(0,1,1)(0,1,1):

Ini adalah model ARIMA dengan komponen non-musiman (0,1,1) yang sama seperti sebelumnya. Komponen musiman (0,1,1) menunjukkan bahwa tidak ada komponen autoregressive musiman, tetapi ada differencing musiman pertama dan moving average musiman order 1.

#Diagnosis checking
#1. uji kesignifikanan parameter : uji t
coeftest(fit1)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.308674   0.088950 -3.4702 0.0005201 ***
## sma1 -0.107446   0.082817 -1.2974 0.1944950    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#2. pengujian residual apakah sudah white noise?
Box.test(fit1$residuals, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  fit1$residuals
## X-squared = 0.0020647, df = 1, p-value = 0.9638
#3.pengujian residual apakah sudah berdistribusi normal?
shapiro.test(fit1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.97603, p-value = 0.0125
#shapiro-francia normality test
sf.test(fit1$residuals)
## 
##  Shapiro-Francia normality test
## 
## data:  fit1$residuals
## W = 0.97088, p-value = 0.004843
#ARIMA(0,1,1)(1,1,0)
fit2<-arima(penumpang, order=c(0,1,1), seasonal = list(order=c(1,1,0), periode=12), method = "ML")
fit2
## 
## Call:
## arima(x = penumpang, order = c(0, 1, 1), seasonal = list(order = c(1, 1, 0), 
##     periode = 12), method = "ML")
## 
## Coefficients:
##           ma1     sar1
##       -0.3096  -0.1459
## s.e.   0.0885   0.0974
## 
## sigma^2 estimated as 134.7:  log likelihood = -507.21,  aic = 1020.43

ARIMA(0,1,1)(1,1,0):

Ini adalah model ARIMA dengan komponen non-musiman (0,1,1) yang sama seperti sebelumnya. Komponen musiman (1,1,0) menunjukkan bahwa terdapat komponen autoregressive musiman (order 1), differencing musiman pertama, tetapi tidak ada komponen moving average musiman.

#Diagnosis checking
#1. uji kesignifikanan parameter : uji t
coeftest(fit2)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.309644   0.088520  -3.498 0.0004687 ***
## sar1 -0.145874   0.097445  -1.497 0.1343964    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#2. pengujian residual apakah sudah white noise?
Box.test(fit2$residuals, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  fit2$residuals
## X-squared = 0.0017501, df = 1, p-value = 0.9666
#3.pengujian residual apakah sudah berdistribusi normal?
shapiro.test(fit2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit2$residuals
## W = 0.97477, p-value = 0.009215
#shapiro-francia normality test
sf.test(fit2$residuals)
## 
##  Shapiro-Francia normality test
## 
## data:  fit2$residuals
## W = 0.96971, p-value = 0.003831
#ARIMA(0,1,1)(1,1,1)
fit3<-arima(penumpang, order=c(0,1,1), seasonal = list(order=c(1,1,1), periode=12), method = "ML")
fit3
## 
## Call:
## arima(x = penumpang, order = c(0, 1, 1), seasonal = list(order = c(1, 1, 1), 
##     periode = 12), method = "ML")
## 
## Coefficients:
##           ma1     sar1    sma1
##       -0.3439  -0.9320  0.8459
## s.e.   0.0880   0.2159  0.3289
## 
## sigma^2 estimated as 130.7:  log likelihood = -506.16,  aic = 1020.33

ARIMA(0,1,1)(1,1,1):

Ini adalah model ARIMA dengan komponen non-musiman (0,1,1) yang terdiri dari differencing pertama dan moving average order 1. Komponen musiman (1,1,1) menunjukkan bahwa terdapat komponen autoregressive musiman (order 1), differencing musiman pertama, dan moving average musiman order 1.

#Diagnosis checking
#1. uji kesignifikanan parameter : uji t
coeftest(fit3)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)    
## ma1  -0.34387    0.08804 -3.9058 9.39e-05 ***
## sar1 -0.93199    0.21588 -4.3172 1.58e-05 ***
## sma1  0.84590    0.32893  2.5716  0.01012 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#2. pengujian residual apakah sudah white noise?
Box.test(fit3$residuals, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  fit3$residuals
## X-squared = 0.0022009, df = 1, p-value = 0.9626
#3.pengujian residual apakah sudah berdistribusi normal?
shapiro.test(fit3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit3$residuals
## W = 0.97594, p-value = 0.01223
#shapiro-francia normality test
sf.test(fit3$residuals)
## 
##  Shapiro-Francia normality test
## 
## data:  fit3$residuals
## W = 0.97089, p-value = 0.004854

dari ketiga metode ARIMA tersebut yang paling bagus di gunakan untuk data ini adalah ARIMA ARIMA(0,1,1)(1,1,1) karna p-valuenya > 0.01

#forecasting untuk 12 tahap kedepan berdasarkan model terbaik
forcesting <- forecast(penumpang, model=fit3, h=12)
forcesting
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 145       450.4486 435.7872 465.1100 428.0259 472.8713
## 146       425.6826 408.1481 443.2171 398.8659 452.4993
## 147       460.2926 440.2936 480.2917 429.7068 490.8785
## 148       499.2244 477.0329 521.4159 465.2855 533.1634
## 149       511.2022 487.0162 535.3882 474.2129 548.1915
## 150       569.5118 543.4837 595.5400 529.7052 609.3185
## 151       657.0813 629.3330 684.8296 614.6440 699.5186
## 152       642.4121 613.0443 671.7799 597.4979 687.3263
## 153       547.6304 516.7278 578.5330 500.3690 594.8919
## 154       498.0593 465.6947 530.4240 448.5618 547.5568
## 155       428.8696 395.1061 462.6331 377.2328 480.5064
## 156       472.3987 437.2921 507.5053 418.7077 526.0896
plot(forcesting, main="plot hasil peramalan")