5. Time Series & Forecasting
Carolin
2023-09-30
1 Metode Smoothing
1.1 Data
setwd("D:/Kuliah/Mat/TSA Kominfo/Praktikum")
df <- read.csv('5. Moving Average.csv')
head(df,15)## t xt
## 1 1 86
## 2 2 90
## 3 3 86
## 4 4 102
## 5 5 111
## 6 6 92
## 7 7 96
## 8 8 106
## 9 9 102
## 10 10 104
## 11 11 99
## 12 12 92
## 13 13 87
## 14 14 85
## 15 15 94df.ts.1 = ts(df$xt)plot(df.ts.1, main= "Time Series Plot",
xlab = "Period", ylab = "Yt")
points(df.ts.1)
# terlihat data memiliki pola yang konstanlibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(TTR)
library(graphics)1.2 Single Moving Average
Single Moving Average dengan n = 3
df.ts.sma <-SMA(df.ts.1, n = 3)
ramal.sma <- c(NA, df.ts.sma)
df.sma <- cbind(df.aktual = c(df.ts.1,NA),
pemulusan = c(df.ts.sma,NA), ramal.sma)
data.frame(df.sma)## df.aktual pemulusan ramal.sma
## 1 86 NA NA
## 2 90 NA NA
## 3 86 87.33333 NA
## 4 102 92.66667 87.33333
## 5 111 99.66667 92.66667
## 6 92 101.66667 99.66667
## 7 96 99.66667 101.66667
## 8 106 98.00000 99.66667
## 9 102 101.33333 98.00000
## 10 104 104.00000 101.33333
## 11 99 101.66667 104.00000
## 12 92 98.33333 101.66667
## 13 87 92.66667 98.33333
## 14 85 88.00000 92.66667
## 15 94 88.66667 88.00000
## 16 112 97.00000 88.66667
## 17 128 111.33333 97.00000
## 18 101 113.66667 111.33333
## 19 124 117.66667 113.66667
## 20 124 116.33333 117.66667
## 21 118 122.00000 116.33333
## 22 120 120.66667 122.00000
## 23 124 120.66667 120.66667
## 24 130 124.66667 120.66667
## 25 108 120.66667 124.66667
## 26 109 115.66667 120.66667
## 27 128 115.00000 115.66667
## 28 124 120.33333 115.00000
## 29 106 119.33333 120.33333
## 30 122 117.33333 119.33333
## 31 119 115.66667 117.33333
## 32 125 122.00000 115.66667
## 33 127 123.66667 122.00000
## 34 105 119.00000 123.66667
## 35 103 111.66667 119.00000
## 36 105 104.33333 111.66667
## 37 93 100.33333 104.33333
## 38 107 101.66667 100.33333
## 39 103 101.00000 101.66667
## 40 110 106.66667 101.00000
## 41 94 102.33333 106.66667
## 42 102 102.00000 102.33333
## 43 94 96.66667 102.00000
## 44 96 97.33333 96.66667
## 45 106 98.66667 97.33333
## 46 106 102.66667 98.66667
## 47 118 110.00000 102.66667
## 48 109 111.00000 110.00000
## 49 NA NA 111.00000# Plot
ts.plot(df.ts.1, xlab = "Periode", ylab = "YT", col = "blue", lty = 2)
points(df.ts.1)
lines(df.ts.sma, col = "red", lwd=2)
lines(ramal.sma, col = "black",lwd = 2)
title("Rataan Bergerak Sederhana N = 3",
cex.main = 1, font.main = 2, col.main = "black")
legend("topleft", c("Data Aktual", "Pemulusan SMA", "Ramalan SMA"),
lty = 1:3, col = c("blue", "red", "black"))
#Ukuran Keakuratan
error.sma <- df.sma[,1] - df.sma[,3]
## SSE : SUM SQUARE ERROR
sse.sma<- sum(error.sma^2, na.rm = T) # ketika ada nilai NA pada data
# agar NA dibuang otomatis dari perhitungan
# MSE: MEAN SQUARED ERROR
mse.sma <- mean(error.sma^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse.sma <- sqrt(mse.sma)
# MAD: mean absolute deviation
mad.sma <- mean(abs(error.sma), na.rm = T)
## MAPE
r.error.sma <- (error.sma/df.sma[,1])*100 # relative error
mape.sma<- mean(abs(r.error.sma), na.rm = T)
data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse.sma, mse.sma, mape.sma, rmse.sma,
mad.sma))## Ukuran.Keakuratan Nilai
## 1 SSE 5380.111111
## 2 MSE 119.558025
## 3 MAPE 8.319077
## 4 RMSE 10.934259
## 5 MAD 8.9555561.3 Double Moving Average
Double Moving Average dengan n = 3
df.ts.dma <-SMA(df.ts.sma, n = 3)
At <- 2*df.ts.sma - df.ts.dma
Bt <- df.ts.sma - df.ts.dma
pemulusan.dma <- At + Bt
ramal.dma <- c(NA, pemulusan.dma)
df.dma <- cbind(df.aktual = c(df.ts.1,NA),
pemulusan.dma= c(pemulusan.dma,NA), ramal.dma)
data.frame(df.dma)## df.aktual pemulusan.dma ramal.dma
## 1 86 NA NA
## 2 90 NA NA
## 3 86 NA NA
## 4 102 NA NA
## 5 111 112.55556 NA
## 6 92 109.00000 112.55556
## 7 96 98.33333 109.00000
## 8 106 94.44444 98.33333
## 9 102 104.66667 94.44444
## 10 104 109.77778 104.66667
## 11 99 100.33333 109.77778
## 12 92 92.33333 100.33333
## 13 87 82.88889 92.33333
## 14 85 78.00000 82.88889
## 15 94 86.44444 78.00000
## 16 112 108.55556 86.44444
## 17 128 136.00000 108.55556
## 18 101 126.33333 136.00000
## 19 124 124.55556 126.33333
## 20 124 117.22222 124.55556
## 21 118 128.66667 117.22222
## 22 120 122.66667 128.66667
## 23 124 119.77778 122.66667
## 24 130 130.00000 119.77778
## 25 108 118.00000 130.00000
## 26 109 106.33333 118.00000
## 27 128 110.77778 106.33333
## 28 124 127.00000 110.77778
## 29 106 121.55556 127.00000
## 30 122 114.00000 121.55556
## 31 119 112.11111 114.00000
## 32 125 129.33333 112.11111
## 33 127 130.11111 129.33333
## 34 105 113.88889 130.11111
## 35 103 98.77778 113.88889
## 36 105 89.66667 98.77778
## 37 93 90.11111 89.66667
## 38 107 100.77778 90.11111
## 39 103 101.00000 100.77778
## 40 110 113.77778 101.00000
## 41 94 100.33333 113.77778
## 42 102 98.66667 100.33333
## 43 94 89.33333 98.66667
## 44 96 94.66667 89.33333
## 45 106 100.88889 94.66667
## 46 106 108.88889 100.88889
## 47 118 122.44444 108.88889
## 48 109 117.22222 122.44444
## 49 NA NA 117.22222# Plot Perbandingan SMA dan DMA
ts.plot(df.ts.1, xlab = "Periode", ylab = "YT",
col = "blue", lty = 2, ylim = c(77,138))
points(df.ts.1)
lines(pemulusan.dma, col = "red", lwd=2)
lines(df.ts.sma, col = "black",lwd = 2)
title("Perbandingan SMA dan DMA dengan N = 3",
cex.main = 1, font.main = 2, col.main = "black")
legend("topleft", c("Data Aktual", "Pemulusan SMA", "Pemulusan DMA"),
lty = 1:3, col = c("blue", "red", "black"))
#Ukuran Keakuratan
error.dma <- df.dma[,1] - df.dma[,3]
## SSE : SUM SQUARE ERROR
sse.dma<- sum(error.dma^2, na.rm = T) # ketika ada nilai NA pada data
# agar NA dibuang otomatis dari perhitungan
# MSE: MEAN SQUARED ERROR
mse.dma <- mean(error.dma^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse.dma <- sqrt(mse.dma)
# MAD: mean absolute deviation
mad.dma <- mean(abs(error.dma), na.rm = T)
## MAPE
r.error.dma <- (error.dma/df.dma[,1])*100 # relative error
mape.dma<- mean(abs(r.error.dma), na.rm = T)
# Perbandingan SMA dan DMA
akurasi <- data.frame(
"Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Simple Moving Average N = 3" = c(sse.sma, mse.sma, mape.sma, rmse.sma, mad.sma),
"Double Moving Average N = 3" = c(sse.dma, mse.dma, mape.dma, rmse.dma, mad.dma)
)
akurasi## Ukuran.Keakuratan Simple.Moving.Average.N...3 Double.Moving.Average.N...3
## 1 SSE 5380.111111 7521.925926
## 2 MSE 119.558025 174.928510
## 3 MAPE 8.319077 9.804028
## 4 RMSE 10.934259 13.226054
## 5 MAD 8.955556 10.4392761.4 Exponential
1.4.1 Single Exponential Smoothing (SES)
df2 <- read.csv('5. Exponential Smoothing.csv')
head(df2, 15)## t yt
## 1 1 9.640
## 2 2 7.596
## 3 3 9.684
## 4 4 6.723
## 5 5 7.508
## 6 6 8.135
## 7 7 8.269
## 8 8 9.489
## 9 9 9.986
## 10 10 5.455
## 11 11 5.843
## 12 12 7.939
## 13 13 8.828
## 14 14 8.574
## 15 15 7.291Plot Time Series
df2.ts.1 = ts(df2$yt)
plot(df2.ts.1, main= "Time Series Plot",
xlab = "Period", ylab = "Yt")
points(df2.ts.1)
# terlihat data memiliki pola tren SES dengan alpha = 0.2
df2.ses <- HoltWinters(df2.ts.1, alpha = 0.2, beta = F, gamma = F)
df2.ses## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = df2.ts.1, alpha = 0.2, beta = F, gamma = F)
##
## Smoothing parameters:
## alpha: 0.2
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 2.760377# diperoleha nilai S_0 = yaitu 2.670377df2.ses$SSE # Nilai Keakuratan## [1] 115.5729datases <- data.frame(df2.ts.1, c(NA, df2.ses$fitted[,1]))
colnames(datases) = c("y", "yhat")
head(datases)## y yhat
## 1 9.640 NA
## 2 7.596 9.640000
## 3 9.684 9.231200
## 4 6.723 9.321760
## 5 7.508 8.802008
## 6 8.135 8.543206ramal2.ses <- forecast::forecast(df2.ses, h = 5) # periode ramalan sebanyak 5
(df2.ramal.ses <- data.frame(ramal2.ses))## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## 37 2.760377 0.7729830 4.747771 -0.2790797 5.799833
## 38 2.760377 0.7336249 4.787129 -0.3392728 5.860026
## 39 2.760377 0.6950166 4.825737 -0.3983190 5.919073
## 40 2.760377 0.6571169 4.863637 -0.4562816 5.977035
## 41 2.760377 0.6198882 4.900865 -0.5132180 6.033972plot(df2.ts.1,
xlab = "Period", ylab = "Yt", col = "blue", lty = 2)
points(df2.ts.1)
lines(datases[,2], col = "red", lwd = 2) # nilai dugaan
title("Single Exponential Smoothing Alpha = 0.2", cex.main = 1,
font.main = 4, col.main = "black")
legend("bottomleft", c("Data Aktual","Fitted SES"), lty = 1:3,
col = c("blue", "red"))
plot(ramal2.ses, xlab = "periode", ylab = "Yt")
Gabungan Data Aktual, Pemulusan dan Ramalan
data.ses<-cbind(aktual = c(df2.ts.1, rep(NA, 5)),
pemulusan = c(df2.ses$fitted[,1], as.numeric(df2.ses$coefficients),
rep(NA,5)),
ramalan = c(NA, df2.ses$fitted[,1], df2.ramal.ses$Point.Forecast))
data.ses <-ts(data.ses)
#Plot
plot(data.ses[,1],
xlba = "Time Period",
ylab = "Yt",
main = "FOrecasting Using Simple Exponential Smoothing",
lty = 3, col = "black")## Warning in plot.window(xlim, ylim, log, ...): "xlba" is not a graphical
## parameter## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "xlba" is not a
## graphical parameter## Warning in axis(1, ...): "xlba" is not a graphical parameter## Warning in axis(2, ...): "xlba" is not a graphical parameter## Warning in box(...): "xlba" is not a graphical parameterpoints(data.ses[,1])
lines(data.ses[,2], lwd = 3, lty = 3, col = "blue")
lines(data.ses[,3], lwd = 3, lty = 3, col = "blue")
legend("bottomleft", legend=c("data aktual", "data pemulusan",
"data peramalan"),
lty = c(2,2,3),
col = c("black", "blue", "red"), cex =0.8)
Ukuran Keakuratan
error.ses <- data.ses[,1] - data.ses[,3]
## SSE : SUM SQUARE ERROR
sse.ses <- sum(error.ses^2, na.rm = T)
# MSE: MEAN SQUARED ERROR
mse.ses <- mean(error.ses^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse.ses <- sqrt(mse.ses)
# MAD: mean absolute deviation
mad.ses <- mean(abs(error.ses), na.rm = T)
## MAPE
r.error.ses <- (error.ses/data.ses[,1])*100 # relative error
mape.ses<- mean(abs(r.error.ses), na.rm = T)
data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse.ses, mse.ses, mape.ses, rmse.ses, mad.ses))## Ukuran.Keakuratan Nilai
## 1 SSE 115.572927
## 2 MSE 3.302084
## 3 MAPE 52.462591
## 4 RMSE 1.817164
## 5 MAD 1.4824041.4.2 Double Exponential Smoothing(DES)
# DES alpha = 0.2 dan beta = 0.3
df2.des <- HoltWinters(df2.ts.1, alpha = 0.2, beta = 0.3, gamma = F)
head(df2.des)## $fitted
## Time Series:
## Start = 3
## End = 36
## Frequency = 1
## xhat level trend
## 3 5.552000 7.596000 -2.044000000
## 4 4.582320 6.378400 -1.796080000
## 5 3.342817 5.010456 -1.667639200
## 6 2.758125 4.175853 -1.417728208
## 7 2.738384 3.833500 -1.095115722
## 8 3.081229 3.844508 -0.763278790
## 9 3.983971 4.362783 -0.378812517
## 10 5.165686 5.184376 -0.018690747
## 11 5.222217 5.223549 -0.001331887
## 12 5.382288 5.346373 0.035915115
## 13 6.082949 5.893631 0.189317809
## 14 6.985980 6.631959 0.354020896
## 15 7.752886 7.303584 0.449302112
## 16 8.082098 7.660509 0.421588958
## 17 8.336221 7.953278 0.382943097
## 18 8.344707 8.048177 0.296529823
## 19 7.455453 7.432565 0.022887414
## 20 7.420243 7.410762 0.009480242
## 21 6.874300 6.992994 -0.118694310
## 22 6.758127 6.876240 -0.118112293
## 23 5.824362 6.130702 -0.306339942
## 24 5.316168 5.669090 -0.352921665
## 25 5.471243 5.706934 -0.235691743
## 26 5.248488 5.481194 -0.232706301
## 27 5.184135 5.377990 -0.193855569
## 28 4.392764 4.724508 -0.331743650
## 29 3.931602 4.293211 -0.361609495
## 30 3.055296 3.535681 -0.480385602
## 31 2.270113 2.820837 -0.550723351
## 32 1.448960 2.062091 -0.613130149
## 33 1.635601 2.064168 -0.428567779
## 34 1.953597 2.209880 -0.256283816
## 35 1.812598 2.042277 -0.229679617
## 36 1.750203 1.941278 -0.191075480
##
## $x
## Time Series:
## Start = 1
## End = 36
## Frequency = 1
## [1] 9.640 7.596 9.684 6.723 7.508 8.135 8.269 9.489 9.986 5.455 5.843 7.939
## [13] 8.828 8.574 7.291 7.438 6.896 3.784 7.232 5.284 6.884 3.621 5.048 7.270
## [25] 5.521 5.896 2.886 3.895 1.952 1.883 1.230 4.525 4.507 2.397 2.456 0.578
##
## $alpha
## [1] 0.2
##
## $beta
## [1] 0.3
##
## $gamma
## [1] FALSE
##
## $coefficients
## a b
## 1.5157622 -0.2614076df2.des$SSE## [1] 267.0453datades <- data.frame(df2.ts.1, c(NA, NA, df2.des$fitted[,1]))
colnames(datades) = c("y","yhat")
head(datades)## y yhat
## 1 9.640 NA
## 2 7.596 NA
## 3 9.684 5.552000
## 4 6.723 4.582320
## 5 7.508 3.342817
## 6 8.135 2.758125ramal.des <- forecast::forecast(df2.des, h = 5)
(df2.ramal.des <- data.frame(ramal.des))## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## 37 1.2543545 -2.209531 4.718240 -4.043201 6.551910
## 38 0.9929469 -2.586103 4.571997 -4.480738 6.466632
## 39 0.7315392 -3.015226 4.478304 -4.998643 6.461722
## 40 0.4701316 -3.501119 4.441382 -5.603372 6.543635
## 41 0.2087239 -4.044950 4.462398 -6.296709 6.714157# Plot
ts.plot(datades[,1],xlab = "periode", ylab = "Yt", col = "blue", lty = 3)
points(datades[,1])
lines(datades[,2], col = "red", lwd = 2)
title("Double Exponential Smoothing (Alpha = 0.2 dan beta = 0.3)", cex.main = 1,
font.main = 4, col.main = "black")
legend("topright", c("Data Aktual", "Fitted DES"), lty = c(3,1), col = c("blue","red"))
plot(ramal.des, xlab = "periode waktu", ylab = "Yt")
1.5 Winter’s Method
df3 <- read.csv("5. Winter.csv")
head(df3,15)## t yt
## 1 1 11354
## 2 2 11273
## 3 3 17050
## 4 4 18007
## 5 5 20288
## 6 6 21619
## 7 7 18749
## 8 8 20429
## 9 9 19689
## 10 10 17904
## 11 11 16940
## 12 12 14799
## 13 13 12993
## 14 14 13664
## 15 15 19013plot(df3$t, df3$yt, type = 'b')
df3.ts = ts(df3$yt, frequency = 12)
# perlu ditambahkan pola seasonal itu terjadi setiap periode berapa kali
# time series plot
plot(df3.ts, main = "Time Series Plot", xlab = "Period", ylab = "Yt")
points(df3.ts)
1.5.1 Winter Aditif
Melakukan Winter Aditif Method
aditif.hw <- HoltWinters(df3.ts)## Warning in HoltWinters(df3.ts): optimization difficulties: ERROR:
## ABNORMAL_TERMINATION_IN_LNSRCHhead(aditif.hw)## $fitted
## xhat level trend season
## Jan 2 12421.06 17515.14 68.54283 -5162.62500
## Feb 2 13119.20 17602.26 70.31546 -4553.37500
## Mar 2 18524.14 17690.26 72.00398 761.87500
## Apr 2 17635.79 17778.14 73.51910 -215.87500
## May 2 22265.15 17865.18 74.80907 4325.16667
## Jun 2 21491.14 17951.67 75.92434 3463.54167
## Jul 2 19360.77 18036.69 76.79171 1247.29167
## Aug 2 20937.81 18102.67 75.76036 2759.37500
## Sep 2 20116.85 18200.98 77.91188 1837.95833
## Oct 2 18296.89 18252.23 75.36782 -30.70833
## Nov 2 17323.31 18340.66 76.61408 -1093.95833
## Dec 2 15140.59 18403.91 75.33930 -3338.66667
## Jan 3 13534.49 18475.01 74.93457 -5015.45307
## Feb 3 14260.91 18594.88 79.22248 -4413.18656
## Mar 3 19651.77 18683.94 80.16184 887.66733
## Apr 3 18730.87 18759.89 79.75964 -108.77555
## May 3 23333.26 18836.08 79.41911 4417.76174
## Jun 3 22504.00 18891.33 77.11242 3535.55499
## Jul 3 20242.40 19000.56 80.17764 1161.66371
## Aug 3 22101.75 19083.32 80.42434 2938.00397
## Sep 3 20900.55 19190.81 83.00683 1626.73919
## Oct 3 19418.78 19263.96 82.06603 72.76210
## Nov 3 18184.33 19305.89 78.23596 -1199.79686
## Dec 3 16141.38 19430.94 82.70414 -3372.26898
## Jan 4 14930.91 19508.18 82.18228 -4659.45085
## Feb 4 15310.71 19566.04 79.86117 -4335.19638
## Mar 4 20580.37 19646.20 79.88998 854.27486
## Apr 4 19681.23 19737.32 80.96119 -137.04801
## May 4 24144.48 19835.61 82.61550 4226.24925
## Jun 4 23765.33 19894.90 80.38872 3790.04404
## Jul 4 21240.90 19978.10 80.65733 1182.14589
## Aug 4 23299.80 20066.04 81.35188 3152.41458
## Sep 4 21783.10 20152.62 81.85148 1548.62927
## Oct 4 20050.60 20215.76 80.06596 -245.22840
## Nov 4 19548.97 20297.56 80.23146 -828.82772
## Dec 4 17020.20 20357.50 78.29449 -3415.59615
##
## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1 11354 11273 17050 18007 20288 21619 18749 20429 19689 17904 16940 14799
## 2 12993 13664 19013 18052 22625 21771 19028 21632 19296 18699 16912 15010
## 3 14918 14564 19522 18621 22589 23493 20322 22935 20597 18183 19626 15973
## 4 14182 15320 20926 20215 23426 23852 21465 23461 21207 20104 18924 16885
##
## $alpha
## alpha
## 0.03247517
##
## $beta
## beta
## 0.09543618
##
## $gamma
## gamma
## 0.2659558
##
## $coefficients
## a b s1 s2 s3
## 20431.4039631 77.8754698 -4852.1598952 -4332.8049068 943.2121603
## s4 s5 s6 s7 s8
## 0.3005512 4041.3710965 3812.3459320 1239.8101247 3193.8939442
## s9 s10 s11 s12
## 1400.3871781 -231.4881004 -989.6438250 -3450.3852427# ciri aditif: lebarnya dari waktu ke waktu konstan
# multiplikatif: lebarnya akan mengalami perubahan dari waktu ke waktu aditif.hw$SSE # nilai keakuratan## [1] 13160801range(df3$yt)## [1] 11273 23852data.frame(aditif.hw$fitted)## xhat level trend season
## 1 12421.06 17515.14 68.54283 -5162.62500
## 2 13119.20 17602.26 70.31546 -4553.37500
## 3 18524.14 17690.26 72.00398 761.87500
## 4 17635.79 17778.14 73.51910 -215.87500
## 5 22265.15 17865.18 74.80907 4325.16667
## 6 21491.14 17951.67 75.92434 3463.54167
## 7 19360.77 18036.69 76.79171 1247.29167
## 8 20937.81 18102.67 75.76036 2759.37500
## 9 20116.85 18200.98 77.91188 1837.95833
## 10 18296.89 18252.23 75.36782 -30.70833
## 11 17323.31 18340.66 76.61408 -1093.95833
## 12 15140.59 18403.91 75.33930 -3338.66667
## 13 13534.49 18475.01 74.93457 -5015.45307
## 14 14260.91 18594.88 79.22248 -4413.18656
## 15 19651.77 18683.94 80.16184 887.66733
## 16 18730.87 18759.89 79.75964 -108.77555
## 17 23333.26 18836.08 79.41911 4417.76174
## 18 22504.00 18891.33 77.11242 3535.55499
## 19 20242.40 19000.56 80.17764 1161.66371
## 20 22101.75 19083.32 80.42434 2938.00397
## 21 20900.55 19190.81 83.00683 1626.73919
## 22 19418.78 19263.96 82.06603 72.76210
## 23 18184.33 19305.89 78.23596 -1199.79686
## 24 16141.38 19430.94 82.70414 -3372.26898
## 25 14930.91 19508.18 82.18228 -4659.45085
## 26 15310.71 19566.04 79.86117 -4335.19638
## 27 20580.37 19646.20 79.88998 854.27486
## 28 19681.23 19737.32 80.96119 -137.04801
## 29 24144.48 19835.61 82.61550 4226.24925
## 30 23765.33 19894.90 80.38872 3790.04404
## 31 21240.90 19978.10 80.65733 1182.14589
## 32 23299.80 20066.04 81.35188 3152.41458
## 33 21783.10 20152.62 81.85148 1548.62927
## 34 20050.60 20215.76 80.06596 -245.22840
## 35 19548.97 20297.56 80.23146 -828.82772
## 36 17020.20 20357.50 78.29449 -3415.59615data.winter.aditif <- data.frame(ts(df3$yt[13:length(df3$yt)],
frequency = 12), aditif.hw$fitted[,1])
colnames(data.winter.aditif) = c("y", "yhat")
data.winter.aditif## y yhat
## 1 12993 12421.06
## 2 13664 13119.20
## 3 19013 18524.14
## 4 18052 17635.79
## 5 22625 22265.15
## 6 21771 21491.14
## 7 19028 19360.77
## 8 21632 20937.81
## 9 19296 20116.85
## 10 18699 18296.89
## 11 16912 17323.31
## 12 15010 15140.59
## 13 14918 13534.49
## 14 14564 14260.91
## 15 19522 19651.77
## 16 18621 18730.87
## 17 22589 23333.26
## 18 23493 22504.00
## 19 20322 20242.40
## 20 22935 22101.75
## 21 20597 20900.55
## 22 18183 19418.78
## 23 19626 18184.33
## 24 15973 16141.38
## 25 14182 14930.91
## 26 15320 15310.71
## 27 20926 20580.37
## 28 20215 19681.23
## 29 23426 24144.48
## 30 23852 23765.33
## 31 21465 21240.90
## 32 23461 23299.80
## 33 21207 21783.10
## 34 20104 20050.60
## 35 18924 19548.97
## 36 16885 17020.20ramal.aditif <- forecast::forecast(aditif.hw, h = 12)
head(ramal.aditif)## $method
## [1] "HoltWinters"
##
## $model
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = df3.ts)
##
## Smoothing parameters:
## alpha: 0.03247517
## beta : 0.09543618
## gamma: 0.2659558
##
## Coefficients:
## [,1]
## a 20431.4039631
## b 77.8754698
## s1 -4852.1598952
## s2 -4332.8049068
## s3 943.2121603
## s4 0.3005512
## s5 4041.3710965
## s6 3812.3459320
## s7 1239.8101247
## s8 3193.8939442
## s9 1400.3871781
## s10 -231.4881004
## s11 -989.6438250
## s12 -3450.3852427
##
## $level
## [1] 80 95
##
## $mean
## Jan Feb Mar Apr May Jun Jul Aug
## 5 15657.12 16254.35 21608.24 20743.21 24862.15 24711.00 22216.34 24248.30
## Sep Oct Nov Dec
## 5 22532.67 20978.67 20298.39 17915.52
##
## $lower
## 80% 95%
## Jan 5 14878.82 14466.81
## Feb 5 15475.56 15063.29
## Mar 5 20828.87 20416.29
## Apr 5 19963.15 19550.22
## May 5 24081.32 23667.97
## Jun 5 23929.28 23515.46
## Jul 5 21433.61 21019.25
## Aug 5 23464.43 23049.48
## Sep 5 21747.53 21331.91
## Oct 5 20192.13 19775.76
## Nov 5 19510.30 19093.11
## Dec 5 17125.73 16707.64
##
## $upper
## 80% 95%
## Jan 5 16435.42 16847.43
## Feb 5 17033.14 17445.41
## Mar 5 22387.62 22800.19
## Apr 5 21523.26 21936.19
## May 5 25642.99 26056.33
## Jun 5 25492.73 25906.55
## Jul 5 22999.08 23413.43
## Aug 5 25032.17 25447.13
## Sep 5 23317.81 23733.43
## Oct 5 21765.21 22181.58
## Nov 5 21086.48 21503.67
## Dec 5 18705.32 19123.41plot(aditif.hw, xlab = "Periode Waktu", lwd = 2, lyt = 2)## Warning in plot.window(xlim, ylim, log, ...): "lyt" is not a graphical
## parameter## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "lyt" is not a
## graphical parameter## Warning in axis(1, ...): "lyt" is not a graphical parameter## Warning in axis(2, ...): "lyt" is not a graphical parameter## Warning in box(...): "lyt" is not a graphical parameterpoints(df3.ts)
legend("bottomright", c("Data Aktual", "Fitted Aditif Winter's"), lty = c(2,1),
col = c("black", "maroon"))
plot(ramal.aditif, xlab = "Periode Waktu", ylab = "Yt")
dt <- ts(data.winter.aditif)
error <- dt[,1] - dt[,2]
## SSE : SUM SQUARE ERROR
sse<- sum(error^2, na.rm = T)
# MSE: MEAN SQUARED ERROR
mse<- mean(error^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse<- sqrt(mse)
# MAD: mean absolute deviation
mad<- mean(abs(error), na.rm = T)
## MAPE
r.error<- (error/dt[,1])*100 #relative error
mape<- mean(abs(r.error), na.rm = T)
akurasi.wa <- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse, mse, mape, rmse, mad))
format(akurasi.wa, scientific = F)## Ukuran.Keakuratan Nilai
## 1 SSE 13160801.09707
## 2 MSE 365577.80825
## 3 MAPE 2.59711
## 4 RMSE 604.63031
## 5 MAD 483.133241.5.2 Winter’s Multiplikatif
multi.hw <- HoltWinters(df3.ts, seasonal = "mult")
multi.hw## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = df3.ts, seasonal = "mult")
##
## Smoothing parameters:
## alpha: 0.03486983
## beta : 0.0268336
## gamma: 0.3030974
##
## Coefficients:
## [,1]
## a 2.030983e+04
## b 7.362873e+01
## s1 7.408672e-01
## s2 7.692830e-01
## s3 1.052334e+00
## s4 1.002149e+00
## s5 1.215205e+00
## s6 1.205450e+00
## s7 1.068809e+00
## s8 1.174645e+00
## s9 1.075766e+00
## s10 9.900418e-01
## s11 9.511360e-01
## s12 8.214672e-01multi.hw$SSE # nilai keakuratan## [1] 13910813data.frame(multi.hw$fitted)## xhat level trend season
## 1 12577.93 17515.14 68.54283 0.7153183
## 2 13249.39 17603.91 69.08577 0.7496969
## 3 18494.38 17692.28 69.60324 1.0412394
## 4 17629.77 17779.26 70.06928 0.9876995
## 5 22162.53 17864.23 70.46928 1.2357345
## 6 21416.75 17947.75 70.81945 1.1885929
## 7 19384.79 18028.96 71.09833 1.0709794
## 8 20988.38 18088.44 70.78661 1.1557968
## 9 20121.51 18178.65 71.30766 1.1025513
## 10 18255.91 18223.85 70.60708 0.9978931
## 11 17258.64 18309.94 71.02255 0.9389410
## 12 15037.86 18368.09 70.67711 0.8155568
## 13 13366.93 18437.57 70.64515 0.7222156
## 14 14114.06 18583.11 72.65468 0.7565521
## 15 19682.93 18676.50 73.21116 1.0497725
## 16 18716.08 18744.37 73.06772 0.9946136
## 17 23481.78 18814.10 72.97827 1.2432723
## 18 22614.06 18862.04 72.30637 1.1943409
## 19 20274.13 18960.01 72.99496 1.0652093
## 20 22282.41 19034.57 73.03700 1.1661538
## 21 20915.30 19127.12 73.56062 1.0893002
## 22 19359.56 19190.49 73.28721 1.0049721
## 23 18010.48 19222.95 72.19177 0.9334205
## 24 15837.12 19355.50 73.81120 0.8151148
## 25 14566.11 19435.12 73.96718 0.7466321
## 26 14939.55 19491.15 73.48581 0.7635995
## 27 20584.92 19582.01 73.95200 1.0472609
## 28 19606.05 19667.32 74.25674 0.9931353
## 29 24389.09 19762.95 74.83046 1.2294263
## 30 24018.60 19810.47 74.09748 1.2079018
## 31 21269.57 19879.76 73.96842 1.0659450
## 32 23562.98 19960.12 74.13997 1.1761345
## 33 21803.15 20031.23 74.05883 1.0844481
## 34 19898.96 20086.12 73.54446 0.9870676
## 35 19387.24 20166.91 73.73883 0.9578365
## 36 16585.96 20223.79 73.28631 0.8171601data.winter.multliplikatif <- data.frame(ts(df3$yt[13:length(df3$yt)],
frequency = 12), multi.hw$fitted[,1])
colnames(data.winter.multliplikatif) = c("y", "yhat")
data.winter.multliplikatif## y yhat
## 1 12993 12577.93
## 2 13664 13249.39
## 3 19013 18494.38
## 4 18052 17629.77
## 5 22625 22162.53
## 6 21771 21416.75
## 7 19028 19384.79
## 8 21632 20988.38
## 9 19296 20121.51
## 10 18699 18255.91
## 11 16912 17258.64
## 12 15010 15037.86
## 13 14918 13366.93
## 14 14564 14114.06
## 15 19522 19682.93
## 16 18621 18716.08
## 17 22589 23481.78
## 18 23493 22614.06
## 19 20322 20274.13
## 20 22935 22282.41
## 21 20597 20915.30
## 22 18183 19359.56
## 23 19626 18010.48
## 24 15973 15837.12
## 25 14182 14566.11
## 26 15320 14939.55
## 27 20926 20584.92
## 28 20215 19606.05
## 29 23426 24389.09
## 30 23852 24018.60
## 31 21465 21269.57
## 32 23461 23562.98
## 33 21207 21803.15
## 34 20104 19898.96
## 35 18924 19387.24
## 36 16885 16585.96(ramal.multi <- forecast::forecast(multi.hw, h = 12))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 5 15101.44 14841.53 15361.35 14703.95 15498.93
## Feb 5 15737.29 15475.72 15998.87 15337.24 16137.34
## Mar 5 21605.18 21339.13 21871.22 21198.30 22012.06
## Apr 5 20648.63 20381.48 20915.77 20240.07 21057.18
## May 5 25127.98 24854.79 25401.17 24710.17 25545.79
## Jun 5 25015.03 24739.95 25290.11 24594.33 25435.72
## Jul 5 22258.20 21984.55 22531.85 21839.69 22676.71
## Aug 5 24548.75 24269.80 24827.71 24122.13 24975.38
## Sep 5 22561.50 22283.55 22839.44 22136.42 22986.57
## Oct 5 20836.54 20559.34 21113.75 20412.59 21260.49
## Nov 5 20087.76 19809.77 20365.74 19662.62 20512.89
## Dec 5 17409.67 17319.12 17500.22 17271.18 17548.15# Plot
plot(multi.hw, xlab = "Periode Waktu", ylab = "Yt", lwd = 2, lty = 2)
points(df3.ts)
legend("bottomright", c("Data Aktual", "Fitted Multiplicative Winter's"), lty = c(2,1),
col = c("black", "red"), cex = 0.8)
plot(ramal.multi, xlab = "Periode Waktu", ylab = "Yt")
# Ukuran Keakuratan
dt <- ts(data.winter.multliplikatif)
error <- dt[,1] - dt[,2]
## SSE : SUM SQUARE ERROR
sse<- sum(error^2, na.rm = T)
# MSE: MEAN SQUARED ERROR
mse<- mean(error^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse<- sqrt(mse)
# MAD: mean absolute deviation
mad<- mean(abs(error), na.rm = T)
## MAPE
r.error<- (error/dt[,1])*100 #relative error
mape<- mean(abs(r.error), na.rm = T)
akurasi.wm <- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Winter's Multiplikatif" = c(sse, mse, mape, rmse, mad))
format(akurasi.wa, scientific = F)## Ukuran.Keakuratan Nilai
## 1 SSE 13160801.09707
## 2 MSE 365577.80825
## 3 MAPE 2.59711
## 4 RMSE 604.63031
## 5 MAD 483.133241.5.3 Perbandingan Aditif dan Multiplikatif
ts.plot(df3.ts, xlab = "Periode Waktu", ylab = "Yt", col = "blue",lty = 3, ylim = c(10100, 24050))
points(df3.ts)
lines(data.winter.aditif[,2], col = "red", lwd = 2)
lines(data.winter.multliplikatif[,2], col = "black", lwd = 2)
title("Perbandingan Winters Aditif dan Multiplikatif",
cex.main = 1, font.main = 4, col.main = "black")
legend("bottomright", c("Data Aktual", "Ramalan Aditif", "Ramalan Multlipikatif"),
lty = 1:3, col = c("blue", "red", "black"), cex = 0.8)
# Perbandingan Keakuratan Winter's Aditif dan Multiplikatif
cbind(akurasi.wa, "Winter's Multiplikatif" = akurasi.wm$Winter.s.Multiplikatif)## Ukuran.Keakuratan Nilai Winter's Multiplikatif
## 1 SSE 1.316080e+07 1.391081e+07
## 2 MSE 3.655778e+05 3.864115e+05
## 3 MAPE 2.597110e+00 2.650933e+00
## 4 RMSE 6.046303e+02 6.216200e+02
## 5 MAD 4.831332e+02 4.975390e+021.6 Time Series Decomposition
1.6.1 Dekomosisi Aditif
df.dekomposisi <- decompose(df3.ts, type = "additive")
plot(df.dekomposisi)
1.6.2 Dekomposisi Multiplikatif
df.dekomposisi.multi <- decompose(df3.ts, type = "multiplicative")
plot(df.dekomposisi.multi)
1.7 Data 2
library(readxl)setwd("D:/Kuliah/Mat/TSA Kominfo/Praktikum")
lat <- read_xlsx("5. Pertemuan1_latihan.xlsx", sheet = "Sheet1")
head(lat,15)## # A tibble: 15 × 2
## Period yt
## <dbl> <dbl>
## 1 1 48.7
## 2 2 45.8
## 3 3 46.4
## 4 4 46.2
## 5 5 44
## 6 6 53.8
## 7 7 47.6
## 8 8 47
## 9 9 47.6
## 10 10 51.1
## 11 11 49.1
## 12 12 46.7
## 13 13 47.8
## 14 14 45.8
## 15 15 45.5Membentuk Objek Time Series
lat.ts = ts(lat$yt)Plot Time Series
plot(lat.ts, main= "Time Series Plot",
xlab = "Period", ylab = "Yt")
points(lat.ts)
1.7.1 Single Moving Average
lat.ts.sma <- SMA(lat.ts, n = 3)
ramal.sma <- c(NA, lat.ts.sma)
df.sma <- cbind(df.aktual = c(lat.ts, NA),
pemulusan = c(lat.ts.sma, NA), ramal.sma)
data.frame(df.sma)## df.aktual pemulusan ramal.sma
## 1 48.7 NA NA
## 2 45.8 NA NA
## 3 46.4 46.96667 NA
## 4 46.2 46.13333 46.96667
## 5 44.0 45.53333 46.13333
## 6 53.8 48.00000 45.53333
## 7 47.6 48.46667 48.00000
## 8 47.0 49.46667 48.46667
## 9 47.6 47.40000 49.46667
## 10 51.1 48.56667 47.40000
## 11 49.1 49.26667 48.56667
## 12 46.7 48.96667 49.26667
## 13 47.8 47.86667 48.96667
## 14 45.8 46.76667 47.86667
## 15 45.5 46.36667 46.76667
## 16 49.2 46.83333 46.36667
## 17 54.8 49.83333 46.83333
## 18 44.7 49.56667 49.83333
## 19 51.1 50.20000 49.56667
## 20 47.3 47.70000 50.20000
## 21 45.3 47.90000 47.70000
## 22 43.3 45.30000 47.90000
## 23 44.6 44.40000 45.30000
## 24 47.1 45.00000 44.40000
## 25 53.4 48.36667 45.00000
## 26 44.9 48.46667 48.36667
## 27 50.5 49.60000 48.46667
## 28 48.1 47.83333 49.60000
## 29 45.4 48.00000 47.83333
## 30 51.6 48.36667 48.00000
## 31 50.8 49.26667 48.36667
## 32 46.4 49.60000 49.26667
## 33 52.3 49.83333 49.60000
## 34 50.5 49.73333 49.83333
## 35 53.4 52.06667 49.73333
## 36 53.9 52.60000 52.06667
## 37 52.3 53.20000 52.60000
## 38 53.0 53.06667 53.20000
## 39 48.6 51.30000 53.06667
## 40 52.4 51.33333 51.30000
## 41 47.9 49.63333 51.33333
## 42 49.5 49.93333 49.63333
## 43 44.0 47.13333 49.93333
## 44 53.8 49.10000 47.13333
## 45 52.5 50.10000 49.10000
## 46 52.0 52.76667 50.10000
## 47 50.6 51.70000 52.76667
## 48 48.7 50.43333 51.70000
## 49 51.4 50.23333 50.43333
## 50 47.7 49.26667 50.23333
## 51 NA NA 49.26667Plot
ts.plot(lat.ts, xlab = "Periode", ylab = "Yt",
col = "blue", lty = 2)
points(lat.ts)
lines(lat.ts.sma, col = "red", lwd = 2)
lines(ramal.sma, col = "black", lwd = 2)
title("Rataan Bergerak Sederhana N = 3", cex.main = 1,
font.main = 2, col.main = "black")
legend("topleft", c("Data Aktual", "Pemulusan SMA",
"Ramalan SMA"), lty = 1:3,
col = c("blue", "red","black"))
Ukuran Keakuratan
error.sma <- df.sma[,-1] - df.sma[,-3]
## SSE : SUM SQUARE ERROR
sse.sma <- sum(error.sma^2, na.rm = T)
# MSE: MEAN SQUARED ERROR
mse.sma <- mean(error.sma^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse.sma <- sqrt(mse.sma)
# MAD: mean absolute deviation
mad.sma <- mean(abs(error.sma), na.rm = T)
## MAPE
r.error.sma <- (error.sma/df.sma[,-1])*100
mape.sma <- mean(abs(r.error.sma), na.rm = T)
akurasi <- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse.sma, mse.sma, mape.sma, rmse.sma, mad.sma))
akurasi## Ukuran.Keakuratan Nilai
## 1 SSE 357.902222
## 2 MSE 3.767392
## 3 MAPE 3.013992
## 4 RMSE 1.940977
## 5 MAD 1.4694741.7.2 Double Moving Average (DMA)
lat.ts.dma <-SMA(lat.ts.sma, n = 3)
At <- 2*lat.ts.sma - lat.ts.dma
Bt <- lat.ts.sma - lat.ts.dma
pemulusan.dma <- At + Bt
ramal.dma <- c(NA, pemulusan.dma)
df.dma <- cbind(df.akutal = c(lat.ts, NA),
pemulusan.dma = c(pemulusan.dma, NA), ramal.dma)
data.frame(df.dma)## df.akutal pemulusan.dma ramal.dma
## 1 48.7 NA NA
## 2 45.8 NA NA
## 3 46.4 NA NA
## 4 46.2 NA NA
## 5 44.0 44.17778 NA
## 6 53.8 50.88889 44.17778
## 7 47.6 50.73333 50.88889
## 8 47.0 51.11111 50.73333
## 9 47.6 45.31111 51.11111
## 10 51.1 48.74444 45.31111
## 11 49.1 50.97778 48.74444
## 12 46.7 49.03333 50.97778
## 13 47.8 46.20000 49.03333
## 14 45.8 44.56667 46.20000
## 15 45.5 45.10000 44.56667
## 16 49.2 47.18889 45.10000
## 17 54.8 54.14444 47.18889
## 18 44.7 51.21111 54.14444
## 19 51.1 50.86667 51.21111
## 20 47.3 44.78889 50.86667
## 21 45.3 46.50000 44.78889
## 22 43.3 41.96667 46.50000
## 23 44.6 41.46667 41.96667
## 24 47.1 45.20000 41.46667
## 25 53.4 53.25556 45.20000
## 26 44.9 50.84444 53.25556
## 27 50.5 51.17778 50.84444
## 28 48.1 46.23333 51.17778
## 29 45.4 47.04444 46.23333
## 30 51.6 48.96667 47.04444
## 31 50.8 50.71111 48.96667
## 32 46.4 50.64444 50.71111
## 33 52.3 50.36667 50.64444
## 34 50.5 49.75556 50.36667
## 35 53.4 55.11111 49.75556
## 36 53.9 54.86667 55.11111
## 37 52.3 54.35556 54.86667
## 38 53.0 53.28889 54.35556
## 39 48.6 48.85556 53.28889
## 40 52.4 50.20000 48.85556
## 41 47.9 47.38889 50.20000
## 42 49.5 49.20000 47.38889
## 43 44.0 43.60000 49.20000
## 44 53.8 49.85556 43.60000
## 45 52.5 52.74444 49.85556
## 46 52.0 56.98889 52.74444
## 47 50.6 52.05556 56.98889
## 48 48.7 48.03333 52.05556
## 49 51.4 49.12222 48.03333
## 50 47.7 47.84444 49.12222
## 51 NA NA 47.84444ts.plot(lat.ts, xlab = "Periode", ylab = "Yt",
col = "blue", lty = 2, ylim = c(77, 138))
points(lat.ts)
lines(pemulusan.dma, col = "red", lwd = 2)
lines(ramal.dma, col = "black", lwd = 2)
title("Rataan Bergerak Berganda N = 3", cex.main = 1,
font.main = 2, col.main = "black")
legend("topleft", c("Data Aktual", "Pemulusan",
"Ramalan DMA"), lty = 1:3,
col = c("blue", "red","black"))
Ukuran Keakuratan
error.dma <- df.dma[,1] - df.dma[,3]
## SSE : SUM SQUARE ERROR
sse.dma <- sum(error.dma^2, na.rm = T)
# MSE: MEAN SQUARED ERROR
mse.dma <- mean(error.dma^2, na.rm = T)
# RMSE: Root Mean Square Error
rmse.dma <- sqrt(mse.dma)
# MAD: mean absolute deviation
mad.dma <- mean(abs(error.dma), na.rm = T)
## MAPE
r.error.dma <- (error.dma/df.dma[,1])*100 # relative error
mape.dma <- mean(abs(r.error.dma), na.rm = T)
akurasi <- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse.dma, mse.dma, mape.dma, rmse.dma, mad.dma))
akurasi## Ukuran.Keakuratan Nilai
## 1 SSE 874.330370
## 2 MSE 19.429564
## 3 MAPE 7.126068
## 4 RMSE 4.407898
## 5 MAD 3.5111111.7.3 Perbandingan SMA dan DMA
ts.plot(lat.ts, xlab = "Periode Waktu", ylab = "Yt",
col = "blue", lty = 2, ylim = c(77, 138))
points(lat.ts)
lines(pemulusan.dma, col = "red", lwd = 2)
lines(lat.ts.sma, col = "black", lwd = 2)
title("RPerbandingan SMA dan DMA", cex.main = 1,
font.main = 2, col.main = "black")
legend("topleft", c("Data Aktual", "Pemulusan SMA",
"Pemulusan DMA"), lty = 1:3,
col = c("blue", "red","black"))
Akurasi Perbandingan SMA dan DMA
akurasi <- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Simple Moving Average N = 3" = c(sse.sma,
mse.sma, mape.sma,
rmse.sma, mad.sma),
"Double Moving Average N = 3" = c(sse.dma, mse.dma,
mape.dma, rmse.dma, mad.dma))
akurasi## Ukuran.Keakuratan Simple.Moving.Average.N...3 Double.Moving.Average.N...3
## 1 SSE 357.902222 874.330370
## 2 MSE 3.767392 19.429564
## 3 MAPE 3.013992 7.126068
## 4 RMSE 1.940977 4.407898
## 5 MAD 1.469474 3.5111112 ARIMA
MODEL ARIMA (Auto - Regressive Integrated Moving Average)
library(TTR)
library(tseries)
library(forecast)
library(TSA)## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tar2.1 Deret Stationer
setwd("D:/Kuliah/Mat/TSA Kominfo/Praktikum")
dt.arima <- read.csv('5. Bagi Hasil.csv')
head(dt.arima, 15)## t Xt
## 1 1 120.8399
## 2 2 120.4160
## 3 3 122.3064
## 4 4 121.5694
## 5 5 122.8942
## 6 6 123.3013
## 7 7 124.1932
## 8 8 124.8930
## 9 9 125.8147
## 10 10 125.8902
## 11 11 126.6755
## 12 12 124.1835
## 13 13 125.5775
## 14 14 124.7087
## 15 15 124.9495dt.ar.ts <- ts(dt.arima$Xt)2.1.1 Memeriksa Kestasioneran
# Melalui Plot TS
plot(dt.ar.ts, main = "Plot Bagi Hasil", ylab = expression(Y[t]))
Data ini merupakan data dengan pola tren, namun datanya belum stasioner.
Tidak memiliki pola musiman!!
df.dekomposisi <- decompose(ts(dt.arima$Xt, frequency = 2)) # jika tak ada frequency, maka error
plot(df.dekomposisi)
# Melalui Uji Augmented Dickey-Fuller (ADF)
tseries::adf.test(dt.ar.ts, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: dt.ar.ts
## Dickey-Fuller = -3.8475, Lag order = 4, p-value = 0.0194
## alternative hypothesis: stationaryHipotesis pada uji ADF: H0: Data mempunyai unit root (tidak stasioner) H1: Data tidak mempunyai unit root (stasioner) \(p-value = 0.0194<0.05\) maka data sudah stasioner
# Melalui Plot ACF
acf(dt.ar.ts, main = "Plot ACF Bagi Hasil")
2.1.2 Identifikasi Model
par(mfrow = c(1,2))
Acf(dt.ar.ts)
Pacf(dt.ar.ts)
Kandidat model untuk data ialah AR(1) yang diperoleh dari identifikasi plot PACF yang signifikan pada lag ke-1 dan cut off pada lag ke-2
# Identifikasi Kandidat Model Menggunakan Extended ACF
TSA::eacf(dt.ar.ts)## 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 o 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 o o o o o o o o o o o o
## 5 o x o o o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x o o x o o o o o o o o o o2.1.3 Pendugaan Parameter
fit_model <-Arima(dt.ar.ts, order = c(1,0,0)) #AR(1)
summary(fit_model)## Series: dt.ar.ts
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.7637 124.6449
## s.e. 0.0678 0.4257
##
## sigma^2 = 1.082: log likelihood = -145.26
## AIC=296.53 AICc=296.78 BIC=304.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04255821 1.029761 0.8427568 0.02756793 0.6753622 0.9345546
## ACF1
## Training set 0.02322355lmtest::coeftest(fit_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.763653 0.067761 11.27 < 2.2e-16 ***
## intercept 124.644911 0.425688 292.81 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit_model_2 <- Arima(dt.ar.ts, order = c(1,0,2)) #AR(1)
summary(fit_model_2)## Series: dt.ar.ts
## ARIMA(1,0,2) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 mean
## 0.5522 0.2542 0.3462 124.6899
## s.e. 0.1335 0.1389 0.1175 0.3506
##
## sigma^2 = 1.03: log likelihood = -141.89
## AIC=293.78 AICc=294.42 BIC=306.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0302255 0.9945218 0.7910991 0.01808066 0.6340434 0.8772701
## ACF1
## Training set -0.02296842lmtest::coeftest(fit_model_2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.55225 0.13347 4.1375 3.511e-05 ***
## ma1 0.25421 0.13892 1.8299 0.067268 .
## ma2 0.34618 0.11751 2.9461 0.003218 **
## intercept 124.68988 0.35058 355.6723 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cbind(AR1 = AIC(fit_model), ARIMA2 = AIC(fit_model_2))## AR1 ARIMA2
## [1,] 296.5282 293.7843cbind(AR1 = BIC(fit_model), ARIMA2 = BIC(fit_model_2))## AR1 ARIMA2
## [1,] 304.3437 306.8102karena komponen ma1 pada ARMA(1,2) tidak signifikan, yang dipilih adalah ARMA(1,0)
# Fungsi Pengujian Parameter Mode
printstatarima <- function (x, digits = 4,se=TRUE){
if (length(x$coef) > 0) {
cat("\nCoefficients:\n")
coef <- round(x$coef, digits)
if (se && nrow(x$var.coef)) {
ses <- rep(0, length(coef))
ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits)
coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
coef <- rbind(coef, s.e. = ses)
statt <- coef[1,]/ses
pval <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
coef <- rbind(coef, t=round(statt, digits),sign.=round(pval,digits))
coef <- t(coef)
}
print.default(coef, print.gap = 2)
}
}
printstatarima(fit_model)##
## Coefficients:
## s.e. t sign.
## ar1 0.7637 0.0678 11.2640 0
## intercept 124.6449 0.4257 292.7999 0lmtest::coeftest(fit_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.763653 0.067761 11.27 < 2.2e-16 ***
## intercept 124.644911 0.425688 292.81 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Intercept dan koefisien AR(1) signifikan. Sehingga persamaan AR(1) yang diperoleh adalah \(y_t - 124.6449 = 0.7637(y_{t-1} - 124.6449) + e_t\) atau \(y_t = 29.45948 + 0.7637y_{t-1} + e_t\)
2.1.4 Diagnostik Model
# Plot Sisaan dan Sisaan Terstandarisasi
par(mfrow = c(1,2))
plot(residuals(fit_model), type = "p", ylab = 'sisaan'); abline(h=0)
plot(rstandard(fit_model), ylab = "sisaan terstandarisasi"); abline(h=0)
t.test(residuals(fit_model), mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: residuals(fit_model)
## t = 0.41156, df = 99, p-value = 0.6815
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1626227 0.2477391
## sample estimates:
## mean of x
## 0.04255821Hipotesis yang diuji: \(H_0:\mu=0\) \(H_1 :\mu \ne 0\) Tidak tolak Ho. Tidak cukup bukti bahwa terjadi pelanggaran asumsi rataan tidak 0 (rataan sisaan 0)
# Plot Q-Q dan Histogram
par(mfrow = c(1,2))
h <- hist(residuals(fit_model), main = 'Histogram Residual')
xfit <- seq(min(residuals(fit_model)), max(residuals(fit_model)), length = 100)
yfit <- dnorm (xfit, mean = mean(residuals(fit_model)), sd = sd(residuals(fit_model)))
yfit <- yfit*diff(h$mids[1:2])*length(residuals(fit_model))
lines(xfit, yfit)
qqnorm(residuals(fit_model)); qqline(residuals(fit_model))
# Uji Kenormalan
shapiro.test(residuals(fit_model))##
## Shapiro-Wilk normality test
##
## data: residuals(fit_model)
## W = 0.99173, p-value = 0.8015ks.test(residuals(fit_model),'pnorm')##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuals(fit_model)
## D = 0.068131, p-value = 0.7421
## alternative hypothesis: two-sidedHipotesis yang diuji: H0: sisaan menyebar normal H1: sisaan tidak menyebar normal
Tidak tolak H0. Tidak cukup bukti bahwa terjadi pelanggaran asumsi sisaan tidak menyebar normal (sisaan menyebar normal).
# Plot ACF dan PACF Sisaan Model
par(mfrow = c(1,2))
Acf(residuals(fit_model), main = "Plot ACF Sisaan")
Pacf(residuals(fit_model), main = "Plot PACF Sisaan")
autokorelasi dalam sisaan tidak ada.
# L-JUng Box Test
Box.test(residuals(fit_model), lag = 24, type = "Ljung")##
## Box-Ljung test
##
## data: residuals(fit_model)
## X-squared = 28.26, df = 24, p-value = 0.2492Hipotesis yang diuji: H0: tidak ada autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan Tidak tolak H0. TIdak cukup bukti bahwa terjadi pelanggaran asumsi autokorelasi (sisaan saling bebas)
2.1.5 Overfitting Model
# Overfitting AR(1)
overfit_1 <- Arima(dt.ar.ts, order = c(2,0,0)) # AR(2)
summary(overfit_1)## Series: dt.ar.ts
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 0.8144 -0.0724 124.6730
## s.e. 0.0994 0.1037 0.3924
##
## sigma^2 = 1.088: log likelihood = -145.02
## AIC=298.04 AICc=298.46 BIC=308.46
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03674748 1.027288 0.8360755 0.02288207 0.6701901 0.9271456
## ACF1
## Training set -0.04265301lmtest::coeftest(overfit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.814392 0.099360 8.1964 2.477e-16 ***
## ar2 -0.072357 0.103716 -0.6976 0.4854
## intercept 124.672993 0.392436 317.6899 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1diperoleh ar2 tidak signifikan pada AR(2), maka yang dipilih adalah AR(1)
printstatarima(overfit_1)##
## Coefficients:
## s.e. t sign.
## ar1 0.8144 0.0994 8.1932 0.0000
## ar2 -0.0724 0.1037 -0.6982 0.4867
## intercept 124.6730 0.3924 317.7192 0.0000Koefisien AR(2) tidak signifikan, sehingga model terbaik adalah AR(1)
c(AR1 = AIC(fit_model), Ar2 = AIC(overfit_1))## AR1 Ar2
## 296.5282 298.04252.1.6 Peramalam
#Peramalan menggunakan model terbaik ARIMA (0,1,1)
ramal_fit_model <- forecast::forecast(fit_model, h = 30)
ramal_fit_model## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 101 123.9471 122.6140 125.2802 121.9083 125.9859
## 102 124.1120 122.4347 125.7894 121.5467 126.6773
## 103 124.2380 122.3892 126.0867 121.4105 127.0654
## 104 124.3341 122.3924 126.2759 121.3645 127.3038
## 105 124.4076 122.4136 126.4015 121.3581 127.4571
## 106 124.4637 122.4399 126.4875 121.3686 127.5588
## 107 124.5065 122.4655 126.5475 121.3851 127.6279
## 108 124.5392 122.4883 126.5902 121.4026 127.6759
## 109 124.5642 122.5075 126.6209 121.4187 127.7097
## 110 124.5833 122.5232 126.6434 121.4326 127.7339
## 111 124.5978 122.5358 126.6599 121.4442 127.7515
## 112 124.6090 122.5458 126.6722 121.4536 127.7644
## 113 124.6175 122.5536 126.6813 121.4611 127.7739
## 114 124.6240 122.5597 126.6882 121.4669 127.7810
## 115 124.6289 122.5644 126.6934 121.4716 127.7863
## 116 124.6327 122.5681 126.6973 121.4751 127.7902
## 117 124.6356 122.5709 126.7003 121.4779 127.7932
## 118 124.6378 122.5730 126.7025 121.4800 127.7955
## 119 124.6395 122.5747 126.7042 121.4817 127.7972
## 120 124.6408 122.5760 126.7055 121.4829 127.7986
## 121 124.6417 122.5769 126.7065 121.4839 127.7996
## 122 124.6425 122.5777 126.7073 121.4847 127.8003
## 123 124.6431 122.5783 126.7079 121.4852 127.8009
## 124 124.6435 122.5787 126.7083 121.4857 127.8013
## 125 124.6438 122.5790 126.7086 121.4860 127.8017
## 126 124.6441 122.5793 126.7089 121.4863 127.8019
## 127 124.6443 122.5795 126.7091 121.4864 127.8021
## 128 124.6444 122.5796 126.7092 121.4866 127.8023
## 129 124.6445 122.5797 126.7093 121.4867 127.8024
## 130 124.6446 122.5798 126.7094 121.4868 127.8025plot(ramal_fit_model)
2.2 Deret Tidak Stasioner (Nilai Tengah)
Data merupakan data suhu tanah global yang diperoleh dari Data NCEI.NOAA.GOV
wages = read.csv("https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land/ann/2/1850-2022/data.csv")
head(wages, 15)## Global.Land.January...December.Temperature.Anomalies
## Units: Degrees Celsius
## Base Period: 1901-2000
## Missing: -999
## Year Anomaly
## 1850 -0.38
## 1851 -0.26
## 1852 -0.35
## 1853 -0.34
## 1854 -0.27
## 1855 -0.35
## 1856 -0.55
## 1857 -0.49
## 1858 -0.41
## 1859 -0.27
## 1860 -0.56wages_2 <- as.numeric(wages[35:140,])
data.ts.2 <- ts(wages_2)2.2.1 Memeriksa Kestasioneran
plot(data.ts.2, ylab = expression(Y[t]))
points(data.ts.2)
par(mfrow = c(1,2))
plot(as.numeric(wages[5:dim(wages)[1],]))
plot(as.numeric(wages[5:dim(wages)[1],]), type = 'l')
adf.test(data.ts.2, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: data.ts.2
## Dickey-Fuller = -3.0978, Lag order = 4, p-value = 0.1214
## alternative hypothesis: stationaryKarena p-value > 0.05 maka tidak tolak H0. Artinya deret tidak stationer.
acf(data.ts.2, main = 'Plot ACF Data')
Karena deret tidak stasioner maka deret harus ditranformasi. Berdasarkan plot deret tampak bahwa deret tidak stasioner dalam nilai tengahnya. Sehingga transformasi yang digunakan adalah pembedaan (differencing)
data.ts.diff = diff(data.ts.2)
#Cek kestasioneran data hasil differencing
adf.test(data.ts.diff)## Warning in adf.test(data.ts.diff): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: data.ts.diff
## Dickey-Fuller = -7.7677, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryplot(data.ts.diff, main = 'Hasil Differencing Data')
Setelah dilakukan pembedaan satu kali diperoleh hasil uji ADF p-alue <0.05, sehingga deret sudah stationer.
acf(data.ts.diff, main = 'Plot ACF Hasil Differencing')
2.2.2 Identifikasi Model
par(mfrow = c(1,2))
Acf(data.ts.diff, main = 'Plot ACF')
Pacf(data.ts.diff, main = 'Plot PACF')
Kandidat model untuk data ialah AR(1) yang diperoleh dari identifikasi plot PACF yang signifikan pada lag ke-1 dan cut off pada lag ke-2
# Identifikasi Kandidat Model Menggunakan Extended ACF
TSA::eacf(data.ts.diff)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x x x x x o x x o o o o
## 1 x o x 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 o x o o o o o o o o o o o o
## 4 x x o x o o o o o o o o o o
## 5 o x o x o o o o o o o o o o
## 6 o x o o x o o o o o o o o o
## 7 o x o o x o o o o o o o o oauto.arima(data.ts.2)## Series: data.ts.2
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -0.7506 0.0073
## s.e. 0.0693 0.0038
##
## sigma^2 = 0.02393: log likelihood = 47.58
## AIC=-89.15 AICc=-88.91 BIC=-81.192.2.3 Pendugaan Parameter
model_arma_1 <-Arima(data.ts.2, order = c(3,1,2)) #AR(1)
summary(model_arma_1)## Series: data.ts.2
## ARIMA(3,1,2)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2
## -0.7036 0.0310 -0.2218 0.1573 -0.5708
## s.e. 0.1475 0.1933 0.1188 0.1288 0.1190
##
## sigma^2 = 0.02241: log likelihood = 52.33
## AIC=-92.66 AICc=-91.8 BIC=-76.73
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.0210547 0.1453965 0.1155169 Inf Inf 0.7815257 -0.03028915lmtest::coeftest(model_arma_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.703582 0.147482 -4.7706 1.836e-06 ***
## ar2 0.030963 0.193271 0.1602 0.87272
## ar3 -0.221791 0.118802 -1.8669 0.06192 .
## ma1 0.157290 0.128831 1.2209 0.22213
## ma2 -0.570820 0.118994 -4.7970 1.610e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1printstatarima(model_arma_1)##
## Coefficients:
## s.e. t sign.
## ar1 -0.7036 0.1475 -4.7702 0.0000
## ar2 0.0310 0.1933 0.1604 0.8729
## ar3 -0.2218 0.1188 -1.8670 0.0647
## ma1 0.1573 0.1288 1.2213 0.2247
## ma2 -0.5708 0.1190 -4.7966 0.0000model_arma_2 <- Arima(data.ts.2, order = c(3,1,0)) #AR(1)
summary(model_arma_2)## Series: data.ts.2
## ARIMA(3,1,0)
##
## Coefficients:
## ar1 ar2 ar3
## -0.5514 -0.3188 -0.3685
## s.e. 0.0904 0.1006 0.0912
##
## sigma^2 = 0.02405: log likelihood = 47.89
## AIC=-87.78 AICc=-87.38 BIC=-77.16
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01388612 0.1521235 0.1162104 Inf Inf 0.7862175 -0.03443426lmtest::coeftest(model_arma_2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.551446 0.090440 -6.0974 1.078e-09 ***
## ar2 -0.318837 0.100648 -3.1679 0.001536 **
## ar3 -0.368504 0.091216 -4.0399 5.347e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1printstatarima(model_arma_2)##
## Coefficients:
## s.e. t sign.
## ar1 -0.5514 0.0904 -6.0996 0e+00
## ar2 -0.3188 0.1006 -3.1690 2e-03
## ar3 -0.3685 0.0912 -4.0406 1e-04model_arma_3 <- Arima(data.ts.2, order = c(0,1,1)) #AR(1)
summary(model_arma_3)## Series: data.ts.2
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.6987
## s.e. 0.0703
##
## sigma^2 = 0.0244: log likelihood = 46.12
## AIC=-88.24 AICc=-88.12 BIC=-82.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.0219713 0.1547257 0.121763 Inf Inf 0.823783 0.02320088lmtest::coeftest(model_arma_3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.698673 0.070288 -9.9402 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1printstatarima(model_arma_3)##
## Coefficients:
## s.e. t sign.
## ma1 -0.6987 0.0703 -9.9388 0model_arma_4 <- Arima(data.ts.2, order = c(2,1,2))
summary(model_arma_4)## Series: data.ts.2
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.6774 0.2886 0.0945 -0.6983
## s.e. 0.1436 0.1360 0.1016 0.0796
##
## sigma^2 = 0.02289: log likelihood = 50.79
## AIC=-91.57 AICc=-90.96 BIC=-78.3
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.02226008 0.1476903 0.1169867 Inf Inf 0.7914696 0.006142657lmtest::coeftest(model_arma_4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.677351 0.143596 -4.7171 2.393e-06 ***
## ar2 0.288551 0.135986 2.1219 0.03384 *
## ma1 0.094488 0.101616 0.9298 0.35245
## ma2 -0.698313 0.079597 -8.7731 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1printstatarima(model_arma_4)##
## Coefficients:
## s.e. t sign.
## ar1 -0.6774 0.1436 -4.7173 0.0000
## ar2 0.2886 0.1360 2.1221 0.0362
## ma1 0.0945 0.1016 0.9301 0.3544
## ma2 -0.6983 0.0796 -8.7726 0.0000cbind(arima312 = AIC(model_arma_1), arima310 = AIC(model_arma_2),
arima011 = AIC(model_arma_3), arima111 = AIC(model_arma_4))## arima312 arima310 arima011 arima111
## [1,] -92.65598 -87.77843 -88.23929 -91.57106cbind(arima312 = BIC(model_arma_1), arima310 = BIC(model_arma_2),
arima011 = BIC(model_arma_3), arima111 = BIC(model_arma_4))## arima312 arima310 arima011 arima111
## [1,] -76.73222 -77.16259 -82.93137 -78.30125Dipilih ARIMA (3,1,0) karena semua ar dan ma nya signifikan
2.2.4 Diagnostik Model
# Plot Sisaan dan Sisaan Terstandarisasi
# Plot Tebaran Sisaan
par(mfrow = c(1,2))
plot(residuals(model_arma_2), type = "p", ylab = 'sisaan'); abline(h=0)
plot(rstandard(model_arma_2), ylab = "sisaan terstandarisasi"); abline(h=0)
t.test(residuals(model_arma_2), mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: residuals(model_arma_2)
## t = 0.93928, df = 105, p-value = 0.3497
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.01542734 0.04319958
## sample estimates:
## mean of x
## 0.01388612Hipotesis yang diuji: \(H_0:\mu=0\) \(H_1 :\mu \ne 0\) Karena p-value = 0.3527 > 0.05 Tidak tolak H0. Tidak cukup bukti bahwa terjadi pelanggaran asumsi rataan tidak 0 (rataan sisaan 0)
# Plot Q-Q dan Histogram
par(mfrow = c(1,2))
h <- hist(residuals(model_arma_2), main = 'Histogram Residual')
xfit <- seq(min(residuals(model_arma_2)), max(residuals(model_arma_2)), length = 100)
yfit <- dnorm (xfit, mean = mean(residuals(model_arma_2)), sd = sd(residuals(model_arma_2)))
yfit <- yfit*diff(h$mids[1:2])*length(residuals(model_arma_2))
lines(xfit, yfit)
qqnorm(residuals(model_arma_2)); qqline(residuals(model_arma_2))
shapiro.test(residuals(model_arma_2))##
## Shapiro-Wilk normality test
##
## data: residuals(model_arma_2)
## W = 0.98934, p-value = 0.5699ks.test(residuals(model_arma_2),'pnorm')##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuals(model_arma_2)
## D = 0.37186, p-value = 3.713e-13
## alternative hypothesis: two-sidedJika berdasarkan Shapiro, maka tidak tolak H0. Artinya, Data tidak menyebar Normal Jika berdasarkan Kolmogorov - Smirnov, maka tolak H0, artinya data sudah menyebar normal.
Hipotesis yang diuji: H0: sisaan menyebar normal H1: sisaan tidak menyebar normal
Tidak tolak H0. Tidak cukup bukti bahwa terjadi pelanggaran asumsi sisaan tidak menyebar normal (sisaan tidak menyebar normal).
# Plot ACF dan PACF Sisaan Model
par(mfrow = c(1,2))
Acf(residuals(model_arma_2), main = "Plot ACF Sisaan")
Pacf(residuals(model_arma_2), main = "Plot PACF Sisaan")
# L-JUng Box Test
Box.test(residuals(model_arma_2), lag = 24, type = "Ljung")##
## Box-Ljung test
##
## data: residuals(model_arma_2)
## X-squared = 14.334, df = 24, p-value = 0.9387Hipotesis yang diuji: H0: tidak ada autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan Tidak tolak H0. TIdak cukup bukti bahwa terjadi pelanggaran asumsi autokorelasi (sisaan saling bebas)
2.2.5 Overfitting Model
# Overfitting
overfit_1 <- Arima(data.ts.2, order = c(3,1,0))
lmtest::coeftest(overfit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.551446 0.090440 -6.0974 1.078e-09 ***
## ar2 -0.318837 0.100648 -3.1679 0.001536 **
## ar3 -0.368504 0.091216 -4.0399 5.347e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1diperoleh ketiganya siginifikan.
# Overfitting
overfit_1 <- Arima(data.ts.2, order = c(4,1,0))
lmtest::coeftest(overfit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.574859 0.098645 -5.8276 5.624e-09 ***
## ar2 -0.339506 0.106518 -3.1873 0.0014360 **
## ar3 -0.399901 0.105490 -3.7909 0.0001501 ***
## ar4 -0.060271 0.102339 -0.5889 0.5559086
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1diperoleh ar4 tidak signifikan
# Underfitting
overfit_1 <- Arima(data.ts.2, order = c(2,1,0))
lmtest::coeftest(overfit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.504044 0.096541 -5.2211 1.779e-07 ***
## ar2 -0.137108 0.096931 -1.4145 0.1572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1diperoleh ar2 tidak signifikan. Maka yang dipakai adalah overfitting dengan order (3,1,0).
# Overfitting
overfit_1 <- Arima(data.ts.2, order = c(3,1,0))
lmtest::coeftest(overfit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.551446 0.090440 -6.0974 1.078e-09 ***
## ar2 -0.318837 0.100648 -3.1679 0.001536 **
## ar3 -0.368504 0.091216 -4.0399 5.347e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(overfit_1)## Series: data.ts.2
## ARIMA(3,1,0)
##
## Coefficients:
## ar1 ar2 ar3
## -0.5514 -0.3188 -0.3685
## s.e. 0.0904 0.1006 0.0912
##
## sigma^2 = 0.02405: log likelihood = 47.89
## AIC=-87.78 AICc=-87.38 BIC=-77.16
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01388612 0.1521235 0.1162104 Inf Inf 0.7862175 -0.03443426printstatarima(overfit_1)##
## Coefficients:
## s.e. t sign.
## ar1 -0.5514 0.0904 -6.0996 0e+00
## ar2 -0.3188 0.1006 -3.1690 2e-03
## ar3 -0.3685 0.0912 -4.0406 1e-04c(AR1 = AIC(model_arma_2), Ar2 = AIC(overfit_1))## AR1 Ar2
## -87.77843 -87.778432.2.6 Peramalan
ramal_arima <- predict(model_arma_2, n.ahead = 20)
fitted_pred = ramal_arima$pred
fitted_ci_lo = ramal_arima$pred - 1.96*ramal_arima$se
fitted_ci_up = ramal_arima$pred + 1.96*ramal_arima$se
ramalan_arima = data.frame(prediksi = fitted_pred, Lower = fitted_ci_lo, upper = fitted_ci_up)
ramalan_arima## prediksi Lower upper
## 1 0.1049501 -0.1990021 0.4089023
## 2 0.2264661 -0.1066633 0.5595956
## 3 0.1780019 -0.1802740 0.5362777
## 4 0.1788995 -0.1872992 0.5450983
## 5 0.1490775 -0.2546110 0.5527661
## 6 0.1830958 -0.2433738 0.6095655
## 7 0.1735141 -0.2772644 0.6242927
## 8 0.1789412 -0.2867666 0.6446489
## 9 0.1664675 -0.3204697 0.6534048
## 10 0.1751466 -0.3300040 0.6802973
## 11 0.1723378 -0.3522240 0.6968996
## 12 0.1757161 -0.3648893 0.7163214
## 13 0.1715504 -0.3862061 0.7293069
## 14 0.1738055 -0.3999356 0.7475466
## 15 0.1726452 -0.4174254 0.7627157
## 16 0.1741011 -0.4310986 0.7793008
## 17 0.1728372 -0.4475176 0.7931920
## 18 0.1734975 -0.4614026 0.8083977
## 19 0.1729999 -0.4763937 0.8223934
## 20 0.1735295 -0.4898322 0.8368913# menggabungkan data asli dan ramalan
all_result = data.frame(Data = c(data.ts.2, rep(NA, 20)), Ramalan = c(rep(NA, 106), fitted_pred),
Lower = c(rep(NA, 106), fitted_ci_lo), Upper = c(rep(NA, 106), fitted_ci_up))
plot(all_result$Data, type = 'l', ylim = c(-0.8, 0.8),
ylab = expression(Y[t]), xlab = 'Periode')
polygon(c(1:dim(all_result)[1], rev(1:dim(all_result)[1])),
c(all_result$Lower, rev(all_result$Upper)), col = "grey75", border = F)
lines(all_result$Ramalan)
pred_1<- forecast(model_arma_2, h = 20)
pred_1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 107 0.1049501 -0.093789945 0.3036901 -0.1989966 0.4088967
## 108 0.2264661 0.008648492 0.4442838 -0.1066572 0.5595895
## 109 0.1780019 -0.056257837 0.4122616 -0.1802674 0.5362711
## 110 0.1788995 -0.060540591 0.4183396 -0.1872925 0.5450915
## 111 0.1490775 -0.114875375 0.4130304 -0.2546036 0.5527586
## 112 0.1830958 -0.095752544 0.4619442 -0.2433660 0.6095577
## 113 0.1735141 -0.121228706 0.4682570 -0.2772561 0.6242844
## 114 0.1789412 -0.125563197 0.4834455 -0.2867581 0.6446404
## 115 0.1664675 -0.151917730 0.4848528 -0.3204607 0.6533958
## 116 0.1751466 -0.155147552 0.5054408 -0.3299947 0.6802880
## 117 0.1723378 -0.170648466 0.5153240 -0.3522144 0.6968899
## 118 0.1757161 -0.177760274 0.5291924 -0.3648793 0.7163115
## 119 0.1715504 -0.193140292 0.5362411 -0.3861959 0.7292967
## 120 0.1738055 -0.201336776 0.5489478 -0.3999251 0.7475361
## 121 0.1726452 -0.213174127 0.5584645 -0.4174145 0.7627049
## 122 0.1741011 -0.221610440 0.5698126 -0.4310875 0.7792897
## 123 0.1728372 -0.232783556 0.5784579 -0.4475062 0.7931806
## 124 0.1734975 -0.241633709 0.5886288 -0.4613909 0.8083860
## 125 0.1729999 -0.251607948 0.5976077 -0.4763817 0.8223815
## 126 0.1735295 -0.260211448 0.6072705 -0.4898200 0.83687912.3 Deret Tidak Stasioner (Ragam)
Berikut merupakan Data Deret Waktu Harga Minyak Mentah.
data("oil.price")
head(oil.price, 15)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1986 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11
## 1987 18.65 17.75 18.30plot(oil.price, ylab = 'Yt', main = "Harga Minyak Mentah")
adf.test(oil.price)## Warning in adf.test(oil.price): p-value greater than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: oil.price
## Dickey-Fuller = 0.7045, Lag order = 6, p-value = 0.99
## alternative hypothesis: stationarytrans <- BoxCox.ar(oil.price)## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
trans$mle## [1] 0.1trans$ci## [1] -0.1 0.3log.oil = log(oil.price)plot(log.oil, main = 'Transformasi Logaritma Harga Minyak Mentah', ylab = 'Yt')
adf.test(log.oil)##
## Augmented Dickey-Fuller Test
##
## data: log.oil
## Dickey-Fuller = -1.1119, Lag order = 6, p-value = 0.9189
## alternative hypothesis: stationaryplot(diff(log.oil), main = 'Transformasi DIfferencing Hasil Logaritma Harga Minyak Mentah', ylab = 'Yt')
plot(diff(oil.price), main = 'Transformasi DIfferencing Hasil Logaritma Harga Minyak Mentah', ylab = 'Yt')
Bandingkan dengan Hasil Plot jika data tidak di-logaritma kan terlebih dahulu. Terlihat pada ujung plot, ragamnnya tidak homogen (terjadi lonjakan yang ragamnya berbeda dengan yang sebelumnya).
Akibatnya, yang harus ditangani terlebih dahulu adalah: 1. Data yang tidak stasioner dalam ragam 2. Data tersebut ditransformasikan agar stasioner dalam ragam (ragamnya homogen) 3. Kemudian, dilanjutkan dengan melakukan differencing pada data, agar data menjadi stasioner data niai tengah
acf(diff(log.oil, main = "ACF Plot"))
par(mfrow = c(1,2))
Acf(diff(log.oil), main = 'ACF Plot')
Pacf(diff(log.oil), main = 'PACF Plot')
eacf(diff(log.oil))## 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 x o o o
## 2 o x o o o o o o o o o o o o
## 3 o x o o o o o o o o o o o o
## 4 o x x o o o o o o o o o o o
## 5 o x o x o o o o o o o o o o
## 6 o x o x o o o o o o o o o o
## 7 x x o x o o o o o o o o o oauto.arima(log.oil)## Series: log.oil
## ARIMA(1,1,0)(1,0,0)[12] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.2337 -0.0055 0.0040
## s.e. 0.0661 0.0726 0.0069
##
## sigma^2 = 0.006864: log likelihood = 258.72
## AIC=-509.44 AICc=-509.27 BIC=-495.52mod_fit_1 <- Arima(log.oil, order = c(1,1,0), seasonal = c(1,0,0))
summary(mod_fit_1)## Series: log.oil
## ARIMA(1,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 sar1
## 0.2364 -0.0038
## s.e. 0.0660 0.0726
##
## sigma^2 = 0.006844: log likelihood = 258.55
## AIC=-511.11 AICc=-511.01 BIC=-500.67
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003486933 0.08221294 0.0624257 0.06998619 2.036076 0.2763038
## ACF1
## Training set 0.008493756lmtest::coeftest(mod_fit_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.2363981 0.0659736 3.5832 0.0003394 ***
## sar1 -0.0037812 0.0725971 -0.0521 0.9584611
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1printstatarima(mod_fit_1)##
## Coefficients:
## s.e. t sign.
## ar1 0.2364 0.0660 3.5818 0.0004
## sar1 -0.0038 0.0726 -0.0523 0.9583mod_fit_2 <- Arima(log.oil, order = c(1,1,2))
summary(mod_fit_2)## Series: log.oil
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.8429 -0.5730 -0.3104
## s.e. 0.1548 0.1595 0.0675
##
## sigma^2 = 0.006682: log likelihood = 261.88
## AIC=-515.75 AICc=-515.58 BIC=-501.83
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005441889 0.08106116 0.06170506 0.108913 2.010809 0.2731142
## ACF1
## Training set -0.03202714lmtest::coeftest(mod_fit_2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.842925 0.154810 5.4449 5.183e-08 ***
## ma1 -0.572957 0.159451 -3.5933 0.0003265 ***
## ma2 -0.310443 0.067461 -4.6018 4.188e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_fit_3 <- Arima(log.oil, order = c(0,1,1))
summary(mod_fit_3)## Series: log.oil
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.2956
## s.e. 0.0693
##
## sigma^2 = 0.006717: log likelihood = 260.29
## AIC=-516.58 AICc=-516.53 BIC=-509.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.00349408 0.08161418 0.06196464 0.07021944 2.020688 0.2742631
## ACF1
## Training set -0.04384663lmtest::coeftest(mod_fit_3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.295600 0.069347 4.2626 2.02e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_fit_4 <- Arima(log.oil, order = c(0,1,2))
summary(mod_fit_4)## Series: log.oil
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## 0.269 -0.0941
## s.e. 0.068 0.0746
##
## sigma^2 = 0.0067: log likelihood = 261.08
## AIC=-516.17 AICc=-516.06 BIC=-505.72
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003774913 0.08134217 0.06171914 0.07383129 2.012101 0.2731765
## ACF1
## Training set -0.01820862lmtest::coeftest(mod_fit_4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.268962 0.068035 3.9533 7.708e-05 ***
## ma2 -0.094137 0.074555 -1.2627 0.2067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cbind(Seasonal= AIC(mod_fit_1), arima112 = AIC(mod_fit_2),
arima011 = AIC(mod_fit_3), arima012 = AIC(mod_fit_4))## Seasonal arima112 arima011 arima012
## [1,] -511.1078 -515.7545 -516.5827 -516.1665cbind(Seasonal = BIC(mod_fit_1), arima112 = BIC(mod_fit_2),
arima011 = BIC(mod_fit_3), arima012 = BIC(mod_fit_4))## Seasonal arima112 arima011 arima012
## [1,] -500.6659 -501.832 -509.6214 -505.7246# Plot Sisaan dan Sisaan Terstandarisasi
# Plot Tebaran Sisaan
par(mfrow = c(1,2))
plot(residuals(mod_fit_3), type = "p", ylab = 'sisaan'); abline(h=0)
plot(rstandard(mod_fit_3), ylab = "sisaan terstandarisasi"); abline(h=0)
t.test(residuals(mod_fit_3), mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: residuals(mod_fit_3)
## t = 0.66385, df = 240, p-value = 0.5074
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.006874165 0.013862326
## sample estimates:
## mean of x
## 0.00349408Hipotesis yang diuji: \(H_0:\mu=0\) \(H_1 :\mu \ne 0\) Karena p-value = 0.3527 > 0.05 Tidak tolak H0. Tidak cukup bukti bahwa terjadi pelanggaran asumsi rataan tidak 0 (rataan sisaan 0)
# Plot Q-Q dan Histogram
par(mfrow = c(1,2))
h <- hist(residuals(mod_fit_3), main = 'Histogram Residual')
xfit <- seq(min(residuals(mod_fit_3)), max(residuals(mod_fit_3)), length = 100)
yfit <- dnorm (xfit, mean = mean(residuals(mod_fit_3)), sd = sd(residuals(mod_fit_3)))
yfit <- yfit*diff(h$mids[1:2])*length(residuals(mod_fit_3))
lines(xfit, yfit)
qqnorm(residuals(mod_fit_3)); qqline(residuals(mod_fit_3))
shapiro.test(residuals(mod_fit_3))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_fit_3)
## W = 0.96883, p-value = 3.937e-05ks.test(residuals(mod_fit_3),'pnorm')##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuals(mod_fit_3)
## D = 0.42369, p-value < 2.2e-16
## alternative hypothesis: two-sidedJika berdasarkan Shapiro, maka tidak tolak H0. Artinya, Data menyebar tidak Normal Jika berdasarkan Kolmogorov - Smirnov, maka tidak tolak H0, artinya data menyebar tidak normal.
Hipotesis yang diuji: H0: sisaan menyebar normal H1: sisaan tidak menyebar normal
Karena p-value < 0.05 maka tidak cukup bukti untuk tolak H0
# Plot ACF dan PACF Sisaan Model
par(mfrow = c(1,2))
Acf(residuals(mod_fit_3), main = "Plot ACF Sisaan")
Pacf(residuals(mod_fit_3), main = "Plot PACF Sisaan")
# L-JUng Box Test
Box.test(residuals(mod_fit_3), lag = 36, type = "Ljung")##
## Box-Ljung test
##
## data: residuals(mod_fit_3)
## X-squared = 45.311, df = 36, p-value = 0.1374Hipotesis yang diuji: H0: tidak ada autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan
2.3.1 Overfitting Model
# Overfitting
overfit_3 <- Arima(log.oil, order = c(0,1,2))
lmtest::coeftest(overfit_3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.268962 0.068035 3.9533 7.708e-05 ***
## ma2 -0.094137 0.074555 -1.2627 0.2067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1diperoleh ma2 tidak siginifikan.
2.3.2 Peramalan
prediksi_3 <- forecast(mod_fit_3, h = 20)
plot(prediksi_3)
2.4 Deret Musiman
data_seasonal <- read.csv('5. seasonal.csv')
head(data_seasonal, 15)## Waktu Xt
## 1 1 24.083
## 2 2 26.146
## 3 3 27.341
## 4 4 22.431
## 5 5 29.361
## 6 6 28.511
## 7 7 30.432
## 8 8 26.929
## 9 9 31.235
## 10 10 32.472
## 11 11 35.815
## 12 12 30.104
## 13 13 35.892
## 14 14 37.232
## 15 15 37.512ts_seasonal <-ts(data_seasonal$Xt)
plot(ts_seasonal, main = 'Plot Data', ylab = 'Yt')
ts.plot(ts_seasonal[1:40], main = 'Plot Data', ylab = 'Yt')
points(ts_seasonal[1:40])
ts_season_2 <- ts(data_seasonal$Xt, frequency = 4)Acf(ts_season_2, main = "ACF Data Training", lag.max = 36)
adf.test(ts_season_2)##
## Augmented Dickey-Fuller Test
##
## data: ts_season_2
## Dickey-Fuller = -2.1969, Lag order = 5, p-value = 0.4944
## alternative hypothesis: stationaryTidak Stasioner
diff_season_1 <- diff(ts_season_2)ts.plot(diff_season_1, main = "Plot Hasil Pembedaan Reguler", ylab = 'Yt')
adf.test(diff_season_1)## Warning in adf.test(diff_season_1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff_season_1
## Dickey-Fuller = -4.8533, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarySudah stasioner dalam rataan.
acf_1 <- acf(diff_season_1, lag.max = 35, xaxt = "n", main = "ACF d1 ")
axis(1, at = 0:36/4, labels= 0:36)
tidak stasioner untuk data musimannya.
acf_1$lag <- acf_1$lag * 4
acf_1_1 <- as.data.frame(cbind(acf_1$acf, acf_1$lag))
acf_1_2 <- acf_1_1[which(acf_1_1$V2%%4 == 0),]
barplot(height = acf_1_2$V1, names.arg = acf_1_2$V2, ylab = "ACF", xlab = "Lag", main = "ACF d1")
diff_seasonal <- diff(ts_season_2, lag = 4, differences = 1) # melakukan differencing untuk musiman lag = 4
ts.plot(diff_seasonal, type = 'l', main = "Plot Data Training Xt-d1", xlab = "Waktu(t)", ylab = "Xt d1")
acf2 <- Acf(diff_seasonal, lag_max = 36, xaxt = 'n', main = 'ACF D1')## Warning in plot.window(...): "lag_max" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "lag_max" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "lag_max" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "lag_max" is not a
## graphical parameter## Warning in box(...): "lag_max" is not a graphical parameter## Warning in title(...): "lag_max" is not a graphical parameter
acf2$lag <- acf2$lag *4
acf2.1 <- as.data.frame(cbind(acf2$acf, acf2$lag))
acf2.2 <- acf2.1[which(acf2.1$V2%%4==0),]
barplot(height = acf2.2$V1, names.arg = acf2.2$V2, ylab = "ACF", xlab = "Lag")
adf.test(diff_seasonal)##
## Augmented Dickey-Fuller Test
##
## data: diff_seasonal
## Dickey-Fuller = -3.7262, Lag order = 5, p-value = 0.0246
## alternative hypothesis: stationaryseason_d1D1 <- diff(diff_seasonal) # data hasil differencing seasonal dan reguler
ts.plot(season_d1D1, main = "Plot Data Training Xt-d1D1", xlab = "Waktu(t)", ylab = "Xt d1D1")
Data telah stasioner baik di reguler maupun di seasonal
par(mfrow = c(1,2))
acf3 <- acf(season_d1D1, lag.max = 24, xaxt = "n", main = "ACF d1D1 - Non Musiman")
axis(1, at = 0:24/4, labels = 0:24)
pacf3 <- pacf(season_d1D1, lag.max = 24, xaxt = "n", main = "PACF d1D1 - Non Musiman")
axis(1, at = 0:24/4, labels = 0:24)
Plot ACF cut-off setelah lag 1 (reguler) dan lag 4 (musiman), sementara plot PACF tail-off. Sehingga kandidat modelnya adalah \(ARIMA(0,1,1)(0,1,1)^4\)
auto.arima(ts_season_2)## Series: ts_season_2
## ARIMA(1,1,1)(1,1,1)[4]
##
## Coefficients:
## ar1 ma1 sar1 sma1
## -0.2773 -0.5917 -0.3390 -0.4511
## s.e. 0.1097 0.0942 0.1216 0.1173
##
## sigma^2 = 0.8663: log likelihood = -195.08
## AIC=400.16 AICc=400.59 BIC=415.04model_seasonal <- Arima(ts_season_2, order = c(0,1,1), seasonal = c(0,1,1))
summary(model_seasonal)## Series: ts_season_2
## ARIMA(0,1,1)(0,1,1)[4]
##
## Coefficients:
## ma1 sma1
## -0.7544 -0.6292
## s.e. 0.0517 0.0610
##
## sigma^2 = 0.9201: log likelihood = -200.36
## AIC=406.71 AICc=406.88 BIC=415.64
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01824931 0.9365558 0.7369625 0.01625189 0.9705009 0.1717576
## ACF1
## Training set -0.1015128lmtest::coeftest(model_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.754432 0.051651 -14.606 < 2.2e-16 ***
## sma1 -0.629226 0.060991 -10.317 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Plot Q-Q dan Histogram
par(mfrow = c(1,2))
h <- hist(residuals(model_seasonal), main = 'Histogram Residual')
xfit <- seq(min(residuals(model_seasonal)), max(residuals(model_seasonal)), length = 100)
yfit <- dnorm (xfit, mean = mean(residuals(model_seasonal)), sd = sd(residuals(model_seasonal)))
yfit <- yfit*diff(h$mids[1:2])*length(residuals(model_seasonal))
lines(xfit, yfit)
qqnorm(residuals(model_seasonal)); qqline(residuals(model_seasonal))
shapiro.test(residuals(model_seasonal))##
## Shapiro-Wilk normality test
##
## data: residuals(model_seasonal)
## W = 0.99393, p-value = 0.7847ks.test(residuals(model_seasonal),'pnorm')##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuals(model_seasonal)
## D = 0.049138, p-value = 0.8619
## alternative hypothesis: two-sidedHipotesis yang diuji: H0: sisaan menyebar normal H1: sisaan tidak menyebar normal
par(mfrow = c(1,2))
plot(residuals(model_seasonal), type = "p", ylab = 'sisaan'); abline(h = 0)
plot(rstandard(model_seasonal), ylab = 'sisaan terstandarisasi'); abline(h=0)
par(mfrow = c(1,2))
Acf(residuals(model_seasonal), main = "Plot ACF Sisaan")
Pacf(residuals(model_seasonal), main = "Plot PACF Sisaan")
# L-JUng Box Test
Box.test(residuals(model_seasonal), lag = 36, type = "Ljung")##
## Box-Ljung test
##
## data: residuals(model_seasonal)
## X-squared = 26.847, df = 36, p-value = 0.8657Hipotesis yang diuji: H0: tidak ada autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan
prediksi_season <- forecast(model_seasonal, h = 40)
plot(prediksi_season)
3 ARCH / GARCH
3.1 Data JISDOR
library(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(TTR)
library(fpp2)## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──## ✔ ggplot2 3.4.3 ✔ expsmooth 2.3
## ✔ fma 2.5## library(tseries)
library(astsa)##
## Attaching package: 'astsa'## The following objects are masked from 'package:fma':
##
## chicken, sales## The following object is masked from 'package:fpp2':
##
## oil## The following object is masked from 'package:forecast':
##
## gaslibrary(fpp3)## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.2.1 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.0 ✔ feasts 0.3.1
## ✔ lubridate 1.9.2 ✔ fable 0.3.3
## ✔ tsibble 1.1.3 ✔ fabletools 0.3.4## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()##
## Attaching package: 'fpp3'## The following object is masked from 'package:fpp2':
##
## insurancelibrary(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:forecast':
##
## gghistogramlibrary(fGarch)## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")##
## Attaching package: 'fGarch'## The following object is masked from 'package:TTR':
##
## volatilitylibrary(rugarch)## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:fabletools':
##
## report## The following object is masked from 'package:stats':
##
## sigmalibrary(TSA)
library(aTSA)##
## Attaching package: 'aTSA'## The following objects are masked from 'package:fabletools':
##
## estimate, forecast## The following objects are masked from 'package:tseries':
##
## adf.test, kpss.test, pp.test## The following object is masked from 'package:forecast':
##
## forecast## The following object is masked from 'package:graphics':
##
## identifylibrary(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:tsibble':
##
## index## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:fpp3':
##
## prices## The following objects are masked from 'package:TSA':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legendlibrary(xts)Data JISDOR, yaitu kurs referensi yang merepresentasikan nilai tukar spot dolar AS terhadap rupiah yang dihitung dari seluruh transaksi spot USD/IDR antar bank di pasar valuta asing domestik yang dilakukan oleh seluruh Bank Devisa yang tercatat di Sitem monitoring Transaksi Valuta Asing terhadap rupiah (SISMONTAVAR). Data dari periode 4 Januari 2016 hingga 29 Oktober 2021.
df <- read.csv('5. garch.csv')%>%
mutate(Tanggal = as.Date(Tanggal, format = '%d/%m/%Y'))%>%
as_tsibble(index = Tanggal, regular = F)
head(df, 15)## # A tsibble: 15 x 4 [!]
## NO Tanggal Kurs Return
## <int> <date> <int> <dbl>
## 1 1 2016-01-04 13898 NA
## 2 2 2016-01-05 13931 0.00237
## 3 3 2016-01-06 13863 -0.00489
## 4 4 2016-01-07 13946 0.00597
## 5 5 2016-01-08 13874 -0.00518
## 6 6 2016-01-11 13935 0.00439
## 7 7 2016-01-12 13835 -0.00720
## 8 8 2016-01-13 13861 0.00188
## 9 9 2016-01-14 13877 0.00115
## 10 10 2016-01-15 13886 0.000648
## 11 11 2016-01-18 13931 0.00324
## 12 12 2016-01-19 13921 -0.000718
## 13 13 2016-01-20 13896 -0.00180
## 14 14 2016-01-21 13899 0.000216
## 15 15 2016-01-22 13874 -0.00180head(diff(log(df$Kurs))) # nilai return adalah hasil kurs yang diff dan log-kan sekali## [1] 0.002371628 -0.004893152 0.005969308 -0.005176144 0.004387076
## [6] -0.007202048df %>%
autoplot(Return)+
labs(title = "Plot Data Return - Transaksi Spot USD/IDR (JISDOR)",
subtitle = "04 Januari 2016 - 29 Oktober 2021") +
ylab("Return USD/IDR")+
theme_bw()## Warning: Removed 1 row containing missing values (`geom_line()`).
df_frame = as.data.frame(df)
ggplot(df_frame, aes(Return)) +
geom_histogram(aes(y = ..density..), color = 'steelblue', fill = "grey")+
stat_function(fun = dnorm,
args = list(mean = mean(df_frame$Return, na.rm = T),
sd = sd(df_frame$Return, na.rm = T)))+
theme_bw()## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
daily_ret_xts <- xts(df_frame[,4], order.by = df_frame[,2])
realizedvol <- rollapply(daily_ret_xts, width = 20, FUN= sd.annualized)vol <- data.frame(index(realizedvol),realizedvol)
colnames(vol) <- c("date", "volatility")p3 <- ggplot(vol, aes(x = date, y = volatility))
p3 + geom_line(color = "steelblue") + theme_bw()## Warning: Removed 19 rows containing missing values (`geom_line()`).
acf.raw <-
ACF(df, Return, lag_max = 20) %>%
autoplot() +
labs(title = "Auto Correlation Function",
subtitle = "Data Return Harian") +
ylab("ACF")+
theme_bw()## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.pacf.raw <-
PACF(df, Return, lag_max = 20) %>%
autoplot() +
labs(title = "Partial Auto Correlation Function",
subtitle = "Data Return Harian") +
ylab("PACF")## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.theme_bw()## List of 97
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.75points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.75points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.75points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.2points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'simpleUnit' num 1.2lines
## ..- attr(*, "unit")= int 3
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "unit")= int 1
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'simpleUnit' num 5.5points
## ..- attr(*, "unit")= int 8
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ linewidth : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.clip : chr "inherit"
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.x.bottom : NULL
## $ strip.text.x.top : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.y.right : NULL
## $ strip.switch.pad.grid : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.switch.pad.wrap : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUEggarrange(acf.raw, pacf.raw, ncol = 2)
tseries::adf.test(df$Return%>%na.omit())## Warning in tseries::adf.test(df$Return %>% na.omit()): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: df$Return %>% na.omit()
## Dickey-Fuller = -9.9219, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary3.2 Pemodelan Nilai Tengah (ARIMA)
# splitting
df.train <- slice_head(df, prop = 0.80)
df.test <- df[-index(df.train),]
tail(df.train)## # A tsibble: 6 x 4 [!]
## NO Tanggal Kurs Return
## <int> <date> <int> <dbl>
## 1 1126 2020-08-25 14632 -0.0110
## 2 1127 2020-08-26 14636 0.000273
## 3 1128 2020-08-27 14714 0.00532
## 4 1129 2020-08-28 14702 -0.000816
## 5 1130 2020-08-31 14554 -0.0101
## 6 1131 2020-09-01 14615 0.00418head(df.test)## # A tsibble: 6 x 4 [!]
## NO Tanggal Kurs Return
## <int> <date> <int> <dbl>
## 1 1132 2020-09-02 14804 0.0128
## 2 1133 2020-09-03 14818 0.000945
## 3 1134 2020-09-04 14792 -0.00176
## 4 1135 2020-09-07 14754 -0.00257
## 5 1136 2020-09-08 14798 0.00298
## 6 1137 2020-09-09 14853 0.00371dim(df.train)## [1] 1131 4dim(df.test)## [1] 283 4df <- df[-1,]auto.arima(df.train$Return[-1], trace = F, method = 'ML')## Series: df.train$Return[-1]
## ARIMA(1,0,1) with zero mean
##
## Coefficients:
## ar1 ma1
## 0.7803 -0.6774
## s.e. 0.0972 0.1151
##
## sigma^2 = 1.922e-05: log likelihood = 4533.34
## AIC=-9060.68 AICc=-9060.65 BIC=-9045.59TSA::eacf(df.train$Return[-1])## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x o o o o x o o x o x
## 1 x o o o o o o o x o o x o o
## 2 x x o o o o o o x o o o o o
## 3 x x x o o o o o x o o o o o
## 4 x x x o o o o o x o o o o x
## 5 x x o o o o o o o o o o o x
## 6 x x o o o o o o o o o x o o
## 7 x x o o x o x o o o o o o omod_fit1 <- stats::arima(df.train$Return[-1], order = c(1,0,1),
include.mean = F, method = 'ML')
lmtest::coeftest(mod_fit1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.780338 0.097178 8.030 9.749e-16 ***
## ma1 -0.677439 0.115055 -5.888 3.910e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1checkresiduals(mod_fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with zero mean
## Q* = 16.636, df = 8, p-value = 0.03413
##
## Model df: 2. Total lags used: 10Ragam dari sisaan tidak homogen
3.3 Identifikasi Efek ARCH:
LM Test -> mengecek apakah ada volatilitas dalam sisaan?
res_arma <- ts(residuals(mod_fit1))
plot.ts(res_arma, ylab = 'residual')
plot.ts(res_arma^2, ylab = "Residual Kuadrat")
par(mfrow = c(1,2))
Acf(res_arma^2)
Pacf(res_arma^2)
Masih ada sisaan yg signifikan -> sisaan tidak saling bebas
aTSA::arch.test(mod_fit1)## ARCH heteroscedasticity test for residuals
## alternative: heteroscedastic
##
## Portmanteau-Q test:
## order PQ p.value
## [1,] 4 351 0
## [2,] 8 539 0
## [3,] 12 722 0
## [4,] 16 950 0
## [5,] 20 1025 0
## [6,] 24 1046 0
## Lagrange-Multiplier test:
## order LM p.value
## [1,] 4 711.7 0.00e+00
## [2,] 8 306.5 0.00e+00
## [3,] 12 172.6 0.00e+00
## [4,] 16 113.1 0.00e+00
## [5,] 20 85.6 1.95e-10
## [6,] 24 66.5 4.19e-06
Dilihat dari nilai Lagrange-Multiplier Test (LM Test), nilai p-value uji LM tersebut hingga lag ke-24 adalah kurang dari 0,05. Sehingga dapat disimpulkan bahwa sisaan terdapat heteroskedastisitas yang belum bisa diatasi oleh model ARIMA yang digunakan. Heteroscedastisitas tersebut dapat diatasi dengan ARCH/GARCH. Nilai p-value signifikan hingga lag ke-24, hal ini berarti ada indikasi proses long memory. Adanya indikasi proses long memory ini membuat penggunaan model ARCH menjadi kurang tepat. Sehingga pemodelan untuk ragam bersyarat yang dilakukan adalah pemodelan GARCH.
3.4 Pemilihan Model GARCH
TSA::eacf(res_arma^2)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 x x o o x o x o x o x o o x
## 2 x o x o o o x o x o x o o x
## 3 x x x o o o o o o o x o o o
## 4 x x x o x o o o x o x x o x
## 5 x x x x o o o o o o o o o o
## 6 x o x x x x o o o o o o o o
## 7 x x x x x x o o o o o o o oauto.arima(res_arma^2)## Series: res_arma^2
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.7104 -0.1699
## s.e. 0.0287 0.0293
##
## sigma^2 = 2.563e-09: log likelihood = 9565.31
## AIC=-19124.62 AICc=-19124.6 BIC=-19109.53#ARIMA(1,0,1) + GARCH(1,1)
garchSpec11 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(1,1)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE), #ARIMA (1,0,1)
distribution.model = 'norm')
(garchFit11 <- rugarch::ugarchfit(spec = garchSpec11,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.21241 0.445893 -0.47637 0.63381
## ma1 0.26823 0.437611 0.61293 0.53992
## omega 0.00000 0.000001 0.45671 0.64788
## alpha1 0.16945 0.044504 3.80753 0.00014
## beta1 0.82365 0.036098 22.81684 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.21241 1.403368 -0.151358 0.87969
## ma1 0.26823 1.507873 0.177884 0.85881
## omega 0.00000 0.000023 0.016392 0.98692
## alpha1 0.16945 0.957691 0.176938 0.85956
## beta1 0.82365 0.831176 0.990940 0.32171
##
## LogLikelihood : 4770.09
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.4338
## Bayes -8.4115
## Shibata -8.4338
## Hannan-Quinn -8.4254
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.150 0.1425623
## Lag[2*(p+q)+(p+q)-1][5] 5.511 0.0005274
## Lag[4*(p+q)+(p+q)-1][9] 7.869 0.0643428
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.787 0.1813
## Lag[2*(p+q)+(p+q)-1][5] 3.197 0.3723
## Lag[4*(p+q)+(p+q)-1][9] 4.189 0.5572
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.715 0.500 2.000 0.1904
## ARCH Lag[5] 1.777 1.440 1.667 0.5227
## ARCH Lag[7] 2.577 2.315 1.543 0.5966
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 92.9581
## Individual Statistics:
## ar1 0.04553
## ma1 0.03817
## omega 28.35826
## alpha1 0.10301
## beta1 0.17498
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.30156 0.7630
## Negative Sign Bias 0.06786 0.9459
## Positive Sign Bias 1.17418 0.2406
## Joint Effect 1.48974 0.6846
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 70.74 6.917e-08
## 2 30 80.30 1.047e-06
## 3 40 111.20 7.312e-09
## 4 50 121.59 4.195e-08
##
##
## Elapsed time : 0.2213001Hipotesis yang diuji pada Ljung-Box: H0: tidak ada autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan Karena p-value > 0.05 maka Tidak tolak H0. TIdak cukup bukti bahwa terjadi pelanggaran asumsi autokorelasi (sisaan saling bebas)
Hipotesis yang diuji pada Adjusted Pearson Goodness-of-Fit Test: H0: menyebar normal H1: tidak menyebar normal Tidak tolak H0. p-value < 0.05 artinya tolak H0 -> data sisaan tidak menyebar normal
#ARIMA(1,0,1) + GARCH(1,1)
garchSpec11 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(1,1)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE), #ARIMA (1,0,1)
distribution.model = 'sstd')
(garchFit11 <- rugarch::ugarchfit(spec = garchSpec11,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.28337 0.462260 0.61302 0.539865
## ma1 -0.31258 0.455683 -0.68596 0.492738
## omega 0.00000 0.000002 0.14542 0.884376
## alpha1 0.19875 0.112883 1.76066 0.078295
## beta1 0.80025 0.089348 8.95657 0.000000
## skew 0.94934 0.035124 27.02863 0.000000
## shape 4.91138 0.845903 5.80608 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.28337 1.262012 0.224541 0.822336
## ma1 -0.31258 1.604692 -0.194792 0.845556
## omega 0.00000 0.000077 0.004166 0.996676
## alpha1 0.19875 4.115027 0.048298 0.961479
## beta1 0.80025 3.240768 0.246932 0.804961
## skew 0.94934 0.197179 4.814630 0.000001
## shape 4.91138 36.907326 0.133073 0.894135
##
## LogLikelihood : 4825.992
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.5292
## Bayes -8.4980
## Shibata -8.5293
## Hannan-Quinn -8.5174
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 14.57 1.347e-04
## Lag[2*(p+q)+(p+q)-1][5] 18.12 0.000e+00
## Lag[4*(p+q)+(p+q)-1][9] 20.50 3.871e-08
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.392 0.2381
## Lag[2*(p+q)+(p+q)-1][5] 3.400 0.3389
## Lag[4*(p+q)+(p+q)-1][9] 4.746 0.4676
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.924 0.500 2.000 0.1654
## ARCH Lag[5] 2.162 1.440 1.667 0.4365
## ARCH Lag[7] 3.211 2.315 1.543 0.4746
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 177.3918
## Individual Statistics:
## ar1 0.7288
## ma1 0.7505
## omega 62.9521
## alpha1 0.1479
## beta1 0.1757
## skew 0.1977
## shape 0.2926
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2569 0.7973
## Negative Sign Bias 0.3447 0.7304
## Positive Sign Bias 0.6230 0.5334
## Joint Effect 0.5676 0.9038
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 16.51 0.6228
## 2 30 23.06 0.7737
## 3 40 30.78 0.8235
## 4 50 48.41 0.4971
##
##
## Elapsed time : 0.442096#ARIMA(1,0,1) + GARCH(2,1)
garchSpec21 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(2,1)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE)) #ARIMA (1,0,1)
(garchFit21 <- rugarch::ugarchfit(spec = garchSpec21,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(2,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.199487 0.429860 -0.464075 0.642594
## ma1 0.255809 0.415767 0.615270 0.538377
## omega 0.000000 0.000001 0.300606 0.763715
## alpha1 0.168527 0.034255 4.919812 0.000001
## alpha2 0.000082 0.065135 0.001259 0.998995
## beta1 0.824345 0.066384 12.417740 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.199487 5.538248 -0.036020 0.97127
## ma1 0.255809 6.278054 0.040747 0.96750
## omega 0.000000 0.000051 0.007372 0.99412
## alpha1 0.168527 0.436075 0.386463 0.69915
## alpha2 0.000082 3.211586 0.000026 0.99998
## beta1 0.824345 2.899853 0.284271 0.77620
##
## LogLikelihood : 4770.137
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.4321
## Bayes -8.4054
## Shibata -8.4322
## Hannan-Quinn -8.4220
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.151 0.1424569
## Lag[2*(p+q)+(p+q)-1][5] 5.462 0.0006326
## Lag[4*(p+q)+(p+q)-1][9] 7.812 0.0674596
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.808 0.1788
## Lag[2*(p+q)+(p+q)-1][8] 4.011 0.5045
## Lag[4*(p+q)+(p+q)-1][14] 4.872 0.7852
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.061 0.500 2.000 0.8049
## ARCH Lag[6] 0.609 1.461 1.711 0.8611
## ARCH Lag[8] 1.161 2.368 1.583 0.8993
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 129.0315
## Individual Statistics:
## ar1 0.04742
## ma1 0.03904
## omega 28.48040
## alpha1 0.10804
## alpha2 0.57484
## beta1 0.18298
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.26980 0.7874
## Negative Sign Bias 0.04336 0.9654
## Positive Sign Bias 1.17407 0.2406
## Joint Effect 1.50560 0.6810
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 68.94 1.380e-07
## 2 30 80.25 1.066e-06
## 3 40 106.67 3.329e-08
## 4 50 123.36 2.411e-08
##
##
## Elapsed time : 0.1676159#ARIMA(1,0,1) + GARCH(1,2)
garchSpec12 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(1,2)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE), #ARIMA (1,0,1)
distribution.model = 'norm')
(garchFit12 <- rugarch::ugarchfit(spec = garchSpec12,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,2)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.16262 0.553586 -0.29375 0.768946
## ma1 0.20668 0.546833 0.37796 0.705458
## omega 0.00000 0.000001 0.49591 0.619959
## alpha1 0.22862 0.060198 3.79786 0.000146
## beta1 0.38053 0.144189 2.63912 0.008312
## beta2 0.38527 0.118034 3.26403 0.001098
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.16262 0.975835 -0.166645 0.86765
## ma1 0.20668 0.984524 0.209931 0.83372
## omega 0.00000 0.000021 0.021914 0.98252
## alpha1 0.22862 1.066531 0.214361 0.83027
## beta1 0.38053 1.678360 0.226729 0.82064
## beta2 0.38527 0.762624 0.505186 0.61343
##
## LogLikelihood : 4774.315
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.4395
## Bayes -8.4128
## Shibata -8.4396
## Hannan-Quinn -8.4294
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.903 0.0884144
## Lag[2*(p+q)+(p+q)-1][5] 5.845 0.0001447
## Lag[4*(p+q)+(p+q)-1][9] 8.143 0.0511048
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5134 0.4737
## Lag[2*(p+q)+(p+q)-1][8] 2.4167 0.7904
## Lag[4*(p+q)+(p+q)-1][14] 3.2774 0.9388
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.005228 0.500 2.000 0.9424
## ARCH Lag[6] 0.500247 1.461 1.711 0.8919
## ARCH Lag[8] 1.104644 2.368 1.583 0.9084
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 90.4144
## Individual Statistics:
## ar1 0.05410
## ma1 0.04481
## omega 25.57275
## alpha1 0.10454
## beta1 0.21733
## beta2 0.19320
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2710 0.7864
## Negative Sign Bias 0.4915 0.6232
## Positive Sign Bias 0.7169 0.4736
## Joint Effect 0.9071 0.8237
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 64.90 6.332e-07
## 2 30 84.12 2.840e-07
## 3 40 88.55 1.017e-05
## 4 50 106.28 4.088e-06
##
##
## Elapsed time : 0.144074#ARIMA(1,0,1) + GARCH(1,3)
garchSpec13 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(1,3)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE), #ARIMA (1,0,1)
distribution.model = 'norm')
(garchFit13 <- rugarch::ugarchfit(spec = garchSpec13,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,3)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.071428 0.774732 -0.092197 0.926541
## ma1 0.109332 0.770337 0.141928 0.887137
## omega 0.000000 0.000001 0.537105 0.591195
## alpha1 0.260652 0.055557 4.691609 0.000003
## beta1 0.337492 0.183292 1.841282 0.065580
## beta2 0.000128 0.289248 0.000442 0.999647
## beta3 0.396592 0.161355 2.457884 0.013976
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.071428 1.201113 -0.059468 0.95258
## ma1 0.109332 1.205827 0.090670 0.92775
## omega 0.000000 0.000015 0.030467 0.97569
## alpha1 0.260652 0.700881 0.371893 0.70997
## beta1 0.337492 1.507492 0.223877 0.82285
## beta2 0.000128 2.035775 0.000063 0.99995
## beta3 0.396592 1.115245 0.355610 0.72213
##
## LogLikelihood : 4780.394
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.4485
## Bayes -8.4173
## Shibata -8.4486
## Hannan-Quinn -8.4367
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.940 0.086399
## Lag[2*(p+q)+(p+q)-1][5] 5.488 0.000574
## Lag[4*(p+q)+(p+q)-1][9] 7.622 0.078801
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1482 0.7002
## Lag[2*(p+q)+(p+q)-1][11] 1.1994 0.9931
## Lag[4*(p+q)+(p+q)-1][19] 2.0537 0.9994
## d.o.f=4
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[5] 0.118 0.500 2.000 0.7312
## ARCH Lag[7] 0.991 1.473 1.746 0.7598
## ARCH Lag[9] 1.117 2.402 1.619 0.9146
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 112.3846
## Individual Statistics:
## ar1 0.11083
## ma1 0.09386
## omega 28.59319
## alpha1 0.09608
## beta1 0.24776
## beta2 0.28207
## beta3 0.21003
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2166 0.8286
## Negative Sign Bias 0.7442 0.4569
## Positive Sign Bias 0.4500 0.6528
## Joint Effect 0.9852 0.8048
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 59.20 5.169e-06
## 2 30 74.57 7.017e-06
## 3 40 85.65 2.394e-05
## 4 50 100.00 2.365e-05
##
##
## Elapsed time : 0.1719999#ARIMA(1,0,1) + GARCH(2,2)
garchSpec22 <- rugarch::ugarchspec(variance.model = list(model = 'sGARCH',
garchOrder = c(2,2)),
mean.model = list(armaOrder = c(1,1),
include.mean = FALSE), #ARIMA (1,0,1)
distribution.model = 'norm')
(garchFit22 <- rugarch::ugarchfit(spec = garchSpec22,
data = df.train$Return[-1], solver = 'hybrid'))##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(2,2)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.176804 0.570617 -0.309848 0.756677
## ma1 0.218200 0.563285 0.387371 0.698482
## omega 0.000000 0.000000 1.743703 0.081211
## alpha1 0.228800 0.045863 4.988783 0.000001
## alpha2 0.000082 0.033244 0.002481 0.998020
## beta1 0.381936 0.107150 3.564488 0.000365
## beta2 0.383348 0.003044 125.943151 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.176804 0.652696 -0.270883 0.78648
## ma1 0.218200 0.673947 0.323765 0.74612
## omega 0.000000 0.000001 0.381731 0.70266
## alpha1 0.228800 0.180583 1.267010 0.20515
## alpha2 0.000082 0.758177 0.000109 0.99991
## beta1 0.381936 2.438174 0.156649 0.87552
## beta2 0.383348 1.950154 0.196573 0.84416
##
## LogLikelihood : 4774.316
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.4377
## Bayes -8.4066
## Shibata -8.4378
## Hannan-Quinn -8.4260
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 3.148 0.076026
## Lag[2*(p+q)+(p+q)-1][5] 6.104 0.000051
## Lag[4*(p+q)+(p+q)-1][9] 8.404 0.040772
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5158 0.4727
## Lag[2*(p+q)+(p+q)-1][11] 2.9122 0.8789
## Lag[4*(p+q)+(p+q)-1][19] 3.9022 0.9807
## d.o.f=4
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[5] 0.03841 0.500 2.000 0.8446
## ARCH Lag[7] 1.25265 1.473 1.746 0.6877
## ARCH Lag[9] 1.33987 2.402 1.619 0.8796
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 205.3495
## Individual Statistics:
## ar1 0.04547
## ma1 0.03995
## omega 25.25171
## alpha1 0.10409
## alpha2 0.50466
## beta1 0.21627
## beta2 0.19180
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2689 0.7881
## Negative Sign Bias 0.4875 0.6260
## Positive Sign Bias 0.7152 0.4746
## Joint Effect 0.8999 0.8255
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 66.28 3.772e-07
## 2 30 84.44 2.545e-07
## 3 40 86.92 1.648e-05
## 4 50 103.10 1.005e-05
##
##
## Elapsed time : 0.2005482garchFit11@fit$coef## ar1 ma1 omega alpha1 beta1
## 2.833735e-01 -3.125805e-01 3.194177e-07 1.987486e-01 8.002486e-01
## skew shape
## 9.493416e-01 4.911380e+00df_2 <- df.train[-1,] %>%
mutate(res11 = garchFit11@fit$residuals,
sigma11 = garchFit11@fit$sigma)# plot standar deviasi model garch terpilih
df_2 %>%
autoplot(sigma11) +
labs(title = "Simpangan Baku", subtitle = "Model ARIMA (1,0,1) + GARCH (1,1)") +
xlab("Tanggal") + ylab("Simpangan Baku") +
theme_bw()
Hasil Pemodelan Ragam
3.5 Analisis Sisaan Model
garchFit11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.28337 0.462260 0.61302 0.539865
## ma1 -0.31258 0.455683 -0.68596 0.492738
## omega 0.00000 0.000002 0.14542 0.884376
## alpha1 0.19875 0.112883 1.76066 0.078295
## beta1 0.80025 0.089348 8.95657 0.000000
## skew 0.94934 0.035124 27.02863 0.000000
## shape 4.91138 0.845903 5.80608 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.28337 1.262012 0.224541 0.822336
## ma1 -0.31258 1.604692 -0.194792 0.845556
## omega 0.00000 0.000077 0.004166 0.996676
## alpha1 0.19875 4.115027 0.048298 0.961479
## beta1 0.80025 3.240768 0.246932 0.804961
## skew 0.94934 0.197179 4.814630 0.000001
## shape 4.91138 36.907326 0.133073 0.894135
##
## LogLikelihood : 4825.992
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.5292
## Bayes -8.4980
## Shibata -8.5293
## Hannan-Quinn -8.5174
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 14.57 1.347e-04
## Lag[2*(p+q)+(p+q)-1][5] 18.12 0.000e+00
## Lag[4*(p+q)+(p+q)-1][9] 20.50 3.871e-08
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.392 0.2381
## Lag[2*(p+q)+(p+q)-1][5] 3.400 0.3389
## Lag[4*(p+q)+(p+q)-1][9] 4.746 0.4676
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.924 0.500 2.000 0.1654
## ARCH Lag[5] 2.162 1.440 1.667 0.4365
## ARCH Lag[7] 3.211 2.315 1.543 0.4746
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 177.3918
## Individual Statistics:
## ar1 0.7288
## ma1 0.7505
## omega 62.9521
## alpha1 0.1479
## beta1 0.1757
## skew 0.1977
## shape 0.2926
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2569 0.7973
## Negative Sign Bias 0.3447 0.7304
## Positive Sign Bias 0.6230 0.5334
## Joint Effect 0.5676 0.9038
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 16.51 0.6228
## 2 30 23.06 0.7737
## 3 40 30.78 0.8235
## 4 50 48.41 0.4971
##
##
## Elapsed time : 0.442096Pada tabel Weighted Ljung-Box berisi hasil uji Ljung-Box dengan hipotesis: H0: tidak terdapat autokorelasi dalam sisaan H1: terdapat autokorelasi dalam sisaan
plot(garchFit11, which = 'all')##
## please wait...calculating quantiles...
# GARCH(1,1)
a = ggplot(df_2, aes(Tanggal, res11)) +
geom_line() +
ylab("Residuals") +xlab("Tanggal") +
labs(title = "Residual ARIMA (1,0,1) + GARCH (1,1)")
b = ACF(df_2, res11) %>%
autoplot(alpha = .2, fill = '#FF6666')## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.c = ggplot(df_2, aes(res11)) +
geom_histogram(aes(y = ..density..), color = 'black', fill = "white") +
geom_density() +
xlab("Residuals")
d = ggplot(df_2, aes(sample = res11)) +
stat_qq() +
stat_qq_line()
ggarrange(a, ggarrange(b, ggarrange(c,d)), nrow = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
pred.roll <- ugarchroll(spec = garchSpec11,
data = df$Return,
n.ahead = 1,
forecast.length = dim(df.test)[1],
#refit.window = "recursive"
window.wize = length(df$Return))# extract output
pred <- as.data.frame(pred.roll, which = "density")
plot.ts(pred$Mu, xlab = "Waktu(hari)", ylab = "Return Kurs USD/IDR", main = "Predicted-Value ARIMA(1,0,1) + GARCH(1,1)")
plot.ts(df.test$Return, xlab = "Waktu(hari)", ylab = "Return kurs USD/IDR", main = "Data Testing")
# tambahkan forecast diff dan sigma ke data test
df.test <-
df.test %>%
mutate(diff = pred$Mu,
Kurs.lag = dplyr::lag(Return, 1), #lag 1 peubah Return
sigma = pred$Sigma)
# observasi pertama pada test jadi terakhir pada train
df.test[1,"Kurs.lag"] <- df.train[dim(df.train)[1], "Return"]
df.test <-
df.test %>%
mutate(Forecast = diff + Kurs.lag,
residual = Return - Forecast,
low95 = Forecast - (1.96*sigma),
up95 = Forecast + (1.96*sigma))
#Plot
point.forecast <-
df.test %>%
select(Return, Forecast) %>%
pivot_longer(-Tanggal)
ggplot() +
geom_line(point.forecast, mapping = aes(Tanggal, value, color = point.forecast$name)) +
geom_ribbon(df.test, mapping = aes(Tanggal, ymax = up95, ymin = low95), fill = "grey", alpha = 0.5)+
xlab("Tanggal")+
labs(title = "Nilai Sebenarnya vs Nilai Forecasting",
subtitle = "Return Kurs USD/IDR")+
theme(legend.title = element_blank())## Warning: Use of `point.forecast$name` is discouraged.
## ℹ Use `name` instead.
# Ukuran Keakuratan
error <- df.test$Return - df.test$Forecast
## SSE : SUM SQUARE ERROR
sse<- sum(error^2)
# MSE: MEAN SQUARED ERROR
mse<- mean(error^2)
# RMSE: Root Mean Square Error
rmse<- sqrt(mse)
# MAD: mean absolute deviation
mad<- mean(abs(error))
## MAPE
r.error<- (error/df.test$Return)*100 #relative error
r.error[!is.finite(r.error)] <- NA
mape <- mean(abs(r.error), na.rm = T)
akurasi<- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse, mse, mape, rmse, mad))
format(akurasi, scientific = F)## Ukuran.Keakuratan Nilai
## 1 SSE 0.00538352403
## 2 MSE 0.00001902305
## 3 MAPE 245.15839649140
## 4 RMSE 0.00436154251
## 5 MAD 0.00323074278invers <- function(data1, data2){
idx = index(data1)
kurs <- c()
for (j in 1:(length(idx)-1)){
kurs[j] = exp(data1[j] + log(data2[j+1]))
}
return(kurs)
}pred_kurs = invers(df.test$Forecast, df.test$Kurs)
head(pred_kurs)## [1] 14879.64 14977.52 14765.77 14772.11 14815.98 14914.41# Ukuran Keakuratan
error <- df.test$Kurs[-1] - pred_kurs
## SSE : SUM SQUARE ERROR
sse<- sum(error^2)
# MSE: MEAN SQUARED ERROR
mse<- mean(error^2)
# RMSE: Root Mean Square Error
rmse<- sqrt(mse)
# MAD: mean absolute deviation
mad<- mean(abs(error))
## MAPE
r.error<- (error/df.test$Kurs[-1])*100 #relative error
r.error[!is.finite(r.error)] <- NA
mape <- mean(abs(r.error), na.rm = T)
akurasi<- data.frame("Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Nilai" = c(sse, mse, mape, rmse, mad))
format(akurasi, scientific = F)## Ukuran.Keakuratan Nilai
## 1 SSE 653085.8753615
## 2 MSE 2315.9073594
## 3 MAPE 0.2431013
## 4 RMSE 48.1238751
## 5 MAD 34.9106192ts.plot(ts(pred_kurs), col = "red")
lines(df.test$Kurs[-1])