### Archive

Archive for the ‘Programming’ Category

## Programming an estimation command in Stata: Allowing for sample restrictions and factor variables

I modify the ordinary least-squares (OLS) command discussed in Programming an estimation command in Stata: A better OLS command to allow for sample restrictions, to handle missing values, to allow for factor variables, and to deal with perfectly collinear variables.

Sample restrictions

The myregress4 command described in Programming an estimation command in Stata: A better OLS command has the syntax

myregress4 depvar [indepvars]

where the indepvars may be time-series operated variables. myregress5 allows for sample restrictions and missing values. It has the syntax

myregress5 depvar [indepvars] [if] [in]

A user may optionally specify an if expression or an in range to restrict the sample. I also make myregress5 handle missing values in the user-specified variables.

*! version 5.0.0  22Nov2015
program define myregress5, eclass
version 14

syntax varlist(numeric ts) [if] [in]
marksample touse

gettoken depvar : varlist

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist' if touse'
local p : word count varlist'
local p = p' + 1
matrix xpx'                = zpz'[2..p', 2..p']
matrix xpy'                = zpz'[2..p', 1]
matrix xpxi'               = syminv(xpx')
matrix b'                  = (xpxi'*xpy')'
quietly matrix score double xbhat' = b' if touse'
quietly generate double res'       = (depvar' - xbhat') if touse'
quietly generate double res2'      = (res')^2 if touse'
quietly summarize res2' if touse' , meanonly
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(p'-1))
matrix V'                  = s2'*xpxi'
ereturn post b' V', esample(touse')
ereturn scalar           N  = N'
ereturn local         cmd   "myregress5"
ereturn display
end


The syntax command in line 5 specifies that a user may optionally restrict the sample by specifying an if expression or in range. When the user specifies an if expression, syntax puts it into the local macro if; otherwise, the local macro if is empty. When the user specifies an in range, syntax puts it into the local macro in; otherwise, the local macro in is empty.

We could use the local macros if and in to handle user-specified sample restrictions, but these do not account for missing values in the user-specified variables. The marksample command in line 6 creates a local macro named touse, which contains the name of a temporary variable that is a sample-identification variable. Each observation in the sample-identification variable is either one or zero. It is one if the observation is included in the sample. It is zero if the observation is excluded from the sample. An observation can be excluded by a user-specified if expression, by a user-specified in range, or because there is a missing value in one of the user-specified variables.

Lines 20–23 use the sample-identification variable contained in the local macro touse to enforce these sample restrictions on the OLS calculations.

Line 28 posts the sample-identification variable into e(sample), which is one if the observation was included in the estimation sample and it is zero if the observation was excluded from the estimation sample.

Line 29 stores the number of observations in the sample in e(N).

Example 1 illustrates that myregress5 runs the requested regression on the sample that respects the missing values in rep78 and accounts for an if expression.

Example 1: myregress5 with missing values and an if expression

. sysuse auto
(1978 Automobile Data)

. count if !missing(rep78)
69

. count if !missing(rep78) & mpg < 30
62

. myregress5 price mpg trunk rep78 if mpg < 30
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -376.8591   107.4289    -3.51   0.000    -587.4159   -166.3024
trunk |  -36.39376   102.2139    -0.36   0.722    -236.7294    163.9418
rep78 |   556.3029   378.1101     1.47   0.141    -184.7793    1297.385
_cons |    12569.5   3455.556     3.64   0.000     5796.735    19342.27
------------------------------------------------------------------------------

. ereturn list

scalars:
e(N) =  62

macros:
e(cmd) : "myregress5"
e(properties) : "b V"

matrices:
e(b) :  1 x 4
e(V) :  4 x 4

functions:
e(sample)


Allowing for factor variables

Example 1 includes the number of repairs as a continuous variable, but it might be better treated as a discrete factor. myregress6 accepts factor variables. Factor-variable lists usually imply variable lists that contain perfectly collinear variables, so myregress6 also handles
perfectly collinear variables.

*! version 6.0.0  22Nov2015
program define myregress6, eclass
version 14

syntax varlist(numeric ts fv) [if] [in]
marksample touse

gettoken depvar : varlist
_fv_check_depvar depvar'

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist' if touse'
local p                    = colsof(zpz')
matrix xpx'               = zpz'[2..p', 2..p']
matrix xpy'               = zpz'[2..p', 1]
matrix xpxi'              = syminv(xpx')
local k                    = p' - diag0cnt(xpxi') - 1
matrix b'                 = (xpxi'*xpy')'
quietly matrix score double xbhat' = b' if touse'
quietly generate double res'       = (depvar' - xbhat') if touse'
quietly generate double res2'      = (res')^2 if touse'
quietly summarize res2' if touse' , meanonly
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(k'))
matrix V'                  = s2'*xpxi'
ereturn post b' V', esample(touse')
ereturn scalar N            = N'
ereturn scalar rank         = k'
ereturn local  cmd             "myregress6"
ereturn display
end


The fv in the parentheses after varlist in the syntax command in line 5 modifies the varlist to accept factor variables. Any specified factor variables are stored in the local macro varlist in a canonical form.

Estimation commands do not allow the dependent variable to be a factor variable. The _fv_check_depvar command in line 9 will exit with an error if the local macro depvar contains a factor variable.

Line 15 stores the number of columns in the matrix formed by matrix accum in the local macro p. Line 19 stores the number of linearly independent columns in the local macro k. This calculation uses diag0cnt() to account for the perfectly collinear variables that were dropped. (Each dropped variable puts a zero on the diagonal of the generalized inverse calculated by syminv() and diag0cnt() returns the number of zeros on the diagonal.)

Line 31 stores the number of linearly independent variables in e(rank) for postestimation commands.

Now, I use myregress6 to include rep78 as a factor-operated variable. The base category is dropped because we included a constant term.

Example 2: myregress6 with a factor-operated variable

. myregress6 price mpg trunk i.rep78
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -262.7053   73.49434    -3.57   0.000    -406.7516   -118.6591
trunk |   41.75706    93.9671     0.44   0.657    -142.4151    225.9292
|
rep78 |
2  |   654.7905   2136.246     0.31   0.759    -3532.175    4841.756
3  |   1170.606   2001.739     0.58   0.559     -2752.73    5093.941
4  |   1473.352   2017.138     0.73   0.465    -2480.167     5426.87
5  |   2896.888   2121.206     1.37   0.172    -1260.599    7054.375
|
_cons |   9726.377   2790.009     3.49   0.000      4258.06    15194.69
------------------------------------------------------------------------------


Done and undone

I modified the OLS command discussed in Programming an estimation command in Stata: A better OLS command to allow for sample restrictions, to handle missing values, to allow for factor variables, and to deal with perfectly collinear variables. In the next post, I show how to allow options for robust standard errors and to suppress the constant term.

Categories: Programming Tags:

## Programming an estimation command in Stata: A better OLS command

I use the syntax command to improve the command that implements the ordinary least-squares (OLS) estimator that I discussed in Programming an estimation command in Stata: A first command for OLS. I show how to require that all variables be numeric variables and how to make the command accept time-series operated variables.

Stata syntax and the syntax command

The myregress2 command described in Programming an estimation command in Stata: A first command for OLS has the syntax

myregress2 depvar [indepvars]

This syntax requires that the dependent variable be specified because depvar is not enclosed in square brackets. The independent variables are optional because indepvars is enclosed in square brackets. Type

for an introduction to reading Stata syntax diagrams.

This syntax is implemented by the syntax command in line 5 of myregress2.ado, which I discussed at length in Programming an estimation command in Stata: A first command for OLS. The user must specify a list of variable names because varlist is not enclosed in square brackets. The syntax of the syntax command follows the rules of a syntax diagram.

*! version 2.0.0  26Oct2015
program define myregress2, eclass
version 14

syntax varlist
gettoken depvar : varlist

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist'
local p : word count varlist'
local p = p' + 1
matrix xpx'                = zpz'[2..p', 2..p']
matrix xpy'                = zpz'[2..p', 1]
matrix xpxi'               = syminv(xpx')
matrix b'                  = (xpxi'*xpy')'
quietly matrix score double xbhat' = b'
quietly generate double res'       = (depvar' - xbhat')
quietly generate double res2'      = (res')^2
quietly summarize res2'
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(p'-1))
matrix V'                  = s2'*xpxi'
ereturn post b' V'
ereturn local         cmd   "myregress2"
ereturn display
end


Example 1 illustrates that myregress2 runs the requested regression when I specify a varlist.

Example 1: myregress2 with specified variables

. sysuse auto
(1978 Automobile Data)

. myregress2 price mpg trunk
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   65.59262    -3.36   0.001    -348.7241    -91.6057
trunk |   43.55851   88.71884     0.49   0.623    -130.3272    217.4442
_cons |   10254.95   2349.084     4.37   0.000      5650.83    14859.07
------------------------------------------------------------------------------


Example 2 illustrates that the syntax command displays an error message and stops execution when I do not specify a varlist. I use set trace on to see each line of code and the output it produces.

Example 2: myregress2 without a varlist

. set trace on

. myregress2
--------------------------------------------------------- begin myregress2 --
- version 14
- syntax varlist
varlist required
----------------------------------------------------------- end myregress2 --
r(100);


Example 3 illustrates that the syntax command is checking that the specified variables are in the current dataset. syntax throws an error because DoesNotExist is not a variable in the current dataset.

Example 3: myregress2 with a variable not in this dataset

. set trace on

. myregress2 price mpg trunk DoesNotExist
--------------------------------------------------------- begin myregress2 --
- version 14
- syntax varlist
----------------------------------------------------------- end myregress2 --
r(111);

end of do-file

r(111);


Because the syntax command on line 5 is not restricting the specified variables to be numeric, I get the no observations error in example 4 instead of an error indicating the actual problem, which is the string variable make.

Example 4: myregress2 with a string variable

. describe make

storage   display    value
variable name   type    format     label      variable label
-------------------------------------------------------------------------------
make            str18   %-18s                 Make and Model

. myregress2 price mpg trunk make
no observations
r(2000);

end of do-file

r(2000);


On line 5 of myregress3, I modify varlist to only accept numeric variables This change produces a more informative error message when I try to include a string variable in the regression.

*! version 3.0.0  30Oct2015
program define myregress3, eclass
version 14

syntax varlist(numeric)
gettoken depvar : varlist

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist'
local p : word count varlist'
local p = p' + 1
matrix xpx'                = zpz'[2..p', 2..p']
matrix xpy'                = zpz'[2..p', 1]
matrix xpxi'               = syminv(xpx')
matrix b'                  = (xpxi'*xpy')'
quietly matrix score double xbhat' = b'
quietly generate double res'       = (depvar' - xbhat')
quietly generate double res2'      = (res')^2
quietly summarize res2'
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(p'-1))
matrix V'                  = s2'*xpxi'
ereturn post b' V'
ereturn local         cmd   "myregress3"
ereturn display
end


Example 5: myregress3 with a string variable

. set trace on

. myregress3 price mpg trunk make
--------------------------------------------------------- begin myregress3 --
- version 14
- syntax varlist(numeric)
string variables not allowed in varlist;
make is a string variable
----------------------------------------------------------- end myregress3 --
r(109);

end of do-file

r(109);


On line 5 of myregress4, I modify the varlist to accept time-series (ts) variables. The syntax command puts time-series variables in a canonical form that is stored in the local macro varlist, as illustrated in the display on line 6, whose output appears in example 6.

*! version 4.0.0  31Oct2015
program define myregress4, eclass
version 14

syntax varlist(numeric ts)
display "varlist is varlist'"
gettoken depvar : varlist

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist'
local p : word count varlist'
local p = p' + 1
matrix xpx'                = zpz'[2..p', 2..p']
matrix xpy'                = zpz'[2..p', 1]
matrix xpxi'               = syminv(xpx')
matrix b'                  = (xpxi'*xpy')'
quietly matrix score double xbhat' = b'
quietly generate double res'       = (depvar' - xbhat')
quietly generate double res2'      = (res')^2
quietly summarize res2'
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(p'-1))
matrix V'                  = s2'*xpxi'
ereturn post b' V'
ereturn local         cmd   "myregress4"
ereturn display
end


Example 6: myregress4 with time-series variables

. sysuse gnp96

. myregress4  L(0/3).gnp
varlist is gnp96 L.gnp96 L2.gnp96 L3.gnp96
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
gnp96 |
L1. |   1.277086   .0860652    14.84   0.000     1.108402    1.445771
L2. |   -.135549   .1407719    -0.96   0.336    -.4114568    .1403588
L3. |  -.1368326   .0871645    -1.57   0.116    -.3076719    .0340067
|
_cons |   -2.94825   14.36785    -0.21   0.837    -31.10871    25.21221
------------------------------------------------------------------------------


Done and undone

I used the syntax command to improve how myregress2 handles the variables specified by the user. I showed how to require that all variables be numeric variables and how to make the command accept time-series operated variables. In the next post, I show how to make the command allow for sample restrictions, how to handle missing values, how to allow for factor-operated variables, and how to deal with perfectly collinear variables.

Categories: Programming Tags:

## Programming an estimation command in Stata: A first command for OLS

Programming an estimation command in Stata: A first ado command, in which I introduced some ado-programming concepts. Although I introduce some local macro tricks that I use all the time, I also build on Programing an estimation command in Stata: Where to store your stuff.

Local macro tricks

I use lots of local macro tricks in my ado-files. In this section, I illustrate ones that I use in the commands that I develop in this post. In every ado-file that I write, I ask questions about lists of variable names stored in local macros. I frequently use the extended macro functions and the gettoken command to ask these questions and store the results in a local macro.

The syntax for storing the result of an extended macro function in a local macro is

local localname : extended_fcn

Below, I use the extended macro function word count to count the number of elements in the list and store the result in the local macro count.

Example 1: Storing and extracting the result of an extended macro function

. local count : word count a b c

. display "count contains count'"
count contains 3


There are many extended macro functions, but I illustrate just the one I use in this post; type help extended fcn for a complete list.

A token is an element in a list. I frequently use the gettoken command to split lists apart. The gettoken command has the syntax

gettoken localname1 [localname2] : localname3

gettoken stores the first token in the list stored in the local macro localname3 into the local macro localname1. If the optional localname2 is specified, the remaining tokens are stored in the local macro localname2.

I use gettoken to store the first token stored in mylist into the local macro first, whose contents I subsequently extract and display.

Example 2: Using gettoken to store first token only

. local mylist y x1 x2

. display "mylist contains mylist'"
mylist contains y x1 x2

. gettoken first : mylist

. display "first contains first'"
first contains y


Now, I use gettoken to store the first token stored in mylist into the local macro first and the remaining tokens into the local macro left. I subsequently extract and display the contents of first and left.

Example 3: Using gettoken to store first and remaining tokens

. gettoken first left: mylist

. display "first contains first'"
first contains y

. display "left  contains left'"
left  contains  x1 x2


I frequently want to increase the value of a local macro by some fixed amount, say, $$3$$. I now illustrate a solution that I use.

Example 4: Local macro update

. local p = 1

. local p = p' + 3

. display "p is now p'"
p is now 4


When the update value, also known as the increment value, is $$1$$, we can use the increment operator, as below:

Example 5: Local macro update

. local p = 1

. local ++p

. display "p is now p'"
p is now 2


A first version of myregress

The code in myregress1 implements a version of the OLS formulas. I use myregress1 in example 6. Below example 6, I discuss the code and the output.

*! version 1.0.0  23Oct2015
program define myregress1, eclass
version 14

syntax varlist
display "The syntax command puts the variables specified by the "
display "    user into the local macro varlist"
display "    varlist contains varlist'"

gettoken depvar : varlist
display "The dependent variable is depvar'"

matrix accum zpz = varlist'
display "matrix accum forms Z'Z"
matrix list zpz

local p : word count varlist'
local p = p' + 1

matrix xpx                = zpz[2..p', 2..p']
matrix xpy                = zpz[2..p', 1]
matrix xpxi               = syminv(xpx)
matrix b                  = (xpxi*xpy)'

matrix score double xbhat = b
generate double res       = (depvar' - xbhat)
generate double res2      = res^2
summarize res2
local N                   = r(N)
local sum                 = r(sum)
local s2                  = sum'/(N'-(p'-1))
matrix V                  = s2'*xpxi
ereturn post b V
ereturn local         cmd   "myregress1"
ereturn display
end


Example 6: myregress1 output

. sysuse auto
(1978 Automobile Data)

. myregress1 price mpg trunk
The syntax command puts the variables specified by the
user into the local macro varlist
varlist contains price mpg trunk
The dependent variable is price
(obs=74)
matrix accum forms Z'Z

symmetric zpz[4,4]
price        mpg      trunk      _cons
price  3.448e+09
mpg    9132716      36008
trunk    6565725      20630      15340
_cons     456229       1576       1018         74

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
res2 |         74     6674851    1.30e+07   11.24372   9.43e+07
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   65.59262    -3.36   0.001    -348.7241    -91.6057
trunk |   43.55851   88.71884     0.49   0.623    -130.3272    217.4442
_cons |   10254.95   2349.084     4.37   0.000      5650.83    14859.07
------------------------------------------------------------------------------


Here are my comments on the code and the output in example 6.

• Line 2 specifies that myregress1 is an e-class command that stores its results in e().
• Lines 5–8 illustrate that the syntax command stores the names of the variables specified by the user in the local macro varlist. This behavior is also illustrated in example 6.
• Line 10 uses the gettoken command to store the first variable name stored in the local macro varlist in the local macro depvar. Line 11 displays this name and the usage is illustrated in example 6.
• Line 13 uses matrix accum to put $$(\Xb’\Xb)$$ and $$(\Xb’\yb)$$ into a Stata matrix named zpz, as discussed in Programming an estimation command in Stata: Using Stata matrix commands and functions to compute OLS objects and further illustrated in lines 14–15 and example 6.
• Line 17 stores the number of variables in the local macro varlist into the local macro p.
• Line 18 increments the local macro p by $$1$$ to account for the constant term included by matrix accum by default.
• Lines 20–23 extract $$(\Xb’\Xb)$$ and $$(\Xb’\yb)$$ from zpz and put the vector of point estimates $$\widehat{\betab}$$ into the Stata row vector b.
• Line 25 puts $$\Xb\widehat{\betab}$$ into the variable xbhat.
• Lines 26 and 27 calculate the residuals and the squared residuals, respectively.
• Lines 28–32 calculate the estimated variance-covariance matrix of the estimator (VCE) from the sum of squared residuals.
• Line 33 stores b and V into e(b) and e(V), respectively.
• Line 34 stores the name of the estimation command (myregress1) in e(cmd).
• Line 35 produces a standard Stata output table from the results in e(b) and e(V).

myregress1 contains code to help illustrate how it works, and it uses hard-coded names for global objects like Stata variables and Stata matrices. Users do not want to see the output from the illustration lines, so they must be removed. Users do not want their global Stata matrices overwritten by a command they use, which is what myregress1 would do to a matrix named zpz, xpx, xpxi, b, or V.

The code in myregress2 fixes these problems.

*! version 2.0.0  26Oct2015
program define myregress2, eclass
version 14

syntax varlist
gettoken depvar : varlist

tempname zpz xpx xpy xpxi b V
tempvar  xbhat res res2

quietly matrix accum zpz' = varlist'
local p : word count varlist'
local p = p' + 1
matrix xpx'                = zpz'[2..p', 2..p']
matrix xpy'                = zpz'[2..p', 1]
matrix xpxi'               = syminv(xpx')
matrix b'                  = (xpxi'*xpy')'
quietly matrix score double xbhat' = b'
generate double res'       = (depvar' - xbhat')
generate double res2'      = (res')^2
quietly summarize res2'
local N                     = r(N)
local sum                   = r(sum)
local s2                    = sum'/(N'-(p'-1))
matrix V'                  = s2'*xpxi'
ereturn post b' V'
ereturn local         cmd   "myregress2"
ereturn display
end

• Line 8 uses tempname to put safe names into the local macros zpz, xpx, xpy, xpxi, b, and V.
• Line 9 uses tempvar to put safe names into the local macros xbhat, res, res2.
• Lines 11, 14–18, and 25–26 use the safe names in the local macros created by tempname instead of the hard-coded names for the matrices.
• Lines 18–20 use the safe names in the local macros created by tempvar instead of the hard-coded names for the variables it creates.

The output below shows the output produced by myregress2.

Example 7: myregress2 output

. myregress2 price mpg trunk
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   65.59262    -3.36   0.001    -348.7241    -91.6057
trunk |   43.55851   88.71884     0.49   0.623    -130.3272    217.4442
_cons |   10254.95   2349.084     4.37   0.000      5650.83    14859.07
------------------------------------------------------------------------------


Done and undone

After reviewing some tricks with local macros that I use in most of the ado-files that I write, I discussed two versions of an ado-command that implements the (OLS) estimator. In the next post, I extend this command so that the user may request a robust VCE, or that the constant term be suppressed, or both.

Categories: Programming Tags:

## Programming an estimation command in Stata: Using Stata matrix commands and functions to compute OLS objects


OLS formulas

Recall that the OLS point estimates are given by

$\widehat{\betab} = \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1} \left( \sum_{i=1}^N \xb_i’y_i \right)$

where $${\bf x}_i$$ is the $$1\times k$$ vector of independent variables, $$y_i$$ is the dependent variable for each of the $$N$$ sample observations, and the model for $$y_i$$ is

$y_i = \xb_i\betab’ + \epsilon_i$

If the $$\epsilon_i$$ are independently and identically distributed, we estimate the variance-covariance matrix of the estimator (VCE) by

$\widehat{\Vb} = \widehat{s} \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1}$

where $$\widehat{s} = 1/(N-k)\sum_{i=1}^N e_i^2$$ and $$e_i=y_i-{\bf x}_i\widehat{{\boldsymbol \beta}}$$.

See Cameron and Trivedi (2005), Stock and Watson (2010), or Wooldridge (2015) for introductions to OLS.

Stata matrix implementation

I use the matrix accum command to compute the sum of the products over the observations. Typing

.  matrix accum zpz = z1 z2 z3

puts $$\left( \sum_{i=1}^N {\bf z}_i'{\bf z}_i \right)$$ into the Stata matrix zpz, where $${\bf z}_i=( {\tt z1}_i, {\tt z2}_i, {\tt z3}_i, 1)$$. The $$1$$ appears because matrix accum has included the constant term by default, like almost all estimation commands.

Below, I use matrix accum to compute $$\left( \sum_{i=1}^N {\bf z}_i'{\bf z}_i \right)$$, which contains $$\left( \sum_{i=1}^N {\bf x}_i'{\bf x}_i \right)$$ and $$\left( \sum_{i=1}^N {\bf x}_i’y_i \right)$$.

Example 1: Using matrix accum

. sysuse auto
(1978 Automobile Data)

. matrix accum zpz = price mpg trunk
(obs=74)

. matrix list zpz

symmetric zpz[4,4]
price        mpg      trunk      _cons
price  3.448e+09
mpg    9132716      36008
trunk    6565725      20630      15340
_cons     456229       1576       1018         74


Now, I extract $$\left( \sum_{i=1}^N {\bf x}_i'{\bf x}_i \right)$$ from rows 2–4 and columns 2–4 of zpz and $$\left( \sum_{i=1}^N {\bf x}_i’y_i \right)$$ from rows 2–4 and column 1 of zpz.

Example 2: Extracting submatrices

. matrix xpx       = zpz[2..4, 2..4]

. matrix xpy       = zpz[2..4, 1]

. matrix list xpx

symmetric xpx[3,3]
mpg  trunk  _cons
mpg  36008
trunk  20630  15340
_cons   1576   1018     74

. matrix list xpy

xpy[3,1]
price
mpg  9132716
trunk  6565725
_cons   456229


I now compute $$\widehat{{\boldsymbol \beta}}$$ from the matrices formed in example 2.

Example 3: Computing $$\widehat{\betab}$$

. matrix xpxi      = syminv(xpx)

. matrix b         = xpxi*xpy

. matrix list b

b[3,1]
price
mpg  -220.16488
trunk    43.55851
_cons    10254.95

. matrix b         = b'

. matrix list b

b[1,3]
mpg       trunk       _cons
price  -220.16488    43.55851    10254.95


I transposed b to make it a row vector because point estimates in Stata are stored as row vectors.

Example 3 illustrates that the Stata matrix b contains the estimated coefficients and the names of the variables on which these values are estimated coefficients. To clarify, our model is
$\Eb[{\tt price}|{\tt mpg}, {\tt trunk} ] = {\tt mpg}*\beta_{\tt mpg} + {\tt trunk}*\beta_{\tt trunk} + {\tt \_cons}$

and b contains the information that $$-220.16$$ is the estimated coefficient on mpg, that $$43.56$$ is the estimated coefficient on trunk, and that $$10254.95$$ is the estimated constant. We can compute the linear combination $$\xb_i\widehat{\betab}$$ over the observations using the information in b, because b contains both the value and the name for each coefficient.

I use matrix score to compute this linear combination for each observation, and I use generate to reiterate what this linear combination is.

Example 4: Using matrix score to compute $$\xb_i\widehat{\betab}’$$

. matrix score double xbhat1 = b

. generate     double xbhat2 = mpg*(-220.16488) + trunk*(43.55851) + 10254.95

. list xbhat1 xbhat2 in 1/4

+-----------------------+
|    xbhat1      xbhat2 |
|-----------------------|
1. | 5890.4661   5890.4663 |
2. | 6991.2905   6991.2907 |
3. | 5934.0246   5934.0248 |
4. | 6548.5884   6548.5886 |
+-----------------------+


I use the predictions for $$\Eb[{\tt price}|{\tt mpg}, {\tt trunk} ]$$ in xbhat1 to compute the residuals and the estimated VCE.

Example 5: Computing the estimated VCE

. generate double res       = (price - xbhat1)

. generate double res2      = res^2

. summarize res2

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
res2 |         74     6674851    1.30e+07   11.24372   9.43e+07

. return list

scalars:
r(N) =  74
r(sum_w) =  74
r(mean) =  6674850.504745401
r(Var) =  168983977867533.1
r(sd) =  12999383.74952956
r(min) =  11.24371634723049
r(max) =  94250157.2111593
r(sum) =  493938937.3511598

. local N                   = r(N)

. local sum                 = r(sum)

. local s2                  = sum'/(N'-3)

. matrix V                  = (s2')*xpxi


(See Programing an estimation command in Stata: Where to store your stuff for discussions of using results from r-class commands and using local macros.)

I verify that my computations for $$\widehat{\betab}$$ and the VCE match those of regress.

Example 6: Comparing against regress

. regress price mpg trunk

Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(2, 71)        =     10.14
Model |   141126459         2  70563229.4   Prob > F        =    0.0001
Residual |   493938937        71  6956886.44   R-squared       =    0.2222
Total |   635065396        73  8699525.97   Root MSE        =    2637.6

------------------------------------------------------------------------------
price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   65.59262    -3.36   0.001    -350.9529    -89.3769
trunk |   43.55851   88.71884     0.49   0.625    -133.3418    220.4589
_cons |   10254.95   2349.084     4.37   0.000      5571.01    14938.89
------------------------------------------------------------------------------

. matrix list e(b)

e(b)[1,3]
mpg       trunk       _cons
y1  -220.16488    43.55851    10254.95

. matrix list b

b[1,3]
mpg       trunk       _cons
price  -220.16488    43.55851    10254.95

. matrix list e(V)

symmetric e(V)[3,3]
mpg       trunk       _cons
mpg   4302.3924
trunk   3384.4186   7871.0326
_cons  -138187.95  -180358.85   5518194.7

. matrix list V

symmetric V[3,3]
mpg       trunk       _cons
mpg   4302.3924
trunk   3384.4186   7871.0326
_cons  -138187.95  -180358.85   5518194.7


Robust standard errors

The frequently used robust estimator of the VCE is given by

$\widehat{V}_{robust}=\frac{N}{N-k} \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1} \Mb \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1}$

where

$\Mb=\sum_{i=1}^N \widehat{e}_i^2\xb_i’\xb_i$

See Cameron and Trivedi (2005), Stock and Watson (2010), or Wooldridge (2015) for derivations and discussions.

matrix accum with weights $$\widehat{e}_i^2$$ computes the formula for $$\Mb$$. Below, I use matrix accum to compute $$\Mb$$ and $$\widehat{V}_{robust}$$

Example 7: A robust VCE

. matrix accum M    = mpg trunk [iweight=res2]
(obs=493938937.4)

. matrix V2         = (N'/(N'-3))*xpxi*M*xpxi


I now verify that my computations match those reported by regress.

Example 8: Comparing computations of robust VCE

. regress price mpg trunk, robust

Linear regression                               Number of obs     =         74
F(2, 71)          =      11.59
Prob > F          =     0.0000
R-squared         =     0.2222
Root MSE          =     2637.6

------------------------------------------------------------------------------
|               Robust
price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   72.45388    -3.04   0.003    -364.6338   -75.69595
trunk |   43.55851    71.4537     0.61   0.544    -98.91613    186.0331
_cons |   10254.95   2430.641     4.22   0.000      5408.39    15101.51
------------------------------------------------------------------------------

. matrix list e(V)

symmetric e(V)[3,3]
mpg       trunk       _cons
mpg   5249.5646
trunk   3569.5316   5105.6316
_cons  -169049.76  -147284.49   5908013.8

. matrix list V2

symmetric V2[3,3]
mpg       trunk       _cons
mpg   5249.5646
trunk   3569.5316   5105.6316
_cons  -169049.76  -147284.49   5908013.8


Cluster-robust standard errors

The cluster-robust estimator of the VCE is frequently used when the data have a panel structure, also known as a longitudinal structure. This VCE accounts for the within-group correlation of the errors, and it is given by

$\widehat{V}_{cluster}=\frac{N}{N-k}\frac{g}{g-1} \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1} \Mb_c \left( \sum_{i=1}^N \xb_i’\xb_i \right)^{-1}$

where

$\Mb_c=\sum_{j=1}^g \Xb_j’ (\widehat{\eb}_j \widehat{\eb}_j’) \Xb_j$

$$\Xb_j$$ is the $$n_j\times k$$ matrix of observations on $$\xb_i$$ in group $$j$$, $$\widehat{\eb}_j$$ is the $$n_j\times 1$$ vector of residuals in group $$j$$, and $$g$$ is the number of groups. See Cameron and Trivedi (2005), Wooldridge (2010), and [R] regress for derivations and discussions.

matrix opaccum computes the formula for $$\Mb_c$$. Below, I create the group variable cvar from rep78 and use matrix opaccum to compute $$\Mb_c$$ and $$\widehat{V}_{cluster}$$

Example 9: A cluster-robust VCE

. generate cvar = cond( missing(rep78), 6, rep78)

. tab cvar

cvar |      Freq.     Percent        Cum.
------------+-----------------------------------
1 |          2        2.70        2.70
2 |          8       10.81       13.51
3 |         30       40.54       54.05
4 |         18       24.32       78.38
5 |         11       14.86       93.24
6 |          5        6.76      100.00
------------+-----------------------------------
Total |         74      100.00

. local Nc = r(r)

. sort cvar

. matrix opaccum M2     = mpg trunk , group(cvar) opvar(res)

. matrix V2          = ((N'-1)/(N'-3))*(Nc'/(Nc'-1))*xpxi*M2*xpxi


I now verify that my computations match those reported by regress.

Example 10: Comparing computations of cluster-robust VCE

. regress price mpg trunk, vce(cluster cvar)

Linear regression                               Number of obs     =         74
F(2, 5)           =       9.54
Prob > F          =     0.0196
R-squared         =     0.2222
Root MSE          =     2637.6

(Std. Err. adjusted for 6 clusters in cvar)
------------------------------------------------------------------------------
|               Robust
price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |  -220.1649   93.28127    -2.36   0.065     -459.952    19.62226
trunk |   43.55851   58.89644     0.74   0.493    -107.8396    194.9566
_cons |   10254.95   2448.547     4.19   0.009     3960.758    16549.14
------------------------------------------------------------------------------

. matrix list e(V)

symmetric e(V)[3,3]
mpg       trunk       _cons
mpg   8701.3957
trunk   4053.5381   3468.7911
_cons     -223021  -124190.97   5995384.3

. matrix list V2

symmetric V2[3,3]
mpg       trunk       _cons
mpg   8701.3957
trunk   4053.5381   3468.7911
_cons     -223021  -124190.97   5995384.3


Done and undone

I reviewed the formulas that underlie the OLS estimator and showed how to compute them using Stata matrix commands and functions. In the next two posts, I write an ado-command that implements these formulas.

References

Cameron, A. C., and P. K. Trivedi. 2005. Microeconometrics: Methods and applications. Cambridge: Cambridge University Press.

Stock, J. H., and M. W. Watson. 2010. Introduction to Econometrics. 3rd ed. Boston, MA: Addison Wesley New York.

Wooldridge, J. M. 2010. Econometric Analysis of Cross Section and Panel Data. 2nd ed. Cambridge, Massachusetts: MIT Press.

Wooldridge, J. M. 2015. Introductory Econometrics: A Modern Approach. 6th ed. Cincinnati, Ohio: South-Western.

Categories: Programming Tags:

## Programming an estimation command in Stata: A first ado command

I discuss the code for a simple estimation command to focus on the details of how to implement an estimation command. The command that I discuss estimates the mean by the sample average. I begin by reviewing the formulas and a do-file that implements them. I subsequently introduce ado-file programming and discuss two versions of the command. Along the way, I illustrate some of the postestimation features that work after the command.

The formulas for our estimator

The formulas for the sample average and its estimated sampling variance, assuming an independently and identically distributed process, are

$\widehat{\mu} = 1/N \sum_{i=1}^N y_i$

$\widehat{Var}(\widehat{\mu}) = 1/[N(N-1)] \sum_{i=1}^N (y_i-\widehat{\mu})^2$

The code mean1.do performs these computations on price from the auto dataset.

Code block 1: mean1.do

// version 1.0.0 20Oct2015
version 14
sysuse auto
quietly summarize price
local sum          = r(sum)
local N            = r(N)
local mu           = (1/N')*sum'
generate double e2 = (price - mu')^2
quietly summarize e2
local V            = (1/((N')*(N'-1)))*r(sum)
display "muhat = " mu'
display "   V  = " V'


mean1.do uses summarize to compute the summations. Lines 5–7 and line 11 store results stored by summarize in r() into local macros that are subsequently used to compute the formulas. I recommend that you use double, instead of the default float, to compute all variables used in your formulas because it is almost always worth taking up the extra memory to gain the extra precision offered by double over float. (Essentially, each variable takes up twice as much space, but you get calculations that are correct to about $$10^{-16}$$ instead of $$10^{-8}$$.)

These calculations yield

Example 1: Computing the average and its sampling variance

. do mean1

. // version 1.0.0 20Oct2015
. version 14

. sysuse auto
(1978 Automobile Data)

. quietly summarize price

. local sum          = r(sum)

. local N            = r(N)

. local mu           = (1/N')*sum'

. generate double e2 = (price - mu')^2

. quietly summarize e2

. local V            = (1/((N')*(N'-1)))*r(sum)

. display "muhat = " mu'
muhat = 6165.2568

. display "   V  = " V'
V  = 117561.16

.
end of do-file


Now I verify that mean produces the same results.

Example 2: Results from mean

. mean price

Mean estimation                   Number of obs   =         74

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
price |   6165.257   342.8719      5481.914      6848.6
--------------------------------------------------------------

. matrix list e(b)

symmetric e(b)[1,1]
price
y1  6165.2568

. matrix list e(V)

symmetric e(V)[1,1]
price
price  117561.16


The code in mymean1.ado performs the same calculations as mean1.do. (The file mymean1.ado is in my current working directory.)

*! version 1.0.0 20Oct2015
program define mymean1
version 14

quietly summarize price
local sum          = r(sum)
local N            = r(N)
local mu           = (1/N')*sum'
capture drop e2				// Drop e2 if it exists
generate double e2 = (price - mu')^2
quietly summarize e2
local V            = (1/((N')*(N'-1)))*r(sum)
display "muhat = " mu'
display "   V  = " V'
end


Line 1 of mymean1.ado specifies that file defines the command mymean1. The command name must be the same as the file name that precedes the suffix .ado. The mymean1 command performs the same computations as the do-file mean1.do.

Example 3: Results from mymean1

. mymean1
muhat = 6165.2568
V  = 117561.16


A slightly better command

We want our command to be reusable; we want it to estimate the mean for any variable in memory, instead of only for price as performed by mymean1.ado. On line 5 of mymean2.ado, we use the syntax command to store the name of the variable specified by the user into the local macro varlist, which we use in the remainder of the computations.

*! version 2.0.0 20Oct2015
program define mymean2
version 14

syntax varlist
display "The local macro varlist contains varlist'"

quietly summarize varlist'
local sum          = r(sum)
local N            = r(N)
local mu           = (1/N')*sum'
capture drop e2				// Drop e2 if it exists
generate double e2 = (varlist' - mu')^2
quietly summarize e2
local V            = (1/((N')*(N'-1)))*r(sum)
display  "The average of varlist' is " mu'
display  "The estimated variance of the average is " V'
end


The extremely powerful syntax command puts the elements of Stata syntax specified by the user into local macros and throws errors when the user makes a mistake. I will discuss syntax in greater detail in subsequent posts.

I begin by illustrating how to replicate the previous results.

Example 4: Results from mymean2 price

. mymean2 price
The local macro varlist contains price
The average of price is 6165.2568
The estimated variance of the average is 117561.16


I now illustrate that it works for another variable.

Example 5: Results from mymean2 trunk

. mymean2 trunk
The local macro varlist contains trunk
The average of trunk is 13.756757
The estimated variance of the average is .24724576

. mean trunk

Mean estimation                   Number of obs   =         74

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
trunk |   13.75676   .4972381      12.76576    14.74775
--------------------------------------------------------------

. display "The variance of the estimator is " (_se[trunk])^2
The variance of the estimator is .24724576


Storing results in e()

mymean2.ado does not save the results that it displays. We fix this problem in mymean3.ado. Line 2 specifies the option e-class on program define to make mymean3 an e-class command. Line 18 uses ereturn post to move the matrix of point estimates (b) and the estimated variance-covariance of the estimator (VCE) into e(b) and e(V). The estimation-postestimation framework uses parameter names for display, hypothesis tests, and other features. In lines 15 and 16, we put those names into the column stripes of the vector of estimates and the estimated VCE. In line 17, we put those names into the row stripe of the estimated VCE.

*! version 3.0.0 20Oct2015
program define mymean3, eclass
version 14

syntax varlist

quietly summarize varlist'
local sum          = r(sum)
local N            = r(N)
matrix b           = (1/N')*sum'
capture drop e2				// Drop e2 if it exists
generate double e2 = (varlist' - b[1,1])^2
quietly summarize e2
matrix V           = (1/((N')*(N'-1)))*r(sum)
matrix colnames b  = varlist'
matrix colnames V  = varlist'
matrix rownames V  = varlist'
ereturn post b V
ereturn display
end


The ereturn display command on line 19 of mymean3.ado easily creates a standard output table using the results now stored in e(b) and e(V).

Example 6: Results from mymean3 trunk

. mymean3 trunk
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trunk |   13.75676   .4972381    27.67   0.000     12.78219    14.73133
------------------------------------------------------------------------------


Estimation-postestimation

test, lincom, testnl, nlcom, and other Wald-based estimation-postestimation features work after mymean3 because all the required information is stored in e(b) and e(V).

To illustrate, I perform a Wald of the null hypothesis that the mean of trunk is $$11$$.

Example 7: test works after mymean3

. test _b[trunk]==11

( 1)  trunk = 11

chi2(  1) =   30.74
Prob > chi2 =    0.0000


The results stored in e() are the glue that holds the estimation-postestimation framework together. We have only stored e(b) and e(V) so far, so not all the standard features are working yet. (But we will get there in the #StataProgramming series.)

Using temporary names for global objects

Stata variables and matrices are global, as discussed in my previous blog post. We need some safe names for global objects. These safe names should not be in use elsewhere, and they should be temporary in that we want Stata to drop the corresponding objects when the command finishes. The tempvar and tempname commands put safe names into local macros and then drop the corresponding objects when the ado-file or do-file finishes. We explicitly dropped e2, if it existed, in line 9 of code block 2, in line 12 of code block 3, and in line 11 of code block 4. We do not need such a line in code block, because we are using temporary variable names.

In line 7 of mymean4.ado, the tempvar command puts a safe name into the local macro e2. In line 8 of mymean4.ado, the tempname command puts safe names into the local macros b and V. I illustrate the format followed by these safe names by displaying them on lines 9–11. The output reveals that a leading pair of underscores is followed by numbers and capital letters. Line 15 illustrates the use of these safe names. Instead of creating the matrix b, we create the matrix whose name is stored in the local macro b. In line 8, the tempname command created the local macro b to hold a safe name.

*! version 4.0.0 20Oct2015
program define mymean4, eclass
version 14

syntax varlist

tempvar e2
tempname b V
display "The safe name in e2 is e2'"
display "The safe name in b is b'"
display "The safe name in V is V'"
quietly summarize varlist'
local sum            = r(sum)
local N              = r(N)
matrix b'           = (1/N')*sum'
generate double e2' = (varlist' - b'[1,1])^2
quietly summarize e2'
matrix V'           = (1/((N')*(N'-1)))*r(sum)
matrix colnames b'  = varlist'
matrix colnames V'  = varlist'
matrix rownames V'  = varlist'
ereturn post b' V'
ereturn display
end


This code produces the output

Example 8: Results from mymean4 trunk

. mymean4 trunk
The safe name in e2 is __000000
The safe name in b is __000001
The safe name in V is __000002
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trunk |   13.75676   .4972381    27.67   0.000     12.78219    14.73133
------------------------------------------------------------------------------


Removing the lines that display the safe names contained in the local macros yields mymean5.ado.

*! version 5.0.0 20Oct2015
program define mymean5, eclass
version 14

syntax varlist

tempvar e2
tempname b V
quietly summarize varlist'
local sum            = r(sum)
local N              = r(N)
matrix b'           = (1/N')*sum'
generate double e2' = (varlist' - b'[1,1])^2
quietly summarize e2'
matrix V'           = (1/((N')*(N'-1)))*r(sum)
matrix colnames b'  = varlist'
matrix colnames V'  = varlist'
matrix rownames V'  = varlist'
ereturn post b' V'
ereturn display
end


This code produces the output

Example 9: Results from mymean5 trunk

. mymean5 trunk
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trunk |   13.75676   .4972381    27.67   0.000     12.78219    14.73133
------------------------------------------------------------------------------


Done and undone

I illustrated some basic ado-file programming techniques by implementing a command that estimates the mean of variable. Even though we have a command that produces correct, easy-to-read output that has some estimation-postestimation features, we have only scratched the surface of what we usually want to do in an estimation command. I dig a little deeper in the next few posts by developing a command that performs ordinary least-squares estimation.

Categories: Programming Tags:

## Programming an estimation command in Stata: Global macros versus local macros

I discuss a pair of examples that illustrate the differences between global macros and local macros. You can view this post as a technical appendix to the previous post in the #StataProgramming series, which introduced global macros and local macros.

In every command I write, I use local macros to store stuff in a workspace that will not alter a user’s data and to make my code easier to read. A good understanding of the differences between global macros and local macros helps me to write better code. The essential differences between global macros and local macros can be summarized in two points.

1. There is only one global macro with a specific name in Stata, and its contents can be accessed or changed by a Stata command executed at any Stata level.
2. In contrast, each Stata level can have a local macro of a specific name, and each one’s contents cannot be accessed or changed by commands executed at other Stata levels.

If you are already comfortable with 1 and 2, skip the remainder of this post.

Global macros are global

The do-files globala.do and globalb.do in code blocks globala and globalb illustrate what it means to be global.

Code block 1: globala.do

*-------------------------------Begin globala.do ---------------
*! globala.do
*  In this do-file we define the global macro vlist, but we
*  do not use it
global vlist var1 var2 var3

do globalb
*-------------------------------End globala.do ---------------


Code block 2: globalb.do

*-------------------------------Begin globalb.do ---------------
*! globalb.do
*  In this do-file, we use the global macro vlist, defined in globala.do

display "The global macro vlist contains $vlist" *-------------------------------End globalb.do ---------------  The easiest way to see what this code does is to execute it; the output is in example 1. Example 1: Output from do globala . do globala . *-------------------------------Begin globala.do --------------- . *! globala.do . * In this do-file we define the global macro vlist, but we . * do not use it . global vlist var1 var2 var3 . . do globalb . *-------------------------------Begin globalb.do --------------- . *! globalb.do . * In this do-file, we use the global macro vlist, defined in globala.do . . display "The global macro vlist contains$vlist"
The global macro vlist contains var1 var2 var3

. *-------------------------------End globalb.do ---------------
.
end of do-file

. *-------------------------------End globala.do ---------------
.
end of do-file


Line 5 of globalb.do can access the contents of vlist created on line 5 of globala.do because vlist is a global macro.

Figure 1 makes this same point graphically: the global macro vlist is in global memory, and a command executed anywhere can access or change the contents of vlist.

Figure 1: A global macro in global memory

Local macros are local

The do-files locala.do and localb.do in code blocks 3 and 4 illustrate what it means to be local.

Code block 3: locala.do

*-------------------------------Begin locala.do ---------------
*! locala.do
local mylist "a b c"
display "mylist contains mylist'"

do localb

display "mylist contains mylist'"
*-------------------------------End locala.do ---------------


Code block 4: localb.do

*-------------------------------Begin localb.do ---------------
*! localb.do
local mylist "x y z"
display "mylist contains mylist'"
*-------------------------------End localb.do ---------------


The easiest way to see what this code does is to execute it; the output is in example 2.

Example 2: Output from do locala

. do locala

. *-------------------------------Begin locala.do ---------------
. *! locala.do
. local mylist "a b c"

. display "mylist contains mylist'"
mylist contains a b c

.
. do localb

. *-------------------------------Begin localb.do ---------------
. *! localb.do
. local mylist "x y z"

. display "mylist contains mylist'"
mylist contains x y z

. *-------------------------------End localb.do ---------------
.
end of do-file

.
. display "mylist contains mylist'"
mylist contains a b c

. *-------------------------------End locala.do ---------------
.
end of do-file


The code in blocks 3 and 4 and the output in example 2 illustrate that a command executed at the level of localb.do cannot change the local macro mylist that is local to locala.do. Line 8 of locala.do displays the contents of the mylist local to locala.do. The contents are still a b c after localb.do finishes because the local macro mylist created on line 3 of locala.do is local to locala.do and it is unaffected by the mylist created on line 3 of localb.do.

Figure 2 makes this point graphically. The contents of the local macro mylist that is local to locala.do can be accessed and changed by commands run in locala.do, but not by commands run in localb.do. Analogously, the contents of the local macro mylist that is local to localb.do can be accessed and changed by commands run in localb.do, but not by commands run in locala.do.

Figure 2: Local macros are local to do-files

Done and Undone

I essentially provided you with a technical appendix to the previous #StataProgramming post. I illustrated that global macros are global and that local macros are local. I use the concepts developed thus far to present an ado-command in the next post.

Categories: Programming Tags:

## Programing an estimation command in Stata: Where to store your stuff

If you tell me “I program in Stata”, it makes me happy, but I do not know what you mean. Do you write scripts to make your research reproducible, or do you write Stata commands that anyone can use and reuse? In the series #StataProgramming, I will show you how to write your own commands, but I start at the beginning. Discussing the difference between scripts and commands here introduces some essential programming concepts and constructions that I use to write scripts and commands.

Scripts versus commands

A script is a program that always performs the same tasks on the same inputs and produces exactly the same results. Scripts in Stata are known as do-files and the files containing them end in .do. For example, I could write a do-file to

1. read in the National Longitudinal Study of Youth (NLSY) dataset,
2. clean the data,
3. form a sample for some population, and
4. run a bunch of regressions on the sample.

This structure is at the heart of reproducible research; produce the same results from the same inputs every time. Do-files have a one-of structure. For example, I could not somehow tell this do-file that I want it to perform the analogous tasks on the Panel Study on Income Dynamics (PSID). Commands are reusable programs that take arguments to perform a task on any data of certain type. For example, regress performs ordinary least squares on the specified variables regardless of whether they come from the NLSY, PSID, or any other dataset. Stata commands are written in the automatic do-file (ado) language; the files containing them end in .ado. Stata commands written in the ado language are known as ado-commands.

An example do-file

The commands in code block 1 are contained in the file doex.do in the current working directory of my computer.

Code block 1: doex.do

// version 1.0.0  04Oct2015 (This line is comment)
version 14                     // version #.# fixes the version of Stata
use http://www.stata.com/data/accident2.dta
summarize accidents tickets


We execute the commands by typing do doex which produces

Example 1: Output from do doex

. do doex

. // version 1.0.0  04Oct2015 (This line is comment)
. version 14                     // version #.# fixes the version of Stata

. use http://www.stata.com/data/accident2.dta

. summarize accidents tickets

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
accidents |        948    .8512658    2.851856          0         20
tickets |        948    1.436709    1.849456          0          7

.
.
end of do-file

1. Line 1 in doex.do is a comment that helps to document the code but is not executed by Stata. The // initiates a comment. Anything following the // on that line is ignored by Stata.
2. In the comment on line 1, I put a version number and the date that I last changed this file. The date and the version help me keep track of the changes that I make as I work on the project. This information also helps me answer questions from others with whom I have shared a version of this file.
3. Line 2 specifies the definition of the Stata language that I use. Stata changes over time. Setting the version ensures that the do-file continues to run and that the results do not change as the Stata language evolves.
4. Line 3 reads in the accident.dta dataset.
5. Line 4 summarizes the variables accidents and tickets.

Storing stuff in Stata

Programming in Stata is like putting stuff into boxes, making Stata change the stuff in the boxes, and getting the changed stuff out of the boxes. For example, code block 2 contains the code for doex2.do, whose output I display in example 2

Code block 2: doex2.do

// version 1.0.0  04Oct2015 (This line is comment)
version 14                     // version #.# fixes the version of Stata
use http://www.stata.com/data/accident2.dta
generate ln_traffic = ln(traffic)
summarize ln_traffic


Example 2: Output from do doex2

. do doex2

. // version 1.0.0  04Oct2015 (This line is comment)
. version 14                     // version #.# fixes the version of Stata

. use http://www.stata.com/data/accident2.dta

. generate ln_traffic = ln(traffic)

. summarize ln_traffic

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
ln_traffic |        948    1.346907    1.004952  -5.261297   2.302408

.
.
end of do-file


In line 4 of code block 2, I generate the new variable ln_traffic which I summarize on line 5. doex2.do uses generate to change what is in the box ln_traffic and uses summarize to get a function of the changed stuff out of the box. Stata variables are the most frequently used box type in Stata, but when you are programming, you will also rely on Stata matrices.

There can only be one variable named traffic in a Stata dataset and its contents can be viewed or changed interactively, by a do-file, or by an ado-file command. Similarly, there can only be one Stata matrix named beta in a Stata session and its contents can be viewed or changed interactively, by a do-file, or by an ado-file command. Stata variables and Stata matrices are global boxes because there can only be one Stata variable or Stata matrix in a Stata session and its contents can be viewed or changed anywhere in a Stata session.

The opposite of global is local. If it is local in Stata, its contents can only be accessed or changed in the interactive session, in a particular do-file, or a in particular ado-file.

Although I am discussing do-files at the moment, remember that we are learning techniques to write commands. It is essential to understand the differences between global boxes and local boxes to program commands in Stata. Global boxes, like variables, could contain data that the users of your command do not want changed. For example, a command you write should never change a user’s variable in a way that was not requested.

Levels of Stata

The notion that there are levels of Stata can help explain the difference between global boxes and local boxes. Suppose that I run 2 do-files or ado-files. Think of the interactive Stata session as level 0 of Stata, and think of each do-file or ado-file as being Stata levels 1 and 2. Global boxes like variables and matrices live in global memory that can be accessed or changed from a Stata command executed in level 0, 1, or 2. Local boxes can only be accessed or changed by a Stata command within a particular level of Stata. (This description is not exactly how Stata works, but the details about how Stata really handles levels are not important here.)

Figure 1 depicts this structure.

Memory by Stata level

Figure 1 clarifies

• that commands executed at all Stata levels can access and change the objects in global memory,
• that only commands executed at Stata level 0 can access and change the objects local to Stata level 0,
• that only commands executed at Stata level 1 can access and change the objects local to Stata level 1, and
• that only commands executed at Stata level 2 can access and change the objects local to Stata level 2.

Global and local macros: Storing and extracting

Macros are Stata boxes that hold information as characters, also known as strings. Stata has both global macros and local macros. Global macros are global and local macros are local. Global macros can be accessed and changed by a command executed at any Stata level. Local macros can be accessed and changed only by a command executed at a specific Stata level.

The easiest way to begin to understand global macros is to put something into a global macro and then to get it back out. Code block 3 contains the code for global1.do which stores and the retrieves information from a global macro.

Code block 3: global1.do

// version 1.0.0  04Oct2015
version 14
global vlist "y x1 x2"
display "vlist contains $vlist"  Example 3: Output from do global1 . do global1 . // version 1.0.0 04Oct2015 . version 14 . global vlist "y x1 x2" . display "vlist contains$vlist"
vlist contains y x1 x2

.
end of do-file


Line 3 of code block 3 puts the string y x1 x2 into the global macro named vlist. To extract what I put into a global macro, I prefix the name of global macro with a $. Line 4 of the code block and its output in example 3 illustrate this usage by extracting and displaying the contents of vlist. Code block 4 contains the code for local1.do and its output is given in example 4. They illustrate how to put something into a local macro and how to extract something from it. Code block 4: local1.do // version 1.0.0 04Oct2015 version 14 local vlist "y x1 x2" display "vlist contains vlist'"  Example 4: Output from do global1 . do local1 . // version 1.0.0 04Oct2015 . version 14 . local vlist "y x1 x2" . display "vlist contains vlist'" vlist contains y x1 x2 . end of do-file  Line 3 of code block 3 puts the string y x1 x2 into the local macro named vlist. To extract what I put into a local macro I enclose the name of the local macro between a single left quote (‘) and a single right quote (’). Line 4 of code block 3 displays what is contained in the local macro vlist and its output in example 4 illustrates this usage. Getting stuff from Stata commands Now that we have boxes, I will show you how to store stuff computed by Stata in these boxes. Analysis commands, like summarize, store their results in r(). Estimation commands, like regress, store their results in e(). Somewhat tautologically, commands that store their results in r() are also known as r-class commands and commands that store their results in e() are also known as e-class commands. I can use return list to see results stored by an r-class command. Below, I list out what summarize has stored in r() and compute the mean from the stored results. Example 5: Getting results from an r-class command . use http://www.stata.com/data/accident2.dta, clear . summarize accidents Variable | Obs Mean Std. Dev. Min Max -------------+--------------------------------------------------------- accidents | 948 .8512658 2.851856 0 20 . return list scalars: r(N) = 948 r(sum_w) = 948 r(mean) = .8512658227848101 r(Var) = 8.133081817331211 r(sd) = 2.851855854935732 r(min) = 0 r(max) = 20 r(sum) = 807 . local sum = r(sum) . local N = r(N) . display "The mean is " sum'/N' The mean is .85126582  Estimation commands are more formal than analysis commands, so they save more stuff. Official Stata estimation commands save lots of stuff, because they follow lots of rules that make postestimation easy for users. Do not be alarmed by the number of things stored by poisson. Below, I list out the results stored by poisson and create a Stata matrix that contains the coefficient estimates. Example 6: Getting results from an e-class command . poisson accidents traffic tickets male Iteration 0: log likelihood = -377.98594 Iteration 1: log likelihood = -370.68001 Iteration 2: log likelihood = -370.66527 Iteration 3: log likelihood = -370.66527 Poisson regression Number of obs = 948 LR chi2(3) = 3357.64 Prob > chi2 = 0.0000 Log likelihood = -370.66527 Pseudo R2 = 0.8191 ------------------------------------------------------------------------------ accidents | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- traffic | .0764399 .0129856 5.89 0.000 .0509887 .1018912 tickets | 1.366614 .0380641 35.90 0.000 1.29201 1.441218 male | 3.228004 .1145458 28.18 0.000 3.003499 3.45251 _cons | -7.434478 .2590086 -28.70 0.000 -7.942126 -6.92683 ------------------------------------------------------------------------------ . ereturn list scalars: e(rank) = 4 e(N) = 948 e(ic) = 3 e(k) = 4 e(k_eq) = 1 e(k_dv) = 1 e(converged) = 1 e(rc) = 0 e(ll) = -370.6652697757637 e(k_eq_model) = 1 e(ll_0) = -2049.485325326086 e(df_m) = 3 e(chi2) = 3357.640111100644 e(p) = 0 e(r2_p) = .8191422669899876 macros: e(cmdline) : "poisson accidents traffic tickets male" e(cmd) : "poisson" e(predict) : "poisso_p" e(estat_cmd) : "poisson_estat" e(chi2type) : "LR" e(opt) : "moptimize" e(vce) : "oim" e(title) : "Poisson regression" e(user) : "poiss_lf" e(ml_method) : "e2" e(technique) : "nr" e(which) : "max" e(depvar) : "accidents" e(properties) : "b V" matrices: e(b) : 1 x 4 e(V) : 4 x 4 e(ilog) : 1 x 20 e(gradient) : 1 x 4 functions: e(sample) . matrix b = e(b) . matrix list b b[1,4] accidents: accidents: accidents: accidents: traffic tickets male _cons y1 .07643992 1.366614 3.2280044 -7.434478  Done and Undone In this second post in the series #StataProgramming, I discussed the difference between scripts and commands, I provided an introduction to the concepts of global and local memory objects, I discussed global macros and local macros, and I showed how to access results stored by other commands. In the next post in the series #StataProgramming, I discuss an example that further illustrates the differences between global macros and local macros. Categories: Programming Tags: ## Programming estimators in Stata: Why you should Distributing a Stata command that implements a statistical method will get that method used by lots of people. They will thank you. And, they will cite you! This post is the first in the series #StataProgramming about programing an estimation command in Stata that uses Mata to do the numerical work. In the process of showing you how to program an estimation command in Stata, I will discuss do-file programming, ado-file programming, and Mata programming. When the series ends, you will be able to write Stata commands. Stata users like its predictable syntax and its estimation-postestimation structure that facilitates hypothesis testing, specification tests, and parameter interpretation. To help you write Stata commands that people want to use, I illustrate how Stata syntax is predictable and give an overview of the estimation-postestimation structure that you will want to emulate in your programs. Stata structure by example I use and describe some simulated data about the number of traffic accidents observed on 948 people. Example 1: Accident data . use http://www.stata.com/data/accident2.dta . describe Contains data from http://www.stata.com/data/accident2.dta obs: 948 vars: 6 23 Sep 2015 13:04 size: 22,752 -------------------------------------------------------------------------------- storage display value variable name type format label variable label -------------------------------------------------------------------------------- kids float %9.0g number of children cvalue float %9.0g car value index tickets float %9.0g number of tickets in last 2 years traffic float %9.0g local traffic index, larger=>worse male float %9.0g 1=>man, 0=>woman accidents float %9.0g number of traffic in last 5 years -------------------------------------------------------------------------------- Sorted by:  Stata’s predictable syntax I estimate the parameters of a Poisson regression model for accidents as a function of traffic conditions (traffic), an indicator for being a male driver (male), and the number of tickets received in the last two years (tickets). Example 2: A Poisson model for accidents . poisson accidents traffic male tickets , vce(robust) Iteration 0: log pseudolikelihood = -377.98594 Iteration 1: log pseudolikelihood = -370.68001 Iteration 2: log pseudolikelihood = -370.66527 Iteration 3: log pseudolikelihood = -370.66527 Poisson regression Number of obs = 948 Wald chi2(3) = 1798.65 Prob > chi2 = 0.0000 Log pseudolikelihood = -370.66527 Pseudo R2 = 0.8191 ------------------------------------------------------------------------------ | Robust accidents | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- traffic | .0764399 .0165119 4.63 0.000 .0440772 .1088027 male | 3.228004 .1232081 26.20 0.000 2.986521 3.469488 tickets | 1.366614 .0328218 41.64 0.000 1.302284 1.430943 _cons | -7.434478 .2413188 -30.81 0.000 -7.907454 -6.961502 ------------------------------------------------------------------------------  I want to focus on the structure in this example so that you can use it to make your commands easier to use. In particular, I want to discuss the structure of the command syntax and to point out that the output is easy to read and interpret because it is a standard Stata output table. For estimators that table almost always reports estimates (often coefficients), standard errors, tests against zero and their$p\$-values, and confidence intervals.

Stata syntax is predictable, which makes it easy to use. Stata users “speak Stata” and do not even notice the details. I highlight some of these details so that we can make the syntax of the commands we write predictable. Here are some of the standard syntax elements illustrated in example 2.

1. The command has four syntactical elements;
1. command name (poisson),
2. list of variable names (accidents traffic male tickets),
3. a comma,
4. an option (vce(robust)).
2. In the list of variable names, the name of the dependent variable is first and it is followed by the names of the independent variables.
3. The job of the comma is to separate the command name and variable list from the option or options.

The output is also structured; it is composed of an iteration log, a header, and a standard output table.

Estimation-postestimation framework

As a Stata user, I could now use the estimation-postestimation framework. For example, I could perform a Wald test of the hypothesis that the coefficient on male is 3.

Example 3: A Wald test of a linear restriction

. test male = 3

( 1)  [accidents]male = 3

chi2(  1) =    3.42
Prob > chi2 =    0.0642


or I could perform a Wald test of the nonlinear hypothesis that the ratio of the coefficient on male to the ratio of the coefficient on tickets is 2.

Example 4: A Wald test of a nonlinear restriction

. testnl _b[male]/_b[tickets] = 2

(1)  _b[male]/_b[tickets] = 2

chi2(1) =       19.65
Prob > chi2 =        0.0000


I could also predict the mean of accidents for each observation and summarize the results.

Example 5: Summarizing the predicted conditional means

. predict nhat
(option n assumed; predicted number of events)

. summarize nhat

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
nhat |        948    .8512658    2.971087   .0006086    29.0763


Finally, I could use margins to estimate conditional or population-averaged parameters that are functions of the parameters in the original model. I use margins to estimate the average number of accidents that would be observed if each individual received 0 tickets, or 1 ticket, or 2 tickets, …, or 7 tickets. See [R] margins, Long and Freese (2006, sec. 4.4.2-4.4.3), and Cameron and Trivedi (2010, 10.5.6{10.6.9) for introductions to estimating functions of the the model parameters by margins.

Example 6: Estimating functions of model parameters

. margins, at(tickets=(0 1 2 3 4 5 6 7))

Predictive margins                              Number of obs     =        948
Model VCE    : Robust

Expression   : Predicted number of events, predict()

1._at        : tickets         =           0

2._at        : tickets         =           1

3._at        : tickets         =           2

4._at        : tickets         =           3

5._at        : tickets         =           4

6._at        : tickets         =           5

7._at        : tickets         =           6

8._at        : tickets         =           7

------------------------------------------------------------------------------
|            Delta-method
|     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1  |   .0097252   .0015387     6.32   0.000     .0067094     .012741
2  |   .0381426   .0048762     7.82   0.000     .0285854    .0476998
3  |   .1495971   .0148157    10.10   0.000      .120559    .1786353
4  |   .5867272   .0432256    13.57   0.000     .5020066    .6714478
5  |   2.301172   .1302033    17.67   0.000     2.045978    2.556366
6  |   9.025308   .5049176    17.87   0.000     8.035688    10.01493
7  |   35.39769   2.555679    13.85   0.000     30.38865    40.40673
8  |   138.8315   13.49606    10.29   0.000     112.3797    165.2832
------------------------------------------------------------------------------


The glue

The estimation results stored in e() are the glue that holds together the estimation-postestimation framework. The poisson command stores lots of stuff in e(). I could use ereturn list to list all this stuff, but there are many stored objects that do not interest you yet.

Most of the estimation-postestimation features that I discussed were implemented using e(b), e(V), and e(predict), which are the vector of point estimates, the estimated VCE, and the name of the command that implements predict after poisson.

I will show how to store what you need in e() in the #StataProgramming series.

Structure of Stata commands

Here is an outline of the tasks performed by a Stata estimation command.

1. Parse the input to the command.
2. Compute results.
3. Store results in e()
4. Display output.

You need to write a predict command to complete the estimation-postestimation framework. After you have stored the estimation results and written the predict command, margins works.

I will explain each of these steps in the #StataProgramming series of posts.

Use this structure to your advantage. To make your command easy to use, design it to have the predictable syntax implemented in other commands and make it work in the estimation-postestimation framework. This task is far easier than it sounds. In fact, it is just plain easy. The Stata language steers you in this direction.

Done and undone

I will teach you how to program an estimation command in Stata in the #StataProgramming series. I will also show you how do the numerical work in Mata. I discussed the following points, in this first post.

1. The predictable structure of Stata syntax makes Stata easy to use. You should emulate this structure, so that your commands are easy to use.
2. The estimation-postestimation framework makes inference and advanced estimation simple. It is easy for you to make your command work with this framework.
3. The estimation results stored in e(), and the predict command, are the glue that holds the estimation-postestimation framework together.

In the next post, I discuss do-file programming tools that I will subsequently use to parse the input to the command.

References

Cameron, A. C., and P. K. Trivedi. 2010. Microeconometrics Using Stata. Revised ed. College Station, Texas: Stata Press.

Long, J. S., and J. Freese. 2014. Regression models for categorical dependent variables using Stata. 3rd ed. College Station, Texas: Stata Press.

Categories: Programming Tags:

## Monte Carlo simulations using Stata

Overview

A Monte Carlo simulation (MCS) of an estimator approximates the sampling distribution of an estimator by simulation methods for a particular data-generating process (DGP) and sample size. I use an MCS to learn how well estimation techniques perform for specific DGPs. In this post, I show how to perform an MCS study of an estimator in Stata and how to interpret the results.

Large-sample theory tells us that the sample average is a good estimator for the mean when the true DGP is a random sample from a $$\chi^2$$ distribution with 1 degree of freedom, denoted by $$\chi^2(1)$$. But a friend of mine claims this estimator will not work well for this DGP because the $$\chi^2(1)$$ distribution will produce outliers. In this post, I use an MCS to see if the large-sample theory works well for this DGP in a sample of 500 observations.

A first pass at an MCS

I begin by showing how to draw a random sample of size 500 from a $$\chi^2(1)$$ distribution and how to estimate the mean and a standard error for the mean.

Example 1: The mean of simulated data

. drop _all
. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345
. generate y = rchi2(1)
. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------


I specified set seed 12345 to set the seed of the random-number generator so that the results will be reproducible. The sample average estimate of the mean from this random sample is $$0.91$$, and the estimated standard error is $$0.055$$.

If I had many estimates, each from an independently drawn random sample, I could estimate the mean and the standard deviation of the sampling distribution of the estimator. To obtain many estimates, I need to repeat the following process many times:

1. Draw from the DGP
2. Compute the estimate
3. Store the estimate.

I need to know how to store the many estimates to proceed with this process. I also need to know how to repeat the process many times and how to access Stata estimates, but I put these details into appendices I and II, respectively, because many readers are already familiar with these topics and I want to focus on how to store the results from many draws.

I want to put the many estimates someplace where they will become part of a dataset that I can subsequently analyze. I use the commands postfile, post, and postclose to store the estimates in memory and write all the stored estimates out to a dataset when I am done. Example 2 illustrates the process, when there are three draws.

Example 2: Estimated means of three draws

. set seed 12345

. postfile buffer mhat using mcs, replace

. forvalues i=1/3 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y])
7. }

. postclose buffer

. use mcs, clear

. list

+----------+
|     mhat |
|----------|
1. | .9107645 |
2. |  1.03821 |
3. | 1.039254 |
+----------+


The command

postfile buffer mhat using mcs, replace


creates a place in memory called buffer in which I can store the results that will eventually be written out to a dataset. mhat is the name of the variable that will hold the estimates in the new dataset called mcs.dta. The keyword using separates the new variable name from the name of the new dataset. I specified the option replace to replace any previous versions of msc.dta with the one created here.

I used

forvalues i=1/3 {


to repeat the process three times. (See appendix I if you want a refresher on this syntax.) The commands

quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y


drop the previous data, draw a sample of size 500 from a $$\chi^2(1)$$ distribution, and estimate the mean. (The quietly before each command suppresses the output.) The command

post buffer (_b[y])


stores the estimated mean for the current draw in buffer for what will be the next observation on mhat. The command

postclose buffer


writes the stuff stored in buffer to the file mcs.dta. The commands

use mcs, clear
list


drop the last $$\chi^2(1)$$ sample from memory, read in the msc dataset, and list out the dataset.

Example 3 below is a modified version of example 2; I increased the number of draws and summarized the results.

Example 3: The mean of 2,000 estimated means

. set seed 12345

. postfile buffer mhat using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y])
7. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256


The average of the $$2,000$$ estimates is an estimator for the mean of the sampling distribution of the estimator, and it is close to the true value of $$1.0$$. The sample standard deviation of the $$2,000$$ estimates is an estimator for the standard deviation of the sampling distribution of the estimator, and it is close to the true value of $$\sqrt{\sigma^2/N}=\sqrt{2/500}\approx 0.0632$$, where $$\sigma^2$$ is the variance of the $$\chi^2(1)$$ random variable.

Including standard errors

The standard error of the estimator reported by mean is an estimate of the standard deviation of the sampling distribution of the estimator. If the large-sample distribution is doing a good job of approximating the sampling distribution of the estimator, the mean of the estimated standard
errors should be close to the sample standard deviation of the many mean estimates.

To compare the standard deviation of the estimates with the mean of the estimated standard errors, I modify example 3 to also store the standard errors.

Example 4: The mean of 2,000 standard errors

. set seed 12345

. postfile buffer mhat sehat using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y]) (_se[y])
7. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256
sehat |      2,000    .0629644    .0051703   .0464698   .0819693


Mechanically, the command

postfile buffer mhat sehat using mcs, replace


makes room in buffer for the new variables mhat and sehat, and

post buffer (_b[y]) (_se[y])


stores each estimated mean in the memory for mhat and each estimated standard error in the memory for sehat. (As in example 3, the command postclose buffer writes what is stored in memory to the new dataset.)

The sample standard deviation of the $$2,000$$ estimates is $$0.0625$$, and it is close to the mean of the $$2,000$$ estimated standard errors, which is $$0.0630$$.

You may be thinking I should have written “very close”, but how close is $$0.0625$$ to $$0.0630$$? Honestly, I cannot tell if these two numbers are sufficiently close to each other because the distance between them does not automatically tell me how reliable the resulting inference will be.

Estimating a rejection rate

In frequentist statistics, we reject a null hypothesis if the p-value is below a specified size. If the large-sample distribution approximates the finite-sample distribution well, the rejection rate of the test against the true null hypothesis should be close to the specified size.

To compare the rejection rate with the size of 5%, I modify example 4 to compute and store an indicator for whether I reject a Wald test against the true null hypothesis. (See appendix III for a discussion of the mechanics.)

Example 5: Estimating the rejection rate

. set seed 12345

. postfile buffer mhat sehat reject using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         quietly test _b[y]=1
7.         local r = (r(p)<.05)
8.         post buffer (_b[y]) (_se[y]) (r')
9. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256
sehat |      2,000    .0629644    .0051703   .0464698   .0819693
reject |      2,000       .0475     .212759          0          1


The rejection rate of $$0.048$$ is very close to the size of $$0.05$$.

Done and undone

In this post, I have shown how to perform an MCS of an estimator in Stata. I discussed the mechanics of using the post commands to store the many estimates and how to interpret the mean of the many estimates and the mean of the many estimated standard errors. I also recommended using an estimated rejection rate to evaluate the usefulness of the large-sample approximation to the sampling distribution of an estimator for a given DGP and sample size.

The example illustrates that the sample average performs as predicted by large-sample theory as an estimator for the mean. This conclusion does not mean that my friend's concerns about outliers were entirely misplaced. Other estimators that are more robust to outliers may have better properties. I plan to illustrate some of the trade-offs in future posts.

Appendix I: Repeating a process many times

This appendix provides a quick introduction to local macros and how to use them to repeat some commands many times; see [P] macro and [P] forvalues for more details.

I can store and access string information in local macros. Below, I store the string hello" in the local macro named value.

local value "hello"


To access the stored information, I adorn the name of the local macro. Specifically, I precede it with the single left quote () and follow it with the single right quote ('). Below, I access and display the value stored in the local macro value.

. display "value'"
hello


I can also store numbers as strings, as follows

. local value "2.134"
. display "value'"
2.134


To repeat some commands many times, I put them in a {\tt forvalues} loop. For example, the code below repeats the display command three times.

. forvalues i=1/3 {
2.    display "i is now i'"
3. }
i is now 1
i is now 2
i is now 3


The above example illustrates that forvalues defines a local macro that takes on each value in the specified list of values. In the above example, the name of the local macro is i, and the specified values are 1/3=$$\{1, 2, 3\}$$.

Appendix II: Accessing estimates

After a Stata estimation command, you can access the point estimate of a parameter named y by typing _b[y], and you can access the estimated standard error by typing _se[y]. The example below illustrates this process.

Example 6: Accessing estimated values

. drop _all

. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345

. generate y = rchi2(1)

. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------

. display  _b[y]
.91076444

. display _se[y]
.05486467


Appendix III: Getting a p-value computed by test

This appendix explains the mechanics of creating an indicator for whether a Wald test rejects the null hypothesis at a specific size.

I begin by generating some data and performing a Wald test against the true null hypothesis.

Example 7: Wald test results

. drop _all

. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345

. generate y = rchi2(1)

. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------

. test _b[y]=1

( 1)  y = 1

F(  1,   499) =    2.65
Prob > F =    0.1045


The results reported by test are stored in r(). Below, I use return list to see them, type help return list for details.

Example 8: Results stored by test

. return list

scalars:
r(drop) =  0
r(df_r) =  499
r(F) =  2.645393485924886
r(df) =  1
r(p) =  .1044817353734439


The p-value reported by test is stored in r(p). Below, I store a 0/1 indicator for whether the p-value is less than $$0.05|0 in the local macro r. (See appendix II for an introduction to local macros.) I complete the illustration by displaying that the local macro contains the value \(0$$.

. local r = (r(p)<.05)
. display "r'"
0

Categories: Programming Tags:

## Retaining an Excel cell’s format when using putexcel

In a previous blog entry, I talked about the new Stata 13 command putexcel and how we could use putexcel with a Stata command’s stored results to create tables in an Excel file.

After the entry was posted, a few users pointed out two features they wanted added to putexcel:

1. Retain a cell’s format after writing numeric data to it.
2. Allow putexcel to format a cell.

In Stata 13.1, we added the new option keepcellformat to putexcel. This option retains a cell’s format after writing numeric data to it. keepcellformat is useful for people who want to automate the updating of a report or paper.

To review, the basic syntax of putexcel is as follows:

putexcel excel_cell=(expression) … using filename[, options]


If you are working with matrices, the syntax is

putexcel excel_cell=matrix(expression) … using filename[, options]


In the previous blog post, we exported a simple table created by the correlate command by using the commands below.

. sysuse auto
(1978 Automobile Data)

. correlate foreign mpg
(obs=74)

|  foreign      mpg
-------------+------------------
foreign |   1.0000
mpg |   0.3934   1.0000

. putexcel A1=matrix(r(C), names) using corr


These commands created the file corr.xlsx, which contained the table below in the first worksheet.

As you can see, this table is not formatted. So, I formatted the table by hand in Excel so that the correlations Read more…

Categories: Programming Tags: