********************************************************************************
PROCEDURE FACTOR invec outvec

********************************************************************************
*
* Factors  1 - a(1)*L - a(2)*L^2 - ... - a(p)*L^p
*
*  into    (1-b(1)*L)*(1-b(2)*L)*...*(1-b(p)*L)
*
*  where  [a(1),...,a(p)]  =  invec
*  and    outvec  =  [b(1),...,b(p)]
*
* [PRINT]/NOPRINT  Print output or not
* BDROP = [0]      Number of entires to drop at the beginning of the
*                     input vector
* EDROP = [0]      Number of entires to drop at the end of the
*                     input vector
*
* The bdrop and edrop command allows one to use the %beta vector, which is
*     output from the LINREG command, as the input vector when there
*     are other variables included with lagged dependent variables.
*
* E.g., for
*   LINREG y
*   #constant y{1 to 4} dum1974 dum1982
* one would use
*   @FACTOR(bdrop=1,edrop=2) %beta roots
*
* By Norman Morin (m1njm00@frb.gov), 7/97
*
*
********************************************************************************
OPTION SWITCH  PRINT 1
OPTION INTEGER BDROP 0
OPTION INTEGER EDROP 0

TYPE VECTOR  invec *outvec

LOCAL INTEGER p
LOCAL VECT newvec
LOCAL RECT F

DIM newvec(%ROWS(invec)-edrop-bdrop)
EWISE newvec(i) = invec(i+bdrop)
COMPUTE p = %ROWS(newvec)

DIM F(p,p)
COMPUTE F = %CONST(0.0)

DO i = 1,p
  DO j = 1,p
    IF j.eq.(i-1)
      COMPUTE F(i,j) = 1.
  END DO
  COMPUTE F(1,i) = newvec(i)
  IF i.ne.1
    COMPUTE F(i,p) = 0.
END DO

EIGEN(cvalues=cvec) F outvec

IF PRINT.eq.1 {
  DISP
  DISP 'Factoring the polynomial:'
  DISP ######## 1
  DO jj = 1,p
     DISP '   -' ###.##### newvec(jj)  "* L^" @-1  jj
  END
  DISP
  DISP '  The roots of the polynomial in the lag operator: '
  DO jj = 1,p
    IF abs(%imag(cvec(jj))).lt.(.00000001)
       DISP ########.##### %real(cvec(jj))
    ELSE IF abs(%real(cvec(jj))).lt.(.00000001)
       DISP %imag(cvec(jj)) @-1 "*i"
    ELSE
       DISP ########.##### %real(cvec(jj)) @-1 '  + ' @-1 $
            ##.##### %imag(cvec(jj)) @-1 "*i" @33  '('  $
            @-1 #.##### (%real(cvec(jj))**2+%imag(cvec(jj))**2)**.5 @-1 ')'
  END DO
  DISP
  DISP '    implying the roots of the difference equation:'
  DO jj = 1,p
    IF abs(%imag(1./cvec(jj))).lt.(.00000001)
       DISP ########.##### %real(1./cvec(jj))
    ELSE IF abs(%real(1./cvec(jj))).lt.(.00000001)
       DISP %imag(1./cvec(jj)) @-1 "*i"
    ELSE
       DISP ########.##### %real(1./cvec(jj)) @-1 '  + ' @-1 $
            ##.##### %imag(1./cvec(jj)) @-1 "*i"  @33  '('  $
            @-1 #.##### (%real(1./cvec(jj))**2+%imag(1./cvec(jj))**2)**.5 @-1
')'
  END DO
  DISP

}

END
