function out=get_errors(qsave,Ssave,N,L,mature,dt,yields,W)
   NFACT=N;
   T=rows(yields);
    bold=squeeze(mean(qsave));
    sigma_cP=squeeze(mean(Ssave));
      %extract Q parameters
  K1Q_X = (bold(1:NFACT));
    rinfQ = bold(NFACT+1);
    omega=sigma_cP*sigma_cP';

cP = yields*W'; % (T+L)*N, cP stands for math caligraphic P.
[BcP, AcP, AX, BX, DUMMY,penalty] = ...
    jszLoadings(W, diag(K1Q_X-1), rinfQ, ...
    omega(1:N,1:N), mature, dt);
[K0Q_cP, K1Q_cP, rho0_cP, rho1_cP] = ...
    jszRotation(W, diag(K1Q_X-1), rinfQ, dt, [], [], BX, AX);
  
  % term premium
 yields_m =  ones(T,1)*AcP + cP*BcP; % (T+1)*J, model-implied yields

out=yields-yields_m;
% square_orthogonal_yield_errors = (out).^2; % T*J, but N-dimensional projection onto W is always 0, so effectively (J-N) dimensional
%  sigma_e = sqrt( sum(square_orthogonal_yield_errors(:))/(T*(J-N)) );
