function output = sims_ident(results,T,s);
%----------------------------------------------------------------------------
% PURPOSE: 1) performs recursive (orthogonal) identification
%             of a structural VAR
%          2) calculates identified impulse response functions
%             at a horizon of T periods
%          3) calculates the forecast error variance decomposition
%             for the s step ahead forecast 
%----------------------------------------------------------------------------
% USAGE:  output = sims_ident(results,T,s) 
% where:    results = a var output 'structure'
%           T = length of desired impulse response horizon
%           s = forecast error variance horizon
%----------------------------------------------------------------------------
% RETURNS a structure
% output.A0 = LHS matrix in a structural VAR estimates of form:
%             A0y(t)=c + A1y(t-1) + ... + Aqy(t-q) + e(t)
% output.sigma = Cov(e(t))
% output.alpha = RHS matrix of structural VAR estimates 
%                alpha = [c A1 .... Aq]
% output.pIRF = 3-tensor of numerical values of positive impulse 
%               responses where cross-section (i,j) corresponds
%               to the the irf of ith variable to the jth 
%               orthogonalized innovation
% output.nIRF = negitive of output.pIRF
% output.FEVD = 3-tensor of numerical values of forecast error variance
%               decompositions where 'page i' corresponds to the contribution
%               to the MSE from the ith orthogonal innovation 
%----------------------------------------------------------------------------
% written by: Gregory Givens
%             3 March 2003
%----------------------------------------------------------------------------


if ~isstruct(results); 
   error('sims identification requires a var output structure');
end;

output.meth = 'recursive (triangular) identification';
S = results.omega; [A,D] = triang_fact(S);
output.A0 = inv(A); output.sigma = D;
output.alpha = output.A0*results.beta;

% calculate impulse response functions
%--------------------------------------
% if VAR contains a constant, interpret 
% the impulse responses as deviations 
% from a constant mean
%--------------------------------------

    % determine if results contains a constant%
    %-----------------------------------------------
q = results.nlags; n = results.neqs;
p = q*n; r = size(results.beta,2);
if p-r == 0;
    gam = results.beta;            % no constant
end;
if p-r == -1;
    gam = results.beta(:,2:end);   % includes constant
end;
    %-----------------------------------------------

Gamma = [gam; eye(n*(q-1)) zeros((n*(q-1)),n)];
Lambda = [A zeros(n,(n*(q-1))); zeros((n*(q-1)),(n*q))];
for j = 0:1:T;
    page = (Gamma^j)*Lambda;
    ppIRF(:,:,j+1) = page(1:n,1:n)*sqrt(D);
end;
output.Hzn = T;
output.pIRF = reshape(ppIRF,n*n,T+1);
output.nIRF = -output.pIRF;
        
% calculate the forecast error variance decompositions
for i = 1:n;
    for j = 1:s;
        junk = Gamma^(j-1); Phij = junk(1:n,1:n);
        fevd(:,:,j) = Phij*D(i,i)*A(:,i)*A(:,i)'*Phij';
    end;
    FEVD(:,:,i) = sum(fevd,3);
end;
MSE1 = sum(FEVD,3); MSE2 = repmat(MSE1,[1 1 n]); 
output.FEVD = FEVD./MSE2; %----------------------------
                          % decompositions expressed as
                          % fraction of total MSE
                          %----------------------------