# delimit ;
clear;
set memory 1g;
set more off;
version 10;
capture log close;
log using yieldRegression.log, replace;

* ###############################################################################################;
* comparison of how various weather variables predict yields and profits (Tables 1, 2, and A1);
* choose what crop to use: for soybeans use "soyb", for corn use "corn", for profit use "prof");
local cropChoice = "soyb";
* choose what lower/upper degree day bound to use;
local ddayGoodLB 8;
local ddayGoodUB 32;
local ddayBadLB  34sqrt;
* chose temporal control (year fixed effects: i.year / state-by-year fixed effects: sst*);
local timeControl sst*;
* ###############################################################################################;

*************************************************************************************************;
* Part 1: build data set;
*************************************************************************************************;
* load DG data set;
use DATA1, clear;
drop sst* statfips state;
* merge in other weather data sets we replicate;
sort fips year;
merge fips year using dataPanelNew;
drop _merge;

* Create variables (following DG's script);
* get profits into code format of yields;
gen prof_planted = fland;
gen prof_prod = (rtfsale-rtotexp)*1000;

* yields;
gen cropArea = `cropChoice'_planted;
gen cropProd = `cropChoice'_prod;
gen y = cropProd/cropArea;
* keep data set consistent;
drop if (cropProd==.) | (cropArea==.) | (cropArea==0) | (dd89==.) | (dday`ddayGoodLB'_`ddayGoodUB' == .);

* created squared terms;
gen dd89_sq = dd89^2;
gen prcp_sq = prcp^2;

* create coefficients that differ by dryland / irrigated;
* DG data;
gen dry_dd89=    (dry==1)*dd89;
gen dry_dd89_sq= (dry==1)*dd89*dd89;
gen dry_prcp=    (dry==1)*prcp;
gen dry_prcp_sq= (dry==1)*prcp*prcp;
gen irr_dd89=    (dry==0)*dd89;
gen irr_dd89_sq= (dry==0)*dd89*dd89;
gen irr_prcp=    (dry==0)*prcp;
gen irr_prcp_sq= (dry==0)*prcp*prcp;

* our replication of DG data using average daily temperature;
gen dry_dd89dm=   (dry==1)*dday8_32dm;
gen dry_dd89dm_sq=(dry==1)*dday8_32dm*dday8_32dm;
gen dry_prec =    (dry==1)*prec;
gen dry_prec_sq=  (dry==1)*prec*prec;
gen irr_dd89dm=   (dry==0)*dday8_32dm;
gen irr_dd89dm_sq=(dry==0)*dday8_32dm*dday8_32dm;
gen irr_prec =    (dry==0)*prec;
gen irr_prec_sq=  (dry==0)*prec*prec;

* weather data set using temperature distribution between tMin and tMax;
gen dry_good =    (dry==1)*dday`ddayGoodLB'_`ddayGoodUB';
gen dry_good_sq = (dry==1)*dday`ddayGoodLB'_`ddayGoodUB'*dday`ddayGoodLB'_`ddayGoodUB';
gen dry_bad =     (dry==1)*dday`ddayBadLB';
gen irr_good =    (dry==0)*dday`ddayGoodLB'_`ddayGoodUB';
gen irr_good_sq = (dry==0)*dday`ddayGoodLB'_`ddayGoodUB'*dday`ddayGoodLB'_`ddayGoodUB';
gen irr_bad =     (dry==0)*dday`ddayBadLB';

* changes in weather variables;
gen d_ddayGoodH2    = d_dday`ddayGoodLB'_`ddayGoodUB'H2;
gen d_ddayGoodH2_sq = d_dday`ddayGoodLB'_`ddayGoodUB'H2_sq;
gen d_ddayBadH2     = d_dday`ddayBadLB'H2;
gen d_ddayGood      = d_dday`ddayGoodLB'_`ddayGoodUB';
gen d_ddayGood_sq   = d_dday`ddayGoodLB'_`ddayGoodUB'_sq;
gen d_ddayBad       = d_dday`ddayBadLB';

* save data for matlab (get Conly's variance-covariance matrix adjusted for spatial correlation);
outsheet using yieldData.csv, comma replace;

* create state-by-year fixed effects;
qui levelsof year, local(yearList);
qui summ year; local yearMin = r(min);
gen stateFIPS = floor(fips/1000);
qui levelsof stateFIPS, local(stateFIPSlist);
foreach s of local stateFIPSlist {;
    foreach t of local yearList {;
        if (`t' > `yearMin') {;
            gen sst_`s'_`t' = (year == `t') & (stateFIPS == `s');
        };
    };
};


*************************************************************************************************;
* Part 2: regression results;
*************************************************************************************************;
* regression with no weather variable;
xi: areg y x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;

* baseline regression: Column 1a of Table 8 in DG;
xi: areg y dry_dd89 dry_dd89_sq dry_prcp dry_prcp_sq irr_dd89 irr_dd89_sq irr_prcp irr_prcp_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
predict yield_DG, xb;

* regression with our reconstructed variable of DG;
xi: areg y dry_dd89dm dry_dd89dm_sq dry_prec dry_prec_sq irr_dd89dm irr_dd89dm_sq irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
predict yield_dm, xb;

* regression using daily temperature range between tMin-tMax (specification REStat 2006);
xi: areg y dry_good dry_good_sq dry_bad dry_prec dry_prec_sq irr_good irr_good_sq irr_bad irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
predict yield_minMax, xb;

*************************************************************************************************;
* Part 3: non-nested J-test;
*************************************************************************************************;
* testing DG model against other models;
* part a: against our construction using daily mean;
xi: areg y yield_dm dry_dd89 dry_dd89_sq dry_prcp dry_prcp_sq irr_dd89 irr_dd89_sq irr_prcp irr_prcp_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* part b: against weather data using daily tMin tMax;
xi: areg y yield_minMax dry_dd89 dry_dd89_sq dry_prcp dry_prcp_sq irr_dd89 irr_dd89_sq irr_prcp irr_prcp_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;

* testing our models against DG;
* part a: our construction using daily mean;
xi: areg y yield_DG dry_dd89dm dry_dd89dm_sq dry_prec dry_prec_sq irr_dd89dm irr_dd89dm_sq irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* part b: our weather using minimum and maximum temperature;
xi: areg y yield_DG dry_good dry_good_sq dry_bad dry_prec dry_prec_sq irr_good irr_good_sq irr_bad irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;


*************************************************************************************************;
* Part 4: climate change impacts;
*************************************************************************************************;
* average historic crop production area * yield;
* note: this is the normalizing constant to get percent impacts;
qui egen hilf1 = total(cropProd);
qui egen hilf2 = total(cropArea);
qui egen hilf = total(cropArea), by(year);
qui summ hilf; local totProd = r(mean)*hilf1[1]/hilf2[1]/100;
drop hilf hilf1 hilf2;

*************************************************************************************************;
* column a: baseline regression: Column 1a of Table 8 in DG;
* derive average climate in sample;
egen m_dry_dd89=mean(dry_dd89), by(fips);
egen m_dry_prcp=mean(dry_prcp), by(fips);
egen m_irr_dd89=mean(irr_dd89), by(fips);
egen m_irr_prcp=mean(irr_prcp), by(fips);

* derive change in climate variables under global warming;
gen d_dry_dd89=    (dry==1)*(dd89_h2_long - dd89_7000);
gen d_dry_dd89_sq= 2*d_dry_dd89*m_dry_dd89 + d_dry_dd89*d_dry_dd89;
gen d_dry_prcp=    (dry==1)*(prcp_h2_long - prcp_7000);
gen d_dry_prcp_sq= 2*d_dry_prcp*m_dry_prcp + d_dry_prcp*d_dry_prcp;
gen d_irr_dd89=    (dry==0)*(dd89_h2_long - dd89_7000);
gen d_irr_dd89_sq= 2*d_irr_dd89*m_irr_dd89 + d_irr_dd89*d_irr_dd89;
gen d_irr_prcp=    (dry==0)*(prcp_h2_long - prcp_7000);
gen d_irr_prcp_sq= 2*d_irr_prcp*m_irr_prcp + d_irr_prcp*d_irr_prcp;

local irrigList dry irr;
local varList dd89 dd89_sq prcp prcp_sq;
foreach i of local irrigList {;
    foreach v of local varList {;
        * dryland variables;
        qui egen hilf = total(d_`i'_`v'*cropArea/`totProd'), by(year);
        qui summ hilf; local diff_`i'_`v' = r(mean); drop hilf;
    };
};

* use regression coefficients to derive climate change impacts;
xi: areg y dry_dd89 dry_dd89_sq dry_prcp dry_prcp_sq irr_dd89 irr_dd89_sq irr_prcp irr_prcp_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* derive climate impact (change in area-weighted average yield/acre);
lincom _b[dry_dd89] * `diff_dry_dd89' + _b[dry_dd89_sq] * `diff_dry_dd89_sq' + _b[dry_prcp] * `diff_dry_prcp' + _b[dry_prcp_sq] * `diff_dry_prcp_sq'
     + _b[irr_dd89] * `diff_irr_dd89' + _b[irr_dd89_sq] * `diff_irr_dd89_sq' + _b[irr_prcp] * `diff_irr_prcp' + _b[irr_prcp_sq] * `diff_irr_prcp_sq' ;


*************************************************************************************************;
* column b1: regression with our reconstructed variable of DG - Hadley II;
* derive change in climate variables under global warming;
gen d_dry_dd89dmH2    = (dry==1)*d_dday8_32dmH2;
gen d_dry_dd89dmH2_sq = (dry==1)*d_dday8_32dmH2_sq;
gen d_dry_precH2      = (dry==1)*d_precH2;
gen d_dry_precH2_sq   = (dry==1)*d_precH2_sq;
gen d_irr_dd89dmH2    = (dry==0)*d_dday8_32dmH2;
gen d_irr_dd89dmH2_sq = (dry==0)*d_dday8_32dmH2_sq;
gen d_irr_precH2      = (dry==0)*d_precH2;
gen d_irr_precH2_sq   = (dry==0)*d_precH2_sq;

local varList dd89dmH2 dd89dmH2_sq precH2 precH2_sq;
foreach i of local irrigList {;
    foreach v of local varList {;
        * dryland variables;
        qui egen hilf = total(d_`i'_`v'*cropArea/`totProd'), by(year);
        qui summ hilf; local diff_`i'_`v' = r(mean); drop hilf;
    };
};

* use regression coefficients to derive climate change impacts;
xi: areg y dry_dd89dm dry_dd89dm_sq dry_prec dry_prec_sq irr_dd89dm irr_dd89dm_sq irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* derive climate impact (change in area-weighted average yield/acre);
lincom _b[dry_dd89dm] * `diff_dry_dd89dmH2' + _b[dry_dd89dm_sq] * `diff_dry_dd89dmH2_sq' + _b[dry_prec] * `diff_dry_precH2' + _b[dry_prec_sq] * `diff_dry_precH2_sq'
     + _b[irr_dd89dm] * `diff_irr_dd89dmH2' + _b[irr_dd89dm_sq] * `diff_irr_dd89dmH2_sq' + _b[irr_prec] * `diff_irr_precH2' + _b[irr_prec_sq] * `diff_irr_precH2_sq' ;


*************************************************************************************************;
* column b2: regression with our reconstructed variable of DG - Hadley III;
* derive change in climate variables under global warming;
gen d_dry_dd89dm    = (dry==1)*d_dday8_32dm;
gen d_dry_dd89dm_sq = (dry==1)*d_dday8_32dm_sq;
gen d_dry_prec      = (dry==1)*d_prec;
gen d_dry_prec_sq   = (dry==1)*d_prec_sq;
gen d_irr_dd89dm    = (dry==0)*d_dday8_32dm;
gen d_irr_dd89dm_sq = (dry==0)*d_dday8_32dm_sq;
gen d_irr_prec      = (dry==0)*d_prec;
gen d_irr_prec_sq   = (dry==0)*d_prec_sq;

local varList dd89dm dd89dm_sq prec prec_sq;
foreach i of local irrigList {;
    foreach v of local varList {;
        * dryland variables;
        qui egen hilf = total(d_`i'_`v'*cropArea/`totProd'), by(year);
        qui summ hilf; local diff_`i'_`v' = r(mean); drop hilf;
    };
};

* use regression coefficients to derive climate change impacts;
xi: areg y dry_dd89dm dry_dd89dm_sq dry_prec dry_prec_sq irr_dd89dm irr_dd89dm_sq irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* derive climate impact (change in area-weighted average yield/acre);
lincom _b[dry_dd89dm] * `diff_dry_dd89dm' + _b[dry_dd89dm_sq] * `diff_dry_dd89dm_sq' + _b[dry_prec] * `diff_dry_prec' + _b[dry_prec_sq] * `diff_dry_prec_sq'
     + _b[irr_dd89dm] * `diff_irr_dd89dm' + _b[irr_dd89dm_sq] * `diff_irr_dd89dm_sq' + _b[irr_prec] * `diff_irr_prec' + _b[irr_prec_sq] * `diff_irr_prec_sq' ;


*************************************************************************************************;
* column c1: regression using daily temperature range between tMin-tMax (specification REStat 2006) - Hadley II;
* derive change in climate variables under global warming;
gen d_dry_goodH2    = (dry==1)*d_ddayGoodH2;
gen d_dry_goodH2_sq = (dry==1)*d_ddayGoodH2_sq;
gen d_dry_badH2     = (dry==1)*d_ddayBadH2;
gen d_irr_goodH2    = (dry==0)*d_ddayGoodH2;
gen d_irr_goodH2_sq = (dry==0)*d_ddayGoodH2_sq;
gen d_irr_badH2     = (dry==0)*d_ddayBadH2;

local varList goodH2 goodH2_sq badH2 precH2 precH2_sq;
foreach i of local irrigList {;
    foreach v of local varList {;
        * dryland variables;
        qui egen hilf = total(d_`i'_`v'*cropArea/`totProd'), by(year);
        qui summ hilf; local diff_`i'_`v' = r(mean); drop hilf;
    };
};

xi: areg y dry_good dry_good_sq dry_bad dry_prec dry_prec_sq irr_good irr_good_sq irr_bad irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* derive climate impact (change in area-weighted average yield/acre);
lincom _b[dry_good] * `diff_dry_goodH2' + _b[dry_good_sq] * `diff_dry_goodH2_sq' + _b[dry_bad] * `diff_dry_badH2' + _b[dry_prec] * `diff_dry_precH2' + _b[dry_prec_sq] * `diff_dry_precH2_sq'
     + _b[irr_good] * `diff_irr_goodH2' + _b[irr_good_sq] * `diff_irr_goodH2_sq' + _b[irr_bad] * `diff_irr_badH2' + _b[irr_prec] * `diff_irr_precH2' + _b[irr_prec_sq] * `diff_irr_precH2_sq' ;


*************************************************************************************************;
* column c2: regression using daily temperature range between tMin-tMax (specification REStat 2006) - Hadley III;
* derive change in climate variables under global warming;
gen d_dry_good    = (dry==1)*d_ddayGood;
gen d_dry_good_sq = (dry==1)*d_ddayGood_sq;
gen d_dry_bad     = (dry==1)*d_ddayBad;
gen d_irr_good    = (dry==0)*d_ddayGood;
gen d_irr_good_sq = (dry==0)*d_ddayGood_sq;
gen d_irr_bad     = (dry==0)*d_ddayBad;

local varList good good_sq bad prec prec_sq;
foreach i of local irrigList {;
    foreach v of local varList {;
        * dryland variables;
        qui egen hilf = total(d_`i'_`v'*cropArea/`totProd'), by(year);
        qui summ hilf; local diff_`i'_`v' = r(mean); drop hilf;
    };
};

* use regression coefficients to derive climate change impacts;
xi: areg y dry_good dry_good_sq dry_bad dry_prec dry_prec_sq irr_good irr_good_sq irr_bad irr_prec irr_prec_sq
           dry x2-x9 `timeControl' [aweight=cropArea], a(fips) robust;
* derive climate impact (change in area-weighted average yield/acre);
lincom _b[dry_good] * `diff_dry_good' + _b[dry_good_sq] * `diff_dry_good_sq' + _b[dry_bad] * `diff_dry_bad' + _b[dry_prec] * `diff_dry_prec' + _b[dry_prec_sq] * `diff_dry_prec_sq'
     + _b[irr_good] * `diff_irr_good' + _b[irr_good_sq] * `diff_irr_good_sq' + _b[irr_bad] * `diff_irr_bad' + _b[irr_prec] * `diff_irr_prec' + _b[irr_prec_sq] * `diff_irr_prec_sq' ;

log close;