Cabai besar merupakan salah satu komoditas hortikultura strategis di Indonesia, termasuk di Provinsi Jawa Timur. Komoditas ini memiliki peran penting dalam memenuhi kebutuhan konsumsi masyarakat serta menjaga stabilitas harga pangan. Menurut Badan Pusat Statistik (BPS), produktivitas cabai besar di Jawa Timur pada tahun 2023 mengalami penurunan sebesar 1,91 persen dibanding tahun sebelumnya, dari 9,04 ton per hektar menjadi 8,86 ton per hektar. Penurunan ini dipengaruhi oleh berbagai faktor seperti perubahan iklim, serangan hama penyakit, serta ketidakseimbangan distribusi pupuk.

Menghadapi kondisi produksi yang tidak stabil, dibutuhkan metode peramalan yang akurat untuk mendukung pengambilan keputusan dalam pengelolaan produksi dan distribusi. Salah satu model yang umum digunakan dalam meramalkan data deret waktu adalah Autoregressive Integrated Moving Average (ARIMA). Namun, model ini sangat sensitif terhadap kehadiran outlier atau data ekstrem yang dapat mengganggu hasil analisis. Sehingga, perlu dilakukan pendeteksian dan penanganan outlier terlebih dahulu sebelum menerapkan model ARIMA agar hasil peramalan menjadi lebih akurat.

Adapun pendekatan berbasis kecerdasan buatan, seperti Neural Network (NN) mulai banyak digunakan karena kemampuannya dalam mengenali pola nonlinier dan hubungan kompleks antar data. Neural Network dinilai lebih fleksibel dalam menghadapi data yang fluktuatif dan tidak terstruktur. Oleh karena itu, untuk mengetahui metode mana yang lebih sesuai dalam menangani dinamika produksi, dilakukan perbandingan antara model ARIMA dengan deteksi outlier dan Neural Network terhadap data produksi cabai besar di Jawa Timur tahun 2017-2023. Adapun analisis dengan ARIMA Outlier dan Neural Network dapat dilihat sebagai berikut dengan taraf signifikansi sebesar 5% = 0,05.

#------------------------------------------------------#
####             ANALISIS ARIMA OUTLIER             ####
#------------------------------------------------------#
# Memanggil Package
library(MASS)
library(car)
## Loading required package: carData
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(quadprog)
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(fracdiff)
library(fUnitRoots)
library(lmtest)
library(Rssa)
## Loading required package: svd
## 
## Attaching package: 'Rssa'
## The following object is masked from 'package:stats':
## 
##     decompose
library(MLmetrics)
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
## 
##     accuracy
library(urca) 
## 
## Attaching package: 'urca'
## The following objects are masked from 'package:fUnitRoots':
## 
##     punitroot, qunitroot, unitrootTable
library(TSA)
## 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(tsoutliers)
# Memanggil Data Harga Cabai Besar di Jawa Timur Periode Januari 2017 hingga Desember 2023
data <- read.table(file.choose(), header = TRUE)
Zt <- ts(data$Zt)
# Pembagian Data Training dan Data Testing dengan Proporsi 92:8
train_size <- 77
test_size <- 7
train <- window(Zt, end = train_size)
test <- window(Zt, start = train_size + 1, end = train_size + test_size)
# Plot Pembagian Data
plot(Zt, main = "Pembagian Data Train-Test", ylab = "Zt")
abline(v = train_size + 0.5, col = "red", lty = 2)
text(train_size/2, max(Zt), paste("Train (", train_size, " obs)"), cex = 0.8)
text(train_size + test_size/2, max(Zt), paste("Test (", test_size, " obs)"), cex = 0.8)

Berdasarkan output di atas, secara visual dapat dilihat bahwa harga cabai besar di Jawa Timur periode Januari 2017 hingga Desember 2023 memiliki pola tren pada periode ke-55 hingga periode ke-80. Adapun analisis ini dilakukan dengan menggunakan data training sebanyak 77 data dan data testing sebanyak 7 data untuk meramalkan harga cabai besar di Jawa Timur 10 periode ke depan.

# Pemeriksaan Stasioneritas Data
# Transformasi Box-Cox
p <- powerTransform(train)
p
## Estimated transformation parameter 
##     train 
## -1.066255
lambda <- (train^-1.066255)
y <- lambda
y
## Time Series:
## Start = 1 
## End = 77 
## Frequency = 1 
##  [1] 8.799784e-05 4.229333e-05 5.058702e-05 7.524575e-05 8.348435e-05
##  [6] 6.285669e-05 5.273539e-05 6.845822e-05 7.047264e-05 9.571089e-05
## [11] 7.543133e-05 6.006223e-05 6.930240e-05 5.041103e-05 5.668885e-05
## [16] 8.588009e-05 8.118957e-05 9.158577e-05 8.523754e-05 9.067031e-05
## [21] 7.826661e-05 5.770002e-05 7.859706e-05 7.063607e-05 6.727248e-05
## [26] 7.377139e-05 5.555375e-05 4.722179e-05 6.503301e-05 7.065534e-05
## [31] 7.267639e-05 8.361775e-05 6.622779e-05 6.112874e-05 3.997376e-05
## [36] 8.379179e-05 8.439971e-05 7.352067e-05 5.788942e-05 4.548054e-05
## [41] 4.609606e-05 5.681494e-05 6.920041e-05 7.800885e-05 8.955825e-05
## [46] 6.532984e-05 8.209522e-05 9.483345e-05 8.225023e-05 4.146768e-05
## [51] 4.996389e-05 6.417399e-05 6.663843e-05 6.277993e-05 6.109236e-05
## [56] 5.956661e-05 5.067293e-05 4.640991e-05 3.069074e-05 3.812638e-05
## [61] 5.089660e-05 4.039683e-05 5.391428e-05 5.220469e-05 6.927456e-05
## [66] 7.329230e-05 7.304473e-05 7.027177e-05 8.353766e-05 3.596277e-05
## [71] 5.029602e-05 6.227029e-05 6.309580e-05 5.528282e-05 6.566275e-05
## [76] 5.731194e-05 6.905257e-05
ts.plot(y, ylab = "Zt Hasil Transformasi Pada Data Train", xlab = "Periode")

# Differencing 1 kali
datadiff <- diff(y, differences = 1)
datadiff
## Time Series:
## Start = 2 
## End = 77 
## Frequency = 1 
##  [1] -4.570451e-05  8.293683e-06  2.465874e-05  8.238595e-06 -2.062766e-05
##  [6] -1.012129e-05  1.572283e-05  2.014420e-06  2.523824e-05 -2.027955e-05
## [11] -1.536911e-05  9.240176e-06 -1.889138e-05  6.277824e-06  2.919124e-05
## [16] -4.690522e-06  1.039620e-05 -6.348228e-06  5.432762e-06 -1.240370e-05
## [21] -2.056659e-05  2.089705e-05 -7.960997e-06 -3.363590e-06  6.498908e-06
## [26] -1.821763e-05 -8.331962e-06  1.781121e-05  5.622336e-06  2.021049e-06
## [31]  1.094136e-05 -1.738996e-05 -5.099050e-06 -2.115498e-05  4.381803e-05
## [36]  6.079150e-07 -1.087904e-05 -1.563125e-05 -1.240888e-05  6.155246e-07
## [41]  1.071888e-05  1.238547e-05  8.808443e-06  1.154940e-05 -2.422841e-05
## [46]  1.676539e-05  1.273823e-05 -1.258322e-05 -4.078255e-05  8.496207e-06
## [51]  1.421010e-05  2.464441e-06 -3.858499e-06 -1.687570e-06 -1.525747e-06
## [56] -8.893686e-06 -4.263012e-06 -1.571918e-05  7.435640e-06  1.277022e-05
## [61] -1.049977e-05  1.351745e-05 -1.709596e-06  1.706987e-05  4.017740e-06
## [66] -2.475725e-07 -2.772955e-06  1.326589e-05 -4.757489e-05  1.433326e-05
## [71]  1.197427e-05  8.255132e-07 -7.812979e-06  1.037992e-05 -8.350804e-06
## [76]  1.174063e-05
summary(ur.df(datadiff))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -4.549e-05 -8.761e-06  2.583e-06  9.669e-06  3.799e-05 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1     -1.4756     0.1692  -8.720 6.95e-13 ***
## z.diff.lag   0.2635     0.1075   2.451   0.0167 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.476e-05 on 72 degrees of freedom
## Multiple R-squared:  0.6135, Adjusted R-squared:  0.6028 
## F-statistic: 57.15 on 2 and 72 DF,  p-value: 1.369e-15
## 
## 
## Value of test-statistic is: -8.7199 
## 
## Critical values for test statistics: 
##      1pct  5pct 10pct
## tau1 -2.6 -1.95 -1.61
adfTest(datadiff)
## Warning in adfTest(datadiff): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -8.7199
##   P VALUE:
##     0.01 
## 
## Description:
##  Tue May 13 01:58:47 2025 by user: hp

Berdasarkan analisis di atas, dapat diketahui jika data sudah stasioner dalam variansi karena secara visual data telah berfluktuasi secara konstan dan p-value sebesar 0,01 < 0,05 yang artinya menolak H0 yaitu data telah stasioner dalam rata-rata, sehingga dapat dibentuk model ARIMA sementara.

# Identifikasi Model ARIMA Sementara dengan Grafik ACF dan PACF
acf(datadiff,main="",ylab="FOK") #q=2
grid()

pacf(datadiff,main="",ylab="FOKP") #p=2
grid()

Berdasarkan gambar di atas, diketahui bahwa terjadi cut off pada grafik ACF (FOK) setelah lag 2 dan terjadi cut off pada grafik PACF (FOKP) setelah lag 2 sehingga model yang terbentuk yakni ARIMA (0,1,1), ARIMA (0,1,2), ARIMA (1,1,0), ARIMA (1,1,1), ARIMA (1,1,2), ARIMA (2,1,0), ARIMA (2,1,1), dan ARIMA (2,1,2)

# Pengujian Signifkansi Model ARIMA Sementara
uji_model_arima_sig <- function(y, p_range = 0:2, d_range = 1, q_range = 0:2) {
  hasil_all <- list()
  hasil_sig <- list()
  aic_summary <- data.frame(
    Model = character(),
    AIC = numeric(),
    Signifikan = logical(),
    Shapiro_p = numeric(),
    Shapiro_Norm = logical(),
    Ljung_p = numeric(),
    White_Noise = logical(),
    stringsAsFactors = FALSE
  )
  
  for (d in d_range) {
    for (p in p_range) {
      for (q in q_range) {
        # Lewati model ARIMA(0,1,0) dan ARIMA(0,2,0)
        if (p == 0 && q == 0) next
        
        model_label <- paste0("ARIMA(", p, ",", d, ",", q, ")")
        cat("\n====================\n")
        cat("Model:", model_label, "\n")
        cat("====================\n")
        
        fit <- tryCatch({
          arima(x = y, order = c(p, d, q))
        }, error = function(e) {
          cat("❌ Gagal membentuk model:", model_label, "\n")
          return(NULL)
        })
        
        if (!is.null(fit)) {
          df <- length(y) - (p + q)
          test <- coeftest(fit, df = df)
          signif <- all(test[, "Pr(>|t|)"] < 0.05)
          
          # Residual
          resid <- residuals(fit)
          
          # Uji Normalitas
          shap <- shapiro.test(resid)
          shapiro_p <- shap$p.value
          shapiro_norm <- shapiro_p > 0.05
          
          # Uji White Noise
          ljung <- Box.test(resid, lag = 10, type = "Ljung-Box", fitdf = p + q)
          ljung_p <- ljung$p.value
          white_noise <- ljung_p > 0.05
          
          # Cetak ringkasan
          print(fit)
          cat("\nUji Signifikansi Koefisien:\n")
          print(test)
          cat("Signifikan:", signif, "\n")
          cat("Shapiro-Wilk p-value:", shapiro_p, "-> Normal:", shapiro_norm, "\n")
          cat("Ljung-Box p-value:", ljung_p, "-> White Noise:", white_noise, "\n")
          
          # Simpan ke list
          hasil_all[[model_label]] <- list(
            model = fit,
            coeftest = test,
            AIC = AIC(fit),
            signifikan = signif,
            shapiro_p = shapiro_p,
            shapiro_norm = shapiro_norm,
            ljung_p = ljung_p,
            white_noise = white_noise
          )
          
          if (signif) {
            hasil_sig[[model_label]] <- hasil_all[[model_label]]
            cat("✅ Model ini SIGNIFIKAN\n")
          } else {
            cat("❌ Model ini TIDAK signifikan\n")
          }
          
          # Simpan ke tabel
          aic_summary <- rbind(aic_summary, data.frame(
            Model = model_label,
            AIC = AIC(fit),
            Signifikan = signif,
            Shapiro_p = round(shapiro_p, 4),
            Shapiro_Norm = shapiro_norm,
            Ljung_p = round(ljung_p, 4),
            White_Noise = white_noise
          ))
        }
      }
    }
  }
  
  cat("\n=============================\n")
  cat("RINGKASAN MODEL ARIMA\n")
  cat("=============================\n")
  print(aic_summary[order(aic_summary$AIC), ])
  
  return(list(
    semua_model = hasil_all,
    model_signifikan = hasil_sig,
    ringkasan = aic_summary
  ))
}
hasil <- uji_model_arima_sig(y) 
## 
## ====================
## Model: ARIMA(0,1,1) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##           ma1
##       -0.7780
## s.e.   0.2471
## 
## sigma^2 estimated as 2.28e-10:  log likelihood = 735.36,  aic = -1468.73
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##     Estimate Std. Error t value Pr(>|t|)   
## ma1 -0.77796    0.24715 -3.1478 0.002351 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: TRUE 
## Shapiro-Wilk p-value: 0.4224902 -> Normal: TRUE 
## Ljung-Box p-value: 0.0529048 -> White Noise: TRUE 
## ✅ Model ini SIGNIFIKAN
## 
## ====================
## Model: ARIMA(0,1,2) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##           ma1      ma2
##       -0.5065  -0.4096
## s.e.   0.1003   0.0976
## 
## sigma^2 estimated as 1.894e-10:  log likelihood = 742.02,  aic = -1480.05
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##      Estimate Std. Error t value  Pr(>|t|)    
## ma1 -0.506451   0.100256 -5.0516 2.992e-06 ***
## ma2 -0.409629   0.097568 -4.1984 7.313e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: TRUE 
## Shapiro-Wilk p-value: 0.5007526 -> Normal: TRUE 
## Ljung-Box p-value: 0.9307643 -> White Noise: TRUE 
## ✅ Model ini SIGNIFIKAN
## 
## ====================
## Model: ARIMA(1,1,0) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##           ar1
##       -0.1836
## s.e.   0.1190
## 
## sigma^2 estimated as 2.502e-10:  log likelihood = 732.28,  aic = -1462.55
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##     Estimate Std. Error t value Pr(>|t|)
## ar1 -0.18362    0.11902 -1.5428    0.127
## 
## Signifikan: FALSE 
## Shapiro-Wilk p-value: 0.03712412 -> Normal: FALSE 
## Ljung-Box p-value: 0.3309713 -> White Noise: TRUE 
## ❌ Model ini TIDAK signifikan
## 
## ====================
## Model: ARIMA(1,1,1) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##          ar1      ma1
##       0.4265  -0.9590
## s.e.  0.1147   0.0415
## 
## sigma^2 estimated as 1.917e-10:  log likelihood = 741.57,  aic = -1479.15
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##      Estimate Std. Error  t value  Pr(>|t|)    
## ar1  0.426468   0.114655   3.7196 0.0003829 ***
## ma1 -0.959028   0.041532 -23.0913 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: TRUE 
## Shapiro-Wilk p-value: 0.4247542 -> Normal: TRUE 
## Ljung-Box p-value: 0.8863643 -> White Noise: TRUE 
## ✅ Model ini SIGNIFIKAN
## 
## ====================
## Model: ARIMA(1,1,2) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.1839  -0.6529  -0.2833
## s.e.  0.2425   0.2346   0.2112
## 
## sigma^2 estimated as 1.881e-10:  log likelihood = 742.29,  aic = -1478.59
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##     Estimate Std. Error t value Pr(>|t|)   
## ar1  0.18387    0.24251  0.7582 0.450742   
## ma1 -0.65294    0.23464 -2.7827 0.006836 **
## ma2 -0.28325    0.21116 -1.3414 0.183899   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: FALSE 
## Shapiro-Wilk p-value: 0.4632257 -> Normal: TRUE 
## Ljung-Box p-value: 0.9540647 -> White Noise: TRUE 
## ❌ Model ini TIDAK signifikan
## 
## ====================
## Model: ARIMA(2,1,0) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##           ar1      ar2
##       -0.2300  -0.2866
## s.e.   0.1162   0.1152
## 
## sigma^2 estimated as 2.309e-10:  log likelihood = 735.23,  aic = -1466.47
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##     Estimate Std. Error t value Pr(>|t|)  
## ar1 -0.23000    0.11623 -1.9788  0.05151 .
## ar2 -0.28659    0.11522 -2.4873  0.01509 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: FALSE 
## Shapiro-Wilk p-value: 0.04261777 -> Normal: FALSE 
## Ljung-Box p-value: 0.3957052 -> White Noise: TRUE 
## ❌ Model ini TIDAK signifikan
## 
## ====================
## Model: ARIMA(2,1,1) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.4810  -0.1622  -0.9434
## s.e.  0.1212   0.1229   0.0457
## 
## sigma^2 estimated as 1.874e-10:  log likelihood = 742.42,  aic = -1478.85
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##      Estimate Std. Error  t value  Pr(>|t|)    
## ar1  0.480975   0.121195   3.9686 0.0001657 ***
## ar2 -0.162203   0.122855  -1.3203 0.1908103    
## ma1 -0.943383   0.045745 -20.6224 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Signifikan: FALSE 
## Shapiro-Wilk p-value: 0.6165737 -> Normal: TRUE 
## Ljung-Box p-value: 0.9725099 -> White Noise: TRUE 
## ❌ Model ini TIDAK signifikan
## 
## ====================
## Model: ARIMA(2,1,2) 
## ====================
## 
## Call:
## arima(x = y, order = c(p, d, q))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.5309  -0.1838  -0.9944  0.0491
## s.e.  0.7053   0.3191   0.7146  0.6853
## 
## sigma^2 estimated as 1.874e-10:  log likelihood = 742.43,  aic = -1476.85
## 
## Uji Signifikansi Koefisien:
## 
## t test of coefficients:
## 
##      Estimate Std. Error t value Pr(>|t|)
## ar1  0.530928   0.705265  0.7528   0.4540
## ar2 -0.183796   0.319144 -0.5759   0.5665
## ma1 -0.994385   0.714607 -1.3915   0.1683
## ma2  0.049054   0.685252  0.0716   0.9431
## 
## Signifikan: FALSE 
## Shapiro-Wilk p-value: 0.6421245 -> Normal: TRUE 
## Ljung-Box p-value: 0.9438976 -> White Noise: TRUE 
## ❌ Model ini TIDAK signifikan
## 
## =============================
## RINGKASAN MODEL ARIMA
## =============================
##          Model       AIC Signifikan Shapiro_p Shapiro_Norm Ljung_p White_Noise
## 2 ARIMA(0,1,2) -1478.046       TRUE    0.5008         TRUE  0.9308        TRUE
## 4 ARIMA(1,1,1) -1477.146       TRUE    0.4248         TRUE  0.8864        TRUE
## 7 ARIMA(2,1,1) -1476.845      FALSE    0.6166         TRUE  0.9725        TRUE
## 5 ARIMA(1,1,2) -1476.590      FALSE    0.4632         TRUE  0.9541        TRUE
## 8 ARIMA(2,1,2) -1474.850      FALSE    0.6421         TRUE  0.9439        TRUE
## 1 ARIMA(0,1,1) -1466.728       TRUE    0.4225         TRUE  0.0529        TRUE
## 6 ARIMA(2,1,0) -1464.469      FALSE    0.0426        FALSE  0.3957        TRUE
## 3 ARIMA(1,1,0) -1460.552      FALSE    0.0371        FALSE  0.3310        TRUE

Berdasarkan output di atas, diketahui bahwa model ARIMA (2,1,0) tidak memenuhi asumsi normalitas residual. Sehingga akan dianalisis menggunakan pendeteksian outlier.

# Model ARIMA (2,1,0) dengan Deteksi Outlier
model_arima210 <- hasil$semua_model[["ARIMA(2,1,0)"]]$model
Prediksi_TFR1 <- fitted(model_arima210)
Hasil_Prediksi1 <- ifelse(Prediksi_TFR1 > 0, Prediksi_TFR1^(1/-1.066255), NA)
# Perhitungan error tanpa NA
valid_idx <- !is.na(Hasil_Prediksi1)
Error_Prediksi1 <- train[valid_idx] - Hasil_Prediksi1[valid_idx]

# Fungsi Z_threshold
Z_threshold <- function(residual, title) {
  time <- seq_along(residual)
  z_score <- (residual - mean(residual)) / sd(residual)
  idx_out <- which(abs(z_score) > 3)
  
  lower_bound <- mean(residual) - 3 * sd(residual)
  upper_bound <- mean(residual) + 3 * sd(residual)
  
  plot(time, residual, type = "o", col = "black", 
       main = sprintf("TS Plot Residual Model ARIMA %s dengan Z-Score", title))
  grid()
  abline(h = lower_bound, col = "darkgreen", lty = 5, lwd = 1)
  abline(h = upper_bound, col = "darkgreen", lty = 5, lwd = 1)
  points(time[idx_out], residual[idx_out], col = "red", pch = 19, cex = 0.8)
  if (length(idx_out) > 0) {
    text(time[idx_out], residual[idx_out], labels = time[idx_out], pos = 2, col = "black")
  }
}
Z_threshold(Error_Prediksi1, "(2,1,0)")

# Deteksi outlier
deteksi_outlier <- function(residual, title = "") {
  hampel_out <- which(abs(residual - median(residual)) > 3 * mad(residual))
  z_out <- which(abs(scale(residual)) > 3)
  
  plot(seq_along(residual), residual, type = "o", 
       main = paste("Outlier Detection:", title),
       ylab = "Residual", xlab = "Time")
  abline(h = median(residual) + c(-3, 3) * mad(residual), col = "red", lty = 2)
  abline(h = mean(residual) + c(-3, 3) * sd(residual), col = "blue", lty = 2)
  points(hampel_out, residual[hampel_out], col = "red", pch = 19)
  points(z_out, residual[z_out], col = "blue", pch = 1)
  legend("topright", legend = c("Hampel Outlier", "Z-Score Outlier"),
         col = c("red", "blue"), pch = c(19, 1))
  
  return(list(
    hampel_outliers = hampel_out,
    zscore_outliers = z_out,
    AO = tso(ts(residual), types = "AO")$outliers
  ))
}

hasil_deteksi <- deteksi_outlier(Error_Prediksi1, "ARIMA(2,1,0)")

print(hasil_deteksi) # DIPEROLEH DENGAN ZCSORE TITIK 70
## $hampel_outliers
## [1]  2 35 50 59 70
## 
## $zscore_outliers
## [1] 70
## 
## $AO
##   type ind time  coefhat    tstat
## 1   AO   2    2 6101.924 4.069742
## 2   AO  35   35 5151.456 3.435817
## 3   AO  50   50 6046.013 4.032451
## 4   AO  59   59 6261.986 4.176497
## 5   AO  70   70 7870.258 5.249150

Berdasarkan output di atas, dapat diketahui bahwa pada model ARIMA (2,1,0) dengan pendeteksian outlier terdapat 1 outlier pada residual mode ARIMA (2,1,0) yakni pada periode ke 70.

# Model ARIMA dengan Outlier
outlier1 <- rep(0, length(y))
outlier1[70] <- 1
ARIMAAO1 <- arima(y, order = c(2,1,0), xreg = outlier1, method = "CSS")
diagnostik1 <- list(
  Model = ARIMAAO1,
  Signifikansi = coeftest(ARIMAAO1),
  Diagnostik = tsdiag(ARIMAAO1, gof.lag = 72),
  Shapiro = shapiro.test(ARIMAAO1$residuals),
  KS = suppressWarnings(ks.test(ARIMAAO1$residuals, "pnorm", mean(ARIMAAO1$residuals), sd(ARIMAAO1$residuals)))
)

print(diagnostik1)
## $Model
## 
## Call:
## arima(x = y, order = c(2, 1, 0), xreg = outlier1, method = "CSS")
## 
## Coefficients:
##           ar1      ar2   xreg
##       -0.1596  -0.2826  0e+00
## s.e.   0.1098   0.1042  1e-04
## 
## sigma^2 estimated as 1.885e-10:  part log likelihood = 743.06
## 
## $Signifikansi
## 
## z test of coefficients:
## 
##         Estimate  Std. Error z value Pr(>|z|)   
## ar1  -1.5957e-01  1.0983e-01 -1.4528 0.146270   
## ar2  -2.8258e-01  1.0425e-01 -2.7106 0.006715 **
## xreg -3.1574e-05  1.1622e-04 -0.2717 0.785871   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## $Diagnostik
## NULL
## 
## $Shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  ARIMAAO1$residuals
## W = 0.99125, p-value = 0.8816
## 
## 
## $KS
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ARIMAAO1$residuals
## D = 0.052876, p-value = 0.9825
## alternative hypothesis: two-sided
# Uji Residual
checkresiduals(model_arima210)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)
## Q* = 8.3967, df = 8, p-value = 0.3957
## 
## Model df: 2.   Total lags used: 10
checkresiduals(ARIMAAO1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)
## Q* = 8.3469, df = 8, p-value = 0.4003
## 
## Model df: 2.   Total lags used: 10

Berdasarkan pengujian diagnostik di atas, dapat diketahui bahwa model ARIMA (2,1,0) dan ARIMA (2,1,0) dengan outlier memenuhi asumsi WN karena nilai p-value berada di atas garis threshold dan memenuhi asumsi normalitas residual karena p-value lebih besar daripada 0,05. Sehingga, kedua model dapat dilanjutkan ke tahap prediksi dan peramalan.

# Forecast untuk 10 periode ke depan
n_forecast_future <- 10

# Prediksi untuk model regular (ARIMA tanpa outlier)
pred_regular <- fitted(model_arima210)^(1/-1.066255)
fore_regular <- predict(model_arima210, n.ahead = test_size + n_forecast_future)
pred_regular_full <- ts(c(pred_regular, fore_regular$pred^(1/-1.066255)),
                        start = start(Zt), frequency = frequency(Zt))

# Prediksi untuk model dengan outlier
pred_outlier <- fitted(ARIMAAO1)^(1/-1.066255)
newxreg <- matrix(0, nrow = test_size + n_forecast_future, ncol = 1)
colnames(newxreg) <- "outlier1"
fore_outlier <- predict(ARIMAAO1, n.ahead = test_size + n_forecast_future, newxreg = newxreg)
pred_outlier_full <- ts(c(pred_outlier, fore_outlier$pred^(1/-1.066255)),
                        start = start(Zt), frequency = frequency(Zt))
pred_outlier_full
## Time Series:
## Start = 1 
## End = 94 
## Frequency = 1 
##  [1]  6361.000 12646.000 10691.000  8810.162  7994.250  7371.044  8597.461
##  [8]  9084.389  8011.464  8365.558  6156.722  7724.778  8073.565  7660.984
## [15] 10646.118  8966.221  7000.929  7512.844  6147.997  6695.905  6126.279
## [22]  7062.426  8514.083  6869.156  8330.714  7877.393  7513.267  9618.275
## [29] 10108.827  8506.845  8484.674  7803.626  6851.201  8341.436  8230.754
## [36] 11983.646  6736.491  7685.973  7381.805  8643.494 10441.360 10910.146
## [43]  9896.143  8550.354  7569.293  6554.814  8338.893  6480.665  6356.564
## [50]  6902.344 10507.211  9091.045  9180.810  8803.377  8741.539  8773.994
## [57]  9069.181 10320.949 10889.723 15347.892 12907.085 11510.117 13823.330
## [64]  9931.128 11088.069  8210.319  8120.169  7683.151  7801.656 10682.499
## [71]  8290.619  9470.241  8422.887  9170.006  9670.826  8305.599  9766.020
## [78]  7932.617  8313.907  8266.521  8163.607  8192.897  8217.200  8204.971
## [85]  8200.055  8204.287  8205.001  8203.691  8203.698  8204.067  8204.006
## [92]  8203.912  8203.944  8203.966
# Plot kombinasi
par(mar = c(5,4,4,2) + 0.1)
plot(Zt, col = "black", type = "o", 
     main = "Perbandingan Prediksi ARIMA(2,1,0) dengan dan tanpa Outlier",
     xlab = "Periode", 
     ylab = "Nilai",
     xlim = c(1, length(Zt) + test_size + n_forecast_future),
     ylim = range(c(Zt, pred_regular_full, pred_outlier_full), na.rm = TRUE))
# Garis prediksi
lines(pred_regular_full, col = "blue", type = "o", pch = 2)
lines(pred_outlier_full, col = "red", type = "o", pch = 3)

# Garis pembatas
abline(v = train_size + 0.5, col = "green", lty = 2)  # Garis pemisah train-test
abline(v = (train_size + test_size) + 0.5, col = "purple", lty = 2) # Garis pemisah test-forecast

# Legenda
legend("topright", legend = c("Aktual", "ARIMA(2,1,0)", "ARIMA(2,1,0)+Outlier", "Train/Test", "Test/Forecast"),
       col = c("black", "blue", "red", "green", "purple"), lty = c(1,1,1,2,2), pch = c(1,2,3,NA,NA),cex = 0.8)

Berdasarkan grafik peramalan di atas, dapat diketahui bahwa harga cabai besar di Jawa Timur untuk 10 periode ke depan yakni mulai dari Januari 2024 hingga Oktober 2024 mengalami sedikit peningkatan pada 10 periode yakni masing-masing per bulan harga cabai besar di Jawa Timur sebesar: Januari: Rp8.200,06 Februari: Rp8.204,29 Maret: Rp8.205,00 April: Rp8.203,69 Mei: Rp8.203,70; Juni: Rp8.204,07 Juli: Rp8.204,01 Agustus: Rp8.203,91 September: Rp8.203,94 Oktober: Rp8.203.97

# Nilai Akurasi
# Evaluasi model ARIMA (2,1,0)
cat("\n============== In Sample =============")
## 
## ============== In Sample =============
cat("\nMAPE ARIMA (2,1,0) In Sample = ",MAPE(pred_regular_full[1:77],train)*100)
## 
## MAPE ARIMA (2,1,0) In Sample =  16.75486
cat("\nRMSE ARIMA (2,1,0) In Sample = ",RMSE(pred_regular_full[1:77],train))
## 
## RMSE ARIMA (2,1,0) In Sample =  2205.591
cat("\nSMAPE ARIMA (2,1,0) In Sample=",smape(train,pred_regular_full[1:77])*100)
## 
## SMAPE ARIMA (2,1,0) In Sample= 17.06214
cat("\n============== Out Sample =============")
## 
## ============== Out Sample =============
cat("\nMAPE ARIMA (2,1,0) Out Sample = ",MAPE(pred_regular_full[78:84],test)*100)
## 
## MAPE ARIMA (2,1,0) Out Sample =  48.66711
cat("\nRMSE ARIMA (2,1,0) Out Sample = ",RMSE(pred_regular_full[78:84],test))
## 
## RMSE ARIMA (2,1,0) Out Sample =  2771.043
cat("\nSMAPE ARIMA (2,1,0) Out Sample=",smape(test,pred_regular_full[78:84])*100)
## 
## SMAPE ARIMA (2,1,0) Out Sample= 36.87163
# Evaluasi model ARIMA (2,1,0) dengan outlier
cat("\n================ In Sample =============")
## 
## ================ In Sample =============
cat("\nMAPE ARIMA (2,1,0) dengan outlier In Sample = ",MAPE(pred_outlier_full[78:84],train)*100)
## 
## MAPE ARIMA (2,1,0) dengan outlier In Sample =  17.85678
cat("\nRMSE ARIMA (2,1,0) dengan outlier In Sample = ",RMSE(pred_outlier_full[78:84],train))
## 
## RMSE ARIMA (2,1,0) dengan outlier In Sample =  2350.136
cat("\nSMAPE ARIMA (2,1,0) dengan outlier In Sample = ",smape(train,pred_outlier_full[78:84])*100)
## 
## SMAPE ARIMA (2,1,0) dengan outlier In Sample =  18.8121
cat("\n================ Out Sample =============")
## 
## ================ Out Sample =============
cat("\nMAPE ARIMA (2,1,0) dengan outlier Out Sample = ",MAPE(pred_outlier_full[78:84],test)*100)
## 
## MAPE ARIMA (2,1,0) dengan outlier Out Sample =  47.6403
cat("\nRMSE ARIMA (2,1,0) dengan outlier Out Sample = ",RMSE(pred_outlier_full[78:84],test))
## 
## RMSE ARIMA (2,1,0) dengan outlier Out Sample =  2723.023
cat("\nSMAPE ARIMA (2,1,0) dengan outlier Out Sample = ",smape(test,pred_outlier_full[78:84])*100)
## 
## SMAPE ARIMA (2,1,0) dengan outlier Out Sample =  36.1875

Berdasarkan output di atas, dapat dilihat bahwa nilai akurasi data out-sample pada model ARIMA (2,1,0) dengan deteksi outlier lebih kecil dibandingkan dengan nilai akurasi pada data out-sample ARIMA (2,1,0). Akan tetapi, jika dibandingkan dengan data in-sample, nilai akurasi data in-sample pada model ARIMA (2,1,0) lebih kecil dibandingkan dengan data in-sample model ARIMA (2,1,0) dengan outlier. Hal ini menunjukkan model ARIMA (2,1,0) pada data in-sample cukup baik dalam mengenali pola historis namun pada data out-sample terjadi nilai kesalahan meningkat, yang menunjukkan bahwa model belum mampu melakukan generalisasi dengan baik terhadap data baru.

#------------------------------------------------------#
####          ANALISIS NEURAL NETWORK (NN)          ####
#------------------------------------------------------#
# Import library
library(MLmetrics)
library(Metrics)
library(neuralnet)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:neuralnet':
## 
##     compute
## The following object is masked from 'package:car':
## 
##     recode
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
import_data = read.table(file.choose(), header = T)
import_data
##       Zt
## 1   6361
## 2  12646
## 3  10691
## 4   7367
## 5   6683
## 6   8721
## 7  10282
## 8   8050
## 9   7834
## 10  5879
## 11  7350
## 12  9101
## 13  7958
## 14 10726
## 15  9608
## 16  6508
## 17  6860
## 18  6127
## 19  6554
## 20  6185
## 21  7100
## 22  9450
## 23  7072
## 24  7817
## 25  8183
## 26  7505
## 27  9792
## 28 11404
## 29  8447
## 30  7815
## 31  7611
## 32  6673
## 33  8304
## 34  8952
## 35 13333
## 36  6660
## 37  6615
## 38  7529
## 39  9421
## 40 11813
## 41 11665
## 42  9588
## 43  7969
## 44  7122
## 45  6257
## 46  8411
## 47  6789
## 48  5930
## 49  6777
## 50 12882
## 51 10816
## 52  8553
## 53  8256
## 54  8731
## 55  8957
## 56  9172
## 57 10674
## 58 11591
## 59 17083
## 60 13938
## 61 10630
## 62 13202
## 63 10071
## 64 10380
## 65  7961
## 66  7551
## 67  7575
## 68  7855
## 69  6679
## 70 14723
## 71 10749
## 72  8798
## 73  8690
## 74  9837
## 75  8371
## 76  9510
## 77  7985
## 78  7902
## 79  5881
## 80  5300
## 81  5158
## 82  4614
## 83  4466
## 84  7078
cabai_besar= ts(data = import_data$Zt)
head(cabai_besar)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]  6361 12646 10691  7367  6683  8721

Data produksi cabai besar di Jawa Timur terdiri dari 84 data observasi, yang kemudian dibagi menjadi dua bagian, yaitu data training dan data testing. Sebanyak 77 data digunakan sebagai data training untuk keperluan analisis dan pembelajaran model, sedangkan 7 data sisanya digunakan sebagai data testing untuk menguji performa dan akurasi model

jumlah_train = round(NROW(cabai_besar)*0.9167,0)#77
jumlah_test = NROW(cabai_besar)-jumlah_train#7
data_train = cabai_besar[1:jumlah_train]
data_test  = cabai_besar[(jumlah_train+1):(jumlah_test+jumlah_train)]

Pada penelitian ini, proses normalisasi data dilakukan menggunakan metode Min-Max Normalization. Metode ini merupakan salah satu teknik normalisasi yang menggunakan transformasi linier untuk mengubah data dari rentang nilai aslinya ke dalam rentang nilai baru, yaitu antara 0 hingga 1. Tujuan dari normalisasi ini adalah untuk menyamakan skala antar variabel, sehingga tidak ada variabel yang mendominasi proses pemodelan hanya karena memiliki skala yang lebih besar.

## Tahap standarisasi
max_asli = max(cabai_besar)
min_asli = min(cabai_besar)
data_standarisasi = 2*((as.numeric(cabai_besar) - min_asli)/(max_asli - min_asli))-1
data_standarisasi
##  [1] -0.699611635  0.296663232 -0.013236110 -0.540144250 -0.648569391
##  [6] -0.325513196 -0.078069272 -0.431877625 -0.466117144 -0.776016486
## [11] -0.542839027 -0.265277007 -0.446461124 -0.007688040 -0.184909249
## [16] -0.676309741 -0.620512008 -0.736704446 -0.669017992 -0.727510502
## [21] -0.582468099 -0.209954823 -0.586906555 -0.468811920 -0.410794959
## [26] -0.518269002 -0.155742253  0.099786003 -0.368946659 -0.469128953
## [31] -0.501466276 -0.650154553 -0.391614488 -0.288895934  0.405563922
## [36] -0.652215265 -0.659348498 -0.514464611 -0.214551795  0.164619165
## [41]  0.141158754 -0.188079575 -0.444717445 -0.578980740 -0.716097329
## [46] -0.374653246 -0.631766664 -0.767932155 -0.633668859  0.334073076
## [51]  0.006578426 -0.352143933 -0.399223270 -0.323928034 -0.288103353
## [56] -0.254022351 -0.015930887  0.129428549  1.000000000  0.501466276
## [61] -0.022905604  0.384798288 -0.111516208 -0.062534675 -0.445985575
## [66] -0.510977253 -0.507172862 -0.462788301 -0.649203456  0.625901561
## [71] -0.004042165 -0.313307442 -0.330427201 -0.148609020 -0.380993897
## [76] -0.200443846 -0.442181184 -0.455338036 -0.775699453 -0.867797416
## [81] -0.890306729 -0.976539589 -1.000000000 -0.585955457
head(data_standarisasi)
## [1] -0.69961164  0.29666323 -0.01323611 -0.54014425 -0.64856939 -0.32551320

Variabel Masukan dipilih berdasarkan pada lag yang signifikan pada grafik PACF. Adapun grafik dapat dilihat pada gambar berikut.

## Grafik PACF Data Standarisasi
pacf(data_standarisasi, main = "PACF Data Standarisasi", lag.max = 80)
grid()

Berdasarkan grafik PACF, diketahui bahwa terdapat 1 lag dengan nilai PACF yang signifikan. Sehingga, variabel masukan yang digunakan pada metode NN adalah \(\varepsilon_{t-1}\) sebagai \(X_{t}\). Selanjutnya, k eluaran yang digunakan merupakan data produksi cabai besar yang telah distandarisasi dan dilambangkan dengan \(Y\). Adapun data masukan dan keluaran dimulai dari \(t=2\).

## Tahap Pembentukan model NN
### Input data
y = data.frame(
  nilai_asli = data_standarisasi
) %>% 
  mutate(
    nilai_lag = dplyr::lag(nilai_asli, n = 1, default = NA)
  )
#### Splitting data
trainy = y[1:77,]
trainy = na.omit(trainy)
testy = y[78:84,]

head(trainy)
##    nilai_asli   nilai_lag
## 2  0.29666323 -0.69961164
## 3 -0.01323611  0.29666323
## 4 -0.54014425 -0.01323611
## 5 -0.64856939 -0.54014425
## 6 -0.32551320 -0.64856939
## 7 -0.07806927 -0.32551320
head(testy)
##    nilai_asli  nilai_lag
## 78 -0.4553380 -0.4421812
## 79 -0.7756995 -0.4553380
## 80 -0.8677974 -0.7756995
## 81 -0.8903067 -0.8677974
## 82 -0.9765396 -0.8903067
## 83 -1.0000000 -0.9765396

Pelatihan model NN dilakukan menggunakan algoritma backpropagation dengan variasi arsitektur satu dan dua hidden layer, untuk melihat pengaruh jumlah neuron dan lapisan tersembunyi terhadap akurasi peramalan data produksi cabai besar di Jawa Timur. Pelatihan menggunakan maksimum iterasi 10.000 dan target error 0,05. Berdasarkan proses pelatihan yang telah dilakukan, diperoleh arsitektur NN sebagai berikut:

  1. NN 1 hidden layer dengan 1 Neuron
attach(trainy)
attach(testy)
## The following objects are masked from trainy:
## 
##     nilai_asli, nilai_lag
set.seed(13)
mlp_1_1 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=1, act.fct="tanh", linear.output = TRUE, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_1_1)
  1. NN 1 hidden layer dengan 2 Neuron
set.seed(13)
mlp_1_2 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=2, act.fct="tanh", linear.output = T, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_1_2)
  1. NN 2 hidden layer dengan masing-masing 1 neuron
set.seed(23)
mlp_2_1_1 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=c(1,1), act.fct="tanh", linear.output = TRUE, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_2_1_1)
  1. NN 2 hidden layer dengan 1 neuron dan 3 neuron
set.seed(23)
mlp_2_1_3 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=c(1,3), act.fct="tanh", linear.output = TRUE, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_2_1_3)
  1. NN 2 hidden layer dengan 2 neuron dan 1 neuron
set.seed(23)
mlp_2_2_1 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=c(2,1), act.fct="tanh", linear.output = TRUE, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_2_2_1)
  1. NN 2 hidden layer dengan 2 neuron dan 3 neuron
set.seed(23)
mlp_2_2_3 = neuralnet(nilai_asli~nilai_lag, data=trainy, hidden=c(2,3), act.fct="tanh", linear.output = TRUE, learningrate = 0.005, stepmax = 10000, startweights = NULL, threshold = 0.005, algorithm = "rprop+")
plot(mlp_2_2_3)

Tahap selanjutnya adalah melakukan destandarisasi dari hasil pelatihan backpropagation

## Tahap Prediksi

### 1 *hidden layer*
#### NN(1)
predict_NN1 = ts(mlp_1_1$net.result[[1]])
kembali_NN1 = ((predict_NN1 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

#### NN(2)
predict_NN2 = ts(mlp_1_2$net.result[[1]])
kembali_NN2 = ((predict_NN2 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

### 2 *hidden layer*
#### NN(1,1)
predict_NN_1_1 = ts(mlp_2_1_1$net.result[[1]])
kembali_NN_1_1 = ((predict_NN_1_1 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

#### NN(1,3)
predict_NN_1_3 = ts(mlp_2_1_3$net.result[[1]])
kembali_NN_1_3 = ((predict_NN_1_3 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

#### NN(2,1)
predict_NN_2_1 = ts(mlp_2_2_1$net.result[[1]])
kembali_NN_2_1 = ((predict_NN_2_1 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

#### NN(2,3)
predict_NN_2_3 = ts(mlp_2_2_3$net.result[[1]])
kembali_NN_2_3 = ((predict_NN_2_3 + 1) / 2) * (max(cabai_besar) - min(cabai_besar)) + min(cabai_besar)

Tahap selanjutnya adalah melakukan pemilihan NN terbaik dengan menggunakan kriteria MAPE,RMSE, dan SMAPE.

## Pemilihan Model NN Terbaik
cat("================ MAPE =============")
## ================ MAPE =============
cat("\nMAPE Neural Network (1) = ",MAPE(ts(kembali_NN1),ts(data_train[2:77]))*100)
## 
## MAPE Neural Network (1) =  16.60254
cat("\nMAPE Neural Network (2) = ",MAPE(ts(kembali_NN2),ts(data_train[2:77]))*100)
## 
## MAPE Neural Network (2) =  16.47924
cat("\nMAPE Neural Network (1,1) = ",MAPE(ts(kembali_NN_1_1),ts(data_train[2:77]))*100)
## 
## MAPE Neural Network (1,1) =  16.6049
cat("\nMAPE Neural Network (1,3) = ",MAPE(ts(kembali_NN_1_3),ts(data_train[2:77]))*100)
## 
## MAPE Neural Network (1,3) =  16.6022
cat("\nMAPE Neural Network (2,1) = ",MAPE(ts(kembali_NN_2_1),ts(data_train[2:77]))*100)#terkecil
## 
## MAPE Neural Network (2,1) =  16.2376
cat("\nMAPE Neural Network (2,3) = ",MAPE(ts(kembali_NN_2_3),ts(data_train[2:77]))*100)
## 
## MAPE Neural Network (2,3) =  16.47754
cat("\n================ RMSE =============")
## 
## ================ RMSE =============
cat("\nRMSE Neural Network (1) = ",RMSE(ts(kembali_NN1),ts(data_train[2:77])))
## 
## RMSE Neural Network (1) =  2015.951
cat("\nRMSE Neural Network (2) = ",RMSE(ts(kembali_NN2),ts(data_train[2:77])))#terkecil
## 
## RMSE Neural Network (2) =  1991.431
cat("\nRMSE Neural Network (1,1) = ",RMSE(ts(kembali_NN_1_1),ts(data_train[2:77])))
## 
## RMSE Neural Network (1,1) =  2016.225
cat("\nRMSE Neural Network (1,3) = ",RMSE(ts(kembali_NN_1_3),ts(data_train[2:77])))
## 
## RMSE Neural Network (1,3) =  2015.913
cat("\nRMSE Neural Network (2,1) = ",RMSE(ts(kembali_NN_2_1),ts(data_train[2:77])))
## 
## RMSE Neural Network (2,1) =  1982.199
cat("\nRMSE Neural Network (2,3) = ",RMSE(ts(kembali_NN_2_3),ts(data_train[2:77])))
## 
## RMSE Neural Network (2,3) =  1934.725
cat("\n================ SMAPE =============")
## 
## ================ SMAPE =============
cat("\nSMAPE Neural Network (1) = ",smape(ts(data_train[2:77]),ts(kembali_NN1))*100)
## 
## SMAPE Neural Network (1) =  16.33906
cat("\nSMAPE Neural Network (2) = ",smape(ts(data_train[2:77]),ts(kembali_NN2))*100)
## 
## SMAPE Neural Network (2) =  16.20107
cat("\nSMAPE Neural Network (1,1) = ",smape(ts(data_train[2:77]),ts(kembali_NN_1_1))*100)
## 
## SMAPE Neural Network (1,1) =  16.34156
cat("\nSMAPE Neural Network (1,3) = ",smape(ts(data_train[2:77]),ts(kembali_NN_1_3))*100)
## 
## SMAPE Neural Network (1,3) =  16.33864
cat("\nSMAPE Neural Network (2,1) = ",smape(ts(data_train[2:77]),ts(kembali_NN_2_1))*100)#terkecil
## 
## SMAPE Neural Network (2,1) =  15.96218
cat("\nSMAPE Neural Network (2,3) = ",smape(ts(data_train[2:77]),ts(kembali_NN_2_3))*100)
## 
## SMAPE Neural Network (2,3) =  16.27737

Berdasarkan hasil yang diperoleh, seluruh arsitektur NN menunjukkan performa peramalan yang baik terhadap data produksi cabai besar di Jawa Timur, ditunjukkan oleh nilai MAPE dan SMAPE yang berada pada kisaran \(10% \leq x < 20%\), serta RMSE dalam rentang \(1000 \leq x < 2000\). Diantara semua arsitektur yang diuji, model NN(2,1) terpilih sebagai model yang terbaik karena menghasilkan nilai MAPE, SMAPE, dan RMSE yang lebih kecil dibandingkan dengan arsitektur NN lainnya.

Hasil peramalan produksi cabai besar di Jawa Timur untuk 12 periode ke depan sebagai berikut.

## Peramalan NN(2,1)
last_lag = tail(trainy$nilai_asli, 1)
forecast_normalized <- numeric(22)
current_lag = last_lag

for (i in 1:22) {
  new_data = data.frame(nilai_lag = current_lag)
  pred = neuralnet::compute(mlp_2_2_1, covariate = new_data)
  forecast_normalized[i] = pred$net.result[1, 1]
  current_lag = forecast_normalized[i]
}

forecast_denormalized = ((forecast_normalized + 1) / 2) * (max_asli - min_asli) + min_asli
forecast_denormalized
##  [1] 8102.874 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899
##  [9] 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899
## [17] 8102.899 8102.899 8102.899 8102.899 8102.899 8102.899

Grafik perbadingan data aktual dan hasil peramalan produksi cabai besar di Jawa Timur menggunakan metode NN (2,1) disajikan sebagai berikut.

## Grafik prediksi dan Peramalan *Neural Newtork*
Hasil_NN = ts(c(NA,kembali_NN_2_1,forecast_denormalized[1:7]))

plot(1:84,cabai_besar, type = "l", main = "Grafik Prediksi dan Peramalan NN", ylab = "Harga Cabai Besar", xlab = "Time") +
  points(cabai_besar, cex = 0.9, pch = 21) +
  lines(1:77,Hasil_NN[1:77], col = "blue") + 
  lines(78:84,Hasil_NN[78:84], col = "red") +
  abline(v = 77, lty = 2) + 
  abline(v = 84, lty = 2)
## integer(0)
legend("topleft", legend = c("Data aktual","Prediksi In sample", "Prediksi Out Sample"), col = c("black","blue","red"), lty = 1, pch = 20, lwd = 3, bty = "o", cex = 0.8)
grid()

Berdasarkan grafik, dapat dilihat bahwa pola data hasil peramalan cenderung mengikuti pola data aktual pada periode in sample. Meskipun demikian, pada periode out sample, hasil peramalan menunjukkan adanya ketidaksesuaian terhadap data aktual. Hal ini menunjukkan bahwa model belum sepenuhnya mampu menangkap pola baru di luar data in sample. Oleh karena itu, diperlukan evaluasi lebih lanjut terhadap performa model agar prediksi yang dihasilkan lebih akurat di masa mendatang.

Berikut adalah perbandingan hasil evaluasi model Neural Network pada data in sample dan data pengujian out sample berdasarkan nilai MAPE, RMSE, dan SMAPE.

## Evaluasi model NN terbaik
cat("================ In Sample =============")
## ================ In Sample =============
cat("\nMAPE Neural In Sample = ",MAPE(ts(kembali_NN_2_1),ts(data_train[2:77]))*100)
## 
## MAPE Neural In Sample =  16.2376
cat("\nRMSE Neural Network In Sample = ",RMSE(ts(kembali_NN_2_1),ts(data_train[2:77])))
## 
## RMSE Neural Network In Sample =  1982.199
cat("\nSMAPE Neural Network In Sample = ",smape(ts(data_train[2:77]),ts(kembali_NN_2_1))*100)
## 
## SMAPE Neural Network In Sample =  15.96218
cat("\n================ Out Sample =============")
## 
## ================ Out Sample =============
cat("\nMAPE Neural Network Out Sample = ",MAPE(forecast_denormalized[1:7],data_test)*100)
## 
## MAPE Neural Network Out Sample =  45.97608
cat("\nRMSE Neural Network Out Sample = ",RMSE(forecast_denormalized[1:7],data_test))
## 
## RMSE Neural Network Out Sample =  2617.406
cat("\nSMAPE Neural Network Out Sample = ",smape(data_test,forecast_denormalized[1:7])*100)
## 
## SMAPE Neural Network Out Sample =  35.25318

Hasil evaluasi model Neural Network menunjukkan bahwa nilai MAPE, RMSE, dan SMAPE pada data in-sample relatif lebih kecil. Evaluasi ini memiliki kesamaan dengan model sebelumnya yakni ARIMA (2,1,0) dengan outlier dimana nilai akurasi pada in-sample lebih kecil dibandingkan dengan out-sample.

KESIMPULAN: Berdasarkan analisis yang telah dilakukan, dapat disimpulkan bahwa dengan melihat nilai akurasi yang diperoleh, model terbaik untuk meramalkan harga cabai besar di Jawa Timur berdasarkan data out-sample yakni model Neural Network dibandingkan model ARIMA (2,1,0) dengan outlier karena memiliki nilai akurasi yang lebih kecil. Namun, jika dibandingkan berdasarkan data in-sample, model ARIMA (2,1,0) memiliki nilai akurasi yang lebih kecil dibandingkan data in-sample model ARIMA (2,1,0) dengan outlier dan data in-sample pada model Neural Network. Hal ini menunjukkan perlunya perbaikan model agar performa prediksi lebih akurat di masa mendatang.