```
/* Returns a H x (1 + 2*Q) matrix of bootstrapped quantiles of
forecast (FC) densities;
'H' is the forecast horizon up to which the bootstrap is done.
The first column always refers to the median (0.5-quantile).
Q is the number of elements in 'quants', which should have values between 0 and 0.5.
The default for 'quants' is {0.005, 0.025, 0.05, 0.1}.
(Corresponding to two-sided 1%, 5%, 10%, 20%. Each value
is mirrored to the right: e.g. 0.1 together with 0.9.)
'exo' (Hx1) provides the future values of any exogenous terms
(other than the constant!) in the underlying equation (X'b, i.e.
already with the coefficients factored in). Otherwise
-- can pass null or {} as arg -- they are set to zero.
The bootstrap is based on resampling.
(Future possibilities:
- wild bootstrap, with flag in the bundle mod,)
*/
## pre-process the quantiles
if !exists(quants)
matrix quants = {} # unified treatment
endif
quants = procquants(quants)
## pre-process the exogenous terms
if !exists(exo)
matrix exo = {}
endif
if nelem(exo)
matrix Xb = vec(exo)
if rows(Xb) != H
funcerr "num of exo values must match horizon H"
endif
else
matrix Xb = zeros(H,1)
endif
## constant?
wconst = inlist(mod.xlist, const)
if wconst
Xb += mod.coeff[1] # const always comes first in $coeff
endif
# process the lag information
bundle mod2 = rhslags(mod)
# lag coefficients
lagorder = max(mod2.ylags)
matrix comptop = zeros(1, lagorder)
loop i = 1..rows(mod2.ylags) -q
whichcf = mod2.ylagspos[i]
comptop[1, mod2.ylags[i]] = mod.coeff[whichcf]
endloop
# initial values for the FC (last 'lagorder' values in-sample)
matrix init = { mod.ylist[1] } # the dep var as vector
init = init[rows(init) - lagorder + 1 : ] # final lagorder values
## the bootstrap
matrix res = { mod.uhat }
matrix FCboot = zeros(H, 0)
loop d = 1..reps -q
# resampling (non-wild so far...)
matrix rndix = mrandgen(i, 1, rows(res), H, 1)
matrix innov = res[rndix]
FCboot ~= varsimul(comptop, Xb + innov, init)[lagorder + 1: ,]
endloop
# evaluate the distribution
matrix out = transp(quantile(FCboot', quants))
colnames(out, sprintf("% g", 100*quants))
rownames(out, sprintf("% d", seq(1,H)))
return out
```

```
## use the uniFCboot() function to get the alpha-quantile
# of the FC density, which is the optimal FC under
# lin-lin loss with alpha (source Elliott & Timmermann p.43)
#
# 'mod' must be a bundle as given by the $model accessor.
# For construction of 'exo' see the uniFCboot function.
#
# Returns a H-row column vector of linlin point forecasts,
# from forecast step 1 (top) to H (bottom).
if !exists(exo)
matrix exo = {}
endif
matrix temp = uniFCboot(mod, H, reps, {alpha}, exo)
# (the 1st col always has the median, the 2nd then our alpha)
return temp[,2]
```

```
# Glue function that will be called from gretl's
# model window and has access to the $model accessor.
# It just passes that on.
# Returns the same bundle, augmented with the (H x (1+2Q)-matrix
# of bootstrapped forecast quantiles ("uniFCboot").
bundle mod = $model # copy the active model (from GUI model context)
/*
# check whether we have access to future data
...
*/
# pre-process optional args
if !exists(quants)
matrix quants = {}
endif
if !exists(exo)
matrix exo = {}
endif
if alpha != 0.5 # we do lin-lin point FC
semiGUI_FClinlin(&mod, alpha, H, reps, exo)
else # intervals for quantiles
matrix temp = uniFCboot(mod,H,reps,quants,exo)
# add result to bundle
matrix mod.uniFCboot = temp
printf "Bootstrapped forecasts for %s, up to step %d,\n", varnames(mod.ylist)[1], H
printf " with %d replications,\n", reps
printf " starting from sample endpoint %s.\n", obslabel($t2)
printf "Chosen quantiles (plus right-tail analogues): % g\n", procquants(quants)
print temp
endif
return mod
```

```
# (private function)
/* This function expects an input bundle of the type
that is provided by gretl's $model accessor.
It puts the information about the names of the
appearing lagged terms (ylagnames, strings array), which lags are used (ylags, matrix), and which is
their respective position in the regressor list
(ylagspos, matrix) into the returned bundle.
This only works if (a) the lagged terms are exactly
named as gretl's lags() function does it, (b) if
the ordering is strictly ascending (e.g. -1 first).
There can be gaps in the lag structure.
*/
bundle out = mod # copy
# create the max potential number of lags
list templags = lags(nelem(mod.xlist), mod.ylist)
# keep those lags that were actually used in $xlist
list sec = mod.xlist && templags
# and record their names
strings out.ylagnames = varnames(sec)
# determine which lags and where they appear
matrix out.ylags = {}
matrix out.ylagspos = {}
loop foreach s sec --quiet
out.ylags |= inlist(templags, sec.$s) # i-th pos in templags is lag
out.ylagspos |= inlist(mod.xlist, sec.$s) # where in $xlist/$coef ?
endloop
return out
```

```
return {0.005, 0.025, 0.05, 0.1}
```

```
if !exists(quants)
matrix quants = {} # for unified treatment
endif
if nelem(quants)
matrix q = sort(vec(quants))'
if min(q) <= 0 || max(q) >= 1
funcerr "invalid input for quants"
endif
else
matrix q = standardquants()
endif
# add the median
q = 0.5 ~ q ~ mreverse(1 - q')'
return q
```

```
# glue function that will be called (but only
# indirectly from gretl's model window.
# Augments the mod bundle with the H-vector
# of lin-lin optimal point forecasts ("FClinlin").
if !exists(exo)
matrix exo = {} # easier to hand down
endif
matrix temp = FClinlin(mod,alpha,H,reps,exo)
# add result to the bundle
matrix mod.FClinlin = temp
mod.alphalinlin = alpha
printf "Point forecasts for %s, up to step %d,\n", varnames(mod.ylist)[1], H
printf " with %d replications,\n", reps
printf " starting from sample endpoint %s,\n", obslabel($t2)
printf " based on lin-lin loss with asymmetry %g:\n", alpha
print temp
```