# Load packages
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(lubridate)
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(fixest)
library(sampleSelection)
## Loading required package: maxLik
## Loading required package: miscTools
##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/
library(modelsummary)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
financial <- read_excel("GISC Financial.xlsx")|>
filter(`(costat) Active/Inactive Status Marker`== "A", `(fic) Current ISO Country Code - Incorporation` == "USA") |>
select(subindustry_code = `(gsubind) GIC Sub-Industries`,
date = `(datadate) Data Date`,
gvkey = `(gvkey) Global Company Key`,
name = `(conm) Company Name`,
total_asset = `(at) Assets - Total`,
current_debt = `(dlc) Debt in Current Liabilities - Total`,
long_term_debt = `(dltt) Long-Term Debt - Total`,
net_income = `(ni) Net Income (Loss)`,
capex = `(capx) Capital Expenditures`) |>
mutate(total_debt = current_debt + long_term_debt,
ROE = 100*net_income/(total_asset-total_debt),
year = year(date),
gvkey = as.character(gvkey)) |>
select(-current_debt, -long_term_debt, -net_income, -date)
ghg_intensity = read_csv("GHG intensity.csv") |>
mutate(year = year(periodenddate), gvkey = as.character(gvkey), ghg_intensity = log(di_319405)) |>
select(gvkey, year, ghg_intensity)
## Rows: 1926 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): gvkey, ticker, companyname, country
## dbl (4): companyid, institutionid, fiscalyear, di_319405
## date (1): periodenddate
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
esg_score <- read_csv("ESG Global Score.csv") |>
filter(
csascoretypename == "Raw",
aspectname %in% c(
"S&P Global ESG Score",
"Environmental Dimension",
"Social Dimension",
"Economic Governance Dimension"
)
) |>
transmute(
scoredate,
year = lubridate::year(scoredate),
gvkey = as.character(gvkey),
score_type = aspectname,
scorevalue
) |>
arrange(gvkey, year, score_type, scoredate) |>
group_by(gvkey, year, score_type) |>
slice_tail(n = 1) |>
ungroup() |>
select(-scoredate) |>
pivot_wider(
names_from = score_type,
values_from = scorevalue
) |>
rename(
ESG_composite_score = `S&P Global ESG Score`,
E_score = `Environmental Dimension`,
S_score = `Social Dimension`,
G_score = `Economic Governance Dimension`
)
## Rows: 5624 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (6): csascoretypename, scoretype, aspectname, ticker, companyname, country
## dbl (3): institutionid, scorevalue, gvkey
## date (1): scoredate
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gas_production <- read_csv("Oil and gas production.csv") |>
select(year = datadate, gvkey, natural_gas_production = ogpngq, oil_production = ogpoilq, country = fic) |>
filter(country == "USA", !is.na(natural_gas_production), !is.na(oil_production)) |>
mutate(year = year(year))|>
group_by(gvkey, year) |>
mutate(annual_total_gas_production = sum(natural_gas_production,na.rm = TRUE),
annual_total_oil_production = sum(oil_production, na.rm = TRUE),
gvkey = as.character(gvkey)) |>
slice(1) |>
mutate(
gas_boe = annual_total_gas_production / 6000,
total_boe = annual_total_oil_production + gas_boe,
gas_percent = 100*gas_boe / total_boe
) |>
select(year, gvkey, gas_percent)
## Rows: 18531 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): datafmt, indfmt, consol, gvkey, conm, tic, fic, conml
## dbl (3): gsubind, ogpngq, ogpoilq
## date (1): datadate
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
brent_oil <- read_csv("DCOILBRENTEU.csv") |>
mutate(year = year(observation_date)) |>
rename("oil_price" = DCOILBRENTEU) |>
filter(!is.na(oil_price)) |>
group_by (year) |>
summarize(oil_price = mean(oil_price))
## Rows: 4221 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (1): DCOILBRENTEU
## date (1): observation_date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
joint_data <- financial |>
left_join(esg_score, by = join_by(gvkey, year)) |>
left_join(ghg_intensity, by = join_by(gvkey, year)) |>
left_join(gas_production, by = join_by(gvkey, year)) |>
left_join(brent_oil, by = join_by(year))
#Regression with Heckman Selection Model using ESG scores
# Create selection indicators
df2 <- joint_data |>
arrange(gvkey, year) |>
group_by(gvkey) |>
mutate(
ESG_obs = if_else(!is.na(ESG_composite_score), 1, 0),
E_obs = if_else(!is.na(E_score), 1, 0),
S_obs = if_else(!is.na(S_score), 1, 0),
G_obs = if_else(!is.na(G_score), 1, 0)
) |>
ungroup() |>
mutate(
log_capex = log(capex + 1),
subindustry_code = as.factor(subindustry_code),
year = as.integer(year)
) |>
filter(!is.na(log_capex), !is.na(year), !is.na(subindustry_code))
# 3) Function to create IMR from first-stage probit
add_imr <- function(data, sel_var, imr_name) {
f <- as.formula(
paste0(sel_var, " ~ log_capex + factor(year) + subindustry_code")
)
sel_mod <- glm(f, data = data, family = binomial(link = "probit"))
xb <- predict(sel_mod, type = "link")
phi <- dnorm(xb)
Phi <- pnorm(xb)
data[[imr_name]] <- ifelse(data[[sel_var]] == 1, phi / pmax(Phi, 1e-8), NA_real_)
data
}
# 4) Add IMRs for each score
df2 <- add_imr(df2, "ESG_obs", "imr_ESG")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
df2 <- add_imr(df2, "E_obs", "imr_E")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
df2 <- add_imr(df2, "S_obs", "imr_S")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
df2 <- add_imr(df2, "G_obs", "imr_G")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# 5) Function to run 4 Heckman-style regressions for ESG variables
run_heckman_set <- function(data, x_var, sel_var, imr_var, label) {
reg_df <- data |>
filter(
.data[[sel_var]] == 1,
!is.na(ROE),
!is.na(.data[[x_var]]),
!is.na(ghg_intensity),
!is.na(gas_percent),
!is.na(oil_price),
!is.na(.data[[imr_var]])
)
f1 <- as.formula(
paste0("ROE ~ ", x_var, " + ", imr_var)
)
f2 <- as.formula(
paste0("ROE ~ ", x_var, " + gas_percent + ", imr_var)
)
f3 <- as.formula(
paste0("ROE ~ ", x_var, " + gas_percent + ghg_intensity + ", imr_var )
)
f4 <- as.formula(
paste0("ROE ~ ", x_var, " + gas_percent + ghg_intensity + oil_price +", imr_var)
)
f5 <- as.formula(
paste0("ROE ~ ", x_var, " + gas_percent + ghg_intensity + ", imr_var, " | gvkey + year")
)
list(
setNames(list(feols(f1, data = reg_df, vcov = ~ gvkey)), paste0(label, " (1)"))[[1]],
setNames(list(feols(f2, data = reg_df, vcov = ~ gvkey)), paste0(label, " (2)"))[[1]],
setNames(list(feols(f3, data = reg_df, vcov = ~ gvkey)), paste0(label, " (3)"))[[1]],
setNames(list(feols(f4, data = reg_df, vcov = ~ gvkey)), paste0(label, " (4)"))[[1]],
setNames(list(feols(f5, data = reg_df, vcov = ~ gvkey)), paste0(label, " Firm+Year FE (5)"))[[1]]
)
}
# 6) Run model sets
mods2_ESG <- run_heckman_set(df2, "ESG_composite_score", "ESG_obs", "imr_ESG", "Composite ESG")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods2_E <- run_heckman_set(df2, "E_score", "E_obs", "imr_E", "E score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods2_S <- run_heckman_set(df2, "S_score", "S_obs", "imr_S", "S score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods2_G <- run_heckman_set(df2, "G_score", "G_obs", "imr_G", "G score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
output_table2 <- modelsummary(
c(mods2_ESG, mods2_E, mods2_S, mods2_G),
stars = TRUE,
statistic = "({std.error})",
gof_omit = "AIC|BIC|Log.Lik|RMSE|Adj|Within|Pseudo"
)
output_table2
| |
(1) |
(2) |
(3) |
(4) |
(5) |
(6) |
(7) |
(8) |
(9) |
(10) |
(11) |
(12) |
(13) |
(14) |
(15) |
(16) |
(17) |
(18) |
(19) |
(20) |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
| (Intercept) |
13.963 |
4.221 |
12.209 |
-58.366+ |
|
22.350 |
12.758 |
19.755 |
-56.196+ |
|
28.349 |
18.187 |
27.865 |
-53.901+ |
|
0.434 |
-8.006 |
1.361 |
-63.668 |
|
|
(27.240) |
(30.819) |
(28.297) |
(32.110) |
|
(25.994) |
(27.943) |
(27.363) |
(29.374) |
|
(25.007) |
(26.567) |
(26.922) |
(28.287) |
|
(22.146) |
(27.128) |
(24.490) |
(37.935) |
|
| ESG_composite_score |
0.701 |
0.710 |
0.928 |
0.415 |
-2.794+ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.657) |
(0.691) |
(0.738) |
(0.695) |
(1.291) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| imr_ESG |
-23.298 |
-20.978 |
-10.427 |
-4.787 |
65.824 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(17.454) |
(17.683) |
(17.409) |
(14.444) |
(39.220) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| gas_percent |
|
0.199 |
0.233 |
0.149 |
2.164 |
|
0.197 |
0.223 |
0.134 |
3.066 |
|
0.206 |
0.237 |
0.140 |
2.380 |
|
0.184 |
0.212 |
0.147 |
2.128 |
|
|
(0.255) |
(0.262) |
(0.227) |
(1.871) |
|
(0.221) |
(0.224) |
(0.201) |
(2.493) |
|
(0.247) |
(0.257) |
(0.206) |
(1.992) |
|
(0.268) |
(0.262) |
(0.246) |
(1.853) |
| ghg_intensity |
|
|
-4.905 |
-3.030 |
8.868 |
|
|
-3.797 |
-1.836 |
5.395 |
|
|
-4.239 |
-2.134 |
9.656 |
|
|
-4.795 |
-3.605 |
0.559 |
|
|
|
(3.883) |
(3.501) |
(6.167) |
|
|
(3.765) |
(3.226) |
(7.160) |
|
|
(3.747) |
(3.316) |
(6.199) |
|
|
(3.209) |
(3.187) |
(7.668) |
| oil_price |
|
|
|
0.861* |
|
|
|
|
0.921* |
|
|
|
|
0.906* |
|
|
|
|
0.802* |
|
|
|
|
|
(0.341) |
|
|
|
|
(0.351) |
|
|
|
|
(0.346) |
|
|
|
|
(0.331) |
|
| E_score |
|
|
|
|
|
0.468 |
0.475 |
0.616 |
0.114 |
-2.893+ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.531) |
(0.556) |
(0.561) |
(0.498) |
(1.492) |
|
|
|
|
|
|
|
|
|
|
| imr_E |
|
|
|
|
|
-23.064 |
-20.742 |
-12.045 |
-8.428 |
41.500 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(17.648) |
(17.893) |
(18.072) |
(14.651) |
(50.972) |
|
|
|
|
|
|
|
|
|
|
| S_score |
|
|
|
|
|
|
|
|
|
|
0.460 |
0.476 |
0.631 |
0.161 |
-2.187 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.507) |
(0.533) |
(0.554) |
(0.506) |
(1.224) |
|
|
|
|
|
| imr_S |
|
|
|
|
|
|
|
|
|
|
-26.949 |
-24.520 |
-16.082 |
-8.664 |
59.453 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(19.387) |
(19.210) |
(19.108) |
(15.528) |
(40.541) |
|
|
|
|
|
| G_score |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.929 |
0.925 |
1.087 |
0.674 |
-0.839 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.651) |
(0.672) |
(0.702) |
(0.681) |
(1.191) |
| imr_G |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
-24.666 |
-22.647 |
-13.613 |
-4.910 |
115.328 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(16.800) |
(17.282) |
(17.040) |
(15.998) |
(66.900) |
| Num.Obs. |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
| R2 |
0.078 |
0.086 |
0.098 |
0.182 |
0.593 |
0.065 |
0.073 |
0.081 |
0.176 |
0.614 |
0.067 |
0.076 |
0.085 |
0.177 |
0.594 |
0.104 |
0.110 |
0.123 |
0.198 |
0.562 |
| Std.Errors |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
| FE: gvkey |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
| FE: year |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
# Run regression without Heckman selection
reg_df3 <- df2 |>
filter(!is.na(ESG_composite_score),
!is.na(ghg_intensity),
!is.na(gas_percent),
!is.na(oil_price))
run_esg_set_fe <- function(data, x_var, label) {
reg_df3 <- data |>
filter(!is.na(ROE),
!is.na(.data[[x_var]]))
mods3 <- list(
feols(
as.formula(paste0("ROE ~ ", x_var)),
data = reg_df3,
vcov = ~gvkey
),
feols(
as.formula(paste0("ROE ~ ", x_var, " + gas_percent")),
data = reg_df3,
vcov = ~gvkey
),
feols(
as.formula(paste0("ROE ~ ", x_var, " + gas_percent + ghg_intensity")),
data = reg_df3,
vcov = ~gvkey
),
feols(
as.formula(paste0("ROE ~ ", x_var, " + gas_percent + ghg_intensity + oil_price")),
data = reg_df3,
vcov = ~gvkey
),
feols(
as.formula(paste0("ROE ~ ", x_var, " + ghg_intensity + gas_percent | gvkey + year")),
data = reg_df3,
vcov = ~gvkey
)
)
names(mods3) <- c(
paste0(label, " (1)"),
paste0(label, " (2)"),
paste0(label, " (3)"),
paste0(label, " (4)"),
paste0(label, " Firm+Year FE (5)")
)
mods3
}
# Run model sets
mods3_ESG <- run_esg_set_fe(reg_df3, "ESG_composite_score", "Composite ESG")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods3_E <- run_esg_set_fe(reg_df3, "E_score", "E Score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods3_S <- run_esg_set_fe(reg_df3, "S_score", "S Score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
mods3_G <- run_esg_set_fe(reg_df3, "G_score", "G Score")
## NOTE: 1/6 fixed-effect singletons were removed (7 observations).
output_table3 <- modelsummary(
c(mods3_ESG, mods3_E, mods3_S, mods3_G),
stars = TRUE,
statistic = "({std.error})",
gof_omit = "AIC|BIC|Log.Lik|RMSE|Adj|Within|Pseudo"
)
output_table3
| |
Composite ESG (1) |
Composite ESG (2) |
Composite ESG (3) |
Composite ESG (4) |
Composite ESG Firm+Year FE (5) |
E Score (1) |
E Score (2) |
E Score (3) |
E Score (4) |
E Score Firm+Year FE (5) |
S Score (1) |
S Score (2) |
S Score (3) |
S Score (4) |
S Score Firm+Year FE (5) |
G Score (1) |
G Score (2) |
G Score (3) |
G Score (4) |
G Score Firm+Year FE (5) |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
| (Intercept) |
-26.935 |
-34.134 |
0.968 |
-64.177 |
|
-17.174 |
-24.208 |
7.487 |
-65.457+ |
|
-15.647 |
-23.877 |
12.694 |
-63.850 |
|
-42.237 |
-48.764 |
-14.018 |
-70.333 |
|
|
(24.705) |
(26.743) |
(23.322) |
(38.201) |
|
(17.642) |
(18.909) |
(21.947) |
(36.695) |
|
(17.381) |
(20.042) |
(18.631) |
(35.971) |
|
(32.985) |
(34.030) |
(27.254) |
(44.445) |
|
| ESG_composite_score |
0.905 |
0.891 |
1.055 |
0.468 |
-3.294+ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.721) |
(0.749) |
(0.793) |
(0.736) |
(1.530) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| gas_percent |
|
0.253 |
0.261 |
0.161 |
2.193 |
|
0.247 |
0.253 |
0.155 |
3.171 |
|
0.273 |
0.285 |
0.164 |
2.441 |
|
0.242 |
0.249 |
0.159 |
2.183 |
|
|
(0.266) |
(0.249) |
(0.211) |
(2.248) |
|
(0.234) |
(0.215) |
(0.185) |
(2.667) |
|
(0.262) |
(0.253) |
(0.190) |
(2.298) |
|
(0.264) |
(0.240) |
(0.227) |
(2.540) |
| ghg_intensity |
|
|
-6.311+ |
-3.649 |
12.289* |
|
|
-5.350 |
-2.899 |
7.004 |
|
|
-6.418+ |
-3.235 |
12.915* |
|
|
-6.501+ |
-4.174 |
4.319 |
|
|
|
(3.515) |
(2.947) |
(4.583) |
|
|
(3.272) |
(2.507) |
(4.890) |
|
|
(3.564) |
(2.888) |
(4.967) |
|
|
(3.024) |
(2.900) |
(6.802) |
| oil_price |
|
|
|
0.870* |
|
|
|
|
0.929* |
|
|
|
|
0.927* |
|
|
|
|
0.819* |
|
|
|
|
|
(0.351) |
|
|
|
|
(0.360) |
|
|
|
|
(0.365) |
|
|
|
|
(0.345) |
|
| E_score |
|
|
|
|
|
0.692 |
0.673 |
0.755 |
0.206 |
-3.171* |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.582) |
(0.598) |
(0.619) |
(0.551) |
(1.440) |
|
|
|
|
|
|
|
|
|
|
| S_score |
|
|
|
|
|
|
|
|
|
|
0.596 |
0.600 |
0.771 |
0.223 |
-2.556+ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.554) |
(0.584) |
(0.625) |
(0.567) |
(1.356) |
|
|
|
|
|
| G_score |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.060 |
1.040 |
1.198 |
0.704 |
-1.080 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(0.749) |
(0.755) |
(0.794) |
(0.729) |
(1.612) |
| Num.Obs. |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
55 |
55 |
55 |
55 |
48 |
| R2 |
0.052 |
0.066 |
0.094 |
0.182 |
0.585 |
0.043 |
0.055 |
0.077 |
0.174 |
0.612 |
0.032 |
0.047 |
0.076 |
0.174 |
0.588 |
0.074 |
0.086 |
0.116 |
0.197 |
0.535 |
| Std.Errors |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
by: gvkey |
| FE: gvkey |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
| FE: year |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
|
|
|
|
X |
#trend aligns with that of the bigger sample using Heckman selection model where there is a significant negative relationship between ROA and E score only
#Check missing data - number of companies
n_distinct(esg_score$gvkey)
## [1] 84
n_distinct(financial$gvkey)
## [1] 210
n_distinct(gas_production$gvkey)
## [1] 223
n_distinct(ghg_intensity$gvkey)
## [1] 306
n_distinct(joint_data$gvkey)
## [1] 210
n_distinct(df2$gvkey)
## [1] 199
n_distinct(reg_df3$gvkey)
## [1] 14
#Check whether Heckman selection is significant - ie. whether selection bias matters
summary(mods2_E[[4]])
## OLS estimation, Dep. Var.: ROE
## Observations: 55
## Standard-errors: Clustered (gvkey)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.196160 29.374147 -1.913116 0.079889 .
## E_score 0.113816 0.497530 0.228763 0.822905
## gas_percent 0.133868 0.200660 0.667138 0.517307
## ghg_intensity -1.835545 3.226436 -0.568908 0.579911
## oil_price 0.920841 0.351357 2.620813 0.022351 *
## imr_E -8.428039 14.650746 -0.575263 0.575741
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 44.7 Adj. R2: 0.091926
summary(mods2_ESG[[4]])
## OLS estimation, Dep. Var.: ROE
## Observations: 55
## Standard-errors: Clustered (gvkey)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -58.366346 32.110058 -1.817697 0.094147 .
## ESG_composite_score 0.415333 0.695196 0.597433 0.561323
## gas_percent 0.149381 0.226752 0.658783 0.522475
## ghg_intensity -3.029769 3.501215 -0.865348 0.403816
## oil_price 0.861109 0.340599 2.528222 0.026506 *
## imr_ESG -4.787412 14.443990 -0.331447 0.746024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 44.5 Adj. R2: 0.09898
#not significant
# Multicollinearity checks
cor(df2[, c("ghg_intensity", "gas_percent", "oil_price")], use = "complete.obs")
## ghg_intensity gas_percent oil_price
## ghg_intensity 1.0000000 -0.15633846 -0.23240825
## gas_percent -0.1563385 1.00000000 0.06177314
## oil_price -0.2324082 0.06177314 1.00000000
# Check VIF of models without Heckman selection
vif(mods3_ESG[[4]])
## GVIF Df GVIF^(1/(2*Df))
## ESG_composite_score 1.823001 0 Inf
## gas_percent 1.823001 0 Inf
## ghg_intensity 1.823001 0 Inf
## oil_price 1.823001 0 Inf
vif(mods3_E[[4]])
## GVIF Df GVIF^(1/(2*Df))
## E_score 1.506291 0 Inf
## gas_percent 1.506291 0 Inf
## ghg_intensity 1.506291 0 Inf
## oil_price 1.506291 0 Inf
vif(mods3_S[[4]])
## GVIF Df GVIF^(1/(2*Df))
## S_score 1.981436 0 Inf
## gas_percent 1.981436 0 Inf
## ghg_intensity 1.981436 0 Inf
## oil_price 1.981436 0 Inf
vif(mods3_G[[4]])
## GVIF Df GVIF^(1/(2*Df))
## G_score 2.000362 0 Inf
## gas_percent 2.000362 0 Inf
## ghg_intensity 2.000362 0 Inf
## oil_price 2.000362 0 Inf
#GVIG values are all within 1-5 range -> acceptable range -> no indication of multicollinearity