% Export forward rates and decompositions
% Bauer, Rudebusch, Wu (2014, AER) Comment on ``Term Premia and Inflation Uncertainty: Empirical Evidence from an International Panel Dataset''

clear;
close all;

est_vec = {'OLS' , 'BC'}; 

N = 5;

% analyze decomposition of forward rates into TP
addpath(genpath('jsz_library'));
addpath('BC');

load 'data/zcbis';
load 'data/macrovars';
cvec={'Canada','Switzerland','Germany','Norway','Sweden','U.K.','U.S.','Australia','New Zealand','Japan'};

% want this order: 'US','Japan','Germany','UK', 'Canada', 'Norway','Sweden','Switzerland', 'Australia', 'NZ'
ki = [7, 10, 3, 6, 1, 4, 5, 2, 8, 9];   % i=1 => k=7 => US

year = zcus(219:end,1);
month = zcus(219:end,2);

dtplot=year+month/12;
dtplot=dtplot(3:3:end); year = year(3:3:end); month = month(3:3:end);


for i=1:10
    k = ki(i);
    
    if k==1; m=zcca; m_m=[yca pica]; end;  if k==2; m=zcch; m_m=[ych pich]; end;
    if k==3; m=zcde; m_m=[yde pide]; end;  if k==4; m=zcno; m_m=[yno pino]; end;
    if k==5; m=zcse; m_m=[yse pise]; end;  if k==6; m=zcuk; m_m=[yuk piuk]; end;
    if k==7; m=zcus; m_m=[yus pius]; end;  if k==8; m=zcaus; m_m=[yaus piaus]; end;
    if k==9; m=zcnz; m_m=[ynz pinz]; end;  if k==10; m=zcjp; m_m=[yjp pijp]; end;
    if size(m,1)>233; m=m(end-232:end,:);  end;
    forw = 2*m(:,42)-m(:,22);
    
    j=size(m,1); 
    m=m(round([(3*((j/3)-floor(j/3)))+1:3:j]),:);  %Go to quarterly data
    forw=forw(round([(3*((j/3)-floor(j/3)))+1:3:j]),:);  %Go to quarterly data
    m_m=m_m(end-size(m,1)+1:end,:);
    
    ns=[1 2 3 4 6 8 12 16 20 24 28 32 36 40];
    dt=m(:,1:2); m=0.01*m(:,ns+2);
    
    [coeff, score, latent, tsquare] = princomp(m);
    score=score(:,1:3);
    bigt=size(score,1);
    
    fhat = zeros(bigt, 2); frn = zeros(bigt, 2); ftp = zeros(bigt, 2); 
    
    for iEst=1:2
        switch iEst
            case 1
                eval(['load param_OLS/param', num2str(k)]);
            case 2
                eval(['load param_bc2/param', num2str(k)]);
        end
    
        x=[ones(bigt,1) score];
        [delta,vcov3]=ols(x,m(:,1)/4);
        delta0=delta(1); delta1=delta(2:4);

        W=coeff(:,1:3)'; mats=ns/4;
        [BcP, AcP, AX, BX] = jszLoadings(W,K1Q_X,rinfQ,ch*ch',mats,1/4);
        [K0Q_cP,K1Q_cP,rho0_cP,rho1_cP] = jszRotation(W,K1Q_X,rinfQ,1/4,[],[],BX,AX);
        phiQ = [K1Q_cP + eye(3), zeros(3,2); zeros(2,5)];

        X = [score m_m];

        % risk factors
        cP=m*W';
        % fitted yields
        mfit=(ones(bigt,1)*AcP)+(cP*BcP);
        % fitted forward rates
        fhat(:,iEst) = 100*(2*mfit(:,end)-mfit(:,end-5));
        frn(:,iEst) = 400*get_exp(X, delta0, delta1, mu, phi);
        ftp(:,iEst) = fhat(:,iEst)-frn(:,iEst);
    end
    
    if k==4 % Norway
        fhat=[nan*ones(32,2);fhat];
        frn=[nan*ones(32,2);frn];
        ftp=[nan*ones(32,2);ftp];
        forw = [nan*ones(32,1);forw];
    elseif k==5 % Sweden
        fhat=[nan*ones(11,2);fhat];
        frn=[nan*ones(11,2);frn];
        ftp=[nan*ones(11,2);ftp];
        forw = [nan*ones(11,1);forw];
    end
    
   data = [year, month, forw, fhat, frn, ftp];

   csvwrite(['export/tp_', num2str(i), '.csv'], data);
   
end;

