Perbandingan Model TSR & ARIMAX Pada Data Luas Panen Padi
Kelompok 6 dengan Anggota: Mujahid Shobrun Jamil, Nisrina Nur Rahman, Novia Ardani, Rifa Joseph
2025-05-09
Pendahuluan
Kita sering membutuhkan analisis prediksi untuk membantu kita dalam membuat keputusan.Salah satu hal penting dalam melakukan prediksi kita memerlukan data lampau untuk memprediksi data masa depan. Prediksi ini sering disebut dengan peramalan. Peramalan banyak dibutuhkan dalam banyak situasi salah satunya pada masalah ketahanan pangan. Pemerintahan Indonesia menjadikan ketahanan pangan sebagai prioritas kebijakan jangka panjang, karena kekurangan gizi dan kerawanan pangan masih menjadi masalah di Indonesia.Dalam konteks ketahanan pangan nasional, padi memegang peranan penting sebagai komoditas pangan utama. Ketersediaan beras, yang dihasilkan dari tanaman padi, sangat bergantung pada luas panen yang tercapai setiap tahunnya. Oleh karena itu, perkembangan luas panen padi menjadi indikator penting dalam menentukan produksi padi nasional yang pada akhirnya turut memengaruhi stabilitas ketahanan pangan di Indonesia.Data yang digunakan adalah data luas panen padi di Indonesia.Data berupa data bulanan periode Januari 2018 hingga Desember 2024.Data dalam analisis ini, terdapat 3 variabel dummy dimana D1 merupakan variabel dummy 1 bulan sebelum musim panen,D2 merupakan variabel dummy saat musim panen, D3 merupakan variabel dummy 1 bulan setelah musim panen dan dalam analisis ini dibagi menjadi dua yaitu data in sampel dan out sampel. Data In sampel diambil mulai dari periode Januari 2018 hingga Agustus 2023 dan data out sampel dimulai dari periode September 2023 hingga Desember 2024.
Import Library
library(MASS)
library(hms)
library(car)## Loading required package: carDatalibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries)
library(timeSeries)## Loading required package: timeDate##
## Attaching package: 'timeSeries'## The following objects are masked from 'package:graphics':
##
## lines, pointslibrary(quadprog)
library(zoo)##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(fracdiff)
library(fUnitRoots)
library(lmtest)
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#Menginput Data
data<-read.csv(file.choose(),header=TRUE)
head(data)## LP t D1 D2 D3
## 1 540 1 0 0 0
## 2 1070 2 1 0 0
## 3 1770 3 0 1 0
## 4 1410 4 0 0 1
## 5 1000 5 0 0 0
## 6 890 6 0 0 080/100*84 #banyaknya pembagian## [1] 67.2#In Sample - Out Sample
Yt=data$LP[1:68]
t=data$t[1:68]
ts.plot(as.ts(Yt),lwd=5,ylab="Zt",col="red")
D1=data$D1[1:68]
D2=data$D2[1:68]
D3=data$D3[1:68]
Ytout=data$LP[69:84]#Pembentukan model TSR
#Regresi Runtun Waktu
model1=lm(Yt~t+D1+D2+D3)
summary(model1)##
## Call:
## lm(formula = Yt ~ t + D1 + D2 + D3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -463.28 -193.87 54.28 152.50 540.02
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 787.7361 65.2180 12.079 < 2e-16 ***
## t -0.9571 1.5588 -0.614 0.541
## D1 36.2238 108.9713 0.332 0.741
## D2 890.5142 108.9147 8.176 1.76e-11 ***
## D3 744.8046 108.8803 6.841 3.79e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 252 on 63 degrees of freedom
## Multiple R-squared: 0.6248, Adjusted R-squared: 0.6009
## F-statistic: 26.22 on 4 and 63 DF, p-value: 8.066e-13#Seleksi Variabel
model2=lm(Yt~D2+D3)
summary(model2)##
## Call:
## lm(formula = Yt ~ D2 + D3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.39 -195.89 41.61 164.11 541.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 758.39 33.28 22.786 < 2e-16 ***
## D2 888.27 106.99 8.302 8.48e-12 ***
## D3 741.61 106.99 6.932 2.29e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 249.1 on 65 degrees of freedom
## Multiple R-squared: 0.6217, Adjusted R-squared: 0.6101
## F-statistic: 53.42 on 2 and 65 DF, p-value: 1.898e-14Model akhir yang terbentuk adalah Yt = D2 + D3
#Prediksi & Peramalana TSR
TSR=758.39+888.27*data$D2+741.61*data$D3
TSR## [1] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [10] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39
## [19] 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39 1646.66
## [28] 1500.00 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39
## [37] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [46] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39
## [55] 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39 1646.66
## [64] 1500.00 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39
## [73] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [82] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39predinTSR=TSR[1:68]
predoutTSR=TSR[69:84]
ramalTSR=TSR[85:90]#Akurasi Model TSR
library(MLmetrics)
#Akurasi Insampel TSR
RMSEINTSR=RMSE(predinTSR,Yt)
RMSEINTSR## [1] 243.5126MAPEINTSR=mean(abs(Yt-(predinTSR))/(Yt))*100
MAPEINTSR## [1] 31.23996smape <- function(actual, forecast) {
return(mean(200 * abs(forecast - actual) / (abs(actual) + abs(forecast)), na.rm = TRUE))
}
SMAPEINTSR=smape(Yt,predinTSR)
SMAPEINTSR## [1] 25.8103#Akurasi Outsampel TSR
RMSEOUTTSR=RMSE(predoutTSR,Ytout)
RMSEOUTTSR## [1] 303.3278MAPEOUTTSR=mean(abs(Ytout-(predoutTSR))/(Ytout))*100
MAPEOUTTSR## [1] 45.17288smape <- function(actual, forecast) {
return(mean(200 * abs(forecast - actual) / (abs(actual) + abs(forecast)), na.rm = TRUE))
}
SMAPEOUTTSR=smape(Ytout,predoutTSR)
SMAPEOUTTSR## [1] 34.43411Pemodelan ARIMAX
#Residual dari model TSR
res=resid(model2)
res## 1 2 3 4 5 6
## -218.392857 311.607143 123.333333 -90.000000 241.607143 131.607143
## 7 8 9 10 11 12
## 311.607143 331.607143 251.607143 -128.392857 -238.392857 -388.392857
## 13 14 15 16 17 18
## -348.392857 -28.392857 73.333333 180.000000 161.607143 121.607143
## 19 20 21 22 23 24
## 191.607143 421.607143 41.607143 -158.392857 -288.392857 -438.392857
## 25 26 27 28 29 30
## -438.392857 -278.392857 -456.666667 360.000000 541.607143 -18.392857
## 31 32 33 34 35 36
## 131.607143 441.607143 261.607143 21.607143 -188.392857 -468.392857
## 37 38 39 40 41 42
## -348.392857 11.607143 143.333333 -40.000000 21.607143 41.607143
## 43 44 45 46 47 48
## 311.607143 91.607143 71.607143 -8.392857 -238.392857 -388.392857
## 49 50 51 52 53 54
## -288.392857 11.607143 113.333333 -80.000000 71.607143 111.607143
## 55 56 57 58 59 60
## 171.607143 51.607143 81.607143 31.607143 -148.392857 -398.392857
## 61 62 63 64 65 66
## -308.392857 181.607143 3.333333 -330.000000 211.607143 191.607143
## 67 68
## 71.607143 101.607143#Identifikasi Apakah Residual dari TSR Sudah WN
acf(res,68,lwd=1,col="blue")
x<-res
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box TSR",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
pv2## [1] 1.810882e-05 7.405289e-05 2.692312e-04 1.112853e-04 9.839038e-07
## [6] 2.605304e-09 4.156342e-12 4.009015e-13 1.149192e-12 1.448064e-12
## [11] 5.662137e-15 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [31] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [41] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [46] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [51] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [56] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [61] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [66] 0.000000e+00 0.000000e+00 NABerdasarkan residual dari model TSR diketahui Residual tidak memenuhi asumsi white noise, maka dilanjutkan dengan model ARIMAX
#Pemodelan ARIMA dari Data Residual TSR
ts.plot(res,ylab="Residual TSR",col="red",lwd=3)
acf(res,68,col="blue",lwd=1) #Membuat ACF
pacf(res,68,col="blue",lwd=1) #Membuat PACF
Berdasarkan grafik runtun waktu dari residual model TSR, diperoleh bahwa data tidak stasioner
#Stasioneritas dalam Variansi
min(res)## [1] -468.3929res=resid(model2) #Residual dari model TSR
res## 1 2 3 4 5 6
## -218.392857 311.607143 123.333333 -90.000000 241.607143 131.607143
## 7 8 9 10 11 12
## 311.607143 331.607143 251.607143 -128.392857 -238.392857 -388.392857
## 13 14 15 16 17 18
## -348.392857 -28.392857 73.333333 180.000000 161.607143 121.607143
## 19 20 21 22 23 24
## 191.607143 421.607143 41.607143 -158.392857 -288.392857 -438.392857
## 25 26 27 28 29 30
## -438.392857 -278.392857 -456.666667 360.000000 541.607143 -18.392857
## 31 32 33 34 35 36
## 131.607143 441.607143 261.607143 21.607143 -188.392857 -468.392857
## 37 38 39 40 41 42
## -348.392857 11.607143 143.333333 -40.000000 21.607143 41.607143
## 43 44 45 46 47 48
## 311.607143 91.607143 71.607143 -8.392857 -238.392857 -388.392857
## 49 50 51 52 53 54
## -288.392857 11.607143 113.333333 -80.000000 71.607143 111.607143
## 55 56 57 58 59 60
## 171.607143 51.607143 81.607143 31.607143 -148.392857 -398.392857
## 61 62 63 64 65 66
## -308.392857 181.607143 3.333333 -330.000000 211.607143 191.607143
## 67 68
## 71.607143 101.607143res=res+550
boxcox(res~1) 
p<-powerTransform(res)
p## Estimated transformation parameter
## res
## 1.006816#Stasioneritas dalam Rata-Rata ADF
library(urca)##
## Attaching package: 'urca'## The following objects are masked from 'package:fUnitRoots':
##
## punitroot, qunitroot, unitrootTabletest=ur.df(res)
summary(test)##
## ###############################################
## # 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
## -470.1 -126.1 36.3 155.8 823.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.08139 0.04855 -1.676 0.0986 .
## z.diff.lag -0.00572 0.12056 -0.047 0.9623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 232.6 on 64 degrees of freedom
## Multiple R-squared: 0.04457, Adjusted R-squared: 0.01471
## F-statistic: 1.493 on 2 and 64 DF, p-value: 0.2325
##
##
## Value of test-statistic is: -1.6763
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.6 -1.95 -1.61diff1 = diff(res, differences = 2)
head(diff1)## 3 4 5 6 7 8
## -718.27381 -25.05952 544.94048 -441.60714 290.00000 -160.00000adf.test(diff1)## Warning in adf.test(diff1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff1
## Dickey-Fuller = -5.9066, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryacf(diff1,68,col="blue",lwd=1) #Membuat ACF
pacf(diff1,68,col="blue",lwd=1) #Membuat PACF
#Pemodelan ARIMA
fit1=arima(x=res,order=c(0,2,1))
coeftest(fit1) #Signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.999994 0.038889 -25.714 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit2=arima(x=res,order=c(0,2,2))
coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.10383 0.20604 -5.3575 8.44e-08 ***
## ma2 0.10383 0.20231 0.5132 0.6078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit3=arima(x=res,order=c(0,2,3))
coeftest(fit3) #Signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.33455 0.10179 -13.1110 < 2.2e-16 ***
## ma2 -0.31829 0.16145 -1.9714 0.04867 *
## ma3 0.66070 0.09019 7.3257 2.376e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit4=arima(x=res,order=c(1,2,0))
coeftest(fit4) #Sigifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.34834 0.11864 -2.936 0.003325 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit5=arima(x=res,order=c(1,2,1))
coeftest(fit5)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.036735 0.127558 -0.288 0.7734
## ma1 -0.999999 0.038979 -25.655 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit6=arima(x=res,order=c(1,2,2))
coeftest(fit6)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.59956 0.26948 -2.2249 0.02609 *
## ma1 -0.24653 0.21512 -1.1460 0.25179
## ma2 -0.75347 0.21362 -3.5272 0.00042 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit7=arima(x=res,order=c(1,2,3))
coeftest(fit7)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.18302 0.18003 1.0166 0.3093
## ma1 -1.42674 0.15248 -9.3570 < 2.2e-16 ***
## ma2 -0.13361 0.28180 -0.4741 0.6354
## ma3 0.56822 0.14458 3.9302 8.488e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit8=arima(x=res,order=c(2,2,0))
coeftest(fit8) #Signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.55115 0.10576 -5.2115 1.873e-07 ***
## ar2 -0.56333 0.10326 -5.4556 4.881e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit9=arima(x=res,order=c(2,2,1))
coeftest(fit9)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.050588 0.119697 -0.4226 0.672560
## ar2 -0.342031 0.118610 -2.8837 0.003931 **
## ma1 -0.999999 0.040542 -24.6658 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit10=arima(x=res,order=c(2,2,2))
coeftest(fit10)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.36683 0.26209 -1.3997 0.161617
## ar2 -0.35343 0.12023 -2.9397 0.003286 **
## ma1 -0.63802 0.26864 -2.3750 0.017548 *
## ma2 -0.36196 0.26604 -1.3605 0.173663
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit11=arima(x=res,order=c(2,2,3))
coeftest(fit11)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.00581 0.15302 -6.5730 4.930e-11 ***
## ar2 -0.85860 0.11696 -7.3408 2.123e-13 ***
## ma1 0.12384 0.19807 0.6252 0.531835
## ma2 -0.36371 0.12224 -2.9754 0.002926 **
## ma3 -0.76013 0.19345 -3.9294 8.516e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Pemeriksaan Diagnostik
res1=resid(fit1)
res3=resid(fit3)
res4=resid(fit4)
res8=resid(fit8)
#White Noise
x<-res1
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
pv2## [1] 0.98019050 0.02937485 0.05221548 0.10242136 0.08719959 0.14225353
## [7] 0.15540206 0.20352561 0.27336266 0.32557213 0.37308150 0.11106943
## [13] 0.14594580 0.18669180 0.21322365 0.11625462 0.14344609 0.18343967
## [19] 0.22717501 0.24545196 0.28956530 0.34191892 0.28388259 0.24675016
## [25] 0.26936904 0.31585886 0.31752569 0.24587944 0.25339709 0.24129817
## [31] 0.27531544 0.17229965 0.17111386 0.13046864 0.07079262 0.08190493
## [37] 0.03227560 0.04011869 0.04766834 0.05155209 0.06020203 0.06820081
## [43] 0.08287040 0.08625002 0.10238309 0.12162089 0.12543109 0.13127874
## [49] 0.14813162 0.16975241 0.19512981 0.21535305 0.24342387 0.27028846
## [55] 0.29256412 0.32554673 0.31568293 0.33007921 0.20841023 0.14400400
## [61] 0.16476916 0.09474894 0.10778288 0.12160327 0.13865254 0.15838475
## [67] 0.17976028 NAx<-res3
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
x<-res4
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
x<-res8
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
#Cek Normalitas
library(nortest)
lillie.test(res1)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: res1
## D = 0.071771, p-value = 0.5211lillie.test(res3)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: res3
## D = 0.065702, p-value = 0.6618lillie.test(res4)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: res4
## D = 0.056773, p-value = 0.8484lillie.test(res8)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: res8
## D = 0.069812, p-value = 0.5662Setelah melakukan pengujian signifikansi parameter dan pengujian diagnostik didapatkan model ARIMA terbaik adalah model ARIMA (2,2,0)
Pemodelan ARIMAX
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:timeDate':
##
## kurtosis, skewness## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarxreg<-data.frame(D2,D3)
Yt=data$LP[1:68]
modelarimax=arima(Yt,order=c(2,2,0),xreg=xreg)
coeftest(modelarimax)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.46618 0.14681 -3.1754 0.001496 **
## ar2 -0.51959 0.11636 -4.4652 8.000e-06 ***
## D2 823.23968 101.58947 8.1036 5.336e-16 ***
## D3 619.48787 96.36059 6.4289 1.286e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Pemodelan ARIMAX (2,2,0)
#Pengecekkan Diagnostik
#WNhite Noise
x<-resid(modelarimax)
pv2<-rep(1,length(x))
for(i in 1:length(x))
pv2[i]=Box.test(x,lag=i,type="Ljung")$p.value
plot(pv2,type="p",main="Ljung-Box",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")
#Normalitas
lillie.test(x)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x
## D = 0.083837, p-value = 0.2778Setelah melakukan pengujian signifikansi parameter dan pengujian diagnostik didapatkan model ARIMAX terbaik adalah model ARIMA (2,2,0)
#Prediksi & Peramalan Model ARIMAX
#Prediksi In Sample
fits.ARIMA=as.ts(fitted(modelarimax))
Pred_in=fits.ARIMA
Pred_in## Time Series:
## Start = 1
## End = 68
## Frequency = 1
## [1] 539.75850 1069.53937 2301.22521 1609.43106 989.06725 1056.13846
## [7] 738.90812 1280.80866 1033.90898 1059.75241 441.81390 440.00668
## [13] 98.35868 382.20859 1643.98656 1608.96182 1305.28790 922.89397
## [19] 951.23467 916.49492 1278.25595 621.23808 633.03540 213.84135
## [25] 142.95254 260.46419 1310.71213 917.25439 1796.11647 1226.25591
## [31] 891.87671 1030.88777 1066.50335 985.29603 822.56899 377.18975
## [37] 27.04523 379.89762 1633.52050 1734.40716 949.66259 856.67403
## [43] 748.31090 1181.62095 728.53314 971.36120 594.05349 391.10283
## [49] 260.64345 411.88704 1670.10601 1681.20515 874.75129 939.66389
## [55] 818.98500 975.21435 763.52133 893.59821 699.35655 532.17091
## [61] 210.17926 417.86865 1890.10636 1406.39400 663.69176 1149.84421
## [67] 773.38624 984.97086#Prediksi Out Sample
Ytnew=data$LP[69:84]
D2new=data$D2[69:84]
D3new=data$D3[69:84]
xregnewbgt=data.frame(D2new,D3new)
fore.ARIMA=predict(modelarimax,newYt=Ytnew,newxreg=xregnewbgt,16)$pred
se.fore.ARIMA=predict(modelarimax,newxreg=xregnewbgt,16)$se
Pred_out=fore.ARIMA
Pred_out## Time Series:
## Start = 69
## End = 84
## Frequency = 1
## [1] 872.0312 814.5010 798.7356 799.6433 771.0779 747.5892 1560.2874
## [8] 1337.3205 695.9338 679.7928 662.3620 642.5408 624.5041 606.8776
## [15] 588.1326 569.6959#Forecast
Ytnewnew=data$Zt[85:90]
D2newnew=data$D2[85:90]
D3newnew=data$D3[85:90]
xregnewbgtbgt=data.frame (D2newnew,D3newnew)
fore.ARIMAA=predict(modelarimax,NewYt=Ytnewnew,newxreg=xregnewbgtbgt,6)$pred
se.fore.ARIMAA=predict(modelarimax,newxreg=xregnewbgtbgt,6)$se
Forecastt=fore.ARIMAA
Forecastt## Time Series:
## Start = 69
## End = 74
## Frequency = 1
## [1] 872.0312 814.5010 1621.9753 1419.1312 771.0779 747.5892se_foree=se.fore.ARIMAA#Akurasi Model ARIMAX
library(MLmetrics)
#Akurasi Insampel ARIMAX
RMSEINARIMAX=RMSE(Pred_in,Yt)
RMSEINARIMAX## [1] 262.1192MAPEINARIMAX=mean(abs(Yt-(Pred_in))/(Yt))*100
MAPEINARIMAX## [1] 25.42295smape <- function(actual, forecast) {
return(mean(200 * abs(forecast - actual) / (abs(actual) + abs(forecast)), na.rm = TRUE))
}
SMAPEINARIMAX=smape(Yt,Pred_in)
SMAPEINARIMAX## [1] 28.13487#Akurasi Outsampel ARIMAX
RMSEOUTARIMAX=RMSE(Pred_out,Ytout)
RMSEOUTARIMAX## [1] 322.8926MAPEOUTARIMAX=mean(abs(Ytout-(Pred_out))/(Ytout))*100
MAPEOUTARIMAX## [1] 44.5743smape <- function(actual, forecast) {
return(mean(200 * abs(forecast - actual) / (abs(actual) + abs(forecast)), na.rm = TRUE))
}
SMAPEOUTARIMAX=smape(Ytout,Pred_out)
SMAPEOUTARIMAX## [1] 36.14227Plot Perbandingan IN SAMPLE, OUT SAMPLE, PREDIKSI ARIMAX, PERAMALAN ARIMAX
#Model ARIMAX Dengan Data Aktual
prediksii=c(Pred_in,Pred_out,Forecastt)
prediksii## [1] 539.75850 1069.53937 2301.22521 1609.43106 989.06725 1056.13846
## [7] 738.90812 1280.80866 1033.90898 1059.75241 441.81390 440.00668
## [13] 98.35868 382.20859 1643.98656 1608.96182 1305.28790 922.89397
## [19] 951.23467 916.49492 1278.25595 621.23808 633.03540 213.84135
## [25] 142.95254 260.46419 1310.71213 917.25439 1796.11647 1226.25591
## [31] 891.87671 1030.88777 1066.50335 985.29603 822.56899 377.18975
## [37] 27.04523 379.89762 1633.52050 1734.40716 949.66259 856.67403
## [43] 748.31090 1181.62095 728.53314 971.36120 594.05349 391.10283
## [49] 260.64345 411.88704 1670.10601 1681.20515 874.75129 939.66389
## [55] 818.98500 975.21435 763.52133 893.59821 699.35655 532.17091
## [61] 210.17926 417.86865 1890.10636 1406.39400 663.69176 1149.84421
## [67] 773.38624 984.97086 872.03120 814.50101 798.73558 799.64334
## [73] 771.07787 747.58916 1560.28736 1337.32045 695.93376 679.79275
## [79] 662.36198 642.54079 624.50411 606.87756 588.13260 569.69593
## [85] 872.03120 814.50101 1621.97526 1419.13121 771.07787 747.58916plot(prediksii,xlim=c(0,90),ylim=c(0,3000),main="Plot Data Padi ARIMAX",xlab= "Waktu",ylab="Jumlah Padi",col="white")
lines(data$LP,col="black",lwd=4)
points(data$LP[1:68],cex=1,col="black",pch=19)
lines(69:84, data$LP[69:84],col="darkgrey",lwd=4)
points(69:84, data$LP[69:84],cex=1,col="darkgrey",pch=19)
lines(prediksii[1:69],col="coral2",lwd=4)
points(prediksii[1:68],cex=1,col="coral2",pch=19)
lines(69:85, prediksii[69:85], col="green", lwd=4, cex=1.5)
points(69:84, prediksii[69:84],cex=1,col="green",pch=19)
points(85:90, Forecastt,cex=1,col="blue",pch=19)
lines(85:90, Forecastt,col="blue",lwd=4,cex=1.5)
grid()
legend("topright",legend=c("In Sample","Out Sample","Prediksi In Sample","Prediksi Out Sample","Forecasting"),cex=0.5,lty=1,col=c("black","darkgrey","coral2","green","blue"))
grid()
#Model TSR Dengan Data Aktual
prediksi2=c(predinTSR,predoutTSR,ramalTSR)
prediksi2## [1] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [10] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39
## [19] 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39 1646.66
## [28] 1500.00 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39
## [37] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [46] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39
## [55] 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39 1646.66
## [64] 1500.00 758.39 758.39 758.39 758.39 758.39 758.39 758.39 758.39
## [73] 758.39 758.39 1646.66 1500.00 758.39 758.39 758.39 758.39 758.39
## [82] 758.39 758.39 758.39 758.39 758.39 1646.66 1500.00 758.39 758.39plot(prediksi2,xlim=c(0,90),ylim=c(0,3000),main="Plot Data Padi TSR",xlab= "Waktu",ylab="Jumlah Padi",col="white")
lines(data$LP,col="black",lwd=4)
points(data$LP[1:69],cex=1,col="black",pch=19)
lines(69:84, data$LP[69:84],col="palegreen3",lwd=4)
points(69:84, data$LP[69:84],cex=1,col="palegreen3",pch=19)
lines(TSR[1:69],col="orange",lwd=4)
points(TSR[1:68],cex=1,col="orange",pch=19)
lines(69:85, TSR[69:85], col="ivory3", lwd=4, cex=1.5)
points(69:84, TSR[69:84],cex=1,col="ivory3",pch=19)
points(85:90, TSR[85:90],cex=1,col="skyblue1",pch=19)
lines(85:90, TSR[85:90],col="skyblue1",lwd=4,cex=1.5)
grid()
legend("topright",legend=c("In Sample","Out Sample","Prediksi In Sample","Prediksi Out Sample","Forecasting"),cex=0.5,lty=1,col=c("black","palegreen3","orange","ivory3","skyblue1"))
grid()
Kesimpulan:
Berdasarkan hasil analisis diketehaui bahwa model ARIMAX lebih baik dibandingkan model TSR (Time Series Reggression). Hal ini dapat dilihat dari perbandingan grafik model ARIMAX dan model TSR dimana pola data prediksi ARIMAX cenderung lebih mengikuti pola aktual daripada pola data prediksi TSR. Selain itu, nilai akurasi model ARIMAX menghasilkan nilai yang lebih kecil dibandingkan nilai akurasi model TSR.