*
* RANDOMIZE.RPF
* RATS Version 8, User's Guide, Example 16.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 draw=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 draw
*
* 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)