Introduction

A concern that automation will lead to a destruction of jobs was present ever since the first industrial revolution. Events like the Great Recession further facilitate the uncertainty and fear of job loss (Messenger (2018)). One particular sector which is heavily affected by the technological change is the energy sector. Inventions like solar panels, wind turbines, and geothermal power plants lead to the wide usage of new means of energy production. With the climate change concerns, the shift toward low-or-no-emission renewable energy has been heavily supported by various regulations, the most notable being the Kyoto protocol. At the same time, a dread that closing carbon mines will lead to massive structural unemployment has been a topic of the public debate, perhaps since Margaret Thatcher rapidly decided to stop the inefficient mining of coal in the UK in the 70. of the XX century. The shift towards green energy is also hindered in places, where labor unions are strong and which have long traditions of coal mining. For example, the consecutive leaders of Poland have been unable to close Polish coal mines, despite the mining being gainless for the last 20 years or so.

Despite the wide fear that the increasing share of green energy will lead to massive unemployment, the existing research shows a positive effect of the shift towards renewable on employment (Costantini, Crespi, and Paglialunga (2018), Fankhauser, Sehlleier, and Stern (2008), Llera et al. (2013), Ortega et al. (2020), Ryan, Convery, and Ferreira (2006)). Through the productivity effect, the introduction of capital-intensive alternative sources of energy production, facilitates economic growth and, thus, leads to creation of new jobs, rather than destruction of the old ones. However, the composition of jobs changes as a result of this process. On one hand, the increased production of green energy leads to production of jobs for high skilled technicians and engineers. On the other hand, it leads to creation of service jobs for the middle educated. Through the income effect, a possible employment loss or gains can also have a significant negative or positive influence on fertility (Guinnane (2011), Shenk et al. (2013)). Thus, a shift towards green energy can also indirectly influence individual childbearing decisions. This matter is also investigated in this assignment.

I investigate the influence of the shift towards renewable energy on employment and fertility in European Union - 26 countries. I build fixed effects panel models, where I regress the total, gender-specific, and education specific employment rates, as well as the total fertility rate on the share of fossil fuels in total energy production and I control for the gdp growth. To check the robustness of my findings, I also compare bootstraped regression results with my initial results.

Methods

The table below presents the variable names and description.

a <- c("et", "ef", "em", "e02", "e34", "e58", "fossil", "gdp", "tfr")
b <- c("total employment rate", "female-specific employment rate", "male-specific employment rate", "low-educated employemnt rate (isced 0-2)", "middle-educated employment rate (isced 3-4)", "highly-educated employment rate (isced 5-8)", "share of fossil fuels in total energy production", "gdp per capita (real prices) growth", "total fertility rate")
c <- c("%", "%", "%", "%", "%", "%", "%", "%", "births per woman")
d <- c("Eurostat","Eurostat","Eurostat","Eurostat","Eurostat","Eurostat","Eurostat","Eurostat","Eurostat")
df <- data.frame(a,b,c,d)
colnames(df) <- c("variable name", "variable description", "unit", "source")
df %>%
  kbl() %>%
  kable_styling()
variable name variable description unit source
et total employment rate % Eurostat
ef female-specific employment rate % Eurostat
em male-specific employment rate % Eurostat
e02 low-educated employemnt rate (isced 0-2) % Eurostat
e34 middle-educated employment rate (isced 3-4) % Eurostat
e58 highly-educated employment rate (isced 5-8) % Eurostat
fossil share of fossil fuels in total energy production % Eurostat
gdp gdp per capita (real prices) growth % Eurostat
tfr total fertility rate births per woman Eurostat

There are 7 different dependent variables, one main regressor, and one control variable. The most important dependent variables are the total employment rate and the total fertility rate. However, I also check the influence of my main covariate on sex-specific and education-specific employment rates, because this can shed more light on the issue, which society groups are most affected by the change in energy production composition. This, in turn, can deliver presumption about, through what channels, fertility can be affected. While it is not entirely clear, what effects the shift towards green energy has on employment in Europe, there is a broad consensus that economic, incl. labour market related, factors affect fertility by affecting family income and partners’ time availability (Guinnane (2011), Shenk et al. (2013)). The direction of this effect is unclear, however. If more jobs are destroyed than produced and men are more exposed to job loss than women, green energy production will likely lead to a decline in fertility. In contrast, if job creation outruns destruction, fertility should theoretically increase, although this also depends on whether the displaced workers are able to find a new job, e.g. in the service sector, and whether the newly created jobs are compatible with childbearing and childcaring.

My main regressor of interest is the share of energy produced by burning fossil fuels. The logic behind using this variable to assess the influence of the transition to green energy on employment and fertility is as follows: if more green energy plants (solar farms, wind turbines, etc.) are introduced, the composite of the energy production changes in favor of the green energy, i.e. the share of fossils declines. Following the literature, I expect that a decline in the share of fossils boosts employment. This will be verified empirically in the next section.

In addition, I control for the GDP per capita (real prices) growth. This is in order to account for the job creation/destruction which might have been due to overall economic growth or decline. I do not include more controls, because, due to the macroeconomic framework that I have, I am concerned that including more covariates such as import competition, for example, would lead rather to multicollinearity issues, than to greater model interpretability. All in all, my estimation equations can be expressed as follows:

\[ et_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (1)\\ ef_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (2)\\ em_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (3)\\ e02_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (4)\\ e34_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (5)\\ e58_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (6)\\ tfr_{it} = \beta_0 + \beta_1 fossil_{i, t-1} + \beta_2 gdp_{i,t-1} +\beta_3 country_i + \beta_4 year_t + \eta_i + \varepsilon_{it} \,\,\,\,\,\,\,\,\, (7) \] where \(i\)-country, \(t\)-year, \(country_i\) is a country fixed effect, \(year_t\) is a time fixed effect, \(\eta_i\) is the individual effect, and \(\varepsilon_{it}\) is a random disturbance.

In order the check the robustness of the coefficients estimated with the within FE panel model, in the second step of the analysis, I also run a bootstrap procedure. Comparing real and bootstrapped coefficients will let me check the accuracy of the coefficients estimated with the FE method.

Dataset Description

The dataset comprises of observations for 26 EU countries for years 2000-2019. Below I load and browse the dataset. I prepared the dataset earlier in a separate script, where I cleaned and merged the raw files from the Eurostat website.

load("fossil_empl_rate_tfr.Rdata")
df <- subset(data_countries, year>1999)
df <- df[,c(1:10,12)]
# keeping only the EU countries
unique <- c("AT","BE","BG", "CY", "CZ", "DE", "DK", "EE", "ES", "FI", "FR", "EL", "HR",
            "HU", "IE", "IT", "LT", "LU", "LV", "NL", "PO", "PT", "RO", "SE", "SI", "SK")
df <- df[df$country %in% unique, ]

# missing data
na <- is.na(df)
sum(na)/sum(1-na)
## [1] 0.06630477
# dataset summary
kable(summarize(df)) %>% kable_styling
N Missing Mean SD Min Q1 Median Q3 Max
ef 497 3 59.27 7.83 39.60 53.30 59.30 65.20 75.90
em 497 3 70.48 5.97 52.70 66.50 70.40 74.80 82.80
et 497 3 64.85 6.25 48.80 60.10 64.80 69.40 78.20
e02 462 38 40.97 12.08 13.30 32.60 42.00 49.20 67.90
e34 462 38 68.64 6.66 46.30 64.60 67.70 73.70 83.00
e58 462 38 83.15 3.79 67.60 81.10 83.80 85.80 90.80
fossil 500 0 75.93 14.62 30.29 68.12 77.89 86.41 98.61
gdp 483 17 2.22 3.73 -14.30 0.50 2.10 4.10 24.00
tfr 298 202 1.57 0.19 1.21 1.42 1.54 1.70 2.06

As can be seen in the calculation above, about 10% of the observations is missing. This concerns mostly alternative specification of the dependent variable - the education specific employment rates. The series for total fertility rate is also shorter as it begins in 2008. Thus, for the models with this variable, a shorter series will be used. It can be also noted that, in general the level of the male-specific employment is higher than the female-specific employment. The same can be said for the employment of the highly educated in contrast to the employment of the low and middle educated. The mean of TFR for 2000-2019 in the UE is very low - only 0.19 child per 1 woman.

Below, I investigate the trend concerning each of the main variables of interest.

et
# total employment rate
plotmeans(et ~ year, main="Mean total employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

em
# male employment rate
plotmeans(em ~ year, main="Mean male employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

ef
# female employment rate
plotmeans(ef ~ year, main="Mean female employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

e02
# low educated (isco 02)
plotmeans(e02 ~ year, main="Mean low educated employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

e34
# middle educated (isco 34)
plotmeans(e34 ~ year, main="Mean middle educated employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

e58
# highly educated (isco 58)
plotmeans(e58 ~ year, main="Mean high educated employment rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

tfr
# total fertility rate
plotmeans(tfr ~ year, main="Mean total fertility rate across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

fossil
# fossil fuels in total energy production
plotmeans(fossil ~ year, main="Mean fossils share in energy production across years", data=df, ylab="%", n.label=F, col="#440154FF", barcol="#FDE725FF")

It can be seen that the total, as well as female- and male-specific employment rates exhibit an overall growing trend, with a slight decline from 2008-2012, years corresponding to the world crisis and slow recovery. The fluctuations of the education-specific employment in time are greater - it can be seen that all three groups were also subject to the world financial crisis, which increased unemployment, but there is no clear trend there. On the next graph, it can be also seen that the mean total fertility rate is declining in Europe, which is not surprising, but which is a huge problem. The share of fossil fuels in energy production is also declining, which is a result of many regulations undertaken by the EU in order to bound the greenhouse gas emissions and curb negative consequences of the climate change.

Finally, below I check the correlation between the variables in the dataset.

df1 <- df[,c(3:11)]
df1 <- df1[complete.cases(df1), ]
data_corr <- cor(df1)

corrplot.mixed(
  corr = data_corr,
  upper = 'color',
  lower = 'number',
  upper.col = viridis(n=20),
  lower.col = 'black',
  number.cex = 0.9,
  tl.cex=0.9, 
  tl.pos='d',
  font=3,
  tl.col = "black"
)

It can be seen that the correlations between different dependent variables and the covariates (fossil, gdp, imports, air) are rather mild - up to around 30%. It provides reasoning that the estimation results won’t be disturbed by multicollinearity issues.

Empirical Analysis

For the analysis, I use a very recent R package Berge and Krantz (2021), which was deployed on CRAN in March 2021. It is said to be incredibly fast and versatile. I incorporate both country and year fixed effects. To account for possible cross-sectional dependence and serial autocorrelation issues, the standard errors are clustered both by country and by year. Below, I run the analyses with respect to different dependent variables. Mind, that, for the reasons of brevity, I conduct the full analysis only with respect to the two most important dependent variables: total employment rate and total fertility rate. In other cases, the analysis is shortened and without very thorough description.

et

First, I run the fixed effects estimation of the influence of fossil fuels of the total employment rate. I start with the regular model without bootstrap.

# define panel data frame
pdf = panel(df, ~ country + year)
# run the first FE model
et1 <- feols(et ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(et1)
## OLS estimation, Dep. Var.: et
## Observations: 466 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error t value Pr(>|t|)    
## l(fossil, 1) 0.214624   0.082002  2.6173 0.017452 *  
## l(gdp, 1)    0.154208   0.103740  1.4865 0.154455    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.3705     Adj. R2: 0.837683
##                Within R2: 0.079975

It can be seen in the table above that both the share of fossils and gdp growth have a positive influence on the total employment rate. However, only the share of fossils is significant at 5% alpha level (as indicated by one asterisk). The overall \(R^2\) is relatively high - about 85%, which suggest that the model explains a great proportion of variability in data. In order to account check the robustness of the coefficients, I also run a bootstrapped replication of the betas.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(et ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##    0.2146240    0.1542082
# run the bootstrap procedure
et2 <- boot(pdf, model_coef, 500)
# see the findings
et2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.2146240 -0.000699348   0.0732076
## t2* 0.1542082  0.007044267   0.1095092

It can be seen that the bias of the original coefficients is considerably low - roughly 2%. To be even more sure, that my coefficients are relatively unbiased, I can also plot the real and bootstrapped coefficients and compare the confidence intervals which I do below.

# OLS
ci.et1<-as.data.frame(confint(et1, level=0.95, type="norm"))
colnames(ci.et1)<-c("lower", "upper")
ci.et1$variables<-rownames(ci.et1)
# Bootstrapped
ci.et2<-as.data.frame(confint(et2, level=0.95, type="norm"))
colnames(ci.et2)<-c("lower", "upper")
rownames(ci.et2) <- rownames(ci.et1)
ci.et2$variables<-rownames(ci.et2)

ggplot() + 
  geom_errorbar(data=ci.et1, mapping=aes(x=variables, ymin=lower, ymax=upper), col="#FDE725FF", size=3) +
  geom_errorbar(data=ci.et2, mapping=aes(x=variables, ymin=lower, ymax=upper, alpha=0.5), col="#440154FF", size=3) +
  labs(title="comparison of simple and bootstrapped standard errors - et") + 
  theme(plot.title = element_text(color="black", size=14, face="bold"))

As expected, it can be seen that the coefficients are practically identical. This provides reasoning that the regression coefficients are robust to a great extent.

ef
# run the first FE model
ef1 <- feols(ef ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(ef1)
## OLS estimation, Dep. Var.: ef
## Observations: 466 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error t value Pr(>|t|)    
## l(fossil, 1) 0.247327   0.076553  3.2308 0.004638 ** 
## l(gdp, 1)    0.079701   0.073255  1.0880 0.290949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.8811     Adj. R2: 0.934592
##                Within R2: 0.128674

From the above presented model it appears that the shift towards green energy has a negative effect on the employment of women and it is stronger than the effect on total employment rate (significant at 1% level vs. at 5% level).

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(ef ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##   0.24732708   0.07970094
# run the bootstrap procedure
ef2 <- boot(pdf, model_coef, 500)
# see the findings
ef2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##       original      bias    std. error
## t1* 0.24732708 0.004798282  0.07347490
## t2* 0.07970094 0.009774875  0.07292461

In this case, the bias of the coefficients is even lower - 0.2%-0.3% of the original coefficient. This suggest that the findings from the regression of the female-specific employment are even more robust than in the case of the general employment rate.

em
# run the first FE model
em1 <- feols(em ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(em1)
## OLS estimation, Dep. Var.: em
## Observations: 466 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error t value Pr(>|t|)    
## l(fossil, 1) 0.180159   0.099814  1.8050 0.087839 .  
## l(gdp, 1)    0.238551   0.138700  1.7199 0.102596    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 3.1634     Adj. R2: 0.681308
##                Within R2: 0.052403

From the above table, I conclude that the share of fossils in energy production has only a weak positive influence on male-specific employment, as the lagged fossil variable is only significant at 10% alpha level.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(em ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##    0.1801593    0.2385509
# run the bootstrap procedure
em2 <- boot(pdf, model_coef, 500)
# see the findings
em2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.1801593 -0.003114445  0.09971055
## t2* 0.2385509  0.006427880  0.13983001

The bias is in this case higher than in the case of women employment, but still considerably low - roughly 0.5% of the real coefficients. Here, the bias of the coefficients is much higher than in other cases - about 2% of the real coefficient. It is still quite a low bias but it suggests that the results of the male-specific employment estimation are less reliable than in the previous cases.

e02
# run the first FE model
e02_1 <- feols(e02 ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(e02_1)
## OLS estimation, Dep. Var.: e02
## Observations: 439 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error  t value Pr(>|t|) 
## l(fossil, 1) 0.035206   0.136508 0.257902 0.799408 
## l(gdp, 1)    0.069520   0.148107 0.469392 0.644425 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 3.4117     Adj. R2: 0.910583
##                Within R2: 0.003032

From the table above, it can be seen that the share of fossil fuels in energy production has no influence on the employment of the low-educated.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(e02 ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##   0.03520574   0.06952045
# run the bootstrap procedure
e02_2 <- boot(pdf, model_coef, 500)
# see the findings
e02_2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##       original     bias    std. error
## t1* 0.03520574 0.01210146   0.1278707
## t2* 0.06952045 0.02304274   0.1544396

Here the bias of the coefficients is considerably low - only up to 0.2% of the real coefficients. Thus, I believe that the results are generally reliable.

e34
# run the first FE model
e34_1 <- feols(e34 ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(e34_1)
## OLS estimation, Dep. Var.: e34
## Observations: 439 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error t value Pr(>|t|)    
## l(fossil, 1) 0.086504   0.076985  1.1237 0.275924    
## l(gdp, 1)    0.231939   0.126988  1.8265 0.084414 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.5462     Adj. R2: 0.835209
##                Within R2: 0.050706

From the table above, it can be seen that the share of fossil fuels in energy production has no influence on the employment of the middle-educated. However, gdp growth is significant at 10% level in this case.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(e34 ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##   0.08650417   0.23193867
# run the bootstrap procedure
e34_2 <- boot(pdf, model_coef, 500)
# see the findings
e34_2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##       original       bias    std. error
## t1* 0.08650417  0.011586535  0.08424948
## t2* 0.23193867 -0.002759513  0.12829530

Here the bias of the coefficients is considerably low - only up to 0.006% of the real coefficients. Thus, I believe that the results are generally reliable.

e58
# run the first FE model
e58_1 <- feols(e58 ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(e58_1)
## OLS estimation, Dep. Var.: e58
## Observations: 439 
## Fixed-effects: country: 25,  year: 19
## Standard-errors: Two-way (country & year) 
##              Estimate Std. Error t value Pr(>|t|)    
## l(fossil, 1) 0.211056   0.063282  3.3352 0.003684 ** 
## l(gdp, 1)    0.146158   0.082353  1.7748 0.092847 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.9579     Adj. R2: 0.701862
##                Within R2: 0.108587

From the table above, I can conclude that the shift from fossil fuels to renewable energy has a strong (1% alpha) positive effect on the highly-educated employment.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(e58 ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##    0.2110558    0.1461583
# run the bootstrap procedure
e58_2 <- boot(pdf, model_coef, 500)
# see the findings
e58_2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.2110558 -0.001366549  0.06160702
## t2* 0.1461583  0.007733708  0.09093758

Here, again the bias is very low - only 0.4% of the original coefficients.

tfr
# run the first FE model
tfr1 <- feols(tfr ~ l(fossil) + l(gdp) | country + year, data=pdf, cluster = c("country", "year"))
# table with model summary
summary(tfr1)
## OLS estimation, Dep. Var.: tfr
## Observations: 295 
## Fixed-effects: country: 25,  year: 12
## Standard-errors: Two-way (country & year) 
##               Estimate Std. Error   t value Pr(>|t|)    
## l(fossil, 1)  0.009844   0.004225  2.330100 0.039863 *  
## l(gdp, 1)    -0.000248   0.003766 -0.065881 0.948655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.07914     Adj. R2: 0.807241
##                 Within R2: 0.067662

From the above estimation, I get that, through the labor market, the share of fossils in energy production has a slightly positive significant (5%) influence on the total fertility rate. An increase in the share of fossils by 1% leads to an increase of TFR by approximately 0.01 child per woman.

# define custom sampling function
model_coef <- function(data, idx){
  coef(feols(tfr ~ l(fossil) + l(gdp) | country + year, data=pdf[idx,], cluster = c("country", "year")))
}
# run the funtion on the panel data frame
model_coef(pdf, 1:500)
## l(fossil, 1)    l(gdp, 1) 
##  0.009844292 -0.000248091
# run the bootstrap procedure
tfr2 <- boot(pdf, model_coef, 500)
# see the findings
tfr2
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pdf, statistic = model_coef, R = 500)
## 
## 
## Bootstrap Statistics :
##         original        bias    std. error
## t1*  0.009844292 -0.0002912428 0.005353538
## t2* -0.000248091  0.0004211195 0.004714498

The bias of the coefficients is, gain, extremely low. Since the TFR is my second, most important dependent variable of interest, it is worth to plot the error bars for the model and bootstrapped coefficients to be 100% sure that the results are reliable.

# OLS
ci.tfr1<-as.data.frame(confint(tfr1, level=0.95, type="norm"))
colnames(ci.tfr1)<-c("lower", "upper")
ci.tfr1$variables<-rownames(ci.tfr1)
# Bootstrapped
ci.tfr2<-as.data.frame(confint(tfr2, level=0.95, type="norm"))
colnames(ci.tfr2)<-c("lower", "upper")
rownames(ci.tfr2) <- rownames(ci.tfr1)
ci.tfr2$variables<-rownames(ci.tfr2)

ggplot() + 
  geom_errorbar(data=ci.tfr1, mapping=aes(x=variables, ymin=lower, ymax=upper), col="#FDE725FF", size=3) +
  geom_errorbar(data=ci.tfr2, mapping=aes(x=variables, ymin=lower, ymax=upper, alpha=0.5), col="#440154FF", size=3) +
  labs(title="comparison of simple and bootstrapped standard errors - tfr1") + 
  theme(plot.title = element_text(color="black", size=14, face="bold"))

The bias of the fossil variable coefficient is about 0.5%. From the plot, it can be seen that there is a slight shift in the real and bootstrapped confidence intervals for the fossil coefficient. Since the bootstrapped CI for the fossil variable is shifted a bit upwards relative to the CI for the real coefficient. This suggests that the real effect of this variable on TFR might be slightly greater than estimated by my FE model.

Conclusion

Finally, I can compare the models with respect to the different dependent variables and draw conclusions about the relation between the shift from fossil fuels to renewable energy and employment, as well as fertility.

models <- list(
  "et" = et1,
  "ef" = ef1,
  "em" = em1,
  "e02" = e02_1,
  "e34" = e34_1,
  "e58" = e58_1,
  "tfr" = tfr1)
modelsummary(models, estimate = "{estimate}{stars}")
et ef em e02 e34 e58 tfr
l(fossil, 1) 0.215** 0.247*** 0.180* 0.035 0.087 0.211*** 0.010**
(0.082) (0.077) (0.100) (0.137) (0.077) (0.063) (0.004)
l(gdp, 1) 0.154 0.080 0.239 0.070 0.232* 0.146* 0.000
(0.104) (0.073) (0.139) (0.148) (0.127) (0.082) (0.004)
Num.Obs. 466 466 466 439 439 439 295
R2 0.853 0.941 0.711 0.920 0.852 0.732 0.831
R2 Adj. 0.838 0.935 0.681 0.911 0.835 0.702 0.807
R2 Within 0.080 0.129 0.052 0.003 0.051 0.109 0.068
R2 Pseudo
AIC 2216.8 2001.3 2485.8 2413.3 2156.4 1925.7 -583.4
BIC 2403.3 2187.8 2672.3 2597.1 2340.2 2109.5 -443.3
Log.Lik. -1063.424 -955.664 -1197.897 -1161.661 -1033.204 -917.870 329.693
Std. Errors Two-way (country & year) Two-way (country & year) Two-way (country & year) Two-way (country & year) Two-way (country & year) Two-way (country & year) Two-way (country & year)
FE: country X X X X X X X
FE: year X X X X X X X

In the table above, it can be seen that the share of fossils has a significant positive impact on the overall employment, the employment of both women and men, and the employment of the highly-educated in Europe. It also has a positive significant influence on the total fertility rate. These findings are contrary to the previous research about the influence of shift towards green energy on employment. My findings suggest that withdrawal from the traditional means of energy production destroys more jobs than it creates, and thus leads to employment decline. It is not surprising, that if the effect on employment is negative, then the effect of fertility is also negative - through the income effect. Having a stable employment and salary are significant contributors to reducing economic uncertainty, and they can influence individual long-term decisions, such as a decision to have a child. In addition, I checked the reliability of my findings with a bootstrap coefficient estimation. In each case, I got a result that the FE model coefficients are not very biased (2% in the worst case). That is why I believe that the presented method leads to valid statistical inference about the causality between the energy production change, and the employment, as well as fertility in Europe.

Below I plot the coefficients with confidence intervals from the original FE regression.

plot_summs(et1, ef1, em1, e02_1, e34_1, e58_1, tfr1, 
           scale = TRUE, 
           robust = list("HC3", "HC3", "HC3", "HC3", "HC3", "HC3", "HC3"),
           model.names = c("et", "ef", "em", "e02", "e34", "e58", "tfr"),
           inner_ci_level = .9)

References

Berge, L., and S. Krantz. 2021. Fixest. https://cran.r-project.org/web/packages/fixest/fixest.pdf.
Costantini, V., F. Crespi, and E. Paglialunga. 2018. “The Employment Impact of Private and Public Actions for Energy Efficiency: Evidence from European Industries.” Energy Policy 119 (C): 250–67.
Fankhauser, S., F. Sehlleier, and N. Stern. 2008. “Climate Change, Innovation and Jobs.” CLIMATE POLICY Earthscan ISSN 80513: 421–29. https://doi.org/10.3763/cpol.2008.0513.
Guinnane, T. W. 2011. “The Historical Fertility Transition: A Guide for Economists.” Journal of Economic Literature 49 (3): 589–614. https://doi.org/10.1257/jel.49.3.589.
Llera, E., S. Scarpellini, A. Aranda, and I. Zabalza. 2013. “Forecasting job creation from renewable energy deployment through a value-chain approach.” Renewable and Sustainable Energy Reviews 21 (C): 262–71. https://doi.org/10.1016/j.rser.2012.12.05.
Messenger, J. 2018. “Working Time and the Future of Work.” Ilo Future of Work Research Paper Series 6.
Ortega, M., P. del Río, P. Ruiz, W. Nijs, and S. Politis. 2020. “Analysing the influence of trade, technology learning and policy on the employment prospects of wind and solar energy deployment: The EU case.” Renewable and Sustainable Energy Reviews 122 (C). https://doi.org/10.1016/j.rser.2019.10965.
Ryan, L., F. Convery, and S. Ferreira. 2006. “Stimulating the Use of Biofuels in the European Union: Implications for Climate Change Policy.” Energy Policy 34 (17): 3184–94.
Shenk, M. K., M. C. Towner, H. C. Kress, and Nurul A. 2013. “A Model Comparison Approach Shows Stronger Support for Economic Models of Fertility Decline.” Proceedings of the National Academy of Sciences of the United States of America 110 (20): 8045–50. https://doi.org/10.1073/pnas.1217029110.