function [fac,lam,eigval,ssr] = pc_factor(x,k)
%    Compute principal components estimates of factor model
%
%    Input:
%     x = txn matrix
%     k = number of factors
%
%     Model
%
%     x(it) = lam(i)'f(t) + u(it)
%
%     or
%
%     x = f*lam' + u
%
%
%     Output:
%
%     f = txk matrix of factors
%     lam = n*k matrix of factor loadings
%     eigval = kx1 vector of eigenvalues of x*x' (or x'x) ...
%     ssr = sum of squared u's
%
%     Normalization:
%      f is normalized so that each column has std dev = 1, thus F'F = t*I(k)
%
%     Calculation note:
%      Calculations are speeded by using the smaller of x*x' or x'x
[t,n]=size(x);

if k > 0
    if n < t
        xx=x'*x;
        [ve,va]=eig(xx);
        va=diag(va);
        
        
        va=flipud(va); %sort eigen value from large to small.
        ve=fliplr(ve); %sort the corresponding eigen vectors.
        eigval=va(1:k); %keep the first k eigen value.
        lam=ve(:,1:k); %keep the first k eigen vectors.
        fac=x*lam;
    else
        xx=x*x';
        [ve,va]=eig(xx);
        va=diag(va);
        
        va=flipud(va); %sort eigen value from large to small.
        ve=fliplr(ve); %sort the corresponding eigen vectors.
        eigval=va(1:k);
        fac=ve(:,1:k);
    end
    
    %@ Normalize @
    sfac=sqrt(mean(fac.^2));
    fac=fac./sfac;
    lam=(x'*fac)/t;  % Note fac'fac = t @
    ssr=sum(va)-sum(eigval);
else
    ssr=sum(sum(x.^2));
    lam=nan;
    eigval=nan;
    fac=nan;
end


end
