function MainGibbsFunc(CurrFileName,GFEVcHorizon,TotalNumGibbsDraws,RunStartDraw,RunEndDraw)

%Initalize
NetworkName = strcat(CurrFileName,'_H',num2str(GFEVcHorizon),'_Draws',num2str(TotalNumGibbsDraws));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Setup parameters and priors if current run begins with the first iteration.  Otherwise, load variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if RunStartDraw == 1
    %Setup model parameters
    NumFactorsK = 2;
    PositiveFactorSignMember = "RA_1";
    RawData = readtable(strcat(CurrFileName,'.csv'));
    SaveFullNetwork = true;
    warning('off'); TickerMetadata = readtable(strcat(CurrFileName,'_Metadata.csv')); warning('on');
    NumSeriesN = size(RawData,2)-1;
    NumT =  size(RawData,1);
    KappaLambda = 0.2;
    KappaSigma = 0.01;
    KappaGamma = 0.2;
    KappaRho = 0.2;
    EpsilonSmall = 1e-8;
    TauDegFree = 3;

    TimerStartValue = tic;

    rng(0);

    %Clear old diary
    if exist(strcat('.\Output\',NetworkName,'_log.log'),'file')==2
        delete(strcat('.\Output\',NetworkName,'_log.log'));
    end

    %Covariance matrix based PCA inclusive of the means
        %Create actual returns matrix
        RtMatrix = table2array(RawData(:,2:end));
        
        %Get PCA factors and loadings
        [VV,DD] = eig(cov(RtMatrix));
        [~,indx] = sort(diag(DD),'descend');
        VVsort = VV(:,indx);
        EigenVect = VVsort(:,1:NumFactorsK); %NxK matrix of the first K eigenvectors of the covariance matrix

        PCA_Factors_Raw = RtMatrix * EigenVect;
        PCA_Loadings_Raw = EigenVect;
        
        %Divide the factors by their standard deviations and multiply their loadings by them
        FactStdDev = sqrt(diag(cov(PCA_Factors_Raw)));
        PCA_Factors_Adj = PCA_Factors_Raw / diag(FactStdDev);
        PCA_Loadings_Adj = PCA_Loadings_Raw * diag(FactStdDev);
        
        %Get location of panel member that needs to load positively on all factors
        PositiveFactorSignMemberColumn = find(string(RawData.Properties.VariableNames) == PositiveFactorSignMember)-1;
        
        %Adjust factor and loadings signs
        PCA_Factors_Adj  = PCA_Factors_Adj  * diag(sign(PCA_Loadings_Adj(PositiveFactorSignMemberColumn,:)));
        PCA_Loadings_Adj = PCA_Loadings_Adj * diag(sign(PCA_Loadings_Adj(PositiveFactorSignMemberColumn,:)));
        
        %Check adjusted factors & loadings
            %Confirm span same space as raw PCA factors & loadings
            PCA_Recover_Raw = PCA_Factors_Raw * PCA_Loadings_Raw';
            PCA_Recover_Adj = PCA_Factors_Adj * PCA_Loadings_Adj';
            CheckVal = max(abs(PCA_Recover_Raw-PCA_Recover_Adj),[],"all");
            if CheckVal > 1e-8
                error(['PCA recovery differences are large after sign and scale adjustment:  ',num2str(CheckVal,"%10.8f")]);
            end
        
            CheckVal = max(abs(cov(PCA_Factors_Adj) - eye(NumFactorsK)),[],"all");
            if CheckVal > 1e-8
                error(['Adjusted factors are not standardaized to have a variance of one & uncorrelated:  ',num2str(CheckVal,"%10.8f")]);
            end
        
            CheckVal = min(sign(PCA_Loadings_Adj(PositiveFactorSignMemberColumn,:)));
            if CheckVal ~= 1
                error(['Adjusted factors do not all have positive signs on target panel member:  ',num2str(CheckVal)]);
            end
        
            CheckVal = max(abs( arrayfun(@(kk) corr(PCA_Factors_Raw(:,kk), PCA_Factors_Adj(:,kk)), 1:NumFactorsK) * diag(sign(PCA_Loadings_Raw(PositiveFactorSignMemberColumn,:)))-1 ));
            if CheckVal > 1e-8
                error(['Correlations between raw and new factors when sign adjusted are not one:  ',num2str(CheckVal,"%10.8f")]);
            end
        
            clear VV DD indx VVsort EigenVect FactStdDev PCA_Recover_Adj PCA_Recover_Raw CheckVal PCA_Loadings_Raw;

    %Get industries for the firms
        %Left join metadata to sample name list so only get data for firms in the node list.  Remove lnRetD name prefix if exists
        NodeNameList = array2table(strrep(RawData.Properties.VariableNames(2:end),"lnRetD","")', 'VariableNames',{'CleanTicker'});
        NodeIndustryList = join(NodeNameList,TickerMetadata,'Keys',{'CleanTicker','CleanTicker'});
        NodeIndustryList = NodeIndustryList.BEA_SECTOR;
        clear NodeNameList TickerMetadata;

    %Setup EstimatedNetwork object to be saved
        %Save factors, RAtoRA & RItoRA to disk
        %Also, save initial raw PCA & adjusted factors for reference if need to debug
        %Save inside loop in case of failue and need to debug
        %NetworkName
        %IndustriesList     -   Industry list to match networks     NumIndustries,1
        %FDraws             -   Common factors                      NumT,NumFactorsK,TotalNumGibbsDraws
        %RAtoRA             -   3D array for industry network       NumIndustries,NumIndustries,TotalNumGibbsDraws
        %RItoRA             -   3D array for industry network       NumIndustries,NumIndustries,TotalNumGibbsDraws
        %PCA_Factors_Raw    -   PCA common factors                  NumT,NumFactorsK
        %PCA_Factors_Adj    -   Initial adjusted common factors     NumT,NumFactorsK
            EstimatedNetwork.NetworkName        = NetworkName;
            EstimatedNetwork.IndustriesList     = string(unique(NodeIndustryList));
            EstimatedNetwork.FDraws             = zeros(NumT,NumFactorsK,TotalNumGibbsDraws);        
            NumIndustries                       = length(EstimatedNetwork.IndustriesList);
            EstimatedNetwork.RAtoRA             = zeros(NumIndustries,NumIndustries,TotalNumGibbsDraws);
            EstimatedNetwork.RItoRA             = zeros(NumIndustries,NumIndustries,TotalNumGibbsDraws);
            EstimatedNetwork.PCA_Factors_Raw    = PCA_Factors_Raw;
            EstimatedNetwork.PCA_Factors_Adj    = PCA_Factors_Adj;

            %Objects to save firm level network if SaveFullNetwork is true
            if SaveFullNetwork
                EstimatedNetwork.RAtoRA_Firm = zeros(NumSeriesN,NumSeriesN,TotalNumGibbsDraws);
                EstimatedNetwork.RIFtoRAF    = zeros(NumSeriesN + NumFactorsK,NumSeriesN + NumFactorsK,TotalNumGibbsDraws);
            end

        clear NumIndustries PCA_Factors_Raw;

    %Setup constants for priors and dummy observations
        %Priors for the common factors and the LBar matrix
        GibbsDrawPrior.F(:,:) = PCA_Factors_Adj;
        LBar = PCA_Loadings_Adj;
        clearvars PCA_Factors_Adj PCA_Loadings_Adj;
        
        %Initial values for the RI series
        GibbsDrawPrior.RIdio(:,:) = RtMatrix - GibbsDrawPrior.F(:,:) * LBar';
        
        %SigmaSqHat from the variances of the RI series
        SigmaSqHat = diag(cov(GibbsDrawPrior.RIdio(:,:)));
        
        %Set priors for periods zero and 1, assuming they are the same
            %Set S1 to a (K + N) × 1 zero vector.
            S1 = zeros(NumFactorsK + NumSeriesN,1);

            %Use initial PCA first period values for the common factors' S1 means
            for ii = 1:1:NumFactorsK
                S1(ii,1) = GibbsDrawPrior.F(1,ii);
            end

            %Set initial S0 draw to S1
            GibbsDrawPrior.S0Draw(:,1) = S1;

            %Set P1 to a diagonal matrix with ones for the first K terms relating to the common factors and the σ^2 hat variances for the idiosyncratic returns' entries.
            P1 = diag([ones(NumFactorsK,1); SigmaSqHat]);

        %The prior covariance matrix for vec(Γ 0 ) is G ≡ κ 2 Γ I K 2
        GG = KappaGamma^2 * eye(NumFactorsK^2);
        
        %Create the dummy observations
        YD = [zeros(NumSeriesN,NumSeriesN);
              zeros(1,NumSeriesN);
              diag(sqrt(SigmaSqHat))];
        XD = [zeros(NumSeriesN,1) diag(sqrt(SigmaSqHat)) / KappaRho;
              EpsilonSmall zeros(1,NumSeriesN);
              zeros(NumSeriesN,NumSeriesN+1)];
        
        %Initialize draw for Tau at Sii
        Sii = KappaSigma*SigmaSqHat;
        GibbsDrawPrior.Tau(:,1) = Sii;
else
    load(strcat('.\Output\',NetworkName,'_IntermediateVariables.mat')); %#ok<LOAD>
    rng(RndState);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Setup logging & number of cores
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Start diary
diary(strcat('.\Output\',NetworkName,'_log.log'));

%Customized code to set number of threads.  On Windows platforms with multiple processor groups, the default maximum number of computational threads is equal to the number of physical cores within the processor group that is running MATLAB.
[~, CompName] = system('hostname');
if strtrim(CompName) == "EGi9d"
    maxNumCompThreads(24);
elseif strtrim(CompName) == "EGi7"
    maxNumCompThreads(6);
else
    maxNumCompThreads('automatic');
end
disp(strcat("Draw start:  ",num2str(RunStartDraw),".  Draw end:  ",num2str(RunEndDraw)));
disp(strcat("Computer name:  ",CompName));
disp(strcat("Number of computational threads:  ",num2str(maxNumCompThreads)));
fprintf('\n');
clearvars CompName;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Run the Gibbs sampling algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Setup iterations so saves data and fails in a way that could be restarted from prior draw
for gDraw = RunStartDraw:1:RunEndDraw
    %Draw status update
    if gDraw==1
        CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
        disp(strcat("Beginning MCMC iterations.  Runtime:  ",CurrTime));
        fprintf('\n');
        clear CurrTime;
    elseif gDraw==100*floor(.01*gDraw)
        CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
        disp(strcat("Now running MCMC iteration ",num2str(gDraw)," (Out of ",num2str(TotalNumGibbsDraws),".  Runtime:  ",CurrTime,")"));
        fprintf('Memory used by Matlab:  %.0f MB\n', memory().MemUsedMATLAB/(1024 * 1024));
        aaa = whos;
        aaa = struct2table(aaa);
        aaa.MB = aaa.bytes/(1024 * 1024);
        [~,idx] = sort(aaa.bytes,'descend');
        disp(aaa(idx,["name" "size" "MB" "bytes" "class"]));
        clearvars aaa idx;
        fprintf('\n');
        clear CurrTime;
        diary off;
        diary(strcat('.\Output\',NetworkName,'_log.log'));
    end


    %Draw Lambda and Tau

    %Get prior iteration information that will reuse several times.
    FPrev        = GibbsDrawPrior.F(:,:);
    Ftm1Prev     = [GibbsDrawPrior.S0Draw(1:NumFactorsK,:)'; FPrev(1:end-1,:)];
    RIdioPrev    = GibbsDrawPrior.RIdio(:,:);
    RIdiotm1Prev = [GibbsDrawPrior.S0Draw(NumFactorsK+1:end,:)'; RIdioPrev(1:end-1,:)];
    TauPrev      = GibbsDrawPrior.Tau(:,1);

    %While loop checks factors are ordered in descending order of the cumulative variance their innovations explain for the Rt equity returns.
    DrawLambdaTau = 1;
    while DrawLambdaTau > 0
        %Initialize Lambda and Tau objects to fill out
        CurrLambda = zeros(NumSeriesN,NumFactorsK);
        CurrTau    = zeros(NumSeriesN,1);

        %Run a loop over the NumSeriesN observation equations.
        for ii = 1:1:NumSeriesN
        %%%parfor ii = 1:1:NumSeriesN
            %Subtract the idiosyncratic returns from the observed ones for firm ii.
            YCii = RtMatrix(:,ii) - RIdioPrev(:,ii);

            %Draw the Λ coefficients for observation equation ii
            LBariiPrime = LBar(ii,:)';
            LVarianceii = (KappaLambda^2 * SigmaSqHat(ii)) * eye(NumFactorsK);
            CoeffVariance = inv(inv(LVarianceii) + (FPrev' * FPrev)/TauPrev(ii));
            CoeffMean     = CoeffVariance * (LVarianceii\LBariiPrime + (FPrev' * YCii)/TauPrev(ii)); %#ok<MINV>
                %Draw the coefficients - Long way is faster than the built in function
                %%%CurrLambda(ii,:) = mvnrnd(CoeffMean,CoeffVariance);
                CurrLambda(ii,:) = CoeffMean+(randn(1,size(CoeffMean,1))*chol(CoeffVariance))';

            %Draw Tauii
            Residii = YCii - FPrev * CurrLambda(ii,:)';
            TauScaleii = (Sii(ii) + Residii' * Residii)/2;
            %Make a draw from an inverse gamma distribution with given degrees of freedom and scale
            %   See BoE textbook p. 11
                NormDrawVect = randn((TauDegFree+NumT),1)/sqrt(2*TauScaleii);
                CurrTau(ii,1) = 1/(NormDrawVect'*NormDrawVect);
        end
        clear ii YCii LBariiPrime LVarianceii CoeffVariance CoeffMean Residii TauScaleii NormDrawVect;

        %Check factors are ordered in descending order of the cumulative variance their innovations explain for the Rt equity returns.
        if issorted(diag(CurrLambda'*CurrLambda),'descend')
            GibbsDrawNew.Lambda(:,:) = CurrLambda;
            GibbsDrawNew.Tau(:,:) = CurrTau;
            DrawLambdaTau = -1;
        elseif DrawLambdaTau > 50000
            %Error message if redraw too many times.
            error(['Too many Lambda and Tau redraws for MCMC iteration ',num2str(gDraw),' (Draw number ',num2str(DrawLambdaTau),')']);
        else
            %Print a message if redraw.
            CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
            disp(strcat("Redrawing Lambda and Tau due to order for MCMC iteration ",num2str(gDraw)," (Draw number ",num2str(DrawLambdaTau),".  Runtime:  ",CurrTime,")"));
            fprintf('\n');
            clear CurrTime;
            DrawLambdaTau = DrawLambdaTau+1;
        end
    end
    clear TauPrev CurrLambda CurrTau DrawLambdaTau;


    %Draw Gamma

    %Follows BoE textbook formulas on p. 34
    KronCalc      = kron(eye(NumFactorsK), Ftm1Prev' * Ftm1Prev);
    CoeffVariance = inv(inv(GG) + KronCalc);
    vecGammaOLS   = reshape( (Ftm1Prev' * Ftm1Prev) \ (Ftm1Prev' * FPrev) ,NumFactorsK^2,1);
    CoeffMean     = CoeffVariance * (KronCalc * vecGammaOLS); %#ok<MINV>

    %While loop checks for stability.
    DrawGamma = 1;
    while DrawGamma > 0
        %Draw the coefficients - Long way is faster than the built in function
            %%%vecGammaMatrix = mvnrnd(CoeffMean,CoeffVariance);
            vecGammaMatrix = CoeffMean+(randn(1,size(CoeffMean,1))*chol(CoeffVariance))';
        GammaMatrix = reshape(vecGammaMatrix,NumFactorsK,NumFactorsK)';

        %Check stability.  A VAR(1) process is stable if the largest eigenvalue of its coefficient matrix has modulus < 1
        if max(abs(eig(GammaMatrix))) < 1
            GibbsDrawNew.Gamma(1:NumFactorsK,1:NumFactorsK) = GammaMatrix;
            DrawGamma = -1;
        elseif DrawGamma > 50000
            %Error message if redraw too many times.
            error(['Too many Gamma redraws for MCMC iteration ',num2str(gDraw),' (Draw number ',num2str(DrawGamma),')']);
        else
            %Print a message if redraw.
            CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
            disp(strcat("Redrawing Gamma due to stability for MCMC iteration ",num2str(gDraw)," (Draw number ",num2str(DrawGamma),".  Runtime:  ",CurrTime,")"));
            fprintf('\n');
            clear CurrTime;
            DrawGamma = DrawGamma+1;
        end
    end
    clear KronCalc CoeffVariance vecGammaOLS CoeffMean vecGammaMatrix GammaMatrix DrawGamma;


    %Direct Monte Carlo sampling from posterior of Σ and then ρ

    %Calculate the necessary parameters
    Xg = [ones(NumT,1) RIdiotm1Prev];
    YgStar = [RIdioPrev' YD']';
    XgStar = [Xg'        XD']';
    RhogStar = (XgStar' * XgStar) \ (XgStar' * YgStar);
    SigmaStar = (YgStar - XgStar * RhogStar)' * (YgStar - XgStar * RhogStar);

    %While loop checks for stability.
    DrawSigmaRho = 1;
    while DrawSigmaRho > 0
        %Draw Σ from its posterior distribution
        %Make a draw from an inverse Wishart distribution with given degrees of freedom and scale
        %   See BoE textbook p. 35
            %Returns a matrix R of n random vectors chosen from the same multivariate normal distribution, with mean vector mu and covariance matrix Sigma.
            ZZ = mvnrnd(zeros(1,size(SigmaStar,1)),  inv(SigmaStar),  NumT + NumSeriesN + 2);
            CurrSigma = inv(ZZ' * ZZ);

        %Draw ρ conditional on Σ
        %https://en.wikipedia.org/wiki/Matrix_normal_distribution
        CurrRho = RhogStar + ( chol(inv(XgStar' * XgStar))' * randn(NumSeriesN+1,NumSeriesN) ) * chol(CurrSigma);

        %Check stability.  A VAR(1) process is stable if the largest eigenvalue of its coefficient matrix has modulus < 1, and the transpose of a matrix has the same eigenvalues
        if max(abs(eig(CurrRho(2:end,:)))) < 1
            GibbsDrawNew.Sigma(1:NumSeriesN,1:NumSeriesN) = CurrSigma;
            %Separate Rho0 and Rho
                GibbsDrawNew.Rho0(1:NumSeriesN,1) = CurrRho(1,:)';
                GibbsDrawNew.Rho(1:NumSeriesN,1:NumSeriesN) = CurrRho(2:end,:)';
            DrawSigmaRho = -1;
        elseif DrawSigmaRho > 50000
            %Error message if redraw too many times.
            error(['Too many Sigma and Rho redraws for MCMC iteration ',num2str(gDraw),' (Draw number ',num2str(DrawSigmaRho),')']);
        else
            %Print a message if redraw.
            CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
            disp(strcat("Redrawing Sigma and Rho due to stability for MCMC iteration ",num2str(gDraw)," (Draw number ",num2str(DrawSigmaRho),".  Runtime:  ",CurrTime,")"));
            fprintf('\n');
            clear CurrTime;
            DrawSigmaRho = DrawSigmaRho+1;
        end
    end
    clear Xg YgStar XgStar RhogStar SigmaStar CurrSigma CurrRho DrawSigmaRho ZZ;


    %Draw the F and RIdio states
    %Use implementation of the Durbin and Koopman (2002) simulation smoother provided in Jarocinski (2015) to draw the states
    % The model is:
    % y(t) = dd + ZZ a(t) + WW eps(t),  eps(t) ~ N(0,I)
    % a(t+1) = cc + TT a(t) + RR eta(t), eta(t) ~ N(0,I)
    % a(1) ~ N(a1,C1*C1')

    DK_dd = zeros(NumSeriesN,1);
    DK_ZZ = [GibbsDrawNew.Lambda(1:NumSeriesN,1:NumFactorsK) eye(NumSeriesN)];
    DK_WW = chol(diag(GibbsDrawNew.Tau(1:NumSeriesN,1)))';
    DK_cc = [zeros(NumFactorsK,1); GibbsDrawNew.Rho0(1:NumSeriesN,1)];
    DK_TT = [GibbsDrawNew.Gamma(1:NumFactorsK,1:NumFactorsK)    zeros(NumFactorsK,NumSeriesN);
             zeros(NumSeriesN,NumFactorsK)                      GibbsDrawNew.Rho(1:NumSeriesN,1:NumSeriesN)];
    DK_RR = chol([eye(NumFactorsK)                              zeros(NumFactorsK,NumSeriesN);
                  zeros(NumSeriesN,NumFactorsK)                 GibbsDrawNew.Sigma(1:NumSeriesN,1:NumSeriesN)])';
    DK_C1 = chol(P1)';

    %Draw the states conditional on y, DKSim is nstates x T
        %Optimized version of Jarocinski (2015), A note on implementing the Durbin and Koopman simulation smoother.
            yplus = nan(NumSeriesN,NumT);
            aplus = nan(NumSeriesN+NumFactorsK,NumT+1);
            aplus(:,1) = DK_C1*randn(NumSeriesN+NumFactorsK,1); % draw the first state with a1=0
            for tt = 1:NumT
                yplus(:,tt) = DK_ZZ*aplus(:,tt) + DK_WW*randn(NumSeriesN,1);
                aplus(:,tt+1) = DK_TT*aplus(:,tt) + DK_RR*randn(NumSeriesN+NumFactorsK,1);
            end
            aplus(:,end) = [];
            yy = RtMatrix' - yplus;
	        clearvars yplus;

        % allocate space
        ZZp_FFFinv_vvv = cell(1,NumT);
        KKK_DK_ZZ = cell(1,NumT);

        % Kalman filter on yy
        % compute frequently used matrices
        HH = DK_WW*DK_WW';
        RQR = DK_RR*DK_RR';
        % initialize Kalman filter
        at = S1; % at|I(t-1)
        Pt = DK_C1*DK_C1'; % Pt|I(t-1)
        for tt = 1:NumT
            vt = yy(:,tt) - DK_dd - DK_ZZ*at;
            Ftinv = inv(DK_ZZ*Pt*DK_ZZ' + HH);%\eye(NumSeriesN); %Update
            TTPt = DK_TT*Pt;                           Kt = TTPt*DK_ZZ'*Ftinv; %#ok<MINV>
            % update at,Pt; from now on their interpretation is "a(t+1),P(t+1)"
            at = DK_cc + DK_TT*at + Kt*vt;
            Kt_DK_ZZ = Kt*DK_ZZ;
            Pt = TTPt*(DK_TT - Kt_DK_ZZ)' + RQR;
            % store the quantities needed later for smoothing
            ZZp_FFFinv_vvv{tt} = DK_ZZ'*(Ftinv*vt); %#ok<MINV>
            KKK_DK_ZZ{tt} = Kt_DK_ZZ; %KKK_DK_ZZ causes the RAM problem between these two for loops
        end
        clearvars yy HH at Pt vt Ftinv TTPt Kt Kt_DK_ZZ;

        % Kalman smoother
        % backwards recursion on r
        rrr = cell(1,NumT); rrr_tp1 = zeros(NumSeriesN+NumFactorsK,1);
        for tt = NumT-1:-1:1
            rrrCurr = ZZp_FFFinv_vvv{tt+1} + (rrr_tp1'*DK_TT - rrr_tp1'*KKK_DK_ZZ{tt+1})'; KKK_DK_ZZ(end) = []; rrr_tp1 = rrrCurr; rrr{tt} = rrrCurr;
        end
        % one more iteration to get r0
        r0 = ZZp_FFFinv_vvv{1} + (rrr{1}'*DK_TT - rrr{1}'*KKK_DK_ZZ{1})';
        clearvars rrrCurr ZZp_FFFinv_vvv KKK_DK_ZZ rrr_tp1;
        % allocate space for smoothed states
        DKSim = nan(NumSeriesN+NumFactorsK,NumT);
        % forwards recursion to compute smoothed states from r - DK(2002),eq.8
        DKSim(:,1) = S1 + DK_C1*(DK_C1'*r0); % initialize the forward recursion
        for tt = 2:NumT
            DKSim(:,tt) = DK_cc + DK_TT*DKSim(:,tt-1) + RQR*(rrr{tt-1});
        end
        DKSim = DKSim + aplus;
        clearvars tt aplus RQR rrr r0;

    GibbsDrawNew.F(1:NumT,1:NumFactorsK) = DKSim(1:NumFactorsK,:)';
        EstimatedNetwork.FDraws(:,:,gDraw) = GibbsDrawNew.F;
    GibbsDrawNew.RIdio(1:NumT,1:NumSeriesN) = DKSim(NumFactorsK+1:end,:)';

    %Draw for S0 - prior mean and variance same as for time period one
    GainK = (P1 * DK_TT') / (DK_TT*P1*DK_TT'+DK_RR * DK_RR');
    Meant = S1 + GainK * (DKSim(:,1) - DK_cc - DK_TT*S1);
    Vart  = P1 - GainK * DK_TT * P1;
    GibbsDrawNew.S0Draw(:,1) = Meant + chol(Vart)'*randn(size(P1,1),1);


    %Cleanup variables & set current draw to prior
    GibbsDrawPrior = GibbsDrawNew;
    clear GibbsDrawNew FPrev Ftm1Prev RIdioPrev RIdiotm1Prev DK_dd DK_ZZ DK_WW DK_cc DK_TT DK_RR DK_C1 DKSim GainK Meant Vart;


    %Calculate the GFEVc networks
        %Calculate A1 & ThetaCoVar matrices
        A1 = [zeros(NumSeriesN,NumSeriesN)      GibbsDrawPrior.Rho              GibbsDrawPrior.Lambda * GibbsDrawPrior.Gamma;
              zeros(NumSeriesN,NumSeriesN)      GibbsDrawPrior.Rho              zeros(NumSeriesN,NumFactorsK);
              zeros(NumFactorsK,NumSeriesN)     zeros(NumFactorsK,NumSeriesN)   GibbsDrawPrior.Gamma];
        ThetaCoVar = [(GibbsDrawPrior.Sigma + GibbsDrawPrior.Lambda * GibbsDrawPrior.Lambda')   GibbsDrawPrior.Sigma            GibbsDrawPrior.Lambda;
                      GibbsDrawPrior.Sigma                                                      GibbsDrawPrior.Sigma            zeros(NumSeriesN,NumFactorsK);
                      GibbsDrawPrior.Lambda'                                                    zeros(NumFactorsK,NumSeriesN)   eye(NumFactorsK)];

        %Calculate combined Ai matrices
            NumVars = 2*NumSeriesN+NumFactorsK;
            Ais = zeros(NumVars, NumVars*(GFEVcHorizon+1));
        %Iterate through terms from 0 to GFEVcHorizon
            for hh = 0:1:GFEVcHorizon
                %Calculate Ah matrix
                if hh == 0
                    Ah = eye(NumVars);
                else
                  Ah = zeros(NumVars,NumVars);
                  for ii = 1:1:min(hh,1)
                    Ah = A1(:,(1+(ii-1)*NumVars):(ii*NumVars)) * Ais(:,(1+(hh-ii)*NumVars):((hh-ii+1)*NumVars)) + Ah;
                  end
                end
                %Fill in Ais with Ah
                Ais(:,(1+hh*NumVars):((hh+1)*NumVars)) = Ah;
            end
        %Calculate GFEVc
            GFEVc = zeros(NumVars, NumVars);
            for hh = 0:1:GFEVcHorizon
                Ah = Ais(:,(1+hh*NumVars):((hh+1)*NumVars));
                GFEVc = (Ah * ThetaCoVar).^2 + GFEVc;
            end
            GFEVc = GFEVc * diag(   1./diag(ThetaCoVar)   );
        %H=1 GFEVc check
        %%%GFEVc1 = (  ThetaCoVar.^2 + (A1 * ThetaCoVar).^2  ) * diag(diag(ThetaCoVar))^-1;
        %%%max(abs(   GFEVc - GFEVc1   ),[],'all')
    clear A1 ThetaCoVar NumVars Ais hh Ah ii;


    %Aggegate networks to industry level
        %Get RAtoRA and RItoRA firm level networks
            RAtoRA_Firm = GFEVc(1:NumSeriesN,1:NumSeriesN);
            RItoRA_Firm = GFEVc(1:NumSeriesN,NumSeriesN+1:2*NumSeriesN);

        %RAtoRA_Firm;
            WorkingTable = array2table(RAtoRA_Firm);
            WorkingTable.Ind = NodeIndustryList;
            WorkingTable = groupsummary(WorkingTable,'Ind',"sum");
            WorkingTable = removevars(WorkingTable,["Ind","GroupCount"]);
            WorkingTable = array2table(table2array(WorkingTable)');
            WorkingTable.Ind = NodeIndustryList;
            WorkingTable = groupsummary(WorkingTable,'Ind',"sum");
            WorkingTable = removevars(WorkingTable,["Ind","GroupCount"]);
            WorkingTable = array2table(table2array(WorkingTable)');
        EstimatedNetwork.RAtoRA(:,:,gDraw) = table2array(WorkingTable);

        %RItoRA_Firm;
            WorkingTable = array2table(RItoRA_Firm);
            WorkingTable.Ind = NodeIndustryList;
            WorkingTable = groupsummary(WorkingTable,'Ind',"sum");
            WorkingTable = removevars(WorkingTable,["Ind","GroupCount"]);
            WorkingTable = array2table(table2array(WorkingTable)');
            WorkingTable.Ind = NodeIndustryList;
            WorkingTable = groupsummary(WorkingTable,'Ind',"sum");
            WorkingTable = removevars(WorkingTable,["Ind","GroupCount"]);
            WorkingTable = array2table(table2array(WorkingTable)');
        EstimatedNetwork.RItoRA(:,:,gDraw) = table2array(WorkingTable);

        %Save firm level network if SaveFullNetwork is true.
        if SaveFullNetwork
            EstimatedNetwork.RAtoRA_Firm(:,:,gDraw) = RAtoRA_Firm;
            EstimatedNetwork.RIFtoRAF(:,:,gDraw)    = [   GFEVc(1:NumSeriesN,                             NumSeriesN+1:2*NumSeriesN)   GFEVc(1:NumSeriesN,                             2*NumSeriesN+1:2*NumSeriesN+NumFactorsK)   ;
                                                          GFEVc(2*NumSeriesN+1:2*NumSeriesN+NumFactorsK,  NumSeriesN+1:2*NumSeriesN)   GFEVc(2*NumSeriesN+1:2*NumSeriesN+NumFactorsK,  2*NumSeriesN+1:2*NumSeriesN+NumFactorsK)   ];
        end

    clear RAtoRA_Firm RItoRA_Firm GFEVc WorkingTable;

    %Save results.  Save inside loop in case of failure and need to debug
        save(strcat('.\Output\',NetworkName,'.mat'),"EstimatedNetwork","-v7.3");
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save variables for next run or empty file to NotifyComplete directory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if RunEndDraw == TotalNumGibbsDraws
    CurrTime=string(seconds(toc(TimerStartValue)),'hh:mm:ss.SSSS');
    fid = fopen(strcat('.\Output\NotifyComplete_',NetworkName,'.txt'),'w');
    fprintf(fid,'%s',strcat(NetworkName," done!  Runtime:  ",CurrTime));
    fclose(fid);
    fprintf('\n');
	fprintf('Memory used by Matlab:  %.0f MB\n', memory().MemUsedMATLAB/(1024 * 1024));
    aaa = whos;
    aaa = struct2table(aaa);
    aaa.MB = aaa.bytes/(1024 * 1024);
    [~,idx] = sort(aaa.bytes,'descend');
    disp(aaa(idx,["name" "size" "MB" "bytes" "class"]));
    clearvars aaa idx;
    fprintf('\n');
	disp(strcat("Done!  Runtime:  ",CurrTime));
    clear CurrTime fid;
else
    RndState = rng;
    save(strcat('.\Output\',NetworkName,'_IntermediateVariables.mat'),"EpsilonSmall","EstimatedNetwork","GG","GibbsDrawPrior","KappaGamma","KappaLambda","KappaRho","KappaSigma","LBar","NodeIndustryList","NumFactorsK","NumSeriesN","NumT","P1","PositiveFactorSignMember","PositiveFactorSignMemberColumn","RawData","RtMatrix","S1","SigmaSqHat","Sii","TauDegFree","XD","YD","TimerStartValue","RndState","SaveFullNetwork","-v7.3");
end

diary off;

end