rng('default')

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

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

beta=zeros(k,1);
beta(1)=1;
beta(2)=1;
gamma=1;

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

for irep=1:nrep
tic
   fprintf('%5i',irep)
    x0=randn(n,k);
    x1=randn(n,k);
    x2=randn(n,k);
    x3=randn(n,k);
    aa=0.5;
    
    x0(:,2:k)=((1-aa)*x0(:,2:k)+aa*repmat(x0(:,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);

    alpha=0*(x0(:,1)+x1(:,1)+x2(:,1)+x3(:,1))/2;
    y0=exp(alpha+x0*beta)./(1+exp(alpha+x0*beta))-rand(n,1)>0;
    y1=exp(alpha+x1*beta+gamma*y0)./(1+exp(alpha+x1*beta+gamma*y0))-rand(n,1)>0;
    y2=exp(alpha+x2*beta+gamma*y1)./(1+exp(alpha+x2*beta+gamma*y1))-rand(n,1)>0;
    y3=exp(alpha+x3*beta+gamma*y2)./(1+exp(alpha+x3*beta+gamma*y2))-rand(n,1)>0;
    theta0=[gamma; beta];
    clear alpha
    
    yy=[y1;y2;y3];
    xx=[ones(n,1) y0 x1;...
        ones(n,1) y1 x2;...
        ones(n,1) y2 x3];
    [th0,s,f,H,V,scores]=estimate_logit(yy,xx);
    thhat0(irep,:)=th0(2:length(th0))';
    
    thstart=th0(2:length(th0));
    
    [m,S]=moments_full_scaled(y0,y1,y2,y3,x0,x1,x2,x3,thstart);
    W0=inv(diag(diag(S)));
    [th1 fval]=estimate_full_scaled(y0,y1,y2,y3,x0,x1,x2,x3,thstart,W0,true);
    thhat1(irep,:)=th1';
    
    xx1=zeros(3*n,k+1);
    yy1=zeros(3*n,1);
    for i=1:n
       xx1((i-1)*3+1,:)=[y0(i) x1(i,:)];
       xx1((i-1)*3+2,:)=[y1(i) x2(i,:)];
       xx1((i-1)*3+3,:)=[y2(i) x3(i,:)];
       yy1((i-1)*3+1,:)=y1(i);
       yy1((i-1)*3+2,:)=y2(i);
       yy1((i-1)*3+3,:)=y3(i);
    end
    nn1=[(1:n)'*3-2 (1:n)'*3];  
   
    [blogit_FE0,slogit_FE0,fval,H,V,scores]=estimate_logit_FE(yy1,xx1,nn1,th0,1,false);
    thhat2(irep,:)=blogit_FE0';

    fprintf('\n')
    

toc
end
   
save(outfile);
