"Fossies" - the Fresh Open Source Software Archive

Member "gretl/functions/extra/extra.gfn" (21 Nov 2020, 30030 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="extra" no-data-ok="true" minver="2016a" lives-in-subdir="true">
    4 <author email="&lt;use the mailing list or trackers&gt;">The gretl team</author>
    5 <version>0.8</version>
    6 <date>2020-09-19</date>
    7 <description>extra functions for hansl scripting</description>
    8 <tags>C88</tags>
    9 <help>
   10 pdfdoc:extra.pdf
   11 </help>
   12 <gretl-function name="gap_filler" type="series">
   13  <params count="2">
   14   <param name="x" type="series">
   15 <description>Series to fill</description>
   16   </param>
   17   <param name="method" type="int" min="0" max="2" default="2">
   18 <description>Fill method</description>
   19 <labels count="3">
   20 "Don't" "Repetition" "Linear interpolation" </labels>
   21   </param>
   22  </params>
   23 <code>string lbl = argname(x)
   24 if method == 0 # do nothing
   25   series ret = x
   26 elif method == 1 # get previous value
   27   genr time
   28   series OK = ok(x)
   29   series ret = x
   30   series tmp = OK ? time : NA
   31   scalar ini = min(tmp)
   32   scalar fin = max(tmp)
   33   smpl ini fin
   34   ret = OK ? x : ret(-1)
   35   string lbl = sprintf(&quot;gap-filled version of %s (with repetition)&quot;, argname(x))
   36   setinfo ret --description=&quot;@lbl&quot;
   37 elif method == 2 # interpolate_linearly
   38   set warnings off
   39   series ret = lin_int(x)
   40   string lbl = sprintf(&quot;gap-filled version of %s (with interpolation)&quot;, argname(x))
   41   setinfo ret --description=&quot;@lbl&quot;
   42 endif
   43 return ret
   44 </code>
   45 </gretl-function>
   46 <gretl-function name="winsor" type="series">
   47  <params count="3">
   48   <param name="x" type="series"/>
   49   <param name="p" type="scalar" min="0" max="1" default="0.05"/>
   50   <param name="phi" type="scalar" min="0" max="1" default="0"/>
   51  </params>
   52 <code># phi (as in p-high) is optional and defines an upper tail mass
   53 #  different from lower (the default phi == 0 means ignore)
   54 #
   55 # This is a rewrite of the function in the winsor.gfn package
   56 # (originally by JoshuaHe2015@163.com)
   57 smpl --no-missing x
   58 # standard symmetric or asymmetric case?
   59 phi = !phi ? 1 - p : phi
   60 # determine bounds
   61 matrix lowhi = quantile({x}, {p, phi})
   62 # lower end (and address non-existing extreme quantiles)
   63 scalar low = ok(lowhi[1]) ? lowhi[1] : min(x)
   64 x = (x &lt; low) ? low : x
   65 # upper end
   66 scalar hi = ok(lowhi[2]) ? lowhi[2] : max(x)
   67 x = (x &gt; hi) ? hi : x
   68 # prepare output
   69 string label = sprintf(&quot;winsorized %s (%g,%g)&quot;, argname(x), p, phi)
   70 setinfo x --description=&quot;@label&quot;
   71 return x
   72 </code>
   73 </gretl-function>
   74 <gretl-function name="nearPSD" type="scalar">
   75  <params count="2">
   76   <param name="m" type="matrixref"/>
   77   <param name="epsilon" type="scalar" min="0" default="0"/>
   78  </params>
   79 <code># TODO: cover the case with some variances == 0
   80 # (use misszero() after scaling, or something)
   81 #
   82 # Forces the matrix m into the positive semi-definite region.
   83 #
   84 # Ported from &quot;DomPazz&quot; in Stackoverflow, apparently
   85 # mimicking the nearPD() function in R.
   86 # Because of re-scaling ( to correlation matrix ), the
   87 # epsilon should implicitly apply to the correlation-based
   88 # eigenvalues.
   89 #
   90 # The return value 0 or 1 indicates whether m was altered or not.
   91 matrix s = sqrt(diag(m)) # std dev
   92 matrix scaling = s * s'
   93 matrix ms = m ./ scaling	# rescaled
   94 matrix eigvec
   95 matrix eigval = eigensym(ms, &amp;eigvec)
   96 matrix val = (eigval .&gt; epsilon) ? eigval : epsilon	# like xmax
   97 if sum(val .&gt; eigval)	# actually something was changed
   98   matrix T = 1 / ( (eigvec .^2) * val )
   99   # transform vector T to a diagonal matrix
  100   matrix temp = I(rows(T))
  101   temp[diag] = sqrt(T)
  102   # and also val
  103   matrix temp2 = I(rows(T))
  104   temp2[diag] = sqrt(val)
  105   matrix B = temp * eigvec * temp2
  106   ms = B * B'
  107   # undo the re-scaling
  108   m = ms .* scaling
  109   return 1
  110 else
  111   return 0
  112 endif
  113 </code>
  114 </gretl-function>
  115 <gretl-function name="zeroifclose" type="scalar">
  116  <params count="2">
  117   <param name="m" type="matrixref"/>
  118   <param name="thresh" type="scalar" min="0" default="1e-12"/>
  119  </params>
  120 <code># Sets elements to zero if they are really close.
  121 # The return value 0 or 1 indicates whether m was altered or not.
  122 # (an older version copied and returned the matrix)
  123 matrix indicator = (abs(m) .&lt; thresh)
  124 ret = sum(indicator) &gt; 0
  125 if ret
  126   m = indicator ? 0 : m
  127 endif
  128 return ret
  129 </code>
  130 </gretl-function>
  131 <gretl-function name="commute" type="matrix">
  132  <params count="4">
  133   <param name="A" type="matrix" const="true"/>
  134   <param name="m" type="int" min="1"/>
  135   <param name="n" type="int" min="0" default="0"/>
  136   <param name="post" type="bool" default="0"/>
  137  </params>
  138 <code># returns A premultiplied by K_mn,
  139 # K_mn being the commutation matrix, see Magnus/Neudecker
  140 # (more efficient than explicit pre-multiplication,
  141 # just reshuffles stuff)
  142 # If post != 0, does post-multiplication:
  143 # if n == 0 (default), then it's understood to be equal to m,
  144 # as per Magnus/Neudecker convention: K_nn = K_n
  145 n = n ? n : m
  146 if !post &amp;&amp; (rows(A) != n * m)
  147   funcerr &quot;m*n doesn't match row dim of input&quot;
  148 elif post &amp;&amp; (cols(A) != n * m)
  149   funcerr &quot;m*n doesn't match column dimension of input&quot;
  150 endif
  151 matrix e = vec(mshape(seq(1, m*n), m, n)')
  152 return post ? A[,e] : A[e,]
  153 </code>
  154 </gretl-function>
  155 <gretl-function name="eliminate" type="matrix">
  156  <params count="1">
  157   <param name="vecA" type="matrix" const="true"/>
  158  </params>
  159 <code># Each column of the input vecA is assumed to
  160 # come from the operation vec(A) on a square
  161 # matrix, thus rows(vecA) must be a square number.
  162 # Returns vech(A), which is the result of pre-
  163 # multiplying vec(A) with the &quot;elimination&quot;
  164 # matrix L_m.
  165 # If vecA has several columns, each column is
  166 # treated separately as described above
  167 # (and the results stacked side-by-side).
  168 r = sqrt(rows(vecA))
  169 if r != round(r)
  170   funcerr &quot;Input must have a square number of rows&quot;
  171 endif
  172 e = vech(mshape(seq(1, r^2), r, r)')
  173 return vecA[e,]
  174 </code>
  175 </gretl-function>
  176 <gretl-function name="duplicate" type="matrix">
  177  <params count="1">
  178   <param name="vechA" type="matrix" const="true"/>
  179  </params>
  180 <code># The input is a vector assumed to come from
  181 # an operation like vech(A).
  182 # Returns vec(A), which is the result of pre-
  183 # multiplying vech(A) with the &quot;duplication&quot;
  184 # matrix D_m.
  185 # If vechA has several columns, each column is
  186 # treated separately as described above
  187 # (and the results stacked side-by-side).
  188 e = vec(unvech(seq(1,rows(vechA))'))
  189 return vechA[e,]
  190 </code>
  191 </gretl-function>
  192 <gretl-function name="truncnorm" type="matrix">
  193  <params count="5">
  194   <param name="n" type="int" min="1">
  195 <description>Number of random variables</description>
  196   </param>
  197   <param name="m" type="scalar" default="0" const="true">
  198 <description>Mean value</description>
  199   </param>
  200   <param name="sigma" type="scalar" min="0" default="1" const="true">
  201 <description>Std. deviation</description>
  202   </param>
  203   <param name="below" type="scalar" const="true">
  204 <description>left truncation value</description>
  205   </param>
  206   <param name="above" type="scalar" const="true">
  207 <description>right truncation value</description>
  208   </param>
  209  </params>
  210 <code>/* Generates truncated normal random values. Set 'below' and/or 'above'
  211 to NA to skip. Returns a column vector. */
  212 if !ok(m)
  213   funcerr &quot;You must pass a valid (non-nan) value for the mean.&quot;
  214 endif
  215 scalar l = ok(below) ? cnorm((below - m)/sigma) : 0
  216 scalar r = ok(above) ? cnorm((above - m)/sigma) : 1
  217 matrix u = l + (r - l) .* muniform(n, 1)
  218 return invcdf(z, u) .* sigma + m
  219 </code>
  220 </gretl-function>
  221 <gretl-function name="scores2x2" type="matrix">
  222  <params count="2">
  223   <param name="in" type="matrix"/>
  224   <param name="verbose" type="bool" default="0"/>
  225  </params>
  226 <code>/*
  227 Computes some standard score measures for a 2x2
  228 contingency table of the form:
  229 Observed
  230 1      0
  231 --------------
  232 Predicted  1 | h(its)  f(alse)
  233 0 | m(iss)  z(eros)
  234 and n = h + f + m + z (total observations)
  235 1) POD / prob of detection = h / (h + m)
  236 2) POFD / prob of false detection = f / (f + z)
  237 3) HR / hit rate = (h + z) / n
  238 4) FAR / false alarm rate = f / (h + f)
  239 5) CSI / critical success index = h / (h + f + m)
  240 6) OR / odds ratio = h*z / (f*m)
  241 7) BIAS / bias score = (h + f) / (h + m)
  242 8) TSS / true skill stat = POD - POFD
  243 9) HSS / Heidke skill score = 2 * (h*z - f*m) /
  244 ( (h+m)*(m+z) + (h+f)*(f+z) )
  245 10) ETS / equitable threat score = (h*z - f*m) /
  246 ( (f+m)*n + (h*z - f*m) )
  247 11) PRC / precision = h / (h + f)
  248 12) FSC / F-Score = 2 * (PRC*POD) / (PRC + POD)
  249 The TSS is also known as the Hanssen-Kuipers score,and is = h / (h+m) - f / (f+z).
  250 The F-Score can also be expressed as 2 * h / (1 + h + m).
  251 The input is always sanitized by taking the upper 2x2 part, using absolute values, and integer-ization.
  252 Warnings are issued if verbose == 1.
  253 */
  254 # (Skip the checks for efficiency if not verbose)
  255 if verbose
  256   if rows(in) != 2 || cols(in) != 2
  257     print &quot;Warning: Discarding input beyond 2x2&quot;
  258   endif
  259   if minr(minc(in)) &lt; 0
  260     print &quot;Warning: Converting negative entries&quot;
  261   endif
  262   if sum(abs( in - int(in) )) &gt; 1e-6
  263     print &quot;Warning: Rounding non-integers&quot;
  264   endif
  265 endif
  266 in = int(abs( in[1:2, 1:2] ))
  267 scalar h = in[1,1]
  268 scalar m = in[2,1]
  269 scalar f = in[1,2]
  270 scalar z = in[2,2]
  271 scalar n = sum(in)
  272 h_m = h + m
  273 f_z = f + z
  274 h_z = h + z
  275 h_f = h + f
  276 m_z = m + z
  277 f_m = f + m
  278 hz = h * z
  279 fm = f * m
  280 hzMfm = hz - fm
  281 POD = h / h_m
  282 POFD = f / f_z
  283 HR = h_z / n
  284 FAR = f / h_f
  285 CSI = h / (h_f + m)
  286 OR = hz / fm
  287 BIAS = h_f / h_m
  288 TSS = POD - POFD
  289 HSS = 2 * hzMfm / ( h_m * m_z + h_f * f_z )
  290 ETS = hzMfm / ( f_m * n + hzMfm )
  291 PRC = h / h_f
  292 FSC = 2 * h / (1 + h_m)
  293 matrix out = {POD; POFD; HR; FAR; CSI; OR; BIAS; TSS; HSS; ETS; PRC; FSC}
  294 rownames(out, &quot;POD POFD HR FAR CSI OR BIAS TSS HSS ETS PRC FSC&quot;)
  295 return out
  296 </code>
  297 </gretl-function>
  298 <gretl-function name="WSRcritical" type="matrix">
  299  <params count="3">
  300   <param name="n" type="int" min="4">
  301 <description>number of trials</description>
  302   </param>
  303   <param name="prob" type="scalar" min="0" max="1" default="0.05">
  304 <description>two-sided prob mass</description>
  305   </param>
  306   <param name="forcenorm" type="bool" default="0">
  307 <description>always use normal approx</description>
  308   </param>
  309  </params>
  310 <code>/* Tries to find the critical values (low/hi) where the
  311 two-sided area to the outside is as close as possible
  312 to the given prob.
  313 (Note that &quot;outside&quot; means including the critical values
  314 themselves in the exact/discrete case.)
  315 If we end up in the interior region not covered by the
  316 exact table (for prob far away from 0 and also from 1), we fall back to the normal approx.
  317 Returned is col vector {low; hi; epv}, where epv is the actual probability mass
  318 (close to prob but not equal in general for small samples).
  319 'low' and 'hi' can be non-integers in the normal
  320 approximation case.
  321 The normal approximation instead of the exact table values can be
  322 enforced with the 'forcenorm' argument.
  323 */
  324 approxtol = 0.05
  325 if n &lt;= 12 &amp;&amp; !forcenorm
  326   matrix lohiP = WPtable(n)
  327   lohiP[3, ] *= 2	# from one-sided to two-sided pv
  328   # find the min deviation from the p-value
  329   scalar ix = iminr(abs( lohiP[3, ] - prob ))
  330   # now check if we are too far in the interior region
  331   # where the exact table doesn't apply and the result
  332   # would be misleading;
  333   if abs(lohiP[3, ix] - prob ) &lt;= approxtol
  334     return lohiP[, ix]
  335   else
  336     forcenorm = 1 # fall back to normal approx.
  337   endif
  338 endif
  339 if n &gt; 12 || forcenorm # normal approx.
  340   Wmean = n * (n + 1) / 4
  341   Wsigma = sqrt(Wmean * (2 * n + 1) / 6)
  342   cv = critical(N, prob/2) # upper critical value
  343   hi = cv * Wsigma + Wmean
  344   lo = Wmean - (hi - Wmean) # symmetric lower
  345   return {lo; hi; prob}
  346 endif
  347 </code>
  348 </gretl-function>
  349 <gretl-function name="WSRpvalue" type="scalar">
  350  <params count="3">
  351   <param name="n" type="int" min="4">
  352 <description>number of trials</description>
  353   </param>
  354   <param name="W" type="scalar" min="0">
  355 <description>W test stat</description>
  356   </param>
  357   <param name="forcenorm" type="bool" default="0">
  358 <description>always use normal approx</description>
  359   </param>
  360  </params>
  361 <code>/* We return P(X &gt;= W), _not_ strict inequality!
  362 (in contrast to an earlier version!)
  363 In the interior region not covered by the exact table, the true value is
  364 &gt;= 12.5% according to WPtable, so typically based on such a p-value
  365 H0 would not be rejected.
  366 We fall back to the normal approximation in this region.
  367 In the extreme outer regions not explicitly covered by the table, the deviation
  368 from 0 or 1 will be smaller than 0.5% = 0.005. We return
  369 values 1e-3 or 1 - 1e-3 as an approximation here.
  370 The normal approximation instead of the exact table values can be
  371 enforced with the 'forcenorm' argument.
  372 Source of the table: Wilfrid J Dixon and Frank J. Massey, Jr., Introduction to Statistical Analysis, 2nd ed. (New York: McGraw-Hill, 1957), pp. 443-444.
  373 */
  374 # input check
  375 if W &gt; n * (n+1) / 2
  376   printf &quot;Warning: stat %d out of range for %d trials!\n&quot;, W, n
  377   return NA
  378 endif
  379 if W != int(W)
  380   # for a non-integer input (only possible for bindings?)
  381   # we always fall back to the normal approx
  382   forcenorm = 1
  383 else
  384   W = int(W)	# should be redundant, but who knows (numerically)
  385 endif
  386 approxdiff = 1e-3
  387 if n &lt;= 12 &amp;&amp; !forcenorm
  388   if W == 0	# trivial but possible case
  389     return 1
  390   endif
  391   matrix lohiP = WPtable(n)
  392   if W &gt; lohiP[1, cols(lohiP)]  &amp;&amp;  W &lt; lohiP[2, cols(lohiP)]
  393     # (&gt; max left, or &lt; min right) no exact values possible,
  394     # fall back to normal approx
  395     forcenorm = 1
  396   elif W &lt;= lohiP[1, 1]	# extreme lower tail
  397     return 1 - approxdiff
  398   elif W &gt; lohiP[2, 1]	# extreme upper tail
  399     return approxdiff
  400   elif W &lt;= lohiP[1, cols(lohiP)] # covered in lower tail
  401     # get the &lt;= W-1 prob
  402     scalar P = selifc( lohiP[3,], lohiP[1,] .= (W - 1) )
  403     # convert to &gt; (W-1) ==: &gt;= W prob
  404     return 1 - P
  405   elif W &gt;= lohiP[2, cols(lohiP)]	# covered in upper tail
  406     scalar P = selifc( lohiP[3,], lohiP[2,] .= W )
  407     return P
  408   endif
  409 endif # exact values
  410 if n &gt; 12 || forcenorm # normal approx.
  411   Wmean = n * (n + 1) / 4
  412   Wsigma = sqrt(Wmean * (2 * n + 1) / 6)
  413   Wstar = (W - Wmean) / Wsigma
  414   return pvalue(N, Wstar)
  415 endif
  416 </code>
  417 </gretl-function>
  418 <gretl-function name="powerset" type="strings">
  419  <params count="1">
  420   <param name="S" type="strings"/>
  421  </params>
  422 <code>/* Computes the powerset of the input S, i.e. all possible cominations
  423 of the string elements in S. (Including the empty set / empty
  424 string &quot;&quot;.) Each combination yields one string in the output array.
  425 Being a set, the ordering is not defined. */
  426 scalar l = nelem(S)
  427 scalar N = 2^l
  428 matrix P = zeros(N, l)
  429 matrix s = seq(0, N-1)'
  430 strings PS = array(N)
  431 loop i=1..l --quiet
  432   matrix a = s % 2
  433   P[,i] = a
  434   s = (s - a) / 2
  435 endloop
  436 loop i = 1..N --quiet
  437   string c = &quot;&quot;
  438   loop j = 1..l --quiet
  439     if P[i, j]
  440       c = c ~ &quot; &quot; ~ S[j]
  441     endif
  442   endloop
  443   PS[i] = c
  444 endloop
  445 return PS
  446 </code>
  447 </gretl-function>
  448 <gretl-function name="onemode" type="matrix">
  449  <params count="1">
  450   <param name="v" type="matrix">
  451 <description>Vector of values</description>
  452   </param>
  453  </params>
  454 <code>/* Returns one mode (modal value) of the input data v. If that is multi-modal,
  455 details of internal computer arithmetic can influence which of the modes
  456 is actually found. Returns a 2-element column vector with the modal value
  457 and its absolute frequency. Returns a 2-element column vector with the
  458 modal value and its absolute frequency. If v is an empty matrix (comprises
  459 only nan values) a 1x1 matrix with nan is returned. */
  460 matrix ret = mshape(NA, 1, 1)
  461 matrix v = vec(v)
  462 if sum(ok(v)) == 0
  463   printf &quot;\nError: No valid input data.\n&quot;
  464 else
  465   matrix E = ecdf(v)
  466   matrix howmuch = diff(0 | E[,2])[2:]    # make sure 1st is also diffed
  467   /* if multi-modal selection is not clear (finite precision digital
  468   arithmetic) */
  469   matrix where = imaxc(howmuch)
  470   matrix ret = E[where, 1] | howmuch[where]
  471   rnameset(ret, &quot;Mode frequency&quot;)
  472 endif
  473 return ret
  474 </code>
  475 </gretl-function>
  476 <gretl-function name="drill" type="matrix">
  477  <params count="3">
  478   <param name="x" type="matrices" const="true"/>
  479   <param name="rowspec" type="matrix" optional="true"/>
  480   <param name="colspec" type="matrix" optional="true"/>
  481  </params>
  482 <code># This function &quot;drills through&quot; a matrix array and returns a matrix;
  483 # for example, drill(x, 2, 3) returns a vector with the [2,3] elements
  484 # of all matrices in the x array. &quot;0&quot; means &quot;all&quot;. Returns an empty
  485 # 1x1 matrix in case of any error.
  486 # NOTA BENE: all matrices must be the same size
  487 # (might perhaps be relaxed in the future?)
  488 matrix ret = {}
  489 scalar n = nelem(x)
  490 scalar same_dim = 1
  491 if n == 0
  492   return ret
  493 endif
  494 ### check sizes
  495 nr = rows(x[1])
  496 nc = cols(x[1])
  497 loop i = 2..n --quiet
  498   same_dim = same_dim &amp;&amp; (rows(x[i]) == nr) &amp;&amp; (cols(x[i]) == nc)
  499   if !same_dim
  500     printf &quot;Matrix number %d with different dimensions\n&quot;, i
  501     funcerr &quot;Not all matrices are the same size.&quot;
  502   endif
  503 endloop
  504 ### process specs
  505 matrix rs = seq(1, nr)'	# default
  506 if exists(rowspec)
  507   rs = (rowspec[1] == 0) ? rs : vec(rowspec) # force to column
  508 endif
  509 matrix cs = seq(1, nc)' # default
  510 if exists(colspec)
  511   cs = (colspec[1] == 0) ? cs : vec(colspec) # force to column
  512 endif
  513 ### check for multiple or illegal specs
  514 scalar nrspec = rows(rs)
  515 scalar ncspec = rows(cs)
  516 if xmin(nrspec, ncspec) &gt; 1
  517   funcerr &quot;Cannot have multiple row and column specifications&quot;
  518 elif minc(rs|cs) &lt; 0
  519   funcerr &quot;Negative specification not allowed&quot;
  520 elif maxc(rs) &gt; nr
  521   printf &quot;(matrices have %d rows, but %d wanted)\n&quot;, nr, maxc(rs)
  522   funcerr &quot;Incorrect row specification&quot;
  523 elif maxc(cs) &gt; nc
  524   printf &quot;(matrices have %d columns, but %d wanted)\n&quot;, nc, maxc(cs)
  525   funcerr &quot;Incorrect column specification&quot;
  526 endif
  527 ### do the actual drilling
  528 if nrspec == 1
  529   ret = flatten(x)[rs,]
  530   ret = transp(mshape(ret, nc, n))
  531   ret = ret[,cs]
  532 elif ncspec == 1
  533   ret = flatten(x,1)[,cs]
  534   ret = mshape(ret, nr, n)
  535   ret = ret[rs,]
  536 endif
  537 return ret
  538 </code>
  539 </gretl-function>
  540 <gretl-function name="splitfname" type="strings">
  541  <params count="1">
  542   <param name="fn" type="string" const="true"/>
  543  </params>
  544 <code>/* The idea is to take a file name or full path and extract
  545 up to 3 components:
  546 1. The path prefix (may be empty; without the trailing / or \ )
  547 2. The &quot;base&quot; component of the file name, without the
  548 extension and without the path prefix
  549 3. The file extension (without the dot; may be empty)
  550 Example:
  551 Input string: &quot;/what/on/earth/isthisfile.gdt&quot;
  552 Output:
  553 defarray(&quot;/what/on/earth&quot;, &quot;isthisfile&quot;, &quot;gdt&quot;)
  554 (To separate the path prefix we just look for the last / _OR_
  555 \ (forward or backward slash) and also think of the possibility of //.)
  556 We use the \t character as internal separator because \n doesn't
  557 work properly with strsplit (bug fixed in gretl Dec 17th, 2019, for 2019d).
  558 */
  559 hasslash = instring(fn, &quot;/&quot;) || instring(fn, sprintf(&quot;\\&quot;))
  560 # Test for a dot in a position where it signals an extension
  561 # (not the Unixy other meanings)
  562 hasext = 0
  563 if regsub(fn, &quot;.*[\w\s]+\.[\w\s]+&quot;, &quot;!&quot;) == &quot;!&quot; &amp;&amp; (fn != &quot;!&quot;)
  564   hasext = 1
  565 endif
  566 if hasext &amp;&amp; hasslash    # the full monty
  567   string sepa = regsub(fn, &quot;(.*)[/\\]+([^/\\]*)\.([^\./\\]*)&quot;, &quot;\1\t\2\t\3&quot;)
  568   strings out = strsplit(sepa, &quot;\t&quot;)
  569   splitfname_check(out, 3)
  570 elif hasext           # only base file name and ext, no prefix
  571   string sepa = regsub(fn, &quot;(.*)\.([^\.]*)&quot;, &quot;\1\t\2&quot;)
  572   strings parts = strsplit(sepa, &quot;\t&quot;)
  573   splitfname_check(parts, 2)
  574   strings out = defarray(&quot;&quot;) + parts # empty prefix first
  575 elif hasslash        # no extension
  576   string sepa = regsub(fn, &quot;(.*)[/\\]+([^/\\]*)&quot;, &quot;\1\t\2&quot;)
  577   strings parts = strsplit(sepa, &quot;\t&quot;)
  578   splitfname_check(parts, 2)
  579   strings out = parts + defarray(&quot;&quot;) # empty ext last
  580 else                 # no slash, no ext, just simple filename
  581   strings out = defarray(&quot;&quot;, fn, &quot;&quot;)
  582 endif
  583 return out
  584 </code>
  585 </gretl-function>
  586 <gretl-function name="multi_instrings" type="matrix">
  587  <params count="2">
  588   <param name="lookinhere" type="strings"/>
  589   <param name="tofind" type="strings"/>
  590  </params>
  591 <code># Returns the positions (indices) in 'lookinhere' where any of the
  592 # strings from 'tofind' occur.
  593 # If there are duplicates in 'tofind' then the output may also
  594 # contain duplicate indices. Use uniq() or values() afterwards
  595 # if needed.
  596 if $version &lt; 20200
  597   funcerr &quot;Sorry, this function needs at least gretl 2020a.&quot;
  598   # (because of instrings)
  599 endif
  600 matrix which = {}
  601 # The case !nelem(tofind) used to be (in SVAR) some sort of default
  602 # and returned all indices instead of none
  603 # (don't remember why)
  604 # which = seq(1, nelem(lookinhere))
  605 if nelem(tofind)
  606   loop n = 1..nelem(tofind) -q
  607     which |= instrings(lookinhere, tofind[n])
  608   endloop
  609   # which = values(which) # this was done in the old names2indices
  610 endif
  611 return which # may still be empty
  612 </code>
  613 </gretl-function>
  614 <gretl-function name="correspondence" type="scalar">
  615  <params count="2">
  616   <param name="x" type="series"/>
  617   <param name="y" type="series"/>
  618  </params>
  619 <code># This function takes two series and establishes if there's a
  620 # 1-to-1 relationship between them, in which case it returns 2.
  621 # If there's a 1-to-n relationship, it returns 1. If there's
  622 # no relationship, it returns 0.
  623 # recode values so as to ensure we only get integer values
  624 matrix v = values(y)
  625 y = replace(y, v, seq(1,rows(v)))
  626 matrix v = values(x)
  627 x = replace(x, v, seq(1,rows(v)))
  628 matrix H = mxtab(x, y) .&gt; 0
  629 is_a_function = maxr(sumc(H)) == 1
  630 is_1_to_1 = maxc(sumr(H)) == 1
  631 if is_a_function
  632   ret = is_1_to_1 ? 2 : 1
  633 else
  634   ret = 0
  635 endif
  636 return ret
  637 </code>
  638 </gretl-function>
  639 <gretl-function name="fracorder" type="matrix">
  640  <params count="3">
  641   <param name="x" type="series" const="true"/>
  642   <param name="order" type="int" min="0" default="0"/>
  643   <param name="verbosity" type="bool" default="0"/>
  644  </params>
  645 <code># order == 0 means to leave gretl's default;
  646 # Applies all available estimators in gretl:
  647 # - Local Whittle
  648 # - Geweke Porter-Hudak
  649 # - Hurst (exponent minus 0.5)
  650 # (also note the generalized Hurst contributed function package)
  651 if $version &lt; 20202 # 2020c
  652   funcerr &quot;The fracorder function from 'extra' requires at least gretl 2020c.&quot;
  653 endif
  654 # translate the wanted order and whether to have printouts
  655 string sorder = !order    ? &quot;&quot; : sprintf(&quot;%d&quot;, order)
  656 string sverb  = verbosity ? &quot;&quot; : &quot;--quiet&quot;
  657 # Local Whittle and Geweke Porter-Hudak
  658 fractint x @sorder @sverb
  659 # point, se, teststat, pv
  660 matrix LocWhit    = vec($result)' ~ $test ~ $pvalue
  661 fractint x @sorder --gph @sverb
  662 matrix GewPorthud = vec($result)' ~ $test ~ $pvalue
  663 # Hurst
  664 matrix Hurst = mshape(NA, 1, 4)
  665 if $nobs &lt; 128 &amp;&amp; verbosity
  666   print &quot;Warning: not enough obs for Hurst, leaving result row as NA.&quot;
  667 elif $nobs &gt;= 128
  668   if verbosity
  669     hurst x --plot=none
  670   else
  671     outfile null	# workaround for the missing quiet option of hurst
  672       hurst x --plot=none
  673     end outfile 	# This is in &gt;=2018b, so OK here.
  674   endif
  675   # shifted to be comparable to d
  676   Hurst[1:2] = ($result[1] - 0.5) ~ $result[2]
  677 endif
  678 matrix ret = LocWhit | GewPorthud | Hurst
  679 cnameset(ret, &quot;estim SE teststat pval&quot;) # is in &gt;= 2018a
  680 rnameset(ret, &quot;LocalWh GPH Hurst&quot;)
  681 return ret
  682 </code>
  683 </gretl-function>
  684 <gretl-function name="splitfname_check" type="void" private="1">
  685  <params count="2">
  686   <param name="out" type="strings"/>
  687   <param name="num" type="int"/>
  688  </params>
  689 <code># Just a stupid debug checker and error catcher
  690 # for splitfname.
  691 if nelem(out) != num
  692   print out
  693   funcerr &quot;Shouldn't happen (weird input?)&quot;
  694 endif
  695 </code>
  696 </gretl-function>
  697 <gretl-function name="lin_int" type="series" private="1">
  698  <params count="1">
  699   <param name="y" type="series" const="true"/>
  700  </params>
  701 <code># originally from Jack's yahoo_get function
  702 # (helper for gap_filler)
  703 series DOK = diff(ok(y))
  704 series Beg = DOK == -1
  705 series End = DOK(+1) == 1
  706 series splen = 0
  707 splen = DOK == -1 ? 1 : DOK==1 ? 0 : (splen(-1) == 0 ? 0 : splen(-1) + 1)
  708 series y0 = NA
  709 series y0 = Beg ? y(-1) : y0(-1)
  710 series y1 = NA
  711 series y1 = End ? y(+1) : y1(-1)
  712 set skip_missing off
  713 matrix A = {y, y0, y1, splen}
  714 scalar t = lastobs(y)
  715 loop while t &gt; firstobs(y) --quiet
  716   if ok(End[t]) &amp;&amp; (End[t] == 1)
  717     scalar l = A[t, 4]
  718     dy = (A[t,3] - A[t,2]) / (l + 1)
  719     patch = A[t,2] + dy * seq(1,l)'
  720     A[t-l+1:t,1] = patch
  721     t -= l
  722   else
  723     t--
  724   endif
  725 endloop
  726 return A[,1]
  727 </code>
  728 </gretl-function>
  729 <gretl-function name="WPtable" type="matrix" private="1">
  730  <params count="1">
  731   <param name="n" type="int" min="4" default="12">
  732 <description>number of trials</description>
  733   </param>
  734  </params>
  735 <code>/*
  736 Upper and Lower Percentiles of the Wilcoxon Signed
  737 Rank Statistic W
  738 Returns a 3 x M matrix of critical values:
  739 (M differs for each n.)
  740 1st row: lower crit. vals, always ascending
  741 2nd row: upper crit. vals, always descending
  742 3rd row: corresponding P(X &lt;= lower crit.val)
  743 ( equal to P(X &gt;= upper crit.val) )
  744 */
  745 if n == 4
  746   return {0, 1} | {10, 9} | {0.062, 0.125}
  747 elif n == 5
  748   return seq(0, 3) | seq(15, 12) | {0.031,0.062,0.094,0.156}
  749 elif n == 6
  750   return seq(0, 5) | seq(21, 16) | {0.016,0.031,0.047,0.078,0.109,0.156}
  751 elif n == 7
  752   return seq(0, 7) | seq(28, 21) | {0.008,0.016,0.023,0.039,0.055,0.078,0.109,0.148}
  753 elif n == 8
  754   return seq(0, 9) | seq(36, 27) | {0.004,0.008,0.012,0.020,0.027,0.039,0.055,0.074,0.098,0.125}
  755 elif n == 9
  756   return seq(1, 12) | seq(44, 33) | {0.004,0.006,0.010,0.014,0.020,0.027,0.037,0.049,0.064,0.082,0.102,0.125}
  757 elif n == 10
  758   return seq(3, 16) | seq(52, 39) | {0.005,0.007,0.010,0.014,0.019,0.024,0.032,0.042, 0.053,0.065,0.080,0.097,0.116,0.138}
  759 elif n == 11
  760   return seq(5, 20) | seq(61, 46) | {0.005,0.007,0.009,0.012,0.016,0.021,0.027,0.034, 0.042,0.051,0.062,0.074,0.087,0.103,0.120,0.139}
  761 elif n == 12
  762   return seq(7, 24) | seq(71, 54) | {0.005,0.006,0.008,0.010,0.013,0.017,0.021,0.026, 0.032,0.039,0.046,0.055,0.065,0.076,0.088,0.102,0.117,0.133}
  763 endif
  764 </code>
  765 </gretl-function>
  766 <sample-script>
  767 # extra_sample.inp
  768 # Sample script for the extra.gfn package
  769 
  770 include extra.gfn --force
  771 set verbose off
  772 
  773 ######## Examples for functions not needing a dataset ##########
  774 
  775 print &quot;---- nearPSD -----&quot;
  776 matrix mc = {1, -1.2; -1.2, 1}	# eigenvalues -0.2 and 2.2
  777 eval nearPSD(&amp;mc)
  778 print mc
  779 
  780 print &quot;---- zeroifclose -----&quot;
  781 mc = {1, 1e-12}
  782 eval zeroifclose(&amp;mc, 1e-10)
  783 print mc
  784 
  785 print &quot;---- truncnorm -----&quot;
  786 m = 0.5
  787 s = 1
  788 matrix left = truncnorm(20, m, s, NA, 1)
  789 matrix right = truncnorm(20, m, s, -1, NA)
  790 matrix both = truncnorm(20, m, s, -1, 1)
  791 print left right both
  792 
  793 print &quot;---- commute -----&quot; ##
  794 # vec application
  795 eval commute( vec( I(2) ~ ones(2,1) ), 2, 3)
  796 # commutation of the Kronecker product
  797 matrix A = mnormal(4,1)
  798 matrix B = mnormal(1,5)
  799 eval B**A
  800 matrix temp = commute(A**B, rows(B), rows(A))
  801 eval commute(temp, cols(B), cols(A), 1) # post-multipl.; should give B**A
  802 
  803 print &quot;---- eliminate -----&quot;
  804 eval eliminate(vec(I(3)))
  805 
  806 print &quot;---- duplicate -----&quot;
  807 eval duplicate(vech(I(3)))
  808 
  809 print &quot;---- scores2x2 -----&quot;##
  810 eval scores2x2({24, 4; 12, 9})
  811 
  812 print &quot;---- WSRpvalue -----&quot;
  813 print &quot;=============================&quot;
  814 print &quot;Some examples for WSRpvalue()&quot;
  815 loop N=11..14 -q            # loop over n
  816   loop i = 6..7 -q        # loop over W
  817     printf &quot;n=$N, W=$i\n&quot;
  818     printf &quot;try exact: %g, approx: %g\n&quot;, WSRpvalue(N, i), WSRpvalue(N, i, 1)
  819   endloop
  820 
  821   Wmax = N * (N + 1) / 2
  822   loop i = (Wmax - 7)..(Wmax - 6) -q
  823     printf &quot;n=$N, W=$i\n&quot;
  824     printf &quot;try exact: %g, approx: %g\n&quot;, WSRpvalue(N, i), WSRpvalue(N, i, 1)
  825   endloop
  826 
  827 endloop
  828 
  829 print &quot;---- WSRcritical -----&quot;
  830 print &quot;----------------------------------------&quot;
  831 print &quot;Some examples for WSRcritical()&quot;
  832 loop N=11..13 -q            # loop over n
  833   Wmax = N * (N + 1) / 2
  834   loop i=3..4 -q        # loop over W
  835     printf &quot;n = $N, Pr = %g\n&quot;, i/100
  836     printf &quot;try exact:&quot;
  837     eval WSRcritical(N, i/100)'
  838     printf &quot;approx:&quot;
  839     eval WSRcritical(N, i/100, 1)'
  840   endloop
  841 endloop
  842 
  843 print &quot;---- powerset -----&quot;
  844 strings pp =  powerset(defarray(&quot;a&quot;, &quot;b&quot;, &quot;c&quot;, &quot;x&quot;, &quot;y&quot;, &quot;z&quot;))
  845 loop i=1..nelem(pp) -q
  846   eval pp[i]
  847 endloop
  848 
  849 print &quot;---- onemode -----&quot;
  850 eval onemode({1,2,3,1,1,2})
  851 
  852 print &quot;---- drill -----&quot;
  853 matrix A = mshape(seq(1,9),3,3)
  854 matrix B = A
  855 matrices x = defarray(A, B)
  856 eval drill(x, 2)			# one row
  857 eval drill(x, , 1)			# one column
  858 eval drill(x, 2, 1)			# one row &amp; one column
  859 eval drill(x, {1,3}, 2)		# multiple rows
  860 eval drill(x, 2, {1,3})		# multiple columns
  861 
  862 print &quot;---- splitfname -----&quot;  ##
  863 eval splitfname(sprintf(&quot;\\do\\.you\\believe.input&quot;))
  864 
  865 print &quot;---- multi_instrings -----&quot;  ##
  866 strings AA = defarray(&quot;an&quot;,&quot;apple&quot;,&quot;and&quot;,&quot;an&quot;,&quot;egg&quot;)
  867 strings BB = defarray(&quot;an&quot;,&quot;dog&quot;,&quot;egg&quot;)
  868 eval multi_instrings(AA, BB)
  869 
  870 ######## Examples that require a dataset ############
  871 
  872 open denmark	# get some data
  873 
  874 print &quot;---- gap_filler -----&quot;  #
  875 series play = LRM
  876 play[10] = NA	# insert missing
  877 series playout = gap_filler(play, 2)
  878 print playout
  879 
  880 print &quot;---- winsor -----&quot;  #
  881 series winDM = winsor(ldiff(LRM))
  882 series winDM_asy = winsor(ldiff(LRM), 0.01, 0.8)
  883 print winDM_asy
  884 
  885 print &quot;---- correspondence -----&quot;  ##
  886 
  887 x = $obsminor == 4
  888 C = seasonals(1,1)
  889 
  890 printf &quot;correspondence(x, S4) = %d\n\n&quot;, correspondence(x, S4)
  891 
  892 print &quot;---- fracorder ------&quot;
  893 open djclose.gdt --quiet	# need longer data here
  894 matrix M = fracorder(ldiff(djclose), 20)
  895 print M
  896 </sample-script>
  897 </gretl-function-package>
  898 </gretl-functions>