set verbose off

function void my_arima(string spec,
                       list y,
                       list xlist[null],
                       string opt[null])
    list biglist = y xlist
    smpl biglist --contiguous
    cmdline = "arima "~spec~"; y"
    cmdline2 = cmdline
    nx = 0
    matrix mxs = {}
    matrix sdxs = {}
    list zlist = null

    # standardize exogenous regressors    
    if nelem(xlist) > 0
        hascst = inlist(xlist, const)
        xlist = xlist - 0
        nx = nelem(xlist)
        matrix mxs = zeros(1, nx)
        matrix sdxs = mxs
        loop i=1..nx -q
            mi = mean(xlist[i])
            sdi = sd(xlist[i])
            vni = sprintf("oleh%d",i)
            zlist += genseries(vni, (xlist[i] - mi)/sdi)
            mxs[i] = mi
            sdxs[i] = sdi
        endloop
        if hascst
            zlist = 0 zlist
            xlist = 0 xlist
        endif
        cmdline = cmdline ~ " zlist"
        cmdline2 = cmdline2 ~ " xlist"
    endif

    if !isnull(opt)
        opt = " " ~ strstrip(opt)
        cmdline = cmdline ~ opt
        cmdline2 = cmdline2 ~ opt
    endif
    
    # libgretl arima command
    catch @cmdline2
    err = $error
    if err
        em = errmsg(err)
        printf "\n"
        if strlen(em) == 0
            print "libgretl arima failed with empty errmsg()"
        else
            eval errmsg(err)
        endif
        printf "\n"
    else
        matrix b0 = $coeff
        matrix se0 = $stderr
	lnl0 = $lnl # loglik from libgretl arima
    endif
    
    # arima command for transformed data
    @cmdline
    matrix b1 = $coeff
    lnl1 = $lnl
    if nx > 0
        # restore coeff and vcv for original data
        nc = $ncoeff
        difmat = I(nc)
        b1[(nc-nx+1):] = b1[(nc-nx+1):]./sdxs'
        loop i=1..nx -q
            difmat[(nc-nx+i),(nc-nx+i)] = 1/sdxs[i]
        endloop
        if hascst
            matrix vec1 = b1[(nc-nx+1):]
            b1[1] -= mxs * vec1
            difmat[1,(nc-nx+1):] = -mxs./sdxs
        endif
        # print difmat
        # eval qform(difmat,$vcv)
        matrix se1 = sqrt(diag(qform(difmat, $vcv)))
    endif

    # print results for comparison
    printf "transformed (oleh): estimates and standard errors:\n\n%#13.5g\n",
      b1 ~ se1
    printf "    lnl = %g\n", lnl1
    if !err
        printf "\nlibgretl: estimates and standard errors:\n\n%#13.5g\n",
	  b0 ~ se0
        printf "    lnl = %g\n", lnl0
        printf "\nDifference in lnl (oleh - libgretl): %g\n", lnl1 - lnl0
    else
        printf "\nln (oleh) = %g\n", lnl1
    endif
end function

open bad_data.gdt

matrix plus = {2, 4, 1000, 10000}
series x = 0
list zlic10 = 0 y_one x

loop j=1..nelem(plus) -q
  printf "\n*** exogenous x = y_two + %d ***\n\n", plus[j]
  x = y_two + plus[j]
  my_arima("3 0 0; 1 0 0", diff_series, zlic10, "-q")
  arima 3 0 0; 1 0 0; diff_series zlic10 --x-12-arima -q
  printf "X-13-ARIMA lnl = %g\n", $lnl
endloop