Dynamic models

Links to sources:

https://cran.r-project.org/web/packages/pdynmc/vignettes/pdynmc-introLong.pdf

https://cran.r-project.org/web/packages/pdynmc/pdynmc.pdf

The following research is conducted by Maria R. Koldasheva and Nikolai A. Popov.

Libraries

library(dplyr) 
library(psych)
library(readxl)
library(writexl)
library(kableExtra)

# Econometrics
library(tidyverse)
library(plm) # panel models
library(sandwich) #covariance matrixes
library(lmtest) # tests
library(xtable)# latex tables
library(stargazer) # latex regression tables
library(ggpubr) # correlation test


# Dynamic models
library(pdynmc)

Downloading the data

# Downloading
long_data <- read_csv('C:/Users/Kolya/Documents/New_studies/Diploma/Data/Data_use/Long_united_data.csv',
                      show_col_types = FALSE)

long_data

Model construction

Specifications, Panel models

Theil index specification

Theil <- log(GRP) ~ lag(log(GRP)) + lag(Div) + Edu + Pop + lag(Inv) + lag(FDI) + log(Imp)

EM and IM specification

Margins <- log(GRP) ~ lag(log(GRP)) + lag(EM) + lag(IM) + Edu + Pop + lag(Inv) +lag(FDI) + log(Imp)

Manually making logarithms for “pdynmc”

# new dataset
long_data_log = data.frame(long_data)

# Check if memory addresses are the same
tracemem(long_data_log) == tracemem(long_data)
# FALSE expected

# logarithm of appropriate columns
long_data_log$GRP=log(long_data_log$GRP)
long_data_log$Imp=log(long_data_log$Imp)
long_data_log$GRPL=log(long_data_log$GRPL)
colnames(long_data_log)

Model setting

Estimation: system GMM

Theil

Theil <- log(GRP) ~ lag(log(GRP)) + lag(Div) + Edu + Pop + lag(Inv) + lag(FDI) + log(Imp)

# Theil
Theil_d <- pdynmc(dat = long_data_log, varname.i = "region", varname.t = "year",
use.mc.diff = TRUE, use.mc.lev = FALSE, use.mc.nonlin = FALSE, inst.stata = TRUE,
include.y = TRUE, varname.y = "GRP", lagTerms.y = 1,
fur.con = TRUE, fur.con.diff = TRUE, fur.con.lev = FALSE,
varname.reg.fur = c("DivL", "Edu", "Pop", "InvL", "FDIL", "Imp"),
lagTerms.reg.fur = c(0,0,0,0,0,0),
include.dum = TRUE, dum.diff = TRUE, dum.lev = FALSE, varname.dum = "year",
w.mat.stata = TRUE, std.err = "dbl.corrected", # Hwang et al. (2021)
estimation = "onestep", opt.meth = "none")

## tracemem[0x0000016abbc10140 -> 0x0000016abbc15000]: pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous> 
## tracemem[0x0000016abbc15000 -> 0x0000016abbc135c0]: $<-.data.frame $<- pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous>

# w.mat = "zero.cov", 

Theil_d_sum <- summary(Theil_d)
Theil_d_sum

## 
## Dynamic linear panel estimation (onestep)
## Estimation steps: 1
## 
## Coefficients:
##           Estimate    Std.Err z-value Pr(>|z|)    
## L1.GRP   0.2351833  0.5351577   0.439  0.66066    
## L0.DivL -0.0024502  0.0160275  -0.153  0.87840    
## L0.Edu  -0.0001886  0.0004601  -0.410  0.68181    
## L0.Pop  -0.0021090  0.0022362  -0.943  0.34568    
## L0.InvL  0.2002416  0.2001088   1.001  0.31683    
## L0.FDIL -0.0464920  0.0634722  -0.732  0.46417    
## L0.Imp  -0.0102304  0.0153727  -0.665  0.50605    
## 2016     0.0788630  0.0251630   3.134  0.00172 ** 
## 2017     0.1028654  0.0278625   3.692  0.00022 ***
## 2018     0.1564529  0.0484014   3.232  0.00123 ** 
## 2019     0.1725278  0.0722114   2.389  0.01689 *  
## 2020     0.1372527  0.0825678   1.662  0.09651 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  26 total instruments are employed to estimate 12 parameters
##  15 linear (DIF) 
##  6 further controls (DIF) 
##  5 time dummies (DIF) 
##  
## J-Test (overid restrictions):  19.89 with 14 DF, pvalue: 0.1337
## F-Statistic (slope coeff):  0.65 with 7 DF, pvalue: 0.9987
## F-Statistic (time dummies):  15.66 with 5 DF, pvalue: 0.0079

J-Test: null hypothesis is not rejected - > our set of IVs is jointly zero F-Statistic (slope coeff): null hypothesis is not rejected - > our betas are jointly zero (bad model) F-Statistic (time dummies): null hypothesis is rejected - > our dummies betas are not jointly zero (good dummies model)

Margins

Margins <- log(GRP) ~ lag(log(GRP)) + lag(EM) + lag(IM) + Edu + Pop + lag(Inv) +lag(FDI)

# Margins
Margins_d <- pdynmc(dat = long_data_log, varname.i = "region", varname.t = "year",
use.mc.diff = TRUE, use.mc.lev = FALSE, use.mc.nonlin = FALSE, inst.stata = TRUE,
include.y = TRUE, varname.y = "GRP", lagTerms.y = 1,
fur.con = TRUE, fur.con.diff = TRUE, fur.con.lev = FALSE,
varname.reg.fur = c("EML", "IML", "Edu", "Pop", "InvL", "FDIL", "Imp"),
lagTerms.reg.fur = c(0,0,0,0,0,0,0),
include.dum = TRUE, dum.diff = TRUE, dum.lev = FALSE, varname.dum = "year",
w.mat.stata = TRUE, std.err = "dbl.corrected", # Hwang et al. (2021)
estimation = "onestep", opt.meth = "none")

## tracemem[0x0000016abbc10140 -> 0x0000016abbc08580]: pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous> 
## tracemem[0x0000016abbc08580 -> 0x0000016abbc09300]: $<-.data.frame $<- pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous>

# w.mat = "zero.cov",

Margins_d_sum <- summary(Margins_d)
Margins_d_sum

## 
## Dynamic linear panel estimation (onestep)
## Estimation steps: 1
## 
## Coefficients:
##           Estimate    Std.Err z-value Pr(>|z|)    
## L1.GRP   0.2659010  0.4944643   0.538   0.5906    
## L0.EML  -0.0141927  0.0264852  -0.536   0.5920    
## L0.IML  -0.0001670  0.0014348  -0.116   0.9076    
## L0.Edu  -0.0002269  0.0004417  -0.514   0.6072    
## L0.Pop  -0.0025740  0.0020895  -1.232   0.2180    
## L0.InvL  0.2122109  0.1811832   1.171   0.2416    
## L0.FDIL -0.0368087  0.0714718  -0.515   0.6066    
## L0.Imp  -0.0088937  0.0146467  -0.607   0.5438    
## 2016     0.0759298  0.0190095   3.994    6e-05 ***
## 2017     0.0960856  0.0362373   2.652   0.0080 ** 
## 2018     0.1480032  0.0581099   2.547   0.0109 *  
## 2019     0.1611877  0.0839598   1.920   0.0549 .  
## 2020     0.1234621  0.0957047   1.290   0.1971    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  27 total instruments are employed to estimate 13 parameters
##  15 linear (DIF) 
##  7 further controls (DIF) 
##  5 time dummies (DIF) 
##  
## J-Test (overid restrictions):  19.73 with 14 DF, pvalue: 0.1389
## F-Statistic (slope coeff):  0.77 with 8 DF, pvalue: 0.9993
## F-Statistic (time dummies):  16.17 with 5 DF, pvalue: 0.0064

J-Test: null hypothesis is not rejected - > our set of IVs is jointly zero F-Statistic (slope coeff): null hypothesis is not rejected - > our betas are jointly zero (bad model) F-Statistic (time dummies): null hypothesis is rejected - > our dummies betas are not jointly zero (good dummies model) # United table creating

# Table for Theil
united_table_dynamic_panels_Theil <-
cbind(Theil_d_sum$coefficients[c(1:(nrow(Theil_d_sum$coefficients)-5)), c(1, 2, 4)])


# Table for Margins
united_table_dynamic_panels_Margins <- 
cbind(Margins_d_sum$coefficients[c(1:(nrow(Margins_d_sum$coefficients)-5)), c(1, 2, 4)])

# matrixes into dataframe
united_table_dynamic_panels_Theil <- as.data.frame(united_table_dynamic_panels_Theil)

united_table_dynamic_panels_Margins <- as.data.frame(united_table_dynamic_panels_Margins)

# renaming
names(united_table_dynamic_panels_Theil) <- c('Est_Theil', 'SE_Theil', 'Pvalue_Theil')

names(united_table_dynamic_panels_Margins) <- c('Est_Margins', 'SE_Margins', 'Pvalue_Margins')

united_table_dynamic_panels_Theil

united_table_dynamic_panels_Margins

Table

library(modelsummary)

## 
## Присоединяю пакет: 'modelsummary'

## Следующий объект скрыт от 'package:psych':
## 
##     SD

ti_Theil <- data.frame(
  term = c('L1.GRP', 'L1.Div', 'L0.Edu', 'L0.Pop', 'L1.Inv', 'L1.FDI', 'L0.Imp'),
  estimate = united_table_dynamic_panels_Theil$Est_Theil,
  std.error = united_table_dynamic_panels_Theil$SE_Theil,
  p.value = united_table_dynamic_panels_Theil$Pvalue_Theil)

ti_Margins <- data.frame(
  term = c('L1.GRP', 'L1.IM', 'L1.EM', 'L0.Edu', 'L0.Pop', 'L1.Inv', 'L1.FDI', 'L0.Imp'),
  estimate = united_table_dynamic_panels_Margins$Est_Margins,
  std.error = united_table_dynamic_panels_Margins$SE_Margins,
  p.value = united_table_dynamic_panels_Margins$Pvalue_Margins)

gl <- data.frame(
  Observations = "546",
  Model = "System GMM")

mod_Theil <- list(tidy = ti_Theil, glance = gl)
mod_Margins <- list(tidy = ti_Margins, glance = gl)

# Model class
class(mod_Theil) <- "modelsummary_list"
class(mod_Margins) <- "modelsummary_list"


# created named list
models <- list()
models[['log(GRP) (1)']] <- mod_Theil
models[['log(GRP) (2)']] <- mod_Margins

# Model
modelsummary(models, stars = T, title = 'Dynamic model', fmt = 4, coef_map = c('L1.GRP', 'L1.Div',
                                                                               'L1.EM', 'L1.IM',
                                                                               'L0.Edu', 'L0.Pop',
                                                                               'L1.Inv','L1.FDI',
                                                                               'L0.Imp'))

Dynamic model

	log(GRP) (1)	log(GRP) (2)
L1.GRP	0.2352	0.2659
	(0.5352)	(0.4945)
L1.Div	-0.0025
	(0.0160)
L1.EM		-0.0002
		(0.0014)
L1.IM		-0.0142
		(0.0265)
L0.Edu	-0.0002	-0.0002
	(0.0005)	(0.0004)
L0.Pop	-0.0021	-0.0026
	(0.0022)	(0.0021)
L1.Inv	0.2002	0.2122
	(0.2001)	(0.1812)
L1.FDI	-0.0465	-0.0368
	(0.0635)	(0.0715)
L0.Imp	-0.0102	-0.0089
	(0.0154)	(0.0146)
Observations	546	546
Model	System GMM	System GMM
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

#output = 'latex'
# modelsummary(mod)

SAME CODE, BUT FOR ROBUST DATASET (NO MSK, SPB)

Downloading the data, robust

# Downloading
long_data_robust <- read_csv(
  'C:/Users/Kolya/Documents/New_studies/Diploma/Data/Data_use/Long_united_data_robust.csv',
  show_col_types = FALSE)

long_data_robust

Model construction, robust

Manually making logarithms for “pdynmc”, robust

# new dataset
long_data_robust_log = data.frame(long_data_robust)

# Check if memory addresses are the same
tracemem(long_data_robust_log) == tracemem(long_data_robust)
# FALSE expected

# logarithm of appropriate columns
long_data_robust_log$GRP=log(long_data_robust_log$GRP)
long_data_robust_log$Imp=log(long_data_robust_log$Imp)
long_data_robust_log$GRPL=log(long_data_robust_log$GRPL)
colnames(long_data_robust_log)

Model setting, robust

Theil_robust, robust

Theil_robust <- log(GRP) ~ lag(log(GRP)) + lag(Div) + Edu + Pop + lag(Inv) + lag(FDI) + log(Imp)

# Theil_robust
Theil_robust_d <- pdynmc(dat = long_data_robust_log, varname.i = "region", varname.t = "year",
use.mc.diff = TRUE, use.mc.lev = FALSE, use.mc.nonlin = FALSE, inst.stata = TRUE,
include.y = TRUE, varname.y = "GRP", lagTerms.y = 1,
fur.con = TRUE, fur.con.diff = TRUE, fur.con.lev = FALSE,
varname.reg.fur = c("DivL", "Edu", "Pop", "InvL", "FDIL", "Imp"),
lagTerms.reg.fur = c(0,0,0,0,0,0),
include.dum = TRUE, dum.diff = TRUE, dum.lev = FALSE, varname.dum = "year",
w.mat.stata = TRUE, std.err = "dbl.corrected", # Hwang et al. (2021)
estimation = "onestep", opt.meth = "none")

## tracemem[0x0000016abbc0bb80 -> 0x0000016abbc0bc40]: pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous> 
## tracemem[0x0000016abbc0bc40 -> 0x0000016abbc0ca80]: $<-.data.frame $<- pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous>

# w.mat = "zero.cov", 

Theil_robust_d_sum <- summary(Theil_robust_d)
Theil_robust_d_sum

## 
## Dynamic linear panel estimation (onestep)
## Estimation steps: 1
## 
## Coefficients:
##           Estimate    Std.Err z-value Pr(>|z|)    
## L1.GRP   0.2449851  0.5455467   0.449  0.65343    
## L0.DivL -0.0049161  0.0167257  -0.294  0.76876    
## L0.Edu  -0.0002579  0.0004903  -0.526  0.59889    
## L0.Pop  -0.0008184  0.0006804  -1.203  0.22898    
## L0.InvL  0.2187808  0.2039544   1.073  0.28327    
## L0.FDIL -0.0469446  0.0651843  -0.720  0.47152    
## L0.Imp  -0.0101418  0.0163275  -0.621  0.53460    
## 2016     0.0779362  0.0271762   2.868  0.00413 ** 
## 2017     0.1018038  0.0282300   3.606  0.00031 ***
## 2018     0.1544932  0.0494634   3.123  0.00179 ** 
## 2019     0.1701953  0.0732561   2.323  0.02018 *  
## 2020     0.1344938  0.0834511   1.612  0.10696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  26 total instruments are employed to estimate 12 parameters
##  15 linear (DIF) 
##  6 further controls (DIF) 
##  5 time dummies (DIF) 
##  
## J-Test (overid restrictions):  20.72 with 14 DF, pvalue: 0.1091
## F-Statistic (slope coeff):  0.6 with 7 DF, pvalue: 0.999
## F-Statistic (time dummies):  13.67 with 5 DF, pvalue: 0.0178

J-Test: null hypothesis is not rejected - > our set of IVs is jointly zero F-Statistic (slope coeff): null hypothesis is not rejected - > our betas are jointly zero (bad model) F-Statistic (time dummies): null hypothesis is rejected - > our dummies betas are not jointly zero (good dummies model)

Margins_robust

Margins_robust <- log(GRP) ~ lag(log(GRP)) + lag(EM) + lag(IM) + Edu + Pop + lag(Inv) +lag(FDI)

# Margins_robust
Margins_robust_d <- pdynmc(dat = long_data_robust_log, varname.i = "region", varname.t = "year",
use.mc.diff = TRUE, use.mc.lev = FALSE, use.mc.nonlin = FALSE, inst.stata = TRUE,
include.y = TRUE, varname.y = "GRP", lagTerms.y = 1,
fur.con = TRUE, fur.con.diff = TRUE, fur.con.lev = FALSE,
varname.reg.fur = c("EML", "IML", "Edu", "Pop", "InvL", "FDIL", 'L0.Imp'),
lagTerms.reg.fur = c(0,0,0,0,0,0,0),
include.dum = TRUE, dum.diff = TRUE, dum.lev = FALSE, varname.dum = "year",
w.mat.stata = TRUE, std.err = "dbl.corrected", # Hwang et al. (2021)
estimation = "onestep", opt.meth = "none")

## tracemem[0x0000016abbc0bb80 -> 0x0000016abbc08b80]: pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous> 
## tracemem[0x0000016abbc08b80 -> 0x0000016abbc08580]: $<-.data.frame $<- pdynmc eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers process_file <Anonymous> <Anonymous>

## Warning in `[<-.data.frame`(`*tmp*`, , varname.reg.estParam.fur, value = list(:
## предложено 532 переменные для замены 7 переменных

# w.mat = "zero.cov",

Margins_robust_d_sum <- summary(Margins_robust_d)
Margins_robust_d_sum

## 
## Dynamic linear panel estimation (onestep)
## Estimation steps: 1
## 
## Coefficients:
##            Estimate   Std.Err z-value Pr(>|z|)    
## L1.GRP     0.172093  0.357110   0.482  0.62981    
## L0.EML    -0.018148  0.004711  -3.852  0.00012 ***
## L0.IML    -0.075302  0.017611  -4.276    2e-05 ***
## L0.Edu     0.004280  0.002587   1.655  0.09792 .  
## L0.Pop    -0.062490  0.025074  -2.492  0.01270 *  
## L0.InvL   -0.023595  0.005798  -4.070    5e-05 ***
## L0.FDIL   -0.014612  0.006374  -2.293  0.02185 *  
## L0.L0.Imp -0.025875  0.005959  -4.343    1e-05 ***
## 2016       0.021200  0.021516   0.985  0.32462    
## 2017      -0.029708  0.015084  -1.970  0.04884 *  
## 2018       0.044911  0.009638   4.660  < 2e-16 ***
## 2019       0.042518  0.019248   2.209  0.02717 *  
## 2020       0.027102  0.026553   1.021  0.30725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  27 total instruments are employed to estimate 13 parameters
##  15 linear (DIF) 
##  7 further controls (DIF) 
##  5 time dummies (DIF) 
##  
## J-Test (overid restrictions):  22.08 with 14 DF, pvalue: 0.0769
## F-Statistic (slope coeff):  106.35 with 8 DF, pvalue: <0.001
## F-Statistic (time dummies):  53.27 with 5 DF, pvalue: <0.001

J-Test: null hypothesis is rejected - > our set of IVs is not jointly zero F-Statistic (slope coeff): null hypothesis is rejected - > our betas are not jointly zero (good model) F-Statistic (time dummies): null hypothesis is rejected - > our dummies betas are not jointly zero (good dummies model) # United table creating, robust

# Table for Theil_robust
united_table_dynamic_panels_Theil_robust <-
cbind(Theil_robust_d_sum$coefficients[c(1:(nrow(Theil_robust_d_sum$coefficients)-5)), c(1, 2, 4)])


# Table for Margins_robust
united_table_dynamic_panels_Margins_robust <- 
cbind(Margins_robust_d_sum$coefficients[c(1:(nrow(Margins_robust_d_sum$coefficients)-5)), c(1, 2, 4)])

# matrixes into dataframe
united_table_dynamic_panels_Theil_robust <- as.data.frame(united_table_dynamic_panels_Theil_robust)

united_table_dynamic_panels_Margins_robust <- as.data.frame(united_table_dynamic_panels_Margins_robust)

# renaming
names(united_table_dynamic_panels_Theil_robust) <- c('Est_Theil_robust', 'SE_Theil_robust', 'Pvalue_Theil_robust')

names(united_table_dynamic_panels_Margins_robust) <- c('Est_Margins_robust', 'SE_Margins_robust', 'Pvalue_Margins_robust')

united_table_dynamic_panels_Theil_robust

united_table_dynamic_panels_Margins_robust

Table, robust

library(modelsummary)
ti_Theil_robust <- data.frame(
  term = c('L1.GRP', 'L1.Div', 'L0.Edu', 'L0.Pop', 'L1.Inv', 'L1.FDI', 'L0.Imp'),
  estimate = united_table_dynamic_panels_Theil_robust$Est_Theil_robust,
  std.error = united_table_dynamic_panels_Theil_robust$SE_Theil_robust,
  p.value = united_table_dynamic_panels_Theil_robust$Pvalue_Theil_robust)

ti_Margins_robust <- data.frame(
  term = c('L1.GRP', 'L1.IM', 'L1.EM', 'L0.Edu', 'L0.Pop', 'L1.Inv', 'L1.FDI', 'L0.Imp'),
  estimate = united_table_dynamic_panels_Margins_robust$Est_Margins_robust,
  std.error = united_table_dynamic_panels_Margins_robust$SE_Margins_robust,
  p.value = united_table_dynamic_panels_Margins_robust$Pvalue_Margins_robust)

gl <- data.frame(
  Observations = "532",
  Model = "System GMM")

mod_Theil_robust <- list(tidy = ti_Theil_robust, glance = gl)
mod_Margins_robust <- list(tidy = ti_Margins_robust, glance = gl)

# Model class
class(mod_Theil_robust) <- "modelsummary_list"
class(mod_Margins_robust) <- "modelsummary_list"


# created named list
models <- list()
models[['log(GRP) (1)']] <- mod_Theil_robust
models[['log(GRP) (2)']] <- mod_Margins_robust

# Model
modelsummary(models, stars = T, title = 'Dynamic model', fmt = 4, coef_map = c('L1.GRP', 'L1.Div',
                                                                               'L1.EM', 'L1.IM',
                                                                               'L0.Edu', 'L0.Pop',
                                                                               'L1.Inv','L1.FDI',
                                                                               'L0.Imp'))

Dynamic model

	log(GRP) (1)	log(GRP) (2)
L1.GRP	0.2450	0.1721
	(0.5455)	(0.3571)
L1.Div	-0.0049
	(0.0167)
L1.EM		-0.0753***
		(0.0176)
L1.IM		-0.0181***
		(0.0047)
L0.Edu	-0.0003	0.0043+
	(0.0005)	(0.0026)
L0.Pop	-0.0008	-0.0625*
	(0.0007)	(0.0251)
L1.Inv	0.2188	-0.0236***
	(0.2040)	(0.0058)
L1.FDI	-0.0469	-0.0146*
	(0.0652)	(0.0064)
L0.Imp	-0.0101	-0.0259***
	(0.0163)	(0.0060)
Observations	532	532
Model	System GMM	System GMM
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

#output = 'latex'
# modelsummary(mod)

Final table

# created named list
models <- list()
models[['log(GRP) (1)']] <- mod_Theil
models[['log(GRP) (2)']] <- mod_Theil_robust
models[['log(GRP) (3)']] <- mod_Margins
models[['log(GRP) (4)']] <- mod_Margins_robust
# Model
modelsummary(models, stars = T, title = 'Dynamic models', fmt = 4, coef_map = c('L1.GRP', 'L1.Div',
                                                                               'L1.EM', 'L1.IM',
                                                                               'L0.Edu', 'L0.Pop',
                                                                               'L1.Inv','L1.FDI',
                                                                               'L0.Imp'))

Dynamic models

	log(GRP) (1)	log(GRP) (2)	log(GRP) (3)	log(GRP) (4)
L1.GRP	0.2352	0.2450	0.2659	0.1721
	(0.5352)	(0.5455)	(0.4945)	(0.3571)
L1.Div	-0.0025	-0.0049
	(0.0160)	(0.0167)
L1.EM			-0.0002	-0.0753***
			(0.0014)	(0.0176)
L1.IM			-0.0142	-0.0181***
			(0.0265)	(0.0047)
L0.Edu	-0.0002	-0.0003	-0.0002	0.0043+
	(0.0005)	(0.0005)	(0.0004)	(0.0026)
L0.Pop	-0.0021	-0.0008	-0.0026	-0.0625*
	(0.0022)	(0.0007)	(0.0021)	(0.0251)
L1.Inv	0.2002	0.2188	0.2122	-0.0236***
	(0.2001)	(0.2040)	(0.1812)	(0.0058)
L1.FDI	-0.0465	-0.0469	-0.0368	-0.0146*
	(0.0635)	(0.0652)	(0.0715)	(0.0064)
L0.Imp	-0.0102	-0.0101	-0.0089	-0.0259***
	(0.0154)	(0.0163)	(0.0146)	(0.0060)
Observations	546	532	546	532
Model	System GMM	System GMM	System GMM	System GMM
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001