/*
       gauss program for computing Saikkonen's test for nonlinear cointegration constant
*/
library gauss, pgraph;
clearg yt, labs, dy, z, ystar;
p       = 4; /* lag length */
ofile   = "nonlin1.out";
r       = 3;
n       = 4;

#include procs.src;
/* computes trace of a matrix */
proc (1) = tr(x);
    retp(sumc(diag(x)));
endp;

proc (0) = nltest(omega, alpha, res, iNi );
       local         left, nobs, vars, W, i, selitt, omegahat, omegatilde,
                     resids, S, df;
       
       Left = Alpha'inv(Omega)*Res';
       Left = Left';
       Nobs = rows(Res);
       Vars = cols(Res);
       /* creating w(t) */
       W = seqa(0,1/nobs,nobs);
/* print "W: " rows(w) " x " cols(w); */
       i=2;
       do while i<=iNi;
              W = W~(W.^i);
              i = i+1;
       endo;
       selitt = W~ystar~z;
       resids = left - selitt*inv(selitt'selitt)*selitt'left;
        omegahat = resids'resids/rows(resids);
/* print "Omega hat: " omegahat; */
        omegatilde = alpha'inv(omega)*alpha;
/* print "Omega tilde: " omegatilde; */
        S = nobs*tr(inv(omegatilde)*(omegatilde - omegahat));
       df = ini*cols(alpha);
    format /m1 /rd 2,0;
       print "Saikkonen's nonlinearity test: chi^2(" df ;;
    format /m1 /rd 12,6;
       print ")= " S "    (" cdfchic(S,df) ")";
/*
format /re 20,15;
print "Test statistic: " S;
*/
endp;
/*
    loading data file and data labels
*/
clearg dat, datlab;
dat     = loadd("igauss");
datlab  = getname("igauss");
clearg yt, labs, dy, z, ystar, y1, y2, toinen;
toinen  = 4;
yt 		= dat[.,toinen 8 11 3];	 /* 11 = i3m and 4 = iown */
labs    = datlab[toinen 8 11 3];
clearg dum86;
dum86m8 = zeros(rows(yt),1);
dum86m8[80] = 1;
dum86m8 = demean(dum86m8,ones(rows(dum86m8),1));
dum86m9 = zeros(rows(yt),1);
dum86m9[81] = 1;
dum86m9 = demean(dum86m9,ones(rows(dum86m9),1));
yt = deseas(yt,dum86m8);

output file=^ofile reset;
print "FILE: " $ofile "  DATE: " datestr(date()) "  TIME: " timestr(time());
print "by Antti Ripatti, August 19, 1997";
print;
print "=============================================================================";
print "        TESTING FOR CONTINUOUS STRUCTURAL CHANGES IN INTERCEPT TERMS";
print "=============================================================================";
print;
format /m0 /lD 8,0;
print "Variables: " $labs';
print;
format /m0 /lD 1,0;
print "Lag length in VECM: " p "     Cointegration rank: " r;
print;
print "Number of observations: " rows(dat)-p;
print;
output file=^ofile off;
format /m1 /rd 12,6;


dy = createZ0(yt,p);
z = createZ1(yt,p,0);
ystar = createZk(yt,p,1);
clearg l, beta, alpha, tracet;
{l, beta, alpha, tracet} = urrr(dy,z,ystar,n);
beta = beta[.,1:r];
alpha = alpha[.,1:r];
clearg selitt, resids, omega;
selitt = (ystar*beta)~z;
resids = dy - selitt*inv(selitt'selitt)*selitt'dy;
omega = resids'resids/rows(resids);
output file=^ofile on;
print "Test statistics computed using auxiliary regressions_______________";
print;
nltest(omega, alpha, resids, 1);
nltest(omega, alpha, resids, 2);
nltest(omega, alpha, resids, 3);
output file=^ofile off;
/*
	computing another LM type test
*/
proc (1)= analgrad(alpha,beta,y,zstar,pp);
	local 	A, ualpha, rx,  PPi, i, z, F1t, F2t, F3t, F1, F2, F3, F, wt, j, u;
        

	rx		= rows(y);
	A		= beta*inv(beta[2:n,.]);
	A		= A./A[1,.];
	PPi		= alpha*beta';
	ualpha  = PPi*A*inv(A'A);
	y2		= y[.,2:n]';
	u 		= y*A;
	z		= u'|zstar';
	j = 1;
	wt		= seqa(0,1/rx,rx);
	for j (2,pp,1);
		wt = wt~(seqa(0,1/rx,rx)^j);
	endfor;
	wt = wt';
	i = 1;
	F1t = 	(y2[.,i].*.ualpha')|
			(alpha');
    F2t = -(z[.,i].*.eye(n));
    F3t	= wt[.,i].*.ualpha';
	F = F1t|F2t|F3t;
	i = 2;
	for i (2,rx,1);
		F1t = 	(y2[.,i].*.ualpha')|
				(alpha');
	    F2t = -(z[.,i].*.eye(n));
    	F3t	= wt[.,i].*.ualpha';
		F = F~(F1t|F2t|F3t);
	endfor;
/* "F" dims(F); */
	retp(F');
endp;

proc (1) = lmstat(F,Res);
	local n, Ft, omega, invomega, nobs, score, scoret, cova, covat, testst, i;

	n	= cols(res);
	nobs = rows(res);
	Ft = F[1:n,.];
	omega = moment(res,0)/nobs;
	invomega = inv(omega);
	covat = Ft'invomega*Ft;
	cova = covat;
	i = 1;
	scoret = Ft'invomega*res[i,.]';
	score = scoret;
	for i (2,nobs,1);
		Ft = F[n*(i-1)+1:n*(i-1)+n,.];
		covat = Ft'invomega*Ft;
		cova = cova+covat;
		scoret = Ft'invomega*res[i,.]';
		score = score + scoret;
	endfor;
	testst = score'inv(cova)*score;
	retp(testst);
endp;

clearg F, pp, df;
pp = 1;
F = analgrad(alpha,beta,ystar,z,pp);
testst = lmstat(F,resids);
df = r*pp;
output file=^ofile on;
print "Test statistics computed using quadratic forms_______________";
print;
format /m1 /rd 2,0;
print "Saikkonen's nonlinearity test: chi^2("  df;;
format /m1 /rd 12,6;
print ")= " testst "    (" cdfchic(testst,df) ")";
print;
/*
format /re 20,15;	   
print "Test statistic: " testst;
*/
output file=^ofile off;
pp = 2;
F = analgrad(alpha,beta,ystar,z,pp);
testst = lmstat(F,resids);
df = r*pp;
output file=^ofile on;
print;
format /m1 /rd 2,0;
print "Saikkonen's nonlinearity test: chi^2("  df;;
format /m1 /rd 12,6;
print ")= " testst "    (" cdfchic(testst,df) ")";
print;
/*
format /re 20,15;	   
print "Test statistic: " testst;
*/
output file=^ofile off;

end;