%  Time periods are B, A, 1, 2, 3, 4

rng('default')

delete(logfile_name);
diary(logfile_name);
diary on

n=a(1);
k=a(2);
nrep=a(3);

n=10*n

beta=zeros(k,1);
beta(1)=1;
beta(2)=1;
gamma=[1.0; 0.5]; 
% gamma=0*gamma;

thhat1=zeros(nrep,2+k);
thhat2=zeros(nrep,2+k);
thhat3=zeros(nrep,2+k);
fvals=zeros(nrep,1);

ms=zeros(nrep,4*(k*3+4));
ms=zeros(nrep,4);
pvals=zeros(nrep,5);
should_scale=1;

Results(nrep) = struct('S', []);

for irep=1:nrep
tic
   fprintf('%5i',irep)
    xB=randn(n,k);
    xA=randn(n,k);
    x1=randn(n,k);
    x2=randn(n,k);
    x3=randn(n,k);
    x4=randn(n,k);
    aa=0.5;
    
    xB(:,2:k)=((1-aa)*xB(:,2:k)+aa*repmat(xB(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
    xA(:,2:k)=((1-aa)*xA(:,2:k)+aa*repmat(xA(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
    x1(:,2:k)=((1-aa)*x1(:,2:k)+aa*repmat(x1(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
    x2(:,2:k)=((1-aa)*x2(:,2:k)+aa*repmat(x2(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
    x3(:,2:k)=((1-aa)*x3(:,2:k)+aa*repmat(x3(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
    x4(:,2:k)=((1-aa)*x4(:,2:k)+aa*repmat(x4(:,1),1,k-1))/sqrt((1-aa)^2+aa^2);
%     
    

    alpha=1*(xB(:,1)+xA(:,1)+x1(:,1)+x2(:,1)+x3(:,1)+x4(:,1))/2;
    zer=zeros(n,1);
    theta0=[gamma; beta];
    
    yB=generatelogit([alpha zer zer xB],[1;theta0]);
    yA=generatelogit([alpha yB zer xA],[1;theta0]);
%     yB=yB*0+1;
%     yA=yA*0+1;
    
    
    
    y1=generatelogit([alpha yA yB x1],[1;theta0]);
    y2=generatelogit([alpha y1 yA x2],[1;theta0]);
    y3=generatelogit([alpha y2 y1 x3],[1;theta0]);
    y4=generatelogit([alpha y3 y2 x4],[1;theta0]);
 
    clear alpha
    
    mean([yB yA y1 y2 y3 y4])
    
%     for iB=0:1
%         for iA=0:1
            for i1=0:1
                for i2=0:1
                    for i3=0:1
                        for i4=0:1
%                             aa=mean((yB==iB).*(yA==iA).*(y1==i1).*(y2==i2).*(y3==i3).*(y4==i4));
                            aa=mean((y1==i1).*(y2==i2).*(y3==i3).*(y4==i4));
%                             fprintf('%1i',iB)
%                             fprintf('%1i',iA)
                            fprintf('%1i',i1)
                            fprintf('%1i',i2)
                            fprintf('%1i',i3)
                            fprintf('%1i',i4)
                            fprintf('& %6.3f \\%%\n',100*aa);
                        end
                    end
                end
            end
%         end
%     end
    
    stop
    
    yy=[y1;y2;y3;y4];
    xx=[ones(n,1) yA yB x1;...
        ones(n,1) y1 yA x2;...
        ones(n,1) y2 y1 x3;...
        ones(n,1) y3 y2 x4];
    [th0,s,f,H,V,scores]=estimate_logit(yy,xx);
    thhat1(irep,:)=th0(2:length(th0))';

    
%     [m,S]=moments_small(yB,yA,y1,y2,y3,y4,xB,xA,x1,x2,x3,x4,theta0,should_scale);
%     Results(irep).S=S/n;
%     ms(irep,:)=m;
% 
%     chi=n*(m/S)*m';
%     chi=n* m*inv(S)*m';
%     fprintf('\nall Chi2 =  %6.3f',chi)
%     fprintf('    pval =   %6.3f\n\n',1-cdf('Chisquare',chi,length(m)))
%     pvals(irep,1)=1-cdf('Chisquare',chi,length(m));
%     
%     
%     for i=1:length(m)
%         t=sqrt(n)*m(i)/sqrt(S(i,i));
%         fprintf('t= %6.3f\n',t)
%         pvals(irep,1+i)=1-cdf('Chisquare',t^2,1);
%     end
%     stop

%     for i=1:4
%         ii=(1+(i-1)*13):(i*13);
%         mm=m(ii);
%         SS=S(ii,ii);
%         chi=n*(mm/SS)*mm';
%         prob=1-cdf('Chisquare',chi,length(mm));
%          pvals(irep,1+i)=prob; 
%         fprintf(' i = %2i',i)
%          fprintf('    Chi2 =  %6.3f',chi)
%     fprintf('    pval =   %6.3f\n\n',prob)
%     end
%        stop mean(pvals)


    
    
    thstart=th0(2:length(th0));
    [m,S]=moments(yB,yA,y1,y2,y3,y4,xB,xA,x1,x2,x3,x4,thstart,should_scale);

    W0=inv(diag(diag(S)));
%     W0=inv(S);
%     W0=eye(size(S,1));
    [th1 fval]=estimate_GMM(yB,yA,y1,y2,y3,y4,xB,xA,x1,x2,x3,x4,thstart,W0,true,should_scale);
    thhat2(irep,:)=th1';
% 
% %     
    xx1=zeros(4*n,k+2);
    yy1=zeros(4*n,1);
    for i=1:n
       xx1((i-1)*4+1,:)=[yA(i) yB(i) x1(i,:)];
       xx1((i-1)*4+2,:)=[y1(i) yA(i) x2(i,:)];
       xx1((i-1)*4+3,:)=[y2(i) y1(i) x3(i,:)];
       xx1((i-1)*4+4,:)=[y3(i) y2(i) x3(i,:)];
       
       yy1((i-1)*4+1,:)=y1(i);
       yy1((i-1)*4+2,:)=y2(i);
       yy1((i-1)*4+3,:)=y3(i);
       yy1((i-1)*4+4,:)=y4(i);
    end
    nn1=[(1:n)'*4-3 (1:n)'*4];  
   
    [blogit_FE0,slogit_FE0,fval,H,V,scores]=estimate_logit_FE(yy1,xx1,nn1,th0,1,false);
    thhat3(irep,:)=blogit_FE0';

    fprintf('\n')
    

toc
end

pri([median(thhat1); median(thhat2); median(thhat3)])

% SS=zeros(size(ms,2),size(ms,2));
% for irep=1:nrep
%     SS=SS+Results(irep).S;
% end
% 
% S=SS/nrep;
% m=mean(ms);
% chi=nrep* m*inv(S)*m';
% fprintf('\nall Chi2 =  %6.3f',chi)
% fprintf('    pval =   %6.3f\n\n',1-cdf('Chisquare',chi,length(m)))
% 
% mean(pvals)
   
save(outfile);
