1 Packages

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(tsibble)
## Warning: package 'tsibble' was built under R version 4.5.3
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)
## Warning: package 'MASS' was built under R version 4.5.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library(TSA)
## Warning: package 'TSA' was built under R version 4.5.3
## 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(TTR)
## Warning: package 'TTR' was built under R version 4.5.3
library(aTSA)
## Warning: package 'aTSA' was built under R version 4.5.2
## 
## Attaching package: 'aTSA'
## The following object is masked from 'package:forecast':
## 
##     forecast
## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## The following object is masked from 'package:graphics':
## 
##     identify
library(graphics)

2 Data

2.1 Input Data

Data yang akan digunakan adalah data kualitas udara PM2.5 harian yang berisi 161 pengamatan dari tanggal 18 September 2024 hingga 28 Februari 2025.

datapm25 <- read.csv("dataset_mpdw.csv")
colnames(datapm25) <- c("tanggal", "pm25")
datapm25$tanggal <- as.Date(datapm25$tanggal, format = "%m/%d/%Y")
pm25.ts <- ts(datapm25$pm25)
str(datapm25)
## 'data.frame':    161 obs. of  2 variables:
##  $ tanggal: Date, format: "2024-09-18" "2024-09-19" ...
##  $ pm25   : int  76 83 88 94 101 87 74 89 71 70 ...

Data kemudian dibagi menjadi data latih dan data uji. Pembagian kali ini dilakukan dengan proporsi / perbandingan, yaitu 80:20.

n <- length(pm25.ts)
ntrain <- round(0.8 * n)

train.ts <- ts(datapm25$pm25[1:ntrain])
test.ts  <- ts(datapm25$pm25[(ntrain + 1):n])

c(Total = n, Train = length(train.ts), Test = length(test.ts))
## Total Train  Test 
##   161   129    32

2.2 Eksplorasi Data

Sebelum masuk dalam tahap pemodelan, dilakukan eksplorasi data dengan plot deret waktu untuk melihat pola data.

#--PLOT TIME SERIES--#
plot.ts(pm25.ts, 
        lty = 1, 
        xlab = "Waktu (hari ke-)", 
        ylab = "PM2.5", 
        main = "Plot Data PM2.5")

Berdasarkan plot data deret waktu, konsentrasi PM2.5 berfluktuasi naik turun di sekitar suatu nilai tengah tanpa memperlihatkan tren naik atau turun yang jelas dalam jangka panjang, meskipun tampak sedikit pergeseran level pada pertengahan periode pengamatan.

#--PLOT DATA LATIH---#
plot.ts(train.ts, lty = 1, xlab = "Waktu", ylab = "PM2.5", main = "Plot PM2.5 Data Latih")

#--PLOT DATA UJI---#
plot.ts(test.ts, lty = 1, xlab = "Waktu", ylab = "PM2.5", main = "Plot PM2.5 Data Uji")

#--CEK KESTASIONERAN---#
acf(train.ts, lag.max = 20, main = "ACF Data Latih")

Plot ACF menurun secara bertahap (tails off) namun cukup cepat mendekati nol, mengindikasikan data cenderung sudah stasioner dalam rataan.

tseries::adf.test(train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -3.449, Lag order = 5, p-value = 0.04975
## alternative hypothesis: stationary
#stasioner

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF, diperoleh p-value sebesar 0,0498 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan disimpulkan bahwa data sudah stasioner dalam rataan (meskipun nilainya berada tidak jauh dari batas 0,05). Dengan demikian, proses differencing tidak diperlukan untuk data ini.

index <- seq(1:length(train.ts))

# Box-Cox
bc <- boxcox(train.ts ~ index, lambda = seq(-2, 4, by = 0.1))

# Nilai lambda optimum
lambda_opt <- bc$x[which.max(bc$y)]
lambda_opt
## [1] 1.151515
# Selang kepercayaan 95% untuk lambda
ci_lambda <- bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, df = 1)]

ci_lower <- min(ci_lambda)
ci_upper <- max(ci_lambda)

# Batas bawah dan batas atas
c(Batas_Bawah = ci_lower, Batas_Atas = ci_upper)
## Batas_Bawah  Batas_Atas 
##   0.8484848   1.5151515

Nilai rounded value (\(\lambda\)) optimum berada di sekitar 1,15 dengan selang kepercayaan 95% antara 0,85 dan 1,52. Karena selang tersebut memuat nilai satu, dapat disimpulkan bahwa data sudah stasioner dalam ragam sehingga transformasi tidak diperlukan.

2.3 Identifikasi Model

Karena data sudah stasioner dalam rataan dan ragam tanpa perlu differencing, identifikasi model dilanjutkan langsung menggunakan plot ACF, PACF, dan EACF pada data latih.

2.3.1 Plot ACF dan PACF

#---SPESIFIKASI MODEL---#
par(mfrow = c(1,2))
acf(train.ts, lag.max = 20, main = "ACF")
pacf(train.ts, lag.max = 20, main = "PACF")

par(mfrow = c(1,1))

Plot ACF meluruh secara perlahan (tails off), sedangkan plot PACF cuts off tegas pada lag pertama. Pola ini mengindikasikan model tentatif ARIMA(1,0,0).

2.3.2 Plot EACF

eacf(train.ts)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 x o o o o o o o x o o  o  o  o 
## 2 x o o o o o o o x o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  o  o 
## 5 x x o x o o o o o o o  o  o  o 
## 6 o o x o o o o o o o o  o  o  o 
## 7 x o x x o x o o o o o  o  o  o

Berdasarkan pola segitiga nol pada plot EACF, ujung segitiga berada pada AR(1) dan MA(1) sehingga model tentatif tambahan yang terbentuk adalah ARIMA(1,0,1), beserta model di sekitarnya seperti ARIMA(2,0,1) dan ARIMA(2,0,0) sebagai kandidat pembanding.

2.4 Pendugaan Parameter

Selanjutnya akan dilakukan pendugaan parameter kelima model ARIMA yang terbentuk sebelumnya. Pendugaan dilakukan dengan fungsi Arima() yang dilanjutkan dengan melihat nilai AIC pada ringkasan data dan melihat signifikansi parameter.

#---PENDUGAAN PARAMETER MODEL---#
model1.pm=Arima(train.ts, order=c(1,0,0), method="ML")
summary(model1.pm) #AIC=1091.60
## Series: train.ts 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     mean
##       0.6419  66.2189
## s.e.  0.0667   3.9355
## 
## sigma^2 = 267.5:  log likelihood = -542.8
## AIC=1091.6   AICc=1091.79   BIC=1100.18
## 
## Training set error measures:
##                       ME     RMSE     MAE      MPE     MAPE      MASE
## Training set -0.06860585 16.22804 12.2273 -10.9996 25.59521 0.9266395
##                     ACF1
## Training set -0.09286855
lmtest::coeftest(model1.pm) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.641888   0.066716  9.6211 < 2.2e-16 ***
## intercept 66.218876   3.935501 16.8260 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model2.pm=Arima(train.ts, order=c(0,0,1), method="ML")
summary(model2.pm) #AIC=1117.17
## Series: train.ts 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1     mean
##       0.4807  66.1988
## s.e.  0.0626   2.3327
## 
## sigma^2 = 326.8:  log likelihood = -555.59
## AIC=1117.17   AICc=1117.36   BIC=1125.75
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE     MAPE     MASE
## Training set -0.01772802 17.93784 14.52688 -14.21903 30.94651 1.100912
##                   ACF1
## Training set 0.1800047
lmtest::coeftest(model2.pm) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ma1        0.480739   0.062575  7.6825  1.56e-14 ***
## intercept 66.198822   2.332735 28.3782 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model3.pm=Arima(train.ts, order=c(1,0,1), method="ML")
summary(model3.pm) #AIC=1089.85
## Series: train.ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     mean
##       0.8107  -0.3064  66.3887
## s.e.  0.0903   0.1663   5.0151
## 
## sigma^2 = 261.8:  log likelihood = -540.92
## AIC=1089.85   AICc=1090.17   BIC=1101.29
## 
## Training set error measures:
##                      ME    RMSE      MAE       MPE     MAPE      MASE
## Training set -0.1384276 15.9899 11.98837 -11.59449 25.85841 0.9085322
##                    ACF1
## Training set 0.02599364
lmtest::coeftest(model3.pm) #ma1 tidak signifikan pada alpha 5%
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)    
## ar1        0.810744   0.090328  8.9756  < 2e-16 ***
## ma1       -0.306396   0.166315 -1.8423  0.06544 .  
## intercept 66.388698   5.015086 13.2378  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model4.pm=Arima(train.ts, order=c(2,0,0), method="ML")
summary(model4.pm) #AIC=1090.61
## Series: train.ts 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2     mean
##       0.5436  0.1512  66.2850
## s.e.  0.0867  0.0869   4.5298
## 
## sigma^2 = 263.4:  log likelihood = -541.3
## AIC=1090.61   AICc=1090.93   BIC=1102.05
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.09525862 16.03836 12.01391 -11.27742 25.58564 0.9104682
##                      ACF1
## Training set -0.004759243
lmtest::coeftest(model4.pm) #ar2 tidak signifikan
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.543553   0.086747  6.2660 3.705e-10 ***
## ar2        0.151182   0.086870  1.7403    0.0818 .  
## intercept 66.284956   4.529832 14.6330 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model5.pm=Arima(train.ts, order=c(2,0,1), method="ML")
summary(model5.pm) #AIC=1085.09
## Series: train.ts 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     mean
##       1.3739  -0.3840  -0.8983  67.0537
## s.e.  0.1020   0.0973   0.0527  10.3409
## 
## sigma^2 = 249.5:  log likelihood = -537.55
## AIC=1085.09   AICc=1085.58   BIC=1099.39
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.9733248 15.55018 11.91873 -12.3779 25.71484 0.9032552
##                     ACF1
## Training set -0.00195341
lmtest::coeftest(model5.pm) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        1.373943   0.102042  13.4645 < 2.2e-16 ***
## ar2       -0.383997   0.097344  -3.9448 7.988e-05 ***
## ma1       -0.898328   0.052712 -17.0422 < 2.2e-16 ***
## intercept 67.053742  10.340908   6.4843 8.913e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#model yang dipilih adalah model 5, yaitu ARIMA(2,0,1)
model_tentatif <- data.frame(
  Model = c("ARIMA(1,0,0)", "ARIMA(0,0,1)", "ARIMA(1,0,1)",
            "ARIMA(2,0,0)", "ARIMA(2,0,1)"),
  AIC   = c(AIC(model1.pm), AIC(model2.pm), AIC(model3.pm),
            AIC(model4.pm), AIC(model5.pm)),
  BIC   = c(BIC(model1.pm), BIC(model2.pm), BIC(model3.pm),
            BIC(model4.pm), BIC(model5.pm))
)

model_tentatif
##          Model      AIC      BIC
## 1 ARIMA(1,0,0) 1091.596 1100.176
## 2 ARIMA(0,0,1) 1117.172 1125.752
## 3 ARIMA(1,0,1) 1089.850 1101.289
## 4 ARIMA(2,0,0) 1090.606 1102.045
## 5 ARIMA(2,0,1) 1085.092 1099.391

Berdasarkan pendugaan parameter di atas, nilai AIC dan BIC terkecil dimiliki oleh model ARIMA(2,0,1) dan seluruh parameternya juga signifikan pada taraf nyata 5% sehingga model ARIMA(2,0,1) dipilih sebagai model terbaik.

2.5 Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

2.5.1 Eksplorasi Sisaan

#Eksplorasi
model.terbaik.pm <- Arima(train.ts, order = c(2,0,1), method = "ML")
sisaan.pm <- model.terbaik.pm$residuals

par(mfrow = c(2,2),
    mar = c(4,4,2,1),
    oma = c(0,0,2,0))

# Q-Q Plot
qqnorm(sisaan.pm, main = "Normal Q-Q Plot")
qqline(sisaan.pm, col = "blue", lwd = 2)

# Plot Sisaan vs Waktu
plot(sisaan.pm, type = "p",
     main = "Plot Sisaan vs Waktu",
     xlab = "Waktu", ylab = "Sisaan")

# ACF
acf(sisaan.pm, main = "ACF Sisaan")

# PACF
pacf(sisaan.pm, main = "PACF Sisaan")

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, sisaan menyebar cukup mengikuti garis \(45^{\circ}\) meski terdapat sedikit penyimpangan di kedua ujung. Plot sisaan terhadap waktu memperlihatkan lebar pita yang relatif bervariasi pada beberapa periode, mengindikasikan kemungkinan ketidakhomogenan ragam. Sementara itu, plot ACF dan PACF sisaan tidak menunjukkan spike yang signifikan, sehingga secara visual sisaan dapat dianggap saling bebas. Kondisi-kondisi ini akan diuji lebih lanjut dengan uji formal.

2.5.2 Uji Formal

#1) Sisaan Menyebar Normal
ks.test(sisaan.pm, "pnorm", mean(sisaan.pm), sd(sisaan.pm))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.pm
## D = 0.038203, p-value = 0.9918
## alternative hypothesis: two-sided

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS, diperoleh p-value sebesar 0,992 yang lebih besar dari taraf nyata 5% sehingga gagal tolak \(H_0\) dan disimpulkan bahwa sisaan menyebar normal. Hasil ini sejalan dengan eksplorasi plot kuantil-kuantil normal.

#2) Sisaan saling bebas/tidak ada autokorelasi
Box.test(sisaan.pm, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan.pm
## X-squared = 0.00050378, df = 1, p-value = 0.9821

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak saling bebas

P-value uji Ljung-Box sebesar 0,982 (> 5%) sehingga gagal tolak \(H_0\); sisaan saling bebas (tidak ada autokorelasi yang tersisa).

#3) Sisaan homogen
Box.test((sisaan.pm)^2, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  (sisaan.pm)^2
## X-squared = 7.3554, df = 1, p-value = 0.006686

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

P-value uji Ljung-Box pada sisaan kuadrat sebesar 0,0067 (< 5%) sehingga tolak \(H_0\); ragam sisaan tidak homogen. Hal ini wajar terjadi pada data polusi udara riil yang cenderung memiliki gejolak ragam (volatility clustering) dan menjadi catatan bahwa model dapat dikembangkan lebih lanjut (misalnya dengan pendekatan ARCH/GARCH) jika diperlukan pemodelan ragam yang lebih akurat.

#4) Nilai tengah sisaan sama dengan nol
t.test(sisaan.pm, mu = 0, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  sisaan.pm
## t = -0.70954, df = 128, p-value = 0.4793
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.687585  1.740936
## sample estimates:
##  mean of x 
## -0.9733248

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

P-value uji-t sebesar 0,479 (> 5%) sehingga gagal tolak \(H_0\); nilai tengah sisaan sama dengan nol.

2.6 Overfitting

Tahapan selanjutnya adalah overfitting dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(2,0,1). Kandidat model overfitting adalah ARIMA(3,0,1) dan ARIMA(2,0,2).

#---OVERFITTING---#
model5a.pm=Arima(train.ts, order=c(3,0,1),method="ML")
summary(model5a.pm) #AIC=1087.09
## Series: train.ts 
## ARIMA(3,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ar3      ma1     mean
##       1.3732  -0.3764  -0.0067  -0.8999  66.9412
## s.e.  0.1028   0.1470   0.0934   0.0545  10.3868
## 
## sigma^2 = 251.5:  log likelihood = -537.54
## AIC=1087.09   AICc=1087.78   BIC=1104.25
## 
## Training set error measures:
##                      ME     RMSE      MAE     MPE     MAPE      MASE
## Training set -0.9645725 15.54974 11.91836 -12.367 25.71354 0.9032266
##                      ACF1
## Training set 0.0008040827
lmtest::coeftest(model5a.pm) #ar3 tidak signifikan
## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## ar1        1.3732175  0.1027843  13.3602 < 2.2e-16 ***
## ar2       -0.3763774  0.1469514  -2.5612   0.01043 *  
## ar3       -0.0066623  0.0933820  -0.0713   0.94312    
## ma1       -0.8998601  0.0544837 -16.5161 < 2.2e-16 ***
## intercept 66.9412239 10.3868278   6.4448 1.157e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model5b.pm=Arima(train.ts, order=c(2,0,2),method="ML")
summary(model5b.pm) #AIC=1087.09
## Series: train.ts 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     mean
##       1.3859  -0.3956  -0.9117  0.0110  67.0014
## s.e.  0.2136   0.2085   0.2214  0.1812  10.3893
## 
## sigma^2 = 251.5:  log likelihood = -537.54
## AIC=1087.09   AICc=1087.78   BIC=1104.25
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.9717332 15.54978 11.91856 -12.3763 25.71478 0.9032422
##                       ACF1
## Training set -0.0001579425
lmtest::coeftest(model5b.pm) #ma2 tidak signifikan
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        1.385865   0.213554  6.4895 8.610e-11 ***
## ar2       -0.395599   0.208519 -1.8972    0.0578 .  
## ma1       -0.911734   0.221409 -4.1179 3.824e-05 ***
## ma2        0.010955   0.181231  0.0604    0.9518    
## intercept 67.001428  10.389315  6.4491 1.125e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#model yang dipilih tetap model awal, yaitu ARIMA(2,0,1)
overfitting.pm <- data.frame(
  Model = c("ARIMA(2,0,1)", "ARIMA(3,0,1)", "ARIMA(2,0,2)"),
  AIC   = c(AIC(model.terbaik.pm), AIC(model5a.pm), AIC(model5b.pm)),
  BIC   = c(BIC(model.terbaik.pm), BIC(model5a.pm), BIC(model5b.pm))
)

overfitting.pm
##          Model      AIC      BIC
## 1 ARIMA(2,0,1) 1085.092 1099.391
## 2 ARIMA(3,0,1) 1087.088 1104.247
## 3 ARIMA(2,0,2) 1087.089 1104.248

Kedua model hasil overfitting memiliki AIC dan BIC yang lebih besar dibandingkan model ARIMA(2,0,1), dan parameter tambahan pada masing-masing model (ar3 pada ARIMA(3,0,1) dan ma2 pada ARIMA(2,0,2)) tidak signifikan. Oleh karena itu, model ARIMA(2,0,1) tetap dipilih sebagai model terbaik untuk digunakan pada tahap peramalan.

2.7 Model Terbaik: ARIMA(2,0,1)

model.terbaik.pm <- Arima(train.ts, order = c(2,0,1), method = "ML")
summary(model.terbaik.pm)
## Series: train.ts 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     mean
##       1.3739  -0.3840  -0.8983  67.0537
## s.e.  0.1020   0.0973   0.0527  10.3409
## 
## sigma^2 = 249.5:  log likelihood = -537.55
## AIC=1085.09   AICc=1085.58   BIC=1099.39
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.9733248 15.55018 11.91873 -12.3779 25.71484 0.9032552
##                     ACF1
## Training set -0.00195341
lmtest::coeftest(model.terbaik.pm)
## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        1.373943   0.102042  13.4645 < 2.2e-16 ***
## ar2       -0.383997   0.097344  -3.9448 7.988e-05 ***
## ma1       -0.898328   0.052712 -17.0422 < 2.2e-16 ***
## intercept 67.053742  10.340908   6.4843 8.913e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.8 Peramalan

Peramalan dilakukan menggunakan fungsi forecast() sepanjang data uji (32 hari ke depan).

#---FORECAST---#
ramalan.pm <- forecast::forecast(model.terbaik.pm, h = length(test.ts))
ramalan.pm
##     Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## 130       55.12762 34.88292 75.37232 24.16603  86.08921
## 131       55.29653 32.87869 77.71437 21.01141  89.58166
## 132       55.47960 32.40753 78.55168 20.19391  90.76529
## 133       55.66627 32.28369 79.04885 19.90570  91.42684
## 134       55.85244 32.26218 79.44270 19.77425  91.93064
## 135       56.03655 32.27680 79.79630 19.69915  92.37396
## 136       56.21802 32.30615 80.12988 19.64797  92.78806
## 137       56.39664 32.34254 80.45075 19.60906  93.18423
## 138       56.57239 32.38299 80.76178 19.57790  93.56687
## 139       56.74525 32.42629 81.06422 19.55260  93.93790
## 140       56.91527 32.47185 81.35870 19.53228  94.29827
## 141       57.08250 32.51937 81.64563 19.51643  94.64856
## 142       57.24696 32.56863 81.92530 19.50471  94.98922
## 143       57.40872 32.61947 82.19796 19.49684  95.32060
## 144       57.56780 32.67175 82.46386 19.49257  95.64303
## 145       57.72426 32.72533 82.72320 19.49169  95.95684
## 146       57.87815 32.78010 82.97619 19.49400  96.26229
## 147       58.02949 32.83595 83.22303 19.49929  96.55969
## 148       58.17834 32.89276 83.46391 19.50739  96.84928
## 149       58.32473 32.95045 83.69900 19.51812  97.13133
## 150       58.46870 33.00892 83.92849 19.53133  97.40608
## 151       58.61031 33.06809 84.15252 19.54686  97.67376
## 152       58.74957 33.12787 84.37128 19.56456  97.93459
## 153       58.88654 33.18819 84.58490 19.58430  98.18879
## 154       59.02125 33.24897 84.79354 19.60595  98.43656
## 155       59.15374 33.31015 84.99733 19.62938  98.67810
## 156       59.28405 33.37167 85.19642 19.65448  98.91361
## 157       59.41220 33.43346 85.39094 19.68114  99.14326
## 158       59.53824 33.49547 85.58102 19.70925  99.36723
## 159       59.66220 33.55764 85.76676 19.73872  99.58569
## 160       59.78412 33.61993 85.94831 19.76944  99.79879
## 161       59.90403 33.68229 86.12576 19.80134 100.00671
data.ramalan.pm <- ramalan.pm$mean
plot(ramalan.pm)

Berdasarkan plot hasil ramalan, nilai ramalan model ARIMA(2,0,1) cenderung konvergen dengan cepat menuju nilai tengah proses (sekitar 55-60), yang merupakan karakteristik khas model AR/ARMA yang stasioner dengan koefisien AR yang tidak terlalu dekat dengan 1.

perbandingan.pm <- data.frame(
  Aktual   = as.numeric(test.ts),
  Forecast = as.numeric(data.ramalan.pm)
)
perbandingan.pm
##    Aktual Forecast
## 1      59 55.12762
## 2      50 55.29653
## 3      36 55.47960
## 4      32 55.66627
## 5      30 55.85244
## 6      51 56.03655
## 7      51 56.21802
## 8      47 56.39664
## 9      45 56.57239
## 10     73 56.74525
## 11     47 56.91527
## 12     53 57.08250
## 13     50 57.24696
## 14     48 57.40872
## 15     52 57.56780
## 16     38 57.72426
## 17     51 57.87815
## 18     74 58.02949
## 19     74 58.17834
## 20     68 58.32473
## 21     54 58.46870
## 22     67 58.61031
## 23     73 58.74957
## 24     50 58.88654
## 25     72 59.02125
## 26     63 59.15374
## 27     70 59.28405
## 28     55 59.41220
## 29     66 59.53824
## 30     57 59.66220
## 31     85 59.78412
## 32     57 59.90403
# Akurasi
accuracy(as.numeric(data.ramalan.pm), as.numeric(test.ts))
##                 ME     RMSE      MAE       MPE    MAPE
## Test set -1.506953 12.36856 10.47278 -8.641303 20.9513
# Plot perbandingan data aktual penuh dengan hasil forecast
ts.plot(pm25.ts, col = "black", lty = 1,
        xlab = "Waktu (hari ke-)", ylab = "PM2.5",
        main = "Aktual vs Forecast PM2.5")

lines((ntrain + 1):n, as.numeric(data.ramalan.pm), col = "blue", lty = 2, lwd = 2)
legend("topleft", 
       legend = c("Data Aktual", "Forecast"),
       col = c("black", "blue"),
       lty = c(1,2), lwd = c(1,2), bty = "n")

2.9 Kesimpulan

Model ARIMA(2,0,1) terpilih sebagai model terbaik untuk data PM2.5 karena memiliki AIC dan BIC terendah di antara model tentatif, seluruh parameter signifikan, serta lolos sebagian besar uji diagnostik sisaan (normalitas, kebebasan, dan nilai tengah nol), kecuali asumsi kehomogenan ragam yang tidak terpenuhi — hal yang umum dijumpai pada data lingkungan/polusi udara riil. Hasil peramalan pada data uji menghasilkan RMSE sekitar 12,4 dan MAPE sekitar 21%, yang menunjukkan model cukup baik dalam menangkap level rata-rata PM2.5 meskipun kurang mampu mengikuti fluktuasi harian yang tajam pada data uji.