function [F,C,Ao]=hs2QRW(p);
d=p(1); h1=p(2); h2=p(3); theta=p(4); aa=p(5);
rho0=p(7); rho1=p(8); rho2=p(9); delta=p(10);
lambda=p(6); sigma=p(11); beta=p(12);
varm=p(13:17);
coefm=p(18:22);
load /home/erm/data/hsdata3



% 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];
DY=Y(2:T,:)-Y(1:T-1,:)*diag(coefm);
T=T-1;
Omega=diag(varm);
G=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
  [ux,b,q]=hsu4(x,p(1:5));
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;
Wp=W';

% 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,p(1:4));
C=[0,0,0,1,0;
   -F;
   tem2'*(D-E*F)+[tem1-tem2'*x,zeros(1,ns-1)]];


Cbar=C*Ao-G*C;
Obar=Omega+C*Sigma*C';
[K,S]=kfilter(Ao,Cbar,Sigma,Obar,Sigma*C');
AA=Ao-K*Cbar;
Sigmaz=Obar+Cbar*S*Cbar';
Sigmazi=inv(Sigmaz);
t0=3;
Xt(t0,:)=X0';
innov(t0,:)=DY(t0,:)-X0'*Cbar';
sum1=innov(t0,:)'*innov(t0,:);
sum2=Xt(t0,:)'*innov(t0,:);
sum6=Y(t0,:)'*innov(t0,:);
for i=t0+1:T;
  Xt(i,:)=Xt(i-1,:)*Ao'+innov(i-1,:)*K';
  innov(i,:)=DY(i,:)-Xt(i,:)*Cbar';
  sum1=sum1+innov(i,:)'*innov(i,:);
  sum2=sum2+Xt(i,:)'*innov(i,:);
  sum6=sum6+Y(i,:)'*innov(i,:);
end;
sum1=sum1/(T-t0+1);
sum2=sum2/(T-t0+1);
sum6=sum6/(T-t0+1);

M1=Sigmazi-Sigmazi*sum1*Sigmazi;
L= (T-t0+1)*(log(det(Sigmaz))+trace(Sigmazi*sum1)) + ..
   1.5*log(lambda)*(T-t0+1)*(T-t0+2);
%L=L/1100;


