@Program name: reg_res.prg @ 
@Purpose: construct a bootstrap estimated sampling distribution 
of the slope coefficient of a bivariate regression @ 
@Used for: 1996 course in non-parametric inference, Essex Summer School @ 
@Author: Chris Mooney @ 
@Date: June 20, 1996 @ 

new; 
@rndseed 34534;@ 
t1=hsec;@for seting the timer@ 
output file = reg_res.out; @Defining output file@ 
output reset; @Overwrite previous output file @ 
library pgraph; 
graphset; 

@Set parameters @ 
resamp = 2; @set re-sampling technique: 1 = residuals, 2 = cases@ 
n=40; @sample size @ 
b=10; @number of bootstrap re-samples @ 
e_weight = 25; @error weight@ 
tbeta = 2.0; tcons = 1.0; @pseudo-pop parameters for generating regression data@ 

i=1; 
bb2 = zeros(B,1); @Setting up empty vector to hold bootstrapped beta-hat's@ 

@Setting up the pseudo-population data @ 
x1=ones(N,1); @Generating vector of ones for constant term in x matrix @ 
x2=seqa(1,1,N); @Generating a fixed X2@ 
@x2=rndn(n,1) * n;@ @Generating random X2@ 
x=x1~x2; @Setting up X data matrix@ 
e = rndn(n,1); @Generating random error @ 
y=(tbeta*x2) + tcons+ (e*e_weight) ; @Generating the dependent variable in 
the sample@ 

{bet,sd,r2}=regress(x,y); @Regress full sample Y on X@ 
resid= y - (bet[1,.] + bet[2,.]*x2); @ save residuals for resamples@ 



@ RESAMPLING@ 
do while i<=b; @Starting the re-sampling loop@ 
cindex=ceil(rndu(N,1)*N); @Setting up cases re-sampling index@ 

xc=x[cindex,.];yc=y[cindex,1]; @Re-sampling cases @ 
re=submat(resid,cindex,0); @Re-sampling residuals@ 
yr=x*bet+re; @Setting up re-sampled Y vector for residuals re-sampling@ 

if resamp==1; @Responding to re-sampling flag@ 
{b11,se,r2}= regress(x,yr); @Residuals re-samplng regression@ 
else; @Responding to re-sampling flag@ 
{b11,se,r2}=regress(xc,yc); @Cases re-sampling regression@ 
endif; @end of loop choosing which type of re-sampling@ 

bsb2=b11[2,.]; @Identify the regression coefficient@ 
bb2[i,1] = bsb2; @Setting up vector of bootstrapped beta-hats@ 

i=i+1; @increment bootstrap looping index@ 
endo; @end of resampling loop@ 

boot_b=meanc(bb2); @Gives the average of the beta-hat*'s, 
the bootstrap estimate of @ 
sd_b=stdc(bb2); @Gives the standard deviation of the beta-hat*'s, 
the bootstrap estimate of @ 
bb2=sorthc(bb2,1); @Sorting the vector of *'s @ 
ll = ceil((b*0.05)/2); @Setting up lower and upper 
percentile points at .025 of b @ 
ul= ceil((b-ll)+1); 
low_b = bb2[ll,.]; @Gives the value of beta-hat* at the lower 
percentile point set by ll @ 
up_b = bb2[ul,.]; @Gives the value of beta-hat* at the upper 
percentile point set by ul @ 

print; 
print "sample size =" n; 
print "number of bootstrap re-samples =" b; 
print "re-sampling method (1=residuals, 2=cases) =" resamp; 
print "pseudo-population slope =" tbeta; 
print; 
print "R-squared =" r2; 
print "OLS b =" bet[2,.]; 
print "OLS sd(b) =" sd[2,.]; 
print; 
print "Mean b* = " boot_b; 
print "Sd of b* =" sd_b; 
print "2.5th percentile point of b* = " low_b; 
print "97.5th percentile point of b* = " up_b; 
print; 
print "Minutes elapsed = " (hsec-t1)/6000; 
print; 

@REGRESS proc included to perform OLS@ 

proc(3) = regress (x,y); 
local xxi,beta,e,sse,sd,t,r2,A,il,n1; 
n1=rows(x); 
XXi=inv(X'X); 
beta=XXi*(X'Y); 
e=Y-X*beta; 
sse=e'e/(rows(X)-cols(x)); 
sd=sqrt(diag(sse*XXi)); 
t=beta./sd ; 
il=ones(n1,1); 
a=(eye(N1)-il*il'/n1); 
r2= (1-e'e/(y'*A*y)); 
retp(beta,sd,R2); 
output off; 
endp;