% fama bliss cochrane piazzesi problem set 2 problem

clear all
close all

maxdate = 2011; % used for plots sets max axis. Increase as you update data
%enddate = 20061229; % use 0 for full sample. If you use a date, make sure it exists. 
enddate = 0; % use 0 for full sample. If you use a date, make sure it exists. 


x = load('..\data\famabliss.txt'); % quarterly updates
beg = find(x(:,1)==19640131);
if enddate == 0; 
    x = x(beg:end,:);
else 
    en = find(x(:,1)==enddate);
    x = x(beg:en,:);
end
fbdate = x(:,1);
T = size(x,1);
T2 = find(x(:,1)==20061229); % used below 
plotdate = (1964+1/12:1/12:1964+T/12)';
p = log(x(:,2:end)/100);
y = -p.*(ones(T,1)*[1 1/2 1/3 1/4 1/5]);
f = [y(:,1) p(:,1:end-1)-p(:,2:end)];
r = [p(13:end,1:end-1)]-[p(1:end-12,2:end)]; 
rx = r - y(1:end-12,1)*ones(1,4); 

figure
subplot(3,1,1); 
plot(plotdate,100*y);
axis([1964 maxdate -inf inf]); 
title('yields')
subplot(3,1,2); 
plot(plotdate,100*f); 
title('forwards'); 
axis([1964 maxdate -inf inf]); 
subplot(3,1,3); 
plot(plotdate(13:T),100*rx); 
title('excess returns'); 
axis([1964 maxdate -inf inf]); 
print -depsc2 yields_plot.eps

% eigenvalue decompositions

[Qy,Ly] = eig(cov(y));
QLy = Qy*Ly.^0.5;
y_factors = y*Qy(:,end-2:end);  % save for later


[Qd,Ld] = eig(cov(diff(y))); 
QLd = Qd*Ld.^0.5;

[Qf,Lf] = eig(cov(f));
QLf = Qf*Lf.^0.5;
f_factors = f*Qf(:,end-2:end);  % save for later

[Qr,Lr] = eig(cov(rx)); 
QLr = Qr*Lr.^0.5;

fprintf('variance of factors -- eigenvalues  \n'); 
fprintf('  var y %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*diag(Ly)/sum(diag(Ly))); 
fprintf('   sd y %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*(diag(Ly).^0.5)); 
fprintf(' var Dy %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*diag(Ld)/sum(diag(Ld))); 
fprintf('  sd Dy %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*(diag(Ld).^0.5)); 
fprintf(' var f  %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*diag(Lf)/sum(diag(Lf))); 
fprintf('  sd f  %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*(diag(Lf).^0.5)); 
fprintf(' var rx %8.2f %8.2f %8.2f %8.2f \n',100*diag(Lr)/sum(diag(Lr))); 
fprintf('  sd rx %8.2f %8.2f %8.2f %8.2f \n',100*(diag(Lr).^0.5)); 

figure
subplot(2,2,1); 
plot(QLy(:,end:-1:1),'-v','LineWidth',2); % plotting columns backwards means the colors stay the same across y and rx which has fewer
hold on
plot((1:5)',0*(1:5)','--k')
title('y'); 

subplot(2,2,2); 
plot(QLd(:,end:-1:1),'-v','LineWidth',2); 
hold on
plot((1:5)',0*(1:5)','--k')
title('\Deltay'); 

subplot(2,2,3); 
plot(QLf(:,end:-1:1),'-v','LineWidth',2); 
hold on
plot((1:5)',0*(1:5)','--k')
title('f'); 

subplot(2,2,4); 
plot(QLr(:,end:-1:1),'-v','LineWidth',2); 
hold on
plot((1:4)',0*(1:4)','--k')
title('rx'); 

print -depsc2 eig1.eps

figure
subplot(2,2,1); 
plot(Qy(:,end:-1:end-2),'-v','LineWidth',2);
hold on
plot((1:5)',0*(1:5)','--k')
title('y'); 

subplot(2,2,2); 
plot(Qd(:,end:-1:end-2),'-v','LineWidth',2); 
hold on
plot((1:5)',0*(1:5)','--k')
title('\Deltay'); 

subplot(2,2,3); 
plot(Qf(:,end:-1:end-2),'-v','LineWidth',2); 
hold on
plot((1:5)',0*(1:5)','--k')
title('f'); 

subplot(2,2,4); 
plot(Qr(:,end:-1:end-2),'-v','LineWidth',2); 
hold on
plot((1:4)',0*(1:4)','--k')
title('rx'); 
print -depsc2 eig2.eps

% three factor model approximations


if 1
    ax1 = [ 1998 maxdate 0 7];
    axr = [1998 maxdate -10 12];
else
    ax1 = [ 1964 maxdate -inf inf];
    axr = [ 1964 maxdate -inf inf];
end

figure;
subplot(3,1,1); 
plot(plotdate,100*y,'-b');
hold on
plot(plotdate,100*y*Qy(:,end-2:end)*Qy(:,end-2:end)','-r')
title('y'); 
axis(ax1); 
title('yields')
subplot(3,1,2); 
plot(plotdate,100*f,'-b')
hold on
plot(plotdate,100*f*Qf(:,end-3:end)*Qf(:,end-3:end)','-r'); 
title('forwards'); 
axis(ax1); 
subplot(3,1,3);
plot(plotdate(13:T),100*rx,'-b');
hold on
plot(plotdate(13:T),100*rx*Qr(:,end-2:end)*Qr(:,end-2:end)','-r'); 
title('excess returns'); 
axis(axr); 

print -depsc2 threef1.eps

% --------------
% now fama bliss
% --------------

fprintf('fama-bliss \n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'n','a','b','t','R2')
for n = 2:5;  
    lhv = rx(:,n-1); 
    rhv = [ones(T-12,1) f(1:T-12,n)-y(1:T-12,1)]; 
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n', n,100*bv(1),bv(2),bv(2)/sebv(2),R2v);
end;

% find the complementary regressions
% rx(n)_t+1 = p(n-1)_t+1 - p(n)_t =  
% 

fprintf('fama-bliss complemetary  \n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'n','b','1-b','t','R2')
for n = 2:5;  
    lhv = (n-1)*(y(13:T,n-1)-y(1:T-12,n-1)); 
    rhv = [ones(T-12,1) f(1:T-12,n)-y(1:T-12,1)]; 
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n', n,bv(2),1-bv(2),bv(2)/sebv(2),R2v);
end;

% long run interest rate

fprintf('fama-bliss long-run interest rate  \n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'n','a','b','t','R2')
for n = 1:4;  
    lhv = (y(1+12*n:T,1)-y(1:T-12*n,1)); 
    rhv = [ones(T-12*n,1) f(1:T-12*n,n+1)-y(1:T-12*n,1)]; 
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n', n,bv(1),bv(2),bv(2)/sebv(2),R2v);
end;


% ----------------------------
% what CP should have done
% ---------------------------

lhv = rx; 
rhv = [ones(T-12,1) f(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
beta = bv; 
fprintf(' CP unconstrained regressions \n'); 
fprintf('       %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n','c','y1','f2','f3','f4','f5','R2'); 
fprintf(' rx(2) %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',bv(:,1),R2v(1)); 
fprintf(' rx(3) %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',bv(:,2),R2v(2)); 
fprintf(' rx(4) %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',bv(:,3),R2v(3)); 
fprintf(' rx(5) %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',bv(:,4),R2v(4));
fit = rhv*bv;


[Q,L] = eig(cov(fit)); 
figure
subplot(1,2,1)
plot(Q(:,end:-1:end-2),'-v','LineWidth',2);
title('Q'); 
subplot(1,2,2);
Q2 = 100*Q*L.^0.5;
plot(Q2(:,end:-1:1),'-v','LineWidth',2);
title('Q*L^{0.5}'); 
print -depsc2 eigs-er.eps;

fprintf('fraction of variance, eigenvalues of expected returns\n'); 
fprintf('%8.2f %8.2f %8.2f %8.2f\n', diag(L)/sum(diag(L))*100)

%er = Q*factors
%factors = Q'*Er
%factor1 = Q(:,end)'Er
%Er(model) = Q(:end)*Q(:,end)'Er
%Er(model) = Q(:,end)*Q(:,end)'*(a + beta*f)

x = fit*Q(:,end);

figure; 
subplot(2,1,1)
plot(plotdate(1:T-12),100*fit)
title('Expected excess bond returns -- unconstrained'); 
axis([1964 maxdate -inf inf])
subplot(2,1,2)
plot(plotdate(1:T-12),100*fit*Q(:,end)*Q(:,end)');
title('Constrained'); 
axis([1964 maxdate -inf inf])
print -depsc2 CP1.eps

beta_constrained = (Q(:,end)*Q(:,end)'*beta(2:end,:)')'; 
figure; 
subplot(2,1,1); 
plot(beta(2:end,:),'-v','LineWidth',2); 
title('Unconstrained')
subplot(2,1,2);
plot(beta_constrained,'-v','LineWidth',2);
title('Constrained');
if enddate == 0
    print -depsc2 CP2.eps
else
    print -depsc2 CP2_to_2006.eps
end

% portfolio
lhv = rx*Q(:,end);
rhv = [ones(T-12,1) f(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);

figure; 
subplot(2,1,1); 
plot((Q(:,end)'*beta(2:end,:)')','-v','LineWidth',2); 
title('Q(:,end)'' \beta'); 
subplot(2,1,2); 
plot(bv(2:end),'-v','LineWidth',2); 
title('regression')

% factor model

lhv = f(1:T-12,:);
rhv = [ones(T-12,1) x]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
figure
plot(bv(2:end,:),'-v','LineWidth',2)
title('Response of f to x'); 
xlabel('Maturity')
ylabel('Response to x'); 

% now, do we need the new factor? 
% 1. forecast average (across maturity) returns


lhv = mean(rx')'; 
   
rhv = [ones(T-12,1) f_factors(1:T-12,end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on 3 forward factors, full sample\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','curve','slope','level','R2')
    fprintf('b %8.2f %8.2s %8.2s %8.2f %8.2f  \n',bv(1), ' ',' ',bv(2),R2v);
    fprintf('t %8.2f %8.2s %8.2s %8.2f        \n',bv(1)./sebv(1),' ',' ',bv(2)/sebv(2));


rhv = [ones(T-12,1) f_factors(1:T-12,end-1:end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2s %8.2f %8.2f %8.2f  \n',bv(1), ' ',bv(2:end),R2v);
    fprintf('t %8.2f %8.2s %8.2f %8.2f        \n',bv(1)/sebv(1),' ',bv(2:end)./sebv(2:end));
    
rhv = [ones(T-12,1) f_factors(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

rhv = [ones(T-12,1) x ]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x \n') 
fprintf('  %8.8s %8.8s %8.8s   \n', 'const','x','R2')
    fprintf('b %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f         \n',bv./sebv);

rhv = [ones(T-12,1) x f_factors(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x and 3 forward factors\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s  \n','const','x','curve','slope','level','R2')
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

% ------------- ssample to 2006 
    
lhv = mean(rx(1:T2-12,:)')';
  
rhv = [ones(T2-12,1) f_factors(1:T2-12,end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on 3 forward factors, sample to 2006\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','curve','slope','level','R2')
    fprintf('b %8.2f %8.2s %8.2s %8.2f %8.2f  \n',bv(1), ' ',' ',bv(2),R2v);
    fprintf('t %8.2f %8.2s %8.2s %8.2f        \n',bv(1)./sebv(1),' ',' ',bv(2)/sebv(2));


rhv = [ones(T2-12,1) f_factors(1:T2-12,end-1:end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2s %8.2f %8.2f %8.2f  \n',bv(1), ' ',bv(2:end),R2v);
    fprintf('t %8.2f %8.2s %8.2f %8.2f        \n',bv(1)/sebv(1),' ',bv(2:end)./sebv(2:end));
    
rhv = [ones(T2-12,1) f_factors(1:T2-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

rhv = [ones(T2-12,1) x(1:T2-12) ]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x \n') 
fprintf('  %8.8s %8.8s %8.8s   \n', 'const','x','R2')
    fprintf('b %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f         \n',bv./sebv);

rhv = [ones(T2-12,1) x(1:T2-12) f_factors(1:T2-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x and 3 forward factors\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s  \n','const','x','curve','slope','level','R2')
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);


% yield factors 
    
lhv = mean(rx')'; 
   
rhv = [ones(T-12,1) y_factors(1:T-12,end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on 3 yield factors, full sample\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','curve','slope','level','R2')
    fprintf('b %8.2f %8.2s %8.2s %8.2f %8.2f  \n',bv(1), ' ',' ',bv(2),R2v);
    fprintf('t %8.2f %8.2s %8.2s %8.2f        \n',bv(1)./sebv(1),' ',' ',bv(2)/sebv(2));


rhv = [ones(T-12,1) y_factors(1:T-12,end-1:end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2s %8.2f %8.2f %8.2f  \n',bv(1), ' ',bv(2:end),R2v);
    fprintf('t %8.2f %8.2s %8.2f %8.2f        \n',bv(1)/sebv(1),' ',bv(2:end)./sebv(2:end));
    
rhv = [ones(T-12,1) y_factors(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

rhv = [ones(T-12,1) x ]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x \n') 
fprintf('  %8.8s %8.8s %8.8s   \n', 'const','x','R2')
    fprintf('b %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f         \n',bv./sebv);

rhv = [ones(T-12,1) x y_factors(1:T-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x and 3 yield factors\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s  \n','const','x','curve','slope','level','R2')
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

% ------------- yield factors, sample to 2006 
    
lhv = mean(rx(1:T2-12,:)')';
  
rhv = [ones(T2-12,1) y_factors(1:T2-12,end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on 3 yield factors, sample to 2006\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','curve','slope','level','R2')
    fprintf('b %8.2f %8.2s %8.2s %8.2f %8.2f  \n',bv(1), ' ',' ',bv(2),R2v);
    fprintf('t %8.2f %8.2s %8.2s %8.2f        \n',bv(1)./sebv(1),' ',' ',bv(2)/sebv(2));


rhv = [ones(T2-12,1) y_factors(1:T2-12,end-1:end)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2s %8.2f %8.2f %8.2f  \n',bv(1), ' ',bv(2:end),R2v);
    fprintf('t %8.2f %8.2s %8.2f %8.2f        \n',bv(1)/sebv(1),' ',bv(2:end)./sebv(2:end));
    
rhv = [ones(T2-12,1) y_factors(1:T2-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);

rhv = [ones(T2-12,1) x(1:T2-12) ]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x \n') 
fprintf('  %8.8s %8.8s %8.8s   \n', 'const','x','R2')
    fprintf('b %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f         \n',bv./sebv);

rhv = [ones(T2-12,1) x(1:T2-12) y_factors(1:T2-12,:)]; 
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
fprintf('regression of average returns on x and 3 yield factors\n') 
fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s  \n','const','x','curve','slope','level','R2')
    fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',bv,R2v);
    fprintf('t %8.2f %8.2f %8.2f %8.2f %8.2f        \n',bv./sebv);


    
% 2. same excercise, here forecast individual bond returns
if 0
    lhv = rx;
    rhv = [ones(T-12,1) x f_factors(1:T-12,:)];
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('regression of returns on x and 3 forward factors\n')
    fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'n','const','x','curve','slope','level','R2')
    for n = 2:5;
        fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',n,bv(:,n-1),R2v(n-1));
        fprintf('t %8.8s %8.2f %8.2f %8.2f %8.2f %8.2f        \n',' ',bv(:,n-1)./sebv(:,n-1));
    end;
    
    lhv = rx;
    rhv = [ones(T-12,1) x ];
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('regression of returns on x \n')
    fprintf('  %8.8s %8.8s %8.8s %8.8s \n', 'n','const','x','R2')
    for n = 2:5;
        fprintf('b %8.2f %8.2f %8.2f %8.2f \n',n,bv(:,n-1),R2v(n-1));
        fprintf('t %8.8s %8.2f %8.2f        \n',' ',bv(:,n-1)./sebv(:,n-1));
    end;
    
    lhv = rx;
    rhv = [ones(T-12,1) f_factors(1:T-12,:)];
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,12,0);
    fprintf('regression of returns on 3 forward factors\n')
    fprintf('  %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'n','const','curve','slope','level','R2')
    for n = 2:5;
        fprintf('b %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',n,bv(:,n-1),R2v(n-1));
        fprintf('t %8.8s %8.2f %8.2f %8.2f %8.2f        \n',' ',bv(:,n-1)./sebv(:,n-1));
    end;
end