-
Notifications
You must be signed in to change notification settings - Fork 3
/
pmt2.ado
executable file
·178 lines (147 loc) · 8.16 KB
/
pmt2.ado
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
program define pmt2, rclass
// July 26 2010
// July 31 2010 - this routine uses the pmt_eligible_train_test.ado program. I've split this into two, incase one wishes to evaluate the performance of non-regression models
// August 03 2010 : I want to replace the sum command with one that uses survey settings. This works well. Next one will add option to use betas from
// a subsample to predict consumption for the rest of the sample
// August 04 2010 : Adding option to use betas from a subsample. STILL NEED TO VERIFY THAT IT IS WORKING PROPERLY. Also removing reduncancy for pmt_eligible in having to specifiy the Quantiles variable AND the number of Quantiles.
// August 06 2010 : For threshold, will use bottom X percent of ACTUAL (not predicted) consumption of the ENTIRE sample
// August 27: add option to specify whether the the subsamples are a filter - i.e. we toss out observations that are not in the subsample, when calculating leakage, undercoverage, coverage of the difference quantiles, etc.
// March 23 2012 : The cutoff is specified as a percentile of the predicted distribution. I will add the option to specify the cutoff as an absolute number. User should specify either C or CAB
syntax varlist(min=2 ts) [if] [pw aw iw fw], Poor(varname numeric) quantilec(varname numeric) [Cutoffs(numlist asc min=1 max=20 >0 <=70 integer) CABsolute(numlist asc min=0 max=20 >=0) Graphme(real -1) logpline(real 0) Usesubsamplebetas(varname numeric) Filter GENerate(name)]
version 9.1
marksample touse
// Need to get the weights, so that I can use this with the _pctile command below.
qui svyset
local svyweight = r(wtype)
local svyexp = r(wexp)
di "[`svyweight'`svyexp']"
// Was the option to use betas from a subsample specified?
if ( "`usesubsamplebetas'"=="") {
// di "usesubsamplebetas NOT specified"
tempvar usesubsamplebetas
qui gen `usesubsamplebetas' = 1
qui sum `usesubsamplebetas'
}
tempvar logpccd_predicted logpccd_predicted_threshold
qui spearman `varlist' `if'
return scalar rho = r(rho)
qui svy : reg `varlist' if `touse' & `usesubsamplebetas'==1
local R2 = e(r2)
local N = e(N)
return scalar r2 = `R2'
return scalar N = `N'
// predict for the entire sample... may as well
predict `logpccd_predicted'
if ( "`generate'"!="") {
cap drop `generate'
predict `generate' if `touse'
}
qui count if `logpccd_predicted' ~= .
qui local count_predicted = r(N)
// should we loop over percentile or absoulute cutoff?
// check which option user chose
di "cutoffs = `cutoffs'"
if ("`cutoffs'"=="") {
di "log cutoffs not specified - assume absolute cutoff specified"
local cutoffs_to_loopover `cabsolute'
}
di "cabsolute = `cabsolute'"
if ("`cabsolute'"=="") {
di "cabsolute not specified - assume log cutoff specified"
local cutoffs_to_loopover `cutoffs'
}
if ("`cutoffs'"=="" & "`cabsolute'"=="") {
stop
}
// foreach x of numlist `cutoffs' { // loop over all the cutoffs to be used
foreach x of numlist `cutoffs_to_loopover' {
if ("`cutoffs'"~="") {
// Getting percentiles of the ENTIRE population. This part is only for setting the cutoff.
_pctile `1' [`svyweight'`svyexp'], n(100)
// return list
local logcutoff = r(r`x')
di "cutoff: `logcutoff'"
}
else {
local logcutoff = ln(`x')
di "cutoff: `logcutoff'"
}
// for return values, need to fix `x'
// local x2 = substr("`x'",1,5)
local x2 = round(`x')
tempvar eligible nonpoor_ineligible poor_ineligible nonpoor_eligible poor_eligible // Quantilec
qui gen `eligible' =.
qui gen `nonpoor_eligible' = .
qui gen `poor_ineligible' = .
qui replace `eligible' = 1 if `logpccd_predicted' < `logcutoff' & `logpccd_predicted' ~= . & `touse'
qui replace `eligible' = 0 if `logpccd_predicted' >= `logcutoff' & `logpccd_predicted' ~= . & `touse'
// Make ineligible those who are filtered out, if the filter option has been selected
if ("`filter'"=="filter") {
replace `eligible' = 0 if `usesubsamplebetas'==0
}
qui replace `nonpoor_eligible' = 1 if `eligible' == 1 & `poor' ~= 1 & `poor' ~= . & `touse'
qui replace `nonpoor_eligible' = 0 if (`eligible' ~= 1 | `poor' == 1) & `touse'
qui replace `poor_ineligible' = 1 if `eligible' == 0 & `poor' == 1 & `poor' ~= . & `touse'
qui replace `poor_ineligible' = 0 if (`eligible' == 1 | `poor' == 0) & `poor' ~= . & `touse'
qui gen `nonpoor_ineligible' =.
replace `nonpoor_ineligible' = 1 if `poor'==0 & `eligible'==0 & `touse'
replace `nonpoor_ineligible' = 0 if (`poor'==1 | `eligible'==1) & `touse'
qui gen `poor_eligible' =.
replace `poor_eligible' = 1 if `poor'==1 & `eligible'==1 & `touse'
replace `poor_eligible' = 0 if (`poor'==0 | `eligible'==0) & `touse'
// what percent of the population is covered?
qui svy : mean `eligible' if `touse'
local percentofpop_covered = _coef[`eligible']
return scalar percentofpop_covered_`x2' = `percentofpop_covered'
// First generate quantile indicators
// Here the classification of quantile is done based on the entire sample. This is probably what most people want.
/*
qui xtile `Quantilec' = `1' [`svyweight'`svyexp'], n(`quantiles')
*/
// need to get the number of quantiles
svy: tab `quantilec'
local quantiles = e(r)
//
// [`weight'`exp'] is not used, but the darn thing doesn't work unless I include it!
pmt_eligible `eligible' [`svyweight'`svyexp'], p(`poor') qu(`quantilec')
// No need to make adjustments for filters here
// Need to store the results somewhere, somehow
return scalar leakage_`x2' = r(leakage)
return scalar undercoverage_`x2' = r(undercoverage)
return scalar targeting_accuracy_`x2' = r(targeting_accuracy)
return scalar targeting_accuracy2_`x2' = r(targeting_accuracy2)
return scalar coverage_rate_targetgroup_`x2' = r(coverage_rate_targetgroup)
return scalar fractiontotal_targetgroup_`x2' = r(fractionoftotal_in_targetgroup)
return scalar fraction_covered_`x2' = r(fraction_covered)
return scalar inclusion_error_rate_`x2' = r(inclusion_error_rate)
return scalar exclusion_error_rate_`x2' = r(exclusion_error_rate)
forv q = 1/`quantiles' {
return scalar coverage_cutoff`x2'_quantile`q' = r(coverage_cutoff_quantile`q')
local temp = r(coverage_cutoff_quantile`q')
// di "scalar coverage_cutoff`x'_quantile`q' = `temp'"
}
tempname xlimlower xlimupper ylimlower ylimupper
local `xlimlower' = 5
local `xlimupper' = 7
local `ylimlower' = 5
local `ylimupper' = 7
if (`graphme'~= -1) {
if ("`filter'"=="filter") {
twoway (scatter `logpccd_predicted' `1' if `poor_ineligible'==1, xline(`logpline', lcolor(0)) yline(`logcutoff') mc(red) m(x) ) /*
*/ (scatter `logpccd_predicted' `1' if `nonpoor_ineligible' == 1, mc(green) m(x)) /*
*/ (scatter `logpccd_predicted' `1' if `nonpoor_eligible'==1, mc(black) m(x)) /* (scatter `logpccd_predicted' `1' if `Quantilec'==10 & logpccd_predicted < `logcutoff' & `1' ~= . & `logpccd_predicted' ~= ., mc(blue) m(Oh))
*/ (scatter `logpccd_predicted' `1' if `poor_eligible' == 1, mc(blue) m(x) xlabel(8(1)12) ylabel(8(1)12) ysc(r(8 12)) xsc(r(8 12)) ) /*
*/ (scatter `logpccd_predicted' `1' if `usesubsamplebetas' == 0, mc(purple) m(oh)) /*
*/ , title("Cutoff `x' pctile (TJK 2009)") legend(lab(1 "Errors of Exclusion") lab(3 "Errors of Inclusion") lab(2 "Nonpoor, Ineligible") lab(4 "Poor Eligible") lab(5 "Filtered Out") )
}
else {
twoway (scatter `logpccd_predicted' `1' if `poor_ineligible'==1, xline(`logpline', lcolor(0)) yline(`logcutoff') mc(red) m(x) ) /*
*/ (scatter `logpccd_predicted' `1' if `nonpoor_ineligible' == 1, mc(green) m(x)) /*
*/ (scatter `logpccd_predicted' `1' if `nonpoor_eligible'==1, mc(black) m(x)) /* (scatter `logpccd_predicted' `1' if `Quantilec'==10 & logpccd_predicted < `logcutoff' & `1' ~= . & `logpccd_predicted' ~= ., mc(blue) m(Oh))
*/ /*
*/ (scatter `logpccd_predicted' `1' if `poor_eligible' == 1, mc(blue) m(x) xlabel(``xlimlower''(1)``xlimupper'') ylabel(``ylimlower''(1)``ylimupper'') ysc(r(``ylimlower'' ``ylimupper'')) xsc(r(``xlimlower'' ``xlimupper'')) ) /*
*/ , title("Cutoff `x'") legend(lab(1 "Errors of Exclusion") lab(3 "Errors of Inclusion") lab(2 "Nonpoor, Ineligible") lab(4 "Poor Eligible") )
}
}
}
end program