BISMILLAH SIC

Data dan Analisis Eksplorasi

Metodologi

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(readxl)
library(forecast)

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

library(TTR)
library(graphics)
require(smooth)

## Loading required package: smooth

## Warning: package 'smooth' was built under R version 4.2.3

## Loading required package: greybox

## Warning: package 'greybox' was built under R version 4.2.3

## Package "greybox", v1.0.8 loaded.

## This is package "smooth", v3.2.1

## 
## Attaching package: 'smooth'

## The following object is masked from 'package:TTR':
## 
##     lags

require(Mcomp)

## Loading required package: Mcomp

## Warning: package 'Mcomp' was built under R version 4.2.3

library(tseries)
library(readr)
library(TSA)

## 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(ggplot2)
library(cowplot)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(gtable)
library(grid)
library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

Input Data

#Input Data
data<- read_xlsx("C:/SEM6/SIC/Jumlah Korban Mati Kecelakaan 1992-2021.xlsx")
head(data)

## # A tibble: 6 × 2
##   Tahun               `Korban Mati`
##   <dttm>                      <dbl>
## 1 1992-01-01 00:00:00          9819
## 2 1993-01-01 00:00:00         10038
## 3 1994-01-01 00:00:00         11004
## 4 1995-01-01 00:00:00         10990
## 5 1996-01-01 00:00:00         10869
## 6 1997-01-01 00:00:00         12308

str(data)

## tibble [30 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Tahun      : POSIXct[1:30], format: "1992-01-01" "1993-01-01" ...
##  $ Korban Mati: num [1:30] 9819 10038 11004 10990 10869 ...

Yt <- data$`Korban Mati`

Eksplorasi

plot(x = data$Tahun,
     y = Yt,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Korban Mati",
     main = "Time Series Plot of Harga Gula Dunia")

summary(data)

##      Tahun                      Korban Mati   
##  Min.   :1992-01-01 00:00:00   Min.   : 8762  
##  1st Qu.:1999-04-02 06:00:00   1st Qu.:10899  
##  Median :2006-07-02 12:00:00   Median :16535  
##  Mean   :2006-07-02 12:00:00   Mean   :18330  
##  3rd Qu.:2013-10-01 18:00:00   3rd Qu.:25545  
##  Max.   :2021-01-01 00:00:00   Max.   :31262

Splitting Data

train <- data[1:23, 2]
test <- data[24:30, 2]
data.ts <- ts(data)
train.ts <- ts(train)
test.ts <- ts(test)

MOVING AVERAGE

# Parameter Optimum
#Fungsi Mencari Parameter Optimum
pemulusan <- function(x,min,max,method=c("SMA","DMA"),
                      metrik=c("SSE","MSE","RMSE","MAD","MAPE")){
  df.master <- data.frame()
  if(method=="SMA"){
    for(i in min:max) {
      sma <- TTR::SMA(x,i)
      ramal <- c(NA,sma)
      df <- cbind(aktual=c(x,NA),mulus=c(sma,NA),ramal)
      error <- df[,1]-df[,3]
      sse <- sum(error^2,na.rm=T)
      mse <- mean(error^2,na.rm=T)
      rmse <- sqrt(mse)
      mad <- mean(abs(error),na.rm=T)
      rerror <- error/df[,1]*100
      mape <- mean(abs(rerror),na.rm=T)
      ak <- data.frame(n=i,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
      df.master <- rbind(df.master,ak)
    }
  }
else if(method=="DMA") {
    for(i in min:max) {
      df_ts_sma <- TTR::SMA(x,  n=i)
      df_ts_dma <- TTR::SMA(df_ts_sma, n=i)
      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(x,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
      error.dma <- df_dma[, 1] - df_dma[, 3]
      SSE.dma <- sum(error.dma^2, na.rm = T)
      MSE.dma <- mean(error.dma^2, na.rm = T)
      RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))
      MAD.dma <- mean(abs(error.dma), na.rm = T)
      r.error.dma <- (error.dma/df_dma[, 1])*100 
      MAPE.dma <- mean(abs(r.error.dma), na.rm = T)
      ak <- data.frame(n=i,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
                       MAD=MAD.dma,MAPE=MAPE.dma)
      df.master <- rbind(df.master,ak)
    }
}
  opt <- df.master$n[which.min(df.master[,metrik])]
  return(list(Akurasi=df.master,n.optimum=paste("n optimum adalah",opt)))
}

pemulusan(train.ts, min = 2, max = 10, method = "DMA", metrik = "RMSE")

## $Akurasi
##    n       SSE      MSE     RMSE      MAD     MAPE
## 1  2 268926428 13446321 3666.923 2334.075 12.34473
## 2  3 216076690 12004261 3464.717 2334.185 12.88494
## 3  4 216056591 13503537 3674.716 2779.602 15.76814
## 4  5 229315379 16379670 4047.180 3273.729 17.79996
## 5  6 230880655 19240055 4386.349 3767.181 19.14636
## 6  7 241097032 24109703 4910.163 4385.727 20.65743
## 7  8 203074367 25384296 5038.283 4254.410 17.74540
## 8  9 172477644 28746274 5361.555 4375.720 15.92533
## 9 10 162851960 40712990 6380.673 5796.140 19.77973
## 
## $n.optimum
## [1] "n optimum adalah 3"

# SMA 
df_ts_sma <- SMA(train.ts,  n=3)

# DMA
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(train.ts,NA), pemulusan_dma = c(pemulusan_dma,NA), ramal_dma)
head(df_dma)

##      df_aktual pemulusan_dma ramal_dma
## [1,]      9819            NA        NA
## [2,]     10038            NA        NA
## [3,]     11004            NA        NA
## [4,]     10990            NA        NA
## [5,]     10869      11583.89        NA
## [6,]     12308      12153.22  11583.89

# Akurasi DMA
error.dma <- df_dma[, 1] - df_dma[, 3]
sse.dma <- sum(error.dma^2, na.rm = T)
mse.dma <- mean(error.dma^2, na.rm = T)
rmse.dma <- sqrt(mean(error.dma^2, na.rm = T))
mad.dma <- mean(abs(error.dma), na.rm = T)
r.error.dma <- (error.dma/df_dma[, 1])*100 
mape.dma <- mean(abs(r.error.dma), na.rm = T)

#Plot
ts.plot(train.ts, xlab = "periode waktu", ylab="Yt", col="blue", lty=3)
points(train.ts)
lines(pemulusan_dma , col = "red",lwd = 1.5)
lines(ramal_dma,col = "black",lwd = 2)
title("Rataan Bergerak  Berganda N=3", cex.main = 1, font.main=4 ,col.main = "black")
legend("topleft", c("Data Aktual","Pemulusan","Ramalan"),lty=1:3,col=c ("blue","red","black"))

EXPONENTIAL SMOOTHING

# Mencari Parameter Optimum
pemulusan2 <- function(x,alfa=NULL,beta=NULL,method=c("SES","DES"),
                       metrik=c("SSE","MSE","RMSE")){
  df.master <- data.frame()
  if(method=="SES"){
    for(i in alfa){
      df_ses <- HoltWinters(x, alpha = i, beta=F, gamma=F)
      sse <- df_ses$SSE 
      mse <- sse/length(x)
      rmse <-sqrt(mse)
      ak <- data.frame(Alpha=i,SSE=sse,MSE=mse,RMSE=rmse)
      df.master <- rbind(df.master,ak)
      datases <- data.frame(x, c(NA, df_ses$fitted[,1]))
      colnames(datases) = c("y","yhat")
    }
  }
  else if(method=="DES"){
    for (i in alfa) {
      for (j in beta) {
        df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse)
        df.master <- rbind(df.master,ak)
        datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
        colnames(datades) = c("y","yhat")
      }
    }
  }
  if(method=="SES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt)))
  }
  else if(method=="DES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    opt2 <- df.master$Beta[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt,
                                                  ",Beta optimum adalah",opt2)))
  }
}

# Double Exponential Smoothing
pemulusan2(train.ts,alfa=c(seq(0.1, 0.9, by = 0.1)),
           beta=c(seq(0.1, 0.9, by = 0.1)),method = "DES",metrik = "RMSE")

## $Akurasi
##    Alpha Beta       SSE      MSE     RMSE
## 1    0.1  0.1 604790486 26295239 5127.888
## 2    0.1  0.2 555592324 24156188 4914.895
## 3    0.1  0.3 510847958 22210781 4712.832
## 4    0.1  0.4 470727526 20466414 4523.982
## 5    0.1  0.5 435672634 18942288 4352.274
## 6    0.1  0.6 406091405 17656148 4201.922
## 7    0.1  0.7 382174495 16616282 4076.307
## 8    0.1  0.8 363810469 15817846 3977.166
## 9    0.1  0.9 350578049 15242524 3904.168
## 10   0.2  0.1 374096937 16265084 4032.999
## 11   0.2  0.2 326832544 14210111 3769.630
## 12   0.2  0.3 293732227 12770966 3573.649
## 13   0.2  0.4 273026481 11870717 3445.391
## 14   0.2  0.5 261627621 11375114 3372.701
## 15   0.2  0.6 255883661 11125377 3335.472
## 16   0.2  0.7 252515693 10978943 3313.449
## 17   0.2  0.8 249243031 10836654 3291.907
## 18   0.2  0.9 244963356 10650581 3263.523
## 19   0.3  0.1 266012919 11565779 3400.850
## 20   0.3  0.2 236723491 10292326 3208.165
## 21   0.3  0.3 222236976  9662477 3108.453
## 22   0.3  0.4 216004158  9391485 3064.553
## 23   0.3  0.5 212881387  9255712 3042.320
## 24   0.3  0.6 210252310  9141405 3023.476
## 25   0.3  0.7 207666707  9028987 3004.827
## 26   0.3  0.8 205801676  8947899 2991.304
## 27   0.3  0.9 205531314  8936144 2989.338
## 28   0.4  0.1 215979826  9390427 3064.380
## 29   0.4  0.2 200824663  8731507 2954.912
## 30   0.4  0.3 196063384  8524495 2919.674
## 31   0.4  0.4 195230387  8488278 2913.465
## 32   0.4  0.5 195689745  8508250 2916.890
## 33   0.4  0.6 197288395  8577756 2928.781
## 34   0.4  0.7 200618186  8722530 2953.393
## 35   0.4  0.8 206041293  8958317 2993.045
## 36   0.4  0.9 213461302  9280926 3046.461
## 37   0.5  0.1 192622304  8374883 2893.939
## 38   0.5  0.2 186489552  8108241 2847.497
## 39   0.5  0.3 187230758  8140468 2853.150
## 40   0.5  0.4 190393961  8277998 2877.151
## 41   0.5  0.5 195035675  8479812 2912.012
## 42   0.5  0.6 201330272  8753490 2958.630
## 43   0.5  0.7 209302962  9100129 3016.642
## 44   0.5  0.8 218593182  9504051 3082.864
## 45   0.5  0.9 228622911  9940127 3152.797
## 46   0.6  0.1 182297260  7925968 2815.310
## 47   0.6  0.2 182107011  7917696 2813.840
## 48   0.6  0.3 186744084  8119308 2849.440
## 49   0.6  0.4 193301235  8404402 2899.035
## 50   0.6  0.5 201318044  8752958 2958.540
## 51   0.6  0.6 210694682  9160638 3026.655
## 52   0.6  0.7 221113310  9613622 3100.584
## 53   0.6  0.8 232145025 10093262 3176.989
## 54   0.6  0.9 243456988 10585086 3253.473
## 55   0.7  0.1 179002583  7782721 2789.753
## 56   0.7  0.2 183057048  7959002 2821.170
## 57   0.7  0.3 190794762  8295424 2880.178
## 58   0.7  0.4 200215987  8705043 2950.431
## 59   0.7  0.5 210981363  9173103 3028.713
## 60   0.7  0.6 222918215  9692096 3113.213
## 61   0.7  0.7 235821298 10253100 3202.046
## 62   0.7  0.8 249592484 10851847 3294.214
## 63   0.7  0.9 264327015 11492479 3390.056
## 64   0.8  0.1 179990863  7825690 2797.443
## 65   0.8  0.2 187415088  8148482 2854.555
## 66   0.8  0.3 197953702  8606683 2933.715
## 67   0.8  0.4 210184554  9138459 3022.988
## 68   0.8  0.5 223904988  9734999 3120.096
## 69   0.8  0.6 239085761 10395033 3224.133
## 70   0.8  0.7 255781653 11120941 3334.808
## 71   0.8  0.8 274177870 11920777 3452.648
## 72   0.8  0.9 294562007 12807044 3578.693
## 73   0.9  0.1 183916399  7996365 2827.785
## 74   0.9  0.2 194339136  8449528 2906.807
## 75   0.9  0.3 207721951  9031389 3005.227
## 76   0.9  0.4 223055826  9698079 3114.174
## 77   0.9  0.5 240288453 10447324 3232.232
## 78   0.9  0.6 259530600 11283939 3359.158
## 79   0.9  0.7 280951176 12215269 3495.035
## 80   0.9  0.8 304734659 13249333 3639.963
## 81   0.9  0.9 331002001 14391391 3793.599
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.7 ,Beta optimum adalah 0.1"

df_des <- HoltWinters(train.ts, alpha = 0.7, beta = 0.1, gamma = F)
df_des

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = train.ts, alpha = 0.7, beta = 0.1, gamma = F)
## 
## Smoothing parameters:
##  alpha: 0.7
##  beta : 0.1
##  gamma: FALSE
## 
## Coefficients:
##        [,1]
## a 28483.315
## b  1084.889

datades <- data.frame(train.ts, c(NA,NA, df_des$fitted[,1]))
colnames(datades) <- c("y","yhat")
head(datades)

##       y     yhat
## 1  9819       NA
## 2 10038       NA
## 3 11004 10257.00
## 4 10990 11051.19
## 5 10869 11275.36
## 6 12308 11229.47

fore.des <- forecast(df_des, h = 5)

AKURASI DESSSS!!!

# Menghitung SSE (Sum of Squared Errors)
sse.des <- sum((datades$y - datades$yhat)^2)

# Menghitung MSE (Mean Squared Error)
mse.des <- mean((datades$y - datades$yhat)^2)

# Menghitung RMSE (Root Mean Squared Error)
rmse.des <- sqrt(mean((datades$y - datades$yhat)^2))

# Plot
plot(train.ts, xlab="periode  waktu", ylab="Yt",  col="blue", lty=3)
points(train.ts)
lines(datades[,2], col = "red",lwd = 1.5, lty = 1)
title("Double Exponential Smoothing (Alpha = 0.7 dan Beta = 0.1)", cex.main=1, font.main=4 
      ,col.main="black")
legend("topleft", c("Data aktual","Fitted DES"), lty=1:3,col=c ("blue","red"))

tabel1 <- matrix(c(sse.dma, mse.dma, rmse.dma, sse.des, mse.des, rmse.des), ncol = 3, byrow = T)
colnames(tabel1) <- c("SSE","MSE","RMSE")
rownames(tabel1) <- c("DMA","DES")

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

kable(tabel1)

	SSE	MSE	RMSE
DMA	216076690	12004261	3464.717
DES	NA	NA	NA

ACF & PACF

acf(data.ts[,2], lag.max = 35)

adf.test(data.ts[,2])

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.ts[, 2]
## Dickey-Fuller = -1.7871, Lag order = 3, p-value = 0.6545
## alternative hypothesis: stationary

Uji Statistik menghasilkan p-value yang lebih besar dari 0.05 yang berarti pada taraf nyata 5% dapat dikatakan data tidak stasioner sehingga diperlukan penanganan ketidakstasioneran data sebelum meentukan model.

Differencing

data.diff <- diff(data.ts[,2], differences = 2)

# cek  Stasioner Eksploratif
plot(data.diff, type = "l")

# Uji Statistik
adf.test(data.diff)

## Warning in adf.test(data.diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff
## Dickey-Fuller = -4.86, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

acf(data.diff, lag.max = 20)

pacf(data.diff, lag.max = 20)

ARIMA 1,2,1 ARIMA 1,2,2 ARIMA 1,2,4

ARIMA 1,2,1

model1 <- Arima(data.diff, order = c(1,2,1))

summary(model1)

## Series: data.diff 
## ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.6282  -1.0000
## s.e.   0.1425   0.1008
## 
## sigma^2 = 45403517:  log likelihood = -267.43
## AIC=540.86   AICc=541.95   BIC=544.63
## 
## Training set error measures:
##                    ME     RMSE    MAE      MPE     MAPE     MASE       ACF1
## Training set 253.6855 6238.373 4384.8 101.8552 187.0319 0.733435 -0.4254442

ARIMA 1,2,2

model2 <- Arima(data.diff, order = c(1,2,2))

summary(model2)

## Series: data.diff 
## ARIMA(1,2,2) 
## 
## Coefficients:
##           ar1      ma1     ma2
##       -0.4806  -1.9722  1.0000
## s.e.   0.1672   0.1539  0.1526
## 
## sigma^2 = 20582572:  log likelihood = -259.51
## AIC=527.03   AICc=528.93   BIC=532.06
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set -47.77289 4111.826 3034.041 30.69171 207.9522 0.5074967 -0.264613

ARIMA(1,2,4)

model3 <- Arima(data.diff, order = c(1,2,4))

summary(model3)

## Series: data.diff 
## ARIMA(1,2,4) 
## 
## Coefficients:
##           ar1      ma1     ma2     ma3      ma4
##       -0.8164  -1.7957  0.0001  1.7955  -0.9999
## s.e.   0.1246   0.2238  0.3390  0.3385   0.2245
## 
## sigma^2 = 13805548:  log likelihood = -255.01
## AIC=522.01   AICc=526.43   BIC=529.56
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
## Training set -10.92311 3217.788 2182.375 -71.16775 202.228 0.3650406 -0.2483874

MODEL DENGAN PARAMETER OPTIMUM

model4 <- auto.arima(data.diff)

summary(model4)

## Series: data.diff 
## ARIMA(2,0,0) with zero mean 
## 
## Coefficients:
##           ar1      ar2
##       -0.7864  -0.4590
## s.e.   0.1669   0.1621
## 
## sigma^2 = 13466115:  log likelihood = -268.92
## AIC=543.84   AICc=544.84   BIC=547.84
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -92.58338 3536.135 2482.085 53.60804 183.6667 0.4151723
##                     ACF1
## Training set -0.09479376

aic1 <- model1$aic
aic2 <- model2$aic
aic3 <- model3$aic
aic4 <- model4$aic

bic1 <- model1$bic
bic2 <- model2$bic
bic3 <- model3$bic
bic4 <- model4$bic

tabel <- matrix(c(aic1,aic2,aic3,aic4,bic1,bic2,bic3,bic4), ncol = 4, byrow = T)
colnames(tabel) <- c("Model 1","Model 2","Model 3","Model 4")
rownames(tabel) <- c("AIC","BIC")

kable(tabel)

	Model 1	Model 2	Model 3	Model 4
AIC	540.8560	527.0280	522.0126	543.8412
BIC	544.6303	532.0604	529.5611	547.8378

## Forecasting
ramalan <- forecast(model3, h = 5)
ramalan

##    Point Forecast     Lo 80    Hi 80      Lo 95     Hi 95
## 31     -3177.4413 -8223.201 1868.319 -10894.265  4539.383
## 32       665.6866 -5098.505 6429.878  -8149.883  9481.256
## 33       214.8347 -5550.493 5980.163  -8602.474  9032.143
## 34       667.9969 -5510.594 6846.588  -8781.343 10117.337
## 35       383.1550 -5802.209 6568.519  -9076.543  9842.853

data.ramalan <- ramalan$mean

plot(ramalan)