%% Simulate a diff in diff data set
% Draws a sample of N states with T time periods from original data orig.  
% assigns half the states to be treated with a policy effect of gamma.
% 
% The data structure is meant to be a
% generic diff in diff setup with 2 level clustering.  The fields of data
% are:
%  
% * 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
%
% Input arguments are:
%
% * orig = original data structure
% * N = number of states to draw in sample
% * T = number of periods to draw in sample
% * gamma = effect of policy
%
% States are drawn with replacement.  Time is drawn as a continuous block
% of length T.  The policy is randomly chosen to affect half of the
% states beginning at t ~ U(0.2*T,0.8*T).  When the policy is in effect,
% y is changed to (1+gamma)*y
%
% See also OLSDD, FGLSDD, AGGDATA
function simd = simulate(orig,N,T,gamma)
  states = unique(orig.s);
  times = unique(orig.t);
  
  if ~isequal((1:length(times))',sort(times))
    error('simulate requires t go from 1 to T in the original data');
  end

  % ks = indices of states to keep in sample
  ks = ceil(rand(N,1)*length(states));
  % kt = indices of times to keep in sample
  kt = ceil(rand()*(length(times)-T));
  kt = kt:(kt+T-1);
  % tuse = 1 if orig.t = any(kt), else 0
  tuse = zeros(size(orig.y)); 
  for t=1:length(kt)
    tuse = tuse | orig.t==kt(t);
  end
  use = []; % = 1 if want observation
  ss = [];  % simulated state index, because we're sampling with
            % replacement, ks might not consist of unique values, ss will 
  for s=1:length(ks)
    f = find(orig.s==states(ks(s)) & tuse);
    % notice that use and ss are growing within this loop.  this
    % adversely affects Matlab's performance 
    use = [use; f]; %#ok (tells mlint to ignore)
    ss = [ss; s*ones(size(f))]; %#ok
  end
  
  % create the sampled data
  fname = fields(orig); 
  for f=1:length(fname)
    simd.(fname{f}) = orig.(fname{f})(use);
  end
  % make sure sample state indices are unique and go from 1 to N
  simd.s = ss;
  assert(isequal((1:N)',sort(unique(simd.s))));
  % make time start from 1
  simd.t = simd.t - min(simd.t) + 1;

  % generate the policy change
  t1 = min(simd.t);
  t2 = max(simd.t);
  t = round(T*0.2 + rand()*T*0.6);
  if (t>=t2) 
    t = t2-1;
  elseif (t<=t1) 
    t = t1+1;
  end
  
  % pick which states are affected
  p = rand(N,1);
  mp = median(p);
  paff = zeros(size(simd.s));
  for s=1:N
    if (p(s)>mp) 
      paff(simd.s==s) = 1;
    end
  end
  
  % assign policy indicator
  simd.p(simd.t<t) = 0;
  simd.p(simd.t>=t & paff) = 1;
  % change y when policy in effect
  simd.y = simd.y + gamma*simd.p;
  if (max(simd.p)==0) 
    error('simd.p = 0 everywhere');
  end
  
  % create time and state dummies
  u = unique(simd.t);
  simd.T = zeros(length(simd.t),length(u)-1);
  for i=2:length(u) % omit first time dummy
    simd.T(:,i-1) = simd.t==u(i);
  end
  u = unique(simd.s);
  simd.S = zeros(length(simd.s),length(u));
  for i=1:length(u)
    simd.S(:,i) = simd.s==u(i);
  end
  clear u;

end
