    function Phi_tilde = shrink_phi(Phi_tilde, Phi_hat, ev_restr, trace)
    if nargin<4; trace = 1; end;
    % stationarity adjustment
    % see Kilian (1998, REStat) "Small-sample confidence intervals for
    % impulse response functions"
        % impose eigenvalue restriction
        ev_bc = max(abs(eig(Phi_tilde)));   % get scalar largest absolute eigenvalues
        ev_ols = max(abs(eig(Phi_hat)));
        
        % impose restriction:  if ev_bc is larger than ev_restr
        if (ev_bc>ev_restr)
            if trace>0;
                fprintf('*** eigenvalue constraint binding ***\n');
                fprintf('*** constraint: %8.6f \n', ev_restr);
            end
            if (ev_ols<ev_restr)   % and if Phi_hat actually satisfies eigenvalue restr
                % then shrink Phi_tilde to Phi_hat
                Phi_diff = Phi_hat - Phi_tilde;   % difference/bias adjustment -- this is what we shrink to zero
                delta = 1 ;
                % as long as restriction is not satisfied, go through the following loop
                while (ev_bc>ev_restr)
                    delta = delta - .01;      % reduce scaling factor
                    Phi_diff = delta*Phi_diff;     % shrink Phi_diff to zero
                    Phi_tilde = Phi_hat - Phi_diff;    % try this value of Phi_diff
                    ev_bc = max(abs(eig(Phi_tilde)));   % calculate largest eigenvalue
                end
            else
                % OLS estimates do not satisfy restriction
                %warning('ev_ols > ev_restr -- setting Phi_bc equal to Phi_ols'); %#ok<WNTAG>
                Phi_tilde = Phi_hat;
            end
            if trace>0;
                fprintf('largest absolute eigenvalue after imposing eigenvalue restriction:  %8.6f \n', ev_bc);
            end;
        end
    end
