' two-sample permutation test
' (Johnston-DiNardo 11.2)
' matrix permutation function (fast)
' for version 4
' 11/3/99 h
' revised and last checked 9/19/2000 h

' create workfile
workfile cardkrug u 1 1

' input data in matrix (first 33 rows from N.J.)
matrix(40) data
data.fill -20, -17.5, -13, -12.5, -4.5, -4, -3.5, -2, -1.5, -1, -0.5, -0.5, 0, 0, 0.5, 0.5, 1.5, 2, 2, 2.25, 3, 4.5, 4.5, 5.5, 6, 6.25, 8.25, 9, 10, 10.5, 12, 14.75, 34, -7, -6, -2.5, -0.5, 4, 4.5, 4.5

' extract submatrices for each sample
matrix sub0 = @subextract(data,1,1,33,1)
matrix sub1 = @subextract(data,34,1)

' compute difference in means for actual sample
scalar mdiff = @mean(sub0) - @mean(sub1)

' set number of replications
!reps = 1000

' set seed of random number generator
rndseed 123456

' declare storage matrices
matrix pdata				' permuted data
vector(!reps) pdiff			' difference in means from pdata

' permutation loop
for !i=1 to !reps
	' permute data
	pdata = @permute(data)
	' extract submatrices from permuted data
	sub0 = @subextract(pdata,1,1,33,1)
	sub1 = @subextract(pdata,34,1)
	' store difference in means
	pdiff(!i) = @mean(sub0) - @mean(sub1) 
next

' 95% interval (percentile method)
scalar lower = @quantile(pdiff,0.025) 
scalar upper = @quantile(pdiff,0.975)

' display output in table
table out
setcolwidth(out,1,30)
out(1,1) = "Sample difference in means:"
out(1,2) = mdiff

out(2,1) = "Permutation confidence interval:"
out(2,2) = "[" + @str(lower) + "," + @str(upper) + "]" 

out(3,1) = "Number of permutations:"
out(3,2) = @str(!reps)
show out

