Sumber Data
library(readxl)
library(forecast)
library(urca)
library(tseries)
library(dynlm)
library(lmtest)
library(sandwich)
library(knitr)
Hasil Pemodelan
Membaca dan menyiapkan data
file_data <- "C:/TUGAS DAN MATERI STIS/Tingkat 2/Semester 4/Analisis Runtun WaktU/Project UAS/Studi Kasus 3/pdb_triwulanan.xlsx"
raw <- read_excel(file_data, col_names = FALSE, .name_repair = "minimal")
header_triwulan <- as.character(unlist(raw[3, ]))
kolom_tw <- which(grepl("Triwulan", header_triwulan)) # buang kolom "Tahunan"
ambil <- function(pola) {
baris <- which(grepl(pola, raw[[1]]))
suppressWarnings(as.numeric(unlist(raw[baris, kolom_tw]))) # "-" jadi NA
}
C <- ambil("Konsumsi Rumahtangga")
Y <- ambil("PRODUK DOMESTIK BRUTO")
ok <- !is.na(C) & !is.na(Y)
C <- ts(C[ok], start = c(2010, 1), frequency = 4)
Y <- ts(Y[ok], start = c(2010, 1), frequency = 4)
cat("Jumlah observasi:", length(C), "| periode:",
paste(start(C), collapse = "Q"), "s.d.", paste(end(C), collapse = "Q"), "\n")
## Jumlah observasi: 66 | periode: 2010Q1 s.d. 2026Q2
kable(head(data.frame(Periode = as.character(zoo::as.yearqtr(time(C))),
Konsumsi_RT = as.numeric(C), PDB = as.numeric(Y)), 8),
format.args = list(big.mark = "."), caption = "Cuplikan data (miliar rupiah, ADHK 2010)")
Cuplikan data (miliar rupiah, ADHK 2010)
| 2010 Q1 |
926.097.5 |
1.642.356 |
| 2010 Q2 |
937.226.8 |
1.709.132 |
| 2010 Q3 |
959.649.9 |
1.775.110 |
| 2010 Q4 |
963.088.7 |
1.737.535 |
| 2011 Q1 |
964.261.9 |
1.748.731 |
| 2011 Q2 |
980.156.6 |
1.816.268 |
| 2011 Q3 |
1.016.155.4 |
1.881.850 |
| 2011 Q4 |
1.016.714.7 |
1.840.786 |
Uji stasioneritas
uji_ar <- function(x, nama, tipe_adf, tipe_kpss) {
adf <- ur.df(x, type = tipe_adf, lags = 4, selectlags = "AIC")
stat <- adf@teststat[1]; cv5 <- adf@cval[1, "5pct"]
kp <- kpss.test(x, null = tipe_kpss)
data.frame(Variabel = nama, `ADF stat` = round(stat, 3), `CV 5%` = cv5,
`Keputusan ADF` = ifelse(stat < cv5, "Stasioner", "Tidak stasioner"),
`KPSS p-value` = round(kp$p.value, 3),
`Keputusan KPSS` = ifelse(kp$p.value < .05, "Tidak stasioner", "Stasioner"),
check.names = FALSE)
}
tabel_ur <- rbind(
uji_ar(lnC, "ln C (level)", "trend", "Trend"),
uji_ar(lnY, "ln Y (level)", "trend", "Trend"),
uji_ar(diff(lnC), "Δ ln C (diff. 1)", "drift", "Level"),
uji_ar(diff(lnY), "Δ ln Y (diff. 1)", "drift", "Level"))
kable(tabel_ur, caption = "Uji ADF dan KPSS (catatan: p-value KPSS dibatasi 0,01–0,10)")
Uji ADF dan KPSS (catatan: p-value KPSS dibatasi
0,01–0,10)
| ln C (level) |
-1.942 |
-3.45 |
Tidak stasioner |
0.01 |
Tidak stasioner |
| ln Y (level) |
-2.122 |
-3.45 |
Tidak stasioner |
0.01 |
Tidak stasioner |
| Δ ln C (diff. 1) |
-4.941 |
-2.89 |
Stasioner |
0.10 |
Stasioner |
| Δ ln Y (diff. 1) |
-4.925 |
-2.89 |
Stasioner |
0.10 |
Stasioner |
Persamaan jangka panjang (regresi kointegrasi)
waktu <- time(lnC)
kandidat <- seq(2019, 2022.75, by = 0.25)
ssr <- sapply(kandidat, function(b) {
D <- as.numeric(waktu >= b - 1e-6); sum(resid(lm(lnC ~ lnY + D))^2)
})
tgl_break <- kandidat[which.min(ssr)]
plot(kandidat, ssr * 1e4, type = "b", pch = 16, xlab = "Awal dummy", ylab = "SSR (x 10^-4)",
main = "Pemilihan tanggal patahan struktural")
abline(v = tgl_break, col = "red", lty = 2)
cat("Tanggal patahan terpilih:", as.character(zoo::as.yearqtr(tgl_break)), "\n")
## Tanggal patahan terpilih: 2021 Q1
D <- ts(as.numeric(waktu >= tgl_break - 1e-6), start = c(2010, 1), frequency = 4)
lr <- dynlm(lnC ~ lnY + D)
summary(lr)
##
## Time series regression with "ts" data:
## Start = 2010(1), End = 2026(2)
##
## Call:
## dynlm(formula = lnC ~ lnY + D)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0087579 -0.0022956 0.0001153 0.0017527 0.0169457
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.196455 0.054828 -3.583 0.000662 ***
## lnY 0.971736 0.003751 259.039 < 0.0000000000000002 ***
## D -0.020003 0.001646 -12.156 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004049 on 63 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 7.559e+04 on 2 and 63 DF, p-value: < 0.00000000000000022
Uji kointegrasi Engle–Granger
ect <- ts(resid(lr), start = c(2010, 1), frequency = 4)
eg <- ur.df(ect, type = "none", lags = 4, selectlags = "BIC")
summary(eg)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0080400 -0.0011218 0.0003518 0.0012444 0.0096004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.5285 0.1397 -3.782 0.000365 ***
## z.diff.lag -0.1137 0.1269 -0.896 0.373977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002805 on 59 degrees of freedom
## Multiple R-squared: 0.3139, Adjusted R-squared: 0.2906
## F-statistic: 13.5 on 2 and 59 DF, p-value: 0.00001492
##
##
## Value of test-statistic is: -3.7822
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.6 -1.95 -1.61
plot(ect, type = "o", pch = 16, cex = .6, col = "purple",
main = "Residual jangka panjang (ECT): simpangan dari keseimbangan", ylab = "u_t", xlab = "")
abline(h = 0, lty = 2)
tau <- eg@teststat[1]
cat("Statistik ADF residual =", round(tau, 3),
"| nilai kritis MacKinnon 5% ≈ -3,37 ->",
ifelse(tau < -3.37, "TOLAK H0: residual stasioner, C dan Y BERKOINTEGRASI",
"Gagal tolak H0"), "\n")
## Statistik ADF residual = -3.782 | nilai kritis MacKinnon 5% ≈ -3,37 -> TOLAK H0: residual stasioner, C dan Y BERKOINTEGRASI
Estimasi Error Correction Model
ecm <- dynlm(d(lnC) ~ d(lnY) + L(ect, 1))
summary(ecm)
##
## Time series regression with "ts" data:
## Start = 2010(2), End = 2026(2)
##
## Call:
## dynlm(formula = d(lnC) ~ d(lnY) + L(ect, 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0086942 -0.0004715 0.0006710 0.0016902 0.0042509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001378 0.000531 -2.596 0.0118 *
## d(lnY) 1.042131 0.034840 29.912 < 0.0000000000000002 ***
## L(ect, 1) -0.391369 0.087135 -4.492 0.0000314 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002782 on 62 degrees of freedom
## Multiple R-squared: 0.9392, Adjusted R-squared: 0.9373
## F-statistic: 479 on 2 and 62 DF, p-value: < 0.00000000000000022
coeftest(ecm, vcov = NeweyWest(ecm))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0013782 0.0005392 -2.5560 0.013054 *
## d(lnY) 1.0421307 0.0272499 38.2435 < 0.00000000000000022 ***
## L(ect, 1) -0.3913686 0.1271170 -3.0788 0.003095 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t_ect <- summary(ecm)$coefficients["L(ect, 1)", "t value"]
cat("t-stat ECT =", round(t_ect, 3), "->",
ifelse(t_ect < -3.19, "signifikan (mendukung adanya kointegrasi)", "tidak signifikan"), "\n")
## t-stat ECT = -4.492 -> signifikan (mendukung adanya kointegrasi)
Uji diagnostik
diag_tab <- data.frame(
Uji = c("Breusch-Godfrey (lag 1)", "Breusch-Godfrey (lag 4)",
"Breusch-Pagan", "Jarque-Bera", "Ramsey RESET"),
H0 = c("Tidak ada autokorelasi", "Tidak ada autokorelasi", "Homoskedastis",
"Residual normal", "Spesifikasi linear tepat"),
`p-value` = c(bgtest(ecm, 1)$p.value, bgtest(ecm, 4)$p.value, bptest(ecm)$p.value,
jarque.bera.test(resid(ecm))$p.value, resettest(ecm)$p.value),
check.names = FALSE)
diag_tab$Keputusan <- ifelse(diag_tab$`p-value` < 0.05, "Tolak H0", "Gagal tolak H0")
diag_tab$`p-value` <- format(signif(diag_tab$`p-value`, 3), scientific = FALSE, drop0trailing = TRUE)
kable(diag_tab, caption = "Uji asumsi residual ECM (α = 5%)")
Uji asumsi residual ECM (α = 5%)
| Breusch-Godfrey (lag 1) |
Tidak ada autokorelasi |
0.756 |
Gagal tolak H0 |
| Breusch-Godfrey (lag 4) |
Tidak ada autokorelasi |
0.522 |
Gagal tolak H0 |
| Breusch-Pagan |
Homoskedastis |
0.41 |
Gagal tolak H0 |
| Jarque-Bera |
Residual normal |
0.000000387 |
Tolak H0 |
| Ramsey RESET |
Spesifikasi linear tepat |
0.938 |
Gagal tolak H0 |
par(mfrow = c(2, 2), mar = c(4, 4, 2.5, 1))
plot(resid(ecm), type = "o", pch = 16, cex = .5, main = "Residual ECM", ylab = "", xlab = "")
abline(h = 0, lty = 2)
acf(as.numeric(resid(ecm)), main = "ACF residual", lag.max = 12)
hist(resid(ecm), breaks = 15, col = "lightblue", main = "Histogram residual", xlab = "")
qqnorm(resid(ecm)); qqline(resid(ecm), col = "red")
r <- resid(ecm)
kable(data.frame(Periode = as.character(zoo::as.yearqtr(time(r)))[order(r)[1:3]],
Residual = round(sort(as.numeric(r))[1:3], 4)),
caption = "Tiga residual paling ekstrem")
Tiga residual paling ekstrem
| 2010 Q3 |
-0.0087 |
| 2021 Q1 |
-0.0076 |
| 2020 Q4 |
-0.0076 |
Kecepatan penyesuaian
lambda <- coef(ecm)["L(ect, 1)"]
half_life <- log(0.5) / log(1 + lambda)
t90 <- log(0.1) / log(1 + lambda)
cat("Koefisien ECT (lambda) :", round(lambda, 4), "\n",
"Porsi ketidakseimbangan dikoreksi per triwulan:", round(-lambda * 100, 1), "%\n",
"Half-life :", round(half_life, 2), "triwulan\n",
"Waktu sampai 90% terkoreksi :", round(t90, 2), "triwulan\n")
## Koefisien ECT (lambda) : -0.3914
## Porsi ketidakseimbangan dikoreksi per triwulan: 39.1 %
## Half-life : 1.4 triwulan
## Waktu sampai 90% terkoreksi : 4.64 triwulan
# Simulasi: kalau hari ini konsumsi 1% di atas keseimbangan, bagaimana simpangannya menyusut?
h <- 0:10
simpangan <- 1 * (1 + lambda)^h
barplot(simpangan, names.arg = paste0("t+", h), col = "steelblue",
main = "Sisa simpangan (%) dari keseimbangan setelah guncangan 1%", ylab = "%")
fit <- fitted(ecm)
ts.plot(diff(lnC), fit,
col = c("black", "red"), lwd = c(1.5, 1.5), lty = c(1, 2),
main = "Δ ln C aktual (hitam) vs fitted ECM (merah putus-putus)", ylab = "Δ ln C")
ringkas <- data.frame(
Komponen = c("Elastisitas jangka panjang (β1)", "Dummy pasca-pandemi (δ)",
"Elastisitas jangka pendek (γ)", "Koefisien ECT (λ)", "R² ECM"),
Estimasi = round(c(coef(lr)[2], coef(lr)[3], coef(ecm)[2], coef(ecm)[3],
summary(ecm)$r.squared), 4))
kable(ringkas, row.names = FALSE, caption = "Ringkasan hasil")
Ringkasan hasil
| Elastisitas jangka panjang (β1) |
0.9717 |
| Dummy pasca-pandemi (δ) |
-0.0200 |
| Elastisitas jangka pendek (γ) |
1.0421 |
| Koefisien ECT (λ) |
-0.3914 |
| R² ECM |
0.9392 |
sessionInfo()
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Asia/Jakarta
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.51 sandwich_3.1-3 lmtest_0.9-40 dynlm_0.3-6
## [5] zoo_1.8-15 tseries_0.10-61 urca_1.3-4 forecast_9.0.2
## [9] readxl_1.5.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 lattice_0.22-9 digest_0.6.39
## [5] magrittr_2.0.5 evaluate_1.0.5 grid_4.6.0 RColorBrewer_1.1-3
## [9] fastmap_1.2.0 cellranger_1.1.0 jsonlite_2.0.0 Formula_1.2-6
## [13] scales_1.4.0 jquerylib_0.1.4 abind_1.4-8 cli_3.6.6
## [17] rlang_1.2.0 cachem_1.1.0 yaml_2.3.12 otel_0.2.0
## [21] tools_4.6.0 parallel_4.6.0 dplyr_1.2.1 colorspace_2.1-2
## [25] ggplot2_4.0.3 curl_7.1.0 vctrs_0.7.3 R6_2.6.1
## [29] stats4_4.6.0 lifecycle_1.0.5 car_3.1-5 pkgconfig_2.0.3
## [33] pillar_1.11.1 bslib_0.11.0 gtable_0.3.6 glue_1.8.1
## [37] quantmod_0.4.28 Rcpp_1.1.1-1.1 xfun_0.57 tibble_3.3.1
## [41] tidyselect_1.2.1 rstudioapi_0.18.0 farver_2.1.2 htmltools_0.5.9
## [45] nlme_3.1-169 carData_3.0-6 rmarkdown_2.31 xts_0.14.2
## [49] timeDate_4052.112 fracdiff_1.5-4 compiler_4.6.0 S7_0.2.2
## [53] quadprog_1.5-8 TTR_0.24.4
