function [res,stats,Xtran,partran]=steady(X,param)
%COMPSTEADY  computes steady state for McGrattan-Waddle paper.
%
%            The model assumes multinationals with technology (ome) maximize
%            dividends subject to contraints on capital:
%
%               max (1-taud) sum_t p_t D_t(ome)
%
%               st  D(ome) =sum_i {(1-taup_i)[P_i(ome)(Y_i(ome)-delT K_i(ome)
%                                - X_Ii(ome)-chi_i(ome) X_M(ome)]
%                                - W_iL_i(ome)-P_i(ome)(K'_Ti(ome)-K_Ti(ome))}
%
%                   K'_Ti(ome) = (1-delT) K_Ti(ome)+X_Ti(ome)+adj costs
%                   K'_Ii(ome) = (1-delI) K_Ii(ome)+X_Ii(ome)+adj costs
%                   M'(ome)    = (1-delM) M(ome)   +X_M(ome) +adj costs
%                   Y_i(ome)   = A_i sig_i(ome) (N_i M(ome))^phi Z_i(ome)
%                   Z_i(ome)   = K_Ti(ome)^alpT K_Ii(ome)^alpI L_i(ome)^alpL
%
%
%            Households in country i maximixe
%
%               max sum_t beta^t (log C_it + psi log (1-L_it)) N_it
%              
%               st  sum_t p_t {(1+tauc_i) P_it [t] C_it + sum_ome V_t(ome)*
%                                     (S_it+1(ome)-S_it(ome) + B_it+1-B_it}
%                      <= sum_t p_t {(1-taul_i) W_it L_it +(1-taud) 
%                                     sum_ome S_it(ome) D_t(ome)
%                                     +rb_t B_it + Kap_it
%
%                   C_i      = CES(C_i(ome))
%                   C_i(ome) = CES(C^fdi_i(ome),C^trade_i(ome)) 
%                              if multinational with ome not from i
%    
%

%            Ellen McGrattan, 9-7-06
%            Revised, 7-30-18


%
% Fixed parameters:
%

I          = param(1);
alphI      = param(2);
alphT      = param(3);
beta       = param(4);
deltI      = param(5);
deltM      = param(6);
deltT      = param(7);
gamN       = param(8); 
gamY       = param(9);
phi        = param(10);
psi        = param(11);
rho        = param(12);
varrho     = param(13);
A          = param(14:13+I);
N          = param(14+I:13+2*I);
b          = param(14+2*I:13+3*I);
tauc       = param(14+3*I:13+4*I);
taul       = param(14+4*I:13+5*I);
taud       = param(14+5*I:13+6*I);
taup       = param(14+6*I:13+7*I);
Tran       = param(14+7*I:13+8*I);
lnb        = param(14+8*I:13+9*I);
xnby       = param(14+9*I:13+10*I);
ynby       = param(14+10*I:13+11*I);
chiji      = zeros(I);
shrji      = zeros(I);
icej2i     = zeros(I);
sigji      = zeros(I);
cjiwgt     = zeros(I);
cjifwgt    = zeros(I);
cjitwgt    = zeros(I);
chiji(:)   = param(14+11*I:13+11*I+I*I);
shrji(:)   = param(14+11*I+I*I:13+11*I+2*I*I);
icej2i(:)  = abs(param(14+11*I+2*I*I:13+11*I+3*I*I));
sigji(:)   = param(14+11*I+3*I*I:13+11*I+4*I*I);
cjiwgt(:)  = param(14+11*I+4*I*I:13+11*I+5*I*I);
cjifwgt(:) = param(14+11*I+5*I*I:13+11*I+6*I*I);
cjitwgt(:) = param(14+11*I+6*I*I:13+11*I+7*I*I);
adjpar     = 0;
if length(param)>13+11*I+7*I*I;
  adjpar   = param(14+11*I+7*I*I);
end;

alphL      = 1-alphT-alphI;
gamy       = (1+gamY)/(1+gamN)-1;
betg       = beta/(1+gamy);
irate      = (1+gamy)/beta-1;
phi1       = 1-phi;
rho1       = 1-rho;
varrho1    = 1-varrho;
o          = ones(I,1);

%
% Unknowns being computed
%
ci         = X(1:I);
li         = X(I+1:2*I);
Mj         = X(2*I+1:3*I);
pji        = zeros(I);
pji(:)     = [1;X(3*I+1:I*I+3*I-1)];

%
% Intermediate variables
%

%   * prices
%
pjifdi     = pji;
pjitr      = (diag(pji)*o').*(1+icej2i);
pjiast     = (cjifwgt.^varrho.*pjifdi.^varrho1+cjitwgt.^varrho.*pjitr.^ ...
              varrho1).^(1/varrho1);
tem        = cjiwgt.^rho.*pjiast.^rho1;
tem        = tem-diag(diag(tem));
pi         = (diag(cjiwgt).^rho.*diag(pji).^rho1 + sum(tem)').^(1/rho1);

%   * consumptions
%

cjiast     = (pjiast./(cjiwgt.*(ones(I,1)*pi'))).^(-rho).*(ones(I,1)*ci');
cjifdi     = (pjifdi./(cjifwgt.*pjiast)).^(-varrho).*cjiast;
cjitr      = (pjitr./(cjitwgt.*pjiast)).^(-varrho).*cjiast;
cii        = (diag(pji)./(diag(cjiwgt).*pi)).^(-rho).*ci;

%   * utilities
%
uc         = 1./ci;
ul         = psi*o./(1-li-lnb);


%   * interest rates
%
rb         = irate;
rT         = rb*o./(1-taup)+deltT;
rI         = rb*o+deltI;

%   * wage rates
%
wi         = (1+tauc).*pi.*ul./((1-taul).*uc);

%   * productivities
%
prji       = (o*wi')./(phi1*alphL*pji);

%   * capital output ratios
% 
%
kTyji      = phi1*alphT*(o*(o./rT)');
kIyji      = phi1*alphI*(o*(o./rI)');

%   * outputs
%
Aji        = sigji.*(o*A');
yji        = (Aji.*kTyji.^(phi1*alphT).*(kIyji).^(phi1*alphI).* ...
             prji.^(-phi1*alphL)).^(1/phi).*(Mj*o');

%
%   * labor inputs
%
lji        = yji./prji;

%   * investments
%
kTji       = yji.*kTyji;
kIji       = yji.*kIyji;
xTji       = (gamY+deltT)*kTji;
xIji       = (gamY+deltI)*kIji;
xMj        = (gamY+deltM)*Mj;

%   * dividends
%
dji        = (o*(1-taup')).*((o*N').*(pji.*(yji-deltT*kTji-xIji)- ...
             (o*wi').*lji)-chiji.*pji.*(xMj*o'))-gamY*(o*N').*pji.*kTji;
dj         = sum(dji')';         

%   * transfers
%
gdpi       = sum(pji.*(yji-xIji-chiji.*(xMj*(o./N)')))'./(1-ynby);
ynbi       = ynby.*gdpi;
xnbi       = xnby.*gdpi;
kapi       = tauc.*pi.*ci+taul.*wi.*li+taud.*sum(shrji.*(dj*o'))'./N+ ...
             taup.*(sum(pji.*(yji-deltT*kTji-xIji))'-wi.*li- ...
             sum(chiji.*pji.*(xMj*o'))'./N)+ynbi-xnbi-Tran./N;
             %sum(chiji.*pji.*(xMj*o'))'./N)+diag(pji).*(ynbi-xnbi)-Tran./N;

%   * returns to technology capital
%
rMj        = sum(((o*(1-taup')).*pji.*(phi*yji.*((o./Mj)*N')+(1-deltM)* ...
             chiji))')'./sum(((o*(1-taup')).*chiji.*pji)')'-1;

%
% 9 residual equations 
%
ctrade     = sum(((o*N').*(1+icej2i).*cjitr./(N*o'))')'- ...
             diag((1+icej2i).*cjitr);
rcown      = diag(yji-xTji-xIji)-cii-ctrade-xMj./N+(ynbi-xnbi)./diag(pji);
rc         = yji-cjifdi-xTji-xIji;
for i=1:I;
  rc(i,i)  = rcown(i);
end;
rc         = rc(:);
res        = [(rb-gamY)*b+(1-taul).*wi.*li+(1-taud).* ...
                  sum(shrji.*(dj*o'))'./N+kapi-(1+tauc).*pi.*ci;
              sum(lji)'-li;
              rMj-rb;
              rc(2:I*I)];
Nji        = o*N';
GDPji      = pji.*(Nji.*(yji-xIji)-chiji.*(xMj*o')) + diag(N.*ynbi);
GDPi       = sum(GDPji)';
NXi        = GDPi-N.*(sum(pji.*xTji)'+pi.*ci+xnbi);
EQPi       = sum((shrji.*(dj*o'))')'-diag(shrji.*(dj*o'));
EQRi       = sum(shrji.*(dj*o'))'-diag(shrji.*(dj*o'));
DIArei     = gamY*(sum((pji.*kTji.*Nji)')'-diag(pji.*kTji.*Nji));
FDIrei     = gamY*(sum(pji.*kTji.*Nji)'-diag(pji.*kTji.*Nji));
DIAdivi    = sum(dji')'-diag(dji);
FDIdivi    = sum(dji)'-diag(dji);
Bi         = N.*b;
DBRi       = max(rb*Bi,0);
DBPi       = max(-rb*Bi,0);
EXi        = sum((pjitr.*(o*N').*cjitr)')'-diag(pjitr.*(o*N').*cjitr);
IMi        = sum(pjitr.*(o*N').*cjitr)'-diag(pjitr.*(o*N').*cjitr);

NFRi       = DIAdivi+DIArei+EQRi+DBRi;
NFPi       = FDIdivi+FDIrei+EQPi+DBPi;
NXi        = EXi-IMi;
CAi        = NXi+NFRi-NFPi-Tran;
FAout      = gamY*(sum((kTji.*(o*N'))')'-diag(kTji.*(o*N')));
FAin       = gamY*(sum(kTji.*(o*N'))'-diag(kTji.*(o*N')));
FAi        = FAout-FAin;

%
% Stats for setting A, sig, ice:
%
IMji       = pjitr.*(o*N').*cjitr-diag(diag(pjitr.*(o*N').*cjitr));
FDIIji     = gamY*(Nji.*pji.*kTji-diag(diag(Nji.*pji.*kTji)));
RGDPPC     = (GDPi./N)/(GDPi(1)/N(1))*100;
FDIIGDPji  = FDIIji./(o*GDPi') *100;
IMGDPji    = IMji./(o*GDPi') *100;

stats      = [RGDPPC;
              FDIIGDPji(:);
              IMGDPji(:)];
if (nargout>2);
  T        = 100;
  partran  = ones(T,1)*[A;N;taup;shrji(:);sigji(:);icej2i(:)]';
  partran  = [I;T;alphI;alphT;beta;deltI;deltM;deltT;gamN;gamY;
              phi;psi;rho;varrho;kIji(1,:)'./kTji(1,:)';chiji(:);cjiwgt(:);
              cjifwgt(:);cjitwgt(:);kTji(:);Mj;b;tauc;taul;taud;Tran;
              lnb;xnby;ynby;partran(:);adjpar];
  tem      = pji(:);
  Xtran    = ones(T,1)*[ci;li;Mj;kTji(:);tem(2:I*I);b(2:I)]';
  if (adjpar>0);
    Xtran  = ones(T,1)*[ci;li;xMj;xTji(:);xIji(:);tem(2:I*I);b(2:I)]';
  end;
end;

