program define fitpower
version 2.1
*Fit a power curve (log-log) to untransformed data
capture drop h_1
capture drop h_2
capture drop h_3
capture drop h_4
capture drop resid
mac def _varlist "req ex "
mac def _if "optional"
mac def _options "*"
parse "%_*"
parse "%_varlist", parse(" ")
di " "
di "Enter an X-axis label (<20 characters),"
di "     OR press <Enter> to use " in red "   %_2  " in white " as the label."
di "X-label..." _r(xn)
di " "
if "%xn"=="" {mac def xn "%_2"}
mac def xnm "b2("%xn")"
di " "
di "Enter a Y-axis label (<20 characters),"
di "     OR press <Enter> to use " in red "   %_1   " in white " as the label."
di "Y-label..." _r(yn)
di " "
if "%yn"=="" {mac def yn "%_1"}
macro define ynm "l1(" ") l2("%yn")"
quietly{
* Calculate SS Total of raw data
              summarize %_1
              mac def sst=_result(4)*(_result(1)-1)
* Log-transform and regress
      gen h_1=log(%_1)
      gen h_2=log(%_2)
      reg h_1 h_2 %_if
      predict h_3
      gen h_4=exp(h_3)
* Calculate R-square
                  gen resid=(%_1-h_4)^2
                   replace resid=sum(resid)
                  mac def ssres=resid[_N]
                  mac def ssreg=%sst-%ssres
                  mac def r2=%ssreg/%sst
                  mac def R2=1-((1-%r2)*(_N-1)/_result(5))
                  mac def F=((%ssreg)/_result(3))/(%ssres/_result(5))
* di "sst=" %sst
* di "ssresidual=" %ssres
* di "ssreg=" %ssreg
* di "r2="%r2
* di "R-square=" %R2
* di "F=" %F
* more
          macro define a=exp(_coef[_cons])
            mac def a = substr(string(%a),1,5)
          macro define b=_coef[h_2]
            mac def b= substr(string(%b),1,5)
  macro define d "t1("%yn = %a %xn ^ %b") b1("Power Curve Regression")"
macro define p=fprob(_result(3),_result(5),%F)
  if %p<.001 { mac def p = .0001 }
    mac def p = substr(string(%p),1,4)
    mac def R2 = substr(string(%R2),1,4)
macro define e "t2("Adjusted R-Squared= %R2 <[][]> p= %p")"
gr %_1 h_4 %_2 %_if, s(oi) sort c(.s) xlab ylab %e %d %xnm %ynm %_options
if "%_options"=="" {
 display" "
  display "   To print graph, run fitpower again specifying a filename, e.g.,"
  display "         fitpower %_1 %_2, saving(jt6)"
drop h_1 h_2 h_3 h_4
end