%MEHRA     computes the equity premium for the model in Mehra and
%          Prescott, JME 1985.



Alpha  = [0.0001,.5:.5:2]';
Beta   = linspace(.9,.98,10);
na     = length(Alpha);
nb     = length(Beta);

mu     = .018;
delta  = .036;
Phi    = [.43,.57;
          .57,.43];
lam    = [1+mu+delta;
          1+mu-delta];

Phi    = [ 0,1,0;
          .05,.43-.05,.57;
          .00,.57,.43-.0];
lam    = [.1;
          1+mu+delta;
          1+mu-delta];


%[d,Phi]= tauchen2(.14,.012,2);
%lam    = 1+mu+d(:);

%lam    = [1+mu+.036;
%          1+mu-.036;
%          1+mu+.036;
%          1+mu-.036;
%          1+mu+.036;
%          1+mu-.036];
%Phi    = [Phi*.98,Phi*.02,Phi*0;
%          Phi*.0,Phi*.98,Phi*.02
%          Phi*.001,Phi*.0,Phi*.999];

Pi     = ergodic(Phi);
n      = length(Phi);
e      = ones(n,1);
k      = 0;
for j=1:nb;
  beta = Beta(j);
  for i=1:na;
    k     = k+1;
    alpha = Alpha(i);

    A     = eye(n)-beta*Phi.*(e*lam').^(1-alpha);
    b     = beta*Phi*lam.^(1-alpha);
    w     = A\b;
    re    = (e./w) * (lam.*(w+1))'-1;
    pf    = beta *Phi*lam.^(-alpha);
    rf    = e./pf-1;

    Re(k) = sum((Phi.*re)')*Pi;
    Rf(k) = Pi'*rf;

  end;
end;          
  
Rf = Rf(:)*100;
Re = Re(:)*100;

EP = Re-Rf;

Alpha = Alpha*ones(1,nb);
Alpha = Alpha(:);
Beta  = ones(na,1)*Beta;
Beta  = Beta(:);

figure(1)
plot(Rf,Re-Rf,'+'),axis([0,5,0,2])
title('Mehra and Prescott (1985) -- Figure 4')
ylabel('Average Risk Premia R^e-R^f (percent)')
xlabel('Average Risk Free Rate R^f (percent)')

figure(2)
plot(Rf,Re-Rf,'+')
title('Mehra and Prescott (1985) -- Figure 4 -- all points')
ylabel('Average Risk Premia R^e-R^f (percent)')
xlabel('Average Risk Free Rate R^f (percent)')
v=axis;

figure(3)
subplot(221)
ii = find(Alpha<5);
plot(Rf(ii),Re(ii)-Rf(ii),'+')
title('\alpha < 5')
axis(v)
subplot(222)
ii = find(Alpha>5);
plot(Rf(ii),Re(ii)-Rf(ii),'+')
title('\alpha > 5')
axis(v)
subplot(223)
ii = find(Beta<.9);
plot(Rf(ii),Re(ii)-Rf(ii),'+')
title('\beta < .9')
axis(v)
subplot(224)
ii = find(Beta>.9);
plot(Rf(ii),Re(ii)-Rf(ii),'+')
title('\beta > .9')
axis(v)
figure(4)
ii = find(Alpha>.99 & Alpha<3.01 & Beta>.95 & Beta<1);
jj = find(Alpha>.99 & Alpha<1.01 & Beta>.977 & Beta<.983);
plot(Rf(ii),Re(ii)-Rf(ii),'+',Rf(jj),Re(jj)-Rf(jj),'*')
title('.99 < \alpha <= 3 and .95 < \beta < 1 with * for \alpha=1,\beta=.98')
ylabel('Average Risk Premia R^e-R^f (percent)')
xlabel('Average Risk Free Rate R^f (percent)')
