function [res,V] = stdfoc(x,param);
% STDFOC  
%
%   max E sum  beta^t lam {log(Cu[t]/Nu[t])+psi*log(1-Hu[t]/Nu[t])}Nu[t] 
%          t     + (1-lam){log(Cr[t]/Nr[t])+psi*log(1-Hr[t]/Nr[t])}Nr[t]
% 
%   s.t.  Cu[t]+Cr[t]+Xu[t]+Xr[t] = Au[t] F(Ku[t],Hu[t])+ Ar[t] F(Kr[t],Hr[t])
%         Kr[t+1] = (1-delta)Kr[t]+Xr[t]
%         Ku[t+1] = (1-delta)Ku[t]+Xu[t]
%         Ku[0], Kr[0] given
%

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

%
% Parameters
%

alpha = param(1);
beta  = param(2);
delta = param(3);
lam   = param(4);
psi   = param(5);
Ku0   = param(6);
Kr0   = 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,4);
X(:)  = x;
Kun   = X(:,1);
Krn   = X(:,2);
Hu    = X(:,3);
Hr    = X(:,4);

%
% 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)];
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;
cun   = Cun./Nun;
cr    = Cr./Nr;
crn   = Crn./Nrn;
mu    = cu./cun;
mr    = cr./crn;
res   = [Yu+Yr-Xu-Xr-Cu-Cr;
         lam*Nu.*Cr-(1-lam)*Cu.*Nr;
         1-beta*cu./cun.*(alpha*Yun./Kun+1-delta);
         1-beta*cr./crn.*(alpha*Yrn./Krn+1-delta)];
