clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code produces Table 1 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sample size
n = 250;

% Number of Instruments
l = 4;

% Weak Instrument (strong R2=0.2 vs weak R2=0.02)
R2 = 0.2;
pi = (R2/(1-R2))^0.5*(1/l^0.5);

% Degree of Endogeneity
rho = 0.5;

% True value = Pseudo-true value
b = 1;

% Number of Monte Carlo repetition
r = 1000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

del = zeros(l,1);
for j = 1:l
    del(j) = (-1)^(j-1);
end
if mod(l,2) == 1
    del(l) = 0;
end

for a = 0:0.2:1

bhat = zeros(r,1);
s = zeros(r,1);
sw = zeros(r,1);
s0 = zeros(r,1);
J = zeros(r,1);
pv = zeros(r,1);
numit = zeros(r,1);

Size_t = zeros(3,1);
Size_tw = zeros(3,1);
Size_t0 = zeros(3,1);

for iter = 1:r

    z = randn(n,l);
    err = mvnrnd([0;0],[1 rho;rho 1],n);
    u = err(:,2);
    e = err(:,1);

    x = pi*sum(z,2) + u;
    y = x*b + a*(z*del) + e;
    
    xhat = z/(z'*z)*(z'*x);
    R2 = 1 - sum((x-xhat).^2)/sum((x-mean(x)).^2);
    
    [bhat(iter),s(iter),V,sw(iter),Vw,s0(iter),V0,J(iter),pv(iter),numit(iter)] = iterated_gmm_cluster(y,x,z,(1:n)');
    
    t = (bhat(iter)-b)/s(iter);
    tw = (bhat(iter)-b)/sw(iter);
    t0 = (bhat(iter)-b)/s0(iter);

    if isnan(bhat(iter))==0
        Size_t = Size_t + (abs(t)>norminv([0.995 0.975 0.95]',0,1));
        Size_tw = Size_tw + (abs(tw)>norminv([0.995 0.975 0.95]',0,1));
        Size_t0 = Size_t0 + (abs(t0)>norminv([0.995 0.975 0.95]',0,1));
    end

end

Bias_bhat = nanmean(bhat)-b;
SD_bhat = nanstd(bhat);
Mean_s = nanmean(s);
Mean_sw = nanmean(sw);
Mean_s0 = nanmean(s0);
Size_t = Size_t/r;
Size_tw = Size_tw/r;
Size_t0 = Size_t0/r;
Ratio = Mean_s/SD_bhat;
Ratiow = Mean_sw/SD_bhat;
Ratio0 = Mean_s0/SD_bhat;

Median_numiter = median(numit);
Freq_NonConv = sum(isnan(bhat));

Reject_J = nanmean(pv < 0.05);

tb = 1;
fprintf(tb, '& %8.4f & %8.4f & %8.4f & %8.4f & %8.4f & %8.4f & %8.4f & %8.4f & %8.4f\\\\ \n', a,Ratio,Ratiow,Ratio0,Size_t(2),Size_tw(2),Size_t0(2),Reject_J,Median_numiter);

end




