frisch_waugh_lovell.py

"""
Frisch-Waugh-Lovell Theorem.

In a multiple regression, the effect of any single covariate can be obtained 
by removing the effects of the other covariates from both the dependent and 
independent variables.
"""
import numpy as np
import statsmodels.api as sm

def gen_data(nobs, nvars, seed=None):
    np.random.seed(seed)
    X = np.random.random((nobs, nvars))
    X = sm.add_constant(X)
    true_beta = np.random.uniform(1,nvars, size=nvars+1)
    y = np.dot(X, true_beta) + np.random.randn(nobs)
    return y, X, true_beta
 
def gen_data_trend(nobs, seed=None):
    np.random.seed(seed)
    X = np.random.random((nobs, 1))
    X = sm.tsa.add_trend(X, 'ct', prepend=True)
    true_beta = np.random.uniform(1, 5, size=3)
    y = np.dot(X, true_beta) + np.random.randn(nobs)
    return y, X, true_beta



if __name__ == '__main__':
    y, X, true_beta = gen_data(100, 10, seed=0)
    X1 = X[:,:-1]
    Xk = X[:,-1]
     
    # remove the effects of x1 from y
    fit1 = sm.OLS(y, X1).fit()
    gamma1 = fit1.params
    f = fit1.resid
     
    # remove the effects of x1 from Xk
    fit2 = sm.OLS(Xk, X1).fit()
    gamma2 = fit2.params
    g = fit2.resid
     
    gamma3 = np.dot(g,f)/np.dot(g,g)
     
    beta_fwl = np.r_[gamma1 - gamma2*gamma3, gamma3]
    beta = sm.OLS(y, X).fit().params
     
    np.testing.assert_array_almost_equal(beta_fwl, beta)

    print sm.OLS(f, g).fit().summary()

    print sm.OLS(y, X).fit().summary()

    # ---------
    # As detrending.

    yy, XX, bb = gen_data_trend(100)

    ols_wtrend = sm.OLS(yy, XX).fit()

    f1 = sm.OLS(XX[:,-1], XX[:, :2]).fit().resid
    g1 = sm.OLS(yy, XX[:, :2]).fit().resid

    ols_fwl = sm.OLS(g1, f1).fit()

    print ols_wtrend.summary()
    print ols_fwl.summary()