% programs for problem set 5 -- consumption model 

close all; 
clear all; 

shortprint = 0; % print out short or long answers

% ----------------------
% consumption over time problem
% -------------------------

% interest rate graph
%dc = kron([ 2 1 -1 0 2]', ones(4,1));
%dc = [NaN;dc]; 

tim = (0:20)';
c = [ 0 2 4 6 8 9 10 11 12 11 10 9 8 8 8 8 8 10 12 14 16 ]';
figure
plot(tim, c,'-v','linewidth',2);
hold on
plot([1;1]*(4:4:20),[-2 16],':k');
plot([0;20],[1;1]*(0:2:14),':k');
legend('c_t',2); 
print -depsc2 consumption_rates_problem_1.eps;

dc = [NaN; c(2:end)-c(1:end-1)];
figure
plot(tim, c,'-v','linewidth',2);
hold on
plot(tim, dc,'-*r','linewidth',2);
plot(tim(1:end-1), 2+2*dc(2:end),'-ok','linewidth',2)
plot([1;1]*(4:4:20),[-2 16],':k');
plot([0;20],[1;1]*(0:2:14),':k');
legend('c_t','\Deltac_t','r^f_t = 2+2*\Deltac_{t+1}',2); 
print -depsc2 consumption_rates_problem_2.eps


% -------------------------
% plot of c-x utility
% -----------------------
c = (1.01:0.01:3)';
c2 = (0.01:0.01:3)'; 
x = 1;
figure
plot(c,((c-x).^-1)/-1,'linewidth',2);
hold on
plot(c2,(c2.^-1)/-1,'-r','linewidth',2);
plot([1 1],[-20 20],'--k'); 
plot([0 3], [0 0], '-k');
text(1.02,0,'x');
axis([ 0 3 -10 1]); 
xlabel('c');
ylabel('u(c)');
print -depsc2 ps4_habit.eps

figure
plot(c,2*(c./x) ./(c./x-1),'linewidth',2);
hold on
plot(c2,2+0*c2,'-r','linewidth',2);
plot([1 1],[-20 20],'--k'); 
text(1.02,1,'x');
axis([ 0 3 0 10]); 
xlabel('c');
ylabel('risk aversion');
print -depsc2 ps4_habit_risk_aversion.eps



% ----------------------------------------
% plot of p/c vs. g for various gamma
% ----------------------------------------

gamma = [0.5 0.8 1 2 5]'; % vertical vector of gammas to try
beta = 0.95; 
g = (0.95:0.001:1.2)'; %horizontal vector of g to try

bg = g*ones(1,size(gamma,1));   % bg for big g. Before using .* they must be the same shape. 
                                % This repeats g side by side
bgam = ones(size(g,1),1)*gamma'; % similarly, repeats gamma up and down. 
pc = beta*bg.^(1-bgam)./(1-beta*bg.^(1-bgam)); % now this produces a table where gammas change 
                                                % up and down, and g changes side to side

figure; 
plot((g-1)*100,pc); 
legend('\gamma = 0.5','\gamma=0.8','\gamma=1','\gamma=2','\gamma=5'); 
axis([-inf inf 0 100]); 
xlabel('growth in percent 100(g-1)'); 
ylabel('p/c'); 
print -depsc2 ps4_pc.eps; 

% plot of p/c for various gammas, continuous time version 

figure
delta = -log(beta);
blg = log(bg); % big log g
pc2 = 1./(delta-(1-bgam).*blg);
plot(100*log(g),pc2);
xlabel('growth in percent 100(g-1)'); 
ylabel('p/c');
title('continuous time p/c')
legend('\gamma = 0.5','\gamma=0.8','\gamma=1','\gamma=2','\gamma=5'); 
axis([-inf inf 0 100]); 


stop

% -------------------------
% Ignore what follows. IT's a great consumption question that I did not ask
% this year. 
%--------------------------

% ----------------------------------------
% consumption model

% Sudents who don't use matlab: The main thing you need to know to read the
% program is that the syntax A = [ 1 2 ; 3 4] creates a matrix 
% A = [ 1 2 ] 
%     [ 3 4 ]; 
% then you can access matrix elements with () notation. For example
% A(1,2)=2; A(2,2)=4. 
% * is matrix multiplication. .* and ./ are element by element
% multiplication, division, etc. 
% This should show you the computations to make, and translate to excel. 
% ----------------------------------------

u = 1.06; 
d = 0.98; 
beta = 0.99; 
gamma = 5; 
x = -0.3; 
piuu = (1+x)/2;
piud = (1-x)/2; 
pidu = (1-x)/2; 
pidd = (1+x)/2; 

% notice all terms are beta * piuu * u^-gamma etc. Call these buu etc. 
buu = beta*piuu*u^(-gamma);
bud = beta*piud*d^(-gamma);
bdu = beta*pidu*u^(-gamma);
bdd = beta*pidd*d^(-gamma);

A = [ 1-buu  -bud;  % the syntax [..  ;..  ] forms a 2 x 2 matrix A with the indicated elments
      -bdu  1-bdd ]; 

% perpetuity
pp = inv(A)*[buu+bud; ...% the [ buu+bud; bdu+bdd] forms a 2x1 matrix on the right hand side. 
             bdu+bdd];   % the result pp on the left hand side is a 2x1 matrix with the price 
                         % in the u state on top of the price in the d
                         % state. 

% interest rate
Rf = 1./[buu+bud; bdu+bdd]; 

% stock p/c ratio

cuu = beta*piuu*u^(1-gamma);
cud = beta*piud*d^(1-gamma);
cdu = beta*pidu*u^(1-gamma);
cdd = beta*pidd*d^(1-gamma);

A = [ 1-cuu  -cud;
      -cdu  1-cdd ]; 
pc = inv(A)*[cuu+cud;...
             cdu+cdd]; 

% expected return 

ER = [piuu*(pc(1)+1)*u/pc(1)+piud*(pc(2)+1)*d/pc(1); % Expected stock return 
      pidu*(pc(1)+1)*u/pc(2)+pidd*(pc(2)+1)*d/pc(2)];
ERp = [piuu*(pp(1)+1)/pp(1)+piud*(pp(2)+1)/pp(1); % Expected bond return. Result is Er in u state on top of ER in d state
       pidu*(pp(1)+1)/pp(2)+pidd*(pp(2)+1)/pp(2)]; 

Edc = [piuu*u+piud*d;  % expected consumption growth 
       pidu*u+pidd*d]; 
phi = (Edc(2)-Edc(1))/(d-u); 

% actual returns. From on rows, to on columns

% R_t+1 = (P_t+1+D_t+1)/P_t 
% = ( P_t+1/D_t+1 + 1)*D_t+1/D_t /(P_t/D_t)
   
pc1 = [1;1]*pc';  
dc1 = [1;1]*[ u d]; 
pc0 = pc*[1 1]; 
Rs = (pc1+1).*dc1./pc0; 
Rs=(Rs-1)*100;


% construct a 2 x 2 matrix with actual returns. 
% actual returns depend on which state you start from 
% as well as which state you go to. Thus organize it as
%   [ u -> u   u -> d]
%   [ d -> u   d -> d]
pp1 = [1;1]*pp';  % creates tomorrow's price [pp(u)  pp(d)]
                  %                          [pp(u)  pp(d)] 
pp0 = pp*[1 1];   % creates today's price    [pp(u)  pp(u)]
                  %                          [pp(d)  pp(d)]
Rp = (pp1+1)./pp0; % return = tomorrow's price plus dividend / today's price
                   % note element by element division. Result is 
                   %                         [(pp(u)+1)/pp(u)  (pp(d)+1)/pp(d)]
                  %                          [(pp(u)+1)/pp(d)  (pp(d)+1)/pp(d)] 
Rp = (Rp-1)*100;  % translate units to percent returns. 

% display results
% note: disp is much easier to use than fprintf. I use fprintf to make nice
% tables for your solutions, but disp plus cut and paste may be much faster
% for you. Try either
% format short
% or
% format compact
% then
% disp('gamma and x'); 
% disp([gamma x]); 
% for quick output. 

if shortprint; 
fprintf('Results for  %6s %6s \n', 'gamma', 'x'); 
fprintf('             %6.2f %6.2f \n',[gamma x]); 

fprintf('             %6s %6s %6s %6s %6s  \n','Bond p','p/c','Rf','ER-Rf','ERp-Rf');
fprintf('good state u %6.2f %6.2f %6.2f %6.2f %6.2f  \n',[ pp(1) pc(1) 100*(Rf(1)-1) 100*(ER(1)-Rf(1)) 100*(ERp(1)-Rf(1))]); 
fprintf('bad  state d %6.2f %6.2f %6.2f %6.2f %6.2f  \n',[ pp(2) pc(2) 100*(Rf(2)-1) 100*(ER(2)-Rf(2)) 100*(ERp(2)-Rf(2))]);

fprintf('             %12s  %12s  \n','stock excess','bond excess');
fprintf(' To state:  %6s %6s %6s %6s \n','u','d','u','d'); 
fprintf('   From u:  %6.2f %6.2f %6.2f %6.2f \n',  (Rs(1,1:2)-100*(Rf(1)-1)),(Rp(1,1:2)-100*(Rf(1)-1)));
fprintf('   From d:  %6.2f %6.2f %6.2f %6.2f  \n', (Rs(2,1:2)-100*(Rf(2)-1)),(Rp(2,1:2)-100*(Rf(2)-1)));

else; 
fprintf('Results for %6s %6s %6s %6s %6s %6s \n','u', 'd', 'beta', 'gamma', 'x','phi'); 
fprintf('            %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n',[u d beta gamma x phi]); 

fprintf('             %6s %6s %6s %6s %6s %6s %6s %6s %6s \n','Dc','E dC','Bond p','p/c','Rf','ER','ERp','ER-Rf','ERp-Rf');
fprintf('good state u %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n',[ (u-1)*100 (Edc(1)-1)*100 pp(1) pc(1) 100*(Rf(1)-1) 100*(ER(1)-1) 100*(ERp(1)-1) 100*(ER(1)-Rf(1)) 100*(ERp(1)-Rf(1))]); 
fprintf('bad  state d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n',[ (d-1)*100 (Edc(2)-1)*100 pp(2) pc(2) 100*(Rf(2)-1) 100*(ER(2)-1) 100*(ERp(2)-1) 100*(ER(2)-Rf(2)) 100*(ERp(2)-Rf(2))]);

fprintf('             %12s  %12s  %12s  %12s \n','Stock Return','stock excess','Bond Return','bond excess');
fprintf(' To state: %6s %6s %6s %6s %6s %6s %6s %6s \n','u','d','u','d','u','d','u','d'); 
fprintf('   From u:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n',  Rs(1,1:2),(Rs(1,1:2)-100*(Rf(1)-1)),Rp(1,1:2),(Rp(1,1:2)-100*(Rf(1)-1)));
fprintf('   From d:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f  \n', Rs(2,1:2),(Rs(2,1:2)-100*(Rf(2)-1)),Rp(2,1:2),(Rp(2,1:2)-100*(Rf(2)-1)));
end; 

% -------------------------------------------------------------------------
% The rest is not part of the homework. I wanted to try a simulation of the model
% to watch how stocks and bonds worked over time. It wasn't pretty enough
% to make you do it. 
% -------------------------------------------------------------------------

% simulate  

rand('state',0);  % initializes random number generator so it's the same each time
                  % that helps fine tune graphs etc. 
                  
% start in state 1=u ; (state 2 = d) 
st = 1; 
for t = 2:50; 
    error = rand(1,1); % now go up or down according to probability 
    if st(t-1) == 1; 
        if error < piuu; 
            st(t) = 1; 
        else; 
            st(t) = 2; 
        end; 
    else; 
        if error < pidu; 
            st(t) = 1; 
        else; 
            st(t) = 2; 
        end; 
    end; 
end; 

% now look at prices etc through time
dc = [u;d];   % creates vertical vector with u over d
dct = dc(st);  % creaates a time series of consumption growth rates. 
               % st is a vector like [1 2 2 2 1 1 2 ...]'. If you use it as an index,
               % it creates a vector  
ppt = pp(st); 
pct = pc(st); 
Rft = Rf(st);

% R_t+1 = (P_t+1+D_t+1)/P_t 
% = ( P_t+1/D_t+1 + 1)*D_t+1/D_t /(P_t/D_t)
Rt = (pct(2:end)+1).*dct(2:end)./pct(1:end-1); 

% Rp_t+1 = (pp_t+1 + 1)/pp_t
Rpt = (ppt(2:end)+1)./ppt(1:end-1); 

% cumulate consumption growth from 1 to get consumption level
ct = cumprod(dct); 
pt = pct.*ct; % stock price level

% cut off the prices so they all can be plotted together
dct = dct(2:end);  
ppt = ppt(2:end); 
pct = pct(2:end); 
Rft = Rft(2:end);
pt = pt(2:end); 
ct = ct(2:end); 

figure; 
plot([ppt pct]); 
legend('perpetuity price','stock p/c'); 

figure; 
plot([log(ct) log(pt)]); 
legend('log consumption','log stock price level'); 

figure; 
plot([(Rt-1)*100 (Rft-1)*100 (Rt-Rft)*100 (dct-1)*100]); 
legend('Stock return','Interest rate','excess return','Consumption'); 
title('Percent returns'); 
axis([0 20 -inf inf]); 