"Fossies" - the Fresh Open Source Software Archive

Member "gretl/functions/SVAR/SVAR.gfn" (21 Nov 2020, 154574 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="SVAR" needs-time-series-data="true" minver="2017a" lives-in-subdir="true">
    4 <author email="r.lucchetti@univpm.it">Riccardo &quot;Jack&quot; Lucchetti and Sven Schreiber</author>
    5 <version>1.92</version>
    6 <date>2020-09-17</date>
    7 <description>Structural VARs</description>
    8 <tags>C32</tags>
    9 <label>Structural VARs</label>
   10 <help>
   11 pdfdoc:SVAR.pdf
   12 </help>
   13 <data-files count="1">
   14 examples </data-files>
   15 <gretl-function name="SVAR_setup" type="bundle">
   16  <params count="5">
   17   <param name="type_string" type="string"/>
   18   <param name="lY" type="list"/>
   19   <param name="lX" type="list" optional="true"/>
   20   <param name="varorder" type="int" min="1" default="1"/>
   21   <param name="checkid_storeSRirfs" type="bool" default="1"/>
   22  </params>
   23 <code>/*
   24 This creates a bundle holding the info on the SVAR; this will be
   25 filled in 3 steps: This function inserts the initial info: sample
   26 size, data and so on. Then, more stuff will have to be added later:
   27 the VAR estimates in packed form (see below), then the SVAR estimates.
   28 The scalar &quot;step&quot; keeps track of the stage we're at.
   29 (maybe not needed any more)
   30 The switch 'checkid_storeSRirf' is now &quot;overloaded&quot; with different
   31 meanings, depending on the model type:
   32 - For traditional types it takes the meaning of 'checkident'.
   33 (Default 1 = yes, do it.)
   34 - For set-identified (sign restriction) models checking the
   35 identification separately makes no sense, so here this becomes
   36 'storeSRirf', serving to decide whether the full IRF data of all
   37 accepted draws should be stored in the model bundle. This can
   38 inflate the bundle quite a bit, but in general shouldn't be a
   39 problem on modern computers.
   40 (Default 1 = yes, store all IRF data.)
   41 In both cases it is just copied into the model bundle and still needs
   42 to be processed later (e.g. within SVAR_estimate).
   43 */
   44 # Drop collinear exogenous terms
   45 lX = dropcoll(lX)
   46 scalar type = modeltype(type_string)
   47 if type == 4	# SVEC
   48   # SVEC treats restricted deterministic terms specially as per Johansen,
   49   # so we'll drop them if present,
   50   # any other exogenous regressors are under the user's responsibility
   51   # (eg centred dummies)
   52   # the following currently doesn't work in gretl w.r.t. the time trend,
   53   # but fails silently (but const works -- still current in Jan 2019)
   54   list lX -= const time
   55   # therefore check for remaining collinearity
   56   genr time
   57   list lcheck = lX const time
   58   matrix mcheck = { lcheck }
   59   if rank(mcheck'mcheck) &lt; nelem(lcheck)
   60     funcerr &quot;Remove constant and trend terms from X list!&quot;
   61   endif
   62 endif
   63 scalar n = nelem(lY)
   64 scalar k = nelem(lX)
   65 scalar n2 = n*n
   66 bundle ret = null
   67 /*-----------------------------------------------------------------------*/
   68 /* general info                                                          */
   69 /*-----------------------------------------------------------------------*/
   70 /* type goes from 1 to 4 (plain, &quot;C&quot;, &quot;AB&quot; or &quot;SVEC&quot;), plus 10 (&quot;SR&quot;) */
   71 ret.type = type
   72 /* the step we're at: 0 = no estimation done, 1 = VAR only, 2 = SVAR */
   73 ret.step = 0 # maybe not needed any more
   74 matrix mreg = set_default_dimensions(&amp;ret, lY, lX, varorder)
   75 /* the actual data */
   76 matrix ret.Y = mreg[,1:n]
   77 matrix ret.X = (k&gt;0) ? mreg[, n+1 : n+k] : {}
   78 # variable names
   79 strings ret.Ynames = strsplit(strsub(varname(lY), &quot;,&quot; , &quot; &quot;)) # note: string array
   80 strings ret.Xnames = strsplit(strsub(varname(lX), &quot;,&quot; , &quot; &quot;)) # needed?
   81 /* Names for shocks */
   82 # per default, we borrow variable names, following tradition
   83 strings ret.snames = ret.Ynames
   84 /* Don't calculate the long-run matrix by default.
   85 (This does not apply to the case of long-run restrictions.) */
   86 ret.calc_lr = 0
   87 /* Optimisation (not needed for SR, but is expected elsewhere...) */
   88 ret.optmeth = 4 # default = scoring
   89 /* Information for cumulating/normalizing IRFs */
   90 ret.ncumul = 0
   91 matrix ret.cumul = {}
   92 ret.normalize = 0
   93 /* Other initializations */
   94 ret.nboot = 0   # Bootstrap doesn't make sense for SR, but nboot must exist
   95 ret.bestdraw = 0	# only for SR, but must exist, initialize as undefined
   96 ret.checkident = 0  # needs to exist, override below
   97 /* ----------------------------------------------------------------
   98 type-specific stuff
   99 ----------------------------------------------------------------*/
  100 if type != 10    # (traditional, no sign restrictions)
  101   /*
  102   The constraint matrices.
  103   &quot;Rd1&quot; contains short-run constraints on B (and therefore C in non-AB models);
  104   &quot;Rd1l&quot; contains long-run constraints on C (not supported in AB models);
  105   &quot;Rd0&quot; contains short-run constraints on A in AB;
  106   note that &quot;Rd1l&quot; and &quot;Rd0&quot; were both &quot;aux&quot; in previous versions.
  107   Initially, they are empty. Except for the &quot;plain&quot; model, it's up to
  108   the user to fill them up later, via SVAR_restrict() or by hand
  109   */
  110   matrix ret.Rd1 = (type==1) ? cholRd(n) : {}
  111   matrix ret.Rd1l = {}
  112   if type == 3
  113     matrix ret.Rd0 = {}
  114   else
  115     matrix ret.fullRd = {}	# short- and long-run restrictions together later
  116     if type == 4
  117       /* SVEC model: cointegration stuff */
  118       matrix ret.jalpha = {}
  119       matrix ret.jbeta = {}
  120       ret.jcase = 0
  121     endif
  122   endif
  123   ## Other settings for traditional SVARs
  124   /* Copy the id check choice */
  125   ret.checkident = checkid_storeSRirfs
  126   /* Bootstrap */
  127   ret.boot_alpha = -1
  128   matrix ret.bootdata = {}
  129   ret.biascorr = 0
  130   ret.boottype = 1  # standard resid resampl. (choice new 1.5)
  131 else #  10 (SR)
  132   # matrix for holding sign restrictions
  133   matrix ret.SRest = {}
  134   /* Since for sign restrictions the traditional estimation step doesn't apply,
  135   but we need the reduced-form coefficients, we already do it here auto-
  136   matically. */
  137   ret.storeSRirfs = checkid_storeSRirfs
  138   SVAR_estimate(&amp;ret, 0)
  139 endif
  140 return ret
  141 </code>
  142 </gretl-function>
  143 <gretl-function name="SVAR_restrict" type="scalar">
  144  <params count="5">
  145   <param name="b" type="bundleref"/>
  146   <param name="code" type="string"/>
  147   <param name="r" type="int"/>
  148   <param name="c" type="int" default="0"/>
  149   <param name="d" type="scalar" default="0"/>
  150  </params>
  151 <code># c gets a default so that it can be omitted with Adiag, Bdiag (?)
  152 # the d default is also natural
  153 type = b.type
  154 n = b.n
  155 # check input for implemented restriction code
  156 if strstr(&quot;C lrC A B Adiag Bdiag&quot;, code) == &quot;&quot; # code unknown
  157   print &quot;Unknown code in SVAR_restrict.&quot;
  158   return 2
  159   # check for input mismatch
  160 elif (code==&quot;C&quot; || code==&quot;lrC&quot;) &amp;&amp; !(type==1 || type==2 || type==4)
  161   print &quot;C type restriction but not a C model.&quot;
  162   return 1
  163 elif (strstr(&quot;A B Adiag Bdiag&quot;, code) != &quot;&quot;) &amp;&amp; (type!=3)
  164   print &quot;AB type restriction but not an AB model.&quot;
  165   return 1
  166   # check for unsupported case
  167 elif code == &quot;lrC&quot; &amp;&amp; type == 3
  168   print &quot;Long-run restrictions only supported in C model.&quot;
  169   return 3	# (Which code to return here? &lt;Sven&gt;)
  170   # check for bogus long-run constraint
  171 elif code == &quot;lrC&quot; &amp;&amp; type == 4
  172   if c &gt; n - b.crank
  173     print &quot;Long-run constraints only make sense on permanent shocks.&quot;
  174     printf &quot; (First %d shocks.)\n&quot;, n - b.crank
  175     # (p.29, section 7 of the doc)
  176     return 4
  177   endif
  178 elif type == 10
  179   print &quot;Use specialized functions (SVAR_SRplain etc.) for sign restrictions.&quot;
  180   return 5
  181 endif
  182 # if no input error, proceed with this:
  183 err = 0
  184 if code == &quot;C&quot;
  185   matrix Rd = b.Rd1
  186   err = add_constraint(&amp;Rd, n, r, c, d)
  187   if !err
  188     b.Rd1 = Rd
  189   endif
  190 elif code == &quot;lrC&quot;
  191   matrix Rd = b.Rd1l
  192   err = add_constraint(&amp;Rd, n, r, c, d)
  193   if !err
  194     b.Rd1l = Rd
  195   endif
  196 elif code == &quot;A&quot;
  197   matrix Rd = b.Rd0
  198   err = add_constraint(&amp;Rd, n, r, c, d)
  199   if !err
  200     b.Rd0 = Rd
  201   endif
  202 elif code == &quot;B&quot;
  203   matrix Rd = b.Rd1
  204   err = add_constraint(&amp;Rd, n, r, c, d)
  205   if !err
  206     b.Rd1 = Rd
  207   endif
  208 elif code == &quot;Adiag&quot;
  209   matrix Rd = b.Rd0
  210   if ok(r)
  211     b.Rd0 = Rd | diag_Rd(n, r)
  212   else
  213     b.Rd0 = Rd | free_diag_Rd(n)
  214   endif
  215 elif code == &quot;Bdiag&quot;
  216   matrix Rd = b.Rd1
  217   if ok(r)
  218     b.Rd1 = Rd | diag_Rd(n, r)
  219   else
  220     b.Rd1 = Rd | free_diag_Rd(n)
  221   endif
  222 endif
  223 ## Inform the user if the restriction failed.
  224 /*
  225 At this point it should hold that:
  226 err == 0 : add_constraint worked ok, -- or Adiag/Bdiag: is it conceivable that this happens
  227 together with other A/B restrictions? Then it probably should
  228 also be checked in principle (but doesn't happen here).
  229 */
  230 if err == -1
  231   printf &quot;Imposing restriction failed, bad input to &quot;
  232   printf &quot;add_constraint.\n&quot;
  233 elif err == 10
  234   printf &quot;Imposing restriction failed, conflicting with &quot;
  235   printf &quot;earlier restrictions.\n&quot;
  236 elif err == 20
  237   printf &quot;Imposing restriction failed, redundant.\n&quot;
  238 endif
  239 if err
  240   printf &quot;(Code %s, &quot;, code
  241   if ok(r)
  242     printf &quot;element %d,%d restricted to %f.)\n&quot;, r,c,d
  243   else
  244     printf &quot;no manual restriction.)\n&quot;
  245   endif
  246 endif
  247 return err	# 0, -1, 10, or 20
  248 </code>
  249 </gretl-function>
  250 <gretl-function name="SVAR_ident" type="scalar">
  251  <params count="2">
  252   <param name="b" type="bundleref"/>
  253   <param name="verbose" type="int" default="0"/>
  254  </params>
  255 <code>/* Apparently returns a non-zero number if identification holds.
  256 (Coming from ident(); Which is perhaps against convention, usually 0 for &quot;no error&quot;?)
  257 */
  258 if verbose
  259   print &quot;Check for identification&quot;
  260   print &quot;------------------------&quot;
  261   print &quot; (It may happen that in special circumstances under-&quot;
  262   print &quot;  identification goes unnoticed. If estimation fails,&quot;
  263   print &quot;  perhaps try adding restrictions.)&quot;
  264 endif
  265 n2 = b.n * b.n
  266 if verbose
  267   printf &quot;\nConstraints in implicit form:\n&quot;
  268   if b.type == 3
  269     printf &quot; Ra: %4.0f&quot;, b.Rd0[,1:n2]
  270     printf &quot; da_T: %4.0f&quot;, b.Rd0[,n2+1]'
  271     printf &quot; Rb: %4.0f&quot;, b.Rd1[,1:n2]
  272     printf &quot; db_T: %4.0f&quot;, b.Rd1[,n2+1]'
  273   else    # The A restrictions are trivial and confusing here.
  274     matrix b.fullRd = get_full_Rd(&amp;b, 0)
  275     # (maybe change get_full_Rd so that it saves its result directly in the model
  276     #  bundle, just like it does with C1...)
  277     printf &quot; Rc: %4.0f&quot;, b.fullRd[,1:n2]
  278     printf &quot; dc_T: %4.0f&quot;, b.fullRd[,n2+1]'
  279   endif
  280   printf &quot;\n&quot;
  281 endif
  282 ret = ident(&amp;b, verbose) # was: &amp;Ra, &amp;da, &amp;Rb, &amp;db, verbose)
  283 return ret
  284 </code>
  285 </gretl-function>
  286 <gretl-function name="SVAR_estimate" type="scalar">
  287  <params count="2">
  288   <param name="obj" type="bundleref"/>
  289   <param name="verbosity" type="int" default="1"/>
  290  </params>
  291 <code>/*
  292 this function fills the bundle with the estimated structural
  293 matrices and the covariance matrix of their free elements; it also
  294 calls do_IRF at the end so that the structural VMA is stored into
  295 the bundle
  296 */
  297 scalar type = obj.type
  298 scalar meth = obj.optmeth
  299 scalar n = obj.n
  300 scalar T = obj.T
  301 matrix vcv
  302 scalar errcode = 0
  303 if type == 4 # do VECM
  304   if inbundle(obj, &quot;cointsetup&quot;)
  305     if obj.cointsetup
  306       vecm_est(&amp;obj)
  307     else
  308       # Not sure if this case (obj.cointsetup==0)can currently happen.
  309       funcerr &quot;No valid cointegration setup found.&quot;
  310     endif
  311   else
  312     funcerr &quot;Need to do cointegration setup first for 'SVEC'.&quot;
  313   endif
  314 else # estimate ordinary VAR
  315   base_est(&amp;obj)
  316 endif
  317 if obj.checkident
  318   id = SVAR_ident(&amp;obj, (verbosity &gt; 1)) # output for verbosity &gt;= 2
  319   if !id
  320     funcerr &quot;Identification check failed.&quot;
  321   endif
  322 endif
  323 # grab the instantaneous covariance matrix
  324 matrix Sigma = obj.Sigma
  325 scalar obj.LL0 = VARloglik(obj.T, Sigma)
  326 if type == 1
  327   # plain model: just Cholesky decomposition
  328   matrix C = cholesky(Sigma)
  329   matrix param = vech(C')
  330   # compute the covariance matrix for C
  331   matrix Ss = imp2exp(obj.Rd1)
  332   vcv = coeffVCV(Ss[, 1: cols(Ss)-1], &amp;C)
  333 elif (type==2) || (type == 4)
  334   # C models in a broad sense (including SVEC)
  335   # Maybe redundant restrictions have already been detected
  336   # by SVAR_ident and then removed; in this case we mustn't
  337   # use the plain vanilla fullRd.
  338   if inbundle(obj, &quot;cleanfullRd&quot;)
  339     matrix fullRd = obj.cleanfullRd
  340     # (obj.C1 and obj.fullRd already created by SVAR_ident then)
  341   elif inbundle(obj, &quot;fullRd&quot;) &amp;&amp; inbundle(obj, &quot;C1&quot;)
  342     matrix fullRd = obj.fullRd
  343   else
  344     matrix fullRd = get_full_Rd(&amp;obj, verbosity)
  345     # (This also sets obj.C1.)
  346     matrix obj.fullRd = fullRd
  347   endif
  348   # try to set some &quot;sensible&quot; initial values
  349   matrix param = init_C(Sigma, fullRd)
  350   # do estimation; note that vcv is estimated inside &quot;estC&quot;
  351   matrix C = estC(&amp;param, Sigma, fullRd, &amp;vcv, &amp;errcode, meth, verbosity)
  352 elif type == 3
  353   # AB-model
  354   matrix E = obj.E     # grab the VAR residuals (needed for initialisation)
  355   matrix bRd = obj.Rd1 # restrictions on B
  356   matrix aRd = obj.Rd0 # restrictions on A
  357   # try to set some &quot;sensible&quot; initial values
  358   # (substitute out call to (former) init_AB)
  359   matrix param = is_standard_AB(aRd, bRd) ? stdAB_init(E, aRd, bRd) : nonstdAB_init(E, aRd, bRd)
  360   # do estimation; note that vcv is estimated inside &quot;estAB&quot;
  361   # (no it actually wasn't, now it is &lt;Sven&gt;)
  362   matrices transfer
  363   matrix C = estAB(&amp;param, Sigma, aRd, bRd, &amp;vcv, &amp;errcode, meth, verbosity, &amp;transfer)
  364 elif type == 10    # SR
  365   # just enter some dummy stuff, only needed for compatibility right now
  366   matrix C = cholesky(Sigma)
  367   matrix obj.theta = vech(C')
  368   # the covariance matrix for C (just dummy!)
  369   matrix Ss = imp2exp(cholRd(n))
  370   matrix obj.vcv = coeffVCV(Ss[, 1: cols(Ss)-1], &amp;C)
  371   matrix obj.C = C
  372 endif	# which type
  373 /*
  374 Post-estimation; transfers, long-run matrix, over-id test
  375 */
  376 if errcode
  377   funcerr &quot;Estimation failed&quot;
  378 elif type != 10	# estimation ran fine, and we don't have a SR model
  379   # Copy stuff into the bundle
  380   if type == 3
  381     matrix obj.S1 = transfer[1]	# A
  382     matrix obj.S2 = transfer[2]	# B
  383     scalar obj.ka = transfer[3]
  384     scalar obj.kb = transfer[4]
  385   endif
  386   /*
  387   We now also store C in the bundle because its computation
  388   was repeated several times elsewhere
  389   (This was obj.S1 for C models instead of obj.C, no idea why.) */
  390   matrix obj.C = C
  391   matrix obj.theta = param
  392   matrix obj.vcv = vcv
  393   # Long-run matrix
  394   # Jan 2018: add the type 4 possibility
  395   if ( (type &lt; 3) &amp;&amp; ( rows(obj.Rd1l) || obj.calc_lr ) ) || type == 4
  396     # a plain or C model with long-run constraints, or user switch;
  397     # here we now (Oct 2017) calc and save the long-run matrix
  398     # re-use the C1 matrix from above if possible
  399     matrix C1 = (type == 2 || type == 4) ? obj.C1 : C1mat(obj.VARpar)
  400     matrix obj.lrmat = C1 * obj.C
  401   endif
  402   # Indicate that estimation is done (still needed?)
  403   obj.step = 2
  404   # Store IRFs into the bundle
  405   doIRF(&amp;obj)
  406   # store the log-likelihood into the bundle
  407   scalar obj.LL1 = VARloglik(obj.T, obj.Sigma, &amp;C)
  408   # calculate the over-id test in any case (not just for verbosity)
  409   # (C should hopefully be correctly depending on type)
  410   overid = (n * (n+1) / 2 - rows(obj.theta))
  411   if overid &gt; 0
  412     LR = 2 * (obj.LL0 - obj.LL1)
  413     matrix obj.LRoid = {LR; overid; pvalue(X, overid, LR)}
  414     # this is: (stat| dof| pv)
  415   endif
  416   if (verbosity &gt; 0)
  417     SVAR_est_printout(&amp;obj)
  418   endif
  419 endif
  420 return errcode
  421 </code>
  422 </gretl-function>
  423 <gretl-function name="SVAR_cumulate" type="scalar">
  424  <params count="2">
  425   <param name="b" type="bundleref"/>
  426   <param name="nv" type="int"/>
  427  </params>
  428 <code>err = (nv&gt;b.n || nv&lt;1) # was nv&lt;0, but 0 makes no sense (?)
  429 if !err
  430   vn = b.Ynames
  431   printf &quot;Variable %s cumulated\n&quot;,  vn[nv]
  432   b.cumul = b.cumul | {nv}
  433   b.ncumul = b.ncumul + 1
  434   # The following code was in doIRF but doesn't depend
  435   # on the data, so slightly inefficient for the bootstrap.
  436   # Instead we introduce b.cumsel which can be reused later
  437   tmp = zeros(b.n, b.n)
  438   tmp[b.cumul,] = 1
  439   matrix b.cumsel = selifr(transp(seq(1, b.n * b.n)), vec(tmp))
  440 endif
  441 return err
  442 </code>
  443 </gretl-function>
  444 <gretl-function name="SVAR_boot" type="scalar">
  445  <params count="6">
  446   <param name="obj" type="bundleref"/>
  447   <param name="rep" type="int" min="0" default="2000">
  448 <description>bootstrap iterations</description>
  449   </param>
  450   <param name="alpha" type="scalar" min="0" max="1" default="0.9">
  451 <description>CI coverage</description>
  452   </param>
  453   <param name="quiet" type="bool" default="1"/>
  454   <param name="btypestr" type="string" optional="true">
  455 <description>bootstrap type</description>
  456   </param>
  457   <param name="biascorr" type="int" min="-1" max="2" default="-1">
  458 <description>bias correction (non-SVEC)</description>
  459   </param>
  460  </params>
  461 <code># btypestr: can be &quot;resample&quot; / &quot;resampling&quot;,
  462 #  &quot;wildN&quot;/&quot;wild&quot;, &quot;wildR&quot;, &quot;wildM&quot;
  463 # The default value for biascorr of -1 means:
  464 # Do not override the previous setting.
  465 ## Copy some params and choices
  466 loop foreach i n k T p type --quiet
  467   scalar $i = obj.$i
  468 endloop
  469 obj.nboot = rep         # record bootstrap details
  470 obj.boot_alpha = alpha  # into original model
  471 if type == 10
  472   funcerr &quot;Wrong turn: Set-ID not for bootstrapping...&quot;
  473 endif
  474 ## Bootstrap type choice (if different from default)
  475 if exists(btypestr)
  476   boottypechoice(&amp;obj, btypestr)
  477 endif
  478 # Copy optional block length choice
  479 bundle bparams = null
  480 if obj.boottype == 5 &amp;&amp; inbundle(obj, &quot;moveblocklen&quot;)
  481   bparams.moveblocklen = obj.moveblocklen
  482 endif
  483 # Bias correction choice, leave or update?
  484 obj.biascorr = (biascorr == -1) ? obj.biascorr : biascorr
  485 # define default for bias correction iterations
  486 if obj.biascorr &amp;&amp; !inbundle(obj, &quot;BCiter&quot;)
  487   obj.BCiter = 1024
  488 endif
  489 ## Various needed stuff
  490 if type == 3 # AB
  491   matrix bmA bmB # needed as memory for transfer
  492 elif type == 4
  493   matrix J = zeros(n - obj.crank, obj.crank) | I(obj.crank)
  494 endif
  495 matrix start = obj.Y[1:p, ] # Y0
  496 # disentangle determ/exog:
  497 calc_bmu(&amp;obj)	# adds obj.bmu
  498 # store a copy of the model for bootstrap
  499 bundle bobj = obj
  500 # Do the bias correction pre-step if applicable
  501 maybe_do_biascorr(&amp;bobj, bparams)
  502 printf &quot;\nBootstrapping model (%d iterations)\n&quot;, rep
  503 printf &quot;Bootstrap type: %s\n&quot;, btypestring(obj.boottype)
  504 printf &quot;Bias correction: %s\n\n&quot;, BCstring(obj.biascorr)
  505 flush
  506 ## Actual bootstrap simulation
  507 matrices bootout = SVAR_boot_innerloop(&amp;bobj, obj, bparams)
  508 /*
  509 &quot;bootirfs&quot;:
  510 each bootstrap replication on one row; each row contains
  511 the vectorisation of the complete IRF matrix
  512 */
  513 matrix bootirfs = bootout[1] # zeros(rep, (h+1) * n2)
  514 # Spar_mat is probably somewhat redundant..., only for the printout
  515 matrix Spar_mat = bootout[2] # zeros(rep, n2) or zeros(rep, 2*n2)
  516 scalar failed   = bootout[3]
  517 if !quiet
  518   boot_printout(type, n, rep, failed, Spar_mat)
  519 endif
  520 # quantiles of bootstrapped IRFs used in graphs
  521 q_alpha = 0.5 * (1 - alpha)	# changed in v1.5
  522 matrix locb = quantile(bootirfs, q_alpha)
  523 matrix hicb = quantile(bootirfs, 1 - q_alpha)
  524 matrix mdn  = quantile(bootirfs, 0.5)
  525 bundle bootdata = null
  526 bootdata.rep   = rep                # no of replications
  527 bootdata.alpha = alpha              # alpha
  528 bootdata.biascorr  = obj.biascorr   # type of bias correction
  529 scalar h = obj.horizon
  530 scalar n2 = n*n
  531 matrix bootdata.lo_cb = mshape(locb, h+1, n2) # lower bounds
  532 matrix bootdata.hi_cb = mshape(hicb, h+1, n2) # upper bounds
  533 matrix bootdata.mdns  = mshape(mdn, h+1, n2)  # medians
  534 bundle obj.bootdata = bootdata
  535 return failed
  536 </code>
  537 </gretl-function>
  538 <gretl-function name="SVAR_hd" type="list">
  539  <params count="3">
  540   <param name="Mod" type="bundleref"/>
  541   <param name="nv" type="int" min="1" default="1"/>
  542   <param name="drawix" type="int" min="0" default="0"/>
  543  </params>
  544 <code># historical decomposition
  545 # (drawix only meant for the set id case (type 10))
  546 list ret = null
  547 loop foreach i n p t1 t2 type T k --quiet
  548   scalar $i = Mod.$i
  549 endloop
  550 matrix B10 = {} # to be filled for type 10
  551 if nv &gt; n
  552   printf &quot;Hm. There are %d variables in the model. &quot;, n
  553   printf &quot;Chosen shock index: %d\n&quot;, nv
  554   funcerr &quot;Shock index out of range&quot;
  555   ## (further range check for SR partial id below) ##
  556 endif
  557 # Prepare the set id case
  558 if type == 10
  559   errchkSRhisto(&amp;Mod, drawix)
  560   # allow drawix to override the setting in the bundle
  561   whichdraw = drawix ? drawix : Mod.bestdraw
  562   bundle pickdraw = Mod.acc_draws[whichdraw]
  563   matrix B10 = pickdraw.B
  564 endif
  565 # The following might be redrawn in type10/Bayesian
  566 matrices muVARparE = muVARparE_mayberedr(Mod, B10)
  567 # compute the exogenous part
  568 if type &lt; 4
  569   matrix m = Mod.X * Mod.mu
  570 elif type == 4
  571   # here we have to take into account the &quot;5 cases&quot;
  572   dcase = Mod.jcase
  573   # T     = Mod.T
  574   matrix mreg = (dcase == 1) ? {} : ones(T,1)
  575   if dcase &gt; 3
  576     mreg ~= seq(1,T)'
  577   endif
  578   matrix m = (mreg ~ Mod.X) * Mod.mu
  579 elif type == 10
  580   matrix m = Mod.X * muVARparE[1]
  581 endif
  582 # grab the C matrix
  583 if (type == 1) || (type == 2) || (type == 4)
  584   matrix C = Mod.C
  585 elif type == 3
  586   if inbundle(Mod, &quot;C&quot;)
  587     matrix C = Mod.C
  588   else
  589     matrix C = Mod.S1 \ Mod.S2
  590   endif
  591 elif type == 10
  592   matrix C = pickdraw.irfs[1] # impact effect is C
  593   if cols(C) &lt; Mod.n
  594     funcerr &quot;partial id not supported for historical decomp&quot;
  595   endif
  596 endif
  597 matrix iC = inv(C)
  598 strings Ynames = Mod.Ynames
  599 strings snames = Mod.snames
  600 string yn = Ynames[nv]
  601 smpl t1 t2
  602 if cols(m)&gt;0
  603   Xdet = varsimul(muVARparE[2], m[p+1:,], Mod.Y[1:p,]) # was VARpar
  604 else
  605   Xdet = varsimul(muVARparE[2], zeros(T-p, n), Mod.Y[1:p,]) # was Mod.T
  606 endif
  607 ret += genseries( sprintf(&quot;hd_%s_det&quot;, yn), Xdet[,nv])
  608 # the structural shocks
  609 matrix U = muVARparE[3] * iC' # was E
  610 rotVARpar = iC * muVARparE[2] * (I(p) ** C)
  611 loop i = 1..n --quiet
  612   a = (seq(1,n) .= i)
  613   W = varsimul(rotVARpar, U .* a, zeros(p,n)) * C'
  614   ret += genseries(sprintf(&quot;hd_%s_%s&quot;, yn,    fixname(snames[i])), W[,nv])
  615 endloop
  616 return ret
  617 </code>
  618 </gretl-function>
  619 <gretl-function name="SVAR_coint" type="scalar">
  620  <params count="6">
  621   <param name="SVARobj" type="bundleref"/>
  622   <param name="dcase" type="int" min="1" max="5" default="3"/>
  623   <param name="jbeta" type="matrix"/>
  624   <param name="jalpha" type="matrix" optional="true"/>
  625   <param name="verbose" type="bool" default="0"/>
  626   <param name="rexo" type="list" optional="true" const="true"/>
  627  </params>
  628 <code>/*
  629 This function doesn't do very much, except setting
  630 up the model for subsequent VECM estimation; &quot;dcase&quot; tells you
  631 which of the &quot;five cases&quot; we want (no constant, restricted
  632 constant, etc), jbeta is simply checked for dimensions and then
  633 copied into the object.
  634 As for jalpha, if it's an empty matrix, that means &quot;just estimate it unrestrictedly&quot;, and we set up a
  635 flag accordingly. Otherwise, it's taken to be pre-set to some
  636 fixed value; the intermediate case (contraints on alpha) is not
  637 handled, and I doubt it will ever be.
  638 While we're at it, we also label the structural shocks as
  639 &quot;Perm_1&quot;, &quot;Perm_2&quot;, &quot;Trans_1&quot;, &quot;Trans_2&quot;, etc.
  640 */
  641 if SVARobj.type != 4 &amp;&amp; verbose
  642   print &quot;Putting cointegration information into a non-SVEC model --&quot;
  643   print &quot; are you sure you know what you're doing?&quot;
  644 endif
  645 scalar n = SVARobj.n
  646 # define default
  647 if !exists(jalpha)
  648   matrix jalpha = {}
  649 endif
  650 # syntax check
  651 err = (dcase&lt;1) || (dcase&gt;5)
  652 if err
  653   errmsg = sprintf(&quot;Invalid dcase value %d&quot;, dcase)
  654   funcerr errmsg
  655 endif
  656 # dimensions check
  657 if dcase%2 # nice, huh?
  658   err = err || (rows(jbeta) != n)
  659 else
  660   err = err || (rows(jbeta) != n+1)
  661 endif
  662 if err
  663   errmsg = sprintf(&quot;jbeta: has %d rows, should have %d&quot;, rows(jbeta), dcase%2 ? n : n+1)
  664   funcerr errmsg
  665 endif
  666 r = cols(jbeta)
  667 err = err || (n &lt; r) # should this be &lt;=? hmm.
  668 # rank check
  669 err = err || (rank(jbeta) &lt; r)
  670 # now check if alpha is ok
  671 d = rows(jalpha)
  672 # d==0 is ok, we'll estimate alpha later
  673 free_a = (d==0)
  674 if !free_a
  675   err = err || (d != n) || (cols(jalpha) != r)
  676 endif
  677 # if anything goes wrong, return
  678 if err
  679   # (currently err==1 here always, Sven believes)
  680   outfile stderr --write --quiet
  681   printf &quot;SVAR_coint: returning on err = %d\n&quot;, err
  682   outfile --close
  683   return err
  684 endif
  685 # fill up the object with the info
  686 SVARobj.crank = r
  687 SVARobj.jcase = dcase
  688 SVARobj.jbeta = jbeta
  689 SVARobj.jalpha = jalpha
  690 SVARobj.free_a = free_a
  691 if verbose
  692   if dcase == 1
  693     printf &quot;No constant, &quot;
  694   elif dcase == 2
  695     printf &quot;Restricted constant, &quot;
  696   elif dcase == 3
  697     printf &quot;Unrestricted constant, &quot;
  698   elif dcase == 4
  699     printf &quot;Restricted trend, &quot;
  700   elif dcase == 5
  701     printf &quot;Unrestricted trend, &quot;
  702   endif
  703   printf &quot;beta =\n%11.5f\n&quot;, jbeta
  704   if free_a
  705     printf &quot;alpha is unrestricted\n&quot;
  706   else
  707     printf &quot;alpha =\n%9.5f\n&quot;, jalpha
  708     printf &quot;PI =\n%9.5f\n&quot;, jalpha * jbeta'
  709   endif
  710 endif
  711 # relabel transitory structural shocks
  712 strings sn = SVARobj.snames
  713 if n-r == 1
  714   sn[1] = sprintf(&quot;Permanent&quot;)
  715   if r == 1
  716     sn[2] = &quot;Transitory&quot;
  717   else
  718     loop i=2..n --quiet
  719       sn[i] = sprintf(&quot;Transitory_%d&quot;, i-1)
  720     endloop
  721   endif
  722 else
  723   loop i=1..n --quiet
  724     sn[i] = sprintf(&quot;Permanent_%d&quot;, i)
  725   endloop
  726   if r == 1
  727     sn[n] = &quot;Transitory&quot;
  728   else
  729     loop i=1..r --quiet
  730       sn[n-r+i] = sprintf(&quot;Transitory_%d&quot;, i)
  731     endloop
  732   endif
  733 endif
  734 SVARobj.snames = sn
  735 if !err
  736   SVARobj.cointsetup = 1	# flag success
  737 endif
  738 return err
  739 </code>
  740 </gretl-function>
  741 <gretl-function name="GetShock" type="series">
  742  <params count="3">
  743   <param name="SVARobj" type="bundleref"/>
  744   <param name="i" type="int" min="1" default="1"/>
  745   <param name="drawix" type="int" min="0" default="0"/>
  746  </params>
  747 <code>/*
  748 Produces the series corresponding to the historical shock
  749 realizations associated with the point estimates of the model
  750 (and IRFs).
  751 For set identification (sign restrictions) there is no point
  752 estimate; however, we support that
  753 the user picks one of the accepted draws and then the shock series
  754 is based on that particular model draw.
  755 # (drawix only meant for the set id case (type 10))
  756 */
  757 series ret = NA
  758 type = SVARobj.type
  759 matrix B10 = {}	# to be filled in type 10
  760 ## some error checks ##
  761 if type &gt; 4 &amp;&amp; type != 10
  762   printf &quot;Given type %d\n&quot;, type
  763   funcerr &quot;Unknown model type&quot;
  764 elif i &gt; SVARobj.n
  765   printf &quot;Chosen shock index: %d\n&quot;, i
  766   funcerr &quot;Shock index out of range&quot;
  767 elif type != 10 &amp;&amp; drawix &gt; 0
  768   print &quot;Warning: 'drawix' arg meaningless for standard SVAR&quot;
  769 elif type == 10
  770   errchkSRhisto(&amp;SVARobj, drawix)
  771 endif
  772 ## get the C matrix (and then the inv) ##
  773 if (type == 1) || (type == 2) || (type == 4)
  774   matrix C = SVARobj.C
  775 elif type == 3
  776   if inbundle(SVARobj, &quot;C&quot;)	# maybe not yet computed
  777     matrix C = SVARobj.C
  778   else
  779     matrix C = SVARobj.S1 \ SVARobj.S2
  780   endif
  781 elif type == 10 # set id
  782   # allow drawix to override the setting in the bundle
  783   whichdraw = drawix ? drawix : SVARobj.bestdraw
  784   bundle pickdraw = SVARobj.acc_draws[whichdraw]
  785   matrix C = pickdraw.irfs[1] # impact effect is C
  786   B10 = pickdraw.B
  787   if cols(C) &lt; SVARobj.n
  788     funcerr &quot;partial id not supported for shock retrieval&quot;
  789   endif
  790 endif
  791 matrix iC = inv(C')
  792 matrix resids = muVARparE_mayberedr(SVARobj, B10)[3]
  793 ## construct the wanted series ##
  794 extra = $nobs - rows(resids)
  795 matrix tmp = {}
  796 if extra &gt; 0
  797   set warnings off
  798   tmp = ones(extra,1) .* NA
  799 endif
  800 tmp |= resids * iC[,i]
  801 ret = tmp
  802 snames = SVARobj.snames # strings array?
  803 string vlab = snames[i]
  804 setinfo ret --description=&quot;@vlab&quot;
  805 return ret
  806 </code>
  807 </gretl-function>
  808 <gretl-function name="IRFplot" type="void">
  809  <params count="5">
  810   <param name="obj" type="bundleref"/>
  811   <param name="snum" type="int" default="1"/>
  812   <param name="vnum" type="int" default="1"/>
  813   <param name="keypos" type="int" min="0" max="2" default="1"/>
  814   <param name="drawix" type="int" min="0" default="0"/>
  815  </params>
  816 <code>IRFsave(&quot;display&quot;, &amp;obj, snum, vnum, keypos, drawix)
  817 </code>
  818 </gretl-function>
  819 <gretl-function name="IRFsave" type="void">
  820  <params count="6">
  821   <param name="outfilename" type="string"/>
  822   <param name="obj" type="bundleref"/>
  823   <param name="snum" type="int"/>
  824   <param name="vnum" type="int"/>
  825   <param name="keypos" type="int" min="0" max="2" default="1"/>
  826   <param name="drawix" type="int" min="0" default="0"/>
  827  </params>
  828 <code># negative snum is allowed and means to flip the shock
  829 # copy and prepare some input
  830 string tmpfile, tmpout
  831 n = obj.n
  832 whichdraw = drawix ? drawix : obj.bestdraw
  833 if obj.type == 10	# SR
  834   errchkSRhisto(&amp;obj, drawix)
  835   putIrf_to_accdraw(&amp;obj, whichdraw)
  836   matrix IRFmat = obj.acc_draws[whichdraw].IRFs
  837 else
  838   matrix IRFmat = obj.IRFs
  839 endif
  840 scale = 1	# possibly changed later
  841 if obj.normalize == 1
  842   matrix tmp = mshape(IRFmat[1,], n, n)
  843 endif
  844 if obj.nboot
  845   boot = obj.bootdata	# bundle?
  846   bc = boot.biascorr
  847 endif
  848 scalar sfrom sto vfrom vto
  849 is_srange = range(snum, n, &amp;sfrom, &amp;sto)
  850 is_vrange = range(vnum, n, &amp;vfrom, &amp;vto)
  851 # cycle through all (selected) shocks and variables
  852 loop snum = sfrom..sto -q
  853   # do checks
  854   err = check_bounds(snum, 1, n)
  855   if err == 1
  856     printf &quot;Shock number %d out of bounds\n&quot;, abs(snum)
  857     # (abs because of flipping)
  858     return
  859   endif
  860   string sn = obj.snames[abs(snum)]
  861   # normalization / scaling
  862   if obj.normalize == 1
  863     scale = tmp[abs(snum), abs(snum)]
  864   endif
  865   loop vnum = vfrom..vto -q
  866     # do checks
  867     err = check_bounds(1, vnum, n)
  868     if err == 2
  869       printf &quot;Variable number %d out of bounds\n&quot;, vnum
  870       return
  871     endif
  872     string yn = obj.Ynames[vnum]
  873     # produce plots
  874     if obj.ncumul == 0
  875       if obj.nboot == 0
  876         tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, null, null, , whichdraw)
  877       else
  878         tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, null, &amp;boot, bc, whichdraw)
  879       endif
  880     else
  881       matrix cumul = obj.cumul
  882       if obj.nboot == 0
  883         tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, &amp;cumul, null, , whichdraw)
  884       else
  885         tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, &amp;cumul, &amp;boot, bc, whichdraw)
  886       endif
  887     endif
  888     if (outfilename == &quot;display&quot;) || (sfrom == sto &amp;&amp; vfrom == vto)
  889       # (single plot, no indices)
  890       tmpout = outfilename
  891     else
  892       strings be = basename(outfilename)	# gives 2-elem array
  893       if is_vrange	 # several v indices
  894         if is_srange # several s indices
  895           tmpout = sprintf(&quot;%s_%d%d.%s&quot;, be[1], snum, vnum, be[2])
  896         else
  897           tmpout = sprintf(&quot;%s_%d.%s&quot;, be[1], vnum, be[2])
  898         endif
  899       elif is_srange 	# several s indices
  900         tmpout = sprintf(&quot;%s_%d.%s&quot;, be[1], snum, be[2])
  901       else
  902         funcerr &quot;shouldn't happen&quot;
  903       endif
  904     endif
  905     gnuplot --input=&quot;@tmpfile&quot; --output=&quot;@tmpout&quot;
  906   endloop
  907 endloop
  908 </code>
  909 </gretl-function>
  910 <gretl-function name="FEVD" type="matrix">
  911  <params count="2">
  912   <param name="SVARobj" type="bundleref"/>
  913   <param name="drawix" type="int" min="0" default="0"/>
  914  </params>
  915 <code># (drawix only meant for the set id case (type 10))
  916 n = SVARobj.n
  917 h = SVARobj.horizon + 1
  918 if SVARobj.type == 10
  919   # In the set id case in each accepted draw the impulse
  920   # responses are already stored as &quot;irfs&quot;; however, the format
  921   # there is an array of matrices.
  922   errchkSRhisto(&amp;SVARobj, drawix)
  923   # allow drawix to override the setting in the bundle
  924   whichdraw = drawix ? drawix : SVARobj.bestdraw
  925   bundle pickdraw = SVARobj.acc_draws[whichdraw]
  926   if cols(pickdraw.irfs[1]) != n
  927     funcerr &quot;partial id not supported for FEVD&quot;
  928   elif h != nelem(pickdraw.irfs)
  929     funcerr &quot;horizon mismatch&quot;
  930   endif
  931   if !inbundle(pickdraw, &quot;IRFs&quot;) # maybe have already been added there
  932     putIrf_to_accdraw(&amp;SVARobj, whichdraw)
  933     # matrix IRFs = zeros(h, n*n)
  934     # loop ix = 1..h -q
  935     #     IRFs[ix, ] = vec(pickdraw.irfs[ix])'
  936     # endloop
  937     ## copy to origin
  938     # matrix SVARobj.acc_draws[whichdraw].IRFs = IRFs
  939     # else
  940     #     matrix IRFs = pickdraw.IRFs
  941   endif
  942   matrix IRFs = SVARobj.acc_draws[whichdraw].IRFs
  943 else # standard non-SR model
  944   if drawix &gt; 0
  945     print &quot;Warning: 'drawix' arg meaningless for standard SVAR, ignoring&quot;
  946   endif
  947   matrix IRFs = SVARobj.IRFs
  948 endif
  949 matrix ret = zeros(h, n*n)
  950 ctmp = cum(IRFs .* IRFs)
  951 loop i = 1..h --quiet
  952   tmp = mshape(ctmp[i,],n,n)'
  953   ret[i,] = vec(tmp ./ sumc(tmp))'
  954 endloop
  955 return ret
  956 </code>
  957 </gretl-function>
  958 <gretl-function name="GUI_SVAR" type="bundle" pkg-role="gui-main">
  959  <params count="17">
  960   <param name="type" type="int" min="1" max="3" default="1">
  961 <description>Model type</description>
  962 <labels count="3">
  963 "plain (Cholesky)" "C-model" "AB-model" </labels>
  964   </param>
  965   <param name="Y" type="list">
  966 <description>VAR variables</description>
  967   </param>
  968   <param name="X" type="list" optional="true">
  969 <description>Exogenous regressors</description>
  970   </param>
  971   <param name="hasconst" type="bool" default="1">
  972 <description>Constant</description>
  973   </param>
  974   <param name="hastrend" type="bool" default="0">
  975 <description>Time trend</description>
  976   </param>
  977   <param name="hasseas" type="bool" default="0">
  978 <description>Seasonal dummies</description>
  979   </param>
  980   <param name="l" type="int" min="1" default="1">
  981 <description>Lags</description>
  982   </param>
  983   <param name="h" type="int" min="0">
  984 <description>Horizon</description>
  985   </param>
  986   <param name="R1" type="matrixref" optional="true">
  987 <description>Restriction pattern (short-run C or B)</description>
  988   </param>
  989   <param name="R2" type="matrixref" optional="true">
  990 <description>Restriction pattern (long-run C or A)</description>
  991   </param>
  992   <param name="b" type="int" min="0">
  993 <description>Bootstrap replications</description>
  994   </param>
  995   <param name="alpha" type="scalar" min="0" max="1" default="0.9">
  996 <description>Bootstrap alpha</description>
  997   </param>
  998   <param name="biascorr" type="int" min="0" max="2" default="0">
  999 <description>Bias correction</description>
 1000 <labels count="3">
 1001 "None" "Partial" "Full" </labels>
 1002   </param>
 1003   <param name="boottype" type="int" min="1" max="5" default="1">
 1004 <description>Bootstrap type</description>
 1005 <labels count="5">
 1006 "residual resampling" "wild/Normal" "wild/Rademacher" "wild/Mammen" "moving blocks" </labels>
 1007   </param>
 1008   <param name="checkident" type="bool" default="0">
 1009 <description>Check identification</description>
 1010   </param>
 1011   <param name="cumix" type="matrixref" optional="true">
 1012 <description>Indices of responses to cumulate</description>
 1013   </param>
 1014   <param name="optmeth" type="int" min="0" max="4" default="4">
 1015 <description>Optimization method</description>
 1016 <labels count="5">
 1017 "BFGS (numerical score)" "BFGS (analytical score)" "Newton-Raphson (numerical score)" "Newton-Raphson (analytical score)" "Scoring algorithm" </labels>
 1018   </param>
 1019  </params>
 1020 <code>n = nelem(Y)
 1021 # stick together deterministics and other exog.
 1022 list lX = dropcoll(determ(X, hasconst, hastrend, hasseas))
 1023 # initialize the model bundle
 1024 bundle m = SVAR_setup(modelstring(type), Y, lX, l, checkident)
 1025 if h &gt; 0
 1026   m.horizon = h
 1027 endif
 1028 # copy options and/or overwrite defaults
 1029 m.biascorr = biascorr
 1030 m.optmeth = optmeth
 1031 m.boottype = boottype
 1032 ## implement the cumulation spec
 1033 if exists(cumix)
 1034   # ensure column vector
 1035   cumix = vec(cumix)
 1036   # input checks (numbers out of bounds)
 1037   if max(cumix) &gt; nelem(Y) || min(cumix) &lt; 1
 1038     print &quot;Invalid cumulation specification!&quot;
 1039     print &quot;(No responses will be cumulated.)&quot;
 1040   else	# sensible cumulation spec
 1041     loop i=1..rows(cumix) -q
 1042       SVAR_cumulate(&amp;m, cumix[i])
 1043     endloop
 1044   endif
 1045 endif
 1046 ## process restrictions
 1047 if type == 1
 1048   if exists(R1) || exists(R2)
 1049     print &quot;Estimating plain model. Discarding provided restrictions.&quot;
 1050   endif
 1051   if checkident
 1052     print &quot;(Identification trivially given in plain model.)&quot;
 1053   endif
 1054 else # C or AB model
 1055   # input check
 1056   if !exists(R1) &amp;&amp; !exists(R2)
 1057     funcerr &quot;Must provide some restrictions for C and AB models!&quot;
 1058   endif
 1059   if type == 3 &amp;&amp; ( !exists(R1) || !exists(R2) )
 1060     funcerr &quot;Must provide restrictions on A and B for AB model!&quot;
 1061   endif
 1062   # transform the R1-matrix to SVAR-style restrictions
 1063   if exists(R1)
 1064     r = rows(R1)
 1065     c = cols(R1)
 1066     if (r != n || c != n)
 1067       funcerr &quot;wrong R1 dimensions&quot;
 1068     endif
 1069     string sBorC = type == 3 ? &quot;B&quot; : &quot;C&quot; # new by Sven 1.0.2
 1070     loop i=1..n -q
 1071       loop j=1..n -q
 1072         scalar rij = R1[i,j]
 1073         if ok(rij) # valid number = restricted element
 1074           SVAR_restrict(&amp;m, sBorC, i, j, rij)
 1075         endif
 1076       endloop
 1077     endloop
 1078   endif
 1079   # still need to consider the A or longrun-C matrix
 1080   # transform the R2-matrix to SVAR-style restrictions
 1081   if exists(R2)
 1082     r2 = rows(R2)
 1083     c2 = cols(R2)
 1084     if (r2 != n || c2 != n)
 1085       funcerr &quot;wrong R2 dimension&quot;
 1086     endif
 1087     string sAorlrC = type == 3 ? &quot;A&quot; : &quot;lrC&quot;
 1088     loop i=1..n -q
 1089       loop j=1..n -q
 1090         scalar rij = R2[i,j]
 1091         if ok(rij) # valid number = restricted element
 1092           SVAR_restrict(&amp;m, sAorlrC, i, j, rij)
 1093         endif
 1094       endloop
 1095     endloop
 1096   endif
 1097 endif
 1098 # do an explicit ID check (sven 1.0.2)
 1099 scalar id_ok = 1
 1100 if checkident
 1101   if (type == 2  &amp;&amp; exists(R2) ) # longrun C
 1102     print &quot;FIXME: not yet implemented for models with long-run restrictions&quot;
 1103   else
 1104     print &quot;Check identification:&quot;
 1105     id_ok = SVAR_ident(&amp;m, 1)	# request verbosity==1 to get messages
 1106   endif
 1107 endif
 1108 # and of course estimate
 1109 if !id_ok
 1110   return m
 1111 else
 1112   SVAR_estimate(&amp;m)
 1113   if b &gt; 0
 1114     SVAR_boot(&amp;m, b, alpha, 0)
 1115   endif
 1116 endif
 1117 return m
 1118 </code>
 1119 </gretl-function>
 1120 <gretl-function name="GUI_plot" type="void" pkg-role="bundle-plot">
 1121  <params count="2">
 1122   <param name="b" type="bundleref"/>
 1123   <param name="ptype" type="int" min="0" max="2" default="0">
 1124 <description>Plot type</description>
 1125 <labels count="3">
 1126 "IRF" "FEVD" "Historical decomposition" </labels>
 1127   </param>
 1128  </params>
 1129 <code>if ptype == 0
 1130   IRFplot(&amp;b, 0, 0, (b.nboot ? 2 : 0)) # all in one go / no key if no bands
 1131 elif ptype == 1
 1132   FEVDplot(&amp;b, 0)
 1133 elif ptype == 2
 1134   HDplot(&amp;b, 0)
 1135 endif
 1136 </code>
 1137 </gretl-function>
 1138 <gretl-function name="FEVDplot" type="void">
 1139  <params count="4">
 1140   <param name="obj" type="bundleref"/>
 1141   <param name="vnum" type="int" default="0"/>
 1142   <param name="keypos" type="int" min="0" max="2" default="1"/>
 1143   <param name="drawix" type="int" min="0" default="0"/>
 1144  </params>
 1145 <code>FEVDsave(&quot;display&quot;, &amp;obj, vnum, keypos, drawix)
 1146 </code>
 1147 </gretl-function>
 1148 <gretl-function name="FEVDsave" type="void">
 1149  <params count="5">
 1150   <param name="outfilename" type="string"/>
 1151   <param name="obj" type="bundleref"/>
 1152   <param name="vnum" type="int" default="0"/>
 1153   <param name="keypos" type="int" min="0" max="2" default="1"/>
 1154   <param name="drawix" type="int" min="0" default="0"/>
 1155  </params>
 1156 <code>scalar n = obj.n
 1157 scalar vfrom vto
 1158 is_vrange = range(vnum, n, &amp;vfrom, &amp;vto)
 1159 matrix Fmat = FEVD(&amp;obj, drawix)
 1160 # prepare title amendment in SR case
 1161 string titletail = &quot;&quot;
 1162 if obj.type == 10
 1163   whichdraw = drawix ? drawix : obj.bestdraw
 1164   titletail = sprintf(&quot; (draw %d)&quot;, whichdraw)
 1165 endif
 1166 loop vnum = vfrom..vto -q
 1167   err = check_bounds(1, vnum, n)
 1168   if err == 2
 1169     printf &quot;Variable number %d out of bounds\n&quot;, vnum
 1170     return
 1171   endif
 1172   string tmpout = (outfilename == &quot;display&quot;) ? &quot;display&quot; : sprintf(&quot;%s_%d&quot;, outfilename, vnum)
 1173   string tmpfile = FEVDgrph(Fmat, vnum, obj.Ynames[vnum], obj.snames, keypos, titletail)
 1174   gnuplot --input=&quot;@tmpfile&quot; --output=&quot;@tmpout&quot;
 1175 endloop
 1176 </code>
 1177 </gretl-function>
 1178 <gretl-function name="HDplot" type="void">
 1179  <params count="3">
 1180   <param name="obj" type="bundleref"/>
 1181   <param name="vnum" type="int" default="0"/>
 1182   <param name="drawix" type="int" min="0" default="0"/>
 1183  </params>
 1184 <code>HDsave(&quot;display&quot;, &amp;obj, vnum, drawix)
 1185 </code>
 1186 </gretl-function>
 1187 <gretl-function name="HDsave" type="void">
 1188  <params count="4">
 1189   <param name="outfilename" type="string"/>
 1190   <param name="obj" type="bundleref"/>
 1191   <param name="vnum" type="int" default="0"/>
 1192   <param name="drawix" type="int" min="0" default="0"/>
 1193  </params>
 1194 <code># more fashionable (=Dynare-like) style
 1195 # Interpret vnum==0 as meaning &quot;all 1..n&quot;
 1196 scalar n = obj.n
 1197 scalar vfrom vto
 1198 is_vrange = range(vnum, n, &amp;vfrom, &amp;vto)
 1199 loop vnum = vfrom..vto -q
 1200   err = check_bounds(1, vnum, n)
 1201   if err == 2
 1202     printf &quot;Variable number %d out of bounds\n&quot;, vnum
 1203     return
 1204   endif
 1205   list tmp = SVAR_hd(&amp;obj, vnum, drawix)
 1206   # prepare extended plot title, only for SR
 1207   whichdraw = drawix ? drawix : obj.bestdraw
 1208   string drawinfo = obj.type == 10 ? sprintf(&quot; (draw %d)&quot;, whichdraw) : &quot;&quot;
 1209   tmp -= tmp[1] # take away deterministic component
 1210   loop i=1..nelem(tmp) --quiet
 1211     string sn = obj.snames[i]
 1212     setinfo tmp[i] --graph-name=&quot;@sn&quot;
 1213   endloop
 1214   series tmpvar = sum(tmp)
 1215   string gnam = sprintf(&quot;%s (stoch. component)&quot;, obj.Ynames[vnum])
 1216   setinfo tmpvar --graph-name=&quot;@gnam&quot;
 1217   tmp = tmp tmpvar
 1218   # (put the line with the target var (stoch.comp) _last_
 1219   # such that the shock colors match those in FEVD
 1220   string tmpout = (outfilename == &quot;display&quot;) ? &quot;display&quot; : sprintf(&quot;%s_%d&quot;, outfilename, vnum)
 1221   plot tmp
 1222     option time-series
 1223     option single-yaxis
 1224     option with-boxes
 1225     option with-lines=tmpvar
 1226     printf &quot;set title \&quot;HD for %s%s\&quot;&quot;, obj.Ynames[vnum], drawinfo
 1227   end plot --output=&quot;@tmpout&quot;
 1228 endloop
 1229 </code>
 1230 </gretl-function>
 1231 <gretl-function name="SVAR_bundle_print" type="void" pkg-role="bundle-print">
 1232  <params count="1">
 1233   <param name="b" type="bundleref"/>
 1234  </params>
 1235 <code># Some specification echoing (sven 1.0.2)
 1236 # (not sure whether this should go into SVAR_est_printout() instead...)
 1237 loop foreach i type n k p --quiet
 1238   scalar $i = b.$i
 1239 endloop
 1240 printf &quot;Model type: %s\n&quot;, modelstring(type)
 1241 strings Ynames = b.Ynames
 1242 print &quot;Endogenous variables:&quot;
 1243 loop i=1..n --quiet
 1244   printf &quot;%s&quot;, Ynames[i]
 1245   if (i==n)
 1246     printf &quot;\n&quot;
 1247   else
 1248     printf &quot;, &quot;
 1249   endif
 1250 endloop
 1251 printf &quot;\n&quot;
 1252 if k&gt;0
 1253   print &quot;Exogenous variables:&quot;
 1254   strings Xnames = b.Xnames
 1255   loop i=1..k --quiet
 1256     printf &quot;%s&quot;, Xnames[i]
 1257     if (i==k)
 1258       printf &quot;\n&quot;
 1259     else
 1260       printf &quot;, &quot;
 1261     endif
 1262   endloop
 1263 else
 1264   print &quot;(No exogenous variables.)&quot;
 1265 endif
 1266 printf &quot;\n&quot;
 1267 printf &quot;Restriction patterns:\n\n&quot;
 1268 if type == 1
 1269   print &quot;Lower-triangular C matrix (Choleski decomposition)\n&quot;
 1270 elif type == 2 # C-model
 1271   if rows(b.Rd1)&gt;0
 1272     printf &quot;Short-run restrictions:\n&quot;
 1273     printf &quot;%3.0f\n&quot;, b.Rd1
 1274   else
 1275     print &quot;(No short-run restrictions)&quot;
 1276   endif
 1277   if rows(b.Rd1l)&gt;0
 1278     print &quot;Long-run restrictions:&quot;
 1279     printf &quot;%3.0f\n&quot;, b.Rd1l
 1280   else
 1281     print &quot;(No long-run restrictions.)&quot;
 1282   endif
 1283 elif type == 3 # AB-model
 1284   if rows(b.Rd0)
 1285     print &quot;Restrictions on A:&quot;
 1286     printf &quot;%3.0f\n&quot;, b.Rd0
 1287   else
 1288     print &quot;(No restrictions on A.)&quot;
 1289   endif
 1290   if rows(b.Rd1)
 1291     print &quot;Restrictions on B:&quot;
 1292     printf &quot;%3.0f\n&quot;, b.Rd1
 1293   else
 1294     print &quot;(No restrictions on B.)&quot;
 1295   endif
 1296 endif
 1297 # only print this if actual estimation took place
 1298 if inbundle(b, &quot;Sigma&quot;)
 1299   printf &quot;Sigma = \n%10.6f&quot;, b.Sigma
 1300   SVAR_est_printout(&amp;b)
 1301 endif
 1302 </code>
 1303 </gretl-function>
 1304 <gretl-function name="SVAR_SRplain" type="void">
 1305  <params count="6">
 1306   <param name="mod" type="bundleref"/>
 1307   <param name="yname" type="string"/>
 1308   <param name="sname" type="string"/>
 1309   <param name="what" type="string"/>
 1310   <param name="length" type="int" min="0" default="0"/>
 1311   <param name="ini" type="int" min="0" default="0"/>
 1312  </params>
 1313 <code># shorthand for sign restrictions in narrow sense
 1314 # (Model type is checked downstream in SVAR_SRfull)
 1315 scalar hi, lo
 1316 if what == &quot;+&quot;
 1317   hi = NA
 1318   lo = 0
 1319 elif what == &quot;-&quot;
 1320   hi = 0
 1321   lo = NA
 1322 else
 1323   funcerr &quot;invalid restriction&quot;
 1324 endif
 1325 SVAR_SRfull(&amp;mod, yname, sname, lo, hi, ini, ini+length)
 1326 </code>
 1327 </gretl-function>
 1328 <gretl-function name="IRF_plotdata" type="bundle">
 1329  <params count="2">
 1330   <param name="bs" type="bundles"/>
 1331   <param name="coveralpha" type="scalar" default="0.9"/>
 1332  </params>
 1333 <code># here we store the data for plotting the irfs in a matrix, conceptually
 1334 # similar to the &quot;bands&quot; matrix we use in SVAR_boot, that is with H rows
 1335 # (=accepted draws, so that descriptive statistics are easy to compute)
 1336 # and ns * nv * horizon columns
 1337 # (But the outputs will be reshaped after the H is aggregated out;
 1338 #  with the horizon in rows and numvariables * numshocks columns.)
 1339 # get dimensions
 1340 H = nelem(bs)    # number of available (accepted) draws
 1341 nvXns = rows(bs[1].Sigma) * cols(bs[1].irfs[1])
 1342 h = nelem(bs[1].irfs)
 1343 matrix A = zeros(H, nvXns * h)
 1344 loop i = 1..H --quiet
 1345   matrices irfs = bs[i].irfs
 1346   matrix sel = seq(1, nvXns * h, h)
 1347   loop t = 1..h --quiet
 1348     # compute the relevant IRFs
 1349     # and store them in the appropriate column
 1350     A[i, sel] = vec(irfs[t])'     # was: b.irfs
 1351     sel += 1
 1352   endloop
 1353 endloop
 1354 matrix m = meanc(A)	# point-wise means
 1355 matrix s = sdc(A)
 1356 # output bundle (H: num of [accepted] replications)
 1357 # (leave out ret.raw = A)
 1358 bundle ret = defbundle(&quot;rep&quot;,H, &quot;coveralpha&quot;,coveralpha, &quot;biascorr&quot;,0)
 1359 matrix ret.irfSRmeans  = mshape(m, h, nvXns)
 1360 matrix ret.irfSRserrs  = mshape(s, h, nvXns)
 1361 # point-wise pseudo confidence bands (ripped from SVAR_boot)
 1362 q_alpha = 0.5 * (1 - coveralpha)
 1363 matrix locb = quantile(A, q_alpha)
 1364 matrix hicb = quantile(A, 1 - q_alpha)
 1365 matrix mdn  = quantile(A, 0.5)
 1366 matrix ret.lo_cb     = mshape(locb, h, nvXns)
 1367 matrix ret.hi_cb     = mshape(hicb, h, nvXns)
 1368 matrix ret.irfSRmeds = mshape(mdn, h, nvXns)
 1369 return ret
 1370 </code>
 1371 </gretl-function>
 1372 <gretl-function name="SVAR_SRexotic" type="void">
 1373  <params count="6">
 1374   <param name="mod" type="bundleref"/>
 1375   <param name="chkstr" type="string"/>
 1376   <param name="shocks" type="strings"/>
 1377   <param name="length" type="int" min="0" default="0"/>
 1378   <param name="ini" type="int" min="0" default="0"/>
 1379   <param name="needs_model" type="bool" default="0"/>
 1380  </params>
 1381 <code>/* Allows to specify more &quot;exotic&quot; restrictions than just an interval
 1382 for one IRF over some horizon. For example, a difference or ratio
 1383 of two IRFs. But the format is necessarily more free-floating here, see the documentation.
 1384 chkstr: a string with a valid numerical evaluation of some function
 1385 of the IRFs at a given horizon tau. The IRF_tau matrix must be
 1386 hardcoded as &quot;M&quot;. Cross-horizon restrictions are impossible.
 1387 length: length of the horizon span over which the restriction should
 1388 hold
 1389 ini: first horizon after shock impact to consider
 1390 (the span is then from ini to ini+length, where impact timing is 0;
 1391 so the default settings just concern the impact effect)
 1392 shocks: strings array with the names of the shocks that enter this
 1393 restriction (it's the user's responsibility to get this right;
 1394 in the futur maybe we can verify or check this)
 1395 needs_model: If yes, this becomes a &quot;super-exotic&quot; restriction, because it needs more input than just the boundaries; for example
 1396 an estimated coefficient (or just the observed variables).
 1397 */
 1398 if mod.type != 10
 1399   funcerr &quot;Wrong model type for set identification restrictions&quot;
 1400 elif $version &lt; 20201
 1401   # (because of feval and errorif() etc., and maybe nested arrays)
 1402   funcerr &quot;Sorry, need at least gretl 2020b for exotic restrictions&quot;
 1403 endif
 1404 if !inbundle(mod, &quot;exoticSR&quot;)
 1405   strings checks = null
 1406   bundle mod.exoticSR = defbundle(&quot;checks&quot;, checks, &quot;spans&quot;, {}, &quot;super&quot;, {}, &quot;eshocks&quot;, {})
 1407 endif
 1408 # store the relevant shock numbers
 1409 loop s = 1..nelem(shocks) -q
 1410   matrix snum = strpos_allin(mod.snames, shocks[s]) # can be several!
 1411   errorif(!nelem(snum), sprintf(&quot;No shock named %s found\n&quot;, shocks[s]))
 1412   mod.exoticSR.eshocks |= snum
 1413 endloop
 1414 finp1 = ini + length + 1 # the indexing convention here matches SRfull, +1
 1415 if mod.horizon &lt;= finp1
 1416   mod.horizon = finp1 + 1 # on the safe side
 1417 endif
 1418 mod.exoticSR.checks += chkstr
 1419 mod.exoticSR.spans = mod.exoticSR.spans | { ini + 1, finp1 }
 1420 mod.exoticSR.super = mod.exoticSR.super | { needs_model }
 1421 </code>
 1422 </gretl-function>
 1423 <gretl-function name="SVAR_spagplot" type="void">
 1424  <params count="3">
 1425   <param name="mod" type="bundle" const="true"/>
 1426   <param name="vname" type="string"/>
 1427   <param name="sname" type="string"/>
 1428  </params>
 1429 <code>if !inbundle(mod, &quot;acc_draws&quot;)
 1430   funcerr &quot;Need collection acc_draws as model bundle member&quot;
 1431 endif
 1432 # check whether it's one of the SR-id'ed shocks
 1433 s_inSR = nelem(strpos_allin(mod.SRid_snames, sname)) ? 1 : 0
 1434 # get shock number
 1435 matrix snum = strpos_allin(mod.snames, sname)
 1436 errmsgshockmatch(snum, sname)
 1437 # get var number
 1438 matrix vnum = strpos_allin(mod.Ynames, vname)
 1439 if !nelem(vnum)
 1440   printf &quot;Regarding variable name %s:\n&quot;, vname
 1441   funcerr &quot;Variable unknown&quot;
 1442 endif
 1443 rep = nelem(mod.acc_draws)  # was bs
 1444 numh = mod.horizon + 1
 1445 matrix IRFs = zeros(numh, rep) ~ seq(0, mod.horizon)'
 1446 loop i = 1..rep --quiet
 1447   IRFs[,i] = drill(mod.acc_draws[i].irfs, vnum[1], snum[1])
 1448 endloop
 1449 plot IRFs
 1450   option with-lines
 1451   #    literal set linetype cycle 4
 1452   literal set linetype 1 lc rgb &quot;#000000&quot; lw 0.1
 1453   literal set linetype 2 lc rgb &quot;#000000&quot; lw 0.1
 1454   literal set linetype 3 lc rgb &quot;#000000&quot; lw 0.1
 1455   literal set linetype 4 lc rgb &quot;#000000&quot; lw 0.1
 1456   literal set linetype 5 lc rgb &quot;#000000&quot; lw 0.1
 1457   literal set linetype 6 lc rgb &quot;#000000&quot; lw 0.1
 1458   literal set linetype 7 lc rgb &quot;#000000&quot; lw 0.1
 1459   literal set linetype 8 lc rgb &quot;#000000&quot; lw 0.1
 1460   literal set xlabel 'lags'
 1461   literal set nokey
 1462   printf &quot;set title '%s -&gt; %s'\n&quot;, sname, vname
 1463 end plot --output=display
 1464 </code>
 1465 </gretl-function>
 1466 <gretl-function name="SVAR_SRfull" type="void">
 1467  <params count="7">
 1468   <param name="mod" type="bundleref"/>
 1469   <param name="yname" type="string"/>
 1470   <param name="sname" type="string"/>
 1471   <param name="lo" type="scalar" default="NA"/>
 1472   <param name="hi" type="scalar" default="NA"/>
 1473   <param name="ini" type="int" min="0" default="0"/>
 1474   <param name="fin" type="int" min="0" default="0"/>
 1475  </params>
 1476 <code># this function will add to the appropriate element of
 1477 # the bundle &quot;mod&quot; a set of sign restrictions relative
 1478 # to one variable and one shock, organised as a row
 1479 # vector with the following elements; input values are:
 1480 #
 1481 # mod      : the VAR model bundle (in pointer form)
 1482 # yname    : a string with the name of the observable to examine
 1483 # sname    : a string with the shock name
 1484 # lo, hi   : bounds (NA for +- infty)
 1485 # ini, fin : IRF interval to check (0-based)
 1486 #
 1487 # for example, if we assume that a monetary policy shock will
 1488 # impact prices with a negative sign on lags from 0 to 4, then
 1489 # the function call should be
 1490 #
 1491 # SVAR_SRfull(mod, &quot;p&quot;, 1, NA, 0, 0, 4)
 1492 #
 1493 # assuming the series for prices is called &quot;p&quot; and that we want
 1494 # our monetary policy shock to come 1st in the ordering (useful for
 1495 # models with sign restrictions on more than one shock)
 1496 # (Switched the yname/sname ordering to match SVAR_SRplain; Sven)
 1497 if mod.type != 10
 1498   funcerr &quot;Wrong model type for set identification restrictions&quot;
 1499 endif
 1500 if !inbundle(mod, &quot;SRest&quot;)
 1501   matrix mod.SRest = {}
 1502 endif
 1503 matrix snum = strpos_allin(mod.snames, sname)
 1504 errmsgshockmatch(snum, sname)
 1505 matrix k = strpos_allin(mod.Ynames, yname)
 1506 if nelem(k) != 1
 1507   printf &quot;Variable %s is not in the model\n&quot;, yname
 1508   print &quot; (or several ones with that name)&quot;
 1509   funcerr &quot;Specify correct and unique variable names&quot;
 1510 endif
 1511 if missing(lo) &amp;&amp; missing(hi)
 1512   funcerr &quot;No bounds specified!\n&quot;
 1513 endif
 1514 lo = ok(lo) ? lo : -$huge
 1515 hi = ok(hi) ? hi :  $huge
 1516 if lo &gt; hi
 1517   printf &quot;Invalid bound specified! [%g : %g]\n&quot;, lo, hi
 1518   funcerr &quot;Provide valid interval&quot;
 1519 endif
 1520 if mod.horizon &lt;= fin
 1521   mod.horizon = fin + 1 # to be on the safe side
 1522 endif
 1523 # Attention: Here the ini and fin inputs are transformed
 1524 # to 1-based matrix indexing form, not keeping the actual
 1525 # horizon meaning
 1526 mod.SRest |= {k, lo, hi, ini+1, fin+1, snum}
 1527 </code>
 1528 </gretl-function>
 1529 <gretl-function name="SVAR_SRirf" type="void">
 1530  <params count="4">
 1531   <param name="mod" type="bundle" const="true">
 1532 <description>the model bundle</description>
 1533   </param>
 1534   <param name="whichvars" type="strings" optional="true" const="true">
 1535 <description>which variables</description>
 1536   </param>
 1537   <param name="whichshocks" type="strings" optional="true" const="true">
 1538 <description>which shocks</description>
 1539   </param>
 1540   <param name="meanormed" type="bool" default="0">
 1541 <description>choose 1 for medians</description>
 1542   </param>
 1543  </params>
 1544 <code># TODO: make the function return the plot code as string
 1545 # TODO: check if it makes sense to use the existing IRFgrph and
 1546 #       IRFsave functions internally
 1547 #       (perhaps the code here is actually better (more modern))
 1548 # Omitting the strings array inputs means: do all.
 1549 if !inbundle(mod, &quot;SRid_snames&quot;)
 1550   funcerr &quot;No SR (sign restriction) specs found&quot;
 1551 elif !inbundle(mod, &quot;SRirfmeans&quot;)
 1552   funcerr &quot;No SR IRF results found - did you run SVAR_SRdraw?&quot;
 1553 endif
 1554 # bundle IrfData = IRF_plotdata(SR_IRFs, coveralpha)
 1555 # Which variables:
 1556 if !exists(whichvars)
 1557   strings whichvars = array(0) # default induces &quot;all&quot;
 1558 endif
 1559 matrix wvarix = names2indices(mod.Ynames, whichvars)
 1560 if !nelem(wvarix)
 1561   funcerr &quot;No matching variable names found&quot;
 1562 endif
 1563 # Which shocks:
 1564 if !exists(whichshocks)
 1565   strings whichshocks = array(0)
 1566 endif
 1567 matrix wshoix = names2indices(mod.SRid_snames, whichshocks)
 1568 if !nelem(wshoix)
 1569   funcerr &quot;No matching shock names found&quot;
 1570 endif
 1571 ## The actual plotting
 1572 loop metaj = 1..nelem(wshoix) -q # for the shocks
 1573   j = wshoix[metaj]
 1574   loop metai = 1..nelem(wvarix) -q # for the target variables
 1575     i = wvarix[metai]
 1576     # pick the right position in the matrix
 1577     k = (j-1) * nelem(mod.Ynames) + i
 1578     matrix center = meanormed ? mod.SRirfmeds : mod.SRirfmeans
 1579     matrix grph = center[,k] ~ seq(0, mod.horizon)' # IrfData.irfSRmeans
 1580     matrix band = mod.SRlo_cb[,k] ~ mod.SRhi_cb[,k] # IrfData.lo_cb, IrfData.hi_cb
 1581     band = 0.5 * band * {1, 1; 1, -1}
 1582     plot grph
 1583       options with-lines fit=none
 1584       options band=band,1 band-style=fill
 1585       printf &quot;set title '%s shock -&gt; %s'&quot;, mod.SRid_snames[j], mod.Ynames[i]
 1586     end plot --output=display
 1587   endloop
 1588 endloop
 1589 # return the gnuplot plotting code here eventually
 1590 </code>
 1591 </gretl-function>
 1592 <gretl-function name="SVAR_SRdraw" type="bundles">
 1593  <params count="5">
 1594   <param name="mod" type="bundleref"/>
 1595   <param name="rep" type="int"/>
 1596   <param name="DO_BAYES" type="bool" default="0"/>
 1597   <param name="coveralpha" type="scalar" min="0" max="1" default="0.9"/>
 1598   <param name="maxiter" type="int" default="10000"/>
 1599  </params>
 1600 <code># This function draws a random rotations of the Cholesky model
 1601 # until &quot;rep&quot; draws satisfying the sign restrictions have come up.
 1602 # If the DO_BAYES flag is on, the VAR parameters (including Sigma)
 1603 # are resampled too; otherwise, they are kept fixed at the ols
 1604 # estimates.
 1605 if !inbundle(mod, &quot;SRest&quot;) &amp;&amp; !inbundle(mod, &quot;exoticSR&quot;)
 1606   funcerr &quot;No set id restrictions found&quot;
 1607 endif
 1608 #
 1609 # Recover VAR parameters from the model bundle
 1610 #
 1611 matrix B = mod.mu | mod.VARpar'
 1612 scalar numh = mod.horizon + 1
 1613 DO_NORM = 1 # used later, when we check if a draw has to be kept or not
 1614 matrix iXX = {}
 1615 if DO_BAYES     # this is only needed for sampling from the posterior
 1616   matrix Data = (mod.X ~ lags(mod.p, mod.Y, 1))[mod.p + 1 : ,]
 1617   iXX = invpd(Data'Data)
 1618   matrix mod.mreg = Data # needed to re-calc the residuals E later
 1619 endif
 1620 if 0
 1621   # print a few parameters so we don't screw up with lags etc.
 1622   safetycheck1(mod.Sigma, iXX, B)
 1623 endif
 1624 #
 1625 # prepare the output array and other auxiliary stuff
 1626 #
 1627 bundles drawn = array(0)
 1628 good = 0
 1629 iter = 1
 1630 # which shocks are to be identified via sign restrictions?
 1631 matrix shocks = values(mod.SRest[,6])
 1632 matrix allshocks = shocks
 1633 # do we have exotic restrictions?
 1634 n_exotic = get_n_exotic(mod)
 1635 # add the ones from exotic restrictions
 1636 matrix allshocks = n_exotic ? values(shocks | mod.exoticSR.eshocks) : shocks
 1637 scalar nSetIdShocks = rows(shocks)
 1638 # Now extract the corresponding shock names
 1639 strings mod.SRid_snames = getstrings_byindex(mod.snames, allshocks)
 1640 # this array contains as many elements as set-identified shocks
 1641 # basically, it's a per-shock split of the SRest bundle element
 1642 matrices Ais = prepAis(mod, shocks)
 1643 # responses to be cumulated?
 1644 matrix to_cum = mod.ncumul ? mod.cumul : {}
 1645 #
 1646 # the main loop
 1647 #
 1648 bundle b = defbundle(&quot;B&quot;,B, &quot;Sigma&quot;,mod.Sigma, &quot;exoterms&quot;,rows(mod.mu), &quot;df&quot;, mod.T - rows(B), &quot;T&quot;, mod.T) # df was: mod.T - nelem(B)
 1649 # for the future: set a matrix for keeping some shocks from being rotated
 1650 # eventually, this could probably be inferred from the model bundle
 1651 rot_only = {}
 1652 loop while good&lt;rep &amp;&amp; iter&lt;=maxiter --quiet
 1653   # first just using non-exotic restrictions
 1654   # generate the random rotation matrix and, optionally, resample
 1655   rot_redraw(&amp;b, mod, DO_BAYES, B, iXX, rot_only) 	# b.rot, b.B, b.Sigma
 1656   matrices irfs = gen_SRirfs(b, numh, to_cum)
 1657   matrix id_matrix = get_id_mat(mod, shocks, irfs, Ais, DO_NORM)
 1658   chk = DO_NORM ? (rows(id_matrix) == rows(shocks)) : 0 # this check 1st!
 1659   chk = chk &amp;&amp; check_id(id_matrix)
 1660   if chk
 1661     if mod.n == nSetIdShocks
 1662       # complete system
 1663       b.Q = b.rot * id_matrix'
 1664     endif
 1665     # perform the &quot;reflection-permutation&quot; step
 1666     loop i = 1..nelem(irfs) --quiet
 1667       irfs[i] = irfs[i] * id_matrix'
 1668     endloop
 1669     # --- now check for exotic restrictions if necessary -------
 1670     if n_exotic
 1671       loop i = 1..n_exotic --quiet
 1672         out = exotic_inner_check(&amp;chk, i, mod, irfs)
 1673         if out
 1674           break
 1675         endif
 1676       endloop
 1677     endif
 1678     if chk
 1679       good++
 1680       b.irfs = irfs
 1681       # copy the results of accepted draw to output bundles array
 1682       drawn += b
 1683     endif
 1684   endif
 1685   factor = (DO_BAYES ? 10 : 50) * xmax(1,floor(sqrt(maxiter/1000)))
 1686   if iter % factor == 0
 1687     printf &quot;draw %5d done (good so far = %5d)\n&quot;, iter, good
 1688     flush
 1689   endif
 1690   iter++
 1691 endloop
 1692 if good &lt; rep
 1693   printf &quot;\n\nWARNING: COULDN'T ACHIEVE THE DESIRED NUMBER OF DRAWS&quot;
 1694   printf &quot; (only %d out of %d)\n&quot;, good, rep
 1695 else
 1696   printf &quot;draw %5d done (good = %5d)\n&quot;, iter-1, good
 1697 endif
 1698 printf &quot; (acceptance ratio: %g)\n&quot;, good/(iter-1)
 1699 # copy info into model
 1700 mod.SRiter = iter - 1
 1701 mod.SRacc  = good
 1702 mod.DO_BAYES = DO_BAYES
 1703 mod.SRcoveralpha = coveralpha
 1704 # integrate some summary results into the model bundle
 1705 if good
 1706   bundle IrfData = IRF_plotdata(drawn, coveralpha)
 1707   matrix mod.SRirfmeans = IrfData.irfSRmeans
 1708   matrix mod.SRirfserrs = IrfData.irfSRserrs
 1709   matrix mod.SRlo_cb    = IrfData.lo_cb
 1710   matrix mod.SRhi_cb    = IrfData.hi_cb
 1711   matrix mod.SRirfmeds  = IrfData.irfSRmeds # medians
 1712   # also store everything in the main bundle if wanted
 1713   if inbundle(mod, storeSRirfs) &amp;&amp; mod.storeSRirfs
 1714     # (idiom possible because for SR we require recent gretl)
 1715     mod.acc_draws = drawn
 1716   endif
 1717 else
 1718   print &quot;Warning: No accepted draws, plotting impossible&quot;
 1719 endif
 1720 return drawn	# may be 0-element
 1721 </code>
 1722 </gretl-function>
 1723 <gretl-function name="SRgetbest" type="scalar">
 1724  <params count="7">
 1725   <param name="mod" type="bundleref"/>
 1726   <param name="vname" type="string"/>
 1727   <param name="sname" type="string"/>
 1728   <param name="ini" type="int" min="0" default="0"/>
 1729   <param name="length" type="int" min="0" default="0"/>
 1730   <param name="disttarg" type="int" min="0" max="1" default="0">
 1731 <description>measure to target</description>
 1732 <labels count="2">
 1733 "mean" "median" </labels>
 1734   </param>
 1735   <param name="loss" type="string" optional="true"/>
 1736  </params>
 1737 <code># Function to find the set-id accepted draw which best fits a
 1738 # a certain IRF center
 1739 # vname: e.g. &quot;Y&quot;
 1740 # sname: e.g. &quot;supply&quot; (naming must exist in the bundle)
 1741 # ini: horizon range start (0 = impact)
 1742 # length: horizon end is ini + length
 1743 # disttarg: which part of the per-horizon distribution should be
 1744 # targeted;
 1745 # loss: how to evaluate deviations from the target;
 1746 #      possible values: &quot;abs&quot;, &quot;quad&quot; (default), ...
 1747 ## implement defaults and error checks:
 1748 if !exists(loss)
 1749   string loss = &quot;quad&quot;
 1750 elif tolower(loss) != &quot;abs&quot; &amp;&amp; tolower(loss) != &quot;quad&quot;
 1751   print &quot;Warning, unrecognized loss, using 'quad'&quot;
 1752   loss = &quot;quad&quot;
 1753 else
 1754   loss = tolower(loss)
 1755 endif
 1756 if ini + length &gt; mod.horizon
 1757   funcerr &quot;ini + length too large (increase horizon)&quot;
 1758 endif
 1759 if !inbundle(mod, &quot;SRacc&quot;)
 1760   funcerr &quot;No info on accepted draws found&quot;
 1761 elif !inbundle(mod, &quot;acc_draws&quot;)
 1762   funcerr &quot;Model must be set up to store accepted draws in its bundle&quot;
 1763 endif
 1764 startix = ini + 1	# offset for 1-based
 1765 finix = startix + length
 1766 matrix center = !disttarg ? mod.SRirfmeans : mod.SRirfmeds
 1767 # pick the right position in the matrix:
 1768 # (taking into account partial identification, i.e.
 1769 #  less id'ed shocks than variables)
 1770 scalar v_ix = names2indices(mod.Ynames, defarray(vname))
 1771 scalar s_ix = names2indices(mod.SRid_snames, defarray(sname))
 1772 k = (s_ix - 1) * nelem(mod.Ynames) + v_ix
 1773 ## debug:
 1774 /*
 1775 print s_ix
 1776 print v_ix
 1777 print k
 1778 eval cols(center)
 1779 eval nelem(mod.SRid_snames)
 1780 eval nelem(mod.snames)
 1781 print mod
 1782 matrices mm = mod.acc_draws[1].irfs
 1783 matrix mtemp = mm[1]
 1784 eval rows(mtemp)
 1785 eval cols(mtemp)
 1786 */
 1787 center = center[,k]
 1788 best = $huge
 1789 bestdraw = 0
 1790 loop i = 1..mod.SRacc -q
 1791   matrix IRFvec = drill(mod.acc_draws[i].irfs, v_ix, s_ix)
 1792   matrix distance = (IRFvec - center)[startix : finix]
 1793   if loss == &quot;abs&quot;
 1794     scalar objective = sum(abs(distance))
 1795   elif loss == &quot;quad&quot;
 1796     scalar objective = sum(distance.^2)
 1797   endif
 1798   if objective &lt; best
 1799     best = objective
 1800     bestdraw = i
 1801   endif
 1802 endloop
 1803 if !bestdraw
 1804   funcerr &quot;shouldn't happen&quot;
 1805 endif
 1806 # store for further use
 1807 mod.bestdraw = bestdraw
 1808 return bestdraw
 1809 </code>
 1810 </gretl-function>
 1811 <gretl-function name="SVAR_getshock" type="series">
 1812  <params count="3">
 1813   <param name="mod" type="bundleref"/>
 1814   <param name="sname" type="string" optional="true"/>
 1815   <param name="drawix" type="int" min="0" default="0"/>
 1816  </params>
 1817 <code># This is a wrapper to provide a nicer interface, using the
 1818 # name of the shock instead of the number.
 1819 # Default (as in GetShock) is to use the first shock.
 1820 s_ix = !exists(sname) ? 1 : strpos_allin(mod.snames, sname)
 1821 return GetShock(&amp;mod, s_ix, drawix)
 1822 </code>
 1823 </gretl-function>
 1824 <gretl-function name="SVAR_HD" type="list">
 1825  <params count="3">
 1826   <param name="mod" type="bundleref"/>
 1827   <param name="vname" type="string" optional="true"/>
 1828   <param name="drawix" type="int" min="0" default="0"/>
 1829  </params>
 1830 <code># wrapper around SVAR_hd to use a string interface for the
 1831 # variable
 1832 v_ix = !exists(vname) ? 1 : strpos_allin(mod.Ynames, vname)
 1833 return SVAR_hd(&amp;mod, v_ix, drawix)
 1834 </code>
 1835 </gretl-function>
 1836 <gretl-function name="getstrings_byindex" type="strings" private="1">
 1837  <params count="2">
 1838   <param name="S" type="strings" const="true"/>
 1839   <param name="indices" type="matrix" const="true"/>
 1840  </params>
 1841 <code># collects the (in general non-contiguous) sub-array
 1842 # indexed by indices (imitate fancy indexing)
 1843 strings out = null
 1844 matrix m = vec(indices)
 1845 loop i = 1..nelem(m) -q
 1846   out += S[m[i]]
 1847 endloop
 1848 return out
 1849 </code>
 1850 </gretl-function>
 1851 <gretl-function name="drawnormwis" type="matrices" private="1">
 1852  <params count="4">
 1853   <param name="iXX" type="matrix" const="true"/>
 1854   <param name="B" type="matrix" const="true"/>
 1855   <param name="Sigma" type="matrix" const="true"/>
 1856   <param name="T" type="int" min="1"/>
 1857  </params>
 1858 <code># Draw from a standard Normal-(inverse)-Wishart prior
 1859 # for a multi-equation regression model.
 1860 #
 1861 # iXX:  K x K matrix (X'X)^{-1} (same for all equations!)
 1862 #       (in a VAR context gretl provides this as $xtxinv)
 1863 # B:    matrix of purely data-based estimates (max-lik)
 1864 # Sigma: cross-equation covariance matrix of innovations
 1865 # T:    number of observations
 1866 K = rows(B)
 1867 N = cols(B) # how many equations
 1868 # some checks
 1869 if K != rows(iXX)
 1870   funcerr &quot;Coeff and data matrix dims don't match&quot;
 1871 elif N != rows(Sigma)
 1872   funcerr &quot;Coeff and Cov matrix dims don't match&quot;
 1873 endif
 1874 matrix Sigma_draw = iwishart(Sigma*T, T)
 1875 matrix V = Sigma_draw ** iXX
 1876 matrix C = cholesky(V)
 1877 matrix B_draw = vec(B)
 1878 B_draw += C * mnormal(K * N, 1)
 1879 B_draw = mshape(B_draw, K, N)
 1880 return defarray(B_draw, Sigma_draw)
 1881 </code>
 1882 </gretl-function>
 1883 <gretl-function name="strpos_allin" type="matrix" private="1">
 1884  <params count="2">
 1885   <param name="S" type="strings"/>
 1886   <param name="search_pattern" type="string"/>
 1887  </params>
 1888 <code># Function to determine the numerical position(s) of a
 1889 # string in an array
 1890 matrix ret = {}
 1891 loop i=1..nelem(S) -q
 1892   if S[i] == search_pattern
 1893     ret |= i
 1894   endif
 1895 endloop
 1896 return ret
 1897 </code>
 1898 </gretl-function>
 1899 <gretl-function name="basename" type="strings" private="1">
 1900  <params count="1">
 1901   <param name="fn" type="string"/>
 1902  </params>
 1903 <code>string base = regsub(fn, &quot;(.*)\.([^\.]*)&quot;, &quot;\1&quot;)
 1904 string ext = fn + (strlen(base) + 1)
 1905 return defarray(base, ext)
 1906 </code>
 1907 </gretl-function>
 1908 <gretl-function name="drill" type="matrix" private="1">
 1909  <params count="3">
 1910   <param name="x" type="matrices" const="true"/>
 1911   <param name="rowspec" type="matrix" optional="true"/>
 1912   <param name="colspec" type="matrix" optional="true"/>
 1913  </params>
 1914 <code># This function &quot;drills through&quot; a matrix array and returns a matrix;
 1915 # for example, drill(x, 2, 3) returns a vector with the [2,3] elements
 1916 # of all matrices in the x array. &quot;0&quot; means &quot;all&quot;.
 1917 #
 1918 # NOTA BENE: all matrices must be the same size
 1919 matrix ret = {}
 1920 n = nelem(x)
 1921 if n == 0
 1922   return ret
 1923 endif
 1924 ### check sizes
 1925 nr = rows(x[1])
 1926 nc = cols(x[1])
 1927 same_dim = 1
 1928 loop i = 2 .. n --quiet
 1929   same_dim = same_dim &amp;&amp; (rows(x[i]) == nr) &amp;&amp; (cols(x[i]) == nc)
 1930   if !same_dim
 1931     printf &quot;Error: Inconsistent dimensions (not all matrices are the same size)\n&quot;
 1932     return ret
 1933   endif
 1934 endloop
 1935 ### process specs
 1936 if !exists(rowspec)
 1937   matrix rs = seq(1, nr)'
 1938 else
 1939   if rowspec[1] == 0
 1940     matrix rs = seq(1, nr)'
 1941   else
 1942     matrix rs = vec(rowspec) # force to column
 1943   endif
 1944 endif
 1945 if !exists(colspec)
 1946   matrix cs = seq(1, nc)'
 1947 else
 1948   if colspec[1] == 0
 1949     matrix cs = seq(1, nc)'
 1950   else
 1951     matrix cs = vec(colspec) # force to column
 1952   endif
 1953 endif
 1954 ### check for multiple or illegal specs
 1955 scalar nrspec = rows(rs)
 1956 scalar ncspec = rows(cs)
 1957 if xmin(nrspec, ncspec) &gt; 1
 1958   printf &quot;Error: you can’t have multiple row and column specs\n&quot;
 1959   return ret
 1960 endif
 1961 if minc(rs|cs) &lt; 0
 1962   printf &quot;Error: negative spec not allowed\n&quot;
 1963   return ret
 1964 endif
 1965 if maxc(rs) &gt; nr
 1966   printf &quot;Error: incorrect row spec (matrices have %d rows, but %d wanted)\n&quot;, nr, maxc(rs)
 1967   return ret
 1968 endif
 1969 if maxc(cs) &gt; nc
 1970   printf &quot;Error: incorrect col spec (matrices have %d columns, but %d wanted)\n&quot;, nc, maxc(cs)
 1971   return ret
 1972 endif
 1973 ### do the actual drilling
 1974 if nrspec == 1
 1975   ret = flatten(x)[rs,]
 1976   ret = transp(mshape(ret, nc, n))
 1977   ret = ret[,cs]
 1978 elif ncspec == 1
 1979   ret = flatten(x,1)[,cs]
 1980   ret = mshape(ret, nr, n)
 1981   ret = ret[rs,]
 1982 endif
 1983 return ret
 1984 </code>
 1985 </gretl-function>
 1986 <gretl-function name="names2indices" type="matrix" private="1">
 1987  <params count="2">
 1988   <param name="fromwhich" type="strings"/>
 1989   <param name="names" type="strings"/>
 1990  </params>
 1991 <code>matrix which = {}
 1992 if !nelem(names) # default, all
 1993   which = seq(1, nelem(fromwhich))
 1994 else
 1995   loop n = 1..nelem(names) -q
 1996     which |= strpos_allin(fromwhich, names[n])
 1997   endloop
 1998   which = values(which)
 1999 endif
 2000 return which # may still be empty
 2001 </code>
 2002 </gretl-function>
 2003 <gretl-function name="range" type="scalar" private="1">
 2004  <params count="4">
 2005   <param name="a" type="int"/>
 2006   <param name="n" type="int"/>
 2007   <param name="n0" type="scalarref"/>
 2008   <param name="n1" type="scalarref"/>
 2009  </params>
 2010 <code># this is used in the plotting function to have
 2011 # 0 as a synonym for &quot;all&quot;
 2012 #
 2013 # The flipping of the shock by passing a negative number
 2014 # is then not possible to specify; users have to do it explicitly
 2015 # in their own loop then.
 2016 #
 2017 # (n0 and n1 are integers, too, but gretl doesn't accept
 2018 #   pointerized integers)
 2019 if a != 0
 2020   ret = 0
 2021   n0 = a
 2022   n1 = a
 2023 else
 2024   ret = 1
 2025   n0 = 1
 2026   n1 = n
 2027 endif
 2028 return ret
 2029 </code>
 2030 </gretl-function>
 2031 <gretl-function name="max_eval" type="scalar" private="1">
 2032  <params count="1">
 2033   <param name="A" type="matrix" const="true"/>
 2034  </params>
 2035 <code>n = rows(A)
 2036 p = cols(A) / n
 2037 matrix compan = (p==1) ? A : A | (I(n*(p-1)) ~ zeros(n*(p-1), n))
 2038 # TODO: eventually replace this with &quot;eigen&quot; (but needs 2020b or so):
 2039 matrix lambda = eigengen(compan) # deprecated
 2040 scalar maxmod = maxc(sumr(lambda.^2))
 2041 return maxmod
 2042 </code>
 2043 </gretl-function>
 2044 <gretl-function name="determ" type="list" private="1">
 2045  <params count="4">
 2046   <param name="X" type="list" const="true"/>
 2047   <param name="cnst" type="bool"/>
 2048   <param name="trnd" type="bool"/>
 2049   <param name="seas" type="bool"/>
 2050  </params>
 2051 <code>list ret = cnst ? const : null
 2052 if trnd
 2053   ret += time
 2054 endif
 2055 if seas
 2056   ret += seasonals(0, 1) # centered
 2057 endif
 2058 # stick together deterministics and other exog.
 2059 ret = ret || X
 2060 return ret
 2061 </code>
 2062 </gretl-function>
 2063 <gretl-function name="vecm_det" type="matrix" private="1">
 2064  <params count="2">
 2065   <param name="T" type="int"/>
 2066   <param name="dcase" type="int"/>
 2067  </params>
 2068 <code># build the deterministic matrix for the VECM; if alpha is
 2069 # empty, it will be estimated via ols later
 2070 # deterministics
 2071 # note that in the &quot;even cases&quot; (restr. const or trend)
 2072 # the restricted term is _not_ included in x, since its
 2073 # parameter will be recovered later via alpha*beta
 2074 matrix mreg = (dcase &lt; 3) ? {} : ones(T,1)
 2075 if dcase == 5
 2076   matrix mreg ~= seq(1,T)'
 2077 endif
 2078 return mreg
 2079 </code>
 2080 </gretl-function>
 2081 <gretl-function name="N2_ify" type="matrix" private="1">
 2082  <params count="1">
 2083   <param name="A" type="matrix" const="true"/>
 2084  </params>
 2085 <code>n2 = rows(A)
 2086 n = int(sqrt(n2))
 2087 k = cols(A)
 2088 matrix e = vec(transp(mshape(seq(1,n2),n,n)))
 2089 return A + A[e,1:k]
 2090 </code>
 2091 </gretl-function>
 2092 <gretl-function name="has_unit_diag" type="scalar" private="1">
 2093  <params count="1">
 2094   <param name="Rd" type="matrix" const="true"/>
 2095  </params>
 2096 <code>ret = 0
 2097 n2 = cols(Rd) - 1
 2098 n = sqrt(n2)
 2099 matrix test = transp(seq(0,n2-1)) .= seq(0,n2-1,n+1)
 2100 test |= ones(1,n)
 2101 matrix e
 2102 mols(test, Rd', &amp;e)
 2103 return maxr(sumc(e.*e)) &lt; 1.0e-12 # ret
 2104 </code>
 2105 </gretl-function>
 2106 <gretl-function name="has_free_diag" type="scalar" private="1">
 2107  <params count="1">
 2108   <param name="Rd" type="matrix" const="true"/>
 2109  </params>
 2110 <code>n2 = cols(Rd) - 1
 2111 n = sqrt(n2)
 2112 matrix e = 1 + seq(0,n2-1,n+1)
 2113 matrix test = Rd[,e]
 2114 return ( sumc(abs(test)) &lt; 1.0e-12 ) # ret
 2115 </code>
 2116 </gretl-function>
 2117 <gretl-function name="mat_exp" type="matrix" private="1">
 2118  <params count="2">
 2119   <param name="theta" type="matrix" const="true"/>
 2120   <param name="Ss" type="matrix" const="true"/>
 2121  </params>
 2122 <code># we don't check for conformability, but
 2123 # cols(Ss) should be equal to rows(theta)+1
 2124 n2 = rows(Ss)
 2125 n = int(sqrt(n2))
 2126 k = cols(Ss) - 1
 2127 matrix C = (k&gt;0) ? ( Ss[,1:k]*theta + Ss[,k+1] ) : Ss
 2128 return mshape(C,n,n)
 2129 </code>
 2130 </gretl-function>
 2131 <gretl-function name="unscramble_dat" type="void" private="1">
 2132  <params count="3">
 2133   <param name="dat" type="matrixref"/>
 2134   <param name="Sigma" type="matrixref"/>
 2135   <param name="Ss" type="matrixref"/>
 2136  </params>
 2137 <code>n2 = rows(dat)
 2138 n = int(sqrt(n2))
 2139 k = cols(dat)
 2140 Sigma = mshape(dat[,1],n,n)
 2141 Ss = dat[,2:k]
 2142 </code>
 2143 </gretl-function>
 2144 <gretl-function name="maybe_flip_columns" type="void" private="1">
 2145  <params count="2">
 2146   <param name="C" type="matrix" const="true"/>
 2147   <param name="X" type="matrixref"/>
 2148  </params>
 2149 <code>/*
 2150 the objective here is to make X as similar as possible to C
 2151 by flipping the sign of the columns. Used for bootstrapping, to have IRFs with comparable signs.
 2152 */
 2153 n = rows(C)
 2154 matrix sel = seq(1,n)
 2155 matrix plus = sumc((C + X).^2)
 2156 matrix minus = sumc((C - X).^2)
 2157 matrix flip = plus .&lt; minus
 2158 if sumr(flip) &gt; 0
 2159   sel = selifc(sel, flip)
 2160   X[,sel] = -X[,sel]
 2161 endif
 2162 </code>
 2163 </gretl-function>
 2164 <gretl-function name="printStrMat" type="void" private="1">
 2165  <params count="3">
 2166   <param name="X" type="matrix" const="true"/>
 2167   <param name="V" type="matrix"/>
 2168   <param name="name" type="string"/>
 2169  </params>
 2170 <code>n = rows(X)
 2171 matrix x = vec(X)
 2172 matrix cfse = vec(X)
 2173 matrix se = sqrt(diag(V))
 2174 matrix numzero = selifr(seq(1,rows(se))', (se .&lt; 1.0e-15))
 2175 if rows(numzero) &gt; 1
 2176   se[numzero] = 0.0
 2177 endif
 2178 cfse ~= se
 2179 string parnames = &quot;&quot;
 2180 loop j=1..n --quiet
 2181   loop i=1..n --quiet
 2182     # sprintf tmpstr &quot;%s[%2d;%2d]&quot;, name, i, j
 2183     # parnames += tmpstr
 2184     parnames += sprintf(&quot;%s[%2d;%2d]&quot;, name, i, j)
 2185     if (j&lt;n) || (i&lt;n)
 2186       parnames += &quot;,&quot;
 2187     endif
 2188   endloop
 2189 endloop
 2190 modprint cfse parnames
 2191 </code>
 2192 </gretl-function>
 2193 <gretl-function name="add_and_smash" type="scalar" private="1">
 2194  <params count="2">
 2195   <param name="A" type="matrixref"/>
 2196   <param name="Psi" type="matrix" const="true"/>
 2197  </params>
 2198 <code># Note that Psi here seems to correspond to -\Psi
 2199 # in Kilian (1998)! That's why it's added.
 2200 matrix Ab = A + Psi
 2201 # now check stationarity
 2202 scalar maxmod = max_eval(Ab)
 2203 h = 0.99
 2204 H = 1
 2205 maxiter = 1000
 2206 iter = 0
 2207 loop while (maxmod &gt; 0.9999) &amp;&amp; (iter &lt; maxiter) --quiet
 2208   iter++
 2209   H *= h
 2210   Ab = A + H .* Psi
 2211   maxmod = max_eval(Ab)
 2212 endloop
 2213 A = Ab
 2214 return (iter &gt;= maxiter) ? NA : H
 2215 </code>
 2216 </gretl-function>
 2217 <gretl-function name="transshockcheck" type="scalar" private="1">
 2218  <params count="1">
 2219   <param name="b" type="bundle" const="true"/>
 2220  </params>
 2221 <code># only for SVEC case
 2222 # check the last r columns of C (trans shocks),
 2223 # shouldn't have more than r-1 zeros in each of them
 2224 # according to Luetkepohl (2008, EL)
 2225 # (for unrestricted alpha)
 2226 # (passing this check doesn't mean everything is OK)
 2227 if b.type != 4
 2228   print &quot;Transitory shocks only for SVEC models&quot;
 2229   return 1
 2230 elif nelem(b.Rd1)	# else there are no short-run restr.
 2231   # check for weakly exogenous variables
 2232   weakexo = 0
 2233   if inbundle(b, &quot;jalpha&quot;)
 2234     if min(sumr(abs(b.jalpha))) &lt; 1e-8
 2235       weakexo = 1
 2236     endif
 2237   endif
 2238   if !weakexo
 2239     loop i = (b.n - b.crank)..(b.n - 1) -q
 2240       matrix part = b.Rd1[, i * b.n + 1: (i+1) * b.n]
 2241       # check for exclusion restrictions (RHS zero)
 2242       matrix mc = (sumr(part) .= 1) &amp;&amp; (b.Rd1[,cols(b.Rd1)] .= 0)
 2243       if sum(mc) &gt; b.crank - 1
 2244         printf &quot;Concerning shock %d\n&quot;, i+1
 2245         funcerr &quot;Too many zero restrictions&quot;
 2246       endif
 2247     endloop
 2248   endif
 2249 else
 2250   print &quot;No short-run restrictions found (Rd1)&quot;
 2251 endif
 2252 return 0
 2253 </code>
 2254 </gretl-function>
 2255 <gretl-function name="CheckNormalizeRd" type="scalar" private="1">
 2256  <params count="2">
 2257   <param name="R" type="matrixref"/>
 2258   <param name="d" type="matrixref"/>
 2259  </params>
 2260 <code>/*
 2261 Checks that
 2262 (1) the constraints are consistent
 2263 (2) the constraints are non-contradictory (??? non-redundant ? - Sven)
 2264 if (1) fails, an error message is printed and R and d are replaced by
 2265   empty matrices; if (2) fails, redundant rows in R and d are dropped.
 2266   */
 2267   p = rows(R)
 2268   r = rank(R)
 2269   ret = 0
 2270   if r &lt; p
 2271     matrix Rd = R ~ d
 2272     if r &lt; rank(Rd) # contradictory
 2273       R = {}
 2274       d = {}
 2275       ret = 1
 2276     else # redundant
 2277       matrix RR
 2278       matrix QQ = qrdecomp(Rd', &amp;RR)
 2279       matrix e = abs(diag(RR)) .&gt; $macheps
 2280       QQ = selifc(QQ, e')
 2281       RR = selifr(selifc(RR, e'), e)
 2282       Rd = QQ * RR
 2283       R = Rd[1:rows(Rd)-1,]'
 2284       d = Rd[rows(Rd),]'
 2285       ret = 2
 2286     endif
 2287   endif
 2288   return ret
 2289 </code>
 2290 </gretl-function>
 2291 <gretl-function name="add_constraint" type="scalar" private="1">
 2292  <params count="5">
 2293   <param name="Rd" type="matrixref"/>
 2294   <param name="n" type="int"/>
 2295   <param name="i" type="int"/>
 2296   <param name="j" type="int"/>
 2297   <param name="d" type="scalar"/>
 2298  </params>
 2299 <code>err = (i&gt;n) || (j&gt;n) || (i&lt;1) || (j&lt;1)
 2300 if !err
 2301   n2 = n*n
 2302   matrix tmp = zeros(1,n2 + 1)
 2303   k = n*(j-1) + i
 2304   tmp[1,k] = 1
 2305   tmp[1,n2+1] = d
 2306   # check for consistency/redundancy
 2307   if rows(Rd) &gt; 0
 2308     matrix newR = Rd[,1:n2] | tmp[1:n2]
 2309     matrix newd = Rd[,n2+1] | d
 2310   else
 2311     matrix newR = tmp[1:n2]
 2312     matrix newd = d
 2313   endif
 2314   err2 = CheckNormalizeRd(&amp;newR, &amp;newd)
 2315   if err2 == 0
 2316     Rd = Rd | tmp
 2317   elif err2 == 1
 2318     printf &quot;The restriction conflicts with previous ones and was ignored\n&quot;
 2319   elif err2 == 2
 2320     printf &quot;The restriction is redundant and was ignored\n&quot;
 2321   endif
 2322 endif
 2323 return err ? -err : 10*err2 # -1 for bad input, 10 or 20 upstream error
 2324 </code>
 2325 </gretl-function>
 2326 <gretl-function name="cholRd" type="matrix" private="1">
 2327  <params count="1">
 2328   <param name="n" type="int"/>
 2329  </params>
 2330 <code>n2 = n*n
 2331 k = 1
 2332 matrix ret = {}
 2333 loop i=1..n -q
 2334   loop j=1..n -q
 2335     matrix tmp = zeros(1,n2+1)
 2336     if i &gt; j
 2337       tmp[k] = 1
 2338       ret = ret | tmp
 2339     endif
 2340     k++
 2341   endloop
 2342 endloop
 2343 return ret
 2344 </code>
 2345 </gretl-function>
 2346 <gretl-function name="diag_Rd" type="matrix" private="1">
 2347  <params count="2">
 2348   <param name="n" type="int"/>
 2349   <param name="x" type="scalar"/>
 2350  </params>
 2351 <code>return selifr(I(n*n), vec(I(n))) ~ (x*ones(n,1))
 2352 </code>
 2353 </gretl-function>
 2354 <gretl-function name="free_diag_Rd" type="matrix" private="1">
 2355  <params count="1">
 2356   <param name="n" type="int"/>
 2357  </params>
 2358 <code>n2 = n*n
 2359 return selifr(I(n2)~zeros(n2,1), !vec(I(n)))
 2360 </code>
 2361 </gretl-function>
 2362 <gretl-function name="imp2exp" type="matrix" private="1">
 2363  <params count="1">
 2364   <param name="Rd" type="matrix"/>
 2365  </params>
 2366 <code>/*
 2367 Given the constraints in implicit form, returns the matrix [ S | s ]
 2368 of the constraints in explicit form
 2369 */
 2370 if !nelem(Rd)
 2371   funcerr &quot;No restrictions given.&quot;
 2372 endif
 2373 p = cols(Rd)
 2374 matrix R = Rd[, 1:(p-1)]
 2375 matrix d = Rd[,p]
 2376 catch matrix RRi = invpd(R*R')
 2377 err = $error
 2378 if err
 2379   print &quot;Processing of restrictions failed aborting.&quot;
 2380   print &quot;Perhaps redundant or incompatible restrictions?&quot;
 2381   printf &quot;%s\n&quot;, errmsg(err)
 2382   funcerr &quot;aborting&quot;
 2383 endif
 2384 matrix s = R'RRi*d
 2385 matrix S = nullspace(R | s')
 2386 return S ~ s
 2387 </code>
 2388 </gretl-function>
 2389 <gretl-function name="check_const" type="void" private="1">
 2390  <params count="2">
 2391   <param name="K" type="matrix"/>
 2392   <param name="fullRd" type="matrix"/>
 2393  </params>
 2394 <code>k = cols(fullRd) - 1
 2395 matrix R = fullRd[,1:k]
 2396 matrix d = fullRd[,k+1]
 2397 printf &quot;Constraint matrix:\n&quot;
 2398 print R
 2399 matrix tmp = R*vec(K)
 2400 print tmp
 2401 printf &quot;Should be:\n&quot;
 2402 print d
 2403 </code>
 2404 </gretl-function>
 2405 <gretl-function name="lrConstr" type="matrix" private="1">
 2406  <params count="2">
 2407   <param name="C1" type="matrix"/>
 2408   <param name="lrRd" type="matrix"/>
 2409  </params>
 2410 <code>if rows(lrRd) == 0
 2411   return {} # nothing to do
 2412 else
 2413   n = rows(C1)
 2414   k = cols(lrRd) - 1
 2415   matrix d = lrRd[,k+1]
 2416   return ( lrRd[,1:k] * (I(n) ** C1) ) ~ d
 2417 endif
 2418 return ret
 2419 </code>
 2420 </gretl-function>
 2421 <gretl-function name="get_full_Rd" type="matrix" private="1">
 2422  <params count="2">
 2423   <param name="obj" type="bundleref"/>
 2424   <param name="verbosity" type="int" default="0"/>
 2425  </params>
 2426 <code>scalar type = obj.type
 2427 scalar n = obj.n
 2428 matrix fullRd = obj.Rd1 # grab sr-restrictions first
 2429 # (for type 1 or type 2 with short-run only: nothing else to do)
 2430 lr_constr = rows(obj.Rd1l) ? 1 : 0	# further long-run constraints?
 2431 if type == 3
 2432   funcerr &quot;This func not usable for AB models (yet)&quot;
 2433   # because long-run restrictions aren't implemented for them
 2434 elif (type == 4) || lr_constr || obj.calc_lr
 2435   # generic C-model (includes SVEC) with some kind of long-run restr
 2436   # before ML estimation, we need to take into account the
 2437   # permanent/temporary shock classification in SVEC
 2438   if type == 4
 2439     matrix locjalpha = obj.jalpha
 2440     # trim beta from restr. exo terms if needed:
 2441     matrix beta = (obj.jcase % 2 == 0) ? obj.jbeta[1:n,] : obj.jbeta
 2442     r = cols(beta)
 2443     matrix C1 = C1mat(obj.VARpar, 1, locjalpha, beta)
 2444     matrix J = zeros(n-r, r) | I(r)
 2445     matrix permatransRd = (J ** nullspace(locjalpha'))' ~ zeros(r*(n-r), 1)
 2446     fullRd = fullRd | permatransRd  # put the constraints together
 2447   else # can be type 2 (or 1 if obj.calc_lr)
 2448     matrix C1 = C1mat(obj.VARpar, 0) # compute the lr matrix
 2449   endif
 2450   if lr_constr
 2451     # TODO: Here we could also check whether in a SVEC model
 2452     # the additional long-run constraints are really related
 2453     # to the permanent shocks.
 2454     # (Otherwise they will probably be either redundant or
 2455     # impossible. Or should we do that in the ident check?)
 2456     fullRd = fullRd | lrConstr(C1, obj.Rd1l) # put the constraints together
 2457   endif
 2458   if verbosity &gt; 1
 2459     matrix lrSigma = qform(C1, obj.Sigma) # compute the lr cov matrix
 2460     printf &quot;Long-run matrix (C1): \n%8.3f\n&quot;, C1
 2461     printf &quot;Long-run Sigma: \n%8.3f\n&quot;, lrSigma
 2462   endif
 2463   # store the C1 matrix for possible future use
 2464   matrix obj.C1 = C1
 2465 endif
 2466 return fullRd
 2467 </code>
 2468 </gretl-function>
 2469 <gretl-function name="SampleMatrices" type="matrices" private="1">
 2470  <params count="4">
 2471   <param name="Ra" type="matrix"/>
 2472   <param name="da" type="matrix"/>
 2473   <param name="Rb" type="matrix"/>
 2474   <param name="db" type="matrix"/>
 2475  </params>
 2476 <code>/*
 2477 Returns an array with A and B evaluated at a random point
 2478 in the parameter space (useful for checking the constraint
 2479 matrices).
 2480 Parameters: Ra, da, Rb, db = constraint matrices;
 2481 */
 2482 n = sqrt(cols(Ra | Rb))
 2483 matrices ret
 2484 matrix tmp = imp2exp(Ra ~ da)
 2485 if cols(tmp) &gt; 1
 2486   tmp *= muniform(cols(tmp)-1,1) | 1
 2487 endif
 2488 ret += mshape(tmp,n,n)
 2489 matrix tmp = imp2exp(Rb ~ db)
 2490 if cols(tmp) &gt; 1
 2491   tmp *= muniform(cols(tmp)-1,1) | 1
 2492 endif
 2493 ret += mshape(tmp,n,n)
 2494 return ret
 2495 </code>
 2496 </gretl-function>
 2497 <gretl-function name="Dtn" type="matrix" private="1">
 2498  <params count="1">
 2499   <param name="n" type="int"/>
 2500  </params>
 2501 <code>/*
 2502 Creates the matrix \tilde{D}_n.
 2503 Output has n^2 rows and n*(n-1)/2 columns; any (n x n) skewsymmetric
 2504 matrix has a vectorised form which lies in the space spanned by the
 2505 columns of \tilde{D}_n
 2506 */
 2507 p = round(n*(n-1)/2)
 2508 matrix A = zeros(1,n) | (lower(unvech(seq(1,p)')) ~ zeros(n-1,1))
 2509 matrix B = zeros(n^2,p)
 2510 matrix C
 2511 loop i=1..p --quiet
 2512   C = A.=i
 2513   B[,i] = vec(C .- C')
 2514 endloop
 2515 return B
 2516 </code>
 2517 </gretl-function>
 2518 <gretl-function name="Umat" type="matrix" private="1">
 2519  <params count="4">
 2520   <param name="n" type="int"/>
 2521   <param name="R" type="matrix"/>
 2522   <param name="d" type="matrix"/>
 2523   <param name="S" type="matrix"/>
 2524  </params>
 2525 <code># See Lucchetti(2006) -- left-hand side of matrix T;
 2526 # see eq. 26 on page 248
 2527 p = cols(S)
 2528 matrix ret = {}
 2529 loop i=1..p --quiet
 2530   ret |= R * (mshape(S[,i],n,n)' ** I(n))
 2531 endloop
 2532 return ret
 2533 </code>
 2534 </gretl-function>
 2535 <gretl-function name="Tmat" type="matrix" private="1">
 2536  <params count="4">
 2537   <param name="n" type="int"/>
 2538   <param name="R" type="matrix"/>
 2539   <param name="d" type="matrix"/>
 2540   <param name="S" type="matrix"/>
 2541  </params>
 2542 <code># See Lucchetti(2006) -- right-hand side of matrix T;
 2543 # see eq. 26 on page 248
 2544 p = cols(S)
 2545 matrix ret = {}
 2546 loop i=1..p --quiet
 2547   ret |= R * (I(n) ** mshape(S[,i],n,n))
 2548 endloop
 2549 return ret
 2550 </code>
 2551 </gretl-function>
 2552 <gretl-function name="strucond" type="scalar" private="1">
 2553  <params count="6">
 2554   <param name="n" type="int"/>
 2555   <param name="Ra" type="matrix"/>
 2556   <param name="da" type="matrix"/>
 2557   <param name="Rb" type="matrix"/>
 2558   <param name="db" type="matrix"/>
 2559   <param name="verbose" type="int" default="0"/>
 2560  </params>
 2561 <code>/*
 2562 Checks the structure condition and optionally
 2563 prints out some of the relevant matrices:
 2564 Parameters: Ra, da, Rb, db = constraint matrices;
 2565 iprint = output mode:
 2566 */
 2567 matrix Sa = imp2exp(Ra~da)
 2568 matrix Sb = imp2exp(Rb~db)
 2569 if verbose &gt; 1
 2570   print Ra da Rb db Sa Sb
 2571 endif
 2572 matrix Ua = Umat(n,Ra,da,Sa)
 2573 matrix Ub = Umat(n,Rb,db,Sb)
 2574 matrix Tb = Tmat(n,Rb,db,Sb)
 2575 if verbose &gt; 1
 2576   print Ua Ub Tb
 2577 endif
 2578 matrix C = Ua ~ zeros(rows(Ua),(n*(n-1)/2))
 2579 C |= Ub ~ Tb * Dtn(n)
 2580 if verbose &gt; 1
 2581   print C
 2582 endif
 2583 /* purge zero rows from C */
 2584 matrix e = maxr(abs(C)) .&gt; 0
 2585 C = selifr(C,e)
 2586 if verbose &gt; 1
 2587   printf &quot;After filtering zero rows, C is %d x %d\n&quot;, rows(C), cols(C)
 2588 endif
 2589 matrix CC = C'C
 2590 if verbose &gt; 1
 2591   print CC
 2592 endif
 2593 d = det(CC)
 2594 if d == 0
 2595   matrix nspace = nullspace(CC)
 2596   u_rank = cols(nspace)
 2597   if verbose
 2598     loop i=1..u_rank --quiet
 2599       printf &quot;Q_%d = \n&quot;, i
 2600       printf &quot;%6.1f&quot;, mshape(nspace[1:n*n,i],n,n)
 2601       printf &quot;H_%d = \n&quot;, i
 2602       printf &quot;%6.1f&quot;, mshape(Dtn(n) * nspace[n*n+1:,i],n,n)
 2603     endloop
 2604   endif
 2605 endif
 2606 return (d &gt; 0)
 2607 </code>
 2608 </gretl-function>
 2609 <gretl-function name="ordercond" type="scalar" private="1">
 2610  <params count="4">
 2611   <param name="n" type="int"/>
 2612   <param name="Ra" type="matrix"/>
 2613   <param name="Rb" type="matrix"/>
 2614   <param name="verbose" type="int" default="0"/>
 2615  </params>
 2616 <code>/*
 2617 Checks the order condition and optionally
 2618 prints out some of the relevant matrices:
 2619 Parameters: Ra, Rb = constraint matrices;
 2620 iprint = output mode:
 2621 */
 2622 p = 2*n*n - (n + 1)*n/2
 2623 if verbose &gt; 0
 2624   print  &quot;Checking order condition:&quot;
 2625   printf &quot; no. of constraints on A:\t%d\n&quot;, rows(Ra)
 2626   printf &quot; no. of constraints on B /or C:\t%d\n&quot;, rows(Rb)
 2627   printf &quot; no. of total constraints:\t%d\n&quot;, rows(Ra) + rows(Rb)
 2628   printf &quot; no. of necessary constraints:\t%d\n&quot;, p
 2629 endif
 2630 return rank(Ra)+rank(Rb) &gt;= p
 2631 </code>
 2632 </gretl-function>
 2633 <gretl-function name="rankcond" type="scalar" private="1">
 2634  <params count="6">
 2635   <param name="n" type="int"/>
 2636   <param name="Ra" type="matrix"/>
 2637   <param name="da" type="matrix"/>
 2638   <param name="Rb" type="matrix"/>
 2639   <param name="db" type="matrix"/>
 2640   <param name="verbose" type="int" default="0"/>
 2641  </params>
 2642 <code>/*
 2643 Checks the rank condition as per Amisano-Giannini and optionally
 2644 prints out some of the relevant matrices:
 2645 Parameters: Ra, da, Rb db, = constraint matrices;
 2646 Note that in fact we check for the existence of solutions to
 2647 the system given in eq. (9), chapter 4. The condition discussed later
 2648 (matrix Q) is, sadly, wrong.
 2649 */
 2650 matrices AB = SampleMatrices(Ra, da, Rb, db)
 2651 matrix A = AB[1]
 2652 matrix B = AB[2]
 2653 matrix BB = B*B'
 2654 matrix Q11 =  Ra * (A' ** BB)
 2655 matrix Q21 = -Rb * (B' ** BB)
 2656 matrix Q22 = Q21 * Dtn(n)
 2657 matrix Q = (Q11 ~ zeros(rows(Ra), n*(n-1)/2)) | (Q21 ~ Q22)
 2658 scalar r = rank(Q)
 2659 if verbose &gt; 1
 2660   loop foreach m Q11 Q21 Q22 Q --quiet
 2661     printf &quot;\n$m:\n%7.2f&quot;, $m
 2662   endloop
 2663 endif
 2664 if verbose &gt; 0
 2665   print  &quot;Checking rank condition:&quot;
 2666   printf &quot; r = %d, cols(Q) = %d\n&quot;, r, cols(Q)
 2667 endif
 2668 return r == cols(Q)
 2669 </code>
 2670 </gretl-function>
 2671 <gretl-function name="ident" type="scalar" private="1">
 2672  <params count="2">
 2673   <param name="mod" type="bundleref"/>
 2674   <param name="verbose" type="int" default="0"/>
 2675  </params>
 2676 <code># signature was: matrix *Ra, matrix *da, matrix *Rb, matrix *db, int verbose[0]
 2677 /*
 2678 Main function for checking identification.
 2679 The function also rewrites the 'fullRd' matrix to the mod if empty, and adds or overwrites 'cleanfullRd' or 'Rd0' and 'Rd1'
 2680 if redundant restrictions are found.
 2681   The algorithm is described fully in Lucchetti (2006), &quot;Identification Of Covariance Structures&quot;, Econometric
 2682   Theory, Cambridge University Press, vol. 22(02), p 235-257.
 2683   (currently partly unused due to some issues, work in progress)
 2684   The Ra, da, Rb, db constraint matrices are constructed from mod.Rd0, mod.Rd1 and oher ingredients if long-run constraints play a role
 2685   (in C and SVEC models).
 2686   Return value 0 means: no identification (check failed).
 2687   */
 2688   n2 = mod.n * mod.n
 2689   if mod.type == 3 # AB - Model
 2690     matrix Ra = mod.Rd0[, 1:n2] # restrictions on A
 2691     matrix da = mod.Rd0[, n2+1]
 2692     matrix Rb = mod.Rd1[, 1:n2] # restrictions on B
 2693     matrix db = mod.Rd1[, n2+1]
 2694     # check the constraints on A for inconsistencies
 2695     # and/or redundancies
 2696     # (for B further down, along with C)
 2697     err = CheckNormalizeRd(&amp;Ra, &amp;da)
 2698     if err == 1
 2699       printf &quot;--Contradictory constraints on A!--\n&quot;
 2700       return 0
 2701     elif err == 2
 2702       if verbose
 2703         print &quot;Warning: redundant constraints on A, dropping some.&quot;
 2704       endif
 2705       mod.Rd0 = Ra ~ da
 2706     endif
 2707   else             # plain, C-model or SVEC
 2708     if !nelem(mod.fullRd)
 2709       matrix mod.fullRd = get_full_Rd(&amp;mod, 0)
 2710       if !nelem(mod.fullRd) 	# possible only for C model?
 2711         funcerr &quot;No restrictions specified!&quot;
 2712       endif
 2713     endif
 2714     matrix Rb = mod.fullRd[, 1:n2]
 2715     matrix db = mod.fullRd[, n2+1]
 2716     matrix Ra = I(n2)	# just needed as a dummy below
 2717     matrix da = vec(I(mod.n))
 2718     if mod.type == 4    # special SVEC check
 2719       if verbose
 2720         print &quot; (The identification check for SVEC models is work in&quot;
 2721         print &quot;  progress. Be careful especially when loadings alpha&quot;
 2722         print &quot;  are restricted.)&quot;
 2723       endif
 2724       transshockcheck(mod)  # based on Luetkepohl 2008
 2725     endif
 2726   endif
 2727   # check the constraints on B and C (respectively) for inconsistencies
 2728   # and/or redundancies
 2729   err = CheckNormalizeRd(&amp;Rb, &amp;db)
 2730   if err == 1
 2731     printf &quot;--Contradictory constraints on B (or C, respectively)!--\n&quot;
 2732     return 0
 2733   elif err == 2
 2734     if verbose
 2735       print &quot;Warning: redundant constraints on B (or C), dropping some.&quot;
 2736     endif
 2737     if mod.type == 3
 2738       mod.Rd1 = Rb ~ db
 2739     else
 2740       matrix mod.cleanfullRd = Rb ~ db
 2741     endif
 2742   endif
 2743   # removed redundant calc: n = round(sqrt(cols(Ra | Rb)))
 2744   /* check for the order condition */
 2745   scalar id_o = ordercond(mod.n, Ra, Rb, verbose)
 2746   if verbose
 2747     string printout = id_o ? &quot;OK&quot; : &quot;fails!&quot;
 2748     printf &quot;=&gt; Order condition %s\n&quot;, printout
 2749   endif
 2750   # /* check for the structure condition */
 2751   # scalar id_s = strucond(n, locRa, locda, locRb, locdb, verbose)
 2752   # string printout = id_s ? &quot;OK&quot; : &quot;fails!&quot;
 2753   # printf &quot;Structure condition %s\n&quot;, printout
 2754   /* check for the rank condition */
 2755   scalar id_r = rankcond(mod.n, Ra, da, Rb, db, verbose)
 2756   if verbose
 2757     string printout = id_r ? &quot;OK&quot; : &quot;fails!&quot;
 2758     printf &quot;=&gt; Rank condition %s\n\n&quot;, printout
 2759   endif
 2760   return (id_o &amp;&amp; id_r)
 2761 </code>
 2762 </gretl-function>
 2763 <gretl-function name="scoreC" type="void" private="1">
 2764  <params count="3">
 2765   <param name="ret" type="matrixref"/>
 2766   <param name="theta" type="matrixref"/>
 2767   <param name="dat" type="matrixref"/>
 2768  </params>
 2769 <code>matrix Sigma Ss
 2770 unscramble_dat(&amp;dat, &amp;Sigma, &amp;Ss)
 2771 scalar n = rows(Sigma)
 2772 scalar npar = cols(Ss) - 1
 2773 matrix C = mat_exp(theta, Ss) # was: , 0)
 2774 matrix S = Ss[,1:npar]
 2775 matrix iC = inv(C)
 2776 matrix ret = iC' (qform(iC,Sigma) - I(n))
 2777 ret = vec(ret)'S
 2778 </code>
 2779 </gretl-function>
 2780 <gretl-function name="C1mat" type="matrix" private="1">
 2781  <params count="4">
 2782   <param name="A" type="matrix" const="true"/>
 2783   <param name="VECM" type="bool" default="0"/>
 2784   <param name="jalpha" type="matrix" optional="true" const="true"/>
 2785   <param name="jbeta" type="matrix" optional="true" const="true"/>
 2786  </params>
 2787 <code># (Jan 2018: used to be pointer args, coming from a time where it was necessary
 2788 # for optional args; now changed to non-pointers)
 2789 /*
 2790 computes C(1) out of the autoregressive matrices; jalpha and jbeta
 2791 are alpha and beta, which are used only if VECM is nonzero, that
 2792 is under cointegration
 2793 (Remark Sven: here jbeta was assumed to be (n,r)
 2794 without coeffs for restr. terms
 2795 Different from jbeta (n+1,r) in other contexts (?).
 2796 We now (Jan 2018) make sure only rows 1:n are used. */
 2797 scalar n = rows(A)
 2798 scalar p = cols(A) / n
 2799 matrix tmp = mshape(vec(A), n*n, p)
 2800 if VECM == 0
 2801   tmp = mshape(sumr(tmp), n, n)
 2802   matrix ret = inv(I(n) - tmp)
 2803 else # cointegrated
 2804   if !exists(jalpha) || !exists(jbeta)
 2805     funcerr &quot;Need cointegration params for C1 in VECM!&quot;
 2806   endif
 2807   matrix aperp = nullspace(jalpha')
 2808   matrix bperp = nullspace(jbeta[1:n, ]')
 2809   tmp = mshape(tmp*seq(1,p)', n, n)
 2810   matrix ret = bperp * inv(aperp'tmp*bperp) * aperp'
 2811 endif
 2812 return ret
 2813 </code>
 2814 </gretl-function>
 2815 <gretl-function name="init_C" type="matrix" private="1">
 2816  <params count="2">
 2817   <param name="Sigma" type="matrix"/>
 2818   <param name="Rd" type="matrix"/>
 2819  </params>
 2820 <code>matrix Ss = imp2exp(Rd)
 2821 scalar k = cols(Ss)
 2822 if k == 1
 2823   # nothing to estimate
 2824   ret = {}
 2825 else
 2826   scalar n = rows(Sigma)
 2827   matrix S = (k&gt;1) ? Ss[,1:k-1] : {}
 2828   matrix s = Ss[,k]
 2829   matrix K = cholesky(Sigma)
 2830   matrix bigmat = (K ** I(n))
 2831   matrix ret = mols(vec(Sigma) - bigmat*s, bigmat*S)
 2832 endif
 2833 return ret
 2834 </code>
 2835 </gretl-function>
 2836 <gretl-function name="PseudoHessC" type="void" private="1">
 2837  <params count="3">
 2838   <param name="H" type="matrixref"/>
 2839   <param name="theta" type="matrixref"/>
 2840   <param name="dat" type="matrixref"/>
 2841  </params>
 2842 <code>matrix Sigma Ss
 2843 unscramble_dat(&amp;dat, &amp;Sigma, &amp;Ss)
 2844 scalar n = rows(Sigma)
 2845 scalar npar = cols(Ss) - 1
 2846 matrix S = Ss[,1:npar]
 2847 # printf &quot;PseudoHessC\n%10.4f\n%4.1f\n&quot;, theta, Ss
 2848 matrix C = mat_exp(theta, Ss) # was: , 0)
 2849 H = InfoMat(C, S) # was InfoMatC
 2850 </code>
 2851 </gretl-function>
 2852 <gretl-function name="estC" type="matrix" private="1">
 2853  <params count="7">
 2854   <param name="theta" type="matrixref"/>
 2855   <param name="Sigma" type="matrix"/>
 2856   <param name="Rd" type="matrix"/>
 2857   <param name="vcv" type="matrixref" optional="true"/>
 2858   <param name="errcode" type="scalarref"/>
 2859   <param name="method" type="int"/>
 2860   <param name="verbose" type="scalar" default="1"/>
 2861  </params>
 2862 <code>matrix Ss = imp2exp(Rd)
 2863 scalar n = rows(Sigma)
 2864 if cols(Ss) == 1
 2865   # nothing to estimate
 2866   matrix C = mshape(Ss, n, n)
 2867   if exists(vcv)
 2868     vcv = zeros(n*n,n*n)
 2869   endif
 2870   errcode = 0
 2871   return C
 2872 endif
 2873 scalar SCALE = 0 # EXPERIMENTAL
 2874 scalar npar = rows(theta)
 2875 if SCALE == 1
 2876   printf &quot;Scale!\n&quot;
 2877   matrix s = sqrt(diag(Sigma))
 2878   matrix sSig = Sigma ./ (s*s')
 2879 else
 2880   matrix sSig = Sigma
 2881 endif
 2882 # printf &quot;estC\n%10.4f\n%4.1f\n&quot;, theta, Ss
 2883 matrix tmp = mat_exp(theta, Ss) # was: , 1)
 2884 if verbose &gt; 2
 2885   /* obsolete ? */
 2886   printf &quot;check within estC -- before estimation\n&quot;
 2887   check_const(tmp, Rd)
 2888 endif
 2889 matrix dat = vec(sSig) ~ Ss
 2890 matrix g H
 2891 if verbose &gt; 1
 2892   set max_verbose 1
 2893 else
 2894   set max_verbose 0
 2895 endif
 2896 err = 1
 2897 iters = 0
 2898 #set bfgs_toler 1.0e-03
 2899 matrix theta0 = theta
 2900 errcode = 0
 2901 loop while (err==1 &amp;&amp; iters&lt;100) --quiet
 2902   if method == 0
 2903     catch scalar ll = BFGSmax(theta, &quot;loglik(&amp;theta, &amp;dat, -1)&quot;)
 2904     errcode = $error
 2905     scoreC(&amp;g, &amp;theta, &amp;dat)
 2906   elif method == 1
 2907     catch scalar ll = BFGSmax(theta, &quot;loglik(&amp;theta, &amp;dat, -1)&quot;, &quot;scoreC(&amp;g, &amp;theta, &amp;dat)&quot;)
 2908     errcode = $error
 2909   elif method == 2
 2910     catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, -1)&quot;)
 2911     errcode = $error
 2912     scoreC(&amp;g, &amp;theta, &amp;dat)
 2913   elif method == 3
 2914     catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, -1)&quot;, &quot;scoreC(&amp;g, &amp;theta, &amp;dat)&quot;)
 2915     errcode = $error
 2916   elif method == 4
 2917     catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, -1)&quot;, &quot;scoreC(&amp;g, &amp;theta, &amp;dat)&quot;, &quot;PseudoHessC(&amp;H, &amp;theta, &amp;dat)&quot;)
 2918     errcode = $error
 2919   endif
 2920   if errcode&gt;0
 2921     printf &quot;errcode = %d\n&quot;, errcode
 2922   endif
 2923   scalar crit = maxr(abs(g))
 2924   err = (crit &gt; 1.0e-4)
 2925   if err==1
 2926     iters++
 2927     theta = 0.1*mnormal(npar,1) + theta0
 2928     if verbose&gt;1
 2929       printf &quot;Iter %3d: Restarting... ll = %12.7f, crit = %16.10f\n&quot;, iters, ll, crit
 2930       printf &quot;theta = %10.5f grad = %10.5f\n&quot;, theta', g
 2931     endif
 2932   endif
 2933 endloop
 2934 scalar crit = maxr(abs(g))
 2935 warn = (crit &gt; 1.0e-1)
 2936 if !err
 2937   if (iters &gt; 0) &amp;&amp; (verbose &gt; 1)
 2938     printf &quot;Converged after %d restarts\n&quot;, iters
 2939   endif
 2940   # printf &quot;estC(2)\n%10.4f\n%4.1f\n&quot;, theta, Ss
 2941   matrix C = mat_exp(theta, Ss) # was: , 1)
 2942   if SCALE == 1
 2943     C = C .* s'
 2944   endif
 2945   if verbose &gt; 1
 2946     printf &quot;Estimated C matrix:\n%12.5f&quot;, C
 2947   endif
 2948   if exists(vcv)
 2949     vcv = coeffVCV(Ss[,1:npar], &amp;C)
 2950   endif
 2951   if verbose &gt; 1
 2952     matrix Info = InfoMat(C, Ss[,1:npar])	# was InfoMatC
 2953     printf &quot;estC : Info = \n%14.10f\n&quot;, Info
 2954     printf &quot;estC : score = \n%14.10f\n&quot;, g
 2955   endif
 2956 else
 2957   if verbose &gt; 1
 2958     printf &quot;No convergence! :-( \n%12.6f\n&quot;, theta' | g
 2959   endif
 2960   matrix C = {}
 2961   errcode = 1
 2962 endif
 2963 return C
 2964 </code>
 2965 </gretl-function>
 2966 <gretl-function name="is_standard_AB" type="scalar" private="1">
 2967  <params count="2">
 2968   <param name="aRd" type="matrix"/>
 2969   <param name="bRd" type="matrix"/>
 2970  </params>
 2971 <code>test = has_unit_diag(aRd) &amp;&amp; has_free_diag(bRd)
 2972 return test
 2973 </code>
 2974 </gretl-function>
 2975 <gretl-function name="scoreAB" type="void" private="1">
 2976  <params count="4">
 2977   <param name="ret" type="matrixref"/>
 2978   <param name="theta" type="matrixref"/>
 2979   <param name="dat" type="matrixref"/>
 2980   <param name="p1" type="int"/>
 2981  </params>
 2982 <code># FIXME: This func might benefit from getting the pre-computed C matrix
 2983 matrix Sigma Ss A B
 2984 unscramble_dat(&amp;dat, &amp;Sigma, &amp;Ss)
 2985 p2 = rows(theta)
 2986 matrix aSs = Ss[,1:p1+1]
 2987 matrix bSs = Ss[,p1+2:p2+2]
 2988 n = rows(Sigma)
 2989 ABmat_exp(theta, aSs, bSs, &amp;A, &amp;B)
 2990 matrix iA = inv(A)
 2991 matrix C = iA*B
 2992 if p1&gt;0
 2993   matrix S = aSs[,1:p1]
 2994   S = -((C') ** iA) * S
 2995 else
 2996   matrix S = {}
 2997 endif
 2998 if p2&gt;p1
 2999   S ~= bSs[,1:cols(bSs)-1]
 3000 endif
 3001 matrix iC = inv(C)
 3002 matrix ret = iC' (qform(iC,Sigma) - I(n))
 3003 ret = vec(ret)'S
 3004 </code>
 3005 </gretl-function>
 3006 <gretl-function name="stdAB_init" type="matrix" private="1">
 3007  <params count="3">
 3008   <param name="mX" type="matrix"/>
 3009   <param name="aRd" type="matrix"/>
 3010   <param name="bRd" type="matrix"/>
 3011  </params>
 3012 <code># mX should contain the VAR residuals
 3013 matrix aSs = imp2exp(aRd)
 3014 matrix bSs = imp2exp(bRd)
 3015 n = cols(mX)
 3016 p = cols(aSs) - 1
 3017 if p &gt; 0
 3018   matrix S = aSs[,1:p]
 3019   matrix s = aSs[,p+1]
 3020   matrix sel = vec(transp(mshape(seq(1,n*n),n,n)))
 3021   matrix dep = vec(mX)
 3022   matrix reg = (I(n) ** mX)*S[sel,]
 3023   matrix e
 3024   matrix b = -mols(dep, reg, &amp;e)
 3025   e = mshape(e, rows(mX), n)
 3026   e = sqrt(meanc(e.^2))'
 3027   S = bSs[,1:cols(bSs)-1]
 3028   sel = 1 + seq(0,n-1)*(n+1)
 3029   matrix a = zeros(n*n,1)
 3030   a[sel] = e
 3031   e = S'a
 3032   matrix ret = b | e
 3033 else
 3034   matrix Sigma = (mX'mX) / rows(mX)
 3035   matrix ret = init_C(Sigma, bRd)
 3036 endif
 3037 return ret
 3038 </code>
 3039 </gretl-function>
 3040 <gretl-function name="nonstdAB_init" type="matrix" private="1">
 3041  <params count="3">
 3042   <param name="E" type="matrix"/>
 3043   <param name="aRd" type="matrix"/>
 3044   <param name="bRd" type="matrix"/>
 3045  </params>
 3046 <code>matrix aSs = imp2exp(aRd)
 3047 matrix bSs = imp2exp(bRd)
 3048 T = rows(E)
 3049 n = cols(E)
 3050 matrix Sigma = (E'E)./T
 3051 matrix C = cholesky(Sigma)
 3052 startA = (rows(aRd) &gt; rows(bRd))
 3053 #startA = 1 # provisional
 3054 ka = cols(aSs) - 1
 3055 if ka&gt;0
 3056   matrix Sa = aSs[,1:ka]
 3057 endif
 3058 matrix sa = aSs[,ka+1]
 3059 kb = cols(bSs) - 1
 3060 if kb&gt;0
 3061   matrix Sb = bSs[,1:kb]
 3062 endif
 3063 matrix sb = bSs[,kb+1]
 3064 if startA == 1
 3065   matrix giSb = invpd(Sb'Sb)*Sb'
 3066   matrix tmp = giSb*(C' ** I(n))
 3067   if ka&gt;0
 3068     matrix giSa = invpd(Sa'Sa)*Sa'
 3069     matrix gama = 0.1*ones(ka,1)
 3070     matrix gamb = tmp * (Sa*gama + sa) - giSb*sb
 3071     tmp = giSa*(inv(C)' ** I(n))
 3072     gama = tmp * (Sb*gamb + sb) - giSa*sa
 3073   else
 3074     matrix gama = {}
 3075     matrix gamb = tmp*sa - Sb'sb
 3076   endif
 3077 else
 3078   matrix giSa = invpd(Sa'Sa)*Sa'
 3079   matrix tmp = giSa*(inv(C)' ** I(n))
 3080   if kb&gt;0
 3081     matrix giSb = invpd(Sb'Sb)*Sb'
 3082     matrix gamb = 0.1*ones(kb,1)
 3083     matrix gama = tmp * (Sb*gamb + sb) - giSa*sa
 3084     tmp = giSb*(C' ** I(n))
 3085     gamb = tmp * (Sa*gama + sa) - giSb*sb
 3086   else
 3087     matrix gama = tmp*sb - Sa'sa
 3088     matrix gamb = {}
 3089   endif
 3090 endif
 3091 if kb &gt; 0
 3092   scale = tr(Sigma) / tr(C*C')
 3093   printf &quot;nonstdAB_init: Scale = %g\n&quot;, scale
 3094   gamb = sqrt(scale) .* gamb
 3095   /*
 3096 elif ka&gt;0
 3097   scalar scale = tr(Sigma) / tr(C*C')
 3098   printf &quot;nonstdAB_init: Scale = %g\n&quot;, scale
 3099   gama = sqrt(scale) ./ gama
 3100   */
 3101 endif
 3102 matrix ret = gama | gamb
 3103 return ret
 3104 </code>
 3105 </gretl-function>
 3106 <gretl-function name="ABmat_exp" type="void" private="1">
 3107  <params count="5">
 3108   <param name="theta" type="matrix"/>
 3109   <param name="aSs" type="matrix"/>
 3110   <param name="bSs" type="matrix"/>
 3111   <param name="A" type="matrixref"/>
 3112   <param name="B" type="matrixref"/>
 3113  </params>
 3114 <code>scalar p1 = cols(aSs) - 1
 3115 scalar p2 = rows(theta)
 3116 matrix theta_a = (p1&gt;0) ? theta[1:p1]: {}
 3117 matrix A = mat_exp(theta_a, aSs)
 3118 matrix theta_b = (p2&gt;p1) ? theta[p1+1:p2] : {}
 3119 matrix B = mat_exp(theta_b, bSs) # was: , 1)
 3120 </code>
 3121 </gretl-function>
 3122 <gretl-function name="PseudoHessAB" type="void" private="1">
 3123  <params count="4">
 3124   <param name="H" type="matrixref"/>
 3125   <param name="theta" type="matrixref"/>
 3126   <param name="dat" type="matrixref"/>
 3127   <param name="p1" type="int"/>
 3128  </params>
 3129 <code>matrix Sigma Ss A B
 3130 unscramble_dat(&amp;dat, &amp;Sigma, &amp;Ss)
 3131 scalar p2 = rows(theta)
 3132 matrix aSs = Ss[,1:p1+1]
 3133 matrix bSs = Ss[,p1+2:p2+2]
 3134 scalar n = rows(Sigma)
 3135 scalar n2 = n*n
 3136 ABmat_exp(theta, aSs, bSs, &amp;A, &amp;B)
 3137 scalar ka = cols(aSs) - 1
 3138 scalar kb = cols(bSs) - 1
 3139 matrix S = zeros(2*n2, ka+kb)
 3140 if ka&gt;0
 3141   S[1:n2,1:ka] = aSs[,1:ka]
 3142 endif
 3143 if kb&gt;0
 3144   S[n2+1:,ka+1:ka+kb] = bSs[,1:kb]
 3145 endif
 3146 H = InfoMat(B, S, &amp;A)
 3147 </code>
 3148 </gretl-function>
 3149 <gretl-function name="estAB" type="matrix" private="1">
 3150  <params count="9">
 3151   <param name="theta" type="matrixref"/>
 3152   <param name="Sigma" type="matrix"/>
 3153   <param name="aRd" type="matrix"/>
 3154   <param name="bRd" type="matrix"/>
 3155   <param name="vcv" type="matrixref" optional="true"/>
 3156   <param name="errcode" type="scalarref"/>
 3157   <param name="method" type="int"/>
 3158   <param name="verbose" type="int" default="1"/>
 3159   <param name="transfABkakb" type="matricesref" optional="true"/>
 3160  </params>
 3161 <code># new optional arg transfABkakb added by Sven for transfer/
 3162 # to avoid calc duplication
 3163 matrix aSs  = imp2exp(aRd)
 3164 matrix bSs  = imp2exp(bRd)
 3165 scalar npar = rows(theta)
 3166 scalar n    = rows(Sigma)
 3167 scalar p1 = cols(aSs) - 1
 3168 scalar p2 = p1 + cols(bSs) - 1
 3169 matrix dat = vec(Sigma) ~ aSs ~ bSs
 3170 matrix g H
 3171 if verbose &gt; 1
 3172   set max_verbose 1
 3173 else
 3174   set max_verbose 0
 3175 endif
 3176 scalar err = 0
 3177 if method == 0
 3178   catch scalar ll = BFGSmax(theta, &quot;loglik(&amp;theta, &amp;dat, p1)&quot;)
 3179   err = $error
 3180   scoreAB(&amp;g, &amp;theta, &amp;dat, p1)
 3181 elif method == 1
 3182   catch scalar ll = BFGSmax(theta, &quot;loglik(&amp;theta, &amp;dat, p1)&quot;, &quot;scoreAB(&amp;g, &amp;theta, &amp;dat, p1)&quot;)
 3183   err = $error
 3184 elif method == 2
 3185   catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, p1)&quot;)
 3186   err = $error
 3187   scoreAB(&amp;g, &amp;theta, &amp;dat, p1)
 3188 elif method == 3
 3189   catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, p1)&quot;, &quot;scoreAB(&amp;g, &amp;theta, &amp;dat, p1)&quot;)
 3190   err = $error
 3191 elif method == 4
 3192   catch scalar ll = NRmax(theta, &quot;loglik(&amp;theta, &amp;dat, p1)&quot;, &quot;scoreAB(&amp;g, &amp;theta, &amp;dat, p1)&quot;, &quot;PseudoHessAB(&amp;H, &amp;theta, &amp;dat, p1)&quot; )
 3193   err = $error
 3194 endif
 3195 scalar crit = maxr(abs(g))
 3196 scalar warn = (crit &gt; 1.0e-1)
 3197 if err || warn
 3198   matrix C = {}
 3199   errcode = err + 10*warn
 3200   outfile stderr --write
 3201   printf &quot;err = %d; warn = %d; Gradient: %10.6f&quot;, err, warn, g
 3202   outfile --close
 3203   return C
 3204 endif
 3205 matrix A B
 3206 ABmat_exp(theta, aSs, bSs, &amp;A, &amp;B)
 3207 # use matlab-style &quot;matrix division&quot; for speed and accuracy
 3208 matrix C = A\B
 3209 if verbose &gt; 1
 3210   printf &quot;Estimated A matrix:\n%12.5f&quot;, A
 3211   printf &quot;Estimated B matrix:\n%12.5f&quot;, B
 3212   printf &quot;Estimated C matrix:\n%12.5f&quot;, C
 3213   printf &quot;Check:\n%12.5f&quot;, cholesky(Sigma)
 3214 endif
 3215 # new transplanted block from SVAR_estimate():
 3216 n2 = n*n
 3217 ka = cols(aSs) - 1
 3218 kb = cols(bSs) - 1
 3219 matrix S = zeros(2*n2, ka+kb)
 3220 if ka &gt; 0
 3221   S[1:n2,1:ka] = aSs[,1:ka]
 3222 endif
 3223 if kb &gt; 0
 3224   S[n2+1:2*n2,ka+1:ka+kb] = bSs[,1:kb]
 3225 endif
 3226 if exists(vcv)
 3227   vcv = coeffVCV(S, &amp;B, &amp;A)  # was coeffVCV(S, &amp;mB, &amp;mA)
 3228 endif
 3229 # transfer back the A, B, ka, kb results (ugly hack &lt;Sven&gt;)
 3230 if exists(transfABkakb)
 3231   transfABkakb = defarray(A, B, {ka}, {kb})
 3232 endif
 3233 return C
 3234 </code>
 3235 </gretl-function>
 3236 <gretl-function name="doIRF" type="void" private="1">
 3237  <params count="1">
 3238   <param name="SVARobj" type="bundleref"/>
 3239  </params>
 3240 <code>/*
 3241 constructs the structural VMA representation. Note
 3242 that the companion matrix is never used explicitly;
 3243 The output is not returned by the function, but rather
 3244 put into the bundle under the &quot;IRFs&quot; key.
 3245 */
 3246 scalar type = SVARobj.type
 3247 matrix varA = SVARobj.VARpar
 3248 scalar H = SVARobj.horizon + 1
 3249 scalar n = SVARobj.n
 3250 if (type == 1) || (type == 2) || (type == 4)
 3251   matrix C = SVARobj.C	# was: SVARobj.S1
 3252 elif type == 3
 3253   if inbundle(SVARobj, &quot;C&quot;)	# maybe not yet computed
 3254     matrix C = SVARobj.C
 3255   else
 3256     matrix C = SVARobj.S1 \ SVARobj.S2
 3257   endif
 3258 endif
 3259 matrix ret = zeros(H,n*n)
 3260 scalar np = SVARobj.p * n
 3261 matrix tmp = I(np)
 3262 matrix prd = zeros(np,np)
 3263 loop i=1..H --quiet
 3264   ret[i,] = vec(tmp[1:n,1:n] * C)'
 3265   if (np&gt;n)
 3266     prd[n+1:np, ] = tmp[1:np-n, ]
 3267   endif
 3268   prd[1:n,] = varA * tmp
 3269   tmp = prd
 3270 endloop
 3271 if SVARobj.ncumul &gt; 0
 3272   # The following code is now done in SVAR_cumulate
 3273   # once and for all:
 3274   # matrix to_cum = SVARobj.cumul
 3275   # tmp = zeros(n,n)
 3276   # tmp[to_cum,] = 1
 3277   # sel = selifr(transp(seq(1,n*n)), vec(tmp))
 3278   ret[, SVARobj.cumsel] = cum(ret[, SVARobj.cumsel]) # .cumsel was sel
 3279 endif
 3280 matrix SVARobj.IRFs = ret
 3281 </code>
 3282 </gretl-function>
 3283 <gretl-function name="check_bounds" type="scalar" private="1">
 3284  <params count="3">
 3285   <param name="s" type="int"/>
 3286   <param name="v" type="int"/>
 3287   <param name="n" type="int"/>
 3288  </params>
 3289 <code># negative s is allowed and means to flip the shock
 3290 ret = 0
 3291 if abs(s) &gt; n || (s == 0)
 3292   ret = 1
 3293 elif (v &gt; n) || (v &lt; 1)
 3294   ret = 2
 3295 endif
 3296 return ret
 3297 </code>
 3298 </gretl-function>
 3299 <gretl-function name="IRFgrph" type="string" private="1">
 3300  <params count="11">
 3301   <param name="IRFmat" type="matrix" const="true"/>
 3302   <param name="snum" type="int"/>
 3303   <param name="vnum" type="int"/>
 3304   <param name="scale" type="scalar"/>
 3305   <param name="sname" type="string"/>
 3306   <param name="vname" type="string"/>
 3307   <param name="keypos" type="int" min="0" max="2" default="1"/>
 3308   <param name="cumul" type="matrixref" optional="true"/>
 3309   <param name="boot" type="bundleref" optional="true"/>
 3310   <param name="bc" type="int" min="0" max="2" default="0"/>
 3311   <param name="whichdraw" type="int" min="0" default="0"/>
 3312  </params>
 3313 <code># bc: bias correction choice
 3314 # FIXME: We should get rid of sname/vname args and get them automatically
 3315 # via snum/vnum
 3316 # (the number of args here is horrible...)
 3317 if whichdraw &amp;&amp; (exists(boot) || bc)
 3318   funcerr &quot;cannot have SR/whichdraw and bootstrap options together&quot;
 3319 endif
 3320 flip = (snum &lt; 0)
 3321 snum = abs(snum)
 3322 n = round(sqrt(cols(IRFmat)))
 3323 # bootrep = exists(boot) ? boot.rep : 0 # doesn't work due to gretl bug+
 3324 bootrep = 0
 3325 if exists(boot)
 3326   bootrep = boot.rep
 3327 endif
 3328 # anything to cumulate?
 3329 cumulate = exists(cumul) ? sumr(sumc(cumul .= vnum)) &gt; 0 : 0
 3330 tmpfname = int(1000 * muniform(1,1))
 3331 string tmpfile = $windows ? sprintf(&quot;@dotdir\\irf%d.gp&quot;, tmpfname) : sprintf(&quot;@dotdir/irf%d.gp&quot;, tmpfname)
 3332 k = (snum-1) * n + vnum
 3333 h = rows(IRFmat)
 3334 matrix x = flip ? -IRFmat[,k]./scale : IRFmat[,k]./scale	# matrix?
 3335 if bootrep &gt; 0
 3336   matrix hicb = flip ? -boot.lo_cb : boot.hi_cb	# matrix?
 3337   matrix locb = flip ? -boot.hi_cb : boot.lo_cb
 3338   matrix mdn  = flip ? -boot.mdns  : boot.mdns
 3339   matrix locb = locb[,k] ./ scale
 3340   matrix hicb = hicb[,k] ./ scale
 3341   matrix mdn  = mdn[,k]  ./ scale
 3342   scalar miny = minc((locb | x))
 3343   scalar maxy = maxc((hicb | x))
 3344 else
 3345   scalar miny = minc(x)
 3346   scalar maxy = maxc(x)
 3347 endif
 3348 miny = (miny&gt;0) ? 0 : miny
 3349 maxy = (maxy&lt;0) ? 0 : maxy
 3350 set force_decpoint on
 3351 ## write the raw gnuplot code
 3352 # TODO: replace outfile --write ... outfile --close
 3353 # with outfile ... end outfile eventually (but needs gretl 2018b)
 3354 outfile &quot;@tmpfile&quot; --write
 3355 printf &quot;set yrange [%g:%g]\n&quot;, miny*1.05, maxy*1.05
 3356 printf &quot;set xrange [%g:%g]\n&quot;, 0, (h-1)*1.05
 3357 printf &quot;set xzeroaxis\n&quot;
 3358 string l_sname = sname==&quot;&quot; ? sprintf(&quot;%d&quot;, snum) : sname
 3359 string l_vname = vname==&quot;&quot; ? sprintf(&quot;%d&quot;, vnum) : vname
 3360 printf &quot;set title \&quot;IRF: %s -&gt; %s&quot;, l_sname, l_vname
 3361 if bc == 1
 3362   printf &quot;; bias-correction = partial&quot;
 3363 elif bc == 2
 3364   printf &quot;; bias-correction = full&quot;
 3365 elif whichdraw
 3366   printf &quot; (draw %d)&quot;, whichdraw
 3367 endif
 3368 if cumulate &gt; 0
 3369   printf &quot; (cumulated)&quot;
 3370 endif
 3371 printf &quot;\&quot;\n&quot;
 3372 if keypos == 0
 3373   printf &quot;set key off\n&quot;
 3374 elif keypos == 1
 3375   printf &quot;set key outside\n&quot;
 3376 elif keypos == 2
 3377   printf &quot;set key below\n&quot;
 3378 endif
 3379 if bootrep &gt; 0
 3380   printf &quot;set style fill solid 0.125\n&quot;
 3381   printf &quot;plot '-' using 1:2:3 w filledcurve t &quot;
 3382   printf &quot;'Bstrap %d%% CI', \\\n&quot;, floor(100 * boot.alpha)
 3383   printf &quot;'-' w l lw 1 lc 1 t 'Bstrap median', \\\n&quot;
 3384   printf &quot;'-' w l lw 2 lc -1 t 'IRF'\n&quot;
 3385   loop i=1..h -q
 3386     printf &quot;%d\t%g\t%g\n&quot;, i-1, locb[i,], hicb[i,]
 3387   endloop
 3388   printf &quot;e\n&quot;
 3389   loop i=1..h -q
 3390     printf &quot;%d\t%g\n&quot;, i-1, mdn[i,]
 3391   endloop
 3392   printf &quot;e\n&quot;
 3393 else
 3394   printf &quot;plot '-' w l lw 2\n&quot;
 3395 endif
 3396 loop i=1..h -q
 3397   printf &quot;%d\t%g\n&quot;, i-1, x[i]
 3398 endloop
 3399 printf &quot;e\n&quot;
 3400 outfile --close
 3401 set force_decpoint off
 3402 return tmpfile
 3403 </code>
 3404 </gretl-function>
 3405 <gretl-function name="FEVDgrph" type="string" private="1">
 3406  <params count="6">
 3407   <param name="FEVDmat" type="matrix" const="true"/>
 3408   <param name="v" type="scalar"/>
 3409   <param name="vname" type="string"/>
 3410   <param name="snames" type="strings"/>
 3411   <param name="keypos" type="int" min="0" max="2" default="1"/>
 3412   <param name="titletail" type="string" optional="true"/>
 3413  </params>
 3414 <code># Final arg 'titletail' is a hack originally meant to put the draw number
 3415 # in the output in the SR case. (In principle it can be anything.)
 3416 if !exists(titletail)
 3417   string titletail = &quot;&quot;
 3418 endif
 3419 n = round(sqrt(cols(FEVDmat)))
 3420 h = rows(FEVDmat) - 1
 3421 scalar tmpfname = int(10000*muniform(1,1))
 3422 if $windows
 3423   string datfile = sprintf(&quot;@dotdir\\fevd%d.txt&quot;, tmpfname)
 3424   # sprintf datfile &quot;@dotdir\\fevd%d.txt&quot;, tmpfname
 3425   string tmpfile = sprintf(&quot;@dotdir\\fevd%d.gp&quot;, tmpfname)
 3426   # sprintf tmpfile &quot;@dotdir\\fevd%d.gp&quot;, tmpfname
 3427 else
 3428   string datfile = sprintf(&quot;@dotdir/fevd%d.txt&quot;, tmpfname)
 3429   # sprintf datfile &quot;@dotdir/fevd%d.txt&quot;, tmpfname
 3430   string tmpfile = sprintf(&quot;@dotdir/fevd%d.gp&quot;, tmpfname)
 3431   # sprintf tmpfile &quot;@dotdir/fevd%d.gp&quot;, tmpfname
 3432 endif
 3433 matrix sel = (v-1)*n + seq(1,n)
 3434 set force_decpoint on
 3435 outfile &quot;@tmpfile&quot; --write
 3436 if   keypos == 0
 3437   printf &quot;set key off\n&quot;
 3438 elif keypos == 1
 3439   printf &quot;set key outside\n&quot;
 3440 elif keypos == 2
 3441   printf &quot;set key below\n&quot;
 3442 endif
 3443 printf &quot;set yrange [0:100]\n&quot;
 3444 printf &quot;set xrange [%g:%g]\n&quot;, 0, h
 3445 printf &quot;set xzeroaxis\n&quot;
 3446 printf &quot;set title \&quot;FEVD for %s%s\&quot;\n&quot;, vname, titletail
 3447 printf &quot;set style fill solid 0.25\n&quot;
 3448 printf &quot;set style histogram rowstacked\n&quot;
 3449 printf &quot;set style data histogram\n&quot;
 3450 loop i=1..n --quiet
 3451   string sname = snames[i]
 3452   if i == 1
 3453     printf &quot;plot '%s' using 2 t '%s', \\\n&quot;, datfile, sname
 3454   elif i == n
 3455     printf &quot;\t'' using %d t '%s'\n&quot;, i+1, sname
 3456   else
 3457     printf &quot;\t'' using %d t '%s', \\\n&quot;, i+1, sname
 3458   endif
 3459 endloop
 3460 outfile --close
 3461 outfile &quot;@datfile&quot; --write
 3462 printf &quot;%12.4f\n&quot;, seq(0, h)' ~ 100*FEVDmat[,sel]
 3463 outfile --close
 3464 set force_decpoint off
 3465 return tmpfile
 3466 </code>
 3467 </gretl-function>
 3468 <gretl-function name="boot_printout" type="void" private="1">
 3469  <params count="5">
 3470   <param name="type" type="int"/>
 3471   <param name="n" type="int"/>
 3472   <param name="rep" type="int"/>
 3473   <param name="failed" type="int"/>
 3474   <param name="Spar_mat" type="matrix" const="true"/>
 3475  </params>
 3476 <code># Sep 2020: change from pointerized Spar_mat to const
 3477 matrix bm = meanc(Spar_mat)
 3478 matrix bV = mcov(Spar_mat)
 3479 scalar nc = cols(Spar_mat)
 3480 scalar n2 = n*n
 3481 # force numerical zeros to real zeros
 3482 e = diag(bV) .&lt; 1.0e-12
 3483 if maxc(e)
 3484   e = selifc(seq(1, nc), e')
 3485   bV[e,] = 0
 3486   bV[,e] = 0
 3487 endif
 3488 printf &quot;Bootstrap results (%d replications, %d failed)\n&quot;, rep + failed, failed
 3489 if (type != 3) &amp;&amp; ( nc == 2 * n2 )	# was cols(Spar_mat)
 3490   # Long-run matrix at the end exists!
 3491   # And so this is a model with long-run restrictions
 3492   # (Bl-Quah-style or SVEC),
 3493   # or the user forced it via calc_lr.
 3494   matrix bK = mshape( bm[1:n2], n, n)
 3495   printStrMat(bK, bV[1:n2, 1:n2], &quot;C&quot;)
 3496   matrix bL = mshape( bm[n2+1:], n, n)
 3497   printStrMat(bL, bV[1+n2: , 1+n2: ], &quot;LongRun&quot;)
 3498 elif type != 3
 3499   # C model without printing the long-run matrix
 3500   # (SVEC / type 4 should not happen here, because it has
 3501   # long-run constraints by construction)
 3502   matrix bK = mshape(bm,n,n)
 3503   printStrMat(bK, bV, &quot;C&quot;)
 3504 elif type == 3	# AB model
 3505   matrix bmA = mshape( bm[1:n2], n, n)
 3506   printStrMat(bmA, bV[1:n2, 1:n2], &quot;A&quot;)
 3507   matrix bmB = mshape( bm[n2+1:], n, n)
 3508   printStrMat(bmB, bV[1+n2:,1+n2:], &quot;B&quot;)
 3509 else
 3510   funcerr &quot;shouldn't happen&quot;
 3511 endif
 3512 </code>
 3513 </gretl-function>
 3514 <gretl-function name="bias_correction" type="matrix" private="1">
 3515  <params count="4">
 3516   <param name="b" type="bundle" const="true"/>
 3517   <param name="Y0" type="matrix" const="true"/>
 3518   <param name="BC" type="matrixref"/>
 3519   <param name="bparams" type="bundle" optional="true" const="true"/>
 3520  </params>
 3521 <code>/* This function implements a bias correction for
 3522 the estimate of the VAR parameters as per Kilian, REStat (1998).
 3523 (The SVEC case is not allowed here, but it must
 3524 be checked on the outside.)
 3525 Sep 2019: re-did interface to use bundle; the new
 3526 'boottype' member means:
 3527 1: standard resampling (default, as before)
 3528 2-4: wild bootstrap (3 variants)
 3529 5: moving-blocks bootstrap
 3530 (So the check below for what we need in the bundle.)
 3531 */
 3532 if !inbundle(b,&quot;BCiter&quot;) || !inbundle(b,&quot;VARpar&quot;) || !inbundle(b,&quot;E&quot;) || !inbundle(b,&quot;X&quot;) || !inbundle(b,&quot;boottype&quot;) || !inbundle(b,&quot;bmu&quot;)
 3533   funcerr &quot;needed input missing in bundle arg&quot;
 3534 endif
 3535 # check for stationarity first
 3536 scalar maxmod = max_eval(b.VARpar)
 3537 if maxmod &lt; 0.9999
 3538   matrix innov = prepres(b.E, b.boottype)
 3539   bundle bootstuff = defbundle(&quot;Y0&quot;,Y0, &quot;innov&quot;,innov)
 3540   if exists(bparams)
 3541     bootstuff.bparams = bparams
 3542   endif
 3543   matrix Ab = avg_VARpar_boot(b, bootstuff)[1]
 3544   matrix BC = b.VARpar - Ab
 3545   add_and_smash(&amp;Ab, BC) 	 # was H = ..., unused
 3546 else	# not stationary
 3547   matrix Ab = b.VARpar
 3548 endif
 3549 return Ab
 3550 </code>
 3551 </gretl-function>
 3552 <gretl-function name="avg_VARpar_boot" type="matrices" private="1">
 3553  <params count="2">
 3554   <param name="b" type="bundle" const="true"/>
 3555   <param name="stuff" type="bundle" const="true"/>
 3556  </params>
 3557 <code># In principle we could call base_est() below in the loop,
 3558 # because that's what it's for, but since we only need
 3559 # the reduced-form AR coefficients, it's not worth it...
 3560 # This function returns an array (matrices) just in case in the
 3561 # future we want to apply the parallel_func MPI approach to it.
 3562 # (This is a requirement there.)
 3563 # Of course in essence it's still just a single matrix.
 3564 bundle bparams = inbundle(stuff, &quot;bparams&quot;) ? stuff.bparams : null
 3565 # simulate the average reduced-form coeffs
 3566 matrix Absum = zeros(b.n, b.n*b.p)
 3567 matrix lagseq = seq(1,b.p)
 3568 matrix Uinit = zeros(b.p, b.n)
 3569 if rows(stuff.Y0) != b.p
 3570   funcerr &quot;wrong assumption about length of Y initvals here!&quot;
 3571 endif
 3572 loop i=1 .. b.BCiter --quiet
 3573   matrix U = Uinit | drawbootres(stuff.innov, b.boottype, bparams)
 3574   if cols(b.bmu) &gt; 0	# was mu
 3575     U += b.bmu
 3576   endif
 3577   matrix bY  = varsimul(b.VARpar, U[b.p + 1: ,], stuff.Y0) # was rows(stuff.Y0)+1
 3578   matrix reg = b.X ~ mlag(bY, lagseq)
 3579   matrix Pi  = mols(bY[b.p + 1: ,], reg[b.p + 1: ,])
 3580   Absum += transp(Pi[b.k + 1: b.k + b.n*b.p, ])
 3581 endloop
 3582 return defarray(Absum ./ b.BCiter)
 3583 </code>
 3584 </gretl-function>
 3585 <gretl-function name="calc_bmu" type="void" private="1">
 3586  <params count="1">
 3587   <param name="obj" type="bundleref"/>
 3588  </params>
 3589 <code># disentangle determ/exog
 3590 # Sep 2020: add the result directly to the (pointerized) bundle
 3591 matrix bmu = zeros(obj.T, obj.n)
 3592 if obj.k &amp;&amp; obj.type == 4 # SVEC,
 3593   # (with some unrestr. exo apart from const/trend)
 3594   if obj.jcase == 1
 3595     bmu = obj.X * obj.mu
 3596   elif obj.jcase == 2 || obj.jcase == 3
 3597     bmu = obj.X * obj.mu[2:, ]
 3598     # need to add restr. or unrestr. const later
 3599   elif obj.jcase == 4 || obj.jcase == 5
 3600     bmu = obj.X * obj.mu[3:, ]
 3601     # need to add restr./unr. const &amp; trend later
 3602   endif
 3603 elif obj.k    # no SVEC
 3604   bmu = obj.X * obj.mu    # this was the pre-1.4 handled case
 3605 endif
 3606 # more special treatment of SVEC
 3607 if obj.type == 4
 3608   # add constant
 3609   if obj.jcase &gt; 2      # unrestricted
 3610     bmu = bmu .+  obj.mu[1, ]  # (use broadcasting)
 3611   elif obj.jcase == 2   # restricted
 3612     bmu = bmu .+ (obj.jbeta[obj.n + 1, ] * obj.jalpha')
 3613   endif
 3614   # add trend
 3615   if obj.jcase == 4	# restricted
 3616     bmu += seq(1, obj.T)' obj.jbeta[obj.n + 1, ] * obj.jalpha'
 3617   elif obj.jcase == 5 # unrestricted
 3618     bmu += seq(1, obj.T)' obj.mu[2, ]
 3619   endif
 3620 endif
 3621 matrix obj.bmu = bmu
 3622 </code>
 3623 </gretl-function>
 3624 <gretl-function name="prepres" type="matrix" private="1">
 3625  <params count="2">
 3626   <param name="E" type="matrix"/>
 3627   <param name="boottype" type="int" min="1" max="5" default="1">
 3628 <description>bootstrap type</description>
 3629 <labels count="5">
 3630 "resampling" "wildN" "wildR" "wildM" "moving blocks" </labels>
 3631   </param>
 3632  </params>
 3633 <code># The input is demeaned (column-wise) only for
 3634 # the resampling (type 1). It isn't done together with the
 3635 # new draws to avoid repeating it many times.
 3636 # (Cannot have 'const matrix E' as gretl versions &lt;2019d would
 3637 #  crash when doing 'return E'.)
 3638 # (Copying the E matrix w/o doing anything shouldn't hurt
 3639 #  performance too much because it is only done once.)
 3640 if boottype == 1
 3641   return cdemean(E)
 3642 else
 3643   return E
 3644 endif
 3645 </code>
 3646 </gretl-function>
 3647 <gretl-function name="drawbootres" type="matrix" private="1">
 3648  <params count="3">
 3649   <param name="E" type="matrix" const="true"/>
 3650   <param name="boottype" type="int" min="1" max="5" default="1">
 3651 <description>bootstrap type</description>
 3652 <labels count="5">
 3653 "resampling" "wildN" "wildR" "wildM" "moving blocks" </labels>
 3654   </param>
 3655   <param name="bparams" type="bundle" optional="true" const="true">
 3656 <description>further inputs</description>
 3657   </param>
 3658  </params>
 3659 <code>/*
 3660 Construct a new draw of residuals for bootstrapping, where E
 3661 is a Txn matrix.
 3662 E can be original residuals or can be some pre-processed
 3663 input. (See the prepres() function - the pre-processing is
 3664 not done here to avoid doing it repeatedly.)
 3665 Currently boottype can go from 1 to 5:
 3666 1: traditional residual resampling (using gretl's resample())
 3667 2-4: wild bootstrap with the help of a standard normal, Rademacher, Mammen distributions respectively
 3668 5: RBMBB Brüggemann/Jentsch/Trenkler 2016. A block length
 3669 must be chosen.
 3670 Other bootstrap types may be added in the future.
 3671 (The E input could be empty then, e.g. for a purely parametric-
 3672 distribution bootstrap. The bundle bparams is intended to be
 3673 used then to specify parameters, in a way that is to be determined.)
 3674 */
 3675 # process extra input
 3676 bl = 0 	# signals default data-based choice
 3677 if exists(bparams) &amp;&amp; boottype == 5
 3678   if inbundle(bparams, &quot;moveblocklen&quot;)
 3679     bl = bparams.moveblocklen
 3680   endif
 3681 endif
 3682 ## standard
 3683 if boottype == 1
 3684   return resample(E)
 3685 endif
 3686 T = rows(E)
 3687 ## wild
 3688 if boottype &lt;= 4
 3689   if boottype == 2
 3690     # Normal
 3691     matrix w = mnormal(T,1)
 3692   elif boottype == 3
 3693     # Rademacher
 3694     matrix w = muniform(T,1) .&lt; 0.5 ? 1 : -1
 3695   elif boottype == 4
 3696     # Mammen
 3697     scalar s5 = sqrt(5)
 3698     scalar p = (0.5/s5) * (s5 + 1)
 3699     matrix w = 0.5 + (muniform(T,1) .&lt; p ? -s5/2 : s5/2)
 3700   endif
 3701   return E .* w
 3702   ## residual-based moving blocks
 3703 elif boottype == 5
 3704   if bl &gt; T
 3705     print &quot;Warning: block length too large, falling back to auto&quot;
 3706     bl = 0
 3707   endif
 3708   if bl == 0	# use default
 3709     bl = xmax(2, floor(T/10))
 3710   endif
 3711   # necessary number of blocks
 3712   s = ceil(T/bl)
 3713   # blocks starting points
 3714   matrix c = mrandgen(i, 1, T-bl+1, 1, s)
 3715   # convert to indices of all needed obs
 3716   matrix ndx = vec(c .+ seq(0, bl-1)')
 3717   ## recentring
 3718   matrix m = mshape(NA, bl, cols(E))
 3719   # calculate respective averages (bl different ones)
 3720   loop i = 1..bl -q
 3721     m[i,] = meanc(E[i: T-bl+i,])
 3722   endloop
 3723   ndx = ndx[1:T]	# cut off &quot;overhanging tail&quot;
 3724   return E[ndx,] - (ones(s,1) ** m)[1:T, ]
 3725 else
 3726   funcerr &quot;Shouldn't happen&quot;
 3727 endif
 3728 </code>
 3729 </gretl-function>
 3730 <gretl-function name="boottypechoice" type="void" private="1">
 3731  <params count="2">
 3732   <param name="mod" type="bundleref"/>
 3733   <param name="btypestr" type="string"/>
 3734  </params>
 3735 <code>temp = boottypecode(btypestr)
 3736 if temp == 0
 3737   print &quot;Warning: ignoring unrecognized SVAR bootstrap type.&quot;
 3738 else
 3739   mod.boottype = temp
 3740 endif
 3741 </code>
 3742 </gretl-function>
 3743 <gretl-function name="maybe_do_biascorr" type="void" private="1">
 3744  <params count="2">
 3745   <param name="bobj" type="bundleref"/>
 3746   <param name="bparams" type="bundle" const="true"/>
 3747  </params>
 3748 <code># this function adds &quot;ABCorA&quot; and perhaps &quot;Psi&quot; to bobj
 3749 if bobj.biascorr &amp;&amp; (bobj.type != 4)  # not available w/unit roots # BIASCORR
 3750   matrix Psi = {}
 3751   matrix start = bobj.Y[1: bobj.p, ]
 3752   matrix bobj.ABCorA = bias_correction(bobj, start, &amp;Psi, bparams) # ABC # moved up from inside the loop
 3753   matrix bobj.Psi = Psi
 3754 else
 3755   matrix bobj.ABCorA = bobj.VARpar
 3756 endif
 3757 </code>
 3758 </gretl-function>
 3759 <gretl-function name="SVAR_boot_innerloop" type="matrices" private="1">
 3760  <params count="3">
 3761   <param name="bobj" type="bundleref"/>
 3762   <param name="obj" type="bundle" const="true"/>
 3763   <param name="bparams" type="bundle" const="true"/>
 3764  </params>
 3765 <code># will return a three-matrix array:
 3766 # 1: bootirfs
 3767 # 2: Spar_mat
 3768 # 3: number of failures (1x1 scalar)
 3769 # copy
 3770 type = obj.type
 3771 n2   = obj.n * obj.n
 3772 # Prepare the residuals
 3773 matrix innov = prepres(obj.E, obj.boottype)
 3774 ## prepare output
 3775 matrix bootirfs = zeros(obj.nboot, (obj.horizon + 1) * n2)
 3776 # Spar_mat: the result matrix
 3777 # (-- type==10 (set-ID) cannot/should not happen in here --)
 3778 # need more cols if saving either A,B (for type 3) or C and the long-run matrix
 3779 numcols = (type &gt; 2) || (rows(obj.Rd1l) || obj.calc_lr) ? 2*n2 : n2
 3780 matrix Spar_mat = zeros(obj.nboot, numcols)
 3781 i = 1
 3782 failed = 0
 3783 set loop_maxiter 16384
 3784 loop while i &lt;= obj.nboot --quiet
 3785   # clear previous bootstrap indicator
 3786   bobj.step = 0
 3787   # generate bootstrap disturbances (bmu may be zero)
 3788   matrix U = obj.bmu[obj.p + 1:, ] + drawbootres(innov, obj.boottype, bparams)
 3789   # generate bootstrap data and store it in bootstrap object
 3790   matrix bobj.Y = varsimul(bobj.ABCorA, U, obj.Y[1: obj.p, ] )
 3791   # estimate VAR parameters, special treatment VECM/SVEC
 3792   if type == 4
 3793     vecm_est(&amp;bobj)
 3794   else
 3795     base_est(&amp;bobj)
 3796   endif
 3797   matrix bA = bobj.VARpar  # estimates (first n rows of companion mat)
 3798   matrix bSigma = bobj.Sigma
 3799   matrix theta = obj.theta # init original SVAR params
 3800   # (C/A&amp;B apparently, in suitable form...)
 3801   errcode = 0
 3802   ## Full bias correction
 3803   /* (The bc-ed VARpar need to be done before the new C matrix is
 3804   calculated at least if there are long-run constraints, because
 3805   then they enter via C1 through fullRd into C.
 3806   Otherwise the new C only depends on Sigma. - BTW, we do not update
 3807   Sigma in the full biascorr case. (In theory we could, by
 3808   re-calculating the residuals using the bc-ed VARpar.
 3809   We shouldn't, should we?)
 3810   */
 3811   if obj.biascorr == 2 &amp;&amp; type != 4 # only for non-SVEC
 3812     scalar H = add_and_smash(&amp;bA, bobj.Psi)
 3813     if ok(H)
 3814       bobj.VARpar = bA
 3815     else
 3816       errcode = 101
 3817     endif
 3818   endif
 3819   /* now re-estimate C, according to model type */
 3820   if type == 1 # Cholesky
 3821     matrix K = cholesky(bSigma)
 3822   elif type == 2 || type == 4 # consolidate the SVEC case here
 3823     /* Watch out for the possibility of constraints that were
 3824     partly redundant originally.
 3825     In the AB-model the cleaned restrictions are inherited
 3826     from the original model, apart from the fact that they
 3827     should already be caught at specification time (in SVAR_restrict).
 3828     But with long-run constraints (C, and especially SVEC with
 3829     weakly exog.)
 3830     the restrictions depend on the simulated bootstrap data.
 3831     So the cleaning of the restrictions must be re-done here.
 3832     */
 3833     if type != 3 &amp;&amp; inbundle(obj, &quot;cleanfullRd&quot;)
 3834       bobj.fullRd = {}	 # reset
 3835       id = ident(&amp;bobj, 0) # re-creates fullRd and possibly cleanfullRd
 3836       if !id
 3837         printf &quot;Ident problem in bootstrap draw %d\n&quot;, i
 3838         funcerr &quot;Unexpected ID problem in bootstrap&quot;
 3839       endif
 3840     elif type == 4 || nelem(obj.Rd1l) # long-run constraints
 3841       bobj.fullRd = get_full_Rd(&amp;bobj, 0)	# update fullRd
 3842     endif
 3843     if inbundle(bobj, &quot;cleanfullRd&quot;)
 3844       matrix fullRd = bobj.cleanfullRd
 3845     else
 3846       matrix fullRd = bobj.fullRd
 3847     endif
 3848     matrix K = estC(&amp;theta, bSigma, fullRd, null, &amp;errcode, obj.optmeth, 0)
 3849   elif type == 3 # &quot;AB&quot;
 3850     matrices transferAB = array(2)
 3851     # matrix Rd2 = type==3 ? obj.Rd0 : obj.Rd1l
 3852     # new: get A,B instead of re-calc'ing it below
 3853     # (obj.Rd0 was Rd2 before, but redundant...)
 3854     matrix K = estAB(&amp;theta, bSigma, obj.Rd0, obj.Rd1, null, &amp;errcode, obj.optmeth, 0, &amp;transferAB)
 3855   endif
 3856   ## Process and store the simulated C results
 3857   if !errcode &amp;&amp; rows(K) == obj.n
 3858     bobj.step = 2
 3859     bobj.theta = theta
 3860     # we don't treat the AB-model specially here (no reason to)
 3861     maybe_flip_columns(obj.C, &amp;K)
 3862     if (type == 1) || (type == 2) || (type == 4)
 3863       bobj.C = K	# is used in doIRF()
 3864       Spar_mat[i, 1: n2] = vec(K)'
 3865       /* New Oct 2017: Also bootstrap the long-run matrix if wanted
 3866       Jan 2018: add type 4 */
 3867       if ( type &lt; 3 &amp;&amp; ( rows(bobj.Rd1l) || bobj.calc_lr ) ) || type == 4
 3868         # (a plain or C model w/ long-run constr, or user switch)
 3869         # long-run mat (C1 comes from get_full_Rd() above
 3870         # (except type 1)):
 3871         matrix C1 = (type == 2 || type == 4) ? bobj.C1 : C1mat(bobj.VARpar)
 3872         matrix bobj.lrmat = C1 * bobj.C
 3873         # attach it to the other bootstrap result
 3874         Spar_mat[i, n2+1 : 2*n2] = vec(bobj.lrmat)'
 3875       endif
 3876     elif type == 3
 3877       # (Sven): the following stuff comes from estAB above
 3878       bobj.S1 = transferAB[1]
 3879       bobj.S2 = transferAB[2]
 3880       Spar_mat[i,] = vec(bobj.S1)' ~ vec(bobj.S2)'
 3881     endif
 3882   endif
 3883   if !errcode &amp;&amp; rows(K) == obj.n
 3884     doIRF(&amp;bobj)
 3885     bootirfs[i,] = vec(bobj.IRFs)'
 3886     i++
 3887   else
 3888     failed++
 3889     outfile stderr --write
 3890     printf &quot;Iter %4d failed (error code = %d)\n&quot;, i, errcode
 3891     outfile --close	# TODO: should be replaced with 'end outfile' once we can require 2018b or so
 3892   endif
 3893 endloop
 3894 return defarray(bootirfs, Spar_mat, {failed})
 3895 </code>
 3896 </gretl-function>
 3897 <gretl-function name="modelstring" type="string" private="1">
 3898  <params count="1">
 3899   <param name="type" type="int"/>
 3900  </params>
 3901 <code>/*
 3902 The idea about the number 10 for sign restrictions is:
 3903 - First it's very different from the other types (haha)
 3904 - Secondly if later on we really support mixed models, like some sign restrictions and some zero restrictions, perhaps this could be expressed as the sum of the types.
 3905 (E.g. 12 == 2 + 10 would be zero restrictions on the C matrix
 3906 coupled with sign stuff.)
 3907 */
 3908 string ret = &quot;&quot;
 3909 if type == 1
 3910   ret = &quot;plain&quot;
 3911 elif type == 2
 3912   ret = &quot;C&quot;
 3913 elif type == 3
 3914   ret = &quot;AB&quot;
 3915 elif type == 4
 3916   ret = &quot;SVEC&quot;
 3917 elif type == 10
 3918   ret = &quot;SR&quot; # sign restrictions
 3919 endif
 3920 return ret
 3921 </code>
 3922 </gretl-function>
 3923 <gretl-function name="modeltype" type="scalar" private="1">
 3924  <params count="1">
 3925   <param name="s" type="string"/>
 3926  </params>
 3927 <code>scalar ret = 0
 3928 string s2 = toupper(s)	# case insensitive
 3929 if s2 == &quot;PLAIN&quot;
 3930   ret = 1
 3931 elif s2 == &quot;C&quot;
 3932   ret = 2
 3933 elif s2 == &quot;AB&quot;
 3934   ret = 3
 3935 elif s2 == &quot;SVEC&quot;
 3936   ret = 4
 3937 elif s2 == &quot;KPSW&quot;
 3938   print &quot;Please use 'SVEC' for cointegrated SVARs,&quot;
 3939   print &quot;the code 'KPSW' is deprecated (but still works for now).&quot;
 3940   ret = 4
 3941 elif s2 == &quot;SR&quot;	# sign restrictions
 3942   ret = 10
 3943 endif
 3944 return ret
 3945 </code>
 3946 </gretl-function>
 3947 <gretl-function name="optimstring" type="string" private="1">
 3948  <params count="1">
 3949   <param name="type" type="int"/>
 3950  </params>
 3951 <code>string ret = &quot;&quot;
 3952 if type == 0
 3953   ret = &quot;BFGS (numerical)&quot;
 3954 elif type == 1
 3955   ret = &quot;BFGS (analytical)&quot;
 3956 elif type == 2
 3957   ret = &quot;Newton-Raphson (numerical)&quot;
 3958 elif type == 3
 3959   ret = &quot;Newton-Raphson (analytical score)&quot;
 3960 elif type == 4
 3961   ret = &quot;Scoring algorithm&quot;
 3962 endif
 3963 return ret
 3964 </code>
 3965 </gretl-function>
 3966 <gretl-function name="btypestring" type="string" private="1">
 3967  <params count="1">
 3968   <param name="type" type="int"/>
 3969  </params>
 3970 <code>string ret = &quot;&quot;
 3971 if type == 1
 3972   ret = &quot;resampling&quot;
 3973 elif type == 2
 3974   ret = &quot;wild (Normal)&quot;
 3975 elif type == 3
 3976   ret = &quot;wild (Rademacher)&quot;
 3977 elif type == 4
 3978   ret = &quot;wild (Mammen)&quot;
 3979 elif type == 5
 3980   ret = &quot;MBB (moving blocks)&quot;
 3981 endif
 3982 return ret
 3983 </code>
 3984 </gretl-function>
 3985 <gretl-function name="boottypecode" type="scalar" private="1">
 3986  <params count="1">
 3987   <param name="btypestr" type="string"/>
 3988  </params>
 3989 <code>ret = 0
 3990 string s = tolower(btypestr) # case insensitive!
 3991 if s == &quot;resampling&quot; || s == &quot;resample&quot;
 3992   ret = 1
 3993 elif s == &quot;wildn&quot; || s == &quot;wild&quot;
 3994   ret = 2
 3995 elif s == &quot;wildr&quot;
 3996   ret = 3
 3997 elif s == &quot;wildm&quot;
 3998   ret = 4
 3999 elif s == &quot;mbb&quot; || s == &quot;moving blocks&quot;
 4000   ret = 5
 4001 endif
 4002 return ret
 4003 </code>
 4004 </gretl-function>
 4005 <gretl-function name="BCstring" type="string" private="1">
 4006  <params count="1">
 4007   <param name="level" type="int"/>
 4008  </params>
 4009 <code>string ret = &quot;&quot;
 4010 if level == 0
 4011   ret = &quot;none&quot;
 4012 elif level == 1
 4013   ret = &quot;partial&quot;
 4014 elif level == 2
 4015   ret = &quot;full&quot;
 4016 endif
 4017 return ret
 4018 </code>
 4019 </gretl-function>
 4020 <gretl-function name="set_default_dimensions" type="matrix" private="1">
 4021  <params count="4">
 4022   <param name="mod" type="bundleref"/>
 4023   <param name="lY" type="list"/>
 4024   <param name="lX" type="list" optional="true"/>
 4025   <param name="varorder" type="int"/>
 4026  </params>
 4027 <code># trivial stuff: set up dimensions etc.
 4028 n = nelem(lY)
 4029 k = nelem(lX)
 4030 /* no of endogenous variables */
 4031 mod.n = n
 4032 /* no of exogenous variables */
 4033 mod.k = k
 4034 /* VAR order */
 4035 mod.p = varorder
 4036 /* horizon for IRFs etc */
 4037 mod.horizon = 10
 4038 if $pd == 4
 4039   /* quarterly: try 5 years */
 4040   mod.horizon = 20
 4041 elif $pd == 12
 4042   /* monthly: two years */
 4043   mod.horizon = 24
 4044 endif
 4045 /* sample size */
 4046 list everything = lY lX
 4047 smpl everything --no-missing
 4048 matrix mreg = { everything }
 4049 mod.T = rows(mreg)
 4050 mod.t1 = $t1
 4051 mod.t2 = $t2
 4052 return mreg
 4053 </code>
 4054 </gretl-function>
 4055 <gretl-function name="base_est" type="scalar" private="1">
 4056  <params count="1">
 4057   <param name="SVARobj" type="bundleref"/>
 4058  </params>
 4059 <code>/*
 4060 takes a SVAR object and adds the basic VAR estimates
 4061 to the bundle; returns an errcode (virginity check).
 4062 */
 4063 scalar step = SVARobj.step
 4064 err = (step&gt;0)
 4065 if err
 4066   printf &quot;Base estimation done already!\n&quot;
 4067 else
 4068   matrix mY  = SVARobj.Y
 4069   scalar p   = SVARobj.p
 4070   scalar n   = SVARobj.n
 4071   scalar k   = SVARobj.k
 4072   scalar T   = SVARobj.T
 4073   matrix E
 4074   matrix mreg = SVARobj.X ~ mlag(mY, seq(1,p))
 4075   matrix olspar = mols(mY[p+1:T,], mreg[p+1:T,], &amp;E)
 4076   matrix Sigma = (E'E) ./ (T - (n*p + k)) # was:  / df
 4077   matrix SVARobj.VARpar = transp(olspar[k+1:k+n*p,])
 4078   matrix SVARobj.mu = (k&gt;0) ? olspar[1:k,] : {} # was: d
 4079   matrix SVARobj.Sigma = Sigma
 4080   matrix SVARobj.E = E
 4081   SVARobj.step = 1
 4082   SVARobj.LL0 = T * (n* log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
 4083 endif
 4084 return err
 4085 </code>
 4086 </gretl-function>
 4087 <gretl-function name="vecm_est" type="scalar" private="1">
 4088  <params count="1">
 4089   <param name="SVARobj" type="bundleref"/>
 4090  </params>
 4091 <code>/*
 4092 We can't afford to be too flexible here, and the intended usage
 4093 is as follows.
 4094 We assume the model has already been declared as a SVEC (type=4)
 4095 model. We also assume that the cointegration properties (beta mat plus
 4096 deterministics) have already been set up via the SVAR_coint() function, so that we already have beta and alpha available as &quot;jbeta&quot; and
 4097 &quot;jalpha&quot;, respectively. Finally, we take care of proper treatment of
 4098 deterministics, via the scalar &quot;jcase&quot; (1 to 5).
 4099 */
 4100 # --- preliminary checks
 4101 err = (SVARobj.step &gt; 0)
 4102 if err
 4103   printf &quot;Base estimation done already!\n&quot;
 4104   return err
 4105 endif
 4106 err = (SVARobj.type != 4)
 4107 if err
 4108   printf &quot;Wrong model type!\n&quot;
 4109   return err
 4110 endif
 4111 # --- grab stuff from the bundle
 4112 matrix mY     = SVARobj.Y
 4113 matrix jbeta   = SVARobj.jbeta
 4114 matrix jalpha  = SVARobj.jalpha
 4115 scalar p	  = SVARobj.p
 4116 scalar n	  = SVARobj.n
 4117 scalar k	  = SVARobj.k
 4118 scalar r      = cols(jbeta)
 4119 scalar dcase   = SVARobj.jcase
 4120 scalar ols_al = rows(jalpha) == 0
 4121 scalar T      = SVARobj.T
 4122 # --- first thing: do we have a pre-set value for alpha?
 4123 matrix dY = diff(mY)
 4124 matrix dep = dY[p+1:T,]
 4125 matrix E = {}
 4126 ng = n * (p-1) # number coming from Gammas /lagged diffs (p&gt;0 for SVEC)
 4127 # deterministics
 4128 matrix mreg = vecm_det(T, dcase)
 4129 # ECM terms
 4130 matrix ECM = mlag(mY * jbeta[1:n,], 1)
 4131 if dcase == 2
 4132   ECM = ECM .+ jbeta[n+1, ]
 4133 elif dcase == 4
 4134   ECM += seq(1,T)'jbeta[n+1, ]
 4135 endif
 4136 if ols_al
 4137   # alpha must be estimated together with the rest of the VECM
 4138   matrix mreg ~= SVARobj.X ~ ECM
 4139 else
 4140   matrix dep -= (ECM[p+1:T,] * jalpha')
 4141   matrix mreg ~= SVARobj.X
 4142 endif
 4143 # extra lags
 4144 if p &gt; 1
 4145   mreg ~= mlag(dY, seq(1,p-1))
 4146 endif
 4147 # trim the first rows
 4148 if rows(mreg)&gt;0
 4149   mreg = mreg[p+1:T,]
 4150   matrix olspar = mols(dep, mreg, &amp;E)
 4151 else
 4152   matrix olspar = {}
 4153   E = dep
 4154 endif
 4155 df = T - (n*p + k - (n-r)) # FIXME: aren't the unrestr. determ. terms missing?
 4156 # (but apparently df is never used again...)
 4157 matrix Sigma = (E'E)./T
 4158 # --- construct the various matrices required later
 4159 rp = rows(olspar)
 4160 # alpha first (should be before the gammas, if estimated)
 4161 if ols_al
 4162   jalpha = olspar[rp-ng-r+1 : rp-ng,]'
 4163 endif
 4164 # exogenous
 4165 if dcase == 1
 4166   matrix mu = {}
 4167   scalar nd = 0
 4168 elif dcase == 2
 4169   matrix mu = jbeta[n+1,] * jalpha'
 4170   scalar nd = 0
 4171 elif dcase == 3
 4172   matrix mu = olspar[1,]
 4173   scalar nd = 1
 4174 elif dcase == 4
 4175   matrix mu =  olspar[1,] |(jbeta[n+1,] * jalpha')
 4176   scalar nd = 1
 4177 elif dcase == 5
 4178   matrix mu = olspar[1:2,]
 4179   scalar nd = 2
 4180 endif
 4181 if k &gt; 0
 4182   mu = mu | olspar[nd + 1 : nd + k,] # was: mu ~ olspar[sel,]
 4183 endif
 4184 /*
 4185 companion form in levels
 4186 */
 4187 Pi = jalpha * jbeta[1:n,]'
 4188 if p &gt; 1
 4189   # the Gammas are always at the back:
 4190   ini = rp - ng + 1
 4191   matrix A = olspar[ini: ,] | zeros(n,n)	# was: ini:fin
 4192   A += I(n) + Pi' | -olspar[ini: ,]
 4193 else
 4194   matrix A = I(n) + Pi'
 4195 endif
 4196 matrix SVARobj.VARpar = A'
 4197 matrix SVARobj.mu = mu
 4198 matrix SVARobj.Sigma = Sigma
 4199 matrix SVARobj.E = E
 4200 matrix SVARobj.jalpha = jalpha
 4201 SVARobj.step = 1
 4202 SVARobj.LL0 = T*(n*log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
 4203 return err
 4204 </code>
 4205 </gretl-function>
 4206 <gretl-function name="VARloglik" type="scalar" private="1">
 4207  <params count="3">
 4208   <param name="T" type="scalar"/>
 4209   <param name="Sigma" type="matrix"/>
 4210   <param name="C" type="matrixref" optional="true"/>
 4211  </params>
 4212 <code># computes the (concentrated) loglikelihood for a VAR
 4213 # the matrix C (such that in theory CC' = Sigma) may be null,
 4214 # for just-identified models
 4215 n = rows(Sigma)
 4216 ll = n * 1.83787706640934548355 # ln(2*$pi)
 4217 if !exists(C) # just-identified model
 4218   ll += log(det(Sigma)) + n
 4219 else
 4220   matrix KK = invpd(C*C')
 4221   ll += -log(det(KK)) + tr(KK*Sigma)
 4222 endif
 4223 return -(T/2) * ll
 4224 </code>
 4225 </gretl-function>
 4226 <gretl-function name="loglik" type="scalar" private="1">
 4227  <params count="3">
 4228   <param name="theta" type="matrixref"/>
 4229   <param name="dat" type="matrixref"/>
 4230   <param name="modeltype" type="scalar"/>
 4231  </params>
 4232 <code># modeltype &lt; 0-&gt;C; modeltype&gt;=0: AB (contains free elements of a)
 4233 matrix Sigma = {}
 4234 matrix Ss = {}
 4235 unscramble_dat(&amp;dat, &amp;Sigma, &amp;Ss)
 4236 if modeltype &lt; 0
 4237   matrix C = mat_exp(theta, Ss) # was: , 0)
 4238 else
 4239   p1 = modeltype
 4240   p2 = rows(theta)
 4241   matrix aSs = Ss[,1:p1+1]
 4242   matrix bSs = Ss[,p1+2:p2+2]
 4243   matrix A B
 4244   ABmat_exp(theta, aSs, bSs, &amp;A, &amp;B)
 4245   # was: matrix C = B/A
 4246   matrix C = A\B
 4247 endif
 4248 matrix KK = invpd(C*C')
 4249 ll = det(KK) # should always be positive
 4250 ll = (ll&lt;=0) ? NA : -0.5 * (tr(KK*Sigma) - log(ll))
 4251 return ll
 4252 </code>
 4253 </gretl-function>
 4254 <gretl-function name="InfoMat" type="matrix" private="1">
 4255  <params count="3">
 4256   <param name="CorB" type="matrix"/>
 4257   <param name="S" type="matrix"/>
 4258   <param name="A" type="matrixref" optional="true"/>
 4259  </params>
 4260 <code>/*
 4261 merged from InfoMatC and InfoMatAB
 4262 First case C model (A is null), latter AB model.
 4263 */
 4264 matrix tmp = !exists(A) ? I(rows(CorB)) : (A\CorB) | -I(rows(A))
 4265 tmp = S' (tmp ** inv(CorB'))
 4266 return tmp * N2_ify(tmp')
 4267 </code>
 4268 </gretl-function>
 4269 <gretl-function name="coeffVCV" type="matrix" private="1">
 4270  <params count="3">
 4271   <param name="S" type="matrix"/>
 4272   <param name="rhs" type="matrixref"/>
 4273   <param name="lhs" type="matrixref" optional="true"/>
 4274  </params>
 4275 <code># C or AB model
 4276 matrix IM = !exists(lhs) ? InfoMat(rhs, S) : InfoMat(rhs, S, &amp;lhs)
 4277 # (should be correct with new InfoMat)
 4278 # quick-dirty check for singularity
 4279 if rcond(IM)&gt;1.0e-10
 4280   matrix iIM = invpd(IM)/$nobs
 4281 else
 4282   matrix evec
 4283   l = eigensym(IM, &amp;evec)
 4284   printf &quot;\n\nInformation matrix is not pd!!\n\n%12.5f\n&quot;, IM
 4285   printf &quot;Eigenvalues:\n%12.5f\n&quot;, l
 4286   matrix e = (l .&lt; 1.0e-07)
 4287   printf &quot;Troublesome eigenvectors:\n%12.5f\n&quot;, selifc(evec, e')
 4288   printf &quot;S:\n%4.0f\n&quot;, S
 4289   matrix iIM = zeros(rows(IM), cols(IM))
 4290 endif
 4291 return qform(S,iIM)
 4292 </code>
 4293 </gretl-function>
 4294 <gretl-function name="SVAR_est_printout" type="void" private="1">
 4295  <params count="1">
 4296   <param name="mod" type="bundleref"/>
 4297  </params>
 4298 <code># scalar type = mod.type
 4299 printf &quot;Optimization method = &quot;
 4300 strings os = defarray(&quot;BFGS (numerical)&quot;, &quot;BFGS (analytical)&quot;, &quot;Newton-Raphson (numerical)&quot;, &quot;Newton-Raphson (analytical score)&quot;, &quot;Scoring algorithm&quot;)
 4301 printf &quot;%s\n&quot;, os[mod.optmeth + 1]	# Because optmeth starts at 0 here
 4302 printf &quot;Unconstrained Sigma:\n%12.5f\n&quot;, mod.Sigma
 4303 if mod.type &lt; 3 || mod.type == 4 # plain or C-model, or (Jan 2018) SVEC
 4304   # Standard C matrix
 4305   printStrMat(mod.C, mod.vcv, &quot;C&quot;)	# was mod.S1
 4306   # Long-run matrix
 4307   if rows(mod.Rd1l) || mod.calc_lr || mod.type == 4
 4308     # either long-run restr., or user wish, or SVEC
 4309     # (this could be much embellished, probably)
 4310     printf &quot;Estimated long-run matrix&quot;
 4311     if rows(mod.Rd1l)
 4312       printf &quot; (restricted)&quot;
 4313     endif
 4314     printf &quot;\n&quot;
 4315     matrix longrun = mod.lrmat
 4316     # round to zero for printout (also do it in general?)
 4317     longrun = (abs(longrun) .&lt; 1e-15) ? 0.0 : longrun
 4318     print longrun
 4319   endif
 4320 elif mod.type == 3	# AB, new &lt;Sven&gt;
 4321   n2 = (mod.n)^2
 4322   if mod.ka &gt; 0
 4323     printStrMat(mod.S1, mod.vcv[1:n2, 1:n2], &quot;A&quot;)
 4324   endif
 4325   if mod.kb &gt; 0
 4326     printStrMat(mod.S2, mod.vcv[n2+1 : 2*n2, n2+1 : 2*n2], &quot;B&quot;)
 4327   endif
 4328 endif
 4329 printf &quot;  Log-likelihood = %g\n&quot;, mod.LL1
 4330 if inbundle(mod, &quot;LRoid&quot;)	# over-id test was done
 4331   printf &quot;  Overidentification LR test = %g (%d df, pval = %g)\n\n&quot;, mod.LRoid[1], mod.LRoid[2], mod.LRoid[3]
 4332 endif
 4333 </code>
 4334 </gretl-function>
 4335 <gretl-function name="putIrf_to_accdraw" type="void" private="1">
 4336  <params count="2">
 4337   <param name="SVARobj" type="bundleref"/>
 4338   <param name="whichdraw" type="int" min="0" default="0"/>
 4339  </params>
 4340 <code># Only for SR (type 10)
 4341 # to transfer the IRFs to the old format
 4342 if SVARobj.type != 10
 4343   funcerr &quot;this function only for SR type&quot;
 4344 elif !inbundle(SVARobj, &quot;acc_draws&quot;)
 4345   funcerr &quot;need accepted draws (acc_draws)&quot;
 4346 endif
 4347 n = SVARobj.n
 4348 h = SVARobj.horizon + 1
 4349 bundle pickdraw = SVARobj.acc_draws[whichdraw]
 4350 matrix IRFs = zeros(h, n*n)
 4351 loop ix = 1..h -q
 4352   IRFs[ix, ] = vec(pickdraw.irfs[ix])'
 4353 endloop
 4354 # copy to origin
 4355 matrix SVARobj.acc_draws[whichdraw].IRFs = IRFs
 4356 </code>
 4357 </gretl-function>
 4358 <gretl-function name="muVARparE_mayberedr" type="matrices" private="1">
 4359  <params count="2">
 4360   <param name="Smod" type="bundle" const="true"/>
 4361   <param name="Bofonedraw" type="matrix" optional="true"/>
 4362  </params>
 4363 <code># Retrieves all the quantities that perhaps are redrawn
 4364 # in the Bayesian set id case:
 4365 # - mu: exog. coeffs, B[1:k , ]
 4366 # - VARpar: autoregr. coeffs, B[k+1: ,]
 4367 # - E: residuals
 4368 redrawn = 0
 4369 if Smod.type == 10
 4370   if Smod.DO_BAYES &amp;&amp; !exists(Bofonedraw)
 4371     funcerr &quot;Need redrawn matrix input for Bayesian set id case&quot;
 4372   elif Smod.DO_BAYES
 4373     redrawn = 1
 4374   endif
 4375 endif
 4376 if !redrawn
 4377   # old standard cases
 4378   matrix mu = Smod.mu
 4379   matrix VARpar = Smod.VARpar
 4380   matrix resids = Smod.E    # original residuals
 4381 else
 4382   # need to use the resids associated with re-drawn coefficients
 4383   matrix mu = Bofonedraw[1 : Smod.k ,]
 4384   matrix VARpar = Bofonedraw[Smod.k + 1 : ,]
 4385   matrix resids = Smod.Y[Smod.p + 1 : ,] - Smod.mreg * Bofonedraw
 4386 endif
 4387 return defarray(mu, VARpar, resids)
 4388 </code>
 4389 </gretl-function>
 4390 <gretl-function name="errchkSRhisto" type="void" private="1">
 4391  <params count="2">
 4392   <param name="Smod" type="bundleref"/>
 4393   <param name="drawix" type="int"/>
 4394  </params>
 4395 <code>if !inbundle(Smod, &quot;bestdraw&quot;)
 4396   funcerr &quot;Model not properly initialized for SR type&quot;
 4397 elif !Smod.bestdraw &amp;&amp; !drawix
 4398   funcerr &quot;Need to pick one particular draw (run SRgetbest?)&quot;
 4399 elif !Smod.storeSRirfs
 4400   funcerr &quot;Need accepted draws content (acc_draws) in set-id case&quot;
 4401 elif drawix &gt; nelem(Smod.acc_draws)
 4402   printf &quot;Only %d accepted draws exist\n&quot;, nelem(Smod.acc_draws)
 4403   funcerr &quot;Draw index out of range&quot;
 4404 endif
 4405 </code>
 4406 </gretl-function>
 4407 <gretl-function name="errmsgshockmatch" type="void" private="1">
 4408  <params count="2">
 4409   <param name="pos" type="matrix" const="true"/>
 4410   <param name="sname" type="string"/>
 4411  </params>
 4412 <code>if nelem(pos) != 1
 4413   printf &quot;You don't have a shock named %s in your model.\n&quot;, sname
 4414   print &quot; (or several ones)&quot;
 4415   funcerr &quot;Specify correct and unique shock names&quot;
 4416 endif
 4417 </code>
 4418 </gretl-function>
 4419 <gretl-function name="safetycheck1" type="void" private="1">
 4420  <params count="3">
 4421   <param name="Sigma" type="matrix" const="true"/>
 4422   <param name="iXX" type="matrix" const="true"/>
 4423   <param name="B" type="matrix" const="true"/>
 4424  </params>
 4425 <code>hhh = 10
 4426 matrix V = Sigma ** iXX
 4427 chk = B[1:hhh,1] ~ sqrt(diag(V[1:hhh, 1:hhh]))
 4428 strings pnames = array(hhh)
 4429 loop i = 1..hhh -q
 4430   pnames[i] = sprintf(&quot;chk%d&quot;, i)
 4431 endloop
 4432 modprint chk pnames
 4433 </code>
 4434 </gretl-function>
 4435 <gretl-function name="get_n_exotic" type="scalar" private="1">
 4436  <params count="1">
 4437   <param name="mod" type="bundle" const="true"/>
 4438  </params>
 4439 <code>n_exotic = 0
 4440 if inbundle(mod, &quot;exoticSR&quot;)
 4441   n_exotic = nelem(mod.exoticSR.checks)
 4442 endif
 4443 return n_exotic
 4444 </code>
 4445 </gretl-function>
 4446 <gretl-function name="rand_rotation" type="matrix" private="1">
 4447  <params count="2">
 4448   <param name="n" type="scalar"/>
 4449   <param name="sel" type="matrix" const="true"/>
 4450  </params>
 4451 <code>scalar nk = cols(sel)
 4452 if nk == 0
 4453   matrix ret = qrdecomp(mnormal(n, n))
 4454 else
 4455   matrix ret = I(n)
 4456   ret[sel, sel] = qrdecomp(mnormal(nk, nk))
 4457 endif
 4458 return ret
 4459 </code>
 4460 </gretl-function>
 4461 <gretl-function name="rot_redraw" type="void" private="1">
 4462  <params count="6">
 4463   <param name="b" type="bundleref"/>
 4464   <param name="mod" type="bundle" const="true"/>
 4465   <param name="DO_BAYES" type="bool"/>
 4466   <param name="B" type="matrix" const="true"/>
 4467   <param name="iXX" type="matrix" const="true"/>
 4468   <param name="rotate_only" type="matrix" const="true"/>
 4469  </params>
 4470 <code>if DO_BAYES
 4471   matrices Mats = drawnormwis(iXX, B, mod.Sigma, b.df) # changed df at source
 4472   matrix b.B = Mats[1]
 4473   matrix b.Sigma = Mats[2]
 4474 endif
 4475 matrix b.rot = rand_rotation(mod.n, rotate_only)
 4476 </code>
 4477 </gretl-function>
 4478 <gretl-function name="prepAis" type="matrices" private="1">
 4479  <params count="2">
 4480   <param name="mod" type="bundle" const="true"/>
 4481   <param name="shocks" type="matrix" const="true"/>
 4482  </params>
 4483 <code>nshocks = rows(shocks)
 4484 matrices Ais = array(nshocks)
 4485 loop i = 1..nshocks --quiet
 4486   matrix Ais[i] = selifr(mod.SRest, mod.SRest[,6] .= shocks[i])
 4487 endloop
 4488 return Ais
 4489 </code>
 4490 </gretl-function>
 4491 <gretl-function name="get_id_mat" type="matrix" private="1">
 4492  <params count="5">
 4493   <param name="mod" type="bundle" const="true"/>
 4494   <param name="shocks" type="matrix" const="true"/>
 4495   <param name="irfs" type="matrices" const="true"/>
 4496   <param name="Ais" type="matrices" const="true"/>
 4497   <param name="DO_NORM" type="bool"/>
 4498  </params>
 4499 <code># (removed redundant nshocks arg)
 4500 nshocks = rows(shocks)
 4501 matrix id_matrix = zeros(nshocks, mod.n)
 4502 loop i = 1..nshocks --quiet
 4503   id_matrix[i,] = check_irfs(irfs, Ais[i])
 4504 endloop
 4505 if DO_NORM
 4506   # printf &quot;Before normalization:\n%6.0f\n&quot;, id_matrix
 4507   id_matrix = normalize(id_matrix)
 4508   # printf &quot;After normalization:\n%6.0f\n&quot;, id_matrix
 4509 endif
 4510 return id_matrix
 4511 </code>
 4512 </gretl-function>
 4513 <gretl-function name="correspondence" type="matrix" private="1">
 4514  <params count="1">
 4515   <param name="A" type="matrix" const="true"/>
 4516  </params>
 4517 <code># Tag: sign restriction apparatus
 4518 # this function checks for a possible way to couple
 4519 # shocks with observables, given the &quot;candidate&quot; matrix A
 4520 nr = rows(A)
 4521 nc = cols(A)
 4522 if nr == 0
 4523   return {}
 4524 endif
 4525 r = 1
 4526 c = 0
 4527 x = 0
 4528 loop while r &lt;= nr --quiet
 4529   matrix candidate = A[r,]
 4530   matrix z = !(candidate .= 0)
 4531   if sumr(z) == 1
 4532     # candidate is ok
 4533     scalar x = selifc(candidate, z)
 4534     scalar c = selifc(seq(1,nc), z)
 4535     break
 4536   endif
 4537   r++
 4538 endloop
 4539 if c&gt;0
 4540   matrix rest = correspondence(A[-r, -c])
 4541   if rows(rest) &gt; 0
 4542     rest[,1] += rest[,1] .&gt;= r
 4543     rest[,2] += rest[,2] .&gt;= c
 4544   endif
 4545 else
 4546   matrix rest = {}
 4547 endif
 4548 return {r, c, x} | rest
 4549 </code>
 4550 </gretl-function>
 4551 <gretl-function name="normalize" type="matrix" private="1">
 4552  <params count="1">
 4553   <param name="A" type="matrix" const="true"/>
 4554  </params>
 4555 <code># Tag: sign restriction apparatus
 4556 # this function is for normalizing the &quot;candidate&quot; matrix A
 4557 # (-- TODO: maybe rename this function, to avoid confusion
 4558 # with the normalize option in the SVAR model bundle!?)
 4559 matrix U = correspondence(A)
 4560 matrix ret = zeros(rows(A), cols(A))
 4561 loop i = 1..rows(U) --quiet
 4562   r = U[i,1]
 4563   c = U[i,2]
 4564   if c != 0
 4565     x = U[i,3]
 4566     ret[r,c] = x
 4567   endif
 4568 endloop
 4569 return selifr(ret, maxr(abs(ret)) .&gt; 0)
 4570 </code>
 4571 </gretl-function>
 4572 <gretl-function name="check_one_irf" type="matrix" private="1">
 4573  <params count="3">
 4574   <param name="M" type="matrices"/>
 4575   <param name="spec" type="matrix"/>
 4576   <param name="verbose" type="int" default="0"/>
 4577  </params>
 4578 <code># Tag: sign restriction apparatus
 4579 # This function will return a row vector with elements
 4580 # that are either:
 4581 #  1 if the corresponding column satifies the bounds
 4582 # -1 if column satifies the bounds after a sign change
 4583 #  0 if the bounds condition isn't met
 4584 n = cols(M[1])
 4585 h = nelem(M)
 4586 vndx = spec[1] # always matrix?
 4587 # reshape IRFs so as to have columns as candidates
 4588 matrix tmp = mshape(flatten(M)[vndx,], n, h)'
 4589 # consider only selected rows
 4590 matrix sel = seq(spec[4], spec[5])
 4591 tmp = tmp[sel,]
 4592 Max = maxc(tmp)
 4593 Min = minc(tmp)
 4594 scalar lo = spec[2]
 4595 scalar hi = spec[3]
 4596 noflip = (Min .&gt; lo) &amp;&amp; (Max .&lt; hi) # shocks that respect the bounds
 4597 flip = (Max .&lt; -lo) &amp;&amp; (Min .&gt; -hi) # shocks that need to be flipped
 4598 ret = noflip - flip
 4599 if verbose
 4600   print tmp
 4601   if verbose &gt; 1
 4602     print lo hi
 4603     print Min Max
 4604   endif
 4605   print noflip flip ret
 4606 endif
 4607 return ret
 4608 </code>
 4609 </gretl-function>
 4610 <gretl-function name="check_irfs" type="matrix" private="1">
 4611  <params count="3">
 4612   <param name="M" type="matrices"/>
 4613   <param name="A" type="matrix"/>
 4614   <param name="verbose" type="int" default="0"/>
 4615  </params>
 4616 <code># Tag: sign restriction apparatus
 4617 # Here we examine the constraints on the IRFs
 4618 # We return a vector
 4619 # indicating if all the constraints contained in the matrix A are
 4620 # satisfied; each row of A contains
 4621 #
 4622 # 1    : ordinal no of variable to check
 4623 # 2, 3 : bounds (NA for +-infty)
 4624 # 4, 5 : IRF interval to check (0-based)
 4625 #
 4626 # It may happen that the IRFs are ok up to a sign swap; we flag this
 4627 # case by setting the corresponding column to -1
 4628 n = cols(M[1])
 4629 rA = rows(A)
 4630 # check consistency
 4631 matrix chks = zeros(rA, n)
 4632 loop i = 1..rA --quiet
 4633   chks[i,] = check_one_irf(M, A[i,], verbose)
 4634 endloop
 4635 if verbose
 4636   print chks
 4637 endif
 4638 matrix positives = minc(chks) .= 1
 4639 matrix negatives = maxc(chks) .= -1
 4640 return positives - negatives
 4641 </code>
 4642 </gretl-function>
 4643 <gretl-function name="check_id" type="scalar" private="1">
 4644  <params count="1">
 4645   <param name="id" type="matrix"/>
 4646  </params>
 4647 <code># Tag: sign restriction apparatus
 4648 # Here we check if the sign restrictions allow us to
 4649 # form a coherent set of IRFs for the problem at hand
 4650 #
 4651 # the input matrix has r rows and c columns, where r is the
 4652 # number of shocks we want to identify and c is the number
 4653 # of variables in the system; note that r can be smaller
 4654 # than c for partially identified models; for example, in Uhlig
 4655 # (2005) KME, c = 6 but r = 1 (the monetary policy shock)
 4656 #
 4657 # id can only contain 0s, 1s and -1s; therefore, for completely
 4658 # identified models, we just check if &quot;id&quot; is a proper orthogonal
 4659 # matrix; for partially identified models, we need to make sure
 4660 # that we have only one non-zero entry for each row, and that
 4661 # the rank of id is r.
 4662 r = rows(id)
 4663 c = cols(id)
 4664 if r &lt; c
 4665   # partial identification
 4666   ret = (rank(id) == r) &amp;&amp; minc(sumr(abs(id)) .= 1) == 1
 4667 else
 4668   ret = maxc(abs((id'id) - I(r))) &lt; 1.0e-15
 4669 endif
 4670 return ret
 4671 </code>
 4672 </gretl-function>
 4673 <gretl-function name="gen_SRirfs" type="matrices" private="1">
 4674  <params count="3">
 4675   <param name="b" type="bundle"/>
 4676   <param name="numh" type="int"/>
 4677   <param name="to_cum" type="matrix"/>
 4678  </params>
 4679 <code># This calculates results for 0..mod.horizons, which is actually
 4680 # like in the earlier SVAR code.
 4681 # TODO: Maybe merge this function with earlier ones from SVAR?
 4682 #
 4683 # So numh == mod.horizon + 1.
 4684 # (Memo: the doc has been (up to 2020) slightly wrong in this
 4685 #  regard; it claimed h=mod.horizons rows in the IRFs matrix,
 4686 #  one less than true)
 4687 matrix B = b.B
 4688 k = b.exoterms    # just a number?
 4689 # On the outside, this cumulation spec is taken from the mod.
 4690 # (Empty matrix to_cum if doesn't apply.)
 4691 any_cumul = nelem(to_cum) &gt; 0
 4692 n = cols(B)
 4693 p = (rows(B) - k) / n # VAR lags
 4694 matrix comp = B[k+1:,]'
 4695 if p &gt; 1
 4696   comp = comp | I(n*(p-1), n*p)
 4697 endif
 4698 matrix K = cholesky(b.Sigma) * b.rot
 4699 matrices ret = array(numh)
 4700 ret[1] = K	# impact effect
 4701 matrix tmp = I(rows(comp))
 4702 loop i = 2 .. numh --quiet
 4703   tmp = tmp * comp
 4704   matrix irf_i = tmp[1:n,1:n] * K
 4705   if any_cumul
 4706     irf_i[to_cum,] = irf_i[to_cum,] + ret[i-1][to_cum,]
 4707   endif
 4708   ret[i] = irf_i
 4709 endloop
 4710 return ret
 4711 </code>
 4712 </gretl-function>
 4713 <gretl-function name="exotic_inner_check" type="scalar" private="1">
 4714  <params count="4">
 4715   <param name="chk" type="scalarref"/>
 4716   <param name="i" type="int"/>
 4717   <param name="mod" type="bundle" const="true"/>
 4718   <param name="irfs" type="matrices" const="true"/>
 4719  </params>
 4720 <code># (chk is a boolean, but gretl doesn't allow pointerized int or bool)
 4721 # TODO: Check whether this setup is really compatible with a situation
 4722 # where we have only partial identification from the &quot;standard&quot;
 4723 # (non-exotic) restrictions, but the exotic restrictions then refer
 4724 # to additional shocks not covered before. (Somehow I'm worrying
 4725 # that the previous multiplication with the non-square id_matrix in the
 4726 # partial id case removes the needed irfs that would still have to be
 4727 # checked against the exotic restrictions. - Sven)
 4728 string exo_restrict = mod.exoticSR.checks[i]
 4729 scalar inip1 = mod.exoticSR.spans[i,1]
 4730 scalar finp1 = mod.exoticSR.spans[i,2]
 4731 scalar needs_model = mod.exoticSR.super[i]
 4732 out = 0
 4733 loop t = inip1..finp1 --quiet
 4734   matrix M = irfs[t] # hardcoded &quot;M&quot; for the API!
 4735   if needs_model
 4736     # FIXME This part looks very broken...!?
 4737     # or very outdated...
 4738     scalar exocheck = feval(exo_restrict, M, mod)
 4739   else
 4740     scalar exocheck = @exo_restrict
 4741   endif
 4742   if !exocheck
 4743     chk = 0
 4744     out = 1
 4745     if 0
 4746       printf &quot;Exotic check %s failed on M[%d]:\n&quot;, exo_restrict, t-1
 4747       printf &quot;%7.4f\n&quot;, M
 4748     endif
 4749     break
 4750   endif
 4751 endloop
 4752 return out
 4753 </code>
 4754 </gretl-function>
 4755 <sample-script filename="examples/simple_C.inp">
 4756 set verbose off
 4757 include SVAR.gfn
 4758 open sw_ch14.gdt
 4759 
 4760 series infl = 400*ldiff(PUNEW)
 4761 rename LHUR unemp
 4762 
 4763 list X = unemp infl
 4764 list Z = const
 4765 
 4766 Mod = SVAR_setup(&quot;C&quot;, X, Z, 3)
 4767 Mod.horizon = 36
 4768 SVAR_restrict(&amp;Mod, &quot;C&quot;, 1, 2)
 4769 
 4770 set stopwatch
 4771 SVAR_estimate(&amp;Mod)
 4772 printf &quot;Time (Cmodel) = %g\n&quot;, $stopwatch
 4773 
 4774 fevdmat = FEVD(&amp;Mod)
 4775 print fevdmat
 4776 #IRFplot(&amp;Mod, 1, 1)
 4777 IRFsave(&quot;simple_C_11_noboot.pdf&quot;, &amp;Mod, 1, 1)
 4778 
 4779 set stopwatch
 4780 bfail = SVAR_boot(&amp;Mod, 1024, 0.90)
 4781 printf &quot;Failed = %d, Time (bootstrap) = %g\n&quot;, bfail, $stopwatch
 4782 
 4783 IRFsave(&quot;simpleC.pdf&quot;, &amp;Mod, 0, 0)	# 0 means &quot;all&quot;
 4784 </sample-script>
 4785 </gretl-function-package>
 4786 </gretl-functions>