# Leave-one-out mean removal stringency histogram

I am trying to construct the kind of histogram described here (https://blogs.worldbank.org/impactevaluations/judge-leniency-iv-designs-now-not-just-crime-studies). I don't really feel like my code is accomplishing this because my output seems off. Is what I have yielding the correct histogram?


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from sklearn.utils import resample

# 2. Calculate the residuals after accounting for Day * Hospital fixed effects
X = pd.get_dummies(df[['weekday', 'hospital_id']], drop_first=True)
model = sm.OLS(df['removal_stringency_instrument'], X).fit()
df['residuals'] = model.resid

# 3. Plot the histogram of residuals
plt.figure(figsize=(10, 6))
sns.histplot(df['residuals'], kde=False, bins=30, color='lightblue')

# 4. Calculate and plot non-parametric regression using bootstrap for 95% CI
def bootstrap_ci(data, n_iterations=1000):
"""Function to generate bootstrap CI for non-parametric regression."""
sample_size = len(data)
bootstrapped_means = np.zeros(n_iterations)

for i in range(n_iterations):
sample = resample(data)
bootstrapped_means[i] = np.mean(sample['TREATMENT'])

lower = np.percentile(bootstrapped_means, 2.5)
upper = np.percentile(bootstrapped_means, 97.5)

return lower, upper

# Grouping by residuals and getting means of treatment
grouped_means = df.groupby('residuals').agg({'TREATMENT': 'mean'}).reset_index()
sns.lineplot(data=grouped_means, x='residuals', y='TREATMENT', color='black')

# Highlighting the 95% CI using bootstrap
lower, upper = bootstrap_ci(df[['residuals', 'TREATMENT']])
plt.axhline(y=lower, color='red', linestyle='--')
plt.axhline(y=upper, color='red', linestyle='--')

plt.title("Histogram with Non-parametric Regression and 95% CI")
plt.show()



Thanks.