function a = istationary(phi,sig)
%ISTATIONARY is the inverse map of stationary.m which maps a general square
%            matrix to one with eigenvalues inside the unit circle.
%
%            A = istationary(Phi,Sig) takes the matrix Phi, which has
%            eigenvalues inside the unit circle, and Sig, which is a
%            symmetric positive definite matrix and computes:
%
%                  T = chol(V)',  V = Phi *V* Phi' + Sig
%                  P = inv(T)* Phi *T
%                  B = inv( chol( I-P*P' )' )
%                  A = B * P
%
                                                                                
%            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(phi);
i = eye(n);
v = dlyap(phi,sig);
t = chol(v)';
p = inv(t)*phi*t;
b = inv(chol(i-p*p')');
a = b*p;


