%Plot figures in paper
clc
clear all

warning off all;
 load('thetamode');
 load('THETAMCMC');

load('YFB');

Z=YFB;
kbar=25;
yc=[3,12,24,36,48,60;];
yn=length(yc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot one step ahead fit (FIGURE 2 in paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[EE,M,MH,N,NN,R,D1,D2,D3,P,F1]=termsolve(thetamode);
DD1=D1(3:end,:);
DD2=D2(3:end,:);
DD3=[0,0;D3;];

n=length(DD1(:,1)); %dimension of observables;

R1=R(3:end,1:end);R1(4:5,:)=0;
Zadj=Z(:,3:end)-[F1(1,:);D3*F1;]*[Z(1,2:end-1);Z(1,1:end-2);];
  
[ZTILDE]=KalmanLagFIT(M,N,DD1,DD2,R1,Zadj);
ZFIT=Z(:,3:end)-ZTILDE;

date=[(1964+(3/12)):(1/12):2008];
figure(2)
for j=1:3
    subplot(3,1,j);
    if j==1;
        plot(date,Z(1,3:end),'linewidth',2,'linestyle','-','color','black');
        hold on;
        plot(date,Z(1,3:end)-ZTILDE(1,:),'linewidth',2,'linestyle','--','color','blue');        
        xlabel({'Fed Funds rate'},'fontsize',16);
    end
    if j==2;
        plot(date,Z(5,3:end),'linewidth',2,'linestyle','-','color','black');
        hold on;
        plot(date,Z(5,3:end)-ZTILDE(5,:),'linewidth',2,'linestyle','--','color','blue');     
        xlabel({'3 year yield'},'fontsize',16);
    end
    if j==3;
        plot(date,Z(7,3:end),'linewidth',2,'linestyle','-','color','black');
        hold on;
        plot(date,Z(7,3:end)-ZTILDE(7,:),'linewidth',2,'linestyle','--','color','blue');     
        xlabel({'5 year yield'},'fontsize',16);
    end
    
    
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct state distribution and distributions of functions of the state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=100; % Number of draws; 100 is usually more than enough
ra = max(size(THETAMCMC));
f = ceil(ra.*rand(S,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Use this to plot distributions conditional on mode
% x=thetamode;
% [EE,M,MH,N,NN,R,D1,D2,D3,P,F1]=termsolve(thetamode);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for s=1:S;
    s
    x=THETAMCMC(:,f(s,1));
    THETA(:,s)=x;
    [EE,M,MH,N,NN,R,D1,D2,D3,P,F1]=termsolve(x);  
    
    Zadj=Z(:,3:end)-[F1(1,:);D3*F1;]*[Z(1,2:end-1);Z(1,1:end-2);];
    
    [Xss]=KalmanSimSmooth(M,N,D1(3:end,:),D2(3:end,:),R(3:end,:),Zadj);
    
    SV=dlyap(M,N*N');
    STATE(:,:,s)=Xss;
    E=[ones(1,2),zeros(1,kbar*2-2);];                                   %selector vector summing first two elements
    H=[zeros(kbar*2,2),eye(kbar*2);];H=H(1:kbar*2,1:kbar*2);        %defines "+1 order" operator Q
     
    
    MH=M*H;
    b=zeros(yn+1,kbar*2);
    bb=zeros(yn+1,kbar*2);
    b(1,:)=E;
    bb(1,:)=E;
    B=zeros(yn+1,kbar*2);
    
    for j=2:60;
        b(j,:)=E*MH^(j-1);
        B(j,:)=sum(b(1:j,:));
        bb(j,:)=-(1/j)*sum(b(1:j,:));
    end;
    bnn=b;
    bbnn=bb;
    B12=b(12,:);
    B36=b(36,:);
    B60=b(60,:);
    
    
    
    b=zeros(yn+1,kbar*2);
    bb=zeros(yn+1,kbar*2);
    b(1,:)=E;
    bb(1,:)=E;
    
    for j=2:60;
        b(j,:)=E*M^(j-1)*H;
        bb(j,:)=-(1/j)*sum(b(1:j,:));
    end;
    beh=b;
    bbeh=bb;
    BEH12=b(12,:);
    BEH36=b(36,:);
    BEH60=b(60,:);
    
    fwdvar=bnn*SV*bnn';
    specvar=(beh-bnn)*SV*(beh-bnn)';
    yldvar=bbnn*SV*bbnn';
    specvaryld=(bbeh-bbnn)*SV*(bbeh-bbnn)';    
    ratvarfwr=((diag(specvar)).^.5)./((diag(fwdvar)).^.5);
    ratvaryld=((diag(specvaryld)).^.5)./((diag(yldvar)).^.5);
    
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute dispersion of state estimate and short rate forecasts across traders
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [K,P]= KalmanLagSS(M,N,D1,D2,R);
    LOM=(eye(kbar*2)-K*D1)*M;
    ISV=dlyap(LOM,K(:,1:2)*R(1:2,:)*R(1:2,:)'*K(:,1:2)');
    for j=2:60;
        b(j,:)=E*M^(j-1);
        bb(j,:)=-(1/j)*sum(b(1:j,:));
    end;
    
    DISP(:,:,s)=b(:,:)*ISV*b(:,:)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
 
    SPEC12(:,:,s)=(B12-BEH12)*Xss;
    SPEC36(:,:,s)=(B36-BEH36)*Xss;
    SPEC60(:,:,s)=(B60-BEH60)*Xss;
    VARFRACFWR(:,:,s)=ratvarfwr;
    VARFRACYLD(:,:,s)=ratvaryld;
    
    
    LOAD121(:,s)=D1(5,1:2:50);
    LOAD122(:,s)=D1(5,2:2:50);
    LOAD361(:,s)=D1(7,1:2:50);
    LOAD362(:,s)=D1(7,2:2:50);
    LOAD601(:,s)=D1(9,1:2:50);
    LOAD602(:,s)=D1(9,2:2:50);
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct posterior probability intervals for theta
% Used in Table 1 of paper
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lower=ceil(S*.025);
median=ceil(S*.5);
upper=ceil(S*.975);
THETAsort=sort(THETA,2);
THETAlow=THETAsort(:,lower);
THETAmed=THETAsort(:,median);
THETAupp=THETAsort(:,upper);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1)
% STATEsort=sort(STATE,3);
% for j=1:10;
%     subplot(5,2,j);
%     plot(STATEsort(j,1:end-1,lower),':');
%     hold on
%     plot(STATEsort(j,1:end-1,median),'k');
%     hold on
%     plot(STATEsort(j,1:end-1,upper),':');
% end
date=[(1964+(3/12)):(1/12):2008];
SPEC12sort=sort(SPEC12,3);
SPEC36sort=sort(SPEC36,3);
SPEC60sort=sort(SPEC60,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Historical Speculation (FIGURE 3 in paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)
subplot(2,1,1);
plot(date,SPEC12sort(1,1:end,lower),'linestyle',':','linewidth',1,'color','black');
hold on
plot(date,SPEC12sort(1,1:end,median),'linestyle','-','linewidth',2,'color','black');
hold on
plot(date,SPEC12sort(1,1:end,upper),'linestyle',':','linewidth',1,'color','black');

subplot(2,1,2);
plot(date,SPEC60sort(1,1:end,lower),'linestyle',':','linewidth',1,'color','black');
hold on
plot(date,SPEC60sort(1,1:end,median),'linestyle','-','linewidth',2,'color','black');
hold on
plot(date,SPEC60sort(1,1:end,upper),'linestyle',':','linewidth',1,'color','black');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot dispersion of first order expectations of short rate across horizons
% (Used in Table 2 in paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7)
DISPsort=sort(DISP,3);
for j=1:60;
    DispLow(j)=(DISPsort(j,j,lower));
    DispMed(j)=(DISPsort(j,j,median));
    DispUpp(j)=(DISPsort(j,j,upper));
end
plot(DispLow.^.5,'linestyle','--','linewidth',1,'color','black');
hold on
plot(DispMed.^.5,'linestyle','-','linewidth',2,'color','black');
hold on
plot(DispUpp.^.5,'linestyle','--','linewidth',1,'color','black');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot relative standard deviation of speculative term in forward rates and
% yields 
% (FIGURE 4 in paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4)
subplot(1,2,1);
VARFRACFsort=sort(VARFRACFWR,3);
plot(VARFRACFsort(:,1,upper),'linestyle','--','linewidth',1,'color','black');
hold on;
plot(VARFRACFsort(:,1,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(VARFRACFsort(:,1,lower),'linestyle','--','linewidth',1,'color','black');


figure(4)
subplot(1,2,2);
VARFRACYsort=sort(VARFRACYLD,3);
plot(VARFRACYsort(:,1,upper),'linestyle','--','linewidth',1,'color','black');
hold on;
plot(VARFRACYsort(:,1,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(VARFRACYsort(:,1,lower),'linestyle','--','linewidth',1,'color','black');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot distribution of loadings onto factors and higher order expectations
% of factor
% (Used in Figure 7 in the Appendix of paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(7)
subplot(3,2,1);
LOAD121sort=sort(LOAD121,2);
plot(LOAD121sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD121sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD121sort(:,lower),'linestyle',':','linewidth',1,'color','black');

subplot(3,2,2);
LOAD122sort=sort(LOAD122,2);
plot(LOAD122sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD122sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD122sort(:,lower),'linestyle',':','linewidth',1,'color','black');

subplot(3,2,3);
LOAD361sort=sort(LOAD361,2);
plot(LOAD361sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD361sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD361sort(:,lower),'linestyle',':','linewidth',1,'color','black');

subplot(3,2,4);
LOAD362sort=sort(LOAD362,2);
plot(LOAD362sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD362sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD362sort(:,lower),'linestyle',':','linewidth',1,'color','black');


subplot(3,2,5);
LOAD601sort=sort(LOAD601,2);
plot(LOAD601sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD601sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD601sort(:,lower),'linestyle',':','linewidth',1,'color','black');

subplot(3,2,6);
LOAD602sort=sort(LOAD602,2);
plot(LOAD602sort(:,upper),'linestyle',':','linewidth',1,'color','black');
hold on;
plot(LOAD602sort(:,median),'linestyle','-','linewidth',2,'color','black');
hold on;
plot(LOAD602sort(:,lower),'linestyle',':','linewidth',1,'color','black');

