for gg=1:NoGroups
    
    Y_b(:,t_idx,1==IndexM(:,gg)) = reshape( A_IV{gg}' * squeeze(X_b(1:px,t_idx,1==IndexM(:,gg))) , [p,1,sum(IndexM(:,gg))]) ...
        + reshape(DATA.MSA_FE_res(:, 1==IndexM(:,gg)), [p, 1,sum(IndexM(:,gg))]) ...
        + repmat(DATA.Year_FE_res(:, t_idx,gg), [1, 1,sum(IndexM(:,gg))]) ...
        + reshape(ErrorTerm(1:p,1==IndexM(:,gg)), [p, 1,sum(IndexM(:,gg))]) ; % error term includes spatial correlation


    Z_b(:,t_idx,1==IndexM(:,gg)) = ...
        reshape(DATA.MSA_FE_Z(:,1==IndexM(:,gg)), [pz, 1,sum(IndexM(:,gg))])  ...
        + repmat(DATA.Year_FE_Z(:, t_idx,gg), [1, 1,sum(IndexM(:,gg))])  ...
        + reshape(ErrorTerm(p+(1:pz),1==IndexM(:,gg)), [pz, 1,sum(IndexM(:,gg))]) ;

    % Challenge: Periphery RHS depends on current Y_b. Thus cannot update
    % one period in advance as for core VAR
    if t_idx > 1
        for p_idx=1:pp
            pX_b(1:k, t_idx,:,p_idx)=[pY_b(p_idx,t_idx-1,:) ; pX_b(1:k-1, t_idx-1,:,p_idx)]; % lags of periphery variable p_idx
            % cannot update current 
            pX_b(k+(1:p*(1+k)), t_idx,:,p_idx)=[Y_b(:,t_idx,:) ; pX_b(k+(1:p*k), t_idx-1,:,p_idx)]; % current value and lags of core variables
            % SANITY CHECK: 
            % squeeze(DATA.pX(k+(1:p*(1+k)), t_idx,15:16,p_idx)-[DATA.Y(:,t_idx,15:16) ;DATA.pX(k+p+(1:p*(k)), t_idx,15:16,p_idx)])
            if DO_TEST
                squeeze(pX_b(:, t_idx,15,p_idx)-pX(:, t_idx,15,p_idx))'
            end;
        end;
    end
    
    for p_idx=1:pp
        pY_b(p_idx,t_idx,1==IndexM(:,gg)) = reshape( Ap_IV{gg,p_idx}' * squeeze(pX_b(:,t_idx,1==IndexM(:,gg),p_idx)) , [1,1,sum(IndexM(:,gg))]) ...
            + reshape(DATA.MSA_FE_pres(p_idx, 1==IndexM(:,gg)), [1, 1,sum(IndexM(:,gg))]) ...
            + repmat(DATA.Year_FE_pres(p_idx, t_idx,gg), [1, 1,sum(IndexM(:,gg))]) ...
            + reshape(pErrorTerm(p_idx,1==IndexM(:,gg)), [1, 1,sum(IndexM(:,gg))]) ;
    end;
    % % Quality check!            
    % %             Y_b(:,t_idx,89)-DATA.Y(:,t_idx,89)
    % Y_b(:,:,89)-DATA.Y(:,:,89)
    % %             Z_b(:,t_idx,89)-DATA.Z(:,t_idx,89) - reshape(DATA.MSA_FE_Z(:, 89), [pz, 1, 1]) - repmat(DATA.Year_FE_Z(:, t_idx), [1, 1, 1]) 
    % Z_b(:,3,89)-DATA.Z(:,3,89) - reshape(DATA.MSA_FE_Z(:, 89), [pz, 1, 1]) - repmat(DATA.Year_FE_Z(:, 3), [1, 1, 1]) 
    % Z_b(:,13,89)-DATA.Z(:,13,89) - reshape(DATA.MSA_FE_Z(:, 89), [pz, 1, 1]) - repmat(DATA.Year_FE_Z(:, 13), [1, 1, 1]) 
    
    if DO_TEST
        squeeze(Y_b(:,t_idx,87)-DATA.Y(:,t_idx,87))'
        squeeze(pY_b(:,t_idx,87)-DATA.pY(:,t_idx,87))'
    end;
end;

% update the RHS for the next time period
if t_idx < T-MissingObs-k
    X_b(:, t_idx+1,:)=[Y_b(:,t_idx,:) ; X_b(1:p*(k-1), t_idx,:)];
    % % Quality check!            
    % DATA.X(:,2,1) - [ DATA.Y(1:p,1,1); DATA.X(1:p,1,1)] == 0;
end;
