Tugas 67 MPDW

Packages

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(tsibble)

## Warning: package 'tsibble' was built under R version 4.2.3

## 
## Attaching package: 'tsibble'

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(MASS)

## Warning: package 'MASS' was built under R version 4.2.3

library(forecast)

## Warning: package 'forecast' was built under R version 4.2.3

library(TSA)

## Warning: package 'TSA' was built under R version 4.2.3

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(TTR)

## Warning: package 'TTR' was built under R version 4.2.3

library(aTSA)

## 
## Attaching package: 'aTSA'

## The following object is masked from 'package:forecast':
## 
##     forecast

## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test

## The following object is masked from 'package:graphics':
## 
##     identify

library(graphics)

Impor Data

datalat67 <- rio::import("https://raw.githubusercontent.com/afhkaniase/praktikum-mpdw/main/Data/Data%20Temperatur%20Portugal%202023.csv")
str(datalat67)

## 'data.frame':    212 obs. of  2 variables:
##  $ Periode    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Temperature: num  15.9 11.1 10.2 9.9 9.9 11 14.8 16.2 15.4 14.4 ...

head(datalat67)

##   Periode Temperature
## 1       1        15.9
## 2       2        11.1
## 3       3        10.2
## 4       4         9.9
## 5       5         9.9
## 6       6        11.0

Mengubah data agar terbaca sebagai data deret waktu dengan fungsi `ts()` .

datalat67.ts <- ts(datalat67$Temperature)
datalat67.ts

## Time Series:
## Start = 1 
## End = 212 
## Frequency = 1 
##   [1] 15.9 11.1 10.2  9.9  9.9 11.0 14.8 16.2 15.4 14.4 13.9 11.5 11.6 10.5 11.6
##  [16] 12.9 12.8  9.4  9.9 14.0 13.7 10.2  8.2  8.3  7.6  8.8  9.3  8.3  7.5  7.5
##  [31]  8.5  8.3  8.3 10.0 11.6 11.7 11.0  9.7  9.5 12.1 10.4  9.6 11.4 12.3 12.3
##  [46] 14.4 13.7 13.9 14.4 13.9 13.8 13.5 12.5  9.8  7.9  8.8 10.6  9.9  8.5  8.1
##  [61]  8.2  9.5 10.5 11.7 13.3 14.2 15.8 16.5 15.4 16.9 15.2 15.2 14.4 14.3 13.2
##  [76] 14.7 13.7 14.5 14.7 14.3 13.9 14.1 16.0 14.2 15.6 15.7 17.4 17.0 16.8 17.0
##  [91] 15.8 16.3 15.1 17.7 19.0 19.0 17.5 17.1 18.2 16.1 16.7 15.3 13.2 14.7 18.1
## [106] 20.2 22.9 21.4 19.2 16.2 14.4 15.8 16.7 18.7 18.6 20.2 23.7 21.6 18.8 19.0
## [121] 21.0 21.7 20.4 18.7 18.2 18.9 19.2 19.4 19.5 17.4 17.6 16.7 17.7 18.1 17.4
## [136] 19.8 21.6 20.0 19.5 19.3 18.2 18.1 19.2 20.2 21.8 17.7 18.8 18.5 18.2 19.0
## [151] 18.9 19.4 20.0 20.1 20.6 21.6 22.0 21.4 20.0 20.0 20.9 20.1 20.6 20.0 21.0
## [166] 24.0 24.4 22.1 21.3 20.5 20.4 19.9 21.7 26.0 27.4 25.1 24.7 23.7 23.6 21.5
## [181] 21.9 23.7 24.6 22.2 21.5 21.4 21.5 21.1 21.9 21.5 23.3 22.5 22.8 23.9 22.7
## [196] 21.7 20.3 22.0 22.2 21.8 21.3 21.2 21.6 20.9 21.4 20.6 22.4 20.5 22.2 23.3
## [211] 22.6 21.9

Eksplorasi Data

Sebelum masuk dalam tahap pemodelan, dilakukan eksplorasi data dengan plot deret waktu untuk melihat pola data.

#--PLOT TIME SERIES--#
plot(datalat67.ts,
     col = "navyblue",
     lwd = 1,
     type = "o",
     xlab = "Periode",
     ylab = "Temperature",
     main = "Time Series Plot")

Berdasarkan plot data deret waktu penuh di atas, terlihat bahwa data memiliki pola trend sehingga pembagian data latih dan data uji ditetapkan dengan proporsi 80%:20%.

Plot Data Latih

ma2.train <- datalat67$Temperature[1:170]
train.ts<-ts(ma2.train)
plot.ts(train.ts, lty=1, xlab="Periode", ylab="Temperature", main="Plot Temperature Portugal Train")

Berdasarkan plot data deret waktu pada data latih, terlihat bahwa data cenderung memiliki trend yang naik. Hal ini mengindikasikan bahwa data tidak stasioner dalam rataan.

Plot Data Uji

ma2.test <- datalat67$Temperature[171:212]
test.ts<-ts(ma2.test)
plot.ts(test.ts, lty=1, xlab="Periode", ylab="Temperature", main="Plot Temperature Portugal Train")

Berdasarkan plot data deret waktu pada data uji di atas, terlihat pola data musiman yang tidak stasioner dalam rataan dan ragam. Data tidak stasioner dalam rataan karena tidak menyebar/bergerak di sekitar nilai tengahnya (0) dan dikatakan tidak stasioner dalam ragam karena memiliki lebar pita yang tidak sama. Selain dengan plot data deret waktu, akan dilakukan pengecekan stasioneritas data dengan plot ACF dan uji ADF.

Uji Kestasioneran Data

Plot ACF

acf(ma2.train, main="ACF", lag.max=20)

Berdasarkan plot ACF di atas, pada data tersebut menurun secara perlahan (tails off slowly) yang menandakan data tidak stasioner dalam rataan.

Uji ADF

adf.test(ma2.train)

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag      ADF p.value
## [1,]   0 -0.32341   0.551
## [2,]   1 -0.21153   0.583
## [3,]   2 -0.00109   0.643
## [4,]   3  0.21213   0.704
## [5,]   4  0.39258   0.756
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -2.11   0.285
## [2,]   1 -2.58   0.102
## [3,]   2 -2.07   0.301
## [4,]   3 -1.79   0.412
## [5,]   4 -1.55   0.504
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -4.98  0.0100
## [2,]   1 -5.32  0.0100
## [3,]   2 -4.65  0.0100
## [4,]   3 -4.19  0.0100
## [5,]   4 -3.83  0.0194
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

#stasioner

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.0194 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.

Plot Box-Cox

index <- seq(1:170)
bc = boxcox(train.ts~index, lambda = seq(-1,5,by=0.01))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] 1.44

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

##  [1] 1.04 1.05 1.06 1.07 1.08 1.09 1.10 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18
## [16] 1.19 1.20 1.21 1.22 1.23 1.24 1.25 1.26 1.27 1.28 1.29 1.30 1.31 1.32 1.33
## [31] 1.34 1.35 1.36 1.37 1.38 1.39 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48
## [46] 1.49 1.50 1.51 1.52 1.53 1.54 1.55 1.56 1.57 1.58 1.59 1.60 1.61 1.62 1.63
## [61] 1.64 1.65 1.66 1.67 1.68 1.69 1.70 1.71 1.72 1.73 1.74 1.75 1.76 1.77 1.78
## [76] 1.79 1.80 1.81 1.82 1.83

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 1.44 dan pada selang kepercayaan 95% nilai memiliki batas bawah 1.04 dan batas atas 1.83. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data yang digunakan tidak stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

train.diff1<-diff(train.ts,differences = 1) 
plot.ts(train.diff1, lty=1, xlab="Periode", ylab="Data Difference 1 Temperature", main="Plot Difference Temperature Portugal")

Berdasarkan plot data deret waktu, terlihat bahwa data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data).

Plot ACF

acf(train.diff1)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cuts off pada lag ke 2. Hal ini menandakan data sudah stasioner dalam rataan dan ketidakstasioneran data telah berhasil tertangani.

Uji ADF

tseries::adf.test(train.diff1)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff1
## Dickey-Fuller = -7.3987, 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\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan

Identifikasi Model

Plot ACF

acf(train.diff1)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 2, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,2).

Plot PACF

pacf(train.diff1)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 2, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(2,1,0).

Jika baik plot ACF maupun plot PACF keduanya dianggap tails of, maka model yang terbentuk adalah ARIMA(2,1,2)

Plot EACF

eacf(train.diff1)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o x x 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 o o o o o o o o o  o  o  o 
## 5 x x x x o o o o o o o  o  o  o 
## 6 x o x x o o o o o o o  o  o  o 
## 7 x o o x o x o o o o o  o  o  o

Identifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,3), ARIMA(1,1,2), dan ARIMA(2,1,1)

Pendugaan Parameter Model Tentatif

ARIMA(0,1,3)

model1.da=Arima(train.diff1, order=c(0,1,3),method="ML")
summary(model1.da) #AIC=612.02

## Series: train.diff1 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1      ma2     ma3
##       -0.8891  -0.3024  0.1916
## s.e.   0.1012   0.0937  0.1135
## 
## sigma^2 = 2.103:  log likelihood = -302.01
## AIC=612.02   AICc=612.27   BIC=624.52
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.1606816 1.432961 1.129997 Inf  Inf 0.7479884 0.06039065

lmtest::coeftest(model1.da) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.889142   0.101247 -8.7819 < 2.2e-16 ***
## ma2 -0.302430   0.093668 -3.2287  0.001243 ** 
## ma3  0.191572   0.113543  1.6872  0.091561 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(0,1,2)

model2.da=Arima(train.diff1, order=c(0,1,2),method="ML")
summary(model2.da) #AIC=613.21

## Series: train.diff1 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.7872  -0.2128
## s.e.   0.0867   0.0851
## 
## sigma^2 = 2.138:  log likelihood = -303.61
## AIC=613.21   AICc=613.36   BIC=622.58
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE         ACF1
## Training set 0.1528848 1.449196 1.141547 NaN  Inf 0.7556337 -0.001529718

lmtest::coeftest(model2.da) #seluruh parameter signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)    
## ma1 -0.787224   0.086709 -9.0789  < 2e-16 ***
## ma2 -0.212776   0.085145 -2.4990  0.01246 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(0,1,1)

model3.da=Arima(train.diff1, order=c(0,1,1),method="ML")
summary(model3.da) #AIC=616.67

## Series: train.diff1 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.0172
## 
## sigma^2 = 2.191:  log likelihood = -306.34
## AIC=616.67   AICc=616.74   BIC=622.92
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE      ACF1
## Training set 0.1609741 1.471454 1.150394 NaN  Inf 0.7614898 0.1712817

lmtest::coeftest(model3.da) #seluruh parameter signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -1.000000   0.017245 -57.988 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,1)

model4.da=Arima(train.diff1, order=c(1,1,1),method="ML")
summary(model4.da) #AIC=614.92

## Series: train.diff1 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.1539  -1.0000
## s.e.  0.0791   0.0168
## 
## sigma^2 = 2.16:  log likelihood = -304.46
## AIC=614.92   AICc=615.06   BIC=624.29
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.1545637 1.456447 1.145219 NaN  Inf 0.7580645 0.05146633

lmtest::coeftest(model4.da) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.153893   0.079086   1.9459  0.05167 .  
## ma1 -1.000000   0.016840 -59.3822  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,2)

model5.da=Arima(train.diff1, order=c(1,1,2),method="ML")
summary(model5.da) #AIC=614.24

## Series: train.diff1 
## ARIMA(1,1,2) 
## 
## Coefficients:
##           ar1      ma1      ma2
##       -0.2656  -0.5386  -0.4614
## s.e.   0.2801   0.2543   0.2539
## 
## sigma^2 = 2.137:  log likelihood = -303.12
## AIC=614.24   AICc=614.48   BIC=626.73
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 0.154219 1.444554 1.138806 Inf  Inf 0.7538196 0.007661311

lmtest::coeftest(model5.da) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1 -0.26560    0.28006 -0.9484  0.34294  
## ma1 -0.53861    0.25433 -2.1178  0.03419 *
## ma2 -0.46139    0.25388 -1.8173  0.06917 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,1)

model6.da=Arima(train.diff1, order=c(2,1,1),method="ML")
summary(model6.da) #AIC=610.49

## Series: train.diff1 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.1788  -0.1998  -1.0000
## s.e.  0.0781   0.0779   0.0174
## 
## sigma^2 = 2.086:  log likelihood = -301.25
## AIC=610.49   AICc=610.74   BIC=622.99
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE          ACF1
## Training set 0.1579057 1.426942 1.124374 Inf  Inf 0.7442668 -0.0006643083

lmtest::coeftest(model6.da) #seluruh parameter signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.178835   0.078092   2.2901  0.02202 *  
## ar2 -0.199801   0.077945  -2.5633  0.01037 *  
## ma1 -0.999999   0.017390 -57.5034  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan pendugaan parameter di atas, memiliki nilai AIC terkecil dan parameter model ARIMA(2,1,1) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(2,1,1).

Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

Eksplorasi Sisaan

#Eksplorasi 
sisaan.da <- model6.da$residuals 
par(mfrow=c(2,2)) 
qqnorm(sisaan.da) 
qqline(sisaan.da, col = "blue", lwd = 2) 
plot(c(1:length(sisaan.da)),sisaan.da) 
acf(sisaan.da) 
pacf(sisaan.da)

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan menyebar normal ditandai dengan titik titik yang cenderung mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(2,1,1) juga tidak signifikan pada 22 lag awal yang menandakan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

#1) Sisaan Menyebar Normal 
tseries::jarque.bera.test(sisaan.da)  #tak tolak H0 > sisaan menyebar normal

## 
##  Jarque Bera Test
## 
## data:  sisaan.da
## X-squared = 0.00099931, df = 2, p-value = 0.9995

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Jarque Bera . Hipotesis pada uji Jarque Bera adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 0.9995 yang lebih dari taraf nyata 5% sehingga tidak tolak \(H_0\) dan menandakan bahwa sisaan menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.

#2) Sisaan saling bebas/tidak ada autokorelasi 
Box.test(sisaan.da, type = "Ljung")  #tak tolak H0 > sisaan saling bebas

## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 7.5912e-05, df = 1, p-value = 0.993

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.993 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini sesuai dengan hasil eksplorasi.

#3) Sisaan homogen 
Box.test((sisaan.da)^2, type = "Ljung")  #tak tolak H0 > sisaan homogen

## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 2.8891, df = 1, p-value = 0.08918

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.08918 yang lebih 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.da, mu = 0, conf.level = 0.95)  #tak tolak h0 > nilai tengah sisaan sama dengan 0

## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = 1.4432, df = 168, p-value = 0.1508
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.05809939  0.37391088
## sample estimates:
## mean of x 
## 0.1579057

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.1508 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini sama dengan hasil eksplorasi.

Overfitting

Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(2,1,1) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,1,2) dan ARIMA(0,1,3).

Model ARIMA(1,1,2)

model.overfit1=Arima(train.diff1, order=c(3,1,1),method="ML")
summary(model.overfit1) #AIC=609.82

## Series: train.diff1 
## ARIMA(3,1,1) 
## 
## Coefficients:
##          ar1      ar2      ar3     ma1
##       0.1498  -0.1803  -0.1307  -1.000
## s.e.  0.0793   0.0781   0.0796   0.018
## 
## sigma^2 = 2.061:  log likelihood = -299.91
## AIC=609.82   AICc=610.19   BIC=625.44
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.1611987 1.414338 1.096872 Inf  Inf 0.7260619 0.01115869

lmtest::coeftest(model.overfit1) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.149813   0.079331   1.8885  0.05896 .  
## ar2 -0.180252   0.078131  -2.3071  0.02105 *  
## ar3 -0.130749   0.079551  -1.6436  0.10026    
## ma1 -1.000000   0.017955 -55.6958  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(0,1,3)

model.overfit2=Arima(train.diff1, order=c(2,1,3),method="ML")
summary(model.overfit2) #AIC=599.08

## Series: train.diff1 
## ARIMA(2,1,3) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     ma3
##       0.8675  -0.1635  -1.8204  0.6555  0.1728
## s.e.  0.2530   0.2120   0.3329  0.5589  0.2438
## 
## sigma^2 = 1.877:  log likelihood = -293.54
## AIC=599.08   AICc=599.6   BIC=617.82
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.1375622 1.345325 1.062468 NaN  Inf 0.7032887 0.02391455

lmtest::coeftest(model.overfit2) ##terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.86754    0.25303  3.4285 0.0006069 ***
## ar2 -0.16349    0.21201 -0.7712 0.4406063    
## ma1 -1.82038    0.33292 -5.4680 4.553e-08 ***
## ma2  0.65553    0.55886  1.1730 0.2408016    
## ma3  0.17280    0.24383  0.7087 0.4785114    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#model yang dipilih adalah model awal, yaitu ARIMA(2,1,1)

Berdasarkan kedua model hasil overfitting di atas, model ARIMA(3,1,1) dan ARIMA(2,1,3) memiliki AIC yang lebih kecil dibandingkan dengan model ARIMA(2,1,1) dan parameter kedua model ARIMA(1,0,2) dan ARIMA(0,0,3) tidak seluruhnya signifikan. Oleh karena itu, model ARIMA(2,1,1) akan tetap digunakan untuk melakukan peramalan.

Peramalan

Peramalan dilakukan menggunakan fungsi forecast() . Contoh peramalan berikut ini dilakukan untuk 30 periode ke depan.

#---FORECAST---#
ramalan.da <- forecast::forecast(model6.da, h = 42) 
ramalan.da

##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 171   0.0473447536 -1.808869 1.903558 -2.791489 2.886179
## 172   0.1988798540 -1.688709 2.086469 -2.687938 3.085698
## 173   0.0566795845 -1.854338 1.967697 -2.865970 2.979329
## 174   0.0009723446 -1.913233 1.915177 -2.926551 2.928496
## 175   0.0194216422 -1.895428 1.934272 -2.909089 2.947932
## 176   0.0338513729 -1.881434 1.949137 -2.895325 2.963028
## 177   0.0327457345 -1.882528 1.948019 -2.896412 2.961904
## 178   0.0296649369 -1.885582 1.944912 -2.899453 2.958783
## 179   0.0293348891 -1.885909 1.944578 -2.899777 2.958447
## 180   0.0298914104 -1.885359 1.945142 -2.899232 2.959015
## 181   0.0300568798 -1.885196 1.945310 -2.899069 2.959183
## 182   0.0299752782 -1.885276 1.945227 -2.899150 2.959100
## 183   0.0299276241 -1.885324 1.945179 -2.899196 2.959052
## 184   0.0299354059 -1.885316 1.945187 -2.899189 2.959059
## 185   0.0299463189 -1.885305 1.945198 -2.899178 2.959071
## 186   0.0299467157 -1.885305 1.945198 -2.899178 2.959071
## 187   0.0299446063 -1.885307 1.945196 -2.899180 2.959069
## 188   0.0299441497 -1.885307 1.945196 -2.899180 2.959068
## 189   0.0299444896 -1.885307 1.945196 -2.899180 2.959069
## 190   0.0299446415 -1.885307 1.945196 -2.899180 2.959069
## 191   0.0299446008 -1.885307 1.945196 -2.899180 2.959069
## 192   0.0299445632 -1.885307 1.945196 -2.899180 2.959069
## 193   0.0299445646 -1.885307 1.945196 -2.899180 2.959069
## 194   0.0299445724 -1.885307 1.945196 -2.899180 2.959069
## 195   0.0299445735 -1.885307 1.945196 -2.899180 2.959069
## 196   0.0299445721 -1.885307 1.945196 -2.899180 2.959069
## 197   0.0299445717 -1.885307 1.945196 -2.899180 2.959069
## 198   0.0299445718 -1.885307 1.945196 -2.899180 2.959069
## 199   0.0299445720 -1.885307 1.945196 -2.899180 2.959069
## 200   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 201   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 202   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 203   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 204   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 205   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 206   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 207   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 208   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 209   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 210   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 211   0.0299445719 -1.885307 1.945196 -2.899180 2.959069
## 212   0.0299445719 -1.885307 1.945196 -2.899180 2.959069

data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(0,1,2) cenderung stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.

pt_1 <- train.ts[170] #nilai akhir data latih
hasil.forc.Diff <- data.ramalan.da
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1

#has.1 sama hasilnta dengan: cumsum(c(pt_1,hasil.forc.Diff))
ts.plot(train.ts,hasil)

perbandingan.da<-matrix(data=c(head(test.ts, n=42), hasil[-1]),
                     nrow = 42, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da

##       Aktual Hasil Forecast
##  [1,]   20.4       20.54734
##  [2,]   19.9       20.74622
##  [3,]   21.7       20.80290
##  [4,]   26.0       20.80388
##  [5,]   27.4       20.82330
##  [6,]   25.1       20.85715
##  [7,]   24.7       20.88990
##  [8,]   23.7       20.91956
##  [9,]   23.6       20.94890
## [10,]   21.5       20.97879
## [11,]   21.9       21.00884
## [12,]   23.7       21.03882
## [13,]   24.6       21.06875
## [14,]   22.2       21.09868
## [15,]   21.5       21.12863
## [16,]   21.4       21.15857
## [17,]   21.5       21.18852
## [18,]   21.1       21.21846
## [19,]   21.9       21.24841
## [20,]   21.5       21.27835
## [21,]   23.3       21.30830
## [22,]   22.5       21.33824
## [23,]   22.8       21.36819
## [24,]   23.9       21.39813
## [25,]   22.7       21.42808
## [26,]   21.7       21.45802
## [27,]   20.3       21.48796
## [28,]   22.0       21.51791
## [29,]   22.2       21.54785
## [30,]   21.8       21.57780
## [31,]   21.3       21.60774
## [32,]   21.2       21.63769
## [33,]   21.6       21.66763
## [34,]   20.9       21.69758
## [35,]   21.4       21.72752
## [36,]   20.6       21.75747
## [37,]   22.4       21.78741
## [38,]   20.5       21.81735
## [39,]   22.2       21.84730
## [40,]   23.3       21.87724
## [41,]   22.6       21.90719
## [42,]   21.9       21.93713

accuracy(ts(hasil[-1]), head(test.ts, n=42))

##                ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
## Test set 1.022531 1.989376 1.343965 4.114618 5.675175 0.6913528  1.409724

Model tersebut merupakan model yang sangat baik dengan nilai MAPE yang kurang dari 10%.

Tugas 67 MPDW

Arfiah Kania Sektiaruni

2023-10-02

Packages

Impor Data

Mengubah data agar terbaca sebagai data deret waktu dengan fungsi `ts()` .

Eksplorasi Data

Plot Data Latih

Plot Data Uji

Uji Kestasioneran Data

Plot ACF

Uji ADF

Plot Box-Cox

Penanganan Ketidakstasioneran Data

Plot ACF

Uji ADF

Identifikasi Model

Plot ACF

Plot PACF

Plot EACF

Pendugaan Parameter Model Tentatif

ARIMA(0,1,3)

ARIMA(0,1,2)

ARIMA(0,1,1)

ARIMA(1,1,1)

ARIMA(1,1,2)

ARIMA(2,1,1)

Analisis Sisaan

Eksplorasi Sisaan

Uji Formal

Overfitting

Model ARIMA(1,1,2)

Model ARIMA(0,1,3)

Peramalan

Tugas 67 MPDW

Arfiah Kania Sektiaruni

2023-10-02

Packages

Impor Data

Mengubah data agar terbaca sebagai data deret waktu dengan fungsi ts() .

Eksplorasi Data

Plot Data Latih

Plot Data Uji

Uji Kestasioneran Data

Plot ACF

Uji ADF

Plot Box-Cox

Penanganan Ketidakstasioneran Data

Plot ACF

Uji ADF

Identifikasi Model

Plot ACF

Plot PACF

Plot EACF

Pendugaan Parameter Model Tentatif

ARIMA(0,1,3)

ARIMA(0,1,2)

ARIMA(0,1,1)

ARIMA(1,1,1)

ARIMA(1,1,2)

ARIMA(2,1,1)

Analisis Sisaan

Eksplorasi Sisaan

Uji Formal

Overfitting

Model ARIMA(1,1,2)

Model ARIMA(0,1,3)

Peramalan

Mengubah data agar terbaca sebagai data deret waktu dengan fungsi `ts()` .