function NumericalExperiments(model, parameter,filename,R,n,m,psi,gamma)

warning('off','all');

% memory preallocation
RPMLE = zeros(R,1); RSE_PMLE = RPMLE; RCI95_PMLE = RPMLE; RCI90_PMLE = RPMLE;
RGMM1 = zeros(R,1); RSE_GMM1 = RGMM1; RCI95_GMM1 = RGMM1; RCI90_GMM1 = RGMM1;
RGMM2 = zeros(R,1); RSE_GMM2 = RGMM2; RCI95_GMM2 = RGMM2; RCI90_GMM2 = RGMM2; 
RGMM3 = zeros(R,1); RSE_GMM3 = RGMM3; RCI95_GMM3 = RGMM3; RCI90_GMM3 = RGMM3; 



for r=1:R
      % data generation
      [Y X Z] = DataGeneration(model,n,m,psi,gamma,parameter);
      % estimation
      [pmle pmlebc se0] = windmeijeriv(Y,X,X,psi);  
      [gmm1 gmm3   se1] =    mullahyiv(Y,X,Z,psi); se3 = se1;
      [gmm2        se2] =       expiv1(Y,X,Z,psi); 
      
      %[pmle gmm1 gmm3 gmm2],
      
  
      % inference
      % confidence bounds at alpha = .05
      ci95_pmle = pmle-abs(norminv(.025))*se0 <= psi & psi <= pmle+abs(norminv(.025))*se0; 
      ci95_gmm1 = gmm1-abs(norminv(.025))*se1 <= psi & psi <= gmm1+abs(norminv(.025))*se1;
      ci95_gmm3 = gmm3-abs(norminv(.025))*se3 <= psi & psi <= gmm3+abs(norminv(.025))*se3;
      ci95_gmm2 = gmm2-abs(norminv(.025))*se2 <= psi & psi <= gmm2+abs(norminv(.025))*se2;
      % confidence bounds at alpha = .10
      ci90_pmle = pmle-abs(norminv(.050))*se0 <= psi & psi <= pmle+abs(norminv(.050))*se0; 
      ci90_gmm1 = gmm1-abs(norminv(.050))*se1 <= psi & psi <= gmm1+abs(norminv(.050))*se1;
      ci90_gmm3 = gmm3-abs(norminv(.050))*se3 <= psi & psi <= gmm3+abs(norminv(.050))*se3;
      ci90_gmm2 = gmm2-abs(norminv(.050))*se2 <= psi & psi <= gmm2+abs(norminv(.050))*se2;
      % writing away results
      RPMLE(r,:) = pmle'; RSE_PMLE(r,:) = se0'; RCI95_PMLE(r,:) = ci95_pmle'; RCI90_PMLE(r,:) = ci90_pmle';
      RGMM1(r,:) = gmm1'; RSE_GMM1(r,:) = se1'; RCI95_GMM1(r,:) = ci95_gmm1'; RCI90_GMM1(r,:) = ci90_gmm1';
      RGMM3(r,:) = gmm3'; RSE_GMM3(r,:) = se3'; RCI95_GMM3(r,:) = ci95_gmm3'; RCI90_GMM3(r,:) = ci90_gmm3';
      RGMM2(r,:) = gmm2'; RSE_GMM2(r,:) = se2'; RCI95_GMM2(r,:) = ci95_gmm2'; RCI90_GMM2(r,:) = ci90_gmm2';  
end

% saving output
save(filename);

median_r  = [   median(RPMLE);    median(RGMM1);    median(RGMM2);    median(RGMM3)];
stddev_r  = [      std(RPMLE);       std(RGMM1);       std(RGMM2);       std(RGMM3)]; 
ci95_r    = [mean(RCI95_PMLE); mean(RCI95_GMM1); mean(RCI95_GMM2); mean(RCI95_GMM3)];
ci90_r    = [mean(RCI90_PMLE); mean(RCI90_GMM1); mean(RCI90_GMM2); mean(RCI90_GMM3)];
display('--------------------------------------------------------------');
output = [median_r, stddev_r, ci90_r, ci95_r];   display(output);
display('--------------------------------------------------------------');

