## ============================================================
## UAS ANALISIS RUNTUN WAKTU TERAPAN — STUDI KASUS 3
## ============================================================
library(urca)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
df <- read.csv("data_pdb.csv")
df <- df[order(df$Tahun, df$Q), ]
n <- nrow(df); n # 62 triwulan (2011Q1-2026Q2)
## [1] 62
lnC <- log(df$Konsumsi_RT)
lnY <- log(df$PDB)
C_ts <- ts(cbind(lnC = lnC, lnY = lnY), start = c(df$Tahun[1], df$Q[1]), frequency = 4)
# Dummy musiman (Q1 = referensi), pasca-COVID, dan COVID
Q2 <- as.numeric(df$Q == 2); Q3 <- as.numeric(df$Q == 3); Q4 <- as.numeric(df$Q == 4)
post <- as.numeric(df$Tahun > 2020 | (df$Tahun == 2020 & df$Q >= 2)) # mulai 2020Q2
pQ2 <- post * Q2; pQ3 <- post * Q3; pQ4 <- post * Q4 # pola musiman pasca-COVID
covid_i <- as.numeric(df$Triwulan == "2020Q2") # puncak kontraksi
plot(C_ts, plot.type = "single", col = c("red", "navy"), lwd = 2,
ylab = "Logaritma natural", xlab = "",
main = "Ln Konsumsi Rumah Tangga dan Ln PDB (harga konstan 2010)")
legend("topleft", c("ln Konsumsi RT", "ln PDB"), col = c("red", "navy"), lwd = 2, bty = "n")
porsi <- 100 * df$Konsumsi_RT / df$PDB
summary(df[, c("Konsumsi_RT", "PDB")])
## Konsumsi_RT PDB
## Min. : 964262 Min. :1748731
## 1st Qu.:1187064 1st Qu.:2173000
## Median :1417987 Median :2614517
## Mean :1391180 Mean :2593884
## 3rd Qu.:1549032 3rd Qu.:2952280
## Max. :1887683 Max. :3576209
round(c(porsi_rata2 = mean(porsi), porsi_min = min(porsi), porsi_max = max(porsi)), 2)
## porsi_rata2 porsi_min porsi_max
## 53.76 51.87 55.23
# Porsi konsumsi per triwulan, sebelum dan sesudah COVID (pola musiman berubah)
round(rbind(pra_covid = tapply(porsi[post == 0], df$Q[post == 0], mean),
pasca_covid = tapply(porsi[post == 1], df$Q[post == 1], mean)), 2)
## 1 2 3 4
## pra_covid 54.99 53.66 53.77 54.84
## pasca_covid 53.53 53.03 52.18 52.86
round(c(pra_covid = mean(porsi[post == 0]), pasca_covid = mean(porsi[post == 1])), 2)
## pra_covid pasca_covid
## 54.33 52.91
# ADF dengan dummy musiman, lag 0-4 dipilih BIC pada sampel yang sama;
# p-value MacKinnon (urca::punitroot).
adf_musiman <- function(x, q, trend = TRUE, maxlag = 4) {
dx <- diff(x); N <- length(x)
idx <- (maxlag + 2):N # periode t yang dipakai
best <- NULL
for (p in 0:maxlag) {
X <- data.frame(xlag = x[idx - 1])
if (p > 0) for (j in 1:p) X[[paste0("dl", j)]] <- dx[idx - 1 - j]
if (trend) X$tren <- idx
X$S2 <- as.numeric(q[idx] == 2); X$S3 <- as.numeric(q[idx] == 3); X$S4 <- as.numeric(q[idx] == 4)
m <- lm(dx[idx - 1] ~ ., data = X)
if (is.null(best) || BIC(m) < best$bic)
best <- list(bic = BIC(m), tau = coef(summary(m))["xlag", "t value"], p = p)
}
pv <- punitroot(best$tau, N = length(idx), trend = if (trend) "ct" else "c", statistic = "t")
c(tau = round(best$tau, 3), lag = best$p, p_value = round(pv, 4))
}
q <- df$Q
tabel_akar_unit <- rbind(
`ln C (level, + tren)` = adf_musiman(lnC, q, trend = TRUE),
`ln Y (level, + tren)` = adf_musiman(lnY, q, trend = TRUE),
`Δ ln C (selisih)` = adf_musiman(diff(lnC), q[-1], trend = FALSE),
`Δ ln Y (selisih)` = adf_musiman(diff(lnY), q[-1], trend = FALSE))
tabel_akar_unit
## tau lag p_value
## ln C (level, + tren) -2.409 0 0.3710
## ln Y (level, + tren) -2.439 0 0.3566
## Δ ln C (selisih) -9.047 0 0.0000
## Δ ln Y (selisih) -8.535 0 0.0000
# Pembanding: Phillips-Perron (non-parametrik, tidak perlu memilih lag)
round(c(PP_lnC = pp.test(lnC)$p.value, PP_dlnC = pp.test(diff(lnC))$p.value,
PP_lnY = pp.test(lnY)$p.value, PP_dlnY = pp.test(diff(lnY))$p.value), 4)
## Warning in pp.test(diff(lnC)): p-value smaller than printed p-value
## Warning in pp.test(diff(lnY)): p-value smaller than printed p-value
## PP_lnC PP_dlnC PP_lnY PP_dlnY
## 0.5114 0.0100 0.0439 0.0100
dl <- data.frame(lnC, lnY, Q2, Q3, Q4, post)
# Nilai kritis Engle-Granger untuk ADF residual (MacKinnon, 2 variabel, dengan konstanta)
cv_eg <- function(Tn) {
b <- rbind(`1%` = c(-3.89644, -10.9519, -22.527),
`5%` = c(-3.33613, -6.1101, -6.823),
`10%`= c(-3.04445, -4.2412, -2.720))
round(b[, 1] + b[, 2] / Tn + b[, 3] / Tn^2, 3)
}
cv <- cv_eg(n); cv
## 1% 5% 10%
## -4.079 -3.436 -3.114
# Trial and error persamaan jangka panjang
spes_lr <- list(
LR1 = lnC ~ lnY,
LR2 = lnC ~ lnY + Q2 + Q3 + Q4,
LR3 = lnC ~ lnY + Q2 + Q3 + Q4 + post)
hasil_lr <- do.call(rbind, lapply(names(spes_lr), function(nm) {
m <- lm(spes_lr[[nm]], data = dl); u <- residuals(m)
tau <- ur.df(u, type = "none", lags = 4, selectlags = "BIC")@teststat[1]
data.frame(Model = nm, beta_lnY = round(coef(m)["lnY"], 4),
R2 = round(summary(m)$r.squared, 4), BIC = round(BIC(m), 2),
tau_ADF_resid = round(tau, 3), CV_5pct = cv["5%"],
Tolak_H0 = tau < cv["5%"])
}))
rownames(hasil_lr) <- NULL
hasil_lr
## Model beta_lnY R2 BIC tau_ADF_resid CV_5pct Tolak_H0
## 1 LR1 0.9361 0.9952 -356.07 -1.979 -3.436 FALSE
## 2 LR2 0.9392 0.9976 -387.25 -1.867 -3.436 FALSE
## 3 LR3 0.9792 0.9987 -423.12 -2.639 -3.436 FALSE
# Tolak H0 (residual stasioner) -> terkointegrasi.
pilih_lr <- "LR3" # BIC terkecil: dummy musiman + pergeseran level pasca-COVID
model_lr <- lm(spes_lr[[pilih_lr]], data = dl)
summary(model_lr)
##
## Call:
## lm(formula = spes_lr[[pilih_lr]], data = dl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.012996 -0.004953 0.001177 0.004012 0.017219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.295068 0.104716 -2.818 0.00667 **
## lnY 0.979241 0.007165 136.663 < 2e-16 ***
## Q2 -0.017503 0.002355 -7.431 6.72e-10 ***
## Q3 -0.022966 0.002400 -9.567 2.22e-13 ***
## Q4 -0.006173 0.002395 -2.578 0.01261 *
## post -0.019767 0.002775 -7.123 2.17e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006649 on 56 degrees of freedom
## Multiple R-squared: 0.9987, Adjusted R-squared: 0.9986
## F-statistic: 8939 on 5 and 56 DF, p-value: < 2.2e-16
beta_lr <- coef(model_lr)["lnY"]
ect <- residuals(model_lr) # ECT_t = ln C - ln C keseimbangan
round(100 * (exp(coef(model_lr)["post"]) - 1), 2) # % pergeseran level konsumsi pasca-COVID
## post
## -1.96
ect_ts <- ts(100 * ect, start = c(df$Tahun[1], df$Q[1]), frequency = 4)
plot(ect_ts, type = "l", lwd = 2, col = "navy", xlab = "",
ylab = "% di atas/bawah keseimbangan",
main = "Error correction term: deviasi konsumsi dari keseimbangan jangka panjang")
abline(h = 0, lty = 2, col = "grey40")
# Sampel disamakan untuk semua model: mulai 2012Q2 (karena ada lag musiman t-4)
dC <- diff(lnC); dYv <- diff(lnY) # dC[k] = lnC[k+1] - lnC[k]
t_idx <- 6:n
ecm_df <- data.frame(
Triwulan = df$Triwulan[t_idx],
dlnC = dC[t_idx - 1], dlnY = dYv[t_idx - 1],
dlnC_1 = dC[t_idx - 2], dlnY_1 = dYv[t_idx - 2],
dlnC_4 = dC[t_idx - 5], dlnY_4 = dYv[t_idx - 5],
ect1 = ect[t_idx - 1], # ECT periode t-1
Q2 = Q2[t_idx], Q3 = Q3[t_idx], Q4 = Q4[t_idx],
pQ2 = pQ2[t_idx], pQ3 = pQ3[t_idx], pQ4 = pQ4[t_idx],
covid_i = covid_i[t_idx])
nrow(ecm_df)
## [1] 57
spes_sr <- list(
S1 = dlnC ~ dlnY + ect1,
S2 = dlnC ~ dlnY + ect1 + Q2 + Q3 + Q4,
S3 = dlnC ~ dlnY + ect1 + Q2 + Q3 + Q4 + pQ2 + pQ3 + pQ4,
S4 = dlnC ~ dlnY + ect1 + Q2 + Q3 + Q4 + pQ2 + pQ3 + pQ4 + dlnC_4 + dlnY_4,
S5 = dlnC ~ dlnY + ect1 + Q2 + Q3 + Q4 + pQ2 + pQ3 + pQ4 + dlnC_4 + dlnY_4 + covid_i,
S6 = dlnC ~ dlnY + ect1 + Q2 + Q3 + Q4 + pQ2 + pQ3 + pQ4 + dlnC_1 + dlnY_1 + dlnC_4 + dlnY_4 + covid_i)
ket_sr <- c(S1 = "dasar", S2 = "+ musiman", S3 = "+ musiman pasca-COVID",
S4 = "S3 + lag musiman (t-4)", S5 = "S4 + dummy COVID 2020Q2",
S6 = "S5 + lag 1")
hasil_sr <- do.call(rbind, lapply(names(spes_sr), function(nm) {
m <- lm(spes_sr[[nm]], data = ecm_df)
ct <- summary(m)$coefficients
data.frame(Model = nm, Spesifikasi = ket_sr[nm],
lambda = round(ct["ect1", 1], 4), t_lambda = round(ct["ect1", 3], 2),
p_lambda = round(ct["ect1", 4], 4),
gamma_dlnY = round(ct["dlnY", 1], 4),
adjR2 = round(summary(m)$adj.r.squared, 4),
BIC = round(BIC(m), 2),
BG4_p = round(bgtest(m, order = 4)$p.value, 4),
BP_p = round(bptest(m)$p.value, 4),
JB_p = round(jarque.bera.test(residuals(m))$p.value, 4))
}))
rownames(hasil_sr) <- NULL
hasil_sr$Layak <- with(hasil_sr, lambda < 0 & p_lambda < 0.05 & BG4_p > 0.05 & BP_p > 0.05 & JB_p > 0.05)
hasil_sr
## Model Spesifikasi lambda t_lambda p_lambda gamma_dlnY adjR2
## 1 S1 dasar -0.8859 -4.00 2e-04 0.5457 0.6432
## 2 S2 + musiman -0.6819 -5.25 0e+00 1.0127 0.8812
## 3 S3 + musiman pasca-COVID -0.4747 -5.50 0e+00 1.1707 0.9645
## 4 S4 S3 + lag musiman (t-4) -0.3085 -4.83 0e+00 1.1103 0.9837
## 5 S5 S4 + dummy COVID 2020Q2 -0.3242 -5.13 0e+00 1.1954 0.9843
## 6 S6 S5 + lag 1 -0.4049 -4.98 0e+00 1.1767 0.9846
## BIC BG4_p BP_p JB_p Layak
## 1 -342.55 0.0000 0.1359 0.0000 FALSE
## 2 -396.36 0.0000 0.0035 0.3721 FALSE
## 3 -456.60 0.0054 0.0120 0.0000 FALSE
## 4 -495.19 0.0606 0.0896 0.0005 FALSE
## 5 -494.74 0.1008 0.1023 0.1012 TRUE
## 6 -490.16 0.5402 0.0929 0.1068 TRUE
# Layak: lambda negatif & signifikan, lolos uji autokorelasi (BG lag 4),
# heteroskedastisitas (BP), dan normalitas (JB). Dipilih: layak dengan BIC terkecil.
kand_sr <- hasil_sr[hasil_sr$Layak, ]
pilih_sr <- kand_sr$Model[which.min(kand_sr$BIC)]
pilih_sr
## [1] "S5"
ecm_model <- lm(spes_sr[[pilih_sr]], data = ecm_df)
summary(ecm_model)
##
## Call:
## lm(formula = spes_sr[[pilih_sr]], data = ecm_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0045820 -0.0011578 0.0001019 0.0009233 0.0065905
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0052211 0.0008054 6.483 5.98e-08 ***
## dlnY 1.1954490 0.0568770 21.018 < 2e-16 ***
## ect1 -0.3241555 0.0632456 -5.125 6.05e-06 ***
## Q2 -0.0226658 0.0035830 -6.326 1.02e-07 ***
## Q3 -0.0107636 0.0024242 -4.440 5.78e-05 ***
## Q4 0.0068673 0.0016594 4.138 0.000151 ***
## pQ2 0.0075936 0.0015383 4.937 1.14e-05 ***
## pQ3 -0.0043811 0.0017535 -2.498 0.016195 *
## pQ4 -0.0084894 0.0019060 -4.454 5.52e-05 ***
## dlnC_4 0.5281887 0.0665873 7.932 4.33e-10 ***
## dlnY_4 -0.5818918 0.0785679 -7.406 2.56e-09 ***
## covid_i 0.0091706 0.0053575 1.712 0.093830 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002239 on 45 degrees of freedom
## Multiple R-squared: 0.9874, Adjusted R-squared: 0.9843
## F-statistic: 320.9 on 11 and 45 DF, p-value: < 2.2e-16
bgtest(ecm_model, order = 4) # autokorelasi
##
## Breusch-Godfrey test for serial correlation of order up to 4
##
## data: ecm_model
## LM test = 7.7595, df = 4, p-value = 0.1008
bptest(ecm_model) # heteroskedastisitas
##
## studentized Breusch-Pagan test
##
## data: ecm_model
## BP = 17.193, df = 11, p-value = 0.1023
jarque.bera.test(residuals(ecm_model)) # normalitas
##
## Jarque Bera Test
##
## data: residuals(ecm_model)
## X-squared = 4.5811, df = 2, p-value = 0.1012
fit_ts <- ts(cbind(Aktual = ecm_df$dlnC, Fitted = fitted(ecm_model)) * 100,
start = c(2012, 2), frequency = 4)
plot(fit_ts, plot.type = "single", col = c("black", "red"), lty = c(1, 2), lwd = 2,
ylab = "% per triwulan", xlab = "", main = "Pertumbuhan konsumsi (q-to-q): aktual vs fitted")
legend("topleft", c("Aktual", "Fitted"), col = c("black", "red"), lty = 1:2, lwd = 2, bty = "n")
lambda <- coef(ecm_model)["ect1"]
gamma0 <- coef(ecm_model)["dlnY"]
porsi_C <- mean(df$Konsumsi_RT / df$PDB)
ringkas <- c(
elastisitas_jangka_panjang = unname(beta_lr),
elastisitas_jangka_pendek = unname(gamma0),
lambda_ECT = unname(lambda),
persen_koreksi_per_triwulan = unname(-100 * lambda),
half_life_triwulan = unname(log(0.5) / log(1 + lambda)),
MPC_jangka_panjang = unname(beta_lr * porsi_C))
round(ringkas, 4)
## elastisitas_jangka_panjang elastisitas_jangka_pendek
## 0.9792 1.1954
## lambda_ECT persen_koreksi_per_triwulan
## -0.3242 32.4155
## half_life_triwulan MPC_jangka_panjang
## 1.7692 0.5264
# ECT pada periode penting (% deviasi dari keseimbangan)
penting <- c("2019Q4", "2020Q1", "2020Q2", "2020Q3", "2020Q4", "2021Q1", "2021Q2",
"2021Q3", "2021Q4", "2022Q4", tail(df$Triwulan, 1))
round(setNames(100 * ect[match(penting, df$Triwulan)], penting), 2)
## 2019Q4 2020Q1 2020Q2 2020Q3 2020Q4 2021Q1 2021Q2 2021Q3 2021Q4 2022Q4 2026Q2
## 0.48 0.24 1.41 1.72 0.93 0.65 0.50 -0.66 -0.39 -0.77 0.66