function [k1o,k1e,so,se] = kf2stage(Ao,Ae,Co,Ce,V1o,V1e,V2o,V2e,V12o,V12e)
%  KF2STAGE calculates the kalman gains, ki=(Ai si Ci'+V12i) (Ci si Ci'+
%  V2i)^(-1), i=o,e, and the stationary covariance matrices si=
%  cov(x(t)| y(t-1),y(t-2),...;t=i,i in [o,e]) using the kalman filter for:
%  
%		x[t+1] = A[t]x[t] + Bu[t] + w1[t+1]   
%               y[t] = C[t]x[t] + Du[t] + w2[t]           
%               A[t]= { Ao, t odd
%                     { Ae, t even
%               C[t]= { Co, t odd
%                     { Ce, t even
%        
%                                        { [V1o   V12o;    t odd
%               E [w1(t+1)] [w1(t+1)]' = {  V12o' V2o ]
%                 [ w2(t) ] [ w2(t) ]    { [V1e   V12e;    t even
%                                        {  V12e' V2e ]
%
%  where x is the mx1 vector of states, u is the nx1 vector of controls, y 
%  is the px1 vector of observables, A[t] is mxm, C[t] is pxm, V1[t] is mxm, 
%  V2[t] is pxp, V12[t] is mxp.
%
%  INPUTS: Ao,Ae,Co,Ce,V1o,V1e,V2o,V2e,V12o,V12e, (or just Ao,Ae,Co,Ce,V1o,
%          V1e,V2o,V2e if there are no cross products)
%
%  OUTPUTS: 
%    (i) Kalman gains: ko, ke
%   (ii) Covariance matrices: so, se
%
%  NOTES: see also KFILTER
%
m=max(size(Ao));
[rc,cc]=size(Co);
if nargin==8; V12o=zeros(m,rc); V12e=V12o; end;
so=.01*eye(m);
se=.01*eye(m);
k0o= (Ao*so*Co'+V12o)/(V2o+Co*so*Co');
k0e= (Ae*se*Ce'+V12e)/(V2e+Ce*se*Ce');
dd=1;
it=1;
maxit=1000;
tol=1e-8;                             % reset tol if necessary
while (dd>tol & it<=maxit);
  so= Ae*se*(Ae'-Ce'*k0e') + V1e -V12e*k0e';
  se= Ao*so*(Ao'-Co'*k0o') + V1o -V12o*k0o';
  k1o= (Ao*so*Co'+V12o)/(V2o+Co*so*Co');
  k1e= (Ae*se*Ce'+V12e)/(V2e+Ce*se*Ce');
  dd=max(max(abs( [k1o-k0o;k1e-k0e]./(1+abs([k0o;k0e])) )));
  it=it+1;
  k0o=k1o;
  k0e=k1e;
end;
if it>=maxit; disp('WARNING: Iteration limit of 1000 reached in KF2STAGE'); end;
