/***************************************************
    ** ESTIMATION OF SYMPLE DYNAMIC PANEL DATA MODELS **
    **              WRITTEN BY S. C. AHN              **
    ***************************************************/

@ CLEAN MEMMORY @
    new  ;

@ OPEN OUTPUT FILE @
    output file = dynam.out reset ;

@ FORMAT OUTPUT FILE @
    format /rd 10,3 ;

@ LOADING DATA @
	load dat[500,6] = dynam.db ;

@ Define yy @
	yy = dat;

@ Define N and T @
	nn = 500;
	tt = 5;


/* From here, do not make any change */

@ CREATING DIFFERENCED DATA: DYY INCLUDES (Y0-Y1),...,(YT-1-YT) @
    dyy = yy[.,2]-yy[.,1] ;
    j  = 2  ; do while j <= tt ;
        dyy = dyy~(yy[.,j+1]-yy[.,j]) ;
    j  = j+1; endo ;

@ STACKING DATA @
    depy  = vec(yy[.,2:tt+1]')   ; @ Y @
    lagy  = vec(yy[.,1:tt]')     ; @ Y_{-1} @
    pdepy = meanc(yy[.,2:tt+1]') .*. ones(tt,1) ; @ P_v Y @
    plagy = meanc(yy[.,1:tt]' ) .*. ones(tt,1) ;  @ P_v Y_{-1} @
    qdepy = depy - pdepy ;  @ Q_v Y @
    qlagy = lagy - plagy ;  @ Q_v Y_{-1} @
    ddepy = vec(dyy[.,2:tt]');
    dlagy = vec(dyy[.,1:tt-1]');
    abwdy = ddepy ;
    abwly = vec(yy[.,2:tt]') ;
    abwdly= dlagy            ;

@ CREATING ARELLANO AND BOND, ARELLANO-BOVER IVS @

@ CREATING L(T,P) MATRIX: NEEDED TO CREATE AB IVS @
    proc ll(t,p) ;
    local i,u    ;

        u = zeros(t,p) ;
        i = 1 ; do while i <= p ;
            u[t-p+i-1,i] = - 1  ; u[t-p+i,i] = 1  ;
        i = i + 1 ; endo        ;

    retp(u) ; endp ;

@ CREATING AHLL(T,P) MATRIX: NEEDED TO CREATE DIAGONAL ANDERSON-HSIAO IVS @
    proc ahll(t,p) ;
   local i,u ;

        u = zeros(t,1) ;
        u[p+2] = - 1  ; u[p+1] = 1  ;

    retp(u) ; endp ;

@ CREATING AB INSTRUMENTS @
    abiv = zeros(tt*nn,1) ;
    j = 1     ; do while j <= tt-1 ;
        abiv = abiv~( yy[.,j] .*. ll(tt,tt-j) ) ;
    j = j + 1 ; endo ;

    pp   = tt*(tt-1)/2    ;
    abiv = abiv[.,2:pp+1] ;

@ CREATING SUB-AB (Diagonal form of Andersen-Hsiao) INSTRUMENTS @
    ahiv = zeros(tt*nn,1) ;
    j = 1     ; do while j <= tt-1 ;
        ahiv = ahiv~( yy[.,j] .*. ahll(tt,j-1) ) ;
    j = j + 1 ; endo ;

    ahpp=tt-1 ;
    ahiv = ahiv[.,2:ahpp+1]  ;

@ CREATING ABo INSTRUMENTS @
    aboiv = zeros(tt*nn,1) ;
    j = 1     ; do while j <= tt-1 ;
      mmmm = zeros(tt,1) ;
      mmmm[j+1] = 1      ;
      aboiv = aboiv~( dyy[.,j] .*. mmmm ) ;
    j = j + 1 ; endo ;

    aboiv = aboiv[.,2:cols(aboiv)]  ;

@ CREATING AB-Type ABo (AABo) INSTRUMENTS @
    aaboiv = zeros(tt*nn,1) ;
    j = 1     ; do while j <= tt-1 ;
      mmmm = zeros(tt,1) ;
      mmmm[j+1] = 1      ;
      aaboiv = aaboiv~( dyy[.,1:j] .*. mmmm ) ;
    j = j + 1 ; endo ;

    aaboiv = aaboiv[.,2:cols(aaboiv)] ;

@ CREATING DIAGONAL LAGGED Y @
    diagly = zeros(tt*nn,1) ;
    j = 1 ; do while j <= tt ;
        mmmm = zeros(tt,1) ;
        mmmm[j] = 1        ;
        diagly = diagly~(yy[.,j] .*. mmmm) ;
    j = j + 1 ; endo  ;
        g1pp   = cols(diagly) ;
        diagly = diagly[.,2:g1pp] ;

@ # OF AB INSTRUMENTS @
    nabiv  = pp ;

@ # OF AH INSTRUMENTS @
    nahiv  = cols(ahiv) ;

/**********************
** ARELLANO-BOND GMM **
**********************/

@ TSLS ESTIMATOR @

    tscov = invpd( (lagy'abiv)*invpd(abiv'abiv)*(abiv'lagy) ) ;
    tsb   = tscov*(lagy'abiv)*invpd(abiv'abiv)*(abiv'depy)    ;
    qres  = qdepy - qlagy*tsb                                 ;
    tssige2 = (qres'qres/(nn*tt-tt-1)) ;
    tscov   = tssige2*tscov            ;

res  = depy - lagy*tsb  ;

opivab = zeros(1,cols(abiv)) ;
j = 0 ; do while j <= nn-1 ;
    zz = res[tt*j+1:tt*j+tt,.]'abiv[tt*j+1:tt*j+tt,.] ;
    opivab = opivab|zz ;
j = j + 1 ; endo ; opivab = opivab[2:nn+1,.] ;

optwab = invpd(opivab'opivab)    ;


abv = INVPD( (lagy'abiv)*(optwab)*(abiv'lagy) ) ;
abb = abv*(lagy'abiv)*(optwab)*(abiv'depy)      ;

abe  = depy - lagy*abb ;
abj  = (abe'abiv)*optwab*(abiv'abe);
abdf = cols(abiv) - rows(abb);


/***********************
** ARELLANO-BOVER GMM **
***********************/

iv = abiv~aboiv ;

@ TSLS ESTIMATOR @

    tscov = invpd( (lagy'iv)*invpd(iv'iv)*(iv'lagy) ) ;
    tsb   = tscov*(lagy'iv)*invpd(iv'iv)*(iv'depy)    ;
    qres  = qdepy - qlagy*tsb                         ;
    tssige2 = (qres'qres/(nn*tt-tt-1)) ;
    tscov   = tssige2*tscov            ;

res  = depy - lagy*tsb  ;

opivab = zeros(1,cols(iv)) ;
j = 0 ; do while j <= nn-1 ;
    zz = res[tt*j+1:tt*j+tt,.]'iv[tt*j+1:tt*j+tt,.] ;
    opivab = opivab|zz ;
j = j + 1 ; endo ; opivab = opivab[2:nn+1,.] ;

optwab = invpd(opivab'opivab)    ;

abov = invpd( (lagy'iv)*(optwab)*(iv'lagy) ) ;
abob = abov*(lagy'iv)*(optwab)*(iv'depy)    ;

aboe = depy - lagy*abob ;
aboj = (aboe'iv)*optwab*(iv'aboe);
abodf = cols(iv) - rows(abob) ;

/*********************
** AHN-SCHMIDT GMM1 **
*********************/

/*
**  First step GMM
*/

@ define some global variables for gmm1 @

    well = ll(tt,tt-1)                ;
    well = well[1:tt,1:tt-2]          ;
    m11  = abiv'depy ;
    m12  = abiv'lagy ;
    m21  = (yy[.,tt+1].*.well)'depy ;
    m22  = (yy[.,tt].*.well)'depy + (yy[.,tt+1].*.well)'lagy  ;
    m23  = (yy[.,tt].*.well)'lagy ;

library optmum;
#include optmum.ext;
optset ;

@ Finding initial value @

@ step 1 for gmm1 @

    m11  = abiv'depy ;
    m12  = abiv'lagy ;
    m21  = (yy[.,tt+1].*.well)'depy ;
    m22  = (yy[.,tt].*.well)'depy + (yy[.,tt+1].*.well)'lagy  ;
    m23  = (yy[.,tt].*.well)'lagy ;
    wei1 = invpd((abiv'abiv)/nn)~zeros(cols(abiv),tt-2) ;
    wei2 = zeros(tt-2,cols(abiv))~eye(tt-2) ;
    wei  = wei1|wei2 ;


                          /* GRADIENT */

proc grad1(b) ;
local gra, m  ;

    m    = (m11-m12*b)|(m21-m22*b+m23*b^2) ;
    gra  = m12|(m22-2*m23*b) ;

retp( - 2*m'wei*gra
    ) ;
endp  ;

                                /* GMM */

proc  f1(b) ;
local m    ;

    m    = (m11-m12*b)|(m21-m22*b+m23*b^2) ;

retp(  m'wei*m
    ) ;
ENDP  ;

_opgdprc = &grad1 ;

@ starting values @

b01 = -1 ;
b02 = 1 ;

_opgtol = 1e-4 ;
_opstmth = "bfgs brent";
__output = 0 ;

{b,func,grad,retcode} = optmum(&f1,b01) ;

g = m12|(m22-2*m23*b) ;
v = invpd(g'g) ;
bb1 = func~b~v ;

{b,func,grad,retcode} = optmum(&f1,b02) ;

g = m12|(m22-2*m23*b) ;
v = invpd(g'g) ;
bb2 = func~b~v  ;

bb = bb1|bb2 ;
bb = sorthc(bb,1) ;

tsbb = bb[1,2] ;

@ define some global variables for gmm1 @

REST     = YY[.,TT+1] - YY[.,TT]*tsbb ;
IVB1     = REST .*. WELL              ;
RES      = DEPY - LAGY*tsbb           ;
AS1IV    = ABIV~IVB1           ;

opivas1 = zeros(1,cols(as1iv)) ;
j = 0 ; do while j <= nn-1 ;
    zz = res[tt*j+1:tt*j+tt,.]'as1iv[tt*j+1:tt*j+tt,.] ;
    opivas1 = opivas1|zz ;
j = j + 1 ; endo ; opivas1 = opivas1[2:nn,.] ;

opwas1   = invpd(opivas1'opivas1) ;


@ step 1 for gmm1 @

                          /* GRADIENT */

proc grad2(b) ;
local gra, m  ;

    m    = (m11-m12*b)|(m21-m22*b+m23*b^2) ;
    gra  = m12|(m22-2*m23*b) ;

retp( - 2*m'opwas1*gra
    ) ;
endp  ;

                                /* GMM */

proc  f3(b) ;
local m    ;

    m    = (m11-m12*b)|(m21-m22*b+m23*b^2) ;

retp(  m'opwas1*m
    ) ;
ENDP  ;

_opgdprc = &grad2 ;

@ starting values @

b03 = -1 ;
b04 = 1 ;
_opgtol = 1e-4 ;
_opstmth = "bfgs brent";
__output = 0 ;
{b,func,grad,retcode} = optmum(&f3,b03) ;

g = m12|(m22-2*m23*b) ;
v = invpd(g'opwas1*g) ;
bb1 = func~b~v        ;

{b,func,grad,retcode} = optmum(&f3,b04) ;

g = m12|(m22-2*m23*b) ;
v = invpd(g'opwas1*g) ;
bb2 = func~b~v  ;

bb = bb1|bb2 ;
bb = sorthc(bb,1) ;

as1j = bb[1,1] ;
as1b = bb[1,2] ;
as1v = bb[1,3] ;
as1df = rows(g) - rows(as1b) ;


"" ;
"AB GMM RESULTS:" ;
"  estimates    st. err.   t-stat" ;
abb~sqrt(abv)~(abb/sqrt(abv));
"";
" J-test statistic, df, pval:"  abj abdf cdfchic(abj,abdf) ;
"";
"ABo GMM RESULTS:" ;
"  estimates    st. err.   t-stat";
abob~sqrt(abov)~(abob/sqrt(abov));
"";
" J-test statistic, df, pval:"  aboj abodf cdfchic(aboj,abodf) ;
"";
"AS1 GMM RESULTS:" ;
"  estimates    st. err.   t-stat";
as1b~sqrt(as1v)~(as1b/sqrt(as1v));
"";
" J-test statistic, df, pval:"  as1j as1df cdfchic(as1j,as1df) ;
"";


OUTPUT OFF