ns=2;           @ specify the no. of states @
 

@ specify the Markov Transition probabilities @
fi11=0.43;
fi12=1-fi11;
fi22=0.43;
fi21=1-fi22;

fi=zeros(ns,ns);
fi[1,1]=fi11; fi[1,2]=fi12;
fi[2,1]=fi21; fi[2,2]=fi22;

@ Here we solve for the stationary probabilities of the states
  we solve the system of equations: pi = fi'*pi st. sum(pi)=1
  where pi is the vector of stationary probabilities. @
pi1=fi';
pi2=pi1[1:ns-1,1:ns-1]-ones(ns-1,1)*pi1[1:ns-1,ns]';
pi3=pi1[1:ns-1,ns];
plr=(inv(eye(ns-1)-pi2))*pi3;
pst=zeros(ns,1);
pst[1:ns-1]=plr; pst[ns,1]=1-sumc(plr);

@ specify the mean growth rates for the ns states of nature @
l1=1.054;
l2=0.982;
l=zeros(ns,1);
l[1,1]=l1; l[2,1]=l2;

@ create a vector with different values of alfa - the
  coefficient of relative risk aversion @
alfa=seqa(0.01,0.025,400);

@ create a vector with different values of beta @
bstar=seqa(0.999999999,-0.000075,10);

rft=zeros(rows(alfa),rows(bstar));
re=zeros(rows(alfa),rows(bstar));
ep=zeros(rows(alfa),rows(bstar));

@ now compute risk-free rate, equity return, and the
  equity premium for all combinations of the risk
  aversion parameter and the discount factor @
i=1;
 do until i>rows(alfa);
  alfa1=alfa[i,1];
  j=1;
   do until j>rows(bstar);
    beta=bstar[j,1];

    @ compute the uncontional expectation of the risk free rate @
    pf=beta*fi*(l^(-alfa1));
    rfs=pf^(-1)-1;
    rf=pst'*rfs;

    @ set up the matrices for the computation of the
      constants w @
    la=l^(1-alfa1);
    lb=eye(ns).*la;
    x=beta*fi*lb;
    f=eye(ns)-x; g=sumc(x');
    w=(inv(f))*g;

    @ compute the unconditional expectation of equity return @
    wa=w+ones(ns,1); wb=(wa./w')';
    la=eye(ns).*l;
    rss=wb*la-ones(ns,ns);
    rs=sumc((fi.*rss)');

    rft[i,j]=rf;
    re[i,j]=pst'*rs;
    ep[i,j]=100*(re[i,j]-rf);

    j=j+1;
   endo;
  i=i+1;
 endo;
 

ep1=zeros(rows(alfa),rows(bstar));
rf1=zeros(rows(alfa),rows(bstar));
 

@ now we eliminate those observations which imply a risk-free
  rate above 4 percent @
i=1;
do until i>rows(rft);

 if rft[i,1]<0.04;

  j=1;
   do until J>cols(rft);

    if rft[i,j]<0.04;
       rf1[i,j]=rft[i,j];
       ep1[i,j]=ep[i,j];
       j=j+1;
    else; goto stop1;
    endif;
   endo;

  stop1:

 endif;
i=i+1;
endo;

@ now we put all the risk free rates and equity premia
  in one long vector @
rf2=reshape(rf1,rows(rf1)*cols(rf1),1);
ep2=reshape(ep1,rows(rf1)*cols(rf1),1);

@ now we sort the observations in terms of increasing risk free
  rates @
ind=sortind(rf2);
rf3=rf2[ind];
ep3=ep2[ind];

@ finally we plot the result @
library pgraph;
title("the Equity premium in the Mehra-Prescott, 1985, economy");
xlabel("the risk-free rate (%)");
ylabel("the equity premium (%)");

bar(100*rf3,ep3);