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



Alpha  = [0.0001,.5:.5:50]';
Alpha  = [.5,1.0001,1.5]';
Beta   = linspace(.8,.9999,100);
Beta   = linspace(.95,.99,5);
na     = length(Alpha);
nb     = length(Beta);

mu     = .018;
delta  = .036/2;
phi    = .43;
lam1   = 1+mu+delta;
lam2   = 1+mu-delta;

k      = 0;
for j=1:nb;
  beta = Beta(j);
  for i=1:na;
    k     = k+1;
    alpha = Alpha(i);

    A     = [1-beta*phi*lam1^(1-alpha),-beta*(1-phi)*lam2^(1-alpha);
             -beta*(1-phi)*lam1^(1-alpha),1-beta*phi*lam2^(1-alpha)];
    b     = [beta*phi*lam1^(1-alpha)+beta*(1-phi)*lam2^(1-alpha);
             beta*(1-phi)*lam1^(1-alpha)+beta*phi*lam2^(1-alpha)];
 
    w     = A\b;

    re    = [lam1*(w(1)+1)/w(1)-1,lam2*(w(2)+1)/w(1)-1;
             lam1*(w(1)+1)/w(2)-1,lam2*(w(2)+1)/w(2)-1];
    pf    = [beta*phi*lam1^(-alpha)+beta*(1-phi)*lam2^(-alpha);
             beta*(1-phi)*lam1^(-alpha)+beta*phi*lam2^(-alpha)];
    rf    = [1;1]./pf-1;

    Re(k) = .5*(phi*re(1,1)+(1-phi)*re(1,2))+.5*((1-phi)*re(2,1)+phi*re(2,2));
    Rf(k) = .5*rf(1)+.5*rf(2);
  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,4,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)
ii1 = find(Alpha<=5 & Beta<=.9);
ii2 = find(Alpha<=5 & Beta>.9);
ii3 = find(Alpha>5 & Beta<=.9);
ii4 = find(Alpha>5 & Beta>.9);
plot(Rf(ii1),Re(ii1)-Rf(ii1),'green+',Rf(ii2),Re(ii2)-Rf(ii2),'yellow+', ...
     Rf(ii3),Re(ii3)-Rf(ii3),'red+', Rf(ii4),Re(ii4)-Rf(ii4),'blue+')
text(-9,.5,'\alpha\leq 5, \beta> .9')
text(32,.5,'\alpha\leq 5, \beta\leq .9')
text(-9,12,'\alpha> 5, \beta> .9')
text(34,12,'\alpha> 5, \beta\leq .9')
title('Equilibrium Risk Premia and Risk-Free Rates in Mehra and Prescott')
ylabel('Average Risk Premia (percent)')
xlabel('Average Risk Free Rate (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)
