* fastreg.ado
* Save this file in your personal ado directory as fastreg.ado
* Type adopath in Stata to find the correct directory
* Please let me know if you find any errors, or if you have made
* any improvements that you would like to share
* Copyright Paul Geertsema 2014
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see http://www.gnu.org/licenses/.
* CODE START
* USAGE: fastreg varlist [in], [CALCVCV] [NOMISSINGTEST] [NOCONST]
* MATA: set up a mata function that calculates OLS estimates
* version 13
mata:
stata("capture mata mata drop m_fastreg()")
void m_fastreg(string scalar depvar, string scalar regs,
            real scalar startob, real scalar endob,
            real scalar calcvcv, real scalar nomissingtest,
            real scalar noconst)
{
   real matrix X, XX, Xy, V
   real vector beta, e, s2
   real scalar N, k, dof
   st_view(y, (startob, endob), tokens(depvar))
   st_view(X, (startob, endob), tokens(regs))
   if (nomissingtest == 0) {
      missing = rowmissing((y, X))
      X = select(X, missing:==0)
      y = select(y, missing:==0)
}
   if (noconst == 0) {
      XX = quadcross(X,1,X,1)
      Xy = quadcross(X,1,y,0)
      N = rows(X)
      k = cols(X)
      if(rank(XX) < k+1) error(2001)
}
else {
      XX = quadcross(X,X)
      Xy = quadcross(X,y)
      N = rows(X)
      k = cols(X)
      if(rank(XX) < k) error(2001)
}
   beta = (invsym(XX)*Xy)'
   st_matrix("r(beta)", beta)
   st_numscalar("r(N)", N)
   if (calcvcv != 0) {
      if (noconst == 0) {
         e = y - (X, J(N,1,1)) * beta'
		 dof = N-k-1
		 }
      else {
         e = y - X * beta'
         dof = N-k
}
      s2 = quadcross(e,e)/(dof)
      V = s2*cholinv(XX)
      st_matrix("r(V)", V)
      st_numscalar("r(dof)", dof)
} }
end
* STATA: add stata wrapper to mata function
capture program drop fastreg
program fastreg, eclass
* version 13
   syntax varlist [in/], [CALCVCV] [NOMISSINGTEST] [NOCONST]
   if "`in'" != "" {
      * extract starting and ending observation from -in- range
      gettoken startob in: in, parse("/")
      gettoken dummy in: in, parse("/")
      local endob "`in'"
}
else {
      * full dataset
      local startob = 1
      local endob = _N
}
   local depvar: word 1 of `varlist'
   local regs: list varlist - depvar
   * nomissing set -> user guarantees there are no missing obs in sample
   local bool_nomissingtest = cond("`nomissingtest'" == "", 0, 1)
   * noconst set -> regression without intercept
   local bool_noconst = cond("`noconst'" == "", 0, 1)
   * calcvcv set -> if specified,
   * VCV is calculated -> allows statistical inference
   local bool_calcvcv = cond("`calcvcv'" == "", 0, 1)
   if `bool_noconst' == 0 local vnames `regs' _cons
   else local vnames `regs'
   * call mata regression routine
   mata: m_fastreg("`depvar'", "`regs'", `startob', `endob', /*
   */ `bool_calcvcv', `bool_nomissingtest', `bool_noconst')
   tempname b
   mat `b' = r(beta)
* if no VCV required...
   if `bool_calcvcv' == 0 {
      matname `b' `vnames', c(.)
      ereturn post `b'
   }
   else {
      tempname V
      tempname dof
      mat `V' = r(V)
      local dof = r(dof)
      local N = r(N)
      matname `V' `vnames'
      matname `b' `vnames', c(.)
      ereturn post `b' `V', depname(`depvar') obs(`N') dof(`dof')
}
end
* CODE END