Tugas Model SARIMA Menggunakan Dataset Walmart

SAYYID RAFIF

Data processing

# Membaca data Walmart
data_walmart <- read.csv("~/topik statistika/Walmart.csv", sep = ";")

# Menggabungkan penjualan dari semua toko berdasarkan Tanggal
data_gabungan <- aggregate(Weekly_Sales ~ Date, data = data_walmart, FUN = sum)

# Mengubah format teks tanggal menjadi objek Date dan mengurutkannya dari yang tertua
data_gabungan$Date <- as.Date(data_gabungan$Date, format = "%d/%m/%Y")
data_gabungan <- data_gabungan[order(data_gabungan$Date), ]

# Mengubah Weekly_Sales menjadi objek time series (ts)
data.ts <- ts(cbind(Yt = data_gabungan$Weekly_Sales), start = c(2010, 5), frequency = 52)

Visualisasi data time series

ts.plot(data.ts[,"Yt"], type="l", ylab="Yt", col="blue")
title(main = "Time Series Plot of Yt", cex.sub = 0.8)
points(data.ts[,"Yt"], pch = 20, col = "blue")

Mengecek seasonality data time series

seasonplot(data.ts[,"Yt"], 5, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 8, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 10, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 12, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

Berdasarkan percobaan plot seasonal di atas dapat diketahui bahwa tidak menunjukkan pola data yang seasonal maka akan dilakukan analisis model ARIMA reguler.

Split data 80% train dan 20% test

train.Yt <- subset(data.ts[,"Yt"], start = 1, end = 114)
test.Yt  <- subset(data.ts[,"Yt"], start = 115, end = 143)

Visualisasi data training yang akan digunakan untuk mencari model terbaik

ts.plot(train.Yt, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Train Time Series Plot of Yt", cex.sub = 0.8)
points(train.Yt, pch = 20, col = "blue")

Visualisasi data testing yang akan digunakan untuk mencari model terbaik

ts.plot(test.Yt, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Test Time Series Plot of Yt", cex.sub = 0.8)
points(test.Yt, pch = 20, col = "blue")

Model SARIMA

1.Cek kestasioneran data

Mengidentifikasi kestasioneran data dilakukan dengan uji ADF. Uji Augmented Dickey-Fuller (ADF) adalah salah satu alat identifikasi untuk menguji kestasioneran data dan memiliki hipotesis sebagai berikut:

\(𝐻_0:\)Data tidak stasioner (unit roots mendekati 1)

\(𝐻_1:\)Data stasioner (unit roots tidak mendekati 1)

cat("\n--- Uji ADF Data Train ---\n")

## 
## --- Uji ADF Data Train ---

print(adf.test(train.Yt))

## Warning in adf.test(train.Yt): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.Yt
## Dickey-Fuller = -5.3396, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5\%\) di atas dapat diketahui bersama bahwa \(p-value = 0.01 < \alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_0\) artinya data yang digunakan merupakan data yang stasioner.

2.Identifikasi kandidat model

Identifikasi kandidat model diperoleh berdasarkan nilai p dan q dimana nilai d=0

acf(train.Yt, lag.max = 20, col = "blue")

Berdasarkan plot ACF terlihat bahwa plot cut-off di lag 5, ARIMA(0,0,5)

pacf(train.Yt, lag.max = 20, col = "blue")

Berdasarkan plot PACF di atas kandidat model yang diperoleh adalah ARIMA(5,0,0)

eacf(train.Yt)

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

Berdasarkan hasil EACF di atas kandidat model yang diperoleh adalah ARIMA(1,0,5), ARIMA(2,0,6) dan ARIMA(3,0,5)

Kandidat Model

Diperoleh beberapa kandidat model yaitu: ARIMA(0,0,5), ARIMA(5,0,0), ARIMA(1,0,5), ARIMA(2,0,6) dan ARIMA(3,0,5).

# ARIMA(0,0,5)
arima005 <- arima(train.Yt, order = c(0,0,5), include.mean = TRUE, method = "ML")
arima005

## 
## Call:
## arima(x = train.Yt, order = c(0, 0, 5), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ma1     ma2     ma3     ma4      ma5   intercept
##       0.4274  0.2952  0.0088  0.4697  -0.2145  47313344.1
## s.e.  0.0992  0.0899  0.1076  0.1134   0.1117    910174.5
## 
## sigma^2 estimated as 2.424e+13:  log likelihood = -1920.31,  aic = 3852.62

# ARIMA(5,0,0)
arima500 <- arima(train.Yt, order = c(5,0,0), include.mean = TRUE, method = "ML")
arima500

## 
## Call:
## arima(x = train.Yt, order = c(5, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2      ar3     ar4      ar5  intercept
##       0.3649  0.0810  -0.0337  0.2640  -0.3829   47260769
## s.e.  0.0859  0.0884   0.0886  0.0882   0.0846     679450
## 
## sigma^2 estimated as 2.591e+13:  log likelihood = -1922.77,  aic = 3857.53

# ARIMA(1,0,5)
arima105 <- arima(train.Yt, order = c(1,0,5), include.mean = TRUE, method = "ML")
arima105

## 
## Call:
## arima(x = train.Yt, order = c(1, 0, 5), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ma1      ma2      ma3     ma4      ma5   intercept
##       0.8107  -0.4163  -0.1595  -0.2186  0.4027  -0.6083  47142331.0
## s.e.  0.0742   0.1286   0.0944   0.0837  0.1803   0.1642    175130.8
## 
## sigma^2 estimated as 2.304e+13:  log likelihood = -1918,  aic = 3850.01

# ARIMA2,0,6)
arima206 <- arima(train.Yt, order = c(2,0,6), include.mean = TRUE, method = "ML")
arima206

## 
## Call:
## arima(x = train.Yt, order = c(2, 0, 6), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2      ma1      ma2      ma3     ma4      ma5     ma6
##       0.8119  0.0289  -0.4377  -0.1523  -0.2652  0.4255  -0.6540  0.0837
## s.e.  0.7244  0.5831   0.7186   0.2834   0.1465  0.2015   0.3229  0.4706
##        intercept
##       47148770.8
## s.e.    186926.9
## 
## sigma^2 estimated as 2.276e+13:  log likelihood = -1917.68,  aic = 3853.37

# ARIMA(3,0,5)
arima305 <- arima(train.Yt, order = c(3,0,5), include.mean = TRUE, method = "ML")
arima305

## 
## Call:
## arima(x = train.Yt, order = c(3, 0, 5), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2     ar3      ma1      ma2      ma3     ma4      ma5
##       0.6361  0.1202  0.0575  -0.2614  -0.1854  -0.3180  0.3618  -0.5970
## s.e.  0.2399  0.1490  0.1921   0.2321   0.0970   0.1613  0.1680   0.1036
##        intercept
##       47150555.9
## s.e.    190352.1
## 
## sigma^2 estimated as 2.273e+13:  log likelihood = -1917.65,  aic = 3853.3

AICKandidatModel <- c(3852.62,3857.53 ,3850.01 ,3853.37 ,3853.3 )
KndidatModelARIMA <- c("ARIMA(0,0,5)", "ARIMA(5,0,0)", "ARIMA(1,0,5)", "ARIMA(2,0,6)", "ARIMA(3,0,5)")
compmodelARIMA <- cbind(KndidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA

##   Kandidat Model Nilai AIC
## 1   ARIMA(0,0,5)   3852.62
## 2   ARIMA(5,0,0)   3857.53
## 3   ARIMA(1,0,5)   3850.01
## 4   ARIMA(2,0,6)   3853.37
## 5   ARIMA(3,0,5)    3853.3

Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu ARIMA(1,0,5)

Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi parameternya, dimana:

\(𝐻_0:\)parameter \(= 0\)

\(𝐻_1:\)parameter \(\neq 0\)

Pengambilan keputusan yaitu Tolak \(H_0\) jika \(|t_{hitung}| > t_{\frac{\alpha}{2},(n-1)}\)

T-hitung parameter dugaan model \(ARIMA(1,0,5)\)

\[t_{hitung} = \frac{parameterestimasi}{SEparameterestimasi}\]

ar_phi1      <- 0.8107 / 0.0742
ar_phi1

## [1] 10.92588

ma_theta1    <- -0.4163 / 0.1286
ma_theta1

## [1] -3.23717

ma_theta2    <- -0.1595 / 0.0944
ma_theta2

## [1] -1.689619

ma_theta3    <- -0.2186 / 0.0837
ma_theta3

## [1] -2.611708

ma_theta4    <- 0.4027 / 0.1803
ma_theta4

## [1] 2.2335

ma_theta5    <- -0.6083 / 0.1642
ma_theta5

## [1] -3.704629

intercept_mu <- 47142331.0 / 175130.8
intercept_mu

## [1] 269.1836

# --- Menggabungkan ke dalam Tabel ---
prmtrdugaanarima105 <- c("ar1", "ma1", "ma2", "ma3", "ma4", "ma5", "intercept")

# Menggunakan fungsi abs() untuk memutlakkan nilai t-hitung
thitungprmtrarima105 <- c(abs(ar_phi1), abs(ma_theta1), abs(ma_theta2), abs(ma_theta3), abs(ma_theta4), abs(ma_theta5), abs(intercept_mu))

# Menentukan nilai T-Tabel
# n = 114 (dari data train Anda), sehingga derajat bebas (df) = 114 - 1 = 113
# Menggunakan taraf signifikansi alpha = 5% (uji dua arah, sehingga 0.975)
ttabel_val <- qt(0.975, df = 113)

# Mengulang nilai t-tabel sebanyak 7 kali agar sejajar dengan jumlah parameter
ttabel <- rep(ttabel_val, 7) 

# Logika pengambilan keputusan: Jika |t-hitung| > t-tabel, maka Signifikan
keputusan <- ifelse(thitungprmtrarima105 > ttabel, "Signifikan", "Tidak Signifikan")

# Membuat Data Frame untuk merapikan hasil
tksignimodelarima105 <- cbind(prmtrdugaanarima105, thitungprmtrarima105, ttabel, keputusan)
colnames(tksignimodelarima105) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodelarima105 <- as.data.frame(tksignimodelarima105)

# Menampilkan tabel hasil pengujian
tksignimodelarima105

##   Parameter dugaan         T-Hitung          T-Tabel        Keputusan
## 1              ar1 10.9258760107817 1.98118035941466       Signifikan
## 2              ma1 3.23716951788491 1.98118035941466       Signifikan
## 3              ma2  1.6896186440678 1.98118035941466 Tidak Signifikan
## 4              ma3 2.61170848267622 1.98118035941466       Signifikan
## 5              ma4 2.23349972268442 1.98118035941466       Signifikan
## 6              ma5 3.70462850182704 1.98118035941466       Signifikan
## 7        intercept 269.183553092888 1.98118035941466       Signifikan

Berdasarkan tabel di atas maka pengambilan keputusan dari paramater dugaan dari model ARIMA(1,0,5) semua parameternya signifikan, kecuali \(\theta_2\)

Diagnostic model

ARIMA105diag <- stats::arima(train.Yt, order = c(1,0,5), method = "ML")
checkresiduals(ARIMA105diag)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,5) with non-zero mean
## Q* = 6.1776, df = 17, p-value = 0.9919
## 
## Model df: 6.   Total lags used: 23

Berdasarkan pl0t di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, ditinjau dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

sisaan <- arima105$residuals

Uji formal normalitas data

\(H_0\): Sisaan berdistribusi normal

\(H_1\): Sisaan tidak berdistribusi normal

# Uji formal normalitas data
ks.test(sisaan,"pnorm")

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

Berdasarkan hasil pengujian, Normalitas Data: Hasil p-value < \(\alpha = 0.05\), maka keputusan yang diambil adalah tolak \(H_0\).Sisaan tidak menyebar normal

# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  sisaan
## t = -0.42282, df = 113, p-value = 0.6732
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1084614.3   703086.5
## sample estimates:
## mean of x 
## -190763.9

Uji Nilai Tengah Sisaan (T-Test) Pengujian ini menggunakan One Sample t-test untuk memastikan keacakan sisaan dengan hipotesis:

\(H_0 : \mu = 0\) (Nilai tengah sisaan sama dengan nol)

\(H_1 : \mu \neq 0\) (Nilai tengah sisaan tidak sama dengan nol)

Berdasarkan hasil uji-t nilai p-value > \(\alpha = 0.05\), maka terima \(H_0\). Ini menunjukkan bahwa secara statistik, nilai harapan (expected value) dari galat terbukti sama dengan nol.

# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 6.1776, df = 23, p-value = 0.9998

Uji Asumsi White Noise (Ljung-Box) Pengujian kebebasan sisaan dilakukan melalui Ljung-Box Test dengan rumusan hipotesis:

\(H_0\): Sisaan saling bebas / tidak terdapat autokorelasi (white noise)

\(H_1\): Sisaan terautokorelasi

Hasil pengujian memberikan nilai p-value > \(\alpha = 0.05\), maka diambil keputusan terima \(H_0\). Hal ini menyimpulkan bahwa sisaan bersifat saling bebas / tidak terdapat autokorelasi (white noise)

Peramalan

Selanjutnya yaitu melakukan peramalan hasil dari pemodelan dengan menggunakan data training, yaitu sebagai berikut:

Fitted model ARIMA(1,0,5)dengan data training

predictARIMA105 <- fitted(arima105)
fittedARIMA105 <- cbind(train.Yt, predictARIMA105)

Visualisasi hasil fitted model ARIMA(1,0,5)dengan data training

ts.plot(train.Yt, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(train.Yt, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA105)

forecasting <- predict(arima105, n.ahead = 20)
hasilforcasting <- as.data.frame(forecasting)
hasilforcasting

##        pred      se
## 1  48702534 4818692
## 2  47092290 5193778
## 3  43874630 5258650
## 4  48353965 5271206
## 5  43434400 5523019
## 6  44136195 5741736
## 7  44705161 5881067
## 8  45166441 5970877
## 9  45540414 6029179
## 10 45843607 6067195
## 11 46089414 6092053
## 12 46288698 6108337
## 13 46450264 6119016
## 14 46581250 6126026
## 15 46687445 6130628
## 16 46773541 6133652
## 17 46843341 6135638
## 18 46899930 6136944
## 19 46945809 6137801
## 20 46983005 6138365

Visualisasi hasil peramalan model ARIMA(1,0,5) 20 waktu ke depan

plot(arima105, n.ahead=20, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(train.Yt, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA105)

Fitting ARIMA(1,0,5) dengan data testing

arima105test <- arima(test.Yt, order = c(1,0,5), include.mean = TRUE, method = "ML")
arima105test

## 
## Call:
## arima(x = test.Yt, order = c(1, 0, 5), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ma1      ma2      ma3     ma4      ma5   intercept
##       0.7117  -0.3100  -0.2644  -0.2574  0.4845  -0.6520  46808300.3
## s.e.  0.2353   0.2898   0.2779   0.2532  0.3476   0.2869    207311.4
## 
## sigma^2 estimated as 1.735e+12:  log likelihood = -452.74,  aic = 919.47

predictARIMA105test <- fitted(arima105test)
fittedARIMA105test <- cbind(test.Yt, predictARIMA105test)

Visualisasi hasil fitted model ARIMA(1,0,5)dengan data testing

ts.plot(test.Yt, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(test.Yt, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA105test)

forecastingtest <- predict(arima105test, n.ahead = 20)
hasilforcastingtest <- as.data.frame(forecastingtest)
hasilforcastingtest

##        pred      se
## 1  47747172 1374482
## 2  47032498 1471239
## 3  46718205 1468902
## 4  47135002 1488194
## 5  47597371 1546940
## 6  47369850 1631194
## 7  47207933 1672247
## 8  47092703 1692659
## 9  47010698 1702903
## 10 46952339 1708068
## 11 46910807 1710678
## 12 46881250 1711999
## 13 46860216 1712667
## 14 46845246 1713005
## 15 46834593 1713177
## 16 46827012 1713263
## 17 46821617 1713307
## 18 46817777 1713330
## 19 46815045 1713341
## 20 46813100 1713347

Visualisasi hasil peramalan model ARIMA(1,0,5) 20 waktu ke depan

plot(arima105test, n.ahead=20, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(test.Yt, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA105test)

# Akurasi hasil peramalan

Akurasi hasil peramalan model ARIMA(1,0,5). AKurasi hasil peramalan model ARIMA(1,0,5)dibandingin dengan data aktual adalah Data testing.

# Menghitung akurasi
akurasi.arima105training <- accuracy(predictARIMA105, train.Yt)
akurasi.arima105testing <- accuracy(predictARIMA105test, test.Yt)

# Mengambil 5 nilai metrik utama dari output accuracy
ModelARIMA105train <- akurasi.arima105training[1, c("ME", "RMSE", "MAE", "MPE", "MAPE")]
ModelARIMA105test <- akurasi.arima105testing[1, c("ME", "RMSE", "MAE", "MPE", "MAPE")]

# Menggabungkan menjadi satu tabel
Akurasiall <- rbind(ModelARIMA105train, ModelARIMA105test)
rownames(Akurasiall) <- c("Training (ARIMA 1,0,5)", "Testing (ARIMA 1,0,5)")

# Menjadikan Data Frame
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall

##                                ME    RMSE     MAE         MPE     MAPE
## Training (ARIMA 1,0,5) -190763.91 4799802 3264247 -1.16731386 6.601918
## Testing (ARIMA 1,0,5)    60569.44 1317059 1008369  0.03963214 2.159604

Berdasarkan hasil fitting model dengan menggunakan metode ARIMA yang diperoleh model terbaik yaitu ARIMA(1,0,5). Pemilihan model terbaik ini didapatkan berdasarkan nilai akurasi AIC yang terkecil dari model tentatif dan diperoleh akurasi dari model berdasarkan nilai MAPE sebesar 6.601918 untuk data train dan MAPE sebesar 2.159604 untuk data testing. Berdasarkan range nilai MAPE bahwa kemampuan model dalam meramal:

rangeaku <- c("<10%", "10-20%", "20-50%", ">50%")
arti <- c("Kemampuan peramalan sangat baik", "Kemampuan peramalan baik", "Kemampuan peramalan layak", "Kemampuan peramalan buruk")
Akurasiall <- cbind(rangeaku, arti)
colnames(Akurasiall) <- c("Range MAPE", "Arti")
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall

##   Range MAPE                            Arti
## 1       <10% Kemampuan peramalan sangat baik
## 2     10-20%        Kemampuan peramalan baik
## 3     20-50%       Kemampuan peramalan layak
## 4       >50%       Kemampuan peramalan buruk

sehingga berdasarkan range tersebut model ARIMA(1,0,5)dengan nilai MAPE sebesar 6.601918 untuk data train dan MAPE sebesar 2.159604 untuk data testing memiliki tingkat peramalan yang layak.