function create_statistics(inputfile)
%INPUTS:
%   inputfile  ... mat file created by mc_ar1
%   outputfile ... name of cls-file being created

  disp([10 '# Open file: ' inputfile])

  load(inputfile);  %load MC data from mat-file, this deletes all existing variables
  
  disp(['    rep   = ' num2str(rep)])
  disp(['    N     = ' num2str(N)])
  disp(['    T     = ' num2str(T)])
  disp(['    R0    = ' num2str(R0)])
  disp(['    Rmax  = ' num2str(Rmax)])
  disp(['    beta0 = ' num2str(beta0)])
  
  outputfile=[outputfile '-stats.csv'];
  disp(['# Create summary statistics and write to: ' outputfile])  
  
%open output file
  fid = fopen(outputfile,'wt+');

%write parameters:
  fprintf(fid, 'MC-STATIC\n\n');
  fprintf(fid, 'rep, %d\n', rep);
  fprintf(fid, 'N, %d\n', N);
  fprintf(fid, 'T, %d\n', T);
  fprintf(fid, 'R0, %d\n', R0);
  fprintf(fid, 'Rmax, %d\n', Rmax);
  fprintf(fid, 'beta0, %f\n', beta0);

%write header:
fprintf(fid, 'estimator, #factors, bias, std, rmse, mean est. std 1, mean est. std 2, \n\n');
for R=0:Rmax
  b=beta_est(:,1+R);
  bias_est=mean(b)-beta0;
  std_est=sqrt(mean((b-mean(b)).^2));
  rmse_est=sqrt(mean((b-beta0).^2));
  fprintf(fid, 'LS-factor, R=%d, %f, %f, %f, %f, %f, \n',R,bias_est,std_est,rmse_est,mean(est_std1(:,1+R)),mean(est_std2(:,1+R)));
end

fprintf(fid, '\n');
for R=1:Rmax
  b=beta_BC1(:,R);
  bias_est=mean(b)-beta0;
  std_est=sqrt(mean((b-mean(b)).^2));
  rmse_est=sqrt(mean((b-beta0).^2));
  fprintf(fid, 'Bias Corrected 1, R=%d, %f, %f, %f, \n',R,bias_est,std_est,rmse_est);
end

fprintf(fid, '\n');
for R=1:Rmax
  b=beta_BC2(:,R);
  bias_est=mean(b)-beta0;
  std_est=sqrt(mean((b-mean(b)).^2));
  rmse_est=sqrt(mean((b-beta0).^2));
  fprintf(fid, 'Bias Corrected 2, R=%d, %f, %f, %f, \n',R,bias_est,std_est,rmse_est);
end

fprintf(fid, '\n');
for i=1:2
  b=beta_CCE_est(:,i);
  bias_est=mean(b)-beta0;
  std_est=sqrt(mean((b-mean(b)).^2));
  rmse_est=sqrt(mean((b-beta0).^2));
  if i==1
    fprintf(fid, 'mean group CCE, , %f, %f, %f, \n',bias_est,std_est,rmse_est);
  else  
    fprintf(fid, 'pooled CCE, , %f, %f, %f, \n',bias_est,std_est,rmse_est);
  end
end



%QUANTILES:

if mod(rep,40)==0 %NUMBER OF REPETITIONS NEEDS TO BE DIVISIBLE BY 40
    
  fprintf(fid, '\n\n QUANTILES of sqrt(NT)*(ESTIMATOR-beta0) for R>=R0\n\n');
  fprintf(fid, 'estimator, #factors, 2.5%%, 5%%, 10%%, 25%%, 50%%, 75%%, 90%%, 95%%, 97.5%%, \n\n');

  for R=R0:Rmax 
    beta_est_norm=sort(sqrt(N*T)*(beta_est(:,R+1)-beta0));
  
    q2=ceil(rep*2.5/100); 
    qq2=(beta_est_norm(q2)+beta_est_norm(q2+1))/2;
    q5=ceil(rep*5/100); 
    qq5=(beta_est_norm(q5)+beta_est_norm(q5+1))/2;
    q10=ceil(rep*10/100); 
    qq10=(beta_est_norm(q10)+beta_est_norm(q10+1))/2;
    q25=ceil(rep*25/100); 
    qq25=(beta_est_norm(q25)+beta_est_norm(q25+1))/2;
    q50=ceil(rep*50/100); 
    qq50=(beta_est_norm(q50)+beta_est_norm(q50+1))/2;
    q75=ceil(rep*75/100); 
    qq75=(beta_est_norm(q75)+beta_est_norm(q75+1))/2;
    q90=ceil(rep*90/100); 
    qq90=(beta_est_norm(q90)+beta_est_norm(q90+1))/2;
    q95=ceil(rep*95/100); 
    qq95=(beta_est_norm(q95)+beta_est_norm(q95+1))/2;
    q97=ceil(rep*97.5/100); 
    qq97=(beta_est_norm(q97)+beta_est_norm(q97+1))/2;

    fprintf(fid, 'LS-factor, R=%d, %f, %f, %f, %f, %f, %f, %f, %f, %f \n',R,qq2,qq5,qq10,qq25,qq50,qq75,qq90,qq95,qq97);
  end
end    


fclose(fid);  
  





return