# ============================================================================
# ANALISIS DAN PROYEKSI EMISI CO2 INDONESIA (1970–2024)
# Metode: Regresi Linear, ARIMA, ARCH/GARCH, dan Forecasting
# Sumber Data: Worldometer – Indonesia CO2 Emissions
# Penulis: [Nama Penulis]
# ============================================================================
# Memuat semua library
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tseries)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(rugarch)
## Loading required package: parallel
##
## Attaching package: 'rugarch'
##
## The following object is masked from 'package:purrr':
##
## reduce
##
## The following object is masked from 'package:stats':
##
## sigma
library(FinTS)
##
## Attaching package: 'FinTS'
##
## The following object is masked from 'package:forecast':
##
## Acf
library(scales)
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor
library(patchwork)
# ==========================
# 1. IMPORT & PRAPEMROSESAN DATA
# ==========================
# Membaca data dari file Excel
raw_data <- read_excel("d:/Training/C02_Emisi.xlsx")
# Menampilkan struktur awal data
cat("=== STRUKTUR DATA ===\n")
## === STRUKTUR DATA ===
str(raw_data)
## tibble [55 × 5] (S3: tbl_df/tbl/data.frame)
## $ Year : num [1:55] 2024 2023 2022 2021 2020 ...
## $ Fossil CO2 Emissions (tons) : num [1:55] 8.12e+08 7.73e+08 7.23e+08 6.23e+08 5.96e+08 ...
## $ YoY Change (%) : num [1:55] 5.13 6.9 16.06 4.43 -6.48 ...
## $ CO2 per Capita (tons/person): num [1:55] 2.88 2.76 2.61 2.27 2.19 2.37 2.24 2.03 1.94 1.98 ...
## $ Global Share (%) : num [1:55] 2.05 1.98 1.87 1.63 1.65 1.68 1.57 1.45 1.39 1.41 ...
cat("\n=== RINGKASAN STATISTIK DESKRIPTIF ===\n")
##
## === RINGKASAN STATISTIK DESKRIPTIF ===
summary(raw_data)
## Year Fossil CO2 Emissions (tons) YoY Change (%)
## Min. :1970 Min. : 31123671 Min. :-6.480
## 1st Qu.:1984 1st Qu.: 97023342 1st Qu.: 1.722
## Median :1997 Median :275520703 Median : 5.695
## Mean :1997 Mean :295431061 Mean : 6.434
## 3rd Qu.:2010 3rd Qu.:459323580 3rd Qu.: 8.773
## Max. :2024 Max. :812204159 Max. :36.690
## NA's :1
## CO2 per Capita (tons/person) Global Share (%)
## Min. :0.270 Min. :0.2000
## 1st Qu.:0.630 1st Qu.:0.4950
## Median :1.350 Median :1.1200
## Mean :1.295 Mean :0.9767
## 3rd Qu.:1.855 3rd Qu.:1.3100
## Max. :2.880 Max. :2.0500
##
# Pra-pemrosesan: rename kolom, urutkan berdasarkan tahun (ascending)
df <- raw_data %>%
rename(
year = `Year`,
co2_tons = `Fossil CO2 Emissions (tons)`,
yoy_change = `YoY Change (%)`,
per_capita = `CO2 per Capita (tons/person)`,
global_pct = `Global Share (%)`
) %>%
arrange(year)
# Konversi ke satuan Juta Ton (Mt) untuk keterbacaan
df <- df %>%
mutate(co2_mt = co2_tons / 1e6)
cat("\n=== DATA SETELAH PRAPEMROSESAN (5 baris pertama) ===\n")
##
## === DATA SETELAH PRAPEMROSESAN (5 baris pertama) ===
print(head(df, 10))
## # A tibble: 10 × 6
## year co2_tons yoy_change per_capita global_pct co2_mt
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1970 31123671 NA 0.27 0.2 31.1
## 2 1971 31496130 1.2 0.27 0.2 31.5
## 3 1972 37263373 18.3 0.31 0.23 37.3
## 4 1973 42033316 12.8 0.34 0.24 42.0
## 5 1974 46956608 11.7 0.37 0.27 47.0
## 6 1975 48794745 3.91 0.37 0.28 48.8
## 7 1976 52527133 7.65 0.39 0.29 52.5
## 8 1977 71798713 36.7 0.52 0.38 71.8
## 9 1978 84466597 17.6 0.6 0.43 84.5
## 10 1979 90919613 7.64 0.63 0.45 90.9
# Membuat objek time series (ts)
# Frekuensi = 1 karena data tahunan, dimulai tahun 1970
co2_ts <- ts(df$co2_mt, start = 1970, frequency = 1)
cat("\n=== OBJEK TIME SERIES ===\n")
##
## === OBJEK TIME SERIES ===
print(co2_ts)
## Time Series:
## Start = 1970
## End = 2024
## Frequency = 1
## [1] 31.12367 31.49613 37.26337 42.03332 46.95661 48.79474 52.52713
## [8] 71.79871 84.46660 90.91961 90.59942 95.65866 96.40196 95.42208
## [15] 97.64472 106.03832 117.10361 120.21502 128.66807 134.48484 161.76870
## [22] 174.65927 182.55245 199.85071 212.74749 239.09898 252.00756 275.52070
## [29] 278.24342 295.25001 299.07024 319.56204 324.93886 353.62302 360.15132
## [36] 362.59705 383.41937 401.79343 397.25537 411.34908 443.68204 499.93621
## [43] 504.25725 474.96512 505.07617 512.39726 507.09651 535.74741 596.46923
## [50] 637.56916 596.28393 622.68684 722.69958 772.56179 812.20416
# ============================================================================
# INTERPRETASI: Data emisi CO2 Indonesia mencakup 55 observasi tahunan
# (1970–2024). Emisi meningkat secara signifikan dari ~31 Mt (1970)
# menjadi ~812 Mt (2024), menunjukkan tren naik (upward trend) yang kuat.
# ============================================================================
# ==========================
# 2. REGRESI LINEAR KLASIK (BASELINE)
# ==========================
cat("\n\n", paste(rep("=", 60), collapse = ""), "\n")
##
##
## ============================================================
cat("BAGIAN 2: REGRESI LINEAR KLASIK\n")
## BAGIAN 2: REGRESI LINEAR KLASIK
cat(paste(rep("=", 60), collapse = ""), "\n\n")
## ============================================================
# Model regresi linear: CO2 (Mt) sebagai fungsi dari tren waktu (Year)
# Y = β0 + β1*Year + ε
model_lm <- lm(co2_mt ~ year, data = df)
cat("=== RINGKASAN MODEL REGRESI LINEAR ===\n")
## === RINGKASAN MODEL REGRESI LINEAR ===
summary(model_lm)
##
## Call:
## lm(formula = co2_mt ~ year, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.78 -35.26 -20.09 28.73 165.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.571e+04 8.417e+02 -30.54 <2e-16 ***
## year 1.302e+01 4.215e-01 30.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.62 on 53 degrees of freedom
## Multiple R-squared: 0.9474, Adjusted R-squared: 0.9464
## F-statistic: 954.2 on 1 and 53 DF, p-value: < 2.2e-16
# UJI ASUMSI KLASIK REGRESI
# 2a. Uji Normalitas Residual (Shapiro-Wilk)
cat("\n--- Uji Normalitas Residual (Shapiro-Wilk) ---\n")
##
## --- Uji Normalitas Residual (Shapiro-Wilk) ---
shapiro_test <- shapiro.test(residuals(model_lm))
print(shapiro_test)
##
## Shapiro-Wilk normality test
##
## data: residuals(model_lm)
## W = 0.84656, p-value = 5.41e-06
cat("Interpretasi: ")
## Interpretasi:
if (shapiro_test$p.value > 0.05) {
cat("p-value > 0.05 → Residual berdistribusi normal (H0 tidak ditolak).\n")
} else {
cat("p-value <= 0.05 → Residual TIDAK berdistribusi normal (H0 ditolak).\n")
}
## p-value <= 0.05 → Residual TIDAK berdistribusi normal (H0 ditolak).
# 2b. Uji Heteroskedastisitas (Breusch-Pagan)
cat("\n--- Uji Heteroskedastisitas (Breusch-Pagan) ---\n")
##
## --- Uji Heteroskedastisitas (Breusch-Pagan) ---
bp_test <- bptest(model_lm)
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: model_lm
## BP = 2.5139, df = 1, p-value = 0.1128
cat("Interpretasi: ")
## Interpretasi:
if (bp_test$p.value > 0.05) {
cat("p-value > 0.05 → Tidak ada heteroskedastisitas (homokedastis).\n")
} else {
cat("p-value <= 0.05 → Terdapat heteroskedastisitas → asumsi OLS dilanggar.\n")
}
## p-value > 0.05 → Tidak ada heteroskedastisitas (homokedastis).
# 2c. Uji Autokorelasi (Durbin-Watson)
cat("\n--- Uji Autokorelasi (Durbin-Watson) ---\n")
##
## --- Uji Autokorelasi (Durbin-Watson) ---
dw_test <- dwtest(model_lm)
print(dw_test)
##
## Durbin-Watson test
##
## data: model_lm
## DW = 0.18417, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
cat("Interpretasi: ")
## Interpretasi:
if (dw_test$p.value > 0.05) {
cat("p-value > 0.05 → Tidak ada autokorelasi pada residual.\n")
} else {
cat("p-value <= 0.05 → Terdapat autokorelasi → OLS tidak efisien.\n")
}
## p-value <= 0.05 → Terdapat autokorelasi → OLS tidak efisien.
# ============================================================================
# INTERPRETASI REGRESI:
# Koefisien β1 (slope) menunjukkan peningkatan rata-rata emisi CO2 per tahun.
# Nilai R² menunjukkan proporsi variasi emisi yang dijelaskan oleh tren waktu.
# Pelanggaran asumsi klasik (terutama autokorelasi) menegaskan perlunya
# model time series (ARIMA) yang lebih tepat untuk data runtun waktu.
# ============================================================================
# ==========================
# 3. MODEL TIME SERIES — ARIMA
# ==========================
cat("\n\n", paste(rep("=", 60), collapse = ""), "\n")
##
##
## ============================================================
cat("BAGIAN 3: PEMODELAN ARIMA\n")
## BAGIAN 3: PEMODELAN ARIMA
cat(paste(rep("=", 60), collapse = ""), "\n\n")
## ============================================================
# 3a. UJI STASIONERITAS — Augmented Dickey-Fuller (ADF)
cat("--- Uji ADF pada Data Level (Original) ---\n")
## --- Uji ADF pada Data Level (Original) ---
adf_level <- adf.test(co2_ts, alternative = "stationary")
## Warning in adf.test(co2_ts, alternative = "stationary"): p-value greater than
## printed p-value
print(adf_level)
##
## Augmented Dickey-Fuller Test
##
## data: co2_ts
## Dickey-Fuller = 1.2364, Lag order = 3, p-value = 0.99
## alternative hypothesis: stationary
cat("Interpretasi: ")
## Interpretasi:
if (adf_level$p.value > 0.05) {
cat("p-value > 0.05 → Data TIDAK stasioner pada level.\n")
cat(" → Diperlukan differencing untuk mencapai stasioneritas.\n")
} else {
cat("p-value <= 0.05 → Data sudah stasioner pada level.\n")
}
## p-value > 0.05 → Data TIDAK stasioner pada level.
## → Diperlukan differencing untuk mencapai stasioneritas.
# Differencing orde pertama (d = 1)
co2_diff1 <- diff(co2_ts, differences = 1)
cat("\n--- Uji ADF pada Data Setelah Differencing (d=1) ---\n")
##
## --- Uji ADF pada Data Setelah Differencing (d=1) ---
adf_diff1 <- adf.test(co2_diff1, alternative = "stationary")
print(adf_diff1)
##
## Augmented Dickey-Fuller Test
##
## data: co2_diff1
## Dickey-Fuller = -2.1618, Lag order = 3, p-value = 0.5094
## alternative hypothesis: stationary
cat("Interpretasi: ")
## Interpretasi:
if (adf_diff1$p.value > 0.05) {
cat("p-value > 0.05 → Masih BELUM stasioner setelah differencing 1x.\n")
cat(" → Coba differencing orde 2 (d=2).\n")
co2_diff2 <- diff(co2_ts, differences = 2)
cat("\n--- Uji ADF pada Data Setelah Differencing (d=2) ---\n")
adf_diff2 <- adf.test(co2_diff2, alternative = "stationary")
print(adf_diff2)
} else {
cat("p-value <= 0.05 → Data STASIONER setelah differencing 1x (d=1).\n")
}
## p-value > 0.05 → Masih BELUM stasioner setelah differencing 1x.
## → Coba differencing orde 2 (d=2).
##
## --- Uji ADF pada Data Setelah Differencing (d=2) ---
## Warning in adf.test(co2_diff2, alternative = "stationary"): p-value smaller
## than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: co2_diff2
## Dickey-Fuller = -4.9452, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
# 3b. PLOT ACF DAN PACF
cat("\n--- Membuat Plot ACF dan PACF ---\n")
##
## --- Membuat Plot ACF dan PACF ---
# Plot ACF & PACF dari data yang sudah di-differencing
# Digunakan untuk identifikasi order p (AR) dan q (MA)
par(mfrow = c(1, 2))
acf(co2_diff1, main = "ACF — Emisi CO2 (Differenced, d=1)", lag.max = 20)
pacf(co2_diff1, main = "PACF — Emisi CO2 (Differenced, d=1)", lag.max = 20)
par(mfrow = c(1, 1))
# ============================================================================
# INTERPRETASI ACF/PACF:
# - ACF: Jika korelasi menurun secara gradual → komponen AR dominan.
# Jika terpotong tajam setelah lag-q → komponen MA order q.
# - PACF: Jika terpotong tajam setelah lag-p → komponen AR order p.
# Pola ini menjadi dasar pemilihan order (p, d, q) secara manual.
# ============================================================================
# 3c. ESTIMASI MODEL ARIMA
# Metode 1: auto.arima (seleksi otomatis berdasarkan AICc)
cat("\n--- Model ARIMA Otomatis (auto.arima) ---\n")
##
## --- Model ARIMA Otomatis (auto.arima) ---
model_auto <- auto.arima(
co2_ts,
seasonal = FALSE, # data tahunan, tidak ada pola musiman
stepwise = FALSE, # pencarian menyeluruh (exhaustive search)
approximation = FALSE, # tanpa aproksimasi untuk akurasi
trace = TRUE # tampilkan semua model yang dievaluasi
)
##
## ARIMA(0,2,0) : 498.6663
## ARIMA(0,2,1) : 478.1594
## ARIMA(0,2,2) : 476.1763
## ARIMA(0,2,3) : 471.4411
## ARIMA(0,2,4) : 472.0159
## ARIMA(0,2,5) : 473.4288
## ARIMA(1,2,0) : 498.73
## ARIMA(1,2,1) : 479.4677
## ARIMA(1,2,2) : 476.4028
## ARIMA(1,2,3) : 473.2929
## ARIMA(1,2,4) : 474.108
## ARIMA(2,2,0) : 474.472
## ARIMA(2,2,1) : 468.2624
## ARIMA(2,2,2) : 470.5711
## ARIMA(2,2,3) : Inf
## ARIMA(3,2,0) : 469.8265
## ARIMA(3,2,1) : 470.5541
## ARIMA(3,2,2) : 473.0648
## ARIMA(4,2,0) : 471.763
## ARIMA(4,2,1) : 473.0536
## ARIMA(5,2,0) : 474.227
##
##
##
## Best model: ARIMA(2,2,1)
cat("\n=== MODEL ARIMA TERPILIH (auto.arima) ===\n")
##
## === MODEL ARIMA TERPILIH (auto.arima) ===
summary(model_auto)
## Series: co2_ts
## ARIMA(2,2,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.0215 -0.5317 -0.6742
## s.e. 0.1455 0.1244 0.1602
##
## sigma^2 = 349.7: log likelihood = -229.71
## AIC=467.43 AICc=468.26 BIC=475.31
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3.378189 17.83073 12.69953 1.181027 4.438248 0.7260615 -0.05919132
# Metode 2: Estimasi manual berdasarkan pola ACF/PACF
# Sebagai perbandingan, dicoba beberapa spesifikasi
cat("\n--- Perbandingan Model ARIMA Manual ---\n")
##
## --- Perbandingan Model ARIMA Manual ---
# Kandidat model berdasarkan inspeksi visual ACF/PACF
candidates <- list(
c(1, 1, 0), # AR(1) dengan d=1
c(1, 1, 1), # ARMA(1,1) dengan d=1
c(2, 1, 0), # AR(2) dengan d=1
c(2, 1, 1), # ARMA(2,1) dengan d=1
c(0, 1, 1), # MA(1) dengan d=1
c(0, 1, 2), # MA(2) dengan d=1
c(2, 1, 2), # ARMA(2,2) dengan d=1
c(1, 2, 1) # ARMA(1,1) dengan d=2 (jika d=2 diperlukan)
)
# Evaluasi AIC & BIC untuk setiap kandidat
comparison_results <- data.frame(
Model = character(),
AIC = numeric(),
BIC = numeric(),
stringsAsFactors = FALSE
)
for (params in candidates) {
model_name <- paste0("ARIMA(", params[1], ",", params[2], ",", params[3], ")")
tryCatch({
fit <- Arima(co2_ts, order = params)
comparison_results <- rbind(comparison_results, data.frame(
Model = model_name,
AIC = round(fit$aic, 2),
BIC = round(fit$bic, 2)
))
}, error = function(e) {
cat(" ⚠ Model", model_name, "gagal diestimasi:", conditionMessage(e), "\n")
})
}
# Urutkan berdasarkan AIC (terkecil = terbaik)
comparison_results <- comparison_results %>% arrange(AIC)
cat("\n=== PERBANDINGAN MODEL ARIMA (diurutkan berdasarkan AIC) ===\n")
##
## === PERBANDINGAN MODEL ARIMA (diurutkan berdasarkan AIC) ===
print(comparison_results)
## Model AIC BIC
## 1 ARIMA(1,2,1) 478.98 484.89
## 2 ARIMA(2,1,2) 486.20 496.14
## 3 ARIMA(0,1,1) 491.01 494.99
## 4 ARIMA(1,1,1) 492.91 498.87
## 5 ARIMA(0,1,2) 492.92 498.89
## 6 ARIMA(1,1,0) 493.85 497.83
## 7 ARIMA(2,1,1) 494.75 502.71
## 8 ARIMA(2,1,0) 495.25 501.21
# Pilih model terbaik: model dengan AIC terkecil
best_order <- candidates[[which.min(sapply(candidates, function(p) {
tryCatch(Arima(co2_ts, order = p)$aic, error = function(e) Inf)
}))]]
model_best <- Arima(co2_ts, order = best_order)
cat("\n=== MODEL ARIMA TERBAIK (Manual) ===\n")
##
## === MODEL ARIMA TERBAIK (Manual) ===
cat("Order:", paste0("ARIMA(", best_order[1], ",", best_order[2], ",", best_order[3], ")\n"))
## Order: ARIMA(1,2,1)
summary(model_best)
## Series: co2_ts
## ARIMA(1,2,1)
##
## Coefficients:
## ar1 ma1
## 0.1445 -0.9028
## s.e. 0.1492 0.0661
##
## sigma^2 = 444.8: log likelihood = -236.49
## AIC=478.98 AICc=479.47 BIC=484.89
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3.766276 20.30994 12.76764 1.350541 4.19482 0.7299557 0.01300726
# Memutuskan model final: bandingkan auto.arima vs manual
if (model_auto$aic <= model_best$aic) {
model_arima <- model_auto
cat("\n✔ Model final dipilih: auto.arima →",
paste0("ARIMA(", model_auto$arma[1], ",", model_auto$arma[6], ",", model_auto$arma[2], ")"),
"| AIC =", round(model_auto$aic, 2), "\n")
} else {
model_arima <- model_best
cat("\n✔ Model final dipilih: Manual →",
paste0("ARIMA(", best_order[1], ",", best_order[2], ",", best_order[3], ")"),
"| AIC =", round(model_best$aic, 2), "\n")
}
##
## ✔ Model final dipilih: auto.arima → ARIMA(2,2,1) | AIC = 467.43
# 3d. DIAGNOSTIK RESIDUAL ARIMA
cat("\n--- Diagnostik Residual Model ARIMA ---\n")
##
## --- Diagnostik Residual Model ARIMA ---
# Ljung-Box test: menguji apakah residual bersifat white noise
cat("\n--- Uji Ljung-Box pada Residual ARIMA ---\n")
##
## --- Uji Ljung-Box pada Residual ARIMA ---
lb_test <- Box.test(residuals(model_arima), lag = 10, type = "Ljung-Box")
print(lb_test)
##
## Box-Ljung test
##
## data: residuals(model_arima)
## X-squared = 10.813, df = 10, p-value = 0.3722
cat("Interpretasi: ")
## Interpretasi:
if (lb_test$p.value > 0.05) {
cat("p-value > 0.05 → Residual bersifat white noise (model memadai).\n")
} else {
cat("p-value <= 0.05 → Residual masih mengandung autokorelasi.\n")
}
## p-value > 0.05 → Residual bersifat white noise (model memadai).
# Plot diagnostik
checkresiduals(
model_arima,
test = FALSE,
main = "Diagnostik Residual Model ARIMA"
)
# ============================================================================
# INTERPRETASI ARIMA:
# Model ARIMA yang dipilih menangkap pola autokorelasi dalam data emisi.
# AIC/BIC yang lebih rendah menunjukkan model yang lebih parsimoni.
# Residual yang white noise mengkonfirmasi bahwa model telah menangkap
# seluruh pola sistematik dalam data.
# ============================================================================
# ==========================
# 4. UJI EFEK ARCH & MODEL GARCH
# ==========================
cat("\n\n", paste(rep("=", 60), collapse = ""), "\n")
##
##
## ============================================================
cat("BAGIAN 4: UJI EFEK ARCH & MODEL GARCH\n")
## BAGIAN 4: UJI EFEK ARCH & MODEL GARCH
cat(paste(rep("=", 60), collapse = ""), "\n\n")
## ============================================================
# 4a. UJI ARCH-LM PADA RESIDUAL ARIMA
residual_arima <- residuals(model_arima)
# Uji ARCH-LM menggunakan paket FinTS
cat("--- Uji ARCH-LM (Engle's LM Test) pada Residual ARIMA ---\n")
## --- Uji ARCH-LM (Engle's LM Test) pada Residual ARIMA ---
arch_test <- ArchTest(residual_arima, lags = 5, demean = TRUE)
print(arch_test)
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: residual_arima
## Chi-squared = 27.932, df = 5, p-value = 3.753e-05
cat("Interpretasi: ")
## Interpretasi:
if (arch_test$p.value > 0.05) {
cat("p-value > 0.05 → TIDAK terdeteksi efek ARCH.\n")
cat(" → Volatilitas residual konstan (homoskedastik).\n")
cat(" → Model ARIMA sudah cukup, GARCH opsional.\n")
arch_detected <- FALSE
} else {
cat("p-value <= 0.05 → Terdeteksi efek ARCH (heteroskedastisitas).\n")
cat(" → Volatilitas residual berubah-ubah seiring waktu.\n")
cat(" → Model GARCH diperlukan untuk menangkap pola volatilitas.\n")
arch_detected <- TRUE
}
## p-value <= 0.05 → Terdeteksi efek ARCH (heteroskedastisitas).
## → Volatilitas residual berubah-ubah seiring waktu.
## → Model GARCH diperlukan untuk menangkap pola volatilitas.
# 4b. ESTIMASI MODEL GARCH (tetap dilakukan untuk kelengkapan analisis)
cat("\n--- Estimasi Model GARCH(1,1) ---\n")
##
## --- Estimasi Model GARCH(1,1) ---
cat("(Model GARCH diestimasi untuk kelengkapan analisis jurnal)\n\n")
## (Model GARCH diestimasi untuk kelengkapan analisis jurnal)
# Spesifikasi GARCH(1,1) dengan mean model ARMA
# Mengambil order AR dan MA dari model ARIMA terpilih
arima_order <- arimaorder(model_arima)
p_ar <- arima_order["p"]
q_ma <- arima_order["q"]
# Spesifikasi model menggunakan rugarch
# sGARCH = standard GARCH; norm = distribusi normal
garch_spec <- ugarchspec(
variance.model = list(
model = "sGARCH", # Model GARCH standar
garchOrder = c(1, 1) # GARCH(1,1): alpha + beta
),
mean.model = list(
armaOrder = c(p_ar, q_ma), # Mean equation: ARMA dari ARIMA
include.mean = TRUE
),
distribution.model = "norm" # Distribusi error: Normal
)
# Estimasi model pada data yang sudah di-differencing sesuai order d
d_order <- arima_order["d"]
if (d_order > 0) {
co2_diff_garch <- diff(co2_ts, differences = d_order)
} else {
co2_diff_garch <- co2_ts
}
# Fit model GARCH
garch_fit <- tryCatch({
ugarchfit(spec = garch_spec, data = co2_diff_garch, solver = "hybrid")
}, error = function(e) {
cat(" ⚠ GARCH(1,1) dengan ARMA(", p_ar, ",", q_ma, ") gagal.\n")
cat(" → Mencoba GARCH(1,1) dengan mean konstan...\n\n")
# Fallback: GARCH(1,1) tanpa komponen ARMA
spec_simple <- ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "norm"
)
ugarchfit(spec = spec_simple, data = co2_diff_garch, solver = "hybrid")
})
## Warning in .sgarchfit(spec = spec, data = data, out.sample = out.sample, :
## ugarchfit-->waring: using less than 100 data
## points for estimation
cat("\n=== RINGKASAN MODEL GARCH ===\n")
##
## === RINGKASAN MODEL GARCH ===
show(garch_fit)
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(2,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.38502 0.084123 4.57687 0.000005
## ar1 0.11880 0.151878 0.78219 0.434103
## ar2 -0.32069 0.161855 -1.98135 0.047552
## ma1 -1.00000 0.110894 -9.01764 0.000000
## omega 0.00000 0.027946 0.00000 1.000000
## alpha1 0.31491 0.139048 2.26477 0.023527
## beta1 0.80756 0.079991 10.09558 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.38502 0.109931 3.50239 0.000461
## ar1 0.11880 0.160394 0.74066 0.458900
## ar2 -0.32069 0.189466 -1.69261 0.090531
## ma1 -1.00000 0.095009 -10.52531 0.000000
## omega 0.00000 0.000029 0.00000 1.000000
## alpha1 0.31491 0.083214 3.78438 0.000154
## beta1 0.80756 0.061888 13.04865 0.000000
##
## LogLikelihood : -214.5257
##
## Information Criteria
## ------------------------------------
##
## Akaike 8.3595
## Bayes 8.6197
## Shibata 8.3297
## Hannan-Quinn 8.4595
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.657 0.4176
## Lag[2*(p+q)+(p+q)-1][8] 3.080 0.9950
## Lag[4*(p+q)+(p+q)-1][14] 6.768 0.6060
## d.o.f=3
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.854 0.09116
## Lag[2*(p+q)+(p+q)-1][5] 4.764 0.17300
## Lag[4*(p+q)+(p+q)-1][9] 7.674 0.14903
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.03437 0.500 2.000 0.8529
## ARCH Lag[5] 1.66013 1.440 1.667 0.5512
## ARCH Lag[7] 3.75252 2.315 1.543 0.3837
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 4.551
## Individual Statistics:
## mu 0.03604
## ar1 0.12107
## ar2 0.36800
## ma1 0.30139
## omega 0.16183
## alpha1 0.11552
## beta1 0.31751
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.1311 0.8963
## Negative Sign Bias 1.0141 0.3156
## Positive Sign Bias 0.3655 0.7164
## Joint Effect 1.4392 0.6964
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 15.30 0.7032
## 2 30 24.55 0.7015
## 3 40 32.28 0.7681
## 4 50 46.06 0.5932
##
##
## Elapsed time : 0.4542799
# Koefisien GARCH
cat("\n--- Koefisien Model GARCH ---\n")
##
## --- Koefisien Model GARCH ---
print(coef(garch_fit))
## mu ar1 ar2 ma1 omega
## 3.850210e-01 1.187975e-01 -3.206920e-01 -1.000000e+00 2.220446e-16
## alpha1 beta1
## 3.149128e-01 8.075602e-01
# Informasi kriteria
cat("\n--- Kriteria Informasi ---\n")
##
## --- Kriteria Informasi ---
print(infocriteria(garch_fit))
##
## Akaike 8.359459
## Bayes 8.619686
## Shibata 8.329709
## Hannan-Quinn 8.459530
# ============================================================================
# INTERPRETASI GARCH:
# - omega (ω): konstanta dalam persamaan varians
# - alpha1 (α₁): koefisien ARCH, mengukur dampak shock periode sebelumnya
# terhadap volatilitas saat ini (short-term shock)
# - beta1 (β₁): koefisien GARCH, mengukur persistensi volatilitas
# (long-term memory)
# - α₁ + β₁ < 1 → volatilitas bersifat mean-reverting (stabil)
# - α₁ + β₁ ≈ 1 → volatilitas sangat persisten (mendekati IGARCH)
#
# Jika efek ARCH tidak signifikan (Bagian 4a), maka model ARIMA murni
# sudah cukup, dan GARCH hanya bersifat komplementer.
# ============================================================================
# ==========================
# 5. FORECASTING & VISUALISASI
# ==========================
cat("\n\n", paste(rep("=", 60), collapse = ""), "\n")
##
##
## ============================================================
cat("BAGIAN 5: PERAMALAN (FORECASTING) 10 TAHUN KE DEPAN\n")
## BAGIAN 5: PERAMALAN (FORECASTING) 10 TAHUN KE DEPAN
cat(paste(rep("=", 60), collapse = ""), "\n\n")
## ============================================================
# 5a. PERAMALAN ARIMA (10 Tahun: 2025–2034)
h_forecast <- 10 # horizon peramalan
fc_arima <- forecast(model_arima, h = h_forecast, level = c(80, 95))
cat("=== HASIL PERAMALAN ARIMA (2025–2034) ===\n")
## === HASIL PERAMALAN ARIMA (2025–2034) ===
print(fc_arima)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2025 851.1820 827.2156 875.1483 814.5286 887.8353
## 2026 895.5798 855.3687 935.7908 834.0823 957.0773
## 2027 940.4473 891.7145 989.1801 865.9169 1014.9776
## 2028 982.4429 924.7868 1040.0989 894.2655 1070.6202
## 2029 1024.1271 953.2523 1095.0018 915.7335 1132.5207
## 2030 1067.3317 982.0085 1152.6548 936.8412 1197.8222
## 2031 1110.7345 1011.8837 1209.5853 959.5553 1261.9137
## 2032 1153.3331 1040.4873 1266.1789 980.7503 1325.9159
## 2033 1195.8091 1067.5053 1324.1128 999.5855 1392.0327
## 2034 1238.7100 1094.1493 1383.2707 1017.6235 1459.7965
# Membuat dataframe untuk plotting
fc_df <- data.frame(
year = 2025:(2024 + h_forecast),
forecast = as.numeric(fc_arima$mean),
lo80 = as.numeric(fc_arima$lower[, 1]),
hi80 = as.numeric(fc_arima$upper[, 1]),
lo95 = as.numeric(fc_arima$lower[, 2]),
hi95 = as.numeric(fc_arima$upper[, 2])
)
cat("\n=== TABEL PERAMALAN ===\n")
##
## === TABEL PERAMALAN ===
print(fc_df %>% mutate(across(where(is.numeric) & !matches("year"), ~round(., 2))))
## year forecast lo80 hi80 lo95 hi95
## 1 2025 851.18 827.22 875.15 814.53 887.84
## 2 2026 895.58 855.37 935.79 834.08 957.08
## 3 2027 940.45 891.71 989.18 865.92 1014.98
## 4 2028 982.44 924.79 1040.10 894.27 1070.62
## 5 2029 1024.13 953.25 1095.00 915.73 1132.52
## 6 2030 1067.33 982.01 1152.65 936.84 1197.82
## 7 2031 1110.73 1011.88 1209.59 959.56 1261.91
## 8 2032 1153.33 1040.49 1266.18 980.75 1325.92
## 9 2033 1195.81 1067.51 1324.11 999.59 1392.03
## 10 2034 1238.71 1094.15 1383.27 1017.62 1459.80
# 5b. PERAMALAN GARCH (Conditional Volatility)
fc_garch <- ugarchforecast(garch_fit, n.ahead = h_forecast)
cat("\n=== PERAMALAN VOLATILITAS GARCH (10 Tahun) ===\n")
##
## === PERAMALAN VOLATILITAS GARCH (10 Tahun) ===
cat("Series (Mean Forecast - differenced scale):\n")
## Series (Mean Forecast - differenced scale):
print(fitted(fc_garch))
## 2024-01-01
## T+1 -21.27355484
## T+2 1.21293112
## T+3 7.42910763
## T+4 0.95633663
## T+5 -1.80609068
## T+6 -0.05849394
## T+7 1.03500463
## T+8 0.60446914
## T+9 0.20264631
## T+10 0.29298008
cat("\nSigma (Conditional Volatility Forecast):\n")
##
## Sigma (Conditional Volatility Forecast):
print(sigma(fc_garch))
## 2024-01-01
## T+1 42.72754
## T+2 45.26847
## T+3 47.96051
## T+4 50.81264
## T+5 53.83438
## T+6 57.03582
## T+7 60.42764
## T+8 64.02116
## T+9 67.82839
## T+10 71.86203
# 5c. VISUALISASI UTAMA — GRAFIK JURNAL
cat("\n--- Membuat visualisasi untuk jurnal ---\n")
##
## --- Membuat visualisasi untuk jurnal ---
# Warna palette yang konsisten dan profesional
color_hist <- "#2C3E50" # biru tua (data historis)
color_fc <- "#E74C3C" # merah (forecast)
color_ci95 <- "#F5B7B1" # merah muda terang (CI 95%)
color_ci80 <- "#EC7063" # merah muda (CI 80%)
color_grid <- "#ECF0F1" # abu terang (grid)
# Data historis untuk plotting
hist_df <- df %>% select(year, co2_mt)
# PLOT 1: Grafik Utama — Historical + Forecast
p_main <- ggplot() +
# Confidence interval 95%
geom_ribbon(
data = fc_df,
aes(x = year, ymin = lo95, ymax = hi95),
fill = color_ci95, alpha = 0.5
) +
# Confidence interval 80%
geom_ribbon(
data = fc_df,
aes(x = year, ymin = lo80, ymax = hi80),
fill = color_ci80, alpha = 0.5
) +
# Garis data historis
geom_line(
data = hist_df,
aes(x = year, y = co2_mt),
color = color_hist, linewidth = 0.9
) +
geom_point(
data = hist_df,
aes(x = year, y = co2_mt),
color = color_hist, size = 1.2, alpha = 0.6
) +
# Garis peramalan
geom_line(
data = fc_df,
aes(x = year, y = forecast),
color = color_fc, linewidth = 1, linetype = "dashed"
) +
geom_point(
data = fc_df,
aes(x = year, y = forecast),
color = color_fc, size = 2, shape = 17
) +
# Garis pemisah historis vs forecast
geom_vline(
xintercept = 2024.5,
linetype = "dotted", color = "grey40", linewidth = 0.5
) +
# Anotasi
annotate(
"text", x = 2010, y = max(fc_df$hi95) * 0.95,
label = "Data Historis",
color = color_hist, fontface = "bold", size = 3.5, hjust = 0
) +
annotate(
"text", x = 2028, y = max(fc_df$hi95) * 0.95,
label = "Peramalan ARIMA",
color = color_fc, fontface = "bold", size = 3.5, hjust = 0
) +
# Formatting
scale_x_continuous(
breaks = seq(1970, 2035, by = 5),
expand = expansion(mult = c(0.02, 0.02))
) +
scale_y_continuous(
labels = label_comma(suffix = " Mt"),
expand = expansion(mult = c(0.05, 0.1))
) +
labs(
title = "Proyeksi Emisi CO\u2082 Fossil Indonesia (1970\u20132034)",
subtitle = paste0(
"Model: ", paste0(
"ARIMA(", arimaorder(model_arima)["p"], ",",
arimaorder(model_arima)["d"], ",",
arimaorder(model_arima)["q"], ")"
),
" | Horizon: 10 Tahun | CI: 80% & 95%"
),
x = "Tahun",
y = "Emisi CO\u2082 (Juta Ton)",
caption = "Sumber Data: Worldometer | Analisis: R (forecast, rugarch)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14, hjust = 0,
margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 10, hjust = 0,
margin = margin(b = 12)),
plot.caption = element_text(color = "grey50", size = 8, hjust = 1,
margin = margin(t = 10)),
axis.title = element_text(face = "bold", size = 10),
axis.text = element_text(size = 9),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = color_grid, linewidth = 0.3),
plot.margin = margin(15, 15, 10, 10)
)
print(p_main)
# Simpan plot utama (resolusi tinggi untuk jurnal)
ggsave(
filename = "d:/Training/CO2_Forecast_Main.png",
plot = p_main,
width = 10, height = 6, dpi = 300, bg = "white"
)
cat("✔ Plot utama disimpan: d:/Training/CO2_Forecast_Main.png\n")
## ✔ Plot utama disimpan: d:/Training/CO2_Forecast_Main.png
# PLOT 2: ACF/PACF menggunakan ggplot2 (versi estetis untuk jurnal)
# Hitung ACF dan PACF secara manual untuk ggplot
acf_vals <- acf(co2_diff1, plot = FALSE, lag.max = 20)
pacf_vals <- pacf(co2_diff1, plot = FALSE, lag.max = 20)
acf_df <- data.frame(
lag = as.numeric(acf_vals$lag[-1]),
acf = as.numeric(acf_vals$acf[-1])
)
pacf_df <- data.frame(
lag = as.numeric(pacf_vals$lag),
pacf = as.numeric(pacf_vals$acf)
)
# Batas signifikansi (95% CI ≈ ±1.96/√n)
n_obs <- length(co2_diff1)
ci_bound <- 1.96 / sqrt(n_obs)
p_acf <- ggplot(acf_df, aes(x = lag, y = acf)) +
geom_hline(yintercept = 0, color = "grey40", linewidth = 0.4) +
geom_hline(yintercept = c(-ci_bound, ci_bound),
linetype = "dashed", color = "#3498DB", linewidth = 0.4) +
geom_segment(aes(xend = lag, yend = 0), color = color_hist, linewidth = 0.7) +
geom_point(color = color_hist, size = 2) +
labs(title = "ACF — Data Differenced (d=1)", x = "Lag", y = "Autokorelasi") +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0.5),
panel.grid.minor = element_blank()
)
p_pacf <- ggplot(pacf_df, aes(x = lag, y = pacf)) +
geom_hline(yintercept = 0, color = "grey40", linewidth = 0.4) +
geom_hline(yintercept = c(-ci_bound, ci_bound),
linetype = "dashed", color = "#3498DB", linewidth = 0.4) +
geom_segment(aes(xend = lag, yend = 0), color = "#8E44AD", linewidth = 0.7) +
geom_point(color = "#8E44AD", size = 2) +
labs(title = "PACF — Data Differenced (d=1)", x = "Lag", y = "Autokorelasi Parsial") +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0.5),
panel.grid.minor = element_blank()
)
p_acf_pacf <- p_acf + p_pacf +
plot_annotation(
title = "Correlogram Emisi CO\u2082 Indonesia",
caption = "Garis biru putus-putus: batas signifikansi 95%",
theme = theme(
plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
plot.caption = element_text(color = "grey50", size = 8, hjust = 0.5)
)
)
print(p_acf_pacf)
ggsave(
filename = "d:/Training/CO2_ACF_PACF.png",
plot = p_acf_pacf,
width = 10, height = 5, dpi = 300, bg = "white"
)
cat("✔ Plot ACF/PACF disimpan: d:/Training/CO2_ACF_PACF.png\n")
## ✔ Plot ACF/PACF disimpan: d:/Training/CO2_ACF_PACF.png
# PLOT 3: Diagnostik Residual ARIMA (ggplot2)
resid_df <- data.frame(
year = time(residual_arima),
residual = as.numeric(residual_arima)
)
p_resid_ts <- ggplot(resid_df, aes(x = year, y = residual)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_line(color = color_hist, linewidth = 0.5) +
geom_point(color = color_hist, size = 1, alpha = 0.6) +
labs(title = "Residual Time Series", x = "Tahun", y = "Residual (Mt)") +
theme_minimal(base_size = 10) +
theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 10),
panel.grid.minor = element_blank())
p_resid_hist <- ggplot(resid_df, aes(x = residual)) +
geom_histogram(
aes(y = after_stat(density)),
bins = 15, fill = color_hist, alpha = 0.7, color = "white"
) +
geom_density(color = color_fc, linewidth = 0.8) +
labs(title = "Distribusi Residual", x = "Residual (Mt)", y = "Density") +
theme_minimal(base_size = 10) +
theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 10),
panel.grid.minor = element_blank())
p_resid_qq <- ggplot(resid_df, aes(sample = residual)) +
stat_qq(color = color_hist, size = 1.5, alpha = 0.7) +
stat_qq_line(color = color_fc, linewidth = 0.7) +
labs(title = "Q-Q Plot Residual", x = "Theoretical Quantiles", y = "Sample Quantiles") +
theme_minimal(base_size = 10) +
theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 10),
panel.grid.minor = element_blank())
p_diagnostics <- p_resid_ts + p_resid_hist + p_resid_qq +
plot_layout(ncol = 3) +
plot_annotation(
title = "Diagnostik Residual Model ARIMA",
theme = theme(plot.title = element_text(face = "bold", size = 13, hjust = 0.5))
)
print(p_diagnostics)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
ggsave(
filename = "d:/Training/CO2_Residual_Diagnostics.png",
plot = p_diagnostics,
width = 12, height = 4.5, dpi = 300, bg = "white"
)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
cat("✔ Plot diagnostik disimpan: d:/Training/CO2_Residual_Diagnostics.png\n")
## ✔ Plot diagnostik disimpan: d:/Training/CO2_Residual_Diagnostics.png
# PLOT 4: Conditional Volatility (dari GARCH)
sigma_hist <- sigma(garch_fit)
vol_df <- data.frame(
year = time(co2_diff_garch),
volatility = as.numeric(sigma_hist)
)
p_vol <- ggplot(vol_df, aes(x = year, y = volatility)) +
geom_area(fill = "#F39C12", alpha = 0.3) +
geom_line(color = "#E67E22", linewidth = 0.8) +
labs(
title = "Conditional Volatility Emisi CO\u2082 — Model GARCH(1,1)",
subtitle = "Estimasi volatilitas bersyarat dari residual ARIMA",
x = "Tahun", y = "Conditional Std. Dev. (Mt)",
caption = "Sumber: Model GARCH(1,1) via rugarch"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 13, hjust = 0),
plot.subtitle = element_text(color = "grey40", size = 9, hjust = 0),
plot.caption = element_text(color = "grey50", size = 8),
panel.grid.minor = element_blank()
)
print(p_vol)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
ggsave(
filename = "d:/Training/CO2_GARCH_Volatility.png",
plot = p_vol,
width = 10, height = 5, dpi = 300, bg = "white"
)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
cat("✔ Plot volatilitas GARCH disimpan: d:/Training/CO2_GARCH_Volatility.png\n")
## ✔ Plot volatilitas GARCH disimpan: d:/Training/CO2_GARCH_Volatility.png
# ==========================
# 6. RINGKASAN AKHIR
# ==========================
cat("\n\n", paste(rep("=", 60), collapse = ""), "\n")
##
##
## ============================================================
cat("RINGKASAN ANALISIS\n")
## RINGKASAN ANALISIS
cat(paste(rep("=", 60), collapse = ""), "\n\n")
## ============================================================
cat("1. REGRESI LINEAR:\n")
## 1. REGRESI LINEAR:
cat(" - R²:", round(summary(model_lm)$r.squared, 4), "\n")
## - R²: 0.9474
cat(" - Slope (peningkatan per tahun):",
round(coef(model_lm)[2], 2), "Mt/tahun\n")
## - Slope (peningkatan per tahun): 13.02 Mt/tahun
cat(" - Autokorelasi terdeteksi (DW test) → regresi OLS tidak memadai.\n\n")
## - Autokorelasi terdeteksi (DW test) → regresi OLS tidak memadai.
cat("2. MODEL ARIMA:\n")
## 2. MODEL ARIMA:
cat(" - Model terpilih:", paste0(
"ARIMA(", arimaorder(model_arima)["p"], ",",
arimaorder(model_arima)["d"], ",",
arimaorder(model_arima)["q"], ")\n"))
## - Model terpilih: ARIMA(2,2,1)
cat(" - AIC:", round(model_arima$aic, 2), "\n")
## - AIC: 467.43
cat(" - BIC:", round(model_arima$bic, 2), "\n")
## - BIC: 475.31
cat(" - Residual white noise (Ljung-Box p =", round(lb_test$p.value, 4), ")\n\n")
## - Residual white noise (Ljung-Box p = 0.3722 )
cat("3. EFEK ARCH/GARCH:\n")
## 3. EFEK ARCH/GARCH:
cat(" - ARCH-LM Test p-value:", round(arch_test$p.value, 4), "\n")
## - ARCH-LM Test p-value: 0
if (arch_detected) {
cat(" - Efek ARCH TERDETEKSI → Model GARCH(1,1) diestimasi.\n")
} else {
cat(" - Efek ARCH TIDAK terdeteksi → GARCH bersifat komplementer.\n")
}
## - Efek ARCH TERDETEKSI → Model GARCH(1,1) diestimasi.
cat(" - GARCH alpha1 + beta1:",
round(sum(coef(garch_fit)[c("alpha1", "beta1")]), 4), "\n\n")
## - GARCH alpha1 + beta1: 1.1225
cat("4. PERAMALAN 2025–2034:\n")
## 4. PERAMALAN 2025–2034:
for (i in 1:nrow(fc_df)) {
cat(sprintf(" %d: %.2f Mt [95%% CI: %.2f — %.2f]\n",
fc_df$year[i], fc_df$forecast[i], fc_df$lo95[i], fc_df$hi95[i]))
}
## 2025: 851.18 Mt [95% CI: 814.53 — 887.84]
## 2026: 895.58 Mt [95% CI: 834.08 — 957.08]
## 2027: 940.45 Mt [95% CI: 865.92 — 1014.98]
## 2028: 982.44 Mt [95% CI: 894.27 — 1070.62]
## 2029: 1024.13 Mt [95% CI: 915.73 — 1132.52]
## 2030: 1067.33 Mt [95% CI: 936.84 — 1197.82]
## 2031: 1110.73 Mt [95% CI: 959.56 — 1261.91]
## 2032: 1153.33 Mt [95% CI: 980.75 — 1325.92]
## 2033: 1195.81 Mt [95% CI: 999.59 — 1392.03]
## 2034: 1238.71 Mt [95% CI: 1017.62 — 1459.80]
cat("\n5. OUTPUT FILES:\n")
##
## 5. OUTPUT FILES:
cat(" - d:/Training/CO2_Forecast_Main.png\n")
## - d:/Training/CO2_Forecast_Main.png
cat(" - d:/Training/CO2_ACF_PACF.png\n")
## - d:/Training/CO2_ACF_PACF.png
cat(" - d:/Training/CO2_Residual_Diagnostics.png\n")
## - d:/Training/CO2_Residual_Diagnostics.png
cat(" - d:/Training/CO2_GARCH_Volatility.png\n")
## - d:/Training/CO2_GARCH_Volatility.png
cat("\n", paste(rep("=", 60), collapse = ""), "\n")
##
## ============================================================
cat("ANALISIS SELESAI.\n")
## ANALISIS SELESAI.
cat(paste(rep("=", 60), collapse = ""), "\n")
## ============================================================