Hands On SARIMA

Lakukan pemodelan deret waktu untuk data yang diberikan menggunakan model ARIMA Musiman. Gunakan tahapan identifikasi model ARIMA Musiman (plot, identifikasi, pendugaan dan seleksi model).

Data

data <- read.csv("C:/Users/User/OneDrive/Dokumen/SEMESTER 6/Topik Statistika I/Walmart.csv")
head(data)
##   Store       Date Weekly_Sales Holiday_Flag Temperature Fuel_Price      CPI
## 1     1 05-02-2010      1643691            0       42.31      2.572 211.0964
## 2     1 12-02-2010      1641957            1       38.51      2.548 211.2422
## 3     1 19-02-2010      1611968            0       39.93      2.514 211.2891
## 4     1 26-02-2010      1409728            0       46.63      2.561 211.3196
## 5     1 05-03-2010      1554807            0       46.50      2.625 211.3501
## 6     1 12-03-2010      1439542            0       57.79      2.667 211.3806
##   Unemployment
## 1        8.106
## 2        8.106
## 3        8.106
## 4        8.106
## 5        8.106
## 6        8.106
length(data$Weekly_Sales)
## [1] 6435

Forming Data Time Series

Pada tahap ini dilakukan perubahan data Weekly_Sales dari data.frame menjadi data ts atau time series menggunakan fungsi ts.

data.ts <- ts(data, start=c(1,1), frequency=52)

Visualisasi Data Time Series

Data yang digunakan terdiri dari 6435 observasi atau yang diukur dalam 6435 Minggu. Hasil dari visualisasi data tersaji pada plot berikut:

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

Visualisasi per Musim

library(forecast)
seasonplot(data.ts[,"Weekly_Sales"], 52, main="Weekly Sales", ylab="Weekly Sales", year.labels = TRUE, col=rainbow(18))

Visualisasi Deskriptif

monthplot(data.ts[,"Weekly_Sales"],ylab="Weekly Sales", col="blue")

boxplot(data.ts[,"Weekly_Sales"]~data.ts[,"Date"],ylab="Zt",xlab="Week", col="blue")

Boxplot menunjukkan ukuran keragaman (panjang box) dan pemusatan (median) yang berbeda pada masing-masing minggu.

Identifikasi Model

Berdasarkan plot time series di atas bahwa data kemungkinan tidak stasioner terhadap ragam maupun rataan, maka akan dilakukan uji asumsi formal.

Uji Kehomogenan Ragam

Uji asumsi formal terhadap kehomogenan ragam yang digunakan yaitu Fligner-Killen test, dimana: \(H_0\) : Ragam homogen

\(H_1\) : Ragam tidak homogen

library(car)
## Loading required package: carData
fligner.test(Weekly_Sales ~ Date, data=data.ts)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Weekly_Sales by Date
## Fligner-Killeen:med chi-squared = 266.01, df = 142, p-value = 1.4e-09

Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi \(\alpha = 0.05\), nilai \(p−value\) yang diperoleh yaitu \(1.4\times10^{-9}\), maka \(p−value=1.4\times10^{-9}<\alpha=0.05\) sehingga tak tolak \(H_0\) atau dengan kata lain ragam data tidak stasioner. Sehingga perlu dilakukan proses transformasi data.

Transformasi Data

logdata <- log10(data.ts[,"Weekly_Sales"])
ts.plot(logdata, type="l", ylab="Log Weekly Sales", main = "Weekly logged Weekly Sales Data", xlab = "Year", col="blue")

Berdasarkan plot di atas transformasi logaritma tidak berhasil mengatasi ketidakstasioneran dalam ragam dan ketidakstasioneran dalam rataan masih nampak.

Deskriptif Data

monthplot(logdata,ylab="Weekly Sales", col="blue")

boxplot(logdata ~ cycle(logdata), col="blue", xlab = "Week", ylab = "log Weekly Sales", main = "Weekly Sales Data - Boxplot")

Boxplot menunjukkan ukuran keragaman (panjang box) hampir sama akan tetapi ukuran pemusatan (median) yang berbeda pada masing-masing minggu.

Uji Kehomogenan Ragam

\(H_0\) : Ragam homogen

\(H_1\) : Ragam tidak homogen

fligner.test(logdata~ Date, data=data.ts)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  logdata by Date
## Fligner-Killeen:med chi-squared = 24.867, df = 142, p-value = 1

Hasil uji Fligner-Killeen menunjukkan bahwa ragam menjadi stasioner setelah dilakukan transformasi log, karena nilai \(p=value\) lebih dari 5% sehinga terima H0 atau ragam homogen.

Berdasarkan hasil Uji Fligner-Killeen di atas transformasi logaritma berhasil mengatasi ketidakstasioneran dalam ragam sehingga data yang digunakan menggunakan data Weekly Sales yang telah ditransformasikan.

Splitting Data

Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebanyak 80% dan bagian kedua sebagai data testing sebesar 20%.

# jumlah data
n <- length(logdata)

# proporsi
train_size <- floor(0.8 * n)

# split
train <- logdata[1:train_size]
test  <- logdata[(train_size+1):n]

Visual Data Train

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

## Visual Data Test

ts.plot(test, col="blue", ylab = "Weekly Sales", xlab = "Week")
title(main = "Time Series Testing Plot of Weekly Sales", cex.sub = 0.8)
points(test, pch = 20, col = "blue")

# Identifikasi Kestasioneran Data

acf0 <- acf(train, lag.max = 52)

acf0$lag <- acf0$lag * 52
acf0.1 <- as.data.frame(cbind(acf0$acf,acf0$lag))
acf0.2 <- acf0.1[which(acf0.1$V2%%52==0),]
barplot(height = acf0.2$V1, 
names.arg=acf0.2$V2, ylab="ACF", xlab="Lag")

Plot ACF menunjukkan data Weekly Sales tidak stasioner sehingga perlu dilakukan proses differencing.

Differencing Data

d1 <- diff(train,52)
ts.plot(d1, type="l", ylab="d1 Weekly Sales", col="blue")

## Identifikasi Kestasioneran Data Hasil Differencing

acf1 <- acf(d1, lag.max = 52)

acf2 <- acf1$lag <- acf1$lag * 52
acf1.1 <- as.data.frame(cbind(acf1$acf,acf1$lag))
acf1.2 <- acf1.1[which(acf1.1$V2%%52==0),]
barplot(height = acf1.2$V1, names.arg=acf1.2$V2, ylab="ACF", xlab="Lag")

Differencing

Melakukan differencing pertama pada data yang telah dilakukan differencing pada series seasonalya.

d1d1 <- diff(d1)
ts.plot(d1d1, type="l", ylab="d1 d1 Weekly Sales", col="blue")

Identifikasi Model

acf2 <- acf(d1d1,lag.max=52)

pacf2 <- pacf(d1d1,lag.max=52)

acf2$lag <- acf2$lag * 52
acf2.1 <- as.data.frame(cbind(acf2$acf,acf2$lag))
acf2.2 <- acf2.1[which(acf2.1$V2%%52==0),]
barplot(height = acf2.2$V1, names.arg=acf2.2$V2, ylab="ACF", xlab="Lag")

pacf2$lag <- pacf2$lag * 52
pacf2.1 <- as.data.frame(cbind(pacf2$acf,pacf2$lag))
pacf2.2 <- pacf2.1[which(pacf2.1$V2%%52==0),]
barplot(height = pacf2.2$V1, names.arg=pacf2.2$V2, ylab="PACF", xlab="Lag")

## EACF

library(TSA)
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
eacf(d1d1)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x o x o o o  o  o  x 
## 1 x x x x x x o x o o o  o  o  x 
## 2 x x x o x o x x o o o  o  o  x 
## 3 x x x x x x x o o o o  o  o  x 
## 4 x x x x x x x o o o o  o  o  x 
## 5 x x x x x x o o o o o  o  o  x 
## 6 x o x x x x o o o o o  o  o  x 
## 7 x o x x x x o o o o o  x  o  o

Berdasarkan plot di atas diperoleh beberapa kandidat model yaitu \(SARIMA(3,1,7)(1,1,0)_{52}\), \(SARIMA(4,1,7)(1,1,0)_{52}\), \(SARIMA(5,1,6)(1,1,0)_{52}\).

tmodel1 <- Arima(train,order=c(3,1,7),seasonal=c(1,1,0))
summary(tmodel1)
## Series: train 
## ARIMA(3,1,7) 
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     ma2      ma3     ma4      ma5
##       -0.7067  -0.5237  -0.1544  0.4215  0.2849  -0.1192  0.0681  -0.1663
## s.e.   0.1879   0.1586   0.0624  0.1879  0.1406   0.0880  0.0275   0.0351
##           ma6      ma7
##       -0.1196  -0.1207
## s.e.   0.0421   0.0505
## 
## sigma^2 = 0.003241:  log likelihood = 7452.07
## AIC=-14882.14   AICc=-14882.08   BIC=-14810.13
## 
## Training set error measures:
##                         ME       RMSE        MAE         MPE      MAPE
## Training set -0.0002810198 0.05687318 0.03310638 -0.01149598 0.5529331
##                   MASE          ACF1
## Training set 0.9556843 -0.0002092528
tmodel2 <- Arima(train,order=c(4,1,7),seasonal=c(1,1,0))
summary(tmodel2)
## Series: train 
## ARIMA(4,1,7) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4     ma1     ma2      ma3     ma4
##       -0.7674  -0.5670  -0.1713  -0.0223  0.4822  0.3104  -0.1174  0.0774
## s.e.   0.3347   0.2177   0.0966   0.0892  0.3343  0.1605   0.0901  0.0428
##           ma5      ma6      ma7
##       -0.1649  -0.1288  -0.1313
## s.e.   0.0371   0.0668   0.0618
## 
## sigma^2 = 0.003242:  log likelihood = 7452.14
## AIC=-14880.28   AICc=-14880.22   BIC=-14801.72
## 
## Training set error measures:
##                         ME      RMSE        MAE         MPE      MAPE      MASE
## Training set -0.0002802022 0.0568724 0.03310165 -0.01147578 0.5528519 0.9555476
##                       ACF1
## Training set -0.0001327585
tmodel3 <- Arima(train,order=c(5,1,6),seasonal=c(1,1,0))
summary(tmodel3)
## Series: train 
## ARIMA(5,1,6) 
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ar5      ma1     ma2      ma3
##       0.2209  -0.1445  -0.0854  0.0867  -0.0951  -0.5066  0.1697  -0.0454
## s.e.  0.2847   0.0541   0.0713  0.0727   0.0561   0.2846  0.1105   0.0683
##          ma4      ma5     ma6
##       0.1042  -0.1932  0.0188
## s.e.  0.0583   0.0510  0.0750
## 
## sigma^2 = 0.003243:  log likelihood = 7451.4
## AIC=-14878.8   AICc=-14878.74   BIC=-14800.25
## 
## Training set error measures:
##                         ME       RMSE        MAE         MPE      MAPE
## Training set -0.0002737308 0.05688054 0.03307371 -0.01130882 0.5523603
##                   MASE          ACF1
## Training set 0.9547411 -0.0002301638

Model Terbaik

AICKandidatModel <- c(tmodel1$aic, tmodel2$aic, tmodel3$aic)
AICcKandidatModel <- c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc)
BICKandidatModel <- c(tmodel1$bic, tmodel2$bic, tmodel3$bic)
KandidatModelARIMA <- c("ARIMA(3,1,7)(1,1,0)52", "ARIMA(4,1,7)(1,1,0)52", "ARIMA(5,1,6)(1,1,0)52")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##          Kandidat Model         Nilai AIC        Nilai AICc         Nilai BIC
## 1 ARIMA(3,1,7)(1,1,0)52 -14882.1356728813  -14882.084261002 -14810.1278105851
## 2 ARIMA(4,1,7)(1,1,0)52 -14880.2759733634  -14880.215202035 -14801.7219417675
## 3 ARIMA(5,1,6)(1,1,0)52 -14878.8041569861 -14878.7433856577 -14800.2501253902

Model terbaik dipilih berdasarkan nilai AIC dan AICc terkecil dari kandidat model. Oleh karena itu, model terbaik yang diperoleh yaitu \(SARIMA(5,1,6)(1,1,0)_{52}\)

Dibandingan auto.arima

auto.arima(train)
## Series: train 
## ARIMA(5,1,2) 
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ar5      ma1    ma2
##       0.3482  -0.1729  -0.1264  0.1743  -0.2691  -0.6350  0.236
## s.e.  0.0538   0.0556   0.0220  0.0159   0.0148   0.0555  0.060
## 
## sigma^2 = 0.003253:  log likelihood = 7441.15
## AIC=-14866.31   AICc=-14866.28   BIC=-14813.94

Jika berdasarkan fungsi auto.arima model terbaik yang diperoleh yaitu \(ARIMA(5,1,2)(1,1,0)_{52}\).

Model terbaik yang didapatkan adalah: \(ARIMA(5,1,6)(1,1,0)_{52}\) dengan fungsi Arima

\(ARIMA(5,1,2)(1,1,0)_{52}\) dengan fungsi auto.arima

Berdasarkan model terbaik dari fungsi auto.arima yaitu \(ARIMA(5,1,2)(1,1,0)_{52}\) nilai AIC sebesar -14866.31 dan model \(ARIMA(5,1,6)(1,1,0)_{52}\) dengan fungsi Arima nilai AICnya sebesar -14878.8041569861 diperoleh bahwa nilai AIC terkecil yaitu model \(ARIMA(5,1,6)(1,1,0)_{52}\) -14878.8041569861. Oleh karena itu, model yang disarankan untuk dianalisis lebih lanjut yaitu dengan menggunakan model \(ARIMA(5,1,6)(1,1,0)_{52}\) yang mana ini adalah model terbaik (dengan nilai AIC terkecil) dari kandidat model berdasarkan ACF, PACF, dan EACF.

Pengujian Parameter Model

printstatarima <- function (x, digits = 4,se=TRUE){
       if (length(x$coef) > 0) {
         cat("\nCoefficients:\n")
         coef <- round(x$coef, digits = digits)
         if (se && nrow(x$var.coef)) {
           ses <- rep(0, length(coef))
           ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
           coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
           coef <- rbind(coef, s.e. = ses)
           statt <- coef[1,]/ses
           pval  <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
           coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
           coef <- t(coef)
         }
         print.default(coef, print.gap = 2)
       }
     }
printstatarima(tmodel3)
## 
## Coefficients:
##                 s.e.        t   sign.
## ar1   0.2209  0.2847   0.7759  0.4378
## ar2  -0.1445  0.0541  -2.6710  0.0076
## ar3  -0.0854  0.0713  -1.1978  0.2311
## ar4   0.0867  0.0727   1.1926  0.2331
## ar5  -0.0951  0.0561  -1.6952  0.0901
## ma1  -0.5066  0.2846  -1.7800  0.0751
## ma2   0.1697  0.1105   1.5357  0.1247
## ma3  -0.0454  0.0683  -0.6647  0.5063
## ma4   0.1042  0.0583   1.7873  0.0739
## ma5  -0.1932  0.0510  -3.7882  0.0002
## ma6   0.0188  0.0750   0.2507  0.8021

Penduga parameter \(\hat{\phi_1}=0.2209\), \(\hat{\phi_2}=-0.1445\), \(\hat{\phi_3}=-0.0854\), \(\hat{\phi_4}=0.0867\), \(\hat{\phi_5}=-0.0951\), dan \(\hat{\theta_1}=0.5066\), \(\hat{\theta_2}=-0.1697\), \(\hat{\theta_3}=0.0454\), \(\hat{\theta_4}=-0.1042\), \(\hat{\theta_5}=0.1932\), \(\hat{\theta_6}=-0.0188\) dengan nilai dugaan \(\sigma_e^2\) sebesar 0.003243. Model terbaik yang diperoleh yaitu model \(ARIMA(5,1,6)(1,1,0)_{52}\)

Maka \(X_t\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(5,1,6)(1,1,0)_{52}\)

  • \(p=5, d=1, q=6\)
  • \(P=1, D=1, Q=0\)

\[ \begin{aligned} \phi_p(B)\Phi_P(B^s)(1-B)^d(1-B^s)^DX_t&=\theta_q(B)\Theta_Q(B^s)e_t \\ \phi_5(B)\Phi_1(B^{52})(1-B)^1(1-B^{52})^1X_t&=\theta_6(B)\Theta_0(B^{52})e_t \\ \end{aligned} \]

Diagnostik Model

tsdisplay(residuals(tmodel3), lag.max=45, main='ARIMA(5,1,6)(1,1,0)52 Model Residuals', col="blue")

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

library(portes)
ljbtest2 <- LjungBox(residuals(tmodel3), lags = seq(5, 30, 5), fitdf = 12)
ljbtest2
##  lags    statistic df      p-value
##     5  0.007676207  0           NA
##    10  1.613681846  0           NA
##    15 13.566106185  3 3.559394e-03
##    20 30.093504300  8 2.034825e-04
##    25 32.497045853 13 2.026397e-03
##    30 49.744433165 18 8.249884e-05