####-----------------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))
}