%Harrison Fell
%Resources for the Future
%2011

% This code is utilized in the paper
% Fell, Harrison and Alan Haynie.  2012. "Spatial Competition with
% Changing Market Institutions," Journal of Applied Econometrics, 
% forthcoming.

%Description: This program determines which subset of parameters have a
%break at a known time using the technique of Inuoe and Rossi (2010). The
%base estimation estimates an SARAR model via GMM as described in Kapoor,
%Kelejian and Prucha (2007), but updated for unbalanced application.
%Basic model: y = gamma*W*y+gamma2*Wd*y+b1*xc+b2*dx+a(i)+u
%             u = lam1*W2*u+e

function results = sarar_GMM_Wald(y,xc,dx,Wt,W2,T,i_t,info,b_d)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Inputs
%y - dependent variable
%xc - right-hand side variables
%dx - RHS variables interacted with break dummy
%Wt - NxNxT spatial weight matrices 
%W2 - NxNxT spatial weight matrices for the error term
%i_t - NTx1 vector giving the cross-sectional identifier associated with
%      each y value.
%info - an information structure containing:
%info.Nt - Tx1 vector of number of cross-sectional units at each period.
%info.N - scalar giving total number of cross-sectional units
%info.y_spec - NTx2 matrix with time identifier in first column and i_t
%                 in second column.
%b_d - scalar giving location of break period

%Outputs:
%results.sort_i - ordering of parameters from lowest p-value to highest
%p-value based on individual restrictions from Wald test
%results.break_set - set of parameters that are shown to have a break based
%on Inuoe and Rossi (2010) method.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


x = [xc dx];
Nt = info.Nt;
[a b c] = size(Wt);

%Set up Wdt - weight for the break in the SAR term
Wdt = zeros(a,b,c);
db = zeros(length(y),1);
db(sum(Nt(1:b_d-1))+1:end)=1;
for i=b_d:T;
    Wdt(:,:,i) = Wt(:,:,i);
end

% check size of user inputs for comformability
[nobs nvar] = size(x);
[nchk junk] = size(y);
if nchk ~= nobs
error('sar: wrong size vector y or matrix x');
end

%Form W, Wd, M matrices, each nobs x nobs;
for t = 1:T;
    if t==1;
        zeroMat = zeros(sum(Nt(t+1:T)),Nt(t));
        W = [Wt(1:Nt(t),1:Nt(t),t);zeroMat];
        Wd = [Wdt(1:Nt(t),1:Nt(t),t);zeroMat];
        M = [W2(1:Nt(t),1:Nt(t),t);zeroMat];
    else
        zeroMat_u = zeros(sum(Nt(1:t-1)),Nt(t));
        zeroMat_b = zeros(sum(Nt(t+1:T)),Nt(t));
        Wtmp=[zeroMat_u;Wt(1:Nt(t),1:Nt(t),t);zeroMat_b];
        Wdtmp = [zeroMat_u;Wdt(1:Nt(t),1:Nt(t),t);zeroMat_b];
        Mtmp = [zeroMat_u;W2(1:Nt(t),1:Nt(t),t);zeroMat_b];
        W = cat(2,W,Wtmp); Wd = cat(2,Wd,Wdtmp); M = cat(2,M,Mtmp);
        clear zeroMat_u zeroMat_b Wtmp Wdtmp Mtmp;
    end
end
        
%Form Wdy and Wy
Wy = W*y;
Wdy = Wd*y;

%Transformation of the data by demeaning (within estimator)
y_spec_re = info.y_spec;
y_unique = unique(y_spec_re(:,2));
D = eye(length(y_unique));
for i = 1986:2006;
    rs1 = find(y_spec_re(:,1)==i);
    k = 0;
    row_del = [];
    for j = 1:length(y_unique);
        rs2 = find(y_spec_re(rs1,:)==y_unique(j));
        if isempty(rs2)
            k = k+1;
            row_del(k,1) = j;
        end
        clear rs2
    end
    D0 = D;
    if length(row_del)>=1
    D0(row_del,:) = [];
    end
    if i ==1986;
        D1 = D0;
    else
        D1 = cat(1,D1,D0);
    end
    clear rs1 row_del D0
end
NT = length(D1);
P = D1*inv(D1'*D1)*D1';
Q = eye(NT)-P;
ywith = Q*y;
wywith = Q*Wy;
wdywith = Q*Wdy;
xwith = Q*x;

% Create initial instruments [Wx W*W*x]
Wx = W*x;
WWx = W*W*x;
Wdx = Wd*x;
z = [x Wx];

%Remove repeated columns in z
i = 0; while i < size(z,2)
    i = i+1;
    diffs = sum(repmat(z(:,i),1,size(z,2)-i)-z(:,i+1:end));
    cs = find(abs(diffs)<1e-13);
    if ~isempty(cs)
        z(:,i+cs)=[];
    end
    clear cs diffs
end

%Run initial 2sls estimator to get residuals
z = Q*z;
Pz = z/(z'*z);
Pz = Pz*z';
Z = [wywith wdywith xwith];
Zhat = Pz*Z;
parms_2sls = (Zhat'*Z)\(Zhat'*ywith); 
u = ywith-Z*parms_2sls;
results.beta1 = parms_2sls;

%Run first pass of error SAR model: u = lam1*W2*u+e to get lam1 and var. of
%e.
opt = optimset('fmincon');
opt.Algorithm = 'interior-point';
opt.Display = 'notify';
Omega = eye(3);
[parm1,liktmp,exitflag,output] = fmincon('sarar_GMM_minCon',[.25;25],[],[],[],[],[-1;0],[1;100],[],opt,Omega,u,Q,M);
lam1 = parm1(1);
sig2 = parm1(2);
results.u_SAR1 = [lam1;sig2];

%Run second weighting matrix using weight given in eqn (26) of Kapoor,
%Kelejian and Prucha (2007), but updated for unbalanced application. Also
%there is a mistake in KKP eqn. (26) with the Tw matrix.
Tw11 = 2*trace(Q*Q);
Tw12 = 2*trace(Q*M'*Q*M);
Tw22 = 2*trace(M'*Q*M*M'*Q*M);
Tw23 = 2*trace(M'*Q*M*M'*Q);
Tw33 = 2*trace(M'*Q*M'*Q);
Tw = [Tw11 Tw12 0
      Tw12 Tw22 Tw23
      0    Tw23 Tw33];
Mvar = kron(sig2^2,Tw);
Omega = inv(Mvar);
[parm2,liktmp,exitflag,output] = fmincon('sarar_GMM_minCon',parm1,[],[],[],[],[-1;0],[1;100],[],opt,Omega,u,Q,M);
lam1 = (parm2(1));
sig2 = (parm2(2));
results.u_SAR2 = [lam1;sig2];

%Perform Cocharan-Orcutt-type transformation
I_W = eye(nobs)-lam1*M;
y = I_W*y;
Wy = I_W*Wy;
Wdy = I_W*Wdy;
x = I_W*x;
ywith2 = Q*y;
wywith2 = Q*Wy;
wdywith2 = Q*Wdy;
xwith2 = Q*x;

%Create new instruments
Wx = W*x;
WWx = W*W*x;
Wdx = Wd*x;
z = [x Wx];

%Remove repeated columns in z
i = 0; while i < size(z,2)
    i = i+1;
    diffs = sum(repmat(z(:,i),1,size(z,2)-i)-z(:,i+1:end));
    cs = find(abs(diffs)<1e-13);
    if ~isempty(cs)
        z(:,i+cs)=[];
    end
    clear cs diffs
end%}

%Run 2sls on transformed data
z = Q*z;
Pz = z/(z'*z);
Pz = Pz*z';
Z2 = [wywith2  wdywith2 xwith2];
Zhat = Pz*Z2;
parms_2sls = (Zhat'*Z2)\(Zhat'*ywith2); 
results.beta2 = parms_2sls;
u2 = ywith2-Z2*parms_2sls;
    
%Run GMM on tranformed variables
K = size(z,2);
S = zeros(K,K);
S2 =0;
Qxz = 0;
Qzy = 0;

for i = 1:length(y_unique)
    rs = find(i_t==y_unique(i));
    Ti(i,1) = length(rs);
    count = count+Ti(i,1);
    yi(count-Ti(i,1)+1:count,1) = ywith2(rs);
    wyi(count-Ti(i,1)+1:count,1) = wywith2(rs);
    wdyi(count-Ti(i,1)+1:count,1) = wdywith2(rs);
    xi(count-Ti(i,1)+1:count,:) = xwith2(rs,:);
    zi(count-Ti(i,1)+1:count,:) = z(rs,:);
    ei(count-Ti(i,1)+1:count,1) = u2(rs);
    zii = z(rs,:);
    eii = u2(rs);
    xii = xwith2(rs,:);
    wyii = wywith2(rs,:);
    wdyii = wdywith2(rs,:);
    yii = ywith2(rs,:);
    Vari(i,1) = eii'*eii;
    S2 = S2+(zii'*eii*eii'*zii);
    Qxz = Qxz + [wyii wdyii xii]'*zii;
    Qzy = Qzy + zii'*yii;
    clear zii eii rs xii wyii wdyii yii
end
S2 = 1/i*S2;
Xi = [wyi wdyi xi];

%Cameron and Trivedi style GMM
results.B2sgmm2 = (Xi'*zi/S2*zi'*Xi)\(Xi'*zi/S2*zi'*yi);
results.resid2 = ywith - [wywith  wdywith xwith]*results.B2sgmm2;
resid2 = yi-Xi*results.B2sgmm2;
mbar = 1/nobs*zi'*resid2;
nq = sqrt(nobs)*mbar'/S2*sqrt(nobs)*mbar;
results.nq = nq;
results.mbar = mbar;
results.S = S2;
N= i;

%Now test individual restriction
K1 = size(xc,2)+1;
b0 = results.B2sgmm2;
r0 = (K1+1:1:((K1-1)*2+2))';
Q0 = inv(Xi'*zi/(N*S2)*zi'*Xi);
A0 = zeros(length(r0),size(b0,1));
A0(:,K1+1:end) = eye(length(r0));
W0 = (b0(r0)'/(A0*Q0*A0')*b0(r0));
p_prob0 = 1-chis_prb(W0,length(r0));

if p_prob0<0.10||p_prob0>0.900
yme = y-mean(y);
rsqu = 1-resid2'*resid2/(yme'*yme);
for i = 1:K1
    Ar1 = zeros(1,size(b0,1));
    if i ==1
        br1 = b0(2);
        AR1(2) = 1;
    else
        br1 = b0(K1+i);
        Ar1(K1+i) = 1;
    end
    W1(i,1) =(br1'/(Ar1*Q0*Ar1')*br1);
    p_prob(i,1) = 1-chis_prb(W1(i,1),1);
    clear Ar1 
end
[p_sort sort_i] = sort(p_prob);
results.sort_i = sort_i;

for i =1:K1;
    const_set = setdiff((1:K1),sort_i(1:i));
    break_set = sort_i(1:i);
    br2 = b0(const_set+K1);
    Ar2 = zeros(length(const_set),size(b0,1));
    if sum(sort_i(1:i)==1)==1
    for k = 1:length(const_set);
        Ar2(k,const_set(k)+K1)=1;
    end
    else
        Ar2(1,2) = 1;
        const_set2 = const_set(const_set~=1);
        for k = 1:length(const_set2)
            Ar2(k+1,const_set2(k)+K1)=1;
        end
    end
    W2 = (br2'/(Ar2*Q0*Ar2')*br2);
    p_prob2(i,1) = 1-chis_prb(W2,length(const_set));
    if p_prob2(i,1)>0.01 
        break 
    end
    clear const_set br2 Ar2 const_set2
end
if sum(break_set==1)==1
    break_set(break_set==1)=K1+1;
end
break_set = break_set-1;
else
        break_set = 0;
        results.sort_i = zeros(K1,1);
end
results.break_set = break_set;
    
    
