%uhlig.m
%
%function that extract Uhlig impulse matrices and computes variance
%decompositions and impulse responses to
%
%Andre Kurmann, FRB, last modified December 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [output,v]=uhlig(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
Atilde = chol(vmat)';

%compute Rtildes for each l
Rtilde=zeros(nvars,nvars,KUbar+1);
FEV=zeros(nvars,nvars,KUbar+2);
for l=0:KUbar;
    C_l=M^l;
    Rtilde(:,:,l+1)=C_l(1:nvars,1:nvars)*Atilde;
end

%compute S
Eii=zeros(nvars,nvars); Eii(1,1)=1;     % ** target variable must be ordered first **
S=zeros(nvars,nvars);
S_l=zeros(nvars,nvars,KUbar+1);
eY=zeros(nvars,nvars*nlags);
eY(1:nvars,1:nvars)=eye(nvars);
for l=0:KUbar;
    S=S+(KUbar+1-max(KLbar,l))*Rtilde(:,:,l+1)'*Eii*Rtilde(:,:,l+1);
end

[V,D]=eig(S);               %eigenvectors have norm=1
lambda=[diag(D) seqa(1,1,nvars)];
order=sortrows(lambda,-1);      %sort eigenvalues in descending order
Vord=[];
for i=1:nvars;
    Vord=[Vord V(:,order(i,2))];    %reorder eigenvector according to ordered
end                                 %eigenvalues            
V=Vord;
                          
%impulse vectors of 2 largest eigenvalues
q1=V(:,1);
alpha=Atilde*q1;   

v = q1'*inv(Atilde)*res';   %1*T series of the fundamental shock 
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(1)=1;       %selection vector for long rate when long rate is first / spread is second
            h1(1)=1; h1(2)=-1;     %selection vector for short rate
        elseif slope == 1;
            hn(2)=1;                %selection vector for long rate when spread is first / long rate is second
            h1(1)=-1; h1(2)=1;      %selection vector for short rate
        end
        hinf = zeros(1,n); hinf(8)=1;   %ATTENTION: ONLY WORKS IF INFLATION IS IN APPROPRIATE POSITIION IN UHLIG-VAR!!
        
        for k=1:nimp;
            Zk1(k,:)=(M^(k-1)*U1)';     
            Zk2(k,:)=(M^(k-1)*U2)';     
            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(1) = 1;   %long bond yield is first element in Uhlig-VAR if slope=0 
        hspread = zeros(1,n); hspread(1) = 1;   %spread is first element in Uhlig-VAR if slope=1
        hffr = zeros(1,n); hffr(nvars) = 1;     %ffr is usually last element in Uhlig-VAR; ATTENTION: ONLY WORKS IF FFR IS IN APPROPRIATE POSITIION!!
        hinf = zeros(1,n); hinf(nvars-1)=1;     %inflation is usually second to  last in Uhlig-VAR; ATTENTION: ONLY WORKS IF INFLATION IS IN APPROPRIATE POSITIION!!
        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);   %nimp x nvars impulse response matrix to 1st shock
    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 first variable in impulse response vector is positive on impact, otherwise revert
    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
    %-----------------------------------------------------------------------
    sigmaks=zeros(nvars,nvars,nimp);
    sigmak=Rtilde(:,:,1)*Rtilde(:,:,1)';
    sigmaks(:,:,1)=sigmak;
    hh1=Rtilde(:,:,1)*q1*(Rtilde(:,:,1)*q1)';
    vardec1(1,:)=(diag(hh1./sigmak))';
    if use_yields == 2;
        if slope == 0;
            hlong = zeros(1,nvars); hlong(1) = 1;   %long bond yield is first element in Uhlig-VAR
            hffr = zeros(1,nvars); hffr(nvars) = 1;     %ffr is usually last element in Uhlig-VAR; ATTENTION: ONLY WORKS IF FFR IS IN APPROPRIATE POSITIION!!
            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(1) = 1;   %long bond yield is first element in Uhlig-VAR 
            hffr = zeros(1,nvars); hffr(nvars) = 1;     %ffr is usually last element in Uhlig-VAR; ATTENTION: ONLY WORKS IF FFR IS IN APPROPRIATE POSITIION!!
            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+Rtilde(:,:,k+1)*Rtilde(:,:,k+1)';
%         sigmaks(:,:,k+1)=sigmak;
        %fraction explained by the two shocks (eqs. 7 & 8)
        hh1=hh1+Rtilde(:,:,k+1)*q1*(Rtilde(:,:,k+1)*q1)';
        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