Penentuan Model Deret Waktu dan Peramalan Data
Radja Ahmad Nur Fikri (G14190028)
[]
Load library yang dibutuhkan
library("TTR")
library("TSA")
library("graphics")
library("forecast")
library(tidyverse)
library(tseries)
library(lubridate)
library(gridExtra)
library(ggfortify)#Load seluruh data yang akan digunakan
Data yang akan digunakan ada tiga, yaitu suhu minimum harian, produksi listrik bulanan, dan produksi bir bulanan. Ketiga data didapatkan dari bahan kuliah Metode Peramalan Deret Waktu pertemuan ke-7.
temp <- read.csv("daily-minimum-temperatures-in-me.csv")
elec <- read.csv("Electric_Production.csv")
beer <- read.csv("monthly-beer-production-in-austr.csv")Untuk memastikan apakah data sudah ter-import dengan benar, kita dapat melihat lima data pertama dari tiap file yang telah di-import.
str(temp, strict.width="cut")## 'data.frame': 365 obs. of 2 variables:
## $ Date : chr "1/1/1990" "1/2/1990" "1/3/1990" "1/4/199"..
## $ Daily.minimum.temperatures: num 14.8 13.3 15.6 14.5 14.3 15.3 16.4 14.8 17..str(elec, strict.width="cut")## 'data.frame': 241 obs. of 2 variables:
## $ DATE : chr "1/1/1998" "2/1/1998" "3/1/1998" "4/1/1998" ...
## $ IPG2211A2N: num 94.8 87.8 86.6 76.8 78 ...str(beer, strict.width="cut")## 'data.frame': 116 obs. of 2 variables:
## $ Month : chr "1986-01" "1986-02" "1986-03" "1986-04" ...
## $ Monthly.beer.production: num 161 156 142 165 136 ...Mengubah tanggal menjadi format yang sesuai
Dari struktur data tersebut, terlihat bahwa variabel date atau tanggal masih dalam bentuk karakter. Untuk mengubah format peubah date, maka dapat dilakukan dengan syntax sebagai berikut.
temp$Date <- mdy(temp$Date)
elec$DATE <- mdy(elec$DATE)
beer$Month <- ym(beer$Month)Lakukan cek ulang apakah format tanggal sudah sesuai dan sudah disimpan dalam format tanggal yang benar.
str(temp, strict.width="cut")## 'data.frame': 365 obs. of 2 variables:
## $ Date : Date, format: "1990-01-01" "1990-01-02" ...
## $ Daily.minimum.temperatures: num 14.8 13.3 15.6 14.5 14.3 15.3 16.4 14.8 17..str(elec, strict.width="cut")## 'data.frame': 241 obs. of 2 variables:
## $ DATE : Date, format: "1998-01-01" "1998-02-01" ...
## $ IPG2211A2N: num 94.8 87.8 86.6 76.8 78 ...str(beer, strict.width="cut")## 'data.frame': 116 obs. of 2 variables:
## $ Month : Date, format: "1986-01-01" "1986-02-01" ...
## $ Monthly.beer.production: num 161 156 142 165 136 ...Melihat lima data pertama
Setelah itu, mari lihat lima data pertama dari masing-masing data
yang sudah kita masukkan ke dalam R.
head(temp) %>%
knitr::kable(digits = 3, align = "c",
col.names = c("Tanggal", "Suhu Minimum Harian"))| Tanggal | Suhu Minimum Harian |
|---|---|
| 1990-01-01 | 14.8 |
| 1990-01-02 | 13.3 |
| 1990-01-03 | 15.6 |
| 1990-01-04 | 14.5 |
| 1990-01-05 | 14.3 |
| 1990-01-06 | 15.3 |
head(elec) %>% knitr::kable(digits = 3, align = "c",
col.names = c("Tanggal", "Produksi Listrik Bulanan"))| Tanggal | Produksi Listrik Bulanan |
|---|---|
| 1998-01-01 | 94.792 |
| 1998-02-01 | 87.820 |
| 1998-03-01 | 86.555 |
| 1998-04-01 | 76.752 |
| 1998-05-01 | 78.030 |
| 1998-06-01 | 86.458 |
head(beer) %>% knitr::kable(digits = 3, align = "c",
col.names = c("Tanggal", "Produksi Beer Bulanan"))| Tanggal | Produksi Beer Bulanan |
|---|---|
| 1986-01-01 | 161.2 |
| 1986-02-01 | 155.5 |
| 1986-03-01 | 141.9 |
| 1986-04-01 | 164.6 |
| 1986-05-01 | 136.2 |
| 1986-06-01 | 126.8 |
Pisahkan data training dan data testing
Data training adalah 80% dari seluruh data yang tersedia, dan sisanya akan dijadikan sebagai data testing. Pendugaan model data akan menggunakan data testing dari masing-masing dataset yang tersedia. Kemduian hasil dari peramalan data akan dibandingkan dengan data testing yang telah dibuat.
#Data suhu
n_temp <- nrow(temp)
train_temp <- temp[1:(0.8*n_temp),]
test_temp <- temp[(0.8*n_temp+1):n_temp,]
#Data listrik
n_elec <- nrow(elec)
train_elec <- elec[1:(0.8*n_elec),]
test_elec <- elec[(0.8*n_elec+1):nrow(elec),]
#Data beer
n_beer <- nrow(beer)
train_beer <- beer[1:(0.8*n_beer),]
test_beer <- beer[(0.8*n_beer+1):n_beer,]Setelah dataset dibagi menjadi data training dan testing, maka eksplorasi data training dapat dilakukan.
Eksplorasi Data
Eksplorasi data dilakukan untuk melihat bagaimana karakteristik dari data yang akan diolah. Selain itu, eksplorasi data juga memungkinkan untuk melihat apakah terdapat anomali pada data yang harus diperbaiki sebelum melakukan analisis lebih lanjut. Langkah pertama yang dapat dilakukan adalah melihat apakah terdapat missing-value pada data yang akan digunakan.
data.frame("Temperatur Minimum Harian" = sum(is.na(train_temp)),
"Produksi Listrik Bulanan" = sum(is.na(train_elec)),
"Produksi Beer Bulanan" = sum(is.na(train_beer))) %>%
knitr::kable(digits = 3, align = "c")| Temperatur.Minimum.Harian | Produksi.Listrik.Bulanan | Produksi.Beer.Bulanan |
|---|---|---|
| 0 | 0 | 0 |
Dari tabel tersebut terlihat bahwa pada data training tidak ada missing-value yang harus ditangani.
Plot dar masing-masing data dapat dilihat sebagai berikut.
windowsFonts(Times = windowsFont("Times New Roman"))
tsplot_temp <- train_temp %>% ggplot(aes(x = Date, y = Daily.minimum.temperatures)) +
geom_line(width = 0.75) +
geom_point(size = 0.75) +
labs(title = "Plot deret waktu suhu minimum harian \n") +
xlab("\nTanggal")+
ylab("Suhu minimum harian") +
scale_x_date(date_labels = "%m-%Y")+
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 11))
tsplot_elec <- train_elec %>%
ggplot(aes(x = DATE, y = IPG2211A2N)) +
geom_line(width = 0.75) +
geom_point(size = 0.75) +
labs(title = "Plot deret waktu produksi listrik bulanan \n") +
xlab("\nTanggal")+
ylab("Produski listrik bulanan") +
scale_x_date(date_labels = "%m-%Y")+
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 11))
tsplot_beer <- train_beer %>%
ggplot(aes(x = Month, y = Monthly.beer.production)) +
geom_line(width = 0.75) +
geom_point(size = 0.75) +
labs(title = "Plot deret waktu produksi beer bulanan \n") +
xlab("\nTanggal")+
ylab("Produski beer bulanan") +
scale_x_date(date_labels = "%m-%Y")+
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 11))
grid.arrange(tsplot_temp,
tsplot_elec,
tsplot_beer,
nrow = 3)
Dari plot time-series tersebut terlihat bahwa pada data suhu
minimum harian dan produksi beer terlihat ada pola musiman, sehingga
harus didiferensiasi terlebih dahulu.
train_temp_diff1 <- data.frame(
date = train_temp$Date[-1],
temp = diff(train_temp$Daily.minimum.temperatures,1)
)
train_beer_diff1 <- data.frame(
date = train_beer$Month[-1],
prod_beer = diff(train_beer$Monthly.beer.production,1)
)
#Plot suhu
tsplot_temp_diff1 <- train_temp_diff1 %>%
ggplot(aes(x = date, y = temp)) +
geom_line(width = 0.75) +
geom_point(size = 0.75) +
labs(title = "Plot pembedaan 1 deret waktu suhu minimum harian \n") +
xlab("\nTanggal")+
ylab("Suhu minimum harian") +
scale_x_date(date_labels = "%m-%Y")+
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 11))
tsplot_beer_diff1 <- train_beer_diff1 %>%
ggplot(aes(x = date, y = prod_beer)) +
geom_line(width = 0.75) +
geom_point(size = 0.75) +
labs(title = "Plot pembedaan 1 deret waktu produksi beer bulanan \n") +
xlab("\nTanggal")+
ylab("Produksi beer bulanan") +
scale_x_date(date_labels = "%m-%Y")+
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 11))
grid.arrange(tsplot_temp_diff1, tsplot_beer_diff1)
Terlihat bahwa setelah dilakukan pembedaan 1 untuk kedua data, maka data
tersebut sudah terlihat lebih stasioner. Untuk lebih memastikan apakah
seluruh data sudah stasioner atau belum dapat dilakukan dengan melihat
plot ACF, PACF, dan uji formal Augmented Dickey-Fuller.
Plot ACF dan PACF untuk masing-masing data dapat dilihat sebagai berikut.
acf_temp <- autoplot(acf(train_temp_diff1$temp, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Data Suhu Minimum Harian") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
acf_elec <- autoplot(acf(train_elec$IPG2211A2N, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Data Produksi Listrik Bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
acf_beer <- autoplot(acf(train_beer_diff1$prod_beer, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Data Produksi Beer Bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
pacf_temp <- autoplot(pacf(train_temp_diff1$temp, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Data Suhu Minimum Harian") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
pacf_elec <- autoplot(pacf(train_elec$IPG2211A2N, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Data Produksi Listrik Bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
pacf_beer <- autoplot(pacf(train_beer_diff1$prod_beer, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Data Produksi Beer Bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(acf_temp, pacf_temp,
acf_elec, pacf_elec,
acf_beer, pacf_beer)
Dari plot ACF dan PACF masing-masing model, terlihat bahwa seluruh data
sudah stsioner, sehingga untuk memastikan dapat dicek dengan uji
Augmented Dicky-Fuller. Hipotesis yang digunakan pada uji tersebut
adalah:
- \(H_0\): Data tidak stasioner
- \(H_1\): Data stasioner
t.stat <- c(adf.test(train_temp_diff1$temp)$statistic,
adf.test(train_elec$IPG2211A2N)$statistic,
adf.test(train_beer_diff1$prod_beer)$statistic)
t.pval <- c(adf.test(train_temp_diff1$temp)$p.value,
adf.test(train_elec$IPG2211A2N)$p.value,
adf.test(train_beer_diff1$prod_beer)$p.value)
t <- data.frame(t.stat, t.pval)
rownames(t) <- c("Temperatur Minimum Harian", "Produksi Listrik Bulanan",
"Produksi Beer Bulanan")
t %>% knitr::kable(digits = 3)| t.stat | t.pval | |
|---|---|---|
| Temperatur Minimum Harian | -10.051 | 0.01 |
| Produksi Listrik Bulanan | -3.803 | 0.02 |
| Produksi Beer Bulanan | -5.074 | 0.01 |
Karena seluruh data memiliki nilali p-value kurang dari 0.05, maka dapat dikatakan seluruh data sudah stasioner, sehingga dapat lanjut ke tahap selanjutnya.
Pendugaan Model Tentatif
Model dari data-data tersebut dapat diduga dengan melihat plot ACF, PACF, dan plot EACF. Untuk memudahkan maka masing-masing plot ACF dan PACF untuk masing-masing data akan dipanggil kembali beserta dengan plot EACF nya.
Data Suhu Minimum Harian
grid.arrange(acf_temp, pacf_temp, ncol = 2)
eacf(train_temp_diff1$temp)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o o o o o o o o o o o
## 1 x x x o o o o o o o o o o o
## 2 x x x o o o o o o o o o o o
## 3 x x x x o o o o o o o o o o
## 4 x o x o 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 o o o o o o o o o o
## 7 x x x x o x o o o o o o o oDari plot tersebut dapat diduga beberapa model untuk data tersebut:
- Karena ACF terlihat cuts-off di lag-2 dan pada PACF terlihat tails-off maka dapat diduga model mengikuti ARIMA(0,1,2).
- ACF terlihat cuts-off pada lag ke-1 dan PACF terlihat tails-off, maka dapat diduga model mengikuti ARIMA(0,1,1).
- Pada plot EACF, terlihat bahwa salah satu segitiga terbesar pertama yang unsur semuanya 0 memiliki sudut pada \(AR = 1\) dan \(MA = 3\) sehingga model ARIMA(1,1,3) dapat juga diduga sebagai model tersebut.
- Plot EACF juga memungkinkan model ARIMA(0,1,3)
Untuk menduga model terbaik untuk data suhu minimum harian adalah sebagai berikut. Pertama, cari model mana di antara keempat model tersebut yang signifikan untuk semua parameter. Penduga parameter untuk masing-masing model dapat dilihat dibawah.
mod1.suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(0,1,2), method = "ML")
mod2.suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(0,1,1), method = "ML")
mod3.suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(1,1,3), method = "ML")
mod4.suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(0,1,3), method = "ML")
lmtest::coeftest(mod1.suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.440262 0.048582 -9.0623 < 2.2e-16 ***
## ma2 -0.414031 0.046613 -8.8823 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod2.suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.790877 0.069962 -11.304 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod3.suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.354781 0.241677 -1.4680 0.14210
## ma1 -0.026544 0.230773 -0.1150 0.90843
## ma2 -0.506347 0.116886 -4.3320 1.478e-05 ***
## ma3 -0.274600 0.100996 -2.7189 0.00655 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod4.suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.376806 0.062847 -5.9956 2.027e-09 ***
## ma2 -0.369137 0.054981 -6.7139 1.895e-11 ***
## ma3 -0.113177 0.066645 -1.6982 0.08947 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hasil uji signifikansi parameter model ARIMA memperlihatkan bahwa model 1 (ARIMA(0,1,2)), model 2 (ARIMA(0,1,1)), dan ARIMA(0,1,3) yang signifikan untuk semua parameter. Maka ketiga model tersebut akan diinvestigasi untuk melihat model mana yang paling tepat untuk menggambarkan data tersebut.
Nilai AIC dan BIC dapat digunakan untuk melihat model mana di antara kedua model tersebut yang paling tepat mendeskripsikan data tersebut.
b <- data.frame(aic = c(mod1.suhu$aic, mod2.suhu$aic, mod4.suhu$aic),
bic = c(mod1.suhu$bic, mod2.suhu$bic, mod4.suhu$bic))
rownames(b) <- c("ARIMA(0,1,2)", "ARIMA(0,1,1)", "ARIMA(0,1,3)")
b %>% knitr::kable(digits = 3,
row.names = T,
col.names = c("AIC", "BIC"),
align = "c")| AIC | BIC | |
|---|---|---|
| ARIMA(0,1,2) | 1302.050 | 1313.070 |
| ARIMA(0,1,1) | 1357.396 | 1364.743 |
| ARIMA(0,1,3) | 1301.215 | 1315.909 |
Karena AIC dan BIC terkecil dimiliki model ARIMA(0,1,3), maka model yang paling tepat adalah model ARIMA(0,1,3).
Data Produksi Listrik Bulanan
Plot ACF, PACF, dan EACF data produksi listrik bulanan dapat dilihat sebagai berikut.
grid.arrange(acf_elec, pacf_elec)
eacf(train_elec$IPG2211A2N)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x o x x x o x o x x x o
## 1 x o x o x x x x x o x x x x
## 2 x x x x x x x x x x x x x x
## 3 x x o o o o o o o o x x x x
## 4 x x x o o o o o o o o x x x
## 5 x o x x o o o o o o o x o x
## 6 o o x x o o o o o o o x x o
## 7 o x o x x o o o o o o x x oDari plot tersebut dapat diduga beberapa model untuk data tersebut:
- ACF tails-off, sedangkan PACF cuts-off pada lag ke-2, sehingga model tentatif adalah ARIMA(2,0,0).
- Dari EACF juga terlihat bahwa ARIMA(4,0,3) memungkinkan untuk dipilih
- Dari EACF juga terlihat bahwa ARIMA(3,0,3) memungkinkan untuk dipilih
Untuk menduga model terbaik untuk data suhu minimum harian adalah sebagai berikut. Pertama, cari model mana di antara keempat model tersebut yang signifikan untuk semua parameter. Parameter dan uji signifikansi masing-masing parameter dapat dilihat di bawah ini.
mod1.elec <- Arima(train_elec$IPG2211A2N, order = c(2,0,0), method = "ML")
mod2.elec <- Arima(train_elec$IPG2211A2N, order = c(4,0,3), method = "ML")
mod3.elec <- Arima(train_elec$IPG2211A2N, order = c(3,0,3), method = "ML")
lmtest::coeftest(mod1.elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.029667 0.052371 19.661 < 2.2e-16 ***
## ar2 -0.695602 0.052048 -13.365 < 2.2e-16 ***
## intercept 97.193388 0.598926 162.280 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod2.elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.99123877 0.01588819 62.3884 < 2.2e-16 ***
## ar2 -0.00081599 NaN NaN NaN
## ar3 -0.99289497 0.01674712 -59.2875 < 2.2e-16 ***
## ar4 0.99666361 NaN NaN NaN
## ma1 -0.02001827 0.04303209 -0.4652 0.6418
## ma2 -0.46090806 NaN NaN NaN
## ma3 0.55717786 NaN NaN NaN
## intercept 97.85562904 24.58441885 3.9804 6.88e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod3.elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.948647 0.034275 56.8528 < 2.2e-16 ***
## ar2 -1.948420 0.034273 -56.8498 < 2.2e-16 ***
## ar3 0.949356 0.034567 27.4641 < 2.2e-16 ***
## ma1 -1.403077 0.106865 -13.1295 < 2.2e-16 ***
## ma2 1.351583 0.111702 12.0999 < 2.2e-16 ***
## ma3 -0.376015 0.121150 -3.1037 0.001911 **
## intercept 96.983985 2.527012 38.3789 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model dengan seluruh parameternya signifikan adalah model ARIMA(3,0,3) dan ARIMA(2,0,0) sehingga akan dipilih model yang paling tepat digunakan dengan melihat AIC dan BIC nya.
b <- data.frame(aic = c(mod1.elec$aic, mod3.elec$aic),
bic = c(mod1.elec$bic, mod3.elec$bic))
rownames(b) <- c("ARIMA(2,0,0)", "ARIMA(3,0,3)")
b %>% knitr::kable(digits = 3,
row.names = T,
col.names = c("AIC", "BIC"),
align = "c")| AIC | BIC | |
|---|---|---|
| ARIMA(2,0,0) | 1209.984 | 1223.014 |
| ARIMA(3,0,3) | 1035.772 | 1061.832 |
Dari hasil tersebut, terlihat bahwa model paling efektif adalah ARIMA(3,0,3).
Data Produksi Bir Bulanan
Plot ACF, PACF, dan EACF untuk data produksi listrik bulanan dapat dilihat sebagai berikut
grid.arrange(acf_beer, pacf_beer)
eacf(train_beer_diff1$prod_beer)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o x o x o o
## 1 x x o o o o o o o x o x o o
## 2 o x o o o x o o o o o x o o
## 3 o x o o o x o o o o o x x o
## 4 x x o o o o o o o o o x x o
## 5 x x x x o o o o o o o x o o
## 6 x o o x o o o o o o o x o o
## 7 x x o x o o o o o o o x o oTerlihat bahwa pola musiman sudah hilang, sehingga dari plot tersebut dapat diduga beberapa model untuk data tersebut:
- Plot PACF terlihat cuts-off saat lag 2 dan tails off di ACF. Sehingga dapat dikatakan model tersebut adalah ARIMA(2,1,0).
- Plot EACF memperlihatkan ARIMA(2, 1, 2).
- Plot EACF memperlihatkan ARIMA(1, 1, 2).
Untuk menduga model terbaik untuk data produksi beer bulanan adalah sebagai berikut. Pertama, cari model mana di antara ketiga model tersebut yang signifikan untuk semua parameter.
mod1.beer <- Arima(train_beer$Monthly.beer.production, order = c(2,1,0), method = "ML")
mod2.beer <- Arima(train_beer$Monthly.beer.production, order = c(2,1,2), method = "ML")
mod3.beer <- Arima(train_beer$Monthly.beer.production, order = c(1,1,5), method = "ML")
lmtest::coeftest(mod1.beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.30673 0.10116 -3.0322 0.002427 **
## ar2 -0.25323 0.10048 -2.5203 0.011724 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod2.beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -4.1253e-01 1.0372e-01 -3.9772 6.972e-05 ***
## ar2 2.2919e-01 1.0676e-01 2.1467 0.03182 *
## ma1 2.6532e-07 9.5215e-02 0.0000 1.00000
## ma2 -1.0000e+00 9.5215e-02 -10.5025 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod3.beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.158003 0.126434 -1.2497 0.2114149
## ma1 -0.246281 0.078502 -3.1372 0.0017055 **
## ma2 -0.632545 0.116562 -5.4267 5.741e-08 ***
## ma3 -0.256919 0.174737 -1.4703 0.1414761
## ma4 -0.329878 0.098464 -3.3502 0.0008074 ***
## ma5 0.540662 0.060744 8.9006 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hanya model ARIMA(2,1,0) yang memiliki parameter yang signifikan. Untuk itu, model ini dipilih karena satu-satunya model yang nyata di antara model-model lainnya.
b <- data.frame(aic = c(mod1.beer$aic),
bic = c(mod1.beer$bic))
rownames(b) <- c("ARIMA(2,1,0)")
b %>% knitr::kable(digits = 3,
row.names = T,
col.names = c("AIC", "BIC"),
align = "c")| AIC | BIC | |
|---|---|---|
| ARIMA(2,1,0) | 829.812 | 837.345 |
Dari hasil tersebut, terlihat bahwa model paling efektif adalah ARIMA(2,1,0).
Diagnostik Model
Diagnostik model yang akan dilakukan adalah melihat apakah sisaan menyebar normal, sisaan saling bebas atau tidak ada autokorelasi, dan apakah sisaan menyebar di sekitar nol.Diagnostik eksploratif
Uji Eksploratif
Data Suhu Minimum Harian
Model ARIMA yang digunakan untuk data ini adalah ARIMA(0,1,3). Plot sebaran sisaan dari model tersebut adalah sebagai berikut.
resid.suhu <- mod1.suhu$residuals
windowsFonts(Times = windowsFont("Times New Roman"))
#QQ Residual Plot
qq_suhu <- ggplot(resid.suhu, aes(sample = resid.suhu)) +
stat_qq(size = 0.75) +
stat_qq_line(width = 0.75) +
labs(title = "Normal Q-Q PLot") +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(-6, 7, 2)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#Plot sisaan
residplot_suhu <- ggplot(resid.suhu, aes(y = resid.suhu,
x = c(1:length(resid.suhu)))) +
geom_line(width = 0.75)+
geom_point(size = 0.75) +
geom_hline(yintercept = 0, col = "red") +
labs(title = "Plot Sisaan Model ARIMA(0,1,3)") +
xlab("Periode Waktu") +
ylab("Sisaan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#ACF plot
acf_suhu <- autoplot(acf(resid.suhu, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Sisaan Model ARIMA(0,1,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#PACF suhu
pacf_suhu <- autoplot(pacf(resid.suhu, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(0,1,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(qq_suhu, residplot_suhu, acf_suhu,pacf_suhu, nrow = 2)
Dari plot Normal Q-Q Plot, terlihat bahwa data relatif menyebar normal. Pada bagian tengah data terlihat dekat dengan garis 45 derajat. Akan tetapi, pada kedua ujung terlihat bahwa data mulai sedikit melencend dari garis 45 derajat terus. Sehingga peneyelidikan lanjutan harus dilakukan untuk menentukan apakah memang sisaan model tersebut menyebar normal. Plot time series sisaan model menunjukkan bahwa data sudah menyebar secara acak dan berada di sekitar garis 0, sehingga dapat dikatakan stasioner dan tidak ada autokorelasi antar sisaan. Plot ACF dan PACF juga mendukung hal tersebut.
Data Produksi Listrik Bulanan
Model ARIMA yang akan didiagnostik adalah model ARIMA (3,0,3). Plot sebaran sisaan model tersebut adalah sebagai berikut.
resid.elec <- mod3.elec$residuals
windowsFonts(Times = windowsFont("Times New Roman"))
#QQ Residual Plot
qq_elec <- ggplot(resid.elec, aes(sample = resid.elec)) +
stat_qq(size = 0.75) +
stat_qq_line(width = 0.75) +
labs(title = "Normal Q-Q PLot") +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(-6, 7, 2)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#Plot sisaan
residplot_elec <- ggplot(resid.elec, aes(y = resid.elec,
x = c(1:length(resid.elec)))) +
geom_line(width = 0.75)+
geom_point(size = 0.75) +
geom_hline(yintercept = 0, col = "red") +
labs(title = "Plot Sisaan Model ARIMA(3,0,3)") +
xlab("Periode Waktu") +
ylab("Sisaan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#ACF plot
acf_elec <-autoplot(acf(resid.elec, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(3,0,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#PACF suhu
pacf_elec <- autoplot(pacf(resid.elec, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(3,0,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(qq_elec, residplot_elec, acf_elec,pacf_elec, nrow = 2)
Untuk data produksi listrik bulanan, terlihat bahwa secara keseluruhan sisaan menyebar normal, akan tetapi sama seperti sebelumnya kedua ujung data menyimpang relatif jauh sehingga harus dilakukan uji formal untuk memastikan apakah benar-benar sisaan menyebar normal. Pada plot sisaan model terlihat juga bahwa ada beberapa titik, terutama di bagian akhir data, yang menyimpang agak jauh dari garis 0. Hal ini membuat asumsi tidak adanya autokorelasi pada data diragukan, sehingga harus melakukan uji lebih lanjut. Akan tetapi, pada plot ACF dan PACF terlihat bahwa korelasi antar lag data relatif tidak nyata, kecuali untuk lag ke 6 dan ke 11, sehingga asumsi autokorelasi harus dicek kembali dengan uji formal.
Data Produski Beer Bulanan
Model ARIMA yang akan didiagnostik adalah model ARIMA (2,1,0). Plot sebaran sisaan model tersebut adalah sebagai berikut.
resid.beer <- mod1.beer$residuals
windowsFonts(Times = windowsFont("Times New Roman"))
par(mfrow = c(2,2))
#QQ Residual Plot
qq_beer <- ggplot(resid.beer, aes(sample = resid.beer)) +
stat_qq(size = 0.75) +
stat_qq_line(width = 0.75) +
labs(title = "Normal Q-Q PLot") +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(-6, 7, 2)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#Plot sisaan
residplot_beer <- ggplot(resid.beer, aes(y = resid.beer,
x = c(1:length(resid.beer)))) +
geom_line(width = 0.75)+
geom_point(size = 0.75) +
geom_hline(yintercept = 0, col = "red") +
labs(title = "Plot Sisaan Model ARIMA(2,1,0)") +
xlab("Periode Waktu") +
ylab("Sisaan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#ACF plot
acf_beer <-autoplot(acf(resid.beer, plot = F)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Sisaan Model ARIMA(2,1,0)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#PACF suhu
pacf_beer <- autoplot(pacf(resid.beer, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(2,1,0)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(qq_beer, residplot_beer, acf_beer,pacf_beer)
Pada model ARIMA(2,1,0) untuk data produksi bir bulanan, terlihat bahwa sisaan menyebar normal karena mayoritas sisaan berhimpit dengan garis 45 derajat. Kemudian pada plot sisaan model terlihat bahwa sisaan menyebar secara acak di sekitar 0, mesikpun terlihat sedikit pola musiman sehingga asumsi tidak ada autokorelasi mungkin belum terpenuhi. Pada plot ACF dan PACF tidak memperlihatkan ada pola musiman yang terbentuk, tetapi terdapat kejanggalan pada plot ACF lag 12 dan PACF lag 10 dimana korelasi antara data asli dengan masing-masing lag signifikan, sehingga harus diinvestigasi lebih lanjut.
Uji Formal
Uji untuk mengecek apakah benar sisaan tersebut menyebar normal, tidak ada autokorelasi, dan sisaan menyebar di sekitar 0 dapat dilakukan dengan beberapa uji. Untuk menguji apakah data menyebar normal dapat menggunakan uji Shapiro-Wilk. Nilai \(\alpha\) yang digunakan adalah 5\%. Hipotesis yang digunakan adalah sebagai berikut:
- \(H_0\): data menyebar normal
- \(H_1\): data tidak menyebar normal
norm.suhu <- shapiro.test(resid.suhu)
norm.elec <- shapiro.test(resid.elec)
norm.beer <- shapiro.test(resid.beer)
norm.test <- data.frame("Statistik Uji" = c(norm.suhu$statistic,
norm.elec$statistic,
norm.beer$statistic),
"P-value" = c(norm.suhu$p.value,
norm.elec$p.value,
norm.beer$p.value))
rownames(norm.test) <- c("Sisaan suhu minimum harian",
"Sisaan produksi listrik bulanan",
"Sisaan produksi beer bulanan")
norm.test %>% knitr::kable(digits = 3, row.names = T,
col.names = c("Statistik Uji",
"p-value"))| Statistik Uji | p-value | |
|---|---|---|
| Sisaan suhu minimum harian | 0.993 | 0.233 |
| Sisaan produksi listrik bulanan | 0.993 | 0.485 |
| Sisaan produksi beer bulanan | 0.985 | 0.371 |
Karena seluruh nilai p-value > \(\alpha\), maka dapati disimpulkan tak tolak \(H_0\). Data sisaan untuk semua model untuk masing-masing data menyebar dengan normal.
Untuk menguji apakah terdapat autokorelasi pada sebaran sisaan model, dapat menggunakan uji Box-Ljung. Hipotesis yang akan digunakan adalah
- \(H_0\): Tidak ada autokorelasi
- \(H_1\): Terdapat autokorelasi
autocor_suhu <- Box.test(resid.suhu, type = "Ljung-Box")
autocor_elec <- Box.test(resid.elec, type = "Ljung-Box")
autocor_beer <- Box.test(resid.beer, type = "Ljung-Box")
autocor <- data.frame(stat = c(autocor_suhu$statistic,
autocor_elec$statistic,
autocor_beer$statistic),
pval = c(autocor_suhu$p.value,
autocor_elec$p.value,
autocor_beer$p.value))
rownames(autocor) <- c("Sisaan suhu minimum harian",
"Sisaan produksi listrik bulanan",
"Sisaan produksi beer bulanan")
autocor %>% knitr::kable(digits = 3, row.names = T,
col.names = c("Statistik Uji",
"P-Value"))| Statistik Uji | P-Value | |
|---|---|---|
| Sisaan suhu minimum harian | 0.295 | 0.587 |
| Sisaan produksi listrik bulanan | 0.795 | 0.373 |
| Sisaan produksi beer bulanan | 0.001 | 0.970 |
Karena seluruh model untuk masing-masing data memiliki nili p-value > \(\alpha\) maka dapat disimpulkan bahwa sisaan dari seluruh data dengan masing-masing model saling bebas dan tidak ada autokorelasi pada data tersebut.
Untuk memeriksa apakah sisaan masing-masing model untuk masing-masing data bergerak di sekitar 0, maka dapat menggunakan uji-t untuk melihat apakah rata-rata dari sisaan berbeda secara nyata dari 0.
ttest_temp <- t.test(resid.suhu, mu = 0, conf.level = 0.95)
ttest_elec <- t.test(resid.elec, mu = 0, conf.level = 0.95)
ttest_beer <- t.test(resid.beer, mu = 0, conf.level = 0.95)
ttest <- data.frame(stat = c(ttest_temp$statistic, ttest_elec$statistic,
ttest_beer$statistic),
pval = c(ttest_temp$p.value, ttest_elec$p.value,
ttest_beer$p.value))
rownames(ttest) <- c("Suhu Minimum Harian",
"Produksi Listrik Bulanan",
"Produksi Beer Bulanan")
ttest %>% knitr::kable(digits = 3, row.names = T,
col.names = c("Statistik", "P-value"))| Statistik | P-value | |
|---|---|---|
| Suhu Minimum Harian | -0.793 | 0.428 |
| Produksi Listrik Bulanan | 0.591 | 0.555 |
| Produksi Beer Bulanan | -0.178 | 0.859 |
Karena hasil uji-t memperlihatkan nilai p-value untuk semua sisaan data lebih dari 0.05, maka dapat dikatakan bahwa sisaan model menyebar di sekitar 0.
Overfitting untuk menentukan model terbaik
Data Suhu Minimum Harian
Overfitting adalah proses meningkatkan salah satu parameter dari model tentatif untuk melihat apakah terdapat model yang lebih baik di sekitar dugaan model tentatif yang telah ditentukan. Model pertama yang akan dilihat adalah model ARIMA(0,1,3). Model overfit yang dapat diuji adalah model ARIMA(0,1,4) dan ARIMA(0,2,3).
overfit_mod1_suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(0,1,4), method = "ML")
overfit_mod2_suhu <- Arima(train_temp$Daily.minimum.temperatures, order = c(0,2,3), method = "ML")
lmtest::coeftest(overfit_mod1_suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.406540 0.058671 -6.9292 4.233e-12 ***
## ma2 -0.401858 0.063084 -6.3702 1.888e-10 ***
## ma3 -0.163301 0.066017 -2.4736 0.01337 *
## ma4 0.120848 0.057233 2.1115 0.03473 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit_mod2_suhu)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.447737 0.048306 -29.9702 <2e-16 ***
## ma2 0.021110 0.088679 0.2381 0.8118
## ma3 0.430954 0.047528 9.0674 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dari 2 model overfit yang ditetapkan, hanya model ARIMA(0,1,4) yang memiliki parameter signifikan, sehingga model tersebut akan dibandingkan dengan model ARIMA(0,1,3).
crit_temp <- data.frame(aic = c(mod4.suhu$aic, overfit_mod1_suhu$aic),
bic = c(mod4.suhu$bic, overfit_mod1_suhu$bic))
rownames(crit_temp) <- c("ARIMA(0,1,3)", "ARIMA(0,1,4)")
crit_temp %>% knitr::kable(digits = 3, row.names = T,
col.names = c("AIC", "BIC"))| AIC | BIC | |
|---|---|---|
| ARIMA(0,1,3) | 1301.215 | 1315.909 |
| ARIMA(0,1,4) | 1298.810 | 1317.177 |
Karena AIC dari model overfit lebih kecil dibandingkan dengan model aslinya, maka model tentatif yang akan digunakan untuk peramalan adalah model ARIMA(0,1,4). Diagnostik eksploratif model tersebut dapat dilihat sebagai berikut.
resid.suhu.overfit <- overfit_mod1_suhu$residuals
windowsFonts(Times = windowsFont("Times New Roman"))
#QQ Residual Plot
qq_suhu_overfit <- ggplot(resid.suhu.overfit, aes(sample = resid.suhu.overfit)) +
stat_qq(size = 0.75) +
stat_qq_line(width = 0.75) +
labs(title = "Normal Q-Q PLot") +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(-6, 7, 2)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#Plot sisaan
residplot_suhu_overfit <- ggplot(resid.suhu.overfit, aes(y = resid.suhu.overfit,
x = c(1:length(resid.suhu.overfit)))) +
geom_line(width = 0.75)+
geom_point(size = 0.75) +
geom_hline(yintercept = 0, col = "red") +
labs(title = "Plot Sisaan Model ARIMA(0,1,4)") +
xlab("Periode Waktu") +
ylab("Sisaan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#ACF plot
acf_suhu_overfit <- autoplot(acf(resid.suhu.overfit, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Sisaan Model ARIMA(0,1,4)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#PACF suhu
pacf_suhu_overfit <- autoplot(pacf(resid.suhu.overfit, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(0,1,4)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(qq_suhu_overfit, residplot_suhu_overfit, acf_suhu_overfit, pacf_suhu_overfit, nrow = 2)
Dari plot Normal Q-Q Plot, terlihat bahwa sisaan model overfit
relatif menyebar normal. Pada bagian tengah sisaan data terlihat
bersinggungan dengan garis 45 derajat. Akan tetapi, pada kedua ujung
terlihat bahwa data mulai sedikit melenceng dari garis 45 derajat.
Sehingga peneyelidikan lanjutan harus dilakukan untuk menentukan apakah
memang sisaan model tersebut menyebar normal. Plot time series sisaan
model menunjukkan bahwa data sudah menyebar secara acak dan berada di
sekitar garis 0, sehingga dapat dikatakan stasioner dan tidak ada
autokorelasi antar sisaan. Plot ACF dan PACF juga mendukung hal
tersebut.
Untuk uji formal diagnostik model dapat dilihat sebagai berikut.
data.frame(uji = c("Uji kenormalan sisaan", "Uji autokorelasi sisaan",
"Uji nilai tengah sisaan"),
stat = c(shapiro.test(resid.suhu.overfit)$statistic,
Box.test(resid.suhu.overfit, type = "Ljung-Box")$statistic,
t.test(resid.suhu.overfit, mu = 0, conf.level = 0.95)$statistic),
p.val =c(shapiro.test(resid.suhu.overfit)$p.val,
Box.test(resid.suhu.overfit, type = "Ljung-Box")$p.val,
t.test(resid.suhu.overfit, mu = 0, conf.level = 0.95)$p.val)
) %>%
knitr::kable(row.names = F,
col.names = c("Jenis Uji", "Statistik Uji",
"P Value"))| Jenis Uji | Statistik Uji | P Value |
|---|---|---|
| Uji kenormalan sisaan | 0.9951972 | 0.5009195 |
| Uji autokorelasi sisaan | 0.0038747 | 0.9503659 |
| Uji nilai tengah sisaan | -0.7785367 | 0.4368852 |
Dari hasil uji formal, terlihat bahwa untuk uji kenormalan, nilai p-value > 0.05, sehingga dapat dikatakan bahwa sisaan data sudah menyebar normal. Kemudian untuk uji autokorelasi, nilai p-value > 0.05, sehingga sisaan model tersebut sudah saling bebas. Hasil uji t juga memperlihatkan nilai p-value > 0.05 sehingga data sudah terbukti menyebar di sekitar 0.
Data Produksi Listrik Bulanan
Model tentatif yang ditetapkan sebelumnya adalah model ARIMA(3,0,3). Model overfit yang dapat diteliti adalah model ARIMA(4,0,3), ARIMA(3,1,3), dan ARIMA(3,0,4).
overfit_mod1_elec <- Arima(train_elec$IPG2211A2N, order = c(4,0,3), method = "ML")
overfit_mod2_elec <- Arima(train_elec$IPG2211A2N, order = c(3,1,3), method = "ML")
overfit_mod3_elec <- Arima(train_elec$IPG2211A2N, order = c(3,0,4), method = "ML")
lmtest::coeftest(overfit_mod1_elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.99123877 0.01588819 62.3884 < 2.2e-16 ***
## ar2 -0.00081599 NaN NaN NaN
## ar3 -0.99289497 0.01674712 -59.2875 < 2.2e-16 ***
## ar4 0.99666361 NaN NaN NaN
## ma1 -0.02001827 0.04303209 -0.4652 0.6418
## ma2 -0.46090806 NaN NaN NaN
## ma3 0.55717786 NaN NaN NaN
## intercept 97.85562904 24.58441885 3.9804 6.88e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit_mod2_elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.462389 0.079331 18.4339 < 2.2e-16 ***
## ar2 -1.462511 0.079246 -18.4553 < 2.2e-16 ***
## ar3 0.463074 0.079344 5.8363 5.337e-09 ***
## ma1 -1.948752 0.048245 -40.3932 < 2.2e-16 ***
## ma2 1.943657 0.098288 19.7752 < 2.2e-16 ***
## ma3 -0.900815 0.064146 -14.0431 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit_mod3_elec)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.981306 0.015692 126.2626 < 2.2e-16 ***
## ar2 -1.981337 0.015619 -126.8534 < 2.2e-16 ***
## ar3 0.984122 0.015683 62.7491 < 2.2e-16 ***
## ma1 -1.433721 0.089684 -15.9863 < 2.2e-16 ***
## ma2 1.185243 0.147342 8.0441 8.685e-16 ***
## ma3 -0.295850 0.137680 -2.1488 0.03165 *
## ma4 -0.189419 0.082132 -2.3063 0.02110 *
## intercept 97.683529 3.293843 29.6564 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model overfit dengan seluruh parameter signifikan adalah model ARIMA(3,1,3) dan ARIMA(3,0,4). Kedua model tersebut akan dibandingkan dengan model tentatif yang telah ditentukan sebelumnya.
crit_elec <- data.frame(aic = c(mod3.elec$aic, overfit_mod2_elec$aic, overfit_mod3_elec$aic),
bic = c(mod3.elec$bic, overfit_mod2_elec$bic, overfit_mod3_elec$bic))
rownames(crit_elec) <- c("ARIMA(3,0,3)", "ARIMA(3,1,3)", "ARIMA(3,0,4)")
crit_elec %>% knitr::kable(digits = 3, row.names = T,
col.names = c("AIC", "BIC"))| AIC | BIC | |
|---|---|---|
| ARIMA(3,0,3) | 1035.772 | 1061.832 |
| ARIMA(3,1,3) | 1012.259 | 1035.025 |
| ARIMA(3,0,4) | 1037.724 | 1067.041 |
Karena model ARIMA(3,1,3) memiliki nilai AIC dan BIC terkecil, maka model yang akan digunakan untuk peramalan data adalah model ARIMA(3,1,3). Diagnostik eksploratif model tersebut dapat dilihat sebagai berikut,
resid.elec.overfit <- overfit_mod2_elec$residuals
windowsFonts(Times = windowsFont("Times New Roman"))
#QQ Residual Plot
qq_elec_overfit <- ggplot(resid.elec.overfit, aes(sample = resid.elec.overfit)) +
stat_qq(size = 0.75) +
stat_qq_line(width = 0.75) +
labs(title = "Normal Q-Q PLot") +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(-6, 7, 2)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#Plot sisaan
residplot_elec_overfit <- ggplot(resid.elec.overfit, aes(y = resid.elec.overfit,
x = c(1:length(resid.elec.overfit)))) +
geom_line(width = 0.75)+
geom_point(size = 0.75) +
geom_hline(yintercept = 0, col = "red") +
labs(title = "Plot Sisaan Model ARIMA(3,1,3)") +
xlab("Periode Waktu") +
ylab("Sisaan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#ACF plot
acf_elec_overfit <- autoplot(acf(resid.elec.overfit, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot ACF Sisaan Model ARIMA(3,1,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
#PACF suhu
pacf_elec_overfit <- autoplot(pacf(resid.elec.overfit, plot = FALSE)) +
geom_hline(yintercept = 0, col = "black") +
labs(title = "Plot PACF Sisaan Model ARIMA(3,1,3)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(qq_elec_overfit, residplot_elec_overfit, acf_elec_overfit, pacf_elec_overfit, nrow = 2)
Pada plot Q-Q, terlihat bahwa data relatif melenceng dari garis 45 derajat, sehingga harus dilakukan uji formal. Untuk asumsi kebebasan, plot sisaan, acf, dan pacf sudah dapat memperlihatkan bahwa asumsi-asumsi tersebut sudah terpenuhi. Akan tetapi, untuk asumsi persebaran di sekitar 0, plot sisaan memperlihatkan bahwa data relatif menyebar di atas 0 sehingga harus dilakukan uji formal. Untuk pengujian formal dapat dilihat sebagai berikut.
data.frame(uji = c("Uji kenormalan sisaan", "Uji autokorelasi sisaan",
"Uji nilai tengah sisaan"),
stat = c(shapiro.test(resid.elec.overfit)$statistic,
Box.test(resid.elec.overfit, type = "Ljung-Box")$statistic,
t.test(resid.elec.overfit, mu = 0, conf.level = 0.95)$statistic),
p.val =c(shapiro.test(resid.elec.overfit)$p.val,
Box.test(resid.elec.overfit, type = "Ljung-Box")$p.val,
t.test(resid.elec.overfit, mu = 0, conf.level = 0.95)$p.val)
) %>%
knitr::kable(row.names = F,
col.names = c("Jenis Uji", "Statistik Uji",
"P Value"))| Jenis Uji | Statistik Uji | P Value |
|---|---|---|
| Uji kenormalan sisaan | 0.9912713 | 0.2998663 |
| Uji autokorelasi sisaan | 0.0056082 | 0.9403040 |
| Uji nilai tengah sisaan | 2.0256970 | 0.0441872 |
Dari hasil uji formal terlihat bahwa sisaan model sudah menyebar normal dan sudah saling bebas karena pada uji normalitas dan uji kebebeasan memiliki p-value > 0.05. Akan tetapi, untuk uji persebaran data di sekitar 0 terlihat bahwa p-value < 0.05 sehingga dapat dikatakan sisaan tidak menyebar di sekitar 0. Karena pelanggaran asumsi ini, model yang akan digunakan untuk melakukan peramalan adalah model ARIMA(3,0,3).
Data Produksi Beer Bulanan
Model tentatif yang telah ditentukan sebelumnya adalah model ARIMA(2,1,0). Model overfit yang dapat diuji kelayakannya adalah model ARIMA(3,1,0), ARIMA(2,2,0), dan ARIMA(2,1,1).
overfit_mod1_beer <- Arima(train_beer$Monthly.beer.production, order = c(3,1,0), method = "ML")
overfit_mod2_beer <- Arima(train_beer$Monthly.beer.production, order = c(2,2,0), method = "ML")
overfit_mod3_beer <- Arima(train_beer$Monthly.beer.production, order = c(2,1,1), method = "ML")
lmtest::coeftest(overfit_mod1_beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.3074340 0.1046441 -2.9379 0.003304 **
## ar2 -0.2540663 0.1054632 -2.4091 0.015994 *
## ar3 -0.0027093 0.1039090 -0.0261 0.979199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit_mod2_beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.788385 0.090299 -8.7308 < 2.2e-16 ***
## ar2 -0.503485 0.089704 -5.6127 1.992e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit_mod3_beer)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.443318 0.104465 4.2437 2.199e-05 ***
## ar2 -0.068131 0.104939 -0.6492 0.5162
## ma1 -0.999999 0.043411 -23.0356 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hanya model ARIMA(2,2,0) yang seluruh parameternya signifikan. Sehingga model ini akan dibandingkan dengan model tentatif ARIMA(2,1,0).
crit_beer <- data.frame(aic = c(mod1.beer$aic, overfit_mod2_beer$aic),
bic = c(mod1.beer$bic, overfit_mod2_beer$bic))
rownames(crit_beer) <- c("ARIMA(2,1,0)", "ARIMA(2,2,0)")
crit_beer %>% knitr::kable(digits = 3, row.names = T,
col.names = c("AIC", "BIC"))| AIC | BIC | |
|---|---|---|
| ARIMA(2,1,0) | 829.812 | 837.345 |
| ARIMA(2,2,0) | 859.863 | 867.363 |
Uniknya, AIC dan BIC untuk kedua model tersebut sama persis, sehingga model tentatif awal ARIMA(2,1,0) masih layak untuk digunakan.
Peramalan data
Data peramalan dan plot peramalan dapat dilihat sebagai berikut.
forecast_temp <- forecast(overfit_mod1_suhu, h = nrow(test_temp))
forecast_elec <- forecast(mod3.elec, h = nrow(test_elec))
forecast_beer <- forecast(mod1.beer, h = nrow(test_beer))
comb_forecast_temp <- head(cbind(test_temp, forecast_temp$mean),10) %>%
knitr::kable(digits = 3,
col.names = c("Date", "Suhu minimum harian",
"Data forecasting suhu"),
row.names = F)
comb_forecast_elec <- head(cbind(test_elec, forecast_elec$mean), 10) %>%
knitr::kable(digits = 3,
col.names = c("Date", "Produksi listrik bulanan",
"Data forecasting produksi listrik bulanan"),
row.names = F)
comb_forecast_beer <- head(cbind(test_beer, forecast_beer$mean), 10) %>%
knitr::kable(digits = 3,
col.names = c("Date", "Produksi beer bulanan",
"Data forecasting produksi beer bulanan"),
row.names = F)
comb_forecast_temp| Date | Suhu minimum harian | Data forecasting suhu |
|---|---|---|
| 1990-10-20 | 6.4 | 10.601 |
| 1990-10-21 | 10.9 | 10.415 |
| 1990-10-22 | 9.0 | 10.748 |
| 1990-10-23 | 10.9 | 10.677 |
| 1990-10-24 | 12.4 | 10.677 |
| 1990-10-25 | 11.6 | 10.677 |
| 1990-10-26 | 13.3 | 10.677 |
| 1990-10-27 | 14.4 | 10.677 |
| 1990-10-28 | 18.4 | 10.677 |
| 1990-10-29 | 13.6 | 10.677 |
comb_forecast_elec| Date | Produksi listrik bulanan | Data forecasting produksi listrik bulanan |
|---|---|---|
| 2014-01-01 | 124.255 | 118.420 |
| 2014-02-01 | 112.881 | 114.030 |
| 2014-03-01 | 104.763 | 101.760 |
| 2014-04-01 | 90.287 | 93.428 |
| 2014-05-01 | 92.134 | 96.932 |
| 2014-06-01 | 101.878 | 108.345 |
| 2014-07-01 | 108.550 | 115.849 |
| 2014-08-01 | 108.194 | 111.559 |
| 2014-09-01 | 100.417 | 99.414 |
| 2014-10-01 | 92.384 | 91.231 |
comb_forecast_beer| Date | Produksi beer bulanan | Data forecasting produksi beer bulanan |
|---|---|---|
| 1993-09-01 | 143 | 136.572 |
| 1993-10-01 | 151 | 134.585 |
| 1993-11-01 | 177 | 136.063 |
| 1993-12-01 | 184 | 136.113 |
| 1994-01-01 | 151 | 135.723 |
| 1994-02-01 | 134 | 135.830 |
| 1994-03-01 | 164 | 135.896 |
| 1994-04-01 | 126 | 135.849 |
| 1994-05-01 | 131 | 135.846 |
| 1994-06-01 | 125 | 135.859 |
forecastplot_temp <- autoplot(forecast(overfit_mod1_suhu, h = nrow(test_temp))) +
xlab("Periode waktu") +
ylab("Suhu minimum harian") +
labs(title = "Forecasting Suhu Minimum Harian") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
forecastplot_elec <- autoplot(forecast_elec) +
xlab("Periode waktu") +
ylab("Produksi listrik bulanan") +
labs(title = "Forecasting Produksi Listrik bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
forecastplot_beer <- autoplot(forecast_beer) +
xlab("Periode waktu") +
ylab("Produksi beer bulanan") +
labs(title = "Forecasting Produksi Beer bulanan") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
text = element_text(family = "Times", size = 9))
grid.arrange(forecastplot_temp, forecastplot_elec, forecastplot_beer)
Dari plot forecasting tersebut terlihat bahwa data forecast sudah cukup relatif mengikuti pola data training. Akurasi dari data dapat dilihat sebagai berikut.
forecast_temp_data <- forecast_temp$mean
forecast_elec_data <- forecast_elec$mean
forecast_beer_data <- forecast_beer$mean
acc_temp <- accuracy(forecast_temp_data,
test_temp$Daily.minimum.temperatures)
acc_elec <- accuracy(forecast_elec_data,
test_elec$IPG2211A2N)
acc_beer <- accuracy(forecast_beer_data,
test_beer$Monthly.beer.production)
acc_tot <- rbind(acc_temp, acc_elec, acc_beer)
rownames(acc_tot) <- c("Suhu Minimum Harian",
"Produksi Listrik Bulanan",
"Produksi Beer Bulanan")
acc_tot %>% knitr::kable(
digits = 3,
row.names = T
)| ME | RMSE | MAE | MPE | MAPE | |
|---|---|---|---|---|---|
| Suhu Minimum Harian | 2.692 | 3.836 | 3.125 | 16.489 | 21.881 |
| Produksi Listrik Bulanan | 2.235 | 4.893 | 3.934 | 2.064 | 3.754 |
| Produksi Beer Bulanan | 11.285 | 23.089 | 17.181 | 6.037 | 10.748 |
Dari data akurasi tersebut terlihat bahwa galat untuk masing-masing model untuk masing-masing data sudah relatif kecil, sehingga dugaan model tentatif sudah cukup baik untuk menggambarkan masing-masing data.
Kesimpulan
Pada bentuk aslinya, data temperatur minimum harian tidak stasioner dan produksi beer bulanan terlihat belum stasioner, sedangkan untuk data produksi listrik bulanan sudah stasioner. Untuk menangani ketidakstasioneran data temperatur minimimum harian dan produksi beer bulanan dilakukan operasi first-differencing. Dengan proses tersebut sudah cukup untuk menstasionerkan kedua data tak stasioner tersebut.
Pemilihan model terbaik untuk masing-masing data dilakukan dengan menentukan parameter dugaan dari plot ACF, PACF, dan EACF. Kemudian mencari dan menguji signifikansi parameter masing-masing model lalu menghitung AIC dan BIC untuk masing-masing model. Kemudian diakhir dengan overfitting model. Model ARIMA(0,1,4) paling tepat digunakan untuk data suhu minimum harian dibandingkan dengan dugaan model lainnya untuk data tersebut. Model ARIMA(3,0,3) satu-satunya dugaan model yang signifikan pada data produksi listrik bulanan, sehingga model tersebut paling tepat dibandingkan dugaan model lainnya. Model ARIMA(2,1,0) adalah dugaan model paling tepat untuk data produksi bir bulanan. Model-model tersebut menghasilkan peramalan data yang sudah baik dengan galat yang relatif kecil.