function res = PCAFM_fast(rt,gt,pmax,q)
% rt is n by T matrix
% gt is d by T factor proxies
% pmax is the number of latent factors
% q is # of lags used in Newy-West standard errors

%% INITIALIZATION

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

    %% ESTIMATION

    rtbar                 =      rt - repmat(mean(rt,2),1,T);% 
    gtbar                 =      gt - repmat(mean(gt,2),1,T);% 
    
    [U,S,V] = svd(rtbar/sqrt(n)/sqrt(T),'econ');    
    
    vhat              =       T^0.5  * V(:,1:pmax)';
    betahat           =       T^(-1) * rtbar * vhat';
    
    Gammatilde            =       inv([ones(n,1),betahat]'*[ones(n,1),betahat]) * [ones(n,1),betahat]'* mean(rt,2);
    
    Gammatilde_nozero     =      [0; inv(betahat'*betahat) * betahat'* mean(rt,2)];

    etahat                =       gtbar   *  vhat' * inv(vhat*vhat');   
        
        
    Gammahat      =       [Gammatilde(1);etahat*Gammatilde(2:end)];   
    Gammahat_nozero=      [Gammatilde_nozero(1);etahat*Gammatilde_nozero(2:end)];   
    
    res=[Gammahat Gammahat_nozero];
    
  