function [Sigmaz]=hsmle(Theta,dummy);
% 
% The model to be solved is a version of Hansen and Sargent (1987).
%          
%          max E sum{ beta^t (log(c(t)) -a1*(n1(t)-n2(t)) -a2*n2(t)
%          s.t.
%              c(t) + i(t) = xi(t) k^theta {h1*(n1(t)-d/2*(n1(t)-n1(t-1))
%                            ^2)^(1-theta)+h2*n2(t)^(1-theta)} 
%              k(t+1)=(1-delta)/lambda k(t) + (1/lambda) i(t)
%              xi(t+1)=rho0+rho1*xi(t)+rho2*xi(t-1)+epsilon(t+1)
%
%  States:       [1 xi(t) xi(t-1) k(t) n1(t-1)]'
%  Controls:     [n1(t) n2(t) i(t)]'
%  Vars in u(.): [xi(t) k(t) n1(t-1) n1(t) n2(t) i(t)]'
%
Theta=Theta(:);
lb=[ -10;1e-8;1e-8;.2;.01;.0001;.6;-.5;.02;.0075;.96;1;1e-8;1e-8;1;1];
ub=[ 10;.8;.3;.8;10;10;1.1;.3;.03;10;.999;10000;.1;.1;10000;10000];
if any(lb>Theta) | any(ub<Theta) | Theta(2)+Theta(3)>1 | Theta(7)+Theta(8)>1;
  L=-1e+10; dL=-1e+10*ones(Theta);
  return
end;

% default values of the parameters
d=3.08; 
h1=.432; 
h2=.0639; 
theta=.365;
aa=3.62; 
delta=.0241; 
lambda=1.005186;
rho0=5.20; 
rho1=.7; 
rho2=.277;
sigma=1.64;
beta=.995; 
varm=[1000;.001;.001;1000;1000];
coefm=zeros(5,1);

% parameter vector
param=[d h1 h2 theta aa lambda rho0 rho1 rho2 delta sigma beta]';
param=[param;varm;coefm];

% parameters to be estimated
ind=[1:5,7:17];
d=Theta(1); h1=Theta(2); h2=Theta(3); theta=Theta(4); aa=Theta(5);
rho0=Theta(6); rho1=Theta(7); rho2=Theta(8); 
delta=Theta(9); sigma=Theta(10); beta=Theta(11); varm=Theta(12:16);
%coefm=Theta( );
param(ind)=Theta;
np=length(ind);
betain=11;

% A,B,D,E,Sigma,ct
a1=-aa*log(1-h1);
a2=-aa*log(1-h1-h2);
A=[1     0    0   0      0;
   rho0 rho1 rho2 0      0;
   0     1    0   0      0;
   0     0    0 (1-delta)/lambda  0;
   0     0    0   0      0];
B=[0 0 0;0 0 0;0 0 0;0 0 1/lambda;1 0 0];
D=[0 1 0 0 0;
   0 0 0 1 0;
   0 0 0 0 1;
   0 0 0 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];
Sigma=zeros(5);
Sigma(2,2)=sigma^2;
ct=1;


% steady states
lbd=(lambda-beta*(1-delta))/(beta*theta);
ti=1/theta;
tt=1-ti;
x(1,1)=rho0/(1-rho1-rho2);
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(1) x(2) x(3)]';
u=[x(4) x(5) x(6)]';
ns=length(X);
nc=length(u);
nx=length(x);


% data Y(t)=[k(t),n1(t),n2(t),i(t),c(t)]';
X0=X;
%load /home/erm/data/hsdata3
T=length(zxvar);
tem=lambda.^([1:T]');
%Y=[zxvar(:,1)./tem,zxvar(:,2),zxvar(:,2).*(zxvar(:,3)-13.045*80*h1)./ ..
%   (13.045*80*h2),zxvar(:,5)./tem,zxvar(:,4)./tem];
Y=[zxvar(:,1)./tem,zxvar(:,2),zxvar(:,3),zxvar(:,5)./tem,zxvar(:,4)./tem];
Y=Y(2:T,:)-Y(1:T-1,:)*diag(coefm);
T=T-1;
Omega=diag(varm);
Coefm=diag(coefm);
% see C below

% retrieve b=dr(x,theta)/x 
%          q=d^2r(x,theta)/dx^2
% plus     cub=d^3r(x,theta)/dx^3
%          du=dr(x,theta)/dtheta
%          db=db/dtheta 
%          dq=dq/dtheta   if derivatives required
if nargout > 1;
  [ux,b,q,cub,du,db,dq]=hsu4(x,param(1:5));
  du=[du(:);zeros(np-5,1)];
  db=[db,zeros(6,np-5)];
  dq=[dq,zeros(36,np-5)];
else;
  [ux,b,q]=hsu4(x,param(1:5));
end;
b=b(:);

% construct Q,R,W 
Q0=zeros(ns);
Q0(ct,ct)=ux-b'*x+x'*q*x/2;
d=(b'-x'*q)';
Q1=zeros(ns);
Q1(ct,:)=d'*D;
Q1(:,ct)=(d'*D)';
Q1=Q1/2;
Q1(ct,ct)=2*Q1(ct,ct);
Q=Q0+Q1+D'*q*D/2;
R=E'*q*E/2;
Q2=zeros(ns,nc);
Q2(ct,:)=(d'*E)/2;
W=Q2+D'*q*E/2;


% make sure R is invertible and construct undiscounted, no cross-products
% versions of A,B,Q,R 
if rank(R)>=length(R);
  An=sqrt(beta)*(A-B*(R\W'));
  Bn=sqrt(beta)*B;
  Qn=Q-W*(R\W');
else;
  disp('WARNING: R is not invertible at Theta equal:')
  Theta
  return
end;

% compute Riccati matrix and decision rules 
[Fn,P]=double(An',Bn',Qn,R);
Fn=Fn';
Aon=An-Bn*Fn;

% construct matrices for original problem
F=Fn+(R\W');
Ao=A-B*F;

% compute innovations for likelihood function
[tem1,tem2]=hsf4(x,param(1:4));
C=[0,0,0,1,0;
   -F;
   tem2'*(D-E*F)+[tem1-tem2'*x,zeros(1,ns-1)]];


Cbar=C*Ao-Coefm*C;
Obar=Omega+C*Sigma*C';
[K,S]=kfilter(Ao,Cbar,Sigma,Obar,Sigma*C');
Sigmaz=Obar+Cbar*S*Cbar';

