ARIMA Musiman

Library

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

Data 1

Import Data

library(readxl)
data1a<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet1")
data1a<-ts(data1a$Yt,start=c(1903))

data1<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet1")
data1 <- ts(data1,start=c(1903))
ts.plot(data1[,"Yt"], type="l", ylab="Yt")

Berdasarkan hasil plot di atas, terlihat bahwa data Yt tidak memiliki pola musiman, namun data Yt cenderung memiliki pola trend.

Partisi Data

Data Yt dipartisi menjadi dua bagian yaitu data training 80% (1903-1994) dan data testing 20% (1995-2016). Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk validasi model.

train1<-subset(data1[,"Yt"], start = 1, end =92)
head(train1)

## Time Series:
## Start = 1903 
## End = 1908 
## Frequency = 1 
## [1] 221 279 276 111 503 360

test1<-subset(data1[,"Yt"], start = 93, end =114)
head(test1)

## Time Series:
## Start = 1995 
## End = 2000 
## Frequency = 1 
## [1] 2238 2122 1749 1781 1862 1604

Plot ACF dan Porses Differencing Data

Plot ACF Data

acf(train1, lag.max=20)

Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi turun secara perlahan yang mengindikasikan bahwa data tidak stasioner. Karena data tidak stasioner maka selanjutnya dilakukan differencing untuk membuat data stasioner sebelum dilakukan identifikasi model tentatif.

Differencing

# Differencing Ordo 1
train.dif1 <- diff(train1, difference=1)
plot.ts(train.dif1, lty=1, xlab="Year", ylab="Yt")
points(train.dif1)

# Cek Kestasioneran
acf(train.dif1, lag.max=20)

adf.test(train.dif1, alternative = c("stationary"),
         k = trunc((length(train.dif1)-1)^(1/3))) # Augmented Dickey-Fuller

## Warning in adf.test(train.dif1, alternative = c("stationary"), k =
## trunc((length(train.dif1) - : p-value smaller than printed p-value

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

Berdasarkan plot data dan ACF yang tidak turun secara melandai, menunjukkan bahwa data stasioner. Selanjutnya, berdasarkan hasil uji ADF, diperoleh nilai p-value = 0.01 < 0.05 yang berarti \(H_0\) ditolak, sehingga dapat disimpulkan bahwa Data Yt telah stasioner setelah dilakukan differencing ordo 1. Karena data telah stasioner selanjutnya akan dilakukan identifikasi model tentatif.

Identifikasi Model Tentatif

Plot ACF

acf(train.dif1, lag.max=20)

Plot PACF

pacf(train.dif1, lag.max=20)

EACF

eacf(train.dif1)

## 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 o o  o  o  o 
## 1 x x o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x x o o o o o o o o o  o  o  o 
## 6 x x o o o o o o o o o  o  o  o 
## 7 x o x x o o o o o o o  o  o  o

Kandidat Model

Berdasarkan plot ACF (cut off setelah lag ke-1), PACF (cut off setelah lag ke-1), dan EACF, diperoleh kandidat model:

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

Pendugaan Parameter dan Penentuan Model Terbaik

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

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

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

arima011

## 
## Call:
## arima(x = train.dif1, order = c(0, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ma1  intercept
##       -0.7013    19.3922
## s.e.   0.1064     6.0228
## 
## sigma^2 estimated as 35052:  log likelihood = -605.6,  aic = 1215.2

lmtest::coeftest(arima011)

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1       -0.70128    0.10639 -6.5917 4.347e-11 ***
## intercept 19.39223    6.02280  3.2198  0.001283 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima110

## 
## Call:
## arima(x = train.dif1, order = c(1, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ar1  intercept
##       -0.4320    20.9508
## s.e.   0.0943    14.5352
## 
## sigma^2 estimated as 39157:  log likelihood = -610.4,  aic = 1224.81

lmtest::coeftest(arima110)

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1       -0.431960   0.094335 -4.5790 4.672e-06 ***
## intercept 20.950791  14.535170  1.4414    0.1495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima111

## 
## Call:
## arima(x = train.dif1, order = c(1, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.3332  -1.0000    19.5427
## s.e.  0.1007   0.0382     1.0389
## 
## sigma^2 estimated as 32115:  log likelihood = -603.2,  aic = 1212.4

lmtest::coeftest(arima111)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        0.333184   0.100701   3.3086 0.0009375 ***
## ma1       -1.000000   0.038158 -26.2071 < 2.2e-16 ***
## intercept 19.542710   1.038939  18.8103 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aic.arimaA<-data.frame(
  "Model"=c("ARIMA(0,1,1)", "ARIMA(1,1,0)", "ARIMA(1,1,1)"),
  "AIC"=c(arima011$aic, arima110$aic, arima111$aic))

aic.arimaA

##          Model      AIC
## 1 ARIMA(0,1,1) 1215.202
## 2 ARIMA(1,1,0) 1224.810
## 3 ARIMA(1,1,1) 1212.400

Berdasarkan hasil di atas dapat disimpulkan bahwa ARIMA (1,1,1) merupakan model terbaik karena model tersebut memiliki nilai AIC terkecil yaitu 1212.400. Selain itu, apabila dilihat dari hasil uji signifikansi penduga parameter, dengan fungsi coeftest terlihat bahwa seluruh penduga parameter telah signifikan terhadap model.

Diagnostik Model

Diagnostik model dilakukan untuk mengetahui apakah sisaan dari model yang dibangun telah memenuhi asumsi pemodelan. Asumsi sisaan yang harus terpenuhi pada model ARIMA(p,d,q) adalah sebagai berikut:

Sisaan saling bebas (tidak ada autokorelasi) dan identik
Sisaan mengikuti sebaran Normal (0, \(\sigma^2\))

Eksploratif

Berikut dilakukan pemeriksaan asumsi model secara eksploratif melalui Normal Q-Q Plot dan Histogram untuk Normalitas, Plot ACF dan PACF untuk non-autokorelasi.

arima.111y<-arima(train1,   order=c(1,1,1), include.mean = TRUE,method="ML") #ARIMA (1,1,1)
sisaan111 <- arima.111y$residuals

# Eksplorasi
par(mfrow=c(2,2))
qqnorm(sisaan111)
qqline(sisaan111, col = "blue", lwd = 2)
plot(c(1:length(sisaan111)),sisaan111)
acf(sisaan111)
pacf(sisaan111)

#diagnostik model
checkresiduals(arima.111y)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 5.2048, df = 8, p-value = 0.7355
## 
## Model df: 2.   Total lags used: 10

Berdasarkan hasil diatas, secara eksploratif Normal Q-Q Plot menunjukkan bahwa sisaan tidak mengikuti sebaran normal karena banyak titik-titik yang tidak berada di sekitar garis. Selanjutnya, apabila dilihat dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal ini menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal.

Uji Non-Autokorelasi dengan `Ljung-Box test`

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

##  lags statistic df   p-value
##     5  2.756503  5 0.7374639
##    10  5.204847 10 0.8770806
##    15  9.746814 15 0.8353565
##    20 12.944725 20 0.8797430
##    25 20.957567 25 0.6949734
##    30 22.570912 30 0.8324806

#atau

Box.test(sisaan111, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan111
## X-squared = 0.33255, df = 1, p-value = 0.5642

Berdasarkan hasil uji Ljung-Box di atas dapat disimpulkan bahwa \(H_0\) tidak ditolak karena diperoleh p-value > 0.05 yang berarti tidak terdapat gejala autokorelasi pada sisaan dari model ARIMA (1,1,1) atau antar \(e_t\) tidak berkorelasi.

Uji Normalitas dengan `Jarque Bera test`

library(tseries)
jarque.bera.test(residuals(arima.111y))

## 
##  Jarque Bera Test
## 
## data:  residuals(arima.111y)
## X-squared = 15.514, df = 2, p-value = 0.0004278

Berdasarkan hasil uji normalitas di atas dapat disimpulkan bahwa \(H_0\) ditolak karena p-value < 0.05 sehingga sisaan dari model ARIMA (1,1,1) tidak mengikuti sebaran Normal.

Uji Nilai Tengah Sisaan

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

## 
##  One Sample t-test
## 
## data:  sisaan111
## t = 2.2643, df = 91, p-value = 0.02593
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   5.473797 83.724977
## sample estimates:
## mean of x 
##  44.59939

Berdasarkan hasil uji nilai tengah sisaan di atas dapat disimpulkan bahwa \(H_0\) ditolak karena p-value < 0.05 sehingga nilai tengah sisaan tidak sama dengan nol.

Overfitting Model

Overfitting model dilakukan dengan menggunakan:

ARIMA (1,1,2)
ARIMA (2,1,1)

arima.112<-arima(train.dif1,    order=c(1,0,2), include.mean = TRUE,method="ML") #ARIMA (1,1,2)

lmtest::coeftest(arima.112)

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.66912    0.19105  3.5023 0.0004612 ***
## ma1       -1.37569    0.23310 -5.9018 3.597e-09 ***
## ma2        0.37570    0.23106  1.6260 0.1039574    
## intercept 19.56412    1.25814 15.5501 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima.211<-arima(train.dif1,    order=c(2,0,1), include.mean = TRUE,method="ML") #ARIMA (2,1,1)

lmtest::coeftest(arima.211)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        0.284728   0.104733   2.7186  0.006556 ** 
## ar2        0.156935   0.104396   1.5033  0.132770    
## ma1       -0.999991   0.034011 -29.4016 < 2.2e-16 ***
## intercept 19.574237   1.205589  16.2362 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena koefisien penduga parameter pada model ARIMA (1,1,2) dan ARIMA (2,1,1) tidak nyata/tidak signifikan maka model yang dipilih tetap model ARIMA (1,1,1).

Validasi Model ARIMA (1,1,1)

#plot dugaan dengan data asli
dugaan <- fitted(arima.111y)
cbind(train1,dugaan)

## Time Series:
## Start = 1903 
## End = 1994 
## Frequency = 1 
##      train1    dugaan
## 1903    221  220.7790
## 1904    279  228.0713
## 1905    276  255.9685
## 1906    111  262.6980
## 1907    503  193.2543
## 1908    360  338.5124
## 1909    515  341.6347
## 1910    349  419.8851
## 1911    366  383.6755
## 1912    363  377.0500
## 1913    320  371.1043
## 1914    473  348.1597
## 1915    403  406.1836
## 1916    474  401.9845
## 1917    505  434.7460
## 1918    503  465.1319
## 1919    363  480.7407
## 1920    555  426.1936
## 1921    419  487.4653
## 1922    542  453.5006
## 1923    641  495.2318
## 1924    651  559.7055
## 1925    661  597.9458
## 1926    661  624.4844
## 1927    639  639.6150
## 1928    654  638.4551
## 1929    597  645.5134
## 1930    803  623.0662
## 1931    764  706.0991
## 1932    746  728.4864
## 1933    780  735.0028
## 1934    854  755.0467
## 1935    870  799.0936
## 1936    751  829.1327
## 1937    833  791.8616
## 1938    867  812.2814
## 1939    780  836.3535
## 1940    730  809.4234
## 1941    948  774.4564
## 1942   1036  855.3351
## 1943    835  933.8161
## 1944   1011  884.6009
## 1945    886  944.2168
## 1946    931  914.9512
## 1947   1051  923.4526
## 1948   1101  981.2404
## 1949   1220 1032.9212
## 1950   1339 1115.3353
## 1951   1338 1212.9091
## 1952   1403 1264.7007
## 1953    899 1324.6808
## 1954   1148 1127.5596
## 1955   1412 1146.2737
## 1956   1139 1267.2419
## 1957   1270 1202.8716
## 1958   1183 1236.0766
## 1959   1105 1210.5044
## 1960   1267 1163.5784
## 1961   1427 1213.0973
## 1962   1542 1308.3130
## 1963   1459 1409.8751
## 1964    912 1426.8157
## 1965   1266 1190.9914
## 1966    979 1236.6364
## 1967   1196 1118.0741
## 1968   1057 1159.2915
## 1969   1463 1111.1871
## 1970   1280 1273.6684
## 1971   1775 1268.7629
## 1972   1354 1498.8932
## 1973   1156 1421.5342
## 1974   1302 1303.3610
## 1975   1300 1308.8039
## 1976   1845 1305.0736
## 1977   1885 1551.2206
## 1978   1209 1691.1714
## 1979   2001 1463.5662
## 1980   1874 1718.8425
## 1981   1863 1777.9086
## 1982   1924 1812.7146
## 1983   2226 1861.3366
## 1984   2150 2024.8638
## 1985   1953 2073.5885
## 1986   2138 2015.5163
## 1987   2087 2073.8801
## 1988   2070 2077.2182
## 1989   1562 2073.5278
## 1990   1630 1840.6704
## 1991   1690 1756.1745
## 1992   1838 1731.2230
## 1993   1955 1781.5563
## 1994   2209 1858.2382

plot.ts(train1, xlab="Month", ylab="Data")
points(train1)
par(col="red")
lines(dugaan)

par(col="black")

Berdasarkan plot data aktual vs dugaan di atas terlihat bahwa nilai dan pola dugaan yang dihasilkan dari model ARIMA (1,1,1) mendekati nilai aktualnya.

Ukuran keakuratan ramalan

#forecast
(ramalan_arima111<- forecast::forecast(train1,model=arima.111y, h=22))

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1995       2014.030 1765.184 2262.877 1633.453 2394.608
## 1996       2006.009 1732.563 2279.455 1587.809 2424.208
## 1997       2005.679 1711.760 2299.598 1556.169 2455.189
## 1998       2005.665 1692.688 2318.642 1527.008 2484.322
## 1999       2005.665 1674.729 2336.601 1499.542 2511.788
## 2000       2005.665 1657.695 2353.634 1473.491 2537.838
## 2001       2005.665 1641.457 2369.872 1448.658 2562.671
## 2002       2005.665 1625.913 2385.416 1424.885 2586.444
## 2003       2005.665 1610.981 2400.348 1402.048 2609.281
## 2004       2005.665 1596.594 2414.736 1380.044 2631.285
## 2005       2005.665 1582.695 2428.634 1358.789 2652.541
## 2006       2005.665 1569.239 2442.090 1338.209 2673.120
## 2007       2005.665 1556.186 2455.144 1318.246 2693.083
## 2008       2005.665 1543.501 2467.828 1298.846 2712.483
## 2009       2005.665 1531.155 2480.174 1279.965 2731.365
## 2010       2005.665 1519.122 2492.207 1261.562 2749.767
## 2011       2005.665 1507.380 2503.949 1243.604 2767.725
## 2012       2005.665 1495.908 2515.421 1226.059 2785.270
## 2013       2005.665 1484.689 2526.640 1208.901 2802.428
## 2014       2005.665 1473.706 2537.623 1192.104 2819.225
## 2015       2005.665 1462.945 2548.384 1175.647 2835.682
## 2016       2005.665 1452.394 2558.935 1159.511 2851.819

plot(ramalan_arima111)

library(dplyr)
gabungan2<- cbind(test1,ramalan_arima111)
df.gab<- as.data.frame(gabungan2, digits=3)
df.gab %>%
  rename(
    "Data Aktual Yt"=test1,
    Ramalan="ramalan_arima111.Point Forecast",
    "Lo 80"="ramalan_arima111.Lo 80",
    "Hi 80"="ramalan_arima111.Hi 80",
    "Lo 95"="ramalan_arima111.Lo 95",
    "Hi 95"="ramalan_arima111.Hi 95"
        )

##    Data Aktual Yt  Ramalan    Lo 80    Hi 80    Lo 95    Hi 95
## 1            2238 2014.030 1765.184 2262.877 1633.453 2394.608
## 2            2122 2006.009 1732.563 2279.455 1587.809 2424.208
## 3            1749 2005.679 1711.760 2299.598 1556.169 2455.189
## 4            1781 2005.665 1692.688 2318.642 1527.008 2484.322
## 5            1862 2005.665 1674.729 2336.601 1499.542 2511.788
## 6            1604 2005.665 1657.695 2353.634 1473.491 2537.838
## 7            1698 2005.665 1641.457 2369.872 1448.658 2562.671
## 8            1800 2005.665 1625.913 2385.416 1424.885 2586.444
## 9            1599 2005.665 1610.981 2400.348 1402.048 2609.281
## 10           1643 2005.665 1596.594 2414.736 1380.044 2631.285
## 11           1636 2005.665 1582.695 2428.634 1358.789 2652.541
## 12           1884 2005.665 1569.239 2442.090 1338.209 2673.120
## 13           2076 2005.665 1556.186 2455.144 1318.246 2693.083
## 14           2031 2005.665 1543.501 2467.828 1298.846 2712.483
## 15           1808 2005.665 1531.155 2480.174 1279.965 2731.365
## 16           1919 2005.665 1519.122 2492.207 1261.562 2749.767
## 17           1421 2005.665 1507.380 2503.949 1243.604 2767.725
## 18           1680 2005.665 1495.908 2515.421 1226.059 2785.270
## 19           1958 2005.665 1484.689 2526.640 1208.901 2802.428
## 20           1888 2005.665 1473.706 2537.623 1192.104 2819.225
## 21           1323 2005.665 1462.945 2548.384 1175.647 2835.682
## 22           1622 2005.665 1452.394 2558.935 1159.511 2851.819

accuracy(ramalan_arima111, test1) #arima

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set   44.59939 193.1176 138.5356   1.952961 14.60255 0.8942216
## Test set     -217.78845 307.2360 257.3913 -13.906400 15.73451 1.6614138
##                     ACF1 Theil's U
## Training set -0.05915499        NA
## Test set      0.28505257  1.307121

Berdasarkan ukuran keakuratan ramalan, diperoleh nilai MAPE yang tidak berbeda jauh antara model pada data training dan data testing, sehingga dengan nilai MAPE sebesar 15.73% mengindikasikan bahwa hasil ramalan yang diperoleh cukup baik.

Peramalan 6 tahun ke depan

Berikut merupakan hasil ramalan data Yt untuk 6 tahun ke depan:

#forecast
arima111y<-arima(data1a, order = c(1,1,1), include.mean = TRUE, method="ML")
(ramalan_arima111<- forecast::forecast(data1a,model=arima111y, h=6))

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2017       1644.656 1389.338 1899.975 1254.180 2035.133
## 2018       1648.694 1363.936 1933.453 1213.194 2084.195
## 2019       1649.414 1346.579 1952.249 1186.268 2112.561
## 2020       1649.542 1330.949 1968.136 1162.296 2136.789
## 2021       1649.565 1316.173 1982.957 1139.686 2159.444
## 2022       1649.569 1302.045 1997.093 1118.077 2181.061

plot(ramalan_arima111)

Data 2

Import Data

library(readxl)
data2.1<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet2")
data2<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet2")
data2 <- ts(data2,start=c(1979,1), frequency = 12)
ts.plot(data2[,"Xt"], type="l", ylab="Yt")

Berdasarkan hasil plot di atas, terlihat bahwa data bulanan Xt memiliki pola musiman. Sehingga selanjutnya akan dilakukan evaluasi periode musiman.

Evaluasi periode musiman

#Plot tahunan
library(forecast)
seasonplot(data2 [,"Xt"],12,main="Xt", ylab="Xt",year.labels = TRUE, col=rainbow(18))

monthplot(data2 [,"Xt"],ylab="Xt")

boxplot(data2 [,"Xt"]~data2[,"bulan"],ylab="Xt",xlab="Bulan")

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 homogenitas dengan kelompok = tahun

library(car)
fligner.test(`Xt` ~ tahun, data=data2)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Xt by tahun
## Fligner-Killeen:med chi-squared = 6.0862, df = 12, p-value = 0.9117

Berdasarkan hasil uji homogenitas di atas, diperoleh nilai p-value=0.9117 > 0.05 sehingga dapat disimpulkan bahwa data Xt telah homogen.

Partisi Data

Data Xt dipartisi menjadi dua bagian yaitu data training 80% (Januari 1979-Desember 1988) dan data testing 20% (Januari 1989-Juni 1991). Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk validasi model.

#partisi data
train2 <- subset(data2[,"Xt"],start=1,end=120)
head(train2)

##         Jan    Feb    Mar    Apr    May    Jun
## 1979 17.083 19.146 20.341 15.431 22.361 21.511

test2 <- subset(data2[,"Xt"],start=121,end=150)
head(test2)

##          Jan     Feb     Mar     Apr     May     Jun
## 1989 145.554 151.668 148.257 144.962 149.425 156.110

Plot ACF dan Proses Differencing

Plot ACF Xt dan Musiman Xt

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

Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi Xt turun secara perlahan yang mengindikasikan bahwa pola data tidak stasioner.

#ACF Musiman Xt

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

Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi untuk lag 12,24, 36, dan 48 turun secara perlahan yang mengindikasikan bahwa pola data tidak stasioner.

Differencing Xt (Non musiman)

# Plot data differencing non musiman
d1.Xt <- diff(train2)
ts.plot(d1.Xt, type="l", ylab="d1 Xt")

Plot data hasil differencing untuk komponen non musiman di atas menunjukkan bahwa, dengan differencing ordo 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non musiman.

Plot ACF Data Differencing Xt (Non Musiman)

# ACF Data d1.Xt

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

acf1.Xt$lag <- acf1.Xt$lag * 12
acf1.Xt.1 <- as.data.frame(cbind(acf1.Xt$acf,acf1.Xt$lag))
acf1.Xt.2 <- acf1.Xt.1[which(acf1.Xt.1$V2%%12==0),]

barplot(height = acf1.Xt.2$V1, names.arg=acf1.Xt.2$V2, ylab="ACF", xlab="Lag")

Plot ACF untuk data differencing Xt untuk komponen non musiman menunjukkan bahwa komponen non musiman telah stasioner. Namun untuk komponen musiman (autokorelasi lag 12, 24, 36, 48) terlihat masih penunjukkan ketidakstasioneran. Sehingga, selanjutnya dilakukan differencing pada data Xt untuk komponen musiman.

Differencing Xt (Musiman)

#Plot data differencing musiman

D1.Xt <- diff(train2,12)
ts.plot(D1.Xt, type="l", ylab="D1 Xt")

Plot Data Differencing Xt untuk komponen musiman menunjukkan pola yang belum stasioner.

Plot ACF Data Differencing Xt (Musiman)

#Plot ACF D1
acf2.Xt <- acf(D1.Xt,lag.max=48,xaxt="n", main="ACF D1 Xt")

acf2.Xt$lag <- acf2.Xt$lag*12
acf2.Xt.1 <- as.data.frame(cbind(acf2.Xt$acf,acf2.Xt$lag))
acf2.Xt.2 <- acf2.Xt.1[which(acf2.Xt.1$V2%%12==0),]

barplot(height = acf2.Xt.2$V1, names.arg=acf2.Xt.2$V2, ylab="ACF", xlab="Lag")

Plot ACF Differencing Xt untuk komponen musiman menunjukkan bahwa differencing mampu mengatasi ketidakstasioneran dalam rataan pada komponen musiman namun tidak untuk komponen non musimannya. Sehingga, selanjutnya dilakukan differencing lagi untuk data Xt yang telah di differencing pada komponen musimannya (D1.Xt).

Differencing D1.Xt

# Plot data d1D1
d1D1.Xt <- diff(D1.Xt)
ts.plot(d1D1.Xt, type="l", ylab="d1 D1 Xt")

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

adf.test(d1D1.Xt)

## Warning in adf.test(d1D1.Xt): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  d1D1.Xt
## Dickey-Fuller = -6.0411, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan plot data dan ACF d1D1.Xt yang tidak turun secara melandai, menunjukkan bahwa data telah stasioner. Selain itu, berdasarkan adf test diperoleh p-value=0.01 < 0.05 sehingga H0 ditolak yang berarti data d1D1.Xt telah stasioner.

Identifikasi Model (Komponen Non Musiman)

# Correlogram Data d1D1.Xt
acf3.Xt <- acf(d1D1.Xt,lag.max=48,xaxt="n", main="ACF d1D1.Xt")
axis(1, at=0:48/12, labels=0:48)

#PACF
pacf3.Xt <- pacf(d1D1.Xt,lag.max=48,xaxt="n", main="PACF d1D1.Xt")
axis(1, at=0:48/12, labels=0:48)

#EACF
eacf(d1D1.Xt)

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

Berdasarkan plot ACF (cut off setelah lag ke-2), PACF (cut off setelah lag ke-2), dan EACF, diperoleh kandidat model untuk komponen non musiman:

ARIMA (0,1,2)
ARIMA (2,1,0)
ARIMA (2,1,2)
ARIMA (3,1,1)

Identifikasi Model (Komponen Musiman)

# Correlogram Musiman Data d1D1

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

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

Berdasarkan plot ACF (cut off setelah lag 1) dan PACF diperoleh kandidat model untuk komponen non musiman adalah ARIMA (0,1,1)12.

Pendugaan Parameter Model ARIMA Musiman

Model tentatif yang diperoleh berdasarkan hasil identifikasi menggunakan plot ACF, PACF, dan EACF adalah:

ARIMA (0,1,2)x(0,1,1)12
ARIMA (2,1,0)x(0,1,1)12
ARIMA (2,1,2)x(0,1,1)12
ARIMA (3,1,1)x(0,1,1)12

#ARIMA(0,1,2)x(0,1,1)12
model1 <- Arima(train2,order=c(0,1,2),seasonal=c(0,1,1))
summary(model1)

## Series: train2 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     ma2     sma1
##       -0.8570  0.2182  -0.4201
## s.e.   0.0939  0.0851   0.0963
## 
## sigma^2 = 1.335:  log likelihood = -167.3
## AIC=342.6   AICc=342.99   BIC=353.29
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.08552043 1.075466 0.7926163 0.09748527 1.197007 0.06069222
##                     ACF1
## Training set -0.02992964

#ARIMA(2,1,0)x(0,1,1)12
model2 <- Arima(train2,order=c(2,1,0),seasonal=c(0,1,1))
summary(model2)

## Series: train2 
## ARIMA(2,1,0)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     sma1
##       -0.7852  -0.2438  -0.3905
## s.e.   0.0943   0.0937   0.0983
## 
## sigma^2 = 1.392:  log likelihood = -169.29
## AIC=346.59   AICc=346.98   BIC=357.28
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.05500423 1.098231 0.8028162 0.05257322 1.222689 0.06147325
##                     ACF1
## Training set -0.06216435

#ARIMA(2,1,2)x(0,1,1)12
model3 <- Arima(train2,order=c(2,1,2),seasonal=c(0,1,1))
summary(model3)

## Series: train2 
## ARIMA(2,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     ma1      ma2     sma1
##       -1.1697  -0.1758  0.3640  -0.5464  -0.5058
## s.e.   0.1380   0.1363  0.1097   0.0929   0.0943
## 
## sigma^2 = 1.26:  log likelihood = -163.85
## AIC=339.7   AICc=340.54   BIC=355.74
## 
## Training set error measures:
##                     ME     RMSE       MAE       MPE     MAPE       MASE
## Training set 0.1055238 1.034696 0.7605394 0.1129866 1.140696 0.05823603
##                     ACF1
## Training set -0.02840871

#ARIMA(3,1,1)x(0,1,1)12
model4 <- Arima(train2,order=c(3,1,1),seasonal=c(0,1,1))
summary(model4)

## Series: train2 
## ARIMA(3,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     sma1
##       -1.6630  -1.0187  -0.3504  0.9142  -0.5228
## s.e.   0.0998   0.1567   0.0954  0.0557   0.0988
## 
## sigma^2 = 1.287:  log likelihood = -165.11
## AIC=342.22   AICc=343.06   BIC=358.26
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE    MAPE       MASE
## Training set 0.07725587 1.045893 0.7716576 0.07532462 1.14633 0.05908738
##                     ACF1
## Training set -0.07146779

akurasiXt<-data.frame("Model ARIMA"=c("ARIMA (0,1,2)x(0,1,1)12", "ARIMA(2,1,0)x(0,1,1)12", "ARIMA(2,1,2)x(0,1,1)12", "ARIMA (3,1,1)x(0,1,1)12"),
                     "MAD"=c(0.7926163, 0.8028162,0.7605394, 0.7716576),
                     "MSE"=c(1.075466^2, 1.098231^2, 1.034696^2, 1.045893^2),
                     "MAPE"=c(1.197007, 1.222689, 1.140696,1.14633), 
                     "AIC"=c(342.6, 346.59, 339.7, 342.22))

akurasiXt

##               Model.ARIMA       MAD      MSE     MAPE    AIC
## 1 ARIMA (0,1,2)x(0,1,1)12 0.7926163 1.156627 1.197007 342.60
## 2  ARIMA(2,1,0)x(0,1,1)12 0.8028162 1.206111 1.222689 346.59
## 3  ARIMA(2,1,2)x(0,1,1)12 0.7605394 1.070596 1.140696 339.70
## 4 ARIMA (3,1,1)x(0,1,1)12 0.7716576 1.093892 1.146330 342.22

Berdasarkan hasil di atas dapat disimpulkan bahwa ARIMA(2,1,2)x(0,1,1)12 merupakan model terbaik karena model tersebut memiliki ukuran keakuratan ramalan terbaik.

Pengujian Parameter Model

printstatarima.Xt <- 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.Xt(model1)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.8570  0.0939  -9.1267  0.0000
## ma2    0.2182  0.0851   2.5640  0.0116
## sma1  -0.4201  0.0963  -4.3624  0.0000

printstatarima.Xt(model2)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.7852  0.0943  -8.3266  0.0000
## ar2   -0.2438  0.0937  -2.6019  0.0104
## sma1  -0.3905  0.0983  -3.9725  0.0001

printstatarima.Xt(model3)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -1.1697  0.1380  -8.4761  0.0000
## ar2   -0.1758  0.1363  -1.2898  0.1996
## ma1    0.3640  0.1097   3.3181  0.0012
## ma2   -0.5464  0.0929  -5.8816  0.0000
## sma1  -0.5058  0.0943  -5.3637  0.0000

printstatarima.Xt(model4)

## 
## Coefficients:
##                  s.e.         t  sign.
## ar1   -1.6630  0.0998  -16.6633  0e+00
## ar2   -1.0187  0.1567   -6.5010  0e+00
## ar3   -0.3504  0.0954   -3.6730  4e-04
## ma1    0.9142  0.0557   16.4129  0e+00
## sma1  -0.5228  0.0988   -5.2915  0e+00

Berdasarkan hasil pengujian signifikansi parameter model, model ARIMA (0,1,2)x(0,1,1)12, ARIMA (2,1,0)x(0,1,1)12, dan ARIMA (3,1,1)x(0,1,1)12 memiliki seluruh penduga parameter signifikan.

Oleh karena itu, berdasarkan ukuran keakuratan ramalan dan hasil pengujian signifikansi parameter dapat disimpulkan bahwa model terbaik adalah ARIMA (0,1,2)x(0,1,1)12.

Diagnostik Model

Diagnostik model dilakukan untuk mengetahui apakah sisaan dari model yang dibangun telah memenuhi asumsi pemodelan. Asumsi sisaan yang harus terpenuhi pada model ARIMA adalah sebagai berikut:

Sisaan saling bebas (tidak ada autokorelasi) dan identik
Sisaan mengikuti sebaran Normal (0, \(\sigma^2\))

Uji Normalitas dan Autokorelasi

#diagnostik model 1
checkresiduals(model1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 47.84, df = 21, p-value = 0.0007231
## 
## Model df: 3.   Total lags used: 24

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  5.398993  5 0.369148851
##    10 19.634276 10 0.032907871
##    15 28.319784 15 0.019650908
##    20 41.876887 20 0.002869562
##    25 50.534009 25 0.001830393
##    30 55.111368 30 0.003448313

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

## 
##  Jarque Bera Test
## 
## data:  residuals(model1)
## X-squared = 1.7637, df = 2, p-value = 0.414

Secara eksploratif, berdasarkan plot ACF dan PACF terlihat bahwa masih terdapat autokorelasi yang signifikan sehingga terdapat gejala autokorelasi pada sisaan. Karena terdapat gejala autokorelasi selanjutnya dilakukan pengujian menggunakan Ljung-Box Test dan diperoleh hasil p-value = 0.0007231 < \(\alpha\) = 0.05 yang berarti bahwa H0 ditolak (terdapat autokorelasi pada sisaan model).

Selanjutnya, apabila dilihat melalui histogram sisaan, terlihat bahwa sisaaan menyebar normal, selain itu, berdasarkan pengujian menggunakan Jarque bera Test diperoleh p-value = 0.414 > 0.05 yang berarti H0 tidak ditolak sehingga dapat disimpulkan sisaan model menyebar normal.

Model ARIMA (0,1,2)x(0,1,1)12 tidak memenuhi asumsi Non-autokorelasi, sehingga hal tersebut dapat mempengaruhi nilai ramalan yang diperoleh.

Overfitting Model

Selanjutnya dilakukan overfitting model dengan menggunakan:

ARIMA(0,1,2)x(1,1,1)12

#overfitting model
#ARIMA(0,1,2)x(1,1,1)12
model5 <- Arima(train2,order=c(0,1,2),seasonal=c(1,1,1))
summary(model5)

## Series: train2 
## ARIMA(0,1,2)(1,1,1)[12] 
## 
## Coefficients:
##           ma1     ma2     sar1     sma1
##       -0.8574  0.2186  -0.0419  -0.3870
## s.e.   0.0942  0.0854   0.2085   0.1953
## 
## sigma^2 = 1.347:  log likelihood = -167.28
## AIC=344.56   AICc=345.16   BIC=357.93
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.08414932 1.075138 0.7903552 0.09652891 1.194659 0.06051909
##                     ACF1
## Training set -0.02902105

printstatarima.Xt(model5)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.8574  0.0942  -9.1019  0.0000
## ma2    0.2186  0.0854   2.5597  0.0117
## sar1  -0.0419  0.2085  -0.2010  0.8411
## sma1  -0.3870  0.1953  -1.9816  0.0498

Berdasarkan hasil di atas, model ARIMA(0,1,2)x(1,1,1)12 memiliki ukuran keakuratan ramalan yang yang tidak lebih baik dari ARIMA(0,1,2)x(0,1,1)12 serta penduga parameter tambahan tidak nyata/signifikan.

ARIMA(1,1,2)x(0,1,1)12

#overfitting model
#ARIMA(1,1,2)x(0,1,1)12
model6 <- Arima(train2,order=c(1,1,2),seasonal=c(0,1,1))
summary(model6)

## Series: train2 
## ARIMA(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2     sma1
##       -0.9921  0.2700  -0.6104  -0.5137
## s.e.   0.0129  0.0705   0.0671   0.0944
## 
## sigma^2 = 1.265:  log likelihood = -164.66
## AIC=339.32   AICc=339.91   BIC=352.68
## 
## Training set error measures:
##                     ME     RMSE       MAE       MPE     MAPE       MASE
## Training set 0.1128961 1.042073 0.7675287 0.1203949 1.150457 0.05877121
##                   ACF1
## Training set -0.105319

printstatarima.Xt(model6)

## 
## Coefficients:
##                  s.e.         t  sign.
## ar1   -0.9921  0.0129  -76.9070  0e+00
## ma1    0.2700  0.0705    3.8298  2e-04
## ma2   -0.6104  0.0671   -9.0969  0e+00
## sma1  -0.5137  0.0944   -5.4417  0e+00

Berdasarkan hasil di atas, model ARIMA(1,1,2)x(0,1,1)12 memiliki ukuran keakuratan ramalan yang lebih baik dibandingkan dengan model ARIMA(0,1,2)x(0,1,1)12. Selain itu, penduga parameter tambahan nyata/signifikan. Sehingga, dapat disimpulkan bahwa model ARIMA(1,1,2)x(0,1,1)12 lebih sesuai untuk data memodelkan data Xt. Selanjutnya dilakukan pengujian diagnostik untuk mengetahui normalitas serta autokorelasi sisaan model.

#diagnostik model 6
checkresiduals(model6)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,2)(0,1,1)[12]
## Q* = 25.845, df = 20, p-value = 0.171
## 
## Model df: 4.   Total lags used: 24

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

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

##  lags statistic df   p-value
##     5  2.546928  5 0.7694116
##    10 14.231402 10 0.1626959
##    15 20.767656 15 0.1444185
##    20 24.757849 20 0.2108527
##    25 26.302063 25 0.3916225
##    30 29.863937 30 0.4726380

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

## 
##  Jarque Bera Test
## 
## data:  residuals(model6)
## X-squared = 2.8795, df = 2, p-value = 0.237

Secara eksploratif, berdasarkan plot ACF dan PACF terlihat bahwa masih terdapat autokorelasi yang signifikan sehingga terdapat gejala autokorelasi pada sisaan. Karena terdapat gejala autokorelasi selanjutnya dilakukan pengujian menggunakan Ljung-Box Test dan diperoleh hasil p-value = 0.171 > \(\alpha\) = 0.05 yang berarti bahwa H0 tidak ditolak (tidak terdapat autokorelasi pada sisaan model).

Selanjutnya, apabila dilihat melalui histogram sisaan, terlihat bahwa sisaaan menyebar normal, selain itu, berdasarkan pengujian menggunakan Jarque bera Test diperoleh p-value = 0.237 > 0.05 yang berarti H0 tidak ditolak sehingga dapat disimpulkan sisaan model menyebar normal.

Oleh karena itu, dapat disimpulkan bahwa model ARIMA(1,1,2)x(0,1,1)12 lebih sesuai untuk memodelkan data Xt dibandingkan dengan model ARIMA(0,1,2)x(0,1,1)12. Hal tersebut dikarenakan model ARIMA(1,1,2)x(0,1,1)12 memiliki ukuran keakuratan ramalan yang lebih baik, penduga parameter yang nyata/signifikan, dan memenuhi asumsi model ARIMA.

Validasi Model

#validasi model
ramalan_arima = forecast(model6, 30)
accuracy(ramalan_arima,test2)

##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set  0.1128961 1.042073 0.7675287  0.1203949 1.150457 0.05877121
## Test set     -1.7379838 2.197475 1.8058300 -1.0514454 1.096842 0.13827603
##                    ACF1 Theil's U
## Training set -0.1053190        NA
## Test set      0.5271378  0.437112

plot(ramalan_arima)

Berdasarkan ukuran keakuratan ramalan di atas, dapat disimpulkan bahwa model ARIMA(1,1,2)x(0,1,1)12 dapat digunakan untuk meramalkan data Xt dengan baik. Selain itu, apabila dilihat berdasarkan plot ramalan, nilai ramalan yang dihasilkan memiliki pola yang mirip dengan pola data sebelumnya.

forecast_arima <- cbind(ramalan_arima$mean,ramalan_arima$lower,ramalan_arima$upper)
forecast_arima <- cbind(data2.1[121:150,1:2],forecast_arima)
forecast_arima <- as.data.frame(forecast_arima)
plot(forecast_arima[,2],forecast_arima[,1], xlab='Month', ylab='Xt', type='l', col='black', lwd=1, ylim=c(150,200))
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)

Plot di atas menunjukkan garis berwarna merah yang berhimpit dengan garis hitam, yang dapat diartikan bahwa nilai ramalan yang dihasilkan sangat mirip dengan data testing (aktual) yang digunakan.

Peramalan (6 bulan ke depan)

Berikut merupakan hasil ramalan data Xt untuk 6 bulan ke depan:

modelXt<-Arima(data2[,"Xt"],order=c(1,1,2),seasonal=c(0,1,1))
ramalan_arimaXt = forecast(modelXt, 6)
ramalan_arimaXt

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jul 1991       177.3741 175.9664 178.7817 175.2213 179.5269
## Aug 1991       173.2601 171.8003 174.7199 171.0275 175.4926
## Sep 1991       177.8730 176.3219 179.4242 175.5007 180.2454
## Oct 1991       183.9353 182.3363 185.5342 181.4899 186.3807
## Nov 1991       182.0662 180.3838 183.7487 179.4932 184.6393
## Dec 1991       177.5853 175.8585 179.3121 174.9444 180.2263

plot(ramalan_arimaXt)

Data 3

Import Data

library(readxl)
data3.1<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet3")
data3<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/DERET WAKTU/latihan Praktikum 8.xlsx", "Sheet3")
data3 <- ts(data3,start=c(1984,1), frequency = 12)
ts.plot(data3[,"Zt"], type="l", ylab="Zt")

Berdasarkan hasil plot di atas, terlihat bahwa data bulanan Xt memiliki pola musiman. Sehingga selanjutnya akan dilakukan evaluasi periode musiman.

Evaluasi periode musiman

#Plot tahunan
library(forecast)
seasonplot(data3 [,"Zt"],12,main="Zt", ylab="Zt",year.labels = TRUE, col=rainbow(18))

monthplot(data3 [,"Zt"],ylab="Zt")

boxplot(data3 [,"Zt"]~data3[,"bulan"],ylab="Zt",xlab="Bulan")

Berdasarkan hasil plot di atas dapat terlihat bahwa data memiliki pola yang hampir sama dari tahun ke tahun, meskipun ada tahun yang memiliki pola sedikit berbeda, sehingga dapat disimpulkan bahwa data berpola musiman dan periode musimannya adalah 12. Berdasarkan Boxplot terlihat bahwa ragam data dari tahun ke tahun cenderung tidak homogen, sehingga untuk memastikan, selanjutnya dilakukan pengujian dengan fligner test.

Uji Homogenitas Zt

# uji homogenitas dengan kelompok = tahun

library(car)
fligner.test(`Zt` ~ tahun, data=data3)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Zt by tahun
## Fligner-Killeen:med chi-squared = 43.369, df = 24, p-value = 0.009025

Berdasarkan hasil uji homogenitas di atas, diperoleh nilai p-value=0.009025 < 0.05 sehingga dapat disimpulkan bahwa data Zt tidak homogen.

Oleh karena itu, berdasarkan hasil plot data dan pengujian homogenitas dapat disimpulkan bahwa data Zt tidak stasioner dalam ragam dan rataan.

Partisi Data

Data Zt dipartisi menjadi dua bagian yaitu data training 80% (Januari 1984-Desember 2003) dan data testing 20% (Januari 2004-Juni 2008). Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk validasi model.

#partisi data
train3 <- subset(data3[,"Zt"],start=1,end=240)
head(train3)

##      Jan Feb Mar Apr May Jun
## 1984 308 317 316 259 277 308

test3 <- subset(data3[,"Zt"],start=241,end=300)
head(test3)

##      Jan Feb Mar Apr May Jun
## 2004 695 716 729 736 739 698

Plot ACF dan Proses Differencing

Plot ACF Zt dan Musiman Zt

#ACF Zt

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

Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi Zt turun secara perlahan yang mengindikasikan bahwa pola data tidak stasioner.

#ACF Musiman Zt

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

Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi untuk lag 12, 24, 36, dan 48 turun secara perlahan yang mengindikasikan bahwa pola data tidak stasioner.

Differencing Zt (Non musiman)

# Plot data differencing non musiman
d1.Zt <- diff(train3)
ts.plot(d1.Zt, type="l", ylab="d1 Zt")

Plot data hasil differencing untuk komponen non musiman di atas menunjukkan bahwa, dengan differencing ordo 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non musiman.

Plot ACF Data Differencing Zt

# ACF Data d1.Zt

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

acf1.Zt$lag <- acf1.Zt$lag * 12
acf1.Zt.1 <- as.data.frame(cbind(acf1.Zt$acf,acf1.Zt$lag))
acf1.Zt.2 <- acf1.Zt.1[which(acf1.Zt.1$V2%%12==0),]

barplot(height = acf1.Zt.2$V1, names.arg=acf1.Zt.2$V2, ylab="ACF", xlab="Lag")

Plot ACF untuk data differencing Zt untuk komponen non musiman menunjukkan bahwa komponen non musiman telah stasioner. Selain itu, pada komponen musiman (autokorelasi lag 12, 24, 36, 48) juga menunjukkan bahwa data telah stasioner. Sehingga, selanjutnya dapat dilakukan identifikasi model ARIMA Musiman.

Identifikasi Model (Komponen Non Musiman)

# Correlogram Data d1.Zt
acf3.Zt<- acf(d1.Zt,lag.max=48,xaxt="n", main="ACF d1.Zt")
axis(1, at=0:48/12, labels=0:48)

#PACF
pacf3.Zt <- pacf(d1.Zt,lag.max=48,xaxt="n", main="PACF d1.Zt")
axis(1, at=0:48/12, labels=0:48)

#EACF
eacf(d1.Zt)

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

Berdasarkan plot ACF (cut off setelah lag ke-1), PACF (cut off setelah lag ke-3), dan EACF, diperoleh kandidat model untuk komponen non musiman:

ARIMA (0,1,1)
ARIMA (3,1,0)
ARIMA (3,1,1)
ARIMA (1,1,2)

Identifikasi Model (Komponen Musiman)

# Correlogram Data d1.Zt
acf3.Zt<- acf(d1.Zt,lag.max=120,xaxt="n", main="ACF d1.Zt")
axis(1, at=0:120/12, labels=0:120)

#PACF
pacf3.Zt <- pacf(d1.Zt,lag.max=120,xaxt="n", main="PACF d1.Zt")
axis(1, at=0:120/12, labels=0:120)

# Correlogram Musiman Data d1.Zt

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

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

Berdasarkan plot ACF dan PACF diperoleh kandidat model untuk komponen musiman adalah

ARIMA (0,0,2)12
ARIMA (1,0,1)12
ARIMA (1,0,2)12

Pendugaan Parameter Model ARIMA Musiman

Model tentatif yang diperoleh berdasarkan hasil identifikasi menggunakan plot ACF, PACF, dan EACF adalah:

ARIMA (0,1,1)x(0,0,2)12
ARIMA (3,1,0)x(0,0,2)12
ARIMA (3,1,1)x(0,0,2)12
ARIMA (1,1,2)x(0,0,2)12
ARIMA (0,1,1)x(1,0,1)12
ARIMA (3,1,0)x(1,0,1)12
ARIMA (3,1,1)x(1,0,1)12
ARIMA (1,1,2)x(1,0,1)12
ARIMA (0,1,1)x(1,0,2)12
ARIMA (3,1,0)x(1,0,2)12
ARIMA (3,1,1)x(1,0,2)12
ARIMA (1,1,2)x(1,0,2)12

Model untuk komponen musiman ARIMA (0,0,2)12

#ARIMA(0,1,1)x(0,0,2)12
modela <- Arima(train3,order=c(0,1,1),seasonal=c(0,0,2))
summary(modela)

## Series: train3 
## ARIMA(0,1,1)(0,0,2)[12] 
## 
## Coefficients:
##           ma1     sma1     sma2
##       -0.5243  -0.0137  -0.1497
## s.e.   0.0617   0.0686   0.0759
## 
## sigma^2 = 1345:  log likelihood = -1198.99
## AIC=2405.97   AICc=2406.15   BIC=2419.88
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 3.642167 36.37349 27.95632 0.1764506 6.441591 0.5359941 0.0127319

#ARIMA(3,1,0)x(0,0,2)12
modelb <- Arima(train3,order=c(3,1,0),seasonal=c(0,0,2))
summary(modelb)

## Series: train3 
## ARIMA(3,1,0)(0,0,2)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3     sma1     sma2
##       -0.4740  -0.2037  -0.1769  -0.0104  -0.1483
## s.e.   0.0655   0.0696   0.0638   0.0695   0.0785
## 
## sigma^2 = 1355:  log likelihood = -1198.76
## AIC=2409.52   AICc=2409.88   BIC=2430.38
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 3.163084 36.34116 27.93469 0.1002274 6.450431 0.5355793
##                     ACF1
## Training set -0.02186838

#ARIMA(3,1,1)x(0,0,2)12
modelc <- Arima(train3,order=c(3,1,1),seasonal=c(0,0,2))
summary(modelc)

## Series: train3 
## ARIMA(3,1,1)(0,0,2)[12] 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1     sma1     sma2
##       0.0646  0.0399  -0.1087  -0.5628  -0.0133  -0.1553
## s.e.  0.1880  0.1106   0.0791   0.1824   0.0688   0.0763
## 
## sigma^2 = 1342:  log likelihood = -1197.17
## AIC=2408.34   AICc=2408.82   BIC=2432.67
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 4.045868 36.09141 27.76326 0.259694 6.358971 0.5322925 -0.01494826

#ARIMA (1,1,2)x(0,0,2)12
modeld <- Arima(train3,order=c(1,1,2),seasonal=c(0,0,2))
summary(modeld)

## Series: train3 
## ARIMA(1,1,2)(0,0,2)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sma1     sma2
##       0.4102  -0.9053  0.1601  -0.0140  -0.1570
## s.e.  0.4196   0.4241  0.2370   0.0687   0.0753
## 
## sigma^2 = 1351:  log likelihood = -1198.46
## AIC=2408.92   AICc=2409.28   BIC=2429.78
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE      MASE        ACF1
## Training set 4.085148 36.28828 27.84618 0.2609174 6.38023 0.5338823 -0.02277293

# pengujian penduga parameter model

printstatarima.Zt <- 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.Zt(modela)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.5243  0.0617  -8.4976  0.0000
## sma1  -0.0137  0.0686  -0.1997  0.8419
## sma2  -0.1497  0.0759  -1.9723  0.0497

printstatarima.Zt(modelb)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4740  0.0655  -7.2366  0.0000
## ar2   -0.2037  0.0696  -2.9267  0.0038
## ar3   -0.1769  0.0638  -2.7727  0.0060
## sma1  -0.0104  0.0695  -0.1496  0.8812
## sma2  -0.1483  0.0785  -1.8892  0.0601

printstatarima.Zt(modelc)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.0646  0.1880   0.3436  0.7314
## ar2    0.0399  0.1106   0.3608  0.7186
## ar3   -0.1087  0.0791  -1.3742  0.1707
## ma1   -0.5628  0.1824  -3.0855  0.0023
## sma1  -0.0133  0.0688  -0.1933  0.8469
## sma2  -0.1553  0.0763  -2.0354  0.0429

printstatarima.Zt(modeld)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.4102  0.4196   0.9776  0.3293
## ma1   -0.9053  0.4241  -2.1346  0.0338
## ma2    0.1601  0.2370   0.6755  0.5000
## sma1  -0.0140  0.0687  -0.2038  0.8387
## sma2  -0.1570  0.0753  -2.0850  0.0381

akurasiA<-data.frame("Model ARIMA"=c("ARIMA (0,1,1)x(0,0,2)12", "ARIMA(3,1,0)x(0,0,2)12", "ARIMA(3,1,1)x(0,0,2)12", "ARIMA (1,1,2)x(0,0,2)12"),
                     "MAD"=c(27.95632,27.93469, 27.76326,27.84618 ),
                     "MSE"=c(36.37349^2,36.34116^2, 36.09141^2, 36.28828^2),
                     "MAPE"=c(6.441591, 6.450431, 6.358971, 6.38023), 
                     "AIC"=c(2405.97,2409.52,2408.34,2408.92))

akurasiA

##               Model.ARIMA      MAD      MSE     MAPE     AIC
## 1 ARIMA (0,1,1)x(0,0,2)12 27.95632 1323.031 6.441591 2405.97
## 2  ARIMA(3,1,0)x(0,0,2)12 27.93469 1320.680 6.450431 2409.52
## 3  ARIMA(3,1,1)x(0,0,2)12 27.76326 1302.590 6.358971 2408.34
## 4 ARIMA (1,1,2)x(0,0,2)12 27.84618 1316.839 6.380230 2408.92

Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (0,0,2)12 model terbaik adalah ARIMA (0,1,1)x(0,0,2)12.

Model untuk komponen musiman ARIMA (1,0,1)12

#ARIMA(0,1,1)x(1,0,1)12
modele <- Arima(train3,order=c(0,1,1),seasonal=c(1,0,1))
summary(modele)

## Series: train3 
## ARIMA(0,1,1)(1,0,1)[12] 
## 
## Coefficients:
##           ma1     sar1    sma1
##       -0.5501  -0.9102  0.9850
## s.e.   0.0585   0.1786  0.3742
## 
## sigma^2 = 1328:  log likelihood = -1199.58
## AIC=2407.15   AICc=2407.32   BIC=2421.06
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 3.110113 36.13539 27.77592 0.08589515 6.414994 0.5325353 0.0108926

#ARIMA(3,1,0)x(1,0,1)12
modelf<- Arima(train3,order=c(3,1,0),seasonal=c(1,0,1))
summary(modelf)

## Series: train3 
## ARIMA(3,1,0)(1,0,1)[12] 
## 
## Coefficients:
##           ar1      ar2     ar3     sar1    sma1
##       -0.5014  -0.2263  -0.195  -0.9110  1.0000
## s.e.   0.0634   0.0702   0.064   0.0507  0.2885
## 
## sigma^2 = 1317:  log likelihood = -1199.1
## AIC=2410.21   AICc=2410.57   BIC=2431.07
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE    MASE        ACF1
## Training set 2.640456 35.83569 27.55397 0.0141046 6.376895 0.52828 -0.01774199

#ARIMA(3,1,1)x(1,0,1)12
modelg <- Arima(train3,order=c(3,1,1),seasonal=c(1,0,1))
summary(modelg)

## Series: train3 
## ARIMA(3,1,1)(1,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1     sar1    sma1
##       0.0263  0.0302  -0.1138  -0.5536  -0.9044  0.9757
## s.e.  0.2042  0.1231   0.0832   0.2001   0.1787  0.2573
## 
## sigma^2 = 1332:  log likelihood = -1197.85
## AIC=2409.71   AICc=2410.19   BIC=2434.04
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 3.339513 35.95627 27.7045 0.1356012 6.364486 0.5311661
##                      ACF1
## Training set -0.009177149

#ARIMA (1,1,2)x(1,0,1)12
modelh <- Arima(train3,order=c(1,1,2),seasonal=c(1,0,1))
summary(modelh)

## Series: train3 
## ARIMA(1,1,2)(1,0,1)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2     sar1    sma1
##       -0.8182  0.2539  -0.4010  -0.9143  0.9829
## s.e.   0.1315  0.1463   0.0998   0.1784  0.3234
## 
## sigma^2 = 1336:  log likelihood = -1198.97
## AIC=2409.93   AICc=2410.29   BIC=2430.79
## 
## Training set error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 2.985047 36.0885 27.87163 0.0660108 6.453911 0.5343703 0.03519108

# pengujian penduga parameter model

printstatarima.Zt <- 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.Zt(modele)

## 
## Coefficients:
##                  s.e.        t  sign.
## ma1   -0.5501  0.0585  -9.4034  0.000
## sar1  -0.9102  0.1786  -5.0963  0.000
## sma1   0.9850  0.3742   2.6323  0.009

printstatarima.Zt(modelf)

## 
## Coefficients:
##                  s.e.         t   sign.
## ar1   -0.5014  0.0634   -7.9085  0.0000
## ar2   -0.2263  0.0702   -3.2236  0.0014
## ar3   -0.1950  0.0640   -3.0469  0.0026
## sar1  -0.9110  0.0507  -17.9684  0.0000
## sma1   1.0000  0.2885    3.4662  0.0006

printstatarima.Zt(modelg)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.0263  0.2042   0.1288  0.8976
## ar2    0.0302  0.1231   0.2453  0.8064
## ar3   -0.1138  0.0832  -1.3678  0.1727
## ma1   -0.5536  0.2001  -2.7666  0.0061
## sar1  -0.9044  0.1787  -5.0610  0.0000
## sma1   0.9757  0.2573   3.7921  0.0002

printstatarima.Zt(modelh)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.8182  0.1315  -6.2221  0.0000
## ma1    0.2539  0.1463   1.7355  0.0839
## ma2   -0.4010  0.0998  -4.0180  0.0001
## sar1  -0.9143  0.1784  -5.1250  0.0000
## sma1   0.9829  0.3234   3.0393  0.0026

akurasiB<-data.frame("Model ARIMA"=c("ARIMA (0,1,1)x(1,0,1)12","ARIMA(3,1,0)x(1,0,1)12","ARIMA(3,1,1)x(1,0,1)12","ARIMA (1,1,2)x(1,0,1)12"),
                     "MAD"=c(27.77592, 27.55397, 27.7045, 27.87163),
                     "MSE"=c(36.13539^2, 35.83569^2, 35.95627^2, 36.0885^2),
                     "MAPE"=c(6.414994, 6.376895, 6.364486, 6.453911 ), 
                     "AIC"=c(2407.15, 2410.21, 2409.71, 2409.93))

akurasiB

##               Model.ARIMA      MAD      MSE     MAPE     AIC
## 1 ARIMA (0,1,1)x(1,0,1)12 27.77592 1305.766 6.414994 2407.15
## 2  ARIMA(3,1,0)x(1,0,1)12 27.55397 1284.197 6.376895 2410.21
## 3  ARIMA(3,1,1)x(1,0,1)12 27.70450 1292.853 6.364486 2409.71
## 4 ARIMA (1,1,2)x(1,0,1)12 27.87163 1302.380 6.453911 2409.93

Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (1,0,1)12 model terbaik adalah ARIMA (3,1,0)x(1,0,1)12.

Model untuk komponen musiman ARIMA (1,0,2)12

#ARIMA(0,1,1)x(1,0,2)12
modelg <- Arima(train3,order=c(0,1,1),seasonal=c(1,0,2))
summary(modelg)

## Series: train3 
## ARIMA(0,1,1)(1,0,2)[12] 
## 
## Coefficients:
##           ma1     sar1    sma1     sma2
##       -0.5417  -0.7261  0.7329  -0.0962
## s.e.   0.0618   0.4411  0.4628   0.1061
## 
## sigma^2 = 1350:  log likelihood = -1199.13
## AIC=2408.26   AICc=2408.52   BIC=2425.65
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 3.344232 36.35256 27.88188 0.1254645 6.443937 0.5345669 0.01362992

#ARIMA(3,1,0)x(1,0,2)12
modelh<- Arima(train3,order=c(3,1,0),seasonal=c(1,0,2))
summary(modelh)

## Series: train3 
## ARIMA(3,1,0)(1,0,2)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3    sar1     sma1     sma2
##       -0.4800  -0.2031  -0.1879  0.5839  -0.5954  -0.1278
## s.e.   0.0651   0.0697   0.0636  0.2521   0.2558   0.0778
## 
## sigma^2 = 1345:  log likelihood = -1197.64
## AIC=2409.27   AICc=2409.76   BIC=2433.61
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 3.89497 36.13322 27.99569 0.2258393 6.445595 0.5367488 -0.02492367

#ARIMA(3,1,1)x(1,0,2)12
modeli <- Arima(train3,order=c(3,1,1),seasonal=c(1,0,2))
summary(modeli)

## Series: train3 
## ARIMA(3,1,1)(1,0,2)[12] 
## 
## Coefficients:
##           ar1     ar2      ar3      ma1    sar1     sma1     sma2
##       -0.0026  0.0160  -0.1280  -0.4986  0.5200  -0.5347  -0.1305
## s.e.   0.2135  0.1204   0.0799   0.2105  0.3102   0.3136   0.0803
## 
## sigma^2 = 1338:  log likelihood = -1196.5
## AIC=2409.01   AICc=2409.63   BIC=2436.82
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 4.562606 35.96971 27.79506 0.3463615 6.362438 0.5329022
##                     ACF1
## Training set -0.01819143

#ARIMA (1,1,2)x(1,0,2)12
modelj<- Arima(train3,order=c(1,1,2),seasonal=c(1,0,2))
summary(modelj)

## Series: train3 
## ARIMA(1,1,2)(1,0,2)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2     sar1    sma1     sma2
##       -0.8184  0.2618  -0.4021  -0.7777  0.7901  -0.0773
## s.e.   0.1442  0.1590   0.1047   0.3544  0.3770   0.1004
## 
## sigma^2 = 1355:  log likelihood = -1198.67
## AIC=2411.34   AICc=2411.82   BIC=2435.67
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE       ACF1
## Training set 3.186689 36.27613 27.9182 0.0996393 6.471266 0.5352632 0.03406512

# pengujian penduga parameter model
printstatarima.Zt <- 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.Zt(modelg)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.5417  0.0618  -8.7654  0.0000
## sar1  -0.7261  0.4411  -1.6461  0.1011
## sma1   0.7329  0.4628   1.5836  0.1146
## sma2  -0.0962  0.1061  -0.9067  0.3655

printstatarima.Zt(modelh)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4800  0.0651  -7.3733  0.0000
## ar2   -0.2031  0.0697  -2.9139  0.0039
## ar3   -0.1879  0.0636  -2.9544  0.0034
## sar1   0.5839  0.2521   2.3161  0.0214
## sma1  -0.5954  0.2558  -2.3276  0.0208
## sma2  -0.1278  0.0778  -1.6427  0.1018

printstatarima.Zt(modeli)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.0026  0.2135  -0.0122  0.9903
## ar2    0.0160  0.1204   0.1329  0.8944
## ar3   -0.1280  0.0799  -1.6020  0.1105
## ma1   -0.4986  0.2105  -2.3686  0.0186
## sar1   0.5200  0.3102   1.6763  0.0950
## sma1  -0.5347  0.3136  -1.7050  0.0895
## sma2  -0.1305  0.0803  -1.6252  0.1054

printstatarima.Zt(modelj)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.8184  0.1442  -5.6755  0.0000
## ma1    0.2618  0.1590   1.6465  0.1010
## ma2   -0.4021  0.1047  -3.8405  0.0002
## sar1  -0.7777  0.3544  -2.1944  0.0292
## sma1   0.7901  0.3770   2.0958  0.0372
## sma2  -0.0773  0.1004  -0.7699  0.4421

akurasiC<-data.frame("Model ARIMA"=c("ARIMA (0,1,1)x(1,0,2)12","ARIMA(3,1,0)x(1,0,2)12","ARIMA(3,1,1)x(1,0,2)12","ARIMA (1,1,2)x(1,0,2)12"),
                     "MAD"=c(27.88188, 27.99569, 27.79506, 27.9182 ),
                     "MSE"=c(36.35256^2, 36.13322^2, 35.96971^2, 36.27613^2),
                     "MAPE"=c(6.443937, 6.445595, 6.362438, 6.471266), 
                     "AIC"=c(2408.26, 2409.27, 2409.01, 2411.34))

akurasiC

##               Model.ARIMA      MAD      MSE     MAPE     AIC
## 1 ARIMA (0,1,1)x(1,0,2)12 27.88188 1321.509 6.443937 2408.26
## 2  ARIMA(3,1,0)x(1,0,2)12 27.99569 1305.610 6.445595 2409.27
## 3  ARIMA(3,1,1)x(1,0,2)12 27.79506 1293.820 6.362438 2409.01
## 4 ARIMA (1,1,2)x(1,0,2)12 27.91820 1315.958 6.471266 2411.34

Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (1,0,2)12 model terbaik adalah ARIMA (3,1,0)x(1,0,2)12.

Pemilihan Model Terbaik

printstatarima.Zt(modela)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.5243  0.0617  -8.4976  0.0000
## sma1  -0.0137  0.0686  -0.1997  0.8419
## sma2  -0.1497  0.0759  -1.9723  0.0497

printstatarima.Zt(modelf)

## 
## Coefficients:
##                  s.e.         t   sign.
## ar1   -0.5014  0.0634   -7.9085  0.0000
## ar2   -0.2263  0.0702   -3.2236  0.0014
## ar3   -0.1950  0.0640   -3.0469  0.0026
## sar1  -0.9110  0.0507  -17.9684  0.0000
## sma1   1.0000  0.2885    3.4662  0.0006

printstatarima.Zt(modelh)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4800  0.0651  -7.3733  0.0000
## ar2   -0.2031  0.0697  -2.9139  0.0039
## ar3   -0.1879  0.0636  -2.9544  0.0034
## sar1   0.5839  0.2521   2.3161  0.0214
## sma1  -0.5954  0.2558  -2.3276  0.0208
## sma2  -0.1278  0.0778  -1.6427  0.1018

akurasiD<-data.frame("Model ARIMA"=c("ARIMA(0,1,1)(0,0,2)12", "ARIMA(3,1,0)(1,0,1)12", "ARIMA(3,1,0)(1,0,2)12"),
                     "MAD"=c(27.95632, 27.55397, 27.99569),
                     "MSE"=c(36.37349^2, 35.83569^2, 36.13322^2),
                     "MAPE"=c(6.441591, 6.376895, 6.445595), 
                     "AIC"=c(2405.97, 2410.21, 2409.27))

akurasiD

##             Model.ARIMA      MAD      MSE     MAPE     AIC
## 1 ARIMA(0,1,1)(0,0,2)12 27.95632 1323.031 6.441591 2405.97
## 2 ARIMA(3,1,0)(1,0,1)12 27.55397 1284.197 6.376895 2410.21
## 3 ARIMA(3,1,0)(1,0,2)12 27.99569 1305.610 6.445595 2409.27

Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dari tiga model terbaik dari masing-masing model musiman, model terbaik untuk data Zt adalah ARIMA(3,1,0)(1,0,1)12.

Diagnostik Model

#diagnostik model
checkresiduals(modelf)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,0)(1,0,1)[12]
## Q* = 26.466, df = 19, p-value = 0.1177
## 
## Model df: 5.   Total lags used: 24

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

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

##  lags statistic df   p-value
##     5  4.487812  5 0.4815158
##    10  8.347219 10 0.5949611
##    15 14.704577 15 0.4728990
##    20 21.645705 20 0.3600409
##    25 27.039863 25 0.3538960
##    30 38.589987 30 0.1351824

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

## 
##  Jarque Bera Test
## 
## data:  residuals(modelf)
## X-squared = 0.9439, df = 2, p-value = 0.6238

Secara eksploratif, berdasarkan plot ACF dan PACF terlihat bahwa tidak terdapat autokorelasi yang signifikan sehingga tidak ada gejala autokorelasi pada sisaan. Selanjutnya dilakukan pengujian menggunakan Ljung-Box Test dan diperoleh hasil p-value = 0.1177 > \(\alpha\) = 0.05 yang berarti bahwa H0 tidak ditolak (tidak terdapat autokorelasi pada sisaan model).

Selanjutnya, apabila dilihat melalui histogram sisaan, terlihat bahwa sisaaan menyebar normal, selain itu, berdasarkan pengujian menggunakan Jarque bera Test diperoleh p-value = 0.6238 > 0.05 yang berarti bahwa H0 tidak ditolak sehingga dapat disimpulkan sisaan model menyebar normal.

Validasi Model

#validasi model
ramalan_arima3 = forecast(modelf, 60)
ramalan_arima3

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2004       652.2388 605.1608 699.3168 580.2392 724.2384
## Feb 2004       649.7381 597.1349 702.3412 569.2885 730.1876
## Mar 2004       654.7100 596.6177 712.8023 565.8655 743.5545
## Apr 2004       650.0253 588.5145 711.5362 555.9527 744.0980
## May 2004       650.9889 583.9032 718.0746 548.3902 753.5876
## Jun 2004       651.7837 580.4459 723.1215 542.6819 760.8855
## Jul 2004       661.5239 585.9648 737.0831 545.9663 777.0816
## Aug 2004       659.6233 580.3830 738.8636 538.4357 780.8109
## Sep 2004       649.2832 566.2729 732.2934 522.3300 776.2364
## Oct 2004       649.6693 563.1440 736.1945 517.3403 781.9982
## Nov 2004       658.1003 568.1521 748.0485 520.5364 795.6642
## Dec 2004       662.5543 569.3556 755.7529 520.0193 805.0892
## Jan 2005       654.6522 557.4652 751.8392 506.0175 803.2869
## Feb 2005       656.2081 555.6049 756.8113 502.3488 810.0674
## Mar 2005       655.2399 551.3049 759.1749 496.2851 814.1948
## Apr 2005       656.1865 549.1003 763.2728 492.4123 819.9608
## May 2005       656.3089 546.0565 766.5614 487.6924 824.9255
## Jun 2005       655.1406 541.8503 768.4309 481.8780 828.4031
## Jul 2005       646.9115 530.6519 763.1711 469.1078 824.7152
## Aug 2005       648.2256 529.0864 767.3648 466.0179 830.4333
## Sep 2005       657.7952 535.8318 779.7586 471.2683 844.3221
## Oct 2005       657.3371 532.6182 782.0560 466.5960 848.0782
## Nov 2005       649.7575 522.3402 777.1748 454.8895 844.6254
## Dec 2005       645.6443 515.5870 775.7017 446.7388 844.5498
## Jan 2006       652.8687 520.7423 784.9951 450.7988 854.9386
## Feb 2006       651.4313 517.0193 785.8434 445.8658 856.9969
## Mar 2006       652.3284 515.6812 788.9756 443.3445 861.3122
## Apr 2006       651.4580 512.5648 790.3511 439.0393 863.8767
## May 2006       651.3510 510.3160 792.3861 435.6565 867.0455
## Jun 2006       652.4120 509.2413 795.5826 433.4514 871.3726
## Jul 2006       659.9106 514.6423 805.1790 437.7419 882.0794
## Aug 2006       658.7123 511.3658 806.0587 433.3653 884.0593
## Sep 2006       649.9954 500.6084 799.3825 421.5277 878.4632
## Oct 2006       650.4122 499.0086 801.8159 418.8603 881.9642
## Nov 2006       657.3174 503.9254 810.7093 422.7246 891.9101
## Dec 2006       661.0641 505.7076 816.4206 423.4669 898.6613
## Jan 2007       654.4830 496.7725 812.1936 413.2856 895.6805
## Feb 2007       655.7924 495.9703 815.6145 411.3655 900.2192
## Mar 2007       654.9752 493.0591 816.8914 407.3459 902.6046
## Apr 2007       655.7681 491.8224 819.7139 405.0347 906.5016
## May 2007       655.8656 489.8596 821.8716 401.9813 909.7498
## Jun 2007       654.8991 486.8792 822.9190 397.9348 911.8633
## Jul 2007       648.0680 478.0528 818.0832 388.0522 908.0839
## Aug 2007       649.1597 477.1807 821.1386 386.1406 912.1788
## Sep 2007       657.1005 483.1730 831.0280 391.1013 923.0996
## Oct 2007       656.7208 480.8689 832.5726 387.7786 925.6629
## Nov 2007       650.4304 472.6737 828.1871 378.5750 922.2858
## Dec 2007       647.0173 467.3773 826.6573 372.2817 921.7529
## Jan 2008       653.0124 471.8240 834.2008 375.9086 930.1162
## Feb 2008       651.8196 468.9431 834.6962 372.1341 931.5052
## Mar 2008       652.5640 468.0225 837.1056 370.3320 934.7960
## Apr 2008       651.8417 465.6209 838.0625 367.0416 936.6418
## May 2008       651.7530 463.9107 839.5952 364.4730 939.0329
## Jun 2008       652.6334 463.1673 842.0995 362.8700 942.3968
## Jul 2008       658.8563 467.7843 849.9282 366.6369 951.0756
## Aug 2008       657.8618 465.1908 850.5328 363.1969 952.5267
## Sep 2008       650.6280 456.3766 844.8795 353.5461 947.7099
## Oct 2008       650.9739 455.1526 846.7952 351.4910 950.4568
## Nov 2008       656.7042 459.3267 854.0817 354.8413 958.5671
## Dec 2008       659.8135 460.8908 858.7362 355.5875 964.0395

accuracy(ramalan_arima3,test3)

##                     ME      RMSE      MAE       MPE      MAPE     MASE
## Training set  2.640456  35.83569 27.55397 0.0141046  6.376895 0.528280
## Test set     51.355837 106.74260 89.79604 5.5689630 12.295929 1.721619
##                     ACF1 Theil's U
## Training set -0.01774199        NA
## Test set      0.88812529  2.271099

plot(ramalan_arima3)

Berdasarkan ukuran keakuratan ramalan di atas, dapat disimpulkan bahwa model ARIMA(3,1,0)x(1,0,1)12 dapat digunakan untuk meramalkan data Zt dengan baik. Selain itu, apabila dilihat berdasarkan plot ramalan, nilai ramalan yang dihasilkan memiliki nilai yang mirip dengan nilai data sebelumnya.

Peramalan (6 bulan ke depan)

Berikut merupakan hasil ramalan data Zt untuk 6 bulan ke depan:

modelZt<-Arima(data3[,"Zt"],order=c(3,1,0),seasonal=c(1,0,1))
ramalan_arimaZt = forecast(modelZt, 6)
ramalan_arimaZt

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2009       604.6808 556.5817 652.7799 531.1195 678.2421
## Feb 2009       607.0078 553.0057 661.0098 524.4187 689.5968
## Mar 2009       607.4169 547.5811 667.2527 515.9059 698.9279
## Apr 2009       608.4863 544.3451 672.6276 510.3907 706.5819
## May 2009       607.3433 537.5036 677.1831 500.5327 714.1540
## Jun 2009       610.5337 536.0789 684.9885 496.6649 724.4025

plot(ramalan_arimaZt)

ARIMA Musiman

ARIMA Musiman

Data 1

Import Data

Partisi Data

Plot ACF dan Porses Differencing Data

Plot ACF Data

Differencing

Identifikasi Model Tentatif

Plot ACF

Plot PACF

EACF

Kandidat Model

Pendugaan Parameter dan Penentuan Model Terbaik

Diagnostik Model

Eksploratif

Uji Non-Autokorelasi dengan Ljung-Box test

Uji Normalitas dengan Jarque Bera test

Uji Nilai Tengah Sisaan

Overfitting Model

Validasi Model ARIMA (1,1,1)

Ukuran keakuratan ramalan

Peramalan 6 tahun ke depan

Data 2

Import Data

Evaluasi periode musiman

Uji Homogenitas

Partisi Data

Plot ACF dan Proses Differencing

Plot ACF Xt dan Musiman Xt

Differencing Xt (Non musiman)

Plot ACF Data Differencing Xt (Non Musiman)

Differencing Xt (Musiman)

Plot ACF Data Differencing Xt (Musiman)

Differencing D1.Xt

Identifikasi Model (Komponen Non Musiman)

Identifikasi Model (Komponen Musiman)

Pendugaan Parameter Model ARIMA Musiman

Pengujian Parameter Model

Diagnostik Model

Uji Normalitas dan Autokorelasi

Overfitting Model

Validasi Model

Peramalan (6 bulan ke depan)

Data 3

Import Data

Evaluasi periode musiman

Uji Homogenitas Zt

Partisi Data

Plot ACF dan Proses Differencing

Plot ACF Zt dan Musiman Zt

Differencing Zt (Non musiman)

Plot ACF Data Differencing Zt

Identifikasi Model (Komponen Non Musiman)

Identifikasi Model (Komponen Musiman)

Pendugaan Parameter Model ARIMA Musiman

Model untuk komponen musiman ARIMA (0,0,2)12

Model untuk komponen musiman ARIMA (1,0,1)12

Model untuk komponen musiman ARIMA (1,0,2)12

Pemilihan Model Terbaik

Diagnostik Model

Validasi Model

Peramalan (6 bulan ke depan)

Uji Non-Autokorelasi dengan `Ljung-Box test`

Uji Normalitas dengan `Jarque Bera test`