function res = MP(rt,gt,pmax,q)


% Just computes the MP

%% 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

%% OUTPUT

% gtMP             is d by 1 vector of the risk premia using MP
% avarMP           is d by 1 vector of avar of the risk premia using MP

%% INITIALIZATION

    T  =  size(rt,2);
    n  =  size(rt,1);
    d  =  size(gt,1);
    
    alphahat              =      nan(n,pmax); 
    R2F                   =      nan(1,pmax); %
    R2G                   =      nan(d,pmax); %
    Gammahat              =      nan(d+1,pmax); 
    Gammahat_nozero       =      nan(d+1,pmax); 
    avarhat               =      nan(d+1,pmax);
    avarhat_nozero        =      nan(d,pmax);
    gthat                 =      nan(d,T,pmax); 
    chi2test              =      nan(d,pmax); % 
    SR2hat                =      nan(1,pmax);

    %% ESTIMATION

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

    erhat    =  gtbar * rtbar'*inv(rtbar*rtbar');
    zrbar    =  gtbar - erhat*rtbar;
    varr     =  T^(-1)*rtbar*rtbar';
    avarMP   =  zeros(d,1);
    for ll=1:d
        
        vareta   =  inv(varr)*(rtbar*diag(zrbar(ll,:)')*diag(zrbar(ll,:)')* rtbar')*inv(varr)/T;
        avarMP(ll)=  T^(-1)*(rbar'*vareta*rbar+erhat(ll,:)*varr*erhat(ll,:)');
    end
   
   
%% OUTPUT

   
    res.gtMP           =  mean(erhat*rt,2);
    res.avarMP         = avarMP; 
  