*
* Example RANDOMIZ.PRG
* Manual Example 13.7
*
open data states.wks
data(org=obs,format=wks) 1 50 expend pcaid pop pcinc
order(all) pop
set pcexp = expend/pop
linreg pcexp
# constant pcaid pcinc
*
* The first test statistic is the ratio between the sum of the residuals squared over
* the first 22 observations (small states) to that over the last 22 (large states), skipping
* the middle eight. The second is the correlation between the ranks of the population
* (which is just the entry number, since the data set is sorted by population) and the
* squared residual.
*
set ressqr = %resids**2
sstats / ressqr*(t<=22)>>sum1 ressqr*(t>=29)>>sum2
compute refer1=sum1/sum2
set prank = t
order(rank=vrank) ressqr
compute refer2=%corr(prank,vrank)
*
* count1 and count2 are the number of times we get a more extreme value after
* reshuffling. We do 999 shuffles.
*
compute count1=count2=0.0
compute ns=999
do draws=1,ns
*
* Use BOOT with NOREPLACE to come up with a new permutation of the ressqr variable.
* Recompute the test statistic for these.
*
   boot(noreplace) entries 1 50
   set shuffle = ressqr(entries(t))
   sstats / shuffle*(t<=22)>>sum1 shuffle*(t>=29)>>sum2
   compute teststat=sum1/sum2
   compute count1=count1+(teststat>refer1)
   order(rank=vrank) shuffle
   compute teststat=%corr(prank,vrank)
   compute count2=count2+(teststat<refer2)
end do draws
*
* The p-values for the tests are computed by taking (count+1)/(draws+1). The +1's
* are needed because we are, in effect, adding the actual sample in with the
* draws and seeing where it ends up.
*
disp "Test 1, p-value" (count1+1)/(ns+1)
disp "Test 2, p-value" (count2+1)/(ns+1)