library("readxl")
library("tseries")
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library("forecast")
## Warning: package 'forecast' was built under R version 4.5.1
library("TTR")
library("TSA")
## Warning: package 'TSA' was built under R version 4.5.2
## 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("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
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(psych)
## Warning: package 'psych' was built under R version 4.5.3
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
suntik <- read_excel("D:/MAGANG BKKBN/DATA FIKS/SUNTIK.xlsx")
data.table::data.table(suntik)
## PERIODE SUNTIK
## <char> <num>
## 1: 2020M1 988
## 2: 2020M2 1189
## 3: 2020M3 1011
## 4: 2020M4 1031
## 5: 2020M5 1021
## 6: 2020M6 1653
## 7: 2020M7 1102
## 8: 2020M8 1180
## 9: 2020M9 1071
## 10: 2020M10 999
## 11: 2020M11 1310
## 12: 2020M12 929
## 13: 2021M1 894
## 14: 2021M2 884
## 15: 2021M3 963
## 16: 2021M4 821
## 17: 2021M5 1015
## 18: 2021M6 1030
## 19: 2021M7 859
## 20: 2021M8 1009
## 21: 2021M9 1168
## 22: 2021M10 1030
## 23: 2021M11 957
## 24: 2021M12 981
## 25: 2022M1 456
## 26: 2022M2 413
## 27: 2022M3 414
## 28: 2022M4 362
## 29: 2022M5 413
## 30: 2022M6 294
## 31: 2022M7 253
## 32: 2022M8 301
## 33: 2022M9 271
## 34: 2022M10 324
## 35: 2022M11 276
## 36: 2022M12 358
## 37: 2023M1 680
## 38: 2023M2 571
## 39: 2023M3 544
## 40: 2023M4 404
## 41: 2023M5 533
## 42: 2023M6 559
## 43: 2023M7 573
## 44: 2023M8 543
## 45: 2023M9 643
## 46: 2023M10 579
## 47: 2023M11 575
## 48: 2023M12 601
## 49: 2024M1 586
## 50: 2024M2 635
## 51: 2024M3 551
## 52: 2024M4 577
## 53: 2024M5 562
## 54: 2024M6 626
## 55: 2024M7 669
## 56: 2024M8 520
## 57: 2024M9 473
## 58: 2024M10 526
## 59: 2024M11 509
## 60: 2024M12 636
## 61: 2025M1 603
## 62: 2025M2 478
## 63: 2025M3 437
## 64: 2025M4 447
## 65: 2025M5 427
## 66: 2025M6 453
## 67: 2025M7 382
## 68: 2025M8 376
## 69: 2025M9 440
## 70: 2025M10 261
## 71: 2025M11 311
## 72: 2025M12 240
## 73: 2026M1 211
## 74: 2026M2 362
## 75: 2026M3 225
## PERIODE SUNTIK
suntik.ts<-ts(suntik$SUNTIK)
suntik.ts.month<-ts(suntik$SUNTIK,start = c(2020,1),frequency = 12)
#Statistik Deskriptif
describe(suntik$SUNTIK)
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 75 647.44 310.42 571 625.38 309.86 211 1653 1442 0.76 -0.08
## se
## X1 35.84
autoplot(suntik.ts) +
ggtitle("Plot Time Series suntik (2020–2026)") +
xlab("Bulan") +
ylab("Jumlah suntik") +
theme_minimal()
#Partisi Data (Train-Test)
suntik.training<-subset(suntik.ts, start = 1, end =52 )
suntik.training
## Time Series:
## Start = 1
## End = 52
## Frequency = 1
## [1] 988 1189 1011 1031 1021 1653 1102 1180 1071 999 1310 929 894 884 963
## [16] 821 1015 1030 859 1009 1168 1030 957 981 456 413 414 362 413 294
## [31] 253 301 271 324 276 358 680 571 544 404 533 559 573 543 643
## [46] 579 575 601 586 635 551 577
suntik.test<-subset(suntik.ts, start = 53, end = 75)
suntik.test
## Time Series:
## Start = 53
## End = 75
## Frequency = 1
## [1] 562 626 669 520 473 526 509 636 603 478 437 447 427 453 382 376 440 261 311
## [20] 240 211 362 225
# data training
plot.ts(suntik.training, main="Data Training suntik", xlab="Bulan", ylab="Jumlah suntik")
# data testing
plot.ts(suntik.test, main="Data Testing suntik", xlab="Bulan", ylab="Jumlah suntik")
#Stasioneritas Data
acf(suntik.training,lag.max=20, main="Plot ACF Data suntik")
adf.test(suntik.training) # Augmented Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: suntik.training
## Dickey-Fuller = -1.633, Lag order = 3, p-value = 0.722
## alternative hypothesis: stationary
# Differencing
# Differencing Ordo 1
suntik.dif1 <- diff(suntik.training, difference=1)
plot.ts(suntik.dif1, lty=1, xlab="Waktu", ylab="Data Ekspor Diff Ordo 1")
points(suntik.dif1)
# Cek Kestasioneran
acf(suntik.dif1, main="Plot ACF Data suntik setelah Diff 1", lag.max = 20)
adf.test(suntik.dif1) # Augmented Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: suntik.dif1
## Dickey-Fuller = -3.8429, Lag order = 3, p-value = 0.02316
## alternative hypothesis: stationary
#Plot ACF
acf(suntik.dif1, main="Plot ACF Data suntik", lag.max=20)
#Plot PACF
pacf(suntik.dif1, main="Plot PACF Data suntik", lag.max=20)
#EACF
eacf(suntik.dif1)
## 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 o o o o o o o o o o o o o
## 2 o o o o o o o o o o o o o o
## 3 x o o o o o o o o o o o o o
## 4 x x x o o o o o o o o o o o
## 5 x o o o o o o o o o o o o o
## 6 o o o o o o o o o o o o o o
## 7 o o o o o o o o o o o o o o
#Pendugaan Parameter dan Penentuan Model Terbaik
arima110<-arima(suntik.dif1, order=c(1,1,0),include.mean=FALSE, method="ML") #ARIMA (1,1,0)
arima011<-arima(suntik.dif1, order=c(0,1,1),include.mean=FALSE, method="ML") #ARIMA (0,1,1)
arima111<-arima(suntik.dif1, order=c(1,1,1),include.mean=FALSE, method="ML") #ARIMA (1,1,1)
#Uji signifikansi
arima110
##
## Call:
## arima(x = suntik.dif1, order = c(1, 1, 0), include.mean = FALSE, method = "ML")
##
## Coefficients:
## ar1
## -0.6499
## s.e. 0.1065
##
## sigma^2 estimated as 53796: log likelihood = -343.54, aic = 689.09
lmtest::coeftest(arima110)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.64985 0.10650 -6.1019 1.048e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
arima011
##
## Call:
## arima(x = suntik.dif1, order = c(0, 1, 1), include.mean = FALSE, method = "ML")
##
## Coefficients:
## ma1
## -1.0000
## s.e. 0.0546
##
## sigma^2 estimated as 33960: log likelihood = -333.74, aic = 669.47
lmtest::coeftest(arima011)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.999999 0.054558 -18.329 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
arima111
##
## Call:
## arima(x = suntik.dif1, order = c(1, 1, 1), include.mean = FALSE, method = "ML")
##
## Coefficients:
## ar1 ma1
## -0.3816 -1.0000
## s.e. 0.1306 0.0658
##
## sigma^2 estimated as 28642: log likelihood = -329.88, aic = 663.75
lmtest::coeftest(arima111)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.38159 0.13059 -2.9221 0.003477 **
## ma1 -1.00000 0.06578 -15.2021 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#AIC
aic.arima<-data.frame(
"Model"=c("ARIMA(1,1,0)", "ARIMA(0,1,1)", "ARIMA(1,1,1)"),
"AIC"=c(arima110$aic, arima011$aic, arima111$aic))
aic.arima
## Model AIC
## 1 ARIMA(1,1,0) 689.0898
## 2 ARIMA(0,1,1) 669.4728
## 3 ARIMA(1,1,1) 663.7502
auto.arima(suntik.training, trace = TRUE)
##
## ARIMA(2,1,2) with drift : 680.1253
## ARIMA(0,1,0) with drift : 680.0519
## ARIMA(1,1,0) with drift : 673.775
## ARIMA(0,1,1) with drift : 672.9251
## ARIMA(0,1,0) : 677.9829
## ARIMA(1,1,1) with drift : 675.1286
## ARIMA(0,1,2) with drift : 675.1175
## ARIMA(1,1,2) with drift : Inf
## ARIMA(0,1,1) : 671.143
## ARIMA(1,1,1) : 673.2111
## ARIMA(0,1,2) : 673.1881
## ARIMA(1,1,0) : 671.8217
## ARIMA(1,1,2) : Inf
##
## Best model: ARIMA(0,1,1)
## Series: suntik.training
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.4135
## s.e. 0.1188
##
## sigma^2 = 28412: log likelihood = -333.45
## AIC=670.89 AICc=671.14 BIC=674.76
arima.111y<-arima(suntik.training, order=c(1,1,1), include.mean = FALSE,method="ML") #ARIMA (1,1,1)
sisaan111 <- arima.111y$residuals
# Eksplorasi
par(mfrow=c(2,2))
qqnorm(sisaan111)
qqline(sisaan111, col = "blue", lwd = 2)
plot(c(1:length(sisaan111)),sisaan111)
acf(sisaan111,lag.max = 10)
pacf(sisaan111,lag.max = 10)
par(mfrow=c(1,1))
#Uji Non-Autokorelasi dengan Ljung-Box test
Box.test(sisaan111, type = "Ljung")
##
## Box-Ljung test
##
## data: sisaan111
## X-squared = 0.0036575, df = 1, p-value = 0.9518
#Uji Normalitas dengan Kolmogorov-Smirnov test
shapiro.test(sisaan111)
##
## Shapiro-Wilk normality test
##
## data: sisaan111
## W = 0.90075, p-value = 0.0003905
ks.test(sisaan111,"pnorm")
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan111
## D = 0.59615, p-value < 2.2e-16
## alternative hypothesis: two-sided
#plot dugaan dengan data asli
dugaan <- fitted(arima.111y)
cbind(suntik.training,dugaan)
## Time Series:
## Start = 1
## End = 52
## Frequency = 1
## suntik.training dugaan
## 1 988 987.0120
## 2 1189 1005.4067
## 3 1011 1109.6668
## 4 1031 1065.1094
## 5 1021 1038.6846
## 6 1653 1027.7653
## 7 1102 1376.0292
## 8 1180 1260.6962
## 9 1071 1194.0229
## 10 999 1123.4074
## 11 1310 1046.7928
## 12 929 1187.2832
## 13 894 1059.7433
## 14 884 949.3782
## 15 963 905.3260
## 16 821 934.6326
## 17 1015 875.0326
## 18 1030 945.8460
## 19 859 1002.2597
## 20 1009 926.0281
## 21 1168 963.2411
## 22 1030 1083.8256
## 23 957 1065.2248
## 24 981 999.9872
## 25 456 983.5307
## 26 413 688.5684
## 27 414 503.0047
## 28 362 441.0433
## 29 413 393.2170
## 30 294 400.0155
## 31 253 342.5752
## 32 301 285.9356
## 33 271 289.8648
## 34 324 280.8451
## 35 276 303.6061
## 36 358 290.9650
## 37 680 326.3656
## 38 571 528.1718
## 39 544 572.7613
## 40 404 556.4588
## 41 533 469.6167
## 42 559 496.0819
## 43 573 536.2469
## 44 543 559.8707
## 45 643 552.2361
## 46 579 601.6689
## 47 575 594.6357
## 48 601 581.5407
## 49 586 591.5179
## 50 635 589.7272
## 51 551 614.5039
## 52 577 581.8284
plot.ts(suntik.training, xlab="Month", ylab="Data")
points(suntik.training)
par(col="red")
lines(dugaan)
par(col="black")
#forecast
ramalan_arima111<- forecast::forecast(suntik.training,model=arima.111y, h=23)
ramalan_arima111
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 53 574.9347 361.470031 788.3993 248.468671 901.4007
## 54 575.2159 330.716818 819.7149 201.286842 949.1449
## 55 575.1776 297.299043 853.0561 150.199029 1000.1561
## 56 575.1828 268.270954 882.0946 105.801640 1044.5639
## 57 575.1821 241.662506 908.7017 65.107906 1085.2563
## 58 575.1822 217.037953 933.3264 27.447851 1122.9165
## 59 575.1822 193.999249 956.3651 -7.786799 1158.1511
## 60 575.1822 172.275989 978.0884 -41.009660 1191.3740
## 61 575.1822 151.665482 998.6989 -72.530711 1222.8951
## 62 575.1822 132.012477 1018.3519 -102.587387 1252.9517
## 63 575.1822 113.194758 1037.1696 -131.366605 1281.7309
## 64 575.1822 95.114091 1055.2502 -159.018599 1309.3829
## 65 575.1822 77.690108 1072.6742 -185.666282 1336.0306
## 66 575.1822 60.856064 1089.5083 -211.411731 1361.7761
## 67 575.1822 44.555808 1105.8085 -236.340820 1386.7052
## 68 575.1822 28.741572 1121.6228 -260.526608 1410.8909
## 69 575.1822 13.372310 1136.9920 -284.031866 1434.3962
## 70 575.1822 -1.587551 1151.9519 -306.910998 1457.2753
## 71 575.1822 -16.169081 1166.5334 -329.211525 1479.5759
## 72 575.1822 -30.399612 1180.7640 -350.975243 1501.3396
## 73 575.1822 -44.303331 1194.6677 -372.239147 1522.6035
## 74 575.1822 -57.901772 1208.2661 -393.036168 1543.4005
## 75 575.1822 -71.214202 1221.5785 -413.395772 1563.7601
plot(ramalan_arima111)
gabungan<- cbind(suntik.test,ramalan_arima111)
df.gab<- as.data.frame(gabungan, digits=3)
df.gab %>%
rename(
"Data Aktual suntik"=suntik.test,
Ramalan="ramalan_arima111.Point Forecast",
"Lo 80"="ramalan_arima111.Lo 80",
"Hi 80"="ramalan_arima111.Hi 80",
"Lo 95"="ramalan_arima111.Lo 95",
"Hi 95"="ramalan_arima111.Hi 95"
)
## Data Aktual suntik Ramalan Lo 80 Hi 80 Lo 95 Hi 95
## 1 562 574.9347 361.470031 788.3993 248.468671 901.4007
## 2 626 575.2159 330.716818 819.7149 201.286842 949.1449
## 3 669 575.1776 297.299043 853.0561 150.199029 1000.1561
## 4 520 575.1828 268.270954 882.0946 105.801640 1044.5639
## 5 473 575.1821 241.662506 908.7017 65.107906 1085.2563
## 6 526 575.1822 217.037953 933.3264 27.447851 1122.9165
## 7 509 575.1822 193.999249 956.3651 -7.786799 1158.1511
## 8 636 575.1822 172.275989 978.0884 -41.009660 1191.3740
## 9 603 575.1822 151.665482 998.6989 -72.530711 1222.8951
## 10 478 575.1822 132.012477 1018.3519 -102.587387 1252.9517
## 11 437 575.1822 113.194758 1037.1696 -131.366605 1281.7309
## 12 447 575.1822 95.114091 1055.2502 -159.018599 1309.3829
## 13 427 575.1822 77.690108 1072.6742 -185.666282 1336.0306
## 14 453 575.1822 60.856064 1089.5083 -211.411731 1361.7761
## 15 382 575.1822 44.555808 1105.8085 -236.340820 1386.7052
## 16 376 575.1822 28.741572 1121.6228 -260.526608 1410.8909
## 17 440 575.1822 13.372310 1136.9920 -284.031866 1434.3962
## 18 261 575.1822 -1.587551 1151.9519 -306.910998 1457.2753
## 19 311 575.1822 -16.169081 1166.5334 -329.211525 1479.5759
## 20 240 575.1822 -30.399612 1180.7640 -350.975243 1501.3396
## 21 211 575.1822 -44.303331 1194.6677 -372.239147 1522.6035
## 22 362 575.1822 -57.901772 1208.2661 -393.036168 1543.4005
## 23 225 575.1822 -71.214202 1221.5785 -413.395772 1563.7601
accuracy(ramalan_arima111, suntik.test) #arima
## ME RMSE MAE MPE MAPE MASE
## Training set -13.55395 164.9580 109.8078 -5.222276 16.08957 0.9347684
## Test set -132.82487 185.9152 153.1068 -44.749309 47.90692 1.3033629
## ACF1 Theil's U
## Training set -0.008150368 NA
## Test set 0.714925200 2.601475
#forecast
(ramalan_arima111<- forecast::forecast(suntik.ts,model=arima.111y, h=9))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 76 267.2134 53.74873 480.6780 -59.25263 593.6794
## 77 261.4662 16.96720 505.9653 -112.46278 635.3953
## 78 262.2487 -15.62985 540.1272 -162.72987 687.2272
## 79 262.1422 -44.76968 569.0540 -207.23900 731.5233
## 80 262.1567 -71.36292 595.6762 -247.91752 772.2308
## 81 262.1547 -95.98954 620.2989 -285.57964 809.8890
## 82 262.1550 -119.02796 643.3379 -320.81401 845.1239
## 83 262.1549 -140.75126 665.0611 -354.03691 878.3468
## 84 262.1549 -161.36176 685.6716 -385.55796 909.8678
plot(ramalan_arima111)
