####-----------------ARIMAX-----------------####
library(MASS)
library(hms)
library(car)
## Loading required package: carData
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tseries)
library(timeSeries)
## Loading required package: timeDate
##
## Attaching package: 'timeSeries'
## The following objects are masked from 'package:graphics':
##
## lines, points
library(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.numeric
library(fracdiff)
library(fUnitRoots)
library(lmtest)
library(Metrics)
##
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
##
## accuracy
Ayam<-read.csv("Data TM.csv", sep = ";")
Ayam$Minggu <- as.Date(Ayam$Minggu, format="%d/%m/%Y")
head(Ayam);tail(Ayam)
## Minggu Harga D1 D2 D3
## 1 2019-05-06 36650 0 0 0
## 2 2019-05-13 35400 0 0 0
## 3 2019-05-20 34200 0 0 0
## 4 2019-05-27 34650 1 0 0
## 5 2019-06-03 37350 0 1 0
## 6 2019-06-10 37150 0 0 1
## Minggu Harga D1 D2 D3
## 258 2024-04-08 38900 0 1 0
## 259 2024-04-15 38850 0 0 1
## 260 2024-04-22 39400 0 0 0
## 261 2024-04-29 38150 0 0 0
## 262 2024-05-06 38550 0 0 0
## 263 2024-05-13 39150 0 0 0
Zt = as.numeric(Ayam$Harga)
Zt_zoo <- zoo(Zt, order.by = Ayam$Minggu)
plot(Zt_zoo, lwd = 5, col = "black", ylab = "Harga Ayam", xlab = "Tanggal")

#In Sample(210) & Out Sample(53)
Insample=Ayam$Harga[1:210]
t=c(1:210)
Insample_zoo = zoo(Insample, order.by = Ayam$Minggu)
plot(Insample_zoo, lwd = 5, col = "green", ylab = "Insample", xlab = "Tanggal")

D1=Ayam$D1[1:210]
D2=Ayam$D2[1:210]
D3=Ayam$D3[1:210]
####-----------------Time Series Regression-----------------####
model1=lm(Insample~t+D1+D2+D3)
summary(model1)
##
## Call:
## lm(formula = Insample ~ t + D1 + D2 + D3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5623.9 -1336.4 152.1 1155.6 12005.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33309.19 301.73 110.394 < 2e-16 ***
## t 13.74 2.46 5.587 7.3e-08 ***
## D1 -259.57 978.53 -0.265 0.79107
## D2 2290.23 895.80 2.557 0.01129 *
## D3 3860.19 978.54 3.945 0.00011 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2160 on 205 degrees of freedom
## Multiple R-squared: 0.2081, Adjusted R-squared: 0.1927
## F-statistic: 13.47 on 4 and 205 DF, p-value: 9.161e-10
#Seleksi Variabel
model2=lm(Insample~t+D2+D3)
summary(model2)
##
## Call:
## lm(formula = Insample ~ t + D2 + D3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5617.3 -1329.9 158.6 1162.1 12012.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33302.598 300.026 110.999 < 2e-16 ***
## t 13.745 2.454 5.600 6.81e-08 ***
## D2 2296.750 893.437 2.571 0.010855 *
## D3 3866.712 976.023 3.962 0.000103 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2156 on 206 degrees of freedom
## Multiple R-squared: 0.2078, Adjusted R-squared: 0.1963
## F-statistic: 18.02 on 3 and 206 DF, p-value: 2.017e-10
res=resid(model2) #Residual dari model TSR
#Identifikasi Apakah Residual dari TSR Sudah WN
acf(res,100,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")

#Pemodelan ARIMA dari Ayam Residual TSR
resZoo = zoo(res, order.by = Ayam$Minggu)
plot(resZoo,ylab="Residual TSR",col="green",lwd=3)

acf(res,100,col="blue",lwd=1) #Membuat ACF

pacf(res,100,col="red",lwd=1) #Membuat PACF

#Stasioneritas dalam Variansi
min(res)
## [1] -5617.326
res1 =res+min(Zt)
boxcox(res1~1) #belum stasioner dalam varians

p<-powerTransform(res1);p
## Estimated transformation parameter
## res1
## -0.3233274
lambda = p$lambda
res2 = res1^lambda
boxcox(res2~1);powerTransform(res2) #sudah stasioner dalam varians setelah pangkat lambda

## Estimated transformation parameter
## res2
## 0.9999999
#Stasioneritas dalam Rata-Rata ADF
library(urca)
##
## Attaching package: 'urca'
## The following objects are masked from 'package:fUnitRoots':
##
## punitroot, qunitroot, unitrootTable
adfTest(res2)#belum stasioner dalam rata-rata
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: 0.0173
## P VALUE:
## 0.621
##
## Description:
## Mon Apr 21 22:46:56 2025 by user: mrio1
datadiff=diff(res2,differences=1)
adfTest(datadiff)#sudah stasioner dalam rata-rata setelah differencing 1 kali (d=1)
## 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: -10.4663
## P VALUE:
## 0.01
##
## Description:
## Mon Apr 21 22:46:56 2025 by user: mrio1
acf(datadiff,50,col="blue",lwd=1);grid() #Membuat ACF (signifikan pada lag 7 -> q=7)

pacf(datadiff,50,col="red",lwd=1);grid() #Membuat PACF (signifikan pada lag 7 -> p=7))

#Karena p=7, d=1, dan q=7 akan sangat banyak menghasilkan model sementara (tidak parsimoni), maka kita akan lakukan differencing 2 kali
datadiff2 = diff(res2,differences=2)
adfTest(datadiff2)#masih stasioner dalam rata-rata setelah differencing 2 kali (d=2)
## 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: -17.8218
## P VALUE:
## 0.01
##
## Description:
## Mon Apr 21 22:46:56 2025 by user: mrio1
acf(datadiff2,50,col="blue",lwd=1);grid() #Membuat ACF (cut-off setelah lag 1 -> q=1)

pacf(datadiff2,50,col="red",lwd=1);grid() #Membuat PACF (cut-off setelah lag 4 -> p=4))

####-----------------Pemodelan ARIMA-----------------####
fit1=arima(x=res2,order=c(0,0,1))
fit2=arima(x=res2,order=c(1,0,0))
fit3=arima(x=res2,order=c(2,0,0))
fit4=arima(x=res2,order=c(3,0,0))
fit5=arima(x=res2,order=c(4,0,0))
fit6=arima(x=res2,order=c(1,0,1))
fit7=arima(x=res2,order=c(2,0,1))
fit8=arima(x=res2,order=c(3,0,1))
fit9=arima(x=res2,order=c(4,0,1))
fit10=arima(x=res2,order=c(0,1,1))
fit11=arima(x=res2,order=c(1,1,0))
fit12=arima(x=res2,order=c(2,1,0))
fit13=arima(x=res2,order=c(3,1,0))
fit14=arima(x=res2,order=c(4,1,0))
fit15=arima(x=res2,order=c(1,1,1))
fit16=arima(x=res2,order=c(2,1,1))
fit17=arima(x=res2,order=c(3,1,1))
fit18=arima(x=res2,order=c(4,1,1))
fit19=arima(x=res2,order=c(0,2,1))
fit20=arima(x=res2,order=c(1,2,0))
fit21=arima(x=res2,order=c(2,2,0))
fit22=arima(x=res2,order=c(3,2,0))
fit23=arima(x=res2,order=c(4,2,0))
fit24=arima(x=res2,order=c(1,2,1))
fit25=arima(x=res2,order=c(2,2,1))
fit26=arima(x=res2,order=c(3,2,1))
fit27=arima(x=res2,order=c(4,2,1))
coeftest(fit1) #signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 6.4982e-01 3.7949e-02 17.124 < 2.2e-16 ***
## intercept 3.6359e-02 9.9459e-05 365.569 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit2) #signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.83154945 0.03807257 21.841 < 2.2e-16 ***
## intercept 0.03633027 0.00020671 175.755 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.86491508 0.06897401 12.54 <2e-16 ***
## ar2 -0.04008492 0.06910762 -0.58 0.5619
## intercept 0.03633115 0.00019994 181.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit4)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.86306869 0.06882697 12.5397 <2e-16 ***
## ar2 0.01771051 0.09186883 0.1928 0.8471
## ar3 -0.06765379 0.07112805 -0.9512 0.3415
## intercept 0.03633096 0.00018917 192.0569 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit5)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.85692414 0.06872548 12.4688 <2e-16 ***
## ar2 0.01983587 0.09158564 0.2166 0.8285
## ar3 0.00969778 0.09348669 0.1037 0.9174
## ar4 -0.08996008 0.07091313 -1.2686 0.2046
## intercept 0.03633156 0.00017592 206.5285 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit6)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.81861550 0.04655657 17.5832 <2e-16 ***
## ma1 0.04213848 0.07792850 0.5407 0.5887
## intercept 0.03633021 0.00020104 180.7089 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit7) #signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.70523100 0.09279334 18.3767 < 2.2e-16 ***
## ar2 -0.75554163 0.07619300 -9.9162 < 2.2e-16 ***
## ma1 -0.81533559 0.10079059 -8.0894 5.996e-16 ***
## intercept 0.03633521 0.00014025 259.0729 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit8)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.59858710 0.13783280 11.5980 < 2.2e-16 ***
## ar2 -0.58032340 0.17446294 -3.3263 0.0008799 ***
## ar3 -0.08429268 0.07490880 -1.1253 0.2604744
## ma1 -0.75930273 0.12242014 -6.2024 5.56e-10 ***
## intercept 0.03633445 0.00013912 261.1814 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit9)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.50374882 0.17801568 8.4473 < 2.2e-16 ***
## ar2 -0.54088241 0.18804172 -2.8764 0.004022 **
## ar3 0.03627858 0.13026164 0.2785 0.780624
## ar4 -0.08917888 0.08011077 -1.1132 0.265625
## ma1 -0.67024365 0.16740600 -4.0037 6.236e-05 ***
## intercept 0.03633526 0.00013918 261.0589 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit10)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.048463 0.070781 -0.6847 0.4935
coeftest(fit11)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.046246 0.069142 -0.6689 0.5036
coeftest(fit12)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.047755 0.069267 -0.6894 0.4905
## ar2 -0.024658 0.071192 -0.3464 0.7291
coeftest(fit13)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.047905 0.069279 -0.6915 0.4893
## ar2 -0.025024 0.071275 -0.3511 0.7255
## ar3 -0.007571 0.071136 -0.1064 0.9152
coeftest(fit14)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.0477729 0.0691768 -0.6906 0.4898
## ar2 -0.0248412 0.0711303 -0.3492 0.7269
## ar3 -0.0046801 0.0711297 -0.0658 0.9475
## ar4 0.0547956 0.0710108 0.7717 0.4403
coeftest(fit15) #Signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.840395 0.039178 21.451 < 2.2e-16 ***
## ma1 -1.000000 0.013628 -73.379 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit16)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.119727 1.232287 -0.0972 0.9226
## ar2 -0.027029 0.078486 -0.3444 0.7306
## ma1 0.071890 1.232188 0.0583 0.9535
coeftest(fit17)
## Warning in sqrt(diag(se)): NaNs produced
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.0240798 NaN NaN NaN
## ar2 -0.0238348 0.0374529 -0.6364 0.5245
## ar3 -0.0057087 NaN NaN NaN
## ma1 -0.0236599 NaN NaN NaN
coeftest(fit18)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.782758 0.270500 -2.8937 0.003807 **
## ar2 -0.060994 0.090977 -0.6704 0.502578
## ar3 -0.022388 0.090622 -0.2471 0.804869
## ar4 0.057241 0.076599 0.7473 0.454894
## ma1 0.740488 0.263299 2.8123 0.004918 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit19)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.000000 0.012788 -78.201 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit20)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.524227 0.059412 -8.8236 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit21)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.699992 0.065019 -10.7659 < 2.2e-16 ***
## ar2 -0.349991 0.066319 -5.2774 1.311e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit22)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.802484 0.066218 -12.1187 < 2.2e-16 ***
## ar2 -0.562599 0.079245 -7.0995 1.252e-12 ***
## ar3 -0.302637 0.067882 -4.4583 8.263e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit23)#signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.843466 0.068671 -12.2827 < 2.2e-16 ***
## ar2 -0.640147 0.087473 -7.3182 2.512e-13 ***
## ar3 -0.417127 0.088014 -4.7393 2.144e-06 ***
## ar4 -0.141916 0.070512 -2.0127 0.04415 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit24)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.041765 0.069460 -0.6013 0.5477
## ma1 -1.000000 0.012867 -77.7201 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit25)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.043100 0.069601 -0.6192 0.5358
## ar2 -0.019922 0.071558 -0.2784 0.7807
## ma1 -0.999997 0.012884 -77.6126 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit26)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.0431863 0.0696207 -0.6203 0.5351
## ar2 -0.0200722 0.0716615 -0.2801 0.7794
## ar3 -0.0028246 0.0714991 -0.0395 0.9685
## ma1 -0.9999948 0.0128876 -77.5933 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit27)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.04273295 0.06954320 -0.6145 0.5389
## ar2 -0.01954246 0.07154127 -0.2732 0.7847
## ar3 0.00064017 0.07153929 0.0089 0.9929
## ar4 0.05981555 0.07138585 0.8379 0.4021
## ma1 -0.99999992 0.01281326 -78.0441 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Pemeriksaan Diagnostik
resid1=resid(fit1)
resid2=resid(fit2)
resid7=resid(fit7)
resid15=resid(fit15)
resid19=resid(fit19)
resid20=resid(fit20)
resid21=resid(fit21)
resid22=resid(fit22)
resid23=resid(fit23)
#WN
x<-resid1
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 Resid1",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")

x<-resid2
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 Resid2",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

x<-resid7
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 Resid7",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

x<-resid15
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 Resid15",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

x<-resid19
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 Resid19",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

x<-resid20
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 Resid20",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")

x<-resid21
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 Resid21",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red")

x<-resid22
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 Resid22",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

x<-resid23
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 Resid23",xlab="lag",ylab="p value",ylim=c(0,1))
abline(h=0.05,lty=1,col="red") #Lolos

#Cek Normalitas (2,7,15,19,22,23)
ks.test(resid2,"pnorm",mean(resid2),sqrt(var(resid2))) #Normal
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid2
## D = 0.099146, p-value = 0.03221
## alternative hypothesis: two-sided
ks.test(resid7,"pnorm",mean(resid7),sqrt(var(resid7)))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid7
## D = 0.12138, p-value = 0.004108
## alternative hypothesis: two-sided
ks.test(resid15,"pnorm",mean(resid7),sqrt(var(resid15))) #Normal
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid15
## D = 0.10983, p-value = 0.01261
## alternative hypothesis: two-sided
ks.test(resid19,"pnorm",mean(resid19),sqrt(var(resid19))) #Normal
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid19
## D = 0.10045, p-value = 0.02887
## alternative hypothesis: two-sided
ks.test(resid22,"pnorm",mean(resid22),sqrt(var(resid22))) #Normal
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid22
## D = 0.10658, p-value = 0.01694
## alternative hypothesis: two-sided
ks.test(resid23,"pnorm",mean(resid23),sqrt(var(resid23))) #Normal
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: resid23
## D = 0.098813, p-value = 0.03312
## alternative hypothesis: two-sided
##Akurasi (2,15,19,22,23)
#Vektor Predict dan Vektor Akurasi
Predict = cbind(data.frame(fitted(fit2)),
data.frame(fitted(fit15)),
data.frame(fitted(fit19)),
data.frame(fitted(fit22)),
data.frame(fitted(fit23))
)
SMAPE = data.frame(SMAPE=rep(0,5))
MAPE = data.frame(MAPE=rep(0,5))
RMSE = data.frame(RMSE=rep(0,5))
Akurasi = cbind(SMAPE,MAPE,RMSE)
rownames(Akurasi) <- c("ARIMA(1,0,0)",
"ARIMA(1,1,1)",
"ARIMA(0,2,1)",
"ARIMA(3,2,0)",
"ARIMA(4,2,0)")
#SMAPE
for (i in 1:5){
Akurasi$SMAPE[i]=smape(actual = (res2),predicted = as.numeric(unlist(Predict[i])))
};Akurasi
## SMAPE MAPE RMSE
## ARIMA(1,0,0) 0.008719704 0 0
## ARIMA(1,1,1) 0.008684834 0 0
## ARIMA(0,2,1) 0.008763802 0 0
## ARIMA(3,2,0) 0.009940301 0 0
## ARIMA(4,2,0) 0.009800283 0 0
#MAPE
for (i in 1:5){
Akurasi$MAPE[i]=mape(actual = (res2),predicted = as.numeric(unlist(Predict[i])))
};Akurasi
## SMAPE MAPE RMSE
## ARIMA(1,0,0) 0.008719704 0.008736206 0
## ARIMA(1,1,1) 0.008684834 0.008696591 0
## ARIMA(0,2,1) 0.008763802 0.008782099 0
## ARIMA(3,2,0) 0.009940301 0.009949090 0
## ARIMA(4,2,0) 0.009800283 0.009812094 0
#RMSE
for (i in 1:5){
Akurasi$RMSE[i]=rmse(actual = (res2),predicted = as.numeric(unlist(Predict[i])))
};Akurasi
## SMAPE MAPE RMSE
## ARIMA(1,0,0) 0.008719704 0.008736206 0.0004862522
## ARIMA(1,1,1) 0.008684834 0.008696591 0.0004862640
## ARIMA(0,2,1) 0.008763802 0.008782099 0.0005068275
## ARIMA(3,2,0) 0.009940301 0.009949090 0.0005595840
## ARIMA(4,2,0) 0.009800283 0.009812094 0.0005541118
#Didapatkan Model terbaik adalah ARIMA(1,1,1)
####-----------------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':
##
## tar
#Pemodelan ARIMAX (1,1,1)
xreg<-data.frame(t,D2,D3)
modelarimax=arima(Insample,order=c(1,1,1),xreg=xreg)
coeftest(modelarimax)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.832731 0.040046 20.7942 < 2.2e-16 ***
## ma1 -0.999998 0.013794 -72.4942 < 2.2e-16 ***
## t 11.856624 7.891548 1.5024 0.132982
## D2 885.553765 489.746874 1.8082 0.070577 .
## D3 1354.832358 490.848966 2.7602 0.005777 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Cek Diagnostik
#WN
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
ks.test(x,"pnorm",mean(x),sqrt(var(x)))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: x
## D = 0.10428, p-value = 0.02077
## alternative hypothesis: two-sided
#Prediksi In Sample
Pred_in=as.ts(fitted(modelarimax))
Pred_in
## Time Series:
## Start = 1
## End = 210
## Frequency = 1
## [1] 36613.36 36607.94 35465.53 34426.87 35753.48 37769.61 35801.09 34688.38
## [9] 33581.90 32967.73 34553.13 36421.86 35797.42 34580.62 33776.43 34283.93
## [17] 32509.43 31630.60 31058.33 31020.55 32194.94 32180.39 31973.21 32827.62
## [25] 33679.20 34144.61 34142.62 33975.09 34069.58 34706.73 34791.41 34576.47
## [33] 34116.90 34378.62 34929.22 34882.26 34244.71 33688.57 33049.44 32666.80
## [41] 32711.11 33267.79 34318.10 34099.96 33295.13 32874.62 33121.57 32814.69
## [49] 31826.79 30220.26 30159.96 29307.71 29258.22 29723.86 32133.68 34716.32
## [57] 38305.27 37777.57 37957.26 37917.10 36758.74 37656.82 38488.19 36085.89
## [65] 33778.78 33575.53 34401.25 32542.03 30905.72 31480.60 30559.41 31562.10
## [73] 31772.65 32439.56 32266.92 31888.25 32099.64 32770.28 33028.05 34081.29
## [81] 34798.44 35136.90 35097.61 34930.21 34435.12 35237.82 35828.37 35791.19
## [89] 35794.02 35502.05 34999.38 34411.22 33703.82 34293.35 34671.91 34630.59
## [97] 34341.94 34804.75 34848.59 34765.40 34436.77 35244.79 37173.38 36881.24
## [105] 35922.71 36691.13 38058.29 37193.94 36992.88 36705.09 36712.21 37177.20
## [113] 37095.62 36424.84 35167.87 34122.57 33953.34 33447.22 32569.48 32944.40
## [121] 32900.54 32693.59 33449.60 34330.54 34834.35 34964.02 35509.22 35345.73
## [129] 35642.12 35895.59 35814.97 35817.90 35818.70 35318.47 35070.84 35364.68
## [137] 35159.08 35496.35 35670.05 36680.37 37649.07 38366.11 38620.76 38077.89
## [145] 36740.94 35739.99 35532.44 35326.21 35537.91 35751.57 36422.64 36342.08
## [153] 36349.02 37438.20 37191.04 37956.55 38393.63 39767.76 38948.10 45575.67
## [161] 37183.48 36605.36 37655.06 38157.95 37366.04 36742.16 36623.31 37709.17
## [169] 36625.27 36333.68 35831.54 34954.50 34623.58 35417.37 35166.45 34708.61
## [177] 34960.49 35127.72 34876.68 34458.63 34166.72 34250.10 34167.07 34293.16
## [185] 34503.33 34839.33 35257.94 35259.37 35345.62 35850.02 36439.44 37360.91
## [193] 37361.36 36443.14 35691.05 34899.51 35148.29 34231.54 34690.87 34565.06
## [201] 34189.09 34023.78 34610.71 35362.34 35068.59 34150.45 34995.43 35853.23
## [209] 34652.46 35619.60
plot(zoo(Insample, order.by = Ayam$Minggu),main="Data Insample",col="green",lwd=3)
lines(zoo(Pred_in, order.by = Ayam$Minggu),col = "blue", lwd = 3, lty = 2)
legend("topleft", legend = c("Insample", "Pred Insample"),col = c("green", "blue"), lty = c(1, 2), lwd = 3)

#Prediksi Out Sample
Outsample = Ayam$Harga[211:263]
tnew =c(211:263)
D2new =Ayam$D2[211:263]
D3new =Ayam$D3[211:263]
newxreg=data.frame(tnew,D2new,D3new)
Pred_out = predict(modelarimax,NewInsample=Outsample,newxreg=newxreg,n.ahead=53)$pred
se.Pred_out=predict(modelarimax,newxreg=newxreg,n.ahead=53)$se
Pred_out
## Time Series:
## Start = 211
## End = 263
## Frequency = 1
## [1] 36289.32 36282.40 36278.63 36277.47 36278.49 36281.32 36285.66 36291.26
## [9] 36297.90 36305.42 36313.66 36322.51 36331.86 36341.63 36351.75 36362.16
## [17] 36372.81 36383.66 36394.68 36405.84 36417.12 36428.49 36439.95 36451.47
## [25] 36463.05 36474.67 36486.34 36498.03 36509.76 36521.50 36533.26 36545.04
## [33] 36556.84 36568.64 36580.45 36592.27 36604.10 36615.93 36627.76 36639.60
## [41] 36651.44 36663.29 36675.13 36686.98 36698.83 36710.68 36722.53 37619.94
## [49] 38101.07 36758.09 36769.95 36781.80 36793.66
#Forecast
Ayam.fore = rep(0,12)
t.fore = c(264:275)
D2.fore= rep(0,12)
D3.fore= rep(0,12)
xreg.fore= data.frame (t.fore,D2.fore,D3.fore)
Forecastt=predict(modelarimax,NewInsample=Ayam.fore,newxreg=xreg.fore,12)$pred
se.fore.ARIMAA=predict(modelarimax,NewInsample=Ayam.fore,newxreg=xreg.fore,12)$se
Forecastt
## Time Series:
## Start = 211
## End = 222
## Frequency = 1
## [1] 36917.72 36910.80 36907.03 36905.87 36906.89 36909.72 36914.06 36919.66
## [9] 36926.30 36933.82 36942.06 36950.91
se_foree=se.fore.ARIMAA
data.frame(Forecastt)
## Forecastt
## 1 36917.72
## 2 36910.80
## 3 36907.03
## 4 36905.87
## 5 36906.89
## 6 36909.72
## 7 36914.06
## 8 36919.66
## 9 36926.30
## 10 36933.82
## 11 36942.06
## 12 36950.91
#Akurasi
smape(Ayam$Harga,c(Pred_in,Pred_out))
## [1] 0.02498434
mape(Ayam$Harga,c(Pred_in,Pred_out))
## [1] 0.02491483
rmse(Ayam$Harga,c(Pred_in,Pred_out))
## [1] 1397.031
#Plot
{plot(zoo(Zt, order.by = Ayam$Minggu),
main="Data Harga Ayam",col="black",lwd=2,
ylab = "Harga Ayam dalam Rupiah",
xlab = "Periode",
xlim = as.Date(c("2019-05-06", "2024-08-05")))
lines(zoo(Pred_in, order.by = Ayam$Minggu[1:210]),col = "blue", lwd = 3, lty = 2)
lines(zoo(Pred_out, order.by = Ayam$Minggu[211:263]),col = "green", lwd = 3, lty = 2)
lines(zoo(Forecastt, order.by = seq(as.Date("2024-05-20"), by = "week", length.out = 12)),col = "orange", lwd = 3, lty = 6)
abline(v = Ayam$Minggu[210], col = "red", lty = 3, lwd = 2)
text(x = Ayam$Minggu[210], y = 42000, labels = "Batas Insample",srt = 90, pos = 4, col = "red", cex = 0.8)
abline(v = Ayam$Minggu[263], col = "red", lty = 3, lwd = 2)
text(x = Ayam$Minggu[263], y = 42000, labels = "Batas Outsample",srt = 90, pos = 4, col = "red", cex = 0.8)
legend("topleft", lwd = 3,bty = "n",
legend = c("Data Ayam", "Pred Insample","Pred Outsample","Forecast"),
col = c("black","blue", "green","orange"), lty = c(1,2,2,6))}

{plot(zoo(Zt, order.by = Ayam$Minggu),
main="Data Harga Ayam",col="black",lwd=2,
ylab = "Harga Ayam dalam Rupiah",
xlab = "Periode",
xlim = as.Date(c("2022-05-06", "2024-08-05")),
cex = 3)
lines(zoo(Pred_in, order.by = Ayam$Minggu[1:210]) ,col = "blue", lwd = 3, lty = 2)
lines(zoo(Pred_out, order.by = Ayam$Minggu[211:263]) ,col = "green", lwd = 3, lty = 2)
lines(zoo(Forecastt, order.by = as.Date(c("2024-05-20", "2024-08-05"))) ,col = "orange", lwd = 3, lty = 6)
abline(v = Ayam$Minggu[210], col = "red" , lty = 3, lwd = 2)
text(x = Ayam$Minggu[210], y = 42000, labels = "Batas Insample",srt = 90, pos = 4, col = "red", cex = 0.8)
abline(v = Ayam$Minggu[263] , col = "red", lty = 3, lwd = 2)
text(x = Ayam$Minggu[263], y = 42000, labels = "Batas Outsample",srt = 90, pos = 4, col = "red", cex = 0.8)
legend("bottomleft" , lwd = 3,bty = "n",
legend = c("Data Ayam", "Pred Insample","Pred Outsample","Forecast"),
col = c("black","blue", "green","orange"), lty = c(1,2,2,6))
}
