%TWOSS    calculates steady state of model described in twosect.f.

%         Ellen McGrattan, 6-26-97


% Parameters for:
%    U(c,l,g)   = [c*l^psi*g^zeta]^(1-sigma) / (1-sigma)
%    F(k,l,g)   = a*k^alpk *l^alph *g^(1-alpk-alph) 
%    G(k,l,g)   = b*k^thetk*l^theth*g^(1-thetk-theth) 
%    beta        = discount factor
%    delk        = depreciation of physical capital
%    delh        = depreciation of human capital
%    eta         = population growth
%

a     = 1.0;
alpk  = .33;
alph  = .67;
b     = .10;
beta  = .95;
delk  = .06;
delh  = .06;
eta   = 0.;
psi   = .5;
sigma = 2.0001;
thetk = .33;
theth = .67;
zeta  = 0;

%
% Policy variables:
%

tauxk = .4;
tauxh = .0;
tauk1 = .0;
tauk2 = .0;
tauh1 = .0;
tauh2 = .0;
tauc  = .0;
z1    = 1e-8;
z2    = 1e-8;
z3    = 1e-8;

par   = [a;alpk;alph;b;beta;delk;delh;eta;psi;sigma;thetk;theth;zeta;
         tauxk;tauxh;tauk1;tauk2;tauh1;tauh2;tauc;z1;z2;z3];

% Guesses for k/y, l, gam, v:

ky    = 1;
l     = 0.25;
gam   = 1;
v     = theth;

% Use guesses of ratios to get initial guess in levels:  

phi   = (1-tauk1)*(1-tauh2)*alpk*theth/((1-tauk2)*(1-tauh1)*alph*thetk);
u     = 1/(1-phi+phi/v);
y     = (a*(ky*v)^alpk*(l*u)^alph*z2^(1-alpk-alph) )^(1/alph);
k     = ky * y;

zi    = [k;l;gam;v];


% Call to fixed point algorithm

z     = secant('twofoc',zi,par,1e-8,100);
z     = secant('twofoc',real(z),par,1e-8,100);

k     = z(1);
l     = z(2);
gam   = z(3);
v     = z(4);
xh    = gam-1+delh;
xk    = (gam-1+delk)*k;
u     = 1/(1-phi+phi/v);
y     = (a*(k*v)^alpk*(l*u)^alph*z2^(1-alpk-alph))^(1/(alpk+alph));
c     = (1-z1-z2-z3)*y-xk;
q     = ((1-tauh1) * alph * y/u )/( (1-tauh2)* theth*xh/(1-u) );
gdp   = y+q*xh;

disp('Steady state values for Two-sector Growth Model')
disp('-----------------------------------------------')
disp(sprintf('   (H''/H-1)*100             =   %7.4f ',(gam-1)*100))
disp(sprintf('   (h''/h-1)*100             =   %7.4f ',(gam/(1+eta/100)-1)*100))
disp(sprintf('   c/y                      =   %7.4f ',c/y))
disp(sprintf('   c/(y+q*xh)               =   %7.4f ',c/gdp))
disp(sprintf('   g/y                      =   %7.4f ',(z1+z2+z3)))
disp(sprintf('   g/(y+q*xh)               =   %7.4f ',(z1+z2+z3)*y/gdp))
disp(sprintf('   l                        =   %7.4f ',l))
disp(sprintf('   xk/y                     =   %7.4f ',xk/y))
disp(sprintf('   xk/(y+q*xh)              =   %7.4f ',xk/gdp))
disp(sprintf('   k/y                      =   %7.4f ',k/y))
disp(sprintf('   k/(y+q*xh)               =   %7.4f ',k/gdp))
disp(sprintf('   xh/y                     =   %7.4f ',xh/y))
disp(sprintf('   xh/(y+q*xh)              =   %7.4f ',xh/gdp))
disp(sprintf('   q*xh/(y+q*xh)            =   %7.4f ',q*xh/gdp))
disp(sprintf('   k/h                      =   %7.4f ',k))
disp(sprintf('   v                        =   %7.4f ',v))
disp(sprintf('   u                        =   %7.4f ',u))
disp('-----------------------------------------------')

