clc;
close all;
clear all;
warning off;
addpath 'TSM_Matlab_Code';
%addpath('/Applications/MATLAB_R2014a.app/toolbox/poi_library')

IndivN=200;   % N

for nn=1:length(IndivN)

TimeT=10;    % 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);

TT_block= kron(eye(N),trend_i(2:end));




J_N   = ones(1,N);
v1    = randn(T,N);
x    = zeros(T,N);

% %% DGP1
% for tt=2:T
%     u(tt,:) = rho*u(tt-1,:)+(1-rho^2)^(1/2)*e(tt,:);
%     %u(tt,:) = rho*u(tt-1,:)+e(tt,:);
%     x(tt,:)= alp_i'.*(J_N-phi)...
%     +beta_i'*trend_i(tt,:).*(J_N-phi)...
%     +phi.*x(tt-1,:)...
%     +(1-phi.^2).^(1/2).*(v1(tt,:).*sqrt(sig(:,1)'));
% end

%% 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_i = (randn(N,1)+1);
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