Stata - Write your own "fast" bootstrap

  * Imagine we would like to bootstrap the standard errors of a command using a bootstrap routine.

  * I created a previous post demonstrating how to write a bootstrap command.  This is a similar post however the bootstrap is much faster than the previous one.
 
  * First let's generate some data.
 
  clear
  set obs 1000
  gen x1 = rnormal()
  gen x2 = 2*rnormal()
  gen u = 5*rnormal()
 
  gen y = 5 + x1 + x2 + u
 
  local dependent_var x1 x2
 
  local command_bs reg y `dependent_var'
 
  * First let's see how the base command works directly.
  `command_bs'
 
  * As a matter of comparison this is the built in bootstrap command.
  bs, rep(100): `command_bs'
 
  * The following code is yet another user written bootstrap alternative.
 
  * I wrote this
 
  * Specify the number of bootstrap iterations
  local bs = 100
 
  * Save the number of observations to be drawn
  local N_obs = _N
 
  * Number of terms plus one for the constant of coefficients to be saved
  local ncol = wordcount("`dependent_var'")+1
 
  mata theta = J(`bs',`ncol',0)
  forv i = 1(1)`bs' {
    * Preserve the initial form of the data
preserve

    * Draw the indicators of the resample
mata draw = ceil(runiform(`N_obs',1):*`N_obs')

* Create an empty matrix to hold the number of items to expand
mata: expander = J(`N_obs',1, 0)

* Count the number of items per observation to generate.
    qui mata: for (i=1 ; i <= `N_obs'; i++) expander[i]=sum(draw:==i) ; "draws complete"

* Pull the expander value into stata
getmata expander=expander

* Drop unnessessary data
qui drop if expander == 0

* Expand data, expand==1 does nothing
    qui expand expander

* Run a regression
    qui `command_bs'

* Send to mata the matrix of results
mata theta[`i',] = st_matrix("e(b)")
   
* Configure the visual display
di _c "."
    if (int(`i'/10)==`i'/10) di _c "|"
    if (int(`i'/50)==`i'/50) di " " `i'

restore
  }
 
  * The estimates of the coefficients have been saved into a theta matrix.
  mata theta
 
  * Now let's calculate the standard deviations.
  mata bs_var = variance(theta)
  mata bs_var
  mata bs_se = diag(bs_var):^.5
  mata bs_se

Matrix Operations in Mata

* This post demonstrates a few methods for how to input matrices into Mata and how to do some basic matrix operations.  Mata is a matrix programming language so basic matrix operations are extremely easy.

mata
// Let's first build some matrices in Mata
// They can be built directly
A = ( 2 , -1, 5 \ 3, 0 , -1 \ 3, 3 , -1)
A

B1 = (0 , 3 ,-1)
B2 = (3 , 0 , 0)
B3 = (0 , 2 , 0)

// Or through a combination of vectors
B = (B1 \ B2 \ B3)

B

// We can also start with an empty matrix and fill in values

C = J(3,2,4)
// The J command creates a matrix with 3 rows, 2 colums, and with default values of 4

C
// We can replace individual elements once the matrix is created
C[1,2] = 3
C[3,1] = 7

// Or entire submatrices
C[2,] = (2,6)
C

// Now let's see how various matrix operations perform in Mata:

// a. AB
A*B

// b. BA
B*A

// c. A+B
A+B

// d. A'B'
A'*B'

// e. B'A'
B'*A'

// f. A'B
A'*B

// g. (AB)'
(A*B)'

// h. ABC
A*B*C

// i. C'AC
C'*A*C

// j. CAC'
C*A*C'
// This does not exist because C*A cannot be multiplied as C is 3x2 and A is 3x3

// k. trace(C'AC)
trace(C'*A*C)

// l. trace(CAC')
trace(C*A*C')
// Likewise the trace does not exist

end

* Matrices can also be input into Mata from data sets.

clear
set obs 5
gen y = 3
replace y = 6  if _n == 2
replace y = 10 if _n == 3
replace y = 8  if _n == 4
replace y = 2  if _n == 5

mata
// The command st_data retrieves data from stata variables to be used in Mata
y = st_data(. , "y")
y

// Once input, matrices can easily be manipulated
y*y'

// The square of the norm of y
norm(y)^2

// This happens to be identical to:
sum(y:^2)

// The norm is the equclidean norm which is the square square root of the sum of all of the squares of a vector.
// Thus the square of it is just the sum of the squares.

// As should be clear, manipulating matrices in Mata is extremely easy.
// Thus Stata is able to pack a powerful Matrix programming language inside an effective high level user language.

// The largest frustration that I have had with Mata is the relative quality of the documentation.
// I find the documentation of Mata much harder to use than that of Stata (at least in version 11, perhaps version 12 has better documentation).
end

Mata speed gains over Stata

* The inclusion of Mata as an available alternative programming language for Stata users was a great move by Stata.

* Mata in general runs much quicker than programming on the surface level in Stata.

* In Stata each loop that runs is compiled (interpretted into machine code) as it runs creating a lot of work for the machine.

* In Mata on the other hand, the entire loop is compiled prior to running.

* Let's see how this works.

* Let's say we want to add up the square of the numbers 1 through 100000

* Method 1: Surface loop

timer clear 1
timer on 1
local x2 = 0

forv i = 1/1000000 {
  local x2 = `x2'+`i'^2
}

di `x2'

timer off 1
timer list 1

* On my laptop, this takes about 13.5 seconds

* Method 2: Mata loop
timer clear 1
timer on 1
mata
  x2=0
  // This command can be read as start i at 1,
  // keep looping so long as i is less than 1000000,
  // the third argument looks a little fishy but it is syntax
  // that has been around for a while (at least since C).
  // It would be identical to writing i=i+1, in other words, add 1 to i.
  // Following the for loop we can immediately place a since line command.
  for (i = 1; i <= 1000000; i++) x2=x2+i^2
  // If there is nothing done with the value x2 then mata displays this value.
  // R handles this identically
  x2
end

timer off 1
timer list 1

* In contrast, my computer completed the loop using mata in .27 seconds, many magnitudes of speed faster.

* However this does not mean you need to learn to use mata (since it has its own limitations and syntax) in order to speed up your commands.

* Method 3: Use Stata's data structure to accomplish vector tasks
timer clear 1
timer on 1

clear
set obs 1000000
gen x2 = _n^2

* The sum command will calculate the mean of x2 which is the same as the sum of x2 divided by it's number of observations.
sum x2
* We can reverse that operation easily.
di r(N)*r(mean)

timer off 1
timer list 1
* Using a little knowledge of how Stata stores post command information this method does the same trick in .2 seconds

* Method 4: The speed gains in 3 was as a result of using the vector structure of data columns.  Mata can do very similar things even easier.

timer clear 1
timer on 1
// This command looks a little fishy, but it is easy to understand.
// Order of operations must be taken into account.
// First the 10^6 is evaluated which equals 1000000
// Then the vector 1..10^6 is made which looks like 1 2 3 ... 1000000
// The .. tells mata to make a count vector.
// If I had written :: then mata would have made a column vector instead.
// Once the vector is made then the command :^2 tells stata to do a piece wise squaring of each term in the vector.
// Finally the sum command adds all of the elements of the vector together to generate the result we were looking for.
mata: sum((1..10^6):^2)
timer off 1
timer list 1
* The result is that this command only took .04 seconds to run through efficient coding in Mata.

# As a matter of comparison, this command
system.time(sum((1:10^6)^2))
# took .04 seconds in R

# And the loop:
x=0
system.time(for(i in 1:10^6) x=x+i^2)
# 1.3 seconds

# Thus Mata in this example is significantly faster than Stata and about the same speed as R.

Mata: program your own IV Estimator

* Mata: program your own IV Estimator

* First let see how we could program an IV estimator without matrices.
version 11

clear
set obs 1000

gen z = rnormal()

gen u= rnormal()

gen x = u + z + rnormal()

gen y = 2 + 2*x + 15*u

reg y x

foreach v in x y z {
  foreach vv in x y z {

  gen `v'`vv' = `v'*`vv'
  qui sum
  local `v'`vv' = r(mean)*r(N)
  }
}

di "IVhat = " (`xz'*(`zz')^-1*`zx')^-1 * (`xz'*(`zz')^-1*`zy')

ivreg y (x=z)
* Looks pretty close.  However, there is some divergence due to the above estimator not allowing for a constant.

ivreg y (x=z), nocon

* Now let's see how to do it in Mata.

mata

 x = st_data(.,("x"))
 z = st_data(.,("z"))
 y = st_data(.,("y"))

 IVhat1 = invsym((x'*z)*invsym(z'*z)*(z'*x))*(x'*z)*invsym(z'*z)*(z'*y)

"IV estimate without constant"
 IVhat1

"Create a constant vector of length x"
 cons = J(rows(x), 1, 1)

 X = (x, cons)
 Z = (z, cons)

 IVhat2 = invsym((X'*Z)*invsym(Z'*Z)*(Z'*X))*(X'*Z)*invsym(Z'*Z)*(Z'*y)

"IV estimate with constant"
 IVhat2

end
* End is necessary with mata do files to indicate when the mata subsystem is completed.

Write your own System OLS in Mata

* This is a method related to Seemingly Unrelated Regression (SUR) referred to in a previous post.

* Imagine that you have three different equations that have different dependent variables but the same explanatory variables.

* Let us think of the example of the demand for capturing and mining an astroid and the research that goes into it. See Tony Cookson's post on the subject

* In general we may think that might have three different relationships that we are trying to uncover but we don't understand the relationship between them.

* In equation 1: we have the price of platnium
* In equation 2: we have millions of dollars of investments in private spacecraft
* In equation 3: we have an index prepresenting the probability of making a launch attempt

* This will specify the means of the three variables to be drawn jointly
matrix m = (0,0,0)

* This will specify the covariance matrix for the errors
matrix C = (12, -5, 5 \ -5, 12, -3 \ 5, -3, 6)
matrix list C

* Clear the memory and tell stata to prepare for 1000 observations
clear
set obs 100

* Draw three variablesb
drawnorm u1 u2 u3, means(m) cov(C)

* Now let's generate our exogenous variables
gen platnum_Q = 10 + rnormal()
  label var platnum_Q "The current quantity of plantinum available on the market"
gen invest_age = 50 + rnormal()
  label var invest_age "Average age of investors"
gen national_news_coverage = rpoisson(10)
  label var national_news_coverage "The current amount of national news coverage of space crafts."

* Now generate the dependent variables
gen platnum_price = 3*platnum_Q-0*invest_age-.1*national_news_coverage+u1*20
gen space_investment = 2*platnum_Q+3.5*invest_age +5.2*national_news_coverage+u2*5
gen launch_index = 6.01*platnum_Q+5.3*invest_age-10*national_news_coverage+u3*10

* send variables to mata.

putmata x1=platnum_Q x2=invest_age x3=national_news_coverage, replace
putmata y1=platnum_price y2=space_investment y3=launch_index, replace

mata
N=rows(x1)
cons = J(N, 1, 1)

// Save the explanatory variables to a matrix
X=(x1 , x2 , x3, cons)

// Save the dependent variables to a matrix
// Notice that at this point this is the same as OLS except that in OLS there is only one Y.
Y=(y1 , y2 , y3)

// alternatively the Sytem OLS estimator is:
B_OLS = invsym(X'*X) * (X'*Y)
B_OLS

end

reg platnum_price platnum_Q invest_age national_news_coverage
reg space_investment platnum_Q invest_age national_news_coverage
reg launch_index platnum_Q invest_age national_news_coverage

* This is the same as using seemingly unrelated regression.
mvreg platnum_price  space_investment launch_index = platnum_Q invest_age national_news_coverage

* However, seemlingl unrelated regression allows for restrictions on coefficients.
sureg (platnum_price=platnum_Q invest_age national_news_coverage) ///
     (space_investment=platnum_Q invest_age national_news_coverage) ///
(launch_index=platnum_Q invest_age national_news_coverage)