“Prediksi dengan Metode SARIMA dan Neural Network pada Data Produksi Perkebunan Kopi di Indonesia”

(Athaya Azzahra Rahman, Diego Christian, Filzah Syakirah, Widya Ayu Andini)

                               May 12th, 2025

Latar Belakang.

Pengembangan perkebunan merupakan bagian penting dalam sektor pertanian yang perlu ditingkatkan karena mampu mengoptimalkan lahan kosong, meningkatkan pendapatan masyarakat, dan menyumbang devisa negara. Berbagai komoditas seperti cengkeh, teh, kelapa sawit, kakao, dan kopi berperan dalam menyediakan bahan baku industri. Kopi sendiri memiliki potensi dalam menciptakan lapangan kerja, mendorong pertumbuhan wilayah, dan mengembangkan sektor agroindustri (Setiadi, 2023). Kopi merupakan salah satu komoditas yang memiliki peran strategis dalam mendukung perekonomian Indonesia. Komoditas ini berkontribusi dalam mendorong pertumbuhan ekonomi serta meningkatkan kesejahteraan para petani (Pratiwi, Fitriani, dan Berliana, 2023). Menurut Kementrian Pertanian Republik Indonesia produksi kopi di Indonesia sangat signifikan, dengan Indonesia menjadi produsen kopi terbesar keempat di dunia setelah Brazil, Vietnam, dan Kolombia. Pada Tahun 2022-2025, total produksi kopi diperkirakan mencapai 789.000 ton per tahun dengan beberapa daerah sebagai penghasil kopi terbesar di Indonesia termasuk Sumatera Selatan, Lampung, dan Sumatera Utara.

Ketahanan pangan adalah kondisi terpenuhinya kebutuhan pangan secara cukup, aman, bergizi, merata, dan terjangkau, yang selaras dengan nilai agama dan budaya, untuk mendukung hidup sehat dan produktif secara berkelanjutan (Badan Pangan Nasional, 2023). Produksi kopi dipengaruhi oleh berbagai faktor, antara lain tanaman kopi yang telah memasuki usia tidak produktif, penerapan teknik budidaya yang kurang optimal, serta pengaruh perubahan iklim yang cenderung ekstrim pada periode tertentu, sehingga menyebabkan terjadinya fluktuasi (Angka dan Dewi, 2021). Jika dilihat dari fluktuasi data produksi perkebunan kopi tahun 2017 hingga 2024 yang diperoleh dari Badan Pusat Statistik (BPS) yaitu mengalami pola musiman setiap 12 periode. Pola musiman yang dimiliki oleh data produksi perkebunan kopi menandakan bahwa data tersebut lebih cocok menggunakan metode yang cocok dengan pola musiman seperti Autoregressive Integrated Moving Average (ARIMA) Musiman. Selain itu, terdapat metode analisis runtun waktu non linear yang dimana tidak memiliki asumsi khusus pada metodenya. Metode non linear dapat menangkap pola musiman pada data produksi perkebunan kopi. Metode non linear yang banyak digunakan salah satunya adalah Neural Network (NN).

Metode SARIMA

##### Packages and Function Libraries 
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(urca)
## Warning: package 'urca' was built under R version 4.4.3
## 
## Attaching package: 'urca'
## The following objects are masked from 'package:fUnitRoots':
## 
##     punitroot, qunitroot, unitrootTable
library(nortest)
library(neuralnet)
## Warning: package 'neuralnet' was built under R version 4.4.3
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.4.3
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
library(magick)
## Warning: package 'magick' was built under R version 4.4.3
## Linking to ImageMagick 6.9.12.98
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
#Function 
akurasi_model <- function(aktual, prediksi) { 
aktual <- as.vector(aktual) 
prediksi <- as.vector(prediksi) 
residual <- aktual-prediksi
rmse <- sqrt(mean(residual^2, na.rm = TRUE)) 
mape <-mean(abs(residual / aktual), na.rm = TRUE) * 100 
smape <- mean(2 *abs(residual) /(abs(aktual) + abs(prediksi)), na.rm = TRUE) * 100
ukuran_akurasi <- data.frame( MAPE = round(mape, 4), RMSE = round(rmse,4), SMAPE = round(smape, 4) )
return(ukuran_akurasi) 
}

Mengimpor data produksi perkebunan kopi di Indonesia periode Januari 2017 - Desember 2024

data<-read.csv(file.choose(),header=TRUE) 
head(data)
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
data
0.03
0.05
0.05
1.14
2.24
3.89

Grafik Runtun Waktu untuk data keseluruhan

#data Keseluruhan
ts.plot(data$data,type="o",xlab="Periode",ylab="Jumlah Produksi",main="Produksi Perkebunan Kopi Periode Januari 2017 - Desember 2024") 
grid()

Terlihat bahwasanya grafik tersebut menunjukkan pola musiman yang konsisten dengan puncak produksi yang terjadi secara berkala, namun disertai dengan tren menurun dari waktu ke waktu. Pada awal periode, produksi kopi mencapai angka tertinggi, namun secara bertahap mengalami penurunan signifikan hingga akhir periode. Pola ini mengindikasikan adanya siklus panen tahunan yang kuat.

Membagi Data Insample dan Outsample dengan Proporsi 90:10

insample <- data[1:86,]
outsample <- data[87:96,]

Zt <- insample 
outsample <- data.frame(zt = outsample)

Grafik Runtun Waktu Data Insample

ts.plot(Zt,type="o",xlab="Periode",ylab="Jumlah Produksi",main="Produksi Perkebunan Kopi (Insample)") 
grid()

Berdasarkan pada grafik di atas, yang mencerminkan data in-sample sebesar 90% dari keseluruhan data. Terlihat bahwa pola produksi menunjukkan siklus musiman tahunan yang kuat dengan fluktuasi yang konsisten. Pada grafik tersebut juga terlihat bahwasanya data insample cenderung belum stasioner baik dalam variansi maupun rata-rata. Sehingga langkah selanjutnya akan dilakukan pengecekan stasioneritas data insample baik dalam variansi maupun dalam rata-rata.

Pengecekan Stasioneritas dalam Variansi

boxcox(Zt~1)

p<-powerTransform(Zt) 
p 
## Estimated transformation parameter 
##           Zt 
## -0.006150835

Berdasarkan grafik boxcox data insample diperoleh bahwasanya nilai lambda belum mendekati angka 1, dan untuk memvalidasi hal tersebut dilakukan pengecekan nilai lambda. Diperoleh nilai lambda sebesar -0,006150835, angka tersebut belum mendekati angka 1. Sehingga akan dilakukan transformasi pada data insample agar data tersebut stasioner dalam variansi dengan memangkatkan data insample dengan lambda sebesar 0,006150835.

Proses Transformasi

y<-(Zt^p$lambda) 
p1<-powerTransform(y)
p1 
## Estimated transformation parameter 
##         y 
## 0.9999994
boxcox(y~1)

Berdasarkan grafik boxcox data hasil transformasi, diperoleh bahwasanya nilai lambda sudah mendekati angka 1 dengan nilailambda sebesar 0,9999994. Hal ini menunjukkan bahwasanya data sudah stasioner dalam variansi.

Grafik Runtun Waktu Data Setelah Transformasi

ts.plot(y,type="o",xlab="Time",ylab="data", main="TS Plot setelah
transformasi") 
grid()

Berdasarkan grafik diatas, terlihat bahwasanya data setelah transformasi telah stasioner dalam variansi. Namun grafik tersebut terlihat belum stasioner dalam rata-rata. Sehingga untuk memvalidasi hal tersebut, akan dilakukan pengecekan stasioneritas dalam rata-rata.

Pengecekan Stasioner Dalam Rata-rata

acf(y)

Berdasarkan hasil plot ACF di atas terlihat bahwasanya, plot tersebut cenderung dies down. Hal ini menandakan bahwasanya data hasil transformasi belum stasioner dalam rata-rata. Untuk memvalidasi hal tersebut dilakukan pengecekan dengan menggunakan uji ADF (Augmented Dickey-Fuller).

Uji ADF Data Transformasi

adfTest(y)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: 0.0016
##   P VALUE:
##     0.613 
## 
## Description:
##  Mon May 12 14:17:12 2025 by user: filza

Diperoleh nilai p-value sebesar 0,613. Nilai p-value yang lebih besar dari taraf signifikansi sebesar 0,05 menunjukkan bahwasanya data tersebut belum stasioner dalam rata-rata. Sehingga selanjutnya akan dilakukan differencing dengan orde 1 agar data tersebut stasioner dalam rata-rata.

datadiff=diff(y,differences=1)

Grafik Runtun Waktu Data Hasil Differencing Orde 1

ts.plot(datadiff,type="o",xlab="Time",ylab="data", main="TS Plot data In
Sample Setelah Transformasi dan Differencing")

Berdasarkan grafik tersebut, secara visual terlihat bahwasanya data setelah differencing orde 1 telah stasioner baik dalam variansi maupun dalam rata-rata. Untuk memvalidasi hal tersebut dilakukan pengecekan dengan uji ADF.

Uji ADF Setelah Differencing Orde 1

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: -4.7719
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon May 12 14:17:12 2025 by user: filza

Diperoleh nilai p-value sebesar 0,01. Nilai p-value yang lebih kecil dari taraf signifikansi sebesar 0,05 menunjukkan bahwasanya data tersebut telah stasioner dalam rata-rata. Sehingga selanjutnya akan dilakukan estimasi model ARIMA non musiman dari data yang telah di differencing.

Estimasi Model ARIMA Non Musiman

###PACF

Pacf(datadiff, length(datadiff))

Berdasarkan plot PACF di atas, diketahui bahwa grafik PACF mengalami cut off setelah lag 1. Sehingga diperoleh model ARIMA sementara yang terbentuk memiliki orde AR yang mungkin adalah 0, dan 1.

ACF

acf(datadiff,length(datadiff)) #cutoff 1 grid()

Berdasarkan plot ACF di atas, diketahui bahwa grafik ACF mengalami cut off setelah lag 1. Sehingga diperoleh model ARIMA sementara yang terbentuk memiliki orde MA yang mungkin adalah 0, dan 1.Terlihat pada plot tersebut, dimana terdapat pola osilasi setiap 12 kali. Sehingga diasumsikan terjadi pola musiman pada data dengan orde 12. Selanjutnya akan dilakukan pembentukan model ARIMA musiman dengan melakukan differencing orde 12 pada data hasil differencing orde 1.

Estimasi Model ARIMA Musiman

Proses Differencing dengan Orde 12

datadiff2=diff(datadiff,differences=12)
adfTest(datadiff2)
## Warning in adfTest(datadiff2): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -28.6246
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon May 12 14:17:13 2025 by user: filza

Diperoleh nilai p-value sebesar 0,01. Nilai p-value yang lebih kecil dari taraf signifikansi sebesar 0,05 menunjukkan bahwasanya data differencing orde 12 telah stasioner dalam rata-rata. Sehingga selanjutnya akan dilakukan estimasi model ARIMA musiman dari data yang telah di differencing dengan orde 12.

PACF

Pacf(datadiff2, length(datadiff2)) #cutoff 5 grid()

Berdasarkan plot PACF, diketahui bahwa grafik pACF mengalami cut off setelah lag 5. Sehingga diperoleh model SARIMA musiman sementara yang terbentuk memiliki orde AR yang mungkin adalah 0, 1, 2, 3, 4, dan 5.

ACF

Acf(datadiff2,  length(datadiff2)) #cutoff 4 grid()

Berdasarkan plot ACF, diketahui bahwa grafik pACF mengalami cut off setelah lag 4. Sehingga diperoleh model SARIMA musiman sementara yang terbentuk memiliki orde MA yang mungkin adalah 0, 1, 2, 3, dan 4.

Selanjutnya dilakukan uji signifikansi terhadap 87 model yang terbentuk berdasarkan model ARIMA (1,1,1)(5,1,4)^12. Diperoleh 9 model yang seluruh parameternya signifikan sebagai berikut.

#Model SARIMA (1,1,1)(5,1,4)^12
fit1=arima(x=y,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
fit1 
## 
## Call:
## arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##           ma1     sma1
##       -0.8797  -0.3913
## s.e.   0.0886   0.1510
## 
## sigma^2 estimated as 3.879e-05:  log likelihood = 265.33,  aic = -524.66
coeftest(fit1) 
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.879685   0.088581 -9.9309 < 2.2e-16 ***
## sma1 -0.391318   0.151045 -2.5907  0.009577 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit5=arima(x=y,order=c(0,1,1),seasonal=list(order=c(1,1,0),period=12))
fit5 
## 
## Call:
## arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(1, 1, 0), period = 12))
## 
## Coefficients:
##           ma1     sar1
##       -0.8977  -0.2165
## s.e.   0.0889   0.1091
## 
## sigma^2 estimated as 4.104e-05:  log likelihood = 263.93,  aic = -521.87
coeftest(fit5) 
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ma1  -0.89768    0.08888 -10.1000   <2e-16 ***
## sar1 -0.21651    0.10910  -1.9845   0.0472 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit10=arima(x=y,order=c(0,1,1),seasonal=list(order=c(2,1,0),period=12))
fit10 
## 
## Call:
## arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(2, 1, 0), period = 12))
## 
## Coefficients:
##           ma1     sar1     sar2
##       -0.8852  -0.2963  -0.3361
## s.e.   0.0820   0.1113   0.1274
## 
## sigma^2 estimated as 3.634e-05:  log likelihood = 266.95,  aic = -525.91
coeftest(fit10) 
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.885165   0.082005 -10.7940 < 2.2e-16 ***
## sar1 -0.296332   0.111297  -2.6625  0.007756 ** 
## sar2 -0.336089   0.127419  -2.6377  0.008348 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit31=arima(x=y,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))
fit31 
## 
## Call:
## arima(x = y, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.3928  -0.9559  -0.4869
## s.e.  0.1262   0.0599   0.1646
## 
## sigma^2 estimated as 3.319e-05:  log likelihood = 270.07,  aic = -532.14
coeftest(fit31) 
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.392797   0.126205   3.1124  0.001856 ** 
## ma1  -0.955904   0.059921 -15.9527 < 2.2e-16 ***
## sma1 -0.486907   0.164635  -2.9575  0.003102 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit35=arima(x=y,order=c(1,1,1),seasonal=list(order=c(1,1,0),period=12))
fit35 
## 
## Call:
## arima(x = y, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 0), period = 12))
## 
## Coefficients:
##          ar1      ma1     sar1
##       0.3637  -0.9692  -0.2585
## s.e.  0.1221   0.0584   0.1118
## 
## sigma^2 estimated as 3.607e-05:  log likelihood = 268.2,  aic = -528.4
coeftest(fit35) 
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.363672   0.122141   2.9775  0.002906 ** 
## ma1  -0.969201   0.058364 -16.6062 < 2.2e-16 ***
## sar1 -0.258540   0.111786  -2.3128  0.020732 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit61=arima(x=y,order=c(1,1,0),seasonal=list(order=c(0,1,1),period=12))
fit61
## 
## Call:
## arima(x = y, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##           ar1     sma1
##       -0.2369  -0.6057
## s.e.   0.1164   0.1666
## 
## sigma^2 estimated as 4.096e-05:  log likelihood = 262.41,  aic = -518.81
coeftest(fit61)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.23688    0.11637 -2.0356 0.0417875 *  
## sma1 -0.60574    0.16658 -3.6363 0.0002766 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit65=arima(x=y,order=c(1,1,0),seasonal=list(order=c(1,1,0),period=12))
fit65
## 
## Call:
## arima(x = y, order = c(1, 1, 0), seasonal = list(order = c(1, 1, 0), period = 12))
## 
## Coefficients:
##           ar1     sar1
##       -0.2724  -0.3091
## s.e.   0.1138   0.1120
## 
## sigma^2 estimated as 4.7e-05:  log likelihood = 259.51,  aic = -513.02
coeftest(fit65)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1  -0.27237    0.11382 -2.3930 0.016713 * 
## sar1 -0.30913    0.11200 -2.7601 0.005778 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit70=arima(x=y,order=c(1,1,0),seasonal=list(order=c(2,1,0),period=12))
fit70 
## 
## Call:
## arima(x = y, order = c(1, 1, 0), seasonal = list(order = c(2, 1, 0), period = 12))
## 
## Coefficients:
##           ar1     sar1     sar2
##       -0.2571  -0.4321  -0.3353
## s.e.   0.1147   0.1211   0.1329
## 
## sigma^2 estimated as 4.182e-05:  log likelihood = 262.29,  aic = -516.57
coeftest(fit70)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.25710    0.11473 -2.2409 0.0250332 *  
## sar1 -0.43212    0.12107 -3.5690 0.0003583 ***
## sar2 -0.33530    0.13285 -2.5239 0.0116059 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit75=arima(x=y,order=c(1,1,0),seasonal=list(order=c(3,1,0),period=12))
fit75 
## 
## Call:
## arima(x = y, order = c(1, 1, 0), seasonal = list(order = c(3, 1, 0), period = 12))
## 
## Coefficients:
##           ar1     sar1     sar2     sar3
##       -0.2841  -0.5565  -0.4934  -0.3367
## s.e.   0.1147   0.1268   0.1426   0.1528
## 
## sigma^2 estimated as 3.721e-05:  log likelihood = 264.27,  aic = -518.54
coeftest(fit75)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.28414    0.11473 -2.4765 0.0132671 *  
## sar1 -0.55652    0.12675 -4.3906  1.13e-05 ***
## sar2 -0.49338    0.14260 -3.4599 0.0005404 ***
## sar3 -0.33672    0.15279 -2.2038 0.0275367 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan sembilan model di atas, selanjutnya akan dilakukan pengujian diagnostik.

Pengujian Diagnostik

Uji Diagnostik yang pertama dilakukan adalah uji normalitas residual dengan menggunakan Uji Kolmogorov Smirnov.

#Uji Normalitas Residual
res1=resid(fit1) 
res5=resid(fit5) 
res10=resid(fit10) 
res31=resid(fit31)
res35=resid(fit35) 
res61=resid(fit61) 
res65=resid(fit65)
res70=resid(fit70) 
res75=resid(fit75)

ks.test(res1,"pnorm",mean(res1),sd(res1)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res1
## D = 0.12622, p-value = 0.1183
## alternative hypothesis: two-sided
ks.test(res5,"pnorm",mean(res5),sd(res5))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res5
## D = 0.15456, p-value = 0.02928
## alternative hypothesis: two-sided
ks.test(res10,"pnorm",mean(res10),sd(res10)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res10
## D = 0.12443, p-value = 0.128
## alternative hypothesis: two-sided
ks.test(res31,"pnorm",mean(res31),sd(res31)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res31
## D = 0.14204, p-value = 0.0562
## alternative hypothesis: two-sided
ks.test(res35,"pnorm",mean(res35),sd(res35)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res35
## D = 0.14384, p-value = 0.05135
## alternative hypothesis: two-sided
ks.test(res61,"pnorm",mean(res61),sd(res61)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res61
## D = 0.12982, p-value = 0.1007
## alternative hypothesis: two-sided
ks.test(res65,"pnorm",mean(res65),sd(res65)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res65
## D = 0.13432, p-value = 0.08173
## alternative hypothesis: two-sided
ks.test(res70,"pnorm",mean(res70),sd(res70)) #Normal
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res70
## D = 0.11492, p-value = 0.1907
## alternative hypothesis: two-sided
ks.test(res75,"pnorm",mean(res75),sd(res75))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res75
## D = 0.14576, p-value = 0.04656
## alternative hypothesis: two-sided

Berdasarkan hasil pengujian normalitas residual, diperoleh 7 model dengan nilai p-value lebih dari taraf signifikansi sebesar 0,05. Sehingga disimpulkan bahwa terdapat 7 model dengan residual yang mengikuti distribusi normal. Sehingga dari 7 model tersebut akan dilakukan pengujian diagnostik selanjutnya yaitu pengujian asumsi white noise.

Uji Asumsi White Noise

#Uji White Noise
tsdiag(fit1, length(y)) 

tsdiag(fit10, length(y)) 

tsdiag(fit31,length(y)) 

tsdiag(fit35, length(y)) 

tsdiag(fit61, length(y))

tsdiag(fit65, length(y)) 

tsdiag(fit70, length(y))

Berdasarkan hasil grafik Ljung Box ketujuh model, diperoleh bahwasanya hanya fit31 yang memenuhi syarat white noise karena nilai p-value setiap lag lebih besar dari 0,05 (secara visual berada di atas garis batas signifikansi 0,05).

Sehingga disimpulkan bahwa fit31 yaitu model ARIMA(1,1,1)(0,1,1)^12 merupakan model terbaik, karena model tersebut memenuhi seluruh asumsi. Maka model tersebut yang akan terpilih menjadi model untuk peramalan. Sebelum peramalan akan dilakukan pengecekan nilai evaluasi pada model terbaik.

Pengecekan Nilai Evaluasi Model

model1=arima(y,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))
summary(model1)
## 
## Call:
## arima(x = y, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.3928  -0.9559  -0.4869
## s.e.  0.1262   0.0599   0.1646
## 
## sigma^2 estimated as 3.319e-05:  log likelihood = 270.07,  aic = -532.14
## 
## Training set error measures:
##                         ME        RMSE         MAE         MPE      MAPE
## Training set -0.0002203287 0.005321582 0.003322583 -0.02267762 0.3272724
##                  MASE       ACF1
## Training set 0.608996 0.02170787

Berdasarkan hasil evaluasi di atas diperoleh nilai MAPE sebesar 32%, hal ini menunjukkan bahwasanya model memiliki akurasi yang baik dalam memprediksi data. Dengan model tersebut, selanjutnya akan dilakukan peramalan.

Prediksi

#Prediksi 
prediksi_sarima <- fitted(model1)^(1/p$lambda)

Peramalan

#Peramalan 
peramalan_sarima <- predict(model1, n.ahead = nrow(outsample))$pred^(1/-0.006150835)

Akurasi Peramalan

#Akurasi Peramalan 
akurasi_sarima <- akurasi_model(outsample$zt,peramalan_sarima)
akurasi_sarima
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
MAPE RMSE SMAPE
31.047 0.2049 37.1082

#Perbandingan Akurasi

#Perbandingan Akurasi
data.frame( Model = c( "SARIMA"), MAPE = c( akurasi_sarima$MAPE),
  RMSE = c( akurasi_sarima$RMSE), SMAPE = c( akurasi_sarima$SMAPE) )
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
Model MAPE RMSE SMAPE
SARIMA 31.047 0.2049 37.1082

Grafik Perbandingan Metode SARIMA

#Grafik Perbandingan 
ts.plot(data$data,type = "o",col="grey22",main="Produksi Kopi 2017 sampai 2024",xlab="Periode",ylab="Jumlah Pencarian",xlim=c(0,nrow(data)),ylim=c(0,10),lwd=1.5) 
lines(c(prediksi_sarima, peramalan_sarima), type = "o",col="green", lwd=1.5) 
abline(v=length(Zt)+0.5, col="orange", lwd=1.5,lty=2)
legend('topleft',legend=c("Aktual","Prediksi"),col=c("grey22", "green"),lty=1, cex = 1, lwd = 2)

Berdasarkan grafik tersebut terlihat bahwa garis hitam merepresentasikan data aktual, sementara garis hijau menunjukkan hasil prediksi model. Secara umum, model mampu menangkap pola musiman tahunan dengan cukup baik, ditunjukkan oleh kemiripan bentuk dan letak puncak antara data aktual dan prediksi. Puncak pencarian yang terjadi secara berkala berhasil direplikasi oleh model, meskipun dengan intensitas yang sedikit lebih rendah pada beberapa titik. Hal ini mengindikasikan bahwa model dapat mengidentifikasi karakteristik musiman dalam data. Setelah titik vertikal oranye (batas data insample), hasil peramalan tetap mempertahankan pola musiman yang serupa, hal ini menunjukkan kemampuan model dalam memproyeksikan tren jangka pendek dengan baik.

NEURAL NETWORKK

Zt<-ts(data$data)
Zt
## Time Series:
## Start = 1 
## End = 96 
## Frequency = 1 
##  [1] 0.03 0.05 0.05 1.14 2.24 3.89 4.75 5.14 4.89 4.48 2.79 0.85 0.02 0.07 0.03
## [16] 0.71 2.33 4.40 7.05 5.94 3.59 2.17 1.25 0.57 0.57 0.01 0.01 0.02 0.19 0.86
## [31] 2.96 2.74 1.78 1.13 0.22 0.13 0.05 0.04 0.05 0.03 0.04 0.33 1.11 1.28 1.47
## [46] 0.54 0.09 0.06 0.01 0.01 0.02 0.07 0.22 0.87 1.60 1.71 0.53 0.12 0.12 0.05
## [61] 0.01 0.01 0.03 0.03 0.32 0.91 1.07 1.04 0.37 0.14 0.06 0.02 0.05 0.05 0.01
## [76] 0.04 0.40 0.85 0.89 0.70 0.22 0.05 0.04 0.02 0.03 0.03 0.01 0.03 0.37 0.90
## [91] 0.99 0.88 0.29 0.09 0.05 0.02
jumlah_train = round(NROW(Zt)*0.90,0)
jumlah_test = NROW(Zt)-jumlah_train
data_train = Zt[1:jumlah_train]
data_test  = Zt[(jumlah_train+1):(jumlah_test+jumlah_train)]

STANDARISASI DATA

max<-max(data_train)
min<-min(data_train)
newmax<-(1)
newmin<-(-1)
Yt<-(((data_train-min)/(max-min))*(newmax-newmin))+newmin
Yt
##  [1] -0.99431818 -0.98863636 -0.98863636 -0.67897727 -0.36647727  0.10227273
##  [7]  0.34659091  0.45738636  0.38636364  0.26988636 -0.21022727 -0.76136364
## [13] -0.99715909 -0.98295455 -0.99431818 -0.80113636 -0.34090909  0.24715909
## [19]  1.00000000  0.68465909  0.01704545 -0.38636364 -0.64772727 -0.84090909
## [25] -0.84090909 -1.00000000 -1.00000000 -0.99715909 -0.94886364 -0.75852273
## [31] -0.16193182 -0.22443182 -0.49715909 -0.68181818 -0.94034091 -0.96590909
## [37] -0.98863636 -0.99147727 -0.98863636 -0.99431818 -0.99147727 -0.90909091
## [43] -0.68750000 -0.63920455 -0.58522727 -0.84943182 -0.97727273 -0.98579545
## [49] -1.00000000 -1.00000000 -0.99715909 -0.98295455 -0.94034091 -0.75568182
## [55] -0.54829545 -0.51704545 -0.85227273 -0.96875000 -0.96875000 -0.98863636
## [61] -1.00000000 -1.00000000 -0.99431818 -0.99431818 -0.91193182 -0.74431818
## [67] -0.69886364 -0.70738636 -0.89772727 -0.96306818 -0.98579545 -0.99715909
## [73] -0.98863636 -0.98863636 -1.00000000 -0.99147727 -0.88920455 -0.76136364
## [79] -0.75000000 -0.80397727 -0.94034091 -0.98863636 -0.99147727 -0.99715909
## [85] -0.99431818 -0.99431818

PREPROCESSING STANDARDIZED

pacf(data_train)

Berdasarkan grafik dapat dilihat bahwa terdapat 2 lag dengan nilai PACF yang signifikan, yaitu lag ke-1 dan lag ke-2. Berdasarkan hal tersebut, variabel masukan yang digunakan pada metode NN adalah error t???1 dan error t???2 sebagai x1 dan x2. Selanjutnya, keluaran yang digunakan merupakan data produksi perkebunan kopi di Indonesia yang telah distandarisasi dan dilambangkan dengan Y, selanjutnya data masukan dan keluaran dimulai dari t = 3.

MENENTUKAN INPUT DAN OUTPUT

X<-Yt[1:84]
Y<-Yt[3:86]
length(X)
## [1] 84
length(Y)
## [1] 84

Kemudian, dilakukan proses pelatihan backpropagation dengan 1 neuron, 2 neuron, dan 3 neuron pada lapisan tersembunyi. Kondisi pemberhentian yang digunakan adalah maksimum iterasi sebesar 10.000 dan target error sebesar 0,05. Berdasarkan proses pelatihan yang dilakukan, diperoleh arsitektur NN sebagai berikut:

METODE NN DENGAN 1 NEURON PADA HIDDEN LAYER

data<-cbind(Y,X)
set.seed(1000)
nn1=neuralnet(Y~X,data=data,hidden=1,act.fct="tanh",threshold=0.05,
              algorithm="backprop",stepmax=10000,startweights=NULL,
              learningrate=0.01,linear.output=T)
nn1
## $call
## neuralnet(formula = Y ~ X, data = data, hidden = 1, threshold = 0.05, 
##     stepmax = 10000, startweights = NULL, learningrate = 0.01, 
##     algorithm = "backprop", act.fct = "tanh", linear.output = T)
## 
## $response
##              Y
## 1  -0.98863636
## 2  -0.67897727
## 3  -0.36647727
## 4   0.10227273
## 5   0.34659091
## 6   0.45738636
## 7   0.38636364
## 8   0.26988636
## 9  -0.21022727
## 10 -0.76136364
## 11 -0.99715909
## 12 -0.98295455
## 13 -0.99431818
## 14 -0.80113636
## 15 -0.34090909
## 16  0.24715909
## 17  1.00000000
## 18  0.68465909
## 19  0.01704545
## 20 -0.38636364
## 21 -0.64772727
## 22 -0.84090909
## 23 -0.84090909
## 24 -1.00000000
## 25 -1.00000000
## 26 -0.99715909
## 27 -0.94886364
## 28 -0.75852273
## 29 -0.16193182
## 30 -0.22443182
## 31 -0.49715909
## 32 -0.68181818
## 33 -0.94034091
## 34 -0.96590909
## 35 -0.98863636
## 36 -0.99147727
## 37 -0.98863636
## 38 -0.99431818
## 39 -0.99147727
## 40 -0.90909091
## 41 -0.68750000
## 42 -0.63920455
## 43 -0.58522727
## 44 -0.84943182
## 45 -0.97727273
## 46 -0.98579545
## 47 -1.00000000
## 48 -1.00000000
## 49 -0.99715909
## 50 -0.98295455
## 51 -0.94034091
## 52 -0.75568182
## 53 -0.54829545
## 54 -0.51704545
## 55 -0.85227273
## 56 -0.96875000
## 57 -0.96875000
## 58 -0.98863636
## 59 -1.00000000
## 60 -1.00000000
## 61 -0.99431818
## 62 -0.99431818
## 63 -0.91193182
## 64 -0.74431818
## 65 -0.69886364
## 66 -0.70738636
## 67 -0.89772727
## 68 -0.96306818
## 69 -0.98579545
## 70 -0.99715909
## 71 -0.98863636
## 72 -0.98863636
## 73 -1.00000000
## 74 -0.99147727
## 75 -0.88920455
## 76 -0.76136364
## 77 -0.75000000
## 78 -0.80397727
## 79 -0.94034091
## 80 -0.98863636
## 81 -0.99147727
## 82 -0.99715909
## 83 -0.99431818
## 84 -0.99431818
## 
## $covariate
##                  
##  [1,] -0.99431818
##  [2,] -0.98863636
##  [3,] -0.98863636
##  [4,] -0.67897727
##  [5,] -0.36647727
##  [6,]  0.10227273
##  [7,]  0.34659091
##  [8,]  0.45738636
##  [9,]  0.38636364
## [10,]  0.26988636
## [11,] -0.21022727
## [12,] -0.76136364
## [13,] -0.99715909
## [14,] -0.98295455
## [15,] -0.99431818
## [16,] -0.80113636
## [17,] -0.34090909
## [18,]  0.24715909
## [19,]  1.00000000
## [20,]  0.68465909
## [21,]  0.01704545
## [22,] -0.38636364
## [23,] -0.64772727
## [24,] -0.84090909
## [25,] -0.84090909
## [26,] -1.00000000
## [27,] -1.00000000
## [28,] -0.99715909
## [29,] -0.94886364
## [30,] -0.75852273
## [31,] -0.16193182
## [32,] -0.22443182
## [33,] -0.49715909
## [34,] -0.68181818
## [35,] -0.94034091
## [36,] -0.96590909
## [37,] -0.98863636
## [38,] -0.99147727
## [39,] -0.98863636
## [40,] -0.99431818
## [41,] -0.99147727
## [42,] -0.90909091
## [43,] -0.68750000
## [44,] -0.63920455
## [45,] -0.58522727
## [46,] -0.84943182
## [47,] -0.97727273
## [48,] -0.98579545
## [49,] -1.00000000
## [50,] -1.00000000
## [51,] -0.99715909
## [52,] -0.98295455
## [53,] -0.94034091
## [54,] -0.75568182
## [55,] -0.54829545
## [56,] -0.51704545
## [57,] -0.85227273
## [58,] -0.96875000
## [59,] -0.96875000
## [60,] -0.98863636
## [61,] -1.00000000
## [62,] -1.00000000
## [63,] -0.99431818
## [64,] -0.99431818
## [65,] -0.91193182
## [66,] -0.74431818
## [67,] -0.69886364
## [68,] -0.70738636
## [69,] -0.89772727
## [70,] -0.96306818
## [71,] -0.98579545
## [72,] -0.99715909
## [73,] -0.98863636
## [74,] -0.98863636
## [75,] -1.00000000
## [76,] -0.99147727
## [77,] -0.88920455
## [78,] -0.76136364
## [79,] -0.75000000
## [80,] -0.80397727
## [81,] -0.94034091
## [82,] -0.98863636
## [83,] -0.99147727
## [84,] -0.99715909
## 
## $model.list
## $model.list$response
## [1] "Y"
## 
## $model.list$variables
## [1] "X"
## 
## 
## $err.fct
## function (x, y) 
## {
##     1/2 * (y - x)^2
## }
## <bytecode: 0x000001d0ad077470>
## <environment: 0x000001d0ad0703f8>
## attr(,"type")
## [1] "sse"
## 
## $act.fct
## function (x) 
## {
##     tanh(x)
## }
## <bytecode: 0x000001d0ad0754a8>
## <environment: 0x000001d0ad076628>
## attr(,"type")
## [1] "tanh"
## 
## $linear.output
## [1] TRUE
## 
## $data
##              Y           X
## 1  -0.98863636 -0.99431818
## 2  -0.67897727 -0.98863636
## 3  -0.36647727 -0.98863636
## 4   0.10227273 -0.67897727
## 5   0.34659091 -0.36647727
## 6   0.45738636  0.10227273
## 7   0.38636364  0.34659091
## 8   0.26988636  0.45738636
## 9  -0.21022727  0.38636364
## 10 -0.76136364  0.26988636
## 11 -0.99715909 -0.21022727
## 12 -0.98295455 -0.76136364
## 13 -0.99431818 -0.99715909
## 14 -0.80113636 -0.98295455
## 15 -0.34090909 -0.99431818
## 16  0.24715909 -0.80113636
## 17  1.00000000 -0.34090909
## 18  0.68465909  0.24715909
## 19  0.01704545  1.00000000
## 20 -0.38636364  0.68465909
## 21 -0.64772727  0.01704545
## 22 -0.84090909 -0.38636364
## 23 -0.84090909 -0.64772727
## 24 -1.00000000 -0.84090909
## 25 -1.00000000 -0.84090909
## 26 -0.99715909 -1.00000000
## 27 -0.94886364 -1.00000000
## 28 -0.75852273 -0.99715909
## 29 -0.16193182 -0.94886364
## 30 -0.22443182 -0.75852273
## 31 -0.49715909 -0.16193182
## 32 -0.68181818 -0.22443182
## 33 -0.94034091 -0.49715909
## 34 -0.96590909 -0.68181818
## 35 -0.98863636 -0.94034091
## 36 -0.99147727 -0.96590909
## 37 -0.98863636 -0.98863636
## 38 -0.99431818 -0.99147727
## 39 -0.99147727 -0.98863636
## 40 -0.90909091 -0.99431818
## 41 -0.68750000 -0.99147727
## 42 -0.63920455 -0.90909091
## 43 -0.58522727 -0.68750000
## 44 -0.84943182 -0.63920455
## 45 -0.97727273 -0.58522727
## 46 -0.98579545 -0.84943182
## 47 -1.00000000 -0.97727273
## 48 -1.00000000 -0.98579545
## 49 -0.99715909 -1.00000000
## 50 -0.98295455 -1.00000000
## 51 -0.94034091 -0.99715909
## 52 -0.75568182 -0.98295455
## 53 -0.54829545 -0.94034091
## 54 -0.51704545 -0.75568182
## 55 -0.85227273 -0.54829545
## 56 -0.96875000 -0.51704545
## 57 -0.96875000 -0.85227273
## 58 -0.98863636 -0.96875000
## 59 -1.00000000 -0.96875000
## 60 -1.00000000 -0.98863636
## 61 -0.99431818 -1.00000000
## 62 -0.99431818 -1.00000000
## 63 -0.91193182 -0.99431818
## 64 -0.74431818 -0.99431818
## 65 -0.69886364 -0.91193182
## 66 -0.70738636 -0.74431818
## 67 -0.89772727 -0.69886364
## 68 -0.96306818 -0.70738636
## 69 -0.98579545 -0.89772727
## 70 -0.99715909 -0.96306818
## 71 -0.98863636 -0.98579545
## 72 -0.98863636 -0.99715909
## 73 -1.00000000 -0.98863636
## 74 -0.99147727 -0.98863636
## 75 -0.88920455 -1.00000000
## 76 -0.76136364 -0.99147727
## 77 -0.75000000 -0.88920455
## 78 -0.80397727 -0.76136364
## 79 -0.94034091 -0.75000000
## 80 -0.98863636 -0.80397727
## 81 -0.99147727 -0.94034091
## 82 -0.99715909 -0.98863636
## 83 -0.99431818 -0.99147727
## 84 -0.99431818 -0.99715909
## 
## $exclude
## NULL
## 
## $net.result
## $net.result[[1]]
##              [,1]
##  [1,] -0.89369584
##  [2,] -0.89027745
##  [3,] -0.89027745
##  [4,] -0.67761705
##  [5,] -0.44024171
##  [6,] -0.15123408
##  [7,] -0.05533150
##  [8,] -0.02314425
##  [9,] -0.04304309
## [10,] -0.08152526
## [11,] -0.32997543
## [12,] -0.73842721
## [13,] -0.89539738
## [14,] -0.88683872
## [15,] -0.89369584
## [16,] -0.76688900
## [17,] -0.42149301
## [18,] -0.08994782
## [19,]  0.06487893
## [20,]  0.02537709
## [21,] -0.19368965
## [22,] -0.45497465
## [23,] -0.65404718
## [24,] -0.79464112
## [25,] -0.79464112
## [26,] -0.89709379
## [27,] -0.89709379
## [28,] -0.89539738
## [29,] -0.86578517
## [30,] -0.73636942
## [31,] -0.29835846
## [32,] -0.33952068
## [33,] -0.53881370
## [34,] -0.67974826
## [35,] -0.86041089
## [36,] -0.87640150
## [37,] -0.89027745
## [38,] -0.89198920
## [39,] -0.89027745
## [40,] -0.89369584
## [41,] -0.89198920
## [42,] -0.84033701
## [43,] -0.68400435
## [44,] -0.64758314
## [45,] -0.60638557
## [46,] -0.80048769
## [47,] -0.88337974
## [48,] -0.88856062
## [49,] -0.89709379
## [50,] -0.89709379
## [51,] -0.89539738
## [52,] -0.88683872
## [53,] -0.86041089
## [54,] -0.73430851
## [55,] -0.57804689
## [56,] -0.55405658
## [57,] -0.80242834
## [58,] -0.87815356
## [59,] -0.87815356
## [60,] -0.89027745
## [61,] -0.89709379
## [62,] -0.89709379
## [63,] -0.89369584
## [64,] -0.89369584
## [65,] -0.84218541
## [66,] -0.72603431
## [67,] -0.69249009
## [68,] -0.69883000
## [69,] -0.83289743
## [70,] -0.87464446
## [71,] -0.88856062
## [72,] -0.89539738
## [73,] -0.89027745
## [74,] -0.89027745
## [75,] -0.89709379
## [76,] -0.89198920
## [77,] -0.82727006
## [78,] -0.73842721
## [79,] -0.73017744
## [80,] -0.76889585
## [81,] -0.86041089
## [82,] -0.89027745
## [83,] -0.89198920
## [84,] -0.89539738
## 
## 
## $weights
## $weights[[1]]
## $weights[[1]][[1]]
##            [,1]
## [1,] -0.6145484
## [2,] -1.1282497
## 
## $weights[[1]][[2]]
##            [,1]
## [1,] -0.5752797
## [2,] -0.6806213
## 
## 
## 
## $generalized.weights
## $generalized.weights[[1]]
##              [,1]
##  [1,]  -0.3544346
##  [2,]  -0.3585718
##  [3,]  -0.3585718
##  [4,]  -0.6602410
##  [5,]  -1.1634374
##  [6,]  -2.6985676
##  [7,]  -5.4760699
##  [8,] -11.0879845
##  [9,]  -6.6449726
## [10,]  -4.1258205
## [11,]  -1.5224963
## [12,]  -0.5638289
## [13,]  -0.3523823
## [14,]  -0.3627526
## [15,]  -0.3544346
## [16,]  -0.5218056
## [17,]  -1.2162346
## [18,]  -3.8500226
## [19,]   1.4602051
## [20,]   6.8669232
## [21,]  -2.2773482
## [22,]  -1.1237841
## [23,]  -0.7003223
## [24,]  -0.4825373
## [25,]  -0.4825373
## [26,]  -0.3503407
## [27,]  -0.3503407
## [28,]  -0.3523823
## [29,]  -0.3887742
## [30,]  -0.5669396
## [31,]  -1.6541834
## [32,]  -1.4858857
## [33,]  -0.9235016
## [34,]  -0.6566960
## [35,]  -0.3955364
## [36,]  -0.3755605
## [37,]  -0.3585718
## [38,]  -0.3564978
## [39,]  -0.3585718
## [40,]  -0.3544346
## [41,]  -0.3564978
## [42,]  -0.4212407
## [43,]  -0.6496547
## [44,]  -0.7116050
## [45,]  -0.7866908
## [46,]  -0.4744681
## [47,]  -0.3669774
## [48,]  -0.3606567
## [49,]  -0.3503407
## [50,]  -0.3503407
## [51,]  -0.3523823
## [52,]  -0.3627526
## [53,]  -0.3955364
## [54,]  -0.5700651
## [55,]  -0.8418227
## [56,]  -0.8909788
## [57,]  -0.4718048
## [58,]  -0.3733979
## [59,]  -0.3733979
## [60,]  -0.3585718
## [61,]  -0.3503407
## [62,]  -0.3503407
## [63,]  -0.3544346
## [64,]  -0.3544346
## [65,]  -0.4188438
## [66,]  -0.5827163
## [67,]  -0.6357649
## [68,]  -0.6255147
## [69,]  -0.4309507
## [70,]  -0.3777343
## [71,]  -0.3606567
## [72,]  -0.3523823
## [73,]  -0.3585718
## [74,]  -0.3585718
## [75,]  -0.3503407
## [76,]  -0.3564978
## [77,]  -0.4383631
## [78,]  -0.5638289
## [79,]  -0.5763608
## [80,]  -0.5189110
## [81,]  -0.3955364
## [82,]  -0.3585718
## [83,]  -0.3564978
## [84,]  -0.3523823
## 
## 
## $startweights
## $startweights[[1]]
## $startweights[[1]][[1]]
##            [,1]
## [1,] -0.4457783
## [2,] -1.2058566
## 
## $startweights[[1]][[2]]
##            [,1]
## [1,] 0.04112631
## [2,] 0.63938841
## 
## 
## 
## $result.matrix
##                               [,1]
## error                   5.15848150
## reached.threshold       0.04981731
## steps                 161.00000000
## Intercept.to.1layhid1  -0.61454842
## X.to.1layhid1          -1.12824969
## Intercept.to.Y         -0.57527967
## 1layhid1.to.Y          -0.68062134
## 
## attr(,"class")
## [1] "nn"
plot(nn1,cex=1)

img1 <- image_read("D:/kuliah/6/ARWW II/RPUBS/NN1.png")

print(img1)
## # A tibble: 1 × 7
##   format width height colorspace matte filesize density
##   <chr>  <int>  <int> <chr>      <lgl>    <int> <chr>  
## 1 PNG      645    645 sRGB       FALSE     5702 72x72

Predict1<-ts(nn1$net.result[[1]])
Predict1
## Time Series:
## Start = 1 
## End = 84 
## Frequency = 1 
##          Series 1
##  [1,] -0.89369584
##  [2,] -0.89027745
##  [3,] -0.89027745
##  [4,] -0.67761705
##  [5,] -0.44024171
##  [6,] -0.15123408
##  [7,] -0.05533150
##  [8,] -0.02314425
##  [9,] -0.04304309
## [10,] -0.08152526
## [11,] -0.32997543
## [12,] -0.73842721
## [13,] -0.89539738
## [14,] -0.88683872
## [15,] -0.89369584
## [16,] -0.76688900
## [17,] -0.42149301
## [18,] -0.08994782
## [19,]  0.06487893
## [20,]  0.02537709
## [21,] -0.19368965
## [22,] -0.45497465
## [23,] -0.65404718
## [24,] -0.79464112
## [25,] -0.79464112
## [26,] -0.89709379
## [27,] -0.89709379
## [28,] -0.89539738
## [29,] -0.86578517
## [30,] -0.73636942
## [31,] -0.29835846
## [32,] -0.33952068
## [33,] -0.53881370
## [34,] -0.67974826
## [35,] -0.86041089
## [36,] -0.87640150
## [37,] -0.89027745
## [38,] -0.89198920
## [39,] -0.89027745
## [40,] -0.89369584
## [41,] -0.89198920
## [42,] -0.84033701
## [43,] -0.68400435
## [44,] -0.64758314
## [45,] -0.60638557
## [46,] -0.80048769
## [47,] -0.88337974
## [48,] -0.88856062
## [49,] -0.89709379
## [50,] -0.89709379
## [51,] -0.89539738
## [52,] -0.88683872
## [53,] -0.86041089
## [54,] -0.73430851
## [55,] -0.57804689
## [56,] -0.55405658
## [57,] -0.80242834
## [58,] -0.87815356
## [59,] -0.87815356
## [60,] -0.89027745
## [61,] -0.89709379
## [62,] -0.89709379
## [63,] -0.89369584
## [64,] -0.89369584
## [65,] -0.84218541
## [66,] -0.72603431
## [67,] -0.69249009
## [68,] -0.69883000
## [69,] -0.83289743
## [70,] -0.87464446
## [71,] -0.88856062
## [72,] -0.89539738
## [73,] -0.89027745
## [74,] -0.89027745
## [75,] -0.89709379
## [76,] -0.89198920
## [77,] -0.82727006
## [78,] -0.73842721
## [79,] -0.73017744
## [80,] -0.76889585
## [81,] -0.86041089
## [82,] -0.89027745
## [83,] -0.89198920
## [84,] -0.89539738

METODE NN DENGAN 2 NEURON PADA HIDDEN LAYER

Data<-cbind(Y,X)
set.seed(1000)
nn2=neuralnet(Y~X,data=data,hidden=c(1,1),act.fct="tanh",threshold=0.05,
              algorithm="backprop",stepmax=10000,startweights=NULL,
              learningrate=0.01,linear.output=T)
plot(nn2,cex=1)

img2 <- image_read("D:/kuliah/6/ARWW II/RPUBS/NN2.png")

print(img2)
## # A tibble: 1 × 7
##   format width height colorspace matte filesize density
##   <chr>  <int>  <int> <chr>      <lgl>    <int> <chr>  
## 1 PNG      645    645 sRGB       FALSE     7890 72x72

Predict2=ts(nn2$net.result[[1]])
Predict2
## Time Series:
## Start = 1 
## End = 84 
## Frequency = 1 
##          Series 1
##  [1,] -0.88022593
##  [2,] -0.87799698
##  [3,] -0.87799698
##  [4,] -0.70711776
##  [5,] -0.44901223
##  [6,] -0.12583758
##  [7,] -0.04191241
##  [8,] -0.01768547
##  [9,] -0.03241753
## [10,] -0.06313913
## [11,] -0.31891183
## [12,] -0.76254593
## [13,] -0.88132957
## [14,] -0.87573887
## [15,] -0.88022593
## [16,] -0.78668344
## [17,] -0.42687524
## [18,] -0.07024428
## [19,]  0.03727533
## [20,]  0.01477652
## [21,] -0.16805325
## [22,] -0.46633614
## [23,] -0.68425876
## [24,] -0.80909943
## [25,] -0.80909943
## [26,] -0.88242601
## [27,] -0.88242601
## [28,] -0.88132957
## [29,] -0.86156351
## [30,] -0.76075578
## [31,] -0.28235354
## [32,] -0.33006489
## [33,] -0.56264186
## [34,] -0.70914759
## [35,] -0.85784773
## [36,] -0.86878704
## [37,] -0.87799698
## [38,] -0.87911508
## [39,] -0.87799698
## [40,] -0.88022593
## [41,] -0.87911508
## [42,] -0.84361452
## [43,] -0.71318259
## [44,] -0.67786021
## [45,] -0.63582591
## [46,] -0.81368109
## [47,] -0.87345136
## [48,] -0.87687158
## [49,] -0.88242601
## [50,] -0.88242601
## [51,] -0.88132957
## [52,] -0.87573887
## [53,] -0.85784773
## [54,] -0.75895685
## [55,] -0.60572730
## [56,] -0.57956302
## [57,] -0.81519108
## [58,] -0.86996436
## [59,] -0.86996436
## [60,] -0.87799698
## [61,] -0.88242601
## [62,] -0.88242601
## [63,] -0.88022593
## [64,] -0.88022593
## [65,] -0.84494859
## [66,] -0.75167346
## [67,] -0.72115293
## [68,] -0.72704252
## [69,] -0.83819656
## [70,] -0.86760218
## [71,] -0.87687158
## [72,] -0.88132957
## [73,] -0.87799698
## [74,] -0.87799698
## [75,] -0.88242601
## [76,] -0.87911508
## [77,] -0.83404664
## [78,] -0.76254593
## [79,] -0.75533267
## [80,] -0.78834150
## [81,] -0.85784773
## [82,] -0.87799698
## [83,] -0.87911508
## [84,] -0.88132957

METODE NN DENGAN 3 NEURON PADA HIDDEN LAYER

Data<-cbind(Y,X)
set.seed(1000)
nn3=neuralnet(Y~X,data=data,hidden=3,act.fct="tanh",threshold=0.05,
              algorithm="backprop",stepmax=10000,startweights=NULL,
              learningrate=0.01,linear.output=T)
plot(nn3,cex=1)

img3 <- image_read("D:/kuliah/6/ARWW II/RPUBS/NN3.png")

print(img3)
## # A tibble: 1 × 7
##   format width height colorspace matte filesize density
##   <chr>  <int>  <int> <chr>      <lgl>    <int> <chr>  
## 1 PNG      645    645 sRGB       FALSE     8954 72x72

Predict3<-ts(nn3$net.result[[1]])
Predict3
## Time Series:
## Start = 1 
## End = 84 
## Frequency = 1 
##           Series 1
##  [1,] -0.882922254
##  [2,] -0.880106080
##  [3,] -0.880106080
##  [4,] -0.696937965
##  [5,] -0.471762651
##  [6,] -0.155967348
##  [7,] -0.037216200
##  [8,]  0.004424309
##  [9,] -0.021420705
## [10,] -0.070448046
## [11,] -0.357716368
## [12,] -0.750938090
## [13,] -0.884322501
## [14,] -0.877269016
## [15,] -0.882922254
## [16,] -0.775756024
## [17,] -0.452851482
## [18,] -0.081005597
## [19,]  0.120871892
## [20,]  0.068589663
## [21,] -0.206048932
## [22,] -0.486494357
## [23,] -0.675637980
## [24,] -0.799682585
## [25,] -0.799682585
## [26,] -0.885717519
## [27,] -0.885717519
## [28,] -0.884322501
## [29,] -0.859809496
## [30,] -0.749132644
## [31,] -0.323670540
## [32,] -0.367870406
## [33,] -0.568294085
## [34,] -0.698853547
## [35,] -0.855328090
## [36,] -0.868632739
## [37,] -0.880106080
## [38,] -0.881516779
## [39,] -0.880106080
## [40,] -0.882922254
## [41,] -0.881516779
## [42,] -0.838501486
## [43,] -0.702673926
## [44,] -0.669759229
## [45,] -0.631904717
## [46,] -0.804689258
## [47,] -0.874411087
## [48,] -0.878690158
## [49,] -0.885717519
## [50,] -0.885717519
## [51,] -0.884322501
## [52,] -0.877269016
## [53,] -0.855328090
## [54,] -0.747322945
## [55,] -0.605464660
## [56,] -0.582815545
## [57,] -0.806348522
## [58,] -0.870085128
## [59,] -0.870085128
## [60,] -0.880106080
## [61,] -0.885717519
## [62,] -0.885717519
## [63,] -0.882922254
## [64,] -0.882922254
## [65,] -0.840056653
## [66,] -0.740042042
## [67,] -0.710270772
## [68,] -0.715929139
## [69,] -0.832230292
## [70,] -0.867175157
## [71,] -0.878690158
## [72,] -0.884322501
## [73,] -0.880106080
## [74,] -0.880106080
## [75,] -0.885717519
## [76,] -0.881516779
## [77,] -0.827474106
## [78,] -0.750938090
## [79,] -0.743690870
## [80,] -0.777495209
## [81,] -0.855328090
## [82,] -0.880106080
## [83,] -0.881516779
## [84,] -0.884322501

DESTANDARISASI PREDIKSI

Prediksi1<-(((Predict1-newmin)/(newmax-newmin))*(max-min))+min
Prediksi2<-(((Predict2-newmin)/(newmax-newmin))*(max-min))+min
Prediksi3<-(((Predict3-newmin)/(newmax-newmin))*(max-min))+min

PEMILIHAN MODEL NN TERBAIK

MAPE_1<-MAPE(Prediksi1[1:84],Yt[3:86])*100
MAPE_2<-MAPE(Prediksi2[1:84],Yt[3:86])*100
MAPE_3<-MAPE(Prediksi3[1:84],Yt[3:86])*100
MAPE_1  
## [1] 527.8965
MAPE_2  #TERBAIK
## [1] 520.1451
MAPE_3 
## [1] 541.4504

Berdasarkan hasil output, dapat dilihat bahwa NN dengan 1 neuron, 2 neuron, dan 3 neuron mempunyai kemampuan peramalan yang buruk untuk meramalkan produksi kopi di Indonesia. Namun NN dengan 2 neuron lebih baik dibandingkan model lainnya, sehingga dipilih model terbaiknya yaitu NN 2 neuron.

Grafik Perbandingan Aktual vs Prediksi

ts.plot(data_train[1:84],col="black",type="o",main="Peramalan Produksi Kopi di Indonesia",xlab="Waktu",ylab="Produksi",
        lwd=1, xlim = c(0, length(data_train)))
lines(Prediksi2,col="blue",type="o",lwd=1)
legend('topright',legend=c("Aktual","Prediksi"),col=c("black","blue"),lty=1, cex = 1.05)

PREDIKSI 11 PERIODE SELANJUTNYA DENGAN NN 2 NEURON

Ytest=c(Y,rep(0,11))
k=11
for (i in 87:97)
{ 
  Xtest=t(matrix(c(Ytest[i-1]),byrow=FALSE))
  Ytest[i]=compute(nn2,covariate=Xtest)$net.result
}
Ytest=ts(Ytest[87:97])
Ytest
## Time Series:
## Start = 1 
## End = 11 
## Frequency = 1 
##  [1] -0.1774272 -0.2938434 -0.3867841 -0.4667031 -0.5364119 -0.5958455
##  [7] -0.6442900 -0.6816864 -0.7090536 -0.7281858 -0.7410938

DESTANDARISASI PERAMALAN

Peramalan<-(((Ytest-newmin)/(newmax-newmin))*(max-min))+min

GRAFIK PERBANDINGAN

Aktual<-Zt
Hasil_NN<-ts(c(NA,Prediksi2,Peramalan))
plot(Hasil_NN,ylim=c(0,10),main="",sub="Januari 2017 - Desember 2024",xlim=c(0,120),
     ylab="Produksi Kopi (Ribu Ton)", xlab="Periode",col="white",lty=1,lwd=2)
points(Hasil_NN[1:96],cex=1,col="orange",pch=19)
lines(Hasil_NN[1:96],col="orange",lwd=3,cex=1.5)
points(Aktual[1:86],cex=1,col="black",pch=10)
lines(Aktual[1:86],col="black",lwd=1.5,cex=1)
lines(87:96, Aktual[87:96], col="blue", cex=1, lwd=1.5)
points(87:96, Aktual[87:96], pch=10, col="blue")
legend("topleft",c("Data Train"),bty="n",cex=1,lty=1, col="black",lwd=3,bg="white")
legend("top",c("Data test"),bty="n",cex=1,lty=1, col="blue",lwd=3,bg="white")
legend("topright",c("NN 2 Neuron"),bty="n",cex=1,lty=1, col="orange",lwd=3,bg="white")
grid()

Berdasarkan gambar diatas dapat dilihat bahwa model NN mampu mengikuti pola musiman tahunan dengan cukup baik, terutama pada bagian data latih. Namun, terlihat bahwa pada beberapa puncak produksi (khususnya sekitar periode ke-20 dan ke-80), model cenderung meremehkan nilai aktual (underfitting). Pada data uji, meskipun model masih dapat mengikuti arah tren secara umum, terdapat deviasi antara nilai aktual dan hasil prediksi, yang menunjukkan bahwa model memiliki keterbatasan dalam menangkap fluktuasi mendadak. Secara keseluruhan, model NN dengan 3 neuron cukup mampu menangkap pola musiman, namun peningkatan jumlah neuron atau teknik lain mungkin diperlukan untuk meningkatkan akurasi prediksi.