Skip to content

Instrumental Variables: Two Stage Least Squares in Python

Last Update: March 24, 2022

Instrumental Variables: Two Stage Least Squares in Python can be done using linearmodels package IV2SLS function found within linearmodels.iv.model module for estimating linear regression with independent variables which are correlated with error term (endogenous). Main parameters within IV2SLS function are dependent with model dependent variable, exog with model exogenous independent variable, endog with model endogenous independent variable and instruments with model instrumental variables.

As example, we can compare estimated coefficients tables and F-statistics from original multiple linear regression of house price explained by its lot size and number of bedrooms and second stage least squares multiple linear regression of house price explained by its lot size first stage multiple linear regression fitted values and number of bedrooms with whether house has a driveway and number of garage places as instrumental variables using data included within AER R package HousePrices object [1].

First, we import packages statsmodels for data downloading and ordinary least squares original model fitting and linearmodels for two stage least squares model fitting [2].

In [1]:
import statsmodels.api as sm
import statsmodels.formula.api as smf
import linearmodels.iv.model as lm

Second, we create houseprices data object using get_rdataset function and display first five rows and first three columns together with sixth and eleventh columns of data using print function and head data frame method to view its structure.

In [2]:
houseprices = sm.datasets.get_rdataset(dataname="HousePrices", package="AER", cache=True).data
print(houseprices.iloc[:, list(range(3)) + [5] + [10]].head())
Out [2]:
     price  lotsize  bedrooms driveway  garage
0  42000.0     5850         3      yes       1
1  38500.0     4000         2      yes       0
2  49500.0     3060         3      yes       0
3  60500.0     6650         3      yes       0
4  61000.0     6360         2      yes       0

Third, we fit original model with ols function using variables within houseprices data object and store outcome within mlr1 object. Within ols function, parameter formula = “price ~ lotsize + bedrooms” fits model where house price is explained by its lot size and number of bedrooms.

In [3]:
mlr1 = smf.ols(formula="price ~ lotsize + bedrooms", data=houseprices).fit()

Fourth, we create mdatac model data object and add a constant column using add_constant function. Within add_constant function, parameters data=houseprices includes houseprices data object and prepend=False includes logical value to add constant at last column of mdatac data object. Then, we fit two stage least squares model with IV2SLS function using variables within mdatac data object and store outcome within mlr2 object. Within IV2SLS function, parameters dependent=mdatac["price"] includes model house price dependent variable, exog=mdatac[["const", "bedrooms"]] includes model number of bedrooms exogenous independent variable, endog=mdatac["lotsize"] includes model lot size endogenous independent variable, instruments=mdatac[["driveway", "garage"]] includes model whether house has a driveway and number of garage places instrumental variables, cov_type="homoskedastic" includes model homoskedastic variance covariance matrix estimation and debiased=True includes logical value to adjust model variance covariance matrix estimation for degrees of freedom. Notice that IV2SLS function parameters cov_type="homoskedastic" and debiased=True were only included as educational examples which can be modified according to your needs. Also, notice that doing stage by stage instead of simultaneous stages estimation of two stage least squares model with ols function would estimate correct coefficients but incorrect standard errors and F-statistic. Additionally, notice that two stage least squares mlr2 model estimation assumes errors are homoskedastic.

In [4]:
mdatac = sm.add_constant(data=houseprices, prepend=False)
mlr2 = lm.IV2SLS(dependent=mdatac["price"], exog=mdatac[["const", "bedrooms"]], endog=mdatac["lotsize"], instruments=mdatac[["driveway", "garage"]]).fit(cov_type="homoskedastic", debiased=True)

Fifth, we can print mlr1 model estimated coefficients table using its summary method and selecting its second tables attribute.

In [5]:
print(mlr1.summary().tables[1])
Out [5]:
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept   5612.5997   4102.819      1.368      0.172   -2446.741    1.37e+04
lotsize        6.0530      0.424     14.265      0.000       5.219       6.887
bedrooms    1.057e+04   1247.676      8.470      0.000    8116.488     1.3e+04
==============================================================================

Sixth, we can print mlr2 model estimated coefficients table using its summary method and selecting its second tables attribute.

In [6]:
print(mlr2.summary.tables[1])
Out [6]:
                             Parameter Estimates                              
==============================================================================
            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI
------------------------------------------------------------------------------
const      -1.913e+04     6540.7    -2.9248     0.0036  -3.198e+04     -6282.0
bedrooms       7680.1     1574.1     4.8791     0.0000      4588.1   1.077e+04
lotsize        12.519     1.2400     10.096     0.0000      10.084      14.955
==============================================================================

Seventh, we can print mlr1 model estimated F-statistic and associated p-value using its fvalue and f_pvalue properties.

In [7]:
print(mlr1.fvalue, mlr1.f_pvalue)
Out [7]:
159.636704936077 2.954866733815849e-55

Eighth, we can print mlr2 model estimated F-statistic and associated p-value using its f_statistic attribute.

In [8]:
print(mlr2.f_statistic)
Out [8]:
Model F-statistic
H0: All parameters ex. constant are zero
Statistic: 91.5197
P-value: 0.0000
Distributed: F(2,543)

Courses

My online courses are hosted at Teachable website.

For more details on this concept, you can view my Linear Regression in Python Course.

References

[1] Data Description: Sales prices of houses sold in the city of Windsor, Canada, during July, August and September, 1987.

Original Source: Anglin, P., and Gencay, R. (1996). Semiparametric Estimation of a Hedonic Price Function. Journal of Applied Econometrics, 11, 633–648.

[2] statsmodels Python package: Seabold, Skipper, and Josef Perktold. (2010). “statsmodels: Econometric and statistical modeling with python”. Proceedings of the 9th Python in Science Conference.

linearmodels Python package: Kevin Sheppard. (2021). “Linear (regression) models for Python. Extends statsmodels with Panel regression, instrumental variable estimators, system estimators and models for estimating asset prices”. Python package version 4.25.

My online courses are closed for enrollment.
+