"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020b/addons/SVAR/SVAR_signrestr.inp" (11 Apr 2020, 21030 Bytes) of package /linux/misc/gretl-2020b.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Charmm source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "SVAR_signrestr.inp": 2020a_vs_2020b.

    1 # For now drill(), strpos_allin() and drawnormwis() come from:
    2 # include SVAR_compatibility_funcs.inp
    3 #
    4 # And the following functions come from include SVAR_signrestr_utils.inp:
    5 #  correspondence, normalize, check_one_irf, check_irfs, check_id;
    6 #  safetycheck1, get_n_exotic, rot_redraw, prepAis, get_id_mat
    7 
    8 ###
    9 ### functions for setting the set constraints
   10 ###
   11 
   12 function void SVAR_SRfull(bundle *mod,
   13                           string yname, string sname,
   14                           scalar lo[NA], scalar hi[NA],
   15                           int ini[0::0], int fin[0::0])
   16 
   17     # this function will add to the appropriate element of
   18     # the bundle "mod" a set of sign restrictions relative
   19     # to one variable and one shock, organised as a row
   20     # vector with the following elements; input values are:
   21     #
   22     # mod      : the VAR model bundle (in pointer form)
   23     # yname    : a string with the name of the observable to examine
   24     # sname    : a string with the shock name
   25     # lo, hi   : bounds (NA for +- infty)
   26     # ini, fin : IRF interval to check (0-based)
   27     #
   28     # for example, if we assume that a monetary policy shock will
   29     # impact prices with a negative sign on lags from 0 to 4, then
   30     # the function call should be
   31     #
   32     # SVAR_SRfull(mod, "p", 1, NA, 0, 0, 4)
   33     #
   34     # assuming the series for prices is called "p" and that we want
   35     # our monetary policy shock to come 1st in the ordering (useful for
   36     # models with sign restrictions on more than one shock)
   37 
   38     # (Switched the yname/sname ordering to match SVAR_SRplain; Sven)
   39 
   40     if mod.type != 10
   41         funcerr "Wrong model type for set identification restrictions"
   42     endif
   43     if !inbundle(mod, "SRest")
   44         matrix mod.SRest = {}
   45     endif
   46 
   47     matrix snum = strpos_allin(mod.snames, sname)
   48     errmsgshockmatch(snum, sname)
   49 
   50     matrix k = strpos_allin(mod.Ynames, yname)
   51     if nelem(k) != 1
   52         printf "Variable %s is not in the model\n", yname
   53         print " (or several ones with that name)"
   54         funcerr "Specify correct and unique variable names"
   55     endif
   56 
   57     if missing(lo) && missing(hi)
   58         funcerr "No bounds specified!\n"
   59     endif
   60 
   61     lo = ok(lo) ? lo : -$huge
   62     hi = ok(hi) ? hi :  $huge
   63 
   64     if lo > hi
   65         printf "Invalid bound specified! [%g : %g]\n", lo, hi
   66         funcerr "Provide valid interval"
   67     endif
   68 
   69     if mod.horizon <= fin
   70         mod.horizon = fin + 1 # to be on the safe side
   71     endif
   72 
   73     # Attention: Here the ini and fin inputs are transformed
   74     # to 1-based matrix indexing form, not keeping the actual
   75     # horizon meaning
   76     mod.SRest |= {k, lo, hi, ini+1, fin+1, snum}
   77 
   78 end function
   79 
   80 function void SVAR_SRplain(bundle *mod,
   81                            string yname,
   82                            string sname,
   83                            string what,
   84                            int length[0::0], int ini[0::0])
   85     # shorthand for sign restrictions in narrow sense
   86     # (Model type is checked downstream in SVAR_SRfull)
   87     scalar hi, lo
   88     if what == "+"
   89         hi = NA
   90         lo = 0
   91     elif what == "-"
   92         hi = 0
   93         lo = NA
   94     else
   95         funcerr "invalid restriction"
   96     endif
   97 
   98     SVAR_SRfull(&mod, yname, sname, lo, hi, ini, ini+length)
   99 end function
  100 
  101 #####################################################################
  102 
  103 function void SVAR_SRexotic(bundle *mod, string chkstr,
  104                             strings shocks,
  105                             int length[0::0], int ini[0::0],
  106                             bool needs_model[0])
  107 
  108     /* Allows to specify more "exotic" restrictions than just an interval
  109        for one IRF over some horizon. For example, a difference or ratio
  110        of two IRFs. But the format is necessarily more free-floating here,
  111        see the documentation.
  112 
  113        chkstr: a string with a valid numerical evaluation of some function
  114        of the IRFs at a given horizon tau. The IRF_tau matrix must be
  115        hardcoded as "M". Cross-horizon restrictions are impossible.
  116 
  117        length: length of the horizon span over which the restriction should
  118        hold
  119        ini: first horizon after shock impact to consider
  120        (the span is then from ini to ini+length, where impact timing is 0;
  121        so the default settings just concern the impact effect)
  122 
  123        shocks: strings array with the names of the shocks that enter this
  124        restriction (it's the user's responsibility to get this right;
  125        in the futur maybe we can verify or check this)
  126 
  127        needs_model: If yes, this becomes a "super-exotic" restriction,
  128        because it needs more input than just the boundaries; for example
  129        an estimated coefficient (or just the observed variables).
  130 
  131     */
  132 
  133     if mod.type != 10
  134         funcerr "Wrong model type for set identification restrictions"
  135     elif $version < 20193
  136         # (because of feval and funcerr() etc., and maybe nested arrays)
  137         funcerr "Sorry, need at least gretl 2019d for exotic restrictions"
  138     endif
  139 
  140     if !inbundle(mod, "exoticSR")
  141         strings checks = null
  142         bundle mod.exoticSR = defbundle("checks", checks, "spans", {},
  143           "super", {}, "eshocks", {})
  144     endif
  145 
  146     # store the relevant shock numbers
  147     loop s = 1..nelem(shocks) -q
  148         matrix snum = strpos_allin(mod.snames, shocks[s]) # can be several!
  149         if !nelem(snum)
  150             funcerr("No shock named %s found\n", shocks[s])
  151         endif
  152         mod.exoticSR.eshocks |= snum
  153     endloop
  154 
  155     finp1 = ini + length + 1 # the indexing convention here matches SRfull, +1
  156     if mod.horizon <= finp1
  157         mod.horizon = finp1 + 1 # on the safe side
  158     endif
  159 
  160     mod.exoticSR.checks += chkstr
  161     mod.exoticSR.spans = mod.exoticSR.spans | { ini + 1, finp1 }
  162     mod.exoticSR.super = mod.exoticSR.super | { needs_model }
  163 
  164 end function
  165 
  166 
  167 #-----------------------------------------------------------------
  168 
  169 function matrices gen_SRirfs(bundle b, int numh, matrix to_cum)
  170 
  171     # This calculates results for 0..mod.horizons, which is actually
  172     # like in the earlier SVAR code.
  173     # TODO: Maybe merge this function with earlier ones from SVAR?
  174     #
  175     # So numh == mod.horizon + 1.
  176     # (Memo: the doc has been (up to 2020) slightly wrong in this
  177     #  regard; it claimed h=mod.horizons rows in the IRFs matrix,
  178     #  one less than true)
  179 
  180     matrix B = b.B
  181     k = b.exoterms    # just a number?
  182 
  183     # On the outside, this cumulation spec is taken from the mod.
  184     # (Empty matrix to_cum if doesn't apply.)
  185     any_cumul = nelem(to_cum) > 0
  186 
  187     n = cols(B)
  188     p = (rows(B) - k) / n # VAR lags
  189     matrix comp = B[k+1:,]'
  190 
  191     if p > 1
  192         comp = comp | I(n*(p-1), n*p)
  193     endif
  194 
  195     matrix K = cholesky(b.Sigma) * b.rot
  196     matrices ret = array(numh)
  197     ret[1] = K  # impact effect
  198 
  199     matrix tmp = I(rows(comp))
  200     loop i = 2 .. numh --quiet
  201         tmp = tmp * comp
  202         matrix irf_i = tmp[1:n,1:n] * K
  203         if any_cumul
  204             irf_i[to_cum,] = irf_i[to_cum,] + ret[i-1][to_cum,]
  205         endif
  206         ret[i] = irf_i
  207     endloop
  208 
  209     return ret
  210 end function
  211 
  212 # -------------------
  213 
  214 # helper function for SVAR_SRdraw
  215 function scalar exotic_inner_check(scalar *chk, int i, const bundle mod,
  216                                    const matrices irfs)
  217     # (chk is a boolean, but gretl doesn't allow pointerized int or bool)
  218     # TODO: Check whether this setup is really compatible with a situation
  219     # where we have only partial identification from the "standard"
  220     # (non-exotic) restrictions, but the exotic restrictions then refer
  221     # to additional shocks not covered before. (Somehow I'm worrying 
  222     # that the previous multiplication with the non-square id_matrix in the
  223     # partial id case removes the needed irfs that would still have to be
  224     # checked against the exotic restrictions. - Sven)
  225 
  226     string exo_restrict = mod.exoticSR.checks[i]
  227     scalar inip1 = mod.exoticSR.spans[i,1]
  228     scalar finp1 = mod.exoticSR.spans[i,2]
  229     scalar needs_model = mod.exoticSR.super[i]
  230 
  231     out = 0
  232     loop t = inip1..finp1 --quiet
  233         matrix M = irfs[t] # hardcoded "M" for the API!
  234         if needs_model
  235             # FIXME This part looks very broken...!?
  236             # or very outdated...
  237             scalar exocheck = feval(exo_restrict, M, mod)
  238         else
  239             scalar exocheck = @exo_restrict
  240         endif
  241 
  242         if !exocheck
  243             chk = 0
  244             out = 1
  245             if 0
  246                 printf "Exotic check %s failed on M[%d]:\n", exo_restrict, t-1
  247                 printf "%7.4f\n", M
  248             endif
  249             break
  250         endif
  251     endloop
  252     return out
  253 end function
  254 
  255 function bundles SVAR_SRdraw(bundle *mod, int rep, bool DO_BAYES[0],
  256                              scalar coveralpha[0:1:0.9], int maxiter[10000])
  257 
  258     # This function draws a random rotations of the Cholesky model
  259     # until "rep" draws satisfying the sign restrictions have come up.
  260     # If the DO_BAYES flag is on, the VAR parameters (including Sigma)
  261     # are resampled too; otherwise, they are kept fixed at the ols
  262     # estimates.
  263 
  264     if !inbundle(mod, "SRest") && !inbundle(mod, "exoticSR")
  265         funcerr "No set id restrictions found"
  266     endif
  267 
  268     #
  269     # Recover VAR parameters from the model bundle
  270     #
  271 
  272     matrix B = mod.mu | mod.VARpar'
  273     scalar numh = mod.horizon + 1
  274 
  275     DO_NORM = 1
  276 
  277     matrix iXX = {}
  278     if DO_BAYES     # this is only needed for sampling from the posterior
  279         matrix Data = (mod.X ~ lags(mod.p, mod.Y, 1))[mod.p + 1 : ,]
  280         iXX = invpd(Data'Data)
  281         matrix mod.mreg = Data # needed to re-calc the residuals E later
  282     endif
  283 
  284     if 0
  285         # print a few parameters so we don't screw up with lags etc.
  286         safetycheck1(mod.Sigma, iXX, B)
  287     endif
  288 
  289     #
  290     # prepare the output array and other auxiliary stuff
  291     #
  292 
  293     bundles drawn = array(0)
  294     good = 0
  295     iter = 1
  296 
  297     # which shocks are to be identified via sign restrictions?
  298     matrix shocks = values(mod.SRest[,6])
  299     matrix allshocks = shocks
  300 
  301     # do we have exotic restrictions?
  302     n_exotic = get_n_exotic(mod)
  303     # add the ones from exotic restrictions
  304     matrix allshocks = n_exotic ? values(shocks | mod.exoticSR.eshocks) : shocks
  305 
  306     # Now extract the corresponding shock names
  307     strings mod.SRid_snames = getstrings_byindex(mod.snames, allshocks)
  308 
  309     matrices Ais = prepAis(mod, shocks)
  310 
  311     # responses to be cumulated?
  312     matrix to_cum = mod.ncumul ? mod.cumul : {}
  313 
  314     #
  315     # the main loop
  316     #
  317 
  318     bundle b = defbundle("B",B, "Sigma",mod.Sigma, "exoterms",rows(mod.mu), \
  319       "df", mod.T - rows(B), \
  320       "T", mod.T) # df was: mod.T - nelem(B)
  321 
  322     loop while good<rep && iter<=maxiter --quiet
  323 
  324         # first just using non-exotic restrictions
  325         rot_redraw(&b, mod, DO_BAYES, B, iXX) # b.rot, b.B, b.Sigma
  326         matrices irfs = gen_SRirfs(b, numh, to_cum)
  327         matrix id_matrix = get_id_mat(mod, shocks, irfs, Ais, DO_NORM)
  328 
  329         chk = DO_NORM ? (rows(id_matrix) == rows(shocks)) : 0 # this check 1st!
  330         chk = chk && check_id(id_matrix)
  331 
  332         if chk
  333             # perform the "reflection-permutation" step
  334             loop i = 1..nelem(irfs) --quiet
  335                 irfs[i] = irfs[i] * id_matrix'
  336             endloop
  337 
  338             # --- now check for exotic restrictions if necessary -------
  339             if n_exotic
  340                 loop i = 1..n_exotic --quiet
  341 
  342                     out = exotic_inner_check(&chk, i, mod, irfs)
  343                     if out
  344                         break
  345                     endif
  346                 endloop
  347             endif
  348 
  349             if chk
  350                 good++
  351                 b.irfs = irfs
  352 
  353                 # copy the results of accepted draw to output bundles array
  354                 drawn += b
  355             endif
  356         endif
  357 
  358         factor = (DO_BAYES ? 10 : 50) * xmax(1,floor(sqrt(maxiter/1000)))
  359 
  360         if iter % factor == 0
  361             printf "draw %5d done (good so far = %5d)\n", iter, good
  362             flush
  363         endif
  364 
  365         iter++
  366     endloop
  367 
  368     if good < rep
  369         printf "\n\nWARNING: COULDN'T ACHIEVE THE DESIRED NUMBER OF DRAWS"
  370         printf " (only %d out of %d)\n", good, rep
  371     else
  372         printf "draw %5d done (good = %5d)\n", iter-1, good
  373     endif
  374     printf " (acceptance ratio: %g)\n", good/(iter-1)
  375 
  376     # copy info into model
  377     mod.SRiter = iter - 1
  378     mod.SRacc  = good
  379     mod.DO_BAYES = DO_BAYES
  380     mod.SRcoveralpha = coveralpha
  381 
  382     # integrate some summary results into the model bundle
  383     if good
  384         bundle IrfData = IRF_plotdata(drawn, coveralpha)    
  385         matrix mod.SRirfmeans = IrfData.irfSRmeans
  386         matrix mod.SRirfserrs = IrfData.irfSRserrs
  387         matrix mod.SRlo_cb    = IrfData.lo_cb
  388         matrix mod.SRhi_cb    = IrfData.hi_cb
  389         matrix mod.SRirfmeds  = IrfData.irfSRmeds # medians
  390 
  391         # also store everything in the main bundle if wanted
  392         if inbundle(mod, storeSRirfs) && mod.storeSRirfs
  393             # (idiom possible because for SR we require recent gretl)
  394             mod.acc_draws = drawn
  395         endif
  396 
  397     else
  398         print "Warning: No accepted draws, plotting impossible"
  399     endif
  400 
  401     return drawn    # may be 0-element
  402 end function
  403 
  404 ###
  405 ### plotting functions
  406 ###
  407 
  408 function void SVAR_spagplot(const bundle mod,
  409                             string vname, string sname)
  410 
  411     if !inbundle(mod, "acc_draws")
  412         funcerr "Need collection acc_draws as model bundle member"
  413     endif
  414 
  415     # check whether it's one of the SR-id'ed shocks
  416     s_inSR = nelem(strpos_allin(mod.SRid_snames, sname)) ? 1 : 0
  417     # get shock number
  418     matrix snum = strpos_allin(mod.snames, sname)
  419     errmsgshockmatch(snum, sname)
  420     # get var number
  421     matrix vnum = strpos_allin(mod.Ynames, vname)
  422     if !nelem(vnum)
  423         printf "Regarding variable name %s:\n", vname
  424         funcerr "Variable unknown"
  425     endif
  426 
  427     rep = nelem(mod.acc_draws)  # was bs
  428     numh = mod.horizon + 1
  429 
  430     matrix IRFs = zeros(numh, rep) ~ seq(0, mod.horizon)'
  431 
  432     loop i = 1..rep --quiet
  433         IRFs[,i] = drill(mod.acc_draws[i].irfs, vnum[1], snum[1])
  434     endloop
  435 
  436     plot IRFs
  437         option with-lines
  438         #    literal set linetype cycle 4
  439         literal set linetype 1 lc rgb "#000000" lw 0.1
  440         literal set linetype 2 lc rgb "#000000" lw 0.1
  441         literal set linetype 3 lc rgb "#000000" lw 0.1
  442         literal set linetype 4 lc rgb "#000000" lw 0.1
  443         literal set linetype 5 lc rgb "#000000" lw 0.1
  444         literal set linetype 6 lc rgb "#000000" lw 0.1
  445         literal set linetype 7 lc rgb "#000000" lw 0.1
  446         literal set linetype 8 lc rgb "#000000" lw 0.1
  447         literal set xlabel 'lags'
  448         literal set nokey
  449         printf "set title '%s -> %s'\n", sname, vname
  450     end plot --output=display
  451 end function
  452 
  453 
  454 function bundle IRF_plotdata(bundles bs, scalar coveralpha[0.9])
  455 
  456     # here we store the data for plotting the irfs in a matrix, conceptually
  457     # similar to the "bands" matrix we use in SVAR_boot, that is with H rows
  458     # (=accepted draws, so that descriptive statistics are easy to compute)
  459     # and ns * nv * horizon columns
  460     # (But the outputs will be reshaped after the H is aggregated out;
  461     #  with the horizon in rows and numvariables * numshocks columns.)
  462 
  463     # get dimensions
  464     H = nelem(bs)    # number of available (accepted) draws
  465     nvXns = rows(bs[1].Sigma) * cols(bs[1].irfs[1])
  466     h = nelem(bs[1].irfs)
  467 
  468     matrix A = zeros(H, nvXns * h)
  469 
  470     loop i = 1..H --quiet
  471         matrices irfs = bs[i].irfs
  472         matrix sel = seq(1, nvXns * h, h)
  473 
  474         loop t = 1..h --quiet
  475             # compute the relevant IRFs
  476             # and store them in the appropriate column
  477             A[i, sel] = vec(irfs[t])'     # was: b.irfs
  478             sel += 1
  479         endloop
  480     endloop
  481 
  482     matrix m = meanc(A) # point-wise means
  483     matrix s = sdc(A)
  484 
  485     # output bundle (H: num of [accepted] replications)
  486     # (leave out ret.raw = A)
  487     bundle ret = defbundle("rep",H, "coveralpha",coveralpha, "biascorr",0)
  488     matrix ret.irfSRmeans  = mshape(m, h, nvXns)
  489     matrix ret.irfSRserrs  = mshape(s, h, nvXns)
  490 
  491     # point-wise pseudo confidence bands (ripped from SVAR_boot)
  492     q_alpha = 0.5 * (1 - coveralpha)
  493     matrix locb = quantile(A, q_alpha)
  494     matrix hicb = quantile(A, 1 - q_alpha)
  495     matrix mdn  = quantile(A, 0.5)
  496 
  497     matrix ret.lo_cb     = mshape(locb, h, nvXns)
  498     matrix ret.hi_cb     = mshape(hicb, h, nvXns)
  499     matrix ret.irfSRmeds = mshape(mdn, h, nvXns)
  500 
  501     return ret
  502 end function
  503 
  504 function void SVAR_SRirf(const bundle mod "the model bundle",
  505                          const strings whichvars[null] "which variables",
  506                          const strings whichshocks[null] "which shocks",
  507                          bool meanormed[0] "choose 1 for medians")
  508 
  509     # TODO: make the function return the plot code as string
  510     # TODO: check if it makes sense to use the existing IRFgrph and
  511     #       IRFsave functions internally
  512     #       (perhaps the code here is actually better (more modern))
  513 
  514     # Omitting the strings array inputs means: do all.
  515 
  516     if !inbundle(mod, "SRid_snames")
  517         funcerr "No SR (sign restriction) specs found"
  518     elif !inbundle(mod, "SRirfmeans")
  519         funcerr "No SR IRF results found - did you run SVAR_SRdraw?"
  520     endif
  521 
  522     # bundle IrfData = IRF_plotdata(SR_IRFs, coveralpha)
  523 
  524     # Which variables:
  525     if !exists(whichvars)
  526         strings whichvars = array(0) # default induces "all"
  527     endif
  528     matrix wvarix = names2indices(mod.Ynames, whichvars)
  529     if !nelem(wvarix)
  530         funcerr "No matching variable names found"
  531     endif
  532 
  533     # Which shocks:
  534     if !exists(whichshocks)
  535         strings whichshocks = array(0)
  536     endif
  537     matrix wshoix = names2indices(mod.SRid_snames, whichshocks)
  538     if !nelem(wshoix)
  539         funcerr "No matching shock names found"
  540     endif
  541 
  542     ## The actual plotting
  543     loop metaj = 1..nelem(wshoix) -q # for the shocks
  544         j = wshoix[metaj]
  545 
  546         loop metai = 1..nelem(wvarix) -q # for the target variables
  547             i = wvarix[metai]
  548 
  549             # pick the right position in the matrix
  550             k = (j-1) * nelem(mod.Ynames) + i
  551 
  552             matrix center = meanormed ? mod.SRirfmeds : mod.SRirfmeans
  553             matrix grph = center[,k] ~ seq(0, mod.horizon)' # IrfData.irfSRmeans
  554             matrix band = mod.SRlo_cb[,k] ~ mod.SRhi_cb[,k] # IrfData.lo_cb, IrfData.hi_cb
  555             band = 0.5 * band * {1, 1; 1, -1}
  556 
  557             plot grph
  558                 options with-lines fit=none
  559                 options band=band,1 band-style=fill
  560                 printf "set title '%s shock -> %s'", mod.SRid_snames[j], \
  561                   mod.Ynames[i]
  562             end plot --output=display
  563         endloop
  564     endloop
  565 
  566     # return the gnuplot plotting code here eventually
  567 end function
  568 
  569 ##############################
  570 
  571 function scalar SRgetbest(bundle *mod,
  572                         string vname, string sname,
  573                         int ini[0::0], int length[0::0],
  574                         int disttarg[0:1:0] "measure to target" \
  575                         {"mean", "median"},
  576                         string loss[null])
  577     # Function to find the set-id accepted draw which best fits a
  578     # a certain IRF center
  579     # vname: e.g. "Y"
  580     # sname: e.g. "supply" (naming must exist in the bundle)
  581     # ini: horizon range start (0 = impact)
  582     # length: horizon end is ini + length
  583     # disttarg: which part of the per-horizon distribution should be
  584     # targeted;
  585     # loss: how to evaluate deviations from the target;
  586     #      possible values: "abs", "quad" (default), ...
  587 
  588     ## implement defaults and error checks:
  589     if !exists(loss)
  590         string loss = "quad"
  591     elif tolower(loss) != "abs" && tolower(loss) != "quad"
  592         print "Warning, unrecognized loss, using 'quad'"
  593         loss = "quad"
  594     else
  595         loss = tolower(loss)
  596     endif
  597 
  598     if ini + length > mod.horizon
  599         funcerr "ini + length too large (increase horizon)"
  600     endif
  601     if !inbundle(mod, "SRacc")
  602         funcerr "No info on accepted draws found"
  603     endif
  604 
  605     scalar v_ix = names2indices(mod.Ynames, defarray(vname))
  606     scalar s_ix = names2indices(mod.snames, defarray(sname))
  607     startix = ini + 1   # offset for 1-based
  608     finix = startix + length
  609     matrix center = !disttarg ? mod.SRirfmeans : mod.SRirfmeds
  610     # pick the right position in the matrix:
  611     k = (s_ix - 1) * nelem(mod.Ynames) + v_ix
  612     center = center[,k]
  613 
  614     best = $huge
  615     bestdraw = 0
  616     loop i = 1..mod.SRacc -q
  617         matrix IRFvec = drill(mod.acc_draws[i].irfs, v_ix, s_ix)
  618         matrix distance = (IRFvec - center)[startix : finix]
  619         if loss == "abs"
  620             scalar objective = sum(abs(distance))
  621         elif loss == "quad"
  622             scalar objective = sum(distance.^2)
  623         endif
  624         if objective < best
  625             best = objective
  626             bestdraw = i
  627         endif
  628     endloop
  629 
  630     if !bestdraw
  631         funcerr "shouldn't happen"
  632     endif
  633 
  634     # store for further use
  635     mod.bestdraw = bestdraw
  636     return bestdraw
  637 end function