fixed appending & summary graph

CER-calculate-AR.do

@@ -34,12 +34,16 @@ GNU General Public License for more details.
 
 global where "/Users/ben/Documents/Work"
 global proot "$where/Projects/ESRC-Transformative-Census2022"
-global rpath "$proot/results/CER-Irish-SM-Trial"
+
 * original files
 global odfiles "$where/Data/Social Science Datatsets/CER Smart Metering Project"
+
 * processed files
-global pdfiles "$proot/data/cer"
+global pdfiles "$proot/data/cer/CER_OctHH_data"
+
 * results path
+global rpath "$proot/results/CER-Irish-SM-Trial"
+
 
 global version "v1"
 
@@ -60,93 +64,90 @@ timer on 1
 * see https://github.com/dataknut/Census2022/blob/master/CER-data-processing.do
 
 * load 1/2 hour data (don't need survey for this)
-use "$pdfiles/CER_OctHH_data/CER_Oct2009HH_30min_survey.dta", clear
+use "$pdfiles/CER_Oct2009HH_30min_survey.dta", clear
 
 local midweek "0 1" // 0 = weekends
-local midweekt1 "Mid-week"
+local midweekt1 "Midweek"
 local midweekt0 "Weekends"
 
+* we have 36 values 'per day' (not 48)
+* so lag 35 = this time tomorrow
+* for mid-week, this time next week = lag 105
+* for weekend, this time next week = lag 70
+* so set max lags = 150 to cover both (will also mean 2 weekends)
+local max_lag = 150
+local midweektl0 xline(35 70, lstyle(refline )) 
+* the labels don't seem to show up?
+local midweektll0 text(35 0.5 "This time tomorrow", place(w)) text(70 0.5 "This time next weekend", place(e)) // weekend
+
+local midweektl1 xline(35 70 105, lstyle(refline )) 
+* the labels don't seem to show up?
+local midweektll1 text(35 0 "This time tomorrow", place(w)) text(70 0 "This time the day after tomorrow", place(w)) text(105 0 "This time next week", place(w)) // midweek
 
-local midweektl0 xline(36 72, lstyle(refline )) 
-local midweektll0 text(36 0 "This time tomorrow") text(72 0 "This time next weekend") // weekend
+local econt1 "In work"
+local econt2 "Unemployed"
+local econt3 "Retired"
+local econt4 "Caring for relative or family"
 
-local midweektl1 xline(36 72 108, lstyle(refline )) 
-local midweektll1 text(36 0 "This time tomorrow") text(72 0 "This time the day after tomorrow") text(108 0 "This time next week") // midweek
+* filter for 
+gen s_hour = hh(s_halfhour)
 
 * this creates a big list of all IDs so we can loop over it - takes a while
 * could just start from min & loop to max - but would then be testing for non-existent households
-qui: levelsof ID, local(ids)
+* qui: levelsof ID, local(ids)
 
-local econt1 "In work"
-local econt2 "Unemployed"
-local econt3 "Retired"
-local econt4 "Caring for relative or family"
+* for testing
+local ids "1002 1003 1004"
 
 foreach m of local midweek {
 	di "****************"
-	di "* Calculating ARs for midweek = `m'"
+	di "* Calculating ARs for midweek = `m' (`midweekt`m'') " 
 	foreach id of local ids {
 		di "* -> ID = `id'"
 
 		preserve
-			keep if ID == `id' & midweek == `m' & ds_halfhour >= 12 & ds_halfhour <= 48 // 06:00 - 24:00 -> 36 hour day
-			* qui: su ds_halfhour
-			* check
-			* di "Halfhours: `r(min)' - `r(max)'"
-			* skip over if fails due to lack of hubid
-			capture noisily {
-				* use s_datetime so ac knows what is missing
-				qui: tsset s_datetime
+			* take out 00:00 - 06:00 as repeat cycle of sleep will swamp other cycles
+			keep if ID == `id' & midweek == `m' & s_hour > 6
+			* need to fool ac into thinking these observations are continuous
+			qui: su ID
+			qui: egen lag = fill(1(1)`r(N)')
+			tsset lag
 			* do this here so local is set to the unique valuye for this hh
 			qui: levelsof ba_nchildren, local(nch)
 			qui: levelsof ba_nadults, local(ba_nadults)
 			qui: levelsof ba_empl, local(econ)
-				di "ID: `id' (`midweekt`m'', N = `ba_nadults' adults, `nch' children, respondent: `econt`econ'')"
-				* we have 36 values 'per day' (not 48)
-				* so lag 36 = this time tomorrow
-				* for mid-week, this time next week = lag 108
-				* for weekend, this time next week = lag 72
-				* so set max lags = 144 to cover both (will also mean 2 weekends)
+			qui: su s_hour
+			di "Halfhours: `r(min)' - `r(max)' for ID: `id' (`midweekt`m'', N = `ba_nadults' adults, `nch' children, respondent: `econt`econ'')"
 			* this draws graph - slow
-				qui: ac kwh, gen(archr) lags(150) name(ac_hubid`id'_`midweekt`m'') ///
+			qui: ac kwh, gen(archr) lags(`max_lag') name(ac_hubid`id'_`midweekt`m'') ///
 				`midweektl`m'' /// draw lines
 				`midweektll`m'' /// draw line labels
-					note("ID: `id' (`midweekt`m'', N = `ba_nadults' adults, `nch' children, respondent: `econt`econ'')")
+				xlabel(0(35)140) ///
+				note("Halfhours: `r(min)' - `r(max)' for ID: `id' (`midweekt`m'', `ba_nadults' adult(s), `nch' children, respondent: `econt`econ'')")
 			graph export "$rpath/graphs/kwh_archr_hubid`id'_`midweekt`m''_$version.png", replace
 			
 			di "* save out arch results for household ID = `id'"
-				qui: su t
-				qui: egen lag = fill(1(1)`r(max)')
 			qui: keep lag archr ID
 			qui: save "$rpath/tmp/archr-`id'_`midweekt`m''_$version.dta", replace
 			di "* Done"
-			}
 		restore
 	}
 }
 
-* reset the id loop values
-qui: su ID
-local id_min = `r(min)'
-local id_max = `r(max)'
-local id = `r(min)'
-
 clear
 
 * pool all the results
 foreach m of local midweek {
 	di "****************"
 	di "* Loading `midweekt`m'' files"
-	while `id' < `id_max' {
-		capture noisily {
+	foreach id of local ids {
 		append using "$rpath/tmp/archr-`id'_`midweekt`m''_$version.dta"
 		erase "$rpath/tmp/archr-`id'_`midweekt`m''_$version.dta"
 	}
 	* distribution of average values for first 70?
-		graph box archr if lag < 144, over(lag)
-		graph export "`rpath'/graphs/box_mean_archr_all_hubids_`midweekt`m''_`version'.png", replace
-		save "`rpath'/archr-all_hubids_`midweekt`m''_`version'.dta", replace
-	}
+	graph box archr if lag <= 105, over(lag)
+	graph export "$rpath/graphs/box_mean_archr_all_hubids_`midweekt`m''_$version.png", replace
+	save "$rpath/archr-all_hubids_`midweekt`m''_$version.dta", replace
 }