function  [y_hat]=make_hstep_forecast_MAIAR(gamma,B0,A,G_inv,sigmas,y,qp,H)

n=size(y,1);
q=qp(1); p=qp(2);
% INPUT :  y is nxT,  H is the horizon
Gammas=zeros(n,n,q);
Gammas(logical(repmat(eye(n),[1,1,q])))=gamma;

% preallocation
x_bul_h_p=ones(n,p,H); 
x_bul_h_q=ones(n,q,H); 
yy=y; %nxT
y_hat=zeros(n,H); 


rnsims=randn(n,H);
u= (G_inv * diag(sigmas))*rnsims;
    
    for h=1:H        
       x_bul_h_p(:,:,h)= yy(:,end:-1:end-p+1);
       x_bul_h_q(:,:,h)= yy(:,end:-1:end-q+1);
       Z_h=vec(B0*x_bul_h_p(:,:,h));
       y_hat(:,h)=A*Z_h+sum(mtimesx(Gammas, reshape(x_bul_h_q(:,:,h),[n,1,q])  ),3)+u(:,h);
       yy=[yy, y_hat(:,h) ];
    end

    


end