% Code replicating main results of "The shocks matter: improving our
% estimates of exchange rate pass-through" by Kristin Forbes, Ida Hjortsoe
% and Tsvetelina Nenova, 2018.

% This code estimates a Bayesian VAR model on UK data, saving all results
% in Matlab and Excel files.

clear all
clc
addpath([pwd '/functions/'])


%% Load and set up data

[num,txt,raw] = xlsread('data_JIE.xlsx','data_base');
dates = datenum(txt(3:end,1),'dd/mm/yyyy');
first = find(dates == datenum('31/03/1993','dd/mm/yyyy'));
last = size(num,1)-find(dates == datenum('31/03/2015','dd/mm/yyyy'));

DateNumber = datenum(txt(2+first:end-last,1),'dd/mm/yyyy');
data = num(first:end-last,:);

% %Data transformation (already done in dataset provided)
% Detrended = [3];
% [Trend,HPData] = hpfilter(data(:,3),Inf); % detrend interest rate
% data(:,Detrended) = HPData;
% LogDifferenced = [1,2,4,5,6];
% data(:,LogDifferenced) =
% 100*(log(data(2:end,LogDifferenced))-log(data(1:end-1,LogDifferenced));

p = 2; % VAR lag order

Y = data(p+1:end,:);
[n,k] = size(Y);
X = [ones(n,1), xlags(data,p)];


%% Priors for Bayesian VAR

%compute standard deviation of each series residual via an ols regression
[s,b0_prior] = se_prior(Y,X,k);

%specify uncertainty parameters for Minnesota prior
lamda1 = 0.2; %own lag
lamda2 = 0.5; %other variables' lags
lamda3 = 1; %decay with lag
lamda4 = 10000; %constant term

%specify prior mean of VAR coefficients
B0 = bmean_prior(b0_prior,k,p,1);

%specify prior variance of VAR coefficients
H = bvar_prior(s,lamda1,lamda2,lamda3,lamda4,k,p);

%prior scale matrix for VAR covariance
S = eye(k);
alpha = k+1; %prior degrees of freedom

%starting values for the Gibbs sampling algorithm
Sigma = eye(k);
betaols = vec(inv(X'*X)*(X'*Y));

Reps = 11000;
burn = 10000;


%% Short run, long run and sign restrictions

% f = 2.k x k matrix contain ones and zeros. The zeros occur where the zero
%   restrictions are on the short and long-run impact matrices.  The top k x
%   k block are the short run restrictions, while the bottom k x k
%   restrictions are the long run restrictions.

% sr = k x k matrix of sign restrictions

%  s,   d,  mp,  er,  s*,  d*
f = [1,1,1,1,1,1; %gdp
    1,1,1,1,1,1; %cpi
    1,1,1,1,1,1; %shadow BR
    1,1,1,1,1,1; %ERI(nom)
    1,1,1,1,1,1; %import prices
    0,0,0,0,1,1; %wxpoil
    %=========
    1,0,0,0,1,0; %gdp
    1,1,1,1,1,1; %cpi
    1,1,1,1,1,1; %shadow BR
    1,1,1,1,1,1; %ERI(nom)
    1,1,1,1,1,1;%import prices
    0,0,0,0,1,1];%wxpoil

%  shocks:s,    d,  mp,  er,  s*,  d*
sr = [   +1,   +1,  -1, nan, nan, nan;  %gdp
    -1,   +1,  -1,  -1, nan, nan;  %cpi
    nan,   +1,  +1,  -1, nan, nan;  %interest rate
    nan,   +1,  +1,  +1, nan, nan;  %exchange rate
    nan,  nan, nan, nan, nan, nan;  %import prices
    nan,  nan, nan, nan,  +1,  +1]; %wxpoil

[Q,index,flag] = findQs(k,f);

if flag == 1
    error('Rank condition not satisfied, the model is overidentified');
end


%% Prepare empty matrices for results

shocks = eye(k);
T = 40; % horizon for impulse responses
R = zeros(k,T,Reps-burn,length(shocks)); % Contains impulse resonse functions
cumR = zeros(k,T,Reps-burn,length(shocks)); % Contains cumulative impulse resonse functions
% 1st dimension = variable
% 2nd dimension = time
% 3rd dimension = draw from posterior
% 4th dimension = shock
shock_pos = logical(shocks); % position of shock
var_pos = [1,2,3,4,6,6];       % position of corresponding variable to shock
% eg positive monetary policy shock should result in an increase in interest
% rates, aggregate demand shock should result in a positive increase in gdp
% etc.

A0big = [];
B0big = [];
HD_big = nan(k*(n+T),Reps-burn,k+1);
U_big = [];
Forecasts_big = [];


%% Main Bayesian loop
tic;

counter = 1;
rng('shuffle','twister');
while counter < Reps+1
    if counter < burn+1
        %draw beta & Sigma from posterior
        [beta, Sigma, e] = drawpost(H,Sigma,S,alpha,X,Y,B0,betaols,n,k,p);
        
        counter = counter + 1
        
    else
        % Draw impact matrix candidate (W)
        chck=-1;
        while chck<0  %check for stability
            %draw beta & Sigma from posterior
            [beta, Sigma, e] = drawpost(H,Sigma,S,alpha,X,Y,B0,betaols,n,k,p);
            
            C1 = chol(Sigma,'lower');
            
            %draw candidate impact matrix W
            C = generateDraw(C1,k); % orthonormal matrix Q such that new C=C1*Q
            betanew = reshape(beta,k*p+1,k);
            P = findP(C,betanew,Q,p,k,index);
            W = C*P; % C1*Q*P
            betanewBIG = [ betanew(2:end,:)';  eye(k*(p-1))  zeros(k*(p-1),k) ]; %companion form matrix for IRFs
            ee = max(abs(eig(betanewBIG)));
            if ee<=1
                chck=1;
            end
        end
        
        % Check signs for each column of W (speeds up algorithm)
        A = nan(k,k);
        
        [VV, cumVV] = irfs(1,W,var_pos,shocks,k,p,T,betanewBIG);
        [tmp1_supply, tmp2_supply] = signchk2(1,VV,sr,k);
        
        if cumVV(4,4)~=0 && all(tmp1_supply==0) && all(tmp2_supply==0) %avoid 4th quarter ERI responses equal to exactly 0 (denominator of rescaled IRFs)
            A(:,1) = W(:,1);
            R(:,:,counter-burn,1) = VV(1:k,:);
            cumR(:,:,counter-burn,1) = cumVV(1:k,:);
            for jj = 2:3 % using columns 2-4 gives the same zero restrictions
                for ii = 2:3 % only applying procedure to demand and monetary policy shocks as it does not allow imposing >= restriction on ir
                    
                    Wnew = W;
                    Wnew(:,ii) = W(:,jj);
                    Wnew(:,jj) = W(:,ii);
                    
                    [VV, cumVV] = irfs(ii,Wnew,var_pos,shocks,k,p,T,betanewBIG);
                    [tmp1, tmp2] = signchk2(ii,VV,sr,k);
                    
                    if cumVV(4,4)~=0 && all(tmp1==0) && all(isnan(A(:,ii))) && all(tmp2==0)
                        A(:,ii) = Wnew(:,ii);
                        R(:,:,counter-burn,ii) = VV(1:k,:);
                        cumR(:,:,counter-burn,ii) = cumVV(1:k,:);
                        break %move on to next column of W
                    end
                end
            end
            
            if jj == 3 && all(all(~isnan(A(:,1:3))))
                [VV, cumVV] = irfs(4,W,var_pos,shocks,k,p,T,betanewBIG);
                
                if cumVV(4,4)~=0 && VV(4,1)>0 && VV(3,1)<=0 && VV(2,1)<0 && VV(4,2)>0 && VV(3,2)<=0 && VV(2,2)<0
                    A(:,4) = W(:,4);
                    R(:,:,counter-burn,4) = VV(1:k,:);
                    cumR(:,:,counter-burn,4) = cumVV(1:k,:);
                    if all(all(~isnan(A(:,1:4))))
                        for pp = 5:6
                            [VV, cumVV] = irfs(pp,W,var_pos,shocks,k,p,T,betanewBIG);
                            [tmp1, tmp2] = signchk2(pp,VV,sr,k);
                            
                            if cumVV(4,4)~=0 && all(tmp1==0) && all(tmp2==0)
                                A(:,pp) = W(:,pp);
                                R(:,:,counter-burn,pp) = VV(1:k,:);
                                cumR(:,:,counter-burn,pp) = cumVV(1:k,:);
                            end
                        end
                        if all(all(~isnan(A)))
                            A0big = [A0big, reshape(A,k*k,1)];
                            B0big = [B0big, reshape(betanew,k*(k*p+1),1)];
                            
                            % Save historical decomposition for each conforming draw
                            us = e*inv(A'); %extract structural shocks
                            U_big = [U_big, reshape(us,k*n,1)];
                            
                            for ii = 1:k+1
                                % Empty matrix to save decomposition
                                xmat = hd_xmat(n,k,T,p,X);
                                
                                if ii==1
                                    % Contribution of base and constant
                                    hd = hd_hd(0,n,k,T,p,xmat,betanew,A,us);
                                    HD_big(:,counter-burn,ii) = reshape(hd,k*(n+T),1);
                                else
                                    % Contribution of shocks
                                    hd = hd_hd(ii-1,n,k,T,p,xmat,betanew,A,us);
                                    HD_big(:,counter-burn,ii) = reshape(hd,k*(n+T),1)-HD_big(:,counter-burn,1);
                                end
                            end
                            
                            counter = counter + 1
                        end
                    end
                end
            end
        end
    end
end

toc;

%% FEVD: average across draws
[FEVD,var1,var2,var3,var4,var5,var6] = fevd(R,k,T,Reps,burn);

%% Historical decomposition
U = reshape(mean(U_big,2),n,k);
HD = hd_mean(HD_big,n,k,T,p,DateNumber);

%% IRF graphs
% Scale impulse responses to deliver 1% sterling appreciation in initial
% period
for kk = 1:k
    Rnew(:,:,:,kk) = R(:,:,:,kk)./median(squeeze(cumR(4,4,:,kk)),1);
    cumRnew(:,:,:,kk) = cumR(:,:,:,kk)./median(squeeze(cumR(4,4,:,kk)),1);
end

irfplot(Rnew,cumRnew,T,k,shocks);

%% ERPT ratios with 68% confidence
[ERtoPM, ERtoCPI] = erpt_ratios(cumRnew,50,16,84);

%% Save Matlab file
save('output\BVAR_results.mat');

%% Save Excel sheets
xlswrite('output\FEVD.xlsx', var1, 'var1');
xlswrite('output\FEVD.xlsx', var2, 'var2');
xlswrite('output\FEVD.xlsx', var3, 'var3');
xlswrite('output\FEVD.xlsx', var4, 'var4');
xlswrite('output\FEVD.xlsx', var5, 'var5');
xlswrite('output\FEVD.xlsx', var6, 'var6');

xlswrite('output\HD.xlsx',HD,'HD');
xlswrite('output\HD.xlsx',Y,'Y');
xlswrite('output\HD.xlsx',U,'Shocks');

xlswrite('output\ERPT.xlsx', ERtoPM, 'ERtoPM');
xlswrite('output\ERPT.xlsx', ERtoCPI, 'ERtoCPI');

