clear all; close all;

% choose whether to use state-by-year fixed effects;
p.stateByYear = true;

% load data;
d = readCSVfileIntoStruct('yieldData.csv');
% sort data;
structNames = fieldnames(d);
[dump,sortIndex] = sortrows([d.fips d.year]);
for s = 1:max(size(structNames))
    eval(['d.' char(structNames(s)) '= d.' char(structNames(s)) '(sortIndex,:);'])
end
groupID = d.fips;
n0(:,1) = find(groupID ~= [groupID(1)-1; groupID(1:size(groupID,1)-1)]);
n1(:,1) = find(groupID ~= [groupID(2:size(groupID,1)); groupID(size(groupID,1),1)+1]);
clear structNames s dump sortIndex groupID;

% average profit in data;
yearUnique = unique(d.year);
for t = 1:max(size(yearUnique))
    hilf(t,1) = sum(d.cropArea(d.year == yearUnique(t)));
end
totProd = mean(hilf)*sum(d.cropProd)/sum(d.cropArea);
clear hilf t;

% weighting variable;
d.weight = sqrt(d.cropArea);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timeControl = [];
if p.stateByYear
    stateUnique = unique(floor(d.fips/1000));
    for s = 1:max(size(stateUnique))
        % need at least two counties in state
        if max(size(unique(d.fips(floor(d.fips/1000)==stateUnique(s))))) > 1 
            yearUniqueState = unique(d.year(floor(d.fips/1000)==stateUnique(s)));
            for t = 2:max(size(yearUniqueState))
                timeControl = [timeControl ((floor(d.fips/1000) == stateUnique(s)) & (d.year == yearUniqueState(t)))];
            end
        end
    end
    clear stateUnique s t;
else
    for t = 2:max(size(yearUnique))
        timeControl = [timeControl (d.year == yearUnique(t))];
    end
    clear t;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% column a;
y = d.y;
X = [d.dry_dd89 d.dry_dd89_sq d.dry_prcp d.dry_prcp_sq d.irr_dd89 d.irr_dd89_sq d.irr_prcp d.irr_prcp_sq d.dry d.x2 d.x3 d.x4 d.x5 d.x6 d.x7 d.x8 d.x9 timeControl];
% demean by county
for i = 1:size(n0,1)
    wf = d.weight(n0(i):n1(i))'/sum(d.weight(n0(i):n1(i)));
    y(n0(i):n1(i),1) = y(n0(i):n1(i))   - wf*y(n0(i):n1(i));
    X(n0(i):n1(i),:) = X(n0(i):n1(i),:) - repmat(wf*X(n0(i):n1(i),:),n1(i)-n0(i)+1,1);
end
y = y.*d.weight;
X = X.*repmat(d.weight,1,size(X,2));
[b,V] = conley(y,X,[d.latitude d.longitude],d.year,[5 5]);
clear i;

% get average by fips;
for i = 1:size(n0,1)
    h.m_dry_dd89(n0(i):n1(i),1) = mean(d.dry_dd89(n0(i):n1(i)));
    h.m_dry_prcp(n0(i):n1(i),1) = mean(d.dry_prcp(n0(i):n1(i)));
    h.m_irr_dd89(n0(i):n1(i),1) = mean(d.irr_dd89(n0(i):n1(i)));
    h.m_irr_prcp(n0(i):n1(i),1) = mean(d.irr_prcp(n0(i):n1(i)));
end
clear i;

hilf(:,1) = (d.dry == 1).*(d.dd89_h2_long - d.dd89_7000);
hilf(:,2) = 2*hilf(:,1).*h.m_dry_dd89 + hilf(:,1).*hilf(:,1);
hilf(:,3) = (d.dry == 1).*(d.prcp_h2_long - d.prcp_7000);
hilf(:,4) = 2*hilf(:,3).*h.m_dry_prcp + hilf(:,3).*hilf(:,3);
hilf(:,5) = (d.dry == 0).*(d.dd89_h2_long - d.dd89_7000);
hilf(:,6) = 2*hilf(:,5).*h.m_irr_dd89 + hilf(:,5).*hilf(:,5);
hilf(:,7) = (d.dry == 0).*(d.prcp_h2_long - d.prcp_7000);
hilf(:,8) = 2*hilf(:,7).*h.m_irr_prcp + hilf(:,7).*hilf(:,7);

for t = 1:max(size(yearUnique))
    Xdiff(t,:) = d.cropArea(d.year == yearUnique(t))'*hilf(d.year == yearUnique(t),:)/totProd;
end
clear t;
Xdiff = mean(Xdiff,1);
disp(['Climate Impacts in column a1: ' sprintf('%6.2f',Xdiff*b(1:8)*100)])
disp(['    Conley standard error is: ' sprintf('%6.2f',sqrt(Xdiff*V(1:8,1:8)*Xdiff')*100)])
disp(' ');
clear h hilf y X wf Xdiff V b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% column b1;
y = d.y;
X = [d.dry_dd89dm d.dry_dd89dm_sq d.dry_prec d.dry_prec_sq d.irr_dd89dm d.irr_dd89dm_sq d.irr_prec d.irr_prec_sq d.dry d.x2 d.x3 d.x4 d.x5 d.x6 d.x7 d.x8 d.x9 timeControl];
% demean by county
for i = 1:size(n0,1)
    wf = d.weight(n0(i):n1(i))'/sum(d.weight(n0(i):n1(i)));
    y(n0(i):n1(i),1) = y(n0(i):n1(i))   - wf*y(n0(i):n1(i));
    X(n0(i):n1(i),:) = X(n0(i):n1(i),:) - repmat(wf*X(n0(i):n1(i),:),n1(i)-n0(i)+1,1);
end
y = y.*d.weight;
X = X.*repmat(d.weight,1,size(X,2));
[b,V] = conley(y,X,[d.latitude d.longitude],d.year,[5 5]);
clear i;

h(:,1) = (d.dry == 1).*(d.d_dday8_32dmH2);
h(:,2) = (d.dry == 1).*(d.d_dday8_32dmH2_sq);
h(:,3) = (d.dry == 1).*(d.d_precH2);
h(:,4) = (d.dry == 1).*(d.d_precH2_sq);
h(:,5) = (d.dry == 0).*(d.d_dday8_32dmH2);
h(:,6) = (d.dry == 0).*(d.d_dday8_32dmH2_sq);
h(:,7) = (d.dry == 0).*(d.d_precH2);
h(:,8) = (d.dry == 0).*(d.d_precH2_sq);

for t = 1:max(size(yearUnique))
    Xdiff(t,:) = d.cropArea(d.year == yearUnique(t))'*h(d.year == yearUnique(t),:)/totProd;
end
clear t;
Xdiff = mean(Xdiff,1);
disp(['Climate Impacts in column b1: ' sprintf('%6.2f',Xdiff*b(1:8)*100)])
disp(['    Conley standard error is: ' sprintf('%6.2f',sqrt(Xdiff*V(1:8,1:8)*Xdiff')*100)])
disp(' ');
clear h Xdiff wf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% column b2;
h(:,1) = (d.dry == 1).*(d.d_dday8_32dm);
h(:,2) = (d.dry == 1).*(d.d_dday8_32dm_sq);
h(:,3) = (d.dry == 1).*(d.d_prec);
h(:,4) = (d.dry == 1).*(d.d_prec_sq);
h(:,5) = (d.dry == 0).*(d.d_dday8_32dm);
h(:,6) = (d.dry == 0).*(d.d_dday8_32dm_sq);
h(:,7) = (d.dry == 0).*(d.d_prec);
h(:,8) = (d.dry == 0).*(d.d_prec_sq);

for t = 1:max(size(yearUnique))
    Xdiff(t,:) = d.cropArea(d.year == yearUnique(t))'*h(d.year == yearUnique(t),:)/totProd;
end
clear t;
Xdiff = mean(Xdiff,1);
disp(['Climate Impacts in column b2: ' sprintf('%6.2f',Xdiff*b(1:8)*100)])
disp(['    Conley standard error is: ' sprintf('%6.2f',sqrt(Xdiff*V(1:8,1:8)*Xdiff')*100)])
disp(' ');
clear h y X Xdiff V b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% column c1;
y = d.y;
X = [d.dry_good d.dry_good_sq d.dry_bad d.dry_prec d.dry_prec_sq d.irr_good d.irr_good_sq d.irr_bad d.irr_prec d.irr_prec_sq d.dry d.x2 d.x3 d.x4 d.x5 d.x6 d.x7 d.x8 d.x9 timeControl];
% demean by county
for i = 1:size(n0,1)
    wf = d.weight(n0(i):n1(i))'/sum(d.weight(n0(i):n1(i)));
    y(n0(i):n1(i),1) = y(n0(i):n1(i))   - wf*y(n0(i):n1(i));
    X(n0(i):n1(i),:) = X(n0(i):n1(i),:) - repmat(wf*X(n0(i):n1(i),:),n1(i)-n0(i)+1,1);
end
y = y.*d.weight;
X = X.*repmat(d.weight,1,size(X,2));
[b,V] = conley(y,X,[d.latitude d.longitude],d.year,[5 5]);
clear i;

h(:, 1) = (d.dry == 1).*(d.d_ddayGoodH2);
h(:, 2) = (d.dry == 1).*(d.d_ddayGoodH2_sq);
h(:, 3) = (d.dry == 1).*(d.d_ddayBadH2);
h(:, 4) = (d.dry == 1).*(d.d_precH2);
h(:, 5) = (d.dry == 1).*(d.d_precH2_sq);
h(:, 6) = (d.dry == 0).*(d.d_ddayGoodH2);
h(:, 7) = (d.dry == 0).*(d.d_ddayGoodH2_sq);
h(:, 8) = (d.dry == 0).*(d.d_ddayBadH2);
h(:, 9) = (d.dry == 0).*(d.d_precH2);
h(:,10) = (d.dry == 0).*(d.d_precH2_sq);

for t = 1:max(size(yearUnique))
    Xdiff(t,:) = d.cropArea(d.year == yearUnique(t))'*h(d.year == yearUnique(t),:)/totProd;
end
clear t;
Xdiff = mean(Xdiff,1);
disp(['Climate Impacts in column c1: ' sprintf('%6.2f',Xdiff*b(1:10)*100)])
disp(['    Conley standard error is: ' sprintf('%6.2f',sqrt(Xdiff*V(1:10,1:10)*Xdiff')*100)])
disp(' ');
clear h Xdiff wf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% column c2;
h(:, 1) = (d.dry == 1).*(d.d_ddayGood);
h(:, 2) = (d.dry == 1).*(d.d_ddayGood_sq);
h(:, 3) = (d.dry == 1).*(d.d_ddayBad);
h(:, 4) = (d.dry == 1).*(d.d_prec);
h(:, 5) = (d.dry == 1).*(d.d_prec_sq);
h(:, 6) = (d.dry == 0).*(d.d_ddayGood);
h(:, 7) = (d.dry == 0).*(d.d_ddayGood_sq);
h(:, 8) = (d.dry == 0).*(d.d_ddayBad);
h(:, 9) = (d.dry == 0).*(d.d_prec);
h(:,10) = (d.dry == 0).*(d.d_prec_sq);

for t = 1:max(size(yearUnique))
    Xdiff(t,:) = d.cropArea(d.year == yearUnique(t))'*h(d.year == yearUnique(t),:)/totProd;
end
clear t;
Xdiff = mean(Xdiff,1);
disp(['Climate Impacts in column c2: ' sprintf('%6.2f',Xdiff*b(1:10)*100)])
disp(['    Conley standard error is: ' sprintf('%6.2f',sqrt(Xdiff*V(1:10,1:10)*Xdiff')*100)])
disp(' ');
clear h y X Xdiff V b;