Oleh Komashko 0.93 2026-07-19 advanced treatment of marginal effects and elasticities C20 C21 pdfdoc:a_eff.pdf

a discrete series

uni = values(y)
n = nelem(uni)
series ret = y>=uni[1]
loop i=2..n -q
  ret = ret+(y>=uni[i])
endloop
return ret

#gn = g_varnumber()
gretl_version = $version

list z = seq(1,$nvars - 1)
gna = strsub(varname(z),","," ")
gnas = strsplit(gna)
ret = ""
ret2 = 0
#string line
loop i=1..($nvars - 1) -q
  isdu = isduname(gnas[i])
  isdi = isdiname(gnas[i])
  isinti = chek_intname(gnas[i])
  #         print isinti
  #         funcerr "14"
  if (isdu||isdi)&&(!isinti)
    isdum = isdummy(z[i])
    if !isdum
      ret = ret~" "~gnas[i]
    endif
    if isdu&&isdum
      znam = getdiname (corr_duname (gnas[i]))
      fifi = findind(znam,gna)
      if fifi==0
        gnasi = gnas[i]
        printf "\n%s is a dummy, the name looks like a level dummy name,\n", gnasi
        printf "but variable %s is absent in the dataset\n", znam
        ret2 = ret2+1
      endif
    endif
    if !(isdu)&&isdi&&isdum
      znam = getdiname (gnas[i])
      fifi = findind(znam,gna)
      if fifi==0
        gnasi = gnas[i]
        printf "\n%s is a dummy, the name looks like a level dummy name,\n", gnasi
        printf "but variable %s is absent in the dataset\n", znam
        ret2 = ret2+1
      endif
    endif
  endif
  if isinti
    #             outfile "@dotdir\inner666.txt" --write
    #                 labels z[i]
    #             outfile --close
    #             string s = readfile("@dotdir\inner666.txt")
    #             remove("@dotdir\inner666.txt")
    #             scalar j = 1
    #             loop while getline(s, line)&&(j<3) -q
    #                 uu=sprintf( "%s\n",  line)
    #                 j++
    #             endloop
    #             bign = strlen(uu)
    #             uuu = strstr(uu,":")
    #             n = strlen(uuu)
    #             vvv = strstrip(substr(uu,1,bign-n))
    #             uuu = substr(uuu,3,n)
    #             uuu = strstrip(uuu)
    bundle bui = getinfo(i)
    fififi = !inbundle(bui,"description")
    vvv = bui.name
    uuu = fififi? "" : bui.description
    ul = uniqvarlist(uuu)
    inti = substr(uuu,1,14)=="interaction of"
    inti2 = isstrstr(uuu,"==")
    intis = strsplit3(uuu,"==")
    nintis = nelem(intis)
    inti3 = (nintis==2)
    inti4 = 0
    if inti3
      intis2 = intis[2]
      nintis2 = strlen(intis2)
      cknin1 = substr(intis2,nintis2,nintis2)==")"
      cknin2 = isnat(strstrip(substr(intis2,1,nintis2-1)))
      inti4 = cknin1*cknin2
    endif
    inti = inti*inti2*inti3*inti4>0
    if inti == 0
      gnasi = gnas[i]
      printf "\nThe name '%s' looks like an interaction term name,\n", gnasi
      printf "but %s is not an interaction term.\n", gnasi
      ret2 = ret2+1
    endif
  endif
endloop
ret = strstrip(ret)
if (ret=="")&&(ret2==0)
  printf "\nAll variables names are OK!\n"
else
  if ret !=""
    printf "\nThe names of the following variables looks as names\n"
    printf "of dummies, but they are not dummies:"
    print_long_string(ret)
  endif
  nn = nelem(strsplit(ret))+ret2
  if nn==1
    printf "\nThere was 1 warning.\n"
  else
    printf "\nThere were %d warnings.\n", nn
  endif
endif

b is $model in most cases the model formula at and wrt vars names, eff type values of 'at vars' for zlist containing models dataset, for gretl before 2018d

gretl_version = $version
isnuatwrt = !exists(at_wrt_type)
if isnuatwrt
  string at_wrt_type = ""
endif
strings all_string_input = fndtype(at_wrt_type)
string link = all_string_input[3]
string wrt = all_string_input[1]
string wrtBACKUP = wrt
string at_nams =  all_string_input[2]
cmd = b.command
cmd2 = cmd
name_of_y = b.depvar
matrix coe = b.coeff
matrix vcov = b.vcv
list xx = b.xlist
exifla = exists(formula)
if !exifla
  formula = smart_formula2(xx)
else
  isdotcomma = instring(formula,";")
  #print isdotcomma
  if !isdotcomma
    formula_backup = formula
    formula_2s = strsplit(formula)
    nf2s = nelem(formula_2s)
    list indep0 = null
    loop i=1..nf2s -q
      jjj = inlist(xx,formula_2s[i])
      if jjj > 0
        indep0 += xx[jjj]
      endif
    endloop
    #eval indep0
    formula = smart_formula2(xx,indep0)
  endif
endif
series sample = b.sample
formulas = strsplit3(formula,";")
newformula = formulas[1]
if !exifla
  newformula = strstrip(smartstrsub(newformula,"const",""))
endif
clistnames = formulas[2]
dilistnames = formulas[3]
list clist = null
list dilist = null
nvar = $nvars
stri = mk_txt_range(nvar - 1)
list da = @stri
if nelem(da) == 0
  list da = seq(1,$nvars)
endif
if clistnames != "_empty"
  clistnamess = strsplit(clistnames)
  nclist = nelem(clistnamess)
  loop i=1..nclist -q
    clnamsi = clistnamess[i]
    clist+=da[inlist(da,clnamsi)]
  endloop
endif
if dilistnames != "_empty"
  dilistnamess = strsplit(dilistnames)
  ndilist = nelem(dilistnamess)
  loop i=1..ndilist -q
    dilnamsi = dilistnamess[i]
    dilist+=da[inlist(da,dilnamsi)]
  endloop
endif
isnumat = !exists(at_vals)
smpl sample --dummy
string orig_names = clistnames == "_empty"? "" : clistnames
orig_names = dilistnames == "_empty"? orig_names : orig_names~" "~dilistnames
if cmd == "heckit"
  scalar grvers = $version
  string at_namsZ = ""
  string wrtZ = ""
  list zlist = b.zlist - b.zlist[1]
endif
if (!isnuatwrt)&&(cmd == "heckit")
  string at_namsALL = at_nams
  at_nams = str_setinrersect(at_nams,orig_names)
  at_namsZ = str_setdiff(at_namsALL,at_nams)
  string wrtALL = wrt
  wrt = str_setinrersect(wrt,orig_names)
  string wrtZ = str_setdiff(wrtALL,wrt)
endif
#smpl sample --dummy
isnw = wrt == ""
#string orig_names = clistnames == "_empty"? "" : clistnames
#orig_names = dilistnames == "_empty"? orig_names : orig_names~" "~dilistnames
if !isnw
  wtwt = str_setdiff(wrt,orig_names)
  wrtchk = strlen(wtwt)
  errtext = sprintf("wrong name(s) in 'wrt' list: %s",wtwt)
  if wrtchk > 0
    funcerr errtext
  endif
endif
if strlen(at_nams) == 0
  bundle bu = defvals (newformula,clist,dilist)
else
  if nelem(strsplit(at_nams))!=nelem(at_vals)
    funcerr "inconsistent names and values for 'at' variables"
  endif
  bundle bu = findvals(at_nams,at_vals,newformula,clist,dilist)
endif
allnames = bu["allnames"]
matrix allvalues = bu["allvalues"]
contvars = bu["contvars"]
if contvars == "_empty"
  contvars = ""
else
  if !isnw
    contvars = str_setinrersect(contvars,wrt)
  endif
endif
dumvars = bu["dumvars"]
if dumvars == "_empty"
  dumvars = ""
else
  if !isnw
    dumvars = str_setinrersect(dumvars,wrt)
  endif
endif
divars = bu["divars"]
if divars == "_empty"
  divars = ""
else
  if !isnw
    divars = str_setinrersect(divars,wrt)
  endif
endif
if divars == ""
  dilevs = ""
  dialll ={}
else
  divarss = strsplit(divars)
  ndis = nelem(divarss)
  dilevs = ""
  matrix dialll=zeros(1,ndis)
  alldinam = strjoin(varnames(dilist))
  matrix vadii={}
  loop i=1..ndis -q
    dili = mk_factlevs2 (formula, divarss[i],allnames)
    indii = findind(divarss[i],alldinam)
    vadii = values(dilist[indii])
    nvadii = nelem(vadii)
    if i==1
      txt1=""
    else
      txt1=","
    endif
    dilevs = dilevs~txt1~dili
    dialll[i] = nvadii
  endloop
endif
lang = $lang
lincheck = (cmd == "ols")||(cmd == "hsk")||(cmd == "wls")||(cmd == "tsls")||(cmd == "arima")
if lincheck
  #if (cmd == "ols") && (link == "exp")
  if link == "exp"
    cmd2 = "poisson"
  else
    cmd2 = "lin"
  endif
elif cmd == "negbin"
  cmd2 = "poisson"
elif  cmd == "logistic"
  cmd2 = "logit"
elif cmd == "logit"
  catch orde = b.ordered
  errord = $error
  if !errord
    cmd2 = "orderedL"
  else
    catch mult = b.multinom
    errmul = $error
    if !errmul
      cmd2 = "multinomial"
    endif
  endif
elif cmd == "probit"
  catch orde = b.ordered
  errord = $error
  if !errord
    cmd2 = "orderedP"
  endif
endif
#print cmd cmd2
lppcheck = (cmd2 == "logit")||(cmd2 == "probit")||(cmd2 == "poisson")
############ IF BY COMMAND
if lincheck
  if link != "exp"
    printf "\n%s marginal effects on %s:\n", cmd, name_of_y
  else
    if (strlen(name_of_y) > 2)&&(substr(name_of_y,1,2)=="l_")
      new_name_y = substr(name_of_y,3,strlen(name_of_y))
      printf "\n%s marginal effects on %s:\n", cmd, new_name_y
    else
      printf "\n%s marginal effects on exp(%s):\n", cmd, name_of_y
    endif
  endif
  if (cmd == "arima")
    armnxx = nelem(xx)
    armnn = b.ncoeff
    if inlist(xx,const)
      armindi = 1~seq(armnn-armnxx+2,armnn)
    else
      armindi = seq(armnn-armnxx+1,armnn)
    endif
    coe = coe[armindi]
    vcov = vcov[armindi,armindi]
  endif
endif
if lppcheck && (!lincheck)
  if cmd == "negbin"
    ncoe = rows(coe)
    matrix vcov = vcov[1:(ncoe-1),1:(ncoe-1)]
    matrix coe = coe[1:(ncoe-1)]
  endif
  if (cmd2 == "poisson")
    printf "\n%s marginal effects on E(%s):\n", cmd, name_of_y
  elif cmd == "logistic"
    printf "\n%s marginal effects on %s:\n", cmd, name_of_y
  else
    printf "\n%s marginal effects on P(%s = 1):\n", cmd, name_of_y
  endif
endif
oneq = lppcheck||lincheck
if oneq
  if  contvars == ""
    matr1 = {}
  else
    matr1 = cont_grads_many(newformula,allnames,allvalues,contvars,coe,cmd2)
  endif
  if  dumvars == ""
    matr2 = {}
  else
    matr2 = dum_grads_many(newformula,allnames,allvalues,dumvars,coe,cmd2)
  endif
  if  divars == ""
    matr3 = {}
  else
    matr3 = {}
    ndss = nelem(divarss)
    dilevelss = strsplit3(dilevs,",")
    loop i=1..ndss -q
      matr3i = discr_grads2_many(newformula, allnames, allvalues, divarss[i],dilevelss[i],mk_alllevs(dialll[i]),coe,cmd2)
      matr3 = matr3|matr3i
    endloop
  endif
endif
if cmd == "heckit"
  name_of_yz = varnames(b.zlist)[1]
  if link != "Ey"
    printf "\n%s marginal effects on E(%s|%s > 0):\n", cmd, name_of_y, name_of_yz
  else
    printf "\n%s marginal effects on E(%s):\n", cmd, name_of_y
  endif
  list xxZ = zlist
  exiflaZ = exists(formulaZ)
  if !exiflaZ
    formulaZ = smart_formula2(xxZ)
    if exifla
      #else
      #print formula
      isdotcomma = instring(formula_backup,";")
      if !isdotcomma
        formula_2s = strsplit(formula_backup)
        nf2s = nelem(formula_2s)
        list indep0 = null
        loop i=1..nf2s -q
          jjj = inlist(xxZ,formula_2s[i])
          #indep0 += xxZ[jjj]
          if jjj > 0
            indep0 += xxZ[jjj]
          endif
        endloop
        formulaZ = smart_formula2(xxZ,indep0)
      endif
    endif
  endif
  # if !exifla
  #         formula = smart_formula2(xx)
  #     else
  #         isdotcomma = instring(formula,";")
  #
  #         if !isdotcomma
  #             formula_2s = strsplit(formula)
  #             nf2s = nelem(formula_2s)
  #             list indep0 = null
  #             loop i=1..nf2s -q
  #                 jjj = inlist(xx,formula_2s[i])
  #                 indep0 += xx[jjj]
  #             endloop
  #             formula = smart_formula2(xx,indep0)
  #         endif
  #     endif
  formulasZ = strsplit3(formulaZ,";")
  newformulaZ = formulasZ[1]
  if !exiflaZ
    newformulaZ = strstrip(smartstrsub(newformulaZ,"const",""))
  endif
  clistnamesZ = formulasZ[2]
  dilistnamesZ = formulasZ[3]
  list clistZ = null
  list dilistZ = null
  if clistnamesZ != "_empty"
    clistnamessZ = strsplit(clistnamesZ)
    nclistZ = nelem(clistnamessZ)
    loop i=1..nclistZ -q
      clnamsiZ = clistnamessZ[i]
      clistZ+=da[inlist(da,clnamsiZ)]
    endloop
  endif
  if dilistnamesZ != "_empty"
    dilistnamessZ = strsplit(dilistnamesZ)
    ndilistZ = nelem(dilistnamessZ)
    loop i=1..ndilistZ -q
      dilnamsiZ = dilistnamessZ[i]
      dilistZ+=da[inlist(da,dilnamsiZ)]
    endloop
  endif
  string orig_namesZ = clistnamesZ == "_empty"? "" : clistnamesZ
  orig_namesZ = dilistnamesZ == "_empty"? orig_namesZ : orig_namesZ~" "~dilistnamesZ
  at_namsZ = str_setinrersect(orig_namesZ,at_namsZ)
  if strlen(at_namsZ) == 0
    bundle buZ = defvals (newformulaZ,clistZ,dilistZ)
  else
    #####?????????????   at_valsZ
    matrix findiZZ = findinds(at_namsZ, at_namsALL)
    #print findiZZ
    matrix at_valsZ = at_vals[findiZZ]
    bundle buZ = findvals(at_namsZ,at_valsZ,newformulaZ,clistZ,dilistZ)
  endif
  allnamesZ = buZ.allnames
  contvarsZ = buZ.contvars
  if contvarsZ == "_empty"
    contvarsZ = ""
  endif
  dumvarsZ = buZ.dumvars
  if dumvarsZ == "_empty"
    dumvarsZ = ""
  endif
  matrix allvaluesZ = buZ.allvalues
  divarsZ = buZ.divars
  if divarsZ == "_empty"
    divarsZ = ""
  endif
  #         if (strlen(wrt) > 0)||(strlen(wrtZ) > 0)
  #             contvars = str_setinrersect(contvars,wrt)
  #             contvarsZ = str_setinrersect(contvarsZ,wrtZ)
  #             dumvars = str_setinrersect(dumvars,wrt)
  #             dumvarsZ = str_setinrersect(dumvarsZ,wrtZ)
  #             divars = str_setinrersect(divars,wrt)
  #             divarsZ = str_setinrersect(divarsZ,wrtZ)
  #         endif
  if (strlen(wrtBACKUP) > 0)
    contvars = str_setinrersect(contvars,wrtBACKUP)
    contvarsZ = str_setinrersect(contvarsZ,wrtBACKUP)
    dumvars = str_setinrersect(dumvars,wrtBACKUP)
    dumvarsZ = str_setinrersect(dumvarsZ,wrtBACKUP)
    divars = str_setinrersect(divars,wrtBACKUP)
    divarsZ = str_setinrersect(divarsZ,wrtBACKUP)
  endif
  #eval "IS ANY CONTWARS?"
  if  (strlen(contvars)+ strlen(contvarsZ)) == 0
    matr1 = {}
  else
    matr1 = all_hk_fun(newformula,allnames,contvars,newformulaZ, allnamesZ, contvarsZ,b.base_coeffs+1,b.ncoeff,allvalues,allvaluesZ,coe,1,link)
    #print matr1
    second = str_setdiff(contvarsZ,contvars)
    contvars = strstrip(contvars~" "~second)
  endif
  #####################
  #print dumvars dumvarsZ
  if  (strlen(dumvars)+ strlen(dumvarsZ)) == 0
    matr2 = {}
  else
    matr2 = all_hk_fun(newformula,allnames,dumvars,newformulaZ, allnamesZ, dumvarsZ,b.base_coeffs+1,b.ncoeff,allvalues,allvaluesZ,coe,2,link)
    # print matr2
    secondd = str_setdiff(dumvarsZ,dumvars)
    dumvars = strstrip(dumvars~" "~secondd)
  endif
  #  eval "divars"
  #                 print divars
  ################
  if divarsZ == ""
    dilevsZ = ""
    dialllZ ={}
  else
    divarssZ = strsplit(divarsZ)
    ndisZ = nelem(divarssZ)
    dilevsZ = ""
    matrix dialllZ=zeros(1,ndisZ)
    alldinamZ = strjoin(varnames(dilistZ))
    matrix vadiiZ={}
    loop i=1..ndisZ -q
      #############?????????????????
      diliZ = mk_factlevs2 (formulaZ, divarssZ[i],allnamesZ)
      indiiZ = findind(divarssZ[i],alldinamZ)
      vadiiZ = values(dilistZ[indiiZ])
      nvadiiZ = nelem(vadiiZ)
      if i==1
        txt1=""
      else
        txt1=","
      endif
      dilevsZ = dilevsZ~txt1~diliZ
      dialllZ[i] = nvadiiZ
    endloop
  endif
  #print dialllZ
  ##################
  if (divars == "")&&(divarsZ == "")
    matr3 = {}
  else
    uniqZdivars = str_setdiff(divarsZ,divars)
    diZind = findinds(uniqZdivars,divarsZ)
    NOUNIQ = nelem(strsplit(divarsZ)) - nelem(strsplit(uniqZdivars))
    ALLDIVARS = strstrip(divars~" "~uniqZdivars)
    # eval "ALLDIVARS"
    # print ALLDIVARS
    if strlen(dilevs)==0
      ALLDILEVS = dilevsZ
      # eval "ALLDILEVS"
      # print ALLDILEVS
    elif NOUNIQ ==0
      if strlen(dilevsZ)==0
        ALLDILEVS = dilevs
      else
        ALLDILEVS = dilevs~","~dilevsZ
      endif
      # eval "ALLDILEVS"
      # print ALLDILEVS
    else
      #print diZind
      if nelem(diZind) == 0
        ALLDILEVS = dilevs
      else
        dilevsZs = strsplit3(dilevsZ,",")
        ndZ = nelem(dilevsZs)
        ALLDILEVS = dilevs
        loop i=1..ndZ -q
          #eval iminc(abs(diZind - i))
          if minc(abs(diZind - i))==0
            #strizi = dilevsZs[i]
            #strizis[1] = strizi
            #DIZ = DIZ+ strizis
            ALLDILEVS = ALLDILEVS~","~dilevsZs[i]
          endif
        endloop
        # eval "ALLDILEVS"
        # print ALLDILEVS
      endif
    endif
    if nelem(diZind)==0
      matrix DIALL = dialll
      # eval "DIALL"
      # print DIALL
    else
      matrix DIALL = dialll~dialllZ[diZind]
      # eval "DIALL"
      # print DIALL
    endif
    cst1 = inlist(b.xlist,const) > 0
    cst2z = inlist(b.zlist,const) > 0
    uniqznam = str_setdiff(allnamesZ,allnames)
    zZindi = findinds(uniqznam,allnamesZ)
    ALLNAMES = strstrip(allnames~" "~uniqznam)
    #eval "ALLNAMES"
    #print ALLNAMES
    matrix ALLVALUES = allvalues~allvaluesZ[zZindi]
    #print ALLVALUES
    matr3 = {}
    divarss = strsplit(ALLDIVARS)
    ndss = nelem(divarss)
    dilevelss = strsplit3(ALLDILEVS,",")
    loop i=1..ndss -q
      alllevelsi = mk_alllevs(DIALL[i])
      matr3i = hk_prefun(newformula,newformulaZ,ALLNAMES,divarss[i], dilevelss[i],alllevelsi,ALLVALUES,b.base_coeffs,b.ncoeff,cst1,cst2z,coe,link)
      ####
      matr3 = matr3|matr3i
      ####
    endloop
    divars = ALLDIVARS
    dilevels = ALLDILEVS
    dialll = DIALL
  endif
  #############3
  #         divars = ALLDIVARS
  #         dilevels = ALLDILEVS
  #         dialll = DIALL
  ############3
endif
if cmd == "tobit"
  vcov = b.full_vcv
  s__ = b.sigma
  catch lo__ = b.llimit
  errtob1 = $error
  if errtob1 > 0
    lo__ = 0
  endif
  catch hi__ = b.rlimit
  errtob2 = $error
  if errtob2 > 0
    hi__ = $huge
    hi22 = NA
  else
    hi22 = hi__
  endif
  if link != "inside"
    if lang == "uk_UA"
      printf "\n%s: граничні ефекти на E(%s)\n", cmd, name_of_y
    else
      printf "\n%s marginal effects on E(%s)\n", cmd, name_of_y
    endif
  else
    tytyt = mk_tob_outtext(lo__,hi22,name_of_y)
    if lang == "uk_UA"
      printf "\n%s: граничні ефекти на E(%s|%s)\n", cmd,name_of_y, tytyt
    else
      printf "\n%s marginal effects on E(%s|%s)\n", cmd,name_of_y, tytyt
    endif
  endif
  if  contvars == ""
    matr1 = {}
  else
    if link != "inside"
      matr1 = tb_cont_grads(newformula,allnames,allvalues,contvars,coe,s__,lo__,hi__)
    else
      matr1 = tb_cont_grads22(newformula,allnames,allvalues,contvars,coe,s__,lo__,hi__)
    endif
  endif
  if  dumvars == ""
    matr2 = {}
  else
    if link != "inside"
      matr2 = tb_dum_grads(newformula,allnames,allvalues,dumvars,coe,s__,lo__,hi__)
    else
      matr2 = tb_dum_grads22(newformula,allnames,allvalues,dumvars,coe,s__,lo__,hi__)
    endif
  endif
  if  divars == ""
    matr3 = {}
  else
    matr3 = {}
    divarss = strsplit(divars)
    ndss = nelem(divarss)
    dilevelss = strsplit3(dilevs,",")
    loop i=1..ndss -q
      if link != "inside"
        matr3i = t_discr_grads2(newformula, allnames, allvalues, divarss[i],dilevelss[i],mk_alllevs(dialll[i]),coe,s__,lo__,hi__)
      else
        matr3i = t_discr_grads222(newformula, allnames, allvalues, divarss[i],dilevelss[i],mk_alllevs(dialll[i]),coe,s__,lo__,hi__)
      endif
      matr3 = matr3|matr3i
    endloop
  endif
endif
if oneq||(cmd == "heckit")
  matrix matr=matr1|matr2|matr3
  matrix mfxs = matr[,1]
  matr = matr[,2:]
  matrix mvc = qform(matr,vcov)
  matrix coeffmat  = mfxs~sqrt(diag(mvc))
elif cmd=="tobit"
  matrix matr=matr1|matr2|matr3
  matrix mfxs = matr[,1]
  matr = matr[,2:]
  nnn = cols(vcov)
  matrix Jm = I(nnn)
  Jm[nnn,nnn] = s__
  matrix Vbig = qform(Jm, vcov)
  matrix mvc = qform(matr,Vbig)
  matrix coeffmat  = mfxs~sqrt(diag(mvc))
endif
if oneq||cmd=="tobit"||cmd=="heckit"
  lang = $lang
  nc = 0
  nd = 0
  if contvars != ""
    nc = nelem(strsplit(contvars))
    matrix coeffmat1 = coeffmat[1:nc,]
    if lang == "uk_UA"
      if nc == 1
        flex = ""
        flex2 = "ої"
      else
        flex = "и"
        flex2 = "их"
      endif
      printf "\nЕфект%s неперервн%s змінн%s: \n",flex,flex2,flex2
      modprint coeffmat1 contvars
    else
      if nc == 1
        flex = ""
      else
        flex = "s"
      endif
      printf "\nEffect%s of continuous variable%s: \n",flex,flex
      modprint coeffmat1 contvars
    endif
  endif
  if dumvars != ""
    nd = nelem(strsplit(dumvars))
    matrix coeffmat2 = coeffmat[(nc+1):(nc+nd),]
    if lang == "uk_UA"
      if nc == 1
        flex = ""
        flex2 = "ої"
      else
        flex = "и"
        flex2 = "их"
      endif
      printf "\nЕфект%s фіктивн%s змінн%s: \n",flex,flex2,flex2
      modprint coeffmat2 dumvars
    else
      if nd == 1
        flex = ""
      else
        flex = "s"
      endif
      printf "\nEffect%s of dummy variable%s: \n",flex,flex
      modprint coeffmat2 dumvars
    endif
  endif
  bigdivars = ""
  if divars != ""
    ndi = nelem(divarss)
    if lang == "uk_UA"
      if ndi == 1
        flex2 = "ої"
      else
        flex2 = "их"
      endif
      printf "\nЕфекти дискретн%s змінн%s: \n",flex2,flex2
    else
      if ndi == 1
        flex = ""
      else
        flex = "s"
      endif
      printf "\nEffects of discrete variable%s.\n",flex
    endif
    counter = nc+nd+1
    matrix matri={}
    loop  i=1..ndi -q
      levi = dilevelss[i]
      levis = strsplit(levi)
      ndii = nelem(levis)
      nami = divarss[i]
      isallev = (ndii==dialll[i])
      if isallev
        basi = "1"
        ndii -=1
        levi = str_setdiff(levi,basi)
      else
        alli = mk_alllevs(dialll[i])
        basi = str_setdiff(alli,levi)
      endif
      basiold = basi
      if (dialll[i] - ndii) > 1
        basis = strsplit(basi)
        basi = basis[1]
      endif
      if lang == "uk_UA"
        printf "\n Ефекти %s (базовий рівень = %s):\n",nami,basi
      else
        printf "\n Effects of %s (base level = %s):\n",nami,basi
      endif
      matri = coeffmat[counter:(counter+ndii-1),]
      namii = mk_nami(nami, levi)
      modprint matri namii
      if (dialll[i] - ndii) > 1
        babasi = strsub(basiold," ",", ")
        if lang == "uk_UA"
          printf "  Увага: поточна модель не розрізняє такі рівні %s: %s.\n",nami, babasi
        else
          printf "  Note: levels %s of %s are indistinguishable in the current model.\n",babasi,nami
        endif
      endif
      counter = counter+ndii
      bigdivars = bigdivars~" "~mk_namlevi(nami,levi)
    endloop
  endif
  bigdivars = strstrip(bigdivars)
  newnames = strstrip(contvars~" "~dumvars~" "~bigdivars)
  bundle bout
  bout["names"] = newnames
  matrix bout["mfx"] = coeffmat[,1]
  matrix bout["se"] = coeffmat[,2]
  matrix bout["vcov"] = mvc
endif
isordered = (cmd2 == "orderedL")||(cmd2 == "orderedP")
if isordered
  list yy = b.ylist
  printf "\nMarginal effects in the ordered %s model.\n", cmd
  nalts = nelem(uniq(yy[1]))
  if  contvars == ""
    matr1 = {}
  else
    matr1 = olp_cont_grads2(newformula,allnames,allvalues,contvars,coe,nalts,cmd)
  endif
  if  dumvars == ""
    matr2 = {}
  else
    matr2 = olp_dum_grads2(newformula,allnames,allvalues,dumvars,coe,nalts,cmd)
  endif
  if  divars == ""
    matr3 = {}
  else
    matr3 = {}
    divarss = strsplit(divars)
    ndss = nelem(divarss)
    dilevelss = strsplit3(dilevs,",")
    loop i=1..ndss -q
      matr3i = olp_discr_grad2(newformula, allnames,allvalues, divarss[i],dilevelss[i],mk_alllevs(dialll[i]),coe,nalts,cmd)
      matr3 = matr3|matr3i
    endloop
  endif
  matrix matr=matr1|matr2|matr3
  rm = rows(matr)
  ra = rm/nalts
  sequ =  seq(0,ra-1)
  bigseq = zeros(rm,1)
  loop i=1..nalts -q
    seqi = i.+ nalts.*sequ
    indsi = ra*(i-1).+sequ.+1
    bigseq[indsi]=seqi'
  endloop
  matr = matr[bigseq,]
  matrix mfxs = matr[,1]
  matr = matr[,2:]
  matrix mvc = qform(matr,vcov)
  matrix uuu = sqrt(diag(mvc))
  matrix bigcoeffmat  = mfxs~sqrt(diag(mvc))
  lang = $lang
  nc = 0
  nd = 0
  matrix coeffmat ={}
  string nam_y = varname(yy)
  loop iii=1..nalts -q
    indsi = ra*(iii-1).+sequ.+1
    coeffmat =  bigcoeffmat[indsi,]
    if lang == "uk_UA"
      printf "\nГраничні ефекти для Pr(%s = %d)\n",nam_y,iii
    else
      printf "\nMarginal effects on Pr(%s = %d)\n",nam_y,iii
    endif
    if contvars != ""
      nc = nelem(strsplit(contvars))
      matrix coeffmat1 = coeffmat[1:nc,]
      if lang == "uk_UA"
        if nc == 1
          flex = ""
          flex2 = "ої"
        else
          flex = "и"
          flex2 = "их"
        endif
        printf "\nЕфект%s неперервн%s змінн%s: \n",flex,flex2,flex2
        modprint coeffmat1 contvars
      else
        if nc == 1
          flex = ""
        else
          flex = "s"
        endif
        printf "\nEffect%s of continuous variable%s: \n",flex,flex
        modprint coeffmat1 contvars
      endif
    endif
    if dumvars != ""
      nd = nelem(strsplit(dumvars))
      matrix coeffmat2 = coeffmat[(nc+1):(nc+nd),]
      if lang == "uk_UA"
        if nc == 1
          flex = ""
          flex2 = "ої"
        else
          flex = "и"
          flex2 = "их"
        endif
        printf "\nЕфект%s фіктивн%s змінн%s: \n",flex,flex2,flex2
        modprint coeffmat2 dumvars
      else
        if nd == 1
          flex = ""
        else
          flex = "s"
        endif
        printf "\nEffect%s of dummy variable%s: \n",flex,flex
        modprint coeffmat2 dumvars
      endif
    endif
    bigdivars = ""
    if divars != ""
      ndi = nelem(divarss)
      if lang == "uk_UA"
        if ndi == 1
          flex2 = "ої"
        else
          flex2 = "их"
        endif
        printf "\nЕфекти дискретн%s змінн%s: \n",flex2,flex2
      else
        if ndi == 1
          flex = ""
        else
          flex = "s"
        endif
        printf "\nEffects of discrete variable%s.\n",flex
      endif
      counter = nc+nd+1
      matrix matri={}
      loop  i=1..ndi -q
        levi = dilevelss[i]
        levis = strsplit(levi)
        ndii = nelem(levis)
        nami = divarss[i]
        isallev = (ndii==dialll[i])
        if isallev
          basi = "1"
          ndii -=1
          levi = str_setdiff(levi,basi)
        else
          alli = mk_alllevs(dialll[i])
          basi = str_setdiff(alli,levi)
        endif
        basiold = basi
        if (dialll[i] - ndii) > 1
          basis = strsplit(basi)
          basi = basis[1]
        endif
        if lang == "uk_UA"
          printf "\n Ефекти %s (базовий рівень = %s):\n",nami,basi
        else
          printf "\n Effects of %s (base level = %s):\n",nami,basi
        endif
        matri = coeffmat[counter:(counter+ndii-1),]
        namii = mk_nami(nami, levi)
        modprint matri namii
        if (dialll[i] - ndii) > 1
          babasi = strsub(basiold," ",", ")
          if lang == "uk_UA"
            printf "  Увага: поточна модель не розрізняє такі рівні %s: %s.\n",nami, babasi
          else
            printf "  Note: levels %s of %s are indistinguishable in the current model.\n",babasi,nami
          endif
        endif
        counter = counter+ndii
        bigdivars = bigdivars~" "~mk_namlevi(nami,levi)
      endloop
    endif
  endloop
  bigdivars = strstrip(bigdivars)
  newnames = strstrip(contvars~" "~dumvars~" "~bigdivars)
  bundle bout
  bout["names"] = newnames
  matrix bout["mfx"] = bigcoeffmat[,1]
  matrix bout["se"] = bigcoeffmat[,2]
  matrix bout["vcov"] = mvc
  bout["n_alts"] = nalts
  bout["nam_y"] = nam_y
  #print bout
endif
if cmd2 == "multinomial"
  list yy = b.ylist
  nnyy = values(yy[1])
  nalts = nelem(nnyy)
  ismucst = inlist(b.xlist,const)
  #print ismucst
  printf "\nMarginal effects after multinomial logit.\n"
  if  contvars == ""
    matr10 = {}
    n_contvars = 0
  else
    matr10 = ml_cont_grads(coe,newformula,allnames,allvalues,contvars,nalts,ismucst)
    n_contvars = rows(matr10)
    #print matr10
  endif
  if  dumvars == ""
    matr20 = {}
    n_dumvars = 0
  else
    matr20 = ml_dum_grads(newformula,allnames,allvalues,dumvars,coe,nalts,ismucst)
    n_dumvars = rows(matr20)
  endif
  if  divars == ""
    matr300 = {}
    bign_divars = 0
    didims = {}
  else
    matr300 = {}
    divarss = strsplit(divars)
    ndss = nelem(divarss)
    matrix didims = zeros(1,ndss)
    dilevelss = strsplit3(dilevs,",")
    loop i=1..ndss -q
      matr300i = ml_discr_grads2(newformula, allnames,allvalues, divarss[i],dilevelss[i],mk_alllevs(dialll[i]),coe,nalts,ismucst)
      didims[i] = rows(matr300i)
      matr300 =  matr300|matr300i
    endloop
    bign_divars = rows(matr300)
    matr300 = matr300[didindi(didims,nalts),]
  endif
  nb = nelem(coe)
  mega_n = n_contvars+n_dumvars+bign_divars
  large_n = mega_n/nalts
  matrix mega_matr = zeros(mega_n,nb+1)
  if n_contvars > 0
    count_seq = seq(1,n_contvars/nalts)
  endif
  if n_dumvars > 0
    dum_seq2 = seq(1,n_dumvars/nalts)
    dum_seq = dum_seq2 + n_contvars/nalts
  endif
  if bign_divars > 0
    di_seq2 = seq(1,bign_divars/nalts)
    di_seq = di_seq2 + (n_contvars + n_dumvars)/nalts
  endif
  loop i = 1..nalts -q
    i_counter = large_n * (i - 1)
    if n_contvars > 0
      mega_matr[i_counter+count_seq,] = matr10 [n_contvars/nalts* (i - 1)+count_seq,]
    endif
    if n_dumvars > 0
      mega_matr[i_counter+dum_seq,] = matr20 [n_dumvars/nalts* (i - 1)+dum_seq2,]
    endif
    if bign_divars > 0
      mega_matr[i_counter+di_seq,] = matr300 [bign_divars/nalts* (i - 1)+di_seq2,]
    endif
  endloop
  matrix all_mfx = mega_matr[,1]
  mega_matr = qform(mega_matr[,2:],vcov)
  matrix all_coeffmat  = all_mfx ~ sqrt(diag(mega_matr))
  #print all_coeffmat
  lang = $lang
  if contvars != ""
    nc = n_contvars/nalts
    if lang == "uk_UA"
      if nc == 1
        flex = ""
        flex2 = "ої"
      else
        flex = "и"
        flex2 = "их"
      endif
      cont_text1 = sprintf( "\nЕфект%s неперервн%s змінн%s: \n",flex,flex2,flex2)
    else
      if nc == 1
        flex = ""
      else
        flex = "s"
      endif
      cont_text1 = sprintf( "\nEffect%s of continuous variable%s: \n",flex,flex)
    endif
  else
    cont_text1 = ""
    nc = 0
  endif
  if dumvars != ""
    nd = n_dumvars/nalts
    if lang == "uk_UA"
      if nc == 1
        flex = ""
        flex2 = "ої"
      else
        flex = "и"
        flex2 = "их"
      endif
      dum_text1 = sprintf("\nЕфект%s фіктивн%s змінн%s: \n",flex,flex2,flex2)
    else
      if nd == 1
        flex = ""
      else
        flex = "s"
      endif
      dum_text1 = sprintf( "\nEffect%s of dummy variable%s: \n",flex,flex)
    endif
  else
    dum_text1 = ""
    nd = 0
  endif
  bigdivars = ""
  if divars != ""
    if lang == "uk_UA"
      if ndis == 1
        flex2 = "ої"
      else
        flex2 = "их"
      endif
      di_text1 = sprintf( "\nЕфекти дискретн%s змінн%s: \n",flex2,flex2)
    else
      if ndis == 1
        flex = ""
      else
        flex = "s"
      endif
      di_text1 = sprintf( "\nEffects of discrete variable%s.\n",flex)
    endif
    strings basiss = array(ndis)
    strings babasiss = array(ndis)
    strings namiiss = array(ndis)
    loop  i=1..ndis -q
      levi = dilevelss[i]
      ndi = didims[i]/nalts
      nami = divarss[i]
      isallev = (ndi==dialll[i])
      if isallev
        basi = "1"
      else
        alli = mk_alllevs(dialll[i])
        basi = str_setdiff(alli,levi)
      endif
      basiold = basi
      if (dialll[i] - ndi) > 1
        basis = strsplit(basi)
        basi = basis[1]
      endif
      basiss[i] = basi
      namii = mk_nami(nami, levi)
      namiiss[i] = namii
      if (dialll[i] - ndi) > 1
        babasiss[i] = strsub(basiold," ",", ")
      else
        babasiss[i] = "_empty"
      endif
      bigdivars = bigdivars~" "~mk_namlevi(nami,levi)
    endloop
  endif
  bigdivars = strstrip(bigdivars)
  newnames = strstrip(contvars~" "~dumvars~" "~bigdivars)
  loop big_i = 1..nalts -q
    printf "\nEffects on P(%s = %d):\n", name_of_y, big_i
    matrix coeffmat = all_coeffmat[(large_n * (big_i -1) + 1): large_n * big_i,]
    if contvars != ""
      matrix coeffmat1 = coeffmat[1:nc,]
      printf "%s", cont_text1
      modprint coeffmat1 contvars
    endif
    if dumvars != ""
      matrix coeffmat2 = coeffmat[(nc+1):(nc+nd),]
      printf "%s", dum_text1
      modprint coeffmat2 dumvars
    endif
    if divars != ""
      printf "%s", di_text1
      counter = nc+nd+1
      loop  i=1..ndis -q
        nami = divarss[i]
        basi = basiss[i]
        if lang == "uk_UA"
          printf "\n Ефекти %s (базовий рівень = %s):\n",nami,basi
        else
          printf "\n Effects of %s (base level = %s):\n",nami,basi
        endif
        matrix matri = coeffmat[counter:(counter+didims[i]/nalts-1),]
        namii = namiiss[i]
        modprint matri namii
        counter = counter+didims[i]/nalts
      endloop
    endif
  endloop
  bigdivars = strstrip(bigdivars)
  newnames = strstrip(contvars~" "~dumvars~" "~bigdivars)
  bundle bout
  bout["names"] = newnames
  matrix bout["mfx"] = all_coeffmat[,1]
  matrix bout["se"] = all_coeffmat[,2]
  matrix bout["vcov"] = mega_matr
  string bout["pkg"] = "a_eff"
  string bout["command"] = b.command
  string bout["cmd2"] = cmd2
  scalar bout["nalt"] = 1
  #         if isordered
  #             bout["nalt"] = nalts
  #         endif
  #print bout
endif
string bout["pkg"] = "a_eff"
string bout["command"] = b.command
string bout["cmd2"] = cmd2
scalar bout["nalt"] = 1
if (cmd2 == "orderedL")||(cmd2 == "orderedP")||(cmd2 == "multinomial")
  bout["nalt"] = nalts
endif
bout["type"] = "mfx"
return bout

output of mfx() z- or t- based intervals conf. level, NA = 0.95

iscl = missing(cl)
if iscl
  cl = 0.95
endif
matrix mafx = b["mfx"]
matrix ses = b["se"]
names = b["names"]
if df==0
  stat = qnorm((1+cl)/2)
else
  stat = critical(t, df, (1-cl)/2)
endif
matrix half = ses.*stat
matrix left = mafx - half
matrix right = mafx + half
matrix ret = mafx~ses~left~right
rnameset(ret,strsplit(names))
proc = round(100*cl)
textproc = sprintf("%d", proc)
leb = textproc~"% l.b."
rib = textproc~"% r.b."
if b.type == "el"
  cna = "dy/dx*x/y,s.e.,"~leb~","~rib
else
  cna = "dy/dx,s.e.,"~leb~","~rib
endif
cnas = strsplit3(cna,",")
cnameset(ret,cnas)
return ret

list of vars with 2-nd order terms list of vars 'as is' whether to include cubes

li1 = li1 - 0
li2 = li2 - 0
list ret = null
namli1 = varnames(li1)
namli110 = cut_names(namli1)
namli2 = varnames(li2)
namli220 = cut_names(namli2)
bign = nelem(namli110)
n2 = nelem(namli220)
list ret = null
string descr
strings S = array(2)
scalar count = 0
loop i=1..bign -q
  loop j = i..bign -q
    loop l = 1..n2 -q
      name = namli110[i]~"_"~namli110[j]~"_"~namli220[l]
      descr = namli1[i]~"*"~namli1[j]~"*"~namli2[l]
      S = corr_d_n_cb2(namli110[i],namli110[j],namli220[l],descr)
      name = S[1]
      descr = S[2]
      if cubes == 0
        uni = uniqvarlist(descr)
        unis = strsplit(uni)
        nnn = nelem(unis)
        if nnn >1
          ret += genseries(name,li1[i]*li1[j]*li2[l])
          text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
          @text2
        endif
      else
        ret += genseries(name,li1[i]*li1[j]*li2[l])
        text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
        @text2
      endif
    endloop
  endloop
endloop
return dropcoll(ret)

$xlist in most cases forced ind. vars

nvar = $nvars
stri = mk_txt_range(nvar - 1)
list da = @stri
vns = varname(da)
allvns = strsub(varname(da),","," ")
xlist = xlist - 0
nx = nelem(xlist)
namss = varnames(xlist)
ret_fla = ""
ret_var = ""
ret_dis = ""
loop i = 1..nx -q
  nam = namss[i]
  list Normals = xlist[i]
  reti = smart_formula200(Normals,indlist,stri,da,allvns)
  retis = strsplit3(reti,";")
  forli = retis[1]
  varsii = retis[2]
  discri = retis[3]
  if retis[2]!="_empty"&&(!inlist(indlist,xlist[i]))
    varsi = strsplit(retis[2])
    nxi = nelem(varsi)
    string normalsi = ""
    loop j=1..nxi -q
      varsij = varsi[j]
      ckj = inlist(da,varsij)
      if ckj
        normalsij = sprintf(" %d",ckj)
        normalsi = normalsi~normalsij
      endif
    endloop
    list dai = @normalsi
    oldvars = varnames(dai)
    retii = smart_formula200(dai,indlist,stri,da,allvns)
    retiis = strsplit3(retii,";")
    newvars = strsplit(retiis[1])
    varsii = retiis[2]
    counter = nelem(oldvars)
    loop j=1..counter -q
      if  oldvars[j]!=newvars[j]
        forli = smartstrsub(forli,oldvars[j],newvars[j])
      endif
    endloop
  endif
  ret_fla = ret_fla~" "~forli
  if varsii != "_empty"
    ret_var =  ret_var~" "~varsii
  endif
  if discri != "_empty"
    ret_dis =  ret_dis ~" "~discri
  endif
endloop
ret_fla = strstrip(ret_fla)
ret_var = strstrip(ret_var)
ret_var = uniqvarlist(ret_var)
ret_dis = uniqvarlist(ret_dis)
ret_dis = strstrip(ret_dis)
bigret = ret_fla~";"~ret_var~";"~ret_dis
return bigret

list of series whether to include cubes

li1 = li1 - 0
list ret = null
namli1 = varnames(li1)
namli110 = cut_names(namli1)
bign = nelem(namli110)
list ret = null
string descr
strings S = array(2)
scalar count = 0
loop i=1..bign -q
  loop j = i..bign -q
    loop l = i..bign -q
      name = namli110[i]~"_"~namli110[j]~"_"~namli110[l]
      descr = namli1[i]~"*"~namli1[j]~"*"~namli1[l]
      S = corr_d_n_cb2(namli110[i],namli110[j],namli110[l],descr)
      name = S[1]
      descr = S[2]
      if cubes == 0
        if (i!=j)||(i!=l)
          ret += genseries(name,li1[i]*li1[j]*li1[l])
          text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
          @text2
        endif
      else
        ret += genseries(name,li1[i]*li1[j]*li1[l])
        text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
        @text2
      endif
    endloop
  endloop
endloop
return dropcoll(ret)

$xlist in most cases forced ind. vars

nvar = $nvars
stri = mk_txt_range(nvar - 1)
list da = @stri
allvns = strsub(varname(da),","," ")
# print allvns
list x = x-0
bignamss = varnames(x)
hugenams = strsub(varname(da),","," ")
#print hugenams
cknull = nelem(indlist)==0
if !cknull
  indnams = strjoin(varnames(indlist))
else
  indnams = ""
endif
string line
dili = ""
coli = ""
dilis = ""
colis = ""
string first
rets = ""
rets_left = ""
rets_right = ""
loop foreach i x -q
  cknull = 0
  bundle bui = getinfo(x[i])
  if !cknull
    fififi = inlist(indlist,x[i])
    if fififi == 0
      fififi = !inbundle(bui,"description")
      #print fififi
    endif
    if fififi
      bignamsi = bignamss[i]
      rets = strstrip(rets~" "~bignamsi)
      colis = strstrip(colis~" "~bignamsi)
    else
      cknull = 1
    endif
  endif
  if cknull
    #bundle bui = getinfo(x[i])
    vvv = bui.name
    uuu = bui.description
    ul = uniqvarlist(uuu)
    isspace = isstrstr(uuu," ")
    if isspace == 0
      isfrla = isfla(uuu)
      isfrla2 = vf(uuu)
      #print isfrla2
      if isfrla&&(isfrla2 > 0)
        ret3 = uuu
        flavars = uniqvarlist(uuu)
      else
        ret3 = vvv
        flavars = vvv
      endif
      nfla = nelem(strsplit(flavars))
      #eval "ret3"
      #print ret3
      #eval "flavars"
      #print flavars
      flavarss = strsplit(flavars)
      bubu_c = ""
      bubu_di = ""
      funij = ret3
      loop ij = 1..nfla -q
        disij =  disuspect(flavarss[ij])
        #eval disij
        #print disij
        if disij != "_empty"
          disijss = strsplit(disij)
          preff = disijss[1]
          nsuff = atof(disijss[2])
          findij = findind(preff,hugenams)
          #print nsuff
          #print findij
          ckij = findij
          if findij > 0
            unqij = values(da[findij])'
            if unqij != seq(1,nelem(unqij))
              ckij = 0
            endif
          endif
          preff2 = substr(preff,2,strlen(preff))
          findij2 = findind(preff2,hugenams)
          #print findij2
          ckij2 = findij2
          if findij2 > 0
            unqij2 = values(da[findij2])'
            #print unqij2
            if unqij2 != seq(1,nelem(unqij2))
              ckij2 = 0
            endif
          endif
          #print ckij ckij2
          if (ckij+ckij2) == 0
            bubu_c = bubu_c~" "~flavarss[ij]
          elif (ckij > 0)&&(ckij2 == 0)
            bubu_di = bubu_di~" "~preff
          elif (ckij == 0)&&(ckij2 > 0)
            bubu_di = bubu_di~" "~preff2
            substij = strstrip(preff2~"_"~disijss[2])
            funij = smartstrsub(funij,flavarss[ij],substij)
          else
            printf "\nAmbigous var. names: both %s and %s are discrete",pref,pref2
            printf "\nYou may need manual correction of formula"
            printf "\nor renaming variables\n"
            bubu_di = bubu_di~" "~preff
          endif
        else
          bubu_c = bubu_c~" "~flavarss[ij]
        endif
      endloop
      bubu_c = strstrip(bubu_c)
      bubu_di = strstrip(bubu_di)
      ret3 = funij
      lili = bubu_c
      didi = bubu_di
      colis = colis~" "~lili
      dilis = dilis~" "~didi
      rets = rets~" "~ret3
    else
      iseq = substr(uuu,1,1)=="="
      if iseq
        issq = substr(vvv,1,3)=="sq_"
        sdpart = substr(vvv,4,strlen(vvv))
        fifi = findind(sdpart,ul)
        issq = issq*fifi>0
        islg = substr(vvv,1,2)=="l_"
        sdpart2 = substr(vvv,3,strlen(vvv))
        fifi2 = findind(sdpart2,ul)
        fifi20l = str_setinrersect(ul,allvns)
        new_condition = sdpart2==fifi20l
        islg = islg*fifi2*new_condition>0
        underscpos = catch_times(vvv, uuu)
        if issq
          ret3 = strstrip(sdpart~"^2")
          flavars = uniqvarlist(ret3)
          lili = flavars
          colis = colis~" "~lili
          rets = rets~" "~ret3
        elif islg
          ret3 = strstrip("log("~sdpart2~")")
          flavars = uniqvarlist(ret3)
          lili = flavars
          colis = colis~" "~lili
          rets = rets~" "~ret3
        elif  underscpos > 0
          nvvv = strlen(vvv)
          part1 = substr(vvv,1,underscpos-1)
          part2 = substr(vvv,underscpos+1,nvvv)
          ret3 = part1~"*"~part2
          rets = rets~" "~ret3
        else
          ret3 = vvv
          rets = rets~" "~ret3
        endif
        flavars = uniqvarlist(ret3)
        lili = flavars
        colis = colis~" "~lili
      else
        inti = substr(uuu,1,14)=="interaction of"
        inti2 = isstrstr(uuu,"==")
        intis = strsplit3(uuu,"==")
        nintis = nelem(intis)
        inti3 = (nintis==2)
        inti4 = 0
        if inti3
          intis2 = intis[2]
          nintis2 = strlen(intis2)
          cknin1 = substr(intis2,nintis2,nintis2)==")"
          cknin2 = isnat(strstrip(substr(intis2,1,nintis2-1)))
          inti4 = cknin1*cknin2
        endif
        inti = inti*inti2*inti3*inti4>0
        if inti
          zx = strsplit3(uuu,"and (")
          zx2 = zx[2]
          zx1 = zx[1]
          namis2 = strsplit3(zx2,"==")
          nami2 = strstrip(namis2[1])
          ni1 = strlen(zx1)
          nami1 = strstrip(substr(zx1,15,ni1))
          didi = getdiname(nami1)
          if strlen(didi) == 0
            first = nami1
          else
            first = null
          endif
          lili = strstrip(didi~" "~nami2)
          ret3 = doret2(uuu,vvv,lili,first)
          if !exists(first)
            first= ""
          endif
          rets = rets~" "~ret3
          dilis = dilis~" "~lili
          colis = colis~" "~first
        else
          doommy = isduname(vvv)
          #eval "doommy"
          #print doommy
          disij =  disuspect(vvv)
          #print disij
          bubu_di = ""
          bubu_c = ""
          funij = vvv
          ij = 1
          flavars = vvv
          flavarss = strsplit(flavars)
          if disij != "_empty"
            disijss = strsplit(disij)
            preff = disijss[1]
            nsuff = atof(disijss[2])
            findij = findind(preff,hugenams)
            #print nsuff
            #print findij
            ckij = findij
            if findij > 0
              unqij = values(da[findij])'
              if unqij != seq(1,nelem(unqij))
                ckij = 0
              endif
            endif
            preff2 = substr(preff,2,strlen(preff))
            findij2 = findind(preff2,hugenams)
            #print findij2
            ckij2 = findij2
            if findij2 > 0
              unqij2 = values(da[findij2])'
              #print unqij2
              if unqij2 != seq(1,nelem(unqij2))
                ckij2 = 0
              endif
            endif
            # eval "ckij ckij2"
            #print ckij ckij2
            if (ckij+ckij2) == 0
              bubu_c = bubu_c~" "~flavarss[ij]
            elif (ckij > 0)&&(ckij2 == 0)
              bubu_di = bubu_di~" "~preff
            elif (ckij == 0)&&(ckij2 > 0)
              bubu_di = bubu_di~" "~preff2
              substij = strstrip(preff2~"_"~disijss[2])
              funij = smartstrsub(funij,flavarss[ij],substij)
            else
              printf "\nAmbigous var. names: both %s and %s are discrete",pref,pref2
              printf "\nYou may need manual correction of formula"
              printf "\nor renaming variables\n"
              bubu_di = bubu_di~" "~preff
            endif
            #eval "bubus"
            #print bubu_c
            #eval "bubudi"
            #print bubu_di
            #eval "funij"
            # print funij
            ret3 = funij
            didi = bubu_di
            dilis = dilis~" "~didi
          else
            bubu_c = bubu_c~" "~flavarss[ij]
            ret3 = vvv
            lili = ret3
            colis = colis~" "~lili
          endif
          rets = rets~" "~ret3
        endif
      endif
    endif
    colis = strstrip(colis)
    dilis = strstrip(dilis)
    rets_left = strstrip(rets_left~" "~vvv)
    rets_right = rets_right~" "~ret3
  endif
endloop
rets = strstrip(rets)
discrsi = uniqvarlist(dilis)
conti = uniqvarlist(colis)
return strstrip(rets)~";"~ strstrip(conti)~";"~strstrip(discrsi)

list of series whether to include squares

li1 = li1 - 0
list ret = null
namli1 = varnames(li1)
namli110 = cut_names(namli1)
bign = nelem(namli110)
list ret = null
string descr
strings S = array(2)
scalar count = 0
loop i=1..bign -q
  if squares==0
    count = i+1
  else
    count = i
  endif
  loop j = count..bign -q
    name = namli110[i]~"_"~namli110[j]
    descr = namli1[i]~"*"~namli1[j]
    uni = uniqvarlist(descr)
    unis = strsplit(uni)
    n = nelem(unis)
    if n == 1
      S = corr_d_n_sq(name,descr)
      name = S[1]
      descr = S[2]
    endif
    ret += genseries(name,li1[i]*li1[j])
    text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
    @text2
  endloop
endloop
return dropcoll(ret)

$xlist in most cases indep. cont. and dummy vars indep discrete vars

gretl_version = $version
xlist = xlist-const
n = nelem(xlist)
string line
vns = varnames(xlist)
nv = nelem(vlist)
ndi = nelem(discrlist)
if (nv + ndi) == 0
  funcerr "no independent variables given"
endif
indlist = vlist || discrlist
alln = strjoin(varnames(indlist))
if nv == 0
  string othn = null
  isno = !exists(othn)
else
  string othn = strjoin(varnames(vlist))
endif
if ndi == 0
  string disn = null
else
  string disn = strjoin(varnames(discrlist))
endif
ret = ""
loop i=1..n -q
  inli = inlist(indlist,xlist[i])
  if inli
    ret3 = vns[i]
    ret = ret~" "~ret3
  else
    outfile mybuf --buffer --quiet
      labels xlist[i]
      outfile --close
        string s = mybuf
        scalar j = 1
        loop while getline(s, line)&&(j<3) -q
          uu=sprintf( "%s\n",  line)
          j++
        endloop
        bign = strlen(uu)
        uuu = strstr(uu,":")
        n = strlen(uuu)
        vvv = strstrip(substr(uu,1,bign-n))
        uuu = substr(uuu,3,n)
        uuu = strstrip(uuu)
        ret3 = doret2(uuu,vvv,disn,othn)
        ret = ret~" "~ret3
      endif
    endloop
    ret = strstrip(ret)
    if nv == 0
      retvl = "_empty"
    else
      retvl = strjoin(varnames(vlist))
    endif
    if ndi == 0
      retdis = "_empty"
    else
      retdis = strjoin(varnames(discrlist))
    endif
    ret = ret~";"~retvl~";"~retdis
    return ret

b is $model in most cases the model formula at and wrt vars names, eff type values of 'at vars' dataset, for gretl before 2018d

gretl_version = $version
cmd = b.command
name_of_y = b.depvar
matrix coe = b.coeff
matrix vcov = b.vcv
#list xx = b.xlist
#funcerr "just stopped"
exifla = exists(formula)
if !exifla
  formula = smart_formula2(b.xlist)
  ##
else
  isdotcomma = instring(formula,";")
  #print isdotcomma
  if !isdotcomma
    formula_backup = formula
    formula_2s = strsplit(formula)
    nf2s = nelem(formula_2s)
    list indep0 = null
    loop i=1..nf2s -q
      jjj = inlist(b.xlist,formula_2s[i])
      if jjj > 0
        indep0 += b.xlist[jjj]
      endif
    endloop
    formula = smart_formula2(b.xlist,indep0)
  endif
  ###
endif
#funcerr "just stopped"
isnuatwrt = !exists(at_wrt_type)
if isnuatwrt
  string at_wrt_type = ""
endif
strings all_string_input = fndtype(at_wrt_type)
string link = all_string_input[3]
string wrt = all_string_input[1]
wrtBACKUP = wrt
string at_nams =  all_string_input[2]
if strlen(wrt)+strlen(at_nams) == 0
  isnuatwrt = 1
endif
series sample = b.sample
formulas = strsplit3(formula,";")
newformula = formulas[1]
clistnames = formulas[2]
dilistnames = formulas[3]
list clist = null
list dilist = null
nvar = $nvars
#stri = mk_txt_range(nvar - 1)
#funcerr "just stopped"
if nelem(da) == 0
  list da = seq(1,nvar - 1)
endif
if clistnames != "_empty"
  clistnamess = strsplit(clistnames)
  nclist = nelem(clistnamess)
  loop i=1..nclist -q
    clnamsi = clistnamess[i]
    clist+=da[inlist(da,clnamsi)]
  endloop
else
  funcerr "no continuous variables in the model"
endif
matrix t_c_indi = {}
string true_cont = ""
loop foreach i clist -q
  isdi = isdummy(clist[i])
  if !isdi
    true_cont = true_cont~" "~varname(clist[i])
    t_c_indi = t_c_indi~i
  endif
endloop
#funcerr "just stopped"
if cmd != "heckit"
  if strlen(true_cont) == 0
    funcerr "no continuous variables in the model"
  endif
endif
isnumat = !exists(at_vals)
if isnuatwrt
  string at_nams = null
  string wrt = true_cont
else
  if strlen(wrt) == 0
    string wrt = true_cont
  else
    checkwrt = str_setdiff(wrt,true_cont)
    if (strlen(checkwrt) > 0)&&(cmd != "heckit")
      funcerr "non-continuous variable(s) in 'wrt' list"
    endif
  endif
endif
if dilistnames != "_empty"
  dilistnamess = strsplit(dilistnames)
  ndilist = nelem(dilistnamess)
  loop i=1..ndilist -q
    dilnamsi = dilistnamess[i]
    dilist+=da[inlist(da,dilnamsi)]
  endloop
endif
if cmd == "heckit"
  string orig_names = clistnames == "_empty"? "" : clistnames
  orig_names = dilistnames == "_empty"? orig_names : orig_names~" "~dilistnames
  scalar grvers = $version
  string at_namsZ = ""
  string wrtZ = ""
  list zlist = b.zlist - b.zlist[1]
  if !isnuatwrt
    string at_namsALL = at_nams
    at_nams = str_setinrersect(at_nams,orig_names)
    at_namsZ = str_setdiff(at_namsALL,at_nams)
    string wrtALL = wrt
    wrt = str_setinrersect(wrt,orig_names)
    string wrtZ = str_setdiff(wrtALL,wrt)
  endif
  isnw = wrt == ""
  if !isnw
    wtwt = str_setdiff(wrt,orig_names)
    wrtchk = strlen(wtwt)
    errtext = sprintf("wrong name(s) in 'wrt' list: %s",wtwt)
    if wrtchk > 0
      funcerr errtext
    endif
  endif
endif
#     print clist
mclist = {clist}
mdilist = {dilist}
naclist = varnames(clist)
nadilist = varnames(dilist)
#funcerr "just stopped"
if strlen(at_nams) == 0
  bundle bu = defvals (formula,clist,dilist)
else
  bundle bu = findvals(at_nams,at_vals,formula,clist,dilist)
endif
#funcerr "just stopped"
allnames = bu["allnames"]
matrix allvalues = bu["allvalues"]
orerr = 1
orerrlo = 1
merrlo = 1
contvars = bu.contvars
contvars = str_setinrersect(contvars,wrt)
#############
#funcerr "just stopped"
if (cmd == "ols")||(cmd == "hsk")||(cmd == "wls")||(cmd == "tsls")
  #printf "\n%s elasticities of %s:\n", cmd, name_of_y
  if link !="exp"
    printf "\n%s elasticities of %s:\n", cmd, name_of_y
    matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,"ols",allnames,vcov)
  else
    yformula = smart_formula2(b.ylist)
    yformulas = strsplit(yformula,";")
    ckylog = substr(yformulas[1],1,4)=="log("
    #             print ckylog
    #             eval uniqvarlist(yformulas[1])
    #if (strlen(name_of_y) > 2)&&(substr(name_of_y,1,2)=="l_")
    if ckylog
      #new_name_y = substr(name_of_y,3,strlen(name_of_y))
      new_name_y = substr(yformulas[1],5,strlen(yformulas[1])-1)
      printf "\n%s elasticities of %s:\n", cmd, new_name_y
    else
      printf "\n%s elasticities of exp(%s):\n", cmd, name_of_y
    endif
    matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,"poisson",allnames,vcov)
  endif
  #matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,"ols",allnames,vcov)
elif cmd == "poisson"
  printf "\n%s elasticities of E(%s):\n", cmd, name_of_y
  matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,cmd,allnames,vcov)
elif cmd == "logistic"
  printf "\n%s elasticities of Pr(%s = 1):\n", cmd, name_of_y
  matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,"logit",allnames,vcov)
elif cmd == "logit"
  catch scalar ordere = b.ordered
  orerrlo = $error
  if orerrlo == 0
    onx = b.nx
    matrices M = op_elast1 (newformula,allnames,allvalues,wrt,b.coeff,vcov,onx,2)
  else
    catch scalar mulitino = b.multinom
    merrlo = $error
    if merrlo == 0
      matrices M = ml_elast22 (newformula,allnames,allvalues,wrt,b.coeff,mulitino +1,vcov)
    else
      printf "\n%s elasticities of Pr(%s = 1):\n", cmd, name_of_y
      matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,cmd,allnames,vcov)
    endif
  endif
elif cmd == "probit"
  catch scalar ordere = b.ordered
  orerr = $error
  if orerr == 0
    onx = b.nx442w
    printf "\n%s elasticities of Pr(%s = 1):\n", cmd, name_of_y
    matrices M = el_coeff002(newformula,wrt,allvalues,b.coeff,cmd,allnames,vcov)
  endif
elif cmd == "heckit"
  name_of_yz = varnames(b.zlist)[1]
  if link != "Ey"
    printf "\n%s elasticities of E(%s|%s > 0):\n", cmd, name_of_y, name_of_yz
  else
    printf "\n%s elasticities of E(%s):\n", cmd, name_of_y
  endif
  list xxZ = zlist
  exiflaZ = exists(formulaZ)
  if !exiflaZ
    formulaZ = smart_formula2(xxZ)
    ###
    if exifla
      isdotcomma = instring(formula_backup,";")
      if !isdotcomma
        formula_2s = strsplit(formula_backup)
        nf2s = nelem(formula_2s)
        list indep0 = null
        loop i=1..nf2s -q
          jjj = inlist(xxZ,formula_2s[i])
          if jjj > 0
            indep0 += xxZ[jjj]
          endif
        endloop
        formulaZ = smart_formula2(xxZ,indep0)
      endif
    endif
    #####
  endif
  formulasZ = strsplit3(formulaZ,";")
  newformulaZ = formulasZ[1]
  if !exiflaZ
    newformulaZ = strstrip(smartstrsub(newformulaZ,"const",""))
  endif
  clistnamesZ = formulasZ[2]
  dilistnamesZ = formulasZ[3]
  list clistZ = null
  list dilistZ = null
  if clistnamesZ != "_empty"
    clistnamessZ = strsplit(clistnamesZ)
    nclistZ = nelem(clistnamessZ)
    loop i=1..nclistZ -q
      clnamsiZ = clistnamessZ[i]
      clistZ+=da[inlist(da,clnamsiZ)]
    endloop
  endif
  if dilistnamesZ != "_empty"
    dilistnamessZ = strsplit(dilistnamesZ)
    ndilistZ = nelem(dilistnamessZ)
    loop i=1..ndilistZ -q
      dilnamsiZ = dilistnamessZ[i]
      dilistZ+=da[inlist(da,dilnamsiZ)]
    endloop
  endif
  string orig_namesZ = clistnamesZ == "_empty"? "" : clistnamesZ
  orig_namesZ = dilistnamesZ == "_empty"? orig_namesZ : orig_namesZ~" "~dilistnamesZ
  at_namsZ = str_setinrersect(orig_namesZ,at_namsZ)
  if strlen(at_namsZ) == 0
    bundle buZ = defvals (newformulaZ,clistZ,dilistZ)
  else
    matrix findiZZ = findinds(at_namsZ, at_namsALL)
    matrix at_valsZ = at_vals[findiZZ]
    bundle buZ = findvals(at_namsZ,at_valsZ,newformulaZ,clistZ,dilistZ)
  endif
  allnamesZ = buZ.allnames
  contvarsZ = buZ.contvars
  if contvarsZ == "_empty"
    contvarsZ = ""
  endif
  dumvars = bu.dumvars
  if dumvars == "_empty"
    dumvars = ""
  endif
  dumvarsZ = buZ.dumvars
  if dumvarsZ == "_empty"
    dumvarsZ = ""
  endif
  matrix allvaluesZ = buZ.allvalues
  divarsZ = buZ.divars
  if divarsZ == "_empty"
    divarsZ = ""
  endif
  divars = bu.divars
  if divars == "_empty"
    divars = ""
  endif
  if (strlen(wrtBACKUP) > 0)
    contvars = str_setinrersect(contvars,wrtBACKUP)
    contvarsZ = str_setinrersect(contvarsZ,wrtBACKUP)
    dumvars = str_setinrersect(dumvars,wrtBACKUP)
    dumvarsZ = str_setinrersect(dumvarsZ,wrtBACKUP)
    divars = str_setinrersect(divars,wrtBACKUP)
    divarsZ = str_setinrersect(divarsZ,wrtBACKUP)
  endif
  if  (strlen(contvars)+ strlen(contvarsZ)) == 0
    funcerr "no continuous variables"
  else
    matrices M = ALL_hk_fun(newformula,allnames,contvars,newformulaZ, allnamesZ, contvarsZ,b.base_coeffs+1,b.ncoeff,allvalues,allvaluesZ,coe,vcov,1,link)
    second = str_setdiff(contvarsZ,contvars)
    contvars = strstrip(contvars~" "~second)
  endif
elif cmd == "tobit"
  tobcond = 0
  if exists(link)
    if link == "inside"
      tobcond = 1
    endif
  endif
  vcov = b.full_vcv
  s__ = b.sigma
  catch lo__ = b.llimit
  errtob1 = $error
  if errtob1 > 0
    lo__ = 0
  endif
  catch hi__ = b.rlimit
  errtob2 = $error
  if errtob2 > 0
    hi__ = $huge
    scalar hi22 = NA
  else
    scalar hi22 = hi__
  endif
  if tobcond == 0
    printf "\n%s elasticities of E(%s):\n", cmd, name_of_y
    matrices M = tb_elast1(newformula,allnames,allvalues,wrt,coe,s__,lo__,hi__,vcov)
  else
    tytyt = mk_tob_outtext(lo__,hi22,name_of_y)
    printf "\n%s elasticities of E(%s|%s):\n", cmd,name_of_y, tytyt
    matrices M = tb_elast2(newformula,allnames,allvalues,wrt,coe,s__,lo__,hi__,vcov)
  endif
else
  errmess = sprintf("not imlemented yet for %s",cmd)
  funcerr errmess
endif
if (orerr == 0)||(orerrlo == 0)||(merrlo == 0)
  scalar nal = M[5]
  matrix coco = M[2]
  matrix soso = M[4]
  if merrlo == 0
    printf "\nmultinomial %s elasticities\n", cmd
  else
    printf "\nordered %s elasticities\n", cmd
  endif
  loop i=1..nal -q
    printf "\nelasticities of Pr(%s = %d):\n",name_of_y, i
    matrix coeffmat  = coco[i,]'~soso[i,]'
    modprint coeffmat contvars
  endloop
  matrix coeffmat  = vec(coco')~vec(soso')
else
  matrix coeffmat  = M[2]~M[4]
  modprint coeffmat contvars
endif
newnames = contvars
bundle bout
bout["names"] = newnames
matrix bout["mfx"] = coeffmat[,1]
matrix bout["se"] = coeffmat[,2]
matrix bout["vcov"] = M[3]
string bout["type"] = "el"
string bout["pkg"] = "a_eff"
string bout["command"] = b.command
string bout["cmd2"] = cmd
scalar bout["nalt"] = 1
if cmd == "probit"
  if orerr == 0
    bout["cmd2"] = "orderedP"
    bout["nalt"] = b.ncoeff - b.nx + 1
  endif
endif
if cmd == "logit"
  if orerrlo == 0
    bout["cmd2"] = "orderedL"
    bout["nalt"] = b.ncoeff - b.nx + 1
  endif
  if merrlo == 0
    bout["cmd2"] = "multinomial"
    bout["nalt"] =  b.ncoeff/b.cblock + 1
  endif
endif
return bout

list of series list of series whether to include squares

li1 = li1 - 0
li2 = li2 - 0
if squares == 0
  li2 -=li1
  eval varname(li2)
endif
list ret = null
namli1 = varnames(li1)
namli110 = cut_names(namli1)
bign = nelem(namli110)
namli2 = varnames(li2)
namli220 = cut_names(namli2)
bign2 = nelem(namli220)
string descr
strings S = array(2)
scalar count = 0
loop i=1..bign -q
  loop j = 1..bign2 -q
    name = namli110[i]~"_"~namli220[j]
    descr = namli1[i]~"*"~namli2[j]
    uni = uniqvarlist(descr)
    unis = strsplit(uni)
    n = nelem(unis)
    if n == 1
      S = corr_d_n_sq(name,descr)
      name = S[1]
      descr = S[2]
    endif
    ret += genseries(name,li1[i]*li2[j])
    text2 = sprintf("setinfo %s --description=\"%s\"" ,name, descr)
    @text2
  endloop
endloop
return dropcoll(ret)

nf = nelem(flas)
if !exists(x)
  x = "x___"
endif
if iscst
  rets = x~sprintf("[%d]", start)
  loop i=1..nf -q
    flsi = flas[i]
    reti = x~sprintf("[%d]*%s",i+start,flsi)
    rets = rets~"+"~reti
  endloop
else
  rets = ""
  loop i=1..nf -q
    flsi = flas[i]
    reti = x~sprintf("[%d]*%s",i+start-1,flsi)
    if i == 1
      rets = reti
    else
      rets = rets~"+"~reti
    endif
  endloop
endif
return rets

strings flas = strsplit(fla)
strings wrtvarss = strsplit(wrtvars)
if !exists(allvars)
  allvars = uniqvarlist(fla)
endif
if !exists(vect_name)
  vect_name = "x___"
endif
if !exists(b_name)
  b_name = "b___"
endif
strings allvarss = strsplit(allvars)
nf = nelem(flas)
nv = nelem(wrtvarss)
matrix indis = findinds(wrtvars,allvars)
indis = sort(indis)
nindi = nelem(indis)
fla1 = fla
loop i = 1..nindi -q
  subi = sprintf("%s[%d]",vect_name,i)
  fla1 = smartstrsub(fla1,allvarss[indis[i]],subi)
endloop
ret1 = xtextfun(strsplit(fla1),nindi+1,iscst)
ret2 = xtextfun(strsplit(fla1),1,iscst,b_name)
strings S = array(2)
S[1] = ret2
S[2] = ret1
return S

if exists(z___)
  strings S___ = strsplit(vars___)
  loop i___=1..nelem(z___) -q
    ddd___ = S___[i___]
    scalar @ddd___ = z___[i___]
  endloop
endif
scalar val0___ = @fun___
sign___ = val0___ > 0
if link___ == "ols"
  fun___ = "log(abs"~fun___~")"
elif link___ == "poisson"
  fun___ = "abs"~fun___
elif link___ == "logit"
  fun___ = "log(1+exp(-("~fun___~")))"
  sign___ = 0
elif link___ == "probit"
  fun___ = "log(cnorm("~fun___~"))"
  sign___ = 0
endif
nx___ = nelem(x___)
x___ = vec(x___)
jac___ = sign___ ? fdjac(x___,@fun___) : -fdjac(x___,@fun___)
jac___ = vec(jac___)
xbcp___ = x___
x___ = x___|vec(b___)
if link___ == "ols"
  fun2___ = "log(abs"~fun2___~")"
elif link___ == "logit"
  fun2___ = "log(1+exp(-("~fun2___~")))"
elif link___ == "probit"
  fun2___ = "log(cnorm("~fun2___~"))"
endif
#if (link___ == "ols")||(link___ == "poisson")
matrix hess___ = numhess(x___,@fun2___)[1:nx___,(nx___+1):]
matrices M___ = defarray({val0___},jac___.*xbcp___,hess___.*xbcp___)
# elif link___ == "poisson"
#     matrix hess___ = numhess(x___,@fun2___)[1:nx___,(nx___+1):]
#     matrices M___ = defarray({val0___},jac___,hess___)
#endif
return M___

nf = nelem(strsplit(fla))
nb = nelem(b)
iscst = nb > nf
if !exists(allvars)
  allvars = uniqvarlist(fla)
endif
indi0 = findinds(wrtvars,allvars)
wrtvars = str_setinrersect(allvars,wrtvars)
if wrtvars == allvars
  allvals = allvals[indi0]
endif
strings work_flass =make_flas_el (fla, wrtvars,allvars,"x___","b___",iscst)
el_fla = "("~work_flass[1]~")"
hess_fla = "("~work_flass[2]~")"
zvars = str_setdiff(allvars,wrtvars)
findix = findinds(wrtvars,allvars)
if strlen(zvars) > 0
  findiz = findinds(zvars,allvars)
  matrix x = allvals[sort(findix)]
  matrix z = allvals[findiz]
  matrices M = super3(el_fla,hess_fla ,zvars,x,b,z,cmd)
else
  matrix x = allvals
  matrices M = super3(el_fla,hess_fla ,null,x,b,null,cmd)
endif
if exists(vcv)
  matrices retms = array(4)
  retms[1] = M[1]
  retms[2] = M[2]
  matrix vcov = qform(M[3],vcv)
  retms[3] = vcov
  matrix ses = sqrt(diag(vcov))
  retms[4] = ses
else
  matrices retms = array(2)
  retms[1] = M[1]
  retms[2] = M[2]
endif
return retms

ni = nelem(inds)
nnn = n - ni
if nnn == 0
  matrix ret = {}
else
  scalar mu
  matrix ret = zeros(1,nnn)
  scalar j = 1
  loop i = 1..n -q
    mu = minc(abs(inds-i))
    if mu > 0
      ret[j] = i
      j++
    endif
  endloop
endif
return ret

bundle bu = defvals (formula,clist,dilist)
matrix matri = bu["allvalues"]
namss = strsplit(nams)
n = nelem(namss)
matrix thelev={}
matrix indi={}
loop i=1..n -q
  i1 = findind(namss[i],bu["contvars"])
  i2 = findind(namss[i],bu["dumvars"])
  if i2
    i24 = minr(abs({0,1} .- vals[i]))
    if i24>0
      msg = "'"~namss[i]~"': wrong level for 'at' value"
      funcerr msg
    endif
  endif
  if (i1+i2)
    i11 = findind(namss[i],bu["allnames"])
    matri[i11]=vals[i]
  else
    i3 = findind(namss[i],bu["divars"])
    alli = values(dilist[i3])
    i4 = minc(abs(alli .- vals[i]))
    if i4>0
      msg = "'"~namss[i]~"': wrong level for 'at' value"
      funcerr msg
    endif
    if i3==0
      msg = "'"~namss[i]~"': wrong name of 'at' variable"
      funcerr msg
    endif
    levi = mk_factlevs2 (formula, namss[i],bu["allnames"])
    dini = discrnams2(namss[i],levi)
    thelev = getatdidvals(levi, vals[i])
    indi = findinds(dini,bu["allnames"])
    matri[indi]=thelev
  endif
endloop
matrix bu["allvalues"]=matri
return bu

ch1=(substr(str1,1,2)=="l_")
ch2 = exists(othernames)
ret = 0
if (ch1 > 0)&&(ch2 > 0)
  tt2 = str_setinrersect(str2,othernames)
  tt2s = strsplit(tt2)
  ttn = nelem(tt2s)
  if ttn==1
    tt21 = tt2s[1]
    ntt21 = strlen(tt21)
    if ntt21 > 0
      ret = 1
    endif
  endif
endif
if ret == 0
  sret = "_empty"
else
  sret = "log("~tt21~")"
endif
return sret

ni = !exists(indeps)
if ni
  indeps = ""
endif
fla0 = fla
vari = uniqvarlist(fla)
varis = strsplit(vari)
nvari = nelem(varis)
contis = ""
diis = ""
loop i=1..nvari -q
  vari = varis[i]
  fifi = findind(vari,indeps)
  if fifi != 0
    contis = contis~" "~vari
  else
    isduvi = isduname(vari)
    if isduvi
      fla0 = smartstrsub(fla0,vari,corr_duname(vari))
      diis = diis~" "~getdiname(corr_duname(vari))
    elif isdiname(vari)
      diis = diis~" "~getdiname(vari)
    else
      contis = contis~" "~vari
    endif
  endif
endloop
contis = strstrip(contis)
diis = strstrip(diis)
ret = fla0~";#"~contis~";#"~diis
rets = strsplit3(ret,";")
return rets

if strlen(strstr(bigst,find))
  varch = "qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM0123456789_"
  varchR = "qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM0123456789_(["
  bigst1 = "#"~strsub(bigst," ", "#")~"#"
  bigst2 = strsub(bigst1,find," "~find~" ")
  bigst2s = strsplit(bigst2)
  bign = nelem(bigst2s)
  ret0 = ""
  loop i = 1..bign -q
    if bigst2s[i] == find
      lengi = strlen(bigst2s[i-1])
      left = strlen(strstr(varch,substr(bigst2s[i-1],lengi,lengi)))
      right = strlen(strstr(varchR,substr(bigst2s[i+1],1,1)))
      if !left&&!right
        bigst2s[i] = repl
      endif
    endif
    ret0 = ret0~bigst2s[i]
  endloop
  return strstrip(strsub(ret0,"#"," "))
else
  return bigst
endif

if !exists(collapce)
  collapce = " "
endif
n = nelem(strs)
ret = strs[1]
if n>1
  loop i = 2..n -q
    ret = ret~collapce~strs[i]
  endloop
endif
ret = strstrip(ret)
return ret

levss = strsplit(levs)
n = nelem(levss)
matrix mat = zeros(1,n)
loop i=1..n -q
  nlevi = atof(levss[i])
  mat[i]=mean(y==nlevi)
endloop
return mat

ret = ""
loop i=1..n -q
  reti = sprintf(" %d",i)
  ret = ret~reti
endloop
ret = strstrip(ret)
return ret

levss = strsplit(levs)
ns = nelem(levss)
txt = ""
loop i=1..ns -q
  if i==1
    sep = ""
  else
    sep = ", "
  endif
  txt = txt~sep~nam~" = "~levss[i]
endloop
txt = strstrip(txt)
if ns==1
  txt=txt~","
endif
return txt

n = strlen(colname(x,1))
ret = ""
if n > 0
  nx = nelem(x)
  loop i=1..nx -q
    ret = ret~" "~colname(x,i)
  endloop
  ret = strstrip(ret)
  rets = strsplit(ret)
  nr = nelem(rets)
  if nr != nx
    ret = ""
  endif
endif
lr = strlen(ret)>0
ret1 = ret
if ret1 == ""
  ret1 = "_empty"
endif
return ret1

ops = "+ - * / ^"
opss = strsplit(ops)
posi = zeros(5,1)
loop i=1..5 -q
  posi[i]=isstrstr(str,opss[i])
endloop
return sum(posi)>0

str0 = str
bign = strlen(str)
if bign < 76
  printf "\n%s\n", str
else
  posi = posstrstr(str0," ")
  if posi[1]==0
    lines = ceil(bign/76)
    printf "\n"
    loop i=1..lines -q
      stri = strstrip(substr(str,(i-1)*76+1,i*76))
      printf "%s\n", stri
    endloop
  else
    j = 1
    str10 = str
    printf "\n"
    loop while 1==1 -q
      n10 = strlen(str10)
      if n10<77
        strprn = str10
      else
        str0 = substr(str10,1,77)
        posii = posstrstr(str10," ")
        break1 = iminr(abs(posii.-76))
        strprn = strstrip(substr(str10,1,posii[break1]))
        str10 = substr(str10,posii[break1]+1,n10)
      endif
      printf "%s\n",strprn
      j++
      if n10<77
        break
      endif
    endloop
  endif
endif

big small

xx = strstr(s1,s2)
nxx = strlen(xx)
n = strlen(s2)
if nxx<(n+1)
  string ret = ""
else
  string ret = substr(xx,n+1,nxx)
endif
return ret

strings strs = strsplit(str)
strings bigstrs = strsplit(bigstr)
scalar bn = nelem(bigstrs)
string ret = ""
scalar fi
loop i = 1..bn -q
  fi = findind(bigstrs[i],str)
  if fi == 0
    ret = ret~" "~bigstrs[i]
  endif
endloop
ret = strstrip(ret)
return ret

dnam = discrnams2(nam,factlevs)
strings dnams = strsplit(dnam)
nf = nelem(dnams)
nv = nelem(allval)
matrix indi = zeros(1,nf)
loop i=1..nf -q
  indi[i]=findind(dnams[i],allnam)
endloop
bign = nelem(strsplit(alllevs))
matrix Z = I(nf)
matrix z ={}
if nf == bign
  matrix ret = zeros(nf,nv)
  loop i=1..nf -q
    z = allval
    z[indi] = Z[i,]
    ret[i,] = z
  endloop
else
  matrix ret = zeros(nf+1,nv)
  Z = zeros(1,nf)|Z
  loop i=1..(nf+1) -q
    z = allval
    z[indi] = Z[i,]
    ret[i,] = z
  endloop
endif
return ret

ret = 0
n = strlen(str)
if n > 0
  matrix ret0 = zeros(1,n)
  loop i=1..n -q
    stri = substr(str,i,i)
    isdi = isdigit (stri)
    if i==1
      isdi=isdi&&(stri!="0")
    endif
    ret0[i] = isdi
  endloop
  ret = (sum(ret0)==n)
endif
return ret

ret = 0
ck1 = substr(name,1,1)=="D"
if ck1
  nn = strlen(name)
  name2 = substr(name,2,nn)
  ret = isdiname(name)
endif
return ret

poseq = posstrstr(nam,"=")
nposeq = nelem(poseq)
ret=""
if nposeq > 1
  nams = strsplit3(nam,";")
  loop i=1..nposeq -q
    namss = strsplit3(nams[i],"=")
    namss1 = strstrip(namss[1])
    if namss1 == "wrt"
      ret = strstrip(namss[2])
    endif
  endloop
endif
if nposeq == 1
  nams = strsplit3(nam,"=")
  ck = strstrip(nams[1]) == "wrt"
  if ck
    ret = strstrip(nams[2])
  endif
endif
return ret

n = nelem(strs)
counter = 1
loop i=1..n -q
  strsi = strs[i]
  if strlen(strsi) > 9
    one = substr(strsi,1,3)
    snumber = sprintf("%d", i)
    sle = strlen(snumber)
    zrs = repchar("0",6-sle)
    two = zrs~snumber
    strs[i] = one~two
  endif
endloop
return strs

if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
smallb = b
b = b|s__
nams = strsplit(nam)
nn = nelem(nams)
nb = nelem(b)
nns = nn
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
loop i=1..nns -q
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
endif
matrix xbackup = x
funs2 = strsplit(fun2)
nf = nelem(funs2)
matrix ret222 = zeros(1,nf)
nnxx = nelem(x)
matrix ret = zeros(nn,nb+1)
loop j = 1..nnxx -q
  x[j] = 0
  matrix grads = zeros(2,nf)
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads[1,i] = @fun3
  endloop
  x[j] = 1
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads[2,i] = @fun3
  endloop
  if nb > nf
    grads = ones(2,1)~grads
  endif
  matrix xbs = grads*smallb
  cnorm1 = cnorm((lo__-(xbs))/s__)
  cnorm2 = cnorm(-(hi__-(xbs))/s__)
  cnorm3 = cnorm((hi__-(xbs))/s__)
  dnorm1 = dnorm((lo__-(xbs))/s__)
  dnorm2 = dnorm((hi__-(xbs))/s__)
  ddd2 = lo__*cnorm1+hi__*cnorm2 +(cnorm3-cnorm1).*((xbs) + s__*(dnorm1-dnorm2)./ (cnorm3-cnorm1))
  ret000 = ddd2[2] - ddd2[1]
  fxbs = cnorm3-cnorm1
  gxbs = dnorm2-dnorm1
  dgxbs = gxbs[1] - gxbs[2]
  matrix retj = ret000 ~ (fxbs[2]*grads[2,] - fxbs[1]*grads[1,]) ~ dgxbs
  x = xbackup
  ret[j,] = retj
endloop
return ret

strs = strsplit(str)
n = nelem(strs)
matrix ret = zeros(1,n)
loop i=1..n -q
  ret[i]=findind(strs[i],bigstr)
endloop
return ret

strlev = sprintf("%d",lev)
n = nelem(strsplit(levs))
ind = findind(strlev,levs)
matrix ret = zeros(1,n)
if (ind>0)&&(ind<=n)
  ret[ind]=1
endif
return ret

ch1=(substr(str1,1,3)=="sq_")
ch2 = exists(othernames)
ret = 0
if (ch1 > 0)&&(ch2 > 0)
  tt2 = str_setinrersect(str2,othernames)
  tt2s = strsplit(tt2)
  ttn = nelem(tt2s)
  if ttn==1
    tt21 = tt2s[1]
    ntt21 = strlen(tt21)
    if ntt21 > 0
      ret = 1
    endif
  endif
endif
if ret == 0
  sret = "_empty"
else
  sret = tt21~"^2"
endif
return sret

fla0 = fla
vari = uniqvarlist(fla)
varis = strsplit(vari)
nvari = nelem(varis)
conti = ""
dii = ""
loop i=1..nvari -q
  vari = varis[i]
  isduvi = isduname(vari)
  if isduvi
    fla0 = smartstrsub(fla0,vari,corr_duname(vari))
  endif
endloop
return fla0

ops = "> < >= <= | || & && = ==  % !="
opss = strsplit(ops)
posi = zeros(nelem(opss),1)
loop i=1..nelem(opss) -q
  posi[i]=isstrstr(str,opss[i])
endloop
return sum(posi)>0

matrix Z = substvals2(allnam,allval,nam,factlevs,alllevs )
rz = rows(Z)
nb = nelem(b)
smallb = b
b = b|s__
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
string NAM = ""
strings  factlevss = strsplit(factlevs)
nns = nelem(factlevss)
strings nams = array(nns)
matrix indi = zeros(nns,1)
string fun2 = mfun
loop i = 1.. nns -q
  nams[i] = nam~"_"~factlevss[i]
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
NAM = strjoin(nams)
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
matrix z = {}
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,NAM)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
matrix xbackup = x
funs2 = strsplit(fun2)
nf = nelem(funs2)
matrix ret222 = zeros(1,nf)
nnxx = nelem(x)
matrix ret = zeros(rz,nb+2)
loop j = 1..rz -q
  x = Z[j,indi]
  matrix grads = zeros(1,nf)
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads[i] = @fun3
  endloop
  if nb > nf
    grads = 1~grads
  endif
  matrix xbs = grads*smallb
  cnorm1 = cnorm((lo__-(xbs))/s__)
  cnorm2 = cnorm(-(hi__-(xbs))/s__)
  cnorm3 = cnorm((hi__-(xbs))/s__)
  dnorm1 = dnorm((lo__-(xbs))/s__)
  dnorm2 = dnorm((hi__-(xbs))/s__)
  ret000 = lo__*cnorm1+hi__*cnorm2 +(cnorm3-cnorm1).*((xbs) + s__*(dnorm1-dnorm2)./ (cnorm3-cnorm1))
  fxbs = cnorm3-cnorm1
  gxbs = dnorm1-dnorm2
  matrix retj = ret000 ~ (fxbs*grads) ~ gxbs
  #x = xbackup
  ret[j,] = retj
endloop
return ret[2:,] .- ret[1,]

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
string fun2 = fun
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
string fun223 = "{"~strsub(fun2," ",",")~"}"
matrix svals = @fun223
if nns < bign
  matrix grads0 = fdjac(x,for_djac(x,fun223,z))'
  if ismucst
    matrix grads0 = zeros(nns,1)~ grads0
  endif
  matrix zzz = ml_cont_grads56 ( b,nalts,&svals,&grads0,ismucst)
  matrix uuu = fdjac(b,ml_cont_grads56( b,nalts,&svals,&grads0,ismucst))
else
  matrix grads0 = fdjac(x,for_djac(x,fun223))'
  if ismucst
    matrix grads0 = zeros(nns,1)~ grads0
  endif
  matrix zzz = ml_cont_grads56 ( b,nalts,&svals,&grads0,ismucst)
  matrix uuu = fdjac(b,ml_cont_grads56( b,nalts,&svals,&grads0,ismucst))
endif
return zzz~uuu

matrix b = vec(b)
matrix Z = substvals2(allnam,allval,nam,factlevs,alllevs )
rz = rows(Z)
nb = nelem(b)
strings  factlevss = strsplit(factlevs)
nns = nelem(factlevss)
strings nams = array(nns)
matrix indi = zeros(nns,1)
string NUM = ""
loop i = 1.. nns -q
  nams[i] = nam~"_"~factlevss[i]
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = mfun
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
NAM = strjoin(nams)
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,NAM)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
string prefun2s = "("~strsub(fun2," ",") (")~")"
string fun23 = "{"~strsub(prefun2s," ",",")~"}"
matrix svals
matrix ret0 = zeros(rz,nalts)
matrix ret20 = zeros(rz,nb*nalts)
matrix seqnal = seq(1,nalts)
loop j=1..rz -q
  x = Z[j,indi]
  svals = @fun23
  ret0[j,] =  mc_p17(b,svals,nalts,ismucst)
  matrix dpij1 = fdjac(b, mc_p17(b,svals,nalts,ismucst))
  ret20[j,] = vec(dpij1')'
endloop
ret0 = vec(ret0[2:,] .- ret0[1,])
ret20 = ret20[2:,] .- ret20[1,]
matrix ret30 = zeros((rz-1)*nalts,nb)
matrix newseqi = seq(1,rz-1)
matrix seqb = seq(1,nb)
loop i = 1..nalts -q
  ret30[(i-1)*(rz-1)+newseqi,] = ret20[,(i-1)*nb+seqb]
endloop
return ret0 ~ ret30

big small

xx = strstr(s1,s2)
nxx = strlen(xx)
ret = (nxx>0)
return ret

big small

xx = cutstrstr(s1,s2)
nxx = strlen(xx)
ww = isstrstr(s1,s2)
if nxx + ww == 0
  ret = {0}
else
  ret = {strlen(s1)-strlen(strstr(s1,s2))+1}
  xx = s1
  loop while nxx>0 -q
    xx = cutstrstr(xx,s2)
    nxx = strlen(xx)
    if nxx>0&&isstrstr(xx,s2)
      ret = ret~(strlen(s1)-strlen(strstr(xx,s2))+1)
    endif
  endloop
endif
return ret

ret = ""
loop i=1..n -q
  ret=ret~char
endloop
return ret

scalar isn = !exists(sep)
if isn||sep == " "
  ret = strsplit(str)
else
  string str2 = strsub(str," ","")
  str2 = strsub(str2,sep," ")
  ret = strsplit(str2)
endif
return ret

unils = strsplit(unil)
n = nelem(unils)
nm = name~"_"
matrix ret = {}
loop i=1..n -q
  isi = isstrstr(unils[i],nm)
  if isi
    str2 = strsub(unils[i],"_"," ")
    nmsi = strsplit(str2)
    nmi = nmsi[2]
    ret = ret~atof(nmi)
  endif
endloop
ret = sort(ret)
n2 = nelem(ret)
ret2 = ""
loop i=1..n2 -q
  reti = sprintf("%d",ret[i])
  ret2 = ret2~" "~reti
endloop
ret2 = strstrip(ret2)
return ret2

matrix inds = findinds(str1,str2)
str1s = strsplit(str1)
n = nelem(str1s)
ret = ""
loop i=1..n -q
  if  inds[i]
    ret = ret~" "~str1s[i]
  endif
endloop
ret = strstrip(ret)
return ret

if !exists(dinames)&&!exists(othernames)
  funcerr "no names given"
endif
iseq = (substr(str,1,1)=="=")
n = strlen(str)
if iseq
  ret3 = "_empty"
  sqr = check_sq (str20,str,othernames)
  if sqr != "_empty"
    ret3 = sqr
  else
    logo = check_log (str20,str,othernames)
    if logo != "_empty"
      ret3 = logo
    else
      times1 = check_times (str20,str,othernames)
      if times1 != "_empty"
        ret3 = times1
      endif
    endif
  endif
  if ret3=="_empty"
    funcerr "incomplete indep. vars. list"
  endif
else
  intera =  substr(str,1,15) == "interaction of "
  if intera
    isdi = !exists(dinames)
    if isdi
      tei1 = "interactions in description "
      tei2 = "but no discrete variables"
      emsg = tei1~tei2
      funcerr emsg
    endif
    ret33 = substr(str,16,n)
    ret33s = strstr(ret33," and ")
    n33 = strlen(ret33)
    n33s = strlen(ret33s)
    ret3a = substr(ret33,1,n33-n33s)
    ret3b = substr(ret33s,7,n33s-1)
    ret3bs = strsplit2(ret3b,"==")
    ret3c = strjoin(ret3bs,"_")
    isund = isstrstr(ret3a,"_")
    if isund
      sret3a = strsplit2(ret3a,"_")
      sret3a1 = sret3a[1]
      isD = substr(sret3a1,1,1)=="D"
      ret3d = ret3a
      if isD
        n3a1 = strlen(sret3a1)
        sret3a2 = substr(sret3a1,2,n3a1)
        iind = findind(sret3a2,dinames)
        if iind
          n3a = strlen(ret3a)
          ret3d = substr(ret3a,2,n3a)
        endif
      endif
    else
      ret3d = ret3a
    endif
    ret3 = ret3d~"*"~ret3c
  else
    isnotdi = !exists(dinames)
    isdu = "_empty"
    if (!isnotdi)
      isdu =  check_dum(str20,str,dinames)
    endif
    if !(isdu=="_empty")
      ret3 = isdu
    else
      ret3=str
      isspace = isstrstr(str," ")
      if isspace
        txt1 = str
        txt2 = "'"~str~"' looks like a description,"
        txt2 = txt2~"it is not a variable name or formula"
        funcerr txt2
      endif
      str2 = uniqvarlist(str)
      isnuoth = !exists(othernames)
      if isnuoth
        othernames = ""
      endif
      fis = findinds(str2,othernames)
      nfis = nelem(fis)
      str2s = strsplit(str2)
      loop i= 1..nfis -q
        ck50 = (fis[i]==0)
        if ck50
          txt60 = "Warning: '"~str2s[i]
          txt62 = "' is not in the list of independent variables"
          txt63 = txt60~txt62
          printf "\n%s\n",txt63
        endif
      endloop
    endif
  endif
endif
return ret3

ret = ""
str = name
if isduname(name)
  str = corr_duname(name)
endif
bign = strlen(str)
matrix mapo = posstrstr(str,"_")
ck2 = mapo!=0
ck3 = 0
if (ck2>0)
  nma = nelem(mapo)
  lapo = mapo[nma]
  level = substr(str,lapo+1,bign)
  ck3 = isnat(level)
  if ck3>0
    ret = substr(str,1,lapo-1)
  endif
endif
return ret

uni = uniqvarlist(des)
unis = strsplit(uni)
n = nelem(unis)
if n == 1
  retname = "cb_"~name1
  retdes = uni~"^3"
elif n == 2
  des2s = strsplit3(des,"*")
  des2= strjoin(des2s)
  matrix fi1 = findind3(unis[1],des2)
  matrix fi2 = findind3(unis[2],des2)
  matrix ma = {nelem(fi1),nelem(fi2)}
  ima = imaxr(ma)
  imi = iminr(ma)
  retdes = unis[ima]~"^2*"~unis[imi]
  smallnames = name1~" "~name2~" "~name3
  matrix fi1 = findind3(name1,smallnames)
  matrix fi2 = findind3(name2,smallnames)
  matrix fi3 = findind3(name3,smallnames)
  matrix ma2 = {nelem(fi1),nelem(fi2),nelem(fi3)}
  ima2 = imaxr(ma2)
  imi2 = iminr(ma2)
  smallnamess = strsplit(smallnames)
  retname = "sq_"~smallnamess[ima2]~"_"~smallnamess[imi2]
else
  retname = name1~"_"~name2~"_"~name3
  retdes = des
endif
strings S = array(2)
S[1]= retname
S[2] = retdes
return S

n = strlen(s)
if n != 1
  ret = 0
else
  string lts = "0123456789"
  ret = isstrstr(lts,s)
endif
return ret

nx = nelem(x)
cx = cols(x)
rx = rows(x)
mrc = xmin(cx,rx)
if mrc > 1
  funcerr "'x' should be a column or row vector"
endif
ninds = nelem(inds)
if ninds > nx
  funcerr "you can't exclude more elements than the length of vector"
endif
inds = vec(inds)
mi = minc(inds)
mai = maxc(inds)
if (mi < 1)||(mai > nx)
  funcerr "subscript(s) is(are) out of range"
endif
if ninds == nx
  matrix ret = {}
else
  matrix eind = indclude(nx,inds)
  matrix ret = x[eind]
endif
return ret

levss = strsplit(factlevs)
n = nelem(levss)
ret = ""
loop i=1..n -q
  reti = nam~"_"~levss[i]
  ret = ret~" "~reti
endloop
ret = strstrip(ret)
return ret

ch2 = exists(othernames)
ret = 0
if (ch2 > 0)
  tt2 = str_setinrersect(str2,othernames)
  tt2s = strsplit(tt2)
  ttn = nelem(tt2s)
  if ttn==2
    tt21 = tt2s[1]
    ntt21 = strlen(tt21)
    tt22 = tt2s[2]
    ntt22 = strlen(tt22)
    if (ntt21 > 0)&& (ntt22 > 0)
      ret = 1
    endif
  endif
endif
if ret == 0
  sret = "_empty"
else
  sret = tt21~"*"~tt22
endif
return sret

ns = strlen(str)
posi = posstrstr(str,"_")
nposi = nelem(posi)
ck1 = (posi[1]!=0)
ret = 0
if ck1
  lp = posi[nposi]
  if lp < ns
    ckck = substr(str,lp+1,ns)
    ret = isnat(ckck)
  endif
endif
if nposi>=2
  ck3 = posi[nposi]-posi[nposi-1]
  if ck3<2
    ret=0
  endif
endif
return ret

matrix posi = posstrstr(vvv,"_")
np = nelem(posi)
if posi[1]==0
  ret = 0
else
  nvvv = strlen(vvv)
  matrix indii = {}
  loop i=1..np -q
    part1 = substr(vvv,1,posi[i]-1)
    part2 = substr(vvv,posi[i]+1,nvvv)
    wordi = part1~" "~part2
    findis = findinds(wordi,uuu)
    timis = sum(findis .>0)==2
    if timis==1
      indii = indii~posi[i]
    endif
  endloop
  if (nelem(indii)==1)&&(indii[1]>0)
    ret = indii[1]
  else
    ret = 0
  endif
endif
return ret

uni = uniqvarlist(des)
unis = strsplit(uni)
n = nelem(unis)
retname = name
retdes = des
if n == 1
  strn = strlen(name)
  smalln = floor(strn/2)
  smallname = substr(name,1,smalln)
  retname = "sq_"~smallname
  retdes = uni~"^2"
endif
strings S = array(2)
S[1]= retname
S[2] = retdes
return S

sprintf ret "%d",n
return ret

strs = strsplit(bigstr)
n = nelem(strs)
scalar ret = 0
loop i=1..n -q
  if smallstr==strs[i]
    ret = i
    break
  endif
endloop
return ret
#     function scalar findind (string smallstr,
#                          string bigstr,
#                          string sep[null])
#     if !exists(sep)
#         sep = " "
#     endif
#     strs = strsplit2(bigstr,sep)
#     n = nelem(strs)
#     matrix ret = {}
#     loop i=1..n -q
#         if smallstr==strs[i]
#             ret = ret~i
#         endif
#     endloop
#     nret = nelem(ret)
#     if nret == 0
#         ret = 0
#     endif
#     if nret > 1
#         funcerr "duplicated names"
#     endif
#     scalar retu = ret[1,1]
#     return retu

#funcerr "just stopped"
ret1 = ""
mret1 = {}
dret1 = ""
mdret1 = {}
n1 = nelem(clist)
if n1>0
  cvns = varnames(clist)
  loop i=1..n1 -q
    isdui = isdummy(clist[i])
    if isdui
      dret1 = dret1~cvns[i]~" "
      mdret1 = mdret1~mean(clist[i])
    else
      ret1 = ret1~cvns[i]~" "
      mret1 = mret1~mean(clist[i])
    endif
  endloop
  dret1 = strstrip(dret1)
  if dret1==""
    dret1a="_empty"
  endif
  ret1 = strstrip(ret1)
  if ret1==""
    ret1a="_empty"
  endif
endif
if dret1==""
  dret1a="_empty"
else
  dret1a=dret1
endif
if ret1==""
  ret1a="_empty"
else
  ret1a=ret1
endif
n2 = nelem(dilist)
ret2 = ""
mret2 = {}
if n2>0
  uli = uniqvarlist(formula)
  divns = varnames(dilist)
  matrix vali = {}
  loop i=1..n2 -q
    levi = mk_factlevs2 (formula, divns[i],uli)
    dini = discrnams2(divns[i],levi)
    vali = getdidvals(levi,dilist[i])
    ret2 = ret2~" "~dini
    mret2 = mret2~vali
  endloop
  ret2 = strstrip(ret2)
endif
if ret2==""
  ret2a="_empty"
else
  ret2a=strjoin(varnames(dilist))
endif
nam1 = ret1~" "~dret1
nam1 = strstrip(nam1)
matrix mat1 = mret1~mdret1
anams = strstrip(nam1~" "~ret2)
matrix avals = mat1~mret2
bundle bu
bu["contvars"]=ret1a
bu["dumvars"]=dret1a
bu["divars"]=ret2a
bu["allnames"]=anams
matrix bu["allvalues"]=avals
return bu

levss = strsplit(levs)
ns = nelem(levss)
txt = ""
loop i=1..ns -q
  txt = txt~" "~nam~"_"~levss[i]
endloop
txt = strstrip(txt)
return txt

matrix posleft = posstrstr(des,"(")
matrix posright = posstrstr(des,")")
ckp1 = nelem(posleft)==nelem(posright)
if (posleft[1]==1)||(posright[nelem(posright)]==nelem(posright))
  ckp1 = 0
endif
ckp2 = 0
if ckp1>0
  ma = posright - posleft
  mima = minr(ma)
  if mima>=0
    ckp2 = 1
  endif
endif
ck3 = 0
if ckp2>0
  uni = uniqvarlist(des)
  ck3 = uni!=des
  des2 = strsub(des,"("," ")
  des2 = strsub(des2,")"," ")
  uni2 = uniqvarlist(des2)
  funs = str_setdiff(uni2,uni)
  if (strlen(funs)==0)&&isop(des)==0
    ck3 = 0
  endif
endif
if islogical(des)
  ck3 = 0
endif
return ck3

ret = 0
posi = posstrstr(name,"_")
nposi = nelem(posi)
#print nposi
bign = strlen(name)
if nposi>1
  loop i=1..nposi -q
    str1 = substr(name,1,posi[i]-1)
    str2 = substr(name,posi[i]+1,bign)
    ck11 = isduname(str1)
    ck12 = isdiname(str1)
    ck21 = isduname(str2)
    ck22 = isdiname(str2)
    #print i ck11 ck12 ck21 ck22
    if (ck11+ck12+ck21+ck22)>0
      ret = 1
    endif
  endloop
endif
return ret

fla2 = "#"~fla~"#"
chars = "qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM"
digits = "0123456789"
lefla = strlen(fla2)
matrix charmap = zeros(lefla,1)
matrix digmap = charmap
matrix bramap = charmap
matrix undmap = charmap
loop i = 1..lefla -q
  chari = substr(fla2,i,i)
  if strlen(strstr(chars,chari))
    charmap[i] =  1
  endif
  if strlen(strstr(digits,chari))
    digmap[i] =  1
  endif
  if chari == "_"
    undmap[i] = 1
  endif
  if (chari == "(")||(chari == "[")
    bramap[i] = 1
  endif
endloop
all_charmap = charmap + digmap + undmap
matrix seql = seq(0,lefla-1)'
matrix begins = seql .*(diff(all_charmap).> 0)
matrix ends = seql.*mlag((diff(all_charmap).< 0),-1)
if  sumc(uniq(begins))
  matrix numbegins = uniq(begins)[2:]
  matrix numends = uniq(ends)[2:]
  smalln = nelem(numends)
  ret = ""
  strings rets
  loop i =1..smalln -q
    summ = 0
    reti = substr(fla,numbegins[i],numends[i])
    nofun = substr(fla,numends[i]+1,numends[i]+1)!="("
    if nofun
      if i==1
        rets +=reti
        ret = ret~" "~reti
      else
        nrets = nelem(rets)
        summ = 0
        loop j=1..nrets -q
          summ += reti == rets[j]
        endloop
        if !(digmap[numbegins[i]+1]+bramap[numends[i]+2]+summ)
          rets +=reti
          ret = ret~" "~reti
        endif
      endif
    endif
  endloop
  return strstrip(ret)
else
  return "_empty"
endif

return x

ret = "_empty"
str1s = strsplit3(str1,"_")
n1 = nelem(str1s)
ck11 = n1 > 1
ck12 = 0
ck13 = 0
if ck11
  str1s1 = str1s[1]
  str1s2 = str1s[n1]
  ck12 = isnat(str1s2)
  ck13 = (substr(str1s1,1,1) == "D")
endif
ck21 = 0
ck22 = 0
ck23 = 0
ck24 = 0
if ck11&&ck12&&ck13
  str2s = strsplit(str2)
  n2 = nelem(str2s)
  ck21 = n2 > 3
  if ck21
    strx = str2s[n2-2]
    ck22 = findind(strx,dinames) > 0
    stry = str2s[n2-1]
    ck23 = (stry == "=")
    strz = str2s[n2]
    ck24 = isnat(strz)
    if ck22&&ck23&&ck24
      ret = strx~"_"~strz
    endif
  endif
endif
return ret

ret = name
if isduname(name)
  bign = strlen(name)
  ret = substr(ret,2,bign)
endif
return ret

poseq = posstrstr(nam,"=")
nposeq = nelem(poseq)
ret=""
if nposeq > 1
  nams = strsplit3(nam,";")
  loop i=1..nposeq -q
    namss = strsplit3(nams[i],"=")
    namss1 = strstrip(namss[1])
    if namss1 == "at"
      ret = strstrip(namss[2])
    endif
  endloop
endif
if nposeq == 1
  #print nam
  nams = strsplit3(nam,"=")
  ck = strstrip(nams[1]) == "at"
  if ck
    ret = strstrip(nams[2])
  endif
endif
#printf "ret: "
#print ret
return ret

if !exists(sep)
  sep = " "
endif
strs = strsplit3(bigstr,sep)
n = nelem(strs)
matrix ret = {}
loop i=1..n -q
  if smallstr==strs[i]
    ret = ret~i
  endif
endloop
#nret = nelem(ret)
return ret

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
string fun2 = fun
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
string prefun2s = "myx("~strsub(fun2," ",") myx(")~")"
strings fun2s = strsplit(prefun2s)
nfs = nelem(fun2s)
matrix svals = zeros(1,nfs)
matrix grads0 = zeros(nns,nfs)
loop i = 1..nfs -q
  funsi = fun2s[i]
  svals[i] = @funsi
  grads0[,i] = fdjac(x,@funsi)'
endloop
if nelem(b) > nfs
  svals = 1~svals
  grads0 = zeros(nns,1)~grads0
endif
matrix dxb_dx =  grads0 * b
scalar  xb = svals*b
scalar  term1 =  (hi__-xb)/s__
scalar term2 = (lo__-xb)/s__
scalar fxb = cnorm(term1)-cnorm(term2)
scalar gxb = dnorm(term1)-dnorm(term2)
matrix ret222 = (-gxb*(dxb_dx .*svals)/s__ + fxb *grads0)~(gxb*dxb_dx*xb/s__^2)
matrix  ret000 = fxb*dxb_dx
return ret000 ~ ret222

strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
string fun2 = fun
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
string prefun2s = "("~strsub(fun2," ",") (")~")"
matrix svals
xbackup = x
matrix ret0 = zeros(nns,nalts)
matrix ret20 = zeros(nns*nalts,nelem(b))
matrix r20seq = nns*seq(0,nalts-1)
string fun23 = "{"~strsub(prefun2s," ",",")~"}"
loop j=1..nns -q
  x[j] = 1
  svals = @fun23
  matrix pij1 = mc_p17(b,svals,nalts,ismucst)
  matrix dpij1 = fdjac(b, mc_p17(b,svals,nalts,ismucst))
  x[j] = 0
  svals = @fun23
  matrix pij0 = mc_p17(b,svals,nalts,ismucst)
  matrix dpij0 = fdjac(b, mc_p17(b,svals,nalts,ismucst))
  x = xbackup
  ret0[j,] =  pij1 - pij0
  ret20[j+r20seq,] = dpij1 - dpij0
endloop
return vec(ret0) ~ret20

scalar isn = !exists(sep)
if isn||sep==" "
  rets = strsplit(str)
else
  bign = strlen(str)
  n = strlen(sep)
  sep1 = substr(sep,1,1)
  isds1 = isdigit(sep1)
  if isds1
    str2 = repchar("a",n)~str
  else
    str2 = repchar("0",n)~str
  endif
  sepn = substr(sep,n,n)
  isdsn = isdigit(sepn)
  if isdsn
    str2 = str2~repchar("a",n)
  else
    str2 = str2~repchar("0",n)
  endif
  matrix posi = posstrstr(str2,sep)
  if posi[1] == 0
    strings rets
    rets+=str
  else
    matrix posi2 = posi .+ (n-1)
    strings rets
    np = nelem(posi)
    pos0 = n+1
    posn = n+bign
    posi0 = posi~(n+bign+1)
    posi20 = n~posi2
    loop i = 1..(np+1) -q
      minp = posi20[i]+1
      maxp = posi0[i]-1
      if minp<=maxp
        rets+=substr(str2,minp,maxp)
      endif
    endloop
  endif
endif
return rets

rx = rows(x)
cx = cols(x)
matrix y = zeros(rx,1)
matrix ret
if pos==1
  ret = y~x
elif pos>cx
  ret = x~y
else
  ret = x[,1:(pos-1)]~y~x[,pos:]
endif
return ret

ind.-specific variables

if ismucst
  X = 1~X
endif
scalar k = cols(X)
matrix b0 = zeros(k,1) ~ mshape(theta[1:(k*(n-1))],k,n-1)
matrix tmp = X*b0
matrix exbadj2 = exp(tmp)
return exbadj2./sumr(exbadj2)

return @fun

alldim = sum(didims)
nvars = nelem(didims)
if nvars > 1
  matrix cudidims = 1|(cum(didims')[1:(nvars-1)]+1)
else
  matrix cudidims = {1}
endif
matrix ret = {}
loop i = 1..nalts -q
  loop j = 1..nvars -q
    qurindj = cudidims[j] + didims[j]*(i-1)/nalts + seq(0,didims[j]/nalts-1)
    ret = ret~ qurindj
  endloop
endloop
return ret

ret = "1"
if n > 1
  loop i=2..n -q
    reti = sprintf(" %d",i)
    ret = ret~reti
  endloop
endif
ret = strstrip(ret)
return ret

$xlist in most cases forced ind. vars

list x = x-0
bignamss = varnames(x)
cknull = nelem(indlist)==0
if !cknull
  indnams = strjoin(varnames(indlist))
else
  indnams = ""
endif
string line
dili = ""
coli = ""
dilis = ""
colis = ""
string first
rets = ""
rets_left = ""
rets_right = ""
loop foreach i x -q
  cknull = 0
  if !cknull
    fififi = inlist(indlist,x[i])
    if fififi
      bignamsi = bignamss[i]
      rets = strstrip(rets~" "~bignamsi)
      colis = strstrip(colis~" "~bignamsi)
    else
      cknull = 1
    endif
  endif
  if cknull
    bundle bui = getinfo(x[i])
    vvv = bui.name
    uuu = bui.description
    ul = uniqvarlist(uuu)
    isspace = isstrstr(uuu," ")
    if isspace == 0
      isfrla = isfla(uuu)
      isfrla2 = vf(uuu)
      #print isfrla2
      if isfrla&&(isfrla2 > 0)
        #if isfrla
        ret3 = uuu
        flavars = uniqvarlist(uuu)
      else
        ret3 = vvv
        flavars = vvv
      endif
      nfla = nelem(strsplit(flavars))
      ret3s = parse_fla3(ret3,indnams)
      ret3 = ret3s[1]
      current_li = ret3s[2]
      current_di = ret3s[3]
      lili = strsub(current_li,"#","")
      didi = strsub(current_di,"#","")
      colis = colis~" "~lili
      dilis = dilis~" "~didi
      rets = rets~" "~ret3
    else
      iseq = substr(uuu,1,1)=="="
      if iseq
        issq = substr(vvv,1,3)=="sq_"
        sdpart = substr(vvv,4,strlen(vvv))
        fifi = findind(sdpart,ul)
        issq = issq*fifi>0
        islg = substr(vvv,1,2)=="l_"
        sdpart2 = substr(vvv,3,strlen(vvv))
        fifi2 = findind(sdpart2,ul)
        fifi20l = str_setinrersect(ul,allvns)
        new_condition = sdpart2==fifi20l
        islg = islg*fifi2*new_condition>0
        underscpos = catch_times(vvv, uuu)
        if issq
          ret3 = strstrip(sdpart~"^2")
          flavars = uniqvarlist(ret3)
          lili = flavars
          colis = colis~" "~lili
          rets = rets~" "~ret3
        elif islg
          ret3 = strstrip("log("~sdpart2~")")
          flavars = uniqvarlist(ret3)
          lili = flavars
          colis = colis~" "~lili
          rets = rets~" "~ret3
        elif  underscpos > 0
          nvvv = strlen(vvv)
          part1 = substr(vvv,1,underscpos-1)
          part2 = substr(vvv,underscpos+1,nvvv)
          ret3 = part1~"*"~part2
          rets = rets~" "~ret3
        else
          ret3 = vvv
          rets = rets~" "~ret3
        endif
        flavars = uniqvarlist(ret3)
        lili = flavars
        colis = colis~" "~lili
      else
        inti = substr(uuu,1,14)=="interaction of"
        inti2 = isstrstr(uuu,"==")
        intis = strsplit3(uuu,"==")
        nintis = nelem(intis)
        inti3 = (nintis==2)
        inti4 = 0
        if inti3
          intis2 = intis[2]
          nintis2 = strlen(intis2)
          cknin1 = substr(intis2,nintis2,nintis2)==")"
          cknin2 = isnat(strstrip(substr(intis2,1,nintis2-1)))
          inti4 = cknin1*cknin2
        endif
        inti = inti*inti2*inti3*inti4>0
        if inti
          zx = strsplit3(uuu,"and (")
          zx2 = zx[2]
          zx1 = zx[1]
          namis2 = strsplit3(zx2,"==")
          nami2 = strstrip(namis2[1])
          ni1 = strlen(zx1)
          nami1 = strstrip(substr(zx1,15,ni1))
          didi = getdiname(nami1)
          if strlen(didi) == 0
            first = nami1
          else
            first = null
          endif
          lili = strstrip(didi~" "~nami2)
          ret3 = doret2(uuu,vvv,lili,first)
          if !exists(first)
            first= ""
          endif
          rets = rets~" "~ret3
          dilis = dilis~" "~lili
          colis = colis~" "~first
        else
          doommy = isduname(vvv)
          if doommy
            ret3 = corr_duname(vvv)
            didi = getdiname(ret3)
            dilis = dilis~" "~didi
          else
            ret3 = vvv
            lili = ret3
            colis = colis~" "~lili
          endif
          rets = rets~" "~ret3
        endif
      endif
    endif
    colis = strstrip(colis)
    dilis = strstrip(dilis)
    rets_left = strstrip(rets_left~" "~vvv)
    rets_right = rets_right~" "~ret3
  endif
endloop
rets = strstrip(rets)
discrsi = uniqvarlist(dilis)
conti = uniqvarlist(colis)
return strstrip(rets)~";"~ strstrip(conti)~";"~strstrip(discrsi)

nams = strsplit(nam)
nn = nelem(nams)
nb = nelem(b)
nns = nn
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = mfun
loop i=1..nns -q
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
xbackup = x
funs2 = strsplit(fun2)
nf = nelem(funs2)
matrix ret = zeros(nn,nf)
matrix ret222 = zeros(1,nf)
nnxx = nelem(x)
matrix ret = zeros(nn,nb+1)
loop j = 1..nnxx -q
  x[j] = 0
  grads = zeros(2,nf)
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads[1,i] = @fun3
  endloop
  x[j] = 1
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads[2,i] = @fun3
  endloop
  if nb > nf
    grads = ones(2,1)~grads
  endif
  matrix mafx = grads*b
  if cmd == "lin"
    matrix logi0 = mafx
    matrix logi = mafx
    matrix ret1 = grads
  elif cmd == "probit"
    matrix logi0 = cnorm(mafx)
    matrix logi = dnorm(mafx)
    matrix ret1 = grads .*logi
  elif cmd == "logit"
    matrix exb = exp(mafx)
    matrix logi0 = exb ./(exb .+ 1)
    matrix logi = logi0 ./(exb .+ 1)
    matrix ret1 = grads .*logi
  elif cmd == "poisson"
    matrix logi0 = exp(mafx)
    matrix ret1 = grads .*logi0
  endif
  ret1 = ret1[2,]-ret1[1,]
  matrix retj = (logi0[2]-logi0[1])~ret1
  x = xbackup
  ret[j,] = retj
endloop
return ret

matrix b = vec(b)
nb = nelem(b)
nn = nelem(strsplit(nam))
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
string fun2 = mfun
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    fun2 = smartstrsub(fun2,znams[i],sprintf("z[%d]",i))
  endloop
endif
string fun223 = "{"~strsub(fun2," ",",")~"}"
if nns < bign
  matrix ret = fdjac(x,for_djac(x,fun223,z))'
else
  matrix ret = fdjac(x,for_djac(x,fun223))'
endif
nf = nelem(strsplit(fun2))
if cmd == "lin"
  if nb > nf
    ret = zeros(nn,1)~ret
  endif
  return (ret*b)~ret
elif cmd == "poisson"
  matrix ret222 = @fun223
  if nb > nf
    ret222 = 1~ret222
    ret = zeros(nn,1)~ret
  endif
  scalar xb = ret222*b
  logixb = exp(xb)
  matrix ret0 = ret*b
  matrix mafx = ret0.*logixb
  matrix ret4 = ret0*ret222
  ret4 = ret + ret4
  return mafx~ (ret4 .*logixb)
elif cmd == "logit"
  matrix ret222 = @fun223
  if nb > nf
    ret222 = 1~ret222
    ret = zeros(nn,1)~ret
  endif
  scalar xb = ret222*b
  matrix xb2 = ret222
  exb = exp(xb)
  logixb = exb/(exb+1)^2
  second = -exb*(exb-1)/((exb+1)^3)
  matrix ret2 = ret.*logixb
  matrix mafx = (ret*b).*logixb
  matrix ret0 = ret*b
  matrix ret4 = ret0*xb2
  ret4 = ret4.*second
  ret = (mafx~(ret4+ret2))
  return ret
elif cmd == "probit"
  matrix ret222 = @fun223
  if nb > nf
    ret222 = 1~ret222
    ret = zeros(nn,1)~ret
  endif
  scalar xb = ret222*b
  matrix xb2 = ret222
  logixb = dnorm(xb)
  second = -logixb*xb
  matrix ret2 = ret.*logixb
  matrix mafx = (ret*b).*logixb
  matrix ret0 = ret*b
  matrix ret4 = ret0*xb2
  ret4 = ret4.*second
  ret = (mafx~(ret4+ret2))
  return ret
endif

strings  factlevss = strsplit(factlevs)
nns = nelem(factlevss)
strings nams = array(nns)
string fun2 = mfun
strings allnamess = strsplit(allnam)
scalar bign = nelem(allnamess)
loop i=1..bign -q
  fun2 = smartstrsub(fun2,allnamess[i],sprintf("x[%d]",i))
endloop
matrix indi = zeros(nns,1)
loop i = 1.. nns -q
  nams[i] = nam~"_"~factlevss[i]
  indi[i] = findind(nams[i],allnam)
endloop
funs2 = strsplit(fun2)
nf = nelem(funs2)
nb = nelem(b)
matrix x = allval
bignns = nelem(strsplit(alllevs))
if bignns > nns
  x[indi] = zeros(1,nns)
  matrix grads2 = zeros(nns+1,nf)
  loop i=1..nf -q
    string fun3 = funs2[i]
    grads2[1,i] = @fun3
  endloop
  loop j = 1..nns -q
    x[indi[j]] = 1
    loop i=1..nf -q
      string fun3 = funs2[i]
      grads2[j+1,i] = @fun3
    endloop
    x[indi] = zeros(1,nns)
  endloop
  if nb > nf
    grads2=ones(nns+1,1)~grads2
  endif
else
  matrix grads2 = zeros(nns,nf)
  x[indi] = zeros(1,nns)
  loop j = 1..nns -q
    x[indi[j]] = 1
    loop i=1..nf -q
      string fun3 = funs2[i]
      grads2[j,i] = @fun3
    endloop
    x[indi] = zeros(1,nns)
  endloop
  if nb > nf
    grads2=ones(nns,1)~grads2
  endif
endif
matrix mafx = grads2*b
if cmd == "logit"
  matrix exb = exp(mafx)
  matrix logi0 = exb ./(exb .+ 1)
  matrix logi = logi0 ./(exb .+ 1)
  matrix ret10 = grads2 .*logi
  ret10 = ret10[2:,] .- ret10[1,]
  matrix ret00 = (logi0[2:] .- logi0[1])~ret10
  return ret00
elif cmd == "poisson"
  matrix logi0 = exp(mafx)
  matrix logi = logi0
  matrix ret10 = grads2 .*logi
  ret10 = ret10[2:,] .- ret10[1,]
  matrix ret00 = (logi0[2:] .- logi0[1])~ret10
  return ret00
elif cmd == "probit"
  matrix logi0 = cnorm(mafx)
  matrix logi = dnorm(mafx)
  matrix ret10 = grads2 .*logi
  ret10 = ret10[2:,] .- ret10[1,]
  matrix ret00 = (logi0[2:] .- logi0[1])~ret10
  return ret00
elif cmd == "lin"
  matrix mafx = grads2*b
  matrix logi0 = mafx
  matrix logi = logi0
  matrix ret10 = grads2
  ret10 = ret10[2:,] .- ret10[1,]
  matrix ret00 = (logi0[2:] .- logi0[1])~ret10
  return ret00
endif

nams = strsplit(nam)
nn = nelem(nams)
nb = nelem(b)
matrix bigb = vec(b)
matrix b = bigb[1:(nb-nalts+1)]
matrix cuts = bigb[(nb-nalts+2):]
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = mfun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
endif
fun2 = "myx("~strsub(fun2," ",") myx(")~")"
funs2 = strsplit(fun2)
nf = nelem(funs2)
matrix ret22 = zeros(nn,nf)
matrix ret222 = zeros(1,nf)
loop i=1..nf -q
  string fun3 = funs2[i]
  ret22[,i] = fdjac(x,@fun3)'
  ret222[i]=@fun3
endloop
matrix xb = ret222*b
matrix xbm = cuts.-xb
if cmd == "logit"
  matrix exbm = exp(xbm)
  matrix logixbm = exbm ./(exbm .+1).^2
  matrix  logixbm1 = 0|logixbm
  matrix  logixbm2 = logixbm|0
  logixbm = logixbm1-logixbm2
  matrix mfxm = vec(logixbm*(ret22*b)')
  ggg = ret222.*ones(nalts-1,1)
  matrix xb2m = ggg~(-I(nalts-1))
  xb100 = xb2m*bigb
  exb100 = exp(xb100)
  logixb100 = exb100./(exb100.+1).^2
  second100 = -exb100.*(exb100.-1)./((exb100.+1).^3)
  ret10 = ret22~zeros(nns,nalts-1)
  matrix ret200 = ret10**logixb100
  matrix ret00 = ret22*b
  matrix ret400 = ret00**xb2m
  ret400 = ret400.*mshape(second100,nns*nelem(second100),1)
  ret600 = ret200+ret400
  ret601 = zeros(nalts*nns,nb)
  ret602 = ret601
  indind = vec(seq(1,nalts-1)'.+ (seq(0,nns-1)*nalts))
  ret601[indind,] = ret600
  ret602[indind+1,] = ret600
  return  mfxm ~(ret601- ret602)
else
  matrix logixbm = dnorm(xbm)
  matrix  logixbm1 = 0|logixbm
  matrix  logixbm2 = logixbm|0
  logixbm = logixbm1-logixbm2
  matrix mfxm = vec(logixbm*(ret22*b)')
  ggg = ret222.*ones(nalts-1,1)
  matrix xb2m = ggg~(-I(nalts-1))
  xb100 = xb2m*bigb
  exb100 = exp(xb100)
  logixb100 = dnorm(xb100)
  second100 = -logixb100.*xb100
  ret10 = ret22~zeros(nns,nalts-1)
  matrix ret200 = ret10**logixb100
  matrix ret00 = ret22*b
  matrix ret400 = ret00**xb2m
  ret400 = ret400.*mshape(second100,nns*nelem(second100),1)
  ret600 = ret200+ret400
  ret601 = zeros(nalts*nns,nb)
  ret602 = ret601
  indind = vec(seq(1,nalts-1)'.+ (seq(0,nns-1)*nalts))
  ret601[indind,] = ret600
  ret602[indind+1,] = ret600
  return  mfxm ~(ret601- ret602)
endif

nams = strsplit(nam)
nn = nelem(nams)
nb = nelem(b)
matrix ret = zeros(nn*nalts,nb+1)
matrix bigb = vec(b)
matrix b = bigb[1:(nb-nalts+1)]
matrix cuts = bigb[(nb-nalts+2):]
strings allnams = strsplit(allnam)
scalar nlas = nelem(allnams)
string fun2 = mfun
loop i=1..nlas -q
  fun2 = smartstrsub(fun2,allnams[i],sprintf("x[%d]",i))
endloop
strings fun2s = strsplit(fun2)
seqna = seq(1,nalts)
loop i=1..nn -q
  indi20 = nalts*(i-1).+ seqna
  ret[indi20,] = olp_dum_grad2(fun2,fun2s,mfun, allnam, &allval, nams[i], &b,&bigb,&cuts,nb, nalts, cmd)
endloop
return ret

scalar indi = findind(nam,allnam)
matrix Z = zeros(2,nelem(allval))
matrix allvi = allval
allvi[indi]=0
Z[1,]= allvi
allvi = allval
allvi[indi]=1
Z[2,]= allvi
rz = 2
scalar nf = nelem(fun2s)
matrix grads = zeros(rz*(nalts-1),nb)
Im = I(nalts-1)
kk = nb-nalts+1
loop i=1..rz -q
  matrix x = Z[i,]
  matrix retX = zeros(1,nf)
  loop iii=1..nf -q
    fun2si = fun2s[iii]
    retX[iii]=@fun2si
  endloop
  jj = seq(1,nalts-1)
  ll = rz*(jj-1)+i
  grads[ll,1:kk] = {1;1} .* retX
  grads[ll,kk+1:nb] = -Im[jj,]
endloop
matrix mfxx = grads*bigb
if cmd == "probit"
  matrix exbx = cnorm(mfxx)
  matrix exbxx =dnorm(mfxx)
  matrix ret300 = grads.*exbxx
  matrix ret400 = zeros(nalts-1,nb)
  matrix mfx0=zeros(nalts-1,1)
  loop i=1..(nalts-1) -q
    mfx0[i]=exbx[2*i-1]-exbx[2*i]
    ret400[i,]=ret300[2*i-1,]-ret300[2*i,]
  endloop
  mfx1 = 0|mfx0
  mfx2 = mfx0|0
  ret11 = mfx2-mfx1
  ret4001 = 0|ret400
  ret4002 = ret400|0
  ret21 = ret4001-ret4002
  return ret11~ret21
else
  matrix exbx0 = exp(mfxx)
  matrix exbx = exbx0./(exbx0.+1)
  matrix exbxx =exbx./(exbx0.+1)
  matrix ret300 = grads.*exbxx
  matrix ret400 = zeros(nalts-1,nb)
  matrix mfx0=zeros(nalts-1,1)
  loop i=1..(nalts-1) -q
    mfx0[i]=exbx[2*i-1]-exbx[2*i]
    ret400[i,]=ret300[2*i-1,]-ret300[2*i,]
  endloop
  mfx1 = 0|mfx0
  mfx2 = mfx0|0
  ret11 = mfx2-mfx1
  ret4001 = 0|ret400
  ret4002 = ret400|0
  ret21 = ret4001-ret4002
  return ret11~ret21
endif

matrix bigb = vec(b)
# sorted by levels of discrete variable
nb = nelem(bigb)
matrix b = bigb[1:(nb-nalts+1)]
matrix cuts = bigb[(nb-nalts+2):]
string fun2 = mfun
strings allnams = strsplit(allnam)
scalar nlas = nelem(allnams)
loop i=1..nlas -q
  fun2 = smartstrsub(fun2,allnams[i],sprintf("x[%d]",i))
endloop
strings fun2s = strsplit(fun2)
scalar nf = nelem(fun2s)
matrix Z = substvals2(allnam,allval,nam,factlevs,alllevs )
rz = rows(Z)
matrix grads = zeros(rz*(nalts-1),nb)
matrix Im = I(nalts-1)
kk = nb-nalts+1
loop i=1..rz -q
  matrix x = Z[i,]
  matrix retX = zeros(1,nf)
  loop iii=1..nf -q
    fun2si = fun2s[iii]
    retX[iii]=@fun2si
  endloop
  jj = seq(1,nalts-1)
  ll = rz*(jj-1)+i
  grads[ll,1:kk] = {1;1} .* retX
  grads[ll,kk+1:nb] = -Im[jj,]
endloop
matrix mfxx = grads*bigb
if cmd == "logit"
  matrix exbx0 = exp(mfxx)
  matrix exbx = exbx0./(exbx0.+1)
  matrix exbxx =exbx./(exbx0.+1)
else
  matrix exbx = cnorm(mfxx)
  matrix exbxx =dnorm(mfxx)
endif
matrix ret300 = grads.*exbxx
matrix ret400 = zeros(nalts-1,nb)
matrix Zi = zeros(rz,1)
matrix indi = Zi
matrix indi2 = zeros(rz-1,1)
matrix Zis = zeros((rz-1)*(nalts-1),1)
matrix ret400 = zeros((rz-1)*(nalts-1),nb)
matrix ret400i = {}
loop i=1..(nalts-1) -q
  indi = rz*(i-1).+seq(1,rz)
  Zi = exbx[indi]
  ret400i = ret300[indi,]
  Zi=Zi[2:,].-Zi[1,]
  ret400i = ret400i[2:,].-ret400i[1,]
  indi2 = i .+ (nalts-1).*seq(0,rz-2)
  Zis[indi2]=Zi
  ret400[indi2,]=ret400i
endloop
matrix mfal = {}
matrix zindi = {}
matrix mret400 = {}
matrix m400i = {}
matrix zi ={}
loop i=1..(rz-1) -q
  zindi =(i-1)*(nalts-1).+ seq(1,nalts-1)
  zi = Zis[zindi]
  m400i = ret400[zindi,]
  mfx1 = 0|zi
  mfx2 = zi|0
  m400i1 = zeros(1,nb)|m400i
  m400i2 = m400i|zeros(1,nb)
  mfal=mfal|(mfx1-mfx2)
  matrix mret400 = mret400|(m400i1-m400i2)
endloop
return mfal~mret400

n = nalts
matrix X = svals
if ismuc
  X = 1~X
endif
scalar k = cols(X)
matrix b0 = mshape(b[1:(k*(n-1))],k,n-1)
b0 = zeros(k,1) ~ b0
matrix tmp = X*b0
matrix exbadj2 = exp(tmp)
matrix pij = exbadj2./sumr(exbadj2)
betam = grads0*b0
alphabar = betam * pij'
mfixes = pij.*(betam .-alphabar)
return vec(mfixes)

xb = "("~xb~")"
cey = "(xb__ + s__*(dnorm((lo__-xb__)/s__)-dnorm((hi__-xb__)/s__))/ (cnorm((hi__-xb__)/s__)-cnorm((lo__-xb__)/s__)))"
cey = smartstrsub(cey,"xb__",xb)
ey = "lo__*cnorm((lo__-xb__)/s__)+hi__*cnorm(-(hi__-xb__)/s__)+(cnorm((hi__-xb__)/s__)-cnorm((lo__-xb__)/s__))*cey__"
ey = smartstrsub(ey,"xb__",xb)
ey = smartstrsub(ey,"cey__",cey)
return ey

if missing(lo__)&&missing(hi__)
  ret = ynam
elif !missing(lo__)
  if missing(hi__)
    ret = sprintf("%s > %g",ynam,lo__)
  else
    ret = sprintf("%g < %s < %g",lo__,ynam,hi__)
  endif
else
  ret = sprintf("%s < %g",ynam,hi)
endif
return ret

xb = "("~xb~")"
Phi1 = "cnorm((lo__-xb__)/s__)"
Phi2 = "cnorm((hi__-xb__)/s__)"
phi1 = "dnorm((lo__-xb__)/s__)"
phi2 = "dnorm((hi__-xb__)/s__)"
num = "("~phi1~"-"~phi2~")"
den = "("~Phi2~"-"~Phi1~")"
ratio = "s__*"~num~"/"~den
ratio = smartstrsub(ratio,"xb__",xb)
ret = xb~"+"~ratio
return ret

funs = strsplit(fun)
nf = nelem(funs)
ret = ""
if nf == n
  loop i = 1..n -q
    term = sprintf("b[%d]*%s", i , funs[i])
    #print term
    if i == 1
      ret = ret~term
    else
      ret = ret~"+"~term
    endif
  endloop
elif n == (nf + 1)
  ret = "b[1]"
  loop i = 2..n -q
    term = sprintf("b[%d]*%s", i , funs[i-1])
    ret = ret~"+"~term
  endloop
else
  funcerr "Wrong lengths of b"
endif
return ret

funs = strsplit(fun)
nf = nelem(funs)
ret = ""
if nf == n
  loop i = 1..n -q
    j = i+nx
    term = sprintf("x[%d]*%s", j , funs[i])
    #print term
    if i == 1
      ret = ret~term
    else
      ret = ret~"+"~term
    endif
  endloop
elif n == (nf + 1)
  j = nx + 1
  ret = sprintf("x[%d]", j)
  loop i = 2..n -q
    j = i+nx
    term = sprintf("x[%d]*%s", j , funs[i-1])
    ret = ret~"+"~term
  endloop
else
  funcerr "Wrong lengths of b"
endif
return ret

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
endif
fun2bkp = fun2
fun2 = text_xb(fun2,nelem(b))
fun200 = text_xb3(fun2bkp,nelem(x),nelem(b))
fun2 = tobit_text(fun2)
fun200 = tobit_text(fun200)
string fun3 = "log("~fun2~")"
string fun300 = "log("~fun200~")"
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
matrix valval = @fun2
matrix ret = fdjac(x,@fun3)
nz = nelem(z)
sub1 = sprintf("z[%d]",(nz+1))
sub2 = sprintf("z[%d]",(nz+2))
z = z~lo__~hi__
nb = nelem(b)
b = b|s__
nx = nelem(x)
sub400 = sprintf("x[%d]",(nx+nb+1))
fun400 = smartstrsub(fun300,"lo__",sub1)
fun400 = smartstrsub(fun400,"hi__",sub2)
fun400 = smartstrsub(fun400,"s__",sub400)
xbkp = x
ret = (ret.*x)'
matrices M = array(4)
M[1] = valval
M[2] = ret
x = x~b'
matrix hess0 = numhess(x,@fun400)[1:nx,(nx+1):]
#print hess0
matrix hess = numhess(x,@fun400)[1:nx,(nx+1):].*xbkp'
nnn = cols(vcv)
matrix Jm = I(nnn)
Jm[nnn,nnn] = s__
print s__
print vcv
matrix Vbig = qform(Jm, vcv)
matrix mvc = qform(hess,Vbig)
M[3] = mvc
M[4] = sqrt(diag(mvc))
return M

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
endif
fun2bkp = fun2
fun2 = text_xb(fun2,nelem(b))
fun200 = text_xb3(fun2bkp,nelem(x),nelem(b))
fun2 = tobit_text2(fun2)
fun200 = tobit_text2(fun200)
string fun3 = "log("~fun2~")"
string fun300 = "log("~fun200~")"
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
matrix valval = @fun2
matrix ret = fdjac(x,@fun3)
ret = (ret.*x)'
nz = nelem(z)
sub1 = sprintf("z[%d]",(nz+1))
sub2 = sprintf("z[%d]",(nz+2))
z = z~lo__~hi__
nb = nelem(b)
b = b|s__
nx = nelem(x)
sub400 = sprintf("x[%d]",(nx+nb+1))
fun400 = smartstrsub(fun300,"lo__",sub1)
fun400 = smartstrsub(fun400,"hi__",sub2)
fun400 = smartstrsub(fun400,"s__",sub400)
xbkp = x
matrices M = array(4)
M[1] = valval
M[2] = ret
x = x~b'

matrix hess = numhess(x,@fun400)[1:nx,(nx+1):].*xbkp'
nnn = cols(vcv)
matrix Jm = I(nnn)
Jm[nnn,nnn] = s__
matrix Vbig = qform(Jm, vcv)
matrix mvc = qform(hess,Vbig)
M[3] = mvc
M[4] = sqrt(diag(mvc))
return M

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
endif
fun2bkp = fun2
fun2 = text_xb(fun2,nelem(b))
fun200 = text_xb3(fun2bkp,nelem(x),nelem(b))
fun2 = tobit_text2(fun2)
fun200 = tobit_text2(fun200)
string fun3 = "("~fun2~")"
string fun300 = "("~fun200~")"
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
matrix valval = @fun2
matrix ret = fdjac(x,@fun3)'
#ret = (ret.*x)'
nz = nelem(z)
sub1 = sprintf("z[%d]",(nz+1))
sub2 = sprintf("z[%d]",(nz+2))
z = z~lo__~hi__
nb = nelem(b)
b = b|s__
nx = nelem(x)
sub400 = sprintf("x[%d]",(nx+nb+1))
fun400 = smartstrsub(fun300,"lo__",sub1)
fun400 = smartstrsub(fun400,"hi__",sub2)
fun400 = smartstrsub(fun400,"s__",sub400)
xbkp = x
x = x~b'
matrix hess = numhess(x,@fun400)[1:nx,(nx+1):]
return ret~hess

if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
smallb = b
b = b|s__
nams = strsplit(nam)
nn = nelem(nams)
nb = nelem(b)
nns = nn
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
loop i=1..nns -q
  fun2 = smartstrsub(fun2,nams[i],sprintf("x[%d]",i))
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z = {}
endif
nz = nelem(z)
sub1 = sprintf("z[%d]",(nz+1))
sub2 = sprintf("z[%d]",(nz+2))
z = z~lo__~hi__
matrix xbackup = x
fun2a = tobit_text2(text_xb(fun2,nb-1))
funs2 = strsplit(fun2)
nf = nelem(funs2)
matrix ret222 = zeros(1,nf)
nnxx = nelem(x)
matrix ret1 = zeros(nn,nb)
matrix ret2 = ret1
matrix retZ1 = zeros(nn,1)
matrix retZ2 = retZ1
sub4 = sprintf("b[%d]",(nb))
fun4 = smartstrsub(fun2a,"lo__",sub1)
fun4 = smartstrsub(fun4,"hi__",sub2)
fun4 = smartstrsub(fun4,"s__",sub4)
fun4 = "("~fun4~")"
loop j = 1..nnxx -q
  x[j] = 0
  retZ1[j] = @fun2a
  ret1[j,] = fdjac(b,@fun4)
  x[j] = 1
  retZ2[j] = @fun2a
  ret2[j,] = fdjac(b,@fun4)
  x = xbackup
endloop
return (retZ2 - retZ1)~(ret2 - ret1)

matrix b = vec(b)
if missing(hi__)
  hi__ = $huge
endif
if missing(lo__)
  lo__ = -$huge
endif
matrix Z = substvals2(allnam,allval,nam,factlevs,alllevs )
rz = rows(Z)
n = nelem(b)
fun3 = text_xb(mfun,n)
fun4 = tobit_text2(fun3)
sub4 = sprintf("b[%d]",(n+1))
fun4 = smartstrsub(fun4,"s__",sub4)
sub1 = sprintf("z[%d]",1)
sub2 = sprintf("z[%d]",2)
fun4 = smartstrsub(fun4,"lo__",sub1)
fun4 = smartstrsub(fun4,"hi__",sub2)
b = b|s__
matrix z = {lo__,hi__}
# print fun4
matrix ret3 = zeros(rz,n+2)
loop i=1..rz -q
  ret3[i,] = t_txtfun(fun4,allnam,Z[i,],z,b)
endloop
ret3 = ret3[2:,] .- ret3[1,]
#return ret3
return ret3

strings allnams = strsplit(allnam)
scalar n = nelem(allnams)
string fun2 = fun
string txti
loop i=1..n -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,allnams[i],txti)
endloop
#print fun2
matrix x = allval
fun3 = "myx("~fun2~")"
matrix ret2 = fdjac(b,@fun3)
matrix ret = @fun2~ret2
return ret

ccut = sprintf("%s[%d]",lette,nume)
if type == 1
  ret = sprintf("cnorm(%s-(%s))",ccut,xb)
else
  ret = sprintf("1/(1+exp(-(%s-(%s))))",ccut,xb)
endif
return ret

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
endif
fun2bkp = fun2
fun2 = text_xb(fun2,numofx)
fun200 = text_xb3(fun2bkp,nns,numofx)
ncuts = nelem(b) - numofx
n_alts = ncuts + 1
ret0100 = zeros(n_alts,nns)
xbkp = x
hessis0100 = zeros(nns*n_alts,nelem(b))
matrix fun2ss = zeros(ncuts,1)
loop i=1..n_alts  -q
  cuf =  cuti(i+numofx,fun2,"b",type)
  cutf = cuti(i+numofx+nns,fun200,"x",type)
  if i == 1
    cuf = "log("~cuf~")"
    cutf = "log("~cutf~")"
  elif i < n_alts
    cuf = "log("~cuf~"-("~cuti(i-1+numofx,fun2,"b",type)~"))"
    cutf = "log("~cutf~"-("~cuti(i-1+numofx+nns,fun200,"x",type)~"))"
  else
    cuf = "log(1-("~cuti(i-1+numofx,fun2,"b",type)~"))"
    cutf = "log(1-("~cuti(i-1+numofx+nns,fun200,"x",type)~"))"
  endif
  ret0100[i,] = fdjac(x,@cuf).*x
  x = x~b'
  hrows = (i-1)*nns.+seq(1,nns)
  hessis0100[hrows,] = numhess(x,@cutf)[1:nns,(nns+1):].*xbkp'
  x = xbkp
endloop
matrix mama = qform(hessis0100,vcv)
matrix ses = sqrt(diag(mama))
matrices M = array(5)
M[1] = fun2ss
M[2] = ret0100
M[3] = mama
M[4] = mshape(ses,nns,n_alts)'
M[5] = {n_alts}
return M

funs = strsplit(fun)
nf = nelem(funs)
ret = ""
if nf == n
  loop i = 1..n -q
    term = sprintf("b[%d]*%s", i+start-1 , funs[i])
    if i == 1
      ret = ret~term
    else
      ret = ret~"+"~term
    endif
  endloop
elif n == (nf + 1)
  ret = sprintf("b[%d]" , start)
  loop i = 2..n -q
    term = sprintf("b[%d]*%s", i+start-1 , funs[i-1])
    ret = ret~"+"~term
  endloop
else
  funcerr "Wrong lengths of b"
endif
return ret

smalln = bign/(nalts-1)
ret0 = ml_text_sum(nalts-1)
loop i = 1..(nalts-1) -q
  find = sprintf("#%d",i)
  subst = ml_text_xb(fun,smalln,smalln*(i-1)+1)
  ret0 = strsub(ret0,find,subst)
endloop
return "1+"~ret0

ret = "exp(#1)"
if k > 1
  loop i = 2..k -q
    reti = sprintf("exp(#%d)", i)
    ret = ret~"+"~reti
  endloop
endif
return ret

den =  ml_text_den(fun,bign,nalts)
num = "1"
if alt > 1
  smalln = bign/(nalts-1)
  num = sprintf("exp(#%d)", alt)
  find = sprintf("#%d",alt)
  subst = ml_text_xb(fun,smalln,smalln*(alt-2)+1)
  num = strsub(num,find,subst)
endif
ret = num~"/("~den~")"
return ret

funre = fun
loop i=1..nb -q
  find = sprintf("b[%d]",i)
  repl = sprintf("x[%d]",i+nx)
  funre = smartstrsub(funre, find, repl)
endloop
return funre

matrix b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
matrix indi = zeros(nns,1)
loop i=1..nns -q
  indi[i] = findind(nams[i],allnam)
endloop
string fun2 = fun
string txti
loop i=1..nns -q
  txti = sprintf("x[%d]",i)
  fun2 = smartstrsub(fun2,nams[i],txti)
endloop
matrix x = allval[indi]
scalar bign = nelem(strsplit(allnam))
if nns < bign
  matrix z = mexclude(allval,indi)
  string znam = str_setdiff(allnam,nam)
  scalar nnz = bign - nns
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  z ={}
endif
matrix ret = zeros(nalts,nns)
hessis0100 = zeros(nns*nalts,nelem(b))
xbkp = x
nnbb = nelem(b)
matrix valval = zeros(nalts,1)
loop i=1..nalts -q
  fun3 = ml_text(fun2,nelem(b),nalts,i)
  valval[i] = @fun3
  fun3 = "log("~fun3~")"
  ret[i,] = fdjac(x,@fun3).*x
  x = x~b'
  hrows = (i-1)*nns.+seq(1,nns)
  fun4 = ml_text_xx(fun3,nnbb,nns)
  hessis0100[hrows,] = numhess(x,@fun4)[1:nns,(nns+1):].*xbkp'
  x = xbkp
endloop
matrix mama = qform(hessis0100,vcv)
matrix ses = sqrt(diag(mama))
matrices M = array(5)
M[1] = valval
M[2] = ret
M[3] = mama
M[4] = mshape(ses,nns,nalts)'
M[5] = {nalts}
return M

strs = strsplit3(str,";")
n = nelem(strs)
rty = ""
rw = ""
rat = ""
if n > 0
  loop i = 1..n -q
    stri = strs[i]
    stri2 = strstrip(strsub(stri," ",""))
    ckt = substr(stri2,1,5) == "type="
    ckw = substr(stri2,1,4) == "wrt="
    cka = substr(stri2,1,3) == "at="
    if ckw
      rw = strstrip(strsplit3(stri,"=")[2])
    elif cka
      rat = strstrip(strsplit3(stri,"=")[2])
    elif ckt
      rty =  strstrip(strsplit3(stri,"=")[2])
    endif
  endloop
endif
strings rets = array(3)
rets[1] = rw
rets[2] = rat
rets[3] = rty
return rets

b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
if nns > 0
  matrix indi = zeros(nns,1)
  loop i=1..nns -q
    indi[i] = findind(nams[i],allnam)
  endloop
  string fun2 = fun
  string txti
  if nns > 0
    loop i=1..nns -q
      txti = sprintf("x[%d]",i)
      fun2 = smartstrsub(fun2,nams[i],txti)
    endloop
  endif
  matrix x = allval[indi]
else
  string fun2 = fun
  matrix x = {}
  matrix indi = {}
endif
scalar bign = nelem(strsplit(allnam))
scalar nnz = bign - nns
matrix z = {}
if (nns < bign)&&(bign > 0)
  if nelem(indi) == 0
    matrix z = allval
  else
    matrix z = mexclude(allval,indi)
  endif
  string znam = str_setdiff(allnam,nam)
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
  string znam = ""
endif
fun2bkp = fun2
fun2 = text_xb(fun2,n1-1)
fun2s = strsplit(fun2bkp)
isscst = (n1-1) > nelem(fun2s)
fun321 = xtextfun(fun2s,1,isscst,"b")
fun321bkp = fun321
maxX = nns
maxZ = nnz
ZNAM = znam
strings nams = strsplit(namZ)
scalar nns = nelem(nams)
matrix indiz = {}
if nns > 0
  matrix indiz = zeros(nns,1)
  loop i=1..nns -q
    indiz[i] = findind(nams[i],allnamZ)
  endloop
endif
string fun2 = funZ
commonX = str_setinrersect(namZ,nam)
commonXs = strsplit(commonX)
ncX = nelem(commonXs)
aloneX = str_setdiff(namZ,nam)
aloneXs = strsplit(aloneX)
naX = nelem(aloneXs)
ALLX = maxX + naX
if ncX > 0
  coindi = findinds(commonX,nam)
  loop i=1..ncX -q
    txti = sprintf("x[%d]",coindi[i])
    fun2 = smartstrsub(fun2,commonXs[i],txti)
  endloop
endif
matrix aindiX = {}
matrix xZZ = {}
matrix aindiZ = {}
matrix zZZ = {}
if naX > 0
  matrix aindiX = findinds(aloneX,namZ)
  matrix xZZ = allvalZ[aindiX]
  loop i=1..naX -q
    txti = sprintf("x[%d]",i+maxX)
    fun2 = smartstrsub(fun2,aloneXs[i],txti)
  endloop
endif
scalar bign = nelem(strsplit(allnamZ))
scalar nnz = bign - nns
if nns < bign
  string znamZ = str_setdiff(allnamZ,namZ)
  strings znamZs = strsplit(znamZ)
  if maxZ == 0
    matrix aindiZ = findinds(znamZ,allnamZ)
    matrix zZZ = allvalZ[aindiZ]
    loop i=1..nnz -q
      txti = sprintf("z[%d]",i)
      fun2 = smartstrsub(fun2,znamZs[i],txti)
    endloop
  else
    commonZ = str_setinrersect(znamZ,ZNAM)
    commonZs = strsplit(commonZ)
    aloneZ = str_setdiff(znamZ,ZNAM)
    aloneZs = strsplit(aloneZ)
    ncZ = nelem(commonZs)
    naZ = nelem(aloneZs)
    if ncZ > 0
      coindiZ = findinds(commonZ,ZNAM)
      loop i=1..ncZ -q
        txti = sprintf("z[%d]",coindiZ[i])
        fun2 = smartstrsub(fun2,commonZs[i],txti)
      endloop
    endif
    ALLZ = maxZ + naZ
    if naZ > 0
      matrix aindiZ = findinds(aloneZ,allnamZ)
      matrix zZZ = allvalZ[aindiZ]
      loop i=1..naZ -q
        txti = sprintf("z[%d]",i+maxZ)
        fun2 = smartstrsub(fun2,aloneZs[i],txti)
      endloop
    endif
  endif
endif
fun2bkpZ = fun2
fun2s = strsplit(fun2bkpZ)
isscst = (n2-n1) > nelem(fun2s)
fun321Z = xtextfun(fun2s,n1+1,isscst,"b")
fun321Zbkp = fun321Z
x = x~xZZ
z = z~zZZ
matrix xBACKUP = x
#print link
if link != "Ey"
  fun321Z = "invmills(-("~fun321Z~"))"
  bzzzz = sprintf("b[%d]",n1)
  fun321Z = bzzzz~"*"~fun321Z
  fun421 = "("~fun321~"+"~fun321Z~")"
  if etype == 1
    dydx = fdjac(x,@fun421)'
    funHe = fun421
    loop i=1..n2 -q
      finD = sprintf("b[%d]",i)
      repL = sprintf("x[%d]",i+ALLX)
      funHe = strsub(funHe,finD,repL)
    endloop
    x = x~b'
    hess =  numhess(x,@funHe)[1:ALLX,(ALLX+1):]
    return dydx~hess
    # endif
    #matrix xBACKUP = x
  elif etype == 2
    retd = zeros(ALLX,n2)
    matrix retdmf = zeros(ALLX,1)
    matrix x = xBACKUP
    loop i = 1..ALLX -q
      x[i] = 0
      nfun1 = @fun421
      ndy1 = fdjac(b,@fun421)
      #print ndy1
      x[i] = 1
      retdmf[i] = @fun421 - nfun1
      retd[i,] = fdjac(b,@fun421) - ndy1
      x = xBACKUP
    endloop
    return retdmf~retd
  endif
else
  fun521Z = "dnorm(("~fun321Zbkp~"))"
  fun521 = "cnorm("~fun321Zbkp~")"
  fun521 = fun521~"*("~fun321~")"
  bzzzz = sprintf("b[%d]",n1)
  fun521Z = bzzzz~"*"~fun521Z
  fun721 = "("~fun521~"+"~fun521Z~")"
  if etype == 1
    dydx = fdjac(x,@fun721)'
    #eval fdjac(x,@fun721)'
    funHe = fun721
    loop i=1..n2 -q
      finD = sprintf("b[%d]",i)
      repL = sprintf("x[%d]",i+ALLX)
      funHe = strsub(funHe,finD,repL)
    endloop
    x = x~b'
    hess =  numhess(x,@funHe)[1:ALLX,(ALLX+1):]
    return dydx~hess
  else
    retd = zeros(ALLX,n2)
    matrix retdmf = zeros(ALLX,1)
    matrix x = xBACKUP
    loop i = 1..ALLX -q
      x[i] = 0
      nfun1 = @fun721
      ndy1 = fdjac(b,@fun721)
      #print ndy1
      x[i] = 1
      retdmf[i] = @fun721 - nfun1
      retd[i,] = fdjac(b,@fun721) - ndy1
      x = xBACKUP
    endloop
    return retdmf~retd
  endif
endif
#return dydx~hess

b = vec(b)
fun11 = xtextfun(strsplit(fla),1,isc1,"b")
fun12 = xtextfun(strsplit(flaZ),n1+2,isc2,"b")
matrix Z = substvals2(allnam,allval,nam,factlevs,alllevs )
allnams = strsplit(allnam)
nna = nelem(allnams)
if link != "Ey"
  fun12 = sprintf("b[%d]*invmills(-(%s))", n1+1,fun12)
  fun3 = sprintf("(%s+%s)",fun11,fun12)
  fun4 = fun3
else
  fun12bkp = fun12
  fun12 = sprintf("b[%d]*dnorm(%s)", n1+1,fun12)
  fun12a = sprintf("cnorm(%s)*(%s)",fun12bkp,fun11)
  fun3 = sprintf("(%s+%s)",fun12a,fun12)
  fun4 = fun3
endif
loop i = 1..nna -q
  findi = allnams[i]
  substi = sprintf("x[%d]",i)
  fun4 = smartstrsub(fun4,findi,substi)
endloop
rz = rows(Z)
matrix ret0 = zeros(rz,1)
matrix ret20 = zeros(rz,nelem(b))
matrix x = allval
matrix xbkp = x
loop i = 1..rz -q
  x = Z[i,]
  ret20[i,]= fdjac(b,@fun4)
  ret0[i] = @fun4
  x = xbkp
endloop
matrix ret01 = (ret0 .-ret0[1])[2:]
matrix ret200 = (ret20 .- ret20[1,])[2:,]
return ret01~ret200

b = vec(b)
strings nams = strsplit(nam)
scalar nns = nelem(nams)
if nns > 0
  matrix indi = zeros(nns,1)
  loop i=1..nns -q
    indi[i] = findind(nams[i],allnam)
  endloop
  string fun2 = fun
  string txti
  if nns > 0
    loop i=1..nns -q
      txti = sprintf("x[%d]",i)
      fun2 = smartstrsub(fun2,nams[i],txti)
    endloop
  endif
  matrix x = allval[indi]
else
  string fun2 = fun
  matrix x = {}
  matrix indi = {}
endif
scalar bign = nelem(strsplit(allnam))
scalar nnz = bign - nns
matrix z = {}
if (nns < bign)&&(bign > 0)
  if nelem(indi) == 0
    matrix z = allval
  else
    matrix z = mexclude(allval,indi)
  endif
  string znam = str_setdiff(allnam,nam)
  strings znams = strsplit(znam)
  loop i=1..nnz -q
    txti = sprintf("z[%d]",i)
    fun2 = smartstrsub(fun2,znams[i],txti)
  endloop
else
  matrix z ={}
  string znam = ""
endif
fun2bkp = fun2
fun2 = text_xb(fun2,n1-1)
fun2s = strsplit(fun2bkp)
isscst = (n1-1) > nelem(fun2s)
fun321 = xtextfun(fun2s,1,isscst,"b")
fun321bkp = fun321
maxX = nns
maxZ = nnz
ZNAM = znam
strings nams = strsplit(namZ)
scalar nns = nelem(nams)
matrix indiz = {}
if nns > 0
  matrix indiz = zeros(nns,1)
  loop i=1..nns -q
    indiz[i] = findind(nams[i],allnamZ)
  endloop
endif
string fun2 = funZ
commonX = str_setinrersect(namZ,nam)
commonXs = strsplit(commonX)
ncX = nelem(commonXs)
aloneX = str_setdiff(namZ,nam)
aloneXs = strsplit(aloneX)
naX = nelem(aloneXs)
ALLX = maxX + naX
if ncX > 0
  coindi = findinds(commonX,nam)
  loop i=1..ncX -q
    txti = sprintf("x[%d]",coindi[i])
    fun2 = smartstrsub(fun2,commonXs[i],txti)
  endloop
endif
matrix aindiX = {}
matrix xZZ = {}
matrix aindiZ = {}
matrix zZZ = {}
if naX > 0
  matrix aindiX = findinds(aloneX,namZ)
  matrix xZZ = allvalZ[aindiX]
  loop i=1..naX -q
    txti = sprintf("x[%d]",i+maxX)
    fun2 = smartstrsub(fun2,aloneXs[i],txti)
  endloop
endif
scalar bign = nelem(strsplit(allnamZ))
scalar nnz = bign - nns
if nns < bign
  string znamZ = str_setdiff(allnamZ,namZ)
  strings znamZs = strsplit(znamZ)
  if maxZ == 0
    matrix aindiZ = findinds(znamZ,allnamZ)
    matrix zZZ = allvalZ[aindiZ]
    loop i=1..nnz -q
      txti = sprintf("z[%d]",i)
      fun2 = smartstrsub(fun2,znamZs[i],txti)
    endloop
  else
    commonZ = str_setinrersect(znamZ,ZNAM)
    commonZs = strsplit(commonZ)
    aloneZ = str_setdiff(znamZ,ZNAM)
    aloneZs = strsplit(aloneZ)
    ncZ = nelem(commonZs)
    naZ = nelem(aloneZs)
    if ncZ > 0
      coindiZ = findinds(commonZ,ZNAM)
      loop i=1..ncZ -q
        txti = sprintf("z[%d]",coindiZ[i])
        fun2 = smartstrsub(fun2,commonZs[i],txti)
      endloop
    endif
    ALLZ = maxZ + naZ
    if naZ > 0
      matrix aindiZ = findinds(aloneZ,allnamZ)
      matrix zZZ = allvalZ[aindiZ]
      loop i=1..naZ -q
        txti = sprintf("z[%d]",i+maxZ)
        fun2 = smartstrsub(fun2,aloneZs[i],txti)
      endloop
    endif
  endif
endif
fun2bkpZ = fun2
fun2s = strsplit(fun2bkpZ)
isscst = (n2-n1) > nelem(fun2s)
fun321Z = xtextfun(fun2s,n1+1,isscst,"b")
fun321Zbkp = fun321Z
x = x~xZZ
z = z~zZZ
matrix xBACKUP = x
#print link
if link != "Ey"
  fun321Z = "invmills(-("~fun321Z~"))"
  bzzzz = sprintf("b[%d]",n1)
  fun321Z = bzzzz~"*"~fun321Z
  fun421 = "log("~fun321~"+"~fun321Z~")"
  if etype == 1
    #print x
    dydx = fdjac(x,@fun421)'.*x'
    #print dydx
    funHe = fun421
    loop i=1..n2 -q
      finD = sprintf("b[%d]",i)
      repL = sprintf("x[%d]",i+ALLX)
      funHe = strsub(funHe,finD,repL)
    endloop
    x = x'|b
    #print x
    hess =  numhess(x,@funHe)[1:ALLX,(ALLX+1):].*xBACKUP'
    #print hess
    matrix VCV = qform(hess,vcv)
    matrix SES = sqrt(diag(VCV))
    matrices M = array(4)
    M[1] = 0
    M[2] = dydx
    M[3] = VCV
    M[4] = SES
    return M
  endif
else
  fun521Z = "dnorm(("~fun321Zbkp~"))"
  fun521 = "cnorm("~fun321Zbkp~")"
  fun521 = fun521~"*("~fun321~")"
  bzzzz = sprintf("b[%d]",n1)
  fun521Z = bzzzz~"*"~fun521Z
  fun721 = "log("~fun521~"+"~fun521Z~")"
  if etype == 1
    dydx = fdjac(x,@fun721)'.*x'
    #eval fdjac(x,@fun721)'
    funHe = fun721
    loop i=1..n2 -q
      finD = sprintf("b[%d]",i)
      repL = sprintf("x[%d]",i+ALLX)
      funHe = strsub(funHe,finD,repL)
    endloop
    x = x'|b
    hess =  numhess(x,@funHe)[1:ALLX,(ALLX+1):].*xBACKUP'
    matrix VCV = qform(hess,vcv)
    matrix SES = sqrt(diag(VCV))
    matrices M = array(4)
    M[1] = 0
    M[2] = dydx
    M[3] = VCV
    M[4] = SES
    return M
  else
  endif
endif
#return dydx~hess

text1 = "abs         acos        acosh       aggregate   argname        array       asin        asinh       atan        atanh        atof        bessel      BFGSmax     BFGSmin     BFGScmax        BFGScmin    bkfilt      boxcox      bread       bwfilt        bwrite      cdemean     cdf         cdiv        cdummify        ceil        cholesky    chowlin     cmult       cnorm        cnumber     colname     colnames    cols        corr        corrgm      cos         cosh        cov         critical        cum         curl        dayspan     defarray    defbundle        deflist     deseas      det         diag        diagcat        diff        digamma     dnorm       dropcoll    dsort        dummify     easterday   ecdf        eigengen    eigensym        eigsolve    epochday    errmsg      exists      exp        fcstats     fdjac       fft         ffti        filter        firstobs    fixname     floor       fracdiff    gammafun        genseries   getenv      getinfo     getline     ghk        gini        ginv        GSSmax      GSSmin      halton        hdprod      hfdiff      hfldiff     hflags      hflist        hpfilt      I           imaxc       imaxr       imhof        iminc       iminr       inbundle    infnorm     inlist        int         inv         invcdf      invmills    invpd        irf         irr         isconst     isdiscrete  isdummy        isnan       isoconv     isodate     iwishart    jsonget        juldate     kdensity    kdsmooth    kfilter     kmeier        kpsscrit    ksetup      ksimul      ksmooth     kurtosis        lags        lastobs     ldet        ldiff       lincomb        linearize   ljungbox    lngamma     loess       log        log10       log2        logistic    lower       lrvar        max         maxc        maxr        mcorr       mcov        mcovg       mean        meanc       meanr       median        mexp        mgradient   min         minc        minr        missing     misszero    mlag        mlincomb    mnormal        mols        monthlen    movavg      mpols       mrandgen        mread       mreverse    mrls        mshape      msortby        muniform    mweights    mwrite      mxtab       naalen        nadarwat    nelem       ngetenv     nlines      NMmax        NMmin       nobs        normal      normtest    npcorr        npv         NRmax       NRmin       nullspace   numhess        obs         obslabel    obsnum      ok          onenorm        ones        orthdev     pdf         pergm       pexpand        pmax        pmean       pmin        pnobs       polroots        polyfit     princomp    prodc       prodr       psd        psdroot     pshrink     psum        pvalue      pxnobs        pxsum       qform       qlrpval     qnorm       qrdecomp        quadtable   quantile    randgen     randgen1    randint        rank        ranking     rcond       readfile    regsub        remove      replace     resample    round       rownames        rows        sd          sdc         sdiff       seasonals        selifc      selifr      seq         setnote     simann        sin         sinh        skewness    sleep       smplspan        sort        sortby      sprintf     sqrt        square        sscanf      sst         stringify   strlen      strncmp        strsplit    strstr      strstrip    strsub      strvals        substr      sum         sumall      sumc        sumr        svd         tan         tanh        toepsolv    tolower        toupper     tr          transp      trimr       typeof        typestr     uniform     uniq        unvech      upper        urcpval     values      var         varname     varnames        varnum      varsimul    vec         vech        weekday        wmean       wsd         wvar        xmax        xmin        xmlget      zeromiss    zeros"
line1 = "acos acosh asin asinh atan atanh cos cosh "
line2 = "digamma exp gammafun lngamma log log10 log2 "
line3 = "logistic sin sinh sqrt tan tanh"
text2 = line1~line2~line3
noar = str_setdiff(text1,text2)
return noar

fla2 = "#"~fla~"#"
chars = "qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM"
digits = "0123456789"
lefla = strlen(fla2)
matrix charmap = zeros(lefla,1)
matrix digmap = charmap
matrix bramap = charmap
matrix undmap = charmap
loop i = 1..lefla -q
  chari = substr(fla2,i,i)
  if strlen(strstr(chars,chari))
    charmap[i] =  1
  endif
  if strlen(strstr(digits,chari))
    digmap[i] =  1
  endif
  if chari == "_"
    undmap[i] = 1
  endif
  #if (chari == "(")||(chari == "[")
  if (chari == "[")
    bramap[i] = 1
  endif
endloop
all_charmap = charmap + digmap + undmap
matrix seql = seq(0,lefla-1)'
matrix begins = seql .*(diff(all_charmap).> 0)
matrix ends = seql.*mlag((diff(all_charmap).< 0),-1)
if  sumc(uniq(begins))
  matrix numbegins = uniq(begins)[2:]
  matrix numends = uniq(ends)[2:]
  smalln = nelem(numends)
  ret = ""
  strings rets
  loop i =1..smalln -q
    summ = 0
    reti = substr(fla,numbegins[i],numends[i])
    #print reti
    nofun = substr(fla,numends[i]+1,numends[i]+1)!="("
    if !nofun
      if i==1
        rets +=reti
        ret = ret~" "~reti
      else
        nrets = nelem(rets)
        summ = 0
        loop j=1..nrets -q
          summ += reti == rets[j]
        endloop
        if !(digmap[numbegins[i]+1]+bramap[numends[i]+2]+summ)
          rets +=reti
          ret = ret~" "~reti
        endif
      endif
    endif
  endloop
  return strstrip(ret)
else
  return "_empty"
endif

nv = $nvars
list li = seq(1,nv-1)
vns = varname(li)
vns = strsub(vns,","," ")
#print vns
descr0 = strsub(descr," ","")
vars = uniqvarlist(descr)
vars0 = uniqvarlist(descr0)
#print vars vars0
ret = 1
funs = ""
if vars != vars0
  return 0
elif vars == "_empty"
  return 0
else
  funs = uniqfunlist(descr)
  #print funs
  noar = no_arith()
  matrix noaris = findinds(funs,noar)
  if nelem(noaris) > 0
    if max(noaris) > 0
      return 0
    endif
  endif
endif
#print vars funs
ck11 = strlen(strstr(descr,"="))
ck12 = strlen(strstr(descr,"<"))
ck13 = strlen(strstr(descr,">"))
if (ck11 + ck12 + ck13) > 0
  ret = 0
endif
#print ret
return ret

firstcheck = strlen(strstr(vname,"_"))
if firstcheck == 0
  return "_empty"
endif
vnn = strlen(vname)
matrix posi = {}
j = 0
loop i = 1..vnn -q
  if substr(vname,i,i) == "_"
    posi = posi~i
    j++
  endif
endloop
pref = substr(vname,1,posi[j]-1)
suff = ""
if posi[j] < vnn
  suff = substr(vname,posi[j]+1,vnn)
endif
if strlen(suff) == 0
  return "_empty"
endif
nsu = atof(suff)
isnasu = missing(nsu)
if isnasu
  return "_empty"
endif
zck = nsu == round(nsu)
if zck == 0
  return "_empty"
endif
if nsu <= 0
  return "_empty"
endif
suff2 = sprintf("%d",nsu)
if suff != suff2
  return "_empty"
else
  return pref~" "~suff
endif

include a_eff.gfn open greene22_2.gdt --quiet gretl_version = $version set echo on list z_small = Z* -Z5 Ymlog = 1+ (Y > 0) + (Y > 1) ### ### marginal effects after multinomial logit ### with group dummies and a nonlinear function ### of a continuous variable # logit Ymlog 0 z_small square(Z2) dummify(Z5) --multinomial bundle b_ml = $model mfx(&b_ml) ### ### marginal effects after tobit ### at different data points for subset of variables # tobit Y 0 z_small square(Z2) dummify(Z5) # # b[sq_Z2] has the expected sign but is insignificant, # multicollinearity? # after tobit Y 0 Z1 Z2 square(Z2) p-value for sq_Z2 is 0.0032 # bundle b_t = $model formula = smart_formula2($xlist) # # note: you can set either one of "at = Z1" or "wrt = Z8", # or the both # mfx(&b_t,formula,"wrt = Z3 Z2; at = Z1 Z2 Z5 Z8", {1,52,2,1}) # # with the same marriage Z2 and Z3 grow at the same rate # so Y grows even at Z1 = 1 and Z2 = 52! # # for E(y|lb<y<rb) # mfx(&b_t,formula,"wrt = Z2; type = inside") # ### elasticities # # for E(y) # eyex(&b_t,formula,"wrt = Z2 Z3; at = Z2 Z3 Z5",{30,10,2}) # # for E(y|lb<y<rb) # eyex(&b_t,formula,"wrt = Z2; type = inside; at = Z5",{2}) ### ### Translog production function: ### tricks with smart_formula2(), double_formula(), ### and 'type' parameter to compute ### different kinds of effects ## creating artifical data for translog production function # set seed 13 K = normal(0,3) K = K - min(K)+3+uniform() L = normal(0,2) L = 0.4*K + K/10*normal(0,1) logs L K list p_fact = l_K l_L list squar = square(p_fact,1) list xxlist = 0 p_fact squar beta = {1,0.946,1.254,-0.108,-0.063,-0.108} l_Q = lincomb(xxlist,beta)+normal(0,0.05) # Cobb -- Douglas # ols l_Q 0 p_fact --simple-print reset --squares-only -q # Translog # ols l_Q xxlist --simple-print bundle b_trl = $model # # here we have 2-stage variables creation: logs then squared logs # so default smart_formula2($xlist) does not work as desirable # eval smart_formula2($xlist) # # then we force l_K and l_L to be independent variables # efla = smart_formula2($xlist,p_fact) eval efla # # the next line wil get elasticities of output wrt. factors # at means of their logs # mfx(&b_trl,efla) # # at sample minima and maxima # bundle bm = mfx(&b_trl,efla,"at = l_L l_K",{min(l_L),min(l_K)}) # # Confidence intervals: # eval eff_confint(&bm) # note that we can test hypotheses in terms of # effects directly: # Copy-paste code region between '/*' and '/*' # into new script editor window to run # (requires waldTest.gfn v. >= 1.1 installed) /* include waldTest.gfn nlwaldtest("b[1]+b[2]=1",&bm) */ mfx(&b_trl,efla,"at = l_L l_K",{max(l_L),max(l_K)}) # double_fla() can deal with 2-stage variables creation: mfla = double_formula($xlist) eval mfla # the next line wil get marginal effects of output wrt. factors # at means mfx(&b_trl,mfla,"type = exp") # ### heckit # open mroz87.gdt --quiet kidsi = ((KL6+K618) == 0)+ 2*((KL6+K618) == 1) +3*((KL6+K618) > 1) heckit WW 0 WE AX; LFP 0 WA WE dummify(kidsi) MTR b_he = $model # # y conditional # mfx(&b_he,null,"wrt = kidsi") # # y expected # mfx(&b_he,null,"wrt = WE; type = Ey") eyex(&b_he,null,"wrt = AX; type = Ey")