Melakukan Instalasi Packages yang diperlukan

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(TTR)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(graphics)
library(forecast)
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ readr::spec()   masks TSA::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rmarkdown)
library(readxl)
dataawal <- read_excel("data penjualan produk farmasi.xlsx")
data1<-dataawal[-111:-120,]
data2<-dataawal[-1:-110,]
data.ts<-ts(data1)
head(data.ts)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
##      Data Penjualan Produk Farmasi
## [1,]                       10618.1
## [2,]                       10537.9
## [3,]                       10209.3
## [4,]                       10553.0
## [5,]                        9934.9
## [6,]                       10534.5

data1 merupakan data yang akan digunakan untuk mengidentifikasi model dan meramal 10 minggu kedepan, sedangkan data2 akan digunakan untuk mengetes tingkat akurasi model.

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

Setelah proses persiapan sudah selesai, langkah selanjutnya adalah mengidentifikasi model yang sesuai.

adf.test(data.ts)
## Warning in adf.test(data.ts): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.ts
## Dickey-Fuller = -5.9291, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
acf(data.ts)

pacf(data.ts, lag.max = 50)

eacf(data.ts)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o x o o o o o o o o  o  o  o 
## 1 o o o o o o o o o o o  o  o  o 
## 2 o x 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 o o o o o o o o o o  o  o  o 
## 5 o o x o o o o o o o o  o  o  o 
## 6 x o o o o o o o o o o  o  o  o 
## 7 x x o o x o o o o o o  o  o  o

Berdasarkan ketiga plot tersebut, berikut kandidat model yang sesuai ACF = Kandidat model 1 : ARIMA(0,0,3) PACF = Kandidat model 2 : ARIMA(3,0,0) EACF = Kandidat model 3 : ARIMA (1,0,1) EACF = Kandidat model 4 : ARIMA (1,0,2) EACF = Kandidat model 4 : ARIMA (3,0,1)

#Penentuan Model Terbaik berdasar AIC

arima(data.ts, order=c(0,0,3), method="ML")
## 
## Call:
## arima(x = data.ts, order = c(0, 0, 3), method = "ML")
## 
## Coefficients:
##          ma1      ma2      ma3   intercept
##       0.0141  -0.0518  -0.2326  10374.4514
## s.e.  0.1002   0.1171   0.0967     14.9337
## 
## sigma^2 estimated as 45079:  log likelihood = -745.56,  aic = 1499.12
arima(data.ts, order=c(3,0,0), method="ML") 
## 
## Call:
## arima(x = data.ts, order = c(3, 0, 0), method = "ML")
## 
## Coefficients:
##          ar1     ar2      ar3   intercept
##       0.0904  0.0086  -0.2312  10374.9540
## s.e.  0.0927  0.0929   0.0926     17.9023
## 
## sigma^2 estimated as 44745:  log likelihood = -745.15,  aic = 1498.3
arima(data.ts, order=c(1,0,1), method="ML")
## 
## Call:
## arima(x = data.ts, order = c(1, 0, 1), method = "ML")
## 
## Coefficients:
##          ar1     ma1   intercept
##       0.0672  0.0246  10375.1957
## s.e.  0.4398  0.4314     22.7734
## 
## sigma^2 estimated as 47350:  log likelihood = -748.18,  aic = 1502.36
arima(data.ts, order=c(1,0,1), method="ML") 
## 
## Call:
## arima(x = data.ts, order = c(1, 0, 1), method = "ML")
## 
## Coefficients:
##          ar1     ma1   intercept
##       0.0672  0.0246  10375.1957
## s.e.  0.4398  0.4314     22.7734
## 
## sigma^2 estimated as 47350:  log likelihood = -748.18,  aic = 1502.36
arima(data.ts, order=c(3,0,1), method="ML") 
## 
## Call:
## arima(x = data.ts, order = c(3, 0, 1), method = "ML")
## 
## Coefficients:
##          ar1      ar2      ar3      ma1   intercept
##       0.6498  -0.0437  -0.2197  -0.6094  10374.0433
## s.e.  0.2868   0.1149   0.1105   0.2950     12.8663
## 
## sigma^2 estimated as 43484:  log likelihood = -743.64,  aic = 1497.28

Berdasarkan nilai AIC, model yang mempunyai AIC terkecil adalah model ARIMA (3,0,1) sehingga model yang akan digunakan untuk peramalan adalah model ARIMA(3,0,1)

ARIMA_data <- Arima(data.ts, order=c(3,0,1), method="ML")

Setelah mendapatkan model terbaik, selanjutnya dilakukan peramalan untuk 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
## 111       10371.02 10097.49 10644.55 9952.692 10789.35
## 112       10361.97 10088.22 10635.72 9943.305 10780.64
## 113       10362.54 10088.75 10636.33 9943.808 10781.27
## 114       10367.76 10086.66 10648.86 9937.849 10797.67
## 115       10373.12 10088.65 10657.58 9938.063 10808.17
## 116       10376.24 10090.72 10661.76 9939.579 10812.91
## 117       10376.89 10091.37 10662.41 9940.230 10813.56
## 118       10376.00 10090.29 10661.72 9939.037 10812.97
## 119       10374.71 10088.73 10660.69 9937.338 10812.08
## 120       10373.76 10087.68 10659.84 9936.242 10811.29
plot(hasil_forecastt)

Garis biru mewakili prediksi titik, area berbayang biru tua menunjukkan interval prediksi 80% dan area berbayang biru muda menunjukkan interval prediksi 95% untuk prediksi titik. 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       10371.02
## 2       10361.97
## 3       10362.54
## 4       10367.76
## 5       10373.12
## 6       10376.24
dataaktual <- as.data.frame(data2)
dataaktual
##    Data Penjualan Produk Farmasi
## 1                        10210.1
## 2                        10352.5
## 3                        10423.8
## 4                        10519.3
## 5                        10596.7
## 6                        10650.0
## 7                        10741.6
## 8                        10246.0
## 9                        10354.4
## 10                       10155.4
akurasi.arima <- accuracy(data_forecast1$data_forecast1,dataaktual$`Sales (in thousands)`)
akurasi.arima
##           ME RMSE MAE MPE MAPE
## Test set NaN  NaN NaN NaN  NaN

Jika nilai MAPE kecil maka kesalahan hasil pendugaannya juga kecil. Peramalan penjualan menggunakan metode ARMA dinilai memiliki kemampuan ramalan yang sangat akurat karena hasil dari MAPE adalah sebesar 1.538676%

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 10393.24
SSE <- data1.ses$SSE
SSE
## [1] 5954375
MSE <- SSE/110
MSE
## [1] 54130.69
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
## 111       10393.24 10093.50 10692.98 9934.830 10851.65
## 112       10393.24 10092.01 10694.48 9932.543 10853.94
## 113       10393.24 10090.52 10695.96 9930.268 10856.22
## 114       10393.24 10089.04 10697.44 9928.004 10858.48
## 115       10393.24 10087.57 10698.92 9925.751 10860.73
## 116       10393.24 10086.10 10700.38 9923.509 10862.97
## 117       10393.24 10084.64 10701.84 9921.278 10865.21
## 118       10393.24 10083.19 10703.29 9919.057 10867.43
## 119       10393.24 10081.74 10704.74 9916.846 10869.64
## 120       10393.24 10080.30 10706.18 9914.645 10871.84
plot(hasil_forecast.sess)

Garis biru mewakili prediksi titik, area berbayang biru tua menunjukkan interval prediksi 80% dan area berbayang biru muda menunjukkan interval prediksi 95% untuk prediksi titik

data_forecast2 <- hasil_forecast.ses$`Point Forecast`
data_forecast2 <- as.data.frame(data_forecast2)
head(data_forecast2)
##   data_forecast2
## 1       10393.24
## 2       10393.24
## 3       10393.24
## 4       10393.24
## 5       10393.24
## 6       10393.24
dataaktual <- as.data.frame(data2)
dataaktual
##    Data Penjualan Produk Farmasi
## 1                        10210.1
## 2                        10352.5
## 3                        10423.8
## 4                        10519.3
## 5                        10596.7
## 6                        10650.0
## 7                        10741.6
## 8                        10246.0
## 9                        10354.4
## 10                       10155.4
akurasi.ses <- accuracy(data_forecast2$data_forecast2,dataaktual$`Sales (in thousands)`)
akurasi.ses
##           ME RMSE MAE MPE MAPE
## Test set NaN  NaN NaN NaN  NaN

Jika nilai MAPE kecil maka kesalahan hasil pendugaannya juga kecil. Peramalan penjualan menggunakan metode ARMA dinilai memiliki kemampuan ramalan yang sangat akurat karena hasil dari MAPE adalah sebesar 1.540696%

Setelah dilakukan peramalan dengan menggunakan dua metode, akan diamati nilai MAPE dari masing-masing metode. Metode dengan nilai MAPE paling kecil akan terpilih menjadi model terbaik.

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

Hasil dari kedua metode tersebut menghasilkan nilai MAPE 1.538676% untuk metode ARIMA dan 1.540696% untuk metode Simple Exponential Smoothing. Dengan demikian, metode yang mempunyai MAPE lebih kecil adalah metode ARIMA(3,0,1) dengan nilai p=3,d=0,q=1

Berikut interval prediksi dari hasil ramalan model ARIMA(3,0,1)

fit <- Arima(data.ts, order=c(3,0,1), method="ML")
fit
## Series: data.ts 
## ARIMA(3,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ar3      ma1        mean
##       0.6498  -0.0437  -0.2197  -0.6094  10374.0433
## s.e.  0.2868   0.1149   0.1105   0.2950     12.8663
## 
## sigma^2 = 45555:  log likelihood = -743.64
## AIC=1499.28   AICc=1500.1   BIC=1515.49
hasil_forecast
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 111       10371.02 10097.49 10644.55 9952.692 10789.35
## 112       10361.97 10088.22 10635.72 9943.305 10780.64
## 113       10362.54 10088.75 10636.33 9943.808 10781.27
## 114       10367.76 10086.66 10648.86 9937.849 10797.67
## 115       10373.12 10088.65 10657.58 9938.063 10808.17
## 116       10376.24 10090.72 10661.76 9939.579 10812.91
## 117       10376.89 10091.37 10662.41 9940.230 10813.56
## 118       10376.00 10090.29 10661.72 9939.037 10812.97
## 119       10374.71 10088.73 10660.69 9937.338 10812.08
## 120       10373.76 10087.68 10659.84 9936.242 10811.29
interval.prediction <- hasil_forecast[,-1]
interval.prediction
##        Lo 80    Hi 80    Lo 95    Hi 95
## 111 10097.49 10644.55 9952.692 10789.35
## 112 10088.22 10635.72 9943.305 10780.64
## 113 10088.75 10636.33 9943.808 10781.27
## 114 10086.66 10648.86 9937.849 10797.67
## 115 10088.65 10657.58 9938.063 10808.17
## 116 10090.72 10661.76 9939.579 10812.91
## 117 10091.37 10662.41 9940.230 10813.56
## 118 10090.29 10661.72 9939.037 10812.97
## 119 10088.73 10660.69 9937.338 10812.08
## 120 10087.68 10659.84 9936.242 10811.29
plot(hasil_forecast$`Point Forecast`, type="n", ylim=range(hasil_forecast$`Lo 80`,hasil_forecast$`Hi 80`))
polygon(c(time(hasil_forecast$`Point Forecast`),rev(time(hasil_forecast$`Point Forecast`))), c(hasil_forecast$`Hi 80`,rev(hasil_forecast$`Lo 80`)), 
   col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecast$`Point Forecast`)
lines(fitted(fit),col='blue')

out <- (hasil_forecast$`Point Forecast`< hasil_forecast$`Lo 80` | hasil_forecast$`Point Forecast` > hasil_forecast$`Hi 80`)
plot(hasil_forecast$`Point Forecast`, type="n", ylim=range(hasil_forecast$`Lo 95`,hasil_forecast$`Hi 95`))
polygon(c(time(hasil_forecast$`Point Forecast`),rev(time(hasil_forecast$`Point Forecast`))), c(hasil_forecast$`Hi 95`,rev(hasil_forecast$`Lo 95`)), 
   col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecast$`Point Forecast`)
lines(fitted(fit),col='blue')

out <- (hasil_forecast$`Point Forecast`< hasil_forecast$`Lo 95` | hasil_forecast$`Point Forecast` > hasil_forecast$`Hi 95`)

Berdasarkan kedua plot yaitu plot interval prediksi 80% dan plot interval prediksi 95% menunjukkan hasil yang hampir sama hanya berbeda di rentang selangnya. Interval prediksi 95% memiliki rentang yang lebih luas dari pada interval prediksi 80%.