function out = weighted_ma(shocks,M)

[dnobs, n_shocks] = size(shocks);

iper=find(isnan(sum(shocks,2))==0);
istart=iper(1);
iend=iper(end);

out = nan(n_shocks,dnobs);
for t = istart:iend
    if t==istart
        out(:,t) =  shocks(t,:)';
    else
        out(:,t) =  M*out(:,t-1) + shocks(t,:)';
    end
end

out = out';

end