										
function[K,S,ST]=doublex(A,C,Q,R,c,sig)
%DOUBLEX.M
%function[K,S,ST]=doublex(A,C,Q,R,c,sig)
%  This program uses the "doubling algorithm" to solve the
%  Riccati matrix difference equations associated with the
%  undiscounted quadratic-Gauusian linear optimal control problems.
%  The control problem has the form
%
%    S(t) = min over u(t) of
%    x(t)' Q x(t) + u(t)' R u(t) + (2/sig)log E(t) exp (sig/2) S(x(t+1)),
%
%  subject to 
%
%     x(t+1) = A' x(t) + C' u(t) + c w(t+1),
% 
%  where w(t+1) is a Gaussian martingale difference sequence 
%  with unit covariance matrix. The program returns the steady
%  state  value function in S.  The optimal control law is
%  u(t) = -K' * x(t).  The program also returns ST, which
%  is the quadratic form in E(t) exp S(x(t+1)).
%  The matrices S and ST should be positive semi-definite

a0=A';
b0=C'*(R\C)-sig*c*c';          %  This is the J of Whittle.
%                       Eventually add a check on this




eds=.0001;    %tolerance
% Check if b0 is positive semi-definite and symmetric
if any(eig(b0) < -eds) | (norm(b0'-b0,1)/norm(b0,1) > eds)
        disp('b0  must be symmetric and positive semi-definite')
end

b0




%[t1,t2]=eig(b0);
%lammax=max(diag(t2));
%if lammax <= 0
%disp('warning: C^T R\C -sig c c^T is not positive definite')
%return
else
end
g0=Q;
tol=1e-15;
dd=1;
ss=max(size(A));
v=eye(ss);
										
while dd>tol
										
a1=a0 *((v+b0*g0)\a0);
b1=b0+a0*((v+b0*g0)\(b0*a0'));
g1=g0+a0'*g0*((v+b0*g0)\a0);
g1t=(v-sig*g1*c*c')\g1;
g0t=(v-sig*g0*c*c')\g0;
k1=A*g1t*C'/(C*g1t*C'+R);
k0=A*g0t*C'/(C*g0t*C'+R);
dd=max(max(abs(k1-k0)));
a0=a1;
b0=b1;
g0=g1;
end
K=k1;S=g1;ST=g1t;
