MPDW Praktikum 67
G1401211094 Hanifa Rahmacindia Nasution
2023-10-16
Packages
library(rio)
library(ggplot2)
library(tsibble)##
## Attaching package: 'tsibble'## The following objects are masked from 'package:base':
##
## intersect, setdiff, unionlibrary(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(MASS)
library(forecast)
library(TSA)## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(TTR)
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':
##
## identifylibrary(graphics)Input Data
gas <- rio::import("https://raw.githubusercontent.com/HanifaNasution/Praktikum-MPDW/main/Data/Data%20Natural%20Gas%20Imputasi.csv")
View(gas)
str(gas)## 'data.frame': 366 obs. of 8 variables:
## $ Period : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Date : chr "9/1/2021" "9/2/2021" "9/3/2021" "9/4/2021" ...
## $ Open : num 4.4 4.64 4.63 NA 4.73 ...
## $ High : num 4.71 4.73 4.72 NA 4.74 ...
## $ Low : num 4.38 4.58 4.59 NA 4.71 ...
## $ Close : num 4.62 4.64 4.71 NA 4.72 ...
## $ Volume : int 241159 176845 110953 NA NA NA 163663 244742 234848 151036 ...
## $ Currency: chr "USD" "USD" "USD" "USD" ...dim(gas)## [1] 366 8Ubah data tanggal ke mode “date”
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:tsibble':
##
## interval## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(stringr)
gas_clean <- gas %>%
mutate(Date = mdy(Date))
str(gas_clean)## 'data.frame': 366 obs. of 8 variables:
## $ Period : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Date : Date, format: "2021-09-01" "2021-09-02" ...
## $ Open : num 4.4 4.64 4.63 NA 4.73 ...
## $ High : num 4.71 4.73 4.72 NA 4.74 ...
## $ Low : num 4.38 4.58 4.59 NA 4.71 ...
## $ Close : num 4.62 4.64 4.71 NA 4.72 ...
## $ Volume : int 241159 176845 110953 NA NA NA 163663 244742 234848 151036 ...
## $ Currency: chr "USD" "USD" "USD" "USD" ...Tampak bahwa kolom tanggal telah menggunakan mode “Date” bukan “chr” lagi
Imputasi Data Hilang
library(imputeTS)##
## Attaching package: 'imputeTS'## The following object is masked from 'package:tseries':
##
## na.removegas.open<- na_interpolation(gas_clean$Open)
ggplot_na_imputations(gas_clean$Open,gas.open)
gas.high<- na_interpolation(gas_clean$High)
ggplot_na_imputations(gas_clean$High,gas.high)
gas.low<- na_interpolation(gas_clean$Low)
ggplot_na_imputations(gas_clean$Low,gas.low)
gas.close <- na_interpolation(gas_clean$Close)
ggplot_na_imputations(gas_clean$Close,gas.close)
gas.volume <- na_interpolation(gas_clean$Volume)
ggplot_na_imputations(gas_clean$Volume,gas.volume)
Data Penuh Pasca Imputasi
gas.imp <- data.frame(gas_clean$Period,gas_clean$Date,cbind(gas.open,gas.high,gas.low,gas.close,gas.volume),gas_clean$Currency)
View(gas.imp)
str(gas.imp)## 'data.frame': 366 obs. of 8 variables:
## $ gas_clean.Period : int 1 2 3 4 5 6 7 8 9 10 ...
## $ gas_clean.Date : Date, format: "2021-09-01" "2021-09-02" ...
## $ gas.open : num 4.4 4.64 4.63 4.68 4.73 ...
## $ gas.high : num 4.71 4.73 4.72 4.73 4.74 ...
## $ gas.low : num 4.38 4.58 4.59 4.65 4.71 ...
## $ gas.close : num 4.62 4.64 4.71 4.72 4.72 ...
## $ gas.volume : num 241159 176845 110953 124131 137308 ...
## $ gas_clean.Currency: chr "USD" "USD" "USD" "USD" ...summary(gas.imp)## gas_clean.Period gas_clean.Date gas.open gas.high
## Min. : 1.00 Min. :2021-09-01 Min. :3.589 Min. :3.794
## 1st Qu.: 92.25 1st Qu.:2021-12-01 1st Qu.:4.661 1st Qu.:4.794
## Median :183.50 Median :2022-03-02 Median :5.484 Median :5.710
## Mean :183.50 Mean :2022-03-02 Mean :5.969 Mean :6.193
## 3rd Qu.:274.75 3rd Qu.:2022-06-01 3rd Qu.:7.306 3rd Qu.:7.678
## Max. :366.00 Max. :2022-09-01 Max. :9.780 Max. :9.987
## gas.low gas.close gas.volume gas_clean.Currency
## Min. :3.536 Min. :3.561 Min. : 1200 Length:366
## 1st Qu.:4.522 1st Qu.:4.684 1st Qu.: 78878 Class :character
## Median :5.349 Median :5.503 Median :115310 Mode :character
## Mean :5.776 Mean :5.981 Mean :119401
## 3rd Qu.:7.045 3rd Qu.:7.384 3rd Qu.:157808
## Max. :9.200 Max. :9.647 Max. :305898Pada peramalan ini hanya digunakan data saham ketika mencapai nilai tertinggi, yaitu data saham kolom high
Eksplorasi Data
Plot Deret Waktu Data Penuh
gas.high.ts <- ts(gas.high)
plot(gas.high.ts, lty=1, xlab="Waktu", ylab="Kurs Gas Alam", main="Plot Data Kurs Gas Alam")
Berdasarkan plot yang dihasilkan data penuh dicoba untuk memotong data dengan ketentuan 80% data latih dan 20% data uji.
Pembagian Data Latih dan Data Uji (Percobaan 1)
Data Latih 80% data awal Data Uji 20% data awal
gas.train <- gas.high[1:293]
gas.test <- gas.high[294:366]
gas.train.ts <- ts(gas.train)
gas.test.ts <- ts(gas.test)Plot Deret waktu Data Latih (Percobaan 1)
plot(gas.train.ts, lty=1, xlab="Waktu", ylab="Kurs Gas Alam", main="Plot Data Latih 1 Kurs Gas Alam")
Berdasarkan plot deret waktu data latih saham gas alam tampak bahwa data cenderung memiliki trend naik dan tidak bergerak dalam satu angka tertentu sehingga data latih tidak stasioner baik dalam rataan maupun ragam Namun, pemotongan ini tidak tepat karena akan memotong data ketika akan naik setelah data turun sehingga disinyalir akan memberikan peramalan yang kurang tepat sehingga harus dicoba pembagian data lainnya
Pembagian Data Latih dan Data Uji (Percobaan 2)
Data Latih 86% data awal Data Uji 14% data awal
gas.train <- gas.high[1:315]
gas.test <- gas.high[316:366]
gas.high.ts <- ts(gas.high)
gas.train.ts <- ts(gas.train)
gas.test.ts <- ts(gas.test)Plot Deret waktu Data Latih (Percobaan 2)
plot(gas.train.ts, lty=1, xlab="Waktu", ylab="Kurs Gas Alam", main="Plot Data Latih 2 Kurs Gas Alam")
Berdasarkan plot deret waktu data latih saham gas alam tampak bahwa data cenderung memiliki trend naik dan tidak bergerak dalam satu angka tertentu sehingga data latih tidak stasioner baik dalam rataan maupun ragam Pemotongan ini dianggap lebih baik karena memotong data ketika data aktual akan naik terus sehingga disinyalir hasil peramalan akan lebih baik nantinya
Plot Deret Waktu Data Uji
plot(gas.test.ts, lty=1, xlab="Waktu", ylab="Kurs Gas Alam", main="Plot Data Uji Kurs Gas Alam")
Uji Stasioneritas Data
Plot ACF
acf(gas.train.ts)
Berdasarkan plot ACF data latih kurs gas alam, tampak bahwa plot ACF tails off slowly. Hal ini menandakan bahwa data latih kurs gas alam tidak stasioner dalam rataan.
Uji ADF
tseries::adf.test(gas.train.ts)##
## Augmented Dickey-Fuller Test
##
## data: gas.train.ts
## Dickey-Fuller = -1.6785, Lag order = 6, p-value = 0.712
## 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.712 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani
Plot Box-Cox
index <- seq(1:315)
bc = boxcox(gas.train.ts~index, lambda = seq(-2,2,by=1))
#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda## [1] -0.3838384#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]## [1] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263 -0.58585859
## [7] -0.54545455 -0.50505051 -0.46464646 -0.42424242 -0.38383838 -0.34343434
## [13] -0.30303030 -0.26262626 -0.22222222 -0.18181818 -0.14141414 -0.10101010
## [19] -0.06060606 -0.02020202 0.02020202Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar -0,3838384 dan pada selang kepercayaan 95% nilai memiliki batas bawah -0,7878789 dan batas atas 0,020202. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data latih kurs gas alam tidak stasioner dalam ragam.
Penanganan Ketidakstasioneran Data
gas.train.diff<-diff(gas.train.ts,differences = 1)
plot.ts(gas.train.diff, lty=1, xlab="Waktu", ylab="Data Kurs Difference 1", main="Plot Kurs Gas Alam Differencing")
Berdasarkan plot data deret waktu, terlihat bahwa data sudah mulai tampak stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data)
Plot ACF
acf(gas.train.diff)
Berdasarkan plot tersebut, terlihat bahwa plot ACF cuts off pada lag ke 1. Hal ini menandakan data sudah stasioner dalam rataan dan ketidakstasioneran data telah berhasil tertangani.
Uji ADF
tseries::adf.test(gas.train.diff)## Warning in tseries::adf.test(gas.train.diff): p-value smaller than printed
## p-value##
## Augmented Dickey-Fuller Test
##
## data: gas.train.diff
## Dickey-Fuller = -7.1758, Lag order = 6, 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 pasca differencing, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan
Identifikasi Model
Plot ACF
acf(gas.train.diff)
Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 1 dan jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,1).
Plot PACF
pacf(gas.train.diff)
Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 1, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(1,1,0).
Jika baik plot ACF maupun plot PACF keduanya dianggap tails of, maka model yang terbentuk adalah ARIMA(1,1,1)
Plot EACF
eacf(gas.train.diff)## 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 o 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 o x x o o o o o o o o o o
## 5 x x o x o o o o o o o o o o
## 6 x x o x x o o o o o o o o o
## 7 x o x x o x o o o o o o o oIdentifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(2,1,2), dan ARIMA(3,1,3).
Berdasarkan model yang didapat dari plot ACF, PACF, dan EACF didapat 5 model tentatif, yaitu ARIMA (1,1,0), ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(2,1,2), dan ARIMA(3,1,3).
Pendugaan Parameter Model Tentatif
ARIMA(0,1,1)
model1=Arima(gas.train.diff, order=c(0,1,1),method="ML")
summary(model1) #AIC=-78.97## Series: gas.train.diff
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -1.0000
## s.e. 0.0143
##
## sigma^2 = 0.04424: log likelihood = 41.48
## AIC=-78.97 AICc=-78.93 BIC=-71.48
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.003126911 0.2096627 0.1433882 -Inf Inf 0.8499932 0.1399336lmtest::coeftest(model1) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.00000 0.01432 -69.834 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,0)
model2=Arima(gas.train.diff, order=c(1,1,0),method="ML")
summary(model2) #AIC=10.7## Series: gas.train.diff
## ARIMA(1,1,0)
##
## Coefficients:
## ar1
## -0.4589
## s.e. 0.0501
##
## sigma^2 = 0.05996: log likelihood = -3.35
## AIC=10.7 AICc=10.73 BIC=18.19
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.0002335151 0.2440907 0.1616892 -Inf Inf 0.9584801 -0.1439141lmtest::coeftest(model2) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.45889 0.05008 -9.163 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,1)
model3=Arima(gas.train.diff, order=c(1,1,1),method="ML")
summary(model3) #AIC=-83.45## Series: gas.train.diff
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.1429 -1.0000
## s.e. 0.0560 0.0114
##
## sigma^2 = 0.04351: log likelihood = 44.72
## AIC=-83.45 AICc=-83.37 BIC=-72.21
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.002663764 0.2075992 0.1408843 -Inf Inf 0.8351508 -0.0104734lmtest::coeftest(model3) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.142868 0.055988 2.5518 0.01072 *
## ma1 -0.999999 0.011378 -87.8871 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(2,1,2)
model4=Arima(gas.train.diff, order=c(2,1,2),method="ML")
summary(model4) #AIC=-80.65## Series: gas.train.diff
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.3584 0.0282 -1.2242 0.2242
## s.e. 0.4425 0.0892 0.4398 0.4397
##
## sigma^2 = 0.04365: log likelihood = 45.32
## AIC=-80.65 AICc=-80.45 BIC=-61.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.002381824 0.2072526 0.1393907 -Inf Inf 0.8262969 -0.003272745lmtest::coeftest(model4) #hanya parameter ma1 signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.358376 0.442540 0.8098 0.418045
## ar2 0.028190 0.089168 0.3161 0.751894
## ma1 -1.224152 0.439834 -2.7832 0.005382 **
## ma2 0.224152 0.439711 0.5098 0.610211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(3,1,3)
model5=Arima(gas.train.diff, order=c(3,1,3),method="ML")
summary(model5) #AIC=-77.67## Series: gas.train.diff
## ARIMA(3,1,3)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## -0.2829 -0.1679 0.113 -0.5821 -0.1432 -0.2747
## s.e. 1.1341 0.6188 0.127 1.1434 0.6977 0.7212
##
## sigma^2 = 0.04379: log likelihood = 45.84
## AIC=-77.67 AICc=-77.31 BIC=-51.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.002368479 0.2069087 0.1393582 -Inf Inf 0.8261039 -0.003150472lmtest::coeftest(model5) #tidak ada parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.28294 1.13411 -0.2495 0.8030
## ar2 -0.16790 0.61876 -0.2713 0.7861
## ar3 0.11304 0.12699 0.8901 0.3734
## ma1 -0.58212 1.14336 -0.5091 0.6107
## ma2 -0.14322 0.69765 -0.2053 0.8374
## ma3 -0.27466 0.72120 -0.3808 0.7033Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(1,1,1) dan parameter model ARIMA(1,1,1) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(1,1,1).
Analisis Sisaan
Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.
Eksplorasi Sisaan
#Eksplorasi
sisaan.model3 <- model3$residuals
par(mfrow=c(2,2))
qqnorm(sisaan.model3)
qqline(sisaan.model3, col = "blue", lwd = 2)
plot(c(1:length(sisaan.model3)),sisaan.model3)
acf(sisaan.model3)
pacf(sisaan.model3)
par(mfrow = c(1,1))Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan tidak menyebar normal karena tidak mengikuti garis \(45^{\circ}\) secara penuh. Kemudian dapat dilihat juga lebar pita sisaan yang tidak sama menandakan bahwa sisaan memiliki ragam yang tidak homogen. Selain itu, plot ACF dan PACF sisaan ARIMA(1,1,1) signifikan pada lag ke-13 sehingga sisaan tidak saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.
Uji Formal
#1) Sisaan Menyebar Normal
ks.test(sisaan.model3,"pnorm") ##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan.model3
## D = 0.3376, p-value < 2.2e-16
## alternative hypothesis: two-sided#tolak H0 > sisaan tidak menyebar normalSelain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.
\(H_0\) : Sisaan menyebar normal
\(H_1\) : Sisaan tidak menyebar normal
Berdasarkan uji KS tersebut, didapat p-value sebesar 2.2e-16 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.
#2) Sisaan saling bebas/tidak ada autokorelasi
Box.test(sisaan.model3, type = "Ljung") ##
## Box-Ljung test
##
## data: sisaan.model3
## X-squared = 0.034773, df = 1, p-value = 0.8521#tak tolak H0 > sisaan saling bebasSelanjutnya 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.8521 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini berbeda dengan eksplorasi.
#3) Sisaan homogen
Box.test((sisaan.model3)^2, type = "Ljung") ##
## Box-Ljung test
##
## data: (sisaan.model3)^2
## X-squared = 2.1754, df = 1, p-value = 0.1402#tak tolak H0 > sisaan homogenHipotesis 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.1402 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen. Hal ini berbeda dengan eksplorasi
#4) Nilai tengah sisaan sama dengan nol
t.test(sisaan.model3, mu = 0, conf.level = 0.95) ##
## One Sample t-test
##
## data: sisaan.model3
## t = -0.22703, df = 313, p-value = 0.8206
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.02574975 0.02042222
## sample estimates:
## mean of x
## -0.002663764#tak tolak h0 > nilai tengah sisaan sama dengan 0Terakhir, 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.8206 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini sesuai dengan eksplorasi.
Overfitting
Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(1,1,1) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(2,1,1) dan ARIMA(1,1,2).
#---OVERFITTING---#
model1a.model3=Arima(gas.train.diff, order=c(2,1,1),method="ML")
summary(model1a.model3) #-82.43## Series: gas.train.diff
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.1355 0.0561 -1.0000
## s.e. 0.0564 0.0564 0.0108
##
## sigma^2 = 0.04353: log likelihood = 45.22
## AIC=-82.43 AICc=-82.3 BIC=-67.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.002463106 0.2073082 0.1396672 -Inf Inf 0.827936 -0.00537233lmtest::coeftest(model1a.model3) #ar2 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.135528 0.056440 2.4013 0.01634 *
## ar2 0.056070 0.056381 0.9945 0.31999
## ma1 -1.000000 0.010754 -92.9881 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model1b.model3=Arima(gas.train.diff, order=c(1,1,2),method="ML")
summary(model1b.model3) #-82.55## Series: gas.train.diff
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.4687 -1.3297 0.3297
## s.e. 0.2330 0.2464 0.2462
##
## sigma^2 = 0.04352: log likelihood = 45.28
## AIC=-82.55 AICc=-82.43 BIC=-67.57
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0023891 0.2072814 0.1395654 -Inf Inf 0.827332 -0.007345763lmtest::coeftest(model1b.model3) #ma2 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.46866 0.23303 2.0112 0.04431 *
## ma1 -1.32973 0.24642 -5.3961 6.811e-08 ***
## ma2 0.32973 0.24622 1.3392 0.18051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model yang dipilih adalah model awal, yaitu ARIMA(1,1,1)Berdasarkan kedua model hasil overfitting di atas, model ARIMA(2,1,1) dan ARIMA(1,1,2) memiliki AIC yang lebih besar dibandingkan dengan model ARIMA(1,1,1) dan parameter kedua model ARIMA(2,1,1) dan ARIMA(1,1,2) tidak seluruhnya signifikan. Oleh karena itu, model ARIMA(1,1,1) akan tetap digunakan untuk melakukan peramalan.
Peramalan
Peramalan dilakukan menggunakan fungsi forecast() .
Contoh peramalan berikut ini dilakukan untuk 51 periode ke depan
(sebanyak data test).
#---FORECAST---#
ramalan.da <- forecast::forecast(sisaan.model3, h = 51)
ramalan.da## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 316 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 317 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 318 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 319 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 320 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 321 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 322 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 323 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 324 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 325 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 326 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 327 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 328 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 329 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 330 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 331 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 332 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 333 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 334 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 335 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 336 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 337 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 338 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 339 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 340 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 341 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 342 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 343 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 344 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 345 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 346 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 347 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 348 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 349 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054607
## 350 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054608
## 351 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054608
## 352 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054608
## 353 -0.002715096 -0.2696069 0.2641767 -0.4108909 0.4054608
## 354 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 355 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 356 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 357 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 358 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 359 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 360 -0.002715096 -0.2696069 0.2641768 -0.4108910 0.4054608
## 361 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608
## 362 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608
## 363 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608
## 364 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608
## 365 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608
## 366 -0.002715096 -0.2696070 0.2641768 -0.4108910 0.4054608data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(1,1,1) cenderung stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.
pt_1 <- gas.train.ts[315] #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(gas.train.ts,hasil)
perbandingan.da<-matrix(data=c(head(gas.test.ts, n=51), hasil[-1]),
nrow = 51, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da## Aktual Hasil Forecast
## [1,] 6.683000 6.644285
## [2,] 6.806000 6.641570
## [3,] 7.047000 6.638855
## [4,] 7.184667 6.636140
## [5,] 7.322333 6.633425
## [6,] 7.460000 6.630709
## [7,] 7.453000 6.627994
## [8,] 7.925000 6.625279
## [9,] 8.045000 6.622564
## [10,] 8.293000 6.619849
## [11,] 8.430000 6.617134
## [12,] 8.567000 6.614419
## [13,] 8.704000 6.611704
## [14,] 9.419000 6.608989
## [15,] 8.990000 6.606274
## [16,] 8.845000 6.603558
## [17,] 8.388000 6.600843
## [18,] 8.358000 6.598128
## [19,] 8.328000 6.595413
## [20,] 8.298000 6.592698
## [21,] 8.233000 6.589983
## [22,] 8.480000 6.587268
## [23,] 8.450000 6.584553
## [24,] 8.248000 6.581838
## [25,] 8.137667 6.579123
## [26,] 8.027333 6.576408
## [27,] 7.917000 6.573692
## [28,] 7.889000 6.570977
## [29,] 8.267000 6.568262
## [30,] 8.994000 6.565547
## [31,] 8.919000 6.562832
## [32,] 8.924667 6.560117
## [33,] 8.930333 6.557402
## [34,] 8.936000 6.554687
## [35,] 9.411000 6.551972
## [36,] 9.677000 6.549257
## [37,] 9.663000 6.546541
## [38,] 9.395000 6.543826
## [39,] 9.580000 6.541111
## [40,] 9.765000 6.538396
## [41,] 9.950000 6.535681
## [42,] 9.987000 6.532966
## [43,] 9.408000 6.530251
## [44,] 9.419000 6.527536
## [45,] 9.668000 6.524821
## [46,] 9.672667 6.522106
## [47,] 9.677333 6.519391
## [48,] 9.682000 6.516675
## [49,] 9.300000 6.513960
## [50,] 9.284000 6.511245
## [51,] 9.394000 6.508530accuracy(ts(hasil[-1]), head(gas.test.ts, n=51))## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 2.086965 2.267502 2.086965 23.26176 23.26176 0.8988544 8.489344Berdasarkan nilai MAPE yang didapat, yaitu sebesar 23% pada perbandingan nilai forecast dan aktual, dapat disimpulkan bahwa model yang dibangun kurang baik dan kurang sesuai terhadap data gas alam yang digunakan, sehingga perlu mencoba model lain (model selain ARIMA).