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

    T  =  size(rt,2);
    n  =  size(rt,1);
    d  =  size(gt,1);
    p  =  size(vt,1);
    
    
    Gammahat              =      nan(d+1,1); 
 %   avarhat               =      nan(d,1);  
    
    rtbar                 =      rt - repmat(mean(rt,2),1,T);% 
    vtbar                 =      vt - repmat(mean(vt,2),1,T);% 
    gtbar                 =      gt - repmat(mean(gt,2),1,T);% 
        
        
        % step 1: TS Regression
   
        betahat = rtbar * vtbar'*inv(vtbar*vtbar');
      
        % step 2: FM        
        Gammatilde  =  inv([ones(n,1),betahat]'*[ones(n,1),betahat]) * [ones(n,1),betahat]'* mean(rt,2);
        
        Gammatilde_nozero  =  inv([betahat]'*[betahat]) * [betahat]'* mean(rt,2);
        
        % step 3: Regression
        etahat  = gtbar * vtbar' * inv(vtbar*vtbar');  
        what    = gtbar  - etahat*vtbar;
        gthat = etahat*vtbar;
        % step 4: Adjustment
        Gammahat(:,1) = [Gammatilde(1);etahat*Gammatilde(2:end)];
        
        Gammahat_nozero(:,1) = [0;etahat*Gammatilde_nozero];
     
        % step 5: R2
        
        %gtbetahat =  rtbar * gtbar'*inv(gtbar*gtbar');
        %gtGamma  =  inv([ones(n,1),gtbetahat]'*[ones(n,1),gtbetahat]) * [ones(n,1),gtbetahat]'* mean(rt,2);
        eps = mean(rt,2) - [ones(n,1),betahat]*Gammatilde;
        R2  = 1 - sum(eps.^2)/sum((mean(rt,2) - mean(mean(rt,2))).^2); 
        
        
        gtbetahat =  rtbar * gthat'*inv(gthat*gthat');
        gtGamma  =  inv([ones(n,1),gtbetahat]'*[ones(n,1),gtbetahat]) * [ones(n,1),gtbetahat]'* mean(rt,2);
        gteps = mean(rt,2) - [ones(n,1),gtbetahat]*gtGamma;
        R2ghat  = 1 - sum(gteps.^2)/sum((mean(rt,2) - mean(mean(rt,2))).^2); 

        
        gtbetahat =  rtbar * gtbar'*inv(gtbar*gtbar');
        gtGamma  =  inv([ones(n,1),gtbetahat]'*[ones(n,1),gtbetahat]) * [ones(n,1),gtbetahat]'* mean(rt,2);
        gteps = mean(rt,2) - [ones(n,1),gtbetahat]*gtGamma;
        R2g  = 1 - sum(gteps.^2)/sum((mean(rt,2) - mean(mean(rt,2))).^2); 
        
        
        %  R2(h,iMC)  = (temp-mean(temp))'* (temp-mean(temp))/((mean(ret,2)-mean(mean(ret)))'*(mean(ret,2)-mean(mean(ret))));
    
        % step 6: Newy-West Estimation of Avar
        Pi11hat  =  zeros(d*p,d*p);
        Pi12hat  =  zeros(d*p,p);
        Pi22hat  =  zeros(p,p);
        
        for t = 1:T
                    
            Pi11hat       = Pi11hat  +  vec(what(:,t) * vtbar(:,t)')*vec(what(:,t) * vtbar(:,t)')'/T;
            Pi12hat       = Pi12hat  +  vec(what(:,t) * vtbar(:,t)')*vtbar(:,t)'/T;
            Pi22hat       = Pi22hat  +  vtbar(:,t)     * vtbar(:,t)'/T;
            for s = 1:min(t-1,q) 
                
                Pi11hat   = Pi11hat  + 1/T*(1-s/(q+1))* (vec(what(:,t) * vtbar(:,t)')*vec(what(:,t-s) * vtbar(:,t-s)')'+vec(what(:,t-s) * vtbar(:,t-s)')*vec(what(:,t) * vtbar(:,t)')');
                Pi12hat   = Pi12hat  + 1/T*(1-s/(q+1))* (vec(what(:,t) * vtbar(:,t)')*vtbar(:,t-s)' + vec(what(:,t-s) * vtbar(:,t-s)')*vtbar(:,t)');
                Pi22hat   = Pi22hat  + 1/T*(1-s/(q+1))* (vtbar(:,t)     * vtbar(:,t-s)' + vtbar(:,t-s) * vtbar(:,t)' );
            
            end;        
        end;
 
        avarhat   = diag(kron(Gammatilde(2:end)',eye(d))*Pi11hat*kron(Gammatilde(2:end),eye(d))/T + ...
             kron(Gammatilde(2:end)',eye(d))*Pi12hat*etahat'/T + (kron(Gammatilde(2:end)',eye(d))*Pi12hat*etahat')'/T + ...
             etahat*Pi22hat*etahat'/T);
        
        
   
%    res.eig      = tempeig;
%    res.vhat     = vhat;
    res.betahat     = betahat;
  %  res.gthat    = gthat;
    res.Gammahat = Gammahat;
    res.Gammahat_nozero = Gammahat_nozero;
    res.avarhat  = avarhat;
    res.R2 = R2;
    res.R2g = R2g;
    res.R2ghat = R2ghat;
    