%Intang      Start with levels for ym, yo, c, xm, xo, g, h, ho, tauc,
% _Steady    and tauh and back out some of the parameters for the
%            3-sector growth model with intangible capital described
%            in the Appendix to
%
%               Unmeasured Investment and the 1990s US Hours Boom
%               by Ellen McGrattan and Ed Prescott
                                                                                
%            Ellen McGrattan, 7-26-05
%            Revised, ERM, 8-12-07

%---------------------------------------------------------------------
% 1. Observed levels

GDPs       = 1;
cs         = .7626*GDPs;
xms        = .2377*GDPs;
xns        = .1325*GDPs;
yns        = .3379*GDPs;
hs         = .2751;
hns        = .0611;
tauc       = .0657;
tauh       = .3109;  
kTs        = 1.74/1.047*GDPs;
NIPAbw     = .464/1.047*GDPs;
NIPAptax   = .025/1.047*GDPs;

%
% 2. Fixed parameters
%

beta       = .98;
eta        = .01;
gamma      = .02;
sigma      = 1.00001;
taup       = .35;
taud       = .15;
taux       = 0;
taug       = 0;
bhat       = beta*(1+gamma)^(-sigma);
grate      = (1+eta)*(1+gamma);
isymm      = 0;   %Symmetric case = 1
q          = 1;

%
% 3. Implied parameters
%

yms        = GDPs;
ybs        = GDPs-yns;
xTs        = xms-xns;
gs         = GDPs-cs-xms;
hbs        = hs-hns;
tauk       = NIPAptax/kTs;

deltT      = xTs/kTs+ 1-grate;
rT         = ((1-bhat*(1-deltT))*(1+taux)-bhat*taup*deltT+ ...
              bhat*(1-taup)*tauk)/ (bhat*(1-taup));
phi1       = input('phi1  = ? ');
deltI      = input('deltI = ? ');
thet2      = input('thet2 = ? ');
phi2       = input('phi2  = ? ');

x0         = [.26;.8];
for i=1:20;
  thet1    = x0(1);
  chi      = x0(2); 
  rI       = q*((1-bhat*(1-deltI))*(taug+(1-chi)*(1-tauh)+ ...
             chi*(1-taup)*(1-taud)))/ (bhat*(1-taud)*(1-taup));
  if isymm==1;
    thet2  = thet1;  
    phi2   = phi1;  
  end;
  lshare1  = 1-thet1-phi1;
  lshare2  = 1-thet2-phi2;
  xIs      = (NIPAbw-lshare1*ybs)/(lshare2*q-(1-chi)*q);
  kIs      = xIs/(grate-1+deltI);
  k1Ts     = thet1*ybs/rT;
  k2Ts     = thet2*q*xIs/rT;
  res      = [kTs-k1Ts-k2Ts;
              rI-(phi1*ybs+phi2*q*xIs)/kIs];

  thet1    = x0(1)+.0001;
  if isymm==1;
    thet2  = thet1;  
    phi2   = phi1;  
  end;
  lshare1  = 1-thet1-phi1;
  lshare2  = 1-thet2-phi2;
  xIs      = (NIPAbw-lshare1*ybs)/(lshare2*q-(1-chi)*q);
  kIs      = xIs/(grate-1+deltI);
  k1Ts     = thet1*ybs/rT;
  k2Ts     = thet2*q*xIs/rT;
  dres(:,1)= ([kTs-k1Ts-k2Ts;
               rI-(phi1*ybs+phi2*q*xIs)/kIs]-res)/.0001;

  thet1    = x0(1);
  chi      = x0(2)+.0001;
  rI       = q*((1-bhat*(1-deltI))*(taug+(1-chi)*(1-tauh)+ ...
             chi*(1-taup)*(1-taud)))/ (bhat*(1-taud)*(1-taup));
  if isymm==1;
    thet2  = thet1;  
    phi2   = phi1;  
  end;
  lshare1  = 1-thet1-phi1;
  lshare2  = 1-thet2-phi2;
  xIs      = (NIPAbw-lshare1*ybs)/(lshare2*q-(1-chi)*q);
  kIs      = xIs/(grate-1+deltI);
  k1Ts     = thet1*ybs/rT;
  k2Ts     = thet2*q*xIs/rT;
  dres(:,2)= ([kTs-k1Ts-k2Ts;
               rI-(phi1*ybs+phi2*q*xIs)/kIs]-res)/.0001;
  x0       = x0-dres\res;
end;

thet1      = x0(1);
chi        = x0(2);
if isymm==1;
  thet2    = thet1;  
  phi2     = phi1;  
end;
lshare1    = 1-thet1-phi1;
lshare2    = 1-thet2-phi2;
xIs        = (NIPAbw-lshare1*ybs)/(lshare2*q-(1-chi)*q);
kIs        = xIs/(grate-1+deltI);
k1Ts       = thet1*ybs/rT;
k2Ts       = thet2*q*xIs/rT;
tem        = lshare1*ybs/(lshare2*q*xIs);
h1s        = tem/(1+tem) *(hs-hns);
psi        = (1-tauh)*lshare1*(1-hs)*ybs/((1+tauc)*cs*h1s);
h2s        = hs-h1s-hns;
A1s        = ybs/(k1Ts^thet1*kIs^phi1*h1s^lshare1);
A2s        = xIs/(k2Ts^thet2*kIs^phi2*h2s^lshare2);
val        = (1-taud)*(1+taux)*kTs+(taug+chi*(1-taud)*(1-taup)+ ...
             (1-chi)*(1-tauh))*q*kIs;

  
% 4. Write out results
%
fid        = fopen('intang_params.m','w');
fprintf(fid,'beta       = %11.8f;\n',beta);
fprintf(fid,'deltI      = %11.8f;\n',deltI);
fprintf(fid,'deltT      = %11.8f;\n',deltT);
fprintf(fid,'eta        = %11.8f;\n',eta);
fprintf(fid,'gamma      = %11.8f;\n',gamma);
fprintf(fid,'phi1       = %11.8f;\n',phi1);
fprintf(fid,'phi2       = %11.8f;\n',phi2);
fprintf(fid,'psi        = %11.8f;\n',psi);
fprintf(fid,'sigma      = %11.8f;\n',sigma);
fprintf(fid,'thet1      = %11.8f;\n',thet1);
fprintf(fid,'thet2      = %11.8f;\n',thet2);
fprintf(fid,'z1s        = %11.8f;\n',A1s^(1/lshare1));
fprintf(fid,'z2s        = %11.8f;\n',A2s^(1/lshare2));
fprintf(fid,'gs         = %11.8f;\n',gs);
fprintf(fid,'taucs      = %11.8f;\n',tauc);
fprintf(fid,'tauhs      = %11.8f;\n',tauh);
fprintf(fid,'taups      = %11.8f;\n',taup);
fprintf(fid,'tauks      = %11.8f;\n',tauk);
fprintf(fid,'tauds      = %11.8f;\n',taud);
fprintf(fid,'tauxs      = %11.8f;\n',taux);
fprintf(fid,'taugs      = %11.8f;\n',taug);
fprintf(fid,'chis       = %11.8f;\n',chi);
fprintf(fid,'kTs        = %11.8f;\n',kTs);
fprintf(fid,'k1Ts       = %11.8f;\n',k1Ts);
fprintf(fid,'k2Ts       = %11.8f;\n',k2Ts);
fprintf(fid,'kIs        = %11.8f;\n',kIs);
fprintf(fid,'xms        = %11.8f;\n',xms);
fprintf(fid,'xns        = %11.8f;\n',xns);
fprintf(fid,'xTs        = %11.8f;\n',xTs);
fprintf(fid,'xIs        = %11.8f;\n',xIs);
fprintf(fid,'hs         = %11.8f;\n',hs);
fprintf(fid,'h1s        = %11.8f;\n',h1s);
fprintf(fid,'h2s        = %11.8f;\n',h2s);
fprintf(fid,'hns        = %11.8f;\n',hns);
fprintf(fid,'ybs        = %11.8f;\n',ybs);
fprintf(fid,'yns        = %11.8f;\n',yns);
fprintf(fid,'gdps       = %11.8f;\n',ybs+yns);
fclose(fid);

