% variance is in tot_var
% variables in the y=A_var*y(-1)+ghy*epsilon are in the order lgy2_
% these are then reordered (as in original lgy_) when assigned to oo_.var

oo_=evalin('base','oo_');
M_=evalin('base','M_');
dr_=oo_.dr;

if 0
    lgy_=M_.endo_names;
    nstat=M_.nstatic;
    if isfield(M_,'nspred');% older verisions of Dynare (<4.3.3) had this
        npred=M_.nspred;
    else
        npred=M_.npred+M_.nboth;
    end;
    
    
    listofvariables=str2cell(lgy_);
    shocks=M_.exo_names;
    
    % [i,j]=sort(dr_.order_var);
    lgy2_=lgy_(dr_.order_var,:); % variables as appear in the canonical VAR rapresentation
    [a,b,c]=intersect(deblank(listofvariables),str2cell(lgy2_));
    c=c(:);
    listofvariables=a;
    %position in forecast vector and VMA representation
    for jj=1:length(c);
        eval([deblank(listofvariables{jj}),'_pos=',num2str(length(nstat+[1:npred])+jj),';']);
    end
end
%c is the order in the output; So lgy2_(c,:) should return listofvariables

% Rewrite the law of motion as a state space representation of a VAR
% The VAR matrix includes all the states plus list of variables,
%
%
%
% reordering (from Dynare manual)
% % When reordered, the variables are stored in the following order: static variables, purely predetermined
% % variables (variables that appear only at the current and lagged periods in the model), variables that are
% % both predetermined and forward-looking (variables that appear at the current, future and lagged
% % periods in the model), purely forward-looking variables (variables that appear only at the current and
% % future periods in the model). In each category, the variables are arranged alphabetically.

% The dynamics is fully described by the predetermined variables only.
% so
% % A_var = [dr_.ghx([nstat+[1:npred],c'],:),zeros(size(dr_.ghx([nstat+[1:npred],c'],:),1),length(c))];%[[dr_.ghx(16:37,:);dr_.ghx(4,:);dr_.ghx(11,:);dr_.ghx(38,:)],zeros(size(dr_.ghx,2)+3,3)];
% A_var2 =[zeros(size(lgy_,1),nstat),...
%     dr_.ghx,zeros(size(lgy_,1),...
%     size(lgy_,1)-npred-nstat)];
% % Sigma matrix for VAR representation.
% % Sig_var = dr_.ghu([nstat+[1:npred],c'],:);%[dr_.ghu(16:37,:);dr_.ghu(4,:);dr_.ghu(11,:);dr_.ghu(38,:)] ;
% Sig_var = dr_.ghu;%[dr_.ghu(16:37,:);dr_.ghu(4,:);dr_.ghu(11,:);dr_.ghu(38,:)] ;

[A_var_0,Sig_var]=dynare_VAR_representation;



%% Take deviations from long-run equilibrium (Beveridge-Nelson)
if ~exist('hor','var')
    horizons=1000;
else
    horizons=hor;
end
% Limit distribution (presupposes that max(eig(A_var))<=1;
AT=A_var_0^horizons;

A_var=A_var_0-AT;

% State stpace form
% y=A_var*y(-1)+Sig_var*epsilon;
Sigma_e=M_.Sigma_e;


z=1;
tot_var=zeros(size(A_var));
for uu=1:horizons
    tot_var=A_var*tot_var*A_var'+Sig_var*Sigma_e*Sig_var';
end
% tot_var=tot_var(oo_.dr.inv_order_var,oo_.dr.inv_order_var);
oo_.var=tot_var;
assignin('base','oo_',oo_);
global verbose
if verbose
    mytable(M_.endo_names,{'variance'},[(diag(oo_.var))]);
end