new;
/*
**	gauss program for computing Saikkonen's test for normalization
**	Written by Antti Ripatti, 15 Nov 1996 (email: antti.ripatti@helsinki.fi)
** 	See also http://www.helsinki.fi/~ripatti/
**	
**	DISCLAIMER!
**	This code is initially written for own research purposes. 
**	It is then submitted for public, non-commercial use. 
**	There are no performance guarantees.
**	Please acknowledge this code (and its author and the authors of relating
**	paper) if you find it useful and use it extensively in your work.
**
*/
library pgraph;
p       = 3; /* lag length */
ofile   = "norm.out";
r       = 2;
n       = 4;
/* computes trace of a matrix */
proc (1) = tr(x);
    retp(sumc(diag(x)));
endp;
clearg qtile;
qtile = (10.56~12.39~16.39)|
        (22.95~25.47~30.65)|
        (39.08~42.20~48.59)|
        (58.96~62.61~70.22);
/* Some general purpose procedures */        
/* ADJOBS
**
** Purpose:   returns the last nobs rows of the matrix
**
** Format:    y = ADJOBS(x,nobs);
**
** Input:     x     T x K matrix, time serie
**            nobs  scalar, the desired number of rows of x
**
** Output:    y     T-nobs x K matrix, where the first T-nobs rows are deleted
**            if nobs>T then prints error message and returns scalar -1 
**
** Created:   15.4.96 by AR
** Revised:
**
*/
proc (1) = AdjObs(x,nobs);
    local xrows, retvalue;

    xrows = rows(x);
    if xrows < nobs;
        errorlog "Error in calling procedure ADJOBS. Number of rows desired";
        errorlog "exceeds the number of rows available. Function returns -1";
        retvalue = -1;
    else;
        retvalue = x[xrows-nobs+1:xrows,.];
    endif;
    retp(retvalue);
endp;


proc (1) = ao(x,nobs);

    retp(adjobs(x,nobs));
endp;
        
/* DIFFE */
/* Calculates first differeces
** >> Adds a miss. value in place of the first observation
*/
FN diffe(X) = miss(zeros(1,cols(X)),0) | trimr(X,1,0)-trimr(X,0,1);

/*
    loading data file and data labels
*/
clearg dat, datlab;
dat     = loadd("igauss");
datlab  = getname("igauss");
clearg xt, labs, dxt, xt1, zt, i;
xt      = dat[.,4 3 11 8];
labs    = datlab[4 3 11 8];
dXt     = diffe(xt);
xt1     = lagn(xt,1);
i = 2;
zt      = lagn(dxt,1);
do while i<=(p-1);
    zt = zt~lagn(dxt,i);
    i = i+1;
endo;
zt = packr(zt);
clearg nobs;
nobs = rows(zt);
dxt = ao(dxt,nobs);
xt1 = ao(xt1,nobs);
const = ones(nobs,1);
/*
    estimation of eq. (4.4)
    -----------------------
*/
clearg selitt, params;
selitt = const~xt1~zt;
params = inv(selitt'selitt)*selitt'dxt;
/*
    disaggregation of parameter matrix
*/
clearg nu, pii, gamm, epsilon, omega;
nu = params[1,.]';
pii = params[2:5,.]';
gamm = params[6:rows(params),.]';
epsilon = dxt - Selitt*params;
omega = moment(epsilon,0)/(nobs-p);
output file=^ofile reset;
print "FILE: " $ofile "  DATE: " datestr(date()) "  TIME: " timestr(time());
print "by Antti Ripatti, November 15 1996";
print;
print "=============================================================================";
print "            TESTING NORMALIZATION RESTRICTIONS BY SAIKKONEN (1996)";
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: " nobs;
print;
output file=^ofile off;
/*
    Testing normalization restrictions using unrestricted estimator of the constant
*/    
proc (0)=TestNorm(c, pii, omega);
    local alpha_c, omega_z, alpha_c_o, selitt, lambda_W, lambda_LM, lambda_LR, dfkorjaus;

    alpha_c = pii*c;
    alpha_c_o = null(alpha_c'); 
    selitt  = ones(nobs,1)~zt;
    omega_z = moment((dxt - selitt*inv(selitt'selitt)*selitt'dxt),0)/(nobs-p);
    dfkorjaus = (nobs-P-p*N)/(nobs-p);
    /*
        computing the test statistics
    */
    lambda_W  = dfkorjaus*(nobs-p)*tr(inv(alpha_c_o'omega*alpha_c_o)*(alpha_c_o'omega_z*alpha_c_o - alpha_c_o'omega*alpha_c_o));
    lambda_LM = dfkorjaus*(nobs-p)*tr(inv(alpha_c_o'omega_z*alpha_c_o)*(alpha_c_o'omega_z*alpha_c_o - alpha_c_o'omega*alpha_c_o));
    lambda_LR = dfkorjaus*(nobs-p)*(ln(det(alpha_c_o'omega_z*alpha_c_o))-ln(det(alpha_c_o'omega*alpha_c_o)));
    /*
        printout of the results
    */
    output file=^ofile on;
    print;
    format /m1 /rd 12,6;
    print "Unrestricted constant_____________________________";
    print "Lagrange multiplier test statistic : " lambda_LM ;
    print "Likelihood ratio test statistic    : " lambda_LR ;
    print "Wald test statistic                : " lambda_W ;
    print "Quantiles from Table 1 by Reinsel and Ahn (1992):";
    print "      50 %       95 %       99 % ";
    print "      9.41      17.97      22.79";
    print;
    print;
    output file=^ofile off;
endp;
/*
    testing normalization restrictions using GLS estimator of the constant
*/
clearg Xlags;
Xlags = lagn(xt,1);
i = 2;
do while i<=p;
    Xlags = Xlags~lagn(xt,i);
    i  = i+1;
endo;
Xlags = packr(Xlags);
Xlags = ao(Xlags,nobs);
xt    = ao(xt, nobs);
clearg selitt;
Selitt = ones(rows(Xlags),1)~Xlags;
Btilde = inv(selitt'selitt)*selitt'xt;
Btilde = Btilde[2:rows(Btilde),.];
clearg dt, Ct, tmpB, tmpXlags, mtilde, tmp1, tmp2, i, invomega;
dt      = xt - Xlags*Btilde;
Ct = eye(n);
i = 2;
do while i<=nobs;
    if i<=p+1;
        Ct = Ct|(Ct[n*(i-2)+1:n*(i-1),.] - Btilde[n*(i-2)+1:n*(i-1),.]');
    else;
        Ct = Ct|(Ct[n*(i-2)+1:n*(i-1),.]);
    endif;
    i = i+1;
endo;
invomega = inv(omega);
i = 1;
tmp1 = zeros(n,n); tmp2 = zeros(n,1);
do while i < nobs;
    tmp1 = tmp1 + Ct[n*(i-1)+1:n*i,.]'invomega*Ct[n*(i-1)+1:n*i,.];
    tmp2 = tmp2 + Ct[n*(i-1)+1:n*i,.]'invomega*dt[i,.]';
    i = i + 1;
endo;
mtilde = inv(tmp1)*tmp2;

output file=^ofile on;
format /m1 /rd 12,6;
print;
print "GLS estimate of the constant ~m:" mtilde;
print;
output file=^ofile off;

selitt = (xt1-ones(nobs,n)*mtilde)~zt;
params = inv(selitt'selitt)*selitt'dxt;
clearg piihat, gammhat, epsilonhat, omegahat;
nu = params[1,.]';
piihat = params[1:4,.]';
gammhat = params[5:rows(params),.]';
epsilonhat = dxt - Selitt*params;
omegahat = moment(epsilonhat,0)/(nobs-p);

proc (0)=TestNormGLS(c, piihat, omegahat);
local alpha_c, omega_z, alpha_c_o, selitt, lambda_W, lambda_LM, lambda_LR, dfkorjaus;
    
    alpha_c = piihat*c;
    alpha_c_o = null(alpha_c'); 
    selitt  = zt;
    omega_z = moment((dxt - selitt*inv(selitt'selitt)*selitt'dxt),0)/(nobs-p);
    dfkorjaus = (nobs-P-p*N)/(nobs-p);
    /*
        computing the test statistics
    */
    lambda_W  = dfkorjaus*(nobs-p)*tr(inv(alpha_c_o'omegahat*alpha_c_o)*(alpha_c_o'omega_z*alpha_c_o - alpha_c_o'omegahat*alpha_c_o));
    lambda_LM = dfkorjaus*(nobs-p)*tr(inv(alpha_c_o'omega_z*alpha_c_o)*(alpha_c_o'omega_z*alpha_c_o - alpha_c_o'omegahat*alpha_c_o));
    lambda_LR = dfkorjaus*(nobs-p)*(ln(det(alpha_c_o'omega_z*alpha_c_o))-ln(det(alpha_c_o'omegahat*alpha_c_o)));
    /*
        printout of the results
    */
    output file=^ofile on;
    format /m1 /rd 12,6;
    print "Constant estimated using GLS______________________";
    print "Lagrange multiplier test statistic : " lambda_LM ;
    print "Likelihood ratio test statistic    : " lambda_LR ;
    print "Wald test statistic                : " lambda_W ;
    print "Quantiles from Table 15.1 by Johansen (1995):";
    print "      50 %       95 %       99 % ";
    print "      5.47      12.21      16.16";
    print;
    print;
    output file=^ofile off;
endp;

/*
    normalization matrix c
*/
clearg c;
c = zeros(n,r);
c[1,1] = 1;
i = 1;
do while i<= n;
    j = 1;
    c[i,1] = 1;
    do while j<= n;
        c[j,2] = 1;
        output file=^ofile on;
        print "The normalization matrix c";
        format /m1 /rD 4,0;
        print c;
        output file=^ofile off;
        TestNorm(c,pii,omega);
        TestNormGLS(c,piihat,omegahat);
        c[j,2] = 0;
        j = j+1;
    endo;
    c[i,1] = 0;
    i = i+1;
endo;
output file=^ofile on;
print "Pervert c matrix";
c = zeros(n,r);
output file=^ofile off;
TestNorm(c,pii,omega);
TestNormGLS(c,piihat,omegahat);
end;