%% parforloop functions for Bias simulation

function [metricsReg, metricsOLS, estimators, tpfpReg, tpfpOLS, Areg] = parforloopBias(Xa,Za,RegA1,RegSigma1,F,deltas,real,target,H)

[Na, ~, ~] = size(Xa);

% Estimate and regularize coefficients 
[Est.Ahat, Est.Areg, Est.Sighat, Est.SigReg] = RegularizeAll(Xa,Za,RegA1,RegSigma1,F,deltas);

% calculates the forecast error variance decompositions 
D = FEVD(Est,real,H,target);

% calculate the norms to the true FEVD
norms = evaluate(D);

normPlotReg = squeeze(norms.frobReg);
normPlotRegStruct = squeeze(norms.frobRegStruct); 
meanRegStruct = squeeze(norms.meanRegStruct);
frobRegone = squeeze(norms.frobRegone);


[~, I] = min(normPlotReg); % Target variable for selection of "oracle". Note that later it will be compared to the structural FEVD.
if length(I)>1
    I = min(I);
end


%% Display if the best hyperparameter is a border solution
if size(D.Reg,4) > 1
    lengthAlpha = size(D.Reg,4)/8;
    if rem(I,lengthAlpha) == 0
        disp(' Border solution I = max(lambda)')
    elseif rem(I,lengthAlpha) == 1
        disp(' Border solution I = min(lambda)')
    end
    
elseif size(D.Reg,5) > 1 && strcmp(RegSigma1, 'Threshold')
    lengthEta = size(D.Reg,5)/length(deltas);
    if rem(I,lengthEta) == 0
        disp(' Border solution I = max(eta)')
    elseif rem(I,lengthEta) == 1
        disp(' Border solution I = min(eta)')
    elseif I <= lengthEta
        disp(' Border solution I = min(delta)')
    elseif I > size(D.Reg,5)-lengthEta
        disp(' Border solution I = max(delta)')
    end
    
elseif strcmp(RegSigma1, 'GLASSO')
    if rem(I,50) == 0
        disp(' Border solution I = max(GLASSOgrid)')
    elseif rem(I,50) == 1
        disp(' Border solution I = min(GLASSOgrid)')
    end
    
end

[~,I2] = min(abs(meanRegStruct));
%% Calculate the metrics
metricsReg = [meanRegStruct(I) normPlotRegStruct(I) frobRegone(I) meanRegStruct(I2)];
metricsOLS = [norms.meanHstruct  norms.frobHstruct norms.frobHone];

Areg = eye(Na);
if length(size(D.Reg)) == 4
    Areg = Est.Areg(:,:,:,I);
    D.Reg = D.Reg(:,:,:,I);
elseif length(size(D.Reg)) == 5
    D.Reg = D.Reg(:,:,:,:,I);
end

[tpfpReg, tpfpOLS] = tpfpcurves(D);

estimators(:,1) = D.Reg(:);
estimators(:,2) = D.estimated(:);

end