# 1. Estimation

beta_hat=inv(X.T@X)@X.T@yarray([[1.31986603e+02],       [9.49700802e+00],       [5.09767067e-02],       [4.56077098e-01]])

# 2. Construct covariance matrix

iden=np.identity(100,dtype=float)M_x=iden-X@inv(X.T@X)@X.Tu_hat_un=M_x@y
def fvar_ols(u_hat_DGP):   s2=(u_hat_DGP.T@u_hat_DGP)/(n-k)   var_ols=s2*inv(X.T@X)   return var_olsdef fHC0(u_hat_DGP):   u_sq=np.square(u_hat_DGP)   omega0=np.diagflat(u_sq)   HC0=inv(X.T@X)@(X.T@omega0@X)@inv(X.T@X)   return HC0P_x=X@inv(X.T@X)@X.Tht=np.diag(P_x)def fHC2(u_hat_DGP):   omega2_diag=[0 for i in range(n)   u_hat_DGP_sq=np.square(u_hat_DGP)   for i in range(n):      omega2_diag[i]=u_hat_DGP_sq[i]/(1-ht[i])   omega2=np.diagflat(omega2_diag)   HC2=inv(X.T@X)@(X.T@omega2@X)@inv(X.T@X)   return HC2

# 3. Wald Statistic and P values

R=np.matrix([[0,0,1,0],[0,0,0,1]])def coef(X,y):    beta_hat=inv(X.T@X)@X.T@y    return beta_hatfrom scipy.stats import chi2def p(var,X,y):   var_hat=var   beta_hat=coef(X,y)   Wald=(R@beta_hat).T@inv(R@var_hat@R.T)@(R@beta_hat)   pvalue=1-(chi2.cdf(Wald,2))   return np.array((Wald,pvalue)).reshape(2)

# V. Conclusion

test_hat=np.array((p(var_ols,X,y),p(HC0,X,y),p(HC2,X,y)),dtype=float).reshape(3)

--

--