% Conventional and bias-corrected estimation of macro-finance dynamic term structure models for ten countries
% Bauer, Rudebusch, Wu (2014, AER) Comment on ``Term Premia and Inflation Uncertainty: Empirical Evidence from an International Panel Dataset''
%
% This code is based on MATLAB estimation code provided to us by Jonathan
% Wright (tpestaffine2.m).
%
% The model specification is Wright's baseline model, a macro-finance DTSM
% with unspanned macro risks: interest rates and macro variables under
% P-measure, only interest rates under Q-measure 
%
% Two specifications: 
%   OLS - replication of Wright's original estimates
%   BC - bias-corrected estimates

clc;
clear;
addpath(genpath('jsz_library'));
addpath('BC'); % for bias-corrected estimation
addpath('utilities');

% set seed for random stream -- make BC results reproducible
stream = RandStream('mt19937ar','Seed',0);
RandStream.setGlobalStream(stream);

load 'data/zcbis';
load 'data/macrovars';
cvec={'Canada','CH','Germany','Norway','Sweden','UK','US','Australia','NZ','Japan'};

specId = 2;

% filename for TP estimates and path for parameter estimates
if specId == 1  % OLS
    tp_filename = 'tp2_ols.xls';
    parampath = 'param_ols/param';
elseif specId == 2 % BC2 -- indirect inference bias correction, eigenvalue^P < 1
    tp_filename = 'tp2_bc2.xls';
    parampath = 'param_bc2/param';
elseif specId == 3 % BC3 -- bootstrap bias correction, eigenvalue^P < 1
    tp_filename = 'tp2_bc3.xls';
    parampath = 'param_bc3/param';
end

for k=1:10;
    fprintf('*********************\n');
    fprintf('** %s \n', cvec{k});
    fprintf('*********************\n');
    % create data matrices: m - yield data, m_m - macro data
    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; % ensure same starting date
    j=size(m,1); m=m(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]; % mats in quarters
    dt=m(:,1:2); % dates
    m=0.01*m(:,ns+2); % select maturities
    
    % JSZ normalization -- risk factors are first three PCs of yields
    [coeff, score, latent, tsquare] = princomp(m);
    score=score(:,1:3); % PCs
    bigt=size(score,1); % sample size
    W=coeff(:,1:3)';    % loadings
    
    % estimate short rate equation by OLS
    x=[ones(bigt,1) score];
    [delta,vcov3]=ols(x,m(:,1)/4); % regress short yield on PCs
    delta0=delta(1); delta1=delta(2:4);
    mats=ns/4; % mats in years
    
    switch specId
        case 1 % OLS
            [bhat,sigma,vcov1,vcov2,vcov1full]=varest([score m_m],1,1);
            phi=bhat(:,2:end);
            mu=[mean(score) mean(m_m)]'-(phi*[mean(score) mean(m_m)]');
            disp(['largest P-eigenvalue: ', num2str(max(abs(eig(phi))))]);
        case 2 % BC2 -- indirect inference bias correction
            [phi, mu, sigma] = est_bc_var([score m_m], 1, 1, 1, 5000, 1000, 5, 0);
        case 3 % BC3 -- bootstrap bias correction
            bhat=varest([score m_m],1,1);
            phi_ols=bhat(:,2:end);
            disp(['largest eigenvalue OLS: ', num2str(max(abs(eig(phi_ols))))]);
            % bootstrap bias correction
            phi = 2*phi_ols - m_var(phi_ols(:), 5000, 1, [score m_m]);
            disp(['largest eigenvalue BC: ', num2str(max(abs(eig(phi))))]);
            % impose eigenvalue restriction
            phi = shrink_Phi(phi, phi_ols, 1, 1);
            mu = (eye(5) - phi) * mean([score m_m])';
            % residual covariance matrix
            Xdem = [score m_m] - ones(bigt,1)*mean([score m_m]);
            resids = Xdem(2:bigt,:)' - phi * Xdem(1:(bigt-1),:)';
            sigma = resids * resids'/(bigt-1);
    end
    sigma=chol(sigma)';  % Cholesky decomposition of cov matrix (lower triangular)
    
    % numerical optimization for cross-sectional parameters
    if specId == 1
        % OLS: starting values as Wright's code
        paramstart=[-0.01;-0.04;-0.08;0.08;sigma(1:5,1);sigma(2:5,2);sigma(3:5,3);sigma(4:5,4);sigma(5,5)];
    else
        % BC: use OLS estimates as starting values
        eval(['load param_OLS/param' num2str(k) ' param'])
        paramstart=param;
    end
    options=optimset('MaxFunEvals',500000,'MaxIter',500000,'Display','Off','Algorithm','active-set');
    
    % parameter restrictions imposed within jpslikel
    % Wright's restrictions:
    % 1. lamQ must have difference of at least .01 and sorted descending
    % 2. lamQ must be between -.002 and -.6
    % 3. rinfQ must be between 0 and .4
    % here slightly weaker to get better fit:
    % - lamQ's arbitrarily close to each other
    % - lamQ's maximum -.0001
    [param,fval,exitflag,output]=fminsearch(@(param) jpslikel(param,m,W,mats,sigma*sigma'),paramstart, options);
    disp(output.message);
    if (exitflag~=1); error('numerical optimization algorithm not converged'); end;
    
    ch=[param(5) 0 0; param(6) param(10) 0; param(7) param(11) param(14)];
    [BcP, AcP, AX, BX] = jszLoadings(W,diag(param(1:3)),param(4),ch*ch',mats,1/4);
    [K0Q_cP,K1Q_cP,rho0_cP,rho1_cP] = jszRotation(W,diag(param(1:3)),param(4),1/4,[],[],BX,AX);
    cP=m*W'; % risk factors -- only yields
    mfit=(ones(bigt,1)*AcP)+(cP*BcP);
    % fitting errors
    u=mfit-m;
    disp('Root Mean square fitting error (percentage points)');
    disp(sqrt(mean(mean((100*u).^2))));
    
    % save parameters
    K1Q_X=diag(param(1:3)); rinfQ=param(4);
    eval(['save ', parampath, num2str(k) ' param mu phi sigma K1Q_X rinfQ ch'])
    
    % average expected short rate
    avexp = get_exp([score m_m], delta0, delta1, mu, phi);
    % term premium
    tp=[dt (2*mfit(:,end))-mfit(:,end-5)-(4*avexp)]; tp(:,3)=100*tp(:,3);
    
    if k==1; tpca=tp; end; if k==2; tpch=tp; end; if k==3; tpde=tp; end; if k==4; tpno=tp; end;
    if k==5; tpse=tp; end; if k==6; tpuk=tp; end; if k==7; tpus=tp; end; if k==8; tpaus=tp; end;
    if k==9; tpnz=tp; end; if k==10; tpjp=tp; end;
    
    clear g avexp tp dv;
end;

xlswrite(tp_filename,{'Canada'},'Baseline','C1'); xlswrite(tp_filename,tpca(:,[1 2 3]),'Baseline','A2');
xlswrite(tp_filename,{'Switzerland'},'Baseline','F1'); xlswrite(tp_filename,tpch(:,[1 2 3]),'Baseline','D2');
xlswrite(tp_filename,{'Germany'},'Baseline','I1'); xlswrite(tp_filename,tpde(:,[1 2 3]),'Baseline','G2');
xlswrite(tp_filename,{'Norway'},'Baseline','L1'); xlswrite(tp_filename,tpno(:,[1 2 3]),'Baseline','J34');
xlswrite(tp_filename,{'Sweden'},'Baseline','O1'); xlswrite(tp_filename,tpse(:,[1 2 3]),'Baseline','M13');
xlswrite(tp_filename,{'UK'},'Baseline','R1'); xlswrite(tp_filename,tpuk(:,[1 2 3]),'Baseline','P2');
xlswrite(tp_filename,{'US'},'Baseline','U1'); xlswrite(tp_filename,tpus(:,[1 2 3]),'Baseline','S2');
xlswrite(tp_filename,{'Australia'},'Baseline','X1'); xlswrite(tp_filename,tpaus(:,[1 2 3]),'Baseline','V2');
xlswrite(tp_filename,{'New Zealand'},'Baseline','AA1'); xlswrite(tp_filename,tpnz(:,[1 2 3]),'Baseline','Y2');
xlswrite(tp_filename,{'Japan'},'Baseline','AD1'); xlswrite(tp_filename,tpjp(:,[1 2 3]),'Baseline','AB2');
