clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication and Extension of Acemoglu Johnson Robinson Yared 2008 AER
% and Cervellati Jung Sunde Vischer 2014 AER
%
% Calculate Arellano-Bond GMM and Iterated GMM for Dynamic Panel
%
% Produces the results in Tables 3 and 4 of Hansen and Lee (forthcoming at Econometrica)
% "Inference for Iterated GMM Under Misspecification"
%
% For errors and/or comments please contact Seojeong (Jay) Lee at
% jay.lee@unsw.edu.au
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need the following functions in the same folder
%
% arellano_bond.m        calculates the data matrices for the difference GMM
% iterated_gmm_cluster.m calculates the iterated GMM, se, and the J test
%
% Need the following data files in the same folder
%
% fiveyearpanel.mat     five year data
% tenyearpanel.mat      ten year data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Choose a Table to be replicated
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Acemoglu Johnson Robinson Yared 2008
% tc = 24: AJRY Table 2 Column 4
% tc = 28: AJRY Table 2 Column 8

% Cervellati Jung Sunde Vischer 2014
% tc = 41: CJSV Table 4 Column 1 Panel B
% tc = 42: CJSV Table 4 Column 2 Panel B
% tc = 43: CJSV Table 4 Column 3 Panel B 
% tc = 44: CJSV Table 4 Column 4 Panel B
% tc = 45: CJSV Table 4 Column 5 Panel B
% tc = 46: CJSV Table 4 Column 6 Panel B

tc = 24;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct Data Matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if tc==24 || tc==41 || tc==43 || tc==45
    load fiveyearpanel.mat
else
    load tenyearpanel.mat
end

d = fhpolrigaug;
y = lrgdpch;

if tc==24 || tc==28
    x = [];
else
    d = d(colony==1);
    y = y(colony==1);
    year_numeric = year_numeric(colony==1);
    country = country(colony==1);
    if tc==41 || tc==42
        x = xc_y;
    elseif tc==43 || tc==44
        x = y_indep_fl;
    elseif tc==45 || tc==46
        x = y_power;
    end
    x = x(colony==1);
end

[Dd,DX,Z,mem,n,G,ng] = arellano_bond(d,y,x,year_numeric,country);

k = length(DX(1,:));
l = max(ng)*(max(ng)+1)/2+k-1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GMM Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% One-step GMM

W0 = zeros(l,l);
W0i = zeros(l,l,G);

for i = 1:G
    Zi = Z(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)),:);

    h0 = 2*ones(ng(i),1);
    h1 = -1*ones(ng(i)-1,1);
    Hm = diag(h0)+diag(h1,1)+diag(h1,-1);

    W0i(:,:,i) = Zi'*Hm*Zi;
    W0 = W0 + W0i(:,:,i);
end
W0 = W0/n;

% One-step GMM in Columns I, III of Table 3 and Columns I, III, V in Table 4
b1 = ((DX'*Z)/W0*(Z'*DX))\(DX'*Z)/W0*(Z'*Dd);

Q = -(Z'*DX)/n;
e1 = Dd - DX*b1;
mu1 = (Z'*e1)/n;

Om01 = zeros(l,l);
Om1 = zeros(k,k);
for i = 1:G
    Zi = Z(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)),:);
    e1i = e1(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)));
    DXi = DX(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)),:);
    
    Om01 = Om01 + (Zi'*e1i)*(Zi'*e1i)';
    
    psi1 = Q'/W0*Zi'*e1i-DXi'*Zi/W0*mu1-Q'/W0*W0i(:,:,i)/W0*mu1;
    Om1 = Om1 + psi1*psi1';
end

Om01 = Om01/n;
Om1 = Om1/n;

H0 = Q'/W0*Q;

% Conventional Heteroskedasticity-and-Cluster-robust SE
V1 = H0\(Q'/W0*Om01/W0*Q)/H0;
s1 = sqrt(diag(V1/n));

% The new fully robust SE
V1r = H0\Om1/H0;
s1r = sqrt(diag(V1r/n));


%%% Two-step GMM
b2 = ((DX'*Z)/Om01*(Z'*DX))\(DX'*Z)/Om01*(Z'*Dd);

e2 = Dd - DX*b2;
mu2 = (Z'*e2)/n;

Om2 = zeros(l,l);
for i = 1:G
    Zi = Z(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)),:);
    e2i = e2(sum(ng(1:i))-ng(i)+1:sum(ng(1:i)));
    
    Om2 = Om2 + (Zi'*e2i)*(Zi'*e2i)';
end

Om2 = Om2/n;

% Conventional Heteroskedasticity-and-Cluster-robust SE
V2 = inv(Q'/Om2*Q);
s2 = sqrt(diag(V2/n));

% Overidentifying restrictions J test
J2 = n*mu2'/Om2*mu2;
pv2 = chi2cdf(J2,l-k,'upper');


%%% Iterated GMM
% b: Iterated GMM in Columns II, IV of Table 3 and Columns II, IV, VI in Table 4
% s: Misspecification-Robust SE
% s0: Arellano-Bond SE
% sw: Windmeijer-corrected SE
% pv: the P-value of the J test
% iter: Numver of Iterations 
[b,sr,Vr,sw,Vw,s0,V0,J,pv,iter] = iterated_gmm_cluster(Dd,DX,Z,mem);



% Cumulative Income Effect and SE's 
if tc==24 || tc==28
    
    % One-step GMM
    b1_cie = b1(2)/(1-b1(1));
    Rb1 = [b1(2)/(1-b1(1))^2; 1/(1-b1(1)); zeros(k-2,1)];

    V1_cie = Rb1'*V1*Rb1;
    s1_cie = sqrt(diag(V1_cie/n));

    V1r_cie = Rb1'*V1r*Rb1;
    s1r_cie = sqrt(diag(V1r_cie/n));

    % Iterated GMM
    b_cie = b(2)/(1-b(1));
    RbI = [b(2)/(1-b(1))^2; 1/(1-b(1)); zeros(k-2,1)];

    V_cie = RbI'*Vr*RbI;
    s_cie = sqrt(diag(V_cie/n));

    V0_cie = RbI'*V0*RbI;
    s0_cie = sqrt(diag(V0_cie/n));
end

