@This program solves the model in Imrohoroglu (1989) without aggregate
    uncertainty.@
@goto start;@
new;
library pgraph;
output file=teach.out reset;
sigma=1.5;   @Risk aversion parameter.@
beta=.995;    @Discount factor.@
pee=.9565;    @Probability of working when employed last period.@
peu=.5;      @Probability of working when unemployed last period.@
maxast=4;   @Maximum number of assets.@
incast=.1;  @Minimum number of assets.@
theta=.25;     @Income that agent gets when unemployed.@
nsim=50000;   @Number of simulations at end.@

@Tabulate return function@
nasts=round(maxast/incast+1);
util1=-10000*ones(nasts,nasts);  @Utility when employed.@
util2=-10000*ones(nasts,nasts);  @Utility when unemployed.@
a=0;
do until a > maxast+.0001;
     i=round(a/incast+1);
     apr=0;
     do until apr > a+1-.0001 or apr > maxast+.0001;
          j=round(apr/incast+1);
          util1[j,i]=(1+a-apr)^(1-sigma)/(1-sigma);
          apr=apr+incast;
     endo;
     apr=0;
     do until apr > a+theta-.0001 or apr > maxast+.0001;
          j=round(apr/incast+1);
          util2[j,i]=(theta+a-apr)^(1-sigma)/(1-sigma);
          apr=apr+incast;
     endo;
     a=a+incast;
endo;

@Iterate on Belman's equation.@
v=zeros(nasts,2);
test=10;
iter=0;
do until test < 1e-4 or iter > 2000;
     r1=util1+beta*(pee*v[.,1]+(1-pee)*v[.,2]);
     r2=util2+beta*(peu*v[.,1]+(1-peu)*v[.,2]);
     tv=maxc(r1)~maxc(r2);
     test=abs(tv-v);
     v=tv;
     iter=iter+1;
 endo;
 "Dynamic Program took";; iter;; "iterations";
decis=maxindc(r1)~maxindc(r2);
decisr=(decis-1)*incast;
decisr;
x=seqa(0,incast,nasts);
xy(x,decisr);
start:
@Compute invariant distribution.@
i=1;
d1=zeros(nasts,nasts); d2=d1;
do until i > nasts;
     d1[decis[i,1],i]=1;
     d2[decis[i,2],i]=1;
     i=i+1;
endo;

f=ones(nasts,2)/(2*nasts);   @Initial guess.@
diter=0;
test=10;
do until test < .000001 or diter > 500;
     fpr=(pee*d1*f[.,1]+peu*d2*f[.,2])~((1-pee)*d1*f[.,1]+(1-peu)*d2*f[.,2]);
     test=abs(fpr-f);
     f=fpr;
     diter=diter+1;
endo;
xy(x,f);
@Compute average asset holdings using invariant distribution.@
fa=f[.,1]+f[.,2];
assets=seqa(0,incast,nasts);
massets=sumc(assets.*fa);
"Mean Asset Holdings Using Invariant Distribution ";;massets;
@Compute average asset holding by simulating decision rules.
     Simualate nsim periods, and throw away the first 1000.@
a=massets;  @initial asset holding.@
s=1;
state=rndu(nsim,1);
suma=zeros(nsim,1);
i=1;
do until i > nsim;
     ia=round(a/incast+1);
     if s == 1; prob=pee; else; prob=peu; endif;
     if state[i,1] < prob; s=1; apr=decisr[ia,1];
     else; s=2; apr=decisr[ia,2]; endif;
     suma[i,1]=apr;
     a=apr;
     i=i+1;
endo;
masset=meanc(suma[101:nsim]);
"Mean Asset Holding from Simulation";;masset;
stop:
end;