*
* ForcedFactor computes a factorization of an input NxN covariance matrix SIGMA
* which controls either
*
*   (a) a column or block of columns in the factorization itself
*   (b) a row or block of rows in the inverse of the factorization
*
* In the case (a), you provide an RxN matrix A. An RxR matrix PI is computed so
* that PI is upper triangular and A x PI makes up the first R columns in a matrix
* F that factors SIGMA (that is FF'=SIGMA). In particular, if R=1, then the first
* column of F is a scale multiple of A. Because of the structure of PI, if R>1,
* the first column of F is a scale multiple of the first column in A, the second
* column in F will be a linear combination of the first two columns of A, etc.
*
* In the case (b), you provide an NxR matrix A. An RxR matrix PI is computed so
* that PI is lower triangular and PI x A makes up the first R rows of the inverse
* of a matrix F that factors SIGMA. In particular, if R=1, the first row in F**-1
* is a scale multiple of A.
*
* This allows you to either force an orthogonalized component to hit the variables
* in a specific pattern (done by setting a column of the factorization), or to
* force that an orthogonalized component be formed from a particular linear
* combination of innovations (forces a row in the inverse). You set the one
* component, and ForcedFactor builds a factorization around it. With multiple
* input components, you can force the factorization to include linear
* combinations of your inputs.
*
* Syntax:
*
* @ForcedFactor( options ) sigma a f
*
*   Parameters:
*      sigma   SYMMETRIC covariance matrix (input)
*      a       RECTANGULAR (or VECTOR) - one dimension should be the same as sigma (input)
*      f       RECTANGULAR - output factor of sigma
*
*   Options:
*      force=[row]/column indicates whether a scale of "a" is to be a row in the
*        inverse or a column in the factorization itself.
*
* Examples:
*
* In a four variable system where the first two variables are two interest rates,
*
* @ForcedFactor sigma ||1.0,-1.0,0.0,0.0|| f1
* @ForcedFactor(force=column) sigma ||1.0,1.0,0.0,0.0|| f2
*
* f1 will be a factorization where the first orthogonal component is the
* (non-orthogonal) innovation in the difference between the rates. f2 will be a
* factorization where the first orthogonal component loads equally onto the two
* interest rates, and hits none of the other variables contemporaneously.
*
*  Revision Schedule
*   09/02 Written by Estima
*   11/05 Revised to allow multiple rows or columns
*
* %SVDZeros picks out the "r" smallest elements in the vector v, returning a
* VECT[INT] with the subscripts. This is used to locate the zero singular values
* in a singular value decomposition.
*
function %SVDZeros v r
type vect[int] %SVDZeros
type vector v
type integer r
local vector ranks
local integer i j
*
compute ranks=%ranks(v)
dim %SVDZeros(r)
compute j=0
do i=1,%rows(v)
   if ranks(i)<=r {
      compute j=j+1
      compute %SVDZeros(j)=i
   }
end do i
end
*
procedure ForcedFactor sigma a f
type symm sigma
type rect a
type rect *f
*
option choice force 1 row column
*
local integer n r
local rect vv w pi p ff fill
local symm pp
local vect[rect] ss
local vect[int]  zeros
*
compute n=%rows(sigma)
compute r=%imin(%rows(a),%cols(a))
if %imax(%rows(a),%cols(a))<>n {
   disp "@ForcedFactor sigma a f"
   disp "A is dimension" %rows(a) "x" %cols(a)
   disp "SIGMA " n "x" n
   disp "One dimension of A must match SIGMA"
   return
}
if force==2 {
   compute p=inv(%decomp(sigma))
   if %rows(a)==n
      compute vv=p*a
   else
      compute vv=p*tr(a)
   compute w=%blockglue(||vv,%zeros(n,n-r)||)
   compute ss=%svdecomp(w)
   compute zeros=%SVDZeros(ss(2),n-r)
   compute pp=inv(tr(vv)*vv)
   compute pi=%psdfactor(pp,%seq(r,1))
   dim fill(n,n-r)
   ewise fill(i,j)=ss(1)(i,zeros(j))
   compute f=inv(p)*%blockglue(||vv*pi,fill||)
}
else {
   compute p=%decomp(sigma)
   if %cols(a)==n
      compute vv=a*p
   else
      compute vv=tr(a)*p
   compute w=%blockglue(||vv|%zeros(n-r,n)||)
   compute ss=%svdecomp(w)
   compute zeros=%SVDZeros(ss(2),n-r)
   compute pi=inv(%decomp(%outerxx(vv)))
   dim fill(n-r,n)
   ewise fill(i,j)=ss(3)(j,zeros(i))
   compute ff=%blockglue(||pi*vv|fill||)
   compute f=p*inv(ff)
}
end



