%% FGLs for diffs in diffs.  
% Estimates using FGLS and reports whether
% reject 
%
% $$H_0: \beta=0$$ 
%
% at 5% level in three ways: 
%
% # plain FGLS allowing arbitrary autocorrelation
% # Hausman and Kuersteiner bias correction
% # bias corrected with size corrected test
%
% Takes as input a data structure with fields:
%  
% * y = outcome variable, size number of observations by 1
% * x = regressors, size number of observations by number of regressors
% * t = time index, size number of observations by 1
% * s = state index, size number of observations by 1
% * p = policy indicator, size number of observations by 1
%
% As in BDM and HK, we work with data aggregated to state-time cells,
% this makes the computations simpler.  The justification is that there
% are so many within state-time observations that we can ignore
% uncertainty from that.  We require that the input data has already been
% aggregated.  This can be accomplished with aggData()
%
% See also AGGDATA, OLSDD

function [bGLS testGLS vGLS] = fglsDD(data)
  alpha = 0.05; % critical level for tests
  crit = abs(norminv(alpha/2));
  
  % check that data has been collapsed to state-time cells
  T = length(unique(data.t));
  S = length(unique(data.s));
  if (length(data.y)~=T*S) 
    error('data must be collapsed to state-time averages');
  end
  x = data.x;
  y = data.y;

  %% ols
  ixx = inv(x'*x); 
  betaP = ixx*x'*y;
  u = y - x*betaP;
  bGLS(4) = betaP(1);
  V = u'*u/length(u)*ixx;
  vGLS(4) = V(1);
  testGLS = zeros(1,4);
  testGLS(4) = abs(bGLS(4)/sqrt(vGLS(4))) > crit;
  
%% Estimate of weighting matrix
% The presence of fixed effects makes doing FGLS
% somewhat tricky. In particular, the natural estimate of the variance
% matrix, 
% 
% $$ \hat{\Sigma} = \sum_{s=1}^S \hat{e}_s \hat{e}_s' $$
%
% will be rank T-1 instead of rank T.  We can proceed as in Kiefer
% (1980).  
  Omega = zeros(T); 
  M1t = eye(T) - ones(T)/T;
  for s=1:S
    Omega = Omega + M1t*u(data.s==s)*u(data.s==s)'*M1t;
  end
  Omega = Omega/S;
  % need generalized inverse of Omega b/c of non full rank
  iO = zeros(T);
  iO(1:T-1,1:T-1) = inv(Omega(1:T-1,1:T-1));
  W = kron(eye(S),M1t)*kron(eye(S),iO)*kron(eye(S),M1t);
%% FGLS
  xf = x(:,1:T); % already projected out state fixed effects 
  V = inv(xf'*W*xf);  
  betaGLS = V * (xf'*W*y);
  bGLS(1) = betaGLS(1);
  testGLS(1) = abs(bGLS(1)/sqrt(V(1))) > crit;
  vGLS(1) = sqrt(V(1));

%% Bias-corrected FGLS
% see Hausman and Kuersteiner
  yt = reshape(y,T,S);
  v = zeros(S,2);
  % must be careful about defining v to make it full rank
  for s=1:S
    v(s,:) = [1 any(x(data.s==s,1)==1)];
  end
  Mv = eye(S) - v*inv(v'*v)*v';
  Stilde = yt*Mv*yt'/trace(Mv);
  B = eye(T);
  B = B(2:T,:);
  SigHat = B*M1t*Stilde*M1t*B';
  
  tran = kron(eye(S),B*M1t);
  z = tran*xf(:,2:size(xf,2));
  ups = tran*xf(:,1);
  ybc = tran*y;
  iSig = inv(SigHat);
  iO = kron(eye(S),iSig);
  Oz = iO - iO*z*inv(z'*iO*z)*z'*iO;
  V = inv(ups'*Oz*ups);
  bGLS(2:3) = V*ups'*Oz*ybc;
  t1 = bGLS(2)/sqrt(V);
  testGLS(2) = abs(t1) > crit;
  vGLS(2:3) = V;
  
%% size corrected test
% see Hausman and Kuersteiner
  sccrit = crit*(1 + ((1+crit^2) + 2*(T-2))/(2*S));
  testGLS(3) = abs(t1) > sccrit; 
end
