Download the Jupyter Notebook for this section: statsmodels.ipynb

# statsmodels¶

:

import pyhdfe
import numpy as np
import statsmodels.api as sm
from sklearn import datasets

pyhdfe.__version__

:

'0.1.0'


In this tutorial, we’ll use the boston data set from scikit-learn to demonstrate how pyhdfe can be used to absorb fixed effects before running regressions with statsmodels. We’ll also demonstrate how pyhdfe can be used to compute degrees of freedom used by fixed effects.

First, load the data set and create a matrix of fixed effect IDs. We’ll use a dummy for the Charles river and an index of accessibility to radial highways.

:

boston = datasets.load_boston().data
ids = boston[:, [3, 8]]
ids

:

array([[0., 1.],
[0., 2.],
[0., 2.],
...,
[0., 1.],
[0., 1.],
[0., 1.]])


Next, choose our variables: per capita crime rate, proportion of residential land zoned for lots over 25,000 square feet, and proportion of non-retail business acres per town.

:

variables = boston[:, :3]
variables

:

array([[6.3200e-03, 1.8000e+01, 2.3100e+00],
[2.7310e-02, 0.0000e+00, 7.0700e+00],
[2.7290e-02, 0.0000e+00, 7.0700e+00],
...,
[6.0760e-02, 0.0000e+00, 1.1930e+01],
[1.0959e-01, 0.0000e+00, 1.1930e+01],
[4.7410e-02, 0.0000e+00, 1.1930e+01]])


The create function initializes an Algorithm for fixed effect absorption that can residualize matrices with Algorithm.residualize. We’ll use the default algorithm. You may want to try other algorithms if it takes a long time to absorb fixed effects into your data.

:

algorithm = pyhdfe.create(ids)
residualized = algorithm.residualize(variables)
residualized

:

array([[-1.08723516e-01, -2.20167195e+01, -2.65583593e+00],
[-5.59754167e-02, -2.04166667e+01, -2.56083333e+00],
[-5.59954167e-02, -2.04166667e+01, -2.56083333e+00],
...,
[-5.42835164e-02, -4.00167195e+01,  6.96416407e+00],
[-5.45351644e-03, -4.00167195e+01,  6.96416407e+00],
[-6.76335164e-02, -4.00167195e+01,  6.96416407e+00]])


We can now run a regression of per capita crime rate on the other two variables and our fixed effects.

:

y = residualized[:, ]
X = residualized[:, 1:]
ols = sm.OLS(y, X)
result = ols.fit()
result.params

:

array([-6.97058632e-05,  5.53038164e-02])


Standard errors can be adjusted to account for the degrees of freedom that are lost because of the fixed effects. By default, fixed effect degrees of freedom are computed when create initializes an algorithm and are stored in Algorithm.degrees.

:

se = result.HC0_se
se

:

array([0.00109298, 0.00962226])

:

se_adjusted = np.sqrt(np.square(se) * result.df_resid / (result.df_resid - algorithm.degrees))

:

array([0.00110398, 0.00971916])