---
title: "Project"
author: "Thomas Reumiller"
date: "22 3 2023"
output:
  html_document: default
  pdf_document: default
editor_options:
  chunk_output_type: console
---

# Einfuehrung in das Paper "The Economic Effects of the Abolition of Serfdom: Evidence from the Russian Empire"
(Markevich, A., & Zhuravskaya, E. (2018))

https://doi.org/10.1257/aer.20160144

Das Paper befasst sich mit den Auswirkungen der Abschaffung der Leibeigenschaft im kaiserlichen Russland im Jahr 1861. Betrachtet werden einerseits die Veraenderung der Produktivitaet in der Landwirtschaft, als auch die Produktivitaetsaenderung in der Industrie. Als dritter Punkt wird die Gesundheit der landwirtschaftlichen Arbeiter untersucht. 

Die Ergebnisse zeigen, dass es zu signifikanten Verbesserungen in allen drei Punkten infolge der Abschaffung der Leibeigenschaft kam. 


# Code zur Reproduktion einiger Ergebnisse des Papers

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE) #Globales Setup um alle Chunks im Dokument anzuzeigen
```

## Laden aller benoetigten Bibliotheken

```{r,message = FALSE, warning=FALSE}
library(haven)        #Fuer das Einlesen von DTA files
library(dplyr)        #Fuer die Bearbeitung von Dataframes
library(ggplot2)      #Fuer Abbildungen
library(car)          #Fuer den Befehl avPlots
library(fixest)       #Fuer den feols-Befehl
library(modelsummary) #Zur Erstellung von Regressionstabellen
library(kableExtra)   #Zur Bearbeitung von Tabellen
```

## Einlesen der Daten

Lade country-Daten und zeige die ersten Zeilen:
```{r}
data_country = read_dta('replication_country_Serfdom.dta')
head(data_country)
```
Der Datensatz enthaelt die durchschnittliche Getreideproduktion der Jahre 1795-1914, sowie Dummys fuer das Vorhandensein von Produktionsdaten von Weizen und industriellem Output in den anderen Datensaetzen im betreffenden Jahr.


Lade district-Daten und zeige die ersten Zeilen:
```{r}
data_district = read_dta('replication_district_Serfdom.dta')
head(data_district)
```
Daten zu den verschiedenen Distrikten mit Fokus auf die eingezogenen Wehrpflichtigen.


Lade province-Daten und zeige die ersten Zeilen (hauptsaechlich verwendet):
```{r}
data_province = read_dta('replication_province_Serfdom.dta')
head(data_province)
```
Daten zu den verschiedenen Provinzen mit Fokus auf die Weizen- und Industrieproduktion, sowie den Anteil der Leibeigenen und die Verschuldung der Grundeigentuemer.

## Vorbereitungen fuer spaetere Regressionen

Erstelle Instrumentalvariable IVland fuer spaetere Regression. In spaeteren Regressionen wird diese Instrumentalvariable dazu verwendet, um die Umsetzung der Landreform zu modellieren. Die Variable beschreibt eine Interpolation zwischen 1-Verschuldung der Landeigentuemer und 1 im Intervall 1862-1882. Vor 1862 betraegt der Wert 0 und nach 1882 1. In westlichen Provinzen wird die Variable um 1963 direkt auf 1 gestellt, dies ist als Ergebnis der polnischen Revolution zu sehen.

```{r,message = FALSE, warning=FALSE}

#Erstelle Variable fuer maximale Verschuldung je Provinz

data_province = data_province %>% 
  group_by(id) %>% 
  mutate('pcollate1858pcpost1862'=max(pcollate1858pc,na.rm = TRUE)) %>% 
  ungroup()

#Bilde Interpolation zwischen 1- maximaler Verschuldung im Jahr 1862 und 0 im Jahr 1882. 

  data_province = mutate(data_province,'IVland' = case_when(
    
    year>=1863 & west==1 ~ 1,
    id==7 | id==8 | id==9 ~ 0,
    year==1862 ~ (1-pcollate1858pcpost1862),
    year==1863 ~ (1-pcollate1858pcpost1862)+(pcollate1858pcpost1862)/20,
    year==1864 ~ (1-pcollate1858pcpost1862)+2*(pcollate1858pcpost1862)/20,
    year==1865 ~ (1-pcollate1858pcpost1862)+3*(pcollate1858pcpost1862)/20,
    year==1866 ~ (1-pcollate1858pcpost1862)+4*(pcollate1858pcpost1862)/20,
    year==1867 ~ (1-pcollate1858pcpost1862)+5*(pcollate1858pcpost1862)/20,
    year==1868 ~ (1-pcollate1858pcpost1862)+6*(pcollate1858pcpost1862)/20,
    year==1869 ~ (1-pcollate1858pcpost1862)+7*(pcollate1858pcpost1862)/20,
    year==1870 ~ (1-pcollate1858pcpost1862)+8*(pcollate1858pcpost1862)/20,
    year==1871 ~ (1-pcollate1858pcpost1862)+9*(pcollate1858pcpost1862)/20,
    year==1872 ~ (1-pcollate1858pcpost1862)+10*(pcollate1858pcpost1862)/20,
    year==1873 ~ (1-pcollate1858pcpost1862)+11*(pcollate1858pcpost1862)/20,
    year==1874 ~ (1-pcollate1858pcpost1862)+12*(pcollate1858pcpost1862)/20,
    year==1875 ~ (1-pcollate1858pcpost1862)+13*(pcollate1858pcpost1862)/20,
    year==1876 ~ (1-pcollate1858pcpost1862)+14*(pcollate1858pcpost1862)/20,
    year==1877 ~ (1-pcollate1858pcpost1862)+15*(pcollate1858pcpost1862)/20,
    year==1878 ~ (1-pcollate1858pcpost1862)+16*(pcollate1858pcpost1862)/20,
    year==1879 ~ (1-pcollate1858pcpost1862)+17*(pcollate1858pcpost1862)/20,
    year==1880 ~ (1-pcollate1858pcpost1862)+18*(pcollate1858pcpost1862)/20,
    year==1881 ~ (1-pcollate1858pcpost1862)+19*(pcollate1858pcpost1862)/20,
    year>=1882 ~ 1,
    year<=1861 ~ 0,
  ))
```

Um in spaetere Regressionen die provinzspezifischen Trends einfliessen lassen zu koennen, werden die Namen ebenjener im folgenden Codeabschnitt aus den Daten herausgeloest und jeweils mit einem + verbunden, sodass sie an die Formel der ols-Regressionen angehaengt werden koennen.

```{r}
#colnames(data_province) zeigt die Namen, sowie die Position der einzelnen Spalten an

#speichere alle year_Iid-Variablen ab und verbinde sie mit einem + (wird spaeter in Regressionsformeln benoetigt)

writeClipboard(paste0('"',colnames(data_province)[76:141],'"', collapse=", "))

year_Iid = c("year_Iid_10", "year_Iid_11", "year_Iid_113", "year_Iid_114", "year_Iid_115", "year_Iid_12", "year_Iid_129", "year_Iid_13", "year_Iid_132", "year_Iid_135", "year_Iid_14", "year_Iid_141", "year_Iid_145", "year_Iid_147", "year_Iid_148", "year_Iid_15", "year_Iid_16", "year_Iid_17", "year_Iid_18", "year_Iid_19", "year_Iid_2", "year_Iid_20", "year_Iid_21", "year_Iid_22", "year_Iid_23", "year_Iid_24", "year_Iid_25", "year_Iid_26", "year_Iid_27", "year_Iid_28", "year_Iid_29", "year_Iid_3", "year_Iid_30", "year_Iid_31", "year_Iid_32", "year_Iid_33", "year_Iid_34", "year_Iid_35", "year_Iid_36", "year_Iid_37", "year_Iid_38", "year_Iid_39", "year_Iid_4", "year_Iid_40", "year_Iid_41", "year_Iid_42", "year_Iid_43", "year_Iid_44", "year_Iid_45", "year_Iid_46", "year_Iid_47", "year_Iid_48", "year_Iid_49", "year_Iid_5", "year_Iid_50", "year_Iid_51", "year_Iid_53", "year_Iid_54", "year_Iid_55", "year_Iid_56", "year_Iid_57", "year_Iid_58", "year_Iid_6", "year_Iid_7", "year_Iid_8", "year_Iid_9"
)

year_Iid.plus = paste0(year_Iid,collapse="+")
```

## Replikation einiger Abbildungen

Im Folgenden werden einige der Abbildungen aus dem Paper repliziert.

### Figure 1

Abbildung 1 zeigt die durchschnittliche Produktivitaet der Weizenproduktion in Russland ueber einen Zeitraum von 1801 bis 1914. Die vertikale rote Linie symbolisiert den Zeitpunkt der Abschaffung der Leibeigenschaft im Jahr 1861. 

```{r, warning=FALSE}  
data_country %>%                          
  filter(!is.na(rus_av_grain_prod)) %>% #Filtere NA Werte aus Datensatz
ggplot(aes(year,rus_av_grain_prod))+  #zeichne ein X-Y -Diagramm mir x = year und y = durschnittliche Weizenproduktion
  geom_line(size = 1,color = 'blue')+ #verbinde die Datenpunkte mit Linien
  geom_point(size = 2)+ #zeichne die Datenpunkte als Punkte
  geom_segment(aes(x = 1845, y = 5.78, xend = 1859, yend = 5.78), arrow = arrow(length = unit(0.5, "cm")),size = 1)+  #Einfuegen des Pfeils
  labs(x = "", y = "Average grain productivity, Russian empire", title = "")+ #Festlegen der Beschriftungen
  geom_vline(xintercept = 1861, color = "red", size = 1.2, linetype = "solid")+ #Einfuegen der vertikalen roten Linie
   annotate("text", x = 1820, y = 5.8, label = "The abolition of serfdom", 
           color = "black", size = 4, fontface = "bold")+         #Einfuegen des Labels der roten Linie
   scale_x_continuous(limits = c(1795, 1920), breaks = seq(1800, 1920, 20))+  #Skalierung und Limits der y-Achse
  scale_y_continuous(limits = c(2, 7))+ #Skalierung und Limits der y-Achse
   theme_classic() +
  theme(plot.title = element_text(size = 20, face = "bold"), #Textformt des Abbildungstitels
        plot.subtitle = element_text(size = 16, face = "bold"),#Textformt des Abbildungssubtitels
        axis.title = element_text(size = 14),#Textformt der Axenbeschriftung
        axis.text = element_text(size = 12),#Textformt der Axenelemente
        legend.position = "none",#Verbergen der Legende
        panel.grid.major.y = element_line(color = 'gray'), #Format der Gitterlinien entlang der y-Achse
        panel.grid.minor = element_blank(),#Format der Subgitterlinien
        panel.background = element_rect(fill = "white"))#Hintergrundfarbe
```

### Figure 4

In spaeteren Regressionen wird der Anteil der kloesterlichen Leibeigenen vor der Verstaatlichung als Instrument fuer den Anteil an Leibeigenen an der Gesamtbevoelkerung im Jahr 1858 verwendet.  Um den Zusammenhang dieser beiden Werte zu zeigen, eignet sich ein sogenannter added-variable-plot (oder auch conditional scatter plot genannt). Dieser zeigt die Korrelation zwischen einer unabhaengigen Variablen in einer Regression (hier Anteil der kloesterlichen Leibeigenen) und einer abhaengigen Variablen (hier Anteil an Leibeigenen im Jahr 1858), wobei die Effekte aller andern unabhaengigen Variablen herausgerechnet werden. (siehe http://fmwww.bc.edu/repec/bocode/a/avciplot_sj.pdf)

Der folgende added-variable-plot zeigt den Zusammenhang zwischen dem Anteil der kloesterlichen Leibeigenen und dem Anteil an Leibeigenen im Jahr 1858. In der Regression wird auf die Entfernung von Moskau und die Fruchtbarkeit des Landes getestet. Es werden beispielhaft Daten aus dem Jahre 1870 verwendet.

```{r}
#lm-Regression zur Bestimmung der Koeffizienten
reg_model = lm(sh_serfs_postEmancip ~ sh_nat_monasterialpost1861 + lndist_Moscowpost1861 + median_suitabilitypost1861, data = data_province, subset = year == 1870 & nobaltics & id != 44)

coef(reg_model)

#Added-variable-Plot
avPlot(reg_model, variable = 'sh_nat_monasterialpost1861',main = 'Conditional scatter plot and fitted line
Controls: log distance from Moscow and land suitability',xlab = 'Share of monasterial serfs before nationalization', ylab = 'Share of serfs in 1858' , id = FALSE) #zeige einen added-variable-Plot aus der durchgefuehrten Regression mit den gewuenschten Variablen. 

text(x = 0.101, y = 0.3,  labels = 'Coef = -1.25605',
     cex=1,col="black") #Anzeigen des Koeffizienten im Plot
text(x = 0.15, y = 0.25,  labels = '(robust) SE = 0.40275, t = -3.12',
     cex=1,col="black")#Anzeigen der Standardabweichung des Koeffizienten im Plot
```

Der Plot zeigt, dass das gewaehlte Instrument einen starken linearen Zusammenhang mit der abhaengigen Variablen gibt. Dies ist ein wichtiges Kriterium fuer die Auswahl eines passenden Instruments.
### Figure 5

Ein weiteres im Paper verwendetes Instrument ist die weiter oben erzeugte Variable IVland. Diese wird verwendet, um den Anteil an Arbeitern zu schaetzen, welche bereits Uebernahmevertraege fuer das Land unterzeichnet haben. Die folgende Abbildung zeigt wiederum in einem added-variable-plot den Zusammenhang zwischen diesen beiden Groessen. Es wird ebenfalls auf Entfernung von Moskau und Fruchtbarkeit des Landes getestet. Die beispielhaft verwendeten Daten stammen aus dem Jahre 1872

```{r}
#lm-Regression zur Bestimmung der Koeffizienten
reg_model = lm(sh_serfs1858p1861sh_buyout ~ IVland + lndist_Moscowpost1861 + median_suitabilitypost1861, data = data_province, subset = year == 1872 & nobaltics & id != 44)

coef(reg_model)

#Added-variable-Plot
avPlot(reg_model, variable = 'IVland' ,main = 'Conditional scatter plot and fitted line
Controls: log distance from Moscow and land suitability',xlab = 'Non-indebtedness instrument in 1872', ylab = 'Share of peasants with buyout contract signed in 1872' , id = FALSE) #zeige einen added-variable-Plot aus der durchgefuehrten Regression mit den gewuenschten Variablen.

text(x = 0.02, y = -0.38,  labels = 'Coef = 3.19862, (robust) SE = 1.20106, t = 2.66',
     cex=0.9,col="black")#Anzeigen des Koeffizienten, sowie der Standardabweichung des Koeffizienten im Plot
```

Auch hier zeigt sich ein passender linearer Zusammenhang, wobei dieser nicht ganz so klar ist wie bei der anderen Instrumentalvariablen.

### Figure A2

Abbildung A2 aus dem Anhang zeigt eindruecklich, wie sehr der Anteil der Leibeigenen an der Gesamtbevoelkerung abhaengt von der Entfernung zu Moskau. Die alten Provinzen (naeher an der Hauptstadt) besitzen einen wesentlich hoeheren Anteil, waehrend die weiter entfernten einen niedrigeren Anteil aufweisen. Aufgrund dieser starken Abhaengigkeit wird in den meisten durchgefuehrten Regressionen auf dieses Merkmal getestet (lndist_Moscowpost1861).

```{r}
summary(lm(sh_serfs_postEmancip~dist_Moscowpost1861,data = subset(data_province,year==1870 & nobaltics)))#ols-Zusammenhang zwischen dem Anteil der Leibeigenen und der Distanz von Moskau, beispielhaft fuer das Jahr 1970

ggplot(subset(data_province,year==1870 & nobaltics),aes(x = dist_Moscowpost1861, y = sh_serfs_postEmancip))+  #Plot des Zusammenahngs
  geom_point(aes(color = 'Share of serfs'),size = 2) +
  geom_text(aes(label = province_name,color = 'Share of serfs'),color = 'blue',hjust=-0.1, vjust=0)+  #einfuegen der Provinznahmen zu den Datenpunkten
  geom_smooth(aes(color = 'Fitted values'),method = "lm", se = FALSE) + #einfuegen der roten Regressionsgerade
  labs(title = "Geography of serfdom: the share of serfs in 1858 and the distance from Moscow", x = "Distance to Moscow from the centeroid of the province (km)", y = "The share of serfs in 1858", color = NULL) +
   scale_x_continuous(limits = c(0, 1500), breaks = seq(0, 1500, 500))+
  theme_classic() +
  theme(plot.title = element_text(size = 20, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.position = 'bottom',
        panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white"))
```

### Figure 6

Abbildung 6 zeigt die durchnittliche Produktivitaet in der Weizenproduktion in den verschiedenen Provinzen je Dekade. Um den unterschiedlichen Effekt auf Provinzen mit viel oder wenig voremanzipatorischem Anteil an Leibeigenschaft zu erkennen, wurden die Provinzen je nach Anteil der Leibeigenen im Jahr 1861 in die drei dargestellten Terzentiele unterteilt. Die vertikale rote Markierung gibt den Zeitpunkt der Abschaffung der Leibeigenschaft an. Durch diese Darstellung laesst sich ein absoluter Effekt der Abschaffung quantifizieren. 

```{r, warning=FALSE}
#Teile die Provinzen je einem der drei Terzentile zu
data = data_province
terc1 <- data %>%
  filter(grain_prod != "." & nobaltics & year == 1861) %>%  #beachte nur passende Daten im Jahre 1861
  mutate(terc_1 = ntile(sh_serfs1858, 3)) %>% #Teilt die Daten je nach Anteil an Leibeigenen in 3 Buckets auf
  group_by(id) %>%
  summarize(terc = terc_1)  #Gibt jeder Provinz-ID ein fixes Terzentil

#Teile jedem Datenpunkt abhaengig von seiner Provinz-ID ein Terzentil zu
data = data %>% 
  group_by(id) %>% 
  mutate(terc = case_when(
    id %in% filter(terc1, terc==1)$id~1,
    id %in% filter(terc1, terc==2)$id~2,
    id %in% filter(terc1, terc==3)$id~3,
  ))

#Teile jedem Datenpunkt eine Dekade zu
data = mutate(data,intervales=floor(((year-1)/10))*10+5)

#View(terc1)
#View(select(data,terc,intervales))

#Berechne die mittlere Produktivitaet je Terzentil und Dekade

grain_prod_below <- data %>%
  filter(grain_prod != "." & nobaltics & terc == 1) %>%
  group_by(intervales) %>%
  summarize(grain_prod_below = mean(grain_prod))

grain_prod_middle <- data %>%
  filter(grain_prod != "." & nobaltics & terc == 2) %>%
  group_by(intervales) %>%
  summarize(grain_prod_middle = mean(grain_prod))

grain_prod_above <- data %>%
  filter(grain_prod != "." & nobaltics & terc == 3) %>%
  group_by(intervales) %>%
  summarize(grain_prod_above = mean(grain_prod))

#Fuege alle Daten berechneten Daten zusammen
merged_data <- merge(grain_prod_below, grain_prod_middle, by = "intervales") %>%
  merge(grain_prod_above, by = "intervales")

#View(merged_data)

#plotte zusammengefuehrte Daten
merged_data %>%
  ggplot(aes(x = intervales)) +
  geom_line(aes(y = grain_prod_below, color = "1st tercile: range (0.1%−40%), mean 17%"), size = 1, linetype = 'dashed') +
  geom_point(aes(y = grain_prod_below, color = "1st tercile: range (0.1%−40%), mean 17%"), size = 3, shape = 16)+
  geom_line(aes(y = grain_prod_middle, color = "2nd tercile: range (41%−57%), mean 50%"), size = 1) +
  geom_point(aes(y = grain_prod_middle, color = "2nd tercile: range (41%−57%), mean 50%"), size = 4, shape = 18)+
  geom_line(aes(y = grain_prod_above, color = "3rd tercile: range (60%−83%), mean 69%"), size = 1, linetype = 'dashed') +
  geom_point(aes(y = grain_prod_above, color = "3rd tercile: range (60%−83%), mean 69%"), size = 3, shape = 15)+
  geom_vline(xintercept = 1861, color = "red", size = 0.7, linetype = "solid")+
  labs(x = "Decades",
       y = "Grain productivity, decade averages, raw data",
       title = "Productivity in provinces by terciles of the share of serfs and the emancipation", color = NULL) +
  scale_x_continuous(limits = c(1845, 1885), breaks = seq(1845, 1885, 10), labels = c("41-50", "51-60", "61-70", "71-80", "81-90")) +
  scale_y_continuous(limits = c(2.5, 4.5))+
   theme_classic() +
  theme(plot.title = element_text(size = 20, face = "bold"),
        plot.subtitle = element_text(size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.position = 'bottom',
        legend.direction = 'vertical',
        panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white"))
```

### Figure 7

Abbildung 7 dient der Untersuchung einer wichtigen Annahme des DiD-Ansatzes, welcher von den Autoren gewaehlt wurde. Und zwar soll untersucht werden, ob es divergierende Trends vor der Abschaffung in der landwirtschaftlichen Produktivitaet von Provinzen mit hoher und niedriger Leibeigenschaft gibt.

```{r}
data = data_province

# Erzeuge Interaktionsvariablen aus Dummys von 5 Jahresintervallen und dem Anteil Leibeigener
data$sh_serfs1858_1840 <- ifelse(data$year >= 1841 & data$year <= 1850, data$sh_serfs1858, 0)
data$sh_serfs1858_185155 <- ifelse(data$year >= 1851 & data$year <= 1855, data$sh_serfs1858, 0)
data$sh_serfs1858_185660 <- ifelse(data$year >= 1856 & data$year <= 1860, data$sh_serfs1858, 0)
data$sh_serfs1858_186165 <- ifelse(data$year >= 1861 & data$year <= 1865, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_186670 <- ifelse(data$year >= 1866 & data$year <= 1870, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_187175 <- ifelse(data$year >= 1871 & data$year <= 1875, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_187680 <- ifelse(data$year >= 1876 & data$year <= 1880, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_188185 <- ifelse(data$year >= 1881 & data$year <= 1885, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_188690 <- ifelse(data$year >= 1886 & data$year <= 1890, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_189195 <- ifelse(data$year >= 1891 & data$year <= 1895, data$sh_serfs_postEmancip, 0)
data$sh_serfs1858_1896 <- ifelse(data$year >= 1896, data$sh_serfs_postEmancip, 0)

# Ersetzt fehlende Varaibelen mit 0, sofern die urspruengliche Variable nicht fehlte
data$sh_serfs1858_1840[data$sh_serfs1858_1840 == "" & !is.na(data$sh_serfs1858)] <- 0
data$sh_serfs1858_185155[data$sh_serfs1858_185155 == "" & !is.na(data$sh_serfs1858)] <- 0
data$sh_serfs1858_185660[data$sh_serfs1858_185660 == "" & !is.na(data$sh_serfs1858)] <- 0
data$sh_serfs1858_186165[data$sh_serfs1858_186165 == "" & !is.na(data$sh_serfs_postEmancip)] <- 0
data$sh_serfs1858_186670[data$sh_serfs1858_186670 == "" & !is.na(data$sh_serfs_postEmancip)] <- 0
data$sh_serfs1858_187175[data$sh_serfs1858_187175 == "" & !is.na(data$sh_serfs_postEmancip)] <- 0
data$sh_serfs1858_187680[data$sh_serfs1858_187680 == "" & !is.na(data$sh_serfs_postEmancip)] <- 0


#Erstelle die fuer die Regression verwendeten Interaktionsvariablen
data = data %>% 
  mutate(D1=sh_serfs1858_1840) %>%
  mutate(D2=sh_serfs1858_185155) %>%
  mutate(D3=sh_serfs1858_185660) %>%
  mutate(D4=sh_serfs1858_186165) %>%
  mutate(D5=sh_serfs1858_186670) %>%
  mutate(D6=sh_serfs1858_187175) %>%
  mutate(D7=sh_serfs1858_187680) %>%
  mutate(D8=sh_serfs1858_188185) %>%
  mutate(D9=sh_serfs1858_188690) %>%
  mutate(D10=sh_serfs1858_189195) %>%
  mutate(D11=sh_serfs1858_1896)


#Teile die Zeitintervalle ein
data = data %>%
  mutate(time = 0) %>% 
  mutate(time = case_when(
    year>=1830 & year<=1840 ~ 0.5,
    year>=1841 & year<=1850 ~ 1,
    year>=1851 & year<=1855 ~ 2,
    year>=1856 & year<=1860 ~ 3,
    year>=1861 & year<=1865 ~ 4,
    year>=1866 & year<=1870 ~ 5,
    year>=1871 & year<=1875 ~ 6,
    year>=1876 & year<=1880 ~ 7,
    year>=1881 & year<=1885 ~ 8,
    year>=1886 & year<=1890 ~ 9,
    year>=1891 & year<=1895 ~ 10,
    year>=1896 & year<=1900 ~ 11,
    year>=1900 ~ 11.5,
  ))


#Regressionsformel
ols.model5 = "grain_prod ~ D1+D2+D3+D4+D5+D6+D7+D8+D9+D10+D11+sh_stpeasants1858post1866+sh_royal1858post1859+median_suitabilityDemeanpost1861+lndist_MoscowDemeanpost1861+factor(year)+ "

#Fuege Regionale Trends hinzu

#writeClipboard(paste0('"',colnames(data)[142:156],'"', collapse=", ")) 
l5trend = c("l5trendAstOren", "l5trendBaltic", "l5trendBelLit", "l5trendCBER", "l5trendCIRNoMsc", "l5trendLeftBank", "l5trendLowVolgaNA", "l5trendMidVolga", "l5trendNorth", "l5trendNorthwestNoSPb", "l5trendRightBank", "l5trendSouth", "l5trendUralsNoO", "l5trendWest", "l5trendcapitals"

)

l5trend.plus = paste0(l5trend,collapse="+")

ols.model5 = formula(paste0(ols.model5, l5trend.plus,'|id'))

#Fuere Regression aus
ols5 = feols(ols.model5, data = subset(data,nobaltics== TRUE),cluster = 'idref') 

#summary(ols5)

model = ols5

#Fuehre Ergebnisse in Tabelle zusammen und berechne Konfidenzintervall
coef_df = data.frame(time = 1:11)
coef_df = coef_df %>% 
  mutate(coef = coef(model)[1:11],
            low = coef(model)[1:11] - 1.64 * sqrt(diag(vcov(model))[1:11]),
            high = coef(model)[1:11] + 1.64 * sqrt(diag(vcov(model))[1:11]))

# Plotter Ergebnisse
ggplot(coef_df, aes(x = time, y = coef)) +
  geom_line(color = "blue", size = 1.2) +
  geom_point(color = "blue", size = 3) +
  geom_line(aes(y = low),color = "grey", size = 1.2, linetype = 'dashed') +
  geom_line(aes(y = high),color = "grey", size = 1.2, linetype = 'dashed') +
  geom_vline(xintercept = 3.5, color = "red", size = 1.2, linetype = "solid") +
  scale_x_continuous(breaks = 1:11, labels = c("41-50", "51-55", "56-60", "61-65", "66-70", "71-75", "76-80", "81-85", "86-90", "90-95", "post 96")) +
  labs(title = "Grain productivity", subtitle = "(relative to years 1795-1829)", x = "Time", y = "Coefficient") +
  theme_classic() +
  theme(plot.title = element_text(size = 20, face = "bold"),
        plot.subtitle = element_text(size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.position = "none",
        panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white"))

```

In den Jahren vor der Abschaffung laesst sich kein Trend erkennen, dies ist ein gutes Zeichen fuer den DiD-Ansatz. Nach der Emanzipation zeigt sich ein nicht signifikanter positiver Trend. Dieser bedeutet, dass Provinzen mit ehemals grossem Anteil an Leibeigenen etwas zu den anderen Provinzen aufschliessen.

## Replikation der Ergebnistabellen

Im folgenden Abschnitt werden die zur Erstellung der Ergebnistabellen erforderlichen Regressionen durchgefuehrt, sowie die Tabellen erstellt. Hierbei werden die Paneldaten mit Fixed-Effect-Regressionen verarbeitet und die Fehler nach der Provinz-ID (idref) vor bzw. nach der Emanzipation geclustert. Der ssc Input in der Funktion feols dient dazu, die Fehlerberechnung an die Berechnungsart des Stata Befehls xtreg anzupassen. Verwendet werden jeweils nur Provinzen, welche nicht zu den baltischen Provinzen gehoeren.

### Table 2

Tabelle 2 stellt die Ergebnisse der Untersuchungen zur Produktivitaetssteigerung in der Landwirtschaft dar. Geschaetzt wird jeweils der Einfluss der Interaktion zwischen dem Anteil an Leibeigenen und einer Dummyvariable fuer die Abschaffung der Leibeigenschaft (sh_serfs_postEmancip) an der Produktivitaet (grain_prod). Die unterschiedlichen Regressionen beziehen jeweils veraenderte Kontrollvariablen mit ein bzw. verwenden Instrumentalvariablen.

#### Table 2 OLS 1

Regression 1 von Tabelle 2 fuehrt die Basic-Regression fuer die Produktivitaet durch und korrigiert nur auf die fixed-effects.

```{r,message = FALSE, warning=FALSE}
ols.model21 = formula("grain_prod ~ sh_serfs_postEmancip + factor(year)| id") #verwendete Formel der Regression
ols21 = feols(ols.model21, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols21)     #Ausgabe der Regressionsergebnisse

#1-ols21$ssr/ols21$ssr_fe_only    #Beispielhafte Berechnung des R2-Wertes bei normaler OLS
```


#### Table 2 OLS 2

Regression 2 von Tabelle 2 prueft nun zusaetzlich auf provinzspezifische Trends (year_Iid), die logarithmische Distanz von der Hauptstadt Moskau (sieh Abbildung A2) und die Fruchtbarkeit des jeweiligen Landes.

```{r,message = FALSE, warning=FALSE}
ols.model22 = "grain_prod ~ sh_serfs_postEmancip + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model22 = formula(paste0(ols.model22, year_Iid.plus,"|id")) #Hinzufuegen der zu provinzspezifischen Trends

ols22 =  feols(ols.model22, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols22)
```

#### Table 2 OLS 3

Regression 3 von Tabelle 2 schaetzt den bekannten Zusammenhang nun mittels einer Instrumentalvariable fuer "sh_serfs_postEmancip". Hierfuer wird der Anteil an kloesterlichen Leibeigenen verwendet (siehe Abbildung 4).

```{r,message = FALSE, warning=FALSE}
ols.model23 = "grain_prod ~ lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + factor(id)+"
ols.model23 = formula(paste0(ols.model23, year_Iid.plus,'|id|sh_serfs_postEmancip ~ sh_nat_monasterialpost1861' ))

ols23 = feols(ols.model23, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols23)
#1-ols23$ssr/ols23$ssr_null #Beispielhafte Berechnung des R2-Wertes bei  IV-Regression
```

#### Table 2 OLS 4

Regression 4 von Tabelle 2 ist wieder eine normale OLS-Regression. Hier wird jedoch zusaetzlich noch auf die Anteile an staatlichen Bauern nach 1866 (sh_stpeasants1858post1866) und koeniglichen Landarbeitern nach 1859 (sh_royal1858post1859) getestet.

```{r,message = FALSE, warning=FALSE}
ols.model24 = "grain_prod ~ sh_serfs_postEmancip + sh_stpeasants1858post1866 + sh_royal1858post1859 + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model24 = formula(paste0(ols.model24, year_Iid.plus,'|id'))

ols24 =  feols(ols.model24, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols24)
```

#### Table 2 OLS 5

Regression 5 von Tabelle 2 prueft dieses Mal nicht auf den Anteil der Leibeigenen nach der Abschaffung, sondern jeweils den Anteil der Leibeigenen zwischen 1861 und 1870 (sh_serfs_1861_1870), sowie nach 1871 (sh_serfs_1871_1900). Dies ermoeglicht eine seperate Betrachtung des Effekts in der ersten Dekade und dem restlichen 19. Jahrhundert.

```{r,message = FALSE, warning=FALSE}
ols.model25 = "grain_prod ~ sh_serfs_1861_1870 + sh_serfs_1871_1900 + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model25 = formula(paste0(ols.model25, year_Iid.plus,'|id'))

ols25 =  feols(ols.model25, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols25)
```

#### Table 2 OLS 6

Regression 6 von Tabelle 2 verwendet wieder den Anteil der Leibeigenen nach der Abschaffung. Dieses Mal prueft die Regression jedoch auch auf den Anteil der Landarbeiter mit abgeschlossenen Abloesevertraegen (sh_serfs1858p1861sh_buyout). Dies ermoeglicht eine Aufspaltung des Gesamteffektes auf die beiden Teileffekte der Emanzipation und der Landreform.

```{r,message = FALSE, warning=FALSE}
ols.model26 = "grain_prod ~ sh_serfs_postEmancip + sh_serfs1858p1861sh_buyout + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model26 = formula(paste0(ols.model26, year_Iid.plus,'|id'))

ols26 =  feols(ols.model26, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols26)
```

#### Table 2 OLS 7

Regression 7 von Tabelle 2 ist eine IV-Regression, welche wie in OLS 3 fuer den Anteil an Leibeigenen nach der Abschaffung den Anteil der kloesterlichen Leibeigenen als Instrument nutzt und zusaetzlich fuer den Anteil  der Landarbeiter mit abgeschlossenen Abloesevertraegen als Instrument IVland (siehe oben) verwendet.

```{r,message = FALSE, warning=FALSE}
ols.model27 = "grain_prod ~ lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + factor(id)+"
ols.model27 = formula(paste0(ols.model27, year_Iid.plus,'|id|sh_serfs_postEmancip+sh_serfs1858p1861sh_buyout~sh_nat_monasterialpost1861+IVland'))

ols27 = feols(ols.model27, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols27)

#1-ols27$ssr/ols27$ssr_null
```

#### Table 2 Panel A

In den folgenden Codeabschnitten werden die Ergebnisse der verschiedenen Regressionen zusammengefuehrt.

```{r,message = FALSE, warning=FALSE}
models = list(ols21, ols22, ols23, ols24, ols25, ols26, ols27) #alle in die Tabelle aufgenommenen Regressionen

#Reihen in der Tabelle die manuell hinzugefuegt werden
rows = data.frame(
  c('Year and province fixed effects','Province-
specific trends','State and royal peasant reforms','Observations','R2'),
  c('Yes','No','No',ols21$nobs,round(r2(ols21, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols22$nobs,round(r2(ols22, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols23$nobs,round(r2(ols23, type = "r2"), digits = 3)),
  c('Yes','Yes','Yes',ols24$nobs,round(r2(ols24, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols25$nobs,round(r2(ols25, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols26$nobs,round(r2(ols26, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols27$nobs,round(r2(ols27, type = "r2"), digits = 3)))

table2a = modelsummary(models,
             coef_map = c( #Auflistung aller dargestellter Koeffizienten mit entsprechenden Bezeichnungen
'fit_sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation',
'sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation',
'sh_serfs_1861_1870'= 'Share of serfs × 1861-1870',
'sh_serfs_1871_1900'= 'Share of serfs × post-1871',
'fit_sh_serfs1858p1861sh_buyout' = 'Share of peasants with signed buyout contracts',
'sh_serfs1858p1861sh_buyout' = 'Share of peasants with signed buyout contracts',
'lndist_MoscowDemeanpost1861' = 'Demeaned log distance to Moscow x Post-emancipation',
'median_suitabilityDemeanpost1861' = 'Demeaned crop suitability x Post-emancipation'),

             add_rows = rows,  #manuelle hinzufuegung der oben genannten Reihen
             stars = FALSE,    #Verbergen der Sterne zum anzeigen der Signifikanz
             gof_map = "none", #Verbergen der uebrigen statistischen Angaben

fmt = function(x) round(x, digits = 2),  #Festlegung der signifikanten Stellen
title = 'Table 2—The Effect of the Abolition of Serfdom on Productivity in Agriculture, Panel A. Panel data estimation'   #Titel der Tabelle

) %>% 
  add_header_above(c("Model:" = 1, "OLS" = 1, "OLS" = 1, "IV, second
stage" = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1, "IV, second
stage" = 1)) %>%  #Hinzufuegen von Ueberschriften
  add_header_above(c(" " = 1, "Grain productivity" = 7))
  
table2a #Ausgeben der Tabelle

#save_kable(table2a,'Tabel 2, Panel A.png')

```

Die Ergebnisse in Tabelle 2 zeigen über alle Regressionen hinweg einen positiven signifikanten Zusammenhang zwischen der Emanzipation und der Produktivitaet in der Landwirtschaft. Die Ergebnisse aus OLS 5 legen nahe, dass 3/4 des gesamten Effektes der Emanzipation bereits in der ersten Dekade zu Tragen gekommen ist. Ebenso zeigen die Ergebnisse aus OLS 3 vs OLS 7, dass der reine Effekt aus der Emanzipation 2,1 mal so gross wie der gesamte Effekt waere. Dies liegt daran, dass die Landreform einen negativen Einfluss auf die Entwicklung der Produktivitaet hatte.

#### Table 2 Panel B

1st-stage-Ergebnisse der IV-Regressionen zusammengetragen in einer Tabelle.

```{r, message = FALSE, warning=FALSE}
models = list("(3)"=ols23$iv_first_stage$sh_serfs_postEmancip,
              "(7.1)" = ols27$iv_first_stage$sh_serfs_postEmancip,
              "(7.2)" = ols27$iv_first_stage$sh_serfs1858p1861sh_buyout)

rows = data.frame(
  c('Controls as in respective column of panel A'),
  c('Yes'),
  c('Yes'),
  c('Yes'))

attr(rows, 'position') = c(5)  #Position der manuell hinzugefuegten Zeilen

table2b = modelsummary(models,
             coef_map = c(
'sh_nat_monasterialpost1861' =  'Share of nationalized monasterial serfs ×
Post-emancipation',
'IVland' =  'Interpolation between (1-
indebtedness) and 1 in the interval 1862–
1882
Post-emancipation'),

             add_rows = rows,
             stars = FALSE, 
             gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = function(x) format(round(x, 3), big.mark=""))),

fmt = function(x) round(x, digits = 2),
title = 'Table 2—The Effect of the Abolition of Serfdom on Productivity in Agriculture, Panel B. First stages of the corresponding 2SLS panel regressions'

)  %>% 
  add_header_above(c("Model:" = 1, "IV, first stage" = 1, "IV, first stage" = 1, "IV, first stage" = 1)) %>% 
  add_header_above(c(" " = 1, "Share of serfs × Post-
emancipation" = 1, "Share of
serfs
× Postemancipation" = 1, "Grain productivityShare of
peasants with
signed buyout
contracts" = 1)) 
 
  table2b
```


### Table 3

Tabelle 3 versucht einige Mechanismen hinter der Produktivitaetssteigerung in der Landwirtschaft aufzudecken. Konkret wird hierbei zum einen auf Aspekte der Abschaffung der Leibeigenschaft, als auch auf Aspekte bezueglich der Landreform eingegangen. Ausgangsvariablen sind hierbei nicht mehr nur die Getreideproduktivitaet insgesamt, sondern auch der Anteil von Winter- bzw. Fruehlingsgetreide. Dies soll zeigen, dass sowohl Anreizverbesserungen ein Rolle bei der Verbesserung der Produktivitaet spielen, als auch dass sich die Landarbeiter bei eigener Entscheidungsgewalt besser den Marktbedingungen anpassen konnten, um ihre Profite zu maximieren.

#### Table 3 OLS 1

Regression 1 von Tabelle 3 schaetzt wie in den Regressionen zu Tabelle 2 nochmals die Produktivitaet in der Weizenproduktion anhand dem Anteil an Leibeigenen nach der Emanzipation. Es werden wieder auf provinzspezifische Trends getestet, als auch auf den Anteil an abgeschlossenen Abloesevertraegen. Zusaetzlich wir hier noch auf eine Interaktion mit einer Dummyvariable getestet, ob die betreffende Provinz von Umverteilungen betroffen war (shserf1858p1861shbuy_RC).

```{r}
ols.model31 = "grain_prod ~ sh_serfs_postEmancip + sh_serfs1858p1861sh_buyout + shserf1858p1861shbuy_RC + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model31 = formula(paste0(ols.model31, year_Iid.plus,'|id'))

ols31 =  feols(ols.model31, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols31)
```

#### Table 3 OLS 2

Regression 2 von Tabelle 3 verfolgt die selbe Regression wie Regression 1, nur wird hierbei nicht auf Betroffenheit von Umverteilung, sondern auf den Abschluss von impliziten Arbeitsvertraegen getestet (sh_serfs1858p1861Chartp63).

```{r}
ols.model32 = "grain_prod ~ sh_serfs_postEmancip + sh_serfs1858p1861sh_buyout + sh_serfs1858p1861Chartp63 + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model32 = formula(paste0(ols.model32, year_Iid.plus,'|id'))

ols32 =  feols(ols.model32, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols32)
```

#### Table 3 OLS 3

Regression 3 von Tabelle 3 schaetzt den Anteil an Wintergetreide (sh_wgrains) am gesamt angebauten Getreide.Die erklaerenden Variablen sind wiederum der Anteil der Leibeigenen getestet auf provinzspezifische Trends, den Abstand von Moskau, die Fruchtbarkeit, sowie die gemessene Temperatur im Vorjahr (lagtemp_mean_e400kmDemean) und der Interaktionsterm zwischen dieser Temperatur und dem Dummy fuer die Abschaffung (lagtempDemagrsh_pr1858post1861).

```{r}
ols.model33 = "sh_wgrains ~ sh_serfs_postEmancip + lagtempDemagrsh_pr1858post1861 + lagtemp_mean_e400kmDemean + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model33 = formula(paste0(ols.model33, year_Iid.plus,'|id'))

ols33 =  feols(ols.model33, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols33)
```

#### Table 3 OLS 4

Regression 4 von Tabelle 3 verfolgt dieselbe Regression, testet jedoch nicht auf die Temperatur, sondern auf eine Interaktion zwischen dem Anteil der Leibeigenen, dem Dummy zur Abschaffung und dem Verhaeltnis vom Weltmarktpreis von Weizen zu Hafer im Vorjahr (sh_serfs1858p1861_lag_rye_wheat).

```{r}
ols.model34 = "sh_wgrains ~ sh_serfs_postEmancip + sh_serfs1858p1861_lag_rye_wheat + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model34 = formula(paste0(ols.model34, year_Iid.plus,'|id'))

ols34 =  feols(ols.model34, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols34)
```

#### Table 3 OLS 5

Regression 5 von Tabelle 3 verknuepft nun die beiden Ansaetze aus den vorhergegangenen Regressionen.

```{r}
ols.model35 = "sh_wgrains ~ sh_serfs_postEmancip + sh_serfs1858p1861_lag_rye_wheat + lagtempDemagrsh_pr1858post1861 + lagtemp_mean_e400kmDemean + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model35 = formula(paste0(ols.model35, year_Iid.plus,'|id'))

ols35 =  feols(ols.model35, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols35)
```

#### Table 3 OLS 6

Regression 6 von Tabelle 3 testet anstatt auf Temperatur und Weizenpreise auf lokale Haferpreise vor (sh_serfs1858_lag_rye_oats) bzw. nach (sh_serfs1858p1861_lag_rye_oats) der Abschaffung der Leibeigenschaft, sowie auf das Verhaeltnis von Weizen zu Haferpreis (lagryeoatpricerrratiodemean).

```{r}
ols.model36 = "sh_wgrains ~ sh_serfs_postEmancip + sh_serfs1858p1861_lag_rye_oats + lagryeoatpricerrratiodemean + sh_serfs1858_lag_rye_oats + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model36 = formula(paste0(ols.model36, year_Iid.plus,'|id'))

ols36 =  feols(ols.model36, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols36)
```

#### Table 3 OLS 7

Regression 7 von Tabelle 3 fuegt der Regression 6 noch zusaetzlich die Tests auf Temperatur hinzu. 

```{r}
ols.model37 = "sh_wgrains ~ sh_serfs_postEmancip + sh_serfs1858p1861_lag_rye_oats + lagryeoatpricerrratiodemean + sh_serfs1858_lag_rye_oats + lagtempDemagrsh_pr1858post1861 + lagtemp_mean_e400kmDemean + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model37 = formula(paste0(ols.model37, year_Iid.plus,'|id'))

ols37 =  feols(ols.model37, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols37)
```

#### Table 3

In den folgenden Codeabschnitten werden die Ergebnisse der verschiedenen Regressionen zusammengefuehrt.

```{r}
models = list(ols31,ols32,ols33,ols34,ols35,ols36,ols37)
rows = data.frame(c('Year and province fixed effects','Province-
specific trends'),c('Yes','Yes'),c('Yes','Yes'),c('Yes','Yes'),c('Yes','Yes'),c('Yes','Yes'),c('Yes','Yes'),c('Yes','Yes'))
attr(rows, 'position') = c(25:26)

table3 = modelsummary(models,
             coef_map = c(
'sh_serfs_postEmancip'= 'Share of serfs × Post-emancipation',
'sh_serfs1858p1861sh_buyout' = 'Share of peasants with signed buyout contracts',
'shserf1858p1861shbuy_RC' = 'Share of peasants with signed buyout contract x repartition commune',
'sh_serfs1858p1861Chartp63' = 'Share of serfs × Post-emancipation × Implicit contracts',
'lagtemp_mean_e400kmDemean' = 'Demeaned temperature (t−1)',
'lagtempDemagrsh_pr1858post1861' = 'Share of serfs × Post-emancipation × Demeaned temperature (t−1)',
'sh_serfs1858p1861_lag_rye_wheat' = 'Share of serfs × Post-emancipation × Demeaned rye-to-wheat world price ratio (t−1)',
'sh_serfs1858p1861_lag_rye_oats' = 'Share of serfs × Post-emancipation × Demeaned rye-to-oat local price ratio (t−1)',
'sh_serfs1858_lag_rye_oats' = 'Share of serfs × Demeaned rye-to-oat local price ratio (t−1)',
'lagryeoatpricerrratiodemean' = 'Demeaned rye-to-oat local price ratio (t−1)',
'lndist_MoscowDemeanpost1861' = 'Demeaned log distance to Moscow x Post-emancipation',
'median_suitabilityDemeanpost1861' = 'Demeaned crop suitability x Post-emancipation'),
             add_rows = rows,
             stars = FALSE,
             gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = function(x) format(round(x, 3), big.mark="")),
                            list("raw" = "r2.within", "clean" = "R2", "fmt" = function(x) format(round(x, 3), big.mark=""))),
title = 'Table 3—The Mechanisms behind the Effects of the Land Reform and the Emancipation'

) %>% 
  add_header_above(c(" " = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1, "OLS" = 1)) %>% 
  add_header_above(c("Dependent variable:" = 1, "Grain productivity" = 2, "Share of winter crops seeded at t−1 in total winter and spring crops seeded at [t−1;t] production cycle" = 5))

table3

#save_kable(table3,'Tabel 3.png')
```

Die Ergebnisse von Tabelle 3 unterstuetzen die Hypothese, dass ein hoeherer Anreiz und bessere lokale Entscheidungen Grund fuer das Ansteigen der Produktivitaet sein koennten. Ebenso wird die Annahme gestuetzt, dass die Umverteilungskommunen einen negativen Effekt hatten. 

### Table 4

Tabelle 4 beschaeftigt sich mit der Produktivitaet im Industriesektor. Die Regressionsformeln sind im Grunde identisch mit denen in Tabelle 2 nur ist die zu erklaerende Variable nicht mehr die Produktivitaet in der Landwirtschaft, sondern der logarithmierte industrielle Output.

#### Table 4 OLS 1

Regression 1 von Tabelle 4 ist wieder die Standardregression, welche als erklaerende Variable den Anteil der Leibeigenen nach der Abschaffung verwendet. 

```{r}
ols.model41 = "ln_indoutput_adinf ~ sh_serfs_postEmancip + factor(year)"
ols.model41 = formula(paste0(ols.model41,'|id'))

ols41 =  feols(ols.model41, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols41)
```

#### Table 4 OLS 2

Regression 2 von Tabelle 4 bezieht regionale Trends mit ein. Die im Code genannten "l5trend" sind regionale Trends, welche sich auf 15 grosse Regionen beziehen. Es wird wiederum auf den Abstand von Moskau und die Fruchtbarkeit getestet.

```{r}
ols.model42 = "ln_indoutput_adinf ~ sh_serfs_postEmancip + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model42 = formula(paste0(ols.model42,l5trend.plus ,'|id'))

ols42 =  feols(ols.model42, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols42)
```

#### Table 4 OLS 3

Regression 3 von Tabelle 4 schaetzt den Zusammenhang aus Regression 2 nun mit dem bekannten Instrument fuer den Anteil an Leibeigenen, naemlich dem Anteil der kloesterlichen Leibeigenen.

```{r}
ols.model43 = "ln_indoutput_adinf ~ lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + factor(id)+"
ols.model43 = formula(paste0(ols.model43, l5trend.plus,'|id|sh_serfs_postEmancip~sh_nat_monasterialpost1861'))

ols43 = feols(ols.model43, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols43)

#1-ols43$ssr/ols43$ssr_null
```

#### Table 4 OLS 4

Regression 4 von Tabelle 4 schaetzt den gewohnten Zuammenhang, testet jedoch nicht nur auf regionspezifische Trends, sondern auch auf den Anteil der staatlichen Leibeigenen nach 1866, sowie den koeniglichen Leibeigenen nach 1859. 

```{r}
ols.model44 = "ln_indoutput_adinf ~ sh_serfs_postEmancip + sh_stpeasants1858post1866 + sh_royal1858post1859 + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + "
ols.model44 = formula(paste0(ols.model44,l5trend.plus ,'|id'))

ols44 =  feols(ols.model44, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols44)
```

#### Table 4 Panel A

In den folgenden Codeabschnitten werden die Ergebnisse der verschiedenen Regressionen zusammengefuehrt.

```{r}
models = list(ols41, ols42, ols43, ols44)

rows = data.frame(
  c('Year and province fixed effects','Region-
specific trends','State and royal peasant reforms','Observations','R2'),
  c('Yes','No','No',ols41$nobs,round(r2(ols41, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols42$nobs,round(r2(ols42, type = "wr2"), digits = 3)),
  c('Yes','Yes','No',ols43$nobs,round(r2(ols43, type = "r2"), digits = 3)),
  c('Yes','Yes','Yes',ols44$nobs,round(r2(ols44, type = "wr2"), digits = 3)))

table4a = modelsummary(models,
             coef_map = c(
'fit_sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation',
'sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation',
'lndist_MoscowDemeanpost1861' = 'Demeaned log distance to Moscow x Post-emancipation',
'median_suitabilityDemeanpost1861' = 'Demeaned crop suitability x Post-emancipation'),
             add_rows = rows,
             stars = FALSE, 
             gof_map = "none",
fmt = function(x) round(x, digits = 2),
title = 'Table 4—The Abolition of Serfdom and Industrial Development, Panel A. Panel data estimation'

) %>% 
  add_header_above(c(" " = 1, "OLS" = 1, "OLS" = 1, "IV, second
stage" = 1, "OLS" = 1)) %>% 
  add_header_above(c(" " = 1, "ln(industrial output)" = 4))
  
table4a

#save_kable(table4a,'Tabel 4, Panel A.png')
```

Auch hier zeigen sich deutliche positive signifikante Ergebnisse. Jedoch weicht hier das Ergebnis der OLS zur IV Regression stark ab. Dies liegt wohl an heterogenen Effekten der Abschaffung auf industriellen Output.

#### Table 4 Panel B

1st-stage-Ergebnisse der IV-Regression zusammengetragen in einer Tabelle.

```{r}
models = list("(3)"=ols43$iv_first_stage$sh_serfs_postEmancip)

rows = data.frame(
  c('Controls as in respective column of panel A'),
  c('Yes'))

attr(rows, 'position') = c(3)

table4b = modelsummary(models,
             coef_map = c(
'sh_nat_monasterialpost1861' =  'Share of nationalized monasterial serfs ×
Post-emancipation'),
             add_rows = rows,
             stars = FALSE, 
             gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = function(x) format(round(x, 3), big.mark=""))),
fmt = function(x) round(x, digits = 2),
title = 'Table 4—The Abolition of Serfdom and Industrial Development, Panel B. First stages of the corresponding 2SLS panel regressions'

)  %>% 
  add_header_above(c("Model:" = 1, "IV, first stage" = 1)) %>% 
  add_header_above(c(" " = 1, "Share of serfs × Post-
emancipation" = 1)) 

  table4b
```


### Table 5

Tabelle 5 praesentiert schliesslich die Ergebnisse aus der Untersuchung der Gesundheitsentwicklung der Landarbeiter nach der Abschaffung der Leibeigenschaft. Da allgemeine Gesundheit schwer direkt zu messen ist, wird hierfuer die Groesse der wehrpflichtigen Jugend als Schaetzer verwendet. 

#### Table 5 OLS 1

Regression 1 von Tabelle 5 vollzieht die Basisversion der Regression. Hierbei ist die Groesse der Wehrpflichtigen die zu erklaerende Variable und der Anteil der Leibeigenen nach der Abschaffung ebenjener die erklaerende Variable. Getestet wird auf jahrspezifische Trends.

```{r}
ols.model51 = "height ~ sh_serfs_postEmancip + factor(year)"
ols.model51 = formula(paste0(ols.model51,'|id'))

ols51 =  feols(ols.model51, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols51)
```

#### Table 5 OLS 2

Regression 2 von Tabelle 5 testet nun zusätzlich auf die Enfernung von Moskau und Fruchtbarkeit.

```{r}
ols.model52 = "height ~ sh_serfs_postEmancip + lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) "
ols.model52 = formula(paste0(ols.model52 ,'|id'))

ols52 =  feols(ols.model52, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE, fixef.K="none"))

#summary(ols52)
```

#### Table 5 OLS 3

Regression 3 von Tabelle 5 schaetzt den Zusammenhang mit dem bekannten Instrument sh_nat_monasterialpost1861 fuer den Anteil der Leibeigenen.

```{r}
ols.model53 = "height ~ lndist_MoscowDemeanpost1861 + median_suitabilityDemeanpost1861 + factor(year) + factor(id)"
ols.model53 = formula(paste0(ols.model53,'|id|sh_serfs_postEmancip~sh_nat_monasterialpost1861'))

ols53 = feols(ols.model53, data = subset(data_province, nobaltics == TRUE),cluster = 'idref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols53)

#1-ols53$ssr/ols53$ssr_null
```

#### Table 5 OLS 4

Regression 4 von Tabelle 5 bedient sich nun nicht mehr den Daten aus dem Provinzdatensatz, sondern verwendet die kleinraeumigeren Distriktdaten. Hiermit werden die selben Regression wie oben nochmal durchgefuehrt. Einzig die Variablennamen sind etwas veraendert. 
 "height"       steht fuer die Groesse der Wehrpflichtigen
 "sh_allserfspost1861"    steht fuer den Anteil der Leibeigenen an der Gesamtbevoelkerung
 "yearbirth"    steht fuer das Geburtsjahr
 "iddistrict"   ist die Distrikt-ID

```{r}
ols.model54 = "height ~ sh_allserfspost1861 + factor(yearbirth)"
ols.model54 = formula(paste0(ols.model54,'|iddistrict'))

ols54 = feols(ols.model54, data = subset(data_district,nobaltics & !moscow & alldrafted>30),cluster = 'iddistref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols54)

```

#### Table 5 OLS 5

Regression 5 von Tabelle 5 testet nun wiederum auf den Abstand zu Moskau und Fruchtbarkeit.

```{r}
ols.model55 = "height ~ sh_allserfspost1861 + lndistMoscowDemeanpost1861 + median_suitDemeanpost1861 + factor(yearbirth)"
ols.model55 = formula(paste0(ols.model55,'|iddistrict'))

ols55 = feols(ols.model55, data = subset(data_district, (nobaltics==TRUE) & (moscow != TRUE) & alldrafted>30),cluster = 'iddistref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols55)
```

#### Table 5 OLS 6

Regression 6 von Tabelle 5 verwendet als Instrumentalvariable fuer den Anteil der Leibeigenen wieder den Anteil an kloesterlichen Leibeigenen nach der Abschaffung der Leibeigenschaft.

```{r}
ols.model56 = "height ~ lndistMoscowDemeanpost1861 + median_suitDemeanpost1861 + factor(iddistrict) + factor(yearbirth)"
ols.model56 = formula(paste0(ols.model56,'|iddistrict|sh_allserfspost1861~sh_church_1814aver_post1861'))

ols56 = feols(ols.model56, data = subset(data_district, (nobaltics==TRUE) & (moscow != TRUE) & alldrafted>30),cluster = 'iddistref', ssc = ssc(adj = TRUE, cluster.adj = TRUE,fixef.K="none"))

#summary(ols56)

#1-ols56$ssr/ols56$ssr_null
```

#### Table 5 Panel A

In den folgenden Codeabschnitten werden die Ergebnisse der verschiedenen Regressionen zusammengefuehrt.

```{r}
models = list(ols51, ols52, ols53, ols54, ols55, ols56 )

rows = data.frame(
  c('Birth cohort and province or district fixed effects','Province-
specific trends','Observations','R2'),
  c('Yes','No',ols51$nobs,round(r2(ols51, type = "wr2"), digits = 3)),
  c('Yes','No',ols52$nobs,round(r2(ols52, type = "wr2"), digits = 3)),
  c('Yes','No',ols53$nobs,round(r2(ols53, type = "r2"), digits = 3)),
  c('Yes','No',ols54$nobs,round(r2(ols54, type = "wr2"), digits = 3)),
  c('Yes','No',ols55$nobs,round(r2(ols55, type = "wr2"), digits = 3)),
  c('Yes','No',ols56$nobs,round(r2(ols56, type = "r2"), digits = 3))
  )

table5a = modelsummary(models,
             coef_map = c(
'fit_sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation cohorts',             
'sh_serfs_postEmancip' =  'Share of serfs × Post-emancipation cohorts',               
'fit_sh_allserfspost1861' =  'Share of serfs × Post-emancipation cohorts',
'sh_allserfspost1861' =  'Share of serfs × Post-emancipation cohorts',
'lndistMoscowDemeanpost1861' = 'Demeaned log distance to Moscow x Post-emancipation',
'lndist_MoscowDemeanpost1861' = 'Demeaned log distance to Moscow x Post-emancipation',
'median_suitDemeanpost1861' = 'Demeaned crop suitability x Post-emancipation',
'median_suitabilityDemeanpost1861' = 'Demeaned crop suitability x Post-emancipation'),
             add_rows = rows,
             stars = FALSE, 
             gof_map = "none",
fmt = function(x) round(x, digits = 2),
title = 'Table 5—The Abolition of Serfdom and Peasant Living Standards: Draftees’ Height, Panel A. Panel data estimation'

) %>% 
  add_header_above(c(" " = 1, "OLS" = 2, "IV, second
stage" = 1, "OLS" = 2, "IV, second
stage" = 1)) %>% 
  add_header_above(c(" " = 1, "Draftees’ height
(cohorts 1853–1866)
 Province-level data" = 3, "Draftees’ height
(cohorts 1853–1862)
 District-level data" = 3))
  
table5a

#save_kable(table5a,'Tabel 5, Panel A.png')
```

Die Ergebnisse in Tabelle 5 zeigen auch hier einen positiven und signifikanten Zusammenhang zwischen der Groesse von Wehrpflichtigen und der Abschaffung der Leibeigenschaft.

#### Table 5 Panel B

1st-stage-Ergebnisse der IV-Regressionen zusammengetragen in einer Tabelle.

```{r}
models = list("(3)"=ols53$iv_first_stage$sh_serfs_postEmancip,
              "(6)"=ols56$iv_first_stage$sh_allserfspost1861)

rows = data.frame(
  c('Controls as in respective column of panel A'),
  c('Yes'),
  c('Yes'))

attr(rows, 'position') = c(3)

table5b = modelsummary(models,
             coef_map = c(
'sh_nat_monasterialpost1861' =  'Share of nationalized monasterial serfs ×
Post-emancipation cohorts',
'sh_church_1814aver_post1861' =  'Share of nationalized monasterial serfs ×
Post-emancipation cohorts'),
             add_rows = rows,
             stars = FALSE, 
             gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = function(x) format(round(x, 3), big.mark=""))),
fmt = function(x) round(x, digits = 2),
title = 'Table 5—The Abolition of Serfdom and Peasant Living Standards: Draftees’ Height, Panel B. First stages of the corresponding 2SLS panel regressions'

)  %>% 
  add_header_above(c(" " = 1,"Share of serfs × Post-emancipation cohorts
IV, first stage IV" = 1, "Share of serfs × Post-emancipation cohorts
IV, first stage IV" = 1))

  table5b

```

# Fazit

Das Paper legt nahe, dass die Abschaffung der Leibeigenschaft in Russland zu deutlichen Verbesserungen in der Produktivitaet, sowie in der Gesundheit der Landbevoelkerung gekommen ist. Das BIP stieg in Folge der Abschaffung in der zweiten Haelfte des 19. Jahrhundert um ganze 17,7 Prozent. Ein großer Teil dieses Effekts ist den besseren Anreizen der Bauern zu verdanken, die ploetzlich nicht mehr einem Landherren unterstellt waren, sondern fuer sich selbst arbeiteten. Der gesamte Effekt der Abschaffung waere sogar noch wesentlich groesser, haette die Landreform nicht dazu gefuehrt, dass Umverteilungskommunen erstarkten und so einen Teil der Anreize (Anreize fuer Investitionen in Grund und Boden) zunichte machen konnten. 

