Tugas Pertemuan 2 - Metode Pemulusan Deret Waktu

Data PM2.5 (Anggota 2) - STA1341

1 Pendahuluan

Data yang digunakan pada tugas ini adalah data kualitas udara harian, tepatnya variabel PM2.5 (particulate matter 2.5 mikron, ug/m3), yang diambil dari dataset Air Quality Historical. Data asli memiliki panjang 780 periode (1 Januari 2024 - 18 Februari 2026) dan dibagi rata secara berurutan (bukan diacak, agar pola deret waktu tidak rusak) ke 5 anggota kelompok, masing-masing mendapat 156 observasi.

Bagian data yang digunakan pada laporan ini adalah bagian milik Anggota 2, yaitu periode 5 Juni 2024 - 7 November 2024 (156 hari).

Tugas ini akan menerapkan dan membandingkan empat metode pemulusan yang telah dipelajari: Single Moving Average (SMA), Double Moving Average (DMA), Single Exponential Smoothing (SES), dan Double Exponential Smoothing (DES), kemudian menentukan metode terbaik berdasarkan metrik akurasi SSE, MSE, RMSE, dan MAPE pada data uji.

2 Library / Packages

Package R yang digunakan pada laporan ini adalah forecast, TTR, readxl, dan ggplot2. Jika package tersebut belum ada, silakan install terlebih dahulu.

# install.packages("forecast")
# install.packages("TTR")
# install.packages("readxl")
# install.packages("ggplot2")

library(forecast)
library(TTR)
library(readxl)
library(ggplot2)

3 Impor Data

Pastikan file PM25_Anggota_2.xlsx berada pada folder kerja (working directory) yang sama dengan file Rmd ini sebelum melakukan knit.

data <- read_excel("C:/Users/oktav/Downloads/PM25_Anggota_2.xlsx", skip = 3)
colnames(data) <- c("t", "tanggal", "pm25")
data$tanggal <- as.Date(data$tanggal, format = "%Y-%m-%d")

head(data)

## # A tibble: 6 × 3
##       t tanggal     pm25
##   <dbl> <date>     <dbl>
## 1     1 2024-06-05 104. 
## 2     2 2024-06-06  87.1
## 3     3 2024-06-07  72.5
## 4     4 2024-06-08  67.8
## 5     5 2024-06-09  73.4
## 6     6 2024-06-10  83.2

tail(data)

## # A tibble: 6 × 3
##       t tanggal     pm25
##   <dbl> <date>     <dbl>
## 1   151 2024-11-02  69.8
## 2   152 2024-11-03  87.9
## 3   153 2024-11-04  89.1
## 4   154 2024-11-05  86.1
## 5   155 2024-11-06  40.8
## 6   156 2024-11-07  73.2

4 Eksplorasi Data

str(data)

## tibble [156 × 3] (S3: tbl_df/tbl/data.frame)
##  $ t      : num [1:156] 1 2 3 4 5 6 7 8 9 10 ...
##  $ tanggal: Date[1:156], format: "2024-06-05" "2024-06-06" ...
##  $ pm25   : num [1:156] 104.5 87.1 72.5 67.8 73.4 ...

dim(data)

## [1] 156   3

summary(data$pm25)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   26.48   61.80   70.16   70.43   80.40  108.38

Mengubah data menjadi objek deret waktu (ts) dengan fungsi ts().

data.ts <- ts(data$pm25)
data.ts

## Time Series:
## Start = 1 
## End = 156 
## Frequency = 1 
##   [1] 104.4583  87.0750  72.5083  67.8042  73.4083  83.2292  66.2083  50.4083
##   [9]  51.1500  52.0667  70.9083  92.8833  88.7125  83.3250  63.7708  64.5375
##  [17]  74.1917 108.1250  84.3167  64.3333  61.1417  53.0625  41.9708  40.8500
##  [25]  43.9292  77.3875  90.9625  64.3458  87.3583  96.3000  45.9167  26.4750
##  [33]  30.7500  44.4833  46.6542  71.8542 108.3750 108.3792 101.2208  90.3750
##  [41]  76.4000  66.8333  60.0833  77.6042  69.0542  70.5875  86.0417  81.1208
##  [49]  88.4417  79.9042  81.2375  78.0042  64.0792  80.3583  83.6167  70.5250
##  [57]  62.1750  81.5542  88.4792  62.6417  72.5417 100.5917  84.4583  83.8708
##  [65]  76.8167  67.8125  69.1417  79.9250  88.5000  61.9917  58.3125  68.4250
##  [73]  78.3833  51.9833  63.2875  54.8542  57.7000  56.5917  50.0958  70.1375
##  [81]  94.3708  74.5458  62.3125  77.8000  77.9458  63.2708  76.6208  62.2417
##  [89]  60.3917  72.9208  87.6792  72.1083  59.0458  48.4542  49.1333  53.3042
##  [97]  52.0375  63.7458  46.5167  74.9625  83.3042  93.6208  88.0083  79.1583
## [105]  76.4875  58.9583  62.7750  49.5000  57.2458  71.0625  76.6875  80.5292
## [113]  73.2250  63.7208  60.9333  64.1542  76.4042  66.8083  71.3792  70.0375
## [121]  64.3875  75.0167  85.8000  65.3958  64.6708  66.7792  75.2833  70.1875
## [129]  60.1958  61.2375  57.5000  96.3000  89.7833  83.4917  75.6208  71.7708
## [137]  77.1333  42.6667  62.4708  67.2875  56.2333  64.9792  62.2750  53.8458
## [145]  79.7250  65.6042  52.4000  68.1667  66.9458  64.2833  69.7833  87.8958
## [153]  89.0917  86.1125  40.7583  73.1750

Membuat plot data deret waktu secara keseluruhan.

ts.plot(data.ts, xlab = "Periode (hari)", ylab = "PM2.5 (ug/m3)",
        main = "Plot Data Deret Waktu PM2.5 - Anggota 2")
points(data.ts)

Interpretasi awal: dari plot terlihat data PM2.5 berfluktuasi di sekitar nilai rataan tertentu tanpa menunjukkan pola tren naik/turun yang tegas dalam jangka panjang, sehingga secara visual data cenderung stasioner/konstan. Pemeriksaan pola musiman mingguan (hari kerja vs akhir pekan) juga dilakukan sebagai berikut.

data$hari <- weekdays(data$tanggal)
rata_hari <- aggregate(pm25 ~ hari, data = data, FUN = mean)
rata_hari[order(rata_hari$pm25), ]

##        hari     pm25
## 7 Wednesday 67.59021
## 1    Friday 68.09716
## 3  Saturday 68.43011
## 5  Thursday 70.61793
## 6   Tuesday 71.43675
## 4    Sunday 72.18977
## 2    Monday 74.76666

Selisih rataan PM2.5 antar hari dalam seminggu relatif kecil dibanding rataan keseluruhan data, dan tidak membentuk pola musiman yang konsisten sehingga metode Holt-Winters (musiman) tidak digunakan pada data ini; empat metode non-musiman (SMA, DMA, SES, DES) dinilai lebih relevan untuk dibandingkan.

5 Pembagian Data Latih dan Data Uji

Pembagian data latih dan data uji dilakukan dengan perbandingan 80% data latih dan 20% data uji.

n <- nrow(data)
n_train <- floor(0.8 * n)
n_test  <- n - n_train

training <- data[1:n_train, ]
testing  <- data[(n_train + 1):n, ]

train.ts <- ts(training$pm25)
test.ts  <- ts(testing$pm25)

n; n_train; n_test

## [1] 156

## [1] 124

## [1] 32

ggplot() +
  geom_line(data = training, aes(x = t, y = pm25, col = "Data Latih")) +
  geom_line(data = testing, aes(x = t, y = pm25, col = "Data Uji")) +
  labs(x = "Periode (hari)", y = "PM2.5 (ug/m3)", color = "Keterangan") +
  scale_colour_manual(name = "Keterangan:",
                       breaks = c("Data Latih", "Data Uji"),
                       values = c("blue", "red")) +
  theme_bw() +
  theme(legend.position = "bottom")

6 Single Moving Average (SMA) & Double Moving Average (DMA)

6.1 Pemilihan Orde (m) Terbaik

Sebelum menerapkan SMA dan DMA, dicari orde m yang menghasilkan MSE data latih terkecil, dengan mencoba beberapa kandidat m (3 sampai 10).

m_candidates <- 3:10
mse_m <- sapply(m_candidates, function(m) {
  sma_m <- SMA(train.ts, n = m)
  err <- train.ts - c(NA, sma_m[-length(sma_m)])
  mean(err[!is.na(err)]^2)
})

hasil_m <- data.frame(m = m_candidates, MSE_latih = mse_m)
hasil_m

##    m MSE_latih
## 1  3  287.0107
## 2  4  310.1260
## 3  5  333.5234
## 4  6  342.3422
## 5  7  333.2544
## 6  8  314.4147
## 7  9  295.0424
## 8 10  275.2525

m_terbaik <- hasil_m$m[which.min(hasil_m$MSE_latih)]
m_terbaik

## [1] 10

Orde m dengan MSE data latih terkecil akan digunakan pada pemulusan SMA dan DMA berikut ini (m = 10).

6.2 Single Moving Average (SMA)

data.sma <- SMA(train.ts, n = m_terbaik)
data.ramal.sma <- c(NA, data.sma)

# forecast n_test periode ke depan (SMA -> nilai ramalan konstan)
ramalan_sma_test <- rep(data.ramal.sma[length(data.ramal.sma)], n_test)

data.gab.sma <- cbind(
  aktual   = c(train.ts, test.ts),
  pemulusan = c(data.sma, rep(NA, n_test)),
  ramalan   = c(data.ramal.sma, ramalan_sma_test[-1])
)
tail(data.gab.sma, 10)

##         aktual pemulusan  ramalan
## [147,] 52.4000        NA 70.03167
## [148,] 68.1667        NA 70.03167
## [149,] 66.9458        NA 70.03167
## [150,] 64.2833        NA 70.03167
## [151,] 69.7833        NA 70.03167
## [152,] 87.8958        NA 70.03167
## [153,] 89.0917        NA 70.03167
## [154,] 86.1125        NA 70.03167
## [155,] 40.7583        NA 70.03167
## [156,] 73.1750        NA 70.03167

6.2.1 Akurasi SMA - Data Latih

error_train.sma <- train.ts - data.ramal.sma[1:length(train.ts)]
idx <- (m_terbaik + 1):length(train.ts)

SSE_train.sma  <- sum(error_train.sma[idx]^2)
MSE_train.sma  <- mean(error_train.sma[idx]^2)
RMSE_train.sma <- sqrt(MSE_train.sma)
MAPE_train.sma <- mean(abs(error_train.sma[idx] / train.ts[idx]) * 100)

akurasi_train.sma <- data.frame(SSE = SSE_train.sma, MSE = MSE_train.sma,
                                 RMSE = RMSE_train.sma, MAPE = MAPE_train.sma)
akurasi_train.sma

##        SSE      MSE     RMSE     MAPE
## 1 31378.79 275.2525 16.59074 20.18208

6.2.2 Akurasi SMA - Data Uji

error_test.sma <- test.ts - ramalan_sma_test

SSE_test.sma  <- sum(error_test.sma^2)
MSE_test.sma  <- mean(error_test.sma^2)
RMSE_test.sma <- sqrt(MSE_test.sma)
MAPE_test.sma <- mean(abs(error_test.sma / test.ts) * 100)

akurasi_test.sma <- data.frame(SSE = SSE_test.sma, MSE = MSE_test.sma,
                                RMSE = RMSE_test.sma, MAPE = MAPE_test.sma)
akurasi_test.sma

##        SSE      MSE     RMSE     MAPE
## 1 5374.037 167.9387 12.95911 16.19147

6.3 Double Moving Average (DMA)

DMA dilakukan dengan merata-ratakan hasil SMA sekali lagi (moving average dari moving average), untuk mengoreksi bias jika terdapat tren.

Sp  <- SMA(train.ts, n = m_terbaik)          # M't
Spp <- SMA(Sp, n = m_terbaik)                # M''t

at <- 2 * Sp - Spp
bt <- (2 / (m_terbaik - 1)) * (Sp - Spp)

# forecast 1 periode ke depan pada data latih
ramalan_dma_train <- c(NA, at[-length(at)] + bt[-length(bt)])

data.gab.dma <- cbind(aktual = train.ts, Sp = Sp, Spp = Spp,
                       at = at, bt = bt, ramalan = ramalan_dma_train)
tail(data.gab.dma, 10)

## Time Series:
## Start = 115 
## End = 124 
## Frequency = 1 
##      aktual       Sp      Spp       at          bt  ramalan
## 115 60.9333 65.46374 70.26145 60.66603 -1.06615733 62.36571
## 116 64.1542 65.98333 69.69178 62.27488 -0.82410044 59.59987
## 117 76.4042 67.34625 69.15103 65.54147 -0.40106289 61.45078
## 118 66.8083 69.07708 68.92583 69.22833  0.03361222 65.14040
## 119 71.3792 70.49042 68.73466 72.24618  0.39016889 69.26195
## 120 70.0375 70.38792 68.57225 72.20359  0.40348333 72.63635
## 121 64.3875 69.15792 68.35300 69.96284  0.17887178 72.60708
## 122 75.0167 68.60667 68.20954 69.00380  0.08825111 70.14171
## 123 85.8000 69.86417 68.33967 71.38867  0.33877867 69.09205
## 124 65.3958 70.03167 68.64092 71.42242  0.30905622 71.72745

6.3.1 Akurasi DMA - Data Latih

error_train.dma <- train.ts - ramalan_dma_train
idx2 <- which(!is.na(error_train.dma))

SSE_train.dma  <- sum(error_train.dma[idx2]^2)
MSE_train.dma  <- mean(error_train.dma[idx2]^2)
RMSE_train.dma <- sqrt(MSE_train.dma)
MAPE_train.dma <- mean(abs(error_train.dma[idx2] / train.ts[idx2]) * 100)

akurasi_train.dma <- data.frame(SSE = SSE_train.dma, MSE = MSE_train.dma,
                                 RMSE = RMSE_train.dma, MAPE = MAPE_train.dma)
akurasi_train.dma

##        SSE      MSE     RMSE     MAPE
## 1 35215.39 335.3847 18.31351 22.01043

6.3.2 Akurasi DMA - Data Uji

Ramalan untuk h periode ke depan pada DMA menggunakan persamaan \(F_{n+m} = a_n + b_n \cdot m\), dengan \(a_n, b_n\) diambil dari nilai terakhir pada data latih.

at_n <- at[length(at)]
bt_n <- bt[length(bt)]
m_ahead <- 1:n_test
ramalan_dma_test <- at_n + bt_n * m_ahead

error_test.dma <- test.ts - ramalan_dma_test

SSE_test.dma  <- sum(error_test.dma^2)
MSE_test.dma  <- mean(error_test.dma^2)
RMSE_test.dma <- sqrt(MSE_test.dma)
MAPE_test.dma <- mean(abs(error_test.dma / test.ts) * 100)

akurasi_test.dma <- data.frame(SSE = SSE_test.dma, MSE = MSE_test.dma,
                                RMSE = RMSE_test.dma, MAPE = MAPE_test.dma)
akurasi_test.dma

##        SSE      MSE     RMSE     MAPE
## 1 7503.033 234.4698 15.31241 20.88182

7 Single Exponential Smoothing (SES) & Double Exponential Smoothing (DES)

7.1 Single Exponential Smoothing (SES)

SES cocok digunakan untuk data dengan pola stasioner/konstan. Nilai parameter alpha dioptimalkan dengan alpha = NULL sehingga R mencari alpha yang meminimumkan SSE data latih.

ses.opt <- ses(train.ts, h = n_test, alpha = NULL)
ses.opt$model

## Simple exponential smoothing 
## 
## Call:
## ses(y = train.ts, h = n_test, alpha = NULL)
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 78.912 
## 
##   sigma:  14.7142
## 
##      AIC     AICc      BIC 
## 1268.525 1268.725 1276.986

autoplot(ses.opt) +
  autolayer(fitted(ses.opt), series = "Fitted") +
  ylab("PM2.5 (ug/m3)") + xlab("Periode") +
  ggtitle("SES - Data Latih dan Hasil Peramalan")

7.1.1 Akurasi SES - Data Latih

Perhitungan akurasi data latih dilakukan langsung dari nilai residual (ses.opt$residuals = data aktual - nilai fitted) agar tidak bergantung pada nama field internal objek ets yang dapat berbeda antar versi package forecast.

resid_train.ses <- ses.opt$residuals

SSE_train.ses  <- sum(resid_train.ses^2)
MSE_train.ses  <- mean(resid_train.ses^2)
RMSE_train.ses <- sqrt(MSE_train.ses)
MAPE_train.ses <- mean(abs(resid_train.ses / train.ts) * 100)

akurasi_train.ses <- data.frame(SSE = SSE_train.ses, MSE = MSE_train.ses,
                                 RMSE = RMSE_train.ses, MAPE = MAPE_train.ses)
akurasi_train.ses

##        SSE      MSE     RMSE     MAPE
## 1 26414.07 213.0167 14.59509 17.03032

7.1.2 Akurasi SES - Data Uji

ramalan_ses_test <- as.numeric(ses.opt$mean)
error_test.ses <- test.ts - ramalan_ses_test

SSE_test.ses  <- sum(error_test.ses^2)
MSE_test.ses  <- mean(error_test.ses^2)
RMSE_test.ses <- sqrt(MSE_test.ses)
MAPE_test.ses <- mean(abs(error_test.ses / test.ts) * 100)

akurasi_test.ses <- data.frame(SSE = SSE_test.ses, MSE = MSE_test.ses,
                                RMSE = RMSE_test.ses, MAPE = MAPE_test.ses)
akurasi_test.ses

##       SSE     MSE     RMSE     MAPE
## 1 5714.91 178.591 13.36379 14.99509

7.2 Double Exponential Smoothing (DES) - Metode Holt

DES (metode Holt, dua parameter) digunakan untuk memeriksa apakah memasukkan komponen tren dapat memperbaiki hasil peramalan, meskipun dari eksplorasi awal data cenderung tidak menunjukkan tren yang kuat.

des.opt <- HoltWinters(train.ts, gamma = FALSE)
des.opt

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

ramalan.des <- forecast(des.opt, h = n_test)
autoplot(ramalan.des) + ylab("PM2.5 (ug/m3)") + xlab("Periode") +
  ggtitle("DES (Holt) - Hasil Peramalan Data Uji")

7.2.1 Akurasi DES - Data Latih

Objek HoltWinters menyimpan langsung nilai SSE (sum of squared errors) pada elemen $SSE. Nilai MSE dihitung dengan membagi SSE dengan banyaknya titik fitted (karena beberapa observasi awal hilang akibat inisialisasi level & tren).

SSE_train.des <- des.opt$SSE
fitted.des <- fitted(des.opt)[, 1]
actual.des <- train.ts[(length(train.ts) - length(fitted.des) + 1):length(train.ts)]

MSE_train.des  <- SSE_train.des / length(fitted.des)
RMSE_train.des <- sqrt(MSE_train.des)
MAPE_train.des <- mean(abs((actual.des - fitted.des) / actual.des) * 100)

akurasi_train.des <- data.frame(SSE = SSE_train.des, MSE = MSE_train.des,
                                 RMSE = RMSE_train.des, MAPE = MAPE_train.des)
akurasi_train.des

##        SSE      MSE     RMSE     MAPE
## 1 29098.58 238.5129 15.44386 18.12265

7.2.2 Akurasi DES - Data Uji

ramalan_des_test <- as.numeric(ramalan.des$mean)
error_test.des <- test.ts - ramalan_des_test

SSE_test.des  <- sum(error_test.des^2)
MSE_test.des  <- mean(error_test.des^2)
RMSE_test.des <- sqrt(MSE_test.des)
MAPE_test.des <- mean(abs(error_test.des / test.ts) * 100)

akurasi_test.des <- data.frame(SSE = SSE_test.des, MSE = MSE_test.des,
                                RMSE = RMSE_test.des, MAPE = MAPE_test.des)
akurasi_test.des

##        SSE      MSE     RMSE     MAPE
## 1 8796.957 274.9049 16.58026 17.86193

8 Perbandingan Metode & Pemilihan Model Terbaik

Perbandingan dilakukan berdasarkan metrik akurasi pada data uji, karena performa pada data uji mencerminkan kemampuan generalisasi model terhadap data baru yang belum “dilihat” saat pemulusan.

perbandingan <- data.frame(
  Metode = c("SMA", "DMA", "SES", "DES"),
  SSE  = c(SSE_test.sma, SSE_test.dma, SSE_test.ses, SSE_test.des),
  MSE  = c(MSE_test.sma, MSE_test.dma, MSE_test.ses, MSE_test.des),
  RMSE = c(RMSE_test.sma, RMSE_test.dma, RMSE_test.ses, RMSE_test.des),
  MAPE = c(MAPE_test.sma, MAPE_test.dma, MAPE_test.ses, MAPE_test.des)
)
perbandingan[order(perbandingan$MAPE), ]

##   Metode      SSE      MSE     RMSE     MAPE
## 3    SES 5714.910 178.5910 13.36379 14.99509
## 1    SMA 5374.037 167.9387 12.95911 16.19147
## 4    DES 8796.957 274.9049 16.58026 17.86193
## 2    DMA 7503.033 234.4698 15.31241 20.88182

metode_terbaik <- perbandingan$Metode[which.min(perbandingan$MAPE)]
metode_terbaik

## [1] "SES"

ggplot(perbandingan, aes(x = reorder(Metode, MAPE), y = MAPE, fill = Metode)) +
  geom_col() +
  geom_text(aes(label = round(MAPE, 2)), vjust = -0.3) +
  labs(x = "Metode", y = "MAPE Data Uji (%)",
       title = "Perbandingan MAPE Data Uji Antar Metode Pemulusan") +
  theme_bw() + theme(legend.position = "none")

9 Kesimpulan

Berdasarkan hasil perbandingan metrik akurasi pada data uji, metode dengan nilai MAPE terkecil adalah SES, sehingga metode ini dipilih sebagai metode pemulusan terbaik untuk data PM2.5 pada bagian Anggota 2.

Beberapa hal yang mendasari hasil ini:

• Dari eksplorasi awal, data PM2.5 pada rentang 5 Juni - 7 November 2024 cenderung berfluktuasi di sekitar rataan tertentu tanpa tren jangka panjang yang tegas dan tanpa pola musiman mingguan yang konsisten, sehingga metode non-musiman (SMA, DMA, SES, DES) sudah tepat dipilih dan metode Holt-Winters musiman tidak diperlukan.
• Metode SMA dan SES dirancang untuk data dengan pola konstan/stasioner, sedangkan DMA dan DES dirancang untuk data yang mengandung tren. Karena data ini tidak menunjukkan tren yang kuat, secara teori SMA/SES diharapkan memberi hasil yang sebanding atau lebih baik dibanding DMA/DES - hal ini dapat dicek langsung pada tabel perbandingan di atas.
• Nilai MAPE SES pada data uji adalah 15%, yang menurut kriteria umum (Lewis, 1982) dapat dikategorikan sebagai akurasi baik.