Skip to content

Pset on Inference

Download notebook

Goal of the problem set:

  • understand heteroskedasticity and clustered standard errors
  • familiarize yourself with pandas, grouping and formating
  • learn about numba to accelerate your code for large data-sets

some useful references:

We are going to reproduce an exercise similar to the example for the computation of standard error. Start by downloading the CPS data from here. We first load the needed libraries, my solutions wich are hidden to you and the data.

1
2
3
4
5
6
7
8
%cd ..
%load_ext autoreload
%autoreload 2
from solutions.sol_pset2 import *
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from tqdm import tqdm
1
2
3
/Users/thibautlamadon/git
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
1
data = pd.read_stata("~/Downloads/files to share/CPS_2012_micro.dta")

Next generate a fictuous policy that you randomly assigned at the state times gender level. Run the regression and report standard errors given by R for one draw of the poilcy.

1
2
3
4
5
6
7
np.random.seed(60356548) # I fix the seed to make sure the draws are reproducible

# we draw the policy for each state
fpol = { k:np.random.uniform() > 0.5 for k in np.unique(data['statefip']) }
data['fp'] = data['statefip'].map(fpol)

data
year statefip wtsupp age sex yrseduc wage_per_hour lnwage age2 fp
0 2012 Maine 569.43 44 Female 14.0 7.020000 1.948763 1936.0 True
1 2012 Maine 595.47 25 Male 16.0 2.117143 0.750067 625.0 True
2 2012 Maine 635.66 61 Female 16.0 16.672501 2.813761 3721.0 True
3 2012 Maine 635.66 62 Male 16.0 17.784000 2.878299 3844.0 True
4 2012 Maine 513.39 25 Female 12.0 9.633000 2.265195 625.0 True
... ... ... ... ... ... ... ... ... ... ...
65680 2012 Hawaii 282.87 24 Male 16.0 28.899000 3.363807 576.0 False
65681 2012 Hawaii 362.23 49 Male 18.0 18.525000 2.919121 2401.0 False
65682 2012 Hawaii 362.23 44 Female 20.0 18.525000 2.919121 1936.0 False
65683 2012 Hawaii 282.87 28 Male 18.0 24.082500 3.181485 784.0 False
65684 2012 Hawaii 282.87 28 Female 18.0 11.856000 2.472834 784.0 False

65685 rows × 10 columns

Question 1

As an exercise to get acquainted with pandas, here I would like for you to implement the following procedure:

  1. compute the average wage in each state for each age decile (less than 20, betwen 20 and 30, ...)
  2. pivot the table with states in rows and age deciles in columns
  3. using pandas styling, I want you to highlight the lowest wage in each age decile and I want you to format the average wage with onely 2 decimal points.
1
2
# we tabulate the counts by age group and state
tibo_question1(data)
age_grp 20.0 30.0 40.0 50.0 60.0 70.0
statefip
Alabama 8.20 11.33 15.50 19.39 17.38 15.21
Alaska 8.68 14.14 19.92 19.59 20.67 20.46
Arizona 9.76 11.97 16.55 18.63 18.46 21.30
Arkansas 6.52 12.47 14.95 17.36 16.13 14.26
California 8.52 12.82 18.98 20.55 21.52 21.34
Colorado 9.63 12.35 19.19 21.00 21.80 21.29
Connecticut 6.79 13.29 20.98 24.64 25.11 20.69
Delaware 6.05 12.61 17.07 19.70 17.89 20.10
District of Columbia 10.24 17.36 26.19 25.80 25.35 29.54
Florida 7.42 11.93 16.08 18.60 18.95 17.42
Georgia 7.84 12.52 17.03 19.88 19.21 15.75
Hawaii 10.73 11.90 17.18 18.53 18.68 19.63
Idaho 6.17 11.31 16.17 17.45 17.50 13.56
Illinois 6.32 12.82 18.08 20.63 21.26 20.45
Indiana 6.91 12.84 16.55 18.36 18.28 19.17
Iowa 6.82 11.93 16.47 17.85 16.86 18.02
Kansas 7.09 11.30 17.16 17.11 19.50 18.55
Kentucky 8.08 11.88 15.83 16.51 17.61 17.01
Louisiana 6.94 12.73 14.69 16.83 15.77 23.54
Maine 5.33 12.13 16.72 17.50 18.44 15.45
Maryland 7.67 15.10 20.85 22.72 22.84 24.99
Massachusetts 7.83 14.09 22.51 23.45 23.85 22.04
Michigan 6.53 12.10 18.30 19.88 20.40 20.24
Minnesota 6.71 13.07 19.61 20.67 19.35 16.23
Mississippi 6.10 11.67 14.92 16.70 14.97 19.39
Missouri 6.36 11.81 16.48 18.48 19.46 17.78
Montana 5.02 11.13 14.47 16.56 15.90 16.35
Nebraska 6.97 11.62 16.71 16.99 16.65 16.14
Nevada 9.01 12.16 16.87 17.99 17.47 17.67
New Hampshire 8.62 14.04 20.53 20.96 21.03 21.27
New Jersey 8.41 12.60 21.93 22.38 24.25 20.92
New Mexico 6.88 11.87 17.32 20.24 18.94 22.81
New York 7.24 14.76 19.35 20.09 20.53 19.21
North Carolina 7.46 11.89 15.89 18.26 19.69 18.44
North Dakota 6.80 12.60 16.92 17.70 16.56 14.66
Ohio 6.69 11.11 16.18 17.74 19.00 20.15
Oklahoma 6.06 10.85 15.90 17.77 18.41 20.17
Oregon 6.90 12.05 17.71 20.14 19.64 19.14
Pennsylvania 8.93 12.74 18.96 20.05 19.13 15.29
Rhode Island 8.02 12.71 18.33 21.42 22.38 20.13
South Carolina 5.23 10.63 13.82 15.68 15.67 16.60
South Dakota 5.38 10.99 14.96 15.74 16.01 15.31
Tennessee 5.02 11.40 16.19 16.49 15.87 17.85
Texas 7.62 12.44 16.57 16.95 17.83 18.53
Utah 9.07 13.14 16.75 20.00 20.05 19.97
Vermont 8.23 12.71 16.30 17.24 18.89 18.57
Virginia 7.01 14.41 19.78 22.61 22.40 20.67
Washington 6.62 13.34 19.38 19.99 20.53 22.46
West Virginia 6.97 11.81 15.80 16.56 16.36 21.27
Wisconsin 7.12 11.72 16.97 18.75 19.11 17.61
Wyoming 7.19 13.69 18.80 18.08 17.69 15.72

A first regression

1
2
3
4
5
import statsmodels.api as sm
import statsmodels.formula.api as smf

results = smf.ols('lnwage ~ fp', data=data).fit()
print(results.summary())
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 lnwage   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                     11.53
Date:                Mon, 27 Apr 2020   Prob (F-statistic):           0.000687
Time:                        13:45:10   Log-Likelihood:                -62336.
No. Observations:               65685   AIC:                         1.247e+05
Df Residuals:                   65683   BIC:                         1.247e+05
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      2.6687      0.003    836.168      0.000       2.662       2.675
fp[T.True]     0.0168      0.005      3.395      0.001       0.007       0.026
==============================================================================
Omnibus:                        5.904   Durbin-Watson:                   1.582
Prob(Omnibus):                  0.052   Jarque-Bera (JB):                5.915
Skew:                           0.023   Prob(JB):                       0.0519
Kurtosis:                       2.992   Cond. No.                         2.47
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Note We do not control for state specific fixed effect as these would would be perfectly colinear with the policy.

Now this is surprising. We generated fp randomly across states, yet we find a significant effect of the policy on log wages. To gain understanding on what is happening we will generate our own data in a way where we control exactly what is happening.

IID errors

Let's start by reassuring ourselves. Let's use an IID data generating process (DGP), run the regression and check the significance.

  1. compute the variance of lnwage in the sample. This is an estimate of our homoskedastic error.
  2. simulate a fictuous outcome y2 simply equal to normal error with the estimated variance, and truly independent across individuals. Use .assign(y2 = ...) on the pandas dataFrame.
  3. regress this outcome y2 on fp, our fictuous policy and collect the coefficient, also save if the coefficient is significant at 5%.
  4. run steps (2,3) 500 times.

Question 2

Follow the previous steps and report how often we would consider teh coefficient on fpsignificant at 5%. You should find something close to 5% and you should feel better about the theory!

1
2
3
4
5
# I colled tstats across 500 replications and plot the histogram
tstats = tibo_question2(data)
plt.hist(tstats)
plt.axvline(1.96,color='red',linestyle=":")
plt.axvline(-1.96,color='red',linestyle=":")
1
2
3
4
5
6
7
100%|██████████| 500/500 [02:51<00:00,  2.91it/s]





<matplotlib.lines.Line2D at 0x1a1c340ad0>

png

1
2
# I check that the distribution is outside of the 0.025 to 0.0975 quantiles of the normal with probability close to 0.05
(np.abs(tstats)>1.96).mean()
1
0.038

Back to using wages in data

We go back and we do the same exercise, only this time we redraw the policy in each replciations and regress the wages on it, and collect the t-stats. We then plot the histogram of the t-stats and compute the probability of rejecting the null that the coefficient on the policy is 0.

Question 3

Replicate 500 times the very first regression to recreate an histogram of t-stats similar to the following one. Compute the probability to reject the null of effect and comment on its value given the way you simulated the data (2 sentences max).

1
2
3
4
tstats = tibo_question3(data)
plt.hist(tstats)
plt.axvline(1.96,color='red',linestyle=":")
plt.axvline(-1.96,color='red',linestyle=":")
1
2
3
4
5
6
7
100%|██████████| 500/500 [03:21<00:00,  2.48it/s]





<matplotlib.lines.Line2D at 0x1a1c82e810>

png

Heteroskedastic errors

Now we want to simulate heteroskedastic robust standard errors which requires us to use some co-variates. We then want to repeat the previous procedure, but we are going to use a different code for simulation and a new test for the significance. statsmodels can do that for you using the cov_type argument of the fit method.

We want to check this by simulating from a model with heteroskedesatic errors. To do so we are going to use linear model for the variance:

y_{i} = 0 + s(x_i) \cdot \epsilon_i

We are going to use a linear specification for s^2(x_i).

  1. run the following regression lnwage ~ yrseduc + age + age^2 to get the m(x_i) = x_i' \beta where the x_i are education, age and age squared.
  2. extract the residuals \hat{u}_i from the previous regression. Then regress \hat{u}^2_i on yrseduc + age + age^2, to get the s^2(x_i) = x_i' \gamma model.
  3. using s^2(x_i) = x_i' \gamma, construct $ y_{i} = 0 + s(x_i) \cdot \epsilon_i$ where \epsilon_i is drawn iid Normal(0,1).
  4. replicate this 500 times, evaluate the significance of fp using heteroskedastic roduct inference by calling .fit(cov_type='HC0'), also save without using the robust errors.

Question 4

Follow the steps and report the distributions of t-stats, as well as the rejection rates for each of the two variance specifications.

1
tstats = tibo_question4(data,R=500)
1
100%|██████████| 500/500 [05:28<00:00,  1.52it/s]
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
plt.subplot(1,2,1)
plt.hist(tstats[:,0])
plt.axvline(1.96,color='red',linestyle=":")
plt.axvline(-1.96,color='red',linestyle=":")
plt.title('Homoskedastic')

plt.subplot(1,2,2)
plt.hist(tstats[:,1])
plt.axvline(1.96,color='red',linestyle=":")
plt.axvline(-1.96,color='red',linestyle=":")
plt.title('Heteroskedastic')

plt.show()

png

State clustered errors

We are going to simulate corrolated error within state. To do so we continue to draw independent error terms u_i but we add to them a common draw at the state level. In practice you can do that by doing similar to drawing the fp. We then create the outcome as

y_{ig} = (1-\rho) u_i + \rho v_g

Here is a way of doing this:

1
2
3
rho=0.5
data['u'] = np.random.normal(size=len(data))
data['v'] = data.groupby('statefip')['u'].transform(lambda x : (1-rho)*x.to_numpy() + rho*np.random.normal(size=1))

Question 5

First, explain the groupby expression. Next find another way of doing line 2 and 3, but try to minize the number of characters (see python golfing or more generaly wikipedia).

Before we get started estimating the effect of the policy, let's write a function that computes the within state correlation in the error term.

To that end draw 500 pairs of people where the 2 people in a pair belong to the same state. This gives you 2 vector of length 500. Compute the correlation between these 2 vectors. Use the sample and query methods of pandas. Similarly construct a placebo function that just takes random pairs of workers without imposing that they come from the same state.

Question 6

Write your own version of the tibo_pair_cor and tibo_cor_placebo functions, see what my version delivers:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
rho=0.5
data['u'] = np.random.normal(size=len(data))
data['v'] = data.groupby('statefip')['u'].transform(lambda x : (1-rho)*x.to_numpy() + rho*np.random.normal(size=1))

print(tibo_pair_cor(data))
print(tibo_pair_cor_placebo(data))

rho=0.8
data['u'] = np.random.normal(size=len(data))
data['v'] = data.groupby('statefip')['u'].transform(lambda x : (1-rho)*x.to_numpy() + rho*np.random.normal(size=1))

print(tibo_pair_cor(data))
print(tibo_pair_cor_placebo(data))
1
2
3
4
0.536557231009759
-0.010315753596236765
0.9487955581792823
0.01587768220072409

We are then going to replicate the data construction 500 times.

Question 7

For \rho=0.2,0.5,0.8 run 500 replications and report the proportion at each value of \rho for which the coefficient on\text{fp} is significant at 5%. Report the results for three different regression using regular fit, using cov_type='HC0' and using cov_type='cluster'.

1
2
3
4
5
6
7
8
# here is an example of clusetered errors 
rho=0.8
data['u'] = np.random.normal(size=len(data))
data['v'] = data.groupby('statefip')['u'].transform(lambda x : (1-rho)*x.to_numpy() + rho*np.random.normal(size=1))

data['y'] =  data.v
results = smf.ols('y ~ fp', data=data).fit(cov_type='cluster',cov_kwds={"groups":data['statefip']})
print(results.summary())
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.027
Model:                            OLS   Adj. R-squared:                  0.027
Method:                 Least Squares   F-statistic:                     1.337
Date:                Sun, 26 Apr 2020   Prob (F-statistic):              0.253
Time:                        18:25:20   Log-Likelihood:                -73718.
No. Observations:               65685   AIC:                         1.474e+05
Df Residuals:                   65683   BIC:                         1.475e+05
Df Model:                           1                                         
Covariance Type:              cluster                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.1198      0.142     -0.844      0.399      -0.398       0.158
fp[T.True]     0.2495      0.216      1.156      0.248      -0.173       0.672
==============================================================================
Omnibus:                      160.476   Durbin-Watson:                   0.145
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              152.337
Skew:                           0.094   Prob(JB):                     8.33e-34
Kurtosis:                       2.859   Cond. No.                         2.47
==============================================================================

Warnings:
[1] Standard Errors are robust to cluster correlation (cluster)
1
2
3
4
5
6
7
8
# and here the regression without clustered standard errors
rho=0.8
data['u'] = np.random.normal(size=len(data))
data['v'] = data.groupby('statefip')['u'].transform(lambda x : (1-rho)*x.to_numpy() + rho*np.random.normal(size=1))

data['y'] =  data.v
results = smf.ols('y ~ fp', data=data).fit()
print(results.summary())
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                  0.001
Method:                 Least Squares   F-statistic:                     80.51
Date:                Sun, 26 Apr 2020   Prob (F-statistic):           2.96e-19
Time:                        18:25:50   Log-Likelihood:                -69687.
No. Observations:               65685   AIC:                         1.394e+05
Df Residuals:                   65683   BIC:                         1.394e+05
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.2244      0.004     62.874      0.000       0.217       0.231
fp[T.True]    -0.0497      0.006     -8.973      0.000      -0.061      -0.039
==============================================================================
Omnibus:                     1332.865   Durbin-Watson:                   0.167
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1411.728
Skew:                           0.357   Prob(JB):                    2.80e-307
Kurtosis:                       2.919   Cond. No.                         2.47
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
1