function [s1]    = fGCOV(mY,d,sSim)%size mY; d=[Phi_1, Phi_2, ... , Phi_p], matrici riempite per colonna
    H=sSim.H;
    K=sSim.K;
    cN=sSim.cN;
    sample_size=sSim.sample_size;
    iLags=sSim.iLags;
    te = sample_size-iLags;
    eps = zeros(te,cN*K);
    ps=zeros(sample_size,cN);
    mPhi=reshape(d,cN,cN,iLags);
    for t = (1+iLags):sample_size
        S=zeros(cN,1);
        for v=1:iLags
            S=S+mPhi(:,:,v)* mY(t-v, :)';
        end
        ps(t, :) = mY(t, :)' - S;
    end
    ps=ps((iLags+1):end,:);
    eps(:,1:cN) = ps;%
    for i= 1:cN
        for j =2:K  
            eps(:,(j-1)*cN+i)   = (ps(:,i)).^j;%
        end
    end
    meann = mean(eps);
    for s=1:cN*K
        eps(:,s)   = eps(:,s)   - meann(s);    
    end
    corr       = 0;
    vc         = cov(eps);
    invv        =inv(vc);%
    epl        = eps;
    i          = 1;
    while i<(H+1)
        eps    = eps(2:end,:);
        epl    = epl(1:(end-1),:);
        co     = (eps'*epl)/(te-1);
        sco    = co*invv*co'*invv;
        corr   = corr+(sum(diag(sco)))';
        i      = i+1;
    end
    s1 = corr;%
end

