new;
Tvec = {25, 50, 100}; tau = 0.5; /* tau = t0/T*/
betavec = {0.3, 0.6, 0.9}; 
alpha = 0.05; /*nominal size of test*/
nrep = 100; /*increase to 10,000 when program working*/
seed1=161058;
count = zeros(rows(Tvec), rows(betavec)); /*storage for rejection counts*/
it = 1;
do until it > rows(Tvec);
 	T = Tvec[it]; t0 = tau*T;
	ib=1;
	do until ib > rows(betavec);
		beta=betavec[ib];
/*now generate nrep autoregressions of length T, coefficient beta, SP y0*/
		in = 1;
		do until in > nrep;
			eps = rndns(T+1,1,seed1);
			y0=eps[1]/sqrt(1-beta*beta);
			y = recserar( eps,y0,beta);
			xy = y[1:T]~y[2:T+1];
				/*3 alternative Cholesky factorisations of xy'xy */
				R1 = qr(xy[1:t0,.]); /* obs 1 to t0 */
				R2 = qr(xy[t0+1:T,.]); /*obs t0+1 to T */
				R = qr(R1|R2); /*obs 1 to T */
			urss = R1[2,2]^2 + R2[2,2]^2; rrss = r[2,2]^2;
			F = (T-2)*(rrss - urss)/urss;
			count[it,ib] = count[it,ib] + ( cdffc(F,1,T-2) < alpha );
			in = in + 1;
		endo;
		ib = ib +1;
	endo;
	it = it+1;
endo;
format 12,3;
"rejection %" betavec';
Tvec~( 100*count/nrep);
format 20,15; "newseed" seed1; format 14,6;/*return format to default*/