CKC_Rcode
Bastien PATRAS, Chiara PARTON
3/25/2021
library(readr)
library(ggplot2)
library(reshape2)
library(gapminder)
library(plotly)
library(stargazer)
library(xtable)
library(dplyr)
library(MatchIt)
library(readxl)
library(broom)
library(tidyr)
library(plm)
library(rmarkdown)
library(zoo)
library(purrr)
GDP_PER_CAPITA <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/GDP_PER_CAPITA.csv")
GDP <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/gdp-world-regions-stacked-area.csv")
MEAN_YEAR_OF_SCHOOLING <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/mean-years-of-schooling-1.csv")
CAPITA_INTENSITY <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/capital-intensity.csv")
POPULATION <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/projected-population-by-country.csv")
CO2_EMMISSION <- read_csv("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/annual-co2-emissions-per-country.csv")
CLASS <- read_excel("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/CLASS.xls", skip = 3) %>%
select(Economy, Region,`Income group`) %>%
rename(Entity = Economy)
GOV_Expenditure <- read_excel("data/GOV_10A_EXP__custom_9647011621327295154.xlsx",
sheet = "Sheet 1", col_types = c("text",
"text", "numeric", "text", "numeric",
"text", "numeric", "text", "numeric",
"text", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "text", "numeric",
"numeric", "numeric", "numeric",
"numeric", "text", "numeric", "text"),
skip = 10, n_max = 35)\(~\)
# Data prep
dta.matching <- left_join(GDP, POPULATION, by = c("Entity", "Code", "Year")) %>% rename(Population = `Population by country and region, historic and projections (Gapminder, HYDE & UN)`) %>%
na.omit()
dta.matching <- left_join(dta.matching, CO2_EMMISSION, by = c("Entity", "Code", "Year")) %>% rename(CO2 = `Annual CO2 emissions`) %>%
na.omit()
dta.matching <- left_join(dta.matching, CLASS, by = c("Entity"))
dta.matching <- left_join(dta.matching, MEAN_YEAR_OF_SCHOOLING, by = c("Entity", "Code", "Year")) %>% rename(Average.Schooling = `Average Total Years of Schooling for Adult Population (Lee-Lee (2016), Barro-Lee (2018) and UNDP (2018))`)
dta.matching <- left_join(dta.matching, CAPITA_INTENSITY, by = c("Entity", "Code", "Year")) %>%
rename(Capital.intensity = `Capital Intensity (Bergeaud, Cette, and Lecat (2016))`)
dta.matching <- dta.matching %>%
group_by(Entity) %>%
mutate(GDP.sqr = GDP^2,
GDP.cube = GDP^3,
GDP.growth = (GDP-dplyr::lag(GDP))/dplyr::lag(GDP),
Lagged.GDP.1 = dplyr::lag(GDP),
Lagged.GDP.2 = dplyr::lag(GDP, n=2))
## Data government expenditure
Year <- replicate((2019-1995), 0)
for (i in 1:(length(Year)+1)) {
Year[i] <- i+1994
}
Year.filter <- c(as.character(Year))
GOV_Expenditure <- GOV_Expenditure %>%
select(TIME, matches(Year.filter)) %>%
rename(Country = TIME) %>%
pivot_longer(cols = matches(Year.filter), names_to = "TIME", values_to = "Expenditure") %>%
na.omit() %>%
group_by(Country) %>%
mutate(Moving.average.3 = zoo::rollmean(Expenditure, k = 3, fill = NA),
Treatment.3 = ifelse(Expenditure > Moving.average.3, 1, 0),
Country = ifelse(Country == "Germany (until 1990 former territory of the FRG)","Germany", Country))
## Mearging datasets
GOV_Expenditure.1 <- GOV_Expenditure %>% rename(Entity = Country, Year = TIME)%>%
mutate(Year = as.numeric(Year))
dta.matching <- left_join(dta.matching, GOV_Expenditure.1, by = c("Entity", "Year"))\(~\)
1.1 - Stylised facts - C02 Emission
\(~\)
plot.1.1 <- dta.matching %>%
filter(Entity == "France" |
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain" |
Entity == "United Kingdom" |
Entity == "Greece") %>%
ggplot(aes(x=Year, y=CO2, color=Entity)) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray")) +
scale_color_manual(values = c("#D9ED92", "#B5E48C",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))
ggplotly(plot.1.1)\(~\)
1.2 - Stylised facts - GDP Growth
\(~\)
plot.1.2 <- dta.matching %>%
filter(Entity == "France" |
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain" |
Entity == "United Kingdom" |
Entity == "Greece") %>%
ggplot(aes(x=Year, y=GDP, color=Entity)) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray")) +
scale_color_manual(values = c("#D9ED92", "#B5E48C",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))
ggplotly(plot.1.2)\(~\)
2 - Naive estimation
\(~\)
- First we proceed according to our estimation strategy
- Then, as our functional form is third degree polynomial, we compute the inflection points according to the deltas:
\[ f(\texttt{GDP}) : \texttt{CO}_2 = \beta_0 + \beta_1\times \texttt{GDP}_{it} + \beta_2 \times\texttt{GDP}^2_{it}+ \beta_3 \times\texttt{GDP}^3_{it} + \epsilon_{it}\] Differentiating with respect to GDP gives the FOC such that:
\[ \frac{\partial f(\texttt{GDP})}{\partial \texttt{GDP}} = \beta_1 + 2\beta_2\times\texttt{GDP}_{it} + 3\beta_3\times\texttt{GDP}_{it}^2=0 \] Computing the \(\Delta\) by country \(i\) gives:
\[ \Delta_i = (2\times\beta_{2i})^2-12\times \beta_{1i}\times \beta_{3i} \] For the theory to hold we need this \(\Delta\) to be equal to 0 otherwise there is more than one inflection points
\(~\)
# Running OLS regression by countries
dta.matching.estimates.time <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube, data = .))) %>%
unnest(Country.reg)
table.naive <- dta.matching.estimates.time %>%
filter(Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain"|
Entity == "United Kingdom" |
Entity == "Greece")
# Latex output
#xtable(table.naive, digits = -3)
# Output
paged_table(table.naive)\(~\)
# Defining function for computing turning points
turn <- function(GDP.sqr.hat, GDP.cube.hat){
GDP.star <- -GDP.sqr.hat/(3*GDP.cube.hat)
return(GDP.star)
}
# Computing the inflection points by country
Inflection.table <- table.naive %>% select(Entity, term, estimate) %>%
pivot_wider(id_cols = Entity, names_from = term, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr, GDP.cube.hat = GDP.cube) %>%
mutate(Delta_i = (2*GDP.sqr.hat)^2 - 12*GDP.cube.hat*GDP.hat,
`Truning points` = ifelse(Delta_i<0.01, abs(turn(GDP.sqr.hat,GDP.cube.hat)), "NA")) %>%
select(Entity, Delta_i, `Truning points`)
# Latex output
#xtable(Inflection.table, digits = -3)
# Output
paged_table(Inflection.table)\(~\)
2 - Estimation controlling for pop density
\(~\)
# Running OLS regression by countries
dta.matching.estimates <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population, data = .))) %>%
unnest(Country.reg)
dta.matching.estimates.reduced <- dta.matching.estimates %>%
select(Entity, term, estimate) %>%
rename(Regressor = term) %>%
pivot_wider(names_from = Regressor, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr,
GDP.cube.hat = GDP.cube, Population.hat = Population)
# Merging etimates to dta.matching
dta.matching <- left_join(dta.matching, dta.matching.estimates.reduced, by = c("Entity"))
# Computing empirical Kuznets curve by country
dta.matching <- dta.matching %>%
group_by(Entity, Year) %>%
mutate(Fitted.Kuznets = GDP.hat*GDP + GDP.sqr.hat*GDP.sqr + GDP.cube.hat*GDP.cube + Population.hat*Population + `(Intercept)`)
table.pop <- dta.matching.estimates %>%
filter(Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain" |
Entity == "United Kingdom" |
Entity == "Greece")
# Latex output
#xtable(table.pop, digits = -3)
# Output
paged_table(table.pop)# Running Wald Test for France
m1.dta <- dta.matching %>% filter(Entity == "France")
# OLS Estimation
m1.pop.fr <- lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population, data = m1.dta)
m1.fr <- lm(CO2 ~ GDP + GDP.sqr + GDP.cube, data = m1.dta)
# Extracting vcov matrixes
cor.pop <- as.matrix(vcov(m1.pop.fr))
cor <- as.matrix(vcov(m1.fr))
# Formating
dta.test.pop <- tidy(m1.pop.fr) %>% select(term, estimate)
dta.test <- tidy(m1.fr) %>% select(term, estimate)
dta.Wald <- left_join(dta.test.pop, dta.test, by=c("term")) %>%
rename(With.pop = estimate.x, Without.pop = estimate.y) %>%
mutate(Without.pop = ifelse(is.na(Without.pop), 1, Without.pop))
beta.vector <- dta.Wald$Without.pop
r.vector <- dta.Wald$With.pop
R.vector <- diag(c(1, 1, 1, 1, 0))\(~\)
# Defining function for computing turning points
turn <- function(GDP.sqr.hat, GDP.cube.hat){
GDP.star <- -GDP.sqr.hat/(3*GDP.cube.hat)
return(GDP.star)
}
# Computing the inflection points by country
Inflection.table.pop <- table.pop %>% select(Entity, term, estimate) %>%
pivot_wider(id_cols = Entity, names_from = term, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr, GDP.cube.hat = GDP.cube) %>%
mutate(Delta_i = (2*GDP.sqr.hat)^2 - 12*GDP.cube.hat*GDP.hat,
`Truning points` = ifelse(Delta_i<0.01, abs(turn(GDP.sqr.hat,GDP.cube.hat)), "NA")) %>%
select(Entity, Delta_i, `Truning points`)
# Latex output
#xtable(Inflection.table.pop, digits = -3)
# Output
paged_table(Inflection.table.pop)\(~\)
2 - Estimation controlling for density accounting for time trend
\(~\)
# Running OLS regression by countries
dta.matching.estimates.time <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population + Year, data = .))) %>%
unnest(Country.reg)
table.time <- dta.matching.estimates.time %>%
filter(Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain"|
Entity == "United Kingdom" |
Entity == "Greece")
# Latex output
#xtable(table.time, digits = -3)
# Output
paged_table(table.time)\(~\)
# Computing the inflection points by country
Inflection.table.time <- table.time %>% select(Entity, term, estimate) %>%
pivot_wider(id_cols = Entity, names_from = term, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr, GDP.cube.hat = GDP.cube) %>%
mutate(Delta_i = (2*GDP.sqr.hat)^2 - 12*GDP.cube.hat*GDP.hat,
`Truning points` = ifelse(Delta_i<0.01, abs(turn(GDP.sqr.hat,GDP.cube.hat)), "NA")) %>%
select(Entity, Delta_i, `Truning points`)
# Latex output
#xtable(Inflection.table.time, digits = -3)
# Output
paged_table(Inflection.table.time)\(~\)
Raw plotting
\(~\)
dta.matching.sub <- dta.matching %>%
group_by(Entity) %>%
mutate(Fitted.Kuznets = GDP.hat*GDP + GDP.sqr.hat*GDP.sqr +
GDP.cube.hat*GDP.cube + Population.hat*Population + `(Intercept)`) %>%
filter(Entity == "France" |
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain"|
Entity == "United Kingdom" |
Entity == "Greece")
# Plotting the empirical Kuznets curves
plot.2 <- dta.matching.sub %>%
ggplot(aes(x=GDP, y=CO2, color=Entity)) +
stat_smooth(method = "lm", formula = CO2 ~ GDP + I(GDP^2) + I(GDP^3), size = 1) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray")) +
scale_color_manual(values = c("#D9ED92", "#B5E48C",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))
ggplotly(plot.2)\(~\)
# Plotting the empirical Kuznets curves
plot.1 <- dta.matching.sub %>%
ggplot(aes(x=GDP, y=Fitted.Kuznets, color=Entity)) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray")) +
scale_color_manual(values = c("#D9ED92", "#B5E48C",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))
ggplotly(plot.1)\(~\)
2 - Estimation controlling for density accounting for time trend and evironmental policies
\(~\)
# Running OLS regression by countries
dta.matching.estimates.time <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population + Year, data = .))) %>%
unnest(Country.reg)
table.time <- dta.matching.estimates.time %>%
filter(Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "Spain"|
Entity == "United Kingdom" |
Entity == "Greece")
# Latex output
#xtable(table.time, digits = -3)
# Output
paged_table(table.time)\(~\)
2 - Truning points for every countries with time trend
\(~\)
# Running OLS regression by countries
dta.matching.estimates.time <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population + Year, data = .))) %>%
unnest(Country.reg)
table.time.turning <- dta.matching.estimates.time
# Defining function for computing turning points
turn <- function(GDP.sqr.hat, GDP.cube.hat){
GDP.star <- -GDP.sqr.hat/(3*GDP.cube.hat)
return(GDP.star)
}
# Computing the inflection points by country
table.time.turning <- table.time.turning %>% select(Entity, term, estimate) %>%
pivot_wider(id_cols = Entity, names_from = term, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr, GDP.cube.hat = GDP.cube) %>%
mutate(Delta_i = (2*GDP.sqr.hat)^2 - 12*GDP.cube.hat*GDP.hat,
`Turning points` = ifelse(Delta_i<0.01, abs(turn(GDP.sqr.hat,GDP.cube.hat)), "NA")) %>%
select(Entity, Delta_i, `Turning points`)
paged_table(table.time.turning)\(~\)
3 - Plotting maps for GDP_2018 with time trend effects
\(~\)
### Computing distance
# Extracting GDP from 2018
GDP_2018 <- dta.matching %>% filter(Year == 2018) %>%
select(Entity, GDP)
table.time.turning.2018 <- left_join(table.time.turning, GDP_2018, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1)
paged_table(table.time.turning.2018)\(~\)
3 - Plotting maps for GDP_2018 with time trend effects
\(~\)
### Creating blank map
# Geospatial data available at the geojson format
library(geojsonio)
library(mapproj)
library(RColorBrewer)
library(scales)
spdf <- geojson_read("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/world/world.geojson", what = "sp")
# Fortifying data
spdf_fortified <- tidy(spdf, region = "NAME_LONG")
plot <- ggplot() +
geom_polygon(data = spdf_fortified, aes( x = long, y = lat, group = group), fill="white", color="grey") +
theme_void() +
coord_map()
# Merging fortified data with turning point table
table.time.turning.2018 <- table.time.turning.2018 %>% rename(id = Entity)
spdf_fortified.merged.2018 <- left_join(spdf_fortified, table.time.turning.2018, by = c("id")) %>%
mutate(Distance = ifelse(is.na(Distance), 0.001, Distance))
plot.map <- ggplot() +
geom_polygon(data = spdf_fortified.merged.2018, aes(fill = Distance, x = long, y = lat, group = group), size=0, alpha=0.9, col = "white") +
theme_void() +
coord_map() +
scale_fill_gradient(low = "#34A0A4",
high = "#D9ED92") +
theme(legend.position = c(0.1, 0.2))
show(plot.map)
\(~\)
3 - Plotting maps for GDP_2005 with time trend effects
\(~\)
### Computing distance
# Extracting GDP from 2018
GDP_2005 <- dta.matching %>% filter(Year == 2005) %>%
select(Entity, GDP)
table.time.turning.2005 <- left_join(table.time.turning, GDP_2005, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1)
paged_table(table.time.turning.2005)\(~\)
3 - Plotting maps for GDP_2005 with time trend effects
\(~\)
### Creating blank map
# Geospatial data available at the geojson format
library(geojsonio)
library(mapproj)
library(RColorBrewer)
library(scales)
spdf <- geojson_read("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/world/world.geojson", what = "sp")
# Fortifying data
spdf_fortified.2005 <- tidy(spdf, region = "NAME_LONG")
plot <- ggplot() +
geom_polygon(data = spdf_fortified.2005, aes( x = long, y = lat, group = group), fill="white", color="grey") +
theme_void() +
coord_map()
# Merging fortified data with turning point table
table.time.turning.2005 <- table.time.turning.2005 %>% rename(id = Entity)
spdf_fortified.merged.2005 <- left_join(spdf_fortified, table.time.turning.2005, by = c("id")) %>%
mutate(Distance = ifelse(is.na(Distance), 0.001, Distance))
plot.map.2005 <- ggplot() +
geom_polygon(data = spdf_fortified.merged.2005, aes(fill = Distance, x = long, y = lat, group = group), size=0, alpha=0.9, col = "white") +
theme_void() +
coord_map() +
scale_fill_gradient(low = "#34A0A4",
high = "#D9ED92") +
theme(legend.position = c(0.1, 0.2))
show(plot.map.2005)
\(~\)
2 - Truning points for every countries without time trend
\(~\)
# Running OLS regression by countries
dta.matching.estimates.without <- dta.matching %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube + Population, data = .))) %>%
unnest(Country.reg)
table.turning.without <- dta.matching.estimates.without
# Defining function for computing turning points
turn <- function(GDP.sqr.hat, GDP.cube.hat){
GDP.star <- -GDP.sqr.hat/(3*GDP.cube.hat)
return(GDP.star)
}
# Computing the inflection points by country
table.turning.without <- table.turning.without %>% select(Entity, term, estimate) %>%
pivot_wider(id_cols = Entity, names_from = term, values_from = estimate) %>%
rename(GDP.hat = GDP, GDP.sqr.hat = GDP.sqr, GDP.cube.hat = GDP.cube) %>%
mutate(Delta_i = (2*GDP.sqr.hat)^2 - 12*GDP.cube.hat*GDP.hat,
`Turning points` = ifelse(Delta_i<0.01, abs(turn(GDP.sqr.hat,GDP.cube.hat)), "NA")) %>%
select(Entity, Delta_i, `Turning points`)
paged_table(table.turning.without)\(~\)
3 - Plotting maps for GDP_2018 without time trend effects
\(~\)
### Computing distance
# Extracting GDP from 2018
GDP_2018 <- dta.matching %>% filter(Year == 2018) %>%
select(Entity, GDP)
table.turning.without.2018 <- left_join(table.turning.without, GDP_2018, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1)
paged_table(table.turning.without.2018)\(~\)
3 - Plotting maps for GDP_2018 without time trend effects
\(~\)
### Creating blank map
# Geospatial data available at the geojson format
library(geojsonio)
library(mapproj)
library(RColorBrewer)
library(scales)
spdf <- geojson_read("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/world/world.geojson", what = "sp")
# Fortifying data
spdf_fortified.2018.without <- tidy(spdf, region = "NAME_LONG")
plot <- ggplot() +
geom_polygon(data = spdf_fortified.2018.without, aes( x = long, y = lat, group = group), fill="white", color="grey") +
theme_void() +
coord_map()
# Merging fortified data with turning point table
table.turning.without.2018 <- table.turning.without.2018 %>% rename(id = Entity)
spdf_fortified.merged.without.2018 <- left_join(spdf_fortified, table.turning.without.2018, by = c("id")) %>%
mutate(Distance = ifelse(is.na(Distance), 0.001, Distance))
plot.map.without.2018 <- ggplot() +
geom_polygon(data = spdf_fortified.merged.without.2018, aes(fill = Distance, x = long, y = lat, group = group), size=0, alpha=0.9, col = "white") +
theme_void() +
coord_map() +
scale_fill_gradient(low = "#34A0A4",
high = "#D9ED92") +
theme(legend.position = c(0.1, 0.2))
show(plot.map.without.2018)
\(~\)
3 - Plotting maps for GDP_2005 without time trend effects
\(~\)
### Computing distance
# Extracting GDP from 2018
GDP_2005 <- dta.matching %>% filter(Year == 2005) %>%
select(Entity, GDP)
# Computing the relative distance as a percentage
table.without.turning.2005 <- left_join(table.turning.without, GDP_2005, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1)
paged_table(table.without.turning.2005)\(~\)
3 - Plotting maps for GDP_2005 without time trend effects
\(~\)
### Creating blank map
# Geospatial data available at the geojson format
library(geojsonio)
library(mapproj)
library(RColorBrewer)
library(scales)
spdf <- geojson_read("/Users/bastienpatras/Dropbox/Econometric project - Bastien, Chiara/data/world/world.geojson", what = "sp")
# Fortifying data
spdf_fortified.without.2005 <- tidy(spdf, region = "NAME_LONG")
plot <- ggplot() +
geom_polygon(data = spdf_fortified.without.2005, aes( x = long, y = lat, group = group), fill="white", color="grey") +
theme_void() +
coord_map()
# Merging fortified data with turning point table
table.without.turning.2005 <- table.without.turning.2005 %>% rename(id = Entity)
spdf_fortified.merged.without.2005 <- left_join(spdf_fortified, table.without.turning.2005, by = c("id")) %>%
mutate(Distance = ifelse(is.na(Distance), 0.001, Distance))
plot.map.without.2005 <- ggplot() +
geom_polygon(data = spdf_fortified.merged.without.2005, aes(fill = Distance, x = long, y = lat, group = group), size=0, alpha=0.9, col = "white") +
theme_void() +
coord_map() +
scale_fill_gradient(low = "#34A0A4",
high = "#D9ED92") +
theme(legend.position = c(0.1, 0.2))
show(plot.map.without.2005)
\(~\)
3.1 - Effective year of the turning point by country (without time trend) table.time.turning
\(~\)
### Computing distance by country
# Extracting GDP
GDP <- dta.matching %>% select(Entity, Year, GDP)
# Computing the year of the turning point (GDP_t < GDP* --> t*)
table.without.turning.time <- left_join(table.turning.without, GDP, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1, Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "United Kingdom" |
Entity == "Greece",
GDP<`Turning points`)
table.without.turning.time <- table.without.turning.time %>% group_by(Entity) %>%
summarise(Entity = Entity, `Turning points` = `Turning points`,
Year = max(Year)) %>%
distinct()
# Output
paged_table(table.without.turning.time)\(~\)
3.1 - Effective year of the turning point by country (with time trend)
\(~\)
### Computing distance by country
# Extracting GDP
GDP <- dta.matching %>% select(Entity, Year, GDP)
# Computing the year of the turning point (GDP_t < GDP* --> t*)
table.with.turning.time <- left_join(table.time.turning, GDP, by = c("Entity")) %>%
mutate(Distance = (GDP-`Turning points`)/`Turning points`) %>%
filter(Distance<1, Entity == "France"|
Entity == "Italy" |
Entity == "Germany" |
Entity == "United Kingdom" |
Entity == "Greece",
GDP<`Turning points`)
table.with.turning.time <- table.with.turning.time %>% group_by(Entity) %>%
summarise(Entity = Entity, `Turning points` = `Turning points`,
Year = max(Year)) %>%
distinct()
# Output
paged_table(table.with.turning.time)\(~\)
4.1 - 2SLS Estimation on European Countries
\(~\)
4.1.1 - First Stage regression
\(~\)
# Data prep
dta.matching.2SLS <- dta.matching %>%
na.omit()
# First stage
First.stage <- dta.matching.2SLS %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(GDP ~ Expenditure, data = .))) %>%
unnest(Country.reg)
# Latex output
First.stage.latex <- First.stage %>% filter(Entity == "France" |
Entity == "Germany" |
Entity == "Italy" |
Entity == "United Kingdom")
#xtable(First.stage.latex, digits=-3)
# Output
paged_table(First.stage.latex)\(~\)
4.1.2 - Second Stage regression
\(~\)
# Building variables
Fitted.First.Stage.dta <- dta.matching.2SLS %>%
select(Entity, Year, CO2, GDP, Population,
Average.Schooling, Capital.intensity)
First.Stage.dta <- First.stage %>% select(Entity, term, estimate)
Fitted.First.Stage.dta.2 <- left_join(Fitted.First.Stage.dta, First.Stage.dta, by = c("Entity")) %>% pivot_wider(names_from = term, values_from = estimate) %>% group_by(Entity) %>%
mutate(Fitted.GDP = Expenditure*GDP + `(Intercept)`,
Fitted.GDP.sqr = (Fitted.GDP)^2,
Fitted.GDP.cube = (Fitted.GDP)^3)
# Second Stage reg
Second.stage <- Fitted.First.Stage.dta.2 %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ Fitted.GDP+Fitted.GDP.sqr+
Fitted.GDP.cube+Year+Population, data = .))) %>%
unnest(Country.reg)
# Latex output
Second.stage.latex <- Second.stage %>% filter(Entity == "France" |
Entity == "Germany" |
Entity == "Italy" |
Entity == "United Kingdom")
#xtable(Second.stage.latex, digits=-3)
# Output
paged_table(Second.stage)\(~\)
4.2 - Fitted corrected from treatment
\(~\)
# Fitted Kuznets corrected
Second.stage.dta <- Second.stage %>% select(Entity, term, estimate)
Fitted.First.Stage.dta <- left_join(Fitted.First.Stage.dta, Second.stage.dta, by = c("Entity")) %>% pivot_wider(names_from = term, names_repair = "unique", values_from = estimate)
Fitted.Kuznets.corrected <- Fitted.First.Stage.dta %>% group_by(Entity) %>%
mutate(Fitted.Kuznets.corrected = Fitted.GDP*GDP + Fitted.GDP.sqr*(GDP^2) + Fitted.GDP.cube*(GDP^3)+`Year...12`*`Year...2`+`Population...13`*`Population...5`+`(Intercept)`)%>% select(Fitted.Kuznets.corrected, GDP)
# Output
paged_table(Fitted.Kuznets.corrected)\(~\)
4.2 - Plotting IV reg to observe the smoothing
\(~\)
Plot.IV <- Fitted.Kuznets.corrected %>%
filter(Entity == "France" |
Entity == "Germany" |
Entity == "Italy" |
Entity == "Spain") %>%
ggplot(aes(x=GDP, y=Fitted.Kuznets.corrected, color=Entity)) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray"))+
scale_color_manual(values = c("#34A0A4", "#D9ED92",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))+
facet_wrap(~Entity)
ggplotly(Plot.IV)\(~\)
4.3 - Regression with additional covariate Average.Schooling
\(~\)
# Data prep
dta.matching.Edu <- dta.matching %>%
filter(Year > 1990) %>%
select(CO2, GDP, GDP.sqr, GDP.cube,
Average.Schooling, Population, Year) %>%
na.omit()
# Education as additional covariate
reg.edu.time <- dta.matching.Edu %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube +
Average.Schooling + Population + Year, data = .))) %>%
unnest(Country.reg)
reg.edu <- dta.matching.Edu %>%
group_by(Entity) %>%
do(Country.reg = tidy(lm(CO2 ~ GDP + GDP.sqr + GDP.cube +
Average.Schooling + Population , data = .))) %>%
unnest(Country.reg)
# Latex output
reg.edu.time <- reg.edu.time %>% filter(Entity == "France" |
Entity == "Germany" |
Entity == "Italy" |
Entity == "United Kingdom")
reg.edu <- reg.edu %>% filter(Entity == "France" |
Entity == "Germany" |
Entity == "Italy" |
Entity == "United Kingdom")
#xtable(reg.edu.time, digits=-3)
#xtable(reg.edu, digits=-3)
# Output
paged_table(reg.edu.time)\(~\)
5.1 - Building treatment variable
\(~\)
plot.1.3 <- GOV_Expenditure %>%
na.omit()%>%
pivot_longer(cols = c("Expenditure", "Moving.average.3"),
names_to= "Spending", values_to = "values") %>%
rename(`Government spending` = values) %>%
filter(Country == "France" |
Country == "Germany" |
Country == "Italy" |
Country == "Spain" |
Country == "United Kingdom" |
Country == "Greece") %>%
ggplot(aes(x=TIME, y=`Government spending`, group=Country, color=Spending)) +
geom_line() +
geom_point(size=0.1) +
theme(axis.text = element_text(angle = 90),
panel.background = element_rect(fill = "white", colour = "gray"),
axis.line = element_line(size = 1, colour = "gray")) +
facet_wrap(~Country) +
scale_color_manual(values = c("#34A0A4", "#D9ED92",
"#99D98C", "#76C893",
"#52B69A", "#34A0A4"))
ggplotly(plot.1.3)