% Analyze and visualize behavior of risk-neutral rates and term premia across countries
% Bauer, Rudebusch, Wu (2014, AER) Comment on ``Term Premia and Inflation Uncertainty: Empirical Evidence from an International Panel Dataset''
% Table 1, table 2, figure 1, figure 2

clear;
close all;

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

M = 2000;  % bootstrap sample
N = 5;

% analyze decomposition of forward rates into TP
addpath(genpath('jsz_library'));
addpath('export_fig');
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

dtplot=zcus(219:end,1)+(zcus(219:end,2)/12);
dtplot=dtplot(3:3:end);

%Read in ECRI dates
ecri=xlsread('data/Ecridates.xls','Sheet1','B362:J594');
ecri=ecri(:,[9 8 3 4 2 5 7 1 6]);  %order the right way
ecri=ecri([3:3:233]',:); %And make quarterly
T = length(ecri(:,1));
ecri = [ecri(:, 1:5), zeros(T,1), ecri(:, 6:9)];

% difference between early and late part of the sample
d_forw = zeros(10,1); d_fhat = zeros(10,2); d_exp = zeros(10,2); d_ftp = zeros(10,2);
%ind_early = 1:8; % first two years
ind_early = 3:7; % first two years
ind_late = 77-4:77; % last five quarters

% Figure 1
left = [.08, .56];
width = .40;
height = .14;
bot = [.03, .23, .43, .63, .83]; bot = bot(5:-1:1);

figure(1);
set(gcf, 'Units', 'inches')
set(gcf, 'Position', [1, 1, 6.5, 8.5])
set(gcf, 'Color', 'w');

figure(2);
set(gcf, 'Units', 'inches')
set(gcf, 'Position', [1.2, 1.2, 6.5, 8.5])
set(gcf, 'Color', 'w');

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);
        
        % confidence intervals for risk-neutral rates
        if (iEst == 2)
            phi_ols = est_var(X, 1, true, false);
            frn_smpl = zeros(bigt, M);
            ftp_smpl = zeros(bigt, M);
            X_sim = gen_var(phi, M, X, 1); % bootstrap sample
            for b = 1:M
                % estimate VAR on each one
                phi_ols_i = est_var(X_sim(:,:,b)', 1, true, false);
                % bias correction
                phi_i = phi_ols_i - (phi_ols - phi);
                % stationarity adjustment
                phi_i = shrink_phi(phi_i, phi_ols_i, 1, 0);
                mu_i = (eye(N) - phi_i) * mean(X)';
                frn_smpl(:,b) = 400*get_exp(X, delta0, delta1, mu_i, phi_i);
                ftp_smpl(:,b) = fhat(:,iEst) - frn_smpl(:,b);
            end
            frn_lb = zeros(bigt,1); frn_ub = zeros(bigt,1);
            ftp_lb = zeros(bigt,1); ftp_ub = zeros(bigt,1);
            for t=1:bigt;
                frn_lb(t) = quantile(frn_smpl(t,:), .05);
                frn_ub(t) = quantile(frn_smpl(t,:), 1-.05);
                ftp_lb(t) = quantile(ftp_smpl(t,:), .05);
                ftp_ub(t) = quantile(ftp_smpl(t,:), 1-.05);
            end
        end
        
        rmse(iEst) = sqrt(mean(mean((100*mfit-100*m).^2))); %#ok<SAGROW>
        maxeig(iEst) = max(abs(eig(phi))); %#ok<SAGROW>
        maxeigQ(iEst) = max(abs(eig(phiQ(1:3,1:3)))); %#ok<SAGROW>

        hmax = 200; h = 20;
        irf = irf_var1(phi, hmax);
        irf_5y(iEst) = irf(h); %#ok<SAGROW>
        try
            halflife(iEst) = find(irf<.5,1,'first'); %#ok<SAGROW>
        catch ME
            halflife(iEst) = NaN; 
        end

        irfQ = irf_var1(phiQ, hmax);
        irfQ_5y(iEst) = irfQ(h); %#ok<SAGROW>
        try
            halflifeQ(iEst) = find(irfQ<.5,1,'first'); %#ok<SAGROW>
        catch ME
            halflifeQ(iEst) = NaN; 
        end            
        
        vol_f(iEst) = std(fhat(:,iEst)); %#ok<SAGROW>
       % vol_f(iEst) = std(forw); %#ok<SAGROW>
        vol_frn(iEst) = std(frn(:,iEst)); %#ok<SAGROW>
        vol_ftp(iEst) = std(ftp(:,iEst)); %#ok<SAGROW>
   
    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];
        frn_lb = [nan*ones(32,1); frn_lb];
        frn_ub = [nan*ones(32,1); frn_ub];
        ftp_lb = [nan*ones(32,1); ftp_lb];
        ftp_ub = [nan*ones(32,1); ftp_ub];
    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];
        frn_lb = [nan*ones(11,1); frn_lb];
        frn_ub = [nan*ones(11,1); frn_ub];
        ftp_lb = [nan*ones(11,1); ftp_lb];
        ftp_ub = [nan*ones(11,1); ftp_ub];
    end
    tp=fhat-frn;

    % averages over early and late part of the sample
    d_forw(i) = 100*(mean(forw(ind_late)) - mean(forw(ind_early)));
    d_fhat(i,:) = 100*(mean(fhat(ind_late,:)) - mean(fhat(ind_early,:)));
    d_exp(i,:) = 100*(mean(frn(ind_late,:)) - mean(frn(ind_early,:)));
    d_ftp(i,:) = 100*(mean(tp(ind_late,:)) - mean(tp(ind_early,:)));
    
    % obtain start and end dates for recessions
    j = 0; % recession
    inRec = false; % currently in recession
    for t = 1:T
      if (ecri(t,i)==1 && ~inRec)
          % new recession
          j = j+1; inRec = true;
          startRec(j) = dtplot(t); %#ok<*SAGROW>
      elseif (ecri(t,i)==0 && inRec)
          % end of recession
          inRec = false;
          endRec(j) = dtplot(t-1);
      end
    end
    if inRec; endRec(j) = max(dtplot); end
    nRec = j;  % number of recessions
    
    % plot Figure 1 -- risk-neutral rates
    figure(1);
    figcol = mod(i-1,2)+1;
    figrow = ceil(i/2);
    subplot('Position', [left(figcol), bot(figrow), width, height]);
    plot(dtplot, fhat(:,1), 'k-', 'Linewidth',1);
    hold on;
    % shaded confidence intervals
    y = frn(:,2); 
    x = dtplot; 
    x = x(~isnan(y)); y = y(~isnan(y));
    p = [frn(:,2)-frn_lb, frn_ub-frn(:,2)]; p = p(~isnan(p(:,1)),:);
    boundedline(x, y, p, 'b-');
    % fitted
    h1 = plot(dtplot, fhat(:,1), 'k-', 'Linewidth',1);
    % OLS
    h2 = plot(dtplot, frn(:,1), 'r:', 'Linewidth',2);
    % BC
    h3 = plot(dtplot, frn(:,2), 'b-', 'Linewidth',2);
    ylabel('Percent');
    axis([1990 2010 -2 12]);
    line([1990 2010],[0 0],'Color','k');
    title(cvec{k});
    if i==2; legend([h1, h2, h3], [{'Fitted'}, est_vec]); end;
%     if (i==10); 
%         export_fig '../../- Wright comment/figures/frn.eps' -eps
%         export_fig 'Figure_1_frn.pdf' -pdf
%     end

    % plot Figure 2 -- term premia
    figure(2);
    subplot('Position', [left(figcol), bot(figrow), width, height]);
    plot(dtplot, ftp(:,1), 'r:', 'Linewidth',2);
    hold on;
    % shade recessions
    curax=axis;
    min_yrange=-2; max_yrange=6;
    y=[min_yrange max_yrange max_yrange min_yrange];
    colorRec = [.7,.7,.7];
    for j=1:nRec
        x=[startRec(j) startRec(j) endRec(j) endRec(j)];
        h_fill = fill(x,y,colorRec);
        set(h_fill,'EdgeColor',colorRec);
    end
    set(gca, 'Layer', 'top')
    % shaded confidence intervals
    % OLS
    h1 = plot(dtplot, ftp(:,1), 'r:', 'Linewidth',2);
    % BC
    h2 = plot(dtplot, ftp(:,2), 'b-', 'Linewidth',2);
    ylabel('Percent');
    axis([1990 2010 -2 6]);
    line([1990 2010],[0 0],'Color','k');
    title(cvec{k});
    if i==2;  h_legend = legend([h1, h2], {'OLS', 'BC'}); end;
%     if (i==10); 
%         export_fig '../../- Wright comment/figures/ftp.eps' -eps
%         export_fig 'Figure_2_ftp.pdf' -pdf
%     end
    
    % Table 1 -- summary statistics
    fprintf('%s \n', cvec{k});
    for iEst = 1:2
        fprintf('& %4s ', est_vec{iEst});
        fprintf('& %5.3f ', rmse(iEst)); 
        fprintf('& %5.4f ', maxeig(iEst)); 
        fprintf('& %5.4f ', maxeigQ(iEst)); 
        fprintf('& %5.2f ', irf_5y(iEst)); 
        fprintf('& %5.2f ', irfQ_5y(iEst));  
        fprintf('& %5.2f ', vol_f(iEst)); 
        fprintf('& %5.2f ', vol_frn(iEst)); 
        fprintf('& %5.2f ', vol_ftp(iEst)); 
        fprintf('\\\\ \n');      
    end
    fprintf('\\hline \n');
   
end;

%% Table 2 -- downward trend in long-term rates
% -> inflation expectations vs. risk-neutral rates

% Consensus survey: first obs is October 1990, last obs is April 2009

consensus0=xlsread('data/Longtermsurveydata.xls','consensus','C3:AF40'); 
%US, Japan, Germany, UK, Canada, Norway, Sweden, Switzerland, Australia, New Zealand
pi_exp = consensus0(:,((1:10)-1)*3+2);
g_exp = consensus0(:,((1:10)-1)*3+1);

for i=[1:5,9]
    k = ki(i);
     g_diff = 100*(mean(g_exp(end-2:end,i)) - mean(g_exp(1:3,i)));
     pi_diff = 100*(mean(pi_exp(end-2:end,i)) - mean(pi_exp(1:3,i)));
     if i==9 % for Australia first obs is NA
        g_diff = 100*(mean(g_exp(end-2:end,i)) - mean(g_exp(2:3,i)));
        pi_diff = 100*(mean(pi_exp(end-2:end,i)) - mean(pi_exp(2:3,i)));
     end
    
    fprintf('%s \n', cvec{k});
    fprintf('& %5.0f ', [pi_diff, g_diff, d_forw(i)]);
    for iEst = 1:2
        if (iEst>1); fprintf('&&                      '); end
        fprintf('& %4s & %5.0f & %5.0f & %5.0f \\\\\n', est_vec{iEst}, d_fhat(i,iEst), d_exp(i,iEst), d_ftp(i,iEst));
    end    
    fprintf('\\hline\n');
end




