Package

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.2
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.2
## Warning: package 'lubridate' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TSA)
## Warning: package 'TSA' was built under R version 4.3.3
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## 
## The following object is masked from 'package:readr':
## 
##     spec
## 
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## 
## The following object is masked from 'package:utils':
## 
##     tar
library(aTSA)
## Warning: package 'aTSA' was built under R version 4.3.2
## 
## Attaching package: 'aTSA'
## 
## The following object is masked from 'package:forecast':
## 
##     forecast
## 
## The following object is masked from 'package:graphics':
## 
##     identify
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Pendahuluan

Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan pengembangan dari model Autoregressive Integrated Moving Average (ARIMA) pada data deret waktu yang memiliki pola musiman.

Model ARIMA musiman menggabungkan faktor-faktor non-musiman (regular) dan musiman dalam model multiplikasi, dengan notasi \(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, s = time span of repeating seasonal pattern.

Tahapan identifikasi model SARIMA sama halnya seperti yang dilakukan pada model ARIMA regular atau model ARIMA non-seasonal, yaitu :

  1. Plot time series
  2. Identifikasi model
  3. Pendugaan parameter model
  4. Seleksi Model
  5. Melakukan peramalan menggunakan model terbaik

Impor Data

nott <- datasets::nottem
head(nott)
##       Jan  Feb  Mar  Apr  May  Jun
## 1920 40.6 40.8 44.4 46.7 54.1 58.5
class(nott)
## [1] "ts"

Eksplorasi Data

ts.plot(nott, type="l", xlab = "Year", ylab="Temperature F", col="blue")
title(main = "Time Series Plot Monthly Temperatures at Nottingham", cex.sub = 0.8)
points(nott, pch = 20, col = "blue")

dec.nott <- decompose(nott)
plot(dec.nott)

Secara eksplorasi, terlihat adanya kecenderungan data konstan dan perilaku berulang kecenderungan musiman dalam deret tersebut. Kecenderungan musiman dapat dilihat dengan lebih jelas dengan menampilkan deret waktu per tahun.

seasonplot(nott,12,main="Seasonal Plot of Temperature", ylab="Temperature F",
           year.labels = TRUE, col=rainbow(18))

Gambar menunjukkan bahwa temperatur tinggi pada bulan Mei, Juni, Juli, Agustus, September, dan Oktober dan rendah pada bulan Januari, Februari, Maret, April, November, dan Desember. Perilaku tersebut terus berulang dari tahun ke tahun.

monthplot(nott,ylab="Temperature F", col="blue")

library(ggplot2)
#-------------------------------------------------------------------------------
frame<-data.frame(values=as.matrix(nott), date=lubridate::year(zoo::as.Date(nott)))
ggplot(frame,aes(y=values,x=date,group=date))+
  geom_boxplot()

Berdasarkan hasil plot di atas dapat terlihat bahwa data memiliki pola yang hampir sama dari tahun ke tahun sehingga dapat disimpulkan bahwa periode musimannya adalah 12. Selain itu, apabila dilihat dari boxplot, terlihat bahwa data cenderung homogen dari tahun ke tahun. Untuk memastikan bahwa data homogen akan dilakukan uji homogenitas dengan \(fligner.test\).

Uji Homogenitas

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

fligner.test(values ~ date, data=frame)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  values by date
## Fligner-Killeen:med chi-squared = 8.8669, df = 19, p-value = 0.9756

Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi \(\alpha=5\%\) didapatkan p-value sebesar 0.9756. \(p-value=0.9756>\alpha=0.05\) sehingga tak tolak \(H_0\) atau dengan kata lain ragam data sudah stasioner.

Pembagian Data

Pembagian data dilakukan dengan mengambil sekitar 80% data awal (192 observasi) sebagai data latih dan 20% sisanya (48 observasi) sebagai data uji.

train.ts <- subset(nott,start=1,end=192)
test.ts <- subset(nott,start=193,end=240)

Plot Data Latih

autoplot(train.ts) + theme_bw() + xlab("Year") + ylab("Temp (F)")

Plot Data Uji

autoplot(test.ts) + theme_bw() + xlab("Year") + ylab("Temp (F)")

Non-Seasonal ARIMA

Kestasioneran Data

par(mfrow = c(1, 2))

acf0 <- acf(train.ts,main="ACF",lag.max=50,xaxt="n", col="blue")
axis(1, at=0:50/6, labels=0:50)

pacf0 <- pacf(train.ts,main="PACF",lag.max=50,xaxt="n", col="blue")
axis(1, at=0:50/6, labels=0:50)

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")

Berdasarkan plot deret sebelumnya diketahui bahwa perilaku deret berulang setiap 12 bulan, atau dikatakan bahwa deret memiliki periode musiman bulanan, sehingga \(s=12\). Perhatikan nilai fungsi autokorelasi pada lag-lag musiman (lag 12, 24, 36,…) dalam plot ACF contoh di atas. Tampak bahwa nilai autokorelasi pada lag-lag tersebut memiliki hubungan yang kuat.

tseries::adf.test(train.ts)
## Warning in tseries::adf.test(train.ts): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -11.744, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF.

Seasonal ARIMA

D1 <- diff(train.ts,12)
ts.plot(D1, type="l", ylab="D1 Temp", col="blue")

acf2<-acf(D1,lag.max=50,xaxt="n", main="ACF D1", col="blue")
axis(1, at=0:50/12, labels=0:50)

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")

ts.plot(D1, type="l", ylab="D1 Temp", col="blue")

Setelah pembedaan musiman tampak bahwa deret sudah tidak memiliki kecenderungan apapun. Selanjutnya penentuan ordo p, q dan P, Q dapat dilakukan menggunakan plot ACF dan PACF contoh dari deret hasil pembedaan musiman tersebut.

Identifikasi Model

acf3 <- acf(D1,lag.max=50,xaxt="n", main="ACF D1", col="blue")
axis(1, at=0:50/12, labels=0:50)

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

Berdasarkan plot ACF contoh lag 1 signifikan sehingga dipilih ordo q=1 , dan lag 12 serta 24 adalah lag musiman yang signifikan sehingga order Q=2.

pacf3 <- pacf(D1,lag.max=50,xaxt="n", main="PACF D1", col="blue")
axis(1, at=0:50/12, labels=0:50)

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

Plot PACF contoh menunjukkan cuts-off pada lag-1 sehingga ordo p=1, pada pola musimannya menunjukkan cuts-off pada lag-24 sehingga ordo P=2.

Model musiman yang dipilih untuk deret temperature di Nottingham adalah \(ARIMA(1,0,1)\times(2,1,2)_{12}\), \(ARIMA(1,0,1)\times(2,1,2)_{12}\), \(ARIMA(0,0,1)\times(2,1,2)_{12}\), \(ARIMA(1,0,0)\times(2,1,2)_{12}\), \(ARIMA(2,0,1)\times(2,1,2)_{12}\), \(ARIMA(3,0,2)\times(2,1,2)_{12}\). Ingat kembali bahwa model yang digunakan bersifat tentatif dan dapat berubah saat diagnostik model.

EACF

TSA::eacf(D1)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o x o  x  o  o 
## 1 x o o o o o o o o o o  x  x  o 
## 2 o o o o o o o o o o o  x  o  x 
## 3 o x o o o o o o o o o  x  o  x 
## 4 x x x o o o o o o o o  x  x  x 
## 5 x x o o o o o o o o o  x  x  x 
## 6 x x o o x o o o o o o  x  x  x 
## 7 x o o o o o o o o o o  x  o  x

Sehingga model tentatif yang diperoleh adalah: \(ARIMA(1,0,1)\times(2,1,2)_{12}\) \(ARIMA(0,0,1)\times(2,1,2)_{12}\) \(ARIMA(2,0,0)\times(2,1,2)_{12}\) \(ARIMA(2,0,1)\times(2,1,2)_{12}\) \(ARIMA(3,0,2)\times(2,1,2)_{12}\)

Pendugaan Parameter

#ARIMA(1,0,1)x(2,1,2)12
tmodel1 <- Arima(train.ts,order=c(1,0,1),seasonal=c(2,1,2))
summary(tmodel1)
## Series: train.ts 
## ARIMA(1,0,1)(2,1,2)[12] 
## 
## Coefficients:
##          ar1      ma1     sar1    sar2     sma1    sma2
##       0.4789  -0.1641  -0.1731  0.1916  -0.8288  0.0254
## s.e.  0.1735   0.1925   0.4618  0.1853   0.4650  0.3188
## 
## sigma^2 = 5.373:  log likelihood = -411.09
## AIC=836.18   AICc=836.83   BIC=858.53
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 0.05964773 2.206665 1.70528 -0.1873942 3.646091 0.6112115
##                     ACF1
## Training set -0.00253994
lmtest::coeftest(tmodel1)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1   0.478895   0.173453  2.7609 0.005763 **
## ma1  -0.164052   0.192542 -0.8520 0.394197   
## sar1 -0.173105   0.461776 -0.3749 0.707759   
## sar2  0.191585   0.185339  1.0337 0.301277   
## sma1 -0.828817   0.464996 -1.7824 0.074682 . 
## sma2  0.025434   0.318771  0.0798 0.936407   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ARIMA(0,0,1)x(2,1,2)12
tmodel2 <- Arima(train.ts,order=c(0,0,1),seasonal=c(2,1,2))
summary(tmodel2)
## Series: train.ts 
## ARIMA(0,0,1)(2,1,2)[12] 
## 
## Coefficients:
##          ma1     sar1    sar2     sma1     sma2
##       0.2731  -0.2164  0.1774  -0.7713  -0.0120
## s.e.  0.0644   0.4411  0.1818   0.4422   0.2973
## 
## sigma^2 = 5.503:  log likelihood = -413.5
## AIC=839   AICc=839.49   BIC=858.16
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.07727077 2.239553 1.728005 -0.1634677 3.693358 0.6193566
##                    ACF1
## Training set 0.04419382
lmtest::coeftest(tmodel2)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1   0.273069   0.064387  4.2411 2.225e-05 ***
## sar1 -0.216370   0.441058 -0.4906   0.62373    
## sar2  0.177395   0.181826  0.9756   0.32925    
## sma1 -0.771260   0.442176 -1.7442   0.08112 .  
## sma2 -0.012002   0.297305 -0.0404   0.96780    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ARIMA(2,0,0)x(2,1,2)12
tmodel3 <- Arima(train.ts,order=c(2,0,0),seasonal=c(2,1,2))
summary(tmodel3)
## Series: train.ts 
## ARIMA(2,0,0)(2,1,2)[12] 
## 
## Coefficients:
##          ar1     ar2     sar1    sar2     sma1    sma2
##       0.3085  0.0711  -0.1721  0.1893  -0.8302  0.0269
## s.e.  0.0754  0.0750   0.4726  0.1879   0.4758  0.3265
## 
## sigma^2 = 5.367:  log likelihood = -410.99
## AIC=835.99   AICc=836.64   BIC=858.34
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.05970602 2.205427 1.703295 -0.1866621 3.641121 0.6105001
##                   ACF1
## Training set 0.0040457
lmtest::coeftest(tmodel3)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.308460   0.075432  4.0892 4.328e-05 ***
## ar2   0.071149   0.075003  0.9486   0.34282    
## sar1 -0.172142   0.472610 -0.3642   0.71568    
## sar2  0.189308   0.187936  1.0073   0.31379    
## sma1 -0.830182   0.475803 -1.7448   0.08102 .  
## sma2  0.026886   0.326511  0.0823   0.93437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ARIMA(2,0,1)x(2,1,2)12
tmodel4 <- Arima(train.ts,order=c(2,0,1),seasonal=c(2,1,2))
summary(tmodel4)
## Series: train.ts 
## ARIMA(2,0,1)(2,1,2)[12] 
## 
## Coefficients:
##          ar1     ar2     ma1     sar1    sar2     sma1    sma2
##       0.0863  0.1473  0.2224  -0.1676  0.1891  -0.8356  0.0323
## s.e.  0.6055  0.2045  0.6081   0.4760  0.1893   0.4793  0.3289
## 
## sigma^2 = 5.395:  log likelihood = -410.93
## AIC=837.87   AICc=838.71   BIC=863.41
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 0.06082412 2.204834 1.70333 -0.1830591 3.641687 0.6105126
##                     ACF1
## Training set 0.003160852
lmtest::coeftest(tmodel4)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)  
## ar1   0.086344   0.605475  0.1426  0.88660  
## ar2   0.147316   0.204510  0.7203  0.47132  
## ma1   0.222419   0.608122  0.3657  0.71455  
## sar1 -0.167632   0.476025 -0.3521  0.72473  
## sar2  0.189104   0.189321  0.9989  0.31786  
## sma1 -0.835568   0.479292 -1.7433  0.08127 .
## sma2  0.032294   0.328912  0.0982  0.92179  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ARIMA(3,0,2)x(0,1,1)12
tmodel5 <- Arima(train.ts,order=c(3,0,2),seasonal=c(2,1,2))
summary(tmodel5)
## Series: train.ts 
## ARIMA(3,0,2)(2,1,2)[12] 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##          ar1      ar2     ar3      ma1     ma2     sar1     sar2     sma1
##       1.3625  -0.4686  0.0304  -1.0321  0.1649  -0.9599  -0.1914  -0.0509
## s.e.  0.0101      NaN     NaN   0.0095     NaN   0.0094   0.0057   0.0118
##         sma2
##       -0.415
## s.e.     NaN
## 
## sigma^2 = 5.563:  log likelihood = -412.02
## AIC=844.03   AICc=845.34   BIC=875.96
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.06759669 2.225812 1.720394 -0.1695233 3.690417 0.6166286
##                     ACF1
## Training set -0.02405128
lmtest::coeftest(tmodel5)
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##        Estimate Std. Error   z value  Pr(>|z|)    
## ar1   1.3625195  0.0101149  134.7039 < 2.2e-16 ***
## ar2  -0.4686344        NaN       NaN       NaN    
## ar3   0.0303833        NaN       NaN       NaN    
## ma1  -1.0320609  0.0095305 -108.2902 < 2.2e-16 ***
## ma2   0.1648730        NaN       NaN       NaN    
## sar1 -0.9598734  0.0093505 -102.6543 < 2.2e-16 ***
## sar2 -0.1914215  0.0057220  -33.4538 < 2.2e-16 ***
## sma1 -0.0509344  0.0118237   -4.3078 1.649e-05 ***
## sma2 -0.4150439        NaN       NaN       NaN    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICKandidatModel <- c(tmodel1$aic, tmodel2$aic, tmodel3$aic,
                      tmodel4$aic, tmodel5$aic)
AICcKandidatModel <- c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc,
                       tmodel4$aicc, tmodel5$aicc)
BICKandidatModel <- c(tmodel1$bic, tmodel2$bic, tmodel3$bic,
                      tmodel4$bic, tmodel5$bic)
KandidatModelARIMA <- c("ARIMA(1,0,1)(2,1,2)12", "ARIMA(0,0,1)(2,1,2)12",
                        "ARIMA(2,0,0)(2,1,2)12", "ARIMA(2,0,1)(2,1,2)12",
                        "ARIMA(3,0,2)(2,1,2)12")
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(1,0,1)(2,1,2)12   836.1819790342 836.833141824898 858.532676990432
## 2 ARIMA(0,0,1)(2,1,2)12  839.00105489069 839.486604023638 858.158795996031
## 3 ARIMA(2,0,0)(2,1,2)12 835.987909931285 836.639072721983 858.338607887517
## 4 ARIMA(2,0,1)(2,1,2)12 837.868959599341 838.711064862499 863.412614406463
## 5 ARIMA(3,0,2)(2,1,2)12 844.034439786801  845.33621493473 875.964008295703

Model terbaik berdasarkan nilai AIC dan AICc terkecil dari kandidat model yaitu \(ARIMA(2,0,0)\times(2,1,2)_{12}\).

model.auto.arima <- auto.arima(train.ts)
summary(model.auto.arima)
## Series: train.ts 
## ARIMA(1,0,2)(1,1,2)[12] with drift 
## 
## Coefficients:
##          ar1     ma1     ma2     sar1     sma1     sma2   drift
##       0.1970  0.1054  0.1271  -0.5710  -0.4468  -0.2068  0.0051
## s.e.  0.3199  0.3149  0.1193   0.1749   0.1907   0.1721  0.0060
## 
## sigma^2 = 5.389:  log likelihood = -410.74
## AIC=837.49   AICc=838.33   BIC=863.03
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.05607925 2.203592 1.690523 -0.428198 3.631101 0.6059223
##                     ACF1
## Training set 0.003376878
lmtest::coeftest(model.auto.arima)
## 
## z test of coefficients:
## 
##         Estimate Std. Error z value Pr(>|z|)   
## ar1    0.1970332  0.3199319  0.6159 0.537987   
## ma1    0.1053538  0.3149366  0.3345 0.737984   
## ma2    0.1270840  0.1193134  1.0651 0.286818   
## sar1  -0.5709741  0.1748664 -3.2652 0.001094 **
## sma1  -0.4468435  0.1907443 -2.3426 0.019148 * 
## sma2  -0.2068450  0.1721025 -1.2019 0.229413   
## drift  0.0051457  0.0060220  0.8545 0.392835   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model yang didapat jikan menggunakan fungsi \(auto.arima\) adalah \(ARIMA(1,0,2)\times(1,1,2)_{12}\) dengan nilai AIC 837.5 dan BIC 863.03. Tetapi pada model tersebut cukup tidak sederhana dan hanya 2 yang signifikan

Diagnostik Model

tsdisplay(residuals(tmodel3), lag.max=48, 
          main='ARIMA(2,0,0)(2,1,2)12 Model Residuals', col="blue")

#Eksplorasi
sisaan.model3 <- tmodel3$residuals
par(mfrow=c(2,2))
car::qqPlot(sisaan.model3)
## [1] 110  19
plot(c(1:length(sisaan.model3)),sisaan.model3)
acf(sisaan.model3)
pacf(sisaan.model3)

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:

Uji Formal

#1) Sisaan Menyebar Normal
ks.test(sisaan.model3,"pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.model3
## D = 0.18819, p-value = 2.483e-06
## alternative hypothesis: two-sided
#tak tolak H0 > sisaan menyebar normal
shapiro.test(sisaan.model3)
## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan.model3
## W = 0.99104, p-value = 0.2791
nortest::ad.test(sisaan.model3)
## 
##  Anderson-Darling normality test
## 
## data:  sisaan.model3
## A = 0.78115, p-value = 0.04191
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.2
## 
## Attaching package: 'tseries'
## The following objects are masked from 'package:aTSA':
## 
##     adf.test, kpss.test, pp.test
jarque.bera.test(sisaan.model3)
## 
##  Jarque Bera Test
## 
## data:  sisaan.model3
## X-squared = 0.65697, df = 2, p-value = 0.72

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS), Shapiro-Wilk, Anderson-Darling, dan Jarque Bera. Hipotesis pada uji kenormalan adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value jarque bera sebesar 0.72 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan menyebar normal. Namun, hasil ini berbeda dengan ks test dan Anderson-Darling test yang menunjukkan p-value < 5%. Hasil uji formal Shapiro-Wilk dan Jarque bera lebih sesuai dengan plot eksplorasi, juga untuk data time series lebih disarankan menggunakan uji Jarque Bera oleh Bu Yenni.

#2) Sisaan saling bebas/tidak ada autokorelasi
Box.test(sisaan.model3, type = "Ljung") 
## 
##  Box-Ljung test
## 
## data:  sisaan.model3
## X-squared = 0.003192, df = 1, p-value = 0.9549
#tak tolak H0 > sisaan saling bebas

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9549 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas.

#3) Sisaan homogen 
Box.test((sisaan.model3)^2, type = "Ljung")  
## 
##  Box-Ljung test
## 
## data:  (sisaan.model3)^2
## X-squared = 0.044548, df = 1, p-value = 0.8328
#tak tolak H0 > sisaan homogen

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.8328 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.

#4) Nilai tengah sisaan sama dengan nol 
t.test(sisaan.model3, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaan.model3
## t = 0.37428, df = 191, p-value = 0.7086
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2549423  0.3743544
## sample estimates:
##  mean of x 
## 0.05970602
#tak tolak h0 > nilai tengah sisaan sama dengan 0

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.7086 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini berbeda dengan eksplorasi.

Overfitting

Pertama, overfit pada model non-musimannya (p,q)

#ARIMA(3,0,0)x(2,1,2)12
tmodel1.ofq <- Arima(train.ts,order=c(3,0,0),seasonal=c(2,1,2))
summary(tmodel1.ofq)
## Series: train.ts 
## ARIMA(3,0,0)(2,1,2)[12] 
## 
## Coefficients:
##          ar1     ar2      ar3     sar1    sar2     sma1    sma2
##       0.3110  0.0832  -0.0385  -0.1702  0.1860  -0.8337  0.0333
## s.e.  0.0755  0.0785   0.0753   0.4881  0.1941   0.4912  0.3360
## 
## sigma^2 = 5.392:  log likelihood = -410.86
## AIC=837.73   AICc=838.57   BIC=863.27
## 
## Training set error measures:
##                      ME     RMSE   MAE        MPE     MAPE      MASE
## Training set 0.06359505 2.204238 1.702 -0.1746113 3.639055 0.6100359
##                     ACF1
## Training set 0.000283503
lmtest::coeftest(tmodel1.ofq)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.310965   0.075473  4.1202 3.785e-05 ***
## ar2   0.083180   0.078463  1.0601   0.28909    
## ar3  -0.038463   0.075315 -0.5107   0.60956    
## sar1 -0.170163   0.488084 -0.3486   0.72736    
## sar2  0.185954   0.194073  0.9582   0.33798    
## sma1 -0.833704   0.491156 -1.6974   0.08961 .  
## sma2  0.033279   0.336027  0.0990   0.92111    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ARIMA(2,0,1)x(2,1,2)12
tmodel2.ofq <- Arima(train.ts,order=c(2,0,1),seasonal=c(2,1,2))
summary(tmodel1.ofq)
## Series: train.ts 
## ARIMA(3,0,0)(2,1,2)[12] 
## 
## Coefficients:
##          ar1     ar2      ar3     sar1    sar2     sma1    sma2
##       0.3110  0.0832  -0.0385  -0.1702  0.1860  -0.8337  0.0333
## s.e.  0.0755  0.0785   0.0753   0.4881  0.1941   0.4912  0.3360
## 
## sigma^2 = 5.392:  log likelihood = -410.86
## AIC=837.73   AICc=838.57   BIC=863.27
## 
## Training set error measures:
##                      ME     RMSE   MAE        MPE     MAPE      MASE
## Training set 0.06359505 2.204238 1.702 -0.1746113 3.639055 0.6100359
##                     ACF1
## Training set 0.000283503
lmtest::coeftest(tmodel2.ofq)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)  
## ar1   0.086344   0.605475  0.1426  0.88660  
## ar2   0.147316   0.204510  0.7203  0.47132  
## ma1   0.222419   0.608122  0.3657  0.71455  
## sar1 -0.167632   0.476025 -0.3521  0.72473  
## sar2  0.189104   0.189321  0.9989  0.31786  
## sma1 -0.835568   0.479292 -1.7433  0.08127 .
## sma2  0.032294   0.328912  0.0982  0.92179  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pada model musiman, ordo yang dilakukan overfit adalah ordo musiman (P, Q).

#ARIMA(2,0,0)x(2,1,3)12
tmodel1.ofQS <- Arima(train.ts,order=c(2,0,0),seasonal=c(2,1,3))
summary(tmodel1.ofQS)
## Series: train.ts 
## ARIMA(2,0,0)(2,1,3)[12] 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2    sma1     sma2     sma3
##       0.2984  0.0526  -1.3709  -0.6918  0.3827  -0.2886  -0.4475
## s.e.  0.0988  0.0784   0.4127   0.2270  0.3507   0.5100   0.4438
## 
## sigma^2 = 5.322:  log likelihood = -410.17
## AIC=836.34   AICc=837.18   BIC=861.88
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 0.06363773 2.189783 1.703117 -0.1759209 3.644613 0.610436
##                     ACF1
## Training set 0.002575289
lmtest::coeftest(tmodel1.ofQS)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.298426   0.098849  3.0190 0.0025360 ** 
## ar2   0.052589   0.078433  0.6705 0.5025369    
## sar1 -1.370868   0.412667 -3.3220 0.0008938 ***
## sar2 -0.691805   0.227000 -3.0476 0.0023068 ** 
## sma1  0.382689   0.350678  1.0913 0.2751488    
## sma2 -0.288641   0.510022 -0.5659 0.5714354    
## sma3 -0.447494   0.443796 -1.0083 0.3132948    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model overfitting yang dicobakan menghasilkan nilai AIC lebih besar tetapi signifikansi parameter bertambah menjadi 3 parameter dari model awal yaitu hanya 1 parameter saja. Oleh karena itu, model yang digunakan adalah model overvtting ARIMA(2,0,0)x(2,1,3)12.

Peramalan

ramalan_sarima = forecast::forecast(tmodel1.ofQS, 48)
ramalan_sarima
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1936       38.37110 35.41455 41.32764 33.84945 42.89275
## Feb 1936       38.11211 35.02672 41.19750 33.39342 42.83080
## Mar 1936       40.87891 37.76524 43.99259 36.11695 45.64087
## Apr 1936       46.60665 43.48826 49.72504 41.83749 51.37581
## May 1936       52.89342 49.77417 56.01267 48.12294 57.66390
## Jun 1936       58.58610 55.46670 61.70550 53.81539 63.35681
## Jul 1936       63.99749 60.87806 67.11692 59.22674 68.76824
## Aug 1936       60.74319 57.62376 63.86263 55.97244 65.51395
## Sep 1936       57.62007 54.50063 60.73950 52.84931 62.39083
## Oct 1936       49.49907 46.37964 52.61850 44.72831 54.26983
## Nov 1936       43.06528 39.94584 46.18471 38.29452 47.83604
## Dec 1936       42.04614 38.92671 45.16557 37.27538 46.81689
## Jan 1937       39.88613 36.76652 43.00574 35.11510 44.65716
## Feb 1937       39.72745 36.60782 42.84708 34.95639 44.49851
## Mar 1937       42.27254 39.15291 45.39217 37.50147 47.04361
## Apr 1937       46.48570 43.36606 49.60533 41.71463 51.25677
## May 1937       51.43645 48.31682 54.55609 46.66539 56.20752
## Jun 1937       59.18381 56.06418 62.30344 54.41274 63.95488
## Jul 1937       63.50609 60.38646 66.62573 58.73502 68.27716
## Aug 1937       62.65458 59.53495 65.77422 57.88352 67.42565
## Sep 1937       57.15117 54.03154 60.27081 52.38011 61.92224
## Oct 1937       49.14542 46.02578 52.26505 44.37435 53.91648
## Nov 1937       43.33886 40.21923 46.45849 38.56780 48.10993
## Dec 1937       38.97740 35.85778 42.09703 34.20635 43.74846
## Jan 1938       38.51719 35.19537 41.83902 33.43690 43.59748
## Feb 1938       39.03844 35.69919 42.37768 33.93150 44.14537
## Mar 1938       41.75686 38.41371 45.10002 36.64395 46.86978
## Apr 1938       47.03595 43.69214 50.37977 41.92204 52.14987
## May 1938       52.67880 49.33487 56.02273 47.56470 57.79290
## Jun 1938       59.08157 55.73762 62.42552 53.96744 64.19570
## Jul 1938       63.82958 60.48562 67.17354 58.71544 68.94372
## Aug 1938       61.54110 58.19715 64.88506 56.42696 66.65524
## Sep 1938       57.46732 54.12336 60.81128 52.35318 62.58146
## Oct 1938       49.48387 46.13991 52.82783 44.36973 54.59801
## Nov 1938       42.98029 39.63634 46.32425 37.86615 48.09443
## Dec 1938       39.93622 36.59227 43.28018 34.82209 45.05036
## Jan 1939       39.50167 36.14217 42.86118 34.36375 44.63959
## Feb 1939       38.94663 35.58574 42.30751 33.80659 44.08666
## Mar 1939       41.53209 38.17089 44.89329 36.39158 46.67260
## Apr 1939       46.37924 43.01799 49.74049 41.23865 51.51983
## May 1939       51.98951 48.62825 55.35077 46.84890 57.13011
## Jun 1939       58.81071 55.44945 62.17198 53.67011 63.95132
## Jul 1939       63.72712 60.36586 67.08838 58.58651 68.86773
## Aug 1939       61.74567 58.38441 65.10693 56.60506 66.88628
## Sep 1939       57.35850 53.99723 60.71976 52.21789 62.49910
## Oct 1939       49.26463 45.90337 52.62590 44.12403 54.40524
## Nov 1939       43.28261 39.92135 46.64387 38.14200 48.42321
## Dec 1939       40.74479 37.38353 44.10605 35.60419 45.88539
autoplot(ramalan_sarima, col="blue")

accuracy(ramalan_sarima,test.ts)
##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  0.06363773 2.189783 1.703117 -0.1759209 3.644613 0.6104360
## Test set     -0.09199005 2.168364 1.666182 -0.1893507 3.587790 0.5971979
##                      ACF1 Theil's U
## Training set  0.002575289        NA
## Test set     -0.071653424 0.4374942