%% OLS for diff-in-diff
% Estimates diff in diff by OLS.  Calculates variance and reports whether
% reject 
%
% $$H_0: \beta=0$$ 
%
% at 5% level in three ways: 
%
% # usual homoskedastic OLS variances
% # clustered on data.t,data.s 
% # clustered on data.s 
% 
% Returns the variance and test results in vectors test and var
%
% 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
%
% See also FGLSDD
function [b test var]=olsDD(data);

  alpha = 0.05; % critical level for tests
  crit = abs(norminv(alpha/2));

  % form X
  X = [data.p data.x data.T data.S];

%% 
% Note that inv should be avoided whenever possible.  
% See
% <http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/
% purpose of inv> regarding the perils of using inv.  Here we need
% inv(X'*X) to compute the variance of beta, but maybe there's a better
% way. 
  iXX = inv(X'*X); 
  if (any(isnan(iXX))) 
    err = 1;
  end
  % compute b
  beta = iXX*X'*data.y;
  b = beta(1);
  % compute residuals
  e = data.y - X*beta;

  % ols test
  n = length(data.y);
  k = length(beta);
  V = e'*e*iXX/(n - k);  
  test(1) = abs(b/sqrt(V(1)))>crit;
  var(1) = V(1);
%% Clustered on state, year
%
% $$ V = (X'X)^{-1} \sum_{s,t} X_{s,t}'\hat{e}_{s,t} \hat{e}_{s,t}'
% X_{s,t} $$
% 
  state=unique(data.s);
  time=unique(data.t);
  V = zeros(size(iXX));
  for s=1:length(state)
    for t=1:length(time)
      g = data.s==state(s) & data.t==time(t);
      V = V + X(g,:)'*e(g)*e(g)'*X(g,:);
    end
  end
  V = iXX*V*iXX;
  test(2) = abs(b/sqrt(V(1)))>crit;
  var(2) = V(1);
%% Clustered on state
%
% $$ V = (X'X)^{-1} \sum_{s} X_{s}'\hat{e}_{s} \hat{e}_{s}'
% X_{s} $$
% 
  V = zeros(size(iXX));
  for s=1:length(state)
    g = data.s==state(s);
    V = V + X(g,:)'*e(g)*e(g)'*X(g,:);
  end
  V = iXX*V*iXX;
  test(3) = abs(b/sqrt(V(1))) > crit;
  var(3) = V(1);
end
