# II. Resampling

#construct H0 modelZ=np.column_stack((iota,x_1))gamma_hat=inv(Z.T@Z)@Z.T@y #run H0M_z=iden-Z@inv(Z.T@Z)@Z.Tu_hat_r=M_z@y# This function takes any residual vector!!def resample_DGP(residual):    u_boot=np.random.choice(residual.flatten(),n,replace=True)  .reshape(n,1)     y_resample=np.add(Z@gamma_hat,u_boot)     return y_resample
def 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)var_ols=fvar_ols(u_hat_un)p(var_ols,X,y)
test_hat=np.array((p(var_ols,X,y),p(HC0,X,y),p(HC2,X,y)),dtype=float).reshape(3)
def test(var_hat,X,y):    beta_hat=coef(X,y)    t=(R@beta_hat).T@inv(R@var_hat@R.T)@(R@beta_hat)    return t
def resample_ols(B):    p=0    for i in range(B):        y_resample=resample_DGP(u_hat_r)        w=test(fvar_ols(M_x@y_resample),X,y_resample)        I=bool(w>test_hat)        p=p+I    return p/B

# III. Wild Bootstrap

from math import sqrtet=[-((sqrt(5)-1)/2), (sqrt(5)+1)/2]prob=(sqrt(5)+1)/(2*sqrt(5))print(et*prob+et*(1-prob), (et**2)*prob+(et**2)*(1-prob)) # expectation zero variance 1def wild_DGP(residual):    st=np.random.choice(et,n,p=[prob,1-prob]).reshape(n,1)    u_boot=np.multiply(st,residual)    y_wild=np.add(Z@gamma_hat,u_boot) #gamma_hat from running H0    return y_wild
def wild_ols(B):    p=0    for i in range(B):        y_wild=wild_DGP(u_hat_r)        w=test(fvar_ols(M_x@y_wild),X,y_wild)        I=bool(w>test_hat)        p=p+I    return p/B

--

--