/*
       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 gauss, pgraph, cml;
cmlset;
_nlp       = 4; /* lag length */
ofile   = "mlagnr.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;

/* 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;
/**********************************************************************
**
** + - - - - - - - - - -   ------------------+
** | A n a l y t i c a l   gradient function |
** + - - - - - - - - - -   ------------------+
**
**********************************************************************/
proc analgrad(par,data);
	local 	res, dy, alpha, y1, A, y2, gtheta, gdelta, ggamma, gtau, 
			pgamma, rx, z, i, gmu, F1t, F2t, F1, F2, F, e1, ecm;
        

    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;
	ecm = eye(_nlr)*y1 - A*y2 - gmu';	
	z = data[.,_iz]~(ecm');
	z = z';
	e1 = 1|0|0;
	i = 1;
    F2t = -(z[.,i].*.eye(4));
	_nobs = rx;
	F1t = 	(y2[.,i]*alpha')|
			(alpha')|
			(gx2(ggamma,gtau,i)*alpha')|
			(dgx2dgamma(ggamma,gtau,i)*gdelta'alpha')|
			(dgx2dtau(ggamma,gtau,i)*gdelta'alpha');
	F = F2t[_nln*(_nln-1)*_nlp+1:rows(F2t),.]|F1t|F2t[1:_nln*(_nln-1)*_nlp,.];
	i = 2;
	do while i<=_nobs;
		F1t = 	(y2[.,i]*alpha')|
				(alpha')|
				(gx2(ggamma,gtau,i)*alpha')|
				(dgx2dgamma(ggamma,gtau,i)*gdelta'alpha')|
				(dgx2dtau(ggamma,gtau,i)*gdelta'alpha');
    	F2t = -(z[.,i].*.eye(4));
		F = F~(F2t[_nln*(_nln-1)*_nlp+1:rows(F2t),.]|F1t|F2t[1:_nln*(_nln-1)*_nlp,.]);
		i=i+1;
	endo;
"F" dims(F);
	retp(F');
endp;

proc lmtest(par,data);
	local 	res, dy, alpha, y1, A, y2, gtheta, gdelta, ggamma, gtau, 
			pgamma, rx, z, i, gmu, F1t, F2t, F1, F2, F, e1, ecm;
        

    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;
	ecm = eye(_nlr)*y1 - A*y2 - gmu';	
	z = data[.,_iz]~(ecm');
	z = z';
	e1 = 1|0|0;
	i = 1;
    F2t = -(z[.,i].*.eye(4));
	_nobs = rx;
	F1t = 	(y2[.,i]*alpha')|
			(alpha')|
			(gx2(ggamma,gtau,i)*alpha')|
			(dgx2dgamma(ggamma,gtau,i)*gdelta'alpha')|
			(dgx2dtau(ggamma,gtau,i)*gdelta'alpha')|
			i.*alpha'./rx|
			i.*i.*alpha'./rx; 
	F = F2t[_nln*(_nln-1)*_nlp+1:rows(F2t),.]|F1t|F2t[1:_nln*(_nln-1)*_nlp,.];
	i = 2;
	do while i<=_nobs;
		F1t = 	(y2[.,i]*alpha')|
				(alpha')|
				(gx2(ggamma,gtau,i)*alpha')|
				(dgx2dgamma(ggamma,gtau,i)*gdelta'alpha')|
				(dgx2dtau(ggamma,gtau,i)*gdelta'alpha')|
				i.*alpha'./rx|
				i.*i.*alpha'./rx;
    	F2t = -(z[.,i].*.eye(4));
		F = F~(F2t[_nln*(_nln-1)*_nlp+1:rows(F2t),.]|F1t|F2t[1:_nln*(_nln-1)*_nlp,.]);
		i=i+1;
	endo;
"F" dims(F);
	retp(F');
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);
/*****************************************************************************

    E s t i m a t i o n

*****************************************************************************/
clearg mlagnr;
load mlagnr;
/*
 _gpar0 = mlagnr[1:15]|(0.07|0.07|0.07)|mlagnr[16:71];
*/
_gpar0 = mlagnr;
/*
_gpar0[13:15] = ones(3,1)./_gpar0[13:15];
_gpar0 = mlagnr[1:19 23:rows(mlagnr)];
_gpar0[16:18] = _gpar0[16:18]*10; 
*/
/*
load parfiml[] = parfiml.txt;
_gpar0 = parfiml;
	some reasonable initial values
_gpar0[13:15] = ones(3,1);
_gpar0[16] = 0.04;
_gpar0[19] = -0.04;
_gpar0[22] = 90;
_gpar0[23] = 0.44;
*/

clearg hessu, gradi, cova, retcode;
/*
restrictions in the transition function
_cml_Active = ones(rows(_gpar0),1);
/*
*/
_cml_Active[20] = 0;
_cml_Active[21] = 0;
/*
_cml_Active[23] = 0;
*/
_gpar0[20] = 0;
_gpar0[21] = 0;
/*
_gpar0[23] = 0.46;
*/
*/
/*
_cml_Active[15:25] = zeros(11,1);
*/
/*
	to compute std err we set up linear constraints of the 
	troublesome parameters.
_eka = 20; 
_toka = 21;
_cml_B = zeros(rows(_gpar0),1);
_cml_B[_eka:_toka] = _gpar0[_eka:_toka];
_cml_B[_eka:_toka] = _gpar0[_eka:_toka];
_cml_A = zeros(rows(_gpar0),rows(_gpar0));
_cml_A[_eka:_toka,_eka:_toka] = eye(_toka-_eka+1);
*/

/* or alternatively 
_cml_B = _gpar0[_eka];
_cml_A = zeros(1,rows(_gpar0)); _cml_A[_eka] = 1;
*/

_cml_Diagnostics = 2;
_cml_Options = { bfgs stepbt info central screen };
_cmlshess=1;     /* 1= compute hessian in the beginning; 0 = I as hessian */
_cml_MaxIters = 1000000;
_cml_DirTol = 1e-8;
_cml_UserSearch = 1;
/* definitions related to analytical gradient
_cml_GradProc = &analgrad;
_cml_GradCheckTol = 1e-3;					  */
/*_cml_UserGrad = &argradre;   hopefully better gradient procedure */
__title = "ML estimation of continuous structural changes in intercept terms";
__row=0;
clearg mlagnr, hessu, gradi, cova, retcode;
output file=^ofile on;
{mlagnr,hessu,gradi,cova,retcode}=cmlprt(cml(data,seqa(1,1,cols(data)),&suf,_gpar0));
/* call DisplayResults(Cova, mlagnr, _cml_ParNames, rows(data)); */
output file=^ofile off;
clearg hh; hh = _cml_finalhess;
save mlagnr;
save cova;
format /m1 /rz 15,5;
/**********************************************************************

	likelihood profile

_cml_Select = {16, 17, 18, 19, 22, 23};
{mlagnr,hessu,gradi,cova,retcode}=cmlProfile(data,seqa(1,1,cols(data)),&suf,mlagnr);
end;
**********************************************************************/
/**********************************************************************

	g r a p h i n g  transition function and cointegration vectors

**********************************************************************/
clearg gtheta, gdelta, ggamma, gtau, gmu, A, y1, y2
alpha  = reshape(mlagnr[_ialpha],_nln,_nlr);
A      = mlagnr[_iA];
gtheta = mlagnr[_itheta];
gdelta = mlagnr[_idelta];
ggamma = mlagnr[_iggamma];
gtau   = mlagnr[_itau];
pgamma = reshape(mlagnr[_ipgamma],_nln,_nln*(_nlp-1));
/*
rearranging the data
*/
dy = data[.,_idy]';
y1 = data[.,_iy1]';
y2 = data[.,_iy2]';     /* */
rx = rows(data);
i=1;
do while i<=_nlr;
   	gmu = gx(gtheta[i], gdelta[i], ggamma, gtau, seqa(0,1/rx,rx));
	graphset;
	titext = "Nonlinear Constant in beta Vector: " $+ ftos(i,"%*.*lf",1,0);
	prtstr = "-C=1 -CS=F -CM=0.1,0.1,0.1,0.1 -CF=gmu" $+ ftos(i,"%*.*lf",1,0) $+ ".ps";
	_pscreen = 1; /* 0 do not display graph; 1 do */
	_pltype  = 6;
	_pcolor  = 0;
	graphprt(prtstr);
	title(titext);
	ts(crmiss(_nlp,1)|gmu,80,12);
	i=i+1;
endo;
/* computing cointegration vector */
clearg gmu;
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';
/* save gmu; */
ecm = eye(_nlr)*y1 - mlagnr[_iA]*y2 - gmu;
ecm = ecm';
clearg res;
res = residuals(mlagnr,data);
/*
res = dy  - pgamma * z';
res = dy - alpha*(ones(_nlr,1)*y1 - A'y2 - gmu);
res = alpha*(ones(_nlr,1)*y1 - A'y2 - gmu);
*/
i=1;
do while i<=_nlr;
   	gmu = ecm[.,i];
	graphset;
	titext = "Cointegration Vector: " $+ ftos(i,"%*.*lf",1,0);
	prtstr = "-C=1 -CS=F -CM=0.1,0.1,0.1,0.1 -CF=ecm" $+ ftos(i,"%*.*lf",1,0) $+ ".ps";
	_pscreen = 1; /* 0 do not display graph; 1 do */
	_pltype  = 6;
	_pcolor  = 0;
	ylabel("%/100");
	graphprt(prtstr);
	title(titext);
	ts(crmiss(_nlp,1)|gmu,80,12);
	i=i+1;
endo;
/* call saved(ecm,"ecm","ecm");  Saving the ecm term for later analysis
print seqa(1,1,rows(ecm))~ecm; */
/**********************************************************************

	residual diagnostics

**********************************************************************/
/* computing parameter covariance matrix */
clearg resmom, invresmom, covat, mlagnrcova, Ft, i, F;
_nobs = rx;
F = analgrad(mlagnr,data);
resmom = moment(res,0)/_nobs;
invresmom = inv(resmom);
i = 1;
Ft = F[1:4,.];
covat = Ft'invresmom*Ft;
mlagnrcova = covat;
i=2;
do while i<=_nobs;
	Ft = F[_nln*(i-1)+1:_nln*(i-1)+4,.];
	covat = Ft'invresmom*Ft;
	mlagnrcova = mlagnrcova+covat;
	i = i+1;
endo;
/*
	delete zero delta's
mlagnrcova = mlagnrcova[1:19 22:rows(mlagnrcova),1:19 22:rows(mlagnrcova)];
mlagnr = mlagnr[1:19 22:rows(mlagnr)];
_cml_parnames = _cml_parnames[1:19 22:rows(_cml_parnames)];
*/
mlagnrcova = inv(mlagnrcova);
/* then outputting */
output file=^ofile on;
call DisplayResults(mlagnrCova, mlagnr, _cml_ParNames, rows(data));
call DResDiag(res, rows(yt), _nlp*_nln+_nlr, _nln, labs);
output file=^ofile off;
/*
	computation of LM test for further nonlinearities
*/
clearg Fstar, Ftstar, _score, Fcova, mlagnrcova, i, covat, _scoret, lmstat, df;
Fstar = lmtest(mlagnr,data);
Ftstar = Fstar[1:4,.];
covat = Ftstar'invresmom*Ftstar;
mlagnrcova = covat;
i = 1;
_scoret = Ftstar'invresmom*res[i,.]';
_score = _scoret;
i=2;
do while i<=_nobs;
	Ftstar = Fstar[_nln*(i-1)+1:_nln*(i-1)+4,.];
	covat = Ftstar'invresmom*Ftstar;
	mlagnrcova = mlagnrcova+covat;
	_scoret = Ftstar'invresmom*res[i,.]';
	_score = _score + _scoret;
	i = i+1;
endo;
mlagnrcova = inv(mlagnrcova);
format /rz 12,5;
lmstat = _score'mlagnrcova*_score;
df = _nlr;
df = _nlr*2;
output file=^ofile on;
print;
format /m1 /rd 2,0;
print "Saikkonen's nonlinearity test: chi^2("  df;;
format /m1 /rd 12,6;
print ")= " lmstat "    (" cdfchic(lmstat,df) ")";
print;
output file=^ofile off;
clearg i, tofile;
graphset;
i = 1;
do while i<=cols(res);
	tofile = labs[i] $+ ".ps";
	title(labs[i]);
	resgraph(res[.,i],80,_nln+1,12,13,tofile);
	i = i+1;
endo;
hh = _cml_finalhess;
end;