% JCPLOTS
% makes plots of daily interest rate data
% defines shocks by change in interest rates along with target change
% starts with data loading copied from policy.m

clear; 
close all;

pub = 0; % 0 to make pretty colored graphs for appendix 
if pub; % B& W for publication 
    set(0,'DefaultAxesColorOrder',[0,0,0])
end; 

%--------------target--------------------------%
% target 
load ..\data\target.dat;
% target starts Sept 13, 74
% stop in 20011130
target=target(1:end-20,:); 

tgt_day  = target(:,2);
tgt_mo = target(:,1);
tgt_yr = target(:,3);
target = target(:,4);


%--------------treasuries----------------------%
% 19620201 is the start date
% 20011130 is the end date
load ..\data\dgs1.dat;
load ..\data\dgs3.dat;
load ..\data\dgs5.dat;
load ..\data\dgs10.dat;
load ..\data\dtb3.dat; 
load ..\data\dtb6.dat;

yields=[dgs10(:,1) dtb3(:,2) dtb6(:,2) dgs1(:,2) dgs3(:,2) dgs5(:,2) dgs10(:,2)];
% take average of past observations for missing ones
% different than Monika since there are serveral places with 
% multiple zeros. 

T=length(yields);
M=size(yields,2);
for t=1:T;
    for m=2:M
        if yields(t,m)< 1E-4;
            yields(t,m)=yields(t-1,m);
        end
    end
end

dates=dgs10(:,1);
yield_yr = floor(dates/10000);
yield_mo = dates-10000*yield_yr;
yield_mo = floor(yield_mo/100);
yield_day = dates-10000*yield_yr-100*yield_mo;

yield_date_graph = yield_yr + (yield_mo-1)/12 + (yield_day-1)/365;
%date for graphing. It would be better to recognize different days in each mol.


% --------------- eurodollars ---------------- %
% starts o1041971
% ends 11302001

load ..\data\ed1m.dat;
load ..\data\ed3m.dat;
load ..\data\ed6m.dat;

euro_yields = [ed1m(:,1) ed1m(:,2) ed3m(:,2) ed6m(:,2)];
T=length(euro_yields);
M=size(euro_yields,2);
for t=1:T;
    for m=2:M
        if euro_yields(t,m)< 1E-4;
            euro_yields(t,m)=euro_yields(t-1,m);
        end
    end
end

% shift euro yields back a day to account for 9:30 AM reporting 
euro_yields = [euro_yields(:,1) [euro_yields(2:end,2:end); euro_yields(end,2:end)]];

%--------------fed funds rate ------------------%
% fed funds rate starts in 1954 ends 2001 12 09
% start in 19620201 like the others
load dff.dat;
dff=dff(2773:end,:);

% fed funds rate has weekends, 
% take those out
ff=dff(1:2,:); 
for i=5:7:length(dff);
    j=(i:i+4)';
    ff=[ff;dff(j,:)];
end;
% stop in 20011130
ff=ff(1:end-5,:);


% ------------------------------------------------%
% target index 1         to 1309
%        13 Sept 74         19 Sept 79
% assume that we therefore can construct 
% policy shocks from Sept 74 to Sept 79

% target index 2449      to end
%        1 February 1984 
% have data from 3292       to 3292+1308
%                3292+2448  to end
% we certainly have policy shocks from 
% February 1984 onwards

dates4=dates(3292+2448:end,:);
yields4 = yields(3292+2448:end,:);
target4 = target(2449:end,:);
yield_date_graph4 = yield_date_graph(3292+2448:end,:);
yield_yr4 = yield_yr(3292+2448:end,:);
yield_mo4 = yield_mo(3292+2448:end,:);
yield_day4 = yield_day(3292+2448:end,:);

euro_yields4 = euro_yields(3413:end,:);
ff4 = ff(3292+2448:end,:);

% select those days at which target moves happen
changes=diff(target4); %  first difference of target
ind=find(changes); %  indices of when first difference is nonzero
ind = ind(2:end);  % ignore huge 1.5 target change in 198401. 
ind=ind+1;         %  add 1 because first differ

% make plots

% -------------------------
% generic plot of all rates for our information, zooming. Not for paper. 
% -------------------------
if 0; 
figure;
plot(yield_date_graph4(1:2328), [target4(1:2328) yields4(1:2328,3:end)]);
legend('target','6 mo','1 yr','3 yr','5 yr','10 yr');

figure;
plot(yield_date_graph4(2328:end), [target4(2328:end) yields4(2328:end,3:end)]);
legend('target','6 mo','1 yr','3 yr','5 yr','10 yr');
end; 

% ------------------------
% graph of 2001 for paper 
% ------------------------

disp('dates of changes in 2001');
disp([yield_yr4(ind(end-9:end)) yield_mo4(ind(end-9:end)) yield_day4(ind(end-9:end))]);

disp('interest rates around target changes -- use for text accompanying figure'); 
for indx = 9:-1:0;
   disp('month') 
   disp(yield_mo4(ind(end-indx))); 
   disp(' date euros 1m yields 1 3 5 10');
   disp([  yield_day4(ind(end-indx)-5:ind(end-indx)+4) ...    
           target4(ind(end-indx)-5:ind(end-indx)+4) ...
           euro_yields4(ind(end-indx)-5:ind(end-indx)+4,2)...
           yields4(ind(end-indx)-5:ind(end-indx)+4,4:end)]);
end;

if not(pub); % color and appendix version 

 figure;
 subplot(2,1,1);
 plot(yield_date_graph4(4375:end), [target4(4375:end) euro_yields4(4375:end,2)] );
 legend('target','1 mo Euro');
 axis([yield_date_graph4(4375) yield_date_graph4(end) 2 7 ]);
 text(yield_date_graph4(ind(end-9)+2),target4(ind(end-9))+0.25,'Jan 3*');
 text(yield_date_graph4(ind(end-8)+2),target4(ind(end-8))+0.25,'Jan 31');
 text(yield_date_graph4(ind(end-7)+2),target4(ind(end-7))+0.25,'Mar 20');
 text(yield_date_graph4(ind(end-6)+2),target4(ind(end-6))+0.25,'Apr 18*');
 text(yield_date_graph4(ind(end-5)+2),target4(ind(end-5))+0.25,'May 15');
 text(yield_date_graph4(ind(end-4)+2),target4(ind(end-4))+0.25,'Jun 27');
 text(yield_date_graph4(ind(end-3)+2),target4(ind(end-3))+0.25,'Aug 21');
 text(yield_date_graph4(ind(end-2)+2),target4(ind(end-2))+0.25,'Sept 17*');
 text(yield_date_graph4(ind(end-1)+2),target4(ind(end-1))+0.25,'Oct 2');
 text(yield_date_graph4(ind(end-0)+2),target4(ind(end-0))+0.25,'Nov 6');

 subplot(2,1,2);
 plot(yield_date_graph4(4375:end), [target4(4375:end) yields4(4375:end,4:end)]);
 legend('target','1 yr','3 yr','5 yr','10 yr');
 axis([yield_date_graph4(4375) yield_date_graph4(end) 2 7 ]);
 text(yield_date_graph4(ind(end-9)+2),target4(ind(end-9))+0.25,'Jan 3*');
 text(yield_date_graph4(ind(end-8)+2),target4(ind(end-8))+0.25,'Jan 31');
 text(yield_date_graph4(ind(end-7)+2),target4(ind(end-7))+0.25,'Mar 20');
 text(yield_date_graph4(ind(end-6)+2),target4(ind(end-6))+0.25,'Apr 18*');
 text(yield_date_graph4(ind(end-5)+2),target4(ind(end-5))+0.25,'May 15');
 text(yield_date_graph4(ind(end-4)+2),target4(ind(end-4))+0.25,'Jun 27');
 text(yield_date_graph4(ind(end-3)+2),target4(ind(end-3))+0.25,'Aug 21');
 text(yield_date_graph4(ind(end-2)+2),target4(ind(end-2))+0.25,'Sept 17*');
 text(yield_date_graph4(ind(end-1)+2),target4(ind(end-1))+0.25,'Oct 2');
 text(yield_date_graph4(ind(end-0)+2),target4(ind(end-0))+0.25,'Nov 6'); 

 print -depsc2 history1.eps;
end; 
if pub; % simpler version for published versoin
 figure;
 plot(yield_date_graph4(4375:end), [target4(4375:end) euro_yields4(4375:end,2) ] );
 axis([yield_date_graph4(4375) yield_date_graph4(end) 2 7 ]);
 set(gca,'XTick',2000+10/12:1/12:2002);
 set(gca,'XTickLabel',{'N','D','J','F','M','A','M','J','J','A','S','O','N','D'},'FontSize',16);
 set(gca,'ytick',[2 3 4 5 6 7],'FontSize',16);
 text(yield_date_graph4(ind(end-9)+2),target4(ind(end-9))+0.25,'Jan 3*','FontSize',16);
 text(yield_date_graph4(ind(end-8)+2),target4(ind(end-8))+0.25,'Jan 31','FontSize',16);
 text(yield_date_graph4(ind(end-7)+2),target4(ind(end-7))+0.25,'Mar 20','FontSize',16);
 text(yield_date_graph4(ind(end-6)+2),target4(ind(end-6))+0.25,'Apr 18*','FontSize',16);
 text(yield_date_graph4(ind(end-5)+2),target4(ind(end-5))+0.25,'May 15','FontSize',16);
 text(yield_date_graph4(ind(end-4)+2),target4(ind(end-4))+0.25,'Jun 27','FontSize',16);
 text(yield_date_graph4(ind(end-3)+2),target4(ind(end-3))+0.25,'Aug 21','FontSize',16);
 text(yield_date_graph4(ind(end-2)+2),target4(ind(end-2))+0.25,'Sep 17*','FontSize',16);
 text(yield_date_graph4(ind(end-1)+2),target4(ind(end-1))+0.25,'Oct 2','FontSize',16);
 text(yield_date_graph4(ind(end-0)+2),target4(ind(end-0))+0.25,'Nov 6','FontSize',16);


 print -deps2 hist1a.eps
 figure; 
 plot(yield_date_graph4(4375:end), [target4(4375:end) yields4(4375:end,4:end)]);
 axis([yield_date_graph4(4375) yield_date_graph4(end) 2 7 ]);
 set(gca,'XTick',2000+10/12:1/12:2002);
 set(gca,'XTickLabel',{'N','D','J','F','M','A','M','J','J','A','S','O','N','D'},'FontSize',16);
 set(gca,'ytick',[2 3 4 5 6 7],'FontSize',16);
 print -deps2 hist1b.eps;
end; 

% --------------
% plot of 1987 for paper
% -------------

if not(pub); 
 figure;
 subplot(2,1,1);
 plot(yield_date_graph4(780:1050), [target4(780:1050) euro_yields4(780:1050,2)] );
 legend('target','1 mo Euro');
 axis([yield_date_graph4(780) yield_date_graph4(1050) 5.5 10.5 ]); 

 subplot(2,1,2);
 plot(yield_date_graph4(780:1050), [target4(780:1050) yields4(780:1050,4:end)]);
 legend('target','1 yr','3 yr','5 yr','10 yr');
 axis([yield_date_graph4(780) yield_date_graph4(1050) 5.5 10.5 ]);

 print -depsc2 history3.eps;
end; 

% -----------------------------------------------------
% identify shocks by change in 1 month rate over shock. 
% ------------------------------------------------

tgtchange = changes(ind-1);  % amount of target changes 
lags = 3; % how far back to look for initial rate. ind is the 1st day of new rate
   
   %if change on date ind=4 (i.e. 3 is last day of old rate) then look 
   % at change from 1 to 5. In many episodes it moves first. 

   eurochange  = euro_yields4(ind+1,2:end)-euro_yields4(ind-lags,2:end);
   yieldchange = yields4(ind+1,2:end)-yields4(ind-lags,2:end);
   
   [bv,sebv,R2v,v] = olsgmm([eurochange yieldchange],[ones(size(tgtchange,1),1) tgtchange],0,0);
   disp('OLS regressions of yield changes on target change  -- b, t, R2');
   disp(' Euro         Treasury ');
   disp('  1mo     3mo    6 mo     3 mo       6 mo         1y        3y         5y         10y');
   disp(bv(2,:));
   disp(bv(2,:)./sebv(2,:)); 
   disp(R2v'); 
   
   [bve,sebve,R2ve,v] = olsgmm([eurochange(:,2:end) yieldchange],[ones(size(eurochange,1),1) eurochange(:,1)],0,0);
   disp('OLS regressions of yield changes on 1 mo Euro change around target -- b, t, R2');
   disp(' Euro         Treasury ');
   disp('    3mo    6 mo     3 mo       6 mo         1y        3y         5y         10y');
   disp(bve(2,:));
   disp(bve(2,:)./sebve(2,:)); 
   disp(R2ve'); 
     
    % time series plot of target change, one month rate, and shocks. 
    eurochange4 = zeros(size(euro_yields4,1),size(euro_yields4,2)-1);
    eurochange4(ind,:)=eurochange;
    
    figure; 
    plot(yield_date_graph4, [target4 euro_yields4(:,2) eurochange4(:,1)+4]);
    
    % scatter plots of change i on 1 mo target  
   if 0 ; % commented out. Really a silly waste of space -- r2 caputres it 
    allsame = 1; % all on same axis -- may miss some outliers
    allax = [-1 1 -0.75 0.75];  
 	figure ; 
	 subplot(2,3,1);
    plot(tgtchange,eurochange(:,2),'.');
    if allsame; 
       axis(allax);
    end; 
    hold on;
    plot(tgtchange,bv(1,1)+bv(2,1)*tgtchange,'-k');
	 title('1 mo euro');
	 xlabel('target');
	 ylabel('yield');
 	
	 subplot(2,3,2);
    plot(tgtchange,eurochange(:,3),'.');
       if allsame; 
       axis(allax);
    end; 
    hold on;
    plot(tgtchange,bv(1,2)+bv(2,2)*tgtchange,'-k');
	 title('3 mo euro');
	 xlabel('target');
	 ylabel('yield');  

  	subplot(2,3,3); 
  	plot(tgtchange,yieldchange(:,2),'.');
      if allsame; 
      axis(allax);
   end; 
    hold on;
    plot(tgtchange,bv(1,4)+bv(2,4)*tgtchange,'-k');
    title('3 mo Tresury');
   	xlabel('target'); 
	 ylabel('yield');

	subplot(2,3,4);
	plot(tgtchange,yieldchange(:,4),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(tgtchange,bv(1,6)+bv(2,6)*tgtchange,'-k');
   title('1 y Treaury');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,5);
	plot(tgtchange,yieldchange(:,5),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(tgtchange,bv(1,7)+bv(2,7)*tgtchange,'-k');
   title('3 y Treaury');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,6);
	plot(tgtchange,yieldchange(:,6),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(tgtchange,bv(1,8)+bv(2,8)*tgtchange,'-k');
   title('5 y Treaury');
	xlabel('target');
	ylabel('yield');
   
   print -depsc2 changes1.eps;
   
   figure;
	subplot(2,3,1);
	plot(eurochange(:,2),eurochange(:,2),'.');
   if allsame; 
      axis(allax);
   end; 
   title('1 mo euro');
	xlabel('target');
	ylabel('yield');
   
	subplot(2,3,2);
	plot(eurochange(:,2),eurochange(:,3),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(eurochange(:,2),bve(1,1)+bve(2,1)*eurochange(:,2),'-k');
   title('3 mo euro');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,3);
	plot(eurochange(:,2),yieldchange(:,2),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(eurochange(:,2),bve(1,3)+bve(2,3)*eurochange(:,2),'-k');
   title('3 mo Tresury');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,4);
	plot(eurochange(:,2),yieldchange(:,4),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(eurochange(:,2),bve(1,5)+bve(2,5)*eurochange(:,2),'-k');
   title('1 y Treaury');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,5);
	plot(eurochange(:,2),yieldchange(:,5),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(eurochange(:,2),bve(1,6)+bve(2,6)*eurochange(:,2),'-k');
   title('3 y Treaury');
	xlabel('target');
	ylabel('yield');

	subplot(2,3,6);
	plot(eurochange(:,2),yieldchange(:,6),'.');
   if allsame; 
      axis(allax);
   end; 
   hold on;
   plot(eurochange(:,2),bve(1,7)+bve(2,7)*eurochange(:,2),'-k');
   title('5 y Treaury');
	xlabel('target');
   ylabel('yield');
   
   print -depsc2 changes2.eps;
end; % ends scatterplots    
   % ---------------------
   % plot of 3 mo euro vs. 3 mo treasury -- why the slope?
   % ---------------------
   if not(pub);
    figure;
    plot(euro_yields4(:,3),[yields4(:,2) euro_yields4(:,3)]);
    title('why not a 45o line?');
end;    
   % ----------------
   % plots of some changes to see dating
   % -----------------
   
if 0; 
   for indx=size(ind,1):-1:size(ind,1)-10;
      figure;
      plot(yield_date_graph4(ind(indx)-5:ind(indx)+4),...
         [target4(ind(indx)-5:ind(indx)+4) ...
             euro_yields4(ind(indx)-5:ind(indx)+4,2:3) ...
             yields4(ind(indx)-5:ind(indx)+4,3:end)],'-+');
                 legend('target','1 mo Euro','3 mo Euro','6 mo','1 yr','3 yr','5 yr','10 yr');
   end;
   
    for indx=46:56;
      figure;
      plot(yield_date_graph4(ind(indx)-5:ind(indx)+4),...
         [target4(ind(indx)-5:ind(indx)+4) ...
             euro_yields4(ind(indx)-5:ind(indx)+4,2:3) ...
             yields4(ind(indx)-5:ind(indx)+4,3:end)],'-+');
                 legend('target','1 mo Euro','3 mo Euro','6 mo','1 yr','3 yr','5 yr','10 yr');
   end;
end;