clear all
close all
disp('#### Application for "the Great Moderation through the lense of sectoral spillovers" ####')
curdict=cd;
addpath([curdict '/glmnet_matlab'])
addpath([curdict '/Graphical Lasso'])
addpath([curdict '/Functions'])
addpath([curdict '/InputData'])
if isempty(gcp('nocreate'))==1
    myCluster = parcluster('local');
    delete(myCluster.Jobs);
    parpool([4 28],'IdleTimeout',60);
end
% dbstop if error
pctRunOnAll warning('off','all')
%% Load and cut data
NAIClevel = 3;
disp(['Read in data for NAIC level ' num2str(NAIClevel) ' and fill missing values.'])
file_name = 'IP98updated_formulas.xlsx';
StartColumn = 1;

RawTable = readtable(file_name,'Sheet',['L' num2str(NAIClevel)],'Range',[NAIClevel+3 StartColumn]);
RawTableWeights = readtable(file_name,'Sheet',['L' num2str(NAIClevel) '-W'],'Range',[NAIClevel+3 StartColumn]);
RawTableFactors = readtable(file_name,'Sheet',['AggregateShocks'],'Range',[3 1]);

RawData = table2array(RawTable(:,2:end));
RawData = fillmissing(RawData,'linear',1,'EndValues','nearest');

weights = table2array(RawTableWeights(:,2:end));
weights = fillmissing(weights,'linear',1,'EndValues','nearest')';

FactorData = table2array(RawTableFactors(:,2:end));
FactorData = fillmissing(FactorData,'linear',1,'EndValues',0)';

names = RawTable.Properties.VariableNames(2:end);
if NAIClevel==3 
    transformNames;
end

Inventories = table2array(readtable(['Inventories.xlsx']));
Shipments = table2array(readtable(['Shipments.xlsx']));

% align factors with data
timeVectorIP = datetime(RawTable.x_DESC);
timeVectorFactor = datetime(RawTableFactors.x_DESC);

month_date = month(timeVectorFactor(2:end));
year_date = year(timeVectorFactor(2:end));

DataOffset = find(timeVectorFactor == timeVectorIP(1)) -1; % Set where the dependent variable starts
Data(:,1+DataOffset:size(FactorData,2)) = RawData';

clearvars -except Data FactorData month_date year_date Inventories Shipments names weights NAIClevel DataOffset
%% Prepare IP growth data and factors
Y = 1200*(log(Data(:,2:end))-log(Data(:,1:end-1)));

% Factors, to be log-growth if they are in levels
Commodity = (log(FactorData(1,2:end)) - log(FactorData(1,1:end-1)));    
TFP = FactorData(2,2:end);
MPx = FactorData(3,2:end);
MPx2 = FactorData(4,2:end);

%% Estimation
% Specifications
RegA = 'elastic-net';                   % regularization for AR coefficient
RegSigma = 'allCV_AIC';                 % regularization for Covariance matrix. Put "allCV_AIC" if it should Cross validate the hyperparameters and select the regularization method with AIC
BreakTests = 'all-dynamic';             % selection of dynamic. Either VAR-dynamic, Sigma-dynamic, all-static or all-dynamic

[N, T] = size(Y); deltas=0:0.0025:3;
p = 1;                              % Set number of VAR lags
RW = 144;                           % Rolling-window Bandwidth
ts = DataOffset + [145:RW:589];     % window end-dates
FactNumLags = 12;                   % number of Laggs for factors
Kfold = 12;                         % folds for cross validation

firstTic = tic;
disp(['  Onset of the Great Moderation has been set to ' num2str(month_date(ts(1)+1)) '/' num2str(year_date(ts(1)+1))])

switch BreakTests
    case 'VAR-dynamic'
        % make dynamic VAR matrix, but static Covariance matrix.
        % this section estimates all windows seperately and concatinates
        % the residual series to estimate one Cov-matrix.
        application_staticSigma
        return
        
    case {'Sigma-dynamic','all-static'}
        % make static VAR matrix.
        % this section estimates one big window and splits
        % the residual series in subwindows to estimate the Cov-matrix.
        application_staticA
        return
        
    case 'all-dynamic'
        % % % % % % % % % % % % % % Main estimation % % % % % % % % % % % % % %
        %  estimates VAR and Covariance dyanmic.
        k = 0;
        for t = ts
            k = k+1;
            
            X = [Y(:,t-RW+1:t)];
            % Add factors to specification
            Factors = [];
            for jj = 0:FactNumLags
                timeFrame = t-RW-jj + 1:t-jj;
                if k < 4
                    Factors =[Factors;
                        [ Commodity(:, timeFrame); MPx(:, timeFrame); TFP(:, timeFrame)]];
                else
                    Factors =[Factors;
                        [ Commodity(:, timeFrame); MPx2(:, timeFrame); TFP(:, timeFrame)]];
                end
            end
            Z = [ones(1,size(X,2));  Factors; Y(:,t-RW:t-1)];

            factors = size(Z,1) - size(X,1);
            Ta = size(X,2);
            [Est.Ahat, Est.Areg, Est.Sighat, Est.SigReg, RegInfo(t)] =...
                RegularizeApp(X,Z,RegA,RegSigma,deltas,Kfold,[],[]);
            
            % If VAR matrix isnt squared yet, it will be squared in order to
            % calculate FEVDs. Est2 is the version needed for the FEVD
            % function
            if size(Est.Areg,1) == size(Est.Areg,2)
                Est2 = Est;
            else
                [Est2.Areg, Est2.Ahat, Est2.Sighat, Est2.SigReg] = deal(zeros(p*N,p*N));
                Est2.Areg(1:N,:) = Est.Areg(1:end,1+factors:end);
                Est2.Areg(N+1:end,1:N*(p-1)) = eye(N*(p-1));
                Est2.Ahat(1:N,:) = Est.Ahat(1:end,1+factors:end);
                Est2.Ahat(N+1:end,1:N*(p-1)) = eye(N*(p-1));
                Est2.Sighat(1:N,1:N) =  Est.Sighat(1:end,1:end);
                Est2.SigReg(1:N,1:N) =  Est.SigReg(1:end,1:end);
            end
            
            % Calculate regression results
            Nonzero(:,t) = RegInfo(t).dfA;
            NonzeroSig(:,t) = RegInfo(t).dfSig;
            averageShrink(t) = mean(mean(Est2.Areg(1:N,1:N*p)./Est2.Ahat(1:N,1:N*p)));
            averageShrinkSig(t) = mean(mean(Est2.SigReg(1:N,1:N)./Est2.Sighat(1:N,1:N)));
            
            % Ledoit-Wolf test for orthogonality
            FinalSigmaReg(:,:,t) = Est2.SigReg(1:N,1:N);
            [KullbackLeibner(1,t) ,pLWtest(1,t)] = LWtest(FinalSigmaReg(:,:,t),N,Ta);
            [KullbackLeibner(2,t) ,pLWtest(2,t)] = LWtest(Est2.Sighat(1:N,1:N),N,Ta);
            
            % save estimation results for later calculations and plots
            AICcCoef(:,t) = RegInfo(t).AICcCoef;
            AICcSig(:,t) = RegInfo(t).AICcSig;
            MSE(:,:,t) = RegInfo(t).CVlossCoef;
            LL(:,:,t) = RegInfo(t).CVlossSigma;
            eta(t) = RegInfo(t).eta;
            EstSave2(t) = Est2;
            EstSave1(t) = Est;
            FactorsSave{t} = Factors;
        end
end
disp(['Time elapsed in minutes:                ' num2str(toc(firstTic)/60)])

%% calculate FEVDs
clear FEVD DfinalRegIRFfactor
H = 3;
for t = ts
    Factors = FactorsSave{t};
    Est2 =  EstSave2(t) ;
    Est = EstSave1(t);
    
    % FEVD calculation for network analysis
    D = FEVD(Est2,0,H,'FEVD');
    DfinalOLS(:,:,t) = diag(sum(D.estimated(1:N,1:N),2))\D.estimated(1:N,1:N);
    DfinalReg(:,:,t) = diag(sum(D.Reg(1:N,1:N),2))\squeeze(D.Reg(1:N,1:N));
    
    delta = logspace(-1.5,-0.5,100);
    parfor f=1:size(delta,2)
        [loglikelihood(f),~] = KfoldCVsigma(Factors, 12,delta(f),0:1,'normalized GLASSO');
    end
    %     covFactor = cov(Factors');%shrinkLedoit(cov(Factors'),Factors);
    [~,I] = min(loglikelihood);
    
    deltaMatr = delta(I)*(ones(size(Factors,1))-eye(size(Factors,1)));
    SaveStd = std(Factors');
    resTrain = SaveStd'.\Factors;
    SigReg = GraphicalLasso(resTrain',deltaMatr);%, 0, 0, 0, 0, 0);
    covFactor(:,:) = diag(SaveStd)*SigReg*diag(SaveStd);
    
    % Factor spillover estimation
    D = FEVD(Est2,0,12,'IRV');
    DfinalRegIRF(:,:,t) = D.Reg(1:N,1:N);
        
    DfinalRegIRFfactor{t} = (Est.Areg(1:end,(2:size(Factors,1)+1))*(covFactor))*Est.Areg(1:end,(2:size(Factors,1)+1))';Est2.SigReg; 
    DfinalRegIRFfactorPlusSpill = (DfinalRegIRF(:,:,t)/(FinalSigmaReg(:,:,t)))*DfinalRegIRFfactor{t};%D.Reg(1:N,1:N);%
    DfinalRegIRFfactorSpillovers{t} = DfinalRegIRFfactorPlusSpill - DfinalRegIRFfactor{t};
    
    avConnectednessOLS(t) = sum(mean(DfinalOLS(:,:,t)-diag(diag(DfinalOLS(:,:,t)))));
    avConnectednessReg(t) = sum(mean(DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t)))));
end

%% Mean weight calculation and partial R2
k = 0;
for t = ts - DataOffset
    k = k+1;
    meanWeight(:,k)=mean(weights(:,t-RW+1:t),2);
    % partial R2 are weighted by the meanweight per period
    partR2(k) = meanWeight(:,k)'*RegInfo(t+DataOffset).partR2/100;
    partR2spillover(k) = meanWeight(:,k)'*RegInfo(t+DataOffset).partR2spillover/100;
    partR2OLS(k) = meanWeight(:,k)'*RegInfo(t+DataOffset).partR2OLS/100;
    partR2spilloverOLS(k) = meanWeight(:,k)'*RegInfo(t+DataOffset).partR2spilloverOLS/100;
end

%% Output to show in command window
disp(['average Connectedness per period (ols): ' num2str(avConnectednessOLS(ts))])
disp(['average Connectedness per period (reg): ' num2str(avConnectednessReg(ts))])
disp(['Degrees of freedom per period (coef):   ' num2str(Nonzero(2,ts))])
disp(['Degrees of freedom per period (sig):    ' num2str(NonzeroSig(ts))])
disp(['Average Shrinkage (coef):               ' num2str(averageShrink(ts))])
disp(['Average Shrinkage (sig):                ' num2str(averageShrinkSig(ts))])
disp(['LW test of orthogonality:               ' num2str(pLWtest(1,ts))])
disp(['Partial R2 for factors:                 ' num2str(partR2)])
disp(['Partial R2 for spillovers:              ' num2str(partR2spillover)])

AdaptiveEnetOLS = [AICcCoef(1,ts)'; sum(AICcCoef(1,ts))]; AdaptiveEnetEnet = [AICcCoef(2,ts)'; sum(AICcCoef(2,ts))];
Enet = [AICcCoef(3,ts)'; sum(AICcCoef(3,ts))]; OLS =[inf(1,length(ts))'; inf];
Period = ["1972 to 1983"; "1984 to 1995"; "1996 to 2007"; "2008 to 2019" ; "Sum"];
T_aic = table(Period,AdaptiveEnetOLS,AdaptiveEnetEnet,Enet,OLS); disp(T_aic)

Threshold = [AICcSig(1,ts)'; sum(AICcSig(1,ts))]; GLASSO = [AICcSig(2,ts)'; sum(AICcSig(2,ts))];
LedoitWolf = [AICcSig(3,ts)'; sum(AICcSig(3,ts))]; SampleCov =[AICcSig(4,ts)'; sum(AICcSig(4,ts))];
T_aic = table(Period,Threshold,GLASSO,LedoitWolf,SampleCov); disp(T_aic)
%% Export and save data

save(['./Output/ApplicationRun_'+string(month(now))+'_'+ string(day(now)) ])
for t = ts
    table2export = array2table(DfinalReg(:,:,t),'VariableNames',names);
    table2export = [ table(names', 'VariableNames', {'Sectors'})  table2export];
    writetable(table2export,['./Output/FEVD_H' num2str(H) '_from' num2str(year_date(t-144),'%02.f') 'to' num2str(year_date(t),'%02.f') '.csv']);
end

%% generate estimation table
AvCon(1,:) = 100*avConnectednessOLS(ts); AvCon(2,:) = 100*avConnectednessReg(ts);
dfTableA(1,:) = ones(1,length(ts))*(sum(Est.Ahat~=0,[1 2]) ); dfTableA(2,:) = Nonzero(2,ts) ;
dfTableSig(1,:) = ones(1,length(ts))*( N*(N+1)/2); dfTableSig(2,:) = NonzeroSig(ts);
LWtable(1,:) = pLWtest(2,ts); LWtable(2,:) = pLWtest(1,ts);
R2fact(1,:) = 100*partR2OLS; R2fact(2,:) = 100*partR2;
R2spill(1,:) = 100*partR2spilloverOLS; R2spill(2,:) = 100*partR2spillover;

latexTableEstimationOverview(NAIClevel,N,AvCon,dfTableA,dfTableSig,LWtable,R2fact,R2spill)

%% generate AIC table

OutputTable = latexTableAIC(MSE,LL,AICcCoef,AICcSig,ts);

%% Sample Variance
Yq=[ones(N,2) 400*(log(Data(:,4:end))-log(Data(:,1:end-3)))];       %Quarterly equivalent
k = 0;
for t=ts
    k= k+1;
    weightComponent = (meanWeight(:,k)/sum(meanWeight(:,k)));
    sampleVarQuarterly(k) = var(weightComponent'*Yq(:,t-RW+1:t));
    sampleVar(k) = var(weightComponent'*Y(:,t-RW+1:t));
end

%% Decomposition Table
% creates a latex table for the decomposition of effects
decompositionTable(ts,DfinalRegIRF,DfinalRegIRFfactor,DfinalRegIRFfactorSpillovers,FinalSigmaReg,sampleVar,meanWeight)

%% Plots

application_plots;