function [ betaStar, alphaHat, OmegaHat ]...
    = GetParamsSwitching(beta0, S00, S01, S11, T, p, r, EstnOptions)

% function [ betaStar, alphaHat, OmegaHat ]
%       = GetParamsSwitching(beta0, S00, S01, S11, T, p, r, EstnOptions)
% Lee Morin & Morten Nielsen
% August 26, 2011
%
% 
% GetParamsSwitching(beta0, S00, S01, S11, T, p, r, EstnOptions) returns the 
%		parameter estimates of the fractionally cointegrated VAR model using 
%		regression and reduced rank regression, by way of Johansen's switching algorithm.
% 
% input = matrix beta0 of starting values in the switching algorithm.
%       matrices S00, S01 and S11 of inner products, see Johansen (1996).
%       scalar T is the number of observations used for estimation.
%       scalar p is the dimension of the system.
%       scalar r is the cointegrating rank.
%       cell array of options EstnOptions, see DefaultEstnOptions(db0, k, r, p), 
%			RestrictEstnOptions(db0, k, r, p, defaultFCVARoptions), 
%			and AdjustEstnOptions(db0, k, r, p, initialFCVARoptions).
% 
% output = matrices betaStar, alphaHat, OmegaHat of parameter estimates.
%
% No function dependencies.
%______________________________________________________



% Determine desired set of options.
deterministics = EstnOptions(3);
% Define p1 to be number of rows of betaStar.
p1 = p + strcmp(deterministics, 'restricted constant');

% Determine the restrictions on the cointegrating relations.
% EstnOptions = { N, M, deterministics, R_Gamma, r_Gamma, [], [], [],
%                   nCols, A, H, a, b, [], Abar, Aperp, replaceCols, keepCols,
%                   [], [], [], print2screen, printGammas, printRoots,
%                   restrictFminOptions }
nCols = EstnOptions{9};
nRestns = length(nCols);    % Note: nRestns = 0 when nCols = [].
A = EstnOptions{10};
H = EstnOptions{11};
a = EstnOptions{12};
b = EstnOptions{13};
Abar = EstnOptions{15};
Aperp = EstnOptions{16};
replaceCols = EstnOptions{17};
keepCols = EstnOptions{18};


% Determine the optimization settings for the switching algorithm.
restrictFminOptions = EstnOptions{end};
MaxFunEvals = restrictFminOptions.MaxFunEvals;
TolX = restrictFminOptions.TolX;
TolFun = restrictFminOptions.TolFun;


% Initialize betaStar.
betaStar = beta0;
for j = 1:nRestns
    if ~isempty(b{j}) 
        % This(These) beta column(s) is (are) known.
        betaR = b{j}; 
    elseif size(H{j}, 2) < p1
        % Estimate current beta (Note: Uses entire unrestricted beta0 matrix).
        [V,D] = eig( beta0'*H{j}*inv(H{j}'*H{j})*H{j}'*beta0, beta0'*beta0 );
        VpD = [ V' diag(D) ];
        p11 = size(VpD, 1);
        V = sortrows( VpD, p11 + 1 )';
        V = V( 1 : p11 , :);
        % Assemble updated beta0 matrix.
        vj = V( 1:p11, p11 : -1 : p11-nCols(j)+1 );
        betaR = beta0*vj;
    else
        % No initialization necessary when H{j} = I and b{j} is empty.
        betaR = beta0(:, replaceCols{j} );
    end
    betaStar(:, replaceCols{j} ) = betaR;
end

% Initialize alphaHat.
alphaHat = S01*betaStar*inv(betaStar'*S11*betaStar);
for j = 1:nRestns
    if ~isempty(a{j})
        alphaHat(:, replaceCols{j} ) = a{j};
    end
end


% Initialize stopping rules for switching algorithm.
i = 0;
% betaStar = beta0; % Already set in initialization loop.
betaStarLast = betaStar + 10*ones(size(betaStar));
logLike = 0;
logLikeLast = 10;

% Iterate estimation of beta and alpha.
while (  (i < MaxFunEvals)... 
         && (norm(betaStar - betaStarLast) > TolX)...
         && (abs(logLike - logLikeLast) > TolFun)...
      )
    
    % Update values used as stopping rules.
    i = i+1;
    betaStarLast = betaStar;
    logLikeLast = logLike;
    
    for j = 1:nRestns
    
        % Reset the inner product matrices.
        S00ab = S00;
        S01ab = S01;
        S11ab = S11;
        
        % Create Aj, a temporary A{j}, that will be replaced if:
        %       (i)  There are other submatrices of alpha and beta that will
        %            be considered known (then Aj = alphaKperp).
        %  AND  (ii) There is NO linear restriction specified in A{j}.
        %  AND  (iii)This submatrix of alpha is NOT specified as known in a{j}.
        Aj = A{j};
        

        % Adjust Sij matrices for other columns known assumed known, if any.
        if nRestns > 1
            % There are multiple submatrices of alpha and beta.

            % Create matrix of other known betas, called betaK.
            betaK = betaStar(:, keepCols{j} ); % Only used if nRestns > 1.
            % This is the matrix of other betas assumed known for now.
            % Same for alpha.
            alphaK = alphaHat(:, keepCols{j});
            
            % Adjust Sij's to account for the parameters assumed known.
            S00temp = S00ab;
            S01temp = S01ab;
            S00ab = S00temp - S01temp*betaK*alphaK' - alphaK*betaK'*S01temp' + alphaK*betaK'*S11ab*betaK*alphaK';
            S01ab = S01temp - alphaK*betaK'*S11ab;
            
            % Update Aj matrix to be alphaKperp if there are no
            %    restrictions on alpha{j} or beta{j}.
            if  (   ( size(H{j}, 2) == p1 ) && ( isempty(b{j}) )...
                &&  ( size(A{j}, 2) == p ) && ( isempty(a{j}) )  )
                
                Aj = null( alphaK' ); % = alphaKperp.
                Abar{j} = Aj*inv(Aj'*Aj);
                Aperp{j} = null( Aj' );
                
            end
            
        end
        
        
        % Correct Sij matrices for linear restrictions or known values for current alphaj.
        % In any other case, will have AjPerp = [].
        if ~isempty(Aperp{j})
            % Linear restrictions or known values for alpha (FWL out the extra regressors).
            S00temp = S00ab;
            S01temp = S01ab;
            S11temp = S11ab;
            S00ab = S00temp - S00temp*Aperp{j}*inv(Aperp{j}'*S00temp*Aperp{j})*Aperp{j}'*S00temp;
            S01ab = S01temp - S00temp*Aperp{j}*inv(Aperp{j}'*S00temp*Aperp{j})*Aperp{j}'*S01temp;
            S11ab = S11temp - S01temp'*Aperp{j}*inv(Aperp{j}'*S00temp*Aperp{j})*Aperp{j}'*S01temp;
        end

        % Matrix of cross terms is symmetric.
        S10ab = S01ab'; 
        
        
        
        % Finally, it's time for estimation.
        if ( isempty(b{j}) && isempty(a{j}) ) 
            % Beta and alpha estimated without further restriction.
            
            % Ensure exact symmetry of the matrix B in eig(A, B).
            HSH = H{j}'* S11ab * H{j};
            HSH = ( HSH + HSH' )/2;
            
            [V,D] = eig( H{j}'* S10ab * Abar{j}...
                            *inv(Abar{j}' *S00ab * Abar{j})...
                                * Abar{j}' *S01ab * H{j},...
                                    HSH); 
            VpD = [ V' diag(D) ];
            p11 = size(VpD, 1);
            V = sortrows( VpD, p11 + 1 )';
            V = V( 1 : p11 , :);
            
            % Assemble updated beta0 matrix.
            phij = V( 1:p11, p11 : -1 : p11-nCols(j)+1 );
            betaR = H{j}*phij;
            betaStar(:, replaceCols{j} ) = H{j}*phij; % = betaR.
            
            alphaHat(:, replaceCols{j} ) = Aj * Abar{j}'...
                                            * S01ab * betaR * inv(betaR' * S11ab * betaR);
            
        elseif ( ~isempty(b{j}) && isempty(a{j}) )  
            % This (These) beta column(s) is (are) known (but alpha unknown).
            
            betaR = b{j};
            % Do nothing to betaStar since this known beta is already initialized correctly.
            alphaHat(:, replaceCols{j} ) = Aj * Abar{j}'...
                                            * S01ab * betaR * inv(betaR' * S11ab * betaR);
            % Note: Note that Abar{j} accounts for the possibility that a
            %   linear restriction is imposed on this alphaR (and Abar{j} = I if A{j} = I).
            
        elseif ( isempty(b{j}) && ~isempty(a{j}) )  
            % This (These) alpha column(s) is (are) known (but beta unknown).
            
            betaStar(:, replaceCols{j} ) = H{j} * inv(H{j}'*S11ab*H{j}) * H{j}'*S10ab * Abar{j};
            
            % Do nothing to alphaHat since this known alpha is already initialized correctly.
        
        end
        % Else both beta and alpha are known and are already set from initialization.
        
    end
    
    
    % Obtain estimate of Omega.
    OmegaHat = S00 - S01*betaStar*alphaHat' - alphaHat*betaStar'*S01' + alphaHat*betaStar'*S11*betaStar*alphaHat';
    
    % Calculate current value of likelihood.
    logLike = - T*p/2*( log(2*pi) + 1)  - T/2*log(det(OmegaHat));
    
end


% end


