PERCOBAAN

PROPOSAL

Kontak	: \(\downarrow\)
Email	clara.evania@student.matanauniversity.ac.id
Instagram	https://www.instagram.com/claraevania/
RPubs	https://rpubs.com/claradellaevania/

library(forecast)
library(ggplot2)
library(tidyverse)
library(tidyr)
library(plotly)
library(tseries)
library(lmtest)
library(xgboost)
#library(tidymodels)
library(modeltime)
library(tidyverse)
library(timetk)
library(dplyr)

# This toggles plots from plotly (interactive) to ggplot (static)
interactive <- FALSE

Prediksi pertama menggunakan data Januari 2009-Januari 2020

data <- read.csv("C:/Users/brigi/Downloads/dataawal.csv", sep = ";")
data

timeseries = ts(data$Jumlah, start = c(2009,1), frequency = 12)
plot(timeseries)

timeseries = ts(data$Jumlah, start = c(2009,1), frequency = 12)
plot(timeseries)

decompose(timeseries) %>%
  autoplot()

hw.add = HoltWinters(timeseries, alpha = NULL, beta = NULL, gamma = NULL, seasonal = "additive")
hw.add

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = timeseries, alpha = NULL, beta = NULL, gamma = NULL,     seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 0.8306758
##  beta : 0
##  gamma: 0.5967294
## 
## Coefficients:
##              [,1]
## a   546866.068550
## b     1386.296620
## s1   -8514.529094
## s2  -16677.609654
## s3  -12691.774384
## s4  -14154.661548
## s5   15072.346969
## s6   54801.147318
## s7   24519.387280
## s8   19459.136027
## s9      -5.922192
## s10 -40886.672799
## s11  11672.850034
## s12 -10634.149381

decompose(timeseries) %>%
  autoplot()

prediksi = forecast(hw.add, h= 12)
prediksi

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2020       539737.8 498985.7 580490.0 477412.8 602062.8
## Mar 2020       532961.1 479982.9 585939.2 451938.0 613984.1
## Apr 2020       538333.2 475463.3 601203.1 442181.9 634484.4
## May 2020       538256.6 466852.3 609660.9 429053.2 647460.0
## Jun 2020       568869.9 489847.6 647892.2 448015.8 689724.0
## Jul 2020       609985.0 524017.2 695952.8 478508.6 741461.4
## Aug 2020       581089.5 488696.8 673482.2 439787.1 722392.0
## Sep 2020       577415.6 479016.6 675814.5 426927.4 727903.8
## Oct 2020       559336.8 455277.7 663395.9 400192.1 718481.5
## Nov 2020       519842.4 410415.5 629269.3 352488.4 687196.3
## Dec 2020       573788.2 459244.8 688331.6 398609.2 748967.2
## Jan 2021       552867.5 433426.6 672308.4 370198.3 735536.6

accuracy(hw.add$fitted[,1], timeseries)

##                ME     RMSE     MAE       MPE     MAPE       ACF1 Theil's U
## Test set 1796.146 31718.28 22078.1 0.1032714 6.123455 0.02872675 0.8097799

train = timeseries %>% head(-12)
test = timeseries %>% tail(12)

autoplot(prediksi)+autolayer(test)

mse.add = hw.add$SSE/frequency(hw.add$fitted)
rmse.add = sqrt(mse.add)
mape.add = mean(abs(timeseries-hw.add$fitted[,1])/timeseries)*100
mse.add

## [1] 10144332749

rmse.add

## [1] 100719.1

mape.add

## [1] 6.123455

Prediksi kedua menggunakan data Januari 2009-Januari 2021

data1 <- read.csv("C:/Users/brigi/Downloads/dataaforecast1.csv", sep = ";")
data1

timeseries1 = ts(data1$Jumlah, start = c(2009,1), frequency = 12)
plot(timeseries1)

hw.add1 = HoltWinters(timeseries1, alpha = NULL, beta = NULL, gamma = NULL, seasonal = "additive")
hw.add1

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = timeseries1, alpha = NULL, beta = NULL, gamma = NULL,     seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 0.8306761
##  beta : 0
##  gamma: 0.5967353
## 
## Coefficients:
##              [,1]
## a   563501.962465
## b     1386.296620
## s1   -8514.413214
## s2  -16677.578821
## s3  -12691.700948
## s4  -14154.583937
## s5   15072.330305
## s6   54801.238058
## s7   24519.441066
## s8   19459.129228
## s9      -6.017413
## s10 -40886.893524
## s11  11672.813265
## s12 -10634.001623

prediksi1 = forecast(hw.add1, h= 12)
prediksi1

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2021       556373.8 517512.6 595235.1 496940.7 615807.0
## Mar 2021       549597.0 499077.0 600117.0 472333.3 626860.6
## Apr 2021       554969.2 495016.3 614922.0 463279.2 646659.1
## May 2021       554892.6 486801.4 622983.8 450756.0 659029.1
## Jun 2021       585505.8 510150.1 660861.5 470259.1 700752.4
## Jul 2021       626621.0 544642.0 708600.0 501244.9 751997.0
## Aug 2021       597725.5 509619.7 685831.3 462979.3 732471.6
## Sep 2021       594051.5 500218.1 687884.8 450545.7 737557.2
## Oct 2021       575972.6 476741.7 675203.5 424212.1 727733.2
## Nov 2021       536478.0 432128.4 640827.7 376889.1 696067.0
## Dec 2021       590424.0 481195.3 699652.8 423373.1 757475.0
## Jan 2022       569503.5 455604.5 683402.5 395310.0 743697.0

accuracy(hw.add1$fitted[,1], timeseries1)

##               ME     RMSE      MAE        MPE     MAPE       ACF1 Theil's U
## Test set 1634.09 30253.56 20086.12 0.09395409 5.570968 0.02891877 0.8013517

train1 = timeseries1 %>% head(-12)
test1 = timeseries1 %>% tail(12)

autoplot(prediksi1)+autolayer(test1)

Prediksi kedua menggunakan data Januari 2009-Januari 2022

data2 <- read.csv("C:/Users/brigi/Downloads/dataaforecast2.csv", sep = ";")
data2

timeseries2 = ts(data2$Jumlah, start = c(2009,1), frequency = 12)
plot(timeseries2)

hw.add2 = HoltWinters(timeseries2, alpha = NULL, beta = NULL, gamma = NULL, seasonal = "additive")
hw.add2

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = timeseries2, alpha = NULL, beta = NULL, gamma = NULL,     seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 0.8306753
##  beta : 0
##  gamma: 0.5967365
## 
## Coefficients:
##              [,1]
## a   580137.885813
## b     1386.296620
## s1   -8514.314173
## s2  -16677.533245
## s3  -12691.642171
## s4  -14154.494075
## s5   15072.330018
## s6   54801.292766
## s7   24519.513228
## s8   19459.137957
## s9      -6.082491
## s10 -40887.066621
## s11  11672.741778
## s12 -10633.912212

prediksi2 = forecast(hw.add2, h= 12)
prediksi2

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2022       573009.9 535798.5 610221.2 516100.0 629919.7
## Mar 2022       566232.9 517857.9 614608.0 492249.7 640216.2
## Apr 2022       571605.1 514197.8 629012.5 483808.1 659402.1
## May 2022       571528.6 506328.4 636728.8 471813.4 671243.7
## Jun 2022       602141.7 529985.4 674298.0 491788.2 712495.2
## Jul 2022       643257.0 564758.6 721755.3 523204.1 763309.8
## Aug 2022       614361.5 529996.5 698726.5 485336.3 743386.6
## Sep 2022       610687.4 520838.0 700536.8 473274.6 748100.2
## Oct 2022       592608.5 497590.7 687626.3 447291.3 737925.6
## Nov 2022       553113.8 453194.6 653033.0 400300.6 705927.0
## Dec 2022       607059.9 502468.8 711651.0 447101.6 767018.2
## Jan 2023       586139.5 477076.4 695202.7 419341.9 752937.2

accuracy(hw.add2$fitted[,1], timeseries2)

##               ME     RMSE      MAE        MPE     MAPE       ACF1 Theil's U
## Test set 1498.86 28974.66 18423.84 0.08617929 5.109925 0.02914809 0.7936465

train2 = timeseries2 %>% head(-12)
test2 = timeseries2 %>% tail(12)

autoplot(prediksi2)+autolayer(test2)

datagabung <- read.csv("C:/Users/brigi/Downloads/datagabung.csv", sep = ";")
datagabung

timeseries3 = ts(datagabung$Jumlah, start = c(2009,1), frequency = 12)
plot(timeseries3)

timeseries_df <- data.frame(
  Date = as.Date(time(timeseries3)),
  Jumlah = as.numeric(timeseries3))

library(ggplot2)
p <- ggplot(timeseries_df, aes(x = Date, y = Jumlah)) +
  geom_line() +
  labs(title = "Time Series Plot", x = "Date", y = "Jumlah") +
  theme_minimal()
ggplotly(p)

HOLTS WINTERS

hw_gabung = HoltWinters(timeseries3, alpha = NULL, beta = NULL, gamma = NULL, seasonal = "additive")
hw_gabung

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = timeseries3, alpha = NULL, beta = NULL, gamma = NULL,     seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 0.8450977
##  beta : 0
##  gamma: 0.5693698
## 
## Coefficients:
##           [,1]
## a   436872.706
## b     1386.297
## s1  -33467.893
## s2  -16633.352
## s3  -10840.553
## s4  -11878.237
## s5   15981.202
## s6   55564.657
## s7   25150.936
## s8   18658.547
## s9   -1784.541
## s10 -40696.785
## s11  14841.902
## s12 -14478.615

forecast_hw = forecast(hw_gabung, h= 12)
forecast_hw

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2024       404791.1 360658.5 448923.7 337296.1 472286.1
## Mar 2024       423011.9 365230.4 480793.5 334642.7 511381.2
## Apr 2024       430191.0 361418.1 498964.0 325011.8 535370.2
## May 2024       430539.7 352304.5 508774.8 310889.3 550190.0
## Jun 2024       459785.4 373115.0 546455.8 327234.4 592336.3
## Jul 2024       500755.1 406400.6 595109.7 356452.3 645058.0
## Aug 2024       471727.7 370269.4 573186.1 316560.6 626894.9
## Sep 2024       466621.6 358525.3 574717.9 301302.6 631940.7
## Oct 2024       447564.8 333215.2 561914.4 272682.2 622447.5
## Nov 2024       410038.9 289760.7 530317.1 226089.2 593988.6
## Dec 2024       466963.9 341035.8 592891.9 274373.6 659554.2
## Jan 2025       439029.7 307694.6 570364.7 238170.1 639889.2

accuracy(hw_gabung$fitted[,1], timeseries3)

##                ME     RMSE      MAE        MPE     MAPE       ACF1 Theil's U
## Test set 27.99458 34334.83 19149.18 -0.3471361 5.263091 0.00258918 0.7994985

train3 = timeseries3 %>% head(-12)
test3 = timeseries3 %>% tail(12)

Plot_hw= autoplot(forecast_hw)+autolayer(test3)
Plot_hw

SARIMA

library(ggplot2)
p <- ggplot(timeseries_df, aes(x = Date, y = Jumlah)) +
  geom_line() +
  labs(title = "Time Series Plot", x = "Date", y = "Jumlah") +
  theme_minimal()
ggplotly(p)

diff_data <- diff(timeseries3)
plot(diff_data)

adf.test(diff_data)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_data
## Dickey-Fuller = -7.7484, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

transdiff_data <- diff(log(timeseries3))
plot(transdiff_data)

adf.test(transdiff_data)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  transdiff_data
## Dickey-Fuller = -8.1381, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(transdiff_data, lag.max=37)
Pacf(transdiff_data, lag.max=37)

sarima_model1 <- Arima(log(timeseries3), order = c(1,1,1), seasonal = list(order = c(1,0,0), period = 12), include.mean = FALSE)
summary(sarima_model1)

## Series: log(timeseries3) 
## ARIMA(1,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1
##       0.5446  -0.8712  0.5468
## s.e.  0.0915   0.0505  0.0733
## 
## sigma^2 = 0.008168:  log likelihood = 176.53
## AIC=-345.06   AICc=-344.83   BIC=-332.29
## 
## Training set error measures:
##                      ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.00646496 0.08937008 0.05934942 0.04886972 0.4664839 0.5063927
##                     ACF1
## Training set -0.02439862

modelsarim <- auto.arima(log(timeseries3), seasonal=TRUE)
summary(modelsarim)

## Series: log(timeseries3) 
## ARIMA(2,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1
##       0.5398  0.0238  -0.8773  0.5471
## s.e.  0.0919  0.0830   0.0530  0.0734
## 
## sigma^2 = 0.008211:  log likelihood = 176.57
## AIC=-343.14   AICc=-342.8   BIC=-327.18
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.00655587 0.0893552 0.0591729 0.04957923 0.4650883 0.5048866
##                     ACF1
## Training set -0.01364196

model <- auto.arima((timeseries3), seasonal=TRUE)
summary(model)

## Series: (timeseries3) 
## ARIMA(1,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1
##       0.6521  -0.8939  0.5735
## s.e.  0.0883   0.0501  0.0712
## 
## sigma^2 = 1.192e+09:  log likelihood = -2137.25
## AIC=4282.5   AICc=4282.72   BIC=4295.27
## 
## Training set error measures:
##                   ME     RMSE     MAE         MPE     MAPE      MASE       ACF1
## Training set 1170.36 34143.98 20883.2 -0.02171558 5.892033 0.4584264 0.01946932

checkresiduals(model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,0,0)[12]
## Q* = 16.532, df = 21, p-value = 0.7391
## 
## Model df: 3.   Total lags used: 24

residuals_sarima <- residuals(model)
Box.test(residuals_sarima, lag=24, fitdf=3, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_sarima
## X-squared = 16.532, df = 21, p-value = 0.7391

# Forecasting
forecast_data <- forecast(model, h=12)  # h adalah horizon waktu untuk prediksi
forecast_data

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2024       287846.0 243597.1 332095.0 220173.1 355519.0
## Mar 2024       326912.2 271382.6 382441.9 241987.0 411837.5
## Apr 2024       358147.7 296587.5 419707.9 263999.6 452295.9
## May 2024       379290.1 313909.0 444671.2 279298.4 479281.8
## Jun 2024       404825.2 336723.5 472927.0 300672.6 508977.8
## Jul 2024       443212.6 372991.4 513433.8 335818.5 550606.7
## Aug 2024       433611.1 361620.9 505601.3 323511.6 543710.6
## Sep 2024       426625.7 353081.6 500169.9 314149.6 539101.9
## Oct 2024       400327.6 325367.2 475288.0 285685.6 514969.6
## Nov 2024       367295.3 291010.7 443580.0 250628.0 483962.6
## Dec 2024       412562.7 335017.6 490107.7 293967.8 531157.6
## Jan 2025       377399.7 298640.6 456158.9 256948.0 497851.5

plot(forecast_data)

plotly_forecast <- plot_ly() %>%
  add_lines(x = index(forecast_data$mean), y = forecast_data$mean, name = "Forecast Mean") %>%
  add_ribbons(x = index(forecast_data$mean), ymin = forecast_data$lower[, "95%"], ymax = forecast_data$upper[, "95%"], 
              name = "95% Prediction Interval", line = list(color = "rgba(0,0,0,0)"), fillcolor = "rgba(0,100,80,0.2)")

plotly_forecast <- plotly_forecast %>%
  layout(title = "Forecast Plot using Plotly",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Forecast Value"))
plotly_forecast

accuracy(forecast_data)

##                   ME     RMSE     MAE         MPE     MAPE      MASE       ACF1
## Training set 1170.36 34143.98 20883.2 -0.02171558 5.892033 0.4584264 0.01946932

PERBANDINGAN HW DAN SARIMA

par(mfrow=c(1,2))
autoplot(forecast_hw)+autolayer(test3)

plot(forecast(model, h=12))

subplot(Plot_hw, plot(forecast_data), nrows = 1)

timeseries3

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2009 174541 147704 168205 188776 190803 200566 235198 232255 218443 221282
## 2010 179273 191926 192579 184907 203388 228045 254907 243154 240947 229904
## 2011 209093 207195 207907 224704 209058 245652 283524 258377 258440 247565
## 2012 253286 225993 230957 225488 220700 244080 271512 254079 257363 255021
## 2013 232935 241868 252210 242369 247972 275667 297878 309219 305629 266562
## 2014 279257 275795 276573 280096 286033 330396 361066 336763 354762 341651
## 2015 301748 338991 305272 313763 295973 359702 382683 303621 389060 369447
## 2016 350592 375744 364113 380767 394557 405835 484231 438135 445716 432215
## 2017 460824 453985 425499 477464 489376 504141 592046 601884 550520 465085
## 2018 358065 452423 492678 516777 528512 544550 624366 573766 555903 517889
## 2019 455570 437456 449569 477069 486602 549516 604323 606412 590398 567967
## 2020 536611 539738 532961 538333 538257 568870 609985 581090 577416 559337
## 2021 552868 556374 549597 554969 554893 585506 626621 597726 594052 575973
## 2022 569504 573010 566233 571605 571529 602142 643257 614362 610687 592609
## 2023 586140 323510 370695 411510 439475 478198 541353 522141 508350 461441
## 2024 420037                                                               
##         Nov    Dec
## 2009 184803 222546
## 2010 199861 227251
## 2011 221603 253591
## 2012 242781 268072
## 2013 307276 299013
## 2014 296876 347370
## 2015 270935 370640
## 2016 413232 442800
## 2017 361006 315909
## 2018 406725 498819
## 2019 497925 552403
## 2020 519842 573788
## 2021 536478 590424
## 2022 553114 607060
## 2023 403154 481646
## 2024

dftimeseries <-  as.data.frame(timeseries3)
dftimeseries <- tibble::rownames_to_column(dftimeseries, "date")
dftimeseries

# Model 2: arima_boost ----
model_fit_arima_boosted <- arima_boost(
    min_n = 2,
    learn_rate = 0.015
) %>%
    set_engine(engine = "auto_arima_xgboost") %>%
    fit(value ~ date + as.numeric(date) + factor(month(date, label = TRUE), ordered = F),
        data = training(timeseries))

datagabung4 <- read.csv("C:/Users/brigi/Downloads/datagabung4.csv", sep = ";")
datagabung4

data <- datagabung4 %>%
  mutate(Tahun = parse_date_time(Tahun, orders = "my", locale = "id_ID")) %>%
  mutate(Tahun = format(Tahun, "%B %Y", locale = "id_ID"))
data

# Split Data 80/20
splits <- initial_time_split(data, prop = 0.9)

# Model 1: auto_arima ----
model_fit_arima_no_boost <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(value ~ date, data = training(splits))
#> frequency = 12 observations per 1 year

class(training(splits))
str(training(splits))
training_data <- as.data.frame(training(splits))
names(training_data)
splits <- initial_time_split(data, prop = 0.9)
training_data <- training(splits)
splits

model_fit_arima_no_boost <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(value ~ date, data = training(splits))
#> frequency = 12 observations per 1 year

# Fit the ARIMA model using `auto_arima`
model_fit_arima_no_boost <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(value ~ date, data = training(splits))

# Print model summary
print(model_fit_arima_no_boost)