function output =  EB_vs_BQML(T,n,n_com_factors,q,mu,load_type,dist_type,delta, tau, n_blocks,factor_per_block, block_type,theta, share,k);



%%% Inputs storing
output.T = T;
output.n = n;
output.r_c = n_com_factors;
output.q = q;
output.mu = mu;
output.load_dist = load_type;
output.idio_dist = dist_type;
output.idio_auto = delta;
output.idio_cross = tau;
output.n_blocks = n_blocks;
output.r_b = factor_per_block;
output.block_dist = block_type;
output.signalt_to_noise = theta;
output.NA_share = share;
output.MA_smooth = k;

output

%%%% Generate and store blocks properties
%%% Generate the indices for each block according to the chosen type
%%% (homogeneous/heterogeneuos)
if (block_type == 1)
    min_n = n/n_blocks - 1;
    max_n = n/n_blocks;
else
    min_n = max(n/100,factor_per_block*2); 
    max_n = 0.8*n;
end
block_generator = round(randfixedsum(n_blocks, 1, n, min_n, max_n));        % Simulate a series of number summing to the number of series
diff = sum(block_generator)- n;                                             % Calculate mismatch
[v, ind] = max(block_generator);                                            
block_generator(ind) = block_generator(ind) - diff;                         % Fix mismatch 

blocks = {};                                                                % Store the block info
n_cum = 0;
for b = 1:n_blocks
     label_block = strcat('B',string(b));                                   
     blocks{b,1} = label_block;                                             % Define block labels
     blocks{b,2} = n_cum+1:n_cum + block_generator(b);                      % Store reference index
     blocks{b,3} = factor_per_block;                                        % Store n_facto per block
     blocks{b,4} = 1;                                                       % Store autoregressive order
     n_cum = n_cum + block_generator(b);
end    


%%%  Simulate common factors  %%%
r = n_com_factors;
u=randn(T+10,q);                                                        	% common shocks 
uu=custom_Standardize(u);                                                       % standardize common shocks 
            % ------------------------------------------------------------- % VAR Common Factors
A=diag(custom_uniform(q,.5,.8));                                                % diagonal elements of A(L)
A(~eye(q))=custom_uniform(q*(q-1),.1,.4);                                       % extra diagonal elements of A(L)
mut=max(eig(A));                                                            % Max Root of the polynomial
A=mu*A/(mut);                                                               % Fix max eigenvalue
            % ------------------------------------------------------------- % 
ff0=uu; for tt=2:T+10; ff0(tt,:)=(A*ff0(tt-1,:)')'+uu(tt,:); end            % Generate common factors
ff1=custom_lag(ff0,(r/q)-1,0);                                                  % take lags of the common factors
ff1=ff1(end-T+1:end,:);                                                     % keeps only the last T observations 


%%%  Simulate block-factors   %%%
factors = {};

for b = 1:n_blocks

r = blocks{b,3};
u=randn(T+10,q);                                                        	% common shocks 
uu=custom_Standardize(u);                                                       % standardize common shocks 
            % ------------------------------------------------------------- % VAR Common Factors
A=diag(custom_uniform(q,.5,.8));                                                % diagonal elements of A(L)
A(~eye(q))=custom_uniform(q*(q-1),.1,.4);                                       % extra diagonal elements of A(L)
mut=max(eig(A));                                                            % Max Root of the polynomial
A=mu*A/(mut);                                                               % Fix max eigenvalue
            % ------------------------------------------------------------- % 
ff0_b =uu; for tt=2:T+10; ff0_b(tt,:)=(A*ff0_b(tt-1,:)')'+uu(tt,:); end     % Generate common factors
ff1_b=custom_lag(ff0_b,(r/q)-1,0);                                              % take lags of the common factors
ff1_b=ff1_b(end-T+1:end,:);                                                 % keeps only the last T observations

factors{b} = ff1_b;


end


factor_matrix = [ff1,cell2mat(factors)];


%%%  Simulate idiosyncratic components  %%%

if delta==0     % --------------------------------------------------------- % Autocorrelation idiosyncratic component
    d=zeros(n,1);                                                           % d = 0
else; d=custom_uniform(n,delta,1-2*delta);                                      % d ~ U(delta,1-2*delta)
end             % --------------------------------------------------------- %
TT1=diag(rand(n,1)+0.5);                                                    % variance idiosyncratic shocks
if tau==0; TT2=eye(n); else; TT2=toeplitz(tau.^(0:n-1)); end                % covariance matrix idiosyncratic shocks
TT2=TT2-triu(TT2,11)-tril(TT2,-11);                                         % keep only 10 correlations

switch dist_type
    case 'normal'
        v = randn(T+10,n)*sqrt(TT1)*cholcov(TT2);                                % idiosyncratic shocks 
    case 'laplace'
        v = randl(T+10,n)*sqrt(TT1)*cholcov(TT2);
    case 'student_3'
        v = trnd(3,T+10,n)*sqrt(TT1)*cholcov(TT2);
    case 'student_1'
        v = trnd(1,T+10,n)*sqrt(TT1)*cholcov(TT2);
end


e=v; for ii=1:n; e(:,ii)=filter(1,[1 -d(ii)],v(:,ii)); end                  % Idiosyncratic component
e(1:10,:)=[];                                                               % eliminate the first 10 observations
e=custom_center(e);                                                             % center idiosyncratic component


%%% ================================= %%%
%%%  Common components and variables  %%%
%%% ================================= %%%
loadings = {};

switch load_type
    case 'unif'
          LL1_c = -1 + 2*rand(n,n_com_factors);                            % factor loadings
          for b =1:n_blocks
            loadings{b} = -1+2*rand(length(blocks{b,2}),blocks{b,3});
          end
    case 'normal'
          LL1_c = randn(n,n_com_factors);                                  % factor loadings
          for b =1:n_blocks
              loadings{b} = randn(length(blocks{b,2}),blocks{b,3});
          end                                                              % idiosyncratic shocks 
    case 'laplace'
          LL1_c = randl(n,n_com_factors);                                        % factor loadings
          for b =1:n_blocks
              loadings{b} = randl(length(blocks{b,2}),blocks{b,3});
          end    
    case 'student_3'
          LL1_c = trnd(3,n,n_com_factors);                                       % factor loadings
          for b =1:n_blocks
              loadings{b} = trnd(3,length(blocks{b,2}),blocks{b,3});
          end
    case 'student_1'
          LL1_c = trnd(1,n,n_com_factors);                                       % factor loadings
          for b =1:n_blocks
             loadings{b} = trnd(1,length(blocks{b,2}),blocks{b,3});
          end
end

LL1 = [LL1_c,blkdiag(loadings{:})];


chi = custom_center(factor_matrix*LL1(1:n,:)' );                                % common components
% ------------------------------------------------------------------------- % fix signal to noise ratio
temp=var(chi)./var(e);                                                      % ratio of variances
theta2=theta+custom_center(custom_uniform(n,theta-.25,theta+.25));                  % draw signal to noise ratio
for ii=1:n; chi(:,ii)=(1/sqrt(theta2(ii)))*chi(:,ii)/sqrt(temp(ii)); end    % impose signal to noise ratio
% ------------------------------------------------------------------------- %
Y=chi+e;                                                                    % variables

idx = randsample(size(Y,1)*size(Y,2),floor(share*size(Y,1)*size(Y,2))) ;
Y(idx) = NaN;


P.r = n_com_factors;% number of factors
P.p = 1; % number of lags in the factor VAR
P.max_iter = 100; % maximal number of EM iterations
P.thresh = 1e-4;
optNaN.method = 1; 
optNaN.k = k;

try
startTime = clock;
res_EB = EB_est(Y,P,optNaN,blocks);
elapsedTime_EB = etime(clock, startTime);

%%%% Model score
output.trace_stat_C_EB = trace_stat(LL1,res_EB.lambda);
output.trace_stat_F_EB = trace_stat(factor_matrix,res_EB.factors);
output.MSE_EB = mean(mean((LL1*factor_matrix'-res_EB.lambda*res_EB.factors').^2));

%%%% Model performance
output.elapsedTime_EB = elapsedTime_EB;

catch 
    
output.trace_stat_C_EB = nan;
output.trace_stat_F_EB = nan;
output.MSE_EB = nan;

output.elapsedTime_EB = nan;
    
end    

try 
startTime = clock;
res_BQML = BDFM_estimator(Y,P,blocks);
elapsedTime_BQML = etime(clock, startTime);

%%%% Model score
output.trace_stat_C_BQML = trace_stat(LL1,res_BQML.C);   
output.trace_stat_F_BQML = trace_stat(factor_matrix,res_BQML.F);
output.MSE_BQML = mean(mean((LL1*factor_matrix'-res_BQML.C*res_BQML.F').^2));

%%%% Model performance
output.elapsedTime_BQML = elapsedTime_BQML;

catch
%%%% Model score

output.trace_stat_C_BQML = nan;   
output.trace_stat_F_BQML = nan;
output.MSE_BQML = nan;

%%%% Model performance
output.elapsedTime_BQML = nan;
end
    





output;


