new;
Tvec = {25, 50, 100}; 						@vector of T values	  	 @
tau = 0.5;  								@break at t0, tau = t0/T @
betavec = {0.3, 0.6, 0.9}; 					@vector of beta values	 @
alpha = 0.05;  								@nominal size of test	 @
nrep = 100;					 @increase to 10,000 when program working@
seed1=161058;				   @initial seed for random no. generator@
count = zeros(rows(Tvec), rows(betavec)); 	@storage for rejection counts@

it = 1;										@loop on T values		 @
do until it > rows(Tvec);
 	T = Tvec[it]; t0 = int(tau*T);				@extract T, set break t0 @
	t1 = seqa(1,1,t0); t2 = seqa(t0+1,1,T-t0);  @time periods@
	ib=1;									@loop on beta values	 @
	do until ib > rows(betavec);
		beta=betavec[ib];					@extact beta			 @
/*now generate nrep autoregressions of length T, coefficient beta, SP y0*/
		in = 1;
		do until in > nrep;					@replication loop		 @
			eps = rndns(T+1,1,seed1);		@random error vector	 @
			y0=eps[1]/sqrt(1-beta*beta);	@ y0, correct variance	 @
			y = recserar( eps,y0,beta);		@ y(t) = beta*y(t-1)+eps(t)@
			yt = y[2:T+1]; yt_1=y[1:T];		@ y(t), y(t-1)			 @


/* regress y(t) on y(t-1) for whole period, 2 sub-periods, find rss  */
			e  = yt     - yt_1*(yt/yt_1);             rss = e'e; 
			e1 = yt[t1] - yt_1[t1]*(yt[t1]/yt_1[t1]); rss1 = e1'e1;
			e2 = yt[t2] - yt_1[t2]*(yt[t2]/yt_1[t2]); rss2 = e2'e2;
/*construct F test */
			urss = rss1 + rss2 ;
			F = (T-2)*(rss - urss)/urss;
/*add to count matrix if tail probability < nominal*/
			count[it,ib] = count[it,ib] + ( cdffc(F,1,T-2) < alpha );
			in = in + 1;
		endo;								@end of replication loop@
		ib = ib +1;
	endo;									@end of beta loop		@
	it = it+1;
endo;										@end of T loop			@


format 12,3;
"rejection %";
"            Beta";
"         T "   betavec';					@ print results heading	 @
Tvec~( 100*count/nrep);						@ print results table 	 @
format 20,15; "newseed" seed1; 				@save seed				 @
format 14,6;								@return format to default@