function mc_static(rep,N,T,R0,Rmax,beta0,outputfile)
%INPUTS:
%   rep   ... number of repetitions
%   N,T   ... panel dimensions
%   R0    ... true number of factors
%   Rmax  ... maximal number of factors in estimation
%   beta0 ... true regressions parameter
%   outputfile ... file name, where to write results

disp([10 '# Running MC-STATIC with the following parameters:'])
disp(['    rep   = ' num2str(rep)])
disp(['    N     = ' num2str(N)])
disp(['    T     = ' num2str(T)])
disp(['    R0    = ' num2str(R0)])
disp(['    Rmax  = ' num2str(Rmax)])
disp(['    beta0 = ' num2str(beta0)])
disp(['# Starting MC runs:'])

dof=5;      %degrees of freedom for t-distribution of errors
MAtheta=1;  %MA-parameter of errors


for r=1:rep
  if mod(r,min(50,ceil(rep/20)))==0
    disp(['    ' num2str(r) ' / ' num2str(rep)]);
  end
  
  %Define DGP:
  
  lambda=1+randn(N,R0);           %factor loading
  lambdaX=1+randn(N,R0);   %factor loading that enters into regressors
  ff = randn(T+1,R0);        
  f = ff(2:T+1,:);                  %factors
  f1 = ff(1:T,:);               %lagged factors
  fX = f1;                %factors that enter into regressors
  
  v = tinv(rand(N,T+1),dof);   %shocks are t-distributed with degrees of freedom = dof
  v0 = v(:,2:T+1); 
  v1 = v(:,1:T); 
  e = (1 + MAtheta^2)^(-0.5) * (v0 + MAtheta*v1);   %MA1 process for errors
  
  X = 1 + randn(N,T) + (lambda+lambdaX)*(f+fX)';    %regressors
  Y = beta0 * X + lambda*f' + e;     %outcome variable
  XX(1,:,:)=X;   %only one regressors

  %calculate OLS estimator:
  
  beta_est(r,1)=inv(trace(X*X'))*trace(X*Y');
  res=Y-beta_est(r,1)*X;
  est_std1(r,1)=sqrt(var(res(:))*inv(trace(X*X')));
  est_std2(r,1)=sqrt( inv(trace(X*X')) * trace((res.*X)*(res.*X)')  * inv(trace(X*X')) );
  
  %calculate QMLE with different R:
  
  for R=1:Rmax
    start_beta=0;
    repMIN=3;
    repMAX=10;
    [beta,exitflag,lambda,f,Vbeta1,Vbeta2,Vbeta3,bcorr1,bcorr2,bcorr3]=LS_factor(Y,XX,R,'silent',10^-8,'m1',start_beta,repMIN,repMAX,2,2);
    if exitflag<=0
      disp('QMLE OPTIMIZATION PROBLEM, may need to increase repMIN or repMAX or optimset-parameters in function qmle')  
    end 
    beta_est(r,1+R)=beta;        %LS-estimator using R factors in estimation
    beta_BC1(r,R)=beta-bcorr2-bcorr3; %bias corrected estimator (only using static biases)
    beta_BC2(r,R)=beta-bcorr1-bcorr2-bcorr3; %bias corrected estimator (using all bias terms)
    est_std1(r,1+R)=sqrt(Vbeta1); %estimated standard error, assuming homoscedasticity
    est_std2(r,1+R)=sqrt(Vbeta2); %estimated standard error, assuming heteroscedasticity
  end

  %calculate common correlated effets estimator (CCE):
  
  w=[mean(Y)',mean(X)'];
  Mw=eye(T)-w*inv(w'*w)*w';

  for i=1:N
    b(i)=inv(X(i,:)*Mw*X(i,:)')*X(i,:)*Mw*Y(i,:)';  %individual slope CCE
  end
  beta_CCE_est(r,1)=mean(b);                             %mean group CCE (CCEMG)
  beta_CCE_est(r,2)=inv(trace(X*Mw*X'))*trace(X*Mw*Y');  %pooled CCE (CCEP)
end  

%delete all variables, except those needed to create MC tables:
clearvars -except rep N T R0 Rmax beta0 dof MAthete outputfile start_beta repMIN repMAX beta_est beta_BC1 beta_BC2 est_std1 est_std2 beta_CCE_est beta_CCE_est

%save results:
disp(['# Save result to file: ' outputfile])
save([outputfile '.mat']);
clear all

return

%alternative way to name output file:
%save(['results_rep=' num2str(rep) '_N=' num2str(N) '_T=' num2str(T) '_R0=' num2str(R0) '_Rmax=' num2str(Rmax)  '_beta0=' num2str(beta0) '_rhof=' num2str(rhof) '.mat'])




