% variance decomposition
function [varY_]=var_dec_symb(PHI,Pi,Msh,Sigma,nsh);

szphi=size(PHI,1);
szY=size(Pi,1);

hor=150;
N=PHI(1:nsh,1:nsh);
PHI1=PHI(nsh+1:end,nsh+1:end);
PHI2=PHI(nsh+1:end,1:nsh);
M=[Msh;
    zeros(szphi-nsh,nsh)];
varS=zeros(szphi,szphi);
varSC=zeros(szphi,szphi,7);
varY=zeros(szY,7);
zz=1;
disp(upper('Variance decomposition computed at horizons: 1q,2q,4q,8q,12q,150q'));
for jj=1:hor
   % if jj==1
    varS=PHI*varS*PHI'+M*Sigma*M';
   % else
   %      varS=PHI*varS*PHI';
   % end
        if jj==1|jj==2|jj==4|jj==8|jj==12|jj==hor-1|jj==hor
        varSC(:,:,zz)=varS;

        varY(:,zz)=diag(Pi*varSC(:,:,zz)*Pi');
                zz=zz+1;
        end
end

varY(varY==0)=NaN;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for uu=1:nsh
eval(['varS_',num2str(uu),'=zeros(szphi,szphi);']);
eval(['varSC_',num2str(uu),'=zeros(szphi,szphi,7);']);
varY_=zeros(szY,7,nsh);
end

for uu=1:nsh
    zz=1;
    Sigma_0=zeros(nsh,nsh);
    Sigma_0(uu,uu)=Sigma(uu,uu);
for jj=1:hor
    %if jj==1
    eval(['varS_',num2str(uu),'=PHI*varS_',num2str(uu),'*PHI''+M*Sigma_0*M'';']);
    %else
    %eval(['varS_',num2str(uu),'=PHI*varS_',num2str(uu),'*PHI'';']);    
    %end
        if jj==1|jj==2|jj==4|jj==8|jj==12|jj==hor-1|jj==hor
         
        eval(['varSC_',num2str(uu),'(:,:,zz)=varS_',num2str(uu),';']);
        eval(['tmp=varSC_',num2str(uu),'(:,:,zz);']);
        varY_(:,zz,uu)=diag(Pi*tmp*Pi')./varY(:,zz);
        zz=zz+1;
        end
end
end

if abs(1-max(max(abs(sum(varY_,3)))))>1e-5;
    disp('Variance decomposition does not add up to one')
end
return
%tot variance;
