# Discrepancy in 2SLS Estimation Results

I've been exploring the 2SLS (Two-Stage Least Squares) estimation method to analyze a model involving endogeneity and instrumental variables. To better understand the process, I performed manual calculations using NumPy to implement the two-stage procedure. However, I noticed a discrepancy when comparing my results with those obtained using the linearmodels library in Python.

import numpy as np
import pandas as pd
from linearmodels.system import SUR

data = pd.DataFrame({
'Y1': [4, 1, 1, 2],
'Y2': [1, 2, 1, 5],
'X1': [1, 1, 3, 2],
'X2': [2, 5, 2, 7]
})

equations = {
'eq1': 'Y1 ~ 1 + Y2 + X2',
'eq2': 'Y2 ~ 1 + Y1 + X1'
}

model = SUR.from_formula(equations, data=data)

results = model.fit()

print(results)


Here's what I did in numpy:

being the Matrix of Cross Products this:

and the equations:

I carefully set up the model equations and organized the data into matrices for Y1, Y2, X1, and X2. I conducted the first-stage regressions to obtain the fitted values of Y1 and Y2 using NumPy. Then, I performed the second-stage regressions by replacing the fitted values in the original equations. However, the parameter estimates I obtained from my NumPy calculations do not match the results I got when using the linearmodels library, which implements the same 2SLS method. I ensured that the steps and matrices are consistent in both implementations.

import numpy as np

data = np.array([
[4, 1, 1, 2],
[1, 2, 1, 5],
[1, 1, 3, 2],
[2, 5, 2, 7]
])

Y1 = data[:, 0]
Y2 = data[:, 1]
X1 = data[:, 2]
X2 = data[:, 3]

# First-stage regressions
X_1st = np.column_stack((np.ones(len(X1)), X1))
X_2nd = np.column_stack((np.ones(len(X2)), X2))

coeff_1st_Y1 = np.linalg.inv(X_1st.T @ X_1st) @ X_1st.T @ Y1
coeff_1st_Y2 = np.linalg.inv(X_2nd.T @ X_2nd) @ X_2nd.T @ Y2

Y1_fitted = X_1st @ coeff_1st_Y1
Y2_fitted = X_2nd @ coeff_1st_Y2

# Second-stage regressions
X_eq1 = np.column_stack((np.ones(len(Y1_fitted)), Y2_fitted, X2))
X_eq2 = np.column_stack((np.ones(len(Y2_fitted)), Y1_fitted, X1))

coeff_eq1 = np.linalg.inv(X_eq1.T @ X_eq1) @ X_eq1.T @ Y1
coeff_eq2 = np.linalg.inv(X_eq2.T @ X_eq2) @ X_eq2.T @ Y2

print("Equation 1:")
print("Intercept:", coeff_eq1[0])
print("Y2 coefficient:", coeff_eq1[1])
print("X2 coefficient:", coeff_eq1[2])

print("\nEquation 2:")
print("Intercept:", coeff_eq2[0])
print("Y1 coefficient:", coeff_eq2[1])
print("X1 coefficient:", coeff_eq2[2])


Can anyone shed light on why there might be a discrepancy between my manual calculations with NumPy and the results from the linearmodels library? Are there any potential pitfalls or nuances in the 2SLS procedure that could explain this difference?