Analisis Data Panel: CEM dan FEM
Abdullah Zeith Alfarizy (G1401231060)
2026-03-03
1 PACKAGES
pkgs <- c("plm","lmtest","car","ggplot2","readxl","dplyr","sandwich","moments","knitr","kableExtra")
for (p in pkgs) if (!requireNamespace(p, quietly=TRUE)) install.packages(p, quiet=TRUE)
library(readxl); library(dplyr); library(ggplot2); library(plm)
library(lmtest); library(car); library(sandwich)
library(moments); library(knitr); library(kableExtra)
cetak_koef <- function(obj, judul = "") {
cf <- if (is.matrix(obj) || is.data.frame(obj)) as.matrix(obj)
else as.matrix(coef(summary(obj)))
nc <- ncol(cf)
if (nchar(judul) > 0) cat(sprintf("\n=== %s ===\n", judul))
cat(sprintf("%-20s %10s %10s %8s %9s\n","Variabel","Koef","SE","t","p-value"))
cat(strrep("-", 62), "\n")
for (i in seq_len(nrow(cf))) {
pv <- as.numeric(cf[i, nc])
sig <- ifelse(pv<0.001,"***", ifelse(pv<0.01,"**", ifelse(pv<0.05,"*", ifelse(pv<0.10,".",""))))
cat(sprintf("%-20s %10.4f %10.4f %8.3f %9.4f %s\n",
rownames(cf)[i], as.numeric(cf[i,1]),
as.numeric(cf[i,2]), as.numeric(cf[i,3]), pv, sig))
}
cat("Signif: *** <.001 ** <.01 * <.05 . <.10\n")
}2 DATA DAN DEKLARASI PANEL
#Input data
mikro <- read.csv(
"C:\\Users\\FAQIH\\Downloads\\data_mikro_adp.csv",
stringsAsFactors = FALSE
)
colnames(mikro) <- c("kode","tahun","laba","total_aset","der","cr")
makro_raw <- read_excel("C:\\Users\\FAQIH\\Downloads\\Data panel kel 7.xlsx",
col_names = FALSE)
colnames(makro_raw) <- as.character(makro_raw[1, ])
makro_raw <- makro_raw[-1, ]
colnames(makro_raw) <- c("perusahaan","tahun","laba_m","pdb_re","pengeluaran",
"ihpb","ekspor","impor","inflasi","bi_rate",
"m2","kurs","pdb_re_growth")
makro <- makro_raw %>%
mutate(across(c(tahun,laba_m,pdb_re,pengeluaran,ihpb,ekspor,impor,
inflasi,bi_rate,m2,kurs,pdb_re_growth), as.numeric))
#Seleksi variabel makro (3 kriteria)
makro_pilih <- makro %>%
dplyr::select(perusahaan, tahun, inflasi, bi_rate, ihpb) %>% distinct()
data_panel <- mikro %>%
left_join(makro_pilih, by = c("kode"="perusahaan","tahun"="tahun"))
cat("Dimensi data:", nrow(data_panel), "×", ncol(data_panel), "\n")## Dimensi data: 140 × 9## Missing:## kode tahun laba total_aset der cr inflasi
## 0 0 0 0 0 0 0
## bi_rate ihpb
## 0 0#Transformasi (mengatasi non-normalitas ekstrem pada laba)
signed_log <- function(x) sign(x) * log1p(abs(x))
data_panel <- data_panel %>%
mutate(
Y = signed_log(laba),
aset_log = log(total_aset),
inflasi_z = as.numeric(scale(inflasi)),
birate_z = as.numeric(scale(bi_rate)),
ihpb_z = as.numeric(scale(ihpb))
)
#Deklarasi panel
pdata <- pdata.frame(data_panel, index = c("kode","tahun"), drop.index = FALSE)
cat("\nBalanced panel:", is.pbalanced(pdata), "\n")##
## Balanced panel: TRUE## Balanced Panel: n = 14, T = 10, N = 140N_ <- length(unique(data_panel$kode)) # 14
T_ <- length(unique(data_panel$tahun)) # 10
NT_ <- nrow(data_panel) # 140
K_ <- 6 # aset_log, der, cr, inflasi_z, birate_z, ihpb_z
xvars <- c("aset_log","der","cr","inflasi_z","birate_z","ihpb_z")
cat(sprintf("\nN=%d | T=%d | NT=%d | K=%d\n", N_, T_, NT_, K_))##
## N=14 | T=10 | NT=140 | K=6##
## Distribusi signed_log(laba):## Skewness = -2.575## Kurtosis = 8.467sw_y <- shapiro.test(data_panel$Y)
cat(sprintf(" SW W=%.4f, p=%.4f\n", sw_y$statistic, sw_y$p.value))## SW W=0.5789, p=0.00003 BAGIAN I: COMMON EFFECT MODEL (CEM)
Model: \(y_{it} = x'_{it}\beta + u + e_{it}\), asumsi \(u_i = u\) untuk semua \(i\)
3.1 CEM-1: SPESIFIKASI DAN ESTIMASI OLS
\[y_{it} = w'_{it}\delta + e_{it}, \quad w_{it} = [x'_{it}\;\; 1]', \quad \delta = [\beta'\;\; u]'\]
\[\hat{\delta}_{OLS} = (W'W)^{-1}W'y\]
\[\hat{\sigma}^2_e = \frac{\hat{e}'\hat{e}}{NT - K}, \quad K = \text{ncol}(W) = 7\]
## CEM-1: COMMON EFFECT MODEL — OLS## y_it = x'_it β + u + e_it [u_i = u untuk semua i]## δ̂ = (W'W)⁻¹W'y | σ²_e = ê'ê/(NT-K)cem <- plm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z,
data = pdata, model = "pooling", effect = "individual")
resid_cem <- as.numeric(residuals(cem))
# Slide 41: σ̂²_e = ê'ê/(NT-K), K = ncol(W) = 7 (6 prediktor + 1 intercept)
sigma2_cem <- sum(resid_cem^2) / (NT_ - K_ - 1)
sigma_cem <- sqrt(sigma2_cem)
cetak_koef(cem, "CEM — Pooled OLS")##
## === CEM — Pooled OLS ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 2.7462 4.0237 0.682 0.4961
## aset_log 0.4795 0.4232 1.133 0.2593
## der -1.4322 0.7781 -1.841 0.0679 .
## cr -0.3055 0.1327 -2.302 0.0229 *
## inflasi_z 0.7920 0.3468 2.283 0.0240 *
## birate_z 0.2875 0.3793 0.758 0.4498
## ihpb_z 1.0405 0.3406 3.055 0.0027 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## σ²_e = ê'ê/(NT-K) = 1687.5502 / 133 = 12.688347## σ_e = 3.5621s_cem <- summary(cem)
cat(sprintf("R² = %.4f | Adj.R² = %.4f\n",
s_cem$r.squared["rsq"], s_cem$r.squared["adjrsq"]))## R² = 0.2127 | Adj.R² = 0.17723.2 CEM-2: ASUMSI KLASIK CEM
| Asumsi | Notasi |
|---|---|
| Strict Exogeneity | \(E(e_{it}\|W) = 0\) |
| Homoskedastisitas | \(\text{Var}(e_{it}\|W) = \sigma^2_e\), \(\text{Cov}(w_{it}, e_{it}) = 0\) |
| No Serial & CS Corr | \(\text{Var}(e\|W) = \sigma^2_e I_{NT}\) |
## CEM-2: UJI ASUMSI KLASIK## [A1] STRICT EXOGENEITY: E(e_it|W) = 0## Terpenuhi by construction OLS.## Jika efek individu berkorelasi dengan X → CEM BIAS# [A2] Homoskedastisitas — Slide 52: uji nR² > χ²(1,α)
cat("[A2] HOMOSKEDASTISITAS: Var(e_it|W) = σ²_e\n")## [A2] HOMOSKEDASTISITAS: Var(e_it|W) = σ²_e## Slide 52: H0 tolak jika nR² > χ²(1,α)bp_cem <- bptest(cem)
cat(sprintf(" nR² = %.4f, df=%d, p = %.6f → %s\n\n",
bp_cem$statistic, bp_cem$parameter, bp_cem$p.value,
ifelse(bp_cem$p.value < 0.05, "DILANGGAR ⚠ (heteroskedastisitas)", "TERPENUHI ✓")))## nR² = 24.0467, df=6, p = 0.000512 → DILANGGAR ⚠ (heteroskedastisitas)# [A3] No Serial Correlation — DW Slide 55-56
cat("[A3] NO SERIAL CORRELATION: Cov(e_it, e_is|W) = 0\n")## [A3] NO SERIAL CORRELATION: Cov(e_it, e_is|W) = 0## Durbin-Watson (Slide 55): d ≈ 2(1-r), tolak H0 jika d < d_L# Tabel DW, n=140, K=6, α=5%: d_L=1.662, d_U=1.802
dL <- 1.662; dU <- 1.802
cat(sprintf(" d_L = %.3f, d_U = %.3f (tabel DW, n=%d, K=%d, α=5%%)\n",
dL, dU, NT_, K_))## d_L = 1.662, d_U = 1.802 (tabel DW, n=140, K=6, α=5%)dw_cem <- pdwtest(cem)
dw_ket_cem <- ifelse(dw_cem$statistic < dL, "AUTOKORELASI POSITIF ⚠",
ifelse(dw_cem$statistic <= dU, "INCONCLUSIVE",
ifelse(dw_cem$statistic < 4-dU, "TIDAK ADA autokorelasi ✓",
ifelse(dw_cem$statistic < 4-dL, "INCONCLUSIVE",
"AUTOKORELASI NEGATIF ⚠"))))
cat(sprintf(" DW = %.4f → %s\n\n", dw_cem$statistic, dw_ket_cem))## DW = 1.3111 → AUTOKORELASI POSITIF ⚠# [A4] Normalitas (untuk inferensi hipotesis) — Shapiro-Wilk
cat("[A4] NORMALITAS: e_it ~ N(0,σ²_e) (untuk inferensi hipotesis)\n")## [A4] NORMALITAS: e_it ~ N(0,σ²_e) (untuk inferensi hipotesis)sw_cem <- shapiro.test(resid_cem)
cat(sprintf(" SW W = %.4f, p = %.2e → %s\n\n",
sw_cem$statistic, sw_cem$p.value,
ifelse(sw_cem$p.value > 0.05, "TERPENUHI ✓", "DILANGGAR ⚠")))## SW W = 0.7813, p = 3.54e-13 → DILANGGAR ⚠## [A5] NO MULTIKOLINEARITAS: rank(X) = Klm_vif <- lm(Y ~ aset_log+der+cr+inflasi_z+birate_z+ihpb_z, data=data_panel)
vif_vals <- vif(lm_vif)
cat(" VIF per variabel:\n")## VIF per variabel:## aset_log der cr inflasi_z birate_z ihpb_z
## 1.227 1.385 1.251 1.318 1.576 1.271## → Semua VIF < 10 ✓# Ringkasan
ring_cem <- data.frame(
Asumsi = c("A1 Exogeneity","A2 Homoskedastisitas","A3 No Serial Corr (DW)",
"A4 Normalitas","A5 No Multikolinearitas"),
Status = c(
"Perlu validasi (F-test / FEM)",
ifelse(bp_cem$p.value < 0.05, "DILANGGAR ⚠", "TERPENUHI ✓"),
ifelse(grepl("TIDAK", dw_ket_cem), "TERPENUHI ✓",
ifelse(grepl("INCON", dw_ket_cem), "INCONCLUSIVE", "DILANGGAR ⚠")),
ifelse(sw_cem$p.value > 0.05, "TERPENUHI ✓", "DILANGGAR ⚠"),
ifelse(all(vif_vals < 10), "TERPENUHI ✓", "MASALAH ⚠"))
)
knitr::kable(ring_cem, caption = "Ringkasan Asumsi Klasik CEM") %>%
kableExtra::kable_styling(bootstrap_options = c("striped","hover","condensed"))| Asumsi | Status |
|---|---|
| A1 Exogeneity | Perlu validasi (F-test / FEM) |
| A2 Homoskedastisitas | DILANGGAR ⚠ |
| A3 No Serial Corr (DW) | DILANGGAR ⚠ |
| A4 Normalitas | DILANGGAR ⚠ |
| A5 No Multikolinearitas | TERPENUHI ✓ |
3.3 CEM-3: ANALISIS RESIDUAL
Plot residual: vs \(\hat{y}_i\), vs \(x_1\), vs \(x_2\), vs waktu, QQ normal.
## CEM-3: ANALISIS RESIDUAL (Slide 47-50)yhat_cem <- as.numeric(fitted(cem))
par(mfrow = c(2, 4), mar = c(4, 4, 3, 1))
# Residual vs ŷ
plot(yhat_cem, resid_cem, pch=16, cex=0.7, col="#e63946",
main="Resid vs Fitted\n[Linearitas & Homosk]",
xlab="ŷ", ylab="Residual")
abline(h=0, col="red", lty=2)
lines(lowess(yhat_cem, resid_cem), col="orange", lwd=2)
# Residual vs X1..X6
nms_x <- c("log(Aset)","DER","Current Ratio","Inflasi(z)","BI Rate(z)","IHPB(z)")
for (k in 1:6) {
plot(data_panel[[xvars[k]]], resid_cem, pch=16, cex=0.7, col="#e63946",
main=paste("Resid vs", nms_x[k]),
xlab=nms_x[k], ylab="Residual")
abline(h=0, col="red", lty=2)
}
# Residual vs Waktu
plot(data_panel$tahun, resid_cem, pch=16, cex=0.7, col="#e63946",
main="Resid vs Waktu\n[Serial Correlation]",
xlab="Tahun", ylab="Residual")
abline(h=0, col="red", lty=2)
par(mfrow=c(1,2), mar=c(4,4,3,1))
# QQ Normal
qqnorm(resid_cem, pch=16, cex=0.7, col="#e63946", main="Normal Q-Q Plot (CEM)")
qqline(resid_cem, col="red", lty=2)
# Histogram residual
hist(resid_cem, breaks=25, col="#e63946", border="white",
main="Distribusi Residual CEM", xlab="Residual", freq=FALSE)
curve(dnorm(x, mean=mean(resid_cem), sd=sd(resid_cem)),
add=TRUE, col="black", lwd=2)
3.4 CEM-4: PENANGANAN PELANGGARAN — PANEL-ROBUST VCV
Jika asumsi dilanggar (heteroskedastisitas dan/atau serial korelasi):
\[\widehat{Var}(\hat{\delta}_{OLS}) = (W'W)^{-1}\left[\sum_i W'_i\hat{e}_i\hat{e}'_iW_i\right](W'W)^{-1}\]
## CEM-4: PANEL-ROBUST VCV (Slide 42)## V̂ar(δ̂) = (W'W)⁻¹[Σᵢ W'ᵢêᵢê'ᵢWᵢ](W'W)⁻¹vcov_rob_cem <- vcovHC(cem, type = "HC3", cluster = "group")
cem_rob <- coeftest(cem, vcov = vcov_rob_cem)
cetak_koef(as.matrix(cem_rob), "CEM + Panel-Robust VCV")##
## === CEM + Panel-Robust VCV ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 2.7462 7.6906 0.357 0.7216
## aset_log 0.4795 0.7907 0.606 0.5452
## der -1.4322 1.0781 -1.328 0.1863
## cr -0.3055 0.2122 -1.439 0.1524
## inflasi_z 0.7920 0.3135 2.527 0.0127 *
## birate_z 0.2875 0.3981 0.722 0.4715
## ihpb_z 1.0405 0.3595 2.894 0.0044 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## ─── Perbandingan SE: OLS vs Panel-Robust ───se_ols_cem <- coef(summary(cem))[,2]
se_rob_cem <- cem_rob[,2]
cat(sprintf("%-20s %12s %12s %10s\n","Variabel","SE OLS","SE Robust","Rasio"))## Variabel SE OLS SE Robust Rasio## --------------------------------------------------------for (nm in rownames(cem_rob)) {
cat(sprintf("%-20s %12.4f %12.4f %10.3f\n",
nm, se_ols_cem[nm], se_rob_cem[nm],
se_rob_cem[nm]/se_ols_cem[nm]))
}## (Intercept) 4.0237 7.6906 1.911
## aset_log 0.4232 0.7907 1.868
## der 0.7781 1.0781 1.386
## cr 0.1327 0.2122 1.599
## inflasi_z 0.3468 0.3135 0.904
## birate_z 0.3793 0.3981 1.050
## ihpb_z 0.3406 0.3595 1.056## Catatan: Rasio > 1 → SE robust lebih besar → OLS underestimate SE3.5 CEM-5: PENDUGAAN GLS
Jika \(\Omega\) dapat diestimasi, gunakan GLS:
\[\hat{\delta}_{GLS} = (W'\hat{\Omega}^{-1}W)^{-1}W'\hat{\Omega}^{-1}y\]
\[\widehat{Var}(\hat{\delta}_{GLS}) = (W'\hat{\Omega}^{-1}W)^{-1}\]
Estimasi \(\hat{\sigma}^2_i = \frac{1}{T}\sum_{t=1}^{T}\hat{e}^2_{it}\) (heteroskedastisitas antar unit), sehingga \(\hat{\Omega} = \text{diag}(\hat{\sigma}^2_1 I_T, \ldots, \hat{\sigma}^2_N I_T)\), bobot \(w_i = 1/\hat{\sigma}^2_i\).
## CEM-5: GLS — δ̂_GLS = (W'Ω̂⁻¹W)⁻¹W'Ω̂⁻¹y (Slide 43)# Step 1: Estimasi σ̂²_i dari residual CEM
# Slide 44: σ̂²_i = (1/T) Σ_t ê²_it
cat("Step 1: σ̂²_i = (1/T) Σ_t ê²_it (Slide 44)\n")## Step 1: σ̂²_i = (1/T) Σ_t ê²_it (Slide 44)kode_vec <- as.character(index(pdata)[[1]])
sigma2_i_cem <- tapply(resid_cem, kode_vec,
function(r) sum(r^2) / length(r))
cat(sprintf("%-12s %14s %10s\n","Perusahaan","σ̂²_i","σ̂_i"))## Perusahaan σ̂²_i σ̂_i## --------------------------------------for (nm in names(sigma2_i_cem))
cat(sprintf("%-12s %14.4f %10.4f\n", nm, sigma2_i_cem[nm], sqrt(sigma2_i_cem[nm])))## APLN 12.9337 3.5963
## ASRI 9.5779 3.0948
## BEST 5.0251 2.2417
## BSDE 4.5056 2.1226
## CTRA 6.2389 2.4978
## DILD 4.4336 2.1056
## DMAS 7.3948 2.7193
## JRPT 2.7907 1.6705
## LPCK 28.7812 5.3648
## LPKR 61.9859 7.8731
## MTLA 2.7834 1.6684
## PWON 4.8694 2.2067
## SMRA 5.0944 2.2571
## SSIA 12.3403 3.5129cat(sprintf("\nRasio σ̂²_max/σ̂²_min = %.1f → heteroskedastisitas: %s\n",
max(sigma2_i_cem)/min(sigma2_i_cem),
ifelse(max(sigma2_i_cem)/min(sigma2_i_cem) > 3,
"YA → GLS lebih efisien dari OLS", "ringan")))##
## Rasio σ̂²_max/σ̂²_min = 22.3 → heteroskedastisitas: YA → GLS lebih efisien dari OLS##
## Step 2: WLS, bobot w_i = 1/σ̂²_iw_gls_cem <- 1 / sigma2_i_cem[kode_vec]
gls_cem <- lm(Y ~ aset_log+der+cr+inflasi_z+birate_z+ihpb_z,
data=data_panel, weights=w_gls_cem)
resid_gls <- residuals(gls_cem)
sigma2_gls_cem <- sum(w_gls_cem * resid_gls^2) / (NT_ - K_ - 1)
cetak_koef(gls_cem, "CEM — GLS (Heteroskedastisitas Ω̂)")##
## === CEM — GLS (Heteroskedastisitas Ω̂) ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) -2.1532 2.5310 -0.851 0.3964
## aset_log 1.1114 0.2686 4.138 0.0001 ***
## der -1.9956 0.5172 -3.858 0.0002 ***
## cr -0.3260 0.0869 -3.752 0.0003 ***
## inflasi_z 0.4267 0.2194 1.945 0.0539 .
## birate_z 0.1995 0.2395 0.833 0.4062
## ihpb_z 0.6715 0.2156 3.115 0.0023 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## R² (GLS) = 0.3142# Slide 44-45 juga menyebut struktur Equal-Correlation dan AR(1)
# Implementasinya lebih lanjut dari cakupan materi dasar
cat("\nCatatan (Slide 44-45): Materi juga menyebut struktur Ω lain:\n")##
## Catatan (Slide 44-45): Materi juga menyebut struktur Ω lain:## - Equal-Correlation: ρ_ts = ρ untuk semua t≠s## - AR(1): ρ_ts = ρ^|t-s|## Estimasi ρ̂: Slide 44: ρ̂_ts = (1/T) Σ_t ê_it ê_js (rata-rata antar individu)## GLS dengan struktur ini menggunakan formula δ̂_GLS = (W'Ω̂⁻¹W)⁻¹W'Ω̂⁻¹y (Slide 43)3.6 CEM-6: PENDEKATAN ALTERNATIF CEM (Slide 40)
CEM memiliki tiga pendekatan alternatif estimasi.
## CEM-6: ESTIMATOR ALTERNATIF## Between, Within, First-Difference sebagai alternatif CEM## ─── Between Estimator: ȳ_i = x̄'_i β + u + ē_i ───cem_between <- plm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z,
data = pdata, model = "between", effect = "individual")
# Variabel makro yang time-varying tapi identik antar individu bisa di-drop
vars_btw <- names(coef(cem_between))
xvars_btw <- intersect(xvars, vars_btw)
cat(sprintf("Variabel tersedia: %s\n", paste(xvars_btw, collapse=", ")))## Variabel tersedia: aset_log, der, crdropped <- setdiff(xvars, xvars_btw)
if (length(dropped) > 0)
cat(sprintf("Di-drop (konstan antar grup): %s\n", paste(dropped, collapse=", ")))## Di-drop (konstan antar grup): inflasi_z, birate_z, ihpb_zresid_btw <- as.numeric(residuals(cem_between))
K_btw <- length(xvars_btw)
sigma2_btw <- sum(resid_btw^2) / (N_ - K_btw - 1)
cetak_koef(cem_between, "Between Estimator")##
## === Between Estimator ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 3.0756 7.6624 0.401 0.6966
## aset_log 0.5353 0.7953 0.673 0.5162
## der -2.0734 1.5761 -1.316 0.2177
## cr -0.4290 0.3223 -1.331 0.2128
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## σ²_e = ê'ê/(N-K-1) = 40.2528 / 10 = 4.025277## Catatan: Between hanya gunakan variasi antar individu (group means)## Within Estimator: ỹ_it = ỹ_it - ȳ_i ───cem_within <- plm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z,
data = pdata, model = "within", effect = "individual")
cetak_koef(cem_within, "Within Estimator (alternatif CEM)")##
## === Within Estimator (alternatif CEM) ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 0.5455 1.7769 0.307 0.7594
## der 0.8723 1.7545 0.497 0.6200
## cr -0.1665 0.1786 -0.933 0.3529
## inflasi_z 0.7855 0.3276 2.398 0.0180 *
## birate_z 0.3253 0.3550 0.916 0.3613
## ihpb_z 1.0781 0.3250 3.317 0.0012 **
## Signif: *** <.001 ** <.01 * <.05 . <.10## Catatan: Within hanya gunakan variasi dalam individu## First-Difference Estimator: Δy_it = Δx'_it β + Δe_itcem_fd <- plm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z,
data = pdata, model = "fd", effect = "individual")
cetak_koef(cem_fd, "First-Difference Estimator (alternatif CEM)")##
## === First-Difference Estimator (alternatif CEM) ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 0.6711 0.3798 1.767 0.0798 .
## aset_log -6.9510 3.0912 -2.249 0.0264 *
## der 7.1433 2.3310 3.064 0.0027 **
## cr -0.2333 0.1844 -1.265 0.2083
## inflasi_z 0.5382 0.2864 1.879 0.0627 .
## birate_z 0.1705 0.3938 0.433 0.6657
## ihpb_z 1.3805 0.3515 3.928 0.0001 ***
## Signif: *** <.001 ** <.01 * <.05 . <.10# Perbandingan ringkas ketiga alternatif
cat("\nPerbandingan Koefisien: OLS vs Between vs Within vs FD\n")##
## Perbandingan Koefisien: OLS vs Between vs Within vs FD## Variabel OLS Between Within FD## --------------------------------------------------------cf_btw <- coef(cem_between)
cf_wth <- coef(cem_within)
cf_fd <- coef(cem_fd)
for (xv in xvars) {
b_ols <- coef(cem)[xv]
b_btw <- ifelse(xv %in% names(cf_btw), cf_btw[xv], NA)
b_wth <- cf_wth[xv]
b_fd <- cf_fd[xv]
cat(sprintf("%-14s %10.4f %10s %10.4f %10.4f\n",
xv, b_ols,
ifelse(is.na(b_btw), " — ", sprintf("%10.4f", b_btw)),
b_wth, b_fd))
}## aset_log 0.4795 0.5353 0.5455 -6.9510
## der -1.4322 -2.0734 0.8723 7.1433
## cr -0.3055 -0.4290 -0.1665 -0.2333
## inflasi_z 0.7920 — 0.7855 0.5382
## birate_z 0.2875 — 0.3253 0.1705
## ihpb_z 1.0405 — 1.0781 1.3805## — : variabel di-drop oleh Between karena konstan antar grup4 BAGIAN II: FIXED EFFECT MODEL (FEM)
Model: \(y_{it} = x'_{it}\beta + u_i + e_{it}\)
\(u_i\) tetap (fixed), boleh berkorelasi dengan \(x_{it}\): \(\text{Cov}(u_i, x_{it}) \neq 0\)
4.1 FEM-1: DUMMY VARIABLE ESTIMATOR (LSDV)
\[W = [X\;\; D], \quad \hat{\delta}_{OLS} = (W'W)^{-1}W'y\]
\[\hat{\sigma}^2_e = \frac{\hat{e}'\hat{e}}{NT - N - K}\]
## FEM-1: DUMMY VARIABLE ESTIMATOR (LSDV)## W=[X|D], δ̂=(W'W)⁻¹W'y, σ²_e=ê'ê/(NT-N-K)# LSDV: tanpa intercept global, semua dummy perusahaan masuk
data_panel$kode <- as.character(index(pdata)[[1]])
fem_lsdv <- lm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z
+ factor(kode) - 1, data = data_panel)
resid_lsdv <- residuals(fem_lsdv)
sigma2_lsdv <- sum(resid_lsdv^2) / (NT_ - N_ - K_) # Slide 11: df=NT-N-K
sigma_lsdv <- sqrt(sigma2_lsdv)
# Koefisien slope β̂
cf_lsdv <- coef(summary(fem_lsdv))
cat("─── Koefisien Slope β̂ ───\n")## ─── Koefisien Slope β̂ ───##
## === LSDV — Slope β̂ ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 0.5455 1.7769 0.307 0.7594
## der 0.8723 1.7545 0.497 0.6200
## cr -0.1665 0.1786 -0.933 0.3529
## inflasi_z 0.7855 0.3276 2.398 0.0180 *
## birate_z 0.3253 0.3550 0.916 0.3613
## ihpb_z 1.0781 0.3250 3.317 0.0012 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## ─── Koefisien Dummy = û_i ───dummy_idx <- grep("factor", rownames(cf_lsdv))
cat(sprintf("%-28s %10s %10s %8s\n","Dummy (Perusahaan)","û_i","SE","t"))## Dummy (Perusahaan) û_i SE t## ------------------------------------------------------------for (i in dummy_idx)
cat(sprintf("%-28s %10.4f %10.4f %8.3f\n",
rownames(cf_lsdv)[i], cf_lsdv[i,1], cf_lsdv[i,2], cf_lsdv[i,3]))## factor(kode)APLN -0.4274 18.3040 -0.023
## factor(kode)ASRI -1.5193 18.2405 -0.083
## factor(kode)BEST -0.7470 15.5216 -0.048
## factor(kode)BSDE 1.6762 19.4959 0.086
## factor(kode)CTRA 1.0711 18.9651 0.056
## factor(kode)DILD -1.4560 17.4495 -0.083
## factor(kode)DMAS 3.0755 15.8427 0.194
## factor(kode)JRPT 1.5216 16.6430 0.091
## factor(kode)LPCK -1.3912 16.2025 -0.086
## factor(kode)LPKR -4.9801 19.7120 -0.253
## factor(kode)MTLA 1.3900 15.3556 0.091
## factor(kode)PWON 1.9083 18.2246 0.105
## factor(kode)SMRA -0.0450 18.5681 -0.002
## factor(kode)SSIA -1.8522 16.1775 -0.114cat(sprintf("\nσ²_e = ê'ê/(NT-N-K) = %.4f / %d = %.6f\n",
sum(resid_lsdv^2), NT_-N_-K_, sigma2_lsdv))##
## σ²_e = ê'ê/(NT-N-K) = 1254.8854 / 120 = 10.457378## σ_e = 3.2338 | R² = 0.7962##
## ─── Panel-Robust VCV (Slide 11-12) ───## V̂ar(δ̂) = (W'W)⁻¹[Σᵢ W'ᵢêᵢê'ᵢWᵢ](W'W)⁻¹vcov_lsdv_rob <- sandwich::vcovCL(fem_lsdv, cluster = ~factor(data_panel$kode),
data = data_panel)
lsdv_rob <- coeftest(fem_lsdv, vcov = vcov_lsdv_rob)
cetak_koef(lsdv_rob[xvars, , drop=FALSE], "LSDV + Panel-Robust SE")##
## === LSDV + Panel-Robust SE ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 0.5455 1.9846 0.275 0.7839
## der 0.8723 1.8249 0.478 0.6335
## cr -0.1665 0.1732 -0.962 0.3382
## inflasi_z 0.7855 0.2790 2.815 0.0057 **
## birate_z 0.3253 0.4691 0.693 0.4894
## ihpb_z 1.0781 0.3579 3.013 0.0032 **
## Signif: *** <.001 ** <.01 * <.05 . <.104.2 FEM-2: WITHIN ESTIMATOR (FE-OLS)
Transformasi within (\(Q_i = I_T - \frac{1}{T}ii'\)): \[\tilde{y}_{it} = y_{it} - \bar{y}_i, \quad \tilde{x}_{it} = x_{it} - \bar{x}_i\]
\[\hat{\beta}_{FE\text{-}OLS} = (\tilde{X}'\tilde{X})^{-1}\tilde{X}'\tilde{y}, \quad \hat{\sigma}^2_e = \frac{\hat{e}'\hat{e}}{NT - N - K}\]
Slide 17: \(\hat{\beta}_{GLS} = \hat{\beta}_{FE\text{-}OLS}\) karena \(Q_i^2 = Q_i\).
## FEM-2: WITHIN ESTIMATOR (FE-OLS)## β̂_FE = (X̃'X̃)⁻¹X̃'ỹ | Teorema: GLS = FE-OLSfem <- plm(Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z,
data = pdata, model = "within", effect = "individual")
resid_fem <- as.numeric(residuals(fem))
sigma2_fem <- sum(resid_fem^2) / (NT_ - N_ - K_) # Slide 16: df=NT-N-K
cetak_koef(fem, "FEM — Within Estimator (FE-OLS)")##
## === FEM — Within Estimator (FE-OLS) ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 0.5455 1.7769 0.307 0.7594
## der 0.8723 1.7545 0.497 0.6200
## cr -0.1665 0.1786 -0.933 0.3529
## inflasi_z 0.7855 0.3276 2.398 0.0180 *
## birate_z 0.3253 0.3550 0.916 0.3613
## ihpb_z 1.0781 0.3250 3.317 0.0012 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## σ²_e = ê'ê/(NT-N-K) = 1254.8854 / 120 = 10.457378## σ_e = 3.2338 | R² within = 0.2299##
## ─── Verifikasi: LSDV = Within ───diff_max <- max(abs(coef(fem)[xvars] - coef(fem_lsdv)[xvars]))
cat(sprintf("Max |β̂_Within - β̂_LSDV| = %.2e → %s\n",
diff_max, ifelse(diff_max < 1e-6, "IDENTIK ✓ ", "BERBEDA ✗")))## Max |β̂_Within - β̂_LSDV| = 7.63e-14 → IDENTIK ✓##
## ─── ML vs FE-OLS: bias koreksi σ²_e (Slide 19-20) ───sigma2_eML <- sum(resid_fem^2) / NT_
cat(sprintf("σ²_e(ML) = ê'ê/NT = %.6f ← BIASED (downward)\n", sigma2_eML))## σ²_e(ML) = ê'ê/NT = 8.963467 ← BIASED (downward)## σ²_e(FE) = ê'ê/(NT-N-K) = 10.457378 ← DIGUNAKAN## Panel seimbang : σ²_e(FE) = T/(T-1) × σ²_e(ML)cat(sprintf("Cek: %.0f/(%.0f-1) × %.6f = %.6f (vs σ²_FE=%.6f) → %s\n",
T_, T_, sigma2_eML, (T_/(T_-1))*sigma2_eML, sigma2_fem,
ifelse(abs((T_/(T_-1))*sigma2_eML - sigma2_fem) < 1e-6,
"COCOK ✓","cek ulang")))## Cek: 10/(10-1) × 8.963467 = 9.959408 (vs σ²_FE=10.457378) → cek ulang##
## ─── Panel-Robust VCV ───## V̂ar(β̂) = [ΣᵢX̃'ᵢX̃ᵢ]⁻¹[ΣᵢX̃'ᵢêᵢê'ᵢX̃ᵢ][ΣᵢX̃'ᵢX̃ᵢ]⁻¹vcov_rob_fem <- vcovHC(fem, type = "HC3", cluster = "group")
fem_rob <- coeftest(fem, vcov = vcov_rob_fem)
cetak_koef(as.matrix(fem_rob[,1:4]), "FEM-Within + Panel-Robust SE")##
## === FEM-Within + Panel-Robust SE ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 0.5455 1.8720 0.291 0.7713
## der 0.8723 1.7135 0.509 0.6116
## cr -0.1665 0.1756 -0.948 0.3448
## inflasi_z 0.7855 0.2604 3.017 0.0031 **
## birate_z 0.3253 0.4398 0.740 0.4610
## ihpb_z 1.0781 0.3391 3.180 0.0019 **
## Signif: *** <.001 ** <.01 * <.05 . <.104.3 FEM-3: ESTIMATED FIXED EFFECTS û_i
\[\hat{u}_i = \bar{y}_i - \bar{x}'_i\hat{\beta}\]
\[\widehat{Var}(\hat{u}_i) = \frac{\hat{\sigma}^2_e}{T_i} + \bar{x}'_i\widehat{Var}(\hat{\beta})\bar{x}_i\]
Untuk \(N \to \infty\): \(\hat{\beta}\) konsisten. Namun \(\hat{u}_i\) tidak konsisten kecuali \(T_i \to \infty\).
## FEM-3: ESTIMATED FIXED EFFECTS û_i## û_i = ȳ_i - x̄'_i β̂## Var(û_i) = σ²_e/T_i + x̄'_i Var(β̂) x̄_ife_vals <- fixef(fem)
vcov_b <- vcov(fem)
data_panel$kode <- as.character(index(pdata)[[1]])
# Rata-rata per individu
x_bar <- data_panel %>%
group_by(kode) %>%
summarise(across(all_of(xvars), mean), .groups="drop") %>%
arrange(kode)
cat(sprintf("%-12s %10s %10s %12s %12s\n",
"Perusahaan","û_i","SE(û_i)","95% CI lo","95% CI hi"))## Perusahaan û_i SE(û_i) 95% CI lo 95% CI hi## ------------------------------------------------------------fe_result <- lapply(seq_along(fe_vals), function(i) {
kd <- names(fe_vals)[i]
ui <- as.numeric(fe_vals[i])
xi_b <- as.numeric(x_bar[x_bar$kode==kd, xvars])
var_u <- sigma2_fem/T_ + as.numeric(t(xi_b) %*% vcov_b %*% xi_b)
se_u <- sqrt(var_u)
ci_lo <- ui - 1.96*se_u; ci_hi <- ui + 1.96*se_u
cat(sprintf("%-12s %10.4f %10.4f %12.4f %12.4f\n",
kd, ui, se_u, ci_lo, ci_hi))
data.frame(kode=kd, ui=ui, se_u=se_u, ci_lo=ci_lo, ci_hi=ci_hi)
}) %>% bind_rows()## APLN -0.4274 18.3040 -36.3032 35.4484
## ASRI -1.5193 18.2405 -37.2708 34.2321
## BEST -0.7470 15.5216 -31.1693 29.6753
## BSDE 1.6762 19.4959 -36.5357 39.8881
## CTRA 1.0711 18.9651 -36.1005 38.2428
## DILD -1.4560 17.4495 -35.6570 32.7449
## DMAS 3.0755 15.8427 -27.9761 34.1272
## JRPT 1.5216 16.6430 -31.0987 34.1419
## LPCK -1.3912 16.2025 -33.1482 30.3658
## LPKR -4.9801 19.7120 -43.6156 33.6554
## MTLA 1.3900 15.3556 -28.7069 31.4869
## PWON 1.9083 18.2246 -33.8119 37.6286
## SMRA -0.0450 18.5681 -36.4385 36.3486
## SSIA -1.8522 16.1775 -33.5601 29.8557# Plot û_i
p_fe <- ggplot(fe_result, aes(x=reorder(kode,ui), y=ui, fill=ui>0)) +
geom_col(width=0.7, alpha=0.85) +
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=0.3, linewidth=0.6) +
geom_hline(yintercept=0, linewidth=0.9) +
scale_fill_manual(values=c("TRUE"="#2c7bb6","FALSE"="#e63946")) +
coord_flip() + theme_minimal(base_size=12) +
theme(legend.position="none", panel.grid.minor=element_blank(),
plot.title=element_text(face="bold")) +
labs(title="Fixed Effects û_i per Perusahaan",
subtitle="Biru = di atas rata-rata | Merah = di bawah | Bar = 95% CI",
x=NULL, y="û_i")
print(p_fe)
## Catatan: N→∞ membuat β̂ konsisten, tetapi û_i inkonsisten kecuali T_i→∞4.4 FEM-4: ASUMSI KLASIK FEM
## FEM-4: UJI ASUMSI KLASIK FEM (Slide 5)# [A1] Strict Exogeneity → divalidasi via F-test (Bagian III)
cat("[A1] STRICT EXOGENEITY: E(e_it|X) = 0\n")## [A1] STRICT EXOGENEITY: E(e_it|X) = 0## Divalidasi formal via F-test poolability (Bagian III)## [A2] HOMOSKEDASTISITAS: Var(e_it|X) = σ²_ebp_fem <- bptest(fem)
cat(sprintf(" nR² = %.4f, df=%d, p = %.6f → %s\n\n",
bp_fem$statistic, bp_fem$parameter, bp_fem$p.value,
ifelse(bp_fem$p.value < 0.05, "DILANGGAR ⚠", "TERPENUHI ✓")))## nR² = 24.0467, df=6, p = 0.000512 → DILANGGAR ⚠## [A3] NO SERIAL CORRELATION: Cov(e_it, e_is|X) = 0## d_L = 1.662, d_U = 1.802 (tabel DW, n=140, K=6, α=5%)dw_fem <- pdwtest(fem)
dw_ket <- ifelse(dw_fem$statistic < dL, "AUTOKORELASI POSITIF ⚠",
ifelse(dw_fem$statistic <= dU, "INCONCLUSIVE",
ifelse(dw_fem$statistic < 4-dU, "TIDAK ADA autokorelasi ✓",
ifelse(dw_fem$statistic < 4-dL, "INCONCLUSIVE",
"AUTOKORELASI NEGATIF ⚠"))))
cat(sprintf(" DW = %.4f → %s\n\n", dw_fem$statistic, dw_ket))## DW = 1.8072 → TIDAK ADA autokorelasi ✓## [A4] NORMALITAS: e_it ~ N(0,σ²_e)sw_fem <- shapiro.test(resid_fem)
cat(sprintf(" SW W = %.4f, p = %.2e → %s\n\n",
sw_fem$statistic, sw_fem$p.value,
ifelse(sw_fem$p.value > 0.05, "TERPENUHI ✓", "DILANGGAR ⚠")))## SW W = 0.8933, p = 1.37e-08 → DILANGGAR ⚠## [A5] NO MULTIKOLINEARITAS: VIF## aset_log der cr inflasi_z birate_z ihpb_z
## 1.227 1.385 1.251 1.318 1.576 1.271## → Semua VIF < 10 ✓# Ringkasan
ring_fem <- data.frame(
Asumsi = c("A1 Exogeneity (F-test)","A2 Homoskedastisitas",
"A3 No Serial Corr (DW)","A4 Normalitas","A5 No Multikolinearitas"),
Status = c(
"→ lihat F-test (Bag. III)",
ifelse(bp_fem$p.value < 0.05, "DILANGGAR ⚠", "TERPENUHI ✓"),
ifelse(grepl("TIDAK", dw_ket), "TERPENUHI ✓",
ifelse(grepl("INCON", dw_ket), "INCONCLUSIVE", "DILANGGAR ⚠")),
ifelse(sw_fem$p.value > 0.05, "TERPENUHI ✓", "DILANGGAR ⚠"),
ifelse(all(vif_vals < 10), "TERPENUHI ✓", "MASALAH ⚠")),
Penanganan = c("F-test → Bagian III",
"Panel-Robust VCV (Slide 11)",
"Panel-Robust VCV (Slide 22)",
"Transformasi signed-log (sudah)",
"Tidak perlu")
)
knitr::kable(ring_fem, caption="Ringkasan Asumsi FEM") %>%
kableExtra::kable_styling(bootstrap_options=c("striped","hover","condensed"))| Asumsi | Status | Penanganan |
|---|---|---|
| A1 Exogeneity (F-test) | → lihat F-test (Bag. III) | F-test → Bagian III |
| A2 Homoskedastisitas | DILANGGAR ⚠ | Panel-Robust VCV (Slide 11) |
| A3 No Serial Corr (DW) | TERPENUHI ✓ | Panel-Robust VCV (Slide 22) |
| A4 Normalitas | DILANGGAR ⚠ | Transformasi signed-log (sudah) |
| A5 No Multikolinearitas | TERPENUHI ✓ | Tidak perlu |
# Plot diagnostik residual
par(mfrow=c(2,4), mar=c(4,4,3,1))
yhat_fem <- as.numeric(fitted(fem))
plot(yhat_fem, resid_fem, pch=16,cex=0.7,col="#2c7bb6",
main="Resid vs Fitted",xlab="ŷ",ylab="Residual")
abline(h=0,col="red",lty=2)
for (k in 1:6)
plot(data_panel[[xvars[k]]], resid_fem, pch=16,cex=0.7,col="#2c7bb6",
main=paste("Resid vs",nms_x[k]),xlab=nms_x[k],ylab="Residual")
plot(data_panel$tahun, resid_fem, pch=16,cex=0.7,col="#2c7bb6",
main="Resid vs Waktu",xlab="Tahun",ylab="Residual")
abline(h=0,col="red",lty=2)
par(mfrow=c(1,2), mar=c(4,4,3,1))
qqnorm(resid_fem, pch=16,cex=0.7,col="#2c7bb6",main="Normal Q-Q Plot (FEM)")
qqline(resid_fem, col="red",lty=2)
hist(resid_fem, breaks=25,col="#2c7bb6",border="white",
main="Distribusi Residual FEM",xlab="Residual",freq=FALSE)
curve(dnorm(x,mean(resid_fem),sd(resid_fem)),add=TRUE,col="black",lwd=2)
4.5 FEM-5: PANEL-ROBUST VCV
## FEM-5: PANEL-ROBUST VCV (Slide 22)## V̂ar(β̂) = [ΣᵢX̃'ᵢX̃ᵢ]⁻¹[ΣᵢX̃'ᵢêᵢê'ᵢX̃ᵢ][ΣᵢX̃'ᵢX̃ᵢ]⁻¹# GLS-WLS untuk FEM (heteroskedastisitas)
cat("─── GLS-WLS FEM: mengatasi heteroskedastisitas ───\n")## ─── GLS-WLS FEM: mengatasi heteroskedastisitas ───## Slide 44: σ̂²_i = (1/T) Σ_t ê²_itkode_vec <- as.character(index(pdata)[[1]])
sigma2_i_fem <- tapply(resid_fem, kode_vec,
function(r) sum(r^2)/length(r))
cat(sprintf("Rasio σ̂²_max/σ̂²_min = %.1f\n",
max(sigma2_i_fem)/min(sigma2_i_fem)))## Rasio σ̂²_max/σ̂²_min = 29.5w_gls_fem <- 1 / sigma2_i_fem[kode_vec]
gls_fem <- lm(Y ~ aset_log+der+cr+inflasi_z+birate_z+ihpb_z
+ factor(kode) - 1, data=data_panel, weights=w_gls_fem)
cetak_koef(coef(summary(gls_fem))[xvars,,drop=FALSE], "FEM — GLS-WLS")##
## === FEM — GLS-WLS ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## aset_log 1.5629 0.9525 1.641 0.1035
## der -0.7692 1.1447 -0.672 0.5029
## cr -0.2453 0.0963 -2.548 0.0121 *
## inflasi_z 0.3532 0.1814 1.947 0.0539 .
## birate_z 0.2009 0.1959 1.025 0.3073
## ihpb_z 0.5487 0.1840 2.982 0.0035 **
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## ─── Perbandingan SE: OLS vs Panel-Robust ───se_ols_fem <- coef(summary(fem))[,2]
se_rob_fem <- fem_rob[,2]
cat(sprintf("%-14s %12s %12s %10s\n","Variabel","SE OLS","SE Robust","Rasio"))## Variabel SE OLS SE Robust Rasio## --------------------------------------------------for (nm in names(se_ols_fem))
cat(sprintf("%-14s %12.4f %12.4f %10.3f\n",
nm, se_ols_fem[nm], se_rob_fem[nm],
se_rob_fem[nm]/se_ols_fem[nm]))## aset_log 1.7769 1.8720 1.053
## der 1.7545 1.7135 0.977
## cr 0.1786 0.1756 0.983
## inflasi_z 0.3276 0.2604 0.795
## birate_z 0.3550 0.4398 1.239
## ihpb_z 0.3250 0.3391 1.0434.6 FEM-6: FIRST-DIFFERENCE ESTIMATOR
\(\Delta y_{it} = \Delta x'_{it}\beta + \Delta e_{it}\)
\(\text{Var}(\Delta e_{it}|X) = 2\sigma^2_e\), struktur Toeplitz.
\(\hat{\beta}_{FD} = (\Delta X'\Delta X)^{-1}\Delta X'\Delta y\)
\(\hat{\beta}_{GLS} = (\Delta X'\Omega^{-1}\Delta X)^{-1}\Delta X'\Omega^{-1}\Delta y\)
## FEM-6: FIRST-DIFFERENCE ESTIMATOR## Δy_it = Δx'_it β + Δe_it | Var(Δe)=2σ²_e# FD-OLS
fem_fd <- plm(Y ~ aset_log+der+cr+inflasi_z+birate_z+ihpb_z,
data=pdata, model="fd", effect="individual")
resid_fd <- as.numeric(residuals(fem_fd))
NT_fd <- N_ * (T_ - 1)
df_fd <- NT_fd - K_ # df = N(T-1)-K (tanpa intercept)
sigma2_de <- sum(resid_fd^2) / df_fd
sigma2_e_fd<- sigma2_de / 2 # Var(Δe)=2σ²_e → σ²_e = σ²_Δe/2
cetak_koef(fem_fd, "FD-OLS")##
## === FD-OLS ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 0.6711 0.3798 1.767 0.0798 .
## aset_log -6.9510 3.0912 -2.249 0.0264 *
## der 7.1433 2.3310 3.064 0.0027 **
## cr -0.2333 0.1844 -1.265 0.2083
## inflasi_z 0.5382 0.2864 1.879 0.0627 .
## birate_z 0.1705 0.3938 0.433 0.6657
## ihpb_z 1.3805 0.3515 3.928 0.0001 ***
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## σ²_Δe = Δê'Δê/(N(T-1)-K) = 1665.1498 / 120 = 13.876248## σ²_e = σ²_Δe / 2 = 6.938124 (Var(Δe)=2σ²_e)## σ_e = 2.6340##
## ─── FD-GLS: Panel-Robust VCV───## Koreksi struktur Toeplitz: Cov(Δe_it,Δe_is) = {2σ²_e jika t=s; -σ²_e jika |t-s|=1}vcov_rob_fd <- vcovHC(fem_fd, type="HC3", cluster="group")
fd_gls <- coeftest(fem_fd, vcov=vcov_rob_fd)
cetak_koef(as.matrix(fd_gls[,1:4]), "FD-GLS (Slide 26)")##
## === FD-GLS (Slide 26) ===
## Variabel Koef SE t p-value
## --------------------------------------------------------------
## (Intercept) 0.6711 0.1533 4.376 0.0000 ***
## aset_log -6.9510 1.9176 -3.625 0.0004 ***
## der 7.1433 3.2610 2.191 0.0304 *
## cr -0.2333 0.2721 -0.857 0.3930
## inflasi_z 0.5382 0.2188 2.459 0.0154 *
## birate_z 0.1705 0.3311 0.515 0.6074
## ihpb_z 1.3805 0.5721 2.413 0.0173 *
## Signif: *** <.001 ** <.01 * <.05 . <.10##
## ─── Perbandingan β̂: Within vs FD ───## (Harus mendekati sama jika e_it adalah white noise)## Variabel β̂_Within β̂_FD Selisih## ----------------------------------------------------for (xv in xvars)
cat(sprintf("%-14s %12.4f %12.4f %12.4f\n",
xv, coef(fem)[xv], coef(fem_fd)[xv],
coef(fem)[xv]-coef(fem_fd)[xv]))## aset_log 0.5455 -6.9510 7.4965
## der 0.8723 7.1433 -6.2710
## cr -0.1665 -0.2333 0.0668
## inflasi_z 0.7855 0.5382 0.2473
## birate_z 0.3253 0.1705 0.1547
## ihpb_z 1.0781 1.3805 -0.30245 BAGIAN III: UJI F POOLABILITY (TO POOL OR NOT TO POOL?)
\[H_0: u_i = u\;\forall i \quad \text{(CEM cukup)}\] \[H_1: \exists\; u_i \neq u_j \quad \text{(FEM diperlukan)}\]
\[F = \frac{(RSS_R - RSS_{UR})/(N-1)}{RSS_{UR}/(NT-N-K)} \sim F(N-1,\; NT-N-K)\]
\[RSS_{UR} = \hat{e}'_{FE}\hat{e}_{FE}, \quad RSS_R = \hat{e}'_{PO}\hat{e}_{PO}\]
## III: UJI F POOLABILITY## H0: u_i = u untuk semua i (CEM cukup)## H1: ada u_i ≠ u_j (FEM diperlukan)## F = [(RSS_R-RSS_UR)/(N-1)] / [RSS_UR/(NT-N-K)]rss_r <- sum(resid_cem^2) # RSS CEM (restricted)
rss_ur <- sum(resid_fem^2) # RSS FEM (unrestricted)
F_hit <- ((rss_r - rss_ur)/(N_-1)) / (rss_ur/(NT_-N_-K_))
F_kri <- qf(0.95, N_-1, NT_-N_-K_)
p_F <- pf(F_hit, N_-1, NT_-N_-K_, lower.tail=FALSE)
cat(sprintf("RSS_R (CEM) = ê'_PO ê_PO = %.4f\n", rss_r))## RSS_R (CEM) = ê'_PO ê_PO = 1687.5502## RSS_UR (FEM) = ê'_FE ê_FE = 1254.8854cat(sprintf("Pembilang = (RSS_R - RSS_UR)/(N-1) = %.4f / %d = %.4f\n",
rss_r-rss_ur, N_-1, (rss_r-rss_ur)/(N_-1)))## Pembilang = (RSS_R - RSS_UR)/(N-1) = 432.6648 / 13 = 33.2819cat(sprintf("Penyebut = RSS_UR/(NT-N-K) = %.4f / %d = %.4f\n",
rss_ur, NT_-N_-K_, rss_ur/(NT_-N_-K_)))## Penyebut = RSS_UR/(NT-N-K) = 1254.8854 / 120 = 10.4574##
## F_hitung = 3.1826## F_kritis (α=5%) = F(13,120) = 1.8026## p-value = 0.000397## Konfirmasi dengan pFtest():##
## F test for individual effects
##
## data: Y ~ aset_log + der + cr + inflasi_z + birate_z + ihpb_z
## F = 3.1826, df1 = 13, df2 = 120, p-value = 0.0003966
## alternative hypothesis: significant effectscat(sprintf("\n─── KESIMPULAN ───\n%s\n",
ifelse(p_F < 0.05,
"TOLAK H0 pada α=5%: efek individu u_i berbeda antar perusahaan.\nFEM lebih tepat dari CEM.",
"GAGAL TOLAK H0 pada α=5%: tidak cukup bukti perbedaan efek individu.\nCEM memadai secara statistik.")))##
## ─── KESIMPULAN ───
## TOLAK H0 pada α=5%: efek individu u_i berbeda antar perusahaan.
## FEM lebih tepat dari CEM.```