% current_tfp.m
%
% extracts current TFP shock from a Cholesky decomposition
% i.e. implicit assumption is that only two shocks can affect TFP and that
% news shock has no contemporaneous effect on TFP
%
% Andre Kurmann, Federal Reserve Board; last modified December 2012
%-------------------------------------------------------------------------


function [output,v]=current_tfp(b,res,vmat,nvars,nlags,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 
%------------------------

%unique lower triangular matrix = Cholesky decomp of vmat
Gtilde = chol(vmat)';
gamma = zeros(nvars,1); gamma(1) = 1;  %'selection vector' of same format than in sims extraction
                                       %Note: TFP must be ordered first
alpha = Gtilde*gamma;  %extracts first column of cholesky

% 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;  
       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 on impact, otherwise revert
    %ATTENTION: NEED TO CORRECTLY SPECIFY COLUMN OF RELEVANT VARIABLE !!
    %TFP is target for FEV maximization and TFP is 1st in VAR
    if Zk1(1,1) < 0   
       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
    %-----------------------------------------------------------------------
    B = zeros(nvars,nvars,nimp+1);
    for l=0:nimp;
        C=M^l;
        B(:,:,l+1) = C(1:nvars,1:nvars);
    end
    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;   %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_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
