Time Series dengan menggunakan ARIMA

A. Menyiapkan packages yang diperlukan

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

B. Plot Data

ts.plot(data.ts, xlab="Time Period ", ylab= "Sales", main= "Time Series Plot Data Sales")
points(data.ts)

C. Partisi Data

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 Data Partisi

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

Penentuan Model Terbaik Berdasarkan AIC

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")

Peramalan 10 Waktu Kedepan

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

Simple Exponential Smoothing

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

Perbandingan model ARIMA dan Exponential Smoothing Forecast

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