function f = trfoc(Y,param);

%
% Endogenous variables
%

nper   = param(2);
X      = zeros(33,nper);
X(:)   = Y;

l1fp   = X(1,:)';
l1mp   = X(2,:)';
l1fs   = X(3,:)';
l1ms   = X(4,:)';
l2fp   = X(5,:)';
l2fs   = X(6,:)';
l2ms   = X(7,:)';
c1fp   = X(8,:)';
c1mp   = X(9,:)';
c1fs   = X(10,:)';
c1ms   = X(11,:)';
c2fp   = X(12,:)';
c2mp   = X(13,:)';
k1pn   = X(14,:)';
k1fsn  = X(15,:)';
k1msn  = X(16,:)';
k2pn   = X(17,:)'; 
k2fsn  = X(18,:)';
k2msn  = X(19,:)';
hfpn   = X(20,:)';
hmpn   = X(21,:)';
hfsn   = X(22,:)'; 
hmsn   = X(23,:)';
efpn   = X(24,:)';
empn   = X(25,:)';
efsn   = X(26,:)'; 
emsn   = X(27,:)';
T      = X(28,:)';
Td     = X(29,:)';
mufp   = X(30,:)';
mump   = X(31,:)';
mufs   = X(32,:)';
mums   = X(33,:)';

%
% Parameters
%

np     = param(3);
nfs    = param(4);
nms    = 1-2*np-nfs;
beta   = param(5);
delh   = param(6);
dele   = param(7);
delk   = param(8);
kap1   = param(9);
zet1   = param(10);
kap2   = param(11);
zet2   = param(12);
b1f    = param(13);
b2f    = param(14);
b1m    = param(15);
b2m    = param(16);
lam    = param(17);
theta1 = param(18);
theta2 = param(20);
util   = param(22);
upar   = param(23:26);
gpar   = param(27:29);
g      = param(30);
tauk   = param(31);
Tp     = param(36)*T;
Tfs    = param(37)*T;
Tms    = T-Tp-Tfs;
Tdp    = param(38)*Td;
Tdfs   = param(39)*Td;
Tdms   = Td-Tdp-Tdfs;
gama   = param(40);
gamp   = param(41);
beth   = beta*(1+gama)^((1-upar(3))*(1-upar(4))-1);

mtech  = param(46+0*nper:45+1*nper);
htech  = param(46+1*nper:45+2*nper);
taulfp = param(46+2*nper:45+3*nper);
taulmp = param(46+3*nper:45+4*nper);
taulfs = param(46+4*nper:45+5*nper);
taulms = param(46+5*nper:45+6*nper);
qk1    = param(46+6*nper:45+7*nper);
qk2    = param(46+7*nper:45+8*nper);
qh     = param(46+8*nper:45+9*nper);
qe     = param(46+9*nper:45+10*nper);
grate  = ((1+gama)*(1+gamp)).^[0:nper-1]';
gratei = ones(nper,1)./grate;
csfp   = param(46+10*nper)*gratei;
csmp   = param(47+10*nper)*gratei;
csfs   = param(48+10*nper)*gratei;
csms   = param(49+10*nper)*gratei;
tsfp   = min(param(50+10*nper),nper+1);
tsmp   = min(param(51+10*nper),nper+1);
tsfs   = min(param(52+10*nper),nper+1);
tsms   = min(param(53+10*nper),nper+1);
cap0   = param(54+10*nper:67+10*nper);

k1p    = [cap0(1); k1pn(1:nper-1)];
k1fs   = [cap0(2); k1fsn(1:nper-1)];
k1ms   = [cap0(3); k1msn(1:nper-1)];
k2p    = [cap0(4); k2pn(1:nper-1)];
k2fs   = [cap0(5); k2fsn(1:nper-1)];
k2ms   = [cap0(6); k2msn(1:nper-1)];
hfp    = [cap0(7); hfpn(1:nper-1)];
hmp    = [cap0(8); hmpn(1:nper-1)];
hfs    = [cap0(9); hfsn(1:nper-1)];
hms    = [cap0(10);hmsn(1:nper-1)];
efp    = [cap0(11);efpn(1:nper-1)];
emp    = [cap0(12);empn(1:nper-1)];
efs    = [cap0(13);efsn(1:nper-1)];
ems    = [cap0(14);emsn(1:nper-1)];

%
% Construct other endogenous variables
%

z1fp                = b1f*hfp.^kap1.*efp.^zet1.*l1fp;    % = phi1fp
z2fp                = b2f*hfp.^kap2.*efp.^zet2.*l2fp;    % = phi2fp
z1mp                = b1m*hmp.^kap1.*emp.^zet1.*l1mp;    % = phi1mp
z1fs                = b1f*hfs.^kap1.*efs.^zet1.*l1fs;    % = phi1fs
z2fs                = b2f*hfs.^kap2.*efs.^zet2.*l2fs;    % = phi2fs
z1ms                = b1m*hms.^kap1.*ems.^zet1.*l1ms;    % = phi1ms
z2ms                = b2m*hms.^kap2.*ems.^zet2.*l2ms;    % = phi2ms

phi1lfp             = z1fp./l1fp;
phi1hfp             = kap1*z1fp./hfp;
phi1efp             = zet1*z1fp./efp;
phi2lfp             = z2fp./l2fp;
phi2hfp             = kap2*z2fp./hfp;
phi2efp             = zet2*z2fp./efp;

phi1lmp             = z1mp./l1mp;
phi1hmp             = kap1*z1mp./hmp;
phi1emp             = zet1*z1mp./emp;

phi1lfs             = z1fs./l1fs;
phi1hfs             = kap1*z1fs./hfs;
phi1efs             = zet1*z1fs./efs;
phi2lfs             = z2fs./l2fs;
phi2hfs             = kap2*z2fs./hfs;
phi2efs             = zet2*z2fs./efs;

phi1lms             = z1ms./l1ms;
phi1hms             = kap1*z1ms./hms;
phi1ems             = zet1*z1ms./ems;
phi2lms             = z2ms./l2ms;
phi2hms             = kap2*z2ms./hms;
phi2ems             = zet2*z2ms./ems;

zbar                = nfs*z1fs+nms*z1ms+np*z1fp+np*z1mp;
kbar                = nfs*k1fs+nms*k1ms+np*k1p;

[F1,F1z,F1k]        = prodfn(zbar,kbar,theta1,mtech);
[F2p,F2zp,F2kp]     = prodfn(z2fp,k2p, theta2,htech);
[F2fs,F2zfs,F2kfs]  = prodfn(z2fs,k2fs,theta2,htech);
[F2ms,F2zms,F2kms]  = prodfn(z2ms,k2ms,theta2,htech);
[g1p,g1fp,g1mp]     = consfn(c1fp,c1mp,gpar);
[g2p,g2fp,g2mp]     = consfn(c2fp,c2mp,gpar);
c2fs                = F2fs;
c2ms                = F2ms;
c1p                 = c1fp+c1mp;
r                   = F1k;
w                   = F1z;
xk1p                = (1+gama)*(1+gamp)*k1pn-(1-delk)*k1p;
xk1fs               = (1+gama)*(1+gamp)*k1fsn-(1-delk)*k1fs;
xk1ms               = (1+gama)*(1+gamp)*k1msn-(1-delk)*k1ms;
xk2p                = (1+gama)*(1+gamp)*k2pn-(1-delk)*k2p;
xk2fs               = (1+gama)*(1+gamp)*k2fsn-(1-delk)*k2fs;
xk2ms               = (1+gama)*(1+gamp)*k2msn-(1-delk)*k2ms;
xhfp                = (1+gama)*(1+gamp)*hfpn-(1-delh)*hfp;
xefp                = (1+gama)*(1+gamp)*efpn-(1-dele)*efp;
xhmp                = (1+gama)*(1+gamp)*hmpn-(1-delh)*hmp;
xemp                = (1+gama)*(1+gamp)*empn-(1-dele)*emp;
xhfs                = (1+gama)*(1+gamp)*hfsn-(1-delh)*hfs;
xefs                = (1+gama)*(1+gamp)*efsn-(1-dele)*efs;
xhms                = (1+gama)*(1+gamp)*hmsn-(1-delh)*hms;
xems                = (1+gama)*(1+gamp)*emsn-(1-dele)*ems;

if util==1; 
  [Ufp,Uf1p,Uf2p,Uflp] = util1(c1fp,c2fp,1-l1fp-l2fp,upar);
  [Ump,Um1p,Um2p,Umlp] = util1(c1mp,c2mp,1-l1mp,upar);
  [Ufs,Uf1s,Uf2s,Ufls] = util1(c1fs,c2fs,1-l1fs-l2fs,upar);
  [Ums,Um1s,Um2s,Umls] = util1(c1ms,c2ms,1-l1ms-l2ms,upar);
end;
if util==2;
  [Ufp,Uf1p,Uf2p,Uflp] = util2(c1fp,c2fp,1-l1fp-l2fp,upar);
  [Ump,Um1p,Um2p,Umlp] = util2(c1mp,c2mp,1-l1mp,upar);
  [Ufs,Uf1s,Uf2s,Ufls] = util2(c1fs,c2fs,1-l1fs-l2fs,upar);
  [Ums,Um1s,Um2s,Umls] = util2(c1ms,c2ms,1-l1ms-l2ms,upar);
end;
if util==3;
  [Ufp,Uf1p,Uf2p,Uflp] = util3(c1fp,c2fp,1-l1fp-l2fp,upar);
  [Ump,Um1p,Um2p,Umlp] = util3(c1mp,c2mp,1-l1mp,upar);
  [Ufs,Uf1s,Uf2s,Ufls] = util3(c1fs,c2fs,1-l1fs-l2fs,upar);
  [Ums,Um1s,Um2s,Umls] = util3(c1ms,c2ms,1-l1ms-l2ms,upar);
end;

rhs1 = beth*((1-tauk)*r+delk*tauk+(1-delk)*qk1);
rhs2 = beth*((1-delh)*qh+(1-taulfp).*w.*(phi1hfp+(phi1lfp.*phi2hfp)./phi2lfp));
rhs3 = beth*((1-dele)*qe+(1-taulfp).*w.*(phi1efp+(phi1lfp.*phi2efp)./phi2lfp));
rhs4 = beth*((1-delh)*qh+(1-taulmp).*w.*phi1hmp);
rhs5 = beth*((1-dele)*qe+(1-taulmp).*w.*phi1emp);
rhs6 = beth*((1-delh)*qh+(1-taulfs).*w.*(phi1hfs+(phi1lfs.*phi2hfs)./phi2lfs));
rhs7 = beth*((1-dele)*qe+(1-taulfs).*w.*(phi1efs+(phi1lfs.*phi2efs)./phi2lfs));
rhs8 = beth*((1-delh)*qh+(1-taulms).*w.*(phi1hms+(phi1lms.*phi2hms)./phi2lms));
rhs9 = beth*((1-dele)*qe+(1-taulms).*w.*(phi1ems+(phi1lms.*phi2ems)./phi2lms));

Uf1p = Uf1p./g1fp;
Um1p = Um1p./g1mp;
psi2 = upar(2);
psi3 = upar(3);
sigma= upar(4);
mufp = mufp.*grate.^((1-psi3)*(1-sigma)-psi2)*10^10;
mump = mump.*grate.^((1-psi3)*(1-sigma)-psi2)*10^10;
mufs = mufs.*grate.^((1-psi3)*(1-sigma)-psi2)*10^10;
mums = mums.*grate.^((1-psi3)*(1-sigma)-psi2)*10^10;
Uf2p = (Uf2p-mufp)./g2fp;
Um2p = (Um2p-mump)./g2mp;
Uf2s = Uf2s-mufs;
Um2s = Um2s-mums;


f    = [Uflp./Uf1p-w.*(1-taulfp).*phi1lfp, ...
        Umlp./Um1p-w.*(1-taulmp).*phi1lmp, ...
        Uflp./Uf2p-F2zp.*phi2lfp, ...
        Uf1p./Um1p-(1-lam)/lam*ones(nper,1), ...
        Uf2p./Um2p-(1-lam)/lam*ones(nper,1), ...
        g2p-F2p, ...
        Ufls./Uf1s-w.*(1-taulfs).*phi1lfs, ...
        Ufls./Uf2s-F2zfs.*phi2lfs, ...
        Umls./Um1s-(1-taulms).*w.*phi1lms, ...
        Umls./Um2s-F2zms.*phi2lms, ...
        (1-tauk)*r.*k1p+tauk*delk*k1p+(1-taulfp).*w.*z1fp+ ...
          (1-taulmp).*w.*z1mp+(Tdp+Tp)/np-(c1p+qk1.*xk1p+qk2.*xk2p+ ...
          qh.*xhfp+qe.*xefp+qh.*xhmp+qe.*xemp), ...
        (1-tauk)*r.*k1fs+tauk*delk*k1fs+(1-taulfs).*w.*z1fs+ ...
          (Tdfs+Tfs)/nfs-(c1fs+qk1.*xk1fs+qk2.*xk2fs+qh.*xhfs+qe.*xefs), ...
        (1-tauk)*r.*k1ms+tauk*delk*k1ms+(1-taulms).*w.*z1ms+ ...
          (Tdms+Tms)/nms-(c1ms+qk1.*xk1ms+qk2.*xk2ms+qh.*xhms+qe.*xems), ...
        nfs*(c1fs+qk1.*xk1fs+qk2.*xk2fs+qh.*xhfs+qe.*xefs)+ ...
           nms*(c1ms+qk1.*xk1ms+qk2.*xk2ms+qh.*xhms+qe.*xems)+ ...
           np*(c1p+qk1.*xk1p+qk2.*xk2p+qh.*xhfp+qe.*xefp+qh.*xhmp+ ...
               qe.*xemp)-F1*(1-g), ...
        Tdp+Tdfs+Tdms+np*(taulms-taulmp).*w.*z1mp- ...
          np*(taulfp-taulms).*w.*z1fp-nfs*(taulfs-taulms).*w.*z1fs, ...
        [Uf1p(2:nper)./Uf1p(1:nper-1)-Uf1s(2:nper)./Uf1s(1:nper-1), ...
         Uf1p(2:nper)./Uf1p(1:nper-1)-Um1s(2:nper)./Um1s(1:nper-1), ...
         qk1(1:nper-1)-rhs1(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qh(1:nper-1)-rhs2(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qe(1:nper-1)-rhs3(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qh(1:nper-1)-rhs4(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qe(1:nper-1)-rhs5(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qh(1:nper-1)-rhs6(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qe(1:nper-1)-rhs7(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qh(1:nper-1)-rhs8(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qe(1:nper-1)-rhs9(2:nper).*Uf1p(2:nper)./Uf1p(1:nper-1), ...
         qk2(1:nper-1)./qk1(1:nper-1)- ...
           (F2kp(2:nper).*Uf2p(2:nper)+(1-delk)*qk2(2:nper).* ...
           Uf1p(2:nper))./(Uf1p(2:nper).*((1-tauk)*r(2:nper)+tauk*delk+ ...
           (1-delk)*qk1(2:nper))), ...
         qk2(1:nper-1)./qk1(1:nper-1)- ...
           (F2kfs(2:nper).*Uf2s(2:nper)+(1-delk)*qk2(2:nper).* ...
           Uf1s(2:nper))./(Uf1s(2:nper).*((1-tauk)*r(2:nper)+tauk*delk+ ...
           (1-delk)*qk1(2:nper))), ...
         qk2(1:nper-1)./qk1(1:nper-1)- ...
           (F2kms(2:nper).*Um2s(2:nper)+(1-delk)*qk2(2:nper).* ...
           Um1s(2:nper))./(Um1s(2:nper).*((1-tauk)*r(2:nper)+tauk*delk+ ...
           (1-delk)*qk1(2:nper)));
         [k1pn(nper),k2pn(nper),k1fsn(nper),k2fsn(nper), ...
          k1msn(nper),k2msn(nper),hfpn(nper),efpn(nper), ...
          hmpn(nper),empn(nper),hfsn(nper),efsn(nper), ...
          hmsn(nper),emsn(nper)]- ...
         [k1pn(nper-1),k2pn(nper-1),k1fsn(nper-1),k2fsn(nper-1), ...
          k1msn(nper-1),k2msn(nper-1),hfpn(nper-1),efpn(nper-1), ...
          hmpn(nper-1),empn(nper-1),hfsn(nper-1),efsn(nper-1), ...
          hmsn(nper-1),emsn(nper-1)]], ...
         [mufp(1:tsfp-1);c2fp(tsfp:nper)-csfp(tsfp:nper)], ...
         [mump(1:tsmp-1);c2mp(tsmp:nper)-csmp(tsmp:nper)], ...
         [mufs(1:tsfs-1);c2fs(tsfs:nper)-csfs(tsfs:nper)], ...
         [mums(1:tsms-1);c2ms(tsms:nper)-csms(tsms:nper)]];
         
f    = f';
f    = f(:);

