function [res,V]  = fdifoc(X0,param)
%Equations for 
%
%   max sum beta^t (log c[i,t] + psi log (1-l[i,t])) N[i,t]
%    
%   subj. to 
%
%     sum_t p[t]{(1+tauc) C[i,t] + sum_j V[j,t] (S[j,i,t+1]-S[j,i,t])
%         + B[i,t+1]-B[i,t]} <= sum_t p[t]{(1-taul) W[i,t] L[i,t]
%                                       +(1-taud[t]) sum_j S[j,i,t] D[j,t]
%                                       +rb[t] B[i,t] + Kap[i,t]
%
%     KT[j,i,t+1]   = (1-deltM) KT[j,i,t] + XT[j,i,t]
%     KI[j,i,t+1]   = (1-deltI) KI[j,i,t] + XI[j,i,t]
%     M[j,t+1]      = (1-deltM) M[j,t]    + XM[j,t]
%
%
%   i   = country index (u=US, r=ROW)
%   j   = multinational index (d=Dell or domestic, f=Fujitsu or foreign)
%
%   Technologies: Y   = A Z^(1-phi) M^phi
%                 Z   = KT^alphT KI^alphI L^(1-alphT-alphI)
%

%   Ellen McGrattan, 9-7-06
%   Revised, 4-9-08

%
% Fixed parameters:
%

beta    = param(1);
chidu   = param(2);
chifu   = param(3);
deltT   = param(4);
deltM   = param(5);
deltI   = param(6);
gamN    = param(7); 
gamY    = param(8);
phi     = param(9);
psi     = param(10);
tauc    = param(11);
taul    = param(12);
alphT   = param(13);
alphI   = param(14);
b0      = param(15);
sdu0    = param(16);
sfu0    = param(17);
T       = param(18);
bn      = param(19:18+T);
taud    = param(19+T:18+2*T);
taupu   = param(19+2*T:18+3*T);
taupr   = param(19+3*T:18+4*T);
Tran    = param(19+4*T:18+5*T);
lnbu    = param(19+5*T:18+6*T);
lnbr    = param(19+6*T:18+7*T);
xnbuy   = param(19+7*T:18+8*T);
ynbuy   = param(19+8*T:18+9*T);
xnbry   = param(19+9*T:18+10*T);
ynbry   = param(19+10*T:18+11*T);
Nr      = param(19+11*T:18+12*T);
Nu      = param(19+12*T:18+13*T);
sigr    = param(19+13*T:18+14*T);
sigu    = param(19+14*T:18+15*T);
sfun    = param(19+15*T:18+16*T);
rsize   = param(19+16*T:18+17*T);
Adu0    = param(19+17*T);
kTdu0   = param(20+17*T);
kTdr0   = param(21+17*T);
kTfu0   = param(22+17*T);
kTfr0   = param(23+17*T);
Md0     = param(24+17*T);
Mf0     = param(25+17*T);
kIdivkT = param(26+17*T);
kIdu0   = kIdivkT*kTdu0;
kIdr0   = kIdivkT*kTdr0;
kIfu0   = kIdivkT*kTfu0;
kIfr0   = kIdivkT*kTfr0;
Adu     = Adu0*ones(T,1);
Afr     = Adu.*(rsize.*Nu./Nr).^(1-(alphT+alphI)*(1-phi));
Adr     = sigr.*Afr;
Afu     = sigu.*Adu;
sfu     = [sfu0;sfun(1:T-1)];
taupun  = [taupu(2:T);taupu(T)];
tauprn  = [taupr(2:T);taupr(T)];
Nrn     = [Nr(2:T);Nr(T)];
Nun     = [Nu(2:T);Nu(T)];
chidr   = 1-chidu;
chifr   = 1-chifu;
alphL   = 1-alphT-alphI;
gamy    = (1+gamY)/(1+gamN)-1;
betg    = beta/(1+gamy);
irate   = (1+gamy)/beta-1;
phi1    = 1-phi;

%
% Unknowns being computed
%
X       = zeros(T,11);
X(:)    = X0;
cu      = X(:,1);
cr      = X(:,2);
lu      = X(:,3);
lr      = X(:,4);
sdun    = X(:,5);
kTdun   = X(:,6);
kTdrn   = X(:,7);
kTfun   = X(:,8);
kTfrn   = X(:,9);
Mdn     = X(:,10);
Mfn     = X(:,11);

%
% Intermediate variables
%

%   * utilities
%
o       = ones(T,1);
ucu     = o./cu;
ucr     = o./cr;
ulu     = psi*o./(1-lu-lnbu);
ulr     = psi*o./(1-lr-lnbr);
ucun    = [ucu(2:T);ucu(T)];
ucrn    = [ucr(2:T);ucr(T)];

%   * interest rates
%
rbn     = (ucu./ucun/betg)-1;
rb      = [irate;rbn(1:T-1)];
rTun    = rbn./(1-taupun)+deltT;
rTrn    = rbn./(1-tauprn)+deltT;
rIun    = (1-taupu).*(1+rbn)./(1-taupun)+deltI-1;
rIrn    = (1-taupr).*(1+rbn)./(1-tauprn)+deltI-1;

%   * intangible capitals
%
uratio  = alphI*rTun./(alphT*rIun);
rratio  = alphI*rTrn./(alphT*rIrn);
kIdun   = uratio.*kTdun;
kIdrn   = rratio.*kTdrn;
kIfun   = uratio.*kTfun;
kIfrn   = rratio.*kTfrn;

%   * this period stocks
%
kTdu    = [kTdu0;kTdun(1:T-1)];
kTdr    = [kTdr0;kTdrn(1:T-1)];
kTfu    = [kTfu0;kTfun(1:T-1)];
kTfr    = [kTfr0;kTfrn(1:T-1)];
kIdu    = [kIdu0;kIdun(1:T-1)];
kIdr    = [kIdr0;kIdrn(1:T-1)];
kIfu    = [kIfu0;kIfun(1:T-1)];
kIfr    = [kIfr0;kIfrn(1:T-1)];
Md      = [Md0;Mdn(1:T-1)];
Mf      = [Mf0;Mfn(1:T-1)];
sdu     = [sdu0;sdun(1:T-1)];
sfu     = [sfu0;sfun(1:T-1)];
b       = [b0;bn(1:T-1)];
kIdun(T)= kIdu(T);
kIdrn(T)= kIdr(T);
kIfun(T)= kIfu(T);
kIfrn(T)= kIfr(T);

%   * wage rates
%
wu      = (1+tauc)/(1-taul)*ulu./ucu;
wr      = (1+tauc)/(1-taul)*ulr./ucr;

%   * productivities
%
prdu    = wu/(phi1*alphL);
prdr    = wr/(phi1*alphL);

%   * outputs
%
ydu     = (Adu.*kTdu.^(phi1*alphT).*kIdu.^(phi1*alphI).* ...
           prdu.^(-phi1*alphL).*Md.^phi).^(1/(1-phi1*alphL));
ydr     = (Adr.*kTdr.^(phi1*alphT).*kIdr.^(phi1*alphI).* ...
           prdr.^(-phi1*alphL).*Md.^phi).^(1/(1-phi1*alphL));
yfu     = (Afu.*kTfu.^(phi1*alphT).*kIfu.^(phi1*alphI).* ...
           prdu.^(-phi1*alphL).*Mf.^phi).^(1/(1-phi1*alphL));
yfr     = (Afr.*kTfr.^(phi1*alphT).*kIfr.^(phi1*alphI).* ...
           prdr.^(-phi1*alphL).*Mf.^phi).^(1/(1-phi1*alphL));
ydun    = [ydu(2:T);ydu(T)];
ydrn    = [ydr(2:T);ydr(T)];
yfun    = [yfu(2:T);yfu(T)];
yfrn    = [yfr(2:T);yfr(T)];

%
%   * labor inputs
%
ldu     = ydu./prdu;
ldr     = ydr./prdr;
lfu     = yfu./prdu;
lfr     = yfr./prdr;

%   * investments
%

xTdu    = (1+gamY)*Nun.*kTdun./Nu-(1-deltT)*kTdu;
xTdr    = (1+gamY)*Nrn.*kTdrn./Nr-(1-deltT)*kTdr;
xTfu    = (1+gamY)*Nun.*kTfun./Nu-(1-deltT)*kTfu;
xTfr    = (1+gamY)*Nrn.*kTfrn./Nr-(1-deltT)*kTfr;
xIdu    = (1+gamY)*Nun.*kIdun./Nu-(1-deltI)*kIdu;
xIdr    = (1+gamY)*Nrn.*kIdrn./Nr-(1-deltI)*kIdr;
xIfu    = (1+gamY)*Nun.*kIfun./Nu-(1-deltI)*kIfu;
xIfr    = (1+gamY)*Nrn.*kIfrn./Nr-(1-deltI)*kIfr;
xMd     = (1+gamY)*Mdn -(1-deltM)*Md;
xMf     = (1+gamY)*Mfn -(1-deltM)*Mf;

%   * other shares
%
sdrn    = 1-sdun;
sfrn    = 1-sfun;
sdr     = 1-sdu;
sfr     = 1-sfu;

%   * dividends
%
dd      = (1-taupu).*(Nu.*(ydu-wu.*ldu-deltT*kTdu-xIdu)-chidu*xMd)- ...
           ((1+gamY)*Nun.*kTdun-Nu.*kTdu)+ ...
          (1-taupr).*(Nr.*(ydr-wr.*ldr-deltT*kTdr-xIdr)-chidr*xMd)- ...
           ((1+gamY)*Nrn.*kTdrn-Nr.*kTdr);
df      = (1-taupu).*(Nu.*(yfu-wu.*lfu-deltT*kTfu-xIfu)-chifu*xMf)- ...
           ((1+gamY)*Nun.*kTfun-Nu.*kTfu)+ ...
          (1-taupr).*(Nr.*(yfr-wr.*lfr-deltT*kTfr-xIfr)-chifr*xMf)- ...
           ((1+gamY)*Nrn.*kTfrn-Nr.*kTfr);

%   * values (assuming interiority)
%
vd      = (1-taud).*(Nu.*(kTdun+(1-taupu).*kIdun)+ ...
                     Nr.*(kTdrn+(1-taupr).*kIdrn)+ ...
                     ((1-taupu)*chidu+(1-taupr)*chidr).*Mdn);
vf      = (1-taud).*(Nu.*(kTfun+(1-taupu).*kIfun)+ ...
                     Nr.*(kTfrn+(1-taupr).*kIfrn)+ ...
                     ((1-taupu)*chifu+(1-taupr)*chifr).*Mfn);

%   * transfers
%

gdpu    = (ydu+yfu-xIdu-xIfu-chidu*xMd./Nu-chifu*xMf./Nu)./(1-ynbuy);
gdpr    = (ydr+yfr-xIdr-xIfr-chidr*xMd./Nr-chifr*xMf./Nr)./(1-ynbry);
ynbu    = ynbuy.*gdpu;
xnbu    = xnbuy.*gdpu;
ynbr    = ynbry.*gdpr;
xnbr    = xnbry.*gdpr;
kapu    = tauc*cu+taul*wu.*lu+taud.*(sdu.*dd+sfu.*df)./Nu+ ...
          taupu.*(ydu+yfu-wu.*lu-deltT*(kTdu+kTfu)-xIdu-xIfu- ...
          (chidu*xMd+chifu*xMf)./Nu)+ynbu-xnbu-Tran./Nu;
kapr    = tauc*cr+taul*wr.*lr+taud.*(sdr.*dd+sfr.*df)./Nr+ ...
          taupr.*(ydr+yfr-wr.*lr-deltT*(kTdr+kTfr)-xIdr-xIfr- ...
          (chidr*xMd+chifr*xMf)./Nr)+ynbr-xnbr+Tran./Nr;


%   * returns to capital
%
rkTdun  = phi1*alphT*ydun./kTdun;
rkTdrn  = phi1*alphT*ydrn./kTdrn;
rkTfun  = phi1*alphT*yfun./kTfun;
rkTfrn  = phi1*alphT*yfrn./kTfrn;
rMdn    = ((1-taupun).*(phi*Nun.*ydun./Mdn+chidu*(1-deltM))+ ...
           (1-tauprn).*(phi*Nrn.*ydrn./Mdn+chidr*(1-deltM)))./ ...
          ((1-taupu)*chidu+(1-taupr)*chidr)-1;
rMfn    = ((1-taupun).*(phi*Nun.*yfun./Mfn+chifu*(1-deltM))+ ...
           (1-tauprn).*(phi*Nrn.*yfrn./Mfn+chifr*(1-deltM)))./ ...
          ((1-taupu)*chifu+(1-taupr)*chifr)-1;


%
% 11 residual equations 
%
res     = [(1+rb).*b-(1+gamY)*bn.*Nun./Nu+(1-taul)*wu.*lu+ ...
             (1-taud).*(sdu.*dd+sfu.*df)./Nu+ ...
             kapu-(1+tauc)*cu-(vd.*(sdun-sdu)+vf.*(sfun-sfu))./Nu;
          (-(1+rb).*b+(1+gamY)*bn.*Nun./Nu).*Nu./Nr+(1-taul)*wr.*lr+ ...
             (1-taud).*(sdr.*dd+sfr.*df)./Nr+ ...
             kapr-(1+tauc)*cr-(vd.*(sdrn-sdr)+vf.*(sfrn-sfr))./Nr;
           ldu+lfu-lu;
           ldr+lfr-lr;
           [ucu(1:T-1)./ucun(1:T-1)-ucr(1:T-1)./ucrn(1:T-1);sdun(T)-sdu(T)];  
           [rkTdun(1:T-1)-rTun(1:T-1);kTdun(T)-kTdu(T)];
           [rkTdrn(1:T-1)-rTrn(1:T-1);kTdrn(T)-kTdr(T)];
           [rkTfun(1:T-1)-rTun(1:T-1);kTfun(T)-kTfu(T)];
           [rkTfrn(1:T-1)-rTrn(1:T-1);kTfrn(T)-kTfr(T)];
           [rMdn(1:T-1)-rbn(1:T-1);Mdn(T)-Md(T)];
           [rMfn(1:T-1)-rbn(1:T-1);Mfn(T)-Mf(T)]];

%
% Compute series for the constraints
%

Cu      = Nu.*cu;
Cr      = Nr.*cr;
Yu      = Nu.*(ydu+yfu);
Yr      = Nr.*(ydr+yfr);
XTu     = Nu.*(xTdu+xTfu);
XTr     = Nr.*(xTdr+xTfr);
KTu     = Nu.*(kTdu+kTfu);
KTr     = Nr.*(kTdr+kTfr);
KTun    = Nun.*(kTdun+kTfun);
KTrn    = Nrn.*(kTdrn+kTfrn);
XIu     = Nu.*(xIdu+xIfu);
XIr     = Nr.*(xIdr+xIfr);
KIu     = Nu.*(kIdu+kIfu);
KIr     = Nr.*(kIdr+kIfr);
KIun    = Nun.*(kIdun+kIfun);
KIrn    = Nrn.*(kIdrn+kIfrn);
XMu     = chidu*xMd+chifu*xMf;
XMr     = chidr*xMd+chifr*xMf;
Mu      = chidu*Md+chifu*Mf;
Mr      = chidr*Md+chifr*Mf;
Mun     = chidu*Mdn+chifu*Mfn;
Mrn     = chidr*Mdn+chifr*Mfn;
GDPu    = Yu-XIu-chidu*xMd-chifu*xMf+Nu.*ynbu;
GDPr    = Yr-XIr-chidr*xMd-chifr*xMf+Nr.*ynbr;
NXu     = GDPu-XTu-Cu-Nu.*xnbu;
NXr     = GDPr-XTr-Cr-Nr.*xnbr;
Lu      = Nu.*lu;
Lr      = Nr.*lr;
Taxu    = taupu.*(Yu-wu.*Lu-XIu-XMu-deltT*KTu);
Taxr    = taupr.*(Yr-wr.*Lr-XIr-XMr-deltT*KTr);
Divu    = Yu-wu.*Lu-XIu-XMu-XTu-Taxu;
Divr    = Yr-wr.*Lr-XIr-XMr-XTr-Taxr;
REu     = (1+gamY)*KTun-KTu;
REr     = (1+gamY)*KTrn-KTr;
Prfu    = Divu+Taxu+REu;
Prfr    = Divr+Taxr+REr;
Depu    = deltT*KTu;
Depr    = deltT*KTr;
DIAdivu = Nr.*(ydr-wr.*ldr-xTdr-xIdr)-chidr*xMd- ...
          taupr.*(Nr.*(ydr-wr.*ldr-deltT*kTdr-xIdr)-chidr*xMd);
DIAreu  = ((1+gamY)*Nrn.*kTdrn-Nr.*kTdr);
EQRu    = sfu.*df;
FDIdivu = Nu.*(yfu-wu.*lfu-xTfu-xIfu)-chifu*xMf- ...
          taupu.*(Nu.*(yfu-wu.*lfu-deltT*kTfu-xIfu)-chifu*xMf);
FDIreu  = ((1+gamY)*Nun.*kTfun-Nu.*kTfu);
EQPu    = sdr.*dd;
DBRu    = max(Nu.*rb.*b,0);
DBPu    = max(-Nu.*rb.*b,0);
NFRu    = DIAdivu+DIAreu+EQRu+DBRu;
NFPu    = FDIdivu+FDIreu+EQPu+DBPu;
GNPu    = GDPu+NFRu-NFPu;
CAmod   = [(DIAdivu+DIAreu)./GNPu*100, ...
           (FDIdivu+FDIreu)./GNPu*100, ...
            NXu./GNPu*100, ...
           (EQRu+DBRu-EQPu-DBPu)./GNPu*100];
V       = [CAmod(:,3)
           GDPu(1)/(GDPu(1)+GDPr(1));
           xTdu(2)*xTdu(2)/(xTdu(1)*xTdu(3));
           xTdr(2)*xTdr(2)/(xTdr(1)*xTdr(3));
           xTfu(2)*xTfu(2)/(xTfu(1)*xTfu(3));
           xTfr(2)*xTfr(2)/(xTfr(1)*xTfr(3));
           xMd(2)*xMd(2)/(xMd(1)*xMd(3));
           xMf(2)*xMf(2)/(xMf(1)*xMf(3));
           GDPu(1)];

