/*GAUSS CODE FOR GMM ESTIMATION OF THE NUMBER OF LATENT FACTORS
BY MARCOS FABRICIO PEREZ  OCTUBER 20 2008
*/


new;
cls;
seed=1;
rndseed seed;
fabtime1=date;


/*INITIAL DATA*/

bandw = 0 ;      @ if bandw = 0, White-Correction, if auto =0 set bandw =0 @
defw=0;          @if defw=1 andrews, if  defw = 0 newey west@
numero2=100;     @number of different partitions used in estimation@
changeornot=1;   @1 if we do change partitions@


nameout="put_a_output_name_here" ;  



/* LOAD DATA HERE AND CALL IT "companybase"
THIS IS AN t BY n MATRIX OF THE VARAIBLES TO BE ANALIZED
*/


 companybase=rndns(300,10,seed); // for example with random numbers we should have 0 factors//





//DO NOT CHANGE FROM HERE

/*START ITERATIONS FOR CHANGE OF BASES*/

base=cols(companybase);
tt=rows(companybase);

company=companybase;
numinstr=floor(base/2);

iterchan=1;
do while iterchan<=numero2;
if changeornot==1; 
alfa=seqa(1,1,base);
phil=rndi(base,1); 
thao=alfa~phil;
lista = sortc(thao,2);
lista1=lista[.,1]';
if iterchan==1; basmat=lista1; 
else; basmat=basmat|lista1; endif;
company=company[.,lista1];
endif;

//  TIMER 
cls;
format /rd 12,4;

locate 4,1;
locate 4,30;
"iterchan" iterchan;
fabtime2=ethsec(fabtime1,date);
locate 6,1;
"time running   "     etstr(fabtime2);
if iterchan>5;
fabtime3=(fabtime2/(iterchan-1))*numero2;
locate 8,1; "time left     "     etstr(fabtime3-fabtime2);
locate 10,1; "total time    "     etstr(fabtime3);
endif;
 
/*BAI'S METHOD*/


proc factorn(x,n);

local a,b,c,d,factor,loadings;

{a,b}=eigv(x'x);
c=a~b';
d=rev(sortc(c,1));
loadings=sqrt(n)*d[.,2:n+1]';
factor=x*loadings/n;

retp(factor);

endp;

proc factort(x,t);

local a,b,c,d,factor;

{a,b}=eigv(x*x');
c=a~b';
d=rev(sortc(c,1));
factor=sqrt(t)*d[.,2:t+1]';

retp(factor);

endp;

/* load data FOR BAI'S METHOD*/

data=company;

kmax=8;    @maximum amount of factors to be tested@

/* # of Observations*/

T=rows(data);
N=rows(data');

/* Factor calculation*/

if N<T;
factores=factorn(data,N); else;
factores=factort(data,T);
endif;

/* Rescaled Factors */

rfact=zeros(rows(factores),kmax);

for i(1,kmax,1);
rfact[.,i]=factores[.,i]*sqrt(factores[.,i]'factores[.,i]/T);
endfor;

/* Squared errors */

V=zeros(kmax,1);
err=zeros(n,1);
for j(1,kmax,1);
ff=rfact[.,1:j];
for i(1,n,1);
lambda=invpd(ff'ff)*ff'data[.,i];
err[i,1]=(data[.,i]-ff*lambda)'(data[.,i]-ff*lambda)/(n*t);
endfor;
V[j,1]=sumc(err);
endfor;

/* PC criteria for choosing the number of factors */

pc=zeros(kmax,3);
cn=sqrt(n)|sqrt(t);
cn=minc(cn);

for i(1,kmax,1);
pc[i,1]=v[i]+i*v[kmax]*((n+t)/(n*t))*ln((n*t)/(n+t));
endfor;

for i(1,kmax,1);
pc[i,2]=v[i]+i*v[kmax]*((n+t)/(n*t))*ln(cn^2);
endfor;

for i(1,kmax,1);
pc[i,3]=v[i]+i*v[kmax]*ln(cn^2)/cn^2;
endfor;

/* IC criteria for choosing the number of factors */

ic=zeros(kmax,3);

for i(1,kmax,1);
ic[i,1]=ln(v[i])+i*((n+t)/(n*t))*ln((n*t)/(n+t));
endfor;

for i(1,kmax,1);
ic[i,2]=ln(v[i])+i*((n+t)/(n*t))*ln(cn^2);
endfor;

for i(1,kmax,1);
ic[i,3]=ln(v[i])+i*ln(cn^2)/cn^2;
endfor;

criic=minindc(ic+100000);
cripc=minindc(pc+100000);

baiIC1=criic[1]; 
baiIC2=criic[2];
baiIC3=criic[3];
baiPC1=cripc[1];
baiPC2=cripc[2];
baiPC3=cripc[3];

/*AHN AND PEREZ GMM METHOD  */


/*creation fo groups*/

company1=company[.,1:base-numinstr];
company3=ones(tt,1)~company[.,base-numinstr+1:base];  /* intrument group, change here if you have observable instrumnents*/


nn=cols(company1);
inst=cols(company3);


/*estimation*/

p=0;   
do while p<inst-1 and p<nn;

@creating yy@

ystar=company1[.,1:nn-p];
yy=vec(ystar);


@creating zz regresors@

ident=eye(nn-p);
if P==0; RR=ONES(TT,1); else;
rr=ones(tt,1)~company1[.,nn-p+1:nn]; endif;
 
zz=ident.*.rr;
xx=ident .*. company3 ;
@creating intruments@

x=company3;
z= rr ;
xz = x'z;
ixx=invpd(x'x);


@first 2SLS@

/*delta= (invpd( (z'x)*invpd(x'x)*(x'z) )*(z'x)*inv(x'x)*x')*ystar;*/


delta= invpd( xz'*ixx*xz )*(xz'*ixx)*x'*ystar;
e=ystar-(z*delta);
ee=vec(e);



/*second round 3SIV*/


gg=zeros(tt,(nn-p)*inst);

jj=0;
do while jj<=nn-p-1;

gg[.,1+(jj*inst):(jj*inst)+inst]=company3.*ee[1+(tt*jj):(tt*jj)+tt];

jj=jj+1;
endo;


 bb=ones(1+(nn-p)*p,1);

    
/*   Procedure for Newey-West   */

proc nw(uv,l) ;
local j,omeo,n,w,ome  ;

    n    = rows(uv)  ;
    omeo = uv'uv/n   ;

    j = 1; do while j <= l  ;

        w    = 1 - j/(l+1)       ;
        ome  = uv[j+1:n,.]'uv[1:n-j,.]/n ;
        omeo = omeo + w*(ome+ome')       ;

    j = j + 1 ; endo ;

retp(omeo) ;
endp       ;

/*   Procedure for Andrews (1991) */

proc quadk(z) ;
local rr ;

  rr = ( 3/(6*pi*z/5)^2 )
        *( 5*sin(6*pi*z/5)/(6*pi*z) - cos(6*pi*z/5) ) ;

retp(rr) ;
endp     ;

proc andrews(uv,l) ;
local j,omeo,n,w,ome  ;

    n    = rows(uv)  ;
    omeo = uv'uv/n   ;

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

        w    = j/(l+1)                    ;
        ome  = uv[j+1:n,.]'uv[1:n-j,.]/n  ;
        omeo = omeo + quadk(w)*(ome+ome') ;

    j = j + 1 ; endo ;

retp((n/(n-rows(bb)))*omeo) ;
endp       ;



if defw ==1; w= inv(tt*andrews(gg,bandw));
else; w= invPD(tt*nw(gg,bandw)); endif;

delta1=invPD(zz'xx*w*xx'zz)*(zz'xx*w*xx'yy);


/*Jtest*/


jtest1=(yy-zz*delta1)'xx*w*xx'(yy-zz*delta1);


momencond=(nn-p)*inst;
parame=(nn-p)*(p+1);

df=momencond-parame; 

JTOD=JTEST1;

if p==0; dfsto=df; paramesto=parame; jtodas=jtod;
else; dfsto=dfsto|df; paramesto=paramesto|parame; jtodas=jtodas|jtod;
endif;



p=p+1;
endo;


/*save the results for the change of bases and the number of factros*/
/* PENALTIES OF MSC*/

penaltybic=jtodas/ln(tt)-dfsto;
penaltybic2=jtodas/log(tt)-dfsto;
penaltyaic=jtodas - 2*dfsto;
penaltyhqq=jtodas/ln(ln(tt))-DFSTO;
penaltyfabri=jtodas/ln(tt)-(dfsto+paramesto);
auxilio=rows(jtodas);


pbic=minindc(penaltybic)-1;
pbic2=minindc(penaltybic2)-1;
paic=minindc(penaltyaic)-1;
phqq=minindc(penaltyhqq)-1;
pfabri=minindc(penaltyfabri)-1;

pvaluet=cdfchic(jtodas,dfsto);
SHT_a=seqa(0,1,P)~pvaluet;
SHT1=SELIF(SHT_a, SHT_a[.,2] .gt 0.05); 
SHT2=MINC(SHT1[.,1]);

if ITERchan==1; 
FAC_SHT=sht2;
pbic_t=pbic;
pbic2_t=pbic2;
paic_t=paic;
phqq_t=phqq;
pfabri_t=pfabri;
baiIC1_t=baiIC1; 
baiIC2_t=baiIC2;
baiIC3_t=baiIC3;
baiPC1_t=baiPC1;
baiPC2_t=baiPC2;
baiPC3_t=baiPC3;
jtsto=jtodas;


else; 
FAC_SHT=FAC_SHT|sht2;
pbic_t=pbic_t|pbic;
pbic2_t=pbic2_t|pbic2;
paic_t=paic_t|paic;
phqq_t=phqq_t|phqq;
pfabri_t=pfabri_t|pfabri;
baiIC1_t=baiIC1_t|baiIC1; 
baiIC2_t=baiIC2_t|baiIC2;
baiIC3_t=baiIC3_t|baiIC3;
baiPC1_t=baiPC1_t|baiPC1;
baiPC2_t=baiPC2_t|baiPC2;
baiPC3_t=baiPC3_t|baiPC3;
jtsto=jtodas~jtsto;


endif;

iterchan=iterchan+1;
endo;
cls;

/*SAVE DATA FOR CHANGE OF BASES*/

vv=seqa(0,1,P);
vvbai=seqa(1,1,kmax);
fBIC=100*counts(pBIC_T,vv)./numero2; 
fBIC2=100*counts(pBIC2_T,vv)./numero2; 
fAIC=100*counts(pAIC_T,vv)./numero2; 
fHQQ=100*counts(pHQQ_T,vv)./numero2; 
fFABRI=100*counts(pFABRI_T,vv)./numero2; 
fSHT=100*counts(FAC_SHT,vv)./numero2; 
Fbai1=100*counts(baiIC1_t,vvbai)./numero2; 
Fbai2=100*counts(baiIC2_t,vvbai)./numero2;
Fbai3=100*counts(baiIC3_t,vvbai)./numero2;
Fbai4=100*counts(baiPC1_t,vvbai)./numero2;
Fbai5=100*counts(baiPC2_t,vvbai)./numero2;
Fbai6=100*counts(baiPC3_t,vvbai)./numero2;


porcentaje1=fSHT~fbic~fbic2~ffabri~Faic~fhqq;
porcentaje2=Fbai1~fbai2~fbai3~fbai4~fbai5~fbai6;



/*results*/

cls;
outwidth 120;
nombre=nameout$+"_march6.index";
output file =^nombre reset; 


"
tt" tt;
"base" base;
"bandw" bandw             "if bandw = 0, White-Correction, if auto =0 set bandw =0 ";
"defw" defw               "if defw=1 andrews, if  defw = 0 newey west";
"numero" numero2           "iterations change bases";
"numinstr" inst;

"NUMBER OF FACTORS AHN PEREZ";
"                   SHT         BIC          BIC2           BIC3        AIC            HQQ";

averages=vv'*porcentaje1;
aaaa=1/0;

vv~porcentaje1;
aaaa~averages;

"NUMBER OF FACTORS BAI";
"                   IC1         IC2          IC3           PC1        PC2            PC3";

averagesbai=vvbai'*porcentaje2/100;
aaaa=1/0;

vvbai~porcentaje2/100;
aaaa~averagesbai;
OUTPUT OFF;