function [Phi_new, Phi_sample] = m_var(theta, M, p, Y, flag_mean)
        % find mean/median of OLS when DGP is VAR(p)
        % works for VAR(1) and VAR(2) only
        % inputs:
        %  theta - vector of parameters determining Phi_1 (and Phi_2)
        %  M - number of Monte Carlo replications
        %  p - order of the VAR
        %  Y - data, T x k matrix
        %  flag_mean - flag whether mean (TRUE) or median (FALSE) is to be returned
        
        if nargin<5; flag_mean = true; end;
        
        [T, k] = size(Y);
        
        if (p==1)
            % get coefficient matrix
            Phi_tilde =  reshape(theta, k, k);
            
            % simulate M datasets
            X_sim = gen_var(Phi_tilde, M, Y, p);
            
            % estimate VAR(1) on each series
            theta_new_i = zeros(M, k^2);
            for m=1:M
                Phi_new = est_var(X_sim(:,:,m)', 1, true, false);
                theta_new_i(m,:) = Phi_new(:);
            end
            
            if (flag_mean)
                % find mean of OLS estimates
                Phi_new = reshape( mean(theta_new_i), k, k );
            else
                % find median of OLS estimates
                Phi_new = reshape( median(theta_new_i), k, k );
            end
            
            Phi_sample = theta_new_i;
        elseif (p==2)
            % get coefficient matrix
            Phi_tilde =  reshape(theta, k, 2*k);
            
            % simulate M datasets
            X_sim = gen_var(Phi_tilde, M, Y, p);
            
            % estimate VAR(2) on each series
            theta_new_i = zeros(M, 2*k^2);
            for m=1:M
                F_new = est_var(X_sim(:,:,m)', 2, true, false);
                tmp = F_new(1:k,:);
                theta_new_i(m,:) = tmp(:);
            end
            
            if (flag_mean)
                % find mean of OLS estimates
                Phi_new = reshape( mean(theta_new_i), k, 2*k );
            else
                % find median of OLS estimates
                Phi_new = reshape( median(theta_new_i), k, 2*k );
            end
        else
            error('not yet implemented for order>2');
        end
        
    end
