## ============================================================
## STUDI KASUS 3 — ECM: Konsumsi Rumah Tangga dan PDB
## UAS Analisis Runtun Waktu Terapan
## ============================================================
## Metode utama: Engle-Granger dua tahap
## Tahap 1: regresi jangka panjang ln C terhadap ln Y + uji kointegrasi
## Tahap 2: ECM jangka pendek dengan ECT(t-1)
## Uji kointegrasi pembanding: uji-t ECM (Banerjee-Dolado-Mestre) & ARDL bounds.
## Data: data_pdb.csv (BPS, PDB triwulanan harga konstan 2010, miliar Rp)
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
library(sandwich)
library(strucchange)
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("firebrick", "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("firebrick", "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, pQ2, pQ3, pQ4)
# Nilai kritis Engle-Granger (MacKinnon 2010, 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)
}
# Dengan pergeseran level: nilai kritis Gregory-Hansen (1996), model C, m = 1 (lebih ketat)
cv_gh <- c(`1%` = -5.13, `5%` = -4.61, `10%` = -4.34)
spes_lr <- list(
LR1 = list(f = lnC ~ lnY, cv = cv_eg(n), ket = "tanpa dummy"),
LR2 = list(f = lnC ~ lnY + Q2 + Q3 + Q4, cv = cv_eg(n), ket = "+ dummy musiman"),
LR3 = list(f = lnC ~ lnY + Q2 + Q3 + Q4 + post, cv = cv_gh, ket = "+ musiman + level shift pasca-COVID"),
LR4 = list(f = lnC ~ lnY + Q2 + Q3 + Q4 + post + pQ2 + pQ3 + pQ4, cv = cv_gh, ket = "LR3 + musiman pasca-COVID"))
tau_resid <- function(u, sel) ur.df(u, type = "none", lags = 4, selectlags = sel)@teststat[1]
hasil_lr <- do.call(rbind, lapply(names(spes_lr), function(nm) {
s <- spes_lr[[nm]]
m <- lm(s$f, data = dl); u <- residuals(m)
t_bic <- tau_resid(u, "BIC"); t_aic <- tau_resid(u, "AIC")
data.frame(Model = nm, Spesifikasi = s$ket,
beta_lnY = round(coef(m)["lnY"], 4),
sd_resid_pct = round(100 * sd(u), 3),
BIC = round(BIC(m), 2),
tau_BIC = round(t_bic, 3), tau_AIC = round(t_aic, 3),
CV_5pct = s$cv["5%"], CV_1pct = s$cv["1%"],
Kointegrasi_BIC = t_bic < s$cv["5%"], Kointegrasi_AIC = t_aic < s$cv["5%"])
}))
rownames(hasil_lr) <- NULL
hasil_lr
## Model Spesifikasi beta_lnY sd_resid_pct BIC
## 1 LR1 tanpa dummy 0.9361 1.250 -356.07
## 2 LR2 + dummy musiman 0.9392 0.880 -387.25
## 3 LR3 + musiman + level shift pasca-COVID 0.9792 0.637 -423.12
## 4 LR4 LR3 + musiman pasca-COVID 0.9780 0.437 -457.38
## tau_BIC tau_AIC CV_5pct CV_1pct Kointegrasi_BIC Kointegrasi_AIC
## 1 -1.979 -1.979 -3.436 -4.079 FALSE FALSE
## 2 -1.867 -1.867 -3.436 -4.079 FALSE FALSE
## 3 -2.639 -2.639 -4.610 -5.130 FALSE FALSE
## 4 -1.804 -2.132 -4.610 -5.130 FALSE FALSE
# Lag ADF residual dipilih dengan BIC (konsisten, menghindari lag berlebih);
# AIC ditampilkan sebagai pembanding kepekaan.
pilih_lr <- "LR3" # terkointegrasi (BIC-lag) dengan nilai kritis Gregory-Hansen yang ketat
model_lr <- lm(spes_lr[[pilih_lr]]$f, data = dl)
summary(model_lr)
##
## Call:
## lm(formula = spes_lr[[pilih_lr]]$f, 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
# Elastisitas jangka panjang dengan Dynamic OLS (standar error valid, HAC)
dY <- c(NA, diff(lnY))
dols_df <- cbind(dl, dY_lead1 = c(dY[-1], NA), dY0 = dY, dY_lag1 = c(NA, dY[-n]))
model_dols <- lm(update(spes_lr[[pilih_lr]]$f, . ~ . + dY_lead1 + dY0 + dY_lag1), data = dols_df)
V_d <- NeweyWest(model_dols, lag = 4, prewhite = FALSE)
b_d <- coef(model_dols)["lnY"]; se_d <- sqrt(V_d["lnY", "lnY"])
round(c(beta_DOLS = b_d, se_HAC = se_d, t_beta_sama_1 = (b_d - 1) / se_d,
p_beta_sama_1 = 2 * pt(-abs((b_d - 1) / se_d), df = model_dols$df.residual)), 4)
## beta_DOLS.lnY se_HAC t_beta_sama_1.lnY p_beta_sama_1.lnY
## 0.9746 0.0076 -3.3280 0.0016
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),
AIC = round(AIC(m), 2), 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),
RESET_p = round(resettest(m, power = 2:3)$p.value, 4))
}))
rownames(hasil_sr) <- NULL
hasil_sr$Layak <- with(hasil_sr, lambda < 0 & lambda > -2 & p_lambda < 0.05 & BG4_p > 0.05)
hasil_sr$Lolos_semua <- with(hasil_sr, Layak & BP_p > 0.05 & JB_p > 0.05 & RESET_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
## AIC BIC BG4_p BP_p JB_p RESET_p Layak Lolos_semua
## 1 -350.72 -342.55 0.0000 0.1359 0.0000 0.0046 FALSE FALSE
## 2 -410.66 -396.36 0.0000 0.0035 0.3721 0.0246 FALSE FALSE
## 3 -477.03 -456.60 0.0054 0.0120 0.0000 0.0039 FALSE FALSE
## 4 -519.70 -495.19 0.0606 0.0896 0.0005 0.0124 TRUE FALSE
## 5 -521.30 -494.74 0.1008 0.1023 0.1012 0.0212 TRUE FALSE
## 6 -520.81 -490.16 0.5402 0.0929 0.1068 0.0446 TRUE FALSE
# Layak: lambda negatif & signifikan, -2 < lambda < 0, bebas autokorelasi (BG lag 4).
# Dipilih: lolos semua diagnostik dengan BIC terkecil (jika tidak ada: layak dengan BIC terkecil).
kand_sr <- hasil_sr[hasil_sr$Lolos_semua, ]
if (nrow(kand_sr) == 0) kand_sr <- hasil_sr[hasil_sr$Layak, ]
if (nrow(kand_sr) == 0) kand_sr <- hasil_sr
pilih_sr <- kand_sr$Model[which.min(kand_sr$BIC)]
pilih_sr # <-- boleh diganti manual
## [1] "S4"
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.0052447 -0.0011155 0.0001525 0.0007488 0.0071537
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0048569 0.0007929 6.126 1.88e-07 ***
## dlnY 1.1102956 0.0281452 39.449 < 2e-16 ***
## ect1 -0.3084738 0.0638779 -4.829 1.56e-05 ***
## Q2 -0.0200528 0.0033089 -6.060 2.35e-07 ***
## Q3 -0.0082154 0.0019530 -4.207 0.000119 ***
## Q4 0.0062197 0.0016493 3.771 0.000462 ***
## pQ2 0.0081094 0.0015398 5.267 3.57e-06 ***
## pQ3 -0.0055309 0.0016534 -3.345 0.001644 **
## pQ4 -0.0068545 0.0016837 -4.071 0.000182 ***
## dlnC_4 0.5151890 0.0675263 7.629 1.05e-09 ***
## dlnY_4 -0.5492042 0.0777943 -7.060 7.42e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002286 on 46 degrees of freedom
## Multiple R-squared: 0.9866, Adjusted R-squared: 0.9837
## F-statistic: 338.5 on 10 and 46 DF, p-value: < 2.2e-16
coeftest(ecm_model, vcov = NeweyWest(ecm_model, lag = 4, prewhite = FALSE)) # SE robust
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0048569 0.0009457 5.1358 5.565e-06 ***
## dlnY 1.1102956 0.0245653 45.1978 < 2.2e-16 ***
## ect1 -0.3084738 0.0724130 -4.2599 0.0001000 ***
## Q2 -0.0200528 0.0041746 -4.8035 1.694e-05 ***
## Q3 -0.0082154 0.0019782 -4.1530 0.0001407 ***
## Q4 0.0062197 0.0016383 3.7964 0.0004282 ***
## pQ2 0.0081094 0.0015268 5.3113 3.071e-06 ***
## pQ3 -0.0055309 0.0017273 -3.2021 0.0024764 **
## pQ4 -0.0068545 0.0012132 -5.6500 9.659e-07 ***
## dlnC_4 0.5151890 0.1011756 5.0920 6.451e-06 ***
## dlnY_4 -0.5492042 0.1018558 -5.3920 2.334e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji kointegrasi berbasis ECM (Banerjee, Dolado & Mestre, 1998):
# H0 tidak ada kointegrasi ditolak jika t(lambda) < nilai kritis
# (k = 1 regresor, dengan konstanta, T = 50: 1% = -3,94; 5% = -3,28; 10% = -2,93)
round(coef(summary(ecm_model))["ect1", "t value"], 3)
## [1] -4.829
bgtest(ecm_model, order = 4)
##
## Breusch-Godfrey test for serial correlation of order up to 4
##
## data: ecm_model
## LM test = 9.0212, df = 4, p-value = 0.06057
bptest(ecm_model)
##
## studentized Breusch-Pagan test
##
## data: ecm_model
## BP = 16.368, df = 10, p-value = 0.08957
jarque.bera.test(residuals(ecm_model))
##
## Jarque Bera Test
##
## data: residuals(ecm_model)
## X-squared = 15.294, df = 2, p-value = 0.0004775
resettest(ecm_model, power = 2:3)
##
## RESET test
##
## data: ecm_model
## RESET = 4.8609, df1 = 2, df2 = 44, p-value = 0.01238
cusum <- efp(spes_sr[[pilih_sr]], data = ecm_df, type = "OLS-CUSUM")
sctest(cusum)
##
## OLS-based CUSUM test
##
## data: cusum
## S0 = 0.52709, p-value = 0.9439
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
plot(cusum, main = "OLS-CUSUM model ECM (batas 5%)")
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", "firebrick"), lty = c(1, 2), lwd = 2,
ylab = "% per triwulan", xlab = "", main = "Pertumbuhan konsumsi (q-to-q): aktual vs fitted")
legend("bottomleft", c("Aktual", "Fitted"), col = c("black", "firebrick"), lty = 1:2, lwd = 2, bty = "n")
par(mfrow = c(1, 1))
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_DOLS = unname(b_d),
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),
MPC_jangka_pendek = unname(gamma0 * porsi_C))
round(ringkas, 4)
## elastisitas_jangka_panjang elastisitas_DOLS
## 0.9792 0.9746
## elastisitas_jangka_pendek lambda_ECT
## 1.1103 -0.3085
## persen_koreksi_per_triwulan half_life_triwulan
## 30.8474 1.8792
## MPC_jangka_panjang MPC_jangka_pendek
## 0.5264 0.5969
# 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
# Komponen deterministik sama dengan model utama (musiman, pasca-COVID, musiman pasca-COVID).
if (requireNamespace("ARDL", quietly = TRUE)) {
library(ARDL)
ardl_df <- data.frame(lnC, lnY, Q2, Q3, Q4, post, pQ2, pQ3, pQ4)
cari <- auto_ardl(lnC ~ lnY | Q2 + Q3 + Q4 + post + pQ2 + pQ3 + pQ4,
data = ardl_df, max_order = 4, selection = "BIC")
print(cari$best_order)
ardl_best <- cari$best_model
print(bounds_f_test(ardl_best, case = 3))
print(multipliers(ardl_best)) # elastisitas jangka panjang versi ARDL
print(summary(recm(ardl_best, case = 3))) # ECM versi ARDL (koefisien ect)
} else {
cat("Package ARDL belum terpasang: install.packages('ARDL') untuk uji bounds.\n")
}
## To cite the ARDL package in publications:
##
## Use this reference to refer to the validity of the ARDL package.
##
## Natsiopoulos, Kleanthis, and Tzeremes, Nickolaos G. (2022). ARDL
## bounds test for cointegration: Replicating the Pesaran et al. (2001)
## results for the UK earnings equation using R. Journal of Applied
## Econometrics, 37(5), 1079-1090. https://doi.org/10.1002/jae.2919
##
## Use this reference to cite this specific version of the ARDL package.
##
## Kleanthis Natsiopoulos, Nickolaos Tzeremes and Daniel Finnan (2026).
## ARDL: ARDL, ECM and Bounds-Test for Cointegration. R package version
## 0.2.5. https://CRAN.R-project.org/package=ARDL
## lnC lnY
## 4 4
##
## Bounds F-test (Wald) for no cointegration
##
## data: d(lnC) ~ L(lnC, 1) + L(lnY, 1) + d(L(lnC, 1)) + d(L(lnC, 2)) + d(L(lnC, 3)) + d(lnY) + d(L(lnY, 1)) + d(L(lnY, 2)) + d(L(lnY, 3)) + Q2 + Q3 + Q4 + post + pQ2 + pQ3 + pQ4
## F = 57.389, p-value = 1e-06
## alternative hypothesis: Possible cointegration
## null values:
## k T
## 1 1000
##
## Term Estimate Std. Error t value Pr(>|t|)
## 1 (Intercept) -0.4911415 0.04577066 -10.73049 1.779719e-13
## 2 lnY 0.9925231 0.00308150 322.09085 2.110665e-71
##
## Time series regression with "zooreg" data:
## Start = 5, End = 62
##
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start,
## end = end)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0040222 -0.0006789 0.0002306 0.0009082 0.0028270
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3616369 0.0321685 -11.242 3.04e-14 ***
## d(L(lnC, 1)) -0.1205509 0.0652438 -1.848 0.071697 .
## d(L(lnC, 2)) -0.1479191 0.0644767 -2.294 0.026850 *
## d(L(lnC, 3)) -0.4003809 0.0629973 -6.356 1.22e-07 ***
## d(lnY) 0.9586966 0.0251353 38.142 < 2e-16 ***
## d(L(lnY, 1)) -0.1082965 0.0821727 -1.318 0.194676
## d(L(lnY, 2)) -0.0700864 0.0837690 -0.837 0.407518
## d(L(lnY, 3)) 0.3348965 0.0840043 3.987 0.000262 ***
## Q2 -0.0199389 0.0037688 -5.291 4.12e-06 ***
## Q3 -0.0002477 0.0064653 -0.038 0.969622
## Q4 0.0169155 0.0043294 3.907 0.000334 ***
## post -0.0124282 0.0033125 -3.752 0.000532 ***
## pQ2 0.0035685 0.0032979 1.082 0.285403
## pQ3 -0.0176438 0.0040190 -4.390 7.50e-05 ***
## pQ4 -0.0171383 0.0027205 -6.300 1.47e-07 ***
## ect -0.7363192 0.0679052 -10.843 9.47e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001498 on 42 degrees of freedom
## (0 observations deleted due to missingness)
## Multiple R-squared: 0.9948, Adjusted R-squared: 0.9929
## F-statistic: 532.4 on 15 and 42 DF, p-value: < 2.2e-16