clc;
close all;
clear all;
warning off;

% addpath('/Applications/MATLAB_R2014a.app/toolbox/poi_library')

IndivN=100;%100:100:500; % Choose n? of Ind

for nn=1:length(IndivN)
nn
TimeT=10;               % Choose Time periods
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);
%% Seeting time series coeff
 eta   = randn(N,3);
 sig   = eta.*eta+0.25;
 phi   = 0.05;%0.05 + (0.95-0.05).*rand(1,N); % parameter of lag of dependent variable

rho = 0; % for time dependensices in errors

%% Create autocorrelated errors (with AR(1))
u     = zeros(T,N);
e     = randn(T,N);
u(1,:) = e(1,:);

alp_i = (randn(N,1)+1);
beta_i =0;%((0.08-(-0.08)).*rand(N,1) + (-0.08)); %time trends

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);
% v2    = randn(T,N);
% v3    = randn(T,N);
x1    = zeros(T,N);
% x2    = zeros(T,N);
% x3    = zeros(T,N);

%% Main DGP
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,:);
    x1(tt,:)= alp_i'.*(J_N-phi)+beta_i'*trend_i(tt,:).*(J_N-phi)+phi.*x1(tt-1,:)+(1-phi.^2).^(1/2).*(v1(tt,:).*sqrt(sig(:,1)'));
 end

%% Beta / null / alternative
% mu    = [1,1,1];
% sigma = [0.03,0,0; 0,0.02, 0; 0,0,0.01];
% beta_i= mvnrnd(mu,sigma,N);
beta_i= ones(N,1);

Y     = zeros(T,N);
X     = zeros(T*N,1);
for ii = 1:N
    X(T*(ii-1)+1:ii*T,:) = [x1(:,ii)];
    Y(:,ii)              = alp_i(ii)+X(T*(ii-1)+1:ii*T,:)*beta_i(ii,:)' + u(:,ii);    
end
Y   = reshape(x1(2:end,:),N*(T-1),1);
Y1  = reshape(x1(1:end-1,:),N*(T-1),1);
TT=reshape(Trend(2:end,:),N*(T-1),1);

%% Model
[bw,sef,sw,ew,br_1,ser_1,sr,er2,hw,FIX,RND,EFIX,ERND,R2w,R2r]=panel(Y,[Y1 TT_block],N); % Fixed Effects and Trend
 

[bw1,sef,sw,ew,br_2,ser_1,sr,er2,hw,FIX,RND,EFIX,ERND,R2w,R2r]=panel(Y,Y1,N); % Fixed Effects


Z=[Y Y1];                       %No Deterministics
[ b, seb, ty, e1] = OLS( Z );


Output(Rep,:)=[bw(1) bw1(1) b];

end

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

kkkk=nn+3;
% BaseName='Output_';
% FileName=[BaseName,num2str(kkkk)];
%  
% %xlwrite(FileName,results);
% xlswrite(FileName,results)


end
results