/********************************************************
Author: Ralf Becker <r.becker AT qut.edu.au>
Date: 12/29/2003
Caveat: free for noncommercial use.
this procedure calculates the Nyblom (1992)
test
input:	(i) y, vector of depenedent var
		(ii) x, matrix of indep. var
				MAX DIM OF 4 !!!
output: (i)	3/2/1/0, sig at 1%/%5/10%/not sig
********************************************************/

proc nyblom(y,x);
        local e,sigsq_hat,xe,s,nyb,i,cv,sig;
e = y - x*(y/x);
sigsq_hat = meanc(e.*e);
xe = x.*e;
s = x'x/t;
nyb = zeros(cols(x),cols(x));

for i(1,t,1);
  nyb = nyb + sumc(xe[i:t,.])*sumc(xe[i:t,.])';
endfor;

nyb = sumc(diag(inv(s)*nyb))/(sigsq_hat*t^2);

			/* critical values
				col 1/2/3: 10%/5%/1%
				row 1/2/3/4: 1/2/3/4 d.o.f */
cv = {0.347 0.461 0.743, 0.607 0.748 1.074,
		0.841 1.000 1.359, 1.063 1.237 1.623};

cv = cv[cols(x),.];
sig = sumc(nyb.>cv');

retp(sig);
endp;