*! version 0.2.6 04July2017
*! Copyright (C) World Bank 2016-17 - Minh Cong Nguyen & Paul Andres Corral Rodas
*! Minh Cong Nguyen - mnguyen3@worldbank.org
*! Paul Andres Corral Rodas - pcorralrodas@worldbank.org
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
cap program drop povmap
program define povmap, eclass byable(recall)
version 11, missing
mata st_local("StataVersion", lsae_povmapStataVersion()); st_local("CodeVersion", lsae_povmapVersion())
if `StataVersion' != c(stata_version) | "`CodeVersion'" < "00.05.00" {
cap findfile "lsae_povmap.mlib"
while !_rc {
erase "`r(fn)'"
cap findfile "lsae_povmap.mlib"
}
qui findfile "lsae_povmap.mata"
run "`r(fn)'"
}
syntax varlist(min=2 numeric fv) [if] [in] [aw pw fw], area(varname numeric) ///
stage(string) VARest(string) [UNIQid(varname numeric) ETA(string) EPSilon(string) ///
PSU(varname numeric) Zvar(varlist numeric fv) yhat(varlist numeric fv) ///
yhat2(varlist numeric fv) lny plinevar(string) PLINEs(numlist sort) ///
rep(integer 1) seed(integer 123456789) matin(string) matbeta(string) ///
PWcensus(string) ydump(string) ydumpdta(string) numfiles(integer 1) prefix(string) ///
saveold RESults(string) addvars(string) BOOTstrap EBest allmata ///
COLprocess(integer 1) NOOUTput vce(string) INDicators(string) aggids(numlist sort) alfatest(string)]
if c(more)=="on" set more off
local version : di "version " string(_caller()) ", missing:"
local cmdline: copy local 0
set seed `seed'
local colprocess = 1 //temporary code
//if ("`stage'" == "") local stage second
//housekeeping
if "`matin'"!="" { //Check if census mata file exists
cap confirm file "`matin'"
if _rc!=0 {
display as error "`matin' not found, please verify"
error 198
}
}
local epsdecomp `varest'
local process col
local hhid `uniqid'
foreach vs in stage epsilon eta ebest epsdecomp vce glsmethod alfatest {
local `vs' = lower("``vs''")
}
//join addvars and plinevar
if "`plinevar'"!="" & "`plines'"!=""{
dis as error "You must specify only one option: plinevar or plines"
error 198
}
local plinevar: list uniq plinevar
local addvars : list addvars | plinevar
local addvars : list uniq addvars
if ("`plinevar'"!="" & "`allmata'"=="") {
if `: list sizeof plinevar' > 1 {
dis as error "You must specify only one plinevar"
error 198
}
}
if "`indicators'"~="" local indicators `=upper("`indicators'")'
//missing observation check
marksample touse
local flist `"`varlist' `zvar' `yhat' `yhat2' `by' `grvar'"'
markout `touse' `flist'
gettoken lhs varlist: varlist
_fv_check_depvar `lhs'
//Selection of X vars
_rmcoll `varlist' if `touse', forcedrop
local okvarlist `r(varlist)'
//fvexpand `okvarlist'
//local okvarlist `r(varlist)'
local okvarlist : list sort okvarlist
//local okvarlist `varlist' `cons'
//local okvarlist `varlist'
//CHECK: varinmodel put all variables used in the model here to mask
*local varinmodel "`okvarlist' `pwcensus' `area' `addvars'"
local varinmodel "`okvarlist' `pwcensus' `area'"
local varinmodel : list uniq varinmodel
local addvars: list uniq addvars
//Model setting
//hheff household effect alpha model
if "`zvar'"=="" & "`yhat'"=="" & "`yhat2'"=="" local hheffs 0
else local hheffs 1
//Change to ensure that if H3 is chosen only bootstrapped drawing of Beta and alpha is done
if "`bootstrap'"=="" local boots 0
else local boots 1
if "`epsdecomp'"=="h3" & `boots'==0 {
local boots 1
local bootstrap bootstrap
display in yellow "You chose H3, parameters must be obtained via bootstrap I changed it for you."
}
if "`psu'"!="" local psu1 = 1
else local psu1 = 0
//Stage
if ("`ebest'"!="") local ebest 1
else local ebest 0
if "`stage'"=="" | "`stage'"=="second" {
if ("`uniqid'"=="") | ("`eta'"=="") | ("`epsilon'"=="") | ("`pwcensus'"=="") {
display as error "Invalid syntax: uniqid, eta, pwcensus, and epsilon are required when running simulations"
error 198
}
if "`epsilon'"!="normal" & "`epsilon'"!="nonnormal" {
dis as error "You must specify drawing method for epsilon as either epsilon(normal) or epsilon(nonnormal)"
error 198
}
if "`epsilon'"=="normal" local epsnormal = 1
else local epsnormal = 0
if "`eta'"!="normal" & "`eta'"!="nonnormal" {
dis as error "You must specify drawing method for eta as either eta(normal) or eta(nonnormal)"
error 198
}
if "`eta'"=="normal" local etanormal = 1
else local etanormal = 0
if ("`psu'"=="" & "`bootstrap'"!="") {
display in gr "You have not specified a PSU for bootstrapping, cluster level is used by default"
}
//this ebest below is not correct.
if ("`ebest'"=="ebest" & "`eta'"!="normal") {
dis as error "Semi-parametric drawing of eta is not allowed under EB"
error 198
}
if ("`epsdecomp'"=="h3" & "`bootstrap'"!="bootstrap") {
dis as error "Henderson epsilon decomposition is only allowed if bootstrap distribution drawing method is chosen"
error 198
}
} //end stage
if ("`plines'"=="" & "`stage'"=="") {
display in gr "You have not specified a poverty line, all predicted incomes were saved in `ydump'"
}
if ("`epsdecomp'"!="h3" & "`epsdecomp'"!="ell") {
dis as error "Either H3 or ELL method must be specified"
error 198
}
//why this?
if ("`stage'"=="" & "`ydump'"=="" & "`results'"=="") {
dis as error "You have not specified where to save your results. If you'd like to only see the first stage estimates use the stage(first) option."
error 198
}
if ("`epsdecomp'"=="h3") local h3 1
else local h3 0
//local nolocation //what is this
if ("`lny'"=="lny") {
dis in gr _n "Note: Dependent variable in logarithmic form"
local lny 1
}
else {
local lny 0
}
if "`allmata'"=="allmata" local matay 1
else local matay 0
//VCE option: none, robust, cluster
if "`vce'"=="" local vce ell
if "`vce'"~="robust" & "`vce'"~="cluster" & "`vce'"~="none" & "`vce'"~="ell" {
dis in gr "Either vce(robust), vce(cluster), or vce(ell) can be specified; or no vce option can be specified"
error 198
}
if "`vce'"=="none" local vceopt 0
if "`vce'"=="robust" local vceopt 1
if "`vce'"=="cluster" local vceopt 2
if "`vce'"=="ell" local vceopt 3
//GLS method: old ELL or new weighted GLS
if "`glsmethod'"=="" local glsmethod wgls
if "`glsmethod'"~="ell" & "`glsmethod'"~="wgls" {
dis in gr "Either gls(ell) or gls(wgls) option can be specified"
error 198
}
//Weights
local wvar : word 2 of `exp'
if "`wvar'"=="" {
tempvar w
gen `w' = 1
local wvar `w'
}
//Area options
if ("`area'"=="") {
tempvar grv
qui gen `grv' = 1
local grvar `grv'
}
else {
if (`:word count `area''==1) { // one group variable
local grvar `area'
}
else { // more than one group variable
dis as error "Please provide only one hierarchical ID variable."
error 198
}
}
//Remove areas with less than 3 obs per area
tempvar cue
qui egen `cue'=count(_n), by(`grvar')
qui count if `cue'<3
display in yellow "WARNING: `r(N)' observations removed due to less than 3 observations in the cluster."
qui drop if `cue'<3
//Setup the variables for the alpha model
tempname one
qui gen __mz1_000 = .
qui gen __myh_000 = .
qui gen __myh2_000 = .
qui gen `one' = 1
if "`zvar'"=="" local zvarn __mz1_000
else local zvarn : list sort zvar
if "`yhat'"=="" {
local yhatn __myh_000
}
else if "`yhat'"=="_one" {
local yhatn "`one'"
local yhnames yhat
}
else {
local yhatn : list sort yhat
local yhnames
foreach va of local yhatn {
local yhnames "`yhnames' `va'_yhat"
}
}
if "`yhat2'"=="" {
local yhat2n __myh2_000
}
else if "`yhat2'"=="_one" {
local yhat2n "`one'"
local yh2names yhat2
}
else {
local yhat2n : list sort yhat2
local yh2names
foreach va of local yhat2n {
local yh2names "`yh2names' `va'_yhat2"
}
}
//name matrix
local bnames `okvarlist' _cons
local anames `zvarn' `yhnames' `yh2names' _cons
//ID check
if ("`hhid'"=="") {
if lower("`stage'")=="second" {
dis as error "Note: Identifier is not specified. You need to specify the unique identifier."
error 198
exit
}
if lower("`stage'")=="first" {
tempvar _hHiD
gen `_hHiD' = _n
local hhid1 `_hHiD'
}
}
else {
cap isid `grvar' `hhid'
if (_rc!=0) {
dis as error "Warning: Identifier is not unique at cluster level, please provide the unique identifier."
error 198
exit
}
else {
capture confirm numeric variable `hhid'
if (_rc==0) {
local hhid1 `hhid'
}
else {
display as error "Warning: Your household ID in the survey data set is not numeric, please provide numeric IDs."
error 198
exit
}
}
}
qui if ("`stage'"=="first" & "`alfatest'"!="" & "`zvar'"!=""){
gen `alfatest'_alfa =.
lab var `alfatest'_alfa "Dependent variable of alfa model, after `epsdecomp' decomposition"
gen `alfatest'_xb =.
lab var `alfatest'_xb "Yhat prediction from first stage ols"
local alfaerr `alfatest'_alfa
local exbeta `alfatest'_xb
local alfatest=1
}
sort `grvar' `hhid1'
//Estimates from the modeling parts
mata: _s2cs_base0("`lhs'","`okvarlist'", "`zvarn'", "`yhatn'", "`yhat2n'", "`grvar'", "`wvar'", "`hhid1'", "`touse'")
//OLS
mat rownames _vols = `bnames'
mat colnames _vols = `bnames'
mat rownames _bols = `lhs'
mat colnames _bols = `bnames'
mat vols = _vols
mat bols = _bols
di in gr _n "OLS model:"
ereturn post _bols _vols, depname(`lhs')
ereturn local cmd "s2sc"
ereturn display
est store bOLS
if `hheffs'==1 {
//Alpha model
mat rownames _vzols = `anames'
mat colnames _vzols = `anames'
mat rownames _bzols = Residual
mat colnames _bzols = `anames'
mat vzols = _vzols
mat bzols = _bzols
di in gr _n "Alpha model:"
ereturn post _bzols _vzols, depname(Residual)
ereturn local cmd "s2sc"
ereturn display
est store balphaOLS
}
//GLS model
mat rownames _vgls = `bnames'
mat colnames _vgls = `bnames'
mat rownames _bgls = `lhs'
mat colnames _bgls = `bnames'
mat bgls = _bgls
mat vgls = _vgls
di in gr _n "GLS model:"
ereturn post _bgls _vgls, depname(`lhs')
ereturn local cmd "s2sc"
ereturn display
est store bGLS
di in gr _n "Comparison between OLS and GLS models:"
est table bOLS bGLS
//return stuff after estimation
ereturn matrix b_beta = bols
ereturn matrix V_beta = vols
ereturn scalar rmse_beta = _rmse_beta[1,1]
ereturn scalar r2_beta = _r2_beta[1,1]
ereturn scalar F_beta = _F_beta[1,1]
ereturn scalar r2a_beta = _r2a_beta[1,1]
ereturn scalar N_beta = _N_beta[1,1]
ereturn scalar dfm_beta = _dfm_beta[1,1]
ereturn scalar eta_ratio = _v_ratio[1,1]
ereturn scalar eta_var = _v_eta[1,1]
ereturn scalar eps_var = _v_eps[1,1]
if ("`epsdecomp'"=="ell") ereturn scalar var_eta_var = _fvar_sigeta2[1,1]
if `hheffs'==1 {
ereturn matrix b_alpha = bzols
ereturn matrix V_alpha = vzols
ereturn scalar rmse_alpha = _rmse_alpha[1,1]
ereturn scalar r2_alpha = _r2_alpha[1,1]
ereturn scalar F_alpha = _F_alpha[1,1]
ereturn scalar r2a_alpha = _r2a_alpha[1,1]
ereturn scalar N_alpha = _N_alpha[1,1]
ereturn scalar dfm_alpha = _dfm_alpha[1,1]
}
ereturn matrix b_gls = bgls
ereturn matrix V_gls = vgls
display _newline in ye "Model settings"
di as text "{hline 61}"
display in gr "Error decomposition" _col(50) in ye upper("`epsdecomp'")
if (lower("`stage'")=="second") {
if ("`bootstrap'"=="") display in gr "Beta drawing" _col(50) in ye "Parametric"
else display in gr "Beta drawing" _col(50) in ye "Bootstrapped"
if (`ebest'==1) { //Always normal when EB
display in gr "Eta drawing method" _col (50) in ye lower("normal")
display in gr "Epsilon drawing method" _col (50) in ye lower("normal")
display in gr "Empirical best methods" _col (50) in ye "Yes"
}
else {
display in gr "Eta drawing method" _col (50) in ye lower("`eta'")
display in gr "Epsilon drawing method" _col (50) in ye lower("`epsilon'")
display in gr "Empirical best method" _col (50) in ye "No"
}
}
display _newline in ye "Beta model diagnostics"
di as text "{hline 61}"
display in gr "Number of observations" _col(45) in gr "=" _col(50) in ye e(N_beta)
display in gr "Adjusted R-squared" _col(45) in gr "=" _col(50) in ye e(r2a_beta)
display in gr "R-squared" _col(45) in gr "=" _col(50) in ye e(r2_beta)
display in gr "Root MSE" _col(45) in gr "=" _col(50) in ye e(rmse_beta)
display in gr "F-stat" _col(45) in gr "=" _col(50) in ye e(F_beta)
if `hheffs'==1 {
display _newline in ye "Alpha model diagnostics"
di as text "{hline 61}"
display in gr "Number of observations" _col(45) in gr "=" _col(50) in ye e(N_alpha)
display in gr "Adjusted R-squared" _col(45) in gr "=" _col(50) in ye e(r2a_alpha)
display in gr "R-squared" _col(45) in gr "=" _col(50) in ye e(r2_alpha)
display in gr "Root MSE" _col(45) in gr "=" _col(50) in ye e(rmse_alpha)
display in gr "F-stat" _col(45) in gr "=" _col(50) in ye e(F_alpha)
}
display _newline in ye "Model parameters"
di as text "{hline 61}"
display in gr "Sigma ETA sq. = " _col(50) in ye e(eta_var)
display in gr "Ratio of sigma eta sq over MSE = " _col(50) in ye e(eta_ratio)
display in gr "Variance of epsilon = " _col(50) in ye e(eps_var)
if ("`epsdecomp'"=="ell") display in gr "Sampling variance of Sigma eta sq. = " _col(50) in ye e(var_eta_var)
di as text "{hline 61}"
display in ye _col(20) ""
local negeta = e(eta_var)
local negeps = e(eps_var)
if (`negeta'<0){
noi dis as error "Your Sigma ETA sq. value is negative, please revise the model or cluster level"
error 121
}
if (`negeps'<0){
noi dis as error "Your variance of epsilon value is negative, please revise the model"
error 121
}
cap drop __mz1_000 __myh_000 __myh2_000
//Processing the census data
if "`stage'"=="" | "`stage'"=="second" {
if (!("`indicators'"!="" & "`aggids'"!="" & ("`plinevar'"!="" | "`plines'"!="")) & ("`ydump'"=="")) {
noi dis as error "You have requested the second stage to run, but forgot to specify"
noi dis as error "where to save your ydump, or options for processing" _n
noi dis as error "Options include: indicators(); aggids(); and plinevar() or plines()" _n
error 121
}
if "`allmata'"~="" {
//Indicator checklis
if "`indicators'"=="" local indicators fgt0
local indicators = lower("`indicators'")
local fgtlist
local gelist
foreach ind of local indicators {
if "`ind'"=="fgt0" local fgtlist "`fgtlist' `ind'"
if "`ind'"=="fgt1" local fgtlist "`fgtlist' `ind'"
if "`ind'"=="fgt2" local fgtlist "`fgtlist' `ind'"
if "`ind'"=="ge0" local gelist "`gelist' `ind'"
if "`ind'"=="ge1" local gelist "`gelist' `ind'"
if "`ind'"=="ge2" local gelist "`gelist' `ind'"
}
}
if "`ydump'"=="" & "`allmata'"=="" {
tempfile ydumpdata
local ydump `ydumpdata'
}
else { //check we can write or path valid
//some code here to check if we can write the files before going to mata
}
//yhatlist for ydump matadata
local yhatlist
forv v=1(1)`rep' {
local yhatlist "`yhatlist' _YHAT`v'"
}
if "`=lower("`process'")'"=="col" { //column-wise calculation
local _itran
if `=mod(`rep',`colprocess')'~=0 {
noi dis as error "Number of simulation should be divisiable by `colprocess'. Please change the number in the colprocess() option."
error 121
}
tempfile estout logout
di as text "{hline 61}"
di in gr _n "Initializing the Second Stage, this may take a while..." _n
sae_closefiles
mata: _s2sc_sim_cols2("`okvarlist'", "`zvarn'", "`yhatn'", "`yhat2n'", "`grvar'", "`plines'", "`pwcensus'", "`touse'", "`hhid1'", "`matin'")
sae_closefiles
dis in gr _n "Finished running the Second Stage"
di as text "{hline 61}"
//allmata or plugin
if "`allmata'"~="" {
if "`ydump'"=="" dis in gr _n "Results are saved into Stata memory"
else {
if "`plines'"!="" sae_proc_inds, matasource(`ydump') aggids(`aggids') ind(`indicators') plines(`plines')
if "`plinevar'"!="" sae_proc_inds, matasource(`ydump') aggids(`aggids') ind(`indicators') plinevar(`plinevar')
}
}
else { //plugin
//Process ydump
if "`_itran'"=="0" {
local doydump = ("`indicators'"!="" & "`aggids'"!="" & ("`plinevar'"!="" | "`plines'"!=""))
if `doydump'==1 {
local plinesused : list plinevar | plines
dis in yellow _n "Opening plugin"
cap program define povc, plugin using("PovcPlugin`=cond(strpos(`"`=c(machine_type)'"',"64"),64,32)'.dll")
cap plugin call povc, 1 "`ydump'" "`estout'" "`plinesused'" "`aggids'" "`indicators'" "`logout'"
if _rc==0 {
if "`results'"~="" {
cap copy "`estout'" "`results'", replace
if _rc==0 noi dis as text "Output bas been saved in `results'. Please see the results."
else error 691
}
else {
cap insheet using "`estout'", clear double case
if _rc==0 noi dis as text "Output bas been loaded into Stata. Please see the results."
else noi dis as text "Unable to load the result file into Stata."
}
*cap erase "`ydump'"
*cap erase "`estout'"
}
else {
sae_closefiles
cap erase "`ydump'"
cap erase "`estout'"
noi di "`_return_string'"
error 1002
} //rc plugin
} //doydump
else { //dont do ydump
dis in yellow "No indicators were procssed, if you want indicators please provide all necessary options"
dis in yellow "Options include: indicators(); aggids(); and plinevar() or plines()"
}
//copy mata ydump if needed
if "`ydumpdta'"~="" noi sae_data_export, matasource(`ydump') numfiles(`numfiles') prefix(`prefix') datasave(`ydumpdta') `saveold'
}
else {
sae_closefiles
}
} //end plugin
}
//if "`=lower("`process'")'"=="block" { //block-wise calculation
//}
}
end