"""Case 005: same-sample sign reversal, NumPy + SciPy.
Self-contained deterministic dataset; no download needed.
"""
import itertools
import numpy as np
import scipy
from scipy.stats import t

rows=np.array(list(itertools.product([-1.,1.],repeat=3)))
a,b,e=rows.T
x=a
z=a+b
y=10-x+2*z+e

def ols(X):
    beta=np.linalg.lstsq(X,y,rcond=None)[0]
    df=len(y)-X.shape[1]
    cov=np.sum((y-X@beta)**2)/df*np.linalg.inv(X.T@X)
    se=np.sqrt(np.diag(cov))
    ci=np.column_stack([beta-t.ppf(.975,df)*se,beta+t.ppf(.975,df)*se])
    return beta,se,ci,2*t.sf(np.abs(beta/se),df)

simple=ols(np.column_stack([np.ones(len(y)),x]))
adjusted=ols(np.column_stack([np.ones(len(y)),x,z]))
print('NumPy',np.__version__,'SciPy',scipy.__version__)
for label,result in [('unadjusted',simple),('adjusted',adjusted)]:
    beta,se,ci,p=result
    print(label,'x estimate:',beta[1],'SE:',se[1],'95% CI:',ci[1],'p:',p[1])
Z=np.column_stack([np.ones(len(y)),z])
x_res=x-Z@np.linalg.lstsq(Z,x,rcond=None)[0]
y_res=y-Z@np.linalg.lstsq(Z,y,rcond=None)[0]
partial=x_res@y_res/(x_res@x_res)
correlation=np.corrcoef(x,z)[0,1]
vif=1/(1-correlation**2)
print('Residualized slope:',partial,'correlation:',correlation,'VIF:',vif)
print('Both x intervals include zero; this is an algebra illustration.')
assert np.isclose(simple[0][1],1) and np.isclose(adjusted[0][1],-1)
assert np.isclose(adjusted[0][2],2) and np.isclose(partial,-1)
assert np.isclose(vif,2)
assert simple[2][1,0]<0<simple[2][1,1]
assert adjusted[2][1,0]<0<adjusted[2][1,1]
print('All checks passed.')