/* THIS IS FREESHARE GAUSS CODE: THE USUAL CAVEATS APPLY                                        **
** RESEARCHERS USING THIS CODE SHOULD QUOTE THE ORIGINAL PAPER:                                 **
** HARRIS, M.N. and MATYAS, L. (2004),                                                          **
** "A Comparative Analysis of Different IV and GMM Estimators of Dynamic Panel Data Models",    **
** International Statistical Review, vol. 72, 3, pp. 397-408.                                   **
*/
 

library co;

            /* DPD.PRC: DYNAMIC PANEL DATA ESTIMATORS */

/*
** Estimates in levels: OLS, Within, FGLS,
** Estimates in first differences: OLS, Anderson-Hsiao, Arellano,
**                                 Arellano-Bond
** Estimates unspecified transformation: Wansbeek-Bekker
** Estimates in levels: non-linear GMM
**
** Takes form:
**
**    {ols,R2_ols,NTols,
**    within,R2_wn,
**    fgls,  R2_fgls,
**    dols,  R2_dols,NTdols
**    aha,   R2_aha, ahaid, AR2aha,NTah,
**    ahb,   R2_ahb, ahbid, AR2ahb,
**    ahc,   R2_ahc, ahcid, AR2ahc,
**    ara,   R2_ara, araid, AR2ara,
**    arb,   R2_arb, arbid, AR2arb,
**    arc,   R2_arc, arcid, AR2arc,
**    aba,   R2_aba, abaid, AR2aba,
**    abc,   R2_abc, abcid, AR2abc,
**    bb,    R2_bb, bbid, AR2bb,
**    wba,   R2_wba, wbaid,
**    wbb,   R2_wbb, wbbid,
**    gmmi,  R2_gmmi,pwaldi,
**    gmm,   R2_gmm, pwald}
**
**        = dpd("data_set",ivtype,timev,xvars,endvars,gmm_con,gmm6,gmm7,gmm_se);
**
** Input:   dataset - string, first column is y
**                      (DON'T SUPPLY LAGGED DEPENDENT VARIABLE)
**                      second column must be id vector and 3rd time vector.
**                      Rest of data set has K* columns consisting of both
**                      explanatory and instrumental variables.
**          note    - id vector, MUST BE CONTIGUOUS, but can start at any #
**                  - time vector, MUST BE CONTIGUOUS, but can start at any #
**                  - panel must be BALANCED
**                  - min T = 4
**          type    - "linear" or "non" - linear gives standard panel
**                    estimators and linear IV's;
**                    "wb" gives standard panel and Wanbeek-Bekker estimators
**                    "gmm" gives standard panel and GMM estimators
**          timev   - 1 by K* vector - has a 1 corresponding
**                    to time-varying X's (and Z's) and zeros elsewhere
**                    (first element corresponding
**                    to constant is therefore 0, and if no
**                    explanatory vars timev = 0 ie scalar)
**          xvars   - K* row vector - has a 1 corresponding to
**                    explanatory variables and zeros otherwise (first element
**                    should be 1, corresponding to the constant)
**          endvars - K* row vector - has a 1 corresponding to endogenous
**                    explanatory variables and zeros otherwise (thus the
**                    instruments will be everything in XZ not identified as
**                    endogenous NOTE THESE IMPLIED INSTRUMENTS ARE FOR
**                    ENDOGENOUS X'S NOT THE LAGGED DEPENDENT VARIABLE
**          gmm_con - row vector takes 1 if condition is to be used. Order:
**                    1. E(vit = a_i + u_it)        = 0
**                    2: V(v_it)                    = sig_a^2 + sig_u^2
**                    3: E(vit)(v_is)               = sig_a^2
**                    4: E([(v_it)-(v_it-1)]y_it-l) = 0 l>=2
**                    5: E(v_it)(x_it)              = 0
**                      Note this uses 1 - endvars
**                    6: E(v_it)(x_i0)              = 0
**                    7: E(v_it)(x_is)              = 0
**          gmm6    - additional info for gmm_con: defines which vars to use
**                    GMM condition 6 with
**          gmm7    - additional info for gmm_con: defines which vars to use
**                    GMM condition 7 with
**          gmm_se  - 1 = numerical standard errors (recommended),
**                    for optimal gmmm only
**                    0 = analytical standard errors
**
**
** Output:  -   all outputs give simple panel estimates, plus the selected
**              consistent estimators in form: coefficients~se's~t-stats,
**              Adj Rsqr. The GMM estimators also include a Sargan test for over-identifying 
**              restrictions: test~df's~p-value
**              NT is also reported, OLS, FGLS, WITHIN, WB, GMM, all have same sample size, 
**              as do DOLS, AR and AB
**
*/

cls;
trap 1;

proc 59 = dpd(dset,ivtype,timev,xvars,endvars,gmm_con,gmm6,gmm7,gmm_se);
local infile1, fh, varsx, idvec, tvec, tvec1, z, x, k, n, tmax, xcheck;
local xrep, t_i, col, ivs, exog, exogtv, xvarstv;
local wn, wnabv, bn, bnabv;
local irep, ylag, dy, dylag, e1, e2, xlag, dx, dxz, dxlag, dylag2, ylag2;
local sel, ytemp, xtemp, abrep, abtrep, abicols, abadd, abz, y_i, trep;
local xolstmp, ols, uols, seols, tols, sig2_ols, Rols, R2_ols, wnxtmp1;
local wny, wnylag, wnx, wnxall, uwn, within, sewn, twn, sig2_wn, Rwn, R2_wn;
local wnytmp, wnxtmp, selwn, wny1tmp;
local ybar, xbar, ylagbar, uhat, sigv2, sigmu2, theta, yfgls, xfgls, ylagfgls;
local xallfgls,xfglstmp, fgls,sefgls,tfgls,R2_fgls,ufgls,sig2_gls,Rfgls,idvec1;
local xbartmp, zbn, pzbn, bn1, ubn1, sigu, talphu, omeginv,lamda, p;
local dxalltmp, dxall, dols, sedols, tdols, R2_dols, Rdols, udols, sig2_dls;
local xahtmp, xah, zahatmp,zaha, aha, uaha, sig2_aha, seaha, taha, Raha, R2_aha;
local phi1, phi, phirep, ahb, zahb, pzahb, seahb, tahb, ahc, seahc, tahc;
local pzaha, zahbtmp;
local uahb, sig2_ahb, Rahb, R2_ahb, colsah, rowsah, zhatahc, z1i, uahbi;
local pzahc, uahc, sig2_ahc, Rahc, R2_ahc;
local xar, zara, ara, uara, sig2_ara, seara, tara, Rara, R2_ara;
local phii, phii1, zarb, pzarb;
local arb, uarb, sig2_arb, searb, tarb, Rarb, R2_arb;
local zhatarc, pzarc, arc, uarc, sig2_arc, searc, tarc, Rarc, R2_arc;
local colsar, rowsar, uarbi;
local abz1, pzaba, aba, seaba, taba, R2_aba, uaba, sig2_aba, Raba;
local colsab1, rowsab1, zhatabc, irepst, irepfn, uabai;
local pzabc, abc, seabc, tabc, R2_abc, uabc, sig2_abc, Rabc;
local ai, arep, aivec, ct, crep, eyerep, ctemp1, ctemp2, ctemp3, ctemp4;
local ctemp5, c1, rt, b, ywb;
local ahat, f, q, qw, g, retcode, aia, aawb, zwba, pzwba, wba, uwba;
local mtype, zhataba, yab, xab, tvec2, dx1tmp, tnew, dynew, dyab;
local colswb, rowswb, xwbtmp, zhatwba, uwbai, sewba, twba, Rwba, R2_wba;
local wbb, sewbb, twbb, R2_wbb, aib, abwb, zwbb, pzwbb, uwbb, zhatwbb;
local uwbbi, Rwbb;
local gmmi, segmmi, tgmmi, ugmmi, Rgmmi, R2_gmmi, sig2_gmi;
local gmm, segmm, tgmm, ugmm, Rgmm, R2_gmm;
local xartmp, zaratmp, pzara, xarbtmp, xarb, zarbtmp;
local tvx, tvexogx, tvexogx1, tvx1, dz, tvz, zlag, dzlag, ztemp;
local start, vit, vitmat, vitmatt, m5itmp, m13i, m14i, yt0, xt0;
local m1a2i, sel_con, m3i, m5irep, vtreplf, vtreprg, m5i;
local m6itmp, m6i, colsm7, colrep, xtmp, m7i, colsm11, xt0tmp, m11i;
local x13, colsm13, rep5, rep6, mi, checkall, start2;
local dum, avy, avytmp, avx, avdx, avxtmp, avu;
local m1a2g, m3g, capM, capG, vargmmi, m5g, m6g, m7g, m11g, m13g;
local avzgmtmp, avzgmm, ztmp;
local jgmm, dfgmm, pwald, pwaldi;
local m1a2c, m3c, m5c, m6c, m7c, m11c, m13c, m14c, whichgmm, xt0m11;
local vargmm, sig2_gmm;
local xm7, avzm7, cdo, varchk, nofgls, overidi, overid, qi;
local infomat, check;
local delta_v, ymat, m6rep, deltav_i;
local ahasarg, ahadfsrg, ahapsarg, ahbsarg, ahbdfsrg, ahbpsarg, ahcsarg, ahcdfsrg, ahcpsarg;
local arasarg, aradfsrg, arapsarg, arbsarg, arbdfsrg, arbpsarg, arcsarg, arcdfsrg, arcpsarg;
local abasarg, abadfsrg, abapsarg, abcsarg, abcdfsrg, abcpsarg;
local wbasarg, wbbsarg, wbadfsrg, wbbdfsrg, wbapsarg, wbbpsarg;
local v, utmp;
local ahaid, ahbid, ahcid, araid, arbid, arcid, abaid, abcid, wbaid, wbbid;
local ustar, u2, utop, AR2ara, AR2arb, AR2arc, sel1, selu2, part1, part2;
local AR2aha, AR2ahb, AR2ahc, AR2aba, AR2abc, ubot;
local varaha, varahb, varahc, varara, vararb, vararc, varaba, varabc;
local NTols, NTdols,NTah;
local bb, R2_bb, bbid, AR2bb, sebb, tbb, ytild, xbb, xtild, bbz, pzbb, z_l, bbsarg, bbdfsrg, bbpsarg;
local bbrep, ubb, colsbb, rowsbb, ubbi, zhatbb, varbb, Rbb, idbb, pzhatbb;

clearg xols, yols, xwb, y, zwbatmp, xfull, e, xgmm, xt0, what;


cdo         = cdir(0);
load path   = ^cdo;
save path   = ^cdo;


/* READ IN FULL DATA SET */

infile1 = dset;
open fh = ^infile1 varindxi; 
varsx   = getname(infile1);
varsx   = varsx[2:rows(varsx)];

z   = readr(fh,rowsf(fh));
fh  = close(fh);


ivs     = (1-xvars)';
exog    = (1-endvars)';
y       = z[.,1];
idvec   = z[.,2];

/* CREAT CONTIGUOUS ID VECTOR */

idvec1   = 1;

for irep (2,rows(z),1);

    if idvec[irep] .eq idvec[irep-1];
        idvec1 = idvec1|idvec1[irep-1];
    else;
        idvec1 = idvec1|(idvec1[irep-1]+1);
    endif;               

endfor;

idvec = idvec1;
clear idvec1;

tvec    = z[.,3];
x       = z[.,4:cols(z)];
k       = cols(x);


/* CHECK DIMENSIONS */

if cols(timev) ne k;
    cls;"ERROR IN DIMENSIONS OF TIME VARYING INDICATOR - program aborted";;stop;
endif;

if cols(xvars) ne k;
    cls;"ERROR IN DIMENSIONS OF X VARIABLE INDICATOR - program aborted";;stop;
endif;

if cols(endvars) ne k;
    cls;"ERROR IN DIMENSIONS OF Z VARIABLE INDICATOR - program aborted";;stop;
endif;


tvx     = timev .* xvars;
tvexogx = tvx .* exog';
tvexogx1= selif(tvexogx',tvx')';   /* indicator for time varying
                                       exogenous variables, consistent with
                                       dimensions of delta X */

gmm6 = gmm6 .* tvexogx;            /* remove time invariant x's */
gmm7 = gmm7 .* tvexogx;            /* remove time invariant x's */
tvz     = timev .* ivs';


if sumc(tvexogx') ge 1;
    xfull   = selif(x',tvexogx')';
else;
    xfull = 0;
endif;



            /* DEFINE N AND T ETC */

tvec    = tvec - minc(tvec) + 1;
tmax    = maxc(tvec);
idvec   = idvec - minc(idvec) + 1;
n       = maxc(idvec);
yt0     = selif(y,tvec .eq minc(tvec));
tvec1   = delif(tvec,tvec .eq minc(tvec));



save n, tmax, y, xfull, timev, xvars, ivs, gmm_con, yt0, gmm6, gmm7, tvexogx;



/* CHECK FOR SUFFICIENT VARIATION IN X & Z */

if k gt 1;
    xcheck  = stdc(x[.,2:k]);
    xcheck = (xcheck .le 0.1);

    for xrep (1,cols(x)-1,1);
            if xcheck[xrep] eq 1;
                    "Potentially Insufficient Variation in:";$varsx[xrep+3];
            endif;
    endfor;
endif;



/* GENERATE DELTA VARIABLES AND LAGGED DEPENDENT VARIABLE &
        DELETE OBS IF LESS THAN 4 OBS */

ylag    	= 0;
dy      	= 0;
dylag   	= 0;
e       	= 0;
e1      	= 0;
e2      	= 0;


if sumc(tvx') ge 1;
    xlag    = zeros(1,sumc(tvx'));
    dx      = zeros(1,sumc(tvx'));
    dxlag   = zeros(1,sumc(tvx'));
    wnx     = zeros(1,sumc(tvx'));
endif;

if sumc(tvz') ge 1;
    zlag    = zeros(1,sumc(tvz'));
    dz      = zeros(1,sumc(tvz'));
    dzlag   = zeros(1,sumc(tvz'));
endif;


dylag2  = 0;
ylag2   = 0;
wny     = 0;
wnylag  = 0;
ybar    = 0;
ylagbar = 0;
xbar    = zeros(1,k);


for irep (1,n,1);
    sel     = idvec .eq irep;
    ytemp   = selif(y,sel);
        if rows(sel) lt 4;
            y   = delif(y,sel);
            x   = delif(x,sel);
            idvec   = delif(idvec,sel);
            idvec   = ((idvec - 1) .* (idvec .gt irep)) +
                        (idvec .* (idvec lt irep));
            tvec    = delif(tvec,sel);
            goto nexti;
        endif;


    if sumc(tvx') ge 1;
        xtemp   = selif(x',tvx');
        xtemp   = selif(xtemp',sel);
    endif;


    if sumc(tvz') ge 1;
        ztemp   = selif(x',tvz');
        ztemp   = selif(ztemp',sel);
    endif;

    ylag    = ylag|-999|(ytemp[1:rows(ytemp)-1]);
    ylag2   = ylag2|(ones(2,1).*-999)|ytemp[1:rows(ytemp)-2];

    if cols(timev) ge 2;

        if sumc(tvx') ge 1;

           xlag    = xlag|(ones(1,sumc(tvx')).*-999)|xtemp[1:rows(xtemp)-1,.];
           wnxtmp  = (selif(x',tvx'))';
           wnxtmp  = selif(wnxtmp,sel);
           wnxtmp  = wnxtmp[2:rows(wnxtmp),.];
           wnx     = wnx|(wnxtmp - (meanc(wnxtmp)' .*. ones(sumc(sel)-1,1)));
           dx      = dx|(ones(1,sumc(tvx')).*-999)|
                     (xtemp[2:rows(xtemp),.]- xtemp[1:rows(xtemp)-1,.]);
           dxlag   = dxlag|(ones(2,sumc(tvx')).*-999)|
                     (xtemp[2:sumc(sel)-1,.] - xtemp[1:sumc(sel)-2,.]);
        endif;

        if sumc(tvz') ge 1;

           zlag    = zlag|(ones(1,sumc(tvz')).*-999)|ztemp[1:rows(ztemp)-1,.];
           dz      = dz|(ones(1,sumc(tvz')).*-999)|
                     (ztemp[2:rows(ztemp),.] - ztemp[1:rows(ztemp)-1,.]);
           dzlag   = dzlag|(ones(2,sumc(tvz')).*-999)|
                     (ztemp[2:sumc(sel)-1,.] - ztemp[1:sumc(sel)-2,.]);
        endif;

    endif;
    e       = e|1|zeros(sumc(sel)-1,1);
    e1      = e1|ones(2,1)|zeros(sumc(sel)-2,1);
    e2      = e2|ones(3,1)|zeros(sumc(sel)-3,1);

    wnytmp  = selif(y,sel);
    wny1tmp = wnytmp[1:rows(wnytmp)-1];
    wnytmp  = wnytmp[2:rows(wnytmp)];

    xbartmp = selif(x,sel);
    xbartmp = xbartmp[2:rows(xbartmp),.];
    xbar    = xbar|(meanc(xbartmp)'.*. ones(sumc(sel)-1,1));


    wny     = wny|(wnytmp - meanc(wnytmp));
    wnylag  = wnylag|(wny1tmp - meanc(wny1tmp));

    ybar    = ybar|(meanc(wnytmp).*. ones(sumc(sel)-1,1));
    ylagbar = ylagbar|(meanc(wny1tmp).*.ones(sumc(sel)-1,1));

    dy      = dy|-999|(ytemp[2:rows(ytemp)]-ytemp[1:rows(ytemp)-1]);
    dylag   = dylag|(ones(2,1).*-999) |
                (ytemp[2:sumc(sel)-1] - ytemp[1:sumc(sel)-2]);
    dylag2  = dylag2|(ones(3,1).*-999)|
                (ytemp[2:sumc(sel)-2] - ytemp[1:sumc(sel)-3]);

    nexti:
endfor;


/* DEFINE N AND NT */

n       = maxc(idvec);
ylag    = ylag[2:rows(ylag)];

if cols(timev) ge 2;

    if sumc(tvx') ge 1;
        xlag    = xlag[2:rows(xlag),.];
        dx      = dx[2:rows(dx),.];     @ NOTE DX DOESN'T CONTAIN A CONSTANT @
        dxlag   = dxlag[2:rows(dxlag),.];
        wnx     = wnx[2:rows(wnx),.];
    endif;

    if sumc(tvz') ge 1;
        zlag    = zlag[2:rows(zlag),.];
        dz      = dz[2:rows(dz),.];     @ NOTE DZ DOESN'T CONTAIN A CONSTANT @
        dzlag   = dzlag[2:rows(dzlag),.];
    endif;

endif;

dy      = dy[2:rows(dy)];
e       = e[2:rows(e)];
e1      = e1[2:rows(e1)];
e2      = e2[2:rows(e2)];
dylag   = dylag[2:rows(dylag)];
dylag2  = dylag2[2:rows(dylag2)];
ylag2   = ylag2[2:rows(ylag2)];
wny     = wny[2:rows(wny)];
wnylag  = wnylag[2:rows(wnylag)];
ybar    = ybar[2:rows(ybar)];
xbar    = xbar[2:rows(xbar),.];
ylagbar = ylagbar[2:rows(ylagbar)];

save e, exog;


/* GENERATE INSTRUMENT SET FOR ARELLANO-BOND ESTIMATOR */

/* CALCULATE TOTAL # OF COLUMNS: IF T > 7,
    ONLY MOST RECENT 5 PAST Y'S ARE USED */

abrep   = 1;
abtrep  = 3;
abicols = 1;

abadd   = 2;

do while abtrep lt tmax;
    abicols = abicols + abadd;
    if abadd lt 5;
        abadd = abadd + 1;
    else;
        abadd = 5;
    endif;
    abtrep = abtrep + 1;
endo;

abz = zeros(1,abicols);

irep  = 1;

do while irep le n;

    sel = (idvec .eq irep);
    y_i = selif(y,sel);
    t_i = sumc(sel);

    abrep   = 1;
    abtrep  = 3;
    col     = 1;
    abadd   = 2;
    start   = 1;

    for trep (1,t_i-2,1);

        if trep eq 1;
            abz = abz|(y_i[start]~zeros(1,abicols-1));
            goto nextt;
        endif;

        if trep eq rows(y_i)-2;

            if abadd lt 6;
                abz = abz|zeros(1,col)~y_i[start:abadd]';
            else;
               abz = abz|zeros(1,cols(abz)-5)~y_i[start:start+4]';
            endif;

        else;

            if abadd lt 6;

                abz = abz|zeros(1,col)~
                        y_i[start:abadd]'~zeros(1,abicols-col-abadd-start+1);

            else;

                abz = abz|zeros(1,col)~
                        y_i[start:abadd]'~zeros(1,cols(abz)-col-5);
            
            endif;

        endif;

        if abadd lt 6;
            col     = col + abadd + start - 1;
        else;
            col     = col + 5;
        endif;        

        abadd   = abadd + 1;
 
        if abadd ge 6;
            start = start + 1;
        endif;
     
    nextt:
    endfor;

irep=irep+1;
endo;



abz = abz[2:rows(abz),.];   @ NOTE ABZ DELETES FIRST 2 OBS PER INDIVID @


        /* DEFINE EXOGENOUS VARIABLES FOR GMM: CAN BE TIME INVARIANT */

xt0     = selif(x,e);
xgmm    = delif(x,e);


save xt0;


       /* ESTIMATE USUAL LAGGED DEPENDENT VARIABLE MODEL */

            /* OLS */

xolstmp = (selif(x',xvars'))';


if sumc(xvars') gt 1;
    xols    = delif(xolstmp[.,1]~ylag~xolstmp[.,2:cols(xolstmp)],e);
else;
    xols    = delif(ones(rows(x),1)~ylag,e);
endif;

if ismiss(xols) eq 1;
	cls;
	"WARNING - MISSING VALUES IN LEVELS X'S (SIMPLE OLS ESTIMATOR)";
endif;


yols    = delif(y,e);
NTols   = rows(yols);                    
ols     = inv(xols'*xols)*xols'*yols;
if ismiss(inv(xols'*xols)) eq 1;
    ols     = zeros(cols(xols),1);
    seols   = ols;
    tols    = ols;
    R2_ols  = -999;
    nofgls  = 1;
else;
    nofgls      = 0;
    uols        = delif(y,e) - (xols*ols);
    sig2_ols    = (uols'*uols)/(rows(uols)-rows(ols));
    seols       = diag(sig2_ols*inv(xols'*xols))^0.5;
    tols        = abs(ols ./ seols);
    Rols        = 1 - ((uols'uols) / sumc((yols - meanc(yols))^2));
    R2_ols      = 1 - (((rows(uols)-1) / (rows(uols)-rows(ols)))
                    * (1 - Rols));
endif;


/* FIXED EFFECTS ESTIMATOR (WITHIN) */

if cols(timev) eq 1;

    wnxall  = wnylag;

else;

    if sumc(timev[2:cols(timev)]' .* xvars[2:cols(xvars)]') ge 1;
        wnxall  = wnylag~wnx;
    else;
        wnxall  = wnylag;
    endif;

endif;

if ismiss(wnxall) eq 1;
	print;
	"WARNING - MISSING VALUES IN LEVELS X'S (WITHIN)";
endif;


if ismiss(inv(wnxall'*wnxall)) eq 1;
    within  = zeros(cols(wnxall),1);
    sewn    = within;
    twn     = within;
    R2_wn   = -999;
else;
    within  = inv(wnxall'*wnxall)*wnxall'*wny;
    uwn         = wny - (wnxall*within);
    sig2_wn     = (uwn'*uwn)/(rows(uwn)-n-rows(ols)+1);
    sewn        = diag(sig2_wn*inv(wnxall'*wnxall))^0.5;
    twn         = abs(within ./ sewn);
    Rwn         = 1 - ((uwn'uwn) / sumc((yols - meanc(yols))^2));
    R2_wn       = 1 - (((rows(uwn)-1) / (rows(uwn)-rows(ols)-n+1))
                    * (1 - Rwn));
endif;

clear wnxall, uwn, sig2_wn, Rwn, xcheck;



/* FGLS ESTIMATOR BASED UPON OLS RESIDS */

if nofgls eq 0;

    idvec1  = delif(idvec,e);
    uhat = 0;
    i = 1;
    do while i le n;
        sel = idvec1 .eq i;
        uhat= uhat|((sumc(selif(uols,sel))) ./ sumc(sel)).^2;
    i=i+1;
    endo;

    uhat = uhat[2:rows(uhat)];

    sigv2   = ((sumc(uols.^2)./ (rows(uols)-rows(ols))) -
                (sumc(uhat) ./ n - rows(ols))) .* (tmax/(tmax-1));
    sigmu2  = (sumc(uols.^2) - sigv2) ./ (rows(yols) - rows(ols));

    theta   = sigv2 / (sigv2 + (tmax*sigmu2));


    if cols(xbar) gt 1;
        xbar    = (selif(xbar',xvars'))';
    endif;


    yfgls       = yols - ((1-sqrt(theta)).*ybar);
    xfgls       = (delif(xolstmp,e)) - ((1-sqrt(theta)).*xbar);
    ylagfgls    = delif(ylag,e) - ((1-sqrt(theta)).*ylagbar);


    if sumc(xvars') gt 1;
        xallfgls    = ylagfgls~xfgls;
    else;
        xallfgls    = ylagfgls~xfgls[.,1];
    endif;

	if ismiss(xallfgls) eq 1;
		print;
		"WARNING - MISSING VALUES IN LEVELS X'S (FGLS)";			
	endif;

    fgls = inv(xallfgls'xallfgls)*xallfgls'yfgls;

    if ismiss(inv(xallfgls'xallfgls)) eq 1;
        fgls    = zeros(cols(wnxall),1);
        sefgls  = fgls;
        tfgls   = fgls;
        R2_fgls = -999;
    else;
        ufgls       = yfgls - (xallfgls*fgls);
        sig2_gls    = (ufgls'*ufgls)/(rows(ufgls)-rows(fgls));
        sefgls      = abs(diag(sig2_gls*inv(xallfgls'*xallfgls)))^0.5;
        tfgls       = abs(fgls ./ sefgls);
        Rfgls       = 1 - ((ufgls'ufgls) / sumc((yfgls - meanc(yfgls))^2));
        R2_fgls     = 1 - (((rows(ufgls)-1) / (rows(ufgls)-rows(fgls)))
                    * (1 - Rfgls));
    endif;

else;

    fgls    = zeros(cols(wnxall),1);
    sefgls  = fgls;
    tfgls   = fgls;
    R2_fgls = -999;

endif;

clear idvec1, sel, uhat, sigv2, sigmu2, theta, yfgls,
        ylagfgls, xfgls, yfgls, xallfgls, ufgls, sig2_gls, Rfgls;
clear uols, Rols;


                /* FIRST DIFFERENCE MODEL ESTIMATORS */

/* USUAL OLS */

if sumc(tvx') ge 1;

    if sumc(timev[2:cols(timev)]' .* xvars[2:cols(xvars)]') ge 1;
        dxall   = delif(dylag~dx,e1);
    else;
        dxall   = delif(dylag,e1);
    endif;

else;

    dxall   = delif(dylag,e1);

endif;

if ismiss(dxall) eq 1;
	print;
	"WARNING - MISSING VALUES IN DELTA X'S (DELTA OLS)";
endif;

NTdols = rows(dxall);

dols    = inv(dxall'*dxall)*dxall'*delif(dy,e1);
if ismiss(inv(dxall'*dxall)) eq 1;

    dols    = zeros(cols(dxall),1);
    sedols  = dols;
    tdols   = dols;
    R2_dols = -999;

else;

    udols       = delif(dy,e1) - (dxall*dols);
    sig2_dls    = (udols'*udols)/(rows(udols)-rows(dols));
    sedols      = diag(sig2_dls*inv(dxall'*dxall))^0.5;
    tdols       = abs(dols ./ sedols);
    Rdols       = 1 - ((udols'udols) /
                    sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
    R2_dols     = 1 - (((rows(udols)-1) / (rows(udols)-rows(dols)))
                    * (1 - Rdols));
endif;

clear udols, sig2_dls, Rdols;


/* Which estimators? */

let mtype = linear wb gmm;

mtype = indcv(ivtype,mtype);


if mtype eq 2 or mtype eq 3;
        clear aha, seaha, taha, R2_aha, ahasarg, ahadfsrg, ahapsarg, ahaid, AR2aha, NTah,
              ahb,  seahb, tahb, R2_ahb, ahbsarg, ahbdfsrg, ahbpsarg, ahbid, AR2ahb, 
              ahc,  seahc, tahc, R2_ahc, ahcsarg, ahcdfsrg, ahcpsarg, ahcid, AR2ahc,
              ara,  seara, tara, R2_ara, arasarg, aradfsrg, arapsarg, araid, AR2ara, 
              arb,  searb, tarb, R2_arb, arbsarg, arbdfsrg, arbpsarg, arbid, AR2arb,
              arc,  searc, tarc, R2_arc, arcsarg, arcdfsrg, arcpsarg, arcid, AR2arc,
              aba,  seaba, taba, R2_aba, abasarg, abadfsrg, abapsarg, abaid, AR2aba, 
              abc,  seabc, tabc, R2_abc, abcsarg, abcdfsrg, abcpsarg, abcid, AR2abc,
              bb,   sebb,  tbb,  R2_bb,  bbsarg,  bbdfsrg,  bbpsarg,  bbid,  AR2bb,
              bb,   sebb,  tbb, R2_bb, bbsarg, bbdfsrg, bbpsarg, bbid, AR2bb,
              gmmi, segmmi, tgmmi, R2_gmmi, overidi,
              gmm,  segmm, tgmm, R2_gmm, overid;
endif;

if mtype eq 2;
    goto wbetc;
elseif mtype eq 3;
    clear wba,sewba,twba,R2_wba,wbb,sewbb,twbb,R2_wbb, 
        wbasarg, wbbsarg, wbadfsrg, wbbdfsrg, wbapsarg, wbbpsarg, wbaid, wbbid;
    goto gmmproc;
endif;


/* SIMPLE ANDERSON-HSIAO SSIV: AHA */

exogtv = selif(exog,timev');    @ re-define exog to be consistent with timev @
xvarstv= selif(xvars',timev');  @ re-define xvars to be consistent with timev @


if sumc(tvx') ge 1;

    xahtmp  = dx;

    if cols(timev) ge 2;
        xah     = delif(dylag~xahtmp,e2);
    else;
        xah     = delif(dylag,e2);
    endif;

    if sumc(tvz') ge 1;
        zahatmp = (selif(dx',tvexogx1')')~dz;
    else;
        zahatmp = xahtmp;
    endif;


    if cols(timev) ge 2;
        zaha    = delif(dylag2~zahatmp,e2);
    else;
        zaha    = delif(dylag2,e2);
    endif;

else;

    xah     = delif(dylag,e2);

    if sumc(tvz') ge 1;
        zaha    = delif(dylag2~dz,e2);
    else;
        zaha    = delif(dylag2,e2);
    endif;

endif;

NTah = rows(xah);


/* DEFINE DELTA (NEEDED FOR COVARIANCE MATRIX) */

phi1 = eye(tmax-3) * 2;
phirep = 2;
for phirep (2,rows(phi1),1);
    phi1[phirep,phirep-1] = -1;
    phi1[phirep-1,phirep] = -1;
endfor;
phi = eye(n) .*. phi1;


if ismiss(zaha) eq 1;
	print;
	"WARNING - MISSING VALUES IN A-H Z'S";
endif;

if cols(zaha) eq cols(xah);

    aha     = inv(zaha'xah)*zaha'*delif(dy,e2);

    if ismiss(inv(zaha'*xah)) eq 1;

        aha     = zeros(cols(xah),1);
        seaha   = aha;
        taha    = aha;
        R2_aha  = -999;

        ahasarg = -999;
        ahadfsrg= -999;
        ahapsarg= -999;

        ahaid   = ahasarg~ahadfsrg~ahapsarg;
        AR2aha  = -999;

    else;
    
        uaha        = delif(dy,e2) - (xah*aha);
        sig2_aha    = ((uaha'*uaha)/(rows(uaha)-rows(aha)))/2;
        varaha      = (inv(zaha'xah)*zaha'phi*zaha*inv(xah'zaha)) .* sig2_aha;
        seaha       = (diag(varaha))^0.5;
        taha        = abs(aha ./ seaha);
        Raha        = 1 - ((uaha'uaha) /
                        sumc((delif(dy,e2) - meanc(delif(dy,e2)))^2));
        R2_aha      = 1 - (((rows(uaha)-1) / (rows(uaha)-rows(aha)))
                        * (1 - Raha));
                                
        if (maxc(tvec)-minc(tvec)+1) ge 6;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-4,1));
       
            u2          = delif(lagn(uaha,2),sel1);
            ustar       = delif(uaha,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2aha = -999;
            
        endif;    

        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e2) .eq irep;
            utmp    = selif(uaha,sel);	            			
            ztmp    = selif(zaha,sel);               
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 6;
                                
                selu2   = delif(delif(idvec,e2),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
            
         endfor;

         ahasarg     = uaha'zaha*inv(v)*zaha'uaha;
         ahadfsrg    = cols(zaha)-cols(xah);
         ahapsarg    = cdfchic( ahasarg,ahadfsrg );
            
         ahaid       = ahasarg~ahadfsrg~ahapsarg;   

         if (maxc(tvec)-minc(tvec)+1) ge 6;
         
            ubot = part1 - (2*u2'delif(xah,sel1)*inv(zaha'xah)*
                            xah'zaha*
                            part2) +
                            (u2'delif(xah,sel1)*varaha*delif(xah,sel1)'u2);
         
            AR2aha = utop/(ubot)^0.5;
            
         else;               
            
            AR2aha = -999;
            
         endif;
            
    endif;

elseif cols(zaha) gt cols(xah);

    pzaha   = zaha*inv(zaha'*phi*zaha)*zaha';
    aha     = inv(xah'*pzaha*xah)*xah'*pzaha*delif(dy,e2);

    if ismiss(inv(xah'*pzaha*xah)) eq 1;

        aha     = zeros(cols(xah),1);
        seaha   = aha;
        taha    = aha;
        R2_aha  = -999;

        ahaid   = ones(1,3) .* -999;

        AR2aha  = -999;

        ahc     = zeros(cols(xah),1);
        seahc   = ahc;
        tahc    = ahc;
        R2_ahc  = -999;

        ahcid   = ones(1,3) .* -999;

        AR2ahc  = -999;

    else;

        uaha        = delif(dy,e2) - (xah*aha);
        sig2_aha    = ((uaha'*uaha)/(rows(uaha)-rows(aha)))/2;
        varaha      = inv(xah'pzaha*xah) .* sig2_aha;
        seaha       = (diag(varaha))^0.5;
        taha        = abs(aha ./ seaha);
        Raha        = 1 - ((uaha'uaha) /
                        sumc((delif(dy,e2) - meanc(delif(dy,e2)))^2));
        R2_aha      = 1 - (((rows(uaha)-1) / (rows(uaha)-rows(aha)))
                        * (1 - Raha));

       if (maxc(tvec)-minc(tvec)+1) ge 6;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-4,1));
       
            u2          = delif(lagn(uaha,2),sel1);
            ustar       = delif(uaha,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2aha = -999;
            
        endif;    


        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e2) .eq irep;
            utmp    = selif(uaha,sel);				
            ztmp    = selif(zaha,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);
  
              if (maxc(tvec)-minc(tvec)+1) ge 6;
                                
                selu2   = delif(delif(idvec,e2),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;

            
         endfor;

         ahasarg     = uaha'zaha*inv(v)*zaha'uaha;
         ahadfsrg    = cols(zaha)-cols(xah);
         ahapsarg    = cdfchic( ahasarg,ahadfsrg );

         ahaid       = ahasarg~ahadfsrg~ahapsarg;

         if (maxc(tvec)-minc(tvec)+1) ge 6;
         
            ubot = part1 - (2*u2'delif(xah,sel1)*
                            inv(xah'*pzaha*xah)*
                            xah'zaha*inv(zaha'*phi*zaha)
                            *part2) +
                            (u2'delif(xah,sel1)*varaha*delif(xah,sel1)'u2);
         
            AR2aha = utop/(ubot)^0.5;
            
         else;               
            
            AR2aha = -999;
            
         endif;


        /* ANDERSON-HSIAO PLUS OLSIV HAT: AHC */

        zahb    = zaha;
        uahb    = uaha;
        colsah  = cols(zahb);
        rowsah  = rows(zahb);
        zhatahc = zeros(cols(zahb),cols(zahb));

        irepst = 1;
        do while irepst le rowsah;
            irepfn = irepst + (rowsah/n) - 1;
            z1i = zahb[irepst:irepfn,.];
            uahbi = uahb[irepst:irepfn,1];
            if irepst eq 1;
                zhatahc = z1i'*uahbi*uahbi'*z1i;
            elseif irepst gt 1;
                zhatahc = zhatahc + (z1i'*uahbi*uahbi'*z1i);
            endif;
            irepst = irepfn + 1;
        endo;

        pzahc = zahb*inv(zhatahc)*zahb';
        ahc = inv(xah'*pzahc*xah)*xah'*pzahc*delif(dy,e2);

        if ismiss(inv(xah'*pzahc*xah)) eq 1;

            ahc     = zeros(cols(xah),1);
            seahc   = ahc;
            tahc    = ahc;
            R2_ahc  = -999;

            ahcid   = ones(1,3) .* -999;

            AR2ahc  = -999;

        else;

            uahc    = delif(dy,e2) - (xah*ahc);
            varahc  = inv(xah'pzahc*xah);
            seahc   = (diag(varahc))^0.5;
            tahc    = abs(ahc ./ seahc);
            Rahc    = 1 - ((uahc'uahc) /
                        sumc((delif(dy,e2) - meanc(delif(dy,e2)))^2));
            R2_ahc  = 1 - (((rows(uahc)-1) / (rows(uahc)-rows(ahc)))
                        * (1 - Rahc));

            if (maxc(tvec)-minc(tvec)+1) ge 6;

                sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-4,1));
       
                u2          = delif(lagn(uahc,2),sel1);
                ustar       = delif(uahc,sel1);
                
                utop        = u2'ustar;
                part1       = 0; 
                part2       = 0;
                        
            else;
        
                AR2ahc = -999;
            
            endif;    


            clear v;
            for irep (1, n, 1);

                sel     = delif(idvec,e2) .eq irep;
                utmp    = selif(uahc,sel);				
                ztmp    = selif(zahb,sel);   
                v       = v + (ztmp'utmp*utmp'ztmp);

                if (maxc(tvec)-minc(tvec)+1) ge 6;
                                
                    selu2   = delif(delif(idvec,e2),sel1) .eq irep;
                    part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                    part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

                endif;
           
            endfor;

            ahcsarg     = uahc'zahb*inv(v)*zahb'uahc;
            ahcdfsrg    = cols(zahb)-cols(xah);
            ahcpsarg    = cdfchic( ahcsarg,ahcdfsrg );
            
            ahcid       = ahcsarg~ahcdfsrg~ahcpsarg;

            if (maxc(tvec)-minc(tvec)+1) ge 6;
         
                ubot =  part1 - (2*u2'delif(xah,sel1)*
                        inv(xah'*pzahc*xah)*
				        xah'zahb*inv(zhatahc)
				        *part2) +
                        (u2'delif(xah,sel1)*varahc*delif(xah,sel1)'u2);
         
                AR2ahc = utop/(ubot)^0.5;
            
            else;               
            
                AR2ahc = -999;
            
            endif;

        endif;

    endif;

    ahb     = aha;
    seahb   = seaha;
    tahb    = taha;
    R2_ahb  = R2_aha;

    ahbid   = ahaid;

    AR2ahb  = AR2aha;

    goto arellano;

endif;


/* ANDERSON-HSIAO PLUS OLSIV: AHB (ie expanded instrument set) */

if sumc(tvx') ge 1;

    if sumc(tvz') ge 1;
        zahbtmp = (selif(dxlag',tvexogx1')')~(selif(dx',tvexogx1')')~dz;
    else;
        zahbtmp = xahtmp;
    endif;

    if cols(timev) ge 2;
        zahb    = delif(dylag2~zahbtmp,e2);
    else;
        zahb    = delif(dylag2,e2);
    endif;

	if ismiss(zahb) eq 1;
		print;
		"WARNING - MISSING VALUES IN A-H B & C Z'S";			
	endif;


    pzahb   = zahb*inv(zahb'*phi*zahb)*zahb';
    ahb     = inv(xah'*pzahb*xah)*xah'*pzahb*delif(dy,e2);

    if ismiss(inv(xah'*pzahb*xah)) eq 1;

        ahb     = zeros(cols(xah),1);
        seahb   = ahb;
        tahb    = ahb;
        R2_ahb  = -999;

        ahbid   = ones(3,1) .* -999;
        AR2ahb  = -999;

        ahc     = zeros(cols(xah),1);
        seahc   = ahc;
        tahc    = ahc;
        R2_ahc  = -999;

        ahcid   = ones(3,1) .* -999;

        AR2ahc  = -999;

    else;

        uahb        = delif(dy,e2) - (xah*ahb);
        sig2_ahb    = ((uahb'*uahb)/(rows(uahb)-rows(ahb)))/2;
        varahb      = inv(xah'pzahb*xah) .* sig2_ahb;
        seahb       = (diag(varahb))^0.5;
        tahb        = abs(ahb ./ seahb);
        Rahb        = 1 - ((uahb'uahb) /
                        sumc((delif(dy,e2) - meanc(delif(dy,e2)))^2));
        R2_ahb      = 1 - (((rows(uahb)-1) / (rows(uahb)-rows(ahb)))
                        * (1 - Rahb));

       if (maxc(tvec)-minc(tvec)+1) ge 6;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-4,1));
       
            u2          = delif(lagn(uahb,2),sel1);
            ustar       = delif(uahb,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2ahb = -999;
            
        endif;    

        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e2) .eq irep;
            utmp    = selif(uahb,sel);				
            ztmp    = selif(zahb,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 6;
                                
                selu2   = delif(delif(idvec,e2),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;

            
        endfor;

        ahbsarg     = uahb'zahb*inv(v)*zahb'uahb;
        ahbdfsrg    = cols(zahb)-cols(xah);
        ahbpsarg    = cdfchic( ahbsarg,ahbdfsrg );

        ahbid       = ahbsarg~ahbdfsrg~ahbpsarg;

        if (maxc(tvec)-minc(tvec)+1) ge 6;
         
            ubot = part1 - (2*u2'delif(xah,sel1)*
                            inv(xah'*pzahb*xah)*
                            xah'zahb*inv(zahb'*phi*zahb)*
                            part2) +
                            (u2'delif(xah,sel1)*varahb*delif(xah,sel1)'u2);
         
            AR2ahb = utop/(ubot)^0.5;
            
         else;               
            
            AR2ahb = -999;
            
         endif;


        /* ANDERSON-HSIAO PLUS OLSIV HAT: AHC */

        colsah = cols(zahb);
        rowsah = rows(zahb);
        zhatahc = zeros(cols(zahb),cols(zahb));

        irepst = 1;
        do while irepst le rowsah;
            irepfn = irepst + (rowsah/n) - 1;
            z1i = zahb[irepst:irepfn,.];
            uahbi = uahb[irepst:irepfn,1];

            if irepst eq 1;
                zhatahc = z1i'*uahbi*uahbi'*z1i;
            elseif irepst gt 1;
                zhatahc = zhatahc + (z1i'*uahbi*uahbi'*z1i);
            endif;

            irepst = irepfn + 1;
        endo;

        pzahc = zahb*inv(zhatahc)*zahb';
        ahc = inv(xah'*pzahc*xah)*xah'*pzahc*delif(dy,e2);
        if ismiss(inv(xah'*pzahc*xah)) eq 1;

            ahc     = zeros(cols(xah),1);
            seahc   = ahc;
            tahc    = ahc;
            R2_ahc  = -999;

            ahcid   = ones(3,1) .* -999;

            AR2ahc  = -999;

        else;

            uahc        = delif(dy,e2) - (xah*ahc);
            varahc      = inv(xah'pzahb*xah);
            seahc       = (diag(varahc))^0.5;
            tahc        = abs(ahc ./ seahc);
            Rahc        = 1 - ((uahc'uahc) /
                            sumc((delif(dy,e2) - meanc(delif(dy,e2)))^2));
            R2_ahc      = 1 - (((rows(uahc)-1) / (rows(uahc)-rows(ahc)))
                            * (1 - Rahc));

            if (maxc(tvec)-minc(tvec)+1) ge 6;

                sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-4,1));
       
                u2          = delif(lagn(uahc,2),sel1);
                ustar       = delif(uahc,sel1);
                
                utop        = u2'ustar;
                part1       = 0; 
                part2       = 0;
                        
            else;
        
                AR2ahc = -999;
            
            endif;    


            clear v;
            for irep (1, n, 1);

                sel     = delif(idvec,e2) .eq irep;
                utmp    = selif(uahc,sel);				
                ztmp    = selif(zahb,sel);   
                v       = v + (ztmp'utmp*utmp'ztmp);

                if (maxc(tvec)-minc(tvec)+1) ge 6;
                                
                    selu2   = delif(delif(idvec,e2),sel1) .eq irep;
                    part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                    part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

                endif;
            
            endfor;

            ahcsarg     = uahc'zahb*inv(v)*zahb'uahc;
            ahcdfsrg    = cols(zahb)-cols(xah);
            ahcpsarg    = cdfchic( ahcsarg,ahcdfsrg );

            ahcid       = ahcsarg~ahcdfsrg~ahcpsarg;

            if (maxc(tvec)-minc(tvec)+1) ge 6;
         
                ubot = part1 - (2*u2'delif(xah,sel1)*
					   inv(xah'*pzahc*xah)*
					   xah'zahb*inv(zhatahc)*
					   part2) +
                                (u2'delif(xah,sel1)*varahc*delif(xah,sel1)'u2);
         
                AR2ahc = utop/(ubot)^0.5;
            
            else;               
            
                AR2ahc = -999;
            
            endif;

        endif;

    endif;

elseif sumc(tvx') eq 0;

    ahb     = zeros(cols(xah),1);
    seahb   = ahb;
    tahb    = ahb;
    R2_ahb  = -999;

    ahbid   = ones(3,1) .* -999;

    AR2ahb  = -999;

    ahc     = zeros(cols(xah),1);
    seahc   = ahc;
    tahc    = ahc;
    R2_ahc  = -999;

    ahcid   = ones(3,1) .* -999;

    AR2ahc  = -999;

endif;


arellano:


clear xah, zaha, uaha, sig2_aha, R2_aha;
clear phi1, phi, phirep, zaha, zahb, pzahb;
clear uahb, sig2_ahb, Rahb, colsah, rowsah, zhatahc, z1i, uahbi;
clear pzahc, uahc, sig2_ahc, Rahc;
clear ahasarg, ahadfsrg, ahapsarg, ahbsarg, ahbdfsrg, ahbpsarg, ahcsarg, ahcdfsrg, ahcpsarg;
clear varaha, varahb, varahc;

/* SIMPLE ARELLANO SSIV: ARA */

if sumc(tvx') ge 1;

    xartmp  = dx;

    if cols(timev) ge 2;
        xar     = delif(dylag~xartmp,e1);
    else;
        xar     = delif(dylag,e1);
    endif;

    if sumc(tvz') ge 1;
         zaratmp = (selif(dx',tvexogx1')')~dz;
    else;
         zaratmp = xartmp;
    endif;

    if cols(timev) ge 2;
        zara = delif(ylag2~zaratmp,e1);
    else;
        zara = delif(ylag2,e1);
    endif;

else;

    xar  = delif(dylag,e1);

    if sumc(tvz') ge 1;
        zara = delif(ylag2~dz,e1);
    else;
        zara = delif(dylag2,e1);
    endif;

endif;



/* CONSTRUCT DELTA COVARIANCE MATRIX */

phi1 = eye(tmax-2) * 2;
for phirep (2,rows(phi1),1);
    phi1[phirep,phirep-1] = -1;
    phi1[phirep-1,phirep] = -1;
endfor;
phi = eye(n) .*. phi1;

if ismiss(zara) eq 1;
	print;
	"WARNING - MISSING VALUES IN AR A'S Z'S";			
endif;


if cols(zara) eq cols(xar);

    ara = inv(zara'*xar)*zara'*delif(dy,e1);

    if ismiss(inv(zara'*xar)) eq 1;

        ara     = zeros(cols(xar),1);
        seara   = ara;
        tara    = ara;
        R2_ara  = -999;

        araid   = ones(3,1) .* -999;

        AR2ara  = -999;

    else;

        uara        = delif(dy,e1) - (xar*ara);
        sig2_ara    = ((uara'*uara)/(rows(uara)-rows(ara)))/2;

        varara      = (inv(zara'xar)*zara'phi*zara*inv(xar'zara)) .* sig2_ara;
        seara       = (diag(varara))^0.5;
        tara        = abs(ara ./ seara);
        Rara        = 1 - ((uara'uara) /
                        sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
        R2_ara      = 1 - (((rows(uara)-1) / (rows(uara)-rows(ara)))
                        * (1 - Rara));

        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(uara,2),sel1);
            ustar       = delif(uara,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2ara = -999;
            
        endif;    


        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e1) .eq irep;
            utmp    = selif(uara,sel);			
            ztmp    = selif(zara,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
       
            
        endfor;

        arasarg     = uara'zara*inv(v)*zara'uara;
        aradfsrg    = cols(zara)-cols(xar);
        arapsarg    = cdfchic( arasarg,aradfsrg );

        araid       = arasarg~aradfsrg~arapsarg;


        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xar,sel1)*
					inv(zara'xar)*
					xar'zara*
					part2) +
                            (u2'delif(xar,sel1)*varara*delif(xar,sel1)'u2);
         
            AR2ara = utop/(ubot)^0.5;
            
         else;               
            
            AR2ara = -999;
            
         endif;

    endif;

elseif cols(zara) gt cols(xar);

    pzara   = zara*inv(zara'*phi*zara)*zara';
    ara     = inv(xar'*pzara*xar)*xar'*pzara*delif(dy,e1);

    if ismiss(inv(xar'pzara*xar)) eq 1;

        ara     = zeros(cols(xar),1);
        seara   = ara;
        tara    = ara;
        R2_ara  = -999;

        araid   = ones(3,1) .* -999;

        AR2ara  = -999;

        arc     = zeros(cols(xar),1);
        searc   = arc;
        tarc    = arc;
        R2_arc  = -999;

        arcid   = ones(3,1) .* -999;

        AR2ara  = -999;

    else;

        uara        = delif(dy,e1) - (xar*ara);
        sig2_ara    = ((uara'*uara)/(rows(uara)-rows(ara)))/2;
        varara      = (inv(xar'pzara*xar)) .*sig2_ara;
        seara       = (diag(varara))^0.5;
        tara        = abs(ara ./ seara);
        Rara        = 1 - ((uara'uara) /
                        sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
        R2_ara      = 1 - (((rows(uara)-1) / (rows(uara)-rows(ara)))
                        * (1 - Rara));

        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(uara,2),sel1);
            ustar       = delif(uara,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2ara = -999;
            
        endif;    


        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e1) .eq irep;
            utmp    = selif(uara,sel);			
            ztmp    = selif(zara,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
            
        endfor;

        arasarg     = uara'zara*inv(v)*zara'uara;
        aradfsrg    = cols(zara)-cols(xar);
        arapsarg    = cdfchic( arasarg,aradfsrg );

        araid       = arasarg~aradfsrg~arapsarg;

        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xar,sel1)*
					inv(xar'*pzara*xar)*
					xar'zara*inv(zara'*phi*zara)*
					part2) +
                            (u2'delif(xar,sel1)*varara*delif(xar,sel1)'u2);
         
            AR2ara = utop/(ubot)^0.5;
            
         else;               
            
            AR2ara = -999;
            
         endif;


        /* ARELLANO PLUS OLSIV HAT: ARC */

        zarb    = zara;
        uarb    = uara;
        colsar  = cols(zarb);
        rowsar  = rows(zarb);
        zhatarc = zeros(cols(zarb),cols(zarb));

        irepst = 1;
        do while irepst le rowsar;
            irepfn  = irepst + (rowsar/n) - 1;
            uarbi   = uarb[irepst:irepfn];
            z1i     = zarb[irepst:irepfn,.];
            if irepst eq 1;
                zhatarc = z1i'*uarbi*uarbi'*z1i;
            elseif irepst gt 1;
                zhatarc = zhatarc + (z1i'*uarbi*uarbi'*z1i);
            endif;
        irepst = irepfn + 1;
        endo;

        pzarc = zarb*inv(zhatarc)*zarb';
        arc = inv(xar'*pzarc*xar)*xar'*pzarc*delif(dy,e1);

        if ismiss(inv(xar'*pzarc*xar)) eq 1;

            arc     = zeros(cols(xar),1);
            searc   = arc;
            tarc    = arc;
            R2_arc  = -999;

            arcid   = ones(3,1) .* -999;

            AR2arc  = -999;

        else;

            uarc        = delif(dy,e1) - (xar*arc);
            vararc      = inv(xar'pzarc*xar);
            searc       = (diag(vararc))^0.5;
            tarc        = abs(arc ./ searc);
            Rarc        = 1 - ((uarc'uarc) /
                            sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
            R2_arc      = 1 - (((rows(uarc)-1) / (rows(uarc)-rows(arc)))
                            * (1 - Rarc));

            if (maxc(tvec)-minc(tvec)+1) ge 5;

                sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
                u2          = delif(lagn(uarc,2),sel1);
                ustar       = delif(uarc,sel1);
                
                utop        = u2'ustar;
                part1       = 0; 
                part2       = 0;
                        
            else;
        
                AR2arc      = -999;
            
            endif;    


            clear v;
            for irep (1, n, 1);

                sel     = delif(idvec,e1) .eq irep;
                utmp    = selif(uarc,sel);			
                ztmp    = selif(zarb,sel);   
                v       = v + (ztmp'utmp*utmp'ztmp);

                if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                    selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                    part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                    part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

                endif;

            
            endfor;

            arcsarg     = uarc'zarb*inv(v)*zarb'uarc;
            arcdfsrg    = cols(zarb)-cols(xar);
            arcpsarg    = cdfchic( arcsarg,arcdfsrg );

            arcid       = arcsarg~arcdfsrg~arcpsarg;

            if (maxc(tvec)-minc(tvec)+1) ge 5;
         
                ubot = part1 - (2*u2'delif(xar,sel1)*
					   inv(xar'*pzarc*xar)*
					   xar'zarb*inv(zhatarc)*
					   part2) +
                                (u2'delif(xar,sel1)*vararc*delif(xar,sel1)'u2);
         
                AR2arc = utop/(ubot)^0.5;
            
            else;               
            
                AR2arc = -999;
            
            endif;


        endif;


    endif;

    arb         = ara;
    searb       = seara;
    tarb        = tara;
    R2_arb      = R2_ara;

    arbid       = araid;

    AR2arb      = AR2ara;
    
    goto a_bond;

endif;


/* ARELLANO PLUS OLSIV: ARB */


if sumc(tvx') ge 1;

    if sumc(tvz') ge 1;
        zaratmp = (selif(dxlag',tvexogx1')')~(selif(dx',tvexogx1')')~dz;
    else;
        zaratmp = xartmp;
    endif;

    if cols(timev) ge 2;
        zarb    = delif(ylag2~zaratmp,e1);
    else;
        zarb    = delif(ylag2,e1);
    endif;

    pzarb   = zarb*inv(zarb'*phi*zarb)*zarb';

    arb     = inv(xar'*pzarb*xar)*xar'*pzarb*delif(dy,e1);

	if ismiss(zarb) eq 1;
		print;
		"WARNING - MISSING VALUES IN AR'S B & C Z'S";			
	endif;


    if ismiss(inv(xar'pzarb*xar)) eq 1;

        arb     = zeros(cols(xar),1);
        searb   = arb;
        tarb    = arb;
        R2_arb  = -999;

        arbid   = ones(3,1) .* -999;

        AR2arb  = -999;

        arc     = zeros(cols(xar),1);
        searc   = arc;
        tarc    = arc;
        R2_arc  = -999;

        arcid   = ones(3,1) .* -999;

        AR2arc  = -999;
        
    else;

        uarb        = delif(dy,e1) - (xar*arb);
        sig2_arb    = ((uarb'*uarb)/(rows(uarb)-rows(arb)))/2;
        vararb      = (inv(xar'pzarb*xar)) .* sig2_arb;
        searb       = (diag(vararb))^0.5;
        tarb        = abs(arb ./ searb);
        Rarb        = 1 - ((uarb'uarb) /
                        sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
        R2_arb      = 1 - (((rows(uarb)-1) / (rows(uarb)-rows(arb)))
                        * (1 - Rarb));
                        
        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(uarb,2),sel1);
            ustar       = delif(uarb,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2arb = -999;
            
        endif;    


        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e1) .eq irep;
            utmp    = selif(uarb,sel);			
            ztmp    = selif(zarb,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
            
        endfor;

        arbsarg     = uarb'zarb*inv(v)*zarb'uarb;
        arbdfsrg    = cols(zarb)-cols(xar);
        arbpsarg    = cdfchic( arbsarg,arbdfsrg );

        arbid       = arbsarg~arbdfsrg~arbpsarg;

        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xar,sel1)*
					inv(xar'*pzarb*xar)*
					xar'zarb*inv(zarb'*phi*zarb)*
					part2) +
                            (u2'delif(xar,sel1)*vararb*delif(xar,sel1)'u2);
         
            AR2arb = utop/(ubot)^0.5;
            
         else;               
            
            AR2arb = -999;
            
         endif;



        /* ARELLANO PLUS OLSIV HAT: ARC */

        colsar  = cols(zarb);
        rowsar  = rows(zarb);
        zhatarc = zeros(cols(zarb),cols(zarb));

        irepst = 1;
        do while irepst le rowsar;
            irepfn  = irepst + (rowsar/n) - 1;
            uarbi   = uarb[irepst:irepfn];
            z1i     = zarb[irepst:irepfn,.];
            if irepst eq 1;
                zhatarc = z1i'*uarbi*uarbi'*z1i;
            elseif irepst gt 1;
                zhatarc = zhatarc + (z1i'*uarbi*uarbi'*z1i);
            endif;
        irepst = irepfn + 1;
        endo;

        pzarc = zarb*inv(zhatarc)*zarb';
        arc = inv(xar'*pzarc*xar)*xar'*pzarc*delif(dy,e1);

        if ismiss(inv(xar'*pzarc*xar)) eq 1;

            arc     = zeros(cols(xar),1);
            searc   = arc;
            tarc    = arc;
            R2_arc  = -999;

            arcid   = ones(3,1) .* -999;

            AR2arc  = -999;

        else;
        
            uarc        = delif(dy,e1) - (xar*arc);
            vararc      = inv(xar'pzarc*xar);
            searc       = (diag(vararc))^0.5;
            tarc        = abs(arc ./ searc);
            Rarc        = 1 - ((uarc'uarc) /
                            sumc((delif(dy,e1) - meanc(delif(dy,e1)))^2));
            R2_arc     = 1 - (((rows(uarc)-1) / (rows(uarc)-rows(arc)))
                        * (1 - Rarc));

            if (maxc(tvec)-minc(tvec)+1) ge 5;

                sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
                u2          = delif(lagn(uarc,2),sel1);
                ustar       = delif(uarc,sel1);
                
                utop        = u2'ustar;
                part1       = 0; 
                part2       = 0;
                        
            else;
        
                AR2arc = -999;
            
            endif;    


            clear v;
            for irep (1, n, 1);

                sel     = delif(idvec,e1) .eq irep;
                utmp    = selif(uarc,sel);			
                ztmp    = selif(zarb,sel);   
                v       = v + (ztmp'utmp*utmp'ztmp);

                if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                    selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                    part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                    part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

                endif;
            
            endfor;

            arcsarg     = uarc'zarb*inv(v)*zarb'uarc;
            arcdfsrg    = cols(zarb)-cols(xar);
            arcpsarg    = cdfchic( arcsarg,arcdfsrg );

            arcid       = arcsarg~arcdfsrg~arcpsarg;

            if (maxc(tvec)-minc(tvec)+1) ge 5;
         
                ubot = part1 - (2*u2'delif(xar,sel1)*
					   inv(xar'*pzarc*xar)*
					   xar'zarb*inv(zhatarc)*
					   part2) +
                                (u2'delif(xar,sel1)*vararc*delif(xar,sel1)'u2);
         
                AR2arc = utop/(ubot)^0.5;
            
            else;               
            
                AR2arc = -999;
            
            endif;


        endif;

    endif;

elseif sumc(tvx') eq 0;

    arb     = zeros(cols(xar),1);
    searb   = arb;
    tarb    = arb;
    R2_arb  = -999;

    arbid   = ones(3,1) .* -999;
    
    AR2arb  = -999;
    
    arc     = zeros(cols(xar),1);
    searc   = arc;
    tarc    = arc;
    R2_arc  = -999;

    arcid   = ones(3,1) .* -999;

    AR2arc  = -999;
    
endif;


clear xar, zara, uara, sig2_ara, R2_ara, phii, phii1, zarb, pzarb, uarb, sig2_arb, R2_arb,
zhatarc, pzarc, uarc, sig2_arc, R2_arc, colsar, rowsar, uarbi, arasarg, aradfsrg, arapsarg, 
arbsarg, arbdfsrg, arbpsarg, arcsarg, arcdfsrg, arcpsarg, varara, vararb, vararc, e, xols,
xfull, xgmm, yols;


a_bond:


                /* ARELLANO-BOND OLSIV: ABA */


if sumc(tvx') ge 1;

    xarbtmp = dx;

    if cols(timev) ge 2;
        xarb    = delif(dylag~xarbtmp,e1);
    else;
        xarb    = delif(dylag,e1);
    endif;

    if sumc(tvz') ge 1;
        zarbtmp = (selif(dx',tvexogx1')')~dz;
    else;
        zarbtmp = xarbtmp;
    endif;

    if cols(timev) ge 2;
        zarbtmp = delif(zarbtmp,e1);
        abz1    = abz~zarbtmp;
    else;
        abz1    = abz;
    endif;

else;

    xarb    = delif(dylag,e1);

    if sumc(tvz') ge 1;
        abz1    = abz~delif(dz,e1);
    else;
        abz1    = abz;
    endif;

endif;

dyab = delif(dy,e1);


if ismiss(abz1) eq 1;
	print;
	"WARNING - MISSING VALUES IN AB Z'S";
endif;

pzaba   = abz1*inv(abz1'phi*abz1)*abz1';
aba     = inv(xarb'pzaba*xarb)*xarb'pzaba*dyab;

if ismiss(inv(xarb'*pzaba*xarb)) eq 1;

    aba     = zeros(cols(xarb),1);
    seaba   = aba;
    taba    = aba;
    R2_aba  = -999;

    abaid   = ones(3,1) .* -999;

    AR2aba  = -999;

    abc     = zeros(cols(xarb),1);
    seabc   = abc;
    tabc    = abc;
    R2_abc  = -999;

    abcid   = ones(3,1) .* -999;

    AR2abc  = -999;

else;

    uaba        = dyab - (xarb*aba);
    sig2_aba    = ((uaba'*uaba)/(rows(uaba)-rows(aba)))/2;

    varaba      = (inv(xarb'abz1*inv(abz1'*(phi)*abz1)*abz1'xarb)) .* sig2_aba;

    seaba       = (diag(varaba))^0.5;
    taba        = abs(aba ./ seaba);
    Raba        = 1 - ((uaba'uaba) /
                    sumc((dyab - meanc(dyab))^2));
    R2_aba      = 1 - (((rows(uaba)-1) / (rows(uaba)-rows(aba)))
                    * (1 - Raba));

    if (maxc(tvec)-minc(tvec)+1) ge 5;

        sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
        u2          = delif(lagn(uaba,2),sel1);
        ustar       = delif(uaba,sel1);
                
        utop        = u2'ustar;
        part1       = 0; 
        part2       = 0;
                        
    else;
        
        AR2aba = -999;
            
    endif;    


    clear v;
    for irep (1, n, 1);

        sel     = delif(idvec,e1) .eq irep;
        utmp    = selif(uaba,sel);			
        ztmp    = selif(abz1,sel);   
        v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
            
    endfor;

        abasarg     = uaba'abz1*inv(v)*abz1'uaba;
        abadfsrg    = cols(abz1)-cols(xarb);
        abapsarg    = cdfchic( abasarg,abadfsrg );

        abaid       = abasarg~abadfsrg~abapsarg;

        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xarb,sel1)*
					inv(xarb'pzaba*xarb)*
					xarb'abz1*inv(abz1'phi*abz1)*
					part2) +
                            (u2'delif(xarb,sel1)*varaba*delif(xarb,sel1)'u2);
         
            AR2aba = utop/(ubot)^0.5;
            
         else;               
            
            AR2aba = -999;
            
         endif;


           /* ARELLANO-BOND GMM: ABC */

    colsab1 = cols(abz1);
    rowsab1 = rows(abz1);

    irepst = 1;
    do while irepst le rowsab1;
        irepfn  = irepst + (rowsab1/n) - 1;
        uabai   = uaba[irepst:irepfn];
        z1i     = abz1[irepst:irepfn,.];
        if irepst eq 1;
            zhatabc = z1i'uabai*uabai'z1i;
        elseif irepst gt 1;
            zhatabc = zhatabc + (z1i'uabai*uabai'z1i);
        endif;
    irepst = irepfn + 1;
    endo;

    pzabc = abz1*inv(zhatabc)*abz1';

    abc = inv(xarb'*pzabc*xarb)*xarb'pzabc*dyab;

    if ismiss(inv(xarb'*pzabc*xarb)) eq 1;

        abc     = zeros(cols(xarb),1);
        seabc   = abc;
        tabc    = abc;
        R2_abc  = -999;

        abcid   = ones(1,3) .* -999;

    else;

        uabc        = dyab - (xarb*abc);
        sig2_abc    = ((uabc'*uabc)/(rows(uabc)-rows(abc)))/2;
        varabc      = inv(xarb'pzabc*xarb);
        seabc       = (diag(varabc))^0.5;
        tabc        = abs(abc ./ seabc);
        Rabc        = 1 - ((uabc'uabc) /
                        sumc((dyab - meanc(dyab))^2));
        R2_abc      = 1 - (((rows(uabc)-1) / (rows(uabc)-rows(abc)))
                        * (1 - Rabc));

        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec),1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(uabc,2),sel1);
            ustar       = delif(uabc,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2abc = -999;
            
        endif;    


        clear v;
        for irep (1, n, 1);

            sel     = delif(idvec,e1) .eq irep;
            utmp    = selif(uabc,sel);			
            ztmp    = selif(abz1,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);


            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(delif(idvec,e1),sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
          
        endfor;
        
        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xarb,sel1)*
					inv(xarb'*pzabc*xarb)*
					xarb'abz1*inv(zhatabc)*
					part2) +
                            (u2'delif(xarb,sel1)*varabc*delif(xarb,sel1)'u2);
         
            AR2abc = utop/(ubot)^0.5;
         
         else;               
            
            AR2abc = -999;
            
         endif;

        abcsarg     = uabc'abz1*inv(v)*abz1'uabc;
        abcdfsrg    = cols(abz1)-cols(xarb);
        abcpsarg    = cdfchic( abcsarg,abcdfsrg );

        abcid       = abcsarg~abcdfsrg~abcpsarg;


    endif;

endif;


                    /* BLUNDEL-BOND ESTIMATOR */

idbb    = delif(idvec,e1)|delif(idvec,e1);
ytild   = dyab|delif(y,e1);            
xbb     = zeros(rows(xarb),1)~xarb[.,1];            @ need to pad out delta x's for any time invariant x's @ 

bbrep = 2;

for xrep (2,cols(xvars),1);

    if xvars[xrep] eq 1;
        xbb     = xbb~xarb[.,bbrep];
        bbrep   = bbrep + 1; 
    else;
        xbb     = xbb~zeros(rows(xbb),1);
    endif;

endfor;

xbb = xbb[.,1:sumc(xvars')+1];

if sumc(xvars') gt 1;
    xtild   = xbb|(delif(xolstmp[.,1]~ylag~xolstmp[.,2:cols(xolstmp)],e1));
else;
    xtild   = xbb|(delif(ones(rows(x),1)~ylag,e1));
endif;

z_l = zeros(1,tmax-2);

for irep (1, n, 1);
    sel     = (idvec .eq irep) .* (1-e1);
    z_l     = z_l|(eye(tmax-2) .* selif(dylag,sel));
endfor;

z_l = z_l[2:rows(z_l),.];

if sumc(ivs) ge 1;
    z_l = z_l~delif(selif(x',ivs)',e1);
endif;

z_l = z_l~delif(selif(x',exog .* xvars')',e1);

bbz     = (abz1~zeros(rows(abz1),cols(z_l)))|(zeros(rows(z_l),cols(abz1))~z_l);

pzbb    = bbz*inv(bbz'bbz)*bbz';

bb      = inv(xtild'pzbb*xtild)*xtild'pzbb*ytild;


if ismiss(inv(xtild'pzbb*xtild)) eq 1;

    bb      = zeros(cols(xtild),1);
    sebb    = bb;
    tbb     = bb;
    R2_bb   = -999;

    bbid    = ones(3,1) .* -999;

    AR2bb   = -999;

else;

    /* TWO-STEP BLUNDELL-BOND GMM: BB */

    ubb        = ytild - (xtild*bb);

    zhatbb = 0;
    for irep (1, n, 1);

        sel     = idbb .eq irep;
        utmp    = selif(ubb,sel);			
        ztmp    = selif(bbz,sel);   
        zhatbb  = zhatbb + (ztmp'utmp*utmp'ztmp);
            
    endfor;

    zhatbb  = zhatbb ./ n; 

    pzhatbb = bbz*inv(zhatbb)*bbz';

    bb = inv(xtild'*pzhatbb*xtild)*xtild'pzhatbb*ytild;


    if ismiss(inv(xtild'*pzhatbb*xtild)) eq 1;

        bb      = zeros(cols(xtild),1);
        sebb    = bb;
        tbb     = bb;
        R2_bb   = -999;
        bbid    = ones(1,3) .* -999;

    else;

        ubb         = ytild - (xtild*bb);

        varbb       = inv(xtild'pzhatbb*xtild);
        sebb        = (diag(varbb))^0.5;

        if ismiss(sebb) eq 1; 
            varbb   = inv(xtild'pzbb*xtild);
            sebb    = (diag(varbb))^0.5;
        endif;

        tbb         = abs(bb ./ sebb);
        Rbb         = 1 - ((ubb'ubb) /
                        sumc((ytild - meanc(ytild))^2));
        R2_bb       = 1 - (((rows(ubb)-1) / (rows(ubb)-rows(bb)))
                        * (1 - Rbb));

        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec)*2,1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(ubb,2),sel1);
            ustar       = delif(ubb,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2bb       = -999;
            
        endif;    

        varbb   = inv(xtild'pzhatbb*xtild);
        sebb    = (diag(varbb))^0.5;

        if ismiss(sebb) eq 1; 
            varbb   = inv(xtild'pzbb*xtild);
            sebb    = (diag(varbb))^0.5;
        endif;


        tbb     = abs(bb ./ sebb);
        Rbb     = 1 - ((ubb'ubb) /
                    sumc((ytild - meanc(ytild))^2));
        R2_bb   = 1 - (((rows(ubb)-1) / (rows(ubb)-rows(bb)))
                    * (1 - Rbb));

        if (maxc(tvec)-minc(tvec)+1) ge 5;

            sel1        = ones(maxc(idvec)*2,1) .*. (ones(2,1)|zeros(maxc(tvec)-minc(tvec)-3,1));
            u2          = delif(lagn(ubb,2),sel1);
            ustar       = delif(ubb,sel1);
                
            utop        = u2'ustar;
            part1       = 0; 
            part2       = 0;
                        
        else;
        
            AR2bb       = -999;
            
        endif;    

        clear v;
        for irep (1, n, 1);

            sel     = idbb .eq irep;
            utmp    = selif(ubb,sel);			
            ztmp    = selif(bbz,sel);   
            v       = v + (ztmp'utmp*utmp'ztmp);

            if (maxc(tvec)-minc(tvec)+1) ge 5;
                                
                selu2   = delif(idbb,sel1) .eq irep;
                part1   = part1 + (selif(u2,selu2)'selif(ustar,selu2)*selif(ustar,selu2)'selif(u2,selu2));
                part2   = part2 + (ztmp'utmp*selif(ustar,selu2)'selif(u2,selu2));

            endif;
            
        endfor;

        bbsarg      = abs(ubb'bbz*inv(v)*bbz'ubb);

        bbdfsrg     = cols(bbz)-cols(xtild);
        bbpsarg     = cdfchic( bbsarg,bbdfsrg );

        bbid        = bbsarg~bbdfsrg~bbpsarg;

        if (maxc(tvec)-minc(tvec)+1) ge 5;
         
            ubot = part1 - (2*u2'delif(xtild,sel1)*
					inv(xtild'*pzhatbb*xtild)*
					xtild'bbz*inv(zhatbb)*
					part2) +
                            (u2'delif(xtild,sel1)*varbb*delif(xtild,sel1)'u2);
         
            AR2bb = utop/(ubot)^0.5;
            
         else;               
            
            AR2bb = -999;
            
         endif;

    endif;

endif;


clear dxall, abz1, pzaba, uaba, sig2_aba, Raba, colsab1, rowsab1, zhatabc;
clear irepst, irepfn, uabai, z1i, pzabc, uabc, Rabc;

if mtype eq 1;
    clear wba,  sewba, twba, R2_wba, wbb,  sewbb, twbb, R2_wbb;
    clear wbasarg, wbbsarg, wbadfsrg, wbbdfsrg, wbapsarg, wbbpsarg, wbaid, wbbid;
    clear gmmi, segmmi,tgmmi,R2_gmmi, overidi;
    clear gmm,  segmm, tgmm, R2_gmm, overid;
    goto finish;
endif;

wbetc:

        /* WANSBEEK-BEKKER ESTIMATOR: A'Y+ AS INSTRUMENT: WBA */

ai      = zeros(tmax,tmax-1);
arep    = 1;
do while arep lt tmax-1;
    ai[arep,arep]   = -1;
    ai[arep,arep+1] = 1;
    arep = arep + 1;
endo;

aivec = vec(ai);

ct      = zeros(tmax-1,tmax*(tmax-1));
eyerep  = tmax-1;

for crep (1,tmax-1,1);
    if crep eq 1;
        ct[crep,.] = (vec((zeros(tmax-1,crep)~eye(eyerep))'))';
    else;
        ctemp1 = zeros(tmax-1,crep);
        ctemp2 = eye(eyerep);
        ctemp3 = zeros(tmax-1-eyerep,eyerep);
        ctemp4 = ctemp2|ctemp3;
        ctemp5 = ctemp1~ctemp4;
        ct[crep,.] = (vec((ctemp5')))';
    endif;
    eyerep = eyerep - 1;
endfor;

c1  = ones(1,tmax-1) .*. eye(tmax);
rt  = c1|ct;
b   = zeros(cols(rt'),1);


coset;
_co_dirTol      = 0.0001;
_co_maxiters    = 50;
_co_A           = rt;
_co_B           = b;
_co_Algorithm   = 1;
__output        = 0;
__title         = "Wansbeek-Bekker A";


if sumc(tvx') ge 1;

    xwbtmp  = selif(x',tvx')';
    xwb     = delif(ylag~xwbtmp,e); @ no constant in XWB @

    if sumc(ivs .* timev') ge 1;
        zwbatmp = selif(x',exog)';
        zwbatmp = zwbatmp[.,2:cols(zwbatmp)];     @ zwbtmp are the exog X's @
        zwbatmp = delif(zwbatmp,e);
    else;
        zwbatmp = xwb;
    endif;


elseif sumc(tvx') eq 0;

    xwb     = delif(ylag,e);
    zwbatmp = 0;

endif;

ywb = delif(y,e);

save xwb, zwbatmp;

if ismiss(zwbatmp) eq 1;
	print;
	"WARNING - MISSING VALUES IN W-B Z'S";
endif;

{ ahat,q,g,retcode } = co(&wbiva,aivec);


aia     = (reshape(ahat,tmax-1,tmax))';
aawb    = eye(n) .*. aia;

if sumc(timev' .* xvars') ge 1;

    if sumc(ivs .* timev') ge 1;
        zwbatmp = selif(x',timev' .* exog)';
    else;
        zwbatmp = xwbtmp;
    endif;

    zwba    = (aawb'y)~delif(zwbatmp,e);

else;

    zwba    = (aawb'y);

endif;



/* CONSTRUCT ESTIMATOR */

pzwba   = zwba*inv(zwba'zwba)*zwba';
wba     = inv(xwb'pzwba*xwb)*xwb'pzwba*ywb;



if ismiss(inv(xwb'pzwba*xwb)) eq 1;

    wba     = zeros(cols(xwb,1));
    sewba   = wba;
    twba    = wba;
    R2_wba  = 0;

    wbasarg = -999;
    wbadfsrg= -999;
    wbapsarg= -999;

    wbaid   = wbasarg~wbadfsrg~wbapsarg;

else;

    uwba    = ywb-(xwb*wba);


    /* ESTIMATE COVARIANCE MATRIX */

    colswb  = cols(zwba);
    rowswb  = rows(zwba);
    zhatwba = zeros(colswb,colswb);

    irepst = 1;
    do while irepst le rowswb;
        irepfn  = irepst + (rowswb/n) - 1;
        uwbai   = uwba[irepst:irepfn];
        z1i     = zwba[irepst:irepfn,.];
        if irepst eq 1;
            zhatwba = z1i'*uwbai*uwbai'*z1i;
        elseif irepst gt 1;
            zhatwba = zhatwba + (z1i'*uwbai*uwbai'*z1i);
        endif;
        irepst = irepfn + 1;
    endo;


    if cols(zwba) eq cols(xwb);
        sewba = (abs(diag(inv(zwba'xwb)*zhatwba*inv(xwb'zwba))))^0.5;
    else;
        sewba = abs((diag(inv(xwb'zwba*inv(zhatwba)*zwba'xwb))))^0.5;
    endif;

    twba        = wba ./ sewba;
    Rwba        = 1 - ((uwba'uwba) /
                    sumc((ywb - meanc(ywb))^2));
    R2_wba     = 1 - (((rows(uwba)-1) / (rows(uwba)-rows(wba)))
                    * (1 - Rwba));

    clear v;
    for irep (1, n, 1);

        sel     = delif(idvec,e) .eq irep;
        utmp    = selif(uwba,sel);			
        ztmp    = selif(zwba,sel);   
        v       = v + (ztmp'utmp*utmp'ztmp);
            
    endfor;

    wbasarg     = uwba'zwba*inv(v)*zwba'uwba;
    wbadfsrg    = cols(zwba)-cols(xwb);
    wbapsarg    = cdfchic( wbasarg,wbadfsrg );

    wbaid       = wbasarg~wbadfsrg~wbapsarg;

endif;


/* WANSBEEK-BEKKER ESTIMATOR: A'Y+ & A'X AS INSTRUMENTS: WBB */

if sumc(timev' .* exog) eq 0;

    wbb     = zeros(cols(xwb),1);
    sewbb   = wbb;
    twbb    = wbb;
    R2_wbb  = 0;

    wbbsarg = -999;
    wbbdfsrg= -999;
    wbbpsarg= -999;

    wbbid   = wbbsarg~wbbdfsrg~wbbpsarg;

    goto finish;

endif;


coset;
_co_dirTol      = 0.0001;
_co_maxiters    = 50;
_co_A           = rt;
_co_B           = b;
_co_Algorithm   = 1;
__output        = 0;
__title         = "Wansbeek-Bekker B";


{ ahat,f,g,retcode } = co(&wbivb,aivec);

aib     = (reshape(ahat,tmax-1,tmax))';
abwb    = eye(n) .*. aia;


if sumc(timev' .* xvars') ge 1;

    if sumc(ivs .* timev') ge 1;
        zwbatmp = selif(x',timev' .* exog)';
    else;
        zwbatmp = xwbtmp;
    endif;

    zwba    = (aawb'y)~delif(zwbatmp,e);

else;

    zwba    = (aawb'y);

endif;

zwbb = (abwb'*y)~(abwb'*xfull)~xwb[.,2:cols(xwb)];

if ismiss(zwbb) eq 1;
	print;
	"WARNING - MISSING VALUES IN W-B B'S Z'S";
endif;


/* CONSTRUCT ESTIMATOR */

pzwbb = zwbb*inv(zwbb'*zwbb)*zwbb';
wbb = inv(xwb'*pzwbb*xwb)*xwb'*pzwbb*ywb;

if ismiss(inv(xwb'pzwbb*xwb)) eq 1;

    wbb     = zeros(cols(xwb,1));
    sewbb   = wbb;
    twbb    = wbb;
    R2_wbb  = 0;

    wbbsarg = -999;
    wbbdfsrg= -999;
    wbbpsarg= -999;

    wbbid   = wbbsarg~wbbdfsrg~wbbpsarg;

else;

    uwbb    = ywb-(xwb*wbb);

    /* ESTIMATE COVARIANCE MATRIX */

    colswb  = cols(zwbb);
    rowswb  = rows(zwbb);
    zhatwbb = zeros(colswb,colswb);

    irepst = 1;
    do while irepst le rowswb;
        irepfn  = irepst + (rowswb/n) - 1;
        uwbbi   = uwbb[irepst:irepfn];
        z1i     = zwbb[irepst:irepfn,.];
        if irepst eq 1;
            zhatwbb = z1i'*uwbbi*uwbbi'*z1i;
        elseif irepst gt 1;
            zhatwbb = zhatwbb + (z1i'*uwbbi*uwbbi'*z1i);
        endif;
        irepst = irepfn + 1;
    endo;

    sewbb = abs((diag(inv(xwb'zwbb*inv(zhatwbb)*zwbb'xwb))))^0.5;
    twbb        = wbb ./ sewbb;
    Rwbb        = 1 - ((uwbb'uwbb) /
                    sumc((ywb - meanc(ywb))^2));
    R2_wbb     = 1 - (((rows(uwbb)-1) / (rows(uwbb)-rows(wbb)))
                    * (1 - Rwbb));

    clear v;
    for irep (1, n, 1);

        sel     = delif(idvec,e) .eq irep;
        utmp    = selif(uwbb,sel);			
        ztmp    = selif(zwbb,sel);   
        v       = v + (ztmp'utmp*utmp'ztmp);
            
    endfor;

    wbbsarg     = uwbb'zwbb*inv(v)*zwbb'uwbb;
    wbbdfsrg    = cols(zwbb)-cols(xwb);
    wbbpsarg    = cdfchic( wbbsarg,wbbdfsrg );

    wbbid       = wbbsarg~wbbdfsrg~wbbpsarg;

endif;

goto finish;

gmmproc:

            /* GMM ESTIMATORS: STARTING VALUES */

if ismiss(xt0) eq 1;
	print;
	"WARNING - MISSING VALUES IN GMM INITIAL VALUES OF X";		
endif;


jgmm        = sumc(gmm_con); 
whichgmm    = 0;
save whichgmm;

            /* GMM: IDENTITY WEIGHTING MATRIX */

coset;
_co_dirTol = 0.0001;
_co_maxiters = 50;
_co_Bounds = { -100 100};
/* _co_LineSearch = 1; */
__title         = "GMM Identity";
__output = 1;


{ gmmi,qi,g,retcode } = co(&_gmm,ols);

if qi lt 0;
   print;"IDENTITY GMM FUNCTION FAILED";;stop;
endif;

if retcode eq 6;
   print;"IDENTITY GMM FUNCTION FAILED";;stop;
endif;

/* DEFINE COMPOSITE RESIDUAL */

vit = _vit(gmmi);

vitmat = reshape(vit,n,tmax-1)';    /* vitmat is (T-1)*N */
vitmatt = vitmat';              /* vitmatat is N*(T-1) */

                    /* ESTIMATE W */

/* 1ST CONDITION: E(V_IT) = 0 */

sel_con = 0;


if gmm_con[1,1] eq 1;
    m1a2i   = vitmat;
    sel_con = sel_con|ones(rows(m1a2i),1);

    trap 1;
    m1a2c   = ismiss(inv((m1a2i*m1a2i')./n));
    cls;
    if m1a2c eq 1;
        ?;"ERROR - FIRST CONDITION APPEARS INAPPROPRIATE - DELETE";?;
    endif;

    trap 0;

else;

    m1a2i   = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 2ND CONDITION: E(V_IT^2) = SIG^2_A + SIG^2_U */


if gmm_con[1,2] eq 1;

    m3i     = _m3i(vitmat,tmax,n);
    sel_con = sel_con|ones(rows(m3i),1);

    trap 1;
    m3c   = ismiss(inv((m3i*m3i')./n));
    cls;
    if m3c eq 1;
        ?;"ERROR - SECOND CONDITION APPEARS INAPPROPRIATE - DELETE";?;
    endif;

    trap 0;

else;

    m3i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 3RD CONDITION: E(V_IT^2)(V_IS) = SIG^2_A */

if gmm_con[1,3] eq 1;

    m5i = _m5i(vitmat,tmax,n);
    sel_con = sel_con|ones(rows(m5i),1);

    trap 1;
    m5c   = ismiss(inv((m5i*m5i')./n));
    cls;
    if m5c eq 1;
        ?;"ERROR - THIRD CONDITION APPEARS INAPPROPRIATE - DELETE";?;
    endif;

    trap 0;

else;

    m5i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


/* 4TH CONDIT: E[(V_IT - V_IT-2)(Y_I,T-L)] = 0 */

if gmm_con[1,4] eq 1;

    delta_v = vitmat[2:tmax-1,1:n] - vitmat[1:tmax-2,1:n];
    ymat    = reshape(y,n,tmax)';

    m6rep = 1;
    do while m6rep le tmax-2;

        y_i         = ymat[m6rep,.] .*. ones(tmax-m6rep-1,1);
        deltav_i    = delta_v[m6rep:tmax-2,.];

        if m6rep eq 1;
            m6i         = _m6i(deltav_i,y_i);
        else;
            m6i = m6i|_m6i(deltav_i,y_i);
        endif;

    m6rep = m6rep + 1;
    endo;

    sel_con = sel_con|ones(rows(m6i),1);

    trap 1;
    m6c   = ismiss(inv((m6i*m6i')./n));
    cls;
    if m6c eq 1;
        ?;"ERROR - FOURTH CONDITION APPEARS INAPPROPRIATE - DELETE";?;
    endif;

    trap 0;

else;

    m6i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 5TH CONDITION: E(V_IT)(X_IT) = 0 */

if gmm_con[1,5] eq 1;

    colsm7 = sumc(tvexogx');

    if colsm7 ge 1;

        colrep  = 1;

        xm7     = selif(xgmm',tvexogx')';

        do while colrep le colsm7;

            xtmp    = xm7[.,colrep];       @ DOES NOT INCLUDE LAGGED Y @

            if colrep eq 1;
                m7i  = _m7i(vitmat,xtmp,tmax,n);
            elseif colrep gt 1;
                m7i  = m7i|_m7i(vitmat,xtmp,tmax,n);
            endif;


        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m7i),1);

        trap 1;
        m7c   = ismiss(inv((m7i*m7i')./n));

        if m7c eq 1;
            ?;"ERROR - FIFTH CONDITION APPEARS INAPPROPRIATE - DELETE";?;
        endif;

        trap 0;

    else;

        m7i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m7i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 6TH CONDITION: E(V_IT)(X_I0) = 0 */

if gmm_con[1,6] eq 1;

    load xt0;

    xt0m11 = selif(xt0',gmm6')';

    colsm11 = sumc(gmm6');

    if colsm11 ge 1;

        colrep = 1;
        do while colrep le colsm11;

            xt0tmp = xt0m11[.,colrep];

            if colrep eq 1;
                m11i = _m11i(vitmat,xt0tmp,tmax);
            elseif colrep gt 1;
                m11i = m11i|_m11i(vitmat,xt0tmp,tmax);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m11i),1);

        trap 1;
        m11c   = ismiss(inv((m11i*m11i')./n));
        if m11c eq 1;
            ?;"ERROR - SIXTH CONDITION APPEARS INAPPROPRIATE - DELETE";?;
        endif;

        trap 0;

    else;

        m11i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m11i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 7TH CONDITION: E(V_IT)(X_IS) = 0 */

if gmm_con[1,7] eq 1;

    x13     = selif(xgmm',gmm7')';
    colsm13 = sumc(gmm7');

    if colsm13 ge 1;
        colrep = 1;
        do while colrep le colsm13;

            xtmp    = x13[.,colrep];

            if colrep eq 1;
                m13i = _m13i(vitmat,xtmp,tmax,n);
            elseif colrep gt 1;
                m13i = m13i|_m13i(vitmat,xtmp,tmax,n);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m13i),1);

        trap 1;
        m13c   = ismiss(inv((m13i*m13i')./n));
        if m13c eq 1;
            ?;"ERROR - SEVENTH CONDITION APPEARS INAPPROPRIATE - DELETE";?;
        endif;

        trap 0;

    else;

        m13i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m13i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

sel_con = sel_con[2:rows(sel_con)];

mi = m1a2i|m3i|m5i|m6i|m7i|m11i|m13i;
mi = selif(mi,sel_con);


what = (mi*mi')./n;

trap 1;
checkall = ismiss(invpd(what));
if checkall eq 1;
    ?;"ERROR - OPTIMAL GMM COVARIANCE MATRIX NOT POSITIVE DEFINITE";
    "RE-SPECIFY GMM CONDITIONS";stop;
endif;

trap 0;

                /* DEFINE COVARIANCE MATRIX */

/* DEFINE TIME AVERAGES OF VARIABLES */

avy     = 0;
avx     = zeros(1,cols(xols));
avzgmm  = zeros(1,cols(xgmm));
trep    = 2;
do while trep le tmax;
    dum         = (tvec1 .eq trep);
    avytmp      = selif(yols,dum);
    avy         = avy|meanc(avytmp);
    avxtmp      = selif(xols,dum);
    avx         = avx|meanc(avxtmp)';
    avzgmtmp    = selif(xgmm,dum);
    avzgmm      = avzgmm|meanc(avzgmtmp)';
trep = trep + 1;
endo;

avy     = avy[2:rows(avy)];
avx     = avx[2:rows(avx),.];
avu     = avy - (avx*gmmi);
avzgmm  = avzgmm[2:rows(avzgmm),.];

/* CONDITION 1 BLOCK: E(V_IT) = 0 */

sel_con = 0;

if gmm_con[1,1] eq 1;
    m1a2g   = _m1a2g(avx);
    sel_con = sel_con|ones(rows(m1a2g),1);
else;
    m1a2g = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;


/* 2ND CONDITION: E(V_IT^2) = SIG^2_A + SIG^2_U */

if gmm_con[1,2] eq 1;
    m3g = _m3g(avu,avx,tmax);
    sel_con = sel_con|ones(rows(m3g),1);
else;
    m3g = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;

/* 3RD CONDITION: E(V_IT^2)(V_IS) = SIG^2_A */

if gmm_con[1,3] eq 1;
    m5g = _m5g(avu,avx,tmax);
    sel_con = sel_con|ones(rows(m5g),1);
else;
    m5g     = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;

/* 4TH CONDIT: E[(V_IT - V_I,T-1)(Y_I,T-L)] = 0   */

if gmm_con[1,4] eq 1;
    m6g = _m6g(avy,avx,tmax);
    sel_con = sel_con|ones(rows(m6g),1);
else;
    m6g     = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;

/* 5TH CONDITION: E(V_IT)(Z_IT) = 0 */

if gmm_con[1,5] eq 1;

    colsm7  = sumc(tvexogx');

    if colsm7 ge 1;

        avzm7   = selif(avzgmm',tvexogx')';

        colrep = 1;
        do while colrep le colsm7;

            ztmp    = avzm7[.,colrep];

            if colrep eq 1;
                m7g  = _m7g(avx,ztmp);
            elseif colrep gt 1;
                m7g  = m7g|_m7g(avx,ztmp);
            endif;
        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m7g),1);

    else;

        m7g     = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m7g     = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;

/* 6TH CONDITION: E(V_IT)(X_I0) = 0 */

if gmm_con[1,6] eq 1;

    load xt0;
    xt0m11 = selif(xt0',gmm6')';

    colsm11 = sumc(gmm6');

    if colsm11 ge 1;

        colrep = 1;
        do while colrep le colsm11;

            xt0tmp = xt0m11[.,colrep];

            if colrep eq 1;
                m11g = _m11g(avx,xt0tmp);
            elseif colrep gt 1;
                m11g = m11g|_m11g(avx,xt0tmp);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m11g),1);

    else;

        m11g    = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m11g    = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;

/* 7TH CONDITION: E(V_IT)(X_IS) = 0 */

if gmm_con[1,7] eq 1;

    x13     = selif(avzgmm',gmm7')';
    colsm13 = sumc(gmm7');

    if colsm13 ge 1;

        colrep = 1;
        do while colrep le colsm13;

            ztmp    = x13[.,colrep];

            if colrep eq 1;
                m13g = _m13g(avx,ztmp,tmax);
            elseif colrep gt 1;
                m13g = m13g|_m13g(avx,ztmp,tmax);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m13g),1);

    else;

        m13g     = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m13g     = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;


        /* FINALLY CONSTRUCT COVARIANCE MATRIX */

sel_con = sel_con[2:rows(sel_con)];

capM    = m1a2g|m3g|m5g|m6g|m7g|m11g|m13g;
capG    = selif(capM,sel_con);


trap 1;
vargmmi = inv(capG'eye(rows(capG))*capG) *
            (capG'eye(rows(capG))*what*eye(rows(capG))*capG) *
            inv(capG'eye(rows(capG))*capG);

if ismiss(vargmmi) eq 1;
    vargmmi = ones(rows(gmmi),1) .* -999;
    segmmi  = vargmmi;
    tgmmi   = vargmmi;
else;               
    vargmmi = vargmmi ./n;
    segmmi  = (abs(diag(vargmmi)))^0.5;
    tgmmi   = gmmi./segmmi;
endif;    

            /* CONSTRUCT R-SQUARED FOR GMM_I */


if ismiss(gmmi~segmmi~tgmmi) eq 1;
    gmmi    = zeros(cols(xols),1);
    segmmi  = gmmi;
    tgmmi   = gmmi;
    R2_gmmi = -999;
else;
    ugmmi       = _vit(gmmi);
    sig2_gmi    = (ugmmi'*ugmmi)/(rows(ugmmi)-rows(gmmi));
    Rgmmi       = 1 - ((ugmmi'ugmmi) / sumc((yols - meanc(yols))^2));
    R2_gmmi     = 1 - (((rows(ugmmi)-1) / (rows(ugmmi)-rows(gmmi)))
                    * (1 - Rgmmi));
endif;


/* CONSTRUCT HAUSMAN TEST FOR OVER-IDENTIFYING RESTRICTIONS */

load jgmm;

dfgmm   = jgmm - rows(gmmi);
pwaldi  = 1-cdfchic(qi,dfgmm);
overidi = qi~dfgmm~pwaldi;

                /* GMM ON ALL CONDITIONS */

whichgmm    = 1;
save whichgmm;

coset;
_co_dirTol      = 0.001;
_co_maxiters    = 75;
_co_Bounds      = { -100 100};
/* _co_LineSearch = 1; */
__output        = 1;
__title         = "GMM Optimal";



{ gmm,qw,g,retcode } = co(&_gmm,gmmi);


if qw lt 0;
    print;"OPTIMAL GMM FUNCTION FAILED";;stop;
endif;

if retcode eq 6;
    print;"OPTIMAL GMM FUNCTION FAILED";;stop;
endif;

if gmm_se eq 1;
    infomat = _co_FinalHess;
    infomat = inv(infomat);
    segmm  = sqrt(abs(diag(infomat)));
    tgmm   = gmm./segmm;
    goto skipse;
endif;

                    /* RE-ESTIMATE W */

/* DEFINE COMPOSITE RESIDUAL */

vit = _vit(gmm);

vitmat = reshape(vit,n,tmax-1)';    /* vitmat is (T-1)*N */
vitmatt = vitmat';              /* vitmatat is N*(T-1) */


/* 1ST CONDITION: E(V_IT) = 0 */

sel_con = 0;


if gmm_con[1,1] eq 1;
    m1a2i   = vitmat;
    sel_con = sel_con|ones(rows(m1a2i),1);

    trap 1;
    m1a2c   = ismiss(inv((m1a2i*m1a2i')./n));
    cls;
    if m1a2c eq 1;
        "ERROR - FIRST CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
    endif;

    trap 0;

else;

    m1a2i   = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 2ND CONDITION: E(V_IT^2) = SIG^2_A + SIG^2_U */

if gmm_con[1,2] eq 1;

    m3i     = _m3i(vitmat,tmax,n);
    sel_con = sel_con|ones(rows(m3i),1);

    trap 1;
    m3c   = ismiss(inv((m3i*m3i')./n));
    cls;
    if m3c eq 1;
        "ERROR - SECOND CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
    endif;

    trap 0;

else;

    m3i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


/* 3RD CONDITION: E(V_IT^2)(V_IS) = SIG^2_A */

if gmm_con[1,3] eq 1;

    m5i = _m5i(vitmat,tmax,n);
    sel_con = sel_con|ones(rows(m5i),1);

    trap 1;
    m5c   = ismiss(inv((m5i*m5i')./n));
    cls;
    if m5c eq 1;
        "ERROR - THIRD CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
    endif;

    trap 0;

else;

    m5i     = zeros(1,n);
    sel_con = sel_con|0;

endif;

/* 4TH CONDIT: E(V_IT)(Y_I0) = CONSTANT   */

if gmm_con[1,4] eq 1;


    delta_v = vitmat[2:tmax-1,1:n] - vitmat[1:tmax-2,1:n];
    ymat    = reshape(y,n,tmax)';

    m6rep = 1;
    do while m6rep le tmax-2;

        y_i         = ymat[m6rep,.] .*. ones(tmax-m6rep-1,1);
        deltav_i    = delta_v[m6rep:tmax-2,.];

        if m6rep eq 1;
            m6i         = _m6i(deltav_i,y_i);
        else;
            m6i = m6i|_m6i(deltav_i,y_i);
        endif;

    m6rep = m6rep + 1;
    endo;

    sel_con = sel_con|ones(rows(m6i),1);

    trap 1;
    m6c   = ismiss(inv((m6i*m6i')./n));
    cls;
    if m6c eq 1;
        "ERROR - FOURTH CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
    endif;

    trap 0;

else;

    m6i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


/* 5TH CONDITION: E(V_IT)(X_IT) = 0 */

if gmm_con[1,5] eq 1;

    colsm7 = sumc(tvexogx');

    if colsm7 ge 1;

        xm7     = selif(xgmm',tvexogx')';

        colrep  = 1;

        do while colrep le colsm7;

            xtmp    = xm7[.,colrep];       @ DOES NOT INCLUDE LAGGED Y @

            if colrep eq 1;
                m7i  = _m7i(vitmat,xtmp,tmax,n);
            elseif colrep gt 1;
                m7i  = m7i|_m7i(vitmat,xtmp,tmax,n);
            endif;


        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m7i),1);

        trap 1;
        m7c   = ismiss(inv((m7i*m7i')./n));
        if m7c eq 1;
            cls;"ERROR - FIFTH CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
        endif;

        trap 0;

    else;

        m7i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m7i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


/* 6TH CONDITION: E(V_IT)(X_I0) = 0 */

if gmm_con[1,6] eq 1;

    load xt0;

    xt0m11 = selif(xt0',gmm6')';

    colsm11 = sumc(gmm6');

    if colsm11 ge 1;

        colrep = 1;
        do while colrep le colsm11;

            xt0tmp = xt0m11[.,colrep];

            if colrep eq 1;
                m11i = _m11i(vitmat,xt0tmp,tmax);
            elseif colrep gt 1;
                m11i = m11i|_m11i(vitmat,xt0tmp,tmax);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m11i),1);

        trap 1;
        m11c   = ismiss(inv((m11i*m11i')./n));
        if m11c eq 1;
            cls;"ERROR - SIXTH CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
        endif;

        trap 0;

    else;

        m11i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m11i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


/* 7TH CONDITION: E(V_IT)(X_IS) = 0 */

if gmm_con[1,7] eq 1;

    x13     = selif(xgmm',gmm7')';
    colsm13 = sumc(gmm7');

    if colsm13 ge 1;

        colrep = 1;
        do while colrep le colsm13;

            xtmp    = x13[.,colrep];

            if colrep eq 1;
                m13i = _m13i(vitmat,xtmp,tmax,n);
            elseif colrep gt 1;
                m13i = m13i|_m13i(vitmat,xtmp,tmax,n);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m13i),1);

        trap 1;
        m13c   = ismiss(inv((m13i*m13i')./n));
        if m13c eq 1;
            cls;"ERROR - SEVENTH CONDITION APPEARS INAPPROPRIATE - DELETE";stop;
        endif;

        trap 0;

    else;

        m13i     = zeros(1,n);
        sel_con = sel_con|0;

    endif;

else;

    m13i     = zeros(1,n);
    sel_con = sel_con|0;

endif;


sel_con = sel_con[2:rows(sel_con)];

mi = m1a2i|m3i|m5i|m6i|m7i|m11i|m13i;
mi = selif(mi,sel_con);


what = (mi*mi') ./n;

trap 1;
checkall = ismiss(invpd(what));
cls;
if checkall eq 1;

    "ERROR - OPTIMAL GMM COVARIANCE MATRIX NOT POSITIVE DEFINITE";
    "RE-SPECIFY GMM CONDITIONS";stop;
endif;

trap 0;


        /* CONSTRUCT COVARIANCE MATRIX */

avu     = avy - (avx*gmm);


/* CONDITION 1 BLOCK: E(V_IT) = 0 */

sel_con = 0;

if gmm_con[1,1] eq 1;
    m1a2g   = _m1a2g(avx);
    sel_con = sel_con|ones(rows(m1a2g),1);
else;
    m1a2g   = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;


/* 2ND CONDITION: E(V_IT^2) = SIG^2_A + SIG^2_U */

if gmm_con[1,2] eq 1;
    m3g     = _m3g(avu,avx,tmax);
    sel_con = sel_con|ones(rows(m3g),1);
else;
    m3g     = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;

/* 3RD CONDITION: E(V_IT^2)(V_IS) = SIG^2_A */

if gmm_con[1,3] eq 1;
    m5g     = _m5g(avu,avx,tmax);
    sel_con = sel_con|ones(rows(m5g),1);
else;
    m5g     = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;


/* 4TH CONDIT: E[(V_IT - V_I,T_1)(Y_I,T-L)] = 0  */

if gmm_con[1,4] eq 1;
    m6g     = _m6g(avy,avx,tmax);
    sel_con = sel_con|ones(rows(m6g),1);
else;
    m6g     = zeros(1,cols(avx));
    sel_con = sel_con|0;
endif;


/* 5TH CONDITION: E(V_IT)(X_IT) = 0 */

if gmm_con[1,5] eq 1;

    colsm7  = sumc(tvexogx');

    if colsm7 ge 1;

        avzm7   = selif(avzgmm',tvexogx')';

        colrep = 1;
        do while colrep le colsm7;

            ztmp    = avzm7[.,colrep];

            if colrep eq 1;
                m7g  = _m7g(avx,ztmp);
            elseif colrep gt 1;
                m7g  = m7g|_m7g(avx,ztmp);
            endif;
        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m7g),1);

    else;

        m7g     = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m7g     = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;


/* 6TH CONDITION: E(V_IT)(X_I0) = 0 */

if gmm_con[1,6] eq 1;

    load xt0;
    xt0m11 = selif(xt0',gmm6')';

    colsm11 = sumc(gmm6');

    if colsm11 ge 1;

        colrep = 1;
        do while colrep le colsm11;

            xt0tmp = xt0m11[.,colrep];

            if colrep eq 1;
                m11g = _m11g(avx,xt0tmp);
            elseif colrep gt 1;
                m11g = m11g|_m11g(avx,xt0tmp);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m11g),1);

    else;

        m11g    = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m11g    = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;

/* 7TH CONDITION: E(V_IT)(X_IS) = 0 */

if gmm_con[1,7] eq 1;

    x13     = selif(avzgmm',gmm7')';
    colsm13 = sumc(gmm7');

    if colsm13 ge 1;

        colrep = 1;
        do while colrep le colsm13;

            ztmp    = x13[.,colrep];

            if colrep eq 1;
                m13g = _m13g(avx,ztmp,tmax);
            elseif colrep gt 1;
                m13g = m13g|_m13g(avx,ztmp,tmax);
            endif;

        colrep = colrep + 1;
        endo;

        sel_con = sel_con|ones(rows(m13g),1);

    else;

        m13g     = zeros(1,cols(avx));
        sel_con = sel_con|0;

    endif;

else;

    m13g     = zeros(1,cols(avx));
    sel_con = sel_con|0;

endif;

        /* FINALLY CONSTRUCT COVARIANCE MATRIX FOR OPTIMAL GMM */

sel_con = sel_con[2:rows(sel_con)];

capM    = m1a2g|m3g|m5g|m6g|m7g|m11g|m13g;
capG    = selif(capM,sel_con);

trap 1;
/* vargmm = (inv(capG'inv(what)*capG)) ./(n^2); */
vargmm = (inv(capG'inv(what)*capG));
varchk = ismiss(vargmm);
if varchk eq 1;
    print;
    "ERROR - GMM VARAIANCE MATRIX NOT PD - REDEFINE CONDITIONS?";
    vargmm = inv(capG'inv(what)*capG);
endif;
trap 0;

vargmm = vargmm ./n;


segmm  = (diag(vargmm))^0.5;
tgmm   = abs(gmm./segmm);


skipse:

            /* CONSTRUCT R-SQUARED FOR GMM */


if ismiss(gmm~segmm~tgmm) eq 1;
    gmm     = zeros(cols(xols),1);
    segmm   = gmm;
    tgmm    = gmm;
    R2_gmmi = -999;
else;
    ugmm        = _vit(gmm);
    sig2_gmm    = (ugmm'*ugmm)/(rows(ugmm)-rows(gmm));
    Rgmm        = 1 - ((ugmm'ugmm) / sumc((yols - meanc(yols))^2));
    R2_gmm      = 1 - (((rows(ugmm)-1) / (rows(ugmm)-rows(gmm)))
                    * (1 - Rgmm));
endif;

/* CONSTRUCT HAUSMAN TEST FOR OVER-IDENTIFYING RESTRICTIONS */

load jgmm;

dfgmm   = jgmm - rows(gmm);
pwald   = 1-cdfchic(qw,dfgmm);
overid  = qw~dfgmm~pwald;


finish:

/*
within;
sewn;
twn;
R2_wn;
fgls;
sefgls;
tfgls;
R2_fgls;
dols;
sedols;
tdols;
R2_dols;
NTdols;
aha;
seaha;
taha;
R2_aha;
ahaid;
AR2aha;
NTah;
ahb;
seahb;
tahb;
R2_ahb;
ahbid;
AR2ahb;
*/
  
              /* RETURN OUTPUTS */

retp(   ols~seols~tols,R2_ols,NTols,
        within~sewn~twn,R2_wn,
        fgls~sefgls~tfgls,R2_fgls,
        dols~sedols~tdols,R2_dols,NTdols,
        aha~seaha~taha,R2_aha,ahaid,AR2aha,NTah,
        ahb~seahb~tahb,R2_ahb,ahbid,AR2ahb,
        ahc~seahc~tahc,R2_ahc,ahcid,AR2ahc,
        ara~seara~tara,R2_ara,araid,AR2ara,
        arb~searb~tarb,R2_arb,arbid,AR2arb,
        arc~searc~tarc,R2_arc,arcid,AR2arc,
        aba~seaba~taba,R2_aba,abaid,AR2aba,
        abc~seabc~tabc,R2_abc,abcid,AR2abc,
        bb~sebb~tbb,R2_bb,bbid,AR2bb,
        wba~sewba~twba,R2_wba,wbaid,
        wbb~sewbb~twbb,R2_wbb,wbbid,
        gmmi~segmmi~tgmmi,R2_gmmi,overidi,
        gmm~segmm~tgmm,R2_gmm,overid
);
endp;
            /* OPTIMISATION PROCEDURES */


/* WANSBEEK-BEKKER "A" */

proc (1) = wbiva(aivec);
local tmax, n, fnctwba, aawb, aia, zwba, pzwba, timev, xvars, ivs, e; 
local exog, check;
load tmax, n;

    aia     = (reshape(aivec,tmax-1,tmax))';
    aawb    = eye(n) .*. aia;

    if rows(zwbatmp) gt 1;
        zwba    = (aawb'*y)~zwbatmp;
        pzwba   = zwba*inv(zwba'*zwba)*zwba';
        fnctwba = sumc(diag(inv(xwb'*pzwba*xwb)));
    else;
        zwba    = (aawb'*y);
        fnctwba = (inv(zwba'xwb)*zwba'zwba*inv(xwb'zwba));
    endif;


retp(fnctwba);
endp;


/* WANSBEEK-BEKKER "B" */

proc (1) = wbivb(aivec);
local tmax, n, fnctwbb, abwb, aib, zwbb, pzwbb;
load tmax, n;


    aib     = (reshape(aivec,tmax-1,tmax))';
    abwb    = eye(n) .*. aib;


    zwbb    = (abwb'y)~(abwb'xfull)~delif(zwbatmp,e);
    pzwbb   = zwbb*inv(zwbb'*zwbb)*zwbb';
    fnctwbb = sumc(diag(inv(xwb'*pzwbb*xwb)));

retp(fnctwbb);
endp;


        /* VARIOUS SMALL PROCEDURES FOR GMM ESTIMATORS */

proc _vit(beta);
retp(yols - (xols*beta));
endp;

proc _m1a2i(vitmat);
retp(vitmat);
endp;

proc _m3i(vitmat,tmax,n);
retp((vitmat[1:tmax-2,1:n]^2) - (vitmat[2:tmax-1,1:n]^2));
endp;


proc _m5i(vitmat,tmax,n);
local m5itmp, m5irep, vtreplf, vtreprg, m5i;

    m5itmp  = zeros(1,n);
    m5irep  = 1;
    vtreplf = 1;
    vtreprg = 2;
    do while vtreplf le tmax-2;
        do while vtreprg le tmax-1;
            m5itmp   = m5itmp|
                    ((vitmat[vtreplf,1:n]).*(vitmat[vtreprg,1:n]));
            m5irep   = m5irep + 1;
            vtreprg  = vtreprg + 1;
        endo;
        vtreplf  = vtreplf + 1;
        vtreprg  = vtreplf + 1;
    endo;

    m5itmp  = m5itmp[2:rows(m5itmp),.];
    m5i     = m5itmp[1:rows(m5itmp)-1,1:n] - m5itmp[2:rows(m5itmp),1:n];

retp(m5i);
endp;


proc _m6i(deltav_i,y_i);
local m6i;


   m6i = (y_i .* deltav_i);


retp(m6i);
endp;

proc _m7i(vitmat,xtmp,tmax,n);
local m7i;

    xtmp    = reshape(xtmp,n,tmax-1)';
    m7i     = vitmat .* xtmp;

retp(m7i);
endp;

proc _m11i(vitmat,xt0tmp,tmax);
local m11i;

    m11i =  (xt0tmp' .*. ones(tmax-1,1)) .* vitmat;

retp(m11i);
endp;


proc _m13i(vitmat,xtmp,tmax,n);
local m13i, rep5, rep6;

    m13i    = zeros(1,n);
    xtmp    = reshape(xtmp,n,tmax-1)';

    rep5 = 1;
    rep6 = 2;

    do while rep5 le tmax-2;
        do while rep6 le tmax-1;
            m13i = m13i|((vitmat[rep5,1:n]) .* (xtmp[rep6,1:n]));
            rep6 = rep6 + 1;
        endo;
        rep5 = rep5 + 1;
        rep6 = rep5 + 1;
    endo;

    m13i = m13i[2:rows(m13i),.];

retp(m13i);
endp;

proc _m1a2g(avx);
retp(-avx);
endp;

proc _m3g(avu,avx,tmax);
local m3g;

    m3g = (2*((avu[1:tmax-2] .*. ones(1,cols(avx))) .* -avx[1:tmax-2,.])) -
            (2*((avu[2:tmax-1] .*. ones(1,cols(avx))) .* -avx[2:tmax-1,.]));
retp(m3g);
endp;


proc _m5g(avu,avx,tmax);
local m5gtmp, vtreplf, vtreprg, m5g;

    m5gtmp  = zeros(1,cols(avx));
    vtreplf = 1;
    vtreprg = 2;
    do while vtreplf le tmax-2;
        do while vtreprg le tmax-1;
            m5gtmp   = m5gtmp|
                    (((avu[vtreplf]) .* (-avx[vtreprg,.])) +
                    ((avu[vtreprg]) .* (-avx[vtreplf,.])));
            vtreprg  = vtreprg + 1;
        endo;
        vtreplf  = vtreplf + 1;
        vtreprg  = vtreplf + 1;
    endo;

    m5gtmp  = m5gtmp[2:rows(m5gtmp),.];
    m5g     = m5gtmp[1:rows(m5gtmp)-1,.] - m5gtmp[2:rows(m5gtmp),.];

retp(m5g);
endp;


proc _m6g(avy,avx,tmax);
local delta_x, m6rep, m6gtmp, m6g;

    delta_x = avx[2:rows(avx),.] - avx[1:rows(avx)-1,.];

    m6rep = 1;
    do while m6rep le tmax-2;

        m6gtmp = -(avy[m6rep] .*. ones(tmax-m6rep-1,1))
                .* delta_x[m6rep:tmax-2,.];

        if m6rep eq 1;
            m6g = m6gtmp;
        else;
            m6g = m6g|m6gtmp;
        endif;

    m6rep = m6rep + 1;
    endo;

retp(m6g);
endp;

proc _m7g(avx,ztmp);
local m7g;

    m7g = -(avx .* (ztmp .*. ones(1,cols(avx))));

retp(m7g);
endp;

proc _m11g(avx,xt0tmp);
local m11g;

    xt0tmp  = meanc(xt0tmp);
    m11g    =  -(avx .*xt0tmp);

retp(m11g);
endp;

proc _m13g(avx,ztmp,tmax);
local m13g, rep5, rep6;

    m13g    = zeros(1,cols(avx));

    rep5 = 1;
    rep6 = 2;

    do while rep5 le tmax-2;
        do while rep6 le tmax-1;
            m13g = m13g|((avx[rep5,.]) .* (ztmp[rep6,.]));
            rep6 = rep6 + 1;
        endo;
        rep5 = rep5 + 1;
        rep6 = rep5 + 1;
    endo;

    m13g = m13g[2:rows(m13g),.];


retp(m13g);
endp;


/* GMM PROCEDURE */

proc (1) = _gmm(beta);
local fnct,m,m1a2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14;
local m1a2i,m3i,m4i,m5itmp,m5i,m6i,m7i,m8i,m9i,m10i,m11i,m12i,m13i;
local m3irep,vtreplf,vtreprg,vit,vitmat,vitmatt,vi0,vi0sq;
local m3irows,m5irows,m5irep,rep1,rep2,rep3,rep4,rep5,rep6;
local m6irows,m6itmp,m6irep,m13irep;
local tmax, n, sel_con, yt0, colsm7, xm7, exog, colrep, xtmp;
local xt0, xt0m11, colsm11, xt0tmp, gmm6, gmm7, colsm13, x13, jgmm, whichgmm;
local gmm_con, delta_v, ymat, m6rep, y_i, deltav_i, tvexogx;


load tmax, n, gmm_con, yt0, exog, xt0, gmm6, gmm7, whichgmm, tvexogx;

    sel_con = 0;

    /* DEFINE COMPOSITE RESIDUAL */

    vit = _vit(beta);

    vitmat = reshape(vit,n,tmax-1)';    /* vitmat is (T-1)*N */
    vitmatt = vitmat';              /* vitmatat is N*(T-1) */


    /* 1ST CONDITION: E(V_IT) = 0 */

    if gmm_con[1,1] eq 1;
        m1a2i   = _m1a2i(vitmat);
        m1a2    = meanc(m1a2i');
        sel_con = sel_con|ones(tmax-1,1);
    else;
        m1a2    = 0;
        sel_con = sel_con|0;
    endif;


    /* 2ND CONDITION: E(V_IT^2) = SIG^2_A + SIG^2_U */

    if gmm_con[1,2] eq 1;
        m3i = _m3i(vitmat,tmax,n);
        m3  = meanc(m3i');
        sel_con = sel_con|ones(tmax-2,1);
    else;
        m3      = 0;
        sel_con = sel_con|0;
    endif;


    /* 3RD CONDITION: E(V_IT)(V_IS) = SIG^2_A */

    if gmm_con[1,3] eq 1;
        m5i     = _m5i(vitmat,tmax,n);
        m5      = meanc(m5i');
        sel_con = sel_con|ones(rows(m5),1);
    else;
        m5      = 0;
        sel_con = sel_con|0;
    endif;


    /* 4TH CONDIT: E[(V_IT) - (V_IT-1)] Y_IT-L] = 0; l>=2*/

    if gmm_con[1,4] eq 1;

        delta_v = vitmat[2:tmax-1,1:n] - vitmat[1:tmax-2,1:n];
        ymat    = reshape(y,n,tmax)';

        m6rep = 1;
        do while m6rep le tmax-2;

            y_i         = ymat[m6rep,.] .*. ones(tmax-m6rep-1,1);
            deltav_i    = delta_v[m6rep:tmax-2,.];
            m6i         = _m6i(deltav_i,y_i);

            if m6rep eq 1;
                m6 = meanc(m6i');
            else;
                m6 = m6|meanc(m6i');
            endif;

        m6rep = m6rep + 1;
        endo;

        sel_con = sel_con|ones(rows(m6),1);

    else;

        m6      = 0;
        sel_con = sel_con|0;

    endif;


        /* 5TH CONDITION: E(V_IT)(X_IT) = 0 */

    if gmm_con[1,5] eq 1;

        colsm7 = sumc(tvexogx');

        if colsm7 ge 1;

            xm7    = selif(xgmm',tvexogx')';

            colrep = 1;
            do while colrep le colsm7;

                xtmp    = xm7[.,colrep];

                m7i  = _m7i(vitmat,xtmp,tmax,n);

                if colrep eq 1;
                    m7  = meanc(m7i');
                elseif colrep gt 1;
                    m7  = m7|meanc(m7i');
                endif;

            colrep = colrep + 1;
            endo;

            sel_con = sel_con|ones(rows(m7),1);

        else;

            m7      = 0;
            sel_con = sel_con|0;

        endif;

    else;

        m7      = 0;
        sel_con = sel_con|0;

    endif;


        /* 6TH CONDITION: E(V_IT)(X_I0) = 0 */

    if gmm_con[1,6] eq 1;

        colsm11 = sumc(gmm6');

        if colsm11 ge 1;
            colrep = 1;

            xt0m11 = selif(xt0',gmm6')';

            do while colrep le colsm11;

                xt0tmp = xt0m11[.,colrep];

                m11i = _m11i(vitmat,xt0tmp,tmax);

                if colrep eq 1;
                    m11  = meanc(m11i');
                elseif colrep gt 1;
                    m11  = m11|meanc(m11i');
                endif;

            colrep = colrep + 1;
            endo;

            sel_con = sel_con|ones(rows(m11),1);

        else;

            m11     = 0;
            sel_con = sel_con|zeros(rows(m11),1);

        endif;

    else;

        m11     = 0;
        sel_con = sel_con|zeros(rows(m11),1);

    endif;


        /* 7TH CONDITION: E(V_IT)(X_IS) = 0 */

    if gmm_con[1,7] eq 1;

        x13     = selif(xgmm',gmm7')';
        colsm13 = sumc(gmm7');

        if colsm13 ge 1;

            colrep = 1;
            do while colrep le colsm13;

                xtmp = x13[.,colrep];
                m13i = _m13i(vitmat,xtmp,tmax,n);

                if colrep eq 1;
                    m13 = meanc(m13i');
                else;
                    m13 = m13|meanc(m13i');
                endif;

            colrep = colrep + 1;
            endo;

            sel_con = sel_con|ones(rows(m13),1);

        else;

            m13     = 0;
            sel_con = sel_con|zeros(rows(m13),1);

        endif;

    else;

        m13     = 0;
        sel_con = sel_con|zeros(rows(m13),1);

    endif;


sel_con = sel_con[2:rows(sel_con)];

m = m1a2|m3|m5|m6|m7|m11|m13;
m = selif(m,sel_con);

jgmm = rows(m);
save jgmm;

if jgmm le rows(beta);
    cls;"GMM SYSTEM EXACTLY/UNDERIDENTIFIED - ADD CONDITIONS";
    "# OF PARAMETERS";;rows(beta);
    "# OF MOMENT CONDITIONS";;jgmm;stop;
endif;

if whichgmm eq 0;
    fnct = n*(m'*eye(rows(m))*m);
elseif whichgmm eq 1;
    fnct = n*(m'*invpd(what)*m);
endif;

retp(fnct);
endp;


