function like=ofn(xj)
%OFN: to use this file follow instructions listed under Sections I-IV
%
%  This file contains the objective function to minimized by mini.m.
%  It has been made easier to use by those interested in solving: *
%     
%  (P)   max E sum{ beta^t u(c1(t),c2(t),..,cN(t),v2~(t),param2) }
%        s.t.
%        c1(t) + i1(t)*IC(:,1) <= f1(v11(t),param11)
%          :         :                     :
%        cN(t) + iN(t)*IC(:,N) <= fN(v1N(t),param1N)    
%        X(t+1) = A*X(t)+B*u(t)+g*epsilon(t),  epsilon=N(0,Sigma)
%        param=[param1; param2],  
%        param1=[param11; param12; ... param1N]
%        y(t)=C*X(t)+w1(t)         agents' observables in Hidden state case
%        y1(t)=C1*X(t)+w1(t)  }    agents' observables in Multi-stage 
%        y2(t)=C2*X(t)+w2(t)  }    decision case
%        x(t)=[v11(t);v12(t);..;v1N(t);i1(t);i2(t);..;iN(t);v2~(t)]
%        x(t)=D*X(t)+E*u(t)    Note: v2~ is subset of v2 not in v11,..,v1N
%  * See Plspec.m for more details on variables, parameters, etc.
%
%  From (P), we get the following system:
%       
%        W(t)=AA*W(t-1)+eta(t),  where W,AA,eta depend on the agent's info. set
%                                eta=N(0,vcov)
%  Add:   
%        z(t)=CE*X(t)+v(t)  z(t): econometrician's observables at t
%                           v(t): measurement errors at t
%        v(t+1)=DD*v(t)+nu(t),  nu=N(0,noise)
%        zz(t)=z(t+1)-DD*z(t)=(CE*AA-DD*CE)*W(t)+nu(t)+CE*eta(t)= CC*W(t)+nnu(t)
%        
z=zxvar;
eval456=eval456+1;
diary off
%
% ----------------------------------------------------------------------------
%                     I.  SPECIFICATION OF PARAMETERS:
% Instructions: set "info", and all parameters as required in Plspec.m,
% instructions 2 and 3. Any parameters not used (e.g. if info=1, G is not) 
% should be set to 0. xj is the vector of parameters to be estimated.
% ----------------------------------------------------------------------------
%
% Parameters used by Kydland and Prescott (1982). 
%
a0=.5; beta=.99; delta=xj(1,1); eta=.1; gamma=-.5; phi1=.25; phi2=.25;
phi3=.25; phi4=.25; rho=xj(2,1); theta=.64; sigma=.0000028344535; v=4;
vxi1=.009^2; vxi2=.0018^2; vxi3=.009^2;
%
A=[1,0,0,0,0,0,0,0,0;
   0,rho,0,0,0,0,0,0,0;
   0,0,0,0,0,0,0,0,0;
   0,0,0,1-delta,0,0,1,0,0;
   0,0,0,0,0,0,0,0,0;
   0,0,0,0,0,1-eta,0,0,0;
   0,0,0,0,0,0,0,1,0;
   0,0,0,0,0,0,0,0,1;
   0,0,0,0,0,0,0,0,0];
B=[0,0,0;0,0,0;0,0,0;0,0,0;0,0,1;1,0,0;0,0,0;0,0,0;0,1,0];
D=[1 1 1 0 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0;
   0 0 0 0 0 0 0 0 0;
   0 0 0 0 -1 0 phi1 phi2 phi3;
   0 0 0 0 0 1 0 0 0];
E=[0 0 0;
   0 0 0;
   0 0 0;
   1 0 0;
   0 phi4 1;
   0 0 0];
g=zeros(9);
g(2,2)=1;
g(3,3)=1;
%
N=1;
IC=1;
param1=[sigma theta v]';
param11=param1;
param2=[a0 eta gamma]';
param=[param1; param2];
Sigma=zeros(9);
Sigma(2,2)=vxi1;
Sigma(3,3)=vxi2;
nov1=4;
%
info=3;
ytknown=0;
C=0;
V2=0;
V3=0;
G=[1 0;0 1;0 0];
H=[0 0 1]';
C1=[0 1 1 0 0 0 0 0 0];
C2=[0 1 1 0 0 0 0 0 0];
V21=vxi3;
V22=0;
V31=zeros(9,1);
V32=zeros(9,1);
hidden=[2 3]';
% ---------------------------------------------------------------------------
%           II.  SPECIFICATION OF STEADY STATES FOR x, X, u
% Instructions: for speed, give analytic solution of steady-state of non-
% stochastic version of model. If you can't, use steadyst.m.
% ---------------------------------------------------------------------------
%   Analytic solution of steady-states for Kydland-Prescott (1982)
%
int=(1-beta)/beta;
pk=phi1+phi2*(1+int)+phi3*(1+int)^2+phi4*(1+int)^3;
b1=.1;
sigma=(1+pk*(int+delta)/(int*b1^(v+1)))^(-1);
b2=1-sigma+sigma*b1^(-v);
b3=(((1-theta)*(1-sigma)*b2^(-(1-theta+v)/v))/(pk*(int+delta)))^(1/theta);
b4=b2^(-(1-theta)/v)*b3^(1-theta);
x(1)=1;
x(4)=(1+2*(1-(delta*b3/b4))*(a0*int+eta)/(theta*(int+eta)))^(-1);
x(2)=b3*x(4);
x(3)=b1*x(2);
x(5)=delta*x(2);
x(6)=(1/eta)*x(4);
x=x';
X=[1 0 0 x(2) x(3) x(6) x(5) x(5) x(5)]';
u=[x(4) x(5) x(3)]';
%
% ---------------------------------------------------------------------------
%                III.  RESTRICTIONS ON PARAMETER SPACE
% Instructions: set bounds on xj if there are any.  If xj goes out of the
% parameter space, set like=999999
% ---------------------------------------------------------------------------
pp = ( xj(1:2,1) <= zeros(2,1) );
pp=sum(pp);
if pp > .5
  like = 999999;
  diary on;
else;
  %
  % -------------------------------------------------------------------------
  %                        IV. SOLVE (P) 
  % -------------------------------------------------------------------------
  [AA,FF,vcov,XI]=mcplan(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));
  nh=length(hidden);
  lx=length(x);
  DEF=[D zeros(lx,3*nh)]-E*FF;
  [f df]=ffn1(x(1:nov1),param11);
  CE=[DEF(2:5,:);df'*DEF(1:nov1,:)+[f-df'*x(1:nov1); zeros(ns+3*nh-1,1)]' ];
  % -------------------------------------------------------------------------
  %                  VI. ADDING MEASUREMENT ERRORS  
  % ------------------------------------------------------------------------ 
  [RCE,col]=size(CE);
  noise=zeros(RCE);
  for i=1:RCE;
    noise(i,i)=.0001;
  end;
  DD=zeros(RCE);
  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 = dlsim(AA-K*CC,K,-CC,eye(cz),zz,XI);
  isut=inv(S)*U';
  ut=U';
  %
  like = .5*(rz-2)*log(det(S)) + .5*ut(:)'*isut(:);
  %
  diary on
end;
