An Ordered Probit PMT Model

A proxy means test (PMT) predicts a household's welfare from things a programme officer can see or ask at the door: the house, its water and sanitation, the assets, who lives there. The welfare measure itself, per-capita consumption expenditure, is only known from a survey. The PMT is fitted on survey data and then applied to households outside the survey.

In this assignment the survey is Susenas March, four rounds. You will:

• read the four rounds as BPS delivers them (dBase files, several files per round) and build one household-level analysis file per round;
• fit two linear PMT models on the pooled rounds 2022–2024 with district fixed effects: the full BPS/World Bank specification of thirty-nine proxies, and a parsimonious specification of eight underlying variables;
• fit the same two specifications as ordered probit models on the welfare decile;
• use each of the four models to predict the 2025 round, select the predicted bottom half of the population, and compute the inclusion error of that rule (equal, here, to its exclusion error).

The data are synthetic: no record belongs to a real household. Treat them as you would treat real microdata; the skills are the same.

The four rounds are split in two. Every model is estimated on 2022, 2023 and 2024 pooled. This assignment calls those three rounds the training sample. No model is estimated on 2025; that round is used only to measure how well the predictions work, and is called the test sample. The words come from the machine learning literature. In econometric terms the training sample is the estimation sample, and everything measured on the test sample is out-of-sample prediction. The do files use a flag called train, so the words are kept here.

The split exists because a PMT predicts consumption for households whose consumption has not been measured, so its performance has to be measured on households the model did not use. An R² computed on the estimation sample cannot do that: it never falls when a proxy is added, and part of it comes from noise in that sample which will not repeat.

The round held out is the latest one rather than a random subsample of the four, because a programme estimates the formula once and then uses it for several years. Over those years real welfare changes, the questionnaire changes, and the relation between the proxies and consumption drifts. A random split puts the same years on both sides and hides all of it; holding out 2025 measures it. For the same reason the 2025 decile is computed within 2025 alone (section 4.6): if 2025 were ranked together with the training rounds, information from the test period would enter the training ranking.

From step 2 on, the two samples sit in one file with the flag train: train == 1 on the 2022–2024 rows, train == 0 on the 2025 rows. Estimate every model with if train == 1, predict for every row including the test rows, and measure every error on train == 0 only.

One folder per round, laid out as a Silastik download from BPS:

2022_maret/  2023_maret/  2024_maret/  2025_maret/
    KOR/            the household file (RT) and the individual file(s) (IND), dBase
    KP/             consumption: BLOK43 (household totals) and BLOK42 (three items only)
    kuesioner/      the questionnaires of the round
    metadata_bps/   the BPS variable lists
    Metadata_Synthetic_Susenas_YYYYMM.xlsx
                    the variable list of this release, with value labels
poverty_line_2022_2025.csv
                    the published BPS poverty line by year, province and urban/rural

2.1 Read the questionnaire and the metadata before the data

Every variable in these files is an answer to a question on a form. The habit this exercise builds is to go to the form first, and to the data second. In each round's folder:

• kuesioner/ holds the two questionnaires of the round: VSENyy.K, the core (KOR) questionnaire with the household and individual blocks, and VSENyy.KP, the consumption questionnaire. Find the question you need — the block (Roman numeral) and the item number within it. A field is named after its question: r1804 is block XVIII, item 4 ("floor area of the dwelling"); r407 is block IV, item 7 (age). When a block moves between rounds, every field in it is renamed.
• The metadata workbook Metadata_Synthetic_Susenas_YYYYMM.xlsx lists what the files of that round actually contain: the sheets variabel ruta (household file), variabel individu, Blok 43 and Blok 42 give each field's name, label, type and width; the sheet value label gives the answer codes (which value means "owned", "PLN metered", "yes"). metadata_bps/ holds the corresponding original BPS variable lists.

Working order for every field you use: the question in the questionnaire → its field name and codes in the workbook of that round → a tabulation of the field in Stata to see that the codes are what the workbook says → only then the code. Variable definitions in section 5 name the 2022 fields; for the other rounds the workbook is the authority, not this page.

The four files of a round carry different units of observation:

Where to get the data

The package is one zip, 323 MB, on the course repository's release page:

• synthetic_susenas_2022-2025_stata_release.zip — the four rounds as dBase files, questionnaires, metadata, poverty linezip · 323 MB

Unzip it anywhere. The folder that holds the four YYYY_maret folders and poverty_line_2022_2025.csv is the path to put in the pkg line at the top of s0_import.do.

Optional: the Stata program that generates the synthetic data is here. Not needed for the assignment.

• synthetic_susenas_mc_generator_student.zipzip · 29 MB

Every file is complete for the first case and stops where the pattern repeats. The places where your code goes are marked >>> write ... here <<<. There is no master file: run the four files one after another, in order, each with do. Each file opens with the one or two path lines to edit and writes its own log into the work folder.

• s0_import.doEdit two path lines: pkg, where the data package was unzipped, and work, a folder for the Stata files. Reads every dBase file with import dbase, appends the files BPS splits by province range, saves one .dta per file, then reads the poverty line.given March 2022, the poverty line  ·  yours March 2023, 2024, 2025
• s1_build.doReads the deflator base from the poverty-line table, then one household record per round: collapse the person file into counts, reshape BLOK42 from long to wide, merge the four files, check the result, build the model variables. Then rank the deciles.given March 2022, the deciles  ·  yours March 2023, 2024, 2025
• s2_pmt_ols.doModel 1, the PMT with district fixed effects; model 2, the parsimonious model with district fixed effects. Both predict every household including 2025.given the sample, the recodes  ·  yours models 1 and 2
• s3_pmt_ordered.doModel 3, ordered probit on the decile with the PMT proxies; model 4, the same with the parsimonious set; section 5, the inclusion error (= exclusion error) at the bottom 50 per cent for all four models.given the load and the section comments  ·  yours models 3 and 4, all of section 5

In steps 0 and 1 the given 2022 block is the pattern: follow it line by line for the other rounds and change only what the data require. In steps 2 and 3 no model is given; section 6 specifies each one and section 7 the targeting rule, precisely enough to write the code from. Keep every variable name as specified: later steps read them.

Each round becomes one file, analysis_YYYY.dta, with one row per household. The steps, in the order the given 2022 block does them:

4.1 Person counts

From the individual file, build 0/1 indicators per person and collapse (sum) them by household. The head of the household is the person with relationship code 1; the head's age and certificate are taken with (max) after keeping them only on the head's row. Two questions asked of every person, phone ownership and internet use, are carried as household characteristics by taking the first listed member. Count the members with (count); this count is checked later against the household file.

4.2 BLOK42 from long to wide

BLOK42 has one row per household per item. Keep the three items, put the quantity of electricity and the values of phone credit and internet packages into three variables, and collapse (max) by household. Households that bought none of the three items still have one all-zero row in the module; keep that row, so that these households come out with zeros rather than disappear.

4.3 Merge

Start from the household file, rename the housing and asset fields to the common names of section 5, then merge one-to-one on urut: BLOK43 (keep food nonfood expend kapita wert weind), the person counts, the BLOK42 variables. Then merge the poverty line many-to-one on year, province and urban/rural. Every merge must match every household; the given code keeps matched rows only and the checks below catch a failed merge.

4.4 Check before trusting

Household size can be got three ways: round(expend / kapita), the household file's r301, and the number of person rows. All three agree in these data. The given code asserts this, and that food + nonfood = expend and weind = wert × hhsize. If an assert fails in a round you wrote, the error is in that round's merge or collapse.

4.5 The welfare measure

Monthly per-capita expenditure is deflated by the poverty line of the household's year, province and urban/rural cell, and re-expressed in the prices of one cell: DKI Jakarta urban, March 2022, whose poverty line is Rp 738 955 (the do file reads that value from the poverty-line table; it is not typed in). One decile ranking then means the same living standard in every province and every year.

4.6 Deciles

Two rules. The training rounds 2022, 2023 and 2024 are appended and ranked together into ten weighted deciles of rpcexp, so that decile 1 means the same in each year. The 2025 round is ranked on its own: it is the test round, and nothing from it may enter the training ranking. Deciles are weighted with weind, the individual weight, so that a decile is a tenth of people, not of households. weind is rounded to an integer first because xtile takes frequency weights.

Field names are those of the March 2022 files. The same variable may have another name in another round; the definition does not change. Categorical variables enter the models as factor variables (i.) with their original codes.

5.1 Household file (RT)

5.2 Counts from the individual file (IND)

5.3 Consumption (BLOK43, BLOK42)

5.4 Transformed terms

The +0.1 inside each logarithm keeps households with a zero in the sample; without it ln(0) is missing and the household drops out of every model.

5.5 Structural missing values

Four variables are blank for a reason, not at random: h_closet and h_sewage when the household has no toilet, h_epower when lighting is not PLN metered, and hd_edu for heads with no certificate recorded. A factor variable with a blank drops the household from the regression without a message. Step 2 recodes all four to the explicit category 99 and prints how many households each recode keeps; those counts are part of the reference output.

This section is the specification you write the code from. Nothing in steps 2 and 3 beyond the sample assembly and the recodes is given.

6.1 Sample and variables common to all four models

• Estimation sample: the pooled training rounds 2022, 2023 and 2024 — the households with train == 1 in the file that step 2 assembles: 1 024 865 households. The 2025 round (train == 0, 343 460 households) is never used in estimation.
• Prediction: every model predicts for every household in the file, the 2025 round included. In Stata, predict after an estimation command does this by default; do not restrict it to the estimation sample.
• Urban dummy: urban enters every model as a regressor.
• District fixed effects: every model carries a separate intercept for each of the 514 districts, entered as the factor variable i.kode_kab (one dummy per district, the first as the base). They are part of the model and of the printed output.
• Factor variables: the categorical housing variables (h_house h_roof h_wall h_floor h_toilet h_closet h_sewage h_water h_light h_fuel h_epower) enter as factor variables, i., with their original codes and the recoded 99 as a category of its own. All other regressors enter as they are.
• Output: let Stata print its full output for every model, as it does when the estimation command is not silenced. The reference output shows it in full, district coefficients included.

6.2 The two regressor sets

The PMT — thirty-nine proxies (the global $baseline_39, given in the do file). Four continuous terms, eleven factor variables, thirteen asset dummies, six household counts, four schooling and marital counts, and three head characteristics:

ln_hhsize ln_hhsize_2 ln_h_luaslantai ln_h_luaslantai_2
i.h_house i.h_roof i.h_wall i.h_floor i.h_toilet i.h_closet i.h_sewage
i.h_water i.h_light i.h_fuel i.h_epower
h_asset_a h_asset_b h_asset_c h_asset_d h_asset_e h_asset_f h_asset_g
h_asset_h h_asset_i h_asset_j h_asset_k h_asset_l h_asset_m
n_age04 n_age519 n_age2064 n_age65up n_female n_worker
h_married h_widowed h_neverschool h_stillschool
head_female hd_age hd_edu

The parsimonious model — fifteen terms from eight underlying variables (the global $predictors_7): household size, electricity per head and telecom spending per head each with its square; the two meter-capacity dummies and their interactions with electricity; three assets; floor area with its square:

ln_hhsize ln_elec_kwh_pc ln_telecom_cons_pc
ln_hhsize_2 ln_elec_kwh_pc_2 ln_telecom_cons_pc_2
h_epower2 h_epower3 ln_elec_kwh_pc_h_epower2 ln_elec_kwh_pc_h_epower3
h_asset_car h_asset_computer h_asset_land
ln_h_luaslantai ln_h_luaslantai_2

Every variable is defined in section 5 and built in step 1; the globals are the same lists, written out.

Model 1 — the PMT, linear, district fixed effects

• Estimator: ordinary least squares (Stata's regress).
• Dependent variable: lnrpcexp, the log of real per-capita expenditure of section 4.5.
• Regressors: X = the thirty-nine proxies, plus urban, plus the district dummies i.kode_kab.
• Sample: train == 1.
• Prediction: the linear prediction X_hβ̂ + α̂_d(h) + γ̂·urban_h for every household, stored as a double-precision variable named yhat_pmt. Step 3 reads that name.

Model 2 — the parsimonious model, linear, district fixed effects

The same equation, estimator, sample and district dummies, with X = the fifteen parsimonious terms. Prediction stored as yhat_pars.

Step 2 ends by saving the file with both predictions; the given save line does that.

Model 3 — ordered probit on the decile, the PMT, district fixed effects

• Estimator: ordered probit by maximum likelihood (Stata's oprobit). One slope vector and nine cutpoints κ₁ … κ₉ are estimated; Stata prints the cutpoints at the foot of the coefficient table.
• Dependent variable: decile, the training decile of section 4.6 (1 = poorest tenth of people, 10 = richest).
• Regressors and sample: exactly as model 1 — the thirty-nine proxies, urban, i.kode_kab, train == 1.
• Predicted probabilities: for every household, the ten predicted probabilities Pr(decile = k), k = 1 … 10, as ten double-precision variables named op1 … op10 (Stata's predict with the probability option writes all ten at once when it is given ten new variable names).
• The expected decile: one number per household, the probability-weighted mean decile, stored as op_expected:
  op_expected_h = Σ_k=1..10 k · Pr(decile_h = k | X_h) = 1·op1 + 2·op2 + … + 10·op10
  Summarize it on the 2025 test round (train == 0); that summary (Obs, Mean, Std. dev., Min, Max) is in the reference output. This is the model's score; its predicted decile is made in section 7 by ranking this score with xtile, the same way as the linear prediction — not by rounding it and not by taking the most likely decile.
• Runtime: the model carries the proxies and 514 district dummies; it converges in four or five iterations, about three minutes on a four-core machine. Do not interrupt it.

Model 4 — ordered probit on the decile, the parsimonious model, district fixed effects

The same as model 3 with X = the fifteen parsimonious terms. Probabilities named ops1 … ops10, expected decile ops_expected, summarized on the test round.

Section 5 of step 3 applies one targeting rule to each of the four predictions and reports its error. Write it for all four scores: yhat_pmt, yhat_pars, op_expected, ops_expected. The procedure below is complete; the numbers in it are checkpoints from the reference run, so you can verify each step before the next.

Step 0 — where you are

You are in step 3, in the file that step 2 saved (step2_predictions.dta): the four rounds in one file with the flag train, the two linear predictions from step 2, and, after models 3 and 4, the twenty probabilities and the two expected deciles. Everything in this section is computed on the test round only, train == 0: 343 460 households, who stand for 287 771 246 people (the sum of weind over those households).

Step 1 — the true decile (already built)

The true decile is the variable decile that step 1 built for the 2025 round: real per-capita expenditure rpcexp ranked into ten weighted deciles within 2025, with Stata's xtile, ten quantiles, and the individual weight weind as a frequency weight (weind was rounded to an integer for this). Decile 1 is the poorest tenth of people, decile 10 the richest; each decile holds about 28.8 million people but a different number of households, because poor households are larger.

The eligible households are those with decile from 1 to 5: the true bottom half of the population. They stand for E = 143 887 916 people — half of 287.8 million, as they must.

Step 2 — the score of each model

A score is one number per household; a higher score means a better-off household.

• Linear models (1 and 2). The score is the linear prediction that predict stores after regress: the fitted log expenditure X_hβ̂ + the district effect + the urban effect, for every household. You stored it as yhat_pmt (model 1) and yhat_pars (model 2) in step 2. Check: on the test round, yhat_pmt has mean 14.387 and standard deviation 0.475.
• Ordered probit models (3 and 4). After oprobit, predict with the probability option and ten new variable names writes the ten predicted probabilities Pr(decile = 1) … Pr(decile = 10) per household; they sum to 1 for every household. The score is the expected decile, computed with generate as a double: op_expected = 1·op1 + 2·op2 + 3·op3 + … + 10·op10 (and ops_expected from ops1 … ops10 for model 4). It lies between 1 and 10. Check: on the test round op_expected has mean 6.153, standard deviation 2.033, minimum 1.057, maximum 10.000 — the summarize line in the reference output. The most likely decile (the k with the largest probability) is not used: it piles households into the middle deciles and cannot select exactly half the population.

Step 3 — the predicted decile

Rank each score the same way the true decile was ranked. With xtile: a new variable, equal to the score, with weind as a frequency weight in square brackets, restricted by an if to train == 0, and the option nquantiles(10). The if matters: it ranks the 2025 households among themselves and leaves the new variable missing for the training rows. The result is the predicted decile of that model, 1 to 10, for every 2025 household. Make one such variable per score: four in all.

Check, model 1: tabulating the predicted decile on the test round gives 23 119 households in predicted decile 1, rising to 61 711 in predicted decile 10 — unequal in households, equal in people (about 28.8 million each), which is what the weight does.

The rule selects the households with predicted decile 1 to 5: the predicted bottom half of the population.

Step 4 — count people, not households

"People" means the sum of weind over a set of 2025 households. In Stata, summarize weind with the right if condition and the meanonly option leaves that sum in r(sum); store it in a scalar before the next command overwrites it. For each model compute:

• S = people selected: the sum of weind over 2025 households with predicted decile ≤ 5. Model 1: S = 143 886 559 — the same as E to within a few thousand, because both are five deciles of the same 287.8 million people; the small gap is the rounding of the decile cut.
• W = people wrongly selected: the sum of weind over 2025 households with predicted decile ≤ 5 and true decile > 5. Model 1: W = 35 325 508.

Step 5 — the error

The inclusion error: of the people the rule selects, the share that are not eligible. Model 1: 100 × 35 325 508 / 143 886 559 = 24.55 per cent. Store it as a scalar, one per model (ie_pmt, ie_pars, ie_op, ie_ops).

The exclusion error — of the eligible people, the share the rule misses — is the same number, so it is not computed separately. The selected are the bottom five predicted deciles and the eligible the bottom five true deciles of the same population, so the two groups are the same size (S ≈ E); every person wrongly selected takes the place of one eligible person who is wrongly missed, so the number missed equals W. Check it once if you like: the sum of weind over households with predicted decile > 5 and true decile ≤ 5 is 35 326 865 for model 1, and 100 × 35 326 865 / 143 887 916 = 24.55 — equal to the second decimal.

What the number means: the inclusion error is what a programme pays for people who should not be in it; the equal exclusion error is what the poor lose when they are missed.

Step 6 — the table

Finish with one table on the 2025 test round, printed with display: a title, one header line, then one line per model in the order 1 to 4 — the model's name in a left-aligned column and its IE₅₀ (= EE₅₀) in a right-aligned column with two decimals. The reference output shows the exact layout and the four values: 24.55, 21.05, 24.52, 21.09.

reference_output.txt is what the four completed do files print, one after another, when run on this package. It is the answer: every number in your log must equal it. It is ordered as your work is ordered, so a mismatch tells you where to look.

Two row counts you will meet early: the 2022 household file has 339 584 rows and its individual file 1 237 768. The pooled training sample has 1 024 865 households; the 2025 test round has 343 460.

• Before you write a round, read its questionnaire and its metadata workbook (section 2.1). A field name copied from another round that no longer exists does not always produce an error; a code that changed produces wrong numbers with no error at all.
• Write every round and every model out in full, as the given blocks are. Do not use a loop over rounds or models, and do not write a program: the reader of your file, including you in six months, should see which variable belongs to which round.
• Use forward slashes in every path. Stata accepts them on Windows, and a backslash immediately before a macro name is an escape character that fails without a message.
• Keep the given variable names, file names and globals. Steps read one another's files.
• Run the four files in order, one at a time, after editing the path lines at the top of each. Runtime on a four-core machine: about two minutes for the import, four for the build, four for the linear models, six for the ordered probits.
• Submit as a group, through the form in section 10: the four completed do files and the four logs they write (s0_import.log … s3_pmt_ordered.log in the work folder), in one zip.

Two questions to answer in a short note with your submission: model 2 has a higher R² than model 1 — why might a programme still prefer model 1? And modelling the decile directly (models 3 and 4) does not target better than modelling expenditure and cutting it afterwards — what information does the ordered model throw away before it starts?

This is a group assignment, run as a competition. Form groups of five. The goal of every group is the same: four do files that run from start to end and reproduce the reference output exactly. What differs is when a group gets there.

• Groups of five students. One submission per group; every member's name and NPM go on the form.
• How the mark is set. A group whose four logs reproduce the reference output exactly scores at least 70, and the first group to do so scores 100. The others fall between 70 and 100 according to when their submission arrives. A submission that does not replicate exactly scores below 70, by how far the do files get.
• The time that counts is the server time when the form is sent, recorded automatically and shown on the receipt and the board. A group may submit again; the time of its first replicating submission counts.
• What to submit: one zip with the four completed do files and the four logs they wrote, unedited.

The submission form and the live board are on the assignment site: sdgcenter.unpad.ac.id/aay/miktrik/pmtprobit.html.

• The assignment package: four do files, the reference output (txt and PDF), and this page zip
• Reference output: the Stata logs of all four do files, the four models and the inclusion error table PDF
• The reference output as plain text txt
• The data package: four rounds of synthetic Susenas in dBase format zip · 323 MB
• Optional: the Stata program that generates the synthetic data zip · 29 MB

The two large files are on the repository’s release page; the rest are in materials/stata.