function res = PCA_CK(rt,gt,pmax,q,w)
%% INPUT
% rt          is n by T matrix
% gt          is d by T factor proxies
% pmax        is the upperbound of the number of latent factors
% q           is # of lags used in Newy-West standard errors
% w           is a scaler that determines the weight of the return component
%% OUTPUT

% Gammahat    is (d+1) by pmax matrix of risk premia estimates


%% INITIALIZATION

    T  =  size(rt,2);
    n  =  size(rt,1);
    d  =  size(gt,1);

    Gammahat              =      nan(d+1,pmax); 

%% ESTIMATION

    rtbar                 =      rt - repmat(mean(rt,2),1,T);% 
    rbar                  =      mean(rt,2);
    gtbar                 =      gt - repmat(mean(gt,2),1,T);% 
    

        [tempV,tempD]         =      eig(rtbar*rtbar'/n/T+rbar*rbar'/n*w);  
        [tempeig,tempidx]     =      sort(real(diag(tempD)),'descend');
        eigvec                =      tempV(:,tempidx);    

       
   
    for phat  = 1 : pmax   % at least one factor up to pmax
        

            betahat           =       n^0.5  * eigvec(:,1:phat);%n^0.5  * eigvec(:,(n-phat+1):n);
            vhat              =       n^(-1) *  betahat'*rtbar;
      
       
    % step 2: FM with PCs
    
            Gammatilde            =       [0;inv(betahat'*betahat) * betahat'* mean(rt,2)];

    % step 3: TS Regression
    
            etahat                =       gtbar   *  vhat' * inv(vhat*vhat');   
       
        
    % step 4: Combine TS and CS estimates
        
            Gammahat(:,phat)      =       [Gammatilde(1);etahat*Gammatilde(2:end)];   
    

    end
    

%% OUTPUT


    res.Gammahat       = Gammahat;

    
  