% tfpnews.m
%
% function that extract TFP news shock as in Barsky and Sims (JME 2011);
% then computes variance decompositions and impulse responses.
%
% Andre Kurmann, FRB; Last modified December 2012
%-----------------------------------------------------------------------

function [output,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields)

%companion matrix of demeaned VAR (=> don't include constant term)
M=zeros(nvars*nlags,nvars*nlags);			
M(1:nvars,:)=b(1:nvars*nlags,:)';
M(nvars+1:nvars*nlags,1:nvars*nlags-nvars)=eye(nvars*nlags-nvars);

% Extract impulse vectors 
%------------------------

%Lower triangular matrix s.t. vcv of fundamental shocks is I matrix = Cholesky decomp of vmat
Gtilde = chol(vmat)';

% compute VMA coeffs Bi
ei = zeros(nvars,1); 
ei(1,1) = 1;      %selection matrix; TFP needs to be ordered first in VAR

B = zeros(nvars,nvars,KUbar+1);
Bi = zeros(1,nvars,KUbar+1);
for l=0:KUbar;
    C=M^l;
    B(:,:,l+1) = C(1:nvars,1:nvars);
    Bi(:,:,l+1) = ei'*C(1:nvars,1:nvars); % Bi are VMA coeffs on how ith variable responds to different shocks l+1 periods back
                                       % Bi(:,:,l+1) is the coef on epsilon(t-l)
end

%compute V
V = zeros(nvars,nvars);
for l=0:KUbar;
    V = V + (KUbar+1-max(KLbar,l))*(Bi(:,:,l+1)*Gtilde)'*(Bi(:,:,l+1)*Gtilde);
end
    
VHat = V(2:nvars,2:nvars);

[eigenVector,eigenValue]=eig(VHat);               %eigenvectors have norm=1
lambda=[diag(eigenValue) seqa(1,1,nvars-1)];
order=sortrows(lambda,-1);     %sort eigenvalues in descending order
Vord=[];
for i=1:nvars-1;
    Vord=[Vord eigenVector(:,order(i,2))];    %reorder eigenvector according to ordered
end                                 %eigenvalues            
                          
%impulse vectors of the largest eigenvalues
gamma = zeros(nvars,1);
gammaHat=Vord(:,1); %VHat(:,1);  % use directely Vord matrix (AK)
gamma = [0;gammaHat];   %first element of gamma = 0 so as to implement assumption that news shock has no contemp effect on TFP (ordered 1st)
alpha=Gtilde*gamma;   %nvars x 1 matrix with alpha in the column 

% Back out time series of structural shock
v = gamma'*inv(Gtilde)*res';    %1xT vector of structural shocks
v = v';

    % Compute impulse responses
    %----------------------------------------------------------------------
    U1=[alpha; zeros(nvars*nlags-nvars,1)];
    n = size(U1,1); 
    nlong = 20;   %ATTENTION: specify number of quarters for long rate: for 5-year = 20

    if use_yields == 0 | use_yields == 3;
        for k=1:nimp;
            Zk1(k,:)=(M^(k-1)*U1)';
        end
        
    elseif use_yields == 1
        hn = zeros(1,n); h1 = zeros(1,n);
        if slope == 0;
            hn(nvars-1)=1;       %selection vector for long rate; long rate is second-to-last of nvars elements if slope=0
            h1(nvars-1)=1; h1(nvars)=-1;     %short rate = long rate - spread; spread is last of nvars elements if slope=0
        elseif slope == 1;
            hn(nvars)=1;       %selection vector for long rate; long rate is last of nvars elements if slope=1
            h1(nvars)=1; h1(nvars-1)=-1;     %short rate = long rate - spread; spread is second to last of nvars elements if slope=1
        end
        hinf = zeros(1,n); hinf(nvars-3)=1;   %inflation; is usually fourth to last in this specification; ATTENTION: ONLY WORKS IF INFLATION IS IN APPROPRIATE POSITIION IN SIMS-VAR!!
    
        for k=1:nimp;
            Zk1(k,:)=(M^(k-1)*U1)';
            LongEH(k,1) = 1/nlong*(h1*inv(eye(n)-M)*(eye(n)-M^nlong)*Zk1(k,:)')'; %long-rate reaction to 1st shock according to Expectations Hypothesis
            TP(k,1) = hn*Zk1(k,:)' - LongEH(k); %term-premium reaction to 1st shock
            SpreadEH(k,1) = LongEH(k,1) - h1*Zk1(k,:)';
        end
        
    elseif use_yields == 2
            hlong = zeros(1,n); hlong(nvars) = 1;   %long bond yield is last element in sims-var if slope=0 
            hspread = zeros(1,n); hspread(nvars) = 1; %long bond yield is last element in sims-var if slope=1
            hffr = zeros(1,n);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!
            hinf = zeros(1,n); hinf(nvars-2)=1;   %selection vector for inflation; ATTENTION: ONLY WORKS IF INFLATION IS IN THIRD-TO-LAST POSITIION IN SIMS-VAR!!
            for k=1:nimp;
                Zk1(k,:)=(M^(k-1)*U1)';
                if slope == 0;
                    Spread(k,1)=(hlong-hffr)*Zk1(k,:)';     %Spread computed from long - ffr
                    Long(k,1)=hlong*Zk1(k,:)';
                elseif slope == 1;
                    Long(k,1) = (hspread+hffr)*Zk1(k,:)';   %Long rate computed from spread + ffr
                end
                LongEH(k,1) = 1/nlong*(hffr*inv(eye(n)-M)*(eye(n)-M^nlong)*Zk1(k,:)')'; %long-rate reaction to 1st shock according to Expectations Hypothesis
                TP(k,1) = Long(k,1) - LongEH(k);    %term-premium reaction to 1st shock
                SpreadEH(k,1) = LongEH(k,1) - hffr*Zk1(k,:)';
            end
            
    end
    
    impulse1(:,:)=Zk1(:,1:nvars);
    if use_yields == 1 | use_yields == 2
       impulse_add1(:,1) = LongEH;  %stock variables computed from VAR expectations in separate matrix
       impulse_add1(:,2) = TP;
       impulse_add1(:,3) = SpreadEH;
       if use_yields == 2
            if slope == 0
                   impulse_add1(:,4) = Spread;
            elseif slope == 1;
                   impulse_add1(:,4) = Long;
            end
       end
    end 
    
    %test whether TFP impulse response is positive, otherwise revert
    %ATTENTION: NEED TO CORRECTLY SPECIFY COLUMN OF RELEVANT VARIABLE !!
    %TFP is target for FEV maximization and thus 1st in VAR
    if Zk1(5,1) < 0   %in order to avert 180 degree rotations, condition on response after 5 quarters to be positive; otherwise mirror;
       alpha = - alpha; 
       v = -v;
       impulse1 = -impulse1;
       if use_yields == 1 | use_yields == 2;
            impulse_add1 = -impulse_add1;
       end
    end


    % Compute fraction of VD explained at different horizons
    %-----------------------------------------------------------------------
    sigmak=B(:,:,1)*vmat*B(:,:,1)';
    hh1=B(:,:,1)*alpha*(B(:,:,1)*alpha)';
    vardec1(1,:)=(diag(hh1./sigmak))';
    if use_yields == 2;
        if slope == 0;
            hlong = zeros(1,nvars); hlong(nvars) = 1;   %long bond yield is last element in Sims-VAR
            hffr = zeros(1,nvars);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!
            sigmak_spread = (hlong-hffr)*sigmak*(hlong-hffr)';    %FEV of spread computed from FEV(long-ffr)
            hh1_spread = (hlong-hffr)*hh1*(hlong-hffr)';
            vardec1_add(1,1) = hh1_spread/sigmak_spread; 
        elseif slope == 1;
            hspread = zeros(1,nvars); hspread(nvars) = 1;   %spread is last element in Sims-VAR 
            hffr = zeros(1,nvars);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!          
            sigmak_long = (hspread+hffr)*sigmak*(hspread+hffr)';    %FEV of long rate computed from FEV(spread+ffr)
            hh1_long = (hspread+hffr)*hh1*(hspread+hffr)';
            vardec1_add(1,1) = hh1_long/sigmak_long;
        end
    end

    for k=1:nimp-1;
        %add square of k-step ahead forecast error to build k-ahead variance-covariance (eq. 6)
        sigmak=sigmak+B(:,:,k+1)*vmat*B(:,:,k+1)';
        %fraction explained by the news shock (eq 8)    % XX change in descripion (AK)
        hh1=hh1+B(:,:,k+1)*alpha*(B(:,:,k+1)*alpha)';
        vardec1(k+1,:)=(diag(hh1./sigmak))';
        if use_yields == 2;
            if slope == 0;
                sigmak_spread = (hlong-hffr)*sigmak*(hlong-hffr)';    %FEV of spread computed from FEV(long-ffr)
                hh1_spread = (hlong-hffr)*hh1*(hlong-hffr)';
                vardec1_add(k+1,1) = hh1_spread/sigmak_spread;    
            elseif slope == 1;
                sigmak_long = (hspread+hffr)*sigmak*(hspread+hffr)';    %FEV of long rate computed from FEV(spread+ffr)
                hh1_long = (hspread+hffr)*hh1*(hspread+hffr)';
                vardec1_add(k+1,1) = hh1_long/sigmak_long;
            end
        end
    end      
    
if use_yields == 0 | use_yields == 3   
    output=[vardec1 impulse1];
elseif use_yields == 1
    output=[vardec1 impulse1 impulse_add1];  
elseif use_yields == 2;
    output=[vardec1 impulse1 impulse_add1 vardec1_add];
end