*
* @PRINCOMP is a simple procedure for extracting principal components
* from a set of series. Use @PRINFACTORS if you need more detailed
* information.
*
*  @PRINCOMP( options )  start  end  pcomp
*  #   list of variables (regression format) (if no EQUATION option)
*
*     start  end  : Range to use
*     pcomp       : Use this to supply a name for the VECT[SERIES]
*                   that will be created by the procedure to hold
*                   the extracted principal components
*
*  Options:
*    EQUATION = equation providing variables
*    NCOMPS = number of components [number of variables]
*    CENTER/[NOCENTER] subtract means from data
*    CORR/[NOCORR]     convert to matrix of correlations
*    STANDARDIZE/[NOSTANDARDIZE] same as CORR
*    PRINT/[NOPRINT]
*
*  Revision Schedule
*
*  05/2001 Modified to use VECTOR of SERIES for the output,
*          rather than a block of numbered series.
*  01/2006 Added CORR and CENTER options for analyzing data in correlation
*          or centered form. (Requires version 6.20).
*  08/2006 Corrections to the comments describing the
*          procedure syntax.
*  10/2008 Corrected to have the number of components default
*          to the number of variables, rather than zero.
*  01/2010 Add EQUATION option
*  08/2017 Add STANDARDIZE option
*
procedure princomp start end pcomp
type integer start end
type vect[series] *pcomp
*
option integer  ncomps   0
option switch   center   0
option switch   corr     0
option switch   standard 0
option switch   print    0
option equation equation
*
local equation stemp
local rect     eigenvec
local symm     cmom
local vect     eigenval princoef scalef
local int      n
local int      icomp nvars
local int      startl endl
*
* Get the list of variables
*
if %defined(equation)
   compute stemp=equation
else
   equation stemp *
inquire(equation=stemp) startl<<start endl<<end
*
cmom(center=(center.or.corr.or.standard),equation=stemp) startl endl
if corr.or.standard {
   compute cmom=%cvtocorr(%cmom)
   compute scalef=%sqrt(%xdiag(%cmom))
}
else {
   compute cmom=%cmom
   compute scalef=%fill(%ncmom,1,1.0)
   if .not.center
      compute %meanv=%const(0.0)
}

compute nvars = %rows(cmom)

eigen cmom eigenval eigenvec

if ncomps>nvars
{
   display "###@PRINCOMP. NCOMPS must be less than or equal to the number of variables."
   return
 }

if ncomps==0
   compute n = nvars
else
   compute n = ncomps

if %defined(pcomp)
   dim pcomp(n)
if print
   report(action=define)
do icomp=1,n
   compute princoef=%xcol(eigenvec,icomp)
   if print
      report(col=new,atrow=1,fillby=cols) eigenval(icomp) princoef
   if %defined(pcomp)
      set pcomp(icomp) startl endl = %dot((%eqnxvector(stemp,t)-%meanv)./scalef,princoef)
end do icomp
if print
  report(action=show)
end
