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;