%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 estimates a SARAR model with GMM as described in
%Kapoor, Kelejian and Prucha (2007), but updated for unbalanced application
%and allowing a break in the SAR parameter in the primary model. The basic
%form of the model is:
% y = rho1*W*y+rho2*Wd*y+b1*xc+b2*dx+a(i)+u
% u = lam1*W2*u+e

function results = sarar_GMMD(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:
%beta1 - initial estimates from 2sls
%u_SAR2 = [lam1;sig2]; lam1 is SAR term in error SAR model, sig2 is
%         variance of e
%beta2 - 2sls estimate after Cochrane-Orcutt transform
%B2sgmm2 - parameter estimates from GMM
%rho - parameters on W*y and Wd*y
%Bvar - variance of parameters
%tstats - tstats of parameters
%ais - fixed effects
%yhat - predicted y values
%yhat2 - predicted y values assuming rho2 = 0
%BIC - BIC from model
%RSS - residual sum of squares
%corr2 - correlation squared measure
%rsqr - R-squared value
%adj_rsqr - Adjusted R-squared
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%Set up Wdt matrices
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;
    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;
    D0(row_del,:) = [];
    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
if info.model==1
z = Q*z;
else
for i = 1:length(y_unique);
    ry = find(y_spec_re(:,2)==y_unique(i));
    Ty =  y_spec_re(ry(end),1);
    for j = 1:length(ry);
        ty = y_spec_re(ry(j),1);
        if Ty>ty
        zstar(ry(j),:) = sqrt((Ty-ty)/(Ty-ty+1))*(z(ry(j),:)-1/(Ty-ty)*sum(z(ry(j:end),:)));
        else
            zstar(ry(j),:) = sqrt((Ty-ty)/(Ty-ty+1))*(z(ry(j),:)-sum(z(ry(j:end),:)));
        end
    end
    clear ry Ty
end
clear z
z = zstar;
end
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 GMM with weight Omega = eye(3) to get SAR compenents of error
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;
ystar = I_W*y;
wystar = I_W*Wy;
wdystar = I_W*Wdy;
xstar = I_W*x;


%Demean transformed data;
ywith2 = Q*ystar;
wywith2 = Q*wystar;
wdywith2 = Q*wdystar;
xwith2 = Q*xstar;

%Create new instruments
Wx = W*xstar;
WWx = W*W*xstar;
Wdx = Wd*xstar;
z = [xstar 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;
    count = 0;

%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);
A = zeros(length(results.B2sgmm2),1);
A(1,1) = 1; A(end,1) = 1;
A = A*ones(1,2);
A(1,2) = -1;A(end,2) = -1;
results.resid2 = ywith2 - [wywith2 wdywith2 xwith2]*results.B2sgmm2;

%Other Results
N = i;
results.Bvar = sqrt(diag(inv(Xi'*zi/(N*S2)*zi'*Xi)));
results.rho = [results.B2sgmm2(1); results.B2sgmm2(2)];
results.beta = results.B2sgmm2(3:end);
results.ais = (y-[Wy Wdy x]*results.B2sgmm2)-(ywith2-[wywith wdywith xwith]*results.B2sgmm2);
results.yhat = (eye(nobs)-results.rho(1)*W-results.rho(2)*Wd)\(x*results.beta+results.ais);
results.yhat2 = (eye(nobs)-results.rho(1)*W)\(x*results.beta+results.ais);
RSS1 = results.resid2'*results.resid2;
BIC1 = log(RSS1/nobs)+(1+size(x,2)+1)*log(nobs)/nobs;
results.RSS = RSS1;
results.BIC = BIC1;
tstat = results.B2sgmm2./results.Bvar;
results.tstats = tstat;
yme=y-mean(y);
rsqr2=yme'*yme;
rsqr1 = results.resid2'*results.resid2;
results.rsqr=1.0-rsqr1/rsqr2; %rsquared
results.adj_rsqr = 1-(nobs-1)/(nobs-length([results.beta;results.rho]))*...
                   (1-results.rsqr);
ywithhat = (eye(nobs)-results.rho(1)*W-results.rho(2)*Wd)\(xwith2*results.B2sgmm2(3:end));
res1=ywith2-mean(ywith2);
res2=ywithhat-mean(ywith2);
rsq1=res1'*res2;
rsq2=res1'*res1;
rsq3=res2'*res2;
results.corr2=rsq1^2/(rsq2*rsq3); %corr2
