function [x,C,e,Sgm,sing,sample,f] = pc(y,crit,method)
% pc  Principal components of input series.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

%**************************************************************************

% crit = [power,maxnumber].
if length(crit) == 1
   crit = [crit,Inf];
end

[ny,nper] = size(y);

if strcmpi(method,'auto')
   if nper > ny
      method = 1;
   else
      method = 2;
   end
end

% Non-NaN sample.
nanindex = isnan(y);
sample = all(~nanindex,1);
nobs = sum(sample);
if nobs == 2
   FAVAR.error(1);
end

% Covariance matrix of input series.
% The matrix is needed whatever method.
ycov = y(:,sample)*y(:,sample)'/nobs;
if any(any(isinf(ycov) | isnan(ycov)))
   FAVAR.error(2);
end

if method == 1
   Q = ycov;
   n = ny;
else
   Q = y(:,sample)'*y(:,sample)/nobs;
   n = nobs;
end

[U,sing] = svd(Q);
sing = diag(sing(1:n,1:n));
cumsing = cumsum(sing);
cumsing = cumsing/cumsing(end);
r = min([find(cumsing >= crit(1),1),crit(2)]);
x = nan(r,nper);
e = nan(ny,nper);
f = nan(ny,nper,r);

% y = C*x + e;
if method == 1
   C = U(:,1:r);
   x(:,sample) = C.'*y(:,sample);

   % Normalise stdevs of the factors to 1.
   a = sqrt(sing(1:r));
   C = C*diag(a);
   x(:,sample) = diag(1./a)*x(:,sample);

else
   x(:,sample) = U(:,1:r).';
   C = y(:,sample)*x(:,sample).';

  % Normalise stdevs of the factors to 1.
   C = C/sqrt(nobs);
   x = sqrt(nobs)*x;
end
 
% Contributions of input data to the estimated factors.
repeat = ones([1,nobs]);
for i = 1 : r
   Ci = C(:,i*repeat);
   f(:,sample,i) = Ci.*y(:,sample) / sing(i);
end

% The number of factors equals the number of input series.
if r == ny
   e(:,sample) = 0;
   Sgm = zeros(ny);
else
   e(:,sample) = y(:,sample) - C*x(:,sample);
   Sgm = ycov - C*C';
end

end
