第6章　パネルデータ分析

df<-read.csv("data/06_第6章/yamaguchi.csv")
df
glimpse(df)
dim(df)
null_ratio<-colMeans(is.na(df)) #列ごとの欠損値の平均を確認する
null_ratio #欠損なし
#iとtを確認する
library(dplyr)
t_size<-length(unique(df$year))
t_size
i_size<-length(unique(df$pref))
i_size
size<-t_size*i_size
size
#モデル構築
library(fixest) #固定効果モデルに使うモデル
ols<-lm(emp.rate~cap.rate,data=df)
summary(ols)
model1<- lm_robust(emp.rate~cap.rate, cluster = pref, data = df)
summary(model1)
model2<-feols(emp.rate~cap.rate|pref,data=df) #|以下で固定効果と時間効果を指定する
summary(model2,se="cluster")
model3<-feols(emp.rate~cap.rate|year,data=df)
summary(model3,se="cluster")
model4<-feols(emp.rate~cap.rate|pref+year,data=df)
summary(model4,se="cluster")
model5<-feols(emp.rate~cap.rate+age+age.hus+emp.rate.hus+urate|year,data=df)
summary(model5)
model6<-feols(emp.rate~cap.rate+age+age.hus+emp.rate.hus+urate|pref+year,data=df)
summary(model6)
#表にまとめる
msummary(list(ols,model1, model2,model3,model4,model5,model6),
         stars=TRUE)

固定効果モデルについては、{fixest}というパッケージを用いた。下記のサイトを参考にさせていただいた。https://keita43a.hatenablog.com/entry/2020/12/01/202122
表の作成は{msummury}というパッケージを用いた。下記のサイトを参考にさせていただいた。 https://keita43a.hatenablog.com/entry/2020/12/01/202122

第7章　操作変数法

df<-read.csv("data/07_第7章/ipehd_qje2009_master.csv")
df
glimpse(df)
dim(df)
null_ratio<-colMeans(is.na(df)) #列ごとの欠損値の平均を確認する
null_ratio #inctaxpcが5%欠損している


# IV法 ---------------------------------------------------------------------
library(AER) #IV法のためのパッケージ
#ols
model1<-lm(f_rw~f_prot,data=df)
summary(model1)
#IV
model2<-ivreg(f_rw~f_prot+f_young+f_jew+
                f_fem+f_ortsgeb+f_pruss+
                hhsize+lnpop+gpop+
                f_blind+f_deaf+f_dumb+f_miss|
                kmwittenberg+f_young+f_jew+
                f_fem+f_ortsgeb+f_pruss+
                hhsize+lnpop+gpop+
                f_blind+f_deaf+f_dumb+f_miss,data=df
                )
summary(model2)
msummary(list(model1,model2),stars = TRUE,title = "Table6-1 2SLS Estimation \n of the effect of Protestantism ")

#タイガーストックのF検定
F_model<-lm(f_prot~kmwittenberg+f_young+f_jew+
              f_fem+f_ortsgeb+f_pruss+
              hhsize+lnpop+gpop+
              f_blind+f_deaf+f_dumb+f_miss,data=df)
summary(F_model)
# 過剰識別検定
residuals <- residuals(model2)
df<-df |> mutate(residual=residuals)
model_residual<-lm(residual~kmwittenberg+f_young+f_jew+
                      f_fem+f_ortsgeb+f_pruss+
                     hhsize+lnpop+gpop+
                      f_blind+f_deaf+f_dumb+f_miss,data=df)
library(lmtest)
wald_test_result <- waldtest(model_residual)
wald_test_result