new;
nrep = 100; /*increase to 10,000 in production*/
seed1 = 150343;
Tvec = {25, 50, 100, 200};
alpha = 0.05;
K = 4; K2=K+2;
B = 40; @increase, say to 400, in production@
@one could experiment with different K, setting up a Kvec and a loop on ik@
format 6,0; "nrep = " nrep "K " K ; format 20,15; "seed" seed1; format 14,6;
count=zeros(rows(Tvec),1);
tim = hsec; @to time the experiment@
it = 1;
do until it > rows(Tvec);
	T = Tvec[it];
	X = ones(T,1)~rndus(T,K-1,seed1);
	irep = 1;
	do until irep > nrep;
		eps = rndn(T,1);
		e = eps - X*(eps/X); @GAUSS's odd way of calculating b@
		r = qr( X~(0|e[1:T-1])~e); @Cholesky again@
		F = (r[K2-1,K2]/r[K2,K2])^2; @can be shown to give the F=t^2 on e(t-1)@
		Fvec=zeros(B,1);
		ib = 1;
		do until ib > B;
			irow = round(T*rndus(T,1,seed1)+0.5);
			estar=e[irow];
			Xstar=X[irow,.];
			rstar = qr( Xstar~(0|estar[1:T-1])~estar); @Cholesky again@
			Fstar = (rstar[K2-1,K2]/rstar[K2,K2])^2;
			Fvec[ib] = Fstar;
			ib = ib+1;
		endo;
		pstar = 1 - sumc(F .> Fvec)/B;
		count[it] = count[it] + (pstar < alpha);
		irep = irep+1;
		endo;
	it = it+1;
endo;
tim = hsec - tim;
"sizes";
Tvec~(count/nrep);
"time" 0.01*(tim);