Tugas Minggu 8

Import Library

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.2

## Warning: package 'ggplot2' was built under R version 4.3.3

## Warning: package 'tibble' was built under R version 4.3.2

## Warning: package 'tidyr' was built under R version 4.3.2

## Warning: package 'readr' was built under R version 4.3.2

## Warning: package 'purrr' was built under R version 4.3.2

## Warning: package 'dplyr' was built under R version 4.3.2

## Warning: package 'forcats' was built under R version 4.3.2

## Warning: package 'lubridate' was built under R version 4.3.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(forecast)

## Warning: package 'forecast' was built under R version 4.3.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(TSA)

## Warning: package 'TSA' was built under R version 4.3.3

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## 
## The following object is masked from 'package:readr':
## 
##     spec
## 
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## 
## The following object is masked from 'package:utils':
## 
##     tar

library(aTSA)

## Warning: package 'aTSA' was built under R version 4.3.2

## 
## Attaching package: 'aTSA'
## 
## The following object is masked from 'package:forecast':
## 
##     forecast
## 
## The following object is masked from 'package:graphics':
## 
##     identify

library(car)

## Warning: package 'car' was built under R version 4.3.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.2

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.2

## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(readxl)

## Warning: package 'readxl' was built under R version 4.3.2

Deksripsi Data

Data yang diambil adalah Data Konsumsi Listrik dari Tahun 1985 Hingga 2018. Pengunaan Listrik ini diukur dengan satuan KW/H ( KilloWatt Hour) yang merupakan agregasi. Dataset ini bersumber dari Kaggle. Pada Tugas Kali ini saya akan memakai 150 periode dari dataset ini sebagai data dalam melakukan Forecasting seasonalnya.

Data<- read.csv("C:/Users/Admin/Downloads/Electric Production.csv")
Data

##           DATE    Value
## 1   01-01-1985  72.5052
## 2   02-01-1985  70.6720
## 3   03-01-1985  62.4502
## 4   04-01-1985  57.4714
## 5   05-01-1985  55.3151
## 6   06-01-1985  58.0904
## 7   07-01-1985  62.6202
## 8   08-01-1985  63.2485
## 9   09-01-1985  60.5846
## 10  10-01-1985  56.3154
## 11  11-01-1985  58.0005
## 12  12-01-1985  68.7145
## 13  01-01-1986  73.3057
## 14  02-01-1986  67.9869
## 15  03-01-1986  62.2221
## 16  04-01-1986  57.0329
## 17  05-01-1986  55.8137
## 18  06-01-1986  59.9005
## 19  07-01-1986  65.7655
## 20  08-01-1986  64.4816
## 21  09-01-1986  61.0005
## 22  10-01-1986  57.5322
## 23  11-01-1986  59.3417
## 24  12-01-1986  68.1354
## 25  01-01-1987  73.8152
## 26  02-01-1987  70.0620
## 27  03-01-1987  65.6100
## 28  04-01-1987  60.1586
## 29  05-01-1987  58.8734
## 30  06-01-1987  63.8918
## 31  07-01-1987  68.8694
## 32  08-01-1987  70.0669
## 33  09-01-1987  64.1151
## 34  10-01-1987  60.3789
## 35  11-01-1987  62.4643
## 36  12-01-1987  70.5777
## 37  01-01-1988  79.8703
## 38  02-01-1988  76.1622
## 39  03-01-1988  70.2928
## 40  04-01-1988  63.2384
## 41  05-01-1988  61.4065
## 42  06-01-1988  67.1097
## 43  07-01-1988  72.9816
## 44  08-01-1988  75.7655
## 45  09-01-1988  67.5152
## 46  10-01-1988  63.2832
## 47  11-01-1988  65.1078
## 48  12-01-1988  73.8631
## 49  01-01-1989  77.9188
## 50  02-01-1989  76.6822
## 51  03-01-1989  73.3523
## 52  04-01-1989  65.1081
## 53  05-01-1989  63.6892
## 54  06-01-1989  68.4722
## 55  07-01-1989  74.0301
## 56  08-01-1989  75.0448
## 57  09-01-1989  69.3053
## 58  10-01-1989  65.8735
## 59  11-01-1989  69.0706
## 60  12-01-1989  84.1949
## 61  01-01-1990  84.3598
## 62  02-01-1990  77.1726
## 63  03-01-1990  73.1964
## 64  04-01-1990  67.2781
## 65  05-01-1990  65.8218
## 66  06-01-1990  71.4654
## 67  07-01-1990  76.6140
## 68  08-01-1990  77.1052
## 69  09-01-1990  73.0610
## 70  10-01-1990  67.4365
## 71  11-01-1990  68.5665
## 72  12-01-1990  77.6839
## 73  01-01-1991  86.0214
## 74  02-01-1991  77.5573
## 75  03-01-1991  73.3650
## 76  04-01-1991  67.1500
## 77  05-01-1991  68.8162
## 78  06-01-1991  74.8448
## 79  07-01-1991  80.0928
## 80  08-01-1991  79.1606
## 81  09-01-1991  73.5743
## 82  10-01-1991  68.7538
## 83  11-01-1991  72.5166
## 84  12-01-1991  79.4894
## 85  01-01-1992  85.2855
## 86  02-01-1992  80.1643
## 87  03-01-1992  74.5275
## 88  04-01-1992  69.6441
## 89  05-01-1992  67.1784
## 90  06-01-1992  71.2078
## 91  07-01-1992  77.5081
## 92  08-01-1992  76.5374
## 93  09-01-1992  72.3541
## 94  10-01-1992  69.0286
## 95  11-01-1992  73.4992
## 96  12-01-1992  84.5159
## 97  01-01-1993  87.9464
## 98  02-01-1993  84.5561
## 99  03-01-1993  79.4747
## 100 04-01-1993  71.0578
## 101 05-01-1993  67.6762
## 102 06-01-1993  74.3297
## 103 07-01-1993  82.1048
## 104 08-01-1993  82.0605
## 105 09-01-1993  74.6031
## 106 10-01-1993  69.6810
## 107 11-01-1993  74.4292
## 108 12-01-1993  84.2284
## 109 01-01-1994  94.1386
## 110 02-01-1994  87.1607
## 111 03-01-1994  79.2456
## 112 04-01-1994  70.9749
## 113 05-01-1994  69.3844
## 114 06-01-1994  77.9831
## 115 07-01-1994  83.2770
## 116 08-01-1994  81.8872
## 117 09-01-1994  75.6826
## 118 10-01-1994  71.2661
## 119 11-01-1994  75.2458
## 120 12-01-1994  84.8147
## 121 01-01-1995  92.4532
## 122 02-01-1995  87.4033
## 123 03-01-1995  81.2661
## 124 04-01-1995  73.8167
## 125 05-01-1995  73.2682
## 126 06-01-1995  78.3026
## 127 07-01-1995  85.9841
## 128 08-01-1995  89.5467
## 129 09-01-1995  78.5035
## 130 10-01-1995  73.7066
## 131 11-01-1995  79.6543
## 132 12-01-1995  90.8251
## 133 01-01-1996  98.9732
## 134 02-01-1996  92.8883
## 135 03-01-1996  86.9356
## 136 04-01-1996  77.2214
## 137 05-01-1996  76.6826
## 138 06-01-1996  81.9306
## 139 07-01-1996  85.9606
## 140 08-01-1996  86.5562
## 141 09-01-1996  79.1919
## 142 10-01-1996  74.6891
## 143 11-01-1996  81.0740
## 144 12-01-1996  90.4855
## 145 01-01-1997  98.4613
## 146 02-01-1997  89.7795
## 147 03-01-1997  83.0125
## 148 04-01-1997  76.1476
## 149 05-01-1997  73.8471
## 150 06-01-1997  79.7645
## 151 07-01-1997  88.4519
## 152 08-01-1997  87.7828
## 153 09-01-1997  81.9386
## 154 10-01-1997  77.5027
## 155 11-01-1997  82.0448
## 156 12-01-1997  92.1010
## 157 01-01-1998  94.7920
## 158 02-01-1998  87.8200
## 159 03-01-1998  86.5549
## 160 04-01-1998  76.7521
## 161 05-01-1998  78.0303
## 162 06-01-1998  86.4579
## 163 07-01-1998  93.8379
## 164 08-01-1998  93.5310
## 165 09-01-1998  87.5414
## 166 10-01-1998  80.0924
## 167 11-01-1998  81.4349
## 168 12-01-1998  91.6841
## 169 01-01-1999 102.1348
## 170 02-01-1999  91.1829
## 171 03-01-1999  90.7381
## 172 04-01-1999  80.5176
## 173 05-01-1999  79.3887
## 174 06-01-1999  87.8431
## 175 07-01-1999  97.4903
## 176 08-01-1999  96.4157
## 177 09-01-1999  87.2248
## 178 10-01-1999  80.6409
## 179 11-01-1999  82.2025
## 180 12-01-1999  94.5113
## 181 01-01-2000 102.2301
## 182 02-01-2000  94.2989
## 183 03-01-2000  88.0927
## 184 04-01-2000  81.4425
## 185 05-01-2000  84.4552
## 186 06-01-2000  91.0406
## 187 07-01-2000  95.9957
## 188 08-01-2000  99.3704
## 189 09-01-2000  90.9178
## 190 10-01-2000  83.1408
## 191 11-01-2000  88.0410
## 192 12-01-2000 102.4558
## 193 01-01-2001 109.1081
## 194 02-01-2001  97.1717
## 195 03-01-2001  92.8283
## 196 04-01-2001  82.9150
## 197 05-01-2001  82.5465
## 198 06-01-2001  90.3955
## 199 07-01-2001  96.0740
## 200 08-01-2001  99.5534
## 201 09-01-2001  88.2810
## 202 10-01-2001  82.6860
## 203 11-01-2001  82.9319
## 204 12-01-2001  93.0381
## 205 01-01-2002 102.9955
## 206 02-01-2002  95.2075
## 207 03-01-2002  93.2556
## 208 04-01-2002  85.7950
## 209 05-01-2002  85.2351
## 210 06-01-2002  93.1896
## 211 07-01-2002 102.3930
## 212 08-01-2002 101.6293
## 213 09-01-2002  93.3089
## 214 10-01-2002  86.9002
## 215 11-01-2002  88.5749
## 216 12-01-2002 100.8003
## 217 01-01-2003 110.1807
## 218 02-01-2003 103.8413
## 219 03-01-2003  94.5532
## 220 04-01-2003  85.0620
## 221 05-01-2003  85.4653
## 222 06-01-2003  91.0761
## 223 07-01-2003 102.2200
## 224 08-01-2003 104.4682
## 225 09-01-2003  92.9135
## 226 10-01-2003  86.5047
## 227 11-01-2003  88.5735
## 228 12-01-2003 103.5428
## 229 01-01-2004 113.7226
## 230 02-01-2004 106.1590
## 231 03-01-2004  95.4029
## 232 04-01-2004  86.7233
## 233 05-01-2004  89.0302
## 234 06-01-2004  95.5045
## 235 07-01-2004 101.7948
## 236 08-01-2004 100.2025
## 237 09-01-2004  94.0240
## 238 10-01-2004  87.5262
## 239 11-01-2004  89.6144
## 240 12-01-2004 105.7263
## 241 01-01-2005 111.1614
## 242 02-01-2005 101.7795
## 243 03-01-2005  98.9565
## 244 04-01-2005  86.4776
## 245 05-01-2005  87.2234
## 246 06-01-2005  99.5076
## 247 07-01-2005 108.3501
## 248 08-01-2005 109.4862
## 249 09-01-2005  99.1155
## 250 10-01-2005  89.7567
## 251 11-01-2005  90.4587
## 252 12-01-2005 108.2257
## 253 01-01-2006 104.4724
## 254 02-01-2006 101.5196
## 255 03-01-2006  98.4017
## 256 04-01-2006  87.5093
## 257 05-01-2006  90.0222
## 258 06-01-2006 100.5244
## 259 07-01-2006 110.9503
## 260 08-01-2006 111.5192
## 261 09-01-2006  95.7632
## 262 10-01-2006  90.3738
## 263 11-01-2006  92.3566
## 264 12-01-2006 103.0660
## 265 01-01-2007 112.0576
## 266 02-01-2007 111.8399
## 267 03-01-2007  99.1925
## 268 04-01-2007  90.8177
## 269 05-01-2007  92.0587
## 270 06-01-2007 100.9676
## 271 07-01-2007 107.5686
## 272 08-01-2007 114.1036
## 273 09-01-2007 101.5316
## 274 10-01-2007  93.0068
## 275 11-01-2007  93.9126
## 276 12-01-2007 106.7528
## 277 01-01-2008 114.8331
## 278 02-01-2008 108.2353
## 279 03-01-2008 100.4386
## 280 04-01-2008  90.9944
## 281 05-01-2008  91.2348
## 282 06-01-2008 103.9581
## 283 07-01-2008 110.7631
## 284 08-01-2008 107.5665
## 285 09-01-2008  97.7183
## 286 10-01-2008  90.9979
## 287 11-01-2008  93.8057
## 288 12-01-2008 109.4221
## 289 01-01-2009 116.8316
## 290 02-01-2009 104.4202
## 291 03-01-2009  97.8529
## 292 04-01-2009  88.1973
## 293 05-01-2009  87.5366
## 294 06-01-2009  97.2387
## 295 07-01-2009 103.9086
## 296 08-01-2009 105.7486
## 297 09-01-2009  94.8823
## 298 10-01-2009  89.2977
## 299 11-01-2009  89.3585
## 300 12-01-2009 110.6844
## 301 01-01-2010 119.0166
## 302 02-01-2010 110.5330
## 303 03-01-2010  98.2672
## 304 04-01-2010  86.3000
## 305 05-01-2010  90.8364
## 306 06-01-2010 104.3538
## 307 07-01-2010 112.8066
## 308 08-01-2010 112.9014
## 309 09-01-2010 100.1209
## 310 10-01-2010  88.9251
## 311 11-01-2010  92.7750
## 312 12-01-2010 114.3266
## 313 01-01-2011 119.4880
## 314 02-01-2011 107.3753
## 315 03-01-2011  99.1028
## 316 04-01-2011  89.3583
## 317 05-01-2011  90.0698
## 318 06-01-2011 102.8204
## 319 07-01-2011 114.7068
## 320 08-01-2011 113.5958
## 321 09-01-2011  99.4712
## 322 10-01-2011  90.3566
## 323 11-01-2011  93.8095
## 324 12-01-2011 107.3312
## 325 01-01-2012 111.9646
## 326 02-01-2012 103.3679
## 327 03-01-2012  93.5772
## 328 04-01-2012  87.5566
## 329 05-01-2012  92.7603
## 330 06-01-2012 101.1400
## 331 07-01-2012 113.0357
## 332 08-01-2012 109.8601
## 333 09-01-2012  96.7431
## 334 10-01-2012  90.3805
## 335 11-01-2012  94.3417
## 336 12-01-2012 105.2722
## 337 01-01-2013 115.5010
## 338 02-01-2013 106.7340
## 339 03-01-2013 102.9948
## 340 04-01-2013  91.0092
## 341 05-01-2013  90.9634
## 342 06-01-2013 100.6957
## 343 07-01-2013 110.1480
## 344 08-01-2013 108.1756
## 345 09-01-2013  99.2809
## 346 10-01-2013  91.7871
## 347 11-01-2013  97.2853
## 348 12-01-2013 113.4732
## 349 01-01-2014 124.2549
## 350 02-01-2014 112.8811
## 351 03-01-2014 104.7631
## 352 04-01-2014  90.2867
## 353 05-01-2014  92.1340
## 354 06-01-2014 101.8780
## 355 07-01-2014 108.5497
## 356 08-01-2014 108.1940
## 357 09-01-2014 100.4172
## 358 10-01-2014  92.3837
## 359 11-01-2014  99.7033
## 360 12-01-2014 109.3477
## 361 01-01-2015 120.2696
## 362 02-01-2015 116.3788
## 363 03-01-2015 104.4706
## 364 04-01-2015  89.7461
## 365 05-01-2015  91.0930
## 366 06-01-2015 102.6495
## 367 07-01-2015 111.6354
## 368 08-01-2015 110.5925
## 369 09-01-2015 101.9204
## 370 10-01-2015  91.5959
## 371 11-01-2015  93.0628
## 372 12-01-2015 103.2203
## 373 01-01-2016 117.0837
## 374 02-01-2016 106.6688
## 375 03-01-2016  95.3548
## 376 04-01-2016  89.3254
## 377 05-01-2016  90.7369
## 378 06-01-2016 104.0375
## 379 07-01-2016 114.5397
## 380 08-01-2016 115.5159
## 381 09-01-2016 102.7637
## 382 10-01-2016  91.4867
## 383 11-01-2016  92.8900
## 384 12-01-2016 112.7694
## 385 01-01-2017 114.8505
## 386 02-01-2017  99.4901
## 387 03-01-2017 101.0396
## 388 04-01-2017  88.3530
## 389 05-01-2017  92.0805
## 390 06-01-2017 102.1532
## 391 07-01-2017 112.1538
## 392 08-01-2017 108.9312
## 393 09-01-2017  98.6154
## 394 10-01-2017  93.6137
## 395 11-01-2017  97.3359
## 396 12-01-2017 114.7212
## 397 01-01-2018 129.4048

Subsetting Data

# Ambil indeks data ke 241 Hingga selesai
Datalistrik <- Data[241:397,]
head(Datalistrik)

##           DATE    Value
## 241 01-01-2005 111.1614
## 242 02-01-2005 101.7795
## 243 03-01-2005  98.9565
## 244 04-01-2005  86.4776
## 245 05-01-2005  87.2234
## 246 06-01-2005  99.5076

class(Datalistrik)

## [1] "data.frame"

#Mengubah menjadi Time series
Datalistrik <- ts(Datalistrik$Value, start = c(2005, 1), end = c(2018, 12), frequency = 12)

Visualisasi Data

ts.plot(Datalistrik, type="l", xlab = "DATE", ylab="Value", col="blue")
title(main = "Time Series Plot of CO2 (ppm)", cex.sub = 0.8)
points(co2, pch = 20, col = "blue")

Bila dilihat dari Pola datanya pola musiman yang kuat dengan fluktuasi yang berulang secara teratur setiap tahun. Terlihat bahwa setiap tahun memiliki puncak dan lembah yang konsisten pada periode waktu yang hampir sama, mengindikasikan adanya siklus tahunan yang berulang

Plot Decompose

dec.Listrik<- decompose(Datalistrik)
plot(dec.Listrik)

Secara eksplorasi, terlihat adanya kecenderungan data memiliki tren Stasioner (Kosntan) dan perilaku berulang kecenderungan musiman) dalam deret tersebut. Kecenderungan musiman dapat dilihat dengan lebih jelas dengan menampilkan deret waktu per tahun.

Plot Tahunan

seasonplot(Datalistrik,12,main="Seasonal Plot of Electric Consumption", ylab="Year",
           year.labels = TRUE, col=rainbow(18))

Pada Plot Seasonal, Dapat dilihat bahwa dari setiap Tahun terdapat pola data yang sama. Tingkat konsumsi listrik cenderung meningkat pada bulan-bulan tertentu dan menurun pada bulan-bulan tertentu. Dimulai dari Awal Tahun, konsumsi listrik cenderung menurun hingga ke Bulan April, lalu kemudian mulai meningkat lagi Dari Mei Ke September, dan kemudian Menurun lagi Hingga ke November. Konsumsi listrik kemudian kembali meningkat pada akhir tahun (Desember).

Plot Bulanan

monthplot(Datalistrik ,ylab="Value", col="blue")

Pada Plot Bulanan, Dapat dilihat bahwa dari setiap Bulan terdapat pola fluktuasi yang acak dan tidak teratur. Konsumsi listrik cenderung meningkat pada bulan-bulan tertentu dan menurun pada bulan-bulan tertentu. Dapat terlihat juga pada pertengahan Bulan terdapat penurunan penggunaan listrik yang cukup signifikan.

frame<-data.frame(values=as.matrix(Datalistrik), date=lubridate::year(zoo::as.Date(Datalistrik)))

library(ggplot2)
ggplot(frame,aes(y=values,x=date,group=date))+
  geom_boxplot()

Berdasarkan hasil plot di atas dapat terlihat bahwa data memiliki pola yang hampir sama dari tahun ke tahun sehingga dapat disimpulkan bahwa periode musimannya adalah 12. Selain itu, apabila dilihat dari boxplot, terlihat bahwa data cenderung homogen dari tahun ke tahun. Untuk memastikan bahwa data homogen akan dilakukan uji homogenitas dengan fligner.test.

Homogenitas Data

library(car)
fligner.test(values ~ date, data=frame)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  values by date
## Fligner-Killeen:med chi-squared = 7.0754, df = 13, p-value = 0.8982

Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi α=5% didapatkan p-value sebesar 0.8982. p−value=0.9885>α=0.05 sehingga tak tolak H0H0atau dengan kata lain ragam data sudah stasioner.

Splitting Data

# Bagi Data Latih dan Data Uji 70 * 30  dari data Listrik

length(Datalistrik)

## [1] 168

train.ts<- subset(Datalistrik,start=1,end=118)
test.ts<- subset(Datalistrik,start=119,end=168)

autoplot(train.ts) + theme_bw() + xlab("Year") + ylab("Electric Consumption")

autoplot(test.ts) + theme_bw() + xlab("Year") + ylab("Electric Consumption")

Bila dilihat secara Plot, Menunjukkan bahwa data sudah stasioner dalam Rataan. dan Data menunjukan pola musiman tanpa tren / konstan, Untuk Mendukung analisa ini, kemudian akan di cek dengan plot ACF nya.

Non Seasonal ARIMA

Stasioner Data

acf0 <- acf(train.ts,main="ACF",lag.max=48,xaxt="n", col="blue")
axis(1, at=0:48/12, labels=0:48)

Dari plot ACF terlihat adanya pola musiman yang jelas dengan periode 12 bulan, ditunjukkan oleh spike-spike signifikan yang muncul pada lag 12, 24, 36, dan 48. Pola ACF menunjukkan penurunan yang lambat (dies down slowly) dengan alternasi antara nilai positif dan negatif, di mana beberapa lag melewati batas signifikansi (ditandai dengan garis putus-putus).

acf0$lag <- acf0$lag * 12
acf0.1 <- as.data.frame(cbind(acf0$acf,acf0$lag))
acf0.2 <- acf0.1[which(acf0.1$V2%%12==0),]
barplot(height = acf0.2$V1, 
names.arg=acf0.2$V2, ylab="ACF", xlab="Lag")

tseries::adf.test(train.ts)

## Warning in tseries::adf.test(train.ts): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -4.2138, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

H0H0 : Data tidak stasioner dalam rataan

H1H1 : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.01 yang lebih kecil dari taraf nyata 5% sehingga tak tolak H0H0 dan menandakan bahwa data sudah stasioner dalam rataan.

Seasonal ARIMA

Pada Plot ACF Non- Seasonal , dapat dilihat bahwa kestasioneran komponen non-seasonal (namun perhatikan lag 12,24, dst), pada series seasonal belum stasioner. Hal ini menunjukkan adanya kecenderungan musiman

Diferencing Komponen Musiman

D1 <- diff(train.ts,12)
ts.plot(D1, type="l", ylab="D1 Xt", col="blue")

acf2<-acf(D1,lag.max=48,xaxt="n", main="ACF D1", col="blue")

acf2$lag <- acf2$lag * 12
acf2.1 <- as.data.frame(cbind(acf2$acf,acf2$lag))
acf2.2 <- acf2.1[which(acf2.1$V2%%12==0),]
barplot(height = acf2.2$V1, names.arg=acf2.2$V2, ylab="ACF", xlab="Lag")

Non-seasonal differencing D = 12 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen seasonalnya (namun tidak untuk komponen non-seasonalnya). Untuk menghilangkan kecenderungan musiman dilakukan pembedaan musiman terhadap deret hasil pembedaan pertama.

d1D1 <- diff(D1)
ts.plot(d1D1, type="l", ylab="d1 D1 Xt", col="blue")

Setelah pembedaan pertama dan pembedaan musiman tampak bahwa deret sudah tidak memiliki kecenderungan apapun. Selanjutnya penentuan ordo p, q dan P, Q dapat dilakukan menggunakan plot ACF dan PACF contoh dari deret hasil pembedaan pertama dan pembedaan musiman tersebut.

Identifikasi Model

acf3 <- acf(d1D1,lag.max=48,xaxt="n", main="ACF d1D1", col="blue")
axis(1, at=0:48/12, labels=0:48)

Berdasarkan plot ACF contoh lag 2 signifikan sehingga dipilih ordo q=2 , dan lag 12 adalah satu-satunya lag musiman yang signifikan sehingga order Q=1.

pacf3 <- pacf(d1D1,lag.max=48,xaxt="n", main="PACF d1D1", col="blue")
axis(1, at=0:48/12, labels=0:48)

Plot PACF contoh menunjukkan cuts-off pada lag-2 sehingga ordo p=2 sementara pada pola musimannya tidak terlihat model AR yang terbentuk karena tidak adanya lag yang signifikan pada pola musimannya.

Model musiman yang dipilih untuk deret konsentrasi karbon dioksida adalah

ARIMA(0,0,2)×(0,1,1)12

ARIMA(2,0,0)×(0,1,1)12

ARIMA(2,0,2)×(0,1,1)12.

Ingat kembali bahwa model yang digunakan bersifat tentatif dan dapat berubah saat diagnostik model.

TSA::eacf(d1D1)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o x o o o o o o o o o  x  o  o 
## 1 x x o o o o o o o o o  x  o  o 
## 2 o x o x o o o o o o o  o  o  o 
## 3 o o x o o o o o o o o  x  o  o 
## 4 x x x o o o o o o o o  x  o  o 
## 5 x x x x o o o o o o o  x  o  o 
## 6 x x o o o o o o o o o  o  o  o 
## 7 x o o o o o o x o o o  o  o  o

Dari Plot EACFnya akan dicoba pula Model Tentantif Non-Seasonalnya ARIMA (3,0,1) dan (1,0,2), Sehingga Model yang akan diuji adalah :

ARIMA ( 0,0,2) x (0,1,1)12 = Model 1

ARIMA (2,0,0) x (0,1,1)12 = Model 2

ARIMA (2,0,2) x (0,1,1)12= Model 3

ARIMA(3,0,1) x (0,1,1)12= Model 4

ARIMA (1,0,2) x (0,1,1)12= Model 5

Pendugaan Parameter

ARIMA (0,0,2)x(0,0,1)12

tmodel1 <- Arima(train.ts,order=c(0,0,2),seasonal=c(0,1,1))
summary(tmodel1)

## Series: train.ts 
## ARIMA(0,0,2)(0,1,1)[12] 
## 
## Coefficients:
##          ma1    ma2     sma1
##       0.5593  0.131  -0.7165
## s.e.  0.0950  0.082   0.1141
## 
## sigma^2 = 7.883:  log likelihood = -262.77
## AIC=533.55   AICc=533.94   BIC=544.2
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.4848395 2.623072 2.000508 0.4284446 1.933467 0.6761762
##                     ACF1
## Training set -0.01038441

lmtest::coeftest(tmodel1)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1   0.559263   0.094997  5.8872 3.928e-09 ***
## ma2   0.130975   0.081956  1.5981      0.11    
## sma1 -0.716474   0.114075 -6.2807 3.370e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (2,0,0)x(0,0,1)12

tmodel2 <- Arima(train.ts,order=c(2,0,0),seasonal=c(0,1,1))
summary(tmodel2)

## Series: train.ts 
## ARIMA(2,0,0)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2     sma1
##       0.5519  -0.0412  -0.7372
## s.e.  0.0966   0.0979   0.1135
## 
## sigma^2 = 7.704:  log likelihood = -261.95
## AIC=531.9   AICc=532.29   BIC=542.55
## 
## Training set error measures:
##                     ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 0.4158064 2.59326 1.956568 0.3617349 1.892038 0.6613245
##                     ACF1
## Training set -0.01151213

lmtest::coeftest(tmodel2)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.551858   0.096587  5.7136 1.106e-08 ***
## ar2  -0.041172   0.097936 -0.4204    0.6742    
## sma1 -0.737232   0.113478 -6.4967 8.209e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (2,0,2)x(0,0,1)12

tmodel3 <- Arima(train.ts,order=c(2,0,2),seasonal=c(0,1,1))
summary(tmodel3)

## Series: train.ts 
## ARIMA(2,0,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     ma1      ma2     sma1
##       0.3659  0.1303  0.1964  -0.1138  -0.7128
## s.e.  0.6864  0.3044  0.6838   0.1982   0.1177
## 
## sigma^2 = 7.865:  log likelihood = -261.57
## AIC=535.14   AICc=535.99   BIC=551.13
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.3863274 2.594652 1.967072 0.3318968 1.901874 0.6648749
##                     ACF1
## Training set -0.01375998

lmtest::coeftest(tmodel3)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.36586    0.68640  0.5330    0.5940    
## ar2   0.13032    0.30441  0.4281    0.6686    
## ma1   0.19641    0.68377  0.2872    0.7739    
## ma2  -0.11378    0.19822 -0.5740    0.5660    
## sma1 -0.71279    0.11772 -6.0550 1.404e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (3,0,1)x(0,0,1)12

tmodel4 <- Arima(train.ts,order=c(3,0,1),seasonal=c(0,1,1))
summary(tmodel4)

## Series: train.ts 
## ARIMA(3,0,1)(0,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3     ma1     sma1
##       0.2346  0.0781  0.0830  0.3249  -0.7044
## s.e.  0.4140  0.2487  0.1047  0.4070   0.1191
## 
## sigma^2 = 7.867:  log likelihood = -261.44
## AIC=534.89   AICc=535.74   BIC=550.87
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.3758059 2.594845 1.969906 0.3211041 1.904784 0.6658327
##                      ACF1
## Training set -0.007181832

lmtest::coeftest(tmodel4)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.234639   0.414038  0.5667    0.5709    
## ar2   0.078123   0.248658  0.3142    0.7534    
## ar3   0.083006   0.104710  0.7927    0.4279    
## ma1   0.324913   0.406995  0.7983    0.4247    
## sma1 -0.704393   0.119055 -5.9165 3.288e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (1,0,2)x(0,0,1)12

tmodel5 <- Arima(train.ts,order=c(1,0,2),seasonal=c(0,1,1))
summary(tmodel5)

## Series: train.ts 
## ARIMA(1,0,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1      ma2     sma1
##       0.6376  -0.0721  -0.1403  -0.7195
## s.e.  0.2967   0.3191   0.2089   0.1160
## 
## sigma^2 = 7.78:  log likelihood = -261.63
## AIC=533.27   AICc=533.87   BIC=546.58
## 
## Training set error measures:
##                     ME     RMSE     MAE       MPE     MAPE      MASE       ACF1
## Training set 0.3849569 2.593315 1.96731 0.3303347 1.902007 0.6649552 -0.0157086

lmtest::coeftest(tmodel5)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.637556   0.296741  2.1485   0.03167 *  
## ma1  -0.072138   0.319056 -0.2261   0.82113    
## ma2  -0.140293   0.208856 -0.6717   0.50176    
## sma1 -0.719539   0.115969 -6.2046 5.484e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Auto-Arima

model.auto.arima <- auto.arima(train.ts)
summary(model.auto.arima)

## Series: train.ts 
## ARIMA(1,0,0)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     sma1   drift
##       0.4856  -0.8122  0.0306
## s.e.  0.0843   0.1344  0.0148
## 
## sigma^2 = 7.24:  log likelihood = -260.29
## AIC=528.59   AICc=528.98   BIC=539.24
## 
## Training set error measures:
##                      ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.02309323 2.513855 1.906158 -0.03193381 1.844289 0.6442858
##                   ACF1
## Training set 0.0439364

lmtest::coeftest(model.auto.arima)

## 
## z test of coefficients:
## 
##        Estimate Std. Error z value  Pr(>|z|)    
## ar1    0.485571   0.084255  5.7631 8.258e-09 ***
## sma1  -0.812187   0.134434 -6.0415 1.527e-09 ***
## drift  0.030588   0.014776  2.0701   0.03844 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model kemudian dicek dengan menggunakan fungsi Auto-Arima didapatkan Model ARIMA (1,0,0) (0,1,1) 12

Matrix Pendugaan Parameter

AICKandidatModel <- c(tmodel1$aic, tmodel2$aic, tmodel3$aic,
                      tmodel4$aic, tmodel5$aic, model.auto.arima$aic)
AICcKandidatModel <- c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc,
                       tmodel4$aicc, tmodel5$aicc,model.auto.arima$aicc)
BICKandidatModel <- c(tmodel1$bic, tmodel2$bic, tmodel3$bic,
                      tmodel4$bic, tmodel5$bic,model.auto.arima$bic)
KandidatModelARIMA <- c("ARIMA(0,0,2)(0,1,1)12", "ARIMA(2,0,0)(0,1,1)12",
                        "ARIMA(2,0,2)(0,1,1)12", "ARIMA(3,0,1)(0,1,1)12",
                        "ARIMA(1,0,2)(0,1,1)12","Model Auto Arima")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel,
                        AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", 
                              "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA

##          Kandidat Model        Nilai AIC       Nilai AICc        Nilai BIC
## 1 ARIMA(0,0,2)(0,1,1)12 533.548304634723 533.944344238683 544.202061011171
## 2 ARIMA(2,0,0)(0,1,1)12  531.89776647055  532.29380607451 542.551522846998
## 3 ARIMA(2,0,2)(0,1,1)12 535.144915323903 535.993400172388 551.125549888576
## 4 ARIMA(3,0,1)(0,1,1)12 534.889154243281 535.737639091766 550.869788807953
## 5 ARIMA(1,0,2)(0,1,1)12 533.266271446165 533.866271446165 546.583466916725
## 6      Model Auto Arima 528.588177845735 528.984217449696 539.241934222183

Model Terbaik berdasarkan Nilai AIC , AICc, dan BIC dari kandidat model yaitu Model Auto ARIMA, dimana mendapatkan nilai AIC , AICC dan BIC yang paling kecil dari semua perbandingan model. Namun Model ini adalah model yang bertentangan dengan hasil eksplorasi dimana baik dari ACF hingga EACF, sehingga untuk selanjutnya, model yang tetap digunakan adalah model ARIMA (2,0,0)(0,1,1)12. Secara keseluruhan Tidak ada model yang seluruh Parameternya signifikan, Namun Model 1 ARIMA (0,0,2)(0,1,1)12dan Model 2 ARIMA (2,0,0) (0,1,1)12 adalah model yang paling banyak memiliki parameter yang signifikan sehingga setelah dibandingkan dengan metrik diatas, terpilih model 2.

Diagnostik Model

tsdisplay(residuals(tmodel2), lag.max=45, 
          main='ARIMA(2,0,0)(0,1,1)12 Model Residuals', col="blue")

#Eksplorasi
sisaan.model2 <- tmodel2$residuals
par(mfrow=c(2,2))
car::qqPlot(sisaan.model2)

## [1] 24 26

plot(c(1:length(sisaan.model2)),sisaan.model2)
acf(sisaan.model2)
pacf(sisaan.model2)

Berdasarkan plot di atas terlihat bahwa sisaan tidak mengikuti sebaran normal karena ada beberapa data yang keluar dari treshold . Selanjutnya, ditinjau dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa kemungkinan tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

Uji Formal

Normalitas

ks.test(sisaan.model2,"pnorm")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.model2
## D = 0.25761, p-value = 3.157e-07
## alternative hypothesis: two-sided

shapiro.test(sisaan.model2)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan.model2
## W = 0.98693, p-value = 0.3151

nortest::ad.test(sisaan.model2)

## 
##  Anderson-Darling normality test
## 
## data:  sisaan.model2
## A = 0.64162, p-value = 0.09205

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS), Shapiro-Wilk, dan Anderson-Darling. Hipotesis pada uji kenormalan adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 0.000 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Namun, hasil ini berbeda dengan Shapiro-Wilk test dan Anderson-Darling test yang menunjukkan p-value > 5%.

Autokorelasi

Box.test(sisaan.model2, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan.model2
## X-squared = 0.016039, df = 1, p-value = 0.8992

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.8992 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini berbeda dengan eksplorasi.

Homogen

Box.test((sisaan.model2)^2, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  (sisaan.model2)^2
## X-squared = 10.292, df = 1, p-value = 0.001336

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.00013 yang lebih kecil dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak homogen.

Nilai Tengah Sisaan 0

t.test(sisaan.model2, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  sisaan.model2
## t = 1.7571, df = 117, p-value = 0.08152
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.05285675  0.88446951
## sample estimates:
## mean of x 
## 0.4158064

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.08 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini berbeda dengan eksplorasi.

Untuk masalah Kehomogenan ragam, bisa diatasi dengan Transformasi Box-Cox , sehingga Ragam dapat homogen, bisa pula memodelkan sisaanya dengan Metode ARCH/GARCH, sehingga model yang dihasilkan adalah Model Gabungan.

Overfitting

Overfit Non-Seasonal

tmodel1.ofq <- Arima(train.ts,order=c(3,0,0),seasonal=c(0,1,1))
summary(tmodel1.ofq)

## Series: train.ts 
## ARIMA(3,0,0)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2     ar3     sma1
##       0.5552  -0.0872  0.0791  -0.7213
## s.e.  0.0964   0.1137  0.0996   0.1152
## 
## sigma^2 = 7.776:  log likelihood = -261.63
## AIC=533.27   AICc=533.87   BIC=546.59
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.3800029 2.592571 1.968747 0.3247123 1.903638 0.665441
##                      ACF1
## Training set -0.002927216

lmtest::coeftest(tmodel1.ofq)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.555188   0.096437  5.7570 8.561e-09 ***
## ar2  -0.087245   0.113651 -0.7677    0.4427    
## ar3   0.079149   0.099617  0.7945    0.4269    
## sma1 -0.721286   0.115246 -6.2587 3.883e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pada model musiman, ordo yang dilakukan overfit adalah ordo musiman (P, Q).

tmodel1.ofP <- Arima(train.ts,order=c(2,0,0),seasonal=c(1,1,1))
summary(tmodel1.ofP)

## Series: train.ts 
## ARIMA(2,0,0)(1,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2    sar1     sma1
##       0.5545  -0.0344  0.2225  -0.9022
## s.e.  0.0967   0.0980  0.1474   0.2457
## 
## sigma^2 = 7.286:  log likelihood = -260.74
## AIC=531.48   AICc=532.08   BIC=544.79
## 
## Training set error measures:
##                     ME    RMSE     MAE       MPE     MAPE      MASE        ACF1
## Training set 0.4051204 2.50951 1.89434 0.3528974 1.834228 0.6402911 -0.01641311

lmtest::coeftest(tmodel1.ofP)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.554491   0.096674  5.7357 9.711e-09 ***
## ar2  -0.034359   0.097958 -0.3508 0.7257737    
## sar1  0.222464   0.147373  1.5095 0.1311623    
## sma1 -0.902202   0.245682 -3.6722 0.0002404 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tmodel1.ofQ <- Arima(train.ts,order=c(2,0,0),seasonal=c(0,1,2))
summary(tmodel1.ofQ)

## Series: train.ts 
## ARIMA(2,0,0)(0,1,2)[12] 
## 
## Coefficients:
##          ar1      ar2     sma1     sma2
##       0.5513  -0.0426  -0.6456  -0.2586
## s.e.  0.0966   0.0979   0.2457   0.1681
## 
## sigma^2 = 7.059:  log likelihood = -260.1
## AIC=530.2   AICc=530.8   BIC=543.51
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.4036941 2.470266 1.857008 0.3527502 1.799395 0.6276727
##                     ACF1
## Training set -0.01969153

lmtest::coeftest(tmodel1.ofQ)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.551270   0.096589  5.7074 1.147e-08 ***
## ar2  -0.042618   0.097878 -0.4354  0.663258    
## sma1 -0.645566   0.245659 -2.6279  0.008592 ** 
## sma2 -0.258621   0.168116 -1.5384  0.123962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model overfitting yang dicobakan menghasilkan nilai AIC dan signifikansi parameter yang tidak lebih baik dari model awal. Dimana parameter-parameter overfit tidak signifikan. Maka karena itu, model yang digunakan tetap model awal.

Forecasting

ramalan_sarima = forecast::forecast(tmodel2, 50)
ramalan_sarima

##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## Nov 2014       95.02788  91.46742  98.58835  89.58262 100.47314
## Dec 2014      109.94648 105.88022 114.01275 103.72767 116.16530
## Jan 2015      117.64080 113.46913 121.81247 111.26078 124.02082
## Feb 2015      108.14415 103.94973 112.33857 101.72933 114.55896
## Mar 2015      100.49586  96.29657 104.69515  94.07360 106.91812
## Apr 2015       89.35557  85.15524  93.55590  82.93172  95.77942
## May 2015       91.13972  86.93917  95.34028  84.71553  97.56392
## Jun 2015      101.51229  97.31169 105.71290  95.08802 107.93656
## Jul 2015      110.37246 106.17185 114.57308 103.94818 116.79675
## Aug 2015      109.58197 105.38135 113.78259 103.15768 116.00626
## Sep 2015       98.92617  94.72556 103.12679  92.50189 105.35046
## Oct 2015       91.10612  86.90552  95.30672  84.68185  97.53039
## Nov 2015       94.38423  90.07960  98.68885  87.80087 100.96758
## Dec 2015      109.64388 105.30843 113.97932 103.01338 116.27437
## Jan 2016      117.50030 113.15905 121.84156 110.86093 124.13968
## Feb 2016      108.07907 103.73629 112.42185 101.43737 114.72078
## Mar 2016      100.46573  96.12263 104.80884  93.82353 107.10794
## Apr 2016       89.34162  84.99845  93.68480  82.69931  95.98394
## May 2016       91.13327  86.79008  95.47646  84.49093  97.77560
## Jun 2016      101.50930  97.16611 105.85250  94.86696 108.15164
## Jul 2016      110.37108 106.02789 114.71427 103.72874 117.01342
## Aug 2016      109.58133 105.23814 113.92452 102.93899 116.22367
## Sep 2016       98.92588  94.58268 103.26907  92.28354 105.56821
## Oct 2016       91.10598  86.76280  95.44916  84.46366  97.74830
## Nov 2016       94.38416  89.94036  98.82796  87.58796 101.18037
## Dec 2016      109.64385 105.17020 114.11749 102.80200 116.48570
## Jan 2017      117.50029 113.02102 121.97956 110.64983 124.35075
## Feb 2017      108.07907 103.59832 112.55981 101.22635 114.93178
## Mar 2017      100.46573  95.98467 104.94679  93.61253 107.31893
## Apr 2017       89.34162  84.86049  93.82275  82.48832  96.19492
## May 2017       91.13326  86.65212  95.61441  84.27994  97.98659
## Jun 2017      101.50930  97.02815 105.99045  94.65598 108.36263
## Jul 2017      110.37108 105.88993 114.85223 103.51775 117.22441
## Aug 2017      109.58133 105.10018 114.06248 102.72800 116.43466
## Sep 2017       98.92588  94.44473 103.40702  92.07255 105.77920
## Oct 2017       91.10598  86.62484  95.58712  84.25267  97.95929
## Nov 2017       94.38416  89.80544  98.96289  87.38161 101.38672
## Dec 2017      109.64385 105.03615 114.25154 102.59698 116.69071
## Jan 2018      117.50029 112.88713 122.11345 110.44507 124.55551
## Feb 2018      108.07907 103.46447 112.69366 101.02165 115.13648
## Mar 2018      100.46573  95.85083 105.08063  93.40785 107.52361
## Apr 2018       89.34162  84.72666  93.95659  82.28364  96.39960
## May 2018       91.13326  86.51828  95.74824  84.07526  98.19127
## Jun 2018      101.50930  96.89432 106.12429  94.45130 108.56731
## Jul 2018      110.37108 105.75610 114.98606 103.31307 117.42909
## Aug 2018      109.58133 104.96635 114.19631 102.52332 116.63934
## Sep 2018       98.92588  94.31089 103.54086  91.86787 105.98388
## Oct 2018       91.10598  86.49101  95.72095  84.04799  98.16397
## Nov 2018       94.38416  89.67437  99.09395  87.18116 101.58716
## Dec 2018      109.64385 104.90589 114.38180 102.39777 116.88993

autoplot(ramalan_sarima, col="blue")

accuracy(ramalan_sarima,test.ts)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.4158064 2.593260 1.956568 0.3617349 1.892038 0.6613245
## Test set     0.4094448 5.598299 3.820322 0.1422125 3.761979 1.2912774
##                     ACF1 Theil's U
## Training set -0.01151213        NA
## Test set      0.36811045 0.5961036

Didapatkan MAPE yang tidak terlau jauh berbeda yang mengindikasikan bahwa tidak terjadi Overfitting. Nilai MPAE yang kecil ini ini ( Under 10%) mengindikasikan bahwa model dapat meramalkan data dengan akurasi yang sangat baik.