/* Multi-site count data model */

#delimit ; 

set mem 30m;
set matsize 80;
set logtype text;
set more off; /* means: produce all output at once */

log using c:\klaus\apec64\spring_07\stata\logs\tc_multiple,replace;
use c:\klaus\apec64\spring_07\stata\data\multi_sites3,replace;

/* STEP 1: variable definitions and sample statistics*/

describe;
sum;
table site, c(sum visits) row; /* get total visits by site plus a row total*/
sort site;
by site: tab visits;
by site: sum visits;

/* STEP 2: run multi-site Poisson model and examine marginal effects*/

poisson visits s3 s4 s5 gender age fishyrs flyfish tc3 tc4 tc5 inc0000, nocon;

/* note: "nocon" means "run the model without a general constant term (intercept); since we already have site-specific intercepts, this would cause a "dummy-trap" error */

/* Step 3a: for each site, predict trips at current situation (status quo), and examine the distribution of predictions */
/* note: because we need to use different coefficients for S and tc for each site, we can't just use "predict" anymore... */

gen lhat3 = exp(_b[s3]*s3+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc3]*tc3);
/* note: this created lhat3 for all observations, but only those for site=3 are meaningful & will be used below */
gen lhat4 = exp(_b[s4]*s4+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc4]*tc4);/* same for site=4*/
gen lhat5 = exp(_b[s5]*s5+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc5]*tc5);/*same for site=5*/

sum lhat3 if site==3; /* so now we're only considering site3-stuff */
sum lhat4 if site==4; /* same for site 4 */
sum lhat5 if site==5; /* same for site 5 */

/*STEP 3b: summarize trips for entire system*/
gen lhat =.; /*first, generate just missing values */
replace lhat=lhat3 if site==3;/* then replace them with site-specific lhat values*/
replace lhat=lhat4 if site==4;
replace lhat=lhat5 if site==5;

sum lhat;


/*STEP 4a: derive & summarize seasonal welfare at status quo for each site*/
/* note: focus on CS since income effect is negative, which complicates use of CV, EV */
gen cs3 =-(1/_b[tc3])*lhat3; /* this creates cs3 for all observations, but only those corresponding for site 3 are meaningful & will be used below */
gen cs4 =-(1/_b[tc4])*lhat4; /* same for site 4*/
gen cs5 =-(1/_b[tc5])*lhat5; /* same for site 5*/

sum cs3 if site==3; /* so now we're only considering site3-stuff */
sum cs4 if site==4; /* same for site 4 */
sum cs5 if site==5; /* same for site 5 */

/*STEP 4b: derive & summarize seasonal welfare for entire system*/
gen cs =.; /*first, generate just missing values */
replace cs=cs3 if site==3;/* then replace them with site-specific cs values*/
replace cs=cs4 if site==4;
replace cs=cs5 if site==5;

sum cs;


/* ************************************************************************** */
/* Assume the policy scenario is "implement access fees of $5 at all 3 sites" */
/* ************************************************************************* */

/*STEP 5a: Create new travel cost variables and predict visits under these fees*/
gen tc3_new = tc3+5; /* add $5 to everybody's travel cost for each site*/
gen tc4_new = tc4+5;
gen tc5_new = tc5+5;

gen lhat3_new = exp(_b[s3]*s3+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc3]*tc3_new);
gen lhat4_new = exp(_b[s4]*s4+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc4]*tc4_new);
gen lhat5_new = exp(_b[s5]*s5+ _b[gender]*gender +_b[age]*age + _b[fishyrs]*fishyrs + _b[flyfish]*flyfish + _b[inc0000]*inc0000 + _b[tc5]*tc5_new);

sum lhat3_new if site==3; /* so now we're only considering site3-stuff */
sum lhat4_new if site==4; /* same for site 4 */
sum lhat5_new if site==5; /* same for site 5 */

/*STEP 5b: summarize trips for entire system under new fees*/
gen lhat_new =.; /*first, generate just missing values */
replace lhat_new=lhat3_new if site==3;/* then replace them with site-specific lhat values*/
replace lhat_new=lhat4_new if site==4;
replace lhat_new=lhat5_new if site==5;

sum lhat_new;



/* STEP 6a: Compute the per-person welfare loss under the new fee for each site*/
gen cs3_new =-(1/_b[tc3])*(lhat3-lhat3_new);
gen cs4_new =-(1/_b[tc4])*(lhat4-lhat4_new); /* same for site 4*/
gen cs5_new =-(1/_b[tc5])*(lhat5-lhat5_new); /* same for site 5*/

sum cs3_new if site==3; 
sum cs4_new if site==4; 
sum cs5_new if site==5; 

/*STEP 6b: derive & summarize seasonal welfare for entire system under the new fees*/
gen cs_new =.; /*first, generate just missing values */
replace cs_new = cs3_new if site==3;/* then replace them with site-specific cs values*/
replace cs_new=cs4_new if site==4;
replace cs_new=cs5_new if site==5;

sum cs_new;

log close;