nulldata 50
set seed 54321
series x = normal()
series y = 10 + x + 2*normal()
ols y 0 x
# the reported t-stat
t0 = abs($coeff[2] / $stderr[2])
# save the residuals
series u = $uhat
scalar ybar = mean(y)
# number of replications for bootstrap
scalar B = 1000
scalar tcount = 0
series ysim
loop B
  # generate simulated y by resampling
  ysim = ybar + resample(u)
  ols ysim 0 x --quiet
  scalar tsim = abs($coeff[2] / $stderr[2])
  tcount += (tsim > t0)
endloop      
printf "proportion of cases with |t| > %.3f = %g\n", t0, tcount / B