Pendahuluan

ARIMA Musiman atau Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan pengembangan dari model Autoregressive Integrated Moving Average (ARIMA) pada data deret waktu yang memiliki pola musiman. Metode ini dipopulerkan oleh George Box dan Gwilym Jenskins sekitar tahun 1970-an.

Model ARIMA musiman menggabungkan faktor-faktor non-musiman (regular) dan musiman dalam model multiplikasi, dengan notasi sebagai berikut : \[ARIMA (p,d,q)×(P,D,Q)_s\] dengan :

p = non-seasonal AR order,
d = non-seasonal differencing,
q = non-seasonal MA order,
P = seasonal AR order,
D = seasonal differencing,
Q = seasonal MA order, dan
s = time span of repeating seasonal pattern.

Contoh :

\[ARIMA (0,0,1)×(0,0,1)_{12}\] Model ini terdiri dari non-musiman MA(1), musiman MA(1), tanpa pembedaan d=0,D=0 dan s=12

Tahapan Identifikasi Model ARIMA Musiman sama halnya seperti yang dilakukan pada model ARIMA regular atau model ARIMA tanpa musiman, yaitu :

Plot time series
Identifikasi model
Pendugaan parameter model
Seleksi Model
Melakukan peramalan menggunakan model terbaik

0. Pra Analisis

Install Library yang dibutuhkan

Silahkan install package-package berikut, jika Anda belum pernah melakukannya sebelumnya.

install.packages("tseries")
     install.packages("forecast")
     install.packages("TTR")
     install.packages("TSA")
     install.packages("astsa")
     install.packages("car")
     install.packages("portes")

Load Library yang digunakan

Jangan lupa, untuk Load library-library tersebut.

library("tseries")
     library("forecast")
     library("TTR")
     library("TSA")
     library("graphics")
     library("astsa")
     library("car")
     library("portes")

Input Data

USpass<-read.csv("data/P8_passenger.csv", sep=",")
     class(USpass)

## [1] "data.frame"

head(USpass)

##   Air.Passenger.Miles Bulan Tahun
     ## 1                2.42     1  1960
     ## 2                2.14     2  1960
     ## 3                2.28     3  1960
     ## 4                2.50     4  1960
     ## 5                2.44     5  1960
     ## 6                2.72     6  1960

Merubah class

Di sini dilakukan perubahan data USpass dari data.frame menjadi data ts atau Time Series

USpass.ts <- ts(USpass, frequency = 12, start = c(1960,1), end = c(1977,12))
     class(USpass.ts)

## [1] "mts"    "ts"     "matrix"

1. Plot Time Series

Membuat Plot Time Series

plot.ts(USpass.ts[,1], main = "Monthly U.S. Air Passenger Data", xlab = "Year",ylab = "Air Passenger Miles")

Plot data asal memperlihatkan pola musiman dengan s = 12, serta terdapat perilaku nonstasioner baik dalam rataan maupun ragam.

Plot per Musim

seasonplot(USpass.ts[,1],12,main="US Air Passenger", ylab="Air Passenger (miles)",year.labels = TRUE, col=rainbow(18))

Deskriptif USpass

monthplot(USpass.ts[,1],ylab="Air Passenger (miles)")

boxplot(USpass.ts[,1] ~ cycle(USpass.ts[,1]), xlab = "Month", ylab = "APM", main = "Monthly U.S. Air Passenger Data - Boxplot")

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

2. Identifikasi Model

Transformasi Logaritma

lnAPM=log(USpass.ts[,1])

Plot Hasil Transformasi Logaritma

plot.ts(lnAPM, main = "Monthly logged U.S. Air Passenger Data", xlab = "Year", ylab = "ln(APM)")

Transformasi logaritma berhasil mengatasi ketidakstasioneran dalam ragam meskipun ketidakstasioneran dalam rataan masih nampak.

Deskriptif ln US Pass

monthplot(lnAPM,ylab="Air Passenger (miles)")

boxplot(lnAPM ~ cycle(lnAPM), xlab = "Month", ylab = "ln(APM)", main = "Monthly U.S. Air Passenger Data - Boxplot")

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

Tes Kehomogenan Ragam

fligner.test(USpass.ts[,1]~ Tahun, data=USpass)

## 
     ##  Fligner-Killeen test of homogeneity of variances
     ## 
     ## data:  USpass.ts[, 1] by Tahun
     ## Fligner-Killeen:med chi-squared = 48.557, df = 17, p-value = 7.056e-05

Nilai-p kurang dari 5% sehinga tolak H0 dengan kata lain ragam tidak homogen

fligner.test(log(USpass.ts[,1])~ Tahun, data=USpass)

## 
     ##  Fligner-Killeen test of homogeneity of variances
     ## 
     ## data:  log(USpass.ts[, 1]) by Tahun
     ## Fligner-Killeen:med chi-squared = 7.6488, df = 17, p-value = 0.9735

Nilai-p lebih dari 5% sehingga tidak tolak H0 dengan kata lain ragam data setelah di ln kan homogen

Pembagian data lnAPM menjadi training-testing

Data setahun terakhir digunakan untuk data testing sementara data sebelumnya untuk membangun model

lpass<- subset(lnAPM,start=1,end=204)
     test <- subset(lnAPM,start=205,end=216)

Identifikasi kestasioneran data ln US Pass

acf0<-acf(lpass,main="ACF",lag.max=48,xaxt="n")
     axis(1, at=0:48/12, labels=0:48)

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

Plot ACF menunjukkan data ln US Pass tidak stasioner baik di series non musiman maupun musiman

Melakukan proses pembedaan pertama pada series non musiman data lnAPM

diff1.lpass=diff(lpass)
     plot(diff1.lpass, main = "Time series plot of lpas d=1")

Differencing non musiman d = 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non musimannya.

Identifikasi kestasioneran data ln US Pass

acf1 <- acf(diff1.lpass,lag.max=48,xaxt="n", main="ACF d1")
     axis(1, at=0:48/12, labels=0:48)

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

Plot ACF data nonseasonal differencing d = 1, mengkonfirmasi kestasioneran komponen non musiman (namun perhatikan lag 12,24, dst), pada series musiman belum stasioner

Melakukan proses pembedaan pada series musiman

diff12.lpass <- diff(lpass,12)
     plot(diff12.lpass, main = "Time series plot of ln(APM) D=12")
     axis(1, at=0:48/12, labels=0:48)

acf2<- acf(diff12.lpass,lag.max=48,xaxt="n", main="ACF d1D1")

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

Nonseasonal differencing D = 12 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen seasonalnya (namun tidak untuk komponen non musimannya).

Melakukan pembedaan pertama pada data yang telah dilakukan pembedaan pada series musimannya

diff12.1lpass <- diff(diff12.lpass,1)
     plot(diff12.1lpass, main = "Time series plot of ln(APM) d=1, D=12")

Identifikasi Model

acf2(diff12.1lpass,48)

##       [,1]  [,2]  [,3]  [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12]
     ## ACF  -0.18 -0.16 -0.07 -0.01  0.00 0.06 -0.06 -0.03 0.06  0.08  0.06 -0.50
     ## PACF -0.18 -0.20 -0.15 -0.10 -0.07 0.02 -0.06 -0.05 0.04  0.09  0.13 -0.47
     ##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
     ## ACF   0.08  0.12 -0.01 -0.03  0.14 -0.10 -0.01  0.02  0.05 -0.09  0.07  0.04
     ## PACF -0.12 -0.07 -0.13 -0.14  0.11 -0.01 -0.07 -0.05  0.13 -0.03  0.10 -0.23
     ##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
     ## ACF  -0.02 -0.03  0.04 -0.03 -0.13  0.14  0.05 -0.01 -0.10  0.08 -0.02 -0.01
     ## PACF -0.04 -0.07 -0.05 -0.14 -0.04  0.07  0.03  0.02  0.05  0.03  0.13 -0.14
     ##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
     ## ACF  -0.04  0.08 -0.03  0.12  0.05 -0.18 -0.08  0.09  0.08 -0.07 -0.02 -0.02
     ## PACF -0.12  0.03  0.01  0.04  0.05  0.02 -0.01  0.03  0.10  0.00  0.05 -0.13

Kedua komponen telah stasioner. Identifikasi komponen non musiman adalah \(ARIMA(0,1,2)\). Identifikasi komponen musiman adalah \(ARIMA(0,1,1)_{12}\), sehingga model tentatif adalah \(ARIMA(0,1,2)×(0,1,1)_{12}\).

eacf(diff12.1lpass)

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

auto.arima(lpass)

## Series: lpass 
     ## ARIMA(2,1,1)(1,1,1)[12] 
     ## 
     ## Coefficients:
     ##          ar1      ar2      ma1     sar1     sma1
     ##       0.3826  -0.1191  -0.6580  -0.1196  -0.6758
     ## s.e.  0.1506   0.0900   0.1393   0.1011   0.0842
     ## 
     ## sigma^2 estimated as 0.00375:  log likelihood=260.17
     ## AIC=-508.34   AICc=-507.88   BIC=-488.82

3. Pendugaan Parameter Model

\(ARIMA(0,1,2)×(0,1,1)_{12}\)

model1 <- Arima(lpass,order=c(0,1,2),seasonal=c(0,1,1))
     summary(model1)

## Series: lpass 
     ## ARIMA(0,1,2)(0,1,1)[12] 
     ## 
     ## Coefficients:
     ##           ma1      ma2     sma1
     ##       -0.2842  -0.2340  -0.7415
     ## s.e.   0.0694   0.0692   0.0548
     ## 
     ## sigma^2 estimated as 0.003744:  log likelihood=259.2
     ## AIC=-510.4   AICc=-510.18   BIC=-497.39
     ## 
     ## Training set error measures:
     ##                       ME       RMSE        MAE       MPE     MAPE      MASE
     ## Training set 0.001656771 0.05873648 0.03782158 0.1868852 2.270385 0.3036671
     ##                     ACF1
     ## Training set 0.008383106

\(ARIMA(2,1,1)×(1,1,1)_{12}\)

model2<-Arima(lpass,order=c(2,1,1),seasonal=c(1,1,1))
     summary(model2)

## Series: lpass 
     ## ARIMA(2,1,1)(1,1,1)[12] 
     ## 
     ## Coefficients:
     ##          ar1      ar2      ma1     sar1     sma1
     ##       0.3826  -0.1191  -0.6580  -0.1196  -0.6758
     ## s.e.  0.1506   0.0900   0.1393   0.1011   0.0842
     ## 
     ## sigma^2 estimated as 0.00375:  log likelihood=260.17
     ## AIC=-508.34   AICc=-507.88   BIC=-488.82
     ## 
     ## Training set error measures:
     ##                       ME       RMSE        MAE     MPE    MAPE     MASE
     ## Training set 0.001446393 0.05847019 0.03761942 0.18651 2.25693 0.302044
     ##                      ACF1
     ## Training set -0.004576573

\(ARIMA(1,1,2)×(1,1,1)_{12}\)

model3<-Arima(lpass,order=c(1,1,2),seasonal=c(1,1,1))
     summary(model3)

## Series: lpass 
     ## ARIMA(1,1,2)(1,1,1)[12] 
     ## 
     ## Coefficients:
     ##          ar1      ma1      ma2     sar1     sma1
     ##       0.1809  -0.4579  -0.1704  -0.1186  -0.6748
     ## s.e.  0.2803   0.2773   0.1327   0.1012   0.0843
     ## 
     ## sigma^2 estimated as 0.003754:  log likelihood=260.07
     ## AIC=-508.15   AICc=-507.69   BIC=-488.63
     ## 
     ## Training set error measures:
     ##                       ME       RMSE        MAE       MPE     MAPE    MASE
     ## Training set 0.001455749 0.05850619 0.03764383 0.1891525 2.256741 0.30224
     ##                      ACF1
     ## Training set -0.001596888

4. Seleksi Model

Pengujian Parameter Model

Di sini digunakan User Defined Function printstatarima

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(model1)

## 
     ## Coefficients:
     ##                  s.e.         t  sign.
     ## ma1   -0.2842  0.0694   -4.0951  1e-04
     ## ma2   -0.2340  0.0692   -3.3815  9e-04
     ## sma1  -0.7415  0.0548  -13.5310  0e+00

printstatarima(model2)

## 
     ## Coefficients:
     ##                  s.e.        t   sign.
     ## ar1    0.3826  0.1506   2.5405  0.0118
     ## ar2   -0.1191  0.0900  -1.3233  0.1872
     ## ma1   -0.6580  0.1393  -4.7236  0.0000
     ## sar1  -0.1196  0.1011  -1.1830  0.2382
     ## sma1  -0.6758  0.0842  -8.0261  0.0000

printstatarima(model3)

## 
     ## Coefficients:
     ##                  s.e.        t   sign.
     ## ar1    0.1809  0.2803   0.6454  0.5194
     ## ma1   -0.4579  0.2773  -1.6513  0.1002
     ## ma2   -0.1704  0.1327  -1.2841  0.2006
     ## sar1  -0.1186  0.1012  -1.1719  0.2426
     ## sma1  -0.6748  0.0843  -8.0047  0.0000

Model terbaik adalah model 1 karena semua dugaan parameternya berpengaruh nyata. Jika ada beberapa model yang semua dugaan parameternya nyata maka yang dipilih adalah yang nilai keakuratannya paling tinggi

Diagnostik Model

#checkresiduals(model1,lag=c(10,30))

Residual telah menyebar normal dan tidak terdapat autokorelasi

Alternatif Tampilan

tsdisplay(residuals(model1), lag.max=45, main='ARIMA(0,1,2)(0,1,1)12 Model Residuals')

library(portes)
     ljbtest <- LjungBox(residuals(model1),lags=seq(5,30,5),order=0,season=1,squared.residuals=FALSE)
     ljbtest

##  lags statistic df   p-value
     ##     5  1.032695  5 0.9598858
     ##    10  4.720549 10 0.9090470
     ##    15  5.980311 15 0.9800766
     ##    20  9.361220 20 0.9783332
     ##    25 12.929254 25 0.9773566
     ##    30 15.428446 30 0.9870838

Tidak terdapat autokorelasi pada sisaan

library(tseries)
     jarque.bera.test(residuals(model1))

## 
     ##  Jarque Bera Test
     ## 
     ## data:  residuals(model1)
     ## X-squared = 612.69, df = 2, p-value < 2.2e-16

Sisaan tidak menyebar normal

5. Melakukan Peramalan menggunakan Model Terbaik

Validasi Model

ramalan_arima = forecast(model1, 12)
     accuracy(ramalan_arima,test)

##                       ME       RMSE        MAE       MPE      MAPE      MASE
     ## Training set 0.001656771 0.05873648 0.03782158 0.1868852 2.2703855 0.3036671
     ## Test set     0.005948736 0.02820536 0.02200860 0.2135300 0.7963059 0.1767057
     ##                     ACF1 Theil's U
     ## Training set 0.008383106        NA
     ## Test set     0.551763208 0.2494055

plot(ramalan_arima)

forecast_arima <- cbind(ramalan_arima$mean,ramalan_arima$lower,ramalan_arima$upper)

Transformasi balik

forecast_arima <- cbind(USpass[205:216,1:2],exp(forecast_arima))
     forecast_arima <- as.data.frame(forecast_arima)

Plot ramalan perbulan

plot(forecast_arima[,2],forecast_arima[,1], xlab='Month', ylab='Air Passenger (miles)', type='l', col='black', lwd=1, ylim=c(10,25))
     lines(forecast_arima[,2],forecast_arima[,3], col='red', lwd=1)
     lines(forecast_arima[,2],forecast_arima[,5], col='blue', lwd=1)
     lines(forecast_arima[,2],forecast_arima[,7], col='green',lwd=1)
     legend("topleft", c("Actual","Forecast","Lower 95%","Upper 95%"), 
     col=c("black","red","blue","green"), lty=1,lwd=1,cex=0.6, inset=0.02, box.lty=0)

TUGAS

Lakukan pemodelan deret waktu untuk data musiman yang terdapat pada file latihan Praktikum 8.xlsx (Ada 3 gugus data)

Catatan :

Gunakan Tahapan Identifikasi Model ARIMA Musiman (plot, identifikasi, pendugaan dan seleksi model)
Gunakan 80% data yang ada sebagai data training, sisanya 20 % data untuk testing
Tugas dikumpulkan pada 24 April 2021 melalui newlms

IPB University - akbar.ritzki@apps.ipb.ac.id ↩︎

VIII : SARIMA-Seasonal ARIMA

Akbar Rizki¹

4/19/2021

Pendahuluan

0. Pra Analisis

Install Library yang dibutuhkan

Load Library yang digunakan

Input Data

Merubah class

1. Plot Time Series

Membuat Plot Time Series

Plot per Musim

Deskriptif USpass

2. Identifikasi Model

Transformasi Logaritma

Plot Hasil Transformasi Logaritma

Deskriptif ln US Pass

Tes Kehomogenan Ragam

Pembagian data lnAPM menjadi training-testing

Identifikasi kestasioneran data ln US Pass

Melakukan proses pembedaan pertama pada series non musiman data lnAPM

Identifikasi kestasioneran data ln US Pass

Melakukan proses pembedaan pada series musiman

Melakukan pembedaan pertama pada data yang telah dilakukan pembedaan pada series musimannya

Identifikasi Model

3. Pendugaan Parameter Model

\(ARIMA(0,1,2)×(0,1,1)_{12}\)

\(ARIMA(2,1,1)×(1,1,1)_{12}\)

\(ARIMA(1,1,2)×(1,1,1)_{12}\)

4. Seleksi Model

Pengujian Parameter Model

Diagnostik Model

Alternatif Tampilan

5. Melakukan Peramalan menggunakan Model Terbaik

Validasi Model

Transformasi balik

Plot ramalan perbulan

TUGAS

VIII : SARIMA-Seasonal ARIMA

Akbar Rizki1

4/19/2021

Pendahuluan

0. Pra Analisis

Install Library yang dibutuhkan

Load Library yang digunakan

Input Data

Merubah class

1. Plot Time Series

Membuat Plot Time Series

Plot per Musim

Deskriptif USpass

2. Identifikasi Model

Transformasi Logaritma

Plot Hasil Transformasi Logaritma

Deskriptif ln US Pass

Tes Kehomogenan Ragam

Pembagian data lnAPM menjadi training-testing

Identifikasi kestasioneran data ln US Pass

Melakukan proses pembedaan pertama pada series non musiman data lnAPM

Identifikasi kestasioneran data ln US Pass

Melakukan proses pembedaan pada series musiman

Melakukan pembedaan pertama pada data yang telah dilakukan pembedaan pada series musimannya

Identifikasi Model

3. Pendugaan Parameter Model

\(ARIMA(0,1,2)×(0,1,1)_{12}\)

\(ARIMA(2,1,1)×(1,1,1)_{12}\)

\(ARIMA(1,1,2)×(1,1,1)_{12}\)

4. Seleksi Model

Pengujian Parameter Model

Diagnostik Model

Alternatif Tampilan

5. Melakukan Peramalan menggunakan Model Terbaik

Validasi Model

Transformasi balik

Plot ramalan perbulan

TUGAS

Akbar Rizki¹