"Fossies" - the Fresh Open Source Software Archive

Member "gretl/functions/gig/gig.gfn" (21 Nov 2020, 48383 Bytes) of package /windows/misc/gretl-2020e-win32.zip:


As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 <?xml version="1.0" encoding="UTF-8"?>
    2 <gretl-functions>
    3 <gretl-function-package name="gig" needs-time-series-data="true" minver="2016b" lives-in-subdir="true">
    4 <author email="r.lucchetti@univpm.it">Riccardo &quot;Jack&quot; Lucchetti and Stefano Balietti</author>
    5 <version>2.24</version>
    6 <date>2020-11-10</date>
    7 <description>An assortment of univariate GARCH models</description>
    8 <tags>C22</tags>
    9 <label>gig</label>
   10 <help>
   11 pdfdoc:gig.pdf
   12 </help>
   13 <data-files count="1">
   14 examples </data-files>
   15 <gretl-function name="GUI_gig" type="bundle" pkg-role="gui-main">
   16  <params count="11">
   17   <param name="y" type="series">
   18 <description>Dependent Variable</description>
   19   </param>
   20   <param name="type" type="int" min="1" max="7" default="1">
   21 <description>Model type</description>
   22 <labels count="7">
   23 "GARCH" "Taylor/Schwert GARCH" "GJR" "TARCH" "NARCH" "APARCH" "EGARCH" </labels>
   24   </param>
   25   <param name="p" type="int" default="1">
   26 <description>GARCH</description>
   27   </param>
   28   <param name="q" type="int" default="1">
   29 <description>ARCH</description>
   30   </param>
   31   <param name="X" type="list" optional="true">
   32 <description>Mean regressors</description>
   33   </param>
   34   <param name="hasconst" type="bool" default="1">
   35 <description>Constant</description>
   36   </param>
   37   <param name="ARlags" type="int" min="0" default="0">
   38 <description>AR lags</description>
   39   </param>
   40   <param name="Y" type="list" optional="true">
   41 <description>Variance regressors</description>
   42   </param>
   43   <param name="cdist" type="int" min="0" max="4" default="0">
   44 <description>Distribution</description>
   45 <labels count="5">
   46 "Normal" "t" "GED" "Skewed t" "Skewed GED" </labels>
   47   </param>
   48   <param name="vtype" type="int" min="0" max="2" default="0">
   49 <description>Covariance estimator</description>
   50 <labels count="3">
   51 "Sandwich" "Hessian" "OPG" </labels>
   52   </param>
   53   <param name="verb" type="int" min="0" max="2" default="1">
   54 <description>Verbosity</description>
   55   </param>
   56  </params>
   57 <code>if hasconst
   58   list X = const || X
   59 endif
   60 bundle m = gig_setup(y, type, X, Y, ARlags)
   61 # this has to be forced
   62 m.depvarname = argname(y)
   63 gig_set_dist(&amp;m, cdist)
   64 gig_set_pq(&amp;m, p, q)
   65 m.vcvtype = vtype
   66 gig_estimate(&amp;m, verb)
   67 return m
   68 </code>
   69 </gretl-function>
   70 <gretl-function name="gig_setup" type="bundle">
   71  <params count="5">
   72   <param name="y" type="series"/>
   73   <param name="type" type="int" default="1"/>
   74   <param name="X" type="list" optional="true"/>
   75   <param name="Y" type="list" optional="true"/>
   76   <param name="ARlags" type="scalar" default="0"/>
   77  </params>
   78 <code>/*
   79 autoscaling is much, much trickier that one would think: assume we
   80 have a suitable number to scale the dep. var.; what should we
   81 do with regressors? Surely we need to scale the lags accordingly, but we don't have atm a clean way to tell them apart.
   82 At the moment, we leave the regressors alone and leave all the
   83 scale accounting (coeff, vcv modification etc) for post-estimation.
   84 */
   85 envvar = ngetenv(&quot;NOSCALE&quot;)
   86 AUTOSCALE = missing(envvar) || (envvar == 0)
   87 bundle model
   88 model.type = type
   89 model.AR = ARlags
   90 if (ARlags&gt;0)
   91   list fullX = X || lags(ARlags, y)
   92 else
   93   list fullX = X
   94 endif
   95 list everything = y || fullX || Y
   96 smpl everything --no-missing
   97 model.t1 = $t1
   98 model.t2 = $t2
   99 if (AUTOSCALE == 1) &amp;&amp; (type&lt;7) # skip EGARCH for the moment
  100   if nelem(X) &gt; 0
  101     ols y X --quiet
  102     scalar scale = $sigma
  103   else
  104     scale = sd(y)
  105   endif
  106   model.scale = scale
  107   model.y = y/scale
  108 else
  109   model.scale = 1
  110   model.y = y
  111 endif
  112 model.depvarname = argname(y)
  113 model.q = 1
  114 model.p = (type&gt;0)
  115 model.nobs = $nobs
  116 k = nelem(X)
  117 model.mk = k + ARlags
  118 if model.mk == 0
  119   matrix mlistX = {}
  120   matrix mX = {}
  121   model.mXnames = &quot;&quot;
  122 else
  123   if nelem(X)&gt;0
  124     matrix mlistX = X
  125   else
  126     matrix mlistX = {}
  127   endif
  128   matrix mX = { fullX }
  129   string names = varname(X)
  130   if (ARlags&gt;0)
  131     loop i = 1 .. ARlags
  132       names = sprintf(&quot;%s,AR%d&quot;, names, i)
  133     endloop
  134   endif
  135   model.mXnames = names
  136 endif
  137 model.mlistX = mlistX
  138 model.mX = mX
  139 list lY = const || Y
  140 vk = nelem(lY)
  141 model.vk = vk
  142 matrix vlistX = lY
  143 matrix mX = { lY }
  144 model.vXnames = varname(lY)
  145 model.vlistX = vlistX
  146 model.vX = mX
  147 gig_set_vQR(&amp;model, 0) # disable for now
  148 model.cdist = 0
  149 /* initialisation */
  150 coef_init(&amp;model)
  151 set_active_par(&amp;model)
  152 /* robust */
  153 model.vcvtype = 0
  154 return model
  155 </code>
  156 </gretl-function>
  157 <gretl-function name="gig_set_dist" type="void">
  158  <params count="2">
  159   <param name="model" type="bundleref"/>
  160   <param name="cdist" type="scalar"/>
  161  </params>
  162 <code>old_cdist = model.cdist
  163 if (old_cdist != cdist)
  164   model.cdist = cdist
  165   n_old = n_cdist_par(old_cdist)
  166   n_new = n_cdist_par(cdist)
  167   active = model.active
  168   params = model.coeff
  169   noldpar = rows(params)
  170   # zap existing parameters
  171   if n_old &gt; 0
  172     active = trimr(active,0,n_old)
  173     params = trimr(params,0,n_old)
  174   endif
  175   if n_new &gt; 0
  176     active |= (rows(params)+seq(1,n_new)')
  177     params |= cdist_initpar(cdist)
  178   endif
  179   model.active = active
  180   model.coeff = params
  181 endif
  182 </code>
  183 </gretl-function>
  184 <gretl-function name="gig_set_pq" type="void">
  185  <params count="3">
  186   <param name="model" type="bundleref"/>
  187   <param name="p" type="scalar"/>
  188   <param name="q" type="scalar"/>
  189  </params>
  190 <code>scalar old_p = model.p
  191 scalar old_q = model.q
  192 if (old_p != p) || (old_q != q)
  193   model.p = p
  194   model.q = q
  195   theta = model.coeff
  196   mk = model.mk
  197   if (mk&gt;0)
  198     matrix mpar = theta[1:mk]
  199   else
  200     matrix mpar = {}
  201   endif
  202   dpar = distpar(model.cdist, theta)
  203   vpar = var_coef_init(&amp;model, model.s2)
  204   model.coeff = mpar | vpar | dpar
  205   set_active_par(&amp;model)
  206 endif
  207 </code>
  208 </gretl-function>
  209 <gretl-function name="gig_set_vQR" type="void">
  210  <params count="2">
  211   <param name="model" type="bundleref"/>
  212   <param name="on_off" type="bool"/>
  213  </params>
  214 <code>model.vX_QR = on_off
  215 if model.vX_QR
  216   matrix vX_R
  217   matrix QvX = vreg_QR(model.vX, &amp;vX_R)
  218   model.vX_R = vX_R
  219   model.QvX = QvX
  220 endif
  221 </code>
  222 </gretl-function>
  223 <gretl-function name="gig_set_vcvtype" type="void">
  224  <params count="2">
  225   <param name="model" type="bundleref"/>
  226   <param name="s" type="string"/>
  227  </params>
  228 <code>s = tolower(s)
  229 if s == &quot;sandwich&quot;
  230   model.vcvtype = 0
  231 elif s == &quot;hessian&quot;
  232   model.vcvtype = 1
  233 elif s == &quot;opg&quot;
  234   model.vcvtype = 2
  235 endif
  236 </code>
  237 </gretl-function>
  238 <gretl-function name="gig_print" type="void">
  239  <params count="2">
  240   <param name="b" type="bundleref"/>
  241   <param name="verbose" type="scalar" min="0" max="3" default="0"/>
  242  </params>
  243 <code>err = b.errcode
  244 if err == 0
  245   gig_printHeader(&amp;b)
  246   maybe_gig_print(&amp;b, verbose)
  247 else
  248   printf &quot;\nESTIMATION FAILED: %s\n&quot;, errmsg(err)
  249 endif
  250 </code>
  251 </gretl-function>
  252 <gretl-function name="gig_estimate" type="scalar">
  253  <params count="2">
  254   <param name="model" type="bundleref"/>
  255   <param name="verbose" type="int" min="0" max="2" default="1"/>
  256  </params>
  257 <code># verbose = 0 -&gt; quiet
  258 # verbose = 1 -&gt; show output
  259 # verbose = 2 -&gt; show iterations
  260 scalar garchType = model.type
  261 t1 = model.t1
  262 t2 = model.t2
  263 smpl t1 t2
  264 if garchType &lt;= 7 #APARCH &amp; EGARCH
  265   scalar err = do_mle(&amp;model, verbose==2)
  266 elif garchType == 8
  267   # GARCH-in-Mean
  268   # Not yet, Buster!
  269   # if verbose &gt; 1
  270   #     printParList( parIndex, parList, rows(distrInit))
  271   #     printf &quot;Proceeding with GARCH-in-mean by Engle et al. estimation.\n\n&quot;
  272   # endif
  273   # ret = garchmean( depVar, meanX, lagOrders, garchInit, distrType,\
  274   # 		distrInit, verbose, vcv)
  275   printf &quot;Not yet, Buster!\n&quot;
  276 endif
  277 if model.scale != 1
  278   s = model.scale
  279   adjust_for_scale(&amp;model, 1/s)
  280 endif
  281 if (verbose&gt;0)
  282   gig_print(&amp;model)
  283 endif
  284 return err
  285 </code>
  286 </gretl-function>
  287 <gretl-function name="gig_var_fcast" type="matrix">
  288  <params count="3">
  289   <param name="mod" type="bundleref"/>
  290   <param name="horizon" type="scalar"/>
  291   <param name="rep" type="scalar"/>
  292  </params>
  293 <code>DBG = 0
  294 matrix ret = {}
  295 # checks
  296 if mod.type == 7 || mod.vk &gt; 1 # EGARCH or exo in variance
  297   printf &quot;Variance forecasting not available for this model (yet)\n&quot;
  298   return ret
  299 endif
  300 if mod.type == 2 ||  mod.type &gt; 3
  301   printf &quot;FIXME: results may be fishy for delta != 2\n&quot;
  302 endif
  303 # fetch stuff from the bundle first
  304 loop foreach i p q mk
  305   scalar $i = mod.$i
  306 endloop
  307 # matrix e = s2m(mod.stduhat)
  308 matrix e = s2m(mod.uhat)
  309 matrix h = s2m(mod.h)
  310 scalar ini = mk+1
  311 scalar fin = mk+1
  312 scalar omega = mod.coeff[ini:fin]
  313 ini = fin+1
  314 fin += q
  315 matrix alpha = mod.coeff[ini:fin]
  316 ini = fin+1
  317 fin += q
  318 matrix gamma = mod.coeff[ini:fin]
  319 ini = fin+1
  320 fin += p
  321 matrix beta = mod.coeff[ini:fin]
  322 scalar delta = mod.coeff[fin+1]
  323 # --- debug ----------------------------------
  324 if DBG
  325   printf &quot;real_do_fcast:\n&quot;
  326   print omega alpha gamma beta delta
  327 endif
  328 # --------------------------------------------
  329 if horizon &gt; 1
  330   matrix sel = ceil(muniform(horizon-1,rep) * rows(e)) # b'strap selector
  331 endif
  332 scalar back = xmax(p,q)
  333 matrix h0 = h[rows(h)-back+1:]
  334 matrix e0 = e[rows(e)-back+1:]
  335 # --- debug ----------------------------------
  336 if DBG
  337   printf &quot;real_do_fcast:\n&quot;
  338   print e0 h0
  339 endif
  340 # --------------------------------------------
  341 matrix ret = zeros(rep, horizon)
  342 loop i = 1 .. rep
  343   if horizon &gt; 1
  344     matrix ei = e[sel[,i]]
  345     scalar chk = sumc(ei)
  346     if !ok(chk)
  347       printf &quot;%16.8f\n&quot;, ei ~ sel[,i]
  348     endif
  349   else
  350     ei = {}
  351   endif
  352   ret[i,] = aparch_fcast(ei, e0, h0, omega, alpha, gamma, beta, delta)'
  353 endloop
  354 return ret
  355 </code>
  356 </gretl-function>
  357 <gretl-function name="gig_plot" type="void">
  358  <params count="1">
  359   <param name="model" type="bundleref"/>
  360  </params>
  361 <code>string tmpfile = gig_grph(model.depvarname, model.uhat, model.h)
  362 gnuplot --input=&quot;@tmpfile&quot; --output=&quot;display&quot;
  363 </code>
  364 </gretl-function>
  365 <gretl-function name="gig_dplot" type="void">
  366  <params count="1">
  367   <param name="model" type="bundleref"/>
  368  </params>
  369 <code>series eps = model.stduhat
  370 matrix X = kdensity(eps)
  371 scalar dtype = model.cdist
  372 theta = model.coeff
  373 ncoef = rows(theta)
  374 e = X[,1]
  375 if dtype == 0 # Normal
  376   matrix d = pdf(z, e)
  377   string descstr = &quot;Std Normal&quot;
  378 elif dtype == 1 # Student's t
  379   scalar ni = theta[ncoef]
  380   scalar hadj = sqrt(1-2/ni)
  381   matrix d = pdf(t, ni, e/hadj)/hadj
  382   descstr = sprintf(&quot;t(%g)&quot;, ni)
  383 elif dtype == 2 # GED
  384   scalar ni = theta[ncoef]
  385   scalar p = 1/ni
  386   scalar lg1 = lngamma(p)
  387   scalar lg3 = lngamma(3*p)
  388   scalar lC = ln(ni/2) + 0.5*(lg3 - 3*lg1)
  389   scalar k = exp(0.5*(lg1-lg3)) * (0.5^p)
  390   matrix u = abs(e)/k
  391   matrix d = exp(lC - 0.5*(u.^ni))
  392   descstr = sprintf(&quot;GED(%g)&quot;, ni)
  393 elif dtype == 3 # Skewed-T
  394   scalar ni = theta[ncoef-1]
  395   scalar alpha = tanh(theta[ncoef])
  396   matrix d = pdf_skt(&amp;e, ni, alpha)
  397   descstr = sprintf(&quot;Skewed t(%g, %g)&quot;, ni, alpha)
  398 elif dtype == 4 # Skewed-GED
  399   scalar ni = theta[ncoef-1]
  400   scalar alpha = tanh(theta[ncoef])
  401   matrix d = pdf_skged(&amp;e, ni, alpha)
  402   descstr = sprintf(&quot;Skewed GED(%g, %g)&quot;, ni, alpha)
  403 endif
  404 string tmpfile = gig_dgrph(X~d, descstr)
  405 gnuplot --input=&quot;@tmpfile&quot; --output=&quot;display&quot;
  406 </code>
  407 </gretl-function>
  408 <gretl-function name="gig_vfgraph" type="void">
  409  <params count="4">
  410   <param name="f" type="matrix" const="true"/>
  411   <param name="mod" type="bundle"/>
  412   <param name="before" type="scalar"/>
  413   <param name="alpha" type="scalar"/>
  414  </params>
  415 <code>scalar horizon = cols(f)
  416 scalar a = (1-alpha)/2
  417 q0 = quantile(f, a)'
  418 q1 = quantile(f, 1-a)'
  419 m = meanc(f)'
  420 me = quantile(f,0.5)'
  421 matrix h = s2m(mod.h)
  422 matrix hinit = h[rows(h)-before:]
  423 matrix gp1 = ( hinit | me ) ~ seq(0,horizon+before)'
  424 matrix gp2 = ( hinit ~ 0 ) | ((q0+q1) ~ (q1-q0))/2
  425 string ylab = sprintf(&quot;set ylabel 'conditional variance for %s'&quot;,    mod.depvarname)
  426 string titlelab = sprintf(&quot;set title \&quot;%s\&quot;&quot;, s_modeltype(&amp;mod))
  427 plot gp1
  428   option --with-lines
  429   option --fit=none
  430   option --band=gp2
  431   option --band-style=fill
  432   literal @ylab
  433   literal @titlelab
  434   literal unset xlabel
  435 end plot --output=display
  436 </code>
  437 </gretl-function>
  438 <gretl-function name="gig_bundle_print" type="void" pkg-role="bundle-print">
  439  <params count="1">
  440   <param name="b" type="bundleref"/>
  441  </params>
  442 <code>gig_print(&amp;b, 1)
  443 </code>
  444 </gretl-function>
  445 <gretl-function name="GUI_gig_plot" type="void" pkg-role="bundle-plot">
  446  <params count="2">
  447   <param name="model" type="bundleref"/>
  448   <param name="ptype" type="int" min="0" max="2" default="0">
  449 <description>Plot type</description>
  450 <labels count="3">
  451 "Time series" "Density" "Forecast" </labels>
  452   </param>
  453  </params>
  454 <code>if ptype == 0
  455   gig_plot(&amp;model)
  456 elif ptype == 1
  457   gig_dplot(&amp;model)
  458 elif ptype == 2
  459   matrix tmp = gig_var_fcast(&amp;model, 15, 1024)
  460   gig_vfgraph(tmp, model, 60, 0.9)
  461 endif
  462 </code>
  463 </gretl-function>
  464 <gretl-function name="ln_pdf_skt" type="series" private="1">
  465  <params count="3">
  466   <param name="u" type="seriesref"/>
  467   <param name="df" type="scalar"/>
  468   <param name="ht_skew" type="scalar"/>
  469  </params>
  470 <code>series ret = NA
  471 if (df&gt;2)
  472   # sqrt(pi) = 1.77245385090551602729816748334
  473   scalar q = lngamma((df+1)/2) - lngamma(df/2)
  474   scalar c = exp(q)/(sqrt(df-2)*1.77245385090551602729816748334)
  475   scalar a = 4 * ht_skew * c * ((df-2)/(df-1))
  476   scalar b = sqrt(1 + 3*ht_skew^2 - a^2)
  477   series d = (b*u + a)
  478   d = (d&lt;0) ? d/(1-ht_skew) : d/(1+ht_skew)
  479   ret = log(b) + log(c) - ((df+1)/2) * log(1+(d^2/(df-2)))
  480 endif
  481 return ret
  482 </code>
  483 </gretl-function>
  484 <gretl-function name="ln_pdf_skged" type="series" private="1">
  485  <params count="3">
  486   <param name="x" type="seriesref"/>
  487   <param name="ni" type="scalar"/>
  488   <param name="ta" type="scalar"/>
  489  </params>
  490 <code>scalar p  = 1/ni
  491 lgp  = lngamma(p)
  492 lg2p = lngamma(2*p)
  493 lg3p = lngamma(3*p)
  494 tap1 = 1 + ta
  495 tam1 = 1 - ta
  496 scalar beta = 0.5 * exp(lg3p - lgp) * (tap1^3 + tam1^3) - 4*ta^2 * exp(2 * (lg2p - lgp))
  497 beta = sqrt(beta)
  498 # m: mode
  499 scalar m = - 2*ta/beta * exp( lg2p - lgp )
  500 scalar lnorm = log(0.5 * beta) - lngamma(p+1)
  501 series z = (x&lt;m) ? (m-x)*beta/tam1 : (x-m)*beta/tap1
  502 ret = lnorm - (z^ni)
  503 return ret
  504 </code>
  505 </gretl-function>
  506 <gretl-function name="gig_loglik" type="series" private="1">
  507  <params count="7">
  508   <param name="e" type="seriesref"/>
  509   <param name="h" type="seriesref"/>
  510   <param name="distrType" type="scalar"/>
  511   <param name="addpar" type="matrix"/>
  512   <param name="de" type="seriesref"/>
  513   <param name="dh" type="seriesref"/>
  514   <param name="dd" type="matrixref"/>
  515  </params>
  516 <code>series ret = NA
  517 de = NA
  518 dh = NA
  519 dd = {}
  520 if distrType == 0 # Normal
  521   series u = e/h
  522   ret = -.91893853320467274177 - 0.5*(log(h) + e*u)
  523   de = -u
  524   dh = -0.5/h * (1 - e*u)
  525 elif distrType == 1 # Student's t
  526   scalar ni = addpar[1]
  527   if (ni&gt;2)
  528     # series hadj = sqrt(h*(1-2/ni))
  529     # series u = e/hadj
  530     # ret = log(pdf(t, ni, u)/hadj)
  531     series e2 = e*e
  532     scalar K1 = lngamma((ni+1)/2) - lngamma(ni/2) - 0.5*log($pi*(ni-2))
  533     series ret = K1 - 0.5 * log(h) - 0.5*(ni+1) * log(1 + e2/(h*(ni-2)))
  534     series den = e2 + (ni-2)*h
  535     de = - (ni + 1) * e / den
  536     dh = 0.5/h  * ((ni + 1)* e2 / den - 1)
  537     scalar k1 = digamma((ni+1)/2) - digamma(ni/2) - 1/(ni-2)
  538     series s1 = (ni + 1)/(ni - 2) * e2/den
  539     series s2 = log(1 + e2/(h*(ni-2)))
  540     matrix dd = 0.5*(k1 + s1 - s2)
  541   endif
  542 elif distrType == 2 # GED
  543   scalar ni = addpar[1]
  544   if (ni&gt;0)
  545     scalar p = 1/ni
  546     scalar lg1 = lngamma(p)
  547     scalar lg3 = lngamma(3*p)
  548     scalar lC = log(ni/2) + 0.5*(lg3 - 3*lg1)
  549     scalar k = exp(0.5*(lg1-lg3)) * (0.5^p)
  550     series u = abs(e)/(k*sqrt(h))
  551     ret = lC - 0.5*(u^ni + log(h))
  552   endif
  553 elif distrType == 3 # Skewed-T
  554   scalar ni = addpar[1]
  555   if (ni&gt;2)
  556     series u = e/sqrt(h)
  557     alpha    = tanh(addpar[2])
  558     ret      = ln_pdf_skt(&amp;u, ni, alpha) - 0.5*log(h)
  559   endif
  560 elif distrType == 4 # Skewed-GED
  561   scalar ni = addpar[1]
  562   if (ni&gt;0)
  563     series u = e/sqrt(h)
  564     alpha    = tanh(addpar[2])
  565     ret      = ln_pdf_skged(&amp;u, ni, alpha) - 0.5*log(h)
  566   endif
  567 endif
  568 return ret
  569 </code>
  570 </gretl-function>
  571 <gretl-function name="aparchFilter" type="scalar" private="1">
  572  <params count="11">
  573   <param name="depVar" type="seriesref"/>
  574   <param name="h" type="seriesref"/>
  575   <param name="e" type="seriesref"/>
  576   <param name="mReg" type="matrixref"/>
  577   <param name="vX" type="matrixref"/>
  578   <param name="p" type="scalar"/>
  579   <param name="q" type="scalar"/>
  580   <param name="parameters" type="matrixref"/>
  581   <param name="is_asymmetric" type="scalar"/>
  582   <param name="deriv_h" type="matrixref" optional="true"/>
  583   <param name="deriv_e" type="matrixref" optional="true"/>
  584  </params>
  585 <code>nmX = cols(mReg)
  586 nvX = cols(vX)
  587 scalar base = nmX + nvX
  588 a0pos = base + 1
  589 a1pos = a0pos + q -1
  590 g0pos = a1pos + 1
  591 g1pos = g0pos + q - 1
  592 b0pos = g1pos + 1
  593 b1pos = b0pos + p - 1
  594 dpos  = b1pos + 1
  595 matrix avec = {0}
  596 matrix gvec = {0}
  597 matrix bvec = {0}
  598 matrix omegas = parameters[nmX+1:base]
  599 if q &gt; 0
  600   matrix avec = parameters[a0pos:a1pos]
  601   matrix gvec = parameters[g0pos:g1pos]
  602 endif
  603 if p &gt; 0
  604   matrix bvec = parameters[b0pos:b1pos]
  605 endif
  606 delta = parameters[dpos]
  607 err = 0
  608 # Checking
  609 # checks on alpha &amp; beta are disabled
  610 err = err || ((nvX==1) &amp;&amp; omegas[1] &lt; 0)
  611 #  err = err || (sumc(avec) + sumc(bvec)) &gt; 1
  612 #  err = err || minc(avec | bvec) &lt; 0
  613 err = err || delta &lt;= 0
  614 # shape? gamma?
  615 minh = 0
  616 if err == 0
  617   # handle the conditional mean first
  618   if nmX &gt; 0
  619     series e = depVar - (mReg * parameters[1:nmX])
  620   else
  621     series e = depVar
  622   endif
  623   scalar e0 = mean(e^2)
  624   matrix tmp_ae = mlag({abs(e)}, seq(1,q), e0 ^ (1/delta))
  625   if is_asymmetric == 1
  626     elag = mlag({e}, seq(1,q))
  627     tmp_ae -= elag .* gvec'
  628   endif
  629   # Raising a negative value to a non-integer
  630   # produces a complex number
  631   err = (delta != floor(delta)) &amp;&amp; (minc(minr(tmp_ae)) &lt; 0)
  632   if (!err)
  633     Kd = tmp_ae .^ delta
  634     series h = Kd * avec
  635     # Var regressors
  636     if (nvX&gt;1)
  637       series h += vX * omegas
  638     else
  639       series h += omegas[1]
  640     endif
  641     if p&gt;0
  642       h = filter(h, null, bvec, e0)
  643     endif
  644     err = min(h)&lt;0 || max(!ok(h))
  645     # the loglikelihood function needs sigma^2
  646     if !err &amp;&amp; (delta != 2)
  647       h = h^(2/delta)
  648     endif
  649   endif
  650 endif
  651 if (!err &amp;&amp; exists(deriv_h) &amp;&amp; exists(deriv_e))
  652   # FIXME: incomplete and experimental --------------------------------
  653   #
  654   # deriv_e and deriv_h should contain (eventually), the derivatives
  655   # of (doh!) e and h, WITH RESPECT TO THE PARAMETERS
  656   # what we have atm is a rough attempt to have it working in the GARCH
  657   # case; we'll see about generalising it later
  658   # ----------- mean eq. --------------
  659   if nmX &gt; 0
  660     deriv_e = -mReg
  661   else
  662     deriv_e = {}
  663   endif
  664   scalar zcols = nvX + p + q * (1+is_asymmetric)
  665   deriv_e ~= zeros(rows(vX), zcols)
  666   # ----------- var eq. ---------------
  667   # omega comes 2nd from last
  668   # mu comes last
  669   # alphas
  670   me = Kd
  671   matrix eeff = tmp_ae.^(delta-1)
  672   if is_asymmetric == 1 # deltas
  673     mh = -delta .* eeff .* ( elag .* avec' )
  674     me ~= mh
  675   endif
  676   if (p &gt; 0) # betas
  677     mh = mlag({h}, seq(1, p), e0)
  678     me ~= mh
  679   endif
  680   dr = rows(me) - rows(vX)
  681   if (dr&gt;0)
  682     deriv_h = ( zeros(dr, cols(vX)) | vX ) ~ me
  683   else
  684     deriv_h = vX ~ me
  685   endif
  686   if nmX &gt; 0
  687     series tmpser = (e&gt;0) ? -1 : 1
  688     matrix sgn = {tmpser}
  689     matrix mfocs = meanc({e} .* mReg)
  690     matrix dmu = zeros(rows(deriv_h), nmX)
  691     loop for i = 1 .. q
  692       matrix focs = eeff[,i] .* mlag(mReg, i)
  693       matrix tmpmat = (mlag(sgn, i) + gvec[i]) .* focs
  694       dmu += avec[i] .* tmpmat
  695     endloop
  696     dmu = dmu .* delta
  697     dmu[1,] = -2*mfocs*sumc(avec|bvec)
  698     deriv_h = dmu ~ deriv_h
  699   endif
  700   if p &gt; 0
  701     loop for i = 1 .. cols(deriv_h)
  702       series tmpser = deriv_h[,i]
  703       tmpser = filter(tmpser, null, bvec)
  704       deriv_h[,i] = tmpser
  705     endloop
  706   endif
  707   if (delta != 2)
  708     deriv_h =  h^(delta/2) .* deriv_h .* (2/delta)
  709   endif
  710 endif
  711 # -------------------------------------------------------------------
  712 return err
  713 </code>
  714 </gretl-function>
  715 <gretl-function name="egarchFilter" type="scalar" private="1">
  716  <params count="8">
  717   <param name="y" type="seriesref"/>
  718   <param name="h" type="seriesref"/>
  719   <param name="u" type="seriesref"/>
  720   <param name="mReg" type="matrixref"/>
  721   <param name="vReg" type="matrixref"/>
  722   <param name="p" type="scalar"/>
  723   <param name="q" type="scalar"/>
  724   <param name="params" type="matrixref"/>
  725  </params>
  726 <code>scalar n_mX = cols(mReg)
  727 scalar n_vX = cols(vReg)
  728 scalar err = 0
  729 if n_mX == 0
  730   series e = y
  731 else
  732   series e = y - (mReg * params[1:n_mX])
  733 endif
  734 series u = misszero(e) # Reassigns residuals for the next computation
  735 scalar omegaini = n_mX + 1
  736 scalar base = n_mX + n_vX
  737 if n_vX == 1
  738   series omega = params[omegaini]
  739 else
  740   series omega = vReg * params[omegaini:base]
  741 endif
  742 series logh = log(var(e))
  743 # ERRORS
  744 series d = (e&gt;0)
  745 series ae = abs(e)
  746 if (p == 0) &amp;&amp; (q&lt;3)
  747   if q == 1
  748     scalar g = params[base+1]
  749     scalar a = params[base+2]
  750     series tmp = d(-1) ? (g+a) : (g-a)
  751     series logh = omega + tmp*ae(-1)/exp(logh(-1)*0.5)
  752   elif q == 2
  753     scalar g1 = params[base+1]
  754     scalar g2 = params[base+2]
  755     scalar a1 = params[base+3]
  756     scalar a2 = params[base+4]
  757     series tmp1 = d(-1) ? (g1+a1) : (g1-a1)
  758     series tmp2 = d(-2) ? (g2+a2) : (g2-a2)
  759     series logh = omega + tmp1*ae(-1)/exp(logh(-1)*0.5)   + tmp2*ae(-2)/exp(logh(-2)*0.5)
  760   endif
  761 elif (p == 1) &amp;&amp; (q&lt;3)
  762   if q == 1
  763     scalar g = params[base+1]
  764     scalar a = params[base+2]
  765     scalar b = params[base+3]
  766     if (b&lt;1)
  767       series tmp = d(-1) ? (g+a) : (g-a)
  768       series logh = omega + tmp*ae(-1)/exp(logh(-1)*0.5) + b*logh(-1)
  769     else
  770       err = 1
  771     endif
  772   elif q == 2
  773     scalar g1 = params[base+1]
  774     scalar g2 = params[base+2]
  775     scalar a1 = params[base+3]
  776     scalar a2 = params[base+4]
  777     scalar b  = params[base+5]
  778     series tmp1 = d(-1) ? (g1+a1) : (g1-a1)
  779     series tmp2 = d(-2) ? (g2+a2) : (g2-a2)
  780     series logh = omega + tmp1*ae(-1)/exp(logh(-1)*0.5)) + tmp2*ae(-2)/exp(logh(-2)*0.5))   + b*logh(-1)
  781   endif
  782 else
  783   string evalstr = &quot;omega&quot;
  784   loop for i = 1 .. q
  785     scalar a$i = params[base+i]
  786     scalar g$i = params[base+i+p+q]
  787     evalstr += &quot; + a$i*abs(e(-$i)/exp(logh(-$i)*0.5)) + g$i*e(-$i)/exp(logh(-$i)*0.5)&quot;
  788   endloop
  789   loop for i = 1 .. p
  790     scalar b$i = params[base+p+i]
  791     evalstr += &quot; + b$i*logh(-$i)&quot;
  792   endloop
  793   series logh = @evalstr
  794 endif
  795 if (!err)
  796   if max(logh) &gt; 100
  797     #Overflow check
  798     logh = NA
  799   else
  800     h = exp(logh)
  801   endif
  802 endif
  803 return err
  804 </code>
  805 </gretl-function>
  806 <gretl-function name="gfilter" type="scalar" private="1">
  807  <params count="11">
  808   <param name="type" type="scalar"/>
  809   <param name="depVar" type="seriesref"/>
  810   <param name="h" type="seriesref"/>
  811   <param name="e" type="seriesref"/>
  812   <param name="mReg" type="matrixref"/>
  813   <param name="vReg" type="matrixref"/>
  814   <param name="p" type="scalar"/>
  815   <param name="q" type="scalar"/>
  816   <param name="parameters" type="matrixref"/>
  817   <param name="DH" type="matrixref" optional="true"/>
  818   <param name="DE" type="matrixref" optional="true"/>
  819  </params>
  820 <code>ascore = exists(DH) &amp;&amp; exists(DE)
  821 if type &lt; 7 # aparch
  822   if ascore
  823     err = aparchFilter(&amp;depVar, &amp;h, &amp;e, &amp;mReg, &amp;vReg, p, q,   &amp;parameters, has_asymm_fx(type), &amp;DH, &amp;DE)
  824   else
  825     err = aparchFilter(&amp;depVar, &amp;h, &amp;e, &amp;mReg, &amp;vReg, p, q,   &amp;parameters, has_asymm_fx(type))
  826   endif
  827 elif type == 7
  828   err = egarchFilter(&amp;depVar, &amp;h, &amp;e, &amp;mReg, &amp;vReg, p, q, &amp;parameters)
  829 else
  830   err = 1
  831 endif
  832 return err
  833 </code>
  834 </gretl-function>
  835 <gretl-function name="do_score" type="matrix" private="1">
  836  <params count="5">
  837   <param name="de" type="seriesref"/>
  838   <param name="dh" type="seriesref"/>
  839   <param name="DE" type="matrixref"/>
  840   <param name="DH" type="matrixref"/>
  841   <param name="dd" type="matrixref"/>
  842  </params>
  843 <code>ret = {}
  844 matrix mde = misszero(de)
  845 matrix mdh = misszero(dh)
  846 #    printf &quot;rows(mde) = %d, rows(mdh) = %d\n&quot;, rows(mde), rows(mdh)
  847 #    printf &quot;rows(DE)  = %d, rows(DH)  = %d\n&quot;, rows(DE), rows(DH)
  848 #    printf &quot;%16.9f\n&quot;, mde[1:10] ~ DE[1:10,]
  849 #    printf &quot;%16.9f\n&quot;, mdh[1:10] ~ DH[1:10,]
  850 ret = (mde .* DE + mdh .* DH) ~ dd
  851 return ret
  852 </code>
  853 </gretl-function>
  854 <gretl-function name="do_mle" type="scalar" private="1">
  855  <params count="2">
  856   <param name="model" type="bundleref"/>
  857   <param name="verbose" type="bool"/>
  858  </params>
  859 <code>series depVar = model.y
  860 scalar type = model.type
  861 scalar cdist = model.cdist
  862 list mlX = model.mlistX
  863 list vlX = model.vlistX
  864 scalar p = model.p
  865 scalar q = model.q
  866 scalar nmX 	= model.mk
  867 scalar nvX 	= model.vk
  868 scalar err 	= 0
  869 series h, e, ll
  870 mleString = &quot;-&quot;
  871 if verbose
  872   string mleString += &quot;v&quot;
  873 else
  874   string mleString += &quot;q&quot;
  875 endif
  876 if model.vcvtype == 0 # robust
  877   string mleString += &quot;r&quot;
  878 elif model.vcvtype == 1 # hessian
  879   string mleString += &quot;h&quot;
  880   # otherwise opg
  881 endif
  882 if mleString == &quot;-&quot;
  883   mleString = &quot;&quot;
  884 endif
  885 fulltheta = model.coeff
  886 inipar = fulltheta
  887 sel = model.active
  888 theta = fulltheta[sel]
  889 filtpar_end = rows(fulltheta) - n_cdist_par(cdist)
  890 mX = model.mX
  891 if model.vX_QR == 1
  892   vX = model.QvX
  893   ini = model.mk+1
  894   fin = ini+model.vk-1
  895   theta[ini:fin] = model.vX_R * theta[ini:fin]
  896 else
  897   vX = model.vX
  898 endif
  899 set warnings off
  900 # experimental
  901 set bfgs_toler 1.0e-13
  902 matrix DH = {}
  903 matrix DE = {}
  904 series dh = NA
  905 series de = NA
  906 matrix dd = {}
  907 matrix score = NA
  908 if !ascore_ok(type, cdist)
  909   catch mle loglik = ll
  910     # put the newly estimated values into the full params vector
  911     fulltheta[sel] = theta
  912     filtpar = fulltheta[1:filtpar_end]
  913     matrix dpar = distpar(cdist, fulltheta)
  914     # FILTER
  915     err = gfilter(type, &amp;depVar, &amp;h, &amp;e, &amp;mX, &amp;vX, p, q, &amp;filtpar)
  916     ll = err ? NA : gig_loglik(&amp;e, &amp;h, cdist, dpar, &amp;de, &amp;dh, &amp;dd)
  917     params theta
  918   end mle @mleString
  919   err = $error
  920 else
  921   catch mle loglik = ll
  922     # put the newly estimated values into the full params vector
  923     fulltheta[sel] = theta
  924     filtpar = fulltheta[1:filtpar_end]
  925     matrix dpar = distpar(cdist, fulltheta)
  926     # FILTER
  927     err = gfilter(type, &amp;depVar, &amp;h, &amp;e, &amp;mX, &amp;vX, p, q, &amp;filtpar, &amp;DH, &amp;DE)
  928     ll = err ? NA : gig_loglik(&amp;e, &amp;h, cdist, dpar, &amp;de, &amp;dh, &amp;dd)
  929     # SCORE
  930     score = do_score(&amp;de, &amp;dh, &amp;DE, &amp;DH, &amp;dd)
  931     deriv theta = score
  932   end mle @mleString --no-gradient-check
  933   err = $error
  934 endif
  935 # err = gfilter(type, &amp;depVar, &amp;h, &amp;e, &amp;mX, &amp;vX, p, q, &amp;filtpar, &amp;DH, &amp;DE)
  936 # score = do_score(&amp;de, &amp;dh, &amp;DE, &amp;DH, &amp;dd)
  937 # printf &quot;Score:\n%20.10f\n&quot;, sumc(score)
  938 if (err==0)
  939   matrix crit = {$lnl; $aic; $bic; $hqc}
  940   matrix V = $vcv
  941 else
  942   matrix crit = zeros(4,1)
  943   npar = rows(theta)
  944   matrix V = zeros(npar, npar)
  945 endif
  946 gig_packResults(&amp;model, err, theta, &amp;h, &amp;e, inipar, V, crit)
  947 return err
  948 </code>
  949 </gretl-function>
  950 <gretl-function name="n_cdist_par" type="scalar" private="1">
  951  <params count="1">
  952   <param name="cdist" type="int"/>
  953  </params>
  954 <code># number of parameters for given density
  955 ret = NA
  956 if (cdist == 0) # normal
  957   ret = 0
  958 elif ((cdist == 1) || (cdist == 2)) # t and GED
  959   ret = 1
  960 elif ((cdist == 3) || (cdist == 4)) # Skewed t/GED
  961   ret = 2
  962 endif
  963 return ret
  964 </code>
  965 </gretl-function>
  966 <gretl-function name="cdist_initpar" type="matrix" private="1">
  967  <params count="1">
  968   <param name="cdist" type="int"/>
  969  </params>
  970 <code>/*
  971 initial parameters for given density; in the future, this may
  972 be data-based, although it's not trivial
  973 */
  974 ret = {}
  975 if (cdist == 0) # normal
  976   ret = {}
  977 elif (cdist == 1) # t
  978   ret = {10}
  979 elif (cdist == 2) # GED
  980   ret = {2}
  981 elif (cdist == 3) # Skewed t
  982   ret = {10; 0}
  983 elif (cdist == 4) # Skewed GED
  984   ret = {2; 0}
  985 endif
  986 return ret
  987 </code>
  988 </gretl-function>
  989 <gretl-function name="distpar" type="matrix" private="1">
  990  <params count="2">
  991   <param name="type" type="scalar"/>
  992   <param name="theta" type="matrix"/>
  993  </params>
  994 <code># extract the density parameter from the full vector
  995 ret = {}
  996 if (type&gt;2)
  997   ret = theta[rows(theta)-1:]
  998 elif (type&gt;0)
  999   ret = theta[rows(theta)]
 1000 endif
 1001 return ret
 1002 </code>
 1003 </gretl-function>
 1004 <gretl-function name="has_asymm_fx" type="scalar" private="1">
 1005  <params count="1">
 1006   <param name="code" type="int"/>
 1007  </params>
 1008 <code># does the model accommodate asymmetry (eg GJR)?
 1009 scalar ret = NA
 1010 if (code == 0) || (code == 1) || (code == 2) || (code == 5)
 1011   ret = 0
 1012 elif (code == 3) || (code == 4) || (code == 6) || (code == 7)
 1013   ret = 1
 1014 endif
 1015 return ret
 1016 </code>
 1017 </gretl-function>
 1018 <gretl-function name="ascore_ok" type="scalar" private="1">
 1019  <params count="2">
 1020   <param name="code" type="scalar"/>
 1021   <param name="cdist" type="scalar"/>
 1022  </params>
 1023 <code># do we have ascore in a usable state?
 1024 scalar mod_ok = 0
 1025 scalar dist_ok = 0
 1026 if (code == 0) || (code == 1) || (code == 3)
 1027   # other models may be ok at this point too,
 1028   # need to investigate
 1029   mod_ok = 1
 1030 endif
 1031 if (cdist&lt;2)
 1032   # only normal and t so far
 1033   dist_ok = 1
 1034 endif
 1035 force_ascore = ngetenv(&quot;ASCORE&quot;)
 1036 ret = ok(force_ascore) ? force_ascore : (mod_ok &amp;&amp; dist_ok)
 1037 return ret
 1038 </code>
 1039 </gretl-function>
 1040 <gretl-function name="gigDistString" type="string" private="1">
 1041  <params count="1">
 1042   <param name="code" type="int"/>
 1043  </params>
 1044 <code>string dist = &quot;&quot;
 1045 if code == 0
 1046   string garchName = &quot;Normal&quot;
 1047 elif code == 1
 1048   string garchName = &quot;T&quot;
 1049 elif code == 2
 1050   string garchName = &quot;GED&quot;
 1051 elif code == 3
 1052   string garchName = &quot;Skewed T&quot;
 1053 elif code == 4
 1054   string garchName = &quot;Skewed GED&quot;
 1055 endif
 1056 return dist
 1057 </code>
 1058 </gretl-function>
 1059 <gretl-function name="var_coef_init" type="matrix" private="1">
 1060  <params count="2">
 1061   <param name="mod" type="bundleref"/>
 1062   <param name="s2" type="scalar"/>
 1063  </params>
 1064 <code>/*
 1065 build variance parameters in this order:
 1066 1) omega + possible variance regressors
 1067 2) alphas
 1068 3) gammas
 1069 4) betas
 1070 5) delta
 1071 */
 1072 scalar code = mod.type
 1073 scalar p = mod.p
 1074 scalar q = mod.q
 1075 if (p &gt; 0)
 1076   scalar a = 0.1
 1077   scalar b = 0.8
 1078   scalar omega = s2 * (1 - a - b)
 1079 else
 1080   scalar a = 0.9
 1081   scalar omega = s2 * (1 - a)
 1082 endif
 1083 vk = mod.vk # variance regressors
 1084 matrix initv = {omega} | zeros(vk-1,1)
 1085 matrix alphas = (a/q) * ones(q,1)
 1086 initv |= alphas
 1087 initv |= zeros(q,1) #gammas
 1088 if (p &gt; 0)
 1089   initv |= {b} | zeros(p-1,1)
 1090 endif
 1091 # delta
 1092 if (code == 0) || (code == 1) || (code == 3)
 1093   delta = 2
 1094 elif (code == 2) || (code == 4)
 1095   delta = 1
 1096 elif (code == 5) || (code == 6)
 1097   delta = 1.5 # totally made up
 1098 else # eg EGARCH, though not really needed
 1099   delta = 2
 1100 endif
 1101 initv |= {delta}
 1102 return initv
 1103 </code>
 1104 </gretl-function>
 1105 <gretl-function name="coef_init" type="scalar" private="1">
 1106  <params count="1">
 1107   <param name="mod" type="bundleref"/>
 1108  </params>
 1109 <code>/* for the moment, naive; we ought to try hannan-rissanen on
 1110 the squared ols residuals
 1111 */
 1112 series dep = mod.y
 1113 if mod.mk &gt; 0
 1114   matrix e
 1115   matrix Reg = mod.mX
 1116   matrix initm = mols({dep}, Reg, &amp;e)
 1117   scalar s2 = meanc(e.^2)
 1118 else
 1119   matrix initm = {}
 1120   scalar s2 = var(dep)
 1121 endif
 1122 mod.s2 = s2
 1123 initv = var_coef_init(&amp;mod, s2)
 1124 cdist = mod.cdist
 1125 matrix initd = cdist_initpar(cdist)
 1126 mod.coeff = initm | initv | initd
 1127 return 0
 1128 </code>
 1129 </gretl-function>
 1130 <gretl-function name="set_active_par" type="scalar" private="1">
 1131  <params count="1">
 1132   <param name="mod" type="bundleref"/>
 1133  </params>
 1134 <code>scalar modeltype = mod.type
 1135 scalar mk = mod.mk
 1136 scalar vk = mod.vk
 1137 scalar p  = mod.p
 1138 scalar q  = mod.q
 1139 /* mean and variance regressors are always in
 1140 for any type of model
 1141 */
 1142 matrix act = ones(mk+vk, 1)
 1143 act |= ones(q, 1) # alphas
 1144 if has_asymm_fx(modeltype) # gammas
 1145   act |= ones(q, 1)
 1146 else
 1147   act |= zeros(q, 1)
 1148 endif
 1149 if (p&gt;0) # betas
 1150   act |= ones(p, 1)
 1151 endif
 1152 # delta is active only for NARCH and APARCH
 1153 act |= {(modeltype == 5) || (modeltype == 6)}
 1154 # density parameters
 1155 scalar cdist = mod.cdist
 1156 if (cdist == 1) || (cdist == 2)
 1157   act |= {1}
 1158 elif (cdist == 3)
 1159   act |= {1;1}
 1160 endif
 1161 n = rows(act)
 1162 mod.active = selifr(seq(1,n)', act)
 1163 return 0
 1164 </code>
 1165 </gretl-function>
 1166 <gretl-function name="vreg_QR" type="matrix" private="1">
 1167  <params count="2">
 1168   <param name="vX" type="matrix"/>
 1169   <param name="R" type="matrixref"/>
 1170  </params>
 1171 <code>/* do a QR decomp of the variance regressors
 1172 (seems to help)
 1173 */
 1174 scale = 10 * sqrt(rows(vX)) # totally heuristic
 1175 QvX = qrdecomp(vX, &amp;R)
 1176 matrix flip = selifr(seq(1,cols(vX))', (R[diag] .&lt; 0))
 1177 if rows(flip)&gt;0
 1178   QvX[,flip] = -QvX[,flip]
 1179   R[flip,] = -R[flip,]
 1180 endif
 1181 R = R ./ scale
 1182 return QvX .* scale
 1183 </code>
 1184 </gretl-function>
 1185 <gretl-function name="adjust_for_scale" type="void" private="1">
 1186  <params count="2">
 1187   <param name="m" type="bundleref"/>
 1188   <param name="scale" type="scalar"/>
 1189  </params>
 1190 <code>sc2 = scale^2
 1191 mk = m.mk
 1192 vk = m.vk
 1193 est_done = inbundle(m, &quot;vcv&quot;)
 1194 if est_done
 1195   # we have two dimensions here: &quot;coeff&quot; contains all parameters,
 1196   # &quot;vcv&quot; and &quot;stderr&quot; only active ones; shit, what was I thinking?
 1197   J1 = ones(rows(m.coeff),1)
 1198   J2 = ones(rows(m.vcv),1)
 1199   if mk&gt;0
 1200     J1[1:mk] = 1/scale
 1201     J2[1:mk] = 1/scale
 1202   endif
 1203   J1[mk+1:mk+vk] = 1/sc2
 1204   J2[mk+1:mk+vk] = 1/sc2
 1205 endif
 1206 if inbundle(m, &quot;s2&quot;)
 1207   m.s2 = m.s2 / sc2
 1208 endif
 1209 if est_done
 1210   # FIXME: what happens to the variance equation if
 1211   # delta != 2?
 1212   m.coeff = m.coeff .* J1
 1213   m.stderr = m.stderr .* J2
 1214   m.vcv = m.vcv .* (J2*J2')
 1215   m.uhat = m.uhat / scale
 1216   m.h = m.h / sc2
 1217   crit_adj = ln(scale)*m.nobs .* {1; -2; -2; -2}
 1218   m.criteria = m.criteria + crit_adj
 1219 endif
 1220 </code>
 1221 </gretl-function>
 1222 <gretl-function name="gig_packResults" type="void" private="1">
 1223  <params count="8">
 1224   <param name="mod" type="bundleref"/>
 1225   <param name="err" type="scalar"/>
 1226   <param name="thetahat" type="matrix"/>
 1227   <param name="h" type="seriesref"/>
 1228   <param name="e" type="seriesref"/>
 1229   <param name="inipar" type="matrix"/>
 1230   <param name="Sigma" type="matrix"/>
 1231   <param name="crit" type="matrix"/>
 1232  </params>
 1233 <code>mod.errcode = err
 1234 if (err == 0)
 1235   theta = mod.coeff
 1236   sel = mod.active
 1237   theta[sel] = thetahat
 1238   mod.inipar = inipar
 1239   vk = mod.vk
 1240   if (vk&gt;1) &amp;&amp; (mod.vX_QR == 1)
 1241     J = I(rows(Sigma))
 1242     mk = mod.mk
 1243     sel = seq(mk+1, mk+vk)
 1244     invR = inv(mod.vX_R)
 1245     J[sel,sel] = invR
 1246     theta[sel] = invR * theta[sel]
 1247     Sigma = qform(J, Sigma)
 1248   endif
 1249   mod.coeff = theta
 1250   mod.vcv = Sigma
 1251   mod.stderr = sqrt(diag(Sigma))
 1252   mod.h = h
 1253   mod.uhat = e
 1254   mod.stduhat = e/sqrt(h)
 1255   mod.criteria = crit
 1256 endif
 1257 </code>
 1258 </gretl-function>
 1259 <gretl-function name="s_modeltype" type="string" private="1">
 1260  <params count="1">
 1261   <param name="mod" type="bundleref"/>
 1262  </params>
 1263 <code>scalar type = mod.type
 1264 scalar p    = mod.p
 1265 scalar q    = mod.q
 1266 if type ==  0
 1267   return sprintf(&quot;ARCH(%d) [Engle]&quot;, p)
 1268 elif type ==  1
 1269   return sprintf(&quot;GARCH(%d,%d) [Bollerslev]&quot;, p, q)
 1270 elif type ==  2
 1271   return sprintf(&quot;Taylor/Schwert's GARCH(%d,%d)&quot;, p, q)
 1272 elif type ==  3
 1273   return sprintf(&quot;GJR(%d,%d) [Glosten et al.]&quot;, p, q)
 1274 elif type ==  4
 1275   return sprintf(&quot;TARCH(%d,%d) [Zakoian]&quot;, p, q)
 1276 elif type ==  5
 1277   return sprintf(&quot;NARCH(%d,%d) [Higgins and Bera]&quot;, p, q)
 1278 elif type ==  6
 1279   return sprintf(&quot;APARCH(%d,%d) [Ding]&quot;, p, q)
 1280 elif type ==  7
 1281   return sprintf(&quot;EGARCH(%d,%d) [Nelson]&quot;, p, q)
 1282 elif type ==  8
 1283   return sprintf(&quot;GARCH-in-mean(%d,%d) [Engle]&quot;, p, q)
 1284 elif type ==  9
 1285   return sprintf(&quot;free-APARCH(%d,%d)&quot;, p, q)
 1286 endif
 1287 </code>
 1288 </gretl-function>
 1289 <gretl-function name="gig_printHeader" type="void" private="1">
 1290  <params count="1">
 1291   <param name="mod" type="bundleref"/>
 1292  </params>
 1293 <code>ncoeff  = rows(mod.coeff) # rows
 1294 nfitted = rows(mod.active) # cols
 1295 type  = mod.type
 1296 dtype = mod.cdist
 1297 vtype = mod.vcvtype
 1298 p    = mod.p
 1299 q    = mod.q
 1300 t1   = mod.t1
 1301 t2   = mod.t2
 1302 T    = mod.nobs
 1303 printf &quot;\nModel: %s&quot;, s_modeltype(&amp;mod)
 1304 # DISTRIBUTIONS ordered by no. of extra params
 1305 if dtype ==  0
 1306   printf &quot; (Normal)&quot;
 1307 elif dtype ==  1
 1308   printf &quot; (Student's t)&quot;
 1309 elif dtype ==  2
 1310   printf &quot; (GED)&quot;
 1311 elif dtype ==  3
 1312   printf &quot; (Skewed T)&quot;
 1313 elif dtype ==  4
 1314   printf &quot; (Skewed GED)&quot;
 1315 endif
 1316 if ascore_ok(type, dtype)
 1317   printf &quot;*&quot;
 1318 endif
 1319 # OBS used
 1320 printf &quot;\nDependent variable: %s&quot;, mod.depvarname
 1321 printf &quot;\nSample: %s -- %s (T = %d)&quot;, obslabel(t1), obslabel(t2), T
 1322 # VCV
 1323 printf &quot;, VCV method: &quot;
 1324 if vtype ==  0
 1325   printf &quot;Robust&quot;
 1326 elif vtype ==  1
 1327   printf &quot;Hessian&quot;
 1328 else
 1329   printf &quot;OPG&quot;
 1330 endif
 1331 if (mod.vX_QR == 1) &amp;&amp; (mod.vk&gt;1)
 1332   printf &quot;\nQR decomposition used for variance regressors&quot;
 1333 endif
 1334 # if (mod.scale != 1)
 1335 #     printf &quot;\nexperimental: scaling = %g\n&quot;, mod.scale
 1336 # endif
 1337 printf &quot;\n\n&quot;
 1338 </code>
 1339 </gretl-function>
 1340 <gretl-function name="GJR_alt_param" type="void" private="1">
 1341  <params count="7">
 1342   <param name="coeff" type="matrix"/>
 1343   <param name="vcv" type="matrix"/>
 1344   <param name="p" type="scalar"/>
 1345   <param name="q" type="scalar"/>
 1346   <param name="nmX" type="scalar"/>
 1347   <param name="nvX" type="scalar"/>
 1348   <param name="svX" type="string"/>
 1349  </params>
 1350 <code>matrix c_om = coeff[nmX+1:nmX+nvX]
 1351 matrix v_om = vcv[nmX+1:nmX+nvX, nmX+1:nmX+nvX]
 1352 matrix c_ag = coeff[nmX+nvX+1:nmX+nvX+2*q]
 1353 matrix v_ag = vcv[nmX+nvX+1:nmX+nvX+2*q, nmX+nvX+1:nmX+nvX+2*q]
 1354 matrix c_bt = coeff[nmX+nvX+2*q+1:nmX+nvX+2*q+p]
 1355 matrix v_bt = vcv[nmX+nvX+2*q+1:nmX+nvX+2*q+p, nmX+nvX+2*q+1:nmX+nvX+2*q+p]
 1356 if nvX == 1
 1357   parNames = &quot;delta&quot;
 1358 else
 1359   parNames = svX
 1360 endif
 1361 cs = c_om ~ sqrt(diag(v_om))
 1362 if q == 1
 1363   parNames += &quot;,alpha,gamma&quot;
 1364 else
 1365   loop for i=1..q --quiet
 1366     parNames += &quot;,alpha_$i&quot;
 1367   endloop
 1368   loop for i=1..q --quiet
 1369     parNames += &quot;,gamma_$i&quot;
 1370   endloop
 1371 endif
 1372 c_ag2 = zeros(2*q,1)
 1373 J = zeros(2*q,2*q) # Jacobian
 1374 loop i = 1 .. q
 1375   alpha = c_ag[i]
 1376   gamma = c_ag[i+q]
 1377   c_ag2[i] = alpha * (1-gamma)^2
 1378   c_ag2[i+q] = 4 * alpha * gamma
 1379   J[i,i] = (1-gamma)^2
 1380   J[i,i+q] = -2 * alpha * (1-gamma)
 1381   J[i+q,i] = 4 * gamma
 1382   J[i+q,i+q] = 4 * alpha
 1383 endloop
 1384 v_ag2 = qform(J, v_ag)
 1385 cs |= (c_ag2 ~ sqrt(diag(v_ag2)))
 1386 if p == 1
 1387   parNames += &quot;,beta&quot;
 1388 else
 1389   loop for i=1..p --quiet
 1390     parNames += &quot;,beta_$i&quot;
 1391   endloop
 1392 endif
 1393 cs |= (c_bt ~ sqrt(diag(v_bt)))
 1394 printf &quot;\n   (alt. parametrization)\n&quot;
 1395 modprint cs parNames
 1396 </code>
 1397 </gretl-function>
 1398 <gretl-function name="maybe_gig_print" type="void" private="1">
 1399  <params count="2">
 1400   <param name="mod" type="bundleref"/>
 1401   <param name="verbose" type="scalar" min="0" max="3" default="0"/>
 1402  </params>
 1403 <code>ncoeff  = rows(mod.coeff) # rows
 1404 nfitted = rows(mod.active) # cols
 1405 type  = mod.type
 1406 dtype = mod.cdist
 1407 vtype = mod.vcvtype
 1408 p    = mod.p
 1409 q    = mod.q
 1410 nmX  = mod.mk
 1411 nvX  = mod.vk
 1412 crit = mod.criteria
 1413 vcv = mod.vcv
 1414 # Composing print-out results
 1415 coeff  = mod.coeff
 1416 sel    = mod.active
 1417 stderr = mod.stderr
 1418 cs2 = coeff[sel] ~ stderr
 1419 # MEAN REGR
 1420 if nmX &gt; 0
 1421   parNames = mod.mXnames
 1422   cs = cs2[1:nmX,]
 1423   if cols(cs)&gt;0
 1424     printf &quot;    Conditional mean equation\n&quot;
 1425     modprint cs parNames
 1426   endif
 1427 endif
 1428 # VAR REGR
 1429 if nvX == 1
 1430   parNames = &quot;omega&quot;
 1431 else
 1432   parNames = mod.vXnames
 1433 endif
 1434 #ALPHAS
 1435 counter = nmX+nvX
 1436 if q == 1
 1437   parNames += &quot;,alpha&quot;
 1438 else
 1439   loop for i=1..q --quiet
 1440     parNames += &quot;,alpha_$i&quot;
 1441   endloop
 1442 endif
 1443 counter += q
 1444 #GAMMAS
 1445 if has_asymm_fx(type)
 1446   if q == 1
 1447     parNames += &quot;,gamma&quot;
 1448   else
 1449     loop for i=1..q --quiet
 1450       parNames += &quot;,gamma_$i&quot;
 1451     endloop
 1452   endif
 1453 endif
 1454 counter += q
 1455 #BETAS
 1456 if type&gt;0
 1457   if p == 1
 1458     parNames += &quot;,beta&quot;
 1459   else
 1460     loop for i=1..p --quiet
 1461       parNames += &quot;,beta_$i&quot;
 1462     endloop
 1463   endif
 1464 endif
 1465 counter += p
 1466 #DELTAS
 1467 if (type==5) || (type==6)
 1468   parNames += &quot;,delta&quot;
 1469 endif
 1470 if type != 9
 1471   cs = cs2[nmX+1:nfitted,]
 1472   limit = nfitted-nmX
 1473 else
 1474   cs = cs2[nmX+1:,]
 1475   limit = ncoeff-nmX
 1476 endif
 1477 printf &quot;    Conditional variance equation\n&quot;
 1478 if dtype &gt; 0
 1479   ndenspar = n_cdist_par(dtype)
 1480   csdens = cs[limit-ndenspar+1:limit,]
 1481   cs = cs[1:limit-ndenspar,]
 1482 endif
 1483 modprint cs parNames
 1484 if (type == 3) # GJR
 1485   GJR_alt_param(coeff, vcv, p, q, nmX, nvX, mod.vXnames)
 1486 endif
 1487 if dtype &gt; 0
 1488   # DENSITY FUNCTION PART
 1489   /*
 1490   FIXME (or maybe not):
 1491   for asymmetric densities what we're actually printing here
 1492   is not really \lambda, but rather atanh(\lambda). In most cases, the difference is so minuscule to be inconsequential. Besides, I don't think anyone is going to care. However, it's annoying.
 1493   What shall we do? Transform the estimated parameter (and vcv) into
 1494   \lambda or just correct its label? Note that when we plot the
 1495   density (see gig_plot.inp) we do transform the parameter back to
 1496   \lambda, so we should at least be consistent. Hmmm.
 1497   */
 1498   if dtype&lt;3
 1499     parNames = &quot;,ni&quot;
 1500   elif dtype&lt;5
 1501     parNames = &quot;ni,lambda&quot;
 1502   endif
 1503   printf &quot;    Conditional density parameters\n&quot;
 1504   modprint csdens parNames
 1505 endif
 1506 # INFO
 1507 printf &quot;\tLlik: %12.5f\t&quot;, crit[1]
 1508 printf &quot; AIC: %12.5f\n&quot;,   crit[2]
 1509 printf &quot;\tBIC:  %12.5f\t&quot;, crit[3]
 1510 printf &quot; HQC: %12.5f\n\n&quot;, crit[4]
 1511 if verbose &gt; 2
 1512   # Starting Values
 1513   printf &quot;Starting/Fitted Values Comparison Matrix:\n\n&quot;
 1514   printf &quot;       Initial         Final    Difference\n\n&quot;
 1515   fitted  = mod.coeff
 1516   initial = mod.inipar
 1517   loop for i=1..ncoeff -q
 1518     difff = initial[i] - fitted[i]
 1519     printf &quot;%14.7f%14.7f%14.7f\n&quot;, initial[i], fitted[i], difff
 1520   endloop
 1521   printf &quot;\n&quot;
 1522 endif
 1523 </code>
 1524 </gretl-function>
 1525 <gretl-function name="dateton" type="scalar" private="1">
 1526  <params count="1">
 1527   <param name="t" type="scalar"/>
 1528  </params>
 1529 <code>scalar ret dd mm yy
 1530 string s = obslabel(t)
 1531 scalar pd = $pd
 1532 if $version&lt;10913
 1533   string datefmt = &quot;%d/%d/%d&quot;
 1534 else
 1535   string datefmt = &quot;%d-%d-%d&quot;
 1536 endif
 1537 if (pd==5) || (pd==6) || (pd==7)
 1538   sscanf(s,datefmt, yy, mm, dd)
 1539   ret = yy + ok(mm)*(mm-1)/12 + ok(dd)*(dd-1)*12/365
 1540 elif (pd==12) || (pd==4)
 1541   sscanf(s,&quot;%d:%d&quot;, yy, mm)
 1542   ret = yy + ok(mm)*(mm-1)/pd
 1543 else
 1544   ret = t
 1545 endif
 1546 return ret
 1547 </code>
 1548 </gretl-function>
 1549 <gretl-function name="do_xtics" type="string" private="1">
 1550  <params count="5">
 1551   <param name="t1" type="scalar"/>
 1552   <param name="t2" type="scalar"/>
 1553   <param name="d1" type="scalar"/>
 1554   <param name="d2" type="scalar"/>
 1555   <param name="n" type="scalar"/>
 1556  </params>
 1557 <code>string ret = &quot;set xtics(&quot;
 1558 s = seq(1,n)/n - 1/(2*n)
 1559 s_d = d1 + s .* (d2 - d1)
 1560 t = t1 + (t2 - t1)/(2*n)
 1561 loop for i = 1 .. n
 1562   string tmp = sprintf(&quot;'%s' %g&quot;, obslabel(round(t)), s_d[i])
 1563   ret += tmp
 1564   if i==n
 1565     ret += &quot;)&quot;
 1566   else
 1567     ret += &quot;, &quot;
 1568   endif
 1569   t += (t2 - t1)/n
 1570 endloop
 1571 return ret
 1572 </code>
 1573 </gretl-function>
 1574 <gretl-function name="gig_grph" type="string" private="1">
 1575  <params count="3">
 1576   <param name="yname" type="string"/>
 1577   <param name="u" type="series"/>
 1578   <param name="h" type="series"/>
 1579  </params>
 1580 <code>series s = sqrt(h)
 1581 smpl ok(s) --restrict
 1582 scalar tmpfname = int(1000*muniform(1,1))
 1583 if $windows
 1584   tmpfile=sprintf(&quot;@dotdir\\gig%d.gp&quot;, tmpfname)
 1585 else
 1586   tmpfile=sprintf(&quot;@dotdir/gig%d.gp&quot;, tmpfname)
 1587 endif
 1588 scalar t1 = $t1
 1589 scalar t2 = $t2
 1590 ### FIXME 10/11/2020: this can probably be made much smarter by using
 1591 ### more recent calendar functions
 1592 scalar d1 = dateton(t1)
 1593 scalar d2 = dateton(t2)
 1594 scalar incr = (d2-d1)/($nobs-1)
 1595 series xt = d1 + incr*(time - min(time))
 1596 set force_decpoint on
 1597 outfile &quot;@tmpfile&quot; --write
 1598 printf &quot;set nokey\n&quot;
 1599 #    printf &quot;unset xzeroaxis\n&quot;
 1600 printf &quot;set title '%s: Residuals and conditional sd'\n&quot;, yname
 1601 printf &quot;%s\n&quot;, do_xtics(t1, t2, d1, d2, 6)
 1602 printf &quot;plot '-' using 1:2 title 'residual' w lines, \\\n&quot;
 1603 printf &quot;'-' using 1:2 w lines lt 2, \\\n&quot;
 1604 printf &quot;'-' using 1:2 w lines lt 2 \n&quot;
 1605 loop for t = t1 .. t2
 1606   printf &quot;%12.6f %12.6f\n&quot;, xt[t], u[t]
 1607 endloop
 1608 printf &quot;e\n&quot;
 1609 loop for t = t1 .. t2
 1610   printf &quot;%12.6f %12.6f\n&quot;, xt[t], s[t]
 1611 endloop
 1612 printf &quot;e\n&quot;
 1613 loop for t = t1 .. t2
 1614   printf &quot;%12.6f %12.6f\n&quot;, xt[t], -s[t]
 1615 endloop
 1616 printf &quot;e\n&quot;
 1617 outfile --close
 1618 set force_decpoint off
 1619 return tmpfile
 1620 </code>
 1621 </gretl-function>
 1622 <gretl-function name="gig_dgrph" type="string" private="1">
 1623  <params count="2">
 1624   <param name="X" type="matrix"/>
 1625   <param name="desc" type="string"/>
 1626  </params>
 1627 <code>n = rows(X)
 1628 x  = X[,1]
 1629 kd = X[,2]
 1630 td = X[,3]
 1631 scalar tmpfname = int(1000*muniform(1,1))
 1632 if $windows
 1633   tmpfile = sprintf(&quot;@dotdir\\gig%d.gp&quot;, tmpfname)
 1634 else
 1635   tmpfile = sprintf(&quot;@dotdir/gig%d.gp&quot;, tmpfname)
 1636 endif
 1637 set force_decpoint on
 1638 outfile &quot;@tmpfile&quot; --write
 1639 printf &quot;set nokey\n&quot;
 1640 printf &quot;set yzeroaxis\n&quot;
 1641 printf &quot;set title 'std residuals: kernel density vs %s'\n&quot;, desc
 1642 printf &quot;plot '-' using 1:2 title 'kernel' w lines, \\\n&quot;
 1643 printf &quot;'-' using 1:2 title 'theoretical' w lines lt 2 \n&quot;
 1644 loop for i = 1 .. n
 1645   printf &quot;%12.6f %12.6f\n&quot;, x[i], kd[i]
 1646 endloop
 1647 printf &quot;e\n&quot;
 1648 loop for i = 1 .. n
 1649   printf &quot;%12.6f %12.6f\n&quot;, x[i], td[i]
 1650 endloop
 1651 printf &quot;e\n&quot;
 1652 outfile --close
 1653 set force_decpoint off
 1654 return tmpfile
 1655 </code>
 1656 </gretl-function>
 1657 <gretl-function name="pdf_skt" type="matrix" private="1">
 1658  <params count="3">
 1659   <param name="e" type="matrixref"/>
 1660   <param name="df" type="scalar"/>
 1661   <param name="ht_skew" type="scalar"/>
 1662  </params>
 1663 <code>matrix ret = NA
 1664 if (df&gt;2)
 1665   # sqrt(pi) = 1.77245385090551602729816748334
 1666   scalar q = lngamma((df+1)/2) - lngamma(df/2)
 1667   scalar c = exp(q)/(sqrt(df-2)*1.77245385090551602729816748334)
 1668   scalar a = 4 * ht_skew * c * ((df-2)/(df-1))
 1669   scalar b = sqrt(1 + 3*ht_skew^2 - a^2)
 1670   matrix d = (b.* e + a)
 1671   sk = (d.&lt;0) .* (1-ht_skew) + (d.&gt;0) .* (1+ht_skew)
 1672   d = d ./ sk
 1673   ret = b*c / exp( ((df+1)/2) * ln(1+(d.^2/(df-2))) )
 1674 endif
 1675 return ret
 1676 </code>
 1677 </gretl-function>
 1678 <gretl-function name="pdf_skged" type="matrix" private="1">
 1679  <params count="3">
 1680   <param name="x" type="matrixref"/>
 1681   <param name="ni" type="scalar"/>
 1682   <param name="ta" type="scalar"/>
 1683  </params>
 1684 <code>scalar p  = 1/ni
 1685 lgp  = lngamma(p)
 1686 lg2p = lngamma(2*p)
 1687 lg3p = lngamma(3*p)
 1688 tap1 = 1 + ta
 1689 tam1 = 1 - ta
 1690 scalar beta = 0.5 * exp(lg3p - lgp) * (tap1^3 + tam1^3) - 4*ta^2 * exp(2 * (lg2p - lgp))
 1691 beta = sqrt(beta)
 1692 # m: mode
 1693 scalar m = - 2*ta/beta * exp( lg2p - lgp )
 1694 scalar lnorm = ln(0.5 * beta) - lngamma(p+1)
 1695 matrix sk = (x.&lt;m) .* (-tam1) + (x.&gt;m) .* tap1
 1696 matrix z = (x-m) ./ sk
 1697 matrix ret = lnorm - (z.^ni)
 1698 return exp(ret)
 1699 </code>
 1700 </gretl-function>
 1701 <gretl-function name="aparch_fcast" type="matrix" private="1">
 1702  <params count="8">
 1703   <param name="e" type="matrix"/>
 1704   <param name="e0" type="matrix"/>
 1705   <param name="h0" type="matrix"/>
 1706   <param name="omega" type="scalar"/>
 1707   <param name="alpha" type="matrix"/>
 1708   <param name="gamma" type="matrix"/>
 1709   <param name="beta" type="matrix"/>
 1710   <param name="delta" type="scalar"/>
 1711  </params>
 1712 <code># in principle, we could just as well use the aparchFilter
 1713 # function from gig_mle.inp, but this is lighter and more specialised
 1714 q = rows(alpha)
 1715 p = rows(beta)
 1716 head = xmax(p,q)
 1717 scalar hor = rows(e) + 1
 1718 scalar is_asymmetric = maxc(abs(gamma)) .&gt; 1.0e-07
 1719 matrix ee  = e0 | e
 1720 matrix ae = abs(ee)
 1721 h0init = h0.^(delta/2)
 1722 scalar s = head + 1
 1723 matrix elag = mlag(ee, seq(0,q-1))
 1724 if is_asymmetric
 1725   tmp_ae = abs(elag) - elag .* gamma'
 1726 else
 1727   tmp_ae = abs(elag)
 1728 endif
 1729 # print e e0 tmp_ae h0init
 1730 matrix h = omega + (tmp_ae.^delta) * alpha
 1731 matrix ret = p&gt;0 ? filter(h0init | h, 1, beta) : h
 1732 scalar good = rows(ret) - hor
 1733 ret = ret[good+1:]
 1734 if delta != 2
 1735   ret = ret.^(2/delta)
 1736 endif
 1737 return ret
 1738 </code>
 1739 </gretl-function>
 1740 <gretl-function name="s2m" type="matrix" private="1">
 1741  <params count="1">
 1742   <param name="x" type="series"/>
 1743  </params>
 1744 <code>matrix a = ok(x)
 1745 matrix b = misszero(x)
 1746 return selifr(b,a)
 1747 </code>
 1748 </gretl-function>
 1749 <sample-script filename="examples/example1.inp">
 1750 # Example n. 1
 1751 
 1752 # GARCH(1,1) with Normal conditional distribution, 
 1753 # no mean regressors.
 1754   
 1755 
 1756 set verbose off
 1757 
 1758 # Import the gig library.
 1759 include gig.gfn
 1760 open djclose
 1761 rr = 100*ldiff(djclose)
 1762 
 1763 model = gig_setup(rr)
 1764 gig_estimate(&amp;model)
 1765 
 1766 # Compare with native implementation
 1767 garch 1 1 ; rr --nc --robust
 1768 </sample-script>
 1769 </gretl-function-package>
 1770 </gretl-functions>