function phi = stationary(a,sig)
%STATIONARY maps a square matrix to one with eigenvalues inside the unit circle
%
%           Phi = stationary(A,Sig) takes the matrix A, which is a general
%           square matrix, and Sig, which is a symmetric positive definite 
%           matrix and computes:
%
%                P   = inv( chol( I+A*A' )' ) * A
%                B   = inv( chol( I-P*P' )' )
%                T   = chol(Sig)'* B
%                Phi = T* P *inv(T)
%
%           Reference: C.F. Ansley and R. Kohn, "A Note on Reparameterizing
%           a Vector Autoregressive Moving Average Model to Enforce 
%           Staionarity," J. of Statistical Computational Simulations,
%           1986, Vol 24, 99-106. 

%           Ellen McGrattan, 5-13-06
%           Revised, 5-30-06

n   = length(a);
i   = eye(n);
b   = chol(i+a*a')';
p   = inv(b)*a;
t   = chol(sig)'*b;
phi = t*p*inv(t);
