% this code does the calculations in Nelson, Pinter and Theodoridis,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% written by Gabor Pinter
% Last updated july, 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
global n_count
addpath('functions')

%% Load my Data
ddata = xlsread('dataset_NPT_JAE.xls','data','B3:F166');
%% VAR Estimation
ttrend    = 0;        % (1): time trend; other values no time trend
c         = 1;        % number of deterministic variables
nlags     = 4;        % number of lags of the VAR.
nimp      = 40;       % number of impulses to be computed.
ho        = 4;        % number of horizons over which sign restrictions are imposed.
ndraws    = 500;     % number of posterior draws for Bayesian
maxrot    = 50 ;      % how many random matrices (orthonormal)
LAMBDA    = 1;        % tightness of the prior
nlagplot  = 20;       % number of periods over which we plot the IRFs
n_count   = 1;        % counter of A0 matrices for sign restrictions
xmon      = 3;        % the choleski column for monetary policy shocks

% coloring for figures
colorgraydark = [120/255 120/255 120/255];
colorgray     = [176/255 176/255 176/255];
colorbluedark = [90/255 160/255 255/255];
colorblue     = [153/255 204/255 255/255];
% Bands for IRFs
bandtight     = [0.32 0.68];
bandwide      = [0.16  0.84];

%% start the estimation
    [Traw, n ]     = size(ddata); % number of variables in the VAR
    Y1            = gen_y(ddata,Traw,nlags);
    X1            = gen_x_trend(ddata,Traw,n,nlags,ttrend,c);
    [T, K]         = size(X1);              % dimensions of the endogenous matrix
    [ALPHA, SIGMA] = ols_gabor(X1,Y1,Traw,n,nlags);
    %% -----------------Prior hyperparameters for bvar model
    lambdaP  = LAMBDA;
    tauP     = 10*lambdaP;       % this controls the tightness of the priors on sum of coefficients
    epsilonP = 1/1000;  % this controls tightness of the prior on the constant
    muP      = 1*mean(Y1)';
    %%  Robustness for lags
    [ALPHA_star, Y_star, X_star, SIGMA_star] = bvar_prior_haroon(lambdaP,tauP,epsilonP,muP,Y1,X1,n,nlags,ttrend,c);
    [ALPHA_draws, SIGMA_draws] = bvar_post_haroon(ndraws,n,nlags, X_star, Y_star, ALPHA_star, ttrend, c);
    %% Sign restriction - Robustness for lags
     [histall2, sirall2, histall_trendy, errors_str] = sign_NPT_demandshockINTEREST(ALPHA_draws, SIGMA_draws, n, nlags, ttrend, nimp, maxrot, Y1,X1, ndraws, ho,c,n_count);

histall_trend = histall_trendy;
    histall       = histall2;
    sirall        = sirall2;
    %% Choleski identification - Robustness for lags
    [histall_ch, sirall_ch, histall_trendy_chol, errors_chol] = leverage_choleski(ALPHA_draws, SIGMA_draws, n, nlags, ttrend, nimp,Y1,X1, ndraws, c);
    histall_trend_chol = histall_trendy_chol;
    histall_chol       = histall_ch;
    sirall_chol        = sirall_ch;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Impulse responses
% extract Choleski
sir_chol        = sort(sirall_chol,4);
% extract sign restriction
sir        = sort(sirall,4);
xsize      = size(sir,1);

me11 = zeros(nlagplot,xsize);
up11 = zeros(nlagplot,xsize);
lo11 = zeros(nlagplot,xsize);
up11wide = zeros(nlagplot,xsize);
lo11wide = zeros(nlagplot,xsize);
me11_chol = zeros(nlagplot,xsize);
up11_chol = zeros(nlagplot,xsize);
lo11_chol = zeros(nlagplot,xsize);
up11_cholwide = zeros(nlagplot,xsize);
lo11_cholwide = zeros(nlagplot,xsize);

    lab{1}='Output \newline ';lab{2}='Price Level \newline ';lab{3}='Interest Rate \newline ';lab{4}='Commercial \newline     Bank';
    lab{5}='Shadow \newline   Bank \newline '; lab{6}='Spread';

figure;

for i2=1:xsize
    me11(:,i2) = squeeze(sir(i2,1,1:nlagplot,fix(ndraws*0.50)))/(squeeze(sir(xmon,1,1,fix(ndraws*0.50))));          % Sign restriction -- monetary shock is the 3rd column!
    up11(:,i2) = squeeze(sir(i2,1,1:nlagplot,fix(ndraws*bandtight(2))))/(squeeze(sir(xmon,1,1,fix(ndraws*0.50))));
    lo11(:,i2) = squeeze(sir(i2,1,1:nlagplot,fix(ndraws*bandtight(1))))/(squeeze(sir(xmon,1,1,fix(ndraws*0.50))));
    up11wide(:,i2) = squeeze(sir(i2,1,1:nlagplot,fix(ndraws*bandwide(2))))/(squeeze(sir(xmon,1,1,fix(ndraws*0.50))));
    lo11wide(:,i2) = squeeze(sir(i2,1,1:nlagplot,fix(ndraws*bandwide(1))))/(squeeze(sir(xmon,1,1,fix(ndraws*0.50))));
    
    me11_chol(:,i2) = squeeze(sir_chol(i2,xmon,1:nlagplot,fix(ndraws*0.50)))/(squeeze(sir_chol(xmon,xmon,1,fix(ndraws*0.50)))); % choleski monetary shock is the 3rd column!
    up11_chol(:,i2) = squeeze(sir_chol(i2,xmon,1:nlagplot,fix(ndraws*bandtight(2))))/(squeeze(sir_chol(xmon,xmon,1,fix(ndraws*0.50))));
    lo11_chol(:,i2) = squeeze(sir_chol(i2,xmon,1:nlagplot,fix(ndraws*bandtight(1))))/(squeeze(sir_chol(xmon,xmon,1,fix(ndraws*0.50))));
    up11_cholwide(:,i2) = squeeze(sir_chol(i2,xmon,1:nlagplot,fix(ndraws*bandwide(2))))/(squeeze(sir_chol(xmon,xmon,1,fix(ndraws*0.50))));
    lo11_cholwide(:,i2) = squeeze(sir_chol(i2,xmon,1:nlagplot,fix(ndraws*bandwide(1))))/(squeeze(sir_chol(xmon,xmon,1,fix(ndraws*0.50))));
    
    xx = linspace(1,nlagplot,nlagplot);
    plot_chol = subplot(xsize,2,1+(i2-1)*2);
    a_chol = shadedplot(xx, lo11_chol(:,i2)', up11_chol(:,i2)', colorgraydark, colorgraydark); hold on;
    plot(me11_chol(:,i2),'k-o', 'LineWidth',1), axis tight,
    
    a_chol = shadedplot(xx, lo11_cholwide(:,i2)', lo11_chol(:,i2)', colorgray, colorgray); hold on;
    a_chol = shadedplot(xx, up11_chol(:,i2)', up11_cholwide(:,i2)', colorgray, colorgray); axis tight, hold on;
    
    xlim = get(gca,'xlim');  %Get x range
    plot([xlim(1) xlim(2)],[0 0],'k')
    hold off
    
    if i2 == 1
        title('(a) Cholesky','FontSize',24);
    end
    ylabel(lab{i2},'FontSize',18)
    
    set(gca,'FontSize',20)
    plot_sign = subplot(xsize,2,2+(i2-1)*2);
    ha = shadedplot(xx, lo11(:,i2)', up11(:,i2)', colorbluedark, colorbluedark); axis tight,hold on;
    plot(me11(:,i2),'k-o', 'LineWidth',1), axis tight,
    
    a_chol = shadedplot(xx, lo11wide(:,i2)', lo11(:,i2)', colorblue, colorblue); hold on;
    a_chol = shadedplot(xx, up11(:,i2)', up11wide(:,i2)', colorblue, colorblue); axis tight, hold on;
    
    xlim = get(gca,'xlim');  %Get x range
    plot([xlim(1) xlim(2)],[0 0],'k')
    hold off
    
    if i2 == 1
        title('(b) Sign restrictions','FontSize',24);
    end
    set(gca,'FontSize',20)
end

% save sign restrictions for DSGE matching
% save('dsgematch_feb27','me11','up11wide','lo11wide')

%% Counterfactual historical
% this is the path for interest rate (2003Q3 - 2007Q2) I want to achieve
test2 = [0.7525; 1.4103; 1.8454; 2.5230; 2.8374; 2.9626; 3.0310 ; ...
3.1131; 2.50161; 2.0961; 1.6318; 0.75; 0.2; 0; 0; 0];
test2 = 0.5*test2/100;

ng2 = size(test2,1);

err2      = zeros(ng2+1, ng2);

spath2    = zeros(ng2, ng2);
sumpaths2 = zeros(1,ng2);
ssize2    = zeros(1,ng2);
paths2    = zeros(ng2, ng2);

spath2_sign     = zeros(ng2, ng2);
sumpaths2_sign  = zeros(1,ng2);
ssize2_sign     = zeros(1,ng2);
paths2_sign     = zeros(ng2, ng2);

% calculate the sizes of the series of MP shocks needed to achieve
% the interest rate target
for ii = 1:ng2
    % calculate the MP IRF starting ii to ii+20
    spath2(ii,ii:end)      = squeeze(sir_chol(3,3,1:ng2+1-ii,fix(ndraws*0.50)));
    spath2_sign(ii,ii:end) = squeeze(sir(3,1,1:ng2+1-ii,fix(ndraws*0.50)));

    % calculate ratio of the "gap" that needs filling in with
    % the new MP shock
    ssize2(ii)      = (test2(ii)-sumpaths2(ii))/spath2(ii,ii);
    ssize2_sign(ii) = (test2(ii)-sumpaths2_sign(ii))/spath2_sign(ii,ii);

%   ssize(ii) = err(ii,ii)/spath(ii,ii);
    sumpaths2 = [];
    sumpaths2_sign = [];
    % use this "ssize" multiplier to calculate the MP IRF needed to
    % achieve the interest rate level at time ii
    paths2(ii,ii:end)      = spath2(ii,ii:end)      .* ssize2(ii);
    paths2_sign(ii,ii:end) = spath2_sign(ii,ii:end) .* ssize2_sign(ii);

    sumpaths2      = sum(paths2);
    sumpaths2_sign = sum(paths2_sign);

  end

% use the calculate MP shock sizes to construct output paths
varpath2      = zeros(ng2,n,ng2);
varpath2_sign = zeros(ng2,n,ng2);

for ii = 1:ng2
    mpsize2      = ssize2(ii);    
    mpsize2_sign = ssize2_sign(ii);

    for i2=1:n
    varpath2(ii:end,i2,ii)      =  mpsize2*squeeze(sir_chol(i2,3,1:ng2+1-ii,fix(ndraws*0.50)));          % Sign restriction -- monetary shock is the 3rd column!
    varpath2_sign(ii:end,i2,ii) =  mpsize2_sign*squeeze(sir(i2,1,1:ng2+1-ii,fix(ndraws*0.50)));          % Sign restriction -- monetary shock is the 3rd column!
    end
end
mpvars2 = sum(varpath2,3);
mpvars2_sign = sum(varpath2_sign,3);

% construct counterfactual series
t3      = size(mpvars2_sign,1);
cd_data = [ddata(end-(t3-1)-4:end-(t3-1)-1,:) ; ddata(end-(t3-1):end,:) + mpvars2_sign];
rr_data = ddata(end-(t3-1)-4:end,:);
d4_cd   = cd_data(5:end,:)-cd_data(1:end-4,:);
d4_rr   = rr_data(5:end,:)-rr_data(1:end-4,:);

toexcel = [cd_data rr_data];
beencomCF = 100*(exp(cd_data(end-2,4)) / exp(cd_data(end-15,4)) - 1);
beencomRR = 100*(exp(rr_data(end-2,4)) / exp(rr_data(end-15,4)) - 1);

beenshadCF = 100*(exp(cd_data(end-2,5)) / exp(cd_data(end-15,5)) - 1);
beenshadRR = 100*(exp(rr_data(end-2,5)) / exp(rr_data(end-15,5)) - 1);
%% Text for counterfactual - components of the banking sector
display(['Commercial banks that grew by ' num2str(beencomRR,'%10.2f') '% over 2003Q3-2006Q4,']);
display(['would have grown by ' num2str(beencomCF,'%10.2f') '% over the same period.']);
display(['Equally, shadow banks that grew by ' num2str(beenshadRR,'%10.2f') '% over 2003Q3-2006Q4,']);
display(['would have grown by ' num2str(beenshadCF,'%10.2f') '% over the same period.']);

figure(34);
datet = 2003.5:0.25:2007.25; datet=datet';
subplot(221),
plot(datet,100*d4_cd(:,1),'r--x','MarkerSize',10,'LineWidth',2.2),hold on,
plot(datet,100*d4_rr(:,1),'k-o','MarkerSize',8,'LineWidth',1.8),hold on, axis tight,
title('Four-Quarter Output Growth','FontSize',26),
xlim = get(gca,'xlim');  %Get x range
plot([xlim(1) xlim(2)],[0 0],'k')
hold off
ylabel('%','FontSize',22)
set(gca,'FontSize',18)

subplot(222),
plot(datet,100*cd_data(5:end,3),'r--x','MarkerSize',10,'LineWidth',2.2),hold on,
plot(datet,100*rr_data(5:end,3),'k-o','MarkerSize',8,'LineWidth',1.8),hold on, axis tight,
leg = legend('Counterfactual','Data');
set(leg,'FontSize',24);
title('Interest Rate','FontSize',26),
xlim = get(gca,'xlim');  %Get x range
plot([xlim(1) xlim(2)],[0 0],'k')
hold off
ylabel('%','FontSize',22)
set(gca,'FontSize',18)

subplot(223),
plot(datet,100*d4_cd(:,4),'r--x','MarkerSize',10,'LineWidth',2.2),hold on,
plot(datet,100*d4_rr(:,4),'k-o','MarkerSize',8,'LineWidth',1.8),hold on, axis tight,
title('Four-Quarter Commercial Bank Asset Growth','FontSize',26),
xlim = get(gca,'xlim');  %Get x range
plot([xlim(1) xlim(2)],[0 0],'k')
hold off
ylabel('%','FontSize',22)
set(gca,'FontSize',18)

subplot(224),
plot(datet,100*d4_cd(:,5),'r--x','MarkerSize',10,'LineWidth',2.2),hold on,
plot(datet,100*d4_rr(:,5),'k-o','MarkerSize',8,'LineWidth',1.8),hold on, axis tight,
title('Four-Quarter Shadow Bank Asset Growth','FontSize',26),
xlim = get(gca,'xlim');  %Get x range
plot([xlim(1) xlim(2)],[0 0],'k')
hold off
ylabel('%','FontSize',22)
set(gca,'FontSize',18)

