*
* 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