4) Do Protestants Prefer Less Leisure than Catholics?¶
E-mail: econometrics.methods@gmail.com
Last updated: 10-4-2020
Max Weber (1930) argued that the Protestant ethic, especially the variants of Calvinism are more aligned with capitalism than the Catholicism. Weber (1930) observed that Protestant regions in Northern Europe were more developed than the Catholic regions in the Southern Europe. He hypothesized that Protestants work harder, save more, rely more on themselves, and expect less from the Government. All characteristics that would lead to a greater economic prosperity.
Maybe, it is not the religion the cause of a great economic performance. Education is a confound factor. Historically, Protestants have higher level of literacy, because they were incentivized to read the Bible.
The causal effect can be reverse as well. Perhaps, an industrial person is more likely to become Protestant. Religion is a choice variable. People self-select the ideology that confirms their own view of the world.
Let’s open the data from Basten & Betz (2013). Each row represents a municipality in Western Switzerland, the cantons of Vaud and Fribourg.
# Load data from Basten & Betz (2013)
import numpy as np
import pandas as pd
path = "https://github.com/causal-methods/Data/raw/master/"
df = pd.read_stata(path + "finaldata.dta")
df.head(4)
| id | name | district | district_name | foreignpop2000 | prot1980s | cath1980s | noreligion1980s | canton | totalpop2000 | ... | pfl | pfi | pfr | reineink_pc_mean | meanpart | popden2000 | foreignpopshare2000 | sub_hs | super_hs | murten | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2135.0 | Gruyères | 1003.0 | District de la Gruyère | 159 | 0.062251 | 0.907336 | 1.621622 | 10 | 1546.0 | ... | 42.722500 | 52.234196 | 40.143444 | 48.099865 | 40.359428 | 54.455795 | 10.284605 | 85.910339 | 6.770357 | 0.0 |
| 1 | 2128.0 | Châtel-sur-Montsalvens | 1003.0 | District de la Gruyère | 23 | 0.053191 | 0.917526 | 2.061856 | 10 | 205.0 | ... | 49.223751 | 56.793213 | 44.365696 | 42.465569 | 45.434593 | 102.499985 | 11.219512 | 83.832336 | 8.383233 | 0.0 |
| 2 | 2127.0 | Charmey | 1003.0 | District de la Gruyère | 166 | 0.028424 | 0.960818 | 0.255537 | 10 | 1574.0 | ... | 41.087502 | 53.120682 | 39.674942 | 44.451229 | 42.641624 | 20.066286 | 10.546379 | 87.331535 | 6.019766 | 0.0 |
| 3 | 2125.0 | Bulle | 1003.0 | District de la Gruyère | 2863 | 0.053967 | 0.923239 | 1.013825 | 10 | 11149.0 | ... | 47.326248 | 55.033939 | 43.350178 | 50.217991 | 40.885822 | 897.664978 | 25.679434 | 82.203598 | 8.904452 | 0.0 |
4 rows × 157 columns
Switzerland is a very diverse country in terms of geography and institutions. It is not fair to compare a rural Catholic that lives in the Alpes with an urban high educated Protestant in Zurich.
Historically, the cities had different incentives to adopt Protestantism or remain Catholic. Cities with a stronger merchant guild were more likely to adopt the Protestantism; whereas cities governed by aristocrats were more likely to remain Catholic.
There are too much confound factors, if we use the whole country. The analysis will be restricted to the cantons of Vaud (historically Protestant) and Fribourg (historically Catholic). See the map below from Basten & Betz (2013).
This region of 4,883 \(km^2\), that represents 4.5 percent of Switzerland, is institutionally and geographically homogeneous. In 1536, the canton of Vaud didn’t self-select to become Protestant, but it was forced because of a war. Therefore, this is a quasi-experiment setting, where treatment region and control region are similar to each other, because of a historical event.
Source: Basten & Betz (2013)
In the graphic below, we can see that higher the proportion of protestants in a municipality, lower the preference for leisure. The blue dots are historically Catholic municipalities (Fribourg), while the red dots are historically Protestant municipalities (Vaud). It looks like distinct subgroups. How can we figure out if there is evidence of causal effect or it is a mere correlation?
# Create variable "Region" for the graphic
def category(var):
if var == 1:
return "Protestant (Vaud)"
else:
return "Catholic (Fribourg)"
df['Region'] = df["vaud"].apply(category)
# Rename variables with auto-explanatory names
df = df.rename(columns={"prot1980s": "Share_of_Protestants",
"pfl": "Preference_for_Leisure"})
# Scatter plot
import plotly.express as px
leisure = px.scatter(df,
x="Share_of_Protestants",
y="Preference_for_Leisure",
color="Region")
leisure.show()
Let’s refine the analysis. Remember that Regression Discontinuity is the closer technique to an experiment.
In the graphic below, there is a discontinuity in the preference for leisure at border distance = 0. The border distance above 0 encompasses historically Protestant municipalities (Vaud); whereas, the border distance below 0 encompasses historically Catholic municipalities (Fribourg).
The running variable “Border Distance” determines the region, but not the share of protestants. However, the share of protestants increases as function of the distance.
df = df.rename(columns={"borderdis": "Border_Distance_in_Km"})
leisure = px.scatter(df,
x="Border_Distance_in_Km",
y="Preference_for_Leisure",
color="Share_of_Protestants",
title="Discontinuity at Distance = 0")
leisure.show()
As the border is arbitrary, that is, determined by a historic event, the municipalities closer to the border are likely to be more similar to each other than the municipalities far away of the border. Therefore, let’s restrict the analysis to municipalities inside a range of 5 Km.
df5 = df[df['Border_Distance_in_Km'] >= -5]
df5 = df5[df5['Border_Distance_in_Km'] <= 5]
The simple mean comparison shows that the preference for leisure is lower in the Protestant municipalities (39.5%) compared with the Catholic municipalities (48.2%). The difference is -8.7%.
Note that the Protestant region has higher mean income measure in Swiss Franc (47.2K vs 43.7K) and higher inequality captured by Gini index (0.36 vs 0.30).
df5 = df5.rename(columns={"reineink_pc_mean": "Mean_Income_(CHF)",
"Ecoplan_gini" : "Gini_1996"})
outcome = ['Preference_for_Leisure', 'Mean_Income_(CHF)', 'Gini_1996']
The comparison above is not “bad”, considering that only the municipalities inside a 5 Km range are used (49 Catholic and 84 Protestant municipalities). Furthermore, the two regions are similar in terms of share of no religious affiliation in 1980 (1.7% vs 2.9%) and altitude above the sea level (642 vs 639 meters).
However, a more credible approach is to use a regression discontinuity framework with the running variable (\(r_i\)).
control = ['noreligion1980s', 'altitude']
df5.loc[:, control].groupby(df5['Region']).agg([np.mean, np.std]).T
| Region | Catholic (Fribourg) | Protestant (Vaud) | |
|---|---|---|---|
| noreligion1980s | mean | 1.729423 | 2.949534 |
| std | 1.499346 | 2.726086 | |
| altitude | mean | 642.591858 | 639.607117 |
| std | 120.230320 | 113.563847 |
In the graphic “Fuzzy Regression Discontinuity”, we can clearly see that the running variable “Border Distance” is very correlated with the treatment variable “Share of Protestants”. The variable “Border Distance” does not determine the treatment status but increases the probability of being Protestant. Therefore, this is a case of a Fuzzy Regression Discontinuity and not a Sharp Regression Discontinuity.
Let \(D_i\) be the treatment status of unit \(i\). \(P(D_i=1|r_i)\) is a jump in the probability of treatment at cutoff \(r_0\):
where \(f_1(r_i)\) and \(f_0(r_i)\) are functions that can assume any value. In the Sharp Regression Discontinuity, \(f_1(r_i)\) was 1 and \(f_0(r_i)\) was 0.
fuzzy = px.scatter(df5,
x="Border_Distance_in_Km",
y="Share_of_Protestants",
color="Region",
title="Fuzzy Regression Discontinuity")
fuzzy.show()
In the graphic below, the variable share of protestants is simulated to illustrate what would be a case of Sharp Regression Discontinuity.
def dummy(var):
if var >= 0:
return 1
else:
return 0
df5["Simulated_Share_Protestant"] = df5["Border_Distance_in_Km"].apply(dummy)
sharp = px.scatter(df5,
x="Border_Distance_in_Km",
y="Simulated_Share_Protestant",
color="Region",
title="Sharp Regression Discontinuity")
sharp.show()
Let’s assume \(Y\) = Preference_for_Leisure, \(D_r\) = Share_of_Protestants, and \(r\) = Border_Distance. Now, we have a problem to estimate the equation below:
The variable of interest \(D_r\) is not anymore “purified” by \(r\), that is, it is not anymore completely determined by the running variable \(r\). Therefore, \(D_r\) is likely to be correlated with the error term \(\epsilon\):
We can fix this problem using an instrumental variable \(Z\) that is uncorrelated with the error term \(\epsilon\), and correlated with \(D_r\) after controlling for other factors:
The natural candidate for \(Z\) is the variable “vaud”: 1 if it is a municipality in Vaud; and 0 if it is a municipality in Fribourg. There is no reason to believe that this variable is correlated with the error term \(\epsilon\), as the border that divides the region was determined in 1536, when the Republic of Berne conquered Vaud. The second assumption is also valid, as more Protestants live in the Vaud region than Fribourg.
The instrumental variable method consists first in “purifying” \(D_r\) using \(Z\):
Then, we get the fitted values of \(\hat{D}_r\) by running an ordinary least square (OLS) and plug it in the equation:
The logic is that the “purified” \(\hat{D}_r\) is uncorrelated with the error term \(\epsilon\). Now, we can run an ordinary least square (OLS) to get the isolated effect of \(\hat{D}_r\) on \(Y\), that is, \(\rho\) will be the unbiased causal effect.
# Install libray to run Instrumental Variable estimation
!pip install linearmodels
Requirement already satisfied: linearmodels in c:\anaconda\envs\textbook\lib\site-packages (4.17)
Requirement already satisfied: scipy>=1 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.5.3)
Requirement already satisfied: pandas>=0.23 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.1.3)
Requirement already satisfied: property-cached>=1.6.3 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.6.4)
Requirement already satisfied: patsy in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.5.1)
Requirement already satisfied: Cython>=0.29.14 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.29.21)
Requirement already satisfied: statsmodels>=0.9 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (0.10.2)
Requirement already satisfied: numpy>=1.15 in c:\anaconda\envs\textbook\lib\site-packages (from linearmodels) (1.19.2)
WARNING: No metadata found in c:\anaconda\envs\textbook\lib\site-packages
ERROR: Could not install packages due to an EnvironmentError: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA'
The computer automatically run the two stages of Instrumental Variable (IV) procedure. We indicated that the endogenous variable \(D_r\) is “Share of Protestants”, and the instrument variable \(Z\) is “vaud”. We also add the variable “t_dist” that is the interaction between the variables “vaud” and “Border Distance”.
The result is that “Share of Protestants” decreases the preference for leisure in 13.4%.
from linearmodels.iv import IV2SLS
iv = 'Preference_for_Leisure ~ 1 + Border_Distance_in_Km + t_dist + [Share_of_Protestants ~ vaud]'
iv_result = IV2SLS.from_formula(iv, df5).fit(cov_type='robust')
print(iv_result)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:
pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
IV-2SLS Estimation Summary
==================================================================================
Dep. Variable: Preference_for_Leisure R-squared: 0.4706
Estimator: IV-2SLS Adj. R-squared: 0.4583
No. Observations: 133 F-statistic: 99.425
Date: Fri, Nov 06 2020 P-value (F-stat) 0.0000
Time: 21:12:05 Distribution: chi2(3)
Cov. Estimator: robust
Parameter Estimates
=========================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
-----------------------------------------------------------------------------------------
Intercept 50.528 1.8885 26.755 0.0000 46.826 54.229
Border_Distance_in_Km 0.4380 0.6291 0.6963 0.4862 -0.7949 1.6710
t_dist -0.3636 0.7989 -0.4552 0.6490 -1.9294 1.2022
Share_of_Protestants -13.460 3.1132 -4.3235 0.0000 -19.562 -7.3582
=========================================================================================
Endogenous: Share_of_Protestants
Instruments: vaud
Robust Covariance (Heteroskedastic)
Debiased: False
C:\Anaconda\envs\textbook\lib\site-packages\linearmodels\iv\data.py:25: FutureWarning:
is_categorical is deprecated and will be removed in a future version. Use is_categorical_dtype instead
We can also check the first stage to see if the instrumental variable “vaud” is correlated with “Share of Protestants” after controlling for other factors like “Border Distance” and “t_dist”. Vaud increases 67% the share of Protestants. The t-value of “vaud” is 20, that is, statistically significant without any margin of doubt.
Therefore, we are confident that the second stage result is more credible than the simple mean comparison. The Instrumental Variable impact of 13.4% is more credible than the simple mean comparison of 8.7%.
print(iv_result.first_stage)
First Stage Estimation Results
===============================================
Share_of_Protestants
-----------------------------------------------
R-squared 0.9338
Partial R-squared 0.7095
Shea's R-squared 0.7095
Partial F-statistic 403.04
P-value (Partial F-stat) 0.0000
Partial F-stat Distn chi2(1)
========================== ===========
Intercept 0.1336
(7.7957)
Border_Distance_in_Km 0.0169
(2.8397)
t_dist -0.0057
(-0.4691)
vaud 0.6709
(20.076)
-----------------------------------------------
T-stats reported in parentheses
T-stats use same covariance type as original model
The simple mean comparison result of 8.7% is closer to the result of 9% from the naive Sharp Regression Discontinuity (SRD) below. The Vaud region has a 9% less preference for leisure than the Fribourg. We cannot conclude that Protestants have a 9% less preference for leisure than Catholics. The Vaud region is not 100% Protestant. Neither the Fribourg region is 100% Catholic.
The Fuzz Regression Discontinuity (FRD), that uses the Instrumental Variable (IV) estimation, is a correction for the naive comparison. The FRD isolates the impact of Protestants on preference for leisure. Therefore, the most credible estimation is that Protestants have 13.4% less preference for leisure than Catholics.
naive_srd = 'Preference_for_Leisure ~ 1 + vaud + Border_Distance_in_Km + t_dist'
srd = IV2SLS.from_formula(naive_srd, df5).fit(cov_type='robust')
print(srd)
OLS Estimation Summary
==================================================================================
Dep. Variable: Preference_for_Leisure R-squared: 0.3830
Estimator: OLS Adj. R-squared: 0.3686
No. Observations: 133 F-statistic: 91.767
Date: Fri, Nov 06 2020 P-value (F-stat) 0.0000
Time: 21:12:05 Distribution: chi2(3)
Cov. Estimator: robust
Parameter Estimates
=========================================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
-----------------------------------------------------------------------------------------
Intercept 48.730 1.5731 30.976 0.0000 45.646 51.813
vaud -9.0303 2.1821 -4.1383 0.0000 -13.307 -4.7534
Border_Distance_in_Km 0.2107 0.6028 0.3496 0.7266 -0.9707 1.3922
t_dist -0.2874 0.8479 -0.3390 0.7346 -1.9493 1.3744
=========================================================================================
Exercises¶
1| What would be the confound factors in estimating the causal impact of Protestants against Catholics on economic prosperity of Switzerland (whole country)? Explain in detail each confound factor.
2| Somebody could argue that Western Switzerland was a diverse region before 1536, when the Republic of Berne conquered Vaud. In this case, Fribourg (Catholic) would not be a reasonable control group for Vaud (Protestant). What variables could be used to test the homogeneity/diversity of the region before 1536? Indicate if the data exists, or if it is feasible to collect the data.
3| I replicated the main results of the paper Basten and Betz (2013) and I noticed that they don’t control for education. Becker & Woessmann(2009) argue that Protestantism effect on the economic prosperity of Prussia in nineteenth century is due to higher literacy. Basten and Betz (2013) present the table below in the Online Appendix. Do the numbers in the table strengthen or weaken the results of Basten and Betz (2013)? Justify your reasoning.
Source: Basten & Betz (2013)
4| Preference for leisure is a self-reported variable in a survey. Maybe, Protestants have more incentives to declare lower preference for leisure than Catholics. In the real word, Protestants might enjoy leisure as much as Catholics. Declared preference might not match with real behavior. The relevant reseach question is if the religion (Protestant) causes people to be hard-working in actuality. Be creative and propose a way (methods and data) to fix the problem described above.
5| Use the data from Basten and Betz (2013) to investigate if Protestants cause higher mean income in Western Switzerland. Adopt the specifications that you think it is the most credible to recover the unbiased causal effect. Explain and justify each step of your reasoning. Interpret the main results.
Reference¶
Basten, Christoph, and Frank Betz (2013). Beyond Work Ethic: Religion, Individual, and Political Preferences. American Economic Journal: Economic Policy, 5 (3): 67-91. Online Appendix
Becker, Sascha O., and Ludger Woessmann. (2009). Was Weber Wrong? A Human Capital Theory of Protestant Economic History. Quarterly Journal of Economics 124 (2): 531–96.
Weber, Max. (1930). The Protestant Ethic and the Spirit of Capitalism. New York: Scribner, (Orig.pub. 1905).