I will focus exclusively on 5-year panel data, especially the table 5 and Table 6 in Instrumental Variable part.

Load data from Acemoglu et al (2008)

#load library(plm)
#load library(modelsummary) or library(stargazer) - nice way of creating output/plots
#load library(wesanderson) - nice color
X5_yr_panel <- read.csv("5_yr_panel.csv", header=TRUE)

Presetting background for some plots below:

b <- list(geom_vline(xintercept = 0, color = 'red'),
          annotate("rect", alpha = .1,
                   xmin = -.5, xmax = .5, 
                   ymin = -Inf, ymax = Inf),
          geom_point(aes(y = term, x = estimate), alpha = .1, 
                     size = 5, color = 'blue'))
fe2sls <- c("FE 2SLS (4)", "FE 2SLS (5)", "FE 2SLS (7)", "FE 2SLS (8)", "FE 2SLS (9)")

Table 2

table2_m1 <- plm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "pooling",
              effect = "twoways")
table2_m2 <- plm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table2_m4 <- pgmm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch, 2) | lag(lrgdpch, 2),
               data = X5_yr_panel,
               index = c("code_numeric", "year_numeric"),
               effect = "twoways",
               model = "twosteps")
table2_m5 <- plm(fhpolrigaug ~ lag(lrgdpch),
                data=X5_yr_panel,
                index = c("code_numeric", "year_numeric"),
                model = "within",
                effect = "twoways")

#I need to remove Model 4 (GMM) as the stargazer cannot put together with other models
table2 <- list(table2_m1, table2_m2, table2_m5)

#Let's look at Table 2 - Model 4

stargazer(table2_m4, 
          title="Table 2 - Model 4 Arellano-Bond GMM", 
          dep.var.labels = "Democracy",
          covariate.labels = c("Democracy (t-1)", "log Income per Capita (t-2)"), 
          type="html", out.header=TRUE)
Table 2 - Model 4 Arellano-Bond GMM
Dependent variable:
Democracy
Democracy (t-1) -0.172***
(0.048)
log Income per Capita (t-2) 0.050
(0.133)
Observations 211
Note: p<0.1; p<0.05; p<0.01 
stargazer(table2, title="Table 2 - Model 1, 2, 5", align=TRUE,
          dep.var.labels = "Democracy",
          covariate.labels = c("Democracy (t-1)", "log Income per capita (t-1)", "Intercept"), type="html", out.header=TRUE
)
Table 2 - Model 1, 2, 5
Dependent variable:
Democracy
(1) (2) (3)
Democracy (t-1) 0.700*** 0.389***
(0.023) (0.033)
log Income per capita (t-1) 0.078*** 0.001 0.062**
(0.008) (0.025) (0.025)
Intercept -0.462***
(0.058)
Observations 988 988 1,048
R2 0.708 0.147 0.007
Adjusted R2 0.707 -0.018 -0.171
F Statistic 1,192.707*** (df = 2; 985) 71.180*** (df = 2; 827) 5.999** (df = 1; 888)
Note: p<0.1; p<0.05; p<0.01

Let’s take a look at models 1, 2, and 5 from Table 2

table_map <- c("lag(lrgdpch)"="log Income per capita (t-1)",
               "lag(fhpolrigaug)"="Democracy (t-1)")

modelplot(table2,
          coef_omit = "Intercept", 
          coef_map = table_map,
          background = b) +
  labs(x = "95% Confidence Interval of coefficients",
        y = "Coefficients name",
       caption = "Dependent variable is Democracy (Freedom House)") +
    scale_color_manual(values = wes_palette('Darjeeling1'))

Pool OLS shows that there is an effect of income on democracy.
After controlling for unit-specific time-invariant effects (i.e., conditional on the confounders/ close the backdoor path), the effects are gone (Model 2) or substantively small to neglect (Model 5).

Table 5

Table 5 Replication from Model 4 to Model 9 (excluding model 6 - Arellano-Bond GMM model).
From Model 4 to 9, Acemoglu et al (2008) used instrumental variable with two-way fixed effects. For my interest, I will replicate only this part and exclude model 6.

#Table 5 - model 4
table5_4ivA <- plm(fhpolrigaug ~ lag(lrgdpch) | lag(nsave, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table5_4ivB <- plm(lag(lrgdpch) ~ lag(nsave, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# model 5
table5_5ivA <- plm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch) | lag(fhpolrigaug) + lag(nsave, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table5_5ivB <- plm(lag(lrgdpch) ~ lag(fhpolrigaug) + lag(nsave, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 7
table5_7ivA <- plm(fhpolrigaug ~ lag(lrgdpch) + lag(laborshare) | lag(laborshare) + lag(nsave, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table5_7ivB <- plm(lag(lrgdpch) ~ lag(laborshare) + lag(nsave, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 8
table5_8ivA <- plm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch)| lag(fhpolrigaug) + lag(nsave, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table5_8ivB <- plm(lag(lrgdpch) ~ lag(fhpolrigaug) + lag(nsave, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 9
table5_9ivA <- plm(fhpolrigaug ~ lag(lrgdpch) | lag(nsave, 2) + lag(nsave, 3),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table5_9ivB <- plm(lag(lrgdpch) ~ lag(nsave, 2) + lag(nsave, 3),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

table5_sumA <- list(
                "FE 2SLS (4)"=table5_4ivA, 
                "FE 2SLS (5)"=table5_5ivA, 
                "FE 2SLS (7)"=table5_7ivA, 
                "FE 2SLS (8)"=table5_8ivA, 
                "FE 2SLS (9)"=table5_9ivA
                )

table5_sumB <- list(
                "FE 2SLS (4)"=table5_4ivB, 
                "FE 2SLS (5)"=table5_5ivB, 
                "FE 2SLS (7)"=table5_7ivB, 
                "FE 2SLS (8)"=table5_8ivB, 
                "FE 2SLS (9)"=table5_9ivB
                )

#rename coefficients for coef_map
sumAm <- c("lag(laborshare)"="Labor share (t-1)", 
           "lag(fhpolrigaug)"="Democracy (t-1)", 
           "lag(lrgdpch)"="log Income per capita (t-1)")

sumBm <- c("lag(nsave, 3)" = "Saving rate (t-3)",
           "lag(laborshare)"="Labor share (t-1)",
           "lag(fhpolrigaug)"="Democracy (t-1)",
           "lag(nsave, 2)" = "Saving rate (t-2)"
           )

#background
b <- list(geom_vline(xintercept = 0, color = 'red'),
          annotate("rect", alpha = .1,
                   xmin = -.5, xmax = .5, 
                   ymin = -Inf, ymax = Inf),
          geom_point(aes(y = term, x = estimate), alpha = .1, 
                     size = 5, color = 'blue'))

modelplot(table5_sumA, 
          coef_omit = "Intercept", 
          coef_map = sumAm,
          background = b) +
  labs(x = "95% Confidence Interval of coefficients",
        y = "Coefficients name",
       caption = "Dependent variable is Democracy (Freedom House)") +
    scale_color_manual(values = wes_palette('Darjeeling1'))

modelplot(table5_sumB, 
          coef_omit = "Intercept", 
          coef_map = sumBm,
          background = b) +
  labs(x = "95% Confidence Interval of coefficients",
        y = "Coefficients name",
       caption = "Dependent variable is lag of log Income per capita") +
    scale_color_manual(values = wes_palette('Darjeeling2'))

Table 6

Table 6 Replication, from Model 4 to 9 (excluding model 5) Here I use plm package

#Table 6 - model 4
table6_4ivA <- plm(fhpolrigaug ~ lag(lrgdpch) | lag(worldincome),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table6_4ivB <- plm(lag(lrgdpch) ~ lag(worldincome),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")


# model 5
table6_5ivA <- plm(fhpolrigaug ~ lag(fhpolrigaug) + lag(lrgdpch) | lag(fhpolrigaug) + lag(worldincome),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table6_5ivB <- plm(lag(lrgdpch) ~ lag(fhpolrigaug) + lag(worldincome),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 7
table6_7ivA <- plm(fhpolrigaug ~ lag(lrgdpch) + lag(worlddemocracy) | lag(worlddemocracy) + lag(worldincome), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table6_7ivB <- plm(lag(lrgdpch) ~ lag(worldincome) + lag(worlddemocracy), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 8
table6_8ivA <- plm(fhpolrigaug ~ lag(lrgdpch)| lag(worldincome, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table6_8ivB <- plm(lag(lrgdpch) ~ lag(worldincome, 2), 
                  data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")

# Model 9
table6_9ivA <- plm(fhpolrigaug ~ lag(lrgdpch) | lag(worldincome) + lag(worldincome, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")
table6_9ivB <- plm(lag(lrgdpch) ~ lag(worldincome) + lag(worldincome, 2),
              data = X5_yr_panel,
              index = c("code_numeric", "year_numeric"),
              model = "within",
              effect = "twoways")



#Preparing for plotting

table6_sumA <- list(
                "FE 2SLS (4)"=table6_4ivA, 
                "FE 2SLS (5)"=table6_5ivA, 
                "FE 2SLS (7)"=table6_7ivA, 
                "FE 2SLS (8)"=table6_8ivA, 
                "FE 2SLS (9)"=table6_9ivA
                )

table6_sumB <- list(
                "FE 2SLS (4)"=table6_4ivB, 
                "FE 2SLS (5)"=table6_5ivB, 
                "FE 2SLS (7)"=table6_7ivB, 
                "FE 2SLS (8)"=table6_8ivB, 
                "FE 2SLS (9)"=table6_9ivB
                )

#rename coefficients for coef_map
sumAm_table6 <- c("lag(worlddemocracy)"="Trade-weighted democracy (t-1)", 
           "lag(fhpolrigaug)"="Democracy (t-1)", 
           "lag(lrgdpch)"="Income per capita (t-1)")

sumBm_table6 <- c( "lag(worldincome, 2)"="Trade-weighted GDP (t-2)", 
                   "lag(worlddemocracy)"="Trade-weighted democracy (t-1)",
                   "lag(fhpolrigaug)"="Democracy (t-1)", 
                   "lag(worldincome)"="Trade-weighted GDP (t-1)")


#Plot function from here
modelplot(table6_sumA, 
          coef_omit = "Intercept", 
          coef_map = sumAm_table6,
          background = b) +
  labs(x = "95% Confidence Interval of coefficients",
        y = "Coefficients name",
       caption = "Dependent variable is lag of log Income per capita") +
    scale_color_manual(values = wes_palette('Darjeeling1'))

modelplot(table6_sumB, 
          coef_omit = "Intercept", 
          coef_map = sumBm_table6,
          background = b) +
  labs(x = "95% Confidence Interval of coefficients",
        y = "Coefficients name",
       caption = "Dependent variable is lag of log Income per capita") +
    scale_color_manual(values = wes_palette('Darjeeling2'))

Full estimation results

We can see all the estimation results of Table 5 and 6 here.

stargazer(table5_sumA, title="Table 5 - Panel A", align=TRUE,
          dep.var.labels = "Democracy",
          covariate.labels = c("Democracy (t-1)", "log Income per capita (t-1)",  "Labor share (t-1)"), 
          type="html", 
          out.header=TRUE,          
          model.names = FALSE,
          column.labels = fe2sls
)
Table 5 - Panel A
Dependent variable:
Democracy
FE 2SLS (4) FE 2SLS (5) FE 2SLS (7) FE 2SLS (8) FE 2SLS (9)
(1) (2) (3) (4) (5)
Democracy (t-1) 0.360*** 0.360***
(0.036) (0.036)
log Income per capita (t-1) -0.001 -0.025 -0.029 -0.025 0.028
(0.069) (0.071) (0.140) (0.071) (0.077)
Labor share (t-1) 0.281*
(0.163)
Observations 931 898 482 898 811
R2 0.008 0.123 0.005 0.123 0.003
Adjusted R2 -0.171 -0.043 -0.276 -0.043 -0.192
F Statistic 0.0001 107.386*** 3.313 107.386*** 0.129
Note: p<0.1; p<0.05; p<0.01 
stargazer(table5_sumB, title="Table 5 - Panel B", align=TRUE,
          dep.var.labels = "log Income per capita",
          covariate.labels = c("Democracy (t-1)", "Labor share (t-1)", "Saving rate (t-2)", "Saving rate (t-3)"), 
          type="html", 
          out.header=TRUE,
          model.names = FALSE,
          column.labels = fe2sls
)
Table 5 - Panel B
Dependent variable:
log Income per capita
FE 2SLS (4) FE 2SLS (5) FE 2SLS (7) FE 2SLS (8) FE 2SLS (9)
(1) (2) (3) (4) (5)
Democracy (t-1) 0.146*** 0.146***
(0.041) (0.041)
Labor share (t-1) 0.331***
(0.122)
Saving rate (t-2) 1.391*** 1.357*** 1.173*** 1.357*** 0.989***
(0.109) (0.115) (0.136) (0.115) (0.128)
Saving rate (t-3) 0.678***
(0.124)
Observations 952 902 487 902 820
R2 0.167 0.166 0.192 0.166 0.181
Adjusted R2 0.021 0.009 -0.034 0.009 0.023
F Statistic 162.469*** (df = 1; 809) 75.599*** (df = 2; 758) 45.082*** (df = 2; 380) 75.599*** (df = 2; 758) 75.988*** (df = 2; 686)
Note: p<0.1; p<0.05; p<0.01 
stargazer(table6_sumA, title="Table 6 - Panel A", align=TRUE,
          dep.var.labels = "Democracy",
          covariate.labels = c("Democracy (t-1)", "log Income per capita (t-1)", "Trade-weighted democracy (t-1)"), 
          type="html", 
          out.header=TRUE,
          model.names = FALSE,
          column.labels = fe2sls
)
Table 6 - Panel A
Dependent variable:
Democracy
FE 2SLS (4) FE 2SLS (5) FE 2SLS (7) FE 2SLS (8) FE 2SLS (9)
(1) (2) (3) (4) (5)
Democracy (t-1) 0.407***
(0.036)
log Income per capita (t-1) -0.221** -0.153* -0.224** -0.219* -0.209**
(0.110) (0.090) (0.112) (0.120) (0.104)
Trade-weighted democracy (t-1) -0.253
(0.335)
Observations 989 936 989 951 951
R2 0.004 0.118 0.004 0.005 0.005
Adjusted R2 -0.151 -0.029 -0.152 -0.156 -0.156
F Statistic 4.060** 132.852*** 4.312 3.363* 4.028**
Note: p<0.1; p<0.05; p<0.01 
stargazer(table6_sumB, title="Table 6 - Panel B", align=TRUE,
          dep.var.labels = "log Income per capita",
          covariate.labels = c("Democracy (t-1)", "Trade-weighted democracy (t-1)", "Trade-weighted GDP (t-1)", "Trade-weighted GDP (t-2)"), 
          type="html", 
          out.header=TRUE,
          model.names = FALSE,
          column.labels = fe2sls
)
Table 6 - Panel B
Dependent variable:
log Income per capita
FE 2SLS (4) FE 2SLS (5) FE 2SLS (7) FE 2SLS (8) FE 2SLS (9)
(1) (2) (3) (4) (5)
Democracy (t-1) 0.174***
(0.043)
Trade-weighted democracy (t-1) 0.301*** 0.407*** 0.301*** 0.423***
(0.046) (0.046) (0.045) (0.149)
Trade-weighted GDP (t-1) -1.846***
(0.300)
Trade-weighted GDP (t-2) 0.297*** -0.099
(0.046) (0.147)
Observations 1,040 942 1,040 989 989
R2 0.045 0.101 0.083 0.046 0.055
Adjusted R2 -0.096 -0.049 -0.053 -0.101 -0.092
F Statistic 42.349*** (df = 1; 906) 45.191*** (df = 2; 807) 40.989*** (df = 2; 905) 41.597*** (df = 1; 856) 24.985*** (df = 2; 855)
Note: p<0.1; p<0.05; p<0.01

Robust standard errors - clustered by group

#Robust SE clustered by group (Table 5)
table5_4ivAse <- table5_4ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_5ivAse <- table5_5ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_7ivAse <- table5_7ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_8ivAse <- table5_8ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_9ivAse <- table5_9ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_4ivBse <- table5_4ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_5ivBse <- table5_5ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_7ivBse <- table5_7ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_8ivBse <- table5_8ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table5_9ivBse <- table5_9ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))

table5se <- list(table5_4ivAse, table5_4ivBse, table5_5ivAse, table5_5ivBse,table5_7ivAse, table5_7ivBse, table5_8ivAse, table5_8ivBse, table5_9ivAse, table5_9ivBse)

#Robust SE clustered by group (Table 6)
table6_4ivAse <- table6_4ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_5ivAse <- table6_5ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_7ivAse <- table6_7ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_8ivAse <- table6_8ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_9ivAse <- table6_9ivA %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_4ivBse <- table6_4ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_5ivBse <- table6_5ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_7ivBse <- table6_7ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_8ivBse <- table6_8ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))
table6_9ivBse <- table6_9ivB %>% coeftest(.,vcov=vcovHC(., type="sss", cluster="group"))

table6se <- list(table6_4ivAse, table6_4ivBse, table6_5ivAse, table6_5ivBse,table6_7ivAse, table6_7ivBse, table6_8ivAse, table6_8ivBse, table6_9ivAse, table6_9ivBse)
library(fixest) #plot robust SE
## Warning: package 'fixest' was built under R version 4.2.3
coefplot(table5se)

coefplot(table6se)