*
* @BaiNg( options ) x
*
* Estimates the required number of factors in a linear factor model
* using the formulas in Bai and Ng, "Determining the Number of Factors
* in Approximate Factor Models", Econometrica, 2002, vol 70, pp 191-222.
*
* Parameters:
*
*  x = nper x nvar matrix of data. The factors are nper x 1 vectors
*
* Options:
*  [CENTER]/NOCENTER - subtract means from data columns
*  STANDARDIZE/[NOSTANDARDIZE] - standardize data columns to mean zero, unit variance
*  MAX=maximum number of factors to consider
*  TITLE=title for report ["Bai-Ng Factor Determination"]
*
* Revision Schedule:
*  04/2007 Written by Tom Doan, Estima
*  09/2010 Convert report to local. Add title.
*  08/2017 Restrict the maximum number of factors to be less than an exact fit.
*
procedure BaiNg x
type rect x
*
option switch  center      1
option switch  standardize 0
option integer max			5
option string  title
*
local rect    colavg wx eigvect ex
local integer nper nvar nf
local vect    lx sigsq eigval y
local symm    xx mxx
local real    penalty1 penalty2 sigsqhat
local rect    stats
local report  freport
local string  ltitle
local integer lmax
*
if .not.%defined(x) {
   disp "###@BaiNg( options)  >>>X matrix of data<<<"
   return
}
compute nper=%rows(x),nvar=%cols(x)
*
dim wx(nper,nvar)
if center.or.standardize {
   compute colavg=%fill(1,nper,1.0/nper)*x
   ewise wx(i,j)=x(i,j)-colavg(1,j)
}
else
   compute wx=x

if standardize {
   dim lx(nvar)
   ewise lx(j)=sqrt(%normsqr(%xcol(wx,j))/nper)
   ewise wx(i,j)=wx(i,j)/lx(j)
}
*
* Do the eigenvalue calculations in the smaller dimension
*
if nper<nvar
   compute xx=%outerxx(wx)
else
   compute xx=%innerxx(wx)

eigen xx eigval eigvect
if nper>nvar
   compute eigvect=wx*eigvect
*
dim sigsq(nvar)
*
* Restrict the maximum number to 1 less than the possible number of
* factors.
*
compute lmax=%imin(max,%rows(xx)-1)
dim stats(lmax,4)
*
* Compute the two penalty factors
*
compute penalty1=((nvar+nper)/float(nvar*nper))*log((nvar*nper)/float(nvar+nper))
compute penalty2=((nvar+nper)/float(nvar*nper))*log(%min(nvar,nper))
*
* Do the calculations in reverse order so that we have an estimate of
* sigsq available for the PC variants.
*
do nf=lmax,1,-1
   compute ex=%xsubmat(eigvect,1,nper,1,nf)
   compute mxx=ex*inv(tr(ex)*ex)*tr(ex)
   *
   do j=1,nvar
      compute y=%xcol(wx,j)
      compute sigsq(j)=(%normsqr(y)-%qform(mxx,y))/nper
   end do j
   if nf==lmax
      compute sigsqhat=%avg(sigsq)
   compute stats(nf,1)=%avg(sigsq)+nf*sigsqhat*penalty1
   compute stats(nf,2)=%avg(sigsq)+nf*sigsqhat*penalty2
   compute stats(nf,3)=log(%avg(sigsq))+nf*penalty1
   compute stats(nf,4)=log(%avg(sigsq))+nf*penalty2
end do nf
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Bai-Ng Factor Determination"

report(use=freport,action=define,title=ltitle)
report(use=freport,atrow=1,atcol=1,span) ltitle
report(use=freport,atrow=2,atcol=1,align=center) "Factors" "PCP1" "PCP2" "ICP1" "ICP2"
report(use=freport,atrow=3,atcol=1,fillby=cols) %seq(1,lmax)
report(use=freport,atrow=3,atcol=2) stats
report(use=freport,action=format,atcol=2,tocol=2,tag=min,special=onestar)
report(use=freport,action=format,atcol=3,tocol=3,tag=min,special=onestar)
report(use=freport,action=format,atcol=4,tocol=4,tag=min,special=onestar)
report(use=freport,action=format,atcol=5,tocol=5,tag=min,special=onestar)
report(use=freport,action=format,atrow=3,atcol=2,picture="*.#####",align=decimal)
report(use=freport,action=show)
end
