% ch4ss.m
% Based on xssir.m program from EER (93) paper
% Figure 4.8
rand('seed',0)
T=3500;
eps=.25;
sig=.1;
a=.025;
alp=.9;
lam=.5;
k=40;
istar=14.1935;
bet=1;
sub=0;
del=.5;
pi1=.98;
pi2=.995;
capa=.0792046;
ss=0;
n=zeros(T,1);
xe1=1^(1+eps)/alp ;
xe2=2.5^(1+eps)/alp;
q=zeros(T,1);
for j=1:T;
r=rand(1);
if j==2000,
        sub=.05;
end;
if ss==1,
        if r>pi1,
        ss=2;
        end;
else,
        if r>pi2,
        ss=1;
        end;
end;
if j==1,
        ss=2;
end;
ne1=(alp*xe1)^(1/(1+eps));
ne2=(alp*xe2)^(1/(1+eps));
psi1=capa*(max(istar,lam*k*ne1/(1+a*lam*k*ne1)))^bet;
psi2=capa*(max(istar,lam*k*ne2/(1+a*lam*k*ne2)))^bet;
q1=ne1^alp*psi1;
q2=ne2^alp*psi2;
if ss==1,
x(j)=(pi1*q1^(1-sig)+(1-pi1)*q2^(1-sig))/(1-sub);
xe1=xe1+del*(x(j)-xe1);
else,
x(j)=(pi2*q2^(1-sig)+(1-pi2)*q1^(1-sig))/(1-sub);
xe2=xe2+del*(x(j)-xe2);
end;
n(j)=(alp*x(j))^(1/(1+eps));
psi=capa*(max(istar,lam*k*n(j)/(1+a*lam*k*n(j))))^bet;
q(j)=psi*n(j)^alp;
end;
w=zeros(T-1,1);
for j=1:T-1;
w(j)=-n(j)^(1+eps)/(1+eps)+q(j+1)^(1-sig)/(1-sig);
end;
plot(n);
title('Figure 4.8')
xlabel('time')
ylabel('employment')