%% Date: 02/02/2022
%% Panel reg with ...
%% Dates in datetime form
%% Id is an identifier
function [beta,tstat,V_T,numobs,numgroup,rsquared]=panel_reg(dates,id,y,x,typereg,lag)

%% typereg:={1-Pooled,2->FE within}
%% Lag=amount of lags in NW regression

alldates=unique(dates); 
T=size(alldates,1);
secids=unique(id); %% that's how we need how many i's we have
N=size(secids,1);
numgroup=N;
%% To get estimator, (X'X)^{-1}X'y,
%% We stack it according to id. First T are 1, second T are 2,..

q=size(x,2);
X=NaN(T*N,q+1); %% q+1 because we'll add ones
Y=NaN(T*N,1);

%% Averages for the FE's
Xhat=NaN(T*N,q+1);
Yhat=NaN(T*N,1);

ids_matrix=NaN(T*N,1);

date_matrix=repmat(alldates,N,1); %% all dates repeated

for i=1:N
 indices=(i-1)*T+1:i*T;
 ids_matrix(indices)=secids(i);
end

%% Run on all ids

for i=1:N
idx1=ismember(id,secids(i)); 
tempdata=[datenum(dates(idx1,:)),y(idx1,:),x(idx1,:)];

idx2=ismember(datenum(date_matrix),tempdata(:,1));
idx3=ismember(ids_matrix,secids(i));
idx=logical(idx2.*idx3);

c=ones(size(tempdata,1),1);
X(idx,:)=[c,tempdata(:,3:end)];
Y(idx,:)=tempdata(:,2);
l=size(tempdata,1);
Xhat(idx,:)=[repmat(nanmean([c,tempdata(:,3:end)],1),l,1)];
Yhat(idx,:)=repmat(nanmean([tempdata(:,2)],1),l,1);


end

%% For F.E
Xhathat=repmat(nanmean(X,1),T*N,1);
Yhathat=repmat(nanmean(Y,1),T*N,1);

%% get rid of NaN's
idx=~(sum(isnan([Y,X]),2)); %% these are the lines we keep

clean=packr([Y,X]);
h_it=NaN(T*N,q+1);

Y=clean(:,1);
X=clean(:,2:end);
numobs=size(X,1);

Xhat=Xhat(idx,:);
Yhat=Yhat(idx,:);
Xhathat=Xhathat(idx,:);
Yhathat=Yhathat(idx,:);
h_t=NaN(T,q+1);

if typereg==0

beta=(X'*X)\X'*Y;

%% DK standard errors
epsilons=Y-X*beta;
SSE=sum(epsilons.*epsilons);
SST=sum( (Y-mean(Y)).^2);
rsquared=1-SSE/SST;

h_it(idx,:)=X.*epsilons;

for i=1:T
 indices=i:T:i+T*(N-1);
 h_t(i,:)=nansum(h_it(indices,:),1);
end

h_t=h_t';

Omega_0=h_t*h_t';

S_T=Omega_0;

Omegas=NaN(q+1,q+1,lag);
    
for l=1:lag
       Omegas(:,:,l)=h_t(:,l+1:end)*h_t(:,1:end-l)';    
end


for l=1:lag
weight=1-l/(lag+1);
Omega_j=reshape(Omegas(:,:,l),q+1,q+1);
S_T=S_T+weight*(Omega_j+Omega_j');
end

%% Var-Covar matrix

V_T=inv(X'*X)*S_T*inv(X'*X);


else

    %% F.E. estimator
    %% first we need averages of each group
    
    Ytag=Y-Yhat+Yhathat;
    Xtag=X-Xhat+Xhathat;
    beta=(Xtag'*Xtag)\Xtag'*Ytag;

    
%% DK standard errors
epsilons=Ytag-Xtag*beta;
SSE=sum(epsilons.*epsilons);
SST=sum( (Ytag-mean(Ytag)).^2);
rsquared=1-SSE/SST;

h_it(idx,:)=Xtag.*epsilons;

for i=1:T
 indices=i:T:i+T*(N-1);
 h_t(i,:)=nansum(h_it(indices,:),1);
end

h_t=h_t';

Omega_0=h_t*h_t';

S_T=Omega_0;

Omegas=NaN(q+1,q+1,lag);
    
for l=1:lag
       Omegas(:,:,l)=h_t(:,l+1:end)*h_t(:,1:end-l)';    
end


for l=1:lag
weight=1-l/(lag+1);
Omega_j=reshape(Omegas(:,:,l),q+1,q+1);
S_T=S_T+weight*(Omega_j+Omega_j');
end

%% Var-Covar matrix

V_T=inv(Xtag'*Xtag)*S_T*inv(Xtag'*Xtag);


end

tstat=beta./sqrt(diag(V_T));
