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)
% IV estimation with differenced data

if sum(strcmp('WhichChol',fieldnames(Params)))==0
    Params.WhichChol=1;
end
if sum(strcmp('CombineDistances',fieldnames(Params)))==0
    Params.CombineDistances=0;
end
NoGroups= size(IndexM,2);

dY=diff(Y,1,2);

dX=diff(X(1:Params.k*Params.p,:,:),1,2);
dYp=diff(pY,1,2);
dXp=diff(pX(1:(Params.k + Params.p* (1+Params.k)),:,:,:),1,2);
p_IV=Params.NoIVLags*Params.p;
p_IVp=Params.NoIVLags + Params.p* (1+Params.k);

%% (1) Construct data for IV

dY(:,1:(Params.NoIVLags-Params.k),:)=[];
dX(:,1:(Params.NoIVLags-Params.k),:)=[];

dYp(:,1:(Params.NoIVLags-Params.k),:)=[];
dXp(:,1:(Params.NoIVLags-Params.k),:,:)=[];

IV_Y=zeros( p_IV ,Params.T-Params.MissingObs-Params.NoIVLags-1,Params.N);
IV_Yp=zeros( p_IVp ,Params.T-Params.MissingObs-Params.NoIVLags-1,Params.N, Params.pp);
IV_Yp(Params.NoIVLags +1:end,:,:,:)=dXp(Params.k +1:end,:,:,:); % exogenous block of periphery regressors
for l_idx=2:(Params.NoIVLags+1)
    for p_idx=1:max(Params.p,Params.pp)
        if p_idx <= Params.p
            IV_Y((l_idx-2)*Params.p+p_idx,:,:)=Yall(p_idx,(Params.NoIVLags+2:Params.T-Params.MissingObs)-l_idx,:);
        end;
        if p_idx <= Params.pp
            IV_Yp((l_idx-2)+1,:,:,p_idx)=pYall(p_idx,(Params.NoIVLags+2:Params.T-Params.MissingObs)-l_idx,:);
        end;
    end;
end;

%% (2) Take out year fixed effects 
DATA.Year_FE_dY=zeros(Params.p, Params.T-Params.MissingObs-Params.NoIVLags-1,NoGroups);
DATA.Year_FE_dYp=zeros(Params.pp, Params.T-Params.MissingObs-Params.NoIVLags-1,NoGroups);
for t_idx=1:Params.T-Params.MissingObs-Params.NoIVLags-1
    for gg=1:NoGroups
        DATA.Year_FE_dY(:,t_idx, gg)=mean(squeeze(dY(:,t_idx,1==IndexM(:,gg))),2);
        dY(:,t_idx,1==IndexM(:,gg))=dY(:,t_idx,1==IndexM(:,gg))-repmat( DATA.Year_FE_dY(:,t_idx, gg), [1, 1, sum(IndexM(:,gg))]);
        dX(:,t_idx,1==IndexM(:,gg))=dX(:,t_idx,1==IndexM(:,gg))-repmat( mean(squeeze(dX(:,t_idx,1==IndexM(:,gg))),2) , [1, 1, sum(IndexM(:,gg))]);
        IV_Y(:,t_idx,1==IndexM(:,gg))=IV_Y(:,t_idx,1==IndexM(:,gg))-repmat( mean(squeeze(IV_Y(:,t_idx,1==IndexM(:,gg))),2) , [1, 1, sum(IndexM(:,gg))]);

        DATA.Year_FE_dYp(:,t_idx, gg)=mean(reshape(dYp(:,t_idx,1==IndexM(:,gg)), [Params.pp, sum(IndexM(:,gg))]),2);
        dYp(:,t_idx,1==IndexM(:,gg))=dYp(:,t_idx,1==IndexM(:,gg))-repmat( DATA.Year_FE_dYp(:,t_idx, gg), [1, 1, sum(IndexM(:,gg))]);
        dXp(:,t_idx,1==IndexM(:,gg),:)=dXp(:,t_idx,1==IndexM(:,gg),:)-repmat( mean( dXp(:,t_idx,1==IndexM(:,gg),:) ,3) , [1, 1, sum(IndexM(:,gg)),1]);
        IV_Yp(:,t_idx,1==IndexM(:,gg),:)=IV_Yp(:,t_idx,1==IndexM(:,gg),:)-repmat( mean( IV_Yp(:,t_idx,1==IndexM(:,gg),:),3) , [1, 1, sum(IndexM(:,gg)),1]);
    end;
end;

%% (3) IV estimation

A_IV=cell(NoGroups,1);
LHS=cell(NoGroups,1);
Ap_IV=cell(NoGroups,Params.pp);
A_IV_combined=cell(NoGroups,1);
K_IV=cell(NoGroups,1);
K_IV_combined=cell(NoGroups,1);
for gg=1:NoGroups
    gamma_IV=reshape(IV_Y(:,:,1==IndexM(:,gg)),p_IV,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))'...
        \reshape(dX(:,:,1==IndexM(:,gg)),Params.p*Params.k,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))';
    dX_IV=reshape(IV_Y(:,:,1==IndexM(:,gg)),p_IV,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))'*gamma_IV;
    A_IV{gg}=dX_IV\reshape(dY(:,:,1==IndexM(:,gg)),Params.p,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))';
    K_IV{gg}=[A_IV{gg}(1:Params.p*Params.k,1:Params.p)'; eye( (Params.k-1)*Params.p , Params.k*Params.p)] ;

    for p_idx=1:Params.pp
        gamma_IVp=reshape(IV_Yp(:,:,1==IndexM(:,gg),p_idx),p_IVp,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))'...
            \reshape(dXp(:,:,1==IndexM(:,gg),p_idx),Params.p*(1+Params.k)+Params.k,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))';
        dXp_IV=reshape(IV_Yp(:,:,1==IndexM(:,gg),p_idx),p_IVp,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))'*gamma_IVp;
        Ap_IV{gg,p_idx}=dXp_IV\reshape(dYp(p_idx,:,1==IndexM(:,gg)),1,(Params.T-Params.MissingObs-Params.NoIVLags-1)* sum(IndexM(:,gg)))';
    end;
    % note that A_IV_combined is already transposed
    [K_IV_combined{gg}, A_combined{gg}, LHS{gg}] = my_companion_form(A_IV{gg}, cell2mat(Ap_IV(gg,:)), Params);
    if Params.DO_VERBOSE
        fprintf('Eigenvalues converted to quarterly frequency are: \n');
        display(sort(abs(eig(K_IV_combined{gg})).^(1/4)));
    end;
end;
%% (4) Compute residuals 

% (4a) Residuals net of year fixed effects
Res_Lvl=nan(size(Y));
pRes_Lvl=nan(size(pY));
for gg=1:NoGroups
    Res_Lvl(:,:,1==IndexM(:,gg))=reshape( ...
        reshape(Y(:,:,1==IndexM(:,gg)),Params.p,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg)) ) ...
                -A_IV{gg}'*reshape(X(1:Params.p*Params.k,:,1==IndexM(:,gg)), ...
        Params.p*Params.k,(Params.T-Params.MissingObs-Params.k)* sum(IndexM(:,gg))), [Params.p, Params.T-Params.MissingObs-Params.k, sum(IndexM(:,gg)) ]);
    for p_idx = 1:Params.pp
        pRes_Lvl(p_idx,:,1==IndexM(:,gg))=reshape( ...
            reshape(pY(p_idx,:,1==IndexM(:,gg)),1,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg)) ) ...
                        -Ap_IV{gg,p_idx}'*reshape(pX(:,:,1==IndexM(:,gg),p_idx), Params.ppx,(Params.T-Params.MissingObs-Params.k)* sum(IndexM(:,gg))), ...
                    [1, Params.T-Params.MissingObs-Params.k, sum(IndexM(:,gg)) ]);
    end;
end;
DATA.Year_FE_res=zeros(Params.p, Params.T-Params.MissingObs-Params.k, NoGroups);
DATA.Year_FE_pres=zeros(Params.pp, Params.T-Params.MissingObs-Params.k, NoGroups);
DATA.Year_FE_Z=zeros(Params.pz, Params.T-Params.MissingObs-Params.k, NoGroups);
for t_idx=1:Params.T-Params.k-Params.MissingObs
    for gg=1:NoGroups
        DATA.Year_FE_res(:,t_idx, gg) = mean(squeeze(Res_Lvl(:,t_idx,1==IndexM(:,gg))),2);
        DATA.Year_FE_pres(:,t_idx, gg) = mean(reshape(pRes_Lvl(:,t_idx,1==IndexM(:,gg)), [Params.pp, 1, sum(IndexM(:,gg))]),3);
        Res_Lvl(:,t_idx,1==IndexM(:,gg))=Res_Lvl(:,t_idx,1==IndexM(:,gg)) - repmat( DATA.Year_FE_res(:,t_idx,gg) , [1, 1, sum(IndexM(:,gg))]);
        pRes_Lvl(:,t_idx,1==IndexM(:,gg))=pRes_Lvl(:,t_idx,1==IndexM(:,gg)) - repmat( DATA.Year_FE_pres(:,t_idx,gg) , [1, 1, sum(IndexM(:,gg))]);
        DATA.Year_FE_Z(:,t_idx, gg) = mean(reshape(Z(:,t_idx,1==IndexM(:,gg)), [Params.pz, sum(IndexM(:,gg))]),2);
        Z(:,t_idx,1==IndexM(:,gg))=Z(:,t_idx,1==IndexM(:,gg)) - repmat( DATA.Year_FE_Z(:,t_idx, gg) , [1, 1, sum(IndexM(:,gg))]);
    end;
end;

% (4b) Residuals net of MSA fixed effects
EXO=ones(Params.T-Params.k-Params.MissingObs,1);
PROJECTION=EXO*( (EXO'*EXO)\EXO' );
DATA.MSA_FE_res=zeros(Params.p, Params.N);
DATA.MSA_FE_pres=zeros(Params.pp, Params.N);
DATA.MSA_FE_Z=zeros(Params.pz, Params.N);
for n_idx=1:Params.N
    DATA.MSA_FE_res(:,n_idx)=mean(squeeze(Res_Lvl(:,:,n_idx)),2);
    Res_Lvl(:,:,n_idx)=Res_Lvl(:,:,n_idx)-repmat(reshape(DATA.MSA_FE_res(:,n_idx), [Params.p, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

    DATA.MSA_FE_Z(:,n_idx)=mean(squeeze(Z(:,:,n_idx)),2);
    Z(:,:,n_idx)=Z(:,:,n_idx) - repmat(reshape(DATA.MSA_FE_Z(:,n_idx), [Params.pz, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

    DATA.MSA_FE_pres(:,n_idx)=mean(squeeze(pRes_Lvl(:,:,n_idx)),2);
    pRes_Lvl(:,:,n_idx)=pRes_Lvl(:,:,n_idx)-repmat(reshape(DATA.MSA_FE_pres(:,n_idx), [Params.pp, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

end;

%% (5) Compute spatial lag

if Params.DO_VERBOSE
    my_opt_options=optimset('TolX', 1e-10, 'TolFun', 1e-10, 'Display', 'iter', 'MaxFunEval', 5e3);
else
    my_opt_options=optimset('TolX', 1e-10, 'TolFun', 1e-10, 'Display', 'off', 'MaxFunEval', 5e3);
end
DIST_WGT = ones(2,NoGroups);
if Params.CombineDistances==1
    % optimize over optimal combination of distance matrices and spatial correlation
    aux_sol =fminunc(@(x) -my_spatial_LL_split_Periphery_v2(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
        DistanceMatrix1/(1+exp(x(2))) + DistanceMatrix2/(1+exp(-x(2))),  ...
        0.999*ones(Params.p+Params.pp+Params.pz,NoGroups)*tanh(x(1)), ... spatial correlation
        [],IndexM ) /(Params.N*Params.T), [0;0] , my_opt_options);
    cons_sol = aux_sol(1);
    DIST_WGT(2,:) = 1/(1+exp(aux_sol(2)));

    aux_sol =fminunc(@(x) -my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
        my_weighted_Distance_matrix(DistanceMatrix1, DistanceMatrix2, 1./(1+exp(x(end,:))), IndexM, NoGroups), ... distance matrix combination
        0.999*tanh(x(1:end-1,:)), [],IndexM ) /(Params.N*Params.T), ... spatial correlation
        [cons_sol*ones(Params.p+Params.pz+Params.pp, NoGroups);  ones(1,NoGroups)*aux_sol(2)] , my_opt_options);

    sol=aux_sol(1:end-1,:);
    DIST_WGT(1,:) = 1./(1+exp(aux_sol(end,:)));
else
    cons_sol =fminunc(@(x) -my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
        DistanceMatrix1,  ...
        0.999*ones(Params.p+Params.pp+Params.pz,NoGroups)*tanh(x(1)), ... spatial correlation
        [],IndexM ) /(Params.N*Params.T), 0 , my_opt_options);

    sol =fminunc(@(x) -my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
        DistanceMatrix1, ... distance matrix combination
        0.999*tanh(x), [],IndexM ) /(Params.N*Params.T), ... spatial correlation
        cons_sol*ones(Params.p+Params.pz+Params.pp, NoGroups) , my_opt_options);
end;
%% (6) Compute output of interest
[logLL_opt, V_opt, E_opt, pV_opt, pE_opt]=my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
    my_weighted_Distance_matrix(DistanceMatrix1, DistanceMatrix2, DIST_WGT(1,:), IndexM, NoGroups), ...
    0.999*tanh(sol),[], IndexM);

[logLL_cons, V_cons, E_cons, pV_cons, pE_cons]=my_spatial_LL_split_Periphery_jae(Res_Lvl,pRes_Lvl,  Z, Params.T-Params.MissingObs-Params.k, Params.N, ...
    DistanceMatrix1*DIST_WGT(2,1) + DistanceMatrix2*(1-DIST_WGT(2,1)),  ...
    0.999*ones(Params.p+Params.pz+Params.pp,NoGroups)*tanh(cons_sol),[], IndexM ) ;
% no weights here, because the distance matrix is irrelevant
[logLL_iid, V_iid, E_iid, pV_iid, pE_iid]=my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, DistanceMatrix1,  zeros(Params.p+Params.pz+Params.pp,NoGroups),[],IndexM);

IV_IRF_opt=cell(NoGroups,4);
IV_IRF_cons=cell(NoGroups,1);
IV_IRF_iid=cell(NoGroups,1);
for gg=1:NoGroups
    if Params.pv1<Params.pz
        [opt_IRF_IV, opt_IRF_IV_noDir, opt_IRF_IV_noInd, opt_IRF_IV_Neither]=my_IRF_IV(V_opt{gg}(Params.p+(1:Params.pz),:),         V_opt{gg}(1:Params.p,1:Params.p), Params.pv1, -Params.WhichChol);
        IV_IRF_opt{gg,1}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV;
        IV_IRF_opt{gg,2}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_noDir;
        IV_IRF_opt{gg,3}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_noInd;
        IV_IRF_opt{gg,4}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_Neither;
        IV_IRF_cons{gg}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * my_IRF_IV(V_cons{gg}(Params.p+(1:Params.pz),:),         V_cons{gg}(1:Params.p,1:Params.p), Params.pv1, -Params.WhichChol);
        IV_IRF_iid{gg}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * my_IRF_IV(V_iid{gg}(Params.p+(1:Params.pz),:),         V_iid{gg}(1:Params.p,1:Params.p), Params.pv1, -Params.WhichChol);
    else
        [opt_IRF_IV, opt_IRF_IV_noDir, opt_IRF_IV_noInd, opt_IRF_IV_Neither]=my_IRF_IV(V_opt{gg}(Params.p+(1:Params.pz),1:Params.p),V_opt{gg}(1:Params.p,1:Params.p),[],  Params.WhichChol);
        IV_IRF_opt{gg,1}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV;
        IV_IRF_opt{gg,2}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_noDir;
        IV_IRF_opt{gg,3}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_noInd;
        IV_IRF_opt{gg,4}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * opt_IRF_IV_Neither;
        IV_IRF_cons{gg}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * my_IRF_IV(V_cons{gg}(Params.p+(1:Params.pz),1:Params.p),V_cons{gg}(1:Params.p,1:Params.p),[], Params.WhichChol);
        IV_IRF_iid{gg}=[eye(Params.p); -LHS{gg}( Params.p + (1:Params.pp), 1:Params.p )] * my_IRF_IV(V_iid{gg}(Params.p+(1:Params.pz),1:Params.p),V_iid{gg}(1:Params.p,1:Params.p),[], Params.WhichChol);
    end
end

if Params.DO_IRF==2
    error('Params.DO_IRF==2 is disabled');
    sol_w =fminsearch(@(x) -my_spatial_LL_v1(Res_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, DistanceMatrix2, 0.999*tanh(x) ) /(Params.N*Params.T), cons_sol*ones(Params.p+Params.pz,1) , my_opt_options);

    % % Mix two weighting matrices
    sol2 =fminsearch(@(x) -my_spatial_LL_v1(Res_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, DistanceMatrix1, DistanceMatrix2, ...
        0.999*tanh(x(1:Params.p+Params.pz)), ...
        0.999*tanh(x(Params.p+Params.pz+(1:Params.p+Params.pz) )) ) /(Params.N*Params.T),  [sol/2; sol_w/2] , my_opt_options);

    [logLL_alt, V_alt, E_alt]=my_spatial_LL_W2(Res_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, DistanceMatrix1, DistanceMatrix2, ...
        0.999*(1-exp(-sol2(1:Params.p+Params.pz)))./(1+exp(-sol2(1:Params.p+Params.pz))), ...
        0.999*(1-exp(-sol2(Params.p+Params.pz+(1:Params.p+Params.pz) )))./(1+exp(-sol2(Params.p+Params.pz+(1:Params.p+Params.pz) ))));
    IV_IRF_alt=my_IRF_IV(V_alt(Params.p+(1:Params.pz),1:Params.p),V_alt(1:Params.p,1:Params.p),[], 1);

    logLL=[logLL_opt, logLL_cons, logLL_iid, logLL_alt];

    RHOS=0.999*[tanh(sol), ones(Params.p+Params.pz,1)*tanh(cons_sol), ...
        tanh(sol2(1:Params.p+Params.pz)), tanh(sol2(Params.p+Params.pz+(1:Params.p+Params.pz) )) ];
else
    logLL=[logLL_opt, logLL_cons, logLL_iid ];

    RHOS=0.999*[tanh(sol), ones(Params.p+Params.pz+Params.pp,NoGroups)*tanh(cons_sol) ];

    IV_IRF_alt=[];
    E_alt=[];
    pE_alt=[];

end;

%% (7) Output data

DATA.Y=Y;
DATA.dY=dY;
DATA.pY=pY;
DATA.dYp=dYp;
DATA.Yall=Yall;
DATA.pYall=pYall;
DATA.X=X;

DATA.pX=pX;

DATA.IV_Y=IV_Y;
DATA.Z=Z;
DATA.Res_Lvl=Res_Lvl;
DATA.pRes_Lvl=pRes_Lvl;
DATA.V_opt=V_opt;
DATA.V_iid=V_iid;
DATA.V_cons=V_cons;
DATA.pV_opt=pV_opt;
DATA.pV_iid=pV_iid;
DATA.pV_cons=pV_cons;
DATA.A_IV=A_IV;
end