% variance
% [varY]=variance(PHItl,Pi,Msh,Sigma,nsh);
function [varY,covy,covyHP]=variance(PHI,Pi,Msh,Sigma,hor,asyn,HP);
varY=[];covy=[];covyHP=[];
szY=size(Pi,1);
% PHI=PHItl;
szphi=size(PHI,1);
% keyboard
% szphitl=size(PHItl,1);
% szphi=sqrt(szphitl);
M=Msh;
if nargin<7;
    HP=0;
end
if nargin<6;
    if szphi>20;
        asyn=0;
    else
        asyn=1;
    end
end
if HP==0
    if asyn==0
        if nargin<5
            hor=2000;
        end
        %N=PHI(1:nsh,1:nsh);
        % PHI1=PHI(nsh+1:end,nsh+1:end);
        % PHI2=PHI(nsh+1:end,1:nsh);
        
        % [Msh;
        %     zeros(szphi-nsh,nsh)];
        
        varS=zeros(szphi,szphi);
        shock=M*Sigma*M';
        PHIt=PHI';
        % zz=1;
        for jj=1:hor
            % if jj==1
            
            varS=PHI*varS*PHIt+shock;
            % else
            %      varS=PHI*varS*PHI';
            % end
            
            
            
            
        end
        
        
        covy=Pi*varS*Pi';
    else
        % vars_v=(eye(size(PHI,1)^2)-kron_d(PHI,PHI'))\vec(Msh*Sigma
        varS=(eye(szphi^2)-kron_d(PHI,PHI'))\vec(M*Sigma*M');
        % varY=diag(Pi*reshape(varS,szphi,szphi)*Pi');
        covy=Pi*reshape(varS,szphi,szphi)*Pi';
    end
    
    varY=diag(covy);
    
    % variance by simulation if HP filter
else %HP
    
    periods=500;
    drop=round(.6*periods);
    randn('seed',1984);
    
    covyHP=zeros(size(Pi,1),periods-drop,20);
    
    [nzsig,xx]=find(Sigma~=0);
    chS=chol(Sigma(nzsig,nzsig))';
    for uu=1:10 % with parfor the seed gets messed up
        initval=zeros(size(Sigma,1),periods);
        inno=randn(length(nzsig),periods);
        initval(nzsig,:)=chS*(inno-repmat((mean(inno,2)),1,periods));
        vary=zeros(size(PHI,1),periods);
        for jj=2:periods
            vary(:,jj)=PHI*vary(:,jj-1)+M*initval(:,jj);
        end
        c=Pi*vary(:,drop+1:end);
%         ctild=c-c*HPM';
        covyHP(:,:,uu)=(c);
        
    end
    

end
return

%tot variance;
