*
*  PRINCOMP is a simple procedure for extracting principal components
*  from a set of series.
*
*  PRINCOMP  start  end
*  #   list of variables (regression format)
*
*  Options:
*    NCOMPS = number of components [number of variables]
*    PCOMP  = VECT[SERIES] of the extracted principal components
*
*  Revision Schedule
*
*  05/01 Modified to use VECTOR of SERIES for the output,
*        rather than a block of numbered series.
*  01/06 Added CORR and CENTER options for analyzing data in correlation
*        or centered form. (Requires version 6.20).
*
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  print   0
*
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
*
equation stemp *
inquire(equation=stemp) startl<<start endl<<end
*
cmom(center=(center.or.corr),equation=stemp) startl endl
if corr {
   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)
}

eigen cmom eigenval eigenvec
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





