function rra = func_rraff_weightchange(R,G,K,S2)
% using R=a+BG+u to extract K factors
% R: T*N, G: T*L
% output: T*K factors with RRA

[T,L] = size(G);
N = size(R,2);
Z = [ones(T,1) G];
M = L+1;

X = ones(T,1);

W1 = eye(N);   
W2 = eye(M);  
P0 = Z*W2*Z';
P = P0 -P0*X*inv(X'*P0*X)*X'*P0;
Q = G'*P*G/T^2;
A = (Q^(-1/2))'*(G'*P*R/T^2)*W1*(G'*P*R/T^2)'*Q^(-1/2);
[E,v] = eig(A);
[v,ind] = sort(diag(v),'descend');
E = E(:,ind);
Phi = Q^(-1/2)*E(:,1:K);
Gstar=G*Phi;
Beta = (Gstar'*P*Gstar)^(-1)*Gstar'*P*R;
Theta = Phi*Beta;
Alpha = (X'*P0*X)^(-1)*X'*P0*(R-G*Theta);
Alpha = Alpha';
U = R-X*Alpha'-G*Theta;
S1 = diag(diag(U'*U/T));
%S2 = diag(diag(Z'*Z/T));


W1 = inv(S1); 
W2 = inv(S2);
P0 = Z*W2*Z';
P = P0 -P0*X*inv(X'*P0*X)*X'*P0;
Q = G'*P*G/T^2;
A = (Q^(-1/2))'*(G'*P*R/T^2)*W1*(G'*P*R/T^2)'*Q^(-1/2);
[E,v] = eig(A);
[v,ind] = sort(diag(v),'descend');
E = E(:,ind);
Phi = Q^(-1/2)*E(:,1:K);
Gstar=G*Phi;

rra = Gstar;
 

%[~,~,F] = plsregress(G,R,K);
%plsf = F;
 
