function charPolyRoots = CharPolyRoots(b, alpha, beta, Gamma, printRoots)

% function charPolyRoots = CharPolyRoots(b, alpha, beta, Gamma, printRoots)
% Lee Morin and Morten Nielsen
% August 23, 2011
%
% CharPolyRoots calculates the roots of the characteristic polynomial and plots them
%       with the unit circle transformed for the fractional model, see Johansen (2008).
% 
% input = scalar b is the fractional differencing parameter.
%       matrices alpha, beta and Gamma are the parameters for the
%           fractional cointegrated VAR model.
%       matrix Gamma is a p x p*k matrix of coefficients in the lag
%           polynomial.
%       printRoots is an indicator to determine if a graph of the roots
%           will be produced.
% 
% output = complex vector charPolyRoots with the roots of the characteristic polynomial.
% 
% No dependencies.
% 
% Note: The roots are calculated from the companion form of the VAR, where the roots 
%		are given as the eigenvalues of the coefficient matrix.
%______________________________________________________


% Get dimension of system from alpha.
p = size(alpha, 1);
% Get cointegrating rank from alpha.
r = size(alpha, 2);
% Get lag length from Gamma.
k = floor(size(Gamma,2)/p + 0.1);


% First construct the coefficient matrix for the companion form of the VAR. 
PiStar = eye(p);
if r > 0
    PiStar = PiStar + alpha*beta';
end

if k > 0
    Gamma1 = Gamma(:, 1 : p );
    PiStar = PiStar + Gamma1;
end
for i = 2:k
    
    Gammai = Gamma(:, (i-1)*p + 1 : i*p );
    GammaiMinus1 = Gamma(:, (i-2)*p + 1 : (i-1)*p );
    
    PiStar = [ PiStar (Gammai - GammaiMinus1) ];
    
end
if k > 0
    Gammak = Gamma(:, (k-1)*p + 1 : k*p );
    PiStar = [ PiStar  ( - Gammak ) ];
end

% Pad with an identity for the transition of the lagged variables.
PiStar = [  PiStar;...
            eye(p*k) zeros( p*k, p*(k>0) )  ];

        
% The roots are then the eigenvalues of the matrix PiStar.
charPolyRoots = eig(PiStar);
charPolyRoots = sort( charPolyRoots , 'descend');

% Generate graph depending on the indicator printRoots.
if printRoots
    
    % Now calculate the line for the transformed unit circle.
    % First do the negative half.
    unitCircle = ( pi : - 0.001 : 0)';
    psi = - (pi - unitCircle)/2;
    unitCircleX = cos( - unitCircle);
    unitCircleY = sin( - unitCircle);
    transformedUnitCircleX = [ (1 - (2*cos(psi)).^b.*cos(b*psi)) ];
    transformedUnitCircleY = [ (    (2*cos(psi)).^b.*sin(b*psi)) ];
    % Then do the positive half.
    unitCircle = (0 : 0.001 : pi)';
    psi = (pi - unitCircle)/2;
    unitCircleX = [ unitCircleX; cos(unitCircle) ];
    unitCircleY = [ unitCircleY; sin(unitCircle) ];
    transformedUnitCircleX = [ transformedUnitCircleX;   (1 - (2*cos(psi)).^b.*cos(b*psi)) ];
    transformedUnitCircleY = [ transformedUnitCircleY;   (    (2*cos(psi)).^b.*sin(b*psi)) ];
    
    % Plot the unit circle and its image under the mapping
    % along with the roots of the characterisitc polynomial.
    plot(transformedUnitCircleX, transformedUnitCircleY);
    hold on;
    plot(unitCircleX, unitCircleY);
    scatter(real(charPolyRoots), imag(charPolyRoots));
    maxXYaxis = max( [ transformedUnitCircleX; unitCircleX;...
                       transformedUnitCircleY; unitCircleY ] );
    minXYaxis = min( [ transformedUnitCircleX; unitCircleX;...
                       transformedUnitCircleY; unitCircleY ] );
    maxXYaxis = max( maxXYaxis, -minXYaxis );
    axis( 2*[ -maxXYaxis; maxXYaxis; -maxXYaxis; maxXYaxis; ] );
    hold off;
    title('Roots of the Characteristic Polynomial with the Image of the Unit Circle');
    
end


% end