function res = decstdfoc(x,param);
% DECSTDFOC
%
%   max E sum  beta^t {log(Ci[t]/Ni[t])+psi*log(1-Hi[t]/Ni[t])}Ni[t] 
%          t     
% 
%   s.t.  Ci[t]+Xi[t]+Bi[t+1] = Ai[t] F(Ki[t],Hi[t]) + (1+rb[t])Bi[t]
%         Ki[t+1] = (1-delta)Ki[t]+Xi[t]
%
%   Market clearing: 
%         sum_i Bi[t]=0
%         Cu[t]+Cr[t]+Xu[t]+Xr[t]=Au[t]F(Ku[t],Hu[t])+Ar[t]F(Kr[t],Hr[t])
%

%   Ellen McGrattan, 10-30-07
%   Revised, ERM, 10-30-07

%
% Parameters
%

alpha = param(1);
beta  = param(2);
delta = param(3);
psi   = param(4);
Ku0   = param(5);
Kr0   = param(6);
B0    = param(7);
T     = param(8);
Au    = param(9:T+8);
Ar    = param(T+9:2*T+8);
Nu    = param(2*T+9:3*T+8);
Nr    = param(3*T+9:4*T+8);

%
% Variables
%

X     = zeros(T,5);
X(:)  = x;
Kun   = X(:,1);
Krn   = X(:,2);
Hu    = X(:,3);
Hr    = X(:,4);
Bn    = X(:,5);

%
% Intermediate variables and residuals
%

Nun   = [Nu(2:T);Nu(T)];
Nrn   = [Nr(2:T);Nr(T)];
Ku    = [Ku0;Kun(1:T-1)];
Kr    = [Kr0;Krn(1:T-1)];
B     = [B0;Bn(1:T-1)];
Yu    = Au.*Ku.^alpha.*Hu.^(1-alpha);
Yr    = Ar.*Kr.^alpha.*Hr.^(1-alpha);
Cu    = (1-alpha)*Yu.*(Nu-Hu)./(psi*Hu);
Cr    = (1-alpha)*Yr.*(Nr-Hr)./(psi*Hr);
Cun   = [Cu(2:T);Cu(T)];
Crn   = [Cr(2:T);Cr(T)];
Yun   = [Yu(2:T);Yu(T)];
Yrn   = [Yr(2:T);Yr(T)];
Xu    = Kun-(1-delta)*Ku;
Xr    = Krn-(1-delta)*Kr;
cu    = Cu./Nu;
cr    = Cr./Nr;
rb    = [1/beta-1;cu(2:T)./cu(1:T-1)/beta-1];
res   = [Cu+Xu+Bn-(1+rb).*B-Yu;
         Cr+Xr-Bn+(1+rb).*B-Yr;
         [cu(1:T-1)./cu(2:T)-cr(1:T-1)./cr(2:T);Bn(T)-B(T)];
         [1-beta*cu(1:T-1)./cu(2:T).*(alpha*Yu(2:T)./Ku(2:T)+1-delta);
            Kun(T)-Ku(T)];
         [1-beta*cr(1:T-1)./cr(2:T).*(alpha*Yr(2:T)./Kr(2:T)+1-delta);
            Krn(T)-Kr(T)]];

