function [moments,wtmat] = Wtmat(data,optweigtmat_selector,maxq)
% function [momemts,wtmat] = Wtmat(data,optweigtmat_selector,maxq);
% Compute the weighting matrix for SMM/GMM estimation
%
% Inputs:
%   - data                      [matrix]    data matrix
%   - optweigtmat_selector      [scalar]    indicator specifying weighting matrix: 0: Identity matrix 1: diagonal 2:
%   - maxq                      [scalar]    bandwidth selector
%
% Outputs:
%   - moments                   [vector]    vector of moments (includes mean of labor, variances and first and second autocovariances)    
%   - wtmat
% 
% adapted by Johannes Pfeifer and Benjamin Born from code by Francisco Ruge-Murcia 



% Moments
nob = size(data,1); % number of observations 
nvar = size(data,2); % number of variables

vars =data(3:nob,end); %get labor
vars = [vars data(3:nob,:).*data(3:nob,:)]; %add variances
for ii = 2:nvar
     vars = [vars data(3:nob,1).*data(3:nob,ii)];  
end
 
autoc1 = data(3:nob,:).*data(2:nob-1,:);
autoc2 = data(3:nob,:).*data(1:nob-2,:);

moments = [mean(vars)';mean(autoc1)';mean(autoc2)'];

% Weighting matrix
if optweigtmat_selector > 0 % weigthing matrix is the optimal or its diagonal %
    nob = size(vars,1)+1;
	wm = (vars-(ones(size(vars,1),1)*mean(vars)))';
	wm = [wm;(autoc1-(ones(size(autoc1,1),1)*mean(autoc1)))'];
    wm = [wm;(autoc2-(ones(size(autoc2,1),1)*mean(autoc2)))'];
    wtmat = wm*wm';
	wtmat = wtmat/(nob-1);
	for iim = 1:maxq
		kapg = wm(:,iim+1:nob-1)*wm(:,1:nob-1-iim)';
		kapg = kapg + kapg';
		kapg = kapg/(nob-1);
		wtmat = wtmat + (1-iim/(maxq+1))*kapg;
    end
	if optweigtmat_selector == 1 % weigthing matrix is the diagonal of the optimal one %
		wtmat = diag(wtmat);
        wtmat = eye(size(wtmat,1)).*(wtmat*ones(1,size(wtmat,1)));
    end
else % weigthing matrix is the identity matrix %
	wtmat = eye(size(moments,1));
end

