#==========================@Install Packages@=======================#
library(zoo)
library(tseries)
library(forecast)
library(lmtest)
library(aTSA)
library(xts)
library(fGarch)
library(FinTS)
library(rugarch)
library(MASS)

#==========================@Setting Working Directory@=======================#
getwd()
setwd("C:/PAPER RIZAL")
ls()
rm(list=ls(all=TRUE))

#===============================@Read Data@==================================#
data<-read.csv("Indofood_2020-2023.csv", sep=";", header=T)
data

#==============================@Memuat Data ts@==============================#
ts_data<-ts(data$Close)
ts_data
str(ts_data)

#Membagi data train dan test
data_Indofood.ts.train <- ts_data[1:778]
data_Indofood.ts.train
data_Indofood.ts.test <- ts_data[779:974]
data_Indofood.ts.test

#=====================@Pemeriksaan Kestasioneran Data@=======================#
#==============@pemeriksaan kestasioneran melalui plot data@=================#
#uji varians
lambda=BoxCox.lambda(data_Indofood.ts.train)
lambda
library(TSA)
plot(data_Indofood.ts.train,type="l",xlab="Data ke",ylab="Harga Penutupan Saham Indofood (Rp)",main="Plot Data Harian Harga Penutupan Saham Indofood")
###berdasarkan plot, terlihat data sangat berfluktuasi dari waktu ke waktu yang mengindikasikan data tidak stasioner

#=====@pemeriksaan kestasioneran melalui grafik correlogram ACF dan PACF@====#
par(mfrow=c(2,1))
acf(data_Indofood.ts.train,lag.max=50)
pacf(data_Indofood.ts.train,lag.max=50)
###berdasarkan plot ACF dan PACF data, banyak lag yang keluar dari garis bartlet yang mengindikasikan data tidak stasioner
#=============@pemeriksaan kestasioneran melalui uji ADF@====================#
library(tseries)
adf.test(data_Indofood.ts.train)
###berdasarkan adf tes, p-value>0.05 yang menunjukkan data tidak stasioner

#==========================@Differencing Data@===============================#
diff=diff(data_Indofood.ts.train)
diff

#===========@Pemeriksaan kestasioneran data differencing@=====================#
#==============@pemeriksaan kestasioneran melalui plot data@=================#
plot(diff,type="l",xlab="Data ke",ylab="Harga Penutupan Saham Indofood (Rp)",main="Plot Data Differencing Harian Harga Penutupan Saham Indofood",col = "blue")
###berdasarkan plot, terlihat data bergerak dikisaran garis tengah yang mengindikasikan data sudah stasionerr
#=====@pemeriksaan kestasioneran melalui grafik correlogram ACF dan PACF@====#
par(mfrow=c(2,1))
acf(diff,lag.max=50)
pacf(diff,lag.max=50)
###berdasarkan plot ACF dan PACF data, hanya beberapa lag yang keluar dari garis bartlet yang mengindikasikan data  stasioner
#=============@pemeriksaan kestasioneran melalui uji ADF@====================#
adf.test(diff)
###p-value=0.01<alpha(0.05), sehingga data differencing sudah stasioner dan dpt dilanjutkan ke pemodelan

#==========================@Pemodelan ARIMA@=================================#
library(forecast)
arima = auto.arima(data_Indofood.ts.train,trace=T)
summary(arima) 

model1 <- Arima(diff,order = c(0,1,2),method ="ML")
summary(model1)
model2 <- Arima(diff,order = c(0,1,9),method ="ML")
summary(model2)
model3 <- Arima(diff,order = c(2,1,0),method ="ML")
summary(model3)

#======================@Pemilihan Model Terbaik@=============================#
library(kableExtra)
Model <- c("ARIMA (0,1,2)","ARIMA (0,1,9)","ARIMA (2,1,0)")
AIC <- c(model1$aic,model2$aic,model3$aic)
BIC <- c(model1$bic,model2$bic,model3$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)
paste("Model yang terbaik berdasarkan akurasi adalah model",Akurasi$Model[which.min(Akurasi[,"AIC"])])

#======================@Overfitting@=============================#
model4 <- Arima(diff, order=c(1,1,2),method="ML") #intersep signifikan
lmtest::coeftest(model4)
Model <- c("ARIMA (1,1,2)")
AIC <- c(model4$aic)
BIC <- c(model4$bic)
Akurasi <- data.frame(Model,AIC,BIC)
Akurasi

model5 <- Arima(diff, order=c(0,1,3),method="ML") #intersep signifikan
lmtest::coeftest(model5)
Model <- c("ARIMA (0,1,3)") #Model terbaik
AIC <- c(model5$aic)
BIC <- c(model5$bic)
Akurasi <- data.frame(Model,AIC,BIC)
Akurasi

#UJI DIAGNOSTIK
#==========================@Analisis Sisaan@================================#
sisaan <- model5$residuals
tsdiag(model5)
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(sisaan, type="o",ylab = "Error", xlab = "Order", main = "Error Model 5 vs Order")
abline(h = 0, col='red')
acf(sisaan, main="Series Error")
pacf(sisaan, main="Series Error")

#=======================@Sisaan Menyebar Normal@============================#
jarque.bera.test(sisaan)
#Berdasarkan Jarque Bera Test, diperoleh p-value < α (0.05), maka Tolak H0. 
#Artinya, cukup bukti untuk menyatakan bahwa sisaan tidak menyebar normal pada taraf nyata 5%.

#==========@Sisaan saling bebas/tidak ada autokorelasi@=====================#
Box.test(sisaan, type = "Ljung") 
#Hasil : p−value=0.941>α=0.05 yang berarti TERIMA H0, Tidak terdapat gejala autokorelasi pada sisaan

#=================@Nilai tengah sisaan sama dengan nol@=====================#
t.test(sisaan, mu = 0, conf.level = 0.95)
#Hasil : p−value=0.3227>α=0.05 yang berarti Terima H0, Nilai tengah sisaan sama dengan 0
#MAPE ARIMA


# Anggaplah data_Indofood.ts.train dan data_Indofood.ts.test adalah objek deret waktu

# Latih model ARIMA
model <- Arima(data_Indofood.ts.train, order = c(0,1,3), method = "ML", include.drift = TRUE)

# Hitung ramalan menggunakan model yang dilatih

ramalan <- forecast::forecast(model, h = length(data_Indofood.ts.test))
data.ramalan <- as.numeric(ramalan$mean)  # Konversi ramalan menjadi vektor numerik

# Hitung kesalahan absolut
error <- data_Indofood.ts.test - data.ramalan

# Hitung kesalahan relatif
r.error <- (error / data_Indofood.ts.test) * 100

# Hitung Mean Absolute Percentage Error (MAPE)
MAPE.ARIMA <- mean(abs(r.error), na.rm = TRUE)
MAPE.ARIMA

#=========================@Identifikasi Efek ARCH@==========================#
library(FinTS)
bydArchTest1 <- ArchTest(sisaan, lags=1, demean=TRUE)
bydArchTest1

bydArchTest2 <- ArchTest(sisaan, lags=2, demean=TRUE)
bydArchTest2

bydArchTest3 <- ArchTest(sisaan, lags=3, demean=TRUE)
bydArchTest3

bydArchTest4 <- ArchTest(sisaan, lags=4, demean=TRUE)
bydArchTest4

bydArchTest5 <- ArchTest(sisaan, lags=5, demean=TRUE)
bydArchTest5

bydArchTest6 <- ArchTest(sisaan, lags=6, demean=TRUE)
bydArchTest6

#setelah dilakukan efek ARCH untuk 5 lag awal, semua p-value <5% menunjukkan signifikan
#sehingga menggunakan GARCH


for (i in 1:15) {
  ArchTest <- ArchTest(model2$residuals, lags=i, demean=TRUE)
  cat("P Value LM Test lag ke", i,"adalah" , ArchTest$p.value, "\n") }
#diperoleh nilai P−Value sampai lag ke-15 yang diperoleh < dari α = 5%, maka tolak H0
#artinya terdapat unsur heteroskedastisitas pada sisaan model ARIMA (1,1,3)

#Pemodelan GARCH
library(fGarch) 
garch10<-garchFit(~arma(0,1,3)+garch(1,0),data = diff, trace = F)
summary(garch10) 

garch11<-garchFit(~arma(0,1,3)+garch(1,1),data = diff, trace = F)
summary(garch11) 

garch13<-garchFit(~arma(0,1,3)+garch(1,3),data = diff, trace = F) #AIC TERKECIL
summary(garch13) 

garch23<-garchFit(~arma(0,1,3)+garch(2,3),data = diff, trace = F)
summary(garch23)

garch33<-garchFit(~arma(0,1,3)+garch(3,3),data = diff, trace = F)
summary(garch33) 

#AIC terkecil terdapat pada model arima(0,1,3)+garch(1,3)

#Peramalan GARCH
garchSpec <- ugarchspec(
  variance.model=list(model="sGARCH",garchOrder=c(1,3)),
  mean.model=list(armaOrder=c(0,1,3)), distribution.model="std")
garchFitt <- ugarchfit(spec=garchSpec, data=diff)

forc<- ugarchforecast(fitORspec = garchFitt, data = data$Close, n.ahead = 20, n.roll = 0)
forc
plot(forc, which= 1) 

# Perbandingan aktual dan ramalan
pt_1 <- data$Close[778] #nilai akhir data latih
hasil.forc.Diff <- forc@forecast$seriesFor[,1]
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1
length(data_Indofood.ts.test) <- length(hasil[-1])
perbandingan <- data.frame("Aktual"= data_Indofood.ts.test, "Ramalan" = hasil[-1])
perbandingan

#head(perbandingan,10)

#MAPE GARCH
Date <- c(data$Date [1998:2260])
length(Date) <- length(perbandingan)
dataframe <- data.frame(Date, perbandingan) 
T <- nrow(dataframe)

MAPE.GARCH <- 1/T*sum(abs((dataframe$Aktual-dataframe$Ramalan)/dataframe$Aktual)*100)
MAPE.GARCH

#Perbandingan MAPE ARIMA dan ARIMA-GARCH
MAPE <- data.frame("MAPE ARIMA"= MAPE.ARIMA, "MAPE GARCH" = MAPE.GARCH)
MAPE