/*
**  CREATING LT MATRIX
*/

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 AM Variables
*/

proc ammat(x,n,t)  ;
local k, amx, i, j ;

    k = cols(x) ;
    amx  = zeros(n*t,k*t)   ;
    i = 1 ; do while i <= n ;

        j = 1 ; do while j <= t ;

            amx[(i-1)*t+1:i*t,(j-1)*k+1:j*k]
                      = ones(t,1) .*. x[(i-1)*t+j,.] ;

        j = j + 1 ; endo ;

    i = i + 1 ; endo ;

retp(amx) ;
endp      ;


/*
** Creating PV Variables
*/

proc pvmat1(x,n,t)  ;
local k, pvx, i    ;

    k = cols(x) ;
    pvx  = zeros(n*t,k)     ;
    i = 1 ; do while i <= n ;

     pvx[(i-1)*t+1:i*t,.]
        = ones(t,1) .*. meanc(x[(i-1)*t+1:i*t,.])';

    i = i + 1 ; endo ;

retp(pvx) ;
endp      ;

proc pvmat(x,n,t)  ;
local k, pvx, i    ;

    k = cols(x) ;
    pvx  = zeros(n,k)     ;
    i = 1 ; do while i <= n ;

     pvx[i,.]
        = meanc(x[(i-1)*t+1:i*t,.])';

    i = i + 1 ; endo ;

retp(pvx .*. ones(t,1) ) ;
endp      ;


/*
** Creating QV Variables
*/

proc qvmat(x,n,t) ;

retp(x - pvmat1(x,n,t)) ;
endp                    ;

/*
**  Creating  BMS instruments
*/

proc bmsmat(x,n,t) ;
local bmsx,xx,k,i,j   ;

    k    = cols(x)        ;
    xx   = qvmat(x,n,t)   ;
    bmsx = zeros(n*t,k*(t-1)) ;
    i = 1 ; do while i <= n   ;

        j = 1 ; do while j <= t-1 ;

            bmsx[(i-1)*t+1:i*t,(j-1)*k+1:j*k]
                      = ones(t,1) .*. xx[(i-1)*t+j,.] ;

        j = j + 1 ; endo ;

    i = i + 1 ; endo ;

retp(bmsx) ;
endp       ;

/*
**  Creating Arellano-Bover Instruments
*/

proc abmat(x,z,n,t) ;
local abins,i       ;

    abins = zeros(n,t*cols(x)+cols(z)) ;
    i  = 1 ; do while i <= n    ;

       abins[i,.]
         = z[(i-1)*t+1,.]~vec(x[(i-1)*t+1:i*t,.])' ;

    i = i + 1 ; endo ;

retp( abins .*. ll(t,t-1)  ) ;
endp          ;

/*
**  PROCEDURE FOR 2SLS
**  INPUT : REG = REGRESSORS;
**  IV = INSTRUMENTS; Y = DEPENDENT VAR.
**  OUTPUT: ESTIMATOR
*/

proc twsls(reg,iv,y) ;
local a,b,c ;

    a = iv'reg       ;
    b = invpd(iv'iv) ;
    c = iv'y         ;

retp( invpd(a'b*a)*a'b*c ) ;
endp ;

/*
**  2SLS WITH COVARIANCE MATRIX AND HANSEN TEST
**  INPUT : REG = REGRESSORS;
**          IV = INSTRUMENTS; Y = DEPENDENT VAR.
*/

proc (3) =  twslsc(reg,iv,y,n,t) ;
local a,b,c,d,e,ssq,han ;

    a = qvmat(iv,n,t)'qvmat(reg,n,t)         ;
    b = invpd( qvmat(iv,n,t)'qvmat(iv,n,t) ) ;
    c = qvmat(iv,n,t)'qvmat(y,n,t)           ;
    d = invpd(a'b*a)*a'b*c ;
    e = y - reg*d ;
    ssq = e'e/(n*t-n-cols(reg)) ;
    han = (e'iv)*invpd(iv'iv)*(iv'e) ;

retp( d, ssq*invpd(a'b*a), han ) ;
endp ;

/*
**  3SLS ESTIMATION
**  INPUT : REG,IV,Y
**  OUTPUT: ESTIMATES, COV. MATRIX OF ESTIMATES
*/

proc (2) =  thsls(reg,iv,y,n,t) ;
local tb,te,a,b1,b2,b,c,d,j,sig ;

    tb = twsls(reg,iv,y) ;
    te = y - reg*tb      ;

    sig  = zeros(t,t) ;
    j = 1 ; do while j <= n ;
        sig = sig + te[(j-1)*t+1:(j-1)*t+t]
                      *te[(j-1)*t+1:(j-1)*t+t]' ;
    j = j + 1 ; endo ;

    sig  = sig/n ;

    b = zeros(n*t,cols(iv));

    j = 1 ; do while j <= n ;

        b[(j-1)*t+1:j*t,.] = sig*iv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    b  = iv'b ;
    b  = invpd(b) ;
    a = iv'reg ;
    c = iv'y   ;
    d = invpd(a'b*a) ;

retp( d*a'b*c, d ) ;
endp ;

/*
**  GMM ESTIMATION
**  INPUT : REG,IV,Y
**  OUTPUT: ESTIMATES, COV. MAT. OF ESTIMATES, Hansen Test
*/

proc (3) = gmm(reg,iv,y,n,t) ;
local tb,te,a,b,c,d,f,j,test,ive  ;

    tb = twsls(reg,iv,y) ;
    te = y - reg*tb      ;

    a = iv'reg ;
    c = iv'y   ;

    b = zeros(n*t,cols(iv)) ;
    j = 1 ; do while j <= n ;

        b[(j-1)*t+1:j*t,.]
            = te[(j-1)*t+1:j*t]*te[(j-1)*t+1:j*t]'
                *iv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    b = invpd(iv'b) ;
    d = invpd( a'b*a ) ;
    f = d*a'b*c ;

    ive = iv'(y-reg*f) ;
    test = ive'b*ive   ;

retp( d*a'b*c, d, test ) ;
endp ;

/*
**  PROCEDURE FOR MGIV2
**  INPUT : REG,IV1,IV2,Y,SIG,N,T
**  OUTPUT: ESTIMATES, COV. MAT. OF ESTIMATES
*/

proc (3) = mgiv2(reg,iv1,iv2,y,n,t) ;
local tb,te,sig,a,b,c,d,tiv,j,is,m,et,mtiv,f,e,bb  ;

    tb = twsls(reg,iv1~iv2,y) ;
    te = y - reg*tb      ;

    sig  = zeros(t,t) ;
    j = 1 ; do while j <= n ;
        sig = sig + te[(j-1)*t+1:(j-1)*t+t]
                     *te[(j-1)*t+1:(j-1)*t+t]' ;
    j = j + 1 ; endo ;

    sig  = sig/n ;

    is   = invpd(sig) ;
    et   = ones(t,1)  ;
    m    = is - is*et*invpd(et'is*et)*et'is ;

    tiv  = zeros(n*t,cols(iv1)+cols(iv2));
    mtiv = zeros(n*t,cols(iv1)+cols(iv2));
    j = 1 ; do while j <= n ;

        tiv[(j-1)*t+1:j*t,.]
        = (m*iv1[(j-1)*t+1:j*t,.])~(is*iv2[(j-1)*t+1:j*t,.]) ;
        mtiv[(j-1)*t+1:j*t,.] = sig*tiv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    a    = tiv'reg ;
    b    = invpd(tiv'mtiv) ;
    c    = tiv'y   ;
    d = invpd(a'b*a) ;
    f = d*a'b*c ;
    e = y - reg*f ;

    bb = zeros(n*t,cols(tiv)) ;
    j  = 1 ; do while j <= n ;

        bb[(j-1)*t+1:j*t,.]
            = e[(j-1)*t+1:j*t]*e[(j-1)*t+1:j*t]'
                *tiv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    bb = tiv'bb ;

retp( f, d, d*a'b*bb*b*a*d  ) ;
endp ;

/*
**  PROCEDURE FOR MGIV1
**  INPUT : REG,IV1,IV2,Y,SIG,N,T
**  OUTPUT: ESTIMATES, COV. MAT. OF ESTIMATES
*/

proc (3) = mgiv1(reg1,reg2,iv1,iv2,y,n,t) ;
local tb,te,sig,a,b,c,d,tiv,j,is,m,et,mtiv,f,e,bb  ;

    tb = twsls(reg1~reg2,iv1~iv2,y) ;
    te = y - (reg1~reg2)*tb         ;

    sig  = zeros(t,t)       ;
    j = 1 ; do while j <= n ;
        sig = sig + te[(j-1)*t+1:(j-1)*t+t]
                        *te[(j-1)*t+1:(j-1)*t+t]' ;
    j = j + 1 ; endo ;

    sig  = sig/n ;

    is   = invpd(sig) ;
    et   = ones(t,1)  ;
    m    = is - is*et*invpd(et'is*et)*et'is ;

    tiv  = zeros(n*t,cols(iv1));
    mtiv = zeros(n*t,cols(iv1));
    j = 1 ; do while j <= n ;

        tiv[(j-1)*t+1:j*t,.]
                = (m*iv1[(j-1)*t+1:j*t,.]) ;
        mtiv[(j-1)*t+1:j*t,.] = sig*tiv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    a    = tiv'reg1 ;
    b    = invpd(tiv'mtiv) ;
    c    = tiv'y   ;
    d = invpd(a'b*a) ;
    f = d*a'b*c ;
    e = y - reg1*f ;

    bb = zeros(n*t,cols(tiv)) ;
    j  = 1 ; do while j <= n ;

        bb[(j-1)*t+1:j*t,.]
            = e[(j-1)*t+1:j*t]*e[(j-1)*t+1:j*t]'
                *tiv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    bb = tiv'bb ;

retp( f, d, d*a'b*bb*b*a*d  ) ;
endp ;


/*
**  RE-GLS for Time-HET/AUTO errors.
**  INPUT : REG,Y,N,T
**  OUTPUT: ESTIMATES, COV. MAT. OF ESTIMATES
*/

proc (3) = regls(x,y,n,t) ;
local olsb,olse,phi,a,b,c,d,isx,j,is,m,et,mtiv,f,e,bb,qvy,qvx ;

    olsb = invpd(x'x)*(x'y) ;
    olse = y - x*olsb       ;

    phi  = zeros(t,t)       ;
    j = 1 ; do while j <= n ;
        phi = phi
         + olse[(j-1)*t+1:(j-1)*t+t]
            *olse[(j-1)*t+1:(j-1)*t+t]' ;
    j = j + 1 ; endo ;

    phi  = phi/n ;

    is   = invpd(phi) ;

    isx  = zeros(n*t,cols(x));
    j = 1 ; do while j <= n ;

        isx[(j-1)*t+1:j*t,.]
                = is*x[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    b    = invpd(x'isx) ;
    c    = isx'y        ;
    f    = b*c ;
    e    = y - x*f ;

    bb = zeros(n*t,cols(isx)) ;
    j  = 1 ; do while j <= n ;

        bb[(j-1)*t+1:j*t,.]
            = e[(j-1)*t+1:j*t]*e[(j-1)*t+1:j*t]'
                *isx[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    bb = isx'bb ;

retp( f, b, b*bb*b  ) ;
endp ;


/*
**  Kiefer's GLS for Fixed Effects
**  INPUT : REG,Y,N,T
**  OUTPUT: ESTIMATES, COV. MAT. OF ESTIMATES
*/

proc (3) = kiefer(x,y,n,t) ;
local olsb,olse,sig,a,b,c,d,tiv,j,is,m,et,mtiv,f,e,bb,qvy,qvx ;

    qvx = qvmat(x,n,t) ;
    qvy = qvmat(y,n,t) ;
    olsb  = invpd(qvx'qvx)*(qvx'qvy) ;
    olse  = qvy - qvx*olsb           ;

    sig  = zeros(t,t)       ;
    j = 1 ; do while j <= n ;
        sig = sig
         + olse[(j-1)*t+1:(j-1)*t+t]
            *olse[(j-1)*t+1:(j-1)*t+t]' ;
    j = j + 1 ; endo ;

    sig  = sig/n ;

    is   = pinv(sig) ;

    tiv  = zeros(n*t,cols(qvx));
    j = 1 ; do while j <= n ;

        tiv[(j-1)*t+1:j*t,.]
                = is*qvx[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    b    = invpd(qvx'tiv) ;
    c    = tiv'qvy        ;
    f    = b*c ;
    e = qvy - qvx*f ;

    bb = zeros(n*t,cols(tiv)) ;
    j  = 1 ; do while j <= n ;

        bb[(j-1)*t+1:j*t,.]
            = e[(j-1)*t+1:j*t]*e[(j-1)*t+1:j*t]'
                *tiv[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    bb = tiv'bb ;

retp( f, b, b*bb*b  ) ;
endp ;

/*
**  WITHIN with HET/AUTO adjustment
**  INPUT : REG,Y,N,T
**  OUTPUT: ESTIMATES, COV. MAT.
*/

proc (2) = w_ha(x,y,n,t) ;
local olsb,olse,sig,a,b,c,d,tiv,j,is,m,et,mtiv,f,e,bb,qvy,qvx ;

    qvx = qvmat(x,n,t) ;
    qvy = qvmat(y,n,t) ;
    b  = invpd(qvx'qvx)*(qvx'qvy) ;
    e  = qvy - qvx*b           ;

    bb = zeros(n*t,cols(qvx)) ;
    j  = 1 ; do while j <= n ;

        bb[(j-1)*t+1:j*t,.]
            = e[(j-1)*t+1:j*t]*e[(j-1)*t+1:j*t]'
                *qvx[(j-1)*t+1:j*t,.] ;

    j = j + 1 ; endo ;

    bb = qvx'bb ;

retp( b, invpd(qvx'qvx)*bb*invpd(qvx'qvx) ) ;
endp ;