# delimit ;
clear;
set memory 1g;
set more off;
version 10;

capture log close;
log using countyCommoditySales.log, replace;

local commodityListAll corn wheat oats barley sorghum soybeans cotton tobacco;
local commodityList corn wheat soybeans cotton;

* ###############################################################################################;
* script looks at county-level storage;
* regresses sales on value of production;
* ###############################################################################################;

*************************************************************************************************;
* Part 1: build data set;
*************************************************************************************************;
use dataCommoditySales;

* get values on per acre basis;
foreach c of local commodityList {;
    gen `c'ProdValueAcres = `c'ProdValue/`c'Acres;
    gen `c'SalesAcres = `c'Sales/`c'Acres; 
};

* totals over all crops: only sum where sales prodValue and Acres are nonmissing;
foreach v in Sales Acres ProdValue {;
    local cropListString;
    foreach c of local commodityListAll {;
        replace `c'`v' = . if ( (`c'Sales >= .) | (`c'Acres >= .) | (`c'ProdValue >= .) );
        local cropListString `cropListString' `c'`v';
    };
    egen `v' = rowtotal(`cropListString');
};
gen ProdValueAcres = ProdValue/Acres;
gen SalesAcres = Sales/Acres; 

preserve;
*************************************************************************************************;
* Part 2: regression results for commodity crops;
* Use all counties and cropland harvested;
*************************************************************************************************;
* loop over individual crops;
foreach c of local commodityList {;
    * loop over whether livestock farms are includes or exclude;
    forvalues e = 0/1 {;
	restore; preserve;
	qui keep if (`c'ProdValueAcres > 0) & (`c'ProdValueAcres < .) & (`c'SalesAcres > 0) & (`c'SalesAcres < .) &
	            (`c'Acres > 0) & (`c'Acres < .) & (exclFarmsWithLivestock == `e');
		
        * create state by year fixed effects;
	qui levelsof year, local(yearList);
	qui summ year;
	local yearMin = r(min);
	gen state = floor(fips/1000);
	qui levelsof state, local(stateList);
	foreach s of local stateList {;
	    foreach t of local yearList {;
		* exlcude first year to avoid perfect multicolinearity;
		if (`t' > `yearMin') {;
		    qui gen sst_state`s'_year`t' = (year == `t') & (state == `s');
                };
            };
        };
		
		
        * regress sales on yields from the same crop;
	qui areg `c'SalesAcres  `c'ProdValueAcres  sst* [aweight=`c'Acres], a(fips) robust;
	local regB  = round(_b[`c'ProdValueAcres]*1000)/1000;
	local regSE = round(_se[`c'ProdValueAcres]*1000)/1000;
	local regN  = e(N);
	display("Regression for `c' (farms with livestock excluded = `e') - `regN' observations");
	display("    coefficient on value of production is: `regB' with a s.e. of `regSE'");
	display(" ");
    };
};

*************************************************************************************************;
* Part 2: regression results for commodity crops: sum of all crops;
*************************************************************************************************;
* loop over whether livestock farms are includes or exclude;
forvalues e = 0/1 {;
    restore; preserve;
    qui keep if (ProdValueAcres > 0) & (ProdValueAcres < .) & (SalesAcres > 0) & (SalesAcres < .) &
                (Acres > 0) & (Acres < .) & (exclFarmsWithLivestock == `e');

    * create state by year fixed effects;
    qui levelsof year, local(yearList);
    qui summ year;
    local yearMin = r(min);
    gen state = floor(fips/1000);
    qui levelsof state, local(stateList);
    foreach s of local stateList {;
        foreach t of local yearList {;
            * exlcude first year to avoid perfect multicolinearity;
            if (`t' > `yearMin') {;
	        qui gen sst_state`s'_year`t' = (year == `t') & (state == `s');
            };
        };
    };

    * regress sales on yields from the same crop;
    qui areg SalesAcres  ProdValueAcres  sst* [aweight=Acres], a(fips) robust;
    local regB  = round(_b[`c'ProdValueAcres]*1000)/1000;
    local regSE = round(_se[`c'ProdValueAcres]*1000)/1000;
    local regN  = e(N);
    display("Regression for all crops (farms with livestock excluded = `e') - `regN' observations");
    display("    coefficient on value of production is: `regB' with a s.e. of `regSE'");
    display(" ");
};
restore;

log close;