library maxlik;
#include maxlik.ext;
maxset;


/* 
**  Data - Klein's model I
**
**  The data set actually to be analyzed is computed:  identities and lagged
**  variables are are computed and appended.
*/


/*
**  klein.asc    1     2    3      4     5     6     7      8    9    10
**             Year    C    P     W_p    I     X     W_g   G     T     K
*/

load data0[22,10] = klein.asc;


/*   Data set to be analyzed:
**
**                  1      2       3      4      5      6
**  Endogenous      C      I      W_p     X      P      K
**
**                  7     8     9   10    11     12     13
**  Exogenous      W_g    T     G    A    X_1    P_1    K_1  
*/

data = data0[.,2]~data0[.,5]~data0[.,4]~data0[.,6]~data0[.,3]~data0[.,10]~
                data0[.,7]~data0[.,9]~data0[.,8]~(1931-data0[.,1])~
                lag1(data0[.,6])~lag1(data0[.,3])~lag(data0[.,10]);

data = packr(data);  /* strip off first obs because of lagged variables */


/* 
** First, specify the column numbers in the data of the
** endogenous and exogenous variables:
*/

_eq_Endo = { 1, 2, 3, 4, 5, 6 };   /* endogenous variables */
_eq_Exo  = { 7, 8, 9, 10, 11, 12, 13 };   /* exogenous variables  */

/*
**  Next are the proc's that insert the parameters into
**  the model arrays:
*/

proc _eq_Gamma(b);   /* coefficient matrix of endogenous on endogenous */
    local g;
    g = eye(6);   /* diagonal fixed to ones */

    g[3,1] = -b[1];   /* free parameters */
    g[5,1] = -b[2];
    g[5,2] = -b[3];
    g[4,3] = -b[4];

    g[1,4] = -1;      /* fixed parameters to enforce identities */
    g[2,4] = -1;
    g[2,6] = -1;
    g[3,5] = 1;
    g[4,5] = -1;

    retp(g);
endp;

proc _eq_Beta(b);    /* coefficient matrix of endogenous on exogenous */
    local be;
    be = zeros(7,6);

    be[1,1] = b[1];  /* Coefficienf of C on W_g same as */
                     /* that of C on W_p */
    be[6,1] = b[5];
    be[6,2] = b[6];
    be[7,2] = b[7];
    be[4,3] = b[8];
    be[5,3] = b[9];

    be[2,5] = -1;     /* fixed parameters to enforce identities */
    be[3,4] = 1;
    be[7,6] = 1;

    retp(be);
endp;

proc _eq_Beta_0(b);    /* constants */
     local c;
     c = zeros(6,1);   /* one for each equation */

     c[1] = b[10];     /* free constants */
     c[2] = b[11];
     c[3] = b[12];
                       /* constants for identities remain zero */
     retp(c);
endp;



/*  
**  Residual Covariance Specification
**
**  In the Klein Model I, the last three equations are identities, i.e.,
**  their residuals are fixed to zero.  This constraint is enforced by
**  providing a matrix of the same size as the residual matrix with 
**  1's where the residual variance or covariance is a free parameter
**  to be estimated and 0's where it is fixed to zero.
**
**  If there are no constraints on the residuals, set _resid_constraints
**  to a scalar 1.  Otherwise it is required to have rows and columns 
**  equal to the number of endogenous variables.
*/

   _resid_constraints = { 1 1 1 0 0 0, 
                          1 1 1 0 0 0,
                          1 1 1 0 0 0,
                          0 0 0 0 0 0,
                          0 0 0 0 0 0,
                          0 0 0 0 0 0 };

/*
**  This global checks for zeros on diagonal of _resid_constraints for
**  the log-likelihood procedure, lpr, below.
*/
   _resid_dz = packr(miss(seqa(1,1,rows(_resid_constraints)).*
                       diag(_resid_constraints),0));


/*
**  parameter labels
*/

_max_ParNames = { ALPHA3, ALPHA1, BETA1, GAMMA1, ALPHA2, BETA2, BETA3, 
                 GAMMA3, GAMMA2, ALPHA0, BETA0, GAMMA0 };


/*
** start values 
*/

start = { .810, .017, .150, .439, .216, .616, -.158, .130, .147, 
           17.8, 17.2, 1.41 };

/*
**  output file
*/

output file = klein.1.out reset;

__output = 0;
_max_Alpha = .01;  /*  99% confidence limits */

{ b1,f1,g1,cov1,ret1 } = maxlik(data,0,&lpr,start);


__title = "Wald confidence limits";
cl1 = maxtlimits(b1,cov1);
call maxclprt(b1,f1,g1,cl1,ret1);


print;
print "Stability Test";

gamma_hat = _eq_gamma(b1);
beta_hat = _eq_beta(b1);

delta = (zeros(3,3) | beta_hat[5 6 7, 1 2 3]) ~ zeros(6,3);

e = eig(inv(gamma_hat)*delta);

print;
print "abs(eigenvalues) " abs(e);



output off;



/**********************************************************
**  Likelihood function
*/


proc lpr(b,z);
    local u, psi, oldt, ldet0, ipsi, _beta, _gamma, _beta_0, id;

    _gamma = _eq_gamma(b);
    _beta = _eq_Beta(b);
    _beta_0 = _eq_Beta_0(b);

    u = z[.,_eq_Endo] * _gamma - z[.,_eq_Exo] * _beta - _beta_0';
    if not scalmiss(_resid_dz);
        psi = _resid_constraints .* (u'u/rows(u));
        psi = psi[_resid_dz,_resid_dz];
        u = u[.,_resid_dz];
     else;
        psi = _resid_constraints .* (u'u/rows(u));
    endif;

    oldt = trapchk(1);
    trap 1,1;
    ipsi = invpd(psi);

    trap oldt,1;
    if scalmiss(ipsi);
        retp(error(0));
    endif;

    ldet0 = ln(detl);

    retp(-0.5 * rows(_eq_Endo) * ln(2*pi) 
            + ln(abs(det(_gamma))) - 0.5*(ldet0 + sumc(((u*ipsi).*u)')));
endp;