function matrix armax (int ar_order[0::1] "max. AR p",
                       int ma_order[0::0] "max. MA q",
                       series y "Dependent variable",
                       list xlist[null] "Exogenous variable list",
                       bool cons[1] "Include a constant",
                       int method[1:2:1] "Method" {"exact ML", "conditional ML"},
                       bool showBestModel[0] "show Best model based on Info. Criterion below",
                       int crit[1:3:1] "Info. Criterion" {"AIC", "BIC", "HQC"},
                       bool quiet[0] "Quiet")
  # set upper limit lags
  maxp=ar_order
  maxq=ma_order
  # genseries(argname(y), y)
  # estimation method: exact ML is the default in gretl
  if method=2
    string eMethod =" --conditional"
  else
    string eMethod =""
  endif
  if cons
    string sCons = ""
    list tlist = 0
    list xlist = xlist || tlist
  else
    string sCons = " --nc"
  endif

  scalar order_aic=0
  scalar order_sc=0
  scalar order_hq=0
  # --- define decimals in output
  nwidth=12
  ndecimal=4
  # declare a Matrix to store information criteria
  matrix tab1 = zeros((maxp+1)*(maxq+1),5)
  colnames(tab1,"p q AIC BIC HQC")
  # loop through every combination of arma lags
  nrow = 1
  loop p=0..maxp -q
    ## for q=0 to maxq
    loop q = 0..maxq -q
      # estimate model
      if nelem(xlist)=0
        if p=0 && q=0
          catch arima p q ; y -q @eMethod
          err=$error
          if err
            printf "armax: Got error %d (%s) while p=%d,q=%d\n", err, errmsg(err),p,q
          else
            printf "armax: successfully-1:p=%d,q=%d\n",p,q
          endif
        else
          catch arima p q ; y -q @eMethod @sCons
          err=$error
          if err
            printf "armax: Got error %d (%s) while p=%d,q=%d\n", err, errmsg(err),p,q
          else
            # printf "armax: successfully-1:p=%d,q=%d\n",p,q
          endif
        endif
      else
        catch arima p q ; y xlist -q @eMethod @sCons
        err=$error
        if err
          printf "armax: Got error %d (%s) while p=%d,q=%d\n", err, errmsg(err),p,q
        endif
      endif

      tab1[nrow,1] = p
      tab1[nrow,2] = q
      if err>0
        tab1[nrow,3] = 99999.9999
        tab1[nrow,4] = 99999.9999
        tab1[nrow,5] = 99999.9999
      else
        tab1[nrow,3] = $aic
        tab1[nrow,4] = $bic
        tab1[nrow,5] = $hqc
        # test for mininum noise model
        if nrow=1
          min_aic = $aic
          min_sc = $bic
          min_hq = $hqc
          order_aic = nrow
          order_sc = nrow
          order_hq = nrow
        else
          if $aic < min_aic
            min_aic = $aic
            order_aic = nrow
          endif
          if $bic < min_sc
            min_sc = $bic
            order_sc = nrow
          endif
          if $hqc < min_hq
            min_hq = $hqc
            order_hq = nrow
          endif
        endif
      endif
      nrow = nrow+1
    endloop
  endloop
  # --- remember the best model
  scalar naic=order_aic
  scalar nbic=order_sc
  scalar nhqc=order_hq
  if !quiet
    printf "===============================================\n"
    printf " Information Criteria of ARMAX(p,q) for %s\n",argname(y)
    printf "-----------------------------------------------\n"
    printf " p, q %*s %*s %*s\n",nwidth,"AIC",nwidth+1,"BIC",nwidth+1,"HQC"
    printf "-----------------------------------------------\n"
    nrow=1
    loop p=0..maxp -q
      loop q=0..maxq -q
        if nrow==naic
          string star_aic="*"
        else
          string star_aic=""
        endif
        if nrow=nbic
          string star_bic="*"
        else
          string star_bic=""
        endif
        if nrow=nhqc
          string star_hqc="*"
        else
          string star_hqc=""
        endif
        printf "%2d,%2d %*.*f%1s %*.*f%1s %*.*f%1s\n",p,q,nwidth,ndecimal,tab1[nrow,3],star_aic, nwidth,ndecimal,tab1[nrow,4],star_bic, nwidth,ndecimal,tab1[nrow,5],star_hqc
        nrow = nrow + 1
      endloop
    endloop
    printf "===============================================\n"
    printf "* indicates best models.\n"
    printf "'9999.9999' suggests failures to estimate the models.\n"
    #    printf "%10.4f\n",tab1
    if showBestModel
      scalar aic_p=tab1[naic,1]
      scalar aic_q=tab1[naic,2]
      scalar bic_p=tab1[nbic,1]
      scalar bic_q=tab1[nbic,2]
      scalar hqc_p=tab1[nhqc,1]
      scalar hqc_q=tab1[nhqc,2]
      if crit=1
        best_p = aic_p
        best_q = aic_q
      elif crit=2
        best_p = bic_p
        best_q = bic_q
      elif crit=3
        best_p = hqc_p
        best_q = hqc_q
      endif
      # show the best model output
      if best_p==0 && best_q==0
        arima best_p best_q ; y xlist @eMethod
      else
        arima best_p best_q ; y xlist @eMethod @sCons
      endif
    endif
  endif
  return tab1
end function