% problem set 10 programs Fama Bliss

close all;
clear all; 

format short; 
format compact; 
addpath('..\Data'); 
addpath('..\Programs'); 

% -------------------------------
% new interest rate graph -- what do FB and EH say about yields
% clairvoyant yield model
% -------------------------------

 y1 = [   6   6   6 6    4   2 1   1   1   2    4  6   6   6 ]';
 sk = [[  6.2 6.2 6 5.75 4.5 3 2.5 2.5 2.5 3.5  5  6.2 6.2 6.2 ]'...  
       [  6.4 6.4 6 5.5  5   4 3.5 3.5 3.5 4.5 5.5 6.4 6.4 6.4 ]']; 
 
 T=size(y1,1);
 tim = (1:T)';
 feh2 = [y1(2:end); NaN];
 feh3 = [y1(3:end); NaN; NaN];
 feh4 = [y1(4:end); NaN; NaN; NaN];
 feh = [feh2 feh3 feh4 ];
 
 figure;
 plot(tim,y1,'-vk','linewidth',2); 
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)}');
 axis([0 T 0 7]);
 grid on
 print -depsc2 expectations_u.eps
 
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 plot(tim,feh,'--o','linewidth',2); 
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T 0 7]);
 title('Expectations hypothesis')
 grid on
 print -depsc2 expectations_u_EH_answer.eps
  
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 plot(tim,sk,'--o','linewidth',2);
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T 0 7]);
 title('Data')
 grid on
 print -depsc2 expectations_u_fact_answer.eps
 
 p = -y1/100; 
 p = [p p-sk(:,1)/100];
 p = [p p(:,2)-sk(:,2)/100];
 
 r = [NaN NaN; p(2:end,1:2)-p(1:end-1,2:3)];
 rx = 100*(r - [NaN;y1(1:end-1)/100]*ones(1,2));
 
 figure;
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
  plot(tim,sk,'--o','linewidth',1); 
 plot(tim,rx,'linewidth',2);
 plot([0;tim],0*[0;tim],'-k')
 legend('y1','f2','f3','rx2','rx3')
 xlabel('Time, years')
 ylabel('100* logs');
 axis([0 T -inf inf]);
 grid on
 print -depsc2 expectations_r.eps
 
  
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 x=[tim';tim'+1;tim'+2];
 y=[y1'; sk'];
 plot(x(:,1:5),y(:,1:5),'-o','linewidth',2);
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T 0 7]);
 title('Expected path from forwards')
 grid on
 print -depsc2 expectations_interp_facts_setup.eps
 
 
  
 
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 plot([tim';tim'+1;tim'+2],[y1'; sk'],'-o','linewidth',2);
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T 0 7]);
 title('Expected path from forwards')
 grid on
 print -depsc2 expectations_interp_facts.eps
 
 
  % ---------------------------
% resume problem set
% -------------------------

global yields forwards beg T method

% -----------------------------------------
% load fama bliss data
% this is the file from CRSP completely unmodified
% ---------------------------------------------

fama_bliss_data = load('Fama_Bliss_price.txt'); 
enddate = 2015; % used to set axis on plots. Use a bit more than end of data so you can see end of data

% collect the annual yields, form prices, forwards, hpr
T = size(fama_bliss_data,1);

rawdates = fama_bliss_data(:,1); 

dates = 1952+6/12:1/12:1952+5/12+T/12; %used on plots
[yr,mo,day,date_number] = decode_date(rawdates); %proper matlab handling of dates

% yields(t), forwards(t) are yields, forwards at time t
prices = log(fama_bliss_data(:,2:end)/100);
mats = [1 2 3 4 5]'; 
yields = -prices./(ones(T,1)*mats');
forwards = prices(:,1:4)-prices(:,2:5);
fs = forwards-yields(:,1)*ones(1,4);

% hprx(t) is the holding period return over last year
hpr = prices(13:T,1:4)-prices(1:T-12,2:5);
hprx = hpr - yields(1:T-12,1)*ones(1,4);
hpr = [NaN*ones(12,4); hpr]; % pads out the initial values so it's the same length as other series       
hprx = [NaN*ones(12,4); hprx];


beg  = 140;      % set beginning date 140 = 1964

disp('first prices'); 
fprintf(' %10.0f %9.5f %9.5f %9.5f %9.5f %9.5f \n', [rawdates(beg:beg+4) prices(beg:beg+4,:)]'); 
disp('first forwards'); 
fprintf(' %10.0f %9.5f %9.5f %9.5f %9.5f \n', [rawdates(beg:beg+4) forwards(beg:beg+4,:)]'); 
disp('first r'); 
fprintf(' %10.0f %9.5f %9.5f %9.5f %9.5f \n', [rawdates(beg+12:beg+18,:) hpr(beg+12:beg+18,:)]'); 
disp('first rx'); 
fprintf(' %10.0f %9.5f %9.5f %9.5f %9.5f \n', [rawdates(beg+12:beg+18,:) hprx(beg+12:beg+18,:)]'); 

disp(''); 


% --------------------------------------------------------------%
disp('average yields'); 
disp(mean(100*yields(beg:T,:))); 
disp('average yield spread'); 
disp(mean(100*(yields(beg:T,2:end)-yields(beg:T,1)*ones(1,4)))); 
disp('average one-year excess returns');
disp(mean(100*hprx(beg+12:T,:)));
disp('t stat  one-year excess returns');
disp(mean(100*hprx(beg+12:T,:))./(std(100*hprx(beg+12:T,:))/(((T-beg-12)/12)^0.5))); % quick calculation -- use T/12 to correct for overlap
disp('std dev holding period excess returns');
disp(std(100*hprx(beg+12:T,:)));
disp('sharpe ratio');
disp(mean(100*hprx(beg+12:T,:))./std(100*hprx(beg+12:T,:)));

% ***************************
% Fama Bliss Regressions
% ***************************

disp('----------------------------------------------------------');

   fbfcst = zeros(T-beg+1,4); 
   indx = 1; 
   while indx <= 4;
      % regression of holding period returns on forward-spot spread
	  [bi,HHi,R2un(indx),R2(indx),v] = olsgmm(100*hprx(beg+12:T,indx),[ones(T-beg-12+1,1) 100*fs(beg:T-12,indx)],12,0);
         
      bFB(indx,:)=bi';
      HH(indx,1:2) = HHi';
      
      fbfcst(:,indx) = [ones(T-beg+1,1) 100*fs(beg:T,indx)]*bi;
                
      % regression of changes in short rate on forward-spot spread
      rhv = [ones(T-beg-12*indx+1,1) 100*fs(beg:T-12*indx,indx)];
      lhv = 100*yields(beg+12*indx:T,1)-100*yields(beg:T-12*indx,1);
      [bi,HHi,R2run(indx),R2r(indx),v] = olsgmm(lhv,rhv,12*indx,0);
 
      br(indx,:)=bi';
      HHr(indx,1:2) = HHi';
      
	  indx = indx+1;
  end;
  
   fprintf('Fama Bliss Regression Coefficients, %8.0i %8.0i \n',rawdates(beg), rawdates(end) );
   disp('rx(n) (t+1) = a + b [ f (n) (t) - y(1) (t)]')
   fprintf( '%8s %8s %8s %8s %8s %8s %8s \n', 'maturity' , 'a'  ,  'b' ,  's(a)' ,'s(b)'  , 'R2',  'R2unadj');
   fprintf( '%8i %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', [(2:5)' bFB HH R2' R2un']');
   disp('y(1)(t+n) -y(1)(t) = a + b [ f (n) (t) - y(1) (t)]')
   fprintf( '%8i %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', [(1:4)' br HHr R2r' R2run']');
   disp('note: standard errors use Hansen-Hodrick GMM correction for overlapping data');

   
   % plot actual and fitted for lecture 
   
figure; 
plot(dates(beg:T)',fbfcst(:,1),'-r'); 
hold on
plot(dates(beg:T-12),100*hprx(beg+12:T,1),'-b'); % this establishes colors for legend
plot(dates(beg:T)',fbfcst,'-r'); 
plot(dates(beg:T-12),100*hprx(beg+12:T,:),'-b'); 
plot([2003 2003]',[-10 12],'-k'); 
axis([1990 enddate -10 12]); 
legend('forecast','actual',3);
text(2003,-5,'CP Sample ','HorizontalAlignment','right'); 
plot([2003 2003]',[-8 12],'-k'); 
title('Fama-Bliss forecasts') 
print -depsc2 fb_expost.eps;
   
% ******************************************
% regressions of interest rate on forward rate
% *******************************************

levelfcst = zeros(T,4);

for indx = 1:4;
    % regression of interest rate on forward rates
    rhv = [ones(T-beg-12*indx+1,1) forwards(beg:T-12*indx,indx)];
    lhv = yields(beg+12*indx:T,1);
    
    [bi,sei,R2un(indx),R2(indx),v] = olsgmm(lhv,rhv,12*indx,0);
    
    blevel(indx,:)=bi';
    se(indx,1:2) = sei';
    
    levelfcst(:,indx) = [ones(T,1) forwards(:,indx)]*bi;
    
    % regression of interest rate on interet rate
    
    rhv = [ones(T-beg-12*indx+1,1) yields(beg:T-12*indx,1)];
    lhv = yields(beg+12*indx:T,1);
    
    [bi,sei,R2uny(indx),R2y(indx),v] = olsgmm(lhv,rhv,12*indx,0);
    
    blevely(indx,:)=bi';
    sey(indx,1:2) = sei';
    
    levelfcsty(:,indx) = [ones(T,1) yields(:,1)]*bi;
    
    % multiple regression of interest rate on interest rate and forward rate
    
    rhv = [ones(T-beg-12*indx+1,1) yields(beg:T-12*indx,1) forwards(beg:T-12*indx,indx)];
    lhv = yields(beg+12*indx:T,1);
    
    [bmi,semi,R2my(indx),R2my(indx),v] = olsgmm(lhv,rhv,12*indx,0);
    
    bm(indx,:)=bmi';
    sm(indx,:) = semi';
    
end;

   fprintf( '%8s %8s %8s %8s %8s %8s %8s %8s %8s \n', 'maturity' , 'a'  ,  'b' , 'c',  's(a)' ,'s(b)', 's(c)' , 'R2',  'R2unadj');
   disp('y(1)(t+n)  = a + b  f (n-1 -> n) (t) ')
   fprintf( '%8i %8.2f %8.2f %8.0s %8.2f %8.2f %8.0s %8.2f %8.2f \n', [(1:4)' blevel [0;0;0;0] se [0;0;0;0] R2' R2un']');
   disp('y(1)(t+n)  = a + b  y(1)(t) ')
   fprintf( '%8i %8.2f %8.2f %8.0s %8.2f %8.2f %8.0s %8.2f %8.2f \n', [(1:4)' blevely [0;0;0;0] sey [0;0;0;0] R2y' R2uny']');
   disp('y(1)(t+n)  = a + b  y(1)(t) + c f(n-1 ->n) (t) ')
   fprintf( '%8i %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', [(1:4)' bm sm R2y' R2uny']');

   disp('note: standard errors use Hansen-Hodrick GMM correction for overlapping data');

   figure;
        plot(dates(beg:T),levelfcst(beg:T,1),'linewidth',2);
        hold on
        plot(dates(beg:T-12*1),yields(beg+12*1:T,1),'-r','linewidth',2)
        legend('forecast','outcome');
        title('forecast y^{(1)}_{t+1} = a + b f^{(2)}_t and outcome');
        print -depsc2 yflevel.eps;
   figure;
        plot(dates(beg:T),levelfcsty(beg:T,1),'linewidth',2);
        hold on
        plot(dates(beg:T-12*1),yields(beg+12*1:T,1),'-r','linewidth',2)
        legend('forecast','outcome');
        title('forecast y^{(1)}_{t+1} = a + b y^{(1)}_t and outcome');
        print -depsc2 yylevel.eps;
   
        
% *************************
% longer period regressions  
% *************************
   
  if 0; % *******IGNORE NOT USED

  % hprx_n(t) is the holding period excess returns for holding an n period bond n-1 periods, minus the return for holding the n-1 period bond to
  % maturity. These are complementary to the yield forecast regressions. 
  hprx_2 = prices(beg+12:T,1)-prices(beg:T-12,2)+prices(beg:T-12,1);
  hprx_3 = prices(beg+24:T,1)-prices(beg:T-24,3)+prices(beg:T-24,2);
  hprx_4 = prices(beg+36:T,1)-prices(beg:T-36,4)+prices(beg:T-36,3);
  hprx_5 = prices(beg+48:T,1)-prices(beg:T-48,5)+prices(beg:T-48,4);
  
   % regression of changes in short rate on forward-spot spread
  indx = 1; 
  rhv = [ones(length((beg:T-12*indx)),1) 100*fs(beg:T-12*indx,indx)];
  lhv = 100*hprx_2;
  [b_2,HH_2,R2_2,R2r_2,v] = olsgmm(lhv,rhv,12*indx,0);
 
  indx = 2; 
  rhv = [ones(length((beg:T-12*indx)),1) 100*fs(beg:T-12*indx,indx)];
  lhv = 100*hprx_3;
  [b_3,HH_3,R2_3,R2r_3,v] = olsgmm(lhv,rhv,12,0);

  indx = 3; 
  rhv = [ones(length((beg:T-12*indx)),1) 100*fs(beg:T-12*indx,indx)];
  lhv = 100*hprx_4;
  [b_4,HH_4,R2_4,R2r_4,v] = olsgmm(lhv,rhv,12,0);

  indx = 4; 
  rhv = [ones(length((beg:T-12*indx)),1) 100*fs(beg:T-12*indx,indx)];
  lhv = 100*hprx_5;
  [b_5,HH_5,R2_5,R2r_5,v] = olsgmm(lhv,rhv,12,0);

  
  disp('Complementary regressions to long run yield forecasts'); 
  disp('long term excess return (t+n-1) = a + b [ f (n-1 -> n) (t) - y(1) (t)]')
  fprintf('%8s %8s %8s %8s %8s %8s %8s %8s %8s \n', 'n' , 'a', 'b' , 'y coeff' , 'b/(n-1)', 'sigma(a)', 'sigma(b)', 'R2' ,  'R2unadj')
  fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', 2, b_2, br(1,2), b_2(2),  HH_2', R2_2, R2r_2 );
  fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', 3, b_3, br(2,2), 1/2*b_3(2),  HH_3', R2_3, R2r_3 );
  fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', 4, b_4, br(3,2), 1/3*b_4(2),  HH_4', R2_4, R2r_4 );
  fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', 5, b_5, br(4,2), 1/4*b_5(2),  HH_5', R2_5, R2r_5 );
     
   disp('-----------------------------------------------------');
   
  end; % ends if 0 that blocks out this section 
   
   
   % ******************************
   % Cochrane-Piazzesi regression 
   % ******************************
   
FT   = [ones(T-beg-12+1,1) 100*yields(beg:T-12,1) 100*forwards(beg:T-12,:)]; % yeilds and forwards
FT_spread   = [ones(T-beg-12+1,1) 100*(forwards(beg:T-12,:)-yields(beg:T-12,1)*ones(1,4))]; % forward spreads

      
disp(' '); 
disp(' CP Rregressions') 
fprintf('Returns on all forwards, %8.0i %8.0i \n',rawdates(beg), rawdates(end) );
[betas,stbetas,R2,R2adj,v,F_trash] = olsgmm(100*hprx(beg+12:T,:),FT,12,0);  % std errors: HH with 12 lags
fit = [ones(T,1) 100*yields(:,1) 100*forwards]*betas;
[betas_s,stbetas_s,R2_s,R2adj_s,v_s,F_trash] = olsgmm(100*hprx(beg+12:T,:),FT_spread,12,0);  % std errors: HH with 12 lags
fit_s = [ones(T,1) 100*(forwards-yields(:,1)*ones(1,4))]*betas_s;

% make the table

fprintf( '      %8s %8s %8s %8s %8s %8s %8s %8s  \n',...
                'matur', 'const', 'y1', 'f(2)',  'f(3)' , 'f(4)', 'f(5)' ,'R2' );
fprintf('levels %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',[(2:5)' betas' R2]');
fprintf('spread %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n',[(2:5)' betas_s(1,:)' zeros(4,1) betas_s(2:5,:)' R2_s]');

% graph tent
figure
plot((1:5)',betas(2:end,:),'-v','linewidth', 2); 
title('tent')
print -depsc2 tent.eps

figure
plot((2:5)',betas_s(2:end,:),'-v','linewidth', 2); 
title('tent, spreads')


% make the graph 

figure; 
plot(dates(beg:T),fit(beg:T,:));
title('CP fitted values') 
axis([1990 enddate -5 7]); 
print -depsc2 cp_fitted.eps;


% ---------
% Graphs for lecture
% --------

% fitted and actual
figure; 
plot(dates(beg:T),fit(beg:T,1),'-r'); 
hold on
plot(dates(beg:T-12),100*hprx(beg+12:T,1),'-b'); 
plot(dates(beg:T),fit(beg:T,:),'-r'); 
hold on
plot(dates(beg:T-12),100*hprx(beg+12:T,:),'-b'); 
plot([2003 2003]',[-10 12],'-k'); 
axis([1990 enddate -10 12]); 
legend('forecast','actual',3);
text(2003,-5,'CP Sample ','HorizontalAlignment','right'); 
plot([2003 2003]',[-8 12],'-k'); 
title('Unrestricted forecasts using all forwards') 


% *******************
% average across maturity
% *******************

% step 1 find gammas
[betas,stbetas,R2,R2adj,v,F_trash] = olsgmm(100*mean(hprx(beg+12:T,:)')',FT,12,0);  % std errors: HH with 12 lags
gpf = [ones(T,1) 100*yields(:,1) 100*forwards]*betas;
gammas = betas;

[betas_s,stbetas_s,R2,R2adj_s,v_s,F_trash] = olsgmm(100*mean(hprx(beg+12:T,:)')',FT_spread,12,0);  % std errors: HH with 12 lags
gpf_s = [ones(T,1) 100*(forwards-yields(:,1)*ones(1,4))]*betas_s;
gammas_s = betas_s;


%step 2 find b
[b,stb,R2b, R2adjb,v,F_trash] = olsgmm(100*hprx(beg+12:T,:),gpf(beg:T-12),12,0); 
[b_s,stb_s,R2b_s, R2adjb_s,v_s,F_trash] = olsgmm(100*hprx(beg+12:T,:),gpf_s(beg:T-12),12,0); 

% make the tables

disp(' ');
disp(' Cochrane Piazzesi Single Factor Model') 
fprintf( '      %8s %8s %8s %8s %8s %8s %8s %8s  \n',...
                'matur', 'const', 'y(1)', 'f(2)',  'f(3)' , 'f(4)', 'f(5)' ,'R2' );
fprintf('level  %8.2s %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f  \n','avg', betas', R2');
fprintf('spread %8.2s %8.2f %8.2s %8.2f %8.2f %8.2f %8.2f %8.2f  \n','avg', betas_s(1,:)', ' ', betas_s(2:end,:)',R2');

fprintf( '      %8s %8s %8s %8s   \n',...
                'b(2)',  'b(3)' , 'b(4)', 'b(5)' );
fprintf(' level %8.2f %8.2f %8.2f %8.2f  \n', b');
fprintf('spread %8.2f %8.2f %8.2f %8.2f  \n', b_s');

% graph restricted tent
figure
plot((1:5)',b'*gammas(2:end)','-v','linewidth', 2); 
title('tent')
print -depsc2 tent_restricted.eps

figure; 
plot(dates(beg:T),gpf(beg:T,:)*b);
hold on
%plot(dates(beg:T),gpf(beg:T),'-r')
title('CP restricted model fitted values') 
axis([1990 enddate -5 7]); 
print -depsc2 cp_fitted_2.eps;

figure; 
plot(dates(beg:T),gpf_s(beg:T,:)*b_s);
hold on
%plot(dates(beg:T),gpf(beg:T),'-r')
title('CP restricted model fitted values -- spreads') 
axis([1990 enddate -5 7]); 
print -depsc2 cp_fitted_2_spread.eps;


figure; 
plot(dates(beg:T),gpf(beg:T,:)*b(1),'-b');
hold on
plot(dates(beg:T),fit(beg:T,1),'-r')
plot(dates(beg:T),gpf(beg:T,:)*b,'-b');
plot(dates(beg:T),fit(beg:T,:),'-r')
title('CP restricted model fitted values') 
axis([2000 2006 -3 7]); 
legend('restricted','unrestricted')
print -depsc2 cp_fitted_3.eps;

% --------------------------
% CP eigenvalue decompositions
% ------------------

[Q,L] = eig(cov(100*hprx(beg:T,:))); 
figure
plot(Q,'-v','linewidth',2)
title('Loadings for actual excess returns'); 
disp('return factor standard deviations and percent of variance') 
disp(diag(L)'.^0.5);
disp(100*diag(L)'/sum(diag(L)));
print -depsc2 return_eig.eps;

[Q,L] = eig(cov(fit(beg:T,:))); 
figure
plot(Q,'-v','linewidth',2)
title('Loadings for expected excess returns'); 
disp('expected return factor standard deviations and percent of variance') 
disp(diag(L)'.^0.5);
disp(100*diag(L)'/sum(diag(L)));
print -depsc2 expected_return_eig.eps;

q = Q(:,end); 
factor = fit*q;
etrfit = factor*q'; 

% copied plots from above

figure; 
plot(dates(beg:T),etrfit(beg:T,:));
hold on
%plot(dates(beg:T),gpf(beg:T),'-r')
title('Eigenvalue restricted model fitted values') 
axis([1990 enddate -5 7]); 
print -depsc2 cp_fitted_2_eig.eps;


figure; 
plot(dates(beg:T),gpf(beg:T,:)*b(1),'-b');
hold on
plot(dates(beg:T),etrfit(beg:T,1),'-r')
plot(dates(beg:T),gpf(beg:T,:)*b,'-b');
plot(dates(beg:T),etrfit(beg:T,:),'-r')
title('Eigenvalue model vs CP ad hoc') 
axis([2000 2006 -3 7]); 
legend('Eigenvalue ','CP adhoc')
print -depsc2 cp_fitted_3_eig.eps;



% *******************
% plots for lecture
% *******************

if 0
    % abandoned plots, but here for reference in case some file uses them. 
    % these are full sample. following versions with the same file names use real time
    % itneresting real time fits better pre crisis, because crisis change
    % in coefficients doesn't affect earliser sample! 
    % fitted and actual
    fit_avg = [ones(T,1) 100*yields(:,1) 100*forwards]*betas;
    figure;
    plot(dates(beg:T),fit_avg(beg:T,:),'-r','LineWidth',2);
    hold on
    plot(dates(beg:T-12),100*mean(hprx(beg+12:T,:)')','-b','LineWidth',2);
    plot(dates(beg:T),0*fit_avg(beg:T,:),'-k');
    legend('forecast','actual',3);
    text(2003,-5,'CP Sample ','HorizontalAlignment','right');
    plot([2003 2003]',[-10 12],'-k');
    axis([1990 enddate -6 8]);
    title('Forecast of average (across mat) bond returns')
    print -depsc2 cp_expost_avg.eps;
    
    % fitted and actual with modifications
    
    FT_nolev   = [ones(T-beg-12+1,1)  100*forwards(beg:T-12,:)-100*yields(beg:T-12,1)*ones(1,4)]; % yeilds and forwards
    [betas_nolev,stbetas_nolev,R2_nolev,R2adj_nolev,v_nolev,F__nolev] = olsgmm(100*mean(hprx(beg+12:T,:)')',FT_nolev,12,0);  % std errors: HH with 12 lags
    fit_avg_nolev = [ones(T,1)  100*forwards-100*yields(:,1)*ones(1,4)]*betas_nolev;
    figure;
    plot(dates(beg:T),fit_avg_nolev(beg:T,:),'-r','LineWidth',2);
    hold on
    plot(dates(beg:T-12),100*mean(hprx(beg+12:T,:)')','-b','LineWidth',2);
    plot(dates(beg:T),0*fit_avg(beg:T,:),'-k');
    legend('forecast','actual',3);
    text(2003,-5,'CP Sample ','HorizontalAlignment','right');
    plot([2003 2003]',[-10 12],'-k');
    axis([1990 enddate -6 8]);
    title('Forecast average (across mat) returns with spreads')
    print -depsc2 cp_expost_avg_spreads.eps;
end

% real time forecasts -- CP, spreads. 
FT_nolev   = [ones(T-beg-12+1,1)  100*forwards(beg:T-12,:)-100*yields(beg:T-12,1)*ones(1,4)]; % yeilds and forwards
lhv = 100*mean(hprx(beg+12:T,:)')';
t = 100; % point to start; 
[betas_nolev,stbetas_nolev,R2_nolev,R2adj_nolev,v_nolev,F__nolev] = olsgmm(lhv(1:t),FT_nolev(1:t,:),12,0);  % std errors: HH with 12 lags
fit_avg_nolev = [ones(beg+t,1)  100*forwards(1:beg+t,:)-100*yields(1:beg+t,1)*ones(1,4)]*betas_nolev;
for t = 101:size(FT_nolev,1); 
        [betas_nolev,stbetas_nolev,R2_nolev,R2adj_nolev,v_nolev,F__nolev] = olsgmm(lhv(1:t),FT_nolev(1:t,:),12,0);  % std errors: HH with 12 lags
        fit_avg_nolev =[fit_avg_nolev; FT_nolev(t,:)*betas_nolev];
end;
% make forecasts past actuals up to present
fit_avg_nolev =[fit_avg_nolev; [ones(12,1) 100*forwards(T-12+1:T,:)-100*yields(T-12+1:T,1)*ones(1,4)]*betas_nolev];

figure; 
plot(dates(beg:T),fit_avg_nolev(beg:T,:),'-r','LineWidth',2); 
hold on
plot(dates(beg:T-12),100*mean(hprx(beg+12:T,:)')','-b','LineWidth',2); 
plot(dates(beg:T),0*fit_avg_nolev(beg:T,:),'-k'); 
%plot(ffdates_raw,100*ffunds,'-k'); 
legend('forecast','actual',3);
text(2003,-5,'CP Sample ','HorizontalAlignment','right'); 
plot([2003 2003]',[-10 12],'-k'); 
axis([1990 enddate -6 8]); 
title('CP Real-time forecast average (across mat) returns with spreads') 
print -depsc2 cp_expost_avg_spreads.eps;

% real time forecasts -- CP, levels. 
FT   = [ones(T-beg-12+1,1)  100*yields(beg:T-12,1) 100*forwards(beg:T-12,:)]; % yeilds and forwards
lhv = 100*mean(hprx(beg+12:T,:)')';
t = 100; % point to start; 
[betas,stbetas,R2,R2adj,v,F] = olsgmm(lhv(1:t),FT(1:t,:),12,0);  % std errors: HH with 12 lags
fit_avg = [ones(beg+t,1)  100*yields(1:beg+t,1) 100*forwards(1:beg+t,:)]*betas;
for t = 101:size(FT,1); 
        [betas,stbetas,R2,R2adj,v,F] = olsgmm(lhv(1:t),FT(1:t,:),12,0);  % std errors: HH with 12 lags
        fit_avg =[fit_avg; FT(t,:)*betas];
end;
% make forecasts past actuals up to present
fit_avg =[fit_avg; [ones(12,1) 100*yields(T-12+1:T,1) 100*forwards(T-12+1:T,:)]*betas];

figure; 
plot(dates(beg:T),fit_avg(beg:T,:),'-r','LineWidth',2); 
hold on
plot(dates(beg:T-12),100*mean(hprx(beg+12:T,:)')','-b','LineWidth',2); 
plot(dates(beg:T),0*fit_avg(beg:T,:),'-k'); 
%plot(ffdates_raw,100*ffunds,'-k'); 
legend('forecast','actual',3);
text(2003,-5,'CP Sample ','HorizontalAlignment','right'); 
plot([2003 2003]',[-10 12],'-k'); 
axis([1990 enddate -6 8]); 
title('CP Real-time forecast average (across mat) returns') 
print -depsc2 cp_expost_avg_levels.eps;


% real time forecasts -- FB, spreads. 
FT_nolev   = [ones(T-beg-12+1,1)  mean(100*forwards(beg:T-12,:)-100*yields(beg:T-12,1)*ones(1,4),2)]; 
%"FB" uses mean forward spread to forecast mean return. All rest of code is
% identical past here. 
lhv = 100*mean(hprx(beg+12:T,:)')';
t = 100; % point to start; 
[betas_nolev,stbetas_nolev,R2_nolev,R2adj_nolev,v_nolev,F__nolev] = olsgmm(lhv(1:t),FT_nolev(1:t,:),12,0);  % std errors: HH with 12 lags
fit_avg_nolev = [ones(beg+t,1)  mean(100*forwards(1:beg+t,:)-100*yields(1:beg+t,1)*ones(1,4),2)]*betas_nolev;
for t = 101:size(FT_nolev,1); 
        [betas_nolev,stbetas_nolev,R2_nolev,R2adj_nolev,v_nolev,F__nolev] = olsgmm(lhv(1:t),FT_nolev(1:t,:),12,0);  % std errors: HH with 12 lags
        fit_avg_nolev =[fit_avg_nolev; FT_nolev(t,:)*betas_nolev];
end;
% make forecasts past actuals up to present
fit_avg_nolev =[fit_avg_nolev; [ones(12,1) mean(100*forwards(T-12+1:T,:)-100*yields(T-12+1:T,1)*ones(1,4),2)]*betas_nolev];

figure; 
plot(dates(beg:T),fit_avg_nolev(beg:T,:),'-r','LineWidth',2); 
hold on
plot(dates(beg:T-12),100*mean(hprx(beg+12:T,:)')','-b','LineWidth',2); 
plot(dates(beg:T),0*fit_avg_nolev(beg:T,:),'-k'); 
%plot(ffdates_raw,100*ffunds,'-k'); 
legend('forecast','actual',3);
text(2003,-5,'CP Sample ','HorizontalAlignment','right'); 
plot([2003 2003]',[-10 12],'-k'); 
axis([1990 enddate -6 8]); 
title('FB Real-time Forecast average (across mat) returns with spreads') 
print -depsc2 fb_expost_avg_spreads.eps;


% ------------------------------------
% make pictures for lecture
% ------------------------------------

% --------------------
% load fed funds data used in graphs

x = load('fedfunds.txt'); 
ffdates = datenum(x(:,1),x(:,2),x(:,3)); % dates in matlab date counter units
ffdates_raw = x(:,1)+x(:,2)/12; % in the same style as used for fama bliss plots
ffunds = x(:,4); 
ffunds = ffunds/36000; % was bank interest, now daily interest
ffunds = 360*log(1+ffunds);
% ------------------------------------

figure; 
subplot(2,1,1); 
plot(date_number,100*yields);
hold on
plot(ffdates,100*ffunds,'-k'); 
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits') 
title('Yields of 1-5 year zeros and fed funds'); 
subplot(2,1,2); 
plot(date_number,100*(yields-yields(:,1)*ones(1,5)));
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits');% there is a way to specify axis in real dates
title('Yield spreads y^{(n)}-y^{(1)}'); 
print -depsc2 yields_over_time.eps; 


figure; 
subplot(2,1,1); 
plot(date_number,100*[yields(:,1) forwards]);
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits') 
title('1-5 year forwards'); 
subplot(2,1,2); 
plot(date_number,100*([yields(:,1) forwards]-yields(:,1)*ones(1,5)));
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits');% there is a way to specify axis in real dates
title('forward spreads f^{(n)}-y^{(1)}'); 
print -depsc2 forwards_over_time.eps; 



figure; 
subplot(2,1,1); 
plot(date_number,100*hpr);
hold on
plot(date_number, 0*date_number, '-k');
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits') 
title('Annual returns on 1-5 year bonds'); 
subplot(2,1,2); 
plot(date_number,100*hprx);
hold on
plot(date_number, 0*date_number, '-k');
axis([datenum(1952,01,01) datenum(enddate,01,01) -inf inf]);
datetick('keeplimits');% there is a way to specify axis in real dates
title('Annual excess returns on 2-5 yerar bonds'); 
print -depsc2 returns_over_time.eps; 


index_05 = find(rawdates == 20051230);
figure; 
plot((1:5)',100*yields([index_05-24 index_05-12 index_05],:)','--v',(1:5)', 100*[yields([index_05-24 index_05-12 index_05],1) forwards([index_05-24 index_05-12 index_05],:)]','-o');
%legend('yields 12/04','forwards 12/04','yields 12/03','forwards 12/03',4); 
title('yields and forwards'); 
axis([0.5 5 0 5]); 
set(gca,'ytick',0:5); 
set(gca,'xtick',1:5); 
text([1 1 1],100*yields([index_05-24 index_05-12 index_05],1),['12/03  '; '12/04  '; '12/05  '],'HorizontalAlignment','Right');  
xlabel('maturity in years'); 
ylabel('percent'); 
print -depsc2 end_yield_curve_05.eps; 


index_08 = find(rawdates == 20081231);
figure; 
plot((1:5)',100*yields([index_08-24 index_08-12 index_08],:)','--v',(1:5)', 100*[yields([index_08-24 index_08-12 index_08],1) forwards([index_08-24 index_08-12 index_08],:)]','-o');
%legend('yields 12/04','forwards 12/04','yields 12/03','forwards 12/03',4); 
title('yields and forwards in the crash'); 
axis([0.5 5 0 5]); 
set(gca,'ytick',0:5); 
set(gca,'xtick',1:5); 
text([1 1 1],100*yields([index_08-24 index_08-12 index_08],1),['12/06  '; '12/07  '; '12/08  '],'HorizontalAlignment','Right');  
xlabel('maturity in years'); 
ylabel('percent'); 
print -depsc2 end_yield_curve.eps; 

figure; 
plot((1:5)',100*[yields([index_08-6:index_08],1) forwards([index_08-6:index_08],:)]','-o');
%legend('yields 12/04','forwards 12/04','yields 12/03','forwards 12/03',4); 
title('Forwards in the crash'); 
axis([0.5 5 0 6]); 
set(gca,'ytick',0:5); 
set(gca,'xtick',1:5); 
text([1 1 1 1 1 1 1]',100*yields([index_08-6:index_08],1),['06/08 ';'07/08 ';'08/08 ';'09/08 ';'10/08 '; '11/08 '; '12/08 '],'HorizontalAlignment','Right');  
xlabel('maturity in years'); 
ylabel('percent'); 
print -depsc2 end_yield_curve2.eps; 


figure; 
plot((1:5)',100*yields([end-12 end-6 end],:)','--v',(1:5)', 100*[yields([end-12 end-6 end],1) forwards([end-12 end-6 end],:)]','-o');
title('Recent yields and forwards in the crash'); 
axis([0.5 6 0 5]); 
set(gca,'ytick',0:5); 
set(gca,'xtick',1:5); 
text([5.1 5.1 5.1],100*forwards([end-12 end-6 end],end),num2str(rawdates([end-12 end-6 end]),'%8.0i'),'HorizontalAlignment','Left');  
xlabel('maturity in years'); 
ylabel('percent'); 
print -depsc2 recent_yield_curve.eps; 


% *************************************
% END DATA QUESTIONS
% ********************************

% --------------------------------
% example for forward curve question
% -------------------------------

f0 = [0.01 0.02 0.025 0.03 0.0325]'; %forwards including current yield
n = [ 1 2 3 4 5]';
p0 = -cumsum(f0); 
p1 = p0 + f0(1); % these are now the prices of 0 1 2 3 4 year bonds at time 1. All prices rise at one year rate

f1 = [p1(1:end-1)-p1(2:end)]; 
p2 = p1(2:end) + f1(1);

f2 = [p2(1:end-1)-p2(2:end)]; 
p3 = p2(2:end) + f2(1);

f3 = [p3(1:end-1)-p3(2:end)]; 
p4 = p3(2:end) + f3(1);

f4 = [p4(1:end-1)-p4(2:end)]; 
p5 = p4(2:end) + f4(1);

f5 = [-p5(1)]; 

p = [p0'; p1'; [NaN p2']; [NaN NaN p3'] ; [NaN NaN NaN p4']; [NaN NaN NaN NaN p5']]; 
f = [f0'; [f1' NaN]; [f2' NaN NaN ]; [f3' NaN NaN NaN ] ; [ f4' NaN NaN NaN NaN]]; 

if 0
    figure;
    subplot(2,1,1);
    plot((0:5)', p,'-v');
    xlabel('time');
    ylabel('log price');
    title('price of given bonds over time');
    set(gca,'xtick',[0 1 2 3 4 5])
    set(gca,'ytick',[-0.25 -0.20 -0.15 -0.10 -0.05 0]);
    
    subplot(2,1,2);
    plot((0:4)', 100*f,'-v');
    xlabel('time');
    ylabel('log forwards (%)');
    title('forwards of given maturities over time');
    legend('1','2','3','4','5',4);
    set(gca,'xtick',[0 1 2 3 4 5] )
    axis([0 5 -inf inf]);
    %set(gca,'ytick',[-0.25 -0.20 -0.15 -0.10 -0.05 0]);
    print -depsc2 ps9a.eps;
end


figure
plot((1:5)',100*f(1,:)','-v','linewidth',2)
xlabel('maturity')
ylabel('log forwards (%)')
title('Forward curves over time under EH')
axis([ 0 6 0 3.5]);
legend('time 0','time 1','time 2','time 3', 'time 4', 4)
print -depsc2 ps9_forward_example.eps
print -dpng ps9_forward_example_question.png

figure
plot((1:5)',100*f','-v','linewidth',2)
xlabel('maturity')
ylabel('log forwards (%)')
title('Forward curves over time under EH')
axis([ 0 6 0 3.5]);
legend('time 0','time 1','time 2','time 3', 'time 4', 4)
print -depsc2 ps9_forward_example.eps
print -dpng ps9_forward_example.png

stop


%***************
% Unused questions
% **************

% ------------------
% AR(1) term structure model
% -------------------

% --------
% find delta, rho. Time series regression
% this is a regression at annual horizon. I ran overlapping monthly
% observations. You can also run it once per year. 
% --------
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(yields(beg+12:T,1),[ones(T-beg-12+1,1) yields(beg:T-12,1)],12,0);
disp(' ')
disp('AR(1) term structure model') 
disp(' delta and rho estimates'); 
fprintf('           %10.9s %10.9s \n','delta','rho')
fprintf(' TS        %10.2f %10.2f \n',100*bv(1)/(1-bv(2)),bv(2))
delta0 = bv(1)/(1-bv(2));
rho0 = bv(2);
fit_fwd0 = delta0 + (yields(:,1)-delta0)*rho0.^[1 2 3 4];

% minimize mse % 
method = 1; % 1 lsqoln 2 fminsearch
parms0 = [0.06 0.90];
parms = lsqnonlin(@objective,parms0);
delta = parms(1);
rho = parms(2); 
fprintf('  search1   %10.2f %10.2f \n',100*delta,rho)

method = 2; % 1 lsqoln 2 fminsearch
parms0 = [0.06 0.90];
parms = fminsearch(@objective,parms0);
delta = parms(1);
rho = parms(2); 
fprintf('  search2   %10.2f %10.2f \n',100*delta,rho)

% -------------
% now make plots to see if it worked
% --------------

fit_fwd = delta + (yields(:,1)-delta)*rho.^[1 2 3 4];

figure;
plot(dates(beg:T),100*forwards(beg:T,:),'-b');
hold on
plot(dates(beg:T),100*yields(beg:T,1),'-k');
axis([1964 enddate  0 16 ]);
title('Actual forwards')
print -depsc2 ar1_model_1.eps


figure;
plot(dates(beg:T),100*fit_fwd0(beg:T,:),'-b');
hold on
plot(dates(beg:T),100*yields(beg:T,1),'-k');
plot(dates(beg:T),100*delta0*ones(T-beg+1,1),'-k')
axis([1964 enddate  0 16 ]);
title('Using time series estimate')
print -depsc2 ar1_model_2.eps


figure;
plot(dates(beg:T),100*fit_fwd(beg:T,:),'-b');
hold on
plot(dates(beg:T),100*yields(beg:T,1),'-k');
plot(dates(beg:T),100*delta*ones(T-beg+1,1),'-k')
axis([1964 enddate  0 16 ]);
title('Using mse estimate') 
print -depsc2 ar1_model_3.eps




% -------------------------------------
% old example of yield curve rising through time
% -----------------------------------------------
if 0
    
    y0 = [0.01 0.02 0.03 0.04 0.05]';
    n = [ 1 2 3 4 5]';
    p0 = -n.*y0;
    p1 = p0 + y0(1); % these are now the prices of 0 1 2 3 4 year bonds at time 1.
    
    y1 = -p1(2:end)./n(1:end-1);
    p2 = p1(2:end)+y1(1);
    
    y2 = -p2(2:end)./n(1:end-2);
    p3 = p2(2:end)+y2(1);
    
    y3 = -p3(2:end)./n(1:end-3);
    p4 = p3(2:end)+y3(1);
    
    y4 = -p4(2:end)./n(1:end-4);
    p5 = p4(2:end)+y4(1);
    
    p = [p0'; p1'; [NaN p2']; [NaN NaN p3'] ; [NaN NaN NaN p4']; [NaN NaN NaN NaN p5']];
    
    figure;
    subplot(2,1,1);
    plot((0:5)', p,'-v');
    xlabel('time');
    ylabel('log price');
    title('price of given bonds over time');
    set(gca,'xtick',[0 1 2 3 4 5])
    set(gca,'ytick',[-0.25 -0.20 -0.15 -0.10 -0.05 0]);
    
    y = [y0'; [y1' NaN]; [y2' NaN NaN ]; [y3' NaN NaN NaN ] ; [ y4' NaN NaN NaN NaN]];
    subplot(2,1,2);
    plot((0:4)', 100*y,'-v');
    xlabel('time');
    ylabel('log yield (%)');
    title('yields of given maturities over time');
    legend('1','2','3','4','5',4);
    set(gca,'xtick',[0 1 2 3 4 5] )
    axis([0 5 -inf inf]);
    %set(gca,'ytick',[-0.25 -0.20 -0.15 -0.10 -0.05 0]);
    
    print -depsc2 ps9a.eps;
    
end

clear all

% -------------------------------------------------------------------
% exchange rate problem. Not Used.
% load data files. These are created from datastream data by ps9_uip.m
% ------------------------------------------------------------------

x = load('Three_month_rates_international.txt');
yr = x(:,1);
mo = x(:,2);
dy = x(:,3); 
r3 = x(:,4:end); 
lbls = { 'CDN' , 'GER', 'JPN', 'SWZ', 'USA',  'UK'};

x = load('One_month_rates_international.txt');
r1 = x(:,4:end);

x = load('exchange_rates.txt'); 
spot = x(:,4:end);

L = 0; % use 0 for full sample, 12 to stop regressions in T-12
T = size(spot,1);
N = size(spot,2);

% *******************************
% Make pictures to check data, change to right units
% *******************************

figure
plot(datenum(yr,mo,dy),r1); 
text(datenum(yr(end),mo(end),dy(end))*ones(1,N),r1(end,:),lbls)
datetick
title('1 Month Interest Rates'); 

r1 = log(1+r1/100)/12; % convert to logs and monthly rate

figure
plot(datenum(yr,mo,dy),r3); 
text(datenum(yr(end),mo(end),dy(end))*ones(1,N),r3(end,:),lbls)
datetick
title('3 Month Interest Rates');

r3 = log(1+r3/100)/4; % convert to logs and 3 month rate

figure
plot(datenum(yr,mo,dy),r1-r1(:,5)*ones(1,size(r1,2))); 
text(datenum(yr(end),mo(end),dy(end))*ones(1,N),r1(end,:)-r1(end,4),lbls)
datetick
title('1 Month Log Interest Rates Spread'); 

figure
plot(datenum(yr,mo,dy),spot./(ones(T,1)*spot(end,:))); 
text(datenum(yr(1),mo(1),dy(1))*ones(1,N),spot(1,:)./spot(end,:),lbls)
title('exchange rates');
datetick

spot = log(spot);

figure
for i = 1:6, 
    subplot(2,3,i)
    plot(datenum(yr,mo,dy),spot(:,i)-spot(1,i));
    hold on
    plot(datenum(yr,mo,dy),-50*(r3(:,i)-r3(:,5)),'-r');
    legend('exch rate','i spread');
    title(char(lbls{i}))
    axis([datenum(1975,1,1) datenum(enddate,1,1) -inf inf])
    datetick('x','keeplimits')
end;
    

% ******************
% Set up regressions
% *******************

base = 5; % index of base country, US
% regress exchange rate change on matched interest differential. 
N = size(r1,2); 
T = size(r1,1);
b = zeros(2,N);
se = zeros(2,N);

% ****** basic uip regression 

for indx = 1:N;
    if ~(indx == base); 
        lhv = spot(2:T-L,indx)-spot(1:T-1-L,indx);
        rhv = [ones(T-1-L,1) r1(1:T-1-L,indx)-r1(1:T-1-L,base)];
        [b(:,indx),se(:,indx),R2(indx),R2a,v,F] = olsgmm(lhv,rhv,0,0);
    end;
end;
fprintf('forward premium regressions \n');
fprintf('1 month s_(t+1)-s_t on rf_t - r_d t\n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n',' ','a','b','t','R2')
for indx = 1:N; 
    fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n',lbls{indx},b(:,indx),b(2,indx)/se(2,indx),R2(indx)); 
end;

% pooled 

lhv = spot(2:T-L,1)-spot(1:T-1-L,1);
rhv = [r1(1:T-1-L,1)-r1(1:T-1-L,base)];
for indx = 2:N;
    if ~(indx == base); 
        lhv = [lhv; spot(2:T-L,indx)-spot(1:T-1-L,indx)];
        rhv = [rhv;  r1(1:T-1-L,indx)-r1(1:T-1-L,base)];
    end;
end;

[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','Pooled ',0,bpool,NaN,R2pool);        

rhv2 = [ones((N-1)*(T-L-1),1) rhv];
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','Pooled',bpool,NaN,R2pool);        

rhv2 = [kron(eye(N-1),ones(T-L-1,1)) rhv];
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','countryD',NaN,bpool(end),NaN,R2pool);        

rhv2 = [kron(ones(N-1,1),eye(T-L-1)) rhv]; %time dummy 
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','timeD',NaN,bpool(end),NaN,R2pool);        


for indx = 1:N;
    if ~(indx == base); 
        lhv = spot(4:T-L,indx)-spot(1:T-3-L,indx);
        rhv = [ones(T-3-L,1) r3(1:T-3-L,indx)-r3(1:T-3-L,base)];
        [b(:,indx),se(:,indx),R2(indx),R2a,v,F] = olsgmm(lhv,rhv,3,0);
    end;
end;
fprintf('3 month s_(t+1)-s_t on rf_t - r_d t\n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n',' ','a','b','t','R2')
for indx = 1:N; 
    fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n',lbls{indx},b(:,indx),b(2,indx)/se(2,indx),R2(indx)); 
end;


% pooled 

lhv = spot(4:T-L,1)-spot(1:T-3-L,1);
rhv = [r3(1:T-3-L,1)-r3(1:T-3-L,base)];
for indx = 2:N;
    if ~(indx == base); 
       lhv = [lhv; spot(4:T-L,indx)-spot(1:T-3-L,indx)];
       rhv = [rhv;  r3(1:T-3-L,indx)-r3(1:T-3-L,base)];
    end;
end;

[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','Pooled',0,bpool,NaN,R2pool);        


rhv2 = [ones((N-1)*(T-L-3),1) rhv];
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','Pooled',bpool,NaN,R2pool);        

rhv2 = [kron(eye(N-1),ones(T-L-3,1)) rhv];
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','countryD',NaN,bpool(end),NaN,R2pool);        

rhv2 = [kron(ones(N-1,1),eye(T-L-3)) rhv]; %time dummy 
[bpool,sepool,R2pool,R2a,v,F] = olsgmm(lhv,rhv2,0,0);
fprintf(' %8.8s %8.3f %8.3f %8.3f %8.3f\n','timeD',NaN,bpool(end),NaN,R2pool);        

% ----------------------
% fama bliss vs. expectations term structure models
% This is an intereting failure. Follow up some day
% ----------------------

if 0

 y1 = [    6 6 6 6     4   2  1   1   1   2   4   6   6 6 ]' - 6;

 T=size(y1,1);
 tim = (1:T)';
 
rho = 0.95; 
feh2 = y1 + (rho-1)*y1; 
feh3 = y1 + (rho^2-1)*y1; 
feh4 = y1 + (rho^3-1)*y1; 
feh5 = y1 + (rho^4-1)*y1; 
feh = [feh2 feh3 feh4 feh5];

bfb = [ 0.17 0.53 0.84 0.92 ]';  % fama bliss yield forecast coefficients
ffb2 = y1 + (rho-1)*y1/bfb(1); 
ffb3 = y1 + (rho^2-1)*y1/bfb(2); 
ffb4 = y1 + (rho^3-1)*y1/bfb(3); 
ffb5 = y1 + (rho^4-1)*y1/bfb(4); 
ffb = [ffb2 ffb3 ffb4 ffb5];

disp([[ 1-rho 1-rho^2 1-rho^3 1-rho^4];bfb']);
 
 figure;
 plot(tim,y1,'-vk','linewidth',2); 
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T -inf inf]);
 grid on
 print -depsc2 final_expectations_u.eps
 
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 plot(tim,feh,'--o','linewidth',2); 
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T -inf inf]);
 title('Expectations hypothesis')
 grid on
 legend('y1','f2','f3','f4','f5');
 print -depsc2 final_expectations_u_EH_answer.eps
 
 
 figure
 plot(tim,y1,'-vk','linewidth',2); 
 hold on
 plot(tim,ffb,'--o','linewidth',2); 
 xlabel('Time, years')
 ylabel('log yield y_t^{(1)} and fwd. rates');
 axis([0 T -inf inf]);
 title('Fama-Bliss')
 grid on
 legend('y1','f2','f3','f4','f5');
  print -depsc2 final_expectations_u_FB_answer.eps
 
end
