1. Pendahuluan

Salah satu metode yang sering diterapkan dalam analisis runtun waktu yaitu metode ARIMA (Autoregressive Integrated Moving Average). Metode ini mampu mengidentifikasi pola tren dan pergerakan data berdasarkan nilai-nilai terdahulu, sehingga sering digunakan dalam perkiraan di berbagai sektor, termasuk industri energi dan minyak bumi. Dengan memanfaatkan data produksi minyak dari sumur yang telah ada, metode ARIMA diharapkan dapat menghasilkan estimasi produksi yang tepat dan dapat dijadikan acuan dalam perencanaan operasional.

2. Import Library

Library yang diperlukan pada analisis ini adalah sebagai berikut:
readxl : digunakan untuk membaca file Excel (.xls/.xlsx) ke dalam R tanpa perlu membuka Excel.
dplyr : digunakan untuk manipulasi dan pengolahan data seperti menyaring data, menambah atau mengubah variabel, memilih kolom, mengurutkan data, dll.
lubridate : digunakan untuk mengelola data tanggal dan waktu.
forecast : digunakan untuk analisis dan peramalan time series termasuk ARIMA, Box-Cox transformation, forecasting, plot ACF dan PACF.
tseries : digunakan untuk uji statistik time series seperti ADF, Jarque-Bera, dll.
KFAS : digunakan untuk model trate space dan Kalman Filter.
ggplot2 : digunakan untuk visualisasi data yang lebih fleksibel dan estetik.
FinTS : digunakan untuk analisis time series finansial, terutama uji ARCH dan deteksi volatilitas.

library(readxl)
library(dplyr)
library(lubridate)
library(forecast)
library(tseries)
library(KFAS)
library(ggplot2)
library(FinTS)

3. Time Series Plot Analysis

data_raw <- read_excel("D:/KULIAH/UCM/LSE-1174 UCM.xlsx")

colnames(data_raw) <- as.character(unlist(data_raw[2, ]))
data_raw <- data_raw[-c(1, 2, 3), ]
data_raw <- data_raw[-c(1), ]

data <- data_raw %>%
  as_tibble() %>%
  rename(
    Date = Tanggal,
    Oil  = `Oil Display`
  ) %>%
  mutate(
    Date = as.numeric(Date),
    Date = as.Date(Date, origin = "1899-12-30"),
    Oil  = as.numeric(Oil)
  ) %>%
  arrange(Date) %>%
  filter(!is.na(Date), !is.na(Oil))

data_ts <- data %>%
  select(Date, Oil) %>%
  filter(!is.na(Oil))

ggplot(data_ts, aes(x = Date, y = Oil)) +
  geom_line() +
  labs(
    title = "Time Series Plot: Oil",
    x = "Tanggal",
    y = "Oil"
  ) +
  theme_minimal()

Berdasarkan grafik time series, produksi minyak harian pada sumur LDC tersebut menunjukkan tren menurun secara umum. Pada awal periode pengamatan, laju produksi minyak berada pada kisaran yang relatif tinggi, kemudian mengalami penurunan cukup tajam pada fase awal, diikuti oleh fluktuasi dengan kecenderungan menurun yang lebih gradual pada periode selanjutnya.

Selama periode pengamatan, data produksi menunjukkan adanya fluktuasi jangka pendek di sekitar tren utama, yang mengindikasikan variasi produksi harian akibat faktor operasional maupun kondisi sumur. Penurunan produksi terlihat semakin signifikan mendekati bulan September hingga Oktober, sebelum kemudian cenderung stabil pada level yang lebih rendah hingga akhir periode pengamatan. Secara keseluruhan, pola data mencerminkan karakteristik deret waktu produksi sumur minyak yang tidak stasioner, dengan adanya tren dan fluktuasi, sehingga sesuai untuk dianalisis menggunakan metode peramalan berbasis deret waktu seperti Autoregressive Integrated Moving Average (ARIMA).

4. Cek Stasioneritas

4.1 Stasioneritas dalam varian

library(forecast)
lambda <- BoxCox.lambda(
  data_ts$Oil,
  method = "guerrero"
)
lambda
## [1] -0.06952849

Berdasarkan hasil analisis Box-Cox, diperoleh nilai parameter transformasi sebesar λ= -0.07. Karena λ ≠0, maka perlu dilakukan transformasi, dimana transformasi yang digunakan adalah Box-Cox.

oil_var <- BoxCox(data_ts$Oil, lambda)
plot(
  data_ts$Date,
  oil_var,
  type = "l",
  xlab = "Tanggal",
  ylab = "Oil (Box-Cox Transformed)",
  main = paste("Produksi Setelah Transformasi Box-Cox")
)

lambda_trans <- BoxCox.lambda(oil_var, method = "guerrero")
lambda_trans
## [1] 1.187021

4.2 Stasioneritas dalam mean

Hipotesis
H0 : Data tidak stasioner
H1 : Data stasioner

adf_result <- adf.test(oil_var)
adf_result
## 
##  Augmented Dickey-Fuller Test
## 
## data:  oil_var
## Dickey-Fuller = -2.3999, Lag order = 6, p-value = 0.4079
## alternative hypothesis: stationary
adf_result$p.value
## [1] 0.4078697

Berdasarkan hasil uji ADF terhadap data produksi sumur LDC, diperoleh nilai Dickey-Fuller yaitu -2.3999 dengan p-value sebesar 0.4079. Karena nilai p-value lebih besar dari alpha 5%, maka hipotesis nol gagal ditolak. Dengan demikian, data produksi minyak harian Sumur LSE-1174 dinyatakan belum stasioner dalam mean. Oleh karena itu, perlu dilakukan differencing ADF test sebanyak satu kali atau sampai nilai p-value berada dibawah taraf signifikansi.

# Differencing pertama
oil_diff1 <- diff(oil_var, differences = 1)

# Plot hasil differencing
plot(
  data_ts$Date[-1],   # karena 1 data hilang setelah diff
  oil_diff1,
  type = "l",
  xlab = "Tanggal",
  ylab = "Oil (1st Difference)",
  main = "Differencing Orde 1"
)

library(tseries)

adf_diff1 <- adf.test(oil_diff1)
## Warning in adf.test(oil_diff1): p-value smaller than printed p-value
adf_diff1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  oil_diff1
## Dickey-Fuller = -5.8084, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

Setelah dilakukan differencing, didapatkan p-value sebesar 0.01 dimana kurang dari alpha 5% sehingga stasioneritas dalam mean terpenuhi setelah dilakukan differencing satu kali.

5. Identifikasi Orde Model

5.1 Plot ACF dan PACF

par(mfrow = c(1,2))  # supaya 2 plot berdampingan

acf(oil_diff1,
    main = "ACF - Differencing Orde 1")

pacf(oil_diff1,
     main = "PACF - Differencing Orde 1")

par(mfrow = c(1,1))  # reset layout

Berdasarkan plot ACF yang terbentuk, terlihat bahwa hanya lag awal yang menunjukkan nilai autokorelasi yang signifikan, sementara sebagian besar lag selanjutnya berada di dalam batas kepercayaan. Pola ini mengindikasikan bahwa korelasi antar pengamatan menurun dengan cepat dan tidak menunjukkan adanya ketergantungan jangka panjang yang kuat. Kondisi tersebut mengisyaratkan bahwa komponen moving average (MA) yang digunakan dalam pemodelan cenderung berorde rendah.

Sementara itu, plot PACF menunjukkan adanya beberapa spike yang signifikan pada lag awal, kemudian nilai partial autocorrelation berfluktuasi sekitar nol dan berada dalam batas kepercayaan pada lag-lag berikutnya. Pola ini mengindikasikan adanya ketergantungan terhadap nilai masa lalu pada lag tertentu, sehingga mengisyaratkan keberadaan komponen autoregressive (AR) dalam model. Berdasarkan analisis tersebut, maka diperoleh calon model terbaik yaitu ARIMA(0,1,0) dan ARIMA(1,1,0).

5.2 Estimasi dan Pemilihan Model Terbaik

model_110 <- Arima(oil_var, order = c(1,1,0))
model_010 <- Arima(oil_var, order = c(0,1,0))

AIC(model_110, model_010)
##           df       AIC
## model_110  2 -1108.308
## model_010  1 -1110.186

Hasil estimasi menunjukkan bahwa model ARIMA (1,1,0) memiliki nilai AIC yang paling kecil dibandingkan model kandidat lainnya. Oleh karena itu, model ARIMA(1,1,0) dipilih sebagai model terbaik untuk memodelkan dan meramalkan produksi minyak harian Sumur LDC.

6. Uji Diagnostik Residual

6.1 Plot Residual:

resid_110 <- residuals(model_110)

plot(resid_110,
     type = "l",
     main = "Residual ARIMA(1,1,0)",
     ylab = "Residual",
     xlab = "Waktu")

abline(h = 0, col = "red")

6.2 Uji Ljung-Box (White Noise Test)

Hipotesis:
H0 : Residual white noise
H1 : Residual tidak white noise (model kurang baik)

Box.test(resid_110,
         lag = 20,
         type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  resid_110
## X-squared = 23.095, df = 20, p-value = 0.2842

Berdasarkan hasil uji Ljung-Box pada lag ke-20, diperoleh nilai p-value sebesar 0.2842. Karena nilai p-value lebih besar dari 0.05, maka hipotesis nol tidak ditolak, sehingga residual model dianggap bersifat white noise. Dengan demikian, model ARIMA(1,1,0) telah memenuhi asumsi diagnostik dan layak digunakan untuk peramalan.

6.3 Uji Normalitas Residual

Hipotesis:
H0 : Residual berdistribusi normal
H1 : Residual tidak berdistribusi normal

resid_110 <- residuals(model_110)

shapiro.test(resid_110)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_110
## W = 0.77044, p-value < 2.2e-16

Berdasarkan uji Kolmogorov-Smirnov (KS), diketahui bahwa residual belum berdistribusi normal karena p-value kurang dari alpha 5% sehingga H0 ditolak. Namun, ketidaknormalan residual ini untuk sementara diabaikan sebagaimana pada penelitian yang dilakukan oleh Wulansari et al., (2014).

6.3 Uji Homoskedastisitas

Hipotesis:
H0 : Tidak ada efek ARCH (homoskedastisitas terpenuhi)
H1 : Ada efek ARCH (homoskedastisitas tidak terpenuhi)

ArchTest(resid_110, lags = 20)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  resid_110
## Chi-squared = 26.343, df = 20, p-value = 0.1548

Berdasarkan uji ARCH pada lag ke-20, didapatkan p-value sebesar 0.1548 dimana lebih besar dari alpha 5% sehingga H0 gagal ditolak. Maka residual tidak memiliki efek ARCH atau dengan kata lain residual homoskedastisitas.

7. Hasil Peramalan

7.1 Plot Hasil Peramalan

df <- data_ts[!is.na(data_ts$Date) & !is.na(data_ts$Oil), c("Date","Oil")]
df <- df[order(df$Date), ]

y_bc <- oil_var
fit_110 <- Arima(y_bc, order = c(1,1,0))

h <- 7
fc_110 <- forecast(fit_110, h = h) #forecast skala Box-Cox

#Mengembalikan ke skala asli
if (!exists("shift")) shift <- 0

fitted_orig   <- InvBoxCox(fitted(fit_110), lambda) - shift
forecast_orig <- InvBoxCox(as.numeric(fc_110$mean), lambda) - shift

fitted_plot <- rep(NA_real_, nrow(df))
idx_fit <- (length(fitted_plot) - length(fitted_orig) + 1):length(fitted_plot)
fitted_plot[idx_fit] <- fitted_orig

fitted_plot <- rep(NA_real_, nrow(df))
idx_fit <- (length(fitted_plot) - length(fitted_orig) + 1):length(fitted_plot)
fitted_plot[idx_fit] <- fitted_orig

last_date <- max(df$Date, na.rm = TRUE)
future_dates <- seq.Date(from = last_date + 1, by = "day", length.out = h)

# Garis batas mulai forecast (vertical line)
cut_date <- last_date  # garis di tanggal terakhir data aktual

# Atur rentang X supaya muat 7 hari forecast
x_all <- c(df$Date, future_dates)

plot(df$Date, df$Oil,
     type = "l",
     col = "red",
     lwd = 2,
     xlim = range(c(df$Date, future_dates)),
     xlab = "Tanggal",
     ylab = "Oil (BOPD)",
     main = "Perbandingan Produksi Aktual dan Hasil Peramalan (ARIMA 1,1,0)")

points(df$Date, df$Oil, pch = 16, col = "red", cex = 0.5)

lines(df$Date, fitted_plot, col = "blue", lwd = 2)
lines(future_dates, forecast_orig, col = "blue", lwd = 2, lty = 2)

abline(v = max(df$Date), lty = 2, lwd = 2, col = "black")

legend("topright",
       legend = c("AKTUAL",
                  "ARIMA (fitted)",
                  "ARIMA (forecast)"),
       col    = c("red", "blue", "blue"),
       lty    = c(1, 1, 2),
       lwd    = c(2, 2, 2),
       pch    = c(16, NA, NA),
       pt.cex = 0.7,     # ukuran titik lebih kecil
       cex    = 0.8,     # ukuran teks lebih kecil
       y.intersp = 0.7,  # jarak antar baris diperkecil
       x.intersp = 0.5,  # jarak simbol ke teks lebih rapat
       inset  = c(0.01, 0.01),  # makin ke pojok kanan atas
       bty    = "n")

7.2 Tabel Hasil Peramalan

forecast_table <- data.frame(
  Tanggal = future_dates,
  `Hasil Forecast` = round(forecast_orig, 3)
)

forecast_table
##      Tanggal Hasil.Forecast
## 1 2026-01-12         117.61
## 2 2026-01-13         117.61
## 3 2026-01-14         117.61
## 4 2026-01-15         117.61
## 5 2026-01-16         117.61
## 6 2026-01-17         117.61
## 7 2026-01-18         117.61