# Menentukan CRAN mirror
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Paket yang diperlukan
packages <- c("plm", "lmtest", "sandwich",
"car", "stargazer")
# Instalasi otomatis jika paket belum tersedia
for (p in packages) {
if (!requireNamespace(p, quietly = TRUE)) {
install.packages(p, repos = "https://cloud.r-project.org")
}
}
# Memanggil paket
library(plm)
## Warning: package 'plm' was built under R version 4.4.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.4.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sandwich)
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
library(stargazer)
##
## Please cite as:
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
# 2. MEMANGGIL DAN MEMAHAMI DATA
# Memanggil dataset Grunfeld
data("Grunfeld", package = "plm")
# Melihat data
head(Grunfeld)
str(Grunfeld)
## 'data.frame': 200 obs. of 5 variables:
## $ firm : int 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 ...
## $ inv : num 318 392 411 258 331 ...
## $ value : num 3078 4662 5387 2792 4313 ...
## $ capital: num 2.8 52.6 156.9 209.2 203.4 ...
summary(Grunfeld)
## firm year inv value
## Min. : 1.0 Min. :1935 Min. : 0.93 Min. : 58.12
## 1st Qu.: 3.0 1st Qu.:1940 1st Qu.: 33.56 1st Qu.: 199.97
## Median : 5.5 Median :1944 Median : 57.48 Median : 517.95
## Mean : 5.5 Mean :1944 Mean : 145.96 Mean :1081.68
## 3rd Qu.: 8.0 3rd Qu.:1949 3rd Qu.: 138.04 3rd Qu.:1679.85
## Max. :10.0 Max. :1954 Max. :1486.70 Max. :6241.70
## capital
## Min. : 0.80
## 1st Qu.: 79.17
## Median : 205.60
## Mean : 276.02
## 3rd Qu.: 358.10
## Max. :2226.30
# Jumlah observasi dan variabel
dim(Grunfeld)
## [1] 200 5
# Melihat nama variabel
names(Grunfeld)
## [1] "firm" "year" "inv" "value" "capital"
# 3. DEKLARASI STRUKTUR DATA PANEL
# Menentukan unit cross-section dan waktu
pdata <- pdata.frame(
Grunfeld,
index = c("firm", "year")
)
# Memeriksa struktur panel
pdim(pdata)
## Balanced Panel: n = 10, T = 20, N = 200
summary(pdata)
## firm year inv value
## 1 :20 1935 : 10 Min. : 0.93 Min. : 58.12
## 2 :20 1936 : 10 1st Qu.: 33.56 1st Qu.: 199.97
## 3 :20 1937 : 10 Median : 57.48 Median : 517.95
## 4 :20 1938 : 10 Mean : 145.96 Mean :1081.68
## 5 :20 1939 : 10 3rd Qu.: 138.04 3rd Qu.:1679.85
## 6 :20 1940 : 10 Max. :1486.70 Max. :6241.70
## (Other):80 (Other):140
## capital
## Min. : 0.80
## 1st Qu.: 79.17
## Median : 205.60
## Mean : 276.02
## 3rd Qu.: 358.10
## Max. :2226.30
##
# 4. ESTIMASI TIGA MODEL DATA PANEL
# Model regresi
form <- inv ~ value + capital
# Common Effect Model (CEM)
cem <- plm(
form,
data = pdata,
model = "pooling"
)
# Fixed Effect Model (FEM)
fem <- plm(
form,
data = pdata,
model = "within"
)
# Random Effect Model (REM)
rem <- plm(
form,
data = pdata,
model = "random"
)
# Melihat hasil estimasi
summary(cem)
## Pooling Model
##
## Call:
## plm(formula = form, data = pdata, model = "pooling")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -291.6757 -30.0137 5.3033 34.8293 369.4464
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -42.7143694 9.5116760 -4.4907 1.207e-05 ***
## value 0.1155622 0.0058357 19.8026 < 2.2e-16 ***
## capital 0.2306785 0.0254758 9.0548 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 9359900
## Residual Sum of Squares: 1755900
## R-Squared: 0.81241
## Adj. R-Squared: 0.8105
## F-statistic: 426.576 on 2 and 197 DF, p-value: < 2.22e-16
summary(fem)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = form, data = pdata, model = "within")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -184.00857 -17.64316 0.56337 19.19222 250.70974
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.110124 0.011857 9.2879 < 2.2e-16 ***
## capital 0.310065 0.017355 17.8666 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2244400
## Residual Sum of Squares: 523480
## R-Squared: 0.76676
## Adj. R-Squared: 0.75311
## F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16
summary(rem)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = form, data = pdata, model = "random")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Effects:
## var std.dev share
## idiosyncratic 2784.46 52.77 0.282
## individual 7089.80 84.20 0.718
## theta: 0.8612
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -177.6063 -19.7350 4.6851 19.5105 252.8743
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -57.834415 28.898935 -2.0013 0.04536 *
## value 0.109781 0.010493 10.4627 < 2e-16 ***
## capital 0.308113 0.017180 17.9339 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2381400
## Residual Sum of Squares: 548900
## R-Squared: 0.7695
## Adj. R-Squared: 0.76716
## Chisq: 657.674 on 2 DF, p-value: < 2.22e-16
# 5. EFEK INDIVIDU PADA FEM
fixef(fem)
## 1 2 3 4 5 6 7 8
## -70.2967 101.9058 -235.5718 -27.8093 -114.6168 -23.1613 -66.5535 -57.5457
## 9 10
## -87.2223 -6.5678
# 6. PERBANDINGAN HASIL MODEL
stargazer(
cem, fem, rem,
type = "text",
column.labels = c("CEM", "FEM", "REM"),
title = "Hasil Estimasi Model Data Panel",
dep.var.labels = "Investasi"
)
##
## Hasil Estimasi Model Data Panel
## =========================================================================
## Dependent variable:
## ------------------------------------------------------------
## Investasi
## CEM FEM REM
## (1) (2) (3)
## -------------------------------------------------------------------------
## value 0.116*** 0.110*** 0.110***
## (0.006) (0.012) (0.010)
##
## capital 0.231*** 0.310*** 0.308***
## (0.025) (0.017) (0.017)
##
## Constant -42.714*** -57.834**
## (9.512) (28.899)
##
## -------------------------------------------------------------------------
## Observations 200 200 200
## R2 0.812 0.767 0.770
## Adjusted R2 0.811 0.753 0.767
## F Statistic 426.576*** (df = 2; 197) 309.014*** (df = 2; 188) 657.674***
## =========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
# 7. UJI PEMILIHAN MODEL
# Uji Chow (CEM vs FEM)
uji_chow <- pFtest(fem, cem)
uji_chow
##
## F test for individual effects
##
## data: form
## F = 49.177, df1 = 9, df2 = 188, p-value < 2.2e-16
## alternative hypothesis: significant effects
# Uji Breusch-Pagan LM (CEM vs REM)
uji_lm <- plmtest(cem, type = "bp")
uji_lm
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: form
## chisq = 798.16, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
# Uji Hausman (FEM vs REM)
uji_hausman <- phtest(fem, rem)
uji_hausman
##
## Hausman Test
##
## data: form
## chisq = 2.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent
# 8. TWO-WAY FIXED EFFECT
# Model efek individu dan waktu
fem_tw <- plm(
form,
data = pdata,
model = "within",
effect = "twoways"
)
summary(fem_tw)
## Twoways effects Within Model
##
## Call:
## plm(formula = form, data = pdata, effect = "twoways", model = "within")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -162.6094 -19.4710 -1.2669 19.1277 211.8420
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.117716 0.013751 8.5604 6.653e-15 ***
## capital 0.357916 0.022719 15.7540 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1615600
## Residual Sum of Squares: 452150
## R-Squared: 0.72015
## Adj. R-Squared: 0.67047
## F-statistic: 217.442 on 2 and 169 DF, p-value: < 2.22e-16
# Menguji apakah efek waktu diperlukan
pFtest(fem_tw, fem)
##
## F test for twoways effects
##
## data: form
## F = 1.4032, df1 = 19, df2 = 169, p-value = 0.1309
## alternative hypothesis: significant effects
# 9. UJI ASUMSI
# A. Uji heteroskedastisitas
uji_hetero <- bptest(
inv ~ value + capital + factor(firm),
data = Grunfeld
)
uji_hetero
##
## studentized Breusch-Pagan test
##
## data: inv ~ value + capital + factor(firm)
## BP = 85.836, df = 11, p-value = 1.086e-13
# B. Uji autokorelasi
uji_auto <- pbgtest(fem)
uji_auto
##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: form
## chisq = 65.063, df = 20, p-value = 1.14e-06
## alternative hypothesis: serial correlation in idiosyncratic errors
# Uji Wooldridge
uji_wooldridge <- pwartest(fem)
uji_wooldridge
##
## Wooldridge's test for serial correlation in FE panels
##
## data: fem
## F = 76.929, df1 = 1, df2 = 188, p-value = 1.057e-15
## alternative hypothesis: serial correlation
# C. Uji dependensi antarunit
uji_dependensi <- pcdtest(fem, test = "cd")
uji_dependensi
##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: inv ~ value + capital
## z = 4.6612, p-value = 3.144e-06
## alternative hypothesis: cross-sectional dependence
# D. Uji multikolinearitas
model_lm <- lm(
inv ~ value + capital,
data = Grunfeld
)
vif(model_lm)
## value capital
## 1.313773 1.313773
# Korelasi antarvariabel independen
cor(Grunfeld[, c("value", "capital")])
## value capital
## value 1.0000000 0.4887058
## capital 0.4887058 1.0000000
# 10. UJI STASIONERITAS
uji_stasioner <- purtest(
pdata$inv,
test = "levinlin",
exo = "intercept",
lags = "AIC",
pmax = 4
)
## Warning in selectT(l, theTs): the time series is short
uji_stasioner
##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$inv
## z = 3.1416, p-value = 0.9992
## alternative hypothesis: stationarity
# 11. STANDARD ERROR ROBUST
# A. Cluster per unit
robust_cluster <- coeftest(
fem,
vcov = vcovHC(
fem,
type = "HC1",
cluster = "group"
)
)
robust_cluster
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.014414 7.6398 1.076e-12 ***
## capital 0.310065 0.050043 6.1959 3.576e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# B. Driscoll-Kraay
robust_dk <- coeftest(
fem,
vcov = vcovSCC(
fem,
type = "HC1"
)
)
robust_dk
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.017775 6.1954 3.586e-09 ***
## capital 0.310065 0.034996 8.8601 6.017e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 12. KESIMPULAN HASIL
cat("ANALISIS DATA PANEL SELESAI\n")
## ANALISIS DATA PANEL SELESAI
cat("Jumlah observasi:", nrow(Grunfeld), "\n")
## Jumlah observasi: 200
cat("Jumlah perusahaan:",
length(unique(Grunfeld$firm)), "\n")
## Jumlah perusahaan: 10
cat("Periode:",
min(Grunfeld$year), "-",
max(Grunfeld$year), "\n")
## Periode: 1935 - 1954