{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE)

{r setup, include=FALSE} ## Global options knitr::opts_chunk$set(cache = TRUE)

library(readxl)
library(forecast)
library(TTR)
library(imputeTS)
library(tseries)
library(ggplot2)
library(dplyr)
library(graphics)
library(TSA)
library(tidyverse)
library(lubridate)
library(gridExtra)
library(ggfortify)
library(cowplot)
library(graphics)
library(lmtest)
library(stats)
library(MASS)
library(fpp2)
library(hrbrthemes)

setwd("D:/SATRIA DATA")
mpdw <- read_xlsx("data pm.xlsx")
mpdw <- mpdw[,1:2]
mpdw$`Data PM 2.5`<- as.numeric(mpdw$`Data PM 2.5`)
ts.mpdw <- ts(mpdw)
kableExtra::kable(head(mpdw) ,caption = 'Subset Data CO DKI Jakarta 2021-2022')
view(mpdw)

# Train-test Split (80/20)
train <- mpdw[1:363,]
test <- mpdw[364:454,]

datampdw <- ts(mpdw$`Data PM 2.5`)
ts.plot(datampdw, main = "Kadar PM 2.5 DKI Jakarta Tahun 2021 - 2022", ylab = "Harga", lwd = 1.5)
points(datampdw)
legend("topleft", c("Premium"), cex = 0.7,
       col = c("black", "red", "blue"), lty = 1)

plot(x = mpdw$t,
     y = mpdw$`Data PM 2.5`,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Kadar PM 2.5 DKI Jakarta",
     main = "Time Series Plot Kadar PM 2.5 DKI Jakarta")

# Time-Series Object
data.ts <- ts(mpdw[,2], start = 2021, end = 2022, frequency = 91)
train.ts <- ts(train)
test.ts <- ts(test)

acf(mpdw$`Data PM 2.5`, lag = 35, start= 1)

pacf(mpdw$`Data PM 2.5`, lag =35)

adf.test(data.ts)

Data Stasioner

Differencing

data.diff <- diff(data.ts, differences = 1)

plot(data.diff)

adf.test(data.diff)

Stasioner

acf(data.diff, lag.max = 50, start = 1)

pacf(data.diff, lag.max =50)

eacf(data.diff)

ARIMA (0,1,1) ARIMA (1,1,1) ARIMA (3,1,3) ARIMA (5,1,5)

Perbandingan Kebaikan Model Tentatif

model1 <- Arima(data.ts,order = c(0,1,1),method ="ML")
model2 <- Arima(data.ts,order = c(1,1,1),method ="ML")
model3 <- Arima(data.ts,order = c(3,1,3),method ="ML")
model4 <- Arima(data.ts,order = c(5,1,5),method = "ML")
summary(model1)
summary(model2)
summary(model3)
summary(model4)
coeftest(model1)
coeftest(model2)
coeftest(model3)

Model <- c("ARIMA (0,1,1)","ARIMA (1,1,1)","ARIMA (3,1,3)","ARIMA (5,1,5)")
AIC <- c(model1$aic,model2$aic,model3$aic,model4$aic)
BIC <- c(model1$bic,model2$bic,model3$bic,model4$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)

model5 <- auto.arima(data.t)
summary(model5)
coeftest(model6)

paste("Model yang terbaik adalah model",Akurasi$Model[which.min(Akurasi[,"AIC"])])

coeftest(model1)
summary(model1)

over1 <- Arima(train.ts[,2], order = c(0,1,2))
summary(over1)
coeftest(over1)

over2 <- Arima(train.ts[,2], order= c(2,1,0))
summary(over2)
coeftest(over2)

over3 <- Arima(train.ts[,2], order = c(1,1,2))
summary(over3)

over4 <- Arima(train.ts[,2], order = c(2,1,1))
summary(over4)


# Uji Diagnostik
## Diagnostik Model: Eksploratif

```{r}
sisaan <- over3$residuals
par(mfrow=c(2,2))
par(bg = '#141415')
qqnorm(sisaan)
box(col="white",lwd=2)
qqline(sisaan, col = "red", lwd =1, col.lab="white",
       col.axis="white",col.sub = "white")
box(col="white",lwd=2)
plot(c(1:length(sisaan)),sisaan,col="white",col.lab="white",col.axis="white")
box(col="white",lwd=2)
acf(sisaan,col="white",col.sub = "white",col.axis="white", col.lab="white")
box(col="white",lwd=2)
pacf(sisaan,col="white",col.sub = "white",col.axis="white", col.lab="white",col.main="white")
box(col="white",lwd=2)

Berdasarkan hasil eksplorasi menggunakan Q-Q plot, terlihat bahwa sisaan berdistribusi mengikuti garis normal, sehingga dapat dikatakan bahwa sisaan menyebar normal. Kemudian, plot sisaan yang diperoleh menunjukkan bahwa sisaan memiliki pola acak dan tersebar di sekitar nilai nol. Sedangkan pada plot ACF dan PACF, nilai awal amatan tidak melewati garis signifikan, atau dapat dikatakan bahwa sisaan saling bebas.

Diagnostik Model: Uji Formal

1. Sisaan Menyebar Normal

Uji formal ini dilakukan dengan Jarque Bera test dan Shapiro-Wilk test.

Hipotesis yang diuji: H0 : Sisaan Menyebar Normal H1: Sisaan Tidak Menyebar Normal

jarque.bera.test(sisaan)

shapiro.test(sisaan)

2. Sisaan Saling Bebas

Uji formal ini dilakukan dengan LJung-Box test.

Hipotesis yang diuji: H0 : Sisaan antara lag saling bebas H1: Sisaan antara lag tidak saling bebas

Box.test(sisaan, type = "Ljung")

Berdasarkan LJung-Box test, diperoleh p-value (0.909) > α (0.05), maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag saling bebas atau dapat dikatakan tidak ada autokorelasi antara sisaan lag pada taraf nyata 5%.

3. Nilai Tengah Sisaan Sama dengan Nol

Uji formal ini dilakukan dengan t-test.

Hipotesis yang diuji: H0 : Nilai tengah sisaan sama dengan nol H1: Nilai tengah sisaan tidak sama dengan nol

t.test(sisaan, mu = 0, conf.level = 0.95)

FORECASTING

train.ts <- ts(train$`Data PM 2.5`,start = 1)

ramalan <- forecast::forecast(Arima(train.ts, order=c(1,1,2),method="ML",include.drift = TRUE), h=10) 
data.ramalan <- ramalan$mean
plot(ramalan,lwd=2)

accuracy(ramalan)

accuracy(test.ts,data.ramalan)

accuracy(ramalan,test = test.ts)