/*
       gauss program for computing Saikkonen's test for nonlinear cointegration constant
	- dummy in the cointegration space / ar 16.7.1998 / removed 
	- concentrating the 1986 peak before estimation
*/
library pgraph, cml;
cmlset;
_nlp       = 4; /* lag length */
ofile   = "lrtest.out";
_nlr       = 3; /* cointegration rank */
_nln       = 4; /* dimension of the system */
_nle	   = 0; /* number of exogenous variables in cointegration space */
_ggamma    = 78.000369;
_nobs	   = 192-4;
load savdelta;
_gdelta    = savdelta;
/* index vectors for vectorized parameters */
_i	   = 1;
_ialpha    = seqa(_i,1,_nln*_nlr); 		_i = _i + rows(_ialpha);
_iA	   = seqa(_i,1,_nlr); 			_i = _i + rows(_iA);
if _nle >0;
	_iAexo = seqa(_i,1,_nle*_nlr); 		_i = _i + rows(_iAexo);
else;
	_iAexo = -1;
endif;
_itheta    = seqa(_i,1,_nlr); 			_i = _i + rows(_itheta);
_idelta    = seqa(_i,1,_nlr); 			_i = _i + rows(_idelta);
_iggamma   = seqa(_i,1,1); 			_i = _i + rows(_iggamma);
_itau      = seqa(_i,1,1); 			_i = _i + rows(_itau);
_ipgamma   = seqa(_i,1,_nln*_nln*(_nlp-1));

/* index vectors for data */
_i = 1;
_idy 	   = seqa(_i,1,_nln);			_i = _i + rows(_idy);
_iy1 	   = seqa(_i,1,_nlr);			_i = _i + rows(_iy1);
_iy2 	   = seqa(_i,1,1+_nle);		_i = _i + rows(_iy2);
_iz 	   = seqa(_i,1,_nln*(_nlp-1));		_i = _i + rows(_iz);

#include procs.src; /* some auxiliary procedures */

/* computes trace of a matrix */
proc (1) = tr(x);
    retp(sumc(diag(x)));
endp;
/*
    nonlin trend
*/
proc (1)=gx(theta, delta, pgamma, tau, x);
local   rx, yksikkov, retti;

    rx = rows(x);
    yksikkov = ones(rx,1);
    retti = theta*yksikkov + delta ./ (yksikkov + exp(-pgamma * (x - tau)));
    retp(retti);
endp;

/*
	another version of the above function, which is needed in computation
	of analytical gradient
*/	
proc (1)=gx2(pgamma, tau, x);
local   yksikkov, retti;

	yksikkov = 1;
    retti = yksikkov ./ (yksikkov + exp(-pgamma * (x/_nobs - tau)));
    retp(retti);
endp;

/* 
	and the correspondig derivatives
*/
proc (1)=dgx2dgamma(pgamma, tau, x);
local   yksikkov, retti;


    yksikkov = 1;
    retti = (x/_nobs - tau) .* exp(-pgamma * (x/_nobs - tau)) ./ (yksikkov + exp(-pgamma * (x/_nobs - tau)))^2;
    retp(retti);
endp;

proc (1)=dgx2dtau(pgamma, tau, x);
local   yksikkov, retti;

    yksikkov = 1;
    retti = -pgamma .* exp(-pgamma * (x/_nobs - tau)) ./ (yksikkov + exp(-pgamma * (x/_nobs - tau)))^2;
    retp(retti);
endp;
	
  
/*
  	residual function
*/  
proc (1) = residuals(par,data);
local   res, dy, alpha, y1, A, y2, gtheta, gdelta, ggamma, gtau, pgamma, rx,
        z, i, gmu;
    /*
    rearranging the parameters
clearg tmp;
tmp = par;
    */
    alpha  = reshape(par[_ialpha],_nln,_nlr);
    A      = par[_iA];
/*    A	   = diagrv(zeros(_nlr,_nlr),par[_iA]);
    if _nle >0;
    	A      = A|reshape(par[_iAexo],_nle,_nlr);
    endif;
*/    
    gtheta = par[_itheta];
    gdelta = par[_idelta];
    ggamma = par[_iggamma];
    gtau   = par[_itau];
    pgamma = reshape(par[_ipgamma],_nln,_nln*(_nlp-1));
    /*
    rearranging the data
    */
    dy = data[.,_idy]';
    y1 = data[.,_iy1]'; /* this is not very clear to me; it is assumed that r = 1 */
    y2 = data[.,_iy2]';     /* */
    rx = rows(dy');
    gmu = gx(gtheta[1], gdelta[1], ggamma, gtau, seqa(0,1/rx,rx));
	i=2;
	do while i<=_nlr;
    	gmu = gmu~(gx(gtheta[i], gdelta[i], ggamma, gtau, seqa(0,1/rx,rx)));
		i=i+1;
	endo;
    gmu = gmu';
    z = data[.,_iz]';
/*
print;
print "dy     " dims(dy);
print "alpha  " dims(alpha);
print "y1     " dims(y1);
print "A      " dims(A);
print "y2     " dims(y2);
print "gmu    " dims(gmu);
print "pgamma " dims(pgamma);
print "z      " dims(z);
print "gmus   " dims(gmu);
print "eka    " dims(ones(3,1)*y1);
print "ecm    " dims(alpha*(ones(3,1)*y1 - A'y2 - gmu ));
*/
    res = dy - alpha*(eye(3)*y1 - A*y2 - gmu) - pgamma * z;
    retp(res');
endp;
/*
    likelihood function
*/
proc (1)=suf(par, data);
local   res, rx, cova, invcova, detcova, loglik, cs, fll;

    rx = rows(data);
    res = residuals(par,data);
    cova = moment(res,0)/rx;
    detcova = det(cova);
    fll = -rx*ln(detcova)/2 - rx*rows(cova)/2;
    loglik = ones(rx,1).*fll./rx;
    retp(loglik);
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];
/*
	printout of headers
*/
output file=^ofile reset;
print "FILE: " $ofile "  DATE: " datestr(date()) "  TIME: " timestr(time());
print "by Antti Ripatti, September 29, 1997";
print;
print "=============================================================================";
print "        ESTIMATION OF CONTINUOUS STRUCTURAL CHANGES IN INTERCEPT TERMS";
print "=============================================================================";
print;
print "Note that the 1986:8 spike has been removed";
print;
format /m0 /lD 8,0;
print "Variables: " $labs';
print;
format /m0 /lD 1,0;
print "Lag length in VECM: " _nlp-1 "     Cointegration rank: " _nlr;
print;
print "Number of observations: " rows(dat)-_nlp;
print;
print "delta(2) and delta(3) fixed to zero!";
print;
output file=^ofile off;
/*
	concentrating 1986 spike
*/
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);
/*
    creating data
*/
dy = createZ0(yt,_nlp);
z = createZ1(yt,_nlp,0);
ystar = createZk(yt,_nlp,0);
y1 = ystar[.,1 2 3];
y2 = ystar[.,4];
clearg data;
data = dy~y1~y2~z;
/*
    computing initial values
*/
clearg l, beta, alpha, tracet, A, gmus, GGamma, _gpar0, ll;
{l, beta, alpha, tracet} = urrr(dy,z,ystar,_nln);
beta = beta[.,1:_nlr];
alpha = alpha[.,1:_nlr];
A  = diagrv(zeros(_nlr,_nlr),(1|1|1));
/*
A  = A|(0.07~0.07~0.07);
*/
gmus = ones(6,1)|(0.42*ones(2,1));
clearg selitt, resids, omega;
selitt = (ystar*beta)~z;
GGamma = inv(selitt'selitt)*selitt'dy;
GGamma = GGamma[_nlr+1:rows(GGamma),.]';
_gpar0 = vecr(alpha)|(1|1|1)|gmus|vecr(GGamma);
/********************************************************************

    computing the parameter names
    
********************************************************************/
clearg PNalpha, PNA, PNtheta, PNdelta, PNgamma, PNtau, PNGamma, ro, co, Pgmus;
i=1;
PNalpha = alpha;
do while i<=_nln;
    j=1;
    do while j<=_nlr;
        PNalpha[i,j]="à(" $+ ftos(i,"%*.*lf",1,0) $+ "," $+ ftos(j,"%*.*lf",1,0) $+ ")";
        j=j+1;
    endo;
    i=i+1;
endo;
i=1;
PNA = A;
ro = rows(A); co=cols(A);
do while i<=ro;
    j=1;
    do while j<=co;
        PNA[i,j]="A(" $+ ftos(i,"%*.*lf",1,0) $+ "," $+ ftos(j,"%*.*lf",1,0) $+ ")";
        j=j+1;
    endo;
    i=i+1;
endo;
PNA = "A1"|"A2"|"A3";
i=1;
PNGamma = GGamma;
ro = rows(PNGamma); co=cols(PNGamma);
do while i<=ro;
    j=1;
    do while j<=co;
        PNGamma[i,j]="â(" $+ ftos(i,"%*.*lf",1,0) $+ "," $+ ftos(j,"%*.*lf",1,0) $+ ")";
        j=j+1;
    endo;
    i=i+1;
endo;
i=1;
Pgmus = gmus;
do while i<=_nlr;
	Pgmus[i] = "phi(" $+  	ftos(i,"%*.*lf",1,0) $+ ")";
	i=i+1;
endo;
do while i<=2*_nlr;
	Pgmus[i] = "delta(" $+  ftos(i-_nlr,"%*.*lf",1,0) $+ ")";
	i=i+1;
endo;
Pgmus[i] = "gamma";
i=i+1;
Pgmus[i] = "tau";

_cml_ParNames = vecr(PNalpha)|vecr(PNA)|vecr(Pgmus)|vecr(PNGamma);
/*****************************************************************************

    T e s t i n g

*****************************************************************************/
clearg mlagnr, mlag;
load mlagnr;
load mlag;
nobs = rows(data);
sufmlag = sumc(suf(mlag,data));
sufmlagnr = sumc(suf(mlagnr,data));
lrtest = 2*sufmlagnr -2*sufmlag;
df = 2;
output file=^ofile on;
print "-----------------------------------------------------------";
print "Likelihood ratio test statistic for Ho: delta(2)=delta(3)=0";
print "-----------------------------------------------------------";
print;
format /m1 /rd 12,6;
print "Likelihood function (unrestricted):" sufmlagnr;
print "Likelihood function (restricted):" sufmlag;
format /m1 /rd 2,0;
print "LR test:  chi^2("  df;;
format /m1 /rd 12,6;
print ")= " lrtest "    (" cdfchic(lrtest,df) ")";
print;
end;

end;