Stata Programming Tips of Writting Loops

Stata Programming 2 1 Tips on writing loops Now, lets see some tips on how to write a loop. The basic function of a loop is to repeat a command following some pattern. To start off, you need to figure out what that pattern is. One way of doing that is writing out the first couple of runs of the loop. Once you identify the pattern, think of the best type of loop for the task. The “best type of loop” depends on the context and on your taste. It may be the loop that requires the least amount of typing, or the loop that is the easiest to interpret. That’s entirely up to you, the important thing is to be sure you are getting the correct result. First, lets start by writing out the commands. Since this is a pretty short loop, we’ll write out each command. I call this “spelling out” the loop (see next page). Look at the code for a minute, and try to find out patterns in it. • First lets notice what part gets repeated without change. We are using the replace command on the same variable many times. The part that doesn’t change is ”replace line =” and ”if id==”. In almost every case, things that don’t change will go in the loop exactly as they are. • Two things change: the number of the chorus macro and the id number. – chorus0, chorus1, ..., chorus5. It seems that one of our counters could be a sequence from 0 to 5 with increments of 1 (0, 1, 2, 3, 4, 5). Lets call this counter i. – This sequence appears twice. – id==25, id==26, id==27, ... ,id==36. It seems there’s a second counter, which could be a sequence from 25 to 36 with increments of 1 (25, 26, 27, ... , 34, 35, 36). Lets call this counter j. • There is another pattern: the difference between i and j is 25 for the first five lines of code (the first time we loop over all the values i), and then it becomes 31. This means we may want to run two separate loops, or that we will need to introduce some sort of change in the loop for the second run. We’ll try that first, then we’ll look at them and see if we can write them as one loop. 1 Please contact me with any comments (typos, errors, unclear stuff, or other suggestions on how to improve these notes) at mbarron4 [at] ucsc 1 do-file use "week4/doubletrouble.dta" g id = _n local local local local local local chorus0 chorus1 chorus2 chorus3 chorus4 chorus5 "" "[Chorus]" "Double trouble--that?s what my friends all call me" "(Double trouble)" "I said, double trouble" "T-R-O-U-B-L-E" replace replace replace replace replace replace line line line line line line = = = = = = "‘chorus0’" "‘chorus1’" "‘chorus2’" "‘chorus3’" "‘chorus4’" "‘chorus5’" if if if if if if id==25 id==26 id==27 id==28 id==29 id==30 replace replace replace replace replace replace line line line line line line = = = = = = "‘chorus0’" "‘chorus1’" "‘chorus2’" "‘chorus3’" "‘chorus4’" "‘chorus5’" if if if if if if id==31 id==32 id==33 id==34 id==35 id==36 browse Its useful to begin with the loop for the first six lines of code. You could choose to use a while loop for this. do-file local i = 0 while ‘i’<= 5 { local j = ‘i’ + 25 replace line = "‘chorus‘i’’" if id==‘i’ local i = ‘i’ + 1 } The code above adds the chorus for the first six lines of code. In this particular case, a forvalues loop provides more efficient (= less typing, without too much complexity) way of doing the same task. In particular, in a forvalues loop you don’t need to update the counter at the end of each run of the loop. We’ll use two for values loops, one for each half of the lines. 2 do-file forvalues i = 0/5 { local j = ‘i’ + 25 replace line = "‘chorus‘i’’" if id==‘j’ } forvalues i = 0/5 { local j = ‘i’ + 31 replace line = "‘chorus‘i’’" if id==‘j’ } Now, notice that both loops look almost identical, but in the first case the difference between i and j is 25 and in the second case its 31. We can write both loops in one, as follows: do-file forvalues i = 0/5 { local j = ‘i’ + 25 local k = ‘i’ + 31 replace line = "‘chorus‘i’’" if id==‘j’ | id==‘k’ } I know that this exercise may seem tedious, especially since for a loop that replaces only 12 lines of code. Sometimes you’ll want to just write down the commands instead of using loops, and that could be fine in some cases, but if you want to get familiar with loops, you should try and invest some time in learning how to make them save you time. There’s a learning curve, so you may be slow at first, but with practice writing the loop you need will become easier and easier. 2 Accessing Stored Results When you run a command, Stata not only prints out some output, but it also stores some values in memory as scalars or matrices. We can access these values with return list (after non-estimation commands) and eretrun list (after estimation commands) Lets see an example. 2.1 return list 3 output window . summarize price Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------price | 74 6165.257 2949.496 3291 15906 . ret list scalars: r(N) r(sum_w) r(mean) r(Var) r(sd) r(min) r(max) r(sum) = = = = = = = = 74 74 6165.256756756757 8699525.974268788 2949.495884768919 3291 15906 456229 . display r(mean) 6165.2568 . display r(kurtosis) . Running a command with different options may produce a different return list. 4 output window . summarize price, d Price ------------------------------------------------------------Percentiles Smallest 1% 3291 3291 5% 3748 3299 10% 3895 3667 Obs 74 25% 4195 3748 Sum of Wgt. 74 50% 75% 90% 95% 99% 5006.5 6342 11385 13466 15906 Largest 13466 13594 14500 15906 Mean Std. Dev. 6165.257 2949.496 Variance Skewness Kurtosis 8699526 1.653434 4.819188 . ret list scalars: r(N) r(sum_w) r(mean) r(Var) r(sd) r(skewness) r(kurtosis) r(sum) r(min) r(max) r(p1) r(p5) r(p10) r(p25) r(p50) r(p75) r(p90) r(p95) r(p99) = = = = = = = = = = = = = = = = = = = 74 74 6165.256756756757 8699525.97426879 2949.495884768919 1.653433511704859 4.819187528464004 456229 3291 15906 3291 3748 3895 4195 5006.5 6342 11385 13466 15906 . display r(mean) 6165.2568 . display r(kurtosis) 4.8191875 5 These values are stored in memory until you issue another command. For instance, output window . summarize mpg Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------mpg | 74 21.2973 5.785503 12 41 . ret list scalars: r(N) r(sum_w) r(mean) r(Var) r(sd) r(min) r(max) r(sum) = = = = = = = = 74 74 21.2972972972973 33.47204738985561 5.785503209735141 12 41 1576 . display r(mean) 21.297297 . display r(kurtosis) . If we need to use the values at a later stage (e.g. in a loop), we can store them, typically as scalars or macros. Once you store them as scalars, they remain in memory until you drop them (for instance, with clear all). 6 output window . sum price Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------price | 74 6165.257 2949.496 3291 15906 . scalar price_mn = r(mean) . sum mpg Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------mpg | 74 21.2973 5.785503 12 41 . scalar mpg_mn = r(mean) . display price_mn 6165.2568 . display mpg_mn 21.297297 2.2 ereturn list The ereturn list command works pretty much the same as return list. The only difference is that instead of r(name) we refer to the results as e(name). Let’s see an example: output window . regress mpg price weight Source | SS df MS -------------+-----------------------------Model | 1595.93249 2 797.966246 Residual | 847.526967 71 11.9369995 -------------+-----------------------------Total | 2443.45946 73 33.4720474 Number of obs F( 2, 71) Prob > F R-squared Adj R-squared Root MSE = = = = = = 74 66.85 0.0000 0.6531 0.6434 3.455 -----------------------------------------------------------------------------mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------price | -.0000935 .0001627 -0.57 0.567 -.000418 .0002309 weight | -.0058175 .0006175 -9.42 0.000 -.0070489 -.0045862 _cons | 39.43966 1.621563 24.32 0.000 36.20635 42.67296 ------------------------------------------------------------------------------ 7 output window . ereturn list scalars: e(N) e(df_m) e(df_r) e(F) e(r2) e(rmse) e(mss) e(rss) e(r2_a) e(ll) e(ll_0) e(rank) = = = = = = = = = = = = e(cmdline) e(title) e(marginsok) e(vce) e(depvar) e(cmd) e(properties) e(predict) e(model) e(estat_cmd) : : : : : : : : : : 74 2 71 66.84814256414501 .6531446579233134 3.454996314099513 1595.932492798133 847.5269666613265 .6433740849070687 -195.2169813478502 -234.3943376482347 3 macros: "regress mpg price weight" "Linear regression" "XB default" "ols" "mpg" "regress" "b V" "regres_p" "ols" "regress_estat" matrices: e(b) : e(V) : 1 x 3 3 x 3 functions: e(sample) We can access them by their names. Just remember to use e(name) instead of r(name) 8 output window . display e(F) 66.848143 . display r(F) . . matrix list e(b) e(b)[1,3] price y1 -.00009351 weight -.00581754 _cons 39.439656 . matrix list e(V) symmetric e(V)[3,3] price weight price 2.648e-08 weight -5.412e-08 3.813e-07 _cons 1.776e-07 -.00081778 3 _cons 2.6294664 Aside: Matrices in Stata In this course we will not cover matrix algebra in Stata. However, since Stata stores some useful results in matrices, you should learn how to access those results. First, it is easier to create a copy of the matrix and working with that instead of using the e(b) or e(V) syntax. The following code shows how to create a copy of the e( ) matrix and list its contents: 9 output window . matrix B = e(b) . matrix list B B[1,3] y1 price -.00009351 weight -.00581754 _cons 39.439656 . matrix V = e(V) . matrix list V symmetric V[3,3] price price 2.648e-08 weight -5.412e-08 _cons 1.776e-07 weight _cons 3.813e-07 -.00081778 2.6294664 To refer to a particular element in the matrix, we use the row and column number (in that order), inside square brackets. output window * Display the element in row 1, column 2 . display B[1,2] -.00581754 * Since B has only 1 row, this will produce a missing value: . display B[2,1] . We can use these in loops: do-file forvalues i = 1/3 { display B[1,‘i’] scalar B_‘i’ = B[1,‘i’] } 10 output window -.00009351 -.00581754 39.439656 . . scalar list B_3 = 39.439656 B_2 = -.00581754 B_1 = -.00009351 4 Application: Hypothesis Testing After running a regression, you can use the results to test hypotheses about the coefficients in your model. The lincom command provides an easy way of doing this. First, we need to estimate the model output window . reg wage hours tenure ttl_exp age Source | SS df MS -------------+-----------------------------Model | 6218.16546 4 1554.54137 Residual | 67802.6451 2222 30.5142417 -------------+-----------------------------Total | 74020.8106 2226 33.2528349 Number of obs F( 4, 2222) Prob > F R-squared Adj R-squared Root MSE = = = = = = 2227 50.94 0.0000 0.0840 0.0824 5.524 -----------------------------------------------------------------------------wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------hours | .0545902 .0115028 4.75 0.000 .0320329 .0771474 tenure | .0366771 .0260101 1.41 0.159 -.0143296 .0876837 ttl_exp | .2848398 .0317393 8.97 0.000 .222598 .3470817 age | -.1211649 .0385944 -3.14 0.002 -.1968498 -.04548 _cons | 6.713919 1.571873 4.27 0.000 3.631426 9.796413 ------------------------------------------------------------------------------ Now we can test, for example, if the coefficient on hours is statistically equal to zero. The null hypothesis is that the coefficient is zero, and the resulting p-value is the p-value associated with that hypothesis. If the p-value less than 0.05, we reject the hypothesis at the 95% of confidence. 11 output window . lincom hours ( 1) hours = 0 /* <= This is the null hypothesis */ -----------------------------------------------------------------------------wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------(1) | .0545902 .0115028 4.75 0.000 .0320329 .0771474 ------------------------------------------------------------------------------ We can do the same with tenure. output window . lincom tenure ( 1) tenure = 0 -----------------------------------------------------------------------------wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------(1) | .0366771 .0260101 1.41 0.159 -.0143296 .0876837 ------------------------------------------------------------------------------ These first two examples are pretty useless, since the t-tests in the regression table gives the results to those tests. The main use of lincom is to test linear combinations of the coefficients. For example, we can test if the coefficient associated to “hours” is equal to the coefficient associated to “tenure”. When we test linear combinations, we implicitly test if the expression is equal to zero. So instead of testing directly if a coefficient is equal to the other, we test if the difference between them is equal to zero. output window . lincom hours - tenure ( 1) hours - tenure = 0 -----------------------------------------------------------------------------wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------(1) | .0179131 .0288196 0.62 0.534 -.0386031 .0744293 ------------------------------------------------------------------------------ Note that despite the coefficient on hours is significantly different than zero and the coefficient on tenure is not statistically different than zero, these coefficients are not statistically different from each other. Boom. We may want to test for specific values, that are derived from a theoretical model, or that are useful benchmarks. For instance, we may want to test if the elasticity of wage to schooling is equal to one. We can test any linear combination of the coefficients, but only linear combinations. The nlcom command allows testing hypotheses that involve nonlinear 12 combinations of estimators. output window . g lnschooling = ln(grade) (4 missing values generated) . g lnschooling = ln(grade) . reg lnwage lnschooling hours tenure ttl_exp age Source | SS df MS -------------+-----------------------------Model | 205.206226 5 41.0412452 Residual | 522.597137 2217 .23572266 -------------+-----------------------------Total | 727.803363 2222 .327544268 Number of obs F( 5, 2217) Prob > F R-squared Adj R-squared Root MSE = = = = = = 2223 174.11 0.0000 0.2820 0.2803 .48551 -----------------------------------------------------------------------------lnwage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------lnschooling | .9748042 .052895 18.43 0.000 .8710753 1.078533 hours | .0058235 .0010116 5.76 0.000 .0038397 .0078073 tenure | .0117327 .002288 5.13 0.000 .0072457 .0162196 ttl_exp | .02857 .0028258 10.11 0.000 .0230284 .0341115 age | -.0080064 .0034015 -2.35 0.019 -.0146769 -.0013359 _cons | -.9480892 .1944988 -4.87 0.000 -1.329508 -.5666703 ------------------------------------------------------------------------------ output window . lincom lnschooling - 1 ( 1) lnschooling = 1 -----------------------------------------------------------------------------lnwage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------(1) | -.0251958 .052895 -0.48 0.634 -.1289247 .0785331 ------------------------------------------------------------------------------ 13 output window . lincom 2 * hours - (1/3) * tenure ( 1) 2*hours - .3333333*tenure = 0 -----------------------------------------------------------------------------lnwage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------(1) | .0077361 .0021878 3.54 0.000 .0034459 .0120264 -----------------------------------------------------------------------------. lincom 2 * hours - (1/3) * tenure^2 not possible with test r(131); In this week’s assignment you’ll be asked to use test, another command that allows testing linear hypotheses. 5 Adding Scalars to Tables A very useful feature of the return and ereturn values is that we can add them to regression tables with outreg2 do-file reg wage hours tenure ttl_exp age outreg2 using table1, replace /// adds("F-stat", e(F)) do-file reg wage hours tenure ttl_exp age summarize wage if e(sample) estadd scalar mn_wage = r(mean) outreg2 using table1, append /// adds("F-stat", e(F),"Mean Wage", e(mn_wage)) 14 do-file reg wage hours tenure ttl_exp age summarize wage if e(sample) estadd scalar mn_wage = r(mean) lincom hours - tenure estadd scalar diff = r(estimate) estadd scalar se_diff = r(se) estadd scalar t_diff = r(estimate) / r(se) outreg2 using table1, append /// adds("F-stat", e(F), "Mean Wage", e(mn_wage), /// "H0: hours = tenure", e(diff), /// "Difference SE", e(se_diff), "t-stat", e(t_diff)) Here is the resulting table: VARIABLES hours tenure ttl exp age Constant (1) wage (2) wage 0.0546*** 0.0546*** (0.0115) (0.0115) 0.0367 0.0367 (0.0260) (0.0260) 0.285*** 0.285*** (0.0317) (0.0317) -0.121*** -0.121*** (0.0386) (0.0386) 6.714*** 6.714*** (1.572) (1.572) Observations 2,227 2,227 R-squared 0.084 0.084 F-stat 50.94 50.94 Mean Wage 7.800 H0: hours = tenure Difference SE t-stat Standard errors in parentheses *** p<0.01, ** p<0.05, * p<0.1 15 (3) wage 0.0546*** (0.0115) 0.0367 (0.0260) 0.285*** (0.0317) -0.121*** (0.0386) 6.714*** (1.572) 2,227 0.084 50.94 7.800 0.0179 0.0288 0.622