R Notebook
Library
library("forecast")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary("graphics")
library("TTR")
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(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forecast)
library(lmtest) ## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt)
library(HoRM)
library(corrplot)## corrplot 0.92 loadedlibrary(dLagM)## Loading required package: nardl## Loading required package: dynlm##
## Attaching package: 'dLagM'## The following object is masked from 'package:forecast':
##
## forecastlibrary(dynlm)
library(MLmetrics)##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:dLagM':
##
## MAPE## The following object is masked from 'package:base':
##
## Recalllibrary(ggplot2)
library(tsibble)##
## Attaching package: 'tsibble'## The following object is masked from 'package:zoo':
##
## index## The following objects are masked from 'package:base':
##
## intersect, setdiff, unionlibrary(tseries)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(aTSA)##
## Attaching package: 'aTSA'## The following objects are masked from 'package:tseries':
##
## adf.test, kpss.test, pp.test## The following object is masked from 'package:dLagM':
##
## forecast## The following object is masked from 'package:forecast':
##
## forecast## The following object is masked from 'package:graphics':
##
## identifyInput data
dataziq <- read.csv("https://raw.githubusercontent.com/Raziqizzan03/MPDW/main/data/Pertemuan1/data%20suhu%20mekkah.csv")
dataziqUbah data menjadi
dataziq.ts <- ts(dataziq$suhu.mekkah)Membagi data menjadi data latih (80%) dan data uji (20%) dan jadikan time series
dt.latih <- dataziq[1:115,2]
dt.uji <- dataziq[116:144,2]
dt.latih.ts <- ts(dt.latih, frequency = 12)
dt.uji.ts <- ts(dt.uji, frequency = 12)Eksplorasi data
Membuat plot time series
# data full
ts.plot(dataziq.ts, col="purple",main="Plot semua data", type="o",xlab="Time",ylab="Data")
points(dataziq.ts)
# data latih
ts.plot(dt.latih.ts, col="orange",main="Plot data latih", type="o",xlab="Time",ylab="Data")
points(dt.latih.ts)
# data uji
ts.plot(dt.uji.ts, col="green",main="Plot data uji", type="o",xlab="Time",ylab="Data")
points(dt.uji.ts)
dapat dilihat bahwwa data deret waktu pada suhu mekkah membentuk pola aditif musiman tanpa trend
Pemulusan
Single Moving Average
Pemulusan dan peramalan
data.sma <- SMA(dt.latih.ts,n=4)
data.sma## Jan Feb Mar Apr May Jun Jul Aug Sep
## 1 NA NA NA 27.6800 30.6950 33.9625 36.3200 37.7250 37.9125
## 2 27.0050 24.5875 24.0525 25.6675 29.2750 32.4350 35.2625 36.8625 36.9875
## 3 26.3850 24.5475 24.7350 26.4400 29.9475 32.9775 35.6900 37.0925 37.0025
## 4 27.0450 25.5050 25.7875 27.3300 30.3350 32.9200 35.0700 36.5250 36.9550
## 5 26.9325 24.9675 24.8550 27.1250 30.1350 33.3425 35.8125 37.1150 37.4500
## 6 26.7575 24.8875 25.4150 26.9450 30.4900 33.4125 35.4225 37.3950 37.7525
## 7 26.8425 24.9475 25.3075 27.3550 30.8150 33.9550 35.8750 37.3675 37.5850
## 8 27.6825 25.2000 24.9275 26.8250 29.8750 33.8225 36.4900 38.1900 38.5375
## 9 27.2625 25.3775 25.6625 27.3150 30.7225 33.6375 35.6275 37.0650 37.4750
## 10 26.6750 24.6325 24.5550 26.2225 29.2050 32.8600 35.9250
## Oct Nov Dec
## 1 36.8550 34.5625 30.9650
## 2 35.9950 33.0050 29.8625
## 3 35.6475 33.2500 30.1375
## 4 36.0225 33.6125 30.4325
## 5 36.2000 33.5575 30.3850
## 6 36.7825 34.6025 30.6250
## 7 36.3450 34.0525 30.9475
## 8 37.1675 34.7950 31.1050
## 9 36.1725 33.2475 30.0550
## 10data.ramal<-c(NA,data.sma)
data.ramal## [1] NA NA NA NA 27.6800 30.6950 33.9625 36.3200 37.7250
## [10] 37.9125 36.8550 34.5625 30.9650 27.0050 24.5875 24.0525 25.6675 29.2750
## [19] 32.4350 35.2625 36.8625 36.9875 35.9950 33.0050 29.8625 26.3850 24.5475
## [28] 24.7350 26.4400 29.9475 32.9775 35.6900 37.0925 37.0025 35.6475 33.2500
## [37] 30.1375 27.0450 25.5050 25.7875 27.3300 30.3350 32.9200 35.0700 36.5250
## [46] 36.9550 36.0225 33.6125 30.4325 26.9325 24.9675 24.8550 27.1250 30.1350
## [55] 33.3425 35.8125 37.1150 37.4500 36.2000 33.5575 30.3850 26.7575 24.8875
## [64] 25.4150 26.9450 30.4900 33.4125 35.4225 37.3950 37.7525 36.7825 34.6025
## [73] 30.6250 26.8425 24.9475 25.3075 27.3550 30.8150 33.9550 35.8750 37.3675
## [82] 37.5850 36.3450 34.0525 30.9475 27.6825 25.2000 24.9275 26.8250 29.8750
## [91] 33.8225 36.4900 38.1900 38.5375 37.1675 34.7950 31.1050 27.2625 25.3775
## [100] 25.6625 27.3150 30.7225 33.6375 35.6275 37.0650 37.4750 36.1725 33.2475
## [109] 30.0550 26.6750 24.6325 24.5550 26.2225 29.2050 32.8600 35.9250# Forecasting 29 periode
data.gab<-cbind(aktual=c(dt.latih.ts,rep(NA,29)),pemulusan=c(data.sma,rep(NA,29)),ramalan=c(data.ramal,rep(data.ramal[length(data.ramal)],28)))
data.gab## aktual pemulusan ramalan
## [1,] 24.10 NA NA
## [2,] 25.12 NA NA
## [3,] 28.51 NA NA
## [4,] 32.99 27.6800 NA
## [5,] 36.16 30.6950 27.6800
## [6,] 38.19 33.9625 30.6950
## [7,] 37.94 36.3200 33.9625
## [8,] 38.61 37.7250 36.3200
## [9,] 36.91 37.9125 37.7250
## [10,] 33.96 36.8550 37.9125
## [11,] 28.77 34.5625 36.8550
## [12,] 24.22 30.9650 34.5625
## [13,] 21.07 27.0050 30.9650
## [14,] 24.29 24.5875 27.0050
## [15,] 26.63 24.0525 24.5875
## [16,] 30.68 25.6675 24.0525
## [17,] 35.50 29.2750 25.6675
## [18,] 36.93 32.4350 29.2750
## [19,] 37.94 35.2625 32.4350
## [20,] 37.08 36.8625 35.2625
## [21,] 36.00 36.9875 36.8625
## [22,] 32.96 35.9950 36.9875
## [23,] 25.98 33.0050 35.9950
## [24,] 24.51 29.8625 33.0050
## [25,] 22.09 26.3850 29.8625
## [26,] 25.61 24.5475 26.3850
## [27,] 26.73 24.7350 24.5475
## [28,] 31.33 26.4400 24.7350
## [29,] 36.12 29.9475 26.4400
## [30,] 37.73 32.9775 29.9475
## [31,] 37.58 35.6900 32.9775
## [32,] 36.94 37.0925 35.6900
## [33,] 35.76 37.0025 37.0925
## [34,] 32.31 35.6475 37.0025
## [35,] 27.99 33.2500 35.6475
## [36,] 24.49 30.1375 33.2500
## [37,] 23.39 27.0450 30.1375
## [38,] 26.15 25.5050 27.0450
## [39,] 29.12 25.7875 25.5050
## [40,] 30.66 27.3300 25.7875
## [41,] 35.41 30.3350 27.3300
## [42,] 36.49 32.9200 30.3350
## [43,] 37.72 35.0700 32.9200
## [44,] 36.48 36.5250 35.0700
## [45,] 37.13 36.9550 36.5250
## [46,] 32.76 36.0225 36.9550
## [47,] 28.08 33.6125 36.0225
## [48,] 23.76 30.4325 33.6125
## [49,] 23.13 26.9325 30.4325
## [50,] 24.90 24.9675 26.9325
## [51,] 27.63 24.8550 24.9675
## [52,] 32.84 27.1250 24.8550
## [53,] 35.17 30.1350 27.1250
## [54,] 37.73 33.3425 30.1350
## [55,] 37.51 35.8125 33.3425
## [56,] 38.05 37.1150 35.8125
## [57,] 36.51 37.4500 37.1150
## [58,] 32.73 36.2000 37.4500
## [59,] 26.94 33.5575 36.2000
## [60,] 25.36 30.3850 33.5575
## [61,] 22.00 26.7575 30.3850
## [62,] 25.25 24.8875 26.7575
## [63,] 29.05 25.4150 24.8875
## [64,] 31.48 26.9450 25.4150
## [65,] 36.18 30.4900 26.9450
## [66,] 36.94 33.4125 30.4900
## [67,] 37.09 35.4225 33.4125
## [68,] 39.37 37.3950 35.4225
## [69,] 37.61 37.7525 37.3950
## [70,] 33.06 36.7825 37.7525
## [71,] 28.37 34.6025 36.7825
## [72,] 23.46 30.6250 34.6025
## [73,] 22.48 26.8425 30.6250
## [74,] 25.48 24.9475 26.8425
## [75,] 29.81 25.3075 24.9475
## [76,] 31.65 27.3550 25.3075
## [77,] 36.32 30.8150 27.3550
## [78,] 38.04 33.9550 30.8150
## [79,] 37.49 35.8750 33.9550
## [80,] 37.62 37.3675 35.8750
## [81,] 37.19 37.5850 37.3675
## [82,] 33.08 36.3450 37.5850
## [83,] 28.32 34.0525 36.3450
## [84,] 25.20 30.9475 34.0525
## [85,] 24.13 27.6825 30.9475
## [86,] 23.15 25.2000 27.6825
## [87,] 27.23 24.9275 25.2000
## [88,] 32.79 26.8250 24.9275
## [89,] 36.33 29.8750 26.8250
## [90,] 38.94 33.8225 29.8750
## [91,] 37.90 36.4900 33.8225
## [92,] 39.59 38.1900 36.4900
## [93,] 37.72 38.5375 38.1900
## [94,] 33.46 37.1675 38.5375
## [95,] 28.41 34.7950 37.1675
## [96,] 24.83 31.1050 34.7950
## [97,] 22.35 27.2625 31.1050
## [98,] 25.92 25.3775 27.2625
## [99,] 29.55 25.6625 25.3775
## [100,] 31.44 27.3150 25.6625
## [101,] 35.98 30.7225 27.3150
## [102,] 37.58 33.6375 30.7225
## [103,] 37.51 35.6275 33.6375
## [104,] 37.19 37.0650 35.6275
## [105,] 37.62 37.4750 37.0650
## [106,] 32.37 36.1725 37.4750
## [107,] 25.81 33.2475 36.1725
## [108,] 24.42 30.0550 33.2475
## [109,] 24.10 26.6750 30.0550
## [110,] 24.20 24.6325 26.6750
## [111,] 25.50 24.5550 24.6325
## [112,] 31.09 26.2225 24.5550
## [113,] 36.03 29.2050 26.2225
## [114,] 38.82 32.8600 29.2050
## [115,] 37.76 35.9250 32.8600
## [116,] NA NA 35.9250
## [117,] NA NA 35.9250
## [118,] NA NA 35.9250
## [119,] NA NA 35.9250
## [120,] NA NA 35.9250
## [121,] NA NA 35.9250
## [122,] NA NA 35.9250
## [123,] NA NA 35.9250
## [124,] NA NA 35.9250
## [125,] NA NA 35.9250
## [126,] NA NA 35.9250
## [127,] NA NA 35.9250
## [128,] NA NA 35.9250
## [129,] NA NA 35.9250
## [130,] NA NA 35.9250
## [131,] NA NA 35.9250
## [132,] NA NA 35.9250
## [133,] NA NA 35.9250
## [134,] NA NA 35.9250
## [135,] NA NA 35.9250
## [136,] NA NA 35.9250
## [137,] NA NA 35.9250
## [138,] NA NA 35.9250
## [139,] NA NA 35.9250
## [140,] NA NA 35.9250
## [141,] NA NA 35.9250
## [142,] NA NA 35.9250
## [143,] NA NA 35.9250
## [144,] NA NA 35.9250Plot forecasting
ts.plot(dataziq.ts, xlab="Time Period ", ylab="Data", main= "SMA m = 4")
points(dataziq.ts)
lines(data.gab[,2],col="green",lwd=2)
lines(data.gab[,3],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.5)
Akurasi data latih
#Menghitung nilai keakuratan data latih
error_train.sma = dt.latih.ts-data.ramal[1:length(dt.latih.ts)]
SSE_train.sma = sum(error_train.sma[5:length(dt.latih.ts)]^2)
MSE_train.sma = mean(error_train.sma[5:length(dt.latih.ts)]^2)
MAPE_train.sma = mean(abs((error_train.sma[5:length(dt.latih.ts)]/dt.latih.ts[5:length(dt.latih.ts)])*100))
akurasi_train.sma <- matrix(c(SSE_train.sma, MSE_train.sma, MAPE_train.sma))
row.names(akurasi_train.sma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi_train.sma) <- c("Akurasi m = 4")
akurasi_train.sma## Akurasi m = 4
## SSE 4407.87464
## MSE 39.71058
## MAPE 18.45565Dalam hal ini nilai MAPE data latih pada metode pemulusan SMA kurang dari 18%, nilai ini dapat dikategorikan sebagai nilai akurasi yang sangat baik.
Akurasi data uji
#Menghitung nilai keakuratan data uji
error_test.sma = dt.uji.ts-data.gab[116:144,3]
SSE_test.sma = sum(error_test.sma^2)
MSE_test.sma = mean(error_test.sma^2)
MAPE_test.sma = mean(abs((error_test.sma/dt.uji.ts*100)))
akurasi_test.sma <- matrix(c(SSE_test.sma, MSE_test.sma, MAPE_test.sma))
row.names(akurasi_test.sma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi_test.sma) <- c("Akurasi m = 4")
akurasi_test.sma## Akurasi m = 4
## SSE 1460.63012
## MSE 50.36656
## MAPE 20.75195Perhitungan akurasi menggunakan data latih menghasilkan nilai MAPE yang kurang dari 20% sehingga nilai akurasi ini dapat dikategorikan sebagai sangat baik.
Double Moving Average (DMA)
Pemulusan dan Peramalan
dma <- SMA(data.sma, n = 4)
At <- 2*data.sma - dma
Bt <- 2/(4-1)*(data.sma - dma)
data.dma<- At+Bt
# Peramalan
data.ramal2<- c(NA, data.dma)
t = 1:29
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}
data.gab2 <- cbind(aktual = c(dt.latih.ts,rep(NA,29)), pemulusan1 = c(data.sma,rep(NA,29)),pemulusan2 = c(data.dma, rep(NA,29)),At = c(At, rep(NA,29)), Bt = c(Bt,rep(NA,29)),ramalan = c(data.ramal2, f[-1]))
data.gab2## aktual pemulusan1 pemulusan2 At Bt ramalan
## [1,] 24.10 NA NA NA NA NA
## [2,] 25.12 NA NA NA NA NA
## [3,] 28.51 NA NA NA NA NA
## [4,] 32.99 27.6800 NA NA NA NA
## [5,] 36.16 30.6950 NA NA NA NA
## [6,] 38.19 33.9625 NA NA NA NA
## [7,] 37.94 36.3200 43.24604 40.47563 2.77041667 NA
## [8,] 38.61 37.7250 42.80729 40.77437 2.03291667 43.24604
## [9,] 36.91 37.9125 40.30000 39.34500 0.95500000 42.80729
## [10,] 33.96 36.8550 36.27479 36.50688 -0.23208333 40.30000
## [11,] 28.77 34.5625 30.89375 32.36125 -1.46750000 36.27479
## [12,] 24.22 30.9650 24.11708 26.85625 -2.73916667 30.89375
## [13,] 21.07 27.0050 18.10188 21.66313 -3.56125000 24.11708
## [14,] 24.29 24.5875 16.76667 19.89500 -3.12833333 18.10188
## [15,] 26.63 24.0525 19.71917 21.45250 -1.73333333 16.76667
## [16,] 30.68 25.6675 26.23312 26.00687 0.22625000 19.71917
## [17,] 35.50 29.2750 34.90729 32.65437 2.25291667 26.23312
## [18,] 36.93 32.4350 40.06417 37.01250 3.05166667 34.90729
## [19,] 37.94 35.2625 42.93333 39.86500 3.06833333 40.06417
## [20,] 37.08 36.8625 42.53542 40.26625 2.26916667 42.93333
## [21,] 36.00 36.9875 39.65521 38.58812 1.06708333 42.53542
## [22,] 32.96 35.9950 35.52521 35.71312 -0.18791667 39.65521
## [23,] 25.98 33.0050 28.49250 30.29750 -1.80500000 35.52521
## [24,] 24.51 29.8625 23.02917 25.76250 -2.73333333 28.49250
## [25,] 22.09 26.3850 18.17354 21.45812 -3.28458333 23.02917
## [26,] 25.61 24.5475 18.04333 20.64500 -2.60166667 18.17354
## [27,] 26.73 24.7350 21.98917 23.08750 -1.09833333 18.04333
## [28,] 31.33 26.4400 27.96188 27.35313 0.60875000 21.98917
## [29,] 36.12 29.9475 35.83083 33.47750 2.35333333 27.96188
## [30,] 37.73 32.9775 40.39833 37.43000 2.96833333 35.83083
## [31,] 37.58 35.6900 43.06708 40.11625 2.95083333 40.39833
## [32,] 36.94 37.0925 42.36854 40.25812 2.11041667 43.06708
## [33,] 35.76 37.0025 39.18896 38.31437 0.87458333 42.36854
## [34,] 32.31 35.6475 34.46312 34.93687 -0.47375000 39.18896
## [35,] 27.99 33.2500 29.08646 30.75187 -1.66541667 34.46312
## [36,] 24.49 30.1375 23.68437 26.26562 -2.58125000 29.08646
## [37,] 23.39 27.0450 19.58667 22.57000 -2.98333333 23.68437
## [38,] 26.15 25.5050 19.70604 22.02562 -2.31958333 19.58667
## [39,] 29.12 25.7875 23.56875 24.45625 -0.88750000 19.70604
## [40,] 30.66 27.3300 28.85187 28.24312 0.60875000 23.56875
## [41,] 35.41 30.3350 35.49437 33.43062 2.06375000 28.85187
## [42,] 36.49 32.9200 39.29812 36.74687 2.55125000 35.49437
## [43,] 37.72 35.0700 41.16375 38.72625 2.43750000 39.29812
## [44,] 36.48 36.5250 41.21250 39.33750 1.87500000 41.16375
## [45,] 37.13 36.9550 39.60083 38.54250 1.05833333 41.21250
## [46,] 32.76 36.0225 35.82146 35.90187 -0.08041667 39.60083
## [47,] 28.08 33.6125 30.00208 31.44625 -1.44416667 35.82146
## [48,] 23.76 30.4325 24.06062 26.60937 -2.54875000 30.00208
## [49,] 23.13 26.9325 18.90333 22.11500 -3.21166667 24.06062
## [50,] 24.90 24.9675 18.26958 20.94875 -2.67916667 18.90333
## [51,] 27.63 24.8550 21.61854 22.91312 -1.29458333 18.26958
## [52,] 32.84 27.1250 29.05000 28.28000 0.77000000 21.61854
## [53,] 35.17 30.1350 35.74229 33.49937 2.24291667 29.05000
## [54,] 37.73 33.3425 40.80604 37.82062 2.98541667 35.74229
## [55,] 37.51 35.8125 42.82708 40.02125 2.80583333 40.80604
## [56,] 38.05 37.1150 42.13792 40.12875 2.00916667 42.82708
## [57,] 36.51 37.4500 39.98333 38.97000 1.01333333 42.13792
## [58,] 32.73 36.2000 35.45937 35.75562 -0.29625000 39.98333
## [59,] 26.94 33.5575 29.35229 31.03437 -1.68208333 35.45937
## [60,] 25.36 30.3850 23.69646 26.37187 -2.67541667 29.35229
## [61,] 22.00 26.7575 18.47833 21.79000 -3.31166667 23.69646
## [62,] 25.25 24.8875 18.20521 20.87812 -2.67291667 18.47833
## [63,] 29.05 25.4150 23.00458 23.96875 -0.96416667 18.20521
## [64,] 31.48 26.9450 28.51792 27.88875 0.62916667 23.00458
## [65,] 36.18 30.4900 36.41604 34.04562 2.37041667 28.51792
## [66,] 36.94 33.4125 40.65729 37.75937 2.89791667 36.41604
## [67,] 37.09 35.4225 41.84750 39.27750 2.57000000 40.65729
## [68,] 39.37 37.3950 42.75333 40.61000 2.14333333 41.84750
## [69,] 37.61 37.7525 40.68062 39.50937 1.17125000 42.75333
## [70,] 33.06 36.7825 36.68979 36.72687 -0.03708333 40.68062
## [71,] 28.37 34.6025 31.21812 32.57187 -1.35375000 36.68979
## [72,] 23.46 30.6250 23.43229 26.30937 -2.87708333 31.21812
## [73,] 22.48 26.8425 17.89146 21.47187 -3.58041667 23.43229
## [74,] 25.48 24.9475 17.76937 20.64062 -2.87125000 17.89146
## [75,] 29.81 25.3075 22.60229 23.68437 -1.08208333 17.76937
## [76,] 31.65 27.3550 29.42479 28.59687 0.82791667 22.60229
## [77,] 36.32 30.8150 36.99625 34.52375 2.47250000 29.42479
## [78,] 38.04 33.9550 41.61646 38.55187 3.06458333 36.99625
## [79,] 37.49 35.8750 42.33333 39.75000 2.58333333 41.61646
## [80,] 37.62 37.3675 42.14146 40.23187 1.90958333 42.33333
## [81,] 37.19 37.5850 39.90062 38.97437 0.92625000 42.14146
## [82,] 33.08 36.3450 35.59812 35.89687 -0.29875000 39.90062
## [83,] 28.32 34.0525 30.24417 31.76750 -1.52333333 35.59812
## [84,] 25.20 30.9475 24.63917 27.16250 -2.52333333 30.24417
## [85,] 24.13 27.6825 20.05854 23.10812 -3.04958333 24.63917
## [86,] 23.15 25.2000 18.08229 20.92937 -2.84708333 20.05854
## [87,] 27.23 24.9275 21.15771 22.66562 -1.50791667 18.08229
## [88,] 32.79 26.8250 27.93542 27.49125 0.44416667 21.15771
## [89,] 36.33 29.8750 35.15521 33.04312 2.11208333 27.93542
## [90,] 38.94 33.8225 42.08917 38.78250 3.30666667 35.15521
## [91,] 37.90 36.4900 44.38479 41.22687 3.15791667 42.08917
## [92,] 39.59 38.1900 44.18271 41.78562 2.39708333 44.38479
## [93,] 37.72 38.5375 41.50000 40.31500 1.18500000 44.18271
## [94,] 33.46 37.1675 36.45292 36.73875 -0.28583333 41.50000
## [95,] 28.41 34.7950 30.83250 32.41750 -1.58500000 36.45292
## [96,] 24.83 31.1050 23.94458 26.80875 -2.86416667 30.83250
## [97,] 22.35 27.2625 18.39583 21.94250 -3.54666667 23.94458
## [98,] 25.92 25.3775 18.28167 21.12000 -2.83833333 18.39583
## [99,] 29.55 25.6625 22.84687 23.97312 -1.12625000 18.28167
## [100,] 31.44 27.3150 28.83271 28.22562 0.60708333 22.84687
## [101,] 35.98 30.7225 36.47771 34.17562 2.30208333 28.83271
## [102,] 37.58 33.6375 40.80937 37.94062 2.86875000 36.47771
## [103,] 37.51 35.6275 41.96396 39.42937 2.53458333 40.80937
## [104,] 37.19 37.0650 41.73479 39.86687 1.86791667 41.96396
## [105,] 37.62 37.4750 40.01458 38.99875 1.01583333 41.73479
## [106,] 32.37 36.1725 35.48500 35.76000 -0.27500000 40.01458
## [107,] 25.81 33.2475 28.67667 30.50500 -1.82833333 35.48500
## [108,] 24.42 30.0550 23.08417 25.87250 -2.78833333 28.67667
## [109,] 24.10 26.6750 18.57083 21.81250 -3.24166667 23.08417
## [110,] 24.20 24.6325 17.93250 20.61250 -2.68000000 18.57083
## [111,] 25.50 24.5550 21.34771 22.63062 -1.28291667 17.93250
## [112,] 31.09 26.2225 27.39125 26.92375 0.46750000 21.34771
## [113,] 36.03 29.2050 34.29042 32.25625 2.03416667 27.39125
## [114,] 38.82 32.8600 40.60896 37.50938 3.09958333 34.29042
## [115,] 37.76 35.9250 44.04479 40.79688 3.24791667 40.60896
## [116,] NA NA NA NA NA 44.04479
## [117,] NA NA NA NA NA 47.29271
## [118,] NA NA NA NA NA 50.54063
## [119,] NA NA NA NA NA 53.78854
## [120,] NA NA NA NA NA 57.03646
## [121,] NA NA NA NA NA 60.28438
## [122,] NA NA NA NA NA 63.53229
## [123,] NA NA NA NA NA 66.78021
## [124,] NA NA NA NA NA 70.02813
## [125,] NA NA NA NA NA 73.27604
## [126,] NA NA NA NA NA 76.52396
## [127,] NA NA NA NA NA 79.77188
## [128,] NA NA NA NA NA 83.01979
## [129,] NA NA NA NA NA 86.26771
## [130,] NA NA NA NA NA 89.51563
## [131,] NA NA NA NA NA 92.76354
## [132,] NA NA NA NA NA 96.01146
## [133,] NA NA NA NA NA 99.25938
## [134,] NA NA NA NA NA 102.50729
## [135,] NA NA NA NA NA 105.75521
## [136,] NA NA NA NA NA 109.00313
## [137,] NA NA NA NA NA 112.25104
## [138,] NA NA NA NA NA 115.49896
## [139,] NA NA NA NA NA 118.74688
## [140,] NA NA NA NA NA 121.99479
## [141,] NA NA NA NA NA 125.24271
## [142,] NA NA NA NA NA 128.49063
## [143,] NA NA NA NA NA 131.73854
## [144,] NA NA NA NA NA 134.98646Plot Forecasting
ts.plot(dataziq.ts, xlab="Time Period", ylab="Data", main= "DMA N=4")
points(dataziq.ts)
lines(data.gab2[,3],col="green",lwd=2)
lines(data.gab2[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)
Akurasi data latih
error_train.dma = dt.latih.ts-data.ramal2[1:length(dt.latih.ts)]
SSE_train.dma = sum(error_train.dma[8:length(dt.latih.ts)]^2)
MSE_train.dma = mean(error_train.dma[8:length(dt.latih.ts)]^2)
MAPE_train.dma = mean(abs((error_train.dma[8:length(dt.latih.ts)]/dt.latih.ts[8:length(dt.latih.ts)])*100))
akurasi_train.dma <- matrix(c(SSE_train.dma, MSE_train.dma, MAPE_train.dma))
row.names(akurasi_train.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi_train.dma) <- c("Akurasi m = 4")
akurasi_train.dma## Akurasi m = 4
## SSE 4620.34145
## MSE 42.78094
## MAPE 19.28705Perhitungan akurasi pada data latih menggunakan nilai MAPE menghasilkan nilai MAPE yang kurang dari 10% sehingga dikategorikan sangat baik.
Akurasi data uji
error_test.dma = dt.uji.ts-data.gab2[116:144,6]
SSE_test.dma = sum(error_test.dma^2)
MSE_test.dma = mean(error_test.dma^2)
MAPE_test.dma = mean(abs((error_test.dma/dt.uji.ts*100)))
akurasi_test.dma <- matrix(c(SSE_test.dma, MSE_test.dma, MAPE_test.dma))
row.names(akurasi_test.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi_test.dma) <- c("Akurasi m = 4")
akurasi_test.dma## Akurasi m = 4
## SSE 119036.1225
## MSE 4104.6939
## MAPE 192.0865Perhitungan akurasi menggunakan data latih menghasilkan nilai MAPE yang kurang dari 10%
Winter aditif
Pemulusan dan Peramalan
# optimum
winter.opt<- HoltWinters(dt.latih.ts, alpha= NULL, beta = NULL, gamma = NULL, seasonal = "additive")
winter.opt## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = dt.latih.ts, alpha = NULL, beta = NULL, gamma = NULL, seasonal = "additive")
##
## Smoothing parameters:
## alpha: 0.1347328
## beta : 0.07756776
## gamma: 0.2846425
##
## Coefficients:
## [,1]
## a 31.11710557
## b -0.02136399
## s1 6.51612115
## s2 5.70991374
## s3 1.42323689
## s4 -3.92052769
## s5 -6.82342932
## s6 -8.37750095
## s7 -6.64513911
## s8 -3.69226249
## s9 0.21277882
## s10 4.66285490
## s11 6.66189996
## s12 6.18527997winter.opt$fitted## xhat level trend season
## Jan 2 21.75672 32.02042 -0.1025480769 -10.1611458
## Feb 2 24.83822 31.82534 -0.1097249626 -6.8773958
## Mar 2 27.09057 31.64176 -0.1154544056 -4.4357292
## Apr 2 31.03783 31.46425 -0.1202678102 -0.3061458
## May 2 35.84353 31.29577 -0.1240075009 4.6717708
## Jun 2 37.20381 31.12548 -0.1275977147 6.2059375
## Jul 2 36.83480 30.96099 -0.1304593328 6.0042708
## Aug 2 37.69563 30.97943 -0.1189089435 6.8351042
## Sep 2 35.90026 30.77758 -0.1253428291 5.2480208
## Oct 2 33.01398 30.66567 -0.1243004231 2.4726042
## Nov 2 27.81559 30.53410 -0.1248645490 -2.5936458
## Dec 2 22.95423 30.16192 -0.1440481702 -7.0636458
## Jan 3 19.76942 30.22749 -0.1277889274 -10.3302800
## Feb 3 23.29640 30.41236 -0.1035367186 -7.0124188
## Mar 3 25.99202 30.62054 -0.0793574868 -4.5491643
## Apr 3 30.17469 30.64061 -0.0716448793 -0.3942773
## May 3 35.25222 30.72462 -0.0595708171 4.5871620
## Jun 3 36.86997 30.78197 -0.0505016746 6.1384991
## Jul 3 37.08231 30.84735 -0.0415135730 6.2764732
## Aug 3 37.52006 30.87289 -0.0363122030 6.6834799
## Sep 3 35.98864 30.75842 -0.0423743171 5.2725867
## Oct 3 33.09979 30.68524 -0.0447637680 2.4593097
## Nov 3 27.43531 30.53407 -0.0530178008 -3.0457369
## Dec 3 23.82809 30.55579 -0.0472208308 -6.6804722
## Jan 4 20.79870 30.59775 -0.0403032937 -9.7587400
## Feb 4 24.45075 30.90658 -0.0132218472 -6.4425986
## Mar 4 26.75943 31.12230 0.0045368308 -4.3674051
## Apr 4 31.36435 31.44488 0.0292069699 -0.1097338
## May 4 36.20192 31.37919 0.0218458225 4.8008901
## Jun 4 37.65822 31.29434 0.0135694870 6.3503174
## Jul 4 37.55092 31.15051 0.0013604841 6.3990514
## Aug 4 37.71839 31.17465 0.0031275341 6.5406170
## Sep 4 36.21738 31.01092 -0.0098148153 5.2162757
## Oct 4 33.38858 31.12407 -0.0002771453 2.2647909
## Nov 4 28.12313 31.03910 -0.0068464073 -2.9091225
## Dec 4 24.50170 31.02644 -0.0072971723 -6.5174500
## Jan 5 21.78364 30.91921 -0.0150485766 -9.1205248
## Feb 5 25.06050 31.08557 -0.0009778908 -6.0240885
## Mar 5 27.27429 31.06296 -0.0026552491 -3.7860159
## Apr 5 30.82608 31.10823 0.0010622321 -0.2832103
## May 5 36.00859 31.38064 0.0221095004 4.6058457
## Jun 5 37.36570 31.28976 0.0133454524 6.0625937
## Jul 5 37.81004 31.35219 0.0171527368 6.4406946
## Aug 5 37.57854 31.32892 0.0140170881 6.2356110
## Sep 5 36.86644 31.40645 0.0189442329 5.4410453
## Oct 5 33.50257 31.37737 0.0152190693 2.1099763
## Nov 5 28.37590 31.28850 0.0071450114 -2.9197455
## Dec 5 24.39420 31.10218 -0.0078614859 -6.7001235
## Jan 6 22.43775 31.22445 0.0022320225 -8.7889278
## Feb 6 25.10174 31.16770 -0.0023428898 -6.0636180
## Mar 6 27.48613 31.18533 -0.0007934292 -3.6984078
## Apr 6 31.62359 31.39524 0.0155504363 0.2128004
## May 6 35.80480 31.39145 0.0140497396 4.3993077
## Jun 6 37.62634 31.45605 0.0179708753 6.1523181
## Jul 6 37.75914 31.38155 0.0107980185 6.3667983
## Aug 6 37.65772 31.30219 0.0038048546 6.3517266
## Sep 6 36.91165 31.53669 0.0216997493 5.3532562
## Oct 6 33.60118 31.65248 0.0289981463 1.9196989
## Nov 6 28.35851 31.60857 0.0233422969 -3.2733963
## Dec 6 25.19467 31.63346 0.0234623362 -6.4622547
## Jan 7 22.53179 31.42320 0.0053334836 -8.8967425
## Feb 7 25.39925 31.42156 0.0047921799 -6.0271026
## Mar 7 28.12963 31.43723 0.0056361086 -3.3132399
## Apr 7 31.86990 31.66927 0.0231975541 0.1774342
## May 7 36.17545 31.66284 0.0208993951 4.4917152
## Jun 7 37.70890 31.70321 0.0224100472 5.9832789
## Jul 7 37.99810 31.77023 0.0258703285 6.2019938
## Aug 7 38.52165 31.72765 0.0205602488 6.7734469
## Sep 7 37.16311 31.62672 0.0111371459 5.5252540
## Oct 7 33.43931 31.64148 0.0114181207 1.7864103
## Nov 7 28.34159 31.60449 0.0076629834 -3.2705674
## Dec 7 24.72719 31.60924 0.0074373914 -6.8894886
## Jan 8 22.78326 31.68038 0.0123786526 -8.9094991
## Feb 8 25.89345 31.87421 0.0264532814 -6.0072142
## Mar 8 28.62944 31.53103 -0.0022183208 -2.8993779
## Apr 8 31.44670 31.34026 -0.0168437273 0.1232746
## May 8 36.02892 31.50441 -0.0028049595 4.5273160
## Jun 8 37.60734 31.54217 0.0003416116 6.0648257
## Jul 8 37.81319 31.72206 0.0142691780 6.0768537
## Aug 8 38.31458 31.74803 0.0151764556 6.5513772
## Sep 8 37.49543 31.93505 0.0285057352 5.5318756
## Oct 8 33.72258 31.99381 0.0308527207 1.6979148
## Nov 8 28.74151 31.98928 0.0281085482 -3.2758838
## Dec 8 25.22433 31.97273 0.0246439707 -6.7730403
## Jan 9 23.38696 31.94424 0.0205228419 -8.5778092
## Feb 9 25.15184 31.82505 0.0096856948 -6.6829039
## Mar 9 28.71190 31.93824 0.0177137255 -3.2440477
## Apr 9 32.54946 32.06887 0.0264726219 0.4541194
## May 9 36.56221 31.94586 0.0148777314 4.6014697
## Jun 9 38.28414 31.88230 0.0087931185 6.3930499
## Jul 9 37.89589 31.79622 0.0014342169 6.0982351
## Aug 9 38.60856 31.74566 -0.0025986668 6.8655018
## Sep 9 37.12170 31.55194 -0.0174239697 5.5871858
## Oct 9 33.22268 31.60165 -0.0122162446 1.6332442
## Nov 9 28.09589 31.47455 -0.0211275046 -3.3575318
## Dec 9 24.23026 31.14544 -0.0450171604 -6.8701609
## Jan 10 22.24975 31.12598 -0.0430341948 -8.8332030
## Feb 10 24.81483 31.33224 -0.0236973489 -6.4937112
## Mar 10 28.15795 31.22570 -0.0301229045 -3.0376311
## Apr 10 30.96044 30.83747 -0.0579009389 0.1808683
## May 10 35.19855 30.79702 -0.0565468752 4.4580766
## Jun 10 37.02427 30.85250 -0.0478574973 6.2196262
## Jul 10 37.02069 31.04659 -0.0290904547 6.0031941xhat.opt <- winter.opt$fitted[,2]
# Peramalan
forecast.opt <- predict(winter.opt, n.ahead = 29)
forecast.opt## Jan Feb Mar Apr May Jun Jul Aug
## 10 37.61186
## 11 22.61142 24.32242 27.25393 31.13761 35.56632 37.54400 37.04602 37.35549
## 12 22.35505 24.06605 26.99756 30.88124 35.30995 37.28763 36.78965 37.09913
## Sep Oct Nov Dec
## 10 36.78429 32.47625 27.11112 24.18686
## 11 36.52792 32.21988 26.85475 23.93049
## 12 36.27156 31.96351 26.59839 23.67412Plot forecasting
plot(dt.latih.ts,main="Winter optimum",type="l",col="black",
xlim=c(1,25),pch=20)
lines(xhat.opt,type="l",col="green")
lines(forecast.opt,type="l",col="green")
Akurasi data latih
SSE.opt<-winter.opt$SSE
MSE.opt<-winter.opt$SSE/length(dt.latih.ts)
RMSE.opt<-sqrt(MSE.opt)
akurasi.opt <- matrix(c(SSE.opt,MSE.opt,RMSE.opt))
row.names(akurasi.opt)<- c("SSE.opt", "MSE.opt", "RMSE.opt")
colnames(akurasi.opt) <- c("Akurasi")
akurasi.opt## Akurasi
## SSE.opt 116.780067
## MSE.opt 1.015479
## RMSE.opt 1.007710Akurasi data uji
forecast.opt<-data.frame(forecast.opt)
dt.uji.ts <- data.frame(dt.uji.ts)
selisih.opt<-(forecast.opt)-(dt.uji.ts)
SSEuji.opt<-sum(selisih.opt^2)
MSEuji.opt<-SSEuji.opt/length(dt.uji.ts)
akurasi.opt <- matrix(c(SSE.opt,MSE.opt))
row.names(akurasi.opt)<- c("SSEuji.opt", "MSEuji.opt")
colnames(akurasi.opt) <- c("Akurasi")
akurasi.opt## Akurasi
## SSEuji.opt 116.780067
## MSEuji.opt 1.015479Penanganan Autokorelasi
Input data
datasol <- read.csv("https://raw.githubusercontent.com/Raziqizzan03/MPDW/main/data/Pertemuan1/data%20tugas%203.csv")
View(datasol)
str(datasol)## 'data.frame': 200 obs. of 10 variables:
## $ open : num 11.1 10.8 10.9 11.4 11.4 ...
## $ High : num 11.6 11.3 11 11.4 11.7 ...
## $ low : num 11.1 10.8 10.7 10.9 11.3 ...
## $ close : num 11.1 10.8 10.9 11.4 11.4 ...
## $ volume : int 4864322 5520210 7038886 7282602 7075356 5626819 5579050 9052957 4987451 5274492 ...
## $ close.t.1. : num NA 11.1 10.8 10.9 11.4 ...
## $ open.t.1. : num NA 11.1 10.8 10.9 11.4 ...
## $ high.t.1. : num NA 11.6 11.3 11 11.4 ...
## $ low.t.1. : num NA 11.1 10.8 10.7 10.9 ...
## $ volume.t.1.: int NA 4864322 5520210 7038886 7282602 7075356 5626819 5579050 9052957 4987451 ...# Peubah yang digunakan
high <- datasol$High
low <- datasol$low
close <- datasol$close
Volume <- datasol$volume
Open <- datasol$openMengubah data menjadi data deret waktu
high.ts <- ts(high)
low.ts <- ts(low)
close.ts <- ts(close)
open.ts <- ts(Open)
volume.ts <- ts(Volume)Pembuatan regresi time series
Pembuatan model awal
model1 <- lm(Volume~high+low)
summary(model1)##
## Call:
## lm(formula = Volume ~ high + low)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13894618 -5097724 -2745142 541551 51971867
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36129309 13400383 -2.696 0.00763 **
## high 16558970 4992203 3.317 0.00108 **
## low -13398593 4531796 -2.957 0.00349 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9333000 on 196 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.06787, Adjusted R-squared: 0.05836
## F-statistic: 7.135 on 2 and 196 DF, p-value: 0.00102Model yang dihasilkan adalah \[y_i=-36129324+ 16558983x1_t+-13398606x2_t\] Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model. Hasil uji-t parsial parameter regresi, yaitu koefisien regresi juga menunjukkan hal yang sama, yaitu memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%. Selanjutnya dapat dilihat juga nilai \(R^2=0.06787\). Artinya, sebesar 6% keragaman nilai volume dapat dijelaskan oleh peubah high dan low. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.
plot asumsi
#sisaan dan fitted value
sisaan1<- residuals(model1)
fitValue1<- predict(model1)
#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan1)
qqline(sisaan1, col = "steelblue", lwd = 2)
plot(fitValue1, sisaan1, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan1, col = "steelblue")
plot(seq(1,199,1), sisaan1, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,199,1), sisaan1, col = "red")
abline(a = 0, b = 0, lwd = 2)
Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot di atas menunjukkan bahwa sisaan cenderung tidak menyebar normal, dan histogram dari sisaan menunjukkan demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. Plot Sisaan vs Fitted Value dan Plot Sisaan vs Order menunjukkan adanya pola pada sisaan. Untuk lebih lanjut akan digunakan uji formal melihat normalitas sisaan dan plot ACF dan PACF untuk melihat apakah ada autokorelasi atau tidak.
Uji formal Normalitas
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
ks.test(sisaan1, "pnorm", mean=mean(sisaan1), sd=sd(sisaan1))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan1
## D = 0.2401, p-value = 2.172e-10
## alternative hypothesis: two-sidedBerdasarkan uji formal Kolmogorov-Smirnov didapatkan nilai p-value < \(\alpha\) (5%). Artinya, belum cukup bukti untuk menyatakan sisaan berdistribusi normal.
Transformasi karna data tidak normal (1/x dan 1/y)
Volumetrans <- 1/(Volume)
hightrans <- 1/(high)
lowtrans <- 1/(low)
# Model transformasi
modeltrans <- lm(Volumetrans~hightrans+lowtrans)
summary(modeltrans)##
## Call:
## lm(formula = Volumetrans ~ hightrans + lowtrans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.210e-07 -4.025e-08 -2.120e-09 3.496e-08 3.756e-07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.296e-07 8.848e-08 -2.595 0.0102 *
## hightrans 2.662e-05 6.023e-06 4.420 1.63e-05 ***
## lowtrans -2.118e-05 5.199e-06 -4.074 6.70e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.522e-08 on 196 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1021, Adjusted R-squared: 0.0929
## F-statistic: 11.14 on 2 and 196 DF, p-value: 2.618e-05Model transforasi yang dihasilkan adalah \[y_i=-2.296e-07+ 2.662e-05x1_t+-2.118e-05x2_t\] Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model. Hasil uji-t parsial parameter regresi, yaitu koefisien regresi juga menunjukkan hal yang sama, yaitu memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%. Selanjutnya dapat dilihat juga nilai \(R^2=0.1021\). Artinya, sebesar 10% keragaman nilai volume dapat dijelaskan oleh peubah high dan low. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.
cek ulang normalitas
#sisaan dan fitted value
sisaan1<- residuals(modeltrans)
fitValue1<- predict(modeltrans)
#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan1)
qqline(sisaan1, col = "steelblue", lwd = 2)
plot(fitValue1, sisaan1, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan1, col = "steelblue")
plot(seq(1,199,1), sisaan1, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,199,1), sisaan1, col = "red")
abline(a = 0, b = 0, lwd = 2)
# Uji Formal normalitas
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
ks.test(sisaan1, "pnorm", mean=mean(sisaan1), sd=sd(sisaan1))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan1
## D = 0.063155, p-value = 0.4054
## alternative hypothesis: two-sidedDua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot di atas menunjukkan bahwa sisaan cenderung menyebar normal, dan histogram dari sisaan menunjukkan demikian. Berdasarkan uji formal Kolmogorov-Smirnov didapatkan nilai p-value > \(\alpha\) (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.
ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan1)
pacf(sisaan1)
Berdasarkan plot ACF dan PACF, terlihat terdapat data yang berada diluar rentang batas sehingga ada yang signifikan. Namun, untuk lebih memastikan akan dilakukan uji formal dengan uji Durbin Watson.
Uji formal autokorelasi
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(modeltrans)##
## Durbin-Watson test
##
## data: modeltrans
## DW = 1.6646, p-value = 0.006793
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 1.6646\) dan p-value = \(0.006793\). Dengan nilai p-value < 0.05 dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penangan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.
Metode Cochrane-Orcutt
#Penanganan Autokorelasi Cochrane-Orcutt
modelCO1<-cochrane.orcutt(modeltrans)
modelCO1## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = Volumetrans ~ hightrans + lowtrans)
##
## number of interaction: 6
## rho 0.173819
##
## Durbin-Watson statistic
## (original): 1.66458 , p-value: 6.793e-03
## (transformed): 2.03818 , p-value: 5.794e-01
##
## coefficients:
## (Intercept) hightrans lowtrans
## 0.0e+00 2.3e-05 -1.9e-05modelCO1$coefficients## (Intercept) hightrans lowtrans
## -1.916405e-07 2.339588e-05 -1.854018e-05Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-1.916405e-07+2.339589e-05x1_t+-1.854019e-05x2_t\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.03818\) dan \(0.579\). dengan nilai p-value > 0.05, artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.173819\). Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho1<- modelCO1$rho
rho1## [1] 0.1738193Transformasi manual
#Transformasi Manual
volume.trans<- Volumetrans[-1]-Volumetrans[-199]*rho1
high.trans<- hightrans[-1]-hightrans[-199]*rho1
low.trans<- lowtrans[-1]-lowtrans[-199]*rho1
modelCOmanual1<- lm(volume.trans~high.trans+low.trans)
summary(modelCOmanual1)##
## Call:
## lm(formula = volume.trans ~ high.trans + low.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.244e-07 -4.176e-08 -6.260e-09 3.392e-08 3.863e-07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.583e-07 8.394e-08 -1.886 0.060760 .
## high.trans 2.340e-05 6.107e-06 3.831 0.000172 ***
## low.trans -1.854e-05 5.271e-06 -3.517 0.000542 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.439e-08 on 195 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.07585, Adjusted R-squared: 0.06637
## F-statistic: 8.002 on 2 and 195 DF, p-value: 0.0004571#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang1 <- modelCOmanual1$coefficients[1]
b0.1<- b0bintang1/(1-rho1)
b1.1 <- modelCOmanual1$coefficients[2]
b2.1 <- modelCOmanual1$coefficients[3]
b0.1## (Intercept)
## -1.916405e-07b1.1## high.trans
## 2.339588e-05b2.1## low.trans
## -1.854018e-05Metode Hildreth-Lu
hildreth.lu.func<- function(r, model){
x1 <- model.matrix(model)[,c(-1,-3)]
x2 <- model.matrix(model)[,c(-1,-2)]
y <- model.response(model.frame(model))
n <- length(y)
t <- 2:n
y <- y[t]-r*y[t-1]
x1 <- x1[t]-r*x1[t-1]
x2 <- x2[t]-r*x2[t-1]
return(lm(y~x1+x2))
}
#Pencariab rho yang meminimumkan SSE
r1 <- c(seq(0.1,0.9, by= 0.1))
tab1 <- data.frame("rho" = r1, "SSE" = sapply(r1, function(i){deviance(hildreth.lu.func(i, modeltrans))}))
tab1Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.2. Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali \(ρ\) yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar \(ρ\) yang dicari adalah 0.1, kali ini jarak antar \(ρ\) adalah 0.001 dan dilakukan pada selang 0.1 sampai dengan 0.4.
rOpt1<- seq(0.1,0.4, by= 0.001)
tabOpt1 <- data.frame("rho" = rOpt1, "SSE" = sapply(rOpt1, function(i){deviance(hildreth.lu.func(i, modeltrans))}))
head(tabOpt1[order(tabOpt1$SSE),])#Grafik SSE optimum
par(mfrow = c(1,1))
plot(tab1$SSE ~ tab1$rho , type = "l", xlab = "Rho", ylab = "SSE")
abline(v = tabOpt1[tabOpt1$SSE==min(tabOpt1$SSE),"rho"], lty = 2, col="red",lwd=2)
text(x=0.174, y=8.084217e-13 , labels = "rho=0.174", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(ρ=0.174\). Hal tersebut juga ditunjukkan
pada plot. Selanjutnya, model dapat didapatkan dengan mengevaluasi nilai
\(ρ\) ke dalam fungsi
hildreth.lu.func, serta dilanjutkan dengan pengujian
autokorelasi dengan uji Durbin-Watson. Namun, setelah pengecekan
tersebut tidak lupa koefisien regresi tersebut digunakan untuk
transformasi balik. Persamaan hasil transformasi itulah yang menjadi
persamaan sesungguhnya.
#Model terbaik
modelHL1 <- hildreth.lu.func(0.174, modeltrans)
summary(modelHL1)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.244e-07 -4.176e-08 -6.260e-09 3.393e-08 3.863e-07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.583e-07 8.394e-08 -1.886 0.060843 .
## x1 2.339e-05 6.107e-06 3.830 0.000173 ***
## x2 -1.854e-05 5.271e-06 -3.517 0.000543 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.439e-08 on 195 degrees of freedom
## Multiple R-squared: 0.07582, Adjusted R-squared: 0.06635
## F-statistic: 7.999 on 2 and 195 DF, p-value: 0.0004582#Transformasi Balik
cat("y = ", coef(modelHL1)[1]/(1-0.174), "+", coef(modelHL1)[2],"x1","+", coef(modelHL1)[3],"x2", sep = "")## y = -1.916084e-07+2.339269e-05x1+-1.85375e-05x2Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-1.916083e-07+ 2.339e-05x1_t+-1.854e-05x2_t\]
#Deteksi autokorelasi
dwtest(modelHL1)##
## Durbin-Watson test
##
## data: modelHL1
## DW = 2.0386, p-value = 0.5805
## alternative hypothesis: true autocorrelation is greater than 0dwtest(modelCOmanual1)##
## Durbin-Watson test
##
## data: modelCOmanual1
## DW = 2.0382, p-value = 0.5794
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bahwa p-value sebesar \(0.5805\), di mana p-value > \(\alpha=5%\). Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai volume dengan metode Hildreth-Lu pada taraf nyata 5%.
Perbandingan keakuratan
sseModelawal1 <- anova(modeltrans)$`Sum Sq`[c(-1,-2)]
sseModelCO1 <- anova(modelCOmanual1)$`Sum Sq`[c(-1,-2)]
sseModelHL1 <- anova(modelHL1)$`Sum Sq`[c(-1,-2)]
mseModelawal1 <- sseModelawal1/length(Volumetrans)
mseModelCO1 <- sseModelCO1/length(Volumetrans)
mseModelHL1 <- sseModelHL1/length(Volumetrans)
akurasi1 <- matrix(c(sseModelawal1,sseModelCO1,sseModelHL1,
mseModelawal1,mseModelCO1,mseModelHL1),nrow=2,ncol=3,byrow = T)
colnames(akurasi1) <- c("Model Awal", "Model Cochrane-Orcutt", "Model Hildreth-Lu")
row.names(akurasi1) <- c("SSE","MSE")
akurasi1## Model Awal Model Cochrane-Orcutt Model Hildreth-Lu
## SSE 8.33816e-13 8.084217e-13 8.084217e-13
## MSE 4.16908e-15 4.042108e-15 4.042109e-15Simpulan
Autokorelasi yang terdapat pada data Volume terjadi akibat adanya korelasi di antara unsur penyusunnya. Adanya autokorelasi menyebabkan model regresi kurang baik karena akan meingkatkan galatnya. Autokorelasi dapat dideteksi secara eksploratif melalui plot sisaan, ACF, dan PACF, serta dengan uji formal Durbin-Watson. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang sama, artinya keduanya baik untuk digunakan
CEK Kestasioneran data
Input data
dataSmekah<- read.csv("https://raw.githubusercontent.com/Raziqizzan03/MPDW/main/data/Pertemuan1/data%20suhu%20mekkah.csv")Ubah data ke ts
dataSmekah.ts <- ts(dataSmekah$suhu.mekkah)Membuat plot time series
ts.plot(dataSmekah.ts, xlab="Time Period ", ylab="Suhu Mekkah",
main = "Time Series Plot")
points(dataSmekah.ts)
dalam plot terlihat bahwa data stasioner dalam rataan dan ragam
Plot ACF
acf(dataSmekah.ts)
Berdasarkan plot ACF, terlihat bahwa plot ACF pada data tersebut cenderung tails off dan membentuk gelombang sinus.
Uji ADF
tseries::adf.test(dataSmekah.ts)## Warning in tseries::adf.test(dataSmekah.ts): p-value smaller than printed
## p-value##
## Augmented Dickey-Fuller Test
##
## data: dataSmekah.ts
## Dickey-Fuller = -11.127, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary\(H_0\) : Data tidak stasioner dalam rataan
\(H_1\) : Data stasioner dalam rataan
Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF.
Plot Box-Cox
index <- seq(1:144)
bc = boxcox(dataSmekah.ts~index, lambda = seq(0,4,by=0.01))
#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 1.68#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]## [1] 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.70 0.71
## [16] 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86
## [31] 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00 1.01
## [46] 1.02 1.03 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
## [61] 1.17 1.18 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
## [76] 1.32 1.33 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
## [91] 1.47 1.48 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
## [106] 1.62 1.63 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
## [121] 1.77 1.78 1.79 1.80 1.81 1.82 1.83 1.84 1.85 1.86 1.87 1.88 1.89 1.90 1.91
## [136] 1.92 1.93 1.94 1.95 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05 2.06
## [151] 2.07 2.08 2.09 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21
## [166] 2.22 2.23 2.24 2.25 2.26 2.27 2.28 2.29 2.30 2.31 2.32 2.33 2.34 2.35 2.36
## [181] 2.37 2.38 2.39 2.40 2.41 2.42 2.43 2.44 2.45 2.46 2.47 2.48 2.49 2.50 2.51
## [196] 2.52 2.53 2.54 2.55 2.56 2.57 2.58 2.59 2.60 2.61 2.62 2.63 2.64 2.65 2.66
## [211] 2.67 2.68 2.69 2.70 2.71 2.72 2.73 2.74 2.75 2.76 2.77 2.78 2.79 2.80 2.81
## [226] 2.82 2.83Gambar di atas menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 1.68 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0.57 dan batas atas 2.83. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.
Kesimpulan: Maka data tersebut stasioner dalam rataan dan ragam
Analisis seluruhnya
Pembangkitan Data
Data yang akan dibangkitkan adalah data dengan model MA(2) sebagai berikut.
set.seed(99)
ma2 <- arima.sim(list(order = c(0,0,2), ma = c(0.55,0.65)), n = 175)Data kemudian dibagi menjadi data latih dan data uji. Pembagian kali ini dilakukan dengan proporsi / perbandingan, yaitu 80:20.
ma2 <- ma2[-c(1:25)]
ma2.train <- ma2[1:120]
ma2.test <- ma2[121:150]Eksplorasi Data
Sebelum masuk dalam tahap pemodelan, dilakukan eksplorasi data dengan plot deret waktu untuk melihat pola data.
#--PLOT TIME SERIES--#
plot(ma2.train,
col = "navyblue",
lwd = 1,
type = "o",
xlab = "Time",
ylab = "Data")
Berdasarkan plot data deret waktu di atas, terlihat data cenderung stasioner dalam rataan dan ragam. Data stasioner dalam rataan karena menyebar/bergerak di sekitar nilai tengahnya (0) dan dikatakan stasioner dalam ragam karena memiliki lebar pita yang cenderung sama. Selain dengan plot data deret waktu, akan dilakukan pengecekan stasioneritas data dengan plot ACF dan uji ADF.
#--CEK KESTASIONERAN---#
acf(ma2.train, main="ACF", lag.max=20)
Berdasarkan plot ACF di atas, dapat dilihat bahwa plot cuts off pada lag ke-2. Hal ini sesuai dengan proses pembangkitan model MA(2).
adf.test(ma2.train) ## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -4.91 0.01
## [2,] 1 -4.89 0.01
## [3,] 2 -5.65 0.01
## [4,] 3 -3.57 0.01
## [5,] 4 -3.43 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -4.99 0.01
## [2,] 1 -4.95 0.01
## [3,] 2 -5.78 0.01
## [4,] 3 -3.69 0.01
## [5,] 4 -3.59 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -5.10 0.0100
## [2,] 1 -5.04 0.0100
## [3,] 2 -6.02 0.0100
## [4,] 3 -3.99 0.0123
## [5,] 4 -3.96 0.0138
## ----
## 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.01358 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.
Spesifikasi Model
#---SPESIFIKASI MODEL---#
par(mfrow = c(1,2))
acf(ma2.train, main="ACF", lag.max=20) #ARIMA(0,0,2)
pacf(ma2.train, main="PACF", lag.max=20) #ARIMA(1,0,0)
par(mfrow = c(1,1))Berdasarkan Plot ACF, terlihat cuts off pada lag ke-2 sehingga dapat kita asumsikan model yang terbentuk adalah ARIMA(0,0,2). Selanjutnya, berdasarkan plot PACF, terlihat cuts off pada lag pertama sehingga model yang terbentuk adalah ARIMA(1,0,0). Selain dengan plot ACF dan PACF, penentuan spesifikasi model dilakukan dengan extended ACF (EACF) berikut ini.
eacf(ma2.train) ## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 o 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 o o o o o o o o o o o
## 4 o o o x o o o o o o o o o o
## 5 x o o x o o o o o o o o o o
## 6 x o o o o o o o x o o o o o
## 7 o o o o o x o o o o o o o o#ARIMA(0,0,2) #ARIMA(1,0,3) #ARIMA(2,0,3) #ARIMA(3,0,3)
#Terdapat 5 model tentatifMenggunakan plot EACF, dapat diambil beberapa model dengan melihat ujung segitiga yang terbentuk, antara lain ARIMA(0,0,2), ARIMA(1,0,3), ARIMA(2,0,3), dan ARIMA(3,0,3).
Pendugaan Parameter
Selanjutnya akan dilakukan pendugaan parameter kelima model ARIMA
yang terbentuk sebelumnya. Pendugaan dilakukan dengan fungsi
Arima() yang dilanjutkan dengan melihat nilai AIC pada
ringkasan data dan melihat signifikansi parameter.
#---PENDUGAAN PARAMETER MODEL---#
model1.ma2=Arima(ma2.train, order=c(0,0,2),method="ML")
summary(model1.ma2) #AIC=326.87## Series: ma2.train
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## 0.6953 0.5349 -0.2082
## s.e. 0.0829 0.0673 0.1841
##
## sigma^2 = 0.8497: log likelihood = -159.44
## AIC=326.87 AICc=327.22 BIC=338.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004254881 0.910219 0.73195 -767.1756 938.3389 0.852277
## ACF1
## Training set 0.05965759lmtest::coeftest(model1.ma2) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.695312 0.082903 8.3871 < 2.2e-16 ***
## ma2 0.534929 0.067337 7.9441 1.957e-15 ***
## intercept -0.208236 0.184115 -1.1310 0.258
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model2.ma2=Arima(ma2.train, order=c(1,0,0),method="ML")
summary(model2.ma2) #AIC=340.47## Series: ma2.train
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.6428 -0.2185
## s.e. 0.0690 0.2450
##
## sigma^2 = 0.9625: log likelihood = -167.24
## AIC=340.47 AICc=340.68 BIC=348.84
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.002975587 0.9728552 0.7723112 -1335.017 1507.476 0.8992733
## ACF1
## Training set 0.06622426lmtest::coeftest(model2.ma2) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.642823 0.069023 9.3131 <2e-16 ***
## intercept -0.218485 0.245000 -0.8918 0.3725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model3.ma2=Arima(ma2.train, order=c(1,0,3),method="ML")
summary(model3.ma2) #AIC=329.22## Series: ma2.train
## ARIMA(1,0,3) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 ma3 mean
## 0.0610 0.684 0.6137 0.0934 -0.2082
## s.e. 0.4257 0.416 0.3263 0.2714 0.2083
##
## sigma^2 = 0.8518: log likelihood = -158.61
## AIC=329.22 AICc=329.96 BIC=345.94
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004465851 0.9034798 0.725319 -548.3662 713.5294 0.844556
## ACF1
## Training set 0.005156424lmtest::coeftest(model3.ma2) #tidak ada yang signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.060970 0.425686 0.1432 0.88611
## ma1 0.684021 0.415962 1.6444 0.10009
## ma2 0.613727 0.326319 1.8808 0.06001 .
## ma3 0.093411 0.271412 0.3442 0.73072
## intercept -0.208191 0.208338 -0.9993 0.31765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model4.ma2=Arima(ma2.train, order=c(2,0,3),method="ML")
summary(model4.ma2) #AIC=330.6## Series: ma2.train
## ARIMA(2,0,3) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 mean
## 0.2670 -0.1935 0.4749 0.6454 0.0524 -0.2075
## s.e. 0.4227 0.1744 0.4136 0.3050 0.2587 0.1916
##
## sigma^2 = 0.8538: log likelihood = -158.3
## AIC=330.6 AICc=331.6 BIC=350.11
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004511346 0.9006288 0.7299239 -174.2871 339.1512 0.8499178
## ACF1
## Training set 0.006394934lmtest::coeftest(model4.ma2) #hanya ma2 yang signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.267023 0.422731 0.6317 0.52761
## ar2 -0.193470 0.174380 -1.1095 0.26723
## ma1 0.474856 0.413641 1.1480 0.25097
## ma2 0.645418 0.304955 2.1164 0.03431 *
## ma3 0.052372 0.258742 0.2024 0.83960
## intercept -0.207543 0.191625 -1.0831 0.27878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model5.ma2=Arima(ma2.train, order=c(3,0,3),method="ML")
summary(model5.ma2) #AIC=329.87## Series: ma2.train
## ARIMA(3,0,3) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 mean
## -0.3810 0.2176 -0.2170 1.1432 0.7306 0.3543 -0.2053
## s.e. 0.3026 0.1535 0.1541 0.2896 0.3056 0.2426 0.1890
##
## sigma^2 = 0.8418: log likelihood = -156.93
## AIC=329.87 AICc=331.16 BIC=352.17
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004402562 0.8903241 0.70915 -674.9739 835.0008 0.8257288
## ACF1
## Training set 0.01679809lmtest::coeftest(model5.ma2) #hanya ma1 dan ma2 yang signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.38103 0.30264 -1.2590 0.2080
## ar2 0.21757 0.15355 1.4170 0.1565
## ar3 -0.21701 0.15408 -1.4084 0.1590
## ma1 1.14318 0.28962 3.9472 7.909e-05 ***
## ma2 0.73062 0.30555 2.3911 0.0168 *
## ma3 0.35428 0.24258 1.4605 0.1442
## intercept -0.20530 0.18897 -1.0864 0.2773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model yang dipilih adalah model 1, yaitu ARIMA(0,0,2)Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(0,0,2) dan parameter model ARIMA(0,0,2) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(0,0,2).
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.ma2 <- model1.ma2$residuals
par(mfrow=c(2,2))
qqnorm(sisaan.ma2)
qqline(sisaan.ma2, col = "blue", lwd = 2)
plot(c(1:length(sisaan.ma2)),sisaan.ma2)
acf(sisaan.ma2)
pacf(sisaan.ma2)
par(mfrow = c(1,1))Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan menyebar normal mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Akan tetapi, plot ACF dan PACF sisaan ARIMA(0,0,2) signifikan pada lag ke-6 sehingga sisaan tidak saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.
Uji Formal
#1) Sisaan Menyebar Normal
ks.test(sisaan.ma2,"pnorm") ##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan.ma2
## D = 0.043157, p-value = 0.9788
## alternative hypothesis: two-sided#tak tolak H0 > sisaan 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 0.9788 yang lebih besar dari taraf nyata 5% sehingga tak 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.ma2, type = "Ljung") ##
## Box-Ljung test
##
## data: sisaan.ma2
## X-squared = 0.43785, df = 1, p-value = 0.5082#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.5082 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.ma2)^2, type = "Ljung") ##
## Box-Ljung test
##
## data: (sisaan.ma2)^2
## X-squared = 2.4711, df = 1, p-value = 0.116#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.116 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.
#4) Nilai tengah sisaan sama dengan nol
t.test(sisaan.ma2, mu = 0, conf.level = 0.95) ##
## One Sample t-test
##
## data: sisaan.ma2
## t = -0.050994, df = 119, p-value = 0.9594
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1694719 0.1609621
## sample estimates:
## mean of x
## -0.004254881#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.9594 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini berbeda dengan eksplorasi.
Overfitting
Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(0,0,2) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,0,2) dan ARIMA(0,0,3).
#---OVERFITTING---#
model1a.ma2=Arima(ma2.train, order=c(1,0,2),method="ML")
summary(model1a.ma2) #327.31## Series: ma2.train
## ARIMA(1,0,2) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 mean
## 0.1990 0.5451 0.5029 -0.2083
## s.e. 0.1412 0.1239 0.0836 0.2093
##
## sigma^2 = 0.8453: log likelihood = -158.65
## AIC=327.31 AICc=327.84 BIC=341.25
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004448537 0.9039253 0.7240135 -642.4808 808.051 0.8430358
## ACF1
## Training set 0.007803042lmtest::coeftest(model1a.ma2) #ar1 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.199039 0.141222 1.4094 0.1587
## ma1 0.545109 0.123914 4.3991 1.087e-05 ***
## ma2 0.502936 0.083609 6.0153 1.795e-09 ***
## intercept -0.208276 0.209264 -0.9953 0.3196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model1b.ma2=Arima(ma2.train, order=c(0,0,3),method="ML")
summary(model1b.ma2) #327.24## Series: ma2.train
## ARIMA(0,0,3) with non-zero mean
##
## Coefficients:
## ma1 ma2 ma3 mean
## 0.7415 0.656 0.1286 -0.2084
## s.e. 0.0878 0.118 0.0997 0.2067
##
## sigma^2 = 0.8446: log likelihood = -158.62
## AIC=327.24 AICc=327.77 BIC=341.18
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004403296 0.9035497 0.7261616 -518.1374 683.3745 0.8455371
## ACF1
## Training set 0.007516805lmtest::coeftest(model1b.ma2) #ma3 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.741488 0.087816 8.4436 < 2.2e-16 ***
## ma2 0.656016 0.118034 5.5579 2.731e-08 ***
## ma3 0.128647 0.099655 1.2909 0.1967
## intercept -0.208355 0.206722 -1.0079 0.3135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model yang dipilih adalah model awal, yaitu ARIMA(0,0,2)Berdasarkan kedua model hasil overfitting di atas, model ARIMA(1,0,2) dan ARIMA(0,0,3) memiliki AIC yang lebih besar dibandingkan dengan model ARIMA(0,0,2) dan parameter kedua model ARIMA(1,0,2) dan ARIMA(0,0,3) tidak seluruhnya signifikan. Oleh karena itu, model ARIMA(0,0,2) akan tetap digunakan untuk melakukan peramalan.
Peramalan
Peramalan dilakukan menggunakan fungsi forecast() .
Contoh peramalan berikut ini dilakukan untuk 30 hari ke depan.
#---FORECAST---#
ramalan <- forecast::forecast(model1.ma2, h = 30)
ramalan## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 121 -0.5618735 -1.743227 0.6194795 -2.368597 1.244850
## 122 -0.3214588 -1.760315 1.1173975 -2.522000 1.879082
## 123 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 124 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 125 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 126 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 127 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 128 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 129 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 130 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 131 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 132 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 133 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 134 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 135 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 136 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 137 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 138 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 139 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 140 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 141 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 142 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 143 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 144 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 145 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 146 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 147 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 148 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 149 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187
## 150 -0.2082361 -1.779750 1.3632776 -2.611659 2.195187data.ramalan <- ramalan$mean
plot(ramalan)
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(0,0,2) cenderung meningkat di awal periode dan stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.
perbandingan<-matrix(data=c(ma2.test, data.ramalan),
nrow = 30, ncol = 2)
colnames(perbandingan)<-c("Aktual","Hasil Forecast")
perbandingan## Aktual Hasil Forecast
## [1,] -0.28783640 -0.5618735
## [2,] -0.64432058 -0.3214588
## [3,] 1.32792863 -0.2082361
## [4,] -1.28736112 -0.2082361
## [5,] -0.51769396 -0.2082361
## [6,] -0.04258320 -0.2082361
## [7,] 1.61687328 -0.2082361
## [8,] 2.11543361 -0.2082361
## [9,] 0.36342978 -0.2082361
## [10,] 0.02669402 -0.2082361
## [11,] 0.06572631 -0.2082361
## [12,] 2.41852142 -0.2082361
## [13,] 2.45760770 -0.2082361
## [14,] 1.78313790 -0.2082361
## [15,] -0.22320740 -0.2082361
## [16,] 0.73861434 -0.2082361
## [17,] -0.09887757 -0.2082361
## [18,] -0.98072544 -0.2082361
## [19,] -2.10125376 -0.2082361
## [20,] -1.02837992 -0.2082361
## [21,] 0.09943497 -0.2082361
## [22,] 1.53945885 -0.2082361
## [23,] -0.38934482 -0.2082361
## [24,] 1.28958884 -0.2082361
## [25,] -0.74746774 -0.2082361
## [26,] 1.06698088 -0.2082361
## [27,] 0.48742994 -0.2082361
## [28,] -0.09244800 -0.2082361
## [29,] 0.68000318 -0.2082361
## [30,] -1.17136609 -0.2082361accuracy(data.ramalan, ma2.test)## ME RMSE MAE MPE MAPE
## Test set 0.5059314 1.242369 0.9656338 103.712 151.7162Data Asli
Digunakan data kurs yang dalam hal ini hanya digunakan data 500 periode awal
datakurs<-read.csv("https://raw.githubusercontent.com/rizkynurhambali/Praktikum-MPDW-2324/main/Pertemuan%2067/kurs.csv")
datakurs <- datakurs[1:500,]
datakurs.ts<-ts(datakurs)Eksplorasi Data
Plot Data Penuh
plot.ts(datakurs.ts, lty=1, xlab="waktu", ylab="Kurs", main="Plot Data Kurs")
Berdasarkan plot data deret waktu, terlihat bahwa data cenderung memiliki trend yang naik. Berdasarkan pola data, pembagian data latih dan data uji ditetapkan dengan proporsi 86%:14%.
Plot Data Latih
kurstrain<-datakurs[1:430]
train.ts<-ts(kurstrain)
plot.ts(train.ts, lty=1, xlab="waktu", ylab="Kurs", main="Plot Kurs Train")
Berdasarkan plot data deret waktu pada data latih, terlihat bahwa data cenderung memiliki trend yang naik dan cenderung tidak bergerak pada nilai tengah tertentu. Hal ini mengindikasikan bahwa data tidak stasioner dalam rataan.
Plot Data Uji
kurstest<-datakurs[431:500]
test.ts<-ts(kurstest)
plot.ts(test.ts, lty=1, xlab="waktu", ylab="Kurs", main="Plot Kurs Test")
Uji Stasioneritas Data
Plot ACF
acf(train.ts)
Berdasarkan plot ACF, terlihat bahwa plot ACF data menurun secara perlahan (tails of slowly). Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan
Uji ADF
tseries::adf.test(train.ts)##
## Augmented Dickey-Fuller Test
##
## data: train.ts
## Dickey-Fuller = -2.0526, Lag order = 7, p-value = 0.5553
## 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.5553 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:430)
bc = boxcox(train.ts~index, lambda = seq(5,10,by=1))
#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 6.616162#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]## [1] 5.656566 5.707071 5.757576 5.808081 5.858586 5.909091 5.959596 6.010101
## [9] 6.060606 6.111111 6.161616 6.212121 6.262626 6.313131 6.363636 6.414141
## [17] 6.464646 6.515152 6.565657 6.616162 6.666667 6.717172 6.767677 6.818182
## [25] 6.868687 6.919192 6.969697 7.020202 7.070707 7.121212 7.171717 7.222222
## [33] 7.272727 7.323232 7.373737 7.424242 7.474747 7.525253 7.575758 7.626263
## [41] 7.676768Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 6,64 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0,48 dan batas atas 5,27. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.
Penanganan Ketidakstasioneran Data
train.diff<-diff(train.ts,differences = 1)
plot.ts(train.diff, lty=1, xlab="waktu", ylab="Data Difference 1 Kurs", main="Plot Difference Kurs")
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.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(train.diff)## Warning in tseries::adf.test(train.diff): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: train.diff
## Dickey-Fuller = -6.3673, Lag order = 7, 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.diff)
Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 1, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,1).
Plot PACF
pacf(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(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 x o o o o o o o o o o o o
## 2 x x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x x o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o x x o o o o o o o o
## 7 x x o x x x x 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,2), ARIMA(1,1,2), ARIMA(2,1,2), dan ARIMA(3,1,2).
Pendugaan Parameter Model Tentatif
ARIMA(0,1,1)
model1.da=Arima(train.diff, order=c(0,1,1),method="ML")
summary(model1.da) #AIC=4753.18## Series: train.diff
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.9879
## s.e. 0.0148
##
## sigma^2 = 3835: log likelihood = -2374.59
## AIC=4753.18 AICc=4753.21 BIC=4761.3
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.6762224 61.78387 43.22489 NaN Inf 0.7363284 0.112352lmtest::coeftest(model1.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.987941 0.014806 -66.727 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,0)
model2.da=Arima(train.diff, order=c(1,1,0),method="ML")
summary(model2.da) #AIC=4917.41## Series: train.diff
## ARIMA(1,1,0)
##
## Coefficients:
## ar1
## -0.3940
## s.e. 0.0444
##
## sigma^2 = 5676: log likelihood = -2456.7
## AIC=4917.41 AICc=4917.43 BIC=4925.52
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.05787648 75.16292 54.08917 NaN Inf 0.9213996 -0.121856lmtest::coeftest(model2.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.394037 0.044414 -8.8719 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,1)
model3.da=Arima(train.diff, order=c(1,1,1),method="ML")
summary(model3.da) #AIC=4761.39## Series: train.diff
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.1207 -0.9966
## s.e. 0.0491 0.0210
##
## sigma^2 = 3782: log likelihood = -2371.6
## AIC=4749.21 AICc=4749.27 BIC=4761.39
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.082868 61.28528 43.2876 -Inf Inf 0.7373967 0.006544589lmtest::coeftest(model3.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.120655 0.049064 2.4591 0.01393 *
## ma1 -0.996606 0.020970 -47.5253 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(0,1,2)
model4.da=Arima(train.diff, order=c(0,1,2),method="ML")
summary(model4.da) #AIC=4748.3## Series: train.diff
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.8600 -0.1359
## s.e. 0.0512 0.0507
##
## sigma^2 = 3775: log likelihood = -2371.15
## AIC=4748.3 AICc=4748.35 BIC=4760.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.059629 61.22216 43.24997 -Inf Inf 0.7367557 -0.009277841lmtest::coeftest(model4.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.860021 0.051204 -16.796 < 2.2e-16 ***
## ma2 -0.135945 0.050725 -2.680 0.007361 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,2)
model5.da=Arima(train.diff, order=c(1,1,2),method="ML")
summary(model5.da) #AIC=4749.85 ## Series: train.diff
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.1723 -0.6912 -0.3021
## s.e. 0.2443 0.2348 0.2305
##
## sigma^2 = 3782: log likelihood = -2370.92
## AIC=4749.85 AICc=4749.94 BIC=4766.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.798497 61.21036 43.15664 -Inf Inf 0.7351659 -0.006549048lmtest::coeftest(model5.da) #terdapat parameter tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.17233 0.24425 -0.7056 0.480467
## ma1 -0.69120 0.23481 -2.9436 0.003244 **
## ma2 -0.30214 0.23054 -1.3106 0.189995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(2,1,2)
model6.da=Arima(train.diff, order=c(2,1,2),method="ML")
summary(model6.da) #AIC=4749.52## Series: train.diff
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.4635 -0.1308 -1.3295 0.3386
## s.e. 0.3091 0.0565 0.3089 0.3073
##
## sigma^2 = 3777: log likelihood = -2369.76
## AIC=4749.52 AICc=4749.66 BIC=4769.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.4957896 61.09641 43.1911 NaN Inf 0.7357528 -0.004047784lmtest::coeftest(model6.da) #terdapat parameter tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.463451 0.309061 1.4995 0.13373
## ar2 -0.130837 0.056523 -2.3148 0.02063 *
## ma1 -1.329455 0.308856 -4.3045 1.674e-05 ***
## ma2 0.338625 0.307281 1.1020 0.27046
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(0,1,2) dan parameter model ARIMA(0,1,2) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(0,1,2).
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 <- model4.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 tidak menyebar normal ditandai dengan titik titik yang cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung tidak sama menandakan bahwa sisaan memiliki ragam yang heterogen. Plot ACF dan PACF sisaan ARIMA(0,0,2) juga tidak signifikan pada 20 lag awal yang menandakan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.
Uji Formal
#1) Sisaan Menyebar Normal
ks.test(sisaan.da,"pnorm") #tak tolak H0 > sisaan menyebar normal##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan.da
## D = 0.46816, p-value < 2.2e-16
## alternative hypothesis: two-sidedSelain 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 0.00 yang kurang 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.da, type = "Ljung") #tak tolak H0 > sisaan saling bebas##
## Box-Ljung test
##
## data: sisaan.da
## X-squared = 0.037186, df = 1, p-value = 0.8471Selanjutnya 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.8471 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.da)^2, type = "Ljung") #tak tolak H0 > sisaan homogen##
## Box-Ljung test
##
## data: (sisaan.da)^2
## X-squared = 29.906, df = 1, p-value = 4.536e-08Hipotesis 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.000 yang kurang dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak 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 = -0.69638, df = 428, p-value = 0.4866
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.872876 3.753619
## sample estimates:
## mean of x
## -2.059629Terakhir, 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.4866 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini berbeda dengan eksplorasi.
Peramalan
Peramalan dilakukan menggunakan fungsi forecast() .
Contoh peramalan berikut ini dilakukan untuk 30 hari ke depan.
#---FORECAST---#
ramalan.da <- forecast::forecast(model4.da, h = 30)
ramalan.da## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 431 0.3239709 -78.42499 79.07294 -120.1122 120.7601
## 432 6.7730086 -72.74745 86.29346 -114.8430 128.3891
## 433 6.7730086 -72.74808 86.29410 -114.8440 128.3900
## 434 6.7730086 -72.74872 86.29473 -114.8450 128.3910
## 435 6.7730086 -72.74935 86.29537 -114.8459 128.3920
## 436 6.7730086 -72.74998 86.29600 -114.8469 128.3929
## 437 6.7730086 -72.75062 86.29664 -114.8479 128.3939
## 438 6.7730086 -72.75125 86.29727 -114.8489 128.3949
## 439 6.7730086 -72.75189 86.29790 -114.8498 128.3958
## 440 6.7730086 -72.75252 86.29854 -114.8508 128.3968
## 441 6.7730086 -72.75316 86.29917 -114.8518 128.3978
## 442 6.7730086 -72.75379 86.29981 -114.8527 128.3988
## 443 6.7730086 -72.75442 86.30044 -114.8537 128.3997
## 444 6.7730086 -72.75506 86.30108 -114.8547 128.4007
## 445 6.7730086 -72.75569 86.30171 -114.8556 128.4017
## 446 6.7730086 -72.75633 86.30234 -114.8566 128.4026
## 447 6.7730086 -72.75696 86.30298 -114.8576 128.4036
## 448 6.7730086 -72.75760 86.30361 -114.8586 128.4046
## 449 6.7730086 -72.75823 86.30425 -114.8595 128.4055
## 450 6.7730086 -72.75886 86.30488 -114.8605 128.4065
## 451 6.7730086 -72.75950 86.30552 -114.8615 128.4075
## 452 6.7730086 -72.76013 86.30615 -114.8624 128.4085
## 453 6.7730086 -72.76077 86.30678 -114.8634 128.4094
## 454 6.7730086 -72.76140 86.30742 -114.8644 128.4104
## 455 6.7730086 -72.76204 86.30805 -114.8653 128.4114
## 456 6.7730086 -72.76267 86.30869 -114.8663 128.4123
## 457 6.7730086 -72.76330 86.30932 -114.8673 128.4133
## 458 6.7730086 -72.76394 86.30996 -114.8683 128.4143
## 459 6.7730086 -72.76457 86.31059 -114.8692 128.4152
## 460 6.7730086 -72.76521 86.31122 -114.8702 128.4162data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(0,012) cenderung stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.
pt_1 <- train.ts[430] #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=30), hasil[-1]),
nrow = 30, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da## Aktual Hasil Forecast
## [1,] 12866 12813.32
## [2,] 12887 12820.10
## [3,] 12862 12826.87
## [4,] 12863 12833.64
## [5,] 12993 12840.42
## [6,] 12962 12847.19
## [7,] 12963 12853.96
## [8,] 13022 12860.74
## [9,] 12983 12867.51
## [10,] 13047 12874.28
## [11,] 13059 12881.05
## [12,] 13164 12887.83
## [13,] 13176 12894.60
## [14,] 13191 12901.37
## [15,] 13237 12908.15
## [16,] 13209 12914.92
## [17,] 13164 12921.69
## [18,] 13008 12928.47
## [19,] 13075 12935.24
## [20,] 13076 12942.01
## [21,] 12972 12948.78
## [22,] 12932 12955.56
## [23,] 13003 12962.33
## [24,] 13064 12969.10
## [25,] 13086 12975.88
## [26,] 13084 12982.65
## [27,] 13043 12989.42
## [28,] 13000 12996.20
## [29,] 12942 13002.97
## [30,] 12982 13009.74accuracy(ts(hasil[-1]), head(test.ts, n=30))## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 118.9674 156.6908 126.4518 0.9079581 0.9657541 0.7929994 2.646555