# 1. INSTALL PACKAGE
packages <- c(
"plm",
"lmtest",
"car",
"sandwich",
"stargazer"
)
# Mengecek dan menginstal package yang belum tersedia
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
# 2. MEMANGGIL PACKAGE
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(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(sandwich)
## Warning: package 'sandwich' 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
# Mengatur tampilan angka agar tidak menggunakan notasi ilmiah
options(scipen = 999)
# 2. MEMANGGIL DAN MENGENALI DATASET
# Memanggil dataset bawaan package plm
data("Grunfeld", package = "plm")
# Menampilkan 6 observasi pertama
head(Grunfeld)
## firm year inv value capital
## 1 1 1935 317.6 3078.5 2.8
## 2 1 1936 391.8 4661.7 52.6
## 3 1 1937 410.6 5387.1 156.9
## 4 1 1938 257.7 2792.2 209.2
## 5 1 1939 330.8 4313.2 203.4
## 6 1 1940 461.2 4643.9 207.2
# Menampilkan 6 observasi terakhir
tail(Grunfeld)
## firm year inv value capital
## 195 10 1949 4.21 71.24 7.39
## 196 10 1950 3.42 69.05 8.74
## 197 10 1951 4.67 83.04 9.07
## 198 10 1952 6.00 74.42 9.93
## 199 10 1953 6.53 63.51 11.68
## 200 10 1954 5.12 58.12 14.33
# Melihat struktur dataset
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 ...
# Melihat nama variabel
names(Grunfeld)
## [1] "firm" "year" "inv" "value" "capital"
# Melihat dimensi dataset
dim(Grunfeld)
## [1] 200 5
# Melihat statistik deskriptif
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
# 3. DEKLARASI STRUKTUR DATA PANEL
# Menetapkan firm sebagai unit cross-section
# Menetapkan year sebagai unit time series
pdata <- pdata.frame(
Grunfeld,
index = c("firm", "year")
)
# Memeriksa struktur panel
pdim(pdata)
## Balanced Panel: n = 10, T = 20, N = 200
# Melihat indeks panel
index(pdata)
## firm year
## 1 1 1935
## 2 1 1936
## 3 1 1937
## 4 1 1938
## 5 1 1939
## 6 1 1940
## 7 1 1941
## 8 1 1942
## 9 1 1943
## 10 1 1944
## 11 1 1945
## 12 1 1946
## 13 1 1947
## 14 1 1948
## 15 1 1949
## 16 1 1950
## 17 1 1951
## 18 1 1952
## 19 1 1953
## 20 1 1954
## 21 2 1935
## 22 2 1936
## 23 2 1937
## 24 2 1938
## 25 2 1939
## 26 2 1940
## 27 2 1941
## 28 2 1942
## 29 2 1943
## 30 2 1944
## 31 2 1945
## 32 2 1946
## 33 2 1947
## 34 2 1948
## 35 2 1949
## 36 2 1950
## 37 2 1951
## 38 2 1952
## 39 2 1953
## 40 2 1954
## 41 3 1935
## 42 3 1936
## 43 3 1937
## 44 3 1938
## 45 3 1939
## 46 3 1940
## 47 3 1941
## 48 3 1942
## 49 3 1943
## 50 3 1944
## 51 3 1945
## 52 3 1946
## 53 3 1947
## 54 3 1948
## 55 3 1949
## 56 3 1950
## 57 3 1951
## 58 3 1952
## 59 3 1953
## 60 3 1954
## 61 4 1935
## 62 4 1936
## 63 4 1937
## 64 4 1938
## 65 4 1939
## 66 4 1940
## 67 4 1941
## 68 4 1942
## 69 4 1943
## 70 4 1944
## 71 4 1945
## 72 4 1946
## 73 4 1947
## 74 4 1948
## 75 4 1949
## 76 4 1950
## 77 4 1951
## 78 4 1952
## 79 4 1953
## 80 4 1954
## 81 5 1935
## 82 5 1936
## 83 5 1937
## 84 5 1938
## 85 5 1939
## 86 5 1940
## 87 5 1941
## 88 5 1942
## 89 5 1943
## 90 5 1944
## 91 5 1945
## 92 5 1946
## 93 5 1947
## 94 5 1948
## 95 5 1949
## 96 5 1950
## 97 5 1951
## 98 5 1952
## 99 5 1953
## 100 5 1954
## 101 6 1935
## 102 6 1936
## 103 6 1937
## 104 6 1938
## 105 6 1939
## 106 6 1940
## 107 6 1941
## 108 6 1942
## 109 6 1943
## 110 6 1944
## 111 6 1945
## 112 6 1946
## 113 6 1947
## 114 6 1948
## 115 6 1949
## 116 6 1950
## 117 6 1951
## 118 6 1952
## 119 6 1953
## 120 6 1954
## 121 7 1935
## 122 7 1936
## 123 7 1937
## 124 7 1938
## 125 7 1939
## 126 7 1940
## 127 7 1941
## 128 7 1942
## 129 7 1943
## 130 7 1944
## 131 7 1945
## 132 7 1946
## 133 7 1947
## 134 7 1948
## 135 7 1949
## 136 7 1950
## 137 7 1951
## 138 7 1952
## 139 7 1953
## 140 7 1954
## 141 8 1935
## 142 8 1936
## 143 8 1937
## 144 8 1938
## 145 8 1939
## 146 8 1940
## 147 8 1941
## 148 8 1942
## 149 8 1943
## 150 8 1944
## 151 8 1945
## 152 8 1946
## 153 8 1947
## 154 8 1948
## 155 8 1949
## 156 8 1950
## 157 8 1951
## 158 8 1952
## 159 8 1953
## 160 8 1954
## 161 9 1935
## 162 9 1936
## 163 9 1937
## 164 9 1938
## 165 9 1939
## 166 9 1940
## 167 9 1941
## 168 9 1942
## 169 9 1943
## 170 9 1944
## 171 9 1945
## 172 9 1946
## 173 9 1947
## 174 9 1948
## 175 9 1949
## 176 9 1950
## 177 9 1951
## 178 9 1952
## 179 9 1953
## 180 9 1954
## 181 10 1935
## 182 10 1936
## 183 10 1937
## 184 10 1938
## 185 10 1939
## 186 10 1940
## 187 10 1941
## 188 10 1942
## 189 10 1943
## 190 10 1944
## 191 10 1945
## 192 10 1946
## 193 10 1947
## 194 10 1948
## 195 10 1949
## 196 10 1950
## 197 10 1951
## 198 10 1952
## 199 10 1953
## 200 10 1954
# Menampilkan ringkasan data panel
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
##
# Memeriksa apakah panel balanced
is.pbalanced(pdata)
## [1] TRUE
# 4. MENENTUKAN MODEL REGRESI DATA PANEL
# Variabel dependen : inv
# Variabel independen: value dan capital
form <- inv ~ value + capital
# 4.1 COMMON EFFECT MODEL (CEM)
cem <- plm(
formula = form,
data = pdata,
model = "pooling"
)
# Menampilkan hasil estimasi CEM
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 0.00001207 ***
## value 0.1155622 0.0058357 19.8026 < 0.00000000000000022 ***
## capital 0.2306785 0.0254758 9.0548 < 0.00000000000000022 ***
## ---
## 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: < 0.000000000000000222
# 4.2 FIXED EFFECT MODEL (FEM)
fem <- plm(
formula = form,
data = pdata,
model = "within",
effect = "individual"
)
# Menampilkan hasil estimasi FEM
summary(fem)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = form, data = pdata, effect = "individual", 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 < 0.00000000000000022 ***
## capital 0.310065 0.017355 17.8666 < 0.00000000000000022 ***
## ---
## 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: < 0.000000000000000222
# 4.3 RANDOM EFFECT MODEL (REM)
rem <- plm(
formula = form,
data = pdata,
model = "random",
effect = "individual"
)
# Menampilkan hasil estimasi REM
summary(rem)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = form, data = pdata, effect = "individual", 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 < 0.0000000000000002 ***
## capital 0.308113 0.017180 17.9339 < 0.0000000000000002 ***
## ---
## 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: < 0.000000000000000222
# 5. MELIHAT EFEK INDIVIDU PADA FEM
# Menampilkan nilai efek individu setiap perusahaan
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
# Menampilkan jenis efek yang digunakan
pdim(fem$model)
## Balanced Panel: n = 10, T = 20, N = 200
# 6. MEMBANDINGKAN KETIGA MODEL
stargazer(
cem,
fem,
rem,
type = "text",
title = "Perbandingan Estimasi Model Data Panel",
column.labels = c("CEM", "FEM", "REM"),
dep.var.labels = "Investasi (inv)",
covariate.labels = c("Value", "Capital"),
digits = 4,
model.numbers = FALSE,
omit.stat = c("f", "ser")
)
##
## Perbandingan Estimasi Model Data Panel
## =============================================
## Dependent variable:
## --------------------------------
## Investasi (inv)
## CEM FEM REM
## ---------------------------------------------
## Value 0.1156*** 0.1101*** 0.1098***
## (0.0058) (0.0119) (0.0105)
##
## Capital 0.2307*** 0.3101*** 0.3081***
## (0.0255) (0.0174) (0.0172)
##
## Constant -42.7144*** -57.8344**
## (9.5117) (28.8989)
##
## ---------------------------------------------
## Observations 200 200 200
## R2 0.8124 0.7668 0.7695
## Adjusted R2 0.8105 0.7531 0.7672
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
# 7. UJI PEMILIHAN MODEL TERBAIK
# 7.1 UJI CHOW (CEM VS FEM)
uji_chow <- pFtest(fem, cem)
print(uji_chow)
##
## F test for individual effects
##
## data: form
## F = 49.177, df1 = 9, df2 = 188, p-value < 0.00000000000000022
## alternative hypothesis: significant effects
# Hipotesis:
# H0: Common Effect Model
# H1: Fixed Effect Model
#
# Keputusan:
# p-value < 0.05 : FEM
# p-value >= 0.05: CEM
# 7.2 UJI BREUSCH-PAGAN LM (CEM VS REM)
uji_lm <- plmtest(
cem,
type = "bp",
effect = "individual"
)
print(uji_lm)
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: form
## chisq = 798.16, df = 1, p-value < 0.00000000000000022
## alternative hypothesis: significant effects
# Hipotesis:
# H0: Common Effect Model
# H1: Random Effect Model
#
# Keputusan:
# p-value < 0.05 : REM
# p-value >= 0.05: CEM
# 7.3 UJI HAUSMAN (FEM VS REM)
uji_hausman <- phtest(fem, rem)
print(uji_hausman)
##
## Hausman Test
##
## data: form
## chisq = 2.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent
# Hipotesis:
# H0: Random Effect Model
# H1: Fixed Effect Model
#
# Keputusan:
# p-value < 0.05 : FEM
# p-value >= 0.05: REM
# 8. TWO-WAY FIXED EFFECT MODEL
# Mengestimasi FEM dengan efek individu dan waktu
fem_tw <- plm(
formula = form,
data = pdata,
model = "within",
effect = "twoways"
)
# Menampilkan hasil Two-Way FEM
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 0.000000000000006653 ***
## capital 0.357916 0.022719 15.7540 < 0.00000000000000022 ***
## ---
## 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: < 0.000000000000000222
# Menguji apakah efek waktu diperlukan
uji_twoway <- pFtest(fem_tw, fem)
print(uji_twoway)
##
## F test for twoways effects
##
## data: form
## F = 1.4032, df1 = 19, df2 = 169, p-value = 0.1309
## alternative hypothesis: significant effects
# Hipotesis:
# H0: Efek individu saja sudah cukup
# H1: Efek individu dan waktu diperlukan
# 9. UJI ASUMSI MODEL
# 9.1 UJI HETEROSKEDASTISITAS
# BREUSCH-PAGAN TEST
uji_hetero <- bptest(
inv ~ value + capital + factor(firm),
data = Grunfeld
)
print(uji_hetero)
##
## studentized Breusch-Pagan test
##
## data: inv ~ value + capital + factor(firm)
## BP = 85.836, df = 11, p-value = 0.0000000000001086
# Hipotesis:
# H0: Tidak terdapat heteroskedastisitas
# H1: Terdapat heteroskedastisitas
#
# p-value < 0.05 : Heteroskedastisitas
# p-value >= 0.05: Tidak terdapat heteroskedastisitas
# 9.2 UJI AUTOKORELASI BREUSCH-GODFREY
uji_autokorelasi <- pbgtest(fem)
print(uji_autokorelasi)
##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: form
## chisq = 65.063, df = 20, p-value = 0.00000114
## alternative hypothesis: serial correlation in idiosyncratic errors
# Hipotesis:
# H0: Tidak terdapat autokorelasi
# H1: Terdapat autokorelasi
#
# p-value < 0.05 : Terdapat autokorelasi
# p-value >= 0.05: Tidak terdapat autokorelasi
# 9.3 UJI AUTOKORELASI WOOLDRIDGE
uji_wooldridge <- pwartest(fem)
print(uji_wooldridge)
##
## Wooldridge's test for serial correlation in FE panels
##
## data: fem
## F = 76.929, df1 = 1, df2 = 188, p-value = 0.000000000000001057
## alternative hypothesis: serial correlation
# Hipotesis:
# H0: Tidak terdapat autokorelasi
# H1: Terdapat autokorelasi
# 9.4 UJI CROSS-SECTIONAL DEPENDENCE
uji_dependensi <- pcdtest(
fem,
test = "cd"
)
print(uji_dependensi)
##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: inv ~ value + capital
## z = 4.6612, p-value = 0.000003144
## alternative hypothesis: cross-sectional dependence
# Hipotesis:
# H0: Tidak terdapat dependensi antarunit
# H1: Terdapat dependensi antarunit
# 9.5 UJI MULTIKOLINEARITAS
# Membentuk model regresi untuk menghitung VIF
model_vif <- lm(
inv ~ value + capital,
data = Grunfeld
)
# Menghitung nilai VIF
hasil_vif <- vif(model_vif)
print(hasil_vif)
## value capital
## 1.313773 1.313773
# Melihat korelasi antarvariabel independen
hasil_korelasi <- cor(
Grunfeld[, c("value", "capital")]
)
print(hasil_korelasi)
## value capital
## value 1.0000000 0.4887058
## capital 0.4887058 1.0000000
# Kriteria umum VIF:
# VIF < 10 : Tidak ada indikasi multikolinearitas serius
# VIF >= 10 : Indikasi multikolinearitas serius
# 10. UJI AKAR UNIT (STASIONERITAS)
# LEVIN-LIN-CHU TEST
# Menguji stasioneritas variabel investasi (inv)
uji_stasioneritas <- purtest(
pdata$inv,
test = "levinlin",
exo = "intercept",
lags = "AIC",
pmax = 4
)
## Warning in selectT(l, theTs): the time series is short
print(uji_stasioneritas)
##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$inv
## z = 3.1416, p-value = 0.9992
## alternative hypothesis: stationarity
# Hipotesis:
# H0: Data mengandung unit root (tidak stasioner)
# H1: Data tidak mengandung unit root (stasioner)
#
# p-value < 0.05 : Data stasioner
# p-value >= 0.05: Data tidak stasioner
# 11. PENANGANAN PELANGGARAN ASUMSI
# ROBUST STANDARD ERROR
# 11.1 ROBUST STANDARD ERROR CLUSTER PER UNIT
robust_cluster <- coeftest(
fem,
vcov. = vcovHC(
fem,
method = "arellano",
type = "HC1",
cluster = "group"
)
)
print(robust_cluster)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.014414 7.6398 0.000000000001076 ***
## capital 0.310065 0.050043 6.1959 0.000000003576064 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 11.2 DRISCOLL-KRAAY STANDARD ERROR
robust_dk <- coeftest(
fem,
vcov. = vcovSCC(
fem,
type = "HC1"
)
)
print(robust_dk)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.017775 6.1954 0.0000000035861808868 ***
## capital 0.310065 0.034996 8.8601 0.0000000000000006017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 12. SELESAI
cat("\n====================================\n")
##
## ====================================
cat("ANALISIS DATA PANEL TELAH SELESAI\n")
## ANALISIS DATA PANEL TELAH SELESAI
cat("Dataset: Grunfeld\n")
## Dataset: Grunfeld
cat("Jumlah observasi: 200\n")
## Jumlah observasi: 200
cat("Jumlah perusahaan: 10\n")
## Jumlah perusahaan: 10
cat("Periode: 1935-1954\n")
## Periode: 1935-1954
cat("====================================\n")
## ====================================