% Listing 14.8 of the book for SVD commands
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% September 04, 2015
% Freely distributed under no warranty

[N,p] = size(Phi);
[U,S,V] = svd(Phi);  % full factorization
Phiplus = V*pinv(S)*U'	
Phiplus = pinv(Phi); % both give pseudo-inverse
Theta = Phiplus*Y;
pred_err = sum(U(:,p+1:end).*repmat(Y,1,N-p),1); 
pred_err = Y - Phi*Theta; % both give residuals
[U,S,V] = svd(Phi,0); 
[U,S,V] = svd(Phi,'econ'); % try with economy factorization