function [se_var, se_sd] = std_err(id,year,X,w,var)

% Create a balanced matrix of growth rates abd weights, where 0's fill in 
% missing, which is ok b/c will sum over zeros and not count when taking 
% sums

[h1,h2,ii]=unique(id);
[h1,h2,jj]=unique(year);

N=max(ii);
T=max(jj);

q = round(.75*T^(1/3) - 1);

bigX=zeros(N,T);
bigX(sub2ind([N T],ii,jj))= X;
bigO=zeros(N,T);
bigO(sub2ind([N T],ii,jj))=1;
counts=sum(bigO,2);
mX=sum(bigX,2)./sum(bigO,2);
demeanX=(bigX-repmat(mX,1,T)).*bigO;

% Demeaned matrix which is N,T and create a (T_n - 1) to divide by 

Y = reshape(demeanX,N,T);

t = counts - ones(size(counts,1),1);
tt = repmat(t,1,T);
TT = reshape(tt,N,T);

Y = Y./sqrt(TT);

W = zeros(N,T);
W(sub2ind([N T],ii,jj))=w;

% Create inner product of W and Y
Z = W'*Y;

% Difference of innner product and variance, repeated for each variance. Do
% some matrix manipulations so can do all 1...T standard errors in one
% calculation without a loop

% For contemporaneous variance calculations

dif_c = sqrt(T-1)*Z.^2-1/sqrt(T-1)*repmat(var,1,T);
s2 = dif_c.^2;
S2 = sum(s2,2); 

% For autocorrelations calculations

S12 = zeros(T,1);
for k = 1:q
    c = 2*(1 - k/(q+1));
    dif_Tk = sqrt(T-k-1)*Z.^2-1/sqrt(T-k-1)*repmat(var,1,T);
    dif_Tmk = dif_Tk(:,1:end-k);
    dif_Tpk = dif_Tk(:,1+k:end);
    s12 = sum(dif_Tmk.*dif_Tpk,2);
    SS12 = c*s12;
    S12 = S12 + SS12;
end


se_var = sqrt(S2 + S12)/sqrt(T);
se_sd = se_var./(2*sqrt(var));