%STOKEY   This file calculates steady state growth rates for a taxed
%         economy and an untaxed economy as in Stokey and Rebelo (1995).  
%         The model is the following two-sector endogenous growth model:
%
%             max sum beta^t (c[t](1-l[t])^psi)^(1-sigma)/(1-sigma)
%
%                s.t.  c[t] + (1+tauxk[t]) xk[t]
%                           + (1+tauxh[t]) q[t] xh[t]
%                         <= (1-tauk1[t]) r1[t] k[t] v[t]
%                           +(1-tauk2[t]) r2[t] k[t] (1-v[t])
%                           +(1-tauh1[t]) w1[t] l[t] h[t] u[t]
%                           +(1-tauh2[t]) w2[t] l[t] h[t] (1-u[t]) + T[t]
%
%                       (1+n) k[t+1] = (1-delk) k[t] + xk[t]
%                       (1+n) h[t+1] = (1-delh) h[t] + xh[t]
%
%                       where r's, w's are marginal products of F,G where
%  
%                       c[t] + xk[t] = y[t] = F(k[t]v[t],h[t]l[t]u[t])
%                              xh[t]        = G(k[t](1-v[t]),h[t]l[t](1-u[t]))
%
%                       F(kv,hlu)           = A(kv)^alpk * (hlu)^alph
%                       G(k(1-v),hl(1-u))   = B(k(1-v))^thetk * (hl(1-u))^theth

%         Ellen McGrattan, 9-24-97
%         Revised, ERM, 9-11-98



par      = input(['Pick parameters: Jones-Manuelli-Rossi=1,', ...
                  'Kim=2,King-Rebelo=3,Lucas=4,Yours=5 ']);
if par==1;
  load par_jones.dat
  par1   = par_jones(:,1);
  par2   = par_jones(:,2);
  z1     = [.43;.29;1.020;.62];
  z2     = [.57;.38;1.052;.62];
elseif par==2;
  load par_kim.dat
  par1   = par_kim(:,1);
  par2   = par_kim(:,2);
  z1     = [4.04;0.9998;1.015;0.72];
  z2     = [5.76;0.9998;1.024;0.67];
elseif par==3;
  load par_king.dat
  par1   = par_king(:,1);
  par2   = par_king(:,2);
  z1     = [2.8;0.9998;1.020;0.51];
  z2     = [2.8;0.9998;1.053;0.38];
elseif par==4;
  load par_lucas.dat
  par1   = par_lucas(:,1);
  par2   = par_lucas(:,2);
  z1     = [2;0.75;1.015;0.9999];
  z2     = [3;0.84;1.016;0.9999];
else 
  load par_mine.dat
  par1   = par_mine(:,1);
  par2   = par_mine(:,2);
  z1     = [1.15;0.75;1.015;0.9941];
  z2     = [2.77;0.84;1.016;0.9941];
end;

z1       = secant('twofoc',z1,par1);
z2       = secant('twofoc',z2,par2);

param    = par1;
z        = z1;
a        = param(1);
alpk     = param(2);
alph     = param(3);
delk     = param(6);
delh     = param(7);
eta      = param(8);
thetk    = param(11);
theth    = param(12);
tauxk    = param(14);
tauxh    = param(15);
tauk1    = param(16);
tauk2    = param(17);
tauh1    = param(18);
tauh2    = param(19);
tauxk1   = 1+tauxk;
tauxh1   = 1+tauxh;
tauk11   = 1-tauk1;
tauk21   = 1-tauk2;
tauh11   = 1-tauh1;
tauh21   = 1-tauh2;
kss      = z(1);
lss      = z(2);
gamss    = z(3);
vss      = z(4);
phi      = tauk11*tauh21*alpk*theth/(tauk21*tauh11*alph*thetk);
uss      = 1/(1-phi+phi/vss);
yss      = a*(kss*vss)^alpk*(lss*uss)^alph;
xhss     = gamss-1+delh;
xkss     = (gamss-1+delk)*kss;
css      = yss-xkss;

disp(' ')
disp('Steady state values for Two-sector Growth Model w/ Taxes')
disp('----------------------------------------------------------')
disp(sprintf('   (H''/H-1)*100                   =   %7.4f ',(gamss-1)*100))
disp(sprintf('   (h''/h-1)*100                   =   %7.4f ', ...
                                                    (gamss/(1+eta/100)-1)*100))
disp(sprintf('   c/y                            =   %7.4f ',css/yss))
disp(sprintf('   l                              =   %7.4f ',lss))
disp(sprintf('   xk/y                           =   %7.4f ',xkss/yss))
disp(sprintf('   xh/y                           =   %7.4f ',xhss/yss))
disp(sprintf('   k/y                            =   %7.4f ',kss/yss))
disp(sprintf('   k/h                            =   %7.4f ',kss))
disp(sprintf('   v                              =   %7.4f ',vss))
disp(sprintf('   u                              =   %7.4f ',uss))
disp('----------------------------------------------------------')


param    = par2;
z        = z2;
a        = param(1);
alpk     = param(2);
alph     = param(3);
delk     = param(6);
delh     = param(7);
thetk    = param(11);
theth    = param(12);
tauxk    = param(14);
tauxh    = param(15);
tauk1    = param(16);
tauk2    = param(17);
tauh1    = param(18);
tauh2    = param(19);
tauxk1   = 1+tauxk;
tauxh1   = 1+tauxh;
tauk11   = 1-tauk1;
tauk21   = 1-tauk2;
tauh11   = 1-tauh1;
tauh21   = 1-tauh2;
kss      = z(1);
lss      = z(2);
gamss    = z(3);
vss      = z(4);
phi      = tauk11*tauh21*alpk*theth/(tauk21*tauh11*alph*thetk);
uss      = 1/(1-phi+phi/vss);
yss      = a*(kss*vss)^alpk*(lss*uss)^alph;
xhss     = gamss-1+delh;
xkss     = (gamss-1+delk)*kss;
css      = yss-xkss;

disp(' ')
disp('Steady state values for Two-sector Growth Model wo/ Taxes')
disp('-------------------------------------------------------------')
disp(sprintf('   (H''/H-1)*100                   =   %7.4f ',(gamss-1)*100))
disp(sprintf('   (h''/h-1)*100                   =   %7.4f ', ...
                                                    (gamss/(1+eta/100)-1)*100))
disp(sprintf('   c/y                            =   %7.4f ',css/yss))
disp(sprintf('   l                              =   %7.4f ',lss))
disp(sprintf('   xk/y                           =   %7.4f ',xkss/yss))
disp(sprintf('   xh/y                           =   %7.4f ',xhss/yss))
disp(sprintf('   k/y                            =   %7.4f ',kss/yss))
disp(sprintf('   k/h                            =   %7.4f ',kss))
disp(sprintf('   v                              =   %7.4f ',vss))
disp(sprintf('   u                              =   %7.4f ',uss))
disp('-------------------------------------------------------------')

