library("forecast")
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library("TTR")
library("TSA")
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library("graphics")
library("tseries")
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, union
data <- read.csv("latihan time series.csv")
data1 <- data[-199:-397,]
data2 <- data[-1:-198]
data.ts <- ts(data1)
head(data.ts)
## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## DATE IPG2211A2N
## 1 1 72.5052
## 2 66 70.6720
## 3 83 62.4502
## 4 100 57.4714
## 5 117 55.3151
## 6 134 58.0904
ts.plot(data.ts, xlab="Time Period ", ylab= "Sales", main= "Time Series Plot Data Sales")
points(data.ts)
data.training<-subset(data.ts, start = 1, end = 99)
data.training
## Time Series:
## Start = 1
## End = 99
## Frequency = 1
## DATE IPG2211A2N
## 1 1 72.5052
## 2 66 70.6720
## 3 83 62.4502
## 4 100 57.4714
## 5 117 55.3151
## 6 134 58.0904
## 7 151 62.6202
## 8 167 63.2485
## 9 183 60.5846
## 10 18 56.3154
## 11 34 58.0005
## 12 50 68.7145
## 13 2 73.3057
## 14 67 67.9869
## 15 84 62.2221
## 16 101 57.0329
## 17 118 55.8137
## 18 135 59.9005
## 19 152 65.7655
## 20 168 64.4816
## 21 184 61.0005
## 22 19 57.5322
## 23 35 59.3417
## 24 51 68.1354
## 25 3 73.8152
## 26 68 70.0620
## 27 85 65.6100
## 28 102 60.1586
## 29 119 58.8734
## 30 136 63.8918
## 31 153 68.8694
## 32 169 70.0669
## 33 185 64.1151
## 34 20 60.3789
## 35 36 62.4643
## 36 52 70.5777
## 37 4 79.8703
## 38 69 76.1622
## 39 86 70.2928
## 40 103 63.2384
## 41 120 61.4065
## 42 137 67.1097
## 43 154 72.9816
## 44 170 75.7655
## 45 186 67.5152
## 46 21 63.2832
## 47 37 65.1078
## 48 53 73.8631
## 49 5 77.9188
## 50 70 76.6822
## 51 87 73.3523
## 52 104 65.1081
## 53 121 63.6892
## 54 138 68.4722
## 55 155 74.0301
## 56 171 75.0448
## 57 187 69.3053
## 58 22 65.8735
## 59 38 69.0706
## 60 54 84.1949
## 61 6 84.3598
## 62 71 77.1726
## 63 88 73.1964
## 64 105 67.2781
## 65 122 65.8218
## 66 139 71.4654
## 67 156 76.6140
## 68 172 77.1052
## 69 188 73.0610
## 70 23 67.4365
## 71 39 68.5665
## 72 55 77.6839
## 73 7 86.0214
## 74 72 77.5573
## 75 89 73.3650
## 76 106 67.1500
## 77 123 68.8162
## 78 140 74.8448
## 79 157 80.0928
## 80 173 79.1606
## 81 189 73.5743
## 82 24 68.7538
## 83 40 72.5166
## 84 56 79.4894
## 85 8 85.2855
## 86 73 80.1643
## 87 90 74.5275
## 88 107 69.6441
## 89 124 67.1784
## 90 141 71.2078
## 91 158 77.5081
## 92 174 76.5374
## 93 190 72.3541
## 94 25 69.0286
## 95 41 73.4992
## 96 57 84.5159
## 97 9 87.9464
## 98 74 84.5561
## 99 91 79.4747
data.test <- subset(data.ts, start = 100, end = 198)
data.test
## Time Series:
## Start = 100
## End = 198
## Frequency = 1
## DATE IPG2211A2N
## 100 108 71.0578
## 101 125 67.6762
## 102 142 74.3297
## 103 159 82.1048
## 104 175 82.0605
## 105 191 74.6031
## 106 26 69.6810
## 107 42 74.4292
## 108 58 84.2284
## 109 10 94.1386
## 110 75 87.1607
## 111 92 79.2456
## 112 109 70.9749
## 113 126 69.3844
## 114 143 77.9831
## 115 160 83.2770
## 116 176 81.8872
## 117 192 75.6826
## 118 27 71.2661
## 119 43 75.2458
## 120 59 84.8147
## 121 11 92.4532
## 122 76 87.4033
## 123 93 81.2661
## 124 110 73.8167
## 125 127 73.2682
## 126 144 78.3026
## 127 161 85.9841
## 128 177 89.5467
## 129 193 78.5035
## 130 28 73.7066
## 131 44 79.6543
## 132 60 90.8251
## 133 12 98.9732
## 134 77 92.8883
## 135 94 86.9356
## 136 111 77.2214
## 137 128 76.6826
## 138 145 81.9306
## 139 162 85.9606
## 140 178 86.5562
## 141 194 79.1919
## 142 29 74.6891
## 143 45 81.0740
## 144 61 90.4855
## 145 13 98.4613
## 146 78 89.7795
## 147 95 83.0125
## 148 112 76.1476
## 149 129 73.8471
## 150 146 79.7645
## 151 163 88.4519
## 152 179 87.7828
## 153 195 81.9386
## 154 30 77.5027
## 155 46 82.0448
## 156 62 92.1010
## 157 14 94.7920
## 158 79 87.8200
## 159 96 86.5549
## 160 113 76.7521
## 161 130 78.0303
## 162 147 86.4579
## 163 164 93.8379
## 164 180 93.5310
## 165 196 87.5414
## 166 31 80.0924
## 167 47 81.4349
## 168 63 91.6841
## 169 15 102.1348
## 170 80 91.1829
## 171 97 90.7381
## 172 114 80.5176
## 173 131 79.3887
## 174 148 87.8431
## 175 165 97.4903
## 176 181 96.4157
## 177 197 87.2248
## 178 32 80.6409
## 179 48 82.2025
## 180 64 94.5113
## 181 16 102.2301
## 182 81 94.2989
## 183 98 88.0927
## 184 115 81.4425
## 185 132 84.4552
## 186 149 91.0406
## 187 166 95.9957
## 188 182 99.3704
## 189 198 90.9178
## 190 33 83.1408
## 191 49 88.0410
## 192 65 102.4558
## 193 17 109.1081
## 194 82 97.1717
## 195 99 92.8283
## 196 116 82.9150
## 197 133 82.5465
## 198 150 90.3955
plot.ts(data.training, main = "Data Training Produksi Listrik", xlab = "pertanggal", ylab = "Produksi Listrik (MWh)")
plot.ts(data.test, main = "Data Training Produksi Listrik", xlab = "Pertanggal", ylab = "Produksi Listrik (MWh)")
# Model ARIMA
data_vektor <- as.vector(data.ts)
print(data_vektor)
## [1] 1.0000 66.0000 83.0000 100.0000 117.0000 134.0000 151.0000 167.0000
## [9] 183.0000 18.0000 34.0000 50.0000 2.0000 67.0000 84.0000 101.0000
## [17] 118.0000 135.0000 152.0000 168.0000 184.0000 19.0000 35.0000 51.0000
## [25] 3.0000 68.0000 85.0000 102.0000 119.0000 136.0000 153.0000 169.0000
## [33] 185.0000 20.0000 36.0000 52.0000 4.0000 69.0000 86.0000 103.0000
## [41] 120.0000 137.0000 154.0000 170.0000 186.0000 21.0000 37.0000 53.0000
## [49] 5.0000 70.0000 87.0000 104.0000 121.0000 138.0000 155.0000 171.0000
## [57] 187.0000 22.0000 38.0000 54.0000 6.0000 71.0000 88.0000 105.0000
## [65] 122.0000 139.0000 156.0000 172.0000 188.0000 23.0000 39.0000 55.0000
## [73] 7.0000 72.0000 89.0000 106.0000 123.0000 140.0000 157.0000 173.0000
## [81] 189.0000 24.0000 40.0000 56.0000 8.0000 73.0000 90.0000 107.0000
## [89] 124.0000 141.0000 158.0000 174.0000 190.0000 25.0000 41.0000 57.0000
## [97] 9.0000 74.0000 91.0000 108.0000 125.0000 142.0000 159.0000 175.0000
## [105] 191.0000 26.0000 42.0000 58.0000 10.0000 75.0000 92.0000 109.0000
## [113] 126.0000 143.0000 160.0000 176.0000 192.0000 27.0000 43.0000 59.0000
## [121] 11.0000 76.0000 93.0000 110.0000 127.0000 144.0000 161.0000 177.0000
## [129] 193.0000 28.0000 44.0000 60.0000 12.0000 77.0000 94.0000 111.0000
## [137] 128.0000 145.0000 162.0000 178.0000 194.0000 29.0000 45.0000 61.0000
## [145] 13.0000 78.0000 95.0000 112.0000 129.0000 146.0000 163.0000 179.0000
## [153] 195.0000 30.0000 46.0000 62.0000 14.0000 79.0000 96.0000 113.0000
## [161] 130.0000 147.0000 164.0000 180.0000 196.0000 31.0000 47.0000 63.0000
## [169] 15.0000 80.0000 97.0000 114.0000 131.0000 148.0000 165.0000 181.0000
## [177] 197.0000 32.0000 48.0000 64.0000 16.0000 81.0000 98.0000 115.0000
## [185] 132.0000 149.0000 166.0000 182.0000 198.0000 33.0000 49.0000 65.0000
## [193] 17.0000 82.0000 99.0000 116.0000 133.0000 150.0000 72.5052 70.6720
## [201] 62.4502 57.4714 55.3151 58.0904 62.6202 63.2485 60.5846 56.3154
## [209] 58.0005 68.7145 73.3057 67.9869 62.2221 57.0329 55.8137 59.9005
## [217] 65.7655 64.4816 61.0005 57.5322 59.3417 68.1354 73.8152 70.0620
## [225] 65.6100 60.1586 58.8734 63.8918 68.8694 70.0669 64.1151 60.3789
## [233] 62.4643 70.5777 79.8703 76.1622 70.2928 63.2384 61.4065 67.1097
## [241] 72.9816 75.7655 67.5152 63.2832 65.1078 73.8631 77.9188 76.6822
## [249] 73.3523 65.1081 63.6892 68.4722 74.0301 75.0448 69.3053 65.8735
## [257] 69.0706 84.1949 84.3598 77.1726 73.1964 67.2781 65.8218 71.4654
## [265] 76.6140 77.1052 73.0610 67.4365 68.5665 77.6839 86.0214 77.5573
## [273] 73.3650 67.1500 68.8162 74.8448 80.0928 79.1606 73.5743 68.7538
## [281] 72.5166 79.4894 85.2855 80.1643 74.5275 69.6441 67.1784 71.2078
## [289] 77.5081 76.5374 72.3541 69.0286 73.4992 84.5159 87.9464 84.5561
## [297] 79.4747 71.0578 67.6762 74.3297 82.1048 82.0605 74.6031 69.6810
## [305] 74.4292 84.2284 94.1386 87.1607 79.2456 70.9749 69.3844 77.9831
## [313] 83.2770 81.8872 75.6826 71.2661 75.2458 84.8147 92.4532 87.4033
## [321] 81.2661 73.8167 73.2682 78.3026 85.9841 89.5467 78.5035 73.7066
## [329] 79.6543 90.8251 98.9732 92.8883 86.9356 77.2214 76.6826 81.9306
## [337] 85.9606 86.5562 79.1919 74.6891 81.0740 90.4855 98.4613 89.7795
## [345] 83.0125 76.1476 73.8471 79.7645 88.4519 87.7828 81.9386 77.5027
## [353] 82.0448 92.1010 94.7920 87.8200 86.5549 76.7521 78.0303 86.4579
## [361] 93.8379 93.5310 87.5414 80.0924 81.4349 91.6841 102.1348 91.1829
## [369] 90.7381 80.5176 79.3887 87.8431 97.4903 96.4157 87.2248 80.6409
## [377] 82.2025 94.5113 102.2301 94.2989 88.0927 81.4425 84.4552 91.0406
## [385] 95.9957 99.3704 90.9178 83.1408 88.0410 102.4558 109.1081 97.1717
## [393] 92.8283 82.9150 82.5465 90.3955
adf.test(data_vektor)
## Warning in adf.test(data_vektor): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: data_vektor
## Dickey-Fuller = -12.084, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
acf(data.ts)
eacf(data_vektor)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o x x x x x x x x x x x
## 1 x x o x o o x x x o o x x x
## 2 x o x x x o o x x o o x x o
## 3 x o x x x o o x x o o x x o
## 4 x x o x o o o o o o o x x x
## 5 x x o o x x x o o o o x x x
## 6 x x o o x x x o o o o x x x
## 7 x x x o x x x o o o o x x x
str(data.ts)
## Time-Series [1:198, 1:2] from 1 to 198: 1 66 83 100 117 134 151 167 183 18 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "DATE" "IPG2211A2N"
arima(data_vektor, order = c(2,0,2), method = "ML")
##
## Call:
## arima(x = data_vektor, order = c(2, 0, 2), method = "ML")
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## -0.0367 0.3397 0.7497 -0.0859 88.1807
## s.e. 0.0987 0.0559 0.0929 0.0711 4.0007
##
## sigma^2 estimated as 1119: log likelihood = -1952.31, aic = 3914.62
arima(data_vektor, order = c(4,0,5), method = "ML")
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
##
## Call:
## arima(x = data_vektor, order = c(4, 0, 5), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 0.0048 1.0045 -0.0102 -0.9836 0.6244 -1.0620 -0.6329 0.8901
## s.e. 0.0077 0.0094 0.0099 0.0108 0.0487 0.0519 0.0705 0.0224
## ma5 intercept
## 0.5433 88.4926
## s.e. 0.0534 1.6673
##
## sigma^2 estimated as 574.7: log likelihood = -1823.8, aic = 3667.59
arima(data_vektor, order = c(4,0,9), method = "ML")
##
## Call:
## arima(x = data_vektor, order = c(4, 0, 9), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 0.3787 -0.283 0.6891 -0.6427 0.2540 0.5328 -0.6643 -0.1175
## s.e. 0.0448 0.027 0.0274 0.0474 0.0449 0.0369 0.0389 0.0556
## ma5 ma6 ma7 ma8 ma9 intercept
## -0.3458 -0.1908 0.4518 0.4375 0.7083 88.6564
## s.e. 0.0456 0.0448 0.0408 0.0391 0.0314 2.5360
##
## sigma^2 estimated as 446.7: log likelihood = -1778.53, aic = 3585.05
arima(data_vektor, order = c(4,0,10), method = "ML")
##
## Call:
## arima(x = data_vektor, order = c(4, 0, 10), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## -0.6223 -0.1560 0.7842 0.4328 1.3337 0.9349 -0.2971 -1.2775
## s.e. 0.0754 0.0534 0.0534 0.0768 0.0623 0.0626 0.0612 0.0793
## ma5 ma6 ma7 ma8 ma9 ma10 intercept
## -1.3042 -1.1184 -0.1338 0.8353 1.1211 0.7665 88.6699
## s.e. 0.0950 0.0696 0.0580 0.0593 0.0557 0.0269 3.4272
##
## sigma^2 estimated as 443.9: log likelihood = -1774.8, aic = 3579.59
arima(data_vektor, order = c(5,0,9), method = "ML")
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
##
## Call:
## arima(x = data_vektor, order = c(5, 0, 9), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 1.4519 -0.7922 0.9197 -1.5169 0.7234 -1.0495 0.2906 -0.9951
## s.e. 0.0976 0.0729 0.0311 0.0788 0.0961 0.0982 0.0786 0.0791
## ma4 ma5 ma6 ma7 ma8 ma9 intercept
## 1.2821 0.2807 -0.3048 -0.2260 -0.3801 0.6019 88.6791
## s.e. 0.1313 0.1178 0.0932 0.0871 0.0802 0.0471 2.1697
##
## sigma^2 estimated as 350.5: log likelihood = -1731.63, aic = 3493.26
arima(data_vektor, order = c(6,0,10), method = "ML")
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
##
## Call:
## arima(x = data_vektor, order = c(6, 0, 10), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ma1 ma2
## -0.9852 0.018 0.9940 -0.0065 -0.9893 -0.9826 1.9547 1.3810
## s.e. 0.0088 0.010 0.0127 0.0099 0.0082 0.0077 0.0512 0.1141
## ma3 ma4 ma5 ma6 ma7 ma8 ma9 ma10
## -0.6498 -1.2549 0.1664 1.7286 1.4333 0.3403 -0.2062 -0.0762
## s.e. 0.1299 0.1376 0.1394 0.0753 0.1389 0.1217 0.0937 0.0521
## intercept
## 88.5862
## s.e. 1.6505
##
## sigma^2 estimated as 278.1: log likelihood = -1685.58, aic = 3405.17
Berdasarkkan nilai AIC, model yang mempunyai AIC terkecil adalah model ARIMA (6,0,10) sehingga model yang digunakan untuk peramalan adalah model ARIMA (6,0,10)
ARIMA_data <- Arima(data_vektor, order = c(6,0,10), method = "ML")
hasil_forecastt <- forecast(ARIMA_data)
hasil_forecast <- as.data.frame(hasil_forecastt)
hasil_forecast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 397 90.77812 68.93297 112.6233 57.36884 124.1874
## 398 90.79803 60.37322 121.2228 44.26730 137.3288
## 399 87.14854 55.21591 119.0812 38.31180 135.9853
## 400 89.77177 57.79634 121.7472 40.86957 138.6740
## 401 89.65454 57.33527 121.9738 40.22648 139.0826
## 402 91.43834 58.88795 123.9887 41.65681 141.2199
## 403 92.33596 59.67955 124.9924 42.39230 142.2796
## 404 87.65181 54.97219 120.3314 37.67264 137.6310
## 405 91.34001 58.51275 124.1673 41.13505 141.5450
## 406 86.86693 53.96433 119.7695 36.54674 137.1871
plot(hasil_forecastt)
# Penilaian Akurasi dengan MAPE Penilaian akurasi dilakukan dengan
membandingkan data hasil ramalan dengan data aktual 10 amatan
terakhir.
data_forecast1 <- hasil_forecast$`Point Forecast`
data_forecast1 <- as.data.frame(data_forecast1)
head(data_forecast1)
## data_forecast1
## 1 90.77812
## 2 90.79803
## 3 87.14854
## 4 89.77177
## 5 89.65454
## 6 91.43834
dataaktual <- as.data.frame(data2)
dataaktual
## data frame with 0 columns and 397 rows
plot(data.ts, main = "Plot semua data")
points(data.ts)
data1.ses <- HoltWinters(data.ts, alpha = 0.1, beta = FALSE, gamma = FALSE)
data1.ses
## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = data.ts, alpha = 0.1, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.1
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 91.0391
SSE <- data1.ses$SSE
SSE
## [1] 751711.7
MSE <- SSE/199
MSE
## [1] 3777.446
data1.ses <- as.ts(data1.ses)
data1.ses <- ses(data_vektor)
hasil_forecast.sess <- forecast(data1.ses)
hasil_forecast.ses <- as.data.frame(hasil_forecast.sess)
hasil_forecast.ses
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 397 89.03974 39.828876 138.2506 13.77822 164.3013
## 398 89.03974 25.490548 152.5889 -8.15036 186.2298
## 399 89.03974 13.838075 164.2414 -25.97128 204.0508
## 400 89.03974 3.763242 174.3162 -41.37941 219.4589
## 401 89.03974 -5.241073 183.3205 -55.15032 233.2298
## 402 89.03974 -13.457393 191.5369 -67.71610 245.7956
## 403 89.03974 -21.062271 199.1417 -79.34675 257.4262
## 404 89.03974 -28.174779 206.2543 -90.22440 268.3039
## 405 89.03974 -34.879726 212.9592 -100.47873 278.5582
## 406 89.03974 -41.240054 219.3195 -110.20601 288.2855
plot(hasil_forecast)
data_forecast2 <- hasil_forecast.ses$`Point Forecast`
data_forecast2 <- as.data.frame(data_forecast2)
head(data_forecast2)
## data_forecast2
## 1 89.03974
## 2 89.03974
## 3 89.03974
## 4 89.03974
## 5 89.03974
## 6 89.03974
dataaktual <- as.data.frame(data2)
dataaktual
## data frame with 0 columns and 397 rows
akurasi.ses <- accuracy(data_forecast2$data_forecast2, dataaktual[["Sales (in thousand)"]])
akurasi.ses
## ME RMSE MAE MPE MAPE
## Test set NaN NaN NaN NaN NaN
Setelah dilakukan peramalan dengan menggunakan dua metode, akan diamati nilai MAPE dari masing-masing metode. Metode dengan nilai MAPE lebih kecil adalah metode ARIMA(6,0,10) dengan p = 6, d = 0 q = 10
hasil_forecastt <- forecast(ARIMA_data, h=10)
akurasi.arima <- accuracy(hasil_forecastt$mean, dataaktual$`Sales (in thousand)`)
print(akurasi.arima)
## ME RMSE MAE MPE MAPE
## Test set NaN NaN NaN NaN NaN
akurasi.ses
## ME RMSE MAE MPE MAPE
## Test set NaN NaN NaN NaN NaN