# delimit ;
clear;
set memory 1g;
set more off;
version 9;
log using hedonicRobustness.log, replace;

* ###############################################################################################;
* script looks at robustness of hedonic model (Section 3);
* ###############################################################################################;

local yearList 1969 1974 1978 1982 1987 1992 1997 2002;

*************************************************************************************************;
* Part 1: build data set;
*************************************************************************************************;
* load data;
use dataHedonic, clear;

* create variables of interest;
foreach t of local yearList {;
    qui gen logFarmVal`t' = log(farmvalue`t');
    qui gen popDens`t'_sq = popDens`t'^2;
};
gen state = floor(fips/1000);

* climate change predictions;
gen d_dday8_32    = dday8_32B2-dday8_32;
gen d_dday8_32_sq = dday8_32B2_sq-dday8_32_sq;
gen d_dday34_root = dday34B2_root-dday34_root;
gen d_prec        = precB2-prec;
gen d_prec_sq     = precB2_sq-prec_sq;

*************************************************************************************************;
* Part 2: regression results;
*************************************************************************************************;
* loop over years;
foreach t of local yearList {;
    * run regression;
    qui regress logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq, robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impA0_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

foreach t of local yearList {;
    * run regression;
    qui regress logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq [aweight = farmland`t'], robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impA1_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

foreach t of local yearList {;
    * run regression;
    qui regress logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq
                       latitude incomeCapita`t' popDens`t' popDens`t'_sq
                       waterCapacity pctClay minPerm topKfact bestSoil, robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impB0_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

foreach t of local yearList {;
    * run regression;
    qui regress logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq
                       latitude incomeCapita`t' popDens`t' popDens`t'_sq
                       waterCapacity pctClay minPerm topKfact bestSoil [aweight = farmland`t'], robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impB1_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

foreach t of local yearList {;
    * run regression;
    qui areg logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq
                       latitude incomeCapita`t' popDens`t' popDens`t'_sq
                       waterCapacity pctClay minPerm topKfact bestSoil, a(state) robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impC0_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

foreach t of local yearList {;
    * run regression;
    qui areg logFarmVal`t' dday8_32 dday8_32_sq dday34_root prec prec_sq
                       latitude incomeCapita`t' popDens`t' popDens`t'_sq
                       waterCapacity pctClay minPerm topKfact bestSoil [aweight = farmland`t'], a(state) robust;

    * derive farmland-area weighted sum of change in climate variables;
    egen sumLand = sum(farmland`t');
    qui egen hilf = sum(d_dday8_32*farmland`t'/sumLand);
    local diff_dday8_32 = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday8_32_sq*farmland`t'/sumLand);
    local diff_dday8_32_sq = hilf[1]; drop hilf;
    qui egen hilf = sum(d_dday34_root*farmland`t'/sumLand);
    local diff_dday34_root = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec*farmland`t'/sumLand);
    local diff_prec = hilf[1]; drop hilf;
    qui egen hilf = sum(d_prec_sq*farmland`t'/sumLand);
    local diff_prec_sq = hilf[1]; drop hilf;
    drop sumLand;

    * derive impact on county;
    gen impC1_`t' = ((exp(_b[dday8_32]*d_dday8_32 + _b[dday8_32_sq]*d_dday8_32_sq + _b[dday34_root]*d_dday34_root + _b[prec]*d_prec + _b[prec_sq]*d_prec_sq)-1)*100);

    * derive aggregate climate impact;
    disp("################################################################################################################");
    disp("no controls (Panel A0): unweighted results for year `t'");
    nlcom ((exp(_b[dday8_32] * `diff_dday8_32' + _b[dday8_32_sq] * `diff_dday8_32_sq' + _b[dday34_root] * `diff_dday34_root' + _b[prec] * `diff_prec' + _b[prec_sq] * `diff_prec_sq')-1)*100);
};

log close;