clc;
close all;
clear all;
warning off;

IndivN=200;   % N

for nn=1:length(IndivN)

TimeT=5;    % T
results=zeros(length(TimeT),3);

for kk=1:length(TimeT)
    
Replicas=50;
Output=zeros(Replicas,3);

for Rep=1:Replicas

T=TimeT(kk);
N=IndivN(nn);

id=1:N;
id=repmat(id,T,1);
id=reshape(id,[],(T)*N)';


%add years to cluster the standard errors
year=(1:T)';
yearvec=repmat(year,N,1);
trend_i=(1:1:T)';

Trend=repmat(trend_i, 1, N);


%% DGP 2
inovs=randn(T,N);
Mdl = garch('Constant',1,'GARCH',0.9,'ARCH',0.05);
[V,Yn] = simulate(Mdl,T,'NumPaths',N);

%E_vec01=sqrt(V).*inovs;
E_vec01=inovs;
E_vec2=randn(T,N);
%alp_ii = mean(E_vec2);
alp_ii = (randn(1,N));
x=zeros(T,N);
YY=zeros(T,N);

for t=2:T
x(t,:)=0+0.25*x(t-1,:)+E_vec2(t,:);   % dynamics of predictor variable 
YY(t,:)=alp_ii + 0*x(t,:)+E_vec01(t,:);  % Change to other model
end



%% Model
%   - 'po': pooling estimation
%   - 'fe': fixed effects (whitin) estimation
%   - 'be': between estimation
%   - 're': random effects GLS estimation


Y   = reshape(YY(1:end,:),N*T,1);
Y1  = reshape(x(1:end,:),N*T,1);
TT  =reshape(Trend(1:end,:),N*T,1);


Po = panel(id, yearvec, Y, Y1, 'po', 'cluster', yearvec);
Fe = panel(id, yearvec, Y, Y1, 'fe', 'cluster', yearvec);
Re = panel(id, yearvec, Y, Y1, 're', 'cluster', yearvec);


%Output(Rep,:)=[Po.coef(1)-phi Fe.coef(1)-phi Re.coef(1)-phi];
Output(Rep,:)=[Po.coef(1) Fe.coef(1) Re.coef(1)];
end

results(kk,:)=mean(Output);
end
results
end