function L = mleq(Theta,fixedpar);
%MLE     Log-likelihood function for the multi-sector model 
%        with intangible investments.  See notes for details.
%

%        Ellen R. McGrattan, 6-7-15
%

global ZVAR % = observables
[T,nz]        = size(ZVAR);
%---------------------------------------------------------------------
% 1. Fixed parameters 
%
n             = fixedpar(1);
o             = ones(n,1);
alpha         = fixedpar(2);
beta          = fixedpar(3);
cbar          = fixedpar(4);
deli          = fixedpar(5:4+n);
delt          = fixedpar(5+n:4+2*n);
g             = fixedpar(5+2*n:4+3*n);
gamma         = zeros(n);
gamma(:)      = fixedpar(5+3*n:4+n*n+3*n);
gn            = fixedpar(5+n*n+3*n);
gz            = fixedpar(6+n*n+3*n);
nu            = zeros(n);
nu(:)         = fixedpar(7+n*n+3*n:6+2*n*n+3*n);
omega         = fixedpar(7+2*n*n+3*n:6+2*n*n+4*n);
p1            = fixedpar(7+2*n*n+4*n);
phi           = fixedpar(8+2*n*n+4*n:7+2*n*n+5*n);
psi           = fixedpar(8+2*n*n+5*n);
rho           = fixedpar(9+2*n*n+5*n);
tauc          = fixedpar(10+2*n*n+5*n);
taud          = fixedpar(11+2*n*n+5*n);
tauh          = fixedpar(12+2*n*n+5*n);
tauk          = fixedpar(13+2*n*n+5*n);
taup          = fixedpar(14+2*n*n+5*n);
taux          = fixedpar(15+2*n*n+5*n);
theta         = fixedpar(16+2*n*n+5*n:15+2*n*n+6*n);
z1            = fixedpar(16+2*n*n+6*n:15+2*n*n+7*n);
z2            = fixedpar(16+2*n*n+7*n:15+2*n*n+8*n);
zeta          = zeros(n);
zeta(:)       = fixedpar(16+2*n*n+8*n:15+3*n*n+8*n);
z             = fixedpar(16+3*n*n+8*n);
xguess        = fixedpar(17+3*n*n+8*n:16+3*n*n+11*n);
pvec          = fixedpar(17+3*n*n+11*n:16+3*n*n+12*n);
qvec          = fixedpar(17+3*n*n+12*n:16+3*n*n+13*n);
yvec          = fixedpar(17+3*n*n+13*n:16+3*n*n+14*n);
xvec          = fixedpar(17+3*n*n+14*n:16+3*n*n+15*n);
agghss        = fixedpar(17+3*n*n+15*n);

param         = [n;alpha;beta;cbar;deli;delt;g;gamma(:);gn;gz;nu(:);omega;
                 p1;phi;psi;rho;tauc;taud;tauh;tauk;taup;taux;theta;z1;z2;
                 zeta(:);z];

%---------------------------------------------------------------------
% 2. Initialize parameters that may be estimated
%

ns            = 3*n+7;
P             = zeros(ns);
Q             = zeros(ns);
Dd            = zeros(nz);
Rr            = zeros(nz);
ii            = 1:2*n;
P(ii,ii)      = .9*eye(2*n);
ii            = 1:n;
Q(ii,ii)      = .01*eye(n);
ii            = n+1:2*n;
Q(ii,ii)      = .0001*eye(n);


P(1,1)        = .995;
P(2,2)        = .995;
P(4,4)        = .995;
P(5,5)        = .995;
P(6,6)        = .995;
P(8,8)        = .995;
P(20,20)      = .995;
P(26,26)      = .995;
P(ns,ns)      = .995;
Q(ns,ns)      = .0104;
%
% Turn off cyclical variation in taul
%
P(48,48)      = .95;
Q(48,48)      = 0.00001;
wgt           = ones(n,1);

%---------------------------------------------------------------------
% 3. Assign elements of Theta to underlying parameters.  Use the vector
%    ind to indicate which parameters are to be estimated and which
%    are to be fixed in estimation.  For example, if the vector of
%    parameters is [beta,gamma,delta] and gamma is to be fixed during
%    estimation then set:
%       beta  = Theta(1);
%       delta = Theta(2);
%       ind   = [1,0,1];
%
%    For the current example, the ordering of the parameter vector is:
%      diag(P), diag(Q),     diag(Dd),        diag(Rr)
%      1:ns     ns+1:2*ns,   2*ns+1:2*ns+nz   2*ns+nz+1:2*ns+2*nz

ind           = zeros(2*ns+2*nz+4,1);
ii            = [3, 7,9:14, 24,  ...
                 ns+1:ns+14, ns+20,ns+24,ns+26,  ...
                 2*ns+2*nz+1:2*ns+2*nz+3+14];
ind(ii)       = 1;

P(3,3)        = Theta(1);
P(7,7)        = Theta(2);
P(9,9)        = Theta(3);
P(10,10)      = Theta(4);
P(11,11)      = Theta(5);
P(12,12)      = Theta(6);
P(13,13)      = Theta(7);
P(14,14)      = Theta(8);
P(24,24)      = Theta(9);
Q(1,1)        = Theta(10);
Q(2,2)        = Theta(11);
Q(3,3)        = Theta(12);
Q(4,4)        = Theta(13);
Q(5,5)        = Theta(14);
Q(6,6)        = Theta(15);
Q(7,7)        = Theta(16);
Q(8,8)        = Theta(17);
Q(9,9)        = Theta(18);
Q(10,10)      = Theta(19);
Q(11,11)      = Theta(20);
Q(12,12)      = Theta(21);
Q(13,13)      = Theta(22);
Q(14,14)      = Theta(23);
Q(20,20)      = Theta(24);
Q(24,24)      = Theta(25);
Q(26,26)      = Theta(26);
Q(5,20)       = Theta(27);
Q(9,24)       = Theta(28);
Q(11,26)      = Theta(29);
wgt(1)        = Theta(30);
wgt(2)        = Theta(31);
wgt(3)        = Theta(32);
wgt(4)        = Theta(33);
wgt(5)        = Theta(34);
wgt(6)        = Theta(35);
wgt(7)        = Theta(36);
wgt(8)        = Theta(37);
wgt(9)        = Theta(38);
wgt(10)       = Theta(39);
wgt(11)       = Theta(40);
wgt(12)       = Theta(41);
wgt(13)       = Theta(42);
wgt(14)       = Theta(43);


%---------------------------------------------------------------------
% 4. Put bounds on Theta if there are values of the parameters for
%    which L cannot be computed. (Optional)
%---------------------------------------------------------------------
Ub            = [.999*ones(ns,1);1000*ones(ns,1);
                 .999*ones(nz,1);1000*ones(nz,1); 1000*ones(3+n,1)];
Lb            = -Ub;
Ub            = Ub(find(ind));
Lb            = Lb(find(ind));
if (any(Theta>Ub) | any(Theta<Lb));
  L           = 1e+20; dL=1e+20*ones(size(Theta));
  return
end;
penalty       = 500000*max(max(abs(eig(P)))-.995,0).^2;

%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 5.  Compute equilibrium
%---------------------------------------------------------------------
%
% 5a. Steady state:
%

x1            = ss_fullmodel(xguess,param);
o             = ones(n,1);
betg          = beta*(1+gz)^(-alpha);
gam           = sum(gamma)';
lshare        = 1-theta-phi-gam;
thet1         = 1-theta;
aggcnrm       = x1(1);
pvec          = [p1;x1(2:n)];
h2h1          = x1(n+1:2*n);
h             = x1(2*n+1:3*n);
pmvec         = gam.*prod(((pvec*o')./(gamma+1e-15)).^(gamma./(o*gam')))';
ptvec         = prod(((pvec*o')./(zeta+1e-15)).^zeta)';
aggp          = sum(omega.^rho.*pvec.^(1-rho))^(1/(1-rho));
cnrm          = (aggp*omega./pvec).^rho * aggcnrm;
c             = cnrm*cbar;
aggc          = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
aggl          = 1-sum(h);
m1y           = gam.*pvec./pmvec;
h1            = h./(1+h2h1);
h2            = h-h1;
qxipy         = h2h1;
xiy           = qxipy.^(1-phi).*(z2./z1).^lshare;
qp            = qxipy./xiy;
qvec          = qp.*pvec;
pivec         = prod(((qvec*o')./(nu+1e-15)).^nu)';
kiy           = betg*phi.*(qvec.*xiy+pvec)./((1-betg*(1-deli)).*pivec);
ret           = ((1+taux)*ptvec./pvec-betg*(1-(1-taup)*(delt+tauk) + ...
                (1-delt).*((1+taux)*ptvec./pvec-1)))/(betg*(1-taup));
k1ty          = theta./ret;
y             = (k1ty.^theta.*kiy.^phi.*m1y.^gam).^(o./lshare).*z1.*h1*z;
k1t           = k1ty.*y;
k2t           = h2h1.*k1t;
kt            = k1t+k2t;
m1            = m1y.*y;
xi            = xiy.*y;
ki            = kiy.*y;
qxi           = qvec.*xi;
qxip          = qxi./pvec;
prodxt        = ((1+gz)*(1+gn)-1+delt).*kt;
prodxi        = ((1+gz)*(1+gn)-1+deli).*ki;
m1mat         = gamma.*((o./pvec)*pvec').*(o*y');
m2mat         = gamma.*((o./pvec)*pvec').*(o*qxip');
xtmat         = zeta.*((o./pvec)*ptvec').*(o*prodxt');
ximat         = nu.*((o./qvec)*pivec').*(o*prodxi');
Sbar          = [log(z1);log(z2);log(g);tauc;taud;tauh;tauk;taup;taux;log(z)];

%
% 5b. Call subroutine with residuals:
%
m1matplus     = max(m1mat,1e-10);
m2matplus     = max(m2mat,1e-10);
v             = [log(kt);log(ki);log(kt);log(ki);log(kt);log(ki);Sbar;Sbar;
                log(c);log(h1);log(h2);log(m1matplus(:));log(m2matplus(:));
                log(c);log(h1);log(h2);log(m1matplus(:));log(m2matplus(:))];
del           = max(abs(v)*1e-4,1e-8);
nres          = 8*n+4*n*n;
nv            = length(v);
[res,xxs]     = res_fullmodel(v,[param;pvec;qvec;yvec;xvec;1-aggl;h2([5,9,11]);wgt]);
dR            = zeros(nres,nv);
dV            = zeros(length(xxs(:)),nv);

for i=1:nv;
  vp          = v;
  vm          = v;
  vp(i)       = vp(i)+del(i);
  vm(i)       = vm(i)-del(i);
  [rp,xp]     = res_fullmodel(vp,[param;pvec;qvec;yvec;xvec;1-aggl;h2([5,9,11]);wgt]);
  [rm,xm]     = res_fullmodel(vm,[param;pvec;qvec;yvec;xvec;1-aggl;h2([5,9,11]);wgt]);
  dR(:,i)     = (rp-rm)/(2*del(i));
  dV(:,i)     = (xp-xm)/(2*del(i));
end;

i1            = 1:6*n+4*n*n;
i2            = 6*n+4*n*n+1:8*n+4*n*n;
j1            = 1:6*n+2*ns;
j2            = 6*n+2*ns+1:nv;
Coef          = -inv(dR(i1,j2))*dR(i1,j1);
dRnew         = dR(i2,j1)+dR(i2,j2)* Coef;
dVnew         = dV(:,j1)+dV(:,j2)* Coef;


%
% 5c. Solution:  Find deterministic part of solution
%
%         Residuals are A1*[X_{t+1}] + A2*[X_t]  + shock terms    = 0
%                          [Z_{t+1}]      [Z_t]
%
%                           X_t     = [kt[t],ki[t]],
%                           Z_t     = [X_t,h1[t],h2[t],y[t]]
%
%

nX            = 2*n;
nZ            = 2*n;
nXZ           = 4*n;
a1            = dRnew(:,1:2*n);
a2            = dRnew(:,2*n+1:4*n);
a3            = dRnew(:,4*n+1:6*n);
a4            = dRnew(:,6*n+1:6*n+ns);
a5            = dRnew(:,6*n+ns+1:6*n+2*ns);
A1            = [eye(2*n),zeros(2*n);zeros(2*n),a3];
A2            = [zeros(2*n),-eye(2*n);a1,a2];

[V,L]         = eig(A2,-A1);

L             = diag(L);
[D,I]         = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) > 1) ) )
  i      = sum(D<=1);
  disp(sprintf( ...
  'WARNING: program assumes %g eigenvalues inside unit circle, found %g',nX,i));
end;
V11           = V(1:nX,I(1:nX));
V21           = V(nX+1:nXZ,I(1:nX));
W             = real([V11*diag(L(I(1:nX)))/V11;
                      V21/V11]);
A             = W(1:2*n,:);
C             = W(2*n+1:4*n,:);

%
% Step 5d. Solution: Find the stochastic part of the solution (ie., B,D)
%
%         Residuals (k in logs):
%
%           E_t{ a1 k[t]+ a2 k[t+1] + a3 k[t+2] +
%                         a4 S[t] + a5 S[t+1] } = 0,
%
%           where k =[ln(kt);ln(ki)]
%                 S =[ln(z1),ln(z2),ln(g),tauc,taud,tauh,tauk,taup,taux,ln(z)]
%                     (3n+7 x 1)
%
%           ==>  (a1+a2*A+a3*A*A)*k
%              + (a2*B+a3*A*B+a3*B*P+a4+a5*P)*S = 0
%
%           Note: need vec(A+B)=vec(A)+vec(B), vec(ABC)=(C' kron A)*vec(B)
%
%

resk          = a1+a2*A+a3*A*A;

if (any(abs(resk)>1e-8));
  disp('Check deterministic part of the solution')
end;

I             = eye(ns);
tem           =-[kron(I,a2+a3*A)+kron(P',a3)]\vec(a4+a5*P);
B             = zeros(2*n,ns);
B(:)          = tem(1:2*n*ns);
D             = B;

ii            = [1:3*n+2*n*n]; % indices for c[t],h1[t],h2[t],m1mat[t],m2mat[t]
C2            = Coef(ii,1:2*n)+Coef(ii,2*n+1:4*n)*A;
D2            = Coef(ii,6*n+1:6*n+ns)+Coef(ii,2*n+1:4*n)*B;

%
% Variables used in the observer equation:
%

C2b           = dVnew(:,1:2*n)+dVnew(:,2*n+1:4*n)*A;
D2b           = dVnew(:,6*n+1:6*n+ns)+dVnew(:,2*n+1:4*n)*B;

%
% 5e. State-space system:   X[t+1] = G X[t] + H eps[t+1]
%                           Y[t]   = J X[t] + ome[t]
%                           ome[t] = D ome[t-1] + eta[t],  E eta eta' ~ N(0,R)
%

A0            = (eye(2*n)-A)*log([kt;ki])-B*Sbar;
C0            = log([c;h1;h2;m1matplus(:);m2matplus(:)])-C2*log([kt;ki])-D2*Sbar;
C0b           = xxs-C2b*log([kt;ki])-D2b*Sbar;
P0            = (eye(ns)-P)*Sbar;
G             = [A,B,A0;zeros(3*n+7,2*n),P,P0;zeros(1,5*n+7),1];
H             = [zeros(2*n,3*n+7);Q;zeros(1,3*n+7)];
J             = [C2,D2,C0];
X0            = [log(kt);log(ki);Sbar;1];
J2            = [C2b,D2b,C0b];   % variables in residual file


X0=[ ...
    -7.031786900000000e+00
    -4.429751700000000e+00
    -6.094775400000000e+00
    -6.451456100000000e+00
    -4.445784400000000e+00
    -6.438274900000000e+00
    -6.296600500000000e+00
    -6.731027300000000e+00
    -6.124995100000000e+00
    -3.752860100000000e+00
    -6.340284500000000e+00
    -6.413406700000000e+00
    -6.252717400000000e+00
    -6.838779800000000e+00
    -8.205554700000000e+00
    -9.766190399999999e+00
    -8.100303500000001e+00
    -7.703470800000000e+00
    -6.268378500000000e+00
    -1.562797000000000e+00
    -5.424832000000000e+00
    -6.164947700000000e+00
    -7.723203400000000e+00
    -3.721536800000000e+00
    -5.254491500000000e+00
    -4.297137700000000e+00
    -6.651353600000000e+00
    -6.273974000000000e+00
    -7.219698100000000e+00
    -9.191427700000000e+00
     2.340021700000000e+00
     2.133227000000000e+00
     2.876830200000000e+00
     6.208896200000000e-01
    -5.753050199999999e-01
     5.235052100000001e-01
     4.864037600000000e-01
    -9.295717400000000e-02
    -1.757904400000000e+00
    -1.781703800000000e+00
    -4.747424100000000e-01
    -6.385185000000000e-01
    -2.223075400000000e-01
     4.480574400000000e-02
     1.520635000000000e-01
    -1.271630700000000e-10
     2.250315000000000e-07
    -9.283079000000000e-09
    -3.067337600000000e-08
     1.470262200000000e+01
    -2.388825300000000e-07
    -6.071379600000000e-08
    -3.402443500000000e-09
     5.730735300000000e+00
    -1.958195600000000e-05
    -2.586646000000000e+00
    -8.526113200000000e-07
    -5.056565100000000e-09
    -6.910572100000000e-08
     7.614837000000000e-07
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
    -3.453877600000000e+01
     6.500000000000000e-02
     1.440000000000000e-01
     4.465999800000000e-01
     3.000000000000000e-03
     3.300000000000000e-01
                         0
     2.148136700000000e-01
     1.000000000000000e+00];
X0(78)=.45;


%---------------------------------------------------------------------
% 6. Specify observables (Z). 
%
Y             = ZVAR;
Ybar          = Y(2:T,:)-Y(1:T-1,:)*Dd';
T             = T-1;
Jbar          = J2*G-Dd*J2;
Rbar          = Rr+J2*H*H'*J2';
[K,Sigma]     = kfilter(G,Jbar,H*H',Rbar,H*H'*J2');
Omega         = Rbar+Jbar*Sigma*Jbar';
Omegai        = inv(Omega);
Xt(1,:)       = X0';
innov(1,:)    = Ybar(1,:)-X0'*Jbar';
for i=2:T;
  Xt(i,:)     = Xt(i-1,:)*G'+innov(i-1,:)*K';
  innov(i,:)  = Ybar(i,:)-Xt(i,:)*Jbar';
end;
sum1          = innov(1:T,:)'*innov(1:T,:)/T;
L             = 0.5*(T*(log(det(Omega))+trace(Omegai*sum1))+penalty);

