%%%
%%% Rietz-Barro Model
%%%

clear all ; close all; clc ;

beta = .97 ; % discount factor
p = .017 ; % probability of a disaster
b = .4 ; % size of disaster
mu = .025 ; % mean growth of consumption (w/o disasters)
sigma = .02 ; % std dev of consumption growth
rrr = 0.8 ; % repayment of bond if disaster (for question 4)

uuu = exp(-p);
listgamma = linspace(.5,6,400);

ep = zeros(size(listgamma));
rf = zeros(size(listgamma));
ret = zeros(size(listgamma));
rfrepay = zeros(size(listgamma));
eprepay = zeros(size(listgamma));
epnodis = zeros(size(listgamma));
rfnodis = zeros(size(listgamma));

for k=1:length(listgamma)

    gamma = listgamma(k);
    
    %%% equity premium in basic case (question 2,3)
    ep(k) = 100*log( exp(gamma*sigma^2)*(uuu + (1-uuu)*(1-b)^(-gamma))*(uuu+(1-uuu)*(1-b))/(uuu+(1-uuu)*(1-b)^(1-gamma)) ) ;
    %%% risk free rate (question 1)
    rf(k) = 100 / (beta*exp(-gamma*mu+gamma^2*sigma^2/2)*(uuu+(1-uuu)*(1-b)^(-gamma))) - 100 ;
    %%% return on equity (question 2,3)
    ret(k) = exp(ep(k)/100)*(rf(k)/100+1);
    %%% equity premium if no disaster:
    epnodis(k) = 100*gamma*sigma^2;
    %%% risk free rate if no disaster:
    rfnodis(k) = 100 / (beta*exp(-gamma*mu+gamma^2*sigma^2/2)) - 100 ;
    %%% risk free rate if bond is repaid only with some probability rrr in
    %%% disasters (question 5)
    pb =  (beta*exp(-gamma*mu+gamma^2*sigma^2/2)*(uuu+rrr*(1-uuu)*(1-b)^(-gamma))) ;
    rfrepay(k) = (uuu+(1-uuu)*rrr)/pb*100 - 100 ;
    %%% equity premium if bond is repaid only with some probability rrr in
    %%% disasters
    eprepay(k) = 100*log(ret(k)/(rfrepay(k)/100+1));
    
    %%% risk free rate and equity premium if 'disasters' are symetric ie exceptional events
    %%% can go up or down (question 4)
    rfsym(k) = 100 / (beta*exp(-gamma*mu+gamma^2*sigma^2/2)*(uuu+(1-uuu)/2*(1-b)^(-gamma)+(1-uuu)/2*(1-b)^(gamma))) - 100 ;
    epsym(k) = 100*log( exp(gamma*sigma^2)*(uuu + (1-uuu)/2*(1-b)^(-gamma)+ (1-uuu)/2*(1-b)^(gamma))*(uuu+(1-uuu)/2*(1-b)+(1-uuu)/2*1/(1-b))/(uuu+(1-uuu)/2*(1-b)^(1-gamma)+(1-uuu)/2*(1-b)^(gamma-1)) ) ;

end

figure; plot(listgamma,ep,'b-'); title('Equity Premium in the Rietz-Barro model');
hold on ; plot(listgamma,epnodis,'k-.'); legend('disasters','no disasters',2)
xlabel('risk aversion'); ylabel('% per year'); grid on;

figure; plot(listgamma,rf,'b-'); title('Risk Free Rate in the Rietz-Barro model');
hold on ; plot(listgamma,rfnodis,'k-.'); legend('disasters','no disasters',2)
xlabel('risk aversion'); ylabel('% per year') ;grid on

figure; plot(listgamma,ep,'b-'); title('Equity Premium in the Rietz-Barro model');
hold on ; plot(listgamma,epnodis,'k-.');
hold on ; plot(listgamma,epsym,'r--');legend('disasters','no disasters','symetric extreme events',2)
xlabel('risk aversion'); ylabel('% per year'); grid on;

figure; plot(listgamma,ep,'b-'); title('Equity Premium in the Rietz-Barro model');
hold on ; plot(listgamma,eprepay,'k-.'); legend('RF bonds','r=.8',2)
xlabel('risk aversion'); ylabel('% per year'); grid on;


%%%
%%% model with 'recoveries' (question 6)
%%%

pilist = linspace(0,1,100);
neweplist = zeros(size(pilist)) ;

for k=1:length(pilist);
    pi = pilist(k);
    
    %%% change this number: note how it is different for gamma>1 or gamma<1
    gamma = 5 ; 
    %gamma  = 4.67 ;
    
    cst = beta*exp((1-gamma)*mu+(1-gamma)^2*sigma^2/2);

    %%% write the system of equations as A*[phi0,phi1]' = B
    AA = zeros(2); BB = zeros(2,1);
    AA(1,1) = cst*( (1-pi)*uuu + (1-b)^(gamma-1)*pi*uuu) ; %%% COEFF ON PHI(0)
    AA(1,2) = - 1 + cst*(pi*(1-uuu)+(1-b)^(1-gamma)*(1-pi)*(1-uuu)) ; %%% coeff ON PHI(1)
    BB(1) = - (AA(1,1)+AA(1,2)+1);
    AA(2,1) = -1 + cst*uuu;
    AA(2,2) = cst*(1-b)^(1-gamma)*(1-uuu);
    BB(2) = -(1+AA(2,1)+AA(2,2));

    phi = inv(AA)*BB;
    phi0 = phi(1);
    phi1 = phi(2);

    expret1 = exp(mu+sigma^2/2)*( (1-pi)*uuu*(phi0+1) + pi*(1-uuu)*(phi1+1) + pi*uuu/(1-b)*(phi0+1) + (1-pi)*(1-uuu)*(phi1+1)*(1-b) )/phi1 ;
    expret0 = exp(mu+sigma^2/2)*( uuu*(phi0+1) + (1-uuu)*(phi1+1)*(1-b) )/phi0 ;

    ret = expret1*(1-uuu) + uuu*expret0 ;

    rf1 = 1 / (beta*exp(-gamma*mu+gamma^2*sigma^2/2)*( (1-pi)*uuu +pi*(1-uuu)+pi*uuu*(1-b)^gamma+(1-pi)*(1-uuu)*(1-b)^(-gamma) ) );
    rf0 = 1 / (beta*exp(-gamma*mu+gamma^2*sigma^2/2)*( uuu + (1-uuu)*(1-b)^(-gamma) ) );

    meanrf = (1-uuu)*rf1 + uuu*rf0 ;

    meanep = 100*((1-uuu)*(expret1/rf1)+uuu*expret0/rf0) - 100 ;

    neweplist(k) = meanep ;
    neweplist2(k) = 100*expret0/rf0 - 100 ;
    neweplist1(k) = 100*expret1/rf1 - 100 ;
    
end

figure;plot(pilist,neweplist,'b-');
title('Equity Premium in the Rietz-Barro Model with recoveries')
xlabel('probability of a recovery');
ylabel('equity premium, in %')

figure;plot(pilist,neweplist1,'b-');
title('Equity Premium in the Rietz-Barro Model with recoveries, after a disaster')
xlabel('probability of a recovery');
ylabel('equity premium, in %')

figure;plot(pilist,neweplist2,'b-');
title('Equity Premium in the Rietz-Barro Model with recoveries, no disaster in previous period')
xlabel('probability of a recovery');
ylabel('equity premium, in %')