* Typical (easy) power analysis relies upon drawing from known distributions often with differences in means using estimators which are unbiased and for which variances are easy to estimate.

* However, many econometric procedures are not unbiased but rather instead only instead consistent and have likewise difficult to estimate variances.

* In order to estimate power under these circumstances simulations can be useful.

* Andreas Tiemann recently informed me of this phenomena citing an interesting article by Benjamin Arnold, Daniel Hogan, John Colford, and Alan Hubbard 2011(http://www.biomedcentral.com/1471-2288/11/94/abstract).

* This post will estimate the sample size required to reduce the chance of false rejection to 5% (Type I error) and have a minimum chance of correct rejection of 80% (Type II error).

* I will use the example recently posted on Dave Giles blog October 30 on Non-linear Least Squares (http://davegiles.blogspot.com/2012/10/some-properties-of-non-linear-least.html#more)

* This code is similar to code shared with me by Andreas Tiemann and Juan José Matta.

cap program drop NLLS_sim

program define NLLS_sim, rclass

clear

* This first argument is the number of observations

set obs `1'

* I will generate x2 as a uniform

gen x2 = runiform()

* Generate the error

gen e = rnormal()

* Generate the parameters to estimate (these are the second and third arguments following the command name.

local beta1 = `2'

local beta2 = `3'

gen y = `beta1' + (`beta2'*x2)^`beta1' + e

* Run the non-linear least squares estimation. Note: I set the initial value of beta1 to .5 because it was having difficulty converging.

nl (y = {beta1=1} + ({beta2}*x2)^{beta1})

* I don't know any way of directly extracting the p values but they can be recomputed

local p1 = 2 * ttail(e(df_r), abs(_b[/beta1]/_se[/beta1]))

* Return p values and significance levels at 5%

cap return scalar p1 = `p1'

cap return scalar sig1 = `p1'<.05

local p2 = 2 * ttail(e(df_r), abs(_b[/beta2]/_se[/beta2]))

cap return scalar p2 = `p2'

cap return scalar sig2 = `p2'<.05

end

* Now let's test out our program

NLLS_sim 10 1 2

return list

* Looking good. The p values are the same and the significance is properly rejecting

* First I would like to warn that this can take a while.

local sig = 1

local n=2

* First let's estimate the likelihood of falsely rejecting the null for any sample size.

while `sig'>.5 {

local n=`n'+1

* Number of observations:

local nobs = ceil(10^(`n'/3))

simulate sig2=r(sig2), reps(100): NLLS_sim `nobs' 1 0

qui: sum sig2

local sig=r(mean)

di "We reject the null `sig2'% of the time where beta2=0"

}

local sig = 0

* We can let n start at the value from the previous simulation since

* Now let's estimate the sample size required to correctly reject the null 80% of the time.

while `sig'<.8 {

* Number of observations:

local nobs = ceil(10^(`n'/3))

simulate sig2=r(sig2), reps(100): NLLS_sim `nobs' 1 1

qui: sum sig2

local sig=r(mean)

di "We reject the null `sig2'% of the time where beta2=1"

local n=`n'+1

}

## No comments:

## Post a Comment