clear variables structures
% Data source
% load MergedMDPFirmChar
VAR_SET=[13]  %3 is our default from the previous. % 1 includes house prices % 2 & 20 include migration 4 has patents
% 12: small wage VAR. 11: Small house price VAR. 
% 13: migration sample and extra variables in periphery
% 130: As 13, but with more variables in the VAR itself - used for WHICH_IV==5
WHICH_IV=[1]  %default: 1. 5 firm entry rate 
% DO_CLUSTER = -1: block bootstrap
% DO_CLUSTER = -2: spatial block bootstrap
DO_CLUSTER=[-2]
k=[  2 ]
% for DO_OTHER_CHOL=[ -1 0 1  2 ]
% for DO_OTHER_CHOL=[  1  2]
DO_OTHER_CHOL=[   1   ]
DO_WHARTON_LAND=[0 ] % wharton governs whether and how to group MSAs

EndYear = inf;
    
close all

B_REPS=100;

CI_lvls=[5 16 50 84 95];

% use constant base year Bartik weights?
DO_1974_WEIGHTS = false; % true;

DataPath = '../Data/';

DO_BIAS_CORRECTION=1; % 0 or 1. (Not true or false)
CombineDistances=0;

DO_WEIGHTED_PATENTS=1;
% DO_WHARTON_LAND=-1;
if VAR_SET==2 || VAR_SET==20
    EndYear=2012;
end
if VAR_SET==13 || VAR_SET==130  || VAR_SET==131 || VAR_SET==1333
    StartYear=1986;
    if k==3
        StartYear=StartYear+1;
    end
    EndYear=min(EndYear,2013);
    if (DO_WHARTON_LAND)~=1
        load([DataPath 'msa_data_jae']);
    else
        load([DataPath 'msa_data_wharton_jae']);
    end
    DO_WEIGHTED_PATENTS=0;
else
    if (DO_WHARTON_LAND)~=1
        load([DataPath 'msa_data_jae']);
    else
        load([DataPath 'msa_data_wharton_jae']);
    end
%     StartYear=1984;
    StartYear=1980;
%     StartYear=1990;
    DO_WEIGHTED_PATENTS=0;
end

%     StartYear=1980;
if VAR_SET==331
    StartYear=max(StartYear,1980);
end

%     StartYear=1986;
%     EndYear=min(EndYear,2013);

DO_MR_Standardization = false;
POP_ADJUST_Standardization = true;
DO_GOV_EMP=0;

warning('off','optim:fminunc:SwitchingMethod')
h=10;

DO_YEAR_FE=1;
DO_MSA_IRF=1;

USE_TV=0;
DO_WILD=0;

ChoiceOfDistance='geo'; %Options: 'Geo' uses geographical distance
% ChoiceOfDistance='cycle'; %Options: 'cycle' uses Ted Crone style correlation of cycles
% ChoiceOfDistance='mix'; %Options: 'mix' combines geography with state based distance matrix
% ChoiceOfDistance=''; %Options: all other use common state
TruncateDistance=false;

NON_STANDARD_H0=''; %Options: 'cons', 'iid', or default (varying rhos)

Suffix='';

if StartYear>min(vYear);
    Suffix=[Suffix , num2str(StartYear, '_%1.0f')];
end
if EndYear<max(vYear);
    Suffix=[Suffix , num2str(EndYear, '_%1.0f')];
end

if DO_OTHER_CHOL==1
    WhichChol=99;
    Suffix=[Suffix '_Org2'];
elseif DO_OTHER_CHOL==-1
    WhichChol=-99;
    Suffix=[Suffix '_Org1'];
elseif DO_OTHER_CHOL==2
    WhichChol=3;
    Suffix=[Suffix '_UpMR'];
else
    WhichChol=1;
end

Suffix=[Suffix '_split' ];

if DO_WILD
    Suffix=[Suffix '_wild'];
end
%% Config Data
my_config_data_periphery;
%% Beginning of actual code
[T,N]=size(vjob_creation_rate_births(vYear<=EndYear,:));

%% Read distance matrices
my_Read_Distance_Matrices_jae;

%% Split sample into groups based on initial conditions
if (DO_WHARTON_LAND~=1) && (DO_WHARTON_LAND~=-1) && (DO_WHARTON_LAND~=-2)
    NoGroups=1;
   GroupVariable={'log(vdensity)', '_densG', [0 75; 75+1e-10 100]};
elseif DO_WHARTON_LAND==-1
    NoGroups=2;
    GroupVariable={'vpop(1,:)', '_LgP', [0 95;95+1e-10 100]};
elseif DO_WHARTON_LAND==-2
    NoGroups=2;
    GroupVariable={'HDD_m_CDD', '_HeatG', [50+1e-10 100; 0 50]};
else
    NoGroups=2;
    GroupVariable={'v_whartonlandindex(1,:)', '_Wharton', [0 67;67+1e-10 100]};
end

% Adjust Suffix for number of different ex ante groups
if NoGroups==1
    Suffix=[Suffix '_pooled'];
else
    for nn=1:size(GroupVariable,1)
        Suffix=[Suffix GroupVariable{nn,2}];
    end
end

HDD_m_CDD=(HDD-CDD)';

if size(GroupVariable,1)==1
    IndexM = zeros(N,NoGroups);
    grouping_var=eval(GroupVariable{1});
    if NoGroups>1
        for gg=1:NoGroups
            IndexM(:, gg) =  (grouping_var(1,:) >= prctile(grouping_var(1,:), GroupVariable{1,3}(gg,1) )) & ...
                                (grouping_var(1,:) <= prctile(grouping_var(1,:), GroupVariable{1,3}(gg,2) ));
        end
    else
        for gg=1:NoGroups
            IndexM(:, gg) =  (grouping_var(1,:) >= prctile(grouping_var(1,:), 100*(gg-1)/NoGroups)) & ...
                                (grouping_var(1,:) <= prctile(grouping_var(1,:), 100*(gg)/NoGroups));
        end
    end
else
    IndexM = zeros(N,NoGroups^size(GroupVariable,1));
    AuxIndex=cell(size(GroupVariable,1),1);
    
    % loop over the grouping variables
    for vv=1:size(GroupVariable,1)
        % simplify code by having a common grouping variable
        grouping_var=reshape(eval(GroupVariable{vv,1}),1,[]);
        % generate index based on current grouping only
        AuxIndex{vv}=zeros(N, NoGroups);
        for gg=1:NoGroups
            AuxIndex{vv}(:, gg) =  (grouping_var(1,:) >= prctile(grouping_var(1,:), GroupVariable{vv,3}(gg,1) )) & ...
                                (grouping_var(1,:) <= prctile(grouping_var(1,:), GroupVariable{vv,3}(gg,2) ));
        end
    end
    
    % combine into a single index
    % the index consists of zeros and ones
    for vv1=1:size(GroupVariable,1)
        
        % initialize index variable equal to the index for the current variable
        IndexM(:, (vv1-1)*NoGroups+(1:NoGroups)) = AuxIndex{vv1} ;
        
        % loop over all groups except the current group and update with the 
        % other variables
        for vv2=1:(vv1-1) 
            %first, interact index for variable 1 with index for variable 2:
            %Attribute~1 and Attribute~2 both in same group
            IndexM(:, (vv1-1)*NoGroups+(1:NoGroups)) = IndexM(:, (vv1-1)*NoGroups+(1:NoGroups)).*AuxIndex{vv2}; 
        end
        for vv2=vv1+1:size(GroupVariable,1)
            %second, interact index for variable 1 with index for variable 2:
            %Attribute~1 and and NOT Attribute~2 
            IndexM(:, (vv1-1)*NoGroups+(1:NoGroups)) = IndexM(:, (vv1-1)*NoGroups+(1:NoGroups)).*(1-AuxIndex{vv2});
        end
    end
end

fprintf('Overall density')
fprintf([num2str(prctile((vdensity(1,:)), [0 25 50 75 100]), ':: %1.2f') '\n'])
fprintf('Overall HDD-CDD')
fprintf([num2str(prctile(HDD_m_CDD, [0 25 50 75 100]), ':: %1.2f') '\n'])


for nn=1:NoGroups^size(GroupVariable,1)
    fprintf(num2str(nn,'Group %1.0f density'))
    fprintf([num2str(prctile((vdensity(1,IndexM(:,nn)==1)), [0 25 50 75 100]), ':: %1.2f') num2str(sum(IndexM(:,nn)==1), ' :: N=%1.0f') '\n'])
    fprintf(num2str(nn,'Group %1.0f HDD-CDD'))
    fprintf([num2str(prctile((HDD_m_CDD(IndexM(:,nn)==1)), [0 25 50 75 100]), ':: %1.2f') num2str(sum(IndexM(:,nn)==1), ' :: N=%1.0f') '\n'])
    if size(GroupVariable,1)>1
        fprintf(num2str(nn, ['Group %1.0f ' GroupVariable{2,2} ]))
        fprintf([num2str(prctile(HDD_m_CDD(IndexM(:,nn)==1), [0 25 50 75 100]), ':: %1.2f') num2str(sum(IndexM(:,nn)==1), ' :: N=%1.0f')  '\n']);
    end
end

if ~isequal(sum(IndexM(:)), N)
    error('sample split does not work -- lack of a tie-breaking rule?');
end

%%
NoGroups=size(IndexM,2);
if DO_YEAR_FE
    px=p*k;
else
    error('not set up without Year FE');
end

pp = size(PerNames,1);
if DO_YEAR_FE
    ppx=k + p*(k+1);
else
    error('not set up without Year FE');
end

if WHICH_IV==1 || WHICH_IV==101 || WHICH_IV==1010 
        ShockLabels={'entry job creation';'agg labor demand'; };    
%         ShockLabels={'agg labor demand'};    
        if WHICH_IV==1
            Suffix=[Suffix '_Identification'];
        elseif WHICH_IV==101
           Suffix=[Suffix '_op'];
        elseif WHICH_IV==1010
           Suffix=[Suffix '_op10'];
        end
        % Switch out instruments
        if ~isempty(strmatch('vestabs_entry',VarNames))
            if VAR_SET>=10 && VAR_SET<=19
                OriginalPos=strmatch('vestabs_entry',VarNames);
                AltPos=strmatch('vjob_creation',VarNames);
            else
                OriginalPos=strmatch('vestabs_entry_rate',VarNames);
                AltPos=strmatch('vjob_creation_rate_births',VarNames);
            end
            Tmp=VarNames{OriginalPos};
            VarNames{OriginalPos}=VarNames{AltPos};
            VarNames{AltPos}=Tmp;

            Tmp=VarLabels{OriginalPos};
            VarLabels{OriginalPos}=VarLabels{AltPos};
            VarLabels{AltPos}=Tmp;

            Tmp=VarLabelsShort{OriginalPos};
            VarLabelsShort{OriginalPos}=VarLabelsShort{AltPos};
            VarLabelsShort{AltPos}=Tmp;

            Tmp=VarScale(OriginalPos);
            VarScale(OriginalPos)=VarScale(AltPos);
            VarScale(AltPos)=Tmp;
        end
end
    
if WHICH_IV==5 % use with VAR_SET=130
        % swap out alternative instruments for those originally included
        ShockLabels={'entry rate';'agg labor demand'; };    
        Suffix=[Suffix '_IVfirm'];
 
        OriginalPos=strmatch('vjob_creation_rate_births',VarNames);
        AltPos=strmatch('vfirm_entry_rate',VarNames);

        Tmp=VarNames{OriginalPos};
        VarNames{OriginalPos}=VarNames{AltPos};
        VarNames{AltPos}=Tmp;

        Tmp=VarLabels{OriginalPos};
        VarLabels{OriginalPos}=VarLabels{AltPos};
        VarLabels{AltPos}=Tmp;

        Tmp=VarLabelsShort{OriginalPos};
        VarLabelsShort{OriginalPos}=VarLabelsShort{AltPos};
        VarLabelsShort{AltPos}=Tmp;

        Tmp=VarScale(OriginalPos);
        VarScale(OriginalPos)=VarScale(AltPos);
        VarScale(AltPos)=Tmp;

end

if WHICH_IV==-1 
        ShockLabels={'entry job creation';'agg labor demand'};    
%         ShockLabels={'agg labor demand'};    
        Suffix=[Suffix '_JCB'];
        % Switch out instruments
        if VAR_SET>=10 && VAR_SET<=19
            OriginalPos=strmatch('vestabs_entry',VarNames);
            AltPos=strmatch('vjob_creation',VarNames);
        else
            OriginalPos=strmatch('vestabs_entry_rate',VarNames);
            AltPos=strmatch('vjob_creation_rate_births',VarNames);
        end
        Tmp=VarNames{OriginalPos};
        VarNames{OriginalPos}=VarNames{AltPos};
        VarNames{AltPos}=Tmp;
        
        Tmp=VarLabels{OriginalPos};
        VarLabels{OriginalPos}=VarLabels{AltPos};
        VarLabels{AltPos}=Tmp;

        Tmp=VarLabelsShort{OriginalPos};
        VarLabelsShort{OriginalPos}=VarLabelsShort{AltPos};
        VarLabelsShort{AltPos}=Tmp;

        Tmp=VarScale(OriginalPos);
        VarScale(OriginalPos)=VarScale(AltPos);
        VarScale(AltPos)=Tmp;
end


if k~=2
        Suffix=[Suffix '_k' num2str(k)];        
end

if ~isempty(NON_STANDARD_H0)
    Suffix=[Suffix '-' NON_STANDARD_H0];
end

switch upper(ChoiceOfDistance)
    case 'GEO'
        Suffix=[Suffix '_geo'];        
    case 'MIX'
        Suffix=[Suffix '_mix'];        
    case 'POP'
        Suffix=[Suffix '_pop'];                
    case 'CYCLE'
        Suffix=[Suffix '_cyc'];
    otherwise
        Suffix=[Suffix '_other'];
end

if TruncateDistance
    Suffix=[Suffix '_trc'];        
end

if CombineDistances==1
    Suffix=[Suffix '_comb'];    
end

%% Check for NaN data missing (at beginning only)

MissingObs=0;
for p_idx=1:p
    [rr,~]=eval(['find(isnan(' VarNames{p_idx} '(1:T,:)));']);
    MissingObs=max([MissingObs,max(rr)]);
end

StartObs=k+1+MissingObs;
if vYear(StartObs)<StartYear
    StartObs=find(vYear==StartYear);
    MissingObs=StartObs-(k+1);
end

%% Assemble the data

% Instrument(s)
if WHICH_IV==1
%     return
    if ~DO_1974_WEIGHTS
        Z(1,:,:)=vZit_Bartik_jc(StartObs:T,:); %vZit_jcb
            VarNames=[VarNames; 'vZit_Bartik_jc'];
            VarLabels=[VarLabels; 'Bartik: Entrant''s job creation'];
        Z(2,:,:)=vZit_Bartik(StartObs:T,:);
        VarNames=[VarNames; 'vZit_Bartik'];
    else
        Z(1,:,:)=vZit_Bartik_jc1974(StartObs:T,:); %vZit_jcb
        VarNames=[VarNames; 'vZit_Bartik_jc1974'];
        VarLabels=[VarLabels; 'Bartik: Entrant''s job creation'];
        Z(2,:,:)=vZit_Bartik1974(StartObs:T,:);
        VarNames=[VarNames; 'vZit_Bartik1974'];
    end
    
    VarLabels=[VarLabels; 'Bartik: Employment'];
elseif WHICH_IV==5
    if ~DO_1974_WEIGHTS
        Z(1,:,:)=vZit_Bartik_firm(StartObs:T,:); %vZit_jcb
        Z(2,:,:)=vZit_Bartik(StartObs:T,:);
        VarNames=[VarNames; 'vZit_Bartik_firm'; 'vZit_Bartik'];
    else
        Z(1,:,:)=vZit_Bartik_firm1974(StartObs:T,:); %vZit_jcb
        Z(2,:,:)=vZit_Bartik1974(StartObs:T,:);
        VarNames=[VarNames; 'vZit_Bartik_firm1974'; 'vZit_Bartik1974'];
    end
    
    VarLabels=[VarLabels; 'Bartik: Firm entry'];
    VarLabels=[VarLabels; 'Bartik: Employment'];      
else
    error('IV option specified that is currently not set up')
end

pz=size(Z,1);
if ~exist('pv1', 'var') 
    pv1=pz;
end
if pz==1
    WhichChol=1; % with a single instrument
end

% VAR variables
Y=NaN(p,T-k-MissingObs,N);
Yall=NaN(p,T-MissingObs,N);

% Periphery variables
pY=NaN(pp,T-k-MissingObs,N);
pYall=NaN(pp,T-MissingObs,N);

% Core
for p_idx=1:p
    Y(p_idx,:,:)=VarScale(p_idx)*eval(['reshape(' VarNames{p_idx} '(StartObs:T,:),[1,T-k-MissingObs,N]);']);
    Yall(p_idx,:,:)=VarScale(p_idx)*eval(['reshape(' VarNames{p_idx} '(StartObs-k:T,:),[1,T-MissingObs,N]);']);
end

% Periphery
for p_idx=1:pp
    pY(p_idx,:,:)=PerScale(p_idx)*eval(['reshape(' PerNames{p_idx} '(StartObs:T,:),[1,T-k-MissingObs,N]);']);
    pYall(p_idx,:,:)=PerScale(p_idx)*eval(['reshape(' PerNames{p_idx} '(StartObs-k:T,:),[1,T-MissingObs,N]);']);
end


X=NaN(px,T-k-MissingObs,N);
for k_idx=1:k
    for p_idx=1:p
        X((k_idx-1)*p+p_idx,:,:)=VarScale(p_idx)*eval(['reshape(' VarNames{p_idx} '( (StartObs-k_idx):(T-k_idx), :),[1,T-k-MissingObs,N]);']);
    end
end
pX=NaN(ppx,T-k-MissingObs,N,pp); % for every p, the set of regressors differs
for k_idx=1:k+1
    % include core variables
    for p_idx=1:p
        % this has pp duplicates for easier data manipulation
        pX(k +(k_idx-1)*p+p_idx,:,:,:)=VarScale(p_idx)*eval(['reshape( repmat(' VarNames{p_idx} '( (StartObs-(k_idx-1) ):(T-(k_idx-1)), :), [1,1,1,pp]),[1,T-k-MissingObs,N,pp]);']);
    end
    if k_idx<=k
        for p_idx=1:pp
            pX(k_idx ,:,:,p_idx)=PerScale(p_idx)*eval(['reshape(' PerNames{p_idx} '( (StartObs-k_idx ):(T-k_idx), :),[1,T-k-MissingObs,N,1]);']);
        end
    end
end


ScaledCycleCorr=CycleCorr;
ScaledStateCorr=StateCorr;
ScaledDistanceCorr=DistanceCorr;
for n_idx=1:N;
    if sum(CycleCorr(n_idx,:))>0
        if TruncateDistance
            ScaledCycleCorr(n_idx,ScaledCycleCorr(n_idx,:)<prctile(ScaledCycleCorr(n_idx,:),75))=0;
        end
        ScaledCycleCorr(n_idx,:)=ScaledCycleCorr(n_idx,:)/sum(ScaledCycleCorr(n_idx,:));
    end
    if sum(StateCorr(n_idx,:))>0
        if TruncateDistance
            ScaledStateCorr(n_idx,ScaledStateCorr(n_idx,:)<prctile(ScaledStateCorr(n_idx,:),75))=0;
        end
        ScaledStateCorr(n_idx,:)=ScaledStateCorr(n_idx,:)/sum(ScaledStateCorr(n_idx,:));
    end
        
    if sum(DistanceCorr(n_idx,:))>0
        if TruncateDistance
            ScaledDistanceCorr(n_idx,ScaledDistanceCorr(n_idx,:)<prctile(ScaledDistanceCorr(n_idx,:),75))=0;
        end
        ScaledDistanceCorr(n_idx,:)=ScaledDistanceCorr(n_idx,:)/sum(ScaledDistanceCorr(n_idx,:));
    end
end

% return;
%% IV
% display( [max(abs(eig(ScaledStateCorr))), max(abs(eig(ScaledRegionCorr))), max(abs(eig(ScaledWeightCorr))), max(abs(eig(ScaledDistanceCorr)))]-1 );

DO_VERBOSE=1;
DO_IRF=1;
NoIVLags=k+2; %>=max(2,ModelStructure.k+1)
Params=var2struct(p,pv1,pz,px, pp,ppx,DO_YEAR_FE, N,T, k, MissingObs, StartObs, WHICH_IV,DO_VERBOSE,NoIVLags,DO_IRF,WhichChol,CombineDistances);
switch upper(ChoiceOfDistance)
    case 'GEO'
        DistanceMatrix1=ScaledDistanceCorr;
        DistanceMatrix2=ScaledStateCorr;
    case 'POP'
        DistanceMatrix1=ScaledDistancePopCorr;
        DistanceMatrix2=ScaledCycleCorr;
    case 'CYCLE'
        DistanceMatrix1=ScaledCycleCorr;
        DistanceMatrix2=ScaledCycleCorr;
    case 'TMP'
        TMP_DISTANCE=zeros(N,N);
        TGT_ROW = [0 N-1:-1:1];
        TGT_ROW =TGT_ROW /sum(TGT_ROW );
        for nn=1:N
            TMP_DISTANCE(nn, nn:N)=TGT_ROW(1:N-nn+1);
            TMP_DISTANCE(nn, 1:nn-1)=TGT_ROW(N-nn+2:end);
        end
        DistanceMatrix1=ScaledDistanceCorr;
        DistanceMatrix2=ScaledStateCorr;
    otherwise
        DistanceMatrix1=ScaledStateCorr;
        DistanceMatrix2=ScaledCycleCorr;
end

% Call the estimation function
[K_IV, K_IV_combined, IV_IRF_opt, IV_IRF_iid, IV_IRF_cons, IV_IRF_alt, logLL, E_opt, E_iid, E_cons, E_alt, Ap_IV, pE_opt, pE_iid, pE_cons, pE_alt, RHOS, DATA, DIST_WGT] ...
    = my_Spatial_Error_split_Periphery_jae(Y, X, pY, pX, Z, Yall,pYall, DistanceMatrix1, DistanceMatrix2, Params, IndexM);

% Normalize sign on the entire employment response, not the
% employment to population ratio response
if POP_ADJUST_Standardization
   POP_ADJUST=zeros(pv1,1);
   POP_IDX = strcmp(VarNames, 'dlog_pop');
   for p_idx=1:pv1
        if ~isempty(strfind(lower(VarNames{p_idx}), 'pop'))
            POP_ADJUST(p_idx,1)=1;
        end
   end
end


try
    Warn = warning('query','last');
    WarnID = Warn.identifier;
    warning('off',WarnID);
catch
    
end
%% First stage F-stat
for gg=1:NoGroups
    % my_first_stage_F_v2 already adjusts for number of MSAs in group
    Original_IV_F{gg}=my_first_stage_F_jae(E_opt(:,:,1==IndexM(:,gg)), T-k-MissingObs, Params, IV_IRF_opt{gg,1}(1:p,:),DATA.V_opt{gg});
    display(Original_IV_F{gg}');
end


%% Diagnostic for residual correlation

ResCorr=zeros(N,N,p+pz);
close all
T_eff = size(squeeze( E_opt(p_idx,:,:) ),1);
W_State_eff = sum(ScaledStateCorr(:));
IndividualCorr = zeros(N, p+pz,3);
for p_idx=1:(p+pz)
    tmpCov = cov( squeeze( E_opt(p_idx,:,:) )  , 1);    
    for nn=1:N
       IndividualCorr(nn, p_idx, 1) = corr(squeeze(E_opt(p_idx, :, nn))' , squeeze(E_opt(p_idx, :, :))*DistanceMatrix1(nn,:)'); 
       IndividualCorr(nn, p_idx, 2) = corr(squeeze(E_cons(p_idx, :, nn))' , squeeze(E_opt(p_idx, :, :))*DistanceMatrix1(nn,:)'); 
       IndividualCorr(nn, p_idx, 3) = corr(squeeze(E_iid(p_idx, :, nn))' , squeeze(E_opt(p_idx, :, :))*DistanceMatrix1(nn,:)'); 
    end
    
    
%     % pairwise correlations: Across all MSAs, correlation over time
%     tmpCorr=(corr( squeeze( E_opt(p_idx,:,:) )  ))';
%     
%     % Z-score based on StateCorr matrix
%     z_score1=((sum(reshape(tril(atanh(tmpCorr(StateCorr==1))), 2*sum(sum( tril(StateCorr==1 ))),1))/(2*sum(sum( tril(StateCorr==1 ))))));
%     ResCorr(:,:,p_idx)=tmpCorr;
% %     display(mean(mean(tmpCorr(StateCorr==1))));
%     TmpHist=tril( tmpCorr(StateCorr==1) );
% %     sum(sum( tril(StateCorr==1 )))
%     figure;
%     hist(TmpHist(:) ,51); % histogram based on common state correlation
%     hold all;
%     plot( tanh(z_score1)* [1,1], ylim, '-r', 'LineWidth', 3)
%     tmpCorr=(corr( squeeze( E_iid(p_idx,:,:) )  ))';
%     % same z-score, but without spatial filtering
%     z_score3=((sum(reshape(tril(atanh(tmpCorr(StateCorr==1))), 2*sum(sum( tril(StateCorr==1 ))),1))/(2*sum(sum( tril(StateCorr==1 ))))));
%     plot( tanh(z_score3)* [1,1], ylim, '-g', 'LineWidth', 3)
%     tmpCorr=(corr( squeeze( E_cons(p_idx,:,:) )  ))';
%     % same z-score, but with common rho in spatial filtering
%     z_score2=((sum(reshape(tril(atanh(tmpCorr(StateCorr==1))), 2*sum(sum( tril(StateCorr==1 ))),1))/(2*sum(sum( tril(StateCorr==1 ))))));
%     plot( tanh(z_score2)* [1,1], ylim, '-m', 'LineWidth', 3)
%     plot( 0.064* [1,1], ylim, '-k', 'LineWidth', 3)
%     plot( -0.064* [1,1], ylim, '-k', 'LineWidth', 3)
%     title(VarLabels{p_idx});
%     set(findall(gcf,'type','axes'),'fontsize',18);
%     set(findall(gcf,'type','text'),'fontSize',18);    
%     axis tight;
%     axis([-1 1 ylim]);
%     
% %     display(abs([z_score1 z_score2 z_score3]))
%     display(tanh([z_score1 z_score2 z_score3]));
%     
% %     display(sum(reshape(tril(tmpCorr(StateCorr==1)), 2*sum(sum( tril(StateCorr==1 ))),1))/sqrt(2*sum(sum( tril(StateCorr==1 )))));
% %     print('-depsc', ['Graphs/Hist_State_' VarNames{p_idx}  Suffix '.eps']);
end

IndividualCorr_t = IndividualCorr .* sqrt( (T_eff-2) ./ (1-IndividualCorr.^2));
for ii=1:3
    for p_idx = 1:(p+pz)
        disp([ii, prctile(abs(IndividualCorr_t(:, p_idx, ii)), [50 70 80 90 95 97.5]), sum(abs(IndividualCorr_t(:,p_idx, ii))>=tinv([0.95 0.975], T_eff-2))/N]);
    end
end

%% Code for wild bootstrap
% if BS_BLOCK_LENGTH~=1
%     Suffix = [Suffix '_Block' num2str(BS_BLOCK_LENGTH, '%1.0f')];
% end
if DO_CLUSTER==1
    Suffix = [Suffix '_clstr'];
elseif DO_CLUSTER==-1
    Suffix = [Suffix '_blck'];
elseif DO_CLUSTER==-2
    Suffix = [Suffix '_mblck'];
elseif DO_CLUSTER==-3
    Suffix = [Suffix '_nblck'];
end
if DO_1974_WEIGHTS
    Suffix = [Suffix '_74'];
end


DO_TEST=0;
if DO_TEST
    B_REPS=1;
end

if ispc && ~DO_TEST
    BootstrapData=matfile(['BootstrapD' Suffix '.mat'], 'writable', true);
end

rng(1);

IRF_pooled_b=cell(2,NoGroups);
IRF_alt_b=cell(3,NoGroups);
IRF_CF_b=cell(1,NoGroups);
for cc=1:2
    for gg=1:NoGroups
        IRF_pooled_b{cc,gg}=zeros(p+pp,pv1,h+1,B_REPS);
    end
end
for cc=1:3
    for gg=1:NoGroups
        IRF_alt_b{cc,gg}=zeros(p+pp,pv1,h+1,B_REPS);
    end
end
for gg=1:NoGroups
    IRF_CF_b{1,gg}=zeros(p+pp,pv1,h+1,B_REPS);
end

% Opt_RHOS_b=zeros(p+pz,B_REPS);
% Cons_RHOS_b=zeros(1,B_REPS);
MaxEig_b=zeros(NoGroups,B_REPS);
logLL_b=zeros(3,B_REPS);

ResCorr_Region=zeros(p+pz+pp,B_REPS);
ResCorr_State=zeros(p+pz+pp,B_REPS);
% ResCorr_State=zeros(p,B_REPS);
% ResCorr_Region=zeros(p,B_REPS);

% Transform from VAR(1) form
A_IV=cell(3,1);
for gg=1:NoGroups
    A_IV{gg}=zeros(k*p,p);
    for kk = 1:k
        A_IV{gg}( (kk-1)*p + (1:p), :)=K_IV_combined{gg}(1:p,(kk-1)*(pp+p)+(1:p))'; 
    end
end
% Initialize
Y_b=zeros(size(Y));
Z_b=zeros(size(Z));
pY_b=zeros(size(pY));

Params.DO_VERBOSE=0;
p_IV=NoIVLags*p;
p_IVp=NoIVLags + p* (1+k);

% ErrorCorr_opt=zeros(p+pz,T-MissingObs-NoIVLags-1);
% ErrorCorr_cons=zeros(p+pz,T-MissingObs-NoIVLags-1);
% ErrorCorr_iid=zeros(p+pz,T-MissingObs-NoIVLags-1);
% 
% for t_idx=1:T-MissingObs-k
% for p_idx=1:(p+pz)
%     if p_idx<=p
%         ErrorCorr_iid(p_idx,t_idx)=corr(squeeze(DATA.Res_Lvl(p_idx,t_idx,:)),squeeze(E_iid(p_idx,t_idx,:)));
%         ErrorCorr_cons(p_idx,t_idx)=corr(squeeze(DATA.Res_Lvl(p_idx,t_idx,:)),(squeeze(E_cons(p_idx,t_idx,:))'/(speye(N)-RHOS(p_idx,2)*ScaledStateCorr'))');
%         ErrorCorr_opt(p_idx,t_idx)=corr(squeeze(DATA.Res_Lvl(p_idx,t_idx,:)),(squeeze(E_opt(p_idx,t_idx,:))'/(speye(N)-RHOS(p_idx,1)*ScaledStateCorr'))');
%     else
%         ErrorCorr_iid(p_idx,t_idx)=corr(squeeze(DATA.Z(p_idx-p,t_idx,:)),squeeze(E_iid(p_idx,t_idx,:)));
%         ErrorCorr_cons(p_idx,t_idx)=corr(squeeze(DATA.Z(p_idx-p,t_idx,:)),(squeeze(E_cons(p_idx,t_idx,:))'/(speye(N)-RHOS(p_idx,2)*ScaledStateCorr'))');
%         ErrorCorr_opt(p_idx,t_idx)=corr(squeeze(DATA.Z(p_idx-p,t_idx,:)),(squeeze(E_opt(p_idx,t_idx,:))'/(speye(N)-RHOS(p_idx,1)*ScaledStateCorr'))');
%     end
% end
% end

GroupedRHO_cons=RHOS(1:(p+pz+pp),NoGroups+(1:NoGroups));
GroupedRHO_opt =RHOS(1:(p+pz+pp), 1:NoGroups );

% Needed for conversion from within-cell indices to indices on 1, ..., N;
FindIndexM=cell(NoGroups,1);
for gg=1:NoGroups
%     FindIndexM{gg}=repmat( find(IndexM(:,gg)') , T-k-MissingObs, 1) ;
    FindIndexM{gg}= find(IndexM(:,gg)') ;
end

if ~DO_TEST
    if DO_BIAS_CORRECTION==1
        B_REPS_current=floor(B_REPS/2);
    else
        B_REPS_current=B_REPS;
    end
else
    B_REPS_current=1;
end

EffectiveObs = T-k-MissingObs; %ceil( (T-k-MissingObs) / BS_BLOCK_LENGTH );
% BlocksIntoObs = ones(T-k-MissingObs,EffectiveObs);
% BlockToObsMap = @(block_indx) t;
% for bb=1:EffectiveObs 
%     BlocksIntoObs(:,bb)= ceil( (1:T-k-MissingObs)' / BS_BLOCK_LENGTH )==bb;
% end            
switch upper(NON_STANDARD_H0)
    case 'CONS'
        DistanceMatrix_combined = DistanceMatrix1 * DIST_WGT(2) + DistanceMatrix2 * (1-DIST_WGT(2)); 
    otherwise
        DistanceMatrix_combined = my_weighted_Distance_matrix(DistanceMatrix1, DistanceMatrix2, DIST_WGT(1,:), IndexM, NoGroups);
end
% for future reference, compute the closest neighbors for each MSA.
Neighbor_Matrix= zeros(N,N);
for nn=1:N
    [check, idx]=sort(DistanceMatrix_combined(nn,:), 'descend');
    Neighbor_Matrix(nn, :)=idx;
    if isequal(upper(ChoiceOfDistance),'TMP')
        [check, idx]=sort(TMP_DISTANCE(nn,:), 'descend');
        Neighbor_Matrix(nn, :)=idx;
    end
end

for BIAS_CORRECTION=DO_BIAS_CORRECTION:-1:0
    BootstrapData.K_IV_combined=zeros([size(K_IV_combined{1}),NoGroups, B_REPS_current]);
    BootstrapData.A_IV=zeros(p*k,p,NoGroups,B_REPS_current);
    BootstrapData.RHOS_opt=zeros(p+pz+pp,NoGroups,B_REPS_current);
    BootstrapData.RHOS_cons=zeros(1,NoGroups,B_REPS_current);
    BootstrapData.DIST_WGT=zeros(2,NoGroups,B_REPS_current);
    BootstrapData.V_opt=zeros(p,p,NoGroups,B_REPS_current);
    BootstrapData.Ap_IV=zeros(k+(1+k)*p,pp,NoGroups,B_REPS_current);
    BootstrapData.pV_opt=zeros(pp,pp,NoGroups,B_REPS_current);

    BootstrapData.IV_IRF_opt=zeros(p+pp,pv1,NoGroups,B_REPS_current);
    BootstrapData.IV_IRF_opt_noDir=zeros(p+pp,pv1,NoGroups,B_REPS_current);
    BootstrapData.IV_IRF_opt_noInd=zeros(p+pp,pv1,NoGroups,B_REPS_current);
    BootstrapData.IV_IRF_opt_Neither=zeros(p+pp,pv1,NoGroups,B_REPS_current);
    BootstrapData.IV_IRF_iid=zeros(p+pp,pv1,NoGroups,B_REPS_current);
 
    BootstrapData.IV_F_opt=zeros(pz,NoGroups,B_REPS_current);
    
    for b_idx=1:B_REPS_current
        if DO_WILD
            if DO_TEST~=1
                SIGN_FLIP=sign(rand(EffectiveObs,1)-0.5)*ones(1,N);
            else
                SIGN_FLIP=ones( EffectiveObs,N );
            end
        else
            if DO_TEST~=1
                % old: IDX=ceil(rand( (T-k-MissingObs),N) * (T-k-MissingObs) );
                IDX=zeros( EffectiveObs , N);
                switch DO_CLUSTER
                    case 1
                        for gg=1:NoGroups
                            % randomize index on old matrix
                            JJ_smpl= ceil(rand( 1,sum(IndexM(:,gg))) * sum(IndexM(:,gg)) );
                            IDX(:,1==IndexM(:,gg)) = (1:(T-k-MissingObs))' * FindIndexM{gg}(1,JJ_smpl);
    %                         % index is on elements within group gg
    %                         [II,JJ]=ind2sub( [EffectiveObs, sum(IndexM(:,gg))], aux_IDX );
    %                         % transform to index on [1, ..., N]
    %                         JJnew=reshape(FindIndexM{gg}(1,JJ), EffectiveObs, sum(IndexM(:,gg)));
    %                         % transform to linear index on full E matrix
    %                         IDX(:,1==IndexM(:,gg)) = sub2ind([EffectiveObs, N], II, JJnew);
    %     %                     setdiff(unique(FindIndexM{gg}(:))',unique(reshape(IDX(:,1==IndexM(:,gg)),1,[])))
                        end
                    case -1 % block bootstrap
                        BlockLength = 3;
                        aux_IDX = zeros( floor(EffectiveObs/BlockLength) , N);
                        for gg=1:NoGroups
                            % randomize index on old matrix
                            
                            % take one extra draw to when blocks could be
                            % too short
                            aux_IDX = ceil(rand( ceil(EffectiveObs/BlockLength),sum(IndexM(:,gg))) * floor(EffectiveObs/BlockLength)*sum(IndexM(:,gg)) );
                            % index is on elements within group gg
                            [aux_II,aux_JJ]=ind2sub( [floor(EffectiveObs/BlockLength), sum(IndexM(:,gg))], aux_IDX );
                            % last block could get sampled multiple times.
                            % In that 
                            II=nan(EffectiveObs, sum(IndexM(:,gg)));
                            JJ=nan(EffectiveObs, sum(IndexM(:,gg)));
                            for ii=1:sum(IndexM(:,gg))
                                ObsCount = 0;
                                for tt=1:ceil(EffectiveObs/BlockLength)
                                    periods = my_BlockPeriodMap(aux_II(tt,ii),EffectiveObs,BlockLength);
                                    II(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), ii)=periods(1: min(length(periods), EffectiveObs-ObsCount) )';
                                    JJ(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), ii)=aux_JJ(tt,ii);
                                    ObsCount = ObsCount + length(periods);
                                end
                            end
                            % transform to index on [1, ..., N]
                            JJnew=reshape(FindIndexM{gg}(1,JJ), EffectiveObs, sum(IndexM(:,gg)));
                            % transform to linear index on full E matrix
                            IDX(:,1==IndexM(:,gg)) = sub2ind([EffectiveObs, N], II, JJnew);
        %                     setdiff(unique(FindIndexM{gg}(:))',unique(reshape(IDX(:,1==IndexM(:,gg)),1,[])))
                        end
                    case -2 % block bootstrap in space and time
                        
                        % code as in previous case, but randomly pick MSAs
                        % afterward and pick shock block from neighbors
                        %
                        % Note, however, that to avoid problems of small
                        % group, assumption added tht e iid across greoups
                        BlockLength = 3;
                        aux_IDX = zeros( floor(EffectiveObs/BlockLength) , N);
%                         for gg=1:NoGroups
                            % randomize index on old matrix
                            
                            % take one extra draw to when blocks could be
                            % too short
                            aux_IDX = ceil(rand( ceil(EffectiveObs/BlockLength),N) * floor(EffectiveObs/BlockLength)*N );
                            % index is on elements within group gg
                            [aux_II,aux_JJ]=ind2sub( [floor(EffectiveObs/BlockLength), N], aux_IDX );
                            % last block could get sampled multiple times.
                            % In that 
                            II=nan(EffectiveObs, N);
                            JJ=nan(EffectiveObs, N);
                            for ii=1:N
                                ObsCount = 0;
                                for tt=1:ceil(EffectiveObs/BlockLength)
                                    periods = my_BlockPeriodMap(aux_II(tt,ii),EffectiveObs,BlockLength);
                                    II(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), ii)=periods(1: min(length(periods), EffectiveObs-ObsCount) )';
                                    JJ(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), ii)=aux_JJ(tt,ii);
                                    ObsCount = ObsCount + length(periods);
                                end
                            end
                            % transform to index on [1, ..., N]
                            JJnew=reshape( JJ, EffectiveObs, N);
                            JJold=JJnew;
                            IsReplaced=false(size(JJnew));
                            % Now substitute clusters
                            reshuffled_idx = randsample(N, N);
                            % now, for each reshuffled idx, find the
                            % 6=354^(1/3)-1 closest neighbors
                            SpatialBlock=6;
                            used_msas_no=0;
                            distinct_msa_no=0;
                            used_msas_truefalse = false(N,1);
                            while any(~used_msas_truefalse)
                                % pick the closest neighbors that haven't
                                % been used, yet
%                                 
%                                 %first, find the first unused MSA
%                                 % need extra map here
%                                 [sharedVals,idxsIntoA] = intersect(reshuffled_idx,1:N,'stable');
                                % reshuffled_idx serves as its own map
                               current_msa=find(~used_msas_truefalse(reshuffled_idx), 1);
                               % this is not working. Why? I am eliminating
                               % the columns, but they are not in order.
                               idx=intersect(Neighbor_Matrix(reshuffled_idx(current_msa),:), find(~used_msas_truefalse), 'stable');
                               %for all neighbors, replace 
                               LL=min(length(idx),SpatialBlock);
                               JJnew(:,idx(1:LL))=repmat(JJnew(:,reshuffled_idx(current_msa)), 1, LL);
%                                IsReplaced(JJold==idx(jj))=true;
                               % mark as "used"
                               used_msas_truefalse(idx(1:LL))=true;
%                                for jj=1:min(length(idx),6)
%                                    % use same draw as for current msa
%                                end
                               distinct_msa_no=distinct_msa_no+1;
                               % mark own MSA as "used"
                               used_msas_truefalse(reshuffled_idx(current_msa))=true;
                               IsReplaced(JJold==reshuffled_idx(current_msa))=true;
                               % adjust count of reused MSAs
                               used_msas_no=used_msas_no+SpatialBlock+1;
                            end
%                             error('the problem is that not all MSAs get used')
                            % good: not all MSAs have been used. Bad: 

                                % transform to linear index on full E matrix
                            IDX = sub2ind([EffectiveObs, N], II, JJnew);
        %                     setdiff(unique(FindIndexM{gg}(:))',unique(reshape(IDX(:,1==IndexM(:,gg)),1,[])))
                    case -3 % block bootstrap in space and time with variable block size
                        
                        % Step 1: Draw indices of time blocks and spatial
                        % blocks
                        % inner loop: Randomly pick groups of target MSAs
                        % from those not picked, yet
                        % Pick with replacement blocks of tau periods
                        % jointly for X closest neighbors to assign to
                        % target MSAs and year.
                        
                        % idea: target MSA and its X closest neighbors can
                        % assigned the same draws as the origin MSA and its
                        % X closest neighbors
                        BlockLength = 3;
%                         aux_IDX = zeros( floor(EffectiveObs/BlockLength) , N);

                        % randomize index on old matrix

                        % take one extra draw too many in time because
                        % blocks could be too short        
                        aux_IDX = ceil(rand( ceil(EffectiveObs/BlockLength),N) ...
                                            * floor(EffectiveObs/BlockLength)*N );
                        % index is on elements within group gg
                        [aux_II,aux_JJ]=ind2sub( [floor(EffectiveObs/BlockLength), N ], aux_IDX );
                        % Expand along the time dimension
                        II=nan(EffectiveObs, N);
                        JJ=nan(EffectiveObs, N);
                        
                        % Now substitute clusters
                        reshuffled_idx = randsample(N, N);
                        % now, for each reshuffled idx, find the
                        % 6=354^(1/3)-1 closest neighbors
                        SpatialBlock=6;
                        used_msas_no=0;
                        distinct_msa_no=0;
                        used_msas_truefalse = false(N,1);
                        replaced_ct=0;
                        while any(~used_msas_truefalse)
                            ObsCount = 0;
                            
                            replaced_ct=replaced_ct+1;
                            % pick the closest neighbors that haven't
                            % been used, yet
%                                 
%                                 %first, find the first unused MSA
%                                 % need extra map here
%                                 [sharedVals,idxsIntoA] = intersect(reshuffled_idx,1:N,'stable');
                            % reshuffled_idx serves as its own map
                           current_msa=find(~used_msas_truefalse(reshuffled_idx), 1);
                           
                           
                           % this is not working. Why? I am eliminating
                           % the columns, but they are not in order.
                           target_msa_idx=intersect([reshuffled_idx(current_msa), Neighbor_Matrix(reshuffled_idx(current_msa),:)], find(~used_msas_truefalse), 'stable');
                           % check how many MSAs are left to be replaced 
                           LL=min(length(target_msa_idx),SpatialBlock+1);
                           % mark target MSAs as "used"
                           used_msas_truefalse(target_msa_idx(1:LL))=true;
                           % adjust count of used MSAs
                           used_msas_no=used_msas_no+LL;
                           
                           for tt=1:ceil(EffectiveObs/BlockLength)
                               core_origin_msa = ceil(rand(1,1)*N); % draw "core" origin MSA at random with replacement
                               origin_msa_idx=[core_origin_msa, Neighbor_Matrix(core_origin_msa, 1:LL-1)]; %assign closest neighbors to the same block
                               origin_block = ceil(rand(1,1)*floor(EffectiveObs/BlockLength)); % draw "core" origin block at random with replacement

                               % which periods correspond to the chosen block?
                               periods = my_BlockPeriodMap(origin_block ,EffectiveObs,BlockLength);
                               % assign time index, truncating the number
                               % of observations if the bootstrap sample
                               % would otherwise be too long
                               II(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), target_msa_idx(1:LL))=periods(1: min(length(periods), EffectiveObs-ObsCount) )' * ones(1, LL);
                               % assign MSA index
                               JJ(ObsCount + 1 : min(ObsCount+length(periods), EffectiveObs), target_msa_idx(1:LL))=ones(min(length(periods), EffectiveObs-ObsCount), 1) * origin_msa_idx;
                               % keep track of length of bootstrapped sample 
                               ObsCount = ObsCount + length(periods);
                            end
                        end
                       
                        % transform to linear index on full E matrix
                        IDX = sub2ind([EffectiveObs, N], II, JJ);
                    otherwise
                        for gg=1:NoGroups
                            % randomize index on old matrix
                            aux_IDX = ceil(rand( EffectiveObs,sum(IndexM(:,gg))) * EffectiveObs*sum(IndexM(:,gg)) );
                            % index is on elements within group gg
                            [II,JJ]=ind2sub( [EffectiveObs, sum(IndexM(:,gg))], aux_IDX );
                            % transform to index on [1, ..., N]
                            JJnew=reshape(FindIndexM{gg}(1,JJ), EffectiveObs, sum(IndexM(:,gg)));
                            % transform to linear index on full E matrix
                            IDX(:,1==IndexM(:,gg)) = sub2ind([EffectiveObs, N], II, JJnew);
        %                     setdiff(unique(FindIndexM{gg}(:))',unique(reshape(IDX(:,1==IndexM(:,gg)),1,[])))
                        end
                end
            else
                IDX=( 1:(T-k-MissingObs))' * ones(1,N);
            end
        end
    %     Res_Lvl=reshape( ...
    %         reshape(Y,Params.p,(Params.T-Params.MissingObs-Params.k)*Params.N) ...
    %                 -A_IV'*reshape(X(1:Params.p*Params.k,:,:), ...
    %         Params.p*Params.k,(Params.T-Params.MissingObs-Params.k)*Params.N), [Params.p, Params.T-Params.MissingObs-Params.k, Params.N]);
        Y_b=NaN(p,(T-k-MissingObs),N);
        Z_b=NaN(pz,(T-k-MissingObs),N);
        X_b=NaN(px,(T-k-MissingObs),N);
        IV_Y_b=zeros( p_IV ,T-MissingObs-NoIVLags-1,N);    
%         dX_b=NaN(px,(T-k-MissingObs)-1,N);

        pY_b=NaN(pp,(T-k-MissingObs),N);
        pX_b=NaN(ppx,(T-k-MissingObs),N,pp);
        IV_Yp_b=zeros( p_IVp ,T-MissingObs-NoIVLags-1,N);    
%         dXp_b=NaN(ppx,(T-k-MissingObs)-1,N,pp);
        
        
        % Initialize: For the first k observations, the data is just data...
    %     Y_b(:, 1:k, :) = DATA.Y(:, 1:k, :);
        X_b(:, 1, :)  = DATA.X(:, 1, :);
%         dX_b(:, 1:k, :) = DATA.dX(:, 1:k, :);
        pX_b(:, 1, :,:)  = DATA.pX(:, 1, :,:);
%         dXp_b(:, 1:(1+k)*pp, :,:) = DATA.dXp(:, 1:(1+k)*pp, :,:);
    %     dY_b(:, 1:k, :) = DATA.dY(:, 1:k, :);
    %     IV_Y_b(:, 1:NoIVLags, :) = DATA.IV_Y(:, 1:NoIVLags , :);    
        ErrorTerm=zeros(p+pz,N);
        pErrorTerm=zeros(pp,N);
%         ErrorTerm=zeros(p+pz,N);    
        if DO_YEAR_FE
            
            switch upper(NON_STANDARD_H0)
                case 'IID'
                    for t_idx = 1 : EffectiveObs
            %             Y_b(:,t_idx,:) = reshape( A_IV' * squeeze(X_b(1:px,t_idx,:)) , [p,1,N]) ...
            %                 + reshape(DATA.MSA_FE_res(:, :), [p, 1, N]) ...
            %                 + repmat(DATA.Year_FE_res(:, t_idx), [1, 1, N]) ...
            %                 + DATA.Res_Lvl(:, t_idx, :) .* repmat(reshape(SIGN_FLIP(t_idx,:), [1, 1, N]), [p, 1, N]) ;
                        if DO_WILD
                            ErrorTerm=( squeeze(E_iid(:,BlocksIntoObs(:,t_idx),:)) .* repmat(SIGN_FLIP(BlocksIntoObs(:,t_idx),:), [p+pz, 1]) );
                            pErrorTerm=( squeeze(pE_iid(:,BlocksIntoObs(:,t_idx),:)) .* repmat(SIGN_FLIP(BlocksIntoObs(:,t_idx),:), [pp, 1]) );
                        else
                            if DO_TEST
                                ErrorTerm=squeeze( E_iid(:,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) );
                                pErrorTerm=squeeze( pE_iid(:,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) );
                            else
                                ErrorTerm=squeeze( E_iid(:,IDX(t_idx,:)) );
                                pErrorTerm=squeeze( pE_iid(:,IDX(t_idx,:)) );
                                [~,bbE,ccE]=ind2sub( size(E_iid), find(E_iid==ErrorTerm(1,5)) );
                                [bbI,ccI]=ind2sub( [T-MissingObs-k, N], IDX(t_idx,5) );
                %                 [~,bb2,cc2]=ind2sub( size(E_iid), find(E_iid==ErrorTerm_iid(2,5)) )
                                if ~isequal([bbI,ccI] , [bbE,ccE])  %this works
                                    error('indexing doesn''t work');
                                end
                            end
                        end


                        run bootstrap_inner_loop_split_Periphery
                    end
                case 'CONS'
                    for t_idx = 1 : T-MissingObs-k
            %             Y_b(:,t_idx,:) = reshape( A_IV' * squeeze(X_b(1:px,t_idx,:)) , [p,1,N]) ...
            %                 + reshape(DATA.MSA_FE_res(:, :), [p, 1, N]) ...
            %                 + repmat(DATA.Year_FE_res(:, t_idx), [1, 1, N]) ...
            %                 + DATA.Res_Lvl(:, t_idx, :) .* repmat(reshape(SIGN_FLIP(t_idx,:), [1, 1, N]), [p, 1, N]) ;
                        if DO_WILD
%                             ErrorTerm=( squeeze(E_iid(:,t_idx,:)) .* repmat(SIGN_FLIP(t_idx,:), [p+pz, 1]) );
                            ErrorTerm=( squeeze(E_cons(:,t_idx,:)) .* repmat(SIGN_FLIP(t_idx,:), [p+pz, 1]) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                            pErrorTerm=( squeeze(pE_cons(:,t_idx,:)) .* repmat(SIGN_FLIP(t_idx,:), [pp, 1]) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                        else
                            if DO_TEST
                                ErrorTerm=squeeze( E_cons(:,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                                pErrorTerm=squeeze( pE_cons(:,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                            else
                                ErrorTerm=squeeze( E_cons(:,IDX(t_idx,:)) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);                                    
                                pErrorTerm=squeeze( pE_cons(:,IDX(t_idx,:)) ) ...
                                        /(speye(N)-spdiags(IndexM*GroupedRHO_cons(p_idx,:)',0,N,N)*DistanceMatrix_combined);                                    
                            end
                        end


                        run bootstrap_inner_loop_split_Periphery
                    end
                otherwise
                    for t_idx = 1 : T-MissingObs-k
            %             Y_b(:,t_idx,:) = reshape( A_IV' * squeeze(X_b(1:px,t_idx,:)) , [p,1,N]) ...
            %                 + reshape(DATA.MSA_FE_res(:, :), [p, 1, N]) ...
            %                 + repmat(DATA.Year_FE_res(:, t_idx), [1, 1, N]) ...
            %                 + DATA.Res_Lvl(:, t_idx, :) .* repmat(reshape(SIGN_FLIP(t_idx,:), [1, 1, N]), [p, 1, N]) ;
                        if DO_WILD
                            for p_idx=1:(p+pz+pp)
                                for gg=1:NoGroups
                                    if p_idx <= p+pz
                                        ErrorTerm(p_idx, 1==IndexM(:,gg) )=( squeeze(E_opt(p_idx,t_idx,1==IndexM(:,gg)))' .* SIGN_FLIP(t_idx,1==IndexM(:,gg)) )...
                                            /(speye(N)-spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                                    else
                                        pErrorTerm(p_idx-p-pz, 1==IndexM(:,gg) )=( squeeze(pE_opt(p_idx-p-pz,t_idx,1==IndexM(:,gg)))' .* SIGN_FLIP(t_idx,1==IndexM(:,gg)) )...
                                            /(speye(N)-spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                                    end
                                end
                            end
                        else
                            if DO_TEST
                                for p_idx=1:(p+pz+pp)
                                    if p_idx <= p+pz
                                        ErrorTerm(p_idx,:)=squeeze( E_opt(p_idx ,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) )...
                                            /(speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N));
    %                                     norm( full((speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N))/(speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)))-eye(N) )
    %                                     transformed_error(p_idx,:,:)=squeeze(original_error(p_idx,:,:))*(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p_idx,:)',0,N,N)) ;
                                        display(corr(ErrorTerm(p_idx,:)',squeeze(E_iid(p_idx,t_idx,:))));
                                    else
                                        pErrorTerm(p_idx-p-pz,:)=squeeze( pE_opt(p_idx-p-pz ,sub2ind([T-MissingObs-k,N], reshape(IDX(t_idx,:),N,1)', 1:N  )) )...
                                            /(speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N));
    %                                     norm( full((speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N))/(speye(N)-sparse(DistanceMatrix_combined')*spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)))-eye(N) )
    %                                     transformed_error(p_idx,:,:)=squeeze(original_error(p_idx,:,:))*(speye(N)-sparse(W')*spdiags(IndexM*GroupedRHO(p_idx,:)',0,N,N)) ;
                                        display(corr(pErrorTerm(p_idx-p-pz,:)',squeeze(pE_iid(p_idx-p-pz,t_idx,:))));
                                    end
                                end
                            else
                                ErrorTerm=squeeze( E_iid(:,IDX(t_idx,:)) );
                                pErrorTerm=squeeze( pE_iid(:,IDX(t_idx,:)) );
                                [~,bbE,ccE]=ind2sub( size(E_iid), find(E_iid==ErrorTerm(1,5)) ); % pick first variable, 5th cross-sectional observation
                                [bbI,ccI]=ind2sub( [T-MissingObs-k, N], IDX(t_idx,5) );
                                if ~isequal([bbI,ccI] , [bbE,ccE])  %this works
                                    error('indexing doesn''t work');
                                end
                                for p_idx=1:(p+pz+pp)
                                    if p_idx <= p+pz
                                        ErrorTerm(p_idx,:)=squeeze( E_opt(p_idx,IDX(t_idx,:)) ) ...
                                            /(speye(N)-spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                                    else
                                        pErrorTerm(p_idx-p-pz,:)=squeeze( pE_opt(p_idx-p-pz,IDX(t_idx,:)) ) ...
                                            /(speye(N)-spdiags(IndexM*GroupedRHO_opt(p_idx,:)',0,N,N)*DistanceMatrix_combined);
                                    end
                                end
                            end
                        end


                        run bootstrap_inner_loop_split_Periphery
                    end
            end
        else
            error('bootstrap currently not set up for this case');
        end

        Yall_b=[DATA.Yall(:,1:NoIVLags-k,:), Y_b];
        pYall_b=[DATA.pYall(:,1:NoIVLags-k,:), pY_b];

        % Call the estimation function
        [~, K_IV_combined_b, IV_IRF_opt_b, IV_IRF_iid_b, IV_IRF_cons_b, IV_IRF_alt_b, logLL_b(:,b_idx), E_opt_b, ~, ~, ~, Ap_IV_b, pE_opt_b, ~, ~, ~, RHOS_b, DATA_b, DIST_WGT_b] ...
            = my_Spatial_Error_split_Periphery_jae(Y_b, X_b, pY_b, pX_b, Z_b, Yall_b,pYall_b, DistanceMatrix1, DistanceMatrix2, Params, IndexM);
    

        for gg=1:NoGroups
            MaxEig_b(gg,b_idx)=max(abs(eig(K_IV_combined_b{gg})));
        end

        if DO_TEST
            for gg=1:NoGroups
                if norm(K_IV_combined_b{gg}-K_IV_combined{gg}) > 1e-6 || norm(RHOS_b-RHOS) > 1e-4
                    error('check if bootstrap produces the right sample');
                end
            end
        end
        
        BootstrapData.RHOS_opt(:,:,b_idx)=RHOS_b(:,1:NoGroups);
        BootstrapData.RHOS_cons(:,:,b_idx)=RHOS_b(1,NoGroups+(1:NoGroups) );
        BootstrapData.DIST_WGT(:,:, b_idx) = DIST_WGT_b;

    %     Opt_RHOS_b(:,b_idx)=RHOS_b(:,1);
    %     Cons_RHOS_b(1,b_idx)=RHOS_b(1,2);
        for gg=1:NoGroups
%             for gg=1:NoGroups
            BootstrapData.K_IV_combined(:,:,gg,b_idx)=K_IV_combined_b{gg}; 
            for kk = 1:k
                BootstrapData.A_IV( (kk-1)*p + (1:p), :,gg,b_idx)=K_IV_combined_b{gg}(1:p,(kk-1)*(pp+p)+(1:p))'; 
            end
            BootstrapData.Ap_IV(:,:,gg,b_idx)=cell2mat(Ap_IV_b(gg,:));
%             end

            BootstrapData.IV_IRF_opt(:,:,gg,b_idx)=IV_IRF_opt_b{gg,1};
            BootstrapData.IV_IRF_opt_noDir(:,:,gg,b_idx)=IV_IRF_opt_b{gg,2};
            BootstrapData.IV_IRF_opt_noInd(:,:,gg,b_idx)=IV_IRF_opt_b{gg,3};
            BootstrapData.IV_IRF_opt_Neither(:,:,gg,b_idx)=IV_IRF_opt_b{gg,4};
            BootstrapData.IV_IRF_iid(:,:,gg,b_idx)=IV_IRF_iid_b{gg};
            BootstrapData.V_opt(:,:,gg,b_idx)=DATA_b.V_opt{gg}(1:p,1:p);
            BootstrapData.V_opt(:,:,gg,b_idx)=DATA_b.V_opt{gg}(1:p,1:p);

            BootstrapData.pV_opt(:,:,gg,b_idx)=DATA_b.pV_opt{gg};
            BootstrapData.pV_opt(:,:,gg,b_idx)=DATA_b.pV_opt{gg};
        end

        if BIAS_CORRECTION~=1
            
            for gg=1:NoGroups
                BootstrapData.IV_F_opt(1:pz,gg,b_idx)=my_first_stage_F_jae(E_opt_b(:,:,1==IndexM(:,gg)), T-k-MissingObs, Params, IV_IRF_opt_b{gg,1}(1:p,:),DATA_b.V_opt{gg});

                % IRF standardization
                % Potentially adjust the levels impact IRF, rather than the
                % one relative to sign
                % Normalize sign on the entire employment response, not the
                % employment to population ratio response
                if DO_MR_Standardization
                    for ii=1:4
                        IV_IRF_opt_b{gg,ii}=IV_IRF_opt_b{gg,ii}*diag(1./diag(IV_IRF_opt_b{gg,ii}(1:pv1,1:pv1)+POP_ADJUST*IV_IRF_opt_b{gg,ii}(POP_IDX,1:pv1)));
                    end
                    IV_IRF_iid_b{gg}=IV_IRF_iid_b{gg}*diag(1./diag(IV_IRF_iid_b{gg}(1:pv1,1:pv1)+POP_ADJUST*IV_IRF_iid_b{gg}(POP_IDX,1:pv1)));
                else
                    for ii=1:4
                        IV_IRF_opt_b{gg,ii}=IV_IRF_opt_b{gg,ii}*diag(sign(diag(IV_IRF_opt_b{gg,ii}(1:pv1,1:pv1)+POP_ADJUST*IV_IRF_opt_b{gg,ii}(POP_IDX,1:pv1))));
                    end
                    IV_IRF_iid_b{gg}=IV_IRF_iid_b{gg}*diag(sign(diag(IV_IRF_iid_b{gg}(1:pv1,1:pv1)+POP_ADJUST*IV_IRF_iid_b{gg}(POP_IDX,1:pv1))));
                end
                
                % Compute IRFs here
                IRF_pooled_b{1,gg}(:,:,1,b_idx)=IV_IRF_opt_b{gg,1}; 
                IRF_alt_b{1,gg}(:,:,1,b_idx)=IV_IRF_opt_b{gg,2}; 
                IRF_alt_b{2,gg}(:,:,1,b_idx)=IV_IRF_opt_b{gg,3}; 
                IRF_alt_b{3,gg}(:,:,1,b_idx)=IV_IRF_opt_b{gg,4}; 
                IRF_pooled_b{2,gg}(:,:,1,b_idx)=IV_IRF_iid_b{gg}; 
                for hh=1:h
                    IRF_pooled_b{1,gg}(:,:,hh+1, b_idx)=eye(p+pp,(p+pp)*k)*K_IV_combined_b{gg}^hh*[IV_IRF_opt_b{gg,1}; zeros((p+pp)*(k-1),pv1)];
                    IRF_pooled_b{2,gg}(:,:,hh+1, b_idx)=eye(p+pp,(p+pp)*k)*K_IV_combined_b{gg}^hh*[IV_IRF_iid_b{gg}; zeros((p+pp)*(k-1),pv1)];

                    IRF_alt_b{1,gg}(:,:,hh+1, b_idx)=eye(p+pp,(p+pp)*k)*K_IV_combined_b{gg}^hh*[IV_IRF_opt_b{gg,2}; zeros((p+pp)*(k-1),pv1)];
                    IRF_alt_b{2,gg}(:,:,hh+1, b_idx)=eye(p+pp,(p+pp)*k)*K_IV_combined_b{gg}^hh*[IV_IRF_opt_b{gg,3}; zeros((p+pp)*(k-1),pv1)];
                    IRF_alt_b{3,gg}(:,:,hh+1, b_idx)=eye(p+pp,(p+pp)*k)*K_IV_combined_b{gg}^hh*[IV_IRF_opt_b{gg,4}; zeros((p+pp)*(k-1),pv1)];
                end
            end
         end

        fprintf('.');
        pause(1);
        if mod(b_idx,10)==0
            fprintf([num2str(b_idx,'%1.0f ') '\n']);
        end
    end
        if BIAS_CORRECTION==1
            A_IV_original=cell(NoGroups,1);
            Ap_IV_original=cell(NoGroups,pp);
            K_IV_original=cell(NoGroups,1);

            for gg=1:NoGroups
                A_IV_original{gg}=A_IV{gg};
                Ap_IV_original(gg,:)=Ap_IV(gg,:);
                K_IV_original{gg}=K_IV_combined{gg};
                RHOS_original=RHOS;

                A_IV_BS=mean(BootstrapData.A_IV(:,:,gg,:),4);
                Ap_IV_BS=mean(BootstrapData.Ap_IV(:,:,gg,:),4);
                if max(abs(eig(K_IV_original{gg})))<1
                    DownBias_A_IV=A_IV_BS-A_IV_original{gg};
                    DownBias_Ap_IV=Ap_IV_BS-cell2mat(Ap_IV_original(gg,:));
                    
                    f_K_IV_BS = @(BM) my_companion_form(A_IV_original{gg} - BM * DownBias_A_IV, cell2mat(Ap_IV_original(gg,:)) - BM * DownBias_Ap_IV, Params);
                    
                    BiasMultiple=1;
                    while max(abs(eig( f_K_IV_BS(BiasMultiple) )))>=1
                        BiasMultiple=BiasMultiple-0.01;
                    end
                    [K_IV_combined{gg}, A_combined{gg}, LHS{gg}] = my_companion_form(A_IV_original{gg} - BiasMultiple * DownBias_A_IV, cell2mat(Ap_IV_original(gg,:)) - BiasMultiple * DownBias_Ap_IV, Params);
                    
                    A_IV{gg} = A_IV_original{gg} -BiasMultiple*DownBias_A_IV;
                    for p_idx=1:pp
                        Ap_IV{gg,pp}=Ap_IV_original{gg,pp} - BiasMultiple * DownBias_Ap_IV(:,p_idx) ;
                    end
                end

                DownBias_RHOS_opt=mean(BootstrapData.RHOS_opt,3)-RHOS(:,1:NoGroups);
                BiasMultiple=1;
                while max(max(abs((RHOS(:,1:NoGroups)-DownBias_RHOS_opt*BiasMultiple))))>=1
                    BiasMultiple=BiasMultiple-0.01;
                end
                RHOS(:,1:NoGroups)=RHOS_original(:,1:NoGroups)-DownBias_RHOS_opt*BiasMultiple;
                
                DownBias_RHOS_cons=mean(BootstrapData.RHOS_cons,3)-RHOS(1,NoGroups+(1:NoGroups) );
                BiasMultiple=1;
                while abs((RHOS(1,NoGroups+(1:NoGroups) )-DownBias_RHOS_cons*BiasMultiple))>=1
                    BiasMultiple=BiasMultiple-0.01;
                end
                RHOS(:,NoGroups+(1:NoGroups) )=RHOS_original(:,NoGroups+(1:NoGroups) )-ones(p+pz+pp,1)*DownBias_RHOS_cons*BiasMultiple;
            end
            B_REPS_current=B_REPS;
        end
end

%% Compute structural implication
my_VAR_Plots;
my_Variance_Decomposition_jae;
close all;
%% Display spatial correlation and fit

varIDX=[1:p, p+pz+(1:pp), p+1:p+pz]; % Core variables first, then periphery, then instruments
for gg=1:NoGroups
    
    NumberTables = [RHOS_original(:,gg) DownBias_RHOS_opt(:,gg) prctile( squeeze(BootstrapData.RHOS_opt(:,gg,:))', [5 16 50 84 95])'];
    
    TableRHOS=[ char([VarLabels(1:p); PerLabels(:); VarLabels(p+1:p+pz)] ) ...
            num2str(NumberTables(varIDX,:) , '& %1.2f' ) repmat(' \\', p+pz+pp,1)];

    disp(cellstr(TableRHOS));
    CellTable=cellstr(TableRHOS);

    my_fid=fopen(['../Tables/RHOs_' Suffix num2str(gg, '%1.0f') '.tex'],'w');
    fprintf(my_fid, '%s\n',CellTable{:});
    fclose(my_fid);
end

TableLL=[ char({'Same $\rho$ vs. varying $\rho$s'; 'No spatial correlation vs. varying $\rho$s'}) ...
       repmat(' &', 2,1)  num2str([-2*(logLL(2:3)'-logLL(1)) prctile( -2*(logLL_b(2:3,:)-[1;1]*logLL_b(1,:))', [1 5 50 95 99])']/N, '& %1.3f' ) ...
         repmat(' \\', 2,1)];

disp(cellstr(TableLL));
CellTable=cellstr(TableLL);

my_fid=fopen(['../Tables/LL_' Suffix '.tex'],'w');
fprintf(my_fid, '%s\n',CellTable{:});
fclose(my_fid);


for gg=1:NoGroups
    try
        CellTable=cellstr([char(VarLabels(1:pv1)) ...
                num2str([Original_IV_F{gg} prctile(squeeze(BootstrapData.IV_F_opt(1:pv1,gg,:)), CI_lvls,2)], '& %1.1f')    repmat(' \\', pv1,1)]);

        disp(CellTable);
        my_fid=fopen(['../Tables/F_stat_BStrap_' Suffix num2str(gg, '%1.0f')  '.tex'],'w');
        fprintf(my_fid, '%s\n',CellTable{:});
        fclose(my_fid);
    catch
        fprintf(2, 'F-stat bootstrap not working in this case');
        fprintf('\n');
    end
end

%%
save(['results' Suffix])

