function [L,ft]=lqmle(theta,t)
%LQMLE  computes log-likelihood from linear-quadratic control problem
%       [L,ft]=lqmle(theta,t) starts with a LQ problem of the form:
%       
%              max E{   sum  beta^t X[t]' Q X[t] + u[t]' R u[t]  |X[0] }
%            {u[t]}   t=0:inf         
% 
%              subject to   X[t+1] = A X[t] + B u[t] + eps[t+1]
%
%       with the user (by editing this file) defining A,B,Q,R,beta as
%       a function of input theta and supplying dA/dtheta, dB/dtheta,
%       dQ/dtheta, dR/dtheta, dbeta/dtheta.  Also supplied by the
%       user are the steady states for X and u.
% 
%                 
%       L= -T/2*log(|S|)-1/2 sum  v(t)'S^(-1) v(t)
%                           t=1:T 
%




z=zxvar;


% ----------------------------------------------------------------------------
%                     I.  SPECIFICATION OF PARAMETERS:
% ----------------------------------------------------------------------------

theta=theta(:);

beta=.99; delta=.023; lambda=1.005186;
aa=theta(1); h1=theta(2); h2=theta(3); d=theta(4); rho0=theta(5);
a1=-aa*log(1-h1); a2=-aa*log(1-h1-h2);
rho1=theta(6); sigma=theta(7); thet=theta(8); 
noise=diag(theta(9:13));
DD=zeros(5);
if nargin==1; t=0; end;


% ----------------------------------------------------------------------------
%              II. SPECIFICATION OF A,B,D,E, AND DERIVATIVES 
% ----------------------------------------------------------------------------


  A=[1     0    0      0;
     rho0 rho1  0      0;
     0     0   (1-delta)/lambda  0;
     0     0    0      0];
  B=[0 0 0;0 0 0;0 0 1/lambda;1 0 0];
  D=[0 1 0 0;
     0 0 1 0;
     0 0 0 1;
     0 0 0 0;
     0 0 0 0;
     0 0 0 0];
  E=[0 0 0;
     0 0 0;
     0 0 0;
     1 0 0;
     0 1 0;
     0 0 1];

  % ---------------------------------------------------------------------------
  %           III.  SPECIFICATION OF STEADY STATES FOR x, X, u 
  % ---------------------------------------------------------------------------
  lbd=(lambda-beta*(1-delta))/(beta*theta);
  ti=1/theta;
  tt=1-ti;
  x(1,1)=rho0/(1-rho1);
  x(2,1)=((h1^ti*a1^tt+h2^ti*(a2-a1)^tt)/lbd)^(-1/tt)*(1-theta)*x(1)^ ..
         (1/(1-theta))/(lbd-lambda+1-delta);
  x(3,1)=(a1*(lbd-lambda+1-delta)*x(2)^(1-theta)/(x(1)*h1*(1-theta)))^(-ti);
  x(4,1)=x(3);
  x(5,1)=(h2*a1/(h1*(a2-a1)))^ti*x(3);
  x(6,1)=(lambda-1+delta)*x(2);
  X=[1 x(1) x(2) x(3)]';
  u=[x(4) x(5) x(6)]';

  % -------------------------------------------------------------------------
  %                        IV. SOLVE (P) 
  % -------------------------------------------------------------------------
  [AA,FF,vcov,XI]=tmcplan2(A,B,D,E,g,beta,param1,param2,Sigma,nov1,x,X,u, ..
                      info,ytknown,C,C1,C2,V2,V21,V22,V3,V31,V32,G,H,hidden);
  % -------------------------------------------------------------------------
  %                     V. CREATE MATRIX CE
  % Instructions: "construct.doc" gives a full description and all formulas
  % needed to set CE (analytically!).
  % -------------------------------------------------------------------------
  ns=max(size(A));
  lx=length(x);
  DEF=D-E*FF;
  [f,df]=ffn1(x(1:nov1),param11);
  CE=[DEF(2,:);DEF(4:5,:);df'*DEF(1:nov1,:)+[f-df'*x(1:nov1);zeros(ns-1,1)]'- ..
      DEF(6,:);DEF(6,:)];
  %
  %  TRANSFORMATION OF DATA:
  %
  [rz,cz]=size(z);
  AVH=z(:,3);
  ztem=z(:,2).*(z(:,3)-13.045*80*h1)/(13.045*80*h2);
  tem=lambda.^([1:rz]');
  z=[z(:,1)./tem z(:,2) ztem z(:,4)./tem z(:,5)./tem];
  % -------------------------------------------------------------------------
  %                  VI. ADDING MEASUREMENT ERRORS  
  % Instructions: specify matrices "noise" and "DD" --definitions given above  
  % ------------------------------------------------------------------------- 
  [RCE,col]=size(CE);
  CC=CE*AA-DD*CE;
  [K,s]=kfilter(AA,CC,vcov,noise+CE*vcov*CE',vcov*CE'); 
  S=CC*s*CC'+noise+CE*vcov*CE';
  %
  [rz,cz]=size(z);
  zz=z(3:rz,:)-z(2:rz-1,:)*DD';
  [U,Xhat] = dlsim(AA-K*CC,K,-CC,eye(cz),zz,XI);
  isut=S\U';
  ut=U';
  %
  jacob=log( z(3:rz,2)./(13.045*80*h2*AVH(3:rz,1)) );
  like = .5*(rz-2)*log(det(S)) + .5*ut(:)'*isut(:)  ..
          -sum(jacob)+penalty*(h1+h2>1);
  if t>0; 
    ft=.5*log(det(S))+.5*U(t,:)*inv(S)*U(t,:)' - jacob(t,1); 
  else; 
    ft=0; 
  end;

like=-like;
