*
* 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 NxR 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 RxN 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/2002 Written by Estima
* 11/2005 Revised to allow multiple rows or columns
* 10/2008 Revised to require all three parameters
* 10/2009 Add %abs(...) to %SVDZeros to deal with the rare occurrence of a
* negative singular value.
* 11/2009 Correct comments
*
* %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(%abs(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
*
if %defined(sigma)==0.or.%defined(a)==0.or.%defined(f)==0 {
disp "The syntax for this procedure is:"
disp " "
disp "@ForcedFactor sigma a f"
disp " "
disp "sigma SYMMETRIC covariance matrix (input)"
disp "a RECTANGULAR (or VECTOR) - one dimension should be the same as sigma (input)"
disp "f RECTANGULAR - output factor of sigma"
disp " "
disp "All three parameters are required."
disp ""
return
}
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