%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Matlab code for producing Tables 1-4 in Gurkaynak and Wright   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z=xlsread('feds200628.xls','Sheet1','A2538:K12117');    %Read in data: Aug 1971-Dec 2009
g=find(z(1:end-1,2)~=z(2:end,2));                       %Take end-of-month data
z=z(g,:);

%First Table 1
nvec=[2;4;6;8;10]; mvec=[0.25;0.5;1];
for i=1:length(mvec); for k=1:length(nvec);             %Loop over m and n
    m=mvec(i); n=nvec(k);
    h=round(12*m);                                      %Horizon in months
    bigt=length(z);
    for j=1:bigt-h;
        if h==3; ym=z(j,4); end;                        
        if h==6; ym=z(j,5); end;
        if h>6; ym=svenssonzc(z(j,6:11)',m); end;
        x(j,1)=1;                                   
        x(j,2)=(svenssonzc(z(j,6:11)',n)-ym)*m/(n-m);   %Regressors in equation (13)
        y(j,1)=svenssonzc(z(j+h,6:11)',n-m)-svenssonzc(z(j,6:11)',n); %Dependent variable in equation (13)
    end;
    [bhat,vcov,rsq(i,k)]=newey(y,x,h);
    est_intercept(i,k)=bhat(1);                         %Save the output
    se_intercept(i,k)=sqrt(vcov(1,1));
    tstat_intercept(i,k)=est_intercept(i,k)/se_intercept(i,k);
    est_slope(i,k)=bhat(2);      
    se_slope(i,k)=sqrt(vcov(2,2));
    tstat_slope(i,k)=(est_slope(i,k)-1)/se_slope(i,k);
    clear x y;
end; end;
disp('Table 1: Estimates: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',est_intercept)
disp('Table 1: Standard Errors: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',se_intercept)
disp('Table 1: T-statistics: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',tstat_intercept);
disp('Table 1: Estimates: Slope');
fprintf('%6.2f %6.2f %6.2f \n',est_slope)
disp('Table 1: Standard Errors: Slope');
fprintf('%6.2f %6.2f %6.2f \n',se_slope)
disp('Table 1: T-statistics: Slope');
fprintf('%6.2f %6.2f %6.2f \n',tstat_slope);
disp('Table 1: R-squared');
fprintf('%6.2f %6.2f %6.2f \n',rsq);
clear est_intercept se_intercept tstat_intercept est_slope se_slope tstat_slope rsq

%Next Table 2
for i=1:length(mvec); for k=1:length(nvec);             %Loop over m and n
    m=mvec(i); n=nvec(k);
    capk=round(n/m);
    h=round(12*m);                                      %Horizon in months
    bigt=length(z);
    for j=1:bigt-round(12*n);
        pfspread=0;                                     %Computation of perfect foresight spread
        if h==3; 
            ym=z(j,4); 
            for o=1:capk-1; pfspread=pfspread+((1-(o/capk))*(z(j+(o*h),4)-z(j+(o*h)-h,4))); end;
        end;
        if h==6; 
            ym=z(j,5); 
            for o=1:capk-1; pfspread=pfspread+((1-(o/capk))*(z(j+(o*h),5)-z(j+(o*h)-h,5))); end;
        end;
        if h>6;
            ym=svenssonzc(z(j,6:11)',m); 
            for o=1:capk-1; pfspread=pfspread+((1-(o/capk))*(svenssonzc(z(j+(o*h),6:11)',m)-svenssonzc(z(j+(o*h)-h,6:11)',m))); end;
        end;
        x(j,1)=1;                                       
        x(j,2)=(svenssonzc(z(j,6:11)',n)-ym);           %Regressors in equation (14)
        y(j,1)=pfspread;
        end;
    [bhat,vcov,rsq(i,k)]=newey(y,x,h);
    est_intercept(i,k)=bhat(1);                         %Save the output
    se_intercept(i,k)=sqrt(vcov(1,1));
    tstat_intercept(i,k)=est_intercept(i,k)/se_intercept(i,k);
    est_slope(i,k)=bhat(2);      
    se_slope(i,k)=sqrt(vcov(2,2));
    tstat_slope(i,k)=(est_slope(i,k)-1)/se_slope(i,k);
    clear x y;
end; end;
disp('Table 2: Estimates: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',est_intercept)
disp('Table 2: Standard Errors: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',se_intercept)
disp('Table 2: T-statistics: Intercept');
fprintf('%6.2f %6.2f %6.2f \n',tstat_intercept);
disp('Table 2: Estimates: Slope');
fprintf('%6.2f %6.2f %6.2f \n',est_slope)
disp('Table 2: Standard Errors: Slope');
fprintf('%6.2f %6.2f %6.2f \n',se_slope)
disp('Table 2: T-statistics: Slope');
fprintf('%6.2f %6.2f %6.2f \n',tstat_slope);
disp('Table 2: R-squared');
fprintf('%6.2f %6.2f %6.2f \n',rsq);
clear est_intercept se_intercept tstat_intercept est_slope se_slope tstat_slope rsq

    
%Next Table 3: top panel
for i=1:length(nvec);
    n=nvec(i);                                          %Loop over n
    bigt=length(z);
    for j=1:bigt-12;
        %p_present is the price of an n-year bond today
        %p_future is the price of an n+1-year bond in one year
        p_present=exp(-n*0.01*svenssonzc(z(j,6:11)',n));    
        p_future=exp(-(n-1)*0.01*svenssonzc(z(j+12,6:11)',n-1));
        y1=svenssonzc(z(j,6:11)',1);
        y2=svenssonzc(z(j,6:11)',2);
        y3=svenssonzc(z(j,6:11)',3);
        y4=svenssonzc(z(j,6:11)',4);
        y5=svenssonzc(z(j,6:11)',5);                    %Yields
        y(j,1)=(100*log(p_future/p_present))-y1;        %Excess returns
        x(j,1)=1;                                       %Predictors in (15)
        x(j,2)=y1;
        x(j,3)=(3*y3)-(2*y2);
        x(j,4)=(5*y5)-(4*y4);
        end;
    [bhat,vcov,rsq(i)]=newey(y,x,h);
    fstat(i)=bhat(2:end)'*inv(vcov(2:end,2:end))*bhat(2:end);
    est(:,i)=bhat;                                      %Save the output
    se(:,i)=sqrt(diag(vcov));
    tstat(:,i)=est(:,i)./se(:,i);
    clear x y;
end;
disp('Table 3: Estimates');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',est')
disp('Table 3: Standard Errors');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',se')
disp('Table 3: T-statistics');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',tstat')
disp('Table 3: F-statistics');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',fstat)
disp('Table 3: R squared');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',rsq)
clear est se tstat rsq fstat

%Finally repeat truncating the sample in Dec 2006 for Table 4: bottom panel
z=z(1:end-36,:); 
for i=1:length(nvec);
    n=nvec(i);
    bigt=length(z);
    for j=1:bigt-12;
        p_present=exp(-n*0.01*svenssonzc(z(j,6:11)',n));
        p_future=exp(-(n-1)*0.01*svenssonzc(z(j+12,6:11)',n-1));
        y1=svenssonzc(z(j,6:11)',1);
        y2=svenssonzc(z(j,6:11)',2);
        y3=svenssonzc(z(j,6:11)',3);
        y4=svenssonzc(z(j,6:11)',4);
        y5=svenssonzc(z(j,6:11)',5);
        y(j,1)=(100*log(p_future/p_present))-y1;
        x(j,1)=1;
        x(j,2)=y1;
        x(j,3)=(3*y3)-(2*y2);
        x(j,4)=(5*y5)-(4*y4);
        end;
    [bhat,vcov,rsq(i)]=newey(y,x,h);
    fstat(i)=bhat(2:end)'*inv(vcov(2:end,2:end))*bhat(2:end);
    est(:,i)=bhat;
    se(:,i)=sqrt(diag(vcov));
    tstat(:,i)=est(:,i)./se(:,i);
    clear x y;
end;
disp('Table 4: Estimates');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',est')
disp('Table 4: Standard Errors');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',se')
disp('Table 4: T-statistics');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',tstat')
disp('Table 4: F-statistics');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',fstat)
disp('Table 4: R squared');
fprintf('%6.2f %6.2f %6.2f %6.2f %6.2f \n',rsq)
clear est se tstat rsq fstat
