    function Sh_end=getNW(Uh,nT,nwLags,S0)
    Sh_end=S0;
    nwWeights=(nwLags+1-(1:nwLags))./(nwLags+1);
    %HAC correction: posterior scale (mean of posterior IW distribution)
    for l=1:nwLags

        Gammal=(Uh(l+1:nT,:)'*Uh(1:nT-l,:))/(nT-l);
        
        Sh_end=Sh_end+nwWeights(l)*(Gammal+Gammal');

    end