function [Ahat, Areg, Sighat, SigReg]=RegularizeAll(X,Z,RegA,RegSigma,Fa,deltaSig)
[N,T,~] =size(X);
% OLS and samplecov estimation
[Ahat, ~, Sighat]=VARols(X,Z,T);

%% Standardize and Regularize the coefficient matrix first
Xa=(X-mean(X,2))./std(X')';
Za=(Z-mean(Z,2))./std(Z')';
AhatStd=VARols(Xa,Za,T);
switch RegA
    
    case 'OLS'
        Areg(:,:,1,1) = Ahat;
        
    case 'Areg'
        Areg = Fa;
        
    case 'elastic-net'
        
        weights=(abs(AhatStd(:))).^(-1);
        alphas2 = [0 logspace(-5,-.3,5) 0.75 1];
        Fsmaller = 50;
        Areg = [];
        
        kk = 0;
        for alpha = alphas2
            kk=kk+1;
            opts = struct('nlambda',Fsmaller,'alpha', alpha,'intr',0,'penalty_factor',weights,'standardize',0);
            AenetAss= glmnet(sparse(kron(Za',speye(N))), Xa(:),[],opts);
            if Fsmaller-size(AenetAss.beta,2)>0
                AenetAss.beta(:,size(AenetAss.beta,2)+1:Fsmaller) = repmat(AenetAss.beta(:,end),1,Fsmaller-size(AenetAss.beta,2));
            end
            for l=1:size(AenetAss.beta,2)
                AregGrid(:,:,1,l) = diag(std(X'))*reshape(AenetAss.beta(:,l),N,N)/diag(std(Z'));
            end
            
            Areg = cat(4,Areg,AregGrid);
            
        end
    case 'elastic-netV2'
        weights=(abs(AhatStd(:))).^(-1);
        alphas2 = [0 logspace(-5,-.3,5) 0.75 1];
        Fsmaller = 50;
        Areg = [];
        lambda = fliplr(logspace(-3,2,Fsmaller));
        
        kk = 0;
        for alpha = alphas2
            kk=kk+1;
            AregGrid = adaptiveElasticNet(X,Z,lambda,alpha,weights);
            Areg = cat(3,Areg,AregGrid);
        end
        Areg = permute(Areg,[1 2 4 3]);
        
    case 'LASSO'
        a = 1;
        weights=(abs(AhatStd(:))).^(-1);
        opts = struct('nlambda',Fa,'alpha', a,'intr',0,'penalty_factor',weights);
        AenetAss= glmnet(sparse(kron(Z',speye(N))), X(:),[],opts);
        if Fa-size(AenetAss.beta,2)>0
            AenetAss.beta(:,size(AenetAss.beta,2)+1:Fa)=repmat(AenetAss.beta(:,end),1,Fa-size(AenetAss.beta,2));
        end
        for l=1:size(AenetAss.beta,2)
            Areg(:,:,1,l)=reshape(AenetAss.beta(:,l),N,N);
        end
        
    case 'RIDGE'
        opts = struct('nlambda',Fa,'alpha',0,'intr',0);
        AenetAss = glmnet(sparse(kron(Z',speye(N))), X(:),[],opts);
        if Fa-size(AenetAss.beta,2)>0
            AenetAss.beta(:,size(AenetAss.beta,2)+1:Fa)=repmat(AenetAss.beta(:,end),1,Fa-size(AenetAss.beta,2));
        end
        for l=1:size(AenetAss.beta,2)
            Areg(:,:,1,l)=reshape(AenetAss.beta(:,l),N,N);
        end
        
end
%% Regularize the covariance matrix 
switch RegSigma
    
    case 'Ledoit-Wolf'
        for l=1:size(Areg,4)
            res=X-squeeze(Areg(:,:,1,l))*Z;
            res=res-mean(res,2);
            SigenetAssisted=(1/(size(res,2)-1))*(res*res');
            SigReg(:,:,1,l,1)=shrinkLedoit(SigenetAssisted,res);
        end
        
    case 'Threshold'
        etas= [logspace(-4,5,50)];
        for l=1:size(Areg,4)
            res=X-squeeze(Areg(:,:,1,l))*Z;
            res=res-mean(res,2);
            SigenetAssisted=(1/(size(res,2)))*(res*res');
            f=0;
            for jj = 1:length(deltaSig)
                for kk = 1:length(etas)
                    f= f+1;
                    SigReg(:,:,1,l,f) = adaptiveThresh(res,SigenetAssisted,etas(kk),deltaSig(jj));
                end
            end
        end
        
    case 'Sample-Cov'
        for l = 1:size(Areg,4)
            res = X-squeeze(Areg(:,:,1,l))*Z;
            res = res-mean(res,2);
            SigReg(:,:,1,l,1) = (1/(size(res,2)-1))*(res*res');
        end
        
    case 'non-normalized GLASSO'
        for l=1:size(Areg,4)
            res = X-squeeze(Areg(:,:,1,l))*Z;
            res = res-mean(res,2);
            deltaSig = logspace(-6,1,50);
                        
            for f = 1:length(deltaSig)

                SigReg(:,:,1,l,f) = GraphicalLasso(res',deltaSig(f));
            end
        end
        
    case 'GLASSO'
        for l=1:size(Areg,4)
            % calculate residuals
            res = X-squeeze(Areg(:,:,1,l))*Z;
            res = res-mean(res,2);
            deltaSig = [0  logspace(-6,1,99)];
            % standardize residuals
            SaveStd = std(res');
            res = SaveStd'.\res;
                
            w = cov(res');
            theta =inv(w);
            warmInit = 1- (T < 1.1*N); % no warm initialization when dimension is too low. Otherwise the simulation can get stuck.

            for f = 1:length(deltaSig)
                [w, theta, iter(f)] = glasso(N, cov(res'), 0, deltaSig(f).*ones(N), 0, warmInit, 0, 0, 1e-4, 1e4, w, theta);
%                 [w, theta, iter] = glasso(numVars, cov(pop), 0, lambda, approximate, warmInit, verbose, penalDiag, tolThreshold, maxIter, w, theta);
                SigReg(:,:,1,1,f) =  diag(SaveStd)*w*diag(SaveStd);
            end

        end
        
    case 'adaptive GLASSO'
        % Glasso with weights corresponding to the sample-cov
        for l=1:size(Areg,4)
            res = X-squeeze(Areg(:,:,1,l))*Z;
            res = res-mean(res,2);
            deltaSig = logspace(-6,1,50);
            SigenetAssisted=(1/(size(res,2)))*(res*res');
            
            for f = 1:length(deltaSig)
                ThetaHat = (abs(inv(SigenetAssisted))).^(-1);
                deltaSigMat = deltaSig(f)*ThetaHat/mean(mean(ThetaHat));
                SigReg(:,:,1,l,f) = GraphicalLasso(res',deltaSigMat);
            end
        end
        
end

end
